CN113484866A - Multi-target detection tracking method based on passive sonar azimuth history map - Google Patents

Abstract

A multi-target detection tracking method based on a passive sonar azimuth course map belongs to the technical field of underwater multi-target tracking. The problem of current underwater target tracking method tracking performance poor is solved. The method comprises the steps of screening a target initial track by adopting a tracking wave gate to obtain a real target track, and determining a target tracking method according to the number of point-in tracks in each real target tracking wave gate at each sampling moment; if the true target tracks are intersected, an MHT algorithm is adopted to measure in the tracking wave gate to form a hypothesis event, the probability of the hypothesis event is calculated, and a state updating equation of the target is obtained; if the real target tracks do not intersect, the measuring point tracks of each target are respectively processed by adopting a PDA algorithm, the measuring point tracks are weighted by utilizing the association probability of all the measuring point tracks in the wave gate and the target to obtain a state updating equation of the target, and the target position is tracked by utilizing the combination of the corresponding target state updating equation and the MPUKF filtering technology. The underwater multi-target tracking method is suitable for underwater multi-target tracking.

Description

Multi-target detection tracking method based on passive sonar azimuth history map

Technical Field

The invention belongs to the technical field of underwater multi-target tracking, and particularly relates to a multi-target detection tracking method based on a passive sonar azimuth course map.

Background

With the gradual deepening of the understanding of human beings on the ocean, various countries are vigorously developing and exploring the ocean field, wherein the detection and tracking of underwater targets are an important task in ocean defense war. The passive sonar detection system receives noise radiated by a target, extracts information such as target direction and the like after processing, and completes tracking of the target through data processing. Because the target can be associated, positioned and identified under the condition that the azimuth information is accurate, the underwater target pure azimuth tracking technology is particularly important in a passive sonar detection system.

Target detection is a tracking condition, and aims to judge whether a target exists or not and provide a data source for subsequent tracking. The N-P criterion is a detection criterion commonly used in the field of target detection, and the idea is that the false alarm probability P is used_faObtaining the decision threshold under a certain condition to make the detection probability P_DThe maximum value is reached, the N-P criterion can reach the optimal detection performance under the fixed background noise intensity, but the problem that the N-P criterion cannot adapt to the changing environment exists.

Data association is the most complex problem in multi-target tracking and is the basic premise for realizing multi-target tracking. The correctness of data association processing directly affects the tracking performance and the track quality, and the wrong data association can cause the situations of target tracking error, target tracking loss or tracking filtering divergence and the like. Under ideal conditions, MHT is considered to be the optimal method for processing data association and is widely adopted in a tracking system, but MHT generates a large number of low probability hypothesis events and has a large number of sequencing operations.

Disclosure of Invention

The invention aims to solve the problem of poor tracking performance of the existing underwater target tracking method, and provides a multi-target detection tracking method based on a passive sonar azimuth course map.

The invention relates to a multi-target detection tracking method based on a passive sonar azimuth course map, which comprises the following steps:

processing sonar observation signals by adopting a beam forming technology to obtain a multi-target azimuth history map, using the multi-target azimuth history map, taking a maximum value point of an azimuth spectrum of each sampling moment in the azimuth history map as a detection unit, and acquiring a self-adaptive noise background of each detection unit;

detecting the detection units by using the N-P criterion detector and the self-adaptive noise background of each detection unit to obtain a measuring point track, and acquiring a plurality of target initial tracks by using the measuring point tracks at the initial two moments;

step three, screening a plurality of target initial tracks by adopting a tracking wave gate to obtain a plurality of real target tracks, and determining a target tracking method according to the number of measuring point tracks falling into the current sampling point moment in each real target tracking wave gate;

if only one measuring point trace falls into the real target tracking wave gate, filtering the measuring point trace in the tracking wave gate by adopting an MPUKF filter; realizing target azimuth detection and tracking;

if the plurality of tracks fall into the true target tracking wave gate, judging whether the tracking wave gate is intersected, if so, intersecting the true target track, otherwise, not intersecting the true target track;

if the true target tracks are intersected, an MHT algorithm is adopted to measure in the tracking wave gate to form a hypothesis event, the probability of the hypothesis event is calculated, and a state updating equation of the target is obtained; executing the step four;

if the real target tracks are not intersected, respectively processing the measuring point track of each target by adopting a PDA algorithm, weighting the measuring point track by utilizing the association probability of all measuring point tracks in the wave gate and the target, and obtaining a state updating equation of the target by taking the weighted sum as an equivalent result; executing the step four;

fourthly, tracking the target position by combining a corresponding target state update equation and an MPUKF filtering technology according to whether the target tracks are intersected or not;

and step five, in the tracking process, when the tracking wave gate does not have the measuring point track within A moments, judging that the target track is terminated, and finishing the detection and tracking of the target position, wherein A is a positive integer larger than 3.

Further, in the second step of the present invention, a specific method for obtaining multiple target initial tracks by using the measurement point tracks at the initial two moments comprises:

step A1, detecting the detection unit at each sampling time by using the self-adaptive noise background of each detection unit by adopting an N-P standard detector, judging whether the detection unit at each sampling time has a measurement trace, and if so, executing step A2;

and A2, obtaining a target initial track according to the measuring point tracks detected at the initial two moments.

Further, in the present invention, the specific method for determining whether the detecting unit has the measurement trace in step a1 includes:

step A11, calculating the level estimation value of sonar measurement environment background noise power

Step A12, calculating the power of the detecting unit

Step A13, utilizing the estimated value

And a known false alarm probability P_faObtaining a detection threshold D_T；

Step A14, detecting the power of the unit

And a detection threshold D_TBy comparison, when

When there isAnd measuring the trace points.

Further, in the present invention, the specific method for obtaining the target initial track in step a2 is as follows:

step A21, taking all measuring point tracks detected at the initial moment as track starting points, and establishing an initial orientation wave gate by taking each measuring point track as a center;

step a22, by formula:

|z(2)-z(1)|≤2θ_BW

judging whether a measuring point trace falls into each initial azimuth wave gate at the second moment, if yes, acquiring an initial flight path according to the number and the distance of the point traces in the initial azimuth wave gate, and if no, deleting the corresponding flight path initial point, wherein theta is equal to the number of the point traces in the initial azimuth wave gate, and if not, deleting the corresponding flight path initial point_BW2arcsin (lambda/md) is the half-power spot beam width of the main lobe of the beam, lambda is the signal wavelength, m is the number of array elements, d is the distance between the array elements, and z (k) (k is 1,2) is the measured trace detected at the initial time and the second time.

Further, in the present invention, the specific method for obtaining the real target track in the third step is as follows:

