Underwater Network Time Synchronization Method Based on Probabilistic Graphical Models

1. Introduction

Network time synchronization refers to the process of synchronizing time across nodes in a network using network protocols. For surface time synchronization, it usually combines the standard time reference provided by GNSS and propagates time information through the network. Network time synchronization is widely used in data centers, distributed systems, communication networks, and other scenarios requiring time synchronization. It focuses on achieving time synchronization through network protocols in distributed computing environments, involving clock hierarchical architecture and adaptation to different network topologies.

The Network Time Protocol (NTP) is the most widely used network time protocol. It uses a hierarchical structure to organize time servers. Research in recent years has focused mainly on improving NTP security and resilience to interference, such as the work of Burbank et al. [1]. The Precision Time Protocol (PTP) is a typical centralized network protocol in which one device in the network is selected as the Master Clock and other devices act as slave Clocks, obtaining time from the Master Clock. Researchers have been working to improve the accuracy of PTP synchronization using hardware timestamps and clock filtering algorithms [2].

Network time synchronization methods are also widely applied in cloud computing and distributed databases. Research focuses on achieving reliable time synchronization in high-latency and unstable network environments [3]. It includes FTSP [4], RBS [5], and Delay Measurement et al. Time Synchronization for Wireless Sensor Networks (DMTS) [6]. The network structures of these protocols can be referenced for underwater time synchronization, but the synchronization process is not applicable due to the large propagation delay underwater. Compared to surface applications, underwater sensor nodes are usually powered by batteries, and the limited power supply severely restricts the computational and communication capabilities of the sensor nodes. Therefore, the structure of underwater networks is relatively simple and can generally be divided into three types based on network methods: centralized, distributed, and multi-hop, as shown in Figure 1.

In 2003, Ganeriwal et al. [7] proposed the Time Synchronization Protocol for Sensor Networks (TPSN) based on centralized network structures. Similarly to the NTP protocol, it classifies network nodes hierarchically and assigns each node a level number. The root node, which serves as the clock source for the entire network, is designated level 0, while other nodes synchronize with a node from the previous level, ultimately synchronizing all nodes with the root node. However, this method cannot prevent network paralysis caused by root node failure and is unsuitable for networks with highly mobile nodes.

In 2017, Wang Shuo et al. [8] proposed an energy-optimal clustering method for underwater sensor networks, dividing time synchronization into inter and intra cluster synchronization [9], and calculating the optimal number of clusters by solving mathematical expectations. This method is not suitable for highly mobile networks, as significant changes in network topology require recalculation of the optimal number of clusters and re-clustering.

In 2020, Kong Weiquan et al. [10] proposed a dual cluster head time synchronization algorithm for underwater sensor networks based on clustering. By selecting two nodes as primary and secondary cluster heads, it improved the network robustness and improved the accuracy and efficiency of time synchronization. It also considered the impact of node mobility on synchronization accuracy, reducing estimation errors to some extent by introducing advanced node mobility models. However, given the complexity of the mobility of the underwater network, this method is difficult to implement.

For time synchronization in mobile networks, the accuracy of synchronization can be improved by providing the assistance of dynamic nodes [11,12,13], or by employing more accurate models of network node mobility [14,15,16,17]. However, despite the benefits of clustering nodes in reducing the messaging overhead for time synchronization in underwater sensor networks and improving the overall network synchronization accuracy, there are still deficiencies when it comes to dealing with mobile nodes. Furthermore, in small areas, the efficiency gains from centralized networking synchronization are not significant.

Multi-hop networking offers higher robustness compared to clustering models. For example, Sun C et al. [18] proposed Multi-hop Time Synchronization for Underwater Acoustic Networks (MSUAN) in 2012, emphasizing the reduction of synchronization errors through multilevel timestamp exchanges. Subsequently, Wen J et al. [19] proposed an Improved Multi-Hop Time Synchronization for Underwater Acoustic Networks (IMSUAN) in 2013, improving synchronization accuracy by establishing a time synchronization tree and optimizing node listening mechanisms, and using global bias thresholds to filter out inaccurate timestamps.

However, both clustering models and multi-hop network models are more suitable for large-scale network nodes, i.e., when node communication distances cannot cover the entire network. In small-scale networks (where node communication distances can cover the network), both models are not different from centralized networks and are limited by central nodes. Moreover, the actual synchronization algorithms do not differ from the time synchronization algorithms between two points(such as DE-Sync [22]) and do not fully utilize the network information. In most cases, underwater sensor nodes, whether mobile submersibles or stationary underwater beacons, have some ability to acquire location information. Integrating location information with sound velocity information can greatly simplify the time-synchronization process. Higher location information accuracy can lead to higher time synchronization accuracy. However, the computational complexity brought about by multiinformation fusion and the communication pressure from distributed network structures increase the energy consumption burden of underwater networks. How to achieve rapid time synchronization in underwater networks while maintaining a certain level of energy consumption is the research objective of this paper.

