Hydrogen atomic clock frequency drift estimation method and device based on kalman filtering

Abstract

This invention relates to the field of hydrogen atomic clock technology, and more specifically to a method and apparatus for estimating the frequency drift of a hydrogen atomic clock based on Kalman filtering. The method includes: receiving raw frequency data of a hydrogen atomic clock sampled at a preset sampling frequency in real time; performing Kalman filtering processing at a preset processing frequency to estimate the state, wherein the state includes the frequency drift and drift rate of the hydrogen atomic clock, and the preset processing frequency is less than the preset sampling frequency.

Description

Technical Field

[0001] This invention relates to the field of hydrogen atomic clock technology, and more specifically to a method and apparatus for estimating the frequency drift of a hydrogen atomic clock based on Kalman filtering. Background Technology

[0002] Hydrogen atomic clocks possess excellent short-time frequency stability, making them a core reference source for navigation, telemetry, and control, national timekeeping, and precision time and frequency systems. As the demands for time and frequency accuracy continue to rise in high-tech fields, the long-term frequency accuracy and stability of hydrogen atomic clocks have become key indicators limiting system performance.

[0003] Due to its own working mechanism and the influence of the external environment, the hydrogen atomic clock inevitably experiences frequency drift during operation. Factors such as the internal resonant cavity traction effect, aging of the cavity wall coating, atomic interaction losses, as well as external temperature fluctuations, magnetic field interference, and air pressure changes, all cause its output frequency to slowly shift over time. Long-term frequency drift will significantly deteriorate the long-term stability of the atomic clock, reduce the accuracy of time synchronization and ranging, and fail to meet the requirements for maintaining a long-period, high-precision spatiotemporal reference.

[0004] Currently, most existing hydrogen atomic clock frequency drift estimation techniques rely on traditional linear fitting, which cannot accurately characterize the long-term frequency drift evolution of hydrogen atomic clocks. The accuracy of frequency drift estimation and compensation is limited, making it difficult to fully realize the ultra-high precision timing potential of hydrogen atomic clocks. Summary of the Invention

[0005] The purpose of this invention is to provide a method and apparatus for estimating the frequency drift of a hydrogen atomic clock based on Kalman filtering, so as to accurately estimate the frequency drift of a hydrogen atomic clock.

[0006] To achieve the above objectives, the present invention provides a method for estimating the frequency drift of a hydrogen atomic clock based on Kalman filtering, comprising: Real-time reception of raw frequency data of the hydrogen atomic clock obtained by sampling at a preset sampling frequency; Kalman filtering is performed at a preset processing frequency to estimate the state, which includes the frequency drift and drift rate of the hydrogen atomic clock. The preset processing frequency is less than the preset sampling frequency. Specifically, performing Kalman filtering at the preset processing frequency to estimate the state includes: For each processing moment, the mean of all raw frequency data of hydrogen atomic clocks received between the current processing moment and the previous processing moment is calculated and used as the observation value at the current processing moment; If the sequence number of the current processing time is less than or equal to the preset preheating window length, preheating is performed to obtain the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance. If the sequence number of the current processing time is greater than the preheating window length, then Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as the iteration benchmark to obtain the posterior state estimate of the current processing time.

[0007] Optionally, if the sequence number of the current processing time is less than or equal to the preset preheating window length, preheating is performed to obtain the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance, specifically including: If the current processing time number is less than the preheating window length, then accumulate the observations; If the sequence number of the current processing time is equal to the length of the preheating window, then the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance are obtained based on the observations accumulated from the current time and all previous processing times; and the process noise covariance matrix is ​​determined based on the initial measurement noise variance.

