Satellite-borne GNSS-R important feature selection and forest above-ground biomass and canopy height estimation method combining tree model ensemble learning algorithm and SHAP interpreter

Abstract

The invention discloses a satellite-borne GNSS-R important feature selection and forest above-ground biomass and canopy height estimation method in combination with a tree model ensemble learning algorithm and an SHAP interpreter. According to the method, three efficient gradient boosting decision tree algorithms, namely CatBoost, XGBoost and LightGBM, and an SHAP interpreter are combined to respectively select important features for modeling in inversion of forest above-ground biomass and canopy height; and on the basis of three tree model ensemble learning algorithms of CatBoost, XGBoost and LightGBM, a bistatic radar scattering cross section, an effective scattering area and selected important features are fused to construct a global high space-time forest above-ground biomass and canopy height estimation model. According to the method, a tree model ensemble learning algorithm and satellite-borne GNSS-R data and multivariable parameters of an SHAP interpreter are combined for collaborative inversion of the global high space-time forest above-ground biomass and canopy height, and a new forest above-ground biomass and canopy height obtaining means is provided; and a satellite-borne GNSS-R inversion forest above-ground biomass and canopy height method taking multiple influence factors into consideration is realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of measuring forest aboveground biomass and canopy height by using GNSS reflection signals, and specifically relates to a spaceborne GNSS-R important feature selection and forest aboveground biomass and canopy height estimation method combining a tree model ensemble learning algorithm and a SHAP interpreter. Background Art

[0002] Forest above-ground biomass (AGB) is an important climate variable. It plays a vital role in global carbon distribution and carbon storage, and can help mitigate climate change. At the same time, AGB also plays a key role in maintaining biodiversity in forests. Canopy height (CH) refers to the vertical distance from the highest point of the canopy in the forest to the ground. It can reflect the structure and growth status of the forest and is closely related to AGB. Traditional satellite missions for measuring AGB include synthetic aperture radar, but it is limited by saturation; microwave radiation measurement also plays an important role, but is affected by low spatial resolution. The widespread forest coverage area in the world makes it a good choice to use remote sensing technology to measure AGB and CH. The working band (L band) signal of the Global Navigation Satellite System-Reflectometry (GNSS-R) has good penetration ability, which makes the GNSS-R signal have the advantage of being able to penetrate the atmosphere, vegetation, etc., and help to perceive the soil under the vegetation. In addition, the passive space-borne GNSS-R receiver makes the use of this technology to invert AGB and CH have great cost advantages. The attenuation of soil near-specular reflection signals caused by vegetation is the main mechanism by which GNSS-R data are sensitive to biomass. Therefore, surface reflectivity is a key observation value for inverting AGB and can be calculated from the delay-Doppler map (DDM).

[0003] The development of computer technology has promoted the rapid development of machine learning. For related research using machine learning or deep learning, the inversion results of AGB and CH using satellite GNSS-R technology are relatively good. However, in these studies, the importance of features is not explained in depth, and it is often difficult to fully explore and explain the complex relationship between influencing factors, resulting in redundant and complicated feature selection or even inaccurate features, which in turn affects the inversion effect and prediction accuracy of the model. The importance of features is one of the key indicators to understand how the model makes predictions. Therefore, the explanation of feature importance is particularly important. In order to solve the limitations of existing machine learning-based satellite GNSS-R feature selection and forest aboveground biomass and canopy height estimation methods, it is urgent to overcome the following problems: insufficient feature explanation, redundant and complicated features, weak model interpretability, and low inversion efficiency and accuracy. SHAP has an excellent performance in explaining feature importance. In addition, the computational time cost of model training during modeling is also an important indicator for judging the quality of the model. Compared with other ensemble learning algorithms, CatBoost, XGBoost and LightGBM, three tree model algorithms based on boosted decision trees, have significant advantages in model training time. In tree models, feature importance is usually measured by the splitting node of the feature. By observing the contribution of the feature in the model, we can have a deeper understanding of the important features in the dataset.

[0004] In view of this, the present invention is proposed. Summary of the invention

[0005] In order to overcome the shortcomings of the existing satellite GNSS-R feature selection and forest aboveground biomass and canopy height estimation methods based on machine learning, the present invention proposes a satellite GNSS-R important feature selection and forest aboveground biomass and canopy height estimation method combining a tree model ensemble learning algorithm and a SHAP interpreter to address the problems of insufficient feature interpretation, complex feature redundancy, weak model interpretability, and low inversion efficiency and accuracy in the current methods. The method combines three tree model algorithms, CatBoost, XGBoost and LightGBM, and the SHAP interpreter to combine and optimize the input variable feature parameters in the satellite GNSS-R forest aboveground biomass and canopy height estimation model, and uses the selected optimal feature parameter combination to construct a forest aboveground biomass and canopy height estimation model based on the three tree model algorithms. This method can well explain the important features in machine learning, improve the model training speed, and has good inversion accuracy.

[0006] In order to solve the above technical problems, the basic concept of the technical solution adopted by the present invention is:

[0007] A method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter includes the following steps:

[0008] Step S1, download and collect eight CYGNSS satellite L1-level observation data and SMAP, NOAA 1km resolution DEM data, European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDI LIDAR data, and canopy height (CH) data integrating GEDI and Landsat data from relevant websites;

[0009] Step S2, the AGB and CH datasets are gridded at a spatial resolution of 5 km × 5 km and spatially matched with the CYGNSS satellite L1 level observation data. At the same time, the CYGNSS data is quality controlled to filter out observations that seriously affect the reflection signal;

[0010] Step S3, extracting variable parameters from the L1 level observation data of 8 CYGNSS satellites and calculating the surface reflectivity using the CYGNSS satellite GNSS-R data;

[0011] Step S4, divide the data set into two parts, 60% of the total data set as training set and 40% of the total data set as

[0012] As a test set;

[0013] Step S5, combining the CatBoost, XGBoost and LightGBM algorithms and the SHAP interpreter to respectively optimize the input variable characteristic parameters in the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model, and using the selected optimal characteristic parameter combination to construct the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model based on the three tree model algorithms;

[0014] Step S6, using the European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDILIDAR data, and canopy height (CH) data integrating GEDI and Landsat data as references, the optimal variable characteristic parameters selected in the test data set are input as input to the trained space-borne GNSS-R forest aboveground biomass and canopy height estimation model to calculate the GNSS-R AGB and CH estimates, and finally, the GNSS-R AGB and CH estimates are compared with the reference data, and the root mean square error (RMSE), mean absolute error (MAE) and Pearson correlation coefficient (CC) are used as indicators to evaluate the inversion performance of the model to verify the accuracy of the model in estimating AGB and CH.

[0015] Furthermore, the SMAP data used in step S1 specifically include soil moisture variables (SM), surface temperature (SST), surface roughness coefficient (RC), and the spatial resolution and temporal resolution are: 36km×36km and 1 day respectively; the spatial resolution and temporal resolution of the European Space Agency (ESA) CCI biomass (AGB) map used are 100m and 1 year respectively; the spatial resolution and temporal resolution of the LUCID AGB data used are 1km and 1 year respectively; the spatial resolution and temporal resolution of the GEDI canopy height (CH) based on GEDI LIDAR data used are 0.5° and 1 year respectively, and the spatial resolution and temporal resolution of the canopy height (CH) data integrating GEDI and Landsat data are 30m and 1 year respectively.

[0016] Furthermore, in step S2, the AGB and CH data sets are unified onto a 5km x 5km spatial resolution grid and spatially matched with the CYGNSS satellite L1 observation data, and the CYGNSS data are quality controlled by using the quality flag (quality_flags) variable provided in the CYGNSS GNSSS-R L1 data set to filter the data, i.e., setting 12 quality flags (S BandPowered Up, Large Spacecraft Attitude Error, Blackbody DDM, DDMI Reconfigured, Spacewire CRC Invalid, DDM is Test Patten, Channel Idle, Direct Signal in DDM, Low Confidence GPS EIRP Estimate, RFIDetected, SP Non Existent Error, and Blackbody Framing Error) to "Yes". In addition, since the incident angle affects the reflected signal, the observation data with an incident angle of less than 10° and >60° needs to be filtered.

[0017] Furthermore, the variable parameters extracted from the L1-level observation data of the eight CYGNSS satellites in step S3 include respectively extracting the surface reflectivity (reflectivity) of each mirror reflection point and extracting the SNR (ddm_snr) at the mirror reflection point, the mirror reflection point position, the incident angle, the range correction gain (RCG), the spacecraft attitude roll angle (sc_roll), the spacecraft attitude pitch angle (sc_pitch), the GPS satellite PRN code (prn_code) for generating the DDM waveform, the GPS spacecraft number (sv_num), the DDM track ID (track_id), the RCG performance index (FOM) of DDM (prn_fig_of_merit), the mirror reflection point longitude (sp_lon), the mirror reflection point latitude (sp_lat), the mirror reflection point altitude (sp_alt), the mirror reflection point incident angle (sp_inc_angle), and the mirror reflection point in the orbital coordinate system. Altitude angle (sp_theta_orbit), azimuth angle of mirror reflection point in orbital coordinate system (sp_az_orbit), altitude angle of mirror reflection point in satellite body coordinate system (sp_theta_body), azimuth angle of mirror reflection point in satellite body coordinate system (sp_az_body), antenna gain from receiver to mirror reflection point (sp_rx_gain), DDM noise floor (ddm_noise_floor), distance from receiver to mirror point (rx_to_sp_range), distance from transmitter to mirror point (tx_to_sp_range), Fresnel power reflection coefficient at mirror point (fresnel_coeff), normalized BRCS (ddm_nbrcs), leading edge slope (ddm_les), NBRCS scattering area (nbrcs_scatter_area), LES scattering area (les_scatter_area), DDM BRCS uncertainty (ddm_brcs_uncert), UTC time corresponding to the middle moment of DDM acquisition (ddm_timestamp_utc), GPS SV transmit power (gps_tx_power_db_w), GPS SV transmit antenna gain (gps_ant_gain_db_i), GPS off-axis angle (gps_off_boresight_angle_deg) and equivalent isotropic radiated power (gps_eirp) of GPS transmitter observed variables, and bistatic radar cross section (BRCS) and effective scattering area (eff_scatter) were also extracted. The calculation formula for calculating the surface reflectivity using CYGNSS satellite GNSS-R data is as follows:

