Non-stationary gaussian process based maneuvering target tracking method

Abstract

The application discloses a kind of based on non-stationary Gaussian process maneuvering target tracking method, target motion mode is modeled using non-stationary Gaussian process and on-line learning recursive updates target state and model parameter, realize state estimation and change point detection simultaneously;This method is based on Bayesian on-line change point detection using run length to divide time series data into non-overlapping observation set, i.e. using run length to model change point position, so as to determine partition;Non-stationary process under maneuvering target tracking scene is converted into partition stationary process, for each partition stationary process, Gaussian process is used to model the time series stationary process of partition, improve the robustness of maneuvering target tracking;Simulation shows that the method of the present application has the advantage of robust tracking effect compared with the method of existing recursive Gaussian process motion tracker.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of single-sensor single-maneuvering target tracking subsystems in 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 unknown target maneuvering mode. BACKGROUND

[0002] Target tracking is a basic task in radar, sonar, navigation and other sensor applications, that is, in a continuous data sequence, the position relationship of the target to be tracked is established, and the complete motion trajectory of the target is obtained. 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 the extended Kalman filter (EKF), the unscented Kalman filter (UKF) and the interacting multiple model (IMM) filter. However, the above methods all depend on the underlying target dynamic model, that is, a pre-set target dynamic model needs to match the actual motion of the target during maneuvering. Therefore, when the target experiences unknown or mixed maneuvering behavior, there may be significant uncertainty in the motion model, so that the target state cannot be approximated by a proper mathematical model, resulting in a decline in the tracking performance based on the Bayesian filter.

[0004] In recent years, target tracking technology based on Gaussian process online model learning, which learns the motion pattern of the target from data, has become a new research hotspot with 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 the position and high-order position derivatives. The Gaussian Process Motion Tracker (GPMT) includes a training and learning phase, which models an infinite number of target trajectories, but the learning phase in GPMT is not recursive. Based on GPMT, the Recursive Gaussian Process Motion Tracker (RGPMT) performs online model-free point target tracking, and the state estimation and learning of the Gaussian process hyperparameters are performed recursively.

[0005] However, the Recursive Gaussian Process Motion Tracker is mainly used to solve the tracking problem of stationary time series processes, and when the target changes from a steady state, that is, a large maneuvering occurs, the tracking performance of the tracker may decline. SUMMARY

[0006] In order to solve the problem that filtering precision influences filtering effect and the solving process is complex in the existing motor target tracking method, the motor target tracking method based on non-stationary Gaussian process is provided, the state estimation problem is effectively handled, and the motor target tracking performance is improved.

[0007] In order to realize the above task, the technical scheme is adopted in the present application:

[0008] The motor target tracking method based on non-stationary Gaussian process comprises the following steps:

[0009] Step 1, the subject constructs the target motion model of online learning by Gaussian process to describe the possible motion mode of the motor target, constructs the measurement model to represent the measurement data obtained by the sensor, and constructs the 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 parameter to be estimated at the current time is determined;

[0010] Step 2, the subject obtains the measurement data of the target through the sensor, takes the measurement data at the previous multiple times as the initial training data of each running length for online learning, and obtains the initial hyperparameter through Gaussian process regression training, which is used for subsequent recursive estimation;

[0011] Step 3, the subject determines the joint distribution of the measurement data at the previous k times and the running length 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 the k time 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 contains the state estimation mean and state estimation variance of the target at the k time; wherein the state estimation mean contains the position estimation of the motor target required, so as to realize the tracking of the motor target.

[0014] Further, the step 1 specifically comprises:

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

[0016] Step 1.2, the time length between adjacent change points in the measurement data sequence is a running length; in each running length, the target motion model is represented as f r (·) is a state transition function, where k represents the current time, u and u' represent any two different inputs of the Gaussian process covariance kernel function , u k represents the input at the current k time; is the target state variable at the k time with a running length of r k , represents the estimated position in the x direction or the estimated position in the y direction with a running length of r k ; is a Gaussian process with a running length of r k at the k time, whose mean function is 0 and whose 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 raw measurement data obtained by the sensor after coordinate conversion and registration preprocessing obtains the position measurement of the target, so the target measurement model can be expressed as where z k is the measurement data at the k time, representing the position measurement data in the x direction or the position measurement data in the y direction; is the measurement noise of the sensor, i.e., zero-mean Gaussian distributed noise, is the measurement noise variance;

