A sound ray tracing incident angle calculation method based on a weight function

Abstract

The application discloses a sound ray tracking incidence angle calculation method based on a weight function, relates to the technical field of ocean underwater acoustic positioning, and is used for obtaining a sound ray tracking incidence angle in precise positioning of a seabed control point realized by relying on a sound ray tracking algorithm, and comprises the following steps: solving an approximate incidence angle based on coordinates of a shipborne transducer and a seabed transponder; converting an initial incidence angle by a Snell refraction law to obtain an iteration initial value; calculating the initial incidence angle by a weight function to obtain a Newton down-slope factor; constructing a sound velocity profile layer by using a depth gauge and a sound velocity profiler to obtain a time delay calculation formula; constructing a Newton down-slope method iteration formula; and calculating the sound ray tracking incidence angle based on the obtained approximate incidence angle, the iteration initial value, the Newton down-slope factor and the time delay calculation formula. The application improves the traditional Newton iteration method, improves the stability and iteration efficiency of the Newton iteration method for solving the sound ray tracking incidence angle, and solves the problem that the sound ray tracking cannot be performed when the incidence angle is unknown.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of ocean underwater acoustic positioning technology, mainly to the field of sound ray tracing algorithm optimization of underwater acoustic positioning system, and particularly to a sound ray tracing incident angle calculation method based on a weight function. BACKGROUND

[0002] The seabed geodetic network can provide navigation guarantee for various unmanned devices on the water surface and underwater, and can also be used to monitor the dynamic changes of seabed plates and water environments, and is an important infrastructure for marine safety, marine economic development and marine environment monitoring. The GNSS-A (Global Navigation Satellite System-Acoustics) observation technology, also known as global navigation satellite / acoustic positioning technology, is a combination of dynamic positioning and underwater acoustic ranging technology. The position of the acoustic transducer is transmitted to the seabed control point by using the underwater acoustic ranging principle, which is a key means for constructing and maintaining the seabed geodetic network. The accurate solution of the position of the seabed control point requires accurate measurement of the seawater sound speed.

[0003] However, the seawater sound speed is affected by the temperature, salinity and pressure of seawater, and is a function of time and space, having time-varying and space-varying characteristics. Treating the seawater sound speed as a fixed value will introduce system errors in the ranging results. In order to eliminate the influence of the sound speed error, the spatial symmetry and observation synchronism of the transponder can be used to eliminate the influence of the sound speed variation; the single difference method and double difference method are used to eliminate the influence of common errors such as long-period term system errors in the case of sufficient observation data; the error compensation algorithm is used to eliminate the influence of the time-varying error and horizontal heterogeneity of the sound speed profile, and the sound ray tracing algorithm is used to directly correct the slant range. The implementation of these error elimination methods all depends on the sound ray tracing algorithm.

[0004] In the process of positioning the seabed control point by relying on the sound ray tracing algorithm, the propagation direction of the sound ray is determined by the sound speed at the current position and the incident angle, and the geometric constraint of the sound ray propagation has a strong dependence on the incident angle. If the incident angle of the sound ray tracing is unknown, the sound ray tracing cannot be performed. The root solving algorithms for numerical analysis in existing sound ray tracing algorithms mainly focus on the bisection method, chord intersection method, Newton iteration method and their improved algorithms. The main role is to gradually approach the real path of sound ray propagation or related parameters (such as incident angle, propagation time, horizontal distance, etc.) through numerical iteration, so as to solve the problem of solving nonlinear equations in sound ray tracing. Among them, the bisection method divides the search interval into two parts by constantly dividing the search interval into two parts, determines the sub-interval where the root is located according to the function value sign change, and gradually narrows the range until the accuracy requirement is met. The bisection method is not sensitive to the initial interval selection and is suitable for handling coarse positioning problems in complex acoustic models, but it is not suitable for accurate positioning of seabed control points. The chord intersection method replaces the derivative with the slope of the line connecting two points, and approximates the root through the iteration formula, which only requires function value calculation and is suitable for scenarios where the derivative in the acoustic model is difficult to express analytically, but it requires two initial values, and if the function changes sharply near the iteration point, it may cause convergence failure and poor robustness. Newton iteration method uses the tangent line of the function at the iteration point to approximate the root, has fast convergence speed and controllable accuracy, and can meet the high-precision requirements of acoustic positioning by adjusting the number of iterations or tolerance, and has wide applicability. But the Newton iteration method is sensitive to the initial value, and the incident angle needs to be selected correctly.