[0018]

[0019] Where P DDM is the peak value of the simulated scattered power of DDM; P t G t is the equivalent isotropic radiated power of the transmitter; P t is the transmission power of the right-hand circularly polarized navigation signal of the GNSS satellite; G t is the transmitting antenna gain; G r is the receiver antenna gain; λ is the carrier wavelength of the transmitted signal; R t and R r are the distance from the transmitter to the ground launch point and the distance from the receiver to the ground launch point, respectively. N is the DDM noise floor.

[0020] Furthermore, the data set is divided in step S4 by using 60% of the data set after time-space matching as a training set and 40% of the data set as a test set;

[0021] Furthermore, the combination of CatBoost, XGBoost and LightGBM algorithms and the SHAP interpreter described in step S5 respectively optimizes the input variable characteristic parameters in the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model, and uses the selected optimal characteristic parameter combination to construct a spaceborne GNSS-R forest aboveground biomass and canopy height estimation model based on three tree model algorithms;

[0022] Further, the European Space Agency (ESA) CCI biomass (AGB) map, LUCIDAGB, GEDI canopy height (CH) based on GEDI LIDAR data, and canopy height (CH) data integrated with GEDI and Landsat data are used as references in step S6, wherein the AGB reference value finally used for AGB estimation model construction and testing is the average value of the European Space Agency (ESA) CCI biomass (AGB) map and LUCID AGB, and the CH reference value finally used for CH estimation model construction and testing is the average value of GEDI canopy height (CH) based on GEDI LIDAR data and canopy height (CH) integrated with GEDI and Landsat data. The optimal variable characteristic parameters selected from the test data set are used as inputs and input into the trained spaceborne GNSS-R forest aboveground biomass and canopy height estimation model to calculate the GNSS-R AGB and CH estimation values, and finally the GNSS-R AGB and CH estimation values ​​are compared with the reference data to verify the accuracy of the model estimation of AGB and CH. The root mean square error (RMSE), mean absolute error (MAE) and Pearson correlation coefficient (CC) are used as indicators to evaluate the inversion performance of the model. The relevant calculation formulas are as follows:

[0023]

[0024]

[0025] Where m is the number of data samples, y i,M and i,T are the AGB (or CH) estimated by the model and the AGB (or CH) obtained from the reference data, and They are y i,M and i,T The average value of .

[0026] The present invention combines the tree model ensemble learning algorithm and the satellite-borne GNSS-R data of the SHAP interpreter, and provides a new acquisition method through the collaborative inversion of multivariate parameters for global high spatiotemporal resolution of forest aboveground biomass and canopy height, and effectively considers multiple influencing factors. Through more precise feature selection and interpretation mechanisms, the complex relationship between influencing factors is fully explored, the interpretability and prediction accuracy of the model are improved, and thus the estimation efficiency and accuracy of forest aboveground biomass and canopy height are significantly improved.

[0027] After adopting the above technical scheme, the present invention has the following beneficial effects compared with the prior art.

[0028] The method of the present invention combines three efficient gradient boosting decision tree algorithms, CatBoost, XGBoost and LightGBM, and SHAP interpreter to select important features for modeling in biomass and canopy height inversion, and builds a global high spatiotemporal forest aboveground biomass and canopy height estimation model based on three tree model ensemble learning algorithms, CatBoost, XGBoost and LightGBM, integrating bistatic radar cross section (BRCS), effective scattering area (Effective Scattering Area) and selected important features. The present invention is the first to realize the collaborative inversion of global high spatiotemporal forest aboveground biomass and canopy height by combining satellite-borne GNSS-R data and multivariate parameters with tree model ensemble learning algorithm and SHAP interpreter, providing a new means for obtaining forest aboveground biomass and canopy height; and realizes a method for inverting forest aboveground biomass and canopy height by satellite-borne GNSS-R taking into account multiple influencing factors.

