Maneuvering target tracking method based on non-stationary Gaussian process

Abstract

The invention discloses a maneuvering target tracking method based on a non-stationary Gaussian process, and the method comprises the steps: carrying out the modeling of a target motion mode through the non-stationary Gaussian process, carrying out the online learning recursive updating of a target state and model parameters, and achieving the simultaneous state estimation and change point detection; the method comprises the following steps: dividing time series data into a non-overlapping observation set by using an operation length based on Bayesian online change point detection, namely modeling change point positions by using the operation length so as to determine partitions; according to the method, a non-stationary process in a maneuvering target tracking scene is converted into a partition stationary process, and for each partition stationary process, modeling is performed on a time sequence stationary process of a partition by using a Gaussian process, so that the robustness of maneuvering target tracking is improved; simulation shows that compared with an existing recursive Gaussian process motion tracker method, the method has the advantage of being robust in tracking effect.

Description

Technical Field

[0001] The invention belongs to the single-sensor single-maneuvering target tracking subsystem technology in the field of target tracking, and relates to a maneuvering target tracking method based on a non-stationary Gaussian process, which is suitable for a dynamic system with an unknown target maneuvering mode. Background Art

[0002] Target tracking is a basic task in sensor applications such as radar, sonar, and navigation. That is, in a continuous data sequence, the position relationship of the target to be tracked is established to obtain the complete motion trajectory of the target. In the Bayesian framework, target tracking is often solved by recursively updating the posterior probability density function (PDF) of the target state.

[0003] For linear target state and measurement models, the famous Kalman filter (KF) provides a recursive optimal solution. For general nonlinear state space models, typical Bayesian filters include extended Kalman filter (EKF), unscented Kalman filter (UKF) and interacting multiple model (IMM) filter. However, the above methods all rely on the underlying target dynamic model, that is, the preset target dynamic model needs to match the actual motion of the target during maneuvering. Therefore, when the target undergoes unknown or mixed maneuvers, the motion model may have significant uncertainty, making it impossible to approximate the target state through an appropriate mathematical model, resulting in a decrease in the tracking performance based on the Bayesian filter.

[0004] In recent years, target tracking technology based on online model learning of Gaussian process has become a new research hotspot by learning the target's motion pattern from data and has the characteristics of model adaptation. Aftab W et al. proposed a non-recursive model-free filtering method based on Gaussian process (GP), namely Gaussian process motion tracker (GPMT), which uses historical measurement data to predict and filter position and high-order position derivatives. Gaussian process motion tracker (GPMT) includes training and learning stages. It models an infinite number of target trajectories, but the learning stage in GPMT is not recursive. Based on GPMT, recursive Gaussian process motion tracker (RGPMT) is used to perform online model-free point target tracking. The state estimation and learning of Gaussian process hyperparameters are performed recursively.

[0005] However, the recursive Gaussian process motion tracker mainly solves the tracking problem of stationary time series processes. When the target changes from the steady-state condition, that is, a large maneuver occurs, the tracking performance of the tracker may deteriorate. Summary of the invention

[0006] In order to solve the problems in existing maneuvering target tracking methods that filtering accuracy affects filtering effect and the solution process is complicated, the present invention proposes a maneuvering target tracking method based on a non-stationary Gaussian process, which effectively handles the state estimation problem and improves the maneuvering target tracking performance.

[0007] In order to achieve the above tasks, the present invention adopts the following technical solutions:

[0008] A method for tracking a maneuvering target based on a non-stationary Gaussian process comprises the following steps:

[0009] Step 1: The subject constructs an online learning target motion model through Gaussian process to describe the possible motion mode of the maneuvering target, constructs a measurement model to represent the measurement data obtained by the sensor, and constructs a change point model to represent the time when the target motion mode changes; the transition probability of the change point model is given and the parameters to be estimated at the current moment are determined;

[0010] Step 2: The subject obtains the measurement data of the target through the sensor, and uses the measurement data of the previous multiple moments as the initial training data of each running length for online learning. The initial hyperparameters can be obtained through Gaussian process regression training for subsequent recursive estimation;

[0011] Step 3, the subject determines the joint distribution of the measurement data and the run length at the previous k moments according to the change point model;

[0012] Step 4, the subject determines the posterior distribution of the running length of the parameter to be estimated at time k according to the joint distribution;

[0013] Step 5: The subject calculates the posterior distribution of the target state quantity of the parameter to be estimated, which includes the state estimation mean and state estimation variance of the target at time k; the state estimation mean includes the required position estimation of the maneuvering target, thereby realizing the tracking of the maneuvering target.

[0014] Furthermore, the step 1 specifically includes:

[0015] Step 1.1, introduce a non-negative discrete integer variable r k ∈[1, k] as the change point model to represent the number of moments from the last change point of the measured data sequence to time k, that is, the running length at time k; when the target motion mode changes, let r k =1, at this time, the measured data sequence produces a change point; when the target motion mode does not change, r k =r k-1 +1, the measurement data sequence continues to grow;

[0016] Step 1.2: The duration between adjacent change points in the measured data sequence is a run length. In each run length, the target motion model is expressed as f r (·) is the state transition function, Among them, k represents the current moment, u and u′ represent the Gaussian process covariance kernel function Any two different inputs, u k Represents the input at the current k moment; is the kth time in a run of length r k The target state variable when the running length is r k The estimated position in the x direction or the estimated position in the y direction; The running length at time k is r k The Gaussian process at time , its mean function is 0, and the covariance kernel function is The covariance kernel function is usually a squared exponential kernel (SE) function;