b1, initializing the unscented Kalman filter of the modified polar coordinate system by using the target initial track, predicting the target initial track by using the unscented Kalman filter of the modified polar coordinate system after initialization to obtain a predicted value of the next sampling point moment, and setting a tracking wave gate by taking the predicted value as the center;

and step B2, judging whether the measuring point trace at the third to fifth sampling points falls into the corresponding tracking wave gate, if the measuring point trace at least one sampling point falls into the corresponding tracking wave gate, the target initial track is the real target track.

Further, in the third step, if a plurality of traces fall into the real target tracking wave gate, the specific method for judging whether the tracking wave gates intersect is as follows:

judging the measurement predicted value z of the target t at the kth sampling point moment^t(k | k-1) (t ═ 1,2) satisfies:

if so, then the gates intersect, otherwise, the gates do not intersect, where z¹(k | k-1) is the measured and predicted value of target 1 at the time of k sampling point, z²And (k | k-1) is a measured predicted value of the target 2 at the time of k sampling points, and gamma is a wave gate parameter.

Further, in the present invention, if the real target tracks intersect in step two, the specific method for obtaining the state update equation of the target is as follows:

step C1, generating a current sampling point time hypothesis event according to the previous sampling point time hypothesis event and the current measurement set z (k); the hypothetical events are: the number of the measurement point tracks in the tracking wave gate belonging to the real target track is assumed as follows:

the number of measuring point traces belonging to the false flight path in the wave gate is as follows:

φ＝m_k-τ

wherein tau is the number of measuring point tracks belonging to the real target track in the tracking wave gate, and tau_iIs the ith measuring point track belonging to the real target track, phi is the number of the measuring point tracks belonging to the false track in the wave gate, m_kMeasuring the total number of traces for the trace gate internal quantity;

step C2, calculating the probability of the assumption of the establishment of the event in the step C1, and obtaining the probability of each measuring point trace in the tracking wave gate being associated with two crossed targets according to the probability of the assumption of the establishment of the event;

step C3, summing the probabilities of each measuring point trace in the tracking wave gate being associated with the two crossed targets to obtain the associated probability of each measuring point trace and the two crossed targets;

step C4, obtaining a state updating equation of the target corresponding to the tracking wave gate according to the association probability of each measuring point trace and the two crossed targets:

in the formula: v^t(k) Is a combined innovation interconnected to the target t; v_i ^t(k) Is the combined innovation of the ith measurement point trace at the moment k and the target t; x^t(k | k) is the state update value for target t at time k;

for the ith measurement point trace z at time k_i(k) Probability of interconnection with target t; x^t(k | k-1) is the predicted value of target t at the sampling point time of state k; k^t(k) Is the filter gain of the target t.

Further, in the second step of the present invention, if the real target tracks do not intersect, the specific method for obtaining the state update equation of the target is as follows:

d1, calculating the probability that all measuring point tracks in the tracking wave gate belong to the real target track;

ith measurement point trace z at time k sample point_i(k) Probability of belonging to an object beta_i(k) Comprises the following steps:

wherein mu is the space density of the false measurement, namely the false measurement number of unit volume; v_i(k) Is the ith quantity at the time of the k sample pointMeasuring innovation, V, corresponding to trace_i' (k) is V_i(k) S (k) is the innovation covariance at the time of k sample points, S^-1(k) Is the inverse of S (k), P_DThe detection probability P of the target track when the real target track does not cross_GMeasuring the probability of the trace of the point falling into the wave gate when the real target flight path does not intersect;

step D2, obtaining a state updating equation of the target by utilizing the probability that all measuring point tracks in the tracking wave gate belong to the true target track:

wherein the content of the first and second substances,

to combine innovation, V_i(k) Measuring the innovation of the trace for the ith measurement point at the k sampling point moment;

X_i(k|k)＝X(k|k-1)+K(k)·V_i(k)

X_i(k | k) is the time of the sampling point k at event θ_i(k) A state estimate of the target of the condition; x (k | k-1) is the prediction of the target state at the sampling point time k by the state equation at the sampling point time k-1, and when i is 0, the predicted value is used as the estimated value, namely X₀(k | k) ═ X (k | k-1); k (k) is the filter gain at the time of the k sample points.

Further, in the third step, the specific method for tracking the target azimuth by combining the target state update equation and the MPUKF filtering technology is as follows:

e1, establishing a state equation and a measurement equation of the target under a polar coordinate system;

e2, obtaining a predicted value X (k | k-1) of the target state at the moment k and a measured predicted value z (k | k-1) by using a state equation and a measurement equation;

e3, obtaining the innovation and the covariance of the k sampling point by using the measurement predicted value z (k | k-1) at the k sampling point and the measurement value z (k) observed at the k sampling point;

and E4, obtaining an estimated value X (k | k) of the target state at the sampling point moment of k by using the target state predicted value X (k | k-1) and the information at the sampling point moment of k and through a target state updating equation, and realizing the tracking of the target azimuth.

Further, in the present invention, in step E4, obtaining an estimated value X (k | k) of the target state at the time of the k sampling point, and implementing a specific method for tracking the target azimuth includes:

establishing a state equation of a target under a polar coordinate system:

wherein the content of the first and second substances,

state vector at time k; x (k-1) is the state vector at time k-1, f [ ·]Is a non-linear state function;

is the azimuth angle at the time of the k sample point,

is the azimuth angle at the time of the k-1 sampling point,

is the rate of change of the azimuth angle at the time of the k sample point,

the azimuth angle change rate at the sampling point time, r (k) is the distance from the target to the observation station at the sampling point time k,

the distance change rate of k sampling points at the moment; r (k-1) is the distance between the target and the observation station at the sampling point moment of k-1,

the distance change rate of the sampling point of k-1 at the moment; t is a measurement time interval;

establishing a measurement equation by using a state equation of a target under a polar coordinate system:

z(k)＝H(k)X(k)+v(k)

wherein h (k) is a measurement matrix of k sample points in time, and h (k) is [ 001 ]; v (k) is the measured noise at the time of k sampling points, which is 0-mean Gaussian white noise;

obtaining a predicted value X (k | k-1) of the target state at the k sampling point time and a covariance P (k | k-1) between the predicted value X (k | k-1) and a state vector X (k) at the k sampling point time by using a state equation:

wherein, Δ X_i(k|k-1)＝ξ_i(k|k-1)-X(k|k-1)，ΔX′_i(k | k-1) is Δ X_iTranspose of (k | k-1); q (k) is the covariance matrix of the process noise; w_iEstimate xi for a sample point_i(k-1| k-1) corresponding weight, n_xIs the dimension of the state vector;

ξ_i(k|k-1)＝f[ξ_i(k-1|k-1)]

