A method of estimating the depth of a moving sound source

Abstract

The present application relates to a kind of motion sound source depth estimation method, the method comprises: simulation is to be monitored sea area with acoustic propagation model, determine the depth of undetermined target sound source, select frequency;Using normal mode theory, the phase velocity of each order normal mode, eigenfunction and the normal mode energy distribution of each undetermined target sound source depth are calculated under selected frequency;Process the receiving signal of selected frequency in each section observation time period, obtain the power curve of compensated output beam;For any undetermined target sound source depth, by cyclic iteration, the optimal value of the phase velocity of each order normal mode is obtained, and the optimal error of the undetermined target sound source depth is obtained according to the beam power curve reconstructed by normal mode energy distribution and the optimal value of the phase velocity of each order normal mode;The undetermined target sound source depth with minimum optimal error is used as target sound source depth.The method of the present application uses the low-frequency single-frequency signal radiated by the target sound source, and realizes the depth estimation of the moving target in the shallow sea with constant environment.

Description

Technical Field

[0001] This invention relates to the field of marine acoustics, and in particular to a method for estimating the depth of a moving sound source. Background Technology

[0002] Underwater vehicles are crucial to the safety of maritime traffic. Target depth information is a critical parameter. Providing reliable target depth estimates based on signals received by acoustic sensors is an effective way to resolve / estimate the depth of underwater targets. Matched field processing is a classic method for target depth estimation. This method uses a sound propagation model to calculate a copy field vector, then matches the copy field with the measured field. However, due to the lack of prior environmental information, the copy field is often inaccurately constructed, leading to the failure of the matched field processing method. Furthermore, the sound source depth estimated by the matched field method originates from the range-depth ambiguity plane; errors in range estimation can also cause deviations in the sound source depth. Summary of the Invention

[0003] In view of this, the main objective of the present invention is to provide a method for estimating the depth of a moving sound source, which uses the mid-to-low frequency single-frequency signal radiated by the target sound source to achieve depth estimation of a moving target in a shallow sea environment with constant horizontal direction.

[0004] To achieve the above objectives, this application provides a method for estimating the depth of a moving sound source, comprising:

[0005] An acoustic propagation model is used to simulate the sea area to be monitored. The position of the horizontal array of hydrophones is marked in the simulated sea area. The depth range to be detected in the simulated sea area is divided into multiple depth grids. For any depth grid, a depth is selected from that depth grid as the depth of the target sound source. Using normal mode theory, it is determined that there is a frequency band of no less than two normal modes in the simulated sea area. From the intersection of this frequency band and the frequency band of the target sound source, a frequency is randomly selected as the selected frequency. Using normal mode theory, the phase velocity and eigenfunction of each normal mode at the selected frequency in the simulated sea area are calculated. Based on the phase velocity and eigenfunction of the normal modes, the normal mode energy distribution at the depth of each target sound source is calculated.

[0006] The received signals at selected frequencies during each observation period of the hydrophone horizontal array are processed using MVDR to obtain the azimuth of the target source during each observation period. An azimuth is randomly selected from the calculated target source azimuth as a reference azimuth. The phase velocity spectrum scanning range is determined by the minimum and maximum values ​​of the phase velocities of each normal mode calculated through simulation. The phase velocity spectrum scanning range is divided into multiple phase velocity grids. A phase velocity is selected from each phase velocity grid. For any selected phase velocity, the following steps are performed: using this phase velocity as the reference velocity, the phase of the received signal during each observation period is compensated to the reference azimuth. The compensated covariance matrix and the compensated weighting vector are calculated, and the compensated output beam power is calculated to obtain the compensated output beam power curve.

[0007] For any unknown target sound source depth, the following steps are performed: using the simulated phase velocities of each normal mode as the initial values ​​of each normal mode phase velocity, performing an optimal value search for each normal mode phase velocity and updating each normal mode phase velocity as one iteration, and repeating the iteration until the number of iterations is greater than the maximum number of iterations or the error change is less than the convergence threshold, at which point the iteration stops, and the error of the last iteration is recorded as the optimal error for the unknown target sound source depth;