[0017] Step 1.3: In each running length, the original measurement data obtained by the sensor is preprocessed by coordinate transformation and registration to obtain the position measurement of the target. Then the target measurement model can be expressed as Among them, z k is the measurement data at time k, indicating the position measurement data in the x direction or the position measurement data in the y direction; is the measurement noise of the sensor, which is zero-mean Gaussian distribution noise. is the measurement noise variance;

[0018] Step 1.4, for a certain run length r k , based on the motion model representation in step 1.2, at time k-1, the induced point set is introduced in,[·] T The T in the upper right corner represents vector or matrix transpose; represents the dth induction point, where d represents the number of induction point sets; represents the input set of induced point sets at time k-1, represents the input of the d-th induction point at time k-1;

[0019] Initial induction point set is the measured data at d moments before the start of the run length; since the target will maneuver during the motion, the hyperparameter vector Unknown, it is necessary to learn hyperparameters while estimating the target state; among them, represents the hyperparameter of the Gaussian process covariance kernel function at time k-1, σ k-1 represents the measurement noise parameter at time k-1; the target state induction point set at time k-1 Hyperparameter vector and the running length for k time is r k The predicted status at As an augmented state vector ψ k , construct the state space model as , the specific meanings are as follows:

[0020]

[0021] Among them, ψ k is the augmented state vector, the subscript k-1 represents the time, I k and I3 represent N×N and 3×3 identity matrices respectively, O represents an all-zero vector or matrix, is the input at time k; represents the cross covariance matrix between the predicted state at time k and the induced point set at time k-1, represents the variance of the predicted state at time k, represents the autocovariance between the induced point sets at time k-1, yes The transpose of

[0022] Step 1.5, the transition probability of the change point model is set to: run length r k-1 to r k The transition probability is p(r k |r k-1 ), when r k =1, p(r k |r k-1 )=H(r k-1 +1); otherwise r k =r k-1 +1, p(r k |r k-1 )=1-H(r k-1 +1); H(·) represents the risk function of the change point model;

[0023] Step 1.6, select (x k , r k ) as the parameter to be estimated at time k, x k is the target state quantity to be estimated at the current k moment, indicating the position of the target in the x direction or the y direction; r k Indicates the running length at the current k moment.

[0024] Furthermore, the step 2 specifically includes:

[0025] When k>d, the subject performs all runs of length r at each moment according to the target motion model defined in step 1. k Calculate the measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , rk ), the posterior prediction density is calculated with a run length of r k and the measurement data z at the previous k-1 moments 1：(k-1) The conditions are as follows:

[0026] Step 2.1, calculate the run length r from the state space model in step 1.4 k The augmented state vector ψ k The predicted mean and the corresponding predicted covariance Since the state-space equation in step 1.4 is given by the hyperparameter It is linear only when , so it is necessary to use unscented transformation to select the hyperparameter vector The mean of s sigma point; given the i Sigma points Its weight is The predicted state mean corresponding to this sigma point and the corresponding covariance matrix as follows:

[0027]

[0028] in, Respectively represent A in the state space model k , The values ​​of the hyperparameters involved are and They represent the estimated mean and estimated covariance of the state of the induced point set at time k-1 respectively; and denote the estimated mean and estimated covariance of the hyperparameter vector respectively; is the estimated cross covariance between the hyperparameter vector and the set of induced points, for The transpose of ; O represents a zero matrix or a zero vector; the superscript -1 of the parameter represents the matrix inversion, and the superscript T represents the matrix transpose;

[0029] After weighted fusion, the augmented state vector ψ can be calculated k The predicted mean and the predicted covariance for:

[0030]

[0031] Step 2.2, for each run length r k , the augmented state vector ψ k Decompose into observable vectors and the unobservable vector where σ k represents the measurement noise parameter at time k, represents the hyperparameter of the Gaussian process covariance kernel function at time k, represents the set of induced points at time k, then the state vector ψ k The predicted mean and the predicted covariance The specific decomposition form is as follows:

[0032]

[0033] Among them, the lower right corner scripts k and k-1 represent the time; the superscript ~ represents the predicted value; represents the unobservable vector g k-1 The estimated mean of ; Observable vector o k The predicted mean of , which includes the measurement noise parameter σ k-1 The estimated mean of and the run length is r k The predicted status at The predicted mean represents the unobservable vector g k-1 The estimated covariance of Represents the observable vector o k The prediction covariance of , which includes the measurement noise parameter σ k-1 The estimated variance of , prediction status The prediction variance and the measurement noise parameter σ k and predicted status The cross covariance between and for The transpose of represents the unobservable variable g k and observable variable o k The cross covariance between for The transpose of

[0034] Step 2.3, using the measurement model in step 1.3 and the decomposition result in step 2.2, calculate the run length r k The measured prediction mean at and the corresponding predicted covariance

[0035]

[0036] Step 2.4, for all runs of length r k, based on the measurement data z at the previous k-1 moments 1：(k-1) Calculate the measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , r k ):

[0037]

[0038] in, and are the measured prediction mean and the corresponding prediction covariance in step 2.3, respectively. Represents a Gaussian distribution.

[0039] Furthermore, the step 3 specifically includes:

[0040] Step 3.1, joint distribution p(z 1：k , r k ) can be expressed recursively as:

[0041]

