Large-scale forest carbon storage estimation method and system based on active and passive remote sensing technology

Abstract

The application discloses a large-scale forest land carbon storage estimation method and system based on active and passive remote sensing technology, which comprises the following steps: acquiring satellite-borne laser radar data of different sensors, L-band and C-band SAR data, multispectral remote sensing image data, airborne laser radar data, terrain data, and land cover category data containing forest land; obtaining discontinuous joint satellite-borne laser canopy height products by using satellite-borne laser radar data of different sensors; improving a geographic weighted regression model to estimate continuous scale canopy height; extracting polarization parameters by using C-band SAR data, extracting vegetation index and texture index by using multispectral remote sensing image data, and constructing feature variables together with terrain data; screening the feature variables based on the variance reduction criterion; and constructing an above-ground carbon storage inversion model to estimate large-scale forest land carbon storage. The application improves the estimation precision of forest land carbon storage in large-scale broken plots.

Description

Technical Field

[0001] The present invention relates to the technical fields of remote sensing, geographic information and ecology, and in particular to a large-scale forest carbon storage estimation method and system based on active and passive remote sensing technologies. Background Art

[0002] Carbon storage, or the amount of carbon stored, typically refers to the amount of carbon in a carbon reservoir (forests, oceans, land, etc.). It is an important indicator of an ecosystem's ability to accumulate organic matter over a specific period of time and a key parameter for describing its structure and function. With the intensification of global climate change and the increasing prominence of environmental issues, the frequent occurrence of extreme weather events threatens the stability and service functions of ecosystems. Therefore, more accurate estimation of vegetation carbon storage over large regions has become a crucial measure of regional vegetation carbon sequestration capacity and the value of its ecological services, as well as a crucial basis for governments to formulate environmental protection and sustainable development strategies.

[0003] Currently, methods for estimating vegetation carbon storage are mainly divided into direct weighing, allometric growth models based on tree parameters, and remote sensing estimation. Direct weighing directly obtains carbon storage values ​​by sampling vegetation, drying it, and weighing it. Allometric growth models based on tree parameters use empirical formulas to estimate carbon storage using parameters such as tree diameter at breast height and height. Remote sensing estimation uses visible light imagery, synthetic aperture radar (SAR) data, and light detection and ranging (LiDAR) data to calculate relevant spectral indices, texture indices, and canopy indices. Combined with ground-based carbon storage data, carbon storage is inverted using regression models.

[0004] LiDAR data can proactively acquire high-resolution vegetation canopy information and plays a crucial role in carbon stock estimation. LiDAR data can be acquired from three platforms: ground-based, airborne, and satellite-based. Ground-based and airborne LiDARs offer high spatial resolution and data accuracy, but their applicability is typically limited due to sampling range and acquisition costs. The development of satellite-based LiDARs (such as ICESat-2 and GEDI) has overcome the limited coverage of traditional platforms, providing canopy information over large areas and enabling accurate and large-scale estimation of vegetation carbon stocks. However, due to different sensor designs and data processing algorithms, ICESat-2 and GEDI satellite-based LiDARs differ in observation accuracy, coverage, and observation time. Therefore, it is crucial to combine these two data types to leverage the strengths of different satellite-based LiDAR data, calibrate and optimize relevant canopy indices (such as canopy height models), and thus achieve more accurate and large-scale carbon stock estimation. Summary of the Invention

[0005] In view of the shortcomings of the existing technology, the present invention provides a large-scale forest carbon storage estimation method based on active and passive remote sensing technology, which includes the following steps:

[0006] Step 1: Obtain spaceborne lidar data, L-band and C-band SAR data, multispectral remote sensing image data, airborne lidar data, terrain data, and land cover category data including forest land from different sensors in the study area;

[0007] Step 2: preprocess and weighted merge the spaceborne lidar data from different sensors to obtain a discontinuous joint spaceborne lidar canopy height product;

[0008] Step 3: Multiple interaction terms between ground elevation and synthetic aperture radar data were introduced to improve the geographically weighted regression model. The coefficients of the improved geographically weighted regression model were obtained using the least squares method, and then the canopy height of the study area was estimated at a continuous scale.

[0009] Step 4: Use C-band SAR data to extract polarization parameters, and multispectral remote sensing image data to extract vegetation index and texture index. Together with the terrain data, these are used to construct feature variables. The feature variables are trained using a random forest model and screened based on the variance reduction criterion.

