Accurate Multiple Ocean Bottom Seismometer Positioning in Shallow Water Using GNSS/Acoustic Technique

The Global Navigation Satellite System combined with acoustic technique has achieved great economic benefits in positioning of ocean bottom seismometers, with hundreds of underwater transponders attached to seismometers typically being deployed during oil exploration. The previous single transponder positioning method ignored the similar underwater environments between the transponders. Due to the refraction effect of sound, the technique usually showed poor positioning accuracy in shallow water when the incidence angles are large. In this paper, the effect of sound ray bending is analyzed based on the sound ray tracing method in shallow water, and a new piecewise incidence angle model is proposed to improve the positioning accuracy of multiple objects in order to estimate the sound ray bending correction. The parameters of the new model are divided into groups and estimated by sequential least squares method, together with all of the transponders. The observability analysis is discussed in simulation and testing experiments in the South China Sea. The results show that the newly proposed method is able to make full use of the acoustic observation data of hundreds of transponders to accurately estimate the SRB correction, which could also significantly improve the positioning accuracy of multiple transponders.


Introduction
During data acquisition using ocean bottom seismometers (OBS) in shallow water, it is usually necessary to locate a large number of submarine seismometers [1]. There are two traditional methods for locating the position of seismometers, one of which is the first break secondary positioning (FBSP) method [2,3]. The positioning accuracy of this method is generally meters (3~10 m), and higher positioning precision for a large number of seismometers within the sub-meter level is a challenge. The other method combines the Global Navigation Satellite System (GNSS) and the acoustic ranging technique, and provides important positioning measurements for underwater objects [4][5][6][7][8]. This method is able to utilize low-cost underwater transponders, which are small and lightweight, and are easily attached to the seismic cable or OBS. This technique can get the position of OBS more accurately than FBSP in shallow water, and has been proved to satisfy the requirements of oil exploration [6,7].
The Real-Time Kinematic (RTK) GNSS technique is always employed for GNSS/acoustic surveys, which achieve a 2-3 cm positioning accuracy and provide a stable absolute position reference [9]. However, the position of exploration vessel is usually far from land for oil exploration and can only be estimated with an accuracy at the decimeter level utilizing Global Navigation Satellite System (DGNSS) technology (e.g., StarFire, VeriPos and Marinestar) based on the communication data link of Inmarsat [6,7]. In contrast to the electromagnetic ranging technique, the acoustic range is the product of sound speed and travel time in the case of ignoring the sound refraction. As a development of the underwater acoustic ranging technique, time delay errors caused by system hardware are effectively controlled and eliminated. Although the speed of sound in the ocean can be directly measured by a conductivity temperature depth gauge (CTD) or sound speed profiler (SSP) gauge, the errors induced by the inaccuracy of SSP and sound refraction are the main factors affecting underwater acoustic ranging. The acoustic ranging is directly affected by the variation of sound speed along the trajectory [10,11]. When the incidence angle is greater than 85 degrees, the bending error of the sound line can reach more than 0.15 m in 100 m of water.
Quite a few studies have measured and analyzed SSP with different techniques to reduce the errors caused by sound speed and sound ray bending (SRB) [12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Previous methods, such as the epoch difference method [12] and the time model [15,17], are usually used under conditions with long observation times, where the SRB changes more significantly over time, and are not convenient during shallow water oil exploration. During OBS acoustic positioning, all the transponders only obtain a small amount of the observed data, and SRB with a high incidence angle has a great effect on positioning. The effect of an inaccurate stochastic model on the underwater positioning accuracy in shallow water has been analyzed [21], and the incidence angle segmented cosine form was used to improve positioning accuracy. This was an empirical approach, and did not analyze the function relation between the SRB and the incidence angle in detail. From the perspective of position accuracy and scientific research, it is more meaningful to explore the relationship between acoustic bending and incidence angle.
Therefore, the inaccurate sound speed and the SRB become the major factors influencing the accuracy of acoustic OBS positioning. The traditional estimated model without sound speed is usually adopted when the sound speed is inaccurate [17]. However, the OBS acoustic data of single transponders are not enough to correct for inaccurate sound speed. As hundreds of transponders hang on an underwater cable, they can provide an abundance of acoustic observations to model the acoustic refraction and estimate the SRB correction. In this research, the Ocean Bottom Cable (OBC) positioning system is introduced, and a relevant method for positioning single transponder is given. Then we analyze the characteristics of the SRB in shallow water based on the sound ray tracing (SRT) method. With the above analysis, we present a new segment incidence angle (SIA) model in shallow water. In the new method, the SRB correction can be divided into groups based on the incidence angles, and the position parameters and model parameters can be estimated together using the Sequence Least Square (SLS) method. Simulation and testing experiments in the South China Sea are used to verify and evaluate the new method.
This paper is organized as follows. Section 2 introduces the OBC measurement system in shallow water and provides details about the estimation method for survey vessel position and transponders. The SRT method, SRB correction and mathematical formulations used for the estimation of each position of seafloor transponders will be described in Section 3. This section also presents a new SIA model, and describes the new approach to positioning multiple transponders. Section 4 introduces the simulation and experiments. Finally, Section 5 presents the conclusions of the study. Figure 1 illustrates the ocean bottom cable positioning system for shallow water oil exploration. This system includes a sea floor cable with hundreds of transponders attached, a deck unit with a GNSS antenna, and a dunking transducer on a rigid pole. The global differential GNSS technology can be utilized far away from the mainland, and the level of positioning accuracy (Veripos) is better than 30 cm [26]. The relative displacement between the GNSS antennas and the transducer can be calculated using attitude angles. The survey ship sails along the survey line after the cable has been deployed, and the transducer will interrogate the sea floor transponders. A deck unit commands the transducer, which transmits and receives acoustic signals and measures the travel times from the transducer to the seafloor acoustic transponders. The number of acoustic observations for a transponder will typically be less than one hundred, and some of them have large incidence angles.

Ocean Bottom Cable Measurement System
Sensors 2019, 19, x FOR PEER REVIEW  3 of 17 deployed, and the transducer will interrogate the sea floor transponders. A deck unit commands the transducer, which transmits and receives acoustic signals and measures the travel times from the transducer to the seafloor acoustic transponders. The number of acoustic observations for a transponder will typically be less than one hundred, and some of them have large incidence angles. The SSP can be obtained by SSP gauge or CTD, and the mean speed can be calculated using the weighted or equivalent SSP method [22,26]. The sampling rate of kinematic GNSS is generally 1Hz, and the acoustic system is usually near 10 s, so the position of GNSS at the acoustic sampling time is usually obtained through the Lagrange interpolation method. The survey vessel surveys in a circle in the process of exploration and positioning to reduce the ray bending error with parallel and symmetrical observation structures. The acoustic transmission source level is 185 dB and the receiver sensitivity is 110 dB. This can be well applied in the acoustic environment of shallow water, and the ranging accuracy is better than 0.5 m at a range of 100 m [3][4][5][6].

The Positioning of the Survey Vessel
The relationship of coordinates between the GNSS antenna and the reference point can be expressed by [5]: where ( ) , , gps gps gps x y z represents the coordinates of the GNSS antenna, and h , p and r represent the attitude measurements of heading, pitch and rolling angles, respectively. This can be measured by an electrical gyrocompass and an attitude sensor, or by GNSS attitude determination.
( ) , , x y z Δ Δ Δ represents the baseline between the GNSS antenna and the reference point, which is defined as the origin of the body-fixed coordinate frame. ( ) , , x y z represents the coordinates of the reference point.
The position of transponders can also be computed using the following formulation: The SSP can be obtained by SSP gauge or CTD, and the mean speed can be calculated using the weighted or equivalent SSP method [22,26]. The sampling rate of kinematic GNSS is generally 1Hz, and the acoustic system is usually near 10 s, so the position of GNSS at the acoustic sampling time is usually obtained through the Lagrange interpolation method. The survey vessel surveys in a circle in the process of exploration and positioning to reduce the ray bending error with parallel and symmetrical observation structures. The acoustic transmission source level is 185 dB and the receiver sensitivity is 110 dB. This can be well applied in the acoustic environment of shallow water, and the ranging accuracy is better than 0.5 m at a range of 100 m [3][4][5][6].

The Positioning of the Survey Vessel
The relationship of coordinates between the GNSS antenna and the reference point can be expressed by [5]: where x gps , y gps , z gps represents the coordinates of the GNSS antenna, and h, p and r represent the attitude measurements of heading, pitch and rolling angles, respectively. This can be measured by an electrical gyrocompass and an attitude sensor, or by GNSS attitude determination. (∆x 1 , ∆y 1 , ∆z 1 ) represents the baseline between the GNSS antenna and the reference point, which is defined as the origin of the body-fixed coordinate frame. (x, y, z) represents the coordinates of the reference point.
The position of transponders can also be computed using the following formulation: where (x, y, z) T denotes the coordinates of the transducer in the navigation coordinate system with the transformation of ship attitude, and (x, y, z) T 0 represents the original coordinates of the transducer in the ship coordinate system and can be directly measured by a total station instrument while the ship is in dry-dock. The main function of the total station instrument is to survey the relative coordinates.

Calculation of Transponder Positions
Generally, the time delay can be measured and corrected by the instrument manufacturer, and the observation equation can be expressed as [12] where ρ(i, k) represents the acoustic ranging at time i. k represents the transponder mark, and f (X T (k), X(i, k)) represents the geometric distance between the transducer and the transponder. X and X T represent the coordinates of the transducer and the transponder, respectively. δρ v (i, k) represents the systematic error due to the refraction of sound rays, and ε(i, k) represents the Gaussian noise. The mean sound velocity (MSV) can be measured and calculated, or determined based on experience, and the linearized observation equation is as follows: where ρ o (i, k) = C e · t(i, k), C e represents the MSV, and t(i, k) represents the travel time; is the Jacobian matrix of the measurement equation, and represents the first partial derivatives with respect to X, and ε p represents the GNSS antenna error at the position X. If the impact of sound speed variation is small, the δρ v can be ignored. dx represents the position correction of the transponder. The position correction can be estimated by the least squares method. This estimation method will subsequently be represented by LS1 in the experiment. In shallow water, the SRB changes sharply with the incidence angle for one transponder. When the incidence angle is greater than 65 degrees, the SRB curvature is obvious. In practice, the observation data are selected according to the cut incidence angle, which can generally reduce the influence of the bending error. This method will later be represented by LS2 in the experiment.
If the SSP is not measured, we can give an initial value of sound speed and estimate its correction. The observation equation can be expressed as [17] where ρ o (i, k) = C o · t(i, k), and C o can be set to 1500 m/s. This method is usually used in the absence of correct MSV, and it will later be represented by LS3 in the experiment.

The Effect of SRB in Shallow Water
The SRT method with constant gradient is introduced in this section. If the sound speed changes little or uniformly, then we can consider the sound speed gradient to be a constant. As shown in  Figure 2, C j represents the sound speed in water layer j, and z j represents the depth of the water. The travel of sound in water follows Snell's law [10,11] sin θ j /C j = P (8) where θ j represents incidence angle, and P is a constant. In each layer, we assume that the medium changes uniformly. In other words, the sound speed has a constant gradient, and the track of the ray is a continuous arc. The radius R j of the arc can be expressed as [27]: where g i represents the sound speed gradient. The horizontal distance of the sound ray x j can be computed as follows where ∆z j represents the length of layer j. In general, the interval between layers is wide if the variation of sound velocity between layers is small. The length of arc S j in layer j can be computed by [28]

The Effect of SRB in Shallow Water
The SRT method with constant gradient is introduced in this section. If the sound speed changes little or uniformly, then we can consider the sound speed gradient to be a constant. As shown in Figure 2, j C represents the sound speed in water layer j , and j z represents the depth of the water.
The travel of sound in water follows Snell's law [10,11] sin / j j where j θ represents incidence angle, and P is a constant. In each layer, we assume that the medium changes uniformly. In other words, the sound speed has a constant gradient, and the track of the ray is a continuous arc. The radius j R of the arc can be expressed as [27]: where i g represents the sound speed gradient. The horizontal distance of the sound ray j x can be computed as follows where j z Δ represents the length of layer j . In general, the interval between layers is wide if the variation of sound velocity between layers is small. The length of arc j S in layer j can be computed by [28] ( ) Thus, the travel time is estimated by In OBS underwater acoustic positioning, a fixed weighted mean sound velocity (WMSV) is usually used in the calculation, and the WMSV can be calculated by [29]  Thus, the travel time is estimated by In OBS underwater acoustic positioning, a fixed weighted mean sound velocity (WMSV) is usually used in the calculation, and the WMSV can be calculated by [29] where H represents the depth of water; w j is the weight of each water column related to the incidence angle, sound speed structure and other factors. If w j = 1, the C w represents the MSV. The MSV assumes a single sound speed value for the entire localization area. In fact, the calculated MSV is different due to the incidence angles and sound ray, and can cause the SRB error. The SRB correction δρ v can be modeled and absorbed by the MSV in many forms, such as a constant, 1st-or 2nd-degree polynomial with a time series [15] as The SRB correction δρ v can also be modeled directly [17] as where a is the coefficient to be estimated, and t i represents the travel time. t ma represents the minimum travel time among the measured epochs. λ ma represents the incidence angle according to the calculation of the position of the transducer. The above methods are usually used under conditions where there is a long observation time, and the SRB varies with time more significantly. Based on the above methods, we can analyze the SRB in shallow water, which helps us to find some useful models to describe it quantitatively.
The relationship among the horizontal distance (HD), geometrical distance (RD) and sound ray trace (SD) is shown in Figure 3. According to the analysis of the above, the acoustic ranging is usually treated as a straight line, while the real sound track is a curve in inhomogeneous water. In this research, we ignore the complexity of the top surface sound speed, and consider the analysis of the steady acoustic environment in small areas within a short period, as the OBS measurement area is only a few hundred square meters, and a single operation period is not usually more than an hour. We temporarily ignore the complexity of those short sound speed changes, and combine all the observation data to estimate the SRB error in this region. We only consider the commonness between these observation data; for example, under the same incidence, the SIB error is approximately equivalent. We use these characteristics and combine all the observation data to estimate the bending error of the acoustic line in this region.
where H represents the depth of water; j w is the weight of each water column related to the incidence angle, sound speed structure and other factors. If The MSV assumes a single sound speed value for the entire localization area. In fact, the calculated MSV is different due to the incidence angles and sound ray, and can cause the SRB error.
The SRB correction v δρ can be modeled and absorbed by the MSV in many forms, such as a constant, 1st-or 2nd-degree polynomial with a time series [15] as The SRB correction v δρ can also be modeled directly [17] as where a is the coefficient to be estimated, and i t represents the travel time. ma t represents the minimum travel time among the measured epochs. ma λ represents the incidence angle according to the calculation of the position of the transducer. The above methods are usually used under conditions where there is a long observation time, and the SRB varies with time more significantly. Based on the above methods, we can analyze the SRB in shallow water, which helps us to find some useful models to describe it quantitatively.
The relationship among the horizontal distance (HD), geometrical distance (RD) and sound ray trace (SD) is shown in Figure 3. According to the analysis of the above, the acoustic ranging is usually treated as a straight line, while the real sound track is a curve in inhomogeneous water. In this research, we ignore the complexity of the top surface sound speed, and consider the analysis of the steady acoustic environment in small areas within a short period, as the OBS measurement area is only a few hundred square meters, and a single operation period is not usually more than an hour. We temporarily ignore the complexity of those short sound speed changes, and combine all the observation data to estimate the SRB error in this region. We only consider the commonness between these observation data; for example, under the same incidence, the SIB error is approximately equivalent. We use these characteristics and combine all the observation data to estimate the bending error of the acoustic line in this region. Based on the above equation, the SRB can be simulated using the measured SSP. Assuming that the coordinates of the transponder and the transducer are known, the beam incidence angle and Based on the above equation, the SRB can be simulated using the measured SSP. Assuming that the coordinates of the transponder and the transducer are known, the beam incidence angle and travel time, as unknown parameters, can be solved by Newton iteration or other secant methods using the SSP.
As shown in Figure 4, we can get some information from the above figures. The HD, RD, SD − RD and SD − C · t increase significantly with the increase of incidence angle and water depth, and show little change when the incidence angle is smaller than 60 degrees in shallow water. The effect of SRB can reach 0.18 m at a depth of 100 m. The largest HD can reach 450 m when the incidence angle is 80 degrees, and this determines the measured horizontal range.
As shown in Figure 4, we can get some information from the above figures. The HD , RD , SD RD − and SD C t − ⋅ increase significantly with the increase of incidence angle and water depth, and show little change when the incidence angle is smaller than 60 degrees in shallow water. The effect of SRB can reach 0.18 m at a depth of 100 m. The largest HD can reach 450 m when the incidence angle is 80 degrees, and this determines the measured horizontal range. From the above analysis (Section 3.1), the SRB correction can be modeled approximately according to incidence angle at the same depth.

3.2.A Segmented Incidence Angle (SIA) Model
In this research, we assume that the SSP is inaccurate or the CTD has not been calibrated; thereby, the MSV and the SRB become the major factors affecting the accuracy of the acoustic positioning. This problem is also the most common for processing acoustic ranging data. In order to build an accurate fitting model, the parameters of SRB correction are divided into different groups by the incidence angles, and we can find a model through many simulation experiments. A new model is proposed to estimate the SRB correction according to the incidence angles where i θ represents the incidence angle, and M w represents the coefficient of the model. M b represents the constant deviation due to inaccurate MSV. N represents the cosine of the incidence angle to the N. When the incidence angle is less than 60°, the SRB correction can be neglected according to the analysis in the previous section. As shown in Figure 5, the new SIA model can fit well with the SRB correction by selecting the parameters. The new model is used to fit the SRB correction, and N is selected as 3 or 4, corresponding to segments 60-70° or 70-80°, depending on the rising trend. From the above analysis (Section 3.1), the SRB correction can be modeled approximately according to incidence angle at the same depth.

A Segmented Incidence Angle (SIA) Model
In this research, we assume that the SSP is inaccurate or the CTD has not been calibrated; thereby, the MSV and the SRB become the major factors affecting the accuracy of the acoustic positioning. This problem is also the most common for processing acoustic ranging data. In order to build an accurate fitting model, the parameters of SRB correction are divided into different groups by the incidence angles, and we can find a model through many simulation experiments. A new model is proposed to estimate the SRB correction according to the incidence angles where θ i represents the incidence angle, and w M represents the coefficient of the model. b M represents the constant deviation due to inaccurate MSV. N represents the cosine of the incidence angle to the N. When the incidence angle is less than 60 • , the SRB correction can be neglected according to the analysis in the previous section. As shown in Figure 5, the new SIA model can fit well with the SRB correction by selecting the parameters. The new model is used to fit the SRB correction, and N is selected as 3 or 4, corresponding to segments 60-70 • or 70-80 • , depending on the rising trend.
In shallow water OBS acoustic positioning, the area of multiple transponder positioning is normally small (one square kilometer). All transponders are located on the shallow and flat seabed (most shallow sea exploration terrain meets this condition); then, the water depth of the area can even be used as a fixed value. Assuming that the effect on SSP caused by the non-barotropic tidal flow or internal gravitational wave in a short time is small, the observation environment will be similar for all transponders within a short time. There are reasons to believe that the same incidence angles of acoustic ranging have the same SRB, which increases with the increase of the incidence angle. Naturally, the key issue is to select correct parameters, and the next section focuses on this problem. In shallow water OBS acoustic positioning, the area of multiple transponder positioning is normally small (one square kilometer). All transponders are located on the shallow and flat seabed (most shallow sea exploration terrain meets this condition); then, the water depth of the area can even be used as a fixed value. Assuming that the effect on SSP caused by the non-barotropic tidal flow or internal gravitational wave in a short time is small, the observation environment will be similar for all transponders within a short time. There are reasons to believe that the same incidence angles of acoustic ranging have the same SRB, which increases with the increase of the incidence angle. Naturally, the key issue is to select correct parameters, and the next section focuses on this problem. A high-precision underwater positioning algorithm for multiple targets with regard to acoustic line bending error is presented.

Calculation Method of Multiple Transponders with Sequential Least Square
For a single transponder, it is easy to cause ill-posed problems due to the introduction of many model parameters. From the above analysis, if the observations of different transponders have the same incidence angle, the same SRB correction can be estimated together. A large number of observations will place a huge burden on the computer. Based on the Sequential Least Square (SLS) method and the matrix orthogonal principle, a convenient solution is given in this research.
We first use the LS1 method or the LS2 method to calculate the initial transponder coordinates, and the observed values can be categorized according to the threshold values and incidence angles.
Then, the unknown parameters 1 ( ) ( ), ( ), ( ), ,..., , and the new observation equation can be expressed as where Δ represents the sum of random errors and modeling errors, and [ ]   Figure 5. The SRB correction using the SIA model.

Calculation Method of Multiple Transponders with Sequential Least Square
For a single transponder, it is easy to cause ill-posed problems due to the introduction of many model parameters. From the above analysis, if the observations of different transponders have the same incidence angle, the same SRB correction can be estimated together. A large number of observations will place a huge burden on the computer. Based on the Sequential Least Square (SLS) method and the matrix orthogonal principle, a convenient solution is given in this research.
We first use the LS1 method or the LS2 method to calculate the initial transponder coordinates, and the observed values can be categorized according to the threshold values and incidence angles. Then, where ∆ represents the sum of random errors and modeling errors, and B a = [B 1 , B 2 ]. L a represents the linearized observation.
and B 1 represents the Jacobian matrix after linearization. B 2 represents the coefficients of sound correction, and p 1 + p 2 + · · · + p m = n 1 = n − n 2 . m represents the number of groups. n 1 is the number of observations where the incidence angle is larger than 65 degrees, and n is the number of all transponder observations. According to the matrix orthogonal principle, the simplification of (17) is as follows where , and I represents an identity matrix. ∑ l represents the covariance of observation L a , and dc represents the SIA parameters.
Then the SLS method can be used to estimate the SIA parameters [30].
where k represents the transponder ID, and A k = (R G1 B 2 ) k .X k =dc, and its variance-covariance matrix can be expressed bŷ As shown in Figure 6, for the final transponder, SIA parameters and variance-covariance matrix can be calculated by solving Equations (21) and (22). When the final new model parameters of the SRB correction are solved, they are brought into Equation (17). Thus, the coordinates and incidence angles can be estimated at this time using the LS1 method again. The new results are brought into the sequential least squares algorithm and iterated until there are few differences from the previous results.  Figure 6. Flowchart of calculation method of multiple transponders using the SLS method.

Case 1: Simulation
This experiment is designed to verify the accuracy attainable using the new methods for multiple transponders. As shown in Figure 7a  If the effect on SSP caused by the non-barotropic tidal flow or internal gravitational wave in a short time cannot be ignored, the new model parameters can be remodeled and estimated for sections of time. Only some new modeling parameters are added and estimated here. Certainly, we could also use the other models (see Equation (14) or Equation (15)) to reduce the impact of this change on positioning. Real-time and high-precision SSP is difficult to obtain in OBC positioning, and it is usually measured before or after positioning. The method proposed in this research is another way of improving its accuracy based on estimating the parameters of the acoustic bending model.

Case 1: Simulation
This experiment is designed to verify the accuracy attainable using the new methods for multiple transponders. As shown in Figure 7a Figure 7b, since the environments are similar within a short period (for example, 1 h), we will obtain the round-trip times using Snell's law of refraction with 12 months' SSP of the same region. The positioning errors of the transducers are 10 cm in the horizontal direction and 30 cm in the vertical direction, and the random errors in travel time are 10 µs. The incidence angles range from 40 • to 80 • , which correspond to the underwater acoustic data in the process of offshore oil exploration. The speed of the ship is 3 knots, and the sampling period of the transducers are 10 s and 45 s, which correspond to good and bad observation conditions. At the same time, in order to approximate the actual situation, the distribution of the transponder was arranged in the form of cosine curve in the bad case to simulate the asymmetry of observation.  Figure 6. Flowchart of calculation method of multiple transponders using the SLS method.

Case 1: Simulation
This experiment is designed to verify the accuracy attainable using the new methods for multiple transponders. As shown in Figure 7a Figure 7b, since the environments are similar within a short period (for example, 1 h), we will obtain the round-trip times using Snell's law of refraction with 12 months' SSP of the same region. The positioning errors of the transducers are 10 cm in the horizontal direction and 30 cm in the vertical direction, and the random errors in travel time are 10 µs. The incidence angles range from 40° to 80°, which correspond to the underwater acoustic data in the process of offshore oil exploration. The speed of the ship is 3 knots, and the sampling period of the transducers are 10 s and 45 s, which correspond to good and bad observation conditions. At the same time, in order to approximate the actual situation, the distribution of the transponder was arranged in the form of cosine curve in the bad case to simulate the asymmetry of observation. As shown in Figure 8a, the SRB correction is related to the incidence angles, depth and time, and increases sharply with the increase of incidence angle and depth. From the histogram of the incidence angle as shown in Figure 8b, we can find that there are many incidence angles greater than 65°, and the problem of SRB must be taken into consideration. As shown in Figure 8a, the SRB correction is related to the incidence angles, depth and time, and increases sharply with the increase of incidence angle and depth. From the histogram of the incidence angle as shown in Figure 8b, we can find that there are many incidence angles greater than 65 • , and the problem of SRB must be taken into consideration. In particular, the SLS algorithm reordered the data according to the incidence angle for grouping calculation, so the epoch (Figure 9) does not have the significance of the time sequence. As shown in Figure 9, the SRB correction for 20 transponders can be calculated using the SLS method. (.) norm represents modulus operation, n represents the number of transponders. Table 1 shows the MPB of horizontal and vertical coordinates calculated using the above methods under good observation conditions. Each table includes three parts, corresponding to water depths of 100 m, 200 m and 300 m. As a result, even though the errors of SRB in the original ranging measurements are very large, the new method (SLS) is more effective at reducing the effect of the errors of SRB. We find that the new method can improve the results with centimeter-level accuracy in the horizontal and vertical components at depths of 100 m, 200 m and 300 m, respectively. The LS3 method can cancel the effect of SRB by estimating the sound speed, and the average positioning accuracy of transponders is higher than 10 cm in the horizontal direction. However, unlike the other methods, the new method is able to calculate the vertical components with centimeter-level accuracy as well, and the average accuracy of the vertical components is better than 4 cm. The actual observation conditions are worse than the above experiments; for example, the number of observations is insufficient and the observation structure is asymmetric. To investigate the effect of sampling period (good or bad observational conditions) on multiple transponders by different methods, the sampling rate is changed from 15 s to 45 s. In particular, the SLS algorithm reordered the data according to the incidence angle for grouping calculation, so the epoch (Figure 9) does not have the significance of the time sequence. As shown in Figure 9, the SRB correction for 20 transponders can be calculated using the SLS method. The true SRB correction can reach up to 0.35 m at a depth of 100 m. Comparing the estimated correction with the true SRB, the fit validity of the SLS method for estimating SRB correction can be proved. The observation epochs of the 20 transponders are counted and the position bias of the three directions is calculated. The effectiveness of the new algorithm is also verified by the final positioning results. The mean position bias (MPB) is used to evaluate the accuracy of the different methods, and is defined by where (x,ŷ,ẑ) represent the estimated transponder coordinates and (x, y, z) are the true transponder coordinates. norm(.) represents modulus operation, n represents the number of transponders. Table 1 shows the MPB of horizontal and vertical coordinates calculated using the above methods under good observation conditions. Each table includes three parts, corresponding to water depths of 100 m, 200 m and 300 m. As a result, even though the errors of SRB in the original ranging measurements are very large, the new method (SLS) is more effective at reducing the effect of the errors of SRB. We find that the new method can improve the results with centimeter-level accuracy in the horizontal and vertical components at depths of 100 m, 200 m and 300 m, respectively. The LS3 method can cancel the effect of SRB by estimating the sound speed, and the average positioning accuracy of transponders is higher than 10 cm in the horizontal direction. However, unlike the other methods, the new method is able to calculate the vertical components with centimeter-level accuracy as well, and the average accuracy of the vertical components is better than 4 cm. The actual observation conditions are worse than the above experiments; for example, the number of observations is insufficient and the observation structure is asymmetric. To investigate the effect of sampling period (good or bad observational conditions) on multiple transponders by different methods, the sampling rate is changed from 15 s to 45 s. Sensors 2019, 19, x FOR PEER REVIEW 12 of 17 Figure 9. The correction of SRB by the SLS method and the value of SRB (depth: 100 m; month: January). We find that the accuracy of the horizontal and vertical components tends to decrease with the decrease in the sample rate under bad conditions (see Table 2). In this case, the SRB correction is difficult to eliminate using a symmetrical structure, but the new SLS method can also improve the accuracy of the positioning. When the depth of the transponders reaches 300 m, the average accuracy of the horizontal components is better than 10 cm.
For one transponder (depth = 300), the number of observations is about 1600/20 = 80, and most of the data have large incidence angles (see Figure 8b). To analyze the effect of observation structure on positioning accuracy, the mean position dilution of precision (MPDOP) (this is a measure of X, Y, Z position geometry) can be calculated by  We find that the accuracy of the horizontal and vertical components tends to decrease with the decrease in the sample rate under bad conditions (see Table 2). In this case, the SRB correction is difficult to eliminate using a symmetrical structure, but the new SLS method can also improve the accuracy of the positioning. When the depth of the transponders reaches 300 m, the average accuracy of the horizontal components is better than 10 cm. For one transponder (depth = 300), the number of observations is about 1600/20 = 80, and most of the data have large incidence angles (see Figure 8b). To analyze the effect of observation structure on positioning accuracy, the mean position dilution of precision (MPDOP) (this is a measure of X, Y, Z position geometry) can be calculated by where q 2 x i , q 2 y i and q 2 z i are the elements of the covariance matrix of the estimated parameters, and n represents the number of transponders. The MPDOP (depth = 100 m, 200 m and 300 m) is as follows.
As shown in Table 3, LS1 is a more robust algorithm than the other three methods, and the MPDOP is minimal. The LS1 method is better than the LS2 method for taking full advantage of all the data, and the system is more stable. The stability of the SLS method is reduced with increasing model parameters. Although this may lead to a certain amount of instability in the solution, the actual positioning accuracy is improved with the more refined solution of SRB correction.

Case 2: Experiment in the South China Sea
As shown in Figure 10, the ship is equipped with GNSS, a transducer, a transponder SSP, and so on. The experiment area was on the east of Hainan Island in the South China Sea (Figure 11). The position of the ship was measured by a kinematic VeriPos DGNSS system. Once the transponders had been sunk in the sea, the exploration ship surveyed around them twice. The ship's heading was recorded with an electric gyrocompass, and four GNSS antennas fixed around the ship were used to measure the ship's roll and pitch [31], and the positioning of the transceiver on the bottom of the ship can be calculated using ship's attitude data and DGNSS measurements. The mean sound speed measured by SSP is 1534 m/s. The water depth at the experiment site was approximately 100 m, and the acoustic ranging system included a transducer and 30 transponders. The acoustic positioning system included two types, an OBC acoustic positioning system developed by UK Sonardyne Company (Yateley, UK), and the China BPS acoustic positioning system; the sea trial was roughly 40 min. In this experiment, some BPS transponders had an electrical failure. Two trajectories were completed successively in the same day. For the first trail, the OBC acoustic system and the BPS acoustic positioning system were tested together. For the second trail, only the OBC acoustic system was tested. As shown in Table 3, LS1 is a more robust algorithm than the other three methods, and the MPDOP is minimal. The LS1 method is better than the LS2 method for taking full advantage of all the data, and the system is more stable. The stability of the SLS method is reduced with increasing model parameters. Although this may lead to a certain amount of instability in the solution, the actual positioning accuracy is improved with the more refined solution of SRB correction.

Case 2: Experiment in the South China Sea
As shown in Figure 10, the ship is equipped with GNSS, a transducer, a transponder SSP, and so on. The experiment area was on the east of Hainan Island in the South China Sea ( Figure 11). The position of the ship was measured by a kinematic VeriPos DGNSS system. Once the transponders had been sunk in the sea, the exploration ship surveyed around them twice. The ship's heading was recorded with an electric gyrocompass, and four GNSS antennas fixed around the ship were used to measure the ship's roll and pitch [31], and the positioning of the transceiver on the bottom of the ship can be calculated using ship's attitude data and DGNSS measurements. The mean sound speed measured by SSP is 1534 m/s. The water depth at the experiment site was approximately 100 m, and the acoustic ranging system included a transducer and 30 transponders. The acoustic positioning system included two types, an OBC acoustic positioning system developed by UK Sonardyne Company (Yateley, UK), and the China BPS acoustic positioning system; the sea trial was roughly 40 min. In this experiment, some BPS transponders had an electrical failure. Two trajectories were completed successively in the same day. For the first trail, the OBC acoustic system and the BPS acoustic positioning system were tested together. For the second trail, only the OBC acoustic system was tested.  The observation number of each transponder is shown in Figure 12a, and the average number of acoustic measurements of a transponder is about 40. No. 1 and No. 3 transponders responded so frequently because the ship began sailing with a time measurement near the first transponder for a long time. The dotted line of incidence angles with epochs is shown in Figure 12b. Nearly half of the observation data have an incidence angle greater than 65 degrees, and abandoning these observations could cause serious ill-posed problems and damage the positioning accuracy. As the number of observations of each transponder is not enough, acoustic observations with high incidence angle are necessary. To estimate the absolute position of the seafloor transponders, we evaluate the positioning accuracy according to the mutual deviation of the two kinds of transponders. The positioning accuracy will increase, while mutual deviation decreases. Two different mean position biases (MPB1 and MPB2), using above approaches, are defined as is calculated using the OBC acoustic system only.
There are many outliers in the BPS observations, the BPS positioning results are not satisfactory, and only five transponders were chosen for investigation. Table 4 provides the absolute deviation of the four methods using the BPS acoustic system and the OBC acoustic system. Inaccurate sound The observation number of each transponder is shown in Figure 12a, and the average number of acoustic measurements of a transponder is about 40. No. 1 and No. 3 transponders responded so frequently because the ship began sailing with a time measurement near the first transponder for a long time. The dotted line of incidence angles with epochs is shown in Figure 12b. Nearly half of the observation data have an incidence angle greater than 65 degrees, and abandoning these observations could cause serious ill-posed problems and damage the positioning accuracy. As the number of observations of each transponder is not enough, acoustic observations with high incidence angle are necessary. The observation number of each transponder is shown in Figure 12a, and the average number of acoustic measurements of a transponder is about 40. No. 1 and No. 3 transponders responded so frequently because the ship began sailing with a time measurement near the first transponder for a long time. The dotted line of incidence angles with epochs is shown in Figure 12b. Nearly half of the observation data have an incidence angle greater than 65 degrees, and abandoning these observations could cause serious ill-posed problems and damage the positioning accuracy. As the number of observations of each transponder is not enough, acoustic observations with high incidence angle are necessary. To estimate the absolute position of the seafloor transponders, we evaluate the positioning accuracy according to the mutual deviation of the two kinds of transponders. The positioning accuracy will increase, while mutual deviation decreases. Two different mean position biases (MPB1 and MPB2), using above approaches, are defined as is calculated using the OBC acoustic system only.
There are many outliers in the BPS observations, the BPS positioning results are not satisfactory, and only five transponders were chosen for investigation. Table 4 provides the absolute deviation of the four methods using the BPS acoustic system and the OBC acoustic system. Inaccurate sound To estimate the absolute position of the seafloor transponders, we evaluate the positioning accuracy according to the mutual deviation of the two kinds of transponders. The positioning accuracy will increase, while mutual deviation decreases. Two different mean position biases (MPB1 and MPB2), using above approaches, are defined as norm((x,ŷ,ẑ) 1i − (x,ŷ,ẑ) 2i )/n (28) where (x,ŷ,ẑ) are the estimated transponder coordinates and (x, y, z) are the true transponder coordinates. MPB1 represents the first trail result using the OBC and BPS systems, and the MRMS2 is calculated using the OBC acoustic system only.
There are many outliers in the BPS observations, the BPS positioning results are not satisfactory, and only five transponders were chosen for investigation. Table 4 provides the absolute deviation of the four methods using the BPS acoustic system and the OBC acoustic system. Inaccurate sound speed measurements also affect the final positioning accuracy. The accuracy of the SLS is better than the other three methods, and the MPB1 in the horizontal direction is 0.54 m. Two repeated observations of MPB2 on 30 transponders using the four methods can also show the performance of different algorithms. Table 5 provides the MPB2 of the LS2 and LS3 methods with different cut angles. As the cut angle increases, the positioning accuracy of the LS2 solution is improved, with the positioning accuracy reaching as high as 0.52 m and 0.64 m in the horizontal and vertical directions, respectively, when the cut angle was 75 degrees. Since a transponder has fewer observations, the accuracy of LS3 was lower than that of LS2, and the positioning accuracy was 0.64 m. In Table 6, we have listed the horizontal and vertical MPB2 using the four methods. It can be seen that the SLS method can also improve the positioning accuracy more significantly than the other methods. The performance difference between SLS and the other methods is also clearly noticeable in favor of the SLS method, which can also be seen in Table 4. Experimental results show that the MPB2 of SLS method in the horizontal and vertical directions, respectively, are about 30 cm and 25 cm.

Conclusions
The SRB error is an important factor influencing underwater positioning. In this study, the effect of SRB has been taken into account, and we have constructed a new SIA model to estimate the SRB correction. The new method differs from the three conventional methods, and all acoustic data are fully utilized to estimate the parameters of the SRB model. Although high-incidence data will introduce SRB errors, these observations are highly conducive to the stability of the model. In the acoustic data procedure, we adopt a new method that yields beneficial results.
Simulation and experimental tests were carried out in order to evaluate the performance of the new method. The simulation results show that the average horizontal positioning accuracies of the new method at different depths are all better than 10 cm, even in the case of a poor positioning environment. Without being restricted to strict symmetrical observations, the other advantage of the new method is that it can effectively improve the accuracy in the vertical direction, and the vertical positioning accuracy can be substantially improved, from 0.8 m in the LS scheme to 0.03 m. The performance of the new method was also validated by experiments in the South China Sea, and the results show that the positioning accuracy can be substantially improved from 0.5 m in the LS scheme to 0.3, based on