Underwater very-low-frequency line spectrum target positioning method for mainland slope sea area

Abstract

The invention discloses a mainland slope sea area underwater very low frequency line spectrum target positioning method. According to the method, firstly, a three-dimensional sound field model based on space coordinate rotation transformation is constructed, the problem that the array element position of a horizontal receiving array cannot accurately fall on a computational grid point is solved, and meanwhile, rapid forecasting is achieved on an array element receiving signal of the horizontal receiving array. And then estimating the multi-path direction of arrival of each modal by adopting an MVDR algorithm, and realizing accurate positioning of the underwater sound source position through horizontal sound ray anti-tracking in combination with a vertical modal-horizontal ray theory. And finally, the average modal phase velocity is introduced to improve the precision of multi-path direction of arrival estimation, the multi-path direction of arrival of each order of modal is estimated more accurately, and the positioning error of the underwater sound source position is further reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of underwater target acoustic detection, and in particular to a method for locating underwater very low frequency line spectrum targets in continental slope sea areas. Background Art

[0002] In the actual three-dimensional ocean environment, the non-uniform continental slope seabed topography will cause the deflection of the propagation direction of the acoustic energy, resulting in a strong three-dimensional effect. Among them, the horizontal refraction effect has the most serious impact on the positioning of underwater targets. Under the strong horizontal refraction effect, the azimuth spectrum estimation result of the received signal is actually the angle at which the energy reaches the receiving array after deflection. If it is directly used as the target azimuth, it will cause an error of up to tens of degrees.

[0003] Most of the previous studies on the horizontal refraction effect focused on the ASA standard wedge-shaped seabed or directly processed and analyzed experimental data. There are relatively few studies on target positioning simulation for actual three-dimensional seabed terrain and on the mechanism of the influence of horizontal refraction effect on target positioning. The important reasons for this are: the calculation speed of the three-dimensional sound field in the real ocean environment is very slow, and if there is an inclination between the receiving array flow pattern and the calculation coordinate axis, the calculation cost of the array receiving signal will increase several times.

[0004] For the three-dimensional sound field calculation problem, it is necessary to consider that most continental slope environments are shallow water waveguides. Under very low frequency conditions, a simple normal wave parabolic equation model that takes into account both accuracy and computational efficiency is established. However, in the process of implementing this model in the numerical program, it is necessary to divide the computational grid in a rectangular coordinate system. In practical applications, there is usually an inclination between the horizontal receiving array and the coordinate axis, which results in the position of each array element not being able to accurately fall on the computational grid point. Previously, the most widely used method was to re-grid and recalculate the sound field for each single array element position, but this method would incur extremely high computational costs. Summary of the invention

[0005] The purpose of the present invention is to provide a method for locating underwater very low frequency line spectrum targets in continental slope waters, aiming to solve the problem of underwater target positioning under very low frequency conditions in actual continental slope waters.

[0006] In order to solve the above technical problems, the object of the present invention is achieved through the following technical solutions: providing a method for locating underwater very low frequency line spectrum targets in continental slope sea areas, comprising: Obtain the topography of the continental slope sea area and set the sound source parameters under very low frequency conditions; Utilizing the sound source parameters to calculate the eigenfunctions and modal function library required for the three-dimensional sound field; Using the eigenfunction and modal function library, a three-dimensional sound field model based on the normal wave parabolic equation theory is constructed; Performing spatial coordinate rotation transformation on the model parameters of the three-dimensional sound field model to obtain receiving signals of all array elements of the horizontal receiving array; The azimuth spectrum is estimated using the received signals of all elements of the horizontal receiving array to obtain the multipath arrival angle estimation results of each order mode; Taking the center position of the horizontal receiving array as the exit position, taking the multipath angle of arrival estimation result as the exit angle, calculating the horizontal sound rays corresponding to each order mode and taking the intersection of each horizontal sound ray as the estimated sound source position; The multipath angle of arrival estimation result of each order mode is updated by using the average modal phase velocity of the horizontal sound line corresponding to each order mode, so as to update the estimated sound source position.

