An underwater multi-platform multi-target tracking method based on extended multi-dimensional assignment

Abstract

The application relates to an underwater multi-platform multi-target tracking method based on extended multi-dimensional distribution, belongs to the underwater target tracking field, and relates to an underwater multi-platform multi-target tracking method. The application aims to solve the problem that in the prior art, when the bearing angle measurement values of platforms are subjected to multi-dimensional distribution and positioning and tracking, the bearing angle measurement values may be combined due to the limitation of the conventional beam forming angle resolution, thereby affecting the positioning and tracking accuracy. The application provides an underwater multi-platform multi-target tracking method based on extended multi-dimensional distribution. The method combines the normalized intensity of a space spectrum, establishes an evaluation criterion of minimum total correlation cost, can effectively split and combine the measurement values, and improves the positioning and tracking performance of multi-targets according to the distribution result.

Description

Technical Field

[0001] This invention belongs to the field of underwater target tracking and relates to a method for tracking multiple targets on multiple underwater platforms. Background Technology

[0002] In underwater multi-target tracking scenarios, integrating measurement information from multiple platforms can compensate for the limitations of a single platform in terms of maneuverability, observation range, and accuracy, making it a key research focus in recent years. In applied engineering, the following process is typically used to achieve azimuth estimation: first, threshold detection is performed on the azimuth-energy spectrum generated by conventional beamforming; 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, the azimuth angle interval between targets with the same frequency signal in space is smaller than the main lobe beamwidth of the array, leading to the problem of measurement value "merging." Traditional methods for associating the azimuth angles of various platforms use multi-dimensional allocation, with the objective function being the lowest total association cost. However, the missing measurement values ​​resulting from measurement value "merging" degrade positioning and tracking performance. Therefore, it is necessary to study how to reduce the impact of measurement value merging caused by beamforming angular resolution on positioning and tracking performance. Summary of the Invention

[0003] The purpose of this invention is to solve the problem that in existing methods, when the peak value of conventional beamforming is used as the azimuth measurement value and then multidimensionally allocated for positioning and tracking, the azimuth measurement value may "merge" due to the limitation of the angular resolution of conventional beamforming, thus affecting the accuracy of positioning and tracking. Therefore, this invention proposes an underwater multi-platform multi-target tracking method based on extended multidimensional allocation.

[0004] The specific process of an underwater multi-platform multi-target tracking method based on extended multidimensional allocation is as follows:

[0005] Step 1: Set time k = 1;

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

[0007] Step 2: Based on each azimuth measurement set Calculate the target position estimate for each measurement set; based on each azimuth measurement set... Based on the target location estimate, calculate the cost function value for each measurement set. Based on cost function value Construct an objective function and solve it to obtain the allocation result; use triangulation to estimate the target position based on the allocation result; use a joint probability data association method to track the target position.

[0008] Record the azimuth measurements of each platform used for each tracking trajectory and the corresponding spatial spectrum intensity;

[0009] Step 3: Determine 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 the two targets relative to platform s is less than the azimuth interval threshold, then open the observation window. After opening the observation window on platform s, determine whether the azimuth measurement value at time k is within the threshold. Inside;

[0010] If the azimuth measurement at time k is If inside, proceed to step four;

[0011] If the azimuth measurement at time k is not If k = k + 1, then execute step one;

[0012] Where ε θ This refers to the azimuth extension value of the observation window; This represents the azimuth angle of the tracking trajectory of the i-th target relative to the s-th platform at time k-1;

[0013] If the azimuth interval between 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 one.

[0014] Step 4: Calculate the rate of change of the spatial spectrum intensity based on the azimuth measurement obtained in Step 2. The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one;

[0015] Step 5: Split the “suspicious merged measurement values”, recombine the split measurement values, calculate the total association cost of each recombinated measurement value, select the combination with the smallest total association cost as the final allocation result, estimate the target location based on the allocation result, and track the target location.

[0016] Preferably, in step one, each platform undergoes conventional beamforming processing to extract the azimuth measurement value and the corresponding spatial spectrum intensity of the platform; the specific process is as follows:

[0017] Step 11: Assume each platform is equipped with a uniform linear array of M-element hydrophones. The signal vector received by the linear array is represented as...

[0018]

[0019] in,

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

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

[0022] Let m be the signal received by the m-th hydrophone in the linear array, where m = 1, 2, ..., M;

[0023] This is the signal received by the Mth hydrophone in the linear array;

[0024] The signal received by the linear array M-element hydrophone; T represents the transpose; t represents time;

[0025] Steps 1 and 2: Assume the signal received by the first hydrophone in 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, and j is the imaginary unit. 2 =-1; the signal received by the m-th hydrophone is

[0028]

[0029] Where, τ m Let m be the time difference between the m-th hydrophone and the 1st hydrophone receiving the target signal;

[0030] The complex envelope in equation (3) can be approximated as:

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

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

[0033]

[0034] Step 13: Represent the signals received by each hydrophone in the linear array as a unified representation.

[0035]

[0036] in, Let τ be the direction vector, and τ2 be the time difference between the second hydrophone and the first hydrophone receiving the target signal. M Let θ be the time difference between the Mth hydrophone and the first hydrophone receiving the target signal, and let θ be the incident angle of the target.

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

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

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