[0042] Among them, p(z k |z 1：(k-1) , r k ) is the measured posterior prediction probability density in step 2.4; p(r k |r k-1 ) is the run length r in step 1.5 k-1 The transition probability p(z 1：(k-1) , r k-1 ) is the measured data z at the previous k-1 moments 1：(k-1) , run length r k-1 The joint distribution of

[0043] Step 3.2: Replace the run length transition probability p(r k |r k-1 ) into the joint distribution p(z 1：k , r k ) in the expression of 1：k , r k ) is decomposed into the joint probability p(z 1：k , r k = 1) and the joint probability p(z 1：k , r k =r k-1 +1), as follows:

[0044]

[0045] Furthermore, the step 4 specifically includes:

[0046] Step 4.1, based on the joint distribution p(z 1：k , r k ), calculate the measurement data z 1：k The marginal likelihood of is:

[0047]

[0048] Step 4.2: According to Bayes’ theorem, based on the marginal likelihood of the above measurement data and the measurement data z at the previous k moments 1：k , calculate the run length r k The posterior distribution of :

[0049]

[0050] Furthermore, the step 5 specifically includes:

[0051] Step 5.1, through the measurement data z at the current k moment k , calculate the run length as r k The estimated mean of the observable state at Unobservable state estimate mean and the corresponding estimated covariance

[0052] in and is the gain matrix;

[0053] Step 5.2, calculate the run length r k The target state vector The estimated mean of and the estimated variance

[0054]

[0055] where h = [0 1];

[0056] Step 5.3, for each run length r k , the target state vector in step 5.2 is Add to the end of the induction point set, and delete the first induction point in the induction point set, so as to complete the update of the induction point set;

[0057] Step 5.4, pruning: Given a threshold τ, prune the run length r in step 4.2 k The posterior distribution p(r k |z 1：k) Sort from large to small, retaining the first τ run lengths r k The posterior distribution of

[0058] Step 5.5: The state estimation mean μ(x) at time k can be obtained by weighted fusion. k|k ) and estimated covariance

[0059]

[0060] The state estimate mean μ(x k|k ) contains the required position estimation of the maneuvering target in the x direction or y direction, thereby achieving tracking of the maneuvering target; is the weight corresponding to the k-th moment when the running length is r, and the expression is as follows:

[0061]

[0062] p(r k |z 1：k ) is the run length r after pruning in step 5.4 k The posterior distribution of and After pruning in step 5.4, the running length at time k is r k The corresponding state estimated mean and estimated covariance;

[0063] Determine whether the measurement at time k is the last measurement of the sensor: if so, end the tracking; otherwise, jump to step 2.

[0064] A terminal device comprises a processor, a memory and a computer program stored in the memory; when the processor is executed by a computer, the mobile target tracking method based on a non-stationary Gaussian process is implemented.

[0065] A computer-readable storage medium stores a computer program; when the computer program is executed by a processor, the maneuvering target tracking method based on a non-stationary Gaussian process is implemented.

[0066] Compared with the prior art, the present invention has the following technical features:

[0067] The mobile target tracking method provided by the present invention uses a non-stationary Gaussian process to model the target motion pattern and recursively update the target state and model parameters through online learning, so as to realize state estimation and change point detection at the same time; the method is based on Bayesian online change point detection and uses the run length to divide the time series data into non-overlapping observation sets, that is, the run length is used to model the change point position, so as to determine the partition. The non-stationary process in the mobile target tracking scenario is converted into a partitioned stationary process, and for each partitioned stationary process, a Gaussian process is used to model the time series stationary process of the partition, thereby improving the robustness of mobile target tracking; through the simulation scenes of six groups of Benchmark aerial targets in Matlab, it is shown that the method of the present invention has the advantage of robust tracking effect compared with the existing recursive Gaussian process motion tracker (RGP) method. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 It is a flowchart of a method for tracking a maneuvering target based on a non-stationary Gaussian process of the present invention;

[0069] Figure 2-7 They are target motion trajectory diagrams for scenes S1-S6 respectively;

[0070] Figure 8-13 These are the target position RMSE curves for scenes S1-S6 respectively. DETAILED DESCRIPTION

[0071] Aiming at the situation that it is difficult to accurately model the target motion mode with a model in the scenario of mobile single target tracking, the present invention proposes a mobile target tracking method based on a non-stationary Gaussian process. The target trajectory is divided into non-intersecting stationary processes through Bayesian online change point detection. In each stationary process, a data-driven Gaussian process is proposed to model the target behavior, and the target state and parameters related to the Gaussian process are updated online by jointly estimating the augmented state vector.

[0072] Referring to the attached figure, a method for tracking a maneuvering target based on a non-stationary Gaussian process is provided. The subject of the method obtains the measurement data of the maneuvering target in real time through a sensor, and performs the following tracking process to obtain a more accurate estimated position of the target. It is assumed that the coupling between the two-dimensional coordinates of the target position is very weak, and the coordinates of each dimension are related within a certain continuous time. The specific steps are as follows:

[0073] Step 1: The subject constructs an online learning target motion model through Gaussian process to describe the possible motion mode of the maneuvering target, constructs a measurement model to represent the measurement data obtained by the sensor, and constructs a change point model to represent the time when the target motion mode changes; the transition probability of the change point model is given and the parameters to be estimated at the current moment are determined; the details are as follows:

[0074] Step 1.1, introduce a non-negative discrete integer variable rk ∈[1, k] is used as a change point model to represent the number of moments from the last change point of the measured data sequence to time k, that is, the running length at time k; in this scheme, when the target motion mode changes, let r k =1, at this time, the measured data sequence produces a change point; when the target motion mode does not change, r k =r k-1 +1, the measurement data sequence continues to grow.