[0010] Step 5: constructing an aboveground carbon storage inversion model and training it to obtain a trained aboveground carbon storage inversion model;

[0011] Step 6: Input the continuous-scale canopy height data and the selected characteristic variables into the trained aboveground carbon storage inversion model to estimate large-scale forest carbon storage.

[0012] Furthermore, the preprocessing in step 2 includes first aligning the spaceborne lidar data of different sensors to the same coordinate system and removing noise, and then resampling the denoised spaceborne lidar data of different sensors and rasterizing them to In the grid, is the set grid size. The overlapping spaceborne lidar data from different sensors in the grid are weighted according to their photon density and area to obtain discontinuous joint spaceborne lidar canopy height data. The weighting formula is as follows:

[0013] (2)

[0014] (3)

[0015] Where, It is i The weighted canopy height value of each grid;n For the i The total number of photon points in a grid; It is i In the grid j The canopy height value of each photon point is obtained by subtracting the reference ground elevation value from the photon point elevation value; For the i In the grid j The comprehensive weight of each photon point; is the photon spot area; is the point density, i.e. j The onboard data category to which the photon point belongs is i The number of photon points in the grid accounts for i The ratio of the number of all satellite-borne data photon points in the 9 grids centered on the grid.

[0016] Furthermore, in step 3, the discontinuous joint space-borne lidar canopy height product obtained in step 2 is used as the dependent variable, and the L-band SAR data, C-band SAR data, and terrain data are used as explanatory variables. The terrain data includes ground elevation, slope, and aspect. Multiple interaction terms between ground elevation and synthetic aperture radar data are introduced to improve the geographically weighted regression model, namely:

[0017] (4)

[0018] Where, Indicates the i The canopy height value of the discontinuous joint space-borne lidar of the grid, Indicates the i The intercept term of the grid, 、 、 、 、 Indicates the i The regression coefficient of each explanatory variable in each grid, Indicates the i Grid L-band SAR data values, Indicates the i Grid C-band SAR data values, Indicates the i The ground elevation value of each grid, Indicates the i The slope value of each grid, Indicates the i The aspect value of each grid, Indicates the iThe interaction regression coefficient between the ground elevation value of each grid and the L-band SAR data is: Indicates the i The regression coefficient of the interaction term between the ground elevation value of each grid and the C-band SAR data is: Indicates the i The third-order interaction regression coefficient between the ground elevation value of each grid and the L-band SAR data and the C-band SAR data, For the i The error term of the grid.

[0019] Use the least squares method to obtain the intercept term of the improved geographically weighted regression model , regression coefficient 、 、 、 、 and the interaction regression coefficient 、 、 , an improved geographically weighted regression model was used to estimate canopy height at a continuous scale in the study area.

[0020] Furthermore, in step 4, the C-band SAR data obtained in step 1 is used to calculate polarization parameters, and the polarization parameters include VV, VH, VH-VV / VH+VV, VV / VH, alpha, anisotropy, and entropy. The multispectral remote sensing image data obtained in step 1 is used to calculate vegetation index and texture index. The vegetation index includes normalized vegetation index, atmospheric resistance vegetation index, difference vegetation index, green normalized vegetation index, green vegetation index, height-diameter ratio index, MERIS terrestrial chlorophyll index, optimized normalized vegetation index, optimized soil-adjusted vegetation index, ratio vegetation index and red edge vegetation index. The texture index includes contrast, uniformity, roughness, directionality and entropy. The terrain data includes ground elevation, slope, and aspect. The polarization parameters, vegetation index, texture index and terrain data are used together to construct an initial feature set.

[0021] The initial feature set is used to build a random forest model based on the decision tree, and the number of decision trees is set to , the maximum depth of the tree is infinite, and the variance of each tree node s is calculated using the following formula:

[0022] (5)

[0023] Where, represents the variance of node s, Indicates that the node contains a sample set, Indicates the i The feature vector of the sample, Indicates the i The sample value of the sample, represents the regression target value, Represents the mean of the regression target value.

[0024] At node s, select the feature Q with the smallest variance after node splitting. Divide into left nodes and right node , calculate the variance corresponding to the node split based on feature Q , the formula is:

[0025] (6)

[0026] Where, For collection The number of samples in is the left node The number of samples in Right node The number of samples in Represents the left node The variance of Represents the right node The variance of .

[0027] Calculate the mean variance reduction , the formula is:

[0028] (7)