[0008] The depth of the undetermined target sound source with the smallest optimal error is taken as the target sound source depth.

[0009] The error of one iteration is the straight-line distance between the compensated output beam power curve and the reconstructed beam power curve. The reconstructed beam power curve is obtained by reconstructing the compensated weighted vector, the normal mode energy distribution, and the phase velocities of each normal mode after the update of this iteration.

[0010] In one possible implementation, the energy distribution of each normal mode at any depth of a target sound source under a selected frequency is calculated, as expressed by the formula:

[0011]

[0012] in, Let represent the energy of the k-th normal mode at the depth z of the undetermined target sound source; The initial phase velocity of the k-th normal mode obtained from simulation calculation; z r This indicates the depth of the horizontal array of hydrophones. Let f be the eigenfunction of the k-th normal mode of z; f is the selected frequency.

[0013] right Normalized to:

[0014]

[0015] The energy distributions of the normal modes of z are obtained as follows:

[0016]

[0017] Where K is the normal mode order.

[0018] In another possible implementation, the calculated compensated covariance matrix is ​​expressed as follows:

[0019]

[0020] p'(s n ,c l )=p(s n )⊙b(c l ,θ n )

[0021]

[0022] in, Indicated by c l The covariance matrix after compensation for the reference velocity; p(s) n ) represents the s-th frame of data in the nth observation time period; c l This indicates that the phase velocity obtained from the l-th phase velocity grid is used as the reference velocity; θ n This indicates the azimuth of the target signal source during the nth observation time period; θ0 is the reference azimuth; there are a total of N observation time periods, and the received signal during each observation time period is divided into S beats; f is the selected frequency; v(θ n )=-[cosθ n sinθ n ] T It refers to the direction θ n The unit direction vector; δ m Let m be the position of the m-th array element relative to the reference array element.

[0023] In another possible implementation, the formula for calculating the compensated weighted vector is expressed as:

[0024] w(c l )=[w1(c l ),w2(c l ),...,w M (c l )] T

[0025]

[0026] Among them, c l This indicates that the phase velocity obtained from the l-th phase velocity grid is used as the reference velocity; w(c l ) indicates c lThe compensated weighted vector is calculated from the reference velocity; θ0 is the reference orientation; m represents the m-th element; r m =[r xm ,r ym ] T It represents the two-dimensional coordinates of the m-th array element.

[0027] In another possible implementation, completing one iteration for any desired target sound source depth is as follows:

[0028] The search range for the phase velocity of each normal mode is determined as follows:

[0029]

[0030] Where dc is the preset minimum interval of phase velocity; Δc is the maximum search interval; c′ represents the phase velocity of the current k-th normal mode. k Let be the search phase velocity of the current k-th normal mode;

[0031] We construct the cost functions for each normal mode and calculate the optimal value of the phase velocity of the k-th normal mode, expressed as:

[0032]

[0033]

[0034]

[0035]

[0036]

[0037] in, B(c) represents the optimal value of the phase velocity of the k-th normal mode found by the search; B(c) is the output beam power curve after compensation. These are signal energy and noise energy, respectively. This represents the normalized value of the energy of the k-th normal mode at the depth z of the undetermined target sound source; This represents the energy distribution of the normal modes at depth z of the target sound source, excluding the k-th normal mode; Ψ1(c′) k ) represents the beam response curve of the k-th order search normal mode; Ψk represents the steering vector of the k-th search normal mode; Ψ2 represents the beam response curve of the other normal modes besides the k-th normal mode; Ψk * The beam response curve represents the current normal mode; L represents the phase velocity grid number; Γ represents the beam response curve for noise. The steering vector represents the current normal mode; θ0 is the reference azimuth; δm This represents the position of the m-th array element relative to the reference array element.

[0038] Update the current normal wave phase velocity with the optimal values ​​of the phase velocities of each normal wave found, and update the beam response curve of the current normal wave simultaneously; calculate the difference between the current iteration error and the previous iteration error, i.e., the error change, and determine that the error change is not less than the convergence threshold and the number of iterations is not greater than the maximum number of iterations. Use the current iteration error as the previous iteration error and increment the iteration count by 1.