[0018] Step 1.4, for a running length r k , from the motion model expression in step 1.2, at the k-1 time, an induced point set is introduced, where [·] T T in the upper right corner of [·] represents vector or matrix transposition; represents the dth induced point, and d represents the number of induced point sets; represents the input set of the induced point set at the k-1 time, represents the input of the dth induced point at the k-1 time;

[0019] The initial induced point set is the measurement data at the d time before the start of the running length; since the target motion process will occur, the hyperparameter vector is unknown, and needs to be learned while estimating the target state; where represents the hyperparameter of the Gaussian process covariance kernel function k-1 time, σ k-1 represents the measurement noise parameter at the k-1 time; the k-1 time target state induced point set is introduced, the hyperparameter vector and the predicted state of the running length r k at the k time 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 moment, 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 hazard 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 target position 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 time k>d, the subject performs all runs of length r at each time according to the target motion model defined in step 1. k Calculate measurement data z k The posterior prediction probability density p(z k |z 1：(k-1) , rk ), given the measurement data z k at time k-1, as follows: 1：(k-1)

[0026] Step 2.1, compute the predicted mean of the augmented state vector ψ k at time k-r k with the state space model in Step 1.4, as follows: and the corresponding predicted covariance Since the state space equation in Step 1.4 is linear given the hyperparameters , the hyperparameter vector is selected by using unscented transformation, as follows: sigma points of the mean s ; given the i th sigma point with weight , the predicted state mean and the corresponding covariance matrix corresponding to the sigma point are computed as follows:

[0027]

[0028] where denote the values of the hyperparameters involved in A k , in the state space model, respectively; and denote the state estimation mean and estimation covariance of the induced point set at time k-1, respectively; and denote the estimation mean and estimation covariance of the hyperparameter vector, respectively; is the estimation cross-covariance between the hyperparameter vector and the induced point set, is the transpose of ; O represents zero matrix or zero vector; the parameter superscript -1 denotes matrix inversion, and the superscript T denotes matrix transpose;

[0029] The predicted mean and the predicted covariance of the augmented state vector ψ k can be computed by weighted fusion as follows:

[0030]

[0031] Step 2.2, for each run length r k , decompose the augmented state vector ψ k into the observable vector 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 a run length of r k The predicted state 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 , predicted status The prediction variance of 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 variables 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 results 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 1：(k-1) at the previous k-1 time instants k : k 1：(k-1) |z k , r 1：k ):

[0037]

[0038] where and are the measurement prediction mean and the corresponding prediction covariance in step 2.3, respectively, denotes a Gaussian distribution.

[0039] Further, the step 3 specifically comprises:

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

[0041]

[0042] where p(z k |z 1：(k-1) , r k ) is the measurement posterior prediction probability density in step 2.4; p(r k |r k-1 ) is the transition probability of the run length r k-1 in step 1.5; p(z 1：(k-1) , r k-1 ) is the joint distribution of the measurement data z 1：(k-1) and the run length r k-1 at the previous k-1 time instants;

[0043] Step 3.2, the run length transition probability p(r k |r k-1 ) in step 1.5 is substituted into the expression of the joint distribution p(z 1：k , r k ) in step 3.1, and then the above joint distribution p(z 1：k , r k ) is decomposed into the joint probability p(z 1：k , r k = 1) when the change point occurs and the joint probability p(z 1：k , r k = r k-1 + 1) when the run length increases, as follows:

[0044]

[0045] Further, the step 4 specifically comprises:

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

[0047]

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

[0049]

[0050] Further, the step 5 specifically comprises:

[0051] Step 5.1, through the measurement data z k at the current k time, calculate the observable state estimation mean , the unobservable state estimation mean and the corresponding estimation covariance when the run length is r k :

[0052] wherein and are gain matrices;

[0053] Step 5.2, calculate the estimation mean and the estimation variance of the target state vector when the run length is r k :

[0054]

[0055] wherein h = [0 1];

[0056] Step 5.3, for each run length r k , add the target state vector in step 5.2 to the end of the inducing point set, and delete the first inducing point in the inducing point set, thereby completing the update of the inducing point set;

[0057] Step 5.4, pruning: given a threshold τ, delete the posterior distribution p(r k |z k ) of the run length r 1：k in step 4.2, and the corresponding target state vector) from large to small, and keep the first tau run length r k Posterior distribution of r

[0058] Step 5.5, by weighted fusion, the state estimation mean mu(x k|k ) and estimation covariance

[0059]

[0060] Wherein the state estimation mean mu(x k|k ) contains the required position estimation of the maneuvering target x direction or y direction, so as to realize the tracking of the maneuvering target; The weight corresponding to the run length r at k time, the expression is as follows:

[0061]

[0062] P(r k |z 1：k ) is the posterior distribution of the run length r k Pruned in step 5.4; And The state estimation mean and estimation covariance corresponding to the run length r k At k time after pruning in step 5.4;

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

[0064] A terminal device, comprising a processor, a memory and a computer program stored in the memory; the processor is computer executed to realize the maneuvering target tracking method based on the non-stationary Gaussian process.

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

[0066] Compared with the prior art, the present application has the following technical characteristics:

[0067] The application provides a motorized target tracking method, which models target motion mode by using a non-stationary Gaussian process, and recursively updates target state and model parameters on-line to realize simultaneous state estimation and change point detection; the method is based on Bayesian on-line change point detection, uses run length to divide time series data into non-overlapping observation sets, i.e., uses run length to model change point positions, so as to determine partitions. The non-stationary process in the motorized target tracking scene is converted into a partitioned stationary process, and for each partitioned stationary process, a Gaussian process is used to model the partitioned time series stationary process, so that the robustness of motorized target tracking is improved; through simulation scenes of six groups of Benchmark air targets in Matlab, it is shown that the method has the advantage of robust tracking effect compared with the method of the existing recursive Gaussian process motion tracker (RGP). BRIEF DESCRIPTION OF DRAWINGS

[0068] Figure 1 A flow chart of the motorized target tracking method based on the non-stationary Gaussian process of the application;

[0069] Figures 2-7 The target motion trajectory graphs of scenes S1-S6 respectively;

[0070] Figures 8-13 The target position RMSE curve graphs of scenes S1-S6 respectively. DETAILED DESCRIPTION

[0071] The application proposes a motorized target tracking method based on a non-stationary Gaussian process for the case that it is difficult to accurately model target motion mode by using a model in the motorized single target tracking scene, divides the target trajectory into disjoint stationary processes by Bayesian on-line change point detection, adopts a data-driven Gaussian process to model target behavior in each stationary process, and updates target state and parameters related to the Gaussian process on-line by joint estimation of an augmented state vector.

[0072] Referring to the drawings, a motorized target tracking method based on a non-stationary Gaussian process is provided, a subject applying the method acquires motorized target measurement data in real time by a sensor, and performs the following tracking process to obtain a more accurate target estimated position, assuming that the coupling between two-dimensional coordinates of the target position is very weak, and the coordinates in each dimension are related within a certain continuous time; the specific steps are as follows:

[0073] Step 1, the subject constructs an on-line learning target motion model by a Gaussian process to describe possible motion modes of the motorized target, constructs a measurement model to represent 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 time are determined; the details are as follows:

[0074] Step 1.1, a non-negative discrete integer variable r is introducedk ∈[1, k] as a change point model, represents the number of time experienced since the last change point of the measurement data sequence to time k, i.e., the running length at time k; in the present scheme, when the target motion pattern changes, let r k = 1, at this time the measurement data sequence produces a change point; when the target motion pattern does not change, r k = r k-1 + 1, then the measurement data sequence continues to grow.

[0075] Step 1.2, the length of time between adjacent change points in the measurement data sequence is a running length; in each running length, the target motion model can be uniformly represented as f r (·) is a state transition function, wherein k represents the current time, u and u' represent any two different inputs of the Gaussian process covariance kernel function ; the input may be, for example, time, state at the last time, etc. parameter, which is determined according to the design of the target motion model; u k represents the input at the current time k; is the target state variable at time k with a running length of r k , which represents the estimated position in the x direction or the estimated position in the y direction with a running length of r k . is a Gaussian process at time k with a running length of r k , whose mean function is 0 and the covariance kernel function is The covariance kernel function is usually a squared exponential kernel (SE) function, and the specific expression is:

[0076]

[0077] wherein and are the covariance kernel functions at time k with a running length of r k , and ||u-u'| 2 is the square of the Euclidean distance between the input points u and u'.

[0078] Step 1.3, in each running length, it is assumed that the original measurement data obtained by the sensor after coordinate conversion and registration and other pretreatments can obtain the position measurement of the target, so the target measurement model can be represented as wherein z k is the measurement data at time k, representing the position measurement data in the x direction or the position measurement data in the y direction; is the measurement noise of the sensor, i.e., zero-mean Gaussian distributed noise, is the measurement noise variance.

