Forest Emissions Reduction Assessment Using Optical Satellite Imagery and Space LiDAR Fusion for Carbon Stock Estimation

Abstract

Forests offer significant climate mitigation benefits, but existing emissions reduction assessment methodologies in forest-based mitigation activities are not scalable, which limits the development of carbon offset markets. In this study, we propose a measurement method using optical satellite imagery and space LiDAR data fusion to assess forest emissions reduction. Compared with the ALS-based carbon stock density estimation method, our approach presented a strong scalability for mapping 10 m-resolution carbon stock at a large scale. It was observed that dense canopy top height estimated by combining GEDI and Sentinel-2 could accurately predict forest carbon stock measurements estimated by the ALS-based method ($R^2$ = 0.72). By conducting an on-site experiment of an ongoing forest carbon project in China, we found the consistency between the emissions reduction assessed by the data fusion measurement method (589,169 ${\mathrm{tCO}}_2\mathrm{e}$) and the official ex post-monitored emissions reduction in the monitoring report (598,442 ${\mathrm{tCO}}_2\mathrm{e}$). Our results demonstrated that forest carton stock estimation using optical satellite imagery and space LiDAR data fusion is efficient and economical for forest emissions reduction assessment. The acquisition of the data was more efficient over large areas with high frequencies using space-based technology. We further discussed the challenge of building a near-real-time monitoring system for forest-based mitigation activities by utilizing optical satellite imagery and space LiDAR data and pointed out that a quality control framework should be established to help us understand the sources of uncertainty in LiDAR-based models and improve carbon stock estimation from individual trees to forest carbon projects to meet the requirements of carbon standards better.

1. Introduction

Forests, which cover 31 percent of the global land area, play a key role in the global carbon cycle. In the past two decades, global forests were a net carbon sink of 7.6 ${\mathrm{GtCO}}_2{\mathrm{eyr}}^{-1}$ [1]. Previous studies have shown that letting forests regrow naturally has the potential to absorb up to 8.9 billion metric tons of carbon dioxide equivalent from the atmosphere each year through 2050, which accounted for nearly one-quarter of global fossil fuel emissions from the atmosphere every year [2]. Therefore, forests offer significant climate mitigation benefits, and monitoring and managing forests for carbon sequestration offer the opportunity to reduce greenhouse gas emissions and stabilize the Earth’s temperature. Accurate information about forest cover, forest type, and carbon stock are critical ecological indicators for actions aiming at reducing emissions from deforestation and forest degradation (REDD+).

The carbon market is an important tool to reach global climate goals in the short and medium term, in which market regulation incentivizes emissions reductions, carbon offsets, and investment in technologies to reduce emissions. As forestry is the current dominant carbon offset solution, forest carbon credits constitute a large part of the total offset–credit transaction volume in the global or local voluntary carbon markets [3]. The development of forest carbon projects in the voluntary carbon market is governed by a number of major carbon standards: Verified Carbon Standard, Plan Vivo, Gold Standard, Climate Action Reserve, China Certified Emission Reduction, Fujian Forestry Certified Emission Reduction, etc. Following these standards, forest carbon projects, after achieving validation, would periodically conduct monitoring, reporting, and verification (MRV) to assess the consistency between estimated emissions reduction before the project implementation and monitored emissions reduction after the project implementation.

Existing emissions reduction assessment methodologies in carbon standards are costly and laborious, which require forest inventories to obtain sampling biophysical parameters, e.g., diameter at breast height (DBH) and tree height (TH) [4]. These biophysical parameters are then used to estimate biomass from species-specific allometric equations. Carbon stock for assessing emissions reduction are finally calculated from the estimated forest biomass. Meanwhile, only a sample of the target forest is monitored in existing emissions reduction assessment methodologies, which makes forest carbon projects completely reliant on forestry experts’ opinion to select sample plots. Therefore, existing emissions reduction assessment methodologies are not robust enough, and there is room for forest carbon projects to overstate their emissions reductions. To address these problems, [4] presented an airborne-laser-scanning (ALS)-based measurement method to assess forest emissions reduction. They demonstrated the capability of ALS technology to quantify forest biomass and associated emissions reduction for forest-based ${\mathrm{CO}}_2$ mitigation activities. However, this emissions reduction assessment method is not scalable, so it is still costly for most forest carbon projects, because forest carbon project areas usually consist of multiple discrete project area parcels scattered over large areas, which are filtered by complex implementation applicability conditions defined in carbon standards. Hence, there is an urgent need for a highly automated, more economical approach to monitor forest carbon stock density and quantify forest emissions reduction for REDD+ projects.

LiDAR applications for forest aboveground biomass and associated carbon stock quantification have been developed over a range of spatial resolutions, forest types, and ecosystems [5,6,7]. In practice, airborne discrete-return small-footprint LiDAR systems are the most widely used for biomass estimation [8], providing detailed forest structure variables of either individual trees or areas. Although this technology has shown its strong capability of assessing forest emissions reduction, it can only be adopted for small areas because of data availability constraints. Recently, the development of space LiDAR has provided a new solution for high-frequency large-scale carbon stock mapping. Key datasets containing millions of LiDAR footprints are generated by the new generation of space LiDAR sensors, such as the Global Ecosystem Dynamics Investigation LiDAR (GEDI) and ICESat-2, which improve the ability to characterize important carbon cycling processes [9]. Moreover, the new novel multi-sensor Earth observation methods combining spaceborne LiDAR data, SAR data, and optical data using machine learning techniques brings precise measurements of forest canopy height and carbon stock monitoring towards large-scale, more efficient, and actionable forest emissions reduction assessment [10,11,12,13,14,15,16,17,18,19,20]. For example, ref. [11] produced the first wall-to-wall map of aboveground live woody vegetation carbon density at 500 m resolution for the pan-tropics using a combination of MODIS and GLAS remote sensing data. Reference [13] developed an aboveground biomass mapping method using a combination of ground observations, MODIS data, forest cover/gain/loss maps based on Landsat, GLAS forest canopy height, and climatic and terrain data. Reference [19] used GEDI data fused with Sentinel-1 SAR, Sentinel-2 multispectral, elevation, and land cover data to produce biomass maps of Australia and the United States. However, to the best of our knowledge, optical satellite imagery and space LiDAR data fusion have not been tested in REDD+ projects to quantify forest emissions reduction.

This study presents a practical framework to assess forest emissions reduction using optical satellite imagery and space LiDAR data fusion, which provides a method using multiple remote sensing data sources to model carbon stock from individual trees to REDD+ projects. By conducting an on-site experiment of an ongoing forest carbon project, we evaluated the method’s performance compared with the official ex post-monitored emissions reduction in the monitoring report. Our results demonstrated that forest emissions reduction assessment using optical satellite imagery and space LiDAR data fusion is efficient and economical. We also discuss the challenge of building a near-real-time monitoring system for forest-based mitigation activities by utilizing optical satellite imagery and space LiDAR data.

2. Data and Methods

Detailed below is our approach to assess forest emissions reduction using multiple remote sensing data sources (see Figure 1). In this work, we used four major data sources: field measurements, ALS data, Sentinel-2 optical images, and sparse canopy top height collected from GEDI L2A. We developed ALS-based carbon stock estimation models to assess forest emissions reduction from individual trees to individual sampled plots. Using optical satellite imagery and space LiDAR data fusion by a deep fully convolutional neural network (FCN), we estimated the canopy top height for each 10 m pixel within the project boundary. Combined with plot-level carbon stock density values, the project-level canopy top height product was translated into carbon density by learning the ordinary least-squares (OLS) regression models. Finally, the total project carbon stock and the corresponding forest emissions reduction can be easily assessed with the help of the project-level carbon stock density product.

2.1. Research Area

Our analysis focused on the Fujian Carbon Neutrality Forest Project in Shunchang Country, Fujian, China (${117}^{°}{29}^{\prime }$–${118}^{°}{14}^{\prime }$E, ${26}^{°}{38}^{\prime }$–${27}^{°}{12}^{\prime }$N), which is developed under the Fujian Forestry Certified Emission Reduction standard. Shunchang Country’s climate is characterized as humid subtropical, and the main tree species in the project area is Cunninghamia lanceolata. The project covers a total area of 3395.51 ha, and 909 discrete project area parcels scattered over mountainous Shunchang Country make up the whole project, which makes an ALS-based measurement method hardly suitable for assessing forest emissions reduction in this case, because of the high implementation costs. When the project was first monitored in 2022, the estimated total carbon stock from field measurements was 1,121,855 ${\mathrm{tCO}}_2\mathrm{e}$. In 2021, the Paulson prize for sustainability was awarded to the local authority of Shunchang Country for the development of the Forest Eco Bank and Carbon Sink Project.

2.2. Field Data and Biomass Allometry

Field measurement data collection was conducted by the Shunchang state-owned forestry farm on 1 July 2022. Eight field sample plots were randomly selected within eight sampled LiDAR-covered areas (see Figure 2).

In each sample plot, we selected a 20 m × 20 m square. Information on the diameter at breast height (DBH), tree height (TH), and GPS coordinates for each tree in the square was collected. A total number of 290 trees were selected to develop our ALS-based individual tree biomass estimation models. Note that individual tree biomass includes aboveground and belowground biomass. For individual trees, the biomass ${\mathrm{b}}_{\mathrm{i}}$ of species i is estimated by using an allometric equation that combines key tree features [21], including TH, DBH, wood density (WD), the biomass expansion factor (BEF), and the ratio of belowground and aboveground biomass (R) [22,23]. These allometric equations take the form:

$$b_i=f(TH,DBH,W{D}_i,BE{F}_i,R_i),$$

(1)

in which species-specific tree features WD, BEF, and R are provided by the carbon standards. The allometric equation we used is defined by the Fujian Forestry Administration. For Cunninghamia lanceolata (CL), we have

$$b_{CL}=0.00005806\times DB{H}^{1.955335}\times H^{0.894033}\times 307\times 1.634\times (1+0.246),$$

(2)

in which $W{D}_{CL}$ is 307 $\mathrm{kg}·{\mathrm{m}}^{-3}$, $BE{F}_{CL}$ is 1.634, and $R_{CL}$ is 0.246.

After the total biomass was obtained by summing individual biomass values for all trees within the plot, the forest carbon at the plot level was then translated from total biomass using the carbon fraction of Cunninghamia lanceolata (0.52).

2.3. ALS Data

We collected ALS data using an unmanned aerial vehicle (UAV) equipped with a long-range laser scanner (Zenmuse L1) LiDAR system throughout 8 sampled plots in Shunchang in July 2022. These eight areas encompass 30 ha in total. With a nominal flying height of 100 m and a pulse frequency of 240 kHz, the resulting point cloud featured an average pulse density of 482.9 pulses/${\mathrm{m}}^2$. The horizontal and vertical accuracies of the point cloud were, respectively, 10 cm and 5 cm.

2.4. Optical Satellite Imagery and Space LiDAR Data

Twelve bands (including B02, B03, B04, B05, B06, B07, B08A, B09, B01, B10, B11, B12) of the atmospherically corrected Sentinel-2 L2A product were used in our study. We collected the Sentinel-2 L2A images (4 Sentinel-2 tiles) of our study area with the lowest overall cloud coverage in July 2022. Ten-meter-resolution bands for all the existing spectral bands with 20 m and 60 m were created using bilinear upsampling [16]. After upsampling, these twelve bands were merged into a single tiff file. The Sentinel-2 L2A scene classification results were used to get rid of non-vegetated areas, so we could guarantee that the model merely estimates canopy top height for vegetated areas. Note that we did not generate any vegetation indices from Sentinel-2 imagery for canopy top height mapping. Vegetation indices, such as the Normalized Difference Vegetation Index (NDVI), are merely used by the Sentinel-2 Level-2A algorithms to generate a scene classification map, which includes six different classifications for shadows, cloud shadows, vegetation, soils/deserts, water, and snow.

For space LiDAR, all GEDI pulses acquired over Shunchang were collected in March 2022, the closest GEDI time available to the field inventory. Shots flagged as high quality by the GEDI L2A product were filtered to ensure the sampling data were suitable for biomass estimation. To control for GEDI data quality, we filtered all data using a beam sensitivity threshold greater than 0.98 [24]. A total of 9500 footprint-level canopy top heights were then derived from the filtered GEDI waveforms. The GEDI L2A product provides relative height (RH) metrics for each footprint, which can be used as the candidate predictor variables for the biomass estimation models. In particular, RH98, giving the height at which 98% of returned energy is reached relative to the ground, was used as our canopy top height as it is a more stable height metric than RH100 [25]. After rastering the RH98 data to the Sentinel-2 pixel grid, Sentinel-2 and RH98 image pairs of 150 m × 150 m pixels centered at each waveform location were clipped.

2.5. LiDAR Data Pre-Processing and Individual Tree Segmentation

The advent of ALS technology makes directly estimating plot-level biomass by summing the biomass of individual trees efficient, and it also has the potential to substantially reduce the uncertainties associated with tree height variations [26]. Similar to field measurements, the ALS-based LiDAR biomass estimation approach is also tree-centric. In the LiDAR data pre-processing pipeline, we firstly removed noise in the original point clouds using the statistical outliers removal method, then the processed point clouds were classified into ground and non-ground points using the progressive morphological filter algorithm [27]. A high-resolution (0.5 m) digital terrain model and digital surface model were then extracted and used to create a canopy height model (CHM) of the area.

After the LiDAR point clouds were normalized, individual tree segmentation can be conducted to obtain the point clouds of each tree. We used the region-based segmentation approach developed by [28,29] to delineate individual trees, which was implemented by finding the local maximums in a rasterized CHM. The complete approach consisted basically of three main steps: (1) the CHM was pre-filtered by a low-pass filter, which can smooth the surface and reduce the number of local maximums; (2) the local maximums of the CHM were then identified using a circular moving window; (3) the extracted local maximums were labeled as the initial seed points, then the seeded region was grown by adding all the neighbors to the region if their vertical distance from the local maximum was less than the difference threshold until no further pixels were added to the region. Note that the number of individual trees in the area covered by the ALS data was determined by the parameters of the segmentation algorithm. To achieve high segmentation accuracy (95.7%), field measurements should be used to tune these parameters.

2.6. Estimating Plot-Level carbon from ALS Data

The LiDAR derivative metrics were extracted from the tree point clouds after individual tree segmentation. In this study, we considered the following derivative metrics, including TH, crown diameter, crown volume, mean height, and simulated relative height (RH) metrics at the 98th, 90th, 80th, 70th, 60th, 50th, 40th, 30th, 20th, and 10th height above ground-level percentiles (RH98, RH90, RH80, RH70, Rh60, RH50, RH40, Rh30, RH20, and RH10, respectively).

To estimate the tree-centric carbon stock c from a specific LiDAR metric x, we used the OLS regression model with possible transforms, which can be written as

$$g\left(c\right)=\alpha *f\left(x\right)+\beta ,$$

(3)

where $\alpha$ is the regression coefficient and $\beta$ is the constant term. $f(·)$ is a transformation function (identity, square root), and $g(·)$ is a back-transformation function (identity, second power). We can conduct linear regression analysis or log-transformed linear regression analysis on ground-based biomass against LiDAR-derived feature parameters.

Our models were built on a total number of 290 Cunninghamia lanceolata trees. We split the dataset into two subsets, holding out a random set of 87 trees (30 %) for validation. Once the OLS models have been trained, individual tree carbon stock estimation will be easily extrapolated for use in the rest of the project area.

In order to assess the results, the following metrics were applied in the validation dataset: (1) The coefficient of determination ($R^2$) was used for the linear models. Assuming we have an observed dataset of n values marked $y_1$,...,$y_n$, each associated with a modeled value $\widehat{y_1}$,...,$\widehat{y_n}$, the definition of $R^2$ is:

$$R^2=1-\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{{\sum }_{i=1}^n{(y_i-\overline{y_i})}^2},$$

(4)

in which $\overline{y}$ is the mean of the observed data. (2) The $pseudo-R^2$ was used for the machine learning and nonlinear models:

$$pseudo-R^2=1-\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{{\sum }_{i=1}^n{(y_i-\overline{y_i})}^2}.$$

(5)

(3) The root-mean-squared error ($RMSE$) and the relative-root-mean-squared error ($rMSE,\%$):

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{n}},$$

(6)

$$rRMSE=100\times \frac{RMSE}{\overline{y}}.$$

(7)

(4) The mean bias error ($ME$) and the relative mean bias error ($rME,\%$):

$$ME=\frac{{\sum }_{i=1}^n(y_i-\widehat{y_i})}{n},$$

(8)

$$rME=100\times \frac{ME}{\overline{y}}.$$

(9)

Then, we used the OLS model to compute the carbon stock for each tree for 199 ALS-covered areas. After obtaining the sum of the carbon stock from each tree within each 10 m pixel aligned with 200 Sentinel-2 L2 data, we further generated a map of the carbon stock density for each plot.

2.7. Sentinel-2 GEDI Data Fusion for Dense Canopy Top Height Mapping

Because the monitoring cost is still a major obstacle to more widespread and general adoption of ALS, it is impractical to use ALS data merely to create a wall-to-wall carbon stock density map for the whole project area at a higher frequency than traditional field measurement sampling methods. In this study, we attempted to employ an FCN to estimate the canopy top height for every 10 m × 10 m pixel of the Sentinel-2 image and then produced a carbon stock density map.

The FCN is a classic pixel-to-pixel learning system that combines semantic information from a deep, coarse layer with appearance information from a shallow, fine layer to produce segmentations or translate cross-modality information from image to image [30]. Our FCN (see Figure 3), which was trained on the Sentinel-2 L2A image and the rasterized RH98 data pairs, contains 4 residual blocks [31] and ends with a convolutional layer with 1 × 1 filter kernels. The convolutional layers mostly have 3 × 3 filters and 256 output channels. This makes the model learn to extract relevant vegetation indexes from Sentinel-2 images. We did not perform downsampling to make the output of each layer have the same height and width as the input. Once trained, the model can process arbitrarily sized input image patches.

The canopy top height reference data derived from GEDI were sparsely distributed, but the estimated canopy top height by the FCN was dense. Therefore, the FCN was trained with sparse supervision. At training time, only the pixels with a valid GEDI reference height were used to calculate the mean-squared error (MSE) loss; therefore, only at valid reference pixels, the model parameters were optimized. We trained on one GPU for 100k iterations, using the ADAM optimizer [32] with a learning rate of 0.001. We split the dataset into two subsets, holding out a random set of 30% of the data for validation.

2.8. Estimating Carbon Stock Density from Dense Canopy Top Height and ALS Data

A linear regression analysis and a log-transformed linear regression analysis of the ALS-based carbon stock density for every 10 m × 10 m pixel against the canopy top height were conducted to establish the relationship between the carbon stock density and canopy top height for the forest carbon project. A project carbon stock density map was then created.

Note that our method applies a filter between space and airborne LiDAR estimated maximum height to guarantee robust model performance. If there was more than a 5 m difference between space and airborne LiDAR estimated maximum height, the corresponding carbon stock samples estimated from airborne LiDAR data were excluded. We split the dataset (around 3000 pixels) into two subsets for regression analysis, holding out a random set of 30% of the data for validation.

2.9. Project-Level Carbon Stock Assessment

In the traditional field measurement method, it is hardly possible to assess the carbon stock within a sufficiently large area. Therefore, under the guidance of the carbon standards, a project area was partitioned into forest strata that were homogeneous in their carbon stock density. For each forest stratum s, the average carbon stock density $C{D}_s$(Mg/ha) was calculated based on the average carbon stock density of representative sample plots. Each stratum s is usually composed of multiple sample plots J.

$$C{D}_s=\frac{1}{J}{\sum }_1^{J}C{D}_{s,j}.$$

(10)

The total carbon stock within the project boundary is the sum of the total carbon stock of all forest strata S.

$$C={\sum }_1^{S}C{D}_s*A_s,$$

(11)

in which $A_s$ is the area of stratum s.

In our method, we can assess the total project carbon stock by summing all the carbon stock of 10 m × 10 m pixels within the project boundary. The carbon dioxide equivalent for forest emissions reduction is then translated from total project carbon stock.

3. Results

3.1. Performance of Tree-Centric ALS-Based Carbon Stock Estimation

We present summaries of the performances of different ALS feature parameters in estimating individual tree carbon stock. The best-fit ALS metric and model were identified based on the $RMSE$ of the carbon stock (in kilogram). On the validation dataset, forcing TH into models yielded the highest sensitivity to the carbon stock, while RH98 was the second-strongest correlated metric.

The linear regression provided a reliable goodness-of-fit. More than 60% of the variation in the field-measured carbon stock was explained by the ALS-derived tree height, and the $RMSE$ was 23.1241 kg ($R^2$ = 0.6014, Figure 4).

The OLS regression model with the log transformations of the predictor and response variables had a better goodness-of-fit. Figure 5 shows a closer relationship between the ALS-derived tree height and the field-measured carbon stock ($pseudo-R^2$ = 0.6228, $RMSE$ = 20.3396 kg). This model also yielded a lower $rMSE$, $ME$, and $rME$. This phenomenon was caused by the power-law relationship between the carbon stock and canopy top height [33].

Compared with the field measurement method, the ALS-based individual tree carbon stock estimation model had a low relative mean bias error ($rME=3.22\%$). The measurement costs of the ALS-based model were only about one-ninth the costs of the field measurement method. Because of such data acquisition efficiency, ALS-based carbon stock estimation can be applied in plots with more extensive area coverage.

3.2. Performance of Dense Canopy Top Height Mapping

We evaluated the performance of the dense canopy height mapping on the held-out GEDI reference data. Over all unseen GEDI reference samples, although the low canopy top height was overestimated (see Figure 6), our trained FCN yielded a high $R^2$ value of 0.8142 and an $RMSE$ of 3.5467 m. The project-level canopy height product derived from the Sentinel-2 L2A images by the FCN is shown in Figure 7.

We further compared the height predictions from the FCN model with the random forest (RF) model in [34,35], the support vector regression (SVR) in [36,37] and the convolutional neural network (CNN) model in [38]. In

Table 1. Comparison of the error metrics for the different canopy height models: FCN, RF, SVR, and CNN.

Model	pseudo-R2	RMSE (m)	rRMSE (%)	ME (m)	rME (%)
FCN	0.8142	3.5467	16.89	−1.2329	−5.87
RF	0.6279	6.1348	29.21	−0.3176	−1.51
SVR	0.7037	6.3252	30.12	1.5633	7.44
CNN	0.7856	4.6730	22.25	2.2397	10.67

, we can see that the FCN model had the highest $R^2$ and the lowest $RMSE$.

3.3. Performance of Project Emissions Reduction Estimation

ALS-based carbon stock density maps for individual plots were rasterized at a 10 m resolution, which is the same as the resolution of the dense canopy top height mapping. The statistical results of the carbon stock density for each plot are shown in

Table 2. Mean, standard deviation (SD), minimum (Min), and maximum (Max) values of individual plot carbon stock density maps produced using the ALS-based methods.

Plot	$\mathbf{Mean} (\mathbf{Mg}/\mathbf{ha})$	SD (Mg/ha)	Max (Mg/ha)	Min (Mg/ha)
1	83.16	8.47	121.35	56.31
2	93.98	6.51	117.72	70.02
3	23.62	7.18	50.82	7.49
4	49.86	6.16	68.23	20.58
5	92.45	9.82	131.79	76.63
6	73.76	11.96	109.90	57.79
7	137.42	25.15	185.90	95.88
8	118.66	18.92	178.46	92.59

, measured in units of megagrams (Mg) of carbon per hectare (ha). The highest levels of carbon stock density occurred in Plot 7, while the lowest values occurred in Plot 3. The field data in the official monitoring report showed the same tendency.

By conducting a linear and a log-transformed linear regression analysis on the carbon stock density against the canopy top height, we demonstrated that the carbon stock density estimates using optical satellite imagery and space LiDAR fusion had $R^2$ values over 0.72 relative to the ALS-based carbon stock density, over all validation samples (see Figure 8a,b). We observed that the model with the log transformations of the predictor and response variables also exhibited a lower $RMSE$, $rRMSE$, $ME$, and $rME$ and a higher $R^2$ value, indicating that the regression model with the log transformations performed better and the canopy top height derived using the optical satellite imagery and space LiDAR fusion was a key indicator for project-level carbon stock density (biomass) mapping, which are consistent with the results from previous studies [39].

It is clear from

Table 3. Mean, standard deviation (SD), minimum (Min), and maximum (Max) values of individual plot carbon stock density maps produced using the data fusion methods.

Plot	$\mathbf{Mean} (\mathbf{Mg}/\mathbf{ha})$	SD (Mg/ha)	Max (Mg/ha)	Min (Mg/ha)
1	85.47	9.74	105.46	55.63
2	91.70	10.98	105.18	67.52
3	29.12	7.85	47.53	9.40
4	56.92	10.59	71.75	12.92
5	98.17	12.57	126.40	77.71
6	77.90	9.70	101.78	57.67
7	129.40	22.24	173.13	99.34
8	107.63	15.34	161.74	86.17

that the plot-level statistic results produced by the log-transformed linear regression model had an acceptable level of consistency with the ALS-based model. We can observe that the data fusion method slightly overestimated the low carbon stock density and underestimated the carbon stock density.

Using the dense project-level carbon stock density map, we calculated the emissions reduction assessed by our data fusion measurement method (589,169 ${\mathrm{tCO}}_2\mathrm{e}$), which was highly consistent with the values calculated from field measurements in the official monitoring report (598,442 ${\mathrm{tCO}}_2\mathrm{e}$).

4. Discussion

4.1. Benefits of Optical Satellite Imagery and Space LiDAR Data Fusion

Compared with the ALS-based assessment method proposed by [4], using optical satellite imagery and space LiDAR fusion does not add any cost, because Sentinel-2 and GEDI data are free for everyone. Moreover, without the use of fixed-wing aircraft, it is difficult to apply the ALS-based method over large areas of complicated landforms, such as mountainous regions, which are considered to provide the special mixture for optimum forest carbon storage [40]. This is because the flight time of a UAV LiDAR system is limited due to the battery capacity. The Global Navigation Satellite System real-time kinematic (GNSS RTK) receiver also influences the uncertainty of laser measurements markedly if the flight duration is long [41]. The typical industry rate for aerial mapping using UAV-mounted LiDAR is currently around USD 2500 per day. The rate for a fixed-wing aircraft LiDAR system is around USD 30,000 per day. The cost for using a fixed-wing aircraft LIDAR system is higher than the labor cost for the measurement and monitoring using traditional field measurements. Current machine learning techniques provide the opportunity to utilize optical satellite imagery and space LiDAR to extrapolate regional measurements. Therefore, the use of optical satellite imagery and space LiDAR fusion goes a step further to reduce the cost of assessing forest emissions reduction. It makes mapping near-real-time carbon stock for large areas possible.

On the other hand, the method of optical satellite imagery and space LiDAR fusion takes full advantage of the data collected by the ALS from a wider sample coverage. In our study case, the airborne LiDAR covered eight plots, with a total area of 30 ha, while the field measurements only covered a total area of 3 ha. This may reduce the statistical error of sampling and improve the performance of carbon stock assessment.

4.2. Comparison with Other Canopy Height and Biomass Carbon Stock Estimation Models

For most evaluation metrics, the FCN model showed the best canopy height mapping performance for the subtropical Cunninghamia lanceolate forest in China. The other deep-learning-based model (CNN) also showed its advantage for canopy height mapping, which suggests the high potential of deep learning in spatially continuous canopy height estimation. However, some problems still remain when deep learning models are applied to map canopy height. For example, low canopies are overestimated and tall canopies are underestimated by the FCN [42]. More precise deep-learning-based canopy height models should be explored to reduce this effect.

Compared with other biomass carbon stock estimation models [19,43,44], the methodology developed in this study did not combine GEDI data with Sentinel-1 SAR, elevation, land cover data, or ICESat-2 data, because it was demonstrated that using only Sentinel-2 bands for canopy height estimation can lead to relatively good error metrics [16,42]. For REDD+ projects in highly cloudy regions, the proposed methodology should consider Sentinel-1 data, when cloud-free Sentinel-2 data are not available.

4.3. Major Challenges from Field Data Collection to Satellite-Based Project-Level Carbon Stock Mapping

In the stage of field plot biomass estimation from airborne LiDAR data, CHM-based individual tree segmentation methods show reduced performance for multiple layered, mixed species stands of mixed forests [45,46]. Not all of the trees can be recognized by the algorithm. This is because model parameters such as the size of the moving window, the small tree cut-off height, and the percentage and absolute height difference thresholds are difficult to determine for multiple tree species with heterogeneous properties [47]. In our case, even though the dominant tree was Cunninghamia lanceolata, the model parameters varied among plots. To overcome the issues associated with the mismatch between field measurements and individual tree segmentation results, we suggest that the parameters adopted in segmentation models should be carefully tuned, using the geographical and planting information presented in forest inventory reports as a reference. Moreover, the fusion of LiDAR and spectral information, providing a better tree extraction [29,48], should be considered in future work.

When we estimate the total project carbon stock from dense canopy top height maps, a reliable practice is to apply a filter between airborne and space LiDAR estimated maximum height, which yields a more substantial increase in model fitting performance because of the better spatial match between both canopy top height estimates [39]. In our case, for each rasterized pixel, if there was more than a 5 m difference between the space and airborne LiDAR estimated maximum height, the corresponding carbon stock sample was excluded. This filter significantly removed the noise caused by dense canopy top height mapping and precluded GEDI’s 10 m geolocation uncertainty for footprint locations. In future work, we should also consider the uncertainty of deep learning models used to retrieve canopy height from Sentinel-2 images. By integrating Bayesian deep learning techniques [49], the per-pixel probability of uncertainty can be directly used to filter out inaccurate canopy height predictions. These techniques will help us understand the sources of error and pathways to improve carbon stock estimation using optical satellite imagery and space LiDAR fusion.

As the relationships between wall-to-wall remote sensing data and carbon stock are complex and context-dependent [50], most studies have relied on model fitting to identify the most-predictive metrics in a specific given context [49]. To reuse pre-trained carbon stock estimation models, open standard (ground-level) datasets that could be used to cross-validate the carbon stock estimation results from multiple sensors deployed on different platforms with different acquisition parameters should be created in the future. These datasets stratified by plant functional type and geographic region must have undergone rigorous, consistent, and transparent validation. Without enough reference data, it is unlikely that a carbon estimation model can be transferable across forest types and regions to quantify forest emissions reduction for REDD+ projects.

Uncertainty analysis is a major challenge for satellite-based project-level carbon stock mapping. When carbon stock estimation models are first established for individual trees and then applied to predict aggregate plot-level or regional carbon stock by developing a machine learning model with remote sensing data, it is not clear how uncertainties should be estimated [51]. Therefore, using remote sensing data to assess forest emissions reduction hastily may result in sub-optimal or even detrimental policy measures for forest management [52].

In summary, optical satellite imagery and space LiDAR fusion has shown the potential to be an efficient monitoring approach added in the methodology for forest emissions reduction development, because of its lower cost, satisfactory accuracy, and scalability. Building a robust monitoring, reporting, and verification (MRV) system of carbon emissions using these Earth observation techniques allows forest emissions reduction assessment to have a spatial and temporal resolution never seen before. High-quality field data are still needed for both calibration and validation at the current stage. Further work should particularly look into how a quality control framework including improved uncertainty estimates can be established to guide upscale carbon stock estimates from individual trees to forest carbon projects, before Earth observation techniques are integrated into the methodology for forest emissions reduction development.

5. Conclusions

In this study, we proposed a measurement framework using optical satellite imagery and space LiDAR data fusion to assess forest emissions reduction. Compared with the ALS-based carbon stock density estimation method, our approach can be applied for mapping carbon stocks at a 10 m resolution across large areas typically associated with REDD+ projects. By conducting an on-site experiment of an ongoing forest carbon project following the Fujian Forestry Certified Emission Reduction standard, we found an acceptable level of consistency between the emissions reduction assessed by the LiDAR-based measurement method (589,169 ${\mathrm{tCO}}_2\mathrm{e}$) and the official ex post-monitored emissions reduction in the monitoring report (598,442 ${\mathrm{tCO}}_2\mathrm{e}$). Our results demonstrated that forest carton stock estimation using optical satellite imagery and space LiDAR data fusion is efficient for forest emissions reduction assessment because of its strong data acquisition efficiency in terms of monitoring frequency and area coverage. We further discussed the challenge of building a near-real-time monitoring system for forest-based mitigation activities by utilizing optical satellite imagery and space LiDAR data and pointed out that a quality control framework should be established to help us understand the sources of uncertainty in LiDAR-based models and improve carbon stock estimation from individual trees to forest carbon projects to meet the requirements of carbon standards better.