[0005] Among the many incident angle solving algorithms, there is less research on the failure of Newton iteration method at large incident angles. Newton's down-hill method can reduce the influence of initial value on solving stability by reducing the step size, but the traditional Newton's down-hill method usually needs to find the down-hill factor by trial (such as halving method), which increases the calculation amount and reduces the convergence efficiency; while directly selecting the initial down-hill factor according to experience lacks theoretical basis and adaptability.

[0006] Therefore, there is a need for a sound ray tracing incident angle calculation method that can adaptively calculate the Newton down-hill factor, can calculate the sound ray tracing incident angle, and can improve the iteration efficiency as much as possible under the premise of ensuring the success of the incident angle iteration, which is used to optimize the sound ray tracing algorithm in the underwater acoustic positioning system. SUMMARY

[0007] The present application aims at the divergence problem of the existing traditional Newton iteration method at large incident angle and the problem that the sound ray tracing cannot be performed when the incident angle of the sound ray tracing is unknown, and provides a sound ray tracing incident angle calculation method based on a weight function, the sound ray tracing incident angle is obtained by constructing a Newton down-hill method iteration formula, and performing calculation based on a Snell constant iteration initial value, a Newton down-hill factor and a time delay calculation formula of the Snell constant, the Newton iteration method is changed in iteration step by introducing the Newton down-hill factor, the requirement for the initial value of the Newton iteration method in iteration at large incident angle is reduced, the iteration efficiency is improved while ensuring the iteration success.

[0008] A sound ray tracing incident angle calculation method based on a weight function, comprising:

[0009] S1. Obtain the approximate incident angle based on the obtained coordinates by obtaining the approximate coordinates of the seabed transponder and the coordinates of the shipborne transducer at the target time;

[0010] S2. Obtain the Snell constant iteration initial value by converting the approximate incident angle through the Snell refraction law;

[0011] S3. Construct or select a weight function, convert the approximate incident angle through the weight function to obtain a Newton down-hill factor;

[0012] S4. Establish a layer time delay calculation formula of the Snell constant based on the sound ray tracing algorithm;

[0013] S5. Obtain the depth from the shipborne transducer to the seabed transponder, and calculate the sound speed profile layer through which the sound ray passes from the shipborne transducer to the seabed transponder;

[0014] S6. Establish a time delay calculation formula of the Snell constant based on the sound speed profile layer and the layer time delay calculation formula of the Snell constant;

[0015] S7. Construct a Newton down-hill method iteration formula, and calculate based on the Snell constant iteration initial value, the Newton down-hill factor and the time delay calculation formula of the Snell constant to obtain the sound ray tracing incident angle.

[0016] S1 comprises:

[0017] S1.1 Obtain the coordinates of the ship antenna based on the global navigation satellite system (GNSS);

[0018] The coordinates of the ship antenna in the Earth-Centered Earth-Fixed (ECEF) coordinate system are obtained by using the global navigation satellite system (GNSS) receiver loaded on the ship The obtained coordinates of the ship antenna in the Earth-Centered Earth-Fixed (ECEF) coordinate system are expressed as:

[0019] ;

[0020] Wherein, denotes the time of coordinate acquisition, 、 、 denote the projection distances of the ship body antenna in the X-axis, Y-axis, and Z-axis directions in the Earth-Centered Earth-Fixed coordinate system, respectively;

[0021] converts into coordinates in the East-North-Up ENU coordinate system :

[0022] ;

[0023] wherein the coordinates of the phase center of the ship body antenna in the ECEF coordinate system are selected as the coordinates of the reference point, is the reference point longitude, is the reference point latitude, denotes the component of the ship body antenna in the east direction, denotes the component of the ship body antenna in the north direction, denotes the component of the ship body antenna in the vertical direction. S1 further comprises:

[0024] S1.2 calculates the position of the ship-borne transducer in the local ENU coordinate system at the target time

[0025] based on the coordinates of the ship body antenna in the ENU coordinate system in combination with the ship body attitude and the position of the transducer relative to the GNSS antenna in the gyroscope rectangular coordinate system. :

[0026] ;

[0027] ;

[0028] ;

[0029] ;

[0030] ;

[0031] ;

[0032] wherein, is a time-varying rotation matrix, denotes the ship body attitude, 、 、 denote the roll, pitch, and heading attitudes of the ship body, respectively, denotes the position of the ship-borne transducer relative to the GNSS antenna in the gyroscope rectangular coordinate system, ,​ , respectively represent the projection components on the X-axis, Y-axis and Z-axis in the corresponding rectangular coordinate system, , represent the intermediate calculation matrix, represent the matrix transpose.

[0033] S1 further comprises:

[0034] S1.3 obtains the approximate coordinates of the seabed transponder based on the average sound speed method according to the principle of distance intersection;

