Underwater multi-platform multi-target tracking method based on extended multi-dimensional distribution

Abstract

The invention discloses an underwater multi-platform multi-target tracking method based on extended multi-dimensional distribution, belongs to the field of underwater target tracking, and relates to an underwater multi-platform multi-target tracking method. The objective of the invention is to solve the problem that in the prior art, when azimuth angle measurement values of all platforms are subjected to multi-dimensional distribution and positioning and tracking, the azimuth angle measurement values may be'merged 'due to limitation of conventional beam forming angular resolution, so that the positioning and tracking accuracy is affected. According to the underwater multi-platform multi-target tracking method based on extended multi-dimensional distribution, the normalized intensity of the spatial spectrum is combined, the measurement value which can be effectively split and combined with the lowest total correlation cost as the evaluation criterion is established, and the multi-target positioning and tracking performance is improved according to the distribution result.

Description

Technical Field

[0001] The present invention belongs to the field of underwater target tracking and relates to an underwater multi-platform multi-target tracking method. Background Art

[0002] In an underwater multi-target tracking scenario, integrating measurement information from multiple platforms can make up for the limitations of a single platform in terms of maneuverability, observation range, and accuracy, and has become a key research direction in recent years. In application engineering, the following process is usually adopted to achieve azimuth estimation: First, threshold detection is performed on the azimuth-energy spectrum generated by conventional beamforming, and then the azimuth angle corresponding to the energy peak exceeding the threshold is extracted as the azimuth measurement. However, due to the angular resolution limitation of beamforming, when the azimuth angle interval between co-frequency signal targets in space is less than the main lobe beam width of the array, the problem of "merging" of measurement values will occur. The traditional method for associating the azimuth angles of each platform uses the multi-dimensional assignment method, which takes the lowest total association cost as the objective function. However, the lack of assignment results caused by the "merging" of measurement values will lead to a decline in positioning and tracking performance. Therefore, it is necessary to study reducing the influence of measurement value merging caused by beamforming angular resolution on positioning and tracking performance. Summary of the Invention

[0003] The object of the invention is to solve the problem that in the existing method, when the peak value of conventional beamforming is used as the azimuth angle measurement value and then multi-dimensional assignment and positioning and tracking are performed, the azimuth angle measurement value may be "merged" due to the limitation of the conventional beamforming angular resolution, thereby affecting the positioning and tracking accuracy. A method for underwater multi-platform multi-target tracking based on extended multi-dimensional assignment is proposed.

[0004] The specific process of a method for underwater multi-platform multi-target tracking based on extended multi-dimensional assignment is as follows:

[0005] Step 1: Let the time k = 1;

[0006] Each platform performs conventional beamforming processing, extracts the azimuth angle measurement value of the platform and the corresponding spatial spectrum intensity;

[0007] Step 2: Based on each azimuth angle measurement set Calculate the estimated target position value corresponding to each measurement set; Based on each azimuth angle measurement set And the estimated target position value, calculate the cost function value of each measurement set Based on the cost function value Construct an objective function, solve the objective function to obtain the assignment result; Use the triangular measurement method to estimate the target position for the assignment result; Use the joint probabilistic data association method to track the target position;

[0008] Record the azimuth measurement values of each platform used for each tracking trajectory and the corresponding spatial spectral intensities of the azimuth measurement values;

[0009] Step 3: Determine the size relationship between the azimuth interval of the tracking trajectories of any two targets relative to platform s and the azimuth interval threshold. If the azimuth interval of the tracking trajectories of two targets relative to platform s is less than the azimuth interval threshold, then open the observation window. After platform s opens the observation window, determine whether the azimuth measurement value at time k is within ;

[0010] If the azimuth measurement value at time k is within , then proceed to Step 4;

[0011] If the azimuth measurement value at time k is not within , then let k = k + 1 and execute Step 1;

[0012] where ε θ is the azimuth expansion value of the observation window; represents the azimuth of the tracking trajectory of the i-th target relative to the s-th platform at time k - 1;

[0013] If the azimuth interval of the tracking trajectories of two targets relative to platform s is greater than or equal to the azimuth interval threshold, then do not open the observation window, and let k = k + 1 and execute Step 1;

[0014] Step 4: Calculate the change rate of the spatial spectral intensity based on the spatial spectral intensity corresponding to the azimuth measurement value obtained in Step 2 Compare the change rate of the spatial spectral intensity with the intensity change threshold Γ A . If the change rate of the spatial spectral intensity is greater than the intensity change threshold Γ A , then the corresponding measurement value is used as a "suspicious merging measurement value"; if the change rate of the spatial spectral intensity is less than or equal to the intensity change threshold Γ A , then let k = k + 1 and execute Step 1;

[0015] Step 5: Split the "suspicious merging measurement value", recombine the split measurement values, calculate the total association cost of each recombined measurement value, select the combination with the minimum total association cost as the final allocation result, estimate the target position based on the allocation result, and track the target position.

[0016] Preferably, in Step 1, each platform performs conventional beamforming processing to extract the azimuth measurement value of the platform and the corresponding spatial spectral intensity; the specific process is as follows:

[0017] Step 1: Assume that each platform is equipped with a uniform linear array of M hydrophones, and the signal vector received by the linear array is expressed as

[0018]

[0019] where,

[0020] is the signal received by the first hydrophone of the linear array;

[0021] is the signal received by the second hydrophone of the linear array;

[0022] is the signal received by the m-th hydrophone of the linear array, m = 1, 2, …, M;

[0023] is the signal received by the M-th hydrophone of the linear array;

[0024] is the signal received by the M-element hydrophones of the linear array; T represents taking the transpose; t represents time;

[0025] Step 2: Assume that the signal received by the first hydrophone of the linear array is:

[0026]

[0027] where, ω is the angular frequency of the received signal, s(t) is the complex envelope of the received signal, j is the imaginary unit, j 2 = -1; the signal received by the m-th hydrophone is

[0028]

[0029] where, τ m is the time difference between the signals received by the m-th hydrophone and the first hydrophone from the target;

[0030] Approximate the complex envelope in Equation (3) as:

[0031] s(t - τ m ) ≈ s(t) (4)

[0032] Combining Equation (3) and Equation (4), the signal received by the m-th hydrophone is approximately

[0033]

[0034] Step 3: Uniformly represent the signals received by each hydrophone in the linear array as

[0035]

[0036] where, is the direction vector, τ 2 is the time difference between the second hydrophone and the first hydrophone receiving the target signal, τ M is the time difference between the Mth hydrophone and the first hydrophone receiving the target signal, and θ is the incident angle of the target;

[0037] The linear array receiving single target signal model is expressed as in Equation (6):

[0038] y(n) = a(θ)s(n) + v(n) (7)

[0039] where v(n) is the additive noise; s(n) is the received discrete signal; y(n) is the discrete output signal;

[0040] Establish the relationship between the time difference τ m between the mth hydrophone and the first hydrophone receiving the target signal and the incident angle θ of the target, as

[0041]

[0042] where c is the sound speed in water and d is the hydrophone interval;