[0029] Where, K is the number of decision trees in the random forest, represents the variance before the nth decision tree is divided, Represents the variance of the nth decision tree after partitioning by feature Q.

[0030] The average variance reduction calculated according to formula (7) ,according to Sort the feature variables from large to small, remove the feature variable with the smallest average variance reduction, and put the remaining feature variables into the random forest model to repeat the above screening until the last feature variable is characteristic variables, is the set threshold.

[0031] Furthermore, in step 5, a convolutional neural network is used to construct a forest carbon storage inversion model, including a convolution layer, a pooling layer, and an activation layer. The convolution layer is designed with a convolution kernel size of F=3×3. The number of convolution kernels increases with the number of convolution layers. The activation function , the first layer convolution formula is:

[0032] (8)

[0033] Where, For the l After the convolution operation, the and the output value on channel k, For the l The weight matrix of the convolution kernel of the layer, is the first m Channels at position The pixel value at It is l The bias of the convolution kernel of the layer.

[0034] Then the pooling operation is performed to reduce the spatial size of the image while retaining important features. The pooling operation layer formula is:

[0035] (9)

[0036] Where, For the l After the layer pooling operation at position And the output value on channel k, P is the size of the pooling area, For the l -1 layer pooling area k-th channel at position The pixel value at ;

[0037] The activation layer outputs the final prediction result of the model, the formula is:

[0038] (10)

[0039] Where, y is the model output, i.e., the predicted forestland carbon storage value; is the one-dimensional vector obtained by flattening the feature maps output by the convolutional layer and the pooling layer, is the weight matrix, is the bias of the activation layer.

[0040] The airborne lidar data obtained in step 1 is segmented into individual trees to obtain tree height parameters, and then the biomass is obtained according to the allometric equation. Finally, the carbon content corresponding to the tree species is multiplied to obtain the true value of forest carbon storage, which is used as the dependent variable. The land cover category data obtained in step 1 is used to extract the part with the value of tree, and this part is used to compare the continuous canopy height data obtained in step 3 and the selected data in step 4. The characteristic variables are masked with the true carbon storage data obtained by the airborne lidar to construct the forest part dataset. The constructed forest part dataset is organized into a multi-channel matrix and normalized, and then randomly divided into training data and validation data. The forest carbon storage inversion model is trained using the training data. The training adopts the early stopping strategy. When the validation set loss is in the continuous range, the forest land carbon storage inversion model is trained. When there is no significant decrease in the round, the training is terminated and the trained forest carbon storage inversion model is obtained. The coefficient of determination, root mean square error, and bias were used to evaluate the estimated forestland carbon storage results in the study area.

[0041] The present invention also provides a large-scale forest carbon stock estimation system based on active and passive remote sensing technologies, which is used to implement the large-scale forest carbon stock estimation method based on active and passive remote sensing technologies as described above.

[0042] Furthermore, the method includes a processor and a memory, wherein the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the large-scale forest carbon storage estimation method based on active and passive remote sensing technologies as described above.

[0043] Alternatively, it includes a readable storage medium having a computer program stored thereon, and when the computer program is executed, it implements the large-scale forest carbon stock estimation method based on active and passive remote sensing technologies as described above.

[0044] Compared with the prior art, the present invention has the following advantages:

[0045] 1) A method for combining multi-platform lidar and synthetic aperture radar data was proposed, and multi-source spaceborne lidar data were weighted and merged to expand the data sample size and increase credibility; 2) The interaction term between ground elevation values ​​and synthetic aperture radar data was introduced to improve the geographically weighted regression equation, resulting in a high-precision and high-resolution large-scale continuous canopy height model; 3) A convolutional neural network-based aboveground carbon storage inversion model was constructed to improve the accuracy of forest carbon storage estimation in fragmented plots across large areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0047] Figure 1 This is a flow chart of a large-scale forest carbon storage estimation method based on active and passive remote sensing technologies according to an embodiment of the present invention.

[0048] Figure 2 This is a schematic diagram of the denoising and mapping grid for ICESat-2 ATL03 series satellite-borne lidar data according to an embodiment of the present invention. DETAILED DESCRIPTION

[0049] To make the objectives, technical solutions, and advantages of the present invention more clear, the technical solutions of the present invention are further described below with reference to the accompanying drawings and embodiments. It is obvious that the embodiments described are only some of the embodiments of the present invention, not all of them. All other embodiments obtained by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0050] Example 1

