GRACE-FO inter-satellite ranging data processing method based on improved CRN filtering

Abstract

The invention relates to the field of inter-satellite ranging, in particular to a GRACE-FO inter-satellite ranging data processing method based on improved CRN filtering, and the method comprises the steps: constructing an LRI1A data preprocessing process in a low-order gravity satellite detector, reading LRI1A-M and LRI1A-T phase measurement original data, and carrying out the phase unwrapping and time mark correction processing; calculating a gain fluctuation level at an orbit frequency, and selecting a low-pass filter design parameter index; constructing an FIR digital filter based on a window function method, calculating time domain tap coefficients of a CRN filter and an improved filter, and carrying out normalization processing at an orbital frequency; according to an LRI biased distance calculation principle, calculating an inter-satellite biased distance and an inter-satellite distance, a rate and an acceleration after low-pass filtering; zero filling is carried out on the offset distance between satellites, the amplitude-frequency response and phase-frequency response characteristics of the offset distance between the satellites are calculated, and the improved filtering algorithm has great advantages in the aspects of improving passband flatness, enhancing high-frequency noise suppression and shortening the order of the filter. According to the invention, the precision and stability of gravitational field inversion are significantly improved.

Description

Technical Field

[0001] The present invention relates to the field of inter-satellite ranging, and particularly to a method for processing inter-satellite ranging data of GRACE-FO based on an improved CRN filter. Background Art

[0002] As a follow-up project of the GRACE (Gravity Recovery And Climate Experiment) mission, GRACE-FO (GRACE Follow-On) not only carries a K-band ranging system (KBR) but also realizes laser ranging interferometry (LRI) for the first time. This system achieves a nanometer-level measurement accuracy of the inter-satellite distance rate of change through dual-frequency laser phase-locking technology, significantly improving the order and spatial resolution of the time-varying gravity field model.

[0003] In the Level-1 data processing stage, the effective data frequency band for gravity field inversion is below 30 mHz. To meet the gravity field signal bandwidth requirements and reduce the computational complexity, a multi-stage downsampling strategy is implemented: the LRI1A data is downsampled from the original sampling rate of 9.664 Hz to 0.5 Hz after anti-aliasing filtering, and the inter-satellite distance, distance rate of change, and acceleration observables are synchronously extracted. At the same time, since the LRI1A data is interfered by various noises such as time scale error, phase jump, and ranging system error, which seriously affect the measurement accuracy and precision, fine data preprocessing is required to ensure that the extracted observed values have low noise and small amplitude distortion within the signal band.

[0004] The quadratic polynomial fitting method used in the early stage has significant defects in extracting inter-satellite rate and acceleration parameter estimation. This method has a gain of up to 104 times for inter-satellite acceleration noise, resulting in the low-frequency signal being submerged by high-frequency noise; and the amplitude distortion of the inter-satellite rate reaches -28 dB at 0.02 Hz, which cannot meet the accuracy requirements of high-order gravity field inversion. The linear fitting and triangular weighting methods also have problems of noise aliasing and gain ripple and cannot meet the task accuracy requirements.

[0005] The currently widely used CRN (N order self-Convolutions Rectangle) low-pass filtering algorithm is designed by multi-order self-convolution based on a rectangular window function, and its frequency domain response has a simple closed-form expression. This algorithm synchronously extracts the inter-satellite distance, inter-satellite rate, and inter-satellite acceleration through frequency domain differentiation, effectively suppressing high-frequency noise aliasing and controlling the passband gain ripple within Magnitude. However, this design method is not proposed as the best filter, and its sidelobe suppression ability shows an exponential decay trend with the increase of the number of self-convolutions. In addition, as the number of self-convolutions increases, it will also lead to an increase in the time span of the filter, thereby increasing the computational delay and storage requirements, making it inapplicable to real-time processing or satellite platforms with limited resources. Summary of the Invention

[0006] To solve the problems of high filter order, limited stopband suppression ability, and high computational complexity in the high-precision inversion of CRN filters, the present invention provides a GRACE-FO inter-satellite ranging data processing method based on improved CRN filtering. This method mainly includes: S1: Read the original LRI1A-M and LRI1A-T phase measurement data and perform phase unwrapping and time scale correction processing; S2: Calculate the gain fluctuation level at the orbital frequency and select the design parameter indicators of the low-pass filter; S3: According to the design parameter indicators of the low-pass filter, construct an FIR digital filter based on the window function method, calculate the time-domain tap coefficients of the improved filter, and perform normalization processing at the orbital frequency; S4: According to the principle of LRI biased distance calculation, calculate the inter-satellite biased distance, the inter-satellite distance, speed, and acceleration after low-pass filtering through the data obtained in step S1; S5: Zero-fill the inter-satellite biased distance and calculate its amplitude-frequency response, phase-frequency response characteristics, and the impact on the gain fluctuation at the orbital frequency.

[0007] A computer-readable storage medium stores a computer program, which realizes the steps of the above method when executed by a processor.

[0008] The beneficial effects brought by the technical solution provided by the present invention are as follows: By introducing the Kaiser window to replace the self-convolution rectangular window and optimizing the window function parameters, β independently control the stopband attenuation and passband gain levels, and realize flexible regulation of the frequency-domain response characteristics of the filter. Using the Kaiser window to replace the self-convolution rectangular window and combining adjustable parameters β dynamically balance the main lobe width and stopband attenuation, break through the limitation that the stopband suppression ability of CRN filtering decays exponentially with the number of convolutions, and at the same time reduce the passband gain ripple to 1.11×10 -16 magnitude.

[0009] Using the zero-padding method, the original data is extended and its amplitude-frequency response characteristics are calculated to ensure that the discrete Fourier transform can cover all frequency components, improve the spectral resolution, and avoid high-frequency noise aliasing. Using the orbital frequency normalization gain control method, the filter coefficients are gain-normalized at the orbital frequency of 0.176 mHz to ensure that the gain within the passband is 1, avoid unexpected amplification or attenuation of the signal amplitude, and ensure that the gain fluctuation level at the orbit is small enough.

[0010] The present invention verifies the feasibility and scientificity of the CRN low-pass filter proposed by JPL as a concept demonstration in the GRACE-FO satellite mission through LRI1A on-orbit data. While effectively suppressing high-frequency noise, using the principle of frequency-domain differentiation, the inter-satellite distance and its first-order (inter-satellite rate) and second-order derivatives (inter-satellite acceleration) information are synchronously extracted in the frequency domain, that is, low-pass filtering and the extraction of inter-satellite rate and inter-satellite acceleration are realized, avoiding the amplification of high-frequency noise caused by complex differentiation in the time domain, ensuring that the amplitude gain error of the gravity harmonic components is minimized, and solving the inherent limitations of the CRN filter. Through innovative window function design and dynamic parameter optimization, the data processing accuracy and efficiency are significantly improved. This method provides the best performance in terms of gain fluctuation and filter time span. While maintaining the computational efficiency, it provides a more refined input signal for the subsequent identification and removal of phase jumps, provides technical support for generating more accurate LRI1B data products, and significantly improves the accuracy and stability of gravity field inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] The present invention will be further described below in conjunction with the drawings and embodiments. In the drawings: Figure 1 is a flowchart of a method for processing GRACE-FO inter-satellite ranging data based on an improved CRN filter in an embodiment of the present invention; Figure 2 is a schematic diagram of the LRI frequency change in an embodiment of the present invention; Figure 3 is a schematic diagram of the phase delay of the CRN filter in an embodiment of the present invention; Figure 4 is a schematic diagram of the phase delay of the improved filter in an embodiment of the present invention; Figure 5(a) is a schematic diagram of the amplitude-frequency response under the same stopband attenuation in an embodiment of the present invention; Figure 5(b) is a schematic diagram of the partial enlargement of the amplitude-frequency response under the same stopband attenuation in an embodiment of the present invention; Figure 6 is a schematic diagram of the influence of the CRN filter on the LRI ranging signal in an embodiment of the present invention; Figure 7 is a schematic diagram of the influence of the improved filter on the LRI ranging signal under the same stopband attenuation in an embodiment of the present invention; FIG. 8(a) is a schematic diagram of the amplitude-frequency response at the minimum filter order in an embodiment of the present invention; FIG. 8(b) is a schematic diagram of a partially enlarged amplitude-frequency response at the minimum filter order in an embodiment of the present invention; Figure 9 is a schematic diagram of the influence of the improved filter on the LRI ranging signal at the minimum filter order in an embodiment of the present invention; FIG. 10(a) is a schematic diagram of the amplitude-frequency response at the minimum passband ripple in an embodiment of the present invention; FIG. 10(b) is a schematic diagram of a partially enlarged amplitude-frequency response at the minimum passband ripple in an embodiment of the present invention; Figure 11 is a schematic diagram of the influence of the improved filter on the LRI ranging signal at the minimum passband ripple in an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0012] In order to have a clearer understanding of the technical features, objectives, and effects of the present invention, the specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0013] Embodiment 1 Please refer to Figure 1 , Figure 1 is a flowchart of a method for processing GRACE-FO inter-satellite ranging data based on an improved CRN filter in an embodiment of the present invention, specifically including: S1. Construct a data preprocessing flow for LRI1A in a low-order gravity satellite detector, read the original LRI1A-M and LRI1A-T phase measurement data, and perform phase unwrapping and time scale correction processing.

[0014] S2. Calculate the gain fluctuation level requirement at the orbital frequency, and select appropriate low-pass filter design parameter indicators.

[0015] In order to extract the inter-satellite biased distance, velocity, and acceleration at a low rate of 0.5 Hz from the LRI1A data at a high frequency of 9.664 Hz, it is necessary to focus on the high-frequency noise processing before downsampling to avoid aliasing effects. Since the low-frequency gravity harmonic signal in the LRI1A data has a relatively high amplitude near the center frequency of the passband, its gain fluctuation level requirement is extremely strict. Therefore, while implementing the low-pass filtering function, it is necessary to meet the accurate extraction requirements of the inter-satellite distance, rate, and acceleration parameters. By calculating the gain fluctuation level at the orbital frequency, it is ensured that the filter design can balance the gain stability in the passband and the stopband suppression characteristics.

[0016] At the orbital frequency of 0.176 mHz, the amplitude of the gravity harmonic is approximately 3 kmG, the allowable gravity coefficient error is 0.01 cmG (centimeter gal), and the gain fluctuation at this point needs to be less than -150 dB, calculated as follows: (1) (2) For the remaining higher-order harmonics (such as ), the sum of the squares of all the gravitational field coefficients related to the given order value The relationship with the order transformation is expressed as: (3) where n represents the order, the order for spherical harmonic expansion, and the higher-order coefficients usually correspond to finer spatial resolution. represents the sum of the squares of the gravitational harmonic coefficients of the nth order. It can be seen from formula (3) that the sum of the squares of the gravitational field coefficients decays with the order.

[0017] Formula (3) is an empirical formula that describes the trend of the amplitude of the gravitational field coefficients decaying with the increase in order. It provides a basis for the subsequent derivation of the filter gain ripple requirements, that is, the amplitude differences of the gravitational signals of different orders are very large, indicating that the low-order signals have higher requirements for the accuracy of the filter.

[0018] Through the geometric relationship of spherical harmonic expansion, the geoid error scale factor is introduced to convert the dimensionless into the equivalent geoid error: (4) Then the remaining higher-order gain ripple errors can be relaxed to the following level: (5) In the actual filter design, the gain fluctuations of the gravitational harmonics at the orbital frequency must strictly satisfy .

[0019] S3. Based on the design parameter indicators of the low-pass filter, construct an FIR digital filter using the window function method, calculate the time-domain tap coefficients of the improved filter, and perform normalization processing at the orbital frequency.

[0020] The design of the filter parameters is shown in Table 1: Table 1 Low-pass filter parameter indicators for GRACE-FO LRI1A data

[0021] It can be seen from Table 1 that the appropriate design parameter indicators of the low-pass filter are as follows: Cutoff frequency: According to the output sampling rate of 0.5 Hz, the cutoff frequency should be 0.25 Hz; Gain fluctuation level at the orbital frequency: Based on the fact that at the orbital frequency of 0.176 mHz, the amplitude of the gravity harmonic is about 3 kmG and the allowable error of the gravity coefficient is 0.01 cmG, it is obtained through the gain fluctuation calculation formula. Level magnitude; Stopband attenuation: For comparison with CRN, set its stopband attenuation level based on that of CRN, and require that the attenuation of the improved filter in the stopband is not lower than this value. Filter length: Since the minimum filter length is to be obtained, the 747 length of CRN is used as the critical value for discussion.

[0022] S3.1. Calculate the time-domain tap coefficients of the CRN filter The CRN filter designs the N-fold self-convolution of the rectangular window function in the time domain as an ideal low-pass filter. According to the parameters in Table 2, the time-domain tap coefficients of the CRN filter can be calculated. c Table 2 Design specifications of CRN low-pass filter parameters

[0023] Let the window length be N, and its time-domain window function can be expressed as:

[0024] (6) (6) Further obtain its Fourier transform: (7) At this time, the expression of formula (7) corresponds to the Sinc function form in the frequency domain. Therefore, the frequency response of the CRN filter can be expressed as: (8) Among them, B is the filter bandwidth, which is half of the downsampling rate, 0.25 Hz. Perform the discrete inverse Fourier transform on H K to obtain the time-domain tap coefficients (weight function) of the CRN filter: (9) (10) Among them, represents the time-domain tap coefficients of the CRN filter, represents the normalization factor, H K represents the frequency response of the CRN filter, represents the filter length, n represents the position index of the tap coefficient in the time domain, represents the orbital frequency, represents the original data sampling frequency.

[0025] In formula (10), through the normalization factor Pair Normalize at the orbital frequency to ensure that the gain within the passband is 1 and avoid amplification or attenuation of the signal amplitude during the filtering process.

[0026] S3.2. Calculate the time-domain tap coefficients of the improved filter. Introduce the Kaiser window to replace the self-convolution rectangular window. According to the empirical formula, its time-domain expression is: (11) In the formula, denotes the zero-order Bessel function, which can be used to adjust the shape and index of the window function. The β parameter is used to control the sidelobe attenuation. Different β values can obtain different frequency-domain characteristics. β is calculated by the following formula: (12) Among them, denotes the stopband attenuation. Different from the CRN filter that generates the frequency-domain response through multiple self-convolutions and then performs the inverse discrete Fourier transform to obtain the time-domain filter coefficients, the tap coefficients of the filter improved based on the Kaiser window are directly calculated from its time-domain expression. Substitute into formula (13) to obtain the time-domain tap coefficients of the improved filter : (13) Among them, denotes the time-domain expression of the Kaiser window, denotes the cut-off frequency of the improved filter, denotes the original data sampling frequency, and n denotes the position index of the tap coefficient in the time domain.

[0027] Furthermore, normalize it. Define the discrete cosine reference signal generated at the orbital frequency as: (14) Obtain its normalization factor at the orbital frequency as: (15) (16) In the formula, denotes the normalized tap coefficients of the improved filter, G K denotes the gain normalization factor at the orbital frequency and is used to calibrate the response of the filter at this frequency. denotes the filter length, denotes the orbital frequency.

[0028] Calculate the equivalent gain of the filter coefficients at the orbital frequency through time-domain dot product, and scale the filter coefficients according to the gain calibration value to obtain the time-domain tap coefficients of the improved filter.

[0029] S4. According to the LRI biased distance calculation principle, calculate the inter-satellite biased distance ρ, and the inter-satellite distance, rate, and acceleration after low-pass filtering, using the data obtained in step S1. The signal after low-pass filtering can be expressed in the following convolution form: (17) (18) (19) Among them, the length of the filter coefficient is , represents the output index after the convolution of the filter weight function and the original biased distance, and its value range is , N ρ represents the length of the original biased distance data. n is the filter sampling index, represents the time-domain tap coefficient of the CRN filter, represents the first derivative of the time-domain tap coefficient of the CRN filter, represents the second derivative of the time-domain tap coefficient of the CRN filter, represents the original biased distance sequence and the current output index i related to the i-n th data point corresponding to the biased distance value. Each time a new output value is calculated, points need to be input, and the value will increase at a set time interval. , , respectively represent the inter-satellite distance, inter-satellite rate, and inter-satellite acceleration after low-pass filtering.

[0030] In formula (18), the calculation of the first derivative of the time-domain tap coefficient of the CRN filter is expressed as: (20) In formula (19), the calculation of the second derivative of the time-domain tap coefficient of the CRN filter is expressed as: (21) Among them, GRACE-FO adopts a dual-satellite formation configuration, and LRI is carried on both satellites. Satellite C is configured as the main satellite (hereinafter referred to as "M satellite"), and satellite D serves as the responder satellite (hereinafter referred to as "T satellite"). After the laser emitted by the laser on the M satellite is frequency-stabilized, it passes through the fast steering mirror and two laser beam splitters in sequence, and finally a retroparallel beam is formed by the triple reflector and directed to the T satellite. Due to the relative motion between the two satellites, the laser signal received by the T satellite has a frequency shift due to the Doppler effect. A small part of the local reference light on the T satellite interferes with the incident light to complete the tracking of the phase and frequency; at the same time, most of the local reference light passes through phase replication, and the phase and frequency information of the incident light are re-emitted back to the M satellite by the triple reflector. During this process, the returned beam undergoes a Doppler frequency shift again. On the M satellite, the returned laser signal interferes with the initial local reference light. At this time, the beat frequency of the interference signal directly reflects the relative velocity information between the two satellites. By integrating the frequency difference, the biased distance ρ between the satellites can be accurately calculated.

[0031] S4.1 The calculation principle of the biased distance between satellites is as Figure 2 shown in the schematic diagram of the LRI frequency change: Define the absolute laser frequency emitted by the M satellite through the frequency stabilizer as the laser (including laser frequency noise), with an offset frequency of about 10 MHz and a laser wavelength of . After the laser is split by the beam splitter, it forms the local reference laser LO and the transmitted laser TX. Then, the Doppler frequency shift caused by the relative velocity between the two satellites is expressed as: (22) (23) Among them,

[0032] The received frequency reaching the T satellite after the Doppler effect is: (24) (25) Among them, represents the received frequency of the laser reaching the T satellite, represents the received frequency of the laser reaching the M satellite.

[0033] The heterodyne frequency on the two satellites is obtained: (26) (27) Among them, represents the heterodyne frequency of the T satellite, Represents the heterodyne frequency of star M. The beat frequency signal of star M contains ranging information, which is tracked and recorded by the Laser Ranging Processor (LRP) for its phase change as the main observable. The phase measurement data of star T and star M are integrated: (28) (29) Among them, Represents the phase measurement data of star T, Represents the phase measurement data of star M.

[0034] S4.2 The binary star phase information calculated by formulas (28) and (29) is the original phase information of LRI1A - M and LRI1A - T extracted in step S1. After phase unwrapping and time scale correction, the inter - satellite distance signal can be further obtained, expressed as: (30) (31) Among them, Represents the inter - satellite distance at the initial time t0, Represents the change in distance from the initial time t0 to the current time t, Represents the biased inter - satellite distance, Represents the absolute laser frequency emitted by star M after passing through the frequency stabilizer, Represents the offset frequency.

[0035] S5. Zero - pad the biased inter - satellite distance and calculate its amplitude - frequency response, phase - frequency response characteristics, and the impact on the gain fluctuation at the orbital frequency: S5.1 Calculate the amplitude - frequency response: Taking the LRI1A - M star data on March 18, 2019 as an example, there are 368729 data samples on star M. Expand the filter to the original data length, construct a vector of length N, and calculate the amplitude - frequency response as follows: (32) Represents the frequency - domain response of the improved filter. The amplitude - frequency response of the CRN filter is also obtained by solving formula (32).

[0036] S5.2 Further calculate the phase delay: After extracting the real and imaginary parts of the complex signal in formula (30), use the arctangent function to obtain the phase of the frequency response, and then perform phase unwrapping to convert the phase into a continuous and unbroken phase to avoid Phase delay calculation error: (33) And calculate the phase delay corresponding to each frequency point: (34) Among them, represents the segmented characteristic of the phase angle presented in the frequency domain, and f represents the frequency variable in frequency domain analysis. The phase delay corresponding to each frequency point of the CRN filter is also obtained by solving formula (34).

[0037] Such as Figure 3 、 Figure 4 shown in the schematic diagram of the phase delay of different filters, represents the segmented characteristic of the phase angle presented in the frequency domain. In the frequency band less than 0.33 Hz, the constant phase delay is about 38.59 s, which can be used to compensate the influence of phase delay. The oscillation in the high-frequency band, and the signal values in the corresponding frequency band have been suppressed by the filter. Combining with the output frequency of 0.5 Hz of the LRI1B data product, such residual phase perturbations have no substantial impact on the subsequent scientific data products.

[0038] S6. Compare the advantages of the optimized filter and the CRN filter in terms of filter length, passband gain ripple, and stopband attenuation.

[0039] S6.1. Same stopband attenuation The amplitude-frequency responses of the two filters calculated by formula (30) in step S5.1 are as Figure 5(a) 、 5(b) 、6 shown. For the frequency components higher than 0.25 Hz, the stopband attenuation is 170 dB. Improve the filter to match the CRN stopband attenuation. According to formula (12) in step S3.2, β = 13.56 can be calculated, and then the amplitude-frequency response of the improved filter can be obtained.

[0040] It can be seen from Fig. 5(a) and (b) that under the same stopband attenuation, the transition bandwidth of the improved filter is 0.06 Hz (0.31 Hz → 0.25 Hz), which is smaller than the transition bandwidth of the CRN filter of 0.09 H (0.34 Hz → 0.25 Hz). The CRN stopband attenuation finally drops to 350 dB. Although the lowest stopband attenuation of the improved filter is not as good as that of the CRN, it is already sufficient to cope with the actual noise level. That is, under the condition of the same stopband attenuation, the transition band of the improved filter is shorter and steeper, and the stopband attenuation rate is faster.

[0041] Furthermore, draw the influence of different filters on the filtered inter-satellite distance signal. Such as Figure 6 shown in the schematic diagram of the influence of the CRN filter on the LRI ranging signal. The amplitude spectral density diagram of the inter-satellite distance signal after low-pass filtering by the black straight line, corresponding to the left y-axis coordinate. The curve is used to show the influence of the filter on the orbital frequency, corresponding to the right y-axis coordinate. For the frequency band greater than 0.16 Hz, the curve rapidly increases to 1. When When it is close to 1, the M value is close to 0. This indicates that at this frequency, the gain of the system for the input signal is very small, that is, this frequency component is almost completely suppressed. The dashed curve is the result of element-by-element multiplication of the signal represented by the straight line, indicating that at the orbital frequency of 0.176 mHz, the influence of the CRN filter on the ranging signal is very small, only order of magnitude.

[0042] Similarly, as Figure 7 shown in the schematic diagram of the influence of the improved filter on the LRI ranging signal under the same stopband attenuation, the gain fluctuation level at the orbital frequency of the improved filter is order of magnitude. Neither of the two filters causes significant distortion to the gravity signal in the ranging data near the orbital frequency. Although the improved scheme has advantages in stopband characteristics and computational efficiency, its passband gain ripple level is two orders of magnitude higher than that of the CRN low-pass filter.

[0043] S6.2. Minimum-order filter with the same passband gain To reduce the computational load of the on-board processor, under the condition of meeting the gain at the orbital frequency and the stopband attenuation requirements, by gradually increasing the filter order and nesting loops, the minimum feasible filter length and its corresponding β value are globally searched. Iteratively optimized within a certain range, through two-layer nested loops, the outer layer traverses the value, and the inner layer traverses the value of β. Calculate the performance indicators of the filter under each combination and check whether the conditions are met.

[0044] The results are as Figure 8(a) , 8(b) shown in the amplitude-frequency response diagram at the minimum filter order. When the filter order is shortened to , the bandwidth of the improved filter is , and the stopband attenuation is 249.6 dB.

[0045] Similarly, as shown in Figure 10, the gain fluctuation level at the orbital frequency of the improved filter is the same as that of the CRN. Obviously, neither of the two filters causes significant distortion to the gravity signal in the ranging data near the orbital frequency.

[0046] Therefore, it can be concluded that the CRN filter requires 9 convolutions, while the improved filter can achieve the same performance at a lower order by adjusting β. At the same gain fluctuation level at the same orbital frequency as the CRN and the same transition bandwidth, the improved filter has a shorter filtering time, can reduce the storage and computational resource occupancy by 6.3%, and provides a feasible hardware implementation scheme for low-power on-board systems.

[0047] S6.3. Minimum passband gain fluctuation To meet the gain level at the orbital frequency Set the improved filter length to be the same as the CRN filter, which is 747. By traversing and adjusting the parameter β, the gain fluctuation level at the minimum orbital frequency is obtained. Since the stopband attenuation level needs to meet 170 dB ( ), the β value controls the balance between the main lobe width and the sidelobe attenuation. A larger β value will result in a narrower main lobe and higher sidelobe attenuation, but will increase the transition band bandwidth. Therefore, β is set in the range of [12, 23] and traversed and optimized iteratively with a step size of 0.1. At this time, the schematic diagram of the amplitude-frequency response under the minimum passband ripple is as shown in Figure 10(a) 、 10(b) the schematic diagram of the amplitude-frequency response under the minimum passband ripple.

[0048] The results show that when the filter order is the same as that of CRN, which is 747, and β = 22.8, as shown in Figures 5 and 6, the transition bandwidth of the improved filter at this time . It is only 0.005 Hz different from the CRN filter. The stopband attenuation is 241.9 dB, which is far better than that of the CRN filter.

[0049] It can be concluded that the experimental results are as shown in Figure 11 the influence of the minimum passband ripple improved filter on the LRI ranging signal. When , the passband gain fluctuation level presents a minimum value magnitude. Compared with Figure 6 the influence of CRN on the LRI signal, at the same order, the improved filter not only has better stopband attenuation effect, but also has improved ability to suppress passband gain fluctuation, meeting the requirements of the LRI1A data preprocessing for the filter.

[0050] From Figure 5(a) 、 5(b) 、 Figure 6 and Figure 7 it can be seen that at the same stopband attenuation level, compared with the CRN low-pass filter, the transition bandwidth of the improved filter is optimized from 0.09 Hz (0.34 Hz → 0.25 Hz) of CRN to 0.06 Hz (0.31 Hz → 0.25 Hz), the transition bandwidth is reduced by 20%, the stopband attenuation rate is faster, the suppression of out-of-band noise is faster, and the passband gain level can still meet the design requirements; from Figure 8(a) 、 8(b) and Figure 9 it can be seen that the improved filter can achieve the same filtering performance as CRN in a shorter time and is superior to CRN in stopband attenuation; from Figure 10(a) 、 10(b) and Figure 11 it can be seen that at the same filter length, by adjusting βValue, the gain ripple level at the orbital frequency is significantly reduced compared to CRN (6.66×10 -16 →1.11×10 -16 ), effectively reducing the gravity harmonic amplitude distortion, and the stopband suppression is still better than that of CRN filtering, meeting the requirement of less than -150 dB ripple for gain distortion in high-order (200×200) gravity field inversion. At the same gain ripple level as CRN, not only is the transition bandwidth the same, but also it has a shorter filter order, reducing the storage requirement and calculation delay, and improving the calculation efficiency. The present invention uses the dynamic parameter global optimization method to calculate the maximum allowable gain ripple threshold in the passband according to the gravity field inversion error, and optimizes the parameters through global scanning and iteration β , ensuring that the designed parameters meet the accuracy requirements of gravity field inversion.

[0051] Embodiment 2 A computer-readable storage medium stores a computer program, which when executed by a processor, implements the steps of the above method.

[0052] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for processing inter-satellite ranging data of GRACE-FO based on improved CRN filtering, characterized in that Including: S1: Read the original LRI1A - M and LRI1A - T phase measurement data and perform phase unwrapping and time scale correction processing; S2: Calculate the gain fluctuation level at the orbital frequency and select the design parameter indicators of the low - pass filter; S3: According to the design parameter indicators of the low - pass filter, construct a FIR digital filter based on the window function method, calculate the time - domain tap coefficients of the improved filter and perform normalization processing at the orbital frequency; S4: According to the LRI biased - distance calculation principle, calculate the inter - satellite biased distance, the inter - satellite distance, rate and acceleration after low - pass filtering through the data obtained in step S1; S5: Perform zero - padding on the inter - satellite biased distance, and calculate its amplitude - frequency response, phase - frequency response characteristics, and the impact on the gain fluctuation at the orbital frequency.

2. The method for processing inter-satellite ranging data of GRACE-FO based on improved CRN filtering according to claim 1, wherein, In S2, the calculation formula for gain fluctuation is: Among them, represents the high-order gain ripple error, represents the gravity coefficient error, represents the amplitude of the gravity harmonic at the orbital frequency, represents the gain fluctuation.

3. The method for processing the inter-satellite ranging data of GRACE-FO based on the improved CRN filtering according to claim 1, characterized in that In S2, the design parameter specifications of the low-pass filter include the cut-off frequency , the stopband attenuation , the gain fluctuation |1 - M| at the orbital frequency, and the filter length .

4. The method for processing the inter-satellite ranging data of GRACE-FO based on the improved CRN filtering according to claim 1, wherein, In S3, the calculation formula for the time - domain tap coefficients of the improved filter is: Among them, represents the time-domain expression of the Kaiser window, represents the cut-off frequency of the improved filter, represents the original data sampling frequency, and n represents the position index of the tap coefficient in the time domain; Perform normalization processing on the time - domain tap coefficients of the improved filter, and define the discrete cosine reference signal generated at the orbital frequency as: Obtain its normalization factor at the orbital frequency as follows: In the formula, represents the improved filter tap coefficient after normalization, G K represents the gain normalization factor at the orbital frequency for calibrating the response of the filter at this frequency; represents the filter length, n represents the position index of the tap coefficient in the time domain, represents the orbital frequency, represents the original data sampling frequency.

5. A GRACE-FO inter-satellite ranging data processing method based on an improved CRN filter as described in claim 1, characterized in that, In S4, the calculation formula for the inter - satellite biased distance is: Among them, represents the inter-satellite biased distance, represents the inter-satellite distance at the initial time t0, represents the distance change from the initial time t0 to the current time t, represents the absolute laser frequency emitted by satellite M through the frequency stabilizer, represents the offset frequency.

6. The method for processing the ranging data between GRACE-FO satellites based on the improved CRN filtering according to claim 1, wherein, In S4, the calculations of the inter - satellite distance, inter - satellite rate and inter - satellite acceleration are respectively: Among them, , , respectively represent the inter-satellite distance, inter-satellite rate of change, and inter-satellite acceleration after low-pass filtering, represents the output index after convolution of the filter weight function with the original biased distance, n is the filter sampling index, represents the time-domain tap coefficient of the filter, represents the first derivative of the time-domain tap coefficient of the filter, represents the second derivative of the time-domain tap coefficient of the filter, represents the original biased distance sequence and the current output index i related to the i - n biased distance value corresponding to the 7. The method for processing inter-satellite ranging data of GRACE-FO based on improved CRN filtering according to claim 6, wherein, Time-domain tap coefficients of the CRN filter The calculation formula is as follows: wherein, represents a normalization factor, H K represents the frequency response of the CRN filter, represents the filter length, n represents the position index of the tap coefficient in the time domain, represents the orbital frequency, represents the original data sampling frequency.

8. The method for processing inter-satellite ranging data of GRACE-FO based on improved CRN filtering according to claim 1, wherein, In S5, the calculation formula for the amplitude - frequency response is: Among them, represents the frequency-domain response of the filter, and f represents the frequency variable in frequency-domain analysis.

9. A GRACE-FO inter-satellite ranging data processing method based on an improved CRN filter as claimed in claim 1, characterized in that In S5, the phase-frequency response characteristic is phase delay , and the calculation formula is: Among them, represents the segmented characteristics presented by the phase angle in the frequency domain.

10. A computer-readable storage medium, characterized in that, There is a computer program which, when executed by a processor, implements the steps of the GRACE - FO inter - satellite ranging data processing method based on improved CRN filtering according to any one of claims 1 - 9.

Images