A regional groundwater level monitoring method and device integrating multi-source data

Abstract

The present invention relates to a regional groundwater level monitoring method and device that fuses multi-source data. The method comprises the following steps: 1) determining a study area, downloading remote sensing images of the study area and model evapotranspiration products from the Global Land Data Assimilation System, and collecting pan evaporation and groundwater level observation data for the study area; 2) estimating surface evapotranspiration based on the visible and thermal infrared bands of the remote sensing images; 3) performing spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into daily-scale fused evapotranspiration; 4) establishing a surface evapotranspiration and groundwater level linkage model based on pan evaporation, groundwater level observation data, and daily-scale fused evapotranspiration at corresponding locations; and 5) inputting the daily-scale fused evapotranspiration from step 3 into the surface evapotranspiration and groundwater level linkage model obtained from step 4 to estimate the regional groundwater level. The present invention strengthens the physical foundation of groundwater level monitoring through remote sensing evapotranspiration, making it possible to estimate groundwater levels in areas without or lacking data.

Description

Technical Field

[0001] The present invention relates to a regional groundwater level monitoring method, and in particular to a regional groundwater level monitoring method and device that integrates multi-source data. Background Art

[0002] Groundwater, an important component of water resources, is the main source of regional agricultural irrigation, industrial water consumption and domestic water needs. It is the main water source for the vegetation-soil-atmosphere system and an important carrier of mass and energy cycles. It directly affects the regional climate characteristics and vegetation growth and type distribution. Monitoring it is a key guarantee for regional sustainable development.

[0003] Traditional manual field monitoring of groundwater levels is time-consuming and labor-intensive, with sparse observation points and poor spatial representation, making it impossible to dynamically monitor and evaluate groundwater levels over large areas. Remote sensing technology, as an effective means of macroscopic, comprehensive, dynamic, and rapid monitoring and evaluation of natural resources, has become a growing trend in estimating groundwater depth over large areas. Remote sensing groundwater level monitoring primarily utilizes indicators closely related to groundwater levels, such as soil moisture, surface temperature, and vegetation index.

[0004] The principle of groundwater level estimation based on soil moisture is that the amount of capillary water that replenishes surface soil moisture varies with the groundwater level. Therefore, changes in soil moisture at a given depth reflect the groundwater level. The physical basis for estimating soil moisture using remote sensing is that changes in soil moisture affect the spectral characteristics of the surface and vegetation. A remote sensing inversion model for soil moisture is established by linearly and nonlinearly fitting remotely sensed parameters with soil moisture observations. Using groundwater level data, experimental equations are then established that integrate soil moisture and other auxiliary data. Groundwater level distribution is then estimated using remotely sensed soil moisture. In arid regions, soil moisture near the surface is very low, and intense evaporation minimizes soil moisture changes. Therefore, groundwater level estimation based on soil moisture is unsuitable for these regions.

[0005] The principle of groundwater level estimation based on surface temperature is that groundwater enrichment and water level directly affect surface temperature due to the difference in specific heat between water and soil. Groundwater level is indirectly estimated using surface temperature inverted using remote sensing technology and its derived indices. These indices primarily include the Temperature Vegetation Drought Index, the Crop Water Stress Index, the Vegetation Water Supply Index, the Water Deficit Index, and the Apparent Thermal Inertia. The relationship between surface temperature and its derived indices and groundwater level is susceptible to numerous meteorological factors, including the atmosphere-surface energy cycle and surface physical processes, which reduces the accuracy of groundwater level estimation.

[0006] The principle behind estimating groundwater levels based on vegetation indices is that plants require an appropriate underground ecological water level at all stages to maintain their natural growth. Excessively low groundwater levels can lead to plant wilting due to insufficient water demand, while excessively high water levels can cause soil salinization, which can affect plant growth. Therefore, vegetation growth is an indicator of groundwater levels. Vegetation indices calculated using remote sensing technology can indirectly reflect groundwater levels. These indices, such as the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EDVI), Ratio Vegetation Index (RI), and Soil Adjusted Vegetation Index (SAVI), can be combined with groundwater level observations and mathematical modeling to monitor regional groundwater levels. These indices are more suitable for areas with high vegetation cover but are less applicable to arid regions with sparse vegetation cover.

[0007] To overcome the problems of poor universality and weak physical foundation of current remote sensing methods for estimating regional groundwater levels, a regional groundwater level monitoring method was invented by combining remote sensing surface evapotranspiration inversion results with evaporation pan evaporation and groundwater level observation. Summary of the Invention

[0008] In view of the shortcomings of the existing technology, the present invention provides a regional groundwater level monitoring method and device that integrates multi-source data.

[0009] The present invention solves the above-mentioned technical problem with the following technical solution: A method for regional groundwater level monitoring by integrating multi-source data, comprising the following steps:

[0010] Step 1: Determine the study area, download the visible light band and thermal infrared band of remote sensing images of the study area and the model evapotranspiration product of the Global Land Data Assimilation System (GLDAS), crop and project them, and collect evaporation pan and groundwater level observation data of the study area;

[0011] Step 2: Estimate surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of remote sensing images;

[0012] Step 3: Perform spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into fused evapotranspiration at the daily scale;

[0013] Step 4: Build a model linking surface evaporation and groundwater level based on pan evaporation, groundwater level observations, and fused evapotranspiration at the daily scale.

[0014] Step 5: Input the fused evapotranspiration into the surface evapotranspiration and groundwater level linkage model to estimate the regional groundwater level.

[0015] In the method for regional groundwater level monitoring by fusing multi-source data, the specific steps of estimating surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of the remote sensing image in step 2 are as follows:

[0016] Step 2011: Remove cloud-contaminated pixel data using built-in quality control tags in remote sensing images; calibrate the remote sensing data using the gain and offset of the remote sensing data;

[0017] Step 2012: Calculate the Normalized Difference Vegetation Index (NDVI) (dimensionless) using the surface reflectance NIR in the near infrared band and the surface reflectance R in the red band of the remote sensing image data obtained in step 2011. The NDVI calculation formula is as follows:

[0018]

[0019] Step 2013: Determine the maximum NDVI value of the study area based on the soil moisture characteristics and vegetation cover characteristics of the study area max and minimum NDVI value NDVI min ;

[0020] Step 2014: NDVI obtained in step 2012 and NDVI obtained in step 2013 max and NDVI min Calculate the vegetation coverage FVC (dimensionless) of the study area. The FVC calculation formula is as follows:

[0021]

[0022] Step 2015: Determine the emissivity ε of the bare soil coverage of the study area based on the soil type characteristics and vegetation coverage characteristics of the study area soil (dimensionless) and the emissivity ε of full vegetation cover veg (dimensionless);

[0023] Step 2016: Based on the FVC obtained in step 2014 and the ε determined in step 2015 soil and ε veg Calculate the surface emissivity EM (dimensionless) using the following formula:

[0024] EM=FVC*ε veg +(1-FVC)*ε soil (3)

[0025] Step 2017: Based on the brightness temperature data T of the EM and remote sensing images obtained in step 2016 B Calculate the true surface temperature T s (K), calculated as follows:

[0026]

[0027] Step 2018: The blue light band reflectivity ρ of the remote sensing image obtained in step 2011 B(dimensionless), red light band reflectivity ρ R (dimensionless), NIR band reflectivity ρ NIR (dimensionless), mid-wave infrared reflectivity ρ swir (dimensionless), mid-wave infrared reflectivity ρ swir2 (dimensionless) Calculate the albedo α (dimensionless) using the following formula:

[0028] α=(0.356*ρ B +0.13*ρ R +0.373*ρ NIR +0.085*ρ swir +0.072*ρ swir2 -0.018) / 1.016 (5)

[0029] Step 2019: Based on the EM (dimensionless) obtained in step 2016 and the T obtained in step 2017 s (K), α (dimensionless) obtained in step 2018, combined with the solar downward shortwave radiation S d (W / m 2 ), Stefan-Boltzmann constant σ (dimensionless), atmospheric incident longwave radiation (W / m 2 ), calculate the surface net radiation flux R n (W / m 2 ), the calculation formula is as follows:

