Research on Time Synchronization Method under Arbitrary Network Delay in Wireless Sensor Networks

To cope with the arbitrariness of the network delays, a novel method, referred to as the composite particle filter approach based on variational Bayesian (VB-CPF), is proposed herein to estimate the clock skew and clock offset in wireless sensor networks. VB-CPF is an improvement of the Gaussian mixture kalman particle filter (GMKPF) algorithm. In GMKPF, Expectation-Maximization (EM) algorithm needs to determine the number of mixture components in advance, and it is easy to generate overfitting and underfitting. Variational Bayesian EM (VB-EM) algorithm is introduced in this paper to determine the number of mixture components adaptively according to the observations. Moreover, to solve the problem of data packet loss caused by unreliable links, we propose a robust time synchronization (RTS) method in this paper. RTS establishes an autoregressive model for clock skew, and calculates the clock parameters based on the established autoregressive model in case of packet loss. The final simulation results illustrate that VB-CPF yields much more accurate results relative to GMKPF when the network delays are modeled in terms of an asymmetric Gaussian distribution. Moreover, RTS shows good robustness to the continuous and random dropout of time messages.


Introduction
Wireless sensor networks (WSNs) consist of many low-cost sensor nodes capable of onboard sensing, computing and communications. WSNs are gaining importance since its wide applications, such as environment monitoring, object tracking and industrial machines controlling, etc. Most of these applications require the time of nodes to be synchronized to each other. Furthermore, some fundamental operations, such as power management, data fusion and transmission scheduling, etc. Require all the nodes running on a common timescale. However, in WSNs, every individual sensor works independently and maintains a local time measured by its own clock. This makes time synchronization between different nodes a critical piece of infrastructure [Wang, Jiang, Zhou et al. (2017) ;Qiu, Zhang, Qiao et al. (2018); ; Chen, Liu and Han (2018)]. Usually, time synchronization between any two nodes is accomplished through the exchange of time messages. Since deterministic and nondeterministic delay exists in the process of message transmission, the time messages may be arbitrarily delayed. At present, the common solution is to build a distribution model for the nondeterministic network delay, such as, Gaussian, Exponent, Gamma and Weber etc. [Wu, Chaudhari and Serpedin (2011); Wang, Jeske and Serpedin (2015); Noh, Chaudhari and Serpedin (2007)]. However, for sensor networks with complex real environment, various reasons will affect the distribution of network delay in varying degrees. It is difficult to find a network delay distribution model that is in line with the actual environment. Although the rationality and applicability of Gaussian network delay distribution model and Exponential network delay distribution model are verified in reference [Etzlinger, Wymeersch and Springer (2013); Abdel-Ghaffar (2002)], the simulation results in Noh et al. [Noh, Chaudhari and Serpedin (2007)] show that the estimation accuracy of clock parameters is very sensitive to the network delay distribution model. Therefore, it is necessary to research the clock parameter estimation method under arbitrary delay [Kim, Lee, Serpedin et al. (2009) ;Kim, Lee, Serpedin et al. (2011); Guo, Shen, Sun et al. (2015)]. Since Gaussian Mixture Model (GMM) can approximate arbitrary probability density [Anderson and Moore (1979)], Kim et al. [Kim, Lee, Serpedin et al. (2009)] estimated nondeterministic delay distribution using GMM, and proposed two estimation algorithms of clock offset, Gaussian Mixture Kalman Particle Filter (GMKPF), and Iterative Gaussian Mixture Kalman Particle Filter (IGMKPF), respectively. GMKPF combines measurement update steps based on Important Sampling (IS) with Gaussian Sum Filter (GSF) based on Kalman Filter (KF) for time update and proposed distribution generation. Then, Expectation Maximization (EM) algorithm is used to approximate the posterior distribution function of clock parameters by GMM. The introduction of EM algorithm not only alleviates the particle degradation problem caused by the particle filer algorithm, but also avoids the phenomenon that the number of GMM components increases exponentially with the iterations number increases. The simulation results of GMKPF show that when the network delay follows a single non-Gaussian (non-Exponential) distribution or a mixture of arbitrary distributions, it can maintain high synchronization accuracy with fewer message exchanges. However, GMKPF only tracks the clock offset and dose not estimate the clock skew, which will greatly reduce the synchronization period and increase communication overhead. Moreover, when GMKPF uses EM algorithm to approximate posterior distribution function with GMM, it needs to determine the components number of GMM beforehand, which is prone to under-fitting or over-fitting, and the estimation accuracy of parameters depends on the initial values setting. If set improperly, it is likely to converge to the local maximum. On the other hand, the unreliability of sensor network links may lead to the loss of time messages during transmission, while GMKPF does not discuss the algorithm performance in this case. Therefore, we propose a Composite Particle Filter Approach based on Variational Bayesian (VB-CPF) to realize the joint estimation of clock offset and skew in this paper. VB-CPF replaces the EM algorithm with Variational Bayesian EM (VB-EM) algorithm, which is used to estimate the GMM parameters of posterior distribution function. VB-EM algorithm is a deterministic approximate reasoning algorithm. It can determine the number of mixture components while determining the estimated values of mixture model parameters. Moreover, in order to solve the problem of time message loss caused by unreliable links, a robust time synchronization method (RTS) is proposed. RTS establishes an autoregressive model for clock skew, and utilizes the estimated clock skew obtained from each iteration of VB-CPF to estimate the autoregressive model parameters by recursive least squares method. If the time message is not received, the node can estimate the current clock parameters through the established clock skew autoregressive model, thus improving the robustness of the time synchronization method.

