A multipath identification and target distance estimation method for deep-sea horizontal line arrays

Abstract

In order to solve the existing application problems of deep-sea horizontal line arrays such as false targets caused by multipath effects and beam splitting and difficulty in estimating target distances, the present invention proposes a multipath identification and target distance estimation method for deep-sea horizontal line arrays. The method uses a two-dimensional cross-correlation matrix of angle-time delay between the output sequence of the expected target main beam direction and the beam sequences in all directions of the whole space. The angle corresponding to the peak position of the matrix can be used as the secondary beam angle value of the target, and the delay difference corresponding to the peak position is the delay difference between the direct path D and the seabed and sea surface primary reflection path B1S1, which can be used to further estimate the distance of the target. The sea trial data verifies the effectiveness of the method; the method is applicable to broadband targets, has low complexity, is simple and effective, is suitable for engineering application implementation, and can be used to support various horizontal arrays in deep-sea application scenarios to carry out the task of underwater acoustic target positioning.

Description

Technical Field

[0001] The invention belongs to the fields of underwater acoustic detection, underwater acoustic positioning, sonar technology, etc., and relates to a multipath identification and target distance estimation method for a deep-sea horizontal line array. Background Art

[0002] The importance of deep-sea maritime safety is becoming increasingly prominent. In order to improve the detection distance of passive sonar, the sonar operating frequency is getting lower and lower, and the array aperture is getting larger and larger. Large-aperture horizontal line arrays can use beamforming technology to obtain high spatial gain, which is suitable for long-range detection and early warning of weak targets, and is convenient for engineering construction deployment. Therefore, it is one of the important equipment for passive underwater acoustic detection. In the deep-sea direct sound area, the sound rays emitted by the near-sea surface sound source usually reach the receiver near the seabed through the direct path D, the sea surface first reflection path S, the seabed first-sea surface first reflection path B1S1, etc. At this time, due to the existence of deep-sea multipath effects, the performance of the linear array beamformer based on the plane wave model will drop sharply and beam splitting will occur, that is, the spatial spectrum of a single target will usually have two beam extremes: one beam extreme corresponds to the arrival angle of the direct path, and its energy is strong, which is called the main beam direction of the desired target; the other beam extreme corresponds to the arrival angle of the B1S1 path, which has a long propagation distance and interacts with the seabed, resulting in large energy loss, so its corresponding beam peak energy is weaker, which is called the secondary beam direction of the desired target. The beam splitting phenomenon will cause the array gain to decrease and reduce the array detection distance on the one hand, and on the other hand, it will cause false target tracks to appear on the array azimuth history diagram, resulting in false alarms.

[0003] At present, many researchers have conducted extensive research on the causes of linear array beam splitting and direction finding deviation (Kunde Yang, Hui Li, Chuanlin He, Rui Duan. Error analysis on bearing estimation of atowed array to a far-field source in deep water. Acoustic Australia, 44, 429-437, 2016). In order to solve the problem of deteriorated direction finding performance, Wu Junnan et al. selected appropriate phase velocity to correct the direction finding deviation of beam forming under the condition of known sound source distance (Wu Junnan et al. Analysis of direction finding of seabed reflection sound of deep sea horizontal array [C]. Proceedings of the 11th Youth Academic Conference of the Chinese Acoustical Society, 2015.), but the direction finding deviation is difficult to correct when the distance of the passive target is unknown. Regarding the problem of false targets caused by beam splitting, in actual scenarios, there are many target trajectories on the azimuth history diagram of the linear array beam forming, and it is difficult to identify which two target trajectories are formed from the same desired target after arriving through multipath. Moreover, the acoustic signal energy after reaching the array through the B1S1 path is weak, and thus the corresponding beam peak energy is weak, which greatly increases the difficulty of identification.

[0004] Generally speaking, only when the spatial position of the desired ship hydroacoustic target is known can ray acoustics be used to accurately predict which tracks on the azimuth history diagram of the linear array are from the azimuth tracks split from the ship after multipath. However, for hydroacoustic targets with unknown spatial positions, how to identify which tracks on the linear array azimuth history diagram are from the same target is one of the important difficulties currently faced by deep-sea passive linear array sonars. The key to solving the above problems faced by deep-sea linear arrays in direction-finding applications lies in solving the problem of target distance estimation based on linear arrays. Summary of the invention

[0005] Technical issues to be solved

[0006] In order to solve the existing application problems such as false targets caused by multipath effects and beam splitting in deep-sea horizontal line arrays and difficulty in estimating target distances, the present invention proposes a multipath identification and target distance estimation method for deep-sea horizontal line arrays. This method uses the output sequence of the expected target main beam direction and the beam sequence in all directions of the whole space to perform an angle-delay two-dimensional cross-correlation matrix. The angle corresponding to the peak position of the matrix can be used as the secondary beam angle value of the target. The delay difference corresponding to the peak position is the delay difference between the direct path D and the first reflection path B1S1 on the seabed, which can be used to further estimate the distance of the target. The sea trial data verified the effectiveness of this method. This method is suitable for broadband targets, has low complexity, is simple and effective, is suitable for engineering application implementation, and can be used to support various horizontal arrays in deep-sea application scenarios to carry out the task of underwater acoustic target positioning.

[0007] Technical Solution

