Farmland root zone soil humidity prediction method, system and equipment and storage medium

Abstract

The invention provides a farmland root zone soil humidity prediction method, system and device and a storage medium, and belongs to the technical field of soil humidity prediction, and the method comprises the steps: obtaining multi-source remote sensing data of a farmland; downscaling the spatial resolution of the multi-source remote sensing data to obtain a kernel normalized difference vegetation index kNDVI of a first spatial resolution, a surface soil humidity SSM of a second spatial resolution and meteorological data ERA5 of a third spatial resolution; a root zone soil humidity RZSMt at the t moment is given, the kNDVI, the SSM and the ERA5 at the same moment and the SSM at the t + 1 moment are jointly input into a soil humidity prediction model, and the root zone soil humidity RZSMt + 1 at the t + 1 moment is simulated; and taking the current model output as the model input of the next step, and repeating the step to complete the prediction of the soil humidity of the root zone every day. According to the soil humidity prediction model provided by the invention, the problems of sensitivity and complexity in remote sensing inversion can be effectively solved, and relatively high-precision irrigation area scale root zone soil humidity prediction is realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of soil moisture prediction, and in particular relates to a method, system, device and storage medium for predicting soil moisture in the root zone of a farmland. Background Art

[0002] Soil moisture (SM) is a key hydrological state variable that plays a vital role in connecting and regulating the land surface, soil and atmosphere, including evaporation, transpiration, infiltration and runoff. In particular, root zone soil moisture (RZSM) directly affects the water supply of plants, which in turn affects agricultural productivity and ecosystem health. In recent years, the rapid development of remote sensing technology has enabled it to monitor soil moisture quickly, widely and continuously, providing an effective means for obtaining soil moisture information on a large scale.

[0003] At present, the main remote sensing methods include passive microwave remote sensing, active microwave remote sensing, and optical and infrared remote sensing. Passive microwave remote sensing inverts soil moisture by receiving microwave signals naturally radiated from the surface, and has the advantages of wide coverage and all-weather and all-day monitoring. Active microwave remote sensing uses synthetic aperture radar (SAR) to transmit and receive microwave signals, and can provide soil moisture data with high spatial resolution. Optical and infrared remote sensing invert soil moisture through vegetation indices (such as normalized difference vegetation index NDVI) and surface energy balance models (such as SEBS and TSEB, etc.). However, the main limitation of these methods is that their signal penetration depth is limited, and usually only soil moisture information of a few centimeters above the surface can be obtained, resulting in a high inversion accuracy of surface soil moisture, but it is difficult to directly obtain soil moisture information at a deeper level.

[0004] In order to solve the problem of limited signal penetration depth, scientists have developed a variety of methods to simulate root zone soil moisture using surface soil moisture data. Among them, the recursive exponential filter (REF) is a commonly used technique that combines current and historical soil moisture data and uses a recursive algorithm to provide a smoothed and estimated soil profile moisture distribution. The principle of maximum entropy (POME) also plays an important role in soil moisture simulation. This method is based on the concept of entropy and estimates the distribution of soil moisture by inferring the most balanced probability distribution by maximizing the entropy value of the system. In the case of insufficient or incomplete data, POME can provide a reasonable estimate of soil moisture and derive possible moisture distribution using existing data and boundary conditions.

[0005] Ecosystem models also play an important role in simulating root zone soil moisture. Crop growth models (such as DSSAT and APSIM models), watershed hydrological models (such as SWAT and VIC models), and land surface process models (such as NOAH and CLM models) provide comprehensive soil moisture simulations at different scales. Through data assimilation techniques, including Kalman filtering, particle filtering, and variational assimilation, remote sensing observation data and model simulation data are integrated to update and correct soil moisture estimates in real time. These methods can better understand the spatiotemporal changes in soil moisture, thereby improving the reliability of root zone soil moisture estimates.

[0006] Of course, all of the above methods have certain limitations. The advantage of the recursive exponential filter is that it is simple to calculate and easy to implement, but it is sensitive to the selection of initial conditions and filter parameters and requires fine-tuning to ensure accuracy, especially when dealing with drastic changes or soil moisture anomalies. Although POME performs well in dealing with data-deficient conditions, its results depend on the boundary conditions and assumptions set. In practical applications, it is necessary to ensure the rationality of these conditions and assumptions to avoid inaccurate estimates. Data assimilation technology can effectively improve the estimation accuracy, but the computational complexity is high when facing larger-scale scenarios, and the requirements for data and models are also strict, requiring powerful computing resources and efficient algorithm support. Summary of the invention

[0007] In order to solve the problems of the current remote sensing inversion method of root zone soil moisture in initial conditions, model parameter sensitivity and calculation complexity, the present invention provides a method for predicting farmland root zone soil moisture.

[0008] In order to achieve the above object, the present invention provides the following technical solutions:

[0009] A method for predicting soil moisture in a farmland root zone, comprising:

[0010] Obtain multi-source remote sensing data of farmland;

[0011] The spatial resolution of MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data is downscaled to obtain the kernel normalized difference vegetation index kNDVI of the first spatial resolution; the spatial resolution of MODIS LST land surface temperature and ERA5 LST land surface temperature in multi-source remote sensing data is downscaled to obtain the surface soil moisture SSM of the second spatial resolution; the spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data is downscaled to obtain the meteorological data ERA5 of the third spatial resolution;

[0012] The features in the soil water balance equation are replaced with various remote sensing data to construct a soil moisture prediction model; given the root zone soil moisture RZSM_t at time t, the kNDVI, SSM and ERA5 at the same time and the SSM at time t+1 are input into the soil moisture prediction model to simulate the root zone soil moisture RZSM_t+1 at time t+1; the current model output will be used as the model input for the next step, and this step will be repeated to complete the prediction of daily root zone soil moisture.

[0013] Preferably, the downscaling operation specifically includes:

[0014] The kernel normalized difference vegetation index (kNDVI) of the same vegetation was calculated using MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data, and the two calculated kNDVIs were fused to generate the kernel normalized difference vegetation index (kNDVI_20m_Daily) with a daily spatial resolution of 20m.

[0015] The MODIS LST land surface temperature and ERA5 LST land surface temperature in multi-source remote sensing data are fused to obtain the spatiotemporal continuous land surface temperature LST_1km_Daily with a daily spatial resolution of 1km; the vegetation temperature condition index VTCI_1km_Daily with a daily spatial resolution of 1km is generated by combining LST_1km_Daily with the kNDVI generated by MODIS multispectral images; the surface soil moisture SMAP SSM in multi-source remote sensing data is spatially downscaled using VTCI_1km_Daily to obtain the surface soil moisture SSM_1km_Daily with a daily spatial resolution of 1km;

[0016] The spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data is downscaled to obtain the daily meteorological data ERA5_20m_Daily with a spatial resolution of 20m.

[0017] Preferably, before fusing the two calculated kNDVIs to generate the kNDVI_20m_Daily image, it also includes using cubic spline interpolation to fill in the missing pixels of the kNDVI calculated using the Sentinel-2 multispectral image, and performing SG filtering on the kNDVI to complete the reconstruction of the kNDVI;

[0018] The reconstructed kNDVI was fused with the kNDVI calculated using MODIS multispectral images to generate the kernel normalized difference vegetation index kNDVI_20m_Daily with a daily spatial resolution of 20m.

[0019] Preferably, the two calculated kNDVIs are fused using the ESTARFM algorithm to generate a kernel normalized difference vegetation index kNDVI_20m_Daily with a daily spatial resolution of 20 m, and the MODIS LST land surface temperature and the ERA5 LST land surface temperature in the multi-source remote sensing data are fused to obtain the spatiotemporal continuous surface temperature LST_1km_Daily with a daily spatial resolution of 1 km.

[0020] Preferably, the surface soil moisture SMAP SSM in the multi-source remote sensing data is spatially downscaled using VTCI_1km_Daily to obtain the surface soil moisture SSM with a daily spatial resolution of 1km, specifically:

[0021]

[0022] Where SSM refers to the surface soil moisture at a resolution of 1 km after downscaling. is the original 9km resolution surface soil moisture SMAP SSM, VTCI is the vegetation temperature condition index VTCI_1km_Daily with a daily spatial resolution of 1km, is the average value of VTCI within a SMAP 9km grid box;

[0023]

[0024] Where m and n represent the mth row and nth column of the 1 km grid in the SMAP 9 km grid frame. The calculation formula of VTCI is as follows:

[0025]

[0026] In the formula, They respectively represent that in the study area, when kNDVI i The maximum and minimum values ​​of LST when it is equal to a certain value; Indicates that the kNDVI value of a pixel is kNDVI i a, b, a', b' are unknown coefficients, which are approximated by the scatter plot of kNDVI and LST in the farmland area.

[0027] Preferably, the root zone soil moisture RZSM_t at the given time t, and the kNDVI_20m_Daily, SSM_1km_Daily and ERA5_20m_Daily at the same time and SSM_1km_Daily at the time t+1 are input into the soil moisture prediction model to simulate the root zone soil moisture RZSM_t+1 at the time t+1. The specific calculation formula is:

[0028] AWF t+1=f(AWF t ,P t ,SSM t+1 ,SSM t ,RH t ,PET t ,kNDVI t );

[0029] Among them, AWF t+1 represents the available water in the root zone at time t+1, AWF t represents the available water in the root zone at time t, P t represents the rainfall at time t, SSM t+1 represents SSM_1km_Daily at time t+1, SSM t represents SSM_1km_Daily at time t, RH t Relative humidity in ERA5_20m_Daily at time t, PET t represents the potential evapotranspiration at time t, KNDVI t represents kNDVI_20m_Daily at time t;

[0030] RZS t+1 =AWF t+1 ×(FC-WP)+WP;

[0031] Among them, RZSM t+1 is the soil moisture in the root zone at time t+1, FC is the field water holding capacity, and WP is the wilting coefficient.

[0032] The present invention also provides a farmland root zone soil moisture prediction system, comprising:

[0033] Data acquisition module, used to obtain multi-source remote sensing data of farmland;

[0034] The downscaling module is used to downscale the spatial resolution of MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data to obtain the kernel normalized difference vegetation index kNDVI of the first spatial resolution; downscale the spatial resolution of MODIS LST land surface temperature and ERA5LST land surface temperature in multi-source remote sensing data to obtain the surface soil moisture SSM of the second spatial resolution; downscale the spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data to obtain the meteorological data ERA5 of the third spatial resolution;

[0035] The humidity prediction module is used to replace the features in the soil water balance equation with various remote sensing data to build a soil humidity prediction model; given the root zone soil humidity RZSM_t at time t, the kNDVI, SSM and ERA5 at the same time and SSM_1km_Daily at time t+1 are input into the soil humidity prediction model to simulate the root zone soil humidity RZSM_t+1 at time t+1; the current model output will be used as the model input for the next step, and this step will be repeated to complete the prediction of daily root zone soil humidity.

[0036] The present invention also provides a computer device, comprising a memory, a processor and a computer program stored in the memory, wherein the processor executes the computer program to implement any one of the steps in the method for predicting soil moisture in the farmland root zone.

[0037] The present invention also provides a computer-readable storage medium, on which a computer program is stored. When the computer program is loaded by a processor, it can execute any one of the steps in the method for predicting soil moisture in the farmland root zone.

[0038] The method for predicting soil moisture in the root zone of farmland provided by the present invention has the following beneficial effects:

[0039] The present invention is based on the soil water balance equation, replaces the features in the soil water balance equation with various remote sensing data, and constructs a soil moisture prediction model, so that the model requires fewer parameters and is insensitive to initial conditions, and can complete the spatiotemporal continuous prediction of the soil moisture in the root zone of the study area, reducing the complexity of the prediction and the requirements for resources, and improving the prediction efficiency. By downscaling multi-source remote sensing data, the present invention does not rely on the set boundary conditions and assumptions. By fusing multi-source data based on the water balance equation, the sensitivity and complexity problems in remote sensing inversion can be effectively solved, and a relatively high-precision root zone soil moisture prediction at the irrigation area scale can be achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the embodiment of the present invention and its design scheme, the following briefly introduces the drawings required for this embodiment. The drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0041] Figure 1 This is a flow chart of a method for predicting soil moisture in the root zone of a farmland according to Embodiment 1 of the present invention;

[0042] Figure 2 The schematic diagram of the root zone soil moisture prediction model is shown in Figure 2. Figure 2 (a) is the technical route for data processing and model building. Figure 2(b) is a diagram of water dynamics in the soil-plant-atmosphere continuum;

[0043] Figure 3 This is a result diagram of verifying the interpolated relative humidity and potential evapotranspiration using the in-situ observation data of site A; Figure 3 (a) is the relative humidity in 2018, Figure 3 (b) is the potential evapotranspiration in 2018, Figure 3 (c) is the relative humidity in 2019, Figure 3 (d) is the potential evapotranspiration in 2019;

[0044] Figure 4 The spatial distribution of kNDVI and LST at different time points and the verification diagram of downscaling results; Figure 4 (a)-(d) are the spatial distribution of 2018.8.1 and 2019.7.1 and the verification diagram of the downscaling results;

[0045] Figure 5 This is the relevant performance diagram of SMAP SSM in the study area; among them, Figure 5 (a)-(c) are the measured surface soil moisture at site A in 2018, site C in 2018 and 2019, and the downscaled SMAP surface soil moisture. Figure 5 (d) and (e) are the daily increments of SMAP surface soil moisture and precipitation at site A in 2018 and 2019, respectively. Figure 5 (f) and (g) are the local rainfall data of ERA5 products and the increment of surface soil moisture of SMAP in the study area on June 21, 2019, respectively. Figure 5 (h) is a scatter plot of the two;

[0046] Figure 6 is a comparison chart of soil moisture prediction and observation; Figure 6 (a)-(c) are the time series changes of soil moisture at sites A, C, and D in 2019, respectively. Figure 6 (d) and (e) are the comparisons between the predicted and observed soil moisture values ​​at 19 discrete observation points in the study area in 2018 and 2019, respectively;

[0047] Figure 7 The spatial and temporal distribution characteristics of available water in the root zone of the study area in 2018 and 2019; Figure 7 (a) is the spatial distribution of the mean value and coefficient of variation of AWF during the period of interest in 2018. Figure 7 (b) is an overview of the probability density function of the mean AWF in 2018. Figure 7(c) is an overview of the probability density function of the AWF coefficient of variation in 2018. Figure 7 (d) is the spatial distribution of the mean value and coefficient of variation of AWF during the period of interest in 2019. Figure 7 (e) is an overview of the probability density function of the mean AWF in 2019. Figure 7 (f) is an overview of the probability density function of the AWF coefficient of variation in 2019. Figure 7 (g) is the temporal variation of the mean and 90% confidence interval of daily AWF in the study area in 2018. Figure 7 (h) is the temporal variation of the mean and 90% confidence interval of daily AWF in the study area in 2019;

[0048] Figure 8 is the number of days for the model simulation results to converge after setting different AWF initial values; Figure 8 (a) is the time series of the model simulation values ​​when different AWF initial values ​​are set. Taking the F site in 2019 as an example, Figure 8 (b) and (c) are the spatial distribution and probability distribution of AWF convergence days in the study area in 2019;

[0049] Fig. 9 This is the attribution analysis of the impact of model parameters on soil moisture changes; Fig. 9 (a) is the feature density map based on the SHAP value. Fig. 9 (b) is the feature density map after stretching the SHAP value. DETAILED DESCRIPTION

[0050] In order to enable those skilled in the art to better understand the technical solution of the present invention and implement it, the present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and cannot be used to limit the scope of protection of the present invention.

[0051] Example 1

[0052] Based on the water balance equation, the present invention proposes a root zone soil moisture simulation strategy that is insensitive to initial conditions, has fewer model parameters, and is easy to use over a large spatial range. The water balance equation covers all major processes in the hydrological cycle (precipitation, evapotranspiration, runoff, and storage changes, etc.), providing a comprehensive perspective to understand the dynamics of water resources. The equation form is simple and easy to understand, and can be applied to systems of different scales (such as watersheds, irrigation areas, groundwater, etc.), and is easy to implement and calculate, which helps to improve the lightweight and generalization capabilities of the target model. Some remote sensing meteorological data sets, such as ERA5, MERRA-2, JRA-55, CFSR, and GFS Analysis, provide global climate reanalysis data covering multiple variables such as atmospheric composition, wind speed, temperature, and humidity. Projects such as SMAP, SMOS, and AMSR-E provide soil moisture products with global coverage through satellite remote sensing technology. Remote sensing ecological indices such as NDVI (Normalized Difference Vegetation Index), EVI (Enhanced Vegetation Index), SAVI (Soil Adjusted Vegetation Index), etc. provide information about vegetation coverage, health status and growth dynamics by analyzing the spectral characteristics of the surface. Although these products all have the contradiction of mismatch in temporal and spatial resolution, they still greatly promote scientific research and technological progress in related fields. If these remote sensing products are introduced into the water balance equation, it may be possible to simulate root zone soil moisture over a large spatial range. The above aspects will provide the core theoretical basis and construction ideas for the root zone soil moisture prediction model in the present invention.

[0053] Based on the soil water balance equation, the present invention combines multi-source remote sensing data and downscaling algorithms to construct a root zone soil moisture prediction model suitable for irrigation area scale. The goal is achieved through two main steps: first, multi-source remote sensing data is processed by downscaling to solve the contradiction between temporal and spatial resolution and provide suitable agents for each feature in the equation. The ESTARFM algorithm is used to generate kNDVI and LST data with high temporal and spatial resolution, and daily VTCI is calculated to improve the spatial resolution of SMAP SSM products, thereby capturing irrigation and runoff events, and using kNDVI to monitor crop growth. The Kriging interpolation method is used to generate daily meteorological data with a resolution of 20m. Secondly, a time series iteration model of root zone soil moisture is constructed. The model is calibrated using historical ground observation data, and the root zone soil moisture and multi-source remote sensing data at the current moment are used as input to predict the root zone soil moisture at the next moment, and the daily soil moisture of the study area is generated in a cyclic iteration.

[0054] The study area of ​​this example is a certain irrigation area (106°43'~107°27'E, 40°34'~41°14'N), with an altitude of 1032~1050m and a coverage area of ​​approximately 2.157*10 5 Hectares, including 1.543*10 5Hectares of arable land, the soil is mainly silt loam, loam and clay loam. The study area is located in arid, semi-arid and semi-desert steppe areas, with sufficient sunshine, an average annual sunshine duration of about 3180 hours, and an average frost-free period of 130 to 168 days. The climate here is mainly temperate continental, with cold and little snow in winter and hot and dry summers. The average annual temperature is about 6.3℃~7.7℃, the average precipitation is 138.2mm, and the average evaporation is 2096.4mm. Regional agricultural production is highly dependent on irrigation. Irrigation in the irrigation area is mainly based on irrigation with water from the Yellow River, with an annual water diversion volume of 1 billion to 1.2 billion cubic meters. Timely observation of soil moisture in the root zone will have a very important guiding significance for efficient water resource management in local agriculture.

[0055] This embodiment involves the soil water balance equation, which is widely used to describe the water dynamics in the soil-plant-atmosphere continuum and becomes the main basis for selecting remote sensing data sets in the present invention. Figure 2 As shown in (b), the root layer can be regarded as a container of water storage. Events such as irrigation, rainfall, capillary rise of groundwater and water vapor condensation have a supplementary effect on the soil water storage in the root zone. Events such as deep infiltration, soil evaporation and plant transpiration will consume the soil water storage in the root zone (soil moisture in the root zone). Runoff describes the flow of surface water in two-dimensional space after irrigation or rainfall events. The amount of water that rises from the groundwater surface to the root layer through the soil capillaries depends on the soil type, the depth of the water table and the soil moisture. According to the recommendations of FAO56, when the distance from the bottom of the root layer to the groundwater table is greater than 1 meter, the capillary rise water is assumed to be zero. The groundwater depth in the area of ​​the present invention is usually greater than 2.5 meters. The main crop grown is sunflower. The taproot generally grows about 0.3-1.2 meters deep in the soil, of which about 60% of the roots are distributed in the 0-40 cm soil layer. Therefore, the capillary rise water can be ignored. At this time, the soil water balance equation can be expressed as:

[0056] SM t+1 H=SM t H+P t +I t +N t -RO t -DP t -ET t (1);

[0057] Among them, SM t is the soil moisture on day t (m 3 / m 3 ), SM t+1 is the soil moisture on day t+1 (m 3 / m 3), H is the soil depth (mm), P is the precipitation (mm), I is the irrigation volume (mm), N is the amount of water vapor condensation (mm), RO is the runoff (mm), DP is the deep infiltration (mm), and ET is the actual evapotranspiration (mm).

[0058] The water balance equation needs to provide various features as described above, which are replaced by various remote sensing data in the present invention. This is why the present invention selects so many remote sensing data and performs downscaling.

[0059] Based on this, Figure 1 and Figure 2 As shown, the method for predicting soil moisture in the root zone of farmland provided by the present invention comprises the following steps:

[0060] Step 1, obtain multi-source remote sensing data of farmland. As shown in Table 1, the multi-source remote sensing data in this embodiment includes seven types of data sets. The seven types of data sets are: (1) ERA5-Land Daily Aggregated reanalysis data set, (2) MCD43A4.061 Nadir Bidirectional Reflectance DistributionFunctionAdjusted Reflectance (NBAR) product, (3) Sentinel-2L1C atmospheric apparent reflectance product, (4) SPL4SMGP.007 surface soil moisture product, (5) MOD11A1.061 land surface temperature product, (6) SoilGrids global gridded soil information, (7) Ground-based data.

[0061] Table 1 Datasets used to develop soil moisture downscaling methods

[0062]

[0063] (1)ERA5-Land DailyAggregated

[0064] ERA5-land is a reanalysis dataset that provides a consistent view of the evolution of land variables over several decades at a higher resolution than ERA5. ERA5-land is produced by replaying the land component of the ECMWF ERA5 climate reanalysis. The reanalysis uses the laws of physics to combine model data with observations from around the world into a globally complete and consistent dataset. The reanalysis produces data going back several decades, providing an accurate description of past climate. ERA5-Land Daily Aggregated is a daily summation of the ECMWF ERA5 land hourly assets, including both flowing and non-flowing bands. Flow bands are obtained by calculating the sum of the data for one day, while non-flow bands are formed by averaging all the hourly data for the day, and flow bands are marked with the "_sum" identifier. The dataset includes daily dynamics of 50 indices characterizing temperature, lakes, snow cover, soil moisture, radiation and heat, evapotranspiration and runoff, wind speed, air pressure, precipitation, and vegetation since 1950. According to the soil water balance equation, the present invention selects six indicators, namely Dewpointtemperature 2m, Temperature2m, Potential evaporation sum, Total precipitationsum and Skin temperature, from the ERA5-Land DailyAggregated dataset. Dewpoint temperature 2m and Temperature 2m are used to calculate relative humidity, attempting to serve as a proxy indicator for water vapor condensation. Skin temperature will assume the role of using ESTARFM to generate a coarse-resolution image at the daily scale of 1km LST. The data is obtained from the Google Earth Engine platform.

[0065] (2)MCD43A4.061

[0066] As shown in the soil water balance equation (1), actual evapotranspiration is an important component in the soil water balance equation. Its conversion with potential evapotranspiration is closely related to vegetation growth and soil water potential. Therefore, it is necessary to characterize vegetation growth and soil water potential. Multispectral data has a unique advantage in vegetation remote sensing. The MCD43A4 V6 Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (NBAR) product provides 500-meter reflectance data for MODIS "Land" bands 1-7. These data are adjusted using a bidirectional reflectance distribution function to simulate the values ​​collected from a nadir perspective. This dataset is used to calculate the daily continuous kNDVI of the study area and as a coarse-resolution input image for the ESTARFM algorithm. The calculation and download of kNDVI rely on the GEE platform.

[0067] (3) Sentinel-2L1C

[0068] The Sentinel-2 satellite is composed of two twin satellites, Sentinel-2A and Sentinel-2B. The revisit period of a single satellite is 10 days, which can be reduced to 5 days by the complementarity of the two satellites. The revisit period in the area of ​​the present invention is 2-3 days. The Sentinel-2 satellite carries a multispectral imager with 13 spectral bands, a width of 290km, and a ground resolution ranging from 10m to 60m depending on the band. Since the Sentinel-2L2A (S2L2A) surface reflectance product of the present invention area in 2018 is missing on the GEE platform, the Sentinel-2L1C (S2L1C) product is selected, and the atmospheric correction is performed to obtain the S2L2A product. At the same time, the clear sky images with cloud cover less than 20% are screened and the QA60 band is used for cloud removal. The use of this data set is to calculate the kernel normalized difference vegetation index kNDVI with a spatial resolution of 20m in the study area and as a fine resolution input image for the ESTARFM algorithm. Atmospheric correction, cloud removal, kNDVI calculation and download all rely on the GEE platform.

[0069] (4)SPL4SMGP.007

[0070] It is difficult to obtain the irrigation system of specific farmland within the irrigation area scale. However, as an important means of replenishing soil moisture, soil moisture after irrigation will increase to a certain extent. Therefore, the present invention attempts to capture the occurrence of regional irrigation events through surface soil moisture. SPL4SMGP.007 uses data assimilation technology to integrate satellite-borne L-band brightness temperature measurements and precipitation observations into the basin land surface process model, providing a number of terrestrial ecological indicators including surface soil moisture (0-5 cm vertical average). The present invention uses the GEE platform to obtain the SMAP Level-4 surface soil moisture product with a spatial resolution of 9 km and a temporal resolution set to output on a daily basis.

[0071] (5)MOD11A1.061

[0072] Farmland fragmentation and surface runoff make the amount of water brought by irrigation and rainfall events heterogeneous in spatial distribution. The feature space formed by VIS-LST is closely related to the surface vegetation cover and soil moisture conditions. This phenomenon provides an entry point for the spatial downscaling of surface soil moisture products. The MOD11A1.061 product provides daily 1km surface temperature. The MODIS LST land surface temperature dataset under clear sky conditions is screened and obtained on the GEE platform. It is then used as a fine-resolution input image for the ESTARFM algorithm to generate a daily-scale 1km LST image, which serves as the basis for the spatial downscaling of the SMAP SSM product.

[0073] (6)SoilGrids-global gridded soil information

[0074] The same soil moisture represents different soil water potential for soils of different textures, and soil water potential is an important factor affecting actual evapotranspiration. For the irrigation area scale, the spatial heterogeneity of soil cannot be ignored, which means that in the subsequent modeling process, the spatial heterogeneity of soil needs to be considered and soil moisture cannot be directly applied. SoilGrids is a global digital soil mapping system that uses the most advanced machine learning methods to map the spatial distribution of global soil properties. The prediction model is fitted using more than 230,000 soil profile observations from the WoSIS database and a series of environmental covariates. Covariates are selected from more than 400 environmental layers derived from earth observation products and other environmental information (including climate, land cover and terrain morphology). The output of SoilGrids is a global soil property map with six standard depth intervals (0-5cm, 5-15cm, 15-30cm, 30-60cm, 60-100cm and 100-200cm) with a spatial resolution of 250 meters. The present invention selects two indicators, Field Capacity and Wilting Point, in the first five depth intervals in SoilGrids to characterize the spatial differences in soil in the study area. Taking into account the spatial heterogeneity of soil, the present invention converts the root zone soil moisture RZSM into the root zone available water AWF through the field water holding capacity FC and the wilting coefficient WP to avoid affecting the model training results.

[0075]

[0076] (7) Ground-based data

[0077] The weather station at site A can provide long-term daily-scale meteorological data including daily maximum temperature, daily minimum temperature, relative humidity, wind speed at 2m, sunshine hours, and precipitation, and can calculate daily reference evapotranspiration according to the Penman formula. This part of the data is used to simply evaluate the performance of the interpolated ERA5-Land in the study area.

[0078] Step 2: Downscaling the multi-source remote sensing data of farmland, including:

[0079] Step 2.1, downscaling the spatial resolution of MODIS multispectral images and Sentinel-2 multispectral images in the multi-source remote sensing data to obtain the kernel normalized difference vegetation index kNDVI of the first spatial resolution.

[0080] Specifically, MODIS and Sentinel-2 multispectral images are used to calculate the kernel normalized difference vegetation index (kNDVI) of the same vegetation, and the two calculated kNDVIs are fused using the ESTARFM algorithm to generate the kernel normalized difference vegetation index (kNDVI_20m_Daily) with a daily spatial resolution of 20m.

[0081] Each parameter is described in time and space. Remote sensing data has the properties of spatial and temporal resolution. Spatial resolution can be simply understood as pixels. The more pixels a photo of the same scene has, the clearer it is, and the more detailed the depiction of space is. Temporal resolution means how often a photo is taken of the same scene. The higher the temporal resolution, the more photos there are, and it is easier to describe the temporal changes of a certain scene. The daily 20m spatial resolution in the present invention means that there is a photo with a pixel size of 20m*20m every day. Since the present invention is intended to generate daily root zone soil moisture, the temporal and spatial resolution of each feature (parameter) is particularly important. In addition, there are many types of images describing each feature, and different images have different temporal and spatial resolutions. It is usually difficult to have high temporal and spatial resolutions at the same time. Therefore, a downscaling method is needed to improve the temporal and spatial resolution of the image. MODIS images have a resolution of 500m per day, and Sentinel-2 has a resolution of 20m every five days. Therefore, the present invention uses the ESTARFM method to merge the two to obtain a daily 20m spatial resolution image. kNDVI is a feature extracted from the image.

[0082] When plants are under water stress, actual evapotranspiration is equal to the product of potential evapotranspiration, crop coefficient and water stress coefficient. Multispectral data is used to characterize the growth status of crops and attempts to serve as a proxy for crop coefficients to a certain extent. The kernel normalized difference vegetation index (kNDVI) is a vegetation remote sensing index used to enhance vegetation information extraction. Similar to the traditional NDVI (normalized difference vegetation index), both are used to measure the growth status and density of vegetation, but kNDVI can better handle nonlinear characteristics and vegetation information under complex backgrounds by introducing kernel functions. It has stronger distinguishing ability in complex scenes and can improve the robustness of vegetation indices, especially when processing high-resolution remote sensing data.

[0083]

[0084] Where NIR, red∈R are the reflectances of the near infrared and red spectral bands respectively, and the kernel function k measures the similarity between the two reflectances. The radial basis function (RBF) kernel function k(NIR, red) = exp(-(NIR-red) 2 / (2σ 2 )), the parameter σ is introduced to control the distance between the near-infrared band and the red light band, thereby adjusting the sensitivity to canopy density. Although it is possible to estimate the best σ in space or time, it is reasonable to set σ to the average value of 0.5. This simplifies the calculation formula of kNDVI to:

[0085] kNDVI=tanh(NDVI 2 ) (4);

[0086] Due to the cloud removal operation, there are local blanks in the Sentinel-2 kNDVI image, which will also cause missing values ​​in the local kNDVI image in the corresponding time period after fusion. At the same time, considering that image acquisition will inevitably bring noise, kNDVI will have slight fluctuations in time series. In the time dimension, the present invention uses cubic spline interpolation to fill the missing pixels and performs SG filtering to complete the reconstruction of kNDVI.

[0087] Step 2.2: Downscale the spatial resolution of MODIS LST land surface temperature and ERA5 LST land surface temperature in multi-source remote sensing data to obtain the surface soil moisture SSM of the second spatial resolution. Specifically include:

[0088] Step 2.2.1: Also use the ESTARFM algorithm to fuse the sunny day data of MODIS LST land surface temperature with the ERA5LST land surface temperature to obtain the daily surface temperature LST_1km_Daily with a spatial resolution of 1km.

[0089] In view of the "temporal and spatial contradiction" of remote sensing image resolution, many fusion methods have been successfully applied to generate various surface variables at reasonable temporal and spatial resolutions. ESTARFM fully considers the mixed pixel effect in similar pixel screening, similar pixel weight calculation, and coarse and fine resolution conversion coefficient calculation, and has good applicability in heterogeneous areas.

[0090] In order to further apply the ESTARFM algorithm, four main implementation steps are required: (1) determine similar pixels with similar spectral characteristics to the central pixel within a predetermined local search window; (2) calculate the weights of similar pixels; (3) determine the conversion coefficients between coarse pixels and fine pixels; (4) calculate the weights and conversion coefficients involved in predicting the fine-resolution reflectance of the central pixel on a predetermined date.

[0091] The present invention utilizes the algorithm to fuse the sunny day data (1km×1km) of MODIS LST and ERA5 LST (9km×9km) to obtain a 1km×1km temporally and spatially continuous surface temperature LST.

[0092] Step 2.2.2: Use the spatiotemporal continuous land surface temperature LST obtained in step 2.2 and the kNDVI generated by MODIS multispectral images to generate the vegetation temperature condition index VTCI_1km_Daily with a daily spatial resolution of 1km; use VTCI_1km_Daily to spatially downscale the SMAP SSM to obtain the surface soil moisture SSM_1km_Daily with a daily spatial resolution of 1km.

[0093] VTCI is constructed based on the feature space of NDVI and LST. In the present invention, NDVI in VTCI is replaced by kNDVI. The specific downscaling method is expressed as:

[0094]

[0095] The SSM here refers to the surface soil moisture at a resolution of 1 km after downscaling. is the original 9km resolution SMAP SSM, VTCI is the spatial scaling factor calculated by resampling the reflectivity of the MCD43A4.061 product to 1km resolution and combining it with the LST generated in step 2. is the average value of VTCI within a SMAP 9km grid box and can be calculated as follows:

[0096]

[0097] Where m and n represent the mth row and nth column of the 1km grid in the SMAP 9km grid box. The calculation formula of VTCI is as follows:

[0098]

[0099] In the formula, They respectively represent that in the study area, when kNDVI i The maximum and minimum values ​​of LST when it is equal to a certain value; Indicates that the kNDVI value of a pixel is kNDVI i a, a', b' are unknown coefficients, which are approximated by the scatter plot of kNDVI and LST in the study area.

[0100] The above method is used to obtain SSM_1km_Daily with a daily spatial resolution of 1km. Considering that in plain areas, the farm canal scale is usually used as an irrigation unit, 1km 2 The farmlands in the same area have the same irrigation system, so there is no need to increase the true spatial resolution of the SSM product to 20m. It is only necessary to resample SSM_1km_Daily to keep the grid size consistent with the other model input variables. Step 3 is implemented based on MATLAB 2023a.

[0101] Step 2.3, downscaling the spatial resolution of the ERA5-Land meteorological data in the multi-source remote sensing data to obtain the meteorological data ERA5 of the third spatial resolution.

[0102] Specifically: The ERA5-Land meteorological data is downscaled to the daily meteorological data ERA5_20m_Daily with a spatial resolution of 20m through Kriging interpolation.

[0103] The ERA5-Land dataset can provide meteorological data with a spatial resolution of 9 km, which does not match the target resolution of 20 m. Therefore, spatial interpolation of such datasets is required to meet the input requirements of the model.

[0104] Kriging is a regression algorithm that performs spatial modeling and prediction (interpolation) of random processes / random fields based on covariance functions. In a specific random process, Kriging can give the optimal linear unbiased estimate, so it is also called the spatial optimal unbiased estimator in geostatistics and is widely used in the fields of geography, environmental science, atmospheric science, etc.

[0105] The present invention uses Kriging method to interpolate ERA5-Land data based on MATLAB 2023a to achieve a spatial resolution of 20 m. The specific operation steps are: (1) converting ERA5-Land data to grid points; (2) calculating the semivariance value to measure the spatial correlation between points; (3) fitting the theoretical variogram model to construct the relationship between distance and semivariance; (4) using the variogram to estimate the attribute values ​​of unknown points; (5) exporting the estimation results as raster data types.

[0106] Step 3: Construct a soil moisture time series iteration model, also called a soil moisture prediction model.

[0107] Relative humidity (RH) indicates the saturation of water vapor in the air and is also the actual water vapor pressure (e) at the same temperature (T). a ) and saturated water vapor pressure (e 0 A ratio of (T)):

[0108]

[0109] The saturated water vapor pressure is a function of temperature and can be directly calculated from the temperature, see the Tetens empirical formula:

[0110]

[0111] Where, T is the air temperature (unit: °C), e 0 (T) is the saturated water vapor pressure at temperature T (unit: kPa).

[0112] The dew point temperature is the temperature at which the air is cooled so that the water vapor content in the air reaches saturation. a ) at dew point temperature (T dew ) is the saturated water vapor pressure:

[0113]

[0114] The soil water balance equation can be summarized as follows:

[0115] RZS t+1 =f(RZSM t ,P t ,I t ,N t ,RO t ,DP t ,ET t ) (10);

[0116] Where RZSM is the root zone soil moisture at time t, I is the irrigation amount (mm), N is the water vapor condensation amount (mm), RO is the runoff (mm), DP is the deep infiltration (mm), and ET is the actual evapotranspiration (mm).

[0117] Rainfall and irrigation events will replenish the root zone soil water volume. If the root zone soil can store water (compared to the current RZSM t and FC) plus the water consumption of farmland ET cannot fully account for the net water inflow (P t +I t -RO t ) is digested, the water content in the root zone may exceed the field water holding capacity. At this time, deep water seepage will occur in the root zone soil, and the amount of water seepage is equal to the difference between the two. Otherwise, deep seepage will not occur. Therefore, deep seepage is ET t , I t , RZSM t And the function of FC:

[0118] DP t =P t +I t -RO t -ET t -1000(FC-RZSM t );

[0119] The most easily observable evidence of irrigation, precipitation, and runoff events is the change in surface soil moisture, so the occurrence of irrigation and runoff events should be captured using surface soil moisture of the two days before and after the precipitation:

[0120] (I t ,RO t )=f(SSM t+1 ,SSM t ,P t );

[0121] Relative humidity is used as a proxy for the amount of water condensed:

[0122] N t =f(RH t );

[0123] Actual evapotranspiration is equal to potential evapotranspiration PET and crop coefficient K c and water stress coefficient K s The product of the three:

[0124] ET t =K c K s PET t ;

[0125] kNDVI is used to characterize the growth status of crops and attempts to serve as a proxy for crop coefficients to a certain extent:

[0126] K c=f(kNDVI t );

[0127] Taking into account the spatial heterogeneity of soil, the present invention converts soil moisture into available water in the root zone AWF through field water holding capacity FC and wilting coefficient WP to avoid affecting the model training results.

[0128]

[0129] Water stress is closely related to AWF:

[0130] K s =f(AWF t );

[0131] The soil moisture prediction model can be obtained by combining the above formulas. The specific model is:

[0132] AWF t+1 =f(AWF t ,P t ,SSM t+1 ,SSM t ,RH t ,PET t ,kNDVI t ) (12);

[0133] Among them, AWF t+1 represents the available water in the root zone at time t+1, AWF t represents the available water in the root zone at time t, P t represents the rainfall at time t, SSM t+1 represents SSM_1km_Daily at time t+1, SSM t represents SSM_1km_Daily at time t, RH t Relative humidity in ERA5_20m_Daily at time t, PET t represents the potential evapotranspiration at time t, KNDVI t Represents kNDVI_20m_Daily at time t.

[0134] It can be seen that the operation of the model requires the setting of AWF0 at the initial moment, and its value range is 0-1. The present invention sets its initial value to 0 in the entire study area, and runs this model with a step size of 1 day. The current model output will be used as the model input for the next step. At the same time, SoilGrids is used to convert the available water in the root zone into daily soil moisture in the root zone.

[0135] Three long-term observation stations can provide daily soil volumetric moisture content at different depths of 10 cm, 40 cm and 100 cm. The present invention extracts soil moisture data automatically measured by the stations from April 30 to September 30, 2018 and the same period in 2019, and calculates the root zone soil moisture to train and evaluate the soil moisture prediction model proposed by the present invention.

[0136] In addition to the in-situ data of soil surface moisture continuously observed at fixed points (these data refer to meteorological data and continuous soil moisture data, and in-situ data can be understood as data measured on the spot), a total of 14 soil volume moisture content measurements were carried out at 19 common observation points in the study area from April to September in 2018 and 2019. These ground observations will also be used to evaluate the performance of the soil moisture prediction model proposed in this invention.

[0137] Figure 2 The schematic diagram of the root zone soil moisture prediction model is shown in Figure 2. Figure 2 (a) is the technical route for data processing and model building. Figure 2 (b) is a diagram of water dynamics in the soil-plant-atmosphere continuum.

[0138] Step 3: Given the root zone soil moisture RZSM_t at time t, the kNDVI_20m_Daily, SSM_1km_Daily and ERA5_20m_Daily at the same time and SSM_1km_Daily at time t+1 are input into the soil moisture prediction model to simulate the root zone soil moisture RZSM_t+1 at time t+1, and this step is repeated; SoilGrids is used to convert the root zone available water into daily root zone soil moisture, and finally the daily root zone soil moisture prediction is completed.

[0139] The root zone soil moisture simulation model proposed in the present invention integrates multi-source remote sensing data based on a physical model. The required land surface parameters are only the field water holding capacity and wilting coefficient of the soil, which provides the model with a certain theoretical basis and improves the generalization ability of the model.

[0140] Attribution Analysis:

[0141] SHAP (SHapleyAdditive exPlanation) is a method for post-model explanation. Its core idea is to calculate the marginal contribution of features to the model output, and then explain the "black box model" from both the global and local levels. SHAP constructs an additive explanation model, and all features are regarded as "contributors". For each predicted sample, the model generates a predicted value, and the SHAP value is the value assigned to each feature in the sample. Its basic idea is to calculate the marginal contribution of a feature when it is added to the model, and then take into account the different marginal contributions of the feature in all feature sequences, and take the average, which is the SHAP value of the feature.

[0142]

[0143] Where SHAP_value i represents the SHAP value of the i-th feature, and E[f(x)] is the expectation of all samples f(x).

[0144] Results and Analysis:

[0145] First, evaluation of downscaling of remote sensing data

[0146] In order to ensure the scale consistency of various input variables, downscaling operations were performed on kNDVI, LST, SMAP, and ERA5-Land data. This will affect the accuracy of the final soil moisture prediction to a certain extent, so the accuracy of the downscaling results needs to be evaluated.

[0147] (1) Image data interpolation evaluation

[0148] Figure 3 This is a result diagram of verifying the interpolated relative humidity and potential evapotranspiration using the in-situ observation data of site A; Figure 3 (a) is the relative humidity in 2018, Figure 3 (b) is the potential evapotranspiration in 2018, Figure 3 (c) is the relative humidity in 2019, Figure 3 (d) is the potential evapotranspiration in 2019. The relative humidity (RH) and potential evapotranspiration (PET) of station A are used to evaluate the corresponding meteorological data after interpolation. Figure 3It can be seen that the interpolated relative humidity has a high correlation with the actual relative humidity in both 2018 and 2019, and can better represent the actual local relative humidity, with R values ​​of 0.881 and 0.829, and RMSE values ​​of 12.561 and 16.604, respectively. Similarly, the interpolated potential evapotranspiration can also satisfactorily represent the actual potential evapotranspiration, with the statistical characteristics of R values ​​of 0.770 and 0.769, and RMSE values ​​of 3.182 and 3.614, respectively. Since precipitation events do not occur frequently and the sample size is small, no definite conclusion can be drawn. The present invention pays more attention to the consistency between precipitation events and surface soil moisture data, as detailed in the following SMAP downscaling and evaluation of its ability to capture water events.

[0149] (2) kNDVI and LST downscaling evaluation

[0150] Figure 4 The figure below is the spatial distribution of kNDVI and LST and the verification diagram of downscaling results. The kNDVI of the local area is magnified, while LST shows the spatial distribution of the entire study area. Figure 4 As shown in the figure, a real fine image of kNDVI and LST was randomly selected in 2018 and 2019 for comparison with the downscaled image. The first column shows the MODIS-kNDVI and ERA5-LST images, which represent the coarse resolution images before downscaling. The second column shows the Sentinel-2 NDVI and MODIS LST images, which represent the real fine resolution images. The third column shows the kNDVI and LST images after downscaling. The fourth column shows the downscaling accuracy of kNDVI and LST in the form of a scatter plot.

[0151] from Figure 4 It can be seen that the downscaling results are highly consistent with the real images, which can better characterize the spatial heterogeneity of kNDVI and LST, and achieved satisfactory results. In 2018 and 2019, the correlation coefficient R between the downscaling results of kNDVI and the true value reached 0.969 and 0.991, respectively, and the RMSE was 0.051 and 0.031, respectively; the correlation coefficient R between the downscaling results of LST and the true value reached 0.869 and 0.842, respectively, and the RMSE was 4.739 and 2.950, respectively. This positive result is conducive to the subsequent construction of soil moisture prediction models.

[0152] (3) Evaluation of SMAP downscaling and its ability to capture water inflow events

[0153] The surface soil moisture at site B in 2018 and 2019 and site A in 2018 were selected to evaluate the downscaling results of SMAP. Figure 5This is the relevant performance diagram of SMAP SSM in the study area; among them, Figure 5 (a)-(c) are the measured surface soil moisture and downscaled SMAP surface soil moisture at site A in 2018, site C in 2018 and 2019, respectively. Figure 5 (d) and (e) are the daily increments of SMAP surface soil moisture and precipitation at site A in 2018 and 2019, respectively. Figure 5 (f) and (g) are the local rainfall data of ERA5 products and the increment of surface soil moisture of SMAP in the study area on June 21, 2019, respectively. Figure 5 (h) is a scatter plot of the two, with the horizontal axis being the increment of SMAP surface soil moisture and the vertical axis being rainfall. The R, Bias, MAE and RMSE of the downscaled SMAP SSM and the Linehe 2018 measured data are 0.669, -0.104, 0.104 and 0.110, respectively; the R, Bias, MAE and RMSE of the downscaled SMAP SSM and the measured data of site B in 2018 are 0.792, -0.092, 0.092 and 0.096, respectively; the R, Bias, MAE and RMSE of the measured data of site B in 2019 are 0.573, -0.079, 0.079 and 0.083, respectively. The overall high correlation with the measured data indicates that the downscaled SMAP SSM can well reflect the relative changes in soil moisture at different times and locations, but the high RMSE and Bias less than zero mean that the SMAP SSM has a non-negligible overestimation of the surface soil moisture in the study area. This is also reflected in the figure. The two have similar changing trends, but the measured surface soil moisture is significantly higher than SMAP SSM.

[0154] Due to the lack of actual local irrigation system, it is difficult to directly evaluate the ability of SMAP SSM to capture irrigation events. However, like precipitation, after an irrigation event occurs, the result on farmland soil is usually an increase in soil moisture. Therefore, precipitation events can be regarded as natural irrigation. In other words, the response of SMAP SSM to precipitation can be evaluated as a reference for its ability to capture irrigation events.

[0155] The present invention randomly plots the daily increment of SMAP SSM and precipitation events at a certain point in the study area ( Figure 5 (d, e)). It can be seen that when precipitation events occur, SMAP SSM is basically in an increasing state (ΔSSM>0), and SMAP SSM effectively responds to the occurrence of precipitation events. The increase in SMAP SSM when there is no precipitation should be caused by irrigation events. In addition, the increment of SMAP SSM also reflects the difference in precipitation in space, such as Figure 5As shown in (f)-(h), the increment of SMAP SSM is large where the precipitation is large, and vice versa. The spatial correlation coefficient R between the two reaches 0.482. This shows that although SMAP SSM underestimates the actual surface soil moisture in the study area, the high correlation between the two ensures that SMAP SSM can capture water events in the study area.

[0156] It should be added that in addition to water events, changes in surface soil moisture are also affected by previous soil moisture conditions, topography, and environmental factors. When soil moisture is high on the previous day, surface runoff and deep infiltration are more likely to occur, and the amount of water brought by precipitation or irrigation cannot be effectively absorbed by the soil, so soil moisture will not increase significantly. Figure 5 In (d), a large precipitation event occurred on July 22, 2018, and soil moisture also increased significantly, but the precipitation on July 23, 2018 did not increase soil moisture. In addition, when the precipitation is less than the evapotranspiration of the day, the soil moisture will not increase. The slope and tillage pattern of the farmland affect the distribution of water in space by affecting the runoff. At the same time, the water storage capacity of the soil has an upper limit, so that the soil moisture will be saturated when the water inflow reaches a certain level, such as Figure 5 The existence and combined effect of the above factors make it impossible to simply apply the water balance model when using remote sensing data to predict soil moisture. Instead, it is necessary to combine it with a machine learning algorithm to explore the relationship between the various factors in order to better predict soil moisture.

[0157] Second, soil moisture prediction results

[0158] (1) Performance of soil moisture prediction model

[0159] Figure 6 The figure is a comparison chart of soil moisture prediction and observation; (ac) are the time series changes of soil moisture at sites A, C and D in 2019, and (de) are the comparisons of soil moisture prediction and observation at 19 discrete observation points in the study area in 2018 and 2019. The 2019 soil moisture observations at sites A, C and D were used to verify the prediction results of soil moisture in terms of time series changes, as shown in Figure 2. Figure 6 As shown in (a)-(c) of Figure 2, the measured soil moisture data of 19 common observation points in 2018 and 2019 were used to verify the predicted results of soil moisture in terms of spatial distribution. The verification results are shown in Figure 6 As shown in (d) and (e).

[0160] For time series changes, the R of sites A, C and D reached 0.721, 0.538 and 0.433 respectively. The Bias was -0.001. The MAE was 0.005, 0.005 and 0.008 respectively. The RMSE was 0.007, 0.007 and 0.010 respectively. These results show that the changes in soil moisture in time series can be well predicted.

[0161] For spatial distribution, the R values ​​in 2018 and 2019 reached 0.659 and 0.535 respectively. The Bias values ​​were -0.050 and -0.017 respectively. The MAE values ​​were 0.008 and 0.017 respectively. The RMSE values ​​were 0.010 and 0.005 respectively. This result shows that the spatial distribution of soil moisture can also be well predicted.

[0162] The present invention sets the initial value of AWF to 0, which is inconsistent with the actual situation. Figure 6 It can be seen from (a)-(c) that there is a relatively large error between the predicted and observed values ​​of soil moisture in the early stage, which is particularly obvious at sites C and D before May 10, 2019. This shows that the assumption that the initial value of AWF is 0 will have a certain impact on the predicted results of soil moisture. However, in the middle and late stages, the predicted results of soil moisture do not seem to be affected, so the overall prediction accuracy is relatively good.

[0163] (2) Temporal and spatial distribution characteristics of soil moisture

[0164] Figure 7 The spatial and temporal distribution characteristics of available water in the root zone of the study area in 2018 and 2019; Figure 7 (a) is the spatial distribution of the mean value and coefficient of variation of AWF during the period of interest in 2018. Figure 7 (b) is an overview of the probability density function of the mean AWF in 2018. Figure 7 (c) is an overview of the probability density function of the AWF coefficient of variation in 2018. Figure 7 (d) is the spatial distribution of the mean value and coefficient of variation of AWF during the period of interest in 2019. Figure 7 (e) is an overview of the probability density function of the mean AWF in 2019. Figure 7 (f) is an overview of the probability density function of the AWF coefficient of variation in 2019. Figure 7 (g) is the temporal variation of the mean and 90% confidence interval of daily AWF in the study area in 2018. Figure 7(h) is the temporal variation of the mean and 90% confidence interval of daily AWF in the study area in 2019. The AWF in the focus period in 2018 was mainly concentrated around 0.12, and the difference in the study area was not obvious, such as Figure 7 As shown in (b), the AWF in the northeast is slightly higher than that in other regions. Figure 7 The same results were also shown in 2019, as shown in (a). Figure 7 The temporal coefficient of variation of AWF in 2018 is mainly concentrated around 0.14, but the coefficient of variation of AWF in a considerable number of regions is still much larger than this, such as Figure 7 As shown in (c), this part of the area is mainly distributed in the central and southern parts of the study area and at the foot of the mountains, such as Figure 7 In comparison, the temporal variability of AWF in 2019 is greater, with the coefficient of variation mainly concentrated around 0.18 and more homogeneous in spatial distribution, as shown in (b). Figure 7 As shown in (d) and (f).

[0165] It can be found that the two banks of the river have higher AWF and lower coefficient of variation, such as Figure 7 As shown in (a) and (d), it is shown that this area has always maintained a high soil moisture, which is consistent with the objective fact that the river water body infiltrates the soil in the riverside area. This phenomenon also verifies the accuracy of the soil moisture prediction of the present invention from the side.

[0166] In order to avoid the influence of the assumption that the initial value of AWF is 0, the description of the temporal changes of AWF starts from May 1st of each year. The AWF in 2018 decreased slowly from May 1st to late July, and showed a slow upward trend in mid-to-late September. Figure 7 The changes in AWF in 2019 during these two periods are almost the same as those in 2018. The differences between the two years are mainly concentrated in late July to early September. Figure 7 During this period, the AWF in 2018 first increased and then decreased, while the AWF in 2019 remained relatively stable, which was the result of the combined effects of multiple factors such as irrigation and meteorology.

[0167] (1) Influence of initial soil moisture value

[0168] Figure 8 is the number of days for the model simulation results to converge after setting different AWF initial values; Figure 8 (a) is the time series of the model simulation values ​​when different AWF initial values ​​are set. Taking the F site in 2019 as an example, Figure 8(b) and (c) are the spatial distribution and probability distribution of AWF convergence days in the study area in 2019. In order to explore the impact of the initial AWF value on the simulation results, 11 initial AWF values ​​from 0 to 1 were set with a step size of 0.1, and the AWF in the entire period of interest was simulated. Taking Hangjinhou Station in 2019 as an example, Figure 8 As shown in (a), it can be found that after 40 days, no matter what the initial value of AWF is set to, the simulation results of the model converge to the same value. The present invention defines the time spent in this process as the number of days of convergence. There may be two reasons for this result.

[0169] First, in the AWF prediction model constructed by the present invention, in addition to the AWF of the previous day, there are multiple factors that act together on the AWF at the current moment, and there are interactions between different factors, which makes it possible for input factors of different levels to produce the same simulation results. It can be observed that in the early stage of the simulation, the simulation results quickly converged in a certain range, and then the convergence speed slowed down, or even tended to be parallel. This may be because the value of AWF is directly related to the current evapotranspiration. The higher the AWF, the weaker the water stress on the vegetation, the greater the evapotranspiration, and the faster the water consumption. Conversely, the water consumption is slow. This interaction prompts the model to converge quickly, causing AWF to tend to a more stable range. When AWF is in this range, its weight affecting soil moisture changes is reduced, and the convergence speed of the model simulation results begins to slow down. At the same time, the present invention observed that a rainfall event occurred on the 40th day, that is, around the convergence date. In addition, the simulation results of the model with the initial AWF values ​​set to 0 and 0.1 respectively have converged around the 10th day. At almost the same time, the simulation results of the model with the initial AWF values ​​set to 0.3 to 1 have also converged. On this day, a rainfall event happened to occur. This shows that the rainfall event has accelerated the convergence of the model to a great extent. As an important supplementary source of soil moisture, this result is reasonable. Of course, rainfall is only easier to be discussed by the present invention than other factors, and it is not the only factor that affects the number of days to converge. For example, the simulation result with the initial AWF value of 0.2 did not converge with other results on the 10th day.

[0170] Second, the SSM, which is the downscaled SMAP product, may accidentally take on the role of calibrating the model in addition to being used to capture irrigation and runoff events. This role is similar to the principle of data assimilation in many studies: combining the surface soil moisture observed by remote sensing, updating the state and internal parameters of the model, avoiding the accumulation of errors, and simulating soil moisture more accurately. The differences between different AWF initial values ​​are essentially no different from the accumulated errors of the model. The use of SMAP SSM products constrains the simulation results of the model and improves the model's resistance to accumulated errors.

[0171] Figure 8 (b) and (c) show the convergence days and their probability distribution in the study area in 2019. It can be found that the convergence days in the study area are relatively fragmented in spatial distribution, and there is no obvious pattern, but they are not high overall, mostly concentrated below 40 days, which is not long compared to the period of interest, and can be ignored if it is in annual units. In addition, the present invention almost exhausts all possible initial values ​​of AWF, and can actually narrow the range of initial values ​​based on local experience, so that the model results can converge faster.

[0172] The above results mean that setting the initial AWF value to 0 does not affect the long-term soil moisture simulation.

[0173] (2) Attribution analysis of soil moisture changes

[0174] Fig. 9 This is the attribution analysis of the impact of model parameters on soil moisture changes; Fig. 9 (a) is the feature density map based on the SHAP value. Fig. 9 (b) is the feature density diagram after the SHAP value is stretched. In order to better understand the contribution of model parameters to soil moisture changes, the present invention uses SHAP to perform attribution analysis on the first-order difference of AWF (i.e., the difference between the two days before and after, representing the change in soil moisture). The larger the SHAP value, the greater the contribution of the feature to soil moisture changes. A negative value means a negative impact, and a positive value means a positive gain. Based on the absolute average value of the SHAP value, among the seven features used in the present invention, the contribution to soil moisture changes from large to small is as follows: AWF t-1 , SMAP t-1 , RH, kNDVI, SMAP t , Precipitation and PET, such as Fig. 9 As shown in (a). Among them, AWF t-1 and SMAP t-1 The contribution to soil moisture variation is much greater than that of other characteristics.

[0175] In order to more clearly observe the contribution of each feature, the SHAP value is stretched to -1 to 1, such as Fig. 9As shown in (b). It can be found that when AWFt-1 is at a lower value, it has a positive gain on the change of soil moisture, otherwise, it has a negative effect. When AWFt-1 is in a relatively middle range, its contribution to the change of soil moisture is almost 0. This confirms the interaction between AWF and evapotranspiration mentioned in the influence of the initial value of soil moisture. On the other hand, it may also be because climate and management measures limit the fluctuation range of soil moisture. As a feature of characterizing surface soil moisture, the functional similarity makes the performance of SMAPt-1 and AWFt-1 almost the same. According to the basis for establishing the model according to the present invention, SMAPt-1 and SMAPt should be analyzed together. Interestingly, the contribution of SMAPt to soil moisture change is completely opposite to that of SMAPt-1. When SMAPt is at a lower value, it has a negative impact on the change of soil moisture, otherwise, it has a positive gain. In the evaluation of SMAP downscaling and the ability to capture water events in the result analysis, it is mentioned that the difference between the two is used to capture water events. That is to say, when SMAPt is larger and SMAPt-1 is smaller, the water inflow on that day is larger, and the soil moisture increases more. Therefore, the opposite performance of the two in the attribution analysis results confirms that the model of the present invention has learned this law. Moreover, the positive or negative contribution of SMAPt to soil moisture change depends more on the size of SMAPt itself, which is likely to be the role of SMAPt in the calibration model as mentioned in the above-mentioned influence of the initial value of soil moisture. The larger the RH, the greater the possibility of water vapor condensation. The results of the attribution analysis also show that when the RH is higher, the greater its contribution to soil moisture change, and vice versa. The soil moisture in the root zone of farmland is mainly used for crop growth and development. Under the same conditions, the more luxuriant the crop growth, the greater the transpiration, the more water the root system extracts, and the faster the soil moisture is consumed. It can be found that when kNDVI is larger, the SHAP value is lower, and vice versa, the SHAP value is higher. It shows that it is reasonable for the present invention to use kNDVI as a proxy for the crop parameter Kc. The SHAP value of Precipitation is not high overall, mainly because it does not rain most of the time, and the absolute average algorithm is not favorable for ranking the contribution of Precipitation. However, once a rainfall event occurs, the contribution of Precipitation to the increase in soil moisture is obvious. It is worth noting that the performance of PET seems to be contrary to the expectations of the present invention. When PET is high, it actually has a positive gain on the change in soil moisture. In other cases, the contribution of PET does not show dependence on the characteristic value. This may be because irrigation managers will make irrigation decisions in response to higher PET, thereby increasing soil moisture.In addition, PET is not actual evapotranspiration, and its specific impact on soil moisture depends on the current soil moisture and crop growth conditions, which has been expressed through the contribution of AWFt-1 and kNDVI to soil moisture changes, which is why PET has the lowest contribution.

[0176] The results show that (1) multi-source remote sensing datasets can effectively replace relevant features and provide root zone soil moisture simulation with high temporal and spatial resolution and good simulation accuracy; (2) the mean values ​​of historical root zone available water AWF for two consecutive years are not significantly different in the study area, mainly concentrated around 0.12. The difference in the temporal change of the mean value of AWF is mainly concentrated from late July to early September; (3) the initial value of soil moisture has limited impact on the long-term simulation results; (4) the contributions to soil moisture change are: AWFt-1, SMAPt-1, RH, kNDVI, SMAPt, Precipitation and PET. The study shows that the soil moisture prediction model proposed in this paper can effectively solve the sensitivity and complexity problems in remote sensing inversion and achieve relatively high-precision root zone soil moisture prediction at the irrigation area scale.

[0177] Based on the soil water balance equation and integrating multi-source remote sensing data, this paper constructs a time series iterative model to simulate the soil moisture in the root zone of farmland, realizes the prediction of daily soil moisture in the root zone of the study area, and explores the influence of different characteristics on soil moisture changes.

[0178] Based on the same inventive concept, the present invention also proposes a farmland root zone soil moisture prediction system, comprising:

[0179] The data acquisition module is used to obtain multi-source remote sensing data of farmland.

[0180] The downscaling module is used to downscale the spatial resolution of MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data to obtain the kernel normalized difference vegetation index kNDVI of the first spatial resolution; downscale the spatial resolution of MODIS LST land surface temperature and ERA5LST land surface temperature in multi-source remote sensing data to obtain the surface soil moisture SSM of the second spatial resolution; downscale the spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data to obtain meteorological data ERA5 of the third spatial resolution.

[0181] The humidity prediction module is used to replace the features in the soil water balance equation with various remote sensing data to build a soil humidity prediction model; given the root zone soil humidity RZSM_t at time t, the kNDVI, SSM and ERA5 at the same time and SSM_1km_Daily at time t+1 are input into the soil humidity prediction model to simulate the root zone soil humidity RZSM_t+1 at time t+1; the current model output will be used as the model input for the next step, and this step will be repeated to complete the prediction of daily root zone soil humidity.

