Sound ray tracking incident angle calculation method based on weighting function

Abstract

The invention discloses a sound ray tracking incident angle calculation method based on a weighting function, relates to the technical field of ocean underwater acoustic positioning, is used for obtaining a sound ray tracking incident angle in the process of realizing accurate positioning of a seabed control point by depending on a sound ray tracking algorithm, and comprises the following steps: solving an approximate incident angle based on coordinates of a shipborne transducer and a seabed transponder; converting the initial incident angle through the Snell refraction law to obtain an iteration initial value; calculating the initial incidence angle through a weighting function to obtain a Newton downhill factor; constructing a sound velocity profile layer through a depth meter and a sound velocity profiler to obtain a time delay calculation formula; and constructing a Newton downhill method iterative formula, and calculating to obtain a sound ray tracking incident angle based on the obtained approximate incident angle, an iterative initial value, a Newton downhill factor and a time delay calculation formula. According to the method, a traditional Newton iteration method is improved, the stability and iteration efficiency of solving the sound ray tracking incident angle through the Newton iteration method are improved, and the problem that sound ray tracking cannot be carried out when the incident angle is unknown is solved.

Description

Technical Field

[0001] This invention relates to the field of marine underwater acoustic positioning technology, mainly to the field of sound ray tracking algorithm optimization technology for underwater acoustic positioning systems, specifically to a sound ray tracking incident angle calculation method based on a weighted function. Background Technology

[0002] Submarine geodetic reference networks provide navigational support for various unmanned surface and underwater devices, and can also be used to monitor dynamic changes in seabed tectonic plates and the aquatic environment. They are crucial infrastructure for marine security, marine economic development, and marine environmental monitoring. 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. It uses the principle of underwater acoustic ranging to transmit the position of acoustic transducers to seabed control points, and is a key means of constructing and maintaining seabed geodetic reference networks. Accurate determination of the location of seabed control points requires precise measurement of the sound velocity in seawater.

[0003] However, the speed of sound in seawater is affected by seawater temperature, salinity, and pressure, and is a function of time and space, exhibiting time-varying and space-varying characteristics. Treating the speed of sound as a fixed value introduces systematic errors into the ranging results. To eliminate the influence of sound speed errors, the spatial symmetry and observation synchronization of the transponder can be utilized to eliminate the impact of sound speed variations; drawing on the differential technology in GNSS (Global Navigation Satellite System), single-difference and double-difference methods can be used to eliminate the influence of common errors such as long-period systematic errors when sufficient observation data is available; error compensation algorithms can be used to eliminate the influence of time-varying errors in the sound speed profile and horizontal heterogeneity; and ray-tracking algorithms can be used to directly correct the slant range. The implementation of these error elimination methods all rely on ray-tracking algorithms.

[0004] In the process of locating seabed control points using ray tracing algorithms, the propagation direction of sound rays is determined by both the sound speed and the angle of incidence at the current location. Furthermore, the geometric constraints and physical laws governing sound ray propagation are strongly dependent on the angle of incidence; if the angle of incidence is unknown, ray tracing cannot be performed. Existing ray tracing algorithms for numerical analysis mainly focus on the bisection method, the secant method, Newton's iteration method, and their improved versions. Their main function is to gradually approximate the true path or relevant parameters (such as the angle of incidence, propagation time, and horizontal distance) of sound ray propagation through numerical iteration, thereby solving the problem of solving nonlinear equations in ray tracing. Among these, the bisection method continuously divides the search interval into two, determining the sub-interval containing the root based on the change in the function value's sign, gradually narrowing the range until the accuracy requirement is met. The bisection method is insensitive to the initial interval selection and is suitable for coarse localization problems in complex acoustic models, but not for precise localization of seabed control points. The secant method uses the slope of the line connecting two points to approximate the derivative, approximating the root through an iterative formula. It only requires function value calculation, making it suitable for scenarios where derivatives are difficult to express analytically in acoustic models. However, it requires two initial values, and drastic changes in the function near the iteration point may lead to convergence failure, resulting in poor robustness. Newton's iteration method approximates the root using the tangent line to the function at the iteration point. It has fast convergence speed and controllable accuracy. The high-precision requirements of acoustic positioning can be met by adjusting the number of iterations or the tolerance, making it widely applicable. However, Newton's iteration method is sensitive to initial values ​​and requires correct selection of the incident angle.

