CN110706265A - Maneuvering target tracking method for improving SRCKF strong tracking filtering - Google Patents

Maneuvering target tracking method for improving SRCKF strong tracking filtering Download PDF

Info

Publication number
CN110706265A
CN110706265A CN201911072800.3A CN201911072800A CN110706265A CN 110706265 A CN110706265 A CN 110706265A CN 201911072800 A CN201911072800 A CN 201911072800A CN 110706265 A CN110706265 A CN 110706265A
Authority
CN
China
Prior art keywords
time
target
matrix
representing
filtering
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911072800.3A
Other languages
Chinese (zh)
Other versions
CN110706265B (en
Inventor
盛卫东
王雪莹
曾瑶源
安玮
程煜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201911072800.3A priority Critical patent/CN110706265B/en
Publication of CN110706265A publication Critical patent/CN110706265A/en
Application granted granted Critical
Publication of CN110706265B publication Critical patent/CN110706265B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/344—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Operations Research (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention relates to a maneuvering target tracking method for improving SRCKF strong tracking filtering, and belongs to the technical field of data association and tracking. The method comprises the following steps: s1: establishing a motion model and an observation model according to the motion condition of the target and the observation condition of the satellite; s2: initializing filtering; s3: carrying out filtering prediction; s4: performing a first measurement update; s5: calculating an fading factor with a time factor and carrying out threshold judgment; s6: and performing measurement updating again according to the threshold judgment result of the S5. The invention has the following beneficial effects: firstly, a more accurate target motion model does not need to be established, and the tracking effect on the maneuvering target can be achieved only by using the motion model when the target normally moves; and secondly, compared with a common adaptive filtering method, the tracking effect of the method is equivalent to that of a common filtering method without considering the maneuver when the target is not maneuvered, the noise interference is less, and the tracking result has more robustness.

Description

Maneuvering target tracking method for improving SRCKF strong tracking filtering
Technical Field
The invention relates to a maneuvering target tracking method for improving strong tracking filtering of a Square-root cubature Kalman Filter (SRCKF) aiming at optical satellite observation, belonging to the technical field of data association and tracking.
Background
In recent years, optical image satellites have been the focus of research in various countries because of their advantages such as wide observation regions, being not easily limited by national boundaries, and being not easily affected by electromagnetic interference, and in optical satellite observation, state estimation and situation awareness of spatial targets are also important points of research.
However, under the existing conditions, the motion pattern of the spatial target is increasingly complex, and acts such as bait throwing, track changing, diversion and the like are generated, so that the motion pattern deviates from the original motion pattern, a certain maneuvering phenomenon is generated, and further great difficulty is caused to the estimation and situation perception of the target state.
Meanwhile, the existing maneuvering target tracking method is generally divided into two types, the first type is to model the target motion, and then track the target by using an Interactive multi-mode (IMM) method, but this type of method needs to establish a more accurate target motion model, and often needs to consume a large amount of time in the model conversion process, the spacial target throwing bait, orbital transfer, ramification and the like are often short-time behaviors, and the motion model is relatively difficult to establish, so this type of method is not suitable for the tracking of the spacial target. The second type is a strong tracking filtering method based on a square root cubature Kalman filter, which is one of the methods, and the method does not change a motion model but changes filtering to generate adaptability to the motion of a target, but the method easily causes larger error at the initial stage of target observation, has poor robustness to observation noise and is not suitable for the whole process of space target tracking.
Disclosure of Invention
The invention aims to solve the technical problem of providing a maneuvering target tracking method for improving SRCKF strong tracking filtering aiming at the defects of space maneuvering target tracking in the prior art.
In order to solve the technical problems, the technical scheme adopted by the invention is as follows: a maneuvering target tracking method for improving SRCKF strong tracking filtering comprises the following steps:
s1, establishing a motion model and an observation model according to the motion condition of the target and the observation condition of the satellite;
s1.1, establishing a motion model according to the motion condition of a target, and considering that J is used as a missile target which moves in a free section and is only influenced by gravity under the non-maneuvering condition when the target is a missile target which moves in a free section2Modeling the target motion by the perturbation model;
the state vector of the target in the ECI coordinate system is X (k) ═ x (k), y (k), z (k), v (k)x(k),vy(k),vz(k)]TWherein x (k), y (k), z (k), vx(k),vy(k),vz(k) Respectively representing the position and velocity of the target at time k, with target J2The perturbation model is as follows:
Figure BDA0002261481400000021
wherein mue＝3.98604418m3/sec2Represents the constant of the gravity of the earth,
Figure BDA0002261481400000023
Re6378137m represents the radius of the earth, J20.001082626075 represents perturbation coefficient;
since the above equation is difficult to directly express as an equation form about x (k), under the condition that the target state at the time k is known as x (k), a fourth-order Runge-Kutta numerical integration method is considered to obtain the value at the time k +1, where:
X(k+1)＝f(X(k))+q(k)
wherein q (k) is a zero-mean gaussian vector, and T represents a time interval during discretization of a target state;
s1.2, an observation model is established according to the observation condition of a satellite, the satellite is in passive observation, only the azimuth angle and the pitch angle of a target can be measured, and under an ECI coordinate system, the state vector of the satellite at the moment k is assumed as follows:
Sa(k)＝[xs(k),ys(k),zs(k),vsx(k),vsy(k),vsz(k)]T
wherein xs(k),ys(k),zs(k)、vsx(k),vsy(k),vsz(k) Respectively representing the position and the velocity of the satellite at the time k;
with the direction of satellite velocity as XbAxis, taking the line connecting the satellite to the center of the earth as ZbAxis, giving Y according to the right-hand screw rulebAxis, definition XbAxis, YbAxis, ZbAxial vectors of the axes are SXb，SYb，SZbGiving a star coordinate system, αsRepresenting the azimuth of the target, betasRepresenting the target pitch angle and SX the satellite view, when alpha can be calculated as followss、βs：
SXb＝[vx(k),vy(k),vz(k)]T
SZb＝-[x(k),y(k),z(k)]T
SX＝[xs(k)-x(k),ys(k)-y(k),zs(k)-z(k)]T
To fit the filtering model, the above equation is generally written as:
wherein Z (k) represents a measurement of time k of the satellite, d βsAnd d αsRepresenting sensor measurement noise, r (k) is a zero-mean gaussian vector;
s2: in the filtering initialization, because the SRCKF method is used, in the filtering process, the iteration amount is not the covariance matrix but the square root matrix of the covariance matrix, so that the square root decomposition of the initial covariance matrix is needed in the initial value, that is, the filtering initialization process is:
Figure BDA0002261481400000031
wherein the content of the first and second substances,
Figure BDA0002261481400000032
representing the state estimation of the target initial time, E (-) representing the expectation of the matrix, P (0) representing the covariance matrix of the state estimation value of the target initial time, S (0) representing the square root of the covariance matrix of the target initial time, and chol (-) representing the lower triangular matrix obtained by performing Cholesky decomposition on the matrix;
meanwhile, because the change times of the fading factors need to be counted in S5, the change time flag is defined to be 0 at the initial moment;
s3: the process of performing filtering prediction, that is, extrapolating the system state at the time k +1 when the system state at the time k is known, still follows the framework of the SRCKF method, and includes:
for k ∈ {0, …, ∞ }:
calculating a volume point:
Figure BDA0002261481400000033
Figure BDA0002261481400000034
wherein
Figure BDA0002261481400000035
Representsk time volume point matrix vkColumn i, N represents the dimension of the state vector,represents the ith volume point and the ith volume point,
Figure BDA0002261481400000037
represents the estimated value of the system state at time k, S (k) represents the square root of the covariance matrix of the estimated value of the system state at time k, [1 ]]iRepresentative matrix [1 ]]The ith element in (1)]As follows:
Figure BDA0002261481400000038
the prediction process is as follows:
υk+1|k＝f(υk)
Figure BDA0002261481400000039
Se(k+1|k)＝qr([Vk+1|kSQ]T)
wherein
Figure BDA00022614814000000311
Representing the prediction of a k +1 time value volume point matrix upsilon by using k time volume point valuesk+1|kColumn i, Vk+1|kMatrix of representative volume points υk+1|kPredicted value at time k +1The deviation of (2), QR (·), is a function representing transposing the upper triangular matrix obtained by QR decomposition of the matrix, SQChol (q (k)) which represents the square root of the covariance matrix of the motion model noise in the prediction process, q (k)TRepresenting co-ordination of motion model noise during prediction of time kA variance matrix;
namely, the predicted value without fading factor at the moment of k +1 is obtained
Figure BDA0002261481400000041
And the square root of its covariance Se(k+1|k)；
S4: the invention is based on the improvement of a strong tracking filtering method, because the strong tracking filtering method is provided under the framework of expanding Kalman filtering, in the process of calculating the fading factors, a Jacobian matrix of a motion equation and an observation equation needs to be provided, in order to expand the Jacobian matrix under the framework of an SRCKF method without providing the Jacobian matrix, a measurement updating process needs to be added firstly, and the specific steps are as follows:
propagation volume point:
Figure BDA0002261481400000042
Figure BDA0002261481400000043
Figure BDA0002261481400000044
wherein
Figure BDA0002261481400000045
Representing a volume point matrix derived using the predicted value at time k +1 and the square root of the covariance
Figure BDA0002261481400000046
The (c) th column of (a),
Figure BDA0002261481400000047
matrix of representation versus volume points
Figure BDA0002261481400000048
Observing according to the observation equation of S1.2,
Figure BDA0002261481400000049
representative pair
Figure BDA00022614814000000410
Volume point matrix psi for observationk+1|kThe ith column;
the measurement update process without the fading factor is as follows:
Figure BDA00022614814000000411
Szz e(k+1)＝qr([Zk+1SR]T)
Figure BDA00022614814000000412
Pxz e(k+1)＝Xk+1Zk+1 T
Figure BDA00022614814000000413
Figure BDA00022614814000000414
wherein Zk+1Volume point matrix psi representing the time k +1k+1|kAnd a measured value of the prediction
Figure BDA00022614814000000415
Deviation of (A), Xk+1Volume point matrix representing k +1 time
Figure BDA00022614814000000416
Predicted value at time k +1
Figure BDA00022614814000000417
Deviation of (A), SRChol (R (k +1)) represents the square root of the covariance matrix resulting from observation equation noise at time k +1, R (k +1) ═ R (k +1)TThat is toThe method comprises the steps that a covariance matrix generated by observation equation noise at the moment k +1, z (k +1) represents a real measurement value at the moment k +1, Γ (k +1) represents an estimated value of innovation increment covariance, ρ is a forgetting factor, ρ is more than 0 and less than or equal to 1, and the value is generally ρ is 0.95.
I.e. obtaining the predicted measurement value without adding the fading factor
Figure BDA00022614814000000418
And the square root of its covariance Szz e(k +1) and a cross-covariance matrix P of predicted values and predicted metrology valuesxz e(k +1), and an innovation increment ε (k + 1);
s5: calculating an fading factor with a time factor and carrying out threshold judgment, wherein S obtained by S3 and S4 is requirede(k+1|k)、Szz e(k +1) and Pxz e(k +1), in this case:
hx(k+1)＝Pxz e(k+1)T/Se(k+1|k)/Se(k+1|k)T
Figure BDA0002261481400000051
Figure BDA0002261481400000052
Figure BDA0002261481400000053
where Norm (-) represents the sum of squares of the individual elements of the matrix, hx(k +1) represents the gain matrix of the square root of the covariance matrix of the noise of the motion model in the prediction process in the observation process at the moment k +1, and l (k +1) represents the time factor term of the moment k +1, wherein t is the time factor termbRepresenting the sampling interval, T, of apparent time measurementsRepresenting the convergence time, λ, of the satellite observations0Represents the initial calculation value of the fading factor, and lambda (k +1) represents the final value of the fading factor at the moment of k + 1;
in this case, if λ (k +1) > 1, the number of changes is countedAdding 1 to the flag, and comparing the flag with a preset threshold value nλComparison, nλSet to 5 if flag is greater than or equal to nλIf k is equal to k-nλI.e. is a forward trace back nλAt each moment, the change times flag is reset to 0, all the quantities containing k in the filtering process are changed into the value represented by the new moment, and the order is given
Figure BDA0002261481400000054
If flag is less than nλThen let S (k +1| k) ═ Se(k +1| k), continuing the normal filtering process; if λ (k +1) ═ 1, the number of changes flag is set to 0 as it is, and S (k +1| k) is made to be S at the same timee(k +1| k), continuing the normal filtering process;
s6: performing measurement updating again according to the threshold judgment result of S5;
updating the real measurement obtained at the time k +1 according to the value of S (k +1| k) obtained in S5:
propagation volume point:
Figure BDA0002261481400000055
Figure BDA0002261481400000056
Figure BDA0002261481400000057
the measurement update process is as follows:
Figure BDA0002261481400000058
Szz(k+1)＝qr([Zk+1SR]T)
Figure BDA0002261481400000059
Pxz(k+1)＝Xk+1Zk+1 T
W(k+1)＝Pxz(k+1)/Szz(k+1)/Szz(k+1)T
Figure BDA0002261481400000061
S(k+1)＝qr([Xk+1-W(k+1)Zk+1W(k+1)SR]T)
wherein
Figure BDA0002261481400000062
Representing the predicted measurement value at time k +1, Szz(k +1) represents
Figure BDA0002261481400000063
Square root of covariance, Pxz(k +1) represents a cross covariance matrix of the predicted value and the predicted measurement value at the time k +1, and W (k +1) represents a filter gain matrix at the time k + 1.
Namely, the state estimation vector at the moment of k +1 is obtained
Figure BDA0002261481400000064
And the square root of its covariance matrix S (k +1) so that iteration can continue to derive state estimates at all times, the algorithm flow diagram of the filtering process is shown in fig. 3.
Compared with the prior art, the invention has the beneficial effects that: firstly, a more accurate target motion model does not need to be established, and the tracking effect on the maneuvering target can be achieved only by using the motion model when the target normally moves; and secondly, compared with a common adaptive filtering method, the tracking effect of the method is equivalent to that of a common filtering method without considering the maneuver when the target is not maneuvered, the noise interference is less, and the tracking result has more robustness.
Drawings
FIG. 1 is a satellite observation model;
FIG. 2 is a flow chart of an implementation of a maneuvering target tracking method with improved SRCKF strong tracking filtering;
FIG. 3 is a flow chart of an implementation of a filtering process of a maneuvering target tracking method for improving SRCKF strong tracking filtering;
FIG. 4 is a comparison of target no maneuver position velocity RMSE: (a) position error (b) velocity error;
FIG. 5 is a comparison of position velocity RMSE for a long maneuver of interest: (a) position error (b) velocity error;
FIG. 6 is a comparison graph of the RMSE with localized amplification of the position velocity for long small maneuvers of the target: (a) position error (b) velocity error;
FIG. 7 is a comparison of position velocity RMSE for a short maneuver of interest: (a) position error (b) velocity error;
FIG. 8 is a comparison graph of the RMSE for a localized magnification of position velocity for a short duration maneuver of the target: (a) position error (b) velocity error;
Detailed Description
The following further describes embodiments of the present invention with reference to the drawings.
The technical scheme of the invention is as follows: a maneuvering target tracking method for improving SRCKF strong tracking filtering specifically comprises the following steps:
s1: establishing a motion model and an observation model according to the motion condition of the target and the observation condition of the satellite;
s2: initializing filtering;
s3: carrying out filtering prediction;
s4: performing a first measurement update;
s5: calculating an fading factor with a time factor and carrying out threshold judgment;
s6: and performing measurement updating again according to the threshold judgment result of the S5.
In order to verify the effectiveness of the method, three scenes are selected, namely no maneuvering of the target occurs, long-time small maneuvering occurs and short-time large maneuvering occurs, the method is compared with a UKF filtering method and an SRCKF strong tracking filtering method which do not consider maneuvering, the Monte Carlo simulation times are 100 times, Root Mean Square Error (RMSE) is taken as an effect measurement standard, and a position speed RMSE comparison graph of different algorithms is given.
Fig. 4 is a comparison graph of the position velocity RMSE of a target without maneuver, and it can be seen from the graph that the velocity and position convergence results are the same as those of the UKF filtering method under the condition that the target does not have maneuver, while the SRCKF strong tracking filtering method determines the initial time to be in a moving state by mistake because the initial time estimation value is inaccurate, so that the convergence velocity is very slow at the beginning, and the final convergence effect is not as good as those of the other two, and meanwhile, the SRCKF strong tracking filtering method is easily affected by noise during the convergence process, and the result still fluctuates greatly under the condition of 100 monte carlo.
FIG. 5 and FIG. 6 are a comparison diagram of the position speed RMSE of the target long-time small maneuver and a partially enlarged diagram thereof, wherein the default maneuver time is 400s, the maneuver time is 100s, and the maneuver acceleration isAs can be seen from the figure, compared with the UKF filtering method without considering the maneuver, the method and the SRCKF strong tracking filtering method of the present invention can both generate a response to the maneuver of the target in a short time and can converge to a smaller value after the maneuver of the target stops; compared with the SRCKF strong tracking filtering, the method disclosed by the invention has the advantages that the divergence speed is lower after the target maneuvers, the tracking precision of the target maneuvers is better, the SRCKF strong tracking filtering generates volatility in the tracking process, and the method is more robust.
FIG. 7 and FIG. 8 are a comparison diagram of the position velocity RMSE and a partially enlarged view thereof, respectively, of a target short-time large maneuver, the default maneuver time being 400s, the maneuver duration being 20s, and the maneuver acceleration being
Figure BDA0002261481400000072
As can be seen from the figure, compared with the UKF filtering method without considering the maneuvering, the method and the SRCKF strong tracking filtering method can both respond to the maneuvering of the target in a short time and can converge to the ratio after the maneuvering of the target stopsA smaller value; compared with the SRCKF strong tracking filtering, the method has the advantages that the divergence speed is low after the target maneuvers, the tracking precision of the target maneuvers is better, the SRCKF strong tracking filtering generates volatility in the tracking process, and the method is more robust.
The experimental result shows that the method is not only suitable for the maneuvering target tracking, but also has better estimation on the state of the unmovable target, which indicates that the method is effective.

