A soil moisture downscaling method considering vegetation memory and fine-scale vegetation indices

Abstract

A soil moisture downscaling method considering vegetation memory and fine-scale vegetation indices, comprising: 1) generating fine-scale vegetation indices with a spatial resolution of 30 m; 2) obtaining NDVI data with hysteresis from soil moisture data and designing vegetation memory data; 3) obtaining MODIS data, SRTM data, SoilGrids data, ERA5-Land data and Landsat data of the target area, and performing spatio-temporal matching with SMAP data; 4) combining the vegetation memory data, the fine-scale vegetation indices and the data obtained in step 3) to construct a soil moisture downscaling model to obtain high-spatial-resolution soil moisture products. The high-spatial-resolution soil moisture data generated by the present invention has stronger spatial expression ability in the results, can reflect detailed information of soil moisture, and can effectively improve the application ability of soil moisture data in a small-scale range.

Description

Technical Field

[0001] The present invention relates to the field of soil moisture downscaling, and specifically to a soil moisture downscaling method considering vegetation memory and fine-scale vegetation indices. Background Art

[0002] Soil moisture plays an important role in processes such as water resource management, hydrological cycle, and drought development. As an important variable for the heat and water exchange between the earth's surface and the atmosphere, it is crucial for guiding agricultural production and ensuring food security. For regions with complex terrain, diverse land cover types, and heterogeneous soil types, there is an urgent need for high spatio-temporal resolution soil moisture data to improve the accuracy and ability of drought monitoring. However, current global soil moisture products generally have a low spatial resolution and cannot reflect the significant spatial heterogeneity of soil moisture caused by the combined effects of various land covers, different climate conditions, and complex terrain. Therefore, how to fully explore the complex relationships between soil moisture data and existing auxiliary data to obtain high spatio-temporal resolution soil moisture data is a very important issue.

[0003] There are many soil moisture products. Among them, SMAP (Soil Moisture Active and Passive) simultaneously uses L-band radar and radiometer to conduct parallel and synchronous measurements, combining the advantages of active (radar) and passive (radiometer) microwave remote sensing, and providing a requirement that can simultaneously meet high temporal resolution, wide spatial coverage, optimal sensing depth, and high-precision inversion of soil moisture under medium vegetation conditions, and is widely used in global soil moisture monitoring research. The SMAP project produces three soil moisture products, namely SM_P, SM_A, and SM_AP, with spatial resolutions of 36 km, 3 km, and 9 km respectively. Among them, due to the failure of the L-band radar in July 2015, currently only soil moisture data with spatial resolutions of 36 km and 9 km can be obtained. The lower spatial resolution remote sensing pixels contain heterogeneous land cover types and complex terrain, and cannot provide detailed changes in soil moisture, making it difficult to accurately reflect drought information at the local regional scale.

[0004] The spatial distribution of soil moisture interacts with many factors such as climate factors, topography, soil properties, and vegetation. Vegetation affects soil moisture by influencing solar incident radiation and soil evapotranspiration; surface temperature affects soil moisture by influencing soil water evapotranspiration, that is, the higher the surface temperature, the greater the soil water evapotranspiration, and the drier the surface appears; elevation affects soil moisture by influencing surface temperature, and elevation will also affect the magnitude of runoff and the redistribution of rainfall. The spatio-temporal variability of soil moisture distribution is the result of the combined action of multiple factors, and these factors should be applied to the soil moisture downscaling model simultaneously. At the same time, optical sensors are often affected by cloud contamination, resulting in a relatively small amount of available remote sensing image data in the study area. However, soil moisture downscaling requires the participation of multi-spectral data. How to fill in the missing spatial and temporal data of the available data in the study area is an issue that needs to be focused on. In addition, the existing spatial resolution of soil moisture is relatively low. How to provide more fine-scale information about soil moisture and improve the spatial resolution of soil moisture is also a problem that needs to be considered in downscaling research. On the other hand, there is a certain time lag between vegetation and soil moisture. The current state of soil moisture needs to be reflected on vegetation after a certain period of time, usually one month later than soil moisture. Existing soil moisture downscaling models rarely consider the impact of vegetation hysteresis on soil moisture. Summary of the Invention

[0005] Aiming at the problems in the prior art that rarely consider the hysteresis between soil moisture and vegetation, the connection between soil moisture and fine-scale vegetation information, the relationship between soil moisture data and multi-feature auxiliary data, and the problem that optical data cannot be effectively applied to the soil moisture downscaling model due to cloud and fog effects, the present invention proposes a combined multi-feature soil moisture downscaling method considering vegetation memory and fine-scale vegetation data, providing more fine-scale information and vegetation memory information about soil moisture for the soil moisture downscaling method, while filling in the spatial and temporal gaps of the available multi-spectral data in the study area, fully exploring the non-linear relationship between soil moisture and multi-feature data, and obtaining soil moisture data with a spatial resolution of 1 km.

[0006] To achieve the above object, according to one aspect of the present invention, a combined multi-feature soil moisture downscaling method is provided, including the following steps:

[0007] (1) Generate a fine-scale vegetation index with a spatial resolution of 30 m using the enhanced spatio-temporal adaptive fusion algorithm;

[0008] (2) Use the MOD13A3 monthly-scale NDVI to obtain NDVI data with hysteresis with respect to soil moisture data, and design vegetation memory data, one of the soil moisture downscaling auxiliary factors, to ensure that the vegetation memory data has a time lag of about one month with respect to the soil moisture data;

[0009] (3) Obtain SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data of the target area. Preprocess the MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data to have the same projection method, the same image coverage, and a spatial resolution consistent with the SMAP soil moisture, and then perform spatio-temporal matching with the SMAP data to serve as the input data for the soil moisture downscaling model;

[0010] (4) Combine the vegetation memory data, fine-scale vegetation index, SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data to construct a soil moisture downscaling model to obtain a high-spatial-resolution soil moisture product with a resolution of 1 km.

[0011] Further, step (1) includes:

[0012] (1.1) Prepare data: Provide two pairs of MODIS and Landsat reflectance images that are close to the prediction date and ensure that the two pairs of data have the same time and coverage, as well as a set of MODIS images for different time series to be predicted;

[0013] (1.2) Data preprocessing: Resample all MODIS data using ENVI software. The resampling method selects bilinear interpolation to reduce the influence of georeferencing errors. Georectify the resampled MODIS to obtain the same georegistration as Landsat. Then, perform cropping operations on all data to ensure that the images have exactly the same coverage area. All data are preprocessed to obtain surface reflectance. The preprocessed MODIS and Landsat images have the same projection, spatial resolution, and consistent coverage;

[0014] (1.3) Search for neighboring pixels: Use two high-resolution images to search for similar pixels of the central pixel in the local window of the prediction period image. The similar pixels are obtained by using the sliding window method. The adjacent pixel values need to meet the conditions of formula (1), that is, the standard deviation of the two pixels is small to ensure that pixels with high spectral similarity to the central pixel are obtained within the search window range:

[0015] |L(x i ,y i ,t k ,B)-L(x w / 2 ,y w / 2 ,t k ,B)|≤σ(B)·2 / m (10)

[0016] In the formula: L is the Landsat image, (x i , y i ) is the position of the i-th similar pixel; (x w / 2 , y w / 2 ) is the position of the central pixel at the prediction time, B represents the image band, t k represents the image time, L(x i , y i , t k , B) is the reflectance of the B band of the Landsat image at the pixel (x k , y i , y i ) at time t w / 2 , y w / 2 , t p , B) is the reflectance of the central pixel of the search window of the B band of the Landsat image at time t k ; σ(B) is the standard deviation of the reflectance value of band B; m represents the number of estimated categories;

[0017] (1.4) Calculate the weights of similar pixels: The formula for the spectral similarity between each similar pixel and its corresponding high-resolution and low-resolution pixels is shown in (2):

[0018]

[0019] Where:

[0020] L i ={L(x i , y i , t m , B1),..., L(x i , y i , t m , B n ), L(x i , y i , t n , B1),..., L(x i , y i , t n , B n}}(12)

[0021] M i ={M(x i , y i , t m , B1),..., M(x i , y i , t m , B n ), M(x i , yi , t n , B1), ..., M(x i , y i , t n , B n )}(13)

[0022] R i is the spectral correlation coefficient between the high - resolution pixel and the low - resolution pixel that describes the similar pixel i; L i and M i respectively represent the reflectance sets of similar pixels in each band of the high - spatial - resolution and low - spatial - resolution data during the time periods t m and t n ; E() represents the expected value; D(L i ) and D(M i ) are the variances of L i and M i respectively, and di is the geographical distance of the similar pixel i;

[0023] The geographical distance d between the i - th similar pixel and the central pixel i is as shown in formula (5):

[0024]

[0025] Combining the spectral similarity and geographical distance of pixel i, the index D i is calculated as:

[0026] D i = (1 - R i ) × d i (15)

[0027] According to the fact that the contribution of the similar pixel with a larger D i value to the calculation result of the central pixel is smaller, so the calculation formula of the similar - pixel weight W i is as shown in (7), that is, calculating the normalized reciprocal of D i , and the range of W i is 0 - 1, and the total weight of all similar pixels is 1:

[0028]

[0029] The time weight T k is calculated based on the difference in the reflectance of the MODIS image between time t k (k = m, n) and the prediction time t p , as shown in formula (10).

[0030]

[0031] (1.5) Calculate the conversion coefficient: The calculation formula for the conversion coefficient is shown in (9):

[0032]

[0033] v(x, y) is the ratio of the reflectance change between the high-resolution image and the low-resolution image, t m and t n are two different times, (x, y) represents the pixel position, B represents the band, L(x, y, t m , B) and L(x, y, t n , B) respectively represent the pixel values at the pixel position (x, y) of the B band of the Landsat data at t m and t n times, M(x, y, t m , B) and M(x, y, t n , B) respectively represent the pixel values at the pixel position (x, y) of the B band of the MODIS data at t m and t n times;

[0034] (1.6) Calculate the time weight: After calculating the similar pixel weight and the conversion coefficient, use the MODIS data at t m and t n two time points and the MODIS data at the prediction time t p and substitute them into formula (10) respectively to obtain the time weighting, and set a larger time weight for the fine-resolution reflectance to ensure that the fine-resolution data closer to the prediction date has similar reflectance values:

[0035]

[0036] In the formula: m and n represent different times; W represents the search window size, (x i , y i ) represents the pixel position, B represents the band of the image, M(x i , y i , t k , B) represents the MODIS data at the pixel position (x, y) of the B band at t k time, M(x i , y i ) represents the MODIS data at the pixel position (x, y) of the B band at the prediction time t i , y i , t p , B) represents the MODIS data at the pixel position (x, y) of the B band at the prediction time t p , T i , y i ) represents the time weight; k ;

[0037] (1.7) Calculate the central pixel value of the prediction period: According to the time weights T m and T n , for the high spatial resolution prediction period t p , the central pixel value calculation formula is as follows:

[0038] L(x w / 2 , y w / 2 , t p , B) = T m ×L m (x w / 2 , y w / 2 , t p , B)+T n ×L n (x w / 2 , y w / 2 , t p , B) (18)

[0039] In the formula: T m and T n are time weights; w represents the search window size, (x w / 2 , y w / 2 ) represents the central pixel, t p represents the prediction time, B represents the prediction band, L m (x w / 2 , y w / 2 , t p , B) and L n (x w / 2 , y w / 2 , t p , B) are the results of the prediction date calculated based on the central pixels of the B band of the high-resolution image at the reference days of t m and t n moments respectively;

[0040] (1.8) According to the obtained results of the prediction date, that is, the 30m spatial resolution multispectral data fused by the ESTARFM algorithm, obtain the fine-scale vegetation index NDVI fused , and the specific calculation formula is as follows:

[0041] NDVI fused =(NIR - R) / (NIR + R)

[0042] In the formula: NIR represents the near-infrared band of the 30m spatial resolution multispectral data fused, and R represents the red band of the 30m spatial resolution multispectral data fused.

[0043] Further, in step (2), the vegetation memory data with a time lag of about 1 month from the soil moisture data selects MODIS13A2, which is monthly-scale NDVI data. For the predicted SMAP soil moisture data in month b, the selected MODIS13A2 data is the data of month b + 1 in the same year, that is, NDVI lagged 。

[0044] Further, the SMAP data in step (3) is selected as the 9-km daily SPL3SMP_E data, specifically selecting soil moisture and brightness temperature data. The brightness temperature data includes horizontal polarization brightness temperature and vertical polarization brightness temperature;

[0045] MODIS data includes vegetation, temperature, evapotranspiration, primary productivity, surface albedo, land cover type. Vegetation data includes normalized difference vegetation index (NDVI, MOD13A3), enhanced vegetation index (EVI, MOD13A2), and leaf area index (LAI, MCD15A3). Temperature data is land surface temperature (LST, MOD11A1). Evapotranspiration data includes evapotranspiration (ET, MOD16A2) and potential evapotranspiration (PET, MOD16A2). Primary productivity data is gross primary productivity (GPP, MOD17A2). Surface albedo data is (Albedo, MCD43A3). Land cover type is (LandCover, MCD12Q1);

[0046] SRTM data includes elevation. Slope, aspect, and hillshade are calculated on the Google Earth Engine platform using the elevation data;

[0047] SoilGrids data selects the average soil sand content, average silt content, average clay content, and average soil pH value at a depth of 0 - 5 cm with a spatial resolution of 250 m;

[0048] ERA5-Land data downloads the ERA5-Land hourly dataset with a spatial resolution of 11132 m from the public data archive of Google Earth Engine (GEE), specifically obtaining the 0 - 7 cm surface soil moisture in the dataset;

[0049] Landsat data was obtained using the GEE cloud platform. The dataset is the Landsat8 OLI / TIRS sensor, including 6 visible and near-infrared bands and 2 shortwave infrared bands that have been orthorectified to surface reflectance, 2 thermal infrared bands and 1 panchromatic band that have been processed into orthorectified brightness temperature data. All data has been preprocessed for radiometric correction, terrain correction, and geometric correction on the GEE platform, and atmospheric correction has been performed using the LaSRC method, including masks for clouds, shadows, water, and snow generated using the CFMASK algorithm and a saturation mask for each pixel. The projection has been transformed to EPSG:32649, and the output format is GeoTIFF.

[0050] Furthermore, the preprocessing of MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data in step (3) specifically includes: performing projection transformation, cropping, and resampling on SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data on the GEE platform to ensure that the above data has the same projection method, the same image coverage, and the same spatial resolution as SMAP soil moisture.

[0051] Furthermore, the projection coordinate system used for projection transformation is the WGS84 geographic coordinate, the vector data of the study area is used for cropping, and the resampling method is the nearest neighbor method.

[0052] Furthermore, the spatio-temporal matching with SMAP data in step (3) to be used as input data for the soil moisture downscaling model specifically includes:

[0053] The temporal resolution of the horizontal polarization brightness temperature (TBh) and vertical polarization brightness temperature (TBv) of SMAP data is the same as that of SMAP soil moisture, and this data only undergoes spatial downscaling;

[0054] The temporal resolution of MOD13A3 is monthly. Considering that there is a lag of about 1 month between soil moisture and the vegetation index, the vegetation index of the next month after the month of the soil moisture data is selected as the vegetation lag data;

[0055] The temporal resolution of MCD15A3 is 4 days. The nearest data within the 4 days before or after the current day of the soil moisture data is selected, and the data within the 4 days after is preferred;

[0056] The temporal resolution of MOD13A2 is 16 days. The nearest data within the 16 days before or after the current day of the soil moisture data is selected, and the data within the 16 days after is preferred;

[0057] The time resolution of MCD15A3 is 4 days. Select the nearest data within the previous 4 days or the next 4 days of the current day for soil moisture, and give priority to the data of the next 4 days.

[0058] The time resolution of MOD11A1 is daily, which is consistent with the SMAP soil moisture. This data only undergoes spatial downscaling.

[0059] The time resolution of MOD16A2 is 8 days. Select the nearest data within the previous 8 days or the next 8 days of the current day for soil moisture, and give priority to the data of the next 8 days.

[0060] The time resolution of MOD17A2 is 8 days. Select the nearest data within the previous 8 days or the next 8 days of the current day for soil moisture, and give priority to the data of the next 8 days.

[0061] The time resolution of MCD43A3 is daily, which is consistent with the SMAP soil moisture. This data only undergoes spatial downscaling.

[0062] The SRTM data is terrain data. Considering the stability of the terrain, the same SRTM data is selected for all soil moisture downscaled data.

[0063] The SoilGrids data is soil texture data. Considering the stability of the terrain, the same SoilGrids data is selected for all soil moisture downscaled data.

[0064] The ERA5-Land data is coarse-scale soil moisture data with a time resolution of hours. Select the data at 9 am that is consistent with the SMAP soil moisture transit time.

[0065] The MODIS and Landsat data are processed by the ESTARFM spatio-temporal fusion algorithm to obtain 30-meter reflectance data. Calculate the NDVI together with the Landsat data to form a fine-scale vegetation index. The time resolution is less than 16 days. Select the nearest data within the previous 16 days or the next 16 days of the current day for soil moisture, and give priority to the data of the next 16 days.

[0066] Furthermore, the specific content of step (4) includes:

[0067] (4.1) Use the SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data that are spatio-temporally matched and have the same spatial resolution as the SMAP soil moisture as input data. Adopt the random forest as the training model for downscaling, and the output data is the soil moisture with low spatial resolution, thus obtaining the soil moisture downscaling model.

[0068] (4.2) Resample the original SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data after spatio-temporal matching to a spatial resolution of 1 km to obtain corresponding data with high spatial resolution;

[0069] (4.3) Input the resampled SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data with high spatial resolution into the soil moisture downscaling model to obtain soil moisture data with high spatial resolution.

[0070] Furthermore, the soil moisture downscaling model is as follows:

[0071]

[0072] where SM is soil moisture, RF is the soil moisture downscaling model based on random forest, NDVI fused is the NDVI data fused by the ESTARFM algorithm, NDVI lagged is the NDVI data with a time scale lag of 1 month from the time of soil moisture downscaling, EVI is the enhanced vegetation index, LAI is the leaf area index, LST is the land surface temperature data, ET is the evapotranspiration data, PET is the potential evapotranspiration data, GPP is the gross primary productivity, Albedo is the land surface albedo, LandCover is the land cover type, Elevation is the elevation data, Slope is the slope data, Aspect is the aspect data, HillShade is the hill shade data, Sand is the percentage of soil sand content, Silt is the percentage of soil silt content, Clay is the percentage of soil clay content, PH is the PH value in soil moisture, ERA5 selects the surface soil moisture data of 0 - 7 cm, Tbh is the SMAP horizontal polarization brightness temperature, and Tbv is the SMAP vertical polarization brightness temperature; where NDVI iused is the multi-spectral data with a spatial resolution of 30 m fused by the ESTARFM algorithm, and the spectral bands include the near-infrared band NIR and the red band R used to calculate NDVI, NDVI fused The specific calculation formula is as follows:

[0073] NDVI fused =(NIR - R) / (NIR + R).

[0074] Generally speaking, compared with the prior art through the above technical solutions conceived by the present invention, the following beneficial effects can be achieved:

[0075] 1. A soil moisture downscaling model that combines multi-source remote sensing data with multiple features is constructed to achieve spatial downscaling of low-spatial-resolution soil moisture data: Considering that the optical sensor is affected by cloud contamination, resulting in a reduction in the amount of available data, in order to provide more fine-scale information about soil moisture, the present invention generates fine-scale vegetation index NDVIfused data by means of the spatio-temporal fusion algorithm ESTARFM to make up for the spatial and temporal missing data in the study area. Considering that there is a certain time lag between vegetation and soil moisture, the state of current soil moisture can only be reflected on vegetation after a certain period of time, usually one month later than soil moisture. The present invention designs NDVIlagged data with a time lag of one month corresponding to the downscaled soil moisture data to solve the problem that the lag relationship between soil moisture and vegetation affects the accuracy of the soil moisture downscaling model. The present invention also uses vegetation data and temperature data commonly used in soil moisture downscaling. Surface vegetation affects soil moisture by influencing solar incident radiation and soil evapotranspiration. Surface temperature affects soil moisture by influencing soil water evapotranspiration, that is, the higher the surface temperature, the greater the soil water evapotranspiration, and the drier the surface appears. At the same time, the present invention also considers the influence of albedo, evapotranspiration, terrain, brightness temperature, soil texture, land use, and primary productivity factors on downscaling. In addition, the model incorporates coarse-scale soil moisture data, and low-spatial-resolution soil moisture can be used as the average soil moisture content at the coarse scale to promote the generation of high-spatial-resolution and reliable soil moisture.

[0076] 2. The present invention uses the downscaled high-spatial-resolution soil moisture data to achieve regional-scale soil moisture drought monitoring: Soil moisture is usually retrieved from P-band, L-band, C-band, and X-band data. There is less research on P-band soil moisture, and there are only a few studies based on remote sensing modeling and airborne flight observations. The L-band is more sensitive to soil moisture than the C-band and X-band, and L-band soil moisture data is more widely used. Commonly used L-band soil moisture data include SMAP and SMOS data, and the spatial resolutions of the two data are above 9 km and 25 km respectively. The low-spatial-resolution remote sensing pixels contain heterogeneous land cover types and complex terrains, which cannot provide detailed changes in soil moisture and are difficult to accurately reflect drought information at the local regional scale. Soil moisture plays a decisive role in the growth of crops and natural vegetation because the dynamics of soil moisture determine the available water resources in the agricultural ecosystem. The high-spatial-resolution soil moisture data obtained by the present invention shows in more detail the spatial distribution of soil moisture at the regional scale, can provide the soil moisture status on a daily continuous time scale, and provides refined soil moisture data for regional-scale drought monitoring, especially for complex underlying surfaces. Description of the Drawings

[0077] Figure 1 It is a flowchart of a soil moisture downscaling method considering vegetation memory and fine-scale vegetation index provided by an embodiment of the present invention;

[0078] Figure 2 It is a schematic diagram of the implementation process of an embodiment of the present invention;

[0079] Figure 3 It is a land use map of the target area in an embodiment of the present invention;

[0080] Figure 4 It is a result map of the fine-scale vegetation index after spatio-temporal fusion provided by an embodiment of the present invention;

[0081] Figure 5 It is a comparison map of the soil moisture in the study area before and after downscaling provided by an embodiment of the present invention;

[0082] Figure 6 It is a scatter plot of soil moisture before and after downscaling and measured sites. Detailed implementation manners

[0083] The present invention is applicable to downscaling low-spatial-resolution soil moisture data using multi-source and heterogeneous multi-feature remote sensing data, meeting the current demand that the spatial resolution of microwave soil moisture data is relatively low and its application in small-scale regional scales is insufficient. It is a general soil moisture downscaling method. In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0084] Soil moisture data, as one of the important variables for agricultural drought monitoring, plays a crucial role in drought monitoring. High-spatial-resolution soil moisture data has stronger spatial expression ability and can reflect detailed information of soil moisture. For study areas with complex terrain, diverse land cover types, and heterogeneous soil types, it can effectively improve the application ability of soil moisture data at a small scale. Accordingly, the present invention provides a soil moisture downscaling method considering vegetation memory and fine-scale vegetation indices. The enhanced spatio-temporal adaptive reflectance fusion model (ESTARFM) algorithm is used to generate fine-scale vegetation indices with a spatial resolution of 30 m. The monthly-scale NDVI of MOD13A3 is used to design vegetation memory data with a lag with respect to the soil moisture data. By combining MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data, a soil moisture calculation model with higher accuracy, richer spatial details, and more complex underlying surface is constructed to obtain consistent and reliable high-spatial-resolution (1 km) soil moisture downscaled data, providing important data support for agricultural drought monitoring and forecasting. The specific data used is shown in Table 1.

[0085] Table 1 Auxiliary factors for the downscaling method

[0086]

[0087] As Figure 1 and Figure 2 shown, the embodiments of the present invention provide a soil moisture downscaling method considering vegetation memory and fine-scale vegetation indices, which specifically includes the following steps:

[0088] Step 1: Use the enhanced spatio-temporal adaptive reflectance fusion model (ESTARFM) algorithm to generate fine-scale vegetation indices with a spatial resolution of 30 m;

[0089] Among them, according to the requirements of drought monitoring, in the embodiments of the present invention, the low-spatial-resolution reflectance data required by the ESTARFM algorithm is selected as MODIS11GA, and the high-spatial-resolution reflectance data is selected as Landsat:

[0090] In the embodiments of the present invention, the implementation manner of Step 1 is specifically as follows:

[0091] Step 1.1, Prepare data. The embodiments of the present invention provide two pairs of data close to the prediction date for the prediction date, and ensure that the two pairs of data have the same time and coverage range of MODIS and Landsat reflectance images, as well as a set of MODIS images for different time series to be predicted;

[0092] Step 1.2, Data preprocessing. All MODIS data are resampled using ENVI software to ensure the same spatial resolution as the Landsat images. The resampling method selects bilinear interpolation to reduce the influence of georeference errors. Geometric correction is performed on the resampled MODIS to obtain the same georegistration as Landsat. After that, cropping operations are performed on all data to ensure that the images have exactly the same coverage area. All data have been preprocessed to obtain surface reflectance. The preprocessed MODIS and Landsat images have the same projection, spatial resolution, and consistent coverage.

[0093] Step 1.3, Search for neighboring pixels. Two high-resolution images are used to search for similar pixels of the central pixel in the local window of the predicted image. The similar pixels are obtained by using the sliding window method. The adjacent pixel values need to meet the conditions of Equation (1), that is, the standard deviation of the two pixels is small to ensure that pixels with high spectral similarity characteristics to the central pixel are obtained within the search window range.

[0094] |L(x i ,y i ,t k ,B)-L(x w / 2 ,y w / 2 ,t k ,B)|≤σ(B)·2 / m (19)

[0095] In the formula: L is the Landsat image, (x i ,y i ) is the position of the i-th similar pixel; (x w / 2 ,y w / 2 ) is the position of the central pixel at the predicted time, B represents the image band, and t k represents the image time. L(x i ,y i ,t k ,B) is the reflectance of the Landsat image in band B at pixel (x k ,y i ,y i ) at time t w / 2 ,y w / 2 ,t p ,B) is the reflectance of the central pixel in the search window of the Landsat image in band B at time t k ; σ(B) is the standard deviation of the reflectance values of band B; m represents the number of estimated categories.

[0096] Step 1.4, calculate the weights of similar pixels. The more similar a similar pixel is to the predicted pixel, the greater the weight value of the similar pixel. Conversely, the less similar the two are, the smaller the weight value of the similar pixel. The spectral similarity calculation formula between each similar pixel and its corresponding high-resolution and low-resolution pixels is shown in (2):

[0097]

[0098] L i ={L(x i , y i , t m , B1),..., L(x i , y i , t m , B n ), L(x i , y i , t n , B1),..., L(x i , y i , t n , B n )}(21)

[0099]

[0100] R i is the spectral correlation coefficient between the high-resolution pixel and the low-resolution pixel describing the similar pixel i; L i and M i respectively represent the reflectance sets of similar pixels in each band of high-spatial-resolution and low-spatial-resolution data during the time periods t m and t n ; E() represents the expected value; D(L i ) and D(M i ) are the variances of L i and M i respectively, and di is the geographical distance of the similar pixel i.

[0101] The geographical distance d i between the i-th similar pixel and the central pixel is shown in formula (5):

[0102]

[0103] Combining the spectral similarity and geographical distance of pixel i, the index D i can be calculated as follows:

[0104] D i =(1 - R i )×d i (24)

[0105] According to D i Similar pixels with larger D values contribute less to the calculation result of the central pixel. Therefore, the weight W i is calculated as shown in formula (7), that is, calculate the i normalized reciprocal of D. According to the formula, it can be known that the range of W i is 0 - 1, and the total weight of all similar pixels is 1.

[0106]

[0107] The time weight T k is calculated based on the difference in the MODIS image reflectance at time t k (k = m, n) and the predicted time t p as shown in formula (10).

[0108]

[0109] Step 1.5: Calculate the conversion coefficient. The conversion coefficient can be obtained by performing linear regression analysis on each similar pixel in the high - resolution and low - resolution data within the moving window. Considering that the pixel geometric positions in different data sources cannot be ensured to coincide through pre - processing operations, calculating the conversion coefficient for only a single similar pixel will introduce large errors in the results. The ESTARFM algorithm calculates the conversion coefficient by searching for similar pixels within a certain range of the central pixel. The calculation formula of the conversion coefficient is as shown in (9).

[0110]

[0111] v(x, y) is the ratio of the reflectance change between the high - resolution image and the low - resolution image. The algorithm assumes that the pixel reflectance of MODIS and Landsat remote sensing data has a linear change characteristic within different time periods. Therefore, the value of v(x, y) remains constant.

[0112] Step 1.6: Calculate the central pixel value for the prediction period. According to the time weights T m and T n , the calculation formula for the central pixel value of the high - spatial - resolution prediction period t p is as follows:

[0113] L(x w / 2 , y w / 2 , t p , B) = T m ×L m (x w / 2 , y w / 2 , t p , B) + T n ×L n (x w / 2, y w / 2 , t p , B) (28)

[0114] Where: T m and T n are time weights; L m (x w / 2 , y w / 2 , t p , B) and L n (x w / 2 , y w / 2 , t p , B) are the results of the predicted date calculated based on the high-resolution image at time t m and t n as the reference date.

[0115] Step 1.7, according to the obtained results of the predicted date, that is, the 30m spatial resolution multispectral data fused by the ESTARFM algorithm, obtain the fine-scale vegetation index NDVI fused , the fine-scale vegetation index results after spatio-temporal fusion obtained by the present invention are as Figure 4 shown, and the specific calculation formula is as follows:

[0116] NDVI iused =(NIR - R) / (NIR + R)

[0117] Where: NIR represents the near-infrared band after fusion, and R represents the red band after fusion.

[0118] Step 2, use the MOD13A3 monthly-scale NDVI to obtain NDVI data with a lag with respect to the soil moisture data, and design the vegetation memory data, one of the auxiliary factors for downscaling soil moisture, to ensure that the data lags behind the soil moisture data by about 1 month. The specific steps are as follows;

[0119] The vegetation data with a lag of about 1 month with respect to the soil moisture data selects MODIS13A2, which is monthly-scale NDVI data. For the predicted SMAP soil moisture data in month b, the selected MODIS13A2 data is the data in month b + 1 of the same year, that is, NDVI lagged .

[0120] Step 3, obtain the MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data of the target area (as Figure 3 shown), and perform spatio-temporal matching with the SMAP data. The specific steps are as follows:

[0121] Step 3.1, The SMAP data includes three types of products, namely SM_P, SM_A, and SM_AP, with spatial resolutions of 36 km, 3 km, and 9 km respectively. In the embodiment of the present invention, the SPL3SMP_E data with a spatial resolution of 9 km and a temporal resolution of daily is selected for downscaling. Specifically, soil moisture and brightness temperature data are selected. The brightness temperature data includes horizontal polarization brightness temperature (TBh) and vertical polarization brightness temperature (TBv).

[0122] Step 3.2, For the MODIS data, the normalized difference vegetation index (NDVI, MOD13A3), enhanced vegetation index (EVI, MOD13A2), and leaf area index (LAI, MCD15A3) are specifically selected. The temperature data is the land surface temperature (LST, MOD11A1). The evapotranspiration data includes evapotranspiration (ET, MOD16A2) and potential evapotranspiration (PET, MOD16A2). The primary productivity data is the gross primary productivity (GPP, MOD17A2). The land surface albedo data is (Albedo, MCD43A3). The land cover type is (LandCover, MCD12Q1).

[0123] Step 3.3, For the SRTM data, the elevation (DEM) is selected, and the slope, aspect, and hill shade are calculated using the DEM data on the Google Earth Engine (GEE) platform.

[0124] Step 3.4, For the SoilGrids data, the average soil sand content, average silt content, average clay content, and average soil pH value at a depth of 0 - 5 cm with a spatial resolution of 250 m are selected.

[0125] Step 3.5, For the ERA5-Land data, the surface soil moisture (0 - 7 cm) at 9 am with a spatial resolution of 11132 m is selected.

[0126] Step 3.6, The Landsat data is obtained using the GEE cloud platform. The dataset is the Landsat8 OLI / TIRS sensor, which includes 6 visible and near-infrared bands processed as surface reflectance after orthorectification, 2 shortwave infrared bands, 2 thermal infrared bands processed as orthorectified brightness temperature data, and 1 panchromatic band. All the data has been preprocessed on the GEE platform, including radiometric correction, topographic correction, geometric correction, etc., and atmospheric correction has been performed using the LaSRC method, including the masks of clouds, shadows, water, and snow generated using the CFMASK algorithm and the saturation mask for each pixel. The projection is transformed to EPSG:32649, and the output format is GeoTIFF.

[0127] Step 3.7, perform projection transformation, cropping, and resampling on SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data on the GEE platform to ensure that the above-mentioned data have the same projection method, the same image coverage, and the same spatial resolution as SMAP soil moisture;

[0128] Specifically, the projection coordinate system is WGS84 geographic coordinates, and the vector data of the study area is used for cropping; the resampling method is the nearest neighbor method;

[0129] Step 3.8, the time resolution of SMAP TBh and TBv is the same as that of SMAP soil moisture, and this data only undergoes spatial downscaling;

[0130] Step 3.9, the time resolution of MOD13A3 is monthly. Considering that there is a lag of about 1 month between soil moisture and the vegetation index, select the vegetation index of the next month after the month of the soil moisture data as the vegetation lag data;

[0131] Step 3.10, the time resolution of MCD15A3 is 4 days. Select the nearest data within the 4 days before or after the current day of the soil moisture, and give priority to the data of the 4 days after;

[0132] Step 3.11, the time resolution of MOD13A2 is 16 days. Select the nearest data within the 16 days before or after the current day of the soil moisture, and give priority to the data of the 16 days after;

[0133] Step 3.12, the time resolution of MCD15A3 is 4 days. Select the nearest data within the 4 days before or after the current day of the soil moisture, and give priority to the data of the 4 days after;

[0134] Step 3.13, the time resolution of MOD11A1 is daily, which is the same as that of SMAP soil moisture, and this data only undergoes spatial downscaling;

[0135] Step 3.14, the time resolution of MOD16A2 is 8 days. Select the nearest data within the 8 days before or after the current day of the soil moisture, and give priority to the data of the 8 days after;

[0136] Step 3.15, the time resolution of MOD17A2 is 8 days. Select the nearest data within the 8 days before or after the current day of the soil moisture, and give priority to the data of the 8 days after;

[0137] Step 3.16, the time resolution of MCD43A3 is daily, which is the same as that of SMAP soil moisture, and this data only undergoes spatial downscaling;

[0138] Step 3.17: The SRTM data is topographic data. Considering the stability of the terrain, the same SRTM data is selected for all soil moisture downscaling data.

[0139] Step 3.18: The SoilGrids data is soil texture data. Considering the stability of the terrain, the same SoilGrids data is selected for all soil moisture downscaling data.

[0140] Step 3.19: The ERA5-Land data is the coarse-scale soil moisture data with a time resolution of hours. The data at 9:00 am that is consistent with the SMAP soil moisture transit time is selected.

[0141] Step 3.20: The MODIS and Landsat data are processed by the ESTARFM spatio-temporal fusion algorithm to obtain 30-meter reflectance data. The NDVI is calculated together with the Landsat data to form a fine-scale vegetation index with a time resolution of less than 16 days. The data closest to the current day within the previous 16 days or the next 16 days of soil moisture is selected, and the data of the next 16 days is preferred.

[0142] Step 4: Combine the vegetation memory data, fine-scale vegetation index, MODIS data, SRTM data, SoilGrids data, and ERA5-Land data to construct a soil moisture downscaling model to obtain high-spatial-resolution soil moisture products. The specific steps are as follows:

[0143] Step 4.1: Use the random forest as the training model for downscaling. The spatio-temporally matched SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data with a 9-km spatial resolution are used as input data, and the soil moisture with a 9-km low spatial resolution is used as the output data to obtain the soil moisture downscaling model.

[0144] Step 4.2: Project and clip the original SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data to ensure that the above data have the same projection method and the same image coverage, and resample them to obtain the corresponding data with a high spatial resolution (1 km).

[0145] Step 4.3: Match the resampled SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data in time according to Steps 3.8 - 3.20.

[0146] Step 4.4, input the spatially and temporally matched SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data with a spatial resolution of 1 km in height into the soil moisture downscaling model to obtain soil moisture data with a spatial resolution of 1 km in height. The soil moisture downscaling model is as follows:

[0147]

[0148] where SM is soil moisture, RF is the soil moisture downscaling model based on random forest, NDVI fused is the NDVI data fused by the ESTARFM algorithm, NDVI lagged is the NDVI data with a time scale lag of 1 month from the soil moisture downscaling time, EVI is the enhanced vegetation index, LAI is the leaf area index, LST is the land surface temperature data, ET is the evapotranspiration data, PET is the potential evapotranspiration data, GPP is the gross primary productivity, Albedo is the land surface albedo, LandCover is the land cover type, Elevation is the elevation data, Slope is the slope data, Aspect is the aspect data, HillShade is the hill shade data, Sand is the percentage of soil sand content, Silt is the percentage of soil silt content, Clay is the percentage of soil clay content, PH is the PH value in soil moisture, ERA5 selects the surface soil moisture (0 - 7 cm) data, Tbh is the SMAP horizontal polarization brightness temperature, and Tbv is the SMAP vertical polarization brightness temperature; where NDVI fused is the multi-spectral data with a spatial resolution of 30 m fused by the ESTARFM algorithm. The spectral bands include the near-infrared band (NIR) and the red band (R) used to calculate NDVI. The specific calculation formula of NDVI is as follows:

[0149] NDVI = (NIR - R) / (NIR + R)

[0150] Thus, the soil moisture with high precision and high spatial resolution in the study area range is obtained. The comparison of the soil moisture in the study area before and after downscaling provided by the present invention is Figure 5 as shown. The data before downscaling has a low spatial resolution, while the soil moisture data after downscaling can contain more detailed information. Figure 6 It shows that the correlation between the soil moisture after downscaling and the measured sites is above 0.4 m3 / m3, while the correlation between the soil moisture before downscaling and the measured sites is below 0.4 m3 / m3. There is a significant improvement in the correlation, and the deviation between the soil moisture after downscaling and the measured sites is small. The model shows significantly higher performance for basically all sites, and the performance of the downscaling result is better.

