A Combined Ray Tracing Method for Improving the Precision of the USBL Positioning System in Smart Ocean

The ultra-short baseline positioning system (USBL) has the advantages of flexible application and easy installation, and it plays an extremely important role in the underwater positioning and communication. The error of the USBL in underwater positioning is mainly caused by a ranging error due to ray tracing, a phase difference error of the USBL, and acoustic noise in the underwater communication. Most of these errors are related to the changes in the sound speed during its propagation through the ocean. Therefore, when using the USBL for underwater detection, it is necessary to correct the sound speed profile in the detection area and optimize the ray tracing. Taking into account the actual conditions, this paper aims at correcting the model of underwater sound speed propagation and improving the tracking method of sound lines when the marine environment in the shallow sea area changes. This paper proposes a combined ray tracing method that can adaptively determine whether to use the constant sound speed ray tracing method or the equal gradient ray tracing method. The theoretical analysis and simulation results show that the proposed method can effectively reduce the error of slant distance in USBL compared with the traditional acoustic tracking method and the constant sound speed ray tracing method. The proposed sound ray correction algorithm solves the contradiction between the number of iterations and the reduction of positioning error and has engineering application value.


Introduction
The ultra-short baseline positioning system (USBL) is based on the installation of an ultra-short baseline array on a ship bottom and a transponder mounted on an underwater target. The transmitter emits an acoustic signal, and the underwater target receives the transmitted acoustic signal and sends back an answer signal. After receiving the signal, the USBL calculates the azimuth and distance [1] of the target using the time-delay difference or phase difference between each receiving hydrophone [2], and calculates the position coordinates of the underwater target. The ultra-short baseline positioning system (USBL) is generally used for underwater positioning in shallow sea area [1].
Due to the complexity of the actual marine environment, the accuracy of an ultra-short baseline positioning system is limited mainly due to the existence of the ranging error, phase difference error, transmission and reception delay estimation error, and acoustic noise in the communication. The main where d is the spacing between the array elements, 12  and 22  are the phase difference between the second and the other two hydrophones, and H is the height of the target array element, X is the distance of the target along the x-axis direction, Y is the distance of the target along the y-axis direction, l is the slant distance of the target from the USBL, and triple (X, Y, H) denotes the threedimensional coordinates of the target relative to the array. Sound propagation speed has a very little effect on the phase difference of the signals received by the hydrophones in the transceiver head. One of the main factor affecting the positioning accuracy of the USBL is the measuring error caused by sound propagation in the slant range. Therefore, it is particularly necessary to take into account the form of the sound propagation path to improve the slant range accuracy of the target and the array.

Empirical Orthogonal Function
The sound underwater propagation path and attenuation are closely related to the change of underwater sound speed. The change of the sound speed profile (SSP) is one of the important factors affecting the sound propagation in the ocean [15]. The sound speed profile is a structural distribution of the path and velocity variations of an acoustic wave propagating in the seawater. Due to the marine environment misalignment caused by the environmental factors such as the non-uniformity of seawater medium, and acoustic propagation channel interference [16], climate change [17][18][19], and other various reasons such as failure of the equipment for sound speed measuring, the sound speed profile measured in real time has diversity and errors. The empirical orthogonal function is a very effective method for inversion of the sound velocity profile. The EOF can obtain accurate enough sound speed profiles with the history sound speed data. The research of Davis et al. [20] showed that the order of EOF is associated with the ocean physics processes, and usually, only the first few orders are needed to represent the sound speed profile accurately. Sound propagation speed has a very little effect on the phase difference of the signals received by the hydrophones in the transceiver head. One of the main factor affecting the positioning accuracy of the USBL is the measuring error caused by sound propagation in the slant range. Therefore, it is particularly necessary to take into account the form of the sound propagation path to improve the slant range accuracy of the target and the array.