[0039] In another possible implementation, for any unknown target sound source depth, the error formula for the first iteration is expressed as:

[0040]

[0041]

[0042]

[0043] Among them, G new This represents the error of that single iteration; B(c) represents the current k-th normal wave phase velocity after updating the phase velocities of each normal wave in this iteration; B(c) represents the compensated output beam power curve. These are signal energy and noise energy, respectively. Ψ1 represents the beam response curve of the current k-th normal mode after the phase velocity of each normal mode is updated in this iteration; Ψ2 represents the beam response curve of the other normal modes except the current k-th normal mode after the phase velocity of each normal mode is updated in this iteration. This represents the normalized value of the energy of the k-th normal mode at the depth z of the undetermined target sound source; Ψ represents the energy distribution of the normal modes at depth z of the target sound source, excluding the k-th normal mode; * Γ represents the beam response curve of the current normal mode after updating the phase velocities of each normal mode in this iteration; Γ represents the beam response curve of the noise. This represents the steering vector of the current normal mode after the phase velocities of each normal mode are updated in this iteration.

[0044] Based on the above, the method for estimating the depth of a moving sound source provided by the present invention has the following advantages and features:

[0045] 1. By minimizing the phase velocity spectrum residual, the depth of the sound source is estimated, which to some extent overcomes the energy leakage of normal modes caused by insufficient array aperture. It can be applied to the depth estimation of shallow sea motion sound sources with shorter aperture arrays and higher frequencies, and overcomes the problem of depth estimation of shallow sea target sound sources when the aperture of the horizontal array is relatively small.

[0046] 2. By combining modal filtering to estimate the normal mode energy with modal matching to estimate the source depth, the overfitting caused by too many unknowns in the reconstruction process is overcome, making the source depth estimation results more robust. Attached Figure Description

[0047] Figure 1 This is a flowchart illustrating a method for estimating the depth of a moving sound source according to an embodiment of the present invention.

[0048] Figure 2 This is a schematic diagram of the sound velocity profile of seawater and sedimentary parameters in a certain sea area.

[0049] Figure 3 This is a schematic diagram showing the location of a horizontal array of hydrophones deployed on the seabed.

[0050] Figure 4 This is a schematic diagram comparing the estimation results of sound sources at different depths using the present invention and traditional modal filtering methods;

[0051] Figure 5 This is a schematic diagram comparing the target's estimated azimuth trajectory and the actual target range as presented in this invention.

[0052] Figure 6 This invention provides a comparison between the estimated depths of different sound sources and the actual depths. Detailed Implementation

[0053] Specifically, the flowchart of a method for estimating the depth of a moving sound source according to an embodiment of the present invention is as follows: Figure 1 As shown, steps 101 to 104 are included.

[0054] Step 101: Simulate the sea area to be monitored using an acoustic propagation model. Mark the position of the horizontal array of the hydrophone in the simulated sea area. Divide the depth range to be detected in the simulated sea area into multiple depth grids. For any depth grid, select a depth from that depth grid as the depth of the target sound source to be determined. Using normal mode theory, determine the frequency band range of at least two normal modes in the simulated sea area. From the intersection of this frequency band range and the target sound source frequency band, arbitrarily select a frequency as the selected frequency. Using normal mode theory, calculate the phase velocity and eigenfunction of each normal mode at the selected frequency in the simulated sea area, and calculate the normal mode energy distribution at the depth of each target sound source to be determined based on the normal mode phase velocity and eigenfunction.

