Simulation residual determination method based on limited observation data

Abstract

The application discloses a kind of based on limited observation data's simulation residual determination method, comprising: pulsar observation data is obtained after processing timing residual, timing residual includes time value and timing residual value;Timing residual value is arranged from big to small, corresponding time value is also arranged;The number of sample timing residual value in different interval range is counted, and the frequency distribution diagram of sample timing residual value in different interval is drawn;Determine the probability statistical distribution model to which frequency distribution belongs;Specific parameters of probability statistical distribution model are solved;The number of simulation residual value is determined;According to probability statistical distribution model, simulation is obtained simulation residual value, set corresponding time value for simulation residual value, simulation residual value is added to sample observation data.The application establishes probability statistical model based on actual data according to the actual residual data obtained, and then simulates to obtain simulation residual value according to the model, which is closer to the original sample value, and improves the accuracy of pulsar time.

Description

Technical Field

[0001] This invention belongs to the field of pulsar timing technology and relates to a method for determining simulation residuals based on limited observation data. Background Technology

[0002] In today's era of rapid advancements in science and technology, precise time systems are a crucial foundation for fields such as atomic energy, aerospace technology, and high-energy physics. Existing time systems primarily rely on the standard second provided by atomic time for maintenance and operation. Within a specific region or even globally, atomic time enables high-precision timekeeping, time synchronization, and time and frequency measurement services, achieving and maintaining a unified time system. However, because atomic clocks gradually age and are susceptible to frequency deviations caused by external interference, their stability gradually declines over long timescales. Therefore, other high-precision time systems are needed to provide references for the calibration of atomic time, allowing for its correction.

[0003] Pulsar time is a time system based on the highly stable rotation frequency of pulsars, and its long-term stability can reach 10⁻⁶. -18 In summary, by employing appropriate analytical methods to process long-term timing observation data from multiple millisecond pulsars, a comprehensive pulsar time can be established. This time exhibits higher stability than that of a single pulsar time. Furthermore, since pulsars, as natural celestial bodies, have physical mechanisms completely different from atomic clocks, they can provide an effective reference for the correction of atomic time. This timeframe allows pulsars to complement and calibrate with the atomic time system over long timescales, thus contributing to the correction of atomic time.

[0004] Constructing a composite pulsar timeline and performing stability analysis requires a large amount of observational data, typically obtained from pulsar observations using radio telescopes and space X-ray detectors. However, the numerous observational tasks of radio telescopes and X-ray detectors make it difficult to ensure timely and quantitative observations of all pulsars. This results in inconsistent observation times for different pulsars, and a relatively small number of observation samples for the same pulsar. Therefore, in the composite pulsar timeline data processing, several pulsars must first be selected. Based on the measured residual data of these pulsars, timing residuals at equal intervals and times are generated through simulation. Existing techniques generally use interpolation to obtain simulated residual values. However, due to the non-uniform distribution of observation times, when two observation times are far apart, the timing residual obtained through interpolation at the intermediate time has already deviated significantly from the actual residual, and the error becomes very significant. Consequently, the calculated pulsar time stability also differs greatly from the true value. The problem this invention aims to solve is: addressing the issues caused by non-uniform sampling, abandoning traditional interpolation methods, and obtaining a probability distribution model by analyzing the observed data using statistical methods, thereby obtaining more accurate simulated residual values.

[0005] Currently, when assessing the stability of composite pulsars, given the limited number of actual observation samples and the uneven distribution of time values, interpolation methods are generally used to address this issue. The core idea of ​​this method is to use the interpolation function with the smallest error across different intervals, thereby obtaining simulated residual values ​​with uniform time distribution and small errors. Commonly used interpolation methods include linear interpolation, nearest neighbor interpolation, and cubic spline interpolation. Specifically, this method sets the time values ​​as the X-axis and the residual values ​​as the Y-axis. When the difference between two adjacent time values ​​is 10... 4 Time intervals of 1 second or longer are considered as discontinuities. These discontinuities divide the time axis into sub-intervals. The time values ​​and corresponding residual values ​​within each sub-interval are considered as a subsequence. Appropriate interpolation functions are used in different sub-sequences to obtain the simulated residual values ​​corresponding to the time values ​​within the discontinuities. The obtained simulated residual values ​​are then filtered to remove those that do not meet the requirements.

[0006] Since pulsar observation data mainly relies on radio telescopes, and radio telescopes are heavily loaded with observations every day, they cannot provide enough pulsar observation data. Furthermore, the traditional methods for generating analog timing residuals have the following drawbacks:

[0007] (1) The simulation residual value obtained by direct interpolation has a large error. This error is mainly caused by the interpolation method. Due to the characteristics of pulsars, the arrival time of the pulses received by radio telescopes is discontinuous, and there are many cases where the interval between adjacent time moments is large in the data. This leads to a large error when using the interpolation method to generate simulation residuals.

[0008] (2) The process of generating residuals is quite complicated. If you want to reduce the error of the residual values ​​obtained by interpolation, you need to segment the time series and give the segmentation criteria in advance, that is, the critical value of the interval between adjacent time points. If the interval between adjacent time points is greater than this critical value, then segmentation is performed. After segmentation, different interpolation methods may be used for different subsequences. Finally, data screening is required to remove residual data with large errors.

[0009] (3) Traditional interpolation methods also need to discard some actual observation samples. When the number of samples is already small, the observation data is wasted, which makes the overall error larger. Summary of the Invention

[0010] To address the aforementioned issues, this invention provides a method for determining simulated residuals based on limited observation data. A probabilistic statistical model based on the actual residual data is established, and then the simulated residual values ​​are obtained through simulation using this model. This method more closely approximates the original sample values, improving the accuracy of pulsar timekeeping and resolving the problems existing in the prior art.

[0011] The technical solution adopted in this invention is a method for determining simulation residuals based on finite observation data, comprising the following steps:

[0012] S1. Pulsar observation data is processed to obtain timing residuals, which include time values ​​and timing residual values. The timing residual values ​​are arranged from largest to smallest, and the corresponding time values ​​are also arranged. The timing residual values ​​are denoted as x1, x2, ..., x n The corresponding time values ​​are set as t1, t2, ..., t n ;

[0013] S2, [x1, x n The interval is divided into equally spaced small interval blocks. The number of sample timing residual values ​​within different interval ranges is counted, and the frequency distribution of sample timing residual values ​​in different intervals is plotted.

[0014] S3, determine the probability and statistical distribution model to which the frequency distribution belongs;

[0015] S4. Based on the timing residual values ​​and the number of residuals in the frequency distribution diagram, the specific parameters of the probability statistical distribution model are obtained, and the probability statistical distribution model of the sample timing residual values ​​is established.

[0016] S5, determine the number of simulation residuals; the number of simulation residuals is the minimum number of residuals for each pulsar when assembling the composite pulsar;

[0017] S6. Based on the probability statistical distribution model, simulated residual values ​​are obtained. Corresponding time values ​​are set for the simulated residual values. The simulated residual values ​​are added to the sample observation data. The time values ​​of the simulated residuals correspond to the time values ​​of the pulsar observation data.

[0018] Furthermore, in step S1, the pulsar observation data is processed by Tempo2 software to obtain timing residuals with a processing accuracy of up to 1 ns.

[0019] Furthermore, in step S2, any interval block is represented as: [x1+(i-1)×d,x1+i×d] (i=1,2,3...,n). Statistics are performed within each small interval block. If the timing residual value x... j If (j=1,2,...,n) is within this interval, then the timing residual value x within this interval block is... j The corresponding quantity y j The value is 1 if the interval is 1 and 0 otherwise, so the corresponding vertical coordinate value Y for this interval on the frequency distribution map is... j The expression for (j = 1, 2, ..., n) is:

[0020]

[0021] The number of sample timing residuals in different interval ranges is determined by equation (4-2).

[0022] Furthermore, the frequency distribution of the sample timing residual values ​​in different intervals is represented by a histogram, with the horizontal axis representing the distribution interval of the timing residual values, each node being an endpoint of the interval, and the vertical axis representing the number of sample timing residual values ​​in the corresponding interval.

[0023] Furthermore, step S3 specifically involves: observing the histogram to preliminarily determine which probability distribution model it belongs to, which can be any one of Gaussian distribution, uniform distribution, geometric distribution, or Poisson distribution; different models correspond to different testing methods, and the test results are used to determine whether the preliminary determination is correct, and finally the probability statistical distribution model corresponding to the timing residual data of the sample can be determined.

[0024] Furthermore, in step S4, the horizontal axis of the histogram is actually a series of small interval blocks, not specific numerical values. Therefore, it is necessary to first convert the interval blocks into specific values. Let any interval block be [x1+(i-1)×d,x1+i×d] (i=1,2,...,n), and the corresponding Y... i Let Z be the number of residual values ​​within this range; for the i-th interval block, let Z be the specific value corresponding to this interval block. i for:

[0025] Z i =[(x1+(i-1)×d)+(x1+i×d)] / 2 (4-6)

[0026] Let there be a set of data (Z) i ,Y i (i = 1, 2, ..., n), each Z i Corresponding to a Y i Then, the parameters can be substituted into the function of the determined probability and statistical distribution model for fitting, and the specific parameters of the probability and statistical distribution model can be obtained.

[0027] Furthermore, the probability and statistical distribution model is a Gaussian distribution, and the functional form of the Gaussian distribution is:

[0028]

[0029] Each Z i Corresponding to a Y i We substitute the function into the Gaussian distribution and fit it. The fitting function is:

[0030]

[0031] Y in the formula maxμ, σ, and μ' are the peak value, mean, and standard deviation of the Gaussian curve, respectively; taking the natural logarithm of both sides of equation (4-7):

[0032]

[0033] make:

[0034]

[0035] Since there is more than one data point, the above formula is represented in matrix form:

[0036]

[0037] The above formula can be described in vector form as follows:

[0038] K = ZB (4-11)

[0039] The vectors in the formula are represented as follows:

[0040]

[0041] According to the least squares method, the parameter B to be solved can be expressed as:

[0042] B = (Z) T Z) -1 Z T K (4-13)

[0043] After finding b0, b1, and b2, substitute them into equation (4-9) to obtain the parameter values ​​of the Gaussian model.

[0044] Furthermore, in step S5, let the number of residual values ​​for each pulsar be n1, n2, n3, n4, ..., n k Then the number of simulated residuals, n, is:

[0045] n = min[n1, n2, n3, n4, ..., n k ].

[0046] The beneficial effects of this invention are:

[0047] 1. The present invention demonstrates that by establishing an accurate statistical model, the simulated residual values ​​obtained from the simulation are all based on the probability statistical distribution model. Therefore, the residual values ​​obtained from the present invention highly overlap with the sample values, thus avoiding the simulation errors present in traditional interpolation methods.

[0048] 2. The embodiments of the present invention only require generating a residual distribution histogram, establishing a probability and statistics type, obtaining model parameter values, and performing simulation to generate simulated residuals. Each step has very mature functions available for use, which makes the method very convenient and overcomes the problem of cumbersome operation process in traditional interpolation methods.

[0049] 3. In order to obtain a more accurate statistical model, the embodiments of the present invention will make full use of every observation data point and will not discard the already scarce observation data, thus fully reflecting the value of each observation data point and improving the accuracy of pulsar timing. Therefore, the problem solved by the embodiments of the present invention is the most common problem in practical engineering applications and is also the bottleneck problem restricting the extraction of comprehensive pulsar timing, and has important engineering application value.

[0050] 4. In practical applications, the embodiments of the present invention only require the establishment of a probability model once, and the simulation residuals can be generated based on the model. In contrast, traditional methods require different interpolation methods in different intervals. Therefore, compared with traditional methods, the embodiments of the present invention can greatly reduce simulation time and have practical engineering significance. Attached Figure Description

[0051] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0052] Figure 1 This is a flowchart of an embodiment of the present invention.

[0053] Figure 2 This is a statistical distribution chart of residual values ​​in an embodiment of the present invention.

[0054] Figure 3 This is a comparison chart of simulated residual values ​​and actual residual values ​​in an embodiment of the present invention.

[0055] Figure 4 It is a comparison chart of the values ​​obtained by spline interpolation and the original values.

[0056] Figure 5 This is a comparison chart of the values ​​obtained through linear interpolation and the original values.

[0057] Figure 6 This is a comparison chart of the values ​​obtained by the nearest interpolation and the original values. Detailed Implementation

[0058] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0059] Example

[0060] A method for determining simulated residuals based on limited observation data abandons the traditional interpolation approach and instead obtains simulated residual values ​​through probability statistical distribution analysis of the sample observation data, such as... Figure 1 As shown, taking pulsar PSR_B1855+09 as an example, the steps include:

[0061] S1, the timing residuals are generally obtained using Tempo2 software, which is widely used in astronomy. This software has powerful data processing capabilities and can quickly fit the raw data to obtain the timing residuals with a processing accuracy of up to 1 ns. The data obtained after fitting by this software is divided into two columns: the first column is the time value, and the second column is the timing residual value corresponding to the time. To facilitate subsequent calculations, the timing residual values ​​are arranged from largest to smallest, and their corresponding time values ​​are also arranged, denoted as x1, x2, ..., x... n The corresponding time values ​​are set as t1, t2, ..., t n .

[0062] S2, after obtaining the timing residuals, the distribution of the residual values ​​is visually displayed in the form of a histogram, showing [x1, x...]. n The interval is divided into equally spaced small interval blocks. Let the number of interval blocks be n (in some embodiments, if the number of timing residual values ​​is greater than 10000, then n is set to 1000; otherwise, n is set to 500). Then the length d of the interval block is as follows:

[0063] d=(x n -x1) / n (4-1)

[0064] Therefore, any interval block can be represented as: [x1+(i-1)×d,x1+i×d] (i=1,2,3...,n). Statistics are performed within each small interval block. If the timing residual value x... j If (j=1,2,...,n) is within this interval, then the timing residual value x within this interval block is... j The corresponding quantity y j The value is 1 if the interval is 1 and 0 otherwise, so the corresponding ordinate value Y on the histogram is 1. j The expression for (j = 1, 2, ..., n) is:

[0065]

[0066] A Y value can be obtained for each interval block (i.e., the sum of all timing residual values ​​within the corresponding interval block). After all intervals are statistically analyzed, a distribution chart of timing residual values ​​is obtained, also called a histogram; for example... Figure 2 The figure shows a probability distribution histogram obtained by statistically analyzing the residual values ​​of the PSR_B1855+09 pulsar. The horizontal axis represents the distribution interval of the timing residuals, with each node being an endpoint of the interval. The vertical axis represents the number of timing residual values ​​within that interval.

[0067] S3. After drawing the histogram, the statistical model of the sample observation data is also established. The histogram counts the number of sample residuals in different intervals. Through this histogram, we can see the frequency of sample residual values ​​in different intervals.

[0068] After the histogram is drawn, it is necessary to analyze which probability distribution model it belongs to. The model analysis is divided into two steps. The first step is to observe the histogram and preliminarily determine which probability distribution model it belongs to. Common probability distribution models include Gaussian distribution, uniform distribution, geometric distribution, Poisson distribution, etc. The second step is to test the preliminary determination. Different models correspond to different testing methods. The test results are used to judge whether the preliminary determination is correct. Finally, the probability distribution model corresponding to the sample residual data can be determined.

[0069] Taking pulsar PSR_B1855+09 as an example, Figure 2 A histogram of the timing residuals of this pulsar is plotted. By observing this plot, it is preliminarily determined that it follows a Gaussian distribution model. Subsequent testing is then conducted to determine if the conjecture is correct. There are many methods for testing whether a Gaussian distribution exists, such as drawing a QQ plot or the Shapiro-Wilk test. This section details the Shapiro-Wilk test, a correlation-based algorithm proposed by Shapiro and Wilke in 1965. Specifically, let the sample residuals be x1, x2, ..., x... n The calculation formula is as follows:

[0070]

[0071] Where, x (i) Let x represent the i-th order statistic, i.e., the i-th smallest number in the sample; i Let x represent the residual value of the i-th sample, and x be the mean of the sample residuals; constant a i Calculated using equation (4-4);

[0072]

[0073] Where m=(m1,...,m n ) T ,m1,...,m n Let W represent the expected value of ordered, independent, and identically distributed statistics sampled from a standard normally distributed random variable, and V be the covariance of these ordered statistics. Equations (4-3) and (4-4) are the standard formulas for the Shapiro-Wilk test. The closer the value of the test statistic W is to 1, the more normally it conforms to the distribution.

[0074] S4. After determining the type of distribution model, to obtain the timing residuals using that model for simulation, the specific parameters of the model need to be determined. Continuing from the example in the previous step, after verifying that the statistical model for the pulsar timing residuals is a Gaussian distribution model, its mean and variance need to be determined to truly establish the model and then simulate to obtain the required residuals. The functional form of the Gaussian distribution is:

[0075]

[0076] Where A is the amplitude, μ is the mean, and σ is the standard deviation. To find these three unknowns, we need to substitute the residuals and the number of residuals in the histogram. However, the histogram counts the number of sample residuals within different intervals, and its horizontal axis is actually a series of small interval blocks, not specific values. Therefore, we need to first convert the interval blocks into specific values. Let any interval block be [x1+(i-1)×d,x1+i×d] (i=1,2,...,n), and its corresponding Y i Let Z be the number of residuals within this range, and for the i-th interval block, let the specific value corresponding to this interval block be Z. i for:

[0077] Z i =[(x1+(i-1)×d)+(x1+i×d)] / 2 (4-6)

[0078] Let there be a set of data (Z) i ,Y i (i = 1, 2, ..., n), each Z i Corresponding to a Y i This can be substituted into the function for fitting, and the fitting function is:

[0079]

[0080] Y in the formula max μ, σ, and μ' represent the peak value, mean, and standard deviation of the Gaussian curve, respectively. Taking the natural logarithm of both sides of equation (4-7), the equation becomes:

[0081]

[0082] make:

[0083]

[0084] Since there is more than one data point, the above formula is represented in matrix form:

[0085]

[0086] The above formula can be described in vector form as follows:

[0087] K = ZB (4-11)

[0088] The vectors in the formula are represented as follows:

[0089]

[0090] According to the least squares method, the parameter B to be solved can be expressed as:

[0091] B = (Z) T Z) -1 Z T K (4-13)

[0092] After finding b0, b1, and b2, substitute them into equation (4-9) to obtain the parameter values ​​of the Gaussian model. Multiple experiments have proven that the residual distribution of pulsars basically satisfies the Gaussian distribution model.

[0093] S5. After determining the model and its parameter values, the simulated residual values ​​can be obtained through simulation. However, it is first necessary to determine how many simulated residual values ​​to generate. Since the simulated residual values ​​are generally used to calculate the composite pulsar, the number of simulated residual values ​​is related to the number of residuals of each pulsar when composing the composite pulsar. Let the number of residual values ​​for each pulsar be n1, n2, n3, n4, ..., n k Then the number of simulated residuals, n, is:

[0094] n = min[n1, n2, n3, n4, ..., n k ]

[0095] S6. After determining the number of simulation residuals, perform simulation of the residual values. Figure 3 The image shows the simulation results for the PSR_B1855+09 pulsar. It clearly shows that the simulated residual values ​​closely match the original sample values, avoiding the problem of large fluctuations in the simulated data samples introduced by interpolation. After obtaining the residual values ​​from the simulation, it is necessary to consider how to add them to the sample observation data. The observation data contains not only timing residuals but also the observation times corresponding to the residuals. This is because subsequent calculations using σ... zWhen calculating the stability of the method, a cubic polynomial fitting is required between the two inputs. Therefore, the simulated residual values ​​must be identical to the actual residual values, and corresponding time values ​​must be set. Simultaneously, during the construction of the integrated pulsar time, the time span and observation time values ​​of the residual data from different pulsars must remain consistent. Therefore, the time values ​​set for the simulated residuals must correspond to the time values ​​of other pulsars.

[0096] This invention establishes a probability distribution model by statistically analyzing sample observation data, and then simulates and generates simulated residual values ​​that meet the requirements. The simulated residual values ​​generated by this method do not introduce new errors and can effectively match the real residual data. This effectively solves the problems of high-precision calculation and extraction required for pulsar timekeeping, making the stability calculation of pulsar timekeeping more accurate. Compared with traditional interpolation methods, this invention does not require segmenting the time series, considering the optimal segmentation interval, or searching for the best interpolation function in each interval. Furthermore, it eliminates the need for segmenting and interpolating the residual data, so the obtained simulated residual values ​​do not require filtering or removing the already limited actual observation values. After obtaining the statistical histogram, this invention, in order to ensure that the simulated residual values ​​meet the expected requirements, not only makes a rough guess about the model but also conducts a specific verification of the model, thereby making the simulated residual values ​​obtained more accurate. This solves the problem that insufficient observation samples lead to excessive fitting residual deviations and low stability of the synthesized pulsar timekeeping data during data processing. The conceptual challenge of this invention lies in abandoning the traditional interpolation approach and instead starting with the mathematical model of the residual data to obtain its distribution model. This not only simplifies the simulation process but also eliminates the additional errors introduced by the simulation.

