CN104237853B - A kind of for the particle filter method of trace point mark sequence before multi frame detection - Google Patents

Abstract

The present invention provides a kind of for the particle filter method of trace point mark sequence before multi frame detection.Before the domestic invention of this law rice utilizes multi frame detection, tracking processes the short flight path that obtains and is filtered dbjective state processing, and can improve the precision to Target state estimator；Tracking before multi frame detection is processed the short flight path obtained be further processed, it is achieved that the long-time tracking to target, the problem that before solving multi frame detection, tracking is not provided that target complete flight path information；Can dbjective state is estimated in real time with iteration, and do not introduce too much amount of storage and amount of calculation；Predicted by particle state, the missing inspection problem being likely to occur can be followed the tracks of before effectively solving low signal-to-noise ratio dim target detection, it is ensured that the seriality of flight path.

Description

A kind of for the particle filter method of trace point mark sequence before multi frame detection

Technical field

The invention belongs to radar data processing technology field, it is particularly to radar weak target tracking technique field.

Background technology

The complexity day by day of radar detection environment is increasingly mature with target stealth technology so that radar timely and reliably monitors and floods Weak target in strong clutter becomes more difficult.Improve radar and the detecting and tracking performance of weak target is had important theory It is worth and practical significance.

Before detection, tracking technique is a kind of a kind of technology that weak target carries out in low signal-to-noise ratio environment detecting and tracking.With tradition Detection method different, follow the tracks of not single frames before detection and announce measurement, but after many frame data are processed, announce mesh simultaneously Target measurement and flight path.Before realizing detection, the method for tracking technique includes particle filter, dynamic programming, Hough transformation etc.. Before detection based on particle filter, the track initiation ability of tracking is not as the method such as dynamic programming and Hough transformation, and the party Method needs calculating probability, computationally intensive.But, before the multi frame detection such as dynamic programming, tracking processes and once can only obtain target Short flight path, it is impossible to obtain the complete flight path of target.And, before some multi frame detections, tracking obtains at discrete space State estimation result, there is certain loss in estimated accuracy.

Summary of the invention

The technical problem to be solved is, for improving the estimated accuracy to dbjective state, it is provided that one is used for processing multiframe The filtering method of the short flight path obtained is followed the tracks of before detection.

The present invention solves that above-mentioned technical problem be employed technical scheme comprise that, a kind of for trace point mark sequence before multi frame detection Particle filter method, comprises the following steps:

Step 1) initialization of variable:

Before initializing total simulation time K, one-time detection, tracking processes frame number N；Before detection, tracking processes thresholding V_T, detection before with Before track measurement possible range γ, a multi frame detection, tracking processes the detection probability of n-th frameRadar scanning Cycle T, dbjective state transfer matrix F, observing matrix H, process noise covariance Q, population N_s, x and y direction Distance range [x_min,x_max] and [y_min,y_max], the velocity interval in x and y directionWith

Step 2) particle collectionInitialize:

Initialization time variable k=1 and i-th particle stateIn position, x and y directionAnd speed state Represent transposition；

After initialization particle collection has initialized, enter step 4)；

Step 3) status predication:

The particle set pair k moment utilizing the k-1 moment is predicted:Obtain predicting grain Subset Being the process noise of process noise covariance Q, F is dbjective state transfer matrix；

When in i-th particle state, maximum magnitude value more than its maximum magnitude, is then assigned to this element, and arranges this grain by element The weights that son is corresponding are set to 0；

Step 4) before multi frame detection tracking process:

Kth-N+1, k-N+2 is read from radar receiver ..., k frame echo data, as k≤N, then from radar receiver Read the 1,2nd ..., k frame echo data；In tracking process before kth time detection, there is measurement, then enter step 5), no Then enter step 7)；There is measurement and i.e. represent that measurement exceedes thresholding；

Step 5) right value update:

K≤N, then measurement isThe most more new particle weights:

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=1,2,...,N_s;$

Work as k > N, measurement isStart more new particle weights:

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_{k-N+1}})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_{k-N+2}})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=1,2,...,N_s;$

