CN104750998B

CN104750998B - Target tracking method of passive multi-sensor based on density filter

Info

Publication number: CN104750998B
Application number: CN201510166369.4A
Authority: CN
Inventors: 李翠芸; 代标; 姬红兵; 李宁
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2015-04-09
Filing date: 2015-04-09
Publication date: 2017-08-25
Anticipated expiration: 2035-04-09
Also published as: CN104750998A

Abstract

, can not be the problem of the noise intensity unknown accurate tracking multiple target with the newborn unknown situation progress of target to solve prior art the invention discloses a kind of passive multi-sensor multi-object tracking method based on density filter.Its implementation is：Corresponding coordinate under passive multi-sensor bearingonly measurement is obtained by the least square Cross Location Method based on Newton iteration first, and the coordinate at former and later two moment is associated, newborn object set is obtained；Secondly, newborn target sample collection is produced according to the distribution of the state of newborn target；Then multiple target tracking is carried out using density filter.The present invention can not only quickly obtain newborn target sample collection, reduce target sample number, improve target following efficiency, and noise intensity can be updated in real time, available for multiple target tracking while dbjective state is estimated.

Description

Target tracking method of passive multi-sensor based on density filter

Technical field

The invention belongs to technical field of guidance, it is related to unknown newborn multiple target tracking, specifically a kind of passive multi-sensor Method for tracking target, available for systems such as infrared guidances.

Background technology

The basic object of target following is the number according to the measurement set estimation current goal for mixing clutter and noise And state.Due to the appearance of target, disappearance and the presence of target derivatization process, the number of each moment target may all occur Change.Passive multi-sensor can not obtain the range information of target, and its target information obtained is usually the azimuth of target With pitching angle information, passive multi-sensor target following is substantially a Nonlinear Tracking problem, how by nonlinear filtering It is always the focus and difficult point of academia and engineer applied area research applied to passive multi-sensor target following.

The probability hypothesis density filtering and density filter realized based on particle method are solving Nonlinear Multiobjective tracking There is preferably application in problem.Although probability hypothesis density filtering can be good at being applied to passive multi-sensor target following In problem, but its newborn dbjective state is known a priori, and can not estimate to measure noise intensity in real time.Density filter is used Continuously-discrete space represents new life and the extinction of target by the mutual conversion of continuous space and discrete space to Target Modeling, Estimation noise intensity can be measured while target is tracked.Although density filter can obtain new by region uniform sampling Raw target information, but population can be sharply increased with the expansion of target area, the missing inspection situation when particle density is relatively low Than more serious, when particle density is big compared with high hour operation quantity, track less efficient.

The content of the invention

It is an object of the invention to for problem above, propose a kind of passive multi-sensor target based on density filter Tracking, target is accurately tracked by realize under the unknown newborn unknown situation with target of noise intensity.

Realizing the key of the technology of the present invention is：, will be by using many base least square cross bearingses based on Newton iteration The azimuth that dynamic multisensor is obtained is converted into coordinate；The velocity interval setpoint distance threshold value moved according to target, by continuous two Coordinate points of the individual moment distance in threshold range are considered as newborn target point.

According to above-mentioned thinking, implementation steps of the invention include as follows：

(1) initial time t=0 is made, is that sensor produces N according to the initial distribution of dbjective state_t|tIndividual initial sample structure Into target sample collection：Wherein,I-th of target sample of t is represented,Represent t i-th The intensity of target sample；

If initial time is without newborn target, i.e., newborn target sample collectionFor empty set, wherein,Table Show i-th of newborn target sample of t,Represent the intensity of i-th of newborn target sample of t, J_tRepresent that t is newborn Target sample number；

(2) the target sample collection of t is mergedWith newborn target sample collectionObtain New set：WhereinWithI-th of sample after merging and its intensity are represented respectively；

(3) whether be empty, if the sample set after merging is sky, perform step if judging the sample set after t merging (9) step (4), otherwise, is performed；

(4) the merging sample set obtained according to (2)It is pre- to the t+1 moment that t is calculated with state equation Survey sample setWherein,T is represented to i-th of forecast sample of t+1 moment,Represent The intensity of i-th of forecast sample, N_t+1|t=N_t|t+J_tRepresent forecast sample number；