Empirical Orthogonal Function
The sound underwater propagation path and attenuation are closely related to the change of underwater sound speed. The change of the sound speed profile (SSP) is one of the important factors affecting the sound propagation in the ocean [15]. The sound speed profile is a structural distribution of the path and velocity variations of an acoustic wave propagating in the seawater. Due to the marine environment misalignment caused by the environmental factors such as the non-uniformity of seawater medium, and acoustic propagation channel interference [16], climate change [17][18][19], and other various reasons such as failure of the equipment for sound speed measuring, the sound speed profile measured in real time has diversity and errors. The empirical orthogonal function is a very effective method for inversion of the sound velocity profile. The EOF can obtain accurate enough sound speed profiles with the history sound speed data. The research of Davis et al. [20] showed that the order of EOF is associated with the ocean physics processes, and usually, only the first few orders are needed to represent the sound speed profile accurately.
The EOF utilizes the spatiotemporal correlation of sound velocity profile data. After feature vector extraction and modal decomposition of the group of the sound speed profiles with the same characteristics, a common feature vector is obtained and combined with the sampled data, so that an accurate representation of the continuous sound velocity profile data can be achieved [21].
Assuming there are M measured sound velocity profiles in the measurement area, M sound speed profile samples can be obtained, which are expressed as where N is the number of discrete points in the depth direction. These samples can be expressed as [22] The average of M speed profiles is obtained by In the sound speed matrix C, M sound speed profiles are subtracted from the average sound speed profile c respectively to obtain the perturbation of each sound speed profile with respect to the average sound speed profile; namely, the perturbation matrix ∆C is defined as where ∆c j (z i ) = c j (z i −c(z i ) i = 1, 2, · · · , N, j = 1, 2, · · · , M. The covariance matrix R of the sound speed perturbation matrix is defined by The covariance matrix R represents the uncertainty of the fluctuation in the sound speed in the measurement area. The characteristic decomposition of R results in In Equation (9), D = [λ 1 λ 2 · · · λ N ] is the eigenvalue matrix of the covariance matrix R, and F is the matrix of the eigenvectors of R, expressed as The component f i (z) corresponding to the eigenvalue λ i is the ith EOF mode. The above N feature values λ 1 , λ 2 , . . . , λ N are arranged in ascending order. The sound velocity profile at any point in the measurement region can be represented by the first K-order EOF approximation as In Equation (11), α i is the EOF coefficient of the sound speed profile. The determination method of α i is (11) can be used to represent the values of each sound velocity profile retrieved in the SSP. Normally, the sound velocity profile can be well inverted by selecting 3 to 6 orders of the EOF [23].

Basic Theory of Ray Acoustics in Layered Media
Ray acoustic is an intuitive sound field analysis method, which has the good robustness to the complex media and boundary inversion. Moreover, it is an important method to analyze the sound field. The basic assumption of the ray acoustics is that acoustic waves are transmitted along a certain direction, and the trajectory of acoustic waves is the sound line, and the ray lines are perpendicular to the isosceles (wavefronts) [24]. As illustrated in Figure 2, in a layered media (where there is a border between two media or the properties of the same media are changed), the acoustic wave will be reflected and refracted, and the propagation path of acoustic wave will always bend toward the area with a lower sound speed. The propagation law of sound lines cross over the medium satisfies the Snell's law, which is expressed as In Equation (11), i  is the EOF coefficient of the sound speed profile. The determination (11) can be used to represent the values of each sound velocity profile retrieved in the SSP. Normally, the sound velocity profile can be well inverted by selecting 3 to 6 orders of the EOF [23].

Basic Theory of Ray Acoustics in Layered Media
Ray acoustic is an intuitive sound field analysis method, which has the good robustness to the complex media and boundary inversion. Moreover, it is an important method to analyze the sound field. The basic assumption of the ray acoustics is that acoustic waves are transmitted along a certain direction, and the trajectory of acoustic waves is the sound line, and the ray lines are perpendicular to the isosceles (wavefronts) [24]. As illustrated in Figure 2, in a layered media (where there is a border between two media or the properties of the same media are changed), the acoustic wave will be reflected and refracted, and the propagation path of acoustic wave will always bend toward the area with a lower sound speed. The propagation law of sound lines cross over the medium satisfies the Snell's law, which is expressed as   In Equation (13), θ i is the angle between the sound wave of the ith layer and the normal to the horizontal plane, and it is called the sonic incident angle of the ith layer, c i is the sound speed of the ith layer, and p is the Snell constant.
In different layers of a propagation media, the sound speed changes, so the propagation path of the sound is not a straight line [25]. Using the Snell's law, we can develop a method for determining the sound path. Presently, two methods are commonly used for acoustic path tracking [26], the constant sound speed ray tracing method, where it is assumed that the sound speed in the layer is constant, and the equal gradient ray tracing, where it is assumed that there is a gradient change in the sound speed in the layer.

Constant Sound Speed Ray Tracing Method
Assume that the sound beam propagates in an underwater area composed of n layers with different sound speeds. Also, assume that a sound wave propagates along a straight line at a constant velocity in one layer, while the sound speed in each layer is different. Then, the sound wave's propagation trajectory can be subdivided into several straight lines [27], as shown in Figure 3. The path propagates in the ith layer beam is S i , the angle between the incident beam and the horizontal plane is θ i , c i is the propagation speed of a sound wave in each layer, t i is the wave propagation time through the ith layer, and z i is the water layer thickness of the ith layer. In different layers of a propagation media, the sound speed changes, so the propagation path of the sound is not a straight line [25]. Using the Snell's law, we can develop a method for determining the sound path. Presently, two methods are commonly used for acoustic path tracking [26], the constant sound speed ray tracing method, where it is assumed that the sound speed in the layer is constant, and the equal gradient ray tracing, where it is assumed that there is a gradient change in the sound speed in the layer.

Constant Sound Speed Ray Tracing Method
Assume that the sound beam propagates in an underwater area composed of n layers with different sound speeds. Also, assume that a sound wave propagates along a straight line at a constant velocity in one layer, while the sound speed in each layer is different. Then, the sound wave's propagation trajectory can be subdivided into several straight lines [27], as shown in Figure 3. The path propagates in the ith layer beam is i S , the angle between the incident beam and the horizontal plane c is the propagation speed of a sound wave in each layer, i t is the wave propagation time through the ith layer, and i z is the water layer thickness of the ith layer. According to the physical and geometrical functions, the trigonometric function of the acoustic wave at the incident angle i The horizontal displacement i x and time i t of the acoustic wave passing through the ith layer are defined by According to the physical and geometrical functions, the trigonometric function of the acoustic wave at the incident angle θ i of each medium is The horizontal displacement x i and time t i of the acoustic wave passing through the ith layer are defined by The total horizontal displacement X of the signal received from the underwater target is defined by This method is simple and iterative. Besides, there is a small amount of calculation iterations and fast ray tracing can be achieved. However, the error made by constant sound speed ray tracing method is larger than the equal gradient ray tracing method; especially in the ocean areas where the sound speed changes greatly, so an accurate tracing cannot be achieved.

Equal Gradient Ray Tracing Method
In the equal gradient ray tracing method, it is assumed that the relationship between the velocity of sound and the depth is linear, and the trajectory of a sound wave in each layer is considered as an arc segment with the radius r [28], as shown in Figure 4. In Figure 4, S i is the path of the arc segment through the ith layer, and r i is the radius of the circle corresponding to the radius of the ith layer; P is the Snell constant which is defined by Equation (13).
The total horizontal displacement X of the signal received from the underwater target is defined by This method is simple and iterative. Besides, there is a small amount of calculation iterations and fast ray tracing can be achieved. However, the error made by constant sound speed ray tracing method is larger than the equal gradient ray tracing method; especially in the ocean areas where the sound speed changes greatly, so an accurate tracing cannot be achieved.

Equal Gradient Ray Tracing Method
In the equal gradient ray tracing method, it is assumed that the relationship between the velocity of sound and the depth is linear, and the trajectory of a sound wave in each layer is considered as an arc segment with the radius r [28], as shown in Figure 4. In Figure 4, i S is the path of the arc segment through the ith layer, and i r is the radius of the circle corresponding to the radius of the ith layer; P is the Snell constant which is defined by Equation (13). The sound speed gradient i g in the ith layer can be expressed as where i g is the sound velocity gradient of the ith layer.
The trajectory of a sound wave in each layer can be regarded as a circular arc with a constant curvature; the radius i r and the length i S of the arc corresponding to the ith layer arc can be expressed as The sound speed gradient g i in the ith layer can be expressed as where g i is the sound velocity gradient of the ith layer. The trajectory of a sound wave in each layer can be regarded as a circular arc with a constant curvature; the radius r i and the length S i of the arc corresponding to the ith layer arc can be expressed as The horizontal displacement x i and time t i of the acoustic wave passing through the ith layer are The total horizontal displacement X of the signal received by the transmitter from the underwater target is The ray path calculated by this method is in line with the actual ray path, and the calculation accuracy is higher than that of the constant sound speed ray tracing method [29]. However, the number of iterations is large, making the calculation process time-consuming.

Combined Ray Tracing Method
The USBL is used in the shallow sea environment where the sound speed varies with depth, so the sound speed can no longer be regarded as a fixed value. Using a straight sound line to track and locate the target directly will lead to inaccurate positioning. The equal gradient ray tracing method can be used to improve the positioning accuracy. However, using the equal gradient ray tracing method throughout the entire propagation, although providing high accuracy, requires a large amount of calculation of iterative processing. In the process of moving underwater targets, the deepening of the underwater depth will cause changes in the sound propagation speed. Each change needs to be recalculated, which is quite time consuming, so it is difficult to meet the demand for rapid real-time positioning of the ultra-short baselines.
In Figure 5, the sound speed profile shows a negative gradient trend. The sound velocity is basically constant within 15 m, and the sound velocity gradient changes significantly beyond 15 m. The sound velocity tends to be stable at depth about 100 m, with a slightly negative gradient after the depth of 200 m. In areas where the sound speed is stable, the error generated by the ray tracing is relatively small, and tracking with the constant sound speed ray tracing method still achieves good positioning accuracy. The main error in ray tracing is generated in the area where the sound speed changes greatly. In such an area, the equal gradient ray tracing method can better fit the actual situation.
In this paper, the constant sound velocity tracking method is combined with the equal gradient ray tracing method, and an improved sound ray tracing method is proposed. The principle of the proposed method is as follows. In the area where the sound speed is constant with depth, the constant sound speed tracking method is used for tracking, and in the area where the sound speed changes rapidly, the equal-gradient sound tracking method is used. By increasing the number of calculation iterations in the local area, the overall tracking accuracy can be improved.
The specific process is as follows. The gradient g of the sound speed profile between the ultrashort baseline system and the transponder is used as a judgment condition, and the gradient change threshold  is set according to the actual situation. When   g , the algorithm considers that the sound speed is greatly disturbed, so the equal gradient ray tracing method is used; when , the sound speed is considered stable, so the constant sound speed tracking method is used. The steps of   g In areas where the sound speed is stable, the error generated by the ray tracing is relatively small, and tracking with the constant sound speed ray tracing method still achieves good positioning accuracy. The main error in ray tracing is generated in the area where the sound speed changes greatly. In such an area, the equal gradient ray tracing method can better fit the actual situation.
In this paper, the constant sound velocity tracking method is combined with the equal gradient ray tracing method, and an improved sound ray tracing method is proposed. The principle of the proposed method is as follows. In the area where the sound speed is constant with depth, the constant sound speed tracking method is used for tracking, and in the area where the sound speed changes rapidly, the equal-gradient sound tracking method is used. By increasing the number of calculation iterations in the local area, the overall tracking accuracy can be improved.
The specific process is as follows. The gradient g of the sound speed profile between the ultra-short baseline system and the transponder is used as a judgment condition, and the gradient change threshold δ is set according to the actual situation. When g ≥ δ, the algorithm considers that the sound speed is greatly disturbed, so the equal gradient ray tracing method is used; when g < δ, the sound speed is considered stable, so the constant sound speed tracking method is used. The steps of the proposed method are as follows: 1.
Acquire the gradient of each layer of the sound speed profile: g i , i = 1, 2, · · · , n; 2.
If the consecutive levels n1 i (i = 1, 2, · · · , n) satisfy that in each layer, g i < δ, then n1 i layers will be merged into one layer n1, which will be defined as a constant velocity region, and it will be considered that the sound speed value in this region is constant; accordingly, the constant sound speed ray tracing method will be used: In Equation (24), n1 represents the number of consecutive layers in the constant-sound-speed region. Substituting the average sound speed into the constant sound speed ray tracing method, the horizontal displacement and depth of the sound wave are obtained.

3.
If consecutive levels n2 i (i = 1, 2, · · · , n) satisfy that in each layer, g i ≥ δ, then n2 i layers will be defined as an equal gradient sound-tracking area, and the equal gradient ray tracing method will be used in layer by layer in this area to calculate the horizontal displacement and depth of the sound wave in each layer.
The horizontal displacement of the corresponding layer is accumulated to obtain the total horizontal displacement at the corresponding depth, and after that, the slope is calculated. The above steps of the combined ray tracing method are shown by a flow chart, as shown in Figure 6.
Due to the random uncertainties in the acquisition process of the actual sound speed profile, during the implementation of the above algorithm, if there is a small number of consecutive layers in the equal gradient sound-tracking region but they satisfy g i ≤ δ, then, it is unreasonable to use the constant sound ray tracing method. In practical applications, there are some transient sound speed stable layers in areas with the large changes in sound speed, and the intention of the proposed algorithm is to use the equal gradient ray tracing method for the whole area where the sound speed fluctuates to a greater extent. The small number of consecutive layers also need to be contained within an equal-gradient area, and the equal gradient ray tracing method is used for this area.
Therefore, the correction of an abnormal judgment is added to the algorithm. In the process of subdividing the sound speed profile, it is stipulated that if there are obvious abnormalities in the gradient of the sound level between the contiguous layers (≤3 layers), which leads to the erroneous sound tracking, then, these layers are defined as a determination abnormality. In that case, the algorithm will re-determine the abnormal layer based on the ray tracing method used in the neighboring area to achieve the corrective effect. The horizontal displacement of the corresponding layer is accumulated to obtain the total horizontal displacement at the corresponding depth, and after that, the slope is calculated. The above steps of the combined ray tracing method are shown by a flow chart, as shown in Figure 6. Due to the random uncertainties in the acquisition process of the actual sound speed profile, during the implementation of the above algorithm, if there is a small number of consecutive layers in the equal gradient sound-tracking region but they satisfy   i g , then, it is unreasonable to use the constant sound ray tracing method. In practical applications, there are some transient sound speed stable layers in areas with the large changes in sound speed, and the intention of the proposed algorithm is to use the equal gradient ray tracing method for the whole area where the sound speed fluctuates to a greater extent. The small number of consecutive layers also need to be contained within an equal-gradient area, and the equal gradient ray tracing method is used for this area. Therefore, the correction of an abnormal judgment is added to the algorithm. In the process of subdividing the sound speed profile, it is stipulated that if there are obvious abnormalities in the gradient of the sound level between the contiguous layers ( 3  layers), which leads to the erroneous sound tracking, then, these layers are defined as a determination abnormality. In that case, the algorithm will re-determine the abnormal layer based on the ray tracing method used in the neighboring area to achieve the corrective effect.

Traditional Sound Ray Tracing Model and Slant Distance Correction
In the USBL positioning system, the slant distance from the underwater target to the ultra-short baseline system is calculated by using the signal transmitting and receiving technique. The three-dimensional position coordinates of the underwater target are obtained by using the trigonometric functions to calculate the slant distance and azimuth angle. This work mainly describes the effect of correcting the slant distance to optimize the positioning of the USBL.
In underwater acoustic positioning, the traditional method to find the slant distance is to use the propagation time T and propagation speed c of acoustic signal to obtain a direct distance l between the transponder and transmitter. The traditional method considers the sound speed is a constant value and the sound wave is set to propagate at an incident angle θ (which is also considered as the initial angle in the method below), then, it can be written that where H is the vertical height of the underwater target from the transmitter. The model of the traditional acoustic correction method is shown in Figure 7.
considered as the initial angle in the method below), then, it can be written that cos H l c T     (25) where H is the vertical height of the underwater target from the transmitter. The model of the traditional acoustic correction method is shown in Figure 7. In actual situations, a sound wave does not propagate along a straight line in the shallow sea area, and its propagation path is always bent towards the direction in which the speed of sound decreases. Therefore, the ray tracing method described above can be used to approximate the actual propagation path of a sound wave using the section of a polyline or an arc. The horizontal distance X from transceiver head to the underwater target is obtained by summing up the layer by layer using Equation (13), and the total horizontal distance X using Equation (26). In the USBL, the depth difference H between the transceiver head and the underwater target can be measured by the depth sensor, and the Pythagorean theorem can be used to obtain a more realistic slant distance between transceiver head and the underwater target lˆ, where