[0008] The present invention proposes a multipath identification and target distance estimation method for deep-sea horizontal line arrays. By cross-correlating the main beam direction output sequence of the desired target with the beam sequences in all directions of the whole space, a two-dimensional cross-correlation matrix of angle-delay is obtained; the relationship between the cross-correlation function and the spatial angle and multipath delay difference can be obtained by integration of different physical dimensions, and the multipath arrival angle and multipath delay difference can be extracted respectively. When the multipath arrival angle coincides with a target trajectory on the array azimuth history diagram, it means that the target comes from a false target; the angle corresponding to the peak position of the trajectory and the multipath delay matrix can be used as the secondary beam angle value of the target, and the delay corresponding to the peak position is the delay difference between the direct path D and the seabed-sea surface reflection path B1S1, which can be used to further estimate the distance of the target.

[0009] The specific steps of this method are as follows:

[0010] S1: Establish a three-dimensional space coordinate system and obtain the sound pressure signal of the linear array:

[0011] Let the position where the first primitive of the linear array is projected onto the sea surface be the origin of the coordinates, the direction of the linear array be the x-axis, and the direction perpendicular to the linear array be the y-axis. Establish a rectangular coordinate system, and the coordinate of the nth primitive be l n =[(n-1)d,0,H] T ,n=1,…,N, where N is the total number of primitives; d is the primitive spacing; H is the known sea depth, which can usually be obtained from actual measurements of the array deployment area or from a nautical chart. The sound source depth is z s The horizontal distance and horizontal azimuth to the head of the linear array are r s and α s , where the horizontal azimuth angle α sIt is defined as the angle between the target projected onto the xoy plane and the array's transverse direction, ranging from -90 degrees to 90 degrees. The horizontal distance between the target and the head of the linear array is usually considered to be no more than 5 times the sea depth, that is, r s ≤5H.

[0012] Assume that the spectrum of the sound pressure signal received by the nth element at time t is p n (f,t), where f represents the signal frequency, f L ≤f≤f H , f L 、f H are the lowest and highest frequencies of the signal frequency band of interest, respectively. Then the sound pressure signals received by N primitives can be expressed in vector form as an N×1-dimensional column vector:

[0013] p(f,t)=[p1(f,t),...,p N (f,t)] T (1)

[0014] The superscript T represents the transpose of a vector.

[0015] S2: Carry out plane wave based beamforming:

[0016] The sound pressure signal p(f, t) received by the linear array is processed by conventional beamforming based on plane waves to obtain a beam domain signal B(θ, f, t) pointing to θ:

[0017] B(θ,f,t)=w H (θ)p(θ,f,t) (2)

[0018] where w(θ)=[1,exp(-ikdsinθ),...,exp[-ik(N-1)dsinθ]] T ; The superscript H represents the conjugate transpose; θ is the guide angle; k represents the wave number, k = ω / c, and c is the speed of sound.

[0019] The azimuth history diagram of the broadband beamforming of the beam domain signal B(θ,f,t) pointing to θ is:

[0020]

[0021] In a deep-sea environment, when the horizontal distance between the receiver on the seabed and the sound source near the sea surface is within 5 times the sea depth (that is, the shallow source depth reception condition is met), the intrinsic sound rays received by the receiver usually come from the three paths of direct path D, sea surface first reflection path S, and seabed first-sea surface first reflection path B1S1. The high-order sound rays that interact with the sea surface and seabed more times are difficult to detect due to large propagation losses (Li Zhenglin, Yang Yixin, Qin Jixing, Yan Wei. Deep Sea Acoustics and Detection Technology [M]. Shanghai Science and Technology Press, 2020). For large-aperture horizontal line arrays, when the target is located near the end-fire direction of the array, the difference in the pitch angles of the two sound rays reaching the array through path D and path S can be ignored, while the pitch angle of path B1S1 is significantly different from the pitch angles of paths D and S. Therefore, the azimuth history diagram E of the array bb (θ,t) will have extreme values ​​in the directions of the two eigenvalues ​​propagating through path D and path B1S1 (Yanqun Wu, et al. Directional response of ahorizontal linear array to an acoustic source at close range in deepwater. Acoustics Australia, 50, 91-103, 2022), that is:

[0022]

[0023] in is the solid angle between the sound line and the linear array, φ l is the elevation angle of the eigenvalues ​​of different paths, l = D, B1S1. The phenomenon described above is the beam splitting phenomenon of the deep-sea line array.

[0024] S3: In the position process diagram E bb (θ, t), manually select the initial angle of the target of interest from the target trajectory in the figure as θ0, and search for the angle θ corresponding to the maximum energy in the angle range θ0±b at this moment s1 (0) and record it as the initial angle of the next moment. b represents the search angle range, which is usually any positive number not greater than 5 degrees. Then, at each subsequent moment, the angle corresponding to the maximum value is recorded as θ according to this method. s1 (1),θ s1 (2),θ s1 (3), ..., until the angle sequence θ of all moments is obtained s1 (t). At this time θ s1 (t) is the main beam direction sequence of the target of interest, and the beam domain signal at this direction is B(θ s1 (t),f,t);

[0025] S4: The beam domain signal B(θ s1 (t),f,t) is cross-correlated with the beam domain signal B(θ,f,t) of any steering angle θ in the whole space. The specific steps are as follows:

[0026] S4.1, B(θ s1 (t),f,t) and B(θ,f,t) are cross-spectral G(θ,f,t), that is,

[0027] G(θ,f,t)=B(θ,f,t)B * (θ s1 (t),f,t) (5)

[0028] S4.2, calculate the generalized cross-correlation function c(θ,τ,t):

[0029]

[0030] Where H(θ,f,t) is a filter, and its selection can be found in the literature "Wei Wenliang, Mao Yulong. Research on a GCC delay estimation method with improved weighting function[J]. Computer Simulation, 2021, 38(5):16-20,84";

[0031] For each time t, the generalized cross-correlation function c(θ,τ,t) is a two-dimensional cross-correlation matrix of angle θ-time delay τ. Obviously, the maximum value of this matrix is ​​located at θ = θ s1 (t) and τ = 0, which is equivalent to the beam domain signal B(θ s1 (t),f,t). In the case of multipath, B(θ s1 (t),f,t) is still strongly correlated with other multipath signals arriving at different angles from the same target, so c(θ,τ,t) will have a maximum value at other angles and time delay positions. At the same time, since the radiation noise signals of different ships are usually uncorrelated, the correlation coefficient of c(θ,τ,t) at the corresponding angle of the interfering ship is very small.

[0032] S5: Obtain the multipath arrival angle trajectory and the multipath delay difference between the direct path and other multipaths:

[0033] Calculate the following integral

[0034]

[0035] where τ min , τ max They correspond to the minimum and maximum values ​​of the delay respectively.

[0036] The function g(θ, t) represents the sum of the mutual correlation coefficients of the signals with different guidance angles θ and the signals arriving at the expected target D path at time t. In order to facilitate the extraction of extreme values, g(θ, t) is normalized using its maximum value. f(τ, t) represents the sum of the mutual correlation coefficients with the signals arriving at the expected target D path in the full spatial orientation at time t with a certain time delay τ. In order to facilitate the subsequent extraction of extreme values, f(τ, t) is normalized using its maximum value. The larger the value of the sum of the mutual correlation coefficients, the more likely it is that the peak at that position comes from the signal of the expected target arriving at the array through multipath. As can be seen from step S4.2, the maximum value of g(θ, t) is located at θ = θ s1 (t), the maximum value of f(τ,t) is at τ=0.

[0037] For the purpose of multipath identification, given a time t, extract the angle sequence θ corresponding to the second maximum value of g(θ,t) B1S1 (t), which is the angle at which the target reaches the array via the B1S1 path. This information can be used for multipath identification. Given a time t, extract the delay corresponding to the second maximum value of f(τ,t) This delay is the propagation time difference between the direct path D and the path B1S1 of the target of interest to reach the array, that is, the propagation delay between the two paths. This information can be used to estimate the target distance.

[0038] In order to achieve the purpose of multipath identification, θ B1S1 The results of (t) are superimposed and plotted on the azimuth history diagram E bb When the trajectory overlaps with the target trajectory in the azimuth history diagram, it means that the trajectory comes from the multipath signal of the expected target. If there is no trajectory overlap, it means that the trajectory of the azimuth history diagram comes from other ship targets.

[0039] S6: Utilizing Delay Calculate the target distance:

[0040] In the direct sound area, the sound line bending can be ignored at close distances, and the slant distances of the sound source along the direct path D and the path B1S1 to the head of the linear array are:

[0041]

[0042] where z r is the array depth.

[0043] For deep seabed arrays, z r =H, and the sound source depth does not exceed 100 meters, so it can be ignored relative to the sea depth, so we have

[0044]

[0045] From this, we can get the theoretical delay value of the sound source reaching the head of the linear array along the path D and the path B1S1:

[0046]

[0047] Where c is the sound velocity at the depth where the linear array is located.

[0048] In order to obtain the above theoretical delay value, the traversal method can be used to convert the sound source distance r s From the minimum value 0 to the maximum value 5H every Δr s Establish the distance grid point, i.e. r sj =jΔr s , where Δr s is the grid resolution, j=1,2,...,J, J=5H / Δr s is the number of grid points; traverse the distance r corresponding to each grid point sj Substituting into formula (11) and (12) we can calculate R D , R B1S1 Then, according to formula (13), the theoretical delay value τ is obtained D-B1S1 (r sj ).

[0049] Construct the following cost function, that is, minimize the measured delay value and the theoretical delay value τ D-B1S1 (r sj ) to get the distance estimate

[0050]

[0051] Cost function C(r sj ) reaches the minimum, the corresponding grid distance point is the estimated value of the first distance rs of the expected target to the array Since formulas (11)-(14) are simple to operate and the amount of calculation is very small, Δr s Can be made very small to ensure accuracy, for example, Δr s =1m.

[0052] Beneficial Effects

