Underground water estimation method based on ensemble Kalman filtering

Abstract

The invention relates to a groundwater estimation method based on ensemble Kalman filtering, and the method comprises the steps: presetting the data distribution of each attribute based on a target region, firstly obtaining the water level value and boundary information of each water body region corresponding to each preset historical time point, and then constructing the surface water reserve grid distribution of the target region corresponding to each preset historical time point; and finally, on the basis of ensemble Kalman filtering, calculating groundwater reserve estimation values of all grid points corresponding to each preset historical time point in the target area, and realizing estimation of groundwater reserves in the target area. The design scheme can accurately estimate groundwater reserve changes, is helpful to better understand the climate change rate and the spatial-temporal characteristics of groundwater reserves in water resource shortage areas under the influence of human activities, can provide support for sustainable groundwater management, water resource planning and distribution decision, and is helpful to determine the groundwater resource planning and distribution decision caused by climate change or over-mining. The potential area may face underground water resource risks.

Description

Technical Field

[0001] The invention relates to a groundwater estimation method based on ensemble Kalman filtering, and belongs to the technical field of hydrogeology. Background Art

[0002] In areas with water shortages, the distribution and utilization of water resources face severe challenges. Groundwater resources are not only the key to supporting the survival and development of local societies, but also an important pillar to alleviate water shortages, promote agricultural production and economic development. However, with the increasing impact of global climate change, coupled with population growth, accelerated industrialization, and unreasonable water resource management, groundwater systems in water-scarce areas are under unprecedented pressure. Changes in precipitation patterns, rising temperatures, and frequent extreme weather events caused by climate change have directly affected groundwater recharge and storage. At the same time, excessive groundwater extraction to meet the growing population demand and agricultural production has led to a decline in groundwater levels and an increasing risk of water depletion. This situation not only threatens the safety of local residents' domestic water supply, but also poses a major risk to agricultural output, ecosystem health, and regional economic development. However, traditional groundwater reserve estimation methods still have many limitations, such as the inability to accurately reflect the overall regional hydrological dynamics and the temporal characteristics of groundwater recharge. The lack of accurate groundwater status information also greatly restricts water resource planning and management. Summary of the invention

[0003] The technical problem to be solved by the present invention is to provide a groundwater estimation method based on ensemble Kalman filtering, which can accurately estimate the temporal and spatial changes of groundwater reserves and provide support for water resources management decisions.

[0004] In order to solve the above technical problems, the present invention adopts the following technical solution: The present invention designs a groundwater estimation method based on ensemble Kalman filtering to achieve the estimation of groundwater reserves in the target area, comprising the following steps:

[0005] Step A. Obtain the grid distribution of terrestrial water storage, soil moisture content, snow water equivalent, canopy water storage, precipitation, and evapotranspiration at preset historical time points in the preset historical time period corresponding to the preset target spatial grid resolution in the target area, and obtain the water level and boundary information of each water body area in the target area corresponding to each preset historical time point in the preset historical time period, and then proceed to step B;

[0006] Step B. Based on the boundary information of each water body area in the target area corresponding to each preset historical time point in the preset historical time period, according to the area of ​​each water body area in each grid corresponding to each preset historical time point at the target spatial grid resolution, combined with the water level value of each water body area corresponding to each preset historical time point, calculate and obtain the surface water reserves of each water body area in each grid corresponding to each preset historical time point at the target spatial grid resolution, and then define that the surface water reserves of other grids at the target spatial grid resolution corresponding to each preset historical time point are equal to 0, construct the surface water reserves grid distribution of the target area at each preset historical time point in the preset historical time period corresponding to the target spatial grid resolution, and then enter step C;

[0007] Step C. Based on the ensemble Kalman filter, the estimated groundwater reserves of all grid points in the target area corresponding to the target spatial grid resolution at each preset historical time point in the preset historical time period are calculated to estimate the groundwater reserves in the target area.

[0008] As a preferred technical solution of the present invention: it also includes steps D to F as follows, after executing step C, proceed to step D;

[0009] Step D. Estimated groundwater storage values ​​of grid points at each preset historical time point t in the preset historical time period according to the target area with respect to the target spatial grid resolution According to the following formula:

[0010]

[0011] Obtain the estimated groundwater storage fluctuation value of the target area at the grid point corresponding to the historical time point t compared with the previous historical time point t-1 at the target spatial grid resolution Among them, 1,…,i,…,I, I represents the number of grids in the target area with respect to the target space grid resolution, It represents the estimated value of underground water storage at the grid point of the target area at the historical time point t corresponding to the i-th grid of the target spatial grid resolution. Represents the estimated value of underground water storage at the grid point of the target area at the target spatial grid resolution corresponding to the historical time point t-1, and then enters step E;

[0012] Step E. Based on the historical analysis time period L, the following formula is used:

[0013]

[0014] Obtain the potential recharge R of groundwater storage in the target area corresponding to the historical analysis period L L ;

[0015] At the same time, according to the following formula:

[0016]

[0017] Obtain the potential discharge D of underground water storage in the target area corresponding to the historical analysis period L L ,|| means taking the absolute value;

[0018] Then proceed to step F;

[0019] Step F. According to the following formula:

[0020] R net,L =R L -D L

[0021] Obtain the net recharge R of groundwater storage in the target area corresponding to the historical analysis period L net,L .

[0022] As a preferred technical solution of the present invention: the step A includes the following steps A1 to A4;

[0023] Step A1. Obtain the grid distribution of terrestrial water storage, soil moisture content, snow water equivalent, canopy water storage, precipitation, and evapotranspiration at each preset historical time point in the preset historical time period corresponding to the target area, and obtain the water level and boundary information of each water body area in the target area at each preset historical time point in the preset historical time period, and then proceed to step A2;

[0024] Step A2. Determine whether there are missing water level values ​​at each preset historical time point in each water body area in the target area corresponding to the preset historical time period, and if so, fill in the missing water level values; otherwise, proceed to step A3;

[0025] Step A3. Update and unify the grid distribution of soil moisture content, snow water equivalent, canopy water storage, precipitation, evapotranspiration, and terrestrial water storage at each preset historical time point in the preset historical time period corresponding to the target area to the preset target spatial grid resolution, and then proceed to step A4;

[0026] Step A4. Determine whether there is a gap in the grid distribution of land water reserves at each preset historical time point in the target area corresponding to the preset historical time period. If so, fill in the missing land water reserves; otherwise, proceed to step B.

[0027] As a preferred technical solution of the present invention: in the step A2, when it is determined that there are missing water level values ​​corresponding to each preset historical time point in each water body area in the target area within the preset historical time period, the following steps A2-1 to A2-4 are executed to fill in the missing water level values;

[0028] Step A2-1. For each non-missing water level value, perform the following steps:

[0029] According to the boundary information of the water area corresponding to the non-missing water level value at the corresponding historical time point, each grid under the precipitation grid distribution involved in the water area based on the boundary information is obtained, and the average value of the precipitation of each grid at the historical time point corresponding to the non-missing water level value is obtained, so as to constitute the precipitation of the local area of ​​the target area where the water area corresponding to the non-missing water level value is located at the corresponding historical time point;

[0030] At the same time, according to the boundary information of the water area corresponding to the non-missing water level value at the corresponding historical time point, each grid under the evapotranspiration grid distribution involved in the water area based on the boundary information is obtained, and the average value of the evapotranspiration of each grid at the historical time point corresponding to the non-missing water level value is obtained, so as to constitute the evapotranspiration of the local area of ​​the target area where the water area corresponding to the non-missing water level value is located at the corresponding historical time point;

[0031] That is, the precipitation and evapotranspiration corresponding to the non-missing water level value are obtained, and then the precipitation and evapotranspiration corresponding to each non-missing water level value are obtained, and the non-missing water level value and the corresponding precipitation and evapotranspiration constitute a first sample, and each first sample is obtained, and then step A2-2 is entered;

[0032] Step A2-2. Based on each first sample, with the precipitation and evapotranspiration in the first sample as input and the non-missing water level value in the first sample as output, train the LSTM neural network to obtain a water level value neural network model, and then proceed to step A2-3;

[0033] Step A2-3. For each missing water level value, perform the following steps:

[0034] According to the boundary information of the water area corresponding to the missing water level value at the corresponding historical time point, each grid under the precipitation grid distribution involved in the water area based on the boundary information is obtained, and the average value of the precipitation of each grid at the historical time point corresponding to the missing water level value is obtained, so as to constitute the precipitation of the local area of ​​the target area where the water area corresponding to the missing water level value is located at the corresponding historical time point;

[0035] At the same time, according to the boundary information of the water area corresponding to the missing water level value at the corresponding historical time point, each grid under the evapotranspiration grid distribution involved in the water area based on the boundary information is obtained, and the average value of the evapotranspiration of each grid at the historical time point corresponding to the missing water level value is obtained, so as to constitute the evapotranspiration of the local area of ​​the target area where the water area corresponding to the missing water level value is located at the corresponding historical time point;

[0036] That is, the precipitation and evapotranspiration corresponding to the missing water level value are obtained, and then the precipitation and evapotranspiration corresponding to each missing water level value are obtained, and then step A2-4 is entered;

[0037] Step A2-4. For each missing water level value, take the precipitation and evapotranspiration corresponding to the missing water level value as input, apply the water level value neural network model for processing, obtain the corresponding water level value, fill the missing water level value, and then fill each missing water level value.

[0038] As a preferred technical solution of the present invention: in the step A3, spherical harmonic expansion is performed on the soil moisture grid distribution, snow water equivalent grid distribution, canopy water storage grid distribution, precipitation grid distribution, evapotranspiration grid distribution, and terrestrial water storage grid distribution at each preset historical time point within the preset historical time period corresponding to the target area, to obtain the spherical harmonic model coefficients, and then the spatial resolution of the grid distribution is adjusted and updated to be unified to the preset target spatial grid resolution.

[0039] As a preferred technical solution of the present invention: in the step A3, during the spherical harmonic expansion, a sparse least squares method is used to minimize the residual between the original field and the simulated field, thereby weakening the Gibbs effect in the spherical harmonic expansion.

[0040] As a preferred technical solution of the present invention: in the step A4, when it is determined that there is a missing grid distribution of land water storage at each preset historical time point in the preset historical time period corresponding to the target area, the following steps A4-1 to A4-4 are executed to fill in the missing land water storage grid values;

[0041] Step A4-1. For each non-missing terrestrial water storage grid value, the non-missing terrestrial water storage grid value is combined with the soil moisture grid value, precipitation grid value, and evapotranspiration grid value corresponding to the corresponding historical time point of the grid where the grid is located to form a second sample, obtain each second sample, and then enter step A4-2;

[0042] Step A4-2. Based on each second sample, the soil moisture grid value, precipitation grid value, and evapotranspiration grid value in the second sample are used as input, and the non-missing land water storage grid value in the second sample is used as output, and the LSTM neural network is trained to obtain a land water storage grid value neural network model, and then proceed to step A4-3;

[0043] Step A4-3. For each missing land water storage grid value, take the soil moisture grid value, precipitation grid value and evapotranspiration grid value of the grid where the missing land water storage grid value is located at the corresponding historical time point as input, apply the land water storage grid value neural network model to process, obtain the corresponding land water storage grid value, fill in the missing land water storage grid value, and then fill in each missing land water storage grid value.

[0044] As a preferred technical solution of the present invention: the step C includes the following steps C1 to C8;

[0045] Step C1. According to the following formula:

[0046]

[0047]

[0048] Define the state vector x of the target area corresponding to the historical time point t t , and the observation vector y of the target area corresponding to the historical time point t t , where 1,…,i,…,I, I represents the number of grids in the target area with respect to the target space grid resolution, It represents the underground water storage grid value calculated by combining the soil moisture grid value, snow water equivalent grid value, canopy water storage grid value, surface water storage grid value, and land water storage grid value at the target area corresponding to the historical time point t at the i-th grid of the target spatial grid resolution. represents the grid value of land water storage in the target area at the historical time point t corresponding to the i-th grid of the target spatial grid resolution, [·] T Indicates transposition, and then goes to step C2;

[0049] Step C2. Construct the state vector x t With the observation vector y t The correlation equation between them is as follows:

[0050] y t =H t x t +v t

[0051] Among them, H t Represents the corresponding historical time point t and state vector xt Coefficient matrix of the same dimension, v t Represents the corresponding historical time point t and state vector x t Observation noise of the same dimension, and then enter step C3;

[0052] Step C3. Based on the preset E different disturbances, the state vector x is calculated for the combined soil moisture content, snow water equivalent, canopy water storage, surface water storage, and terrestrial water storage. t , add each disturbance separately and construct the set X corresponding to the historical time point t t as follows:

[0053]

[0054] Among them, 1,…,e,…,E, Represents the state vector x t Add the initial amount of the e-th disturbance, and then go to step C4;

[0055] Step C4. For the state vector x t-1 Add the initial amount of each disturbance Execute time forward propagation to obtain each The corresponding updated values Then proceed to step C5;

[0056] Step C5. According to the following formula:

[0057]

[0058]

[0059] Construct the collection matrix X t|t-1 and the state covariance update value matrix P t|t-1 ,in, Represents the set matrix X t|t-1 The average value of each element in , and then go to step C6;

[0060] Step C6. According to the following formula:

[0061]

[0062] Calculate the gain matrix K corresponding to the historical time point t t , where R t represents the covariance matrix of the observation vector corresponding to the historical time point t, and then proceeds to step C7;

[0063] Step C7. According to the following formula:

[0064]

[0065] Get each Correction value in, Represents the observation vector y t Add the initial amount of the e-th disturbance, and then enter step C8;

[0066] Step C8. According to the following formula:

[0067]

[0068] Obtain the estimated value of groundwater storage corresponding to historical time point t That is, the estimated groundwater reserves of all grid points in the target area at each preset historical time point t within the preset historical time period corresponding to the target spatial grid resolution are used to estimate the groundwater reserves in the target area.

[0069] The groundwater estimation method based on ensemble Kalman filtering described in the present invention has the following technical effects compared with the prior art by using the above technical solution:

[0070] (1) The present invention designs a groundwater estimation method based on ensemble Kalman filtering. Based on the preset attribute data distribution of the target area, the water level value and boundary information of each water body area corresponding to each preset historical time point are first obtained, and then the surface water storage grid distribution corresponding to each preset historical time point of the target area is constructed. Finally, based on the ensemble Kalman filtering, the estimated groundwater storage value of all grid points in the target area corresponding to each preset historical time point is calculated to estimate the groundwater storage of the target area. The design scheme can accurately estimate the change of groundwater storage, which helps to better understand the spatiotemporal characteristics of groundwater storage in water-scarce areas under the influence of climate variability and human activities, and can provide support for sustainable groundwater management, water resources planning and allocation decisions, and is conducive to identifying potential areas that may face groundwater resource risks due to climate change or overexploitation.