Experimental Results and Analysis
In this paper, the experimental data and the result analysis are all simulation analysis. In order to make the simulation environment closer to the experimental environment, the sound velocity profile used in the simulation is the measured multiple sound speed profiles. In actual situations, a sound wave does not propagate along a straight line in the shallow sea area, and its propagation path is always bent towards the direction in which the speed of sound decreases. Therefore, the ray tracing method described above can be used to approximate the actual propagation path of a sound wave using the section of a polyline or an arc. The horizontal distance X from transceiver head to the underwater target is obtained by summing up the layer by layer First, it is needed to calculate θ i (i = 1, 2, · · · , n) using Equation (13), and the total horizontal distance X using Equation (26). In the USBL, the depth difference H between the transceiver head and the underwater target can be measured by the depth sensor, and the Pythagorean theorem can be used to obtain a more realistic slant distance between transceiver head and the underwater targetl, wherê

Experimental Results and Analysis
In this paper, the experimental data and the result analysis are all simulation analysis. In order to make the simulation environment closer to the experimental environment, the sound velocity profile used in the simulation is the measured multiple sound speed profiles.

Inversion of Measured Speed Profiles
To verify the validity of the EOF, the measured velocity distribution samples were compared to the samples inversion by EOF. The measured speed profiles are shown in Figure 8. Due to equipment failures and environmental factors, there was an error in the measured sound speed profiles. As it can be seen in Figure 8, the sound speed fluctuated at the depth in the range 15-100 m and tended to be stable after the depth reached 100 m. The samples were measured at a maximum depth of about 200 m.
In the simulation experiment, we performed the first three steps of the EOF. The eight sound speed profiles was extracted for fitting, and after the inversion, the newly obtained profile was compared with the measured sound speed profile, Figure 9.
To verify the validity of the EOF, the measured velocity distribution samples were compared to the samples inversion by EOF. The measured speed profiles are shown in Figure 8. Due to equipment failures and environmental factors, there was an error in the measured sound speed profiles. As it can be seen in Figure 8, the sound speed fluctuated at the depth in the range 15-100 m and tended to be stable after the depth reached 100 m. The samples were measured at a maximum depth of about 200 m. In the simulation experiment, we performed the first three steps of the EOF. The eight sound speed profiles was extracted for fitting, and after the inversion, the newly obtained profile was compared with the measured sound speed profile, Figure 9. The mean square error (MSE) of the sound speed profile is presented in Figure 10. In Figure 10, it can be seen that the MSE was below 0.5 m/s [30].  In the simulation experiment, we performed the first three steps of the EOF. The eight sound speed profiles was extracted for fitting, and after the inversion, the newly obtained profile was compared with the measured sound speed profile, Figure 9. The mean square error (MSE) of the sound speed profile is presented in Figure 10. In Figure 10, it can be seen that the MSE was below 0.5 m/s [30]. The mean square error (MSE) of the sound speed profile is presented in Figure 10. In Figure 10, it can be seen that the MSE was below 0.5 m/s [30].
The MSE can evaluate the degree of change of the data. The smaller the value of the MSE, the better the accuracy of the prediction model describing the experimental data. At depths in the range 0-100 m, the MSE fluctuations were more frequent, and at depths larger than 100 m, the MSE was reduced to 0.2 m/s. Accordingly, it can be concluded that the EOF inversion of sound speed profile met the experimental requirements. The MSE can evaluate the degree of change of the data. The smaller the value of the MSE, the better the accuracy of the prediction model describing the experimental data. At depths in the range 0-100 m, the MSE fluctuations were more frequent, and at depths larger than 100 m, the MSE was reduced to 0.2 m/s. Accordingly, it can be concluded that the EOF inversion of sound speed profile met the experimental requirements.

Experiment with Combined Ray Tracing Method
The combined ray tracing method automatically judges whether to employ the constant sound speed ray tracing method or the equal gradient ray tracing method based on the variation in sound speed gradient.
Using the EOF inversion of the sound speed profile, the depth interval between every two layers of the sound speed profile was set to be between 2 m and 3 m; thus, the sound speed profile was divided into 100 layers, and the variation in the sound speed profile gradient i g was obtained.
In Figure 11, it can be seen that the gradient i g of the sound speed profile was stable at the 0th layer, and the sound speed was kept constant between the 0th layer and the 10th layer. Between the 10th layer and the 50th layer, significantly floated, and the sound speed changed more frequently. After the 50th layer the gradient of the sound speed tended to stabilize at zero.
We selected different speed gradient thresholds to stratify the profile. Choosing a too small or a too large threshold could cause an inadequate division of the sound speed profile, so the threshold was adjusted based on the experience and considering the actual conditions. The thresholds of 0.01, 0.06, and 0.1 were respectively used in the experiment to ensure the multiple possibilities of a regional division.
In the case of a smaller threshold, such as that of 0.01, when , the sound speed profile was divided into equal-speed-gradient areas, as in Figure 12. So, the areas were interspersed with the equal gradient ray tracing method, which resulted in the unclear distinction between the divisions of the area and increased the iterations.

Experiment with Combined Ray Tracing Method
The combined ray tracing method automatically judges whether to employ the constant sound speed ray tracing method or the equal gradient ray tracing method based on the variation in sound speed gradient.
Using the EOF inversion of the sound speed profile, the depth interval between every two layers of the sound speed profile was set to be between 2 m and 3 m; thus, the sound speed profile was divided into 100 layers, and the variation in the sound speed profile gradient g i was obtained.
In Figure 11, it can be seen that the gradient g i of the sound speed profile was stable at the 0th layer, and the sound speed was kept constant between the 0th layer and the 10th layer. Between the 10th layer and the 50th layer, g i significantly floated, and the sound speed changed more frequently. After the 50th layer the gradient of the sound speed g i tended to stabilize at zero.  We selected different speed gradient thresholds to stratify the profile. Choosing a too small or a too large threshold could cause an inadequate division of the sound speed profile, so the threshold was adjusted based on the experience and considering the actual conditions. The thresholds of 0.01, 0.06, and 0.1 were respectively used in the experiment to ensure the multiple possibilities of a regional division.
In the case of a smaller threshold, such as that of 0.01, when g i ≥ 0.01, the sound speed profile was divided into equal-speed-gradient areas, as in Figure 12. So, the areas were interspersed with the equal gradient ray tracing method, which resulted in the unclear distinction between the divisions of the area and increased the iterations.  When the threshold value was 0.06, the sound speed transform area could be distinguished well from the sound speed stable area, as shown in Figure 13. The entire sound speed profile was divided into three sections. The first section sound speed was stable, and the constant sound speed ray tracing method reduced the amount of iterations and did not produce too much error. In the second section, the equal gradient ray tracing method improved the accuracy of the complex sound speeds. In the last section, the speed of sound tended to be stable, and tracking was performed using the constant sound speed ray tracing method. By increasing the amount of calculation in the local area, the overall tracking accuracy could be improved. When the threshold value was 0.06, the sound speed transform area could be distinguished well from the sound speed stable area, as shown in Figure 13. The entire sound speed profile was divided into three sections. The first section sound speed was stable, and the constant sound speed ray tracing method reduced the amount of iterations and did not produce too much error. In the second section, the equal gradient ray tracing method improved the accuracy of the complex sound speeds. In the last section, the speed of sound tended to be stable, and tracking was performed using the constant sound speed ray tracing method. By increasing the amount of calculation in the local area, the overall tracking accuracy could be improved. When the threshold value was too large, such as that of 0.1 presented in Figure 14, in the region where the sound speed changed frequently, was greater than 0.01, and G was less than or equal to A, which led to the division in too many layers, and the calculation method was switched back and forth, resulting in inaccurate calculation results. When the threshold value was too large, such as that of 0.1 presented in Figure 14, in the region where the sound speed changed frequently, g i was greater than 0.01, and G was less than or equal to A, which led to the division in too many layers, and the calculation method was switched back and forth, resulting in inaccurate calculation results.
When the threshold value was too large, such as that of 0.1 presented in Figure 14, in the region where the sound speed changed frequently, was greater than 0.01, and G was less than or equal to A, which led to the division in too many layers, and the calculation method was switched back and forth, resulting in inaccurate calculation results. Hence, it can be concluded that when the threshold was too small or too large, the combined ray tracing method was inaccurate in judgment, and the sound speed profile was layered and disordered, resulting in errors in the sound propagation path correction. In this simulation experiment, at the threshold of 0.06, the level using constant sound speed ray tracing and the level using the equal gradient ray tracing can be better preserved, and enter the calculation of the next slant range correction. Hence, it can be concluded that when the threshold was too small or too large, the combined ray tracing method was inaccurate in judgment, and the sound speed profile was layered and disordered, resulting in errors in the sound propagation path correction. In this simulation experiment, at the threshold of 0.06, the level using constant sound speed ray tracing and the level using the equal gradient ray tracing can be better preserved, and enter the calculation of the next slant range correction.