[0007] The beneficial effects of the embodiments of the present invention are as follows: a method for underwater target positioning using horizontal multipath information under very low frequency conditions is proposed. First, a three-dimensional sound field model based on spatial coordinate rotation transformation is constructed to solve the problem that the position of the elements of the horizontal receiving array cannot accurately fall on the calculation grid points, and at the same time, a rapid prediction of the received signals of the elements of the horizontal receiving array is achieved. Then, the MVDR algorithm is used to estimate the multipath arrival angle of each order mode, and combined with the vertical mode-horizontal ray theory, the horizontal sound line is reversed to achieve accurate positioning of the underwater sound source position. Finally, the average modal phase velocity is introduced to improve the accuracy of the multipath arrival angle estimation, and the multipath arrival angle of each order mode is re-estimated more accurately, further reducing the positioning error of the underwater sound source position. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0009] Figure 1 A schematic diagram of the flow chart of a method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 2 A schematic diagram of a sub-process of a method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 3 A schematic diagram of another sub-process of the method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 4 A schematic diagram of another sub-process of the method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 5A schematic diagram of another sub-process of the method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 6 A schematic diagram of another sub-process of the method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 7 A schematic diagram of another sub-process of the method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention; Figure 8 A schematic diagram of the continental slope seafloor topography provided by an embodiment of the present invention; Fig. 9 A comparison diagram of the seabed topography range before and after the spatial coordinate rotation transformation provided by an embodiment of the present invention; Fig.10 A schematic diagram of the sound field distribution of each mode provided by an embodiment of the present invention; Fig.11 A schematic diagram of the multipath angle of arrival estimation result when the horizontal receiving array is 21 km away from the sound source provided by an embodiment of the present invention; Fig.12 A schematic diagram of a horizontal sound line reverse tracking positioning result using the speed of sound in water provided by an embodiment of the present invention; Fig.13 A schematic diagram of horizontal sound ray back-tracking positioning results using average modal phase velocity provided in an embodiment of the present invention. DETAILED DESCRIPTION

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

[0011] It should be understood that when used in this specification and the appended claims, the terms "include" and "comprises" indicate the presence of described features, integers, steps, operations, elements and / or components, but do not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or combinations thereof.

[0012] It should also be understood that the terms used in this specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. As used in the specification of the present invention and the appended claims, the singular forms "a", "an" and "the" are intended to include plural forms unless the context clearly indicates otherwise.

[0013] It should be further understood that the term "and / or" used in the present description and the appended claims refers to and includes any and all possible combinations of one or more of the associated listed items.

[0014] See also Figure 1 , Figure 1 A schematic flow chart of a method for locating underwater very low frequency line spectrum targets in continental slope waters provided by an embodiment of the present invention.

[0015] like Figure 1 As shown, the method includes steps S101 to S107:

[0016] S101, obtaining the topography of the continental slope sea area and setting the sound source parameters under very low frequency conditions; S102, using sound source parameters to calculate the eigenfunction and modal function library required for the three-dimensional sound field; S103, constructing a three-dimensional sound field model based on the theory of simple normal wave parabolic equation using the eigenfunction and modal function library; S104, performing spatial coordinate rotation transformation on the model parameters of the three-dimensional sound field model to obtain receiving signals of all array elements of the horizontal receiving array; S105, performing azimuth spectrum estimation using the received signals of all array elements of the horizontal receiving array to obtain multipath arrival angle estimation results of each order mode; S106, taking the center position of the horizontal receiving array as the emission position, taking the multipath angle of arrival estimation result as the emission angle, calculating the horizontal sound rays corresponding to each order mode and taking the intersection of each horizontal sound ray as the estimated sound source position; S107, using the average modal phase velocity of the horizontal sound line corresponding to each order mode to update the multipath arrival angle estimation result of each order mode, so as to update the estimated sound source position.

[0017] In this embodiment, a three-dimensional sound field model based on spatial coordinate rotation transformation is first constructed to solve the problem that the position of the elements of the horizontal receiving array cannot accurately fall on the calculation grid points, and at the same time, a rapid prediction of the received signals of the elements of the horizontal receiving array is realized. Then, the MVDR algorithm is used to estimate the multipath arrival angle of each order mode, and combined with the vertical mode-horizontal ray theory, the horizontal sound line is reversed to achieve accurate positioning of the underwater sound source. Finally, the average modal phase velocity is introduced to improve the accuracy of the multipath arrival angle estimation, and the multipath arrival angle of each order mode is re-estimated more accurately, further reducing the positioning error of the underwater sound source.