[0055] Step 102: Use MVDR to process the received signals of selected frequencies for each observation period of the hydrophone horizontal array to obtain the azimuth of the target source for each observation period. Randomly select one azimuth from the calculated target source azimuth as the reference azimuth. Determine the phase velocity spectrum scanning range based on the minimum and maximum values ​​of the phase velocities of each normal mode calculated by simulation, and divide the phase velocity spectrum scanning range into multiple phase velocity grids. Select one phase velocity from each phase velocity grid, and perform the following for any selected phase velocity: using the phase velocity as the reference velocity, compensate the phase of the received signal for each observation period to the reference azimuth, calculate the compensated covariance matrix and the compensated weighting vector, and calculate the compensated output beam power to obtain the compensated output beam power curve.

[0056] Step 103: For any target sound source depth to be determined, perform the following: Use the phase velocities of each normal mode calculated by simulation as the initial values ​​of each normal mode phase velocity, perform an optimal value search for each normal mode phase velocity, and update each normal mode phase velocity as one iteration. Iterate repeatedly until the number of iterations is greater than the maximum number of iterations or the error change is less than the convergence threshold, and stop iterating. Record the error of the last iteration as the optimal error of the target sound source depth to be determined.

[0057] Step 104: Take the depth of the undetermined target sound source with the smallest optimal error as the target sound source depth.

[0058] The error of one iteration is the straight-line distance between the compensated output beam power curve and the reconstructed beam power curve. The reconstructed beam power curve is obtained by reconstructing the compensated weighted vector, the normal mode energy distribution, and the phase velocities of each normal mode after the update of this iteration.

[0059] Here, in step 101, the simulation of the sea area to be monitored using an acoustic propagation model is based on the waveguide environmental parameters of the sea area to be monitored; the waveguide environmental parameters include: water depth, sound velocity profile and seabed acoustic parameters.

[0060] The depth range to be detected is determined based on the maximum diving depth of the target, which is the maximum diving depth from sea level to the target. The depth grid ranges from 1 to 2 meters. Selecting a depth from each depth grid as the target sound source depth can be done using the upper or lower boundary values, or an intermediate value; here, the upper boundary value is used. The selected frequency is no higher than 200 Hz.

[0061] The formula for calculating the energy distribution of each normal mode wave at the depth of any undetermined target sound source at a selected frequency is expressed as follows:

[0062]

[0063] in, Let represent the energy of the k-th normal mode at the depth z of the undetermined target sound source; The initial phase velocity of the k-th normal mode obtained from simulation calculation; z r This indicates the depth of the horizontal array of hydrophones. Let f be the eigenfunction of the k-th normal mode of z; f is the selected frequency.

[0064] right Normalized to:

[0065]

[0066] The energy distributions of the normal modes of z are obtained as follows:

[0067]

[0068] Where K is the normal mode order.

[0069] In step 102, the location of the target information source obtained in each observation time period is:

[0070] Perform the following on the received signal for any given observation period:

[0071] The received signal during the observation period is divided into S equal parts. Fourier transforms are performed on the s-th part of the data to obtain the s-th beat frequency domain signal, expressed as:

[0072] p (s) =[p s1 ,p s2 ,...,p sM ] T

[0073] Where M is the number of array elements, (·) T This indicates transpose.

[0074] The covariance matrix is ​​obtained from the frequency domain signal, and the formula is expressed as:

[0075]

[0076] in,(·) H This indicates the complex conjugate transpose.

[0077] The weighted vector formula is expressed as:

[0078] w(θ)=[w1(θ),w2(θ),...,w M (θ)] T

[0079]

[0080] Where θ represents the array scanning angle; f is the selected frequency; τm (θ) represents the time delay of the m-th array element relative to the reference array element, given by τ. m (θ)=v T (θ)r m / c ref Calculated, c ref The reference speed of sound is v(θ) = -[cosθ,sinθ]. T It is the unit direction vector pointing to the azimuth θ, r m =[r xm ,r ym ] T These are the two-dimensional coordinates of the array element.

[0081] Therefore, by using MVDR to process the received signal during this observation period, the power spectrum of the output beam is:

[0082]

[0083] The θ corresponding to max(B(θ)) is the azimuth of the target source during that observation period.