[0075] Step 1.2: The duration between adjacent change points in the measured data sequence is a run length. In each run length, the target motion model can be uniformly expressed as f r (·) is the state transition function, Among them, k represents the current moment, u and u′ represent the Gaussian process covariance kernel function Any two different inputs; the input may be, for example, time, state at the previous moment, and other parameters, which are specifically determined according to the design of the target motion model; u k Represents the input at the current k moment; is the kth time in a run of length r k The target state variable when the running length is r k The estimated position in the x direction or the estimated position in the y direction; The running length at time k is r k The Gaussian process at time , its mean function is 0, and the covariance kernel function is The covariance kernel function is usually a square exponential kernel (SE) function, and the specific expression is:

[0076]

[0077] in and is the covariance kernel function The run length at time k is r k The hyperparameters when ||uu′|| 2 is the square of the Euclidean distance between input points u and u′.

[0078] Step 1.3, in each running length, assuming that the original measurement data obtained by the sensor can obtain the target position measurement after preprocessing such as coordinate transformation and registration, then the target measurement model can be expressed as Among them, z k is the measurement data at time k, indicating the position measurement data in the x direction or the position measurement data in the y direction; is the measurement noise of the sensor, which is zero-mean Gaussian distribution noise. is the measurement noise variance.

[0079] Step 1.4, in order to reduce the computational complexity, for a certain run length r k , based on the motion model representation in step 1.2, at time k-1, the induced point set is introduced in,[·] T The T in the upper right corner represents vector or matrix transpose; represents the dth induction point, where d represents the number of induction point sets; represents the input set of induced points at time k-1, represents the input of the d-th induction point at time k-1; the initial induction point set is the measured data at d moments before the start of the run length; since the target will maneuver during the motion, the hyperparameter vector Unknown, it is necessary to learn hyperparameters while estimating the target state; among them, represents the hyperparameter of the Gaussian process covariance kernel function at time k-1, σ k-1 represents the measurement noise parameter at time k-1. Therefore, the target state induction point set at time k-1 is Hyperparameter vector and the running length for k time is r k The predicted status at As an augmented state vector ψ k , construct the state space model as The specific meanings are as follows:

[0080]

[0081] Among them, ψ k is the augmented state vector, the subscript k-1 represents the time, I N and I3 represent N×N and 3×3 identity matrices respectively, O represents an all-zero vector or matrix, is the input at time k; represents the cross covariance matrix between the predicted state at time k and the induced point set at time k-1, represents the variance of the predicted state at time k, represents the autocovariance between the induced point sets at time k-1, yes The transpose of .

[0082] Step 1.5, determine the run length transition probability p(r k |r k-1 ).

[0083] The transition probability of the change point model is set to: run length r k-1 to r k The transition probability is p(rk |r k-1 ), when r k =1, p(r k |r k-1 )=H(r k-1 +1); otherwise r k =r k-1 +1, p(r k |r k-1 )=1-H(r k-1 +1); H(·) represents the hazard function of the change point model.

[0084] Step 1.6, select (x k , r k ) as the parameter to be estimated at time k, x k is the target state quantity to be estimated at the current k moment, indicating the position of the target in the x direction or the y direction; r k Indicates the running length at the current k moment.

[0085] Step 2: The subject obtains the target’s measurement data through the sensor; the measurement data of the previous d moments is used as the length r of each run. k The initial training data {z k-d+1 , …, z k}(k=d) or {z k-d , …, z k-1}(k>d) for online learning; where z k Represents the measurement data at time k; the initial hyperparameters can be obtained through Gaussian process regression training for subsequent recursive estimation.

[0086] When k>d, the subject performs all runs of length r at each moment according to the target motion model defined in step 1. k Calculate the measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , r k ), the posterior prediction density is calculated with a run length of r k and the measurement data z at the previous k-1 moments 1：(k-1) The conditions are as follows:

[0087] Step 2.1, calculate the run length r from the state space model in step 1.4 k The augmented state vector ψ k The predicted mean and the corresponding predicted covariance Since the state-space equation in step 1.4 is given by the hyperparameter It is linear only when , so it is necessary to use unscented transformation to select the hyperparameter vector The mean s sigma points; given the i-th sigma point Its weight is The predicted state mean corresponding to this sigma point and the corresponding covariance matrix as follows:

[0088]

[0089] in, Respectively represent A in the state space model k , The values ​​of the hyperparameters involved are and They represent the estimated mean and estimated covariance of the state of the induced point set at time k-1 respectively; and denote the estimated mean and estimated covariance of the hyperparameter vector respectively; is the estimated cross covariance between the hyperparameter vector and the set of induced points, for The transpose of ; O represents a zero matrix or a zero vector, as long as the dimensions match; the superscript -1 indicates matrix inversion, and the superscript T indicates matrix transpose.

[0090] Therefore, the augmented state vector ψ can be calculated through weighted fusion k The predicted mean and the predicted covariance for:

[0091]

[0092] Step 2.2, for each run length r k , the augmented state vector ψ k Decompose into observable vectors and the unobservable vector where σ k represents the measurement noise parameter at time k, represents the hyperparameter of the Gaussian process covariance kernel function at time k, represents the set of induced points at time k, then the state vector ψ k The predicted mean and the predicted covariance The specific decomposition form is as follows:

[0093]