(5) according to forecast sample collectionWith t without goal hypothesis intensity f_t|tWhen (φ) predicts t+1 Carve without goal hypothesis intensity f_t+1|t(φ)；According to forecast sample collectionWith the measurement collection at t+1 momentPredict the measurement intensity at t+1 momentWherein,Represent that the t+1 moment measures for j-th, N_z,t+1Represent t The measurement number at+1 moment；

(6) by measurement intensityWith measurement likelihoodCalculate the probability W that target is present_j, work as W_j＞ τ When, it is believed that target is present, and estimates dbjective stateWherein, τ is the given threshold that target is present；

(7) prediction is updated respectively without goal hypothesis intensity f_t+1|t(φ) and forecast sample collectionObtain t+ 1 moment without goal hypothesis intensity f_t+1|t+1(φ) and sample set

(8) following parameter is obtained according to the result of step (7)：

8a) to the sample set after renewalIn intensity summation, obtain the number of targets that the t+1 moment estimates Desired value η_t+1；

8b) using without goal hypothesis intensity f_t+1|t+1(φ) and forecast sample number N_t+1|t, calculate resampling sample number N_t+1|t+1；

8c) using without goal hypothesis intensity f_t+1|t+1(φ), without goal hypothesis transition probability ψ_t+1(φ | φ) and pre- test sample This number N_t+1|tCalculate newborn sample number J_t+1；

8d) from sample setMiddle resampling N_t+1|t+1Individual sample, the sample set after being sampled：

(9) newborn coordinates of targets and newborn target velocity are calculated：

Many base least square cross bearingses 9a) are utilized, all measurements at t and t+1 moment are converted into coordinate, point It is not expressed asWithFilter out all satisfactionsCoordinate, structure Into newborn coordinates of targets collection：Wherein,WithRepresent that t and t+1 moment measure for j-th respectively The coordinate of conversion,WithJ-th of newborn coordinates of targets of t and t+1 moment, ε are represented respectively_minAnd ε_maxDifference table Show apart from lower limit and apart from the upper limit, N_zFor all coordinate numbers for meeting condition；

The speed of newborn target 9b) is calculated according to the coordinate of newborn target：Wherein, Δ T Represent time interval；

(10) according to newborn coordinates of targets collectionWith speed collectionConstitute newborn object setWhether be empty, if it is empty, then in target region [a if judging newborn object set_x,b_x]*[a_y,b_y] in uniformly adopt Sample J_t+1Individual target sample；Otherwise, each newborn target is distributed as according to the state of newborn target to produceIndividual target sample, Constitute newborn target sample collectionWherein,Represent newborn target strength；

(11) t=t+1 is made, return to step (2) continues to track target.

The present invention has advantages below：

1) bearingonly measurement is converted into coordinate, utilizes seat by the present invention due to using many base least square cross bearingses Information is marked, all coordinate points met in threshold range are judged as newborn target point, and newborn target velocity is estimated, because And can be in the case where newborn target is unknown, the newborn target of optional position in tracing area.

2) present invention is due to using the estimation to newborn coordinates of targets, reducing sample used in object tracking process Number, reduces operand, improves operational efficiency.

Brief description of the drawings

Fig. 1 be the present invention realize general flow chart；

Fig. 2 is the least square cross bearing sub-process figure based on Newton iteration in the present invention；

Fig. 3 is the scene graph that the present invention is emulated；

Fig. 4 (a) is the trajectory diagram with present invention tracking target；

Fig. 4 (b) is the target number figure estimated with the present invention；

Fig. 5 (a) is the dbjective state OSPA Error Graphs estimated with the present invention；

Fig. 5 (b) is the target velocity OSPA Error Graphs estimated with the present invention；

Fig. 6 (a) is with target sample number figure used in present invention tracking object procedure；

Fig. 6 (b) is the run time figure of the present invention.

Embodiment

First, basic theory introduction

1. system equation

Under cartesian coordinate system, system mode takes x, the position in y directions and speed, can set up following Nonlinear Dynamic System model：

x_t+1=Fx_t+Ge_t 1)

z_t=h (x_t)+w_t 2)