[0005] Among the many algorithms for solving the incident angle, there is relatively little research on the failure of the Newton iteration method when the incident angle is large. The Newton downhill method can reduce the impact of the initial value of the iteration on the solution stability by reducing the step size. However, the traditional Newton downhill method usually requires the search for the downhill factor through trial calculation (such as the halving method), which increases the amount of computation and reduces the convergence efficiency. On the other hand, directly selecting the initial downhill factor based on experience lacks theoretical basis and adaptability.

[0006] Therefore, a method for calculating the ray-tracking incident angle that can adaptively adapt to the Newton downhill factor and calculate the ray-tracking incident angle while ensuring successful iteration of the incident angle is needed to optimize the ray-tracking algorithm in underwater acoustic positioning systems. Summary of the Invention

[0007] This invention addresses the divergence problem of the existing traditional Newton's iteration method at large incident angles and the inability to perform ray tracking when the incident angle is unknown. It provides a method for calculating the incident angle of ray tracking based on a weighted function. By constructing an iterative formula for Newton's downhill method, the method calculates the ray tracking incident angle based on the initial value of the Snell constant, the Newton's downhill factor, and the time delay formula of the Snell constant. By introducing the Newton's downhill factor, the iteration step size is changed, reducing the initial value requirement of the Newton's iteration method at large incident angles, thus improving iteration efficiency while ensuring successful iteration.

[0008] A method for calculating the incident angle of ray tracking based on a weighted function includes: S1. Obtain the approximate coordinates of the seabed transponder and the coordinates of the shipborne transducer at the target time, and calculate the approximate incident angle based on the obtained coordinates; S2. The approximate incident angle is converted using Snell's law of refraction to obtain the initial value of the Snell constant for iteration; S3. Construct or select a weighting function, and use the weighting function to transform the approximate incident angle to obtain the Newton descent factor; S4. Establish an intra-layer delay calculation formula based on the ray tracing algorithm using the Snell constant; S5. Obtain the depth from the shipborne transducer to the seabed transponder, and calculate the sound velocity profile layers that the sound rays pass through from the shipborne transducer to the seabed transponder. S6. Based on the intra-layer time delay calculation formula of the sound velocity profile layer and the Snell constant, establish the time delay calculation formula of the Snell constant; S7. Construct the iterative formula for Newton's downhill method, and calculate the incident angle of the sound ray tracing based on the initial value of the Snell constant iteration, the Newton's downhill factor, and the time delay calculation formula of the Snell constant.

[0009] S1 includes: S1.1 Obtain the coordinates of the ship's antennas based on the Global Navigation Satellite System (GNSS); The coordinates of the ship's antenna in the geocentric-ground-fixed coordinate system ECEF are obtained using the GNSS receiver mounted on the ship. The coordinates of the ship's antenna in the geocentric coordinate system are expressed as follows: ; in, Indicates the time when the coordinates were acquired. , , These represent the projected distances of the ship's antenna in the X, Y, and Z axes, respectively, in the geocentric coordinate system. Will Convert to ENU coordinate system (Northeast Celestial) : ; Among them, the coordinates of the phase center of the ship's antenna in the ECEF coordinate system are selected. The coordinates of the reference point, For reference point longitude, For reference point latitude, This indicates the eastward component of the ship's antenna. This indicates the northward component of the ship's antenna. This indicates the vertical component of the ship's antenna.

[0010] S1 also includes: S1.2 Based on the coordinates of the ship's antenna in the ENU coordinate system By combining the ship's attitude and the transducer's position relative to the GNSS antenna in the gyroscope's Cartesian coordinate system, the target's time is calculated. Position of the shipborne transducer in the local ENU coordinate system : ; ; ; ; ; ; in, It is a time-varying rotation matrix. Indicates the ship's attitude. , , These represent the ship's roll, pitch, and heading attitudes, respectively. This indicates the position of the shipborne transducer relative to the GNSS antenna in the gyroscope's Cartesian coordinate system. , , They represent The projection components on the X, Y, and Z axes in the corresponding Cartesian coordinate system , This represents the intermediate calculation matrix. This indicates the matrix transpose.