[0094] Among them, the lower right corner scripts k and k-1 represent the time; the superscript ~ represents the predicted value; represents the unobservable vector g k-1The estimated mean of Observable vector o k The predicted mean of , which includes the measurement noise parameter σ k-1 The estimated mean of and the run length is r k The predicted status at The predicted mean represents the unobservable vector g k-1 The estimated covariance of Represents the observable vector o k The prediction covariance of , which includes the measurement noise parameter σ k-1 The estimated variance of Prediction Status The prediction variance and the measurement noise parameter σ k and predicted status The cross covariance between and for The transpose of represents the unobservable variable g k and observable variable o k The cross covariance between for The transpose of .

[0095] Step 2.3, using the measurement model in step 1.3 and the decomposition result in step 2.2, calculate the run length r k The measured prediction mean at and the corresponding predicted covariance

[0096]

[0097] Step 2.4, for all runs of length r k , based on the measurement data z at the previous k-1 moments 1：(k-1) Calculate the measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , r k ):

[0098]

[0099] in, and are the measured prediction mean and the corresponding prediction covariance in step 2.3, respectively. Represents a Gaussian distribution.

[0100] Step 3: The subject determines the measurement data z at the previous k moments according to the change point model in step 1. 1：k , run length r k The joint distribution p(z 1：k , r k ), as follows:

[0101] Step 3.1, joint distribution p(z 1：k , r k ) can be expressed recursively as:

[0102]

[0103] Among them, p(z k |z 1：(k-1) , r k ) is the measured posterior prediction probability density in step 2.4; p(r k |r k-1 ) is the run length r in step 1.5 k-1 The transition probability p(z 1：(k-1) , r k-1 ) is the measured data z at the previous k-1 moments 1：(k-1) , run length r k-1 The joint distribution of .

[0104] Step 3.2: Replace the run length transition probability p(r k |r k-1 ) into the joint distribution p(z 1：k , r k ) in the expression of 1：k , r k ) is decomposed into the joint probability p(z 1：k , r k = 1) and the joint probability p(z 1：k , r k =r k-1 +1), as follows:

[0105] p(z 1：k , r k =r k-1 +1)

[0106] =p(z k |z 1：(k-1) , r k )p(z 1：(k-1) , r k-1 )(1-H(r k-1 +1))

[0107]

[0108] Step 4: The subject calculates the joint distribution p(z 1：k , r k ), determine the running length r of the parameter to be estimated at time k k The posterior distribution of is as follows:

[0109] Step 4.1, based on the joint distribution p(z 1：k , r k ), calculate the measurement data z 1：k The marginal likelihood of is:

[0110]

[0111] Step 4.2: According to Bayes’ theorem, based on the marginal likelihood of the above measurement data and the measurement data z at the previous k moments 1：k , calculate the run length r k The posterior distribution of :

[0112]

[0113] Step 5: The subject calculates the target state x of the parameter to be estimated k The posterior distribution of where μ(x k|k ) and ∑(x k|k ) are the state estimation mean and state estimation variance of the target at time k, respectively, as follows:

[0114] Step 5.1, through the measurement data z at the current k moment k , calculate the run length r k The estimated mean of the observable state at Unobservable state estimate mean and the corresponding estimated covariance

[0115] in and is the gain matrix.

[0116] Step 5.2, calculate the run length r k The target state vector The estimated mean of and the estimated variance

[0117]

[0118] Where h = [0 1].

[0119] Step 5.3, for each run length r k, the target state vector in step 5.2 is Add it to the end of the induction point set and delete the first induction point in the induction point set to complete the update of the induction point set.

[0120] Step 5.4, pruning: Given a threshold τ, prune the run length r in step 4.2 k The posterior distribution p(r k |z 1：k ) Sort from large to small, retaining the first τ run lengths r k To reduce the computational complexity.

[0121] Step 5.5: The state estimation mean μ(x) at time k can be obtained by weighted fusion. k|k ) and estimated covariance

[0122]

[0123] The state estimate mean μ(x k|k ) contains the required position estimation of the maneuvering target in the x direction or y direction, thereby achieving tracking of the maneuvering target; is the weight corresponding to the k-th moment when the running length is r, and the expression is as follows:

[0124]

[0125] p(r k |z 1：k ) is the run length r after pruning in step 5.4 k The posterior distribution of ; similarly, and After pruning in step 5.4, the running length at time k is r k The corresponding state estimated mean and estimated covariance.

[0126] Determine whether the measurement at time k is the last measurement of the sensor: if so, end the tracking; otherwise, jump to step 2.

[0127] Example:

[0128] In this specific example, six different scenarios are given to simulate the present invention, as follows:

[0129] 1. Scene Design

[0130] Scenario 1 (S1): The first target trajectory represents a large aircraft, such as a military cargo aircraft. The target initially flies at a constant speed of 290 m / s in the k∈[1,60) time period; at time k∈[60,80) and k∈[110,130), the target makes turns with accelerations of 2g and 3g (g is the acceleration of gravity) respectively; the target continues to fly at a constant speed in the other time periods.

[0131] Scenario 2 (S2): The second target trajectory represents a small and flexible aircraft, such as a jet or other similar high-performance commercial aircraft. The target initially flies at a constant speed of 305 m / s for a period of time k∈[1,30); at time k∈[30,50), the target makes a 90° turn with an acceleration of 2.5 g (g is the acceleration of gravity); at time k∈[100,115), the target makes a turn with an acceleration of 4 g (g is the acceleration of gravity); the target continues to fly at a constant speed for the rest of the time.