To address these issues, this paper proposes a time synchronization method based on a probability graph model for underwater networks. This method solves the marginal probability density function of the system clock offset parameters and performs binary simplification to quickly calculate the clock offset parameters of various sensors within the local network, achieving dynamic time synchronization in distributed networks. Compared to conventional network time synchronization methods, the probability graph model-based time synchronization method employs distributed networking, reducing dependence on central nodes and improving network robustness. This method improves multi-user protocols, integrates location and sound velocity information, significantly enhancing synchronization efficiency, and reduces computational requirements for nodes through binary simplification of the algorithm.

2. Application Scenarios

The underwater network time synchronization methods discussed in this paper is applied in a scenario where sensor networks periodically synchronize their time reference. When node A receives a signal from node B, it records the reception time $T_{\mathrm{AB}}$ and its own position $P_{\mathrm{A}}$ Similarly, each synchronization round can obtain time values: $T_{\mathrm{AB}}$, $T_{\mathrm{BA}}$, $T_{\mathrm{BC}}$, $T_{\mathrm{CB}}$, $T_{\mathrm{AC}}$, $T_{\mathrm{CA}}$ and node positions: $P_{\mathrm{A}}$ , $P_{\mathrm{B}}$ , $P_{\mathrm{C}}$ , which constitute one cycle of time synchronization for the sensor network. In addition,the network must measure the average sound velocity C in the water layer where it is located. The above describes a section of a research paper.

Figure 2. Schematic diagram of underwater network time synchronization scenario.

Figure 2. Schematic diagram of underwater network time synchronization scenario.

Table 1. Measurement parameters table for network time synchronization Methods.

Parameters	Explanations
$T_{ij}\left[k\right]$	The timestamp at which node i receives the synchronization message from node j during the $k_{th}$ interaction
$P_i\left[k\right]$	The position of node in at the $k_{th}$ interaction
C	The average sound velocity in the water layer where the network is located

¹ Tables may have a footer.

Table 1. Measurement parameters table for network time synchronization Methods.

Parameters	Explanations
$T_{ij}\left[k\right]$	The timestamp at which node i receives the synchronization message from node j during the $k_{th}$ interaction
$P_i\left[k\right]$	The position of node in at the $k_{th}$ interaction
C	The average sound velocity in the water layer where the network is located

¹ Tables may have a footer.

3. Principles of Underwater Network Time Synchronization Methods

To simplify the existing network time synchronization protocol and increase synchronization efficiency, additional information is clearly needed as support. Information fusion would increase the computational burden on nodes, while reducing the central node’s constraints on the network. Adopting a distributed network would require further network communication resources. This paper attempts to utilize factor graph models to simplify the entire network time synchronization process, thereby reducing the communication and computational demands of network synchronization.

3.1. Network Configuration Methods

The probabilistic graph method adopts a distributed hierarchical network approach. Taking the underwater network in Figure 3 as an example (it is generally assumed that the communication range of the cluster head nodes is greater than that of ordinary nodes), the network can be divided into 6 networks based on the communication range of the nodes. Each network forms its own distributed network, while the five cluster head nodes form a larger distributed network. In this case, each node is not a central node for its local network, and the cluster head nodes only serve as time-reference standards for the local network rather than initiators of all synchronization activities. When a cluster head node fails, the distributed network can use the time of any node within the network as a reference, achieving unification of the local network.

The scope of hierarchical distributed networking is limited by the communication range of the cluster head nodes. For larger-scale networking, the communication distance of the cluster head nodes cannot cover the entire underwater acoustic network, which requires multi-hop-assisted networking. As shown in Figure 4, distributed networking is used within clusters, while multi-hop communication is employed between clusters.

3.2. Message Passing Scheme

Cluster-based networking and multihop networking can significantly improve synchronization efficiency in large-scale networks. However, these improvements are primarily based on optimizations in topological structure and offer no enhancement in local areas (within the communication coverage of a single node) compared to inter-platform time-synchronization methods.

The intracluster information transmission scheme (multiuser protocol) is shown in Figure 5. It is evident that existing network time-synchronization methods, while needing to consider information differentiation among multiple users, are fundamentally not different from one-way or two-way time-synchronization methods. They do not fully utilize the additional information available in sensor networks, thus missing the opportunity to simplify the time-synchronization process in underwater sensor networks.

The probability graph method for time synchronization, as illustrated in Figure 6, differs from previously used methods that require multiple interaction cycles to complete a single time synchronization. In contrast, the probability graph method can complete a time synchronization in each interaction cycle ($\Delta T$). In addition, multiple sets of interactions can further improve the accuracy of the synchronization.

3.3. Factor Graph Models

3.3.1. Single-Group Interaction Time Synchronization Model

This model is generally applied to underwater sensor networks (clusters, reference networks, etc.) that require high refresh rates. In this scenario, the reference points have certain initial position values and the unification of the time reference can be achieved through relative ranging. Due to the ability to achieve real-time time reference alignment, there is no need to correct clock drift.