Represent that the last frame moment is t_kDetection before tracking process the l frame measurement that obtains,RepresentCorresponding The state of father's particle, as l≤k-1, usesRepresentThe state of corresponding father's particle；P represents likelihood probability；

Representing and round downwards, H is observing matrix, measuring range variable △_xAnd △_ySpan beAnd △_xAnd △_yAll it is not equal to zero；The meter of likelihood probability when γ=1 Without middle entry in formula；

Calculate update after particle weights andAgain to particle weights normalizedEnter Enter step 6)

Step 6) to the particle collection after right value updateCarry out system resampling；

6.1) initialize

6.2) for all i=2,3 ..., N_s, calculate

6.3) random number u is produced₁, u₁ObeyIn the range of be uniformly distributed；

6.4) i=1 is made, for all j=1,2 ..., N_s, perform calculating:

Work as u_j>c_i, then update i=i+1, return step 6.4) until u_j≤c_i；Make againAnd remember Record the father particle i of this particle^j=i；

Particle set representations after resampling isEnter step 7)

Step 7) Target state estimator:

To the particle state after all resamplingsWeighted averageObtain the Target state estimator result in kth moment；

Step 8) update k=k+1, as k≤K, return step 3), otherwise N frame particle filter terminates.

The present invention passes through above step, it is possible to processes, according to tracking before detection, the short flight path obtained and estimates the most in real time Dbjective state, obtains the complete flight path of target.

The invention has the beneficial effects as follows, before the present invention utilizes multi frame detection, tracking processes the short flight path that obtains and filters dbjective state Ripple processes, and can improve the precision to Target state estimator；Tracking before multi frame detection is processed the short flight path obtained do further Process, it is achieved that the long-time tracking to target, before solving multi frame detection, tracking is not provided that target complete flight path information Problem；Can dbjective state is estimated in real time with iteration, and do not introduce too much amount of storage and amount of calculation；Pass through Particle state is predicted, can follow the tracks of the missing inspection problem being likely to occur before effectively solving low signal-to-noise ratio dim target detection, it is ensured that boat The seriality of mark.

Accompanying drawing explanation

Fig. 1 is the flow chart of the present invention.

Fig. 2 is that in the embodiment of the present invention, the true flight path of target follows the tracks of, with the application present invention, the targetpath obtained.

Fig. 3 is to apply the inventive method to estimate the target location root-mean-square error of dbjective state in the embodiment of the present invention.

Fig. 4 is to apply the inventive method to estimate the target velocity root-mean-square error of dbjective state in the embodiment of the present invention.

Detailed description of the invention

As it is shown in figure 1, flow process of the present invention is as follows:

Step 1) initialization of variable:

Initialize total simulation time K；Before one-time detection, tracking processes frame number N；Target signal to noise ratio SNR；Before detection, tracking processes Thresholding V_T；Follow the tracks of before detection and measure possible range γ, γ ∈ { δ²；δ=1,2,3,4 ... }, δ represents detector unit number, before detection with Track measures the scope that not can exceed that between measurement and the real goal of tracking before possible range is expressed as detecting；Multiframe inspection Before surveying, tracking processes the detection probability of n-th frameDetection probability is known when following the tracks of before detecting；Thunder Reach scan period T；Dbjective state transfer matrix F；Observing matrix H；Process noise covariance Q；Population N_s；X and y side To distance range [x_min,x_max] and [y_min,y_max]；The velocity interval in x and y directionWith

Step 2) particle collectionInitialize:

Initialization time variable k=1 and i-th particle stateIn position, x and y directionAnd speed state

$x_k^i=x_{min}+(x_{max}-x_{min})\×rand\left(1\right),i=1,2,...,N_s$

${\stackrel{\·}{x}}_k^i={\stackrel{\·}{x}}_{min}+({\stackrel{\·}{x}}_{max}-{\stackrel{\·}{x}}_{min})\×rand\left(1\right),i=1,2,...,N_s$

$y_k^i=y_{min}+(y_{max}-y_{min})\×rand\left(1\right),i=1,2,...,N_s$

${\stackrel{\·}{y}}_k^i={\stackrel{\·}{y}}_{min}+({\stackrel{\·}{y}}_{max}-{\stackrel{\·}{y}}_{min})\×rand\left(1\right),i=1,2,...,N_s$

