High-Resolution Mapping of Shrub Biomass in Arid and Semi-Arid Environments by Integrating Remote Sensing Observations across Multiple Spatial Scales

1. Introduction

The need to protect the ecological environment and restore ecosystems has grown in recent years due to the effects of human activity and global climate change [1,2,3,4]. As an important vegetation type, shrubs play a crucial role in ecosystems [5,6,7]. They not only maintain soil moisture balance and prevent soil erosion but also provide habitats and fodder [8,9,10], which are essential for preserving biodiversity and maintaining ecological balance [11]. Among their attributes, shrub biomass is a crucial measurement indicator for reflecting shrub growth status and overall vegetation productivity [12,13]. Therefore, accurately estimating shrub biomass is a key task for understanding vegetation growth status, ecological functions, and ecosystem carbon storage [14,15,16].

Conventional methods of field surveying have limitations in accurately estimating the biomass of shrubs [17]. Firstly, these methods necessitate a substantial allocation of human resources and financial investment and are time-consuming, particularly in expansive areas with intricate topography [18,19]. Secondly, field surveys can only yield localized sample data, making it challenging to capture the spatial distribution and temporal patterns of shrubs fully. Conversely, the use of destructive sampling techniques to quantify shrub biomass may potentially impact the surrounding ecological system [20]. Furthermore, the biomass of shrubs is influenced by various factors such as vegetation structure, growth environment, and climatic conditions [21,22]. Consequently, it is challenging to account for these factors using traditional methods comprehensively.

To address the limitations of conventional methods, remote sensing technology has emerged as an effective tool for estimating shrub biomass [23,24,25]. This technology enables the acquisition of extensive and continuous data with high spatial and temporal resolution, facilitating large-scale estimation of shrub biomass [26,27]. Recent advancements in satellite remote sensing technology, along with the availability of data from satellites such as Landsat and Sentinel series, have expanded the possibilities in the field of biomass estimation [28,29]. However, utilizing only the spectral bands of remote sensing images in shrub biomass estimation has its constraints [30]. For example, the visible and near-infrared bands have a narrow wavelength range [31], making it challenging to capture the differences in factors closely related to shrub biomass, such as vegetation structure, leaf area index, and chlorophyll content. To overcome these limitations, many studies have used vegetation indices as important indicators for biomass estimation [32]. Vegetation indices are numerical values calculated from spectral bands that can reflect key information about vegetation growth status, chlorophyll content, leaf area index, and other related factors [33,34]. However, existing vegetation indices may encounter saturation issues and struggle to differentiate higher biomass shrubs in densely vegetated areas [35]. Therefore, it is imperative to evaluate which vegetation indices or spectral bands contribute most significantly to the estimation of shrub biomass, emphasizing the need to identify the optimal feature combination for accurate estimation.

Unmanned aerial vehicles (UAVs) are also essential in the estimation of shrub biomass due to their ability to provide detailed information and significant advantages [36,37]. Equipped with high-resolution sensors and advanced image processing techniques, UAVs can capture high-quality remote sensing imagery, offering comprehensive data on shrub vegetation [38,39]. They are capable of capturing extensive areas of shrub coverage from a high-altitude overhead perspective, enabling the rapid and accurate identification and extraction of shrub objects [37,40,41,42]. Additionally, UAVs have higher spatial resolution and flexibility compared to traditional aerial remote sensing techniques, allowing for the capture of finer details of shrub structure and characteristics [43]. This capability is crucial for biomass estimation, as the spatial heterogeneity and subtle variations of shrubs can significantly impact the accuracy of estimations. Furthermore, UAVs can enhance workflow efficiency, reduce workforce, and save time and costs compared to traditional field surveys and plot measurements [44]. However, there is limited research on integrating UAV-derived shrub biomass estimation at the satellite scale, highlighting the need for further investigation in this area.

Traditionally, biomass estimation models relied on statistical regression techniques to establish empirical relationships between field-measured biomass and various biophysical or remote sensing variables [45]. Over time, several modeling approaches have been developed to enhance the precision of estimating shrub biomass [37,46,47]. In recent times, machine learning methods have emerged as potent tools in this field, demonstrating superior performance compared to conventional modeling methods [48]. However, these models had limitations in capturing complex nonlinear relationships and accounting for spatial heterogeneity [49]. In contrast, machine learning methods, such as random forests, support vector machines, and neural networks, have gained traction due to their capability to handle intricate and nonlinear relationships [50,51]. Machine learning methods offer several advantages over traditional modeling techniques in shrub biomass estimation, particularly in effectively handling large volumes of remote sensing data, including multispectral and hyperspectral imagery, which provide rich spectral and spatial information for accurate biomass estimation [48,52,53].

This work aimed to propose a framework to estimate shrubland biomass at 10m over an arid and semi-arid mountain region based on multi-scale data from field measurements, UAV, Sentinel-1, Sentinel-2 and Landsat observations. The Helan Mountains in Ningxia province, China, were selected as the study area. Firstly, a land cover classification was conducted to identify the shrublands and other land cover types. Secondly, the prediction model of shrub biomass was developed in a Random Forest Regression (RFR) approach driven by different predicted variable datasets based on field measurements, UAV, and satellite images. The field measurement data was used to establish the allometric growth equation between shrub biomass and shrub structure parameters. Using the allometric equation, the shrub biomass was determined in UAV data. Using the UAV-based shrub biomass as inputs, the optimal satellite-based biomass estimation model was developed by comparing different predictor variables from Landsat, Sentinel-1, and Sentinel-2 satellites. Thirdly, with the best model, we created a map of the biomass distribution of shrubland in the Helan Mountains. The accuracy of the resultant map was evaluated over various ranges of shrub biomass or shrub coverage. Finally, we assessed the spatial characteristics of shrubland biomass based on the resultant biomass map and the auxiliary datasets of climate and topography.

2. Material and methods

2.1. Study area

The Helan Mountains are a remnant range of the Kunlun Mountains and are situated at the intersection of the Ningxia Hui Autonomous Region and the Inner Mongolia Autonomous Region. The main peak of the Helan Mountains in Ningxia has an elevation of 3,556 meters and covers an area of approximately 2,100 square kilometers. The mean annual temperature in the study area is around –0.8°C, with a mean annual precipitation of about 430 millimeters.

The study area represents a typical temperate arid mountain region with diverse vegetation types. It is mainly covered by evergreen coniferous forests, deciduous broad-leaved forests, evergreen coniferous shrubs, deciduous broad-leaved shrubs, typical grasslands, and desert grasslands, with slope direction and height influencing distribution patterns. Furthermore, the region has several shrub species, including Amygdalus pedunculata, Caragana korshinskii, and Mongolian almond.

Figure 1. (a, b) Locations of the Helan Mountain in China and Ningxia province, and (c) the distribution of ground truth samples from field measurements, UAV, and visual interpretation.

Figure 1. (a, b) Locations of the Helan Mountain in China and Ningxia province, and (c) the distribution of ground truth samples from field measurements, UAV, and visual interpretation.

2.2. Date

2.2.1. Landsat and Sentinel-1/2 data and pre-processing

Our study utilized all available Landsat 7/8/9, Sentinel-2, and Sentinel-1 satellite images from late July to early August 2023. These satellite datasets were conveniently accessed through Google Earth Engine. This study utilized Level-2 Landsat-7 Enhanced Thematic Mapper Plus (ETM+), Landsat-8 Operational Land Imager (OLI), and Landsat-9 OLI-2 data, which include visible, infrared, and shortwave infrared bands with a spatial resolution of 30m. This study also used Sentinel-2 multispectral instrument (MSI) Level-2A products, which provide red-edge band data additionally at 10m compared to Landsat. Sentinel-1 provides dual-polarization C-band Synthetic Aperture Radar (SAR) observations operating at 5.405 GHz (C band), including vertical receive (VV) and horizontal receive (VH) bands. Thus, our dataset includes three visible bands, three red-edge bands, one near-infrared band, two shortwave infrared bands, and two radar bands.

To identify the quality of the satellite images, we performed image processing using the quality assurance (QA) band and CFMask algorithm within the images [54]. Additionally, we conducted calibration processes to harmonize the Landsat and Sentinel-2 images using the approach of ordinary least squares (OLS) regression at 10m [55].

Based on the harmonized Landsat/Sentinel-2 datasets, various vegetation indices sensitive to shrub biomass [56] were calculated to select the optimal predictor variable sets for constructing the estimated model of shrub biomass. A total number of 28 vegetation indices were selected. These vegetation indices were computed at the GEE platform.

Table 1. Widely used vegetation indices (VIs) for the estimates of shrub biomass. R = Red, G = Green, B = Blue, RE = RedEdge1, NIR = Near-infrared, SWIR1 = Shortwave infrared1 and SWIR2 = Shortwave infrared2. ρ represents the surface reflectance of MODIS bands 2 (841–876 nm), σ = 0.5 × (NIR + Red).

No.	Index	Formula
1	Normalized Difference Vegetation Index [57]	NDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}$
2	Enhanced Vegetation Index 1 [58]	EVI = $\frac{2.5\times (\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+6\times \mathrm{R}-7.5\times \mathrm{B}+1}$
3	Land Surface Water Index [59]	LSWI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}$
4	Difference Vegetation Index [60]	DVI = NIR − R
5	Green Normalized Difference Vegetation Index [61]	GNDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{G}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{G}}$
6	Vegetation Index green [62]	VIgreen = $\frac{\mathrm{G}-\mathrm{R}}{\mathrm{G}+\mathrm{R}}$
7	Infrared Simple Ratio [63]	ISR = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}$
8	Moisture Stress Index [64]	MSI = $\frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{N}\mathrm{I}\mathrm{R}}$
9	Ratio Vegetation Index [65]	RVI = $\frac{\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}}$
10	Simple Ratio [66]	SR = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}$
11	Enhanced Vegetation Index 2 [67]	EVI2 = $\frac{2.5\times (\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+2.4\times \mathrm{R}+1}$
12	Modified Simple Ratio [68]	MSR = $\frac{\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}-1}{\sqrt{\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}+1}}$
13	Optimized Soil-Adjusted Vegetation Index [69]	OSAVI = (1+L) × $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}+\mathrm{L}}$, L was set to 0.16
14	Renormalized Difference Vegetation Index [70]	RDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\sqrt{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}}$
15	Soil Adjusted Vegetation Index [71]	SAVI = (1+L) ×$\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}+\mathrm{L}}$, L=0.5
16	Soil Adjusted Vegetation Index2 [72]	SAVI2 = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\raisebox{1ex}{$b$}\!\left/ \!\raisebox{-1ex}{$a$}\right.}$, b was set to 0.025 and a to 1.25
17	Stress-related Vegetation Index 1 [73]	STVI1 = $\frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1\times \mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}}$
18	Stress-related Vegetation Index 2 [73]	STVI2 = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2\times \mathrm{R}}$
19	Stress-related Vegetation Index 3 [73]	STVI3 = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1\times \mathrm{R}}$
20	Red Edge Normalized Difference Vegetation Index [74]	RENDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}}$
21	Anthocyanin Reflectance Index [75]	ARI = $\frac{1}{\mathrm{G}}$ - $\frac{1}{\mathrm{R}\mathrm{E}}$
22	Vogelmann Red Edge Index [76]	VREI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}\mathrm{E}}$
23	Radar ratio vegetation index	Ratio = $\frac{\mathrm{V}\mathrm{V}}{\mathrm{V}\mathrm{H}}$
24	Radar Difference Vegetation Index	Difference = VV - VH
25	Radar normalized difference vegetation index RNDVI	RNDVI = $\frac{\mathrm{V}\mathrm{V}-\mathrm{V}\mathrm{H}}{\mathrm{V}\mathrm{V}+\mathrm{V}\mathrm{H}}$
26	Near-infrared reflectance of vegetation [77]	NIRv = (NDVI – C) × ρ (C = 0.08)
27	kernel NDVI [78]	kNDVI = tanh(($\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{2\sigma }$)2)
28	Normalized Difference Phenology Index [79]	NDPI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-(0.74\times \mathrm{R}\mathrm{e}\mathrm{d}+0.26\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+(0.74\times \mathrm{R}\mathrm{e}\mathrm{d}+0.26\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R})}$

Table 1. Widely used vegetation indices (VIs) for the estimates of shrub biomass. R = Red, G = Green, B = Blue, RE = RedEdge1, NIR = Near-infrared, SWIR1 = Shortwave infrared1 and SWIR2 = Shortwave infrared2. ρ represents the surface reflectance of MODIS bands 2 (841–876 nm), σ = 0.5 × (NIR + Red).

No.	Index	Formula
1	Normalized Difference Vegetation Index [57]	NDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}$
2	Enhanced Vegetation Index 1 [58]	EVI = $\frac{2.5\times (\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+6\times \mathrm{R}-7.5\times \mathrm{B}+1}$
3	Land Surface Water Index [59]	LSWI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}$
4	Difference Vegetation Index [60]	DVI = NIR − R
5	Green Normalized Difference Vegetation Index [61]	GNDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{G}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{G}}$
6	Vegetation Index green [62]	VIgreen = $\frac{\mathrm{G}-\mathrm{R}}{\mathrm{G}+\mathrm{R}}$
7	Infrared Simple Ratio [63]	ISR = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}$
8	Moisture Stress Index [64]	MSI = $\frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{N}\mathrm{I}\mathrm{R}}$
9	Ratio Vegetation Index [65]	RVI = $\frac{\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}}$
10	Simple Ratio [66]	SR = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}$
11	Enhanced Vegetation Index 2 [67]	EVI2 = $\frac{2.5\times (\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+2.4\times \mathrm{R}+1}$
12	Modified Simple Ratio [68]	MSR = $\frac{\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}-1}{\sqrt{\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}+1}}$
13	Optimized Soil-Adjusted Vegetation Index [69]	OSAVI = (1+L) × $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}+\mathrm{L}}$, L was set to 0.16
14	Renormalized Difference Vegetation Index [70]	RDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\sqrt{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}}$
15	Soil Adjusted Vegetation Index [71]	SAVI = (1+L) ×$\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}+\mathrm{L}}$, L=0.5
16	Soil Adjusted Vegetation Index2 [72]	SAVI2 = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\raisebox{1ex}{$b$}\!\left/ \!\raisebox{-1ex}{$a$}\right.}$, b was set to 0.025 and a to 1.25
17	Stress-related Vegetation Index 1 [73]	STVI1 = $\frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1\times \mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}}$
18	Stress-related Vegetation Index 2 [73]	STVI2 = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2\times \mathrm{R}}$
19	Stress-related Vegetation Index 3 [73]	STVI3 = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1\times \mathrm{R}}$
20	Red Edge Normalized Difference Vegetation Index [74]	RENDVI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}}$
21	Anthocyanin Reflectance Index [75]	ARI = $\frac{1}{\mathrm{G}}$ - $\frac{1}{\mathrm{R}\mathrm{E}}$
22	Vogelmann Red Edge Index [76]	VREI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}\mathrm{E}}$
23	Radar ratio vegetation index	Ratio = $\frac{\mathrm{V}\mathrm{V}}{\mathrm{V}\mathrm{H}}$
24	Radar Difference Vegetation Index	Difference = VV - VH
25	Radar normalized difference vegetation index RNDVI	RNDVI = $\frac{\mathrm{V}\mathrm{V}-\mathrm{V}\mathrm{H}}{\mathrm{V}\mathrm{V}+\mathrm{V}\mathrm{H}}$
26	Near-infrared reflectance of vegetation [77]	NIRv = (NDVI – C) × ρ (C = 0.08)
27	kernel NDVI [78]	kNDVI = tanh(($\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{2\sigma }$)2)
28	Normalized Difference Phenology Index [79]	NDPI = $\frac{\mathrm{N}\mathrm{I}\mathrm{R}-(0.74\times \mathrm{R}\mathrm{e}\mathrm{d}+0.26\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+(0.74\times \mathrm{R}\mathrm{e}\mathrm{d}+0.26\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R})}$

2.2.2. Land cover data

The European Space Agency (ESA) WorldCover dataset provides land cover types at 10m global, which has been generated based on Sentinel-1/2 images with 11 land cover types [80]. We generated 500 samples randomly and acquired the land cover type for each sample point. Subsequently, we examined the classification accuracy of each sample by visual interpretation using high-resolution images from Google Earth and ground in-situ photos. Finally, a total of 354 samples have been obtained, including 97 shrublands, 86 grasslands, 112 barren lands, and 59 forests. According to the MODIS land cover dataset, the land with a shrub cover of more than 10% was classified as shrubland for this study [81]. These samples were separated as training and validation points randomly with a ratio of 7:3 used for shrubland mapping.

2.2.3. Ground-based measurements and UAV field images

We conducted two comprehensive field surveys in the Helan Mountains during late April 2023 and from late July to early August in 2023. Thirteen shrubs in the study area were subjected to destructive sampling and structural parameter measurements including length (CL), width (CW), height (CH), crown area (CA), and volume (CV). These measurements were taken to develop the optimal allometric growth equation for estimating shrub biomass. Furthermore, we established 24 sampling plots that were each 10m × 10m in size, where we measured the structural parameters of every shrub. The shrub samples were later processed in the laboratory for drying and weighing to obtain the biomass data.

Table 2. Field samples, including structural parameters and biomass data of each shrub.

ID	CL (m)	CW (m)	CH (m)	CA (m2)	CV (m3）	Number of branches	Single branch biomass (g)	Single plant biomass (g)
1	1.62	1.02	1.34	1.3	1.74	6	487.81	2926.88
2	1.74	1.55	1.83	2.12	3.88	13	243.91	3170.79
3	1.26	1.02	0.73	1.01	0.74	33	40.95	1351.46
4	0.72	0.58	0.5	0.33	0.16	1	183.43	183.43
5	0.51	0.48	0.45	0.19	0.09	1	513.86	513.86
6	1.15	0.96	0.84	0.87	0.73	10	134.43	1344.3
7	0.42	0.42	0.6	0.14	0.08	1	185.91	185.91
8	1.63	1.58	0.74	2.02	1.5	9	226.38	2037.45
9	0.88	0.99	1	0.68	0.68	9	121.95	1097.58
10	1.96	2.08	1.32	3.2	4.23	40	38.48	1539.07
11	1.32	1.47	1.11	1.52	1.69	15	144.95	2174.3
12	1.05	0.98	0.59	0.81	0.48	3	495.86	1487.58
13	0.89	0.68	0.75	0.48	0.36	5	51.95	259.77

Table 2. Field samples, including structural parameters and biomass data of each shrub.

ID	CL (m)	CW (m)	CH (m)	CA (m2)	CV (m3）	Number of branches	Single branch biomass (g)	Single plant biomass (g)
1	1.62	1.02	1.34	1.3	1.74	6	487.81	2926.88
2	1.74	1.55	1.83	2.12	3.88	13	243.91	3170.79
3	1.26	1.02	0.73	1.01	0.74	33	40.95	1351.46
4	0.72	0.58	0.5	0.33	0.16	1	183.43	183.43
5	0.51	0.48	0.45	0.19	0.09	1	513.86	513.86
6	1.15	0.96	0.84	0.87	0.73	10	134.43	1344.3
7	0.42	0.42	0.6	0.14	0.08	1	185.91	185.91
8	1.63	1.58	0.74	2.02	1.5	9	226.38	2037.45
9	0.88	0.99	1	0.68	0.68	9	121.95	1097.58
10	1.96	2.08	1.32	3.2	4.23	40	38.48	1539.07
11	1.32	1.47	1.11	1.52	1.69	15	144.95	2174.3
12	1.05	0.98	0.59	0.81	0.48	3	495.86	1487.58
13	0.89	0.68	0.75	0.48	0.36	5	51.95	259.77

In late April 2023, we acquired a UAV image of the Helan Mountain region with the help of a DJI Phantom 4 Pro V2.0. The image covered an expansive area of 588m ×196m, with a spatial resolution of 0.06m. The Phantom 4 Pro V2.0 UAV provides six bands, including true color (RGB), blue, green, red, red edge, and near-infrared (NIR).

2.2.4. Auxiliary datasets

Relevant climate and terrain datasets, including precipitation, temperature, aridity index, and elevation datasets for our subsequent analyses. The precipitation data was obtained using the Global Precipitation Measurement (GPM) v6 dataset (Precipitation Processing System (PPS) at NASA GSFC, 2019), with a spatial resolution of 11,132m. It provides global rainfall and snowfall observations every three hours by integrating data from all sensors within the GPM system through calibration, merging, and interpolation. The temperature data used was the ECMWF Climate Reanalysis dataset [83], which provides daily observation data at a spatial resolution of 11,132m. The aridity index data was the Global Aridity Index and Potential Evapotranspiration (ET0), which calculates the aridity index based on the mean annual precipitation and evapotranspiration [84], with lower values indicating more severe drought conditions. The spatial resolution of the aridity index dataset is 30 arc-seconds. The elevation data was obtained from the 30m NASA SRTM Digital Elevation dataset [85]. These datasets were obtained through the Google Earth Engine platform and resampled to a spatial resolution of 1,000m to ensure data compatibility.

2.3. Method

Figure 2 shows an overview of the workflow for biomass estimation in the shrublands. It included four main procedures. Firstly, a sample dataset was constructed based on field surveys, ground measurements, UAVs, and Google Earth images. The dataset included samples of land cover types, shrub structure parameters, and biomass. Secondly, a land cover classification was performed to extract the distribution of shrubland. Afterward, the optimal model for shrub biomass estimates was built by comparing several models using different predictor variables derived from individual spectral bands, spectral vegetation indices, and backscatter values of VV and VH. Finally, a comprehensive analysis of the spatial distribution of biomass for the shrubland was conducted. The subsequent paragraphs provide detailed information for each of these steps.

2.3.1. Mapping of shrubland and other land cover types

Random Forest is a well-liked machine-learning method that combines the characteristics of decision trees and randomness [86,87]. It demonstrates good skills in handling large-scale and high-dimensional datasets and is also able to evaluate each feature’s significance in terms of selecting important features.

For our research, we utilized GEE to gather 11 bands of data to classify various land cover types. These bands consist of three visible bands, three infrared bands, one near-infrared band, two shortwave infrared bands, and two radar bands. The band information of each sample point and the associated land cover categorization results were then entered into the Random Forest model. 70% of the sample points were utilized to train the model, while the remaining sample points were reserved for subsequent accuracy assessment. Our evaluation metrics included accuracy, recall, and F1 score, all calculated using Python (version 3.7) and the sklearn package (version 0.0 post2) for constructing the Random Forest model and computing accuracy metrics.

$\mathrm{Accuracy}=\frac{(\mathrm{True}\mathrm{Positives}+\mathrm{True}\mathrm{Negatives})}{(\mathrm{True}\mathrm{Positives}+\mathrm{True}\mathrm{Negatives}+\mathrm{False}\mathrm{Positives}+\mathrm{False}\mathrm{Negatives})}$

$\mathrm{Recall}=\frac{\mathrm{True}\mathrm{Positives}}{(\mathrm{True}\mathrm{Positives}+\mathrm{False}\mathrm{Negatives})}$

$\mathrm{F}1\mathrm{score}=\frac{2\times (\mathrm{Precision}\times \mathrm{Recall})}{(\mathrm{Precision}+\mathrm{Recall})}$

2.3.2. Calculation of biomass at the shrub level

The structural parameters of shrubs closely reflect their growth status and biomass accumulation, making them important indicators for estimating shrub biomass. By conducting field measurements of shrub length, width, and height, we can calculate the shrub's ground area and volume, both of which are positively correlated with shrub biomass [56]. Generally, shrub area and volume directly reflect biomass values. The ground area represents the horizontal structure of the shrub, while the volume represents its horizontal and vertical structural characteristics [56]. Studies have shown a linear link between the biomass of shrubs and CA [56].

To account for the diverse range of shrub species found in different regions, our study produced allometric equations for shrub development, with a specific focus on the Helan Mountains. We utilized seven samples in each iteration to fit the canopy area (CA) and crown volume (CV), with the remaining six samples being used for validation. We conducted a total of $5\times C_{13}^7$ iterations, using the coefficient of determination (R²), mean absolute error (MAE), and root mean square error (RMSE) as accuracy assessment criteria. The final selected allometric equations for shrub biomass estimation based on area and volume are as follows: AGB = $e^{(7.3+0.41\times \mathrm{l}\mathrm{n}(CV)}$(R² = 0.89), AGB = 1060.93×CA (R² = 0.9).

2.3.3. Calculation of biomass at the UAV plot level

UAV provides an approach to obtain the shrub biomass at the plot level conveniently, which can well alleviate the estimation uncertainty caused by the spatial mismatch between ground sampling and satellite observations [56]. In this study, the shrub was identified firstly by a random forest classification algorithm based on the true color UAV image (Figure 3a). The shrub (number=106) and non-shrub (number=229) samples were collected by visual interpretation of the UAV imagery using a random sampling approach. The RGB band values of each sample in the UAV imagery were extracted as input features for building the random forest classification model. Thus, a binary image of shrub and non-shrub was produced at the UAV level. The model fitness showed that the accuracy, recall, and F1 score of the UAV-based shrub map all exceeded 0.9.

Subsequently, the grids with a size of 10m × 10m were created using the UAV-based shrub and non-shrub map to match the spatial resolution of satellite images. This process was conducted using the “create fishnet” function in ArcMap 10.4. The grids were perfectly aligned with the pixel boundaries of the satellite image (Figure 3c). Within each grid, the crown area (CA) of shrubs was calculated using the shrub pixels. Finally, the total shrub biomass within each grid was estimated by applying a biomass allometric equation based on the CA [56]. To reduce the impacts caused by rugged terrain and clouds, the grids with black or bright backgrounds, steep slopes, or clouds were removed. A total of 86 UAV-based biomass samples were obtained for further analysis.

2.3.4. Satellite-based modeling of shrubland biomass at the regional level

For biomass estimation and accuracy assessment, we selected the Random Forest Regression (RFR) model and employed the leave-one-out cross-validation (LOOCV) method. The RFR model calculates the final score as the mean score of all decision trees rather than through voting, resulting in better predictive performance and higher accuracy than the RF model [88]. We trained the model n times (n=89), where n-1 samples were used for training, and the remaining sample was used for testing. The model with the highest accuracy was selected as the final output, allowing us to make full use of the samples [56].

The model's performance is not positively impacted by every feature [89]. Certain characteristics may cause little variations in samples, which might lower the classification accuracy [90]. For feature selection, we so used the Python RFECV and KFold libraries. To find the optimal feature combination, the Random Forest Regression model's features are combined using the RFECV package. To ensure optimal accuracy for the final feature combination and minimize the possibility of accidental selection, the KFold package is utilized for cross-validation to validate the chosen feature combinations. We then constructed three models: one with only VV and VH and spectral bands (SB), another with only vegetation indices (VI), and a third with all available variables (SBVI). To identify the ideal feature set for each model, feature combination selection was done. Finally, by comparing the accuracy of different models, we identified the optimal feature combination for constructing the biomass model.

2.3.5. Regional implementation of the satellited-based shrubland biomass model

The RFR model with the best predictor variables was selected to estimate the biomass of shrubland in the study area based on the Landsat, Sentinel-1, and Sentinel-2 images. Through the rule of LOOCV, we generated a series of shrubland biomass maps (number=89). Based on these maps, we obtained the resultant biomass map of shrubland by calculating the mean value of each pixel. Meanwhile, the relevant map of standard deviation (SD) for shrubland biomass was also produced to present the potential uncertainty. Subsequently, we analyzed the uncertainty of the resultant biomass map along with the magnitudes of shrubland biomass and coverage. To obtain the shrubland coverage map, we used the shrubland and non-shrubland classification map generated in Section 2.3.1. This 10m binary map was resampled into a shrubland coverage map at a spatial resolution of 1,000m by the mean method. Then, the relationship between the SD map and the shrubland coverage map was examined at 1,000m by a spatial overlay analysis.

2.3.6. Characteristics of the shrubland biomass

To examine the geographical patterns of shrubland biomass, we overlayed the resultant biomass map with different datasets of precipitation, temperature, aridity index, and elevation (Table 3). These datasets were obtained through the GEE platform and were resampled to a consistent spatial resolution of 1,000m. The mean shrubland biomass was calculated within different ranges of the controlling factors.

3. Results

3.1. Mapping of shrubland and other land-cover types

We created a land cover map for the Helan Mountains in 2023 using an RF model constructed with 11 bands. We have identified four land cover types: forest, barren land, grasslands, and shrublands. Figure 4 displays the result map, where only the shrublands were used for biomass estimation. It can be observed that shrublands are predominantly distributed in the southeastern plains and northeastern parts, while forests dominate the southwestern region. Additionally, there is relatively less grassland and more barren land in the study area, which can be attributed to the overall climate type being a typical temperate arid climate with low rainfall. The accuracy assessment on the shrubland and non-shrubland map was shown in Table 4, presenting an overall classification accuracy of 0.91, a recall of 0.92, and an F1 score of 0.92.

3.2. Selection of the best model and variable importance analysis

We conducted three models driven by different variable sets: VV, VH, and spectral bands (SB), another with only vegetation indices (VI), and a third with all available variables (SBVI). The RFECV and KFold packages were used to determine the optimal variables for each model. The best variables for each model are shown in Table 5. The R², MAE, and RMSE for each of the three models are shown in Figure 5. The results suggested that the SBVI model achieved the best performance, followed by the SB and VI models. Therefore, we selected the SBVI model as our final model for estimating shrub biomass. It also showed that a significant number of vegetation indices related to the near-infrared (NIR) band are incorporated into the model, indicating the high explanatory power of NIR-based vegetation indices for shrub biomass.

The accuracy was assessed based on shrub biomass samples across different groups of shrub coverages and biomass. As shown in Figure 6a, the ratio of RMSE to mean shrub cover was calculated to present the error of one percent coverage of shrub (EOPC). EOPC decreases with increasing shrub cover. Similarly, Figure 6b shows the ratio of RMSE to mean biomass, presenting the error of one unit biomass (EOUB). EOUB decreases as the shrub biomass increases. In terms of R², the lowest values were observed within the ranges of 20-30% shrub cover and 20-30 kg/m² shrub biomass. These results indicate good predictive capabilities in high and low shrub cover and shrub biomass ranges, while the estimation accuracy is relatively lower in the moderate shrub cover and biomass ranges.

To determine the contributions of four features to the model, we conducted a sensitivity analysis (Figure 6c). It can be observed that SAVI2 has the most significant impact on model accuracy. Removing SAVI2 resulted in a decrease of R² by 0.13 and an increase of RMSE by 9.53 kg/m². In addition, VH had the least impact on the model, as removing VH only led to a decrease of 0.03 in R² and an increase of 1.93 kg/m² in RMSE. The findings showed that the four characteristics' contributions—SAVI2, DVI, VV, and VH—ranked from high to low.

3.3. Biomass mapping in the shrubland

We used the SBVI model to estimate the shrubland biomass in the Helan Mountains (Figure 7). The results showed the minimum value of shrub biomass was 10.5 kg/m², while the maximum value was 54.73kg/m². Despite having more bushes than other vegetation, the biomass is often lower in the southeast. This is because the shrubs in this area are predominantly low-lying shrubs. By contrast, the central Helan Mountains have greater heights and are mostly characterized by tall shrub, which contributes to usually larger levels of shrub biomass.

3.4. Distribution of shrub biomass under different factors

The mean biomass of shrubs was determined by dividing temperature, precipitation, aridity index, and elevation into several categories. In Figure 8a, as shrub biomass grows with increased precipitation up to 350 mm, then it begins to decline. Figure 8b displays the distribution of shrub biomass across different ranges of mean annual temperature, revealing no significant trend of shrub biomass increase or decrease with rising temperatures. As can be seen in Figure 8c, shrub biomass and the aridity index show a positive link up until the index hits 0.2, at which point the two show a negative correlation. Nevertheless, no discernible pattern of shrub biomass fluctuation with elevation variations is seen in Figure 8d.

4. Discussion

4.1. Algorithms

The integration of unmanned aerial vehicle (UAV) derived shrub biomass estimates into satellite imagery at a larger scale has significant implications for ecological research and environmental monitoring. By establishing a grid with a spatial resolution that matches that of satellite imagery, we were able to estimate shrub biomass within each grid cell. Furthermore, the utilization of UAV imagery provided a cost-effective and efficient method for estimating shrub biomass, which resulted in time and resource savings, as well as improved spatial coverage and accuracy of shrub biomass estimation. However, it is important to acknowledge that there are certain limitations and considerations associated with the integration of UAV and satellite data. The accuracy and reliability of the estimated shrub biomass at the satellite scale are influenced by factors such as the spatial and spectral resolution of the satellite imagery, atmospheric conditions, and the transferability of the allometric equations. Additionally, the temporal difference between the UAV and satellite imagery acquisitions may introduce uncertainties in capturing the dynamic nature of shrub biomass. Consequently, further research is necessary to explore the potential of integrating multi-scale remote sensing data to achieve more comprehensive and accurate estimation and monitoring of shrub biomass.

In this study, we conducted land cover classification in the Helan Mountains region using satellite data. Additionally, we developed three models and selected the model with the highest accuracy for shrub biomass estimation in the shrublands. The selection of the model with the highest accuracy for shrub biomass estimation is a crucial outcome of our study. The precise estimation of shrub biomass enabled by the selected model has several potential applications. Firstly, it can contribute to ecological research and conservation efforts by providing insights into the growth status and productivity of shrub populations [91]. Understanding the spatial distribution and temporal trends of shrub biomass is crucial for assessing ecosystem health, biodiversity, and carbon storage [92,93]. Secondly, our findings have practical implications for land management and restoration projects. Accurate shrub biomass estimation can assist in decision-making processes related to land-use planning, habitat restoration, and the preservation of biodiversity hotspots.

We discovered that errors decreased as vegetation coverage and estimated biomass increased, suggesting that predicting shrub biomass in areas with dense vegetation and lower biomass levels presents challenges. In contrast, our model demonstrated higher predictive accuracy in areas with higher vegetation coverage and biomass. These findings are significant for comprehending the dynamics of shrub ecosystems and can contribute to more effective management and conservation strategies.

4.2. Spatial characteristics of biomass in shrubland

The study found that as precipitation levels increased to 350 mm, shrub biomass growth also increased, suggesting an optimal range for precipitation in the area. However, excessive rainfall beyond a certain threshold may lead to a shift in dominant species from shrubs to other plant types within the region, thereby impacting the growth environment of shrubs and resulting in a reduction in their abundance [94,95,96]. In terms of temperature, there was no discernible pattern of shrub biomass increase or decrease as temperature rose. This suggests that mean annual temperature may not be the primary factor affecting shrub biomass variation in the Helan Mountains. On the other hand, an intriguing relationship was observed between shrub biomass and the aridity index. When the index was below 0.2, there was a positive correlation between shrub biomass and moderately dry conditions. However, beyond an aridity index of 0.2, a negative correlation was observed, indicating that drought conditions are more conducive to shrub growth [96,97]. As for elevation, no significant trend of shrub biomass variation was observed with elevation changes, indicating that elevation may not be a primary factor influencing shrub biomass distribution in the Helan Mountains. Previous studies have suggested a correlation between shrub biomass distribution in the southwestern region of China and the four selected influencing factors in our research [97]. Our study did not find significant effects of temperature and elevation on shrubs, which may be attributed to regional variations. Therefore, future research could further investigate the interactions between these factors to unravel the complex relationships between elevation, environmental conditions, and shrub biomass dynamics.

Although our research yielded promising outcomes, it is essential to acknowledge certain limitations. Firstly, the precision of estimating shrub biomass may be impacted by factors such as soil moisture, topographical variances, and the presence of understory vegetation [8,98]. Future investigations could explore the incorporation of supplementary environmental factors to enhance model efficacy. Secondly, the generalization of our findings to other regions or diverse shrub species should be approached cautiously, as the associations between vegetation indices and shrub biomass may differ. Expanding the study area could provide insights into the applicability of our model. Lastly, our results demonstrated strong performance in depicting the spatial distribution of shrub biomass in the Helan Mountains, contributing significantly to various aspects including vegetation management, shrub surveys, and ecological conservation in the area.

5. Conclusion

In conclusion, this study effectively utilized satellite data to classify shrublands and other land cover types in the Helan Mountains region. Additionally, we georeferenced the shrub biomass estimates computed from the UAV imagery onto satellite imagery. Three models were developed, and the most accurate model was chosen for shrub biomass estimation in the shrublands. We have created a biomass distribution map of shrubs in the Helan Mountains using the most accurate model. Following that, we conducted an accuracy evaluation on various levels of shrub biomass and coverage based on the samples and the shrub biomass distribution map. Finally, an analysis of the regional distribution of shrub biomass concerning temperature, aridity index, elevation, and precipitation was done. This information is essential for evaluating ecosystem health, carbon storage, and biodiversity conservation.