Wherein, x_t=[x_t,v_x,t,y_t,v_y,t]^T, T represents transposition, x_t,y_tRepresent target on x directions and y directions respectively Coordinate, v_x,t,v_y,tSpeed of the target on x directions and y directions is represented respectively, and subscript t ∈ [1, N] represent time, state-noise e_tObeying variance isZero-mean gaussian distribution, diag [@] represent diagonalization element therein, σ_xAnd σ_y Represent that location criteria is poor, F, G is respectively state-transition matrix and input matrix, and h is nonlinear function, measure noise w_tObedience side Difference is distributed for R zero-mean gaussian, e_tWith w_tIt is separate, z_tFor the measuring value of sensor.

In the present invention assume passive sensor can only observed object azimuth information, therefore nonlinear function h definition such as Under：

Wherein, arctan () represents arc tangent, x_s,y_sFor the position of sensor.

2. more than base least square cross bearing

Assuming that many base station s₁,s₂,…,s_nCan obtain the orientation angle measurements of target at each moment, the measurement of t to Measure Z=[θ_t1,θ_t2,…,θ_tn], the actual position X=[x of target_t,y_t], the position of many base stations is (x_s1,y_s1),(x_s2, y_s2),…,(x_sn,y_sn)。

By formula 2) measurement model that can obtain alignment system is：

Z=h (X)+w_θ 4)

Wherein nonlinear function is：

w_θFor observation noise vector, w_θ=[w_θ1,w_θ2,…,w_θn]^T, w_θiFor the observation noise of i-th of base station；

According to the definition w of least square_θ ^TPw_θ=min, then state vector X estimateMeet：

Wherein, min represents minimum, and P is covariance matrix, due to Z^TZ be a constant, therefore formula 6) be equivalent to target letter Number：

Formula 7) be many base least square cross bearing systems least-squares estimation model.

Formula 7) a certain approximate solution can be approached by Newton iterative method recursion so that

Wherein, → represent infinite approach, the approximate solution obtainedCan serve as formula 7) numerical solution.WillIn kth Secondary iterative state X^kNearby utilize second order Taylor series expansion：

Wherein：

Wherein,Represent in kth time iterative state X^kThe value at place, formula 9) inMake formula 9) to Δ X^k's Single order local derviation is zero, is obtained：

g^k+(ΔX^k)^TG^k=0 12)

Work as G^kWhen nonsingular

(ΔX^k)^T=-(G^k)^-1(g^k)^T 13)

Wherein, ()^-1Represent inversion operation.

Then iterative formula can be written as：

X^k+1=X^k+ΔX^k=X^k-(G^k)^-1(g^k)^T 14)

Formula 14) be exactly Newton iteration fundamental formular, stopping criterion for iteration is：

R(X^k+1)≡R(X^k) 15)

Wherein, ≡ represents to be constantly equal to.

Two：The tracking of the present invention

Reference picture 1, tracking step of the invention is as follows：

Step 1. initialized target sample

Initial time t=0 is made, is that sensor produces N according to the initial distribution of dbjective state_t|tIndividual initial sample constitutes mesh Mark sample set：Wherein,I-th of target sample of t is represented,Represent i-th of target of t The intensity of sample；

Step 2. merges target sample

Merge the target sample collection of tWith newborn target sample collectionObtain new Set：WhereinWithI-th of sample after merging and its intensity are represented respectively；

Step 3. judges whether the sample set after t merging is empty, if the sample set after merging is sky, performs step 8, otherwise, perform step 4；

Step 4. predicts target sample collection, without goal hypothesis intensity and measurement intensity.

(4.1) the merging sample set obtained according to tT is obtained to the prediction mesh at t+1 moment Mark sample setWherein,T is represented to i-th of forecast sample of t+1 moment, F is represented State-transition matrix, G represents input matrix, e_tRepresent state-noise,The intensity of i-th of forecast sample is represented, N_t+1|t=N_t|t+J_tRepresent forecast sample number；

(4.2) predict t to the t+1 moment without goal hypothesis intensity f_t+1|t(φ)：

Wherein,T is represented to the non-existent intensity of t+1 moment targets,Represent t to the t+1 moment The intensity of target termination, its calculation formula is as follows：