[0035] Set the approximate coordinates of the seabed transponder as , select observation data of observation epochs, and based on the observation data of each observation epoch, use the least square method to obtain the approximate coordinates of the seabed transponder, and the specific expression is as follows:

[0036] ;

[0037] wherein, is the average sound speed calculated according to the sound speed profile, is the time delay corresponding to each observation epoch, , , represent the coordinates of the observation points of each observation epoch.

[0038] S1 further comprises:

[0039] S1.4 calculates the approximate incidence angle based on the approximate coordinates of the seabed transponder and the position of the ship-borne transducer in the local ENU coordinate system at the time The approximate incidence angle is:

[0040] ;

[0041] wherein, represents the approximate incidence angle.

[0042] An iterative formula related to the Snell constant is established, the approximate incidence angle is converted according to the Snell refraction law, and the initial value of the Snell constant iteration is obtained:

[0043] ;

[0044] wherein, is the sound speed of the incident layer;

[0045] A weight function is constructed or selected, the obtained approximate incidence angle is brought into the weight function, and the Newton descent factor is solved.

[0046] Based on the formula for calculating intra-layer delay related to the incident angle in the ray tracing algorithm, the formula for calculating intra-layer delay related to the Snell constant is established as follows:

[0047] ;

[0048] in, Indicates the first Intra-layer delay of the Snell constant of the layer, For the first The sound velocity gradient of the layer; and These are the sound speed profiles. Layers and The speed of sound in the layer, Snell's constant;

[0049] The depth between the shipborne transducer and the seabed transponder is obtained by a depth gauge, and the sound velocity profile of the waters where the shipborne transducer and the seabed transponder are located is obtained by a sound velocity profiler. Based on the obtained depth and sound velocity profile from the shipborne transducer to the seabed transponder, a sound velocity profile layer is constructed that the sound ray passes through from the shipborne transducer to the seabed transponder.

[0050] Based on the sound velocity profile layers that the sound rays pass through from the shipborne transducer to the seabed transponder and the formula for calculating the intra-layer time delay of the Snell constant, a formula for calculating the time delay of the Snell constant is established:

[0051] ;

[0052] Converting the time delay calculation formula into a form that solves for the zeros of a function is equivalent to obtaining the functional expression of the time delay calculation formula and then calculating the first derivative of that functional expression:

[0053] ;

[0054] ;

[0055] in, This indicates the number of sound velocity profile layers that the sound ray passes through from the shipboard transducer to the seabed transponder. Here is the functional expression for the time delay calculation formula. for The first derivative.

[0056] An adaptive Newton's downhill method iterative formula is constructed by substituting the obtained initial value of the Snell constant, the Newton's downhill factor, the function expression of the time delay calculation formula, and the first derivative of the function expression of the time delay calculation formula into the adaptive Newton's downhill method iterative formula. The adaptive Newton's downhill method iterative formula is as follows:

[0057] ;

[0058] At this time, , is the iteration output value;

[0059] iterative operation is carried out based on the adaptive Newton downhill method iterative formula to obtain the iteration final value ;

[0060] If the iteration output value meets the iteration termination condition, then ; otherwise, let and repeat the iteration until the iteration termination condition is met.

[0061] When the selected initial Newton down factor causes iteration failure, at this time , the iteration is restarted, and if it still fails, then is halved until the iteration is successful.

[0062] The Snell refraction law is used to convert into , and the sound ray tracing incidence angle is obtained.

[0063] Compared with the prior art, the present application has the following beneficial effects:

[0064] The present application provides a sound ray tracing incidence angle calculation method based on a weight function, constructs a Newton down method iterative formula, constructs an adaptive Newton down method iterative formula based on a Snell constant iteration initial value, an adaptive Newton down factor and a time delay calculation formula of the Snell constant, carries out iterative calculation based on the adaptive Newton down method iterative formula, obtains the sound ray tracing incidence angle, and solves the problem that if the incidence angle of the sound ray is unknown, the sound ray tracing cannot be carried out in the process of realizing the positioning of the seabed control point by relying on the sound ray tracing algorithm.

[0065] The present application provides a sound ray tracing incidence angle calculation method based on a weight function, improves the traditional Newton iteration method, introduces the weight into the Newton down factor on the basis of the traditional Newton iteration method by analogy with the weight in the adjustment random model, realizes the adaptive adjustment of the size of the Newton down factor, reduces the iteration step length when the incidence angle is large, and improves the iteration efficiency on the premise of ensuring the iteration success of the incidence angle. BRIEF DESCRIPTION OF DRAWINGS

