CN111948657A - Maneuvering weak target tracking-before-detection method based on multimode particle filtering - Google Patents

Maneuvering weak target tracking-before-detection method based on multimode particle filtering Download PDF

Info

Publication number
CN111948657A
CN111948657A CN202010738626.8A CN202010738626A CN111948657A CN 111948657 A CN111948657 A CN 111948657A CN 202010738626 A CN202010738626 A CN 202010738626A CN 111948657 A CN111948657 A CN 111948657A
Authority
CN
China
Prior art keywords
target
state
particle
tracking
following formula
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010738626.8A
Other languages
Chinese (zh)
Other versions
CN111948657B (en
Inventor
齐滨
王燕
付进
梁国龙
李想
王逸林
张光普
邹男
王晋晋
邱龙皓
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN202010738626.8A priority Critical patent/CN111948657B/en
Publication of CN111948657A publication Critical patent/CN111948657A/en
Application granted granted Critical
Publication of CN111948657B publication Critical patent/CN111948657B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/66—Sonar tracking systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/539—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention relates to a maneuvering weak target tracking-before-detection method based on multimode particle filtering. The invention belongs to the technical field of underwater target tracking, and the method comprises the steps of carrying out parameter initialization processing and determining a receiving signal of a passive sonar array; processing by adopting a broadband conventional beam forming algorithm according to the receiving signals of the passive sonar array to obtain a space spectrum, and taking the space spectrum as measurement data; according to the measured data, noise balance judges the suspicious target at the current moment; partitioning the target state space according to the measurement data; sampling each particle state of each target according to the target state space partition result, and calculating a weight; resampling particles of the same target independently; estimating the state of the target according to the sampling result; and when the target duration exceeds the joint observation frame number, performing joint judgment on the target, and deleting target information which does not pass the judgment. The invention realizes the real-time tracking of a plurality of maneuvering targets and realizes the detection and tracking of maneuvering weak targets in a passive sonar scene.

Description

Maneuvering weak target tracking-before-detection method based on multimode particle filtering
Technical Field
The invention relates to the technical field of underwater target tracking, in particular to a maneuvering weak target tracking-before-detection method based on multimode particle filtering.
Background
Due to the self characteristics of the passive sonar, a received signal often has the characteristic of low signal to noise ratio, most of the traditional detection-before-tracking technology uses the threshold point track direction as measurement information, but the target is easy to miss detection under the condition of low signal to noise ratio due to threshold processing, so the tracking effect of the low signal to noise ratio is poor. In addition, for non-cooperative targets, the motion state of the non-cooperative targets can be changed at any time according to the battle mission, and the non-cooperative targets have certain mobility, so that the difficulty of detection and tracking of the passive sonar is further increased. Therefore, the passive sonar can effectively detect and track the maneuvering target with low signal-to-noise ratio, and has important significance for improving the combat efficiency and the survival capability of the passive sonar.
The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering can effectively realize the detection and tracking of the maneuvering target in a scene with a low signal-to-noise ratio, but at present, related researches are mostly based on a radar system and limited in a single-target scene, different targets are seriously interfered with each other in a multi-target scene, and particularly when newly generated target particles are sampled to be nearby other targets, particles which carry wrong target information and have large weights are generated, so that the tracking result is greatly influenced.
Disclosure of Invention
The invention provides a maneuvering weak target pre-detection tracking method based on multimode particle filtering, which aims to solve the tracking problem of maneuvering targets with low signal-to-noise ratios in a passive sonar scene, and realize effective detection and tracking of a plurality of maneuvering targets with low signal-to-noise ratios in the passive sonar scene by limiting the newly generated particle state space and weakening the problem of target interference during sampling, and the invention provides the following technical scheme:
a maneuvering weak target tracking-before-detection method based on multimode particle filtering comprises the following steps:
step 1: performing parameter initialization processing to determine a receiving signal of the passive sonar array;
step 2: processing by adopting a broadband conventional beam forming algorithm according to the receiving signals of the passive sonar array to obtain a space spectrum, and taking the space spectrum as measurement data;
and step 3: according to the measured data, noise balance judges the suspicious target at the current moment;
and 4, step 4: partitioning the target state space according to the measurement data;
and 5: sampling each particle state of each target according to the target state space partitioning result, and determining a weight;
step 6: resampling particles of the same target independently;
and 7: estimating the state of the target according to the sampling result;
and 8: when the target duration time exceeds the joint observation frame number, performing joint judgment on the target, and deleting target information which does not pass the judgment;
and step 9: and repeating the steps 2 to 8 for real-time tracking.
Preferably, the step 1 specifically comprises: when the number of the simulation scene targets is T, the initial state X of the target is represented by the following formula0：
X0＝[x1,x2,…,xT]
Figure BDA0002606004020000021
Wherein x istRepresents the initial state of the target t, thetatWhich represents the orientation of the target t,
Figure BDA0002606004020000022
the angular velocity of the target t is represented,
Figure BDA0002606004020000023
representing the angular acceleration of the target t, EtA presence state variable representing a target t, rtA motion model variable representing the target t;
the passive sonar array type is a uniform line array, the array model is initialized, the array model comprises array elements M and array element spacing d, and the receiving signals of the passive sonar array are initialized, and the receiving signals comprise sampling frequency fsSignal band fbandThe reception signal y (t) of the passive sonar array is expressed by the following equation:
y(t)＝[y1(t),y2(t),…,yM(t)]T
wherein, ym(t) denotes a received signal of an array element M, where M is 1,2, …, M.
Preferably, the step 2 specifically comprises: processing the received signal of the passive sonar array by utilizing a broadband conventional beam forming algorithm, and taking the obtained space spectrum as measurement data Z ═ Z1,z2,…,zB]And dividing the entire observation range into [ theta ]1,θ2,…,θB]，zbIs an azimuth thetabCorresponding spatial spectral energy amplitude, where B is 1,2, …, B.
Preferably, the step 3 specifically comprises: setting an orientation threshold LambdaθWhen the absolute value of the difference between the suspicious target at the moment k and the predicted azimuth of the target existing at the moment k-1 is smaller than the azimuth threshold lambdaθIf the target is judged to be the existing target, otherwise, the target is judged to be the new target, the prediction direction is calculated by the formula, and the noise balance method comprises the following steps:
when the spatial spectrum at a certain time is measured as z ═ { z (1), z (2), …, z (b) }, corresponding to [ theta ] respectively1,θ1,…,θB]For an angle thetaiSelecting a data window with the window length of 2K +1, wherein the data in the window is as follows:
{ z (i-K), z (i-K +1), …, z (i), …, z (i + K-1), z (i + K) }, ordering data within a window from small to large: { v (1), v (2), … v (2K +1) }, median filtering is performed on the measurement values z (i) to be processed, the median filtered measurement values being represented by:
Figure BDA0002606004020000024
wherein z is0(i) Is the filtered angle thetaiMeasured value of ΛzFor the threshold value, the threshold value is calculated by the following formula:
Figure BDA0002606004020000031
wherein α is a threshold adjustment parameter, and a suspicious target existence sequence is represented by the following formula:
Figure BDA0002606004020000032
wherein e (i) ═ 1 represents an angle θiWith suspicious object, E (i) 0 represents angle θiNo suspicious object.
Preferably, the step 4 specifically includes:
step 4.1: predicting the position of the target j at the current time by using the motion state of the target j at the previous time, and expressing the position of the target j at the current time by the following formula
Figure BDA0002606004020000033
Figure BDA0002606004020000034
A＝[10000]
Wherein the content of the first and second substances,
Figure BDA0002606004020000035
representing the state transition matrix corresponding to target j at time k-1,
Figure BDA0002606004020000036
and (3) representing motion model variables of a plurality of particles of the target j at the moment of k-1, and A representing a target orientation extraction matrix.
Step 4.2: spatial spectrum peak information after azimuth prediction and noise balance screening and target j historical multi-frame amplitude information are used as prior information to realize the mutual correlation of the spatial spectrum peak and the target, limit the state space of the target new particle azimuth at the moment k, and express the state space at the spatial spectrum peak through the following formula
Figure BDA0002606004020000037
Figure BDA0002606004020000038
Wherein the content of the first and second substances,
Figure BDA0002606004020000039
is the spatial spectrum peak orientation, θ, corresponding to target jbandThe size of the target azimuth state space is determined for the azimuth bandwidth.
Preferably, the step 5 specifically comprises:
step 5.1: according to the target j, j is 1,2, …, Tk，TkIf j is 1 for the suspicious target number at the current time, and the current time k is the new time of the target j, initializing the new particle state, and expressing the state initialization process by the following formula:
let i equal 1, the particle orientation initialization is represented by:
Figure BDA00026060040200000310
wherein rand represents the generation of a random number over the [0,1] interval;
performing particle angular velocity initialization, wherein the particle angular velocity initialization is represented by the following formula:
Figure BDA0002606004020000041
performing particle angular acceleration initialization, which is expressed by the following formula:
Figure BDA0002606004020000042
according to the target initial birth probability U0Performing target presence state variable initialization, which is represented by the following formula:
Figure BDA0002606004020000043
initializing target motion model variables according to the target initial motion model probability, and expressing the initialization of the target motion model variables by the following formula:
Figure BDA0002606004020000044
calculating the weight of the particles, and expressing the weight of the particles by the following formula:
Figure BDA0002606004020000045
the likelihood function model of the target and the clutter is obtained by fitting a measurement sample, and fitting different space spectrum energy amplitudes and corresponding likelihood model parameters to obtain a function of likelihood function model parameters related to the space spectrum energy amplitudes, and the likelihood function model parameters are adjusted according to the space spectrum energy amplitudes by using the function related to the space spectrum energy amplitudes to distinguish the target space spectrum and the clutter space spectrum within a certain energy amplitude range;
making i equal to i +1, and when i is less than or equal to N, carrying out particle state initialization again;
step 5.2: and if the current time k is greater than the new generation time of the target j, predicting the particle state, and setting i to be 1 to obtain a particle existing state variable at the current time: when in use
Figure BDA0002606004020000046
And rand is less than or equal to U0Then, then
Figure BDA0002606004020000047
Otherwise
Figure BDA0002606004020000048
When in use
Figure BDA0002606004020000049
And rand is less than or equal to 1-U0Then, then
Figure BDA00026060040200000410
Otherwise
Figure BDA00026060040200000411
For all the current time
Figure BDA00026060040200000412
The particles of (a) are considered in two cases, when the particles are nascent particles,
Figure BDA00026060040200000413
then the particle state initialization is carried out again; when the particles are persistent particles
Figure BDA00026060040200000414
Then pass through ΠmAnd the previous time particle motion model variable
Figure BDA00026060040200000415
Obtaining the current particle motion model variable
Figure BDA00026060040200000416
According to IImRandomly selecting and keeping the current motion model according to the set initial transition probability and the probability of keeping the current state
Figure BDA00026060040200000417
And predicting the state of the particles at the current moment by using a corresponding motion equation, and when the model is a uniform angular velocity motion model, expressing a state transition matrix of the uniform angular velocity motion model by using the following formula:
Figure BDA0002606004020000051
when the model is a uniform angular acceleration motion model, a state transition matrix of the uniform angular acceleration motion model is represented by the following formula:
Figure BDA0002606004020000052
is determined by
Figure BDA0002606004020000053
The corresponding equation of motion:
Figure BDA0002606004020000054
wherein v iskThe noise covariance matrix is Q for process noise;
calculating a particle weight; let i equal to i +1, when i is less than or equal to N, then the particle state is predicted again.
Preferably, the target state estimation result is expressed by the following formula
Figure BDA0002606004020000055
Figure BDA0002606004020000056
Preferably, the step 8 specifically includes:
when the combined observation frame number N is exceeded from the target start to the current momentbAnd then, judging the target by utilizing the joint likelihood ratio, and representing a judgment process by the following formula:
Figure BDA0002606004020000057
wherein, ΛlRepresenting a likelihood ratio decision threshold, E k1 represents the presence of an object, E k0 represents the absence of the target, E1And E0Respectively representing the hypothesis that the suspicious target is a real target and the hypothesis that the suspicious target is a clutter, namely judging that the target exists when the suspicious target exceeds a threshold, otherwise judging that the target does not exist, and deleting the target information which does not pass the judgment;
when j < TkLet j equal j +1 and return to step 5.
Preferably, the joint observation frame number is set to 3-4 frames.
The invention has the following beneficial effects:
the invention directly uses the space spectrum as the measurement of the algorithm, and can fully utilize the data information. In practice, the signal-to-noise ratio of a passive sonar target is unknown, so that a function of a likelihood function parameter and a space spectrum energy amplitude value is obtained through fitting, a likelihood function model parameter is adjusted according to the space spectrum energy amplitude value to enable the likelihood function model parameter to be capable of distinguishing a target space spectrum and a clutter space spectrum within a certain energy amplitude value range, then tracking of a maneuvering target is achieved by matching with a multi-mode particle filtering algorithm, a state space partition sampling method is provided, the influence of sampling to the vicinity of other targets on a tracking result is weakened, real-time tracking of a plurality of maneuvering targets is achieved, and finally starting and stopping of the target are achieved through a joint likelihood ratio, so that detection and tracking of the maneuvering weak target under a passive sonar scene can be achieved through the algorithm under the condition.
Drawings
FIG. 1 is a flow chart of a method for tracking a maneuvering weak target before detection based on multi-mode particle filtering;
FIG. 2 is a diagram of the movement situation of the target and the matrix;
FIG. 3 is a graph of tracking results;
FIG. 4 is a diagram of tracking results of a conventional tracking method;
FIG. 5 is a graph of the estimated target number at different times for a condition of 100 Monte Carlo times;
FIG. 6 is a graph comparing the detection rates of the present method and the conventional tracking method under different SNR conditions of the received signal.
Detailed Description
The present invention will be described in detail with reference to specific examples.
The first embodiment is as follows:
as shown in fig. 1, the present invention provides a method for tracking a maneuvering weak target before detection based on multimode particle filtering, which specifically comprises:
a maneuvering weak target tracking-before-detection method based on multimode particle filtering comprises the following steps:
step 1: performing parameter initialization processing to determine a receiving signal of the passive sonar array;
the step 1 specifically comprises the following steps: when the number of the simulation scene targets is T, the initial state X of the target is represented by the following formula0：
X0＝[x1,x2,…,xT]
Figure BDA0002606004020000061
Wherein x istRepresents the initial state of the target t, thetatWhich represents the orientation of the target t,
Figure BDA0002606004020000062
the angular velocity of the target t is represented,
Figure BDA0002606004020000063
representing the angular acceleration of the target t, EtRepresenting the presence state variable of the target t, rt representing the motion model variable of the target t;
the passive sonar matrix is a uniform line arrayInitializing an array model, including an array element M and an array element distance d, initializing a receiving signal of the passive sonar array, including a sampling frequency fsSignal band fbandThe reception signal y (t) of the passive sonar array is expressed by the following equation:
y(t)＝[y1(t),y2(t),…,yM(t)]T
wherein, ym(t) denotes a received signal of an array element M, where M is 1,2, …, M.
Step 2: processing by adopting a broadband conventional beam forming algorithm according to the receiving signals of the passive sonar array to obtain a space spectrum, and taking the space spectrum as measurement data;
the step 2 specifically comprises the following steps: processing the received signal of the passive sonar array by utilizing a broadband conventional beam forming algorithm, and taking the obtained space spectrum as measurement data Z ═ Z1,z2,…,zB]And dividing the entire observation range into [ theta ]1,θ2,…,θB]，zbIs an azimuth thetabCorresponding spatial spectral energy amplitude, where B is 1,2, …, B.
And step 3: according to the measured data, noise balance judges the suspicious target at the current moment;
the step 3 specifically comprises the following steps: setting an orientation threshold LambdaθWhen the absolute value of the difference between the suspicious target at the moment k and the predicted azimuth of the target existing at the moment k-1 is smaller than the azimuth threshold lambdaθIf the target is judged to be the existing target, otherwise, the target is judged to be the new target, the prediction direction is calculated by the formula, and the noise balance method comprises the following steps:
when the spatial spectrum at a certain time is measured as z ═ { z (1), z (2), …, z (b) }, corresponding to [ theta ] respectively1,θ1,…,θB]For an angle thetaiSelecting a data window with the window length of 2K +1, wherein the data in the window is as follows:
{ z (i-K), z (i-K +1), …, z (i), …, z (i + K-1), z (i + K) }, ordering data within a window from small to large: { v (1), v (2), … v (2K +1) }, median filtering is performed on the measurement values z (i) to be processed, the median filtered measurement values being represented by:
Figure BDA0002606004020000071
wherein z is0(i) Is the filtered angle thetaiMeasured value of ΛzFor the threshold value, the threshold value is calculated by the following formula:
Figure BDA0002606004020000072
wherein α is a threshold adjustment parameter, and a suspicious target existence sequence is represented by the following formula:
Figure BDA0002606004020000073
wherein e (i) ═ 1 represents an angle θiWith suspicious object, E (i) 0 represents angle θiNo suspicious object.
And 4, step 4: partitioning the target state space according to the measurement data;
the step 4 specifically comprises the following steps:
step 4.1: predicting the position of the target j at the current time by using the motion state of the target j at the previous time, and expressing the position of the target j at the current time by the following formula
Figure BDA0002606004020000081
Figure BDA0002606004020000082
A＝[10000]
Wherein the content of the first and second substances,
Figure BDA0002606004020000083
representing the state transition matrix corresponding to target j at time k-1,
Figure BDA0002606004020000084
indicating the k-1 timeThe motion model variables of a plurality of particles of the target j, A represents a target orientation extraction matrix;
step 4.2: spatial spectrum peak information after azimuth prediction and noise balance screening and target j historical multi-frame amplitude information are used as prior information to realize the mutual correlation of the spatial spectrum peak and the target, limit the state space of the target new particle azimuth at the moment k, and express the state space near the spatial spectrum peak through the following formula
Figure BDA0002606004020000085
Figure BDA0002606004020000086
Wherein the content of the first and second substances,
Figure BDA0002606004020000087
is the spatial spectrum peak orientation, θ, corresponding to target jbandThe size of the target azimuth state space is determined for the azimuth bandwidth.
And 5: sampling each particle state of each target according to the target state space partitioning result, and determining a weight;
the step 5 specifically comprises the following steps:
step 5.1: according to the target j, j is 1,2, …, Tk，TkIf j is 1 for the suspicious target number at the current time, and the current time k is the new time of the target j, initializing the new particle state, and expressing the state initialization process by the following formula:
let i equal 1, the particle orientation initialization is represented by:
Figure BDA0002606004020000088
wherein rand represents the generation of a random number over the [0,1] interval;
performing particle angular velocity initialization, wherein the particle angular velocity initialization is represented by the following formula:
Figure BDA0002606004020000089
performing particle angular acceleration initialization, which is expressed by the following formula:
Figure BDA00026060040200000810
according to the target initial birth probability U0Performing target presence state variable initialization, which is represented by the following formula:
Figure BDA0002606004020000091
initializing target motion model variables according to the target initial motion model probability, and expressing the initialization of the target motion model variables by the following formula:
Figure BDA0002606004020000092
calculating the weight of the particles, and expressing the weight of the particles by the following formula:
Figure BDA0002606004020000093
the likelihood function model of the target and the clutter is obtained by fitting a measurement sample, and fitting different space spectrum energy amplitudes and corresponding likelihood model parameters to obtain a function of likelihood function model parameters related to the space spectrum energy amplitudes, and the likelihood function model parameters are adjusted according to the space spectrum energy amplitudes by using the function related to the space spectrum energy amplitudes to distinguish the target space spectrum and the clutter space spectrum within a certain energy amplitude range;
making i equal to i +1, and when i is less than or equal to N, carrying out particle state initialization again;
step 5.2: the current time k is larger than the targetAnd marking j new time, predicting the particle state, setting i to be 1, and acquiring a particle existing state variable at the current time: when in use
Figure BDA0002606004020000094
And rand is less than or equal to U0Then, then
Figure BDA0002606004020000095
Otherwise
Figure BDA0002606004020000096
When in use
Figure BDA0002606004020000097
And rand is less than or equal to 1-U0Then, then
Figure BDA0002606004020000098
Otherwise
Figure BDA0002606004020000099
For all the current time
Figure BDA00026060040200000910
The particles of (a) are considered in two cases, when the particles are nascent particles,
Figure BDA00026060040200000911
then the particle state initialization is carried out again; when the particles are persistent particles
Figure BDA00026060040200000912
Then pass through ΠmAnd the previous time particle motion model variable
Figure BDA00026060040200000913
Obtaining the current particle motion model variable
Figure BDA00026060040200000914
According to IImRandomly selecting a probability of initial transition and a probability of preserving a current state of a set preserving a current motion model or transitionOther motion models, and based on
Figure BDA00026060040200000915
And predicting the state of the particles at the current moment by using a corresponding motion equation, and when the model is a uniform angular velocity motion model, expressing a state transition matrix of the uniform angular velocity motion model by using the following formula:
Figure BDA00026060040200000916
when the model is a uniform angular acceleration motion model, a state transition matrix of the uniform angular acceleration motion model is represented by the following formula:
Figure BDA0002606004020000101
is determined by
Figure BDA0002606004020000102
The corresponding equation of motion:
Figure BDA0002606004020000103
wherein v iskThe noise covariance matrix is Q for process noise;
calculating a particle weight; let i equal to i +1, when i is less than or equal to N, then the particle state is predicted again.
Step 6: resampling particles of the same target independently;
and 7: estimating the state of the target according to the sampling result;
the target state estimation result is expressed by the following formula
Figure BDA0002606004020000104
Figure BDA0002606004020000105
And 8: when the target duration time exceeds the joint observation frame number, performing joint judgment on the target, and deleting target information which does not pass the judgment;
the step 8 specifically comprises the following steps:
when the combined observation frame number N is exceeded from the target start to the current momentbAnd then, judging the target by utilizing the joint likelihood ratio, and representing a judgment process by the following formula:
Figure BDA0002606004020000106
wherein, ΛlRepresenting a likelihood ratio decision threshold, E k1 represents the presence of an object, E k0 represents the absence of the target, E1And E0Respectively representing the hypothesis that the suspicious target is a real target and the hypothesis that the suspicious target is a clutter, namely judging that the target exists when the suspicious target exceeds a threshold, otherwise judging that the target does not exist, and deleting the target information which does not pass the judgment;
when j < TkLet j equal j +1 and return to step 5.
Setting a combined observation frame number NbThe method aims to weaken the influence of false alarm caused by clutter with large intensity at a certain moment or missed detection caused by low target intensity, but the algorithm is slow due to the fact that the combined observation frame number is too long, generally 3-4 frames are set, and the combined observation frame number is set to be 3-4 frames.
And step 9: and repeating the steps 2 to 8 for real-time tracking.
The second embodiment is as follows:
firstly, constructing motion tracks of a target and a base array, wherein the base array and the target motion tracks are shown in figure 2, and the base array does uniform linear motion; the target makes a turn motion at a uniform angular velocity for 9-14 seconds, then makes a turn motion at a uniform angular acceleration for 15-22 seconds, and finally makes a turn motion at a uniform angular velocity for 23-42 seconds, and then disappears, and the initial angular velocity omega12 °/s, angular acceleration a1＝0.2°/s2(ii) a The target makes a turning motion at a uniform angular velocity for two 3-9 seconds, then makes a turning motion at a uniform angular acceleration for 10-17 seconds, and finally makes a turning motion at a uniform angular velocity for 18-28 secondsTurning movement, and subsequent disappearance, initial angular velocity ω2Angular acceleration a of-1 °/s2＝-0.3°/s2(ii) a The target makes a turning motion at a uniform angular velocity for three 5-7 seconds, then makes a turning motion at a uniform angular acceleration for 8-15 seconds, finally makes a turning motion at a uniform angular velocity for 16-25 seconds, and then disappears, and the initial angular velocity omega32.1 °/s, angular acceleration a3＝0.05°/s2(ii) a The target makes a uniform angular velocity turning motion for 11-19 seconds, then makes a uniform angular acceleration turning motion for 20-25 seconds, finally makes a uniform angular velocity turning motion for 26-45 seconds, and then disappears, the initial angular velocity omega4Angular acceleration a of 0 DEG/s4＝0.1°/s2(ii) a The target five makes a uniform angular velocity turning motion for 11-25 seconds, then makes a uniform angular acceleration turning motion for 26-30 seconds, finally makes a uniform angular velocity turning motion for 31-43 seconds, and then disappears, the initial angular velocity omega52 °/s, angular acceleration a5＝-1°/s2(ii) a The whole observation time length K is 50 seconds; the observation period dt is 1 second; the number of particles N is 500; target joint observation frame number Nb4 seconds; maximum value of target azimuth thetamax180 deg. and minimum value thetamin0 °; maximum value of angular velocity ωmax5 °/s, minimum value ωmin-5 °/s; maximum value of angular acceleration amax＝3°/s2Minimum value amin＝-3°/s2(ii) a Probability of particle birth U00.2; probability of target initial motion model
Figure BDA0002606004020000111
Target existing state probability transition matrix pie＝[0.8,0.2；0.2,0.8](ii) a Probability transfer matrix pi of target motion modelm＝[0.8,0.2；0.2,0.8](ii) a The corresponding orientations of the spatial spectrum measured values are 0.5 degrees, 1 degree and … 180 degrees; array element number M is 40; the signal-to-noise ratio SNR of the target receiving signal is-23 dB;
finally, the tracking results of the method and the traditional tracking algorithm under the conditions are respectively given, as shown in fig. 3 and fig. 4, respectively, it can be seen that the method can realize maneuvering target tracking under the condition of low signal-to-noise ratio, and the flight path is complete; the traditional tracking method can track targets ( targets 1, 3 and 4) with weak maneuverability, mainly because a motion equation has process noise, the particle cloud can be more dispersed by increasing the process noise, so that part of particles are sampled to be close to a real target, but effective tracking cannot be realized for the targets (targets 2 and 5) with strong maneuverability, and the targets are easy to miss detection due to low signal-to-noise ratio and threshold detection, so that part of tracks are lost. Fig. 5 shows the estimated target number curve of the present method under the monte carlo 100 times condition. Fig. 6 shows the comparison of the detection rates of the method and the conventional tracking method under different received signal-to-noise ratios, which shows that the method has better performance.
In conclusion, the method can realize the detection and tracking of the maneuvering weak target in the passive sonar scene under the condition that the number of the targets is unknown, and has better tracking performance.
The above description is only a preferred embodiment of the method for tracking the maneuvering weak target before detection based on the multimode particle filtering, and the protection range of the method for tracking the maneuvering weak target before detection based on the multimode particle filtering is not limited to the above embodiments, and all technical solutions belonging to the idea belong to the protection range of the invention. It should be noted that modifications and variations which do not depart from the gist of the invention will be those skilled in the art to which the invention pertains and which are intended to be within the scope of the invention.

