Clock synchronization method for adaptive Kalman filtering based on state extension, computer program product and electronic equipment

Abstract

The invention discloses an adaptive Kalman filtering clock synchronization method based on state extension, a computer program product and electronic equipment, and belongs to the technical field of computer software, and the method comprises the steps: building a discrete time clock synchronization state model based on a precision time protocol PTP; establishing a state space model for estimating clock skew and clock drift; and carrying out clock state estimation based on adaptive Kalman filtering of state expansion. According to the method, the asymmetry of delay is described by adopting a discrete linear difference equation, and input lambda k is added in a state vector to obtain an extended state vector, so that an extended state space model is constructed, the influence of the asymmetry of delay on clock skew is reflected more truly, and on the basis, the clock skew is estimated by using a self-adaptive Kalman filtering algorithm, so that the accuracy of clock skew estimation is improved. According to the method, the optimal clock skew estimation processing can be provided while the low algorithm complexity is kept, so that the accuracy of clock parameter estimation in a real-time system is further improved.

Description

Technical Field

[0001] The present invention relates to the technical field of computer software, and in particular to a clock synchronization method, a computer program product, and an electronic device based on state-expanded adaptive Kalman filtering. Background Art

[0002] Clock synchronization plays a crucial role in modern computing and communication systems, especially for real-time systems with strict time constraints. The Precision Time Protocol (PTP) under the IEEE 1588 standard is widely adopted due to its efficient implementation, low resource occupancy, and high-precision time synchronization characteristics. The reliability of PTP is reflected in two aspects: synchronization accuracy and fault tolerance. In the PTP protocol, the slave device calculates the clock offset with the master device by exchanging packets with timestamps, where the clock offset refers to the time difference between different clocks. The effectiveness of this process is based on a key assumption that the network path delay is symmetric. In addition, the synchronization accuracy provided by PTP is directly affected by the accuracy of event timestamps. However, in the actual environment, the end-to-end path delay is often not symmetric, and the timestamps may be affected by factors such as retransmission due to packet loss, local oscillator failure, etc., resulting in uncertainty, which will cause significant time offset and time drift, where time drift refers to the cumulative difference between clocks. Therefore, in order to accurately estimate the clock parameters, the design of PTP must take into account the delay asymmetry and the uncertainty of timestamps. Some models using probabilistic methods, such as the Kalman filter, are suitable for estimating the change of clock offset. In particular, the Adaptive Kalman Filter (AKF) can provide excellent clock offset estimation processing while maintaining a low algorithm complexity, but the AFK still does not fully consider the delay asymmetry, so there is an urgent need to further improve the accuracy of clock parameter estimation. Summary of the Invention

[0003] The object of the present invention is to overcome the problems of the prior art, and provide a clock synchronization method, a computer program product, and an electronic device based on state-expanded adaptive Kalman filtering.

[0004] The object of the present invention is achieved by the following technical solutions: A clock synchronization method based on state-expanded adaptive Kalman filtering, the method comprising the following steps:

[0005] Establish a discrete-time clock synchronization state model based on the Precision Time Protocol (PTP), including: modeling the clock offset estimation and clock drift estimation between the master and slave clocks during the synchronization process according to the PTP message exchange process; describing the change of delay in the clock offset estimation through a discrete linear difference equation to obtain the discrete linear difference equation of delay asymmetry;

[0006] Establish a state space model for estimating clock offset and clock drift, including: rewriting the measurement equations for clock offset estimation and clock drift estimation according to the clock offset estimation and clock drift estimation models of the master and slave clocks, and giving the representation of the measurement noise random variables related to clock offset and clock drift; deriving the recursive state equations for clock offset and clock drift based on stochastic differential equations; deriving the vector-matrix representation of the state space model of PTP according to the measurement equations and recursive state equations;

[0007] Perform clock state estimation based on state-expanded adaptive Kalman filtering, including: adding the input delay asymmetry to the state vector to obtain an expanded state vector; using the PTP timestamp measurement to numerically describe the recurrence of delay asymmetry according to the discrete linear difference equation of delay asymmetry to obtain a delay asymmetry supplementary equation; establishing the state equation and measurement equation of the expanded state space model; deriving the adaptive Kalman filtering algorithm for estimating the expanded state vector according to the expanded state space model, and using the adaptive Kalman filtering algorithm to estimate the clock offset during the PTP clock synchronization process.

[0008] In one example, the discrete linear difference equation of the delay asymmetry is:

[0009] λ k =λ k-1 +[1, 1]α k +[-1, -1]β k ;

[0010] α k =[Δt k , Δt′ k ;

[0011] β k =[ΔC(t k +d ms,k ), ΔC(t′ k -d sm,k )];

[0012] where λ k is the delay asymmetry at the k-th time point; λ k-1 is the delay asymmetry at the (k - 1)-th time point; Δt kDenote the change in the timestamp of the synchronization message sent by the master clock between the k-th time point and the (k - 1)-th time point; Δt′ k Denote the change in the received time difference recorded after receiving the delay request message from the master clock between the k-th time point and the (k - 1)-th time point; ΔC(t k +d ms,k ) Denote the change in the received time difference recorded after receiving the synchronization message from the slave clock at the k-th time point and delaying it by time d ms,k ; ΔC(t′ k -d sm,k ) Denote the change in the timestamp of the response to the delay request message after recording the reception time C(t k +d ms,k ) at the k-th time point and the (k - 1)-th time point from the slave clock.

[0013] In one example, after the steps of the measurement equations for the rewritten clock offset estimation and clock drift estimation, it further includes:

[0014] Introduce a Bernoulli random process to optimize the clock offset measurement equation, and the optimized clock offset measurement equation is:

[0015]

[0016] Where, Denote the estimated value of the clock offset at the k-th time point, θ k Denote the actual clock offset at the k-th time point; λ k Is the delay asymmetry at the k-th time point; n B，k Denote the Bernoulli random process; v θ，k Is a random variable representing the measurement noise related to the clock offset.

[0017] In one example, the recursive state equations for the clock offset and clock drift are:

[0018]

[0019] Where, θ k Denote the actual clock offset at the k-th time point, θ k-1 Denote the actual clock offset at the (k - 1)-th time point; Δτ represents the time point; γ k Denote the actual clock drift at the k-th time point; γ k Denote the actual clock drift at the (k - 1)-th time point; ω θ，k-1 , ω γ，k-1 Are two uncorrelated Gaussian random noise processes with zero mean.

[0020] In one example, the vector-matrix representation of the state space model of the PTP is:

[0021]

[0022] w k-1 = [ω θ，k-1 , ω γ，k-1 T ;

[0023] v k = [v θ，k , ν γ,k T ;

[0024] x k represents the clock state vector at the k-th time point; x k-1 represents the clock state vector at the (k - 1)-th time point; w k , v k are defined as Gaussian white noises with zero mean and uncorrelated with each other, w k-1 represents the Gaussian white noise at the (k - 1)-th time point; z k represents the measurement vector; λ k is the delay asymmetry at the k-th time point k; n B，k represents a Bernoulli random process; ω θ，k-1 , ω γ，k-1 are two uncorrelated Gaussian random noise processes with zero mean; T represents vector transpose; v θ，k , v γ，k are Gaussian random noise processes with zero mean and finite variance; A and C are both state transition matrices; D, E, H, and G are all measurement matrices.

[0025] In an example, the state equation and measurement equation of the extended state space model are:

[0026]

[0027] where, represents the extended state vector at the k-th time point; represents the extended state vector at the (k - 1)-th time point; A a , C a both represent extended state transition matrices; represents the extended state Gaussian white noise at the (k - 1)-th time point; z k represents the extended state measurement vector at the k-th time point; D a represents the extended state measurement matrix; H, G are both measurement matrices; n B，k represents a Bernoulli random process; v k , w k-1 , ω λ，k-1 ​​They are all Gaussian white noises with zero mean and uncorrelated with each other; T represents vector transpose.

[0028] In one example, the adaptive Kalman filtering algorithm for estimating the extended state vector derived according to the extended state space model is used to estimate the clock offset in the PTP clock synchronization process, including:

[0029] Initialization processing: Define the prior estimation covariance matrix of the extended state vector and the posterior estimation covariance matrix of the extended state vector;

[0030] Prediction processing: Perform prior state estimation according to the state equation and update the prior estimation error covariance matrix;

[0031] Adaptive adjustment: Adjust the measurement noise of the extended state space model according to the weighting coefficient, and correct the difference between the measured value and the predicted value based on the prior state estimation according to the residual;

[0032] Update processing: Calculate the Kalman gain, and obtain the posterior estimation and the posterior estimation error covariance matrix by updating the prior estimation and the measured value;

[0033] Clock offset estimation processing: Compensate the slave clock according to the sum of the clock offset, clock drift and delay asymmetry in the posterior estimation.

[0034] In one example, after the step of using the adaptive Kalman filtering algorithm to estimate the clock offset in the PTP clock synchronization process, it further includes detecting the time stamp by using residual awareness, including the following sub-steps:

[0035] Calculate the squared Mahalanobis distance of the residual, compare the squared Mahalanobis distance with the threshold related to the expected confidence level. If the squared Mahalanobis distance is greater than the threshold related to the expected confidence level, mark the time stamp as an outlier, otherwise, the time stamp passes the verification and is used for posterior estimation.

[0036] It should be further noted that the technical features corresponding to the above examples can be combined or replaced with each other to form a new technical solution.

[0037] The present invention also includes a computer program product, including a computer program, which when executed by a processor implements the steps of the above-mentioned clock synchronization method based on state expansion adaptive Kalman filtering formed by any one example or a combination of multiple examples.

[0038] The present invention further includes an electronic device, comprising a memory and a processor, wherein computer instructions that can run on the processor are stored on the memory, and when the processor runs the computer instructions, the steps of the clock synchronization method based on state-expanded adaptive Kalman filtering formed by any one or more of the above examples are executed.

[0039] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0040] 1. In one example, a discrete linear difference equation is used to describe the delay asymmetry, and the input λ is added to the state vector k to obtain an extended state vector, and then an extended state space model is constructed, which more realistically reflects the influence of delay symmetry on clock offset. On this basis, the adaptive Kalman filtering algorithm is used to estimate the clock offset, which can provide optimal clock offset estimation processing while maintaining a low algorithm complexity, thereby further improving the accuracy of clock parameter estimation in a real-time system.

[0041] 2. In one example, compared with using white noise with a finite variance to represent the uncertainty of timestamps in PTP measurements, the present invention gives a better PTP measurement equation through a Bernoulli random process, and then derives a state space model that more conforms to the actual situation, thereby improving the prediction performance of the model.

[0042] 3. In one example, by introducing residual awareness, detection conditions for timestamp verification are proposed, which can effectively detect and preprocess unreliable timestamps, improving the robustness and clock synchronization accuracy of the entire clock synchronization method; using the verified timestamps for posterior estimation can avoid estimation errors caused by abnormal data, further improving the accuracy of clock synchronization. Description of the Drawings

[0043] The following further details the specific embodiments of the present invention in conjunction with the drawings. The drawings described herein are used to provide a further understanding of the present application and form a part of the present application. The same reference numerals are used to represent the same or similar parts in these drawings. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation of the present application.

[0044] Figure 1 It is a flowchart of the method provided by one example of the present invention;

[0045] Figure 2 It is a schematic diagram of the PTP protocol message exchange process based on the IEEE 1588 standard provided by one example of the present invention. Detailed implementation mode

[0046] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0047] In the description of the present invention, it should be noted that the directions or positional relationships indicated by terms such as "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", etc. are based on the directions or positional relationships shown in the accompanying drawings. It is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus cannot be construed as a limitation to the present invention. In addition, the use of ordinal numbers (for example, "first and second", "first to fourth", etc.) is for distinguishing objects and is not limited to this order, and cannot be construed as indicating or implying relative importance.

[0048] In the description of the present invention, it should be noted that unless otherwise clearly specified and limited, the terms "installation", "connection", and "connection" should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be directly connected or indirectly connected through an intermediate medium, and it can be the communication inside two elements. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific situations.

[0049] In addition, the technical features involved in different embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0050] In one example, as Figure 1 shown, a clock synchronization method based on state-expanded adaptive Kalman filtering can be applied to clock synchronization systems in automotive driving, aerospace, industrial automation, etc. In the clock synchronization system of automotive driving, the method includes the following steps:

[0051] S1: Establish a discrete-time clock synchronization state model based on the Precision Time Protocol (PTP), including the following sub-steps:

[0052] S11: Model the clock offset estimation and clock drift estimation between the master and slave clocks during the synchronization process according to the PTP message exchange process;

[0053] S12: Describe the change of delay in the clock offset estimation through a discrete linear difference equation to obtain the discrete linear difference equation of delay asymmetry.

[0054] Step S1 is the basic model of the whole method. According to the clock synchronization message exchange process of the PTP protocol, it describes and models the clock offset and drift of PTP clock synchronization, uses discrete linear difference equations to describe the change of delay, linearizes the random relationship between sequential delay asymmetries, and abstracts the clock offset problem to deeply study clock offset estimation, delay asymmetry analysis and related issues, providing the necessary basic equations for subsequent analysis.

[0055] S2: Establish a state space model for estimating clock offset and clock drift, including the following sub-steps:

[0056] S21: Rewrite the measurement equations for clock offset estimation and clock drift estimation according to the clock offset estimation and clock drift estimation models of the master and slave clocks, and give the representation of the measurement noise random variables related to clock offset and clock drift;

[0057] S22: Derive the recursive state equations for clock offset and clock drift based on stochastic differential equations;

[0058] S23: Derive the vector-matrix representation of the state space model of PTP according to the measurement equations and the recursive state equations.

[0059] The state space model for estimating clock offset and drift in step S2 is implemented based on step S1, providing the state equations and measurement equations that can be used for Kalman filtering. That is, the process of estimating the clock offset of the slave clock by the state space model can be modeled as a Kalman filtering problem, which includes specific state equations and measurement equations, respectively defining the state transition between different time points of the clock offset and the relationship between the observation and the state, and can effectively predict the next state of the clock offset and correct the predicted state in Kalman filtering.

[0060] S3: Perform clock state estimation based on adaptive Kalman filtering with state extension, including the following sub-steps:

[0061] S31: Add the input delay asymmetry to the state vector to obtain the extended state vector.

[0062] S32: According to the discrete linear difference equation of the delay asymmetry, use the PTP timestamp measurement to numerically describe the recurrence of the delay asymmetry to obtain the delay asymmetry supplementary equation.

[0063] S33: Establish the state equations and measurement equations of the extended state space model, where the parameters can be assumed to be Gaussian white noise with zero mean and uncorrelated with each other, and give the noise and its related noise covariance matrix.

[0064] S34: Derive an adaptive Kalman filtering algorithm for estimating the extended state vector according to the extended state space model, and use the adaptive Kalman filtering algorithm to estimate the clock offset during the PTP clock synchronization process. Among them, the adaptive Kalman filtering algorithm can estimate the clock offset during the PTP clock synchronization process according to the general Kalman filtering process.

[0065] The state space model of PTP is a set of linear state and measurement equations. Using the usual Kalman filter is not sufficient to estimate the clock offset caused by the delay asymmetry λ k Even if λ can be described by PTP timestamp measurements, accurately modeling its dynamics is still challenging due to unreliable time measurements. k Therefore, step S3 considers the delay asymmetry and uses an adaptive Kalman filtering algorithm for state estimation to optimize the clock synchronization accuracy. Specifically, in order to handle the delay asymmetry during the clock synchronization process, step S3 of the present invention extends the state vector by using the delay asymmetry analysis based on the state space model in step S2, and derives the state equation and measurement equation of the extended state space model. On this basis, an adaptive Kalman filtering algorithm for estimating the extended state vector is given, and then the clock state estimation is carried out to improve the clock synchronization accuracy.

[0066] The model in step S1 of the present invention is the basis of the whole method, providing the necessary basic equations and models; the model in step S2 is analyzed based on the model in step S1 to obtain the state and measurement equations that can be used for Kalman filtering; the clock state estimation method in step S3 gives an adaptive Kalman filtering method for state extension based on steps S1 and S2, and uses the optimal estimation of the clock offset.

[0067] In an example, model the clock offset estimation and clock drift estimation between the master and slave clocks during the synchronization process according to the PTP message exchange process, including the following processing procedures:

[0068] The Precision Time Protocol (PTP) is an accurate time synchronization protocol based on the IEEE 1588 standard. It calculates the clock offset between the master and slave clocks through regular message exchanges between the master and slave clocks, and compensates the slave clock to achieve the purpose of clock synchronization. The PTP protocol message exchange process based on the IEEE 1588 standard is as Figure 2 shown. In the kth synchronization iteration, the master clock starts the two-way message exchange by sending a synchronization (sync) message containing the local time t k . After sending the sync message, the master clock then sends a supplementary follow-up message, which contains the timestamp t k . Once the slave clock receives the sync message, after a delay d ms,k , it records the reception time C(tk +d ms,k )。The slave clock then responds with a delay - req message, which contains C(t′ k -d sm,k )。After receiving the delay - req message, the master clock records the reception time t′ k , and sends it back to the slave clock via a delay - resp message. Finally, after receiving the delay - resp message, all the required timestamps {t k , C(t k +d ms,k ), C(t′ k -d sm,k ), t′ k} of the slave clock are collected, and the four timestamps are used for clock parameter estimation. By using these four timestamps, the clock offset estimation and clock drift estimation between the master and slave clocks during the synchronization process can be modeled as:

[0069]

[0070] where, represents the estimated value of the clock offset at the k - th time point, represents the estimated value of the clock offset at the (k - 1)-th time point; represents the estimated value of the clock drift at the current time step k; t k represents the timestamp when the master clock sends the synchronization message, C(t k +d ms,k ) represents the reception time recorded by the slave clock after receiving the synchronization message and delaying for a time d ms，k , C(t′ k -d sm,k ) represents the timestamp when the slave clock responds to the delay - req message after recording the reception time C(t k +d ms,k ), and t′ k represents the reception time recorded by the master clock after receiving the delay - req message.

[0071] Assume that the actual clock offset in the k - th synchronization iteration is θ k , and further assume that it remains unchanged within a sufficiently short time interval. At this time, there is:

[0072] C(t k +d) = t k +d+θ k ;

[0073] where, d represents a short time interval. Combining the above formulas, we get:

[0074]

[0075] where λ k represents the delay asymmetry in the k-th synchronization iteration. It is assumed here that the maximum allowable values of d ms，k and d sm，k are finite and, in general, three times the synchronization interval according to the IEEE 1588 PTP standard. Obviously, this asymmetry causes a deviation in the convergence point of the estimated clock offset. Therefore, the synchronization accuracy is very sensitive to the delay asymmetry.

[0076] Since the delay asymmetry is caused by different factors of timestamp implementation, physical layer hardware, and delay components, it is not feasible to directly handle the delay asymmetry. To quantify this asymmetry, the present invention describes the change in delay through a discrete linear difference equation (LDE). Specifically, the delay asymmetry Δλ k at the k-th time step is defined as a function of the delay difference between the current time point and the previous time point, i.e., Δλ k = λ k - λ k-1 , and can be obtained from PTP timestamp measurements. In one example, the present invention constructs the discrete linear difference equations of d ms，k and d sm，k by defining the change in delay asymmetry of PTP measurements through the derivative of the timestamp with respect to the time step:

[0077] Δd ms，k = d ms，k - d ms，k-1 = ΔC(t k + d ms,k ) - Δt k ;

[0078] Δd sm,k = d sm,k - d sm,k-1 = ΔC(t′ k - d sm,k ) - Δt′ k ;

[0079] where Δ represents numerical differentiation. Further derivation gives the discrete linear difference equation LDE of the delay asymmetry as:

[0080] λ k = λ k-1 + [1, 1]α k + [-1, -1]β k ;

[0081] where α k = [Δt k , Δt′ k , βk = [ΔC(t k + d ms,k ), ΔC(t′ k - d sm,k )]. λ k is the delay asymmetry at the k-th time point; λ k-1 is the delay asymmetry at the (k - 1)-th time point; Δt k represents the change in the timestamp of the synchronization message sent by the master clock between the k-th time point and the (k - 1)-th time point; Δt′ k represents the change in the received time difference recorded after receiving the delay request message from the master clock between the k-th time point and the (k - 1)-th time point; ΔC(t k + d ms,k ) represents the change in the received time difference recorded after receiving the synchronization message from the slave clock at the k-th time point and the (k - 1)-th time point and delaying it by the time d ms，k ; ΔC(t′ k - d sm,k ) represents the change in the timestamp of the response to the delay request message after recording the received time C(t k + d ms,k ) from the slave clock at the k-th time point and the (k - 1)-th time point. This example uses a discrete linear difference equation to describe the change in delay, and the random relationship between sequential delay asymmetries is linearized through the given timestamp measurements.

[0082] In one example, according to the clock offset estimation and clock drift estimation models of the master and slave clocks, the measurement equations for clock offset estimation and clock drift estimation are rewritten, including:

[0083] The process of the slave clock estimating the clock offset can be modeled as a Kalman filtering problem, where state and measurement equations need to be designed to construct a state space model. According to the offset and drift models of the master and slave clocks, the measurement equations of PTP can be rewritten as:

[0084]

[0085] where θ k represents the actual clock offset at the k-th time point; γ k represents the actual clock drift at the current time step k; v θ，k and vvγ,k are random variables, representing the measurement noise related to clock offset and clock drift respectively. Assume that νθ,k and vvγ,k follow a Gaussian process with a mean of zero and have finite variances and Considering the uncertainties of the master clock and the slave clock, the variance of the measurement noise related to the clock offset can be calculated as where, and respectively represent the timestamp uncertainties of the master clock and the slave clock.

[0086] The noise in PTP measurements mainly consists of timestamp uncertainties that originate from multiple different factors such as finite clock resolution, clock phase noise, delay timestamps, etc. In an ideal situation, the timestamp inaccuracies can be well approximated by white noise with a finite variance. However, in practice, this is usually not the case, and a better PTP measurement equation is obtained by introducing a Bernoulli random process n B，k is given. In one example, after step S21, it further includes:

[0087] Introduce a Bernoulli random process to optimize the clock offset measurement equation, and the optimized clock offset measurement equation is:

[0088]

[0089] where the additional term n B，k is a random value with a given probability density function, generated with probability p. By introducing a Bernoulli random process to represent the timestamp uncertainty and optimizing the measurement equation, a vector-matrix form representation of the PTP state space model is derived, providing a basic state space model for step S3.

[0090] In one example, in step S22, in the design of the state equation, the recursive state equations for clock offset and clock drift are:

[0091]

[0092] where, θ k-1 represents the actual clock offset at the (k - 1)-th time point; Δτ represents the time point; γ k represents the actual clock drift at the (k - 1)-th time point; ω θ，k and ω γ，k are two uncorrelated Gaussian random noise processes with zero mean, used to describe the random variations of clock offset and clock offset; the variances are respectively and ω θ，k-1 and ω γ，k-1 respectively represent the random noise processes of clock offset and clock offset at the (k - 1)-th time point.

[0093] In one example, in step S23, let x k = [θ k , γ k T represent the clock state vector at the k-th time point, and define ​For the measurement vector, according to the measurement equation and the recursive state equation, the state space model of PTP can be represented by a vector-matrix as follows:

[0094]

[0095] where w k-1 =[ω θ，k-1 , ω γ，k-1 T , v k =[v θ，k , v γ，k T . x k represents the clock state vector at the k-th time point; x k-1 represents the clock state vector at the (k - 1)-th time point; w k , v k are defined as Gaussian white noises with zero mean and uncorrelated with each other. w k-1 represents the Gaussian white noise at the (k - 1)-th time point; T represents vector transpose; v θ，k , v γ，k are Gaussian random noise processes with zero mean and finite variance; the state transition matrices A and C are:

[0096]

[0097] The measurement matrices D, E, H, and G are:

[0098]

[0099] In one example, in step S31, by adding the input λ k =[θ k , γ k T to the state vector x k , the extended state vector is represented as:

[0100]

[0101] In one example, in step S32, according to the discrete linear difference equation LDE of the delay asymmetry, the recurrence of λ k can be numerically described using PTP timestamp measurements. Considering the influence of noise, the delay asymmetry supplementary equation is obtained as:

[0102] λ k =λ k-1 +η k-1 +ω λ，k-1 ;

[0103] where η k-1 ​​​= [1, 1]α k + [-1, -1]β k represents the random noise in the time synchronization repetition process, and the relevant covariance matrix is Q λ .

[0104] In one example, the state equation and measurement equation of the extended state space model in step S33 are as follows:

[0105]

[0106] Among them, A a , C a , and D a are:

[0107] D a = [D E];

[0108] Among them, A a , C a both represent the extended state transition matrix; represents the extended state Gaussian white noise at the (k - 1)-th time point; z k represents the extended state measurement vector at the k-th time point; D a represents the extended state measurement matrix; H and G are both measurement matrices.

[0109] For the extended state space model, W k , ω λ，k and ν k are assumed to be Gaussian white noises with zero mean and uncorrelated with each other, and the noise and its relevant noise covariance matrix are given by the following formula:

[0110] W k ~ N(0, Q), ν k ~ N(0, R), ω λ，k ~ N(0, S)

[0111]

[0112] Among them, Q represents the process noise covariance matrix; R represents the observation noise covariance matrix; S represents the delay asymmetry noise covariance matrix. So far, the extended state space model has been obtained, and the adaptive Kalman filtering algorithm can be used to estimate the extended state vector.

[0113] In one example, in step S34, according to the extended state space model, an adaptive Kalman filtering algorithm for estimating the extended state vector is derived, and the adaptive Kalman filtering algorithm is used to estimate the clock offset in the PTP clock synchronization process, including:

[0114] Initialization process: Define the prior estimation covariance matrix of the extended state vector and the posterior estimation covariance matrix of the extended state vector:

[0115]

[0116] Among them, is the prior estimation of ; represents the posterior estimation of ; P k|k-1 represents the covariance matrix of the estimation error, and P k|k is defined as the covariance matrix of the estimation error.

[0117] Prediction process: Perform prior state estimation according to the state equation and update the prior estimation error covariance matrix. Among them, the prior state estimation and the estimation error covariance matrix P k|k-1 are expressed as:

[0118]

[0119] P k|k-1 = A a P k|k (A a ) T + Q a ;

[0120] Since the process noise w k is extended by ω λ , the process noise covariance matrix is

[0121] Adaptive adjustment: Adjust the measurement noise of the extended state space model according to the weighting coefficient and correct the difference between the measured value and the predicted value based on the prior state estimation according to the residual. Among them, the expression of the weighting coefficient d k is:

[0122]

[0123] In the formula, b is the forgetting factor, taking 0.9 to 0.95, which determines the value of the weighting coefficient d k . At the beginning, d k is close to 1, making the measurement noise mainly depend on the difference between the measured value and the predicted value; as the number of iterations increases, d k gradually decreases, making the measurement noise and variance gradually depend on the historical values and converge to a relatively stable value, so as to achieve the purpose of adaptively adjusting the measurement noise.

[0124] Furthermore, the residual r k, defined as the difference between the measured value and its predicted value, introducing d k can be expressed as:

[0125]

[0126] Update process: Calculate the Kalman gain, and obtain the posterior estimate and the posterior estimate error covariance matrix by updating the prior estimate and the measured value. Among them, the Kalman gain K k 、The calculation expressions for the posterior estimate and its error covariance matrix are:

[0127] K k = P k|k-1 (A a ) T (A a P k|k-1 (A a ) T + R k ) -1 ;

[0128]

[0129] P k|k = (I - K k A a )(P k|k-1 ) -1 ;

[0130] Among them, the measurement noise covariance matrix R k is updated in the following way:

[0131]

[0132] Clock offset estimation process: According to the clock offset θ in the posterior estimate k compensate the slave clock to complete clock synchronization. Specifically, in each iteration, the posterior estimate is an update of all elements in the extended state vector. Through the posterior estimate, the optimal estimate value of the clock offset in the extended state vector can be obtained, including clock offset, clock drift, and delay asymmetry. Using the sum of the three can achieve clock synchronization compensation.

[0133] In an example, after the step of estimating the clock offset in the PTP clock synchronization process using the adaptive Kalman filter algorithm, it further includes detecting the time stamp using residual awareness, including the following sub-steps:

[0134] Calculate the squared Mahalanobis distance of the residual, compare the squared Mahalanobis distance with the threshold related to the expected confidence level. If the squared Mahalanobis distance is greater than the threshold related to the expected confidence level, mark the time stamp as an outlier. Otherwise, the time stamp passes the verification and is used for the posterior estimate.

[0135] Specifically, when the clock synchronization system is not affected by outliers, the residual is a white noise process with a mean of zero, and its covariance matrix Γ k is:

[0136]

[0137] When the process noise is negligible compared to the measurement noise, Γ k depends only on the assumed measurement uncertainty. To detect potential outliers in the observed timestamps, the corresponding measurement vector z k and the resulting residual vector r k as well as its covariance matrix Γ k are used to calculate the squared Mahalanobis distance

[0138]

[0139] Therefore, the detection condition for timestamp verification is: where η a is the threshold related to the desired confidence level (1 - κ). r k can be shown to be a chi 2 random variable with two degrees of freedom because r k can be assumed to be approximately two-dimensional Gaussian distributed. Therefore, the threshold η a corresponding to the confidence level (1 - κ) is:

[0140] η κ = -2lnκ.

[0141] Whenever a new timestamp is observed, if the detection condition inequality is not satisfied, then z k is marked as an outlier; otherwise, the timestamp is verified and accepted, and at this time, the posterior estimate is updated as follows:

[0142]

[0143] This example obtains an accurate and robust clock offset estimation method based on the state space model. This method uses an adaptive Kalman filter with state extension for clock offset estimation and a residual-aware timestamp verification method to detect and preprocess unreliable timestamps. The specific implementation process is as follows:

[0144] S100: Establish a discrete-time clock synchronization state model based on the Precision Time Protocol (PTP);

[0145] S200: Establish a state space model for estimating clock offset and clock drift;

[0146] S300: Calculate the current measurement value z according to the extended state space model k , calculate the current prior estimate result according to the formula in the prediction stage and the previous posterior estimate result, and calculate the weighting coefficient d k and the residual r k ;

[0147] S400: Detect the timestamp using residual awareness, update based on the verified timestamp, and calculate the posterior estimate and the quantities required for the next iteration (posterior estimate error covariance matrix and Kalman gain);

[0148] S500: According to the posterior estimate of the clock offset θ k + clock drift γ k + delay asymmetry λ k compensate the slave clock, thereby completing the current clock synchronization.

[0149] The present invention also provides a computer program product, including a computer program, which when executed by a processor implements the steps of the above-mentioned clock synchronization method based on state expansion adaptive Kalman filtering formed by any one example or a combination of multiple examples. Wherein, the processor can be a single-core or multi-core central processing unit or a specific integrated circuit, or an integrated circuit configured to implement one or more of the present invention.

[0150] The present invention also provides an electronic device, which has the same inventive concept as any one example or a combination of multiple examples corresponding to the above-mentioned clock synchronization method based on state expansion adaptive Kalman filtering, including a memory and a processor, and a computer instruction that can run on the processor is stored on the memory, and when the processor runs the computer instruction, it executes the steps of the above-mentioned clock synchronization method based on state expansion adaptive Kalman filtering. The processor can be a single-core or multi-core central processing unit or a specific integrated circuit, or an integrated circuit configured to implement one or more of the present invention.

[0151] In one example, the electronic device is presented in the form of a general computing device, and the components of the electronic device may include, but are not limited to: the above-mentioned at least one processing unit (processor), the above-mentioned at least one storage unit, and a bus connecting different system components (including the storage unit and the processing unit).

[0152] Among them, the storage unit stores program codes, which can be executed by the processing unit, so that the processing unit executes the steps according to various exemplary embodiments of the present invention described in the above "Exemplary Method" section of this specification. For example, the processing unit can execute the above clock synchronization method based on state-expanded adaptive Kalman filtering.

[0153] The storage unit may include a readable medium in the form of a volatile storage unit, such as a random access storage unit (RAM) 3201 and / or a cache storage unit, and may further include a read-only storage unit (ROM).

[0154] The storage unit may also include a program / utility with a set (at least one) of program modules. Such program modules include, but are not limited to: an operating system, one or more application programs, other program modules, and program data. Each or some combination of these examples may include the implementation of a network environment.

[0155] The bus may represent one or more of several types of bus structures, including a storage unit bus or a storage unit controller, a peripheral bus, a graphics acceleration port, a processing unit, or a local bus using any bus structure in a variety of bus structures.

[0156] The electronic device can also communicate with one or more external devices (such as a keyboard, a pointing device, a Bluetooth device, etc.), can also communicate with one or more devices that enable a user to interact with the electronic device, and / or communicate with any device that enables the electronic device to communicate with one or more other computing devices (such as a router, a modem, etc.). Such communication can be carried out through an input / output (I / O) interface. And, the electronic device can also communicate with one or more networks (such as a local area network (LAN), a wide area network (WAN), and / or a public network, such as the Internet) through a network adapter. The network adapter communicates with other modules of the electronic device through the bus. It should be understood that other hardware and / or software modules can be used in combination with the electronic device, including but not limited to: microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data backup storage systems, etc.

[0157] Through the above description, those skilled in the art can easily understand that the exemplary embodiments described herein can be implemented by software or by a combination of software and necessary hardware. Therefore, the technical solution according to this exemplary embodiment can be embodied in the form of a software product, which can be stored in a non-volatile storage medium (which can be a CD-ROM, a USB flash drive, a mobile hard disk, etc.) or on a network, and includes several instructions to enable a computing device (which can be a personal computer, a server, an electronic device, or a network device, etc.) to execute the method of the exemplary embodiment of the present application.

[0158] The above specific embodiments are detailed descriptions of the present invention. It cannot be determined that the specific embodiments of the present invention are only limited to these descriptions. For those of ordinary skill in the technical field to which the present invention pertains, without departing from the concept of the present invention, several simple deductions and substitutions can still be made, and all should be regarded as belonging to the protection scope of the present invention.

Claims

1. A clock synchronization method based on state expansion adaptive Kalman filtering, characterized in that: The following steps are involved: A discrete time clock synchronization state model based on the precision time protocol PTP is established, including: modeling the clock offset estimation and clock drift estimation between the master and slave clocks in the synchronization process according to the PTP message exchange process; describing the change of delay in the clock offset estimation through discrete linear difference equations, and obtaining discrete linear difference equations for delay asymmetry; Establish a state space model for estimating clock offset and clock drift, including: rewriting the measurement equations for clock offset estimation and clock drift estimation based on the clock offset estimation and clock drift estimation models of the master and slave clocks, and providing a random variable representation of measurement noise related to clock offset and clock drift; deriving the recursive state equations for clock offset and clock drift based on stochastic differential equations; deriving the vector-matrix representation of the state space model of PTP based on the measurement equations and the recursive state equations; The invention relates to clock state estimation based on state extension adaptive Kalman filtering, including: adding input delay asymmetry to a state vector to obtain an extended state vector; using PTP timestamp measurement to numerically describe the recursion of delay asymmetry based on a discrete linear difference equation of delay asymmetry to obtain a delay asymmetry supplementary equation; establishing a state equation and a measurement equation of an extended state space model; deriving an adaptive Kalman filtering algorithm for estimating the extended state vector based on the extended state space model, and using the adaptive Kalman filtering algorithm to estimate the clock offset in the PTP clock synchronization process.

2. The clock synchronization method based on state expansion adaptive Kalman filtering according to claim 1 is characterized in that: The discrete linear difference equation for the delay asymmetry is: l k =λ k-1 +[1,1]a k +[-1,-1]b k ; a k =[Δt k ,Δt' k ]; β k =[ΔC(t k +d ms,k ),ΔC(t′ k -d sm,k )]; Among them, λ k is the delay asymmetry at the kth time point; k-1 is the delay asymmetry at the k-1th time point; Δt k Indicates the timestamp change between the synchronization message sent by the master clock at the kth time point and the k-1th time point; Δt' k represents the change in the receiving time difference recorded after the master clock receives the delay request message at the kth time point and the k-1th time point; ΔC(t k +d ms,k ) represents the synchronization message received from the clock at the kth time point and the k-1th time point and delayed by time d ms，k The change in the receiving time difference recorded later; ΔC(t′ k -d sm,k ) represents the reception time C(t k +d ms,k ) after responding to the delayed request message.

3. The clock synchronization method based on state expansion adaptive Kalman filtering according to claim 1, characterized in that: After the step of rewriting the measurement equations for clock offset estimation and clock drift estimation, the method further includes: The Bernoulli random process is introduced to optimize the clock offset measurement equation. The optimized clock offset measurement equation is: in, represents the estimated clock offset at the kth time point, θ k represents the actual clock offset at the kth time point; λ k is the delay asymmetry at the kth time point; n B,k represents the Bernoulli random process; v θ,k is a random variable representing the measurement noise associated with the clock offset.

4. The clock synchronization method based on state expansion adaptive Kalman filtering according to claim 1, characterized in that: The recursive state equations of the clock offset and clock drift are: Among them, θ k represents the actual clock offset at the kth time point, θ k-1 represents the actual clock offset at the k-1th time point; Δτ represents the time point; γ k represents the actual clock drift at the kth time point; γ k represents the actual clock drift at the k-1th time point; ω θ,k-1 ,ω γ,k-1 are two uncorrelated Gaussian random noise processes with zero mean.

5. The clock synchronization method based on state expansion adaptive Kalman filtering according to claim 1, characterized in that: The vector-matrix representation of the state-space model of the PTP is: w k-1 =[ω θ,k-1 ,oh γ,k-1 ] T ; n k =[n θ,k ,v γ,k ] T ; x k represents the clock state vector at the kth time point; x k-1 represents the clock state vector at the k-1th time point; w k , ν k Defined as Gaussian white noise with zero mean and no correlation, w k-1 represents the Gaussian white noise at the k-1th time point; z k represents the measurement vector; λ k is the delay asymmetry at the kth time point k; n B,k represents a Bernoulli random process; ω θ,k-1 ,ω γ,k-1 are two unrelated Gaussian random noise processes with zero mean; T represents vector transposition; ν θ,k , ν γ,k is a Gaussian random noise process with zero mean and finite variance; A and C are state transfer matrices; D, E, H and G are measurement matrices.

6. The clock synchronization method based on state expansion adaptive Kalman filtering according to claim 1, characterized in that: The state equation and measurement equation of the extended state space model are: in, represents the extended state vector at the kth time point; A represents the extended state vector at the k-1th time point; a , C a All represent extended state transfer matrices; represents the extended state Gaussian white noise at the k-1th time point; z k represents the extended state measurement vector at the kth time point; D a represents the extended state measurement matrix; H and G are both measurement matrices; n B,k represents the Bernoulli random process; ν k 、w k-1 ,ω λ,k-1 They are all Gaussian white noise with zero mean and no correlation with each other; T represents vector transpose.

7. The clock synchronization method based on state expansion adaptive Kalman filtering according to claim 1, characterized in that: The method of deriving an adaptive Kalman filter algorithm for estimating an extended state vector according to an extended state space model and using the adaptive Kalman filter algorithm to estimate a clock offset in a PTP clock synchronization process includes: Initialization processing: define the a priori estimated covariance matrix of the extended state vector and the a posteriori estimated covariance matrix of the extended state vector; Prediction processing: perform a priori state estimation based on the state equation and update the a priori estimation error covariance matrix; Adaptive adjustment: The extended state space model measurement noise is adjusted according to the weighting coefficients, and the difference between the measured value and the predicted value based on the prior state estimate is corrected according to the residual; Update processing: calculate the Kalman gain, and obtain the posterior estimate and the posterior estimate error covariance matrix by updating the prior estimate and the measurement value; Clock offset estimation processing: Compensate the slave clock based on the sum of clock offset, clock drift and delay asymmetry in the a posteriori estimate.

8. The clock synchronization method based on state extension adaptive Kalman filtering according to claim 1 or 7, characterized in that: After the step of estimating the clock offset in the PTP clock synchronization process using the adaptive Kalman filter algorithm, the method further includes detecting the timestamp using residual perception, including the following sub-steps: The squared Mahalanobis distance of the residual is calculated and compared with a threshold associated with the desired confidence. If the squared Mahalanobis distance is greater than the threshold associated with the desired confidence, the timestamp is marked as an outlier. Otherwise, the timestamp is validated and used for the posterior estimation.

9. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the clock synchronization method based on state extension adaptive Kalman filtering described in any one of claims 1 to 8 are implemented.

10. An electronic device comprising a memory and a processor, wherein the memory stores computer instructions that can be executed on the processor, wherein: When the processor runs the computer instructions, the processor performs the steps of the clock synchronization method based on state extension adaptive Kalman filtering as described in any one of claims 1-8.

Images