[0018] In one embodiment, regarding the topography of the continental slope sea area obtained in step S101, the fineness of the topography will affect the accuracy of the three-dimensional sound field calculation and the time for solving the horizontal rays using the Runge-Kutta method; therefore, it can be obtained from the existing ETOP01 data set, the data accuracy of which is 1′, and the actual continental slope seabed topography extracted can meet the application requirements. Figure 8 For example, Figure 8 (a) is the topographic plan, Figure 8 (b) is a topographic stereogram. In addition, the sound source parameters under very low frequency conditions are set, including: the sound source frequency of the very low frequency is 25Hz, the number of elements of the horizontal receiving array is 36, the array element spacing is the corresponding half-wavelength spacing of 30m; the sound source depth is 60m, the receiving depth is 50m; the speed of sound in water is c w =1500m / s, density ρ w =1.0g / cm 3 ; Sound speed at the bottom of the ocean c b =1700m / s, density ρ b =1.5g / cm 3 , and the attenuation coefficient α b =0.5dB / λ.

[0019] See also Figure 2 In one embodiment, step S102 includes: S201, calculating the eigenfunction of the sound source frequency in the sound source parameters at different sea depths; S202, calculating the modal function of the sound source depth in the sound source parameters at different sea depths; S203, calculating the modal function of the receiving depth in the sound source parameters at different sea depths.

[0020] In this embodiment, the calculation of the three-dimensional sound field requires an eigenfunction and a modal function library. The sea depth range of the seabed topography extracted in step S101 is 0-880m, and the sea depth calculation interval is set to 1m. The CAMBLA method is used to calculate the eigenfunction of the sound source frequency at different sea depths in the sound source parameters, the modal function corresponding to the sound source depth, and the modal function corresponding to the receiving depth (the two modal functions constitute the modal function library), and the calculation results are saved as a database so that they can be called at any time during the subsequent three-dimensional sound field calculation. In the example of this embodiment, this embodiment uses the CAMBLA method to calculate the database in only 63s, while the KrakenC algorithm, which is currently widely used, takes 821s, and if the sea depth range of the seabed topography is larger and the sea depth calculation interval is smaller, then the efficiency difference between the two methods of calculating the database will be more obvious. The specific reason for the efficiency difference is that the KrakenC algorithm only supports the calculation of the eigenvalue and modal function corresponding to one seabed depth each time. If the actual continental slope seabed topography selected is deeper and the depth grid division is finer, then calling the KrakenC algorithm multiple times will greatly increase the calculation time. The present invention uses the CAMBLA method to calculate the database. Although the above two models are based on the finite difference method for numerical calculation, the CAMBLA algorithm can support the configuration of multiple seabed depths at the same time, and introduces a large number of parallel operations in the program implementation process, which greatly reduces the calculation time of the database.

[0021] See also Figure 3 In one embodiment, step S103 includes: S301, setting the position and amplitude of the continuous wave point source, and constructing the non-uniform Helmholtz equation; S302. Based on the theory of simple normal wave parabolic equation, the non-uniform Helmholtz equation is solved using the eigenfunction and modal function library to obtain a three-dimensional sound field model and output the sound field of each mode.

[0022] In this embodiment, most of the continental slope waveguide environments are shallow water waveguides, and the horizontal refraction effect is obvious under very low frequency conditions. For non-uniform waveguide environments, a simple normal wave parabolic equation three-dimensional acoustic field model that best balances calculation accuracy and calculation efficiency can be established: Set the CW point source at ( x 0, y 0, z 0), the amplitude is S ( ω ), then the sound pressure P The non-uniform Helmholtz equation satisfied is: (1); in, k is the wave number, k=ω / c ( x,y,z ), ωis the angular frequency, ω = 2 π f , f is the sound source frequency, c is the speed of sound, ρ is the density, δ represents the Dirac function (unit impulse function).