[0043] Since ω = 2πf = 2πc / λ, and then substituting the direction vector a(θ) into Equation (8), the direction vector in the linear array receiving single target signal model is finally obtained as shown in Equation (9):

[0044] a(θ) = [1, e -jφ , …, e -j(m-1)φ , …, e -j(M-1)φ T , φ = 2πdsinθ / λ (9)

[0045] where f is the signal frequency, λ is the wavelength, and φ is the phase;

[0046] Step 14. Assume that the incident angles of N targets relative to the linear array are [θ 1 , θ 2 , …, θ N , and rewrite the direction vector as

[0047]

[0048] where A(θ) is the direction matrix, a(θ 1 ) is the direction vector of the first target, a(θ 2 ) is the direction vector of the second target, a(θ N ) is the direction vector of the Nth target;

[0049] ω 1 is the angular frequency of receiving the first target, ω​2 To receive the angular frequency of the second target, ω N is the angular frequency for receiving the Nth target;

[0050] Combined with Equation (7), the multi-target signal model received by the linear array is obtained, as shown in Equation (11):

[0051] y(n) = A(θ)s(n) + v(n) (11)

[0052] Step 15: Extract the azimuth measurement value of the platform and the spatial spectral intensity of the measurement value from each output signal y(n); the specific process is as follows:

[0053] For S' platforms equipped with passive sonars, assume the platform numbers are s = 1, 2,..., S';

[0054] The azimuth measurement value number of platform 1 is i 1 = 0, 1,... n 1 ;

[0055] The azimuth measurement value number of platform 2 is i 2 = 0, 1,... n 2 ;

[0056] The azimuth measurement value number of platform 3 is i 3 = 0, 1,... n 3 ;

[0057] The azimuth measurement value number of platform s is i s = 0, 1,..., n s ;

[0058] The azimuth measurement value number of platform S' is i S′ = 0, 1,..., n S′ ;

[0059] n s is the number of measurement values received by platform s, and 0 indicates a missed detection, that is, the situation where platform s does not detect a target;

[0060] represents the ith s azimuth measurement value of platform s;

[0061] At time k, one azimuth measurement value is taken from each platform to form an azimuth measurement set

[0062] At time k, the spatial spectral intensities of one measurement value are taken from each platform to form a spatial spectral intensity set

[0063] Among them,

[0064] Denote the $i$-th azimuth measurement value of platform 1 1 ; Denote the $i$-th azimuth measurement value of platform 2 2 ; Denote the $i$-th azimuth measurement value of platform $s$ s ; Denote the $i$-th azimuth measurement value of platform $S'$ S′ ;

[0065] Denote the spatial spectrum intensity of the $i$-th azimuth measurement value of platform 1 1 ; Denote the spatial spectrum intensity of the $i$-th azimuth measurement value of platform 2 2 ; Denote the spatial spectrum intensity of the $i$-th azimuth measurement value of platform $s$ s ; Denote the spatial spectrum intensity of the $i$-th azimuth measurement value of platform $S'$ S′ ;

[0066] Preferably, in the second step, based on each azimuth measurement set calculate the target position estimation value corresponding to each measurement set; based on each azimuth measurement set and the target position estimation value, calculate the cost function value of each measurement set Based on the cost function value construct an objective function, solve the objective function to obtain the allocation result; use the triangulation measurement method to estimate the target position for the allocation result; use the joint probabilistic data association method to track the target position;

[0067] Record the azimuth measurement values of each platform used for each tracking trajectory and the spatial spectrum intensity corresponding to the azimuth measurement values;

[0068] The specific process is as follows:

[0069] Step 2-1: Based on each azimuth measurement set calculate the target position estimation value corresponding to each measurement set; the specific process is as follows:

[0070] Since the target position $p$ is unknown, use the triangulation measurement method to estimate the target position;

[0071]

[0072] wherein, denotes the target position estimation value; $p$ denotes the target position;

[0073] denotes the azimuth measurement set the probability density function of the target $p$ from a known position;

[0074] Step 22: Based on each azimuth measurement set and the target position estimate calculate the cost function value of each measurement set

[0075] Step 23: Based on the cost function values in Step 22 construct an objective function and solve the objective function to obtain the allocation result;

[0076] Step 24: Use the triangulation measurement method to estimate the target position for the allocation result;

[0077] Use the joint probabilistic data association method to track the target position;

[0078] Assume that at the (k - 1)-th moment in a system with S' platforms, N targets are tracked. At the (k - 1)-th moment, the azimuth angle of the n-th target tracked relative to the s-th platform is

[0079] Record the mean of the spatial spectral intensity of the measurement values used by platform s from the start to the (k - 1)-th moment when tracking the n-th target, denoted as where represents the time window length.

[0080] Preferably, in Step 22, based on each azimuth measurement set and the target position estimate calculate the cost function value of each measurement set The specific process is as follows:

[0081] 1). The probability density function of the azimuth measurement set coming from the target position estimate is:

[0082]

[0083] where

[0084] is the probability density function of the measurement set coming from the target position estimate ;

[0085] is the detection probability of platform s;

[0086] u(i s ) is the indicator function; expressed as:

[0087]

[0088] is the measured value from the target probability density function; expressed as:

[0089]

[0090] wherein,

[0091] represents the i-th s measured value of platform s;

[0092] σ s represents the noise of platform s;

[0093] represents the measurement set estimated target position;

[0094] 2), assuming that the clutter is uniformly distributed in space, then the measurement set probability density function from the clutter in space is:

[0095]

[0096] wherein,

[0097] represents the measurement set probability density function from the clutter in space;

[0098] Φ represents the clutter;

[0099] represents the i-th s measured value probability from the clutter;

[0100] V s represents the observation area of platform s, that is is the clutter density;

[0101] 3), calculate the cost function value of each measurement set given by the negative log-likelihood ratio:

[0102]

[0103] Substitute Equation (13) and Equation (16) into Equation (17) to obtain the cost function value expressed as:

[0104]

[0105] Preferably, in step 23, based on the cost function value of step 22 Construct the objective function and solve the objective function to obtain the allocation result. The specific process is as follows:

[0106] The objective function is:

[0107]

[0108] Constraints:

[0109]

[0110] In the formula, is a binary variable;

[0111] When S' measurement sets are not associated with a certain target position estimate

[0112] When S' measurement sets are associated with a certain target position estimate

[0113] The corresponding of the minimum value of the objective function that satisfies the constraint conditions is the allocation result.

[0114] Preferably, in the third step, determine the size of the azimuth angle interval between the tracking trajectories of any two targets relative to the platform s and the azimuth angle interval threshold. If the azimuth angle interval between the tracking trajectories of the two targets relative to the platform s is less than the azimuth angle interval threshold, then open the observation window. After the platform s opens the observation window, determine whether the azimuth angle measurement value at time k is within ;

[0115] If the azimuth angle measurement value at time k is within , then enter the fourth step;

[0116] If the azimuth angle measurement value at time k is not within , then let k = k + 1 and execute the first step;

[0117] where ε θ is the azimuth angle expansion value of the observation window; represents the azimuth angle of the tracking trajectory of the i-th target relative to the s-th platform at time k - 1;