[0030]

[0031] Among them, the solar downward shortwave radiation S d The calculation formula is as follows:

[0032] S d =S0*cosZ*cosZ / (1.085*cosZ+ea*(2.7+cosZ)*0.001+b) (7)

[0033] Among them, S0 is the solar radiation constant at the top of the atmosphere, which is 1367W / m 2 , Z is the solar zenith angle obtained from the image data, in radians, ea is the actual atmospheric water vapor pressure, obtained from meteorological observations, in hPa, and b = 0.2 is the correction factor;

[0034] Where α is the albedo, dimensionless, σ is the Stefan-Boltzmann constant, which is 0.000000056697, EM is the emissivity of the surface, dimensionless, and T s is the surface temperature in K, ε sky is the sky emissivity, which can be approximately equal to 1, T skyis the sky temperature in K;

[0035] Step 2020: Based on the FVC obtained in step 2014 and the net surface radiation flux R obtained in step 2019 n , when the soil is completely bare, the surface heat flux G accounts for R n The ratio of G to R is 0.315. n The ratio is 0.05, and the surface heat flux G (W / m 2 ), the calculation formula is as follows:

[0036] G=R n [0.05+0.31*(1-FVC)] (8)

[0037] Step 2021: Based on the FVC obtained in step 2014 and the T obtained in step 2017 s , construct a model with FVC as the horizontal axis and T s The two-dimensional scattered point space of the vertical axis is divided into several equal parts, and the maximum T in each FVC partition is counted. s Tmaxv and minimum T s Tminv, with several equally divided FVC values as independent variables and the corresponding Tmaxv and Tminv as dependent variables, the following regression equation is obtained:

[0038] Tmaxv=ar+br*FVC (9)

[0039] Tminv=aw+bw*FVC (10)

[0040] Where ar and br are the slope and intercept of the dry side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively; aw and bw are the slope and intercept of the wet side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively.

[0041] Step 2022: Based on the FVC obtained in step 2014 and the T obtained in step 2017 s , Tmaxv and Tminv obtained in step 2021 are used to calculate the Bowen ratio β (dimensionless) of each pixel. The calculation formula is as follows:

[0042] β=(T s -Tminv) / (Tmaxv-T s ) (11)

[0043] Step 2023: R obtained according to step 2019 n , G obtained in step 2020, β obtained in step 2022, calculate the surface evapotranspiration ET (W / m 2 ), the calculation formula is as follows:

[0044] ET=(R n -G) / (1+β) (12)

[0045] In the method for regional groundwater level monitoring by integrating multi-source data, the specific steps of step 3 are as follows:

[0046] Step 3011: Prepare the corresponding remote sensing evapotranspiration ET (W / m 2 ) GLDAS model evapotranspiration ET r (kg / m 2 / s), corresponding to the GLDAS model evaporation on the observation day of evaporation pan and groundwater level {ET Oi (kg / m 2 / s)} 1≤i≤8 ;

[0047] Step 3012: Convert remote sensing evapotranspiration ET (W / m 2 ) is converted to kg / m 2 / s;

[0048] Step 3013: Convert the dimension of remote sensing evapotranspiration ET (kg / m 2 / s), model evapotranspiration ET r (kg / m 2 / s), model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 Perform spatiotemporal fusion to obtain the corresponding evaporation of the evaporation pan, the observation day of the groundwater level, and the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 , the fusion method is:

[0049]

[0050] Where W is the size of the small window moving on the remote sensing evapotranspiration and model evapotranspiration, k, l are the positions of the pixels in the small window, and w kl is the weight at pixel (k, l), and satisfies 0≤w≤1,∑w=1;

[0051] Step 3014: Combine the evaporation Ri (kg / m 2 / s)} 1≤i≤8 Dimension conversion to integrated evapotranspiration ET R (mm / d);

[0052] In the method for regional groundwater level monitoring by integrating multi-source data, the specific steps of step 4 are as follows:

[0053] Step 4011: Observe the evaporation E00 (mm / d) using the evaporation dish and use the inverse S-shaped curve to establish the fused evapotranspiration ET R The relationship model between (mm / d) and groundwater level H (m) is as follows:

[0054]

[0055] Step 4012: Convert equation (14) into the following form:

[0056]

[0057] Step 4013: Using the observation data of E00 and H and the fused evapotranspiration ET R Perform a linear fit on equation (15), and the fitting steps are as follows:

[0058] a) Using E00 and ET R The dependent variable is calculated based on the data Calculate the independent variable x=lnH and use least squares regression to establish the following fitting equation:

[0059] y=kx+m (16)

[0060] b) Using the fitting results from step a), calculate the parameters c and d in equation (14):

[0061]

[0062] c) Substitute d and c obtained in step b) into equation (14) to obtain the relationship model between surface evaporation and groundwater level:

[0063]

[0064] The method for regional groundwater level monitoring by fusing multi-source data, wherein in step 5, the fused evapotranspiration ET obtained in step 3 is used as the R Substitute the relationship model between surface evapotranspiration and groundwater level obtained in step 4 to calculate the groundwater level. The calculation formula is as follows:

[0065]

[0066] The technical solution of the present invention to solve the above technical problems is: a device for regional groundwater level monitoring method integrating multi-source data, including a data preparation module, a remote sensing evapotranspiration estimation module, a remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module, a surface evapotranspiration and groundwater level modeling module, and a groundwater level estimation module;

[0067] The data preparation module is used to determine the study area, download the visible light band and thermal infrared band of remote sensing images covering the study area, and the GLDAS model evapotranspiration product, and perform cropping and projection conversion, and collect evaporation pan and groundwater level observation data in the study area;

[0068] The remote sensing evapotranspiration estimation module is used to estimate surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of the remote sensing image;

[0069] The remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module is used to perform spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into fused evapotranspiration at a day scale;

[0070] The surface evapotranspiration and groundwater level modeling module is used to establish a connection model between surface evapotranspiration and groundwater level based on the observation data of evaporation pan and groundwater level and the fused evapotranspiration at the day scale;

[0071] The groundwater level estimation module is used to input the fused evapotranspiration into the connection model of surface evapotranspiration and groundwater level to estimate the regional groundwater level;

[0072] The device for regional groundwater level monitoring method integrating multi-source data, wherein the remote sensing evapotranspiration estimation module includes a vegetation index calculation unit, a vegetation coverage calculation unit, an emissivity calculation unit, a surface temperature calculation unit, an albedo calculation unit, a surface net radiation flux calculation unit, a surface heat flux calculation unit, a surface temperature and vegetation coverage two-dimensional scattered point space construction unit, a Bowen ratio calculation unit, and a surface evapotranspiration remote sensing estimation unit;

[0073] The vegetation index calculation unit is used to remove cloud-contaminated pixel data through built-in quality control tags of remote sensing images; and calibrate the remote sensing data through the gain and offset of the remote sensing data;

[0074] The normalized difference vegetation index (NDVI) (dimensionless) is calculated using the surface reflectance NIR of the near infrared band and the surface reflectance R of the red light band of the remote sensing image data. The NDVI calculation formula is as follows:

[0075]

[0076] The vegetation coverage calculation unit is used to determine the maximum NDVI value NDVI of the study area based on the soil moisture characteristics and vegetation coverage characteristics of the study area. max and minimum NDVI value NDVI min ;

[0077] NDVI obtained by calculating the vegetation index unit and the determined NDVI max and NDVI minCalculate the vegetation coverage FVC (dimensionless) of the study area. The FVC calculation formula is as follows:

[0078]

[0079] The emissivity calculation unit is used to determine the emissivity ε of the bare soil full coverage of the study area according to the soil type characteristics and vegetation coverage characteristics of the study area. soil (dimensionless) and the emissivity ε of full vegetation cover veg (dimensionless);

[0080] According to the vegetation coverage calculation unit obtained FVC and the determined ε soil and ε veg Calculate the surface emissivity EM (dimensionless) using the following formula:

[0081] EM=FVC*ε veg +(1-FVC)*ε soil (twenty two)

[0082] The surface temperature calculation unit is used to calculate the brightness temperature data T of the EM and remote sensing image obtained by the emissivity calculation unit. B Calculate the true surface temperature T s (K), calculated as follows:

[0083]

[0084] The albedo calculation unit is used to calculate the blue light band reflectivity ρ of the remote sensing image obtained by the vegetation index calculation unit. B (dimensionless), red light band reflectivity ρ R (dimensionless), NIR band reflectivity ρ NIR (dimensionless), mid-wave infrared reflectivity ρ swir (dimensionless), mid-wave infrared reflectivity ρ swir2 (dimensionless) Calculate the albedo α (dimensionless) using the following formula:

[0085] α=(0.356*ρ B +0.13*ρ R +0.373*ρ NIR +0.085*ρ swir +0.072*ρ swir2 -0.018) / 1.016 (24)

[0086] The surface net radiation flux calculation unit is used to calculate the EM (dimensionless) obtained by the emissivity calculation unit and the T obtained by the surface temperature calculation unit. s(K), α (dimensionless) obtained by the albedo calculation unit, combined with the solar downward shortwave radiation S d (W / m 2 ), Stefan-Boltzmann constant σ (dimensionless), atmospheric incident longwave radiation (W / m 2 ), calculate the surface net radiation flux R n (W / m 2 ), the calculation formula is as follows:

[0087]

[0088] Among them, the solar downward shortwave radiation S d The calculation formula is as follows:

[0089] S d =S0*cosZ*cosZ / (1.085*cosZ+ea*(2.7+cosZ)*0.001+b) (26)

[0090] Among them, S0 is the solar radiation constant at the top of the atmosphere, which is 1367W / m 2 , Z is the solar zenith angle obtained from the image data, in radians, ea is the actual atmospheric water vapor pressure, obtained from meteorological observations, in hPa, and b = 0.2 is the correction factor;

[0091] Where α is the albedo, dimensionless, the Stefan-Boltzmann constant σ is 0.000000056697, EM is the emissivity of the surface, dimensionless, Ts is the surface temperature in K, ε sky is the sky emissivity, which can be approximately equal to 1, T sky is the sky temperature in K;

[0092] The surface heat flux calculation unit is used to calculate the surface net radiation flux R obtained by the vegetation coverage calculation unit and the surface heat flux calculation unit. n , when the soil is completely bare, the surface heat flux G accounts for R n The ratio of G to R is 0.315. n The ratio is 0.05, and the surface heat flux G (W / m 2 ), the calculation formula is as follows:

[0093] G=R n [0.05+0.31*(1-FVC)] (27)

[0094] The surface temperature and vegetation coverage two-dimensional scattered point space construction unit is used to calculate the FVC obtained by the vegetation coverage calculation unit and the T obtained by the surface temperature calculation unit. s, construct a model with FVC as the horizontal axis and T s The two-dimensional scattered point space of the vertical axis is divided into several equal parts, and the maximum T in each FVC partition is counted. s Tmaxv and minimum T s Tminv, with several equally divided FVC values as independent variables and the corresponding Tmaxv and Tminv as dependent variables, the following regression equation is obtained:

[0095] Tmaxv=ar+br*FVC (28)

[0096] Tminv=aw+bw*FVC (29)

[0097] Where ar and br are the slope and intercept of the dry side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively; aw and bw are the slope and intercept of the wet side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively.

[0098] The Bowen ratio calculation unit is used to calculate the FVC obtained by the vegetation coverage calculation unit and the T obtained by the surface temperature calculation unit. s The Tmaxv and Tminv obtained by constructing the unit in the two-dimensional scattered space of surface temperature and vegetation cover are used to calculate the Bowen ratio β (dimensionless) of each pixel. The calculation formula is as follows:

[0099] β=(T s -Tminv) / (Tmaxv-T s ) (30)

[0100] The surface evapotranspiration remote sensing estimation unit is used to calculate the R n , G obtained by the surface heat flux calculation unit, β obtained by the Bowen ratio calculation unit, and the surface evapotranspiration ET (W / m 2 ), the calculation formula is as follows:

[0101] ET=(R n -G) / (1+β) (31)

[0102] The device of the regional groundwater level monitoring method of fusing multi-source data, wherein the remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module includes a fusion data preparation unit, a remote sensing evapotranspiration dimension conversion unit, a spatiotemporal fusion unit, and a fusion evapotranspiration dimension conversion unit;

[0103] The fusion data preparation unit is used to prepare the corresponding remote sensing evapotranspiration ET (W / m 2 ) GLDAS model evapotranspiration ET r (kg / m 2 / s), corresponding to the GLDAS model evaporation on the observation day of evaporation pan and groundwater level {ET Oi (kg / m 2 / s)} 1≤i≤8 ;

[0104] The remote sensing evapotranspiration dimension conversion unit is used to convert the remote sensing evapotranspiration ET (W / m 2 ) is converted to kg / m 2 / s;

[0105] The space-time fusion unit is used to convert the dimension-converted remote sensing evapotranspiration ET (kg / m 2 / s), model evapotranspiration ET r (kg / m 2 / s), model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 Perform spatiotemporal fusion to obtain the corresponding evaporation of the evaporation pan, the observation day of the groundwater level, and the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 , the fusion method is:

[0106]

[0107] Where W is the size of the small window moving on the remote sensing evapotranspiration and model evapotranspiration, k, l are the positions of the pixels in the small window, and w kl is the weight at pixel (k, l), and satisfies 0≤w≤1,∑w=1;

[0108] The fused evapotranspiration dimension conversion unit is used to convert the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 Dimension conversion to integrated evapotranspiration ET R (mm / d).

[0109] The device for regional groundwater level monitoring method integrating multi-source data, wherein the surface evaporation and groundwater level modeling module includes an evaporation pan, groundwater level observation data and fusion evaporation conversion unit, a conversion data regression unit, and a surface evaporation and groundwater level modeling unit;

[0110] The evaporation pan evaporation, groundwater level observation data and fusion evaporation conversion unit are used to use the evaporation pan evaporation observation E00 (mm / d) to establish the fusion evaporation ET using the inverse S-shaped curve. R The relationship model between (mm / d) and groundwater level H (m) is as follows:

[0111]

[0112] The transformed data regression unit is used to transform equation (33) into the following form:

[0113]

[0114] The surface evapotranspiration and groundwater level modeling unit is used to use the observation data of E00 and H and the fusion evapotranspiration ET R Perform a linear fit on equation (34), and the fitting steps are as follows:

[0115] a) Using E00 and ET R The dependent variable is calculated based on the data Calculate the independent variable x=lnH and use least squares regression to establish the following fitting equation:

[0116] y=kx+m (35)

[0117] b) Using the fitting results from step a), calculate the parameters c and d in equation (33):

[0118]

[0119] c) Substitute d and c obtained in step b) into equation (33) to obtain the relationship model between surface evapotranspiration and groundwater level:

[0120]

[0121] The device of the regional groundwater level monitoring method of fusing multi-source data, wherein the groundwater level estimation module is used to combine the fused evapotranspiration ET obtained by the remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module R Substitute the relationship model between surface evapotranspiration and groundwater level obtained in the surface evapotranspiration and groundwater level modeling module to calculate the groundwater level. The calculation formula is as follows:

[0122]

[0123] The beneficial effects of the present invention are as follows: compared with the traditional remote sensing index groundwater level estimation method, the present invention strengthens the physical basis of groundwater level monitoring through remote sensing evaporation, reduces the impact of the complexity of land and air interactions on groundwater level estimation, and makes it possible to estimate groundwater levels in data-free areas and data-deficient areas. It has broad application prospects in research fields such as ecology, hydrology, and climate change. BRIEF DESCRIPTION OF THE DRAWINGS

[0124] Figure 1 This is a flow chart of a regional groundwater level monitoring method integrating multi-source data according to the present invention;