Slant Distance Correction
In underwater acoustic positioning, the large azimuth angle of the sound ray will lead to the underwater sound waves diverge during the incident process, so the azimuth angle control under 70 • [27][28][29][30]. The experiments here are based on the environment of the simulation experiment. In the simulation experiment, the azimuth angle and elevation angle of the USBL transmitter are both set to 30 • , 40 • , 45 • . Therefore, we placed the target respectively at a depth of 10, 50, 100, 150, and 200 m from the USBL transmitter. For simulation purposes, we calculate the true slant range of the target at different depths by Equation (29) in simulation experiment and compare the true slant distance with the corrected slant range. We used Equation (25) for the traditional method iterations, and Equation (17) for the constant sound speed ray tracing method iterations, and compare the obtained results with the results of the proposed combined ray tracing method. The results are shown in Table 1.
where l true represents the slant distance of the underwater target in simulation experiment. In Table 1, the calculation iterations parameter denotes the number of times the horizontal displacement of each layer was added due to the division of the sound speed profile layers. Figure 15 shows the simulated sound ray tracing model at 45 • . It can be seen that the sound ray tracing model under the traditional method is completely a straight line. At the beginning of the constant sound ray tracing method and the combined ray tracing method, the sound ray tracing almost overlaps, but gradually becomes curved later.  Table 2, where M1 represents the constant sound ray tracing method, and M2 represents the combined ray tracing method.  As it can be seen in Table 2, the combined ray correction (method M2) exhibits an error that is always smaller than the constant sound speed ray tracing method (method M1) and using a fixed sound speed (traditional method). The worst error observed for M2 (7.4% for 40°, 10m) is approximately one-third of the error obtained by M1 and the traditional method.
Calculating the ratio error(M2)/error(M1), we can objectively find the reduction of the relative It can be seen from Figure 15 that in the sound line iteration, the model of sound line propagation under the traditional method has always been a straight line. However, the M1 and M2 methods have different degrees of bending.
The error analysis was used to evaluate the comparison of positioning error between different ray tracing methods. The error was defined as where l true represents the true slant distance, l new represents the new slant distance obtained by the three methods. The comparison of the accuracy of different ray tracing methods, calculated using Equation (30), is given in Table 2, where M1 represents the constant sound ray tracing method, and M2 represents the combined ray tracing method. As it can be seen in Table 2, the combined ray correction (method M2) exhibits an error that is always smaller than the constant sound speed ray tracing method (method M1) and using a fixed sound speed (traditional method). The worst error observed for M2 (7.4% for 40 • , 10 m) is approximately one-third of the error obtained by M1 and the traditional method.
Calculating the ratio error(M2)/error(M1), we can objectively find the reduction of the relative error when using M2 instead of M1.
From Table 3 we can see that the reduction of the relative error reduced at least 20% when using the combined ray tracing method, and the maximum reduction is about 70%. From the simulation results, we can find that the combined ray tracing method reduces the relative error of the slant distance in USBL without increasing the amount of calculation iterations, which has research significance and engineering application value.

