An optimization method for building a synthetic pulsar based on Wiener filtering

Abstract

The present invention provides an optimization method for establishing a comprehensive pulsar time based on Wiener filtering, wherein the reference time is TT (BIPM), and the model is fitted using existing data to obtain accurate model parameters; the reference time is changed to the monitored reference time, the timing residual before fitting is obtained and averaged within a set time length to obtain a number of residual points; all the mutual power spectra are calculated and adjusted according to the upper limit of the mutual power spectrum, and the two groups of residuals that determine the upper limit of the mutual power spectrum are filtered using the calculated mutual power spectra, and the remaining residuals are filtered using the arithmetic average of all the mutual power spectra; the correlation coefficients of each group of residuals before and after adding a tiny signal are calculated, and the normalized value of the non-correlation coefficient is used as the weight of each star; the filtering results of each star are weighted to finally obtain the comprehensive pulsar time. The present invention makes the fluctuation magnitude of the comprehensive pulsar time relative to the reference time in the correct range, and has good sensitivity to the change of the reference time.

Description

Technical Field

[0001] The present invention belongs to the field of pulsar time scale establishment, relates to the establishment of a pulsar timing model, and in particular to a method for establishing a comprehensive pulsar time scale. Background Art

[0002] Time is the most accurate basic physical quantity currently measured. Accurate time measurement is the foundation of many scientific and engineering fields. High-precision time standards currently rely mainly on International Atomic Time (TAI). TAI is calculated by the International Bureau of Weights and Measures BIPM based on data from the world's timekeeping laboratories and is released monthly. The Earth Time TT (BIPM) derived from TAI is currently the most accurate time scale based on atomic frequency standards. The International Bureau of Weights and Measures generally releases its clock difference data with the previous year's TAI at the beginning of the next year. For departments with high-precision time requirements, the reliance on foreign data releases reduces the reliability of high-precision time services. Therefore, it is necessary to find an independent high-precision time maintenance method.

[0003] Pulsars are compact stars formed at the end of stellar evolution, with rotation and radiation characteristics. Millisecond pulsars are the evolutionary products of pulsars after accelerated rotation, with a rotation period of milliseconds and very high rotation stability. Studies have shown that the rotation stability of some millisecond pulsars is better than atomic time on a 10-year scale. However, the pulse arrival signal obtained by direct observation cannot be used as a time signal, and it is also necessary to rely on a timing model, that is, a series of mathematical equations that describe the emission, propagation delay and time correction of pulsar signals to deduct the influence of various physical processes on signal propagation. In order to obtain accurate model parameters, it generally takes several years or even more than 10 years of observation. Once an accurate timing model is obtained, the pulsar time can be established.

[0004] As an astronomical signal, the pulses radiated by pulsars are affected by many factors during propagation, resulting in fluctuations in the arrival time of the pulses determined by this method, and making the established pulsar time noisy. The noise of pulsar time is mainly white noise and low-frequency noise, among which the level of white noise is currently higher than that of atomic clock signals and their derived time scales. In order to improve the stability of pulsar time, some methods need to be used to reduce the noise level, which is generally achieved by integrating the observation data of multiple stars. The basic principle of integrated pulsar time is based on the irrelevance of the noise generation process of different pulsars, and certain algorithms are used to eliminate irrelevant noise, so as to obtain a common signal under a certain time scale, and then the noise level is further reduced by weighted averaging. According to existing studies, the noise-reduced integrated pulsar time can successfully detect the fluctuation of TAI, reflecting its important application value.

[0005] At present, the basic methods for building synthetic pulsars include Wiener filtering and Bayesian methods. The calculation of Wiener filtering is relatively simple and easy to implement quickly. Its principle is to estimate the statistical characteristics of the detected signal by performing correlation operations on the timing residuals of different pulsars, and then use it to filter the residuals of each pulsar, and finally perform weighted averaging on the filtering results. The problem with Wiener filtering is that when several pulsars have high noise levels, the synthetic pulsar established will also have high noise, which will directly drown out the detected signal. Improving the signal energy estimation method and taking reasonable weights are the key to using Wiener filtering to build synthetic pulsars. Summary of the invention

[0006] In order to overcome the shortcomings of the prior art, the present invention provides an optimization method for establishing a synthetic pulsar time based on Wiener filtering. According to the estimation of signal energy by the cross-correlation between low-noise pulsar residuals, the magnitude of the cross-correlation (cross-power spectrum) of high-noise pulsar residuals is adjusted to ensure that the fluctuation level of each group of residual filtering results is no higher than the fluctuation level of the reference time signal. At the same time, the weights are determined according to the sensitivity of different pulsar residuals to tiny signal changes, so that the synthetic pulsar time can correctly follow the changes in the reference time.

[0007] The technical solution adopted by the present invention to solve the technical problem comprises the following steps:

[0008] 1) Using the initial model, change the reference time to TT (BIPM), and fit the model using the existing data until the root mean square of the residuals before and after fitting are equal to obtain accurate model parameters; change the reference time to the monitored reference time to obtain the timing residuals before fitting;

[0009] 2) The pre-fitting timing residuals obtained in step 1) are averaged over a continuous set time period to obtain a number of residual points;

[0010] 3) Fit the residual points obtained in step 2) and remove the quadratic terms;

[0011] 4) estimating the upper bound of the cross power spectrum; calculating and adjusting all cross power spectra according to the upper bound of the cross power spectrum so that all cross power spectra do not exceed the upper bound in the specified low frequency interval; after all cross power spectrum estimates are given, filtering the two groups of residuals that determine the upper bound of the cross power spectrum using the calculated cross power spectrum, and filtering the remaining residuals using the arithmetic average of all cross power spectra;

[0012] 5) Add a small signal to each group of residuals. The amplitude of the signal is taken as the magnitude of the signal fluctuation in the reference time. Calculate the correlation coefficient of each group of residuals before and after adding the small signal, and use the normalized value of the non-correlation coefficient as the weight of each star;

[0013] 6) Using the obtained weights, weight the filtering results of each star and finally obtain the comprehensive pulsar time.

[0014] The average residual point in step 2) is located at the midpoint of each interval, and a number of residual points distributed at equal intervals are obtained.

[0015] The step 4) uses an AR model or a periodogram method to estimate the upper bound of the cross-power spectrum.

[0016] In the step 4), if the upper bounds of the cross-power spectra given by the two methods are consistent, or the upper bounds of the cross-power spectra given by the two methods are inconsistent but there is no obvious low-frequency periodic component in the residual, that is, the spectrum shapes in the low-frequency interval are not all peaked, then the AR model method is used; if the estimates given by the two methods are inconsistent, and there are obvious low-frequency periodic components in the residual, that is, the spectrum shapes in the low-frequency interval are all peaked, then the periodogram method is used.

[0017] In the step 4), the low-frequency interval is determined according to the upper limit of the cross-power spectrum and the spectral shape of each calculated cross-power spectrum, and is taken as the 0 frequency point to the turning frequency, where the turning frequency refers to the frequency beyond which the upper limit of the cross-power spectrum becomes a white noise spectrum, or is taken as the half-high frequency interval of the main peak of the upper limit of the cross-power spectrum, so that the difference between the value of the finally calculated cross-power spectrum and the upper limit in the low-frequency interval is less than the set value.

[0018] The filtering formula in step 4) is: where modulus{s k} represents the residual x k The recovered signal k The modulus estimate of is a unit complex number, and the residual for each star is a constant.

[0019] The step 5) described above adds a unit period sine wave to each group of residuals.

[0020] The beneficial effects of the present invention are: reducing the high noise phenomenon of the integrated pulsar time in the traditional Wiener filtering method, making the fluctuation magnitude of the integrated pulsar time relative to the reference time within the correct range, and having good sensitivity to the change of the reference time. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 It is a schematic diagram of the impact of different energy noises on the residual phase;

[0022] Figure 2 This is the flow chart of the improved Wiener filter algorithm to establish a comprehensive pulsar;

[0023] Figure 3 is the filtering result obtained using simulated data. DETAILED DESCRIPTION

[0024] The present invention proposes a method for establishing a comprehensive pulsar time based on the existing frequency domain Wiener filtering method. The core idea of ​​this method is to estimate the signal energy based on the cross-correlation between low-noise pulsar residuals, adjust the cross-correlation (cross-power spectrum) level of high-noise pulsar residuals, and ensure that the fluctuation level of each group of residual filtering results is not higher than the fluctuation level of the reference time signal; at the same time, the weight is determined according to the sensitivity of different pulsar residuals to tiny signal changes, so that the comprehensive pulsar time can correctly follow the changes of the reference time.

[0025] The basic quantity to be determined for Wiener filtering in the frequency domain is the frequency response, which is expressed as

[0026]

[0027] Where x is the observed signal, s is the expected signal, and S represents the power spectrum. The Wiener filter method for building a synthetic pulsar is described as

[0028]

[0029] The subscript k represents the kth star, x k is the residual of the kth star, s k is the filtering result of the kth star, F is the Fourier transform operator, is the arithmetic mean of the cross-power spectra of all residuals. The final integrated pulsar time is expressed as

[0030]

[0031] μ k is the weight of the residual filtering result of the kth star in the weighting, which is generally taken as the normalized value of the inverse square of the root mean square of the filtering result.

[0032] Theoretical analysis and practical calculations show that the energy estimate of the reference time signal made by taking the arithmetic mean of the cross power spectrum will be too large due to the introduction of high-noise pulsars. Specifically, in the traditional Wiener filtering method, the estimate of the reference time signal power spectrum is generally expressed as

[0033]

[0034] N is the number of residuals.

[0035]

[0036]

[0037] According to the relationship between mutual spectrum and self spectrum

[0038]

[0039] It can be seen that the cross power spectrum estimate given by equation (4) is actually the upper bound of the cross power spectrum modulus value. When the noise is high, it will be seriously higher than the modulus of the real cross spectrum, resulting in If the cross power spectrum estimation caused by noise is too high, it must be adjusted.

[0040] On the other hand, according to the relationship between the residual spectrum phase and the signal spectrum phase (see Figure 1 ), we can see that when the noise ε is low, the phase of the residual x can roughly reflect the original phase of the signal; when the noise is very high, the residual phase is the superposition of the signal phase and a random change. At this time, even if the power spectrum is accurate, the signal recovered from the residual is not the desired reference time fluctuation, but the result after adding random disturbance.

[0041] Considering the above two points, we can rewrite equation (2) as

[0042]

[0043] where modulus{s k} represents the residual x k The recovered signal k The modulus estimate of is a unit complex number, and the residual for each star is a constant. k}Estimation method.

[0044] In order to accurately estimate the power spectrum of the clock error signal, the following assumption must be made based on the cross-correlation method: the noise level of at least two pulsar data is not higher than that of the reference time signal. In fact, if the reference time signal is very small compared to the residual noise, then equation (4) must be greater than the power spectrum of the reference time signal. Under the above assumption,

[0045] The real part of will converge to the power spectrum of the reference time signal, but it may be a negative number in actual calculations, so formula (4) is usually used instead. This is a conservative estimate that is larger than the upper limit of the reference time signal power spectrum. For other high-noise pulsars, the estimate given by formula (4) will be higher than the estimate given by the above two low-noise pulsars, so it is necessary to lower the calculated high-level power spectrum, at least to make it no higher than the value of formula (4) given by the above two low-noise pulsars. In this way, all cross-power spectra are limited to the upper limit of the reference time signal power spectrum. The modulus of the Fourier transform is calculated from the power spectrum, and then filtered by formula (8), and the fluctuations of the recovered signal are all limited to the fluctuation level of the reference time signal.

[0046] The above-mentioned restriction of the cross-power spectrum not exceeding the upper bound refers to a certain frequency range. Specifically, since the clock error fluctuations of concern are all low-frequency signals, the restriction on the power spectrum is also for a certain low-frequency interval. This low-frequency interval can be taken as from the 0 frequency point to a certain turning frequency. The power spectrum exceeding the turning frequency can be considered as a white noise spectrum; it can also be taken as the half-high frequency interval of the low-frequency peak of the power spectrum. In actual calculations, no matter which frequency interval is taken, it is necessary to ensure that the value of the reference time power spectrum calculated by the two sets of residuals that determine the upper bound of the cross-power spectrum is close to the upper bound in this interval, so as to ensure that the signal recovered by these two sets of residuals has the correct fluctuation magnitude. In specific problems, the case where the reference time power spectrum estimation result is a white noise spectrum is not ruled out. In order to prevent the high-frequency energy of the recovered signal from being too large, the power spectrum above the high-frequency point in the turning frequency or half-high frequency of the power spectrum can be set to zero, and then the signal modulus value is calculated.

[0047] For power spectrum estimation, in addition to the traditional periodogram method mentioned above, the AR model method can also be used. AR model power spectrum estimation is a high-resolution estimation method that can overcome the defect of large variance of the periodogram method and has advantages in power spectrum adjustment. The power spectrum estimation of the AR model method depends on the selection of the model. For the case where the residual contains periodic signals, the estimated power spectrum may have problems such as spectral peak offset and inaccurate spectral peak. Therefore, in practical applications, two methods should be used to estimate the upper bound of the cross-power spectrum and compare them. If the upper bounds of the cross-power spectrum given by the two methods are basically the same, or there are differences, but it shows that there is no obvious periodic component in the residual, the AR model can be used to estimate the power spectrum.

[0048] When the AR model is used to estimate the cross power spectrum, there are two estimation methods according to the frequency characteristics of the model / signal.

[0049] Specifically, the AR model

[0050]

[0051] The power spectrum of the random sequence described is

[0052]

[0053] Where A(z) is a polynomial in the complex variable z, p is the model order, and a j is the model coefficient (real number), X t represents the time series describing the residual, σ 2 is a white noise process ε t The variance of , i is the imaginary unit, ω is the circular frequency (the range is [0,π]). Model parameters p, a j With σ 2A variety of methods in AR model theory can be used for estimation, and the rationality of parameter estimation can be verified by white noise test. The first method of estimating the cross power spectrum using the AR model is similar to the periodogram method, that is, based on formula (7), the cross power spectrum is given by formula (11). In addition, the Wold coefficient of the stationary part of the AR process can also be used to directly calculate the cross-correlation function and obtain the cross power spectrum. The principle of which method to use is similar to that of the periodogram method, that is, it is necessary to ensure that the value of the cross power spectrum calculated by the two sets of residuals that determine the upper bound of the cross power spectrum is close to the upper bound in this interval. When filtering, the cross power spectrum calculated using the second method must eventually take the possible complex values ​​as real numbers. The specific method is to directly modulo the complex values ​​or take the real part of the complex values, and ensure the non-negativity of the power spectrum.

[0054] When using equation (8) to filter and build a synthetic pulsar, the residuals of different pulsars play different roles in essence. For high-noise pulsars, the cross-power spectrum estimate given by cross-correlation reflects the noise characteristics in the residuals, rather than the reference time itself. Only the cross-power spectrum given by the two low-noise pulsars that determine the upper bound of all cross-power spectra represents the characteristics of the reference time signal to the greatest extent, and the residual phases of these two stars are also closest to the reference time signal. Therefore, modulus {s k For the two pulsars that determine the upper bound of the cross-power spectrum, the best approach is to use the cross-power spectrum estimate given by these two stars; for the remaining high-noise pulsars, the calculated cross-power spectrum reflects the noise characteristics, and the residual phase is also a random variable under the superposition of the signal phase and random influences. Therefore, the mean of all cross-power spectra can be simply used for filtering to obtain a filtering result under an average power spectrum. Probabilistically speaking, when the number of pulsars involved in the establishment of a synthetic pulsar increases infinitely, the noise components in the signals recovered from the residuals of different pulsars will cancel each other out, allowing the true reference time signal to appear.

[0055] After filtering in the above manner, the weighting method mentioned in formula (3) can no longer be used when weighting the filtering results of each star, because the fluctuation level of the filtering results can no longer reflect the noise level of the original residual. The present invention proposes a new weighting method: weighting according to the sensitivity of each star to the detected signal. Specifically, a small signal is added to the residual of each star (without loss of generality, the small signal can be a unit period sine wave of a certain amplitude, and the amplitude of the sine wave is taken as the magnitude of the signal fluctuation at the reference time) to produce a contrast residual, and the correlation coefficient between the original residual and the contrast residual of each star is calculated. For pulsars with large correlation coefficients, lower weights should naturally be used, because high correlation indicates that the noise level of the star is high and is insensitive to changes in the reference time. Based on this, the "non-correlation coefficient" is defined

[0056] w=1-ρ, (12)

[0057] where ρ is the autocorrelation coefficient

[0058]

[0059] The normalized w is the weight μ in formula (3).

[0060] The method for calculating the cross power spectrum using the Wold coefficient under the AR model theory is as follows:

[0061] According to the AR model theory, the random process X described by equation (10) t The stationary part of

[0062]

[0063] In the formula For X t The Wold coefficient can be calculated by the following formula

[0064]

[0065] And has the nature

[0066]

[0067] therefore Absolutely possible.

[0068] Suppose there are two random processes X t , Y t The Wold coefficient and white noise are {ψ j} and {ε s},{ε t ′}, it can be proved

[0069]

[0070] In AR model theory, the white noise of two random processes must be uncorrelated when si≠tj. In practical applications, when the correlation function of two random processes is not 0, the white noise must be correlated when si=tj, otherwise it will lead to E(X s Y t )=0. Definition

[0071] σ c 2 =E(ε t ε t ′), (18)

[0072] Let ts = τ, and we get

[0073]

[0074] Since the Wold coefficient converges to 0 with negative exponential order, t , Y t Under stable conditions, equation (19) can obtain high accuracy by performing only a small number of finite summation operations. According to the Wiener-Khinchin theorem, the Fourier transform of equation (19) is the cross power spectrum.

[0075] (19) where σ c 2 is an unknown quantity. In actual calculation, we only need to adjust its size according to the upper bound of the cross power spectrum to obtain an estimate of the cross power spectrum.

[0076] In summary, the technical solution of the present invention includes the following steps:

[0077] 1) Using the initial model, change the reference time to TT (BIPM), and fit the model using the existing data until the root mean square (rms) of the residuals before and after fitting are equal to obtain accurate model parameters; change the reference time to the monitored reference time to obtain the timing residuals before fitting;

[0078] 2) The pre-fitting timing residuals obtained in step 1) are averaged over a continuous set time period, with the average residual point located at the midpoint of each interval, to obtain a number of residual points distributed at equal intervals;

[0079] 3) Fit the residual points obtained in step 2) and remove the quadratic terms;

[0080] 4) Use the AR model or the periodogram method to estimate the upper bound of the cross-power spectrum, which is determined by the two groups of residuals with the lowest noise level; calculate and adjust all cross-power spectra based on the upper bound of the cross-power spectrum so that all cross-power spectra do not exceed the upper bound within a specific low-frequency interval; the low-frequency interval is determined based on the upper bound of the cross-power spectrum and the spectral shape of each calculated cross-power spectrum, and can be taken as the 0 frequency point to the turning frequency (the turning frequency refers to the frequency beyond which the upper bound of the cross-power spectrum can be regarded as a white noise spectrum), or taken as the half-high frequency interval of the main peak of the upper bound of the cross-power spectrum. In specific problems, the final calculated cross-power spectrum (especially the cross-power spectrum calculated by the two groups of residuals that determine the upper bound of the cross-power spectrum) should be within the low-frequency interval. The value is similar to the upper bound. In the application, the judgment criteria for using the AR model or the periodogram method to estimate the power spectrum are as follows: if the upper bounds of the cross power spectrum given by the two methods are consistent, or the upper bounds of the cross power spectrum given by the two methods are inconsistent but there is no obvious periodic component in the residual, then the AR model method is used; if the estimates given by the two methods are inconsistent, and there is an obvious low-frequency periodic component in the residual (the spectrum in the low-frequency range is all peak-shaped), then the periodogram method is used; after all the cross power spectrum estimates are given, the two groups of residuals that determine the upper bound of the cross power spectrum are filtered using the calculated cross power spectrum, and the remaining residuals are filtered using the arithmetic average of all cross power spectra. The filtering formula is formula (8);

[0081] 5) Add a small signal (such as a unit period sine wave) to each group of residuals. The amplitude of the signal is taken as the magnitude of the signal fluctuation in the reference time. Calculate the correlation coefficient of each group of residuals before and after adding the small signal (Formula (13)), and use the normalized value of the "non-correlation coefficient" (Formula (12)) as the weight of each star;

[0082] 6) Using the obtained weights, weight the filtering results of each star and finally obtain the comprehensive pulsar time.

[0083] The present invention is further described below in conjunction with the accompanying drawings and embodiments. The present invention includes but is not limited to the following embodiments.

[0084] 1) Using the pulsar timing analysis software tempo2, refer to TT (BIPM15) to simulate 6 sets of data located between MJD50188 and 55604, with an observation interval of 10 days, and the white noise levels are 100ns (two sets), 200ns, 300ns, and 1μs (two sets). The first set of data with a white noise of 100ns is added with a transition frequency f c =0.067yr -1 , spectral index α = 2, magnitude amp = 1.86 × 10 -40 yr 3 Low frequency noise; white noise of 200ns is added to the data f c =0.07yr -1, α=4,amp=3×10 -25 yr 3 Low frequency noise; white noise of 300ns is added to the data f c =0.06yr -1 , α=2,amp=5×10 -25 yr 3 Change the reference time to TT(TAI) and obtain and output the timing residual before fitting of the reference TT(TAI).

[0085] 2) The residuals obtained in step 1) are averaged over consecutive 30-day intervals, and the average residual point is taken at the midpoint of each interval, resulting in 180 equally spaced residual points located between MJD50206 and 55576.

[0086] 3) Fit the residuals obtained in step 2) and remove the second-order terms, which is equivalent to obtaining the fitted residuals of the fitted rotation frequency and its first-order derivative. This removes the main non-stationary components in the residuals.

[0087] 4) The upper bound of the cross-power spectrum of the residual with a white noise level of 100ns is calculated using the periodogram method and the AR model. It is found that the upper bound of the cross-power spectrum given by the two methods reflects that the residual contains obvious low-frequency periodic components, and there are differences in magnitude and peak drift. Therefore, the periodogram method is used to estimate the upper bound of the cross-power spectrum and perform cross-power spectrum calculation and adjustment. The upper bound of the cross-power spectrum is determined by the two groups of residuals with the lowest noise level (two groups of residuals with a white noise level of 100ns), and the half-high frequency interval of the spectrum peak is taken as the frequency interval limited by the power spectrum. After calculating and adjusting the cross-power spectrum, each group of residuals is filtered according to the method proposed by the present invention.

[0088] 5) A unit period sine wave with an amplitude of 200 ns is added to each set of simulation residuals, the correlation coefficients of the six sets of residuals before and after the sine wave is added are calculated, and the weight of each star is determined. The results are listed in Table 1.

[0089] Table 1 Correlation coefficient of residuals before and after adding sine wave to simulated data and determined weights

[0090]

[0091] 6) Using the obtained weights, the filtering results of each star are weighted, and the final comparison between the obtained integrated pulsar time and TT(BIPM15)-TT(TAI) is shown in Figure 3 .

Claims

1. An optimization method for building a synthetic pulsar based on Wiener filtering, It is characterized in that The following steps are involved: 1) Using the initial model, change the reference time to TT (BIPM), and fit the model using the existing data until the root mean square of the residuals before and after fitting are equal to obtain accurate model parameters; change the reference time to the monitored reference time to obtain the timing residuals before fitting; 2) The pre-fitting timing residuals obtained in step 1) are averaged over a continuous set time period to obtain a number of residual points; the average residual point is located at the midpoint of each interval to obtain a number of residual points distributed at equal intervals; 3) Fit the residual points obtained in step 2) and remove the quadratic terms; 4) Estimate the upper bound of the cross power spectrum; calculate and adjust all cross power spectra according to the upper bound of the cross power spectrum so that all cross power spectra The upper bound is not exceeded within the specified low-frequency interval; after all cross-power spectrum estimates are given, the two groups of residuals that determine the upper bound of the cross-power spectrum are filtered using the cross-power spectra calculated therefrom, and the remaining residuals are filtered using the arithmetic average of all cross-power spectra; The filtering formula is: , where modulus{s k } represents the residual x k The recovered signal k The modulus estimate of is a unit complex number, and the residual for each star is a constant; F is the Fourier transform operator; 5) Add a small signal to each group of residuals, and take the amplitude of the signal as the magnitude estimate of the signal fluctuation at the reference time. Calculate the correlation coefficient of each group of residuals before and after adding the tiny signal, and use the normalized value of the non-correlation coefficient as the weight of each star; 6) Using the obtained weights, weight the filtering results of each star and finally obtain the comprehensive pulsar time.

2. The optimization method for establishing a synthetic pulsar based on Wiener filtering according to claim 1, It is characterized in that The step 4) uses an AR model or a periodogram method to estimate the upper bound of the cross-power spectrum.

3. The optimization method for establishing a synthetic pulsar based on Wiener filtering according to claim 1, It is characterized in that In the step 4), if the upper bounds of the cross-power spectra given by the two methods are consistent, or the upper bounds of the cross-power spectra given by the two methods are inconsistent but there is no obvious low-frequency periodic component in the residual, that is, the spectrum shapes in the low-frequency interval are not all peaked, then the AR model method is used; if the estimates given by the two methods are inconsistent, and there are obvious low-frequency periodic components in the residual, that is, the spectrum shapes in the low-frequency interval are all peaked, then the periodogram method is used.

4. The optimization method for establishing a synthetic pulsar based on Wiener filtering according to claim 1, It is characterized in that In the step 4), the low-frequency interval is determined according to the upper limit of the cross-power spectrum and the spectral shape of each calculated cross-power spectrum, and is taken as the 0 frequency point to the turning frequency, where the turning frequency refers to the frequency beyond which the upper limit of the cross-power spectrum becomes a white noise spectrum, or is taken as the half-high frequency interval of the main peak of the upper limit of the cross-power spectrum, so that the difference between the value of the finally calculated cross-power spectrum and the upper limit in the low-frequency interval is less than the set value.

5. The optimization method for establishing a synthetic pulsar based on Wiener filtering according to claim 1, It is characterized in that The step 5) described above adds a unit period sine wave to each group of residuals.

Images