Mountain area GNSS-R soil humidity inversion method fusing terrain humidity index

Abstract

The invention discloses a mountain area GNSS-R soil humidity inversion method fused with a terrain humidity index. The method comprises the following steps: acquiring remote sensing data of other areas except a water body area in a target mountain area and soil humidity data of each observation station; the altitude, the gradient, the slope direction, the longitude and latitude, the vegetation index, the clay content and the signal-to-noise ratio are selected as influence factors of mountainous area soil humidity inversion; calculating a corrected signal-to-noise ratio and a terrain humidity index; and constructing a mountainous area soil humidity inversion model, training and verifying, inputting the test set into the trained model, obtaining 1km-resolution soil humidity prediction day data, and averaging monthly soil humidity prediction day data to obtain soil humidity prediction monthly data. According to the method, the terrain humidity index and the CYGNSS data are fused, the advantages of terrain factors and the remote sensing technology can be fully utilized, the precision of soil humidity distribution and the prediction capacity of temporal and spatial variation are further improved, and the dynamic change of soil humidity in complex terrain areas such as mountainous areas is better revealed.

Description

Technical Field

[0001] The present invention relates to a method for retrieving soil moisture in mountainous areas by fusing terrain humidity index, and belongs to the technical field of surface soil moisture retrieval. Background Art

[0002] Soil moisture, as a key variable describing the water and energy exchange between the earth's surface and the atmosphere, is one of the key parameters affecting hydrological processes, ecological processes, and biogeochemical processes, and is a necessary influencing factor in research fields such as climate prediction, water resources management, agricultural meteorology, and flood disasters. Traditional soil moisture monitoring methods have spatio-temporal limitations, and the acquisition process and data are relatively cumbersome and complex. In contrast, professional remote sensing satellite methods have better comprehensiveness, timeliness, and regional continuity. With the progress of satellite navigation technology, the Global Navigation Satellite System (GNSS), represented by the Global Positioning System (GPS) of the United States, the GLONASS system of Russia, the GALILEO system of the European Union, and the Beidou satellite navigation system of China, has provided more than a hundred satellite resources for society. On this basis, the spaceborne Global Navigation Satellite System Reflectometry (GNSS-R) technology, which obtains the state of the target object by receiving and processing the surface reflection signal, has emerged. As a newly emerging microwave remote sensing technology in recent years, spaceborne GNSS-R has the characteristics of passive detection, low power consumption, low deployment cost, and wide range by setting the receiving device on a low-earth orbit satellite. The revisit period can be as low as several hours, and the time resolution is high. At the same time, the GNSS signal belongs to the L band, has good penetration, and can achieve all-weather and all-day observation. Among them, the Cyclone Global Navigation Satellite System (CYGNSS) launched by the United States in 2016, as a member of GNSS-R, is the most popular in the research field of soil moisture.

[0003] Currently, there are still many problems and challenges in the algorithm for retrieving soil moisture using CYGNSS data in practical applications: (1) Mountainous areas are mostly rain-fed agricultural areas, and water resource transportation is difficult. The development of mountainous agriculture urgently needs timely and accurate soil moisture data. The effective acquisition of soil moisture in mountainous areas will provide a more reliable scientific decision-making basis for fields such as water resources management, drought warning, and agricultural irrigation. However, the existing GNSS-R soil moisture retrieval algorithms are mainly designed for large-scale areas, lacking optimization for specific terrain conditions in mountainous areas, resulting in a weakened application effect of the constructed algorithm model in mountainous areas compared to the overall situation. Most of the current algorithm studies only consider conventional terrain elements such as elevation and slope, and simply use them as input features of the machine learning model, which cannot well solve the problem of GNSS-R soil moisture retrieval in mountainous areas.

[0004] (2)In the field of soil moisture research, the Topographic Wetness Index (TWI) is widely used to characterize the spatial variation characteristics of soil moisture distribution in mountainous areas. TWI takes into account topographic factors such as slope, aspect, and the hydrological processes of the watershed, so it can effectively reflect the influence of topography on soil moisture. Especially in mountainous areas where the terrain is complex and the spatial distribution of soil moisture often has significant non-uniformity, TWI provides an important quantitative tool for analyzing and predicting soil moisture in this area. However, current research only simply takes the topographic wetness index as an influencing factor input for machine learning, and does not truly achieve the integration with GNSS-R data. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a method for inverting soil moisture of GNSS-R in mountainous areas by integrating the topographic wetness index, which effectively improves the inversion accuracy of spaceborne GNSS-R in mountainous areas.