System model 2.1 Discrete clock model
The clock model of sensor node A is shown in Formula (1).
where A β and A θ represents the clock rate and initial clock phase, respectively. Since A β is time-varying, the above model can be expressed as an integral form, as shown in Time synchronization between nodes is usually achieved by the exchange of time stamps, which can be regarded as discrete samples of continuous time. Assuming that the sampling period is τ 0 , the discrete clock model of node A is 0 1 where ( ) A v n denotes the cumulative clock offset, ( ) A k β denotes the instantaneous clock skew at the k-th sampling. According to the analysis of Luo [Luo (2014)], the timevarying clock skew ( ) A n β can be described by the Gauss-Markov model, that is On the other hand, according to the definition of cumulative clock offset in formula (3), the recursive form of ( ) A v n can be written as Substituting formula (4) into formula (5), and obtain , according to formulas (4) and (6), the clock parameter evolution model of node A can be written in the following matrix form.

Local timestamp measurement model
In order to establish the clock relationship between two neighbor nodes, a two-way timestamps exchange mechanism is adopted in this paper. As shown in Fig. 1, in the th i round of message exchange, R sends a synchronization message to A at 1 i t , embedding its clock reading The above procedure can be modeled as where d represents the deterministic delay in message transmission between two nodes, i X and i Y represent the nondeterministic delays. Figure 1: Two-way time-stamps exchange between two nodes R and A For the time to complete a round of message exchange is very short, we assume that the clock parameters remain unchanged during a round of message exchange. Firstly, the clock model (2) is represented by reference clock and accumulative clock, as follows: To simplify symbolic representation ， let 1 1, , then formulas (8) and (9) can be simplified as Add formula (11) and formula (12), and let , then, a discrete local time measurement model is obtained by sampling, as shown in formula (13).
is the observation vector, n V is the observation noise. It is easy to observe that formulas (7) and (13) transform the estimation of clock parameters into Gauss-Markov estimation with unknown states.

KF based time update and proposed distribution generation
Firstly, it should be noted that all probability density function can be approximated by the GMM shown in formula (14) [Anderson and Moore (1979)].
where G is the number of mixture components in GMM, ( ) g ϕ represents the mixture weights and satisfies ( ) To simplify symbolic representation，let ( ) It can be known from formula (4), n u is a Gaussian noise with zero mean and variance 2 u σ , so the state noise in formula (7) obey At time n , the distribution of the current state can be predicted based on previous observations and state information, that is 1: 1 n n n n n n n G g g g n n n n n n n g G g g g n n n n n n n n g p z Theorem 1 [Anderson and Moore (1979)] If the state equation of the system is the form of can be expressed in form of Gauss sum as shown in formula (15), then the one-step predictive density function is uniformly close to Gauss sum According to Theorem 1, (18) can be written as follows: The mixture model parameters can be obtained from a set of parallel KF.
Then, according to the latest observation value n z obtained from the system observation Eq. (13), combined with the density functions n n n n n n n n n n n n n n where 1 1: 1 Theorem 2 [Anderson and Moore (1979)] If the observation equation of the system is the form of ( )  (19), then the updated posterior density function 1: Then update its weight . The posterior probability density 1: will be used as the proposed distribution function of the measurement update step, which is based on Importance Sampling (IS). It is easy to find that the number of mixture components in GMM of posterior probability density 1: increases from G to ' G , and with the increase of iterations, the number of mixture components increases exponentially, which will greatly increase the computational complexity. So, the number of mixture components must be reduced by adopting corresponding schemes.

