Groundwater spatial distribution detection method, device and medium based on timing constraints
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- ANHUI & HUAI RIVER WATER RESOURCES RES INST
- Filing Date
- 2026-04-03
- Publication Date
- 2026-06-12
Smart Images
Figure CN122198249A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of groundwater distribution detection technology, specifically a method, equipment, and medium for detecting the spatial distribution of groundwater based on time-series constraints. Background Technology
[0002] With the intensification of global climate change and human activities, accurate and dynamic monitoring and assessment of groundwater resources are crucial for achieving sustainable water resource management. Traditional technologies use global navigation satellite systems to observe vertical crustal deformation to invert changes in groundwater storage, and utilize the bi-basic geometric relationship between direct satellite signal and surface reflection to capture hydrologically relevant physical parameters. The observation data is one of the main signal sources for inverting groundwater quality migration. However, field measurements show that the deployment of traditional observation stations is not adapted to the spatial distribution characteristics of aquifers, and insufficient observation point density leads to uneven spatial coverage, making it difficult to characterize the details of groundwater spatial distribution at the watershed or aquifer scale. At the same time, the observation strategy uses a fixed elevation angle range and sampling frequency, and the time resolution is usually monthly, which cannot capture the deformation signals corresponding to the dynamic changes of groundwater within a short period.
[0003] To overcome the aforementioned spatial coverage and resolution issues, existing technologies have developed methods for multi-source remote sensing data fusion and joint inversion. By adjusting the spatial layout of observation stations to match the distribution of hydrogeological units, introducing a dynamic elevation angle range selection method to eliminate interference signals from non-hydrological factors, and strengthening the collaborative processing of satellite reflection signals and direct signals, theoretically, the extraction accuracy and spatial representativeness of groundwater-related deformation signals can be improved.
[0004] However, in practical use, it still has some shortcomings. For example, when integrating GPS time series data, the existing methods fail to optimize the sampling frequency by combining the hydrological periodic characteristics of groundwater dynamic changes, resulting in incomplete separation of multi-scale hydrological signals from non-hydrological interference signals. Due to the failure to optimize and quantify the uncertainty of GPS time series data for hydrological signal extraction, the optimization of the layout of observation points in areas with complex hydrogeological structures is not targeted, making it difficult to accurately invert the spatial heterogeneity characteristics of groundwater changes. Summary of the Invention
[0005] This invention provides a method, device, and medium for detecting the spatial distribution of groundwater based on time-series constraints. By acquiring vertical displacement data from observation stations, extracting groundwater deformation components through multi-scale decomposition, optimizing the network layout of observation stations, and finally constraining and solving the distribution field of groundwater storage changes, the invention addresses the problems in the background art.
[0006] To achieve the above objectives, the technical solution of the present invention is as follows: A time-constrained method for detecting the spatial distribution of groundwater involves the following steps performed using computer equipment: S1: Obtain vertical displacement data of observation stations in the target area to form first time series data, and preprocess the first time series data to generate first time series features aligned with hydrological features; S2: Perform multi-scale decomposition on the first time series features, and use non-hydrological deformation signals as constraints to identify and extract the pure hydrological load deformation components caused by groundwater changes, so as to obtain the second time series features; S3: Based on the extraction process of the second time series feature, the uncertainty of groundwater changes reflected by each observation station is quantitatively assessed and summarized as the overall uncertainty; S4: With the goal of minimizing overall uncertainty, the network layout of each observation station is simulated and optimized to generate the optimized observation network and its corresponding spatial constraint weights; S5: Using the second temporal feature and the spatial constraint weight as spatial constraint conditions, solve the groundwater storage change distribution field to generate a spatiotemporally continuous third temporal feature as the result of groundwater spatial distribution detection.
[0007] Preferably, the first temporal feature in S1 specifically includes: analyzing the prior hydrological periodic characteristics of groundwater dynamics at the observation stations in the target area; The prior hydrological periodic characteristics include at least interannual trends, seasonal cycles, and short-cycle fluctuations caused by artificial irrigation activities; wherein, during the active period of the short-cycle fluctuations, high-frequency deformation signals related to groundwater changes are captured. Based on the prior hydrological periodic characteristics, a resampling strategy and filtering parameters matching different hydrological activity periods are dynamically determined to generate a first time-series feature aligned with the hydrological features.
[0008] Preferably, in S2, non-hydrological deformation signals are used as constraints to identify and extract the second time series features, specifically including: using the non-hydrological deformation time series obtained from the hydrological physical model simulation as the identification constraint; The hydrophysical model includes at least a surface load model for simulating atmospheric and ocean load effects, and a terrestrial hydrological model for simulating changes in soil water and snowmelt. The non-hydrological deformation time series is correlated with each intrinsic mode component obtained by multi-scale decomposition and phase comparison is performed to identify and remove the mode components associated with the non-hydrological driving source; the remaining mode components are reconstructed to obtain the second time series feature, namely the pure hydrological load deformation component.
[0009] Preferably, the non-hydrological deformation time series is subjected to correlation analysis and phase comparison with each intrinsic mode component obtained through multi-scale decomposition, specifically including: S201: For the i-th intrinsic mode component obtained from the target temporal decomposition Calculate the Pearson correlation coefficient between it and the corresponding component of the non-hydrological deformation time series. ; S202: To Perform a Fourier transform on the model signal to calculate its amplitude squared coherence near the preset hydrological cycle frequency. The model signal is related to The k-th model intrinsic mode component with the closest time scale is denoted as In calculating amplitude squared coherence At the same time, calculate and Phase difference at hydrological periodic frequency points ; S203: Pre-set strong correlation threshold and weak correlation threshold ; S204: Perform identification and removal: like and At key frequencies > 0.8: At the same time, if | If | < 30°, then the two are considered to be in phase. If it is confirmed to be a non-hydrological mode, it will be removed; if | If the angle is greater than 150°, then the two are determined to be out of phase, and the operation is switched to AND. The same frequency domain adaptive filtering process is used to process the cases; like Then determine For mixed modes of hydrological and non-hydrological signals, frequency domain adaptive filtering is used for attenuation: like Then determine Hydrological signals, in particular, should be preserved in their entirety. S205: Perform the reconstruction: Directly superimpose all intrinsic mode components and trend residuals that were determined to be fully retained or retained in attenuated form as determined in S204 to generate the second time series feature.
[0010] Preferably, in S204, the determination is made. For mixed modes of hydrological and non-hydrological signals, frequency-domain adaptive filtering is used for attenuation, specifically including: right and model signals Perform Fourier transforms on each; In the frequency domain, for each discrete frequency point If the coherence between the two at that point Then for Complex amplitude at this frequency point Attenuation:
[0011] in, Represented as a preset attenuation intensity coefficient; for the attenuated spectrum Perform an inverse Fourier transform to obtain the filtered time-domain signal.
[0012] Preferably, in S3, a quantization matrix is constructed for each observation station gpsi. This is used to assess the uncertainty of groundwater changes reflected by each observation station, specifically expressed as: For a time series with t time points It is a t×t matrix; The elements on the diagonal of the matrix represent each time point. Total error variance It consists of the sum of the variances of three parts of the error:
[0013] in, This is represented as the inherent error between observation and solution. This represents the error introduced during the signal separation process. This is represented as hydrological characterization error; Off-diagonal elements This indicates a point in time. and The error covariance between the two is used to quantify the correlation strength of errors at different time points.
[0014] Preferably, in S4, the network layout of each observation station is simulated and optimized, specifically including: An optimization function is constructed with the objective of minimizing the posterior uncertainty of the inversion results of regional groundwater storage changes. Based on the hydrogeological unit boundaries, fault lines, and spatial distribution map of uncertainties of existing observation stations in the target area, a set of virtual observation candidate points is generated in key heterogeneous areas and uncertainty cavity areas. A sequential greedy algorithm is used to iteratively select virtual observation points from the set of virtual observation candidate points that can reduce the function value of the optimization function the most, and add them to the optimization network; The optimized observation network scheme and its corresponding spatial constraint weight map are output, wherein the spatial constraint weight map represents the quantitative difference in the optimized network's ability to constrain groundwater inversion at different locations within the region.
[0015] Preferably, in S5, the third temporal feature is a time-series, high spatial resolution groundwater storage change grid data, where each grid value represents the change in groundwater storage of the corresponding pixel at a specific time relative to the long-term average.
[0016] In another aspect, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the method described above.
[0017] In another aspect, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the method described above.
[0018] As can be seen from the above technical solution compared with the prior art, the present invention has the following beneficial effects: 1. This invention introduces a hydrological physical model to identify and eliminate non-hydrological interference modal components and reconstruct pure hydrological load deformation components, which significantly alleviates the problem of incomplete separation between multi-scale hydrological signals and non-hydrological interference signals, and realizes optimized analysis of time-series data for hydrological signal extraction. 2. This invention quantifies the uncertainty of time series analysis by constructing a quantification matrix for each observation station, thus overcoming the shortcoming of existing methods that lack uncertainty quantification and providing a precise basis for subsequent observation layout optimization. By constructing an optimization function, it solves the problem of the lack of specificity in the optimization of observation point layout in complex hydrogeological areas and strengthens the targeted coverage of the observation network for groundwater detection. Attached Figure Description
[0019] Figure 1 This is a schematic diagram of the method steps of the present invention; Figure 2 This is a schematic diagram of the timing characteristic change process of the present invention. Detailed Implementation
[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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 some embodiments of the present invention, but not all embodiments.
[0021] The embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are used to illustrate the present invention, but should not be used to limit the scope of the present invention.
[0022] Example: In this embodiment, as Figure 1 The temporal-constrained groundwater spatial distribution detection method shown herein acquires vertical displacement data from observation stations, extracts groundwater deformation components through multi-scale decomposition, optimizes the observation station network layout, and finally solves for the constrained distribution field of groundwater storage changes; specifically, it includes the following steps: S1: Obtain vertical displacement data of observation stations in the target area to form first time series data, and preprocess the first time series data to generate first time series features aligned with hydrological features; S2: Perform multi-scale decomposition on the first time series features, and use non-hydrological deformation signals as constraints to identify and extract the pure hydrological load deformation components caused by groundwater changes, so as to obtain the second time series features; S3: Based on the extraction process of the second time series feature, the uncertainty of groundwater changes reflected by each observation station is quantitatively assessed and summarized as the overall uncertainty; S4: With the goal of minimizing overall uncertainty, the network layout of each observation station is simulated and optimized to generate the optimized observation network and its corresponding spatial constraint weights; S5: Using the second temporal feature and the spatial constraint weight as spatial constraint conditions, solve the groundwater storage change distribution field to generate a spatiotemporally continuous third temporal feature as the result of groundwater spatial distribution detection.
[0023] Implementation begins: like Figure 2 As shown, specifically in S1, observation data and ephemeris from all available GPS continuous observation stations in and around the target area are obtained from authoritative data centers such as the International GNSS Service Organization and the China Crustal Movement Observation Network. The observation data includes at least carrier phase and pseudorange. Then, by using precise single-point positioning technology, daily coordinate time series are generated, and the vertical displacement component is extracted as the first time series data. After preprocessing, the first time series features include at least the displacement value, the effective hydrological activity period label corresponding to each data point, and the resampling rate and filtering parameters used for each data point.
[0024] It should be noted that precise point positioning technology can be achieved using scientific software such as GAMIT / GLOBK, Bernese, or GPSTk. When performing calculations, it is necessary to constrain satellite orbit and clock error products under the International Earth Reference Frame and consistently apply solid tide, ocean tidal load, polar tide correction models, and non-tidal atmospheric and ocean load models to obtain high-precision vertical displacement time series that primarily reflects changes in surface mass load.
[0025] In this embodiment, the selection of observation stations should follow these guidelines: priority should be given to stations with a continuous observation time span covering at least one complete hydrological year to ensure the extraction of reliable hydrological cycle characteristics; stations located on stable bedrock or that have undergone strict environmental error correction should be selected, while stations located in soft soil, landslides, or near large buildings should be avoided to minimize local non-tectonic and non-hydrological deformation interference; and the distribution of stations should cover the main hydrogeological units.
[0026] The generation of the first time-series feature aligned with hydrological characteristics includes: analyzing the prior hydrological cycle characteristics of groundwater dynamics in the target area, which include at least interannual trends, seasonal cycles, and short-cycle fluctuations caused by artificial irrigation activities; based on the prior hydrological cycle characteristics, dynamically determining resampling strategies and filtering parameters that match different hydrological activity periods to generate the first time-series feature aligned with hydrological characteristics; wherein, during the active period of short-cycle fluctuations, a higher sampling frequency is used to capture high-frequency deformation signals related to groundwater changes.
[0027] Before generating the first time-series feature, a hydrological activity period label needs to be defined for each data point on the time axis based on the prior hydrological cycle characteristics. The label generation rule can be based on a preset date range or dynamically triggered by the derivative of real-time hydrological data exceeding a specific threshold. In this embodiment, the regional soil water storage change rate is calculated based on the surface soil moisture content of the Global Land Surface Data Assimilation System (GLDASNoah) model. If the regional soil water storage change rate is positive for three consecutive days and is greater than the arithmetic mean of all data from the start of the current hydrological year to this date, then the period is marked as a rapid change period.
[0028] The resampling strategy and filtering parameters are designed to match the time scale of data processing with the inherent period of hydrological signals. In the preferred example, the 3-day average is designed for events such as irrigation that may cause rapid replenishment and discharge of groundwater within several days, aiming to preserve this high-frequency information. The 15-day moving median filter length is approximately equal to the duration of a typical weather-scale event, which can effectively suppress such non-hydrological noise. The monthly average and the 6-month moving average are suitable for characterizing slow seasonal and interannual variations. Those skilled in the art can adaptively adjust the filter window length to be within one to two times the known duration of hydrological events in the target area; a higher sampling frequency specifically means that the sampling interval should be less than half of the shortest period of the target hydrological signal to ensure that the signal is not lost.
[0029] It should be noted that the interannual trend characteristics are determined by studying long-term precipitation, temperature data and historical hydrological reports of the target area to determine whether the groundwater system is in a long-term deficit, equilibrium or recharge trend. The slope and significance of the change can be quantified by performing linear regression on precipitation and temperature series with a timescale of ten years or more, and calibrated in combination with the long-term simulation results of the groundwater numerical model. In this embodiment, for the North China Plain with severe over-extraction, its characteristics include a significant long-term downward trend; the seasonal cycle characteristics are based on the climate type, precipitation concentration period and snowmelt period of the target area, including the expected phase and amplitude range of the dominant annual seasonal signal. Least square spectral analysis can be used to fit the historical GPS vertical displacement sequence to accurately quantify the amplitude and phase of the annual and semi-annual cycle signals; the short-cycle fluctuation characteristics integrate land use data, irrigation system information and soil moisture satellite products to identify sub-seasonal fluctuation signals with a cycle of several days to several months caused by human and natural events such as high-intensity irrigation, concentrated precipitation or reservoir regulation. Soil moisture anomaly index or irrigation water consumption statistics within the sliding time window can be used. If it exceeds 1.5 times the standard deviation of the historical average for the same period, that is, the set of all data within 15 days before and after the same calendar date corresponding to the current date in at least five complete hydrological years in the past, it is identified as a high-activity period time window.
[0030] In a preferred embodiment, if the current time is within a high-activity period window of short-cycle fluctuations, the resampling is defined as a 3-day mean sequence and a 15-day moving median filter is applied; otherwise, the resampling is defined as a monthly mean sequence and a 6-month moving average filter is applied. Both the moving median filter and the moving average filter use center-aligned sliding windows. For data at the beginning and end of the sequence where the window cannot be centered, a method that keeps the window length unchanged but allows asymmetric forward coverage can be used.
[0031] In this embodiment, compared with the existing technology that uses a fixed frequency to simply smooth GPS time series, this method uses prior hydrological feature analysis to actively adopt a higher sampling rate and a suitable filter during the signal active period, thus preserving the integrity of hydrological signals at all scales to the greatest extent and providing high-quality input for the subsequent accurate separation in S2. At the same time, its time series features are a standardized data for hydrological signal analysis, providing crucial contextual basis that is lacking in traditional methods for the uncertainty quantification assessment of each station's data for hydrological applications in S3.
[0032] Specifically, in S2, the first time-series feature is decomposed into multiple time-scale components and a trend residual term by using ensemble empirical mode decomposition.
[0033] It should be noted that the first time-series features of each observation station after S1 processing are decomposed into multiple scales. Empirical mode decomposition is preferred as the decomposition method because it is suitable for non-stationary and nonlinear signals and can effectively alleviate the mode mixing problem. White noise with finite amplitude is added to the input time series, and multiple empirical mode decompositions are performed. Then, all decomposition results are ensembled and averaged to finally obtain a set of intrinsic mode functions arranged from high frequency to low frequency and a residual term representing the overall trend. Each intrinsic mode component represents an oscillation mode at a specific time scale. The termination condition for ensemble empirical mode decomposition follows this rule: decomposition stops when the residual term becomes a monotonic function or the function value has fewer than 3 extreme points. In this embodiment, to eliminate high-frequency pseudo-modes dominated by noise, only intrinsic mode components with a variance contribution rate exceeding 1% of the total variance are selected for subsequent S201-S205 analysis processes. Intrinsic mode IMF1 can correspond to short-term fluctuations of several days to several weeks, intrinsic mode IMF2 and intrinsic mode IMF3 can dominate seasonal changes, while higher-order intrinsic mode components and residuals may reflect interannual or even long-term trends.
[0034] In this embodiment, in the ensemble empirical mode decomposition, the amplitude of the added white noise is usually set to 0.1 to 0.2 times the original time series standard deviation to effectively perturb the signal without excessively distorting its inherent structure; the number of ensemble averaging is set to 100 to 200 times to ensure the statistical stability of the decomposition results, and finally a set of stable intrinsic mode components is obtained.
[0035] In one possible implementation, identifying and extracting the second time-series feature includes: using the non-hydrological deformation time series simulated by the hydrophysical model as an identification constraint; wherein the hydrophysical model includes at least a surface load model for simulating atmospheric and ocean load effects, and a terrestrial hydrological model for simulating soil water and snowmelt changes; performing correlation analysis and phase comparison between the non-hydrological deformation time series and each intrinsic modal component obtained through multi-scale decomposition, identifying and removing modal components associated with non-hydrological driving sources; and reconstructing the remaining modal components to obtain the second time-series feature, namely the pure hydrological load deformation component.
[0036] It should be noted that the hydrophysical model is used to calculate the theoretical vertical displacement of the Earth's surface caused by loads other than hydrology at the same location and within the same time period. The non-hydrological deformation time series obtained through it includes: acquiring high spatiotemporal resolution atmospheric reanalysis data, ocean non-tidal model data, and land surface hydrological model data for the target area; using the Green's function integral method based on the elastic Earth load theory, calculating the mass changes of each load as the theoretical vertical displacement time series at the observation station location; and vector superimposing the displacement time series caused by all loads to generate a comprehensive non-hydrological deformation time series that is perfectly aligned with the GPS observations in terms of timestamps and spatially precisely corresponds to each station.
[0037] In this embodiment, the Green's function integration method uses a table of vertical displacement Green's functions calculated based on the PREM Earth model, which includes, but is not limited to, load types such as atmospheric, oceanic, and terrestrial water. Specifically, recognized Earth load calculation software packages or function libraries such as REAR, GBIS, or LOADDEF are called, and the mass distribution grid data of each load is input. The theoretical displacement time series of each observation station location is obtained through spatial convolution calculation.
[0038] The specific procedures for performing correlation analysis and phase alignment are as follows: S201: For the i-th intrinsic mode component obtained from the target temporal decomposition Calculate the Pearson correlation coefficient between it and the corresponding components of the non-hydrological deformation time series. ; S202: To Perform a Fourier transform on the model signal to calculate its amplitude squared coherence near the preset hydrological cycle frequency. The model signal is obtained by performing the same ensemble empirical mode decomposition on the integrated non-hydrological deformation time series, and is consistent with the observed signal. The k-th model intrinsic mode component with the closest time scale is denoted as In calculating amplitude squared coherence At the same time, calculate and Phase difference at hydrological periodic frequency points ; It should be noted that the hydrological cycle frequency points are derived from the prior hydrological cycle characteristics obtained in S1. Specifically, the typical cycles identified in the prior hydrological cycle characteristics are converted into corresponding frequency values, and a reasonable bandwidth range of ±10% is set as the key frequency band for coherence calculation and analysis, centered on these frequency values. Typical cycles can be represented as an annual cycle of 365 days, a semi-annual cycle of 182.5 days, an irrigation cycle of 30-90 days, etc. The timescales are closest, through comparison of observations. average period With all modes average period The average period can be determined by averaging the instantaneous frequencies obtained through Hilbert transform, and the choice should be made to minimize the absolute value. The smallest model mode as ; In this embodiment, the average period is calculated as follows: For each intrinsic mode component, its analytic signal is obtained by performing a Hilbert transform, and then the instantaneous frequency is obtained by calculating the derivative of its instantaneous phase. Outliers in the instantaneous frequencies are removed, and the reciprocal of the median of the remaining valid instantaneous frequencies is taken as the average period T of the IMF. The quantification standard for outliers is as follows: Let the instantaneous frequency sequence be... Calculate its first quartile. and the third and fourth quartiles ; will satisfy or The values of are considered outliers and are removed. The interquartile range is 1 / 4; the remaining data points constitute the set of effective instantaneous frequencies. S203: Pre-set strong correlation threshold and weak correlation threshold ; In this embodiment, =0.7, =0.3; S204: Perform identification and removal: If and At the key frequency point > 0.8: Meanwhile, if | If | < 30°, then the two are considered to be in phase. If identified as a non-hydrological mode, it will be removed. If | If the angle is greater than 150°, the two modes are considered out of phase, indicating a discrepancy. Therefore, this mode is listed as a mixed mode to be analyzed instead of being directly rejected, and is transferred to the analysis process. The same frequency domain adaptive filtering process is used to process the cases; If 30°≤| If |≤150°, then for the sake of conservatism, it is assumed to be treated as | Handling situations where the angle is >150°; if Then determine For mixed modes of hydrological and non-hydrological signals, frequency domain adaptive filtering is employed: for mixed modes and model signals Perform Fourier transforms separately; in the frequency domain, for each discrete frequency point... If the coherence between the two at that point Then for Complex amplitude at this frequency point Attenuation: ,in This is represented as the preset attenuation intensity coefficient; In this embodiment, the value is set to 0.5; for the attenuated spectrum Perform an inverse Fourier transform to obtain the filtered time-domain signal; if Then determine Hydrological signals, in particular, should be preserved in their entirety. S205: Perform reconstruction: Directly perform linear superposition of all intrinsic mode components and trend residuals that were fully retained or attenuated as determined by S204; as a preferred option, to optimize the signal-to-noise ratio of the final output signal, each retained component can be evaluated based on its variance before superposition. The weighting is performed, and the weights are:
[0039] in, , They are respectively represented as components , Variance over the entire analysis period This represents the total number of intrinsic mode components that are fully retained or retained in attenuated form as determined by S204. The summation applies to all retained components to generate the final second time series feature, which is the vertical displacement time series that corresponds completely to the timestamp of the first time series feature but has been filtered out for non-hydrological interference. Its time series is indexed by station, with each station corresponding to a time series vector, which at least includes the central period, amplitude, variance contribution rate, and filtered metadata of each retained mode, providing a quantitative basis for subsequent analysis.
[0040] It should be noted that the calculations in S202... This refers to the coherence at multiple discrete frequency points within a preset hydrological frequency band. The average value is used for overall correlation assessment; The S204 mid-frequency domain adaptive filtering is based on This refers to the coherence sequence directly calculated after the Fourier transform, which corresponds one-to-one with the frequency f.
[0041] In this embodiment, the method of the present invention was applied to an over-extraction area in the North China Plain and compared with the traditional GPS-GRACE fusion method with unoptimized sampling and simple filtering. The average correlation coefficient between the hydrological signal separated by the traditional method and the measured groundwater well data was 0.65, and significant atmospheric pressure change related signals could still be detected in its residuals. However, the average correlation coefficient between the pure hydrological load deformation component obtained by the present invention and the well data was increased to over 0.82, and the correlation between its residual spectrum and the output of the atmospheric and soil water models was significantly reduced, indicating that non-hydrological interference was effectively removed.
[0042] Specifically, in S3, a quantization matrix is constructed for each observation station gpsi. This is used to assess the uncertainty of groundwater changes reflected by each observation station. Its dimension is: for a time series with t time points, Let be a t×t matrix; its elements are: the elements on the diagonal of the matrix are each time point. Total error variance It consists of the sum of the variances of three parts of the error:
[0043] in, This is represented as the inherent error between observation and solution. This represents the error introduced during the signal separation process. Represented as hydrological characterization error; off-diagonal elements:
[0044] in, Represented as , If expressed as the characteristic time decay constant, then it represents the time point. and The error covariance between the two is used to quantify the correlation strength of errors at different time points.
[0045] It should be noted that in the error model of this embodiment, inherent errors in observation and solution are assumed. Signal separation process error and hydrological characteristic error The three variables are uncorrelated, therefore the total variance is the sum of their variances; inherent errors in observation and solution. variance The main manifestation is high-frequency random noise, which is approximated by calculating the variance of the pure hydrological load deformation component in the high-frequency band of the S2 output; errors introduced during the signal separation process. variance It consists of two parts; One is the variance of the decomposition residual, which is the variance of the remaining residual terms after the set empirical mode decomposition, representing the incompleteness of the decomposition algorithm in interpreting the original signal. The second is the variance of model constraint loss, which is estimated by calculating the variance between the modal components associated with the removed non-hydrological driving sources and the corresponding observed modes, reflecting the inaccuracy of the physical model; The calculation process for model constraint mismatch variance is as follows: The non-hydrological deformation time series from the forward simulation undergoes the same ensemble empirical mode decomposition to obtain a set of model IMFs. For each IMF in the observation time series, its cross-correlation coefficient with all model IMFs is calculated, and the model IMF with the highest cross-correlation coefficient is identified as the corresponding observation mode. The variance of the difference sequence between this observation IMF and its corresponding model IMF is the error variance component introduced by this mode mismatch. For each observation IMF that finds a corresponding model IMF with a cross-correlation coefficient exceeding a preset threshold of 0.5, its mode mismatch error variance is calculated. The mode mismatch error variances of all such observation IMFs are summed to obtain the variance. The constraint mismatch part in the model; Hydrological Characterization Error variance Deformation is interpreted as the physical uncertainty of groundwater changes, and a scaling factor is set based on prior knowledge for estimation. In this embodiment, according to regional hydrogeological studies, the scaling factor between deformation caused by groundwater load and water volume changes is set to have an uncertainty range. Then, this uncertainty is propagated to the deformation sequence using the error propagation law to estimate the resulting variance.
[0046] In this embodiment, the high-frequency band specifically refers to: for each intrinsic mode component obtained by ensemble empirical mode decomposition in S2, a fast Fourier transform is performed on each IMF component to calculate its power spectrum, and the intrinsic mode components with a period of less than 30 days corresponding to the peak of the power spectrum are identified as high-frequency components. That is, it is estimated by calculating the sum of the variances of these intrinsic mode components that are identified as high-frequency components; When applying the error propagation law to propagate to the deformation sequence, let the deformation sequence be S, and the estimated groundwater change be W. The theoretical relationship between the two is S = ks × W, where... Expressed as a proportionality coefficient, under the assumption that the error of W is independent of the error of ks, by Uncertainty The additional variance of the resulting deformation sequence S can be estimated as follows:
[0047] in, Represented as the deformation value at time point ct, ( The uncertainty range of the proportionality coefficient ks in this embodiment is as follows: =±15%, which is based on the summary of the statistical range of this parameter in existing borehole pressure water test reports, indoor geotechnical mechanics test results, or published academic literature on the target area or similar hydrogeological units.
[0048] It should be noted that the characteristic time decay constant The autocorrelation function is obtained by fitting the sample autocorrelation function of the residual sequence of the pure hydrological load deformation component; specifically, after S2, the residual between the observed sequence and the reconstructed pure hydrological sequence is calculated; the autocorrelation coefficient of the residual sequence under different time lags is calculated; and the exponential decay function is used. The calculated autocorrelation coefficient is fitted using nonlinear least squares fitting, and the resulting time constant is... If the residual is close to white noise, then The value approaches 0.
[0049] Furthermore, to support observation network optimization and multi-site data fusion, it is also necessary to construct the spatial covariance between the errors of any two observation sites. In this embodiment, the spatial covariance is assumed to decay with the distance between sites, and an exponential model is used to represent the spatial covariance at the same time. Covariance between the errors of the two stations:
[0050] in, and The sites are gpsi and gpsj respectively. The average value over the entire time series is used to represent the typical error level at that site. Let L represent the spherical distance between the two stations, and let L represent the spatial correlation length. This represents the maximum correlation coefficient at zero distance. In this embodiment, the spatial correlation length L is determined based on the hydrogeological map: for homogeneous porous aquifers with a scale of 50-100 km, L can be set to 20-50 km; for blocks divided by large faults, L is considered to be less than 5 km. Maximum correlation coefficient This reflects the spatial homology of errors and is typically set empirically between 0.7 and 0.9. In this embodiment, it is taken as... =0.8; L can be jointly optimized and determined by fitting the cross-correlation coefficient of the error sequences between different site pairs as a function of distance.
[0051] In this embodiment, taking an inland basin as an example in the simulation experiment, the quantization matrix generated by S3 is used. By optimizing the network, under the constraint of the same number of newly added sites, the posterior mean square error of the spatial details of groundwater changes within the basin obtained by inversion can be reduced by about 20% to 40% compared with the traditional empirical weighted method. In particular, the physical rationality of the inversion results is significantly improved near the hydrogeological boundary.
[0052] Specifically, in S4, the network layout of each observation station is simulated and optimized, including: S401: Construct an optimization function with the objective of minimizing the posterior uncertainty of the inversion results of regional groundwater storage changes; In this embodiment, the optimization function is specifically to minimize the trace of the posterior covariance matrix of the estimated regional groundwater storage changes, derived from Bayesian inversion theory. ,Right now: Where JC represents the sensitivity matrix for groundwater parameters. The covariance matrix of the observed data is represented by the quantization matrix output by S3. It is represented as the prior covariance matrix of groundwater parameters; It should be noted that the construction of the observation sensitivity matrix JC includes: JC is a (jm×jn) dimensional matrix, where jm is the total number of observations at all observation points within the total time series, and jn is the total number of groundwater parameter grids after discretization of the inversion region; the matrix element JC(mi,mj) represents the theoretical response of the mi-th observation to the change in groundwater mass per unit mass in the mj-th grid cell, and its calculation is based on the elastic earth load theory, specifically implemented through the loaded Green's function method: ,in, It is the horizontal distance from the center point of grid mj to the observation point mi. The corresponding vertical displacement Green's function value can be calculated using standard Earth models such as REAR and PREM. SA represents the area of the grid mj, used to extend the point load response to a surface load response; the covariance matrix of the observation data... The construction includes: It is a (cm×cm) dimensional block diagonal matrix. For the existing KN real observation stations, the quantization matrix of the uncertainty of the q-th station output by S3 is used as a sub-block, arranged in the order of the stations. On the diagonal; the prior covariance matrix of the parameters The construction includes: It is a (cn×cn) dimensional matrix, representing the prior uncertainty and correlation between groundwater variations in each grid cell. It is set as follows: It is a diagonal matrix, that is Among them, prior variance It can be set based on statistical data on historical groundwater fluctuations in the study area; In a preferred example, the prior variance The determination can be achieved by: collecting groundwater level change data from representative monitoring wells in the target area over the past 10 years, converting it into equivalent water height change, calculating the variance of this time series data, and multiplying this variance by a safety factor as the basis for the determination. The value of is chosen to reflect a conservative estimate of the uncertainty of groundwater variation based on historical experience; In this embodiment, when implementing the Green's function method, the standard Earth model parameters, the average crustal elastic parameters of the target area, and the geographic coordinates of all grid center points and observation points can be input by calling the geophysical software library. The vertical displacement Green's function value between each pair (observation point, grid) can be calculated in batches. This function value is multiplied by the area SA of the corresponding grid cell to fill in the complete observation sensitivity matrix JC. S402: Based on the hydrogeological unit boundaries, fault lines, and the spatial distribution map of uncertainties of existing observation stations in the target area, a set of virtual observation candidate points is generated in key heterogeneous areas and uncertain cavity areas; It should be noted that the key heterogeneous region refers to the area within a certain buffer zone on both sides of the boundary of the hydrogeological unit; the uncertainty cavity region refers to the area that is more than a certain threshold distance from the existing observation station and whose uncertainty propagation value obtained by Kriging interpolation is higher than the regional average. In this embodiment, the buffer zone of the key heterogeneous region is taken as 5-10 kilometers on both sides of the boundary line of the main aquifer; the specific threshold distance is taken as 1.5 to 2 times the average spacing; S403: Using a sequential greedy algorithm, the virtual observation point that causes the function value of the optimization function to decrease the most is iteratively selected from the set of virtual observation candidate points and added to the optimization network; Furthermore, the sequential greedy algorithm starts from the existing network and executes in each iteration: Suppose the current optimization network contains cs points, and the corresponding observation vector dimension is... The observed covariance matrix is Based on this, the current posterior covariance matrix is calculated:
[0053] and its traces ,in, This is represented as the sensitivity matrix corresponding to the current cs points; for each candidate point cs: create a temporary observation covariance matrix. and the corresponding posterior covariance ; Calculate the reduction in the objective function: ; Choose to make The largest candidate point is moved from the candidate pool to the network; this process is repeated until the number of new sites reaches a preset limit, or the maximum number of new sites added in a single iteration. Less than the set threshold threshold It is a small positive number used to control the stopping of the optimization iteration; S404: Output the optimized observation network scheme and its corresponding spatial constraint weight map, where the spatial constraint weight map represents the quantitative difference in the optimized network's ability to invert groundwater at different locations within the region. It should be noted that in the spatial constraint weight diagram, the weight value of each grid cell Gri is... The calculation is expressed as:
[0054] in, and These are the posterior variances of the estimated groundwater changes in the grid cells (Gri) before and after network optimization, respectively. This can be directly derived from... Extracting the diagonal elements of a matrix.
[0055] In this embodiment, compared with the method of relying solely on experience to deploy or uniformly densify the GPS network, the method in a simulation experiment of a groundwater over-extraction area in Hebei Province, based on the spatiotemporal correlation uncertainty output by S3 and after adding 5 virtual stations through algorithm optimization by S4, significantly reduced the average posterior standard deviation of the groundwater inversion results for the entire area. This is significantly better than the reduction brought about by randomly deploying the same number of stations outside the uncertainty cavity area. This directly proves that the proposed scheme has clear targeting and high efficiency in optimizing the observation network.
[0056] Specifically, in S5, the third temporal feature is a time-series, high spatial resolution grid data of groundwater storage changes, where each grid value represents the change in groundwater storage of that cell relative to the long-term mean at a specific time.
[0057] Furthermore, the spatial constraints include at least: the second temporal characteristics of each observation station, namely, the purified vertical displacement sequence reflecting the changes in groundwater at the station location; a quantization matrix describing the error variance and temporal correlation of the data at different time points of each station; a spatial constraint effectiveness map; and regional total surface constraints: equivalent water high-grid data of monthly land water storage changes in the target area obtained from GRACE / GRACE-FO satellite mission data.
[0058] In one possible implementation, a high-resolution groundwater storage variation field is solved. The third time series feature itself constitutes the joint inversion equation set:
[0059]
[0060]
[0061] in, Represented as a Green's function matrix based on the elastic Earth model, used to... The theoretical vertical displacement mapped to the location of each observation station. The second time series feature is represented by H, which is the spatial averaging operator matrix used to... The spatial resolution is aggregated to a scale that matches the satellite gravity field observation data. Represented as satellite gravity field variation data, LC represents the regularization operator matrix used to stabilize the solution in regions with insufficient data, and λ is the regularization parameter. and It is represented as a weight matrix.
[0062] It should be noted that the construction process of the Green's function matrix G is as follows: Based on the PREM Earth model, the point load Green's function is calculated using the well-known surface load deformation theory in this field; each high-resolution grid rule is discretized into several smaller computational sub-units; for each observation station and each high-resolution grid, the vertical displacement caused by the center point of all sub-units in the grid at the station is calculated, and its arithmetic mean is taken as the average influence of the grid on the station, which is the matrix element; Each element of the spatial averaging operator matrix H represents the area contribution weight of the coarse-resolution GRACE grid corresponding to each high-resolution grid, which is usually determined according to the grid area overlap ratio; the regularization operator matrix LC is usually a first-order or second-order finite difference matrix used to express the smoothness constraint of the solution in space, and λ is the regularization parameter, which can be objectively determined by the L-curve method. The specific operations include: selecting a suitable number of candidate values for λ, solving the minimization problem of the objective function Φ once for each candidate λ; and recording the data fitting residual norm obtained from each solution. Norm of reconciliation constraint terms In a log-log coordinate system, a scatter plot is drawn with the data residual norm as the horizontal axis and the constraint norm as the vertical axis, forming an L-shaped curve. The λ value corresponding to the bend of the curve is selected as the optimal parameter; The construction includes: for each observation station, its basic weights are derived from the quantization matrix output by S3. The inverse ensures that sites with lower uncertainty have higher weights in the inversion, achieving objective and quantitative precision control. The spatial constraint weight map values generated by S4 are used as a spatial modulation factor and multiplied by the corresponding site weights. For sites located at geological boundaries and in key areas revealing spatial heterogeneity, even if their absolute precision is slightly lower, their weights are relatively increased to strengthen the targeted characterization of heterogeneous features. Ultimately, It is a block diagonal matrix composed of weighted submatrices of each site; The construction includes: the weights are mainly determined based on the nominal error variance of the GRACE data itself, while the representativeness error caused by the inconsistency between the data and the GPS observation scale can also be considered.
[0063] Furthermore, the joint inversion equations are solved by minimizing the following objective function: , in, This is represented as a scaling factor balancing the relative weights of GPS and GRACE constraints, initially set based on the prior variance ratio of the two types of data; in this embodiment, the scaling factor... The GPS observation residuals are objectively determined as follows: Before the formal inversion, the GPS observation residuals are calculated using historical or training data of the target area. Compared with GRACE observation residuals variance estimate and ;in, For a known reference field, set This ensures that, in the objective function Φ, the two types of data constraint terms have comparable contributions due to their inherent noise levels; reference field This is achieved by aggregating the daily average groundwater storage changes simulated by the hydrological model of the target area into a monthly average over time, and then sampling it spatially to the same grid resolution as ΔG.
[0064] In another aspect, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the method described above.
[0065] In another aspect, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the method described above.
[0066] In another embodiment provided in this application, a computer program product containing instructions is also provided, which, when run on a computer, causes the computer to execute any of the time-constrained groundwater spatial distribution detection methods described in the above embodiments.
[0067] It is understood that the systems, devices, and storage media provided in the embodiments of the present invention correspond to the methods provided in the embodiments of the present invention, and the explanations, examples, and beneficial effects of the relevant content can be referred to the corresponding parts of the above methods.
[0068] In the above embodiments, implementation can be achieved, in whole or in part, through software, hardware, firmware, or any combination thereof. When implemented in software, it can be implemented, in whole or in part, as a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transferred from one computer-readable storage medium to another.
[0069] For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.). The computer-readable storage medium can be any available medium that a computer can access, or a data storage device such as a server or data center that integrates one or more available media.
[0070] The available media may be magnetic media (e.g., floppy disks, hard disks, magnetic tapes), optical media (e.g., DVDs), or semiconductor media (e.g., solid state disks (SSDs)).
[0071] It should be noted that in this document, relational terms such as first and second are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations.
[0072] Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.
[0073] The various embodiments in this specification are described in a related manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions of the method embodiments.
[0074] The embodiments of the present invention are given for the purposes of illustration and description. Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.
Claims
1. A method for detecting the spatial distribution of groundwater based on time-series constraints, characterized in that, Perform the following steps using a computer device: S1: Obtain vertical displacement data of observation stations in the target area to form first time series data, and preprocess the first time series data to generate first time series features aligned with hydrological features; S2: Perform multi-scale decomposition on the first time series features, and use non-hydrological deformation signals as constraints to identify and extract the pure hydrological load deformation components caused by groundwater changes, so as to obtain the second time series features; S3: Based on the extraction process of the second time series feature, the uncertainty of groundwater changes reflected by each observation station is quantitatively assessed and summarized as the overall uncertainty; S4: With the goal of minimizing overall uncertainty, the network layout of each observation station is simulated and optimized to generate the optimized observation network and its corresponding spatial constraint weights; S5: Using the second temporal feature and the spatial constraint weight as spatial constraint conditions, solve the groundwater storage change distribution field to generate a spatiotemporally continuous third temporal feature as the result of groundwater spatial distribution detection.
2. The groundwater spatial distribution detection method based on time-series constraints as described in claim 1, characterized in that: The first temporal feature in S1 specifically includes: analyzing the prior hydrological periodic characteristics of groundwater dynamics at the observation stations in the target area; The prior hydrological periodic characteristics include at least interannual trends, seasonal cycles, and short-cycle fluctuations caused by artificial irrigation activities; wherein, during the active period of the short-cycle fluctuations, high-frequency deformation signals related to groundwater changes are captured. Based on the prior hydrological periodic characteristics, a resampling strategy and filtering parameters matching different hydrological activity periods are dynamically determined to generate a first time-series feature aligned with the hydrological features.
3. The groundwater spatial distribution detection method based on time-series constraints as described in claim 1, characterized in that: In step S2, non-hydrological deformation signals are used as constraints to identify and extract the second time-series feature, specifically including: The non-hydrological deformation time series obtained from the hydrophysical model simulation are used as identification constraints; The hydrophysical model includes at least a surface load model for simulating atmospheric and ocean load effects, and a terrestrial hydrological model for simulating changes in soil water and snowmelt. The non-hydrological deformation time series is correlated with each intrinsic mode component obtained by multi-scale decomposition and phase comparison is performed to identify and remove the mode components associated with the non-hydrological driving source; the remaining mode components are reconstructed to obtain the second time series feature, namely the pure hydrological load deformation component.
4. The groundwater spatial distribution detection method based on time-series constraints as described in claim 3, characterized in that: The process of performing correlation analysis and phase comparison between the non-hydrological deformation time series and the intrinsic mode components obtained through multi-scale decomposition specifically includes: S201: For the i-th intrinsic mode component obtained from the target temporal decomposition Calculate the Pearson correlation coefficient between it and the corresponding component of the non-hydrological deformation time series. ; S202: To Perform a Fourier transform on the model signal to calculate its amplitude squared coherence near the preset hydrological cycle frequency. The model signal is related to The k-th model intrinsic mode component with the closest time scale is denoted as In calculating amplitude squared coherence At the same time, calculate and Phase difference at hydrological periodic frequency points ; S203: Pre-set strong correlation threshold and weak correlation threshold ; S204: Perform identification and removal: like and At key frequencies > 0.8: At the same time, if | If | < 30°, then the two are considered to be in phase. If it is confirmed to be a non-hydrological mode, it will be removed; if | If the angle is greater than 150°, then the two are determined to be out of phase, and the operation is switched to AND. The same frequency domain adaptive filtering process is used to process the cases; like Then determine For mixed modes of hydrological and non-hydrological signals, frequency domain adaptive filtering is used for attenuation: like Then determine Hydrological signals, in particular, should be preserved in their entirety. S205: Perform the reconstruction: Directly superimpose all intrinsic mode components and trend residuals that were determined to be fully retained or retained in attenuated form as determined in S204 to generate the second time series feature.
5. The method for detecting the spatial distribution of groundwater based on time-series constraints as described in claim 4, characterized in that: In S204, the determination is made. For mixed modes of hydrological and non-hydrological signals, frequency-domain adaptive filtering is used for attenuation, specifically including: right and model signals Perform Fourier transforms on each; In the frequency domain, for each discrete frequency point If the coherence between the two at that point Then for Complex amplitude at this frequency point Attenuation: in, Represented as a preset attenuation intensity coefficient; for the attenuated spectrum Perform an inverse Fourier transform to obtain the filtered time-domain signal.
6. The method for detecting the spatial distribution of groundwater based on time-series constraints as described in claim 5, characterized in that: In S3, a quantization matrix is constructed for each observation station gpsi. This is used to assess the uncertainty of groundwater changes reflected by each observation station, specifically expressed as: For a time series with t time points It is a t×t matrix; The elements on the diagonal of the matrix represent each time point. Total error variance It consists of the sum of the variances of three parts of the error: in, This is represented as the inherent error between observation and solution. This represents the error introduced during the signal separation process. This is represented as hydrological characterization error; Off-diagonal elements This indicates a point in time. and The error covariance between the two is used to quantify the correlation strength of errors at different time points.
7. The method for detecting the spatial distribution of groundwater based on time-series constraints as described in claim 1, characterized in that: The simulation optimization of the network layout of each observation station in S4 specifically includes: An optimization function is constructed with the objective of minimizing the posterior uncertainty of the inversion results of regional groundwater storage changes. Based on the hydrogeological unit boundaries, fault lines, and spatial distribution map of uncertainties of existing observation stations in the target area, a set of virtual observation candidate points is generated in key heterogeneous areas and uncertainty cavity areas. A sequential greedy algorithm is used to iteratively select virtual observation points from the set of virtual observation candidate points that can reduce the function value of the optimization function the most, and add them to the optimization network; The optimized observation network scheme and its corresponding spatial constraint weight map are output, wherein the spatial constraint weight map represents the quantitative difference in the optimized network's ability to constrain groundwater inversion at different locations within the region.
8. The groundwater spatial distribution detection method based on time-series constraints as described in claim 1, characterized in that: The third temporal feature in S5 is a time-series, high spatial resolution grid data of groundwater storage changes, where each grid value represents the change in groundwater storage of the corresponding pixel at a specific time relative to the long-term average.
9. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by a processor, it causes the processor to perform the steps of the method as described in any one of claims 1 to 8.
10. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the computer program is executed by the processor, it causes the processor to perform the steps of the method as described in any one of claims 1 to 8.