[0023] Sound pressure P It can be written as the sum of simple normal waves of each order mode: (2); in, R n ( x,y ) is the amplitude of each order of simple normal wave mode, Φ n (z;x,y) are the eigenfunctions of simple normal waves of various orders, satisfying equations (3) and (4) respectively: (3); (4); The continental slope ocean environment satisfies the adiabatic approximation condition, so the coupling coefficient in the above formula is A mn , B mn , a mn and b mn All should be set to 0; after ignoring the coupling coefficient, equation (3) becomes the following form: (5); in, k m It is m The horizontal wave number of the order mode. Through the above derivation, the three-dimensional non-uniform Helmholtz equation problem in equation (1) is successfully transformed into the two-dimensional Helmholtz equation problem in equation (5), which can be solved by the classical parabolic equation theory. The solution of each order mode amplitude can be obtained by the following formula: (6); R m Indicates m The modal amplitude of the first order simple normal wave, the eigenfunction and modal function values ​​required in the three-dimensional acoustic field model based on the simple normal wave parabolic equation theory can all be obtained by reading the database of steps S201 to S203. Based on the above derivation, this embodiment completes the three-dimensional acoustic field calculation in the continental slope waveguide environment under very low frequency conditions.

[0024] See also Figure 4 In one embodiment, step S104 includes: S401, calculate the horizontal receiving array and x The inclination angle between the axes is obtained according to the inclination angle. z The space coordinate rotation transformation matrix of the axis rotating counterclockwise; S402, calculating the minimum seabed terrain range covering the sound source and the horizontal receiving array; S403, calculating the seafloor topography range required by the simple normal wave parabolic equation model; S404, calculating the seabed terrain range required for the spatial coordinate rotation transformation; S405, performing a spatial coordinate rotation transformation on the calculated seabed terrain range so that the position coordinates of all array elements of the horizontal receiving array fall on the calculation grid points; S406 , performing three-dimensional sound field calculation using the model parameters after the spatial coordinate rotation transformation, to obtain receiving signals of all array elements of the horizontal receiving array.

[0025] In this embodiment, in the process of implementing the program of the three-dimensional sound field model based on the simple normal wave parabolic equation, it is necessary to divide the calculation grid in the Cartesian coordinate system. However, in actual application, there is usually an inclination angle between the horizontal receiving array and the coordinate axis, which causes the position of each array element of the horizontal receiving array to not fall precisely on the calculation grid. Before the present invention, the most common way was to divide the calculation grid once for each array element position, and then perform a sound field calculation for the received signal at each array element. This method greatly increases the calculation cost. To this end, the present invention proposes a three-dimensional sound field fast calculation method based on spatial coordinate rotation transformation. The principle of this method is simple but extremely practical: it solves the problem that the coordinates of the array element position cannot fall precisely on the calculation grid points, and at the same time, the calculation of the received signals of all array elements can be realized in one model call, which greatly improves the calculation efficiency. The theoretical derivation of this method is as follows: Assume that the horizontal receiving array and x The inclination angle between the axes is β , giving the space coordinate rotation transformation matrix for counterclockwise rotation around the z axis T x : (7); Then determine the minimum seafloor terrain range that covers the sound source and the horizontal receiving array: (8); in, L x and L y are the smallest seafloor topography mentioned above. x Axis direction and y Axis range, E represents the position coordinates of the array element,E N ( x,y )and E 1 ( x,y ) are respectively N The position of the array elements and the first array element, S ( x,y ) is the location of the sound source in the horizontal plane.

[0026] Continue to determine the seafloor topography range required for the simple normal wave parabolic equation theory: (9); in, L rx and L ry The seafloor topography required for model calculation is x Axis direction and y Axis direction range.

[0027] Then continue to determine the seabed terrain range required for the spatial coordinate rotation transformation: (10); in, L X and L Y The seafloor topography required for model calculation is x Axis direction and y Axis direction range.

[0028] After determining the seabed topography range, the ETOP01 data set in step S101 is used to extract the seabed topography of the corresponding area. Finally, the sound source position, the position of each element of the horizontal receiving array, and the extracted seabed topography matrix are subjected to spatial coordinate rotation transformation: (11); (12); (13); in, S ( x,y,z ) represents the spatial position coordinates of the sound source, E n ( x,y,z ) represents the spatial position coordinates of each horizontal receiving array element, B ( x,y,z ) represents the seafloor topography matrix extracted from the ETOP01 dataset. The subscripts in (11)-(13) r Represents the coordinates after rotation transformation. After rotation transformation, the coordinates of each element of the horizontal receiving array are parallel to xThe calculation grid line of the axis is then defined according to the distance between the sound source and the horizontal receiving array and the spacing between the array elements, so that the position coordinates of each array element can fall accurately on the calculation grid point.

[0029] Finally, the parameters after rotation transformation are used as input, and the model is called once to perform three-dimensional sound field calculation to obtain the received signals at all array elements. Compared with the traditional method, the three-dimensional sound field fast calculation method based on spatial coordinate rotation transformation proposed by the present invention has a more obvious improvement effect on the calculation efficiency if the number of receiving array elements is larger.

[0030] In a specific example, the horizontal receiving array is x For example, if there is an inclination angle of 30° between the coordinate axes, the spatial coordinate rotation transformation is performed according to the above equations (7) to (13). The effect before and after the transformation can be referred to Fig. 9 As shown in (a) and (b), Fig. 9 (a) shows the seafloor topography range before the spatial coordinate rotation transformation. Fig. 9 (b) shows the seafloor topography range after the spatial coordinates are rotated and transformed; the inner dotted box in (a) and (b) of 9 is the minimum seafloor topography range covering the sound source and the horizontal receiving array, which is consistent with Figure 8 The seafloor topography of (a) and (b) corresponds; Fig. 9 The area demarcated by the outer dotted boxes in (a) and (b) is the seafloor topography range required by the simple normal wave parabolic equation theory (the model calculation requires the seafloor topography to be input in the form of a matrix), and all the colored parts in the two figures show the seafloor topography range required for the rotation transformation. Fig. 9 The seabed topography delineated by the outer dotted box in (b) is extracted, and then the sound source position and the position of each element of the horizontal receiving array are also rotated and transformed. At this time, the calculation grid is delineated according to the position relationship between the sound source position and the horizontal receiving array elements and the array element spacing, which can ensure that all array element positions can fall on the calculation grid points, and the receiving signals at all array elements can be obtained in one sound field calculation. The calculation results of the sound field are shown in Fig.10 As shown, Fig.10 (a) is the total sound field of the first three modes, Fig.10 (b) is the sound field of the first-order mode. Fig.10 (c) is the sound field of the second-order mode, Fig.10 (d) is the sound field of the third-order mode. By comparing with the traditional method, the three-dimensional sound field calculation time of the present invention is 483s when the number of array elements is 36, while the calculation time of the traditional method is 15464s, which can improve the calculation efficiency by more than 32 times. The more array elements there are, the more obvious the improvement effect of the calculation efficiency is.

[0031] See also Figure 5 In one embodiment, step S105 includes: S501, using the MVDR algorithm to perform azimuth spectrum estimation using the received signals of all array elements of the horizontal receiving array to obtain azimuth spectrograms of each order mode; S502, performing peak scanning on the azimuth spectrum of each order mode to obtain the azimuth spectrum peak value of each order mode; S503: Taking the azimuth spectrum peak of each order mode as the corresponding multipath angle of arrival estimation result.

[0032] In this embodiment, under the horizontal refraction effect, the propagation paths of different order modal energies are different, resulting in different multipath arrival angles of each order mode arriving at the horizontal receiving array. However, in some cases, the arrival energy of individual modes is extremely low, and it is impossible to accurately estimate their multipath arrival angles. Because the arrival angle will be used as the exit angle in the subsequent horizontal sound line back-tracking positioning, its estimation accuracy directly affects the positioning accuracy of the final underwater sound source position.

[0033] Therefore, the present invention adopts the MVDR high-resolution algorithm, uses the received signals at each element of the horizontal receiving array, and calculates the peak value according to the following azimuth spectrum formula: P θ Calculation (the azimuth spectrum is estimated using the sound speed in water of 1500m / s): (14); Where a(θ) is the array steering vector of the horizontal receiving array, R is the signal autocorrelation matrix, and the peak scanning of the above azimuth spectrum is obtained. P θ The peak value is the estimated multipath arrival angle. Specifically, the distance between the horizontal receiving array and the sound source is set to 21 km. The multipath arrival angle estimation result at this location is as follows: Fig.11 As shown, Fig.11 The angles of 1.0° and 5.45° corresponding to the two most obvious peaks are used as the estimation results of the two multipath arrival angles.

[0034] See also Figure 6 In one embodiment, step S106 includes: S601, constructing a two-dimensional Helmholtz equation based on the amplitude of each mode; S602. Solve the two-dimensional Helmholtz equation by classical ray theory to obtain a partial differential equation; S603, introducing intermediate variables to transform the partial differential equation to obtain a system of ordinary differential equations; S604, taking the center position of the horizontal receiving array as the emission position, and taking the multipath arrival angle estimation result as the emission angle, so as to confirm the initial conditions of the ordinary differential equation group; S605, under the initial conditions, solving the ordinary differential equations by using the ODE45 function in MATLAB to obtain the horizontal sound lines corresponding to each mode; S606: Taking the intersection of the horizontal sound rays corresponding to each order mode as the estimated sound source position.

[0035] In this embodiment, according to the vertical mode-horizontal ray theory, the multipath corresponding to each order mode can be drawn as a horizontal sound line from the sound source position to the receiving point. The derivation process of the horizontal sound line of each order mode is as follows: In the shallow water waveguide on the continental slope, the amplitude of each mode is A n ( x,y ) only with two horizontal coordinates x and y Related, and satisfies the following two-dimensional Helmholtz equation: (15); in, k n For the n The eigenfunction of the order mode, Equation (15) can be solved by classical ray theory. Based on ray theory, we can introduce a , For a set small amount, φ ( x,y ) is called the eikonal equation and it satisfies the following partial differential equation (PDE): (16); By introducing two intermediate variables and , the partial differential equation (PDE) of formula (16) can be transformed into the following set of ordinary differential equations (ODEs): (17); Here, the subscript s represents the derivative of the arc length. At this time, if the center position of the horizontal receiving array is set as the origin, that is, the emission position of the horizontal sound line, and the multipath arrival angle estimated in steps S501 to S503 is θ As the exit angle, the initial conditions of the ODEs are as follows: (18); The ODEs of formula (17) can be solved by the Runge-Kutta method. The two columns of data in the solution are x and y The present invention uses the function ODE45 in MATLAB to solve the ordinary differential equation system. During the calculation process, the eigenfunction at the recursive position needs to be calculated multiple times. kn , the present invention uses the NM-CT model to solve it. Compared with the KrankenC and CAMBLA models mentioned above, although this model cannot calculate the modal function, it has the advantages of high calculation efficiency, good calculation stability, and convenience for multiple calls in the calculation of intrinsic functions.

[0036] Through the above process, the horizontal sound lines corresponding to each order mode emitted from the receiving horizontal array can be obtained, and then the intersection of each sound line is used as the estimated sound source position to locate the underwater sound source position. Among them, the drawing result of the horizontal sound line is as follows Fig.12 As shown in the figure, the red solid line represents the horizontal sound line corresponding to the first-order mode, the blue solid line represents the horizontal sound line corresponding to the second-order mode, the symbol "O" represents the actual sound source position, and "×" represents the intersection of the horizontal sound lines of each order mode, that is, the sound source position estimated by this embodiment. It can be seen that the error between the estimated sound source position and the actual position is small. The positioning error of this method in this case is 4.14%.

[0037] See also Figure 7 In one embodiment, step S107 includes: S701, calculating the average modal phase velocity of the horizontal sound rays corresponding to each order mode at fixed arc length intervals; S702, updating the multipath angle of arrival estimation result of each order mode by using the average modal phase velocity of the horizontal sound line corresponding to each order mode; S703, updating the horizontal sound rays corresponding to each order mode using the updated multipath angle of arrival estimation results of each order mode; S704: Update the estimated sound source position using the updated intersection points of the horizontal sound rays corresponding to each order mode.

[0038] In this embodiment, according to the simple normal wave theory, the propagation path and propagation speed of each order modal acoustic energy are very different, so the phase velocity of each order mode arriving at the horizontal receiving array is also different. This embodiment proposes to solve the average modal phase velocity based on the horizontal sound line results in steps S601~S606, and then use the average modal phase velocity to re-estimate the multipath arrival angle of each order mode more accurately. The derivation of the average modal phase velocity of each order is as follows: In steps S601 to S606, the horizontal sound lines corresponding to each order mode are known, so the following calculations can be performed at fixed arc length intervals: (19); in, c mean represents the average modal phase velocity, s total represents the total arc length of the horizontal sound line, sm and c m Respectively represent m Segment arc length and its corresponding modal phase velocity Then, the re-estimated multipath arrival angles of each order mode are used as the exit angles for horizontal sound line back-tracking, ultimately achieving more accurate underwater target positioning.