[0182] Each module in the above-mentioned farmland root zone soil moisture prediction system can be implemented in whole or in part by software, hardware and a combination thereof. Each of the above-mentioned modules can be embedded in or independent of a processor in a computer device in the form of hardware, or can be stored in a memory in a computer device in the form of software, so that the processor can call and execute the operations corresponding to each of the above modules.

[0183] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps in the method embodiment for predicting soil moisture in the root zone of a farmland. The specific implementation method can be found in the method embodiment, which will not be described in detail here.

[0184] Furthermore, the present invention also provides a non-temporary computer-readable storage medium containing instructions, and a computer program is stored on the storage medium. For example, a memory containing instructions, the above instructions can be executed by a processor of a computer device to complete the above method. For example, a non-temporary computer-readable storage medium can be a ROM, a random access memory (RAM), a CD-ROM, a magnetic tape, a floppy disk, and an optical data storage device. When the computer program is executed by the processor, the steps in the embodiment of the method for predicting soil moisture in the root zone of farmland can be implemented. The specific implementation method can be found in the method embodiment, which will not be repeated here.

[0185] It will be appreciated by those skilled in the art that embodiments of the present invention may provide methods, systems or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0186] The present invention is based on the soil water balance equation, integrates multi-source remote sensing data and combines multiple downscaling algorithms to characterize each variable in the equation, and constructs a root zone soil moisture prediction model suitable for irrigation district scale. The model requires fewer parameters and is insensitive to initial conditions. It can complete the spatiotemporal continuous prediction of the root zone soil moisture in the study area, reduce the complexity of the prediction and the requirements for resources, and improve the prediction efficiency. The boundary conditions of the present invention are provided by remote sensing data, such as ERA5 data. Through the fusion of multi-source data, the present invention does not rely on the set boundary conditions and assumptions. By fusing multi-source data based on the water balance equation, the sensitivity and complexity problems in remote sensing inversion can be effectively solved, and a relatively high-precision root zone soil moisture prediction at the irrigation district scale can be achieved.

[0187] It should be pointed out that the specific implementation methods described above can enable those skilled in the art to understand the invention more comprehensively, but do not limit the invention in any way. Therefore, although the invention has been described in detail in this specification and embodiments, those skilled in the art should understand that the invention can still be modified or replaced by equivalents; and all technical solutions and improvements that do not deviate from the spirit and scope of the invention are included in the protection scope of the patent for the invention. Any figure mark in the claims should not be regarded as limiting the claims involved. Any simple change or equivalent replacement of the technical solution that can be obviously obtained by any technician familiar with the field within the technical scope disclosed in the present invention belongs to the protection scope of the present invention.

Claims

1. A method for predicting soil moisture in the root zone of a farmland, characterized in that: include: Obtain multi-source remote sensing data of farmland; The spatial resolution of MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data is downscaled to obtain the kernel normalized difference vegetation index kNDVI of the first spatial resolution; the spatial resolution of MODIS LST land surface temperature and ERA5 LST land surface temperature in multi-source remote sensing data is downscaled to obtain the surface soil moisture SSM of the second spatial resolution; the spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data is downscaled to obtain the meteorological data ERA5 of the third spatial resolution; The features in the soil water balance equation are replaced with various remote sensing data to construct a soil moisture prediction model; given the root zone soil moisture RZSM_t at time t, the kNDVI, SSM and ERA5 at the same time and the SSM at time t+1 are input into the soil moisture prediction model to simulate the root zone soil moisture RZSM_t+1 at time t+1; the current model output will be used as the model input for the next step, and this step will be repeated to complete the prediction of daily root zone soil moisture.

2. The method for predicting soil moisture in the root zone of farmland according to claim 1, characterized in that: The downscaling operation specifically includes: The kernel normalized difference vegetation index (kNDVI) of the same vegetation was calculated using MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data, and the two calculated kNDVIs were fused to generate the kernel normalized difference vegetation index (kNDVI_20m_Daily) with a daily spatial resolution of 20m. The MODIS LST land surface temperature and ERA5 LST land surface temperature in multi-source remote sensing data are fused to obtain the spatiotemporal continuous land surface temperature LST_1km_Daily with a daily spatial resolution of 1km; the vegetation temperature condition index VTCI_1km_Daily with a daily spatial resolution of 1km is generated by combining LST_1km_Daily with the kNDVI generated by MODIS multispectral images; the surface soil moisture SMAP SSM in multi-source remote sensing data is spatially downscaled using VTCI_1km_Daily to obtain the surface soil moisture SSM_1km_Daily with a daily spatial resolution of 1km; The spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data is downscaled to obtain the daily meteorological data ERA5_20m_Daily with a spatial resolution of 20m.

3. The method for predicting soil moisture in the root zone of farmland according to claim 2, characterized in that: Before fusing the two calculated kNDVIs to generate the kNDVI_20m_Daily image, it also includes using cubic spline interpolation to fill in the missing pixels of the kNDVI calculated using the Sentinel-2 multispectral image, and performing SG filtering on the kNDVI to complete the reconstruction of the kNDVI; The reconstructed kNDVI was fused with the kNDVI calculated using MODIS multispectral images to generate the kernel normalized difference vegetation index kNDVI_20m_Daily with a daily spatial resolution of 20m.

4. The method for predicting soil moisture in the root zone of farmland according to claim 3, characterized in that: The two calculated kNDVIs were fused using the ESTARFM algorithm to generate the kernel normalized difference vegetation index kNDVI_20m_Daily with a daily spatial resolution of 20 m. The MODIS LST land surface temperature and the ERA5 LST land surface temperature in the multi-source remote sensing data were fused to obtain the spatiotemporal continuous land surface temperature LST_1km_Daily with a daily spatial resolution of 1 km.

5. The method for predicting soil moisture in the root zone of farmland according to claim 2, characterized in that: The VTCI_1km_Daily is used to spatially downscale the surface soil moisture SMAP SSM in the multi-source remote sensing data to obtain the surface soil moisture SSM with a daily spatial resolution of 1km, specifically: Where SSM refers to the surface soil moisture at a resolution of 1 km after downscaling. is the original 9km resolution surface soil moisture SMAP SSM, VTCI is the vegetation temperature condition index VTCI_1km_Daily with a daily spatial resolution of 1km, is the average value of VTCI within a SMAP 9km grid box; Where m and n represent the mth row and nth column of the 1 km grid in the SMAP 9 km grid frame. The calculation formula of VTCI is as follows: In the formula, They respectively represent that in the study area, when kNDVI i The maximum and minimum values ​​of LST when it is equal to a certain value; Indicates that the kNDVI value of a pixel is kNDVI i a, b, a', b' are unknown coefficients, which are approximated by the scatter plot of kNDVI and LST in the farmland area.

6. The method for predicting soil moisture in the root zone of farmland according to claim 5, characterized in that: The root zone soil moisture RZSM_t is given at time t, and the kNDVI_20m_Daily, SSM_1km_Daily and ERA5_20m_Daily at the same time and SSM_1km_Daily at time t+1 are input into the soil moisture prediction model to simulate the root zone soil moisture RZSM_t+1 at time t+1. The specific calculation formula is: AWF t+1 =f(AWF t ,P t ,SSM t+1 ,SSM t ,RH t ,PET t ,kNDVI t ); Among them, AWF t+1 represents the available water in the root zone at time t+1, AWF t represents the available water in the root zone at time t, P t represents the rainfall at time t, SSM t+1 represents SSM_1km_Daily at time t+1, SSM t represents SSM_1km_Daily at time t, RH t Relative humidity in ERA5_20m_Daily at time t, PET t represents the potential evapotranspiration at time t, KNDVI t represents kNDVI_20m_Daily at time t; RZSM t+1 =AWF t+1 ×(FC-WP)+WP; Among them, RZSM t+1 is the soil moisture in the root zone at time t+1, FC is the field water holding capacity, and WP is the wilting coefficient.

7. A farmland root zone soil moisture prediction system, characterized in that: include: Data acquisition module, used to obtain multi-source remote sensing data of farmland; The downscaling module is used to downscale the spatial resolution of MODIS multispectral images and Sentinel-2 multispectral images in multi-source remote sensing data to obtain the kernel normalized difference vegetation index kNDVI of the first spatial resolution; downscale the spatial resolution of MODIS LST land surface temperature and ERA5LST land surface temperature in multi-source remote sensing data to obtain the surface soil moisture SSM of the second spatial resolution; downscale the spatial resolution of ERA5-Land meteorological data in multi-source remote sensing data to obtain the meteorological data ERA5 of the third spatial resolution; The humidity prediction module is used to replace the features in the soil water balance equation with various remote sensing data to build a soil humidity prediction model; given the root zone soil humidity RZSM_t at time t, the kNDVI, SSM and ERA5 at the same time and SSM_1km_Daily at time t+1 are input into the soil humidity prediction model to simulate the root zone soil humidity RZSM_t+1 at time t+1; the current model output will be used as the model input for the next step, and this step will be repeated to complete the prediction of daily root zone soil humidity.

8. A computer device comprising a memory, a processor and a computer program stored in the memory, characterized in that: The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is loaded into a processor, it can execute the steps of the method according to any one of claims 1 to 6.

Images