The fusion parameters include reference position information, reference clock offset, relative ranging information, time delay difference measurement information, and relative ranging signal reception moments. At this point, the global function is represented as:

$$p(\mathsf{\Phi },T,\tau ,D,P)=p(\mathsf{\Phi },T,\tau )p(\tau ,D)p(D,P)$$

(1)

In 1, represents the set of clock offset for each reference beacon, denotes the set of all measurement signal reception times, represents the set of time delay difference measurements, is the set of range measurements, and is the set of underwater reference positions.

Taking the inter-beacon ranging among three beacons A, B, and C (with A as the reference beacon, ${\mathsf{\Phi }}_{\mathrm{A}}=0)$ as an example, the factor graph model as shown in the Figure 7:

$f_1$ represents the reference among signal reception time, time delay and reference clock offset. $f_2$ represents the functional relationship between the information from the range and the time delay difference. $f_3$ represents the functional relationship between the information from the range and the reference position.

Let:

$${\left.\begin{array}{cc}\hfill & f_a\left({\mathsf{\Phi }}_i\mid {\mathsf{\Phi }}_j\right)=\sum _{\sim \left\{{\mathsf{\Phi }}_i\right\}}{f}_1\left({\mathsf{\Phi }}_i,{\mathsf{\Phi }}_j,T_{ij},T_{ji},{\tau }_{ij}\right)\hfill \\ \hfill & f_b\left({\tau }_{ij}\right)=\sum _{\sim \left\{{\tau }_{ij}\right\}}{f}_2\left({\tau }_{ij},D_{ij}\right)\hfill \\ \hfill & f_c\left(D_{ij}\right)=\sum _{\sim \left\{D_{ij}\right\}}{f}_3\left(P_i,P_j,D_{ij}\right)\hfill \end{array}\right|}_{i,j\in \{\mathrm{A},\mathrm{B},\mathrm{C}\}i\ne j}$$

(2)

The expression of the edge function

$$\begin{array}{cc}\hfill & p\left({\mathsf{\Phi }}_{\mathrm{B}}\right)\hfill \\ \hfill & =\left(f_a\left({\mathsf{\Phi }}_{\mathrm{B}}\mid {\mathsf{\Phi }}_{\mathrm{A}}\right)f_b\left({\tau }_{\mathrm{AB}}\right)f_c\left(D_{\mathrm{AB}}\right)\right)\hfill \\ \hfill & \left(f_a\left({\mathsf{\Phi }}_{\mathrm{B}}\mid {\mathsf{\Phi }}_{\mathrm{C}}\right)\left(f_a\left({\mathsf{\Phi }}_{\mathrm{C}}\mid {\mathsf{\Phi }}_{\mathrm{A}}\right)f_b\left({\tau }_{\mathrm{AC}}\right)f_c\left(D_{\mathrm{AC}}\right)\right)f_b\left({\tau }_{\mathrm{BC}}\right)f_c\left(D_{\mathrm{BC}}\right)\right)\hfill \end{array}$$

(3)

The relationship between clock offset $\mathsf{\Phi }$, reception time T, and measured delay $\tau$ is

$${\mathsf{\Phi }}_i=T_{ij}+T_{ji}+2{\epsilon }_T-{\tau }_{ij}-{\mathsf{\Phi }}_j$$

(4)

${\epsilon }_T$ represents the measurement error of time, which follows ${\epsilon }_T\sim \mathrm{N}\left(0,{\sigma }_T^2\right)$. The function $f_1$ can be expressed as:

$$f_1\left({\mathsf{\Phi }}_i\mid {\mathsf{\Phi }}_j\right)={\displaystyle \frac{1}{2\sqrt{2\pi }{\sigma }_T}}exp\left\{-{\displaystyle \frac{1}{8{\sigma }_T^2}}{\left(T_{ij}+T_{ji}-2{\tau }_{ij}-\left({\mathsf{\Phi }}_i+{\mathsf{\Phi }}_j\right)\right)}^2\right\}$$

(5)

The relationship between range information D and time delay $\tau$ is

$${\tau }_{ij}={\displaystyle \frac{D_{ij}}{c+{\epsilon }_c}}$$

(6)

${\epsilon }_c$ represents the measurement error of the sound velocity; it follows that ${\epsilon }_c\sim \mathrm{N}\left(0,{\sigma }_c^2\right)$. Therefore, the function $f_2$ can be expressed as:

$$f_2\left({\tau }_{ij}\right)={\displaystyle \frac{\left|D_{ij}\right|}{\sqrt{2\pi }{\sigma }_c{\tau }_{ij}^2}}exp\left\{-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{D_{ij}}{{\tau }_{ij}}}-c\right)}^2/{\sigma }_c^2\right\}$$

(7)

The relationship between the reference position information P and the range information D :

$$D_{ij}=∥P_i-P_j+2{\epsilon }_p∥$$