Claims (6)

1. A maneuvering target tracking method for improving SRCKF strong tracking filtering is characterized by comprising the following steps:
s1, establishing a motion model and an observation model according to the motion condition of the target and the observation condition of the satellite;
s1.1, establishing a motion model according to the motion condition of a target, and considering that J is used as a missile target which moves in a free section and is only influenced by gravity under the non-maneuvering condition when the target is a missile target which moves in a free section2Modeling the target motion by the perturbation model;
s1.2, an observation model is established according to the observation condition of a satellite, the satellite is in passive observation, only the azimuth angle and the pitch angle of a target can be measured, and under an ECI coordinate system, the state vector of the satellite at the moment k is assumed as follows:
Sa(k)＝[xs(k),ys(k),zs(k),vsx(k),vsy(k),vsz(k)]T
wherein xs(k),ys(k),zs(k)、vsx(k),vsy(k),vsz(k) Respectively representing the position and the velocity of the satellite at the time k;
with the direction of satellite velocity as XbAxis, taking the line connecting the satellite to the center of the earth as ZbAxis, giving Y according to the right-hand screw rulebAxis, definition XbAxis, YbAxis, ZbAxial vectors of the axes are SXb，SYb，SZbGiving a star coordinate system, αsRepresenting the azimuth of the target, betasRepresenting the target elevation angle, SX representing the satellite line of sight, which may be measured as followsCalculating alphas、βs：
SXb＝[vx(k),vy(k),vz(k)]T
SZb＝-[x(k),y(k),z(k)]T
SX＝[xs(k)-x(k),ys(k)-y(k),zs(k)-z(k)]T
Figure FDA0002261481390000011
To fit the filtering model, the above equation is generally written as:
Figure FDA0002261481390000012
wherein Z (k) represents a measurement of time k of the satellite, d βsAnd d αsRepresenting sensor measurement noise, r (k) is a zero-mean gaussian vector;
s2: and (3) carrying out filtering initialization, wherein the filtering initialization process comprises the following steps:
Figure FDA0002261481390000013
wherein the content of the first and second substances,
Figure FDA0002261481390000014
representing the state estimation of the target initial time, E (-) representing the expectation of the matrix, P (0) representing the covariance matrix of the state estimation value of the target initial time, S (0) representing the square root of the covariance matrix of the target initial time, and chol (-) representing the lower triangular matrix obtained by performing Cholesky decomposition on the matrix;
meanwhile, because the change times of the fading factors need to be counted in S5, the change time flag is defined to be 0 at the initial moment;
s3: the process of performing filtering prediction, that is, extrapolating the system state at the time k +1 when the system state at the time k is known, still follows the framework of the SRCKF method, and includes:
for k ∈ {0, …, ∞ }:
calculating a volume point:
Figure FDA0002261481390000022
wherein
Figure FDA0002261481390000023
Volume point matrix v representing k timekColumn i, N represents the dimension of the state vector,
Figure FDA0002261481390000024
represents the ith volume point and the ith volume point,represents the estimated value of the system state at time k, S (k) represents the square root of the covariance matrix of the estimated value of the system state at time k, [1 ]]iRepresentative matrix [1 ]]The ith element in (1)]As follows:
Figure FDA0002261481390000026
the prediction process is as follows:
υk+1|k＝f(υk)
Figure FDA0002261481390000027
Figure FDA0002261481390000028
Se(k+1|k)＝qr([Vk+1|kSQ]T)
wherein
Figure FDA0002261481390000029
Representing the prediction of a k +1 time value volume point matrix upsilon by using k time volume point valuesk+1|kColumn i, Vk+1|kMatrix of representative volume points υk+1|kPredicted value at time k +1
Figure FDA00022614813900000210
The deviation of (2), QR (·), is a function representing transposing the upper triangular matrix obtained by QR decomposition of the matrix, SQChol (q (k)) which represents the square root of the covariance matrix of the motion model noise in the prediction process, q (k)TRepresenting a covariance matrix of motion model noise in a k-moment prediction process;
namely, the predicted value without fading factor at the moment of k +1 is obtained
Figure FDA00022614813900000211
And the square root of its covariance Se(k+1|k)；
S4: the first measurement update is carried out, and the specific steps are as follows:
propagation volume point:
Figure FDA00022614813900000212
Figure FDA00022614813900000213
wherein
Figure FDA00022614813900000215
Representing a volume point matrix derived using the predicted value at time k +1 and the square root of the covariance
Figure FDA00022614813900000216
The (c) th column of (a),
Figure FDA00022614813900000217
matrix of representation versus volume points
Figure FDA00022614813900000218
Observing according to the observation equation of S1.2,representative pair
Figure FDA00022614813900000220
Volume point matrix psi for observationk+1|kThe ith column;
the measurement update process without the fading factor is as follows:
Figure FDA0002261481390000031
Szz e(k+1)＝qr([Zk+1SR]T)
Pxz e(k+1)＝Xk+1Zk+1 T
Figure FDA0002261481390000034
wherein Zk+1Volume point matrix psi representing the time k +1k+1|kAnd a measured value of the prediction
Figure FDA0002261481390000035
Deviation of (A), Xk+1Volume point matrix representing k +1 time
Figure FDA0002261481390000036
Predicted value at time k +1Deviation of (A), SRChol (R (k +1)) represents the square root of the covariance matrix resulting from observation equation noise at time k +1, R (k +1) ═ R (k +1)TNamely a covariance matrix generated by observation equation noise at the moment k +1, z (k +1) represents a real measurement value at the moment k +1, Γ (k +1) represents an estimated value of innovation increment covariance, and ρ is a forgetting factor;
i.e. obtaining the predicted measurement value without adding the fading factor
Figure FDA0002261481390000038
And the square root of its covariance Szz e(k +1) and a cross-covariance matrix P of predicted values and predicted metrology valuesxz e(k +1), and an innovation increment ε (k + 1);
s5: calculating an fading factor with a time factor and carrying out threshold judgment, wherein S obtained by S3 and S4 is requirede(k+1|k)、Szz e(k +1) and Pxz e(k +1), in this case:
hx(k+1)＝Pxz e(k+1)T/Se(k+1|k)/Se(k+1|k)T
Figure FDA0002261481390000039
Figure FDA00022614813900000310
where Norm (-) represents the sum of squares of the individual elements of the matrix, hx(k +1) represents the gain matrix of the square root of the covariance matrix of the noise of the motion model in the prediction process in the observation process at the moment k +1, and l (k +1) represents the time factor term of the moment k +1, wherein t is the time factor termbRepresenting the sampling interval, T, of apparent time measurementsRepresenting the convergence time, λ, of the satellite observations0Represents the initial calculation value of the fading factor, and lambda (k +1) represents the final value of the fading factor at the moment of k + 1;
if lambda (k +1) > 1, adding 1 to the change time flag, and comparing the flag with a preset threshold value nλComparing, if flag is greater than or equal to nλLet k be k-nλI.e. is a forward trace back nλAt each moment, the change times flag is reset to 0, all the quantities containing k in the filtering process are changed into the value represented by the new moment, and the order is given
Figure FDA00022614813900000312
If flag is less than nλThen let S (k +1| k) ═ Se(k +1| k), continuing the normal filtering process; if λ (k +1) ═ 1, the number of changes flag is set to 0 as it is, and S (k +1| k) is made to be S at the same timee(k +1| k), continuing the normal filtering process;
s6: performing measurement updating again according to the threshold judgment result of S5;
updating the real measurement obtained at the time k +1 according to the value of S (k +1| k) obtained in S5:
propagation volume point:
Figure FDA0002261481390000041
Figure FDA0002261481390000042
the measurement update process is as follows:
Figure FDA0002261481390000044
Szz(k+1)＝qr([Zk+1SR]T)
Figure FDA0002261481390000045
Pxz(k+1)＝Xk+1Zk+1 T
Figure FDA0002261481390000046
W(k+1)＝Pxz(k+1)/Szz(k+1)/Szz(k+1)T
Figure FDA0002261481390000047
S(k+1)＝qr([Xk+1-W(k+1)Zk+1W(k+1)SR]T)
wherein
Figure FDA0002261481390000048
Representing the predicted measurement value at time k +1, Szz(k +1) represents
Figure FDA0002261481390000049
Square root of covariance, Pxz(k +1) represents a cross covariance matrix of the predicted value and the predicted measurement value at the time k +1, and W (k +1) represents a filter gain matrix at the time k + 1;
namely, the state estimation vector at the moment of k +1 is obtained
Figure FDA00022614813900000410
And the square root of its covariance matrix S (k +1), so that iteration can continue to derive state estimates at all times.
2. A maneuvering target tracking method for improving SRCKF strong tracking filtering according to claim 1, characterized by: in S1.1, use is made of J2The perturbation model models the motion of the target by the following steps:
the state vector of the target in the ECI coordinate system is X (k) ═ x (k), y (k), z (k), v (k)x(k),vy(k),vz(k)]TWherein x (k), y (k), z (k), vx(k),vy(k),vz(k) Respectively representing the position and velocity of the target at time k, with target J2The perturbation model is as follows:
Figure FDA00022614813900000411
wherein mue＝3.98604418m3/sec2Represents the constant of the gravity of the earth,
Figure FDA0002261481390000051
Figure FDA0002261481390000052
Re6378137m represents the radius of the earth, J20.001082626075 represents perturbation coefficient.
3. A maneuvering target tracking method for improving SRCKF strong tracking filtering according to claim 2, characterized by: under the condition that the target state at the moment k is known to be X (k), a fourth-order Runge-Kutta numerical integration method is used for obtaining a value at the moment k +1, wherein the method comprises the following steps:
X(k+1)＝f(X(k))+q(k)
where q (k) is a zero mean gaussian vector and T represents the time interval during which the target state is discretized.
4. A maneuvering target tracking method for improving SRCKF strong tracking filtering according to claim 1, characterized by: in S4, the forgetting factor rho is greater than 0 and less than or equal to 1.
5. The maneuvering target tracking method for improving SRCKF strong tracking filtering according to claim 4, characterized by: in S4, the forgetting factor ρ is set to 0.95.
6. A maneuvering target tracking method for improving SRCKF strong tracking filtering according to claim 1, characterized by: in S5, a predetermined threshold value nλSet to 5.
CN201911072800.3A 2019-11-05 2019-11-05 Maneuvering target tracking method for improving SRCKF strong tracking filtering Active CN110706265B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911072800.3A CN110706265B (en) 2019-11-05 2019-11-05 Maneuvering target tracking method for improving SRCKF strong tracking filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911072800.3A CN110706265B (en) 2019-11-05 2019-11-05 Maneuvering target tracking method for improving SRCKF strong tracking filtering