[0053] The present invention aims at the problems that are difficult to be solved by the prior art, such as false targets caused by multipath effects and beam splitting after the deep-sea horizontal line array beam forming, and difficulty in estimating the target distance. By utilizing the characteristic that the cross-correlation value of the signal from the same ship target to the array from the direct path and other paths is much greater than the cross-correlation value of the signal from different ship targets, the main beam output sequence of the target of interest and the output sequences in other directions are cross-correlated. The advantage of this technology is that on the one hand, the signal-to-noise ratio of the multipath arrival signal is improved by beam forming, and the azimuth angle of other multipath signals except the direct path is obtained, thereby eliminating false targets on the array azimuth history diagram. On the other hand, the time delay difference between the direct path D and the seabed-sea surface first reflection path B1S1 is obtained, thereby realizing the target distance estimation. This method is suitable for broadband targets, has low complexity, is simple and effective, and is suitable for engineering application implementation. In combination with this method, it can support the deep-sea horizontal line array to carry out the task of target positioning. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 Schematic diagram of the geometric relationship of the sound signal radiated by the moving sound source reaching the deep-sea horizontal linear array via the direct path;

[0055] Figure 2 The sound velocity profile used in the simulation;

[0056] Figure 3 : is a azimuth process diagram of broadband beamforming of the array in the embodiment, in which the dotted line is the delay trajectory θ obtained by extracting the second maximum value of g(θ,t) in step 5 BSR (t).

[0057] Figure 4 In the embodiment, the two-dimensional cross-correlation matrix c(θ, τ, t) and the cross-correlation coefficient g(θ, t) and f(τ, t) at the first moment

[0058] (a) is the two-dimensional cross-correlation matrix of angle θ-time delay τ obtained by cross-correlating the complex beam sequences corresponding to the main beam and the secondary beam;

[0059] (b) is the normalized cross-correlation coefficient and g(θ,t). The three dashed lines in the figure give the theoretical azimuths of the three paths D, B1S1, and B2S2.

[0060] (c) is the normalized cross-correlation coefficient and f(τ, t). The three dashed lines in the figure give the theoretical multipath delay differences between the D path and the D path, between the D path and the B1S1 path, and between the D path and the B1S1 path;

[0061] Figure 5 It is a two-dimensional graph of the mutual correlation coefficient f(τ,t) changing with time;

[0062] Figure 6is a comparison diagram of the target distance estimation result and the actual distance in the embodiment;

[0063] Figure 7 This is the azimuth diagram of the seabed horizontal line array after beamforming in a certain offshore test. The black bold dotted line is the azimuth estimate of the desired target θs1(t), and the pink dotted line is the θ calculated in step 5. B1S1 (t);

[0064] Figure 8 This is the result of the application of the present invention to the offshore test:

[0065] (a) is the multipath delay time series extracted using the extreme value of the cross-correlation matrix;

[0066] (b) is a comparison chart between the target distance estimated using multipath delay and the actual distance. DETAILED DESCRIPTION

[0067] The present invention will be further described below in combination with the embodiments and the accompanying drawings. Through these descriptions, the features and advantages of the present invention will become clearer and more explicit. However, it should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention and are not used to limit the present invention.

[0068] The specific steps of the present invention are as follows:

[0069] S1: Establish a three-dimensional space coordinate system and obtain the sound pressure signal of the linear array:

[0070] Assume that the position where the first primitive of the linear array is projected onto the sea surface is the origin of the coordinate system, the direction of the linear array is the x-axis, and the direction perpendicular to the linear array is the y-axis. Establish a rectangular coordinate system, then the coordinate of the nth primitive is l n =[(n-1)d,0,H] T ,n=1,…,N, where N is the total number of primitives; d is the primitive spacing; H is the known sea depth, which can usually be obtained from actual measurements of the array deployment area or from a nautical chart. The sound source depth is z s The horizontal distance and horizontal azimuth to the head of the linear array are r s and α s , where the horizontal azimuth angle α s It is defined as the angle between the target projected onto the xoy plane and the array's transverse direction, ranging from -90 degrees to 90 degrees. The horizontal distance between the target and the head of the linear array is usually considered to be no more than 5 times the sea depth, that is, r s ≤5H.

[0071] Assume that the spectrum of the sound pressure signal received by the nth element at time t is p n (f,t), where f represents the signal frequency, f L ≤f≤f H , fL 、f H are the lowest and highest frequencies of the signal frequency band of interest, respectively. Then the sound pressure signals received by N primitives can be expressed in vector form as an N×1-dimensional column vector:

[0072] p(f,t)=[p1(f,t),...,p N (f,t)] T (1)

[0073] The superscript T represents the transpose of a vector.

[0074] To illustrate the method of the present invention, in the embodiment, it is assumed that the deep seabed is flat, the sea depth is H = 1245m, the uniform linear array is arranged on the seabed, the array element spacing is d = 6.25m, and the number of primitives is N = 32. The array head is taken as the origin, the array axis is the x-axis, and the direction perpendicular to the array is the y-axis. The sound source transmits a signal frequency band from 40Hz to 200Hz, the depth is 5m, and its depth remains unchanged during the entire movement. The sound source moves in a direction parallel to the x-axis, and the observation time length is T = 20 minutes. Figure 1 A schematic diagram of the geometric relationship between the deep seabed horizontal linear array and the moving sound source is given. In order to obtain the array sound pressure signal p(f, t), the sound pressure field data set from the frequency domain simulation target to different elements of each array is calculated by inputting the underwater acoustic environment parameters into the sound field calculation program KRAKENC, where the simulated underwater acoustic environment parameters: the sound velocity profile is as follows Figure 2 As shown; the seabed parameters are sound speed 1650m / s, density 1.6g / cm 3 , absorption loss is 0.1dB / λ, where λ is the wavelength.