[0006] The present invention adopts the following technical solutions to solve the above technical problems: A method for inverting soil moisture of GNSS-R in mountainous areas by integrating the topographic wetness index, comprising the following steps: Step 1: Obtain remote sensing data of other areas except water body areas in the target mountainous area, including Cyclone Global Navigation Satellite System (GNSS) received data, Digital Elevation Model (DEM) data, Moderate Resolution Imaging Spectroradiometer (MODIS) observation data, and SoilGrids global soil information data; at the same time, obtain soil moisture data of each observation station in the target mountainous area; preprocess the obtained remote sensing data and soil moisture data. Step 2: Select altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and signal-to-noise ratio as influencing factors for inverting soil moisture of GNSS-R in mountainous areas, obtain altitude, slope, aspect, and longitude and latitude data from the preprocessed DEM data, calculate the vegetation index using the preprocessed MODIS observation data, and obtain clay content data from the preprocessed SoilGrids global soil information data. Step 3: Correct the signal-to-noise ratio in the preprocessed Cyclone GNSS received data, and calculate the topographic wetness index using the preprocessed DEM data; spatiotemporally match altitude, slope, aspect, longitude and latitude, vegetation index, clay content, the corrected signal-to-noise ratio, the topographic wetness index, and the preprocessed soil moisture data of the observation station into daily data under a 1 km resolution grid. Step 4: Obtain the daily data of signal-to-noise ratio (SNR) and terrain humidity index corresponding to the 1-km resolution grid with non-null SNR values. Use the SNR as the dependent variable and the terrain humidity index as the independent variable to establish a daily corresponding spatial fitting function model. Utilize the daily corresponding spatial fitting function model to spatialize the SNR of the 1-km resolution grid with null SNR values on the same day, and obtain the SNR of all 1-km resolution grids daily. Step 5: Take the altitude, slope, aspect, longitude and latitude, vegetation index, clay content, corrected SNR, and preprocessed daily soil humidity data of the observation station corresponding to all 1-km resolution grids as a dataset. Sort the dataset in chronological order from the front to the back. Starting from the very front of the dataset, select the daily data corresponding to the 1-km resolution grids with SNR being the corrected value and soil humidity being non-null during a certain period as the training set and the validation set. Take all the daily data corresponding to the remaining time as the test set. The test set includes the daily data corresponding to all 1-km resolution grids with null and non-null SNR values, where the SNR of the 1-km resolution grid with null SNR values is obtained through spatialization in Step 4. Step 6: Construct a mountainous GNSS-R soil moisture inversion model, including a random forest model, an XGBoost model, a support vector machine model, and a Stacking ensemble model. The altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and corrected SNR are all used as inputs for the random forest model, the XGBoost model, and the support vector machine model. Take the soil moisture prediction results output by the random forest model, the XGBoost model, and the support vector machine model as inputs for the Stacking ensemble model. The Stacking ensemble model assigns weights to the soil moisture prediction results output by the random forest model, the XGBoost model, and the support vector machine model, and outputs the final soil moisture prediction result. Use the training set and the validation set to train and validate the mountainous GNSS-R soil moisture inversion model to obtain a trained mountainous GNSS-R soil moisture inversion model. Step 7: Input the test set into the trained mountainous GNSS-R soil moisture inversion model to obtain the predicted daily soil moisture data at the 1-km resolution grid. Take the average of the predicted daily soil moisture data for each month to obtain the predicted monthly soil moisture data at the 1-km resolution grid.