[0071] (2) In the groundwater estimation method based on ensemble Kalman filtering designed by the present invention, spherical harmonic functions are used to unify the spatial resolution of hydrological data, and the sparse least squares method is used to effectively weaken the Gibbs effect. The missing terrestrial water storage and water level data are reconstructed by establishing an LSTM neural network, thereby ensuring the integrity of the groundwater storage estimation results. In addition, the method of estimating groundwater reserves using ensemble Kalman filtering incorporates the uncertainties in the observed data and models, and no longer relies on the prior assumptions and decomposition process of traditional estimation methods, thereby improving the accuracy and reliability of the estimation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 The present invention designs a groundwater estimation method based on ensemble Kalman filtering to realize the principle diagram of groundwater recharge estimation. DETAILED DESCRIPTION

[0073] The specific implementation modes of the present invention will be further described in detail below in conjunction with the accompanying drawings.

[0074] The present invention designs a groundwater estimation method based on ensemble Kalman filtering to realize the estimation of groundwater reserves in the target area. In practical applications, such as Figure 1 As shown, specifically perform the following steps A to C.

[0075] Step A. Based on the GRACE CSR masocn data and the GLDAS-Noah model data, obtain the grid distribution of terrestrial water storage, soil moisture content, snow water equivalent, canopy water storage, precipitation, and evapotranspiration at each preset historical time point in the preset historical time period corresponding to the preset target spatial grid resolution in the target area, and obtain the water level value and boundary information of each water body area in the target area corresponding to each preset historical time point in the preset historical time period, and then enter step B.

[0076] In practical applications, the above step A is specifically designed to execute the following steps A1 to A4.

[0077] Step A1. Obtain the grid distribution of terrestrial water storage, soil moisture content, snow water equivalent, canopy water storage, precipitation, and evapotranspiration at each preset historical time point in the preset historical time period corresponding to the target area, and obtain the water level value and boundary information of each water area in the target area corresponding to each preset historical time point in the preset historical time period, and then proceed to step A2.

[0078] Step A2. Determine whether there are any missing water level values ​​at each preset historical time point in each water area in the target area corresponding to the preset historical time period. If so, execute the following steps A2-1 to A2-4 to fill in the missing water level values; otherwise, proceed to step A3.

[0079] Step A2-1. For each non-missing water level value, perform the following steps:

[0080] According to the boundary information of the water area corresponding to the non-missing water level value at the corresponding historical time point, each grid under the precipitation grid distribution of the water area based on the boundary information is obtained, and the average value of the precipitation of each grid at the historical time point corresponding to the non-missing water level value is obtained, which constitutes the precipitation of the local area of ​​the target area where the water area corresponding to the non-missing water level value is located at the corresponding historical time point.

[0081] At the same time, according to the boundary information of the water area corresponding to the non-missing water level value at the corresponding historical time point, each grid under the grid distribution of the evapotranspiration involved in the water area based on the boundary information is obtained, and the average value of the evapotranspiration of each grid at the historical time point corresponding to the non-missing water level value is obtained, which constitutes the evapotranspiration of the local area of ​​the target area where the water area corresponding to the non-missing water level value is located at the corresponding historical time point.

[0082] That is, the precipitation and evapotranspiration corresponding to the non-missing water level value are obtained, and then the precipitation and evapotranspiration corresponding to each non-missing water level value are obtained. The non-missing water level value and the corresponding precipitation and evapotranspiration constitute the first sample, obtain each first sample, and then enter step A2-2.

[0083] Step A2-2. Based on each first sample, with the precipitation and evapotranspiration in the first sample as input and the non-missing water level value in the first sample as output, train the LSTM neural network to obtain a water level value neural network model, and then enter step A2-3.

[0084] Step A2-3. For each missing water level value, perform the following steps:

[0085] According to the boundary information of the water area corresponding to the missing water level value at the corresponding historical time point, each grid under the precipitation grid distribution of the water area based on the boundary information is obtained, and the average value of the precipitation of each grid at the historical time point corresponding to the missing water level value is obtained, which constitutes the precipitation of the local area of ​​the target area where the water area corresponding to the missing water level value is located at the corresponding historical time point.

[0086] At the same time, according to the boundary information of the water area corresponding to the missing water level value at the corresponding historical time point, each grid under the grid distribution of the evapotranspiration involved in the water area based on the boundary information is obtained, and the average value of the evapotranspiration of each grid at the historical time point corresponding to the missing water level value is obtained, which constitutes the evapotranspiration of the local area of ​​the target area where the water area corresponding to the missing water level value is located at the corresponding historical time point.

[0087] That is, the precipitation and evapotranspiration corresponding to the missing water level value are obtained, and then the precipitation and evapotranspiration corresponding to each missing water level value are obtained, and then step A2-4 is entered.

[0088] Step A2-4. For each missing water level value, take the precipitation and evapotranspiration corresponding to the missing water level value as input, apply the water level value neural network model for processing, obtain the corresponding water level value, fill the missing water level value, and then fill each missing water level value.

[0089] Step A3. Perform spherical harmonic expansion on the grid distribution of soil moisture content, snow water equivalent, canopy water storage, precipitation, evapotranspiration, and terrestrial water storage at each preset historical time point within the preset historical time period corresponding to the target area, obtain the spherical harmonic model coefficients, and then adjust the spatial resolution of the grid distribution, and update it to the preset target spatial grid resolution. In the process of performing the spherical harmonic expansion, the sparse least squares method is used to minimize the residual between the original field and the simulated field, weaken the Gibbs effect in the spherical harmonic expansion, and then enter step A4.

[0090] Step A4. Determine whether there is a gap in the grid distribution of terrestrial water reserves at each preset historical time point within the preset historical time period corresponding to the target area. If so, execute the following steps A4-1 to A4-4 to fill in the missing terrestrial water reserves; otherwise, proceed to step B.

[0091] Step A4-1. For each non-missing terrestrial water storage grid value, a second sample is constructed by combining the non-missing terrestrial water storage grid value with the soil moisture grid value, precipitation grid value, and evapotranspiration grid value corresponding to the corresponding historical time point of the grid, obtaining each second sample, and then entering step A4-2.

[0092] Step A4-2. Based on each second sample, the soil moisture grid value, precipitation grid value, and evapotranspiration grid value in the second sample are used as input, and the non-missing land water storage grid value in the second sample is used as output. The LSTM neural network is trained to obtain a land water storage grid value neural network model, and then enter step A4-3.

[0093] Step A4-3. For each missing land water storage grid value, take the soil moisture grid value, precipitation grid value and evapotranspiration grid value of the grid where the missing land water storage grid value is located at the corresponding historical time point as input, apply the land water storage grid value neural network model to process, obtain the corresponding land water storage grid value, fill in the missing land water storage grid value, and then fill in each missing land water storage grid value.

[0094] Step B. Based on the boundary information of each water body area in the target area corresponding to each preset historical time point in the preset historical time period, according to the area of ​​each water body area in each grid corresponding to each preset historical time point under the target spatial grid resolution, combined with the water level value of each water body area corresponding to each preset historical time point, calculate and obtain the surface water reserves of each water body area in each grid corresponding to each preset historical time point under the target spatial grid resolution, and then define that the surface water reserves of other grids under the target spatial grid resolution corresponding to each preset historical time point are equal to 0, construct the grid distribution of surface water reserves of the target area at each preset historical time point in the preset historical time period corresponding to the target spatial grid resolution, and then enter step C.

[0095] Step C. Based on the ensemble Kalman filter, the estimated groundwater reserves of all grid points in the target area corresponding to the target spatial grid resolution at each preset historical time point in the preset historical time period are calculated to estimate the groundwater reserves in the target area.

[0096] In actual application, the above step C specifically performs the following steps C1 to C8.

[0097] Step C1. According to the following formula:

[0098]

[0099]

[0100] Define the state vector x of the target area corresponding to the historical time point t t , and the observation vector y of the target area corresponding to the historical time point t t , where 1,…,i,…,I, I represents the number of grids in the target area with respect to the target space grid resolution, It represents the underground water storage grid value calculated by combining the soil moisture grid value, snow water equivalent grid value, canopy water storage grid value, surface water storage grid value, and land water storage grid value at the target area corresponding to the historical time point t at the i-th grid of the target spatial grid resolution. represents the grid value of land water storage in the target area at the historical time point t corresponding to the i-th grid of the target spatial grid resolution, [·] T Indicates transposition, and then goes to step C2.

[0101] Step C2. Construct the state vector x t With the observation vector y t The correlation equation between them is as follows:

[0102] y t =H t x t +v t

[0103] Among them, H t Represents the corresponding historical time point t and state vector x t Coefficient matrix of the same dimension, v t Represents the corresponding historical time point t and state vector x t Observation noise of the same dimension, and then enter step C3.

[0104] Step C3. Based on the preset E different disturbances, the state vector x is calculated for the combined soil moisture content, snow water equivalent, canopy water storage, surface water storage, and terrestrial water storage. t , add each disturbance separately and construct the set X corresponding to the historical time point t t as follows:

[0105]

[0106] Among them, 1,…,e,…,E, Represents the state vector x t Add the initial amount of the e-th disturbance, and then go to step C4.

[0107] Step C4. For the state vector x t-1 Add the initial amount of each disturbance Execute time forward propagation to obtain each The corresponding updated values Then go to step C5.

[0108] Step C5. According to the following formula:

[0109]

[0110]

[0111] Construct the collection matrix X t|t-1 and the state covariance update value matrix P t|t-1 ,in, Represents the set matrix X t|t-1 Then go to step C6.

[0112] Step C6. According to the following formula:

[0113]

[0114] Calculate the gain matrix K corresponding to the historical time point t t , where R t Represents the covariance matrix of the observation vector corresponding to the historical time point t, and then enters step C7.

[0115] Step C7. According to the following formula:

[0116]

[0117] Get each Correction value in, Represents the observation vector y t Add the initial amount of the e-th disturbance, and then go to step C8.

[0118] Step C8. According to the following formula:

[0119]

[0120] Obtain the estimated value of groundwater storage corresponding to historical time point t That is, the estimated groundwater reserves of all grid points in the target area at each preset historical time point t within the preset historical time period corresponding to the target spatial grid resolution are used to estimate the groundwater reserves in the target area.

[0121] After the estimation of groundwater reserves in the target area is realized through the above design, further practical applications such as Figure 1 As shown, by executing the following steps D to F, the net recharge of underground water storage in the target area corresponding to the historical analysis time period L can be analyzed.

[0122] Step D. Estimated groundwater storage values ​​of grid points at each preset historical time point t in the preset historical time period according to the target area with respect to the target spatial grid resolution According to the following formula:

[0123]

[0124] Obtain the estimated groundwater storage fluctuation value of the target area at the grid point corresponding to the historical time point t compared with the previous historical time point t-1 at the target spatial grid resolution Among them, 1,…,i,…,I, I represents the number of grids in the target area with respect to the target space grid resolution, It represents the estimated value of underground water storage at the grid point of the target area at the historical time point t corresponding to the i-th grid of the target spatial grid resolution. It represents the estimated value of underground water storage at the grid point of the target area at the historical time point t-1 corresponding to the i-th grid of the target spatial grid resolution, and then enters step E.

[0125] Step E. Based on the historical analysis time period L, the following formula is used:

[0126]

[0127] Obtain the potential recharge R of groundwater storage in the target area corresponding to the historical analysis period LL ;

[0128] At the same time, according to the following formula:

[0129]

[0130] Obtain the potential discharge D of underground water storage in the target area corresponding to the historical analysis period L L ,|| means taking the absolute value;

[0131] Then go to step F.

[0132] Step F. According to the following formula:

[0133] R net,L =R L -D L

[0134] Obtain the net recharge R of groundwater storage in the target area corresponding to the historical analysis period L net,L .

[0135] In the design of the scheme, spherical harmonic functions were used to unify the spatial resolution of hydrological data, and the sparse least squares method was used to effectively weaken the Gibbs effect. The missing terrestrial water storage and water level data were reconstructed by establishing an LSTM neural network to ensure the integrity of the groundwater storage estimation results. In addition, the method of estimating groundwater reserves using ensemble Kalman filtering incorporates the uncertainties in the observed data and models, and no longer relies on the prior assumptions and decomposition processes of traditional estimation methods.

[0136] The above design scheme was applied in practice, and the data of each year in the preset historical period of 21 years corresponding to the preset target spatial grid resolution in the African region were analyzed, which overcame the limitation of missing water reserves in some time periods. Compared with the traditional decomposition method, the accuracy and reliability of ensemble Kalman filtering in groundwater reserve estimation have been significantly improved.

[0137] The embodiments of the present invention are described in detail above with reference to the accompanying drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A groundwater estimation method based on ensemble Kalman filtering to achieve the estimation of groundwater reserves in a target area, characterized in that: The steps include: Step A. Obtain the grid distribution of terrestrial water storage, soil moisture content, snow water equivalent, canopy water storage, precipitation, and evapotranspiration at preset historical time points in the preset historical time period corresponding to the preset target spatial grid resolution in the target area, and obtain the water level and boundary information of each water body area in the target area corresponding to each preset historical time point in the preset historical time period, and then proceed to step B; Step B. Based on the boundary information of each water body area in the target area corresponding to each preset historical time point in the preset historical time period, according to the area of ​​each water body area in each grid corresponding to each preset historical time point at the target spatial grid resolution, combined with the water level value of each water body area corresponding to each preset historical time point, calculate and obtain the surface water reserves of each water body area in each grid corresponding to each preset historical time point at the target spatial grid resolution, and then define that the surface water reserves of other grids at the target spatial grid resolution corresponding to each preset historical time point are equal to 0, construct the surface water reserves grid distribution of the target area at each preset historical time point in the preset historical time period corresponding to the target spatial grid resolution, and then enter step C; Step C. Based on the ensemble Kalman filter, the estimated groundwater reserves of all grid points in the target area corresponding to the target spatial grid resolution at each preset historical time point in the preset historical time period are calculated to estimate the groundwater reserves in the target area.

2. The groundwater estimation method based on ensemble Kalman filtering according to claim 1 is characterized in that: It also includes steps D to F as follows: after executing step C, proceed to step D; Step D. Estimated groundwater storage values ​​of grid points at each preset historical time point t in the preset historical time period according to the target area with respect to the target spatial grid resolution According to the following formula: Obtain the estimated groundwater storage fluctuation value of the target area at the grid point corresponding to the historical time point t compared with the previous historical time point t-1 at the target spatial grid resolution Among them, 1,…,i,…,I, I represents the number of grids in the target area with respect to the target space grid resolution, It represents the estimated value of underground water storage at the grid point of the target area at the historical time point t corresponding to the i-th grid of the target spatial grid resolution. Represents the estimated value of underground water storage at the grid point of the target area at the target spatial grid resolution corresponding to the historical time point t-1, and then enters step E; Step E. Based on the historical analysis time period L, the following formula is used: Obtain the potential recharge R of groundwater storage in the target area corresponding to the historical analysis period L L ; At the same time, according to the following formula: Obtain the potential discharge D of underground water storage in the target area corresponding to the historical analysis period L L ,|| means taking the absolute value; Then proceed to step F; Step F. According to the following formula: R net,L =R L -D L Obtain the net recharge R of groundwater storage in the target area corresponding to the historical analysis period L net,L .

3. The groundwater estimation method based on ensemble Kalman filtering according to claim 1 is characterized in that: The step A includes the following steps A1 to A4; Step A1. Obtain the grid distribution of terrestrial water storage, soil moisture content, snow water equivalent, canopy water storage, precipitation, and evapotranspiration at each preset historical time point in the preset historical time period corresponding to the target area, and obtain the water level and boundary information of each water body area in the target area at each preset historical time point in the preset historical time period, and then proceed to step A2; Step A2. Determine whether there are missing water level values ​​at each preset historical time point in each water body area in the target area corresponding to the preset historical time period, and if so, fill in the missing water level values; otherwise, proceed to step A3; Step A3. Update and unify the grid distribution of soil moisture content, snow water equivalent, canopy water storage, precipitation, evapotranspiration, and terrestrial water storage at each preset historical time point in the preset historical time period corresponding to the target area to the preset target spatial grid resolution, and then proceed to step A4; Step A4. Determine whether there is a gap in the grid distribution of land water reserves at each preset historical time point in the target area corresponding to the preset historical time period, and if so, fill in the gaps in the land water reserves; Otherwise go to step B.

4. The groundwater estimation method based on ensemble Kalman filtering according to claim 3 is characterized by: In step A2, when it is determined that there are missing water level values ​​corresponding to each preset historical time point in each water body area in the target area within the preset historical time period, the following steps A2-1 to A2-4 are executed to fill in the missing water level values; Step A2-1. For each non-missing water level value, perform the following steps: According to the boundary information of the water area corresponding to the non-missing water level value at the corresponding historical time point, each grid under the precipitation grid distribution involved in the water area based on the boundary information is obtained, and the average value of the precipitation of each grid at the historical time point corresponding to the non-missing water level value is obtained, so as to constitute the precipitation of the local area of ​​the target area where the water area corresponding to the non-missing water level value is located at the corresponding historical time point; At the same time, according to the boundary information of the water area corresponding to the non-missing water level value at the corresponding historical time point, each grid under the evapotranspiration grid distribution involved in the water area based on the boundary information is obtained, and the average value of the evapotranspiration of each grid at the historical time point corresponding to the non-missing water level value is obtained, so as to constitute the evapotranspiration of the local area of ​​the target area where the water area corresponding to the non-missing water level value is located at the corresponding historical time point; That is, the precipitation and evapotranspiration corresponding to the non-missing water level value are obtained, and then the precipitation and evapotranspiration corresponding to each non-missing water level value are obtained, and the non-missing water level value and the corresponding precipitation and evapotranspiration constitute a first sample, and each first sample is obtained, and then step A2-2 is entered; Step A2-2. Based on each first sample, with the precipitation and evapotranspiration in the first sample as input and the non-missing water level value in the first sample as output, train the LSTM neural network to obtain a water level value neural network model, and then proceed to step A2-3; Step A2-3. For each missing water level value, execute the following: According to the boundary information of the water area corresponding to the missing water level value at the corresponding historical time point, each grid under the precipitation grid distribution involved in the water area based on the boundary information is obtained, and the average value of the precipitation of each grid at the historical time point corresponding to the missing water level value is obtained, so as to constitute the precipitation of the local area of ​​the target area where the water area corresponding to the missing water level value is located at the corresponding historical time point; At the same time, according to the boundary information of the water area corresponding to the missing water level value at the corresponding historical time point, each grid under the evapotranspiration grid distribution involved in the water area based on the boundary information is obtained, and the average value of the evapotranspiration of each grid at the historical time point corresponding to the missing water level value is obtained, so as to constitute the evapotranspiration of the local area of ​​the target area where the water area corresponding to the missing water level value is located at the corresponding historical time point; That is, the precipitation and evapotranspiration corresponding to the missing water level value are obtained, and then the precipitation and evapotranspiration corresponding to each missing water level value are obtained, and then step A2-4 is entered; Step A2-4. For each missing water level value, take the precipitation and evapotranspiration corresponding to the missing water level value as input, apply the water level value neural network model for processing, obtain the corresponding water level value, fill the missing water level value, and then fill each missing water level value.