[0125] Figure 2 This is a flowchart for the specific implementation of step 1 of the present invention;

[0126] Figure 3 This is a flowchart for the specific implementation of step 2 of the present invention;

[0127] Figure 4 This is a flowchart for the specific implementation of step 3 of the present invention;

[0128] Figure 5 This is a flowchart for the specific implementation of step 4 of the present invention;

[0129] Figure 6 This is a device block diagram of a method for regional groundwater level monitoring that integrates multi-source data according to the present invention;

[0130] Figure 7 This is a structural block diagram of the remote sensing evapotranspiration estimation module of the present invention;

[0131] Figure 8 This is a structural block diagram of the spatiotemporal fusion module of remote sensing evapotranspiration and model evapotranspiration according to the present invention;

[0132] Figure 9 This is a structural block diagram of the surface evapotranspiration and groundwater level modeling module of the present invention;

[0133] In the accompanying drawings, the components represented by the reference numerals are as follows:

[0134] 1. Data preparation module, 2. Remote sensing evapotranspiration estimation module, 3. Remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module, 4. Surface evapotranspiration and groundwater level modeling module, 5. Groundwater level estimation module; 201. Vegetation index calculation unit, 202. Vegetation coverage calculation unit, 203. Emissivity calculation unit, 204. Surface temperature calculation unit, 205. Albedo calculation unit, 206. Surface net radiation flux calculation unit, 207. Surface heat flux calculation unit, 208. 8. Two-dimensional scattered point space construction unit for surface temperature and vegetation cover, 209. Bowen ratio calculation unit, 210. Surface evapotranspiration remote sensing estimation unit; 301. Fusion data preparation unit, 302. Remote sensing evapotranspiration dimension conversion unit, 303. Spatiotemporal fusion unit, 304. Fusion evapotranspiration dimension conversion unit; 401. Pan evaporation, groundwater level observation data and fused evapotranspiration conversion unit, 402. Conversion data regression unit, 403. Surface evaporation and groundwater level modeling unit. DETAILED DESCRIPTION

[0135] The present invention is described in detail below with reference to the accompanying drawings.

[0136] like Figure 1 As shown, a regional groundwater level monitoring method integrating multi-source data includes the following steps:

[0137] Step 1: Determine the study area, download the visible light band and thermal infrared band of remote sensing images of the study area and the model evapotranspiration product of the Global Land Data Assimilation System (GLDAS), crop and project them, and collect evaporation pan and groundwater level observation data of the study area;

[0138] Step 2: Estimate surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of remote sensing images;

[0139] Step 3: Perform spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into fused evapotranspiration at the daily scale;

[0140] Step 4: Build a model linking surface evaporation and groundwater level based on pan evaporation, groundwater level observations, and fused evapotranspiration at the daily scale.

[0141] Step 5: Input the fused evapotranspiration into the surface evapotranspiration and groundwater level linkage model to estimate the regional groundwater level.

[0142] like Figure 2 As shown in Figure 1, in step 1, the required relevant data include the visible light band and thermal infrared band of the remote sensing image of the study area, the model evapotranspiration product of GLDAS, evaporation pan evaporation and groundwater level observation data.

[0143] Remote sensing imagery of the study area in visible and thermal infrared bands can be downloaded from the Geospatial Data Cloud (https: / / www.gscloud.cn / search). Observational data for pan evaporation, vapor pressure, and air temperature in the study area can be downloaded from the National Meteorological Administration (http: / / data.cma.cn / site / index.html). Cloud-affected pixels must be removed from the remote sensing data using quality control bands.

[0144] like Figure 3 As shown, the specific steps of step 2 are:

[0145] Step 2011: Remove cloud-contaminated pixel data using built-in quality control tags in remote sensing images; calibrate the remote sensing data using the gain and offset of the remote sensing data;

[0146] Step 2012: Calculate the Normalized Difference Vegetation Index (NDVI) (dimensionless) using the surface reflectance NIR in the near infrared band and the surface reflectance R in the red band of the remote sensing image data obtained in step 2011. The NDVI calculation formula is as follows:

[0147]

[0148] Step 2013: Determine the maximum NDVI value of the study area based on the soil moisture characteristics and vegetation cover characteristics of the study area max and minimum NDVI value NDVI min ;

[0149] Step 2014: Based on the NDVI obtained in step 2012 and the NDVI determined in step 2013 max and NDVI min Calculate the vegetation coverage FVC (dimensionless) of the study area. The FVC calculation formula is as follows:

[0150]

[0151] Step 2015: Determine the emissivity ε of the bare soil coverage of the study area based on the soil type characteristics and vegetation coverage characteristics of the study area soil (dimensionless) and the emissivity ε of full vegetation cover veg (dimensionless);

[0152] Step 2016: Based on the FVC obtained in step 2014 and the ε obtained in step 2015 soil and ε veg Calculate the surface emissivity EM (dimensionless) using the following formula:

[0153] EM=FVC*ε veg +(1-FVC)*ε soil (41)

[0154] Step 2017: Based on the brightness temperature data T of the EM and remote sensing images obtained in step 2016 B Calculate the true surface temperature T s (K), calculated as follows:

[0155]

[0156] Step 2018: The blue light band reflectivity ρ of the remote sensing image obtained in step 2011 B (dimensionless), red light band reflectivity ρ R (dimensionless), NIR band reflectivity ρ NIR (dimensionless), mid-wave infrared reflectivity ρ swir (dimensionless), mid-wave infrared reflectivity ρ swir2 (dimensionless) Calculate the albedo α (dimensionless) using the following formula:

[0157] α=(0.356*ρ B +0.13*ρ R +0.373*ρ NIR +0.085*ρswir +0.072*ρ swir2 -0.018) / 1.016(43)

[0158] Step 2019: Based on the EM (dimensionless) obtained in step 2016 and the T obtained in step 2017 s (K), α (dimensionless) obtained in step 2018, combined with the solar downward shortwave radiation S d (W / m 2 ), Stefan-Boltzmann constant σ (dimensionless), atmospheric incident longwave radiation (W / m 2 ), calculate the surface net radiation flux R n (W / m 2 ), the calculation formula is as follows:

[0159]

[0160] Among them, the solar downward shortwave radiation S d The calculation formula is as follows:

[0161] S d =S0*cosZ*cosZ / (1.085*cosZ+ea*(2.7+cosZ)*0.001+b) (45)

[0162] Among them, S0 is the solar radiation constant at the top of the atmosphere, which is 1367W / m 2 , Z is the solar zenith angle obtained from the image data, in radians, ea is the actual atmospheric water vapor pressure, obtained from meteorological observations, in hPa, and b = 0.2 is the correction factor;

[0163] Where α is the albedo, dimensionless, σ is the Stefan-Boltzmann constant, which is 0.000000056697, EM is the emissivity of the surface, dimensionless, and T s is the surface temperature in K, ε sky is the sky emissivity, which can be approximately equal to 1, T sky is the sky temperature in K;

[0164] Step 2020: Based on the FVC obtained in step 2014 and the net surface radiation flux R obtained in step 2019 n , when the soil is completely bare, the surface heat flux G accounts for R n The ratio of G to R is 0.315. n The ratio is 0.05, and the surface heat flux G (W / m 2 ), the calculation formula is as follows:

[0165] G=R n[0.05+0.31*(1-FVC)] (46)

[0166] Step 2021: Based on the FVC obtained in step 2014 and the T obtained in step 2017 s , construct a model with FVC as the horizontal axis and T s The two-dimensional scattered point space of the vertical axis is divided into several equal parts, and the maximum T in each FVC partition is counted. s Tmaxv and minimum T s Tminv, with several equally divided FVC values as independent variables and the corresponding Tmaxv and Tminv as dependent variables, the following regression equation is obtained:

[0167] Tmaxv=ar+br*FVC (47)

[0168] Tminv=aw+bw*FVC (48)

[0169] Where ar and br are the slope and intercept of the dry side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively; aw and bw are the slope and intercept of the wet side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively.

[0170] Step 2022: Based on the FVC obtained in step 2014 and the T obtained in step 2017 s , Tmaxv and Tminv obtained in step 2021 are used to calculate the Bowen ratio β (dimensionless) of each pixel. The calculation formula is as follows:

[0171] β=(T s -Tminv) / (Tmaxv-T s ) (49)

[0172] Step 2023: R obtained according to step 2019 n , G obtained in step 2020, β obtained in step 2022, calculate the surface evapotranspiration ET (W / m 2 ), the calculation formula is as follows:

[0173] ET=(R n -G) / (1+β) (50)

[0174] like Figure 4 As shown, the specific steps of step 3 are:

[0175] Step 3011: Prepare the corresponding remote sensing evapotranspiration ET (W / m 2 ) GLDAS model evapotranspiration ET r (kg / m 2 / s), corresponding to the GLDAS model evaporation on the observation day of evaporation pan and groundwater level {ET Oi (kg / m 2 / s)} 1≤i≤8 ;

[0176] Step 3012: Convert remote sensing evapotranspiration ET (W / m 2 ) is converted to kg / m 2 / s;

[0177] Step 3013: Convert the dimension of remote sensing evapotranspiration ET (kg / m 2 / s), model evapotranspiration ET r (kg / m 2 / s), model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 Perform spatiotemporal fusion to obtain the corresponding evaporation of the evaporation pan, the observation day of the groundwater level, and the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 , the fusion method is:

[0178]

[0179] Where W is the size of the small window moving on the remote sensing evapotranspiration and model evapotranspiration, k, l are the positions of the pixels in the small window, and w kl is the weight at pixel (k, l), and satisfies 0≤w≤1,∑w=1;

[0180] Step 3014: Combine the evaporation Ri (kg / m 2 / s)} 1≤i≤8 Dimension conversion to integrated evapotranspiration ET R (mm / d);

[0181] like Figure 5 As shown, the specific steps of step 4 are:

[0182] Step 4011: Observe the evaporation E00 (mm / d) using the evaporation dish and use the inverse S-shaped curve to establish the fused evapotranspiration ET R The relationship model between (mm / d) and groundwater level H (m) is as follows:

[0183]

[0184] Step 4012: Convert equation (52) to the following form:

[0185]

[0186] Step 4013: Using the observation data of E00 and H and the fused evapotranspiration ET R Perform a linear fit on equation (53), and the fitting steps are as follows:

[0187] a) Using E00 and ET R The dependent variable is calculated based on the data Calculate the independent variable x=lnH and use least squares regression to establish the following fitting equation:

[0188] y=kx+m (54)

[0189] b) Using the fitting results from step a), calculate the parameters c and d in equation (52):

[0190]

[0191] c) Substituting d and c obtained in step b) into equation (52), we obtain the relationship model between phreatic evaporation and groundwater level:

[0192]

[0193] In step 5, the fused evapotranspiration ET obtained in step 3 is R Substitute the relationship model between surface evapotranspiration and groundwater level obtained in step 4 to calculate the groundwater level. The calculation formula is as follows:

[0194]

[0195] like Figure 6 As shown, a device for implementing a regional groundwater level monitoring method by fusing multi-source data includes 1. a data preparation module, 2. a remote sensing evapotranspiration estimation module, 3. a remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module, 4. a surface evapotranspiration and groundwater level modeling module, and 5. a groundwater level estimation module;

[0196] The data preparation module 1 is used to determine the study area, download the visible light band and thermal infrared band of remote sensing images covering the study area, and the GLDAS model evapotranspiration product, and perform cropping and projection conversion, and collect evaporation pan and groundwater level observation data of the study area;

[0197] The remote sensing evapotranspiration estimation module 2 is used to estimate surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of the remote sensing image;

[0198] The remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module 3 is used to perform spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into fused evapotranspiration at a day scale;

[0199] The surface evapotranspiration and groundwater level modeling module 4 is used to establish a connection model between surface evapotranspiration and groundwater level based on the observation data of evaporation pan and groundwater level and the fused evapotranspiration at the day scale;

[0200] The groundwater level estimation module 5 is used to input the fused evapotranspiration into the connection model of surface evapotranspiration and groundwater level to estimate the regional groundwater level;

[0201] like Figure 7 As shown, the remote sensing evapotranspiration estimation module 2 includes a vegetation index calculation unit 201, a vegetation coverage calculation unit 202, an emissivity calculation unit 203, a surface temperature calculation unit 204, an albedo calculation unit 205, a surface net radiation flux calculation unit 206, a surface heat flux calculation unit 207, a surface temperature and vegetation coverage two-dimensional scattered point space construction unit 208, a Bowen ratio calculation unit 209, and a surface evapotranspiration remote sensing estimation unit 210;

[0202] The vegetation index calculation unit 201 is used to remove cloud-contaminated pixel data through built-in quality control tags of remote sensing images; and calibrate the remote sensing data through the gain and offset of the remote sensing data;

[0203] The normalized difference vegetation index (NDVI) (dimensionless) is calculated using the surface reflectance NIR of the near infrared band and the surface reflectance R of the red light band of the remote sensing image data. The NDVI calculation formula is as follows:

[0204]

[0205] The vegetation coverage calculation unit 202 is used to determine the maximum NDVI value NDVI of the study area based on the soil moisture characteristics and vegetation coverage characteristics of the study area. max and minimum NDVI value NDVI min ;

[0206] NDVI obtained by calculating the vegetation index unit and the determined NDVI max and NDVI min Calculate the vegetation coverage FVC of the study area. The FVC calculation formula is as follows:

[0207]

[0208] The emissivity calculation unit 203 is used to determine the emissivity ε of the bare soil coverage of the study area according to the soil type characteristics and vegetation coverage characteristics of the study area. soil (dimensionless) and the emissivity ε of full vegetation cover veg (dimensionless);

[0209] According to the vegetation coverage calculation unit obtained FVC and the determined ε soil and ε veg Calculate the surface emissivity EM (dimensionless) using the following formula:

[0210] EM=FVC*ε veg +(1-FVC)*ε soil (60)

[0211] The surface temperature calculation unit 204 is used to calculate the brightness temperature data T of the remote sensing image based on the EM obtained by the emissivity calculation unit. B Calculate the true surface temperature T s (K), calculated as follows:

[0212]

[0213] The albedo calculation unit 205 is used to calculate the blue light band reflectivity ρ of the remote sensing image obtained by the vegetation index calculation unit. B (dimensionless), red light band reflectivity ρ R (dimensionless), NIR band reflectivity ρ NIR (dimensionless), mid-wave infrared reflectivity ρ swir (dimensionless), mid-wave infrared reflectivity ρ swir2 (dimensionless) Calculate the albedo α (dimensionless) using the following formula:

[0214] α=(0.356*ρ B +0.13*ρ R +0.373*ρ NIR +0.085*ρ swir +0.072*ρ swir2 -0.018) / 1.016(62)

[0215] The surface net radiation flux calculation unit 206 is used to calculate the EM (dimensionless) obtained by the emissivity calculation unit and the T obtained by the surface temperature calculation unit. s (K), α (dimensionless) obtained by the albedo calculation unit, combined with the solar downward shortwave radiation S d (W / m 2 ), Stefan-Boltzmann constant σ (dimensionless), atmospheric incident longwave radiation (W / m 2 ), calculate the surface net radiation flux R n (W / m 2 ), the calculation formula is as follows:

[0216]

[0217] Among them, the solar downward shortwave radiation S d The calculation formula is as follows:

[0218] S d =S0*cosZ*cosZ / (1.085*cosZ+ea*(2.7+cosZ)*0.001+b) (64)

[0219] Among them, S0 is the solar radiation constant at the top of the atmosphere, which is 1367W / m 2 , Z is the solar zenith angle obtained from the image data, in radians, ea is the atmospheric water vapor pressure, obtained from meteorological observations, in hPa, and b = 0.2 is the correction factor;

[0220] Where α is the albedo, dimensionless, σ is the Stefan-Boltzmann constant, which is 0.000000056697, EM is the emissivity of the surface, dimensionless, and T s is the surface temperature in K, ε sky is the sky emissivity, which can be approximately equal to 1, T sky is the sky temperature in K;

[0221] The surface heat flux calculation unit 207 is used to calculate the surface net radiation flux R according to the FVC obtained by the vegetation coverage calculation unit and the surface net radiation flux obtained by the surface net radiation flux calculation unit. n , when the soil is completely bare, the surface heat flux G accounts for R n The ratio of G to R is 0.315. n The ratio is 0.05, and the surface heat flux G (W / m 2 ), the calculation formula is as follows:

[0222] G=R n [0.05+0.31*(1-FVC)] (65)

[0223] The surface temperature and vegetation coverage two-dimensional scattered point space construction unit 208 is used to calculate the FVC obtained by the vegetation coverage calculation unit and the T obtained by the surface temperature calculation unit. s , construct a model with FVC as the horizontal axis and T s The two-dimensional scattered point space of the vertical axis is divided into several equal parts, and the maximum T in each FVC partition is counted. s Tmaxv and minimum T s Tminv, with several equally divided FVC values as independent variables and the corresponding Tmaxv and Tminv as dependent variables, the following regression equation is obtained:

[0224] Tmaxv=ar+br*FVC (66)

[0225] Tminv=aw+bw*FVC (67)

[0226] Where ar and br are the slope and intercept of the dry side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively; aw and bw are the slope and intercept of the wet side of the two-dimensional scatter space of surface temperature and vegetation cover, respectively.

[0227] The Bowen ratio calculation unit 209 is used to calculate the FVC obtained by the vegetation coverage calculation unit and the T obtained by the surface temperature calculation unit. s The Tmaxv and Tminv obtained by constructing the unit in the two-dimensional scattered space of surface temperature and vegetation cover are used to calculate the Bowen ratio β (dimensionless) of each pixel. The calculation formula is as follows:

[0228] β=(T s -Tminv) / (Tmaxv-T s ) (68)

[0229] The surface evapotranspiration remote sensing estimation unit 210 is used to calculate the R n , G obtained by the surface heat flux calculation unit, β obtained by the Bowen ratio calculation unit, and the surface evapotranspiration ET (W / m 2 ), the calculation formula is as follows:

[0230] ET=(R n -G) / (1+β) (69)

[0231] like Figure 8 As shown, the remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module 3 includes a fusion data preparation unit 301, a remote sensing evapotranspiration dimension conversion unit 302, a spatiotemporal fusion unit 303, and a fusion evapotranspiration dimension conversion unit 304;

[0232] The fusion data preparation unit 301 is used to prepare the corresponding remote sensing evapotranspiration ET (W / m 2 ) GLDAS model evapotranspiration ET r (kg / m 2 / s), corresponding to the GLDAS model evaporation on the observation day of evaporation pan and groundwater level {ET Oi (kg / m 2 / s)} 1≤i≤8 ;

[0233] The remote sensing evapotranspiration dimension conversion unit 302 is used to convert the remote sensing evapotranspiration ET (W / m 2 ) is converted to kg / m 2 / s;

[0234] The time-space fusion unit 303 is used to convert the dimension-converted remote sensing evapotranspiration ET (kg / m 2 / s), model evapotranspiration ET r (kg / m 2 / s), model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 Perform spatiotemporal fusion to obtain the corresponding evaporation of the evaporation pan, the observation day of the groundwater level, and the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 , the fusion method is:

[0235]

[0236] Where W is the size of the small window moving on the remote sensing evapotranspiration and model evapotranspiration, k, l are the positions of the pixels in the small window, and w kl is the weight at pixel (k, l), and satisfies 0≤w≤1,∑w=1;

[0237] The fused evapotranspiration dimension conversion unit 304 is used to convert the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 Dimension conversion to integrated evapotranspiration ET R (mm / d);

[0238] like Figure 9 As shown, the surface evapotranspiration and groundwater level modeling module 4 includes an evaporation pan, groundwater level observation data and fusion evapotranspiration conversion unit 401 , a conversion data regression unit 402 , and a surface evapotranspiration and groundwater level modeling unit 403 .

[0239] The evaporation pan, groundwater level observation data and fusion evaporation conversion unit 401 is used to use the evaporation pan evaporation observation E00 (mm / d) to establish the fusion evaporation ET using the inverse S-shaped curve. R The relationship model between (mm / d) and groundwater level H (m) is as follows:

[0240]

[0241] The transformed data regression unit 402 is used to transform equation (71) into the following form:

[0242]

[0243] The surface evapotranspiration and groundwater level modeling unit 403 is used to use the observation data of E00 and H and the fusion evapotranspiration ET RPerform a linear fit on equation (72), and the fitting steps are as follows:

[0244] a) Using E00 and ET R The dependent variable is calculated based on the data Calculate the independent variable x=lnH and use least squares regression to establish the following fitting equation:

[0245] y=kx+m (73)

[0246] b) Using the fitting results from step a), calculate the parameters c and d in equation (71):

[0247]

[0248] Substituting d and c obtained in step b) into equation (71), we obtain the relationship model between surface evapotranspiration and groundwater level:

[0249]

[0250] The groundwater level estimation module 5 is used to combine the fusion evapotranspiration ET obtained by the remote sensing evapotranspiration and the model evapotranspiration spatiotemporal fusion module R Substitute the relationship model between surface evapotranspiration and groundwater level obtained in the surface evapotranspiration and groundwater level modeling module to calculate the groundwater level. The calculation formula is as follows:

[0251]

[0252] The above embodiments are only used to illustrate the preferred solutions of the present invention rather than to limit the same. A person skilled in the art may improve, modify or make equivalent replacements for the technical solutions of the present invention, which shall all fall within the scope of protection of the appended claims of the present invention.

Claims

1. A method for regional groundwater level monitoring by integrating multi-source data, characterized in that: The following steps are involved: Step 1: Determine the study area, download the visible and thermal infrared bands of remote sensing images of the study area and the Global Land Data Assimilation System (GLDAS) model evapotranspiration product, perform cropping and projection conversion, and collect evaporation pan and groundwater level observation data for the study area; Step 2: Estimate surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of remote sensing images; Step 3: Perform spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into fused evapotranspiration at the daily scale. The specific steps are as follows: Step 3011: Prepare the corresponding remote sensing evapotranspiration ET (W / m 2 ) GLDAS model evapotranspiration ET r (kg / m 2 / s), corresponding to the evaporation of the evaporation pan and the GLDAS model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 ; Step 3012: Convert remote sensing evapotranspiration ET (W / m 2 ) is converted to kg / m 2 / s; Step 3013: Convert the dimension of remote sensing evapotranspiration ET (kg / m 2 / s), model evapotranspiration ET r (kg / m 2 / s), model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 Perform spatiotemporal fusion to obtain the fusion evaporation {ET Ri (kg / m 2 / s)} 1≤i≤8 , the fusion method is: Where W is the size of the small window moving on the remote sensing evapotranspiration and model evapotranspiration, k, l are the positions of the pixels in the small window, and w kl is the weight at pixel (k, l), and satisfies 0≤w≤1,∑w=1; Step 3014: Combine the evaporation Ri (kg / m 2 / s)} 1≤i≤8 Dimension conversion to integrated evapotranspiration ET R (mm / d); Step 4: Build a model linking surface evaporation and groundwater level based on pan evaporation, groundwater level observations, and fused evapotranspiration at the daily scale. Step 5: Input the fused evapotranspiration into the surface evapotranspiration and groundwater level linkage model to estimate the regional groundwater level.

2. The method for regional groundwater level monitoring by integrating multi-source data according to claim 1, characterized in that: The specific steps of step 2 for estimating surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of the remote sensing image are as follows: Step 2011: Remove cloud-contaminated pixel data using built-in quality control tags in remote sensing images; calibrate the remote sensing data using the gain and offset of the remote sensing data; Step 2012: Calculate the Normalized Difference Vegetation Index (NDVI) (dimensionless) using the surface reflectance NIR in the near infrared band and the surface reflectance R in the red band of the remote sensing image data obtained in step 2011. The NDVI calculation formula is as follows: Step 2013: Determine the maximum NDVI value of the study area based on the soil moisture characteristics and vegetation cover characteristics of the study area max and minimum NDVI value NDVI min ; Step 2014: NDVI obtained in step 2012 and NDVI obtained in step 2013 max and NDVI min Calculate the vegetation coverage FVC (dimensionless) of the study area. The FVC calculation formula is as follows: Step 2015: Determine the emissivity ε of the bare soil coverage of the study area based on the soil type characteristics and vegetation coverage characteristics of the study area soil (dimensionless) and the emissivity ε of full vegetation cover veg (dimensionless); Step 2016: Based on the FVC obtained in step 2014 and the ε obtained in step 2015 soil and ε veg Calculate the surface emissivity EM (dimensionless) using the following formula: EM=FVC*ε veg +(1-FVC)*ε soil (4) Step 2017: Based on the brightness temperature data T of the EM and remote sensing images obtained in step 2016 B Calculate the true surface temperature T s (K), calculated as follows: Step 2018: The blue light band reflectivity ρ of the remote sensing image obtained in step 2011 B (dimensionless), red light band reflectivity ρ R (dimensionless), NIR band reflectivity ρ NIR (dimensionless), mid-wave infrared reflectivity ρ swir (dimensionless), mid-wave infrared reflectivity ρ swir2 (dimensionless) Calculate the albedo α (dimensionless) using the following formula: α=(0.356*ρ B +0.13*p R +0.373*ρ NIR +0.085*r swir +0.072*p swir2 -0.018) / 1.016(6) Step 2019: Based on the EM (dimensionless) obtained in step 2016 and the T obtained in step 2017 s (K), α (dimensionless) obtained in step 2018, combined with the solar downward shortwave radiation S d (W / m 2 ), Stefan-Boltzmann constant σ (dimensionless), atmospheric incident longwave radiation Calculate the surface net radiation flux R n (W / m 2 ), the calculation formula is as follows: Among them, the solar downward shortwave radiation S d The calculation formula is as follows: S d =S0*cosZ*cosZ / (1.085*cosZ+ea*(2.7+cosZ)*0.001+b) (8) Among them, S0 is the solar radiation constant at the top of the atmosphere, which is 1367W / m 2 , Z is the solar zenith angle obtained from the image data, in radians, ea is the actual atmospheric water vapor pressure, obtained from meteorological observations, in hPa, and b = 0.2 is the correction factor; Where α is the albedo, dimensionless, the Stefan-Boltzmann constant σ is 0.000000056697, EM is the emissivity of the surface, dimensionless, Ts is the surface temperature in K, ε sky is the sky emissivity, which can be approximately equal to 1, T sky is the sky temperature in K; Step 2020: Based on the FVC obtained in step 2014 and the net surface radiation flux R obtained in step 2019 n , when the soil is completely bare, the surface heat flux G accounts for R n The ratio of G to R is 0.

315. n The ratio is 0.05, and the surface heat flux G (W / m 2 ), the calculation formula is as follows: G=R n [0.05+0.31*(1-FVC)] (9) Step 2021: Based on the FVC obtained in step 2014 and the T obtained in step 2017 s , construct a model with FVC as the horizontal axis and T s The two-dimensional scattered point space of the vertical axis is divided into several equal parts, and the maximum T in each FVC partition is counted. s Tmaxv and minimum T s Tminv, with several equally divided FVC values as independent variables and the corresponding Tmaxv and Tminv as dependent variables, the following regression equation is obtained: Tmaxv=ar+br*FVC (10) Tminv=aw+bw*FVC (11) Where, ar and br are the slope and intercept of the dry side of the two-dimensional scatter space of surface temperature and vegetation coverage, respectively; aw and bw are the slope and intercept of the wet side of the two-dimensional scatter space of surface temperature and vegetation coverage, respectively; Step 2022: Based on the FVC obtained in step 2014 and the T obtained in step 2017 s , Tmaxv and Tminv obtained in step 2021, calculate the Bowen ratio β (dimensionless) of each pixel, and the calculation formula is as follows: β=(T s -Tminv) / ( Tmaxv-T s ) (12) Step 2023: R obtained according to step 2019 n , G obtained in step 2020, β obtained in step 2022, calculate the surface evapotranspiration ET (W / m 2 ), the calculation formula is as follows: ET=(R n -G) / (1+β) (13).

3. The method for regional groundwater level monitoring by integrating multi-source data according to claim 1, characterized in that: The specific steps of step 4 are: Step 4011: Observe the evaporation E00 (mm / d) using the evaporation dish and use the inverse S-shaped curve to establish the fused evapotranspiration ET R The relationship model between (mm / d) and groundwater level H (m) is as follows: Step 4012: Convert equation (14) into the following form: Step 4013: Using the observation data of E00 and H and the fused evapotranspiration ET R Perform a linear fit on equation (15), and the fitting steps are as follows: a) Using E00 and ET R The dependent variable is calculated based on the data Calculate the independent variable x=lnH and use least squares regression to establish the following fitting equation: y=kx+m(16) b) Using the fitting results from step a), calculate the parameters c and d in equation (14): c) Substitute d and c obtained in step b) into equation (14) to obtain the relationship model between surface evapotranspiration and groundwater level:

4. The method for regional groundwater level monitoring by integrating multi-source data according to claim 1, characterized in that: Step 5 is to combine the fusion evaporation ET obtained in step 3 R Substitute the relationship model between surface evapotranspiration and groundwater level obtained in step 4 to calculate the groundwater level. The calculation formula is as follows:

5. A device for implementing a regional groundwater level monitoring method for integrating multi-source data according to any one of claims 1 to 4, characterized in that: It includes data preparation module, remote sensing evapotranspiration estimation module, remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module, surface evapotranspiration and groundwater level modeling module, and groundwater level estimation module; The data preparation module is used to determine the study area, download the visible light band and thermal infrared band of remote sensing images covering the study area, and the GLDAS model evapotranspiration product, and perform cropping and projection conversion, and collect evaporation pan and groundwater level observation data in the study area; The remote sensing evapotranspiration estimation module is used to estimate surface evapotranspiration pixel by pixel based on the visible light band and thermal infrared band of the remote sensing image; The remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module is used to perform spatiotemporal fusion of remote sensing evapotranspiration and model evapotranspiration to expand remote sensing evapotranspiration into fused evapotranspiration at a day scale; The surface evapotranspiration and groundwater level modeling module is used to establish a connection model between surface evapotranspiration and groundwater level based on the observation data of evaporation pan and groundwater level and the fused evapotranspiration at the day scale; The groundwater level estimation module is used to input the fused evapotranspiration into a connection model between surface evapotranspiration and groundwater level to estimate the regional groundwater level.

6. The device for regional groundwater level monitoring method by fusing multi-source data according to claim 5, characterized in that: The remote sensing evapotranspiration estimation module includes a vegetation index calculation unit, a vegetation coverage calculation unit, an emissivity calculation unit, a surface temperature calculation unit, an albedo calculation unit, a surface net radiation flux calculation unit, a surface heat flux calculation unit, a surface temperature and vegetation coverage two-dimensional scattered point space construction unit, a Bowen ratio calculation unit, and a surface evapotranspiration remote sensing estimation unit; The vegetation index calculation unit is used to remove cloud-contaminated pixel data through built-in quality control tags of remote sensing images; and calibrate the remote sensing data through the gain and offset of the remote sensing data; The normalized difference vegetation index (NDVI) (dimensionless) is calculated using the surface reflectance NIR of the near infrared band and the surface reflectance R of the red light band of the remote sensing image data. The NDVI calculation formula is as follows: The vegetation coverage calculation unit is used to determine the maximum NDVI value NDVI of the study area based on the soil moisture characteristics and vegetation coverage characteristics of the study area. max and minimum NDVI value NDVI min ; NDVI obtained by calculating the vegetation index unit and the determined NDVI max and NDVI min Calculate the vegetation coverage FVC (dimensionless) of the study area. The FVC calculation formula is as follows: The emissivity calculation unit is used to determine the emissivity ε of the bare soil full coverage of the study area according to the soil type characteristics and vegetation coverage characteristics of the study area. soil (dimensionless) and the emissivity ε of full vegetation cover veg (dimensionless); The FVC obtained by the vegetation coverage calculation unit and the determined ε soil and ε veg Calculate the surface emissivity EM (dimensionless) using the following formula: EM=FVC*ε veg +(1-FVC)*ε soil (22) The surface temperature calculation unit is used to calculate the brightness temperature data T of the EM and remote sensing image obtained by the emissivity calculation unit. B Calculate the true surface temperature T s (K), calculated as follows: The albedo calculation unit is used to calculate the blue light band reflectivity ρ of the remote sensing image obtained by the vegetation index calculation unit. B (dimensionless), red light band reflectivity ρ R (dimensionless), NIR band reflectivity ρ NIR (dimensionless), mid-wave infrared reflectivity ρ swir (dimensionless), mid-wave infrared reflectivity ρ swir2 (dimensionless) Calculate the albedo α (dimensionless) using the following formula: α=(0.356*ρ B +0.13*p R +0.373*ρ NIR +0.085*r swir +0.072*p swir2 -0.018) / 1.016 (24) The surface net radiation flux calculation unit is used to calculate the EM (dimensionless) obtained by the emissivity calculation unit and the T obtained by the surface temperature calculation unit. s (K), α (dimensionless) obtained by the albedo calculation unit, combined with the solar downward shortwave radiation S d (W / m 2 ), Stefan-Boltzmann constant σ (dimensionless), atmospheric incident longwave radiation Calculate the surface net radiation flux R n (W / m 2 ), the calculation formula is as follows: Among them, the solar downward shortwave radiation S d The calculation formula is as follows: S d =S0*cosZ*cosZ / (1.085*cosZ+ea*(2.7+cosZ)*0.001+b) (26) Among them, S0 is the solar radiation constant at the top of the atmosphere, which is 1367W / m 2 , Z is the solar zenith angle obtained from the image data, in radians, ea is the actual atmospheric water vapor pressure, obtained from meteorological observations, in hPa, and b = 0.2 is the correction factor; Where α is the albedo, dimensionless, the Stefan-Boltzmann constant σ is 0.000000056697, EM is the emissivity of the surface, dimensionless, T s is the surface temperature in K, ε sky is the sky emissivity, which can be approximately equal to 1, T sky is the sky temperature in K; The surface heat flux calculation unit is used to calculate the surface net radiation flux R obtained by the vegetation coverage calculation unit and the surface heat flux calculation unit. n , when the soil is completely bare, the surface heat flux G accounts for R n The ratio of G to R is 0.

315. n The ratio is 0.05, and the surface heat flux G (W / m 2 ), the calculation formula is as follows: G=R n [0.05+0.31*(1-FVC)] (27) The surface temperature and vegetation coverage two-dimensional scattered point space construction unit is used to calculate the FVC obtained by the vegetation coverage calculation unit and the T obtained by the surface temperature calculation unit. s , construct a model with FVC as the horizontal axis and T s The two-dimensional scattered point space of the vertical axis is divided into several equal parts, and the maximum T in each FVC partition is counted. s Tmaxv and minimum T s Tminv, with several equally divided FVC values as independent variables and the corresponding Tmaxv and Tminv as dependent variables, the following regression equation is obtained: Tmaxv=ar+br*FVC (28) Tminv=aw+bw*FVC (29) Where, ar and br are the slope and intercept of the dry side of the two-dimensional scatter space of surface temperature and vegetation coverage, respectively; aw and bw are the slope and intercept of the wet side of the two-dimensional scatter space of surface temperature and vegetation coverage, respectively; The Bowen ratio calculation unit is used to calculate the FVC obtained by the vegetation coverage calculation unit and the T obtained by the surface temperature calculation unit. s The Tmaxv and Tminv obtained by constructing the unit in the two-dimensional scattered space of surface temperature and vegetation cover are used to calculate the Bowen ratio β (dimensionless) of each pixel. The calculation formula is as follows: β=(T s -Tminv) / ( Tmaxv-T s ) (30) The surface evapotranspiration remote sensing estimation unit is used to calculate the R n , G obtained by the surface heat flux calculation unit, β obtained by the Bowen ratio calculation unit, and the surface evapotranspiration ET (W / m 2 ), the calculation formula is as follows: ET=(R n -G) / (1+β) (31).

7. The device for regional groundwater level monitoring method by fusing multi-source data according to claim 5, characterized in that: The remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module includes a fusion data preparation unit, a remote sensing evapotranspiration dimension conversion unit, a spatiotemporal fusion unit, and a fusion evapotranspiration dimension conversion unit; The fusion data preparation unit is used to prepare the corresponding remote sensing evapotranspiration ET (W / m 2 ) GLDAS model evapotranspiration ET r (kg / m 2 / s), corresponding to the GLDAS model evaporation on the observation day of evaporation pan and groundwater level {ET Oi (kg / m 2 / s)} 1≤i≤8 ; The remote sensing evapotranspiration dimension conversion unit is used to convert the remote sensing evapotranspiration ET (W / m 2 ) is converted to kg / m 2 / s; The space-time fusion unit is used to convert the dimension-converted remote sensing evapotranspiration ET (kg / m 2 / s), model evapotranspiration ET r (kg / m 2 / s), model evapotranspiration {ET Oi (kg / m 2 / s)} 1≤i≤8 Perform spatiotemporal fusion to obtain the corresponding evaporation of the evaporation pan, the observation day of the groundwater level, and the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 , the fusion method is: Where W is the size of the small window moving on the remote sensing evapotranspiration and model evapotranspiration, k, l are the positions of the pixels in the small window, and w kl is the weight at pixel (k, l), and satisfies 0≤w≤1,∑w=1; The fused evapotranspiration dimension conversion unit is used to convert the fused evapotranspiration {ET Ri (kg / m 2 / s)} 1≤i≤8 Dimension conversion to integrated evapotranspiration ET R (mm / d).

8. The device for regional groundwater level monitoring method by fusing multi-source data according to claim 5, characterized in that: The surface evapotranspiration and groundwater level modeling module includes evaporation pan, groundwater level observation data and fusion evapotranspiration conversion unit, conversion data regression unit, surface evapotranspiration and groundwater level modeling unit; The evaporation pan evaporation, groundwater level observation data and fusion evaporation conversion unit are used to use the evaporation pan evaporation observation E00 (mm / d) to establish the fusion evaporation ET using the inverse S-shaped curve. R The relationship model between (mm / d) and groundwater level H (m) is as follows: The transformed data regression unit is used to transform equation (33) into the following form: The surface evapotranspiration and groundwater level modeling unit is used to use the observation data of E00 and H and the fusion evapotranspiration ET R Perform a linear fit on equation (34), and the fitting steps are as follows: a) Using E00 and ET R The dependent variable is calculated based on the data Calculate the independent variable x=lnH and use least squares regression to establish the following fitting equation: y=kx+m(35) b) Using the fitting results from step a), calculate the parameters c and d in equation (33): c) Substitute d and c obtained in step b) into equation (33) to obtain the relationship model between surface evapotranspiration and groundwater level:

9. The device for regional groundwater level monitoring method by fusing multi-source data according to claim 5, characterized in that: The groundwater level estimation module is used to combine the fusion evapotranspiration ET obtained by the remote sensing evapotranspiration and model evapotranspiration spatiotemporal fusion module R Substitute the relationship model between surface evapotranspiration and groundwater level obtained in the surface evapotranspiration and groundwater level modeling module to calculate the groundwater level. The calculation formula is as follows:

Images