[0008] Optionally, Kalman filtering iterations are performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as iterative benchmarks to obtain the posterior state estimate at the current processing time, specifically including: Based on the initial measurement noise variance, determine the measurement noise variance at the current processing moment; The prior state estimate for the current processing time is determined based on the posterior state estimate of the previous processing time, and the prior covariance matrix for the current processing time is determined based on the posterior covariance matrix and the process noise covariance matrix of the previous processing time. The innovation at the current processing time is determined based on the observed values ​​and prior state estimates at the current processing time; the covariance of the innovation at the current processing time is determined based on the prior covariance matrix and measurement noise variance at the current processing time; the innovation at the current processing time is normalized based on the covariance of the innovation at the current processing time to obtain the normalized innovation at the current processing time. Determine if there is a sudden change in the normalized information at the current processing time. If so, reset the prior covariance matrix at the current processing time to the preset value. The Kalman gain at the current processing time is determined based on the prior covariance matrix at the current processing time and the covariance of the innovation at the current processing time. The posterior state estimate at the current processing time is determined based on the prior state estimate and Kalman gain at the current processing time; the posterior covariance matrix at the current processing time is determined based on the prior covariance matrix and Kalman gain at the current processing time.

[0009] Optionally, based on the initial measurement noise variance, the measurement noise variance at the current processing moment is determined, specifically including: The time between the current processing moment and the previous processing moment is taken as the current processing cycle, and the actual amount of raw frequency data of the hydrogen atomic clock received within the current processing cycle is obtained. The ratio of the preset sampling frequency to the preset processing frequency is obtained as the theoretical quantity of raw frequency data of the hydrogen atomic clock received in each processing cycle. The ratio between the actual number and the theoretical number of hydrogen atomic clock frequency data received in the current processing cycle is obtained and used as the weight of the current processing moment. The initial measurement noise variance is multiplied by the weight at the current processing time to obtain the measurement noise variance at the current time.

[0010] Optionally, determining whether there is a sudden change in the normalized information at the current processing time includes: The normalized news at the current time is input into the CUSUM detector to determine whether there is a mutation in the normalized news at the current processing time.

[0011] Another aspect of the present invention provides a hydrogen atomic clock frequency drift estimation device based on Kalman filtering, comprising: The receiving module is used to receive the raw frequency data of the hydrogen atomic clock obtained by sampling at a preset sampling frequency in real time. The processing module is used to perform Kalman filtering at a preset processing frequency to estimate the state, which includes the frequency drift and drift rate of the hydrogen atomic clock, wherein the preset processing frequency is less than the preset sampling frequency; wherein, performing Kalman filtering at the preset processing frequency to estimate the state specifically includes: For each processing moment, the mean of all raw frequency data of hydrogen atomic clocks received between the current processing moment and the previous processing moment is calculated and used as the observation value at the current processing moment; If the sequence number of the current processing time is less than or equal to the preset preheating window length, preheating is performed to obtain the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance. If the sequence number of the current processing time is greater than the preheating window length, then Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as the iteration benchmark to obtain the posterior state estimate of the current processing time.

[0012] Optionally, if the sequence number of the current processing time is less than or equal to the preset preheating window length, preheating is performed to obtain the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance, specifically including: If the current processing time number is less than the preheating window length, then accumulate the observations; If the sequence number of the current processing time is equal to the length of the preheating window, then the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance are obtained based on the observations accumulated from the current time and all previous processing times; and the process noise covariance matrix is ​​determined based on the initial measurement noise variance.

[0013] Optionally, Kalman filtering iterations are performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as iterative benchmarks to obtain the posterior state estimate at the current processing time, specifically including: Based on the initial measurement noise variance, determine the measurement noise variance at the current processing moment; The prior state estimate for the current processing time is determined based on the posterior state estimate of the previous processing time, and the prior covariance matrix for the current processing time is determined based on the posterior covariance matrix and the process noise covariance matrix of the previous processing time. The innovation at the current processing time is determined based on the observed values ​​and prior state estimates at the current processing time; the covariance of the innovation at the current processing time is determined based on the prior covariance matrix and measurement noise variance at the current processing time; the innovation at the current processing time is normalized based on the covariance of the innovation at the current processing time to obtain the normalized innovation at the current processing time. Determine if there is a sudden change in the normalized information at the current processing time. If so, reset the prior covariance matrix at the current processing time to the preset value. The Kalman gain at the current processing time is determined based on the prior covariance matrix at the current processing time and the covariance of the innovation at the current processing time. The posterior state estimate at the current processing time is determined based on the prior state estimate and Kalman gain at the current processing time; the posterior covariance matrix at the current processing time is determined based on the prior covariance matrix and Kalman gain at the current processing time.

[0014] Optionally, based on the initial measurement noise variance, the measurement noise variance at the current processing moment is determined, specifically including: The time between the current processing moment and the previous processing moment is taken as the current processing cycle, and the actual amount of raw frequency data of the hydrogen atomic clock received within the current processing cycle is obtained. The ratio of the preset sampling frequency to the preset processing frequency is obtained as the theoretical quantity of raw frequency data of the hydrogen atomic clock received in each processing cycle. The ratio between the actual number and the theoretical number of hydrogen atomic clock frequency data received in the current processing cycle is obtained and used as the weight of the current processing moment. The initial measurement noise variance is multiplied by the weight at the current processing time to obtain the measurement noise variance at the current time.

[0015] Optionally, determining whether there is a sudden change in the normalized information at the current processing time includes: The normalized news at the current time is input into the CUSUM detector to determine whether there is a mutation in the normalized news at the current processing time. Attached Figure Description

[0016] Figure 1 This is a flowchart of a hydrogen atomic clock frequency drift estimation method based on Kalman filtering according to an embodiment of the present invention; Figure 2 This is a structural block diagram of a hydrogen atomic clock frequency drift estimation device based on Kalman filtering according to an embodiment of the present invention. Detailed Implementation

[0017] The preferred embodiments of the present invention are given below with reference to the accompanying drawings and described in detail.

[0018] like Figure 1 As shown, this embodiment of the invention provides a method for estimating the frequency drift of a hydrogen atomic clock based on Kalman filtering, which includes the following steps S100-S200: S100: Receives raw frequency data of the hydrogen atomic clock in real time, sampled at a preset sampling frequency.

[0019] S200: Perform Kalman filtering at a preset processing frequency to estimate the state, including the frequency drift and drift rate of the hydrogen atomic clock.

[0020] No. k The state vector at each processing moment can be denoted as: X k : (1) in , Let represent the relative frequency deviation and drift rate at time t, respectively. The state transition equation and observation equation of the Kalman filter can be expressed as: (2) (3) in, , , For the first k The observation value at each moment, The processing period is defined by F, the state transition matrix, and H, the observation matrix. This represents the process noise vector, reflecting random disturbances in the actual physical process. Observation noise is used to represent random errors in a single observation.

[0021] Kalman filtering requires determining the process noise vector. covariance matrix With measurement noise variance .in, Used to characterize the model's own uncertainty regarding the fact that the drift rate may slowly wander. Used to characterize the magnitude of random noise in a single window of observation.

[0022] The preset sampling frequency can be N times the preset processing frequency, where N is a positive integer greater than 1. Thus, the processing period (the interval between any two adjacent processing moments) is N times the sampling period (the interval between any two adjacent sampling points). In other words, theoretically, within each processing period, the original frequency data of hydrogen atoms from N sampling points can be received.

[0023] In some embodiments, step S200 specifically includes the following steps S210-S230: S210: For each processing moment, calculate the average of all raw frequency data of the hydrogen atomic clock received between the current processing moment and the previous processing moment, and use this as the observation value for the current processing moment. k The formula for calculating the observations at each processing time is as follows: (4) in, For the first k The processing cycle (i.e., the ) k The processing time and the first k- 1. Processing time between moments. For the first i The raw frequency data of the hydrogen atomic clock at each sampling point For the first k Theoretically, the number of raw frequency data from all hydrogen atomic clocks within a processing cycle... It should be N, but in reality, due to various reasons such as data acquisition system restart, network communication interruption, storage device malfunction, or data acquisition being suspended during instrument maintenance, It could also be less than N, or even 0.

[0024] S220: If the sequence number of the current processing time is less than or equal to the preset preheating window length, preheating is performed to obtain the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance.

[0025] In some embodiments, the preheating window length can be set according to actual needs, for example, 30-50. The preheating stage mainly plays two roles: 1) maintaining the basic stability of the initial state and avoiding drastic fluctuations in estimation due to insufficient initialization at the beginning; 2) adaptively estimating the measurement noise variance based on the preheating data so that the subsequent filter parameters are automatically matched with the current data volume.

[0026] Step S220 specifically includes: S222: If the sequence number of the current processing time is less than the preheating window length, then accumulate the observations; S224: If the sequence number of the current processing time is equal to the preheating window length, then initialize based on the observations of the current time and all previous processing times to obtain the initial posterior state estimate, the initial posterior covariance matrix, the initial measurement noise variance, and the process noise covariance matrix.

[0027] Step S224 specifically includes: Construct a sequence of differences between adjacent observations based on the observations at the current time and all previous processing times, calculate the variance of this sequence, and then take half of it to obtain the initial measurement noise variance. Determine the process noise covariance matrix based on the initial measurement noise variance; Obtain the average of the observations at the current time and all previous processing times as the initial posterior state estimate; The initial posterior covariance matrix is ​​based on the degree of dispersion of the observations at the current time and all previous processing times.

[0028] The initial measurement noise variance is calculated as follows: (5) The process noise covariance matrix can be expressed in the following form: (6) in, This represents the spectral density driven by the continuous-time drift rate, used to control the scale of the noise matrix throughout the process. To avoid directly manipulating the spectral density parameter... Parameter tuning was performed, and a smoothing ratio was further introduced. Defined as: (7) in Represents the process noise matrix The element in the top left corner. top left element The uncertainty of the corresponding frequency shift state process, therefore the smoothing ratio In reality, it describes the relative magnitude between this uncertainty and the variance of the single observation noise. The smaller the value, the smoother and more conservative the filter; conversely, the larger the value, the more sensitive and easier the filter is to adapt to trends. Therefore, it can be written as: (8) For example, the preheating window length is 30. It is 0.001.

[0029] S230: If the sequence number of the current processing time is greater than the preheating window length, then Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as the iteration benchmark to obtain the posterior state estimate of the current processing time.

[0030] In some embodiments, Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as iterative benchmarks to obtain the posterior state estimate at the current processing time. This specifically includes the following steps: S231: Based on the initial measurement noise variance, determine the measurement noise variance at the current processing moment.

[0031] No. k The formula for calculating the measurement noise variance at each processing time is as follows: (9) in N represents the theoretical amount of raw frequency data from the hydrogen atomic clock received in one processing cycle. n k For the first k The actual number of raw frequency data received by the hydrogen atomic clock in a processing cycle. The fewer the number of raw frequency data in a processing cycle, the larger the mean variance, and the lower the observation reliability. As can be seen from formula (9), the corresponding measurement noise variance will also increase. In this way, it is not necessary to discard insufficient processing times in order to maintain the regular time grid, and the difference in observation quality at different processing times can be reflected in the model.

[0032] S232: Determine the prior state estimate of the current processing time based on the posterior state estimate of the previous processing time, and determine the prior covariance matrix of the current processing time based on the posterior covariance matrix and the process noise covariance matrix of the previous processing time.

[0033] The calculation formula is as follows: (10) (11) in, This indicates that the time has not yet been obtained. k Real observation At time, for a specific momentk Prior prediction of the state; This is the prior covariance matrix at the current processing time, describing the uncertainty of the prior prediction.

[0034] S233: Determine the innovation at the current processing time based on the observed value and the prior state estimate at the current processing time; determine the covariance of the innovation at the current processing time based on the prior covariance matrix and the measurement noise variance at the current processing time; normalize the innovation at the current processing time based on the covariance of the innovation at the current processing time to obtain the normalized innovation at the current processing time.

[0035] New This represents the difference between the current observed value and the model's predicted value. If the model perfectly matches the actual system, It should be approximated as a zero-mean random variable; if systematic bias occurs, It will continue to skew in a certain direction, which is the source of the signal for subsequent anomaly detection.

[0036] No. k The formula for calculating the new information at each processing time is as follows: (12) New The covariance is: (13) It combines current prior uncertainty and measurement noise The two parts constitute a complete quantification of the expected fluctuation range of current interest rates. Further normalization of interest rates is then performed: (14) The introduction of normalization has three specific effects: (1) it eliminates the difference in absolute dimensions; (2) when the model matches the actual system well, it theoretically approximately follows a zero-mean unit variance distribution; and (3) it is more suitable as a unified input to the sequential detector.

[0037] S234: Determine if there is a sudden change in the normalized information at the current processing time. If so, reset the prior covariance matrix at the current processing time to the preset value.

[0038] In some embodiments, a CUSUM detector may be used for mutation detection. The CUSUM detector is as follows: (15) (16) in This is a tolerance parameter, representing the maximum allowable offset per step. One standard deviation is not included in the cumulative total, which is used to absorb random fluctuations in measurement noise and prevent a single large innovation from being misjudged as an anomaly; It is the detection threshold, indicating the cumulative threshold exceeding [a certain value]. An alarm is triggered after one standard deviation, corresponding to approximately level.

[0039] when or At that time, it was determined that in the first Structural abrupt changes occur. These structural abrupt changes include, but are not limited to: (1) significant changes in drift rate; (2) frequency jumps; (3) trend reversals; and (4) persistent model mismatch with current observations.

[0040] Specifically, the selected parameters are: , Initialization When a sudden change is detected, the model does not directly modify the state value, but instead resets the prior covariance matrix at the current processing time. (17) The core logic of this design is: state estimation. The frequency offset and drift rate values ​​stored in the old covariance matrix may still be close to the true values ​​under the new trend, so they should not be forced to zero; however, we do no longer trust the uncertainty range described by the old covariance, therefore we reset the covariance matrix. This indicates that the current state and position are highly uncertain.

[0041] S235: Determine the Kalman gain at the current processing time based on the prior covariance matrix at the current processing time and the covariance of the innovation at the current processing time.

[0042] No. k Kalman gain at each moment The calculation formula is as follows: (18) when When the size is large, the filter absorbs the current observation faster; when When the size is small, the filter relies more on prior predictions.

[0043] S236: Determine the posterior state estimate at the current processing time based on the prior state estimate and Kalman gain at the current processing time; determine the posterior covariance matrix at the current processing time based on the prior covariance matrix and Kalman gain at the current processing time.

[0044] Step S236 is the standard Kalman update process, and its specific formula is as follows: (19) (20) in, This indicates that it has been obtained. k The posterior state update obtained after the observation at time step. This represents the posterior covariance at the current processing time.

[0045] The posterior state estimate is the optimal estimate at the current processing time.

[0046] The frequency drift estimation method for hydrogen atomic clocks based on Kalman filtering in this invention uses the mean of the original frequency data within the processing period as the observed value, and estimates the frequency drift and drift rate based on the Kalman filtering algorithm, achieving high accuracy; estimation is performed during the preheating stage. R Then use smoothing ratio generate This avoids the problem of repeated manual parameter tuning; in addition, the dynamic scaling of measurement noise in incomplete windows ensures that the scheme can maintain relatively stable performance under different data quality conditions; through the combination mechanism of normalized innovation CUSUM detector and covariance reset, it maintains smoothness during the stationary period and quickly relearns after abrupt changes, thus simultaneously taking into account robustness and response speed; this invention can estimate the frequency drift of hydrogen atomic clock in real time and output a state estimate based on past processing times at the current processing time. The state estimate can be directly used for real-time monitoring and online control front end.

[0047] like Figure 2 As shown, this embodiment of the invention also provides a hydrogen atomic clock frequency drift estimation device based on Kalman filtering, which includes a receiving module 10 and a processing module 20. The receiving module 10 is used to receive the original frequency data of the hydrogen atomic clock obtained by sampling at a preset sampling frequency in real time. The processing module 20 is used to perform Kalman filtering processing at a preset processing frequency to estimate the state, the state including the frequency drift and drift rate of the hydrogen atomic clock, and the preset processing frequency is less than the preset sampling frequency.

[0048] Specifically, Kalman filtering is performed at a preset processing frequency to estimate the state, including: For each processing moment, the mean of all raw frequency data of hydrogen atomic clocks received between the current processing moment and the previous processing moment is calculated and used as the observation value at the current processing moment; If the sequence number of the current processing time is less than or equal to the preset preheating window length, preheating is performed to obtain the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance. If the sequence number of the current processing time is greater than the preheating window length, then Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, and the initial measurement noise variance as the iteration benchmark to obtain the posterior state estimate of the current processing time.

[0049] The receiving module 10 and the processing module 20 are the functional modules corresponding to steps S100 and S200 in the method embodiment, respectively. Their specific implementation methods can be referred to the description in the method embodiment, and will not be repeated here.

[0050] Another embodiment of the present invention provides a readable storage medium having a computer program stored thereon, which, when executed in a computer, causes the computer to perform the steps of the hydrogen atomic clock frequency drift estimation method based on Kalman filtering in the above embodiments of the present invention.

[0051] Another embodiment of the present invention provides an electronic device, which includes a memory and a processor. The memory stores executable code. When the processor executes the executable code, it performs the steps of the hydrogen atom clock frequency drift estimation method based on Kalman filtering in the above embodiments of the present invention.

[0052] The systems, devices, modules, or units described in the above embodiments can be implemented by computer chips or entities, or by products with certain functions. A typical implementation device is a computer. Specifically, a computer can be, for example, a personal computer, laptop computer, cellular phone, camera phone, smartphone, personal digital assistant, media player, navigation device, email device, game console, tablet computer, wearable device, or any combination of these devices.

[0053] For ease of description, the above apparatus is described in terms of its functions, divided into various units. Of course, in implementing this invention, the functions of each unit can be implemented in one or more software and / or hardware components.

[0054] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0055] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0056] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0057] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0058] In a typical configuration, an electronic device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0059] Memory may include non-persistent storage in computer-readable media, such as random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM. Memory is an example of computer-readable media.

[0060] Computer-readable media includes both permanent and non-permanent, removable and non-removable media that can store information by any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, disk storage or other magnetic storage devices, or any other non-transferable medium that can be used to store information accessible by electronic devices. As defined herein, computer-readable media does not include transient computer-readable media, such as modulated data signals and carrier waves.

[0061] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0062] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0063] This invention can be described in the general context of computer-executable instructions, such as program modules, that are executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, etc., that perform a specific task or implement a specific abstract data type. This invention can also be practiced in distributed computing environments where tasks are performed by remote processing devices connected via a communication network. In distributed computing environments, program modules can reside in local and remote computer storage media, including storage devices.

[0064] The various embodiments in this invention are described in a progressive manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions in the method embodiments.

[0065] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of the invention. Various variations can be made to the above embodiments of the present invention. That is, all simple and equivalent changes and modifications made based on the claims and description of this invention fall within the protection scope of the claims of this patent. All aspects not described in detail in this invention are conventional technical content.

Claims

1. A method for estimating the frequency drift of a hydrogen atomic clock based on Kalman filtering, characterized in that, include: Real-time reception of raw frequency data of the hydrogen atomic clock obtained by sampling at a preset sampling frequency; Kalman filtering is performed at a preset processing frequency to estimate the state, which includes the frequency drift and drift rate of the hydrogen atomic clock. The preset processing frequency is less than the preset sampling frequency. Specifically, performing Kalman filtering at the preset processing frequency to estimate the state includes: For each processing moment, the mean of all raw frequency data of hydrogen atomic clocks received between the current processing moment and the previous processing moment is calculated and used as the observation value at the current processing moment; If the current processing time number is less than or equal to the preset warm-up window length, the observations from each processing time are accumulated; and when the current processing time number equals the warm-up window length, noise parameter self-calibration and filter initialization are performed. The noise parameter self-calibration includes: constructing a difference sequence of adjacent observations based on the accumulated observations, using half of the variance of the difference sequence as the initial measurement noise variance; and using a preset dimensionless smoothing ratio and the product of the initial measurement noise variance to back-calculate the process noise covariance matrix, where the smoothing ratio is the ratio of the element corresponding to the frequency offset state in the process noise covariance matrix to the initial measurement noise variance. The filter initialization includes: determining the initial posterior state estimate and the initial posterior covariance matrix based on the accumulated observations; The initial measurement noise variance and the process noise covariance matrix serve as noise model parameters and are continuously used in the Kalman filtering calculation at each processing time after preheating. If the sequence number of the current processing time is greater than the preheating window length, then Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, the initial measurement noise variance, and the process noise covariance matrix as the iteration benchmark to obtain the posterior state estimate of the current processing time.

2. The hydrogen atom clock frequency drift estimation method based on Kalman filtering according to claim 1, characterized in that, The Kalman filter iteration includes: Based on the initial measurement noise variance and the ratio of the actual to the theoretical number of raw frequency data of the hydrogen atomic clock received between the current processing time and the previous processing time, the measurement noise variance at the current processing time is determined; when the actual number is zero, the Kalman update is not performed at the current processing time, and only the prior estimate is used as the output. Based on the posterior state estimate, the posterior covariance matrix of the previous processing time, and the process noise covariance matrix, the prior state estimate and prior covariance matrix of the current processing time are determined. The innovation is determined based on the observations and prior state estimates at the current processing time. The covariance of the innovation is determined based on the prior covariance matrix at the current processing time and the measurement noise variance. The innovation is then normalized using the covariance of the innovation to obtain the normalized innovation. The normalized information is input into the CUSUM detector to determine whether a structural change has occurred at the current processing time. If a change has occurred, the state estimate at the current processing time is not modified, but the prior covariance matrix at the current processing time is reset to a preset value. The Kalman gain is determined based on the prior covariance matrix at the current processing time and the covariance of the new information, and the posterior state estimate and posterior covariance matrix at the current processing time are determined accordingly.

3. The hydrogen atom clock frequency drift estimation method based on Kalman filtering according to claim 1, characterized in that, The process noise covariance matrix is ​​constructed from a spectral density driven by a random continuous-time drift rate, and the spectral density is determined by inversely calculating the smoothing ratio and the initial measurement noise variance; the smaller the smoothing ratio, the smoother the filter; the larger the smoothing ratio, the more sensitive the filter.

4. The hydrogen atomic clock frequency drift estimation method based on Kalman filtering according to claim 2, characterized in that, Determine the measurement noise variance at the current processing time, specifically including: The time between the current processing moment and the previous processing moment is taken as the current processing cycle, and the actual amount of raw frequency data of the hydrogen atomic clock received within the current processing cycle is obtained. The ratio of the preset sampling frequency to the preset processing frequency is obtained as the theoretical quantity of raw frequency data of the hydrogen atomic clock received in each processing cycle. The ratio between the actual number and the theoretical number of hydrogen atomic clock frequency data received in the current processing cycle is obtained and used as the weight of the current processing moment. The initial measurement noise variance is multiplied by the weight at the current processing time to obtain the measurement noise variance at the current time.

5. The hydrogen atom clock frequency drift estimation method based on Kalman filtering according to claim 2, characterized in that, The CUSUM detector includes a positive accumulator and a negative accumulator; after deducting the tolerance parameter from the normalized information, it is accumulated in both positive and negative directions to obtain the positive accumulation amount and the negative accumulation amount; when the positive accumulation amount or the negative accumulation amount exceeds a preset threshold, it is determined that a structural change has occurred at the current processing time; The length of the preheating window is 30 to 50 mm.

6. A hydrogen atomic clock frequency drift estimation device based on Kalman filtering, characterized in that, include: The receiving module is used to receive the raw frequency data of the hydrogen atomic clock obtained by sampling at a preset sampling frequency in real time. The processing module is used to perform Kalman filtering at a preset processing frequency to estimate the state, which includes the frequency drift and drift rate of the hydrogen atomic clock, wherein the preset processing frequency is less than the preset sampling frequency; wherein, performing Kalman filtering at the preset processing frequency to estimate the state specifically includes: For each processing moment, the mean of all raw frequency data of hydrogen atomic clocks received between the current processing moment and the previous processing moment is calculated and used as the observation value at the current processing moment; If the current processing time number is less than or equal to the preset warm-up window length, the observations from each processing time are accumulated; and when the current processing time number equals the warm-up window length, noise parameter self-calibration and filter initialization are performed. The noise parameter self-calibration includes: constructing a difference sequence of adjacent observations based on the accumulated observations, using half of the variance of the difference sequence as the initial measurement noise variance; and using a preset dimensionless smoothing ratio and the product of the initial measurement noise variance to back-calculate the process noise covariance matrix, where the smoothing ratio is the ratio of the element corresponding to the frequency offset state in the process noise covariance matrix to the initial measurement noise variance. The filter initialization includes: determining the initial posterior state estimate and the initial posterior covariance matrix based on the accumulated observations; The initial measurement noise variance and the process noise covariance matrix serve as noise model parameters and are continuously used in the Kalman filtering calculation at each processing time after preheating. If the sequence number of the current processing time is greater than the preheating window length, then Kalman filtering iteration is performed using the initial posterior state estimate, the initial posterior covariance matrix, the initial measurement noise variance, and the process noise covariance matrix as the iteration benchmark to obtain the posterior state estimate of the current processing time.

7. The hydrogen atomic clock frequency drift estimation device based on Kalman filtering according to claim 6, characterized in that, The Kalman filter iteration includes: Based on the initial measurement noise variance and the ratio of the actual to the theoretical number of raw frequency data of the hydrogen atomic clock received between the current processing time and the previous processing time, the measurement noise variance at the current processing time is determined; when the actual number is zero, the Kalman update is not performed at the current processing time, and only the prior estimate is used as the output. Based on the posterior state estimate, the posterior covariance matrix of the previous processing time, and the process noise covariance matrix, the prior state estimate and prior covariance matrix of the current processing time are determined. The innovation is determined based on the observations and prior state estimates at the current processing time. The covariance of the innovation is determined based on the prior covariance matrix at the current processing time and the measurement noise variance. The innovation is then normalized using the covariance of the innovation to obtain the normalized innovation. The normalized information is input into the CUSUM detector to determine whether a structural change has occurred at the current processing time. If a change has occurred, the state estimate at the current processing time is not modified, but the prior covariance matrix at the current processing time is reset to a preset value. The Kalman gain is determined based on the prior covariance matrix at the current processing time and the covariance of the new information, and the posterior state estimate and posterior covariance matrix at the current processing time are determined accordingly.

8. The hydrogen atomic clock frequency drift estimation device based on Kalman filtering according to claim 6, characterized in that, The process noise covariance matrix is ​​constructed from a spectral density driven by a random continuous-time drift rate, and the spectral density is determined by inversely calculating the smoothing ratio and the initial measurement noise variance; the smaller the smoothing ratio, the smoother the filter; the larger the smoothing ratio, the more sensitive the filter.

9. The hydrogen atomic clock frequency drift estimation device based on Kalman filtering according to claim 7, characterized in that, Determine the measurement noise variance at the current processing time, specifically including: The time between the current processing moment and the previous processing moment is taken as the current processing cycle, and the actual amount of raw frequency data of the hydrogen atomic clock received within the current processing cycle is obtained. The ratio of the preset sampling frequency to the preset processing frequency is obtained as the theoretical quantity of raw frequency data of the hydrogen atomic clock received in each processing cycle. The ratio between the actual number and the theoretical number of hydrogen atomic clock frequency data received in the current processing cycle is obtained and used as the weight of the current processing moment. The initial measurement noise variance is multiplied by the weight at the current processing time to obtain the measurement noise variance at the current time.

10. The hydrogen atom clock frequency drift estimation device based on Kalman filtering according to claim 7, characterized in that, The CUSUM detector includes a positive accumulator and a negative accumulator; after deducting the tolerance parameter from the normalized information, it is accumulated in both positive and negative directions to obtain the positive accumulation amount and the negative accumulation amount; when the positive accumulation amount or the negative accumulation amount exceeds a preset threshold, it is determined that a structural change has occurred at the current processing time; The length of the preheating window is 30 to 50 mm.

Images