[0040] Establish the time difference τ between the m-th hydrophone and the 1st hydrophone in receiving the target signal. m The relationship with the target incident angle θ is as follows:

[0041]

[0042] Where c is the speed of sound in water, and d is the hydrophone spacing;

[0043] Since ω=2πf=2πc / λ, and combining equation (8) with the direction vector a(θ), 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 1, Section 4: Assume the incident angles of N targets relative to the linear array are [θ1, θ2, ..., θ]. N ], 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, and a(θ3) is the direction vector of the third target. N () represents the direction vector of the Nth target;

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

[0050] Combining equation (7), we obtain the linear array receiving multi-target signal model, as shown in equation (11):

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

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

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

[0054] The azimuth measurement values ​​of platform 1 are numbered i1=0,1,…n1;

[0055] The azimuth measurements of platform 2 are numbered i2=0,1,…n2;

[0056] The azimuth measurement values ​​of platform 3 are numbered i3=0,1,…n3;

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

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

[0059] n s This represents the number of measurement values ​​received by platform s. 0 indicates a missed detection, meaning that platform s did not detect the target.

[0060] Indicates the i-th term of platform s s Each azimuth angle measurement value;

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

[0062] The spatial spectral intensity set is composed of one measurement value taken from each platform at time k.

[0063] in,

[0064] This represents the i1th azimuth measurement value of platform 1; This represents the i2th azimuth angle measurement value of platform 2; Indicates the i-th term of platform s s Each azimuth angle measurement value; Indicates the i-th position of platform S′ S′ Each azimuth angle measurement value;

[0065] This represents the spatial spectral intensity of the i1th azimuth angle measurement on platform 1; This represents the spatial spectral intensity of the i2th azimuth angle measurement on platform 2; Indicates the i-th term of platform s s Spatial spectral intensity of each azimuth angle measurement; Indicates the i-th position of platform S′ S′ Spatial spectral intensity of each azimuth angle measurement.

[0066] Preferably, in step two, the measurement is based on each azimuth measurement set. Calculate the target position estimate for each measurement set; based on each azimuth measurement set... Based on the target location estimate, calculate the cost function value for each measurement set. Based on cost function value Construct an objective function and solve it to obtain the allocation result; use triangulation to estimate the target position based on the allocation result; use a joint probability data association method to track the target position.

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

[0068] The specific process is as follows:

[0069] Step 21: Based on each azimuth measurement set Calculate the target location estimate for each measurement set; the specific process is as follows:

[0070] Since the target position p is unknown, triangulation is used to estimate the target position.

[0071]

[0072] in, p represents the estimated target location; p represents the target location.

[0073] Represents the set of azimuth measurements The probability density function of a target p from a known location;

[0074] Step 22: Based on each azimuth measurement set and target location estimate Calculate the cost function value for each measurement set.

[0075] Step Two Three: Based on the cost function value from Step Two Two Construct the objective function and solve the objective function to obtain the allocation result;

[0076] Step 24: Use triangulation to estimate the target location based on the allocation results;

[0077] The target location is tracked using a joint probabilistic data association method;

[0078] Suppose that in a system with S′ platforms, N targets are tracked at time k-1, and the azimuth angle of the nth target tracked at time k-1 relative to the sth platform is...

[0079] Recording platform s The mean of the spatial spectral intensity of the measurements used to track the nth target up to time k-1 is denoted as . This indicates a long time window.

[0080] Preferably, in step two, the measurement set is based on each azimuth angle. and target location estimate Calculate the cost function value for each measurement set. The specific process is as follows:

[0081] 1) Azimuth measurement set From target location estimate The probability density function is:

[0082]

[0083] in,

[0084] For measurement set From target location estimate The probability density function;

[0085] Let be the detection probability of platform s;

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

[0087]

[0088] For measurement value From the target The probability density function is expressed as:

[0089]

[0090] in,

[0091] The i-th term of platform s s Individual measurement values;

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

[0093] Represents the measurement set Estimated target location;

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

[0095]

[0096] in,

[0097] Represents the measurement set The probability density function of clutter originating from space;

[0098] Φ represents clutter;

[0099] The i-th term of platform s s Individual measurement value The probability of originating from clutter;

[0100] V s This represents the area of ​​the observation region of platform s, i.e. clutter density;

[0101] 3) Calculate the cost function value for each measurement set. Given by the negative log-likelihood ratio:

[0102]

[0103] Substituting equations (13) and (16) into equation (17), we obtain the cost function value. Represented as:

[0104]

[0105] Preferably, in steps two and three, the cost function value is based on step two and two. Construct the objective function, 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, It is a binary variable;

[0111] When there are S′ measurement sets Compared with the estimated location of a certain target When no association is generated,

[0112] When there are S′ measurement sets Compared with the estimated location of a certain target When related,

[0113] The minimum value of the objective function that satisfies the constraints corresponds to This represents the allocation result.

[0114] Preferably, in step three, the azimuth interval between the tracking trajectories of any two targets relative to platform s is determined to be greater than or equal to an azimuth interval threshold. If the azimuth interval between the tracking trajectories of the two targets relative to platform s is less than the azimuth interval threshold, an observation window is opened. After the observation window is opened, platform s determines whether the azimuth measurement value at time k is within the threshold range. Inside;

[0115] If the azimuth measurement at time k is If inside, proceed to step four;

