Underwater Network Time Synchronization Method Based on Probabilistic Graphical Models

Abstract

In underwater clustering and benchmark networks, nodes need to reduce the rate and energy consumption of acoustic communication while ensuring synchronization accuracy. In large-scale networks, the improvement in the efficiency of existing network time synchronization often relies on the optimization of topological structures, and the improvement in efficiency within local areas is limited. This paper proposes a method to synchronize underwater time using the probability graph model. The method utilizes the positional and motion status information of sensor networks to construct a factor graph model for distributed network synchronization. By simplifying the marginal probability density function of the system clock difference, it can quickly calculate the clock difference parameters of nodes, thereby effectively improve the synchronization efficiency. The experimental results show that the method can complete global time synchronization within a cycle while achieving a clock difference correction accuracy higher than seconds, which significantly optimized the synchronization cycle and efficiency, and reduced the energy consumption of the acoustic communication.

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 a 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 and involves a 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 their 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 has focused on achieving reliable time synchronization in high-latency and unstable network environments [3]. It includes FTSP [4], RBS [5], and DMTS [6] for wireless sensor networks. The emergence of Time-Sensitive Networking (TSN) [7,8,9] further demonstrates that the future direction of sensor network data lies in reliable and real-time communication, emphasizing efficient and reliable data transmission. A key feature of TSN is the design of redundant paths, which ensures that network traffic can continue smoothly even if one path fails, a capability that is often challenging for most centralized networks to achieve.

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 with 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. [10] 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 not suitable for networks with highly mobile nodes.

In 2017, Wang Shuo et al. [11] proposed an energy-optimal clustering method for underwater sensor networks, dividing time synchronization into inter- and intra-cluster synchronization [12] 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. [13] 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 the time synchronization. It also considered the impact of node mobility on the 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 [14,15,16] or by employing more accurate models of network node mobility [17,18,19,20]. 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 with clustering models. For example, Sun C et al. [21] proposed Multi-Hop Time Synchronization for Underwater Acoustic Networks (MSUAN) in 2012, emphasizing the reduction in synchronization errors through multilevel timestamp exchanges. Wen J et al. [22] proposed an Improved Multi-Hop Time Synchronization for Underwater Acoustic Networks (IMSUAN) in 2013, improving synchronization accuracy by establishing a time synchronization tree, optimizing the node-listening mechanisms, and using global bias thresholds to filter out inaccurate timestamps. In 2021, Shams R et al. [23] proposed an algorithm that simultaneously performs node localization and time synchronization in multi-hop underwater sensor networks using a single anchor node. The algorithm employs gradient techniques to solve unconstrained optimization problems and optimizes node localization and time synchronization through calculations of angle information and propagation delays.

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 [24]) 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 multi-information fusion and the communication pressure from distributed network structures increases the energy consumption burden of underwater networks. How to achieve rapid time synchronization in underwater networks while maintaining a certain level of energy consumption was the research objective of this study.

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 with conventional network time synchronization methods, the probability-graph-model-based time synchronization method employs distributed networking, reducing the dependence on central nodes and improving network robustness. This method improves multi-user protocols, integrates location and sound velocity information, significantly enhances the synchronization efficiency, and reduces the computational requirements for nodes through binary simplification of the algorithm.

This paper briefly introduces the application scenarios and the necessary measurement parameters of the proposed method in Section 2. Section 3 provides a detailed overview of the topologies and limitations of distributed networks suitable for probabilistic graphs, the advantages of the message-passing scheme compared with other underwater network time synchronization methods, the mathematical and physical models when the probabilistic graph method processes single-cycle and multicycle data, the binary simplification method to reduce the model complexity and improve the computational efficiency, the impacts of different parameter errors on the clock offset estimation, and the pseudocode for the calculations of the single-cycle and multi-cycle single-node models. Section 4 presents real field experiment validations that compared the improvements in multicycle data against single-cycle data, as well as performance comparisons with traditional underwater time synchronization methods under multi-cycle data conditions.

2. Application Scenarios

The underwater network time synchronization methods discussed in this document are applied in a scenario, as shown in Figure 2, where the sensor networks periodically synchronize their time reference. As shown in

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 kth interaction
$P_i\left[k\right]$	The position of the node at the kth interaction
C	The average sound velocity in the water layer where the network is located