[0066] Figure 1 A flowchart of a sound ray tracing incidence angle calculation method based on a weight function.

[0067] Figure 2 A sailing path schematic diagram.

[0068] Figure 3 is a schematic diagram of a sound speed profile.

[0069] Figure 4 is a schematic diagram of a weight function.

[0070] Figure 5 is a schematic diagram of a ray tracing algorithm.

[0071] Figure 6 is a schematic diagram of residual variation of different weight models at the 332nd observation epoch.

[0072] Figure 7 is a schematic diagram of iteration number variation of different weight models. DETAILED DESCRIPTION

[0073] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions in the present application are described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

[0074] As shown in Figure 1 , a weight function-based ray tracing incidence angle calculation method comprises the following steps:

[0075] S1. obtaining the approximate coordinates of a seabed transponder and the coordinates of a shipborne transducer at a target time; S2. converting the approximate incidence angle by Snell's law of refraction to obtain an initial value of Snell's constant; S3. constructing or selecting a weight function, and converting the approximate incidence angle by the weight function to obtain a Newton's descent factor;

[0076] S4. establishing a layer internal time delay calculation formula of Snell's constant based on a ray tracing algorithm;

[0077] S5. obtaining the depth from the shipborne transducer to the seabed transponder, and calculating the sound speed profile layer through which the sound ray passes from the shipborne transducer to the seabed transponder;

[0078] S6. establishing a time delay calculation formula of Snell's constant based on the sound speed profile layer and the layer internal time delay calculation formula of Snell's constant;

[0079]

[0080]

[0081] ​​​​​S7. Constructing the iteration formula of Newton's down-hill method, calculating based on the iteration initial value of Snell constant, the calculation formula of Newton's down-hill factor and the time delay of Snell constant, obtaining the sound ray tracing incident angle.

[0082] The approximate coordinates of the seabed transponder represent the static position of the seabed transponder in the earth coordinate system, which is the reference point of acoustic positioning. As a fixed reference point of the acoustic positioning system, it is used to determine the relative position of the ship-borne equipment or target.

[0083] The ship-borne transducer coordinates are the three-dimensional geographic coordinates of the acoustic transducer (device for transmitting / receiving sound waves) installed on the ship body at a certain time, which dynamically changes with the movement of the ship body and reflects the spatial position of the transducer at the measurement moment. Real-time calculation of ship-borne transducer coordinates is realized by acquiring the coordinates of the ship body antenna through GNSS, combining the ship body attitude sensor (IMU) and the transducer installation offset, and real-time solving the transducer coordinates.

[0084] S1 includes:

[0085] S1.1 Acquiring the coordinates of the ship body antenna based on the global navigation satellite system GNSS;

[0086] Firstly, the coordinates of the ship body antenna in the earth-centered earth-fixed coordinate system ECEF are acquired in real time by the global navigation satellite system GNSS receiver loaded on the ship body , the acquired coordinates of the ship body antenna in the earth-centered earth-fixed coordinate system are expressed as:

[0087] ;

[0088] Wherein, represents the time of coordinate acquisition, , , respectively represent the projection distances of the ship body antenna in the X-axis, Y-axis and Z-axis directions in the earth-centered earth-fixed coordinate system;

[0089] Then, the coordinates of the ship body antenna in the earth-centered earth-fixed coordinate system are converted into the coordinates in the east-north-up ENU coordinate system :

[0090] ;

[0091] Wherein, the coordinates of the phase center of the ship body antenna in the ECEF coordinate system are selected as the coordinates of the reference point, is the reference point longitude, is the reference point latitude, represents the east component of the ship body antenna, represents the north component of the ship body antenna, represents the vertical component of the ship body antenna.​

[0092] The GNSS receiver can obtain the position of the surveying ship in real time, and the positions of the acoustic transducers can be calculated based on the coordinate correction by determining the relative position relationship between the ship bottom transducer and the GNSS receiver.

[0093] The ECEF coordinate system takes the center of the earth as the origin, and the coordinate value of the ship body antenna is usually millions to tens of millions of meters. Direct use of these large values for calculation is easy to cause numerical overflow or precision loss. The ENU coordinate system is a local coordinate system, taking the ship body antenna as the origin, and converting the coordinates into local values in meters. By converting the coordinates in the ECEF coordinate system into the coordinates in the ENU coordinate system, the large value operation is converted into a small value operation, which can avoid complex trigonometric function operation in the ECEF coordinate system and significantly improve the calculation stability and precision.

[0094] S1 further comprises:

[0095] S1.2 calculating the coordinates of the ship body antenna in the ENU coordinate system in combination with the ship attitude and the position of the transducer in the gyroscope rectangular coordinate system relative to the GNSS antenna, to calculate the position of the ship-borne transducer in the local ENU coordinate system at the target time