Measurement update based on VB-EM algorithm
IS is a Monte Carlo method, which represents distribution ( ) p x by empirical approximation based on weighted particle (sample) set, that is, is a Dirac delta function and the weighted is obtained from the proposed distribution ( ) q x . The first step in implementing IS is to sample particles from the proposed distribution function 1: (22)), and then calculate their corresponding importance weights.
Then the posterior probability density 1: can be approximated as However, after many iterations, the weights of most particles are negligible, and only a few particles have large weights, thus resulting in the phenomenon of particle weight degradation. Although the introduction of resampling step alleviates the problem of particle degradation to a certain extent [Gordon, Salmond and Smith (1993)], excessive resampling will lead to particle depletion. Therefore, VB-EM algorithm is used to replace the resampling step to avoid particle depletion. Here, the GMM expression for posterior probability density 1: where G denotes the number of mixture components in GMMs, ( ) In order to facilitate the calculation of VB-EM algorithm, the precision matrix ( ) g n Λ is used to replace the covariance matrix ( ) g n Σ , in which the precision matrix is the inverse of the covariance matrix.
According to the above description, particle set  Since the Bayesian estimation is used to solve the GMM parameters, all parameters in the model are considered as random variables. Assume that the prior distribution of the mixture weight n ϕ obeys the Dirichlet distribution [Gorur and Rasmussen (2010)].
The prior distribution of mean n µ and precision matrix n Λ obeys the joint Gaussian-Wishart distribution [Bishop (2016)]. based on the a priori information described above, but it is difficult to directly calculate the posterior distribution with high dimension. Therefore, VB-EM introduces approximate distribution ( , , , ) VB-EM algorithm is a generalized EM algorithm, which is also an iterative algorithm. There are also two iterative steps, namely VB-E step and VB-M step.

VB-E:
Then find the partial derivative of ( ) L p in respect to ( ) lg t g t l g t T g t l g t n n n n n n n g t D g t g t n i l g t T g t l g t n n n n n can be obtained in two ways. The first one is obtained directly from the weighted sum of the particle set of formula (27) before the VB-EM algorithm is executed.
The second is calculated according to the formula (28) by using the weighted particle set fitted by GMM.
The estimation performance of formula (51) is better than that of formula (52). However, because of M G  , the first kind of computation complexity is higher than the second one.
The pseudo code of VB-CPF is as follows. Composite

Robust time synchronization method (RTS)
In practical sensor networks, the links between nodes are unreliable and susceptible to external interference, so the packets containing timestamps may be lost or collided during transmission. In order to solve the problem of data packet loss, the traditional wireless network can solve it by simply retransmitting. However, data packet retransmitting is infeasible for WSNs, because the energy cost of data retransmitting is too large, and the uncertainty delay in the process of data packet retransmitting will also affect the accuracy of time synchronization. So how to ensure the estimation accuracy of clock parameters in the case of data packet loss is the problem to be solved in this section.
According to Section 2.1, ( ) A n β varies with time in the actual environment, and for each sampling, it is not completely independent [Kim (2014); Kim, Ma and Hamilton (2012)]. Moreover, due to the insufficient energy of nodes and the change of temperature, it may change greatly. In this section, the time-varying clock skew is modeled as an Auto-Regressive (AR) process [Tibshirani (2011)], and then the parameters of AR model are estimated by recursive least square method based on the estimated clock skew obtained by VB-CPF. When encountering data packet loss, the clock parameters of the next time can be estimated according to the AR model, which ensures the robustness of VB-CPF method in unreliable link environment. The AR model of clock skew is as follows.
where P is the order of AR model, ( 1, 2,..., ) l l P π = is AR coefficient, ( ) n η is Gauss noise with zero mean and variance the time-varying characteristics of clock skew very well, and its computational complexity is not very high [Kim (2014)]. Next, the process of solving AR coefficients is described.
Assuming that the first five transmitted time messages are received correctly, according to the above estimation method, we can obtain five estimation values of clock skew, i.e., π π π β π β π β π β π β η β β β β η π π       = According to the recursive least squares estimation method, the mean square error matrix of 1 π is calculated first.
Then we can obtain the estimated value of 1 When 6 β is obtained, the AR coefficients are updated in the following recursive form according to the 1 M and 1 π calculated for the first time.
Step 1: calculate the mean square error matrix of 2 π , Step 2: update AR coefficient, 2  Figure 2: The overall method block diagram of RTS

