Method for establishing a composite clock based on an atomic clock and a millisecond pulsar

Abstract

A kind of composite clock establishment method based on atomic clock and millisecond pulsar, reads millisecond pulsar timing observation data, calculates the timing residual of millisecond pulsar;The timing residual of millisecond pulsar is sorted according to the order of reduced julian day;Millisecond pulsar that meets the establishment of millisecond pulsar time is screened out;Whether the timing residual of screened millisecond pulsar is multi-band timing residual is judged;Multi-band timing residual is statistically processed to obtain average timing residual;Average timing residual and timing red noise identification and processing of timing residual that is not multi-band are carried out, to obtain millisecond pulsar time time sequence;Read atomic clock measurement comparison data, carry out frequency and phase jump identification and rough error detection preprocessing, to obtain atomic time time sequence;Millisecond pulsar time and atomic clock time are combined using Vondrak-Cepek combination method to obtain composite clock Q.The composite clock established by the present application has both millisecond pulsar time long-term stability and accuracy, and also has atomic clock time high resolution and medium and short-term stability.

Description

Technical Field

[0001] This invention belongs to the field of atomic time technology, specifically relating to a method for establishing a composite clock based on an atomic clock and a millisecond pulsar. Background Technology

[0002] Since the 13th International Committee for Weights and Measures (CGPM) decided in 1967 to adopt the atomic second as the basic unit of time measurement, atomic time has become the main international time reference standard. Atomic time is typically generated by one or more atomic clocks; International Atomic Time (TAI) is established by approximately 500 atomic clocks in over 80 time laboratories worldwide. Atomic clocks offer nanosecond-level measurement accuracy, high resolution, customizable measurement intervals, and a guaranteed operating environment, enabling them to produce continuous atomic time with excellent short- to medium-term frequency stability. However, due to factors such as the physical structure, operating characteristics, seasonal variations, and annual effects of atomic clocks, the long-term stability of atomic time gradually decreases over time.

[0003] Pulsar time, established using millisecond pulsars as references, exhibits excellent long-term frequency stability, and its accuracy has been proven: Hobbs G. et al., in 2012, used timing data from 19 millisecond pulsars at the Parkes Observatory Millisecond Pulsar Timing Array (PPTA) in Australia to obtain a composite pulsar time comparable to the most accurate Earth time, TT (BIPM10). TT (BIPMXX) is another implementation of Earth time established annually at the beginning of the year by the International Bureau of Weights and Measures (BIPM) using data from all available reference frequencies and atomic clocks, employing a post-processing method. It is currently the most accurate and stable atomic time scale globally, with XX representing the last two digits of the year. However, due to limitations in observation equipment, techniques, and interstellar propagation delays, the timing measurement accuracy of the best-performing millisecond pulsars is only tens of nanoseconds, far lower than the accuracy of atomic clock measurements. Furthermore, the sampling resolution is poor, and current millisecond pulsar timing observations are spaced between days and months apart, with observations being discontinuous.

[0004] Millisecond pulsar timing observations use the atomic clock at the station as a reference. This reference is converted to more accurate and stable TT (BIPMXX) or TT (TAI) when calculating pulsars. Therefore, the performance of the station's atomic clock does not affect the establishment of pulsar time. That is, the establishment of pulsar time depends on one or more millisecond pulsars, which is completely different from the generation mechanism of atomic time that depends on a set of freely operating atomic clocks.

[0005] In summary, while atomic time is currently the most widely used high-precision time scale globally, exhibiting excellent short- to medium-term frequency stability, it also has certain limitations due to various factors. In particular, the lifespan of atomic clocks is currently no more than 20 years. Therefore, atomic time, which relies on one or more freely operating atomic clocks, struggles to maintain good accuracy and long-term frequency stability. Summary of the Invention

[0006] The technical problem to be solved by the present invention is to overcome the shortcomings of existing atomic time and provide a composite clock establishment method based on atomic clock and millisecond pulsar by using pulsar time to constrain and control the phase change of atomic time, thereby improving the long-term stability of atomic time, while maintaining the high sampling rate and excellent short-to-medium-term stability of atomic time.

[0007] The technical solution adopted to solve the above-mentioned technical problems is: a composite clock establishment method based on atomic clocks and millisecond pulsars, including the following steps:

[0008] A1. Read the timing observation data of the millisecond pulsar and use Tempo2 software to calculate the timing residual of the millisecond pulsar;

[0009] A2. The timing residuals of millisecond pulsars are sorted in reduced Julian day order;

[0010] A3. Select millisecond pulsars that meet the initial screening conditions for establishing millisecond pulsars;

[0011] A4. Determine whether the timing residual of the millisecond pulsar selected in step A3 is a multi-band timing residual. If it is a multi-band timing residual, proceed to step A5; otherwise, proceed to step A6.

[0012] A5. Statistical processing of multi-band timing residuals yields the average timing residual;

[0013] A6. Perform timing red noise identification and processing on the timing residuals that are not multi-band in step A4 and the average timing residuals obtained in step A5 to obtain the millisecond pulsar time series y(t).

[0014] y(t) = TT(BIPMXX) - PT(t)

[0015] In the formula, TT(BIPMXX) is the time scale established by the International Bureau of Weights and Measures at the beginning of each year, XX is the last two digits of the year, and PT(t) is the millisecond pulsar time at time t;

[0016] A7. Read the atomic clock measurement and comparison data, perform preprocessing such as frequency and phase jump identification and gross error detection, and obtain the atomic time series Y(t).

[0017] Y(t) = TT(BIPMXX) - AT(t)

[0018] In the formula, AT(t) represents the atomic time at time t;

[0019] The atomic time AT(t) at time t is:

[0020]

[0021] In the formula, p i h represents the weight of the i-th atomic clock. i (t) represents the reading of the i-th atomic clock at time t, a i,j (t k ) represents the time interval t within time period j. k The phase value of the i-th atomic clock at time i, B i,j (t) is [t k The frequency value of the i-th atomic clock in the interval [t], C i,j (t) is [t k The frequency drift of the i-th atomic clock in the interval [t], 1≤i≤M, where M is the number of atomic clocks;

[0022] A8. Based on the atomic time series obtained in step A7, the first derivative sequence of the atomic clock atomic time AT is obtained according to the following formula.

[0023]

[0024] In the formula, Let t be the sequence of first derivatives at the i-th atomic clock. i For the observation epoch of the i-th atomic clock, t i+1 This represents the observation epoch of the (i+1)th atomic clock;

[0025] A9. The Vondrak-Cepek combination method is used to combine pulsar time and atomic time to obtain the composite clock Q.

[0026]

[0027] In the formula, l is the time parameter, and ε is the smoothing factor for the PT time series of pulsars. y′ is the smoothing factor of the first derivative sequence of AT at atomic time. l The measured value of the PT time series at the l-th composite clock time is... q is the measured value of the first derivative sequence of AT at the l-th composite clock time. l The weights of the PT time series measurements at the l-th composite clock time are: y represents the weight of the first derivative sequence measurements of AT at the l-th composite clock time. l The smoothing value at the l-th composite clock time. Let be the smoothed value of the first derivative of the l-th composite clock, N be the total number of smoothed values, and n be the number of pulsar time series. Estimate t1 and t from discrete function smoothing values. N For curves The boundary values ​​of the independent variable;

[0028] The smoothing factor ε of the PT time series during pulsar observation is:

[0029]

[0030] In the formula, f is the frequency value of the periodic signal contained in the pulsar time and atomic time time series, and T is the ratio of the amplitude of the smoothed curve of the function with period P = 1 / f to the amplitude of the observed curve;

[0031] The smoothing factor of the first derivative sequence of atomic time AT for:

[0032]

[0033] In the formula, It is the ratio of the amplitude of the derivative smoothing curve with period P = 1 / f to the amplitude of the observed curve.

[0034] As a preferred technical solution, the method for screening millisecond pulsars in step A3 is as follows:

[0035] Step A31. Starting with the first millisecond pulsar, select millisecond pulsars with a timing residual greater than 10 years;

[0036] Step A32. Further filter out the uncertainties σ that satisfy the pulse arrival time. TOA Millisecond pulsars under conditions <1.5us.

[0037] As a preferred technical solution, the statistical processing method for multi-band timing residuals in step A5 is as follows:

[0038] Step A51. For the same time t, based on the uncertainty of the pulse arrival time of the observation data of different frequency bands, the weights of different frequency bands are obtained according to the following formula;

[0039]

[0040] In the formula, ω s The weight of the s-th frequency band is... Let be the uncertainty of the arrival time of the pulse in the s-th frequency band. Let S be the uncertainty of the arrival time of the pulse in the j-th frequency band, and S be the total number of frequency bands.

[0041] Step 52: Using a weighted average algorithm, obtain the average timing residual of the multi-band timing residuals at time t according to the following formula.

[0042]

[0043] In the formula, u s (t) represents the timing residual of the s-th frequency band.

[0044] As a preferred technical solution, the timing red noise identification and processing method in step A6 is as follows:

[0045] A61. σ is obtained according to the following formula. z (τ),

[0046]

[0047] In the formula, σ z (τ) represents the stability of the pulsar in the time domain, where τ is the sampling interval and c is the coefficient of the cubic term of the polynomial fitting over the τ interval.

[0048] A62. According to log(σ) z The (τ)-log(τ) curve is used to identify timing noise.

[0049] A63. Employ empirical mode decomposition to suppress timing red noise.

[0050] The beneficial effects of this invention are as follows:

[0051] This invention first selects millisecond pulsars with high accuracy in measuring pulse arrival time, long observation span of timing data, and no timing red noise based on the conditions for establishing pulsar time, thus establishing pulsar time with excellent long-term stability. Secondly, it employs the Vondrak-Cepek combination method to use pulsar time to constrain and control the phase changes of atomic time (AT), reducing or eliminating its long-term instability and thereby improving the long-term stability of atomic time. This invention introduces millisecond pulsars into my country's current time reference maintenance system, improving the accuracy and long-term stability of existing time references, and establishing a new, independent, secure, and robust time reference capable of real-time time services. It solves the problem of the long-term stability of atomic time gradually deteriorating over time, thereby addressing the issue of maintaining my country's time reference during periods of international time comparison disruption. This has significant strategic importance for improving my country's ability to maintain its time system independently, as well as for social stability and national defense security during emergencies. Attached Figure Description

[0052] Figure 1 This is a flowchart illustrating the present invention. Detailed Implementation

[0053] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments, but the present invention is not limited to the following embodiments.

[0054] exist Figure 1 The method for establishing a composite clock based on an atomic clock and a millisecond pulsar in this embodiment includes the following steps:

[0055] A1. Read the timing observation data of the millisecond pulsar and use Tempo2 software to calculate the timing residual of the millisecond pulsar;

[0056] A2. The timing residuals of millisecond pulsars are sorted in reduced Julian day order;

[0057] A3. Select millisecond pulsars that meet the initial screening conditions for establishing millisecond pulsars;

[0058] The method for screening millisecond pulsars in this embodiment is as follows:

[0059] Step A31. Starting with the first millisecond pulsar, select millisecond pulsars with a timing residual greater than 10 years;

[0060] Step A32. Further filter out the uncertainties σ that satisfy the pulse arrival time. TOA Millisecond pulsars with a time limit of <1.5µs;

[0061] A4. Determine whether the timing residual of the millisecond pulsar selected in step A3 is a multi-band timing residual. If it is a multi-band timing residual, proceed to step A5; otherwise, proceed to step A6.

[0062] A5. Statistical processing of multi-band timing residuals yields the average timing residual;

[0063] The statistical processing method for multi-band timing residuals is as follows:

[0064] Step A51. For the same time t, based on the uncertainty of the pulse arrival time of the observation data of different frequency bands, the weights of different frequency bands are obtained according to the following formula;

[0065]

[0066] In the formula, ω s The weight of the s-th frequency band is... Let be the uncertainty of the arrival time of the pulse in the s-th frequency band. Let S be the uncertainty of the arrival time of the pulse in the j-th frequency band, and S be the total number of frequency bands.

[0067] Step 52: Using a weighted average algorithm, obtain the average timing residual of the multi-band timing residuals at time t according to the following formula.

[0068]

[0069] In the formula, u s (t) represents the timing residual of the s-th frequency band;

[0070] A6. Perform timing red noise identification and processing on the timing residuals that are not multi-band in step A4 and the average timing residuals obtained in step A5 to obtain the millisecond pulsar time series y(t).

[0071] y(t) = TT(BIPMXX) - PT(t)

[0072] In the formula, TT(BIPMXX) is the time scale established by the International Bureau of Weights and Measures at the beginning of each year, XX is the last two digits of the year, and PT(t) is the millisecond pulsar time at time t;

[0073] The method for identifying and processing timing red noise is as follows:

[0074] A61. The time-domain characterization of pulsar stability σ is obtained according to the following formula. z (τ),

[0075]

[0076] In the formula, τ is the sampling interval time, and c is the coefficient of the cubic term of the polynomial fitting over the τ interval;

[0077] A62. According to log(σ) z The (τ)-log(τ) curve is used to identify timing noise.

[0078] A63. Employ empirical mode decomposition to suppress timing red noise;

[0079] A7. Read the atomic clock measurement and comparison data, perform preprocessing such as frequency and phase jump identification and gross error detection, and obtain the atomic clock time series Y(t).

[0080] Y(t) = TT(BIPMXX) - AT(t)

[0081] In the formula, AT(t) represents the atomic time at time t;

[0082] The atomic time AT(t) at time t is:

[0083]

[0084] In the formula, p i h represents the weight of the i-th atomic clock. i (t) represents the reading of the i-th atomic clock at time t, a i,j (t k ) represents the time interval t within time period j. k The phase value of the i-th atomic clock at time i, Bi,j (t) is [t k The frequency value of the i-th atomic clock in the interval [t], C i,j (t) is [t k The frequency drift of the i-th atomic clock in the interval [t], 1≤i≤M, where M is the number of atomic clocks;

[0085] A8. Based on the atomic time series obtained in step A7, the first derivative sequence of the atomic clock atomic time AT is obtained according to the following formula.

[0086]

[0087] In the formula, Let t be the sequence of the first derivatives of the i-th atomic clock. i For the observation epoch of the i-th atomic clock, t i+1 This represents the observation epoch of the (i+1)th atomic clock;

[0088] A9. The composite clock Q is obtained by combining pulsar time and atomic time using the Vondrak-Cepek combination method.

[0089]

[0090] In the formula, l is the time parameter, and ε is the smoothing factor for the PT time series of pulsars. y′ is the smoothing factor of the first derivative sequence of AT at atomic time. l The measured value of the PT time series at the l-th composite clock time is... q is the measured value of the first derivative sequence of AT at the l-th composite clock time. l The weights of the PT time series measurements at the l-th composite clock time are: y represents the weight of the first derivative sequence measurements of AT at the l-th composite clock time. l The smoothing value at the l-th composite clock time. Let be the smoothed value of the first derivative of the l-th composite clock, N be the total number of smoothed values, and n be the number of pulsar time series. Estimate t1 and t from discrete function smoothing values. N For curves The boundary values ​​of the independent variable;

[0091] when The composite clock time series is obtained.

[0092] The Vondrak-Cepek method in this embodiment is a combined smoothing method proposed by Vondrak and Cepek in 2000. It utilizes derivative information to smooth functions, that is, it incorporates derivative information into the function value, and through smoothing the signal, it obtains a smoothing result that has both the accuracy and long-term stability of the function value, as well as the high time resolution and short-term stability of the derivative value.

[0093] The smoothing factor ε of the pulsar-time PT time series in this embodiment is:

[0094]

[0095] In the formula, f is the frequency value of the periodic signal contained in the pulsar time and atomic time time series, and T is the ratio of the amplitude of the smoothed curve of the function with period P = 1 / f to the amplitude of the observed curve;

[0096] The smoothing factor of the first derivative sequence of atomic time AT in this embodiment for:

[0097]

[0098] In the formula, It is the ratio of the amplitude of the derivative smoothing curve with period P = 1 / f to the amplitude of the observed curve.

[0099] This invention applies the Vondrak-Cepek method to all available measurement data from both millisecond pulsar and atomic clock measurement technologies, establishing a composite clock that possesses both the long-term stability and accuracy of the former and the high resolution and short-to-medium-term stability of the latter.

Claims

1. A method for establishing a composite clock based on an atomic clock and a millisecond pulsar, characterized in that, Includes the following steps: A1. Read the timing observation data of the millisecond pulsar and use Tempo2 software to calculate the timing residual of the millisecond pulsar; A2. The timing residuals of millisecond pulsars are sorted in reduced Julian day order; A3. Select millisecond pulsars that meet the initial screening conditions for establishing millisecond pulsars; A4. Determine whether the timing residual of the millisecond pulsar selected in step A3 is a multi-band timing residual. If it is a multi-band timing residual, proceed to step A5; otherwise, proceed to step A6. A5. Statistical processing of multi-band timing residuals yields the average timing residual; A6. Perform timing red noise identification and processing on the timing residuals that are not multi-band in step A4 and the average timing residuals obtained in step A5 to obtain the millisecond pulsar time series y(t). y(t) = TT(BIPMXX) - PT(t) In the formula, TT(BIPMXX) is the time scale established by the International Bureau of Weights and Measures at the beginning of each year, XX is the last two digits of the year, and PT(t) is the millisecond pulsar time at time t; A7. Read the atomic clock measurement and comparison data, perform preprocessing such as frequency and phase jump identification and gross error detection, and obtain the atomic time series Y(t). Y(t) = TT(BIPMXX) - AT(t) In the formula, AT(t) represents the atomic time at time t; The atomic time AT(t) at time t is: In the formula, p i h represents the weight of the i-th atomic clock. i (t) represents the reading of the i-th atomic clock at time t, a i,j (t k ) represents the time interval t within time period j. k The phase value of the i-th atomic clock at time i, B i,j (t) is [t k The frequency value of the i-th atomic clock in the interval [t], C i,j (t) is [t k The frequency drift of the i-th atomic clock in the interval [t], 1≤i≤M, where M is the number of atomic clocks; A8. Based on the atomic time series obtained in step A7, the first derivative sequence of the atomic clock atomic time AT is obtained according to the following formula. In the formula, Let t be the sequence of first derivatives at the i-th atomic clock. i For the observation epoch of the i-th atomic clock, t i+1 This is the observation epoch of the (i+1)th atomic clock; A9. The Vondrak-Cepek combination method is used to combine the millisecond pulsar time and the atomic clock time to obtain the composite clock Q. In the formula, l is the time parameter, and ε is the smoothing factor for the PT time series of pulsars. y′ is the smoothing factor of the first derivative sequence of AT at atomic time. l The measured value of the PT time series at the l-th composite clock time is... q is the measured value of the first derivative sequence of AT at the l-th composite clock time. l The weights of the PT time series measurements at the l-th composite clock time are: y represents the weight of the first derivative sequence measurements of AT at the l-th composite clock time. l The smoothing value of the l-th composite clock. Let be the smoothed value of the first derivative of the l-th composite clock, N be the total number of smoothed values, and n be the number of pulsar time series. Estimate t1 and t from discrete function smoothing values. N For curves The boundary values ​​of the independent variable; The smoothing factor ε of the PT time series during pulsar observation is: In the formula, f is the frequency value of the periodic signal contained in the pulsar time and atomic time time series, and T is the ratio of the amplitude of the smoothed curve of the function with period P = 1 / f to the amplitude of the observed curve; The smoothing factor of the first derivative sequence of atomic time AT for: In the formula, It is the ratio of the amplitude of the derivative smoothing curve with period P = 1 / f to the amplitude of the observed curve.

2. The method for establishing a composite clock based on an atomic clock and a millisecond pulsar according to claim 1, characterized in that, The method for screening millisecond pulsars in step A3 is as follows: Step A31. Starting with the first millisecond pulsar, select millisecond pulsars with a timing residual greater than 10 years; Step A32. Further filter out the uncertainties σ that satisfy the pulse arrival time. TOA Millisecond pulsars under conditions <1.5us.

3. The method for establishing a composite clock based on an atomic clock and a millisecond pulsar according to claim 1, characterized in that, The statistical processing method for multi-band timing residuals in step A5 is as follows: Step A51. For the same time t, based on the uncertainty of the pulse arrival time of the observation data of different frequency bands, the weights of different frequency bands are obtained according to the following formula; In the formula, ω s The weight of the s-th frequency band is... Let be the uncertainty of the arrival time of the pulse in the s-th frequency band. Let S be the uncertainty of the arrival time of the pulse in the j-th frequency band, and S be the total number of frequency bands. Step 52: Using a weighted average algorithm, obtain the average timing residual of the multi-band timing residuals at time t according to the following formula. In the formula, u s (t) represents the timing residual of the s-th frequency band.

4. The method for establishing a composite clock based on an atomic clock and a millisecond pulsar according to claim 1, characterized in that, The timing red noise identification and processing method in step A6 is as follows: A61. σ is obtained according to the following formula. z (τ), In the formula, σ z (τ) represents the stability of the pulsar in the time domain, where τ is the sampling interval and c is the coefficient of the cubic term of the polynomial fitting over the τ interval. A62. According to log(σ) z The (τ)-log(τ) curve is used to identify timing noise. A63. Employ empirical mode decomposition to suppress timing red noise.

Images