[0116] If the azimuth measurement at time k is not If k = k + 1, then execute step one;

[0117] Where ε θ This refers to the azimuth extension value of the observation window; This 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 interval between 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 one.

[0119] The specific process is as follows:

[0120] The azimuth interval between the tracking trajectories of any two targets relative to the platform s is calculated as follows:

[0121]

[0122] in,

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

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

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

[0126] If the azimuth interval is greater than or equal to the azimuth interval threshold, the observation window will not be opened, and step one will be executed.

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

[0128]

[0129] Among them, Γ θ This is the azimuth angle interval threshold;

[0130] After platform s opens the observation window, it determines whether the azimuth measurement value at time k is within the range. Inside;

[0131] If the azimuth measurement at time k is If inside, proceed to step four;

[0132] If the azimuth measurement at time k is not If k = k + 1, then execute step one;

[0133] Where ε θ This is the azimuth extension value for the observation window.

[0134] Preferably, in step four, the rate of change of the spatial spectrum intensity is calculated based on the spatial spectrum intensity corresponding to the azimuth measurement obtained in step two. The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one;

[0135] The specific process is as follows:

[0136] Calculate the rate of change of the spatial spectrum intensity based on the azimuth measurement obtained in step two. The calculation formula is as follows:

[0137]

[0138] in,

[0139] It is the i-th moment of platform s at time k. s The rate of change of the spatial spectral intensity of the tracking trajectory of the nth target with respect to each measurement value;

[0140] It is the i-th moment of platform s at time k. s Spatial spectral intensity of each measurement value;

[0141] For the time before k The average spatial spectral intensity of the tracking trajectory of the nth target at time n relative to the sth platform;

[0142] The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one;

[0143] Represented as:

[0144]

[0145] Preferably, in step five, the "suspicious merged measurement values" are 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 location is estimated based on the allocation result, and the target location is tracked.

[0146] The specific process is as follows:

[0147] Step 51: Assume that platform s contains J. s One "suspicious combined measurement value", numbered as follows:

[0148] Where, j s,1 This indicates the first "suspicious merged measurement value" number in platform s, j s,2 This indicates the second "suspicious merged measurement value" number on platform s. Indicates platform s Jth s One "suspicious combined measurement value" number;

[0149] Step 52, for J s The "suspicious merged measurement value" is "fissioned" as follows:

[0150]

[0151] in,

[0152] Indicates the measurement value Add to set Z s Z s This represents the original measurement set of platform s;

[0153] For binary decision variables;

[0154] This indicates that the j-th element in platform s will be... s Each measurement value is added to the original measurement set of platform s, which is equivalent to a fission.

[0155] This indicates that fission will not occur;

[0156]

[0157] Step 53: Let the minimum total cost obtained from equation (19) be... Calculate the minimum total cost for each combination, and select the measurement set corresponding to the minimum total cost among all combinations as the final allocation result;

[0158] Step 54: Estimate the target location based on the allocation results and track the target location.

[0159] Preferably, in step five-three, the minimum total cost obtained by equation (19) is set as follows: Calculate the minimum total cost for each combination, and select the measurement set corresponding to the minimum total cost among all combinations as the final allocation result; as shown in the following formula:

[0160]

[0161] Where min represents selecting the minimum value in the group, c min This represents the minimum value in the total cost.

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

[0163] The target location is estimated using triangulation methods based on the allocation results;

[0164] The target location is tracked using a joint probability data association method.

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

[0166] This invention proposes an underwater multi-platform multi-target tracking method based on extended multidimensional allocation. This method combines the normalized intensity of the spatial spectrum and establishes an evaluation criterion based on the lowest correlation cost, which can effectively split and merge measurement values, and improve the tracking and positioning accuracy of multiple targets based on the allocation results. Attached Figure Description

[0167] Figure 1 This is a flowchart of the present invention;

[0168] Figure 2 This is a motion map of the platform's position and the target, using a northeast-sky coordinate system with the X-axis pointing north and the Y-axis pointing east, and a custom origin position.

[0169] Figure 3 This is a location result map based on traditional multidimensional allocation;

[0170] Figure 4 The image shows the tracking results based on traditional multidimensional allocation.

[0171] Figure 5 This is a map showing the localization results based on extended multidimensional allocation;

[0172] Figure 6 This is a graph showing the tracking results based on extended multidimensional allocation. Detailed Implementation

[0173] Specific Implementation Method 1: The specific process of this implementation method for underwater multi-platform multi-target tracking based on extended multidimensional allocation is as follows:

[0174] Step 1: Set time k = 1;

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

[0176] One platform can extract any one azimuth measurement value or not extract it; if not extracted, the value is 0. The set of azimuth measurement values ​​extracted from all platforms is: This represents the i1th azimuth measurement value of platform 1; This represents the i2th azimuth angle measurement value of platform 2; Indicates the i-th term of platform s s Each azimuth angle measurement value; Indicates the i-th position of platform S′ S′ Each azimuth angle measurement value;

[0177] For each platform, one azimuth measurement value can be extracted arbitrarily or not. For all platforms, all possible combinations are taken to obtain the set of azimuth measurements corresponding to all combinations. One combination corresponds to one set of azimuth measurements.

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

[0179] Record the azimuth measurements of each platform used for each tracking trajectory and the corresponding spatial spectrum intensity;