(8)

${\epsilon }_p={\left[\begin{array}{cc}{\epsilon }_{px}\hfill & {\epsilon }_{py}\hfill \end{array}\right]}^{\mathrm{T}}$ represents the positional error, which can be represented by a distribution of ${\epsilon }_{px}\sim \mathrm{N}\left(0,{\sigma }_p^2\right),{\epsilon }_{py}\sim \mathrm{N}\left(0,{\sigma }_p^2\right)$. Then $\Delta x\sim \mathrm{N}\left(P_{ix}-P_{jx},{\sigma }_p^2\right)$, $\Delta y\sim \mathrm{N}\left(P_{ij}-P_{jy},{\sigma }_p^2\right).P_{ix},P_{ij}$, respectively, denote the horizontal and vertical coordinates of the reference point $i,\Delta x$ and $\Delta y$ represent the horizontal and vertical components of the distance $D_{ij}$.

Let $u={\left(\Delta x\right)}^2+{\left(\Delta y\right)}^2,\lambda ={\left(P_{ix}-P_{jx}\right)}^2+{\left(P_{ij}-P_{ij}\right)}^2$ then u follows a noncentral chi-square distribution with 2 degrees of freedom:

$$p\left(u\right)={\displaystyle \frac{1}{2{\sigma }_p^2}}{\left({\displaystyle \frac{u}{\lambda }}\right)}^{{\displaystyle \frac{n}{4}}-{\displaystyle \frac{1}{2}}}exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{u^{\frac{1}{2}}-{\lambda }^{\frac{1}{2}}}{{\sigma }_p}}\right)}^2\right)I_{{\displaystyle \frac{n}{2}}-1}\left({\displaystyle \frac{\sqrt{u\lambda }}{{\sigma }_p^2}}\right)$$

(9)

$D_{ij}=\sqrt{u}$, the function $f_c$ can be expressed as

$$f_3\left(D_{ij}\right)={\displaystyle \frac{\left|D_{ij}\right|}{{\sigma }_p^2}}exp\left(-{\displaystyle \frac{D_{ij}^2+\lambda }{2{\sigma }_p^2}}\right)I_0\left({\displaystyle \frac{D_{ij}\sqrt{\lambda }}{{\sigma }_p^2}}\right)$$

(10)

By analogy to the derivation of the reference B, the marginal probability density function for the clock offset of reference C can be obtained similarly.

3.3.2. Comprehensive Time Synchronization Model for Multiple Group Interactions

Compared to single-cycle interactions, a factor graph model that integrates information from multiple interactions can achieve higher accuracy of clock offset estimation. Taking multiple interactions between reference A and mobile platform B as an example, the single-node time synchronization factor graph model is as follows:

Figure 8. Factor graph model for single node underwater network time synchronization.

Figure 8. Factor graph model for single node underwater network time synchronization.

$f_4$ express the functional relationship between adjacent time interval clock offset. The multicycle interaction model is a combination of the single-cycle model and the single-node model.

Figure 9. Factor graph model for multi-cycle underwater network time synchronization.

Figure 9. Factor graph model for multi-cycle underwater network time synchronization.

Let

$${\left.\begin{array}{cc}\hfill & f_d\left({\mathsf{\Phi }}_{\mathrm{B}i}\right)=\sum _{\sim \left\{{\mathsf{\Phi }}_{\mathrm{B}}\right\}}{f}_1\left({\mathsf{\Phi }}_{\mathrm{B}i},T_i,{\tau }_i\right)\hfill \\ \hfill & f_e\left({\tau }_i\right)=\sum _{\sim \left\{{\tau }_i\right\}}{f}_2\left({\tau }_i,D_i\right)\hfill \\ \hfill & f_f\left(D_i\right)=\sum _{\sim \left\{D_i\right\}}{f}_3\left(P_{\mathrm{A}},P_{\mathrm{Bi}},D_i\right)\hfill \\ \hfill & f_g\left({\mathsf{\Phi }}_{\mathrm{B}i}\mid {\mathsf{\Phi }}_{\mathrm{Bj}j}\right)=\sum _{\sim \left\{{\mathsf{\Phi }}_{\mathrm{B}}\right\}}{f}_4\left({\mathsf{\Phi }}_{\mathrm{B}i},{\mathsf{\Phi }}_{\mathrm{B}j}\right)\hfill \end{array}\right|}_{i,j\in \{1,2,3\cdots n\}i\ne j}$$

(11)

n is the maximum number of interactions. Taking ’${\mathsf{\Phi }}_{\mathrm{B}2}$’ as an example, the expression for its marginal function is:

$$p\left({\mathsf{\Phi }}_{\mathrm{B}2}\right)=\left(f_g\left({\mathsf{\Phi }}_{\mathrm{B}2}\mid {\mathsf{\Phi }}_{\mathrm{B}1}\right)f_g\left({\mathsf{\Phi }}_{\mathrm{B}2}\mid {\mathsf{\Phi }}_{\mathrm{B}3}\right)\right)\prod _{i=1}^3{f}_d\left({\mathsf{\Phi }}_{\mathrm{B}i}\right)f_e\left({\tau }_i\right)f_f\left(D_{\mathrm{i}}\right)$$

