Chinese region soil humidity remote sensing inversion method based on CYGNSS satellite data

Abstract

The invention discloses a Chinese region soil humidity remote sensing inversion method based on CYGNSS satellite data. The method comprises the following steps: S1, preprocessing the collected actually measured soil humidity, delay-Doppler map, normalized vegetation index and land cover type; modeling samples matched with all the data are obtained through space-time matching; s2, taking the soil humidity as a target variable, selecting a CYGNSS trailing edge slope, an incident angle, a ground elevation, a surface roughness and a vegetation optical thickness as characteristic variables, and constructing a soil humidity remote sensing model; the modeling samples are divided into a training set and a test set, the soil humidity remote sensing model is trained, and the optimal soil humidity remote sensing model is stored; and S3, according to the optimal soil humidity remote sensing model, inverting daily soil humidity data from the CYGNSS trailing edge slope and related auxiliary data. According to the invention, the precision of the satellite-borne GNSS-R technology for remotely sensing and monitoring the soil humidity is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of soil moisture prediction, and in particular to a remote sensing inversion method for soil moisture in China based on CYGNSS satellite data. Background Art

[0002] As a key parameter of the water cycle, soil moisture characterizes the evaporation loss of surface water, photosynthesis of plants, and the activities of soil microorganisms. Traditional soil moisture monitoring methods are time-consuming, labor-intensive, inefficient, and have few sites, and can no longer meet current needs. In recent years, with the continuous development of global navigation satellite system reflection measurement technology (GNSS Reflectomtry, GNSS-R), the L-band signal emitted by GNSS satellites is sensitive to soil moisture. After the GNSS satellite signal is reflected by the earth's surface, it will interfere with the GNSS satellite direct signal at the GNSS receiver, and the interference signal will change with the different surface conditions. Therefore, GNSS reflection measurement technology has become a new method for soil moisture monitoring. At present, there are two main ways to measure soil moisture with GNSS-R technology: ground-based and satellite-based. Ground-based GNSS-R measurement can obtain soil moisture at the site with high precision, high observation frequency and long time series, such as the Chinese invention patent "Soil moisture estimation method based on ground-based GNSS-R data based on corrected phase" with publication number CN 115980317 B. However, due to the limited distribution of sites, this method cannot obtain soil moisture in a wide range; satellite-borne GNSS-R can measure surface soil moisture quickly, periodically and over a wide range, which has great advantages in soil moisture monitoring. Since the reflection of the GNSS signal from the surface has two parts, coherent components and incoherent components, the surface roughness will affect the surface reflection characteristics of the GNSS signal. In areas with less surface roughness, the reflection is mainly coherent; as the surface roughness increases, the effective reflection of the coherent component decreases rapidly, while the effective reflection of the incoherent component increases. However, the surface reflectivity used in the current satellite-borne GNSS-R remote sensing monitoring of soil moisture is often calculated under the assumption that the surface is smooth and reflects the coherent component, ignoring the contribution of the incoherent component to the surface reflection. In addition, when the GNSS signal from the navigation satellite is incident on the surface and then reflected by the surface to reach the satellite receiver, it is attenuated by the vegetation covering the surface. Therefore, accurately determining the surface reflection characteristics, surface roughness and vegetation optical thickness are key parameters for high-precision remote sensing monitoring of soil moisture by satellite-borne GNSS-R.

[0003] At present, when spaceborne GNSS-R remote sensing inverts soil moisture, surface roughness and vegetation optical thickness mainly come from SMAP satellite soil moisture products. Even the soil moisture used as the true value in the modeling process also comes from SMAP satellite soil moisture products, such as the Chinese invention patent with publication number CN 114371182 B "A GNSS-R high-precision soil moisture estimation method based on CYGNSS data", the Chinese invention patent with publication number CN 114371182 B "Method for collaborative inversion of surface soil moisture and vegetation water content by integrating spaceborne GNSS-R and multi-source RS data", and the paper "Maria Paola Clarizia et al. Analysis of CYGNSS Data for Moisture Retrieval, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing". Sensing, 2019, 12(7): 2227-2235", Due to the large range of data missing in the SMAP satellite soil moisture product and the limited accuracy of the satellite product data, errors will accumulate when modeling soil moisture remote sensing based on GNSS-R satellite data, making it difficult to achieve high-precision remote sensing inversion of satellite-borne GNSS-R soil moisture. The CYGNSS soil moisture product currently released by NASA is obtained by inverting the SMAP satellite soil moisture product as the true value and the CYGNSS surface reflectivity as the remote sensing variable using a univariate linear regression model. It has not been corrected for surface roughness and vegetation optical thickness, and the accuracy of its soil moisture product needs to be further improved. Summary of the invention

[0004] Purpose of the invention: The purpose of the present invention is to provide a remote sensing inversion method for soil moisture in China based on CYGNSS satellite data, so as to improve the accuracy of modeling in the process of satellite-borne GNSS-R remote sensing monitoring of soil moisture.

[0005] Technical solution: A remote sensing inversion method for soil moisture in China based on CYGNSS satellite data, including the following steps: S1, pre-processing the collected measured soil moisture, delay-Doppler map, normalized vegetation index and land cover type; obtaining modeling samples matching all data through time-space matching; obtaining the CYGNSS trailing edge slope after processing the delay-Doppler map, the normalized vegetation index is used for calculating the vegetation optical thickness, and the land cover type is used for calculating the surface roughness; S2, taking soil moisture as the target variable, by calculating the Pearson correlation coefficient between the target variable and the characteristic variable, selecting CYGNSS trailing slope, incident angle, ground elevation, surface roughness and vegetation optical thickness as characteristic variables, and constructing a soil moisture remote sensing model; dividing the modeling samples into a training set and a test set, and training the soil moisture remote sensing model, and saving the best soil moisture remote sensing model; S3, based on the best soil moisture remote sensing model, invert daily soil moisture data from the CYGNSS trailing edge slope and its related auxiliary data.

[0006] Furthermore, the measured soil moisture data of each meteorological station were screened. First, missing values ​​and outliers were eliminated, and then the mean and standard deviation of the overall data were calculated. Data that exceeded 2 times the standard deviation of the mean were eliminated.

[0007] Furthermore, the steps to obtain the CYGNSS trailing edge slope after processing the delay-Doppler map are as follows: SD1, subtract the noise floor value from the power value corresponding to each delay-Doppler graph, and the noise calculation formula is: , where , are the lower and upper boundaries of the delay chips selected for calculating the noise region, represents the delay chip; , They are the lower and upper boundaries of the frequency selected for calculating the noise region; M is the number of pixels in the selected noise region; Then the calculation formula of the delay-Doppler map after denoising can be expressed as: , where is the delay-Doppler diagram after denoising; DDM is the delay-Doppler diagram; Noise represents noise; SD2, yes Integrate the reflected power in the Doppler dimension to obtain the integrated delay waveform. The calculation formula is: , where is the integrated time delay waveform; N is the number of selected Doppler frequencies; is the kth Doppler frequency; SD3, normalizes the one-dimensional integrated delay waveform, and the calculation formula is: , where is the normalized integrated delay waveform, is the maximum value in the integrated delay waveform, i takes values ​​from 1 to 17; represents the integrated delay waveform of the i-th delay chip; SD4, a linear function is used to perform least squares fitting on the trailing edge of the normalized integrated delay waveform, and the absolute value of the slope is taken as the trailing edge slope.

[0008] Furthermore, the implementation steps for obtaining modeling samples that match all data are as follows: S11, by calculating the distance between each ground reflection point of the CYGNSS satellite and all soil moisture automatic monitoring stations, the point pairs with a distance less than 7 km between the two are taken as spatial matching sites; S12, time matching the CYGNSS data with the soil moisture data of the corresponding station with the same observation date and time; S13, matching the matched CYGNSS-measured soil moisture sample points with auxiliary data, wherein the auxiliary data includes vegetation optical thickness and surface roughness; S14, according to the observation time of the CYGNSS satellite data, select the auxiliary data of the corresponding time period for matching; then correspond the longitude and latitude coordinates of each reflection point of the CYGNSS satellite to the pixel where the planar auxiliary data is located, so that the data of each reflection point of the CYGNSS satellite is matched with the corresponding pixel value of the auxiliary data, and generate a modeling sample with all data matching.

[0009] Furthermore, in step S2, a soil moisture multivariate linear regression remote sensing model, a soil moisture support vector remote sensing model and a soil moisture random forest remote sensing model are constructed respectively, and the three models are trained and tested respectively, and the best remote sensing models are saved respectively; soil moisture predictions are performed on the three best remote sensing models respectively, and the remote sensing model with the highest accuracy is selected as the soil moisture remote sensing model.

[0010] Compared with the prior art, the present invention has the following significant effects: 1. The present invention comprehensively considers the effects of coherent and incoherent reflections of GNSS signals on the ground, and proposes the trailing edge slope (TES) as the main variable for remote sensing monitoring of soil moisture based on the shape characteristics of the time delay-Doppler map (DDM) of the satellite-borne GNSS-R observation quantity. This solves the problem that the surface reflectivity in the previous soil moisture remote sensing only considers the coherent component and ignores the contribution of the incoherent component to the surface reflection, thereby improving the accuracy of remote sensing monitoring of soil moisture using satellite-borne GNSS-R technology; 2. Since the surface roughness will affect the reflection characteristics of the surface GNSS signal and the vegetation coverage will attenuate the GNSS signal, the surface roughness and the vegetation optical thickness are two important factors affecting the satellite-borne GNSS-R remote sensing monitoring of soil moisture. The present invention proposes a method for accurately calculating these two parameters using the MODIS surface product of the Earth Observation System (EOS) satellite, which provides accurate auxiliary variables for the satellite-borne GNSS-R high-precision remote sensing inversion of soil moisture. 3. In response to the problem of inaccurate modeling parameters in the current satellite-borne GNSS-R remote sensing monitoring of soil moisture, the present invention uses the measured soil moisture data in China as the target variable, the CYGNSS satellite trailing edge slope data as the remote sensing variable, and combines the vegetation optical thickness and surface roughness data. It adopts multiple regression and machine learning methods to construct the optimal CYGNSS soil moisture remote sensing inversion model suitable for China, which solves the problem of inaccurate modeling parameters in previous satellite-borne GNSS-R soil moisture remote sensing methods, and obtains daily and monthly soil moisture remote sensing products in China inverted from CYGNSS data, providing an important data source for agricultural climate drought and hydrometeorological research in China. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 is a flow chart of the present invention; Figure 2 (a) is the CYGNSS satellite delay-Doppler diagram, and (b) is the normalized integrated delay waveform diagram; Figure 3 A scatter plot of measured values ​​and predicted values ​​during training and testing of the multivariate linear remote sensing model of soil moisture in an embodiment of the present invention; Figure 4 A scatter plot of measured values ​​and predicted values ​​during training and testing of the soil moisture support vector remote sensing model in an embodiment of the present invention; Figure 5 It is a scatter plot of measured values ​​and predicted values ​​during training and testing of the random forest remote sensing model of soil moisture in an embodiment of the present invention. DETAILED DESCRIPTION

[0012] The present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0013] The present invention provides a method for remote sensing inversion of soil moisture in China based on CYGNSS satellite data. In order to verify the method proposed in the present invention, the data of four typical months in January, April, July and October 2018, representing different seasons, are used to construct the optimal CYGNSS soil moisture remote sensing inversion model suitable for China. The soil moisture data of four typical months in February, May, August and November 2019 are remotely sensed and inverted to generate daily and monthly soil moisture remote sensing products in China. The implementation process of the technical solution is shown in the attached figure. Figure 1 As shown, the following steps are included: Step 1, data input; The data that need to be input in the present invention mainly include: The soil moisture data measured by regional meteorological stations in China contain information such as station latitude and longitude, elevation, and observation date, which are used to construct remote sensing inversion models; CYGNSS satellite data include delay-Doppler map (DDM), specular reflection point latitude (sp_lat), specular reflection point longitude (sp_lon), specular reflection point incidence angle (sp_inc_angle) and other data sets. DDM is used to calculate ground reflectivity as the main variable for soil moisture remote sensing; TERRA or AQUA satellite MODIS vegetation index products, including the normalized difference vegetation index for the calculation of vegetation optical depth; The land cover type data classified by the International Geosphere biosphere Programme (IGBP) of MODIS are used to calculate the surface roughness.

[0014] Delay-Doppler maps (DDM) are the signal power around the mirror reflection point in the ground blaze area obtained by the satellite-borne GNSS-R receiver. It is a two-dimensional distribution diagram of the related power of delay and Doppler frequency shift. Its shape and amplitude can be used to infer the surface reflection characteristics. When the surface is relatively smooth, the signal reflected by the coherent component will show a typical "horseshoe" shape on the DDM; with the influence of surface roughness and soil moisture, the signal reflected by the incoherent component will show various delay wave shapes. Therefore, the leading edge slope (LES) and trailing edge slope (TES) calculated from the DDM normalized delay waveform can well reflect the degree of coherent / incoherent reflection on the surface.

[0015] The data used for soil moisture remote sensing modeling in the embodiment of the present invention are: measured soil moisture data in four typical months of January, April, July and October 2018 in my country, CYGNSS delay-Doppler map (DDM) data and its auxiliary parameters including the incident angle during the same period, MODIS (Moderate-resolution Imaging Spectroradiometer) normalized vegetation index product during the same period and MODIS IGBP global land cover type data.

[0016] Step 2, data preprocessing; (A1) Quality control of measured soil moisture data In order to ensure the validity of the data when modeling soil moisture remote sensing, the measured soil moisture data of each meteorological station were screened. First, missing values ​​and outliers were eliminated, and then the mean and standard deviation of the overall data were calculated. The data that exceeded the range of 2 times the standard deviation of the mean was eliminated to ensure the accuracy of the measured soil moisture data and improve the accuracy of the soil moisture remote sensing model.

[0017] (A2) Processing of CYGNSS satellite data Quality control is performed on CYGNSS satellite data. For data in the CYGNSS data set with a DDM signal-to-noise ratio less than 0db, a receiver antenna gain less than 0dB, and an incident angle greater than 65°, it means that the received signal quality is relatively poor and should be eliminated.

[0018] The CYGNSS satellite GNSS-R receiver has four reflection measurement channels. Each channel observes DDM data of nearly 100,000 reflection points worldwide. The DDM data of each reflection point is a two-dimensional distribution diagram of the relevant power composed of 17 time delay units (Bin) and 11 Doppler frequency shift units (Bin), such as Figure 2 Calculating the CYGNSS trailing edge slope (TES) from DDM data is done in 4 steps: (D1) DDM data denoising Subtract the noise floor from the power value corresponding to each DDM, and the noise calculation formula is: (1) In the formula, , The lower and upper boundaries of the delay chips selected for calculating the noise region respectively; , are the lower and upper boundaries of the frequency selected for calculating the noise area; M is the number of pixels in the selected noise area. The noise area selected by the present invention is the pixels of all frequencies of the first 4 time delays, and M is 44 pixels. The calculation formula of the DDM after denoising can be expressed as: (2) In the formula, is the DDM after denoising, and Noise represents noise.

[0019] (D2) Calculate the integral delay waveform The reflected power of the denoised DDM in the Doppler dimension is integrated to obtain the integrated delay waveform (IDW), which is calculated as follows: (3) In the formula, is the normalized integrated delay waveform; N is the number of selected Doppler frequencies, and the present invention takes all frequencies as 11; is the kth Doppler frequency.

[0020] (D3) Normalization of the integrated delay waveform The one-dimensional integrated delay waveform is normalized, and its calculation formula can be expressed as: (4) In the formula, is the normalized integrated delay waveform, is the maximum value in the integrated delay waveform, i takes a value between 1 and 17; Represents the integrated delay waveform of the i-th delay code chip. Figure 2 (b) is a normalized integrated delay waveform diagram in an embodiment of the present invention.

[0021] (D4) Trailing edge slope calculation The trailing edge of the normalized integrated delay waveform is composed of N sampling points after the waveform peak, and shows a decreasing rule with the delay. The trailing edge slope can be obtained by performing the least squares fitting of the trailing edge with a linear function. ) The last three waveform values Perform linear regression and obtain the absolute value of the slope, which is the CYGNSS trailing edge slope TES.

[0022] (A3) Calculation of vegetation optical depth When the electromagnetic wave signal emitted by the GNSS satellite encounters vegetation, it will propagate from the top of the vegetation canopy to the ground, and then from the ground to the top of the vegetation canopy. During the GNSS signal transmission process, it will experience secondary attenuation of the vegetation. Depends on the vegetation optical thickness and the angle of incidence, the calculation formula is as follows: (5) (6) In the formula, represents the angle of incidence, represents vegetation optical depth, VWC represents vegetation water content, and b represents the scaling factor related to land cover type.

[0023] The proportional factor b is related to the vegetation type, canopy structure, growth condition, and the frequency of the GNSS satellite transmission signal. In this embodiment, the proportional factor b is determined by the lookup table method using the land cover type data classified by MODIS IGBP in 2018 (see Table 1).

[0024] Table 1 Vegetation parameters of different land cover types

[0025] The vegetation water content (VWC) can be measured by the normalized vegetation index (NDVI) in the corresponding period of 2018 and the maximum value of the normalized vegetation index in the past 10 years. and minimum value , and the vegetation stem factor (stem_factor) are calculated according to the following empirical formula (3): (7) In the formula, the vegetation stem factor is determined by the lookup table method based on the 2018 MODIS IGBP land cover type; the maximum value of the normalized difference vegetation index in the past 10 years is and minimum value It can be calculated through programming on the Google Earth Engine (GEE) platform.

[0026] (A4) Determination of surface roughness When the GNSS satellite transmits a signal incident on the smooth surface of the ground, it will show the characteristics of mirror reflection, and the signal received by the GNSS receiver is mainly coherent scattering. As the surface roughness increases, the amount of incoherent scattering on the surface gradually increases, while incoherent scattering is insensitive to polarization characteristics, and the difference between horizontal and vertical polarization is small. For surface roughness, the root mean square height h is usually used in the vertical direction of the surface, which is related to the surface vegetation cover type, terrain undulation, and the polarization mode and frequency of the electromagnetic wave. In this embodiment, the surface roughness is determined by the lookup table method using the 2018 MODIS IGBP land cover type data (Table 1).

[0027] (A5) Spatial and temporal matching of data Since the distance between ground reflection points in the same scanning band of CYGNSS satellite is about 7 km, the position where CYGNSS satellite receives the signal from the ground reflection point changes every day, while the position of the soil moisture automatic monitoring station is fixed. When performing spatial matching of the two types of data, the distance between each ground reflection point of CYGNSS satellite and all soil moisture automatic monitoring stations is calculated every day, and the point pairs with a distance less than 7 km between the two are taken as spatial matching sites.

[0028] As a polar-orbiting satellite, CYGNSS satellite passes through from 00:00 to 23:00 every day, while soil moisture is observed hourly from 00:00 to 23:00 every day. After completing the spatial matching of the two data, the CYGNSS data is temporally matched with the soil moisture data of the corresponding station at the same date and time, thus completing the temporal matching of the two data.

[0029] Since CYGNSS satellite data is point data, and the two auxiliary data (vegetation optical thickness and surface roughness) are surface data, firstly, according to the observation time of CYGNSS satellite data, the auxiliary data of the corresponding period is selected; secondly, according to the longitude and latitude coordinates of each reflection point of CYGNSS satellite, it is matched to the pixel where the surface auxiliary data is located, so that the data of each reflection point of CYGNSS satellite is matched with the corresponding pixel value of the auxiliary data, and the samples matching the two data are generated.

[0030] To obtain modeling samples that match all data, the steps are as follows: Step 11: Calculate the distance between each ground reflection point of the CYGNSS satellite and all soil moisture automatic monitoring stations, and select the point pairs with a distance less than 7 km as the spatial matching sites.

[0031] Step 12, time-matching the CYGNSS data with the soil moisture data of the corresponding station with the same observation date and time, that is, completing the spatiotemporal matching of the two data.

[0032] Step 13, match the matched CYGNSS-measured soil moisture sample points with two auxiliary data (vegetation optical thickness and surface roughness).

[0033] Step 14, according to the observation time of the CYGNSS satellite data, select the auxiliary data of the corresponding time period; then, according to the longitude and latitude coordinates of each reflection point of the CYGNSS satellite, correspond to the pixel where the planar auxiliary data is located, so that the data of each reflection point of the CYGNSS satellite is matched with the corresponding pixel value of the auxiliary data, and generate a modeling sample with matching data.

[0034] Step 3, construction of soil moisture remote sensing model; (B1) Model feature parameter selection The present invention uses soil moisture data as the target variable, and uses CYGNSS trailing edge slope TES, incident angle , Ground elevation Surface roughness and vegetation optical depth As the characteristic variable, the characteristic parameters of the soil moisture remote sensing model are screened. Since the selection of characteristic variables has a great influence on the accuracy of the model, when the correlation between the characteristic variables and the dependent variable is low, it is difficult for the constructed model to have a good fitting effect; when the correlation is too high, the model will fall into a flat state, resulting in low model efficiency. In order to obtain a high-precision and high-efficiency model, it is necessary to screen the characteristic variables of the model, calculate the Pearson correlation coefficient between the characteristic variables and the characteristic variables, between the characteristic variables and the dependent variable, and calculate the variance inflation factor (VIF) between the dependent variable and each characteristic variable. Finally, according to the Pearson correlation coefficient less than 0.8 and the variance inflation factor less than 5, it is used as the optimal characteristic parameter for constructing the soil moisture remote sensing model. Table 2 lists the Pearson correlation coefficient and variance inflation factor between soil moisture and each characteristic variable.

[0035] Table 2 Statistics of soil moisture and characteristic variables

[0036] As shown in Table 2, all correlation coefficients are less than 0.8 and the variance inflation factor is less than 5, indicating that there is no multicollinearity problem between soil moisture and each characteristic variable, and the five characteristic variables can be used for soil moisture remote sensing modeling.

[0037] (B2) Constructing a soil moisture remote sensing model After the model characteristic parameters were selected, 50% of all modeling samples were randomly selected as the training data set for the construction of the soil moisture remote sensing model, and the remaining 50% were used as the test set of the model. Three methods, namely multivariate linear regression, support vector regression and random forest, were selected to construct the soil moisture remote sensing model.

[0038] In the embodiment of the present invention, the number of samples matching the soil moisture of each measured site with the five characteristic variable data is 128,518 pairs in total. 50% of the modeling samples are randomly selected as the training set for establishing the model, with a total of 64,259 pairs of samples, and the remaining 50% of the modeling samples are used as the test set for model prediction, with a total of 64,259 pairs of samples.

[0039] (i) Multiple linear regression model Assume that soil moisture (SM) is related to five characteristic variables: CYGNSS trailing edge slope TES and incident angle , Ground elevation Surface roughness and vegetation optical depth There is a multivariate linear correlation. The multivariate linear regression analysis method is used to establish a soil moisture remote sensing model using training samples. The general form of the model is as follows: (8) Where: is the constant term of the model, is the regression coefficient of the model. The soil moisture remote sensing model was constructed by performing multivariate linear regression on the training set of soil moisture samples using the least squares method.

[0040] For the training samples, the multivariate linear regression analysis method is used to construct the soil moisture remote sensing model:

[0041] Where: SM is the soil volume water content, in cm 3 / cm 3 ; Angle of incidence V z The unit is degree ( ); Ground elevation The unit is kilometer (Km); surface roughness and vegetation optical thickness are dimensionless.

[0042] After the t-test, the p-value of the regression coefficient of each characteristic variable is less than 0.01, indicating that the linear influence of each characteristic variable on the corresponding variable reaches an extremely significant level. The scatter plot of the measured value and the predicted value during the training and testing of the soil moisture multivariate linear regression remote sensing model is shown in Figure 3 .

[0043] (ii) Support Vector Regression Model There is a function and two parameters in the support vector regression (SVR) model that have a great influence on the results of the model, namely the kernel function, the penalty coefficient C and the kernel coefficient gamma. The kernel function determines the accuracy and universality of the SVR model. Commonly used kernel functions include: linear kernel function, polynomial kernel function, Gaussian kernel function and Sigmoid kernel function. Through experiments on four kernel functions, the present invention selects the Gaussian kernel function with better smoothing performance for the construction of soil moisture support vector remote sensing model.

[0044] The default value of the penalty coefficient C is 1.0. The C value controls the complexity of the model and the tolerance to errors. The larger the C value, the higher the degree of fit of the model to the training data, that is, the model is less likely to tolerate errors in the training set, which may cause the model to overfit; on the contrary, a smaller C will cause the model to tolerate more errors, which helps to improve the generalization ability of the model. The present invention adopts a trial number method, and takes the C value as [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0] respectively, tests all training sample sets, and selects the optimal C value by calculating the correlation coefficient R and the root mean square error RMSE between the measured value of soil moisture and the predicted value of the SVR model.

[0045] The kernel coefficient gamma defines the influence range of the kernel function, which affects the fitting of the model on the training data. A larger gamma value will cause the model to only focus on sample points close to the support vector, which may cause the model to overfit the training data; a smaller gamma value will increase the range of the support vector, which may cause the model to underfit. The present invention adopts a trial number method, takes the gamma value as [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], tests all training sample sets, and selects the optimal gamma value by calculating the R and RMSE of the measured value of soil moisture and the predicted value of the SVR model.

[0046] Therefore, a soil moisture support vector remote sensing model with optimal parameters is constructed.

[0047] For the training samples, the support vector regression method is used. The kernel function uses the Gaussian kernel function with good smoothing performance. The penalty coefficient C values ​​are [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0], and the kernel coefficient gamma values ​​are [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]. All training sample sets are tested to determine the optimal parameters in the support vector regression model. The soil moisture support vector regression remote sensing model is established. The scatter plots of the measured values ​​and the predicted values ​​during model training and testing are shown in Figure 2. Figure 4 .

[0048] (iii) Random Forest Model This embodiment focuses on optimizing and selecting two key parameters in the random forest model (the number of decision trees and the feature value of each decision tree).

[0049] For the training samples, the random forest method is used. The number of decision trees is set to {50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}, and the feature value of each decision tree is set to {1, 2, 3, 4, 5}. The 5-fold cross-validation method is used for the training sample data set. The number of decision trees and the feature value of each decision tree are continuously optimized from different combinations to improve the model accuracy. The optimal number of decision trees and the feature value of each decision tree are determined to construct a random forest remote sensing model of soil moisture. The scatter plots of measured values ​​and predicted values ​​during model training and testing are shown in Figure 2. Figure 5 .

[0050] Table 3 lists the statistical parameters of the three soil moisture remote sensing models during training and testing. From the perspective of correlation coefficient R, root mean square error RMSE and mean absolute error MAE, the training accuracy of the soil moisture multivariate linear remote sensing model is equivalent to the prediction accuracy, while the training accuracy of the soil moisture support vector remote sensing model and the random forest remote sensing model is better than the prediction accuracy, especially the training accuracy of the random forest model is very high. Among the three soil moisture remote sensing models, the accuracy of the multivariate linear model is relatively low, the accuracy of the support vector model is second, and the accuracy of the random forest model is the highest. Therefore, the random forest remote sensing model is the best model for soil moisture remote sensing inversion in the embodiment of the present invention.

[0051] Table 3 Statistical parameters of different soil moisture remote sensing models during training and testing

[0052] Step 4, remote sensing inversion of soil moisture; (C1) Inversion data preparation First, the daily trailing edge slope data of the eight small satellites of the CYGNSS constellation are calculated and the corresponding incident angles are extracted. For the convenience of data processing, the data of the eight CYGNSS satellites can be merged every day.

[0053] The embodiment of the present invention uses the constructed optimal soil moisture random forest remote sensing model to invert soil moisture data for four typical months: February, May, August and November 2019. First, the daily trailing edge slope data of the eight small satellites of the CYGNSS constellation in four months of 2019 are calculated, the corresponding angle of incidence is extracted, and the data of the eight CYGNSS satellites are merged every day. Then, according to the longitude and latitude coordinates of the CYGNSS surface reflection point, the elevation value corresponding to the CYGNSS ground reflection point is extracted from the digital elevation model (DEM) of China; the vegetation optical thickness value corresponding to the CYGNSS ground reflection point is extracted from the vegetation optical thickness data calculated by the MODIS global vegetation index product; the surface roughness value corresponding to the CYGNSS ground reflection point is extracted from the surface roughness data calculated by the MODIS IGBP land cover type product. Finally, the daily CYGNSS trailing edge slope and incident angle data, as well as the corresponding ground elevation, surface roughness, and vegetation optical thickness, are merged as independent variables for remote sensing inversion of soil moisture, in preparation for remote sensing inversion of soil moisture.

[0054] (C2) Soil moisture remote sensing product generation According to the best soil moisture remote sensing model, the daily soil moisture data is inverted from the CYGNSS trailing slope and its related auxiliary data. The latitude and longitude coordinates of the CYGNSS surface reflection points are used to generate daily soil moisture point vector files, and then the inverse distance weighted method (IDW) is used to perform spatial interpolation on the soil moisture data of each reflection point to generate soil moisture remote sensing daily product data. The soil moisture daily product data of each month is averaged to generate soil moisture remote sensing monthly product data.

Claims

1. A remote sensing inversion method for soil moisture in China based on CYGNSS satellite data, characterized in that: The steps include: S1, pre-processing the collected measured soil moisture, delay-Doppler map, normalized vegetation index and land cover type; obtaining modeling samples matching all data through time-space matching; obtaining the CYGNSS trailing edge slope after processing the delay-Doppler map, the normalized vegetation index is used for calculating the vegetation optical thickness, and the land cover type is used for calculating the surface roughness; S2, taking soil moisture as the target variable, by calculating the Pearson correlation coefficient between the target variable and the characteristic variable, selecting CYGNSS trailing slope, incident angle, ground elevation, surface roughness and vegetation optical thickness as characteristic variables, and constructing a soil moisture remote sensing model; dividing the modeling samples into a training set and a test set, and training the soil moisture remote sensing model, and saving the best soil moisture remote sensing model; S3, based on the best soil moisture remote sensing model, invert daily soil moisture data from the CYGNSS trailing edge slope and its related auxiliary data.

2. The method for remote sensing inversion of soil moisture in China based on CYGNSS satellite data according to claim 1 is characterized in that: The measured soil moisture data of each meteorological station were screened. First, missing values ​​and outliers were eliminated, and then the mean and standard deviation of the overall data were calculated. Data that exceeded 2 times the standard deviation of the mean were eliminated.

3. The method for remote sensing inversion of soil moisture in China based on CYGNSS satellite data according to claim 1 is characterized in that: The steps to obtain the CYGNSS trailing edge slope after processing the delay-Doppler map are as follows: SD1, subtract the noise floor value from the power value corresponding to each delay-Doppler graph, and the noise calculation formula is: , where , are the lower and upper boundaries of the delay chips selected for calculating the noise region, represents the delay chip; , They are the lower and upper boundaries of the frequency selected for calculating the noise region; M is the number of pixels in the selected noise region; Then the calculation formula of the delay-Doppler map after denoising can be expressed as: , where is the delay-Doppler diagram after denoising; DDM is the delay-Doppler diagram; Noise represents noise; SD2, yes Integrate the reflected power in the Doppler dimension to obtain the integrated delay waveform. The calculation formula is: , where is the integrated time delay waveform; N is the number of selected Doppler frequencies; is the kth Doppler frequency; SD3, normalizes the one-dimensional integrated delay waveform, and the calculation formula is: , where is the normalized integrated delay waveform, is the maximum value in the integrated delay waveform, i takes values ​​from 1 to 17; represents the integrated delay waveform of the i-th delay chip; SD4, a linear function is used to perform least squares fitting on the trailing edge of the normalized integrated delay waveform, and the absolute value of the slope is taken as the trailing edge slope.

4. The method for remote sensing inversion of soil moisture in China based on CYGNSS satellite data according to claim 1 is characterized in that: The implementation steps for obtaining modeling samples that match all data are as follows: S11, by calculating the distance between each ground reflection point of the CYGNSS satellite and all soil moisture automatic monitoring stations, the point pairs with a distance less than 7 km between the two are taken as spatial matching sites; S12, time matching the CYGNSS data with the soil moisture data of the corresponding station with the same observation date and time; S13, matching the matched CYGNSS-measured soil moisture sample points with auxiliary data, wherein the auxiliary data includes vegetation optical thickness and surface roughness; S14, according to the observation time of the CYGNSS satellite data, selecting the auxiliary data of the corresponding period for matching; Then, according to the longitude and latitude coordinates of each reflection point of the CYGNSS satellite, it corresponds to the pixel where the planar auxiliary data is located, so that the data of each reflection point of the CYGNSS satellite is matched with the corresponding pixel value of the auxiliary data to generate a modeling sample with matching data.

5. The method for remote sensing inversion of soil moisture in China based on CYGNSS satellite data according to claim 1 is characterized in that: In step S2, a soil moisture multivariate linear regression remote sensing model, a soil moisture support vector remote sensing model and a soil moisture random forest remote sensing model are constructed respectively, and the three models are trained and tested respectively, and the best remote sensing models are saved respectively; soil moisture prediction is performed for the three best remote sensing models respectively, and the remote sensing model with the highest accuracy is selected as the soil moisture remote sensing model.

Images