[0118] If the azimuth angle interval between the tracking trajectories of the two targets relative to the platform s is greater than or equal to the azimuth angle interval threshold, then do not open the observation window, and let k = k + 1 and execute the first step;

[0119] The specific process is as follows:

[0120] ​​Calculate the azimuth interval of the tracking trajectories of any two targets relative to platform s as shown in the following formula:

[0121]

[0122] where,

[0123] represents the azimuth interval of the tracking trajectories of the i-th target and the j-th target relative to platform s;

[0124] represents the azimuth of the tracking trajectory of the i-th target relative to the s-th platform at time k-1;

[0125] represents the azimuth of the tracking trajectory of the j-th target relative to the s-th platform at time k-1;

[0126] If the azimuth interval is greater than or equal to the azimuth interval threshold, do not open the observation window and execute Step 1;

[0127] If the azimuth interval is less than the azimuth interval threshold, determine that the two targets corresponding to the azimuth interval enter the proximity warning, and platform s opens the observation window:

[0128]

[0129] where, Γ θ is the azimuth interval threshold;

[0130] After platform s opens the observation window, determine whether the azimuth measurement value at time k is within ;

[0131] If the azimuth measurement value at time k is within , enter Step 4;

[0132] If the azimuth measurement value at time k is not within , let k = k + 1 and execute Step 1;

[0133] where ε θ is the azimuth expansion value of the observation window.

[0134] Preferably, in Step 4, calculate the change rate of the spatial spectrum intensity according to the spatial spectrum intensity corresponding to the azimuth measurement value obtained in Step 2 Compare the change rate of the spatial spectrum intensity with the intensity change threshold Γ A . If the change rate of the spatial spectrum intensity is greater than the intensity change threshold Γ A , then the corresponding measurement value is used as a "suspicious merged measurement value"; if the change rate of the spatial spectrum intensity Less than or equal to the intensity change threshold Γ A , then let k = k + 1 and execute Step 1;

[0135] The specific process is as follows:

[0136] Calculate the change rate of the spatial spectrum intensity according to the spatial spectrum intensity corresponding to the azimuth measurement value obtained in Step 2 The calculation formula is as follows,

[0137]

[0138] where,

[0139] is the change rate of the spatial spectrum intensity between the i s -th measurement value of the platform s at time k and the tracking trajectory of the n-th target;

[0140] is the spatial spectrum intensity of the i s -th measurement value of the platform s at time k;

[0141] is the average value of the spatial spectrum intensity of the tracking trajectory of the n-th target relative to the s-th platform in the previous time instants at time k;

[0142] Compare the change rate of the spatial spectrum intensity with the intensity change threshold Γ A . If the change rate of the spatial spectrum intensity is greater than the intensity change threshold Γ A , then the corresponding measurement value is used as a "suspicious merged measurement value"; if the change rate of the spatial spectrum intensity is less than or equal to the intensity change threshold Γ A , then let k = k + 1 and execute Step 1;

[0143] It is expressed as:

[0144]

[0145] Preferably, in Step 5, the "suspicious merged measurement values" are split, the split measurement values are recombined, the total association cost of each recombined measurement value is calculated, and the combination with the minimum total association cost is selected as the final allocation result. Based on the allocation result, the target position is estimated and the target position is tracked;

[0146] The specific process is as follows:

[0147] Step 5-1. Assume that there are J s "suspicious merged measurement values" in the platform s, numbered

[0148] Among them, j s,1 represents the number of the first "suspicious merged measurement value" of platform s, and j s,2 represents the number of the second "suspicious merged measurement value" of platform s, represents the number of the J s th "suspicious merged measurement value" of platform s;

[0149] Step Five Two: Perform "fission" on the J s "suspicious merged measurement values"; specifically as follows:

[0150]

[0151] Among them,

[0152] means adding the measurement value to the set Z s , and Z s represents the original measurement set of platform s;

[0153] is a binary decision variable;

[0154] means adding the j s th measurement value in platform s to the original measurement set of platform s, that is, performing one fission;

[0155] means not performing fission;

[0156]

[0157] Step Five Three: Let the lowest total cost obtained from Equation (19) be Calculate the lowest total cost of each combination, and select the measurement set corresponding to the minimum value among all total costs as the final allocation result;

[0158] Step Five Four: Estimate the target position based on the allocation result and track the target position.

[0159] Preferably, in Step Five Three, let the lowest total cost obtained from Equation (19) be Calculate the lowest total cost of each combination, and select the measurement set corresponding to the minimum value among all total costs as the final allocation result; as shown in the following formula:

[0160]

[0161] Among them, min represents selecting the minimum value in the group of values, and c min represents the minimum value in the total cost.

[0162] Preferably, in step 54, the target position is estimated based on the allocation result, and the target position is tracked; the specific process is as follows:

[0163] The triangulation measurement method is used to estimate the target position from the allocation result;

[0164] The joint probabilistic data association method is used to track the target position.

[0165] The beneficial effects of the present invention are as follows:

[0166] The present invention proposes an underwater multi-platform multi-target tracking method based on extended multi-dimensional assignment. This method combines the normalized intensity of the spatial spectrum and establishes an evaluation criterion with the lowest association cost, which can effectively split and merge measurement values, and improve the tracking and positioning accuracy of multi-targets according to the allocation result. Description of the Drawings

[0167] Figure 1 It is a flow chart of the present invention;

[0168] Figure 2 It is a motion situation diagram of the platform position and the target. The coordinate system is the northeast celestial coordinate system, the X-axis points north, the Y-axis points east, and the origin position is customized;

[0169] Figure 3 It is a positioning result diagram based on traditional multi-dimensional assignment;

[0170] Figure 4 It is a tracking result diagram based on traditional multi-dimensional assignment;

[0171] Figure 5 It is a positioning result diagram based on extended multi-dimensional assignment;

[0172] Figure 6 It is a tracking result diagram based on extended multi-dimensional assignment. Detailed Embodiment

[0173] Detailed Embodiment 1: The specific process of an underwater multi-platform multi-target tracking method based on extended multi-dimensional assignment in this embodiment is as follows:

[0174] Step 1: Let the time k = 1;

[0175] Each platform performs conventional beamforming processing to extract the azimuth measurement value of the platform and the corresponding spatial spectrum intensity;

[0176] Any one platform extracts one azimuth measurement value or does not extract it. If not extracted, the value is 0. The set of azimuth measurement values extracted by all platforms is represents the i-th 1 azimuth measurement value of platform 1; represents the i-th of platform 22 An azimuth measurement value; Indicates the i-th s azimuth measurement value of platform s; Indicates the i-th S′ azimuth measurement value of platform S';

[0177] For each platform, either extract 1 azimuth measurement value or not, and for all platforms, take all possible combinations to obtain the set of azimuth measurement values corresponding to all combination cases. One combination case corresponds to one set of azimuth measurement values;