(12)

The relationship between the clock offset at time i (denoted as ${\mathsf{\Phi }}_{\mathrm{Bi}}$ ) and the clock offset at time j (denoted as $\left.{\mathsf{\Phi }}_{\mathrm{B}j}\right)$ is:

$${\mathsf{\Phi }}_{\mathrm{B}i}={\mathsf{\Phi }}_{\mathrm{B}j}+(i-j)a_{\mathrm{B}}+{\epsilon }_{\varphi }$$

(13)

Among them, $a_{\mathrm{B}}$ represents the skew in clock frequency and ${\epsilon }_{\varphi }$ represents the jitter error in the clock phase, which follows ${\epsilon }_{\varphi }\sim \mathrm{N}\left(0,{\sigma }_{\varphi }^2\right)$.

$$f_4\left({\mathsf{\Phi }}_{\mathrm{B}i}\right|{\mathsf{\Phi }}_{\mathrm{B}j})={\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{\varphi }}}exp\left\{-{\displaystyle \frac{1}{2{\sigma }_{\varphi }}}{\left({\mathsf{\Phi }}_{\mathrm{B}j}+(i-j)a_{\mathrm{B}}-{\mathsf{\Phi }}_{\mathrm{B}i}\right)}^2\right\}$$

(14)

3.4. Binary Simplification

Under normal circumstances, the position, distance, delay, clock offset ${\mathsf{\Phi }}_i$, and other information about the target to be estimated are distributed continuously, and the correlation functions between variables are also distributed continuously. Directly estimating the clock offset by calculating the expectation of the continuous distribution function has a high computational complexity. Therefore, the sampling concept is introduced to discretize the continuous distribution function. The most probable result of the clock offset distribution is divided into grids, each grid representing an estimated clock offset value. The final estimated result of the clock offset is obtained by weighted averaging of the estimated values for each grid. The main steps for calculating the weights are as follows. If the measured data set $\varrho$ is known, fixed-interval scattered points can be used to represent the grid area $(X=x,Y=y).$. The weight $w=f(X=x,Y=y\mid Q)$ of this grid area is represented by the value of the probability density function at the position $x,y$. The calculation of the factor graph mainly focuses on the transmission and updating of messages. To reduce complexity, we adopt a binary mode for message transmission. The binary mode transforms the sum operations in the message transmission process into logical OR operations and product operations into logical AND operations, greatly reducing the computational load.

Based on binary mode, the probability density function of distance information can be simplified to