[0075] S2: Carry out plane wave based beamforming:

[0076] The sound pressure signal p(f, t) received by the linear array is processed by conventional beamforming based on plane waves to obtain a beam domain signal B(θ, f, t) pointing to θ:

[0077] B(θ,f,t)=w H (θ)p(θ,f,t) (2)

[0078] where w(θ)=[1,exp(-ikdsinθ),...,exp[-ik(N-1)dsinθ]] T ; The superscript H represents the conjugate transpose; θ is the guide angle; k represents the wave number, k = ω / c, and c is the speed of sound.

[0079] The azimuth history diagram of the broadband beamforming of the beam domain signal B(θ,f,t) pointing to θ is:

[0080]

[0081] In a deep-sea environment, when the horizontal distance between the receiver on the seabed and the sound source near the sea surface is within 5 times the sea depth (that is, the shallow source depth reception condition is met), the intrinsic sound rays received by the receiver usually come from the three paths of direct path D, sea surface first reflection path S, and seabed first-sea surface first reflection path B1S1. The high-order sound rays that interact with the sea surface and seabed more times are difficult to detect due to large propagation losses (Li Zhenglin, Yang Yixin, Qin Jixing, Yan Wei. Deep Sea Acoustics and Detection Technology [M]. Shanghai Science and Technology Press, 2020). For large-aperture horizontal line arrays, when the target is located near the end-fire direction of the array, the difference in the pitch angles of the two sound rays reaching the array through path D and path S can be ignored, while the pitch angle of path B1S1 is significantly different from the pitch angles of paths D and S. Therefore, the azimuth history diagram E of the array bb (θ,t) will have extreme values ​​in the directions of the two eigenvalues ​​propagating through path D and path B1S1 (Yanqun Wu, et al. Directional response of ahorizontal linear array to an acoustic source at close range in deepwater. Acoustics Australia, 50, 91-103, 2022), that is:

[0082]

[0083] in is the solid angle between the sound line and the linear array, φ l is the elevation angle of the eigenvalues ​​of different paths (l = D, B1S1). The phenomenon described above is the beam splitting phenomenon of the deep-sea line array.

[0084] S3: In the position process diagram E bb On (θ, t), given the initial angle of the target of interest as θ0, the angle corresponding to the maximum energy in the search angle range θ0±b at that moment (b represents the search angle range, which can be any positive number not greater than 5 degrees) is θ s1 (0) and record it as the initial angle of the next moment. Then, at each subsequent moment, the angle corresponding to the maximum value is recorded as θ according to this method. s1 (1),θ s1 (2),θ s1 (3), ..., until the angle sequence θ of all moments is obtained s1 (t). At this time θ s1 (t) is the main beam direction sequence of the target of interest, and the beam domain signal at this direction is B(θ s1 (t),f,t);

[0085] In the embodiment, Figure 3 The horizontal array beamforming is performed according to formula (3), and the azimuth course E is obtained. bb (θ,t). Figure 3 It can be seen that after the deep-sea linear array beam is formed, beam splitting occurs, that is, a single target will form two beam tracks, and the main beam points to the direction of arrival of the D path, and the secondary beam points to the direction of arrival of the B1S1 path. When detecting passive targets, it is impossible to distinguish from the azimuth history diagram whether the secondary peak track is from the real target estimation or the false target caused by multipath. Figure 2 In the example, the initial angle θ0 is selected as -60 degrees, and the search angle range b is set to 5 degrees. According to the above method, θ s1 (t).

[0086] S4: The beam domain signal B(θ s1 (t),f,t) is cross-correlated with the beam domain signal B(θ,f,t) of any steering angle θ in the whole space. The specific steps are as follows:

[0087] S4.1, B(θ s1 (t),f,t) and B(θ,f,t) are cross-spectral G(θ,f,t), that is,

[0088] G(θ,f,t)=B(θ,f,t)B * (θ s1 (t),f,t) (5)

[0089] S4.2, calculate the generalized cross-correlation function:

[0090]

[0091] Where H(θ,f,t) is a filter, and its selection can be found in the literature "Wei Wenliang, Mao Yulong. Research on a GCC delay estimation method with improved weighting function[J]. Computer Simulation, 2021, 38(5):16-20,84";

[0092] For each time t, the generalized cross-correlation function c(θ,τ,t) is a two-dimensional cross-correlation matrix of angle θ-time delay τ. Obviously, the maximum value of this matrix is ​​located at θ = θ s1 (t) and τ = 0, which is equivalent to the beam domain signal B(θ s1 (t),f,t). In the case of multipath, B(θ s1(t),f,t) is still strongly correlated with other multipath signals arriving at different angles from the same target, so c(θ,τ,t) will have a maximum value at other angles and time delay positions. At the same time, since the radiation noise signals of different ships are usually uncorrelated, the correlation coefficient of c(θ,τ,t) at the corresponding angle of the interfering ship is very small.