[0097] Comparative example,

[0098] The timing residuals of a single pulsar can be obtained by fitting data using Tempo2 software, outputting two columns of data: the first column shows the time values, and the second column shows the timing residual values. After Tempo2 fitting is complete, the time values ​​and corresponding residual values ​​are sorted chronologically, and then adjacent time values ​​with a difference of 10 are identified. 4 The data is divided into intervals of seconds and longer, and these intervals are further segmented into several sub-intervals (the number of sub-intervals depends on the actual observation data of the pulsar). For each sub-interval, the residual value is set as the Y-axis and the time value as the X-axis. Different interpolation functions are used to obtain the simulated residual values ​​at the discontinuities. The simulated residual values ​​obtained from different interpolation functions are compared to find the most suitable interpolation function. Once the most suitable interpolation function has been found for each sub-interval, residual values ​​with large errors are removed. The remaining simulated residual values ​​are then used together with the actual observed residual values ​​in the calculation of the synthesized pulsar data.

[0099] Interpolation is essentially about adding a continuous function to discrete data so that the continuous curve passes through all the given discrete data points. Commonly used interpolation methods include linear interpolation, nearest neighbor interpolation, and cubic spline interpolation.

[0100] (1) Linear interpolation method

[0101] Let X1 and X2 be two discrete data points, with corresponding values ​​Y1 and Y2. Then, the function is established as follows:

[0102] Y=Y1+(Y2+Y1)×(X-X1) / (X2-X1) (2-1)

[0103] When X0∈[X1,X2], if you want to find the corresponding value Y0 at X0, you only need to substitute X0 into the equation, and the obtained Y value is the value Y0.

[0104] (2) Interpolation between two nearest points

[0105] Let X1 and X2 be two discrete data points, where X1 < X2, and their corresponding Y-axis coordinates be Y1 and Y2, respectively. Let the desired data be Y0, and its corresponding data point be X0.

[0106] Let X1 < X0 < X2, then...

[0107] d1 = X0 - X1 (2-2)

[0108] d2=X2-X0 (2-3)

[0109] If d1 < d2, then Y0 = Y1; if d1 > d2, then Y0 = Y2; if d1 = d2, then Y0 = Y1 or Y0 = Y2.

[0110] (3) Spline interpolation

[0111] Suppose there are n+1 sample points in the interval [a,b], namely x1, x2, ..., xn. n The corresponding function values ​​are y1, y2, ..., y n Divide the x-axis into n sub-intervals, and in each sub-interval [x j ,x j+1 Within the range (j = 1, 2, ..., n-1), the interpolation function is a cubic function, and its expression is:

[0112] S j (x)=a j x 3 +b j x 2 +c j x+d j (2-4)

[0113] The following conditions must be met simultaneously:

[0114]

[0115] Since n unknowns require n equations to obtain a unique solution, boundary constraints need to be added to obtain a unique interpolation function value, as follows:

[0116]

[0117]

[0118] Only with the above two constraints can a unique value be obtained.

[0119] Figure 4 , Figure 5 , Figure 6 This is a comparison chart of the simulated residual values ​​and the actual residual values ​​obtained by three different interpolation methods for the PSR_B1855+09 pulsar. It can be seen that the Spline method has the largest interpolation error among the three methods, while the other two methods have relatively smaller errors, but they still have a very large error compared with the actual sample residual values.

[0120] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention are included within the scope of protection of the present invention.

Claims

1. A method for determining simulated residuals based on limited observation data, characterized in that, The method comprises the following steps: S1, the pulsar observation data is processed to obtain timing residuals, the timing residuals include time values and timing residual values; the timing residual values are arranged from large to small, and the corresponding time values are also arranged; the timing residual values are set as x1, x2,..., x n ; the corresponding time values are set as t1, t2,..., t n ; S2, [x1, x n The interval is divided into equally spaced small interval blocks. The number of sample timing residual values ​​in different interval ranges is counted, and the frequency distribution of sample timing residual values ​​in different intervals is plotted. S3, determining a probability statistical distribution model to which the frequency distribution belongs; S4, obtaining specific parameters of the probability statistical distribution model based on the timing residual values in the frequency distribution diagram and the number of the residual values, and establishing the probability statistical distribution model of the sample timing residual values; S5, determining the number of the simulation residual values; the number of the simulation residual values is the minimum value of the residual value numbers of the respective pulsars in the integrated pulsar time; S6, obtaining the simulation residual values through simulation according to the probability statistical distribution model, setting corresponding time point values for the simulation residual values, and adding the simulation residual values to the sample observation data, wherein the time point values of the simulation residual values correspond to the time point values of the pulsar observation data.

2. The method of claim 1, wherein, In the step S1, the pulsar observation data is processed by the Tempo2 software to obtain the timing residual values, and the processing accuracy is up to 1 ns.

3. The method of claim 1, wherein, The step S2 is expressed as: [x1+(i-1)×d,x1+i×d] (i=1,2,3...,n), and the statistical value is counted in each small interval block. If the timing residual value x j (j=1,2,...,n) is in the interval, the timing residual value x j The corresponding number y j is 1, otherwise 0, so the vertical coordinate value Y j (j=1,2,...,n) on the frequency distribution diagram is expressed as: The number of the sample timing residual values in different interval ranges is determined through formula (4-2).

4. The method of claim 1, wherein, The frequency distribution diagram of the sample timing residual values in different intervals is represented in the form of a histogram, wherein the abscissa represents the distribution interval of the timing residual values, each node is an end point of the interval, and the ordinate represents the number of the sample timing residual values in the corresponding interval.

5. The method of claim 1, wherein, The step S3 specifically comprises the following steps: observing the histogram to preliminarily determine which probability distribution model belongs to, the probability distribution model being any one of a Gaussian distribution, a uniform distribution, a geometric distribution or a Poisson distribution; different models correspond to different test methods, whether the preliminary determination is correct is determined through the test result, and finally the probability statistical distribution model corresponding to the sample timing residual data can be determined.

6. The method of claim 4, wherein, In the step S4, the abscissa of the histogram is actually a plurality of cell blocks, rather than specific values, so it is necessary to first change the cell blocks into specific values. Assuming that any one cell block is [x1+(i-1)×d, x1+i×d] (i=1, 2, …, n), the Y corresponding to the cell block is i the number of residual values in the range; for the i-th cell block, assuming that the specific value Z corresponding to the cell block is i is: Z i = [(xi + (i - 1) x d) + (xi + i x d)] / 2 (4-6) Let a set of data (Z i ,Y i )(i=1,2,...,n) be given, each Z i corresponding to a Y i , which can be brought into the function of the determined probability and statistics distribution model for fitting to obtain the specific parameters of the probability and statistics distribution model.

7. The method of claim 6, wherein, The probability statistical distribution model is a Gaussian distribution, and the function form of the Gaussian distribution is as follows: Each Z i Corresponding to one Y i , into the function of Gaussian distribution for fitting, the fitting function is: Y in the formula max μ, σ are the peak value, mean value and standard deviation of the Gaussian curve, respectively; and taking the natural logarithm of both sides of equation (4-7): Let: Since there is more than one data, the above formula is represented in a matrix form: The above formula is described in a vector form as follows: K=ZB (4-11) The vectors in the formula are respectively as follows: According to the solution of the least square method, the to-be-solved parameter B can be represented as follows: B = (Z T Z) -1 Z T K (4-13) After b0, b1 and b2 are obtained, they are substituted into formula (4-9), and the parameter value of the Gaussian model can be obtained.

8. The method of claim 1, wherein, In the step S5, the number of residual values of each pulsar is respectively n1, n2, n3, n4,..., n k The number n of the simulation residual values is: n = min[n1, n2, n3, n4,..., n k ].

Images