Underground water reserve abrupt change point detection method combining signal decomposition and sequential test method

Abstract

The invention discloses an underground water reserve abrupt change point detection method combining signal decomposition and a sequential test method, which comprises the following steps: firstly, acquiring a gravity satellite spherical harmonic coefficient product and GLDAS hydrological mode data of a to-be-detected area, and performing inversion to obtain underground water reserve abnormal data of the to-be-detected area; missing data interpolation and signal decomposition are carried out on groundwater reserve abnormal data obtained through inversion, and a non-linear trend term, a seasonal period term and a signal decomposition residual term which do not contain data missing are obtained; and finally, performing underground water reserve abrupt change point identification and detection by using the nonlinear trend term and the signal residual term of the signal decomposition. According to the method, residual terms of signal decomposition are utilized to construct the break variable test threshold value, robust mutation point recognition and detection are carried out, in the actual application process, the break variable test threshold value well eliminates false alarm points of mutation point recognition, and the accuracy and robustness of mutation point detection are guaranteed.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of gravity satellite groundwater storage inversion, and particularly relates to a groundwater storage mutation point detection method combining signal decomposition and sequential testing method. BACKGROUND

[0002] Groundwater is one of the most important water resources in China. Evaluating groundwater storage (GWS) and exploring its temporal variation law have important practical and scientific significance for water resources management, climate change monitoring, and food security protection. Current gravity satellite programs, mainly GRACE (Gravity Recovery and Climate Experiment) and GRACE-FO (Gravity Recovery and Climate Experiment Follow On), are the only technical means to obtain large-area GWS time series. Domestic and foreign scholars have conducted research on the variation law of GWS in the North China Plain based on the gravity field products of GRACE and GRACE-FO, and have made a series of beneficial results in the overall change rate. However, the groundwater system is a complex unity, and the complex coupling relationship with various factors leads to significant segmental mutation characteristics of the linear change rate. Unfortunately, current researches mostly focus on the overall change rate in the entire research period, and there is no sufficient consideration of the change rate estimation results considering the segmental mutation characteristics of the linear change. In view of the important value of the linear change rate of GWS to regional water resources management and social sustainable development, it is of important theoretical research and practical application significance to study the mutation point detection method of GWS time series to improve the estimation accuracy of the linear change rate. SUMMARY

[0003] The present application aims to at least partially solve the technical problems in the related art.

[0004] The present application aims to at least partially solve the technical problems in the related art.

[0005] In order to achieve the above-mentioned purpose, the present application provides a groundwater storage mutation point detection method combining signal decomposition and sequential testing method, comprising the following steps:

[0006] S1, obtaining gravity satellite spherical harmonic coefficient products and GLDAS hydrological model data of a region to be detected, and inversely obtaining groundwater storage anomaly data of the region to be detected ;

[0007] S2, performing missing data interpolation and signal decomposition on the inversely obtained groundwater storage anomaly data to obtain a nonlinear trend item without data missing Seasonal periodic terms and signal decomposition residues ;

[0008] S3. Utilizing the nonlinear trend term of signal decomposition and signal residuals Identification and detection of abrupt changes in groundwater storage.

[0009] As a preferred option, step S1 specifically involves:

[0010] S11. Land water storage data of the area to be tested is obtained by inversion based on gravity satellite spherical harmonic coefficient products acquired through multiple channels. and its outliers ;

[0011] S12. Obtain the vegetation canopy water storage in the area to be tested using GLDAS NOAH hydrological model data. Snowmelt water reserves and soil water storage And calculate its corresponding outlier. , and ;

[0012] S13. Based on the above data, the groundwater storage anomaly in the area to be detected is calculated. Time series.

[0013] Preferably, step S11 is as follows:

[0014] S111. Based on spherical harmonic coefficient products obtained from multiple channels, the terrestrial water storage of the area to be tested is calculated using the mean method. :

[0015]

[0016] in, , and These represent the data obtained from the inversion of three data products: J, G, and C. data;

[0017] S112, Based on the above The data, after being processed to remove the mean, was used to retrieve outliers in terrestrial water storage. data:

[0018]

[0019] in, This indicates the terrestrial water storage in the area to be tested during the research period. average value.

[0020] Preferably, the vegetation canopy water storage in step S12 Snowmelt water reserves and soil water storage The outlier calculation method is as follows:

[0021]

[0022] in, It is the average water storage of the vegetation canopy during the study period. It is the average snowmelt water storage during the study period. , These represent the soil water storage at different depths and the average soil water storage at different depths during the study period, respectively. Indicates different depths.

[0023] Preferably, the groundwater storage anomaly value mentioned in step S13 The time series calculation method is as follows:

[0024] .

[0025] As a preferred option, step S2 specifically involves:

[0026] S21. Using singular spectrum analysis methods to... Data interpolation is performed on the missing points;

[0027] S22. Use STL methods to process the interpolated data. Time series signals are decomposed to obtain nonlinear trend terms. Seasonal periodic terms and residual terms .

[0028] As a preferred option, step S3 specifically involves:

[0029] S31. Using signal decomposition of residual terms Construct an adaptive mutation detection threshold :

[0030]

[0031] in, As a scale factor, The length of the residual term, The first term in the residual term A number;

[0032] S32, Based on groundwater storage anomalies nonlinear trend term The set of abrupt changes in groundwater storage in the area to be tested was calculated. ,in ;

[0033] S33, based on the groundwater reserves mutation point set and nonlinear trend term , calculate the mutation point corresponding to the mutation ;

[0034] S34, in turn, compare the mutation point of the mutation and adaptive mutation variable test threshold size relationship:

[0035] When the mutation variable less than the threshold , eliminate the mutation point;

[0036] When the mutation variable greater than the threshold , then keep the mutation point;

[0037] S35, based on the updated mutation point set of step S34 , step S33-S34, until all the mutation points through the detection.

[0038] As a preferred, step S32 is specifically:

[0039] S321, set the subsequence length to , set the restrictive level of T test to , calculate the difference value of the adjacent two subsequence :

[0040]

[0041] Wherein, T represents the T test value with degree of freedom , significance level , the average value of the standard deviation of each subsequence ; S322, calculate the mean of the previous subsequence of time point t

[0042] and the mean of the next subsequence ,

[0043] ;

[0044] S323, based on the mean of the previous subsequence of time point t and the mean of the next subsequence and the difference value of the adjacent two subsequence , calculate the critical mean of state change and sign index ​

[0045] ;

[0046] S324、According to the critical mean value , the state mutation index of the current moment is calculated

[0047]

[0048] When the state mutation index and the symbol index are the same, the first time point is the mutation point of the trend change; the time of the identified mutation point is defined as , wherein .

[0049] As preferred, step S33 is specifically:

[0050] According to the identified mutation point , the entire nonlinear trend term is divided into different subsequences, one-dimensional linear regression is performed on each subsequence, and then the linear change speed of each subsequence is obtained:

[0051]

[0052] wherein, and are the intercepts of the first and the first subsequences, and are the linear change speeds of the first and the first subsequences, is the kth time point in the entire time sequence;

[0053] When the linear change speeds of the two adjacent subsequences before and after the first mutation point do not change in direction, i.e. , the mutation value of the mutation point is the difference between the average values of the two adjacent subsequences; when the linear change speeds of the two adjacent subsequences before and after the first mutation point change in direction, i.e. , the mutation value of the mutation point is infinite.

[0054] Beneficial effects: firstly, the GWS time series of the North China Plain is obtained based on the GRACE, GRACE-FO gravity field spherical harmonic coefficient products and GLDAS NOAG data inversion; then, the signal decomposition is performed on the GWS time series by the singular spectrum analysis method and the STL method; finally, the mutation point identification and detection are performed by the STARS method and the signal residual term. Considering the complexity of the mutation point detection and the internal correlation between different mutation points, the residual term of the signal decomposition is used to construct the mutation variable test threshold in the present application, so that the robust mutation point identification and detection are realized. In the actual application process, the mutation variable test threshold can well eliminate the false alarm points in the mutation point identification, and the accuracy and robustness of the mutation point detection are ensured. BRIEF DESCRIPTION OF DRAWINGS