[0093] In the embodiment, The advantage of this selection method is that there is no need to estimate the power spectrum of other physical quantities, so it is more convenient in practical applications. Since when the multipath signals are distributed at different horizontal azimuths, their signals are still correlated with the direct sound signals, peaks will appear at the angles where the different paths are distributed and the multipath delay differences. At the target azimuth, c(θ, τ, t) is a two-dimensional cross-correlation matrix of delay τ-time t, and its peak describes how the multipath delay difference changes with time as the sound source position changes. According to formula (6), the generalized cross-correlation function c(θ, τ, t) changes with angle at the first moment t = 1 as follows: Figure 4 As shown in (a), for each time t, c(θ,τ,t) is a two-dimensional cross-correlation matrix of angle θ-time delay τ. Figure 4 There are three significant peaks in (a). The maximum value is located in the autocorrelation between the direct paths D and D, whose corresponding delay value is 0, and the angle corresponding to D is θ = θ s1 (t). The second maximum is the cross-correlation between the D path and the B1S1 path, and its corresponding delay value is τ 12 (t), and the angle is the angle corresponding to the B1S1 path θ=θ s2 (t). The third weakest maximum value is the correlation coefficient between the direct path D and the signal arriving on the B2S2 path (double reflection on the sea surface - double reflection on the seabed). However, in practical applications, the signal arriving on the B2S2 path is too weak and is usually not considered.

[0094] S5: Obtain the multipath arrival angle trajectory and the multipath delay difference between the direct path and other multipaths:

[0095] Calculate the following integral

[0096]

[0097] where τ min , τ max They correspond to the minimum and maximum values ​​of the delay respectively.

[0098] The function g(θ, t) represents the sum of the mutual correlation coefficients of the signals with different guidance angles θ and the signals arriving at the expected target D path at time t. In order to facilitate the extraction of extreme values, g(θ, t) is normalized using its maximum value. f(τ, t) represents the sum of the mutual correlation coefficients with the signals arriving at the expected target D path in the full spatial orientation at time t with a certain time delay τ. In order to facilitate the subsequent extraction of extreme values, f(τ, t) is normalized using its maximum value. The larger the value of the sum of the mutual correlation coefficients, the more likely it is that the peak at that position comes from the signal of the expected target arriving at the array through multipath. As can be seen from step S4.2, the maximum value of g(θ, t) is located at θ = θ s1 (t), the maximum value of f(τ,t) is at τ=0.

[0099] For the purpose of multipath identification, given a time t, extract the angle sequence θ corresponding to the second maximum value of g(θ,t) B1S1 (t), which is the angle at which the target reaches the array via the B1S1 path. This information can be used for multipath identification. Given a time t, extract the delay corresponding to the second maximum value of f(τ,t) This delay is the propagation time difference between the direct path D and the path B1S1 of the target of interest to reach the array, that is, the propagation delay between the two paths. This information can be used to estimate the target distance.

[0100] In order to achieve the purpose of multipath identification, θ B1S1 The results of (t) are superimposed and plotted on the azimuth history diagram E bb When the trajectory overlaps with the target trajectory in the azimuth history diagram, it means that the trajectory comes from the multipath signal of the expected target. If there is no trajectory overlap, it means that the trajectory of the azimuth history diagram comes from other ship targets.

[0101] In the embodiment, according to formulas (7) and (8), g(θ, t) and f(τ, t) at time t=1 are obtained respectively, see Figure 4 (b) Figure 4 (c), the dotted lines in the figure give the theoretical angles of the two paths D and B1S1 and the multipath delay values. It can be seen that the peak value of g(θ,t) coincides with the theoretical angle value, and the peak value of f(τ,t) coincides with the theoretical delay difference. In order to obtain the extreme values ​​of g(θ,t) and f(τ,t), the findpeaks command in the MATLAB 2019 version is used here to obtain the angle estimate θ of the B1S1 path. B1S1 (t)(see Figure 3 Secondary maximum trajectory) and time-delay trajectory (See Figure 4 ).θ B1S1 (t) is superimposed with a black dashed line. Figure 3 Middle position process map Ebb (θ, t), the dotted line overlaps with the target track with weaker energy in the figure, indicating that the track is not another ship, but a multipath arrival signal from the same target as the main beam signal, thus eliminating the false target.

[0102] S6: Utilizing Delay Calculate the target distance:

[0103] In the direct sound area, the sound line bending can be ignored at close distances, and the slant distances of the sound source along the direct path D and the path B1S1 to the head of the linear array are:

[0104]

[0105] where z r is the array depth.

[0106] For deep seabed arrays, z r =H, and the sound source depth does not exceed 100 meters, so it can be ignored relative to the sea depth, so we have

[0107]

[0108] From this, we can get the theoretical delay value of the sound source reaching the head of the linear array along the path D and the path B1S1:

[0109]

[0110] Where c is the sound velocity at the depth where the linear array is located.

[0111] In order to obtain the above theoretical delay value, the traversal method can be used: the sound source distance r s From the minimum value 0 to the maximum value 5H every Δr s Establish the distance grid point, i.e. r sj =jΔr s , where Δr s is the grid accuracy, j=1,2,...,J, J=5H / Δr s is the number of grid points; traverse the distance r corresponding to each grid point sj Substituting into formula (11) and (12) we can calculate R D , R B1S1 Then, according to formula (13), the theoretical delay value τ is obtained D-B1S1 (r sj ).

[0112] Construct the following cost function, that is, minimize the measured delay value and the theoretical delay value τ D-B1S1 (r sj ) to get the distance estimate

[0113]