[0084] The phase velocity spectrum scanning range is determined by using the minimum and maximum values ​​of the simulated normal wave phase velocities as follows: the minimum value minus 20-25 is used as the minimum value of the scanning range, and the maximum value plus 20-25 is used as the maximum value of the scanning range. The phase velocity grid is no greater than half of the minimum phase velocity interval; the phase velocity interval refers to the interval between two adjacent phase velocities in the simulated normal wave phase velocities. Taking a phase velocity from each phase velocity grid can be done using the upper or lower boundary values ​​of each grid, or an intermediate value; here, the intermediate value of each grid is used.

[0085] The formula for compensating the phase of the received signal to the reference azimuth during each observation time period, using any phase velocity as the reference velocity, is expressed as follows:

[0086] p'(s n ,c l )=p(s n )⊙b(c l ,θ n )

[0087]

[0088] Where p'(s) n ,c l p(s) represents the phase compensation data of the s-th frame of the nth observation time period, adjusted to the reference azimuth; n ) represents the s-th frame of data in the nth observation period; θ0 represents the reference azimuth; c l v(θ) represents the phase velocity taken from the l-th phase velocity grid.n )=-[cosθ n sinθ n ] T It refers to the direction θ n The unit direction vector; f is the selected frequency; ⊙ represents the Hadamard product.

[0089] The formula for the compensated covariance matrix is ​​as follows:

[0090]

[0091] in, Indicated by c l The covariance matrix after compensation for the reference velocity; N is the number of observation time periods; S is the number of frames in any observation time period; s N This represents the s-th observation in the Nth observation time period. N shoot.

[0092] After azimuth correction, the weighted vectors of the beamforming all point to the reference azimuth. These compensated weighted vectors are linear and used to calculate the beam response curve and noise response curve, resulting in low computational complexity and saving computational resources. The formula is expressed as:

[0093] w(c l )=[w1(c l ),w2(c l ),...,w M (c l )] T

[0094]

[0095] Among them, c l This indicates that the phase velocity obtained from the l-th phase velocity grid is used as the reference velocity; w(c l ) indicates c l The compensated weighted vector is calculated from the reference velocity; θ0 is the reference orientation; m represents the m-th element; r m =[r xm ,r ym ] T It represents the two-dimensional coordinates of the m-th array element.

[0096] The formula for the compensated output beam power is expressed as follows:

[0097]

[0098] Wherein, c represents l The output beam power after compensation for the reference speed.

[0099] In one possible implementation, in step 103, the maximum number of iterations is 100; the convergence threshold ranges from 6 to 10.

[0100] The process of completing one iteration for any unknown target sound source depth is as follows:

[0101] The search range for the phase velocity of each normal mode is determined as follows:

[0102]

[0103] Where dc is the preset minimum phase velocity interval, and the phase velocity grid ≤ dc ≤ minimum phase velocity interval. Here, the minimum phase velocity interval refers to the minimum value of the interval between two adjacent phase velocities in the normal wave phase velocity obtained by simulation calculation; Δc is the maximum search interval, with a value range of 10 to 15 m. c′ represents the phase velocity of the current k-th normal mode. k Let be the search phase velocity of the current k-th normal mode.

[0104] We construct the cost functions for each normal mode and use convex optimization methods, such as CVX, to calculate the optimal value of the phase velocity of the k-th normal mode, expressed as:

[0105]

[0106]

[0107]

[0108]

[0109]

[0110] in, The optimal value of the phase velocity of the k-th normal mode is found; B(c) = {B(c1), ..., B(c2)} L )} represents the output beam power curve after compensation; These are signal energy and noise energy, respectively. This represents the normalized value of the energy of the k-th normal mode at the depth z of the undetermined target sound source; This represents the energy distribution of the normal modes at depth z of the target sound source, excluding the k-th normal mode; Ψ1(c′) k ) represents the beam response curve of the k-th order search normal mode; Ψk represents the steering vector of the k-th search normal mode; Ψ2 represents the beam response curve of the other normal modes besides the k-th normal mode; Ψk *The beam response curve represents the current normal mode; L represents the phase velocity grid number; Γ represents the beam response curve for noise. The steering vector represents the current normal mode; θ0 is the reference azimuth; δ m This represents the position of the m-th array element relative to the reference array element.

[0111] Update the current normal wave phase velocity with the optimal values ​​of the phase velocities of each normal wave found, and update the beam response curve of the current normal wave simultaneously; calculate the difference between the current iteration error and the previous iteration error, i.e., the error change, and determine that the error change is not less than the convergence threshold and the number of iterations is not greater than the maximum number of iterations. Use the current iteration error as the previous iteration error and increment the iteration count by 1.

[0112] Here, the principle of reconstructing the beam power curve is as follows:

[0113] The signal received by the hydrophone array is

[0114]

[0115] in Q is a function proportional to the sound source spectrum; z s For the target sound source depth, z r The depth of the horizontal array of the hydrophone. If the observation matrix G is known, the vector x, representing the distribution of normal mode energy, can be obtained by solving a system of linear equations. However, the observation matrix G and the sampling distance (r1, r2, ..., r...) are related. M Related to this. Assuming the target is located in the end-fire direction of the horizontal array, and the array aperture is sufficiently small relative to the distance from the target to the array, the following approximate relationship holds:

[0116]

[0117]

[0118] r0 is the distance from the target to the reference element. Therefore,

[0119]

[0120]

[0121]

[0122] Where a k It is the complex amplitude of the normal mode wave. δ m For the position of the m-th array element relative to the reference array element, the new observation matrix... It no longer depends on the sampling distance. Under the assumptions of a constant waveguide environment level, adiabatic approximation, and spatial white noise, when the number of propagating normal modes can be determined, the array covariance matrix can be written as:

[0123]

[0124] in

[0125]

[0126] The off-diagonal elements of matrix Λ are

[0127]

[0128] diagonal element is

[0129]

[0130] in

[0131]

[0132] Under normal circumstances

[0133]

[0134] When the target's movement distance exceeds the normal mode interference span, the off-diagonal terms will gradually approach zero over time. Therefore, when the sound source depth does not change drastically with time and the sound source's movement distance is sufficiently far, the normal mode coherence terms smoothed through the distance dimension approach zero, while the incoherence terms gradually accumulate, achieving decoherence. Matrix Λ can be approximated as a diagonal matrix.

[0135]

[0136] The beam output can be further expressed as

[0137]

[0138] yes The k-th column. The reconstructed beam power curve is

[0139]

[0140] For any unknown target sound source depth, the error formula for the first iteration is expressed as:

[0141]

[0142]

[0143]

[0144]

[0145]

[0146]

[0147] Among them, G new This represents the error of that single iteration; B(c) represents the current k-th normal wave phase velocity after updating the phase velocities of each normal wave in this iteration; B(c) represents the compensated output beam power curve. These are signal energy and noise energy, respectively. Ψ1 represents the beam response curve of the current k-th normal mode after the phase velocity of each normal mode is updated in this iteration; Ψ2 represents the beam response curve of the other normal modes except the current k-th normal mode after the phase velocity of each normal mode is updated in this iteration. This represents the normalized value of the energy of the k-th normal mode at the depth z of the undetermined target sound source; Ψ represents the energy distribution of the normal modes at depth z of the target sound source, excluding the k-th normal mode; * Γ represents the beam response curve of the current normal mode after updating the phase velocities of each normal mode in this iteration; Γ represents the beam response curve of the noise. This represents the steering vector of the current normal mode after updating the phase velocities of each normal mode in this iteration; θ0 is the reference azimuth; δ m Let m be the position of the m-th array element relative to the reference array element.

[0148] Application examples:

[0149] Assuming target detection is being conducted in a certain sea area, the seawater sound velocity profile and sedimentary layer parameters are as follows: Figure 2 As shown, only the parameters of the sedimentary layer are listed, without plotting the corresponding thickness or labeling the parameters. The sound velocity C at the interface of the sedimentary layer... s The sound velocity C at the interface below the sedimentary layer is 1608 m / s. b The velocity is 1800 m / s, and the thickness of the sedimentary layer is H. s The distance is 200m, the sound velocity at the base is 1800m / s, and the density ρ is 1.7g / cm³. 3 The attenuation coefficient is α = 0.517f 1.07 dB / λ.