[0055] Figure 1 is a flow chart of the groundwater storage mutation point detection method of the present application combined with signal decomposition and sequential testing method.

[0056] Figure 2 is the example verification result of the present application; Figure 2 Fig. 2 (a) is the GWS time series of the North China Plain from April 2002 to May 2024 obtained by the present application; Figure 2 Fig. 2 (b) is the entire nonlinear trend item, mutation point and linear change rate of each subsequence of the groundwater storage change in the North China Plain from April 2002 to May 2024 in the embodiment of the present application. DETAILED DESCRIPTION

[0057] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions in the present application will be described clearly and completely below in combination with the drawings in the present application. Obviously, the described embodiments are some embodiments of the present application, but not all the embodiments. They should not be understood as limiting the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application. In the description of the present application, it should be understood that the terms used are only for the purpose of description, and should not be understood as indicating or implying relative importance.

[0058] The present application provides a groundwater storage mutation point detection method combined with signal decomposition and sequential testing method. Figures 1-2 The present application provides a groundwater storage mutation point detection method combined with signal decomposition and sequential testing method.

[0059] Embodiment 1: as shown in the figure, the present embodiment provides a groundwater storage mutation point detection method combined with signal decomposition and sequential testing method, which comprises the following steps: Figure 1

[0060] ​S1, obtain the gravity satellite spherical harmonic coefficient product and GLDAS hydrological model data of the North China Plain region, and obtain the groundwater storage anomaly data of the North China Plain region by inversion ;

[0061] S11, based on the GRACE RL06 spherical harmonic coefficient products provided by JPL, GFZ and CSR, obtain the land water storage data of the North China Plain region of different agencies by inversion 、 and , and convert it into the overall land water storage data by using the mean value operator

[0062] (1)

[0063] wherein, 、 and respectively represent the data obtained by inversion of JPL, GFZ and CSR data products.

[0064] Based on the above data, the formula is used for mean value removal processing to invert the land water storage change

[0065] (2)

[0066] wherein, represents the average value in the research period.

[0067] S12, according to the GLDAS NOAH hydrological model data, extract the canopy water , snowmelt and soil water of the North China Plain region in turn, and use formula for mean value removal processing

[0068] (3)

[0069] wherein, is the average value of the vegetation canopy water storage in the research period, is the average value of the snowmelt water storage in the research period, 、 respectively represent the soil water storage of different depths and the average value of the soil water storage of different depths in the research period, represents different depths.

[0070] ​​S13, based on the water balance equation, the above data are processed by difference to obtain the groundwater storage anomaly data of the North China Plain time series

[0071] (4)

[0072] S2, the missing data interpolation and signal decomposition are performed on the inversion obtained groundwater storage anomaly data to obtain the nonlinear trend item , seasonal cycle item and signal decomposition residual item without data missing;

[0073] S21, the singular spectrum analysis method is used to interpolate the missing points of the time series;

[0074] S22, the STL method is used to perform signal decomposition on the interpolated time series to obtain the nonlinear trend item , seasonal cycle item and residual item .

[0075] S3, the nonlinear trend item and signal residual item of signal decomposition are used for groundwater storage anomaly point identification and detection.

[0076] S31, the signal decomposition residual item is used to construct an adaptive mutation variable test threshold ,

[0077] (5)

[0078] wherein, is a scale factor, is the length of the residual item, is the th value in the residual item.

[0079] S32, the length of the subsequence is set to , the restrictive level of T test is set to , and the deviation value of the adjacent two subsequences is calculated :

[0080] (6)

[0081] wherein, represents the T test value with the degree of freedom and the significance level , and is the average value of the standard deviation of each subsequence .

[0082] For the tth time point, calculate the mean of the previous subsequence of the time point and the mean of the next subsequence ,

[0083] (7)

[0084] On this basis, calculate the critical mean of state change and the sign index

[0085] (8)

[0086] According to the critical mean , calculate the state mutation index of the current time

[0087] (9)

[0088] When the state mutation index and the sign index are the same, the tth time point is the mutation point of trend change; if mutation points are identified, the time of the mutation point is defined as , wherein .

[0089] S33, according to the identified mutation point , divide the entire nonlinear trend item into different subsequences, perform one-dimensional linear regression on each subsequence, and then obtain the linear change speed of each subsequence.

[0090] The fitting method of the two subsequences before and after the tth mutation point (i.e., the tth and the t+1th) is given as follows,

[0091] (10)

[0092] wherein, and are the intercepts of the tth and the t+1th subsequences, and are the linear change speeds of the tth and the t+1th subsequences. ​​​​​​​​​is the kth time point in the whole time series;

[0093] When the linear change speed of the two adjacent sub-sequences before and after the kth mutation point does not change direction , the mutation amount of the mutation point is the difference between the average values of the two adjacent sub-sequences; when the linear change speed of the two adjacent sub-sequences before and after the kth mutation point changes direction, i.e. , the mutation amount of the mutation point is infinite.

[0094] S34, sequentially compare the mutation amount of each mutation point with the size relationship of the adaptive mutation amount detection threshold value: When the mutation amount is less than the threshold value

[0095] , the mutation point is removed; When the mutation amount is greater than the threshold value

[0096] , the mutation point is retained.

[0097] S35, using the updated mutation point set , cycle S33 and S34 until all mutation points pass the detection.

[0098] In order to verify the reliability of the method for detecting mutation points proposed in the present application, the present application uses the spherical harmonic coefficient gravity field products of North China Plain provided by CSR, GFZ and JPL from April 2002 to May 2024, the hydrological model products provided by GLDAS, as shown in (a) of FIG. 8, to obtain the time series of groundwater storage GWS in North China Plain by inversion, and sequentially perform data interpolation, signal decomposition, mutation point identification and mutation point detection by using the method of the present application, wherein the scale factor is 1.2, and the detection result is shown in FIG. 9. Figure 2 Figure 2

[0099] ​​​​​​​​​The signal decomposition result shows that the standard deviation of the decomposition residual term is 2.42 cm, and the mutation variable threshold value constructed therefrom is 2.91 cm; the initial mutation points identified by the first mutation point identification are 8, respectively, 2004-11-01, 2007-03-01, 2009-09-01, 2011-11-01, 2015-11-01, 2016-10-01, 2018-03-01, 2020-03-01, wherein the mutation variables of the 2nd and 5th mutation points are 1.13 cm and 2.68 cm, respectively, which are less than the mutation variable threshold value, so they are excluded; after the second detection, the linear change speeds before and after the remaining 6 mutation points change in direction, so the remaining 6 mutation points are the final mutation points, which are 2004-11-01, 2009-09-01, 2011-11-01, 2016-10-01, 2018-03-01, and 2020-03-01. Figure 2 As shown in (b), from the time distribution result of the GWS nonlinear trend term, the detected mutation points are located at the critical points of the speed change, which embodies the effectiveness and reliability of the method.

[0100] In summary, the groundwater storage mutation point detection method combining signal decomposition and sequential testing method has good effectiveness and reliability, can effectively carry out mutation point detection of groundwater storage time series, and can be used as a feasible method for characteristic analysis of groundwater storage.

[0101] Finally, it should be noted that: the above examples are only used to illustrate the technical solutions of the present application, but not to limit them; although the present application has been described in detail with reference to the foregoing examples, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing examples, or make equivalent replacement for part of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for detecting abrupt changes in groundwater storage using a combined signal decomposition and sequential testing approach, characterized in that, Includes the following steps: S1. Obtain gravity satellite spherical harmonic coefficient products and GLDAS hydrological model data for the area to be monitored, and then retrieve the groundwater storage anomaly data for the area to be monitored. ; S2. Anomaly data of groundwater storage obtained from the inversion. Perform missing data interpolation and signal decomposition to obtain a nonlinear trend term without missing data. Seasonal periodic terms and signal decomposition residues ; S3. Utilizing the nonlinear trend term of signal decomposition and signal residuals Identification and detection of abrupt changes in groundwater storage.

2. The groundwater storage mutation point detection method according to claim 1, characterized in that, Step S1 is as follows: S11. Land water storage data of the area to be tested is obtained by inversion based on gravity satellite spherical harmonic coefficient products acquired through multiple channels. and its outliers ; S12. Obtain the vegetation canopy water storage in the area to be tested using GLDAS NOAH hydrological model data. Snowmelt water reserves and soil water storage And calculate its corresponding outlier. , and ; S13. Based on the above data, the groundwater storage anomaly in the area to be detected is calculated. Time series.

3. The groundwater storage mutation point detection method according to claim 2, characterized in that, Step S11 is as follows: S111. Based on spherical harmonic coefficient products obtained from multiple channels, the terrestrial water storage of the area to be tested is calculated using the mean method. : ； in, , and These represent the data obtained from the inversion of three data products: J, G, and C. data; S112, Based on the above The data, after being processed to remove the mean, was used to retrieve outliers in terrestrial water storage. data: ； in, This indicates the terrestrial water storage in the area to be tested during the research period. average value.

4. The groundwater storage mutation point detection method according to claim 3, characterized in that, The vegetation canopy water storage mentioned in step S12 Snowmelt water reserves and soil water storage The outlier calculation method is as follows: ； in, It is the average water storage of the vegetation canopy during the study period. It is the average snowmelt water storage during the study period. , These represent the soil water storage at different depths and the average soil water storage at different depths during the study period, respectively. Indicates different depths.

5. The groundwater storage mutation point detection method according to claim 4, characterized in that, The groundwater storage anomaly mentioned in step S13 The time series calculation method is as follows: 。 6. The groundwater storage mutation point detection method according to claim 5, characterized in that, Step S2 is as follows: S21. Using singular spectrum analysis methods to... Data interpolation is performed on the missing points; S22. Use STL methods to process the interpolated data. Time series signals are decomposed to obtain nonlinear trend terms. Seasonal periodic terms and residual terms .

7. The groundwater storage mutation point detection method according to claim 6, characterized in that, Step S3 is as follows: S31. Using signal decomposition of residual terms Construct an adaptive mutation detection threshold : ； in, As a scale factor, The length of the residual term. For the residual term A number; S32, Based on groundwater storage anomalies nonlinear trend term The set of abrupt changes in groundwater storage in the area to be tested was calculated. ,in ; S33, Based on the set of groundwater storage mutation points and nonlinear trend term Calculate the mutation amount corresponding to the mutation point. ; S34. Compare the mutation amounts at each mutation point sequentially. With adaptive mutation amount test threshold Size relationship: When mutation amount Less than the threshold When the mutation point is removed, it should be discarded. When mutation amount Greater than the threshold In this case, the mutation point is preserved; S35. Based on the updated set of mutation points from step S34 Repeat steps S33-S34 until all mutation points pass the detection.

8. The groundwater storage mutation point detection method according to claim 7, characterized in that, Step S32 is as follows: S321. Set the subsequence length to... Set the restriction level for the T-test to 0. Calculate the difference between two adjacent subsequences. : ； in, Describing the degrees of freedom as The significance level is The T-test value, For each subsequence The mean of the standard deviations; S322. Calculate the mean of the subsequence preceding time point t. and the mean of the next subsequence , ； S323, The mean of the previous subsequence based on time point t and the mean of the next subsequence and the difference between two adjacent subsequences Calculate the critical mean of state change. and symbol index , ； S324. Based on the critical mean Calculate the state change index at the current moment. , ； When the state change index and symbol index When the symbols are the same, the first The time points are The abrupt change point of the trend; the identified The time of each mutation point is defined as... ,in .

9. The groundwater storage mutation point detection method according to claim 8, characterized in that, Step S33 is as follows: Based on the identified mutation points The entire nonlinear trend term Divided into For each of the three distinct subsequences, a one-dimensional linear regression is performed to determine the linear rate of change of each subsequence. ； in, and The first The and the first The intercept of each subsequence, and For the first The and the first The linear rate of change of each subsequence This represents the k-th time point in the entire time series. When the The linear change rate of two adjacent subsequences before and after the mutation point did not change direction, i.e. Then the mutation amount at that mutation point This is the difference between the averages of two adjacent subsequences; when the linear rate of change of the two adjacent subsequences changes direction, i.e. Then the mutation amount at that mutation point It is infinitely large.

Images