is from the sampling point xi at the time of the sampling point k-1_i(k-1| k-1) predicted value, ξ, at time k_i(k-1| k-1) is calculated from the state vector X (k-1| k-1) estimated at the time of the sampling point k-1 and the covariance P (k-1| k-1):

estimate xi of each sample point_iWeight W corresponding to (k-1| k-1)_iComprises the following steps:

in the formula: k is a scale parameter satisfying n_x+κ≠0；

Is (n)_x+ kappa) the ith row or column, n, of the P (k-1| k-1) matrix_xP (k-1| k-1) is the covariance between the estimated value of the target state updated at the time k-1 and the state vector X (k) at the time k sample;

acquiring an estimated value X (k | k) of the target state at the k sampling point time by using a predicted value X (k | k-1) of the target state at the k sampling point time and a covariance P (k | k-1) between the predicted value X (k | k-1) and a state vector X (k) at the k sampling point time:

X(k|k)＝X(k|k-1)+K(k)·V(k)

in the formula:

V(k)＝z(k)-z(k|k-1)

z(k|k-1)＝H(k)·X(k|k-1)

K(k)＝P(k|k-1)·H′(k)·S^-1(k)

S(k)＝H(k)·P(k|k-1)·H′(k)+R(k)

P(k|k)＝P(k|k-1)-K(k)·S(k)·K'(k)

where K (K) is the gain at the time of K sample points, and K' (K) is the transpose of K (K); s (k) is the covariance of innovation; z (k | k-1) is a measured predicted value at the sampling point time of k; h' (k) is the transpose of the measurement matrix H (k); r (k) is a covariance matrix of the measured noise at the time of k sampling points; p (k | k) is the covariance between the estimated value X (k | k) and the true value X (k) of the target state updated at the time of the k sample points.

The invention uses the N-P criterion detector to detect the target azimuth process chart, estimates the local background average power of each detection unit, and further obtains the detection threshold, so that the detection threshold changes along with the change of the environment. The problem that the false alarm probability is increased when the environment is changed due to the fixed detection threshold of the N-P criterion is solved, and the constant false alarm detection is realized. And simultaneously, tracking the targets by using a multi-target pure direction tracking method combining PDA and MHT, in the multi-target tracking process, regarding the targets as a plurality of single targets when no track crossing occurs, tracking and filtering by using a PDA algorithm, and when the track crossing occurs, filtering and tracking the targets with two crossed tracks by using the MHT algorithm. The method improves the convergence speed of tracking.

Drawings

FIG. 1 is a flow chart of a method of the present invention;

FIG. 2 is a multi-objective azimuth history diagram of the sea trial experiment of the present invention;

FIG. 3 shows the results of the sea trial experiment of the present invention involving the use of a detector to detect a target;

FIG. 4 is a diagram illustrating the results of the multi-target tracking method of the present invention for the azimuth tracking of a target;

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

It should be noted that the embodiments and features of the embodiments may be combined with each other without conflict.

The first embodiment is as follows: the following describes this embodiment with reference to fig. 1, and the passive sonar multi-target orientation detection and tracking method according to this embodiment includes:

processing sonar observation signals by adopting a beam forming technology to obtain a multi-target azimuth history map, using the multi-target azimuth history map, taking a maximum value point of an azimuth spectrum of each sampling moment in the azimuth history map as a detection unit, and acquiring a self-adaptive noise background of each detection unit;

detecting the detection units by using the N-P criterion detector and the self-adaptive noise background of each detection unit to obtain a measuring point track, and acquiring a plurality of target initial tracks by using the measuring point tracks at the initial two moments;

step three, screening a plurality of target initial tracks by adopting a tracking wave gate to obtain a plurality of real target tracks, and determining a target tracking method according to the number of measuring point tracks falling into the current sampling point moment in each real target tracking wave gate;

if only one measuring point trace falls into the real target tracking wave gate, filtering the measuring point trace in the tracking wave gate by adopting an MPUKF filter; realizing target azimuth detection and tracking;

if the plurality of tracks fall into the true target tracking wave gate, judging whether the tracking wave gate is intersected, if so, intersecting the true target track, otherwise, not intersecting the true target track;

if the true target tracks are intersected, an MHT algorithm is adopted to measure in the tracking wave gate to form a hypothesis event, the probability of the hypothesis event is calculated, and a state updating equation of the target is obtained; executing the step four;

if the real target tracks are not intersected, respectively processing the measuring point track of each target by adopting a PDA algorithm, weighting the measuring point track by utilizing the association probability of all measuring point tracks in the wave gate and the target, and obtaining a state updating equation of the target by taking the weighted sum as an equivalent result; executing the step four;

fourthly, tracking the target position by combining a corresponding target state update equation and an MPUKF filtering technology according to whether the target tracks are intersected or not;

and step five, in the tracking process, when the tracking wave gate does not have the measuring point track within A moments, judging that the target track is terminated, and finishing the detection and tracking of the target position, wherein A is a positive integer larger than 3. In the tracking process in the embodiment, when the measurement point trace does not appear in the wave gate, the target state is estimated by using the one-step prediction result of the filter to keep tracking the target. If no measuring point trace still appears in the target wave gate after a period of time, the target flight path is really terminated, and the target is tracked.

The specific process for judging whether the target track is terminated is as follows:

setting a threshold number A of continuous non-detection trace moments, wherein the number of the continuous non-detection trace moments is i, and performing the following circulation when i is less than A:

fifthly, setting that no measuring point trace exists at the current time k, and judging whether a measuring point trace exists in a wave gate at the time k + i when i is 1;

if so, filtering and tracking by using the measuring point trace, keeping the target track, and finishing the track ending judgment;

otherwise, estimating the target state by using the one-step prediction value, taking the one-step prediction result as a target tracking filtering result, and meanwhile, the number of continuous non-detection trace moments is i + 1.

Step two, judging whether i is more than A;

if yes, returning to the step four to continue recursion circulation; otherwise, judging that the target track is ended.

Further, in the second embodiment, in the second step, a specific method for obtaining a plurality of target initial tracks by using the measurement point tracks at the initial two times includes:

step A1, detecting the detection unit at each sampling time by using the self-adaptive noise background of each detection unit by adopting an N-P standard detector, judging whether the detection unit at each sampling time has a measurement trace, and if so, executing step A2;

and A2, obtaining a target initial track according to the measuring point tracks detected at the initial two moments. Further, in the present embodiment, the first and second substrates,

the specific method for determining whether the detection unit is a measurement trace in step a1 is as follows:

step A11, calculating the level estimation value of sonar measurement environment background noise power

Step A12, countingCalculating the power of the detection unit

Step A13, utilizing the estimated value

And a known false alarm probability P_faObtaining a detection threshold D_T；

Step A14, detecting the power of the unit

And a detection threshold D_TBy comparison, when

At all times, there is a trace of measurement points.

Further, in the present embodiment, the first and second substrates,

the specific method for obtaining the target initial track in the step a2 is as follows:

step A21, taking all measuring point tracks detected at the initial moment as track starting points, and establishing an initial orientation wave gate by taking each measuring point track as a center;

step a22, by formula:

|z(2)-z(1)|≤2θ_BW

judging whether a measuring point trace falls into each initial azimuth wave gate at the second moment, if yes, acquiring an initial flight path according to the number and the distance of the point traces in the initial azimuth wave gate, and if no, deleting the corresponding starting point of the flight path, wherein theta_BW2arcsin (lambda/md) is the half-power spot beam width of the main lobe of the beam, lambda is the signal wavelength, m is the number of array elements, d is the distance between the array elements, and z (k) (k is 1,2) is the measured trace detected at the initial time and the second time.

If the second moment falls into one trace point in the initial azimuth wave gate, the two trace point at the second moment form an initial track;

and if the second moment falls into the plurality of point tracks in the initial azimuth wave gate, the point track with the closest distance and the starting point are taken to form an initial flight path.

In this embodiment, a multi-target azimuth history map is obtained by using a beam forming technique, the azimuth history map is obtained by recording and displaying target azimuth spectrums at all sampling times, and a maximum point of the azimuth spectrum obtained at each sampling time is used as a detection unit for detection (the detection process is to compare the power of the maximum point with a threshold value). The half-power spot beam width theta of the main lobe is taken as the main lobe width of the formed beam_BWFor detecting the length of the protection unit, i.e. the length of the left and right protection units is theta_BW2, and the left and right lengths are all theta through the protection unit_BWThe reference unit of (2) obtains a level estimate of the background noise power

θ_BWIs represented by the formula:

θ_BW＝2arcsin(λ/md)

and calculating, wherein λ is the wavelength of the received signal, m is the number of array elements, d is the spacing between the array elements, and d is λ/2 for a uniform linear array.

Due to the influence of environmental noise, underwater acoustic channels and beam side lobes, the maximum r is removed when the background of the reference unit is selected₁R of smallest sum₂A noisy background to prevent an uneven background from making the threshold too high or too low. Averaging the residual noise background in the reference unit to obtain the average power of the noise background

The power of the detection unit is the peak power of the power spectrum

The signal-to-noise ratio is:

applying the N-P criterion again by

With a given false alarm probability P_faCan be represented by the formula:

calculating to obtain detection threshold D_TWhere Q (-) is the right tail function of a standard normal distribution.

Will be provided with

And D_TBy comparison, if the formula is satisfied:

if the energy is not the same as the target, the energy is set to be 0.

The SNR can be calculated according to the following equation:

determining the probability of detection P_DIn the formula, H₀Representing the case where the target does not exist, H₁Representing the situation in which the object is present.

In this embodiment, each trace obtained after the detection corresponds to a measurement orientation value,

represents the set of measurement traces measured at time k, where z_i(k) The ith trace point at the time k, and n (k) the number of all trace points at the time k. Regarding all point tracks in z (1) as track starting points, establishing an initial azimuth wave gate by taking each point track as a center, and according to whether the point tracks in z (2) meet the following conditions:

|z(2)-z(1)|≤2θ_BW

judging whether a point track falls into the initial wave gate at the second moment, and if the point track falls into one point track in the initial wave gate at the second moment, forming an initial flight track by twice point tracks; if a plurality of traces fall into the initial wave gate, the closest trace and the initial point are taken to form an initial flight path; if no trace point falls into the initial wave gate, the starting point of the trace is considered as a false measurement, and the starting point is deleted.

Further, in the present embodiment, the first and second substrates,

the specific method for obtaining the real target track in the third step comprises the following steps:

b1, initializing the unscented Kalman filter of the modified polar coordinate system by using the target initial track, predicting the target initial track by using the unscented Kalman filter of the modified polar coordinate system after initialization to obtain a predicted value of the next sampling point moment, and setting a tracking wave gate by taking the predicted value as the center;

and step B2, judging whether the measuring point trace at the third to fifth sampling points falls into the corresponding tracking wave gate, if the measuring point trace at least one sampling point falls into the corresponding tracking wave gate, the target initial track is the real target track.

In this embodiment, the kalman filter is initialized by forming two point traces of the initial trace, a predicted value of the target at the next time is predicted, and a tracking gate is set with the predicted value as the center, that is, whether the measurement z (k) of the target satisfies:

[z(k)-z(k|k-1)]′·S^-1(k)·[z(k)-z(k|k-1)]≤γ

wherein gamma is a wave gate parameter; p (k | k-1) is a predicted value; s (k) is innovation covariance at time k:

S(k)＝H(k)·P(k|k-1)·H′(k)+R(k)

obtaining the change result of the target azimuth measurement information along with the time from the time azimuth process chart, wherein the measurement point trace obtained by detection at a certain moment appears on a linear area, and the tracking wave gate is taken to take the predicted value as the center and the length as the length

And the measured value z (k) is n_z1-dimensional, the measurement dimension n is given in the table below_zWhen 1, the gate probabilities P corresponding to different parameters lambda_G。

TABLE 1n_zProbability of falling into the wave gate at 1_G

The number of starting time of the flight path is determined by the number of targets, relative positions of the targets, detection probability and false alarm probability, and the starting step number N is 5 in specific implementation.

If a point track falls into the wave gate at the next moment, the initial track can be considered as a real target track; if no trace point falls into the wave gate, the one-step predicted value is used for replacing the measured value at the moment to perform prediction again, and judgment is repeated for two or three times. If the trace falls into the wave gate in the two subsequent moments, the trace is also determined as the target track; if no trace point falls into the wave gate at the two subsequent moments, the flight path is considered as a false flight path and is cancelled.

Further, in the present embodiment, the first and second substrates,

in the third step, if a plurality of traces fall into the real target tracking wave gate, the specific method for judging whether the tracking wave gates are intersected is as follows:

judging the measurement predicted value z of the target t at the kth sampling point moment^t(k | k-1) (t ═ 1,2) satisfies:

if so, then the gates intersect, otherwise, the gates do not intersect, where z¹(k | k-1) is the measured and predicted value of target 1 at the time of k sampling point, z²And (k | k-1) is a measured predicted value of the target 2 at the time of k sampling points, and gamma is a wave gate parameter.

Further, in the present embodiment, the first and second substrates,

in the second step, if the true target track intersects, the specific method for obtaining the state update equation of the target is as follows:

step C1, generating a current sampling point time hypothesis event according to the previous sampling point time hypothesis event and the current measurement set z (k); the hypothetical events are: the number of the measurement point tracks in the tracking wave gate belonging to the real target track is assumed as follows:

the number of measuring point traces belonging to the false flight path in the wave gate is as follows:

φ＝m_k-τ

wherein tau is the number of measuring point tracks belonging to the real target track in the tracking wave gate, and tau_iIs the ith measuring point track belonging to the real target track, phi is the number of the measuring point tracks belonging to the false track in the wave gate, m_kMeasuring the total number of traces for the trace gate internal quantity;

step C2, calculating the probability of the assumption of the establishment of the event in the step C1, and obtaining the probability of each measuring point trace in the tracking wave gate being associated with two crossed targets according to the probability of the assumption of the establishment of the event;

step C3, summing the probabilities of each measuring point trace in the tracking wave gate being associated with the two crossed targets to obtain the associated probability of each measuring point trace and the two crossed targets;

step C4, obtaining a state updating equation of the target corresponding to the tracking wave gate according to the association probability of each measuring point trace and the two crossed targets:

in the formula: v^t(k) Is a combined innovation interconnected to the target t;

is the combined innovation of the ith measurement point trace at the moment k and the target t; x^t(k | k) is the state update value for target t at time k;

for the ith measurement point trace z at time k_i(k) Probability of interconnection with target t; x^t(k | k-1) is the predicted value of target t at the sampling point time of state k; k^t(k) Is the filter gain of the target t.

In the invention, only the situation when two target tracks are crossed is processed by MHT calculation, so the point track z_i(k) The assumption of (1) is that the target (1) and (2) tracks are continuous or false, and each target can only be interconnected with one measurement or no measurement regardless of the generation of new targets.

If m hypothetical events Θ occur at the end of time k-1^k-1,s(s-1, …, m), n possible hypotheses are generated at time k, and a total of mn hypothetical events Θ results at the end of time k^k,l(l＝1,…,mn)。

Let event Θ (k) about the current measurement include: tau measurements from established tracks; phi false measurements. For i ═ 1, …, m_kThe following labeled variables relating to event Θ (k) are defined

Then the measurements from the established track:

the false measurement number is as follows:

φ＝m_k-τ

calculating the hypothesis probability:

for any hypothetical event Θ^k,lThe probability and hypothesis event Θ can be derived^k-1,sThe probability calculation formula of the recursive relationship is as follows:

in the formula: c ═ P { z (k) | z^k-1Is a normalization constant factor;

the detection probabilities of the two track targets; n is a radical of_t[z_i(k)]Indicating that the metrology associated with target t follows a gaussian distribution:

N_t[z_i(k)]＝N[z_i(k)；z_t(k|k-1)S_t(k)]

wherein z is_t(k | k-1) is a predicted measure of target t; s_t(k) Is its corresponding innovation covariance.

Trace of measurement points z_i(k) Probability of interconnection with target t (t ═ 1,2)

Is all the assumed events theta at the end of time k^k,lThe ith measurement z_i(k) The sum of the assumed event probabilities from target t.

State update and covariance update:

the state updating and covariance updating of the MHT algorithm are implemented by the following expressions of a state updating equation for a single target t:

in the formula: v^t(k) Is a combined innovation interconnected to the target t; x^t(k | k) is the state estimate update value for target t at time k; x^t(k | k-1) is the state prediction value of target t; k^t(k) Is the filter gain of the target t.

The covariance update equation for target t is:

in the formula: p is a radical of^t(k | k) is the covariance update value of target t at time k;

p^t(k | k-1) is a covariance one-step predicted value of the target t;

probability of a hypothetical event associated with target t for no measurements;

when the MHT algorithm is applied to filter updating, the covariance updating equation and the state updating equation of the target t are adopted for calculation,

and (3) adopting a tolerant method to perform hypothesis deletion to solve the problem of MHT hypothesis event combination explosion, namely reserving P (P > 0) hypotheses with high probability, wherein P is equal to 5.

Further, in the present embodiment, the first and second substrates,

in the third step, if the real target tracks do not intersect, the specific method for obtaining the state update equation of the target is as follows:

d1, calculating the probability that all measuring point tracks in the tracking wave gate belong to the real target track;

k samplingIth measurement point trace z at a point in time_i(k) Probability of belonging to an object beta_i(k) Comprises the following steps:

wherein mu is the space density of the false measurement, namely the false measurement number of unit volume; v_i(k) Is an innovation, V, corresponding to the ith metrology trace at the time of the k sample point_i' (k) is V_i(k) S (k) is the innovation covariance at the time of k sample points, S^-1(k) Is the inverse of S (k), P_DThe detection probability P of the target track when the real target track does not cross_GMeasuring the probability of the trace of the point falling into the wave gate when the real target flight path does not intersect;

step D2, obtaining a state updating equation of the target by utilizing the probability that all measuring point tracks in the tracking wave gate belong to the true target track:

wherein the content of the first and second substances,

to combine innovation, V_i(k) Measuring the innovation of the trace for the ith measurement point at the k sampling point moment;

X_i(k|k)＝X(k|k-1)+K(k)·V_i(k)

X_i(k | k) is the time of the sampling point k at event θ_i(k) A state estimate of the target of the condition; x (k | k-1) is the number of times of k sampling points obtained by the state equation of the k-1 sampling point timeFor the prediction of the target state, when i is 0, the predicted value is used as the estimated value, i.e. X₀(k | k) ═ X (k | k-1); k (k) is the filter gain at the time of the k sample points.

In this embodiment, the process of applying the PDA algorithm to update the filtering of the target state when the target gates are not intersected is as follows:

first, the probability that all the measurement point traces in the wave gate come from the target is calculated:

by using

Representing a set of point traces falling into the wave gate at the moment k, wherein m (k) is the number of the point traces falling into the wave gate; z is a radical of^kRepresenting the cumulative collection of acknowledgment traces up to time k.

Defining an event:

θ_i(k) denotes z_i(k) Is an event from the correct measurement of the target;

θ₀(k) indicating that time k does not have a measurement from the target.

Ith point trace z at time k_i(k) The probability of interconnecting with the target is:

β_i(k)＝P{θ_i(k)|z^k}

in practical situations, when there are more false measurements in the wave gate or the false alarm rate is high, the number of false measurements falling into the wave gate is approximately obeyed by the parameter μ V_kThe interconnection probability at this time is:

in the formula (I), the compound is shown in the specification,

wherein μ is the spatial density of the spurious measurement, i.e. the number of spurious measurements per unit volume, and μ can be estimated by averaging the number of traces falling into the wave gate until k:

wherein, V_kIs the wave gate volume, when n_zWhen the number is equal to 1, the alloy is put into a container,

state update and covariance update:

considering all measurement point traces in the wave gate, assigning a certain association probability to each measurement point trace, and performing weighted summation to obtain the state of the final target at the time k:

in the formula, X_i(k | k) is an event θ_i(k) State estimation for the condition, namely:

X_i(k|k)＝X(k|k-1)+K(k)·V_i(k)

in the formula, V_i(k) Is an innovation corresponding to the measured value

If no measurement is associated with the target, i.e. i is 0

X₀(k|k)＝X(k|k-1)

The available target state update equation:

in the formula (I), the compound is shown in the specification,

to combine innovation.

The corresponding covariance update is:

in the formula:

P_c(k|k)＝P(k|k-1)-P(k|k-1)·K(k)·H(k)

further, in the present embodiment, the first and second substrates,

in the fourth step, the specific method for tracking the target azimuth by combining the target state update equation and the MPUKF filtering technology is as follows:

e1, establishing a state equation and a measurement equation of the target under a polar coordinate system;

e2, obtaining a predicted value X (k | k-1) of the target state at the moment k and a measured predicted value z (k | k-1) by using a state equation and a measurement equation;

e3, obtaining the innovation and the covariance of the k sampling point by using the measurement predicted value z (k | k-1) at the k sampling point and the measurement value z (k) observed at the k sampling point;

and E4, obtaining an estimated value X (k | k) of the target state at the sampling point moment of k by using the target state predicted value X (k | k-1) and the information at the sampling point moment of k and through a target state updating equation, and realizing the tracking of the target azimuth.

Further, in the present embodiment, the first and second substrates,

in step E4, obtaining an estimated value X (k | k) of the target state at the time of the k sampling point, and implementing a specific method for tracking the target azimuth is as follows:

establishing a state equation of a target under a polar coordinate system:

wherein the content of the first and second substances,

state vector at time k; x (k-1) is the state vector at time k-1, f [ ·]Is a non-linear state function;

is the azimuth angle at the time of the k sample point,

is the azimuth angle at the time of the k-1 sampling point,

is the rate of change of the azimuth angle at the time of the k sample point,

the azimuth angle change rate at the sampling point time, r (k) is the distance from the target to the observation station at the sampling point time k,

the distance change rate of k sampling points at the moment; r (k-1) is the distance between the target and the observation station at the sampling point moment of k-1,

the distance change rate of the sampling point of k-1 at the moment; t is a measurement time interval;

establishing a measurement equation by using a state equation of a target under a polar coordinate system:

z(k)＝H(k)X(k)+v(k)

wherein h (k) is a measurement matrix of k sample points in time, and h (k) is [ 001 ]; v (k) is the measured noise at the time of k sampling points, which is 0-mean Gaussian white noise;

obtaining a predicted value X (k | k-1) of the target state at the k sampling point time and a covariance P (k | k-1) between the predicted value X (k | k-1) and a state vector X (k) at the k sampling point time by using a state equation:

wherein, Δ X_i(k|k-1)＝ξ_i(k|k-1)-X(k|k-1)，ΔX′_i(k | k-1) is Δ X_iTranspose of (k | k-1); q (k) is the covariance matrix of the process noise; w_iEstimate xi for a sample point_i(k-1| k-1) corresponding weight, n_xIs the dimension of the state vector;

ξ_i(k|k-1)＝f[ξ_i(k-1|k-1)]

is from the sampling point xi at the time of the sampling point k-1_i(k-1| k-1) predicted value, ξ, at time k_i(k-1| k-1) is calculated from the state vector X (k-1| k-1) estimated at the time of the sampling point k-1 and the covariance P (k-1| k-1):

estimate xi of each sample point_iWeight W corresponding to (k-1| k-1)_iComprises the following steps:

in the formula: k is a scale parameter satisfying n_x+κ≠0；

Is (n)_x+ kappa) the ith row or column, n, of the P (k-1| k-1) matrix_xP (k-1| k-1) is the covariance between the estimated value of the target state updated at the time k-1 and the state vector X (k) at the time k sample;

acquiring an estimated value X (k | k) of the target state at the k sampling point time by using a predicted value X (k | k-1) of the target state at the k sampling point time and a covariance P (k | k-1) between the predicted value X (k | k-1) and a state vector X (k) at the k sampling point time:

X(k|k)＝X(k|k-1)+K(k)·V(k)

in the formula:

V(k)＝z(k)-z(k|k-1)

z(k|k-1)＝H(k)·X(k|k-1)

K(k)＝P(k|k-1)·H′(k)·S^-1(k)

S(k)＝H(k)·P(k|k-1)·H′(k)+R(k)

P(k|k)＝P(k|k-1)-K(k)·S(k)·K'(k)

where K (K) is the gain at the time of K sample points, and K' (K) is the transpose of K (K); s (k) is covariance of innovation; z (k | k-1) is a measured predicted value at the sampling point time of k; h' (k) is the transpose of the measurement matrix H (k); r (k) is a covariance matrix of the measured noise at the time of k sampling points; p (k | k) is the covariance between the estimated value X (k | k) and the true value X (k) of the target state updated at the time of the k sample points.

The specific embodiment is as follows: sea trial experiment: the description is made with reference to fig. 2 to 4;

and (3) selecting 25-minute time domain signals from data acquired by the 32-element uniform linear array, and tracking the target at 750 steps, wherein the time of each step is 2 seconds. The time domain signals are processed by using a beam forming technology to obtain an azimuth history chart as shown in fig. 2.

Then, a detector is used for detecting the target, the false alarm probability is set to be 0.1%, and the lengths of the left protection unit, the right protection unit and the reference unit are both 1 degree. After detection, the energy of the target position is judged to be unchanged, and the energy of the target position is judged to be 0. The energy values are taken in decibels (dB) and normalized to obtain the detection results as shown in fig. 3.

Substituting the signal-to-noise ratio, the number of array elements, the spacing between the array elements, the snapshot number, the approximate direction of the target and other information of the signals in the experimental data into the Clarithromol bound, calculating the estimated Clarithromol bound of the direction, and adjusting the standard deviation of the measured noise in the filter and the gate probability P according to the calculation result_GThe track ending process takes a as 5, and the tracking result is shown in fig. 4.

From the results in the figure, it can be seen that the multi-target tracking method can be used for effectively and purely tracking the multi-target in the environment. The method can better process the situations of multi-target track intersection with equivalent energy, multi-target track intersection with obvious energy difference, target track discontinuity, target track ending and the like.

Although the invention herein has been described with reference to particular embodiments, it is to be understood that these embodiments are merely illustrative of the principles and applications of the present invention. It is therefore to be understood that numerous modifications may be made to the illustrative embodiments and that other arrangements may be devised without departing from the spirit and scope of the present invention as defined by the appended claims. It should be understood that features described in different dependent claims and herein may be combined in ways different from those described in the original claims. It is also to be understood that features described in connection with individual embodiments may be used in other described embodiments.

Claims (10)

1. A multi-target detection tracking method based on a passive sonar azimuth course map is characterized by comprising the following steps:

processing sonar observation signals by adopting a beam forming technology to obtain a multi-target azimuth history map, using the multi-target azimuth history map, taking a maximum value point of an azimuth spectrum of each sampling moment in the azimuth history map as a detection unit, and acquiring a self-adaptive noise background of each detection unit;

detecting the detection units by using the self-adaptive noise background of each detection unit by using an N-P criterion detector to obtain a measuring point track, and acquiring a plurality of target initial tracks by using the measuring point tracks at the initial two moments;

step three, screening a plurality of target initial tracks by adopting a tracking wave gate to obtain a plurality of real target tracks, and determining a target tracking method according to the number of measuring point tracks falling into the current sampling point moment in each real target tracking wave gate;

if only one measuring point trace falls into the real target tracking wave gate, filtering the measuring point trace in the tracking wave gate by adopting an MPUKF filter; realizing target azimuth detection and tracking;

if the plurality of tracks fall into the true target tracking wave gate, judging whether the tracking wave gate is intersected, if so, intersecting the true target track, otherwise, not intersecting the true target track;

if the true target tracks are intersected, an MHT algorithm is adopted to measure in the tracking wave gate to form a hypothesis event, the probability of the hypothesis event is calculated, and a state updating equation of the target is obtained; executing the step four;

if the real target tracks are not intersected, respectively processing the measuring point track of each target by adopting a PDA algorithm, weighting the measuring point track by utilizing the association probability of all measuring point tracks in the wave gate and the target, and obtaining a state updating equation of the target by taking the weighted sum as an equivalent result; executing the step four;

fourthly, tracking the target position by combining a corresponding target state update equation and an MPUKF filtering technology according to whether the target tracks are intersected or not;

and step five, in the tracking process, when the tracking wave gate does not have the measuring point track within A moments, judging that the target track is terminated, and finishing the detection and tracking of the target position, wherein A is a positive integer larger than 3.

2. The multi-target detection tracking method based on the passive sonar azimuth history map according to claim 1, characterized in that in step two, the specific method for obtaining a plurality of target initial tracks by using the measuring point tracks at the initial two moments is:

step A1, detecting the detection unit at each sampling time by using the self-adaptive noise background of each detection unit by adopting an N-P criterion detector, judging whether the detection unit at each sampling time has a measurement trace, and if so, executing step A2;

and A2, obtaining a target initial track according to the measuring point tracks detected at the initial two moments.

3. The multi-target detection tracking method based on the passive sonar azimuth course map according to claim 1, characterized in that the specific method for judging whether the detection unit is a measurement point trace in step a1 is as follows:

step A11, calculating the level estimation value of sonar measurement environment background noise power

Step A12, calculating the power of the detecting unit

Step A13, utilizing the estimated value

And a known false alarm probability P_faObtaining a detection threshold D_T；

Step A14, detecting the power of the unit

And a detection threshold D_TBy comparison, when

At all times, there is a trace of measurement points.

4. The multi-target detection tracking method based on the passive sonar azimuth course map according to claim 2, characterized in that the specific method for obtaining the target initial track in step a2 is as follows:

step A21, taking all measuring point tracks detected at the initial moment as track starting points, and establishing an initial orientation wave gate by taking each measuring point track as a center;

step a22, by formula:

|z(2)-z(1)|≤2θ_BW

judging whether a measuring point trace falls into each initial azimuth wave gate at the second moment, if yes, acquiring an initial flight path according to the number and the distance of the point traces in the initial azimuth wave gate, and if no, deleting the corresponding flight path initial point, wherein theta is equal to the number of the point traces in the initial azimuth wave gate, and if not, deleting the corresponding flight path initial point_BW2arcsin (lambda/md) is the half-power spot beam width of the main lobe of the beam, lambda is the signal wavelength, m is the number of array elements, d is the spacing of the array elements, z: (m)k) (k is 1,2) is the measurement trace detected at the initial time and the second time.

5. The multi-target detection tracking method based on the passive sonar azimuth history map according to claim 2, characterized in that the specific method for obtaining the true target track in the third step is:

b1, initializing the unscented Kalman filter of the modified polar coordinate system by using the target initial track, predicting the target initial track by using the unscented Kalman filter of the modified polar coordinate system after initialization to obtain a predicted value of the next sampling point moment, and setting a tracking wave gate by taking the predicted value as the center;

and step B2, judging whether the measuring point trace at the third to fifth sampling points falls into the corresponding tracking wave gate, if the measuring point trace at least one sampling point falls into the corresponding tracking wave gate, the target initial track is the real target track.

6. The multi-target detection tracking method based on the passive sonar azimuth history map according to claim 5, wherein in the third step, if a plurality of traces fall in a real target tracking wave gate, the specific method for judging whether the tracking wave gate intersects is as follows:

judging the measurement predicted value z of the target t at the kth sampling point moment^t(k | k-1) (t ═ 1,2) satisfies:

if so, then the gates intersect, otherwise, the gates do not intersect, where z¹(k | k-1) is the measured and predicted value of target 1 at the time of k sampling point, z²And (k | k-1) is a measured predicted value of the target 2 at the time of k sampling points, and gamma is a wave gate parameter.

7. The method for detecting and tracking the multiple targets based on the passive sonar azimuth history map according to claim 6, wherein in the third step, if the true target tracks intersect, the specific method for obtaining the state update equation of the target is as follows:

step C1, generating a current sampling point time hypothesis event according to the previous sampling point time hypothesis event and the current measurement set z (k); the hypothetical events are: the number of the measurement point tracks in the tracking wave gate belonging to the real target track is assumed as follows:

the number of measuring point traces belonging to the false flight path in the wave gate is as follows:

φ＝m_k-τ

wherein tau is the number of measuring point tracks belonging to the real target track in the tracking wave gate, and tau_iIs the ith measuring point track belonging to the real target track, phi is the number of the measuring point tracks belonging to the false track in the wave gate, m_kMeasuring the total number of traces for the trace gate internal quantity;

step C2, calculating the probability of the assumption of the establishment of the event in the step C1, and obtaining the probability of each measuring point trace in the tracking wave gate being associated with two crossed targets according to the probability of the assumption of the establishment of the event;

step C3, summing the probabilities of each measuring point trace in the tracking wave gate being associated with the two crossed targets to obtain the associated probability of each measuring point trace and the two crossed targets;

