Studying Forest Species Classification Methods by Combining PolSAR and Vegetation Spectral Indices

Abstract

Tree species are important factors affecting the carbon sequestration capacity of forests and maintaining the stability of ecosystems, but trees are widely distributed spatially and located in complex environments, and there is a lack of large-scale regional tree species classification models for remote sensing imagery. Therefore, many studies aim to solve this problem by combining multivariate remote sensing data and proposing a machine learning model for forest tree species classification. However, satellite-based laser systems find it difficult to meet the needs of regional forest species classification characters, due to their unique footprint sampling method, and SAR data limit the accuracy of species classification, due to the problem of information blending in backscatter coefficients. In this work, we combined Sentinel-1 and Sentinel-2 data to construct a machine learning tree classification model based on optical features, vegetation spectral features, and PolSAR polarization observation features, and propose a forest tree classification feature selection method featuring the Hilbert–Huang transform for the problem of mixed information on the surface of SAR data. The PSO-RF method was used to classify forest species, including four temperate broadleaf forests, namely, aspen (Populus L.), maple (Acer), peach tree (Prunus persica), and apricot tree (Prunus armeniaca L.), and two coniferous forests, namely, Chinese pine (Pinus tabuliformis Carrière) and Mongolian pine (Pinus sylvestris var. mongolica Litv.). In this study, some experiments were conducted using two Sentinel-1 images, four Sentinel-2 images, and 550 measured forest survey sample data points pertaining to the forested area of Fuxin District, Liaoning Province, China. The results show that the fusion model constructed in this study has high accuracy, with a Kappa coefficient of 0.94 and an overall classification accuracy of 95.1%. In addition, this study shows that PolSAR data can play an important role in forest tree species classification. In addition, by applying the Hilbert–Huang transform to PolSAR data, other feature information that interferes with the perceived vertical structure of forests can be suppressed to a certain extent, and its role in the classification of forest species, combined with PolSAR, should not be ignored.

1. Introduction

Tree species surveying is an important link in forest biomass estimation, which can provide data support for forest resource management, and at the same time, it is important for assessments of the ecological service value of vegetation [1,2,3]. The traditional forest species survey methods, represented by sample plot surveys, are difficult to apply to large-scale forest species survey tasks, because of difficulties related to data acquisition and spatial distribution, as well as their high costs [4,5]. With the development of remote sensing technology, remote sensing satellites have the advantages of a wide sensing range, continuous spatial distribution, and high temporal resolution, which facilitate tasks related to the large-scale investigation of tree species [4,5,6]. The spectral reflectance of multispectral images is widely used in forest tree species classification studies due to its ability to express the biochemical information of tree canopies [7,8,9]. However, multispectral imagery lacks the ability to perceive the vertical structure information of trees due to its insufficient transmittance, thus limiting its accuracy in tree species classification. Airborne LIDAR and satellite lasers detect the vertical structure of vegetation and understory topography by actively transmitting laser pulses, and are therefore able to capture differences in horizontal and vertical structural characteristics, such as canopy height and canopy width, between different tree species [10,11,12]. However, due to the small data coverage of the former and the footprint-type spot measurement of the latter, the remote sensing data are still spatially discontinuous, making it difficult to meet the needs of large-scale forest species survey missions over hundreds of square kilometers, or even larger areas. The Synthetic Aperture Radar (SAR) technique has become an effective remote sensing technique for obtaining the biophysical parameters of forests at regional and global scales, due to its imaging capability under all-weather and all-day conditions, and its sensitivity to the geometric and physical properties of the target [13,14,15,16,17]. The method realizes the transmission to the vegetation area, as well as a certain depth from the ground surface, by actively transmitting microwaves, which can realize the perception of the vertical structure information of the vegetation in a large-scale area. Therefore, this method has been widely applied to the classification of crop planting structure [18,19], forest species classification [13,20,21], and other related studies. In recent years, studies examining forest species classification by combining optical remote sensing images and SAR data have achieved notable results due to this method’s ability to combine the perception of the biochemical information of tree canopies and the vertical structure of trees [22,23]. Since the backscattering coefficients of SAR images contain topographic information, canopy information, and other noise information that interferes with the vertical structure perception of trees, the classification accuracy is limited to a certain extent. The Hilbert–Huang transform (HHT), with the advantages of having an adaptive decomposition property, no requirement for basis function selection, and its time–frequency distribution, can analyze smooth electromagnetic wave signals and nonlinear non-smooth electromagnetic wave signals [24,25]. Good research results have been achieved in seismic data analysis [26,27,28], gravity detection [29,30,31], and bridge deformation [32,33,34]. In our study, we propose a tree species classification method for multi-source remote sensing images combined with HHT. The method decomposes the electromagnetic wave information of SAR images into multiple intrinsic modal functions (IMFs) based on HHT, aiming to reduce to some extent the limiting effect of the coupling of the backscatter coefficient and multi-class information of SAR images on the accuracy of tree species classification, so as to emphasize the contribution of SAR data in tree species classification research. On this basis, we realized forest species classification using optical and synthetic aperture radar (SAR) images, to take into account the vertical structure of forest tree species and differences in vegetation spectral indices. The objectives of this study were (1) to investigate the degree of contribution of different orders of eigenmode functions to forest species classification, (2) to select the best eigenmode function of SAR image information and combine the reflectance information of multispectral images and spectral indices of multiple types of vegetation, and (3) to construct a wide-ranging forest species classification model based on machine learning.

2. Study Area and Data

2.1. Overview of the Study Area

The Fuxin region is located in the northwest of Liaoning Province, China, with an area of 10,355 km² and a forest cover of about 2901 km². The northern part of the Fuxin region is the Horqin Sandland, and the eastern part is in the Liaohe Plain. It is a transition zone between the Inner Mongolian steppe and the stony mountainous terrain of North China. The geomorphological pattern of the region is high in the northern part of the west and low in the southern part of the east, as shown in Figure 1. The Fuxin region’s climate is a semi-arid continental monsoon climate characteristic of the north temperate zone; the main feature of the climate is precipitation, but the spatial and temporal distribution is extremely uneven. The rain and heat simultaneously aid the growth of vegetation in the region, providing favorable conditions. The region is characterized by a variety of natural environments and a rich variety of forest species. It is worth noting that the region is at the junction of the semi-humid and semi-arid zones and that the forests of the Fuxin region are also an important part of the windbreak forests in North China, resulting in the distribution of both natural and planted forests, as well as coniferous and broadleaf forests, in the region [35]. Therefore, the classification of forest species in Fuxin is of great significance to the study of carbon storage, the carbon cycle, and desertification control in the region.

2.2. Research Data

Two Sentinel-1B images and four Sentinel-2 images covering the study area in 2021 were selected as data sources for our forest species classification study. Sentinel-1B and Sentinel-2 are multispectral imaging sensors operated by the European Space Agency and provide open, freely accessible data. (https://dataspace.copernicus.eu/, accessed on 29 August 2021). The Sentinel-1B image has a spatial resolution of 10 m, and the data have been processed by polarization decomposition, refined Lee filtering, and terrain radiation correction. The ground truth for the data used in the study was obtained by the researchers from the annual field survey, which took place in August 2021. The survey was based on 20 m × 20 m square plots covering a range of information on species composition, structural attributes, and health, and a total of 550 sample data points were measured for machine learning model training and validation, which were uniformly distributed within the forest cover of the study area and included all dominant tree species in the study area. According to the survey results, there were six dominant tree species in the study area: aspen (Populus L.), Chinese pine (Pinus tabuliformis Carrière), maple (Acer), Mongolian pine (Pinus sylvestris var. mongolica Litv.), peach tree (Prunus persica), and apricot tree (Prunus armeniaca L.).

2.3. Polarized Synthetic Aperture Radar Observations

In the polarization decomposition theory, the scattering matrix is the basis for the study of scattering of ground objects, containing both amplitude and phase information, and is a complex matrix linking the electric field strengths of the incident wave and the scattered return wave. However, the scattering matrix can only describe the scattering signature of pure scatterers, and it cannot describe that of distributed scatterers. The scattering of distributed scatterers can only be described statistically. In order to reduce the effect of noise, the scattering of distributed scatterers is commonly analyzed in dual-polarization SAR systems using the covariance matrix ($C_2$ matrix), which is a second-order scattering information quantity generated by the spatial averaging of the target vector of the dual-polarization scattering matrix, k [36], which is generated by spatial averaging, as shown in Equation (1), and the $C_2$ matrix as shown in Equation (2) [37].

$$k={\left[\begin{array}{cc}〈S_{VV}〉& 〈S_{VH}〉\end{array}\right]}^T$$

(1)

$$C_2=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]=\left[\begin{array}{cc}〈{\left|S_{VV}\right|}^2〉& 〈\left|S_{VV}{S}_{VH}^{*}\right|〉\\ 〈\left|S_{VH}{S}_{VV}^{*}\right|〉& 〈{\left|S_{VH}\right|}^2〉\end{array}\right]$$

(2)

where * denotes the complex conjugate, and according to the scattering symmetry, there is $〈\left|S_{VV}{S}_{VH}^{*}\right|〉=〈\left|S_{VH}{S}_{VV}^{*}\right|〉$ = 0. Equation (2) can be rewritten as Equation (3). In Equation (3), $\sigma =〈{\left|S_{VV}\right|}^2〉$, $\epsilon =〈{\left|S_{VH}\right|}^2〉/〈{\left|S_{VV}\right|}^2〉$.

$$C_2=\sigma \left[\begin{array}{cc}\sigma & 0\\ 0& \epsilon \end{array}\right]$$

(3)

Based on Equation (3), the $C_2$ matrix is a semi-positive definite Hermite matrix that can be decomposed into a weighted sum of two mutually orthogonal matrices, as shown in Equation (4). In Equation (4), λi and ξi denote the real eigenvalues and eigenvectors of the $C_2$ matrix, respectively, and Ci denotes the independent covariance matrices of rank 1, denoting a scattering mechanism and the corresponding eigenvalues λi. The expression of the eigenvalues of the $C_2$ matrix is shown in Equation (5). Based on Equation (5), the dual-polarization decomposition can be expressed as shown in Equation (6).

$$C_2=\sum _{i=1}^2{\lambda }_i{C}_i={\lambda }_1{\xi }_1{\xi }_1^{*}+{\lambda }_2{\xi }_2{\xi }_2^{*}$$

(4)

$$\begin{array}{c}{C}_2{\lambda }_1=\sigma \\ {\lambda }_2=\sigma \epsilon \end{array}$$

(5)

$$C_2={\lambda }_1{C}_1+{\lambda }_2{C}_2=\sigma \left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]+\sigma \epsilon \left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right]$$

(6)

From Equation (6), we can see that $C_1$ is the other scattering, which is the result of the joint action of odd and even scattering, and ${\lambda }_1$ represents the intensity of this other scattering, while $C_2$ is only related to VH and represents the multiple scattering process. ${\lambda }_2$ is used to represent the intensity of multiple scattering. In order to characterize the disorder in the scattering process of the feature, the scattering entropy (H) is introduced, which is defined as Equation (7). It is important to note that $P_i={\displaystyle \frac{{\lambda }_i}{{\lambda }_1+{\lambda }_2}}$.

$$H=\sum _{i=1}^2-P_i{\log }_2{P}_i$$

(7)

The scattering entropy (H) characterizes the degree of randomness of target scattering. When H = 0, it indicates that the system is in a fully polarized state with only one scattering mechanism. When H tends to 0, the system is close to a fully polarized state, with one of the scattering mechanisms dominating; when H approaches 1, the ${\lambda }_1$ and ${\lambda }_2$ eigenvalues are close in size, and the system is close to a completely unpolarized state; when H = 1, the system is in a completely unpolarized state, and the target scattering is completely random noise. The scattering angle represents the average scattering mechanism from surface scattering to dihedral angle scattering and is defined as Equation (8).

$$\overline{a}=P_1{a}_1+P_2{a}_2$$

(8)

where $\overline{a}$ is the scattering angle. When $\overline{a}$ = 0°, it corresponds to isotropic surface scattering, and when $\overline{a}$ = 90°, it corresponds to isotropic dihedral scattering. A Sentinel-1 Single Look Complex (SLC) can generate the $C_2$ matrix, and then based on the $C_2$ matrix, we can choose the H-alpha decomposition. Then, we can finally obtain the three polarization features of the scattering angle, entropy, and anisotropy.

2.4. Spectral Index of Vegetation in Multispectral Images

The multispectral images used in this study were selected from Sentinel-2 L2A level data in July 2021, during the period of lush vegetation in the study area with less than 10% cloud cover. The Sentinel-2 imagery has 13 bands, which can provide data at up to a 10 m spatial resolution. In order to ensure the consistency of the spatial resolution of the images used in this study, we sampled all the bands to 10 m and went through the band fusion and region cropping preprocessing steps. Among the spectral features, the band information can reflect the differences between the absorption and reflectance intensities of the spectra of the leaf blades of different species [38], and several vegetation indices, including chlorophyll content [39], leaf moisture content [40], and vegetative soil conditions [41], have also been proven to be able to characterize the differences in the branching of the branches and leaves of different tree species. In our study, we selected 11 vegetation spectral features, including the Normalized Vegetation Index (NDVI), Normalized Humidity Index (NHI), Optimized Soil-Adjusted Vegetation Index (OSAVI), Chlorophyll Index (CGI), Green Index (GI), Moisture Stress Index (MSI), Red Edge Normalized Vegetation Index (reNDVI), Enhanced Vegetation Index (EVI), and Ratio Vegetation Index (RVI), along with the Chlorophyll Content Index (LCCI) and Normalized Difference Red Edge Index (NDRESWIR), as shown in

Table 1. Vegetation-based spectral variables [7].

Vegetation Spectral Characteristics	Descriptive
Normalized Vegetation Index (NDVI)	$NDVI={\displaystyle \frac{B_8-B_4}{B_8+B_4}}$
Normalized Humidity Index (NHI)	$NHI={\displaystyle \frac{B_{11}-B_3}{B_{11}+B_3}}$
Optimized Soil-Adjusted Vegetation Index (OSAVI)	$OSAVI={\displaystyle \frac{B_8-B_4}{B_8+B_4+0.16}}$
Chlorophyll Index (CGI)	$CGI={\displaystyle \frac{B_8}{B_3+B_4-B_5}}$
Green Index (GI)	$GI={\displaystyle \frac{B_3}{B_4}}$
Moisture Stress Index (MSI)	$MSI={\displaystyle \frac{B_{11}}{B_8}}$
Red Edge Normalized Vegetation Index (reNDVI)	$reNDVI={\displaystyle \frac{B_8-B_4+B_5}{B_8+B_4-B_5}}$
Enhanced Vegetation Index (EVI)	$EVI={\displaystyle \frac{2.5(B_8-B_4)}{B_8+6B_4-7.5B_2+1}}$
Ratio Vegetation Index (RVI)	$RVI={\displaystyle \frac{B_8}{B_4}}$
Chlorophyll Content Index (LCCI)	$LCCI={\displaystyle \frac{B_4-B_7}{B_4-B_5}}$
Normalized Difference Red Edge Index (NDRESWIR)	$NDRESWIR={\displaystyle \frac{B_4-B_6-B_{12}}{B_4-B_6\mp B_{12}}}$
Spectral Band	B1–B12

. B1–B12 represent the spectral values of Sentinel-2 satellite images in band 1 and band 12, respectively; meanwhile, B1–B12 as optical features were also used as input data for forest species classification study.

3. Forest Species Classification Methods

3.1. The Hilbert–Huang Transform

The Hilbert–Huang transform (HHT) is a signal processing method that mainly consists of two parts, empirical modal decomposition (EMD) and the Hilbert transform (HT) [42]. The empirical modal decomposition can decompose the microwave signals of SAR data into several intrinsic mode function (IMF) components, and each IMF component is arranged in order, from high frequency to low frequency. The Hilbert transform (HT) is applied to each IMF component on the basis of EMD to obtain the amplitude and frequency of each IMF component. The HHT is based on the concept of intrinsic mode functions and the idea that any SAR data microwave signal can be decomposed into several IMFs, so that the instantaneous frequency has a real physical meaning. The HHT can adaptively decompose the signal according to its own time scale characteristics, and there are no human factors in the process of decomposition, meaning that the characteristics of the data itself are retained in the decomposition process, and the resulting individual IMF components contain the relevant information of the original signal.

3.1.1. Empirical Modal Decomposition

The empirical modal decomposition method, used to decompose SAR data microwave signals, should be based on the following three assumptions [43]: First, the SAR data microwave signal has at least two extreme points (one extremely large and one extremely small). Second, the feature scale is defined as the time interval between neighboring extreme points. Finally, if the signal has no extreme points but only “curved points”, then before decomposing it, one must first differentiate it one or more times to obtain the extreme points, and then integrate all the results to obtain the corresponding components. For this, the following steps are required [44]: The signal is known as $x\left(t\right)$. To determine all the extreme values of the microwave signal and fit them using the cubic spline interpolation algorithm, the upper envelope is obtained based on the fitted maximum and minimum values, respectively, $x_{max}\left(t\right)$ and lower envelope $x_{min}\left(t\right)$. Then, we calculate the upper and lower mean values of the envelope, denoted as $m_1\left(t\right)$, as shown in Equation (9). Using SAR data microwave signal data series $x\left(t\right)$, subtract the upper and lower mean values of the envelope $m_1\left(t\right)$ to obtain a new sequence $h_1\left(t\right)$, as shown in Equation (10).

$$m_1\left(t\right)=\left[x_{max}\left(t\right)+x_{min}\left(t\right)\right]/2$$

(9)

$$h_1\left(t\right)=x\left(t\right)-m_1\left(t\right)$$

(10)

If the $h_1\left(t\right)$ IMF condition is satisfied, then $h_1\left(t\right)$ is the first IMF component. If the $h_1\left(t\right)$ IMF condition is not satisfied, then treat $h_1\left(t\right)$ as the original signal sequence, and then calculate the $h_1\left(t\right)$ and the mean value of the envelope of m1.1(t), and calculate the $h_1\left(t\right)$ and m1.1(t), and the difference between the two, to obtain $h_{1.1}\left(t\right)=h_1\left(t\right)-m_{1.1}\left(t\right)$. Repeat the above steps k times until $h_{1.k}\left(t\right)$ satisfies the IMF definition and the first eigenmode function is obtained, denoted as $c_1\left(t\right)$. Usually, $c_1\left(t\right)$ represents the high-frequency part of the microwave signal data sequence of the SAR data, and it can be obtained by subtracting the original signal $x\left(t\right)$. Subtract $c_1\left(t\right)$ to obtain the remaining sequence $r_1\left(t\right),$ that is, $r_1\left(t\right)=x\left(t\right)-c_1\left(t\right)$. Take $r_1\left(t\right)$ as a new original sequence and follow the above steps to extract up to the nth eigenmode function, $c_n\left(t\right)$. After extracting the last eigenmode function, the sequence $r_n\left(t\right)$ becomes a monotonic sequence. By integrating the decomposed components, the original sequence is obtained $x\left(t\right)$, as shown in Equation (11).

$$x\left(t\right)=\sum _{i=1}^n{c}_i\left(t\right)+r_n\left(t\right)$$

(11)

3.1.2. Hilbert Transform

For any continuous real signal $X\left(t\right)$, its Hilbert transform is shown in Equation (12), which is defined as a generalized integral over the real number field with respect to $\tau$, the generalized integral, where $\tau$ is the integration variable used to traverse all time points in the integration process. Equation (13) represents the Hilbert inverse transform.

$$Y\left(t\right)={\displaystyle \frac{1}{\pi }}{\int }_{-\infty }^{\infty }{\displaystyle \frac{X\left(\tau \right)}{t-\tau }}d\tau$$

(12)

$$X\left(t\right)={\displaystyle \frac{1}{\pi }}{\int }_{-\infty }^{\infty }{\displaystyle \frac{Y\left(\tau \right)}{t-\tau }}d\tau$$

(13)

$Y\left(t\right)$ and $X\left(t\right)$ can form a set of conjugate complex pairs, which, in turn, leads to the analytic signal $Z\left(t\right)$, as shown in Equation (14), where $a\left(t\right)$ is the instantaneous amplitude, and $\theta \left(t\right)$ is the phase.

$$\begin{array}{c}Z\left(t\right)=X\left(t\right)+Y\left(t\right)i=a\left(t\right)e^{\theta \left(t\right)i}\\ a\left(t\right)=\sqrt{[{X\left(t\right)}^2+{Y\left(t\right)}^2]}\\ \theta \left(t\right)=arctan{\displaystyle \frac{Y\left(t\right)}{X\left(t\right)}}\end{array}$$

(14)

Furthermore, the continuous real signal $X\left(t\right)$ can be expressed as Equation (15), which, in turn, leads to the instantaneous frequency $f\left(t\right)$.

$$\begin{array}{c}X\left(t\right)=Re{\displaystyle \sum _{j=1}^n}{a}_j\left(t\right)e^{\theta \left(t\right)i}=Re{\displaystyle \sum _{j=1}^n}{a}_j\left(t\right)e^{i\int {\omega }_j\left(t\right)dt}\\ f\left(t\right)={\displaystyle \frac{1}{2}}{\displaystyle \frac{d\theta \left(t\right)}{dt}}\end{array}$$

(15)

where Re denotes taking the real part, the expansion is called the Hilbert spectrum, and the residual component is omitted in the $r_n\left(t\right)$ equation. From the above analysis, it can be seen that the HHT method is a combination of EMD and the Hilbert transform.

3.2. Forest Species Classification Models

Considering the complexity and nonlinearity of forest species classification based on multi-source remote sensing data, it is difficult for existing generalized machine learning models to achieve the accurate classification of forest species in a natural state. Therefore, we proposed a multi-source remote sensing forest species classification method based on the above vegetation spectral features and the Hilbert–Huang transform of PolSAR signals. The model consists of two parts, as shown in Figure 2: (1) the feature transformation module—used to transform PolSAR features and multispectral remote sensing reflectance to obtain new features that can more accurately represent tree information, such as the vegetation-based spectral variables, etc.—and (2) the stochastic forest classifier module. Because the number of data samples is limited by the amount of measured forestry survey data in our study, the stochastic forest classifier, with its excellent classification effect under small-sample conditions, was chosen. Additionally, the random forest classifier, which has an excellent classification effect under small-sample conditions, was chosen as the base model, and hyperparameter optimization was based on the particle swarm algorithm [45].

The decision tree method is a decision classification method based on tree architecture, which mainly includes features such as selection, branch generation, and pruning [46], as shown in Figure 2. The random forest classifier section consists of several decision trees. The decision tree process starts from the root node, for partitioning the solution and performing multiple decompositions at the nodes. Different decision trees will have different criteria for dividing the data samples into two subsets. We iterate repeatedly according to the above logic until the termination of splitting is achieved. A wide variety of decision tree models are currently being developed and improved, while the CART decision tree algorithm has been widely used due to its ease of implementation and the relatively small amount of expertise and experience required [47]. Random forest is a parallel algorithm based on CART decision trees; it has a strong learning ability. Random forest trains multiple CART decision trees and then learns by summarizing the resultant experience of each decision tree in a certain combination. A random forest consists of many decision tree classifiers. Each decision tree classifier can be defined as $h(x,{\theta }_i)$, where x represents the input features, and ${\theta }_i$ represents an independent set of identically distributed features. Each decision tree of the random forest votes independently for classification, and the final classification result is obtained via the combined vote [48]. The process of random forest classification is shown in Equation (16).

$$\begin{array}{c}H=\left\{h_1, h_2, h_3, \dots , h_k\right\}\\ h\left(x,{\theta }_i\right):X\to Y,i=1,2,\dots ,k\end{array}$$

(16)

In the random forest classification process, the result with the highest number of votes in the decision tree is used as the classification result, which is calculated as in Equation (17). Due to random forest classification’s random nature, the risk of overfitting can be reduced by averaging the decision trees, making it highly resistant to overfitting. Random forest classification also has strong stability; even in the midst of a large fraction of feature loss, as long as no error occurs in more than half of the decision tree classifiers, accuracy can still be maintained, including in cases where a new data point has been added, as the presence of this new data point in the dataset will likely affect only one decision tree, hardly affecting all decision trees [49]. For unbalanced data, such classification has an ability to balance error. However, because random forests require the construction of multiple decision trees, they are more time-consuming and costly to compute [50]. The combination of the above features and the fact that the trees are independent of each other during the training of the random forest gives this classification process the ability to parallelize power and be relatively fast in remote sensing image dataset classification problems, having some applicability [51,52,53,54]. Thus, we selected random forest classification as the base model for tree classification.

$$y^{*}=argmax\sum _{y\in Y,j\in H}I(f\left(x|t\right)=y)$$

(17)

Due to the complex nonlinear relationship between the vegetation features of the multi-source remote sensing images selected in this study and the actual tree classification, we use the particle swarm algorithm (PSO) to find optimal solutions for key parameters such as n_estimators, max_features, max_depth, etc., on the random forest base model. The method, which uses population fitness data to find the optimal solution of a given problem, can be used for the optimization of nonlinear problems [55,56]. The algorithm treats each solution to the optimization problem as a random particle, and the position and velocity of each particle are random. During the iteration process, each particle records and updates its current position and velocity, recording its own optimal position and the optimal position of the population [57,58]. For example, to find an optimal problem in n dimensional space, there are m particles forming a population; the position of the 1st particle is written as $x_i=(x_1,x_2,\dots ,x_n)$, and the velocity is denoted as $v_i=(v_1,v_2,\dots ,v_n)$, while the iterative formula for updating the velocity and position of each particle in the population is as follows [59,60]:

$$v_i(t+1)=\omega v_i(t)+c_1{r}_1(p_i(t)-x_i(t))+c_2{r}_2(p_g(t)-x_i(t))$$

(18)

$$x_i(t+1)=x_i(t)+v_i(t+1)$$

(19)

where i = 1, 2, 3, …, m; $P_i$ is the individual pole position; $P_g$ is the global pole position; ω is the initial value of inertia weight; $C_1$ and $C_2$ are the acceleration coefficients; and $r_1$ and $r_2$ are random numbers between 0 and 1. For the RF model, the inclusion of the particle swarm algorithm with global convergence can better ensure the rationality of parameter optimization [59].

4. Experiments and Results

4.1. Forest Species Classification Results

In the regional forest species classification study, we accurately classified the forest species in the study area based on multi-source remote sensing forest species classification methods [61,62]. The model input variables were divided into two parts, namely, multi-source remote sensing data and measured forest survey data. Among them, the multi-source remote sensing data included the Normalized Vegetation Index (NDVI), the Normalized Humidity Index (NHI), the Optimized Soil Conditioning Vegetation Index (OSAVI), the Chlorophyll Index (CGI), the Greenness Index (GI), the Moisture Stress Index (MSI), the Red Edge Normalized Vegetation Index (reNDVI), the Enhanced Vegetation Index (EVI), and the Ratio Vegetation Index (RVI), along with the Chlorophyll Content Index (LCCI), the Normalized Difference Red Edge Index (NDRESWIR), and 11 other spectral vegetation features, in addition to the reflectance of multispectral data and second-order eigenvalues of five polarization observation variables, such as the main diagonal elements of the C2 matrix (C11 and C22) [63], and the H-alpha decomposition parameter (alpha, anisotropy, and entropy) modal function (IMF) [64]. Before random forest species classification, we first performed feature selection based on random forest importance ranking, as shown in Figure 3. In Figure 3, reNDVI and MSI were significantly less important than other features in forest species classification, so we removed these two features. We input the tree species information obtained from the 550 measured forest survey sites described in Section 2.2 and the 26 multi-source remote sensing features mentioned above into the tree species classification module at a ratio of 0.7:0.3, before optimizing the random forest model for parameter tuning based on the particle swarm algorithm. Among them, the tree species information is the type of tree species in the measured sample sites, and the multi-source remote sensing information includes the B1–B12 bands, the vegetation-based spectral variables shown in Table 1, and the five PolSAR features mentioned above. Then, we selected the optimal random forest model parameters for making inferences regarding the tree species in the study area. Regarding the model parameters, n_estimators are 100, max_features is all features.

Table 2. A confusion matrix for the classification of tree species for multi-source remote sensing forest species classification methods.

Forecast Category	Reference Category						User Accuracy
	Aspen	Maple	Mongolian Pine	Peach Tree	Chinese Pine	Apricot Tree
Aspen	42	0	1	0	0	0	0.976
Maple	1	27	1	0	0	1	0.900
Mongolian pine	0	0	34	3	0	0	0.918
Peach tree	0	0	0	15	0	0	1.000
Chinese pine	0	0	0	0	22	0	1.000
Apricot tree	0	1	0	0	0	17	0.944
Producer accuracy	0.976	0.964	0.944	0.833	1.000	0.944
Kappa	0.94
Overall classification accuracy	0.951

represents the classification accuracy of multi-source remote sensing forest species classification methods, with overall classification accuracy and the Kappa coefficient as the evaluation criteria. As shown in Table 2, the overall classification accuracy of the model is 0.951, and the Kappa coefficient is 0.94. This shows the effectiveness of multi-source remote sensing forest species classification methods in the task of classifying forest species in the region.

4.2. Spatial Distribution of Forest Tree Species in the Study Area

In our tree species classification study, based on the optimal input feature combinations and training hyperparameters obtained in Section 4.1, as shown in Figure 4, the white-colored part of the Fuxin area is the non-forested area, consisting of towns, waters, and farmland. The non-forested area was not involved in our forest species classification based on expert visual interpretation. In addition, in order to highlight the differences between different tree species, we display the six dominant tree species in different colors in the forest species distribution map of the study area, shown in Figure 3, including aspen, Chinese pine, maple, Mongolian pine, peach tree, and apricot trees.

As analyzed in Figure 5, the forests in the study area are mainly distributed in the southwest, while the northeast has a small amount distributed, and the forests in the rest of the area are sporadically distributed. The above six dominant tree species in the study area are generally mixed forests, among which maple, aspen, and Chinese pine have a wider distribution area. Particularly noteworthy is that, in the northeastern forest area of the study area, Mongolian pine is the dominant tree species, and this area is also the area with the largest distribution of Mongolian pine in the study area. This is due to the fact that this area is a plantation forest, which is an important part of China’s “Three-North” protection forest, where a large number of Mongolian pine are planted in the windbreak forest, forming a unique Mongolian pine forest area in Fuxin.

5. Discussion

In Section 5.1, we discuss our first set of experiments, conducted to assess the extent to which individual remote sensing image features contributed to the forest tree species classification task. We provide a discussion of our second set of experiments, based on the Hilbert–Huang transform, in Section 5.2. We used these experiments to assess the change in the degree of contribution of polarized observational variables to the forest tree species classification task before and after the application of the Hilbert–Huang transform, as well as the interaction of multiple factors. Finally, in Section 5.3, we discuss the limitations of our proposed forest tree species classification method, and future research.

5.1. Ablation Study

In our study, in order to highlight the effects of various types of features on forest tree species classification, we designed a feature ablation experiment to quantitatively characterize each feature’s contribution to classification accuracy. In the feature space constructed by the optimal feature set, all the features were merged into three categories according to the differences in the feature descriptions on the image, specifically optical features, vegetation spectral features, and PolSAR features, where, optical features are the reflectance of the bands B1–B12 included in the Sentinel-2 image. Vegetation spectral features are the eleven vegetation-based spectral variables mentioned in Table 1. PolSAR features are the five polarization-based observational variables mentioned in Section 4.1. Under the same experimental conditions used to construct the model, featuring the experimental parameters and dataset division mentioned in Section 4.1, the three classes of features were cyclically deleted from the feature variables of the constructed model. The objective was to quantify the contribution of three remote sensing image features in forest species classification studies. If the classification accuracy was severely reduced after subtracting a certain type of remote sensing image feature, it indicated that this feature contributed more to the forest species classification task. It is worth noting that in this section of the experiment, no feature selection was performed in order to quantify the degree of contribution of optical features, vegetation spectral features, or PolSAR features. A total of 28 features were involved in the ablation study. The recognition accuracies of the three models and the changes in model accuracy before and after the deletion of features are shown in Figure 6. The model accuracy is most affected when no vegetation spectral features are introduced, as shown in Figure 6e,f. The overall classification accuracy decreased from 93.3% to 88.48%; from Figure 6e, it can be seen that the discriminative abilities of the peach tree, apricot tree, and oil pine are all greatly affected by the vegetation spectral features. This is due to the fact that there is a large gap between the peach tree, apricot tree, and camphor pine in the vegetation spectral indices, and such features can reflect these divergence gaps. The overall classification accuracy of our method is only second to that of when the PolSAR features extracted from SAR imagery are deleted. The decrease in overall classification accuracy after deleting the PolSAR features from the SAR images is only second to the vegetation spectral features, with a decrease of 1.78%. It can be seen that the canopy vertical structure feature extracted from the SAR imagery makes an important contribution to the classification results.

5.2. The Validity of the Hilbert–Huang Transform

The SAR image polarization observation variable contains terrain information, canopy information, soil moisture, and other information, so some noise information that interferes with the perception of the vertical structure of the tree limits the accuracy of the tree species classification to a certain extent. The Hilbert–Huang transform is usually used in the decomposition of time series signals, but in the processing of non-time series electromagnetic wave signals, a series of intrinsic modal functions that retain the original signal characteristics is generated by the decomposition and transformation of original electromagnetic wave signals, as shown in Figure 7. An electromagnetic wave is decomposed and transformed to generate a series of eigenmode functions that retain the original signal characteristics. Furthermore, some specific orders of eigenmode functions of the polarized observed signals will make a better contribution than the original signals in the tree species classification task. The reason for this phenomenon is that the eigenmode functions of specific orders play a role in suppressing the topographic information and soil moisture, etc., while preserving the characteristics of the original signal, which is more conducive to machine learning-based forest tree species classification. In order to investigate the effect of different orders of eigenmode functions in forest species classification, we used different orders of eigenmode functions to conduct forest species classification experiments in the study area of Fuxin comprising mixed coniferous and broad forest, as shown in Figure 7.

As shown in Figure 7, the first layer of each subplot is the original-input PolSAR signal. From the second layer onwards, there is the first–sixth-order intrinsic modal function (IMF), and the last layer is the residual signal. As can be seen from the results, we find that as the order of the intrinsic modal function increases, starting from the third-order intrinsic modal function (IMF), the signal waveforms appear to differ greatly from the original signal. Therefore, in our study, only the first- and second-order intrinsic modal functions were used to compare the original PolSAR signals, to explore the effectiveness of the Hilbert–Huang transform in tree classification studies. In view of this, based on the measured forest survey sample data and the PolSAR data of corresponding spatial locations mentioned in Section 2.2, we performed tree species classification experiments, respectively, and the results are shown in

Table 3. Results of PolSAR based on Hilbert–Huang transform in forest tree species classification.

Experimental Data	Average User Accuracy	Average Producer Accuracy	Kappa Coefficient (i.e., Height of Force)
Original PolSAR	0.564	0.554	0.488
IMF1	0.504	0.567	0.459
IMF2	0.792	0.879	0.779

. From the results, it can be seen that the random forest model based on particle swarm optimization acquires higher tree species classification accuracy than the second-order IMF, which shows the effectiveness of our proposed tree species classification method based on the Hilbert–Huang transform.

5.3. Limitations and Future Research

We selected a particle swarm optimization-based random forest (PSO-RF) model using the Normalized Vegetation Index (NDVI), the Normalized Humidity Index (NHI), the Optimized Soil Conditioning Vegetation Index (OSAVI), the Chlorophyll Index (CGI), the Greenness Index (GI), the Enhanced Vegetation Index (EVI), the Ratio Vegetation Index (RVI), the Chlorophyll Content Index (LCCI), the Normalized Difference Red Edge Index (NDRESWIR), the reflectance of the multispectral data, and the third-order eigenmode function of the diagonal elements of the C2 matrix (C11 and C22), along with the H-alpha decomposition parameters (alpha, anisotropy, and entropy) (IMFs), vegetation spectral features, and polarization observation variables to classify the forest species in the study area. Since different multi-source remote sensing features contribute differently in the task of forest species classification [65], we explored the contribution of different types of remote sensing image features in the study of tree species classification based on ablation experiments. Meanwhile, in this study, we proposed a machine learning-based classification method (multi-source remote sensing forest species classification methods) using Sentinel-1 polarized observation data and Sentinel-2 multispectral data for forest tree species classification. In this study, a forest species classification experiment was conducted in the Fuxin area of Liaoning Province, China, as the study area. Some of our findings, such as the spatial distribution of forests and the species composition of forests in the northeastern part of the Fuxin area can be corroborated in the previously reported study by Wang et al. [35]. In addition, Yan et al. utilized MODIS images from 2012 to achieve forest classification in Northeast China by constructing a vegetation spectral index. In their results, the distribution of coniferous and broadleaf forests in the Fuxin area had some similarities with our study [66]. In a recently reported study, Liu et al. used GF-1 optical and Sentinel-1 images to map the distribution of coniferous forests in the northwest of Liaoning Province, China. Among them, the distribution of coniferous forests in the Fuxin area is very similar to the two coniferous forest distributions in our study [67]. However, the method in this paper still has some limitations. The method relies on a large amount of actual forest survey data and remote sensing data. In our case, although the number of samples is large, the proportions between samples of different tree species are not completely balanced. This may affect the accuracy of model calibration and estimation. In addition, natural conditions such as soil moisture and vegetation foliar water content may also affect the estimation results, and although we suppressed these signals to a certain extent based on the Hilbert–Huang transform, the effects still exist. Future work will focus on testing the method’s ability to classify other vegetation types (e.g., crops, etc.) as well as other types of forests (e.g., tropical rainforests). The use of other radar frequencies (e.g., LT-1 in the L-band and PolSAR data from the airborne radar in the P-band) and high-resolution multispectral data (e.g., GF-series satellites and ZY-series satellites) will also be investigated to test and analyze the method in different locations in future studies.

6. Conclusions

In this study, we selected optical features, vegetation spectral indices, and PolSAR polarization observation variables and applied them to temperate mixed conifer forests in the Fuxin area of Liaoning Province, China, to analyze their contributions to the classification of forest species in this area. Based on this, we realized forest species classification in the study area via the use of a stochastic forest classification model of particle swarm optimization (PSO). Through our study, we found that the PolSAR polarization observation variable can provide valuable information for forest species classification. The polarization decomposition of backscatter coefficients of SAR images to obtain canopy vertical structure information in feature space analysis experiments proved to be an important contribution to the task of forest tree species classification, demonstrating the effectiveness of the PolSAR polarization observation variable in forest species classification. In addition, by applying the Hilbert–Huang transform to PolSAR polarization observation variables, other feature information that plays an interfering role in perceiving the vertical structure of forests can be suppressed to a certain extent, and its role in forest tree species classification, along with that of PolSAR, should not be neglected. The model constructed shows high accuracy in forest species classification experiments, with an overall accuracy of 95.1% and a Kappa coefficient of 0.94. The multi-source remote sensing forest species classification approach integrates multiple sources of remote sensing data and significantly improves the accuracy of tree species identification compared to classification models derived from a single data source. Therefore, the method is suitable for large-scale forest species surveys of temperate mixed conifer forests based on remote sensing images.