Rand (1) represents random according to be uniformly distributed on one [0,1] interval of generation Number；

I-th particle weightsIt is initialized as:

After particle collection has initialized, enter step 4)；

Step 3) status predication:

The particle set pair k moment utilizing the k-1 moment is predicted:Obtain predicting grain Subset It it is the process noise of process noise covariance Q；

When in i-th particle state, maximum magnitude value more than its maximum magnitude, is then assigned to this element, and arranges this grain by element The weights that son is corresponding are set to 0；

Step 4) before multi frame detection tracking process:

Kth-N+1, k-N+2 is read from radar receiver ..., k frame echo data, as k≤N, then from radar receiver Read the 1,2nd ..., k frame echo data；In tracking process before kth time detection, there is measurement, then enter step 5), no Then enter step 7)；

Step 5) right value update:

K≤N, then measurement isThe most more new particle weights:

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=1,2,...,N_s$

Work as k > N, measurement isStart more new particle weights:

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_{k-N+1}})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_{k-N+2}})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=1,2,...,N_s$

Represent that the last frame moment is t_kDetection before tracking process the l frame measurement that obtains,RepresentCorresponding The state of father's particle, as l≤k-1, usesRepresentThe state of corresponding father's particle；P represents

Likelihood probability；

Representing and round downwards, H is observing matrix, measuring range variable △_xAnd △_ySpan beAnd △_xAnd △_yAll it is not equal to zero；

Calculate update after particle weights andAgain to particle weights normalizedEnter Enter step 6)

Step 6) to the particle collection after right value updateCarry out system resampling；

6.1) initialize

6.2) for all i=2,3 ..., N_s, calculate

6.3) random number u is produced₁, u₁ObeyIn the range of be uniformly distributed；

6.4) i=1 is made, for all j=1,2 ..., N_s, perform calculating:

Work as u_j>c_i, then update i=i+1, return step 6.4) until u_j≤c_i；Make againAnd remember Record the father particle i of this particle^j=i；

Particle set representations after resampling isEnter step 7)

Step 7) Target state estimator:

To the particle state after all resamplingsWeighted averageObtain the Target state estimator result in kth moment；

Step 8, renewal k=k+1, as k≤K, return step 3), otherwise N frame particle filter terminates.

Embodiment

The present invention mainly uses the method for Computer Simulation to verify, institute is in steps, conclusion is all on MATLAB-R2010b Checking is correct.It is embodied as step as follows:

Step 1, initialization of variable:

Initializing variable includes: total simulation time K=30；Before one-time detection, tracking processes frame number N=6；Target signal to noise ratio SNR=8dB；Before detection, tracking processes thresholding V_T=13.2847；Follow the tracks of before detection and measure possible range γ=7²；Multiframe inspection Before surveying, tracking processes the 1,2nd ..., the detection probability of N frame is respectively 0.876,0.892,0.918,0.910,0.908,0.838；Radar scanning Cycle T=1s；Dbjective state transfer matrix F

$F=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]$

Observing matrix H

$H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

Process noise covariance Q

$Q=\left[\begin{array}{cccc}0.001\×T^3/3& 0.001\×T^2/2& 0& 0\\ 0.001\×T^2/2& 0.001\×T& 0& 0\\ 0& 0& 0.001\×T^3/3& 0.001\×T^2/2\\ 0& 0& 0.001\×T^2/2& 0.001\×T\end{array}\right]$

Population N_s=8000；The distance range in x and y direction is [1,59]；The velocity interval in x and y direction is [0.6,1.2].

Step 2, particle collection initialize:

Initialization time variable k=1.Particle x and position, y direction and speed state initialize

$x_k^i=1+58\×rand\left(1\right),i=1,2,...,N_s$

${\stackrel{\·}{x}}_k^i=0.6+0.6\×rand\left(1\right),i=1,2,...,N_s$

$y_k^i=1+58\×rand\left(1\right),i=1,2,...,N_s$

${\stackrel{\·}{y}}_k^i=0.6+0.6\×rand\left(1\right),i=1,2,...,N_s$