[0029] The specific implementation modes of the present invention are further described in detail below in conjunction with the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] The accompanying drawings are part of this application and are used to provide a further understanding of the present invention. The illustrative embodiments of the present invention and their descriptions are used to explain the present invention, but do not constitute an improper limitation of the present invention. Obviously, the drawings described below are only some embodiments. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work. In the drawings:

[0031] Figure 1 Schematic diagram of the method flow in an embodiment of the present invention.

[0032] Figure 2 This is a calculation diagram of the importance of forest aboveground biomass (AGB) inversion features in an embodiment of the present invention.

[0033] Figure 3 Calculation diagram of the importance of canopy height (CH) inversion features in an embodiment of the present invention.

[0034] FIG. 4( a ), FIG. 4( b ), and FIG. 4( c ) are relationship diagrams of important features and target values ​​selected in the inversion of forest aboveground biomass (AGB) in an embodiment of the present invention.

[0035] FIG. 5( a ), FIG. 5( b ), and FIG. 5( c ) are relationship diagrams of important features and target values ​​selected in canopy height (CH) inversion in an embodiment of the present invention.

[0036] FIG. 6( a ) and FIG. 6( b ) are comparison diagrams of the importance of inversion features of forest aboveground biomass (AGB) and canopy height (CH) in an embodiment of the present invention.

[0037] Figure 7 This is a SHAP result diagram of forest biomass (AGB) inverted by the CatBoost model in an embodiment of the present invention.

[0038] Figure 8 This is a SHAP result diagram of canopy height (CH) inverted by the CatBoost model in an embodiment of the present invention.

[0039] Fig. 9 Graphs showing the inversion results of the aboveground biomass (AGB) of the forests of various models in the embodiments of the present invention (top: 35 input characteristic parameters, bottom: 8 input characteristic parameters).

[0040] Fig.10 Graphs showing the canopy height (CH) inversion results of various models in the embodiments of the present invention (top: 35 input feature parameters, bottom: 8 input feature parameters).

[0041] It should be noted that these drawings and textual descriptions are not intended to limit the conceptual scope of the present invention in any way, but are intended to illustrate the concept of the present invention for those skilled in the art by referring to specific embodiments. DETAILED DESCRIPTION

[0042] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. The following embodiments are used to illustrate the present invention but are not used to limit the scope of the present invention.

[0043] Example 1

[0044] In order to verify the effectiveness of the method proposed in the present invention, the 2021 CYGNSS LIB data, the SMAP data with a spatial resolution of 36 km×36 km and a time resolution of 1 day, the European Space Agency (ESA) CCI biomass (AGB) map with a spatial resolution of 100 m and a time resolution of 1 year, the LUCID AGB data with a spatial resolution of 1 km and a time resolution of 1 year, the GEDI canopy height (CH) based on GEDI LIDAR data with a spatial resolution of 0.5° and a time resolution of 1 year, and the canopy height (CH) data of integrated GEDI and Landsat data with a spatial resolution of 30 m and a time resolution of 1 year were selected. Among them, the average AGB value of the European Space Agency (ESA) CCI biomass (AGB) map and the LUCID AGB, the GEDI canopy height (CH) based on GEDI LIDAR data, and the average CH value of the canopy height (CH) data integrating GEDI and Landsat data were used as reference data, and the root mean square error (RMSE), mean absolute error (MAE) and Pearson correlation coefficient (CC) were used as indicators to evaluate the inversion performance of the model.

[0045] The basic configuration of the experimental platform of this embodiment is shown in Table 1:

[0046] Table 1 Configuration of the experimental platform

[0047]

[0048] A method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter. The technical solution implementation process is as follows: Figure 1 As shown, the following steps are included:

[0049] Step S1, download and collect eight CYGNSS satellite L1-level observation data and SMAP, NOAA 1km resolution DEM data, European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDI LIDAR data, and canopy height (CH) data integrating GEDI and Landsat data from relevant websites;

[0050] Step S2, the AGB and CH datasets are gridded at a spatial resolution of 5 km × 5 km and spatially matched with the CYGNSS satellite L1 level observation data. At the same time, the CYGNSS data is quality controlled to filter out observations that seriously affect the reflection signal;

[0051] Step S3, extracting variable parameters from the L1 level observation data of 8 CYGNSS satellites and calculating the surface reflectivity using the CYGNSS satellite GNSS-R data;

[0052] Step S4, divide the data set into two parts, 60% of the total data set as training set and 40% of the total data set as

[0053] As a test set;

[0054] Step S5, combining the CatBoost, XGBoost and LightGBM algorithms and the SHAP interpreter to respectively optimize the input variable characteristic parameters in the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model, and using the selected optimal characteristic parameter combination to construct the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model based on the three tree model algorithms;

[0055] Step S6, using the European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDILIDAR data, and canopy height (CH) data integrating GEDI and Landsat data as references, the optimal variable feature parameters selected in the test data set are input as input to the trained space-borne GNSS-R forest aboveground biomass and canopy height estimation model to calculate the GNSS-R AGB and CH estimates, and finally, the GNSS-R AGB and CH estimates are compared with the reference data to verify the accuracy of the model in estimating AGB and CH.

[0056] As an implementation mode of this embodiment, in step S1, the SMAP data used specifically include three auxiliary data of soil moisture (SM) variable, surface temperature (SST), and surface roughness coefficient (RC), and the spatial resolution and temporal resolution are: 36km×36km and 1 day respectively; the spatial resolution and temporal resolution of the European Space Agency (ESA) CCI biomass (AGB) map used are 100m and 1 year respectively; the spatial resolution and temporal resolution of the LUCID AGB data used are 1km and 1 year respectively; the spatial resolution and temporal resolution of the GEDI canopy height (CH) based on GEDI LIDAR data used are 0.5° and 1 year respectively, and the spatial resolution and temporal resolution of the canopy height (CH) data integrating GEDI and Landsat data are 30m and 1 year respectively.

[0057] As an implementation of this embodiment, in step S2, the AGB and CH data sets are unified on a 5km x 5km spatial resolution grid and are spatiotemporally matched with the CYGNSS satellite L1-level observation data, and the CYGNSS data are quality controlled. Specifically, the quality flag (quality_flags) variable provided in the CYGNSS GNSSS-R L1 data set is used for data filtering, that is, 12 quality flags (S Band PoweredUp, LargeSpacecraft Attitude Error, Blackbody DDM, DDMI Reconfigured, Spacewire CRCInvalid, DDM is Test Patten, Channel Idle, Direct Signal in DDM, Low ConfidenceGPS EIRP Estimate, RFIDetected, SP Non Existent Error, and Blackbody Framing Error) are set to "Yes". In addition, since the incident angle affects the reflected signal, the observation data with an incident angle less than 10° and >60° needs to be filtered.

[0058] As an implementation mode of this embodiment, in step S3, the variable parameters extracted from the L1-level observation data of the eight CYGNSS satellites include respectively extracting the surface reflectivity (reflectivity) of each mirror reflection point and extracting the SNR (ddm_snr) at the mirror reflection point, the mirror reflection point position, the incident angle, the range correction gain (RCG), the spacecraft attitude roll angle (sc_roll), the spacecraft attitude pitch angle (sc_pitch), the GPS satellite PRN code (prn_code) for generating the DDM waveform, the GPS spacecraft number (sv_num), the DDM track ID (track_id), the RCG performance index (FOM) of the DDM (prn_fig_of_merit), the longitude of the mirror reflection point (sp_lon), the latitude of the mirror reflection point (sp_lat), the height of the mirror reflection point (sp_alt), the incident angle of the mirror reflection point (sp_inc_angle), and the mirror reflection point on track. The altitude angle in the orbit coordinate system (sp_theta_orbit), the azimuth angle of the mirror reflection point in the orbit coordinate system (sp_az_orbit), the altitude angle of the mirror reflection point in the satellite body coordinate system (sp_theta_body), the azimuth angle of the mirror reflection point in the satellite body coordinate system (sp_az_body), the antenna gain from the receiver to the mirror reflection point (sp_rx_gain), the DDM floor noise (ddm_noise_floor), the distance from the receiver to the mirror point (rx_to_sp_range), the distance from the transmitter to the mirror point (tx_to_sp_range), the Fresnel power reflection coefficient at the mirror point (fresnel_coeff), the normalized BRCS (ddm_nbrcs), the leading edge slope (ddm_les), the NBRCS scattering area (nbrcs_scatter_area), the LES scattering area (les_scatter_area), the DDM BRCS uncertainty (ddm_brcs_uncert), UTC time corresponding to the middle moment of DDM acquisition (ddm_timestamp_utc), GPS SV transmit power (gps_tx_power_db_w), GPS SV transmit antenna gain (gps_ant_gain_db_i), GPS off-axis angle (gps_off_boresight_angle_deg) and GPS transmitter equivalent isotropic radiated power (gps_eirp) observation variables, and bistatic radar cross section (BRCS) and effective scattering area (eff_scatter) are also extracted. Figure 2 and attached Figure 3The calculation of the inversion feature importance of different input variables of the three models LightGBM, XGBoost and CatBoost for AGB and CH is shown respectively. Figures 4 and 5 show the relationship between the important features selected in the inversion of AGB and CH and the target values, respectively. Finally, the comparison of the feature importance of inverting AGB (right) and CH (left) by different inversion models (i.e., XGBoost, LightGBM and CatBoost) is shown in Figure 6. Figure 6 shows the importance ranking of different input variables in the inversion of AGB (right) and CH (left). It can be seen that there are differences in the feature importance of different inversion models, and through the connection between the nodes, it can be intuitively seen how the importance of each feature changes in different models. For AGB, in general, the longitude of the specular reflection point (sp_lon), latitude of the specular reflection point (sp_lat), altitude of the specular reflection point (sp_alt), soil moisture (SM), surface roughness coefficient (RC) and surface reflectivity (reflectivity) are all important parameters in the three inversion models. For CH, compared with AGB, the importance of the features of the three inversion models is more different. The most important features of XGBoost are soil moisture (SM), surface reflectivity (reflectivity), front slope (ddm_les), specular reflection (sp_alt), soil moisture (SM), surface roughness coefficient (RC) and surface reflectivity (reflectivity). The most important features of LightGBM are the longitude of the mirror reflection point (sp_lon), the latitude of the mirror reflection point (sp_lat), the height of the mirror reflection point (sp_alt), the surface roughness coefficient (RC), and the soil moisture (SM); the most important features of CatBoost are the latitude of the mirror reflection point (sp_lat), the soil moisture (SM), the longitude of the mirror reflection point (sp_lon), the height of the mirror reflection point (sp_alt), and the surface reflectivity (reflectivity). The difference in the importance of the features of different inversion models indicates that when different parameters are selected to train the model, there will be different inversion accuracies. Therefore, the selection of feature parameters for different models is very important.

[0059] The calculation formula for calculating the surface reflectivity using CYGNSS satellite GNSS-R data is as follows:

[0060]

[0061] Where P DDM is the peak value of the simulated scattered power of DDM; P t G t is the equivalent isotropic radiated power of the transmitter; P t is the transmission power of the right-hand circularly polarized navigation signal of the GNSS satellite; G tis the transmitting antenna gain; G r is the receiver antenna gain; λ is the carrier wavelength of the transmitted signal; R t and R r are the distance from the transmitter to the ground launch point and the distance from the receiver to the ground launch point, respectively. N is the DDM noise floor.

[0062] As an implementation of this embodiment, in step S4, the data set is divided into 60% of the data set after spatiotemporal matching as a training set and 40% of the data set as a test set;

[0063] As an implementation method of this embodiment, in step S5, the combination of CatBoost, XGBoost and LightGBM algorithms and SHAP interpreter respectively optimizes the input variable feature parameters in the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model, and uses the selected optimal feature parameter combination to construct a spaceborne GNSS-R forest aboveground biomass and canopy height estimation model based on three tree model algorithms. In order to show which features are most important to the model, the present invention is attached Figure 7 and attached Figure 8 The SHAP results of the CatBoost model inversion of forest aboveground biomass (AGB) and canopy height (CH) are shown respectively, where the SHAP values ​​(upper left) and average absolute total contributions (upper middle) of all features of each sample are shown. At the same time, in order to observe the contribution of different features to the output, the SHAP values ​​of AGB and CH are also shown (upper right). The length of the bar represents the contribution of each feature. The expected value of AGB for the CatBoost model is 16.723. The SHAP value for the AGB sample is more negative, and the inversion rate is slightly reduced. Among them, sp_alt is the largest in this sample, indicating that it has a great influence on the inversion. The expected value of CH for the CatBoost model is 6.528. The SHAP value for the CH sample is also more negative, and the inversion rate is also slightly reduced. Among them, sp_lat is the largest in this sample, indicating that it has a great influence on the inversion.

[0064] In addition, Figure 7 and attached Figure 8 The lower part of the figure is a SHAP force diagram, which can visualize the influence of features. In this figure, we show the features that push up the prediction in red and the features that push down the prediction in blue. You can intuitively see the eigenvalue of each feature, and by combining the influence of all features, you can intuitively see how the model arrives at the final prediction result.

[0065] As an implementation method of this embodiment, in step S6, the European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDI LIDAR data, and canopy height (CH) data integrated with GEDI and Landsat data are used as references respectively. Specifically, the AGB reference value finally used for AGB estimation model construction and testing is the average value of the European Space Agency (ESA) CCI biomass (AGB) map and LUCID AGB, and the CH reference value finally used for CH estimation model construction and testing is the average value of GEDI canopy height (CH) based on GEDI LIDAR data and canopy height (CH) integrated with GEDI and Landsat data. The optimal variable characteristic parameters selected from the test data set are used as inputs and input into the trained spaceborne GNSS-R forest aboveground biomass and canopy height estimation model to calculate the GNSS-R AGB and CH estimation values, and finally the GNSS-R AGB and CH estimation values ​​are compared with the reference data to verify the accuracy of the model estimation of AGB and CH. The root mean square error (RMSE), mean absolute error (MAE) and Pearson correlation coefficient (CC) are used as indicators to evaluate the inversion performance of the model. The relevant calculation formulas are as follows:

[0066]

[0067] Where m is the number of data samples, y i,M and i,T are the AGB (or CH) estimated by the model and the AGB (or CH) obtained from the reference data, and They are y i,M and i,T The average value of .

[0068] Attached Fig. 9 and attached Fig.10 The inversion results of forest aboveground biomass (AGB) and canopy height (CH) of each model are shown in Figure 2 (top: input characteristic parameters are 35, bottom: input characteristic parameters are 8). Fig. 9 and attached Fig.10 The following conclusions can be drawn:

[0069] (1) By Fig. 9It can be seen that among the three models of LightGBM, XGBoost and CatBoost, no matter the input feature parameters are 35 or 8, the CatBoost model has the best inversion performance in terms of AGB. When the input feature parameters are 35, the RMSE, MAE and CC are 9.421 ton / ha, 5.128 ton / ha and 0.889 respectively; when the input feature parameters are 8, the RMSE, MAE and CC are 8.341 ton / ha, 4.609 ton / ha and 0.915 respectively.

[0070] (2) By the attached Fig.10 It can be seen that among the three models of LightGBM, XGBoost and CatBoost, whether the input feature parameters are 35 or 8, the CatBoost model has the best inversion performance in terms of CH. When the input feature parameters are 35, the RMSE, MAE and CC are 0.879m, 0.595m and 0.948 respectively; when the input feature parameters are 8, the RMSE, MAE and CC are 0.784m, 0.533m and 0.959 respectively.

[0071] (3) By the attached Fig. 9 and attached Fig.10 It can be seen that the inversion performance of selecting 8 eigenvalues ​​can improve the inversion performance of 35 observations to some extent. In terms of inversion of AGB, when the input feature parameters are reduced from 35 to 8, the RMSE, MAE and CC obtained by the three models of LightGBM, XGBoost and CatBoost have been slightly improved on the basis of the previous ones; in terms of inversion of CH, when the input feature parameters are reduced from 35 to 8, except for the LightGBM model, the RMSE, MAE and CC obtained by the two models of XGBoost and CatBoost have also been slightly improved on the basis of the previous ones. It can be concluded that, to a certain extent, the inversion results of inputting 8 eigenvalues ​​are significantly better than the inversion results of inputting 35 eigenvalues. At the same time, fewer input parameters can speed up the training speed of the model and improve the generalization ability of the model to a certain extent.

[0072] In summary, the inversion performance of the model can be significantly improved by accurately selecting the important features in the model using the implementation scheme of the present invention, especially in combination with tree model ensemble learning and with the help of the SHAP interpreter tool. The advantage of this method is that it can better understand the impact of features on model output, thereby optimizing the feature selection process, reducing redundant features, and enhancing the interpretability and accuracy of the model. Effective feature selection can not only improve the predictive ability of the model, but also speed up the training speed of the model and improve its generalization ability, especially in the inversion tasks of forest aboveground biomass and canopy height, which will directly lead to performance improvement. It also further shows that selecting important modeling features helps to improve the inversion performance of the model.

[0073] The above is only a preferred embodiment of the present invention, and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment, it is not used to limit the present invention. Any technician familiar with this patent can make some changes or modify the technical contents suggested above into equivalent embodiments without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the solution of the present invention.

Claims

1. A method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter, characterized in that: The following steps are involved: Step S1, download and collect eight CYGNSS satellite L1-level observation data and SMAP, NOAA 1km resolution DEM data, European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDI LIDAR data, and canopy height (CH) data integrating GEDI and Landsat data from relevant websites; Step S2, the AGB and CH datasets are gridded at a spatial resolution of 5 km × 5 km and spatially matched with the CYGNSS satellite L1 level observation data. At the same time, the CYGNSS data is quality controlled to filter out observations that seriously affect the reflection signal; Step S3, extracting variable parameters from the L1 level observation data of 8 CYGNSS satellites and calculating the surface reflectivity using the CYGNSS satellite GNSS-R data; Step S4, dividing the data set, 60% of the total data set as a training set, and 40% of the total data set as a test set; Step S5, combining the CatBoost, XGBoost and LightGBM algorithms and the SHAP interpreter to respectively optimize the input variable characteristic parameters in the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model, and using the selected optimal characteristic parameter combination to construct the spaceborne GNSS-R forest aboveground biomass and canopy height estimation model based on the three tree model algorithms; Step S6, using the European Space Agency (ESA) CCI biomass (AGB) map, LUCID AGB, GEDI canopy height (CH) based on GEDI LIDAR data, and canopy height (CH) data integrating GEDI and Landsat data as references, the optimal variable characteristic parameters selected in the test data set are input as input to the trained space-borne GNSS-R forest aboveground biomass and canopy height estimation model to calculate the GNSS-R AGB and CH estimates, and finally, the GNSS-R AGB and CH estimates are compared with the reference data, and the root mean square error (RMSE), mean absolute error (MAE) and Pearson correlation coefficient (CC) are used as indicators to evaluate the inversion performance of the model to verify the accuracy of the model in estimating AGB and CH.

2. The method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter according to claim 1, characterized in that: The SMAP data used in step S1 specifically include three auxiliary data: soil moisture (SM) variables, surface temperature (SST), and surface roughness coefficient (RC), and the spatial resolution and temporal resolution are: 36km×36km and 1 day, respectively; the spatial resolution and temporal resolution of the European Space Agency (ESA) CCI biomass (AGB) map used are 100m and 1 year, respectively; the spatial resolution and temporal resolution of the LUCID AGB data used are 1km and 1 year, respectively; the spatial resolution and temporal resolution of the GEDI canopy height (CH) based on GEDI LIDAR data are 0.5° and 1 year, respectively, and the spatial resolution and temporal resolution of the canopy height (CH) data integrating GEDI and Landsat data are 30m and 1 year, respectively.

3. The method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter according to claim 1, characterized in that: In step S2, the specific method for quality control of CYGNSS data is to use the quality flag (quality_flags) variable provided in the CYGNSS GNSSS-R L1 dataset to filter the data, that is, set 12 quality flags (S Band Powered Up, LargeSpacecraft Attitude Error, Blackbody DDM, DDMI Reconfigured, Spacewire CRCInvalid, DDM is Test Patten, Channel Idle, Direct Signal in DDM, Low ConfidenceGPS EIRP Estimate, RFIDetected, SP Non Existent Error, and Blackbody Framing Error) to "Yes"; in addition, the observation data with incident angles less than 10° and greater than 60° need to be filtered.

4. The method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter according to claim 1, characterized in that: In step S3, the variable parameters extracted from the L1 observation data of the eight CYGNSS satellites include the surface reflectivity (reflectivity) of each mirror reflection point, the SNR (ddm_snr) at the mirror reflection point, the range correction gain (RCG), the spacecraft attitude roll angle (sc_roll), the spacecraft attitude pitch angle (sc_pitch), the GPS satellite PRN code (prn_code) for generating the DDM waveform, the GPS spacecraft number (sv_num), the DDM track ID (track_id), the RCG performance index (FOM) of DDM (prn_fig_of_merit), the longitude of the mirror reflection point (sp_lon), the latitude of the mirror reflection point (sp_lat), the height of the mirror reflection point (sp_alt), the incident angle of the mirror reflection point (sp_inc_angle), and the altitude angle of the mirror reflection point in the orbital coordinate system (sp_th eta_orbit), azimuth of the mirror reflection point in the orbital coordinate system (sp_az_orbit), altitude of the mirror reflection point in the satellite body coordinate system (sp_theta_body), azimuth of the mirror reflection point in the satellite body coordinate system (sp_az_body), antenna gain from the receiver to the mirror reflection point (sp_rx_gain), DDM floor noise (ddm_noise_floor), distance from the receiver to the mirror point (rx_to_sp_range), distance from the transmitter to the mirror point (tx_to_sp_range), Fresnel power reflection coefficient at the mirror point (fresnel_coeff), normalized BRCS (ddm_nbrcs), leading edge slope (ddm_les), NBRCS scattering area (nbrcs_scatter_area), LES scattering area (les_scatter_area), DDM BRCS uncertainty (ddm_brcs_uncert), UTC time corresponding to the middle moment of DDM acquisition (ddm_timestamp_utc), GPS SV transmit power (gps_tx_power_db_w), GPS SV transmit antenna gain (gps_ant_gain_db_i), GPS off-axis angle (gps_off_boresight_angle_deg) and equivalent isotropic radiated power of GPS transmitter (gps_eirp) observation variables. Bistatic radar cross section (BRCS) and effective scattering area (eff_scatter) are also extracted. The calculation formula for calculating the surface reflectivity using CYGNSS satellite GNSS-R data is as follows: Where P DDM is the peak value of the simulated scattered power of DDM; P t G t is the equivalent isotropic radiated power of the transmitter; P t is the transmission power of the right-hand circularly polarized navigation signal of the GNSS satellite; G t is the transmitting antenna gain; G r is the receiver antenna gain; λ is the carrier wavelength of the transmitted signal; R t and R r are the distance from the transmitter to the ground launch point and the distance from the receiver to the ground launch point, respectively. N is the DDM noise floor.

5. The method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter according to claim 1, characterized in that: In step S4, the data set is divided into 60% of the data set after time-space matching as a training set and 40% of the data set as a test set.

6. The method for selecting important features of spaceborne GNSS-R and estimating forest aboveground biomass and canopy height by combining a tree model ensemble learning algorithm and a SHAP interpreter according to claim 1, characterized in that: In step S6, the root mean square error (RMSE), mean absolute error (MAE) and Pearson correlation coefficient (CC) are used as indicators to evaluate the inversion performance of the model. The relevant calculation formulas are as follows: Where m is the number of data samples, y i,M and i,T are the AGB (or CH) estimated by the model and the AGB (or CH) obtained from the reference data, and They are y i,M and i,T The average value of .

Images