Precise time synchronization clock parameter estimation method based on optimal invariant estimation

Abstract

The present invention relates to the technical field of data communication, and relates to a precise time synchronization clock parameter estimation method based on optimal invariant estimation. In the present invention, in a communication scene where unknown random time delay distribution information is present, on the basis of a precise time protocol, timestamp interaction is achieved by means of bidirectional data packet transfer between a master node and a slave node, and a clock relation model between the master node and the slave node is established; considering that a plurality of intermediate devices exist in the data packet interaction process, resulting that random time delay distribution information cannot be accurately obtained, a Gaussian mixture model is used to model a random time delay distribution, a competitive expectation maximization algorithm is used to learn mixture model parameters; and finally, on the basis of the parameters obtained by learning, an optimal invariant estimation method is used to estimate a clock frequency offset and a clock phase offset simultaneously, so as to realize time synchronization between the master node and the slave node. The present invention improves the accuracy of clock frequency offset estimation and clock phase offset estimation, is suitable for a communication scene where unknown random time delay distribution information is present, and enhances the robustness of the synchronization algorithm to a communication time delay.

Description

A precise clock parameter estimation method for time synchronization based on optimal invariant estimation Technical Field

[0001] The present invention belongs to the technical field of data communication and relates to a method for estimating precise time synchronization clock parameters based on optimal invariant estimation. Background Art

[0002] The essence of time synchronization lies in finding methods to estimate two clock parameters: the clock frequency offset and the clock phase offset. Only by simultaneously obtaining estimates of both can complete network time synchronization be achieved. If the synchronization method only estimates the frequency offset, while it can correct the frequency offset between nodes and maintain the same frequency variation between nodes, the initial phase offset between nodes still exists, preventing true clock synchronization. On the other hand, if the synchronization method only estimates the phase offset, nodes can correct the phase offset based on the estimated phase offset. However, since the frequency offset persists, the phase offset between nodes will reappear as the runtime increases. To ensure the quality of clock synchronization service, frequent resynchronization is required, which results in a significant waste of network resources. Therefore, to achieve complete network time synchronization while minimizing the additional energy consumption incurred by synchronization, it is necessary to simultaneously correct both the clock frequency offset and the clock phase offset.

[0003] The Precise Time Protocol (PTP) has become one of the most commonly used protocols in distributed networks due to its limited restrictions, high cost-effectiveness, and ability to ensure efficient distributed interaction. However, PTP is susceptible to packet delay variation, which can reduce synchronization accuracy and thus impact the quality of service. Therefore, finding methods to overcome the effects of packet delay variation is a key issue in ensuring the synchronization performance of PTP. Currently, there are two main solutions for overcoming the effects of packet delay variation in PTP-based synchronization methods: one uses Kalman filtering to address the clock parameter estimation problem for uncertain delay distributions, and the other uses a Gaussian mixture model to model random delay distributions and then uses the expectation-maximization algorithm to learn the model parameters. However, both solutions have certain limitations. The first solution only works for random delays that follow a single distribution and has relatively low robustness. The second solution's learned model parameters are easily affected by initial conditions and lack fault tolerance.

[0004] Therefore, there is an urgent need for a new method to estimate the clock frequency deviation and clock phase deviation parameters based on the precise time synchronization protocol, so that it can still maintain high synchronization accuracy and strong robustness in communication scenarios with unknown random delay distribution information.

[0005] Summary of the Invention

[0006] In view of this, the purpose of the present invention is to provide a Precise Time Protocol (PTP) clock parameter estimation method based on optimal invariant estimation. It is aimed at communication scenarios with unknown random delay distribution information. Based on the Precise Time Protocol (PTP), it overcomes the influence of data packet delay changes, establishes a clock relationship model between master and slave nodes, and simultaneously estimates two clock parameters, clock frequency deviation and phase deviation, to effectively ensure that the synchronization accuracy of the Precise Time Protocol is not affected. At the same time, it is applicable to different types of delay scenarios and has a strong robustness effect.

[0007] In order to achieve the above object, the present invention provides the following technical solutions:

[0008] A precise time synchronization clock parameter estimation method based on optimal invariant estimation is proposed. Timestamp interaction is achieved through bidirectional data packet transmission between the master node and the slave node. A clock relationship model between the master and slave nodes is established. A mixture Gaussian model is used to model the random delay distribution. The adversary-penalized expectation maximization algorithm is used to learn the mixture Gaussian model parameters. Finally, based on the learned parameters, the optimal invariant estimation method is used to simultaneously estimate the clock frequency deviation and phase deviation, thereby achieving time synchronization between the master and slave nodes.

[0009] The method specifically comprises the following steps:

[0010] S1: For communication scenarios with unknown random delay distribution information, a clock and delay relationship model between master and slave nodes is established based on the Precise Time Protocol (PTP).

[0011] S2: When the number of periodic synchronization rounds reaches the predetermined number N, a mixture Gaussian model is used to model the random delay distribution based on the timestamp sequence information recorded from node B. The adversary penalty expectation maximization algorithm is used to learn the model parameters of the mixture Gaussian model.

[0012] S3: Based on the learned Gaussian mixture model parameters and the timestamp sequence information recorded by slave B, the clock frequency and phase offsets of slave B are estimated using the optimal invariant estimation method.

[0013] S4: Based on the estimated values ​​of the clock frequency offset and phase offset, the slave node B corrects the local clock and achieves clock synchronization with the master node A.

[0014] Furthermore, in step S1, the specific steps of establishing the clock and delay relationship model between the master and slave nodes are:

[0015] S11: Master node A sends a synchronization message to slave node B at a fixed period T. The message does not contain a timestamp and records the sending time as t 1,i , and then send the timestamp t 1,i Encapsulated in a follow message and passed to slave node B;

[0016] S12: When the synchronization message sent by the master node A is successfully received, the slave node B records the receiving time as t 2,i At the same time, the slave node B obtains the sending timestamp t in the follow-up message 1,i ;

[0017] S13: Send a delay request message from node B to master node A, and record the message sending time as t 3,i ;

[0018] S14: When the delay request message sent by slave node B is successfully received, master node A records the receiving time as t 4,i Then the master node A delays the response message to send the timestamp t 4,i Send to slave node B;

[0019] S15: Obtain N sets of continuous timestamp sequences from node B through a series of message exchanges in the network According to the communication delay relationship, the clock information and delay relationship of the master and slave nodes can be modeled as follows:

[0020] in, and δ represent the relative frequency offset and relative phase offset of slave node B relative to master node A, respectively. ms and d sm Denote the deterministic delay on the forward and reverse paths, w 1,i and w 2,i denote the random delays for uplink and downlink, respectively.

[0021] Furthermore, in step S2, a mixture Gaussian model is used to model the random delay distribution and learn the model parameters. The specific steps are as follows:

[0022] S21: Assuming random delay w 1,1 ,…,w 1,N ,w 2,1 ,…,w 2,N They are independent of each other and obey the same Gaussian mixture distribution model as follows:

[0023] Where P(w|Θ) represents the Gaussian mixture distribution of random delays,

[0024] represents the weight coefficient of the unknown mixture component, K represents the number of mixture components, G(w|μ,σ 2 ) means the mean is μ and the variance is σ 2 The Gaussian probability density function of Represents a set of unknown parameters;

[0025] S22: According to Bayes’ theorem, any random delay belongs to the jth th The posterior density probability of the mixture component is:

[0026] Among them, h(j|w,Θ) represents the posterior density probability of the j-th mixture component, represents the probability density parameter of the j-th mixture component;

[0027] S23: Based on the Gaussian mixture distribution model of step S21 and the posterior density probability of step S22, a weighted likelihood framework function is constructed. The framework function needs to satisfy the condition that only Θ = Θ * When Θ * Representing the ideal true value of the parameter Θ, the weighted likelihood framework can reach the global maximum, and its formula is:

[0028] Where Q(Θ; W N ) represents the weighted likelihood framework, r(j|w i ,Θ) represents the evaluation function, w i represents the i-th random delay variable;

[0029] The weighted likelihood framework function consists of a generalized maximum expectation cost function and a measure for evaluating the uncertainty of random delay density; according to the hard-cut maximum expectation algorithm and the divergence function, r(j|w i ,Θ) is defined as: r(j|w i ,Θ)=2I(j|w i ,Θ)-h(j|w i ,Θ)