[0114] Cost function C(r sj ) reaches the minimum, the corresponding grid distance point is the first distance r of the expected target to reach the array s Estimated value of Since formulas (11)-(14) are simple to operate and the amount of calculation is very small, Δr s Can be made very small to ensure accuracy, for example, Δr s =1m.

[0115] In the embodiment, the distance is gridded from 0m to 8000m, and the grid points Δr s =1m. Applying formula (14), we can get the comparison result between the target distance estimation result and the real distance in the simulation as follows: Figure 6 As shown in the figure, the two correspond well. To further verify the effectiveness of the present invention, the present invention was applied to a marine test, and the array and the marine environment were the same as the previous simulation parameters. The target was the acoustic signal radiated by the cooperative ship, and the effective frequency range was 40Hz to 200Hz. The azimuth history diagram after the array broadband beamforming is shown in the figure below. Figure 7 As shown, the black bold dotted line is the main beam direction θ of the desired target s1 (t), the pink dashed line is θ calculated in step 5 B1S1 (t). Figure 7 It can be seen that in actual scenarios, due to the serious interference of surface ships, there are many tracks on the azimuth history map, which highlights the necessity of multipath track identification. B1S1 The energy trace in the (t) direction is very weak, which is related to the weak energy of the signal arriving through the B1S1 path. Figure 8 is the result of the application of the present invention to the offshore test, Figure 8 (a) is the multipath delay time series extracted using the extreme value of the cross-correlation matrix; Figure 8 (b) is a comparison chart between the target distance estimated by multipath delay and the actual distance. In the experiment, the closest distance between the cooperative ship and the array was 1 km, and the farthest distance was 6 km. Figure 8 (a) shows the comparison between the target distance estimation result and the actual distance. Figure 8 (b) shows the distance estimation error, which is basically within 15%, illustrating the effectiveness of this method.

[0116] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the scope of protection of the present invention.

Claims

