# The Spatiotemporal Pattern of the Aerosol Optical Depth (AOD) on the Canopies of Various Forest Types in the Exurban National Park: A Case in Ningbo City, Eastern China

Authors: Yufeng Chi; Shudi Zuo; Yin Ren; Kaichao Chen
Journal: Advances in Meteorology
 (2019)

---

## Abstract

The spatiotemporal distribution pattern of the aerosol optical depth (AOD) is influenced by many environmental factors, such as meteorological condition changes, atmospheric pollution, and topographic changes. Understanding the relationship between the vegetation land cover and the AOD would favor the improvement of forest ecosystem services. This quantitative research integrated remote sensing and ground survey data and used spatial statistical methods to explore the drivers that influence the AOD of the exurban national forest park and analyze the differences between various forest types. The driver analysis was carried out in the hot (Z ≥ 1.64) and cold (Z ≤ −1.64) spots of AOD in 2010 and 2017. Our results showed that (1) the forest type was proved to be the main factor contributing to the AOD pattern and (2) from 2010 to 2017, the average growth rate of broad-leaved forest, coniferous forest, bamboo, and shrub in hot spots was significantly higher than that in cold spots, while there was no significant difference in the mixed forest. The average growth rate of biomass densities of bamboo, coniferous forest, and mixed forest were higher than that of the shrub and broad-leaved forest. These findings provided the guidance for the rational allocation of tree species to increase the biomass and improve the ecosystem service values of forest parks.


---

## Body

## 1. Introduction

There are two major sources of atmospheric aerosol: natural and human emissions, which is a system comprising atmospheric medium with mixed solid and liquid particles. The composition is complex and diverse, including various trace metals, inorganic oxides, sulfates, nitrates, and oxygen-containing organic compounds [1–4]. The composition of the urban atmosphere affected by the various pollution sources changes greatly. However, the composition in the exurbs is relatively stable, which is generally related to the local soil composition. As a crucial natural provider of multiple ecosystem services, the quantitative assessment and improvement measures of exurban national forest parks are of increasing concern to urban managers [5, 6]. Human health and vegetative growth are inevitably affected by aerosols. Aerosols have important effects on air quality, visibility, acid deposition, precipitation, atmospheric radiation balance, and chemical reactions in the stratosphere and troposphere. Thus, understanding the spatiotemporal pattern of the aerosol in the exurban national park and its relationship to the forest type on the ground will be helpful to improve ecological system services provided by the forest.AOD is an important index of atmospheric turbidity [7, 8]. It has been widely used in related studies to improve air quality because of certain linear relationships with the PM 2.5 concentration on the ground. The smaller the AOD, the cleaner the atmosphere and vice versa. Furthermore, it is vital to study the mechanisms influencing the distribution pattern of aerosols. Previous studies have shown that AOD may overestimate aerosol concentration when the relative humidity (RH) is higher than 40% and can better reflect dry aerosol concentration when the RH is lower than 40% [9]. Quantitative analysis of influencing factors of aerosol spatial distribution requires the integration of remote sensing images and field surveys to identify influential factors, which have strong effects in clustering regions (hot or cold spots) differentiated by spatial statistical analyses [10–13]. Besides the influence of human activities, the meteorological condition, topography, soil, and forest vegetation also have an interactive impact on aerosol concentration. The relationship between vegetation and aerosol is mainly reflected in the direct and indirect effects of aerosols on plants. The manifestation of direct influence is the aerosol coverage on plant leaves, affecting plant respiration, stomatal conductance, and utilization rate of sunlight. The indirect effect refers to the influences of aerosols on precipitation and temperature of the atmosphere as well as the scattering of sunlight, etc., further resulting in the effects on the utilization of light, water, and heat by plants [14]. Natural and anthropogenic aerosols mixed in the atmosphere mainly have indirect effects on plant growth, but less direct effects. The influence of aerosols on net primary production (NPP) may be positive, neutral, or negative [15]. Although certain concentrations of aerosols inhibit plant growth, the light scattering caused by aerosols disperses direct sunlight to every layer of leaves, including the leaves under canopy. This phenomenon is named as diffusion radiation fertilization effect, which has positive influences on vegetation [16].At present, the research on the relationship between AOD and plant growth is mainly divided into two spatial scales: the sample scale and a large regional scale [17–19]. Most of the recent studies were based on the quantitative experiment of ground-based remote sensing and the small-scale sample points, where the relationships between the AOD and a particular species or some species were detected by a statistical approach, On small-scale sample points, it has been proved that higher AOD concentration significantly promoted plant growth [20]. A part of the studies used global or large-scale model simulations to measure the effects caused by AOD diffuse radiation and showed that the diffuse radiation improved the efficiency of photosynthesis [21, 22]. However, in order to enhance regional ecosystem service capacity through forest management, these study results may bias the application guidance [23, 24]. On the one hand, sample plots may be difficult to reflect the aggregation effects of different tree species in the forest park on the spatial pattern of aerosol at the regional scale. On the other hand, a wide range of meteorological digital model simulations do not reflect the characteristics of vegetation properties inside the forest park. Furthermore, the AOD spatial resolution of existing satellite remote sensing products are 500 m or coarser, making it difficult to fully describe forest property details [25]. The ground forest management planning inventory (FMPI) data required the remote sensing images with high spatial resolution, such as 30 m, to support the mechanism analysis. In this study, we took the exurban Siming Mountain National Forest Park of Ningbo City in the Yangtze River Delta as an example. To conduct the analysis, we firstly integrated the multisource data composed of the FMPI, Landsat remote sensing data, road network data, meteorological data, and demographic data. Then, we used a geographical detector to analyze the heterogeneity of the AOD spatial-temporal distribution pattern in 2010 and 2017 and its main influencing factors. Finally, a statistical analysis was applied to find the relationship of AOD on different types of forest canopies.The purpose of this study is to find out the main influencing factors of aerosol distribution in exurban forest parks and to explore the differences in biomass growth of different forest types in aggregation areas with different concentrations. We tried to answer two questions: (1) in addition to human emission sources, what are the main ecological factors influencing exurban forest park AOD? and (2) do the AOD concentrations cause differences in biomass growth rates among different types of forests in the study area? The research results can provide guidance for improving the ecosystem service value of forest parks and provide reference for the further study of mechanisms.


## 2. Materials and Methods

### 2.1. Overview

We approached this research in the following steps: First, to create a database that integrates multiple types of GIS and remote sensing (RS) data (including Landsat 8 OLI images, Landsat 7-ETM+ images, FMPI, MODIS data, OpenStreetMap data, and Landscan data). Second, to preprocess the database. There were two substeps included: (1) replenishment of missing remote sensing pixel values, caused by cloud occlusion, through the Kriging interpolation which is a geostatistical method for predicting missing spatial information (AOD, MODIS data) [26, 27] and (2) inversion of the AOD by Landsat data. Third, the hot spots and cold spots of AOD in this area were identified by the high Z values (Z ≥ 1.64) and low Z values (Z ≤ −1.64) respectively, which were calculated through global Moran’s I and Getis-Ord Gi∗ [28, 29]. The optimal threshold distance was calculated using the incremental spatial autocorrelation module. Fourth, the factor detector of the Geodetector model was used to quantify the impact contribution of different factors on the AOD. Lastly, we calculated and analyzed the biomass changes of different forest types in cold and hot spots from 2010 to 2017.


### 2.2. Study Area

Siming Mountain National Forest Park is in Ningbo City, Zhejiang Province, China (Figure1). The location is longitude 120°59′20″E–121°25′16″E and latitude 29°31′39″N–29°59′20″N. The park is subject to the subtropical monsoon climate, which is characterised by 1,277 km2, 976 m highest altitude, 77.5% forest coverage, and yellow soil. The annual average temperature is 10ºC, and the annual average rainfall is 1,800 mm.Figure 1
Map of study area. The inset shows the location in China.


### 2.3. Multisource Dataset and Preprocessing

Three types of data in two years, which included the year 2010 and 2017, were collected. The first type of data was the Landsat satellite remote sensing data for AOD retrieval with a spatial resolution of 30 m. To ensure representation of the data, we selected the two nearest meteorological stations in the China Meteorological Data Network (http://data.cma.cn) (Nos. 58467 and 58562, respectively) to screen for wind speed and air pressure that were within a 5% error limitation of the mean value in both 2010 and 2017. Afterwards, we collected eight remote sensing image data for the four seasons in 2010 and 2017 (row number 118, column number 39). The remote sensing images and data from the four phases of the ETM sensor and the four phases of the OLI sensor were chosen based on the cloud cover percentage. The less cloudy ones were selected, which were the Landsat 7 ETM images on January 1, April 23, July 28, and December 3, 2010, and Landsat 8 OLI image on April 2, August 24, and October 27, 2017, and January 15, 2018 (Table 1). The cloud covers were 0.01%, 6.62%, 24.63%, 0.10%, 0.37%, 0.26%, 14.27%, and 1.27%, respectively. The average wind speeds at the two meteorological stations were (2.1 and 3.2 m/s), (4.1 and 3.1 m/s), (2.1 and 4.7 m/s), (3.6 and 3.5 m/s), (4.2 and 4.1 m/s), (2.2 and 2.3 m/s), (5.1 and 6.1 m/s), and (5.5 and 5.5 m/s) with relative humidity (75 and 77%), (88 and 86%), (90 and 91%), (88 and 85%), (79 and 82%), (85 and 90%), (90 and 90%), and (77 and 78%) in the years 2010 and 2017, correspondingly.Table 1
Information about the GIS and RS data in this study.
Type
Data
Time
GIS
FMPI
2010
2016
OSM
2017
LSP
2017
RS
Landsat 8 OLIRow 118/path 39
April 2, 2017
August 24, 2017
October 27, 2017
January 15, 2018
RS
Landsat 7 ETMRow 118/path 39
January 1, 2010
April 23, 2010
July 28, 2010
December 3, 2010
MODIS TPWV
The time of MODIS data set is the same as that of Landsat data set
MODIS LAI
MODIS LST
MODIS NDVI
MOD04_3KThe second type of data was the FMPI, which is obtained from in situ investigation and observation by the national forestry bureau and its affiliations in China. These data provide forest ground survey results on a regional scale. The attribute database contains information on forest patch area, diameter of breast height (DBH), height, canopy, tree density, tree composition, origin, forest management type, volume, forest resource distribution map, etc. The study area includes 20,452 and 53,980 irregular forest patches with the average of 6.499 ± 15.847 (mean ± SD) ha and 2.376 ± 6.784 ha in 2010 and 2017, respectively. The data selected for this research are (1) forest attribute data (age, DBH, dominant species, and species composition), (2) soil data (layer thickness and type), and (3) topographic data (altitude, slope degree, and slope direction). The values of the above data were calculated by averaging values of all units in the forest patches. The tree species of Siming Mountain National Forest Park mainly includePinus massoniana, Cunninghamia lanceolata, Pinus taiwanensis, Pinus thunbergii, Pinus elliottii, Schima superba, Liquidambar formosana, Phyllostachys heterocycla, Dendrocalamopsis oldhami, Myrica rubra, Eriobotrya japonica, Camellia sinensis, Morus alba, and Camellia oleifera. According to the dominant species and tree species composition data recorded in each patch, forest types were divided into six categories: coniferous tree, broad-leaved tree, mixed coniferous and broad-leaved, bamboo, shrub, and nonforest land (Figure 2).Figure 2
Forest types of the Siming Mountain National Forest Park in (a) 2010 and (b) 2017.
(a)
(b)The third type of data was obtained from open-access resources. It consisted of daily Moderate Resolution Imaging Spectroradiometer (MODIS), total precipitable water vapor (TPWV), daily leaf area index (LAI), daily land surface temperature (LST), 16-day normalised difference vegetation index (NDVI), MODIS Terra Aerosol 5-Min L2 Swath 3 km (MOD04_3K) and Landscan population (LSP) data (https://landscan.ornl.gov/) [30], and OpenStreetMap (OSM) data (https://www.openstreetmap.org) [31]. All data were selected in accordance with the shooting time of the Landsat satellite imagery dataset. Due to cloud problems, some missing pixel values of the MODIS data (TPWV, LAI, LST, and NDVI) were filled by Kriging interpolation [32–35]. Then, the spatial resolution of the whole dataset was resampled to 30 m [36].


### 2.4. Inversion of the AOD

The atmospheric top radiation captured by satellite sensors is the result of the interaction of electromagnetic waves and the earth’s atmospheric system. Atmospheric aerosol remote sensing, which is sensitive to aerosol scattering, depends on the characteristics of the short wavelength of visible light combined with Second Simulation of the Satellite Signal in the Solar Spectrum, 6 S atmosphere transmission model, to realise the inversion [37–39]. Assuming the earth surface to be a Lambert surface, the atmospheric level is uniform, and the apparent atmospheric reflectance (upward reflectance) observed by the satellite is as described in the following equation:(1)ρTOAθs,θv,φ=ρ0θs,θv,φ+TθsTθvρsθs,θv,φ1−ρsθs,θv,φ⋅R,where θs is the solar zenith angle, θv is the satellite zenith angle, φ is the azimuth of the scattered radiation from the solar beam, ρ0 is the path radiance, ρs is the angular surface reflectance, R is the atmospheric backscattering ratio, and Tθs is the normalised downward flux for zero surface reflectance; Tθv represents upward total transmission into the satellite’s field of view.In this study, we used ENVI software and IDL language to calculate AOD values at 550 nm. The first step was preprocessing of the Landsat images, including masking, radiation calibration, geometric correction, and the calculation of the apparent reflectance (ρ_TOA) of the atmosphere. The second step was inversion of the AOD by using the dark pixel method because of the high vegetation cover of the study area. The detailed processes include four substeps [39, 40]: (1) the dark pixels, which are usually in the dense vegetation covered areas with short wavelength and dark surface, were identified by the NDVI; (2) since the apparent reflectance of 2.1-2.2 μm wavelength is almost independent of aerosol [41–43], the short infrared band was chosen to calculate the surface reflectance; (3) the optimal AOD was determined by matching the 6 S LookUp Table (LUT) and apparent reflectance and the mean AOD mean values of red and blue band were obtained; (4) the Kriging interpolation was applied to the preliminary inversion results which were resampled to 30 m [44]; and (5)the AOD inversion results in this study were compared with MOD04_3K data to verify the inversion accuracy. In our study, For this substep, we firstly registered Landsat AOD and MOD04_3K data which were taken at the same time (two sets of data at the same time), then calculated the average value of the corresponding mesh of Landsat AOD according to the mesh size of MOD04_3K, counted the number of effective pixels, and finally counted the linear regression equation and R2.


### 2.5. Statistical Method

#### 2.5.1. Spatial Statistical Analysis

Two mean AOD spatial distribution maps were mapped based on the calculation of the average of each pixel of four AOD images in 2010 and 2017, respectively. The Getis-Ord Gi∗ method effectively identified the statistically significant aggregation areas (e.g., hot spots representing the position with the high AOD concentrations (p<0.05) and cold spots representing low AOD concentrations (p<0.05)) [40], which is popular in spatial analysis [41, 42]. Moreover, we evaluated the AOD autocorrelation pattern by the Moran I index. The Getis-Ord Gi∗ demands optimal threshold distances for clustering. Therefore, we used the incremental spatial autocorrelation module to determine the optimal distance threshold, which is required, by increasing the threshold distance from 200 to 4,000 m at the interval of 200 m until it reached its maximum value (2,200 m) [43]. The distance was calculated using the Euclidean distance.


#### 2.5.2. Geodetector Model

The Geodetector is a statistical tool for detecting spatially stratified heterogeneity and revealing the contribution of factors to the heterogeneity [44]. In this study, we used the factor detector (FD) module of this model to assess the explanation powers of the park features (e.g., tree species and NDVI), the topography (e.g., altitude, slope degree, and slope aspect), the soil (e.g., soil type and thickness), the human activity (e.g., road length and population density), and the meteorology (e.g., temperature and evaporation) on the AOD. The contribution power can be expressed using the following equation:(2)q=1−∑h=1LNhσh2Nσ2=1−SSWSST,where q ∈ (0,1) is the impact power of the independent variable on the dependent variable 1, 2, 3, …, L, L is a strata of dependent or independent variables; stratum h is composed of Nh unit; σh2 is the variance of AOD in the same stratum; N is the number of FMPI patches; σ2 is the variance of the AOD in the entire region; N is the number of FMPI patches in the same stratum; SSW is the within sum of squares; and SST is the total sum of squares. The larger the q value is, the greater the similarity between independent factors and AOD.


#### 2.5.3. Statistical Analysis

The independent sampleT test was used to analyze the difference of biomass growth rate of different forest types in cold and hot spots from 2010 to 2017. The system clustering was used to classify the differences of independent variables.To ensure that the analysis of biomass growth rate is not affected by seedling plantation or forest fire, we selected forest patches with biomass growth from 2010 to 2017. Moreover, we also screened the forest patches with the same forest type in both years in order to reduce statistical errors.


### 2.6. Biomass Estimation

The biomass is an important indicator to evaluate the quality of forest park because it is a fundamental factor in evaluating plant growth. The Volume-Biomass Compatibility Model (VBCM) was used to estimate individual tree biomass by tree height, the DBH, and the coefficients [45–48]. The forest biomass density was inferred from the number of trees and the area of forest patches recorded by the FMPI. Although there were errors in estimating regional biomass by inference, it would not significantly influence the analysis result because the growth rate was the ratio of the regional biomass estimates in the two years. The formula is shown in the following equation:(3)W=aDkHj×numtreeA,where W is the biomass density (t/ha); a, k, and j are the coefficients; D is the DBH (cm); H is the tree height (m); numtree is the number of patch plants; and A is the area (ha). Dominant tree species and forest composition were recorded in the FMPI using the coefficients of the dominant tree species (Cunninghamia lanceolate, Phyllostachys heterocycle, Schima superba, and Camellia oleifera) [47, 49, 50] in Table 2 to calculate the biomass.Table 2
Coefficients of different dominant tree species in the VBCM.
Tree species (dominant species)
a
k
j
Coniferous
Cunninghamia lanceolata
0.0811
1.6942
0.8472
Bamboo
Phyllostachys heterocycla
66.9197
2.5500
0.0437
Broad-leaved
Schima superba
0.2993
1.8530
0.2774
Shrub
Camellia oleifera
0.1510
2.0170
0.0000


## 2.1. Overview

We approached this research in the following steps: First, to create a database that integrates multiple types of GIS and remote sensing (RS) data (including Landsat 8 OLI images, Landsat 7-ETM+ images, FMPI, MODIS data, OpenStreetMap data, and Landscan data). Second, to preprocess the database. There were two substeps included: (1) replenishment of missing remote sensing pixel values, caused by cloud occlusion, through the Kriging interpolation which is a geostatistical method for predicting missing spatial information (AOD, MODIS data) [26, 27] and (2) inversion of the AOD by Landsat data. Third, the hot spots and cold spots of AOD in this area were identified by the high Z values (Z ≥ 1.64) and low Z values (Z ≤ −1.64) respectively, which were calculated through global Moran’s I and Getis-Ord Gi∗ [28, 29]. The optimal threshold distance was calculated using the incremental spatial autocorrelation module. Fourth, the factor detector of the Geodetector model was used to quantify the impact contribution of different factors on the AOD. Lastly, we calculated and analyzed the biomass changes of different forest types in cold and hot spots from 2010 to 2017.


## 2.2. Study Area

Siming Mountain National Forest Park is in Ningbo City, Zhejiang Province, China (Figure1). The location is longitude 120°59′20″E–121°25′16″E and latitude 29°31′39″N–29°59′20″N. The park is subject to the subtropical monsoon climate, which is characterised by 1,277 km2, 976 m highest altitude, 77.5% forest coverage, and yellow soil. The annual average temperature is 10ºC, and the annual average rainfall is 1,800 mm.Figure 1
Map of study area. The inset shows the location in China.


## 2.3. Multisource Dataset and Preprocessing

Three types of data in two years, which included the year 2010 and 2017, were collected. The first type of data was the Landsat satellite remote sensing data for AOD retrieval with a spatial resolution of 30 m. To ensure representation of the data, we selected the two nearest meteorological stations in the China Meteorological Data Network (http://data.cma.cn) (Nos. 58467 and 58562, respectively) to screen for wind speed and air pressure that were within a 5% error limitation of the mean value in both 2010 and 2017. Afterwards, we collected eight remote sensing image data for the four seasons in 2010 and 2017 (row number 118, column number 39). The remote sensing images and data from the four phases of the ETM sensor and the four phases of the OLI sensor were chosen based on the cloud cover percentage. The less cloudy ones were selected, which were the Landsat 7 ETM images on January 1, April 23, July 28, and December 3, 2010, and Landsat 8 OLI image on April 2, August 24, and October 27, 2017, and January 15, 2018 (Table 1). The cloud covers were 0.01%, 6.62%, 24.63%, 0.10%, 0.37%, 0.26%, 14.27%, and 1.27%, respectively. The average wind speeds at the two meteorological stations were (2.1 and 3.2 m/s), (4.1 and 3.1 m/s), (2.1 and 4.7 m/s), (3.6 and 3.5 m/s), (4.2 and 4.1 m/s), (2.2 and 2.3 m/s), (5.1 and 6.1 m/s), and (5.5 and 5.5 m/s) with relative humidity (75 and 77%), (88 and 86%), (90 and 91%), (88 and 85%), (79 and 82%), (85 and 90%), (90 and 90%), and (77 and 78%) in the years 2010 and 2017, correspondingly.Table 1
Information about the GIS and RS data in this study.
Type
Data
Time
GIS
FMPI
2010
2016
OSM
2017
LSP
2017
RS
Landsat 8 OLIRow 118/path 39
April 2, 2017
August 24, 2017
October 27, 2017
January 15, 2018
RS
Landsat 7 ETMRow 118/path 39
January 1, 2010
April 23, 2010
July 28, 2010
December 3, 2010
MODIS TPWV
The time of MODIS data set is the same as that of Landsat data set
MODIS LAI
MODIS LST
MODIS NDVI
MOD04_3KThe second type of data was the FMPI, which is obtained from in situ investigation and observation by the national forestry bureau and its affiliations in China. These data provide forest ground survey results on a regional scale. The attribute database contains information on forest patch area, diameter of breast height (DBH), height, canopy, tree density, tree composition, origin, forest management type, volume, forest resource distribution map, etc. The study area includes 20,452 and 53,980 irregular forest patches with the average of 6.499 ± 15.847 (mean ± SD) ha and 2.376 ± 6.784 ha in 2010 and 2017, respectively. The data selected for this research are (1) forest attribute data (age, DBH, dominant species, and species composition), (2) soil data (layer thickness and type), and (3) topographic data (altitude, slope degree, and slope direction). The values of the above data were calculated by averaging values of all units in the forest patches. The tree species of Siming Mountain National Forest Park mainly includePinus massoniana, Cunninghamia lanceolata, Pinus taiwanensis, Pinus thunbergii, Pinus elliottii, Schima superba, Liquidambar formosana, Phyllostachys heterocycla, Dendrocalamopsis oldhami, Myrica rubra, Eriobotrya japonica, Camellia sinensis, Morus alba, and Camellia oleifera. According to the dominant species and tree species composition data recorded in each patch, forest types were divided into six categories: coniferous tree, broad-leaved tree, mixed coniferous and broad-leaved, bamboo, shrub, and nonforest land (Figure 2).Figure 2
Forest types of the Siming Mountain National Forest Park in (a) 2010 and (b) 2017.
(a)
(b)The third type of data was obtained from open-access resources. It consisted of daily Moderate Resolution Imaging Spectroradiometer (MODIS), total precipitable water vapor (TPWV), daily leaf area index (LAI), daily land surface temperature (LST), 16-day normalised difference vegetation index (NDVI), MODIS Terra Aerosol 5-Min L2 Swath 3 km (MOD04_3K) and Landscan population (LSP) data (https://landscan.ornl.gov/) [30], and OpenStreetMap (OSM) data (https://www.openstreetmap.org) [31]. All data were selected in accordance with the shooting time of the Landsat satellite imagery dataset. Due to cloud problems, some missing pixel values of the MODIS data (TPWV, LAI, LST, and NDVI) were filled by Kriging interpolation [32–35]. Then, the spatial resolution of the whole dataset was resampled to 30 m [36].


## 2.4. Inversion of the AOD

The atmospheric top radiation captured by satellite sensors is the result of the interaction of electromagnetic waves and the earth’s atmospheric system. Atmospheric aerosol remote sensing, which is sensitive to aerosol scattering, depends on the characteristics of the short wavelength of visible light combined with Second Simulation of the Satellite Signal in the Solar Spectrum, 6 S atmosphere transmission model, to realise the inversion [37–39]. Assuming the earth surface to be a Lambert surface, the atmospheric level is uniform, and the apparent atmospheric reflectance (upward reflectance) observed by the satellite is as described in the following equation:(1)ρTOAθs,θv,φ=ρ0θs,θv,φ+TθsTθvρsθs,θv,φ1−ρsθs,θv,φ⋅R,where θs is the solar zenith angle, θv is the satellite zenith angle, φ is the azimuth of the scattered radiation from the solar beam, ρ0 is the path radiance, ρs is the angular surface reflectance, R is the atmospheric backscattering ratio, and Tθs is the normalised downward flux for zero surface reflectance; Tθv represents upward total transmission into the satellite’s field of view.In this study, we used ENVI software and IDL language to calculate AOD values at 550 nm. The first step was preprocessing of the Landsat images, including masking, radiation calibration, geometric correction, and the calculation of the apparent reflectance (ρ_TOA) of the atmosphere. The second step was inversion of the AOD by using the dark pixel method because of the high vegetation cover of the study area. The detailed processes include four substeps [39, 40]: (1) the dark pixels, which are usually in the dense vegetation covered areas with short wavelength and dark surface, were identified by the NDVI; (2) since the apparent reflectance of 2.1-2.2 μm wavelength is almost independent of aerosol [41–43], the short infrared band was chosen to calculate the surface reflectance; (3) the optimal AOD was determined by matching the 6 S LookUp Table (LUT) and apparent reflectance and the mean AOD mean values of red and blue band were obtained; (4) the Kriging interpolation was applied to the preliminary inversion results which were resampled to 30 m [44]; and (5)the AOD inversion results in this study were compared with MOD04_3K data to verify the inversion accuracy. In our study, For this substep, we firstly registered Landsat AOD and MOD04_3K data which were taken at the same time (two sets of data at the same time), then calculated the average value of the corresponding mesh of Landsat AOD according to the mesh size of MOD04_3K, counted the number of effective pixels, and finally counted the linear regression equation and R2.


## 2.5. Statistical Method

### 2.5.1. Spatial Statistical Analysis

Two mean AOD spatial distribution maps were mapped based on the calculation of the average of each pixel of four AOD images in 2010 and 2017, respectively. The Getis-Ord Gi∗ method effectively identified the statistically significant aggregation areas (e.g., hot spots representing the position with the high AOD concentrations (p<0.05) and cold spots representing low AOD concentrations (p<0.05)) [40], which is popular in spatial analysis [41, 42]. Moreover, we evaluated the AOD autocorrelation pattern by the Moran I index. The Getis-Ord Gi∗ demands optimal threshold distances for clustering. Therefore, we used the incremental spatial autocorrelation module to determine the optimal distance threshold, which is required, by increasing the threshold distance from 200 to 4,000 m at the interval of 200 m until it reached its maximum value (2,200 m) [43]. The distance was calculated using the Euclidean distance.


### 2.5.2. Geodetector Model

The Geodetector is a statistical tool for detecting spatially stratified heterogeneity and revealing the contribution of factors to the heterogeneity [44]. In this study, we used the factor detector (FD) module of this model to assess the explanation powers of the park features (e.g., tree species and NDVI), the topography (e.g., altitude, slope degree, and slope aspect), the soil (e.g., soil type and thickness), the human activity (e.g., road length and population density), and the meteorology (e.g., temperature and evaporation) on the AOD. The contribution power can be expressed using the following equation:(2)q=1−∑h=1LNhσh2Nσ2=1−SSWSST,where q ∈ (0,1) is the impact power of the independent variable on the dependent variable 1, 2, 3, …, L, L is a strata of dependent or independent variables; stratum h is composed of Nh unit; σh2 is the variance of AOD in the same stratum; N is the number of FMPI patches; σ2 is the variance of the AOD in the entire region; N is the number of FMPI patches in the same stratum; SSW is the within sum of squares; and SST is the total sum of squares. The larger the q value is, the greater the similarity between independent factors and AOD.


### 2.5.3. Statistical Analysis

The independent sampleT test was used to analyze the difference of biomass growth rate of different forest types in cold and hot spots from 2010 to 2017. The system clustering was used to classify the differences of independent variables.To ensure that the analysis of biomass growth rate is not affected by seedling plantation or forest fire, we selected forest patches with biomass growth from 2010 to 2017. Moreover, we also screened the forest patches with the same forest type in both years in order to reduce statistical errors.


## 2.5.1. Spatial Statistical Analysis

Two mean AOD spatial distribution maps were mapped based on the calculation of the average of each pixel of four AOD images in 2010 and 2017, respectively. The Getis-Ord Gi∗ method effectively identified the statistically significant aggregation areas (e.g., hot spots representing the position with the high AOD concentrations (p<0.05) and cold spots representing low AOD concentrations (p<0.05)) [40], which is popular in spatial analysis [41, 42]. Moreover, we evaluated the AOD autocorrelation pattern by the Moran I index. The Getis-Ord Gi∗ demands optimal threshold distances for clustering. Therefore, we used the incremental spatial autocorrelation module to determine the optimal distance threshold, which is required, by increasing the threshold distance from 200 to 4,000 m at the interval of 200 m until it reached its maximum value (2,200 m) [43]. The distance was calculated using the Euclidean distance.


## 2.5.2. Geodetector Model

The Geodetector is a statistical tool for detecting spatially stratified heterogeneity and revealing the contribution of factors to the heterogeneity [44]. In this study, we used the factor detector (FD) module of this model to assess the explanation powers of the park features (e.g., tree species and NDVI), the topography (e.g., altitude, slope degree, and slope aspect), the soil (e.g., soil type and thickness), the human activity (e.g., road length and population density), and the meteorology (e.g., temperature and evaporation) on the AOD. The contribution power can be expressed using the following equation:(2)q=1−∑h=1LNhσh2Nσ2=1−SSWSST,where q ∈ (0,1) is the impact power of the independent variable on the dependent variable 1, 2, 3, …, L, L is a strata of dependent or independent variables; stratum h is composed of Nh unit; σh2 is the variance of AOD in the same stratum; N is the number of FMPI patches; σ2 is the variance of the AOD in the entire region; N is the number of FMPI patches in the same stratum; SSW is the within sum of squares; and SST is the total sum of squares. The larger the q value is, the greater the similarity between independent factors and AOD.


## 2.5.3. Statistical Analysis

The independent sampleT test was used to analyze the difference of biomass growth rate of different forest types in cold and hot spots from 2010 to 2017. The system clustering was used to classify the differences of independent variables.To ensure that the analysis of biomass growth rate is not affected by seedling plantation or forest fire, we selected forest patches with biomass growth from 2010 to 2017. Moreover, we also screened the forest patches with the same forest type in both years in order to reduce statistical errors.


## 2.6. Biomass Estimation

The biomass is an important indicator to evaluate the quality of forest park because it is a fundamental factor in evaluating plant growth. The Volume-Biomass Compatibility Model (VBCM) was used to estimate individual tree biomass by tree height, the DBH, and the coefficients [45–48]. The forest biomass density was inferred from the number of trees and the area of forest patches recorded by the FMPI. Although there were errors in estimating regional biomass by inference, it would not significantly influence the analysis result because the growth rate was the ratio of the regional biomass estimates in the two years. The formula is shown in the following equation:(3)W=aDkHj×numtreeA,where W is the biomass density (t/ha); a, k, and j are the coefficients; D is the DBH (cm); H is the tree height (m); numtree is the number of patch plants; and A is the area (ha). Dominant tree species and forest composition were recorded in the FMPI using the coefficients of the dominant tree species (Cunninghamia lanceolate, Phyllostachys heterocycle, Schima superba, and Camellia oleifera) [47, 49, 50] in Table 2 to calculate the biomass.Table 2
Coefficients of different dominant tree species in the VBCM.
Tree species (dominant species)
a
k
j
Coniferous
Cunninghamia lanceolata
0.0811
1.6942
0.8472
Bamboo
Phyllostachys heterocycla
66.9197
2.5500
0.0437
Broad-leaved
Schima superba
0.2993
1.8530
0.2774
Shrub
Camellia oleifera
0.1510
2.0170
0.0000


## 3. Results

### 3.1. AOD Distribution and Various Forest Types Spatial Statistical Result in Siming Mountain National Forest Park

As shown in Figure3, there were 717 and 718 verification points in the 2010 and 2017 images of Landsat AOD paired with MOD04_3K, respectively, and R2 was 0.752 and 0.724. These outcomes proved that the inversion results of Landsat AOD are reliable.Figure 3
Verification results of Landsat AOD and MOD04_3K data. (a) 2010. (b) 2017. The color bars represent the counts of points. The red solid line represents the regression line.
(a)
(b)From 2010 to 2017, the annual average AOD decreased from 0.567 to 0.292 with a change in tree species distribution. Specifically, from the proportion of tree species distribution aspect, the area proportional to broad-leaved forest increased greatly (from 24.49% to 42.58%). The area proportion of nonforest land, coniferous forest, mixed coniferous, and broad-leaved forest decreased significantly (from 21.38%, 15.57%, and 9.06% to 16.2%, 6.29%, and 4.14%, respectively), while bamboo (from 23.01% to 23.90%) and shrub (from 6.49 to 6.89%) showed a slight increase. The aggregation degree of AOD distribution decreased in terms of Moran’s I, which was 0.348 and 0.177 in the two years, respectively (Table3).Table 3
Results of Moran’s I.
Year
2010
2017
Moran’s I
0.348∗∗
0.177∗∗
Z-score
98.157
52.983
Pattern
Aggregated
Aggregated
Note.∗∗The significant value (p<0.01).The average AOD of the cold and hot spots area in 2010 was significantly different from that in 2017 (Figure4).The cold spot area increased significantly from 18,119.64 ha to 41,924.62 ha with a decrease in the average AOD from 0.39 to 0.16. While the hot spot area decreased from 57,964.31 ha to 26,706.11 ha with a decrease in the average AOD from 0.69 to 0.52. From the five regions marked in Figure 4, we found the five regions with great changes, including the 1st (north), 3rd (west), and 5th (east) regions where the cold spots areas increased significantly. In 2017, the 2nd and the 4th regions where the hot spots areas increased as well, but with a lower AOD than in 2010.Figure 4
Cold and hot spots of the AOD spatial distribution maps in (a) 2010 and (b) 2017. The regions with 1–5 indicate great changes in hot and cold spots.
(a)
(b)The 10,178 forest patches of cold spots and the 10,581 forest patches of hot spots were selected based on the screening principle in the overlapped regions of cold and hot spots in 2010 and 2017. The area proportion of the various forest types in cold and hot spots area varied dramatically (Table4). The broad-leaved forests and bamboo accounted for the largest areas in the cold spots. Specifically, the percentages of broad-leaved forests and bamboo were 18.28% and 15.15% larger than in the hot spots. In the hot spots, the nonforest land was larger than in the cold spots. Moreover, the proportion of the coniferous forest, shrub, mixed forest, and nonforest land was 5.19%, 7.30%, 8.60%, and 13.15% higher than in cold spots, respectively.Table 4
Area proportion of different forest types in the overlapped cold and hot spots in 2010 and 2017.
%
Coniferous forest
Broad-leaved forest
Shrub
Nonforest
Bamboo
Con and bro mixed forest
Total
Cold spots
3.80
50.50
3.00
1.80
39.80
1.10
100.00
Hot spots
9.09
32.22
10.30
14.95
24.75
9.70
100.00
Con and bro mixed forest = coniferous and broad-leaved mixed forest.


### 3.2. 2010 and 2017 Factor Detector Results of Geodetector Model

The forest type played the most important role in both cold and hot spots during the two years. In the year of 2010, temperature (0.0301 and 0.0159) and NDVI (0.0257 and 0.0354) ranked second and third in both cold and hot spots. In the year of 2017, the second and third important factors were the Landscan population density (0.0189) and the OSM road length (0.0195) in hot spots and the NDVI (0.0204) and Landscan population density (0.0201) in cold spots. In the two-year average, the followers were the OSM road length (0.0234) and NDVI (0.0158) in the hot spots and the NDVI (0.0273) and slope (0.0064) in the cold spots. The forest type and other factors were classified into two clusters through the hierarchical cluster analysis (Table5).Table 5
The results of different factors through FD and cluster analysis in 2010, 2017 and the two periods’ average.
Factors
2010
2017
Two periods’ average
Hierarchical cluster
Hot spots
Cold spots
Hot spots
Cold spots
Hot spots
Cold spots
TSY
0.1248∗∗
0.1196∗∗
0.1389∗∗
0.1450∗∗
0.1569∗∗
0.1735∗∗
1
NDVI
0.0257∗∗
0.0354∗∗
0.0146∗∗
0.0204∗∗
0.0158∗∗
0.0273∗∗
2
LST
0.0301∗∗
0.0159∗∗
0.0101∗∗
0.0095∗∗
0.0094∗∗
0.0057∗∗
2
TPWV
0.0009∗∗
0.0008∗∗
0.0011∗∗
0.0012∗∗
0.0005∗∗
0.0027∗∗
2
ST
0.0099∗∗
0.0121∗∗
0.0115∗∗
0.0141∗∗
0.0093∗∗
0.0051∗∗
2
SLT
0.0085∗∗
0.0099∗∗
0.0054∗∗
0.0021∗∗
0.0056∗∗
0.0027∗∗
2
AT
0.0114∗∗
0.0089∗∗
0.0021∗∗
0.0012∗∗
0.0078∗∗
0.007∗∗
2
SL
0.0088∗∗
0.0066∗∗
0.0025∗∗
0.0045∗∗
0.0098∗∗
0.0058∗∗
2
SD
0.0027∗∗
0.0056∗∗
0.0077∗∗
0.0096∗∗
0.0025∗∗
0.0064∗∗
2
LSP
0.0012∗∗
0.0098∗∗
0.0189∗∗
0.0201∗∗
0.0029∗∗
0.0003∗∗
2
OSM
0.0012∗∗
0.0083∗∗
0.0195∗∗
0.0055∗∗
0.0234∗∗
0.0017∗∗
2
Note.p values for all areas are <0.01 in 2010, 2017, and the two-year average. ∗∗High significance (p<0.001). TSY = tree species type; ST = soil type; SLT =  soil layer thickness; AT =  altitude; SL = slope; SD = slope direction; LSP = Landscan population; OSM = OSM road length.


### 3.3. Biomass Variations and Growth Rates of Various Forest Types in 2010 and 2017

Analysis of samples from the broad-leaved forest, shrub, coniferous and broad-leaved mixed forest, coniferous forest, and bamboo were in the number of 385, 214, 31, 196, and 63 patches in the hot spots and 596, 88, 22, 191, and 239 patches in the cold spots, respectively.In 2010, there was no significant difference between the biomass densities of different forest types in the hot and cold spots. The average biomass densities of the broad-leaved forest, mixed forest, coniferous forest, bamboo, and shrub in the hot spots were 142.82 t/ha, 126.57 t/ha, 102.11 t/ha, 51.82 t/ha, and 38.25 t/ha, and in the cold spots, 143.18 t/ha, 128.14 t/ha, 100.68 t/ha, 51.53 t/ha, and 36.26 t/ha, respectively.However, in 2017, the biomass density differences were significant between the hot and cold spots (p<0.01) for all of the forest types except for the mixed forest. The biomass densities of these five forest types were 187.09 t/ha, 158.21 t/ha, 152.14 t/ha, 80.8 3 t/ha, and 49.34 t/ha for the broad-leaved forest, mixed forest, coniferous forest, bamboo, and shrub in the hot spots; 181.84 t/ha, 161.45 t/ha, 141.96 t/ha, 75.23 t/ha, and 43.87 t/ha of them were in the cold spots, respectively.Comparing the growth rate of biomass density in different forest types from 2010 to 2017 (Figure5), the order of average biomass growth rate was bamboo > coniferous forest > broad-leaved forest > shrub > mixed forest in hot spots. In cold spots, the order of average biomass growth rate was bamboo > coniferous forest > mixed forest > broad-leaved forest > shrub. The average ages of the broad-leaved forest, the mixed forest, the coniferous forest, the bamboo, and the shrub increased from 15.42 ± 2.16 to 15.84 ± 2.31 (mean ± SD), 14.88 ± 2.55 to 15.27 ± 2.78, 9.01 ± 1.25 to 9.25 ± 1.12, 3.18 ± 0.97 to 3.51 ± 0.98, and 2.11 ± 0.88 to 2.44 ± 0.85, respectively. Except for the mixed forest, the growth rate of biomass densities of the four other forest types were significantly higher in hot spots than in cold spots (p<0.01). The growth rate of biomass density of the mixed forest in cold spots was higher than that in hot spots, but not significant.Figure 5
Biomass and biomass growth rates of different forest types in 2010 and 2017 (age and SD in 2010). The different letters indicate significant differences in the biomass increase rate between forest types atp<0.05.


## 3.1. AOD Distribution and Various Forest Types Spatial Statistical Result in Siming Mountain National Forest Park

As shown in Figure3, there were 717 and 718 verification points in the 2010 and 2017 images of Landsat AOD paired with MOD04_3K, respectively, and R2 was 0.752 and 0.724. These outcomes proved that the inversion results of Landsat AOD are reliable.Figure 3
Verification results of Landsat AOD and MOD04_3K data. (a) 2010. (b) 2017. The color bars represent the counts of points. The red solid line represents the regression line.
(a)
(b)From 2010 to 2017, the annual average AOD decreased from 0.567 to 0.292 with a change in tree species distribution. Specifically, from the proportion of tree species distribution aspect, the area proportional to broad-leaved forest increased greatly (from 24.49% to 42.58%). The area proportion of nonforest land, coniferous forest, mixed coniferous, and broad-leaved forest decreased significantly (from 21.38%, 15.57%, and 9.06% to 16.2%, 6.29%, and 4.14%, respectively), while bamboo (from 23.01% to 23.90%) and shrub (from 6.49 to 6.89%) showed a slight increase. The aggregation degree of AOD distribution decreased in terms of Moran’s I, which was 0.348 and 0.177 in the two years, respectively (Table3).Table 3
Results of Moran’s I.
Year
2010
2017
Moran’s I
0.348∗∗
0.177∗∗
Z-score
98.157
52.983
Pattern
Aggregated
Aggregated
Note.∗∗The significant value (p<0.01).The average AOD of the cold and hot spots area in 2010 was significantly different from that in 2017 (Figure4).The cold spot area increased significantly from 18,119.64 ha to 41,924.62 ha with a decrease in the average AOD from 0.39 to 0.16. While the hot spot area decreased from 57,964.31 ha to 26,706.11 ha with a decrease in the average AOD from 0.69 to 0.52. From the five regions marked in Figure 4, we found the five regions with great changes, including the 1st (north), 3rd (west), and 5th (east) regions where the cold spots areas increased significantly. In 2017, the 2nd and the 4th regions where the hot spots areas increased as well, but with a lower AOD than in 2010.Figure 4
Cold and hot spots of the AOD spatial distribution maps in (a) 2010 and (b) 2017. The regions with 1–5 indicate great changes in hot and cold spots.
(a)
(b)The 10,178 forest patches of cold spots and the 10,581 forest patches of hot spots were selected based on the screening principle in the overlapped regions of cold and hot spots in 2010 and 2017. The area proportion of the various forest types in cold and hot spots area varied dramatically (Table4). The broad-leaved forests and bamboo accounted for the largest areas in the cold spots. Specifically, the percentages of broad-leaved forests and bamboo were 18.28% and 15.15% larger than in the hot spots. In the hot spots, the nonforest land was larger than in the cold spots. Moreover, the proportion of the coniferous forest, shrub, mixed forest, and nonforest land was 5.19%, 7.30%, 8.60%, and 13.15% higher than in cold spots, respectively.Table 4
Area proportion of different forest types in the overlapped cold and hot spots in 2010 and 2017.
%
Coniferous forest
Broad-leaved forest
Shrub
Nonforest
Bamboo
Con and bro mixed forest
Total
Cold spots
3.80
50.50
3.00
1.80
39.80
1.10
100.00
Hot spots
9.09
32.22
10.30
14.95
24.75
9.70
100.00
Con and bro mixed forest = coniferous and broad-leaved mixed forest.


## 3.2. 2010 and 2017 Factor Detector Results of Geodetector Model

The forest type played the most important role in both cold and hot spots during the two years. In the year of 2010, temperature (0.0301 and 0.0159) and NDVI (0.0257 and 0.0354) ranked second and third in both cold and hot spots. In the year of 2017, the second and third important factors were the Landscan population density (0.0189) and the OSM road length (0.0195) in hot spots and the NDVI (0.0204) and Landscan population density (0.0201) in cold spots. In the two-year average, the followers were the OSM road length (0.0234) and NDVI (0.0158) in the hot spots and the NDVI (0.0273) and slope (0.0064) in the cold spots. The forest type and other factors were classified into two clusters through the hierarchical cluster analysis (Table5).Table 5
The results of different factors through FD and cluster analysis in 2010, 2017 and the two periods’ average.
Factors
2010
2017
Two periods’ average
Hierarchical cluster
Hot spots
Cold spots
Hot spots
Cold spots
Hot spots
Cold spots
TSY
0.1248∗∗
0.1196∗∗
0.1389∗∗
0.1450∗∗
0.1569∗∗
0.1735∗∗
1
NDVI
0.0257∗∗
0.0354∗∗
0.0146∗∗
0.0204∗∗
0.0158∗∗
0.0273∗∗
2
LST
0.0301∗∗
0.0159∗∗
0.0101∗∗
0.0095∗∗
0.0094∗∗
0.0057∗∗
2
TPWV
0.0009∗∗
0.0008∗∗
0.0011∗∗
0.0012∗∗
0.0005∗∗
0.0027∗∗
2
ST
0.0099∗∗
0.0121∗∗
0.0115∗∗
0.0141∗∗
0.0093∗∗
0.0051∗∗
2
SLT
0.0085∗∗
0.0099∗∗
0.0054∗∗
0.0021∗∗
0.0056∗∗
0.0027∗∗
2
AT
0.0114∗∗
0.0089∗∗
0.0021∗∗
0.0012∗∗
0.0078∗∗
0.007∗∗
2
SL
0.0088∗∗
0.0066∗∗
0.0025∗∗
0.0045∗∗
0.0098∗∗
0.0058∗∗
2
SD
0.0027∗∗
0.0056∗∗
0.0077∗∗
0.0096∗∗
0.0025∗∗
0.0064∗∗
2
LSP
0.0012∗∗
0.0098∗∗
0.0189∗∗
0.0201∗∗
0.0029∗∗
0.0003∗∗
2
OSM
0.0012∗∗
0.0083∗∗
0.0195∗∗
0.0055∗∗
0.0234∗∗
0.0017∗∗
2
Note.p values for all areas are <0.01 in 2010, 2017, and the two-year average. ∗∗High significance (p<0.001). TSY = tree species type; ST = soil type; SLT =  soil layer thickness; AT =  altitude; SL = slope; SD = slope direction; LSP = Landscan population; OSM = OSM road length.


## 3.3. Biomass Variations and Growth Rates of Various Forest Types in 2010 and 2017

Analysis of samples from the broad-leaved forest, shrub, coniferous and broad-leaved mixed forest, coniferous forest, and bamboo were in the number of 385, 214, 31, 196, and 63 patches in the hot spots and 596, 88, 22, 191, and 239 patches in the cold spots, respectively.In 2010, there was no significant difference between the biomass densities of different forest types in the hot and cold spots. The average biomass densities of the broad-leaved forest, mixed forest, coniferous forest, bamboo, and shrub in the hot spots were 142.82 t/ha, 126.57 t/ha, 102.11 t/ha, 51.82 t/ha, and 38.25 t/ha, and in the cold spots, 143.18 t/ha, 128.14 t/ha, 100.68 t/ha, 51.53 t/ha, and 36.26 t/ha, respectively.However, in 2017, the biomass density differences were significant between the hot and cold spots (p<0.01) for all of the forest types except for the mixed forest. The biomass densities of these five forest types were 187.09 t/ha, 158.21 t/ha, 152.14 t/ha, 80.8 3 t/ha, and 49.34 t/ha for the broad-leaved forest, mixed forest, coniferous forest, bamboo, and shrub in the hot spots; 181.84 t/ha, 161.45 t/ha, 141.96 t/ha, 75.23 t/ha, and 43.87 t/ha of them were in the cold spots, respectively.Comparing the growth rate of biomass density in different forest types from 2010 to 2017 (Figure5), the order of average biomass growth rate was bamboo > coniferous forest > broad-leaved forest > shrub > mixed forest in hot spots. In cold spots, the order of average biomass growth rate was bamboo > coniferous forest > mixed forest > broad-leaved forest > shrub. The average ages of the broad-leaved forest, the mixed forest, the coniferous forest, the bamboo, and the shrub increased from 15.42 ± 2.16 to 15.84 ± 2.31 (mean ± SD), 14.88 ± 2.55 to 15.27 ± 2.78, 9.01 ± 1.25 to 9.25 ± 1.12, 3.18 ± 0.97 to 3.51 ± 0.98, and 2.11 ± 0.88 to 2.44 ± 0.85, respectively. Except for the mixed forest, the growth rate of biomass densities of the four other forest types were significantly higher in hot spots than in cold spots (p<0.01). The growth rate of biomass density of the mixed forest in cold spots was higher than that in hot spots, but not significant.Figure 5
Biomass and biomass growth rates of different forest types in 2010 and 2017 (age and SD in 2010). The different letters indicate significant differences in the biomass increase rate between forest types atp<0.05.


## 4. Discussion

The exploration of the impact factors of aerosol spatial distribution and the influences of AOD on the canopies of different forest types in the exurbs of cities encourage forest managers to improve the ecosystem service level. This study used multisource data, which included high-resolution remote sensing images and detailed ground investigation data to discover if the type of tree species is the one of the main factors affecting the AOD concentration difference in Siming Mountain Forest Park. We observed the impacts of AOD concentration on the growth of different tree species in the forest park, which supplemented the evidence of the relationship between AOD spatial distribution and the forest on a regional scale.


### 4.1. Relationship between Different Forest Types and AOD Distribution in the National Park

Siming Mountain Forest Park is located near the east coast of China, with a high forest coverage rate far from the urban center. The variation of AOD values may be caused by the combined action of moisture and particulate matter. In order to minimize the influence of water vapor, we selected eight remote sensing images taken during good and stable weather conditions between 10 : 00 and 11 : 00 am to avoid the condensation of water vapor caused by low temperature in the morning, although the average value of AOD in the study area is low. The AOD is the integral of the extinction coefficient of the medium in the vertical direction, which reflects the comprehensive effect of atmospheric scattering and surface reflection. We cannot distinguish the aerosol types. Based on the analysis of aerosol sources in four forest areas in Sichuan Province, China, it was found that the aerosols were mainly from local primary emission sources such as crust dust and biomass combustion, followed by the long distance transmission of aged organic aerosols (aged OOA) [51]. In addition to soil dust, meteorology, and human activities, plant species was also an important factor for the AOD variation in forest areas. The reactivity of plant VOCs with oxidants (such as OH, NO3 free radical, and O3) was higher than that of artificial VOCs, resulting in the generation of secondary organic aerosols (SOA). Zhu et al. (2017) found that the main source of SOA at the top of the exurban Tai Mountain in China was from an artificial source (11.6–42.6%), such as the use of solvents, fuel combustion, and vehicle exhaust as well as the biological source (6.5%) such as secondary oxidation and biomass burning [52]. However, the source of a large proportion of SOA was not identified. Although plants contribute to the generation of aerosols, the underlying surface of forest cover affects the physical processes in the atmospheric boundary layer. The high vegetation coverage increases the friction between the air and the ground, reduces wind speed, and diminishes transmission of pollutants from the ground and over long distances.The results of the Geodetector and cluster analysis in the cold and hot spots of this study in two years showed that the forest type was the main factor affecting AOD distribution in the exurban forest park, followed by NDVI and temperature. In the cold spots, more than half of the area was broad-leaved forest, while the percentage of broad-leaved forest decreased and the proportion of nonforest land area augmented in the hot spots. The result is similar to that of Deng et al. [53]. The maximized dust retention ability of the coniferous leaf for a single tree is higher than the broad-leaved tree species. However, there are many factors that influence the total dust retention quantity of the forest. As to the specific forest type, in addition to the blade characteristics, the factors influencing tree growth, such as temperature, humidity, soil nutrient, precipitation, and illumination, may affect the stand of leaf area index, thus further modifying the function of particulate matter retention. Moreover, the microstructure of the leaf surface is conducive to the deposition of particulate matter [54]. The forest increases the intensity of end flow in the surface atmosphere, thereby accelerating the deposition rate of particulate matter [55, 56]. Meanwhile, vegetation effectively absorbs pollutants such as SO2, NOx, and O3 in the air in favor of reducing the formation of secondary aerosols [51, 57]. In addition, the large amount of the particle matter in the air can be trapped on the leaf surface as compared to the small amount trapped in the leaf epidermal wax layer [58, 59]. The delayed particulate matter on the leaf surface is eluted to the ground with the precipitation leading to the recovery of its retention ability [60].From the perspective of the impact of aerosols on plants, the concentration threshold reversed the impact of aerosols on photosynthesis above which the net radiation (direct + scattering) effect of atmospheric aerosols would reduce the photosynthetic absorption of plants [61]. Xu quantified the AOD threshold of the impact of aerosol net radiation effect on ecosystem NPP in different regions of China and found that when the AOD was less than 1.5 at 550 nm in the southeast region, the aerosol net radiation effect increased NPP [16]. In this study, the average AOD value in hot spots was much lower than the threshold, and the results were in line with expectations. For example, the growth rate of tree species in hot spots was generally faster than that in cold spots. The average growth rates of the coniferous forest, broad-leaved forest, shrubs, and bamboos in cold spots were significantly different than in hot spots. The reason may be the diffusion radiation fertilization effect mentioned in the introduction, which increases the NPP and promotes growth in hot spots. However, there was not a significant difference in the growth rate of the coniferous broad-leaved mixed forest. The reason may be the large age variation between different tree species (large SD), which induced the different effects of diffusion radiation fertilization on tree species of different ages. The results of this study were consistent with other studies and verified the diffusion radiation fertilization effect at the regional forest park scale.


### 4.2. Limits and Prospects

There are two limitations in this study. First, the change of AOD is continuous with time. Landsat remote sensing satellite image data (a scene every 16 days) only obtains instantaneous information of the environment. Although the average value of four specific periods throughout the year was used as the research object in this study, the representativeness of inversion results still needs to be considered. At present, ground-based remote sensing and monitoring technologies represent the distribution of sample points and have the ability to obtain series time resolution data. Statistical analyses or machine learning algorithms could be used to obtain the relationships between the sample plot and the surface data. The scale expansion analysis based on multisource data fusion technology could be a solution to explore the mechanism of the aerosol characteristics in forest regions affected by the microclimate. Second, the research area was in the developed area of southeastern China. During the period of economic development from 2010 to 2013, many local residents engaged in the planting of flower and tree seedlings for economic benefits, resulting in a large amount of deforestation. From 2013 to 2017, the local government invested 150 million Yuan to restore 226,637 ha of forest area in order to restore the forest ecological service level. Large areas of land were affected by human disturbance in a short time, which may affect the formation of aerosols in the air. Although rigorous screening criteria were set up to ensure the accuracy of forest type classification, the flow feature of aerosols cannot exclude the influence of other disturbances.In addition, this study has additional advantages. High-resolution AOD inversion data and detailed comprehensive ground forest resource survey data not only accurately depicted the spatial heterogeneity of the research objects but also explored the impact of different forest attributes on spatial heterogeneity. On the basis of integrating multisource data, the effect of diffusion radiation fertilization on regional forest park scale was verified, and the relationship between forest type and AOD concentration was discussed.


## 4.1. Relationship between Different Forest Types and AOD Distribution in the National Park

Siming Mountain Forest Park is located near the east coast of China, with a high forest coverage rate far from the urban center. The variation of AOD values may be caused by the combined action of moisture and particulate matter. In order to minimize the influence of water vapor, we selected eight remote sensing images taken during good and stable weather conditions between 10 : 00 and 11 : 00 am to avoid the condensation of water vapor caused by low temperature in the morning, although the average value of AOD in the study area is low. The AOD is the integral of the extinction coefficient of the medium in the vertical direction, which reflects the comprehensive effect of atmospheric scattering and surface reflection. We cannot distinguish the aerosol types. Based on the analysis of aerosol sources in four forest areas in Sichuan Province, China, it was found that the aerosols were mainly from local primary emission sources such as crust dust and biomass combustion, followed by the long distance transmission of aged organic aerosols (aged OOA) [51]. In addition to soil dust, meteorology, and human activities, plant species was also an important factor for the AOD variation in forest areas. The reactivity of plant VOCs with oxidants (such as OH, NO3 free radical, and O3) was higher than that of artificial VOCs, resulting in the generation of secondary organic aerosols (SOA). Zhu et al. (2017) found that the main source of SOA at the top of the exurban Tai Mountain in China was from an artificial source (11.6–42.6%), such as the use of solvents, fuel combustion, and vehicle exhaust as well as the biological source (6.5%) such as secondary oxidation and biomass burning [52]. However, the source of a large proportion of SOA was not identified. Although plants contribute to the generation of aerosols, the underlying surface of forest cover affects the physical processes in the atmospheric boundary layer. The high vegetation coverage increases the friction between the air and the ground, reduces wind speed, and diminishes transmission of pollutants from the ground and over long distances.The results of the Geodetector and cluster analysis in the cold and hot spots of this study in two years showed that the forest type was the main factor affecting AOD distribution in the exurban forest park, followed by NDVI and temperature. In the cold spots, more than half of the area was broad-leaved forest, while the percentage of broad-leaved forest decreased and the proportion of nonforest land area augmented in the hot spots. The result is similar to that of Deng et al. [53]. The maximized dust retention ability of the coniferous leaf for a single tree is higher than the broad-leaved tree species. However, there are many factors that influence the total dust retention quantity of the forest. As to the specific forest type, in addition to the blade characteristics, the factors influencing tree growth, such as temperature, humidity, soil nutrient, precipitation, and illumination, may affect the stand of leaf area index, thus further modifying the function of particulate matter retention. Moreover, the microstructure of the leaf surface is conducive to the deposition of particulate matter [54]. The forest increases the intensity of end flow in the surface atmosphere, thereby accelerating the deposition rate of particulate matter [55, 56]. Meanwhile, vegetation effectively absorbs pollutants such as SO2, NOx, and O3 in the air in favor of reducing the formation of secondary aerosols [51, 57]. In addition, the large amount of the particle matter in the air can be trapped on the leaf surface as compared to the small amount trapped in the leaf epidermal wax layer [58, 59]. The delayed particulate matter on the leaf surface is eluted to the ground with the precipitation leading to the recovery of its retention ability [60].From the perspective of the impact of aerosols on plants, the concentration threshold reversed the impact of aerosols on photosynthesis above which the net radiation (direct + scattering) effect of atmospheric aerosols would reduce the photosynthetic absorption of plants [61]. Xu quantified the AOD threshold of the impact of aerosol net radiation effect on ecosystem NPP in different regions of China and found that when the AOD was less than 1.5 at 550 nm in the southeast region, the aerosol net radiation effect increased NPP [16]. In this study, the average AOD value in hot spots was much lower than the threshold, and the results were in line with expectations. For example, the growth rate of tree species in hot spots was generally faster than that in cold spots. The average growth rates of the coniferous forest, broad-leaved forest, shrubs, and bamboos in cold spots were significantly different than in hot spots. The reason may be the diffusion radiation fertilization effect mentioned in the introduction, which increases the NPP and promotes growth in hot spots. However, there was not a significant difference in the growth rate of the coniferous broad-leaved mixed forest. The reason may be the large age variation between different tree species (large SD), which induced the different effects of diffusion radiation fertilization on tree species of different ages. The results of this study were consistent with other studies and verified the diffusion radiation fertilization effect at the regional forest park scale.


## 4.2. Limits and Prospects

There are two limitations in this study. First, the change of AOD is continuous with time. Landsat remote sensing satellite image data (a scene every 16 days) only obtains instantaneous information of the environment. Although the average value of four specific periods throughout the year was used as the research object in this study, the representativeness of inversion results still needs to be considered. At present, ground-based remote sensing and monitoring technologies represent the distribution of sample points and have the ability to obtain series time resolution data. Statistical analyses or machine learning algorithms could be used to obtain the relationships between the sample plot and the surface data. The scale expansion analysis based on multisource data fusion technology could be a solution to explore the mechanism of the aerosol characteristics in forest regions affected by the microclimate. Second, the research area was in the developed area of southeastern China. During the period of economic development from 2010 to 2013, many local residents engaged in the planting of flower and tree seedlings for economic benefits, resulting in a large amount of deforestation. From 2013 to 2017, the local government invested 150 million Yuan to restore 226,637 ha of forest area in order to restore the forest ecological service level. Large areas of land were affected by human disturbance in a short time, which may affect the formation of aerosols in the air. Although rigorous screening criteria were set up to ensure the accuracy of forest type classification, the flow feature of aerosols cannot exclude the influence of other disturbances.In addition, this study has additional advantages. High-resolution AOD inversion data and detailed comprehensive ground forest resource survey data not only accurately depicted the spatial heterogeneity of the research objects but also explored the impact of different forest attributes on spatial heterogeneity. On the basis of integrating multisource data, the effect of diffusion radiation fertilization on regional forest park scale was verified, and the relationship between forest type and AOD concentration was discussed.


## 5. Conclusion

This study overcame the low resolution of AOD remote sensing products and the unavailability problem of the forest vegetation vertical structure information from the optical remote sensing data by the use of high-resolution Landsat and ground survey data of forest resources. We analyzed the AOD spatial-temporal distribution pattern and its main influencing factors and explored the relationship of different types of forest canopies at the regional scale. It was found that the AOD distribution in the study area was clustered, and the forest type was one of the main impact factors. From 2010 to 2017, the average growth rate of broad-leaved forest, coniferous forest, bamboo, and shrub in hot spots was significantly higher than that in cold spots, while there was not a significant difference in mixed forests. The average growth rate of biomass in the cold and hot spots was higher in bamboo, coniferous forest, and mixed forest than in shrub and broad-leaved forest. In summary, at the regional forest park scale, the vegetation type had the closest interaction with AOD. The research results provide the guidance for the rational allocation of tree species to improve biomass and ecosystem service value of the exurban forest park.



---

*Source: 4942827-2019-10-30.xml*