, 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 round of synchronization can obtain the time values $T_{\mathrm{AB}}$, $T_{\mathrm{BA}}$, $T_{\mathrm{BC}}$, $T_{\mathrm{CB}}$, $T_{\mathrm{AC}}$, and $T_{\mathrm{CA}}$, and node positions $P_{\mathrm{A}}$, $P_{\mathrm{B}}$, and $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 this research paper.

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 a 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 the network synchronization.

3.1. Network Topology

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 the ordinary nodes), the network can be divided into six networks based on the communication ranges 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 multi-hop networking can significantly improve the synchronization efficiency in large-scale networks. However, these improvements are primarily based on optimizations in a topological structure and offer no enhancement in local areas (within the communication coverage of a single node) compared with inter-platform time synchronization methods.

The intra-cluster 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 between 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 methods previously used 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

Underwater network time synchronization methods based on probabilistic graphical models are 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. The mathematical model for the clock source of this method is

$$T=a\times t+\Phi .$$

(1)

where T is the local time axis, t is the standard reference time axis, a is the clock drift (clock frequency error), and $\Phi$ is the clock offset (clock phase error). Due to the ability to achieve a real-time reference alignment of time, in most cases, the model does not require a correction for clock drift a.

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(\Phi ,T,\tau ,D,P)=p(\Phi ,T,\tau )p(\tau ,D)p(D,P).$$

(2)

In Equation (2), $\Phi$ represents the set of clock offsets for each reference beacon, T denotes the set of all reception times of the measurement signal, $\tau$ represents the set of measurements of the time delay difference, D is the set of range measurements, and P is the set of underwater reference positions.

3.3.1. Single-Group Interaction Time Synchronization Model

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