[0007] Compared with the prior art, the present invention adopts the above technical solutions and has the following technical effects: The present invention first focuses on the algorithm design for mountainous areas with complex terrains, and proposes an inversion algorithm that integrates the terrain humidity index and GNSS-R data. The terrain humidity index can reflect the influence of terrain on the non-uniform distribution of soil moisture under the action of gravity. Integrating the terrain humidity index with CYGNSS data can make full use of the advantages of terrain factors and remote sensing technology, further improve the accuracy of soil moisture distribution and the prediction ability of spatio-temporal changes, so as to better reveal the dynamic changes of soil moisture in complex terrain areas such as mountainous areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] Figure 1 is the technical flow chart of the present invention; Figure 2 is the technical flow chart of the specific processing involved in the present invention; Figure 3 is the schematic diagram of the distribution of observation stations involved in the embodiments of the present invention; Figure 4 is the schematic diagram of the distribution of influencing factors; Figure 5 is the schematic diagram of the relationship between CYGNSS SNR and soil moisture; Figure 6 is the schematic diagram of the distribution of the terrain humidity index; Figure 7 is the flow chart of the method for spatializing CYGNSS SNR data using the terrain humidity index; Figure 8 is the schematic diagram of the results of the model inversion of soil moisture in the southern mountainous area of Anhui in January, April, July, and October 2022 involved in the embodiments of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0009] The following details the embodiments of the present invention, and the examples of the embodiments are shown in the drawings. The embodiments described below with reference to the drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0010] As Figure 1 and Figure 2 shown, the present invention proposes a method for integrated inversion of GNSS-R soil moisture in mountainous areas that integrates the terrain humidity index, including the following steps: Step S1, obtain the Cyclone Global Navigation Satellite System (CYGNSS) received data, digital elevation model data (DEM), Moderate Resolution Imaging Spectroradiometer (MODIS) observation data, and SoilGrids world soil information data of the target mountainous area; at the same time, obtain the soil moisture data of each observation station in the target mountainous area; CYGNSS satellite data comes from the Level 1 V3.0 data provided by the National Aeronautics and Space Administration (NASA) of the United States; DEM data comes from the Shuttle Radar Topography Mission (SRTM); SoilGrids world soil information data comes from the SoilGrids data developed by the International Soil Reference and Information Centre (ISRIC).

[0011] Taking the mountainous area in southern Anhui as an example, including seven cities of Huangshan, Wuhu, Ma'anshan, Tongling, Xuancheng, Chizhou, and Anqing, the data of 2021 and 2022 are selected to express the technical solution in the present invention in detail and completely.

[0012] Step S2, select four factors that have a profound impact on CYGNSS soil moisture inversion, namely vegetation, terrain, water body, and soil properties, as the parameters for model construction, as shown in Table 1 and Figure 4 shown. Among them, the vegetation factor includes the vegetation index (NDVI), the terrain factor includes altitude, slope, aspect, and longitude and latitude, the soil property includes clay content, and the water body factor includes monthly water body data; the monthly water body data is only used to exclude the water body area of the target area.

[0013] Table 1

[0014] Preprocess the acquired data: 2.1, conduct quality control on the measured soil moisture data. Screen the soil moisture data of each observation station, exclude the missing values and outliers of the data to ensure the accuracy of the data. The distribution of the observation stations is as Figure 3 shown.

[0015] 2.2, conduct quality control on the acquired CYGNSS data: exclude the observation data with an incident angle greater than 65° to initially ensure that the error caused by the terrain height is within a reasonable range; exclude the data with a signal-to-noise ratio less than 2 dB or greater than 14 dB; the receiver gain Gr less than 0 dB; the peak power of the reflected signal greater than -147 dB.

[0016] 2.3, calculate the vegetation index using the preprocessed MOD09Q1 dataset. The preprocessing steps include: format conversion, projection conversion, and data resampling; the specific formula for calculating the vegetation index is as follows: , In the formula, is the reflectance in the red light band of 620 nm to 670 nm, is the reflectance in the near-infrared band of 841 - 875 nm; NDVI is the vegetation index, and its positive value indicates the surface vegetation coverage, that is, the greater the coverage, the larger the NDVI value.

[0017] 2.4, Process the noise and outliers in the DEM data, and use the "Slope" calculation tool in ArcGIS software to calculate the slope based on the DEM data. The formula is as follows: , where, is the slope, is the elevation value on the terrain surface, and represent the change in spatial position in the horizontal direction.

[0018] In ArcGIS, use the Aspect tool in the Spatial Analyst tool to calculate the aspect. The aspect value is represented from 0 to 360 degrees, which is the main orientation of each raster cell, calculated clockwise, and the due north direction is 0 degrees. Then, use the Reclassify tool to reclassify the aspect value into 8 directions.