[0051] like Figure 1 As shown, the embodiment of the present invention provides a large-scale forest carbon storage estimation method based on active and passive remote sensing technology, including the following steps:

[0052] Step 1: Obtain ICESat-2 ATL03 series satellite-borne lidar data, GEDI satellite-borne lidar data, Land Exploration No. 1 satellite-borne L-band SAR data, Sentinel-1 C-band SAR data, Sentinel-2 multispectral remote sensing image data, airborne lidar data, terrain data, and land cover category data including forest land in the study area.

[0053] This example selected a city as the study area. Considering the study area's scope and woodland types, airborne LiDAR data covering 100 square kilometers was acquired, encompassing three types of woodlands: street trees, parkland, and forest. Because the accuracy of existing land cover products varies significantly across different regions, this example compared ESA and ESRI land cover data for the study area, choosing to use ESA's 10m land cover classification data.

[0054] Step 2: The ICEsat-2 ATL03 series space-borne lidar data and the GEDI space-borne lidar data obtained in step 1 are aligned, denoised, and resampled to construct a 5m×5m grid. The resampled ICEsat-2 ATL03 series and GEDI space-borne lidar data are weighted merged to obtain a discontinuous joint space-borne lidar canopy height product.

[0055] The ICEsat-2 ATL03 series satellite-borne lidar data and the GEDI satellite-borne lidar data were registered to the WGS 1984 coordinate system. The photon point data of the ICEsat-2 ATL03 series satellite-borne lidar data were denoised using the improved clustering algorithm (Density-Based Spatial Clustering of Applications with Noise, DBSCAN). The search area of ​​the DBSCAN algorithm is elliptical, and the neighborhood calculation formula is:

[0056] (1)

[0057] Where, l Center point i With the points in the neighborhood j The sum of the distances, is the total number of photon points in the neighborhood, 、 Center points i With the points in the neighborhood j The distance along the photon track, 、 Photon point i 、 j The elevation, a is the major semi-axis of the elliptical search domain, b is the semi-minor axis of the elliptical search domain.

[0058] Gaussian smoothing was used to remove noise from the GEDI satellite-borne lidar data. The canopy height value was obtained by subtracting the reference ground elevation value from the photon point elevation value of the ICEsat-2 ATL03 series and GEDI satellite-borne lidar data. Figure 2 As shown in the figure, the denoised ICEsat-2 ATL03 series spaceborne lidar data was resampled to a photon scale, and the 25m diameter GEDI spot was resampled to a 5m resolution. Both were then rasterized into a 5m×5m grid. The overlapping ICEsat-2 and GEDI data in the 5m×5m grid were weighted according to their photon density and area to obtain discontinuous joint spaceborne lidar canopy height data. The weighting formula is as follows:

[0059] (2)

[0060] (3)

[0061] Where, It is i The weighted canopy height value of each grid; n For the i The total number of photon points in a grid; It is i In the grid j The canopy height value of each photon point; For the i In the grid j The comprehensive weight of each photon point; is the photon spot area; is the point density, i.e. j The onboard data category to which the photon point belongs is i The number of photon points in the grid accounts for i The ratio of the number of all satellite-borne data photon points in the 9 grids centered on the grid.

[0062] Spaceborne lidar technology fills the gap in acquiring highly relevant data over large scales. However, existing products have varying spatiotemporal resolutions and sensor types, resulting in non-overlapping locations within the same study area. This embodiment combines the two to expand the volume of spaceborne lidar data. Furthermore, by constructing a grid for weighting, it fully accounts for data positioning deviations and enables mutual correction. Ultimately, discontinuous combined spaceborne lidar canopy height data for the study area is obtained.

[0063] Step 3: Multiple interaction terms between ground elevation values ​​and synthetic aperture radar data were introduced to improve the geographically weighted regression model. The coefficients of the improved geographically weighted regression model were obtained using the least squares method, and then the canopy height of the study area was estimated at a continuous scale.

[0064] The Land Exploration-1 satellite group 01, consisting of satellites A and B, carries an L-band multi-polarization SAR payload with a maximum spatial resolution of 3 meters. The L-band signal has a wavelength of approximately 23 centimeters, offering advantages such as long wavelength, high penetration, and immunity to cloud and fog. This allows for all-weather ground detection and more accurate surface information. Compared to other bands, the L-band can penetrate tree canopies better, enabling more precise identification of forest floor structure. Using this as an explanatory variable in canopy height estimation models can further improve the accuracy of canopy height estimation.