$f_1$ represents the reference for the 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({\Phi }_i\mid {\Phi }_j\right)=\sum _{\sim \left\{{\Phi }_i\right\}}{f}_1\left({\Phi }_i,{\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},$$

(3)

and the expression of the edge function is

$$\begin{array}{cc}\hfill & p\left({\Phi }_{\mathrm{B}}\right)\hfill \\ \hfill & =\left(f_a\left({\Phi }_{\mathrm{B}}\mid {\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({\Phi }_{\mathrm{B}}\mid {\Phi }_{\mathrm{C}}\right)\left(f_a\left({\Phi }_{\mathrm{C}}\mid {\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}$$

(4)

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

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

(5)

${\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({\Phi }_i\mid {\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({\Phi }_i+{\Phi }_j\right)\right)}^2\right\}.$$

(6)

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

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

(7)

${\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 [25]

$$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\}.$$

(8)

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

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

(9)

${\epsilon }_p={\left[\begin{array}{cc}{\epsilon }_{px}\hfill & {\epsilon }_{py}\hfill \end{array}\right]}^{\mathrm{T}}$ represents the positional error, which can represent the distributions ${\epsilon }_{px}\sim \mathrm{N}\left(0,{\sigma }_p^2\right)$ and ${\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)$ and $\Delta y\sim \mathrm{N}\left(P_{ij}-P_{jy},{\sigma }_p^2\right)$. $P_{ix}$ and $P_{ij}$ denote the horizontal and vertical coordinates of the reference points i, respectively, while $\Delta x$ and $\Delta y$ represent the horizontal and vertical components of the distance $D_{ij}$, respectively.

Let $u={(\Delta x)}^2+{(\Delta y)}^2$ and $\lambda ={\left(P_{ix}-P_{jx}\right)}^2+{\left(P_{ij}-P_{ij}\right)}^2$; then u follows a non-central chi-square distribution with two degrees of freedom [26]:

$$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),$$

(10)

$D_{ij}=\sqrt{u}$, and 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).$$

(11)

By analogy with 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 with single-cycle interactions, a factor graph model that integrates information from multiple interactions can achieve a 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 shown in Figure 8.

$f_4$ expresses the functional relationship between the clock offset in the adjacent time interval. The multi-cycle interaction model, as shown in Figure 9, is a combination of the single-cycle model and the single-node model.

Let

$${\left.\begin{array}{cc}\hfill & f_d\left({\Phi }_{\mathrm{B}i}\right)=\sum _{\sim \left\{{\Phi }_{\mathrm{B}}\right\}}{f}_1\left({\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({\Phi }_{\mathrm{B}i}\mid {\Phi }_{\mathrm{Bj}j}\right)=\sum _{\sim \left\{{\Phi }_{\mathrm{B}}\right\}}{f}_4\left({\Phi }_{\mathrm{B}i},{\Phi }_{\mathrm{B}j}\right)\hfill \end{array}\right|}_{i,j\in \{1,2,3\cdots n\}i\ne j},$$

(12)

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

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

(13)

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

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

(14)

Here, $a_{\mathrm{B}}$ represents the clock drift 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({\Phi }_{\mathrm{B}i}\right|{\Phi }_{\mathrm{B}j})={\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{\varphi }}}exp\left\{-{\displaystyle \frac{1}{2{\sigma }_{\varphi }}}{\left({\Phi }_{\mathrm{B}j}+(i-j)a_{\mathrm{B}}-{\Phi }_{\mathrm{B}i}\right)}^2\right\}$$

(15)

3.4. Binary Simplification

Under normal circumstances, the position, distance, delay, clock offset ${\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, with each grid representing an estimated clock offset value. The final estimated result of the clock offset is obtained by the 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 the complexity, we adopted a binary mode for the 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 the binary mode, the probability density function of the 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..$$

(16)

The probability density function of the 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..$$

(17)

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

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

(18)

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

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

(19)

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

(20)

${\epsilon }_D,{\epsilon }_{\tau }$, and ${\epsilon }_{\varphi }$ represent the standard deviations of the respective probability density functions. The binary representation of the 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}.$$

(21)

3.5. Error Analysis

3.5.1. Impact of Time Measurement Error on Synchronization Accuracy

From Equations (5) and (6), it can be seen that under the condition that the time measurement values follow a normal distribution ${\epsilon }_T\sim \mathrm{N}\left(0,{\sigma }_T^2\right)$, the synchronization accuracy also follows a normal distribution ${\epsilon }_{\varphi }\sim \mathrm{N}\left(0,4{\sigma }_T^2\right)$. If a standard deviation is used as the accuracy of the estimation, then ${\epsilon }_{\varphi }=2{\sigma }_{\mathrm{T}}$.

3.5.2. Impact of Sound Speed Measurement Error on Synchronization Accuracy

The sound speed measurement error is introduced in Equation (7). Similarly, using one standard deviation as the accuracy of the measurement of the sound speed ${\epsilon }_c$, we have ${\epsilon }_c={\sigma }_c$, ${\epsilon }_c\ll c$, given that

$${\epsilon }_{\tau }=-{\displaystyle \frac{D\times {\epsilon }_c}{c^2}}.$$

(22)

Substituting into Equation (5) yields

$${\epsilon }_{\varphi }=-{\displaystyle \frac{D\times {\epsilon }_c}{c^2}}.$$

(23)

3.5.3. Impact of Position Measurement Error on Synchronization Accuracy

From Equation (10), it can be seen that u follows a noncentral chi-square distribution with two degrees of freedom; thus,

$$var\left(u\right)=2n{\sigma }_p^4.$$

(24)

Let ${\sigma }_p$ represent the accuracy of position measurement $D_{ij}=\sqrt{u}$. The error of $D_{ij}$ is a standard deviation. We have

$${\epsilon }_D=std\left(D\right)=\sqrt[4]{2n}{\sigma }_p=\sqrt{2}{\sigma }_p.$$

(25)

Given $\tau ={\displaystyle \frac{D}{c}}$, we can express

$${\epsilon }_{\tau }={\displaystyle \frac{{\epsilon }_D}{c}}={\displaystyle \frac{\sqrt{2}{\epsilon }_p}{c}}.$$

(26)

Substituting into Equation (5) gives

$${\sigma }_{\varphi }={\displaystyle \frac{\sqrt{2}{\sigma }_p}{c}}.$$

(27)

3.5.4. Comprehensive Impact

When considering the input parameter errors ${\epsilon }_p,{\epsilon }_c$, and ${\epsilon }_T$ simultaneously, substituting ${\epsilon }_p$ and ${\epsilon }_c$ into $\tau ={\displaystyle \frac{D}{c}}$ yields

$${\epsilon }_{\tau }={\displaystyle \frac{D+\sqrt{2}{\epsilon }_p}{c+{\epsilon }_c}}-{\displaystyle \frac{D}{c}}={\displaystyle \frac{\sqrt{2}{\epsilon }_p-\tau \times {\epsilon }_c}{c+{\epsilon }_c}}.$$

(28)

Now, substituting ${\epsilon }_{\tau }$ and ${\epsilon }_T$ into Equation (5) gives

$${\epsilon }_{\varphi }={\epsilon }_{\tau }+2{\epsilon }_{\mathrm{T}}={\displaystyle \frac{\sqrt{2}{\epsilon }_p-\tau \times {\epsilon }_c}{c+{\epsilon }_c}}+2{\epsilon }_{\mathrm{T}}.$$

(29)

It can be seen that the main source of error in the estimation of the clock deviation arises from the time measurement error.

3.5.5. Impact of Resolution on Synchronization Precision

The resolution mentioned in this article refers to the minimum unit of clock deviation after binary simplification, which can be expressed as

$${\epsilon }_{{\varphi }_{\mathrm{Binary}\mathrm{Simplification}}}\ge max\left(\eta ,{\epsilon }_{\varphi }\right).$$

(30)

When $\eta >{\epsilon }_{\varphi }$, the simplified deviation estimation error ${\epsilon }_{{\varphi }_{\mathrm{Binary}\mathrm{Simplification}}}$ is greater than ${\eta }_j$; when $\eta <{\epsilon }_{\varphi }$, the simplified deviation estimation error ${\epsilon }_{{\varphi }_{\mathrm{Binary}\mathrm{Simplification}}}$ is greater than or equal to the clock deviation estimation error ${\epsilon }_{\varphi }$. Therefore, the resolution should ideally be less than the accuracy of the time measurement.

3.6. Algorithm Process

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

The flow of the simplified algorithm is given in Algorithm 1.

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,{\mathit{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 ${\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 (16)

5         For all $D_{\mathrm{A}i}$, there exists a ${\tau }_{\mathrm{A}i}$ that satisfies (17)

6         For all ${\tau }_{\mathrm{A}i}$, there exists a ${\Phi }_{\mathrm{i}\mid \mathrm{A}}$ that satisfies (18)

$//$ Calculating the clock offset at node id when ${\Phi }_A=0$, path 1.2

$7{\Phi }_{j\mid \mathrm{A}}\leftarrow$ Same as steps $2\sim 6$

$//$ Calculate the clock offset of node i in set ${\Phi }_{j\mid A}$, path 1.3

$8$For ${\phi }_{k,j\mid \mathrm{A}},{\phi }_{k,j\mid \mathrm{A}}\in {\Phi }_{j\mid \mathrm{A}}$

$9{\Phi }_{k,i\mid \mathrm{j}}$ Same as steps $2\sim 6$

$//$ Calculate the final clock offset of node i

$10{\Phi }_i=\underset{k=1\sim \mathrm{lengh}h\left({\Phi }_j\right)}{AND}\left[{\Phi }_{i\mid A},{\Phi }_{k,i\mid j}\right]$

Output $\left\{{\Phi }_i\right\}$

Repeating for N cycles, the clock offset $\left\{{\Phi }_{i,1},{\Phi }_{i,2},\cdots ,{\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 by Algorithm 2.

Algorithm 2: Comprehensive multi-cycle probability graph time synchronization algorithm.

Procedure Estimated clock offset (

${\Phi }_{i,1},{\Phi }_{i,2},\cdots ,{\Phi }_{i,N}//$Multiple clock offset measurements for node i

${\sigma }_{\varphi }$ //Single-cycle clock offset estimation accuracy

)

$1\left\{{\Phi }_{i,j}\right\}\leftarrow$Select the clock offset at time j for node i $\left\{{\Phi }_i,{\Phi }_{i,2},\cdots ,{\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 ${\Phi }_{l,n}$, there exists ${\Phi }_{l,j|\mathfrak{m}}$ that satisfies 20

$4\mathbf{End}\mathbf{for}$

$//$Calculate the final clock offset at time j

$5{\Phi }_{i,j}=\underset{n=1\sim \mathbb{N}}{\mathrm{AND}}\left[\begin{array}{ccc}{\Phi }_{i,j\mathbb{N}}& \cdots & {\Phi }_{i,j\mathbb{N}}\end{array}\right]$

Output $\left\{{\Phi }_{i,j}\right\}$

4. Field Experiments

4.1. Overview of the Experiment

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

The device deployment process is illustrated in Figure 10, with the complete anchor system structure highlighted by a red box. During the deployment process, it was necessary to periodically collect the sound speed profile of the deployment area and record the coordinates of the experimental vessel as the initial values for the beacons.

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 point [0, 0] was the midpoint of the line that connects the coordinates of J1 and J4.

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

The benchmark interaction cycle was 20 s, with J1, J2, and J4 having forwarding delays of 0 s, 0.8 s, and 2.4 s, respectively. The beacon clock source was an SA.45s [27] (a chip−scale atomic clock) rubidium clock model, with a cumulative drift of less than 1 s over 135 h. Therefore, it can be assumed that the clock offset remained constant during the experiment. Assuming that each beacon used its transmission moment as the local clock’s 0 moment, the clock offset $\Phi$ for each beacon was 0 s, 0.8 s, and 2.4 s, 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, Figure 14 and Figure 15. It can be seen that the fluctuation amplitude of the delay measurements was consistently below $1\times {10}^{-4}$ s. This error was primarily due to variations in the velocity of the sound. This amplitude represented the upper limit of the accuracy of clock-offset measurement and served as the main basis for setting the resolution in discrete calculations. An excessively high resolution would affect the computation speed, while a too low resolution would reduce the synchronization accuracy. For these experimental data, the clock offset resolution was set to $1\times {10}^{-5}$ s.

4.2. Trial Conditions and Efficiency Assessment

(1) Trial conditions

The source level of the underwater transducer during sea trials was 185 dB, with $f_1=2$ kHz and $f_2=4$ kHz, and the noise level within the bandwidth was approximately 83 dB, the signal length was 1 s, the propagation loss was approximately 70 dB, the relative positioning accuracy underwater was approximately 0.1 m, and the average sound speed in water was approximately 1500 m/s.

(2) Efficiency assessment

The time synchronization technology based on probability graphs was of the type of one-way interaction, with a propagation distance of not less than $3000\mathrm{m}$, an equivalent sound speed of not less than $1500\mathrm{m}/\mathrm{s}$, a forward delay $T_r=1\mathrm{s}$, a maximum clock difference $\beta ={\Phi }_{54}=2.4\mathrm{s}$, and the number of nodes to be synchronized Y. The synchronization efficiency $\Delta T$ is given by

$$\Delta T=T_d\ge T_s+\beta \ge 3.4\mathrm{s}$$

(31)

The synchronization efficiency $\Delta T$ of the DE-Sync method under the same scenario is

$$\Delta T=Y\times \left({\displaystyle \frac{2D}{\tilde{c}}}+T_r\right)\ge Y\times 5\mathrm{s}$$

(32)

4.3. Experimental Results

Comparison Between the Single-Cycle Interaction Model and the Multi-Cycle Interaction Model

The experiment collected data from 206 cycles, although data from a single cycle were sufficient to solve the offset of the beacon clock, as shown in Figure 16 and Figure 17. With increasing redundancy of the messages, the width of the peak gradually decreased and the precision of synchronization progressively improved and approached the $1\times {10}^{-5}$ s resolution. The measurement accuracy of J2 was higher than $8\times {10}^{-4}$ s (see the yellow bar annotations in Figure 16), while the measurement accuracy of J4 exceeded $3\times {10}^{-5}$ s (see the yellow bar annotations in Figure 17).

4.4. Comparison of Probability Graph Method and Conventional Time Synchronization Under a Single Node with Multiple Cycles

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 was higher than $2.3\times {10}^{-4}$ s, while the accuracy of the DE-Sync method approached $2.59\times {10}^{-4}$ s. For the J4 beacon, the accuracy of the probabilistic graphical method was higher than $3\times {10}^{-5}$ s, while the accuracy of the DE-Sync method approached $6.93\times {10}^{-5}$ s. The clock offset estimation precision of both methods was of the same order of magnitude. However, the minimum synchronization cycle for a single beacon using the probabilistic graphical method was 3.4 s, whereas the DE method required 5 s and its cycle was limited by the propagation delay, which resulted in a 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 the synchronization efficiency while maintaining synchronization accuracy. This method solves the problem of low efficiency in existing underwater time synchronization algorithms when synchronizing the time for multiple targets. After the 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 was 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 was higher than $8\times {10}^{-4}$ s, its synchronization cycle was only limited by the maximum clock difference in the cluster and could complete the timing of the entire network within a cycle. When the accuracy of the arrival time measurement was not higher than $1\times {10}^{-4}$ s and the resolution was $1\times {10}^{-5}$ s, the estimation of the time synchronization reached the submillisecond order and the network synchronization efficiency was better than 5 s, meeting the general requirements of underwater sensor networks.