[0096] ;

[0097] ;

[0098] ;

[0099] ;

[0100] ;

[0101] ;

[0102] wherein, is a time-varying rotation matrix, represents the ship attitude, 、 、 represent the roll, pitch and heading attitudes of the ship body, respectively, represents the position of the ship-borne transducer in the gyroscope rectangular coordinate system relative to the GNSS antenna, 、 、 represent the projection components of in the X-axis, Y-axis and Z-axis of the corresponding rectangular coordinate system, respectively, 、 represent the intermediate calculation matrix.​​

[0103] The position of the transducer in the local ENU coordinate system is calculated according to the above method As shown in Figure 2 , the GNSS receiver provides real-time surface position information of the surveying ship, and the acoustic sounding device obtains the relative distance information between the ship and the seafloor through interaction with the seafloor transponder. The GNSS receiver can obtain the position information of the surveying ship in real time, and the motion sensor can synchronously measure the attitude information of the ship body. In the ship body coordinate system, based on the position relationship between the GNSS antenna and the transducer and the motion attitude sensing data of the ship body, the position of the acoustic transducer can be determined. Figure 2 The dashed line represents the sailing path of the surveying ship, which is equivalent to the position of the ship-borne transducer in the local ENU coordinate system at different times, that is, the motion path of the ship-borne transducer in the horizontal plane. The pentagram represents the position of the M01 transponder. As an underwater acoustic device, the transponder is usually installed at a fixed position on the seafloor. When the transponder receives an acoustic signal from the ship, the transponder will transmit a response signal to the ship. By calculating the round-trip time of the acoustic signal, the sound speed, and other information, the distance between the ship and the transponder can be calculated, and the position of the ship relative to the transponder can be determined.

[0104] S1 further comprises:

[0105] S1.3 calculates the approximate coordinates of the seafloor transponder based on the distance intersection principle through the average sound speed method;

[0106] The approximate coordinates of the seafloor transponder are set as , and the observation data of observation epochs are selected. Based on the observation data of each observation epoch, the least squares method is used to calculate the approximate coordinates of the seafloor transponder, and the specific expression is as follows:

[0107] ;

[0108] wherein is the average sound speed obtained by averaging the sound speed profile according to the sound speed profile diagram as shown in Figure 3 , is the time delay corresponding to each observation epoch, , , represent the coordinates of each observation epoch.

[0109] Figure 3 is a sound speed profile diagram of the sea area from 0 to 4000 meters, showing the variation of sound speed with water depth, as shown in Figure 3 , the horizontal axis represents the sound speed, with units of meters per second (m / s), and the vertical axis represents the depth, with units of meters (m). From Figure 3It can be seen that the sound speed does not change linearly with depth, but shows a clear layered feature, at a certain depth There is a minimum value (sound channel axis), the sea surface sound speed is greater than the sea bottom sound speed, which is a full deep sea sound channel. Above the sound channel axis, the sound speed increases with the increase of seawater depth, and the sound speed changes rapidly. Below the sound channel axis, the sound speed decreases with the increase of seawater depth, but in the deep sea area (after the seawater depth is greater than 1500 meters), the sound speed tends to be stable with the change of seawater, which is equivalent to the slope of the sound speed profile gradually tends to be constant.

[0110] The observation epoch refers to the time point or time period when the satellite signal at a certain specific time is observed in the GNSS measurement process.

[0111] The average sound speed method may refer to considering the change of sound speed in different water layers, and taking the average value to calculate the distance. In addition, the specific implementation of the average sound speed method may need to consider the sound speed profile, such as dividing the water layer into multiple layers, each layer having different sound speed, and then taking the average.

[0112] S1 also includes:

[0113] S1.4 the approximate coordinates based on the bottom transponder and The position of the ship-borne transducer in the local ENU coordinate system at the moment Calculate the approximate incidence angle, and the approximate incidence angle is:

[0114] ;

[0115] Wherein, The approximate incidence angle is denoted as.

[0116] An iterative formula related to the Snell constant is established, and the approximate incidence angle is converted according to the Snell refraction law to obtain the initial value of the Snell constant iteration :

[0117] ;

[0118] Wherein, The sound speed of the incident layer is denoted as.

[0119] A weight function is constructed or selected, and the obtained approximate incidence angle is brought into the weight function to solve the Newton descent factor .

[0120] The Newton descent factor is an adaptive Newton descent factor, and the Newton descent factor has various weight function expression forms as shown in Table 1:

[0121] Table 1 weight function expression form of Newton descent factor