Wherein rand (1) represents according to being uniformly distributed the random number produced on [0,1] interval.Order

$x_k^i={\[x_k^i,{\stackrel{\·}{x}}_k^i,y_k^i,{\stackrel{\·}{y}}_k^i\]}^T,i=1,2,...,N_s$

Particle weight initialization

$w_k^i=\frac{1}{N_s},i=1,2,...,N_s$

Then, primary collection is obtained

Step 3, status predication:

If k=1, then skip this step.Otherwise the particle collection in k-1 moment is carried out one-step prediction

$x_{k|k-1}^i={\mathrm{Fx}}_{k-1}^i+v_k^i,i=1,2,...,N_s$

Obtain predicting particle collectionWhereinBe average be null covariance matrix be the process noise of Q.

If i-th particle stateOrThen makeOrIfOrThen makeOrAnd weights corresponding for this particle are set to 0。

Before step 4, multi frame detection, tracking processes:

Kth-N+1, k-N+2 is read from radar receiver ..., and k frame echo data (if k≤N, then from radar receiver Middle reading the 1st, 2 ..., N frame echo data), before carrying out multi frame detection, tracking processes.

Step 5, right value update:

OrderRepresent that the last frame moment is t_kDetection before tracking process the l frame measurement that obtains.If this detection Measurement is there is in front tracking in processingIf (k≤N, then measurement is), then update and predict particle weights:

Convenient in order to represent, useRepresentThe state of corresponding father's particle, as l≤k-1, usesRepresentCorresponding father's grain The state of son.As k > N time, the most more new particle weights

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_{k-N+1}})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_{k-N+2}})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=1,2,...,N_s$

Wherein, likelihood probability

△_xAnd △_ySpan beAnd △_xAnd △_yAll it is not equal to zero.

As k≤N, the most more new particle weights

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=1,2,...,N_s$

Calculate particle weights and

$\mathrm{\Ξ}=\underset{i=1}{\overset{N_s}{\Σ}}{w}_k^i$

To particle weights normalization in the case of Ξ is not zero

$x_k^i=\frac{w_k^i}{\mathrm{\Ξ}},i=1,2,...,N_s$

If tracking does not exist measurement in processing before this detection, i.e. do not have state to cross thresholding, then skip this step.

Step 6, resampling:

If tracking does not exist measurement in processing before this detection, i.e. do not have state to cross thresholding, then skip this step.No Then to the particle collection after right value updateCarry out system resampling.

6.1, initialize

6.2, for all i=2,3 ..., N_s, calculate

$c_i=c_{i-1}+w_k^i$

6.3, u is randomly generated₁, it is obeyedIn the range of be uniformly distributed, i.e.

6.4, i=1 is made, for all j=1,2 ..., N_s, perform calculated below:

$u_j=u_1+(j-1)N_s^{-1}$

If u_j>c_i, then circulation performs i:=i+1, until u_j≤c_i.Then makeAnd record this particle Father particle i^j=i.

Particle set representations after resampling is

Step 7, Target state estimator:

To all particle statesWeighted averageObtain the Target state estimator result in kth moment.

Step 8, make k:=k+1, if k≤K, repeat step 3～step 7.

Accompanying drawing 2 gives the restoration result in this simulation example to targetpath, successfully provides the complete flight path of target Information, and there is not missing inspection problem.Accompanying drawing 3 gives the target location root-mean-square error utilizing the present invention to obtain in the present embodiment, Accompanying drawing 4 gives the target velocity root-mean-square error utilizing the present invention to obtain in the present embodiment, is 200 Monte Carlo experiments Statistical result, result shows that the method has higher and stable estimated accuracy.

Claims (2)

1. one kind for the particle filter method of trace point mark sequence before multi frame detection, it is characterised in that comprise the following steps:

Step 1) initialization of variable:

Before initializing total simulation time K, one-time detection, tracking processes frame number N；Before detection, tracking processes thresholding V_T, detection before with Before track measurement possible range γ, a multi frame detection, tracking processes the detection probability of n-th frameN=1,2 ..., N, radar scanning Cycle T, dbjective state transfer matrix F, observing matrix H, process noise covariance Q, population N_s, x and y direction Distance range [x_min,x_max] and [y_min,y_max], the velocity interval in x and y directionWith

