CN110554597A - hydrogen cesium time scale fusion method based on Vondark-Cepek filtering - Google Patents

Abstract

The invention provides a hydrogen-cesium time scale fusion method based on Vondark-Cepek filtering. First, the exponential filter method is used to generate a clock group time scale for the hydrogen clock group and the cesium clock group respectively, and then Vondark-Cepek is calculated according to the principle of least squares. The key parameters of the Cepek combined filter are selected, and then the performance of the time scale of the cesium clock group is enhanced through the difference information of the time scale time series of the hydrogen clock group, so as to obtain the hydrogen-cesium fusion time scale. The present invention can take into account the characteristics of two types of hydrogen and cesium atomic clocks while effectively suppressing noise, and integrates at the clock group scale level to generate a time scale with good long-term and short-term stability.

Description

Hydrogen-Cesium Time Scale Fusion Method Based on Vondark-Cepek Filter

technical field

The invention belongs to the field of time frequency and signal processing, and relates to a time scale fusion method.

Background technique

At present, high-precision time has become a vital parameter in a country's science and technology, economy, military and social life. It is widely used in navigation, electric power, communication, aviation, national defense and other fields. As the most basic physical quantity, it is very important to improve the level of national scientific research. There is also a fundamental role. For this reason, every time-frequency laboratory in the world has its own set of atomic clocks, and each atomic clock can generate a time scale. But keeping time cannot rely on a single atomic clock, because every physical device has the potential for anomalies. The time scale algorithm is to calculate a standard time based on the time of each clock in the clock watch group according to the demand. The purpose is to find a way to minimize the uncertainty of the algorithm and maximize the stability.

Cesium hydrogen atomic clocks are two different types of atomic frequency standards. They have good short-term stability and long-term stability respectively. They are the two most common time-keeping atomic clocks in time laboratories around the world. The hydrogen-cesium time scale fusion method is mainly to study the effective combination of two different atomic clocks to produce a more stable and accurate time scale. What a time laboratory emphasizes is the stability of the time scale.

At present, many chronological laboratories use hydrogen and cesium as cross-reference methods to fuse the time scales of two different types of atomic clocks: using hydrogen clocks, cesium clocks, or cesium clock group time as references, and measuring the clock difference of another type of atomic clock to form After the time series, the short-term fluctuations of cesium clocks in the clock-keeping group are corrected or the rate and frequency drift of hydrogen clocks are deducted, and then the final time scale is obtained by weighted average. In the above-mentioned research method, a single clock of another atomic clock is mainly evaluated and corrected based on hydrogen or cesium, and then all the clocks in the clock-keeping group are weighted and averaged together to calculate the final time scale. The disadvantage of this method is that when the hydrogen atomic clock is used as a reference to measure the cesium atomic clock, the main noise is phase white noise. After filtering by mathematical methods, the short-term stability of its time scale will still be affected by the noise of the cesium atomic clock. However, when the hydrogen clock is corrected with the time scale of the cesium clock group, when the hydrogen clock breaks down, the reliability of the generated time scale cannot be guaranteed.

Contents of the invention

In order to overcome the deficiencies of the prior art, the present invention provides a hydrogen-cesium time scale fusion method based on Vondark-Cepek filtering, which can effectively suppress noise while taking into account the characteristics of two types of hydrogen-cesium atomic clocks, and perform fusion at the clock group scale level to generate A timescale with good long-term and short-term stability.

The technical solution adopted by the present invention to solve its technical problems comprises the following steps:

Step 1, grouping cesium atomic clocks and hydrogen atomic clocks, and generating exponentially filtered cesium clock group time scale TA _Cs and hydrogen clock group time scale TA _H through clock error prediction, frequency estimation and weight estimation steps respectively;

Step 2: Perform a difference on the time scale sequence of the hydrogen clock group to obtain the corresponding frequency data;

Step 3, transform the cesium clock group time scale sequence into the frequency domain by Fourier transform, and determine the signal f or period P to be suppressed by spectrum analysis; calculate the cesium clock group time scale smoothing factor according to the least square theory and hydrogen clock group time scale smoothing factor where T and T' are the frequency responses corresponding to TA _Cs and TA _H , respectively;

In step 4, the cesium clock group time scale and the hydrogen clock group frequency sequence are used as input, and the smoothing factor is used in the Vondark-Cepek combined filtering method to generate time scale fusion results and realize the fusion of the cesium clock group and hydrogen clock group time scales.

Before the step 1, preprocess the atomic clock difference data measured in the laboratory, select the atomic clock with the best long-term and short-term stability in the cesium atomic clock group and the hydrogen atomic clock group as the reference clock, and use the time interval counter to measure the same clock group The clock difference data between other atomic clocks and the reference clock.

The clock error data obtained by the preprocessing adopts the 3σ rule to remove abnormal points, and simultaneously detects the continuity of the clock error data, and uses interpolation to repair the clock error sequence with continuous missing data less than 5 measurement cycles; for the clock error data continuous If there are more than 5 measurement periods missing, the clock will be excluded from the clock-keeping group.

In the described step 1, if the weight of each atomic clock is greater than the upper limit of the set maximum weight, then the weight of this atomic clock is set to the maximum weight; let h' _i (t) be the time correction amount of the atomic clock _Hi at t, then set weights, atomic clock data, h' _i (t) via the formula calculation, where x _ij (t) is the clock difference between atomic clock H _j and atomic clock H _i at time t, N is the number of atomic clocks involved in the calculation, and x _i (t) is the average time scale between a single atomic clock in the laboratory and a group of atomic clocks in the laboratory The difference, w _i is the weight of each atomic clock.

In the step 1, the maximum weight is set to 2.5/N.

In the step 4, the cesium clock group time scale TA _Cs data is defined as m(t _i ), which is a function of time t _i , _i =1,2,3,..., and its derivative m'( t)=[m(t _i+1 )-m(t _i )]/(t _i+1 -t _i ); define that M(t _j ) is the data of hydrogen clock group time scale (TA _H ), and its derivative M′(t)=[M(t _j+1 )-M(t _j )]/(t _j+1 -t _j ), where j=1,2,3,…; when t _i =t _j When, m′(t)=M′(t), m(MJD2)=m(MJD1)+M′(t _j )Δt _j , MJD represents a simplified Julian day; suppose the input observation sequence 1 is expressed as y′ _j , the corresponding time scale is x _j , and the weight is p _j ; the first order derivative of observation sequence 2 is expressed as The corresponding time scale is Weight is The output sequence is the smooth curve y _i to be obtained, and the corresponding time scale is _xi , which defines the smoothness of the curve In the formula, The expression of is unknown, and its third derivative needs to be based on the smoothed data Estimated; a smooth curve between two points [x _i+1 , y _i+1 ] and [x _i+2 , y _i+2 ] is defined by the adjacent four points i, i+1, i A third-order Lagrange polynomial obtained by +2, i+3 its third derivative Defines the fidelity of the smooth curve to the observations Defines the fidelity of the first derivative of the smooth curve with respect to the observations

The beneficial effects of the present invention are: the fusion of hydrogen and cesium time scale levels is realized, and the absolute smoothness and absolute fidelity of the cesium clock group time scale and the absolute fitting of the first order derivative of the hydrogen clock group time scale are discussed. Good claim. On the basis of the exponential filtering algorithm, the filtering of the time scale of the clock group is further realized, which effectively improves the influence of the short-term fluctuation of cesium clocks on the fusion results. At the same time, the fusion of the two time scales further improves the reliability of the weighted average time scale only once. This method can not only effectively suppress the noise, but also make full use of the characteristics of the cesium hydrogen atomic clock, so that the long-term and short-term stability of the fusion time scale is good.

Description of drawings

Fig. 1 is a schematic flow chart of generating fusion time scales in the present invention (taking three cesium clocks and three hydrogen clocks as examples).

Figure 2 is a schematic diagram of the time scale, where (a) is the time scale of the cesium clock group, and (b) is the time scale of the hydrogen clock group.

Fig. 3 is a schematic diagram of time-scale fusion results of the present invention, wherein (a) is the fusion result, and (b) is the Allan deviation.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings and embodiments, and the present invention includes but not limited to the following embodiments.

The invention provides a hydrogen-cesium time scale fusion method based on Vondark-Cepek filtering. Firstly, an exponential filter method is used to generate a clock group time scale for each of the hydrogen clock group and the cesium clock group, and then according to the least squares principle, the Vondark-Cepek The key parameters of the combined filter are selected, and then the performance of the cesium clock group time scale is enhanced through the difference information of the hydrogen clock group time scale time scale, so as to obtain the hydrogen-cesium fusion time scale.

The present invention specifically comprises the following steps:

Step 1: Preprocessing the clock error data of the atomic clock measured in the laboratory. First, the atomic clock with the best long-term and short-term stability in the clock group is selected as the reference clock, and the time interval counter is used to measure the clock difference data between other atomic clocks and the reference clock. The clock difference is measured every hour. The 3σ rule is used to remove outliers, and the continuity of the clock error data is detected at the same time. The interpolation method is used to repair the clock difference sequence with continuous missing data of less than 5 hours; for the case of continuous missing clock difference data of more than 5 hours, the clock is removed from the clock keeping group.

Step 2: Group cesium atomic clocks and hydrogen atomic clocks into groups, and generate exponentially filtered cesium clock group time scale TA _Cs and hydrogen clock group time scale TA _H through clock error prediction, frequency estimation and weight estimation steps (calculated once per hour) respectively. (The following is an example of calculating the time scale of 3 months)

If the weight of each atomic clock is greater than the upper limit of the set maximum weight, the weight of the atomic clock is set to the maximum weight. According to the optimization of this method, our maximum weight is set to 2.5/N. Where N is the number of atomic clocks involved in the calculation. Suppose h' _i (t) is the time correction amount of atomic clock H _i at time t, then the weight, atomic clock data, and h' _i (t) are calculated by the following formula:

Among them, x _ij (t) is the clock difference between atomic clock H _j and atomic clock H _i at time t; x _i (t) is the difference between a single atomic clock in the laboratory and the average time scale of the laboratory atomic clock group, and w _i is each atomic clock the weight of.

Step 3: Perform a difference on the time scale sequence of the hydrogen clock set calculated in Step 2 to obtain the corresponding frequency data.

Step 4: Transform the cesium clock group time scale sequence into the frequency domain by Fourier transform, and determine the signal f or period P to be suppressed by spectrum analysis. Determine the frequency response T and T' corresponding to TA _Cs and TA _H according to the frequency f or period P, experience and actual algorithm, and calculate the time-scale smoothing factor ε of the cesium clock group and the time-scale smoothing factor of the hydrogen clock group according to the least square theory

Step 5: Take the cesium clock group time scale generated by the calculation in step 2 and the hydrogen clock group frequency sequence generated in step 3 as input, and apply the smoothing factor obtained in step 4 to the Vondark-Cepek combined filtering method to generate time scale fusion results and realize Fusion of cesium clock group and hydrogen clock group time scales.

Establish a hydrogen-cesium time scale fusion model: define the cesium clock group time scale (TA _Cs ) data as m(t _i ), which is a function of time t _i , i=1,2,3,..., its derivative at time t _i It is defined as m′(t)=[m(t _i+1 )−m(t _i )]/(t _i+1 −t _i ). M(t _j ) is the data of hydrogen clock group time scale (TA _H ), and its derivative is defined as M′(t)=[M(t _j+1 )-M(t _j )]/(t _{j+1 -tj} ). Among them, j=1, 2, 3, . . . The physical meanings of m'(t) and M'(t) are the same, they are both the rate of the time scale sequence. Therefore, when t _i =t _j , m'(t)=M'(t).

m(MJD2)=m(MJD1)+M′(t _j )Δt _j (3)

where MJD is the simplified Julian day.

Because the time intervals of the time scales of the cesium clock group and the hydrogen clock group are the same, the formula (3) can be used to estimate the clock difference of the cesium clock group at the next moment by using the clock difference of the cesium clock group and the time scale rate and time of the hydrogen clock group at this moment The sum of interval products representation.

Assume that the input observation sequence 1 is expressed as y′ _j , the corresponding time scale is x _j , and the weight is p _j ; the first order derivative of observation sequence 2 is expressed as The corresponding time scale is Weight is The output sequence is the smooth curve to be obtained as y _i , and the corresponding time scale is _xi . Define three quantities:

1) The smoothness of the curve:

In the formula, The expression of is unknown, and its third derivative needs to be based on the smoothed data Make an estimate. A smooth curve between two points [x _i+1 , y _i+1 ] and [x _i+2 , y _i+2 ] is defined by the adjacent four points i, i+1, i+2, A third-order Lagrangian polynomial L _i (x) obtained by i+3:

Its third derivative is:

The third derivative between each pair of data points is set as a constant, formula (1) can be expressed as:

in

This smooth definition says that the ideal smooth function is a quadratic function of time (whose first derivative is a linear function).

2) The fidelity of the smooth curve to the observed values:

3) The fidelity of the smooth curve to the first derivative of the observed value:

First Derivative The smooth value of can be expressed in terms of the smooth function value y _i .

Using the first derivative of the Lagrangian polynomial L _i (x) used to define S, denoted as L′ _i (x):

L' _i (x)＝A _i (x)y _i +B _i (x)y _i+1 +C _i (x)y _i+2 +D _i (x)y _i+3 (10)

in,

is the smooth value of the first derivative represented by the smooth function value y _i There is a certain degree of freedom in the selection of the four points around x _i at each moment. The constraints of the following formula guarantee that is on a smooth curve.

1) In order to calculate Using the values of time x ₁ , x ₂ , x ₃ , x ₄ , that is to say, formula (11) can be obtained

in,

2) For the first half of the input data The calculation of i=2, 3, ... N/2, using time _xi-1 , _xi , _xi+1 , _xi+2 , namely

in,

3) For the second half of the input data The calculation of i=N/2+1, N/2+2,...N-1, using time _xi-2 , _xi-1 , _xi , _xi+1 , namely

in,

4) For The calculation of using instants x _N-3 , x _N-2 , x _N-1 , x _N , namely

in,

but It can be expressed as:

The VC combined filtering method aims to determine the equalization of smoothed values _yi under three different conditions, each different setting leads to different results:

(a) The curve needs to be smoothed (minimize S)

(b) Smoothing requires observations close to the function (minimizing F)

(c) The first derivative of a smooth curve requires observations close to the first derivative (minimizing )

The least squares method is used to minimize the above constraints to adjust each parameter, which is expressed as:

where ε≥0, The degree of equalization of the three conditions can be determined by choosing the values of these two parameters: smoothing coefficient ε and The larger the value of the smoothing coefficient, the greater the weight of the fidelity relative to the observed function value or its first derivative, and the closer the smoothed value is to the observed value.

An embodiment of the present invention provides a method for fusion of real hydrogen and cesium hydrogen time scales by using the Vondark-Cepek combined filter. For the first time, the fusion of the time-scale level of the hydrogen-cesium clock group is proposed, which avoids the problem that the white noise generated by the hydrogen-cesium cross-reference affects the time-scale results. The hydrogen clock group time scale and the cesium clock group time scale have good short-term and long-term stability respectively, so without affecting the long-term stability of the cesium clock group time scale, the excellent short-term stability of the hydrogen clock group time scale Performance enhancements to produce fusion timescales with good long and short stability.

As shown in Figure 1, the method for generating the fusion time scale includes the following steps:

Step 1. Preprocess the atomic clock error data measured in the laboratory, and use the 3σ rule to remove abnormal points.

The gross error elimination is based on the 3σ criterion, the detected atomic clock frequency difference data is {x ₁ , x ₂ , x ₃ ,…, x _N-1 , x _N }, and the sample mean is The standard deviation is The 3σ criterion is that if the measured value Then this value is a singular point and should be eliminated.

Step 2. Group cesium atomic clocks and hydrogen atomic clocks, and use the clock difference data processed in step 1 to calculate the time scale of each clock group by using the exponential filtering method to obtain TA _Cs and TA _H , as shown in Figure 2.

Clock error estimation: the clock error data X _ij (t) between each watch clock and the reference clock, and the measurement interval is τ. The estimated value of clock error at time (t+τ) is calculated by formula (1) based on the clock error and clock velocity at time t. The optimal estimate of the clock error of the reference clock relative to the time scale is calculated by formula (2) based on the principle of weighted average.

Among them, X _i (t), Y _i (t) are respectively the estimation of time difference and frequency difference of clock i relative to a certain reference time scale at time t; is the clock error forecast of clock i at time t.

Frequency estimation: The average frequency difference estimation of clock i at (t+τ) time is obtained by using the first difference of clock difference data, and then calculated by exponential filter, such as formulas (3) and (4).

in is the frequency difference forecast of clock _i in the time interval from time t to t+τ; mi is the average time constant of exponential frequency, which is the error introduced in minimizing the time difference estimation formula (1), and is estimated according to formula (5).

Among them, τ _mini is the most stable time interval of clock i, which is taken on the Allan variance curve The value of τ at the minimum.

Weight estimation: use the difference between the clock error of each clock relative to the paper clock and the calculated value at the current time to determine the weight.

Such as formula (6):

Where K _i is the deviation of the error estimate, calculated according to formula (10). The mean square time error of each clock is determined by the exponential filter of equation (7).

where ε _i (τ) is the cumulative error of time difference estimation of clock i at time interval τ; <ε ² _x (τ)> is the paper time mean square error at time t at time interval τ, and the initial value is generally taken as N _τ is the time constant of the exponential filter, which is used to estimate the current mean square error, and N _τ is generally 20 days.

n is the number of atomic clocks. The weights are applied in Equation (2), and their calculation is shown in Equation (9).

Firstly, the clock difference file is read, and the initial values of _Xi (t) and _Yi (t) are taken. For the selection of the initial value, the clock difference value of the atomic clock relative to UTC at the initial time is selected for convenience. The T Bulletin of BIPM publishes the clock difference between UTC and UTC (NTSC). It is known that the clock difference of the three hydrogen clocks relative to UTC (NTSC) at this time is subtracted from the two differences to obtain the clock difference of the three hydrogen clocks relative to UTC. Difference. The initial values of the clock velocity Y _i (t) and the drift parameter d are obtained by quadratic fitting of the clock error data of the previous ten days. Secondly, according to the formula (1), calculate the predicted clock difference at the next moment then calculate The initial value of is usually taken as Then the Allan variance of the three clocks should be calculated first. Since what we measure is the clock difference between each species and the main clock, we need to use the triangular hat method to solve the Allen variance of each clock, so we need to subtract the known clock difference data two by two, and remove the main clock. Obtain the clock difference between the three hydrogen atomic clocks. Then calculate the Allan variance between two atomic clocks, and solve the Allan variance of each clock according to the triangular hat calculation formula (11), according to the formula calculate initial value. Then calculate the initial value of the weight according to the formula (8-9). Then, according to the formula (2) and the clock data of each atomic clock, calculate the clock error of the main clock relative to the paper clock at the next moment, and perform a weighted average to obtain the final clock error estimate X _j (t+τ), which is the calculated The difference between the mean time and the reference clock at instant (t+τ), ie TA _Cs or TA _H . Then, according to the estimated clock error, the clock speed forecast for the next moment (t+τ) is carried out. Calculate the clock difference between each clock and the paper clock according to X _i =X _j +X _ij , and substitute it into formula (3) to calculate the frequency difference estimate. Then calculate the exponential frequency constant according to the formula (5), and substitute it into the formula (4) to calculate Y _i (t+τ), which is used for the calculation of the clock error forecast at the next time (t+2τ). Finally, the clock error, weight, and frequency estimation of the next time (t+2τ) are performed. At this time, the time difference estimate is the same as the previous calculation method, and the weight calculation is calculated according to the formula (6)-(10). First, the parameter K _i is determined according to the formula (10), and then according to the formula (6)-(9). Take power. Loop in turn, and then carry out the prediction calculation of the next moment.

Step 3. Perform a difference on the hydrogen clock set time scale sequence calculated in step 2 to obtain the corresponding frequency data.

Step 4. Determine the frequency response T/T' according to the frequency f or period P of the noise to be reduced, and calculate the smoothing factor ε and

The cesium clock group time scale sequence TA _Cs is transformed into the frequency domain by Fourier transform, the noise period P or frequency f is determined by spectrum analysis, and the frequency response T/T' is determined at the same time, and finally the smoothing factor is calculated according to the formula.

Step 5. The time scale of the cesium clock group generated by the calculation in step 2 and the frequency sequence of the hydrogen clock group generated in step 3 are used as input, and the smoothing factor obtained in step 4 is used in the Vondark-Cepek combined filtering method to generate a time scale fusion result.

First, the system equation In matrix form it can be expressed as

in,

in,

Since A is a positive definite matrix, Cholesky decomposition can be used to solve the system equation (14). Matrix A can be decomposed by Cholesky:

A = U ^T DU (15)

Let U ^T w = By', the vector w can be obtained. Then according to Uy=D ^-1 w, y can be obtained. which is:

y=A ^-1 By' (16)

Therefore, the time-scale result of hydrogen-cesium fusion can be obtained, that is, the sequence y. The results are shown in Figure 3, (a) is the time-scale fusion result, and (b) is the Allan deviation.

The above implementations are only exemplary embodiments of the present invention, and are not intended to limit the present invention. The protection scope of the present invention is defined by the claims. Those skilled in the art can make various modifications or equivalent replacements to the present invention within the spirit and protection scope of the present invention, and such modifications or equivalent replacements should also be deemed to fall within the protection scope of the present invention.

Claims (6)

1. a hydrogen-cesium time-scale fusion method based on Vondark-Cepek filtering, is characterized in that comprising the following steps: Step 1, grouping cesium atomic clocks and hydrogen atomic clocks, and generating exponentially filtered cesium clock group time scale TA _Cs and hydrogen clock group time scale TA _H through clock error prediction, frequency estimation and weight estimation steps respectively; Step 2: Perform a difference on the time scale sequence of the hydrogen clock group to obtain the corresponding frequency data; Step 3, transform the cesium clock group time scale sequence into the frequency domain by Fourier transform, and determine the signal f or period P to be suppressed by spectrum analysis; calculate the cesium clock group time scale smoothing factor according to the least square theory and hydrogen clock group time scale smoothing factor where T and T' are the frequency responses corresponding to TA _Cs and TA _H , respectively; In step 4, the cesium clock group time scale and the hydrogen clock group frequency sequence are used as input, and the smoothing factor is used in the Vondark-Cepek combined filtering method to generate time scale fusion results and realize the fusion of the cesium clock group and hydrogen clock group time scales. 2. the hydrogen cesium time scale fusion method based on Vondark-Cepek filtering according to claim 1, is characterized in that: before described step 1, the atomic clock difference data that laboratory records is carried out preprocessing, selects cesium atomic clock respectively The atomic clock with the best long-term and short-term stability in the hydrogen atomic clock group and the hydrogen atomic clock group is used as the reference clock, and the clock difference data between other atomic clocks in the same clock group and the reference clock is measured by the time interval counter. 3. the hydrogen-cesium time scale fusion method based on Vondark-Cepek filtering according to claim 2, is characterized in that: the clock difference data that described preprocessing obtains adopts 3σ rule to remove anomalous point, detects the continuity of clock difference data simultaneously For the clock error sequence whose continuous missing data is less than 5 measurement periods, the interpolation method is used to repair; for the clock error data whose continuous missing data is greater than 5 measurement periods, the clock is removed from the clock keeping group. 4. the hydrogen-cesium time scale fusion method based on Vondark-Cepek filtering according to claim 1, is characterized in that: in described step 1, if the weight of each atomic clock is greater than the maximum weight upper limit of setting, then this atomic clock The weight of is set as the maximum weight; let h' _i (t) be the time correction amount of atomic clock H _i at time t, then the weight, atomic clock data, h' _i (t) are passed through the formula calculation, where x _ij (t) is the clock difference between atomic clock H _j and atomic clock H _i at time t, N is the number of atomic clocks involved in the calculation, and x _i (t) is the average time scale between a single atomic clock in the laboratory and a group of atomic clocks in the laboratory The difference, w _i is the weight of each atomic clock. 5. the hydrogen cesium time scale fusion method based on Vondark-Cepek filtering according to claim 4, is characterized in that: in the described step 1, the maximum weight is set to 2.5/N. 6. the hydrogen cesium time scale fusion method based on Vondark-Cepek filtering according to claim 1, is characterized in that: in described step 4, define cesium clock group time scale TA _Cs data as m (t _i ), is The function of time t _i , _i =1, 2, 3, ..., its derivative m'(t)=[m(t _i+1 )-m(t _i )]/(t _{i +1} -t _i ); define M(t _j ) as the data of hydrogen clock group time scale (TA _H ), its derivative M′(t)=[M(t _j+1 )-M(t _j )]/ (t _j+1 - _j ), where j=1, 2, 3,...; when t _i =t _j , m'(t)=M'(t), m(MJD2)=m( MJD1)+M′(t _j )Δt _j , MJD represents a simplified Julian day; suppose the input observation sequence 1 is expressed as y′ _j , the corresponding time scale is x _j , and the weight is p _j ; the first derivative of observation sequence 2 is expressed as The corresponding time scale is Weight is The output sequence is the smooth curve y _i to be obtained, and the corresponding time scale is _xi , which defines the smoothness of the curve In the formula, The expression of is unknown, and its third derivative needs to be based on the smoothed data Estimated; a smooth curve between two points [x _i+1 , y _i+1 ] and [x _i+2 , y _i+2 ] is defined by the adjacent four points i, i+1, i A third-order Lagrange polynomial obtained by +2, i+3 its third derivative Defines the fidelity of the smooth curve to the observations Defines the fidelity of the first derivative of the smooth curve with respect to the observations