Beidou clock synchronization anomaly detection method for smart fusion terminal of transformer area

Abstract

The application belongs to the technical field of anomaly detection, and particularly relates to a Beidou clock synchronization anomaly detection method for a district intelligent fusion terminal, which comprises the following steps: collecting the Beidou time service original timestamp, the carrier-to-noise ratio and the environmental temperature of the intelligent fusion terminal; extracting parameter components through an iterative reweighted least squares method, establishing a clock long-term drift model and obtaining drift coefficients; removing the parameter components to obtain residual sequences; generating a kernel regression initial bandwidth by mapping the drift coefficients, and optimizing the bandwidth by taking the minimum mutual information of the two components as a criterion; superimposing the parameter component prediction value, the non-parameter component estimated value and the offset amount, and constructing a detection boundary; the offset amount changes with the confidence loss amount determined based on the historical window average weight and the non-parameter component curvature norm; and determining the Beidou clock synchronization anomaly and alarming when a new clock bias sample is out of bounds. The application can fuse multi-dimensional state features and adaptively adjust the detection boundary.

Description

Technical Field

[0001] This invention relates to the field of anomaly detection technology. More specifically, this invention relates to a method for detecting BeiDou clock synchronization anomalies in intelligent converged terminals for distribution areas. Background Technology

[0002] Intelligent converged terminals, as crucial devices for edge computing and on-site data acquisition, are widely used in various distributed business systems. In multi-node collaborative operation scenarios, high-precision time synchronization is a key foundation for ensuring data consistency, accurate event sequence recording, and the safe and stable operation of the system. BeiDou time synchronization technology can provide a unified time reference for distributed terminals, enabling the matching of acquired data and control commands on the same timeline.

[0003] Currently, intelligent fusion terminals are mostly deployed in outdoor transformer substations or complex industrial sites, and are constantly facing unstable factors such as temperature changes, electromagnetic interference, attenuation due to obstruction, and component aging. In actual operation, changes in the spatial environment, electromagnetic interference, and local obstruction can cause fluctuations in the quality of BeiDou satellite signals; changes in ambient temperature and aging of internal components can cause cumulative drift and random disturbances in the terminal's local clock.

[0004] However, existing clock anomaly detection mechanisms typically fail to fully integrate signal quality, environmental conditions, and clock evolution characteristics. They struggle to accurately characterize the evolution of clock deviations under varying operating conditions, and they also find it difficult to distinguish between normal environmental fluctuations and natural drift-induced deviations and genuine clock synchronization anomalies. This leads to false alarms or missed alarms, and their adaptive adjustment capabilities to changes in operating conditions are insufficient. Therefore, a BeiDou clock synchronization anomaly detection mechanism for intelligent fusion terminals is needed. This mechanism should integrate multi-dimensional environmental and state characteristics to characterize long-term clock drift patterns and random fluctuations, thereby improving the accuracy of BeiDou clock synchronization anomaly identification and detection. Summary of the Invention

[0005] To address the technical problem that existing clock anomaly detection methods struggle to distinguish between routine deviations and genuine anomalies, this invention proposes a BeiDou clock synchronization anomaly detection method for intelligent fusion terminals in distribution areas. This method can integrate multi-dimensional state features and adaptively adjust the detection boundary.

[0006] In a first aspect, the present invention provides a method for detecting BeiDou clock synchronization anomalies in intelligent fusion terminals in a distribution area, comprising: S1, acquiring the original timestamp of BeiDou timing, satellite signal carrier-to-noise ratio, and ambient temperature data of the intelligent fusion terminal; acquiring weighting coefficients and constructing an adaptive weighted model to generate a clock deviation sequence; calculating adjustment constants and processing the clock deviation sequence using an iterative reweighted least squares method to extract parameter components, obtaining a drift model and drift coefficients representing the long-term trend of the clock, wherein the drift coefficients are used to represent the inherent drift characteristics of the clock accumulated over time; S2, removing parameter components from the clock deviation sequence to obtain a residual sequence; applying kernel regression to the residual sequence to estimate nonparametric components, wherein the nonparametric components are used to represent the clock... The nonlinear random fluctuation characteristics; the initial bandwidth of kernel regression is generated based on the drift coefficient mapping, and fine-tuned to obtain the bandwidth by minimizing the mutual information between nonparametric and parametric components, thereby improving the estimation accuracy; S3, based on the sum of the predicted value of the parametric component and the estimated value of the nonparametric component as the benchmark, the offset is added or subtracted to set the detection boundary; the offset is proportional to the confidence missing amount determined based on the average weight within the historical time window, the confidence missing amount is negatively correlated with the average weight, and the offset is proportional to the curvature norm of the nonparametric component within the window, thereby realizing the adjustment for different states; when the new clock deviation sample point exceeds the detection boundary, it is determined that the current BeiDou clock synchronization is abnormal and the result is output.

[0007] By employing the above technical solutions, multi-dimensional data such as the original timestamp of BeiDou time synchronization, satellite signal carrier-to-noise ratio, and ambient temperature are integrated to construct an adaptive weighted model to generate clock deviation sequences. An iterative reweighted least squares method is used to extract the parametric components and drift coefficients characterizing the long-term clock trend, accurately depicting the cumulative clock drift characteristics. Simultaneously, a kernel regression method is combined to extract nonparametric components to represent nonlinear random fluctuation characteristics. The kernel regression bandwidth is optimized based on minimizing the mutual information between parametric and nonparametric components, effectively separating long-term drift from short-term fluctuations. In the detection phase, a detection benchmark is constructed by integrating the predicted values ​​of parametric components and the estimated values ​​of nonparametric components. The detection boundary is adaptively adjusted based on the historical window average weight and the curvature norm of the nonparametric components, enabling the anomaly judgment criteria to dynamically match the signal quality and clock fluctuation state, thereby improving the accuracy of BeiDou clock synchronization anomaly detection and reducing the risk of false alarms and missed alarms.

[0008] Preferably, the step of constructing an adaptive weighted model to generate a clock deviation sequence includes: subtracting the original BeiDou time stamp from the local reference time stamp to obtain a basic clock deviation sequence; obtaining the weight coefficient for each time point, wherein the weight coefficient has a positive nonlinear relationship with the amplitude of the satellite signal carrier-to-noise ratio at the current time; using the weight coefficient as the confidence weight of the corresponding data point in the basic clock deviation sequence, performing weighted smoothing, weighted robust estimation, or weighted filtering on the basic clock deviation sequence to obtain a signal quality corrected basic clock deviation sequence, wherein the larger the weight coefficient, the greater the contribution of the corresponding data point to the correction result, and the smaller the weight coefficient, the smaller the contribution of the corresponding data point to the correction result; and adding the corresponding elements of the adjusted basic clock deviation sequence and the time delay increment sequence to generate a clock deviation sequence.