[0178] Step 2: Based on each set of azimuth measurements Calculate the estimated target position corresponding to each measurement set; based on each set of azimuth measurements and the estimated target position, calculate the cost function value of each measurement set Based on the cost function value Construct an objective function, solve the objective function to obtain the allocation result; use the triangulation measurement method to estimate the target position for the allocation result (complete positioning); use the joint probabilistic data association method to track the target position;

[0179] Record the azimuth measurement values of each platform used for each tracking trajectory and the spatial spectral intensity corresponding to the azimuth measurement values;

[0180] Step 3: Judge the size of the azimuth interval between the tracking trajectories of any two targets relative to platform s and the azimuth interval threshold. If the azimuth interval between the tracking trajectories of two targets relative to platform s is less than the azimuth interval threshold, then open the observation window. After platform s opens the observation window, judge whether the azimuth measurement value at time k is within ;

[0181] If the azimuth measurement value at time k is within , then enter Step 4;

[0182] If the azimuth measurement value at time k is not within , then let k = k + 1 and execute Step 1;

[0183] where ε θ is the azimuth expansion value of the observation window; represents the azimuth of the tracking trajectory of the i-th target relative to the s-th platform at time k - 1;

[0184] If the azimuth interval between the tracking trajectories of two targets relative to platform s is greater than or equal to the azimuth interval threshold, then do not open the observation window, and let k = k + 1 and execute Step 1;

[0185] Step 4: Calculate the change rate of the spatial spectrum intensity based on the spatial spectrum intensity corresponding to the azimuth measurement value obtained in Step 2. Compare the change rate of the spatial spectrum intensity with the intensity change threshold Γ A If the change rate of the spatial spectrum intensity is greater than the intensity change threshold Γ A , then the corresponding measurement value is used as a "suspicious merging measurement value"; if the change rate of the spatial spectrum intensity is less than or equal to the intensity change threshold Γ A , then let k = k + 1 and execute Step 1.

[0186] Step 5: Split the "suspicious merging measurement value", recombine the split measurement values, calculate the total association cost of each recombined measurement value, select the combination with the minimum total association cost as the final allocation result, estimate the target position based on the allocation result, and track the target position.

[0187] Specific Embodiment 2: The difference between this embodiment and Specific Embodiment 1 is that in Step 1, each platform performs conventional beamforming processing (Equation 1-11) to extract the azimuth measurement value and the corresponding spatial spectrum intensity of the platform; the specific process is as follows:

[0188] Step 11: Assume that each platform is equipped with a uniform linear array of M hydrophones. In the scenario where the uniform linear array receives a far-field target signal, represent the signal vector received by the linear array as

[0189]

[0190] where,

[0191] is the signal received by the first hydrophone of the linear array;

[0192] is the signal received by the second hydrophone of the linear array;

[0193] is the signal received by the m-th hydrophone of the linear array, m = 1, 2,..., M;

[0194] is the signal received by the M-th hydrophone of the linear array; T represents taking the transpose; t represents time;

[0195] is the signal received by the M-element hydrophone of the linear array; T represents taking the transpose; t represents time;

[0196] Step 12: Assume that the signal received by the first hydrophone of the linear array is:

[0197]

[0198] where ω is the angular frequency of the received signal, s(t) is the complex envelope of the received signal, and j is the imaginary unit, where j 2 = -1;

[0199] Since the target satisfies the far-field condition, the deviation of the incident angle of the target signal arriving at each hydrophone of the linear array can be ignored. That is, the angles of the signal arriving at each hydrophone are the same, only the arrival times at each hydrophone are different. Therefore, the signal received by the m-th hydrophone is

[0200]

[0201] where τ m is the time difference between the m-th hydrophone and the first hydrophone receiving the far-field target signal;

[0202] Assuming that the target signal does not change rapidly, the complex envelope in Equation (3) is approximated as:

[0203] s(t - τ m ) ≈ s(t) (4)

[0204] In practical applications, it is generally independent of the received signal itself, that is, the received signal is independent of e jωt independent;

[0205] Combining Equation (3) and Equation (4), the signal received by the m-th hydrophone is approximately

[0206]

[0207] In Equation (5), it is deduced that the difference in the signals received by each hydrophone in the linear array is only the phase difference caused by the time delay;

[0208] Step One, represent the signals received by each hydrophone in the linear array uniformly as

[0209]

[0210] where is the direction vector, τ 2 is the time difference between the second hydrophone and the first hydrophone receiving the far-field target signal, τ M is the time difference between the M-th hydrophone and the first hydrophone receiving the far-field target signal, and θ is the incident angle of the far-field target;

[0211] In practical engineering, the signals collected by the system are discrete signals, and there will be noise interference in actual situations. Therefore, the single-target signal model received by the linear array, as shown in Equation (6), is expressed as:

[0212] y(n) = a(θ)s(n) + v(n) (7)

[0213] where \(v(n)\) is the additive noise; \(s(n)\) is the received discrete signal; \(y(n)\) is the discrete output signal;

[0214] Establish the time difference \(\tau\) between the received target signals of the \(m\)th hydrophone and the 1st hydrophone m and the relationship with the target incident angle \(\theta\), as

[0215]

[0216] where \(c\) is the sound speed in water and \(d\) is the hydrophone interval;

[0217] Equation (8) shows that the time difference of the signals received by the uniform linear array is related to the hydrophone interval \(d\), the sound speed \(c\) in water and the incident angle \(\theta\). However, the hydrophone interval and the sound speed are constant, and the time difference in the array received signal model is converted into the relationship with the incident angle.

[0218] Since \(\omega = 2\pi f = 2\pi c / \lambda\), and then combined with Equation (8) and substituting into the direction vector \(a(\theta)\), the direction vector in the linear array received single target signal model is finally obtained, as shown in Equation (9):

[0219] \(a(\theta)=[1,e\) -jφ ,…,e\) -j(m-1)φ ,…,e\) -j(M-1)φ T ,\(\varphi = 2\pi d\sin\theta / \lambda\) (9)

[0220] where \(f\) is the signal frequency, \(\lambda\) is the wavelength, and \(\varphi\) is the phase;

[0221] Step 14. In the actual scenario, there may be multiple targets. Assume that the incident angles of \(N\) far-field targets relative to the linear array are \([\theta\) 1 ,\(\theta\) 2 ,…,\(\theta\) N , and rewrite the direction vector as

[0222]

[0223] where \(A(\theta)\) is the direction matrix, \(a(\theta\) 1 ) is the direction vector of the 1st far-field target, \(a(\theta\) 2 ) is the direction vector of the 2nd far-field target, \(a\theta\) N ) is the direction vector of the \(N\)th far-field target;

[0224] \(\omega\) 1 is the angular frequency of receiving the 1st target, \(\omega\) 2 is the angular frequency of receiving the 2nd target, \(\omega\) N is the angular frequency of receiving the \(N\)th target;

[0225] ​Combined with Equation (7), the linear array receiving multi-target signal model is obtained, as shown in Equation (11):

[0226] y(n) = A(θ)s(n) + v(n) (11)

[0227] Step 15: Extract the azimuth measurement value of the platform and the spatial spectrum intensity of the measurement value from each output signal y(n); the specific process is as follows:

[0228] Each output signal y(n) corresponds to 1 platform, and only 1 azimuth measurement value of the platform and the spatial spectrum intensity of the measurement value can be extracted from each output signal y(n);

[0229] For S′ platforms equipped with passive sonars, assume the platform numbers are s = 1, 2, …, S′;

[0230] The azimuth measurement value number of platform 1 is i 1 = 0, 1, … n 1 ;

[0231] The azimuth measurement value number of platform 2 is i 2 = 0, 1, … n 2 ;

[0232] The azimuth measurement value number of platform 3 is i 3 = 0, 1, … n 3 ;

[0233] The azimuth measurement value number of platform s is i s = 0, 1, …, n s ;

[0234] The azimuth measurement value number of platform S′ is i S′ = 0, 1, …, n S′ ;

[0235] n s is the number of measurement values received by platform s, 0 indicates a missed detection, that is, the situation where platform s does not detect a target;

[0236] represents the i s th azimuth measurement value of platform s;

[0237] At time k, 1 azimuth measurement value is taken from each platform to form an azimuth measurement set (1 value is taken from 1 platform, it can be not taken, if not taken, the value is 0);

[0238] At time k, 1 spatial spectrum intensity of the measurement value is taken from each platform to form a spatial spectrum intensity set (1 value is taken from 1 platform, it can be not taken, if not taken, the value is 0);

[0239]

[0240] Among them,

[0241] represents the i-th 1 azimuth measurement value of platform 1; represents the i-th 2 azimuth measurement value of platform 2; represents the i-th s azimuth measurement value of platform s; represents the i-th S′ azimuth measurement value of platform S';

[0242] represents the spatial spectrum intensity of the i-th 1 azimuth measurement value of platform 1; represents the spatial spectrum intensity of the i-th 2 azimuth measurement value of platform 2; represents the spatial spectrum intensity of the i-th s azimuth measurement value of platform s; represents the spatial spectrum intensity of the i-th S′ azimuth measurement value of platform S'.

[0243] Other steps and parameters are the same as those in the first specific implementation manner.

[0244] Third specific implementation manner: The difference between this implementation manner and the first or second specific implementation manner is that in step 2, based on each azimuth measurement set calculate the target position estimation value corresponding to each measurement set; based on each azimuth measurement set and the target position estimation value, calculate the cost function value of each measurement set Based on the cost function value construct an objective function, solve the objective function to obtain the allocation result; use the triangulation measurement method to estimate the target position for the allocation result (complete positioning); use the joint probabilistic data association method to track the target position;

[0245] Record the azimuth measurement values of each platform used for each tracking trajectory and the spatial spectrum intensity corresponding to the azimuth measurement values;

[0246] The S-D assignment algorithm enumerates all possible combination methods based on the criterion of minimizing the established cost function, calculates the association cost of each combination, and thus selects the combination with the lowest cost as the optimal association combination;

[0247] The specific process is as follows:

[0248] Step 2-1: Based on each azimuth measurement set Calculate the estimated target position value corresponding to each measurement set; the specific process is as follows:

[0249] Since the target position p is unknown, the least squares estimation is used instead, and the method of triangular measurement is used to estimate the target position;

[0250]

[0251] Among them, represents the estimated target position value; p represents the target position;

[0252] represents the azimuth measurement set The probability density function of the target p from a known position;

[0253] Step two: Based on each azimuth measurement set and the estimated target position value Calculate the cost function value of each measurement set

[0254] Step two-three: Based on the cost function value in step two Construct an objective function and solve the objective function to obtain the allocation result;

[0255] Step two-four: Use the triangular measurement method to estimate the target position for the allocation result (complete positioning);

[0256] Use the joint probability data association method to track the target position;

[0257] Assume that at the (k - 1)-th moment in a system with S' platforms, N targets are tracked. At the (k - 1)-th moment, the azimuth angle of the n-th target tracked relative to the s-th platform is

[0258] Record the mean value of the spatial spectrum intensity of the measurement values used by platform s to track the n-th target from to the (k - 1)-th moment, denoted as represents the time window length.

[0259] Other steps and parameters are the same as those in the first or second specific implementation manner.

[0260] Specific implementation manner four: The difference between this implementation manner and one of the first to third specific implementation manners is that in the above step two, based on each azimuth measurement set and the estimated target position value Calculate the cost function value of each measurement set The specific process is as follows:

[0261] 1). The azimuth measurement set from the estimated target position The probability density function is as follows:

[0262]

[0263] where

[0264] is the measurement set from the estimated target position The probability density function;

[0265] is the detection probability of platform s;

[0266] u(i s ) is the indicator function; expressed as:

[0267]

[0268] is the measurement value from the target The probability density function; expressed as:

[0269]

[0270] where

[0271] represents the i-th s measurement value of platform s;

[0272] σ s represents the noise of platform s;

[0273] represents the measurement set The estimated target position;

[0274] 2), Assume that the clutter is uniformly distributed in space, then the probability density function of the measurement set from the clutter in space is:

[0275]

[0276] where

[0277] represents the probability density function of the measurement set from the clutter in space;

[0278] Φ represents the clutter;

[0279] represents the i-th s measurement value from the clutter;

[0280] V s represents the observed area of platform s, that is is the clutter density;

[0281] 3), calculate the cost function value of each measurement set given by the negative log-likelihood ratio:

[0282]

[0283] Substitute equations (13) and (16) into equation (17) to obtain the cost function value expressed as:

[0284]

[0285] Other steps and parameters are the same as those in one of the specific embodiments 1 to 3.

[0286] Specific embodiment 5: The difference between this embodiment and one of the specific embodiments 1 to 4 is that: in step 23, based on the cost function value of step 22 construct an objective function, solve the objective function to obtain the allocation result; the specific process is:

[0287] The objective function is:

[0288]

[0289] Constraints:

[0290]

[0291] In the formula, is a binary variable;

[0292] When S' measurement sets are not associated with a certain target position estimate (given by equation 12),

[0293] When S' measurement sets are associated with a certain target position estimate (given by equation 12),

[0294] The corresponding of the minimum value of the objective function that satisfies the constraint conditions is the allocation result.

[0295] Other steps and parameters are the same as those in one of the specific embodiments 1 to 4.

[0296] Embodiment Six: The difference between this embodiment and any one of Embodiments One to Five is as follows: In Step 3, it is judged whether the azimuth interval between the tracking trajectories of any two targets relative to platform s is smaller than the azimuth interval threshold. If the azimuth interval between the tracking trajectories of the two targets relative to platform s is smaller than the azimuth interval threshold, the observation window is opened. After platform s opens the observation window, it is judged whether the azimuth measurement value at time k is within ;

[0297] If the azimuth measurement value at time k is within , then go to Step 4;

[0298] If the azimuth measurement value at time k is not within , then let k = k + 1 and execute Step 1;

[0299] where ε θ is the azimuth expansion value of the observation window; represents the azimuth of the tracking trajectory of the i-th target relative to the s-th platform at time k - 1;