[0079] Step 1.4, to reduce the computational complexity, for a certain run length r k , the motion model in step 1.2 is expressed as, at time k-1, the set of inducing points where [·] T T in the top right corner of a matrix represents the transpose of a vector or a matrix; represents the dthinducing point, d represents the number of the set of inducing points; represents the input set of the set of inducing points at time k-1, represents the input of the dthinducing point at time k-1; the initial set of inducing points is the measurement data of the d time before the start of the run length; since the target motion process will occur maneuver, the hyperparameter vector is unknown, and needs to be learned while estimating the target state; where, represents the hyperparameter at time k-1 of the Gaussian process covariance kernel function k-1, and k-1 represents the measurement noise parameter at time k-1. Therefore, the set of target state inducing points at time k-1 is the hyperparameter vector and the predicted state of the run length r k at time k is as an augmented state vector ψ k , the state space model is constructed as The specific meaning is as follows:

[0080]

[0081] where ψ k is the augmented state vector, the subscript k-1 represents the time, I N and I3 represent the N×N, 3×3 identity matrix, O represents the 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 set of inducing points at time k-1, represents the variance of the predicted state at time k, represents the self-covariance between the set of inducing points at time k-1, is 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 as: the transition probability from run length r k-1 to r k 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 target position 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 time k>d, the subject performs all runs of length r at each time according to the target motion model defined in step 1. k Calculate 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 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 when , so it is necessary to use unscented transformation to select the hyperparameter vector the mean of s sigma points; given the i-th sigma point whose weight is the predicted state mean corresponding to the sigma point and the corresponding covariance matrix as follows:

[0088]

[0089] where, denote the A k , the values of the hyperparameters involved in and denote the state estimation mean and estimation covariance of the induced point set at time k-1, respectively; and denote the estimation mean and estimation covariance of the hyperparameter vector, respectively; is the estimation cross-covariance between the hyperparameter vector and the induced point set, is the transpose of ; O represents a zero matrix or a zero vector, which guarantees the dimension matching; the parameter superscript -1 represents the matrix inversion, and the superscript T represents the matrix transpose.

[0090] Therefore, the prediction mean k and the prediction covariance of the augmented state vector ψ can be calculated through the weighted fusion as follows:

[0091]

[0092] Step 2.2, for each running length r k , the augmented state vector ψ k is decomposed into the observable vector 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 induced point set at time k, then the prediction mean and the prediction covariance of the state vector ψ k have the following specific decomposition forms:

[0093]

[0094] where, the lower right subscripts k and k-1 represent the time; the superscript ~ represents the predicted value; represents the unobservable vector g k-1the estimated mean of the observable vector o k k-1 the predicted mean of and the predicted state at time r k the predicted mean of representing the unobservable vector g k-1 the estimated covariance of representing the observable vector o k the predicted covariance of k-1 the estimated variance of the predicted state the predicted variance of and the cross-covariance between the measurement noise parameter k and the predicted state and is the transpose of denotes the cross-covariance between the unobservable variable g k and the observable variable o k is the transpose of

[0095] Step 2.3, based on the decomposition results of Step 2.2, compute the measurement predicted mean k and the corresponding predicted covariance at time r

[0096]

[0097] Step 2.4, for all run lengths r k , compute the posterior predictive probability density p(z k |z 1：(k-1) , r k ) of the measurement data z k based on the measurement data z 1：(k-1) at the previous k - 1 times:

[0098]

[0099] where and are the measurement predicted mean and the corresponding predicted covariance in Step 2.3, respectively, denotes a Gaussian distribution.

[0100] ​​​​​​Step 3, the main body determines the measurement data z of the previous k time according to the change point model in step 1 1：k , the running length r k The joint distribution p(z 1：k , r k ) is as follows:

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

[0102]

[0103] Wherein, p(z k |z 1：(k-1) , r k ) is the measurement posterior predictive probability density in step 2.4; p(r k |r k-1 ) is the transition probability of running length r k-1 In step 1.5; p(z 1：(k-1) , r k-1 ) is the joint distribution of measurement data z 1：(k-1) , running length r k-1 At the previous k-1 time.

[0104] Step 3.2, the running length transition probability p(r k |r k-1 ) in step 1.5 is substituted into the expression of the joint distribution p(z 1：k , r k ) in step 3.1, and then the above joint distribution p(z 1：k , r k ) is decomposed into the joint probability p(z 1：k , r k =1) when the change point occurs and the joint probability p(z 1：k , r k =r k-1 +1) when the running length increases, 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 is distributed according to 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 :

[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, 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 estimated 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 kThe target state vector in step 5.2 is added to the end of the set of induced points, and the first induced point in the set of induced points is deleted, thereby completing the update of the set of induced points.

[0120] Step 5.4, pruning: given a threshold τ, the posterior distribution p(r k |z k ) of the run length r 1：k in step 4.2 is sorted in descending order, and the posterior distribution of the first τ run lengths r k is retained to reduce the computational complexity.

[0121] Step 5.5, the state estimation mean μ(x k|k ) and estimation covariance

[0122]

[0123] wherein the state estimation mean μ(x k|k ) contains the required position estimation of the maneuvering target x in the x direction or the y direction, thereby achieving tracking of the maneuvering target; is the weight corresponding to the run length r at time k, and the expression is as follows:

[0124]

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

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

[0127] Embodiment:

[0128] In a specific example of the present application, six different scenarios are given to simulate the present application, as follows:

[0129] I. Scene design

[0130] ​Scenario 1 (SI): The first target trajectory represents a large aircraft, such as a military cargo plane. The target flies at a constant speed of 290 m / s for the time period k e [1, 60); when time k e [60, 80), k e [110, 130], the target makes a turn with an acceleration of 2g, 3g (g is the acceleration of gravity) respectively; the target continues to fly at a constant speed for other time periods.

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

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

[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 for the time period k e [1, 30); when time k e [30, 40), the target makes a 45° turn with an acceleration of 4g (g is the acceleration of gravity); when time k e [70, 80), the target makes a turn with an acceleration of 6g (g is the acceleration of gravity); the target continues to fly horizontally at a constant speed for other time periods.

[0134] Scenario 5 (S5): The fifth target represents a fighter. The target initially flies at a constant speed; next, at time k e [30, 40), k e [60, 70), k e [115, 170), the target makes a turn with an acceleration of 5g, 7g, 6g (g is the acceleration of gravity) respectively; the target flies at a constant speed for other time periods.

[0135] Scenario 6 (S6): The 6th target also represents a fighter aircraft. The target keeps constant speed with initial speed 426 m / s for 30 s; then in the time period k∈[30, 40), k∈[70, 85), k∈[115, 125), k∈[125, 155), the target makes a turn with acceleration 7g, 6g, 6g, 7g (g is the gravity acceleration) respectively; in other time period, the target keeps constant speed and heading.

[0136] II. Method related configuration

[0137] In order to solve the influence of unknown target maneuver mode on filtering effect in existing maneuvering target tracking method, the application provides a maneuvering target tracking method based on non-stationary Gaussian process. The maneuvering target described herein includes ships, vehicles, aircrafts and other mobile objects that can be detected by sensors (such as radar).

[0138] In this embodiment:

[0139] When the time k=d, or k>d and a new running length is generated due to a change point, the initial value of the related variable distribution is given, and the specific process is as follows:

[0140] The initial value of the running length is given as r k =1; given that the number of induced point sets is d=12, the initial training data of each running length is {z k-d+1 ,…,z k}(k=d) or {z k-d ,…,z k-1}(k>d), and the initial hyperparameters The initial distribution of the hyperparameters The mean and the covariance The specific form is as follows:

[0141]

[0142] The initial induced point set is given as z The measurement values at d time points are as follows:

[0143]

[0144] The initial distribution of the initial induced point set is The covariance is The mean The values are as follows:

[0145]

[0146] The initial cross-covariance of the induced point set and the hyperparameters is given as is a zero matrix; given the initial joint distribution p(z 1：k , r k )=1; given risk function H=1 / 10.

[0147] 3. Assessment of tracking of mobile targets:

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

[0149] In the present method, the filtering method based on a non-stationary Gaussian process is called RGP-BOCPD. The evaluation metric is the root mean square error (RMSE) of the radial distance to the target. A smaller RMSE indicates a better 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 for 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 for both the proposed algorithm (RGP-BOCPD) and the comparative algorithm (RGP) is d = 12. The number of Monte Carlo simulations is 30.

[0152] Table 1 Statistics of average position RMSE of six simulation scenarios under different methods

[0153] Scenario RGP-BOCPD RGP Metrology 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] Figures 2-7 The following are the tracking scene diagrams for scenes 1-6 respectively. Figures 8-13The target position RMSE curve chart is shown in FIG. 6, in which red is the algorithm RGP-BOCPD of the present scheme, and blue is the comparative algorithm RGP. As can be seen from the figure, the RMSE of the algorithm of the present scheme will increase when the target maneuvers, and the tracking accuracy will be reduced to a certain extent, and the subsequent RMSE will quickly decrease. This is because when the target maneuver is detected, the algorithm of the present scheme needs to re-initialize the learning of the motion pattern of the target according to less measurement data, and with the passage of time, the measurement data is constantly increasing, and the model parameters can be constantly corrected. Table 1 is the average position RMSE of different algorithms in six simulation scenarios. As can be seen from the position RMSE, the RMSE value of the algorithm of the present scheme is smaller than that of the comparative algorithm RGP in the six simulation scenarios, indicating that the filtering effect of the algorithm RGP-BOCPD of the present scheme is better than that of the comparative algorithm. The RMSE value of the comparative algorithm is greater than the measurement RMSE value in scenarios S2, S4, S5 and S6. This is because in the scenario with strong maneuverability, the algorithm system parameter updates slowly, and cannot timely describe the motion pattern of the target after the maneuver, resulting in the decline of the target tracking performance and even loss of tracking. The RMSE value of the algorithm of the present scheme is smaller than the measurement RMSE value in all scenarios, indicating that the algorithm of the present scheme can effectively filter in all maneuvering target tracking scenarios. This is because we combine the Bayesian online change point detection, when the target maneuver is detected, the motion pattern of the target is relearned by using the Gaussian process, which can timely and effectively represent the target motion, thereby ensuring the effectiveness and robustness of the filtering of the algorithm of the present scheme.

Claims

1. A maneuvering target tracking method based on non-stationary Gaussian process, characterized in that: The following steps are involved: In step 1, the subject uses a Gaussian process to construct an online learning target motion model to describe the possible motion of the maneuvering target, a measurement model to represent the measurement data obtained by the sensor, and a change point model to represent the time when the target motion 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 target's measurement data through sensors and uses the measurement data from the previous multiple moments as the initial training data for each run 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 run length of the parameter to be estimated at time k based on 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, wherein: The step 1 specifically includes: Step 1.1, introduce a non-negative discrete integer variable r k ∈[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; when the target motion mode changes, let r k =1, at this time the measurement 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 measurement 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 moment in the run with 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 when ; The length of the run at time k is r k The Gaussian process when , its mean function is 0, and the covariance kernel function is The covariance kernel function is usually the squared exponential kernel (SE) function; In step 1.3, in each running length, the original measurement data obtained by the sensor is subjected to coordinate transformation and registration preprocessing to obtain the target position measurement. 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 the induced point set at time k-1, represents the input of the d-th induction point at time k-1; Initial induction point set is the measurement data at d moments before the start of the run length; since the target will maneuver during the movement, the hyperparameter vector is unknown, and hyperparameters need to be learned 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 is Hyperparameter vector and the running length for k moments is r k The predicted state 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 moment, 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 hazard 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 target position 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, wherein: The step 2 specifically includes: When time k>d, the subject performs all runs of length r at each time according to the target motion model defined in step 1. k Calculate 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 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 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 the sigma point and the corresponding covariance matrix as follows: in, Represents 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 zero vector; the superscript -1 of the parameter represents the matrix inversion, and the superscript T represents the matrix transpose; The augmented state vector ψ can be calculated through weighted fusion 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 a run length of r k The predicted state 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 of 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 variables o k The cross covariance between for The transpose of Step 2.3, using the measurement model in step 1.3 and the decomposition results 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 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, Represents a Gaussian distribution.

4. The method for tracking a maneuvering target based on a non-stationary Gaussian process according to claim 3, wherein: The step 3 specifically includes: Step 3.1, joint distribution p(z 1：k , r k ) can be recursively expressed 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 measurement data z at the previous k-1 moments 1：(k-1) , run length r k-1 The joint distribution of Step 3.2, change 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 the joint distribution p(z 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, wherein: 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 : 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 r k The estimated mean of the observable state at Unobservable state estimated 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 it to the end of the induction point set and delete the first induction point in the induction point set, thereby completing 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 k|k ) and estimated covariance ∑(x k|k ): The state estimation mean μ(x k|k ) contains the required position estimate 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 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 method for tracking a maneuvering target based on a non-stationary Gaussian process according to any one of claims 1 to 6 is implemented.

8. A computer-readable storage medium storing a computer program; wherein: 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