$$f_3\left(D\right)\Rightarrow \left\{\begin{array}{cc}1,\hfill & -{\epsilon }_D<D-{\lambda }^{{\displaystyle \frac{1}{2}}}<{\epsilon }_D\hfill \\ 0,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(15)

The probability density function of delay information can be simplified to

$$f_2\left(\tau \right)\Rightarrow \left\{\begin{array}{cc}1,\hfill & -{\epsilon }_{\tau }<{\displaystyle \frac{D}{\tau }}-c<{\epsilon }_{\tau }\hfill \\ 0,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(16)

Under single-group interaction, the probability density function of node clock offset can be simplified to

$$f_1\left({\mathsf{\Phi }}_i\right|{\mathsf{\Phi }}_j)\Rightarrow \left\{\begin{array}{cc}1,\hfill & -{\epsilon }_{\varphi }<T_{ij}+T_{ji}-2{\tau }_{ij}-\left({\mathsf{\Phi }}_i+{\mathsf{\Phi }}_j\right)<{\epsilon }_{\varphi }\hfill \\ 0,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(17)

The probability density function of clock offset between nodes under multiple interactions can be separately simplified to

$$f_1\left({\mathsf{\Phi }}_{\mathrm{B}i}\right)\Rightarrow \left\{\begin{array}{cc}1,\hfill & -{\epsilon }_{\varphi }<T_i-{\tau }_i-{\mathsf{\Phi }}_{\mathrm{B}i}<{\epsilon }_{\varphi }\hfill \\ 0,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(18)

$$f_4\left({\mathsf{\Phi }}_{\mathrm{B}i}\right|{\mathsf{\Phi }}_{\mathrm{B}j})\Rightarrow \left\{\begin{array}{cc}1,\hfill & -{\epsilon }_{\varphi }<{\mathsf{\Phi }}_{\mathrm{B}i}-{\mathsf{\Phi }}_{\mathrm{B}j}-(i-j)a_{\mathrm{B}}<{\epsilon }_{\varphi }\hfill \\ 0,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(19)

${\epsilon }_D,{\epsilon }_{\tau },{\epsilon }_{\varphi }$ represent the standard deviations of the respective probability density functions.The binary representation of message transmission is

$$\begin{array}{cc}& u_{x-f}\left(x\right)=\underset{h\in n\left(x\right)\setminus \left\{f\right\}}{AND}\left[u_{h-x}\left(x\right)\right]\hfill \\ & u_{f-x}\left(x\right)=\underset{\{\sim \}x}{OR}\left[f\left(X\right)\underset{y\in n\left(f\right)\setminus \left\{x\right\}}{AND}\left[u_{y-f}\left(y\right)\right]\right]\hfill \end{array}$$

(20)

3.5. Algorithm Process

This section focuses on the time synchronization scenario between underwater benchmarks. After establishing a factor graph model, the sum-product algorithm is employed to calculate the marginal probability function. Finally, to simplify the computational process, discretization and binary representation are adopted for the calculations.

The simplified algorithm flow is as follows:

Table 2. Pseudocode of Single-Cycle Probability Graph Method

Algorithm 1: Single-cycle probability graph time synchronization algorithm
Procedure Estimated clock offset(
$p_{\mathrm{A}},p_1,\cdots ,p_{\mathrm{n}}//$ Reference node A and n nodes to be synchronized
$c,{\mathbf{T}}_{(n+1)\times n}//$ Sound velocity and reciprocal range observations
${\sigma }_p,{\sigma }_c,{\sigma }_T//$ Position, speed of sound, time instant observed value accuracy
)
$1\left\{p_i,p_j\right\}\leftarrow$ Select any two points $\left\{p_1,\cdots ,p_2\right\}$
$//$ Calculate the clock offset of node i at ${\mathsf{\Phi }}_A=0$,path 1.1
$2{\sigma }_d\leftarrow \int f_c\left(D_{\mathrm{A}i}\right)D_{\mathrm{A}i}^{2}d{D}_{\mathrm{A}i}-\int f_c\left(D_{\mathrm{A}i}\right)D_{\mathrm{A}i}d{D}_{\mathrm{A}i}$
$3{\sigma }_{\tau },{\sigma }_{\oplus }\leftarrow$ Similarly
4         For all $\left\{p_{\mathrm{A}},p_i\right\}$, there exists a $D_{\mathrm{A}i}$ that satisfies 15
5         For all $D_{\mathrm{A}i}$, there exists a ${\tau }_{\mathrm{A}i}$ that satisfies 16
6         For all ${\tau }_{\mathrm{A}i}$, there exists a ${\mathsf{\Phi }}_{\mathrm{i}\mid \mathrm{A}}$ that satisfies 17
$//$ Calculating the clock offset at node id when ${\mathsf{\Phi }}_A=0$, Path 1.2
$7{\mathsf{\Phi }}_{j\mid \mathrm{A}}\leftarrow$ Same as steps $2\sim 6$
$//$ Calculate the clock offset of node i in set ${\mathsf{\Phi }}_{j\mid A}$, path 1.3
$8$For ${\phi }_{k,j\mid \mathrm{A}},{\phi }_{k,j\mid \mathrm{A}}\in {\mathsf{\Phi }}_{j\mid \mathrm{A}}$
$9{\mathsf{\Phi }}_{k,i\mid \mathrm{j}}$ Same as steps $2\sim 6$
$//$ Calculate the final clock offset of node i
$10{\mathsf{\Phi }}_i=\underset{k=1\sim \mathrm{lengh}h\left({\mathsf{\Phi }}_j\right)}{AND}\left[{\mathsf{\Phi }}_{i\mid A},{\mathsf{\Phi }}_{k,i\mid j}\right]$
Output $\left\{{\mathsf{\Phi }}_i\right\}$

Table 2. Pseudocode of Single-Cycle Probability Graph Method

Algorithm 1: Single-cycle probability graph time synchronization algorithm
Procedure Estimated clock offset(
$p_{\mathrm{A}},p_1,\cdots ,p_{\mathrm{n}}//$ Reference node A and n nodes to be synchronized
$c,{\mathbf{T}}_{(n+1)\times n}//$ Sound velocity and reciprocal range observations
${\sigma }_p,{\sigma }_c,{\sigma }_T//$ Position, speed of sound, time instant observed value accuracy
)
$1\left\{p_i,p_j\right\}\leftarrow$ Select any two points $\left\{p_1,\cdots ,p_2\right\}$
$//$ Calculate the clock offset of node i at ${\mathsf{\Phi }}_A=0$,path 1.1
$2{\sigma }_d\leftarrow \int f_c\left(D_{\mathrm{A}i}\right)D_{\mathrm{A}i}^{2}d{D}_{\mathrm{A}i}-\int f_c\left(D_{\mathrm{A}i}\right)D_{\mathrm{A}i}d{D}_{\mathrm{A}i}$
$3{\sigma }_{\tau },{\sigma }_{\oplus }\leftarrow$ Similarly
4         For all $\left\{p_{\mathrm{A}},p_i\right\}$, there exists a $D_{\mathrm{A}i}$ that satisfies 15
5         For all $D_{\mathrm{A}i}$, there exists a ${\tau }_{\mathrm{A}i}$ that satisfies 16
6         For all ${\tau }_{\mathrm{A}i}$, there exists a ${\mathsf{\Phi }}_{\mathrm{i}\mid \mathrm{A}}$ that satisfies 17
$//$ Calculating the clock offset at node id when ${\mathsf{\Phi }}_A=0$, Path 1.2
$7{\mathsf{\Phi }}_{j\mid \mathrm{A}}\leftarrow$ Same as steps $2\sim 6$
$//$ Calculate the clock offset of node i in set ${\mathsf{\Phi }}_{j\mid A}$, path 1.3
$8$For ${\phi }_{k,j\mid \mathrm{A}},{\phi }_{k,j\mid \mathrm{A}}\in {\mathsf{\Phi }}_{j\mid \mathrm{A}}$
$9{\mathsf{\Phi }}_{k,i\mid \mathrm{j}}$ Same as steps $2\sim 6$
$//$ Calculate the final clock offset of node i
$10{\mathsf{\Phi }}_i=\underset{k=1\sim \mathrm{lengh}h\left({\mathsf{\Phi }}_j\right)}{AND}\left[{\mathsf{\Phi }}_{i\mid A},{\mathsf{\Phi }}_{k,i\mid j}\right]$
Output $\left\{{\mathsf{\Phi }}_i\right\}$

Repeating for N cycles, the clock offset $\left\{{\mathsf{\Phi }}_{i,1},{\mathsf{\Phi }}_{i,2},\cdots ,{\mathsf{\Phi }}_{i,N}\right\}$ can be obtained for a single node, and by combining the single-node factor graph model, the estimated value of the comprehensive multicycle clock offset can be obtained.

Table 3. Pseudocode of Multi-Cycle Probability Graph Method

Algorithm 2: Comprehensive multi-cycle probability graph time synchronization algorithm
Procedure Estimated clock offset(
${\mathsf{\Phi }}_{i,1},{\mathsf{\Phi }}_{i,2},\cdots ,{\mathsf{\Phi }}_{i,N}//$Multiple clock offset measurements for node i
${\sigma }_{\varphi }$ //Single-cycle clock offset estimation accuracy
)
$1\left\{{\mathsf{\Phi }}_{i,j}\right\}\leftarrow$Select the clock offset at time j for node i $\left\{{\mathsf{\Phi }}_i,{\mathsf{\Phi }}_{i,2},\cdots ,{\mathsf{\Phi }}_{i,N}\right\}$
$//$ Calculate the clock offset at time j for node i, following path 2.1 to end
$2\mathbf{For}\mathbf{n}=\mathbf{1}\mathbf{to}\mathbf{N}$
3         For all ${\mathsf{\Phi }}_{l,n}$, there exists ${\mathsf{\Phi }}_{l,j|\mathfrak{m}}$ that satisfies the 19
$4\mathbf{End}\mathbf{for}$
$//$Calculate the final clock offset at time j
$5{\mathsf{\Phi }}_{i,j}=\underset{n=1\sim \mathbb{N}}{\mathrm{AND}}\left[\begin{array}{ccc}{\mathsf{\Phi }}_{i,j\mathbb{N}}& \cdots & {\mathsf{\Phi }}_{i,j\mathbb{N}}\end{array}\right]$
Output $\left\{{\mathsf{\Phi }}_{i,j}\right\}$

Table 3. Pseudocode of Multi-Cycle Probability Graph Method

Algorithm 2: Comprehensive multi-cycle probability graph time synchronization algorithm
Procedure Estimated clock offset(
${\mathsf{\Phi }}_{i,1},{\mathsf{\Phi }}_{i,2},\cdots ,{\mathsf{\Phi }}_{i,N}//$Multiple clock offset measurements for node i
${\sigma }_{\varphi }$ //Single-cycle clock offset estimation accuracy
)
$1\left\{{\mathsf{\Phi }}_{i,j}\right\}\leftarrow$Select the clock offset at time j for node i $\left\{{\mathsf{\Phi }}_i,{\mathsf{\Phi }}_{i,2},\cdots ,{\mathsf{\Phi }}_{i,N}\right\}$
$//$ Calculate the clock offset at time j for node i, following path 2.1 to end
$2\mathbf{For}\mathbf{n}=\mathbf{1}\mathbf{to}\mathbf{N}$
3         For all ${\mathsf{\Phi }}_{l,n}$, there exists ${\mathsf{\Phi }}_{l,j|\mathfrak{m}}$ that satisfies the 19
$4\mathbf{End}\mathbf{for}$
$//$Calculate the final clock offset at time j
$5{\mathsf{\Phi }}_{i,j}=\underset{n=1\sim \mathbb{N}}{\mathrm{AND}}\left[\begin{array}{ccc}{\mathsf{\Phi }}_{i,j\mathbb{N}}& \cdots & {\mathsf{\Phi }}_{i,j\mathbb{N}}\end{array}\right]$
Output $\left\{{\mathsf{\Phi }}_{i,j}\right\}$

4. Field Experiments

4.1. Overview of the Experiment

On 3 April 2023, at 10:58 AM, an underwater cluster time synchronization experiment was conducted in the South China Sea at a depth of 3400 meters. The experiment involved the deployment and calibration of four seafloor beacon transponders, with a total duration of 3 hours.

Figure 10. experimental site

Figure 10. experimental site

The frequency band of the communication signal ranged from 2 to 4 kHz. The relative positions of the deployed beacons are shown in Figure 11 (Beacon J3 was excluded from subsequent data processing due to battery depletion and incomplete data collection).

The sound velocity profile was collected on April 3, 2023, at 19:10:09. The results of the acquisition of the sound velocity are shown in Figure 12.

The benchmark inter-interaction cycle is 20s, with J1, J2, and J4 having forwarding delays of 0s, 0.8s, and 2.4s respectively. The beacon clock source is an SA.45s rubidium clock model, with a cumulative drift of less than s over 135 hours. Therefore, it can be assumed that the clock offset remains constant during the experiment. Assuming that each beacon uses its transmission moment as the local clock’s 0 moment, the clock offset $\mathsf{\Phi }$ for each beacon are 0s, 0.8s, and 2.4s, respectively.

The experiment collected 206 cycles of intermeasurement signals with a 10 s cycle. Each beacon recorded the arrival times of the signals from other beacons, using its own transmission moment as a reference, as shown in Figure 13 to Figure 15. It can be seen that the fluctuation amplitude of the delay measurements is consistently below $1\times {10}^{-4}s$. This error is primarily due to variations in the velocity of the sound. This amplitude represents the upper limit of the accuracy of clock-offset measurement and serves as the main basis for setting the resolution in discrete calculations. Excessively high resolution would affect computation speed, while too low resolution would reduce synchronization accuracy. For these experimental data, the clock offset resolution was set to $1\times {10}^{-5}s$.

Figure 13. Measurement timestamp of J1.

Figure 13. Measurement timestamp of J1.

Figure 14. Measurement timestamp of J2.

Figure 14. Measurement timestamp of J2.

Figure 15. Measurement timestamp of J4.

Figure 15. Measurement timestamp of J4.

4.2. Experimental Results

Comparison between the single-cycle interaction model and the multicycle interaction model.

The experiment collected data from 206 cycles, although data from a single cycle is sufficient to solve for the beacon clock offset, as shown in Figure 16 andFigure 17. As redundant messages increase, the width of the peak gradually decreases, and the precision of synchronization progressively improves, approaching the resolution $8\times {10}^{-4}s$. The measurement accuracy of J2 is higher than s, while the measurement accuracy of J4 exceeds s. With increasing redundant messages, the peak width gradually narrows, and the synchronization precision progressively improves, convergent towards the resolution $7\times {10}^{-5}s$.

Comparison of Probability Graph Method and Conventional Time Synchronization under Single Node

The comparison of the estimation results between the probabilistic graphical method and the Conventional Synchronization Methods (taking DE-Sync as an example) at different interaction frequencies is shown in Figure 18 and Figure 19. For the J2 beacon, the accuracy of the probabilistic graphical method is higher than $2.3\times {10}^{-4}s$, while the accuracy of the DE-Sync method approaches $2.59\times {10}^{-4}s$. For the J4 beacon, the accuracy of the probabilistic graphical method is higher than $3\times {10}^{-5}s$, while the accuracy of the DE-Sync method approaches$7.4\times {10}^{-5}s$. The clock offset estimation precision of both methods is of the same order of magnitude. However, the minimum synchronization cycle for a single beacon using the probabilistic graphical method is 3.4s, whereas the DE method requires 5 s and its cycle is limited by propagation delay, resulting in lower synchronization efficiency.

5. Discussion

This paper innovatively proposes an underwater network time synchronization method suitable for cluster targets, making full use of the known physical information between underwater sensor networks and greatly improving synchronization efficiency while maintaining synchronization accuracy. This method solves the problem of low efficiency in existing underwater time-synchronization algorithms when synchronizing time for multiple targets. After completion of the construction of the factor graph model for the known physical information between underwater sensor networks, the edge probability density function of the clock difference is quickly calculated to solve the clock difference parameters between units in the underwater sensor network. The experimental results show that, on the premise that its time synchronization accuracy is higher than $8\times {10}^{-4}s$, its synchronization period is only limited by the maximum clock difference in the cluster and can complete the timing of the entire network within a cycle. When the accuracy of arrival time measurement is not higher than $1\times {10}^{-4}s$ and the resolution is $1\times {10}^{-5}s$, the estimation of time synchronization reaches the submillisecond order and the network synchronization efficiency is better than $5s$, meeting the general requirements of underwater sensor networks.