Publications (2)

Publication Number Publication Date
CN110706265A true CN110706265A (en) 2020-01-17
CN110706265B CN110706265B (en) 2022-02-01

Family

ID=69205071

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911072800.3A Active CN110706265B (en) 2019-11-05 2019-11-05 Maneuvering target tracking method for improving SRCKF strong tracking filtering

Country Status (1)

Country Link
CN (1) CN110706265B (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113125962A (en) * 2021-04-21 2021-07-16 东北大学 Lithium titanate battery state estimation method under temperature and time variation
CN113466848A (en) * 2021-05-22 2021-10-01 中国人民解放军空军工程大学 Angle flicker noise scene-oriented co-location MIMO radar multi-target tracking resource optimal allocation method
CN113744309A (en) * 2021-08-11 2021-12-03 中国科学院空天信息创新研究院 Maneuvering target tracking method based on motion state change perception
CN116738873A (en) * 2023-05-11 2023-09-12 北京科技大学 Three-dimensional target tracking method and device based on double UKF and aerostat state estimation
CN117455960A (en) * 2023-12-22 2024-01-26 中国人民解放军海军工程大学 Passive positioning filtering algorithm for airborne photoelectric system to ground under time-varying observation noise condition

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030200065A1 (en) * 2001-04-20 2003-10-23 Li Luo Wen Maneuvering target tracking method via modifying the interacting multiple model (IMM) and the interacting acceleration compensation (IAC) algorithms
CN103294897A (en) * 2013-05-09 2013-09-11 哈尔滨工程大学 Self-adaptive filtering method used for ship dynamic-positioning position reference system
CN105652297A (en) * 2014-11-15 2016-06-08 航天恒星科技有限公司 Method and system for realizing real-time orbit determination for single satellite navigation positioning system
CN107390199A (en) * 2017-09-20 2017-11-24 哈尔滨工业大学（威海） A kind of radar maneuvering target tracking waveform design method
CN109520503A (en) * 2018-11-27 2019-03-26 南京工业大学 A kind of square root volume Fuzzy Adaptive Kalman Filtering SLAM method

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030200065A1 (en) * 2001-04-20 2003-10-23 Li Luo Wen Maneuvering target tracking method via modifying the interacting multiple model (IMM) and the interacting acceleration compensation (IAC) algorithms
CN103294897A (en) * 2013-05-09 2013-09-11 哈尔滨工程大学 Self-adaptive filtering method used for ship dynamic-positioning position reference system
CN105652297A (en) * 2014-11-15 2016-06-08 航天恒星科技有限公司 Method and system for realizing real-time orbit determination for single satellite navigation positioning system
CN107390199A (en) * 2017-09-20 2017-11-24 哈尔滨工业大学（威海） A kind of radar maneuvering target tracking waveform design method
CN109520503A (en) * 2018-11-27 2019-03-26 南京工业大学 A kind of square root volume Fuzzy Adaptive Kalman Filtering SLAM method

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张恒 等: "改进的强跟踪容积卡尔曼滤波的机动目标跟踪", 《现代防御技术》 *
袁晓波 等: "强跟踪自适应SRCKF的卫星姿态确定算法", 《测绘科学》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113125962A (en) * 2021-04-21 2021-07-16 东北大学 Lithium titanate battery state estimation method under temperature and time variation
CN113466848A (en) * 2021-05-22 2021-10-01 中国人民解放军空军工程大学 Angle flicker noise scene-oriented co-location MIMO radar multi-target tracking resource optimal allocation method
CN113466848B (en) * 2021-05-22 2023-09-19 中国人民解放军空军工程大学 Co-location MIMO radar multi-target tracking resource optimization allocation method for angular flicker noise scene
CN113744309A (en) * 2021-08-11 2021-12-03 中国科学院空天信息创新研究院 Maneuvering target tracking method based on motion state change perception
CN113744309B (en) * 2021-08-11 2024-05-24 中国科学院空天信息创新研究院 Maneuvering target tracking method based on motion state change perception
CN116738873A (en) * 2023-05-11 2023-09-12 北京科技大学 Three-dimensional target tracking method and device based on double UKF and aerostat state estimation
CN116738873B (en) * 2023-05-11 2024-02-06 北京科技大学 Three-dimensional target tracking method and device based on double UKF and aerostat state estimation
CN117455960A (en) * 2023-12-22 2024-01-26 中国人民解放军海军工程大学 Passive positioning filtering algorithm for airborne photoelectric system to ground under time-varying observation noise condition
CN117455960B (en) * 2023-12-22 2024-03-08 中国人民解放军海军工程大学 Passive positioning filtering method for airborne photoelectric system to ground under time-varying observation noise condition

Also Published As

Publication number Publication date
CN110706265B (en) 2022-02-01

Similar Documents

Publication Publication Date Title
CN110706265B (en) Maneuvering target tracking method for improving SRCKF strong tracking filtering
CN111985093B (en) Adaptive unscented Kalman filtering state estimation method with noise estimator
Kleeman Understanding and applying Kalman filtering
CN108225337B (en) Star sensor and gyroscope combined attitude determination method based on SR-UKF filtering
Morelli Real-time parameter estimation in the frequency domain
CN111798491B (en) Maneuvering target tracking method based on Elman neural network
Gross et al. A comparison of extended Kalman filter, sigma-point Kalman filter, and particle filter in GPS/INS sensor fusion
US7277047B1 (en) Reduced state estimation with biased measurements
CN113792411A (en) Spacecraft attitude determination method based on central error entropy criterion unscented Kalman filtering
CN110989341B (en) Constraint auxiliary particle filtering method and target tracking method
CN115204212A (en) Multi-target tracking method based on STM-PMBM filtering algorithm
CN116047498A (en) Maneuvering target tracking method based on maximum correlation entropy extended Kalman filtering
CN111798494A (en) Maneuvering target robust tracking method under generalized correlation entropy criterion
CN113449384A (en) Attitude determination method based on central error entropy criterion extended Kalman filtering
CN113587926A (en) Spacecraft space autonomous rendezvous and docking relative navigation method
He et al. A dynamic enhanced robust cubature Kalman filter for the state estimation of an unmanned autonomous helicopter
CN114445459B (en) Continuous-discrete maximum correlation entropy target tracking method based on variable decibel leaf theory
CN112595319B (en) Model self-adaptive compensation return trajectory estimation method
CN114858166A (en) IMU attitude resolving method based on maximum correlation entropy Kalman filter
Wang et al. An improved robust estimation method for gnss/sins under gnss-challenged environment
CN114132531A (en) Low-orbit space target orbit correction method and device and electronic equipment
Sun et al. A modified marginalized Kalman filter for maneuvering target tracking
Girrbach et al. Adaptive compensation of measurement delays in multi-sensor fusion for inertial motion tracking using moving horizon estimation
Ye et al. Iterative noise estimation-based cubature Kalman filtering for distributed pos in aerial earth observation imaging
Yue et al. INS/VNS fusion based on unscented particle filter

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant