A method for positioning a very low frequency line spectrum target in a continental slope sea area

Abstract

The application discloses a kind of continental slope sea area underwater very low frequency line spectrum target positioning method.The method first constructs the three-dimensional sound field model based on spatial coordinate rotation transformation, solves the problem that the element position of horizontal receiving array cannot accurately fall on the calculation grid point, simultaneously realizes the fast prediction to the element receiving signal of horizontal receiving array.Then the MVDR algorithm is used to estimate the multi-path wave arrival angle of each order mode, and combined with the vertical mode-horizontal ray theory, the underwater sound source position is accurately positioned by horizontal sound line backtracking.Finally, the average mode phase velocity is introduced to improve the accuracy of multi-path wave arrival angle estimation, and the multi-path wave arrival angle of each order mode is more accurately estimated, and the positioning error of 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 a real three-dimensional ocean environment, the uneven topography of the continental slope seafloor deflects the propagation direction of acoustic energy, producing strong three-dimensional effects. The most severe impact on underwater target positioning is horizontal refraction. Under strong horizontal refraction, the azimuth spectrum estimate of the received signal is actually the angle at which the deflected energy reaches the receiving array. Directly interpreting this as the target azimuth can result in errors of up to tens of degrees.

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

[0004] To calculate three-dimensional acoustic fields, it's important to consider the shallow-water waveguides found in most continental slope environments. Under very low frequency conditions, a simple normal wave parabolic equation model was developed that optimally balances accuracy and computational efficiency. However, numerical implementation of this model requires the computational grid to be divided in a rectangular coordinate system. In practical applications, the horizontal receiving array often has an inclination angle with the coordinate axis, preventing the precise placement of each array element on the computational grid. Previously, the most widely used method was to re-mesh and recalculate the acoustic field for each individual array element position, but this approach incurs 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 actual very low frequency conditions in continental slope waters.

[0006] To solve the above technical problems, the present invention aims to achieve the following technical solutions: providing a method for locating underwater very low frequency line spectrum targets in continental slope waters, comprising:

[0007] Obtain the topography of the continental slope sea area and set the sound source parameters under very low frequency conditions;

[0008] Utilizing the sound source parameters to calculate the eigenfunctions and modal function library required for the three-dimensional sound field;

[0009] A three-dimensional sound field model based on the normal wave parabolic equation theory is constructed by using the eigenfunction and the modal function library;

[0010] A spatial coordinate rotation transformation is performed on model parameters of the three-dimensional sound field model to obtain receiving signals of all array elements of a horizontal receiving array;

[0011] Azimuth spectrum estimation is performed by using the receiving signals of all array elements of the horizontal receiving array to obtain a multipath angle of arrival estimation result of each order mode;

[0012] The center position of the horizontal receiving array is taken as an exit position, and the multipath angle of arrival estimation result is taken as an exit angle, and a horizontal sound ray corresponding to each order mode is calculated, and the intersection of the horizontal sound rays is taken as an estimated sound source position;

[0013] The multipath angle of arrival estimation result of each order mode is updated by using the average modal phase velocity of the horizontal sound ray corresponding to each order mode to update the estimated sound source position.

[0014] The embodiment of the present application has the beneficial effects that: a method for positioning an underwater target by using horizontal multipath information under a very low frequency condition is proposed. First, a three-dimensional sound field model based on spatial coordinate rotation transformation is constructed, which solves the problem that the array element positions of the horizontal receiving array cannot be accurately placed on the calculation grid points, and realizes fast prediction of the array element receiving signals of the horizontal receiving array. Then, the MVDR algorithm is used to estimate the multipath angle of arrival of each order mode, and the vertical mode-horizontal ray theory is combined to realize accurate positioning of the underwater sound source position by horizontal sound ray backtracking. Finally, the average modal phase velocity is introduced to improve the accuracy of the multipath angle of arrival estimation, and the multipath angle of arrival of each order mode is more accurately estimated, further reducing the positioning error of the underwater sound source position. BRIEF DESCRIPTION OF DRAWINGS

[0015] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0016] Figure 1 A flowchart of a continental slope sea area underwater very low frequency line spectrum target positioning method provided by the embodiment of the present application is shown in the figure.

[0017] Figure 2 A sub-flowchart of a continental slope sea area underwater very low frequency line spectrum target positioning method provided by the embodiment of the present application is shown in the figure.

[0018] Figure 3A 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;

[0019] 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;

[0020] Figure 5 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;

[0021] 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;

[0022] 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;

[0023] Figure 8 A schematic diagram of the continental slope seabed topography provided by an embodiment of the present invention;

[0024] Figure 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;

[0025] Figure 10 Schematic diagram of the sound field distribution of each mode provided by an embodiment of the present invention;

[0026] Figure 11 A schematic diagram of the multipath angle of arrival estimation results of a horizontal receiving array at a distance of 21 km from the sound source provided by an embodiment of the present invention;

[0027] Figure 12 A schematic diagram of the results of horizontal sound ray back-tracking positioning using the speed of sound in water provided by an embodiment of the present invention;

[0028] Figure 13 A schematic diagram of the horizontal sound ray back-tracking positioning results using the average modal phase velocity provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0029] The following will provide a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort shall fall within the scope of protection of the present invention.

[0030] It will be understood that when used in this specification and the appended claims, the terms “comprises” and “comprising” indicate the presence of described features, integers, steps, operations, elements and / or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or groups thereof.

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

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

[0033] 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.

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

[0035] S101. Obtain the topography of the continental slope sea area and set the sound source parameters under very low frequency conditions;

[0036] S102, using the sound source parameters to calculate the eigenfunction and modal function library required for the three-dimensional sound field;

[0037] S103, using the eigenfunction and modal function library to construct a three-dimensional sound field model based on the simple normal wave parabolic equation theory;

[0038] S104, performing spatial coordinate rotation transformation on the model parameters of the three-dimensional sound field model to obtain reception signals of all array elements of the horizontal receiving array;

[0039] S105, performing azimuth spectrum estimation using the received signals of all elements of the horizontal receiving array to obtain multipath arrival angle estimation results of each order mode;

[0040] S106, using the center position of the horizontal receiving array as the exit position, using the multipath angle of arrival estimation result as the exit angle, calculating the horizontal sound rays corresponding to each mode, and using the intersection of each horizontal sound ray as the estimated sound source position;

[0041] S107 , using the average modal phase velocity of the horizontal sound rays corresponding to each order mode to update the multipath angle of arrival estimation result of each order mode, so as to update the estimated sound source position.

[0042] In this embodiment, a three-dimensional sound field model based on spatial coordinate rotation transformation is first constructed, solving the problem of the horizontal receiving array element position not being able to accurately fall on the calculation grid point, while also achieving rapid prediction of the horizontal receiving array element received signal. The MVDR algorithm is then used to estimate the multipath arrival angle of each mode. Combined with the vertical mode-horizontal ray theory, the horizontal sound ray back-tracing is used to accurately locate the underwater sound source. Finally, the mean modal phase velocity is introduced to improve the accuracy of the multipath arrival angle estimation, re-estimating the multipath arrival angle of each mode order more accurately, further reducing the positioning error of the underwater sound source.

[0043] 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 required to solve the horizontal ray using the Runge-Kutta method; therefore, it can be obtained from the existing ETOP01 data set, which has a data accuracy of 1′. The actual continental slope seabed topography extracted can meet the application requirements. Figure 8 For example, Figure 8 (a) is a topographic plan, Figure 8 (b) is a three-dimensional map of the terrain. In addition, the sound source parameters under the very low frequency condition are set, including: the very low frequency sound source frequency is 25Hz, the number of elements in the horizontal receiving array is 36, the element spacing is 30m corresponding to the half wavelength spacing; the sound source depth is 60m, the receiving depth is 50m; the speed of sound in water is 0. c w =1500m / s, density ρ w =1.0g / cm 3 ; Sound speed at sea bottom c b =1700m / s, density ρ b =1.5g / cm 3 , and the attenuation coefficient α b =0.5dB / λ.

[0044] See also Figure 2 In one embodiment, step S102 includes:

[0045] S201, calculating the eigenfunction of the sound source frequency in the sound source parameters at different sea depths;

[0046] S202, calculating the modal function of the sound source depth in the sound source parameters at different sea depths;

[0047] S203: Calculate the modal function of the receiving depth in the sound source parameters at different sea depths.

[0048] In this embodiment, calculating the three-dimensional sound field requires a library of eigenvalues ​​and modal functions. The seabed topography extracted in step S101 has a depth range of 0-880m, and the depth calculation interval is set to 1m. The CAMBLA method is used to calculate the eigenvalues, modal functions corresponding to the source depth, and modal functions corresponding to the receiving depth (the two modal functions constitute the modal function library) of the sound source frequency at different sea depths. The calculation results are stored in a database for subsequent use in three-dimensional sound field calculations. In this example, the CAMBLA method only takes 63 seconds to calculate the database, while the widely used KrakenC algorithm takes 821 seconds. The efficiency difference between the two methods is even more significant if the seabed topography has a larger depth range and a smaller depth calculation interval. The specific efficiency difference is that the KrakenC algorithm only supports the calculation of eigenvalues ​​and modal functions corresponding to a single seabed depth at a time. If the actual continental slope seabed topography is deep and the depth grid is fine, calling the KrakenC algorithm multiple times will significantly 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 simultaneous configuration of multiple seabed depths and introduces a large number of parallel operations during the program implementation process, which greatly reduces the calculation time of the database.

[0049] See also Figure 3 In one embodiment, step S103 includes:

[0050] S301, setting the position and amplitude of the continuous wave point source and constructing the non-uniform Helmholtz equation;

[0051] S302. Based on the theory of the 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.

[0052] In this embodiment, most continental slope waveguide environments are shallow water waveguides, and the horizontal refraction effect is significant under very low frequency conditions. For non-uniform waveguide environments, a simple normal wave parabolic equation three-dimensional acoustic field model can be established that best balances calculation accuracy and efficiency:

[0053] 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:

[0054] (1);

[0055] 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).

[0056] sound pressure P It can be written as the sum of simple normal waves of each mode:

[0057] (2);

[0058] in, R n ( x,y ) is the amplitude of each order simple normal wave mode, Φ n (z;x,y) are the eigenfunctions of each order simple normal wave, satisfying equations (3) and (4) respectively:

[0059] (3);

[0060] (4);

[0061] 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 should all be set to 0; after ignoring the coupling coefficient, equation (3) becomes the following form:

[0062] (5);

[0063] 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:

[0064] (6);

[0065] R m Indicates the m The modal amplitudes of the first-order simple normal waves, as well as the eigenfunctions and modal function values ​​required for the three-dimensional acoustic field model based on the simple normal wave parabolic equation theory, can all be obtained by reading the database in 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.

[0066] See also Figure 4 In one embodiment, step S104 includes:

[0067] 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 counterclockwise rotation;

[0068] S402, calculating the minimum seabed topography range covering the sound source and the horizontal receiving array;

[0069] S403, calculating the seabed topography range required by the simple normal wave parabolic equation model;

[0070] S404, calculating the seabed topography range required for spatial coordinate rotation transformation;

[0071] S405, performing a spatial coordinate rotation transformation on the calculated seabed topography range so that the position coordinates of all elements of the horizontal receiving array fall on the calculation grid points;

[0072] S406 , performing three-dimensional sound field calculation using the model parameters after the spatial coordinate rotation transformation to obtain reception signals of all array elements of the horizontal receiving array.

[0073] In this embodiment, during the program implementation 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 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 element position, and then perform a sound field calculation for the received signal at each 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 element position cannot fall precisely on the calculation grid points, and at the same time, the calculation of the received signals of all elements can be realized in one model call, which greatly improves the calculation efficiency. The theoretical derivation of this method is as follows:

[0074] Assume that the horizontal receiving array and x The inclination angle between the axes is β , gives the space coordinate rotation transformation matrix for counterclockwise rotation around the z axis T x :

[0075] (7);

[0076] Then determine the minimum seafloor topography that encompasses the sound source and the horizontal receiving array:

[0077] (8);

[0078] 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 the first N The position of the array elements and the first array element, S ( x,y ) is the position of the sound source in the horizontal plane.

[0079] Continue to determine the seabed topography range required by the simple normal wave parabolic equation theory:

[0080] (9);

[0081] in, L rx and L ry Calculate the seabed topography required for the model x Axis direction and y Axis direction range.

[0082] Then continue to determine the seabed terrain range required for spatial coordinate rotation transformation:

[0083] (10);

[0084] in, L X and L Y Calculate the seabed topography required for the model x Axis direction and y Axis direction range.

[0085] After determining the seabed topography range, the ETOP01 dataset 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:

[0086] (11);

[0087] (12);

[0088] (13);

[0089] 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 the rotation transformation, the coordinates of each element of the horizontal receiving array are parallel to x The calculation grid line of the axis is then determined based on 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 be accurately placed on the calculation grid point.

[0090] Finally, the transformed parameters are used as input, and the model is invoked once to perform a three-dimensional sound field calculation to obtain the received signals at all array elements. Compared to traditional methods, the proposed method for rapid three-dimensional sound field calculation based on spatial coordinate rotation transformation significantly improves computational efficiency as the number of receiving array elements increases.

[0091] In a specific example, the horizontal receiving arrays are 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 effects before and after the transformation can be referred to Figure 9 As shown in (a) and (b), Figure 9 (a) shows the seabed topography range before the spatial coordinate rotation transformation. Figure 9 (b) shows the seabed topography range after the spatial coordinates are rotated; the inner dotted box in (a) and (b) of 9 is the minimum seabed topography range covering the sound source and the horizontal receiving array, which is the same as Figure 8 The seabed topography of (a) and (b) corresponds to each other; Figure 9The area demarcated by the outer dotted lines in (a) and (b) is the seabed topography range required by the simple normal wave parabolic equation theory (the model calculation requires the seabed topography to be input in the form of a matrix). All the colored parts in the two figures show the seabed topography range required by the rotation transformation. After the rotation transformation, Figure 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 positional relationship between the sound source position and the position of the horizontal receiving array elements and the array element spacing. It can be ensured that the positions of all array elements 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 as follows: Figure 10 As shown, Figure 10 (a) is the total sound field of the first three modes, Figure 10 (b) is the sound field of the first-order mode, Figure 10 (c) is the sound field of the second-order mode, Figure 10 (d) is the third-order mode sound field. Compared with traditional methods, the present invention calculates the three-dimensional sound field in 483 seconds with 36 array elements, while the traditional method takes 15,464 seconds, improving computational efficiency by more than 32 times. The improvement in computational efficiency is more pronounced with a greater number of array elements.

[0092] See also Figure 5 In one embodiment, step S105 includes:

[0093] S501, using the MVDR algorithm to perform azimuth spectrum estimation using the received signals of all elements of the horizontal receiving array to obtain azimuth spectra of each order mode;

[0094] S502, performing peak scanning on the azimuth spectrum of each mode order to obtain the azimuth spectrum peak value of each mode order;

[0095] S503: Taking the azimuth spectrum peak of each mode as the corresponding multipath angle of arrival estimation result.

[0096] In this embodiment, due to the horizontal refraction effect, the propagation paths of different-order modal energies are different, resulting in different multipath angles of arrival for each mode at the horizontal receiving array. However, in some cases, the arrival energy of individual modes may be extremely low, making it impossible to accurately estimate their multipath angles of arrival. Because the angle of arrival is used as the departure angle in subsequent horizontal sound ray backtracking, the accuracy of its estimation directly affects the final positioning accuracy of the underwater sound source.

[0097] Therefore, the present invention adopts the MVDR high-resolution algorithm, uses the received signal at each element of the horizontal receiving array, and calculates the peak value according to the following azimuth spectrum formula P θCalculation (using the sound speed of 1500m / s in water for azimuth spectrum estimation):

[0098] (14);

[0099] 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 performed, and the 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 21km. The multipath arrival angle estimation result at this location is as follows: Figure 11 As shown, Figure 11 The angles of 1.0° and 5.45° corresponding to the two most obvious peaks are used as the estimated results of the two multipath arrival angles.

[0100] See also Figure 6 In one embodiment, step S106 includes:

[0101] S601. Construct a two-dimensional Helmholtz equation based on the amplitude of each mode.

[0102] S602. Solve the two-dimensional Helmholtz equation using classical ray theory to obtain a partial differential equation;

[0103] S603. Introducing intermediate variables to transform the partial differential equations to obtain a system of ordinary differential equations;

[0104] S604: Using the center position of the horizontal receiving array as the emission position and the multipath angle of arrival estimation result as the emission angle to confirm the initial conditions of the ordinary differential equation system;

[0105] S605. Under the initial conditions, the ordinary differential equations are solved using the ODE45 function in MATLAB to obtain the horizontal sound lines corresponding to each mode;

[0106] S606: Take the intersection of the horizontal sound rays corresponding to each mode as the estimated sound source position.

[0107] In this embodiment, according to the vertical mode-horizontal ray theory, the multipath corresponding to each 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 mode is as follows:

[0108] In the continental slope shallow water waveguide, the amplitude of each mode is A n ( x,y ) with only two horizontal coordinates x and y Related, and satisfy the following two-dimensional Helmholtz equation:

[0109] (15);

[0110] 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, a , For a set small amount, φ ( x,y ) is called the eikonal equation and it satisfies the following partial differential equation (PDE):

[0111] (16);

[0112] By introducing two intermediate variables and , the partial differential equation (PDE) in formula (16) can be transformed into the following set of ordinary differential equations (ODEs):

[0113] (17);

[0114] Wherein, 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 exit 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:

[0115] (18);

[0116] 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 This invention uses the function ODE45 in MATLAB to solve the ordinary differential equations. During the calculation process, the eigenfunction at the recursive position needs to be calculated multiple times. k n , 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 computational efficiency, good computational stability, and convenience for multiple calls in the calculation of intrinsic functions.

[0117] Through the above process, the horizontal sound lines corresponding to each 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 achieve the positioning of the underwater sound source position. Among them, the drawing result of the horizontal sound line is as follows Figure 12As 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%.

[0118] See also Figure 7 In one embodiment, step S107 includes:

[0119] S701, calculating the average modal phase velocity of the horizontal sound rays corresponding to each mode at fixed arc length intervals;

[0120] S702, updating the multipath angle of arrival estimation result of each order mode using the average modal phase velocity of the horizontal sound line corresponding to each order mode;

[0121] S703, using the updated multipath arrival angle estimation results of each order mode to update the horizontal sound rays corresponding to each order mode;

[0122] S704: Update the estimated sound source position using the updated intersection points of the horizontal sound rays corresponding to each mode order.

[0123] In this embodiment, according to simple normal wave theory, the propagation paths and propagation velocities of acoustic energy of various modes are significantly different, so the phase velocities of various modes arriving at the horizontal receiving array are also different. This embodiment proposes to calculate the average modal phase velocity based on the horizontal sound ray results in steps S601-S606, and then use this average modal phase velocity to re-estimate the multipath angle of arrival of each mode with greater accuracy. The derivation of the average modal phase velocity of each order is as follows:

[0124] In steps S601 to S606, the horizontal sound rays corresponding to each mode are known, so the following calculations can be performed at fixed arc length intervals:

[0125] (19);

[0126] 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

[0127] Then, the re-estimated multipath arrival angles of each mode are used as the exit angles for horizontal sound ray back-tracking, ultimately achieving more accurate underwater target positioning.