[0011] S1 also includes: S1.3 Based on the principle of distance intersection, the approximate coordinates of the seabed transponder are obtained by the average speed of sound method; The approximate coordinates of the seabed transponder are set as follows: Select Based on the observation data from each observation epoch, the approximate coordinates of the seabed transponder are obtained using the least squares method. The specific expression is as follows: ; in, The average speed of sound calculated from the sound speed profile. The time delay corresponding to each observation epoch, ( , , () represents the coordinates of the observation point at each observation epoch.

[0012] S1 also includes: S1.4 Based on the approximate coordinates of the seabed transponder and Position of the shipborne transducer in the local ENU coordinate system Calculate the approximate angle of incidence. The approximate angle of incidence is: ; in, This represents the approximate angle of incidence.

[0013] An iterative formula for the Snell constant is established, and the approximate incident angle is transformed according to Snell's law of refraction to obtain the initial value of the Snell constant for iteration. : ; in, The velocity of sound in the incident layer; Construct or select a weighting function, substitute the obtained approximate incident angle into the weighting function, and solve for the Newton descent factor. .

[0014] 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: ; 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; 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.

[0015] 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: ; 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: ; ; 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.

[0016] 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: ; at this time, , For the iterative output value; Iterative calculations are performed based on the adaptive Newton's downhill method iterative formula to obtain the final value of the iteration. ; If the iteration output value satisfies the iteration termination condition, then let Otherwise, The iteration is repeated until the iteration termination condition is met.

[0017] When the selected initial Newton's downhill factor causes the iteration to fail, then let Restart the iteration; if it still fails, halve the number of iterations. Continue until the iteration is successful; Using Snell's law of refraction Convert to The incident angle for ray tracing is obtained. .

[0018] Compared with the prior art, the present invention has the following advantages: This invention provides a method for calculating the incident angle of acoustic ray tracking based on a weighted function. By constructing an iterative formula for Newton's descent method, and based on the initial value of the Snell constant, the adaptive Newton's descent factor, and the time delay calculation formula for the Snell constant, an adaptive Newton's descent method iterative formula is constructed. Iterative calculation is performed based on the adaptive Newton's descent method iterative formula to obtain the incident angle of acoustic ray tracking. This solves the problem that acoustic ray tracking cannot be performed if the incident angle of the acoustic ray is unknown in the process of positioning seabed control points relying on acoustic ray tracking algorithms.

[0019] This invention provides a method for calculating the incident angle of ray tracking based on a weighted function. It improves upon the traditional Newton iteration method by incorporating the weights into the Newton descent factor, analogous to the weights in the adjustment stochastic model. This enables adaptive adjustment of the Newton descent factor, reducing the iteration step size at large incident angles and improving iteration efficiency while ensuring successful iteration at the incident angle. Attached Figure Description

[0020] Figure 1 This is a flowchart of a method for calculating the incident angle of ray tracing based on a weighting function.

[0021] Figure 2 This is a schematic diagram of the navigation path.

[0022] Figure 3 This is a schematic diagram of the sound velocity profile.

[0023] Figure 4 This is a schematic diagram of the weighting function.

[0024] Figure 5 This is a schematic diagram of a voice tracing algorithm.

[0025] Figure 6 This is a schematic diagram of the residual changes of different weighted models at the 332nd observation epoch.

[0026] Figure 7 This is a schematic diagram showing the change in the number of iterations for different weight models. Detailed Implementation

[0027] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention are described clearly and completely below. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0028] like Figure 1As shown, a method for calculating the incident angle of ray tracking based on a weighted function includes: S1. Obtain the approximate coordinates of the seabed transponder. and target time Coordinates of the shipborne transducer Using the approximate coordinates of the seabed transponder and the target time The coordinates of the shipborne transducer at any given time are calculated to obtain the approximate incident angle. S2. Convert the approximate incident angle using Snell's law of refraction to obtain the initial value of the Snell constant for iteration. ; S3. Construct or select a weighting function, and use the weighting function to transform the approximate incident angle to obtain the Newton descent factor; S4. Establish an intra-layer delay calculation formula based on the ray tracing algorithm using the Snell constant; S5. Obtain the depth from the shipborne transducer to the seabed transponder, and calculate the sound velocity profile layers that the sound rays pass through from the shipborne transducer to the seabed transponder. S6. Based on the intra-layer time delay calculation formula of the sound velocity profile layer and the Snell constant, establish the time delay calculation formula of the Snell constant; S7. Construct the iterative formula for Newton's downhill method, and calculate the incident angle of the sound ray tracing based on the initial value of the Snell constant iteration, the Newton's downhill factor, and the time delay calculation formula of the Snell constant.