[0030] Among them, I(j|w i ,Θ) represents the indicator function, which is defined as:

[0031] S24: Use the weighted likelihood framework function as the cost function and maximize the cost function to estimate the unknown parameter set. In the process of iteratively updating the estimated value, add an opponent penalty mechanism to reward the winning mixed component and penalize the competitor mixed component. Continuously iterate the expectation calculation and maximization calculation until the algorithm converges to obtain the estimated value of the mixed model parameter set. The specific calculation formula is as follows:

[0032] Expected calculation:

[0033] Maximize the winner's mixture component:

[0034] Maximize the competitor mixture components:

[0035] in, c represents the winner component. If there is more than one winner component in the expected calculation result, the first one that appears is selected and marked as c, η and η β Represents the learning rate and satisfies η β =η.

[0036] Furthermore, in step S3, the clock frequency offset and phase offset of the slave node B are estimated using the optimal invariant estimation method. The specific steps are as follows:

[0037] S31: In most communication networks, uplink and downlink information will be communicated through the same path to maximize symmetry, so it is assumed that the delay d ms and d sm Equal and known, the timestamps obtained from node B are stacked into a matrix form, and the clock model can be expressed as:

[0038] The clock model is expressed in vector form as:

[0039] in,

[0040] S32: Definition represents the vector representation of clock frequency offset and clock phase offset, χ represents the finite parameter space of parameter θ, and its formula is Since the time interval between information exchanges is usually sufficient to ensure that the random delays are independent of each other, the joint probability density function of the observations is:

[0041] Among them, the coefficient Calculated from the Jacob matrix, d represents the fixed delay, which is assumed to be a known constant, and f1(·) and f2(·) are the probability density functions of the forward and backward paths, respectively;

[0042] S33: Definition The spatial position scale transformation group G is as follows:

[0043] in, g l,m (h) represents the position scale transformation formula, l represents the scale parameter, and m represents the position parameter. Since the probability density function of the observed data is unchanged under the position scale transformation group G, the corresponding scale and displacement parameters are used. The set of induced transformations can be defined as:

[0044] in, represents the induced transformation group, represents the induced transformation formula;

[0045] S34: The tilted normalized square error is used as the loss function of the optimal invariant estimation method. This loss function remains unchanged under the transformation group G. Its formula is:

[0046] in, express The estimated value of represents the estimated value of δ;

[0047] S35: Based on the Gaussian mixture model parameters learned in step S2, combined with the timestamp sequence information recorded from node B, the clock frequency deviation is obtained by Pittman estimation expansion. The optimal invariant estimate of the clock phase offset δ is given by:

[0048] in,

[0049] The beneficial effects of the present invention are as follows: Addressing the issue of the Precision Time Protocol (PTP) synchronization accuracy being easily affected by packet delay variations, the present invention proposes a method for overcoming packet delay variations in communication scenarios with unknown random delay distribution information, improving synchronization accuracy while also addressing the shortcomings of existing synchronization mechanisms. The present invention utilizes an adversary-penalized expectation-maximization algorithm, which is fault-tolerant to the pre-assigned density parameter setting and can adaptively reduce redundant components in density mixtures, while also improving the method's robustness to multiple types of delays.

[0050] Other advantages, objects, and features of the present invention will be described in part in the following description and, in part, will be apparent to those skilled in the art upon examination of the following or may be learned from practice of the present invention. The objects and other advantages of the present invention may be realized and obtained through the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention will be described in detail below with reference to the accompanying drawings, in which:

[0052] FIG1 is a schematic diagram of communication interaction between a master node A and a slave node B according to the present invention;

[0053] FIG2 is a flow chart of a method for estimating clock parameters of precise time synchronization based on optimal invariant estimation according to the present invention; DETAILED DESCRIPTION

[0054] The following describes the embodiments of the present invention by means of specific examples, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present invention, and the following embodiments and features in the embodiments can be combined with each other without conflict.