[0128] In a specific example, the horizontal sound lines in steps S601 to S606 are taken as known conditions, 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 to 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 two average modal phase velocities are used to re-estimate the multipath arrival angles of the first-order mode and the second-order mode. 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 back tracking positioning. The positioning result of the underwater sound source position at this time is as follows: Figure 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.

[0129] For the above-mentioned underwater very low frequency line spectrum target positioning method in continental slope waters, the basic principles and implementation plans 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 average modal phase velocity to perform more accurate multipath wave arrival angle estimation, thereby achieving more precise underwater target positioning, has a significant effect on improving positioning accuracy.

[0130] Those skilled in the art will 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.

[0131] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and such modifications or substitutions are intended to be within the scope of protection of the present invention. Therefore, the scope of protection of the present invention shall be subject to the scope of protection 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 simple normal wave parabolic equation theory is constructed; Performing a spatial coordinate rotation transformation on the model parameters of the three-dimensional sound field model to obtain receiving signals of all elements of the horizontal receiving array; specifically comprising: calculating the inclination angle between the horizontal receiving array and the x-axis, and obtaining a spatial coordinate rotation transformation matrix rotating counterclockwise around the z-axis according to the inclination angle; calculating the minimum seabed terrain range covering the sound source and the horizontal receiving array; calculating the seabed terrain range required by the simple normal wave parabolic equation model; calculating the seabed terrain 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 elements of the horizontal receiving array fall on the calculation grid points; performing a three-dimensional sound field calculation using the model parameters after the spatial coordinate rotation transformation to obtain receiving signals of all 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 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 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 waters 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 using the ETOP01 dataset, and the sound source parameters under very low frequency conditions are set. The sound source parameters include at least the sound source frequency, the number of elements in the horizontal receiving array, the element spacing, the sound source depth, the receiving depth, the water sound speed, the density, the seabed sound speed, and the attenuation coefficient.

3. The method for locating underwater very low frequency line spectrum targets in continental slope waters 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: Calculating 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 waters according to claim 1, characterized in that: The method of constructing a three-dimensional sound field model based on the simple normal wave parabolic equation theory by utilizing the eigenfunction and modal function library 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 waters according to claim 4, characterized in that: The method is based on the theory of the simple normal wave parabolic equation and uses the eigenfunction and modal function library to solve the non-uniform Helmholtz equation, obtain a three-dimensional sound field model and output the sound field of each mode, including: The following formula is used to solve each mode: ; in, 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 waters according to claim 1, characterized in that: The method of 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 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 spectrum of each order mode. Perform peak scanning on the azimuth spectrum of each mode to obtain the azimuth spectrum peak of each mode; The azimuth spectrum peak of each mode is used as the corresponding multipath arrival angle estimation result.

7. The method for locating underwater very low frequency line spectrum targets in continental slope waters according to claim 1, characterized in that: The method of using the center position of the horizontal receiving array as the exit position, using the multipath angle of arrival estimation result as the exit angle, calculating the horizontal sound rays corresponding to each mode, and using the intersection of each horizontal sound ray as the estimated sound source position includes: Based on the amplitude of each mode, a two-dimensional Helmholtz equation is constructed; Solving the two-dimensional Helmholtz equation using classical ray theory yields a partial differential equation; Introducing intermediate variables to transform the partial differential equation to obtain a system of ordinary differential equations; Taking the center position of the horizontal receiving array as the exit position and the multipath angle of arrival estimation result as the exit 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.

8. The method for locating underwater very low frequency line spectrum targets in continental slope waters according to claim 1, characterized in that: The updating of 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 mode are updated using the average modal phase velocity of the horizontal sound line corresponding to each mode. The horizontal sound rays corresponding to each order mode are 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.

9. The underwater VLF line spectrum target positioning method in continental slope waters according to claim 8, wherein the calculation of the average modal phase velocity of the horizontal sound rays corresponding to each mode at fixed arc length intervals 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