A method and system for constructing a dataset of changes in terrestrial water storage, and a storage medium

Abstract

The application provides a land water storage change dataset construction method and system and a storage medium, and belongs to the field of remote sensing and hydrological data processing. The method comprises the following steps: firstly, the missing time series of Mascon products of three sets of CSR, JPL and GSFC datasets are respectively reconstructed and completed by using a singular spectrum iterative interpolation method, so that respective complete space-time data matrices are obtained; secondly, a generalized triangular cap method is used for pixel-by-pixel error analysis of the data matrices, and an error covariance matrix reflecting error variances and covariances is constructed; then, based on the error covariance matrix, the optimal fusion weight of each data matrix is calculated pixel by pixel through a generalized least square criterion; finally, the three sets of data matrices are weighted and fused according to the weight, so that a global land water storage change fusion dataset is generated. The application effectively improves the continuity and precision of the dataset, and provides a reliable data basis for global water cycle research, drought monitoring and water resource management.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing and hydrological data processing, and specifically relates to a method, system and storage medium for constructing a dataset of changes in terrestrial water storage. Background Technology

[0002] The GRACE series of satellite products have been widely used in drought monitoring, glacier mass balance, sea-level change, watershed water resource assessment, and crustal deformation research, providing a unique observational perspective for deeply revealing the mechanisms of mass transfer and energy exchange in the Earth system. Accurate acquisition of data on changes in terrestrial water storage is of great significance for in-depth research on the evolution of regional climate characteristics and for promoting its application in agriculture, ecology, and global change.

[0003] The Gravity Recovery and Climate Experiment (GRACE) and its follow-on mission (GRACE-FO) have achieved high-precision monitoring of Earth's time-varying gravity field by accurately measuring minute changes in the relative distance between binary stars in Earth's gravitational field. Their inverted monthly global gravity field model yields equivalent water height (EWH) data, comprehensively reflecting the quality changes of various water bodies, including surface water, soil water, glacial meltwater, and groundwater, providing an unprecedented technical means for the dynamic monitoring of Terrestrial Water Storage Anomaly (TWSA). There is an approximately 11-month observation gap between the GRACE and GRACE-FO missions, and some months during operation experienced single or double-month data gaps, resulting in multiple data discontinuities in the time series, with a total of 35 data gaps identified. This timing disruption issue is prevalent in various versions of Mascon products released by institutions such as the University of Texas Space Research Center (CSR), the Jet Propulsion Laboratory (JPL), and NASA's Goddard Space Flight Center (GSFC), severely impacting the continuity of long-term time-series analyses.

[0004] To address the issue of missing data, scholars both domestically and internationally have proposed various data reconstruction and interpolation methods, mainly including spectral analysis-based, physical observation-driven, machine learning-based, and multi-source fusion strategies. Spectral analysis methods (such as Singular Spectral Analysis (SSA) and Least Squares Harmonic Estimation (LS-HE)) can effectively separate trend and periodic signals and are robust to missing data. Physical observation-driven methods (such as GNSS vertical deformation-based or Swarm satellite inversion) can recover surface quality changes during periods of missing data using independent observations. Machine learning methods (such as MLP, LSTM, and Transformer) have strong nonlinear fitting capabilities and high reconstruction accuracy, but suffer from overfitting risks and insufficient physical interpretability. In contrast, SSA and its iterative interpolation methods do not require external driving data and can improve temporal continuity and signal stability while maintaining consistency between trends and periods.

[0005] However, interpolation methods can only achieve continuous time series from a single data source. For Mascon products from different institutions such as CSR, JPL, and GSFC, inherent differences in data processing algorithms (e.g., CSR uses Tikhonov regularization, JPL uses continuous Kalman filtering, and GSFC uses iterative solving and autocovariance regularization), background models, and spatial grid designs lead to significant cross-product inconsistencies in the phase, amplitude, period, and spatial details of the retrieved terrestrial water storage change signals. For example, within the same region, CSR time-varying gravity field products typically exhibit a smoother long-term trend, while JPL products, due to less empirical smoothing, may be more sensitive to rapid hydrological extremes, and GSFC products may exhibit different spatial leakage error patterns. Each individual Mascon product has its regional advantages and disadvantages; no single product can demonstrate comprehensive superiority across all regions and time scales globally.

[0006] Researchers can reduce random errors to some extent by performing a simple arithmetic average on multiple Mascon products. However, this method fails to fully consider the complex situation that products such as CSR, JPL, and GSFC have different error magnitudes in different regions and that there may be correlations between errors. Therefore, it is difficult to achieve statistically optimal fusion and cannot quantitatively assess the uncertainty of the fusion result. Summary of the Invention

[0007] To address the aforementioned problems, this invention provides a method, system, and storage medium for constructing a dataset of changes in terrestrial water storage.

[0008] To achieve the above objectives, the present invention provides a method for constructing a dataset of changes in terrestrial water storage, comprising: Based on monthly time resolution, a global land and ocean quality change dataset is obtained over a continuous multi-year period. The global land and ocean quality change dataset includes the CSR dataset, the JPL dataset, and the GSFC dataset.

[0009] The CSR, JPL, and GSFC datasets were segmented, and trajectory matrices were constructed based on the subsequences of each dataset. Singular value decomposition was performed on each trajectory matrix to obtain the temporal principal component matrices of different modes, the singular value matrices of different mode energy magnitudes, and the empirical orthogonal function matrices describing the spatial structure of each mode. A rank-one matrix was constructed using the temporal principal component matrix, the singular value matrix, and the empirical orthogonal function matrix, and the reconstructed components of each mode were obtained based on the rank-one matrix. The reconstructed components were iteratively interpolated and modal-filtered, and then linearly superimposed to obtain the interpolated spatiotemporal data matrix of each dataset.

[0010] Estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices. Construct an error covariance matrix based on the error variance and error covariance matrix. Based on the error covariance matrix, obtain the optimal fusion weight corresponding to each interpolated spatiotemporal data matrix pixel by pixel. Use the optimal fusion weight to weight and fuse each interpolated spatiotemporal data matrix to obtain a fused dataset of global land water storage changes.

[0011] Preferably, the reconstructed components are iteratively interpolated and modal filtered, then linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset, specifically including: Calculate the relative root mean square error of adjacent reconstructed components at the missing position. If the relative root mean square error is greater than a preset convergence threshold, increase the number of modes, reconstruct the trajectory matrix and perform singular value decomposition, and continue iterating. If the relative root mean square error is less than or equal to the preset convergence threshold, the iteration is complete and the mode selection step is entered. For all reconstructed components obtained after iteration, calculate the cumulative power spectrum distribution function of the empirical orthogonal function matrix corresponding to each reconstructed component; if the value of the cumulative power spectrum distribution function at the preset characteristic frequency is greater than the set threshold, then the current reconstructed component is determined to be an effective low-frequency mode and is retained; otherwise, it is determined to be a high-frequency noise mode and is discarded. All the retained reconstructed components are linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset.

[0012] Preferably, the generalized triangular hat method is used to estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices, specifically including: Each time, one of the interpolated spatiotemporal data matrices from the CSR dataset, JPL dataset, and GSFC dataset is selected as a reference benchmark. The difference sequence of the interpolated spatiotemporal data matrices of the other two datasets relative to the current reference benchmark is calculated. Based on the covariance matrix of each difference sequence, the error variance and error covariance of each dataset are solved simultaneously, and three preliminary error covariance estimation matrices are constructed. The median is extracted from the matrix elements representing the same error variance or error covariance in the three preliminary error covariance estimation matrices, and all medians are fused to generate a preliminary fused error covariance matrix; The preliminary fusion error covariance matrix is ​​optimized using the off-diagonal energy minimization method, and the KKT hard constraint of matrix positive definiteness is introduced. The error correlation coefficient in the optimized matrix is ​​thresholded to obtain the error variance and error covariance of different spatial pixels.

[0013] Preferably, based on the error covariance matrix, the optimal fusion weights corresponding to each interpolated spatiotemporal data matrix are obtained pixel-by-pixel using the generalized least squares method. Specifically, this includes: determining whether the error covariance matrix contains valid error covariance; if the error covariances are correlated, the optimal fusion weights are calculated directly based on the inverse of the error covariance matrix using the generalized least squares method; if not, the weights are allocated according to the reciprocal of the error variance of each interpolated spatiotemporal data matrix in the corresponding pixel, with a smaller error variance resulting in a higher weight; if neither of the above two cases applies, equal weights are assigned to the three interpolated spatiotemporal data matrices CSR, JPL, and GSFC, and an equal-weighted average is performed.

[0014] Preferably, the step of segmenting the CSR dataset, JPL dataset, and GSFC dataset respectively includes: firstly, determining the missing length of the three datasets to obtain two scenarios: single / double month missing and consecutive 11 month missing. Based on the singular spectrum iterative interpolation method, a fixed window length and principal component number are empirically set to interpolate and fill in the missing data in each dataset. If there are 11 consecutive months of missing data, a parameter optimization mode is used to handle this. Within the complete time period of each dataset, cross-validation is used to test the interpolation accuracy of different MK combinations for simulated missing data segments. The MK combination with the smallest root mean square error is selected as the fixed window length and principal component number for region adaptation to segment each dataset. Through iterative reconstruction updates and effective mode selection, the retained reconstructed components are linearly superimposed to output a segmented, gap-free continuous sequence. Indicates the window length. This indicates the number of principal components used for reconstruction.

[0015] Preferably, after obtaining the global land water storage change fusion dataset, the method further includes: obtaining a pixel-by-pixel fusion uncertainty field that matches the global land water storage change fusion dataset; the uncertainty field is obtained based on the error covariance matrix and the optimal fusion weights, and a credibility specification is provided for the global land water storage change fusion dataset based on the uncertainty field.

[0016] Preferably, the method further includes: obtaining the simulation results of a preset independent GLDAS-Noah land surface model, and evaluating the accuracy and reliability of the fused dataset by calculating the time series correlation coefficient, root mean square error, and trend consistency index between the fused dataset of global land water storage changes and the simulation results.

[0017] This invention also provides a dataset construction system for changes in terrestrial water storage, comprising: The data acquisition module is used to acquire a global land and ocean quality change dataset over a continuous multi-year period based on a monthly time resolution. The global land and ocean quality change dataset includes the CSR dataset, the JPL dataset, and the GSFC dataset.

[0018] The interpolation module is used to segment the CSR, JPL, and GSFC datasets respectively, and construct trajectory matrices based on the subsequences of each dataset segment. Singular value decomposition is performed on each trajectory matrix to obtain the temporal principal component matrices of different modes, the singular value matrices of different mode energy magnitudes, and the empirical orthogonal function matrices describing the spatial structure of each mode. A rank-one matrix is ​​constructed using the temporal principal component matrix, singular value matrix, and empirical orthogonal function matrix, and the reconstructed components of each mode are obtained based on the rank-one matrix. The reconstructed components are iteratively interpolated and modal-filtered, and then linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset.

[0019] The fusion module is used to estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices. Based on the error variance and error covariance, an error covariance matrix is ​​constructed. Based on the error covariance matrix, the optimal fusion weight corresponding to each interpolated spatiotemporal data matrix is ​​obtained pixel by pixel. The optimal fusion weight is used to weight and fuse each interpolated spatiotemporal data matrix to obtain a fused dataset of global land water storage changes.

[0020] The present invention also provides a computer-readable storage medium storing a computer program that, when loaded by a processor, is capable of executing any of the steps in the method for constructing a dataset of changes in terrestrial water storage.

[0021] The method for constructing a dataset of terrestrial water storage changes provided by this invention has the following beneficial effects: This invention reconstructs missing values ​​using the physical laws inherent in the data itself, independently and intelligently solving the time series discontinuity problem shared by three original datasets without the need for external data, providing a complete and physically consistent foundation for subsequent analysis; Pixel-by-pixel error analysis is performed on the interpolated data, transforming the inconsistencies of each dataset into a quantifiable error covariance matrix, accurately characterizing the spatial reliability differences and correlations of different data sources; Finally, based on this error matrix, appropriate data weights are dynamically assigned to each location globally, achieving a leap from simple averaging to statistical optimization. The generated fused dataset significantly outperforms any single data source in terms of temporal continuity, spatial consistency, and overall accuracy, providing a more reliable data foundation for global water cycle research. Attached Figure Description

[0022] To more clearly illustrate the embodiments and design schemes of the present invention, the accompanying drawings required for this embodiment will be briefly described below. The drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is a flowchart illustrating a method for constructing a dataset of changes in terrestrial water storage according to an embodiment of the present invention. Figure 2 This is a technical roadmap for the singular spectrum iterative interpolation method and data fusion according to embodiments of the present invention; Figure 3 This is a diagram illustrating the missing month information of the GRACE and GRACE-FO time series according to an embodiment of the present invention. Figure 4 This is the optimal reconstruction error probability distribution of the SSA interpolation method in two typical missing scenarios according to an embodiment of the present invention; Figure 4 (a) indicates missing data for single and double months. Figure 4 (b) represents 11 consecutive months of missing data; Figure 5 This is a diagram showing the energy distribution and cumulative contribution rate of the SSA principal components in an embodiment of the present invention. Figure 6 This is a time-series morphological distribution diagram of the first five reconstructed components in an embodiment of the present invention; Figure 7 The CDF test for high-frequency noise removal in this embodiment of the invention; Figure 8 This is a diagram illustrating the GTCH iteration process in an embodiment of the present invention. Figure 8 (a) is the eigenvalue of R. Figure 8 (b) represents the variance of each dataset. Figure 8(c) is a graph showing the relationship between the number of iterations and the objective function. Figure 8 (d) represents the number of iterations and the KKT constraints. Figure 8 (e) is a graph showing the relationship between the number of iterations and energy. Figure 8 f represents the iteration number and the parameter iteration trajectory; Figure 9 This is an uncertainty distribution diagram of an embodiment of the present invention. Detailed Implementation

[0024] To enable those skilled in the art to better understand and implement the technical solutions of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and should not be construed as limiting the scope of protection of the present invention.

[0025] The main objective of this invention is to simultaneously address the temporal discontinuities and cross-institutional product differences in GRACE / GRACE-FO Mascon data from 2002 to 2024. To achieve this, firstly, the Singular Spectrum Iterative Interpolation (SGIT) method is used to reconstruct the GRACE / GRACE-FO time series. Missing monthly signals are recovered through iterative decomposition and principal component screening, maintaining consistency in trend and seasonal components and ensuring temporal continuity. Secondly, the Generalized Triangular Hat (GTCH) method is used to estimate the pixel-level uncertainty and error correlation of different Mascon products, quantitatively characterizing their error differences. Finally, a hierarchical weighted fusion strategy is implemented based on the covariance matrix. When GLS fusion fails, it degenerates sequentially into inverse variance weighting and equal-weighted averaging to achieve robust optimal fusion of multi-source products. The optimal fusion weight for each product at the pixel scale is determined, integrating the advantages of multi-source Mascon products to generate a high-precision dataset of terrestrial water storage changes. The fusion results are further independently validated with the GLDAS model to evaluate its accuracy improvement effect at global and regional scales.

[0026] Figure 2This paper presents the overall framework of a method for constructing a fusion dataset of global land water storage changes based on singular value decomposition (SVD) and GTCH covariance optimization. The method consists of two core modules: the left half is the SSA iterative interpolation module, which processes three Mascon products (CSR, JPL, and GSFC) respectively (time range: April 2002 to December 2024). First, the missing length is determined, distinguishing between single / double month missing data and consecutive 11-month missing data. Then, a trajectory matrix is ​​constructed through time series embedding, and the principal component matrix, singular value matrix, and empirical orthogonal function matrix are obtained through SVD. A rank-one matrix is ​​then constructed, and reconstructed components of different modes are extracted. A fixed-parameter mode is used to handle single / double month missing data, while a parameter optimization mode is used to handle consecutive 11-month missing data. Through iterative reconstruction updates and effective mode selection (based on the power spectrum CDF criterion), the retained reconstructed components are linearly superimposed to output a gap-free continuous sequence. The right half is the multi-source fusion module, which takes three sets of interpolated Mascon products as input. It constructs a difference sequence through pairwise sequences to remove the true signal. It adopts an alternating reference strategy, combined with objective function construction, KKT hard constraints and off-diagonal energy minimization optimization, and obtains the error covariance matrix through reference median synthesis. On this basis, a hierarchical fusion strategy is implemented: it first adopts generalized least squares (GLS) fusion. If the covariance matrix is ​​not positive definite, it degenerates to inverse variance weighting. If it is still unstable, it further degenerates to equal weighted averaging. The robustness of fusion is ensured through the stepwise strategy. Finally, it outputs a time-series continuous fused product and generates a pixel-by-pixel fusion uncertainty field.

[0027] Based on this, the present invention provides a method for constructing a dataset of terrestrial water storage changes, specifically a method for fusing a global terrestrial water storage change dataset based on singular spectrum reconstruction and GTCH covariance optimization, as detailed below. Figure 1 As shown, it includes: S1. Based on monthly time resolution, obtain global land and ocean quality change datasets over consecutive years. The global land and ocean quality change datasets include the CSR dataset, JPL dataset, and GSFC dataset.

[0028] Mascon products greatly simplify the data processing workflow, allowing direct use of gridded equivalent water height data without the need for complex preprocessing of spherical harmonic coefficients. Furthermore, by employing prior geophysical models (such as JPL's continuous Kalman filtering) or introducing regularization constraints (such as CSR's Tikhonov regularization method), the products effectively suppress land-sea leakage and improve the accuracy of mass change estimation. This invention utilizes the GRACE / GRACE-FO Mascon global land and oceanic water storage anomaly (TWSA) products released by three official institutions, covering the period from April 2002 to December 2024. The three Mascon products are from the University of Texas Center for Space Research (CSR), the Jet Propulsion Laboratory (JPL), and NASA's Goddard Space Flight Center (GSFC), respectively, all providing global coverage, monthly temporal resolution, and using centimeters of equivalent water height (cmEWH) as the unit. The fusion scheme of this invention is based on the CSR dataset, JPL dataset and GSFC dataset. The main parameters and processing characteristics of these three products are shown in Table 1.

[0029] Table 1. Comparison of the latest Mascon product parameters from three official organizations. like Figure 3 As shown, due to the handover of satellite missions and observational anomalies, some months are missing, especially the 11-month observation gap from July 2017 to May 2018, which seriously affected the ability to continuously characterize seasonal changes and interannual fluctuations. Figure 3A heatmap using "year × month" format was used to display the missing data in the GRACE / GRACE-FO time series from April 2002 to December 2024: blue represents normal months (GRACE satellite began normal operation in April 2002), blank represents single-month missing data, some years have consecutive double-month missing data, and there is also a 11-month missing data period from June 2017 to May 2018 during the mission handover gap. Throughout the period from April 2002 to December 2024, single-month and double-month missing data scenarios accounted for approximately 4%; the 11-month consecutive missing data scenario accounted for 4%, but this is a missing data pattern more sensitive to sequence continuity, and the total missing data accounted for 12.8% of the complete time series. The study only used the latest publicly available Release data: RL06 series GRACE / GRACE-FO products from three processing centers (GFZ, JPL, and CSR) were selected to ensure version consistency and support for missing data handling and fusion analysis.

[0030] To obtain the complete time series, this invention, in the subsequent method section, independently performs singular spectrum iterative interpolation (SSII) reconstruction on three Mascon products: JPL, CSR, and GSFC. The interpolation only completes the time dimension, while each product retains its original spatial resolution (0.25° for CSR, and 0.5° for JPL and GSFC). The three sets of Mascon data after complete interpolation will be used as input for subsequent uncertainty estimation and multi-source fusion.

[0031] The TWSA results of GRACE / GRACE-FO inversion of terrestrial water storage changes were independently validated using the Noah land surface model dataset GLDAS-Noah v2.1 provided by the Global Land Data Assimilation System (GLDAS), jointly developed by NASA's Goddard Space Flight Center (GSFC) and NOAA's National Environmental Prediction Center (NCEP). The GLDAS-Noah v2.1 product has a grid size of 0.25° × 0.25°, a monthly temporal resolution, and covers the period from April 2002 to December 2024. The variable units are [missing information]. (Equivalent to millimeter water layer thickness). Data source: NASA Goddard Earth Sciences Data and Information Services Center, GESDISC. The terrestrial water storage variables output by the GLDAS-Noah model include the moisture content of four soil layers, snowmelt equivalent, and canopy water. To maintain consistency with the terrestrial water storage changes estimated by GRACE, the moisture content of the four soil layers, snowmelt equivalent, and canopy water are summed to calculate the total shallow surface water storage change simulated by GLDAS, as shown in the following expression:

[0032] Where SWSt represents the total shallow surface water storage. For time index, representing the first Month; SoilMoi t,0-200cm represents the soil moisture content from the surface to a depth of 200cm in the current month, which is obtained by summing the soil moisture content of four layers (0–10cm, 10–40cm, 40–100cm, 100–200cm); SWE t represents the snow cover water equivalent in the current month; CanopInt t represents the canopy water in the current month.

[0033] S101, Physical meaning of each component To eliminate systematic bias and maintain comparability with GRACE data, the GLDAS-Noah data underwent mean-reduction processing based on the 2004–2009 average to generate a water storage anomaly SWS sequence. Spatially, an area-weighted averaging method was used to resample the original grid size from 0.25°×0.25° to a uniform 0.5°×0.5° global grid, maintaining consistency with the spatial resolution of the GRACEMascon product. The processed ΔSWS was compared with the three Mascon products and their fusion results at global and regional scales to verify the temporal consistency and reliability of the changing trends across different data sources.

[0034] S2. Segment the CSR, JPL, and GSFC datasets respectively, and construct trajectory matrices based on the subsequences of each dataset segment. Perform singular value decomposition on each trajectory matrix to obtain the temporal principal component matrices of different modes, the singular value matrices of different mode energy magnitudes, and the empirical orthogonal function matrices describing the spatial structure of each mode. Construct a rank-one matrix using the temporal principal component matrix, singular value matrix, and empirical orthogonal function matrix, and obtain the reconstruction components of each mode based on the rank-one matrix. After iterative interpolation and mode filtering, linearly superimpose the reconstruction components to obtain the interpolated spatiotemporal data matrix of each dataset.

[0035] This invention employs Singular Spectrum Analysis (SSA) iterative interpolation to fill in missing months in GRACE / GRACE-FO monthly time series. This method effectively denoises time series data, fills in missing values, and achieves stable reconstruction of low-frequency signals. The SSA method has been applied in numerous studies and its implementation has been publicly disclosed. This invention strictly adheres to its mathematical framework and implementation logic, extending its application and fine-tuning the method for the dataset used in this invention, including the inclusion of parameter search and parallel computing strategies.

[0036] To address the issue of missing months in GRACE / GRACE-FO terrestrial water storage change data, embedding processing is first performed on its time series. Assume the GRACE / GRACE-FO terrestrial water storage change time series for one pixel is as follows: Missing observations are represented by NaN. Singular spectral analysis requires a given window length and an embedding dimension of [value missing]. time series Embedded into a trajectory matrix Hankel, the columns of the matrix contain The time interval of consecutive elements. Construct the trajectory matrix and the Hankel matrix. Let X be a column of M consecutive line segments, expressed as follows:

[0037] Where N is the length of the time series (number of months); R represents a real-field matrix. Let be the row number of the trajectory matrix. is the column number of the trajectory matrix; the th The corresponding length of the row is continuous subsequences During the iteration process, missing time points are replaced by the current reconstruction result, ensuring that the matrix remains consistent in each iteration. It can be computed, thus enabling subsequent singular value decomposition and reconstruction.

[0038] Based on the method proposed by Broomhead and King in 1986, the trajectory matrix is... Perform Singular Value Decomposition (SVD), a method that transforms the matrix... It is decomposed into the product of three matrices: the principal component matrix of time, the singular value matrix, and the empirical orthogonal function matrix. in, This decomposition shows that the trajectory matrix can be represented as a linear superposition of several rank-1 structures, each structure corresponding to a mutually orthogonal mode. (Left singular vector matrix) The column vectors describe the modes in the sliding window dimension (length). The time coefficient variation on the right singular vector; The column vectors characterize the modalities in the embedding space (window length). The structural form in ) is usually called the empirical orthogonal function (EOF); Let be a singular value matrix, denoted as The singular values ​​are sorted according to σ1≥σ2≥…≥σ_M, and each singular value… Measure the first The importance of each modality is determined by its numerical value. A higher value indicates a greater contribution of that modality to the original data. The corresponding eigenvalue is defined as follows: Its normalized form It can be used to characterize the energy contribution of different modes.

[0039] Before selection Reconstruct the mode, and reconstruct the mode. Each mode, and its corresponding rank-one matrix is ​​defined as follows: : in, , for the first A number of empirical orthogonal functions (EOF) describe the modality in the embedding space (window length). The structural form in ); , for the first Each principal component (PC) represents the mode in the time direction (length). The coefficient of variation on the modality reflects the intensity of the change of the modality over time. With trajectory matrix The same type is the complete representation of the i-th mode in the trajectory matrix space.

[0040] To elucidate the modal physical characteristics and energy distribution patterns of SSA decomposition, Figure 5 The time series decomposition results of typical pixels are shown, and the normalized eigenvalues ​​(blue line) and their cumulative contribution rates (orange line) of each principal component are given. The energy is highly concentrated in the first 2-3 low-order modes. The contribution rate of the first principal component exceeds 35%, and the cumulative contribution rate of the first 3 modes is over 90%. Subsequently, the eigenvalues ​​decay rapidly and tend to level off, indicating that the main information of the original sequence can be characterized by a few dominant modes. Figure 6The temporal morphology of the first five reconstructed components (RC1–RC5) based on the optimal window length M and the number of modes K is presented. RC1 and RC2 exhibit regular periodic oscillations with large amplitudes and orthogonal phases, corresponding to annual and semi-annual seasonal signals, respectively. RC3 represents a slowly changing low-frequency trend, reflecting long-term background changes in water storage. RC4 and RC5 have smaller amplitudes and higher frequencies, mainly reflecting high-frequency disturbances or local anomalies. These results validate the effective separation capability of the SSA method for GRACE / GRACE-FO time-varying gravity field signals, providing a physical basis for subsequent mode selection and iterative reconstruction of missing values ​​based on the power spectrum CDF criterion.

[0041] Since the trajectory matrix is ​​constructed using a sliding window embedding, its elements satisfy... in, For the row index of the trajectory matrix, ; For column indices of the trajectory matrix, ; From the above equation, it can be seen that all matrices satisfying the following condition... All elements correspond to the same time point in the original time series. , Indicates the original time series at the 1st... The value at each time index. Therefore, each anti-diagonal line in the trajectory matrix corresponds to the same time index in the time series.

[0042] Get the first rank-one matrix Next, it needs to be restored from the two-dimensional trajectory matrix space to a one-dimensional time series. This process is achieved by applying diagonal averaging, also known as Hankelization. That is, averaging the matrix along the anti-diagonal line restores it to a one-dimensional time series, yielding the corresponding reconstructed components. The Hankelization operator is then defined. Then the first Each reconstructed component is represented as:

[0043] Its specific form is as follows: in, For the first The number of elements on each anti-diagonal line. The final result is... For length is The time series reconstruction component. This process essentially involves averaging multiple estimates of the same time point across different sliding windows to obtain the final reconstruction result for that time point; The Hankelization operator represents the anti-diagonal mean, which transforms the matrix... Restore to length A one-dimensional sequence; For the first Length of each reconstructed component .

[0044] To further fill in missing values, an iterative strategy of gradually increasing the number of modalities is adopted. In the initial stage, the missing values ​​are initialized, and then... SSA decomposition and reconstruction begin. In each iteration, based on the previous... The reconstructed form of the time series obtained from the reconstructed components is as follows:

[0045] in, Before use The time series obtained from modal reconstruction is at time... The estimated value (that is, the current iteration relative to the true value) (reconstruction) The number of modes participating in the reconstruction (using from...) Gradually increase to the target (strategy) Indicates the first Each reconstructed component at time... The value, it comes from the first rank-one matrix The Hankelization operator is then used. The values ​​at the missing positions are then updated using the current reconstruction result. Let... Given a set of missing time indices, the iterative convergence criterion is defined as:

[0046] Among them, molecules It is the RMS of the missing position "this estimate vs. the previous estimate"; the denominator is the RMS of the missing position "this estimate vs. the previous estimate". It is the RMS of the last estimate at the missing position, which is used for normalization to obtain the relative change. This is the set of missing indexes (which months are NaN). For example, if the missing months are from July 2017 to May 2018, then the corresponding indexes for these months are... All belong to ; For the number of iterations (the th iteration) (nth iteration) Indicates the first After the next iteration, at time point Reconstruction estimation; This represents the estimate from the previous iteration; For the preset threshold ( When the above conditions are met, the iteration is considered to have converged, and the final filling result is obtained. If the conditions are not met, the number of modes continues to increase. This continues until the preset maximum number of modes is reached.

[0047] After iterative convergence, not all modes participate in the final signal reconstruction. To suppress high-frequency oscillations (HFO), this embodiment performs spectral analysis on the empirical orthogonal function (EOF) for each mode. Figure 7 As shown, the cumulative distribution function (CDF) of the power spectrum is used as the screening criterion: the power spectral density of each EOF is calculated and a CDF curve is constructed. A cutoff frequency of 3 times / year (with a period ≥ 4 months) is used as the threshold. If the CDF value of a mode at the cutoff frequency is ≥ 0.9, it is considered a valid low-frequency mode and retained; otherwise, it is considered high-frequency noise and discarded. The red curve in the figure represents the energy accumulation process of each mode. The left side of the curve that rapidly rises above 0.9 represents the real hydrological signal with energy concentrated in the low frequency, while the right side of the curve close to the gray diagonal dashed line exhibits white noise characteristics. This criterion effectively separates the trend term, seasonal cycle, and high-frequency disturbances, incorporating only physically consistent low-frequency modes into the final reconstruction, thus avoiding the introduction of non-physical fluctuations by higher-order oscillation components within long-term gaps.

[0048] When the number of iteration modes reaches the preset maximum number of modes After convergence, frequency domain filtering is performed on each mode. For the ... For each EOF, calculate its power spectral density (PSD) and construct the cumulative distribution function (CDF):

[0049] in, It calculates the cumulative energy proportion of the power spectral density of the i-th EOF below frequency f; the numerator represents the total energy at frequencies less than or equal to f; the denominator represents the total energy across all frequencies; the ratio represents the proportion of low-frequency energy. For the first Each mode at discrete frequencies Power spectral density at; The total number of frequency discrete points; the power spectral density is estimated by the periodogram method.

[0050] According to the criterion proposed by Yi and Sneeuw, if there are 3 cycles per year (i.e., a cycle of ≥4 months), then... If the following conditions are met: If the modal energy is concentrated in the low-frequency range, it is considered an effective mode; otherwise, it is considered high-frequency noise.

[0051] ; in, As a binary discrimination flag (0 / 1), htest(i)=1 indicates that the i-th mode is a valid mode that passes the low-frequency criterion and will be retained; htest(i)=0 indicates that the mode does not pass the high-frequency noise criterion and needs to be discarded.

[0052] Window length The number of principal components reconstructed The method was determined through cross-validation of the complete time series from 2004 to 2015 to minimize the error in known observation years. This process was performed independently at each pixel, ensuring the method's adaptability on a global scale. Specifically, within the complete observation period of 2004-2015, a leave-one-out cross-validation strategy was employed. Each time, observation data for a particular year was manually removed, and the remaining data was used to perform SSA interpolation reconstruction based on different MK combinations. The root mean square error (RMSE) between the reconstructed values ​​and the actual observation values ​​was calculated. After traversing all possible MK combinations, the MK combination that minimizes the RMSE was selected as the optimal parameter combination for that pixel. Because this process was performed independently at each pixel, regional adaptive optimization of the parameters was achieved.

[0053] To improve computational efficiency, a 40×40 spatial block parallel processing strategy is adopted. The `parfor` parallel loop in MATLAB is used to enhance computational efficiency, automatically skipping processed blocks and saving intermediate results, and supporting breakpoint continuation. Each pixel outputs a complete interpolation sequence, error estimates verror1 and verror2 for stages A and B, and its optimal SSA parameters. All blocks are eventually merged to generate a structured 3D TWSA dataset. This method effectively improves the temporal continuity and signal stability of Mascon data, providing fundamental support for subsequent multi-source product fusion and uncertainty modeling.

[0054] To verify the reconstruction accuracy of the SSA interpolation method under different missing scenarios, this embodiment uses cross-validation to perform optimal reconstruction error statistics for two typical scenarios: single / double month missing and consecutive 11 month missing. Figure 4 The optimal reconstruction error probability distribution of the SSA interpolation method is shown in two typical missing scenarios. Figure 4 (a) represents the scenario of sporadic missing data in odd / even months. Figure 4 (b) represents the 11-month missing data scenario during the handover period between the GRACE and GRACE-FO tasks from July 2017 to May 2018. In both scenarios, the error distribution exhibits a sharp peak near 0 cm, indicating that the method has excellent reconstruction accuracy for the vast majority of grids. In contrast, the distribution tail of the 11-month missing data scenario is slightly wider, with a small number of samples showing larger errors of 80-100 cm, suggesting that long-term continuous data gaps can slightly reduce reconstruction stability, but overall, it still maintains a good level of error control.

[0055] S3. Estimate the error variance of each interpolated spatiotemporal data matrix in different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices. Construct an error covariance matrix based on the error variance and error covariance. Based on the error covariance matrix, obtain the optimal fusion weight corresponding to each interpolated spatiotemporal data matrix pixel by pixel. Use the optimal fusion weight to perform weighted fusion of each interpolated spatiotemporal data matrix to obtain a global land water storage change fused dataset.

[0056] After completing singular spectral decomposition and mode filtering based on spectral CDF, the final time series used for physical signal reconstruction consists of a linear superposition of all valid modes that have passed the spectral filtering: ; in, The reconstructed time series of the final physical low-frequency water storage signal consists only of the linear superposition of low-frequency effective modes selected by spectrum filtering. This is a monthly-scale time index sequence. Each month; The length of the time series indicates the specific number of months. The reconstructed time series, used to represent the physical low-frequency water storage signal, is represented by the one obtained after removing high-frequency oscillation noise modes. For the first The reconstructed component of the SSA mode is derived from the first... A one-dimensional time series is obtained by Hankel anti-diagonal averaging of rank-1 matrices; Modal numbering, usually ;in This represents the upper limit of the number of modes output in the algorithm. This is the set of valid modal indexes filtered by spectrum.

[0057] At the locations where the original time series observations are missing, this reconstructed value can serve as the final fill-in result for the SSA-B stage. Because this fill-in value consists entirely of low-frequency, physically consistent modes, it outperforms results based solely on full mode reconstruction in terms of temporal continuity, amplitude stability, and spectral structure, thus avoiding the introduction of non-physical fluctuations by higher-order oscillations within long-term gaps. The final reconstructed sequence not only includes long-term trends, seasonal variations, and interannual low-frequency signals, but also effectively eliminates high-frequency oscillation noise during the process.

[0058] To rationally allocate the fusion weights of the three Mascon products (CSR, JPL, and GSFC), it is necessary to first estimate their pixel-by-pixel uncertainty, i.e., error variance and error covariance. This invention adopts the GTCH method improved by Duanetal in 2024. The observations of the three products are represented as follows: ,in , It can be decomposed into real signals Product error ,satisfy:

[0059] ; The goal is to estimate the error covariance matrix, constructing the error covariance matrix among the three products at each cell, in the form: in, For the error variance, This is the error correlation coefficient.

[0060] This invention employs a rotating reference strategy for three Mascon products (CSR, JPL, and GSFC). In the In the round, define the set of non-reference products. Construct the difference vector: Construct difference sequences for any two products Thus eliminating the real signal By using the differences between the three products to eliminate the true signal, the variance of the difference sequence satisfies:

[0061] And calculate its difference covariance matrix: Where cov() is the covariance operator, If the matrix is ​​a symmetric positive definite matrix (SPD), then the cell is discarded. Then, the parameters to be estimated are introduced. :

[0062] in, The vector of parameters to be optimized; Represents a two-dimensional real vector. Error covariance vector between the reference product and other products (corresponding) , ); The error variance term for the reference product (the element in the bottom right corner of the matrix); It is an all-1 vector. The error covariance matrix is ​​constructed according to the structural relationship of GTCH. .

[0063] in, Free parameters Constructed candidate matrix; for The top left block (the error covariance submatrix between non-reference products); to obtain a stable and feasible The optimization estimation with hard constraints is adopted: the criterion is to minimize the "off-diagonal energy" of the error covariance matrix. ; in, For matrix The Line number Column elements (covariance / variance); It is the sum of squares of all off-diagonal elements in the upper triangle (off-diagonal energy). Normalization coefficient ( ).

[0064] Furthermore, KKT hard constraints (Karush-Kuhn-Tucker) are introduced to guarantee the feasibility and positive definite structure of the solution. Definition ( (This represents the vector used in the constraint). The constraint is written as: ( The KKT hard constraints are applied (to ensure numerical stability). The optimization problem is solved using the interior-point method (fmincon is used in the implementation) to obtain the optimal parameters θ* (θ is the parameter being tried, and θ* is the optimal parameter obtained after final convergence) and to form the error covariance estimate of the reference wheel. The final per-pixel error variance estimate for the three products is obtained by taking their diagonal elements:

[0065] in, For the first Variance estimation for each product.

[0066] Through objective function optimization and KKT hard constraints, the final result is... It reflects the error structure between data sources. The matrix contains the variance of the errors (i.e., the diagonal elements) and the correlation between the errors (i.e., the off-diagonal elements).

[0067] To verify the optimization stability and convergence characteristics of the GTCH algorithm, Figure 8 The iterative optimization process and convergence characteristics of the GTCH generalized triangular hat are demonstrated. Figure 8 (a) shows the changes in the eigenvalues ​​of the error covariance matrix R from index 1 to 3. All three are positive and remain positive, tending to be stable, indicating that the matrix always maintains positive definiteness and the optimization process is numerically stable. Figure 8Figure (b) shows the error variance estimates for the three products, CSR, JPL, and GSFC. JPL has the largest variance, approximately 16, followed by CSR at approximately 11, and GSFC has the smallest variance, approximately 6, quantitatively characterizing the differences in reliability among products from different organizations. Figure 8(c) shows the objective function. The convergence curve of the ) shows that its value decreases stepwise with the number of iterations and eventually approaches zero, reflecting the achievement of the off-diagonal energy minimization objective. Figure 8 (d) represents the trajectory of the violation under the KKT constraint. It approaches zero in the later stages of the iteration, approximately -0.009, indicating that the positive definite constraint condition is satisfied. Figure 8 (e) represents the decay process of the energy of the off-diagonal elements of matrix R, and its continuous decrease verifies the effectiveness of decorrelation. Figure 8 f is the optimization parameter. The iterative trajectory, h1, h2 and reference variance All values ​​tend to stabilize, indicating that the algorithm has reached a stable convergence state. Combining the six subgraphs, the GTCH algorithm achieves optimized estimation of the error covariance matrix through approximately 40 iterations, providing a reliable quantitative description of the error structure for subsequent GLS fusion.

[0068] To reduce the sensitivity of reference product selection to GTCH estimation, this invention applies a different approach to each pixel. Three sets of error covariance estimation matrices were obtained for reference. Subsequently, the median of each of the six independent elements (three variances and three covariances) of the covariance matrix was taken element-wise to help balance the errors between the various products and construct a robust composite error covariance matrix. :

[0069] in, For the estimation of the error variance of the three products, This is for estimating the error covariance.

[0070] Obtain the pixel-level error covariance matrix Subsequently, the three Mascon products were fused using a GLS-based generalized least squares error optimal weight determination method. Since the error correlation of GRACE products is generally weak, the correlation coefficient in the error covariance matrix may be very small in some cases. To suppress numerical noise, this invention performs threshold pruning on extremely small correlation coefficients, setting them to 0. Calculate the correlation coefficient :

[0071] like Then let And set up synchronously (That is, the covariance term is clipped to 0), so that the error structure of unstable pixels degenerates into a diagonal matrix form when necessary, thereby improving the stability of the fusion weight solution.

[0072] Obtain the pixel-level error covariance matrix Then, the error-optimal weights (weights) based on Generalized Least Squares (GLS) are adopted. The method for fusing three Mascon products is determined. Compared to traditional inverse variance weighting based solely on variance, the GLS of this invention directly uses the complete covariance matrix, allowing the covariance of error correlation to participate in the solution of fusion weights. The goal of fusion is to utilize the confidence information in Σ to give higher weights to products with smaller errors, while considering the error correlation among the three products. The goal of fusion is to minimize the error variance of the fusion result, specifically as follows:

[0073] in, ; w represents the fusion weight of the three products. The MATLAB quadprog algorithm is used to solve the problem to ensure the non-negativity and numerical stability of the weights.

[0074] like When the matrix is ​​a general symmetric positive definite matrix (i.e., the error is correlated): The weight calculation formula is: when When the off-diagonal matrix is ​​0 (i.e., the correlation is clipped to 0), The above formula degenerates into inverse variance weighted: When the pixel-level error covariance matrix does not meet the symmetric positive definite condition or when numerical anomalies occur during the weight calculation process, the arithmetic mean of the three products is used as the fusion result for that pixel to ensure the spatial continuity and numerical stability of the fusion result.

[0075] in, These represent CSR, JPL, and GSFC respectively in the month. The observations ensured the continuity and stability of the fused field.

[0076] For any month's CSR(A), JPL(B), and GSFC(C) products, the fusion result is: This weighted fusion method can produce a unified dataset of global land water storage anomalies, thereby reducing the uncertainty of individual products and improving the overall data quality. Figure 9 This demonstrates the improvement in uncertainty of the fused dataset compared to the single product. The horizontal axis represents the uncertainty ratio (r = 1 / 3 of the fused data versus the single product). The vertical axis represents the empirical cumulative distribution function F(r). When r < 1, the fusion uncertainty is lower than that of the single product; when r = 1, the two are comparable; when r > 1, the fusion is worse. The four curves compare the uncertainty differences between the fusion result and CSR (blue), JPL (red), GSFC (yellow), and the single product best (purple). Taking the CSR curve as an example, This indicates that the fusion uncertainty of 50% of the pixels has decreased to below 47% of that of the CSR; This indicates that 100% of the pixels meet the fusion uncertainty requirement, which is no higher than that of CSR. The JPL and GSFC curves have similar shapes, with medians of 0.41 and 0.44 respectively, showing a slightly better improvement than CSR. The purple control curve, using the minimum uncertainty of the three products for each pixel as a benchmark, shifts to the right, indicating that while the fusion result is better than any single product, the improvement potential relative to the theoretically optimal single product is narrowing, with a median of 0.71. This figure quantitatively verifies the uncertainty reduction benefit of GLS fusion globally, providing a basis for the reliability assessment of the fused dataset.

[0077] After obtaining the optimal fusion weights for each pixel, this invention further generates a pixel-by-pixel fusion uncertainty field for the fused dataset. The calculation of this uncertainty field is based on two key elements: first, it is estimated using the GTCH method, which includes the error variances of the three datasets and the error covariance between the datasets; second, it is solved using the GLS method to obtain the optimal fusion weight vector. For each spatial pixel, the optimal fusion weights are combined, i.e., the weight vectors are weighted together to obtain the error variance of the fusion result for that pixel. This error variance is the fusion uncertainty, and its magnitude comprehensively reflects the error contribution of each data source and their interrelationship.

[0078] The resulting global pixel-by-pixel fusion uncertainty field is output synchronously with the global land water storage change fusion dataset, forming a complete credibility specification. Users can judge the reliability of the fusion results in different regions based on the magnitude of the uncertainty field: in plains and tropical regions with low uncertainty, the fusion results are highly reliable and can be directly used for refined hydrological analysis; in high-latitude icy and snowy regions, near coastlines, and areas with drastic topographic changes, where uncertainty is high, users are advised to perform cross-validation with other independent observations or adopt a more conservative interpretation strategy. This quantitative uncertainty characterization mechanism significantly improves the transparency and credibility of the fusion dataset in practical applications such as water resource management and drought monitoring.

[0079] Based on three sets of GRACE / GRACE-FO Mascon products interpolated using SSII, this paper systematically presents their error structure, uncertainty estimation, and spatiotemporal characteristics of TWSA (Time-Wide Average Scale) of terrestrial water storage changes after multi-source fusion. The fusion results are independently validated using GLDAS-Noah model data. Specific validation indicators include three aspects: first, the time series correlation coefficient, used to assess the synchronization between the fusion results and GLDAS-Noah simulations in interannual fluctuations and seasonal cycles; second, the root mean square error, used to quantify the absolute deviations of the fusion results and model simulations in amplitude and phase; and third, the long-term trend consistency index, which assesses the ability of the fusion dataset to capture long-term water storage surplus and deficit signals by comparing the slopes of the linear trends of terrestrial water storage changes between the two datasets from 2002 to 2024. These three indicators systematically evaluate the accuracy and reliability of the fusion dataset from three time scales: short-term fluctuations, seasonal variations, and long-term trends, providing independent accuracy validation support for the application of the dataset in global water cycle research and regional water resource assessment.

[0080] The selection of effective components was achieved using the htest array. Fourier transform and power spectral density (PSD) were used to calculate the cumulative energy of each principal component, and the cumulative energy was then used to determine which principal components were effective. Only principal components with a cumulative energy exceeding 90% were considered effective; others, such as high-frequency noise, were excluded. The overall weak correlation coefficients obtained from the GTCH estimation are shown in Table 2.

[0081] Table 2. Global pixel cropping ratio of 0 The fused TWSA exhibits greater consistency and coherence in spatial distribution, resulting in reduced strip noise, more continuous spatial transitions, and the elimination of single-product dominance in high-noise regions. The overall uncertainty estimated by the GTCH framework is lower globally than that of any single product. In most regions, the standard deviation of the fused TWSA is reduced by an average of approximately 15%–30% compared to the individual results of CSR, JPL, and GSFC. The spatial distribution of uncertainty correlates well with the weight distribution; regions with relatively balanced weights or dominated by multiple products exhibit the lowest fused uncertainty. Conversely, regions where a single product holds a dominant position show lower uncertainty than single-product uncertainty, but the improvement is relatively limited. This indicates that GLS fusion not only relies on a single high-quality product but also effectively mitigates random errors and structural noise through multi-source complementarity. On a temporal scale, the uncertainty of the fused TWSA remains stable during the transition from GRACE to GRACE-FO, without significant jumps, demonstrating the good stability of SSA filling and GLS fusion during the cross-mission transition period.

[0082] For some pixels, due to extremely weak signals or excessively high condition numbers in the GTCH covariance matrix, GLS fusion may not be able to solve stably. To address these pixels, a fallback strategy is introduced: when the GLS conditions are not met or the solution is unstable, a degenerate fusion scheme, such as an equal-weighted averaging strategy, is employed.

[0083] The fallback strategy primarily targets pixels located in extremely weak signal areas, high noise areas, and a small number of land-sea transition zones. These pixels represent a relatively small spatial proportion, less than 3% of the global land area. From a time-series perspective, the TWSA variation of these pixels is inherently low, and their impact on overall watershed or regional scale analyses is limited.

[0084] Comparative analysis shows that the fallback strategy for processing pixel time series outperforms direct removal or simple interpolation methods in terms of continuity and stability, while avoiding non-physical oscillations that may be introduced by GLS solution failures. Therefore, this strategy improves the usability and robustness of the global-scale fused TWSA product while ensuring overall data integrity.

[0085] This invention constructs a global land water storage change (TWSA) dataset based on SSII interpolation, GTCH-KKT uncertainty estimation, and covariance-based GLS weighted fusion. Results show that the proposed method achieves significant improvements in temporal continuity, error structure characterization, and cross-institutional Mascon product fusion. Regarding time series completion, the SSII interpolation method effectively recovers the systematically missing months during the handover between the GRACE and GRACE-FO missions, maintaining the continuity of trend and seasonal components, and sufficiently smoothing the non-physical jumps in the original time series caused by observation interruptions and monthly missing data. Unlike traditional harmonic fitting or machine learning interpolation, SSA relies on principal component decomposition to automatically distinguish between trends, seasonality, and high-frequency noise, thus maintaining the physical interpretability of the time series without relying on external data.

[0086] In terms of uncertainty estimation, the GTCH-KKT method successfully characterized the error differences among the three Mascon products (CSR, JPL, and GSFC) and revealed significant non-uniformity in their error spatial structure. The error variances of the three sets of products exhibited a spatial distribution with significant regional differences. Errors were relatively low in tropical and temperate plains, while errors were significantly higher in high-latitude icy regions, near coastlines, and areas with drastic topographic changes. Simultaneously, the correlation coefficients estimated by GTCH were generally weak, with approximately 75% of pixels being pruned to 0 due to correlation below the threshold. This indicates that the error of Mascon products from different institutions is dominated by different system regularization strategies, resulting in limited error coupling. This phenomenon differs from previous studies that assumed strong error correlation among GRACE products from different institutions, and is more consistent with the reality of complex and regionally specific GRACE error sources. Threshold pruning helps reduce the interference of unstable covariance on the GLS weight solution, making Σ more numerically stable.

[0087] The advantages of multi-source information integration are further demonstrated based on the GLS fusion results. Since the fusion weights are jointly determined by the error variance and covariance, GLS can adaptively select the most reliable data source for different regions. In regions with significant covariance, the full GLS can fully utilize covariance information and improve the effectiveness of weight allocation; while in regions with extremely weak correlation, GLS automatically degenerates into inverse variance weighting, thereby avoiding weight bias introduced by unstable covariance. Compared with traditional arithmetic averages or simple weighting, the GLS fusion of this invention not only significantly improves the spatial consistency of TWSA but also effectively suppresses noise anomalies in single products, making trend terms and seasonality more consistent with physical expectations. The uncertainty of the fusion results is generally lower than that of single products, with improvements exceeding 50% in some high-error regions, indicating that the constructed covariance-driven fusion framework has good robustness.

[0088] This invention also found that the fallback strategy accounts for only 2.63% of the pixels, indicating that the solution process of GTCH and GLS is highly stable globally, with only a few extreme regions (such as the polar edges and high-leakage regions) experiencing issues due to Σ. The result degenerates into an equal-weighted average, further validating the robustness of the proposed method on large-scale grid data. However, this invention still has certain limitations. GTCH is sensitive to time series length, and may underestimate covariance in early GRACE data segments; covariance terms in weakly correlated regions may be over-pruned, thus ignoring potential systematic error information; the fusion process does not incorporate hydrological model priors, therefore it cannot eliminate common systematic errors in Mascon products, such as GIA model errors. Future research needs to further incorporate hydrological models, GNSS observations, or multi-satellite joint inversion into the error modeling framework to improve the physical interpretation and overall accuracy of the fusion results.

[0089] Based on the same inventive concept, this invention also provides a dataset construction system for terrestrial water storage changes, comprising: The data acquisition module is used to acquire a global land and ocean quality change dataset over a continuous multi-year period based on a monthly time resolution. The global land and ocean quality change dataset includes the CSR dataset, the JPL dataset, and the GSFC dataset.

[0090] The interpolation module is used to segment the CSR, JPL, and GSFC datasets respectively, and construct trajectory matrices based on the subsequences of each dataset segment. Singular value decomposition is performed on each trajectory matrix to obtain the temporal principal component matrices of different modes, the singular value matrices of different mode energy magnitudes, and the empirical orthogonal function matrices describing the spatial structure of each mode. A rank-one matrix is ​​constructed using the temporal principal component matrix, singular value matrix, and empirical orthogonal function matrix, and the reconstructed components of each mode are obtained based on the rank-one matrix. The reconstructed components are iteratively interpolated and modal-filtered, and then linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset.

[0091] The fusion module is used to estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices. Based on the error variance and error covariance, an error covariance matrix is ​​constructed. Based on the error covariance matrix, the optimal fusion weight corresponding to each interpolated spatiotemporal data matrix is ​​obtained pixel by pixel. The optimal fusion weight is used to weight and fuse each interpolated spatiotemporal data matrix to obtain a fused dataset of global land water storage changes.

[0092] The present invention also provides a computer-readable storage medium storing a computer program that can be used to execute the above-described method for constructing a dataset of changes in terrestrial water storage.

[0093] For specific limitations on the computational system for constructing the dataset of terrestrial water storage changes, please refer to the limitations on the dataset construction method for terrestrial water storage changes mentioned above, which will not be repeated here. Each module in the aforementioned dataset construction system for terrestrial water storage changes can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device, or stored in the memory of a computer device as software, so that the processor can call and execute the corresponding operations of each module.

[0094] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification. Furthermore, the above embodiments only illustrate several implementation methods of this application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make several modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims.

Claims

1. A method for constructing a dataset of changes in terrestrial water storage, characterized in that, Includes the following steps: Based on a monthly time resolution, a global land and ocean quality change dataset was obtained over a continuous multi-year period. This dataset includes the CSR dataset, the JPL dataset, and the GSFC dataset. The CSR, JPL, and GSFC datasets were segmented, and trajectory matrices were constructed based on the subsequences of each dataset. Singular value decomposition was performed on each trajectory matrix to obtain the temporal principal component matrices of different modes, the singular value matrices of different mode energy magnitudes, and the empirical orthogonal function matrices describing the spatial structure of each mode. A rank-one matrix was constructed using the temporal principal component matrix, singular value matrix, and empirical orthogonal function matrix, and the reconstructed components of each mode were obtained based on the rank-one matrix. The reconstructed components were iteratively interpolated and mode-filtered, and then linearly superimposed to obtain the interpolated spatiotemporal data matrix of each dataset. Estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices. Construct an error covariance matrix based on the error variance and error covariance. Based on the error covariance matrix, obtain the optimal fusion weight corresponding to each interpolated spatiotemporal data matrix pixel by pixel. By using the optimal fusion weights to weight and fuse each interpolated spatiotemporal data matrix, a fused dataset of global land water storage changes is obtained. Also includes: The generalized triangular hat method is used to estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices, specifically including: Each time, one of the interpolated spatiotemporal data matrices from the CSR dataset, JPL dataset, and GSFC dataset is selected as a reference benchmark. The difference sequence of the interpolated spatiotemporal data matrices of the other two datasets relative to the current reference benchmark is calculated. Based on the covariance matrix of each difference sequence, the error variance and error covariance of each dataset are solved simultaneously, and three preliminary error covariance estimation matrices are constructed. The median is extracted from the matrix elements representing the same error variance or error covariance in the three preliminary error covariance estimation matrices, and all medians are fused to generate a preliminary fused error covariance matrix; The preliminary fusion error covariance matrix is ​​optimized using the off-diagonal energy minimization method, and the KKT hard constraint of matrix positive definiteness is introduced. The error correlation coefficient in the optimized matrix is ​​thresholded to obtain the error variance and error covariance of different spatial pixels. Based on the aforementioned error covariance matrix, the optimal fusion weights corresponding to each interpolated spatiotemporal data matrix are obtained pixel-by-pixel using the generalized least squares method. Specifically, this includes: determining whether the error covariance matrix contains valid error covariance; if the error covariances are correlated, the optimal fusion weights are calculated directly based on the inverse of the error covariance matrix using the generalized least squares method; if not, the weights are assigned according to the reciprocal of the error variance of each interpolated spatiotemporal data matrix in the corresponding pixel, with a higher weight for smaller error variances; if neither of the above two cases applies, equal weights are assigned to the three interpolated spatiotemporal data matrices CSR, JPL, and GSFC, and an equal-weighted average is performed.

2. The method for constructing a dataset of terrestrial water storage changes according to claim 1, characterized in that, The reconstructed components are iteratively interpolated and modal filtered, then linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset, specifically including: Calculate the relative root mean square error of adjacent reconstructed components at the missing position. If the relative root mean square error is greater than a preset convergence threshold, increase the number of modes, reconstruct the trajectory matrix and perform singular value decomposition, and continue iterating. If the relative root mean square error is less than or equal to the preset convergence threshold, the iteration is complete and the mode selection step is entered. For all reconstructed components obtained after iteration, calculate the cumulative power spectrum distribution function of the empirical orthogonal function matrix corresponding to each reconstructed component; if the value of the cumulative power spectrum distribution function at the preset characteristic frequency is greater than the set threshold, then the current reconstructed component is determined to be an effective low-frequency mode and is retained; otherwise, it is determined to be a high-frequency noise mode and is discarded. All the retained reconstructed components are linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset.

3. The method for constructing a dataset of terrestrial water storage changes according to claim 1, characterized in that, The process of segmenting the CSR, JPL, and GSFC datasets involves the following steps: First, the missing length of each dataset is determined, resulting in two scenarios: odd / even month missing data and consecutive 11-month missing data. For odd / even month missing data, singular spectrum iterative interpolation is used, employing an empirically set fixed window length and principal component count to interpolate and fill in the missing data in each dataset. For consecutive 11-month missing data, a parameter optimization approach is used. Within the complete timeframe of each dataset, cross-validation is employed to test the interpolation accuracy of different MK combinations for simulated missing data segments. The MK combination with the smallest root mean square error is selected as the region-adaptive fixed window length and principal component count for segmenting each dataset. Through iterative reconstruction and effective modality selection, the retained reconstructed components are linearly superimposed to output a segmented, gap-free continuous sequence. Indicates the window length. This indicates the number of principal components used for reconstruction.

4. The method for constructing a dataset of terrestrial water storage changes according to claim 1, characterized in that, After obtaining the global land water storage change fusion dataset, the method further includes: obtaining a pixel-by-pixel fusion uncertainty field that is compatible with the global land water storage change fusion dataset; the uncertainty field is obtained based on the error covariance matrix and the optimal fusion weights, and a credibility specification is provided for the global land water storage change fusion dataset based on the uncertainty field.

5. The method for constructing a dataset of terrestrial water storage changes according to claim 1, characterized in that, Also includes: The simulation results of the preset independent GLDAS-Noah land surface model are obtained. The accuracy and reliability of the fused dataset are evaluated by calculating the time series correlation coefficient, root mean square error and trend consistency index between the fused dataset of global land water storage changes and the simulation results.

6. A system for constructing a dataset of terrestrial water storage changes as described in claim 1, characterized in that, include: The data acquisition module is used to acquire a global land and ocean quality change dataset over a continuous multi-year period based on a monthly time resolution. This global land and ocean quality change dataset includes the CSR dataset, the JPL dataset, and the GSFC dataset. The interpolation module is used to segment the CSR, JPL, and GSFC datasets respectively, and construct trajectory matrices based on the subsequences of each dataset segment. Singular value decomposition is performed on each trajectory matrix to obtain the temporal principal component matrices of different modes, the singular value matrices of different mode energy magnitudes, and the empirical orthogonal function matrices describing the spatial structure of each mode. A rank-one matrix is ​​constructed using the temporal principal component matrix, singular value matrix, and empirical orthogonal function matrix, and the reconstructed components of each mode are obtained based on the rank-one matrix. The reconstructed components are iteratively interpolated and modal-filtered, and then linearly superimposed to obtain the interpolated spatiotemporal data matrix for each dataset. The fusion module is used to estimate the error variance of each interpolated spatiotemporal data matrix at different spatial pixels and the error covariance between different interpolated spatiotemporal data matrices. Based on the error variance and error covariance, an error covariance matrix is ​​constructed. Based on the error covariance matrix, the optimal fusion weight corresponding to each interpolated spatiotemporal data matrix is ​​obtained pixel by pixel. The optimal fusion weight is used to weight and fuse each interpolated spatiotemporal data matrix to obtain a fused dataset of global land water storage changes.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is loaded by the processor, it is able to perform the steps of the method according to any one of claims 1 to 5.

Images