Method for diagnosing dynamic transformation of land water reserves based on abrupt change time point analysis

Abstract

The invention discloses a land water reserve dynamic change diagnosis method based on sudden change time point analysis. The method comprises the following steps: S1, obtaining TWS abnormal data of a preset time length in a research area; s2, carrying out preprocessing on the TWS abnormal data of each grid unit; s3, according to the preprocessed TWS abnormal data of each grid unit, calculating a main mutation occurrence time point of each grid unit; s4, for each grid unit, calculating change trends of TWS abnormal data in time periods before and after a main sudden change occurrence time point; s5, classifying each grid unit into a preset dynamic transformation category; s6, acquiring regional hydrological climate index data corresponding to the time and space of the TWS abnormal data in the step S1; s7, analyzing the change trend of each hydrological climate index data in time periods before and after the main sudden change occurrence time point; and S8, carrying out correlation analysis on the dynamic transformation category of each grid unit in the step S5 and the change trend of each hydrological climate index data in the step S7, and outputting correlation information between the dynamic transformation category and the change trend.

Description

Technical Field

[0001] This invention relates to the field of changes in terrestrial water storage, and more specifically, to a diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis. Background Technology

[0002] Terrestrial water storage (TWS), as a comprehensive indicator of surface water, soil moisture, and groundwater, is a key element reflecting regional water resource status (Scanlon et al., 2018). Understanding TWS dynamics is crucial for assessing the response of hydrological systems to climate variability, supporting ecosystem health research, and providing foundational data for climate feedback models (Famiglietti, 2014; Rodell et al., 2018; Calvin et al., 2023).

[0003] Historically, comprehensive assessments of large-scale TWS (Telematics and Weather Anomalies) have been limited by fragmented point measurements (Milly et al., 2005; Trenberth, 2011). The Gravity Recovery and Climate Experiment (GRACE) and its follow-up mission (GRACE-FO) have significantly enhanced monitoring capabilities in this field by providing large-scale, continuous observations of TWS anomalies (Tapley et al., 2004). GRACE data have been widely used for understanding hydrological processes, monitoring drought, and assessing water resources (Rodell et al., 2018; Famiglietti et al., 2011; Thomas et al., 2014).

[0004] However, the response of hydrological systems to climate forcing is often nonlinear, and significant mechanistic shifts or abrupt changes may occur at specific points in time (Milly et al., 2008; Scheffer et al., 2012). Traditional TWS analysis methods, such as the widely used linear trend analysis or methods that only focus on static correlations, often assume that the system is stationary. This assumption has limitations when faced with nonlinear and abrupt behavior, making it difficult to accurately capture key mechanistic shifts occurring in TWS dynamics, and also failing to provide a deep understanding of the timing, type, and potential dynamic and non-stationary associations between these shifts and driving factors (Vishwakarma et al., 2021). Although some studies have used abrupt change detection methods to analyze the time series of hydrological elements such as runoff and precipitation, there is a lack of a systematic and universal method to systematically identify and characterize significant mechanistic shifts in large-scale terrestrial water storage dynamics using satellite remote sensing data (such as GRACE), quantify the changes in the timing, type, and stability of these shifts, and diagnose which changes in driving factors these shifts may be associated with. Recent studies have revealed increasing non-stationarity and variability in TWS dynamics across multiple locations globally (Zhong et al., 2025), as well as trend shifts during the GRACE / GRACE-FO record period (Buzzanga et al., 2025), highlighting the importance of developing new methods to diagnose shifts in terrestrial water storage mechanisms. Ignoring these potential mechanistic shifts could lead to misjudgments of water storage status and biased predictions of future water resource scenarios.

[0005] Furthermore, there is a close relationship between TWS dynamics and climate drivers, including regional hydroclimatic conditions (e.g., SPEI) and large-scale climate modes (e.g., ENSO and PDO). However, most existing studies focus on analyzing the static correlation between these factors and TWS, and our understanding of whether these associations change when TWS mechanisms shift—that is, the non-stationarity of the climate-TWS coupling relationship—remains insufficient. The difficulty in distinguishing the impact of internal variability and climate forcing on TWS dynamic transitions further complicates our understanding of complex hydrological processes (Buzzanga et al., 2025).

[0006] North America, with its diverse climate zones, complex topography, and known sensitivity to climate change—manifested as warming, altered precipitation patterns, and an increase in extreme hydroclimate events (Seager et al., 2007; Diffenbaugh et al., 2017; Cook et al., 2015)—is a typical region exhibiting the nonstationarity of global terrestrial water storage and its complex driving mechanisms. Therefore, using North America as a research background and application example can effectively demonstrate the technical value and capability of a method capable of diagnosing the dynamic changes in terrestrial water storage mechanisms and their dynamic correlation with driving factors.

[0007] The present invention aims to provide a universally applicable method that can systematically diagnose significant mechanistic shifts in terrestrial water storage dynamics, identify key time points of these shifts, characterize the types of shifts and their impact on stability, and correlate them with changes in hydro-climate drivers. This overcomes the limitations of existing methods, provides a deeper scientific understanding, and offers more accurate information for water resource management and risk assessment. Summary of the Invention

[0008] This invention provides a diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis, in order to solve the technical problems existing in the prior art.

[0009] To achieve the above objectives, this invention provides a method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis, comprising:

[0010] S1: Divide the study area into multiple grid cells and obtain TWS anomaly data of a preset time length within the study area;

[0011] S2: Preprocess the TWS anomaly data for each grid cell;

[0012] S3: Based on the TWS anomaly data of each preprocessed grid cell, calculate the main mutation occurrence time point of each grid cell. This main mutation occurrence time point is the critical moment when TWS dynamics change significantly.

[0013] S4: For each grid cell, calculate the trend of TWS anomalous data changes before and after the time point of the main mutation;

[0014] S5: Based on the changing trend of TWS, classify each grid cell into a preset dynamic transformation category;

[0015] S6: Obtain the regional hydrological and climate index data corresponding to the TWS abnormal data in step S1 in terms of time and space. The regional hydrological and climate index data characterizes the regional hydrological and climate conditions.

[0016] S7: Analyze the changing trends of various hydro-climate index data before and after the major abrupt changes;

[0017] S8: Perform correlation analysis between the dynamic transformation category of each grid cell in step S5 and the changing trend of each hydrological and climate index data in step S7, and output the correlation information between the two.

[0018] In one embodiment of the present invention, step S2 is as follows:

[0019] S21: In each grid cell, the TWS outlier data is decomposed into seasonal, trend, and residual terms using the least squares fitting method:

[0020] Y(t) = S(t) + T(t) + R(t)

[0021] Where Y(t) is the original TWS outlier data, S(t) is the seasonal term, T(t) is the trend term, and R(t) is the residual term.

[0022] S22: Discard the seasonal term S(t) and treat T(t)+R(t) as TWS abnormal data for subsequent processing.

[0023] In one embodiment of the present invention, in step S3, the major mutation occurrence time point of each grid cell is calculated using the Pettitt method. Step S3 is as follows:

[0024] S31: Set a test statistic U for each grid cell. t t = 1, 2, ..., N, where N is the total length of TWS abnormal data.

[0025]

[0026] Where sgn(x) is the sign function; if x > 0, sgn(x) equals +1; if x < 0, sgn(x) equals -1; if x = 0, sgn(x) equals 0. i and X j These are the i-th and j-th observations in the TWS anomaly data, respectively;

[0027] S32: Calculate U t The time t when the maximum value is reached CPT ;

[0028] S33: Use the significance p-value for time t determined by S32. CPT An evaluation was conducted; if p < 0.05, it indicates that t CPT This was identified as the time point where the major mutation occurred, at which point the TWS dynamics changed significantly.

[0029] In one embodiment of the present invention, step S4 includes:

[0030] S41: The Theil-Sen slope estimation method was used to calculate the TWS change trend before and after the time point of the major mutation.

[0031] β=median((X j -X i ) / (ji)), for all 1≤i <j≤n seg ,

[0032] Where β is the trend slope, X i and X j These are the i-th and j-th observations in the TWS anomaly data, respectively, and n seg The length of the time period;

[0033] S42: Use the Mann-Kendall test to assess the significance of the trend slope β. If p < 0.05, it indicates that the trend of TWS abnormal data during this period is significant.

[0034] In one embodiment of the present invention, step S5 includes: classifying each grid cell into a preset dynamic transformation category based on the direction and significance of the TWS abnormal data change trend before and after the major mutation occurrence time point. The dynamic transformation category includes a negative transformation category and a positive transformation category. The negative transformation category indicates that the trend of TWS abnormal data changes to a more negative direction or changes from positive to negative. The positive transformation category indicates that the trend of TWS abnormal data changes to a more positive direction or changes from negative to positive.

[0035] In one embodiment of the present invention, in step S5, the classification of the preset dynamic transition categories is based on the combination of the significance and direction of the TWS change trend before and after the time point of the major mutation. The preset dynamic transition categories include at least one category that represents an aggravated decline or a change to a decline, and at least one category that represents an accelerated rise or a slowed decline.

[0036] In one embodiment of the present invention, the method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis further includes step S9:

[0037] Using the time point of the main mutation as the time limit, calculate and output the change of the relative stability index of TWS abnormal data before and after the time point of the main mutation, and include it in the output of step S8.

[0038] In one embodiment of the present invention, in step S9, the relative stability index is the relative difference standard deviation (SDRD), and step S9 includes:

[0039] S91: Calculate the numerical values ​​of TWS abnormal data at each time point during this period.t Compared with the average TWS abnormal data during this period, TWS seg_mean Relative deviation RD t ,

[0040] RD t =(TWS) t -TWS seg_mean ) / TWS seg_mean

[0041] S92: Calculate the SDRD for this period as the standard deviation of these relative deviations:

[0042]

[0043] RD mean It is the average relative deviation for that period, n seg is the length of the time period, and n is the number of time points in the time period.

[0044] In one embodiment of the present invention, the method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis further includes step S10: using cosine similarity analysis to determine the directional consistency between the spatial gradient of the main abrupt change time points and the surface elevation field.

[0045] S101: Calculate the spatial gradient vector in each grid cell. and gradient vector gradient vector The gradient vector represents the direction of the steepest increase in space for the years when the major mutations occur. The direction from which the elevation increases most steeply from the DEM;

[0046] S102: Gradient vector The angle is denoted as θ CPT , gradient vector The angle is denoted as θ Z， θ CPT θ Z Convert to unit vector v respectively CPT and v Z :

[0047] v CPT =[cos(θ) CPT ),sin(θ CPT )]

[0048] v Z =[cos(θ) Z ),sin(θ Z )]

[0049] S103: Calculate the cosine similarity of the mutation point locations as a dot product to obtain the calculation result CosineSimilarity:

[0050] CosineSimilarity = v CPT ·v Z

[0051] A CosineSimilarity value close to 1 indicates that the gradient vector... gradient vector The two directions are parallel, and -1 indicates the gradient vector. gradient vector The two directions are antiparallel, and a value close to 0 indicates the gradient vector. gradient vector The two are orthogonal;

[0052] S104: Calculate the regional average GlobalCorr:

[0053]

[0054] GlobalCorr quantizes the gradient vector using the regional average. gradient vector The overall directional trend of both, where M is the number of grid cells in the region and k is the index of the grid cell.

[0055] In one embodiment of the present invention, the hydroclimatic index data is the standardized precipitation evapotranspiration index (SPEI), and step S7 includes:

[0056] S71: Define drought events as consecutive months with SPEI < -0.5 and wet events as consecutive months with SPEI > 0.5. Calculate the duration D1 of drought events and the duration D2 of wet events, and determine the severity S1 of drought events and the severity S2 of wet events.

[0057]

[0058] Where i represents the i-th drought event and j represents the j-th wet event;

[0059] S72: Calculate the intensity I1 of drought events and the intensity I2 of wet events:

[0060] I1 = S1 / D1

[0061] I² = S² / D²

[0062] S73: Calculate the frequency, average duration, average severity, and average intensity of drought and wet events before and after each major mutation point.

[0063] The diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis provided by this invention offers an integrated and repeatable diagnostic tool. It can generate quantitative information on the spatiotemporal distribution, type classification, stability changes, and dynamic correlations with driving factors of TWS mechanism changes. This provides decision support for water resource management, risk assessment, and climate adaptation that surpasses traditional static analysis. In summary, this invention has the following beneficial technical effects:

[0064] (1) It provides a systematic and objective method for quantifying the transformation of TWS non-stationary mechanisms.

[0065] (2) By using dynamic comparison of anchored major mutation points, the diagnostic accuracy of TWS state transition characteristics and time is improved.

[0066] (3) It reveals the dynamic connection between the dynamic transformation of TWS and the changes in the characteristics of driving factors (especially extreme events), and deepens the understanding of climate-hydrology interaction.

[0067] (4) It can generate regional diagnostic information that distinguishes different TWS evolution paths, providing a scientific basis for formulating differentiated and adaptive water resource management strategies.

[0068] (5) The method is reproducible and easy to apply to different regions and datasets. Attached Figure Description

[0069] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, 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.

[0070] Figure 1 A topographic map of the North American continent;

[0071] Figure 2 This refers to the dataset involved in the embodiments of the present invention;

[0072] Figures 3a-3d Spatiotemporal characteristics of North American terrestrial water storage (TWS) trends (2003-2022);

[0073] Figure 4a , Figure 4b Spatial distribution and spatial autocorrelation analysis of the time of change points (CPT) in terrestrial water storage (TWS);

[0074] Figures 5a-5d Spatial patterns and statistical analysis of temporal stability of terrestrial water storage (TWS) before and after identified abrupt change points;

[0075] Figures 6a-6g To investigate the coupling relationship between regional hydroclimate transition and the dynamics of terrestrial water storage anomalies (TWSA) before and after the identified abrupt change point (CPT);

[0076] Figures 7a-7g The temporal evolution of the correlation between the regional mean land water storage anomaly (TWSA) and large-scale climate indices (2003–2022). Detailed Implementation

[0077] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0078] This invention focuses on the North American continent (NAC), defined as the region between 130° and 60° west longitude and 25° and 75° north latitude. Figure 1 The figure shows the boundaries of the study area (brown outline), major drainage systems (blue lines), major mountain ranges (Cordillera, Pacific Coast Range, and Appalachian Range), and the extent of permafrost (cross-shaded areas; data source: Ranet al., 2022). Arrows indicate the dominant atmospheric circulation patterns influencing moisture transport in the region: cold, dry Arctic airflow (blue); warm, moist airflow from the Pacific and Atlantic / Gulf of Mexico (orange); and the North American monsoon (orange dashed line). Figure 1 The topography of the continent is shown (elevation, in meters). The NAC was chosen as the region for studying the hydrological response to climate change because of its climate sensitivity, manifested in recent decades' warming, altered precipitation patterns, and an increase in extreme hydroclimate events. The continent's terrestrial water storage (TWS) exhibits spatiotemporal heterogeneity and complex dynamics. Topographically, the NAC has significant elevation differences, with mountain ranges in the west (Cordillera) and east (Appalachian Mountains) surrounding the central plain. Figure 1 This large-scale topography influences atmospheric circulation and water vapor transport, thus affecting regional TWS patterns. Permafrost exists in high-latitude and high-altitude regions, influencing hydrological processes, particularly in the northern basins. The continent also contains major freshwater systems, including the Great Lakes and large river basins such as the Mississippi and St. Lawrence Rivers. The combination of diverse climate zones, complex topography, significant water bodies, and recorded climate sensitivity makes North America a suitable region for studying TWS dynamics, abrupt changes, and their climate drivers.

[0079] This invention relates to three main datasets: satellite gravity observation data for TWS, drought index for hydro-climatic conditions, and climate teleconnection index, summarized in... Figure 2 .

[0080] Figure 2 The terrestrial water storage shown is based on satellite gravity observation data used for TWS, as well as monthly TWS anomaly data, derived from the GRACE and GRACE-FO RL06 Mascon solutions from the Center for Space Research (CSR) at the University of Texas at Austin.

[0081] Figure 2 The drought index shown: Regional hydroclimatic conditions are characterized using the Standardized Precipitation Evapotranspiration Index (SPEI) from the SPEIbase v2.9 dataset, calculated on a one-month timescale. The SPEI integrates the effects of precipitation and potential evapotranspiration, providing a climate water balance measure related to TWS fluctuations.

[0082] Figure 2 The climate teleconnection indices shown: To assess the impacts of large-scale climate variability, monthly indices of El Niño-Southern Oscillation and Pacific Decadal Oscillation (PDO) were obtained from NOAA Physical Sciences Laboratory. These indices represent known atmospheric and oceanic circulation patterns that influence North American hydroclimate.