[0019] The DEM obtains the longitude and latitude data through the operation of arcgis software. Specifically: Extract points from the DEM, then use the points to extract the xy coordinates of each point, and then export each xy as a raster.

[0020] 2.5, The global gridded soil information (SoilGrids) is used to obtain the clay content. Since the penetration depth of the L-band signal is usually from the surface to 5 cm, only the surface layer (0 - 5 cm) data in SoilGrids is used.

[0021] Step S3, Calculate the signal-to-noise ratio (SNR) data using the CYGNSS data, and calculate the terrain moisture index using the elevation data; 3.1, Calibrate the signal-to-noise ratio (SNR) data using the CYGNSS data. SNR can be used to characterize the strength of the reflected signal. When the soil moisture increases, it will cause an increase in the soil dielectric constant, and then the signal reflection increases, resulting in a change in SNR. The schematic diagram of the relationship is as Figure 5 shown.

[0022] The formula for the calibrated signal-to-noise ratio is as follows: , In the formula, SNR represents the calibrated signal-to-noise ratio data; is the transmission power of the global navigation satellite system (GNSS) signal; is the gain of the transmitting antenna; is the gain of the receiving antenna; is the carrier wavelength of the global navigation satellite system signal; is the distance between the sampling point and the global navigation satellite system transmitter; is the distance between the sampling point and the Cyclone Global Navigation Satellite System (CYGNSS) receiver; DDM_SNR is the original signal-to-noise ratio data provided by CYGNSS.

[0023] 3.2, calculate the terrain humidity index. It can reflect the influence of terrain on the non-uniform distribution of soil humidity under the action of gravity, and the calculation formula is as follows: , In the formula, is the unit catchment area, is the slope of the grid cell. Usually, it is assumed that the effective contour length is consistent with the length of the unit grid. The slope is the maximum slope direction in the downhill direction, and the relevant definitions are as follows: , , In the formula, represents the total upstream area entering the grid cell; is the effective contour length, that is, the length of the unit grid; represents the height difference between adjacent grid cells; represents the length of the center of adjacent grid cells. Among them, can be calculated by using the flow direction tool in the arcgis software to calculate the flow direction of each grid, and then using the flow accumulation tool to calculate. The distribution of the terrain humidity index is as Figure 6 shown.

[0024] 3.3, perform spatio-temporal matching on all data. Resample all data (altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and terrain humidity index) to a 1km resolution by the nearest neighbor method; represent the signal-to-noise ratio data on the CYGNSS sampling points falling into the same grid within the same day (within 24 hours) in the form of an average to represent the signal-to-noise ratio data of the entire grid, and the signal-to-noise ratio of the grid without sampling points is a null value; the soil humidity data collected at the station represents the soil humidity data of the grid where the station is located, and the soil humidity of the grid without a station is a null value.

[0025] Step S4, establish a spatial fitting function between the terrain humidity index and the daily corrected SNR data, and spatialize the corrected SNR data; Only use the 1-km grid where the actual sampling points with valid values are located. Taking the terrain humidity index as the independent variable and the preprocessed daily CYGNSS SNR data with a 1-km resolution as the dependent variable, select the most suitable functional relationship to establish a spatially fitting function model that updates daily: , In this embodiment, linear fitting is preferably used as the spatially fitting function model: , For the fitting result, it is necessary to calculate the coefficient of determination value, set a threshold for the value. When is greater than the threshold, it indicates that the functional relationship is valid. When is less than the threshold, the data of that day is not retained to ensure the effectiveness of the relationship between the SNR data and the terrain humidity index data distribution. Using the established spatially fitting function model, the known terrain humidity index data can be used to spatialize the reasonable daily SNR data of the fitting result, and finally obtain the 1-km daily SNR data distribution. The detailed process of the spatialization part is as Figure 7 shown.

[0026] Step S5: Take the spatialized SNR data and other data (altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and preprocessed daily soil moisture data) as a data set, and divide the data set in the time series. The data corresponding to a part of the dates is used as the training set and the validation set, and only the 1-km resolution grid with the actual correction value of the signal-to-noise ratio and non-empty soil moisture under this period of dates is selected; the data corresponding to another part of the dates is used as the test set. At this time, all the data including the grid corresponding to the spatialized SNR is selected. In this embodiment, the data in 2021 is randomly divided into the training set and the validation set at a ratio of 80% and 20% for the construction of the optimal model; the data in 2022 is used as the test set for testing the model effect.

[0027] Step S6: Use the training set and the validation set to perform integrated learning on three machine learning algorithms, namely support vector machine (SVM), XGBoost, and random forest (RF), through the Stacking integrated model to obtain the optimal integrated inversion model.

[0028] Select the RF, XGBoost, and SVM models as the base learners of the Stacking integration method. Among them, RF can handle high-dimensional and non-linear data, provide high prediction accuracy, and has strong anti-noise ability and automatic feature selection ability; XGBoost is good at dealing with structured data and can optimize performance through the gradient boosting method, especially suitable for non-linear problems; SVM is particularly suitable for high-dimensional data and can find the optimal separation hyperplane in the high-dimensional space by using the kernel trick.

[0029] Take the outputs of the RF, XGBoost, and SVM models as the inputs of the Stacking integration model, and the final result of the site soil moisture prediction as the output of the Stacking integration model. That is, take the training set as the input item and input it into the 3 base models respectively, and then take the output items of the 3 base models as the inputs of the Stacking integration model to train the model.

[0030] The Stacking integration model sets weights for the outputs of the RF, XGBoost, and SVM models, and verifies the prediction results under different weight combinations through methods such as 10-fold cross-validation, and takes the best verified result as the optimal result. In this embodiment, the method of evenly distributing weights, that is, 1 / 3, 1 / 3, 1 / 3, is selected.

[0031] The present invention uses (root mean square error), (correlation coefficient) to characterize the accuracy, and the calculation formulas are as follows: , , Among them, is the number of samples, is the th observed value of the site soil moisture corresponding to the th sample, is the soil moisture value predicted by the model for the th sample, is the average value of all sample predicted values,

[0032] Set that when is greater than 0.7 and is less than 0.1, output the trained model. The verification result finally obtained in this embodiment is: , .

[0033] Step S7, input the test set into the trained model, and take the monthly average of the simulated daily values to obtain the final soil moisture result at the 1km monthly scale.

[0034] In this embodiment, the CYGNSS data and auxiliary data in 2022 are used as the model input. After taking the monthly average of the simulated daily values, the accuracy is tested on a monthly scale with the SMAP remote sensing soil moisture data and the station data as the true values respectively.

[0035] The 1-km soil moisture data inverted by the model is resampled to 9 km consistent with the SMAP data, and the calculated accuracy is , ; There is good consistency between the soil moisture data inverted by the model and the stations. The correlation coefficient reaches 0.93, being .

[0036] The soil moisture distribution in January, April, July, and October 2022 simulated by the model is as shown in Figure 8 .

[0037] Based on the same inventive concept, an embodiment of the present application provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of the above-mentioned method for inverting soil moisture of GNSS-R in mountainous areas by fusing terrain humidity index are implemented.

[0038] Based on the same inventive concept, an embodiment of the present application provides a computer-readable storage medium. The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the above-mentioned method for inverting soil moisture of GNSS-R in mountainous areas by fusing terrain humidity index are implemented.

[0039] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can 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.

[0040] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, as well as the combination of flows and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate for implementing in the processFigure 1 means for the functions specified in one process or multiple processes and / or one block or multiple blocks Figure 1 or multiple blocks.

[0041] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory produce a manufactured article including an instruction device that implements the functions specified in one process Figure 1 or multiple processes and / or one block Figure 1 or multiple blocks.

[0042] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one process Figure 1 or multiple processes and / or one block Figure 1 or multiple blocks.

[0043] The above embodiments are only for illustrating the technical idea of the present invention, and the protection scope of the present invention cannot be limited thereby. Any modification made on the basis of the technical solution according to the technical idea proposed by the present invention falls within the protection scope of the present invention.

Claims

1. A method for retrieving soil moisture in mountainous areas by integrating the terrain humidity index, characterized in that It includes the following steps: Step 1: Obtain remote sensing data of other areas except water body areas in the target mountainous area, including Cyclone Global Navigation Satellite System (CYGNSS) received data, Digital Elevation Model (DEM) data, Moderate Resolution Imaging Spectroradiometer (MODIS) observation data, and SoilGrids world soil information data; meanwhile, obtain soil moisture data of each observation station in the target mountainous area; preprocess the obtained remote sensing data and soil moisture data; Step 2: Select altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and signal-to-noise ratio as influencing factors for GNSS-R soil moisture inversion in the mountainous area. Obtain altitude, slope, aspect, and longitude and latitude data from the preprocessed DEM data, calculate the vegetation index using the preprocessed MODIS observation data, and obtain clay content data from the preprocessed SoilGrids world soil information data; Step 3: Correct the signal-to-noise ratio in the preprocessed CYGNSS received data, and calculate the terrain moisture index using the preprocessed DEM data; spatiotemporally match altitude, slope, aspect, longitude and latitude, vegetation index, clay content, the corrected signal-to-noise ratio, the terrain moisture index, and the preprocessed soil moisture data of the observation station into daily data under a 1-km resolution grid; Step 4: Obtain the daily data of the signal-to-noise ratio and the terrain moisture index corresponding to the 1-km resolution grid where the signal-to-noise ratio is not a null value. Take the signal-to-noise ratio as the dependent variable and the terrain moisture index as the independent variable, establish a daily corresponding spatial fitting function model, and use the daily corresponding spatial fitting function model to spatialize the signal-to-noise ratio of the 1-km resolution grid where the signal-to-noise ratio is a null value on the same day to obtain the signal-to-noise ratio of all 1-km resolution grids on that day; Step 5: Take the altitude, slope, aspect, longitude and latitude, vegetation index, clay content, the corrected signal-to-noise ratio, and the preprocessed soil moisture daily data corresponding to all 1-km resolution grids as a dataset. Sort the dataset in chronological order from the front to the back. Starting from the very front of the dataset, select the daily data corresponding to the 1-km resolution grids where the signal-to-noise ratio is the corrected value and the soil moisture is not a null value for a certain period of time as the training set and the validation set, and take all the daily data corresponding to the remaining time as the test set. The test set includes the daily data corresponding to all 1-km resolution grids where the signal-to-noise ratio is a null value and the signal-to-noise ratio is not a null value. The signal-to-noise ratio of the 1-km resolution grid where the signal-to-noise ratio is a null value is obtained by spatialization in Step 4; Step 6: Construct a mountainous GNSS-R soil moisture inversion model, including a random forest model, an XGBoost model, a support vector machine model, and a Stacking ensemble model. Altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and the corrected signal-to-noise ratio are used as inputs for the random forest model, the XGBoost model, and the support vector machine model. The soil moisture prediction results output by the random forest model, the XGBoost model, and the support vector machine model are used as inputs for the Stacking ensemble model. The Stacking ensemble model assigns weights to the soil moisture prediction results output by the random forest model, the XGBoost model, and the support vector machine model and outputs the final soil moisture prediction result. Use the training set and the validation set to train and validate the mountainous GNSS-R soil moisture inversion model to obtain a trained mountainous GNSS-R soil moisture inversion model. Step 7: Input the test set into the trained mountainous GNSS-R soil moisture inversion model to obtain daily soil moisture prediction data at a 1 km resolution grid. Take the average of the daily soil moisture prediction data for each month to obtain monthly soil moisture prediction data at a 1 km resolution grid.

2. The GNSS-R soil moisture inversion method for mountainous areas integrating terrain humidity index according to claim 1, wherein In step 1, preprocess the acquired remote sensing data and soil moisture data as follows: Preprocess the Cyclone Global Navigation Satellite System received data, including: sequentially removing received data with an incident angle greater than 65°, a signal-to-noise ratio less than 2 dB or greater than 14 dB, a receiver gain less than 0 dB, and a peak power of the reflected signal greater than -147 dB. Preprocess the digital elevation model data, including: removing noise data in the digital elevation model data. Preprocess the Moderate Resolution Imaging Spectroradiometer observation data, including: format conversion, projection conversion, and data resampling. Preprocess the soil moisture data, including: removing missing values and outliers in the soil moisture data observed at each observation station.

3. The GNSS-R soil moisture inversion method for mountainous areas integrating the terrain humidity index according to claim 1, wherein, In step 2, obtain altitude, longitude, and latitude from the preprocessed digital elevation model data. Use the slope calculation tool in ArcGIS software to calculate the slope. The formula for the slope is as follows: ， Among them, is the slope,[ is the elevation value on the terrain surface,[ and represent the change in spatial position in the horizontal direction.[ Use the aspect calculation tool in ArcGIS software to calculate the aspect, calculated clockwise with the due north direction being 0 degrees. The formula for the vegetation index is as follows: ， Among them, is the vegetation index, is the reflectance in the red light band, is the reflectance in the near-infrared band; Select the surface clay content data corresponding to the L band from the preprocessed SoilGrids world soil information data.

4. The GNSS-R soil moisture inversion method for mountainous areas integrating the terrain humidity index according to claim 1, wherein, In step 3, the formula for the corrected signal-to-noise ratio is as follows: ， Among them, is the corrected signal-to-noise ratio, is the transmission power of the global navigation satellite system signal, is the gain of the transmitting antenna, is the distance between the sampling point and the global navigation satellite system transmitter, is the distance between the sampling point and the Cyclone global navigation satellite system receiver, is the carrier wavelength of the global navigation satellite system signal, is the signal-to-noise ratio in the preprocessed Cyclone global navigation satellite system received data; The formula for the terrain moisture index is as follows: ， Among them, is the terrain humidity index, is the unit catchment area, is the slope of the grid cell, and ， ， wherein, is the total upstream area entering the grid cell, and the grid cell is the grid cell of the digital elevation model data; is the effective contour length, that is, the length of the unit grid; is the height difference between adjacent grid cells; is the length of the center of adjacent grid cells.

5. The GNSS-R soil moisture inversion method for mountainous areas integrating the terrain humidity index according to claim 1, characterized in that, In step 3, the spatio-temporal matching is as follows: Resample altitude, slope, aspect, longitude and latitude, vegetation index, clay content, and terrain moisture index to daily data at a 1 km resolution grid using the nearest neighbor method. Convert the signal-to-noise ratio data format to a 1-km resolution grid. That is, for each grid, calculate the average of the signal-to-noise ratios obtained from all sampling points that fall within the grid within 24 hours to obtain the daily signal-to-noise ratio at a 1-km resolution grid. Denote the daily signal-to-noise ratio of grids with no sampling points within 24 hours as null values. Convert the preprocessed soil moisture data to a 1-km resolution grid, that is, use the preprocessed soil moisture of each observation station as the soil moisture of the grid where the observation station is located. Denote the soil moisture of grids without observation stations as null values.

6. The GNSS-R soil moisture inversion method for mountainous areas integrating the terrain humidity index according to claim 1, characterized in that, In step 4, after establishing the corresponding daily spatial fitting function model, calculate the coefficient of determination for the spatial fitting function model. When the coefficient of determination is greater than a preset third threshold, it indicates that the established spatial fitting function model is valid, and retain the signal-to-noise ratio fitted by the spatial fitting function model established on the current day; otherwise, delete the signal-to-noise ratio fitted by the spatial fitting function model established on the current day.

7. The GNSS-R soil moisture inversion method for mountainous areas integrating terrain humidity index according to claim 1, characterized in that In step 6, when training and validating the mountainous GNSS-R soil moisture inversion model, use the root mean square error and the correlation coefficient to characterize the accuracy of the inversion model. When the root mean square error is less than a preset first threshold and the correlation coefficient is greater than a preset second threshold, the training of the mountainous GNSS-R soil moisture inversion model is completed. Among them, the calculation formulas for the root mean square error and the correlation coefficient are as follows: ， ， Among them, is the root mean square error, is the correlation coefficient, is the number of samples, is the soil moisture value predicted by the inversion model for the th sample, is the soil moisture observation value of the corresponding station of the th sample, is the average value of all sample observation values, is the average value of all sample predicted values.

8. A computer device, comprising a memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the mountainous GNSS-R soil moisture inversion method integrating the terrain humidity index according to any one of claims 1 to 7.

9. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the mountainous GNSS-R soil moisture inversion method integrating the terrain humidity index according to any one of claims 1 to 7.

Images