[0132] Scenario 3 (S3): The third target represents a medium-sized bomber flying at high speed with good maneuverability. When time k∈[1,30], the target maintains a straight and horizontal flight at an initial constant speed of 457m / s; when time k∈[30,40), the target makes a 45° turn with an acceleration of 4g (g is the acceleration of gravity); when time k∈[70,110), the target makes a 90° turn with an acceleration of 4g (g is the acceleration of gravity), the heading changes to a horizontal straight line direction, and the speed drops to 274m / s; in other time periods, the target continues to fly at a constant speed on the new heading.

[0133] Scenario 4 (S4): The fourth target also represents a medium-sized bomber flying at high speed with good maneuverability. The target flies at an initial constant speed of 251 m / s within the time k∈[1,30); when time k∈[30,40), the target makes a 45° turn with an acceleration of 4g (g is the acceleration of gravity); when time k∈[70,80), the target makes a turn with an acceleration of 6g (g is the acceleration of gravity); in other time periods, the target maintains horizontal straight-line motion at a constant speed.

[0134] Scenario 5 (S5): The fifth target represents a fighter jet. The target flies at a constant speed at the beginning; then at time k∈[30, 40), k∈[60, 70), k∈[115, 170), the target makes turns with accelerations of 5g, 7g, and 6g (g is the acceleration of gravity) respectively; in other time periods, the target flies at a constant speed.

[0135] Scenario 6 (S6): The sixth target also represents a fighter jet. The target maintains a constant speed of 426 m / s for 30 seconds. Then, at time k∈[30, 40), k∈[70, 85), k∈[115, 125), k∈[125, 155), the target makes turns at accelerations of 7g, 6g, 6g, and 7g (g is the acceleration of gravity). In other time periods, the target maintains a constant speed and heading.

[0136] 2. Method-related configuration

[0137] In order to solve the influence of unknown target maneuvering mode on filtering effect in existing maneuvering target tracking methods, the present invention proposes a maneuvering target tracking method based on non-stationary Gaussian process. The maneuvering target here includes moving objects such as ships, vehicles, and aircraft that can be detected by sensors (such as radar).

[0138] In this embodiment:

[0139] When k=d, or k>d and a change point occurs to generate a new run length, the initial value of the distribution of the relevant variables is given as follows:

[0140] Given an initial value r for the run length k =1; given that the number of induction point sets is d=12, the initial training data for each run length is {z k-d+1 , …, z k}(k=d) or {z k-d , …, z k-1}(k>d), the initial hyperparameters can be obtained through Gaussian process regression training Initial distribution of hyperparameters Mean and covariance The specific form is as follows:

[0141]

[0142] Given an initial induction point set is the measured value at d moments, and its specific value is as follows:

[0143]

[0144] Then, the initial distribution of the initial inducing point set is is the covariance, the mean The values ​​are as follows:

[0145]

[0146] Given the induction point set and the initial cross covariance of the hyperparameters is a zero matrix; given the initial joint distribution of the run length and its corresponding measurement data p(z 1：k , r k )=1; given risk function H=1 / 10.

[0147] 3. Assessment of mobile target tracking:

[0148] The following comparative evaluation algorithm uses the Recursive Gaussian Process motion tracker (RGP).

[0149] In the method of the present invention, the filtering method based on non-stationary Gaussian process is called RGP-BOCPD. The evaluation index is the root mean square error RMSE of the target radial distance. The smaller the RMSE, the better the filtering effect.

[0150] The simulation results and analysis are as follows:

[0151] Referring to the six sets of Benchmark data in MATLAB, six single-target tracking scenarios with variable acceleration maneuvers are considered. The total simulation step length of each tracking scenario is 187s, the sensor sampling period is T=1s, and the measurement noise variance is The number of Gaussian process induced point sets of the proposed algorithm (RGP-BOCPD) and the comparative algorithm (RGP) is d = 12. The number of Monte Carlo simulations is 30.

[0152] Table 1 Average position RMSE statistics of six simulation scenarios under different methods

[0153] Scenario RGP-BOCPD RGP Measurement S1 37.5992 41.6784 43.7770 S2 35.2562 537.8034 44.2693 S3 34.6308 42.2197 44.1759 S4 34.7259 281.1420 44.3870 S5 38.8173 63.6304 44.8457 S6 40.4823 1018.9 44.3301

[0154] Figure 2-7 The following are the tracking scene diagrams for scenes 1-6. Figure 8-13In the figure, they are the target position RMSE curves. In the target position RMSE curve, the red color is the algorithm RGP-BOCPD of this scheme, and the blue color is the comparison algorithm RGP. It can be seen from the figure that the RMSE of the algorithm of this scheme will increase when the target is maneuvering, the tracking accuracy will be reduced to a certain extent, and the subsequent RMSE will drop rapidly. This is because when the target maneuver is detected, the algorithm of this scheme needs to reinitialize the motion mode of the learning target based on less measurement data. As time goes by, the measurement data continues to increase, and the model parameters can be continuously corrected. Table 1 shows the average position RMSE of different algorithms under six groups of simulation scenarios. It can be seen from the position RMSE that the RMSE value of the algorithm of this scheme is smaller than that of the comparison algorithm RGP in the six simulation scenarios, indicating that the filtering effect of the algorithm RGP-BOCPD of this scheme is better than that of the comparison algorithm. The RMSE values ​​of the comparison algorithm in scenes S2, S4, S5, and S6 are larger than the measured RMSE values. This is because in scenes with strong maneuverability, the algorithm system parameters are updated slowly and cannot describe the target motion pattern after the maneuver in time, resulting in a decrease in target tracking performance or even loss of tracking. The RMSE values ​​of the algorithm in this scheme are smaller than the measured RMSE values ​​in all scenes, indicating that the algorithm in this scheme can effectively filter in all maneuvering target tracking scenes. This is because we combine Bayesian online change point detection. When the target is detected to be maneuvering, the Gaussian process is used to relearn the target's motion pattern, which can timely and effectively characterize the target motion, thereby ensuring the effectiveness and robustness of the filtering of the algorithm in this scheme.