[0150] A horizontal array of hydrophones was deployed on the seabed, consisting of 60 elements and an aperture of approximately 540m. The element positions are as follows: Figure 3As shown. The target sound source travels at a speed of 6 m / s, with a travel time of 800 s and a distance of 4.8 km. The center frequency of the narrowband signal is 100 Hz, and the bandwidth is 1 Hz. The depth range to be detected is 0–150 m, i.e., the search depth range is 0–150 m, and the phase velocity spectrum scanning range is 1450–1750 m / s.

[0151] Figure 4 A comparison of the estimation results of the present invention and the traditional modal filtering method for sound sources at different depths is presented. The root mean square error of the present method is 7.08m, while that of the traditional modal filtering method is 25.68m, indicating that the present method has higher accuracy in estimating the depth of the sound source.

[0152] Figure 5 The target's bearing trajectory, processed from actual data received by the seabed acoustic array during another sea trial, is presented. The experimental sea area was approximately 180-240m deep. The target sound source was a UW350 launched from the experimental vessel, which drifted freely with the ocean currents. During the experiment, the sound source traveled only a few hundred meters at the same depth; therefore, decoherence was achieved through frequency averaging with a bandwidth of 5Hz.

[0153] Figure 6 A comparison between the estimated depths of different sound sources and the actual depths provided by this invention is given, demonstrating that this invention can estimate sound source depths with relatively high accuracy. These results further verify the feasibility and reliability of this invention.

[0154] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention.

Claims

1. A method for estimating the depth of a moving sound source, characterized in that, include: The acoustic propagation model is used to simulate the sea area to be monitored. The position of the horizontal array of hydrophones is marked in the simulated sea area, and the depth range to be detected in the simulated sea area is divided into multiple depth grids. For any depth grid, a depth is selected from that depth grid as the depth of the target sound source to be determined. Using normal mode theory, it is determined that there is a frequency band of no less than two normal modes in the simulated sea area. From the intersection of this frequency band and the frequency band of the target sound source, an arbitrary frequency is selected as the selected frequency. Using normal mode theory, the phase velocity and eigenfunction of each normal mode at the selected frequency in the simulated sea area are calculated, and the normal mode energy distribution at the depth of each target sound source to be determined is calculated based on the normal mode phase velocity and eigenfunction. The received signals of selected frequencies in each observation period of the hydrophone horizontal array are processed by MVDR to obtain the azimuth of the target source in each observation period. Then, one azimuth is randomly selected from the calculated azimuth of the target source as the reference azimuth. The phase velocity spectrum scanning range is determined by the minimum and maximum values ​​of the phase velocities of each normal mode obtained by simulation calculation. The phase velocity spectrum scanning range is divided into multiple phase velocity grids. A phase velocity is selected from each phase velocity grid. For any selected phase velocity, the following steps are performed: using the phase velocity as the reference velocity, the phase of the received signal in each observation time period is compensated to the reference azimuth. The compensated covariance matrix and the compensated weighting vector are calculated, and the compensated output beam power is calculated to obtain the compensated output beam power curve. For any unknown target sound source depth, the following steps are performed: using the phase velocities of each normal mode calculated by the simulation as the initial values ​​of each normal mode phase velocity, performing an optimal value search for each normal mode phase velocity and updating each normal mode phase velocity as one iteration, and repeating the iteration until the number of iterations is greater than the maximum number of iterations or the error change is less than the convergence threshold, at which point the iteration stops, and the error of the last iteration is recorded as the optimal error for the unknown target sound source depth; The depth of the undetermined target sound source with the smallest optimal error is taken as the target sound source depth. The error of one iteration is the straight-line distance between the compensated output beam power curve and the reconstructed beam power curve. The reconstructed beam power curve is obtained by reconstructing the compensated weighted vector, the normal mode energy distribution, and the phase velocities of each normal mode after the update of this iteration.