[0180] Step 3: Determine 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 the two targets relative to platform s is less than the azimuth interval threshold, then open the observation window. After opening the observation window on platform s, determine whether the azimuth measurement value at time k is within the threshold. Inside;

[0181] If the azimuth measurement at time k is If inside, proceed to step four;

[0182] If the azimuth measurement at time k is not If k = k + 1, then execute step one;

[0183] Where ε θ This refers to the azimuth extension value of the observation window; This represents the azimuth angle 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 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 one.

[0185] Step 4: Calculate the rate of change of the spatial spectrum intensity based on the azimuth measurement obtained in Step 2. The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one;

[0186] Step 5: Split the “suspicious merged measurement values”, recombine the split measurement values, calculate the total association cost of each recombinated measurement value, select the combination with the smallest total association cost as the final allocation result, estimate the target location based on the allocation result, and track the target location.

[0187] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that: in step one, each platform undergoes conventional beamforming processing (Formula 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: Assuming each platform is equipped with a uniform linear array of M-element hydrophones, in a scenario where the uniform linear array receives a far-field target signal, the signal vector received by the linear array can be represented as:

[0189]

[0190] in,

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

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

[0193] Let m be the signal received by the m-th hydrophone in the linear array, where m = 1, 2, ..., M;

[0194] This is the signal received by the Mth hydrophone in the linear array;

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

[0196] Steps 1 and 2: Assume the signal received by the first hydrophone in 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. 2 =-1;

[0199] Since the target satisfies the far-field condition, the deviation of the incident angle of the target signal reaching each hydrophone in the linear array is negligible. That is, the angle at which the signal reaches each hydrophone is the same; only the arrival time at each hydrophone is different. Therefore, the signal received by the m-th hydrophone is...

[0200]

[0201] Where, τm Let m be the time difference between the m-th hydrophone and the 1st hydrophone receiving the far-field target signal;

[0202] Assuming the target signal does not change rapidly, the complex envelope in equation (3) can be approximated as:

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

[0204] In practical applications, the received signal itself is generally irrelevant; that is, the received signal is unrelated to e. jωt Irrelevant;

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

[0206]

[0207] Equation (5) shows that the difference in the signals received by each hydrophone in the linear array is only due to the phase difference caused by the time delay.

[0208] Step 13: Represent the signals received by each hydrophone in the linear array as a unified representation.

[0209]

[0210] in, Let τ be the direction vector, and τ2 be the time difference between the second hydrophone and the first hydrophone receiving the far-field target signal. M Let θ be the time difference between the Mth hydrophone and the first hydrophone receiving the far-field target signal, and let θ be the incident angle of the far-field target.

[0211] In practical engineering, the signals acquired by the system are discrete signals, and noise interference may exist in reality. Therefore, the linear array receiving single-target signal model is expressed as shown in equation (6):

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

[0213] Where v(n) is additive noise; s(n) is the received discrete signal; and y(n) is the discrete output signal.

[0214] Establish the time difference τ between the m-th hydrophone and the 1st hydrophone in receiving the target signal. m The relationship with the target incident angle θ is as follows:

[0215]

[0216] Where c is the speed of sound in water, and d is the hydrophone spacing;

[0217] Equation (8) shows that the time difference of the signal received by the uniform linear array is related to the hydrophone spacing d, the sound speed in the water c, and the incident angle θ. However, the hydrophone spacing and the sound speed are constant, which transforms the time difference in the array signal receiving model into a relationship with the incident angle.

[0218] Since ω=2πf=2πc / λ, and combining equation (8) with the direction vector a(θ), the direction vector in the linear array receiving single-target signal model is finally obtained, as shown in equation (9):

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

[0220] Where f is the signal frequency, λ is the wavelength, and φ is the phase;

[0221] Step 1, Section 4: In a real-world scenario, there may be multiple targets. Assume that the incident angles of N far-field targets relative to the linear array are [θ1, θ2, ..., θ]. N ], rewrite the direction vector as

[0222]

[0223] Where A(θ) is the direction matrix, a(θ1) is the direction vector of the first far-field target, a(θ2) is the direction vector of the second far-field target, and aθ N () represents the direction vector of the Nth far-field target;

[0224] ω1 is the angular frequency of the first target received, ω2 is the angular frequency of the second target received, ω N To receive the angular frequency of the Nth target;

[0225] Combining equation (7), we obtain the linear array receiving multi-target signal model, as shown in equation (11):

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

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

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

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

[0230] The azimuth measurement values ​​of platform 1 are numbered i1=0,1,…n1;

[0231] The azimuth measurements of platform 2 are numbered i2=0,1,…n2;

[0232] The azimuth measurement values ​​of platform 3 are numbered i3=0,1,…n3;

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

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

[0235] n s This represents the number of measurement values ​​received by platform s. 0 indicates a missed detection, meaning that platform s did not detect the target.

[0236] Indicates the i-th term of platform s s Each azimuth angle measurement value;

[0237] At time k, one azimuth measurement value is taken from each platform to form an azimuth measurement set. (One value is taken from each platform; it is optional, and the value is 0 if no value is taken.)

[0238] The spatial spectral intensity set is composed of one measurement value taken from each platform at time k. (One value is taken from each platform; it is optional, and the value is 0 if no value is taken.)

[0239]

[0240] in,

[0241] This represents the i1th azimuth measurement value of platform 1; This represents the i2th azimuth angle measurement value of platform 2; Indicates the i-th term of platform s s Each azimuth angle measurement value; Indicates the i-th position of platform S′ S′ Each azimuth angle measurement value;

[0242] This represents the spatial spectral intensity of the i1th azimuth angle measurement on platform 1; This represents the spatial spectral intensity of the i2th azimuth angle measurement on platform 2; Indicates the i-th term of platform s s Spatial spectral intensity of each azimuth angle measurement; Indicates the i-th position of platform S′S′ Spatial spectral intensity of each azimuth angle measurement.

[0243] The other steps and parameters are the same as in Specific Implementation Method 1.

[0244] Specific Implementation Method Three: This implementation method differs from Specific Implementation Method One or Two in that: in step two, based on each azimuth measurement set... Calculate the target position estimate for each measurement set; based on each azimuth measurement set... Based on the target location estimate, calculate the cost function value for each measurement set. Based on cost function value Construct an objective function and solve it to obtain the allocation result; use triangulation to estimate the target position (complete positioning) based on the allocation result; use joint probability data association to track the target position;

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

[0246] The SD allocation algorithm uses the minimum established cost function as the criterion, enumerates all possible combinations, calculates the association cost of each combination, and selects the combination with the lowest cost as the optimal association combination.

[0247] The specific process is as follows:

[0248] Step 21: Based on each azimuth measurement set Calculate the target location estimate for each measurement set; the specific process is as follows:

[0249] Since the target position p is unknown, least squares estimation is used instead, and triangulation is used to estimate the target position.

[0250]

[0251] in, p represents the estimated target location; p represents the target location.

[0252] Represents the set of azimuth measurements The probability density function of a target p from a known location;

[0253] Step 22: Based on each azimuth measurement set and target location estimate Calculate the cost function value for each measurement set.

[0254] Step Two Three: Based on the cost function value from Step Two Two Construct the objective function and solve the objective function to obtain the allocation result;

[0255] Step 24: Use triangulation to estimate the target position based on the allocation results (complete positioning);

[0256] The target location is tracked using a joint probabilistic data association method;

[0257] Suppose that in a system with S′ platforms, N targets are tracked at time k-1, and the azimuth angle of the nth target tracked at time k-1 relative to the sth platform is...

[0258] Recording platform s The mean of the spatial spectral intensity of the measurements used to track the nth target up to time k-1 is denoted as . This indicates a long time window.

[0259] Other steps and parameters are the same as in specific implementation method one or two.

[0260] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One to Three in that: in step two-two, based on each azimuth measurement set... and target location estimate Calculate the cost function value for each measurement set. The specific process is as follows:

[0261] 1) Azimuth measurement set From target location estimate The probability density function is:

[0262]

[0263] in,

[0264] For measurement set From target location estimate The probability density function;

[0265] Let be the detection probability of platform s;

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

[0267]

[0268] For measurement value From the target The probability density function is expressed as:

[0269]

[0270] in,

[0271] The i-th term of platform s s Individual measurement values;

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

[0273] Represents the measurement set Estimated target location;

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

[0275]

[0276] in,

[0277] Represents the measurement set The probability density function of clutter originating from space;

[0278] Φ represents clutter;

[0279] The i-th term of platform s s Individual measurement value The probability of originating from clutter;

[0280] V s This represents the area of ​​the observation region of platform s, i.e. clutter density;

[0281] 3) Calculate the cost function value for each measurement set. Given by the negative log-likelihood ratio:

[0282]

[0283] Substituting equations (13) and (16) into equation (17), we obtain the cost function value. Represented as:

[0284]

[0285] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0286] Specific Implementation Method Five: This implementation method differs from Specific Implementation Methods One to Four in that: in steps Two and Three, the cost function value based on step Two and Two is... Construct the objective function, solve the objective function to obtain the allocation result; the specific process is as follows:

[0287] The objective function is:

[0288]

[0289] Constraints:

[0290]

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

[0292] When there are S′ measurement sets Compared with the estimated location of a certain target (As given in Equation 12) when no association occurs,

[0293] When there are S′ measurement sets Compared with the estimated location of a certain target (As given in Equation 12) When associated,

[0294] The minimum value of the objective function that satisfies the constraints corresponds to This represents the allocation result.

[0295] The other steps and parameters are the same as those in one of the specific implementation methods one to four.

[0296] Specific Implementation Method Six: This implementation method differs from Specific Implementation Methods One through Five in that: in step three, the azimuth interval between the tracking trajectories of any two targets relative to platform s is determined to be smaller than the azimuth interval threshold. If the azimuth interval between the tracking trajectories of the two targets relative to platform s is less than the azimuth interval threshold, an observation window is opened. After the observation window is opened, platform s determines whether the azimuth measurement value at time k is within the threshold. Inside;

[0297] If the azimuth measurement at time k is If inside, proceed to step four;

[0298] If the azimuth measurement at time k is not If k = k + 1, then execute step one;

[0299] Where ε θ This refers to the azimuth extension value of the observation window; This represents the azimuth angle 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 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 one.

[0301] The specific process is as follows:

[0302] The azimuth interval between the tracking trajectories of any two targets relative to the platform s is calculated as follows:

[0303]

[0304] in,

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

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

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

[0308] As the azimuth interval becomes smaller, it indicates that the target may be nearby in the azimuth. When it is smaller than the system's azimuth resolution, merging will occur in the spatial spectrum.

[0309] The azimuth interval was used 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 k = k + 1 is set to execute step one;

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

[0312]

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

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

[0315] If the azimuth measurement at time k is If inside, proceed to step four;

[0316] If the azimuth measurement at time k is not If k = k + 1, then execute step one;

[0317] Where ε θ This is the azimuth extension value for the observation window.

[0318] The other steps and parameters are the same as those in one of the specific implementation methods one to five.

[0319] Specific Implementation Method Seven: This implementation method differs from Specific Implementation Methods One through Six in that: in step four, the rate of change of the spatial spectrum intensity is calculated based on the spatial spectrum intensity corresponding to the azimuth measurement value obtained in step two. The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one;

[0320] The specific process is as follows:

[0321] Calculate the rate of change of the spatial spectrum intensity based on the azimuth measurement obtained in step two. The calculation formula is as follows:

[0322]

[0323] in,

[0324] It is the i-th moment of platform s at time k. s The rate of change of the spatial spectral intensity of the tracking trajectory of the nth target with respect to each measurement value;

[0325] It is the i-th moment of platform s at time k. s Spatial spectral intensity of each measurement value;

[0326] For the time before k The average spatial spectral intensity of the tracking trajectory of the nth target at time n relative to the sth platform;

[0327] The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one;

[0328] Represented as:

[0329]

[0330] The other steps and parameters are the same as those in specific implementation methods one through six.

[0331] Specific Implementation Method Eight: This implementation method differs from Specific Implementation Methods One to Six in that: in step five, the "suspicious merged measurement values" are 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 location is estimated based on the allocation result, and the target location is tracked.

[0332] The specific process is as follows:

[0333] Step 51: Assume that platform s contains J. s One "suspicious combined measurement value", numbered as follows:

[0334] Where, j s,1 This indicates the first "suspicious merged measurement value" number in platform s, j s,2 This indicates the second "suspicious merged measurement value" number on platform s. Indicates platform s Jth s One "suspicious combined measurement value" number;

[0335] Step 52, for J s The "suspicious merged measurement value" is "fissioned" as follows:

[0336]

[0337] in, Indicates the measurement value Add to set Z s Z s This represents the original measurement set of platform s (platform s has its own original measurement set);

[0338] For binary decision variables;

[0339] This indicates that the j-th element in platform s will be... s Each measurement value is added to the original measurement set of platform s, which is equivalent to a fission.

[0340] This indicates that fission will not occur;

[0341]

[0342] Step 53: Let the minimum total cost obtained from equation (19) be... Calculate the minimum total cost for each combination, and select the measurement set corresponding to the minimum total cost among all combinations as the final allocation result;

[0343] Step 54: Estimate the target location based on the allocation results and track the target location.

[0344] The other steps and parameters are the same as those in specific implementation methods one through seven.

[0345] Specific Implementation Method Nine: This implementation method differs from Specific Implementation Methods One through Eight in that: in step five-three, the minimum total cost obtained by equation (19) is set as... Calculate the minimum total cost for each combination, and select the measurement set corresponding to the minimum total cost among all combinations as the final allocation result; as shown in the following formula:

[0346]

[0347] Where min represents selecting the minimum value in the group, c min This represents the minimum value in the total cost.

[0348] The other steps and parameters are the same as those in specific implementation methods one through eight.

[0349] Specific Implementation Method Ten: This implementation method differs from Specific Implementation Methods One through Nine in 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 as follows:

[0350] The target location is estimated using triangulation based on the allocation results (location is completed);

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

[0352] The other steps and parameters are the same as those in specific implementation methods one through nine.

[0353] The beneficial effects of the present invention are verified using the following embodiments:

[0354] Example 1:

[0355] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0356] First, the platform position and the target's trajectory are constructed. This example includes 4 targets and 3 targets, as detailed below. Figure 2As shown. The positions of the four platforms are (300m, 350m), (1100m, 50m), (1900m, 50m), and (2700m, 350m), respectively, and the initial positions of the three targets are (1100m, 1700m), (2000m, 1750m), and (2000m, 17500m), respectively. The positioning results based on traditional multidimensional allocation are as follows. Figure 3 As shown, target 2 has a significant number of missing measurement values, and the tracking results are as follows. Figure 4 As shown, target 2 experienced a tracking interruption. The localization result based on extended multidimensional allocation is as follows... Figure 5 As shown, the corresponding tracking results are as follows: Figure 6 As shown, the results indicate that the proposed method has good tracking performance.

[0357] This invention may have other embodiments. Without departing from the spirit and essence of this invention, those skilled in the art can make various corresponding changes and modifications according to this invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A method for underwater multi-platform multi-target tracking based on extended multidimensional allocation, characterized in that: The specific process of the method is as follows: Step 1: Set time k = 1; Each platform undergoes conventional beamforming processing to extract the azimuth measurement value and the corresponding spatial spectrum intensity of the platform; Step 2: Based on each azimuth measurement set Calculate the target location estimate for each measurement set; Based on each azimuth measurement set Based on the target location estimate, calculate the cost function value for each measurement set. Based on cost function value Construct the objective function and solve the objective function to obtain the allocation result; The target location is estimated using triangulation methods based on the allocation results; The target location is tracked using a joint probabilistic data association method; Record the azimuth measurements of each platform used for each tracking trajectory and the corresponding spatial spectrum intensity; Step 3: Determine 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 the two targets relative to platform s is less than the azimuth interval threshold, then open the observation window. After opening the observation window on platform s, determine whether the azimuth measurement value at time k is within the threshold. Inside; If the azimuth measurement at time k is If inside, proceed to step four; If the azimuth measurement at time k is not If k = k + 1, then execute step one; Where ε θ This refers to the azimuth extension value of the observation window; This represents the azimuth angle of the tracking trajectory of the i-th target relative to the s-th platform at time k-1; If the azimuth interval between 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 one. Step 4: Calculate the rate of change of the spatial spectrum intensity based on the azimuth measurement obtained in Step 2. The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one; Step 5: Split the "suspicious merged measurement values", recombine the split measurement values, calculate the total association cost of each recombinated measurement value, select the combination with the smallest total association cost as the final allocation result, estimate the target location based on the allocation result, and track the target location.

2. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 1, characterized in that: In step one, each platform undergoes conventional beamforming processing to extract the azimuth measurement values ​​and corresponding spatial spectrum intensities of the platform; the specific process is as follows: Step 11: Assume each platform is equipped with a uniform linear array of M-element hydrophones. The signal vector received by the linear array is represented as... in, This is the signal received by the first hydrophone in the linear array; This is the signal received by the second hydrophone in the linear array; Let m be the signal received by the m-th hydrophone in the linear array, where m = 1, 2, ..., M; This is the signal received by the Mth hydrophone in the linear array; The signal received by the linear array M-element hydrophone; T represents the transpose; t represents time; Steps 1 and 2: Assume the signal received by the first hydrophone in the linear array is: 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. 2 =-1; The signal received by the m-th hydrophone is Where, τ m Let m be the time difference between the m-th hydrophone and the 1st hydrophone receiving the target signal; The complex envelope in equation (3) can be approximated as: s(t-τ m )≈s(t) (4) Combining equations (3) and (4), the signal received by the m-th hydrophone is approximately: Step 13: Represent the signals received by each hydrophone in the linear array as a unified representation. in, Let τ be the direction vector, and τ2 be the time difference between the second hydrophone and the first hydrophone receiving the target signal. M Let θ be the time difference between the Mth hydrophone and the first hydrophone receiving the target signal, and let θ be the incident angle of the target. The linear array receiving single-target signal model is expressed 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; and y(n) is the discrete output signal. Establish the time difference τ between the m-th hydrophone and the 1st hydrophone in receiving the target signal. m The relationship with the target incident angle θ is as follows: Where c is the speed of sound in water, and d is the hydrophone spacing; Since ω=2πf=2πc / λ, and combining equation (8) with the direction vector a(θ), the direction vector in the linear array receiving single-target signal model is finally obtained, 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 1, Section 4: Assume the incident angles of N targets relative to the linear array are [θ1, θ2, ..., θ]. N ], rewrite the direction vector as 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, and a(θ3) is the direction vector of the third target. N () represents the direction vector of the Nth target; ω1 is the angular frequency of the first target received, ω2 is the angular frequency of the second target received, ω N To receive the angular frequency of the Nth target; Combining equation (7), we obtain 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 azimuth measurement value and spatial spectral intensity of the platform from each output signal y(n); the specific process is as follows: For S′ platforms equipped with passive sonar, assume the platform numbers are s = 1, 2, ..., S′; The azimuth measurement values ​​of platform 1 are numbered i1=0,1,…n1; The azimuth measurements of platform 2 are numbered i2=0,1,…n2; The azimuth measurement values ​​of platform 3 are numbered i3=0,1,…n3; The azimuth measurement value of platform s is numbered i s =0,1,…,n s ; The azimuth measurement value of platform S′ is numbered i S′ =0,1,…,n S′ ; n s This represents the number of measurement values ​​received by platform s. 0 indicates a missed detection, meaning that platform s did not detect the target. Indicates the i-th term of platform s s Each azimuth angle measurement value; At time k, one azimuth measurement value is taken from each platform to form an azimuth measurement set. The spatial spectral intensity set is composed of one measurement value taken from each platform at time k. in, This represents the i1th azimuth measurement value of platform 1; This represents the i2th azimuth angle measurement value of platform 2; Indicates the i-th term of platform s s Each azimuth angle measurement value; Indicates the i-th position of platform S′ S′ Each azimuth angle measurement value; This represents the spatial spectral intensity of the i1th azimuth angle measurement on platform 1; This represents the spatial spectral intensity of the i2th azimuth angle measurement on platform 2; Indicates the i-th term of platform s s Spatial spectral intensity of each azimuth angle measurement; Indicates the i-th position of platform S′ S′ Spatial spectral intensity of each azimuth angle measurement.

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

4. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 3, characterized in that: In step two, the measurement set for each azimuth angle is used. and target location estimate Calculate the cost function value for each measurement set. The specific process is as follows: 1) Azimuth measurement set From target location estimate The probability density function is: in, For measurement set From target location estimate The probability density function; Let be the detection probability of platform s; u(i s ) is an indicator function; expressed as: For measurement value From the target The probability density function is expressed as: in, The i-th term of platform s s Individual measurement values; σ s This represents the noise of platform s; Represents the measurement set Estimated target location; 2) Assuming the clutter is uniformly distributed in space, then the measurement set The probability density function of clutter from space is: in, Represents the measurement set The probability density function of clutter originating from space; Φ represents clutter; The i-th term of platform s s Individual measurement value The probability of originating from clutter; V s This represents the area of ​​the observation region of platform s, i.e. clutter density; 3) Calculate the cost function value for each measurement set. Given by the negative log-likelihood ratio: Substituting equations (13) and (16) into equation (17), we obtain the cost function value. Represented as:

5. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 4, characterized in that: The cost function value in step two and three is based on step two and two. Construct the objective function, solve the objective function to obtain the allocation result; the specific process is as follows: The objective function is: Constraints: In the formula, It is a binary variable; When there are S′ measurement sets Compared with the estimated location of a certain target When no association is generated, When there are S′ measurement sets Compared with the estimated location of a certain target When related, The minimum value of the objective function that satisfies the constraints corresponds to This represents the allocation result.

6. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 5, characterized in that: In step three, the azimuth interval between the tracking trajectories of any two targets relative to platform s is compared to an azimuth interval threshold. If the azimuth interval between the tracking trajectories of the two targets relative to platform s is less than the azimuth interval threshold, an observation window is opened. After the observation window is opened, platform s determines whether the azimuth measurement value at time k is within the threshold. Inside; If the azimuth measurement at time k is If inside, proceed to step four; If the azimuth measurement at time k is not If k = k + 1, then execute step one; Where ε θ This refers to the azimuth extension value of the observation window; This represents the azimuth angle of the tracking trajectory of the i-th target relative to the s-th platform at time k-1; If the azimuth interval between 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 one. The specific process is as follows: The azimuth interval between the tracking trajectories of any two targets relative to the platform s is calculated as follows: in, This represents the azimuth angle interval between the tracking trajectories of the i-th and j-th targets relative to the platform s; This represents the azimuth angle of the tracking trajectory of the i-th target relative to the s-th platform at time k-1; This represents the azimuth angle of the tracking trajectory of the j-th target at time k-1 relative to the s-th 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 one; If the azimuth interval is less than the azimuth interval threshold, then the two targets corresponding to the azimuth interval are determined to enter the proximity warning, and platform s opens the observation window: Among them, Γ θ This is the azimuth angle interval threshold; After platform s opens the observation window, it determines whether the azimuth measurement value at time k is within the range. Inside; If the azimuth measurement at time k is If inside, proceed to step four; If the azimuth measurement at time k is not If k = k + 1, then execute step one; Where ε θ This is the azimuth extension value for the observation window.

7. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 6, characterized in that: In step four, the rate of change of the spatial spectrum intensity is calculated based on the spatial spectrum intensity corresponding to the azimuth measurement obtained in step two. The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one; The specific process is as follows: Calculate the rate of change of the spatial spectrum intensity based on the azimuth measurement obtained in step two. The calculation formula is as follows: in, It is the i-th moment of platform s at time k. s The rate of change of the spatial spectral intensity of the tracking trajectory of the nth target with respect to each measurement value; It is the i-th moment of platform s at time k. s Spatial spectral intensity of each measurement value; For the time before k The average spatial spectral intensity of the tracking trajectory of the nth target at time n relative to the sth platform; The rate of change of spatial spectral intensity With intensity change threshold Γ A Comparison, if the rate of change of spatial spectral intensity Greater than the intensity change threshold Γ A Then the corresponding measurement value is taken as a "suspected merged measurement value"; if the rate of change of spatial spectral intensity Less than or equal to the intensity change threshold Γ A Then let k = k + 1 and execute step one; Represented as:

8. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 7, characterized in that: In step five, the "suspicious merged measurement value" is split, the split measurement value is 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 location is estimated based on the allocation result, and the target location is tracked. The specific process is as follows: Step 51: Assume that platform s contains J. s One "suspicious combined measurement value", numbered as follows: Where, j s,1 This indicates the first "suspicious merged measurement value" number in platform s, j s,2 This indicates the second "suspicious merged measurement value" number on platform s. Indicates platform s Jth s One "suspicious combined measurement value" number; Step 52, for J s The "suspicious merged measurement value" is "fissioned" as follows: in, Indicates the measurement value Add to set Z s Z s This represents the original measurement set of platform s; For binary decision variables; This indicates that the j-th element in platform s will be... s Each measurement value is added to the original measurement set of platform s, which is equivalent to performing a fission. This indicates that fission will not occur; Step 53: Let the minimum total cost obtained from equation (19) be... Calculate the minimum total cost for each combination, and select the measurement set corresponding to the minimum total cost among all combinations as the final allocation result; Step 54: Estimate the target location based on the allocation results and track the target location.

9. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 8, characterized in that: In step five-three, the lowest total cost obtained by equation (19) is: Calculate the minimum total cost for each combination, and select the measurement set corresponding to the minimum total cost among all combinations as the final allocation result; As shown in the following formula: Where min represents selecting the minimum value in the group, c min This represents the minimum value in the total cost.

10. The underwater multi-platform multi-target tracking method based on extended multidimensional allocation according to claim 9, characterized in 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 as follows: The target location is estimated using triangulation methods based on the allocation results; The target location is tracked using a joint probability data association method.

Images