Wherein, ψ_t+1(φ | φ) no goal hypothesis transition probability is indicated,Represent the probability that target disappears, f_t|t (φ) represents the intensity without goal hypothesis of t；

(4.3) according to the measurement collection at t+1 momentPrediction measures intensity

Wherein, j ∈ [1, N_z,t+1], N_z,t+1The measurement number at t+1 moment is represented,Represent measurement source in target-like The intensity of state,Measurement source is represented in the intensity without goal hypothesis, its calculation formula is as follows：

Wherein,Represent to measure the likelihood for coming from no goal hypothesis,Indicate the inspection of no goal hypothesis Survey probability,Represent target detection probability,Represent to measure likelihood,With Calculation formula it is as follows：

Wherein,！Factorial computing is represented, Ω represents covariance matrix.

Step 5. estimates dbjective state.

(5.1) j-th is measuredCalculate and measureFrom the probability of target：

(5.2) W is worked as_jDuring ＞ τ, it is believed that j-th of measurement comes from target, and estimates dbjective state： Wherein,τ is the given threshold that target is present.

Step 6. updates target sample intensity and without goal hypothesis intensity.

(6.1) using following formula to target strengthIt is updated, the target strength after being updated

(6.2) using following formula to after prediction without goal hypothesis intensity f_t+1|t(φ) is updated, the nothing after being updated Goal hypothesis intensity f_t+1|t+1(φ)：

Step 7. is according to the target sample collection after renewalWith the target strength after renewalObtain mesh Mark sample set

(7.1) to formula 25) summation of obtained sample intensity, the number of targets desired value that the t+1 moment estimates is obtained, is expressed as η_t+1：

(7.2) according to above-mentioned number of targets desired value η_t+1With after renewal without goal hypothesis intensity f_t+1|t+1(φ), calculates weight Oversampling ratio p_s：

(7.3) forecast sample number is multiplied with resampling ratio and obtains resampling sample number N_t+1|t+1：

N_t+1|t+1=N_t+1|t·p_s 29)

(7.4) from the target sample collection after renewalMiddle resampling N_t+1|t+1Individual sample, is obtained after sampling Target sample collection：

Step 8. measures conversion

Reference picture 2, this step is implemented as follows：

(8.1) set initial state value as：According to j-th of measurementInitial target functional value R is calculated as follows (X^k):

Wherein, j ∈ [1, N_z,t+1],x_t+1And y_t+1X directions and y directions coordinate are represented respectively, and k is represented Iterations, nonlinear functionIt is calculated as follows：

Wherein, (x_si,y_si) represent i-th of observation station position, i=1 ..., n；

(8.2) Jacobi matrixes g is calculated as follows^k：

(8.3) Hessian matrixes G is calculated as follows^k：

(8.4) an iteration state value is calculated as follows

Wherein, (G^k)^-1Represent Hessian matrixes G^kIt is inverse, T represents to seek the transposition of matrix；

(8.5) an iteration target function value is calculated as follows

(8.6) judgeWithIt is whether equal, ifThen export j-th of measurement The coordinate of conversion：Otherwise, makeK=k+1, return to step (8.2) is carried out next Secondary iteration；Wherein,

(8.7) step (8.1)~(8.6) are pressed, all measurements of t and t+1 moment are converted into coordinate, are expressed as：With

Step 9. estimates newborn target

(9.1) newborn target sample number J is calculated according to following formula_t+1：

Judge J_t+1Whether it is zero or infinitely great, if J_t+1It is zero or infinitely great, then makes J_t+1=N_newNumPar, its Middle N_newNewborn target number is represented, numPar represents the sampled targets sample number of each newborn target；

(9.2) result obtained according to step 8, the coordinate set of conversion is measured from tMeasure and turn with the t+1 moment The coordinate set of changeIn filter out and all meet ε_min＜ d ＜ ε_maxCoordinate constitute newborn coordinates of targets collectionWherein,Denotation coordinationAnd coordinateDistance, N_zFor The number of all coordinates of targets for meeting condition, ε_minAnd ε_maxThe lower and upper limit apart from d are represented respectively；

(9.3) speed of j-th of newborn target of t+1 moment is calculated：Wherein, j ∈ [1, N_z], Δ T represents time interval；