Claims

1. A maneuvering target tracking method based on a non-stationary Gaussian process, characterized in that: The following steps are involved: Step 1: The subject constructs an online learning target motion model through Gaussian process to describe the possible motion of the maneuvering target, constructs a measurement model to represent the measurement data obtained by the sensor, and constructs a change point model to represent the time when the target motion mode changes; Given the transition probability of the change point model, determine the parameters to be estimated at the current moment; Step 2: The subject obtains the measurement data of the target through the sensor, and uses the measurement data of the previous multiple moments as the initial training data of each running length for online learning. The initial hyperparameters can be obtained through Gaussian process regression training for subsequent recursive estimation; Step 3, the subject determines the joint distribution of the measurement data and the run length at the previous k moments according to the change point model; Step 4, the subject determines the posterior distribution of the running length of the parameter to be estimated at time k according to the joint distribution; Step 5: The subject calculates the posterior distribution of the target state of the parameter to be estimated, which contains the state estimation mean and state estimation variance of the target at time k; The state estimation mean contains the required position estimation of the maneuvering target, thereby achieving tracking of the maneuvering target.

2. The method for tracking a maneuvering target based on a non-stationary Gaussian process according to claim 1, characterized in that: The step 1 specifically includes: Step 1.1, introduce a non-negative discrete integer variable r k ∈[1, k] as the change point model to represent the number of moments from the last change point of the measured data sequence to time k, that is, the running length at time k; when the target motion mode changes, let r k =1, at this time, the measured data sequence produces a change point; when the target motion mode does not change, r k =r k-1 +1, the measurement data sequence continues to grow; Step 1.2: The duration between adjacent change points in the measured data sequence is a run length. In each run length, the target motion model is expressed as f r (·) is the state transition function, Among them, k represents the current moment, u and u′ represent the Gaussian process covariance kernel function Any two different inputs, u k Represents the input at the current k moment; is the kth time in the run length r k The target state variable when the running length is r k The estimated position in the x direction or the estimated position in the y direction; The running length at time k is r k The Gaussian process at time , its mean function is 0, and the covariance kernel function is The covariance kernel function is usually a squared exponential kernel (SE) function; Step 1.3: In each running length, the original measurement data obtained by the sensor is preprocessed by coordinate transformation and registration to obtain the position measurement of the target. Then the target measurement model can be expressed as Among them, z k is the measurement data at time k, indicating the position measurement data in the x direction or the position measurement data in the y direction; is the measurement noise of the sensor, which is zero-mean Gaussian distribution noise. is the measurement noise variance; Step 1.4, for a certain run length r k , based on the motion model representation in step 1.2, at time k-1, the induction point set is introduced in,[·] T The T in the upper right corner represents vector or matrix transpose; represents the dth induction point, where d represents the number of induction point sets; represents the input set of induced point sets at time k-1, represents the input of the d-th induction point at time k-1; Initial induction point set is the measured data at d moments before the start of the run length; since the target will maneuver during the motion, the hyperparameter vector Unknown, it is necessary to learn hyperparameters while estimating the target state; among them, represents the hyperparameter of the Gaussian process covariance kernel function at time k-1, σ k-1 represents the measurement noise parameter at time k-1; the target state induction point set at time k-1 Hyperparameter vector and the running length for k time is r k The predicted status at As an augmented state vector ψ k , construct the state space model as The specific meanings are as follows: Among them, ψ k is the augmented state vector, the subscript k-1 represents the time, I N and I3 represent N×N and 3×3 identity matrices respectively, O represents an all-zero vector or matrix, is the input at time k; represents the cross covariance matrix between the predicted state at time k and the induced point set at time k-1, represents the variance of the predicted state at time k, represents the autocovariance between the induced point sets at time k-1, yes The transpose of Step 1.5, the transition probability of the change point model is set to: run length r k-1 to r k The transition probability is p(r k |r k-1 ), when r k =1, p(r k |r k-1 )=1-H(r k-1 +1); otherwise r k =r k-1 +1, p(r k |r k-1 )=1-H(r k-1 +1); H(·) represents the risk function of the change point model; Step 1.6, select (x k , r k ) as the parameter to be estimated at time k, x k is the target state quantity to be estimated at the current k moment, indicating the position of the target in the x direction or the y direction; r k Indicates the running length at the current k moment.