Step 2) particle collectionInitialize:

Initialization time variable k=1 and i-th particle stateIn position, x and y directionAnd speed state [·]^TRepresent transposition；

After initialization particle collection has initialized, enter step 4)；

Step 3) status predication:

The particle set pair k moment utilizing the k-1 moment is predicted:Obtain predicting grain Subset Being the process noise of process noise covariance Q, F is dbjective state transfer matrix；

When in i-th particle state, maximum magnitude value more than its maximum magnitude, is then assigned to this element, and arranges this grain by element The weights that son is corresponding are set to 0；

Step 4) before multi frame detection tracking process:

Work as k > N, from radar receiver, read kth-N+1, k-N+2 ..., k frame echo data, as k≤N, then from thunder Reach and receiver read the 1,2nd ..., N frame echo data；In tracking processes before kth time detection, there is measurement, then enter Step 5), otherwise enter step 7)；There is measurement and i.e. represent that measurement exceedes tracking process thresholding V before detection_T；

Step 5) right value update:

K≤N, then measurement isThe most more new particle weights:

$w_k^i=P\left(z_1^N\right|x_1^{i_1})P\left(z_2^N\right|x_2^{i_2})...P\left(z_k^N\right|x_{k|k-1}^i),i=1,2,...,N_s;$

Work as k > N, measurement isStart more new particle weights:

$w_k^i=P\left(z_{k-N+1}^{t_k}\right|x_{k-N+1}^{i_{k-N+1}})P\left(z_{k-N+2}^{t_k}\right|x_{k-N+2}^{i_{k-N+2}})...P\left(z_k^{t_k}\right|x_{k|k-1}^i),i=1,2,...,N_s;$

Represent that the last frame moment is t_kDetection before tracking process the l frame measurement that obtains,RepresentCorresponding The state of father's particle, as l≤k-1, usesRepresentThe state of corresponding father's particle；P represents likelihood probability；

Representing and round downwards, H is observing matrix, measuring range variable Δ_xAnd Δ_ySpan beAnd Δ_xAnd Δ_yAll it is not equal to zero；The meter of likelihood probability when γ=1 Without middle entry in formula；

Calculate update after particle weights andAgain to particle weights normalizedEnter Enter step 6)

Step 6) to the particle collection after right value updateCarry out system resampling；

6.1) initialize

6.2) for all i=2,3 ..., N_s, calculate

6.3) random number u is produced₁, u₁ObeyIn the range of be uniformly distributed；

6.4) i=1 is made, for all j=1,2 ..., N_s, perform calculating:

Work as u_j>c_i, then update i=i+1, return step 6.4) until u_j≤c_i；Make againAnd remember Record the father particle i of this particle^j=i；

Particle set representations after resampling isEnter step 7)

Step 7) Target state estimator:

To the particle state after all resamplingsWeighted averageObtain the Target state estimator result in kth moment；

Step 8) update k=k+1, as k≤K, return step 3), otherwise N frame particle filter terminates.

It is a kind of for the particle filter method of trace point mark sequence before multi frame detection, it is characterised in that Initialization time variable k=1 and i-th particle stateIn position, x and y directionAnd speed state Method particularly includes:

$x_k^i=x_{min}+(x_{max}-x_{min})\×rand\left(1\right),i=1,2,...,N_s$

${\stackrel{\·}{x}}_k^i={\stackrel{\·}{x}}_{min}+({\stackrel{\·}{x}}_{max}-{\stackrel{\·}{x}}_{min})\×rand\left(1\right),i=1,2,...,N_s$

$y_k^i=y_{\min }+(y_{max}-y_{min})\×rand\left(1\right),i=1,2,...,N_s$

${\stackrel{\·}{y}}_k^i={\stackrel{\·}{y}}_{min}+({\stackrel{\·}{y}}_{max}-{\stackrel{\·}{y}}_{min})\×rand\left(1\right),i=1,2,...,N_s$

Rand (1) represents according to being uniformly distributed the random number produced on [0,1] interval；

Initialize i-th particle weightsConcrete grammar for for: