A Combined Measurement Method for the Seafloor Positioning, Navigation, and Timing Network

Abstract

The idea of constructing the “GNSS-like” seafloor geodetic network for underwater positioning, navigation, and timing (PNT) has been proposed by many countries. Based on this idea, this paper introduces the principle of the seafloor PNT network and provides a combined measurement method, including the absolute positioning and the relative positioning. Experimental results show that the positioning difference between the proposed method and circle-sailing positioning is approximately 10 cm, and the observation efficiency is higher than the existing measurement in the seafloor PNT network. In addition, a model is derived to determine the optimal configuration of the unit network (the basic component of the seafloor PNT network), considering the ranging capability of acoustic beacons, which is helpful to balance the side length of the unit network and the number of observation sessions in the relative positioning. Finally, a chain coordinate transfer strategy in the whole seafloor PNT network is proposed, and the positions where the absolute positioning should be carried out in the whole network are derived based on the conditional adjustment model. Based on this strategy, the design schemes of a seafloor PNT network with centimeter, decimeter, and meter positioning accuracy when the acoustic velocity measuring accuracy is 0.02 m/s and the time measuring accuracy is 10⁻⁵ s, are given in the experimental section.

1. Introduction

The marine geodetic network is the reference frame for all marine activities. It consists of buoy datum stations on the ocean surface, mobile datum stations in the water, and the seafloor geodetic network [1]. The seafloor geodetic network is the most important component of the marine geodetic network owing to its good stability and high positioning accuracy, which plays a very important role in ocean engineering [2]. In recent years, coastal states have paid more and more attention to the construction of seafloor geodetic networks. Some countries or organizations, such as the United States, Japan, China, and the European Union, have established seafloor geodetic networks in some waters, which mainly serve marine environment monitoring and seafloor crustal deformation observation [3,4,5,6]. The seafloor geodetic network for observing seafloor crustal deformation has three main forms, including the direct-path acoustic network, the indirect-path acoustic network, and the Global Navigation Satellite System-Acoustic (GNSS-A) network [7]. In the direct-path acoustic network (Figure 1a), two precision acoustic beacons are placed on faults, and repeated interrogations over months to years measure the fault displacement. To overcome the influence of the pressure on the propagation of acoustic waves in the direct-path network, the indirect-path network was invented (Figure 1b), in which at least three seafloor acoustic beacons are deployed, and a floating unit with a depth sensor is added several hundred meters above the beacons [8]. To overcome the effect of seawater stratification on acoustic ranging, researchers have further improved the arrangement of acoustic beacons on the seafloor; three or more beacons are evenly spaced around the perimeter of a circle with a radius of the water depth, which is called the GNSS-A network (Figure 1c). In the GNSS-A network, the positioning of the seafloor acoustic beacon is realized by combining the coordinates of the onboard transducer provided by the GNSS system with the geometric distance between the onboard transducer and the seafloor beacon [9,10,11].

With the increase in marine activities, the application scope of the seafloor geodetic network has been expanded, such as in underwater positioning, navigation, and timing (PNT) for survey ships on the surface and underwater vehicles. In 2015, the American Defense Advanced Research Projects Agency (DARPA) announced the development of the Positioning System for Deep Ocean Navigation (POSYDON), which aims to provide omnipresent, robust positioning for ships and submarines across ocean basins by establishing the “GNSS-like” underwater navigation system, that is, the long baseline hydroacoustic positioning network (called the seafloor PNT network for short) using low-frequency sonar [12]. In 2016, China launched the national key research and development project named Marine Geodetic Datum and Marine Navigation and Positioning Technology, which aims to construct the large-scale seafloor PNT network [13,14]. The idea of underwater PNT is as shown in Figure 1d. The network is achieved by first deploying sonar beacons (similar to satellites) on the seabed to form a seafloor geodetic network, and determining the coordinates of all seafloor sonar beacons. Vehicles can then obtain their coordinates by range intersection positioning using acoustic signals (similar to satellite signals) from the seabed sonar beacons with known coordinates. However, the limited coverage area and specialized purpose make it difficult to serve large-scale underwater PNT for the direct-path acoustic network, indirect-path acoustic network, and GNSS-A network. Therefore, how to construct the underwater PNT network is a new question to discuss.

Similar to the terrestrial geodetic network, a large-scale seafloor PNT network is generally composed of several unit networks, such as the white triangles and squares in Figure 1d. Therefore, it is necessary to formulate a scheme for the optimization of the whole network before deploying sonar beacons from the cost, reliability, and precision views. Up to now, many scholars have carried out research on the optimization of the geodetic network. In 1974, Grafarend proposed four orders for designing the geodesic network, including designing the best datum and coordinate system of the network, designing the best configuration of the network, designing the suitable precision of the observations and types of observations, and designing the expansion or densification of the network [15]. Based on the four orders, in 1985, Koch mentioned the optimum configuration problem of the geodetic network first and used the quadratic programming and linear complementary problem to obtain the position changes or configuration of the network in an optimal way [16]. In 2007, Eshagh and Kiamehr put forward a strategy for the optimum design of the geodetic network using optimization models from a mathematical point of view [17]. In 2009, Bagherbandi et al. obtained the optimum configuration and observation weights in geodetic networks considering the postulated reliability and precision requirements [18]. In 2017, Horemuž and Jansson investigated the optimum horizontal location of the total station in a simple network in which a total station measures several control points in the process of being established [19]. In the following year, Alizadeh et al. investigated the optimal vertical placement of the total station with the variances criterion that related to the sum of total station coordinates and orientations [20]. The studies mentioned above mainly focused on the terrestrial geodetic network, while only a few scholars have studied the optimization of the seafloor PNT network. In 2020, Yang et al. gave the basic principles of the seafloor PNT network design. They believed that three sonar beacons can determine the 3D position of the underwater vehicle, so a equilateral triangle should be the optimal configuration of the unit network in the underwater PNT network. If the redundancy is taken into consideration, then at least four necessary reference points and one redundant reference point constitute the simplest network type, that is, an equilateral pentagon should be the optimal configuration of the unit network in the underwater PNT network [21]. However, the conclusion comes from experience and has not been proven by strict theory, and at present there is no strict optimization theory and design principle of the seafloor PNT network.

With the rapid development of GNSS positioning technology, the impact of a terrestrial geodetic networks’ configuration on positioning accuracy has become less prominent [22]. Moreover, because the electromagnetic wave is easily attenuated in water while the acoustic wave can travel great distances, the underwater positioning mode mainly depends on the acoustic wave [23]. Therefore, it is not applicable to directly transplant the optimization principle of the terrestrial geodetic network into the seafloor PNT network. Given the present situation, this paper studies the optimal design of the seafloor PNT network on the basis of Grafarend’s four orders for designing the geodesic network. The research tasks are as follows. Firstly, the positioning method of seabed control points in the underwater PNT network will be studied, based on which there are two main questions that need to be explored. One is the best configuration of the unit network, that is, how long should the side length of the unit network be and how many observation sessions should be carried out during the observation period under different positioning precisions. The other is how to ensure the positioning accuracy of seafloor control points during the expansion of the underwater PNT network.

The remainder of this paper is organized as follows. In Section 2, the measurement method and adjustment model of the seafloor PNT network are presented. In Section 3, the best configuration of the unit network in the seafloor PNT network is introduced. In Section 4, a chain coordinate transfer strategy is proposed to determine the distribution of the known control points in the expansion of the seafloor PNT network. Section 5 presents an experiment to prove the applicability of the proposed method in Section 2 and shows the layout schemes of seafloor PNT networks with centimeter, decimeter, and meter level precision on the basis of Section 3 and Section 4. Finally, conclusions are drawn in Section 6.

2. The Combined Measurement Method for the Seafloor PNT Network

It is necessary to first determine the coordinates of all control points in the seafloor PNT network according to the principles of underwater PNT. At present, the most commonly used method to determine the coordinates of underwater control points is circle-sailing positioning. The principle of circle-sailing positioning is as follows. The survey ship follows a circular trajectory, centering on the control point. During the sailing, the coordinates of the onboard transducer are provided by GNSS in real time, and the distance between the transducer and the acoustic beacon installed on the seafloor control point is measured by sound waves. Finally, the position of the control point is determined by the range intersection model constructed by combining the coordinates of the onboard transducer and the acoustic ranging results [24]. This method is highly accurate but time-consuming and laborious. Self-calibrating beacons simplify the task by reducing the surveying task to only one or a few beacons, with the others determining their own positions relative to the beacons that have been carried out by circle-sailing [25]. Due to the limited range capability of acoustic beacons, usually about 10 km [26], the relative positioning between beacons is only carried out within the unit network in the seafloor PNT network.

Therefore, the combined measurement method is proposed. As shown in Figure 2, there are two parts in the combined measurement method. The first part is the absolute positioning of several certain control points by circle-sailing positioning. The second part is the relative positioning among all control points, that is, the relative coordinates among seafloor control points are determined by constructing the free network adjustment model using the mutual acoustic ranging results provided by the acoustic beacons installed on them. The two parts can be carried out simultaneously, and the coordinates of all control points can be obtained by combining the results of the two.

2.1. The Adjustment Model of the Absolute Positioning

The model of the absolute positioning in a unit network is shown in Equation (1).

(1)$S_{tP}^i=\sqrt{{\left(X_t^i-X_P\right)}^2+{\left(Y_t^i-Y_P\right)}^2+{\left(Z_t^i-Z_P\right)}^2}+{\epsilon }^i$

where $S_{tP}^i$ is the geometric distance between the onboard transducer and the acoustic beacons installed on the seafloor control point P at the ith observation time, which is determined by acoustic ray tracing based on the one-way travel time between the two; $X_t=\left(X_t^i,Y_t^i,Z_t^i\right)$ and ${\mathit{X}}_{\mathit{P}}=\left(X_P,Y_P,Z_P\right)$ are the coordinates of the transducer at the ith observation time and the control point P, respectively; ɛⁱ is the observation error.

Linearize Equation (1) to obtain the error equation:

(2)$v_{tP}^i=\frac{X_t^i-X_P^0}{S_{tP}^{i0}}dx+\frac{Y_t^i-Y_P^0}{S_{tP}^{i0}}dy+\frac{Z_t^i-Z_P^0}{S_{tP}^{i0}}dz-\left(S_{tP}^{i0}-S_{tP}^i\right)$

where $S_{tP}^{i0}=\sqrt{{\left(X_t^i-X_P^0\right)}^2+{\left(Y_t^i-Y_P^0\right)}^2+{\left(Z_t^i-Z_P^0\right)}^2}$ and ${\mathit{X}}_{\mathit{P}}^{\mathit{0}}=\left(X_P^0,Y_P^0,Z_P^0\right)$ is the prior coordinates of the control point P. A series of observation equations with the form of Equation (2) can be represented by the matrix notation, as shown in Equation (3).

(3)$\underset{n\times 1}{V_a}=\underset{n\times 3}{B_a}\underset{3\times 1}{x_a}-\underset{n\times 1}{L_a}$

where $\left[V_a{}_{,i}\right]=v_{tP}^i$, $\left[B_{a,}{}_i\right]=\left[\begin{array}{ccc}\frac{X_t^i-X_P^0}{S_{tP}^{i0}}& \frac{Y_t^i-Y_P^0}{S_{tP}^{i0}}& \frac{Z_t^i-Z_P^0}{S_{tP}^{i0}}\end{array}\right]$, $x_a={\left[\begin{array}{ccc}dx& dy& dz\end{array}\right]}^T$, and $\left[L_{a,}{}_i\right]=S_{tP}^{i0}-S_{tP}^i$; the symbol “[A_i]” represents the ith row of the matirx A, the superscript “T” represents the transpose of the matrix, and the subscript “a” represents the absolute positioning.

On the basis of the least-squares criterion, the coordinate corrections ${\widehat{x}}_a$ and the final coordinates X_P of the control point can be acquired through Equation (4).

(4)${\widehat{x}}_a={\left(B_a{}^T{W}_a{B}_a\right)}^{-1}{B}_a{}^T{W}_a{L}_a,X_p={\mathit{X}}_{\mathit{P}}^{\mathit{0}}+{\widehat{x}}_a$

where W_a is the weight matrix of all observations; generally it is a diagonal matrix, and $W_{a,}{}_{ii}=1/\left|L_a{}_{,i}\right|$.

The precision of the absolute positioning is given in Equation (5).

(5)$D_{\widehat{x}\widehat{x}}={\sigma }_0^2{Q}_{\widehat{x}\widehat{x}},{\sigma }_0^2=\frac{V_a{}^T{W}_a{V}_a}{n-3},Q_{\widehat{x}\widehat{x}}={\left(B_a{}^T{W}_a{B}_a\right)}^{-1}$

where $D_{\widehat{x}\widehat{x}}$ is the variance matrix of the control point coordinates, ${\widehat{\sigma }}_0$ is the standard deviation, and $Q_{\widehat{x}\widehat{x}}$ is the covariance matrix of the control point coordinates.

When the absolute positioning is completed, the coordinates of the control point, which has been carried out by the circle-sailing positioning, can be known. Therefore, in this paper, the control point is called “the known point”.

2.2. The 2D Adjustment Model under Depth Constraint in the Relative Positioning

Generally, there is a pressure sensor installed on the seafloor control point that can provide the high-precision depth from the sea surface to the control point. The elevation of the control point can then be obtained by subtracting the depth from the elevation of the sea surface. In this case, the 3D adjustment model in the relative positioning of a unit network can be transformed into a 2D adjustment model. As shown in Figure 3, the geometric distance, S_pq, between control points p and q can be obtained by in-layer constant gradient acoustic ray tracing using the acoustic observed value, then the horizontal distance, D_pq, between the two control points can be calculated based on Equation (6) by combining S_pq with the elevation difference, ΔZ_pq [27].

(6)$D_{pq}=\sqrt{S_{pq}^2-\Delta Z_{pq}^2},\Delta Z_{pq}=\left|Z_p-Z_q\right|$

If there are m control points in the unit network, then the number of measured horizontal distance observations, n, is equal to m(m − 1)/2. Supposing that the horizontal coordinates of two control points, p and q, are $\left(X_P,Y_P\right)$ and $\left(X_q,Y_q\right)$, respectively, then the horizontal distance observation, D_pq, between p and q can be expressed by Equation (7).

(7)$D_{pq}=\sqrt{{\left(X_p-X_q\right)}^2+{\left(Y_p-Y_q\right)}^2}+{\epsilon }_{pq}$

where ɛ_pq is the observation error. Equation (7) is linearized to obtain the error equation:

(8)$v_{pq}=-\frac{X_p-X_q}{D_{pq}^0}d{\widehat{x}}_p-\frac{Y_p-Y_q}{D_{pq}^0}d{\widehat{y}}_p+\frac{X_p-X_q}{D_{pq}^0}d{\widehat{x}}_q+\frac{Y_p-Y_q}{D_{pq}^0}d{\widehat{y}}_q-\left(D_{pq}^0-D_{pq}\right)$

where $D_{pq}^0=\sqrt{{\left(X_p^0-X_q^0\right)}^2+{\left(Y_p^0-Y_q^0\right)}^2}$ and ${\mathit{X}}_{\mathit{p}}^0=\left(X_p^0,Y_q^0\right)$; ${\mathit{X}}_{\mathit{q}}^0=\left(X_q^0,Y_q^0\right)$ are the prior coordinates of control points p and q, respectively.

The error equations of n horizontal distance observations can be written into matrix form, as shown in Equation (9).

(9)$\underset{n\times 1}{V_r}=\underset{n\times 2m}{B_r}\underset{2m\times 1}{x_r}-\underset{n\times 1}{L_r}$

where the subscript “r” represents the relative positioning. V_r, B_r, x_r, and L_r are given in Equation (10).

(10) $\begin{array}{c}\left[{V_{r,}}_i\right]=v_{pq},\left[{L_{r,}}_i\right]=D_{pq}^0-D_{pq},\underset{2m\times 1}{{\widehat{x}}_r}={\left[\begin{array}{ccccccc}d{\widehat{x}}_1& d{\widehat{y}}_1& d{\widehat{x}}_2& d{\widehat{y}}_2& \cdots & d{\widehat{x}}_m& d{\widehat{y}}_m\end{array}\right]}^T\hfill \\ \left[{B_{r,}}_i\right]=\left[\begin{array}{ccccccccccccc}0& \cdots & 0& -\frac{X_p-X_q}{D_{pq}^0}& -\frac{Y_p-Y_q}{D_{pq}^0}& 0& \cdots & 0& \frac{X_p-X_q}{D_{pq}^0}& \frac{Y_p-Y_q}{D_{pq}^0}& 0& \cdots & 0\end{array}\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\stackrel{\uparrow }{{\displaystyle \mathit{column}}{\displaystyle \left(2p-1\right)}} \stackrel{\uparrow }{{\displaystyle \mathit{column}}{\displaystyle \left(2p\right)}}\hspace{1em}\stackrel{\uparrow }{{\displaystyle \mathit{column}}{\displaystyle \left(2q-1\right)}} \stackrel{\uparrow }{{\displaystyle \mathit{column}}{\displaystyle \left(2q\right)}}\end{array}$

As shown in Figure 4, only one known control point and one known azimuth are provided, and the question is, can the 2D unit network be fixed in a geodetic coordinate system? Because there are no known conditions during the relative positioning, design matrix B is rank-deficient and the number of rank deficiency, d, is equal to 3.

The rank defect network adjustment under the barycentric datum is adopted to solve Equation (10), and the additional datum constraint condition is given in Equation (11).

(11)$\underset{d\times 2m}{S^T}\underset{2m\times 2m}{W_x}\underset{2m\times 1}{{\widehat{x}}_r}=0$

where S is the base condition matrix and meets the following conditions: $rank(S)=d$, $rank\left[\begin{array}{l}{B}_r\\ S^T\end{array}\right]=2m$, $B_{r}S=0$; W_x is the base weight matrix and is an identity matrix. The specific expression of S^T is shown in Equation (12).

(12)$\underset{d\times 2m}{S^T}=\left[\begin{array}{ccccccc}1& 0& 1& 0& \cdots & 1& 0\\ 0& 1& 0& 1& \cdots & 0& 1\\ -Y_1^0& X_1^0& -Y_2^0& X_2^0& \cdots & -Y_m^0& X_m^0\end{array}\right]$

Combining Equations (9) and (11), the coordinate corrections, ${\widehat{x}}_r$, and the final coordinates, X, of the control points can be obtained by Equation (13) under the adjustment criterion of $V_r{}^T{W}_r{V}_r=\min$.

(13)${\widehat{x}}_r={\left(B_r{}^T{W}_r{B}_r+W_{x}S{S}^T{W}_x\right)}^{-1}{B}_r{}^T{W}_r{L}_r,X=X_0+{\widehat{x}}_r$

The precision of the relative positioning is given in Equation (14).

(14)$\begin{array}{l}{D}_{\widehat{x}\widehat{x}}={\widehat{\sigma }}_0^2{Q}_{\widehat{x}\widehat{x}},{\widehat{\sigma }}_0^2=\frac{V_r{}^T{W}_r{V}_r}{n-\left(2m-d\right)}\\ Q_{\widehat{x}\widehat{x}}=Q’B_r{}^T{W}_r{B}_{r}Q’,Q’={\left(B_r{}^T{W}_r{B}_r+W_{x}S{S}^T{W}_x\right)}^{-1}\\ GDO{P}_i=\sqrt{Q_{\widehat{x}\widehat{x}}\left(2i-1,2i-1\right)+Q_{\widehat{x}\widehat{x}}\left(2i,2i\right)}\end{array}$

where $D_{\widehat{x}\widehat{x}}$ is the variance matrix of the control point coordinates, ${\widehat{\sigma }}_0$ is the standard deviation, $Q_{\widehat{x}\widehat{x}}$ is the covariance matrix of the control point coordinates, and GDOP is the geometric dilution of precision of the control point.

The relative coordinates between the known point and the other control points can be acquired after the model of the relative positioning is solved, and then the coordinates of all control points can be calculated by combining the coordinates of the known point and the relative coordinates.

3. The Best Configuration of the Unit Network in the Seafloor PNT Network

In the seafloor PNT network, the longer the unit network side length, the lower the control point density and the lower the layout cost. In the relative positioning, the more acoustic ranging observation sessions there are, the higher the ranging accuracy. In order to balance the layout cost and observation efficiency, it is necessary to study the relationship between the side length of the unit network and the number of observation sessions.

It can be proved that the coverage area of the regular polygon is the largest in the same circumcircle. Therefore, in order to achieve the maximum coverage under the same number of control points, the unit network of the seafloor PNT network is generally arranged as a regular polygon. As shown in Figure 5, in order to achieve seamless full coverage of the water area, the unit network is designed as an equilateral triangle (GDOP = 0.808 when the side length of the unit network is 1 km) or square (GDOP = 0.802) in most cases. Additionally, when the navigation and positioning accuracy requirements are undemanding, the unit network may be designed as a regular hexagon (GDOP = 0.837).

As mentioned in Section 2, the ranging capability of the acoustic beacons is limited. If the upper limit, R, of the range of the acoustic beacon is known, the maximum side length, L_max, of the equilateral triangle and the square unit network are equal to R and $\sqrt{2}R/2$, respectively. In the following, the functional relationship between the side length of the unit network and the navigation and positioning accuracy of the control points will be derived.

As shown in Figure 6, there are two known control point in an equilateral triangle unit network, and the distance, S₀, and the azimuth, α₀, between points A and B can be acquired based on the geometric relations. The coordinates (X_P, Y_P) of the unknown control point P can then be obtained by Equation (15).

(15)$X_P=X_0+S\cos \alpha ,Y_P=Y_0+S\sin \alpha$

where S is the distance between point P and A; α is the azimuth between the two points and it meet the condition in Equation (16).

(16)$\alpha ={\alpha }_0-\left({180}^{\circ }-\angle B-\angle P\right)$

Based on the law of cosines, Equation (17) can be given in the triangle ΔABP.

(17)$S^2=S_0^2+{S^{\prime }}^2-2S_0{S}^{\prime }\cos B$

where S′ is the distance between point P and B. S and S′ can be measured by acoustic ranging, and their ranging accuracy are known. By differentiating both sides of Equation (17), we can acquire:

(18)$SdS=\left(S_0-S^{\prime }\cos B\right)d{S}_0+\left(S^{\prime }-S_0\cos B\right)d{S}^{\prime }+\left(S_0{S}^{\prime }\sin BdB\right)/\rho$

where ρ is the conversion parameter between angle and radian and ρ = 180°/π. Because S₀ is known, which means dS₀ is equal to 0, Equation (18) can be simplified as Equation (19).

(19)$dB=\frac{S\rho dS-\left(S^{\prime }-S_0\cos B\right)\rho d{S}^{\prime }}{S_0{S}^{\prime }\sin B}$

By combining Equations (16) and (19), the following equation can be obtained:

(20)$\left[\begin{array}{c}dS\\ d\alpha \end{array}\right]=\left[\begin{array}{cc}1& 0\\ x_1& x_2\end{array}\right]\left[\begin{array}{c}dS\\ dS’\end{array}\right]$

where $x_1=\frac{S\rho }{S_0{S}^{\prime }\sin B}, x_2=\frac{\left(S’-S_0\cos B\right)\rho }{S_0{S}^{\prime }\sin B}$.

In the equilateral triangle unit network, the ranging accuracy of S is equal to that of S′. Supposing the ranging accuracy of S and S′ is σ_S, then the covariance matrix between S and α can be depicted as Equation (21) by the variance–covariance propagation law.

(21)$\left[\begin{array}{cc}{\sigma }_S^2& {\sigma }_{S\alpha }\\ {\sigma }_{\alpha S}& {\sigma }_{\alpha }^2\end{array}\right]=\left[\begin{array}{cc}{\sigma }_S^2& x_1{\sigma }_S^2\\ x_1{\sigma }_S^2& \left(x_1^2+x_2^2\right){\sigma }_S^2\end{array}\right]$

Equation (15) can be differentiated and the result written into matrix form, as shown in Equation (22).

(22)$\left[\begin{array}{c}dX\\ dY\end{array}\right]=\left[\begin{array}{cc}\cos \alpha & -S\sin \alpha /\rho \\ \sin \alpha & S\cos \alpha /\rho \end{array}\right]\left[\begin{array}{c}dS\\ d\alpha \end{array}\right]$

Then, according to the variance–covariance propagation law, Equation (23) can be obtained.

(23)$\left[\begin{array}{cc}{\sigma }_X^2& {\sigma }_{XY}\\ {\sigma }_{YX}& {\sigma }_Y^2\end{array}\right]=\left[\begin{array}{cc}\cos \alpha & -S\sin \alpha /\rho \\ \sin \alpha & S\cos \alpha /\rho \end{array}\right]\left[\begin{array}{cc}{\sigma }_S^2& {\sigma }_{S\alpha }\\ {\sigma }_{\alpha S}& {\sigma }_{\alpha }^2\end{array}\right]\left[\begin{array}{cc}\cos \alpha & \sin \alpha \\ -S\sin \alpha /\rho & S\cos \alpha /\rho \end{array}\right]$

Based on Equation (23), the coordinate accuracy, σ_P, of point P is given by Equation (24).

(24)${\sigma }_P^2={\sigma }_X^2+{\sigma }_Y^2={\sigma }_S^2\left[1+S^2\left(x_1^2+x_2^2\right)/{\rho }^2\right]$

By substituting x₁ and x₂ into Equation (24) and considering the geometric conditions S = S′ = S₀ and ∠B = ∠P = 60° in the equilateral triangle unit network, Equation (25) can be obtained.

(25)${\sigma }_S^2={\sigma }_P^2/2.667$

According to the same principles, the relation between ranging accuracy, σ_S, and the coordinate accuracy, σ_P, of the unknown point in the square unit network can also be derived, as shown in Equation (26).

(26)${\sigma }_S^2={\sigma }_P^2/11$

Because the distance, S, between two control points is determined by the sound velocity, c, and the one-way propagation time, t, of the acoustic wave, we can acquire Equation (27).

(27)$dS=tdc+cdt$

Suppose that the measurement accuracy of sound velocity and time are σ_c and σ_t, respectively, then, based on the law of error propagation, it can be known that the ranging accuracy, σ_S, meets the condition in Equation (28).

(28)${\sigma }_S^2=t^2{\sigma }_c^2+c^2{\sigma }_t^2={\left(S/c\right)}^2{\sigma }_c^2+c^2{\sigma }_t^2$

Substituting Equation (28) into Equations (25) and (26), the function of the side length, S, and the coordinate accuracy of control points, σ_P, in one observation period can be obtained, as shown in Equation (29), where the subscripts “3” and “4” represent the number of edges in a unit network.

(29)$S_3=\sqrt{\frac{c^2{\sigma }_P^2-2.667c^4{\sigma }_t^2}{2.667{\sigma }_c^2}},S_4=\sqrt{\frac{c^2{\sigma }_P^2-11c^4{\sigma }_t^2}{11{\sigma }_c^2}}$

The limiting conditions for the above formula are ${\sigma }_P>\sqrt{2.667}c{\sigma }_t$ and ${\sigma }_P>\sqrt{11}c{\sigma }_t$. If there are N observation sessions, that is, each side in the unit network has been observed for N times, the ranging accuracy, σ_S, satisfies Equation (30).

(30)${\sigma }_S^2=\left[{\left(S/c\right)}^2{\sigma }_c^2+c^2{\sigma }_t^2\right]/N$

The relation between the side length and the coordinate accuracy of control points in N observation sessions can then be given by Equation (31).

(31)$S_3=\sqrt{\frac{N{c}^2{\sigma }_P^2-2.667c^4{\sigma }_t^2}{2.667{\sigma }_c^2}},S_4=\sqrt{\frac{N{c}^2{\sigma }_P^2-11c^4{\sigma }_t^2}{11{\sigma }_c^2}}$

To sum up, in the case that the positioning accuracy of the control points in the seafloor PNT network and the number of observation sessions are given, the side length of the unit network can be determined based on Equation (31). Furthermore, if the positioning precision and the side length of the unit network have been given, the optimal number of observation sessions can be determined.

4. The Distribution of Known Control Points in the Seafloor PNT Network

Due to the large area of the seafloor PNT network, it is time-consuming and high cost to carry out circle-sailing positioning in each unit network in the absolute positioning. To solve this problem, the chain coordinate transfer strategy is proposed, whose principle is as follows.

As shown in Figure 7: Firstly, carry out the absolute positioning and the relative positioning in the initial unit network A to determine the coordinates of all control points in it. Then, utilize the results of acoustic ranging in the unit network B adjacent to A to obtain the angle (Angles 1, 2, 3, and 4) between each edge in B. Finally, determine the coordinates of unknown points in B using the coordinates of common points between A and B, the acoustic ranging results, and the computed angles in B. Once the coordinates of the control points in B are determined, it can be regarded as “the initial unit network” relative to the unit network C. The coordinates of the control points in other unit networks (C, D, E, etc.) can then be determined by repeating the process.

During the coordinate transfer in the seafloor PNT network, the positioning error will increase with the increase of the number of unit networks. Circle-sailing positioning should be carried out at the control point whose positioning accuracy is less than the lower limit of accuracy (called “the weakest point” for short) to calibrate its coordinates. Therefore, determining the position of the weakest point is inevitable.

As shown in Figure 8, suppose that the diagonal length between the starting known point P₁ and the weakest point P_n is L. Due to the error of acoustic ranging, the computed position of P_n is offset to P_n″, the offset in distance is f (P_n P_n″), and the offset in azimuth is dA. The offset f can be decomposed into the longitudinal offset t and the transverse offset u, which are equal to P_nP_n′ and P_n′P_n″, respectively. The offsets f, t, and u meet the condition in Equation (32).

(32)$f^2=t^2+u^2=t^2+d{A}^2{L}^2/{\rho }^2$

Suppose that the positioning accuracy of P_n is σ_P and the mean square error of t, u, and dA are σ_t, σ_u, and σ_A(L), respectively, the relation among them is given by Equation (33).

(33)${\sigma }_P{}^2={\sigma }_t^2+{\sigma }_u^2={\sigma }_t^2+{\sigma }_A{}^2\left(L\right)L^2/{\rho }^2$

If the number of unit networks shown in Figure 8 is odd, namely there are (2n − 1) unit networks with side length S in the seafloor PNT network, and the observed distance between the control points located on L in the ith unit network is S_i, then Equation (34) can be obtained.

(34)$L=S_1+S_2+\cdots +S_n$

Suppose that the length of the starting edge is S₀, then, according to the geometric relations in the equilateral triangle unit network, we can obtain Equation (35).

(35)$L=S_0\frac{\sin \left({\theta }_3\right)}{\sin \left({\theta }_2\right)}+S_0\frac{\sin \left({\theta }_1\right)\sin \left({\theta }_5\right)\sin \left({\theta }_9\right)}{\sin \left({\theta }_2\right)\sin \left({\theta }_4\right)\sin \left({\theta }_8\right)}+\cdots +S_0\frac{\sin \left({\theta }_1\right)\sin \left({\theta }_5\right)\sin \left({\theta }_7\right)\cdots \sin \left({\theta }_{6n-3}\right)}{\sin \left({\theta }_2\right)\sin \left({\theta }_4\right)\sin \left({\theta }_8\right)\cdots \sin \left({\theta }_{6n-4}\right)}$

where θ₁, θ₂, …, θ_6n−4 are the computed angles in the unit networks, which can be calculated based on the law of cosines, and they satisfy the following conditions:

(36)${\tilde{\theta }}_{6n-5}+{\tilde{\theta }}_{6n-4}+{\tilde{\theta }}_{6n-3}-{180}^{\circ }=0$

where ${\tilde{\theta }}_i$ is the true value of θ_i. The equations corresponding to all unit networks in the form of Equation (36) can be written into the function model of conditional adjustment, as shown in Equation (37).

(37)$\underset{n\times \left(6n-3\right)}{\left[\begin{array}{cccccccccc}1& 1& 1& 0& 0& 0& \cdots & 0& 0& 0\\ 0& 0& 0& 1& 1& 1& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & 1& 1& 1\end{array}\right]}\underset{\left(6n-3\right)\times 1}{\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_{6n-3}\end{array}\right]}+\underset{n\times 1}{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_n\end{array}\right]}=0$

where v_i is the deviation between ${\tilde{\theta }}_i$ and its corresponding calculated value θ_i; $w_n=\left({\theta }_{6n-5}+{\theta }_{6n-4}+{\theta }_{6n-3}-{180}^{\circ }\right)$. The matrix form of Equation (37) is given by Equation (38).

(38)$Av+W=0$

Equation (35) can be expressed as a function between L and all computed angles, as shown in Equation (39).

(39)$L=\phi \left[{\tilde{\theta }}_1,{\tilde{\theta }}_2,{\tilde{\theta }}_3,\cdots {\tilde{\theta }}_{6n-3}\right]+\frac{\partial L}{\partial {\theta }_1}{v}_1+\frac{\partial L}{\partial {\theta }_2}{v}_1+\cdots +\frac{\partial L}{\partial {\theta }_{6n-3}}{v}_{6n-3}$

By differentiating Equation (39), the following equation can be gained:

(40)$d\widehat{L}=\frac{\partial L}{\partial {\theta }_1}d{v}_1+\frac{\partial L}{\partial {\theta }_2}d{v}_1+\cdots +\frac{\partial L}{\partial {\theta }_{6n-3}}d{v}_{6n-3}$

where $\frac{\partial L}{\partial {\theta }_i}$ meets the rule in Equation (41).

(41)$\begin{array}{l}\frac{\partial L}{\partial {\theta }_1}=\left(n-1\right)S_{0}c,\frac{\partial L}{\partial {\theta }_2}=-n{S}_{0}c,\frac{\partial L}{\partial {\theta }_3}=\frac{\partial L}{\partial {\theta }_9}=\frac{\partial L}{\partial {\theta }_{15}}=\cdots =\frac{\partial L}{\partial {\theta }_{6n-3}}=S_{0}c,\\ \frac{\partial L}{\partial {\theta }_6}=\frac{\partial L}{\partial {\theta }_{12}}=\frac{\partial L}{\partial {\theta }_{18}}=\cdots =\frac{\partial L}{\partial {\theta }_{6n-6}}=0,\\ \frac{\partial L}{\partial {\theta }_7}=-\frac{\partial L}{\partial {\theta }_{10}}=\frac{\partial L}{\partial {\theta }_{11}}=-\frac{\partial L}{\partial {\theta }_{14}}=\cdots =\left(n-2\right)S_{0}c,c=ctg\left({60}^{\circ }\right)=\sqrt{3}/3\end{array}$

Let $\frac{\partial L}{\partial \theta }={\left[\begin{array}{ccccccc}\frac{\partial L}{\partial {\theta }_1}& \frac{\partial L}{\partial {\theta }_2}& \frac{\partial L}{\partial {\theta }_3}& \cdots & \frac{\partial L}{\partial {\theta }_{6n-5}}& \frac{\partial L}{\partial {\theta }_{6n-4}}& \frac{\partial L}{\partial {\theta }_{6n-3}}\end{array}\right]}^T$, then the reciprocal of the weight of the adjustment function of Equation (39) can be given by Equation (42) according to the variance–covariance propagation law and the principle of condition adjustment.

(42)$\frac{1}{W_{\widehat{L}}}=Q_{\widehat{L}\widehat{L}}={\left(\frac{\partial L}{\partial \theta }\right)}^T{W}^{-1}\frac{\partial L}{\partial \theta }-{\left(A{W}^{-1}\frac{\partial L}{\partial \theta }\right)}^T{\left(A^{T}WA\right)}^{-1}\left(A{W}^{-1}\frac{\partial L}{\partial \theta }\right)$

where W is the weight matrix of computed angles and is an identity matrix; A is the design matrix of the condition adjustment function model shown in Equation (38).

Based on Equation (41), the components in Equation (42) can be calculated.

(43)$\begin{array}{l}{\left(\frac{\partial L}{\partial \theta }\right)}^T\frac{\partial L}{\partial \theta }=S^2{c}^2\left\{n^2+n+4\left[{\left(n-1\right)}^2+{\left(n-2\right)}^2+{\left(n-3\right)}^2+\cdots 2^2+1^2\right]\right\}\\ {\left(A\frac{\partial L}{\partial \theta }\right)}^T{\left(A^{T}A\right)}^{-1}\left(A\frac{\partial L}{\partial \theta }\right)=0\end{array}$

If Equation (43) is substituted into Equation (42), then we can obtain Equation (44).

(44)$\frac{1}{W_{\widehat{L}}}=S^2{c}^2\left(\frac{4n^3-3n^2+5n}{3}\right)$

Substitute c² = 1/3 into the above formula and take into consideration the condition that the truth value of L is equal to n times S₀, that is, $\tilde{L}=n{S}_0$, then the reciprocal of the weight of the diagonal length L between the starting point P₁ and the weakest point P_n after condition adjustment is:

(45)${\left(\frac{1}{W_{\widehat{L}}}\right)}_{odd}={\tilde{L}}^2\left(\frac{4n^2-3n+5}{9n}\right)$

If the number of unit networks shown in Figure 8 is even, namely there are (2n) unit networks in the seafloor PNT network, based on the same derivation process we can acquire:

(46)${\left(\frac{1}{P_{\widehat{L}}}\right)}_{even}={\tilde{L}}^2\left(\frac{4n^2+3n+5}{9n}\right)$

Considering the error of the starting edge, the mean square error of the longitudinal offset ${\sigma }_{t_3}$ in the seafloor PNT network composed of the equilateral triangle unit network can be given by Equation (47).

(47)${\sigma }_{t_3}={\sigma }_{L_3}=\pm \tilde{L}\sqrt{{\left(\frac{{\sigma }_{S_0}}{\tilde{L}}\right)}^2+{\left(\frac{{\sigma }_{\alpha }}{\rho }\right)}^2\left(\frac{4n^2\mp 3n+5}{9n}\right)}$

where ${\sigma }_{S_0}$ is the error of the starting edge and σ_α is the accuracy of the computed angles, which can be calculated based on Equation (21). Based on the same derivation process, the mean square error of the transverse offset ${\sigma }_{u_3}$ can be given by Equation (48).

(48)${\sigma }_{u_3}=\tilde{L}\sqrt{{\left(\frac{{\sigma }_{{\alpha }_0}}{\rho }\right)}^2+{\left(\frac{{\sigma }_{\alpha }}{\rho }\right)}^2\left(\frac{4n^2-3n+5}{9n}\right)}$

where ${\sigma }_{{\alpha }_0}$ is the error of the starting azimuth. Similarly, the mean square error of the longitudinal and transverse offset in the seafloor geodetic PNT composed of the square unit network can also be acquired, as shown in Equation (49).

(49)${\sigma }_{t_4}=\tilde{L}\sqrt{{\left(\frac{{\sigma }_{S_0}}{\tilde{L}}\right)}^2+{\left(\frac{{\sigma }_{\alpha }}{\rho }\right)}^2\frac{n}{3}},{\sigma }_{u_4}=\tilde{L}\sqrt{{\left(\frac{{\sigma }_{{\alpha }_0}}{\rho }\right)}^2+{\left(\frac{{\sigma }_{\alpha }}{\rho }\right)}^2\frac{n}{3}}$

Finally, it can be concluded that the positioning accuracy of the weakest point is ${\sigma }_{P_n}=\sqrt{{\sigma }_t{}^2+{\sigma }_u{}^2}$. If ${\sigma }_{P_n}$ is greater than the lower limit of the positioning precision, σ_max, the equation ${\sigma }_{\max }=\sqrt{{\sigma }_t{}^2+{\sigma }_u{}^2}$ can be solved to obtain n, then the position of the weakest point where the absolute positioning should be carried out can be gained.

5. Experiments and Analysis

5.1. The Applicability of the Combined Measurement Method in the Seafloor PNT Network

To test the applicability of the combined measurement method in Section 2, it is utilized to solve the coordinates of seafloor control points in a unit network laid in Songhua Lake, Jilin Province, China in October 2012. The seabed for laying the control points is flat, with a depth of 63~64 m, 320 m from north to south, and 400 m from east to west. The topographic map of the experimental water area is shown in Figure 9a.

The unit network is composed of four control points, each of which was carried out by circle-sailing positioning and conducting acoustic ranging of 50 epochs between each other. The structure of the unit network is shown in Figure 9b. Using the data collected in the circle-sailing positioning, the coordinates of the four seafloor control points can be calculated according to the adjustment model of the absolute positioning, as shown in Table 1.

The acoustic ranging results of different combinations of control points (Figure 10) were used to calculate the coordinates of control points according to the relative positioning model. The mean and standard deviation of mutual ranging results between control points are given in Table 2. The difference between the relative positioning results and the absolute positioning results is used as an index to evaluate the accuracy of the relative positioning. The results are shown in Table 3.

The following conclusions can be drawn from Table 3.

• In the Combination C₁C₂C₃C₄, the maximum and minimum of the coordinate difference between the relative positioning and the absolute positioning in the X direction are 0.115 m and 0.009 m, respectively; the maximum and minimum of the coordinate difference in the Y direction are 0.130 m and 0.026 m; the positioning difference is approximately 10 cm.

• In the Combinations C₂C₃C₄ and C₁C₂C₄, the positioning differences between the relative positioning and the absolute positioning are very slight, and less than 10 cm in most cases.

• In the Combinations C₁C₂C₃ and C₁C₃C₄, the positioning differences between the relative positioning and the absolute positioning are quite large relative to other combinations, and they are up to the decimeter level. The reason for this is that the standard deviation of the acoustic ranging results between point C1 and point C3 is very large, which means the data quality is poor.

It can be concluded that if the acoustic ranging results are high-quality, the positioning difference between the combined measurement method and circle-sailing is at the centimeter level, which can satisfy the requirement of underwater PNT for the most part. Therefore, it is feasible to use the combined measurement method in the seafloor PNT network.

5.2. The Design of the Seafloor PNT Network under Different Accuracy Levels

In general, the seafloor PNT network is divided into the centimeter-level network, decimeter-level network, and meter-level network, according to different uses. In this section, the sound velocity, c, in Equation (28) is set to 1500 m/s, the time measuring accuracy of acoustic beacons is set to 10⁻⁵ s, and the sound velocity measuring accuracy is set to 0.02 m/s, referring to the Compatt 6+ transponder of Sonardyne. In this case, the highest positioning precision of control points in the seafloor PNT network composed of the equilateral triangle unit network can reach 0.024 m in one observation period, and that of the seafloor PNT network composed of the square unit network can reach 0.05 m. Therefore, two centimeter-level networks (σ_P = 0.025 m and 0.05 m), three decimeter-level networks (σ_P = 0.1 m, 0.2 m, and 0.5 m) and one meter-level network (σ_P = 1 m) are set. Based on Equation (31), the relations between the side length of the unit network (unit: kilometer) and the number of observation sessions at all precision levels, and the relations between the total observation time (unit: minute) in the unit network and observation sessions can be calculated, as shown in Figure 11 and Figure 12.

From Figure 11 and Figure 12, the following conclusions can be drawn.

• In the same observation sessions, the side length of the unit network decreases with the increase of the positioning precision of the control points.

• The total observation time increases exponentially with the increase of observation sessions and the side length of the unit network.

• When the positioning precision of control points in the seafloor PNT network is given, the more observation sessions, and the longer the side length of the unit network, the smaller the density of control points and the lower the layout cost.

• The side length will not increase infinitely with the increase of observation sessions, and it will gradually approach a certain value. Thereafter, increasing observation sessions has little effect on the side length, but the total observation time will still rise, resulting in low efficiency.

In order to reduce the density of control points so as to reduce the layout cost, the side length of the unit network should be as long as possible within the ranging capability. If the ranging capacity, R, of the transducer is known, then the maximum side length of the equilateral triangle unit network can be set to R and the maximum side length of the square unit network can be set to $\sqrt{2}R/2$. According to the results in Figure 11 and Figure 12, the number of observation sessions in the unit network to achieve different accuracy levels when R is set to 5 km, 10 km, 15 km, and 20 km is shown in Table 4, where “Δ” represents the equilateral triangle unit network, “*“ represents the square unit network, and “N(T)” represents the total observation time T of N observation sessions (unit: minutes).

In this paper, given observation efficiency, the total observation time in a unit network is limited to no more than 2 min. The observation sessions meeting this condition are marked in bold in Table 4. It can be concluded that if the total observation time is less than 2 min, then:

• When the side length is 5 km, the equilateral triangle and the square unit network can both meet the accuracy requirements of all levels.

• When the side length is 10 km, the equilateral triangle and the square unit network can meet the accuracy requirements of 0.1 m and below.

• When the side length is 15 km, the equilateral triangle unit network can meet the accuracy requirements of 0.1 m and below, while the square unit network can only meet the accuracy requirements of 0.2 m and below.

• When the side length is 20 km, the equilateral triangle and the square unit network can both meet the accuracy requirements of 0.2 m and below.

Based on Table 4, it can be seen that there are one, three, eight, sixteen, and twenty-two designing schemes for the seafloor PNT network when the positioning accuracy of control points is 0.05 m, 0.1 m, 0.2 m, 0.5 m, and 1 m under different unit network side lengths and different observation sessions. The positions of the weakest points in seafloor geodetic networks at all precision levels corresponding to different schemes are calculated according to Equations (47)–(49). The results are given in Table 5, Table 6 and Table 7, where “N(n)” represents the weakest points in the nth unit network; N is the original calculated value of n.

The TOPSIS method was utilized to evaluate each scheme by comprehensively considering the observation efficiency (the total observation time to the total control area), cost (the ratio of the number of control points to the total control area), and the number to carry out circle-sailing positioning (the number of known points). The comprehensive scores of each scheme are shown in Table 8.

The following conclusions can be drawn from the columns in the table above:

• The equilateral triangle unit network with the side length (S) of 5 km and observation sessions (T_n) of 20 can be selected when the positioning accuracy, σ_P, of the control points in the seafloor PNT network is 0.05 m.

• When σ_P = 0.1 m, the optimal scheme is the equilateral triangle unit network with S equal to 5 km and T_n equal to 5; the alternative schemes are the equilateral triangle unit network with S equal to 5 km and T_n equal to 20 or the square unit network with S equal to 5 km and T_n equal to 11.

• When σ_P = 0.2 m, the optimal scheme is the equilateral triangle unit network with S equal to 10 km and T_n equal to 5; the alternative scheme is the equilateral triangle unit network with S equal to 15 km and T_n equal to 11.

• When σ_P = 0.5 m, the optimal scheme is the equilateral triangle unit network with S equal to 15 km and T_n equal to 3; the alternative schemes are the equilateral triangle unit network with S equal to 20 km and T_n equal to 5, the equilateral triangle unit network with S equal to 10 km and T_n equal to 5 or 2, the square unit network with S equal to 20 km and T_n equal to 10, and the square unit network with S equal to 10 km and T_n equal to 11.

• When σ_P = 1 m, the optimal scheme is the equilateral triangle unit network with S equal to 15 km and T_n equal to 3 or 1; the alternative schemes are the equilateral triangle unit network with S equal to 20 km and T_n equal to 5 or 1, the equilateral triangle unit network with S equal to 10 km and T_n equal to 1 or 2, the square unit network with S equal to 15 km and T_n equal to 6, and the square unit network with S equal to 20 km and T_n equal to 10.

In this section, the design schemes of the seafloor PNT network are given when the velocity measuring accuracy, σ_c, is set as 0.02 m/s and the time measuring accuracy, σ_t, is set as 10⁻⁵ s. In practical scheme design, the corresponding condition can be substituted into the model given in this paper to obtain the specific layout scheme.

6. Conclusions

In this paper, the measurement method and the configuration and coordinate transfer strategy are studied from the conception of the construction of the seafloor PNT network. The main contributions are as follows.

6.1. A Combined Measurement Method That Is Less Time Consuming and Less Costly

In order to reduce the influence of the attitude on the positioning results, the ship speed is slow in the circle-sailing positioning, which leads to a long observation time. Moreover, it takes time, and the cost is very high to carry out circle-sailing positioning for each control point in a large-area seafloor PNT network. The combined measurement method proposed in this paper only needs to conduct circle-sailing positioning for a few control points, and the coordinates of other points can be obtained based on the results of mutual acoustic ranging. This method can achieve the same positioning accuracy as circle-sailing positioning with a shorter observation time and lower cost when the observation quality is good.

6.2. The Optimal Observation Time under Different Side Lengths of the Unit Network

Through the model given in this paper, the relationship between the number of observation sessions during acoustic ranging and the side length of the unit network under different accuracy is obtained. The more observation sessions, the longer the side length of the unit network, the smaller the density of control points, and the lower the layout cost. However, the side length will not increase infinitely with the increase of observation sessions and will gradually approach a certain value. Thereafter, increasing observation sessions has little effect on the side length, but the total observation time will still rise. The proposed model can be used to find the optimal observation time under different side lengths of the unit network, so as to improve the observation efficiency.

6.3. The Position of Carrying Out Circle-Sailing Positioning in the Seafloor PNT Network

Because the seafloor PNT network covers a large area, it is not practical to carry out circle-sailing positioning in each unit network. Based on the chain coordinate transfer strategy given in this paper, the position of the weakest point in the coordinate transfer process can be found. All that is needed is to carry out circle-sailing positioning at the position of the weakest point; the positioning accuracy of control points in the whole network can meet the requirements, which can reduce operational cost and improve observation efficiency.

Author Contributions

Conceptualization, J.M. and J.Z.; funding acquisition, J.Z.; investigation, J.M.; methodology, J.M. and J.Z.; writing—original draft, J.M.; writing—review and editing, J.M. and J.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Access to the data will be considered upon request to the authors.

Acknowledgments

We would like to thank the editor and anonymous reviewers for their valuable comments and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) The direct-path acoustic network; (b) the indirect-path acoustic network; (c) the GNSS-A network; (d) the seafloor PNT network.

Figure 2. The principle of determining the coordinates of control points in the seafloor PNT network.

Figure 3. The geometric relations between control points in the unit network of the seafloor PNT network.

Figure 4. The necessary conditions to constrain a 2D unit network in a geodetic coordinate system.

Figure 5. (a) the seafloor PNT network composed of equilateral triangle unit networks; (b) the seafloor PNT network composed of square unit networks. The number “i” in the figure represents the ith unit network in the seafloor PNT network.

Figure 6. The geometric relations of control points in the equilateral triangle unit network. The yellowish points are known points.

Figure 7. The principle of the chain coordinate transfer strategy. The yellow point is the known point. The green points are the common points between the unit network A and B. The white arrow is the direction of coordinate transfer.

Figure 8. The chain coordinate transfer in the seafloor PNT network composed of the equilateral triangle unit network. The unit network with the yellowish point is the initial network. The white double line represents the starting edge. The golden point is the weakest point.

Figure 9. (a) The topographic map of the experimental water area; the red box indicates the area where the unit network is located (Reprinted with permission from Ref. [27]. 2022, John Wiley and Sons). (b) The structure of the unit network and trajectories of circle-sailing positioning.

Figure 10. The combinations of control points. (a) Points C1, C2, C3, and C4. (b) Points C1, C2, and C3. (c) Points C1, C3, and C4. (d) Points C1, C2, and C4. (e) Points C2, C3, and C4.

Figure 11. (a) The relations between the side length of the equilateral triangle unit network and the number of observation sessions at all precision levels; (b) The relations between the total observation time in the equilateral triangle unit network and observation sessions.

Figure 12. (a) The relations between the side length of the square unit network and the number of observation sessions at all precision levels; (b) The relations between the total observation time in the square unit network and observation sessions.

Table 1

The coordinates of seafloor control points determined by circle-sailing.

Point	X (m)	Y (m)	Z (m)	σx (m)	σy (m)	σz (m)	σP (m)
C1	955.117	690.385	−60.800	0.013	0.013	0.005	0.019
C2	835.095	697.771	−60.146	0.008	0.008	0.004	0.012
C3	865.058	792.475	−60.044	0.007	0.006	0.006	0.011
C4	969.805	787.069	−60.536	0.010	0.010	0.004	0.015

Table 2

The mean and standard deviation of mutual ranging results between control points.

Point 1	Point 2	Ave. (m)	Std. (m)
C1	C2	120.208	0.017
C1	C3	136.603	3.164
C1	C4	97.812	0.003
C2	C3	99.416	0.002
C2	C4	161.528	0.028
C3	C4	104.885	0.007

Table 3

The coordinates of control points determined by the relative positioning.

Combination	Point	X (m)	Y (m)	dx (m)	dy (m)	dP (m)
[Image omitted. Please see PDF.]	C1	955.002	690.412	0.115	0.026	0.118
	C2	835.028	697.822	0.067	0.051	0.084
	C3	865.048	792.606	0.009	0.130	0.131
	C4	969.906	787.032	0.100	0.037	0.107
[Image omitted. Please see PDF.]	C1	955.247	690.098	0.131	0.287	0.315
	C2	835.249	697.863	0.155	0.092	0.180
	C3	864.702	792.823	0.355	0.348	0.497
[Image omitted. Please see PDF.]	C2	835.040	697.778	0.054	0.007	0.055
	C3	864.990	792.581	0.067	0.106	0.125
	C4	969.843	787.104	0.038	0.034	0.052
[Image omitted. Please see PDF.]	C1	955.048	690.448	0.068	0.063	0.093
	C2	835.080	697.731	0.015	0.040	0.042
	C4	969.869	787.079	0.063	0.009	0.064
[Image omitted. Please see PDF.]	C1	955.725	690.029	0.608	0.356	0.705
	C3	864.643	793.227	0.414	0.752	0.859
	C4	969.546	786.812	0.259	0.257	0.365

Table 4

The relations among the acoustic ranging capability, R, the number of observation sessions, N, and the total observation time, T.

R (km)	Unit Type	σP (m)
		0.025	0.05	0.1	0.2	0.5	1
5	Δ	20 (1.111)	5 (0.278)	2 (0.111)	1 (0.056)	1 (0.056)	1 (0.056)
	☐	44 (2.444)	11 (0.611)	3 (0.167)	1 (0.056)	1 (0.056)	1 (0.056)
10	Δ	77 (8.556)	20 (2.222)	5 (0.556)	2 (0.222)	1 (0.111)	1 (0.111)
	☐	161 (17.889)	41 (4.556)	11 (1.222)	3 (0.333)	1 (0.111)	1 (0.111)
15	Δ	172 (19.111)	43 (4.778)	11 (1.222)	3 (0.333)	1 (0.111)	1 (0.111)
	☐	356 (39.556)	89 (9.889)	23 (2.556)	6 (0.667)	1 (0.111)	1 (0.111)
20	Δ	305 (33.889)	77 (8.556)	20 (2.222)	5 (0.556)	1 (0.111)	1 (0.111)
	☐	630 (70.000)	158 (17.556)	40 (4.444)	10 (1.111)	2 (0.222)	1 (0.111)

Table 5

The position of the weakest point under different side lengths and observation sessions when the number of the equilateral triangle unit network is odd.

Side Length (km)	Observation Sessions	σP (m)
		0.05	0.1	0.2	0.5	1
5	20	5.218 (5)	19.516 (19)	76.053 (76)	471.495 (471)	1883.740 (1883)
	5	/	5.218 (5)	19.516 (19)	118.426 (118)	471.495 (471)
	2	/	/	8.128 (8)	47.799 (47)	189.042 (189)
	1	/	/	/	24.236 (24)	94.886 (94)
10	5	/	/	5.337 (5)	30.840 (30)	121.263 (121)
	2	/	/	/	12.704 (12)	48.934 (48)
	1	/	/	/	/	24.804 (24)
15	11	/	/	5.247 (5)	30.302 (30)	119.113 (119)
	3	/	/	/	8.677 (8)	32.996 (32)
	1	/	/	/	/	11.402 (11)
20	5	/	/	/	8.176 (8)	31.023 (31)
	1	/	/	/	/	6.624 (6)

Table 6

The position of the weakest point under different side lengths and observation sessions when the number of the equilateral triangle unit network is even.

Side Length (km)	Observation Sessions	σP (m)
		0.05	0.1	0.2	0.5	1
5	20	4.425 (4)	18.763 (18)	75.303 (75)	470.745 (470)	1882.990 (1882)
	5	/	4.425 (4)	18.763 (18)	117.676 (117)	470.745 (470)
	2	/	/	7.362 (7)	47.048 (47)	188.292 (188)
	1	/	/	/	23.484 (23)	94.136 (94)
10	5	/	/	4.546 (4)	30.089 (30)	120.513 (120)
	2	/	/	/	11.948 (11)	48.184 (48)
	1	/	/	/	/	24.053 (24)
15	11	/	/	4.454 (4)	29.551 (29)	118.363 (118)
	3	/	/	/	7.913 (7)	32.245 (32)
	1	/	/	/	/	10.644 (10)
20	5	/	/	/	7.410 (7)	30.272 (30)
	1	/	/	/	/	5.849 (5)

Table 7

The position of the weakest point under different side lengths and observation sessions in the square unit network.

The Side Length (km)	Observation Sessions	σP (m)
		0.05	0.1	0.2	0.5	1
5	11	/	2.196 (2)	8.785 (8)	57.536 (57)	230.143 (230)
	3	/	/	2.357 (2)	15.692 (15)	62.766 (62)
	1	/	/	/	5.231 (5)	20.922 (20)
10	11	/	/	2.357 (2)	14.731 (14)	58.923 (58)
	3	/	/	/	4.017 (4)	16.070 (16)
	1	/	/	/	/	5.357 (5)
15	6	/	/	/	3.587 (3)	14.348 (14)
	1	/	/	/	/	2.391 (2)
20	10	/	/	/	3.368 (3)	13.473 (13)
	2	/	/	/	/	2.695 (2)

Table 8

The comprehensive scores of design schemes for the seafloor PNT network at different precision levels.

Unit Type	Side Length (km)	Observation Sessions	σP (m)
			0.05	0.1	0.2	0.5	1
Δ	5	20	1.000	0.323	0.100	0.044	0.028
		5	/	0.470	0.139	0.059	0.038
		2	/	/	0.138	0.060	0.039
		1	/	/	/	0.059	0.039
	10	5	/	/	0.168	0.077	0.052
		2	/	/	/	0.078	0.054
		1	/	/	/	/	0.054
	15	11	/	/	0.148	0.070	0.048
		3	/	/	/	0.083	0.059
		1	/	/	/	/	0.059
	20	5	/	/	/	0.082	0.057
		1	/	/	/	/	0.058
☐	5	11	/	0.207	0.108	0.049	0.032
		3	/	/	0.077	0.048	0.033
		1	/	/	/	0.036	0.031
	10	11	/	/	0.122	0.067	0.047
		3	/	/	/	0.057	0.048
		1	/	/	/	/	0.039
	15	6	/	/	/	0.063	0.054
		1	/	/	/	/	0.035
	20	10	/	/	/	0.068	0.054
		2	/	/	/	/	0.042

References

1. Liu, J.; Chen, G.; Zhao, J.; Gao, K.; Liu, Y. Development and Trends of Marine Space-Time Frame Network. Geomat. Form. Sci. Wuhan Univ.; 2019; 44, pp. 17-37.

2. Yang, Y.; Xu, T.; Xue, S. Progresses and Prospects in Developing Marine Geodetic Datum and Marine Navigation of China. Acta Geod. Cartogr. Sin.; 2017; 46, pp. 1-8.

3. Spiess, F.N.; Chadwell, C.D.; Hildebrand, J.A.; Young, L.E.; Purcell, G.H.; Dragert, H. Precise GPS/Acoustic Positioning of Seafloor Reference Points for Tectonic Studies. Phys. Earth Planet. Inter.; 1998; 108, pp. 101-112. [DOI: https://dx.doi.org/10.1016/S0031-9201(98)00089-2]

4. Creed, E.L.; Glenn, S.; Schofield, O.M.; Barrier, H.; Petrecca, R.; Dobarro, J.; McLean, S.; Barnard, A.; Brown, K.; Adams, R. et al. LEO-15 Observatory—The next generation. Proceedings of the OCEANS 2005 MTS/IEEE; Washington, DC, USA, 17–23 September 2005; Volume 1, pp. 657-661.

5. Favali, P.; Beranzoli, L. Seafloor Observatory Science: A Review. Ann. Geophys.; 2006; 49, pp. 515-567. [DOI: https://dx.doi.org/10.4401/ag-3125]

6. Matsumoto, Y.; Ishikawa, T.; Fujita, M.; Sato, M.; Saito, H.; Mochizuki, M.; Yabuki, T.; Asada, A. Weak interpolate coupling beneath the sub-duction zone off Fukushima, NE Japan, inferred from GPS/Acoustic seafloor geodetic observation. Earth Planet Space; 2008; 60, pp. 9-12. [DOI: https://dx.doi.org/10.1186/BF03353114]

7. Chadwell, C.D.; Sweeney, A.D. Acoustic Ray-Trace Equations for Seafloor Geodesy. Mar. Geod.; 2010; 33, pp. 164-186. [DOI: https://dx.doi.org/10.1080/01490419.2010.492283]

8. Sweeney, A.D.; Chadwell, C.D.; Hildebrand, J.A.; Spiess, F.N. Centimeter-level positioning of seafloor acoustic transponders from a deeply towed interrogator. Mar. Geod.; 2005; 28, pp. 39-70. [DOI: https://dx.doi.org/10.1080/01490410590884502]

9. Gagnon, K.; Chadwell, C.; Norabuena, E. Measuring the onset of locking in the Peru-Chile trench with GPS and acoustic measurements. Nature; 2005; 434, 205. [DOI: https://dx.doi.org/10.1038/nature03412] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15758997]

10. Burgmann, Chadwell. Seafllor Geodesy. Annu. Rev. Earth Planet. Sci.; 2014; 42, pp. 509-534.

11. Yokota, Y.; Ishikawa, T. Shallow Slow Slip Events along the Nankai Trough Detected by GNSS-A. Sci. Adv.; 2020; 6, eaay5786. [DOI: https://dx.doi.org/10.1126/sciadv.aay5786] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31998843]

12. Kris, O. DARPA Discovers “GPS-Like” Undersea Drone Connectivity—Emerging DARPA Technology Connects Undersea Drones with Low-Frequency Sonar. Defense System. 2017; Available online: https://www.defenseone.com/defense-systems/2017/02/darpa-discovers-gps-like-undersea-drone-connectivity/192123/ (accessed on 15 September 2022).

13. Zhao, J.; Zou, Y.; Zhang, H.; Fang, S. A new method for absolute datum transfer in seafloor control network measurement. Mar. Sci. Technol.; 2016; 21, pp. 216-226. [DOI: https://dx.doi.org/10.1007/s00773-015-0344-z]

14. Chen, G.; Liu, Y.; Liu, Y.; Tian, Z.; Liu, J.; Li, M. Adjustment of Transceiver Lever Arm Offset and Sound Speed Bias for GNSS-Acoustic Positioning. Remote Sens.; 2019; 11, 1606. [DOI: https://dx.doi.org/10.3390/rs11131606]

15. Grafarend, E.W. Optimization of geodetic networks. Can. Surv.; 1974; 28, pp. 716-723. [DOI: https://dx.doi.org/10.1139/tcs-1974-0120]

16. Koch, K.R. First Order Design: Optimization of the Configuration of a Network by Introducing Small Position Changes; Springer: Berlin/Heidelberg, Germany, 1985.

17. Eshagh, M.; Kiamehr, R. A strategy for optimum designing of the geodetic networks from the cost, reliability and precision views. Acta Geod. Geoph. Hung.; 2007; 42, pp. 297-308. [DOI: https://dx.doi.org/10.1556/AGeod.42.2007.3.4]

18. Bagherbandi, M.; Eshagh, M.; Sjoberg, L.E. Multi-objective versus single-objective models in geodetic network optimiza-tion. Nord. J. Surv. Real Estate Res.; 2009; 6, pp. 7-20.

19. Horemuz, M.; Patrick, J. Optimum Establishment of Total Station. J. Surv. Eng. Asce; 2017; 143, 06016004. [DOI: https://dx.doi.org/10.1061/(ASCE)SU.1943-5428.0000207]

20. Alizadeh-Khameneh, M.A.; Horemuz, M.; Jensen, A.B.; Andersson, J.V. Optimal Vertical Placement of Total Station. J. Surv. Eng. Asce; 2018; 144, 06018001. [DOI: https://dx.doi.org/10.1061/(ASCE)SU.1943-5428.0000255]

21. Yang, Y.; Liu, Y.; Sun, D.; Xu, T.; Xue, S.; Han, Y.; Zeng, A. Seafloor geodetic network establishment and key technologies. Sci. China Earth Sci.; 2020; 63, pp. 1188-1198. [DOI: https://dx.doi.org/10.1007/s11430-019-9602-3]

22. Zhao, J. Modern Marine Surveying and Charting; Wuhan University Press: Wuhan, China, 2008.

23. Torge, W.; Müller, J. Geodesy; Walter de Gruyter: Berlin, Germany, 2012.

24. Chen, H.H.; Wang, C.C. Optimal localization of a seafloor transponder in shallow water using acoustic ranging and GPS observations. Ocean. Eng.; 2007; 34, pp. 2385-2399. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2007.05.005]

25. Leonard, J.J.; Bennett, A.A.; Smith, C.M.; Feder, H.J.S. Autonomous Underwater Vehicle Navigation; Massachusetts Institute of Technology: Cambridge, MA, USA, 1998.

26. Whitcomb, L.L.; Yoerger, D.R.; Singh, H. Combined doppler/lbl based navigation of underwater vehicles. Proceedings of the 11th International Symposium on Unmanned Untethered Submersible Technology, New England Center; Durham, NH, USA, 22–25 August 1999.

27. Chen, X.; Zhang, H.; Zhao, J. Positioning accuracy model of sailing-circle GPS-Acoustic method. Earth Space Sci.; 2021; 8, e2019EA000639. [DOI: https://dx.doi.org/10.1029/2019EA000639]