[0300] If the azimuth interval between the tracking trajectories of the two targets relative to platform s is greater than or equal to the azimuth interval threshold, the observation window is not opened, then let k = k + 1 and execute Step 1;

[0301] The specific process is as follows:

[0302] Calculate the azimuth interval between the tracking trajectories of any two targets relative to platform s, as shown in the following formula:

[0303]

[0304] where,

[0305] represents the azimuth interval between the tracking trajectories of the i-th target and the j-th target relative to platform s;

[0306] represents the azimuth of the tracking trajectory of the i-th target relative to the s-th platform at time k - 1;

[0307] represents the azimuth of the tracking trajectory of the j-th target relative to the s-th platform at time k - 1;

[0308] When the azimuth interval becomes smaller and smaller, it indicates that the targets may be approaching in azimuth. When it is smaller than the system azimuth resolution, they will merge in the spatial spectrum.

[0309] Take the azimuth interval as the test statistic;

[0310] If the azimuth interval is greater than or equal to the azimuth interval threshold, the observation window is not opened, and let k = k + 1 to execute Step 1;

[0311] If the azimuth interval is less than the azimuth interval threshold, it is determined that the two targets corresponding to the azimuth interval enter the proximity warning, and the platform s opens the observation window:

[0312]

[0313] Among them, Γ θ is the azimuth interval threshold;

[0314] After the platform s opens the observation window, it is judged whether the azimuth measurement value at time k is within ;

[0315] If the azimuth measurement value at time k is within , go to Step 4;

[0316] If the azimuth measurement value at time k is not within , let k = k + 1 to execute Step 1;

[0317] Among them, ε θ is the azimuth expansion value of the observation window.

[0318] Other steps and parameters are the same as those in any one of the specific embodiments 1 to 5.

[0319] Specific embodiment 7: The difference between this embodiment and any one of the specific embodiments 1 to 6 is that in Step 4, according to the spatial spectrum intensity corresponding to the azimuth measurement value obtained in Step 2, the change rate of the spatial spectrum intensity is calculated The change rate of the spatial spectrum intensity is compared with the intensity change threshold Γ A . If the change rate of the spatial spectrum intensity is greater than the intensity change threshold Γ A , the corresponding measurement value is used as the "suspicious combined measurement value"; if the change rate of the spatial spectrum intensity is less than or equal to the intensity change threshold Γ A , let k = k + 1 to execute Step 1;

[0320] The specific process is as follows:

[0321] According to the spatial spectrum intensity corresponding to the azimuth measurement value obtained in Step 2, the change rate of the spatial spectrum intensity is calculated The calculation formula is as follows,

[0322]

[0323] Among them,

[0324] is the rate of change of the spatial spectral intensity between the i-th measurement value of platform s at time k and the tracking trajectory of the n-th target; s ;

[0325] is the spatial spectral intensity of the i-th measurement value of platform s at time k; s ;

[0326] is the average value of the spatial spectral intensity of the tracking trajectory of the n-th target relative to the s-th platform at the previous time instants at time k;

[0327] Compare the rate of change of the spatial spectral intensity with the intensity change threshold Γ A . If the rate of change of the spatial spectral intensity is greater than the intensity change threshold Γ A , then the corresponding measurement value is regarded as a "suspicious merged measurement value"; if the rate of change of the spatial spectral intensity is less than or equal to the intensity change threshold Γ A , then let k = k + 1 and execute Step 1;

[0328] It is expressed as:

[0329]

[0330] Other steps and parameters are the same as those in Embodiments 1 to 6.

[0331] Embodiment 8: The difference between this embodiment and Embodiment 1 to 6 is that in Step 5, the "suspicious merged measurement value" is split, the split measurement values are recombined, the total association cost of each recombined measurement value is calculated, and the combination with the minimum total association cost is selected as the final allocation result, and the target position is estimated based on the allocation result to track the target position;

[0332] The specific process is as follows:

[0333] Step 5-1. Assume that there are J s "suspicious merged measurement values" in platform s, numbered

[0334] where j s,1 represents the number of the first "suspicious merged measurement value" of platform s, and j s,2 represents the number of the second "suspicious merged measurement value" of platform s, represents the number of the J s -th "suspicious merged measurement value" of platform s;

[0335] Step 5-2. For J s"Fission" is performed on a "suspicious combined measurement value"; specifically as follows:

[0336]

[0337] Among them, means adding the measurement value to the set Z s , and Z s represents the original measurement set of platform s (platform s has its own original measurement set);

[0338] are binary decision variables;

[0339] means adding the j s th measurement value in platform s to the original measurement set of platform s, that is, performing one fission;

[0340] means not performing fission;

[0341]

[0342] Step Five Three. Let the lowest total cost obtained from Equation (19) be Calculate the lowest total cost of each combination, and select the measurement set corresponding to the minimum value among all the total costs as the final allocation result;

[0343] Step Five Four. Estimate the target position based on the allocation result and track the target position.

[0344] Other steps and parameters are the same as those in the first to seventh specific embodiments.

[0345] Specific Embodiment Nine: The difference between this embodiment and any one of the first to eighth specific embodiments is that in Step Five Three, the lowest total cost obtained from Equation (19) is Calculate the lowest total cost of each combination, and select the measurement set corresponding to the minimum value among all the total costs as the final allocation result; as shown in the following formula:

[0346]

[0347] Among them, min represents selecting the minimum value in the group of values, and c min represents the minimum value among the total costs.

[0348] Other steps and parameters are the same as those in the first to eighth specific embodiments.

[0349] Specific Embodiment Ten: The difference between this embodiment and any one of the first to ninth specific embodiments is that in Step Five Four, the target position is estimated based on the allocation result and the target position is tracked; the specific process is:

[0350] Use triangulation method to estimate the target position (complete positioning) based on the allocation results;

[0351] The target location is tracked using a joint probabilistic data association method.

[0352] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 9.

[0353] The following examples are used to verify the beneficial effects of the present invention:

[0354] Embodiment 1:

[0355] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0356] First, construct the platform position and target motion trajectory. This example contains 4 targets and 3 targets. Figure 2 The positions of the four platforms are (300m, 350m), (1100m, 50m), (1900m, 50m) and (2700m, 350m), and the initial positions of the three targets are (1100m, 1700m), (2000m, 1750m) and (2000m, 17500m). The positioning results based on traditional multi-dimensional allocation are shown in Figure 2. Figure 3 As shown in Figure 2, target 2 has many missing measurements. The tracking results are as follows: Figure 4 As shown in Figure 2, target 2 is tracked out of range. The positioning result based on extended multi-dimensional allocation is shown in Figure 2. Figure 5 The corresponding tracking results are shown in Figure 6 As shown in Figure 3, the results show that this method has good tracking performance.