[0122] .

[0123] in, , , , These represent the coefficients of the corresponding weighting functions. , , , These represent the constant terms of the corresponding weight functions, and the function curves for each weight function are shown below. Figure 4 As shown.

[0124] like Figure 4 As shown, the horizontal axis represents the incident angle in degrees (°), and the vertical axis represents the weight value. When the Newton descent factor is 1, it is equivalent to the traditional Newton-Raphson iteration method, which has a second-order convergence rate. Therefore, 1 is taken as the weight value. Figure 4 The upper limit of the weight value; the downhill factor in the first halving trial is 0.5, that is, the weight value is 0.5. In order to select a suitable downhill factor for direct iteration before performing the halving trial, 0.5 is used as... Figure 4 The lower limit of the weighted value. As the incident angle increases, the difference between the approximate incident angle and the true incident angle increases significantly, making... Compared to The growth rate is faster, leading to and As the ratio increases, the iteration step size increases. To ensure the success rate of iteration, a smaller weight value should be selected, that is, the larger the incident angle, the smaller the weight value.

[0125] This invention addresses the issue of divergence in the traditional Newton-Raphson iteration method at large incident angles. It introduces an adaptive Newton downhill factor to reduce the step size and improve iteration stability. While the selection of a fixed Newton downhill factor is random and subjective, halving the trial calculations wastes computation. This invention introduces an incident angle weighting function, selecting different downhill factors based on the magnitude of the incident angle. This ensures iteration stability while minimizing halving the trial calculations and improving convergence efficiency, enabling rapid and stable acquisition of the incident angle for ray tracking.

[0126] Figure 5 This is a schematic diagram of a voice tracking algorithm, such as... Figure 5 As shown, the speed of sound changes with a constant gradient in each layer (left figure). According to the theory of ray acoustics, the trajectory of the sound ray in each layer is an arc (right figure). The depth at which the sound wave is emitted. Let be the speed of sound at the depth where the sound wave is emitted. Let the be... The depths at the top and bottom interfaces of the layer are respectively and The speeds of sound are respectively and The incident angle within the layer is The angle of departure is The layer thickness is The sound velocity gradient within the layer is The horizontal propagation distance of the sound wave within that layer is The radius of curvature of the circular arc is .

[0127] Based on the intra-layer delay calculation formula related to the incident angle in the constant gradient ray tracing algorithm, the intra-layer delay calculation formula related to the Snell constant is established as follows:

[0128] ;

[0129] in, Indicates the first Intra-layer delay of the layer's Nell constant, For the first The sound velocity gradient of the layer; and These are the sound speed profiles. Layers and The speed of sound in the layer, Snell's constant;

[0130] The Snell constant is a key parameter in acoustics and optics used to describe the path characteristics of waves propagating in different or layered media. Its definition and application vary depending on the characteristics of the medium. The Snell constant p is a parameter that remains unchanged when a wave propagates in a layered medium, and its expression is: ,in, Indicates a depth of The speed of sound at that location Indicates a depth of The angle of incidence at that point.

[0131] The depths of the shipborne transducer and the seabed transponder are obtained by a depth gauge, and the sound velocity profiler is used to obtain the sound velocity profile of the water area. Based on the obtained depth data and sound velocity profile, a sound velocity profile layer is constructed to show the path of the sound ray from the shipborne transducer to the seabed transponder.

[0132] Based on the sound velocity profile layers that sound rays pass through from the shipborne transducer to the seabed transponder and the intra-layer time delay calculation formula of the Snell constant, a single-pass time delay calculation formula for the Snell constant is established:

[0133] ;

[0134] Converting the one-way delay calculation formula into a form that solves for the zeros of a function is equivalent to obtaining the functional expression of the one-way delay calculation formula and then calculating the first derivative of that functional expression:

[0135] ;

[0136] ;

[0137] wherein, represents the number of layers of the sound speed profile layer through which the sound ray passes from the ship-mounted transducer to the seabed transponder, is a form of the function zero point of the one-way time delay calculation formula, is a first derivative of the form of the function zero point of the one-way time delay calculation formula.

[0138] The adaptive Newton downhill method iterative formula is constructed, the obtained Snell constant iterative initial value, the Newton downhill factor, the function expression of the one-way time delay calculation formula and the first derivative of the function expression of the one-way time delay calculation formula are substituted into the adaptive Newton down hill method iterative formula, and iterative operation is carried out based on the adaptive Newton down hill method iterative formula to obtain an iterative final value :

[0139] ;

[0140] At this time, , is an iterative output value;