step C4, obtaining a state updating equation of the target corresponding to the tracking wave gate according to the association probability of each measuring point trace and the two crossed targets:

in the formula: v^t(k) Is interconnected with the meshCombining innovation of mark t; v_i ^t(k) Is the combined innovation of the ith measurement point trace at the moment k and the target t; x^t(k | k) is the state update value for target t at time k;

for the ith measurement point trace z at time k_i(k) Probability of interconnection with target t; x^t(k | k-1) is the predicted value of target t at the sampling point time of state k; k^t(k) Is the filter gain of the target t.

8. The multi-target detection tracking method based on the passive sonar azimuth history map according to claim 7, wherein in the third step, if the true target tracks do not intersect, the specific method for obtaining the state update equation of the target is as follows:

d1, calculating the probability that all measuring point tracks in the tracking wave gate belong to the real target track;

ith measurement point trace z at time k sample point_i(k) Probability of belonging to an object beta_i(k) Comprises the following steps:

wherein mu is the space density of the false measurement, namely the false measurement number of unit volume; v_i(k) Is an innovation, V, corresponding to the ith metrology trace at the time of the k sample point_i' (k) is V_i(k) S (k) is the innovation covariance at the time of k sample points, S^-1(k) Is the inverse of S (k), P_DTarget track when real target track does not crossProbability of detection of, P_GMeasuring the probability of the trace of the point falling into the wave gate when the real target flight path does not intersect;

step D2, obtaining a state updating equation of the target by utilizing the probability that all measuring point tracks in the tracking wave gate belong to the true target track:

wherein the content of the first and second substances,

to combine innovation, V_i(k) Measuring the innovation of the trace for the ith measurement point at the k sampling point moment;

X_i(k|k)＝X(k|k-1)+K(k)·V_i(k)

X_i(k | k) is the time of the sampling point k at event θ_i(k) A state estimate of the target of the condition; x (k | k-1) is the prediction of the target state at the sampling point time k by the state equation at the sampling point time k-1, and when i is 0, the predicted value is used as the estimated value, namely X₀(k | k) ═ X (k | k-1); k (k) is the filter gain at the time of the k sample points.

9. The method for detecting and tracking the multiple targets based on the passive sonar azimuth history map according to claim 8, wherein in the fourth step, the specific method for tracking the azimuth of the target by combining the state update equation of the target and the MPUKF filtering technology is as follows:

e1, establishing a state equation and a measurement equation of the target under a polar coordinate system;

e2, obtaining a predicted value X (k | k-1) of the target state at the moment k and a measured predicted value z (k | k-1) by using a state equation and a measurement equation;

e3, obtaining the innovation and the covariance of the k sampling point by using the measurement predicted value z (k | k-1) at the k sampling point and the measurement value z (k) observed at the k sampling point;

and E4, obtaining an estimated value X (k | k) of the target state at the sampling point moment of k by using the target state predicted value X (k | k-1) and the information at the sampling point moment of k and through a target state updating equation, and realizing the tracking of the target azimuth.

10. The multi-target detection tracking method based on the passive sonar azimuth process map according to claim 9, wherein in step E4, an estimated value X (k | k) of a target state at a sampling point time k is obtained, and a specific method for tracking the target azimuth is as follows:

establishing a state equation of a target under a polar coordinate system:

wherein the content of the first and second substances,

state vector at time k; x (k-1) is the state vector at time k-1, f [ ·]Is a non-linear state function;

is the azimuth angle at the time of the k sample point,

is the azimuth angle at the time of the k-1 sampling point,

is the rate of change of the azimuth angle at the time of the k sample point,

at the moment of the sampling pointThe azimuth change rate, r (k), is the distance of the target from the observation station at the time of k sampling points,

the distance change rate of k sampling points at the moment; r (k-1) is the distance between the target and the observation station at the sampling point moment of k-1,

the distance change rate of the sampling point of k-1 at the moment; t is a measurement time interval;

establishing a measurement equation by using a state equation of a target under a polar coordinate system:

z(k)＝H(k)X(k)+v(k)

wherein h (k) is a measurement matrix of k sample points in time, and h (k) is [ 001 ]; v (k) is the measured noise at the time of k sampling points, which is 0-mean Gaussian white noise;

obtaining a predicted value X (k | k-1) of the target state at the k sampling point time and a covariance P (k | k-1) between the predicted value X (k | k-1) and a state vector X (k) at the k sampling point time by using a state equation:

wherein, Δ X_i(k|k-1)＝ξ_i(k|k-1)-X(k|k-1)，ΔX′_i(k | k-1) is Δ X_iTranspose of (k | k-1); q (k) is the covariance matrix of the process noise; w_iEstimate xi for a sample point_i(k-1| k-1) corresponding weight, n_xIs the dimension of the state vector;

ξ_i(k|k-1)＝f[ξ_i(k-1|k-1)]

is from the sampling point xi at the time of the sampling point k-1_i(k-1| k-1) predicted value, ξ, at time k_i(k-1| k-1) is estimated from the time of the k-1 sample pointThe calculated state vector X (k-1| k-1) and covariance P (k-1| k-1) are calculated as:

estimate xi of each sample point_iWeight W corresponding to (k-1| k-1)_iComprises the following steps:

in the formula: k is a scale parameter satisfying n_x+κ≠0；

Is (n)_x+ kappa) the ith row or column, n, of the P (k-1| k-1) matrix_xP (k-1| k-1) is the covariance between the estimated value of the target state updated at the time k-1 and the state vector X (k) at the time k sample;

acquiring an estimated value X (k | k) of the target state at the k sampling point time by using a predicted value X (k | k-1) of the target state at the k sampling point time and a covariance P (k | k-1) between the predicted value X (k | k-1) and a state vector X (k) at the k sampling point time:

X(k|k)＝X(k|k-1)+K(k)·V(k)

in the formula:

V(k)＝z(k)-z(k|k-1)

z(k|k-1)＝H(k)·X(k|k-1)

K(k)＝P(k|k-1)·H′(k)·S^-1(k)

S(k)＝H(k)·P(k|k-1)·H′(k)+R(k)

P(k|k)＝P(k|k-1)-K(k)·S(k)·K'(k)

where K (K) is the gain at the time of K sample points, and K' (K) is the transpose of K (K); s (k) is covariance of innovation; z (k | k-1) is a measured predicted value at the sampling point time of k; h' (k) is the transpose of the measurement matrix H (k); r (k) is a covariance matrix of the measured noise at the time of k sampling points; p (k | k) is the covariance between the estimated value X (k | k) and the true value X (k) of the target state updated at the time of the k sample points.

Images