[0055] Please refer to Figures 1 and 2. Figure 1 is a schematic diagram of the communication interaction between the master node A and the slave node B of the present invention, where node A represents the master node and provides the reference time. Node B is the slave node and synchronizes its local clock with the master node clock by exchanging messages. The specific interaction process is as follows: master node A sends synchronization messages to slave node B at a fixed period T. During the i-th synchronization cycle, master node A sends synchronization messages to slave node B at a fixed period T. The message does not contain a timestamp and the sending time is recorded as t 1,i , and then send the timestamp t 1,i Encapsulated in the follow-up message and passed to the slave node B. When the synchronization message sent by the master node A is successfully received, the slave node B records the receiving time as t 2,i At the same time, the slave node B obtains the sending timestamp t in the follow-up message 1,i . Send a delay request message from node B to master node A, and record the message sending time as t 3,i When the delay request message sent by slave B is successfully received, master node A records the receiving time as t 4,i Then the master node A delays the response message to send the timestamp t4,i Sent to slave node B. Through the exchange of a series of messages in the network, slave node B obtains N sets of continuous timestamp sequences According to the communication delay relationship, the clock information and delay relationship of the master and slave nodes can be modeled as follows:

[0056] in, and δ represent the relative frequency offset and relative phase offset of slave node B relative to master node A, respectively. ms and d sm Denote the deterministic delay on the forward and reverse paths, w 1,i and w 2,i denote the random delays for uplink and downlink, respectively.

[0057] Assume random delay w 1,1 ,...,w 1,N ,w 2,1 ,…,w 2,N They are independent of each other and obey the same Gaussian mixture distribution model as follows:

[0058] in, represents the weight coefficient of the unknown mixture component, K represents the number of mixture components, G(w|μ,σ 2 ) means the mean is μ and the variance is σ 2 The Gaussian probability density function of . Let represents a set of unknown parameters. According to Bayes’ theorem, any random delay belongs to the jth th The posterior density probability of the mixture components is:

[0059] Based on the Gaussian mixture model and the posterior density probability, a weighted likelihood framework function is constructed, which needs to satisfy the condition that only Θ=Θ * When Θ * Representing the ideal true value of the parameter Θ, the weighted likelihood framework can reach the global maximum, and its formula is:

[0060] The weighted likelihood framework function consists of a generalized maximum expectation cost function and a measure for evaluating the uncertainty of random delay density. According to the hard-cut maximum expectation algorithm and the divergence function, r(j|w i ,Θ) is defined as: r(j|w i ,Θ)=2I(j|w i ,Θ)-h(j|w i ,Θ) (6)

[0061] Among them, I(j|wi ,Θ) represents the indicator function, which is defined as:

[0062] The weighted likelihood framework is used as the cost function, and the maximum cost function is adopted to estimate the unknown parameter set. In the process of iteratively updating the estimated value, an opponent penalty mechanism is added to reward the winning mixture component and penalize the competitor mixture component. The expectation calculation and maximization calculation are continuously iterated until the algorithm converges to obtain the estimated value of the mixture model parameter set. The specific calculation formula is as follows:

[0063] Expected calculation:

[0064] Maximization calculation (winning mixture component):

[0065] Maximization calculation (competitor mixture component):

[0066] in, c represents the winner component. If there is more than one winner component in the expected calculation result, the first one that appears is selected and marked as c, η and η β Represents the learning rate and satisfies η β =η. By continuously iterating until the algorithm converges, the model parameters of the random delay distribution mixture model can be obtained.

[0067] Assume that the delay d is determined ms and d sm Equal and known, the timestamps obtained from node B (Formula (1) and Formula (2)) are stacked into a matrix form, and the clock model can be expressed as:

[0068] The clock model is expressed in vector form as:

[0069] in, definition represents the vector representation of clock frequency offset and clock phase offset, χ represents the finite parameter space of parameter θ, and its formula is Since the time interval between information exchanges is usually sufficient to ensure that the random delays are independent of each other, the joint probability density function of the observations is:

[0070] Among them, the coefficient Calculated from the Jacob matrix. Definition The spatial position scale transformation group G is as follows:

[0071] in, Since the probability density function of the observed data is unchanged under the position scale transformation group G, the corresponding scale and displacement parameters The set of induced transformations can be defined as:

[0072] The tilted normalized square error is used as the loss function of the optimal invariant estimation method. The loss function remains unchanged under the transformation group G, and its formula is:

[0073] Based on the learned Gaussian mixture model parameters, combined with the timestamp sequence information recorded from node B, the clock frequency deviation is obtained by Pittman estimation expansion. The optimal invariant estimate of the clock phase offset δ is given by:

[0074] in,

[0075] FIG2 is a flow chart of the precise time synchronization method based on optimal invariant estimation of the present invention. The present invention provides a precise time synchronization method based on optimal invariant estimation, as shown in FIG2 , and the specific steps are as follows:

[0076] M1: The master-slave node interaction process begins;

[0077] M2: Master node A sends a synchronization message without a timestamp to slave node B, and records the local time of the message as t 1,i ;

[0078] M3: Master node A will send timestamp t 1,i Encapsulate it in a follow message and send it to slave node B;

[0079] M4: After receiving the synchronization message from node B, the local time of message arrival is recorded as t 2,i , then sends a delay request message to the master node A and records the local time of sending as t 3,i ;

[0080] M5: Master node A records the arrival time of the delay request message as t 4,i , and then the timestamp t 4,i Encapsulated in a delay response message and delivered to slave node B;

[0081] M6-M7: Determine whether the synchronization round has reached the set value N. If so, the slave B estimates the clock frequency offset and phase offset based on the obtained timestamp sequence. Otherwise, the synchronization round i = i + 1, and enter M2 to continue timestamp exchange;

[0082] M8: Node B uses a Gaussian mixture model to model the random delay distribution based on the obtained timestamp sequences. It then iterates the adversary penalty expectation maximization algorithm to learn the mixture model parameters until the algorithm converges.

[0083] M9: Based on the timestamp information and the model parameters learned in the previous step, the slave node B jointly estimates the clock frequency offset and phase offset.

[0084] M10: Based on the estimated clock frequency and phase offsets, slave node B corrects its local clock to achieve clock synchronization with master node A.

[0085] M11: The master-slave node interaction process ends.

[0086] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solutions, which should all be included in the scope of the claims of the present invention.

Claims

1. A method for estimating clock parameters for precise time synchronization based on optimal invariant estimation, characterized in that: The method specifically comprises the following steps: S1: For communication scenarios with unknown random delay distribution information, a clock and delay relationship model between master and slave nodes is established based on the precise time synchronization protocol; S2: When the number of periodic synchronization rounds reaches the predetermined number of synchronization rounds N, a mixed Gaussian model is used to model the random delay distribution according to the timestamp sequence information recorded from node B, and the adversary penalty expectation maximization algorithm is used to learn the model parameters of the mixed Gaussian model; S3: Based on the learned Gaussian mixture model parameters and the timestamp sequence information recorded from node B, the clock frequency deviation and phase deviation of node B are estimated using the optimal invariant estimation method; S4: According to the estimated values ​​of the clock frequency deviation and phase deviation, the local clock of the slave node B is corrected to achieve clock synchronization with the master node A.

2. The method for estimating precise time synchronization clock parameters according to claim 1, characterized in that: In step S1, the specific steps of establishing the clock and delay relationship model between the master and slave nodes are: S11: Master node A sends a synchronization message to slave node B at a fixed period T. The message does not contain a timestamp and records the sending time as t 1,i , and then send the timestamp t 1,i Encapsulated in a follow message and transmitted to slave node B; S12: When the synchronization message sent by the master node A is successfully received, the slave node B records the receiving time as t 2,i At the same time, the slave node B obtains the sending timestamp t in the follow-up message 1,i ; S13: Send a delay request message from node B to master node A, and record the message sending time as t 3,i ; S14: When the delay request message sent by slave node B is successfully received, master node A records the receiving time as t 4,i Then the master node A delays the response message to send the timestamp t 4,i Send to slave node B; S15: Obtain N groups of continuous timestamp sequences from node B through a series of message exchanges in the network According to the communication delay relationship, the clock information and delay relationship of the master and slave nodes are modeled as follows: in, and δ represent the relative frequency offset and relative phase offset of slave node B relative to master node A, respectively. ms and d sm denote the deterministic delay on the forward and reverse paths, respectively, and w 1,i and w 2,i denote the random delays for uplink and downlink, respectively.

3. The method for estimating precise time synchronization clock parameters according to claim 2, characterized in that: In step S2, a mixed Gaussian model is used to model the random delay distribution and learn the model parameters. The specific steps are: S21: Assuming random delay w 1,1 ,...,w 1,N ,w 2,1 ,...,w 2,N They are independent of each other and obey the same Gaussian mixture distribution model as follows: Where P(w|Θ) represents the Gaussian mixture distribution of random delays, represents the weight coefficient of the unknown mixture component, K represents the number of mixture components, G(w|μ,σ 2 ) means the mean is μ and the variance is σ 2 The Gaussian probability density function of Represents a set of unknown parameters; S22: According to Bayes’ theorem, any random delay belongs to the jth th The posterior density probability of the mixture component is: Among them, h(j|w,Θ) represents the posterior density probability of the jth mixture component, represents the probability density parameter of the jth mixture component; S23: Based on the Gaussian mixture distribution model of step S21 and the posterior density probability of step S22, a weighted likelihood framework function is constructed, which needs to satisfy the condition that only Θ=Θ * When * Representing the ideal true value of parameter Θ, the weighted likelihood framework reaches the global maximum, and its formula is: Among them, Q(θ; W N ) represents the weighted likelihood framework, r(j|w i ,Θ) represents the evaluation function, w i represents the i-th random delay variable; The weighted likelihood framework function consists of a generalized maximum expectation cost function and a measure for evaluating the uncertainty of random delay density; according to the hard-cut maximum expectation algorithm and the divergence function, r(j|w i ,Θ) is defined as: r(j|w i ,Θ)=2I(j|in i ,Θ)-h(j|w i ,Θ) Among them, I(j|w i ,Θ) represents the indicator function, which is defined as: S24: Use the weighted likelihood framework function as the cost function, and use the maximum cost function to estimate the unknown parameter set. In the process of iteratively updating the estimated value, add the opponent penalty mechanism, reward the winning mixed component, and punish the competitor's mixed component. Continuously iterate the expected calculation and maximization calculation until the algorithm converges to obtain the estimated value of the mixed model parameter set. The specific calculation formula is as follows: Expected calculation: Maximize the winner mixture component: Maximize the competitor mixture components: in, c represents the winner component. If there is more than one winner component in the expected calculation result, the first one that appears is selected and marked as c, η and η β Represents the learning rate and satisfies η β =η.

4. The method for estimating precise time synchronization clock parameters according to claim 3, characterized in that: In step S3, the clock frequency deviation and phase deviation of the slave node B are estimated using the optimal invariant estimation method. The specific steps are as follows: S31: Assume that the delay d is determined ms and d sm Equal and known, the timestamps obtained from node B are stacked into a matrix form, and the clock model is expressed as: The clock model is expressed in vector form as: in, S32: Definition represents the vector representation of clock frequency offset and clock phase offset, χ represents the finite parameter space of parameter θ, and its formula is The joint probability density function of the observed data is: Among them, the coefficient Calculated from the Jacob matrix, d represents the fixed delay, which is assumed to be a known constant, and f1(·) and f2(·) are the probability density functions of the forward and backward paths, respectively; S33: Definition The spatial position scale transformation group G is as follows: in, g l,m (h) represents the position scale transformation formula, l represents the scale parameter, and m represents the position parameter; since the probability density function of the observed data is unchanged under the position scale transformation group G, according to the corresponding scale and displacement parameters The induced transformation group is defined as: in, represents the induced transformation group, represents the induced transformation formula; S34: The tilt normalized square error is used as the loss function of the optimal invariant estimation method. The loss function satisfies the requirement of being invariant under the transformation group G, and its formula is: in, express The estimated value of represents the estimated value of δ; S35: Based on the Gaussian mixture model parameters learned in step S2, combined with the timestamp sequence information recorded from node B, the clock frequency deviation is obtained by Pittman estimation expansion. The optimal invariant estimate of the clock phase offset δ is given by: in,

Images