[0029] The approximate coordinates of the seabed transponder represent its static position in the Earth coordinate system. It serves as the reference point for acoustic positioning and acts as a fixed reference point for the acoustic positioning system to determine the relative position of shipborne equipment or targets.

[0030] Shipborne transducer coordinates are the three-dimensional geographic coordinates of an acoustic transducer (a device that transmits / receives sound waves) installed on the ship's hull at a given moment. These coordinates change dynamically with the ship's movement, reflecting the transducer's spatial position at the instant of measurement. Real-time calculation of the shipborne transducer coordinates involves obtaining the ship's antenna coordinates via GNSS, combining this with the ship's attitude sensor (IMU) and the transducer's installation offset, to calculate the transducer coordinates in real time.

[0031] S1 includes: S1.1 Obtain the coordinates of the ship's antennas based on the Global Navigation Satellite System (GNSS); First, the coordinates of the ship's antenna in the geocentric-ground-fixed coordinate system ECEF are acquired in real time using a GNSS receiver mounted on the ship. The coordinates of the ship's antenna in the geocentric coordinate system are expressed as follows: ; in, Indicates the time when the coordinates were acquired. , , These represent the projected distances of the ship's antenna in the X, Y, and Z axes, respectively, in the geocentric coordinate system. Then Convert to ENU coordinate system (Northeast Celestial) : ; Among them, the coordinates of the phase center of the ship's antenna in the ECEF coordinate system are selected. The coordinates of the reference point, For reference point longitude, For reference point latitude, This indicates the eastward component of the ship's antenna. This indicates the northward component of the ship's antenna. This indicates the vertical component of the ship's antenna.

[0032] The GNSS receiver can acquire the position of the survey vessel in real time. By determining the relative positional relationship between the transducer on the bottom of the vessel and the GNSS receiver, the position of the acoustic transducer can be calculated based on coordinate correction.

[0033] The ECEF coordinate system uses the Earth's center as its origin. The coordinate values ​​of the ship's antenna are typically in the millions to tens of millions of meters. Directly using these large values ​​for calculations can easily lead to numerical overflow or loss of accuracy. The ENU coordinate system, a local coordinate system, uses the ship's antenna as its origin and converts the coordinates into local values ​​in meters. By transforming the coordinates in the ECEF coordinate system into those in the ENU coordinate system, large numerical calculations are converted into smaller numerical calculations. This avoids the complex trigonometric function calculations in the ECEF coordinate system and significantly improves computational stability and accuracy.

[0034] S1 also includes: S1.2 Based on the coordinates of the ship's antenna in the ENU coordinate system By combining the ship's attitude and the transducer's position relative to the GNSS antenna in the gyroscope's Cartesian coordinate system, the target's time is calculated. Position of the shipborne transducer in the local ENU coordinate system : ; ; ; ; ; ; in, It is a time-varying rotation matrix. Indicates the ship's attitude. , , These represent the ship's roll, pitch, and heading attitudes, respectively. This indicates the position of the shipborne transducer relative to the GNSS antenna in the gyroscope's Cartesian coordinate system. , , They represent The projection components on the X, Y, and Z axes in the corresponding Cartesian coordinate system , This represents the intermediate calculation matrix.