1. A multipath identification and target distance estimation method for deep-sea horizontal line array, characterized in that: The specific steps of this method are as follows: S1: Establish a three-dimensional space coordinate system and obtain the sound pressure signal of the linear array: Let the position where the first primitive of the linear array is projected onto the sea surface be the origin of the coordinates, the direction of the linear array be the x-axis, and the direction perpendicular to the linear array be the y-axis. Establish a rectangular coordinate system, and the coordinate of the nth primitive be l n =[(n-1)d,0,H] T ,n=1,…,N, where N is the total number of primitives; d is the primitive spacing; H is the known sea depth; the sound source depth is z s The horizontal distance and horizontal azimuth to the head of the linear array are r s and α s , where the horizontal azimuth angle α s It is defined as the angle between the target projected onto the xoy plane and the array’s transverse direction; the horizontal distance between the target and the head of the linear array does not exceed 5 times the sea depth, i.e. r s ≤5H; Assume that the spectrum of the sound pressure signal received by the nth element at time t is p n (f,t), where f represents the signal frequency, f L ≤f≤f H , f L 、f H are the lowest and highest frequencies of the signal frequency band of interest, respectively. Then the sound pressure signals received by N primitives can be expressed in vector form as an N×1-dimensional column vector: p(f,t)=[p1(f,t),…,p N (f,t)] T (1) The superscript T represents the transpose of the vector; S2: Carry out plane wave based beamforming: The sound pressure signal p(f, t) received by the linear array is processed by conventional beamforming based on plane waves to obtain a beam domain signal B(θ, f, t) pointing to θ: B(θ,f,t)=w H (θ)p(θ,f,t) (2) where w(θ)=[1,exp(-ikdsinθ),...,exp[-ik(N-1)dsinθ]] T ; The superscript H represents the conjugate transpose; θ is the guide angle; i represents the imaginary unit, k represents the wave number, k = ω / c, ω represents the circular frequency, and c is the speed of sound; The azimuth history diagram of the broadband beamforming of the beam domain signal B(θ,f,t) pointing to θ is: In the deep sea environment, when the horizontal distance between the receiver located on the seabed and the sound source near the sea surface is within 5 times the sea depth, the intrinsic sound rays received by the receiver come from the three paths: the direct path D, the sea surface first reflection path S, and the seabed first-sea surface first reflection path B1S1. The high-order sound rays that interact with the sea surface and the seabed more times are difficult to be detected due to the large propagation loss; for the large aperture horizontal line array, when the target is located near the end-fire direction of the array, the difference in the pitch angles of the two sound rays reaching the array through path D and path S can be ignored, while the pitch angle of path B1S1 is quite different from the pitch angles of paths D and S, so the azimuth history diagram E of the array bb (θ, t) will have extreme values ​​in the directions of the two eigenvalues ​​propagating through path D and path B1S1, that is: in is the solid angle between the sound line and the linear array, φ l is the elevation angle of the eigenvalues ​​of different paths, l = D, B1S1; the phenomenon described above is the beam splitting phenomenon of the deep-sea line array; S3: In the position process diagram E bb (θ, t), manually select the initial angle of the target of interest from the target trajectory in the figure as θ0, and search for the angle θ corresponding to the maximum energy in the angle range θ0±b at this moment s1 (0) and record it as the initial angle of the next moment, b represents the search angle range; then at each subsequent moment, the angle corresponding to the maximum value is recorded as θ according to this method s1 (1),θ s1 (2),θ s1 (3), ..., until the angle sequence θ of all moments is obtained s1 (t); at this time θ s1 (t) is the main beam direction sequence of the target of interest, and the beam domain signal at this direction is B(θ s1 (t),f,t); S4: The beam domain signal B(θ s1 (t),f,t) is cross-correlated with the beam domain signal B(θ,f,t) of any steering angle θ in the whole space. The specific steps are as follows: S4.1, B(θ s1 (t),f,t) and B(θ,f,t) are cross-spectral G(θ,f,t), that is, G(θ,f,t)=B(θ,f,t)B * (θ s1 (t),f,t) (5) S4.2, calculate the generalized cross-correlation function c(θ,τ,t): Where H(θ,f,t) is the filter; For each time t, the generalized cross-correlation function c(θ,τ,t) is the two-dimensional cross-correlation matrix of angle θ-time delay τ; obviously, the maximum value of this matrix is ​​located at θ = θ s1 (t) and τ = 0, which is equivalent to the beam domain signal B(θ s1 (t),f,t) autocorrelation function; in the case of multipath, B(θ s1 (t),f,t) is still strongly correlated with other multipath signals arriving from the same target at different angles, so c(θ,τ,t) will have a maximum value at other angles and time delay positions. At the same time, since the radiation noise signals of different ships are usually uncorrelated, the correlation coefficient of c(θ,τ,t) at the corresponding angle of the interfering ship is very small; S5: Obtain the multipath arrival angle trajectory and the multipath delay difference between the direct path and other multipaths: Calculate the following integral where τ min , τ max They correspond to the minimum and maximum values ​​of the delay respectively; Function g(θ,t) represents the sum of the cross-correlation coefficients of the signals with different guiding angles θ and the signals arriving at the desired target path D at time t. In order to facilitate the extraction of extreme values, g(θ,t) is normalized using its maximum value. f(τ,t) represents the sum of the cross-correlation coefficients with the signals arriving at the desired target path D in the full spatial azimuth at time t with a certain time delay τ. In order to facilitate the subsequent extraction of extreme values, f(τ,t) is normalized using its maximum value. The larger the sum of the cross-correlation coefficients, the more likely it is that the peak at that position comes from the signal of the desired target arriving at the array through multipath. As can be seen from step S4.2, the maximum value of g(θ,t) is located at θ=θ s1 (t), the maximum value of f(τ,t) is at τ=0; For the purpose of multipath identification, given a time t, extract the angle sequence θ corresponding to the second maximum value of g(θ,t) B1S1 (t), which is the angle at which the target reaches the array via path B1S1. This information can be used for multipath identification. Given a time t, extract the delay corresponding to the second maximum value of f(τ,t) This delay is the propagation time difference between the direct path D and the path B1S1 of the target of interest to reach the array, that is, the propagation delay between the two paths. This information can be used to estimate the target distance; In order to achieve the purpose of multipath identification, θ B1S1 The results of (t) are superimposed and plotted on the azimuth history diagram E bb (θ, t), when the trajectory overlaps with the target trajectory in the azimuth history graph, it means that the trajectory comes from the multipath signal of the desired target. If there is no trajectory overlap, it means that the trajectory of the azimuth history graph comes from other ship targets. S6: Utilizing Delay Calculate the target distance: In the direct sound area, the sound line bending can be ignored at close distances, and the slant distances of the sound source along the direct path D and the path B1S1 to the head of the linear array are: where z r is the array depth; For deep seabed arrays, z r =H, and the sound source depth does not exceed 100 meters, so it can be ignored relative to the sea depth, so we have From this, we can get the theoretical delay value of the sound source reaching the head of the linear array along the path D and the path B1S1: Where c is the sound velocity at the depth where the linear array is located; In order to obtain the above theoretical delay value, the traversal method is used to convert the sound source distance r s From the minimum value 0 to the maximum value 5H every Δr s Establish the distance grid point, i.e. r sj =jΔr s , where Δr s is the grid resolution, j=1,2,...,J, J=5H / Δr s is the number of grid points; traverse the distance r corresponding to each grid point sj Substituting into formula (11) and (12) we can calculate R D , R B1S1 Then, according to formula (13), the theoretical delay value τ is obtained D-B1S1 (r sj ); Construct the following cost function, that is, minimize the measured delay value and the theoretical delay value τ D-B1S1 (r sj ) to get the distance estimate Cost function C(r sj ) reaches the minimum, the corresponding grid distance point is the first distance r of the expected target to reach the array s Estimated value of 2. A multipath identification and target distance estimation method for deep-sea horizontal line array according to claim 1, characterized in that: In S1, the sea depth H is obtained by actual measurement of the sea area where the array is deployed or from the nautical chart.

3. A multipath identification and target distance estimation method for deep-sea horizontal line array according to claim 1, characterized in that: In S1, the horizontal azimuth angle α s The value range is -90 degrees to 90 degrees.

4. A multipath identification and target distance estimation method for deep-sea horizontal line array according to claim 1, characterized in that: In S3, the search angle range b is any positive number not greater than 5 degrees.

5. A multipath identification and target distance estimation method for deep-sea horizontal line array according to claim 1, characterized in that: In S6, the grid resolution Δr s =1m.

Images