[0065] The discontinuous joint spaceborne lidar canopy height data obtained in step 2 is used as the dependent variable, and the L-band SAR data and C-band SAR data resampled to 5m resolution using the cubic convolution method and terrain data (ground elevation, slope, and aspect) are used as explanatory variables. Multiple interaction terms between ground elevation and synthetic aperture radar data are introduced to improve the geographically weighted regression model, namely:

[0066] (4)

[0067] Where, Indicates the i The canopy height value of the discontinuous joint space-borne lidar of the grid, Indicates thei The intercept term of the grid, 、 、 、 、 Indicates the i The regression coefficient of each explanatory variable in each grid, Indicates the i Grid L-band SAR data values, Indicates the i Grid C-band SAR data values, Indicates the i The ground elevation value of each grid, Indicates the i The slope value of each grid, Indicates the i The aspect value of each grid, Indicates the i The interaction regression coefficient between the ground elevation value of each grid and the L-band SAR data is: Indicates the i The regression coefficient of the interaction term between the ground elevation value of each grid and the C-band SAR data is: Indicates the i The third-order interaction regression coefficient between the ground elevation value of each grid and the L-band SAR data and the C-band SAR data, For the i The error term of each grid;

[0068] Use the least squares method to obtain the intercept term of the improved geographically weighted regression model , regression coefficient 、 、 、 、 and the interaction regression coefficient 、 、 , an improved geographically weighted regression model was used to estimate canopy height at a continuous scale in the study area.

[0069] In step 4, polarization parameters are extracted using the Sentinel-1 C-band SAR data obtained in step 1, and vegetation index and texture index are extracted using the Sentinel-2 multispectral remote sensing image data. Feature variables are constructed together with the terrain data, and the feature variables are trained using the random forest model. The feature variables are then screened based on the variance reduction criterion.

[0070] Polarimetric parameters are calculated using the Sentinel-1 C-band SAR data obtained in step 1. These parameters include VV, VH, (VH-VV) / (VH+VV), VV / VH, alpha, anisotropy, and entropy. VV represents vertical polarization, VH represents vertical-horizontal polarization, alpha represents the scattering angle, anisotropy represents anisotropy, and entropy represents entropy. Vegetation indices and texture indices are calculated using the Sentinel-2 data obtained in step 1. Vegetation indices include 11 types: Normalized Difference Vegetation Index (NDVI), Atmospheric Resistance Vegetation Index (ARVI), Difference Vegetation Index (DVI), Green Normalized Difference Vegetation Index (GNDVI), Green Vegetation Index (GVI), Aspect Ratio (Axis Ratio), MERIS Land Chlorophyll Index (MCVI), Optimized Normalized Difference Vegetation Index (ODVI), Optimized Soil Adjusted Vegetation Index (OSAVI), Ratio Vegetation Index (RIVI), and Red Edge Vegetation Index (REDEVI). Texture indices include contrast, uniformity, roughness, directionality, and entropy. Polarization parameters, vegetation index, texture index and terrain data (elevation, slope, aspect) are combined to construct the initial feature set.

[0071] The initial feature set is used to build a random forest model based on the decision tree. The number of decision trees is set to 500, the maximum depth of the tree is set to infinity (to ensure that all leaf nodes are pure or the minimum number of samples is reached), and the variance of each tree node s is calculated using the following formula:

[0072] (5)

[0073] Where, represents the variance of node s, Indicates that the node contains a sample set, Indicates the i The feature vector of the sample, Indicates the i The sample value of the sample, represents the regression target value, Represents the mean of the regression target value.

[0074] At this node s, select the feature Q with the smallest variance after node splitting. Divide into left nodes and right node , calculate the variance corresponding to the node split based on feature Q , the formula is:

[0075] (6)

[0076] Where, For collection The number of samples in is the left node The number of samples in Right node The number of samples in Represents the left node The variance of Represents the right node The variance of .

[0077] Calculate the mean variance reduction , the formula is:

[0078] (7)

[0079] Where, K is the number of decision trees in the random forest, represents the variance before the nth decision tree is divided, Represents the variance of the nth decision tree after partitioning by feature Q.

[0080] The average variance reduction calculated according to formula (7) ,according to Sort the feature variables from large to small, remove the feature variable with the smallest average variance reduction, and put the remaining feature variables into the random forest model to repeat the above screening until the last 9 feature variables remain.

[0081] Multicollinearity between indices can affect modeling performance. A feature screening method based on variance reduction at each tree node can assess feature importance and select informative feature variables. During decision tree construction, the greater the reduction in average variance after a tree node is partitioned by a particular feature, the stronger that feature's partitioning ability and importance are considered. This random forest feature screening method, based on average variance reduction, can identify the most influential feature variables for carbon stock regression.

[0082] Step 5: Construct an aboveground carbon storage inversion model based on a convolutional neural network. Use the true carbon storage value obtained by the airborne lidar, the continuous-scale canopy height data obtained in step 3, and the characteristic variables selected in step 4 to train the aboveground carbon storage inversion model to obtain a trained aboveground carbon storage inversion model.

[0083] A convolutional neural network is used to construct a forest carbon storage inversion model, which includes convolutional layers, pooling layers, and activation layers. The convolutional layer is designed with a convolution kernel size of F=3×3. The number of convolution kernels increases with the number of convolutional layers (K=32, 64, 128). The activation function , the first layer convolution formula is:

[0084] (8)

[0085] Where, For the lAfter the convolution operation, the and the output value on channel k, For the l The weight matrix of the convolution kernel of the layer, is the first m Channels at position The pixel value at It is l The bias of the convolution kernel of the layer.

[0086] Then the pooling operation is performed to reduce the spatial size of the image while retaining important features. The pooling operation layer formula is:

[0087] (9)

[0088] Where, For the l After the layer pooling operation at position And the output value on channel k, P is the size of the pooling area, For the l -1 layer pooling area k-th channel at position The pixel value at .

[0089] The activation layer outputs the final prediction result of the model, the formula is:

[0090] (10)

[0091] Where, y is the model output, i.e., the predicted forestland carbon storage value; is the one-dimensional vector obtained by flattening the feature maps output by the convolutional layer and the pooling layer, is the weight matrix, is the bias of the activation layer.

[0092] The airborne LiDAR data obtained in step 1 were segmented into individual trees to obtain tree height parameters. Biomass was then calculated using the allometric growth equation. Finally, forestland carbon storage was multiplied by the carbon content corresponding to the tree species and used as the dependent variable. The ESA land cover data obtained in step 1 were used to extract the tree component. This component was used to mask the continuous canopy height data obtained in step 3, the nine characteristic variables selected in step 4, and the ground-truth carbon storage data obtained by the airborne LiDAR to construct the forestland dataset.

[0093] The constructed forestland dataset was organized into a multi-channel matrix and normalized, then randomly divided into 70% training data and 30% validation data. A forestland carbon storage inversion model was trained using the training data. An early stopping strategy was employed to prevent overfitting. Training was terminated when the validation set loss did not significantly decrease after 50 consecutive rounds. The trained forestland carbon storage inversion model was obtained. The estimated forestland carbon storage results for the study area were evaluated using the coefficient of determination, root mean square error, and bias.

[0094] In step 6, the continuous-scale canopy height data obtained in step 3 and the characteristic variables selected in step 4 are input into the aboveground carbon storage inversion model trained in step 5 to estimate large-scale forest carbon storage.

[0095] The continuous-scale canopy height data and the screened characteristic variables were input into the trained aboveground carbon storage inversion model to obtain the estimated results of forest carbon storage in the study area.

[0096] Example 2

[0097] Based on the same inventive concept, the present invention also provides a large-scale forest carbon stock estimation system based on active and passive remote sensing technology, including a processor and a memory, the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute the above-mentioned large-scale forest carbon stock estimation method based on active and passive remote sensing technology.

[0098] Example 3

[0099] Based on the same inventive concept, the present invention also provides a large-scale forest carbon stock estimation system based on active and passive remote sensing technology, including a readable storage medium, on which a computer program is stored. When the computer program is executed, it implements the above-mentioned large-scale forest carbon stock estimation method based on active and passive remote sensing technology.

[0100] In specific implementation, the method proposed in the technical solution of the present invention can be automatically run by those skilled in the art using computer software technology. System devices that implement the method, such as computer-readable storage media that store the corresponding computer program of the technical solution of the present invention and computer equipment that runs the corresponding computer program, should also be within the scope of protection of the present invention.

[0101] The specific embodiments described herein are merely illustrative of the spirit of the present invention. Persons skilled in the art may make various modifications, additions, or substitutions to the described specific embodiments without departing from the spirit of the present invention or exceeding the scope of the appended claims.

Claims

1. A large-scale forestland carbon storage estimation method based on active and passive remote sensing technology, characterized in that: The following steps are involved: Step 1: Obtain spaceborne lidar data, L-band and C-band SAR data, multispectral remote sensing image data, airborne lidar data, terrain data, and land cover category data including forest land from different sensors in the study area; Step 2: preprocess and weighted merge the spaceborne lidar data from different sensors to obtain a discontinuous joint spaceborne lidar canopy height product; Step 3: Multiple interaction terms between ground elevation and synthetic aperture radar data were introduced to improve the geographically weighted regression model. The coefficients of the improved geographically weighted regression model were obtained using the least squares method, and then the canopy height of the study area was estimated at a continuous scale. The discontinuous joint spaceborne lidar canopy height product obtained in step 2 was used as the dependent variable, and the L-band and C-band SAR data and terrain data were used as explanatory variables. The terrain data included ground elevation, slope, and aspect. Multiple interaction terms between ground elevation and synthetic aperture radar data were introduced to improve the geographically weighted regression model. The intercept term, regression coefficient, and interaction term regression coefficient of the improved geographically weighted regression model were obtained using the least squares method. The improved geographically weighted regression model was used to estimate the canopy height at a continuous scale in the study area. Step 4: Use C-band SAR data to extract polarization parameters, and multispectral remote sensing image data to extract vegetation index and texture index. Together with the terrain data, these are used to construct feature variables. The feature variables are trained using a random forest model and screened based on the variance reduction criterion. Step 5: constructing an aboveground carbon storage inversion model and training it to obtain a trained aboveground carbon storage inversion model; A convolutional neural network is used to construct a forest carbon storage inversion model, including convolution layer, pooling layer and activation layer. The convolution layer is designed with a convolution kernel size of F=3×3, and the number of convolution kernels increases with the number of convolution layers. The activation function , the first layer convolution formula is: （8） Where, For the l After the convolution operation, the and channel k The output value on For the l The weight matrix of the convolution kernel of the layer, is the first m Channels at position The pixel value at It is l The bias of the convolution kernel of the layer; Then the pooling operation is performed to reduce the spatial size of the image while retaining important features. The pooling operation layer formula is: （9） Where, For the l After the layer pooling operation at position and channel k The output value on, P is the size of the pooling area, For the l -1 layer pooling area k Channels at position The pixel value at ; The activation layer outputs the final prediction result of the model, the formula is: （10） Where, y is the model output, i.e., the predicted forestland carbon storage value; is the one-dimensional vector obtained by flattening the feature maps output by the convolutional layer and the pooling layer, is the weight matrix, is the bias of the activation layer; Step 6: Input the continuous-scale canopy height data and the selected characteristic variables into the trained aboveground carbon storage inversion model to estimate large-scale forest carbon storage.

2. The large-scale forestland carbon stock estimation method based on active and passive remote sensing technology according to claim 1, characterized in that: The preprocessing in step 2 includes first registering the spaceborne lidar data of different sensors to the same coordinate system and removing noise, then resampling the denoised spaceborne lidar data of different sensors and rasterizing them to the same coordinate system. In the grid, is the set grid size.

3. The large-scale forestland carbon stock estimation method based on active and passive remote sensing technology according to claim 2 is characterized by: In step 2, weighted merging refers to The overlapping spaceborne lidar data from different sensors in the grid are weighted according to their photon density and area to obtain discontinuous joint spaceborne lidar canopy height data. The weighting formula is as follows: （2） （3） Where, It is i The weighted canopy height value of each grid; n For the i The total number of photon points in a grid; It is i In the grid j The canopy height value of each photon point is obtained by subtracting the reference ground elevation value from the photon point elevation value; For the i In the grid j The comprehensive weight of each photon point; is the photon spot area; is the point density, i.e. j The onboard data category to which the photon point belongs is i The number of photon points in the grid accounts for i The ratio of the number of all satellite-borne data photon points in the 9 grids centered on the grid.

4. The method for estimating large-scale forestland carbon stocks based on active and passive remote sensing technologies according to claim 1, wherein: The improved geographically weighted regression model expression in step 3 is: （4） Where, Indicates the i The canopy height value of the discontinuous joint space-borne lidar of the grid, Indicates the i The intercept term of the grid, 、 、 、 、 Indicates the i The regression coefficient of each explanatory variable in each grid, Indicates the i Grid L-band SAR data values, Indicates the i Grid C-band SAR data values, Indicates the i The ground elevation value of each grid, Indicates the i The slope value of each grid, Indicates the i The aspect value of each grid, Indicates the i The interaction regression coefficient between the ground elevation value of each grid and the L-band SAR data is: Indicates the i The regression coefficient of the interaction term between the ground elevation value of each grid and the C-band SAR data is: Indicates the i The third-order interaction regression coefficient between the ground elevation value of each grid and the L-band SAR data and the C-band SAR data, For the i The error term of the grid.

5. The large-scale forestland carbon storage estimation method based on active and passive remote sensing technology according to claim 1 is characterized in that: In step 4, the C-band SAR data obtained in step 1 is used to calculate polarization parameters. Polarization parameters include VV, VH, VH-VV / VH+VV, VV / VH, alpha, anisotropy, and entropy. The multispectral remote sensing image data obtained in step 1 is used to calculate vegetation index and texture index. Vegetation indices include normalized vegetation index, atmospheric resistance vegetation index, difference vegetation index, green normalized vegetation index, green vegetation index, height-diameter ratio index, MERIS terrestrial chlorophyll index, optimized normalized vegetation index, optimized soil-adjusted vegetation index, ratio vegetation index, and red edge vegetation index. Texture indices include contrast, uniformity, roughness, directionality, and entropy. Terrain data include ground elevation, slope, and aspect. The polarization parameters, vegetation index, texture index, and terrain data are used together to construct an initial feature set.

6. The method for estimating large-scale forestland carbon reserves based on active and passive remote sensing technologies according to claim 5, characterized in that: In step 4, the initial feature set is used to build a random forest model based on the decision tree, and the number of decision trees is set to , the maximum depth of the tree is infinite, and the variance of each tree node s is calculated using the following formula: （5） Where, Representation node s The variance of Indicates that the node contains a sample set, Indicates the i The feature vector of the sample, Indicates the i The sample value of the sample, represents the regression target value, Represents the mean of the regression target value; At the node s On the other hand, the feature Q with the smallest variance after node splitting is selected. Divide into left nodes and right node , calculate the variance corresponding to the node split based on feature Q , the formula is: （6） Where, For collection The number of samples in is the left node The number of samples in Right node The number of samples in Represents the left node The variance of Represents the right node variance; Calculate the mean variance reduction , the formula is: （7） Where, K is the number of decision trees in the random forest, Indicates the n The variance before the decision tree is divided, Indicates the n The variance of the decision tree after being divided by feature Q; The average variance reduction calculated according to formula (7) ,according to Sort the feature variables from large to small, remove the feature variable with the smallest average variance reduction, and put the remaining feature variables into the random forest model to repeat the above screening until the last feature variable is characteristic variables, is the set threshold.

7. The method for estimating large-scale forestland carbon storage based on active and passive remote sensing technologies according to claim 1, characterized in that: In step 5, the airborne lidar data obtained in step 1 are segmented into single trees to obtain tree height parameters, and then the biomass is obtained according to the allometric growth equation. Finally, the carbon content corresponding to the tree species is multiplied to obtain the true value of forest carbon storage, which is used as the dependent variable. The land cover category data obtained in step 1 are used to extract the part with the value of tree, and this part is used to calculate the continuous canopy height data obtained in step 3 and the selected data in step 4. The characteristic variables are masked with the true value data of carbon storage obtained by airborne lidar to construct the data set of forest land. The constructed data set of forest land is organized into a multi-channel matrix and normalized, and then randomly divided into training data and validation data. The forest land carbon storage inversion model is trained using the training data. The early stopping strategy is used for training. When the validation set loss is in the continuous range, the model is stopped. When there is no significant decrease in the round, the training is terminated and the trained forest carbon storage inversion model is obtained. is the set number of training sessions; The coefficient of determination, root mean square error and bias were used to evaluate the estimated forest carbon storage results in the study area.

8. A large-scale forest carbon storage estimation system based on active and passive remote sensing technology, characterized by: The method comprises a processor and a memory, wherein the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute a large-scale forest carbon storage estimation method based on active and passive remote sensing technology as described in any one of claims 1 to 7.

9. A large-scale forest carbon storage estimation system based on active and passive remote sensing technology, characterized by: The method comprises a readable storage medium having a computer program stored thereon, and when the computer program is executed, a large-scale forest carbon stock estimation method based on active and passive remote sensing technology as described in any one of claims 1 to 7 is implemented.

Images