[0141] If the iterative output value meets the iterative termination condition, then ; otherwise, let and repeat the iteration until the iterative termination condition is met. The iterative termination condition is that the vertical displacement of the sound ray is the same as the depth value from the ship-mounted transducer to the seabed transponder, and the iteration is stopped.

[0142] When the selected initial Newton downhill factor leads to iteration failure, at this time let , and start the iteration again, if it still fails, then halve until the iteration is successful; convert to by using Snell's law to obtain the sound ray tracking incidence angle .

[0143] In a real sea trial scene, the incidence angles calculated by using different weight functions are approximately the same, and the approximate incidence angle is equivalent to the incidence angle obtained by the average sound speed method, which is often smaller than the real incidence angle. The larger the incidence angle is, the more obvious the difference between the incidence angle calculated by the average sound speed method and the incidence angle calculated by the weight function is, which indicates that the Snell constant iterative initial value obtained based on the approximate incidence angle is farther away from the real value. Therefore, when solving the incidence angle, the situation of too large incidence angle must be considered.

[0144] The calculation residual of the model solved by different weight functions is significantly reduced relative to the average sound velocity method, the unit weight error of the sound ray tracing model established by different incidence angle solving methods is reduced by about 0.7 meters, and the purpose of sound ray correction is achieved, so the adaptive Newton down-hill method using different incidence angle weight functions has higher precision in solving the incidence angle relative to the average sound velocity method. The change of the calculation residual with the observation epoch is related to the size of the incidence angle solved at each observation epoch, and the greater the incidence angle, the greater the absolute value of the residual.

[0145] Figure 6 The residual change of different weight models at the 332th observation epoch is shown in FIG. 1, wherein the curve of the cosine function model is at the top, the absolute value of the residual is the smallest, the solved incidence angle is the most accurate, the cosine function is the most consistent with the change of the iteration step in the selected function, and has the optimal effect, so the Newton down-hill method established according to the cosine weight function has the optimal effect in solving the incidence angle. Figure 6

[0146] The iteration convergence speed of different models is evaluated, and the iteration convergence times of different models are shown in FIG. 2. Figure 7 For the same function model (same color line), the iteration times change with the observation epoch and the size of the incidence angle solved at each epoch, the greater the incidence angle, the more the iteration convergence times, and the slower the convergence. For different function models of the same observation epoch (different color lines), the corresponding function curve of the cosine function model is at the bottom of FIG. 2, so the iteration times of the cosine function model are the least, and the iteration convergence is the fastest, so the Newton down-hill method established according to the cosine weight function has the fastest iteration convergence speed in solving the incidence angle, and the iteration times are the least in most epochs. Figure 7

[0147] The above is only the preferred embodiment of the present application, and it should be noted that those skilled in the art can make several improvements and refinements without departing from the principles of the present application, and these improvements and refinements should also be considered as the protection scope of the present application.​​

Claims

1. A method for computing an angle of incidence of a sound ray based on a weight function, characterized in that, Comprise: S1. Obtain the approximate coordinates of the seabed transponder and the coordinates of the shipborne transducer at the target time, and calculate the approximate incidence angle based on the obtained coordinates; S2. Convert the approximate incidence angle by Snell's law of refraction to obtain the initial value of the Snell constant iteration; S3. Construct or select a weight function, convert the approximate incidence angle by the weight function to obtain the Newton's descent factor; S4. Establish the intralayer time delay calculation formula of the Snell constant based on the sound ray tracing algorithm; S5. Obtain the depth from the shipborne transducer to the seabed transponder, and calculate the sound speed profile layer through which the sound ray passes from the shipborne transducer to the seabed transponder; S6. Establish the time delay calculation formula of the Snell constant based on the sound speed profile layer and the intralayer time delay calculation formula of the Snell constant; S7. Construct the Newton's descent iteration formula, and calculate based on the initial value of the Snell constant iteration, the Newton's descent factor and the time delay calculation formula of the Snell constant to obtain the sound ray tracing incidence angle.

2. The method of claim 1, wherein, S1 comprises: S1.1 Obtain the coordinates of the ship antenna based on the global navigation satellite system (GNSS); acquiring coordinates of a ship antenna in an Earth-Centered Earth-Fixed (ECEF) coordinate system using a ship-mounted Global Navigation Satellite System (GNSS) receiver The acquired coordinates of the ship antenna in the ECEF coordinate system are represented as: ； wherein, denotes the time of coordinate acquisition; , , denote the projection distances of the ship body antenna in the X-axis, Y-axis, and Z-axis directions in the geocentric and geodetic coordinate system, respectively. Converts the input point to the coordinates in the East-North-Up (ENU) coordinate system. Converts the input point to the coordinates in the East-North-Up (ENU) coordinate system. : ； wherein the coordinates of the phase center of the ship antenna in the ECEF coordinate system are selected the coordinates of the reference point, the longitude of the reference point, the latitude of the reference point, denotes the component of the ship antenna in the east direction, denotes the component of the ship antenna in the north direction, denotes the component of the ship antenna in the vertical direction.

