A method and system for time-domain aliasing error mitigation for gravity satellite inversion

Abstract

The application discloses a time-domain aliasing error weakening method and system for gravity satellite inversion, and relates to the field of gravity satellite data processing, which adopts a two-step gravity field inversion strategy: a closed-loop full-scale numerical simulation scene is constructed, a multi-day short-term gravity field is estimated as a stable background field through experience parameterization processing, residual autocovariance weight determination and Wiese type signal separation strategy; a low-order daily scale residual gravity field signal is estimated based on single-day observation data with the short-term gravity field as a constraint, and two-step daily scale gravity field results are obtained through combination; a regularization extension method based on singular value decomposition is introduced, a regularization parameter is determined through generalized cross-validation, and the stable low-order solution is extended to a middle-order coefficient in control. The application effectively weakens the time-domain aliasing error, improves the inversion stability of the low-order and middle-order coefficients of the daily scale gravity field, suppresses the noise amplification caused by weak observation parameters, significantly reduces the spatial noise at high latitudes and polar regions, and balances the time resolution and spatial expression ability of the gravity field.

Description

Technical Field

[0001] This invention relates to the field of gravity satellite data processing technology, and in particular to a method and system for reducing temporal aliasing errors in gravity satellite inversion. Background Technology

[0002] Gravity satellites can dynamically monitor Earth's time-varying gravity field by observing satellite orbital perturbations, inter-satellite distances, and inter-satellite distance variability, enabling dynamic monitoring of Earth system processes such as changes in global land water storage, ice sheet and glacier mass changes, and ocean mass migration. Existing gravity satellite missions, such as Gravity Recovery and Climate Experiment (GRACE) and Gravity Recovery and Climate Experiment-Follow-On (GRACEFollow-On), have been widely used in global mass change research. However, their standard time-varying gravity field products are mostly based on monthly scales, with low temporal resolution, making it difficult to effectively capture rapidly evolving surface mass migration phenomena such as floods, rapid glacial melting, and short-period ocean dynamic processes, thus failing to meet the application requirements for high temporal resolution gravity monitoring.

[0003] With the development of next-generation gravity satellite missions, high temporal resolution gravity field inversion on periodic and even daily scales has become an important research direction. However, existing technologies have multiple technical shortcomings in daily-scale inversion scenarios: First, the sparse global coverage of daily satellite orbits and insufficient observational geometric constraints lead to a significant increase in the ill-conditioned nature of the gravity field inversion equations. If daily observation data is used directly to estimate the full-order spherical harmonic coefficients, the high proportion of weak observation parameters will cause severe noise amplification, resulting in poor stability of the inversion results and significant spatial striping errors, especially in high-latitude and polar regions where the noise amplification effect is more prominent.

[0004] Secondly, the sources of temporal aliasing errors in diurnal scale inversion are complex and difficult to suppress: unmodeled high-frequency non-tidal mass change signals and tidal model errors will alias within a short time window, mapping across scales to the inverted spherical harmonic coefficients and distorting the true diurnal scale time-varying gravity signal; at the same time, the observation noise of inter-satellite distance variability has significant temporal correlation colored characteristics. Traditional inversion methods usually use a diagonal weight matrix based on the root mean square of the residuals to simplify the colored noise into uncorrelated white noise, which cannot accurately characterize the temporal correlation of the noise. This causes the colored observation noise to continue to propagate during the solution of the normal equations, further aggravating the temporal aliasing error.

[0005] Third, existing methods for retrieving gravity fields at daily or short-period scales often employ fixed time windows to directly solve for all-order spherical harmonic coefficients, or estimate residual gravity field signals based on a fixed background field. These methods fail to adequately consider the varying applicability of daily observation data to different orders of spherical harmonic coefficients. Low-order spherical harmonic coefficients correspond to large-scale mass migration and carry the primary daily-scale time-varying gravity information; their unstable estimation significantly reduces the overall accuracy of spatial domain inversion results. Mid- to high-order spherical harmonic coefficients are limited by the geometric constraints of single-day observations; direct full-order solutions easily introduce numerous weak observation parameters, amplifying small singularity noise and further deteriorating the spatial inversion accuracy in polar and high-latitude regions.

[0006] In summary, existing diurnal gravity field inversion techniques cannot simultaneously achieve high temporal resolution, stable low-order signal recovery capability, and controllable mid-to-high-order noise levels. They also cannot systematically suppress temporal aliasing errors caused by colored observation noise, unmodeled high-frequency signals, and ill-conditioned normal equations.

[0007] Therefore, there is an urgent need for a method and system to reduce temporal aliasing errors in gravity satellite inversion. Summary of the Invention

[0008] To address the aforementioned issues, this application proposes a method and system for mitigating temporal aliasing errors in gravity satellite inversion. This method improves the stability and accuracy of diurnal gravity field inversion while maintaining reasonable spatial representation capabilities. The specific details are as follows: On the one hand, this application proposes a method for reducing temporal aliasing errors in gravity satellite inversion, the specific steps of which include: S1. Construct a closed-loop full-scale numerical simulation scenario for the inversion of the daily-scale gravity field, and determine the satellite formation orbit parameters, forward force model, inversion reference force model and instrument noise model within this simulation framework; Error-free observation data containing time-varying gravity signals is generated based on orbital parameters and forward force model. Instrument noise is superimposed on the data to obtain noisy simulated observation data. Meanwhile, reference observation data is generated based on inverted reference force model. By subtracting the simulated observation data from the reference observation data, the orbital observation residuals and the inter-satellite distance variability observation residuals are finally obtained. S2. Estimate the multi-day short-term gravity field based on the orbit observation residuals and inter-satellite distance variability observation residuals, and use it as the short-term background field for the inversion of the daily-scale gravity field; In the estimation of the short-term gravity field, the inter-satellite range-rate residual is first empirically parameterized (K-band range-rate residual empirical parameterization method, KBR method), and then the observation covariance matrix and observation weight matrix are constructed based on the residual autocovariance (Covariance-based weighting method, COV method). In the solution of the normal equation, a Wiese-type signal separation strategy is adopted to reduce the time-domain aliasing error caused by colored observation noise and unmodeled high-frequency signals. S3. Execute the two-step gravity field inversion method (TSM): Use the multi-day short-term gravity field as the background constraint field for the daily-scale gravity field inversion. Based on this, estimate the low-order spherical harmonic coefficients corresponding to the daily-scale residual gravity field based on the single-day observation residuals. In the spherical harmonic domain, combine the multi-day short-term gravity field with the daily-scale low-order residual spherical harmonic coefficients to obtain the daily-scale gravity field result obtained by the two-step method. S4. Implement an optimized two-step gravity field inversion strategy: divide the day-scale spherical harmonic coefficients into stable low-order coefficients and intermediate-order coefficients to be expanded. The stable low-order coefficients are directly estimated using the least squares method. Then, a regularization constraint based on singular value decomposition is established for the extended intermediate-order coefficients. At the same time, the optimal regularization parameter is determined using the generalized cross-validation method. By combining the estimation results of the two types of coefficients, we obtain the optimized day-scale gravity field results after further reducing the time-domain aliasing error.

[0009] Preferably, in S1, the satellite formation adopts a polar-orbiting plus inclined orbit Bender-type binary satellite pair formation mode, including a near-polar orbit satellite pair and a set of inclined orbit satellite pairs. The near-polar orbit satellite pair and the inclined orbit satellite pair are used to provide complementary observation sampling directions to improve the observability of time-varying gravity signals and reduce directional errors.

[0010] In the Bender-type binary satellite pair formation, the near-polar orbit satellite pair has an orbital inclination of 89° and an orbital altitude of 355km, the inclined orbit satellite pair has an orbital inclination of 70° and an orbital altitude of 340km, and the inter-satellite spacing within the satellite pair is 100km.

[0011] The instrument noise model includes an inter-satellite distance variation observation noise model, an accelerometer noise model, a planetarium attitude noise model, a residual non-conservative force noise model, and an orbital position error model; wherein, the orbital position error model is constructed by adding centimeter-level white noise to the position of each satellite.

[0012] The simulation observation data generation adopts a forward and inverse model setting. The forward modeling process uses a static gravity field model, a non-tidal time-varying gravity field model, and a tide model to generate error-free orbit data and error-free inter-satellite distance variability data. Based on the amplitude spectral density of the instrument noise model, instrument noise is added to the error-free orbit data and the error-free inter-satellite distance variability data to obtain noisy simulated observation data; The static gravity field model in the forward modeling process is GOCO05s, the tidal model is EOT11a, and the non-tidal time-varying gravity field model is the AOHIS time-varying gravity model for atmosphere, ocean, hydrology, ice and snow, and solid earth. The static gravity field model in the inversion process is GOCO05s, the tidal model is FES2004, and the AOHIS signal is retained in the observation residuals as the unmodeled diurnal time-varying gravity signal to be recovered.

[0013] The inversion reference force model uses a static gravity field model and a tidal model that differs from the forward modeling process to generate reference orbit data and reference inter-satellite distance variability data, and does not include a non-tidal time-varying gravity field model in the inversion reference force model; The orbit observation residual and the inter-satellite distance variation observation residual are formed by subtracting the noisy simulated observation data from the reference orbit data and the reference inter-satellite distance variation data. The expression for calculating the observation residual is: ; in, Including orbital observation residuals and inter-satellite distance variability observation residuals, The noisy simulated observations generated for the forward modeling process, Reference observations generated for the inversion process.

[0014] Preferably, the multi-day short-term gravity field in S2 is a 7-day short-term gravity field, and the highest spherical harmonic truncation order of the 7-day short-term gravity field is 60th. The 7-day short-term gravity field is used to provide stable short-term spatial background information, so as to avoid ill-conditioned normal equations and noise amplification caused by directly using single-day observation data to estimate all 60th-order spherical harmonic coefficients while maintaining a high spatial representation capability.

[0015] In S2, the process of constructing the observation weight matrix based on the residual autocovariance is as follows: Empirical parameterization modeling is performed on the low-frequency systematic error term in the inter-satellite distance variability residual; The residual autocovariance function is calculated based on the inter-satellite distance variation residual sequence, and a Toeplitz-type observation covariance matrix is ​​constructed based on the residual autocovariance function. A full-array observation weight matrix is ​​constructed from the Toeplitz-type observation covariance matrix to replace the diagonal weight matrix constructed solely from the root mean square of the residuals; The method equations are constructed based on the observation weight array in the full array form to reduce the influence of colored noise in the inter-satellite distance variability residual on the gravity field inversion results.

[0016] It is understandable that the residual autocovariance function is calculated based on the inter-satellite distance variation residual sequence: ; in, Let i be the residual value of the inter-satellite distance variation. The mean of the residuals, The lag order is... The length of the residual sequence. Let the autocovariance be of order k. This represents the (i+k)th inter-satellite distance variation residual value; A Toeplitz-type observation covariance matrix is ​​constructed based on the autocovariance function, and the matrix elements satisfy: ; in, Let be the covariance between the residual values ​​of the i-th and j-th inter-satellite distance variation. For autocovariance; Inverting the observation covariance matrix yields the observation weight matrix in full matrix form: ; in, For the weight matrix, The inverse of the residual covariance matrix; Based on the full-array observation weight array construction equation, the colored noise time-domain correlation characteristics of the inter-satellite distance variation residual are matched.

[0017] Preferably, in S2, the Wiese-type signal separation strategy is a signal separation strategy introduced during the solution of the normal equation. By grouping the spherical harmonic parameters, the long-wavelength high-frequency signals that are prone to aliasing are preferentially mapped to the low-order spherical harmonic coefficients, thereby reducing the propagation of unmodeled high-frequency signals to the mid-to-high-order spherical harmonic coefficients. The multi-day short-term gravity field is a 7-day short-term gravity field, with a highest spherical harmonic cutoff order of 60; the normal equation form corresponding to the weighted least squares is: ; in, The normal equation matrix, For the right-hand side of the normal equation, For observation weight array, To design the matrix, Let be the vector of spherical harmonic coefficients to be estimated.

[0018] Preferably, the two-step gravity field inversion strategy in S3 includes: The multi-day short-term gravity field is added to the reference model of the single-day inversion, so that the reference model contains short-term time-varying gravity information at a multi-day scale; The daily-scale residual spherical harmonic coefficients are estimated based on the daily observation residuals, and the daily-scale residual spherical harmonic coefficients are truncated to a preset low-order range. The multi-day short-term gravity field and the day-scale residual spherical harmonic coefficients are added in the spherical harmonic domain to obtain the day-scale gravity field result that retains the spatial representation ability of the preset highest spherical harmonic order.

[0019] In S3, the calculation expression for the spherical harmonic domain combination of the two-step gravity field inversion strategy is as follows: ; in, For the spherical harmonic coefficients of the multi-day short-term gravitational field, Let be the day-scale residual spherical harmonic coefficients estimated on day d, and The solution is only applied within a preset low-order range; for spherical harmonic coefficients above the low-order cutoff range, the results from the multi-day short-term gravity field are used. The two-step solution for the spherical harmonic coefficients of the gravitational field; The preset low-order range is selected from any one of the 10th, 15th, 20th, 25th, and 30th orders, and is used to determine the low-order cutoff range that can be stably recovered from daily observation data; when the preset low-order range is 15th order, the two-step gravity field inversion strategy achieves a balance between the accuracy of low-order coefficient recovery and noise amplification control.

[0020] Preferably, in S4, the optimized two-step gravity field inversion strategy specifically includes: The day-scale spherical harmonic coefficients are divided into stable low-order coefficients and intermediate-order coefficients to be expanded. The stable low-order coefficients are directly estimated using the least squares method; Subtracting the contribution of the stable low-order coefficients from the right-hand side of the intermediate-order method equation yields the residual method equation corresponding to the intermediate-order coefficients. Singular value decomposition is performed on the residual equations corresponding to the intermediate-order coefficients, and a regularization filter factor is introduced to suppress noise amplification caused by small singular values. The regularization parameter is determined by a generalized cross-validation method with multiple range shrinkages. By combining the directly estimated low-order coefficients with the regularly estimated mid-order coefficients, the optimized diurnal gravity field results are obtained.

[0021] Specifically, construct the day-scale inversion equation in block form: ; in, To stabilize the low-order coefficient vector, The intermediate-order coefficient vector to be expanded. The normal equation matrix corresponding to the lower-order coefficients. The normal equation matrix corresponding to the intermediate-order coefficients. and This is a coupling term between low-order and mid-order coefficients. and These are the right-hand terms of the corresponding normal equation; Directly solve for the stable low-order coefficients: ; in, These are low-order coefficients. The normal equation matrix corresponding to the lower-order coefficients The reverse; After subtracting the lower-order contributions, the intermediate-order residual method equation is obtained: ; Singular value decomposition of the intermediate-order normal equation matrix: ; in, It is a left singular vector matrix. It is a singular value matrix. It is the transpose of the right singular vector matrix; Introducing a regularization filter factor to construct an intermediate-order regularized solution: ; in, The regularization filter factor matrix has the following diagonal elements: ; in, It is a singular value. For regularization parameters, is the regularization filter factor in the direction of the i-th singular value.

[0022] Preferably, the generalized cross-validation method with multiple range shrinkages is a generalized cross-validation method with four range shrinkages, which is used to determine the optimal regularization parameter based on gradually narrowing the search range of the regularization parameter; The optimal regularization parameter is determined using a generalized cross-validation method with multiple range shrinkage steps. The generalized cross-validation function is calculated using the following formula: ; in, To provide the influence matrix corresponding to regularization, This involves matrix trace operations; by progressively narrowing the search interval for the regularization parameter, a value is selected that makes the matrix trace more efficient. The smallest parameter is used as the optimal regularization parameter.

[0023] The stable low-order coefficients are spherical harmonic coefficients up to order 15, and the intermediate-order coefficients to be expanded are spherical harmonic coefficients from order 16 to 25. By performing regularized estimation based on singular value decomposition on the spherical harmonic coefficients from order 16 to 25, more intermediate-order diurnal time-varying gravity signals can be recovered while maintaining the stability of the low-order coefficients.

[0024] The effectiveness of the two-step gravity field inversion strategy and the optimized two-step gravity field inversion strategy is evaluated by the order error, Pearson correlation coefficient, correlation coefficient variance, and root mean square error of the equivalent water height in the spatial domain. If the optimized diurnal gravity field results improve the correlation of the intermediate-order spherical harmonic coefficients and reduce the root mean square error of the equivalent water height in the spatial domain while maintaining the accuracy of the recovery of the low-order spherical harmonic coefficients, then it is determined that the optimized two-step gravity field inversion strategy can further reduce the temporal aliasing error.

[0025] Preferably, the inversion results are evaluated using an index system that combines the spectral and spatial domains, wherein: The first-order error is calculated using the following formula: ; in, The average radius of the Earth Used to characterize the level of spectral error at different orders , To estimate the spherical harmonic coefficients, , The reference spherical harmonic coefficients; The equivalent water height in the spatial domain is calculated using the following formula: ; in, The average radius of the Earth The average density of the Earth, For the density of water, For the first Order load Love number, To fully normalize the associated Legendre function, Latitude Longitude To calculate the maximum spherical harmonic order used, , This represents the difference in spherical harmonic coefficients on day d.

[0026] Preferably, the formula for calculating the root mean square error of the equivalent water height in the spatial domain is: ; in, For the number of days included in the statistics, To estimate the equivalent water height, For reference, the equivalent water level is [high].

[0027] The equivalent root mean square error of the spatial domain water height is a polar or high-latitude weighted equivalent root mean square error of the water height, used to evaluate the spatial striping error and high-latitude spatial noise caused by temporal aliasing error and noise amplification in the daily gravity field results.

[0028] On the other hand, this application proposes a temporal aliasing error reduction system for gravity satellite inversion, comprising: The closed-loop simulation module is used to construct a closed-loop full-scale numerical simulation scenario and generate orbit observation residuals and inter-satellite distance variability observation residuals. The short-term gravity field calculation module is used to estimate the short-term gravity field over multiple days. It has built-in empirical parameterization unit, covariance weighting unit and Wiese-type signal separation unit, and outputs a stable short-term background field. The short-term gravity field calculation module includes: Empirical parameterization unit is used to parameterize the low-frequency systematic error term in the inter-satellite distance variability residual; Covariance weighting unit is used to construct a Toeplitz-type observation covariance matrix and its corresponding full-weight matrix based on the autocovariance function of the inter-satellite distance variation residuals; Wiese-type signal separation unit is used to group low-order and mid-to-high-order spherical harmonic parameters during the solution of normal equations in order to reduce the aliasing propagation of unmodeled high-frequency signals to mid-to-high-order spherical harmonic coefficients. The two-step inversion module is used to estimate the low-order daily-scale residual gravity field and complete the spherical harmonic domain combination against the background of the multi-day short-term gravity field. The two-step inversion module is optimized to distinguish between high and low order spherical harmonic coefficients, enabling low-order direct estimation and regularized extension of mid-order singular value decomposition. The optimized two-step inversion module includes: The coefficient partitioning unit is used to divide the diurnal scale spherical harmonic coefficients into stable low-order coefficients and intermediate-order coefficients to be expanded. A low-order direct estimation unit is used to directly estimate the stable low-order coefficients using the least squares method; The intermediate-order residual method equation construction unit is used to subtract the contribution of the stable low-order coefficients from the right-hand side of the intermediate-order method equation to form the residual method equation corresponding to the intermediate-order coefficients. The SVD regularization unit is used to perform singular value decomposition on the residual equations corresponding to the intermediate-order coefficients and introduce a regularization filter factor to suppress noise amplification caused by small singular values. The GCV parameter selection unit is used to determine the regularization parameters using the generalized cross-validation method. The spherical harmonic coefficient combination unit is used to combine directly estimated low-order coefficients with regularly estimated mid-order coefficients to obtain optimized diurnal gravity field results.

[0029] The results evaluation module is used to evaluate the effect of reducing temporal aliasing error based on spectral and spatial domain indicators.

[0030] In summary, the temporal aliasing error reduction method and system for gravity satellite inversion of the present invention have the following advantages compared with traditional technologies: 1. This invention transforms the daily gravity field inversion from a direct high-order estimation problem into a hierarchical inversion problem of "short-term background field construction - low-order daily-scale correction - intermediate-order regularization extension", which can reduce the impact of sparse daily observation coverage and ill-conditioned normal equations on the recovery of the daily gravity field. 2. This invention constructs a day-scale inversion background field by using a multi-day short-term gravity field, so that the day-scale inversion no longer directly relies on single-day observation data to estimate all high-order spherical harmonic coefficients. This improves the recovery stability of low-order day-scale time-varying gravity signals while maintaining the spatial representation capability of the short-term gravity field. 3. This invention employs empirical parameterization and a full-array weighting method based on residual autocovariance for inter-satellite distance variation residuals, which can more reasonably describe the characteristics of colored observation noise and avoid simplifying time-correlated observation errors into uncorrelated white noise, thereby improving the accuracy of short-term gravity field calculation. 4. In the process of solving the normal equation, this invention introduces a Wiese-type signal separation strategy, which makes long-wavelength high-frequency signals that are prone to aliasing preferentially mapped to low-order spherical harmonics, reducing the risk of unmodeled high-frequency signals propagating to mid-to-high-order spherical harmonics, thereby weakening time-domain aliasing errors. 5. This invention introduces a regularization constraint based on singular value decomposition in the optimization of the two-step inversion and uses a generalized cross-validation method to determine the regularization parameter. This can maintain the stability of low-order coefficient recovery while performing controlled expansion of mid-order spherical harmonic coefficients and suppressing noise amplification caused by small singular values. 6. This invention comprehensively evaluates the daily-scale gravity field results by using the order error, Pearson correlation coefficient, correlation coefficient variance, and root mean square error of the equivalent water height in the spatial domain. It can jointly evaluate the time-domain aliasing error reduction effect from the spectral domain and the spatial domain, and is especially suitable for the generation of future multi-satellite gravity satellite formations and high temporal resolution gravity field products.

[0031] The technical method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0032] Figure 1 This is a diagram illustrating the closed-loop, full-scale numerical simulation framework of the present invention. Figure 2 This is a comparison diagram of the first step gravity field calculation methods of this invention; Figure 3 This is a statistical graph of Pearson correlation coefficients for the two-step gravity field inversion strategy (TSM) of this invention in the range of 2 to 15. Figure 4 This is a statistical graph of Pearson correlation coefficients for the two-step gravity field inversion strategy of this invention in the range of 15 to 30 orders. Figure 5 Pearson correlation coefficient statistics for the optimized two-step gravity field inversion strategy of this invention in the range of 2 to 15 orders; Figure 6 Pearson correlation coefficient statistics for the optimized two-step gravity field inversion strategy of this invention in the range of 15 to 30 orders; Figure 7 This is a flowchart of the two-step strategy of the time-domain aliasing error reduction method and system for gravity satellite inversion according to the present invention. Detailed Implementation

[0033] The technical method of the present invention will be further described below with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps described in these embodiments do not limit the scope of this application.

[0034] The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the scope of this application or its application or use.

[0035] Techniques, systems, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, they should be considered part of the instruction manual.

[0036] In all the examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0037] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning as understood by one of ordinary skill in the art to which this invention pertains.

[0038] This invention addresses the problems of sparse daily orbit coverage, temporal aliasing of unmodeled high-frequency quality signals, colored noise from inter-satellite distance variability observations, and instability of inversion results caused by ill-conditioned normal equations in daily-scale gravity field inversion. It provides a method and system for reducing temporal aliasing errors in gravity satellite inversion. This method addresses the problems of sparse daily orbit coverage, temporal aliasing of unmodeled high-frequency time-varying gravity signals, colored noise from inter-satellite distance variability observations, and instability of inversion results caused by ill-conditioned normal equations in daily-scale time-varying gravity field inversion. It proposes a two-step gravity field inversion strategy consisting of short-term background field construction, daily-scale low-order correction, and intermediate-order regularization extension.

[0039] Example 1 like Figure 1As shown, this embodiment first constructs a closed-loop, full-scale numerical simulation framework, illustrating the relationships between the forward modeling process, the inversion process, the real model, the reference model, error-free observation data, noisy real observation data, reference observation data, and parameter evaluation. This framework includes both forward and inversion processes. The forward modeling process simulates the real observation environment, while the inversion process constructs the reference observation environment required for parameter estimation. By using different background model settings for the forward and inversion processes, background model errors, unmodeled time-varying gravity signals, and observation noise can be explicitly introduced into the closed-loop simulation, thereby evaluating the generation and mitigation effect of time-domain aliasing errors in diurnal gravity field inversion.

[0040] In one optional embodiment, the forward modeling process uses GOCO05s as the static gravity field model, AOHIS as the non-tidal time-varying gravity field model, and EOT11a as the ocean tidal model; the inversion process uses GOCO05s as the static gravity field model and FES2004 as the ocean tidal model, and does not include the AOHIS signal in the inversion reference model. Therefore, the AOHIS signal is retained in the observation residuals as the unmodeled diurnal time-varying gravity signal to be recovered.

[0041] Let the noisy simulation observations generated by the forward modeling process be... The reference observations generated by the inversion process are The observation residuals used for gravitational field inversion are: ; in, This includes orbital observation residuals and inter-satellite distance variability observation residuals. Since the AOHIS signal exists in the forward modeling process but is not included in the inversion reference model, its contribution is retained. In the subsequent gravity field inversion, it is estimated as part of the target time-varying gravity signal.

[0042] In the first step of gravity field calculation, this embodiment estimates the short-term gravity field based on 7-day observation residuals. Let the linearized observation equation be: ; in, To design the matrix, Let be the vector of spherical harmonic coefficients to be estimated. For observation noise, the corresponding weighted least squares equation is: ; in: ; The normal equation matrix, For the right-hand side of the normal equation, For observation weight array.

[0043] For conventional methods, a diagonal weight matrix is ​​directly constructed based on the root mean square of the observation residuals, and then least squares calculation is performed. This method does not specifically address the low-frequency systematic errors and time-dependent colored noise in the inter-satellite distance variability residuals, which can easily lead to noise amplification and striping errors in the mid-to-high-order spherical harmonic coefficients.

[0044] For the KBR+COV scheme, this embodiment first performs empirical parameterization on the low-frequency systematic errors in the inter-satellite distance variation residuals to reduce the impact of low-frequency colored noise on the gravity field solution. Furthermore, the autocovariance function is estimated based on the inter-satellite distance variation residual sequence: ; in, Let i be the residual value of the inter-satellite distance variation. The mean of the residuals, The lag order is... The length of the residual sequence. Let the autocovariance be of order k. Let be the residual value of the (i+k)th inter-satellite distance variation. Construct a Toeplitz-type observation covariance matrix based on the autocovariance function: ; in, Let be the covariance between the residual values ​​of the i-th and j-th inter-satellite distance variation. For autocovariance; The observation weight matrix in full form is then constructed from the covariance matrix: ; in, For the weight matrix, The inverse of the residual covariance matrix; This processing can describe the temporal correlation of inter-satellite distance variability residuals, avoiding the simplification of colored observation noise to uncorrelated white noise.

[0045] For the KBR+COV+Wiese scheme, this embodiment introduces a Wiese-type signal separation strategy during the solution of the normal equations, based on the empirical parameterization of KBR and the full-array weighting of COV. This strategy, through parameter grouping, preferentially maps long-wavelength high-frequency signals prone to aliasing to low-order spherical harmonics, thereby reducing the risk of their propagation to mid- and high-order spherical harmonics.

[0046] like Figure 2As shown, the first-step short-term gravity field solution strategies—the conventional method, the KBR method and the COV method (KBR+COV), and the KBR method, the COV method and the Wiese method (KBR+COV+Wiese)—are illustrated with respect to the AOHIS reference signal, representing the differences in first-step gravity field solution methods. This comparison chart is used to evaluate the error differences among the three strategies in the spectral domain. The evaluation employs the first-step error index. Let the first-step error be... The estimated spherical harmonic coefficients are , The reference spherical harmonic coefficient is , Then the first The order error can be expressed as: ; in, The average radius of the Earth Used to characterize the level of spectral error at different orders. Figure 2 This demonstrates that the KBR+COV+Wiese scheme, compared to conventional methods and the KBR+COV scheme, can further reduce spectral domain errors caused by colored noise due to inter-satellite distance variability and unmodeled high-frequency signals.

[0047] In the second step, this embodiment executes a two-step gravity field inversion strategy. Let the spherical harmonic coefficients of the 7-day short-term gravity field obtained in the first step be... , No. The day-scale residual spherical harmonic coefficients estimated by Tian are: Then the two-step day-scale gravity field solution can be expressed as: ; in, The two-step solution for the spherical harmonic coefficients of the gravitational field. Only the low-order spherical harmonics are estimated within a preset cutoff range; spherical harmonics above this low-order cutoff range retain spatial structure information in the 7-day short-term gravity field. Therefore, the two-step strategy can maintain the spatial representation capability of the short-term gravity field while enabling daily observation data to be primarily used to recover low-order spherical harmonics that exhibit significant diurnal variations and can be stably estimated.

[0048] like Figure 3 and Figure 4 As shown, Figure 3 The Pearson correlation coefficients for the two-step gravity field inversion strategy (TSM) in the range of 2 to 15 are plotted to evaluate the impact of different truncation orders on the stability of low-order diurnal gravity field signal recovery. The truncation orders include 10, 15, 20, 25 and 30, corresponding to TSM-10, TSM-15, TSM-20, TSM-25 and TSM-30 schemes, respectively. Figure 4The Pearson correlation coefficients for the two-step gravity field inversion strategy are plotted in the range of 15 to 30, which are used to evaluate the impact of different truncation orders on the stability of the recovery of the mesoscale gravity field signal. The truncation orders include 10, 15, 20, 25 and 30, corresponding to the TSM-10, TSM-15, TSM-20, TSM-25 and TSM-30 schemes, respectively.

[0049] This embodiment uses the Pearson correlation coefficient to evaluate the recovery stability of the two-step gravity field inversion strategy in the low- and mid-order ranges. The truncation orders include 10, 15, 20, 25, and 30, corresponding to the TSM-10, TSM-15, TSM-20, TSM-25, and TSM-30 schemes, respectively. Figure 3 Corresponding to the order range of 2 to 15, Figure 4 This corresponds to the order range of 15 to 30. The Pearson correlation coefficient can be expressed as: ; in, To estimate the spherical harmonic coefficient sequence, For reference spherical harmonic coefficient sequence, and These are the mean values ​​of the corresponding sequences. Let represent the number of spherical harmonic coefficients included in the statistics. For multi-day results, the variance of the correlation coefficient can be expressed as: ; in, For the first Pearson correlation coefficient for days for Mean correlation coefficient for each day. Figure 3 and Figure 4 This is used to illustrate that the two-step strategy can improve the stability of low-order diurnal signal recovery, but when the second step truncation order is too high, stability is prone to decline in the mid-order range.

[0050] To further extend the recovery stability of low-order coefficients to intermediate-order spherical harmonics, this embodiment constructs an optimized two-step gravity field inversion strategy. The optimization strategy divides the diurnal-scale spherical harmonics into stable low-order coefficients and intermediate-order coefficients to be extended. Let the low-order coefficient vector be... The intermediate-order coefficient vector is The corresponding block-based equation is: ; in, The normal equation matrix corresponding to the lower-order coefficients. The normal equation matrix corresponding to the intermediate-order coefficients. and This is a coupling term between low-order and mid-order coefficients. and This is the right-hand side term of the corresponding normal equation.

[0051] First, conventional least squares estimation is used for the stable low-order coefficients: ; in, These are low-order coefficients. The normal equation matrix corresponding to the lower-order coefficients The reverse; Then, by subtracting the contribution of lower-order solutions from the right-hand side of the intermediate-order equations, the intermediate-order residual method equations are formed: ; Right now: ; Because the intermediate coefficient is poorly observable under single-day observation conditions. This approach is prone to amplifying noise due to small singular values. Therefore, this embodiment performs singular value decomposition on the intermediate-order residual method equations: ; in, It is a left singular vector matrix. It is a singular value matrix. It is the transpose of the right singular vector matrix; And a regularized estimation based on SVD is introduced: ; in, The diagonal elements of the regularization filter factor matrix can be represented as: ; in, It is a singular value. For regularization parameters, Let be the regularization filter factor along the direction of the i-th singular value. To determine the regularization parameter, this embodiment employs a generalized cross-validation method. The GCV function can be expressed as: ; in, This is the influence matrix corresponding to the regularization. In this embodiment, the optimal regularization parameter can be determined using a GCV search method with multiple range shrinkage steps. Finally, the low-order direct estimation result and the mid-order regularization estimation result are combined to obtain the optimized two-step diurnal gravity field result: ; like Figure 5 and Figure 6 As shown, Figure 5To optimize the Pearson correlation coefficient statistics of the two-step gravity field inversion strategy in the range of 2 to 15, and to evaluate whether the SVD regularization extension maintains the recovery stability of the low-order diurnal gravity field signal, the proposed schemes include TSM-10 / 20, TSM-15 / 20, TSM-15 / 25 and TSM-20 / 30. Figure 6 To optimize the Pearson correlation coefficient statistics of the two-step gravity field inversion strategy in the range of orders 15 to 30, and to evaluate the improvement effect of SVD regularization extension on the stability of mid-order diurnal gravity field signal recovery, the proposed schemes include TSM-10 / 20, TSM-15 / 20, TSM-15 / 25, and TSM-20 / 30. This embodiment further uses the Pearson correlation coefficient to evaluate and optimize the two-step gravity field inversion strategy. The schemes include TSM-10 / 20 with m=10 and n=20, TSM-15 / 20 with m=15 and n=20, TSM-15 / 25 with m=15 and n=25, and TSM-20 / 30 with m=20 and n=30. Figure 5 Corresponding to the order range of 2 to 15, it is used to verify whether the SVD regularization extension maintains the stability of low-order day-scale signal recovery; Figure 6 Corresponding to orders 15 to 30, this approach was used to evaluate the recovery stability of mid-order diurnal scale signals. Evaluation results show that optimizing the two-step strategy can improve the recovery stability of mid-order coefficients without weakening the low-order recovery capability.

[0052] To convert the spherical harmonic coefficient results into spatial domain mass variations, this embodiment uses the equivalent water height (EWH) for expression. Let the... The difference between the normalized spherical harmonic coefficients of the gravitational field results and the reference field is: and Then the corresponding spatial location The equivalent water level at a given location can be expressed as: ; in, The average radius of the Earth The average density of the Earth, For the density of water, For the first Order load Love number, To fully normalize the associated Legendre function, Latitude Longitude The maximum spherical harmonic order used in the calculation is given. This formula can convert the difference in spherical harmonic coefficients in the spectral domain into an equivalent change in water height in the spatial domain.

[0053] Let the first The estimated equivalent water level is: The reference equivalent water level is The spatial domain EWH residual is: ; right The root mean square error can be calculated from the results to obtain the spatial domain EWHRMS: ; in, This refers to the number of days included in the statistics.

[0054] From the error distribution perspective, all three methods exhibit some latitudinal differences, but the error amplification is most pronounced in high-latitude regions with the KBR+COV+Wiese method. Its RMS error increases from 2.43 cm in the tropics to 3.23 cm in the mid-latitudes, reaching 4.07 cm in the polar regions. In contrast, TSM-15 and TSM-15-25 significantly suppress errors in mid- and high-latitude regions. TSM-15 reduces errors in mid-latitudes and the polar regions to 2.61 cm and 2.80 cm, respectively, representing reductions of approximately 19.4% and 31.2% compared to KBR+COV+Wiese. TSM-15-25 further enhances the improvement, reducing errors in mid-latitudes and the polar regions to 2.53 cm and 2.51 cm, respectively, corresponding to reductions of approximately 21.8% and 38.3%. Particularly in the polar regions, TSM-15-25 further reduces errors by approximately 10.4% compared to TSM-15, indicating that this method can more effectively mitigate the error amplification problem in high-latitude regions. It should be noted that the TSM method shows a slight increase in error in tropical regions, but overall it achieves a significant improvement in error in mid-to-high latitudes, especially in polar regions.

[0055] like Figure 7 As shown, the relationship between the AOHIS reference signal, closed-loop full-scale numerical simulation, 7-day higher-order gravity field, sky-based lower-order gravity field, and gravity field fusion solution is illustrated. The two-step strategy process of this invention includes: First, inputting the AOHIS reference signal and closed-loop full-scale numerical simulation results into the two-step gravity field inversion strategy; then, in the first step, estimating the 7-day higher-order gravity field based on 7-day observation data to form a stable short-term background field; subsequently, in the second step, estimating the sky-based lower-order gravity field or the regularized extended intermediate-order gravity field based on daily observation data; finally, fusing the 7-day higher-order gravity field with the sky-based lower-order or intermediate-order corrections to obtain the gravity field fusion solution with reduced temporal aliasing error.

[0056] In the system implementation, this invention may include a closed-loop simulation module, a noise simulation module, a short-term gravity field calculation module, a two-step inversion module, an optimized two-step inversion module, and a result evaluation module. The closed-loop simulation module generates forward modeling observation data and inversion reference observation data; the noise simulation module generates instrument noise based on the ASD model; the short-term gravity field calculation module obtains a stable 7-day short-term background field; the two-step inversion module estimates diurnal low-order corrections under the constraint of the 7-day short-term background field; the optimized two-step inversion module estimates intermediate-order coefficients through SVD-GCV regularization; and the result evaluation module outputs the order error, Pearson correlation coefficient, correlation coefficient variance, and EWHRMS results.

[0057] In summary, this invention centers on a two-step gravity field inversion strategy. It weakens the temporal aliasing error in low-order diurnal gravity field inversion by using a 7-day short-term gravity field background constraint, and improves the stability of intermediate-order spherical harmonic coefficient estimation through SVD regularization constraints and GCV parameter selection methods. This method can improve the recovery accuracy of low- and intermediate-order coefficients of the diurnal gravity field while maintaining the spatial representation capability of the short-term gravity field, and reduce the temporal aliasing error caused by colored observation noise, unmodeled high-frequency signals, and ill-conditioned normal equations.

[0058] Finally, it should be noted that the above embodiments are only used to illustrate the technical methods of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical methods of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical methods to deviate from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for reducing temporal aliasing errors in gravity satellite inversion, characterized in that, Includes the following steps: S1. Construct a closed-loop full-scale numerical simulation framework for the inversion of the daily-scale gravity field, and determine the satellite formation orbit parameters, forward force model, inversion reference force model and instrument noise model within the closed-loop full-scale numerical simulation framework. Error-free observation data containing time-varying gravity signals is generated based on orbital parameters and forward force model. Instrument noise is superimposed on the data to obtain noisy simulated observation data. Meanwhile, reference observation data is generated based on inverted reference force model. By subtracting the simulated observation data from the reference observation data, the orbital observation residuals and the inter-satellite distance variability observation residuals are finally obtained. S2. Estimate the multi-day short-term gravity field based on the orbit observation residuals and inter-satellite distance variability observation residuals, and use it as the short-term background field for the inversion of the daily-scale gravity field; In the estimation of the short-term gravity field, the inter-satellite distance variation residuals are first empirically parameterized, and then the observation covariance matrix and observation weight matrix are constructed based on the residual autocovariance. In the solution of the normal equation, a Wiese-type signal separation strategy is adopted to reduce the time-domain aliasing error caused by colored observation noise and unmodeled high-frequency signals. S3. Execute a two-step gravity field inversion strategy: Use the multi-day short-term gravity field as the background constraint field for the daily-scale gravity field inversion. Based on this, estimate the low-order spherical harmonic coefficients corresponding to the daily-scale residual gravity field based on the single-day observation residuals. In the spherical harmonic domain, combine the multi-day short-term gravity field with the daily-scale low-order residual spherical harmonic coefficients to obtain the daily-scale gravity field result inverted by the two-step method. S4. Implement an optimized two-step gravity field inversion strategy: divide the day-scale spherical harmonic coefficients into stable low-order coefficients and intermediate-order coefficients to be expanded. The stable low-order coefficients are directly estimated using the least squares method. Then, a regularization constraint based on singular value decomposition is established for the extended intermediate-order coefficients. At the same time, the optimal regularization parameter is determined using the generalized cross-validation method. By combining the estimation results of the two types of coefficients, we obtain the optimized daily-scale gravity field results after further reducing the temporal aliasing error.

2. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, In S1, the satellite formation adopts a polar-orbiting plus inclined orbit double-satellite pair formation mode, including a pair of near-polar orbit satellites and a pair of inclined orbit satellites, which are used to provide complementary observation sampling directions; The simulation observation data is generated using a forward and inverse model setting. The forward modeling process includes a complete time-varying gravity field model and a tidal model. The inversion reference force model uses different tidal models and does not include a non-tidal time-varying gravity field model; The expression for calculating the observation residual is: ； in, Including orbital observation residuals and inter-satellite distance variability observation residuals, The noisy simulated observations generated for the forward modeling process, Reference observations generated for the inversion process.

3. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, In S2, the process of constructing the observation weight matrix based on the residual autocovariance is as follows: Calculate the residual autocovariance function based on the inter-satellite distance variability residual sequence: ； in, Let i be the residual value of the inter-satellite distance variation. The mean of the residuals, The lag order is... The length of the residual sequence. Let the autocovariance be of order k. This represents the (i+k)th inter-satellite distance variation residual value; A Toeplitz-type observation covariance matrix is ​​constructed based on the autocovariance function, and the matrix elements satisfy: ； in, Let be the covariance between the residual values ​​of the i-th and j-th inter-satellite distance variation. For autocovariance; Inverting the observation covariance matrix yields the observation weight matrix in full matrix form: ； in, For the weight matrix, The inverse of the residual covariance matrix; Based on the observation weight array construction equation in the full array form, the colored noise time-domain correlation characteristics of the inter-satellite distance variation residual are matched.

4. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, In S2, the Wiese-type signal separation strategy groups the spherical harmonic parameters, preferentially mapping easily aliased long-wavelength high-frequency signals to low-order spherical harmonic coefficients, reducing the propagation of unmodeled high-frequency signals to mid-to-high-order spherical harmonic coefficients; the multi-day short-term gravity field is a 7-day short-term gravity field with a maximum spherical harmonic cutoff order of 60; the normal equation form corresponding to the weighted least squares is: ； in, The normal equation matrix, For the right-hand side of the normal equation, For observation weight array, To design the matrix, Let be the vector of spherical harmonic coefficients to be estimated.

5. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, In S3, the calculation expression for the spherical harmonic domain combination of the two-step gravity field inversion strategy is as follows: ； in, For the spherical harmonic coefficients of the multi-day short-term gravitational field, Let be the day-scale residual spherical harmonic coefficients estimated on day d, and The solution is only applied within a preset low-order range; for spherical harmonic coefficients above the low-order cutoff range, the results from the multi-day short-term gravity field are used. The two-step solution for the spherical harmonic coefficients of the gravitational field; The preset low-order range is selected from any one of the following: 10th order, 15th order, 20th order, 25th order, and 30th order.

6. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, In S4, the optimization of the two-step gravity field inversion strategy specifically includes: Constructing a block-based diurnal scale inversion equation: ； in, To stabilize the low-order coefficient vector, The intermediate-order coefficient vector to be expanded. This is the normal equation matrix corresponding to the lower-order coefficients. This is the normal equation matrix corresponding to the intermediate-order coefficients. and This is a coupling term between low-order and mid-order coefficients. and These are the right-hand terms of the corresponding normal equation; Directly solve for the stable low-order coefficients: ;in, These are low-order coefficients. The normal equation matrix corresponding to the lower-order coefficients The reverse; After subtracting the lower-order contributions, the intermediate-order residual method equation is obtained: ; Singular value decomposition of the intermediate-order normal equation matrix: ;in, It is a left singular vector matrix. It is a singular value matrix. It is the transpose of the right singular vector matrix; Introducing a regularization filter factor to construct an intermediate-order regularized solution: ； in, The regularization filter factor matrix has the following diagonal elements: ； in, It is a singular value. For regularization parameters, is the regularization filter factor in the direction of the i-th singular value.

7. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, The optimal regularization parameter is determined using a generalized cross-validation method with multiple range shrinkage steps. The generalized cross-validation function is calculated using the following formula: ； in, To provide the influence matrix corresponding to regularization, This involves matrix trace operations; by progressively narrowing the search interval for the regularization parameter, a value is selected that makes the matrix trace more efficient. The smallest parameter is used as the optimal regularization parameter.

8. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, The inversion results are evaluated using an index system that combines the spectral and spatial domains, where: The first-order error is calculated using the following formula: ； in, The average radius of the Earth Used to characterize the level of spectral error at different orders , To estimate the spherical harmonic coefficients, , The reference spherical harmonic coefficients; The equivalent water height in the spatial domain is calculated using the following formula: ； in, The average radius of the Earth The average density of the Earth, For the density of water, For the first Order load Love number, To fully normalize the associated Legendre function, Latitude Longitude To calculate the maximum spherical harmonic order used, , This represents the difference in spherical harmonic coefficients on day d.

9. The method for reducing temporal aliasing errors in gravity satellite inversion according to claim 1, characterized in that, The formula for calculating the root mean square error of the equivalent water height in the spatial domain is: ； in, For the number of days included in the statistics, To estimate the equivalent water height, For reference, the equivalent water level is [high].

10. A system for reducing temporal aliasing errors in gravity satellite inversion, characterized in that, include: The closed-loop simulation module is used to construct a closed-loop full-scale numerical simulation scenario and generate orbit observation residuals and inter-satellite distance variability observation residuals. The short-term gravity field calculation module is used to estimate the short-term gravity field over multiple days. It has built-in empirical parameterization unit, covariance weighting unit and Wiese-type signal separation unit, and outputs a stable short-term background field. The two-step inversion module is used to estimate the low-order daily-scale residual gravity field and complete the spherical harmonic domain combination against the background of the multi-day short-term gravity field. The two-step inversion module is optimized to distinguish between high and low order spherical harmonic coefficients, enabling low-order direct estimation and regularized extension of mid-order singular value decomposition. The results evaluation module is used to evaluate the effect of reducing temporal aliasing error based on spectral and spatial domain indicators.

Images