Simulation experiment
To validate the performances of the proposed algorithm, simulation results are presented and compared to GMKPF. Since the VB-CPF method proposed in this paper is an improvement to GMKPF, in order to simplify the description of the simulation results, this section only evaluates the performance of the VB-CPF method under the asymmetric Gaussian delay model, in which the variance of the uplink nondeterministic Gaussian delay is 1 0.5 σ = , and the downlink is 2 1 σ = . And the other parameters used in simulations are as follows. Initial clock offsets, clock offsets and deterministic delays are uniformly drawn from ranges where S is the number of times the simulation experiment is executed, let 100 S = . The smaller the value of MSE, the higher the estimation accuracy of the clock parameters. Firstly, the effects of time message exchange times N (number of observations) and particle number M on the performance of VB-CPF are analyzed. As shown in Fig. 3 and Fig. 4, as the number of time message exchanges increases, the MSE of the clock parameter estimate gradually decreases, and when 20 N > , the magnitude of the MSE decrease become smaller. It can also be seen that when the number of particles is large ( 400 M ≥ ), the value of MSE is smaller. Obviously, through 5 times message exchanges, the estimated value of clock parameters obtained by sampling 500 particles per time is less than that obtained by 25 times message exchanges and sampling 100 particles per time. Therefore, the number of message exchanges can be reduced by increasing the number of particles sampled, thus reducing the communication overhead. As shown in Fig. 5 and Fig. 6, when the number of time message exchanges is fixed, the MSE of the estimated clock parameters decreases with the increase of the number of particles. When the number of particles is greater than 600, the MSE of the cumulative clock offset estimates decreases slightly. When the number of particles is greater than 500, the MSE of the clock skew estimates tends to be stable. Therefore, considering the synchronization accuracy and computational complexity, in the following comparative analysis of the performance of VB-CPF and GMKPF, we select the number of particles    Fig. 7 and Fig. 8 are the joint estimated value of cumulative clock offset and clock skew, which is obtained using the GMKPF algorithm. As can be seen from the figure, the MSE of the clock parameters estimated value decreases gradually with the increase of the number of time message exchanges. Comparing the two hybrid filtering methods: VB-CPF and GMKPF, the performance of VB-CPF is better than that of GMKPF. This is because VB-CPF uses VB-EM algorithm for Gaussian mixture fitting of sampled particles. VB-EM algorithm can adaptively determine the number of mixture components, avoiding the problem of over-fitting and under-fitting caused by the EM algorithm to determine the number of mixture components in advance.
Next, we evaluate the performance of VB-CPF and RTS in the case of time message lost. The simulation experiments consider two types of time message dropout: (1) continuous dropout, (2) random dropout. Fig. 9 shows the performance comparison of VB-CPF and RTS when time messages (observations) are lost five times in a row. Obviously, the performance of RTS is only 6 10 − difference from the estimation result without message dropout, which is better than that of VB-CPF. This is because, in each time of message loss, RTS estimates the cumulative clock offset and clock skew based on the preestablished clock skew model, while VB-CPF cannot update the estimated clock parameters. However, as the number of message loss increases, the MSE of the clock parameter estimate by RTS gradually increases. This is because the autoregressive model established for the clock skew is a statistical prediction model that predicts future values based on past estimated clock skew values. Therefore, when the time message is lost for a long time, the performance of the RTS will not be guaranteed. .39 10 − × on average. Obviously, the MSE of the clock parameter estimation obtained by RTS is smaller than that of VB-CPF, so the robustness of RTS in dealing with time message loss is better than that of VB-CPF. According to the analysis of Fig. 9 and Fig. 10, RTS can show better performance regardless of whether the time messages are continuous dropout or random dropout. This shows that RTS can well solve the problem of synchronization accuracy decline caused by the time message dropout.
(a) (b) Figure 10: When time messages are lost randomly, the MSE of (a) the cumulative clock offset estimated value and (b) the clock skew estimated value 6 Conclusions A composite particle filter approach based on variational Bayesian (VB-CPF) is proposed in this paper. VB-CPF replaces the EM algorithm used by GMKPF in estimating the parameters of Gaussian mixture model of the posterior distribution function of clock parameters with VB-EM algorithm, which makes it possible to determine the number of mixture components adaptively, thereby improving the estimation accuracy of clock parameters. At the same time, in order to solve the problem of time message dropout caused by unreliable link, a robust time synchronization method (RTS) is designed. RTS improves the robustness of VB-CPF in unstable network environment by establishing an autoregressive model for clock skew. The simulation results show that the estimation accuracy of clock parameters obtained by VB-CPF is better than that obtained by GMKPF, and RTS is robust to the continuous and random dropout of time messages.