(9.4) according to newborn coordinates of targets collectionWith speed collectionConstitute newborn object setWhether be empty, if it is empty, then in target region [a if judging newborn object set_x,b_x]*[a_y,b_y] in uniformly adopt Sample J_t+1Individual target sample；Otherwise, each newborn target is distributed as according to the state of newborn target to produceIndividual target sample, Constitute newborn target sample collectionWherein,Represent newborn target strength.

Step 10. makes t=t+1, and return to step 2 continues to track target.

The design sketch of the present invention can be further illustrated by following experiment simulation：

1. simulated conditions and parameter

Simulating scenes such as Fig. 3 represents that the time of day of the target of each in simulating scenes is x=[x, v_x,y,v_y]^T, its In, x, y represents coordinate of the target on cartesian coordinate system x directions and y directions, v respectively_x,v_yRepresent target in x directions respectively With the speed on y directions, state equation and the measurement equation such as formula 1 of target) and 2), wherein：

Wherein, T_sSampling time interval is represented, simulation parameter is as shown in table 1：

The experiment simulation parameter of table 1

2. emulation content and interpretation of result

Pure orientation angle tracking to six targets under three sensor conditions carries out emulation experiment, and target is in two dimensional surface It is interior motion real trace as shown in figure 3, emulation in the present invention tracking LSCL-iFilter with using uniform sampling with The conventional strength wave filter U-iFilter method for tracking target of the unknown newborn target of track is contrasted.

Emulation 1：The track of target and target number are emulated with this method, simulation result as shown in figure 4, wherein, Fig. 4 (a) is target following trajectory diagram, and Fig. 4 (b) is target number analogous diagram, as can be seen that passing through from Fig. 4 (a) and Fig. 4 (b) Often there is target detection leakage phenomenon when tracking unknown newborn target in uniform sampling, fails effective tracking target, the present invention in real time Can postpone two moment estimates newborn target, and number of targets estimation is stable, seldom occurs more estimating and estimating phenomenon less.

Emulation 2：Dbjective state error and velocity error are emulated with this method, simulation result as shown in figure 5, its In, Fig. 5 (a) is dbjective state OSPA Error Graphs, and Fig. 5 (b) is target velocity OSPA Error Graphs, can from Fig. 5 (a) and Fig. 5 (b) To find out, error of the invention is consistently lower than the error that unknown newborn goal approach is tracked by uniform sampling.Fig. 5 (a) and Fig. 5 (b) the two larger peak values occurred in are both present in two moment that newborn target is postponed estimation.After newborn target state estimator goes out, Its state error and velocity error can tend to be steady.

Emulation 3：Sample number and run time are emulated with the present invention, simulation result is as shown in fig. 6, wherein, Fig. 6 (a) it is target number figure, Fig. 6 (b) is simulation time figure, as can be seen that the present invention is in tracking mesh from Fig. 6 (a) and Fig. 6 (b) The sample number and run time that timestamp is used are few, and it tracks efficiency far and tracks unknown newborn target higher than using uniform sampling Conventional strength wave filter.

Claims

1. a kind of target tracking method of passive multi-sensor based on density filter, comprises the following steps:

(1) initial time t=0 is made, is that sensor produces N according to the initial distribution of dbjective state_t|tIndividual initial sample constitutes mesh Mark sample set：Wherein,I-th of target sample of t is represented,Represent i-th of target sample of t This intensity；

If initial time is without newborn target, i.e., newborn target sample collectionFor empty set, wherein,When representing t I-th of newborn target sample is carved,Represent the intensity of i-th of newborn target sample of t, J_tRepresent the newborn target sample of t This number；

(2) the target sample collection of t is mergedWith newborn target sample collectionObtain new collection Close：WhereinWithI-th of sample after merging and its intensity are represented respectively；

(3) whether be empty, if the sample set after merging is sky, execution step (9) is no if judging the sample set after t merging Then, step (4) is performed；

(4) the merging sample set obtained according to (2)T is calculated to the pre- test sample at t+1 moment with state equation This collectionWherein,T is represented to i-th of forecast sample of t+1 moment,Represent i-th The intensity of individual forecast sample, N_t+1|t=N_t|t+J_tRepresent forecast sample number；