Claims (9)

1. A maneuvering weak target tracking method before detection based on multimode particle filtering is characterized by comprising the following steps: the method comprises the following steps:
step 1: performing parameter initialization processing to determine a receiving signal of the passive sonar array;
step 2: processing by adopting a broadband conventional beam forming algorithm according to the receiving signals of the passive sonar array to obtain a space spectrum, and taking the space spectrum as measurement data;
and step 3: according to the measured data, noise balance judges the suspicious target at the current moment;
and 4, step 4: partitioning the target state space according to the measurement data;
and 5: sampling each particle state of each target according to the target state space partitioning result, and determining a weight;
step 6: resampling particles of the same target independently;
and 7: estimating the state of the target according to the sampling result;
and 8: when the target duration time exceeds the joint observation frame number, performing joint judgment on the target, and deleting target information which does not pass the judgment;
and step 9: and repeating the steps 2 to 8 for real-time tracking.
2. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the step 1 specifically comprises the following steps: when the number of the simulation scene targets is T, the initial state X of the target is represented by the following formula0：
X0＝[x1,x2,…,xT]
Figure FDA0002606004010000011
Wherein x istRepresents the initial state of the target t, thetatWhich represents the orientation of the target t,
Figure FDA0002606004010000012
the angular velocity of the target t is represented,
Figure FDA0002606004010000013
representing the angular acceleration of the target t, EtA presence state variable representing a target t, rtA motion model variable representing the target t;
the passive sonar array type is a uniform line array, the array model is initialized, the array model comprises array elements M and array element spacing d, and the receiving signals of the passive sonar array are initialized, and the receiving signals comprise sampling frequency fsSignal band fbandThe reception signal y (t) of the passive sonar array is expressed by the following equation:
y(t)＝[y1(t),y2(t),…,yM(t)]T
wherein, ym(t) denotes a received signal of an array element M, where M is 1,2, …, M.
3. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the step 2 specifically comprises the following steps: processing the received signal of the passive sonar array by utilizing a broadband conventional beam forming algorithm, and taking the obtained space spectrum as measurement data Z ═ Z1,z2,…,zB]And dividing the entire observation range into [ theta ]1,θ2,…,θB]，zbIs an azimuth thetabCorresponding spatial spectral energy amplitude, where B is 1,2, …, B.
4. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the step 3 specifically comprises the following steps: setting an orientation threshold LambdaθWhen the absolute value of the difference between the suspicious target at the moment k and the predicted azimuth of the target existing at the moment k-1 is smaller than the azimuth threshold lambdaθIf the target is judged to be the existing target, otherwise, the target is judged to be the new target, the prediction direction is calculated by the formula, and the noise balance method comprises the following steps:
when the spatial spectrum at a certain time is measured as z ═ { z (1), z (2), …, z (b) }, corresponding to [ theta ] respectively1,θ1,…,θB]For an angle thetaiSelecting a data window with the window length of 2K +1, wherein the data in the window is as follows:
{ z (i-K), z (i-K +1), …, z (i), …, z (i + K-1), z (i + K) }, ordering data within a window from small to large: { v (1), v (2), … v (2K +1) }, median filtering is performed on the measurement values z (i) to be processed, the median filtered measurement values being represented by:
Figure FDA0002606004010000021
wherein z is0(i) Is the filtered angle thetaiMeasured value of ΛzFor the threshold value, the threshold value is calculated by the following formula:
Figure FDA0002606004010000022
wherein α is a threshold adjustment parameter, and a suspicious target existence sequence is represented by the following formula:
Figure FDA0002606004010000023
wherein e (i) ═ 1 represents an angle θiWith suspicious object, E (i) 0 represents angle θiNo suspicious object.
5. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the step 4 specifically comprises the following steps:
step 4.1: predicting the position of the target j at the current time by using the motion state of the target j at the previous time, and expressing the position of the target j at the current time by the following formula
Figure FDA0002606004010000024
Figure FDA0002606004010000025
A＝[1 0 0 0 0]
Wherein the content of the first and second substances,
Figure FDA0002606004010000026
representing the state transition matrix corresponding to target j at time k-1,
Figure FDA0002606004010000027
representing motion model variables of a plurality of particles of the target j at the moment of k-1, wherein A represents a target orientation extraction matrix;
step 4.2: spatial spectrum peak information after azimuth prediction and noise balance screening and target j historical multi-frame amplitude information are used as prior information to realize the mutual correlation of the spatial spectrum peak and the target, limit the state space of the target new particle azimuth at the moment k, and express the state space at the spatial spectrum peak through the following formula
Figure FDA0002606004010000031
Figure FDA0002606004010000032
Wherein the content of the first and second substances,
Figure FDA0002606004010000033
is the spatial spectrum peak orientation, θ, corresponding to target jbandThe size of the target azimuth state space is determined for the azimuth bandwidth.
6. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the step 5 specifically comprises the following steps:
step 5.1: according to the target j, j is 1,2, …, Tk，TkIf j is 1 for the suspicious target number at the current time, and the current time k is the new time of the target j, initializing the new particle state, and expressing the state initialization process by the following formula:
let i equal 1, the particle orientation initialization is represented by:
Figure FDA0002606004010000034
wherein rand represents the generation of a random number over the [0,1] interval;
performing particle angular velocity initialization, wherein the particle angular velocity initialization is represented by the following formula:
Figure FDA0002606004010000035
performing particle angular acceleration initialization, which is expressed by the following formula:
Figure FDA0002606004010000036
according to the target initial birth probability U0Performing target presence state variable initialization, which is represented by the following formula:
Figure FDA0002606004010000037
initializing target motion model variables according to the target initial motion model probability, and expressing the initialization of the target motion model variables by the following formula:
Figure FDA0002606004010000038
calculating the weight of the particles, and expressing the weight of the particles by the following formula:
Figure FDA0002606004010000039
the likelihood function model of the target and the clutter is obtained by fitting a measurement sample, and fitting different space spectrum energy amplitudes and corresponding likelihood model parameters to obtain a function of likelihood function model parameters related to the space spectrum energy amplitudes, and the likelihood function model parameters are adjusted according to the space spectrum energy amplitudes by using the function related to the space spectrum energy amplitudes to distinguish the target space spectrum and the clutter space spectrum within a certain energy amplitude range;
making i equal to i +1, and when i is less than or equal to N, carrying out particle state initialization again;
step 5.2: and if the current time k is greater than the new generation time of the target j, predicting the particle state, and setting i to be 1 to obtain a particle existing state variable at the current time: when in use
Figure FDA0002606004010000041
And rand is less than or equal to U0Then, then
Figure FDA0002606004010000042
Otherwise
Figure FDA0002606004010000043
When in use
Figure FDA0002606004010000044
And rand is less than or equal to 1-U0Then, then
Figure FDA0002606004010000045
Otherwise
Figure FDA0002606004010000046
For all the current time
Figure FDA0002606004010000047
The particles of (a) are considered in two cases, when the particles are nascent particles,
Figure FDA0002606004010000048
then the particle state initialization is carried out again; when the particles are persistent particles
Figure FDA0002606004010000049
Then pass through ΠmAnd the previous time particle motion model variable
Figure FDA00026060040100000410
Obtaining the current particle motion model variable
Figure FDA00026060040100000411
According to IImRandomly selecting and keeping the current motion model according to the set initial transition probability and the probability of keeping the current state
Figure FDA00026060040100000412
And predicting the state of the particles at the current moment by using a corresponding motion equation, and when the model is a uniform angular velocity motion model, expressing a state transition matrix of the uniform angular velocity motion model by using the following formula:
Figure FDA00026060040100000413
when the model is a uniform angular acceleration motion model, a state transition matrix of the uniform angular acceleration motion model is represented by the following formula:
Figure FDA00026060040100000414
is determined by
Figure FDA00026060040100000415
The corresponding equation of motion:
Figure FDA00026060040100000416
wherein v iskThe noise covariance matrix is Q for process noise;
calculating a particle weight; let i equal to i +1, when i is less than or equal to N, then the particle state is predicted again.
7. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the target state estimation result is expressed by the following formula
Figure FDA0002606004010000051
Figure FDA0002606004010000052
8. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 1, characterized in that: the step 8 specifically comprises the following steps:
when the combined observation frame number N is exceeded from the target start to the current momentbAnd then, judging the target by utilizing the joint likelihood ratio, and representing a judgment process by the following formula:
Figure FDA0002606004010000053
wherein, ΛlRepresenting a likelihood ratio decision threshold, Ek1 represents the presence of an object, Ek0 represents the absence of the target, E1And E0Respectively representing the hypothesis that the suspicious target is a real target and the hypothesis that the suspicious target is a clutter, namely judging that the target exists when the suspicious target exceeds a threshold, otherwise judging that the target does not exist, and deleting the target information which does not pass the judgment;
when j < TkLet j equal j +1 and return to step 5.
9. The method for tracking the maneuvering weak target before detection based on the multi-mode particle filtering as claimed in claim 8, characterized by: and the combined observation frame number is set to be 3-4 frames.
CN202010738626.8A 2020-07-28 2020-07-28 Maneuvering weak target tracking-before-detection method based on multimode particle filtering Active CN111948657B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010738626.8A CN111948657B (en) 2020-07-28 2020-07-28 Maneuvering weak target tracking-before-detection method based on multimode particle filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010738626.8A CN111948657B (en) 2020-07-28 2020-07-28 Maneuvering weak target tracking-before-detection method based on multimode particle filtering