3. The method of claim 2, wherein, S1 further comprises: S1.2 coordinates of the hull antenna in the ENU coordinate system , in combination with the hull attitude and the position of the transducer relative to the GNSS antenna in the gyro Cartesian coordinate system, to calculate the target time position of the shipboard transducer in the local ENU coordinate system at the time : ； ； ； ； ； ； wherein is a time-varying rotation matrix, denotes the ship's attitude, , , denote the ship's roll, pitch and heading attitude, respectively, denotes the position of the shipboard transducer relative to the GNSS antenna in the gyroscope rectangular coordinate system, , , denote the projection components on the X, Y, Z axes in the corresponding rectangular coordinate system, , denotes an intermediate calculation matrix.

4. The method of claim 3, wherein, S1 further comprises: S1.3 Obtain the approximate coordinates of the seabed transponder based on the principle of distance intersection by average sound speed method; The approximate coordinates of the seabed transponder are set as , observation data of observation epochs are selected, and the approximate coordinates of the seabed transponder are calculated based on the observation data of each observation epoch by using a least square method, and the specific expression is as follows: ； wherein is the average sound speed calculated from the sound speed profile, is the time delay corresponding to each observation epoch, , , represents the coordinates of the observation point for each observation epoch.

5. The method of claim 4, wherein, S1 further comprises: S1.4 approximate coordinates based on subsea transponder and the position of the ship-borne transducer at the time instant in a local ENU coordinate system calculating an approximate angle of incidence, the approximate angle of incidence being: ； wherein denotes the approximate angle of incidence.

6. The method of claim 5, wherein the weight function is defined as Establishes the iterative formula about Snell constant, transforms the approximate incident angle according to Snell's law of refraction, and obtains the initial value of Snell constant iteration : ； wherein, is the incident layer sound speed; The weight function is constructed or selected, the obtained approximate incident angle is brought into the weight function, and the Newton's descent factor is solved .

7. The method of claim 6, wherein the weight function is defined as: ###0001### where is the angle of incidence of the sound ray, and is the angle of incidence of the sound ray at the previous time step. The intralayer time delay calculation formula related to the Snell constant is established according to the intralayer time delay calculation formula related to the incidence angle in the sound ray tracing algorithm as follows: ； wherein denotes the layer's Snell constant, is the layer's sound speed gradient; and are the sound speeds of the layer and layer, respectively, is the Snell constant; The depth between the shipborne transducer and the seabed transponder is obtained by the depth gauge, and the sound speed profile of the water area where the shipborne transducer and the seabed transponder are located is obtained by the sound speed profiler. Based on the obtained depth from the shipborne transducer to the seabed transponder and the sound speed profile, the sound speed profile layer through which the sound ray passes from the shipborne transducer to the seabed transponder is constructed.

8. The method of claim 7, wherein, Based on the sound speed profile layer through which the sound ray passes from the shipborne transducer to the seabed transponder and the intralayer time delay calculation formula of the Snell constant, the time delay calculation formula related to the Snell constant is established: ； Convert the time delay calculation formula into the form of solving the zero point of the function, which is equivalent to obtaining the function expression of the time delay calculation formula, and calculate the first derivative of the function expression: ； ； wherein, represents the number of layers of the sound speed profile profile through which the acoustic ray travels from the shipboard transducer to the seafloor transponder, is a function expression of the time delay calculation formula, is the first derivative of 9. The method of claim 8, wherein, Construct the adaptive Newton's descent iteration formula, and substitute the obtained initial value of the Snell constant iteration, the Newton's descent factor, the function expression of the time delay calculation formula and the first derivative of the function expression of the time delay calculation formula into the adaptive Newton's descent iteration formula, and the adaptive Newton's descent iteration formula is constructed as follows: ； At this time, , is the iteration output value; Based on the adaptive Newton iterative formula for iterative operation, the iterative final value of the Snell constant is obtained ; If the iteration output value meets the iteration termination condition, let ; Else let and repeat the iteration until the iteration termination condition is met.

10. The method of claim 9, wherein, When the chosen initial Newton down-hill factor leads to iteration failure, let , restart iteration with a new factor, and if still failure, halve , until iteration succeeds; Using Snell's law of refraction, we have Transforming to , we get the sound ray tracing incidence angle .

Images