[0039] In a specific example, the horizontal sound line in steps S601-S606 is taken as a known condition, and the arc length interval is fixed to 50m. The average modal phase velocity of the horizontal sound lines of the first-order mode and the second-order mode in steps S601-S606 along their respective propagation paths is calculated by formula (19). After calculation, the average modal phase velocity of the first-order mode is 1509m / s, and the average modal phase velocity of the second-order mode is 1540m / s. The multipath arrival angles of the first-order mode and the second-order mode are re-estimated using the two average modal phase velocities. At this time, the multipath arrival angle of the first-order mode is updated from 1.0° to 1.01°, and the multipath arrival angle of the second-order mode is updated from 5.45° to 5.60°. Then, the updated multipath arrival angle is used as the exit angle for horizontal sound line reverse tracking positioning. The positioning result of the underwater sound source position at this time is as follows: Fig.13 As shown in the figure, it can be seen that the positioning result using the average modal phase velocity is better than the positioning result using the water sound speed. The positioning error is reduced from the aforementioned 4.14% to 1.91%, and the positioning accuracy is significantly improved.

[0040] For the above-mentioned underwater very low frequency line spectrum target positioning method in the continental slope sea area, the basic principle and implementation scheme of the present invention have been verified by computer numerical simulation, and the results show that: the simple normal wave parabolic equation model of spatial coordinate rotation transformation proposed in the present invention can accurately predict the three-dimensional sound field, and the calculation efficiency is significantly improved; the method for underwater target positioning based on horizontal sound line back tracking proposed in the present invention has high positioning accuracy; the method for using the average modal phase velocity to perform more accurate multipath angle of arrival estimation, thereby achieving more precise underwater target positioning, has a significant effect on improving positioning accuracy.

[0041] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working processes of the above-described equipment, devices and units can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0042] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of various equivalent modifications or replacements within the technical scope disclosed by the present invention, and these modifications or replacements should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A method for locating underwater very low frequency line spectrum targets in continental slope waters, characterized in that: include: Obtain the topography of the continental slope sea area and set the sound source parameters under very low frequency conditions; Utilizing the sound source parameters to calculate the eigenfunctions and modal function library required for the three-dimensional sound field; Using the eigenfunction and modal function library, a three-dimensional sound field model based on the normal wave parabolic equation theory is constructed; Performing spatial coordinate rotation transformation on the model parameters of the three-dimensional sound field model to obtain receiving signals of all array elements of the horizontal receiving array; The azimuth spectrum is estimated using the received signals of all elements of the horizontal receiving array to obtain the multipath arrival angle estimation results of each order mode; Taking the center position of the horizontal receiving array as the exit position, taking the multipath angle of arrival estimation result as the exit angle, calculating the horizontal sound rays corresponding to each order mode and taking the intersection of each horizontal sound ray as the estimated sound source position; The multipath angle of arrival estimation result of each order mode is updated by using the average modal phase velocity of the horizontal sound line corresponding to each order mode, so as to update the estimated sound source position.

2. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The obtaining of the topography of the continental slope sea area and setting of the sound source parameters under very low frequency conditions include: The topography of the continental slope sea area is obtained through the ETOP01 data set, and the sound source parameters under very low frequency conditions are set. The sound source parameters at least include the sound source frequency, the number of elements of the horizontal receiving array, the array element spacing, the sound source depth, the receiving depth, the water sound speed, density, the seabed sound speed, and the attenuation coefficient.

3. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The eigenfunction and modal function library required for calculating the three-dimensional sound field using the sound source parameters includes: Calculate the eigenfunction of the sound source frequency in the sound source parameters at different sea depths; Calculating the modal function of the sound source depth in the sound source parameters at different sea depths; The modal function of the receiving depth in the sound source parameters at different sea depths is calculated.

4. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The method of using the eigenfunction and modal function library to construct a three-dimensional sound field model based on the normal wave parabolic equation theory includes: Set the position and amplitude of the continuous wave point source and construct the inhomogeneous Helmholtz equation; Based on the theory of simple normal wave parabolic equation, the non-uniform Helmholtz equation is solved using the eigenfunction and modal function library to obtain a three-dimensional sound field model and output the sound field of each mode.

5. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 4, characterized in that: Based on the theory of simple normal wave parabolic equation, the eigenfunction and modal function library are used to solve the non-uniform Helmholtz equation, obtain a three-dimensional sound field model and output the sound field of each order mode, including: The following formula is used to solve each mode: ； in, k is the wave number, k = ω / c (x, y, z) , ω is the angular frequency, R m For the m The simple normal wave modal amplitude of the order mode.

6. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The performing of spatial coordinate rotation transformation on the model parameters of the three-dimensional sound field model to obtain receiving signals of all array elements of the horizontal receiving array includes: Calculate the horizontal receiving array and x The inclination angle between the axes is obtained according to the inclination angle. z The space coordinate rotation transformation matrix of the axis rotating counterclockwise; Calculate the minimum seafloor topography range that covers the sound source and horizontal receiving array; Calculate the seafloor topography range required for the simple normal wave parabolic equation model; Calculating the seabed topography range required for the spatial coordinate rotation transformation; Performing a spatial coordinate rotation transformation on the calculated seabed terrain range so that the position coordinates of all array elements of the horizontal receiving array fall on the calculation grid points; The three-dimensional sound field is calculated by using the model parameters after the spatial coordinate rotation transformation to obtain the receiving signals of all array elements of the horizontal receiving array.

7. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The method of using the received signals of all array elements of the horizontal receiving array to perform azimuth spectrum estimation to obtain multipath arrival angle estimation results of each order mode includes: The MVDR algorithm is used to estimate the azimuth spectrum using the received signals of all elements of the horizontal receiving array to obtain the azimuth spectra of each order mode. Perform peak scanning on the azimuth spectrum of each mode to obtain the azimuth spectrum peak value of each mode; The azimuth spectrum peak of each mode is taken as the corresponding multipath angle of arrival estimation result.

8. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The method uses the center position of the horizontal receiving array as the emission position, uses the multipath angle of arrival estimation result as the emission angle, calculates the horizontal sound rays corresponding to each order mode, and uses the intersection point of each horizontal sound ray as the estimated sound source position, including: Based on the amplitude of each mode, the two-dimensional Helmholtz equation is constructed; Solving the two-dimensional Helmholtz equation by classical ray theory to obtain a partial differential equation; Introducing intermediate variables to transform the partial differential equation to obtain a system of ordinary differential equations; The center position of the horizontal receiving array is used as the emission position, and the multipath arrival angle estimation result is used as the emission angle to confirm the initial conditions of the ordinary differential equation group; Under the initial conditions, the ordinary differential equations are solved by using the ODE45 function in MATLAB to obtain the horizontal sound lines corresponding to each mode; The points of the horizontal sound line corresponding to each mode are used as the estimated sound source positions.

9. The method for locating underwater very low frequency line spectrum targets in continental slope sea areas according to claim 1, characterized in that: The method of updating the multipath angle of arrival estimation result of each order mode by using the average modal phase velocity of the horizontal sound line corresponding to each order mode to update the estimated sound source position includes: The average modal phase velocity of the horizontal sound lines corresponding to each mode is calculated at fixed arc length intervals; The multipath arrival angle estimation results of each order mode are updated using the average modal phase velocity of the horizontal sound line corresponding to each order mode; The horizontal sound line corresponding to each order mode is updated using the updated multipath arrival angle estimation results of each order mode; The estimated sound source position is updated using the intersection points of the horizontal sound rays corresponding to the updated modes.

10. The underwater very low frequency line spectrum target positioning method in the continental slope sea area according to claim 9, wherein the average modal phase velocity calculation of the horizontal sound lines corresponding to each order mode at a fixed arc length interval comprises: The average modal phase velocity is calculated as follows: ； in, c mean represents the average modal phase velocity, s total represents the total arc length of the horizontal sound line, s m and c m Respectively represent m Segment arc length and its corresponding modal phase velocity.

Images