Conclusions
The development of ocean needs the cooperation of many kinds of ocean equipment, various sensors are especially needed. The sensors work either independently or form an underwater sensor network [31][32][33][34][35] to accomplish the preset mission. Sufficient accurate position information is the primary prerequisite for achieving the goal.
To improve the accuracy of slant distance in the USBL positioning system, the ray tracing method is introduced in the underwater acoustic positioning system to eliminate the influence of the ray bending on the calculation of the slant distance for the positioning accuracy. This paper presents a combined ray tracing method that combines the constant sound speed ray tracing method with the equal gradient ray tracing method. The proposed method can adaptively determine which ray tracing method should be employed; namely, the constant sound speed ray tracing method is used in the regions where the sound velocity is stable, and the equal gradient ray tracing method is used in regions where the sound velocity changes significantly.
In the simulation experiments, the EOF was used to invert the real-time sound speed profile, which mitigated the errors in the measured sound velocity profile caused by the environmental and some other factors. After the sound speed profile was obtained, the combined ray tracing method autonomously judged which ray tracing method to use based on the gradient change in the sound speed profile.
The operational time of the proposed method is very short, so the sound source of the USBL remains relatively stationary with the transponder on the underwater target. To improve the positioning accuracy of the USBL positioning system further, it is also possible to optimize the sound ray correction and continuous positioning between the ultra-short baselines and moving sound sources, which is not presented in this work due to the limited space. Besides, we will conduct more in-depth studies on that subject in the future.