[0151] As described above, it is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed by the present invention should be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. A soil moisture downscaling method considering vegetation memory and fine-scale vegetation indices, characterized in that: It includes the following steps: (1) Generate a fine-scale vegetation index with a 30m spatial resolution using an enhanced spatio-temporal adaptive fusion algorithm; (2) Use the monthly-scale NDVI of MOD13A3 to obtain NDVI data with a lag with respect to soil moisture data, and design vegetation memory data, one of the auxiliary factors for soil moisture downscaling, to ensure that the vegetation memory data lags behind the soil moisture data by about 1 month; (3) Obtain SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data for the target area. Preprocess the MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data to have the same projection method, the same image coverage, and the same spatial resolution as the SMAP soil moisture. Then perform spatio-temporal matching with the SMAP data to be used as the input data for the soil moisture downscaling model; (4) Combine the vegetation memory data, fine-scale vegetation index, SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data to construct a soil moisture downscaling model to obtain a soil moisture product with a high spatial resolution of 1km; Step (1) includes: (1.1) Prepare data: Provide two pairs of MODIS and Landsat reflectance images that are close to the prediction date and ensure that the two pairs of data have the same time and coverage, as well as a set of MODIS images for different time series to be predicted; (1.2) Data preprocessing: Resample all MODIS data using ENVI software. The resampling method selects bilinear interpolation to reduce the influence of georeferencing errors. Georectify the resampled MODIS to obtain the same georegistration as Landsat. Then perform a cropping operation on all data to ensure that the images have exactly the same coverage area. All data are preprocessed to obtain surface reflectance. The preprocessed MODIS and Landsat images have the same projection, spatial resolution, and consistent coverage; (1.3) Search for neighboring pixels: Use two high-resolution images to search for similar pixels of the central pixel in the local window of the prediction period image. The similar pixels are obtained using a sliding window method. The adjacent pixel values need to meet the conditions of Equation (1), that is, the standard deviation of the two pixels is small, to ensure that pixels with high spectral similarity to the central pixel are obtained within the search window range; |L(x i ,y i ,t k ,B)-L(x w / 2 ,y w / 2 ,t k ,B)| ≤ σ(B)·2 / m (1) Where: L is the Landsat image, (x i , y i ) is the position of the i-th similar pixel; (x w / 2 , y w / 2 ) is the position of the central pixel at the prediction time, B represents the image band, t k represents the image time, L(x i , y i , t k , B) is the reflectance of the B band of the Landsat image at the pixel (x k , y i , y i ) at time t w / 2 , y w / 2 , t p , B) is the reflectance of the central pixel of the search window of the B band of the Landsat image at time t k ; σ(B) is the standard deviation of the reflectance value of band B; m represents the number of estimated classes; (1.4) Calculate the weights of similar pixels: The spectral similarity calculation formula between each similar pixel and its corresponding high-resolution and low-resolution pixels is shown in Equation (2): Where: L i = {L(x i , y i , t m , B1),..., L(x i , y i , t m , B n ), L(x i , y i , t n , B1),..., L(x i , y i , t n , B n )} (3) M i = {M(x i , y i , t m , B1),..., M(x i , y i , t m , B n ), M(x i , y i , t n , B1),..., M(x i , y i , t n , B n )}(4) R i is the spectral correlation coefficient between the high - resolution pixel and the low - resolution pixel that describe the similar pixel i; L i and M i respectively represent the reflectance sets of the similar pixels in each band of the high - spatial - resolution and low - spatial - resolution data during the time periods t m and t n ; E() represents the expected value; D(L i ) and D(M i ) are the variances of L i and M i respectively, and di is the geographical distance of the similar pixel i; The geographical distance d between the i-th similar pixel and the central pixel i As shown in formula (5): Calculate the index D by combining the spectral similarity and geographical distance of pixel i i : D i = (1 - R i ) × d i (6) According to D i Similar pixels with larger values contribute less to the calculation result of the central pixel. Therefore, the weight W of similar pixels i The calculation formula is as shown in (7), that is, calculate the normalized reciprocal of D i The range of W i is 0-1, and the total weight of all similar pixels is 1: Time weight T k Based on the reflectance of the MODIS image at time t k (k = m, n) and the predicted time t p The calculation is based on the difference between them, as shown in Equation (8): (1.5) Calculate the conversion coefficient: The calculation formula of the conversion coefficient is shown in Equation (9): v(x,y) is the ratio of the reflectance change between the high-resolution image and the low-resolution image, t m and t n are two different times, (x,y) represents the pixel position, B represents the band, L(x,y,t m ,B) and L(x,y,t n ,B) respectively represent the pixel values at the pixel position (x,y) in the B band of the Landsat data at time t m and t n ; M(x,y,t m ,B) and M(x,y,t n ,V) respectively represent the pixel values at the pixel position (x,y) in the B band of the MODIS data at time t m and t n ; (1.6) Calculate the time weight: After calculating the similarity pixel weight and the conversion coefficient, use the MODIS data at two time points t m and t n and the MODIS data at the predicted time t p and substitute them into Equation (10) respectively to obtain the time weighting. Set a larger time weight for the fine-resolution reflectance to ensure that the fine-resolution data closer to the predicted date have similar reflectance values: m and t n The MODIS data at two time points and the MODIS data at the predicted time t p are respectively substituted into Equation (10) to obtain the time weighting. Set a larger time weight for the fine-resolution reflectance to ensure that the fine-resolution data closer to the predicted date have similar reflectance values: Where: m and n represent different times; W represents the search window size, (x i , y i ) represents the pixel position, B represents the band of the image, M(x i , y i , t k , B) represents the MODIS data of the pixel position (x k , y i , y i ) at time t and band B, M(x i , y i , t p , B) represents the predicted MODIS data of the pixel position (x p , y i , y i ) at time t and band B, T k represents the time weight; (1.7) Calculate the central pixel value in the prediction period: According to the time weights T m and T n , for the high spatial resolution prediction period t p , the central pixel value calculation formula is as follows: L(x w / 2 ,y w / 2 ,t p ,B) = T m ×L m (x w / 2 ,y w / 2 ,t p ,B) + T n ×L n (x w / 2 ,y w / 2 ,t p ,B) (11) Where: T m and T n are time weights; w represents the search window size, (x w / 2 , y w / 2 ) represents the central pixel, t p represents the prediction time, B represents the prediction band, L m (x w / 2 , y w / 2 , t p , B) and L n (x w / 2 , y w / 2 , t p , B) are the results of the prediction dates calculated based on the central pixels of the B band of the high-resolution image at the reference days of t m and t n moments, respectively; (1.8)According to the result of the obtained prediction date, that is, the multispectral data with a 30m spatial resolution fused by the ESTARFM algorithm, obtain the fine-scale vegetation index NDVI fused , and the specific calculation formula is as follows: NDVI fused =(NIR - R) / (NIR + R) In the formula: NIR represents the near-infrared band of the 30m spatial resolution multispectral data obtained by fusion, and R represents the red band of the 30m spatial resolution multispectral data obtained by fusion.

2. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to claim 1, characterized in that: In step (2), the vegetation memory data with a time lag of about 1 month from the soil moisture data selects MODIS13A2, which is monthly-scale NDVI data. For the SMAP soil moisture data predicted in month b, the selected MODIS13A2 data is the data of month b + 1 of the same year, that is, NDVI lagged .

3. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to claim 1, characterized in that: The SMAP data in step (3) is selected as the 9-km daily SPL3SMP_E data. Specifically, soil moisture and brightness temperature data are selected. The brightness temperature data includes horizontal polarization brightness temperature and vertical polarization brightness temperature; The MODIS data includes vegetation, temperature, evapotranspiration, primary productivity, surface albedo, and land cover type. The vegetation data includes the Normalized Difference Vegetation Index (NDVI, MOD13A3), Enhanced Vegetation Index (EVI, MOD13A2), and Leaf Area Index (LAI, MCD15A3). The temperature data is the Land Surface Temperature (LST, MOD11A1). The evapotranspiration data includes evapotranspiration (ET, MOD16A2) and potential evapotranspiration (PET, MOD16A2). The primary productivity data is the Gross Primary Productivity (GPP, MOD17A2). The surface albedo data is (Albedo, MCD43A3). The land cover type is (LandCover, MCD12Q1); The SRTM data includes elevation. Slope, aspect, and hillshade are calculated from the elevation data on the Google Earth Engine platform; The SoilGrids data selects the average soil sand content, average silt content, average clay content, and average soil pH value at a 250-meter spatial resolution for a depth of 0-5 cm; The ERA5-Land data downloads the ERA5-Land hourly dataset with a spatial resolution of 11132 meters from the public data archive of Google Earth Engine (GEE). Specifically, the 0-7 cm surface soil moisture in the dataset is obtained; The Landsat data is obtained using the GEE cloud platform. The dataset is the Landsat8 OLI / TIRS sensor, including 6 visible and near-infrared bands processed to surface reflectance, 2 shortwave infrared bands, 2 thermal infrared bands processed to orthorectified brightness temperature data, and 1 panchromatic band. All data has been preprocessed for radiometric correction, topographic correction, and geometric correction on the GEE platform, and atmospheric correction has been performed using the LaSRC method, including cloud, shadow, water, and snow masks generated using the CFMASK algorithm and a saturation mask for each pixel. The projection is transformed to EPSG:32649, and the output format is GeoTIFF.

4. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to claim 3, characterized in that: The preprocessing of the MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data in step (3) specifically includes: performing projection conversion, cropping, and resampling on the SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data on the GEE platform to ensure that the above data has the same projection method, the same image coverage, and the same spatial resolution as the SMAP soil moisture.

5. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to claim 4, characterized in that: Among them, the projection coordinate system used for projection conversion is the WGS84 geographic coordinate. The vector data of the study area is used for cropping, and the resampling method is the nearest neighbor method.

6. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to claim 3, characterized in that: In step (3), spatio-temporal matching is performed with SMAP data as the input data for the soil moisture downscaling model, which specifically includes: The temporal resolution of the horizontal polarization brightness temperature (TBh) and vertical polarization brightness temperature (TBv) of SMAP data is the same as that of SMAP soil moisture. This data only undergoes spatial downscaling; The temporal resolution of MOD13A3 is monthly. Considering the lag of about 1 month between soil moisture and vegetation index, the vegetation index of the next month after the month of the soil moisture data is selected as the vegetation lag data; The temporal resolution of MCD15A3 is 4 days. The nearest data within the 4 days before or after the current day of the soil moisture is selected, and the data of the 4 days after is preferred; The temporal resolution of MOD13A2 is 16 days. The nearest data within the 16 days before or after the current day of the soil moisture is selected, and the data of the 16 days after is preferred; The temporal resolution of MCD15A3 is 4 days. The nearest data within the 4 days before or after the current day of the soil moisture is selected, and the data of the 4 days after is preferred; The temporal resolution of MOD11A1 is daily, which is the same as that of SMAP soil moisture. This data only undergoes spatial downscaling; The temporal resolution of MOD16A2 is 8 days. The nearest data within the 8 days before or after the current day of the soil moisture is selected, and the data of the 8 days after is preferred; The temporal resolution of MOD17A2 is 8 days. The nearest data within the 8 days before or after the current day of the soil moisture is selected, and the data of the 8 days after is preferred; The temporal resolution of MCD43A3 is daily, which is the same as that of SMAP soil moisture. This data only undergoes spatial downscaling; The SRTM data is terrain data. Considering the stability of the terrain, the same SRTM data is selected for all soil moisture downscaling data; The SoilGrids data is soil texture data. Considering the stability of the terrain, the same SoilGrids data is selected for all soil moisture downscaling data; The ERA5-Land data is coarse-scale soil moisture data with a temporal resolution of hours. The data at 9 am consistent with the SMAP soil moisture transit time is selected; The MODIS and Landsat data undergo the ESTARFM spatio-temporal fusion algorithm to obtain 30-meter reflectance data. Together with the Landsat data, the NDVI is calculated to form a fine-scale vegetation index with a temporal resolution of less than 16 days. The nearest data within the 16 days before or after the current day of the soil moisture is selected, and the data of the 16 days after is preferred.

7. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to claim 1, characterized in that: Step (4) specifically includes: (4.1) Using the SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data that have undergone spatio-temporal matching and have the same spatial resolution as the SMAP soil moisture as input data, and adopting the random forest as the training model for downscaling. The output data is low-spatial-resolution soil moisture, and a soil moisture downscaling model is obtained; (4.2) Resample the original SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data after spatio-temporal matching to a spatial resolution of 1 km to obtain the corresponding data with high spatial resolution; (4.3) Input the resampled SMAP data, MODIS data, SRTM data, SoilGrids data, ERA5-Land data, and Landsat data with high spatial resolution into the soil moisture downscaling model to obtain soil moisture data with high spatial resolution.

8. The soil moisture downscaling method considering vegetation memory and fine-scale vegetation index according to any one of claims 1-7, characterized in that: The soil moisture downscaling model is: Among them, SM is soil moisture, RF is the soil moisture downscaling model based on random forest, NDVI fused is the NDVI data fused by the ESTARFM algorithm, NDVI lagged is the NDVI data with a 1-month time lag from the soil moisture downscaling time scale. EVI is the enhanced vegetation index, LAI is the leaf area index, LST is the land surface temperature data, ET is the evapotranspiration data, PET is the potential evapotranspiration data, GPP is the gross primary productivity, Albedo is the land surface albedo, LandCover is the land cover type, Elevation is the elevation data, Slope is the slope data, Aspect is the aspect data, HillShade is the hillshade data, Sand is the percentage of soil sand content, Silt is the percentage of soil silt content, Clay is the percentage of soil clay content, PH is the PH value in soil moisture, ERA5 selects the surface soil moisture data of 0-7 cm, Tbh is the SMAP horizontal polarization brightness temperature, and Tbv is the SMAP vertical polarization brightness temperature; among them, NDVI fused is the multispectral data with a 30m spatial resolution fused by the ESTARFM algorithm. The spectral bands include the near-infrared band NIR and the red band R used to calculate NDVI, NDVI fused The specific calculation formula is as follows: NDVI fused =(NIR - R) / (NIR + R).

Images