[0357] The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. An underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation, characterized by: The specific process of the method is: Step 1: Set time k=1; Each platform performs conventional beamforming processing to extract the platform's azimuth measurement value and the corresponding spatial spectrum intensity; Step 2: Based on each azimuth measurement set Calculate the target position estimate corresponding to each measurement set; Based on each azimuth measurement set And the target position estimate, calculate the cost function value of each measurement set Based on the cost function value Construct an objective function and solve it to obtain the allocation result; The target location is estimated using triangulation method on the allocation results; The target location is tracked using a joint probabilistic data association method; Record the azimuth angle measurement value of each platform used in each tracking trajectory and the spatial spectrum intensity corresponding to the azimuth angle measurement value; Step 3: Determine the size of the azimuth interval of the tracking trajectories of any two targets relative to platform s and the azimuth interval threshold. If the azimuth interval of the tracking trajectories of the two targets relative to platform s is less than the azimuth interval threshold, open the observation window. After platform s opens the observation window, determine whether the azimuth measurement value at time k is within Inside; If the azimuth angle measurement value at time k is If yes, go to step 4; If the azimuth angle measurement value at time k is not If k is within, let k=k+1 and execute step 1; where ε θ is the observation window azimuth extension value; represents the azimuth of the tracking trajectory of the i-th target at time k-1 relative to the s-th platform; If the azimuth interval of the tracking trajectories of the two targets relative to the platform s is greater than or equal to the azimuth interval threshold, the observation window is not opened, and k=k+1 is set to execute step 1; Step 4: Calculate the rate of change of the spatial spectrum intensity according to the spatial spectrum intensity corresponding to the azimuth measurement value obtained in step 2 The rate of change of spatial spectrum intensity With the intensity change threshold Γ A By comparison, if the rate of change of spatial spectrum intensity Greater than the intensity change threshold Γ A , then the corresponding measurement value is regarded as the "suspicious combined measurement value"; if the rate of change of the spatial spectrum intensity Less than or equal to the intensity change threshold Γ A , then let k=k+1 and execute step 1; Step 5: Split the "suspicious merged measurement values", recombine the split measurement values, calculate the total association cost of each recombination measurement value, select the combination with the smallest total association cost as the final allocation result, estimate the target position based on the allocation result, and track the target position.

2. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 1, characterized in that: In step 1, each platform performs conventional beamforming processing to extract the azimuth angle measurement value of the platform and the corresponding spatial spectrum intensity; the specific process is: Step 1: Assume that each platform is equipped with a uniform linear array of M hydrophones. The signal vector received by the linear array is expressed as in, is the signal received by the first hydrophone of the linear array; is the signal received by the second hydrophone of the linear array; is the signal received by the mth hydrophone of the linear array, m = 1, 2, …, M; is the signal received by the Mth hydrophone of the linear array; is the signal received by the M-element linear array hydrophone; T represents the transpose; t represents time; Step 1 and 2: Assume that the signal received by the first hydrophone of the linear array is: Where ω is the angular frequency of the received signal, s(t) is the complex envelope of the received signal, j is the imaginary unit, and j 2 = -1; The signal received by the mth hydrophone is Among them, τ m is the time difference between the mth hydrophone and the first hydrophone receiving the target signal; The complex envelope in equation (3) is approximated as: s(t-τ m )≈s(t) (4) Combining equations (3) and (4), the signal received by the mth hydrophone is approximately Step 13: The signals received by each hydrophone in the linear array are uniformly expressed as in, is the direction vector, τ2 is the time difference between the second hydrophone and the first hydrophone receiving the target signal, τ M is the time difference between the Mth hydrophone and the first hydrophone receiving the target signal, and θ is the incident angle of the target; The linear array receiving single target signal model is expressed as follows, as shown in equation (6): y(n)=a(θ)s(n)+v(n) (7) Where, v(n) is additive noise; s(n) is the received discrete signal; y(n) is the discrete output signal; Establish the time difference τ between the mth hydrophone and the first hydrophone receiving the target signal m The relationship with the target incident angle θ is as follows: Where c is the speed of sound in water, d is the hydrophone spacing; Since ω=2πf=2πc / λ, and then combined with equation (8) to bring in the direction vector a(θ), we can finally get the direction vector in the linear array receiving single target signal model, as shown in equation (9): a(θ)=[1,e -jφ ,…,e -j(m-1)φ ,…,e -j(M-1)φ ] T ,φ=2πd sinθ / λ (9) Where f is the signal frequency, λ is the wavelength, and φ is the phase; Step 14: Assume that the incident angles of N targets relative to the linear array are [θ1,θ2,…,θ N ], rewrite the direction vector as Among them, A(θ) is the direction matrix, a(θ1) is the direction vector of the first target, a(θ2) is the direction vector of the second target, and a(θ N ) is the direction vector of the Nth target; ω1 is the angular frequency of receiving the first target, ω2 is the angular frequency of receiving the second target, ω N is the angular frequency of receiving the Nth target; Combining with equation (7), we get the linear array receiving multi-target signal model, as shown in equation (11): y(n)=A(θ)s(n)+v(n) (11) Step 15: Extract the platform's azimuth angle measurement value and the spatial spectrum intensity of the measurement value from each output signal y(n); the specific process is: For S′ platforms equipped with passive sonar, assume that the platforms are numbered s=1, 2, …, S′; The azimuth angle measurement values ​​of platform 1 are numbered as i1=0,1,…n1; The azimuth angle measurement values ​​of platform 2 are numbered as i2=0, 1, ... n2; The azimuth angle measurement values ​​of platform 3 are numbered as i3=0,1,…n3; The azimuth measurement value of platform s is numbered as i s =0,1,…,n s ; The azimuth measurement value of platform S′ is numbered as i S′ =0,1,…,n S′ ; n s is the number of measurement values ​​received by platform s, 0 indicates missed detection, that is, platform s does not detect the target; Indicates platform s i s Azimuth measurements; At time k, take one azimuth measurement value on each platform to form an azimuth measurement set The spatial spectrum intensity of one measurement value is taken at each platform at time k to form a spatial spectrum intensity set in, represents the i1th azimuth measurement value of platform 1; represents the i2th azimuth measurement value of platform 2; Indicates platform s i s Azimuth measurements; The i-th platform S′ S′ Azimuth measurements; represents the spatial spectrum intensity of the i1th azimuth measurement value of platform 1; represents the spatial spectrum intensity of the i2th azimuth measurement value of platform 2; Indicates platform s i s The spatial spectrum intensity of the azimuth measurement value; The i-th platform S′ S′ The spatial spectrum intensity of the azimuth measurement value.

3. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 2, characterized in that: In step 2, based on each azimuth measurement set Calculate the target position estimate corresponding to each measurement set; based on each azimuth measurement set And the target position estimate, calculate the cost function value of each measurement set Based on the cost function value Construct an objective function and solve it to obtain the allocation result; The target location is estimated using triangulation method on the allocation results; The target location is tracked using a joint probabilistic data association method; Record the azimuth angle measurement value of each platform used in each tracking trajectory and the spatial spectrum intensity corresponding to the azimuth angle measurement value; The specific process is: Step 21: Based on each azimuth measurement set Calculate the target position estimate corresponding to each measurement set; the specific process is: Since the target position p is unknown, the triangulation method is used to estimate the target position; in, represents the estimated value of the target position; p represents the target position; Represents a set of azimuth measurements The probability density function of the target p from a known position; Step 22: Based on each azimuth measurement set and target position estimate Calculate the cost function value for each measurement set Step 23: Cost function value based on step 22 Construct an objective function and solve it to obtain the allocation result; Step 24: Use triangulation method to estimate the target position based on the allocation results; The target location is tracked using a joint probabilistic data association method; Assume that N targets are tracked at the k-1th time in a system with S′ platforms. The azimuth of the nth target tracked at the k-1th time relative to the sth platform is Recording platforms The mean value of the spatial spectrum intensity of the measurement value used to track the nth target at time k-1 is recorded as Indicates that the time window is long.

4. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 3 is characterized by: In step 22, based on each azimuth measurement set and target position estimate Calculate the cost function value for each measurement set The specific process is: 1) Azimuth measurement set From the target position estimate The probability density function of is: in, For measurement set From the target position estimate The probability density function of is the detection probability of platform s; u(i s ) is the indicator function; it is expressed as: is the measured value From Target The probability density function of ; expressed as: in, represents the i-th platform s s A measurement value; σ s represents the noise of platform s; Represents a measurement set estimated target location; 2) Assuming that the clutter is uniformly distributed in space, the measurement set The probability density function of the clutter from space is: in, Represents a measurement set The probability density function of the clutter in the space; Φ represents clutter; represents the i-th platform s s Measurement value Probability from clutter; V s represents the observation area of ​​platform s, that is, is the clutter density; 3) Calculate the cost function value of each measurement set Given by the negative log-likelihood ratio: Substituting equations (13) and (16) into equation (17), we get the cost function value: It is expressed as:

5. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 4 is characterized in that: The cost function value in step 23 based on step 22 Construct the objective function and solve it to get the allocation result; the specific process is: The objective function is: Constraints: In the formula, is a binary variable; When S′ measurement sets With a target position estimate When there is no association, When S′ measurement sets With a target position estimate When associated, The minimum value of the objective function that satisfies the constraints corresponds to Assignment result.

6. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 5, characterized in that: In the step 3, the azimuth interval of the tracking trajectories of any two targets relative to the platform s and the azimuth interval threshold are determined. If the azimuth interval of the tracking trajectories of the two targets relative to the platform s is less than the azimuth interval threshold, the observation window is opened. After the observation window is opened on the platform s, it is determined whether the azimuth measurement value at time k is within Inside; If the azimuth angle measurement value at time k is If yes, go to step 4; If the azimuth angle measurement value at time k is not If k is within, let k=k+1 and execute step 1; where ε θ is the observation window azimuth extension value; represents the azimuth of the tracking trajectory of the i-th target at time k-1 relative to the s-th platform; If the azimuth interval of the tracking trajectories of the two targets relative to the platform s is greater than or equal to the azimuth interval threshold, the observation window is not opened, and k=k+1 is set to execute step 1; The specific process is: Calculate the azimuth angle interval of the tracking trajectories of any two targets relative to the platform s, as shown in the following formula: in, represents the azimuth interval between the tracking trajectories of the i-th target and the j-th target relative to the platform s; represents the azimuth of the tracking trajectory of the i-th target at time k-1 relative to the s-th platform; represents the azimuth angle of the tracking trajectory of the jth target at time k-1 relative to the sth platform; If the azimuth interval is greater than or equal to the azimuth interval threshold, the observation window is not opened, and k=k+1 is set to execute step 1; If the azimuth interval is less than the azimuth interval threshold, the two targets corresponding to the azimuth interval are judged to enter the proximity warning, and the platform s opens the observation window: Among them, Γ θ is the azimuth interval threshold; After the observation window is opened on platform s, determine whether the azimuth angle measurement value at time k is within Inside; If the azimuth angle measurement value at time k is If yes, go to step 4; If the azimuth angle measurement value at time k is not If k is within, let k=k+1 and execute step 1; where ε θ is the observation window azimuth extension value.

7. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 6, characterized in that: In step 4, the rate of change of the spatial spectrum intensity is calculated based on the spatial spectrum intensity corresponding to the azimuth angle measurement value obtained in step 2. The rate of change of spatial spectrum intensity With the intensity change threshold Γ A By comparison, if the rate of change of spatial spectrum intensity Greater than the intensity change threshold Γ A , then the corresponding measurement value is regarded as the "suspicious combined measurement value"; if the rate of change of the spatial spectrum intensity Less than or equal to the intensity change threshold Γ A , then let k=k+1 and execute step 1; The specific process is: According to the spatial spectrum intensity corresponding to the azimuth measurement value obtained in step 2, calculate the rate of change of the spatial spectrum intensity The calculation formula is as follows, in, is the i-th platform s at time k s The rate of change of the spatial spectrum intensity of the nth measurement value and the tracking trajectory of the nth target; is the i-th platform s at time k s The spatial spectrum intensity of the measured values; is the number before time k The average value of the spatial spectrum intensity of the tracking trajectory of the nth target relative to the sth platform at the moment; The rate of change of spatial spectrum intensity With the intensity change threshold Γ A By comparison, if the rate of change of spatial spectrum intensity Greater than the intensity change threshold Γ A , then the corresponding measurement value is regarded as the "suspicious combined measurement value"; if the rate of change of the spatial spectrum intensity Less than or equal to the intensity change threshold Γ A , then let k=k+1 and execute step 1; It is expressed as:

8. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 7, characterized in that: In the step 5, the "suspicious combined measurement value" is split, the split measurement values ​​are recombined, the total association cost of each recombined measurement value is calculated, the combination with the smallest total association cost is selected as the final allocation result, the target position is estimated based on the allocation result, and the target position is tracked; The specific process is: Step 51: Assume that there is J in platform s s "Suspicious combined measurements", numbered Among them, j s,1 Indicates the first "suspicious combined measurement value" number of platform s, j s,2 Indicates the second "suspicious combined measurement value" number of platform s. Indicates platform s J s A "Suspicious Combined Measurement Value" number; Step 52: J s The "suspicious combined measurement value" is "split"; the details are as follows: in, Indicates that the measured value Add to collection Z s , Z s represents the original measurement set of platform s; is a binary decision variable; Indicates that the jth s The measurement value is added to the original measurement set of platform s, that is, a fission is performed; It means no fission. Step 53: Assume that the minimum total cost obtained by equation (19) is Calculate the minimum total cost of each combination, and select the measurement set corresponding to the minimum value of all total costs as the final allocation result; Step 54: Estimate the target position based on the allocation result and track the target position.

9. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 8, characterized in that: In step 53, the minimum total cost obtained by setting equation (19) is Calculate the minimum total cost of each combination, and select the measurement set corresponding to the minimum value of all total costs as the final allocation result; As shown below: Among them, min means selecting the minimum value in the group value, c min Represents the minimum value in the total cost.

10. The underwater multi-platform multi-target tracking method based on extended multi-dimensional allocation according to claim 9, characterized in that: In step 54, the target position is estimated based on the allocation result and the target position is tracked; the specific process is: The target location is estimated using triangulation method on the allocation results; The target location is tracked using a joint probabilistic data association method.

Images