(5) according to forecast sample collectionWith t without goal hypothesis intensity f_t|t(φ) prediction t+1 moment Without goal hypothesis intensity f_t+1|t(φ)；According to forecast sample collectionWith the measurement collection at t+1 moment Predict the measurement intensity at t+1 momentWherein,Represent that the t+1 moment measures for j-th, N_z,t+1Represent the amount at t+1 moment Survey number；

(6) by measurement intensityWith measurement likelihoodCalculate the probability W that target is present_j, work as W_jDuring ＞ τ, recognize Exist for target, and estimate dbjective stateWherein, τ is the given threshold that target is present；

(7) prediction is updated respectively without goal hypothesis intensity f_t+1|t(φ) and forecast sample collectionObtain the t+1 moment Without goal hypothesis intensity f_t+1|t+1(φ) and sample set

(8) following parameter is obtained according to the result of step (7)：

8a) to the sample set after renewalIn intensity summation, obtain the number of targets desired value that the t+1 moment estimates η_t+1；

8c) using without goal hypothesis intensity f_t+1|t+1(φ), without goal hypothesis transition probability ψ_t+1(φ | φ) and forecast sample number N_t+1|tCalculate newborn sample number J_t+1；

Many base least square cross bearingses 9a) are utilized, all measurements at t and t+1 moment coordinate are converted into, respectively table It is shown asWithFilter out all satisfactionsCoordinate, constitute new Raw coordinates of targets collection：Wherein,WithI-th of the measurement conversion of t and t+1 moment is represented respectively Coordinate,WithJ-th of newborn coordinates of targets of t and t+1 moment is represented respectively, and T represents transposition, ε_minAnd ε_maxRespectively Represent apart from lower limit and apart from the upper limit, N_zFor all coordinate numbers for meeting condition；

The speed of newborn target 9b) is calculated according to the coordinate of newborn target：Wherein, Δ T is represented Time interval；

(10) according to newborn coordinates of targets collectionWith speed collectionConstitute newborn object setSentence Whether disconnected new life object set is empty, if it is empty, then in target region [a_x,b_x]*[a_y,b_y] interior uniform sampling J_t+1Individual target Sample；Otherwise, each newborn target is distributed as according to the state of newborn target to produceIndividual target sample, constitutes newborn target Sample setWherein,Represent newborn target strength；

(11) t=t+1 is made, return to step (2) continues to track target.

2. according to the method described in claim 1, wherein：According to forecast sample collection in the step (5)And t The measurement collection at+1 momentPredict the measurement intensity at t+1 momentCarry out as follows：

<mrow> <msub> <mi>&lambda;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mover> <mi>&lambda;</mi> <mo>~</mo> </mover> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&zeta;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

Wherein, j ∈ [1, N_z,t+1], N_z,t+1The measurement number at t+1 moment is represented,Represent measurement source in dbjective state Intensity,Measurement source is represented in the intensity without goal hypothesis, its calculation formula is as follows：

<mrow> <msub> <mover> <mi>&lambda;</mi> <mo>~</mo> </mover> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <msub> <mi>f</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

<mrow> <msub> <mi>&zeta;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> </munderover> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> <msubsup> <mi>w</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>;</mo> </mrow>

Wherein,Represent to measure the likelihood for coming from no goal hypothesis,Indicate that the detection of no goal hypothesis is general Rate,Represent to measure likelihood,Represent target detection probability,WithMeter Calculate formula as follows：

<mrow> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>&pi;</mi> <mo>|</mo> <mi>&Omega;</mi> <msup> <mo>|</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>-</mo> <mi>h</mi> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msup> <mi>&Omega;</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>(</mo> <mrow> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>-</mo> <mi>h</mi> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

Wherein,！Factorial computing is represented, Ω represents covariance matrix.

3. according to the method described in claim 1, wherein, by measurement intensity in the step (6)With measurement likelihoodCalculate the probability W that target is present_j, its calculation formula is as follows：

<mrow> <msub> <mi>W</mi> <mi>j</mi> </msub> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> </munderover> <mfrac> <mrow> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&lambda;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&CenterDot;</mo> <msubsup> <mi>w</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>,</mo> </mrow>

Wherein,Represent to measure likelihood,Represent target detection probability,Calculating it is public Formula is as follows：

4. according to the method described in claim 1, wherein updating prediction in the step (7) respectively without goal hypothesis intensity f_t+1|t (φ) and forecast sample collectionObtain the t+1 moment without goal hypothesis intensity f_t+1|t+1(φ) and sample setCarry out as follows：

(7.1) no goal hypothesis intensity f is calculated as follows_t+1|t+1(φ)：

<mrow> <msub> <mi>f</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&lsqb;</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> <mo>)</mo> </mrow> <mo>+</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </munderover> <mfrac> <mrow> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&lambda;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&rsqb;</mo> <mo>&CenterDot;</mo> <msub> <mi>f</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>,</mo> </mrow>