5. The groundwater estimation method based on ensemble Kalman filtering according to claim 3 is characterized by: In step A3, spherical harmonic expansion is performed on the soil moisture grid distribution, snow water equivalent grid distribution, canopy water storage grid distribution, precipitation grid distribution, evapotranspiration grid distribution, and terrestrial water storage grid distribution at each preset historical time point within the preset historical time period corresponding to the target area to obtain the spherical harmonic model coefficients, and then the spatial resolution of the grid distribution is adjusted to be updated and unified to the preset target spatial grid resolution.

6. The groundwater estimation method based on ensemble Kalman filtering according to claim 5 is characterized by: In step A3, during the spherical harmonic expansion, a sparse least squares method is used to minimize the residual between the original field and the simulated field, thereby weakening the Gibbs effect in the spherical harmonic expansion.

7. The groundwater estimation method based on ensemble Kalman filtering according to claim 3 is characterized by: In the step A4, when it is determined that there are missing grid values ​​of terrestrial water storage at each preset historical time point within the preset historical time period corresponding to the target area, the following steps A4-1 to A4-4 are executed to fill in the missing grid values ​​of terrestrial water storage; step A4-1. For each non-missing grid value of terrestrial water storage, a second sample is formed by combining the non-missing grid value of terrestrial water storage with the grid value of soil moisture, the grid value of precipitation, and the grid value of evapotranspiration corresponding to the corresponding historical time point of the grid where the grid is located, and each second sample is obtained, and then the step A4-2 is entered; Step A4-2. Based on each second sample, the soil moisture grid value, precipitation grid value, and evapotranspiration grid value in the second sample are used as input, and the non-missing land water storage grid value in the second sample is used as output, and the LSTM neural network is trained to obtain a land water storage grid value neural network model, and then proceed to step A4-3; Step A4-3. For each missing land water storage grid value, take the soil moisture grid value, precipitation grid value and evapotranspiration grid value of the grid where the missing land water storage grid value is located at the corresponding historical time point as input, apply the land water storage grid value neural network model to process, obtain the corresponding land water storage grid value, fill in the missing land water storage grid value, and then fill in each missing land water storage grid value.

8. The groundwater estimation method based on ensemble Kalman filtering according to claim 1, characterized in that: The step C includes the following steps C1 to C8; Step C1. According to the following formula: Define the state vector x of the target area corresponding to the historical time point t t , and the observation vector y of the target area corresponding to the historical time point t t , where 1,…,i,…,I, I represents the number of grids in the target area with respect to the target space grid resolution, It represents the underground water storage grid value calculated by combining the soil moisture grid value, snow water equivalent grid value, canopy water storage grid value, surface water storage grid value, and land water storage grid value at the target area corresponding to the historical time point t at the i-th grid of the target spatial grid resolution. represents the grid value of land water storage in the target area at the historical time point t corresponding to the i-th grid of the target spatial grid resolution, [·] T Indicates transposition, and then goes to step C2; Step C2. Construct the state vector x t With the observation vector y t The correlation equation between them is as follows: y t =H t x t +v t Among them, H t Represents the corresponding historical time point t and state vector x t Coefficient matrix of the same dimension, v t Represents the corresponding historical time point t and state vector x t Observation noise of the same dimension, and then enter step C3; Step C3. Based on the preset E different disturbances, the state vector x is calculated for the combined soil moisture content, snow water equivalent, canopy water storage, surface water storage, and terrestrial water storage. t , add each disturbance separately and construct the set X corresponding to the historical time point t t as follows: Among them, 1,…,e,…,E, Represents the state vector x t Add the initial amount of the e-th disturbance, and then go to step C4; Step C4. For the state vector x t-1 Add the initial amount of each disturbance Execute time forward propagation to obtain each The corresponding updated values Then proceed to step C5; Step C5. According to the following formula: Construct the collection matrix X t|t-1 and the state covariance update value matrix P t|t-1 ,in, Represents the set matrix X t|t-1 The average value of each element in , and then go to step C6; Step C6. According to the following formula: Calculate the gain matrix K corresponding to the historical time point t t , where R t represents the covariance matrix of the observation vector corresponding to the historical time point t, and then proceeds to step C7; Step C7. According to the following formula: Get each Correction value in, Represents the observation vector y t Add the initial amount of the e-th disturbance, and then enter step C8; Step C8. According to the following formula: Obtain the estimated value of groundwater storage corresponding to historical time point t That is, the estimated groundwater reserves of all grid points in the target area at each preset historical time point t within the preset historical time period corresponding to the target spatial grid resolution are used to realize the estimation of groundwater reserves in the target area.

Images