[0009] By adopting the above technical solution, the basic clock deviation sequence is obtained by subtracting the original timestamp of BeiDou time synchronization from the local reference timestamp. Then, using the weighting coefficients that have a positive nonlinear relationship with the satellite signal carrier-to-noise ratio as confidence weights, the basic clock deviation sequence is subjected to weighted smoothing, weighted robust estimation, or weighted filtering. This makes the contribution of high carrier-to-noise ratio samples to the correction result greater and the contribution of low carrier-to-noise ratio samples smaller, thereby effectively suppressing pseudorange noise introduced by signal quality fluctuations and improving the signal quality correction effect of the clock deviation sequence.

[0010] Preferably, obtaining the weighting coefficients includes: obtaining the satellite signal carrier-to-noise ratio amplitude at the current moment; determining the current signal quality characterization value based on the relative relationship between the satellite signal carrier-to-noise ratio amplitude and a preset reference carrier-to-noise ratio amplitude; generating corresponding weighting coefficients according to a preset mapping relationship between carrier-to-noise ratio and weighting coefficients; wherein the mapping relationship is a piecewise mapping relationship, a lookup table mapping relationship, or a monotonically increasing function mapping relationship, and the weighting coefficients increase with the increase of the satellite signal carrier-to-noise ratio amplitude. The weighting coefficients are used to characterize the reliability of the baseline clock deviation data points at the corresponding moment and to determine the contribution of the corresponding data points in the weighted processing.

[0011] By adopting the above technical solution, the current satellite signal carrier-to-noise ratio amplitude is obtained, and the signal quality characterization value is determined according to its relative relationship with the preset reference carrier-to-noise ratio amplitude. Then, the weighting coefficients are generated according to the preset mapping relationship between carrier-to-noise ratio and weighting coefficients, so that the weighting coefficients can be adaptively adjusted as the signal quality changes. By using methods such as segmented mapping, lookup table mapping, or monotonically increasing function mapping, it is ensured that the weighting coefficients accurately reflect the reliability of the basic clock deviation data points at the corresponding time, providing a reliable confidence basis for subsequent weighted processing.

[0012] Preferably, the step of processing the clock deviation sequence using the iterative reweighted least squares method to extract parameter components and obtain a drift model and drift coefficients representing the long-term trend of the clock includes: constructing a long-term trend fitting model containing a constant term, a first-order trend term, and a second-order trend term based on the clock deviation sequence; performing initial estimation of each model parameter of the long-term trend fitting model according to the clock deviation sequence; extracting the converged second-order trend term parameters as the drift coefficients; generating a corresponding trend sequence based on the converged long-term trend fitting model, and using the trend sequence as a parameter component representing the long-term trend of the clock signal.

[0013] By adopting the above technical solution, a long-term trend fitting model containing constant terms, first-order trend terms, and second-order trend terms is constructed. After initial estimation of the model parameters, the converged second-order trend term parameters are extracted as drift coefficients. This can effectively separate the parameter components representing the inherent drift characteristics of the clock over time from the clock deviation sequence. The generated trend sequence can accurately characterize the long-term drift law of the clock, providing a reliable benchmark reference for subsequent residual analysis and anomaly detection.

[0014] Preferably, after initially estimating the various model parameters of the long-term trend fitting model, the method further includes: determining the weight of each data point in the current round based on a bisquared weight function or other stable weight function, and iteratively updating the model parameters using a weighted fitting method until the change in model parameters between two adjacent iterations is lower than a preset convergence threshold.

[0015] By adopting the above technical solution, the weights of each data point are determined by a bisquared weighting function or other stable weighting functions, and the model parameters are iteratively updated until convergence by weighted fitting. This can effectively suppress the interference of gross errors and outliers in the clock deviation sequence on the long-term trend fitting, and improve the fitting robustness and parameter estimation accuracy of the drift model under complex working conditions.

[0016] Preferably, the calculation of the adjustment constant includes: performing subtraction operations on the data points of two adjacent time points in the clock deviation sequence according to the time order to generate a first-order difference sequence; sorting all the values ​​in the first-order difference sequence in ascending order; subtracting the value corresponding to the 25th percentile from the value corresponding to the 75th percentile in the sorted first-order difference sequence to calculate the interquartile range; dividing the interquartile range by 1.349 and multiplying by 4.685 to obtain the adjustment constant of the current bisquared weighting function.

[0017] Preferably, the initial bandwidth of the kernel regression is generated based on the drift coefficient mapping, including: extracting the drift coefficient and calculating its absolute value; calculating the absolute value of the drift coefficient and adding the absolute value to a preset small positive bias to obtain a drift correction amount; dividing a preset first constant by the drift correction amount to obtain a dimensionless inverse factor; wherein the first constant and the small positive bias have the same dimension as the drift coefficient; multiplying a preset scaling factor constant with time dimension by the dimensionless inverse factor, and setting the resulting value as the initial bandwidth required for kernel regression.

[0018] Preferably, the step of fine-tuning the bandwidth based on minimizing the mutual information between nonparametric and parametric components includes: generating multiple candidate bandwidths based on the initial bandwidth and a preset search step size; performing kernel regression estimation on the residual sequence using each candidate bandwidth to obtain candidate nonparametric component sequences corresponding to each candidate bandwidth; calculating the mutual information evaluation value between each candidate nonparametric component sequence and the parametric component sequence according to a preset mutual information evaluation model; selecting the candidate bandwidth with the smallest mutual information evaluation value as the fine-tuned bandwidth; and performing kernel regression estimation on the residual sequence based on the fine-tuned bandwidth to obtain the nonparametric components.

[0019] Preferably, the offset is proportional to the number of missing confidence levels determined based on the average weight within the historical time window, the number of missing confidence levels is negatively correlated with the average weight, and the offset is proportional to the curvature norm of the nonparametric component within the window. This includes: assigning convergence weights to each data point within the historical time window using a statistical iterative reweighted least squares method, and calculating the average weight value; determining the number of missing confidence levels based on the average weight value, where a larger average weight value results in a smaller number of missing confidence levels, and a smaller average weight value results in a larger number of missing confidence levels; obtaining curvature feature values ​​to characterize the intensity of fluctuations in the nonparametric component; and determining the offset based on the number of missing confidence levels and the curvature feature values ​​according to a preset monotonic mapping relationship, a piecewise mapping relationship, or a lookup table mapping relationship; where a larger number of missing confidence levels results in a larger offset; and a larger curvature feature value results in a larger offset.

[0020] Preferably, obtaining the curvature feature value used to characterize the intensity of fluctuations in the nonparametric components includes: performing differential processing on the nonparametric component sequence to obtain the curvature feature value.

[0021] The present invention has the following beneficial effects:

[0022] This invention integrates multi-dimensional data such as the original timestamp of BeiDou time synchronization, satellite signal carrier-to-noise ratio, and ambient temperature to construct an adaptive weighted model to generate a clock deviation sequence. It uses an iterative reweighted least squares method for robust fitting to extract parameter components and drift coefficients, characterizing the cumulative clock drift characteristics. Furthermore, it employs a kernel regression method to optimize bandwidth based on minimizing mutual information, extracting the nonlinear random fluctuation characteristics of the clock. In the anomaly detection stage, it fuses the predicted values ​​of parameter components and the estimated values ​​of nonparametric components to construct a central benchmark, thereby achieving accurate detection of BeiDou clock synchronization anomalies.

[0023] Furthermore, this invention adaptively adjusts the detection boundary based on the confidence level missing quantity and the nonparametric component curvature norm determined by the historical average weight, so that the anomaly judgment criteria dynamically match the changes in signal quality and clock fluctuation state, thereby reducing the risk of false alarms and missed alarms and improving the system's operational stability. Attached Figure Description

[0024] Figure 1 This is a flowchart of a BeiDou clock synchronization anomaly detection method for intelligent converged terminals in a distribution area; Figure 2 This is a graph showing the fitting of long-term drift parameters; Figure 3 It is a bandwidth mutual information optimization graph; Figure 4 These are comparison images from ablation experiments. Detailed Implementation

[0025] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0026] This invention discloses a method for detecting BeiDou clock synchronization anomalies in intelligent converged terminals for distribution areas, referring to... Figure 1 This includes steps S1-S3: S1. Obtain the original timestamp of BeiDou time synchronization, satellite signal carrier-to-noise ratio and ambient temperature data of the intelligent fusion terminal, obtain the weighting coefficients and construct an adaptive weighted model to generate a clock deviation sequence; calculate the adjustment constant and use the iterative reweighted least squares method to process the clock deviation sequence to extract parameter components, and obtain the drift model and drift coefficient representing the long-term trend of the clock. The drift coefficient is used to represent the inherent drift characteristics of the clock accumulated over time.

[0027] The microcontroller's serial communication interface is used to call the underlying serial port read function to parse the data stream of the national universal receiver data protocol, extracting the standard second-level time data, visible satellite numbers, and satellite signal carrier-to-noise ratio information output by the BeiDou positioning module. Simultaneously, the microcontroller's timer input capture interface is used to acquire the rising edge of the 1PPS pulse output by the BeiDou positioning module. This 1PPS pulse rising edge is then hardware-captured and aligned with the local nanosecond-level count value generated by the microprocessor's high-frequency crystal oscillator counter. Subtraction is performed based on pre-calibrated antenna feeder delay, receiver output delay, and input capture channel delay to obtain the fundamental time difference between the nanosecond-level BeiDou timing original timestamp and the local reference timestamp. The standard second-level time data obtained from the serial port parsing is used to determine the absolute second epoch, and the 1PPS hardware capture result is used to determine the nanosecond-level deviation within that second epoch.

[0028] The hardware abstraction layer interface is called to read the ambient temperature value reported by the analog-to-digital converter channel of the motherboard thermistor, as well as the average carrier-to-noise ratio of each visible satellite in the baseband signal processor memory register, as the carrier-to-noise ratio.

[0029] Weighting coefficients are calculated using a weighting function based on a monotonic mapping of the carrier-to-noise ratio (CNR). Weighting coefficients between 0 and 1 are generated according to the relative relationship between the current CNR and a preset reference CNR amplitude; the higher the CNR, the larger the weighting coefficient, and vice versa. Ambient temperature data is not directly used in the weighting coefficient calculation; instead, it is used to generate time delay increments through a quadratic parabolic temperature-time characteristic equation. The weighting coefficients serve as confidence weights for the corresponding data points in weighted smoothing, weighted robust estimation, or weighted filtering processes, controlling the contribution of each data point to the correction result, rather than directly performing division amplification on single-point deviations. This constructs an adaptive weighted model, generating a continuous clock skew sequence.

[0030] When extracting parameter components, a robust scaling method based on first-order difference interquartile range (ICR) is used to calculate the adjustment constant, and the Tukey bisquared weighting function is used as the stable weighting function for the iterative reweighted least squares method. Specifically, first-order difference is performed on the clock skew sequence, the ICR of the first-order difference sequence is calculated, and the ICR is divided by 1.349 and multiplied by 4.685 to obtain the adjustment constant of the bisquared weighting function. Subsequently, second-order polynomial long-term trend fitting is performed on the clock skew sequence until the change in the Euclidean norm of the coefficient vector obtained from two adjacent iterations is less than a preset convergence threshold, at which point the iteration stops. The second-order polynomial coefficient vector output by the optimizer includes a constant term, a first-order trend term, and a second-order trend term, where the converged second-order trend term parameters serve as drift coefficients. This second-order polynomial equation constitutes a drift model representing the long-term trend of the clock.

[0031] In some implementations, an adaptive weighted model is constructed to generate the clock skew sequence, including: Subtracting the original timestamp from the BeiDou time synchronization from the local reference timestamp yields the basic clock offset sequence; The ambient temperature data is input into the preset quadratic parabolic temperature-time characteristic equation to calculate and obtain the time delay increment sequence corresponding to the ambient temperature. The weighting coefficients at each time point are obtained, and the weighting coefficients have a positive nonlinear relationship with the amplitude of the satellite signal carrier-to-noise ratio at the current time. The weighting coefficients are used as confidence weights for the corresponding data points in the basic clock deviation sequence. Weighted smoothing, weighted robust estimation, or weighted filtering are then applied to the basic clock deviation sequence to obtain the signal quality corrected basic clock deviation sequence. The larger the weighting coefficient, the greater the contribution of the corresponding data point to the correction result; the smaller the weighting coefficient, the smaller the contribution of the corresponding data point to the correction result. The adjusted base clock skew sequence is added element-wise to the time delay increment sequence to generate the clock skew sequence.

[0032] When acquiring the timestamp sequence, a 1-second sampling period is used. The BeiDou time synchronization second epoch is determined via serial port messages. The nanosecond-level edge arrival time corresponding to this second epoch is determined by capturing 1PPS hardware input. This time is then subtracted from the reference count timestamp generated by the local temperature-controlled crystal oscillator to obtain the basic clock deviation sequence. The initial deviation value at a single point is typically distributed in the range of -1000 to 1000 nanoseconds. For example, subtracting at a certain moment yields a basic clock deviation of 250.5 nanoseconds.

[0033] Simultaneously, ambient temperature data in degrees Celsius output from the terminal's built-in temperature sensor is collected and input into a preset quadratic parabolic temperature-time characteristic equation. In the middle, the reference reversal temperature The preferred setting is 25°C, with a temperature coefficient of 100°C. The empirical value ranges from 0.01 to 0.1. Between, the example is taken Zero-degree bias It is usually set to 0 nanoseconds. If the current ambient temperature T is 35℃, then the calculated time delay increment for the current moment is 5 nanoseconds.

[0034] By analyzing the GSV statements in the BeiDou NMEA-0183 protocol message, the carrier-to-noise ratio (CNR) amplitude of the satellite signal at the corresponding time is extracted, typically ranging from 30 to 50 dB-Hz. A nonlinear exponential mapping model is then used to calculate weighting coefficients distributed between 0 and 1. For example, when the CNR is 45 dB-Hz and the reference CNR amplitude is also 45 dB-Hz, the relative ratio R is 1. The calculated weighting coefficient is approximately 0.6321.

[0035] When performing sequence numerical correction, the weighting coefficients are used as confidence weights for weighted filtering. Taking a local window of length m as an example, the base clock deviation values ​​within the window are multiplied by their corresponding weights, summed, and then divided by the sum of the weights within the window to obtain the base clock deviation values ​​after signal quality correction. Among these, samples with high carrier-to-noise ratios contribute more to the correction result due to their larger weights, while samples with low carrier-to-noise ratios contribute less to the correction result due to their smaller weights.

[0036] For example, within the current local window, the base clock deviation at a certain moment is 250.5 nanoseconds and the corresponding weight is approximately 0.6321. This sample is used as a high-confidence observation for weighted filtering. After weighted estimation by multiple samples within the window, the base clock deviation value after signal quality correction at that moment is 251.2 nanoseconds.

[0037] The adjusted deviation value is added point by point to the corresponding time delay increment of 5 nanoseconds to generate a clock deviation sequence value of 256.2 nanoseconds that includes dual corrections for temperature and signal quality. This step is repeated second by second to construct a complete continuous clock deviation time series.

[0038] In some implementations, the weight coefficients are obtained, including: Obtain the current satellite signal carrier-to-noise ratio amplitude; The current signal quality characterization value is determined based on the relative relationship between the satellite signal carrier-to-noise ratio amplitude and the preset reference carrier-to-noise ratio amplitude. Generate the corresponding weighting coefficients according to the preset mapping relationship between carrier-to-noise ratio and weighting coefficients; The mapping relationship can be a piecewise mapping relationship, a lookup table mapping relationship, or a monotonically increasing function mapping relationship. The weight coefficient increases with the increase of the satellite signal carrier-to-noise ratio amplitude. The weight coefficient is used to characterize the reliability of the basic clock deviation data point at the corresponding time and to determine the contribution of the corresponding data point in the weighted processing.

[0039] In the weighting coefficient generation stage, a preset reference carrier-to-noise ratio (CNR) amplitude is configured. This reference amplitude represents the nominal signal quality of the BeiDou terminal in an ideal open environment. The preferred parameter range is set between 40 and 50 dB-Hz, typically 45 dB-Hz. When the system operates at a data refresh rate of 1 Hz, the effective satellite signal CNR amplitude at the current moment is read, and the ratio of the current satellite signal CNR amplitude to the preset reference CNR amplitude is calculated to obtain a dimensionless relative ratio R. If the CNR amplitude at a certain observation epoch is 36 dB-Hz, then 36 is divided by the reference amplitude of 45 to calculate the relative ratio R = 0.8. This relative ratio R is used to characterize the degree of attenuation of the current channel state compared to the nominal state.

[0040] Based on this, according to the formula The weighting coefficients are calculated, where W is the calculated weighting coefficient, e is the natural constant, and R is the relative ratio. Substituting R=0.8 into the above nonlinear function, the natural constant e is approximated as 2.71828, and the calculation process is e^(-0.8)≈0.4493. Subtraction yields the weighting coefficient W=0.5507. This function mapping mechanism ensures that the weighting coefficient W lies within the interval (0,1). When the signal is in a good state, such as a current carrier-to-noise ratio of 45dB-Hz, and R=1, the weighting coefficient rises to approximately 0.6321. However, when the signal is affected by obstruction, such as a current carrier-to-noise ratio dropping to 22.5dB-Hz, and R=0.5, the weighting coefficient decreases to 0.3934. This nonlinear transformation reduces the interference of timestamp pseudorange noise caused by signal jumps on the weighted model.

[0041] In some implementations, iterative reweighted least squares is used to process the clock skew sequence to extract parameter components, resulting in a drift model and drift coefficients representing the long-term trend of the clock, including: A long-term trend fitting model containing a constant term, a primary trend term, and a secondary trend term is constructed based on the clock deviation sequence. Initial estimates of the various model parameters of the long-term trend fitting model are made based on the clock deviation sequence; The weights of each data point in the current round are determined based on the bisquared weight function or other stable weight function, and the model parameters are iteratively updated using a weighted fitting method until the change in model parameters between two adjacent iterations is lower than the preset convergence threshold. Extract the parameters of the converged quadratic trend term as the drift coefficient; The corresponding trend sequence is generated based on the converged long-term trend fitting model, and the trend sequence is used as a parameter component representing the long-term trend of the clock signal.

[0042] In the calculation process for extracting parameter components, a clock deviation sequence with a pre-selected length of 3600 sample points is used to form the analysis window, where the time variable t monotonically increases from 1 to 3600 in seconds. The drift model representing the long-term trend of the clock is set as a second-order polynomial equation containing a constant term: Where D is the predicted parametric component, and t is the time independent variable. , and These are the coefficients of the zero-degree, first-degree, and second-degree polynomials, respectively.

[0043] The initial values ​​of the three polynomial coefficients of the second-order polynomial equation are obtained using ordinary least squares (OLS) based on the clock skew sequence. Specifically, a matrix equation of the form Y=XA is established, where Y represents the clock skew observation vector within the analysis window, and X represents the design matrix consisting of time variables. Represents the vector of polynomial coefficients to be solved, and through The initial coefficient values ​​are obtained. From this, the zero-order constant term is obtained. The corresponding initial phase deviation, for example, is 300 nanoseconds; the first term The corresponding frequency offset, such as a value of 0.05 ns / s; and the quadratic term. The corresponding long-term frequency drift rate, such as the value .

[0044] Starting with the initial coefficients mentioned above, the weights of the data points in the current round are determined based on the Tukey bisquared weighting function, and the polynomial coefficients are updated using an iterative reweighted least squares method. For each sample residual in the current k-th iteration... ,when When the value is less than the adjustment constant, the sample is assigned a corresponding weight according to the bisquared weighting function; when... When the weight of the sample point is not less than the adjustment constant, the weight of the sample point is truncated to 0 to suppress the influence of gross errors on the long-term drift fitting results.

[0045] In each iteration update, the data point weight matrix is ​​based on a diagonal matrix. Perform weighted least squares solution Update the polynomial coefficient vector. Calculate the L2 norm of the difference between the coefficient vectors of two consecutive iterations. When the norm result is lower than the preset convergence threshold, the preferred range is defined as follows: to Between, for example, set as Alternatively, the iteration loop will terminate when the maximum number of iterations reaches 50.

[0046] After determining that the model has converged, extract the quadratic term from the coefficient vector at the stopping time. For example, converges to This is designated as the drift coefficient. Further, the coefficients confirmed by convergence are... , , Substitute back to the second-order polynomial equation A parabolic sequence was generated point-by-point for the independent variable t=1,2,...,3600 within the window to establish the parametric components representing the long-term drift phenomenon. The fitting effect of this drift model is as follows: Figure 2 As shown, this reflects the long-term cumulative trend of clock deviation over time, as well as the model's fitting accuracy to the actual observation data.

[0047] In some implementations, the adjustment constant is calculated, including: For the clock skew sequence, subtract the data points of two adjacent moments in chronological order to generate a first-order difference sequence; Sort all values ​​in the first-order difference sequence in ascending order; The interquartile range is calculated by subtracting the value corresponding to the 25th percentile from the value corresponding to the 75th percentile in the sorted first-order difference sequence. Divide the interquartile range by 1.349 and multiply by 4.685 to obtain the adjustment constant of the current bisquared weighting function.

[0048] Load the full clock offset sequence data within the time window. For a sequence of total length N, such as a set of 3600 sampling points, use the clock offset value at the current time point i (e.g., the offset at the 100th second is 310 nanoseconds) and subtract the value at the previous time point i-1 (e.g., the offset at the 99th second is 308 nanoseconds). This generates a first-order difference sequence of total length N-1, or 3599 points. This sequence physically represents the short-term stable frequency offset jump characteristics between adjacent seconds. Then, arrange all floating-point values ​​in this first-order difference sequence in ascending order from smallest to largest.

[0049] Based on this ordered sequence, obtain the lower and upper quartiles, specifically the 25th percentile at the 0.25×(N-1)th position (e.g., 1.2 nanoseconds) and the 75th percentile at the 0.75×(N-1)th position (e.g., 3.6 nanoseconds). Subtract 2.4 nanoseconds to obtain the interquartile range. Divide this interquartile range (2.4 nanoseconds) by 1.349 to calculate the stable standard deviation. Nanosecond.

[0050] Multiplying this standard deviation by the tuning factor 4.685, we calculate the adjustment constant c ≈ 8.33 nanoseconds applicable to the current bisquared weighting function. This adjustment constant serves as the threshold for residual filtering, ensuring that outlier noise weights exceeding this value are reset to zero.

[0051] S2, remove the parametric components from the clock deviation sequence to obtain the residual sequence, apply kernel regression to the residual sequence to estimate the nonparametric components, the nonparametric components are used to represent the nonlinear random fluctuation characteristics of the clock; the initial bandwidth of the kernel regression is generated based on the drift coefficient mapping, and fine-tuned to obtain the bandwidth by minimizing the mutual information between the nonparametric components and the parametric components, thereby improving the estimation accuracy.

[0052] The array subtraction module in the data processing program is called to subtract the predicted value sequence calculated by substituting the time variable into the drift model from the clock deviation sequence generated in the previous step point by point according to the time index. The resulting discrete difference sequence is the residual sequence.

[0053] The Nadalaya-Watson kernel regression algorithm was used to smooth the residual sequence and estimate the nonparametric components. The Nadalaya-Watson kernel regression algorithm is a nonparametric smoothing method based on kernel function weights. By assigning higher weights to neighboring samples and lower weights to distant samples, a weighted average estimate is performed on the residual sequence to extract continuously smoothed nonparametric components. When determining the initial bandwidth of the kernel regression, the second-order polynomial drift coefficients calculated in the preceding sequence were extracted. The absolute value of the drift coefficient is calculated, and a dimensionless inverse factor is generated based on a preset inverse mapping relationship. This dimensionless inverse factor is then multiplied by a scaling factor constant with a time dimension to obtain the initial bandwidth value. This ensures that a larger absolute value of the drift coefficient results in a smaller initial bandwidth, reducing the masking of short-term fluctuations by excessive smoothing under strong drift conditions; conversely, a smaller absolute value of the drift coefficient results in a larger initial bandwidth, enhancing noise suppression capabilities under stable conditions.

[0054] During the bandwidth fine-tuning optimization phase, multiple candidate bandwidths are generated near the initial bandwidth based on the initial bandwidth and a preset search step size. Kernel regression estimation is then performed on the residual sequence using each candidate bandwidth to obtain candidate nonparametric component sequences corresponding to each bandwidth. Based on a preset mutual information evaluation model, the mutual information evaluation value between each candidate nonparametric component sequence and the parametric component sequence is calculated. The candidate bandwidth with the smallest mutual information evaluation value is selected as the fine-tuned bandwidth. The Gaussian kernel function is then reinitialized using this fine-tuned bandwidth, and a weighted average mapping is performed on the residual sequence to complete the estimation of the nonparametric components.

[0055] In some implementations, the initial bandwidth of kernel regression is generated based on a drift coefficient mapping, including: Extract the drift coefficient and calculate its absolute value; Calculate the absolute value of the drift coefficient and add it to a preset small positive bias to obtain the drift correction amount; divide the preset first constant by the drift correction amount to obtain the dimensionless inverse factor; wherein, the first constant and the small positive bias have the same dimensions as the drift coefficient. Multiply the preset scaling factor constant with time dimension by the dimensionless inverse factor, and set the resulting value as the initial bandwidth required for kernel regression.

[0056] In establishing the nonparametric estimation smoothness benchmark, the second-order polynomial drift coefficients of the output are extracted. And use the absolute value function to obtain its value. Assume the extracted drift coefficient... for nanoseconds / square seconds, after absolute value conversion, become Nanoseconds per square second. The absolute value is offset from a predefined, small positive number in the denominator, typically a value of [value missing]. Nanoseconds and square seconds, ensuring they have the same physical dimensions, are added together to obtain the corrected denominator value. Nanoseconds per square second.

[0057] Using a first constant that also has the dimension of drift rate, this constant is usually set at... Magnitude, for example, selecting Nanoseconds per square second. Divide the first constant by the corrected denominator to calculate the quotient. This yields a dimensionless inverse factor of approximately 3.996. The purpose of this inverse factor is to reduce the impact of over-smoothing by outputting a smaller inverse ratio when the clock itself has a large internal drift.

[0058] A scaling factor constant with a time dimension, typically configured in the range of 2 to 10 seconds (e.g., fixed at 3 seconds), is multiplied by the dimensionless inverse factor 3.996. The resulting product, 11.988 seconds, is used as the initial time bandwidth parameter required for kernel regression. In actual calculations, the initial time bandwidth parameter is limited to a preset minimum and maximum bandwidth to avoid either an excessively large drift coefficient leading to insufficient bandwidth or a drift coefficient approaching zero leading to excessive bandwidth.

[0059] In some implementations, bandwidth is obtained by fine-tuning based on minimizing the mutual information between nonparametric and parametric components, including: Based on the initial bandwidth and the preset search step size, multiple candidate bandwidths are generated; Kernel regression estimation is performed on the residual sequence using each candidate bandwidth to obtain the candidate nonparametric component sequence corresponding to each candidate bandwidth; Based on the preset mutual information evaluation model, calculate the mutual information evaluation value between each candidate nonparametric component sequence and the parametric component sequence; The candidate bandwidth with the smallest mutual information evaluation value is selected as the fine-tuned bandwidth. Based on the fine-tuned bandwidth, kernel regression estimation is performed on the residual sequence to obtain nonparametric components.

[0060] During the bandwidth optimization process, assuming an initial time bandwidth of 12 seconds and a discrete search step size of 0.5 seconds, five exploration nodes are equidistantly derived in both the smaller and larger directions, centered at 12 seconds. This generates a one-dimensional arithmetic sequence with bandwidths ranging from 9.5 seconds to 14.5 seconds, totaling 11 candidate bandwidths. For the clock residual sequence after removing the effects of long-term parameter drift, a kernel regression algorithm is sequentially applied with the aforementioned 11 candidate bandwidths of different sizes to perform smooth fitting, generating 11 corresponding sets of candidate nonparametric component sequences.

[0061] For each set of candidate nonparametric component sequences, a two-dimensional integral plane is constructed, and a two-dimensional kernel density estimation algorithm based on the Gaussian kernel function is used to analyze the joint statistical characteristics between the candidate nonparametric component sequence X and the parametric component sequence D, thereby obtaining the joint probability distribution density matrix of the two-dimensional plane. And the independent marginal probability distribution density vector obtained by integrating along a single axis. and Where X represents the candidate nonparametric component sequence, and D represents the corresponding parametric component sequence. Let X and D represent the joint probability density function. and Let each represent an independent marginal probability distribution density function of the two components.

[0062] Subsequently, based on the mutual information calculation formula The mutual information value I under each candidate condition is calculated using numerical integration with fixed spatial step sizes dx and dd. For example, the mutual information value is 0.42Nat for a 9.5-second bandwidth, while the mutual information value calculated for an 11-second bandwidth is the lowest at 0.18Nat.

[0063] During the bandwidth selection process, the mutual information value corresponding to each candidate bandwidth and the fitting error of each candidate nonparametric component sequence relative to the residual sequence are calculated simultaneously. When the fitting error corresponding to a candidate bandwidth does not exceed a preset error threshold, it is considered a valid candidate bandwidth. Among all valid candidate bandwidths, the group of candidate bandwidths with the smallest mutual information value is selected as the bandwidth after fine-tuning and confirmation. For example, in the above example, if the mutual information value corresponding to the 11-second bandwidth is 0.18Nat and the residual fitting error does not exceed the preset error threshold, then 11 seconds is determined as the final bandwidth, and the nonparametric component sequence calculated from this bandwidth is output. The correspondence between bandwidth and mutual information is as follows: Figure 3 As shown, mutual information first decreases and then increases with increasing bandwidth, with the minimum value corresponding to the optimal bandwidth.

[0064] S3 sets the detection boundary by adding or subtracting the offset based on the sum of the predicted values ​​of the parametric components and the estimated values ​​of the non-parametric components. The offset is proportional to the amount of missing confidence determined based on the average weight within the historical time window, and the amount of missing confidence is negatively correlated with the average weight. Furthermore, the offset is proportional to the curvature norm of the non-parametric components within the window, thereby achieving adjustment for different states. When a new clock deviation sample point exceeds the detection boundary, it is determined that the current BeiDou clock synchronization is abnormal and the result is output.

[0065] The one-dimensional array of time-aligned parametric component predictions and the one-dimensional array of non-parametric component estimates are added element-wise to construct a central baseline representing the normal operating reference state of the clock. During real-time detection, the central baseline and detection boundary are generated based on the historical time window preceding the current sample point to be detected. The current sample point to be detected does not participate in the calculation of its own corresponding detection boundary to avoid abnormal samples being absorbed by the model and affecting the out-of-bounds judgment.

[0066] A sliding historical time window with a length of 120 consecutive sampling periods is opened in system memory. The arithmetic mean of the weights assigned to all time nodes within this sliding window by the preceding stable regression algorithm in the last iteration is calculated to obtain the average weight. The average weight is the data point fitting weight generated after the convergence of the iterative reweighted least squares method, used to characterize the fit consistency of the observed data with respect to the long-term drift model; this average weight differs from the weight coefficients obtained in S1 by mapping the satellite signal carrier-to-noise ratio.

[0067] Simultaneously, a second-order central difference operation is performed on the sequence data of nonparametric components within the same historical window to approximate the second-order derivative sequence of the curve. A sampling period Δt is introduced into the second-order central difference operation, i.e., according to... Calculate the second-order difference result so that the curvature norm has the dimension of clock deviation divided by the square of time. Solve for the L2 norm of this second-order derivative sequence as the curvature norm representing the intensity of the fluctuation.

[0068] A fixed constant is preset in the configuration file as the baseline ratio coefficient, and the confidence level missing value is calculated based on the average weight. The preferred confidence level missing value is q=1- ,in The average weight is calculated within the historical time window. The baseline scaling factor, the missing confidence level, and the curvature norm are multiplied by the microprocessor's floating-point unit to obtain the candidate offset value. This candidate offset value is then compared to a preset minimum detection boundary width, and the larger of the two is taken as the actual offset. Therefore, when the average weight is larger and the fitting confidence level of the observed data relative to the long-term drift model is higher, the missing confidence level is smaller, and the detection boundary is relatively narrower; conversely, when the average weight is smaller and the observation consistency is lower, the missing confidence level is larger, and the detection boundary is relatively wider.

[0069] The upper detection boundary array is generated by adding the actual offset to the center baseline. The lower detection boundary array is generated by subtracting the actual offset from the center baseline.

[0070] In the real-time monitoring loop, the central processing unit uses logic comparison instructions to perform real-time threshold judgment on the latest clock deviation sample point obtained by the Beidou time synchronization input acquisition module or the local clock sampling module. If the sample point is greater than the detection upper limit boundary of the corresponding timestamp or less than the detection lower limit boundary of the corresponding timestamp, a high-level hardware interrupt signal of the microcontroller's general input / output interface is immediately triggered, and the system log recording module is called to write the timestamp of the current abnormality and the extreme value of the deviation exceeding the boundary into the non-volatile flash memory, thus completing the abnormality judgment and alarm result output.

[0071] In some implementations, the offset is proportional to the number of missing confidence values ​​determined based on the average weight within the historical time window, the number of missing confidence values ​​is negatively correlated with the average weight, and the offset is proportional to the curvature norm of the nonparametric components within the window, including: The statistical iterative reweighted least squares method assigns convergence weights to each data point within a historical time window and calculates the average weight value. The number of missing confidence levels is determined based on the average weight value. The larger the average weight value, the smaller the number of missing confidence levels, and the smaller the average weight value, the larger the number of missing confidence levels. By performing differential processing on the nonparametric component sequence, curvature eigenvalues ​​are obtained to characterize the intensity of fluctuations in the nonparametric components. Based on the confidence level missing value and curvature feature value, the offset is determined according to the preset monotonic mapping relationship, segmented mapping relationship or lookup table mapping relationship; Among them, the larger the amount of missing confidence, the larger the offset; the larger the curvature feature value, the larger the offset.

[0072] When calculating the offset, a historical time window is retrieved, preferably containing state data within 120 consecutive valid clock tick samples. During the convergence of the iterative reweighted least squares calculation, the weights of the data points assigned to each time point within this historical time window are calculated, and all weight values ​​are summed (e.g., a total of 102), then divided by the total number of samples (120) to obtain the average weight value. =0.85. This average weight represents the confidence level of the fit of recent observations to the long-term drift model and the consistency of observations.

[0073] Based on average weight value Calculate the confidence level missing quantity q=1- In the example above, q = 1 - 0.85 = 0.15. This confidence level missing value represents the degree of missing recent observations relative to a fully credible state. The larger the average weight, the smaller the confidence level missing value; the smaller the average weight, the larger the confidence level missing value.

[0074] Simultaneously, nonparametric component sequences are input, and second-order central difference processing is performed on the nonparametric component sequences. Time scale normalization is then performed in conjunction with the sampling period Δt to construct the second-order difference results within the associated segment. And calculate the L2 norm of the second-order difference result. For example, the result is 0.35 ns / s². Wherein, The norm of the second-order difference result of the nonparametric components is used to characterize the random fluctuation energy characteristics of the clock in the nonparametric space.

[0075] After obtaining the confidence level missing value and curvature norm, according to the formula Establish candidate offsets, where... The term refers to the candidate offset, where K represents the multiplication gain adjustment constant that includes dimension transformation and gives the candidate offset clock bias dimensions, and q is the confidence level missing value determined based on the average weight within the historical time window. Let be the norm of the second-order difference result of the nonparametric components; then, according to the formula... Determine the actual offset, where Δd represents the actual offset. This indicates the preset minimum detection boundary width, used to prevent the detection boundary from degenerating to zero width under high confidence and low fluctuation conditions. The dimensions are At that time, the gain multiplication adjustment constant K has the dimension of the square of seconds, and in this example it is set to 340. This converts the output dimensions to nanoseconds, a measure of clock skew. Substituting the aforementioned parameters, we obtain the candidate offset calculation results. =340×0.15×0.35=17.85 nanoseconds; if the minimum detection boundary width is preset. If the offset is not greater than 17.85 nanoseconds, then the actual offset Δd is taken as 17.85 nanoseconds. If the time is greater than 17.85 nanoseconds, then the actual offset Δd is taken as... The actual offset is used for boundary calculations of the clock skew model and can be adjusted according to the noise environment and clock status of the current observation window, thereby reducing the false alarm rate of the clock.

[0076] The ablation experiment was conducted using a test platform built on a high and low temperature alternating damp heat test chamber and a BeiDou satellite navigation receiver terminal. During the experiment, the ambient temperature was controlled to undergo periodic and drastic changes between 10℃ and 40℃, and artificial obstruction was applied to the outdoor receiving antenna to cause the satellite signal carrier-to-noise ratio to fluctuate randomly within the range of 25 to 48 dB-Hz. The experiment ran continuously for 72 hours with a sampling period of 1 second, collecting a total of 259,200 pairs of original BeiDou time synchronization timestamps and local reference timestamps. This dataset was used to compare and verify the clock deviation sequence generation and prediction performance of different module combinations.

[0077] The experiment setup included a baseline group and three experimental groups for comparison and verification. The baseline group used a conventional, uncorrected clock skew prediction model, with a root mean square error (RMSE) of 45.3 ns and a maximum prediction deviation of 120 ns. Experimental group one removed the temperature and carrier-to-noise ratio (CNR) weighted correction modules from the complete scheme, resulting in an inability to suppress pseudorange noise caused by environmental and signal attenuation, with an RMS of 28.6 ns and a maximum prediction deviation of 85 ns. Experimental group two removed the bandwidth fine-tuning mechanism based on minimum mutual information and adopted fixed-bandwidth kernel regression from the complete scheme, causing correlation between long-term and short-term component information, with an RMS of 18.4 ns and a maximum prediction deviation of 52 ns. Experimental group three used the complete technical scheme presented in this paper, reducing the RMS to 6.2 ns and controlling the maximum prediction deviation to within 15 ns. A comparison of the prediction errors of each group is provided below. Figure 4 As shown, with the gradual utilization of each functional module, both the root mean square error and the maximum error continue to decrease, indicating that the proposed complete scheme has superior prediction performance.

[0078] The complete solution presented in this paper reduces the root mean square error by 86.3% compared to the traditional benchmark group, thus improving the system's timing reliability. Data from Experiment 1 shows that the dual adaptive weighting mechanism of ambient temperature and signal quality can mitigate timestamp jump interference caused by changes in physical conditions. Data from Experiment 2 demonstrates that by finding the optimal bandwidth based on the mutual information minimization criterion, the parametric components representing long-term trends and the non-parametric components representing random fluctuations are decoupled, suppressing over-smoothing during the fitting process and achieving prediction and modeling of clock deviation in the timing system.

[0079] In the description of this specification, "multiple" or "several" means at least two, such as two, three or more, unless otherwise expressly and specifically defined.

Claims

1. A method for detecting BeiDou clock synchronization anomalies in intelligent converged terminals for distribution areas, characterized in that, include: S1: Obtain the original timestamp of BeiDou time synchronization, satellite signal carrier-to-noise ratio and ambient temperature data of the intelligent fusion terminal, obtain the weighting coefficients and construct an adaptive weighted model to generate a clock deviation sequence; The adjustment constant is calculated and the clock deviation sequence is processed by iterative reweighted least squares to extract the parameter components, so as to obtain the drift model and drift coefficient representing the long-term trend of the clock. The drift coefficient is used to represent the inherent drift characteristics of the clock accumulated over time. S2, the parameter components are removed from the clock deviation sequence to obtain the residual sequence. Kernel regression is applied to the residual sequence to estimate the nonparametric components, which are used to represent the nonlinear random fluctuation characteristics of the clock. The initial bandwidth of kernel regression is generated based on the drift coefficient mapping, and fine-tuned to obtain the bandwidth by minimizing the mutual information between nonparametric and parametric components, thereby improving the estimation accuracy. S3, the detection boundary is set by adding or subtracting the offset based on the sum of the predicted values ​​of parametric components and the estimated values ​​of nonparametric components. The offset is proportional to the amount of missing confidence determined based on the average weight within the historical time window, and the amount of missing confidence is negatively correlated with the average weight. Furthermore, the offset is proportional to the curvature norm of the nonparametric components within the window, thereby achieving adjustment for different states. When a new clock deviation sample point exceeds the detection boundary, it is determined that the current BeiDou clock synchronization is abnormal and the result is output.

2. The method according to claim 1, characterized in that, The process of constructing an adaptive weighted model to generate a clock bias sequence includes: subtracting the original timestamp of BeiDou time service from the local reference timestamp to obtain the basic clock bias sequence; and obtaining the weight coefficient for each time point, wherein the weight coefficient has a positive nonlinear relationship with the amplitude of the satellite signal carrier-to-noise ratio at the current time. The weighting coefficients are used as confidence weights for the corresponding data points in the basic clock deviation sequence. Weighted smoothing, weighted robust estimation, or weighted filtering are then applied to the basic clock deviation sequence to obtain the signal quality corrected basic clock deviation sequence. The larger the weighting coefficient, the greater the contribution of the corresponding data point to the correction result; the smaller the weighting coefficient, the smaller the contribution of the corresponding data point to the correction result. The adjusted base clock skew sequence is added element-wise to the time delay increment sequence to generate the clock skew sequence.

3. The method according to claim 2, characterized in that, The acquisition of weight coefficients includes: Obtain the current satellite signal carrier-to-noise ratio amplitude; The current signal quality characterization value is determined based on the relative relationship between the satellite signal carrier-to-noise ratio amplitude and the preset reference carrier-to-noise ratio amplitude. Generate the corresponding weighting coefficients according to the preset mapping relationship between carrier-to-noise ratio and weighting coefficients; The mapping relationship can be a piecewise mapping relationship, a lookup table mapping relationship, or a monotonically increasing function mapping relationship. The weight coefficient increases with the increase of the satellite signal carrier-to-noise ratio amplitude. The weight coefficient is used to characterize the reliability of the basic clock deviation data point at the corresponding time and to determine the contribution of the corresponding data point in the weighted processing.

4. The method according to claim 1, characterized in that, The clock deviation sequence is processed using an iterative reweighted least squares method to extract parameter components, resulting in a drift model and drift coefficients representing the long-term trend of the clock, including: A long-term trend fitting model containing a constant term, a primary trend term, and a secondary trend term is constructed based on the clock deviation sequence. Initial estimates are made for the various model parameters of the long-term trend fitting model based on the clock deviation sequence; The parameters of the converged quadratic trend term are extracted as the drift coefficients; The corresponding trend sequence is generated based on the converged long-term trend fitting model, and the trend sequence is used as a parameter component representing the long-term trend of the clock signal.

5. The method according to claim 4, characterized in that, After initially estimating the various model parameters of the long-term trend fitting model, the process also includes: The weights of each data point in the current round are determined based on a bisquared weighting function or other stable weighting functions, and the model parameters are iteratively updated using a weighted fitting method until the change in model parameters between two adjacent iterations is lower than a preset convergence threshold.

6. The method according to claim 5, characterized in that, The calculation of the adjustment constant includes: For the clock skew sequence, subtract the data points of two adjacent moments in chronological order to generate a first-order difference sequence; Sort all values ​​in the first-order difference sequence in ascending order; The interquartile range is calculated by subtracting the value corresponding to the 25th percentile from the value corresponding to the 75th percentile in the sorted first-order difference sequence. Divide the interquartile range by 1.349 and multiply by 4.685 to obtain the adjustment constant of the current bisquared weighting function.

7. The method according to claim 1, characterized in that, The initial bandwidth of the kernel regression is generated based on the drift coefficient mapping, including: Extract the drift coefficient and calculate its absolute value; Calculate the absolute value of the drift coefficient and add it to a preset small positive bias to obtain the drift correction amount; divide the preset first constant by the drift correction amount to obtain the dimensionless inverse factor; wherein, the first constant and the small positive bias have the same dimensions as the drift coefficient. Multiply the preset scaling factor constant with time dimension by the dimensionless inverse factor, and set the resulting value as the initial bandwidth required for kernel regression.

8. The method according to claim 1, characterized in that, The method of fine-tuning the bandwidth based on minimizing the mutual information between nonparametric and parametric components includes: Based on the initial bandwidth and the preset search step size, multiple candidate bandwidths are generated; Kernel regression estimation is performed on the residual sequence using each candidate bandwidth to obtain the candidate nonparametric component sequence corresponding to each candidate bandwidth; Based on the preset mutual information evaluation model, calculate the mutual information evaluation value between each candidate nonparametric component sequence and the parametric component sequence; The candidate bandwidth with the smallest mutual information evaluation value is selected as the fine-tuned bandwidth. Based on the fine-tuned bandwidth, kernel regression estimation is performed on the residual sequence to obtain nonparametric components.

9. The method according to claim 1, characterized in that, The offset is proportional to the number of missing confidence levels determined based on the average weight within the historical time window. The number of missing confidence levels is negatively correlated with the average weight. Furthermore, the offset is proportional to the curvature norm of the nonparametric component within the window. This includes: assigning convergence weights to each data point within the historical time window using the statistical iterative reweighted least squares method, and calculating the average weight value; determining the number of missing confidence levels based on the average weight value, where a larger average weight value results in a smaller number of missing confidence levels, and vice versa; obtaining curvature feature values ​​to characterize the intensity of fluctuations in the nonparametric component; and determining the offset based on the number of missing confidence levels and the curvature feature values ​​according to a preset monotonic mapping relationship, piecewise mapping relationship, or lookup table mapping relationship. The larger the number of missing confidence levels, the larger the offset; and the larger the curvature feature value, the larger the offset.

10. The method according to claim 9, characterized in that, The step of obtaining curvature feature values ​​to characterize the intensity of fluctuations in nonparametric components includes: performing differential processing on the nonparametric component sequence to obtain the curvature feature values.

Images