A method for evaluating the applicability of groundwater storage anomaly inversion products

Abstract

The present invention discloses a method for evaluating the applicability of groundwater reserve anomaly inversion products, comprising the following steps: Step S1: generating a regional groundwater reserve anomaly inversion product by inversion based on the GRACE / GRACE-FO total water reserve anomaly product and surface water reserve and soil water reserve anomaly products simulated by a hydrological model; calculating the cross-correlation coefficient, phase difference coefficient of variation, and amplitude ratio coefficient of variation between a pixel-scale measured water level time series and a groundwater reserve anomaly time series; and evaluating the applicability of the groundwater reserve anomaly inversion product from multiple perspectives, including trend, phase, and amplitude. The present invention comprehensively evaluates the applicability of the groundwater reserve anomaly inversion product using three evaluation indicators, namely the cross-correlation coefficient, phase difference coefficient of variation, and amplitude ratio coefficient of variation, to improve the accuracy of the evaluation result.

Description

Technical Field

[0001] The present invention relates to the technical field of groundwater reserve remote sensing monitoring, and in particular to a method for evaluating the applicability of groundwater reserve anomaly inversion products. Background Art

[0002] Groundwater is the world's most widespread freshwater resource. Compared to surface water, groundwater is less polluted and safer to drink, making it a crucial source of freshwater for agricultural irrigation, industrial production, and urban life in many regions of the world, particularly in semi-arid and densely populated areas. Therefore, effectively monitoring the dynamic changes in groundwater reserves is beneficial for optimizing the development and protection of groundwater resources, and is of great significance for maintaining the ecological environment and ensuring sustainable social and economic development.

[0003] Groundwater changes are typically monitored through groundwater observation wells. However, due to limitations such as time-consuming and labor-intensive well deployment and maintenance, insufficient well density, and limited monitoring range, groundwater well observations often fail to fully reflect basin-scale spatiotemporal groundwater variations. The Gravity Recovery and Climate Experiment (GRACE) and its subsequent satellite, GRACE Follow-On (GRACE-FO), launched jointly by NASA and the German Aerospace Center in 2002 and 2018, provided effective approaches for large-scale monitoring of total water storage anomalies (TWSA). By deducting surface water storage anomalies and soil moisture anomalies derived from hydrological model simulations, including anomalies in snow water equivalent, lake and reservoir storage, runoff, and canopy water content, from the GRACE / GRACE-FO inverted total water storage anomalies, groundwater storage anomalies can be successfully isolated. However, because different institutions, such as NASA's Jet Propulsion Laboratory (JPL), the University of Texas Center for Space Research (CSR), and NASA's Goddard Space Flight Center (GSFC), employ different processing strategies and computational models when processing raw GRACE gravity satellite observation data, the applicability of total water storage anomalies derived from the spherical harmonic coefficient (SH) products or Mascon products released by these institutions, such as GRACE / GRACE-FO, varies across regions. Furthermore, different assumptions and conceptualizations, coupled with varying model structures and parameterization methods, result in varying applicability of different hydrological models for simulating surface water storage anomalies and soil moisture anomalies in different regions. Consequently, significant differences in groundwater storage anomalies are observed when combining different GRACE / GRACE-FO products with products simulated by different hydrological models. Therefore, evaluating the applicability of combining different GRACE / GRACE-FO products and hydrological model simulation products to invert groundwater storage anomalies in different basins is crucial for conducting research on the temporal and spatial variation analysis of groundwater storage.

[0004] Currently, the applicability evaluation of groundwater storage anomaly inversion products mainly uses groundwater well observed water level data as measured values. By calculating the correlation coefficient between two time series: groundwater storage anomalies and groundwater well observed water level changes, the applicability of different groundwater storage anomaly inversion products is evaluated based on the size of the correlation coefficient. However, the correlation coefficient mainly represents the trend consistency between the two time series, while groundwater storage changes often show certain intra-annual periodic variation characteristics. Evaluating the applicability of groundwater storage anomaly inversion products using only the correlation coefficient ignores the ability of groundwater storage anomaly inversion products to represent the characteristics of intra-annual periodic changes (phase and amplitude) of groundwater, and the applicability evaluation is not comprehensive and accurate. Summary of the Invention

[0005] In order to solve the problems existing in the prior art, the present invention provides a method for evaluating the applicability of groundwater reserve anomaly inversion products. The present invention uses three evaluation indicators, namely the cross-correlation coefficient, the phase difference variation coefficient and the amplitude ratio variation coefficient, to comprehensively evaluate the applicability of groundwater reserve anomaly inversion products to improve the accuracy of the evaluation results.

[0006] The technical solutions of the present invention are as follows:

[0007] A method for evaluating the applicability of groundwater storage anomaly inversion products comprises the following steps:

[0008] Step S1: Generate regional groundwater storage anomaly inversion products based on the GRACE / GRACE-FO total water storage anomaly products and the surface water storage and soil water storage anomaly products simulated by the hydrological model. Calculate the cross-correlation coefficient, phase difference coefficient of variation, and amplitude ratio coefficient of variation between the pixel-scale measured water level time series and the groundwater storage anomaly time series. Evaluate the applicability of the groundwater storage anomaly inversion products from multiple perspectives: trend, phase, and amplitude.

[0009] Step S2: The evaluation sample set is composed of the mutual correlation coefficient, phase difference variation coefficient, and amplitude ratio variation coefficient of all pixels in the evaluation area. The weights of the mutual correlation coefficient, phase difference variation coefficient, and amplitude ratio variation coefficient are determined using the principal component analysis method, and the comprehensive applicability index of the pixel-scale groundwater storage anomaly inversion product is weighted and calculated.

[0010] Step S3: Based on the comprehensive applicability index results of the groundwater storage anomaly inversion product, the applicability levels are graded and evaluated to evaluate the applicability of the regional groundwater storage anomaly inversion product for the analysis of spatiotemporal changes in groundwater storage.

[0011] Preferably, step S1 includes the following sub-steps:

[0012] Sub-step S11: generating regional groundwater storage anomaly inversion products based on GRACE / GRACE-FO and GLDAS total water storage anomaly products and hydrological model simulation surface water storage and soil water storage anomaly products;

[0013] Sub-step S12: obtaining and pre-processing the measured groundwater level data in the area;

[0014] Sub-step S13: Multi-perspective applicability evaluation of groundwater storage anomaly inversion products in the area based on measured groundwater levels.

[0015] Preferably, the expression of groundwater reserve anomaly in sub-step S11 is:

[0016] GWSA=TWSA1-SWSA-SMSA

[0017] Where GWSA represents groundwater storage anomaly, TWAS represents total water storage anomaly, SWSA represents surface water storage anomaly, and SMSA represents soil water storage anomaly.

[0018] Preferably, the multi-view suitability evaluation indexes of sub-step S13 include the cross-correlation coefficient, the phase difference variation coefficient, and the amplitude ratio variation coefficient;

[0019] The formula for calculating the mutual correlation coefficient is:

[0020]

[0021] Where: σ x , σ y are the mean square errors of the time series x(t) and y(t), respectively; are the means of x(t) and y(t), respectively; k is the lag, which indicates the degree of offset between the two time series; C xy (k) and r xy (k) are the cross-covariance and cross-correlation coefficient of time series x(t) and y(t) at time lag k, respectively; n is the length of the time series; the cross-correlation value range is -1 to 1. The closer the cross-correlation coefficient of two time series is to 1, the more similar the changing trends of the two time series are.

[0022] The calculation formula of the phase difference variation coefficient is:

[0023]

[0024] Where: σ PCE is the standard deviation of PCE; μ PCE is the average value of PCE; the closer PCV is to 0, the higher the phase consistency of the two time series;

[0025] The calculation formula of the coefficient of variation of the amplitude ratio is:

[0026]

[0027] Where: σ AR is the standard deviation of AR; μ AR is the average value of AR; the closer ACV is to 0, the higher the phase consistency of the two time series.

[0028] Preferably, step S2 uses principal component analysis to determine the weights of the mutual correlation coefficient, the phase difference variation coefficient, and the amplitude ratio variation coefficient, and weightedly calculates the comprehensive applicability index of the pixel-scale groundwater storage anomaly inversion product. The calculation formula of the comprehensive applicability index of the groundwater storage anomaly inversion product is as follows:

[0029]

[0030] Where: PCV p is the pth principal component; α p is the contribution rate of the corresponding principal component.

[0031] The comprehensive applicability index of the groundwater storage anomaly inversion product is normalized and the calculation formula is:

[0032]

[0033] Where: GSEI norm is the normalized comprehensive applicability index value; GSEI, GSEL max and GSEL min They represent the original value, the maximum value and the minimum value of the improved groundwater storage anomaly inversion product applicability index in all pixels.

[0034] Preferably, the grading and classification of the applicability levels in step S3 is specifically based on the comprehensive applicability index result of the groundwater storage anomaly inversion product, and the applicability index is graded and divided into five levels with 0.2 as an interval: the first-level numerical interval is [0, 0.2], indicating that the applicability of the groundwater storage anomaly inversion product is "poor"; the second-level numerical interval is (0.2, 0.4], indicating that the applicability is "poor"; the third-level numerical interval is (0.4, 0.6], indicating that the applicability is "general"; the fourth-level numerical interval is (0.6, 0.8], indicating that the applicability is "good"; the fifth-level numerical interval is (0.8, 1], indicating that the applicability is "excellent". According to the above grading results, the applicability of the groundwater storage anomaly inversion product is comprehensively evaluated.

[0035] Compared with the prior art, the technical effects of the present invention are as follows:

[0036] 1. The present invention innovatively proposes a method for evaluating the applicability of groundwater reserve anomaly inversion products. Existing methods for evaluating the applicability of groundwater reserve anomaly inversion products have many shortcomings. For example, the evaluation indicators are single. Existing indicators and methods only focus on the trend characteristics of groundwater reserve anomaly inversion products, ignoring the phase and amplitude characteristics involved in the periodic changes of groundwater within the year, resulting in incomplete and inaccurate applicability evaluation. In view of this, the present invention not only measures the consistency of the overall trend of the time series by calculating the correlation coefficient between the measured water level changes and groundwater reserve anomalies, but also uses the phase difference coefficient of variation and the amplitude ratio coefficient of variation to evaluate the applicability of groundwater reserve anomaly inversion products from the perspectives of phase and amplitude. Three representative indicators are used to evaluate the applicability of groundwater reserve anomaly inversion products from multiple perspectives such as trend, phase and amplitude, and the evaluation results are more comprehensive and accurate.

[0037] 2. This invention introduces principal component analysis (PCA) to comprehensively evaluate multiple time series consistency indicators. The weighting of each indicator is determined entirely based on the data's inherent properties. When the cumulative contribution rate of the principal components in PCA reaches 85% or above, the eigenvalue contribution rate is used as the weighted sum to calculate the suitability index. This approach, while ensuring data accuracy, encompasses a greater degree of information in the original data, making the evaluation results more realistic. This invention largely avoids interference from human factors and improves the objectivity of the applicability evaluation of groundwater storage anomaly inversion products. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 It is a flow chart of a method for evaluating the applicability of groundwater reserve anomaly inversion products according to an embodiment of the present invention.

[0039] Figure 2 This is a distribution map of groundwater monitoring wells in the Lazi to Pai Township watershed according to an embodiment of the present invention.

[0040] Figure 3 Schematic diagram of the spatial distribution of evaluation results of different applicability evaluation indicators of the groundwater reserve anomaly inversion product according to an embodiment of the present invention.

[0041] Figure 3 (a) is the spatial distribution diagram of the cross-correlation coefficient CC results.

[0042] Figure 3 (b) is the spatial distribution diagram of the phase difference coefficient of variation PCV results.

[0043] Figure 3 (c) is the spatial distribution diagram of the amplitude ratio coefficient of variation ACV results.

[0044] Figure 3 (d) is the spatial distribution map of the groundwater storage anomaly inversion product applicability index GSEI results.

[0045] Figure 4 This is a time series diagram of groundwater storage anomalies and measured groundwater level changes in pixels 3-11 and pixels 6-9. DETAILED DESCRIPTION

[0046] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0047] The specific embodiments of the present invention are described below to facilitate understanding of the present invention by those skilled in the art. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the appended claims, these changes are obvious, and all inventions and creations utilizing the concepts of the present invention are protected.

[0048] A method for evaluating the applicability of groundwater storage anomaly inversion products comprises the following steps:

[0049] S1. Generate regional groundwater storage anomaly inversion products based on the GRACE / GRACE-FO total water storage anomaly products and the surface water storage and soil water storage anomaly products simulated by the hydrological model. Calculate the cross-correlation coefficient, phase difference coefficient of variation, and amplitude ratio coefficient of variation between the pixel-scale measured water level time series and the groundwater storage anomaly time series. Evaluate the applicability of the groundwater storage anomaly inversion products from multiple perspectives: trend, phase, and amplitude.

[0050] S2. The evaluation sample set is composed of the mutual correlation coefficient, phase difference variation coefficient and amplitude ratio variation coefficient of all pixels in the evaluation area. The weights of the mutual correlation coefficient, phase difference variation coefficient and amplitude ratio variation coefficient are determined by principal component analysis, and the comprehensive applicability index of the pixel-scale groundwater storage anomaly inversion product is calculated by weighted calculation.

[0051] S3. Based on the comprehensive applicability index results of groundwater storage anomaly inversion products, classify the applicability levels and evaluate the applicability of regional groundwater storage anomaly inversion products for analyzing the spatiotemporal changes of groundwater storage.

[0052] Step S1 of this embodiment specifically includes the following sub-steps:

[0053] S11. Production of regional groundwater storage anomaly inversion products based on the GRACE / GRACE-FO and GLDAS hydrological models.

[0054] GRACE / GRACE-FO products include spherical harmonic coefficient (SH) and Mascon products released by three different institutions: CSR, JPL, and GFZ. GLDAS hydrological models include GLDAS-Noah, GLDAS-VIC, and GLDAS-CLSM. GRACE / GRACE-FO products are used to invert total water storage anomaly (TWSA), which is primarily influenced by surface water storage anomaly (SWSA), soil water storage anomaly (SMSA), and groundwater storage anomaly (GWSA). Surface water storage anomaly (SWSA) can be considered to be composed of snow water equivalent anomaly (SWEA), vegetation canopy water anomaly (CWSA), surface runoff anomaly (SRSA), and lake and reservoir water storage anomaly (LWSA). Therefore, groundwater storage anomaly (GWSA) can be obtained using the following formula:

[0055] GWSA=TWSA-SWSA-SMSA (1)

[0056] TWSA is derived from the inversion of GRACE / GRACE-FO products, while SWSA and SMSA are derived from the GLDAS hydrological model. The unit of each water storage component anomaly is equivalent water height, cm.

[0057] S12. Acquisition and preprocessing of regional measured groundwater level data.

[0058] Groundwater level data are collected from the National Groundwater Monitoring Project. The observation interval is one month, and the water level unit is meter. The measured groundwater level data were quality-checked, and outliers and stations with significant sequence missing values ​​were removed. The remaining station data were then averaged to obtain monthly changes. Based on the preprocessed measured groundwater level change data, the mean of the measured groundwater level change in each pixel was calculated to obtain a pixel-scale time series of measured groundwater level changes.

[0059] S13. Multi-perspective applicability evaluation of regional groundwater storage anomaly inversion products based on measured groundwater levels.

[0060] In step S13 of this embodiment, the multi-view suitability evaluation indicators include the cross-correlation coefficient, the phase difference variation coefficient, and the amplitude ratio variation coefficient.

[0061] The formula for calculating the mutual correlation coefficient is as follows:

[0062]

[0063] Where: σ x , σ y are the mean square errors of the time series x(t) and y(t), respectively; are the means of x(t) and y(t), respectively; k is the lag, which indicates the degree of offset between the two time series; Cxy (k) and r xy (k) are the cross-covariance and cross-correlation coefficient of the time series x(t) and y(t), respectively, at lag k; n is the length of the time series. The cross-correlation coefficient ranges from -1 to 1. The closer the cross-correlation coefficient is to 1, the more similar the trends of the two time series are.

[0064] In order to facilitate the comparison between different indicators, they are normalized according to the characteristics of the numerical domain of the mutual relationship to eliminate the dimension difference. The normalization formula is as follows:

[0065]

[0066] Where: CC is the normalized maximum cross-correlation coefficient; a is the upper limit of the cross-correlation range; b is the lower limit of the cross-correlation range; is the correlation coefficient of the two time series under the optimal time lag, that is, the optimal cross-correlation coefficient, indicating that the cross-correlation coefficient of the two time series is the largest under the optimal time lag k0. The optimal time lag can be used to determine the time lag between the two time series and to align the time series.

[0067] The optimal time lag k0>0 means that sequence B lags behind sequence A by k0 time units. In this case, sequence B needs to be adjusted and moved forward by k0 time units to align with sequence A. In the specific operation, due to the limited time series of measured groundwater level changes and groundwater storage anomalies, truncation is used to discard the last k0 data points of sequence B to obtain a new sequence In order to keep the lengths of the two time series consistent, it is necessary to discard the redundant data points of sequence A and obtain a new sequence This aligns sequence A' and sequence B'.

[0068] If the optimal time lag k0 < 0, it means that sequence B is ahead of sequence A by |k0| time units. The specific operation is opposite to that when k0 < 0.

[0069] When the optimal time lag k0 = 0, it means that the two time series are already aligned in time and no additional operations are required.

[0070] The phase difference coefficient of variation was calculated based on the aligned time series.

[0071] First, the original time series is subjected to Hilbert transform. This method can construct an analytical signal corresponding to the original signal and has a good effect on extracting parameters such as the instantaneous amplitude, instantaneous phase and instantaneous frequency of the signal.

[0072]

[0073] Where: Y i (t) is the analytical signal corresponding to the original signal; is the imaginary part; Represents a given time series X i The Hilbert transform of (t) is calculated as follows:

[0074]

[0075] Where: PV represents the Cauchy principal value. Formula (4) can be expressed in polar coordinate form as:

[0076]

[0077] The time signal X can be determined by formula (6): i Amplitude A of (t) i (t) and instantaneous phase φ i (t), where the instantaneous phase is calculated as follows:

[0078]

[0079] Get two time signals X i (t), X j The instantaneous phase φ of (t) i (t) and φ j (t) After that, the phase difference of the two time series can be calculated:

[0080] Δφ(t)=φ i (t)-φ j (t) (8)

[0081] In order to more accurately reflect the phase difference characteristics and facilitate processing and calculation, the phase difference is projected onto the complex exponential plane:

[0082] P CE =e iΔφ (t) (9)

[0083] Thus, the phase difference coefficient of variation of the two time signals is calculated:

[0084]

[0085] Where: σ PCE is the standard deviation of PCE; μ PCE is the average value of PCE. The closer PCV is to 0, the higher the phase consistency of the two time series.

[0086] According to the range characteristics of the phase difference coefficient of variation, reverse normalization is performed to reduce the impact of outliers on the overall data distribution, ensure that the true difference between the original coefficients of variation can be effectively reflected, and eliminate the dimensional differences between different indicators to facilitate comparison between different indicators. The reverse normalization formula is as follows:

[0087]

[0088] The coefficient of variation of the amplitude ratio was calculated from the aligned time series.

[0089] Calculate the time signal X by formula (6) i Amplitude A of (t) i (t):

[0090]

[0091] Get two time signals X i (t), X j The instantaneous amplitude A of (t) i (t) and A j (t) After that, the amplitude spectrum is normalized to calculate the amplitude ratio of the two time series:

[0092]

[0093] Thus, the coefficient of variation of the amplitude ratio of the two time signals can be calculated:

[0094]

[0095] Where: σ AR is the standard deviation of AR; μ AR is the average value of AR. The closer ACV is to 0, the higher the phase consistency of the two time series.

[0096] Based on the range characteristics of the coefficient of variation of the amplitude ratio, reverse normalization is performed to reduce the impact of outliers on the overall data distribution, ensure that the true differences between the original coefficients of variation can be effectively reflected, and eliminate the dimensional differences between different indicators to facilitate comparison between different indicators. The reverse normalization formula is as follows:

[0097]

[0098] Furthermore, the step S2 uses the principal component analysis method to integrate the multi-perspective applicability evaluation indicators to calculate the comprehensive applicability index of the groundwater storage anomaly inversion product, and the calculation formula is as follows:

[0099]

[0100] Where: PCV pis the pth principal component; α p is the contribution rate of the corresponding principal component.

[0101] The comprehensive evaluation results of the applicability of the groundwater storage anomaly inversion products are normalized, and the calculation formula is as follows:

[0102]

[0103] Where: GSEI norm is the normalized comprehensive applicability index value; GSEI, GSEI max and GSEI min They represent the original value, the maximum value and the minimum value of the improved groundwater storage anomaly inversion product applicability index in all pixels.

[0104] In step S3 of this implementation plan: for any study area, based on the comprehensive applicability index results of the groundwater storage anomaly inversion product, the applicability index is graded and divided into five levels with 0.2 as an interval: the first-level numerical interval is [0, 0.2], indicating that the applicability of the groundwater storage anomaly inversion product is "poor"; the second-level numerical interval is (0.2, 0.4], indicating that the applicability is "poor"; the third-level numerical interval is (0.4, 0.6], indicating that the applicability is "general"; the fourth-level numerical interval is (0.6, 0.8], indicating that the applicability is "good"; the fifth-level numerical interval is (0.8, 1], indicating that the applicability is "excellent". Based on the above grading results, the applicability of the groundwater storage anomaly inversion product is comprehensively evaluated.

[0105] When this implementation plan is implemented,

[0106] like Figure 1 As shown in FIG, a method for evaluating the applicability of groundwater storage anomaly inversion products includes the following steps:

[0107] S1. Generate regional groundwater storage anomaly inversion products based on the GRACE / GRACE-FO total water storage anomaly products and surface water storage anomaly products simulated by hydrological models. Calculate the correlation coefficient, phase difference coefficient of variation, and amplitude ratio coefficient of variation between the pixel-scale measured water level time series and the groundwater storage anomaly time series. Evaluate the applicability of the groundwater storage anomaly inversion products from multiple perspectives, including trend, phase, and amplitude.

[0108] S11. Production of regional groundwater storage anomaly inversion products based on the GRACE / GRACE-FO and GLDAS hydrological models.

[0109] Example 1: The Lazi to Pai Township watershed was used as the study area. The JPLRL06.1_v03 Mascon product released by the Jet Propulsion Laboratory of the United States was obtained. The monthly soil water storage, snow water equivalent, canopy water storage, and surface runoff data provided by the Noah total surface model in the Global Land Data Assimilation System (GLDAS) model were obtained. The changes in the water storage capacity of lakes and reservoirs in the study area had a relatively small impact on the long-term total water storage anomaly and could be ignored.

[0110] JPL Mascon is used to invert total water storage anomalies (TWSA). The total water storage anomaly (TWSA) is primarily influenced by the surface water storage anomaly (SWSA) and groundwater storage anomaly (GWSA). The surface water storage anomaly (SWSA) can be considered to be composed of soil water anomaly (SMSA), snow water equivalent anomaly (SWEA), vegetation canopy water anomaly (CWSA), surface runoff anomaly (SRSA), and lake water storage anomaly (LWSA). Therefore, the groundwater storage anomaly (GWSA) can be obtained as follows:

[0111] GWSA=TWSA-SWSA-SMSA (1)

[0112] TWSA is derived from the inversion of JPLMascon products, while SWSA and SMSA are derived from the GLDAS-Noah hydrological model. The unit of each water storage component anomaly is equivalent water height, cm.

[0113] S12. Acquisition and preprocessing of regional measured groundwater level data.

[0114] The measured groundwater data comes from the National Groundwater Monitoring Project. The observation interval is one month and the water level unit is m. The measured groundwater level data was quality checked and sites with outliers and many missing sequences were removed. A total of 110 monitoring sites were selected. The data span is from July 2018 to December 2023. The distribution of monitoring wells is as follows: Figure 2 The remaining station data were de-averaged to obtain monthly change values. Based on the pre-processed measured groundwater level change data, the mean of the measured groundwater level change in each pixel was calculated to obtain the pixel-scale measured groundwater level change time series.

[0115] S13. Multi-perspective applicability evaluation of regional groundwater storage anomaly inversion products based on measured groundwater levels.

[0116] Furthermore, the multi-view suitability evaluation indexes in step S13 include the cross-correlation coefficient, the phase difference variation coefficient, and the amplitude ratio variation coefficient.

[0117] The formula for calculating the cross-correlation coefficient is as follows:

[0118]

[0119] Where: σ x , σ y are the mean square errors of the time series x(t) and y(t), respectively; are the means of x(t) and y(t), respectively; k is the lag, which indicates the degree of offset between the two time series; C xy (k) and r xy (k) are the cross-covariance and cross-correlation coefficient of the time series x(t) and y(t), respectively, at lag k; n is the length of the time series. The cross-correlation coefficient ranges from -1 to 1. The closer the cross-correlation coefficient is to 1, the more similar the trends of the two time series are.

[0120] In order to facilitate the comparison between different evaluation indicators, they are normalized according to the numerical domain characteristics of the mutual relationship to eliminate the dimension difference. The normalization formula is as follows:

[0121]

[0122] Where: CC is the normalized maximum cross-correlation coefficient; a is the upper limit of the cross-correlation range; b is the lower limit of the cross-correlation range; is the correlation coefficient of the two time series under the optimal time lag, that is, the optimal cross-correlation coefficient, indicating that the cross-correlation coefficient of the two time series is the largest under the optimal time lag k0. The optimal time lag can be used to determine the time lag between the two time series and to align the time series.

[0123] The optimal time lag k0>0 means that sequence B lags behind sequence A by k0 time units. In this case, sequence B needs to be adjusted and moved forward by k0 time units to align with sequence A. In the specific operation, due to the limited time series of measured groundwater level changes and groundwater storage anomalies, truncation is used to discard the last k0 data points of sequence B to obtain a new sequence In order to keep the lengths of the two time series consistent, it is necessary to discard the redundant data points of sequence A and obtain a new sequence This aligns sequence A' and sequence B'.

[0124] If the optimal time lag k0 < 0, it means that sequence B is ahead of sequence A by |k0| time units. The specific operation is opposite to that when k0 < 0.

[0125] When the optimal time lag k0 = 0, it means that the two time series are already aligned in time and no additional operations are required.

[0126] The phase difference coefficient of variation was calculated based on the aligned time series.

[0127] First, the original time series is subjected to Hilbert transform. This method can construct an analytical signal corresponding to the original signal and has a good effect on extracting parameters such as the instantaneous amplitude, instantaneous phase and instantaneous frequency of the signal.

[0128]

[0129] Where: Y i (t) is the analytical signal corresponding to the original signal; i is the imaginary part; Represents a given time series X i The Hilbert transform of (t) is calculated as follows:

[0130]

[0131] Where: PV represents the Cauchy principal value. Formula (4) can be expressed in polar coordinate form as:

[0132]

[0133] The time signal X can be determined by formula (6): i Amplitude A of (t) i (t) and instantaneous phase φ i (t), where the instantaneous phase is calculated as follows:

[0134]

[0135] Get two time signals X i (t), X j The instantaneous phase φ of (t) i (t) and φ j (t) After that, the phase difference of the two time series can be calculated:

[0136] Δφ(t)=φ i (t)-φ j (t) (8)

[0137] In order to more accurately reflect the phase difference characteristics and facilitate processing and calculation, the phase difference is projected onto the complex exponential plane:

[0138] PCE=e iΔφ(t) (9)

[0139] Thus, the phase difference coefficient of variation of the two time signals is calculated:

[0140]

[0141] Where: σ PCE is the standard deviation of PCE; μ PCEis the average value of PCE. The closer PCV is to 0, the higher the phase consistency of the two time series.

[0142] According to the range characteristics of the phase difference coefficient of variation, it is reverse normalized to reduce the impact of outliers on the overall data distribution, ensure that it can effectively reflect the true difference between the original coefficients of variation, and eliminate the dimensional differences between different indicators to facilitate comparison between different evaluation indicators. The reverse normalization formula is as follows:

[0143]

[0144] The coefficient of variation of the amplitude ratio was calculated from the aligned time series.

[0145] Calculate the time signal X by formula (6) i Amplitude A of (t) i (t):

[0146]

[0147] Get two time signals X i (t), X j The instantaneous amplitude A of (t) i (t) and A j (t) After that, the amplitude spectrum is normalized to calculate the amplitude ratio of the two time series:

[0148]

[0149] Thus, the coefficient of variation of the amplitude ratio of the two time signals can be calculated:

[0150]

[0151] Where: σ AR is the standard deviation of AR; μ AR is the average value of AR. The closer ACV is to 0, the higher the phase consistency of the two time series.

[0152] According to the range characteristics of the coefficient of variation of the amplitude ratio, it is reverse normalized to reduce the impact of outliers on the overall data distribution, ensure that it can effectively reflect the true differences between the original coefficients of variation, and eliminate the dimensional differences between different indicators to facilitate comparison between different evaluation indicators. The reverse normalization formula is as follows:

[0153]

[0154] The spatial distribution of the consistency results of the time series trend, phase and amplitude of the groundwater storage anomalies monitored by the combined JPL Mascon and GLDAS-Noah hydrological models and the water level changes observed in groundwater wells is shown in the figure below. Figure 3As shown. Among them, Figure 3 (a) represents the correlation coefficient CC between groundwater storage anomaly and groundwater level change, Figure 3 (b) represents the phase difference variation coefficient PCV, Figure 3 (c) represents the coefficient of variation of amplitude ratio ACV.

[0155] S2. The evaluation sample set is composed of the mutual correlation coefficient, phase difference variation coefficient and amplitude ratio variation coefficient of all pixels in the evaluation area. The weights of the mutual correlation coefficient, phase difference variation coefficient and amplitude ratio variation coefficient are determined by principal component analysis, and the comprehensive applicability index of the pixel-scale groundwater storage anomaly inversion product is calculated by weighted calculation.

[0156] The step S2 uses the principal component analysis method to integrate the consistency indicators of each time series and establishes a suitability evaluation method for groundwater storage anomaly inversion products. The calculation formula is as follows:

[0157]

[0158] Where: PCV p is the pth principal component; α p is the contribution rate of the corresponding principal component; z is the number of principal components whose cumulative contribution rate reaches 85% and above.

[0159] The applicability evaluation method of the constructed groundwater storage anomaly inversion product is normalized, and the calculation formula is as follows:

[0160]

[0161] Where: GSEI norm is the normalized applicability index value; GSEI, GSEI max and GSEI min They represent the original value, the maximum value and the minimum value of the improved groundwater storage anomaly inversion product applicability index in all pixels.

[0162] S3. Based on the comprehensive applicability index results of groundwater storage anomaly inversion products, the applicability levels are divided into different levels to evaluate the applicability of regional groundwater storage anomaly inversion products for analyzing the spatiotemporal changes of groundwater storage.

[0163] The classification processing method in step S3 is as follows:

[0164] The applicability index is graded and divided into five levels with 0.2 as an interval: the first level is [0, 0.2], indicating that the applicability of the groundwater storage anomaly inversion product is "poor"; the second level is (0.2, 0.4], indicating that the applicability is "poor"; the third level is (0.4, 0.6], indicating that the applicability is "general"; the fourth level is (0.6, 0.8], indicating that the applicability is "good"; the fifth level is (0.8, 1], indicating that the applicability is "excellent". Based on the above classification results, the applicability of the groundwater storage anomaly inversion product is evaluated and analyzed. The spatial distribution of the applicability index results of the groundwater storage anomaly monitored by JPL Mascon combined with the GLDAS-Noah hydrological model is shown in the figure below. Figure 3 (d) shown.

[0165] Figure 4 The following are time series diagrams of groundwater reserve anomalies and measured groundwater level changes for pixels 3-11 and 6-9. Using the applicability evaluation method from previous studies, the correlation coefficient between the groundwater reserve anomaly and the measured groundwater level change at pixel 3-11 was calculated to be 0.45, indicating general applicability. However, using the groundwater reserve anomaly applicability evaluation method proposed in this invention, the applicability index for pixel 3-11 was 0.16, indicating poor applicability. It can be seen that there is a significant difference between the two results. It can be intuitively seen from the figure that the phase consistency and amplitude consistency of the time series of the groundwater reserve anomaly and the measured groundwater level change at pixel 3-11 are low. The phase difference coefficient of variation (PCV) of the two time series for this pixel is 0.32, and the amplitude ratio coefficient of variation (ACV) is 0.39, which is basically consistent with the situation presented in the time series diagram, fully demonstrating the superiority of the evaluation method proposed in this invention.

[0166] The correlation coefficient between groundwater storage anomalies and measured groundwater level changes at pixels 6-9 is 0.68; using the proposed groundwater storage anomaly applicability evaluation method, the applicability index is 0.78. The figure shows a high degree of consistency in trend, phase, and amplitude between the two time series. The calculated coefficient of variation (PCV) for the phase difference between the two time series for this pixel is 0.67, and the coefficient of variation (ACV) for the amplitude ratio is 0.68, which is generally consistent with the time series plot.

Claims

1. A method for evaluating the applicability of groundwater reserve anomaly inversion products, characterized in that: The following steps are involved: Step S1: Generate regional groundwater storage anomaly inversion products based on the GRACE / GRACE-FO total water storage anomaly products and the surface water storage and soil water storage anomaly products simulated by the hydrological model. Calculate the cross-correlation coefficient, phase difference coefficient of variation, and amplitude ratio coefficient of variation between the pixel-scale measured water level time series and the groundwater storage anomaly time series. Evaluate the applicability of the groundwater storage anomaly inversion products from multiple perspectives: trend, phase, and amplitude. Step S2: The evaluation sample set is composed of the mutual correlation coefficient, phase difference variation coefficient, and amplitude ratio variation coefficient of all pixels in the evaluation area. The weights of the mutual correlation coefficient, phase difference variation coefficient, and amplitude ratio variation coefficient are determined using the principal component analysis method, and the comprehensive applicability index of the pixel-scale groundwater storage anomaly inversion product is weighted and calculated. Step S3: Based on the comprehensive applicability index results of the groundwater storage anomaly inversion product, the applicability levels are graded and evaluated to evaluate the applicability of the regional groundwater storage anomaly inversion product for the analysis of spatiotemporal changes in groundwater storage; The step S1 includes the following sub-steps: Sub-step S11: generating a regional groundwater storage anomaly inversion product based on the GRACE / GRACE-FO total water storage anomaly product and the surface water storage and soil water storage anomaly products simulated by the GLDAS hydrological model; Sub-step S12: obtaining regional measured groundwater level data; Sub-step S13: evaluating the applicability of regional groundwater storage anomaly inversion products based on measured groundwater levels from multiple perspectives; The multi-view suitability evaluation indexes of the sub-step S13 include the cross-correlation coefficient, the phase difference variation coefficient and the amplitude ratio variation coefficient; The cross-correlation coefficient calculation formula is: Where: σ x , σ y are the mean square errors of the time series x(t) and y(t), respectively; are the means of x(t) and y(t), respectively; k is the lag, which indicates the degree of offset between the two time series; C xy (k) and r xy (k) are the cross-covariance and cross-correlation coefficient of time series x(t) and y(t) at time lag k, respectively; n is the length of the time series; the cross-correlation value range is -1 to 1. The closer the cross-correlation coefficient of two time series is to 1, the more similar the changing trends of the two time series are. The calculation formula of the phase difference variation coefficient is: Where: σ PCE is the standard deviation of PCE; μ PCE is the average value of PCE; the closer PCV is to 0, the higher the phase consistency of the two time series; The calculation formula of the amplitude ratio variation coefficient is: Where: σ AR is the standard deviation of AR; μ AR is the average value of AR; the closer ACV is to 0, the higher the amplitude consistency of the two time series.

2. The method for evaluating the applicability of groundwater reserve anomaly inversion products according to claim 1, characterized in that: The expression of the groundwater reserve anomaly in sub-step S11 is: GWSA=TWSA-SWSA-SMSA Where GWSA represents groundwater storage anomaly, TWSA represents total water storage anomaly, SWSA represents surface water storage anomaly, and SMSA represents soil water storage anomaly.

3. The method for evaluating the applicability of groundwater reserve anomaly inversion products according to claim 1, characterized in that: The step S2 uses the principal component analysis method to determine the weights of the mutual correlation coefficient, the phase difference variation coefficient, and the amplitude ratio variation coefficient, and weightedly calculates the comprehensive applicability index of the pixel-scale groundwater storage anomaly inversion product. The calculation formula of the comprehensive applicability index of the groundwater storage anomaly inversion product is as follows: Where: PCV p is the pth principal component; α p is the contribution rate of the corresponding principal component; The comprehensive applicability index of the groundwater storage anomaly inversion product is normalized and the calculation formula is: Where: GSEI norm is the normalized comprehensive applicability index value; GSEI, GSEI max and GSEI min They represent the original value, the maximum value and the minimum value of the comprehensive applicability index of the groundwater storage anomaly inversion product in all pixels, respectively.

4. The method for evaluating the applicability of groundwater reserve anomaly inversion products according to claim 1, characterized in that: The grading and classification of the applicability levels in step S3 is specifically based on the comprehensive applicability index result of the groundwater storage anomaly inversion product, and the applicability index is graded and divided into five levels with 0.2 as an interval: the first-level numerical interval is [0, 0.2], indicating that the applicability of the groundwater storage anomaly inversion product is "poor"; the second-level numerical interval is (0.2, 0.4], indicating that the applicability is "poor"; the third-level numerical interval is (0.4, 0.6], indicating that the applicability is "general"; the fourth-level numerical interval is (0.6, 0.8], indicating that the applicability is "good"; and the fifth-level numerical interval is (0.8, 1], indicating that the applicability is "excellent". According to the above grading results, the applicability of the groundwater storage anomaly inversion product is comprehensively evaluated.

Images