2. The method according to claim 1, characterized in that, The energy distribution of each normal mode at the depth of any undetermined target sound source at a selected frequency is calculated using the following formula: in, To represent the depth of the undetermined target sound source The energy of the kth normal mode; The initial phase velocity of the k-th normal mode obtained from simulation calculation; This indicates the depth of the horizontal array of hydrophones. for The eigenfunctions of the k-th normal mode; Select a frequency; right Normalized to: get Energy distribution of normal modes of each order: Where K is the normal mode order.

3. The method according to claim 1, characterized in that, The calculated compensated covariance matrix is ​​expressed by the following formula: ⊙ in, Indicates The covariance matrix after compensation for the reference velocity; This represents the s-th frame of data in the nth observation period; This indicates that the phase velocity obtained from the l-th phase velocity grid is used as the reference velocity; This indicates the location of the target information source during the nth observation period; For reference azimuth; there are a total of N observation time periods, and the received signal in each observation time period is divided into S frames; To select a frequency; by It indicates direction. The unit direction vector; Let m be the position of the m-th array element relative to the reference array element.

4. The method according to claim 1, characterized in that, The formula for the weighted vector after compensation is expressed as follows: in, This indicates that the phase velocity obtained from the l-th phase velocity grid is used as the reference velocity; Indicates The compensated weighted vector calculated based on the reference speed; For reference orientation; m represents the m-th array element; It represents the two-dimensional coordinates of the m-th array element.

5. The method according to claim 1, characterized in that, The process of completing one iteration for any unknown target sound source depth is as follows: The search range for the phase velocity of each normal mode is determined as follows: in, The minimum interval of phase velocity is preset. The maximum range to search; This represents the phase velocity of the current k-th normal mode. Let be the search phase velocity of the current k-th normal mode; We construct the cost functions for each normal mode and calculate the optimal value of the phase velocity of the k-th normal mode, expressed as: in, This represents the optimal value of the phase velocity of the k-th normal mode found through searching. The output beam power curve after compensation; , These are signal energy and noise energy, respectively. Indicates the depth of the undetermined target sound source. The normalized value of the energy of the kth order normal mode; Indicates the depth of the undetermined target sound source. The energy distribution of normal modes other than the k-th normal mode; The beam response curve represents the k-th order search normal mode. This represents the steering vector of the k-th order search normal mode; This represents the beam response curves of the normal modes other than the k-th normal mode. The beam response curve represents the current normal mode; L represents the number of phase velocity grids. Beam response curve representing noise; This represents the steering vector of the current normal mode. For reference direction; This represents the position of the m-th array element relative to the reference array element. Update the current normal wave phase velocity with the optimal values ​​of the phase velocities of each normal wave found, and update the beam response curve of the current normal wave simultaneously; calculate the difference between the current iteration error and the previous iteration error, i.e., the error change, and determine that the error change is not less than the convergence threshold and the number of iterations is not greater than the maximum number of iterations. Use the current iteration error as the previous iteration error and increment the iteration count by 1.

6. The method according to claim 1, characterized in that, For any unknown target sound source depth, the error formula for the first iteration is expressed as: in, This represents the error of that single iteration; This refers to the current k-th normal mode phase velocity after updating the phase velocities of each normal mode in this iteration. The output beam power curve after compensation; , These are signal energy and noise energy, respectively. This represents the beam response curve of the current k-th normal mode after updating the phase velocities of each normal mode in this iteration. This represents the beam response curves of the remaining normal modes, excluding the current k-th normal mode, after the phase velocities of each normal mode are updated in this iteration. Indicates the depth of the undetermined target sound source. The normalized value of the energy of the kth order normal mode; Indicates the depth of the undetermined target sound source. The energy distribution of normal modes other than the k-th normal mode; This represents the beam response curve of the current normal mode after updating the phase velocities of each normal mode in this iteration; Beam response curve representing noise; This represents the steering vector of the current normal mode after the phase velocities of each normal mode are updated in this iteration.

Images