[0035] The position of the transducer in the local ENU coordinate system is calculated using the method described above. like Figure 2 As shown. The GNSS receiver provides real-time surface position information for the survey vessel, while the acoustic depth sounder acquires the relative distance between the vessel and the seabed through interaction with the seabed transponder. The GNSS receiver can acquire the survey vessel's position information in real time, and the motion sensor can simultaneously measure the vessel's attitude information. In the vessel's coordinate system, the position of the acoustic transducer can be determined based on the positional relationship between the GNSS antenna and the transducer, and the vessel's motion attitude sensing data. Figure 2 The dashed line represents the survey vessel's navigation path, equivalent to the position of the shipborne transducer in the local ENU coordinate system at different times, i.e., the movement path of the shipborne transducer in the horizontal plane. The pentagram represents the position of the M01 transponder. As an underwater acoustic device, the transponder is usually installed in a fixed position on the seabed. When the transponder receives an acoustic signal from the ship, it transmits a response signal back to the ship. By calculating the round-trip time, speed of sound, and other information of the acoustic signal, the distance between the ship and the transponder can be calculated, thereby determining the ship's position relative to the transponder.

[0036] S1 also includes: S1.3 Based on the principle of distance intersection, the approximate coordinates of the seabed transponder are obtained by the average speed of sound method; The approximate coordinates of the seabed transponder are set as follows: Select Based on the observation data from each observation epoch, the approximate coordinates of the seabed transponder are obtained using the least squares method. The specific expression is as follows: ; in, According to Figure 3 The schematic diagram of the sound velocity profile shown is the average sound velocity obtained by averaging the sound velocity profile. The time delay corresponding to each observation epoch, ( , , () represents the coordinates of each observation epoch.

[0037] Figure 3 This is a schematic diagram of the sound speed profile in the sea from 0 to 4000 meters, showing how the sound speed changes with seawater depth. Figure 3 As shown, the horizontal axis represents the speed of sound, in meters per second (m / s), and the vertical axis represents depth, in meters (m). From Figure 3 It can be seen that the change in sound speed with depth is not linear, but rather exhibits a clear stratified characteristic, with sound speed varying at a certain depth. There is a minimum value (sound channel axis), where the sound speed at the sea surface is greater than the sound speed at the seabed, indicating a fully deep-sea sound channel. Above the sound channel axis, the sound speed increases with increasing seawater depth and changes rapidly. Below the sound channel axis, the sound speed decreases with increasing seawater depth, but in deep sea areas (where the seawater depth is greater than 1500 meters), the sound speed tends to stabilize with seawater changes, equivalent to the slope of the sound speed line gradually approaching a constant in a sound speed profile.

[0038] An observation epoch refers to the point in time or period during GNSS measurements when satellite signals are observed at a specific moment.

[0039] The average sound velocity method likely refers to considering the variation of sound velocity in different water layers and taking the average value to calculate distance. Alternatively, a more specific implementation of the average sound velocity method might require considering sound velocity profiles, such as dividing the water layer into multiple layers with different sound velocities, and then taking the average.

[0040] S1 also includes: S1.4 Based on the approximate coordinates of the seabed transponder and Position of the shipborne transducer in the local ENU coordinate system Calculate the approximate angle of incidence. The approximate angle of incidence is: ; in, This represents the approximate angle of incidence.

[0041] An iterative formula for the Snell constant is established, and the approximate incident angle is transformed according to Snell's law of refraction to obtain the initial value of the Snell constant for iteration. : ; in, The velocity of sound in the incident layer; Construct or select a weighting function, substitute the obtained approximate incident angle into the weighting function, and solve for the Newton descent factor. .

[0042] The Newton's downhill factor is an adaptive Newton's downhill factor, as shown in Table 1. The Newton's downhill factor has multiple weight function expressions: Table 1. Weighting function representation of the Newton's downhill factor .

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

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

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

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

[0047] 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: ; 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; 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.

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

[0049] 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: ; 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: ; ; in, This indicates the number of sound velocity profile layers that the sound ray passes through from the shipboard transducer to the seabed transponder. This is the zero-point form of the solution function for the one-way time delay calculation formula. It is the first derivative of the function in the form of the zero point of the formula for calculating one-way time delay.

[0050] An adaptive Newton's downhill method iterative formula is constructed. The obtained initial value of the Snell constant, the Newton's downhill factor, the functional expression of the one-way time delay calculation formula, and the first derivative of the functional expression of the one-way time delay calculation formula are substituted into the adaptive Newton's downhill method iterative formula. Iterative calculations are then performed based on the adaptive Newton's downhill method iterative formula to obtain the final value. : ; at this time, , For the iterative output value; If the iteration output value satisfies the iteration termination condition, then let Otherwise, The iteration is repeated until the iteration termination condition is met. The iteration terminates when the vertical displacement of the sound ray is the same as the depth from the shipborne transducer to the seabed transponder.

[0051] When the selected initial Newton's downhill factor causes the iteration to fail, then let Restart the iteration; if it still fails, halve the number of iterations. Continue until the iteration is successful; use Snell's law to... Convert to The incident angle for ray tracing is obtained. .

[0052] In real sea trial scenarios, the incident angles calculated using different weighting functions are approximately the same. However, the approximate incident angle, which is equivalent to the incident angle obtained by the average sound speed method, is often smaller than the true incident angle. Moreover, the larger the incident angle, the more obvious the difference between the incident angle calculated by the average sound speed method and the incident angle calculated by the weighting function. This indicates that the initial value of the Snell constant obtained based on the approximate incident angle is further away from the true value. Therefore, when solving for the incident angle, the case of an excessively large incident angle must be considered.

[0053] The computational residuals of the models solved using different weighting functions are significantly lower than those of the average speed of sound method. The unit weight error of the ray tracking models established using different incident angle calculation methods is reduced by approximately 0.7 meters, achieving the goal of ray correction and effectively eliminating the influence of ray bending. Therefore, the adaptive Newton's descent method using different incident angle weighting functions exhibits higher accuracy in solving for the incident angle compared to the average speed of sound method. Furthermore, the variation of the computational residuals with observation epochs is related to the magnitude of the incident angle calculated for each observation epoch; the larger the incident angle, the larger the absolute value of the residual.

[0054] Figure 6 This is a schematic diagram illustrating the residual changes of different weighted models at epoch 332, as shown below. Figure 6 As shown, the curve of the cosine function model is at the top, the absolute value of the residual is the smallest, and the obtained incident angle is the most accurate. Among the selected functions, the cosine function best matches the change of the iteration step size and has the best effect. Therefore, the Newton's downhill method based on the cosine weight function obtains the best effect when solving for the incident angle.

[0055] The convergence speed of different models was evaluated, and the number of iterations required for different models was as follows: Figure 7 As shown, for the same function model (lines of the same color), the number of iterations varies with the observed epoch and is related to the size of the incident angle at each epoch. The larger the incident angle, the more iterations are required for convergence, and the slower the convergence. For different function models (lines of different colors) at the same observed epoch, the corresponding function curve of the cosine function model is located at... Figure 7 Therefore, the cosine function model has the fewest iterations and the fastest convergence. Thus, the Newton's downhill method based on the cosine weight function achieves the fastest convergence speed when solving for the incident angle, and has the fewest iterations in most epochs.

[0056] The above are merely preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for calculating the incident angle of ray tracking based on a weighted function, characterized in that, include: S1. Obtain the approximate coordinates of the seabed transponder and the coordinates of the shipborne transducer at the target time, and calculate the approximate incident angle based on the obtained coordinates; S2. The approximate incident angle is converted using Snell's law of refraction to obtain the initial value of the Snell constant for iteration; S3. Construct or select a weighting function, and use the weighting function to transform the approximate incident angle to obtain the Newton descent factor; S4. Establish an intra-layer delay calculation formula based on the ray tracing algorithm using the Snell constant; S5. Obtain the depth from the shipborne transducer to the seabed transponder, and calculate the sound velocity profile layers that the sound rays pass through from the shipborne transducer to the seabed transponder. S6. Based on the intra-layer time delay calculation formula of the sound velocity profile layer and the Snell constant, establish the time delay calculation formula of the Snell constant; S7. Construct the iterative formula for Newton's downhill method, and calculate the incident angle of the sound ray tracing based on the initial value of the Snell constant iteration, the Newton's downhill factor, and the time delay calculation formula of the Snell constant.

2. The method for calculating the incident angle of ray tracking based on a weighted function according to claim 1, characterized in that, S1 includes: S1.1 Obtain the coordinates of the ship's antennas based on the Global Navigation Satellite System (GNSS); The coordinates of the ship's antenna in the geocentric-ground-fixed coordinate system ECEF are obtained using the GNSS receiver mounted on the ship. The coordinates of the ship's antenna in the geocentric coordinate system are expressed as follows: ； in, Indicates the time when the coordinates were acquired; , , These represent the projected distances of the ship's antenna in the X, Y, and Z axes, respectively, in the geocentric coordinate system. Will Convert to ENU coordinate system (Northeast Celestial) : ； Among them, the coordinates of the phase center of the ship's antenna in the ECEF coordinate system are selected. The coordinates of the reference point, For reference point longitude, For reference point latitude, This indicates the eastward component of the ship's antenna. This indicates the northward component of the ship's antenna. This indicates the vertical component of the ship's antenna.

3. The method for calculating the incident angle of ray tracking based on a weighting function according to claim 2, characterized in that, S1 also includes: S1.2 Based on the coordinates of the ship's antenna in the ENU coordinate system By combining the ship's attitude and the transducer's position relative to the GNSS antenna in the gyroscope's Cartesian coordinate system, the target's time is calculated. Position of the shipborne transducer in the local ENU coordinate system : ； ； ； ； ； ； in, It is a time-varying rotation matrix. Indicates the ship's attitude. , , These represent the ship's roll, pitch, and heading attitudes, respectively. This indicates the position of the shipborne transducer relative to the GNSS antenna in the gyroscope's Cartesian coordinate system. , , They represent The projection components on the X, Y, and Z axes in the corresponding Cartesian coordinate system , This represents the intermediate calculation matrix.

4. The method for calculating the incident angle of ray tracking based on a weighting function according to claim 3, characterized in that, S1 also includes: S1.3 Based on the principle of distance intersection, the approximate coordinates of the seabed transponder are obtained by the average speed of sound method; The approximate coordinates of the seabed transponder are set as follows: Select Based on the observation data from each observation epoch, the approximate coordinates of the seabed transponder are obtained using the least squares method. The specific expression is as follows: ； in, The average speed of sound calculated from the sound speed profile. The time delay corresponding to each observation epoch, ( , , () represents the coordinates of the observation point at each observation epoch.

5. The method for calculating the incident angle of ray tracking based on a weighted function according to claim 4, characterized in that, S1 also includes: S1.4 Based on the approximate coordinates of the seabed transponder and Position of the shipborne transducer in the local ENU coordinate system Calculate the approximate angle of incidence. The approximate angle of incidence is: ； in, This represents the approximate angle of incidence.

6. The method for calculating the incident angle of ray tracking based on a weighting function according to claim 5, characterized in that, An iterative formula for the Snell constant is established, and the approximate incident angle is transformed according to Snell's law of refraction to obtain the initial value of the Snell constant for iteration. : ； in, The velocity of sound in the incident layer; Construct or select a weighting function, substitute the obtained approximate incident angle into the weighting function, and solve for the Newton descent factor. .

7. The method for calculating the incident angle of ray tracking based on a weighted function according to claim 6, characterized in that, 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: ； 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; 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.

8. The method for calculating the incident angle of ray tracing based on a weighting function according to claim 7, characterized in that, 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 time delay calculation formula related to the Snell constant is established: ； 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: ； ； 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.

9. The method for calculating the incident angle of ray tracking based on a weighting function according to claim 8, characterized in that, 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: ； at this time, , For the iterative output value; Iterative calculations were performed using the adaptive Newton's downhill method iterative formula to obtain the final value of the Snell constant. ; If the iteration output value satisfies the iteration termination condition, then let ; Otherwise The iteration is repeated until the iteration termination condition is met.

10. The method for calculating the incident angle of ray tracking based on a weighting function according to claim 9, characterized in that, When the selected initial Newton's downhill factor causes the iteration to fail, let Restart the iteration; if it still fails, halve the number of iterations. Continue until the iteration is successful; Using Snell's law of refraction Convert to The incident angle for ray tracing is obtained. .

Images