Wherein,The detection probability of no goal hypothesis is indicated,Represent to measure and come from no goal hypothesis seemingly So, its calculation formula is as follows：

(7.2) intensity after i-th of target update of t+1 moment is calculated as follows

<mrow> <msubsup> <mover> <mi>w</mi> <mo>^</mo> </mover> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>=</mo> <mo>&lsqb;</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> <mo>)</mo> </mrow> <mo>+</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </munderover> <mfrac> <mrow> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> <msubsup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>D</mi> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&lambda;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&rsqb;</mo> <mo>&CenterDot;</mo> <msubsup> <mi>w</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>,</mo> </mrow>

Wherein,Represent target detection probability,Represent to measure likelihood, its calculation formula is as follows：

<mrow> <msub> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>|</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>&pi;</mi> <mo>|</mo> <mi>&Omega;</mi> <msup> <mo>|</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>-</mo> <mi>h</mi> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msup> <mi>&Omega;</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>(</mo> <mrow> <msubsup> <mi>z</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>j</mi> </msubsup> <mo>-</mo> <mi>h</mi> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> </mrow> <mi>i</mi> </msubsup> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

Wherein, Ω represents covariance matrix.

5. according to the method described in claim 1, wherein the step 8b) it is middle using without goal hypothesis intensity f_t+1|t+1(φ) and Forecast sample number N_t+1|t, calculate resampling sample number N_t+1|t+1, carry out as follows：

8b1) according to number of targets desired value η_t+1With without goal hypothesis intensity f_t+1|t+1(φ) calculates resampling ratio p_s：

8b2) according to adopting ratio p again_sWith forecast sample number N_t+1|tCalculate resampling sample number N_t+1|t+1：

N_t+1|t+1=N_t+1|t·p_s。

6. according to the method described in claim 1, wherein the step 8c) it is middle using without goal hypothesis intensity f_t+1|t+1(φ), nothing Goal hypothesis transition probability ψ_t+1(φ | φ) and forecast sample number N_t+1|tCalculate newborn sample number J_t+1, its calculation formula is as follows：

<mrow> <msub> <mi>J</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mo>&lsqb;</mo> <mfrac> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <msub> <mi>&eta;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mfrac> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>&psi;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>(</mo> <mrow> <mi>&phi;</mi> <mo>|</mo> <mi>&phi;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msub> <mi>f</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>&phi;</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>,</mo> </mrow>

Wherein η_t+1The number of targets desired value estimated for the t+1 moment.

7. according to the method described in claim 1, wherein the step 9a) it is middle using many base least square cross bearingses, by t All measurements at moment and t+1 moment are converted into coordinate, carry out as follows：

9a1) set initial state value as：According to j-th of measurementInitial target functional value R (X are calculated as follows^k):

Wherein, j ∈ [1, N_z,t+1],x_t+1And y_t+1X directions and y directions coordinate are represented respectively, and k represents iteration Number of times, nonlinear functionIt is calculated as follows：

Wherein, (x_si,y_si) represent the observation station position, i=1 ..., n；

Jacobi matrixes g 9a2) is calculated as follows^k：

Hessian matrixes G 9a3) is calculated as follows^k：

An iteration state value 9a4) is calculated as follows

An iteration target function value 9a5) is calculated as follows

9a6) judgeWithIt is whether equal, ifThen export j-th of measurementThe seat of conversion Mark：Otherwise, makeReturn to step (9a2), is changed next time Generation；Wherein,

9a7) press step 9a1)~9a6), all measurements of t and t+1 moment are converted into coordinate, are expressed as：With