3. The method for tracking a maneuvering target based on a non-stationary Gaussian process according to claim 2, characterized in that: The step 2 specifically includes: When k>d, the subject performs all runs of length r at each moment according to the target motion model defined in step 1. k Calculate the measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , r k ), the posterior prediction density is calculated with a run length of r k and the measurement data z at the previous k-1 moments 1：(k-1) The conditions are as follows: Step 2.1, calculate the run length r from the state space model in step 1.4 k The augmented state vector ψ k The predicted mean and the corresponding predicted covariance Since the state-space equation in step 1.4 is given by the hyperparameter It is linear only when , so it is necessary to use unscented transformation to select the hyperparameter vector The mean s sigma points; given the i-th sigma point Its weight is The predicted state mean corresponding to this sigma point and the corresponding covariance matrix as follows: in, Respectively represent A in the state space model k , The values ​​of the hyperparameters involved are and They represent the estimated mean and estimated covariance of the state of the induced point set at time k-1 respectively; and denote the estimated mean and estimated covariance of the hyperparameter vector respectively; is the estimated cross covariance between the hyperparameter vector and the set of induced points, for The transpose of ; O represents a zero matrix or a zero vector; the superscript -1 of the parameter represents the matrix inversion, and the superscript T represents the matrix transpose; After weighted fusion, the augmented state vector ψ can be calculated k The predicted mean and the predicted covariance for: Step 2.2, for each run length r k , the augmented state vector ψ k Decompose into observable vectors and the unobservable vector where σ k represents the measurement noise parameter at time k, represents the hyperparameter of the Gaussian process covariance kernel function at time k, represents the set of induced points at time k, then the state vector ψ k The predicted mean and the predicted covariance The specific decomposition form is as follows: Among them, the lower right corner scripts k and k-1 represent the time; the superscript ~ represents the predicted value; represents the unobservable vector g k-1 The estimated mean of Observable vector o k The predicted mean of , which includes the measurement noise parameter σ k-1 The estimated mean of and the run length is r k The predicted status at The predicted mean represents the unobservable vector g k-1 The estimated covariance of Represents the observable vector o k The prediction covariance of , which includes the measurement noise parameter σ k-1 The estimated variance of Prediction Status The prediction variance and the measurement noise parameter σ k and predicted status The cross covariance between and for The transpose of represents the unobservable variable g k and observable variable o k The cross covariance between for The transpose of Step 2.3, using the measurement model in step 1.3 and the decomposition result in step 2.2, calculate the run length r k The measured prediction mean at and the corresponding predicted covariance Step 2.4, for all runs of length r k , based on the measurement data z at the previous k-1 moments 1：(k-1) Calculate the measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , r k ): in, and are the measured prediction mean and the corresponding prediction covariance in step 2.3, respectively. Represents a Gaussian distribution.

4. The method for tracking a maneuvering target based on a non-stationary Gaussian process according to claim 3, characterized in that: The step 3 specifically includes: Step 3.1, joint distribution p(z 1：k , r k ) can be expressed recursively as: Among them, p(z k |z 1：(k-1) , r k ) is the measured posterior prediction probability density in step 2.4; p(r k |r k-1 ) is the run length r in step 1.5 k-1 The transition probability p(z 1：(k-1) , r k-1 ) is the measured data z at the previous k-1 moments 1：(k-1) , run length r k-1 The joint distribution of Step 3.2: Replace the run length transition probability p(r k |r k-1 ) into the joint distribution p(z 1：k , r k ) in the expression of 1：k , r k ) is decomposed into the joint probability p(z 1：k , r k = 1) and the joint probability p(z 1：k , r k =r k-1 +1), as follows:

5. The method for tracking a maneuvering target based on a non-stationary Gaussian process according to claim 4, characterized in that: The step 4 specifically includes: Step 4.1, based on the joint distribution p(z 1：k , r k ), calculate the measurement data z 1：k The marginal likelihood of is: Step 4.2: According to Bayes’ theorem, based on the marginal likelihood of the above measurement data and the measurement data z at the previous k moments 1：k , calculate the run length r k The posterior distribution of :

6. The method for tracking a maneuvering target based on a non-stationary Gaussian process according to claim 5, characterized in that: The step 5 specifically includes: Step 5.1, through the measurement data z at the current k moment k , calculate the run length as r k The estimated mean of the observable state at Unobservable state estimate mean and the corresponding estimated covariance in and is the gain matrix; Step 5.2, calculate the run length r k The target state vector The estimated mean of and the estimated variance where h = [0 1]; Step 5.3, for each run length r k , the target state vector in step 5.2 is Add to the end of the induction point set, and delete the first induction point in the induction point set, so as to complete the update of the induction point set; Step 5.4, pruning: Given a threshold τ, prune the run length r in step 4.2 k The posterior distribution p(r k |z 1：k ) Sort from large to small, retaining the first τ run lengths r k The posterior distribution of ; Step 5.5: The state estimation mean μ(x) at time k can be obtained by weighted fusion. k|k ) and the estimated covariance ∑(x k|k ): The state estimate mean μ(x k|k ) contains the required position estimation of the maneuvering target in the x direction or y direction, thereby achieving tracking of the maneuvering target; is the weight corresponding to the k-th moment when the running length is r, and the expression is as follows: p(r k |z 1：k ) is the run length r after pruning in step 5.4 k The posterior distribution of ; and After pruning in step 5.4, the running length at time k is r k The corresponding state estimated mean and estimated covariance; Determine whether the measurement at time k is the last measurement of the sensor: if so, end the tracking; otherwise, jump to step 2.

7. A terminal device comprising a processor, a memory and a computer program stored in the memory; characterized in that: When the processor is executed by a computer, the maneuvering target tracking method based on a non-stationary Gaussian process according to any one of claims 1 to 6 is implemented.

8. A computer-readable storage medium, wherein a computer program is stored in the medium; characterized in that: When the computer program is executed by a processor, the method for tracking a maneuvering target based on a non-stationary Gaussian process according to any one of claims 1 to 6 is implemented.

Images