Improved Forest Canopy Closure Estimation Using Multispectral Satellite Imagery within Google Earth Engine

Abstract

The large area estimation of forest canopy closure (FCC) using remotely sensed data is of high interest in monitoring forest changes and forest health, as well as in assessing forest ecological services. The accurate estimation of FCC over the regional or global scale is challenging due to the difficulty of sample acquisition and the slow processing efficiency of large amounts of remote sensing data. To address this issue, we developed a novel bounding envelope methodology based on vegetation indices (BEVIs) for determining vegetation and bare soil endmembers using the normalized differences vegetation index (NDVI), modified bare soil index (MBSI), and bare soil index (BSI) derived from Landsat 8 OLI and Sentinel-2 image within the Google Earth Engine (GEE) platform, then combined the NDVI with the dimidiate pixel model (DPM), one of the most commonly used spectral-based unmixing methods, to map the FCC distribution over an area of more than 90,000 km². The key processing was the determination of the threshold parameter in BEVIs that characterizes the spectral boundary of vegetation and soil endmembers. The results demonstrated that when the threshold equals 0.1, the extraction accuracy of vegetation and bare soil endmembers is the highest with the threshold range given as (0, 0.3), and the estimated spatial distribution of FCC using both Landsat 8 and Sentinel-2 images were consistent, that is, the area with high canopy density was mainly distributed in the western mountainous region of Chifeng city. The verification was carried out using independent field plots. The proposed approach yielded reliable results when the Landsat 8 data were used (R² = 0.6, RMSE = 0.13, and 1-rRMSE = 80%), and the accuracy was further improved using Sentinel-2 images with higher spatial resolution (R² = 0.81, RMSE = 0.09, and 1-rRMSE = 86%). The findings demonstrate that the proposed method is portable among sensors with similar spectral wavebands, and can assist in mapping FCC at a regional scale while using multispectral satellite imagery.

1. Introduction

Forests play indispensable ecological functions and act as a carbon sink, and forest canopy closure (FCC) is one of the most important indices for biophysical parameters estimation, climate response monitoring, and forest changes evaluation at the large scale [1,2,3]. Thus, accurate and timely information on forest canopy density can greatly support the decision-making process of forestry management departments. As such, it is necessary to clearly quantify the spatial distribution of FCC for forest inventories and carbon cycle estimation [4,5]. Although time consuming and costly, forest resource monitoring has a long history of being carried out at intervals of several years through conventional forestry inventory [6,7]. Remote sensing technology is cost-effective and has therefore been broadly used in assessing and monitoring forest resources because it enables the periodic monitoring of forest resources using various sensors at both global and regional scales compared to traditional field work [8,9]. Earlier studies involve the use of multisource remotely sensed data for FCC prediction, including hyperspectral, unmanned aerial vehicle (UAV) light detection and ranging (LiDAR), synthetic aperture radar (SAR), and multispectral data [10,11,12,13,14,15].

Hyperspectral data can provide narrow and contiguous spectral curves capable of capturing subtle spectral differences among different targets, which are crucial for distinguishing forests from their surrounding backgrounds [16]. To estimate FCC, Salvador et al. used the normalized difference vegetation index (NDVI) to estimate canopy closure at 10 m intervals via hyperspectral imagery, and compared the data with those estimated by ground-based methods. The results showed that the estimates of FCC using hyperspectral images compare well to those using ground-based methods [10]. The estimated accuracy of FCC using hyperspectral data was about 80% in related studies [2,10]. However, the hyperspectral imagery access is restrictive due to the high cost and low number of sensors in operation. In addition, it is time consuming to process the hyperspectral data due to the large volume, and thus requires high expertise to filter out the most effective information from a large number of highly correlated spectral bands [17,18]. Therefore, hyperspectral remotely sensed data are mainly used for the prediction of canopy density in small areas. Similarly, UAV LiDAR and airborne LiDAR data are also used in conjunction with other remote sensing imagery for the inversion of small- and medium-scale forest canopy closure, in which LiDAR data mainly provide single-point high-precision (better than 80% overall accuracy) FCC estimation to improve the overall estimation accuracy of the whole study area [2,19,20]. SAR is not affected by atmosphere and can work all day and in all kinds of weather. There are studies using the correlation analysis between the radar backscattering coefficient or its derived indices and ground samples to predict canopy density [21,22]. In addition, polarimetric classification has also been combined with SAR data for canopy closure inversion [23,24]. The high speckle noise and limited information of vegetation separability of SAR data results in an accuracy of FCC estimation of between 60% and 70%, which limits its operational application.

Compared to those data mentioned above for FCC prediction, the cost-free and easy-to-use multispectral imagery approach appears to be the best solution to FCC mapping, especially in both regional and global areas [25]. For example, mixed pixels will reduce the monitoring accuracy of remote sensing for forests, especially in heterogeneous forests [26,27]. The GF and ZY series multispectral satellite imagery with high spatial resolution has been adopted for FCC estimation because it can better capture the geometric details of the forest and reduce the influence of mixed pixels in the results [28,29,30,31]. High-resolution imagery covers a small area, so multispectral imagery at moderate spatial resolution is the primary data source for forest monitoring at large regional scales [32,33]. Landsat data at a spatial resolution of 30 m have been used in many typical studies relevant to vegetation due to its abundance of high-quality historical data [34,35]. In the cloudy and rainy areas, the Landsat time-series data can assist in capturing the forest coverage [32,36,37]. The green band of Landsat data was found to be the most effective band for the estimation of canopy characteristics in hardwood forests [38]. The spectral vegetation indices derived from the visible and near-infrared bands of multispectral images are the feature most often used to characterize the FCC distribution [39]. As another access-friendly satellite data source, Sentinel-2 multispectral data at the finest spatial resolution of 10 m have brought new opportunities for the fine monitoring of vegetation owing to its unique red-edge band and excellent spatial and temporal resolution since its launch in 2015 [40]. The great advantages of red-edge Sentinel-2 data for vegetation studies was first discovered in 2016 [41]. Since then, a large number of studies have used this data to estimate forest canopy closure [42,43,44]. The comparison of Landsat 8 and Sentinel-2 in the estimation of FCC has also been carried out [45].

Overall, the methodology used for FCC prediction at the regional scale is composed of a geometric-optical model, classification, regression statistics, and spectral unmixing [29,30,46,47]. However, the dimidiate pixel model (DPM), one of the semi-expert methods, is one of the most commonly used based on spectral vegetation indices over a large area as it reduces the effect of bias, including atmospheric disturbances and vegetative influences [48,49]. For example, the DPM was often combined with MODIS NDVI [50,51,52,53] and Landsat NDVI [54,55,56] data at the provincial and watershed scale.

In this study, we built on previous studies that used Landsat 8 and Sentinel-2 data to describe the spatial distribution of FCC within the Google Earth Engine (GEE) platform. The objective of this study was to develop a new approach based on DPM to recognize vegetation and bare soil endmembers using the multispectral images. The novelty is that the proposed method can realize regional FCC estimation that can be easily scaled for application over larger areas using GEE, and the ground-truth data are not required in the endmember extraction process.

2. Materials and Methods

2.1. Study Area

Chifeng city, centered at 119°22′58.38″E, 60°35′7.2″N, is located in the southeast of Inner Mongolia, China. Measuring approximately 90,000 km², this region consists of four northern districts (Bahrain left, Bahrain right, Linxi, and Arhorchin), four southern districts (municipal district, Aohan, Karqin, and Ningcheng), the western Hexingten, and the central Ongniud (Figure 1). Our study area belongs to the temperate semi-arid continental monsoon climate and has a variety of landforms, with mountains and hills in the north and south, and plains in the east and west. Affected by the monsoon and topography, the average annual temperature gradually rises from northwest to southeast, and the average annual temperature in most areas remains between 0 and 7 °C. Additionally, the average annual precipitation is between 350 and 450 mm, of which 60–70% is concentrated in summer, showing a decreasing trend from southwest to northeast. The dominant tree species in this area are Pinus tabulaeformis, Quercus Mongolia, Betula spp., Populus spp., Larix spp., and Armeniaca sibirica. The phenological state of the forest corresponding to the satellite data acquisition time in this area is depicted in Figure 2.

2.2. Field Survey Plots

A field survey was carried out for the main tree species in the southern four regions of Chifeng city in September 2019 to collect the forest canopy closure (FCC) data for accuracy assessment. First, considering the large scope of the study area, we set up sample plots as far as possible according to the different degrees of FCC. Then, the FCC value of each plot was measured using the line transects method. The sample points were arranged at a certain horizontal spacing, and at each sample point, the tree canopy crown was observed vertically. The FCC of the plot was obtained by dividing all of the points shaded by the canopy by the total number of sample points. Finally, the FCC data of 71 rectangular plots (30 m × 30 m) were obtained, in which the value range of the recorded FCC was (0.22, 0.97). Their spatial distribution based on the WGS84 Geographic Coordinate System is depicted in Figure 1.

2.3. Canopy Closure Estimation Overview

Our goal was to produce a high-resolution canopy closure distribution map using the powerful computing ability of the GEE platform. We proposed and implemented the novel methodology to extract the vegetation and bare soil endmembers based on the traditional DPM within the GEE cloud-computing platform, which was split into four processes, including addressing satellite images, model construction and calibration, endmembers extraction, and estimation and validation. Figure 3 provides an overview of our workflow, which is described in detail in subsequent sections.

2.4. Processing of Satellite Imagery in Google Earth Engine

The Google Earth Engine (GEE) is a platform with a petabyte analysis-ready imagery and powerful computing [57,58,59]. Three kinds of satellite imagery, composed of Landsat 8, Sentinel-2, and SRTM DEM, were used in this study. The details of their parameters are provided in

Table 1. Details of datasets used in the estimation of forest canopy closure.

Data Type	Year	Month	Spatial Resolution	Bands
Landsat 8 surface reflectance	2019	7, 8, 9	30 m	Band 4, band 5, band 6, and band 7
Sentinel-2 surface reflectance	2019	7, 8, 9	10 m	Band 2, band 4, band 8, and band 12
SRTM DEM	2000	-	30 m	-
Field survey plots	2019	9	-	-

Note: SRTM DEM represents the 30 m digital elevation model produced by the Shuttle Radar Topography Mission.

. Landsat 8 and Sentinel-2 images were used for regional forest canopy closure estimation, and SRTM DEM images were used for the preprocessing.

The open-access Landsat 8 Surface Reflectance Tier 1 (L8SR) products obtained from Landsat 8 OLI sensors were generated from specialized software for atmospheric correction using the Land Surface Reflectance Code (LaSRC) Version 3.0. The Sentinel-2 L2 surface reflectance (S2SR) data were computed by sensor2cor. In this study, the images with visible, near-infrared, and short-wave infrared bands spanning July to September in the GEE platform were used for further processing because during this period, none of the forest types (coniferous forests, broad-leaved forests) in the study area had fallen leaves. In remote sensing images, the presence of clouds will block the target objects on the ground, reducing the monitoring precision, so the pixels with cloud should be removed first. Both L8SR and S2SR data have a band identifying clouds, named “pixel_qa” and “QA60”, respectively, in which the cloud pixels were marked. Therefore, we separately used the two bands to detect and remove the pixels with cloud coverage from L8SR and S2SR, leaving only the cloudless pixels per scene to form good quality images covering the whole region.

To correct the topographic effect caused by mountainous terrain in the satellite images, the digital elevation model of 30 m produced by the Shuttle Radar Topography Mission (SRTM DEM) [60] integrated within the GEE platform was used. Additionally, the sun–canopy–sensor (SCS) with C-Correction (SCS + C), proposed by Scott et al. (2007), was implemented and applied in conjunction with SRTM DEM data to every single image scene [61]. The SCS + C correction was developed on the basis of SCS correction and overcame the overcorrection problem of SCS by introducing the C-Correction, which has been shown to be more appropriate, especially in forested areas [61].

Finally, we produced the mosaic image composite covering the whole study area by reducing the topographic-corrected cloudless L8SR image collection from July to September using the median reducer to pick median values in each band, in each pixel, over time. This is beneficial to remove pixels with high (clouds) and low (shadows) value, and to fill in the missing data after outlier removal.

2.5. Construction of Endmembers Extraction Model

According to the previous studies in which the spectral vegetation index [62,63] and surface reflectance [64] were both used for vegetation coverage estimation based on DPM, NDVI is the most commonly used vegetation index in DPM, and then the FCC calculation formula is depicted as in Equation (1):

$${\mathrm{f}}_{\mathrm{c}}=\frac{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{}{{}_{}}_{}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}}{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{v}\mathrm{e}\mathrm{g}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}}$$

(1)

where NDVI is derived from the remote sensing imagery, $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{v}\mathrm{e}\mathrm{g}}$ and $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$ $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}\mathrm{c}\mathrm{i}\mathrm{l}}$ are respectively the NDVI value of pure pixel of green vegetation and bare soil, and ${\mathrm{f}}_{\mathrm{c}}$ represents FCC [62].

As Equation (1) shows, NDVI_veg and NDVI_soil jointly determine the value of FCC. Therefore, the accurate extraction of vegetation and bare soil endmember information is critical to the FCC estimation.

In order to find the vegetation and bare soil pixels of spectral purity, we proposed a bounding envelope method based on vegetation indices (named BEVIs hereafter) based on Landsat 8 multispectral images to extract green vegetation and bare soil endmembers on the regional scale, relying on the powerful computing of the GEE platform. The BEVIs were built on the basis that the higher the vegetation index value, the higher the proportion of vegetation in the pixel, and the higher the bare soil index value, the higher the proportion of soil in the pixel.

In this study, we used the NDVI and MBSI [65] (short for the modified bare soil index proposed in the latest study in 2021) to identify forest and soil information. First, the two indices were derived from the median composite generated from L8SR according to Equations (2) and (3), and then all water pixels were masked by NDVI > 0 since the NDVI value of water bodies was generally negative to eliminate the effect of water on the statistics for pixel values as much as possible.

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}\mathrm{D}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}\mathrm{D}}$$

(2)

$$\mathrm{M}\mathrm{B}\mathrm{S}\mathrm{I}=\frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2-\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2+\mathrm{N}\mathrm{I}\mathrm{R}}+\mathrm{f}$$

(3)

where RED, NIR (i.e., near-infrared), SWIR1, and SWIR2 are band 4, band 5, band 6, and band 7 of Landsat 8 OLI, respectively, and f = 0.5 for Landsat 8 OLI [65].

Additionally, we performed statistical analysis on the NDVI and MBSI images across the entire region, and acquired their respective maximum values and standard deviations, which were denoted as NDVI_max, NDVI_std, MBSI_max, and MBSI_std. Accordingly, the regions of green vegetation and bare soil were delineated by a bounding envelope instead of a single value. The NDVI_max and MBSI_max were regarded as the upper boundary (named UB_veg and UB_soil), and the NDVI_std and MSBI_std together with an additional parameter k were used to determine the lower boundary (named LB_veg and LB_soil).

The lower boundary of vegetation regions ($\mathrm{L}{\mathrm{B}}_{\mathrm{veg}}$) and the vegetation endmembers (named EM_veg) were depicted as:

$$\mathrm{L}{\mathrm{B}}_{\mathrm{veg}}=\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\max }-\mathrm{k}\times \mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}\mathrm{t}\mathrm{d}}$$

(4)

$$\mathrm{E}{\mathrm{M}}_{\mathrm{v}\mathrm{e}\mathrm{g}}=\left[\mathrm{L}{\mathrm{B}}_{\mathrm{v}\mathrm{e}\mathrm{g}},\mathrm{U}{\mathrm{B}}_{\mathrm{v}\mathrm{e}\mathrm{g}}\right]$$

(5)

The lower boundary of soil regions ($\mathrm{L}{\mathrm{B}}_{\mathrm{soil}}$) and the bare soil endmembers (named EM_soil) were determined as:

$$\mathrm{L}{\mathrm{B}}_{\mathrm{soil}}=\mathrm{M}\mathrm{B}\mathrm{S}{\mathrm{I}}_{\max }-\mathrm{k}\times \mathrm{M}\mathrm{B}\mathrm{S}{\mathrm{I}}_{\mathrm{s}\mathrm{t}\mathrm{d}}$$

(6)

$$\mathrm{E}{\mathrm{M}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}=[\mathrm{L}{\mathrm{B}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}},\mathrm{U}{\mathrm{B}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}]$$

(7)

The key k value of the BEVIs was used to identity the information boundary between vegetation and bare soil, and the larger the k, the more raster pixels were recognized as endmembers. Therefore, parameter k is critical to the endmember determination in this study. We determine the optimal k by comparing the changes in the model accuracy when k is set to seven different values in the range (0, 0.3).

2.6. Estimation of Forest Canopy Closure and Validation

To accurately estimate the FCC, we introduced the BEVIs to delineate the vegetation and soil endmembers using Equations (5) and (7). The spectral information used to represent the two compositions in the mixed pixel was characterized by the mean value of NDVI in their respective endmembers, which can be expressed in mathematical form as Equations (8) and (9). As a result, FCC can be calculated according to Equation (1) using the extracted NDVI_veg and NDVI_soil.

$$\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{v}\mathrm{e}\mathrm{g}}=\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{E}{\mathrm{M}}_{\mathrm{v}\mathrm{e}\mathrm{g}}})$$

(8)

$$\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}=\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{E}{\mathrm{M}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}})$$

(9)

where $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{E}{\mathrm{M}}_{\mathrm{v}\mathrm{e}\mathrm{g}}}$ and $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{E}{\mathrm{M}}_{\mathrm{v}\mathrm{e}\mathrm{g}}}$ respectively represent the NDVI value of all vegetation and soil endmembers.

To quantitatively clarify the model performance of the FCC prediction, three metrics consisted of RMSE, rRMSE, and R² that were widely used to assess the results of regression operations and the 71 field survey plots that recorded the ground-truth FCC values were used for accuracy assessment.

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{({\widehat{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}})}^2}}{\mathrm{n}}}$$

(10)

$$\mathrm{r}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\frac{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}{\overline{\mathrm{y}}}$$

(11)

$${\mathrm{R}}^2=1-\frac{{\displaystyle \sum _{\mathrm{i}}{({\mathrm{y}}_{\mathrm{i}}-{\widehat{\mathrm{y}}}_{\mathrm{i}})}^2}}{{\displaystyle \sum _{\mathrm{i}}{({\widehat{\mathrm{y}}}_{\mathrm{i}}-\overline{\mathrm{y}})}^2}}$$

(12)

where ${\widehat{\mathrm{y}}}_{\mathrm{i}}$, ${\mathrm{y}}_{\mathrm{i}}$, and $\overline{\mathrm{y}}$ were respectively the ith measured FCC, the ith predicted FCC, and the average of all of the measured FCC. The parameter $\mathrm{n}$ recorded the total number of the field plots.

In addition, cross validation was used to further analyze the portability of the model by applying the calibrated BEVIs using Landsat 8 to Sentinel-2 images for FCC estimation. Because the parameter f in the calculation formula of MBSI (Equation (3)) cannot be determined in Sentinel-2 images, MBSI was replaced by BSI, which was described as Equation (13). The accuracy was validated using the same metrics and field plots as Landsat 8.

$$\mathrm{B}\mathrm{S}\mathrm{I}=\frac{(\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2+\mathrm{R}\mathrm{E}\mathrm{D})-(\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2+\mathrm{B}\mathrm{L}\mathrm{U}\mathrm{E})}{(\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2+\mathrm{R}\mathrm{E}\mathrm{D})+(\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2+\mathrm{B}\mathrm{L}\mathrm{U}\mathrm{E})}$$

(13)

where BLUE, RED, NIR (i.e., near-infrared), and SWIR2 are the band 2, band 4, band 8, and band 12 of the Sentinel-2 images.

3. Results

The model parameter k of NDVI and MBSI were first determined using Landsat 8 images, and then the k value was simultaneously applied to estimate the FCC based on both L8SR and S2SR images.

3.1. Optimal Parameter Value Determination for Vegetation Indices

The determination of the k value is the premise of the model calculation. A detailed analysis was performed based on Landsat 8 images to stabilize the model parameters. Equations (4) and (6) show the influence of the parameter k value on BEVIs, as depicted in Figure 4. The model precision gradually increases when k increases from 0 to 0.1, and the model precision gradually decreases when k is greater than 0.1.

Corresponding to different k values, the other six parameters of the BEVIs were also calculated, among which LB_veg and LB_soil changed with the k values, while UB_veg and UB_soil were fixed. When k was between 0–1, the NDVI of vegetation and soil were both stable, about 0.99 and 0.04, respectively. When k was between 0.1–0.3, the NDVI of vegetation decreased to 0.96–0.98, and the NDVI of soil increased to 0.18–0.19 (

Table 2. Details of parameter values corresponding to different k values.

k Value	UBveg	LBveg	UBsoil	LBsoil	NDVIveg	NDVIsoil
0.00	0.999	0.999	0.460	0.460	0.999	0.040
0.05		0.990		0.459	0.994	0.040
0.10		0.982		0.456	0.987	0.040
0.15		0.973		0.453	0.985	0.198
0.20		0.965		0.449	0.979	0.198
0.25		0.956		0.446	0.974	0.174
0.30		0.948		0.443	0.964	0.181

). Therefore, when k falls into the range (0, 0.1), the model parameters are most stable. Combined with the results in Figure 4, 0.1 was selected as the optimal parameter k value in this study.

3.2. Forest Canopy Closure Estimation and Validation

We carried out canopy closure inversion at a regional scale with k = 0.1 using Landsat 8 and Sentinel-2 images based on the GEE platform, and the results are respectively shown in Figure 5 and Figure 6. The spatial distribution of the canopy closure estimated by the two satellite images was extremely similar. High-value areas are mainly distributed in southern, western, and northwestern mountainous areas. The maximum value of the canopy closure derived from Landsat 8 in Chifeng city was about 0.91, and the value derived from Sentinel-2 was approximately 0.98. Because the same spectral indices and model were used, the difference in the results was probably caused by the inconsistency in image resolution and spectral range.

The scatter plots (Figure 7 and Figure 8) show the accuracy validation results using the same 71 field plots data, which proved the methods proposed in this study have great potential for application at the regional scale. When the Landsat 8 images were adopted, a reliable result was produced (R² = 0.6, RMSE = 0.13, 1-rRMSE = 80%), and the result was further improved when the higher spatial resolution Sentinel-2 was used (R² = 0.81, RMSE = 0.09, 1-rRMSE = 86%).

4. Discussion

The proposed method was able to directly determine bare soil and vegetation endmembers through statistical analysis without relying on sample points at the regional scale, which can especially alleviate the difficulty of surveying large-scale forest ground truth samples. We combined the BEVIs and the dimidiate pixel model to achieve the regional-scale forest closure estimation using Landsat 8 and Sentinel-2 satellite images. There were two key points: one was the determination of the key parameters of the model, and the other was the model robustness and accuracy verification.

4.1. Model Key Parameters Calibration

The accuracy of the vegetation and bare soil samples is critical for estimating canopy closure with the dimidiate pixel model [66,67], but samples are generally acquired by field sample plots or visual interpretation, which often results in significant time and economic costs and is susceptible to human factors. Instead, we considered using the computing power of GEE to determine endmember information through the BEVIs and satellite-acquired indices. As Figure 4 shows, a large threshold introduced more errors in endmember selection and led to a decline in model accuracy, which indicates that the estimation accuracy was not linearly positively correlated to the number of endmembers. Note that the endmember determination proposed in this study varies with satellite data sources and study areas, because the statistical results on the spectral indices generally change with specific items [67]. The NDVI and MBSI derived from the primary surface reflectance of the visible, near-infrared, and short-wave infrared bands of L8SR in Chifeng city were used to calibrate the k value and identify vegetation and bare soil, because they were proven to be the most effective indices in multispectral images that can depict the spectral information of vegetation and bare soil [68]. In the follow-up study, we will further optimize the key parameters of the model by combining more satellite spectral information such as red-edge bands.

4.2. Model Robustness and Accuracy Verification

We verified the model accuracy using independent field samples and assessed the model robustness through comparisons with different spatial resolution. The evaluation results using field plots indicate that the FCC estimation accuracy using 30 m Landsat 8 images (R² = 0.6, RMSE = 0.13, accuracy = 80%) is comparable to previous relevant studies (R² varied between 0.4 and 0.8) [42,69,70]. In the same study area, Hua et al. used three models, including multiple stepwise regression (MSR), back propagation neural network (BPNN), and a Li–Strahler geometric-optical (Li–Strahler GO) model with Sentinel-2 data for FCC estimation in the southern regions of Chifeng city. The results showed that the relative error values of the three models were 16.97%, 20.76%, and 24.83%, respectively [47]. In this study, the canopy density accuracy is further improved (R² = 0.81, RMSE = 0.09, accuracy = 86%) when 20 m Sentinel-2 images were used (Figure 7 and Figure 8). This shows that the proposed method has a reliable accuracy in the estimation of canopy density at the regional scale, and it has good robustness evidenced by comparing the application of the proposed method with two kinds of sensor images. The better performance of model accuracy in high spatial resolution images may be due to the fact that high spatial resolution images can provide more detailed information about vegetation and bare land.

4.3. Strengths and Limitations

The approach used in this study has the advantages that the modeling process does not require ground sample plots, and possesses high operating efficiency and portability utilizing the computing power of GEE. In this study, the inversion accuracy of FCC was affected by the endmember extraction precision that was determined by the k value of the BEVIs algorithm. Additionally, the estimation accuracy also depends on the acquisition time of the selected satellite image. The images in the growing season should be selected to ensure that all forest types in the study area are not defoliated and the spectral characteristics of the images are prominent. The limitation of the proposed method is that the calibration of the model parameters requires manual determination of the value range. We will further realize the adaptive calibration of model parameters in a wider range of values, continue to explore the influence of more spatial resolution images on the model accuracy, and carry out further comparative experiments in different study areas.

5. Conclusions

In the present study, we proposed a promising approach to select vegetation and soil endmembers using satellite images, and estimated the forest canopy closure at the regional scale based on DPM combined with Landsat 8 and Sentinel-2 Satellite images within the GEE platform. We found that the area with the highest canopy density in Chifeng city was mainly distributed in the western mountainous region, and the modeling method provided reliable performance, as evidenced by the accuracy based on Landsat 8 (R² = 0.6, RMSE = 0.13, and 1-rRMSE = 80%) and Sentinel-2 (R² = 0.81, RMSE = 0.09, and 1-rRMSE = 86%) images. Meanwhile, the optimal threshold value was found to be 0.1 in the range from 0 to 0.3 in intervals of 0.05. Overall, this work can help forestry authorities better understand the detailed spatial distribution of forest resources and realize the estimation of forest canopy density at the regional scale without ground plots, which indicates that it can fully utilize satellite technology to reduce the time and labor costs associated with national forest resource survey and monitoring programs.