[0083] This invention provides a diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis, comprising:

[0084] S1: Divide the study area into multiple grid cells and obtain TWS anomaly data of a preset time length within the study area;

[0085] S2: Preprocess the TWS anomaly data for each grid cell;

[0086] S3: Based on the TWS anomaly data of each preprocessed grid cell, calculate the main mutation occurrence time point of each grid cell. This main mutation occurrence time point is the critical moment when TWS dynamics change significantly.

[0087] S4: For each grid cell, calculate the trend of TWS anomalous data changes before and after the time point of the main mutation;

[0088] S5: Based on the changing trend of TWS, classify each grid cell into a preset dynamic transformation category;

[0089] S6: Obtain the regional hydrological and climate index data corresponding to the TWS abnormal data in step S1 in terms of time and space. The regional hydrological and climate index data characterizes the regional hydrological and climate conditions.

[0090] S7: Analyze the changing trends of various hydro-climate index data before and after the major abrupt changes;

[0091] S8: Perform correlation analysis between the dynamic transformation category of each grid cell in step S5 and the changing trend of each hydrological and climate index data in step S7, and output the correlation information between the two.

[0092] In one embodiment of the present invention, step S2 is as follows:

[0093] S21: In each grid cell, the TWS outlier data is decomposed into seasonal, trend, and residual terms using the least squares fitting method:

[0094] Y(t) = S(t) + T(t) + R(t)

[0095] Where Y(t) is the original TWS outlier data, S(t) is the seasonal term, T(t) is the trend term, and R(t) is the residual term.

[0096] S22: Discard the seasonal term S(t) and treat T(t)+R(t) as TWS abnormal data for subsequent processing.

[0097] It should be noted that a small number of data gaps in the TWS anomalous data (approximately 11 months during the GRACE / GRACE-FO transition period) can be filled using linear interpolation. Given the short duration of the gaps relative to the 20-year study period and the robustness of the nonparametric method used, the impact of interpolation on detecting major long-term mutation points and trends is considered minimal.

[0098] In one embodiment of the present invention, in step S3, the major mutation occurrence time point of each grid cell is calculated by Pettitt test, and step S3 is as follows:

[0099] S31: Set a test statistic U for each grid cell. t t = 1, 2, ..., N, where N is the total length of TWS abnormal data.

[0100]

[0101] Where sgn(x) is the sign function; if x > 0, sgn(x) equals +1; if x < 0, sgn(x) equals -1; if x = 0, sgn(x) equals 0. i and X j These are the i-th and j-th observations in the TWS anomaly data, respectively;

[0102] S32: Calculate U t The time t when the maximum value is reached CPT ;

[0103] S33: Use the significance p-value for time t determined by S32. CPT An evaluation was conducted; if p < 0.05, it indicates that t CPT This was identified as the time point where the major mutation occurred, at which point the TWS dynamics changed significantly.

[0104] Step S3 above applies the nonparametric Pettitt test independently to the deseasoned TWS time series for each grid cell to identify the single most significant transition in the mean. This rank-based test requires no assumptions about the data distribution and is robust to outliers, making it applicable to potentially nonnormal hydrological data. Given the approximately 20-year record length of Grace, identifying a single dominant CPT provides a statistically reliable characterization, avoiding the reduced statistical power that can occur when detecting multiple breakpoints in shorter sequences.

[0105] To characterize the major mutation point (CPT), the corresponding time is denoted as t. CPT The change in the TWS cumulative rate after t can be quantified using the nonparametric Theil-Sen slope estimation before CPT (time from 1 to t). CPT ) and CPT after (time from t) CPT The linear trend from +1 to N). In one embodiment of the invention, step S4 includes:

[0106] S41: The Theil-Sen slope estimation method was used to calculate the TWS change trend before and after the time point of the major mutation.

[0107] β=median((X j -X i ) / (ji)), for all 1≤i <j≤n seg ,

[0108] Where β is the trend slope, X i and X j These are the i-th and j-th observations in the TWS anomaly data, respectively, and n seg The length of the time period;

[0109] S42: Use the Mann-Kendall test to assess the significance of the trend slope β. If p < 0.05, it indicates that the trend of TWS abnormal data during this period is significant.

[0110] In one embodiment of the present invention, step S5 includes: classifying each grid cell into a preset dynamic transformation category based on the direction and significance of the TWS abnormal data change trend before and after the major mutation occurrence time point. The dynamic transformation category includes a negative transformation category and a positive transformation category. The negative transformation category represents the trend of TWS abnormal data changing in a more negative direction or from positive to negative, such as an accelerated decline before the major mutation occurrence time point, an increase before the change to a decline, or an increase before the change to a slowdown. The positive transformation category represents the trend of TWS abnormal data changing in a more positive direction or from negative to positive, such as an increase before the major mutation occurrence time point, an increase before the change to an increase, or an increase before the change to a slowdown.

[0111] In one embodiment of the present invention, in step S5, the classification of the preset dynamic transition categories is based on the combination of the significance and direction of the TWS change trend before and after the time point of the major mutation. The preset dynamic transition categories include at least one category that represents an aggravated decline or a change to a decline, and at least one category that represents an accelerated rise or a slowed decline.

[0112] In one embodiment of the present invention, the method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis further includes step S9:

[0113] Using the time point of the main mutation as the time limit, calculate and output the change of the relative stability index of TWS abnormal data before and after the time point of the main mutation, and include it in the output of step S8.

[0114] In this invention, the relative standard deviation of difference (SDRD) can be used to assess the change in TWS time variability before and after CPT. SDRD quantifies the variability relative to the average TWS level over a period of time. In one embodiment of this invention, in step S9, the relative stability index is the relative standard deviation of difference (SDRD), and step S9 includes:

[0115] S91: Calculate the numerical values ​​of TWS abnormal data at each time point during this period. t Compared with the average TWS abnormal data during this period, TWS seg_mean Relative deviation RD t ,

[0116] RD t =(TWS) t -TWS seg_mean ) / TWS seg_mean

[0117] S92: Calculate the SDRD for this period as the standard deviation of these relative deviations:

[0118]

[0119] RD mean It is the average relative deviation for that period, n seg is the length of the time period, and n is the number of time points in the time period.

[0120] To investigate whether the spatial evolution of CPT is related to topography, this invention can also use cosine similarity analysis to determine the directional consistency between the spatial gradient of the CPT time node field and the surface elevation field. Cosine similarity is a metric sensitive to directional consistency. In one embodiment of this invention, the method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis further includes step S10: using cosine similarity analysis to determine the directional consistency between the spatial gradient of the main abrupt change time node field and the surface elevation field.

[0121] S101: Calculate the spatial gradient vector in each grid cell. and gradient vector gradient vector The gradient vector represents the direction of the steepest increase in space for the years when the major mutations occur. The direction from which the elevation increases most steeply from the DEM;

[0122] S102: Gradient vector The angle is denoted as θ CPT , gradient vector The angle is denoted as θ Z， θ CPT θ Z Convert to unit vector v respectively CPT and v Z :

[0123] v CPT =[cos(θ) CPT ),sin(θ CPT )]

[0124] v Z =[cos(θ) Z ),sin(θ Z )]

[0125] S103: Calculate the cosine similarity of the mutation point locations as a dot product to obtain the calculation result CosineSimilarity:

[0126] CosineSimilarity = v CPT ·v Z

[0127] A CosineSimilarity value close to 1 indicates that the gradient vector... gradient vector The two directions are parallel, and -1 indicates the gradient vector. gradient vector The two directions are antiparallel, and a value close to 0 indicates the gradient vector. gradient vector The two are orthogonal;

[0128] S104: Calculate the regional average GlobalCorr:

[0129]

[0130] GlobalCorr quantizes the gradient vector using the regional average. gradient vector The overall directional trend of both, where M is the number of grid cells in the region and k is the index of the grid cell.

[0131] This invention can also use runs theory applied to regional SPEI time series to study the relationship between TWS dynamics and regional hydro-climate conditions. Here, it applies... Figure 2 The drought index data in the data. In one embodiment of the present invention, the hydroclimate index data is the standardized precipitation evapotranspiration index (SPEI), and step S7 includes:

[0132] S71: Define drought events as consecutive months with SPEI < -0.5 and wet events as consecutive months with SPEI > 0.5. Calculate the duration D1 of drought events and the duration D2 of wet events, and determine the severity S1 of drought events and the severity S2 of wet events.

[0133]

[0134] Where i represents the i-th drought event and j represents the j-th wet event;

[0135] S72: Calculate the intensity I1 of drought events and the intensity I2 of wet events:

[0136] I1 = S1 / D1

[0137] I² = S² / D²

[0138] S73: Calculate the frequency, average duration, average severity, and average intensity of drought and wet events before and after each major mutation point.

[0139] Based on the above, this invention further explores potential climate impacts from two perspectives. First, it qualitatively and quantitatively compares the shifts in TWS characteristics (trend type, temporal stability) with the corresponding shifts in drought and wet event characteristics derived from regional SPEI, to assess the coupling between TWS pattern shifts and regional hydro-climate change. Second, it uses Spearman rank correlation to study large-scale climate modes (… Figure 2 The impact of the climate teleconnection index (ENSO, PDO) is shown. This nonparametric method assesses the monotonic association between monthly regional TWS anomalies and monthly climate indices. To account for potential lagged responses, monthly lagged correlations (TWS lags behind the climate indices) from 0 to 12 months were calculated. The lags that produce the highest statistical significance (p < 0.05) correlations were identified as the regional characteristic response times.

[0140] Output results of this invention embodiment:

[0141] The diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis described in this invention can be applied to different geographical regions. The following will use North America as an example to demonstrate the application process of this method in a specific embodiment and the obtained diagnostic and correlation analysis results.

[0142] I. Overall Spatiotemporal Pattern of Changes in North American Land Water Storage

[0143] Figures 3a-3d Spatiotemporal characteristics of North American terrestrial water storage (TWS) trends (2003-2022), Figure 3a The figure shows the overall linear trend (cm / year) calculated throughout the study period. Figure 3b The figure shows the trend (cm / year) before the identified mutation point (CPT). Figure 3c The trend (cm / year) is shown after the major mutation point. Figures 3a-3c The histogram in the image shows the frequency distribution of the trend values ​​of the grid cells. Figure 3d Spatial distribution of trend changes: Lower developmental types (LD, LD1-LD3) indicate a negative shift in the TWS accumulation rate after the major inflection point (CPT), corresponding to... Figure 3d The T1 to T3 range, while the higher development type (HD, HD1-HD3 type) represents a positive transition, corresponding to... Figure 3d The pie chart shows the percentage distribution of LD and HD classifications across the entire study area, from T4 to T6. Figure 3b and Figure 3c The statistical significance (p<0.05) was assessed using the Mann-Kendall test.

[0144] Figure 4a , Figure 4bSpatial distribution and spatial autocorrelation analysis of the time of change points (CPT) in terrestrial water storage (TWS). Figure 4a The statistically significant CPT years identified in each grid cell (Pettitt test, p<0.05) are shown in the radar chart, which displays the frequency distribution of CPT occurrences by year. Figure 4b The local spatial autocorrelation index (LISA) clustering plot (p<0.05) identifies statistically significant hotspots (late CPT) and cold spots (early CPT) at CPT time. The plot also shows the statistics and significance level of the Global Moran's I index.

[0145] Analysis of GRACE / GRACE-FO data (2002–2022) reveals spatial heterogeneity in North American terrestrial water storage (TWS) trends. Figure 3a Differentiation exists near the 98°W meridian: net TWS decreases (0.3–1.7 cm / year) in the west coast of the United States, the central plains of Canada, and the southwestern United States, while net TWS increases (0.5–3.0 cm / year) in the central plains of the United States, the northeastern United States / Great Lakes region, eastern Canada, and parts of the southeastern United States.

[0146] The 20-year trend masks the temporal shift. The Pettitt test identified statistically significant (p<0.05) abrupt change points (CPTs) in the TWS time series of most of the study area (>99% of the grid cells). This non-stationarity indicates the need for piecewise trend analysis before and after the CPTs.

[0147] Comparison before CPTs identification ( Figure 3b ) and after ( Figure 3c The TWS trend shows widespread changes. Approximately 54% of the study area experienced an accelerated decline or a reversal to decline in TWS after CPT (negative development, LD type); Figure 3d This pattern is evident in the western regions; for example, the rate of TWS decline near the Sierra Nevada and Rocky Mountains increased by 1.5 to 2.5 cm / year after CPT. Conversely, the remaining approximately 46% of the region showed a trend reversal with either an increase in TWS or a slowdown in decline (positive development, HD type). Figure 3d In the Great Lakes region and the East Coast, TWS increased by 0.3 to 1.2 cm / year after CPT. Spatially, LD regions are mainly distributed in the Central Plains, west of the Rocky Mountains, and Mexico, while HD regions are concentrated in the East Coast of the United States, around the Great Lakes, and eastern Canada. Figure 3d ).

[0148] The time points of TWS mutations also exhibit spatial organization. Figure 4aAlthough CPTs occurred throughout the period from 2002 to 2022, a concentration was observed between 2010 and 2014. Figure 4a Local spatial autocorrelation index (LISA) analysis identified statistically significant hotspots (late-stage CPTs) and coldspots (early-stage CPTs) (p<0.05). Figure 4b Global spatial autocorrelation analysis yielded Moran's I of 0.86 (Z-score = 470.70, p < 0.001). Figure 4b The data indicates that the distribution of CPT years is clustered, suggesting a large-scale synergistic effect in the transformation of the TWS pattern.

[0149] II. Temporal stability of terrestrial water storage changes and topographic effects

[0150] Figures 5a-5d To analyze the spatial pattern and statistical significance of temporal stability of terrestrial water storage (TWS) before and after identified abrupt change points, temporal stability was quantified using the standard deviation of relative differences (SDRD). Figure 5a SDRD calculated for the period prior to the mutation point (SDRD) 01 ). Figure 5b SDRD calculated for the period following the mutation point (SDRD) 02 ). Figure 5c The ratio of SDRD after the mutation point to SDRD before the mutation point (SDRD) 02 / SDRD 01 A value >1 indicates decreased stability (increased relative variability) after CPT, while a value <1 indicates increased stability. Figure 5d Display North American SDRD 01 (Red), SDRD 02 The probability density function of the distribution of (blue) and SDRD ratio (gray) values.

[0151] An assessment of the time stability (SDRD) of TWS revealed changes in North America after CPT. While the spatial pattern of relative variability persisted ( Figure 5a , Figure 5b However, the magnitude of SDRD increased after CPT in regions including the south-central plains of the United States (extending to the Gulf Coast), Mexico, and east-central Canada.

[0152] SDRD ratio (SDRD 02 / SDRD 01 ; Figure 5cThis stability shift was quantified. Decreasing stability (ratio > 1) is characteristic of most of the continent, including central and eastern Canada, central and eastern United States, and Mexico. Increasing stability (ratio < 1) is less common and spatially scattered, mainly found in parts of the western Cordillera and isolated Arctic patches. Probability density functions (PDFs) quantify this stability shift. Figure 5d This is supported by: SDRD after CPT 02 Distribution relative to SDRD before CPT 01 The distribution shifts to the right, and the PDF of the SDRD ratio is right-skewed (mode around 1.2–1.3). This indicates that the decrease in TWS temporal stability is the main result across the continent after CPT.

[0153] III. Climate Drivers of Changes in Terrestrial Water Storage

[0154] Figures 6a-6g Time series of the relationship between the regionally standardized precipitation evapotranspiration index (SPEI) and the regional average TWSA (black line) for seven subdivided regions. SPEI values ​​are represented by shaded areas, with dark red / blue indicating significant drought (SPEI < -0.5) / wet (SPEI > +0.5) events, and light / transparent shading indicating milder conditions (-0.5 ≤ SPEI ≤ +0.5). The vertical dashed line represents the area-weighted average CPT year. (h) Stacked bar charts show the relative changes (%) in the frequency, duration, and intensity of drought and wet events after and before CPT. Bars above / below zero indicate increases / decreases in the corresponding characteristic, and different colors distinguish the six event indicators.

[0155] Figures 7a-7g The time series for seven divided regions is shown, displaying: regional average TWSA (black line, left axis, cm); representing the monthly ENSO status. 3.4 Index values ​​(blue / red lines, right axis); and monthly Pacific Decadal Oscillation (PDO) index values ​​(purple lines, right axis). Light red and light blue shading represent officially defined El Niño and La Niña events, respectively.

[0156] Analysis of the characteristics of dry and wet events as defined by SPEI before and after regional abrupt change points (CPTs) shows that shifts in hydroclimate patterns are coupled with observed terrestrial water storage (TWS) trajectories. TWSSA fluctuations are consistent with dry and wet events as defined by SPEI. Figures 6a-6g A decrease in TWS typically corresponds to a dry period, while replenishment corresponds to a wet period. Based on the hydroclimatic evolution following CPT, the regions were classified as follows:

[0157] The drought-intensified regions (B01, B02, B05) tend towards even drier conditions after the CPT, which is consistent with their negative development (LD) TWS trajectory. Figure 6a (6b, 6e). These regions experienced an increase in drought frequency (e.g., 30% increase in region B01) or intensity (e.g., 72% increase in region B05) after the CPT, while wetting events generally diminished (e.g., a 74% decrease in the frequency of wetting events in region B01). Conversely, regions with enhanced wetting (B03, B06) shifted to wetter conditions or experienced a reduction in drought severity after the CPT, consistent with the positive development (HD) TWS trend. Figure 6c This includes an increase in the duration (78%) and intensity (25%) of wet events in region B03, and a decrease in the frequency of droughts (-52%) in region B06. The complex change regions (B04, B07) exhibit a complex hydroclimatic pattern, with both wet and dry extreme events intensifying, which may correspond to greater TWS variability after the CPT. Figure 6d (6g). Area B04 experienced an increase in the frequency of drought (64%) and the frequency of wet events (219%); while area B07 saw an increase in both the frequency of drought (70%) and wet events (approximately 90%).

[0158] To investigate the impact of ENSO and PDO on regional TWS dynamics, the TWSA responses of seven regions were analyzed. Spearman rank correlation coefficients and optimal lag times are shown below. Figures 7a-7g As shown. In most regions, TWS showed a statistically significant correlation with ENSO and PDO (coefficients typically between 0.1 and 0.45), with responses either synchronous or lagging (up to approximately one year). Figures 7a-7g ).

[0159] The characteristics of TWS changes differed during ENSO and PDO events. In southern North America / Mexico (B01), TWS showed a synchronous negative correlation with both ENSO and PDO; during the 2010–2012 La Niña period, TWS increased by 15.5 cm, while during the 2015–2016 El Niño period, it decreased by 10.2 cm. Figure 7a In contrast, TWS in the southwestern United States (B02) was positively correlated with ENSO and PDO (lagging by several months); a TWS deficit of 13.9 cm was observed during the 2020-2023 La Niña period. Figure 7b The Missouri River Basin (B04) showed a lagged positive correlation; the 2009 El Niño corresponded to an increase of 11.3 cm in TWS (with a lag of about one year), while the 2020-2023 La Niña coincided with a decrease of 16.1 cm in TWS (with a lag). Figure 7d ).

[0160] In the Pacific Coast Range (B05), TWS is negatively correlated with ENSO (synchronous) and PDO (lagging by approximately 10 months), accompanied by a long-term downward trend. La Niña events (e.g., 9.3 cm during 2010–2012) coincide with a temporary slowdown in TWS decline, while El Niño events exacerbate the depletion. During 2014–2018 (including the effects of lagged El Niño and positive PDO), the region lost 45.8 cm of TWS. Figure 7e Although TWS in Alaska (B07) is synchronously positively correlated with ENSO / PDO, its response is complex. During the 2015–2017 El Niño and positive PDO period, TWS decreased by 4.9 cm; while during the 2020–2023 La Niña and negative PDO period, it increased by 2.5 cm. Figure 7g These discrepancies suggest that other factors or more complex response mechanisms dominate the TWS dynamics in this region.

[0161] The diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis provided by this invention offers an integrated and repeatable diagnostic tool. It can generate quantitative information on the spatiotemporal distribution, type classification, stability changes, and dynamic correlations with driving factors of TWS mechanism changes. This provides decision support for water resource management, risk assessment, and climate adaptation that surpasses traditional static analysis. In summary, this invention has the following beneficial technical effects:

[0162] (1) It provides a systematic and objective method for quantifying the transformation of TWS non-stationary mechanisms.

[0163] (2) By using dynamic comparison of anchored major mutation points, the diagnostic accuracy of TWS state transition characteristics and time is improved.

[0164] (3) It reveals the dynamic connection between the dynamic transformation of TWS and the changes in the characteristics of driving factors (especially extreme events), and deepens the understanding of climate-hydrology interaction.

[0165] (4) It can generate regional diagnostic information that distinguishes different TWS evolution paths, providing a scientific basis for formulating differentiated and adaptive water resource management strategies.

[0166] (5) The method is reproducible and easy to apply to different regions and datasets.

[0167] Those skilled in the art will understand that the modules in the apparatus of the embodiments can be distributed in the apparatus of the embodiments as described in the embodiments, or they can be located in one or more devices different from this embodiment with corresponding changes. The modules of the above embodiments can be combined into one module, or they can be further divided into multiple sub-modules.

Claims

1. A diagnostic method for dynamic changes in terrestrial water storage based on abrupt change time point analysis, characterized in that, include: S1: Divide the study area into multiple grid cells and obtain TWS anomaly data of a preset time length within the study area; S2: Preprocess the TWS anomaly data for each grid cell; S3: Based on the TWS anomaly data of each preprocessed grid cell, calculate the main mutation occurrence time point of each grid cell. This main mutation occurrence time point is the critical moment when TWS dynamics change significantly. S4: For each grid cell, calculate the trend of TWS anomalous data changes before and after the time point of the main mutation; S5: Based on the changing trend of TWS, classify each grid cell into a preset dynamic transformation category; S6: Obtain the regional hydrological and climate index data corresponding to the TWS abnormal data in step S1 in terms of time and space. The regional hydrological and climate index data characterizes the regional hydrological and climate conditions. S7: Analyze the changing trends of various hydro-climate index data before and after the major abrupt changes; S8: Perform correlation analysis between the dynamic transformation category of each grid cell in step S5 and the changing trend of each hydrological and climate index data in step S7, and output the correlation information between the two.

2. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, Step S2 is as follows: S21: In each grid cell, the TWS outlier data is decomposed into seasonal, trend, and residual terms using the least squares fitting method: Y(t) = S(t) + T(t) + R(t) Where Y(t) is the original TWS outlier data, S(t) is the seasonal term, T(t) is the trend term, and R(t) is the residual term. S22: Discard the seasonal term S(t) and treat T(t)+R(t) as TWS abnormal data for subsequent processing.

3. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, In step S3, the major mutation time points of each grid cell are calculated using the Pettitt method. Step S3 is as follows: S31: Set a test statistic U for each grid cell. t t = 1, 2, ..., N, where N is the total length of TWS abnormal data. Where sgn(x) is the sign function; if x > 0, sgn(x) equals +1; if x < 0, sgn(x) equals -1; if x = 0, sgn(x) equals 0. i and X j These are the i-th and j-th observations in the TWS anomaly data, respectively; S32: Calculate U t The time t when the maximum value is reached CPT ; S33: Use the significance p-value for time t determined by S32. CPT An evaluation was conducted; if p < 0.05, it indicates that t CPT This was identified as the time point where the major mutation occurred, at which point the TWS dynamics changed significantly.

4. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, Step S4 includes: S41: The Theil-Sen slope estimation method was used to calculate the TWS change trend before and after the time point of the major mutation. β=median((X j -X i ) / (ji)), for all 1≤i <j≤n seg , Where β is the trend slope, X i and X j These are the i-th and j-th observations in the TWS anomaly data, respectively, and n seg The length of the time period; S42: Use the Mann-Kendall test to assess the significance of the trend slope β. If p < 0.05, it indicates that the trend of TWS abnormal data during this period is significant.

5. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 4, characterized in that, Step S5 includes: classifying each grid cell into a preset dynamic transition category based on the direction and significance of the TWS abnormal data change trend before and after the main mutation occurrence time point. The dynamic transition category includes a negative transition category and a positive transition category. The negative transition category indicates that the trend of TWS abnormal data changes to a more negative direction or changes from positive to negative. The positive transition category indicates that the trend of TWS abnormal data changes to a more positive direction or changes from negative to positive.

6. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, In step S5, the classification of the preset dynamic transition categories is based on the combination of the significance and direction of the TWS change trend before and after the time point of the main mutation. The preset dynamic transition categories include at least one category that represents an accelerated decline or a change to a decline, and at least one category that represents an accelerated rise or a slowed decline.

7. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, It also includes step S9: Using the time point of the main mutation as the time limit, calculate and output the change of the relative stability index of TWS abnormal data before and after the time point of the main mutation, and include it in the output of step S8.

8. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 7, characterized in that, In step S9, the relative stability index is the relative standard deviation of difference (SDRD). Step S9 includes: S91: Calculate the numerical values ​​of TWS abnormal data at each time point during this period. t Compared with the average TWS abnormal data during this period, TWS seg_mean Relative deviation RD t , RD t =(TWS t -TWS seg_mean ) / TWS seg_mean S92: Calculate the SDRD for this period as the standard deviation of these relative deviations: RD mean It is the average relative deviation for that period, n seg is the length of the time period, and n is the number of time points in the time period.

9. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, It also includes step S10: using cosine similarity analysis to determine the directional consistency between the spatial gradient at the time points of major mutations and the surface elevation field. S101: Calculate the spatial gradient vector in each grid cell. and gradient vector gradient vector The gradient vector represents the direction of the steepest increase in space for the years when the major mutations occur. The direction from which the elevation increases most steeply from the DEM; S102: Gradient vector The angle is denoted as θ CPT , gradient vector The angle is denoted as θ Z， θ CPT θ Z Convert to unit vector v respectively CPT and v Z : v CPT =[cos(θ CPT ),sin(θ CPT )] v Z =[cos(θ Z ),sin(θ Z )] S103: Calculate the cosine similarity of the mutation point locations as a dot product to obtain the calculation result CosineSimilarity: CosineSimilarity=v CPT ·v Z A cosineSimilarity value close to 1 indicates that the gradient vector... gradient vector The two directions are parallel, and -1 indicates the gradient vector. gradient vector The two directions are antiparallel, and a value close to 0 indicates the gradient vector. gradient vector The two are orthogonal; S104: Calculate the regional average GlobalCorr: GlobalCorr quantizes the gradient vector using the regional average. gradient vector The overall directional trend of both, where M is the number of grid cells in the region and k is the index of the grid cell.

10. The method for diagnosing dynamic changes in terrestrial water storage based on abrupt change time point analysis according to claim 1, characterized in that, The hydroclimatic index data is the Standardized Precipitation Evapotranspiration Index (SPEI). Step S7 includes: S71: Define drought events as consecutive months with SPEI < -0.5 and wet events as consecutive months with SPEI > 0.

5. Calculate the duration D1 of drought events and the duration D2 of wet events, and determine the severity S1 of drought events and the severity S2 of wet events. Where i represents the i-th drought event and j represents the j-th wet event; S72: Calculate the intensity I1 of drought events and the intensity I2 of wet events: I1 = S1 / D1 I² = S² / D² S73: Calculate the frequency, average duration, average severity, and average intensity of drought and wet events before and after each major mutation point.

Images