Publications (2)

Publication Number Publication Date
CN111948657A true CN111948657A (en) 2020-11-17
CN111948657B CN111948657B (en) 2022-08-19

Family

ID=73338346

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010738626.8A Active CN111948657B (en) 2020-07-28 2020-07-28 Maneuvering weak target tracking-before-detection method based on multimode particle filtering

Country Status (1)

Country Link
CN (1) CN111948657B (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112630783A (en) * 2020-11-26 2021-04-09 海鹰企业集团有限责任公司 Passive sonar target tracking method
CN113281729A (en) * 2021-05-31 2021-08-20 中国科学院声学研究所 Target automatic detection method and system based on multi-frame spatial spectrum joint processing
CN113534164A (en) * 2021-05-24 2021-10-22 中船海洋探测技术研究院有限公司 Target path tracking method based on active and passive combined sonar array
CN115616602A (en) * 2022-10-14 2023-01-17 哈尔滨工程大学 Observer optimal maneuvering strategy based on passive sonar pure azimuth positioning pre-detection tracking algorithm

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104820993A (en) * 2015-03-27 2015-08-05 浙江大学 Underwater weak target tracking method combining particle filtering with track before detect
CN107202989A (en) * 2017-05-08 2017-09-26 电子科技大学 A kind of complicated Faint target detection and tracking suitable for passive starboard ambiguity of towed linear array sonar
WO2019071003A1 (en) * 2017-10-06 2019-04-11 Airmar Technology Corporation Aft-looking sonar
CN110187335A (en) * 2019-06-25 2019-08-30 电子科技大学 Tracking before being detected for the particle filter with discontinuous characteristic target

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104820993A (en) * 2015-03-27 2015-08-05 浙江大学 Underwater weak target tracking method combining particle filtering with track before detect
CN107202989A (en) * 2017-05-08 2017-09-26 电子科技大学 A kind of complicated Faint target detection and tracking suitable for passive starboard ambiguity of towed linear array sonar
WO2019071003A1 (en) * 2017-10-06 2019-04-11 Airmar Technology Corporation Aft-looking sonar
CN110187335A (en) * 2019-06-25 2019-08-30 电子科技大学 Tracking before being detected for the particle filter with discontinuous characteristic target

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
WEN CHEN等: "Joint Algorithm Based on Interference Suppression and Kalman Filter for Bearing-Only Weak Target Robust Tracking", 《IEEE ACCESS》 *
徐璐霄: "基于被动阵列的水声弱目标检测前跟踪算法研究", 《中国优秀博硕士学位论文全文数据库（硕士） 工程科技II辑》 *
胡秀华等: "一种新的可变数目机动目标联合检测与跟踪算法", 《中南大学学报(自然科学版)》 *
苗媛媛等: "基于粒子滤波的弱目标检测前跟踪算法", 《计算机系统应用》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112630783A (en) * 2020-11-26 2021-04-09 海鹰企业集团有限责任公司 Passive sonar target tracking method
CN113534164A (en) * 2021-05-24 2021-10-22 中船海洋探测技术研究院有限公司 Target path tracking method based on active and passive combined sonar array
CN113534164B (en) * 2021-05-24 2023-12-12 中船海洋探测技术研究院有限公司 Target path tracking method based on active-passive combined sonar array
CN113281729A (en) * 2021-05-31 2021-08-20 中国科学院声学研究所 Target automatic detection method and system based on multi-frame spatial spectrum joint processing
CN115616602A (en) * 2022-10-14 2023-01-17 哈尔滨工程大学 Observer optimal maneuvering strategy based on passive sonar pure azimuth positioning pre-detection tracking algorithm
CN115616602B (en) * 2022-10-14 2023-08-18 哈尔滨工程大学 Course determination method of observer optimal maneuver strategy based on passive sonar pure azimuth positioning detection pre-tracking algorithm

Also Published As

Publication number Publication date
CN111948657B (en) 2022-08-19

Similar Documents

Publication Publication Date Title
CN111948657B (en) Maneuvering weak target tracking-before-detection method based on multimode particle filtering
CN107450055B (en) High-speed maneuvering target detection method based on discrete linear frequency modulation Fourier transform
CN109188385B (en) Method for detecting high-speed weak target under clutter background
CN109324315B (en) Space-time adaptive radar clutter suppression method based on double-layer block sparsity
WO2003079046A9 (en) An adaptive system and method for radar detection
CN107607916B (en) Self-defense type speed and distance joint deception jamming resisting method
KR100852934B1 (en) System for Detection of Noise Jamming Signal Improving the Side Lobe Blanking and its Method
CN109239703B (en) Real-time tracking method for moving target
CN117008064A (en) Improved YOLOv 5-based radar interference adaptive suppression method
Bagheri et al. A new approach to pulse deinterleaving based on adaptive thresholding
CN111007471B (en) Method for judging interference effect of active suppression interference in simulation environment
CN105652256B (en) A kind of high-frequency ground wave radar TBD methods based on polarization information
CN109507654A (en) Phase information calculation method under a kind of complex environment based on LS
CN112799028A (en) False target identification method based on RCS fluctuation statistical characteristic difference
CN110988867A (en) Elliptical cross target positioning method for one-transmitting and double-receiving through-wall radar
CN111563914A (en) Underwater positioning and tracking method and device and readable storage medium
CN110244289A (en) A kind of adaptive particle filter ground wave radar target integrative detection method
CN113126086B (en) Life detection radar weak target detection method based on state prediction accumulation
CN114325599B (en) Automatic threshold detection method for different environments
CN116299208A (en) Anti-interference method based on active/passive radar composite guide head data association
CN113093174B (en) PHD filter radar fluctuation weak multi-target-based pre-detection tracking method
CN112698295B (en) Knowledge-assisted radar detection and tracking integrated method and system
CN111948613A (en) Ship-borne ground wave radar target detection method based on self-adaptive background area selection
CN111025251B (en) Multi-target composite detection method based on dynamic programming
CN113608178B (en) Anti-drag deception jamming method based on dual-band information fusion

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant