Feature-Level Fusion between Gaofen-5 and Sentinel-1A Data for Tea Plantation Mapping

: The accurate mapping of tea plantations is signiﬁcant for government decision-making and environmental protection of tea-producing regions. Hyperspectral and Synthetic Aperture Radar (SAR) data have recently been widely used in land cover classiﬁcation, but e ﬀ ective integration of these data for tea plantation mapping requires further study. This study developed a new feature-level image fusion method called LPPSubFus that combines locality preserving projection and subspace fusion (SubFus) to map tea plantations. Based on hyperspectral and SAR data, we ﬁrst extracted spectral indexes, textures, and backscattering information. Second, this study applied LPPSubFus to tea plantation mapping with di ﬀ erent classiﬁcation algorithms. Finally, we compared the performance of LPPSubFus, SubFus, and pixel-level image fusion in tea plantation mapping. Feature-level image fusion performed better than pixel-level image fusion. An improvement of about 3% was achieved using feature-level image fusion compared to hyperspectral data alone. Regarding feature-level image fusion, LPPSubFus improved the overall accuracy by more than 3% compared to SubFus. In particular, LPPSubFus using neural network algorithms achieved the highest overall accuracy (95%) and over 90% producer and user accuracy for tea plantations and forests. In addition, LPPSubFus was more compatible with di ﬀ erent classiﬁcation algorithms than SubFus. Based on these ﬁndings, it is concluded that LPPSubFus has better and more stable performance in tea plantation mapping than pixel-level image fusion and SubFus. This study demonstrates the potential of integrating hyperspectral and SAR data via LPPSubFus for mapping tea plantations. Our work o ﬀ ers a promising tea plantation mapping method and contributes to the understanding of hyperspectral and SAR data fusion.


Introduction
The tea plant (Camellia sinensis (L.) O. Kuntze), an evergreen broadleaved perennial shrub, originated in China, and tea has gradually become one of the three most popular beverages in the world [1,2]. Tea is widely distributed globally. China, India, and Sri Lanka are the major tea-producing countries in Asia, and Kenya, Malawi, Rwanda, Tanzania, and Uganda are the major tea-producers in Africa [3]. Tea planting and production play a vital role in the agricultural economy and rural development of these developing countries. Driven by the considerable economic benefits of tea cultivation, the area of tea plantations is continually increasing [4]. According to data of the Food and Agriculture Organization (FAO) of the United Nations, the global area of tea plantations reached 4,193,176 hectares in 2018, and the tea plantation area increased by 77.33% between 2000 and 2018 [5]. In China, a global leader in tea cultivation and production, tea is mainly grown in the red soil hilly areas of the southern provinces, such as Fujian and Zhejiang. In 2018, China's tea plantations occupied an area of 2,347,726 hectares, increasing 161% since 2000 [5]. Although the expansion of tea

Field Survey
Field surveys were conducted in June and December 2019. The field survey results showed that, in addition to the widespread distribution of tea plantations, evergreen forests, croplands, waterbodies, built-up areas, and bareness rocks were distributed throughout the study area ( Figure  2). Trees were unevenly scattered in some tea plantations (Figure 2a-c). The winter and summer crops in the study area are shown in Figure 2d

Field Survey
Field surveys were conducted in June and December 2019. The field survey results showed that, in addition to the widespread distribution of tea plantations, evergreen forests, croplands, waterbodies, built-up areas, and bareness rocks were distributed throughout the study area ( Figure 2). Trees were unevenly scattered in some tea plantations (Figure 2a-c). The winter and summer crops in the study area are shown in Figure 2d

Remote Sensing and Auxiliary Data
The HS data used in this study was acquired by a Chinese Gaofen-5 satellite, which covers the spectral ranges from 400 to 2500 nm with 30 m spatial resolution. The SAR data was acquired by Sentinal-1A satellites that were operated at a frequency of 5405 GHz and collected in Interferometric Wide Swath mode, which acquires data with a 250 km swath width at 5 × 20 m spatial resolution, and vertical transmit and horizontal receive (VH) and vertical transmit and vertical receive (VV) polarization. These two datasets were collected in November because other crops in the study area had been harvested by this time, which helped distinguish the tea plantations. The detailed information of HS and SAR data used in this study is provided in Table 1. Auxiliary data included land cover/land use map of the study area and a 30 m digital elevation model (DEM).

Remote Sensing and Auxiliary Data
The HS data used in this study was acquired by a Chinese Gaofen-5 satellite, which covers the spectral ranges from 400 to 2500 nm with 30 m spatial resolution. The SAR data was acquired by Sentinal-1A satellites that were operated at a frequency of 5405 GHz and collected in Interferometric Wide Swath mode, which acquires data with a 250 km swath width at 5 × 20 m spatial resolution, and vertical transmit and horizontal receive (VH) and vertical transmit and vertical receive (VV) polarization. These two datasets were collected in November because other crops in the study area had been harvested by this time, which helped distinguish the tea plantations. The detailed information of HS and SAR data used in this study is provided in Table 1. Auxiliary data included land cover/land use map of the study area and a 30 m digital elevation model (DEM). bands were removed for the spectral overlap between the visible near-infrared and short wave infrared band. Subsequently, 301 bands were chosen for radiometric calibration, Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes atmospheric correction, Savitzky-Golay filtering, and orthographic correction. Regarding Savitzky-Golay filtering, the center point of the filtering core was 8, the derivative order was 0, and the degree of the polynomial was 4. Finally, we performed visual inspection and removed 47 bands due to stripe noise. Hence, 254 bands were outputted to conduct the following study.
Terrain correction was performed for the HS image via the C-correction model based on DEM data to reduce the terrain effect. C-correction is a useful method for normalizing topographic illumination effects [28]. This method introduces an empirically derived parameter to explain the overcorrection of the simpler cosine method in low illumination areas [29].
where L H is the corrected image value, L T is the uncorrected image value, sz is the solar zenith angle, i is the solar incident angle concerning the surface normal direction, and c is the C-correction coefficient derived by regressing L T against cos(i) and taking the quotient of the intercept and slope.

SAR Data Preprocessing
The Single Look Complex data of VV and VH polarization were first multilooked with a range number of two looks in the azimuth and eight in the range, resulting in a ground-range resolution of 29.52 × 27.87 m. The speckle noise was then removed using the refined Lee filter with a 5 × 5 window. The VV and HH images were then geocoded to the World Geodetic System 1984 datum and Universal Transverse Mercator Zone 50 North coordinate system via Range-Doppler terrain correction using the DEM data. Finally, VV and VH polarized data were resampled to a 30 m spatial resolution and georeferenced with the HS data using 20 ground control points with an error of less than 0.5 pixels. After data preprocessing, we obtained processed HS and SAR data with a size of 587 × 377 pixels and a spatial resolution of 30 m.

Feature Extraction
This study extracted four sets of features from HS and SAR data, namely, principal component bands, spectral indexes, polarimetric features, and texture information. These features were extracted based on the Environment for Visualizing Images (ENVI) software program.

Principal Component Analysis
Principal component analysis (PCA) is a commonly used dimensionality reduction method. PCA transforms a set of variables that may be correlated into fewer variables that are linearly independent using orthogonal transformation [30]. We performed PCA on the HS data. The first three principal components, which accounted for most of the data variance, were extracted to participate in the classification as spectral variables.

Spectral Index
Spectral characteristics of six land cover types are shown in Figure 3. First, six bands with wavelengths 792, 865, 1006, 1173, 1266, and 1662 nm were selected because of the big difference between tea plantation and forest in these bands. Second, four vegetation indexes that reflect the canopy structure, pigmentation, and water content of vegetation canopy were chosen to improve the spectral separability of tea (Table 2), namely, the Chlorophyll Index 2 [31], Anthocyanin Reflectance Index 2 [32], Vogelmann Red Edge Index 2 [33], and Water Band Index [34]. canopy structure, pigmentation, and water content of vegetation canopy were chosen to improve the spectral separability of tea (Table 2), namely, the Chlorophyll Index 2 [31], Anthocyanin Reflectance Index 2 [32], Vogelmann Red Edge Index 2 [33], and Water Band Index [34].  We further analyzed the spectral reflectance curves to fully use the rich spectral characteristics of HS data. We found that the variation tendencies of tea and forest spectral curves are similar, with the reflectance of tea higher than that of the forest, especially at wavelengths 800-1600 nm. For example, the spectral reflectance values of wavelengths 865 nm are different, whereas spectral reflectance values of 454 nm are the same, and thus the slope of the line between wavelengths 865 and 454 nm are different (the black line in Figure 3). Moreover, the triangle area enclosed by the bands with wavelengths 454, 1006, and 2032 nm is different between tea and forest (the blue polygon in Figure 3). Therefore, constructing new spectral indexes for the slope and area-dependent may improve tea and forest separability. We created three tea vegetation indexes (TVIs) based on the spectral slope between wavelengths 454 and 865 nm, 1266 and 2032 nm, and the triangle area enclosed by wavelengths 454, 1006, and 2032 nm. The equations of the slopes between wavelength 454 and 865, and 1266 and 2032 nm, are TVI 1 and TVI 2 , respectively: TVI 2 = (ρ2032 nm − ρ1266 nm)/∆λ 2 (3) where ρ is the spectral reflectance value of the corresponding band, ∆ is the difference in wavelength between 454 and 865 nm, and ∆ is the difference in wavelength between 1266 and 2032 nm.
The equation of the triangle area enclosed by 454, 1006, and 2032 nm is TVI 3 : where ∆ is the difference in wavelength between 454 and 1006 nm, ∆ is the difference in wavelength between 1006 and 2032 nm, and ∆ is the difference in wavelength between 454 and 2032 nm.
The results of the three TVIs are shown in Figure 4. The tea plantations showed a high value in TVI 1 and TVI 3 images, but a low value in TVI 2 images. The value of other land cover types on the

Name of Variable Description
Chlorophyll Index 2 ρ762 nm/ρ702 nm − 1 Anthocyanin Reflectance Index 2 ρ800 nm(1/ρ550 nm − 1/ρ700 nm) Vogelmann Red Edge Index 2 (ρ734 nm − ρ747 nm)/(ρ715 nm + ρ726 nm) Water Band Index ρ900 nm/ρ970 nm We further analyzed the spectral reflectance curves to fully use the rich spectral characteristics of HS data. We found that the variation tendencies of tea and forest spectral curves are similar, with the reflectance of tea higher than that of the forest, especially at wavelengths 800-1600 nm. For example, the spectral reflectance values of wavelengths 865 nm are different, whereas spectral reflectance values of 454 nm are the same, and thus the slope of the line between wavelengths 865 and 454 nm are different (the black line in Figure 3). Moreover, the triangle area enclosed by the bands with wavelengths 454, 1006, and 2032 nm is different between tea and forest (the blue polygon in Figure 3). Therefore, constructing new spectral indexes for the slope and area-dependent may improve tea and forest separability. We created three tea vegetation indexes (TVIs) based on the spectral slope between wavelengths 454 and 865 nm, 1266 and 2032 nm, and the triangle area enclosed by wavelengths 454, 1006, and 2032 nm. The equations of the slopes between wavelength 454 and 865, and 1266 and 2032 nm, are TVI 1 and TVI 2 , respectively: where ρ is the spectral reflectance value of the corresponding band, ∆λ 1 is the difference in wavelength between 454 and 865 nm, and ∆λ 2 is the difference in wavelength between 1266 and 2032 nm. The equation of the triangle area enclosed by 454, 1006, and 2032 nm is TVI 3 : where ∆ω 1 is the difference in wavelength between 454 and 1006 nm, ∆ω 2 is the difference in wavelength between 1006 and 2032 nm, and ∆ω 3 is the difference in wavelength between 454 and 2032 nm. The results of the three TVIs are shown in Figure 4. The tea plantations showed a high value in TVI 1 and TVI 3 images, but a low value in TVI 2 images. The value of other land cover types on the TVI images was opposite to that of the tea plantation. Particularly in TVI 1 and TVI 3 images, the difference between tea plantation and forest was evident (the red polygon in Figure 4). These three TVIs all play a role in increasing the difference between tea plantations and other land cover types, thus improving classification accuracy.
Forests 2020, 11, x FOR PEER REVIEW 7 of 22 TVI images was opposite to that of the tea plantation. Particularly in TVI 1 and TVI 3 images, the difference between tea plantation and forest was evident (the red polygon in Figure 4). These three TVIs all play a role in increasing the difference between tea plantations and other land cover types, thus improving classification accuracy.

Polarimetric Feature
SAR data have the potential for land cover classification because they are not affected by clouds and weather [10,[35][36][37]. In addition, SAR data can reflect the geometric and dielectric properties of scattering, improving classification accuracy [38,39]. In this study, four polarimetric features of Sentinel-1A data were extracted. They were the VH and VV intensity, the difference between VV and VH, and the VV and VH intensity ratio ( Figure 5).

Polarimetric Feature
SAR data have the potential for land cover classification because they are not affected by clouds and weather [10,[35][36][37]. In addition, SAR data can reflect the geometric and dielectric properties of scattering, improving classification accuracy [38,39]. In this study, four polarimetric features of Sentinel-1A data were extracted. They were the VH and VV intensity, the difference between VV and VH, and the VV and VH intensity ratio ( Figure 5).
Forests 2020, 11, x FOR PEER REVIEW 7 of 22 TVI images was opposite to that of the tea plantation. Particularly in TVI 1 and TVI 3 images, the difference between tea plantation and forest was evident (the red polygon in Figure 4). These three TVIs all play a role in increasing the difference between tea plantations and other land cover types, thus improving classification accuracy.

Polarimetric Feature
SAR data have the potential for land cover classification because they are not affected by clouds and weather [10,[35][36][37]. In addition, SAR data can reflect the geometric and dielectric properties of scattering, improving classification accuracy [38,39]. In this study, four polarimetric features of Sentinel-1A data were extracted. They were the VH and VV intensity, the difference between VV and VH, and the VV and VH intensity ratio ( Figure 5).

Spatial Feature
Image texture was defined as the spatial arrangement and variation frequency of the value levels of an adjacent set of pixels in the local neighborhood of an image [21]. Adding texture information contributes to classification accuracy when mapping tea plantations [7]. Hence, the Grey-Level Co-occurrence Matrix textures (contrast, dissimilarity, homogeneity, mean, variance, and correlation) [40] were then sequentially calculated via the first three principal components and four polarimetric features using a window of 3 × 3 with a direction of 45 • and a spatial distance of 1 pixel. Therefore, 18 and 24 features were extracted from the HS and the SAR data, respectively.

Image Fusion
The SubFus method proposed by Rasti and Ghamisi [27] is an innovative FLIF method to integrate diverse remote sensing data for land cover classification. This method proposes a novel subspace minimization with its derived solution to fuse features, making it applicable to various sensors with different characteristics. Moreover, SubFus can be easily extended to other applications in addition to remote sensing. The principle of SubFus is as follows: Let H and S be the data from HS and SAR sensors, respectively. The features obtained from HS and SAR data can be formulated as: where matrices H (n × p 1 ) and S (n × p 2 ) contain the vectorized features for HS and SAR data, respectively, in their i-th columns. Then, this method assumes that multi-source data can be modeled using the fused features (F) in different subspaces. The extracted features contain redundant information, and therefore the multi-source features should have low-rank representations. Accordingly, the features from H and S can be modeled in lower-dimensional space as: N 1 and N 2 are n × p matrices representing the noise and model error, F is an (n × r) matrix containing the fused features, and V 1 and V 2 are (p 1 × r and p 2 × r, respectively) subspace projection matrices. For the standard subspace rank r, there is a lower bound for the rank selection. This method suggests the minimum possible rank, i.e., r ≤ min(n, p 1 , p 2 ), to achieve the maximum reduction in the feature space.
Then the constrain penalized cost function to estimate the unknown matrices F simultaneously, V 1 , and V 2 as where: In Equation (8), λ 1 and λ 2 control the weight between the fidelity of SAR data and the smoothness of the subspace features, respectively, with respect to the fidelity of HS data. The constraints enforce the orthogonality on the subspaces. The alternating direction method of multipliers is used to solve Equation (8). The variable F is divided by: Then, the penalty method is used to add the constraint F = M to the minimization problem. Therefore, the augmented minimization problem can be represented as: where B is the Lagrangian multiplier. After the variables are split, the cyclic descent algorithm is adopted to solve one unknown matrix at a time, using five iterative steps. The detailed steps are not repeated here.
It should be noted that the right singular vectors of H and S are used to initialize V 1 and V 2 , respectively, and M 0 and B 0 are set to zero in SubFus. However, the singular eigenvector computed from the original data matrix may ignore the essence of the data, which may not be the best representation of the original data. Therefore, a more appropriate method should be explored to initialize V 1 and V 2 .

Combination of Locality Preserving Projection and SubFus
Locality preserving projection (LPP), a manifold learning algorithm, was introduced to improve the SubFus, and we called this new method LPPSubFus. The basic idea of LPP is to maintain the close and far affinity relationship between each sample pair in the projection and retain the local neighborhood structure of the samples in the space when reducing the dimensionality [41]. That is, LPP aims to optimize projection P which converts the original data X to a lower-dimensional representation Y = P T X [42]. The objective of this study is to obtain the projection P that is used to replace the singular vector to initialize V 1 and V 2 . The LPP introduced into the SubFus has two advantages: (1) the data dimensionality can be initially reduced before SubFus, and (2) LPP can better reflect the essence of data than the singular vector. The principle of solving projection P is as follows [43]: x m ] is set to the original data and lower-dimensional representation of Y = [y 1 , y 2 , . . . , y m ] is a m × d matrix. A reasonable criterion for deciding Y is to minimize the following objective function: The parameter W is obtained by the k-nearest-neighborhood and induces a punishment when the neighbor area x i is far away from x j . The formula of W is as follows: x i and x j are local neighbors 0 x i and x j are not local neighbors (12) where σ is a filtering parameter that can be determined empirically. Suppose p is a transformation vector, that is, Y = p T X, by simple algebraic formulation, the optimization is formulated as follows: where D is a diagonal matrix; its entries are the column (or row, because W is symmetric) sum of W.
Therefore, the constraint is imposed as follows: Subject to p T XDX T p = 1 The transformation vector p that minimizes the objective functions is given by the solution of the minimum eigenvalue to the generalized eigenvalues problem: It is easy to show that the matrices XLX T and XDX T are symmetric and positive semidefinite. Let the column vectors p 0 , p 1 , . . . p d−1 be the solutions of Equation (15), ordered according to their eigenvalues, λ 0 , λ 1 , . . . λ d−1 . Thus, the embedding is as follows: The projection vector P H and P S of HS and SAR features can be obtained, respectively. Then, they are used to initialize V 1 and V 2 in SubFus.
In this study, SubFus and LPPSubFus were implemented in the MATLAB software program. Regarding SubFus, λ 1 , λ 2 , and µ were set to 0.5, 0.1, and 10, respectively. Then, in LPPSubFus, k was set to 5, and σ was set to 5. The remaining parameter values were the same as those in SubFus.

Pixel-Level Image Fusion
Three PLIF methods were also used in this study to fuse HS and SAR data under the MATLAB software program. The PLIF methods mainly include component substitution, multi-scale decomposition, and the hybrid method [25]. Therefore, we chose PCA [23], discrete wavelet transforms (DWT) [22], and wavelet-principal component analysis (W-PCA) [21] to represent these three types. In this study, six bands of HS data were fused with VV polarized image because the co-polarization SAR data often performs better than cross-polarization in vegetation identification [44,45].

Classification Algorithm
The machine learning algorithms show good performance in remote sensing applications [24,[46][47][48]. This study selected two typical and popular machine learning algorithms, namely SVM and neural network (NN) algorithms, to carry out classifications.
The SVM algorithm aims to define an optimum hyperplane to separate two different classes using training data without any assumptions [49]. Although SVM has limitations in parameter selection and computational requirements, it provides superior performance compared to most other classification algorithms [50]. Many studies have used the SVM method to extract tea plantation distribution and achieved good results [8,9,51]. Moreover, Melgani, and Bruzzone [52] demonstrated that SVM has good performance in HS data classification because it is less sensitive to the Hughes phenomenon. The major parameters of SVM algorithms include kernel function, gamma value, and penalty parameter.
The NN algorithm is a computing model composed of essential nodes, including the input layer, hidden layer, and output layer [53]. This method is competent for parallel computing, automatic learning, and correcting errors [54], and has long been used for remote sensing applications, such as land cover classification [55], wetland vegetation mapping [56], and tree species identification [54]. Many have suggested that NN algorithms are superior to traditional statistical classification algorithms because they learn from training samples without assumptions about data distribution [48]. The major parameters of NN algorithms include training threshold contribution, training rate, training momentum, training root mean square error exit criteria, number of hidden layers, and training iteration times.

Classification Scenario and Sample
Six classification scenarios were developed to compare the performance of different image fusion methods in this study (Table 3). Scenario 1 used only the HS data, scenarios 2-4 used the data fused by PLFI, and scenarios 5 and 6 fused features of HS and SAR data using SubFus and LPPSubFus, respectively. In addition, all corresponding features in each classification scenario were normalized to weaken the impact of different attributes of features. Table 3. Classification scenarios using the support vector machine and neural network algorithms.

Scenario
Image Layer Used for Classification In this study, the built-up area and waterbody were classified into one category. Therefore, five categories needed to be classified. The training and test samples for the five categories were randomly chosen from the field survey results and land cover/land use map ( Table 4). The ratio of training samples to test samples was 1 to 2, and all samples were distributed evenly throughout the study area. The sample number of tea plantations and forests were two times as many as other land cover types, which ensured the accuracy assessment result mainly depended on tea plantations and forests.

Parameter Setting and Accuracy Assessment
The classifications were carried out using the ENVI software program. Regarding the SVM-based classifications, the radial basis function was selected as the kernel function because it is suitable for non-linear data [57]. The gamma value was dependent on the inverse of the number of input variables [9], and the default value (100) of the penalty parameter was used for all classifications. For all NN-based classifications, the logistic function was selected as the activation function, the training threshold contribution was set to 0.9, the different learning rate was set to 0.2, training root mean square error was set to 0.1, and momentum value was set to 0.9. The number of nodes in the hidden layer and the number of training iterations were set to 1 and 1000, respectively. After classifications, an accuracy assessment was also performed, including the traditional confusion matrix, Kappa coefficient, producer accuracy (PA), and user accuracy (UA).

Visual Comparison of Thematic Maps
The classifications of different scenarios using SVM and NN algorithms are shown in Figure 6, where scenarios 5 and 6 yielded the highest OA in SubFus-and LPPSubFus-based classifications (scenario 5: the feature number using SVM and NN algorithms was 27 and 12, respectively; scenario 6: the feature number using SVM and NN algorithms was 28 and 22, respectively). In both SVM-and NN-based classifications, scenario 1 overestimated the tea plantation classification results, and some forests were misclassified as tea plantations (Figure 6a,g). The distribution of tea plantations was reduced in the classifications of scenarios 2-4 ( Figure 6b-d,h-j). However, the classifications using scenarios 2-4 retained speckle noise, and these classifications were unsatisfactory. Scenarios 5 and 6 improved the classification with little speckle noise, consistent with the land cover distribution (Figure 6e,f,k,l). Moreover, scenario 6 more accurately classified forests in the central and southwestern parts of the study area than other scenarios.

Visual Comparison of Thematic Maps
The classifications of different scenarios using SVM and NN algorithms are shown in Figure 6, where scenarios 5 and 6 yielded the highest OA in SubFus-and LPPSubFus-based classifications (scenario 5: the feature number using SVM and NN algorithms was 27 and 12, respectively; scenario 6: the feature number using SVM and NN algorithms was 28 and 22, respectively). In both SVM-and NN-based classifications, scenario 1 overestimated the tea plantation classification results, and some forests were misclassified as tea plantations (Figure 6a,g). The distribution of tea plantations was reduced in the classifications of scenarios 2-4 ( Figure 6b-d,h-j). However, the classifications using scenarios 2-4 retained speckle noise, and these classifications were unsatisfactory. Scenarios 5 and 6 improved the classification with little speckle noise, consistent with the land cover distribution (Figure 6e,f,k,l). Moreover, scenario 6 more accurately classified forests in the central and southwestern parts of the study area than other scenarios.  Some rocks in the upper right part of the study area were incorrectly classified as waterbodies in the classifications of scenarios 2-4 due to topography shadow, whereas the rocks were correctly classified by scenarios 5-6. However, the waterbodies in the classification of scenario 6 using the NN algorithms were not well identified (Figure 6l). By comparison, because the primary cultivated land type in the study area is the paddy field, there was confusion between waterbodies and croplands in the classifications of some scenarios, such as in scenario 4 (Figure 6d,j).
The classifications using only HS data with SVM and NN algorithms had considerable speckle noise. Image fusion improved the noise phenomenon, especially in the FLIF-based classifications. In addition, the classifications using NN algorithms were better than those using SVM algorithms. For example, the distribution of tea plantations and forests was similar in the classifications of scenarios 5 and 6. However, some croplands were mistakenly identified as waterbodies and had evident shadows when using the SVM algorithms (Figure 6e,f,k,l). In terms of visual comparisons, all the scenarios using both SVM and NN algorithms could identify the land cover types well in the study area. Among these, the classifications of scenario 1 had considerable speckle noise, and classifications of scenarios 2-4 were similar. Scenarios 5 and 6 achieved better results than other scenarios.

Comparison of Overall Accuracy
The OA and Kappa coefficient of the six scenarios are summarized in Table 5. In SVM-based classifications, scenario 1 achieved an OA of 89.05%. There was an improvement of OA when using scenarios 2 and 4. However, the OA of the classifications with scenario 3 was 1.31% less than that with scenario 1. Scenario 5 achieved the same OA as scenario 1, and the Kappa coefficient was slightly lower than that of scenario 1. Scenario 6 achieved the highest OA of 92.14% in SVM-based classifications. In NN-based classifications, scenario 1 yielded a 92.14% OA, then the OA of scenarios 2-5 decreased by 1.31%-6.07% compared to scenario 1. Scenario 6 achieved the highest OA of 95% and the highest Kappa coefficient of 0.9355 in NN-based classifications, and the OA with scenario 6 was increased by 2.86% over that with scenario 1.  Figure 7 shows the OA curves of SubFus and LPPSubFus with the feature number. For SubFus-based classifications (Figure 7a), the OA using SVM and NN algorithms had a similar tendency. The OA gradually increased with the increase in the feature number, and finally was close to 90%. For LPPSubFus-based classifications (Figure 7b), the OA of SVM-based classifications changed smoothly and ultimately remained at about 90%. The NN-based classifications had high OA with significant variation later, notably when the feature number exceeded 20. Most of the OA values in NN-based classifications were higher than those in SVM-based classifications. In addition, LPPSubFus-based classifications yielded higher OA than SubFus-based classifications, particularly after the feature number was greater than 16. To summarize, the LPPSubFus achieved better classifications than SubFus, and NN algorithms performed better than SVM algorithms in this study.
and ultimately remained at about 90%. The NN-based classifications had high OA with significant variation later, notably when the feature number exceeded 20. Most of the OA values in NN-based classifications were higher than those in SVM-based classifications. In addition, LPPSubFus-based classifications yielded higher OA than SubFus-based classifications, particularly after the feature number was greater than 16. To summarize, the LPPSubFus achieved better classifications than SubFus, and NN algorithms performed better than SVM algorithms in this study.

Comparison of Producer and User Accuracy
This study aimed to improve the separability of tea plantations via HS and SAR data. The tea and forest in the study area are evergreen broadleaved perennial vegetation [9]. Only if the tea plantation and forest are identified well can accurate tea plantation mapping be achieved. The PA and UA of tea plantations and forests using two classification algorithms under different scenarios are shown in Table 6. Both SVM and NN algorithms achieved high PA of tea plantations, and low PA of forests and low UA of tea plantations, when using scenarios 2-4. The PA and UA of tea plantation and forest only in the classifications of scenario 6 were above 90%, indicating that the tea plantations and forests in this classification were identified well with fewer misclassifications and omissions. Generally, the PA and UA of scenario 6 with NN algorithms were higher than those with SVM algorithms.

Comparison of Producer and User Accuracy
This study aimed to improve the separability of tea plantations via HS and SAR data. The tea and forest in the study area are evergreen broadleaved perennial vegetation [9]. Only if the tea plantation and forest are identified well can accurate tea plantation mapping be achieved. The PA and UA of tea plantations and forests using two classification algorithms under different scenarios are shown in Table 6. Both SVM and NN algorithms achieved high PA of tea plantations, and low PA of forests and low UA of tea plantations, when using scenarios 2-4. The PA and UA of tea plantation and forest only in the classifications of scenario 6 were above 90%, indicating that the tea plantations and forests in this classification were identified well with fewer misclassifications and omissions. Generally, the PA and UA of scenario 6 with NN algorithms were higher than those with SVM algorithms. Table 6. Producer accuracy (PA) and user accuracy (UA) of tea plantations and forests with support vector machine (SVM) and neural network (NN) algorithms for different scenarios.  Figure 8 shows the PA and UA curves of tea plantations with the feature number using SubFus and LPPSubFus. In the SubFus-based classifications, after the feature number exceeded 18, the PA of classifications using SVM algorithms was slightly higher than that using NN algorithms, with the exception of when the feature number was 22 (Figure 8a,b). However, the UA of NN-based classifications at the later stage was slightly higher than that of SVM-based classifications. In the LPPSubFus-based classifications, there was no significant difference between the PA and UA with SVM and NN algorithms (Figure 8c,d). However, curves of NN-based classifications were sensitive to changes in the feature number with apparent volatility, particularly the UA of tea plantations in the LPPSubFus-based classifications. In general, LPPSubFus-based classifications achieved higher PA or UA than SubFus-based classifications. In particular, the PA of tea plantations using LPPSubFus was close to 100% later, and the UA was around 95% in most cases.
Forests 2020, 11, x FOR PEER REVIEW 15 of 22 Figure 8 shows the PA and UA curves of tea plantations with the feature number using SubFus and LPPSubFus. In the SubFus-based classifications, after the feature number exceeded 18, the PA of classifications using SVM algorithms was slightly higher than that using NN algorithms, with the exception of when the feature number was 22 (Figure 8a,b). However, the UA of NN-based classifications at the later stage was slightly higher than that of SVM-based classifications. In the LPPSubFus-based classifications, there was no significant difference between the PA and UA with SVM and NN algorithms (Figure 8c,d). However, curves of NN-based classifications were sensitive to changes in the feature number with apparent volatility, particularly the UA of tea plantations in the LPPSubFus-based classifications. In general, LPPSubFus-based classifications achieved higher PA or UA than SubFus-based classifications. In particular, the PA of tea plantations using LPPSubFus was close to 100% later, and the UA was around 95% in most cases.  Figure 9 shows the PA and UA curves of forests with the feature number using SubFus and LPPSubFus. The PA and UA of SubFus-based classifications were relatively stable. The PA remained below 90% later, and the UA was close to 100% (Figure 9a,b). Regarding LPPSubFus-based classifications, the UA curves changed steadily (Figure 9d). However, the PA volatility of SVM-and NN-based classifications was more noticeable, especially for the latter (Figure 9c). When the feature number was between 18 and 22, the PA in LPPSubFus-based classifications using NN algorithms approached 100%, and the UA gradually approached 100%. In addition, LPPSubFus identified forest better than SubFus in this study.  Figure 9 shows the PA and UA curves of forests with the feature number using SubFus and LPPSubFus. The PA and UA of SubFus-based classifications were relatively stable. The PA remained below 90% later, and the UA was close to 100% (Figure 9a,b). Regarding LPPSubFus-based classifications, the UA curves changed steadily (Figure 9d). However, the PA volatility of SVM-and NN-based classifications was more noticeable, especially for the latter (Figure 9c). When the feature number was between 18 and 22, the PA in LPPSubFus-based classifications using NN algorithms approached 100%, and the UA gradually approached 100%. In addition, LPPSubFus identified forest better than SubFus in this study.

Comparison of Pixel-and Feature-Level Image Fusion
The classifications obtained by FLIF were consistent with the real distribution of the tea plantations and forests. In addition, FLIF achieved higher OA than PLIF, with the exception of the classification of SubFus using SVM algorithms, and a 1.31%-8.93% improvement was achieved using FLIF compared to PLIF (Table 5). In particular, LPPSubFus performed well in both SVM-and NNbased classifications, achieving the highest classification accuracy of 95%, and good PA and UA of tea plantations and forests. The above results demonstrate the superior ability of FLIF in mapping tea plantations. In PLIF-based classifications, the classification accuracy using SVM algorithms was not significantly improved compared to the HS data. Even in NN-based classifications, the accuracy dropped considerably. This shows that fusing HS and SAR data cannot guarantee classification accuracy improvement. Therefore, it is crucial to select the appropriate image fusion method for tea plantation mapping. Furthermore, the low OA of scenarios 2-4 in NN-based classifications indicate that FLIF also has the advantage of data dimension reduction and is compatible with more classification algorithms compared to PLIF.
As expected, the FLIF-based classification offered a more generalized visual appearance presentation and a more contiguous depiction of land cover types with fewer speckles than the PLIFbased classification ( Figure 6). The FLIF-based classification was similar to that of the object-based classification. However, the performance of the object-based image analysis mainly depends on the

Comparison of Pixel-and Feature-Level Image Fusion
The classifications obtained by FLIF were consistent with the real distribution of the tea plantations and forests. In addition, FLIF achieved higher OA than PLIF, with the exception of the classification of SubFus using SVM algorithms, and a 1.31%-8.93% improvement was achieved using FLIF compared to PLIF (Table 5). In particular, LPPSubFus performed well in both SVM-and NN-based classifications, achieving the highest classification accuracy of 95%, and good PA and UA of tea plantations and forests. The above results demonstrate the superior ability of FLIF in mapping tea plantations. In PLIF-based classifications, the classification accuracy using SVM algorithms was not significantly improved compared to the HS data. Even in NN-based classifications, the accuracy dropped considerably. This shows that fusing HS and SAR data cannot guarantee classification accuracy improvement. Therefore, it is crucial to select the appropriate image fusion method for tea plantation mapping. Furthermore, the low OA of scenarios 2-4 in NN-based classifications indicate that FLIF also has the advantage of data dimension reduction and is compatible with more classification algorithms compared to PLIF.
As expected, the FLIF-based classification offered a more generalized visual appearance presentation and a more contiguous depiction of land cover types with fewer speckles than the PLIF-based classification ( Figure 6). The FLIF-based classification was similar to that of the object-based classification. However, the performance of the object-based image analysis mainly depends on the complex parameter selection of image segmentation [46]. FLIF can avoid the complex parameter decisions and has an equal effect as object-based image analysis. Therefore, compared with PLIF and object-based image analysis, FLIF is superior in tea plantation mapping.

Comparison of SubFus and LPPSubFus
Based on the SubFus algorithms proposed by Rasti and Ghamisi, this study introduced LPP into SubFus and proposed the LPPSubFus algorithms. The main difference between LPPSubFus and SubFus is that LPPSubFus uses the projection matrix obtained by LPP to initialize variables instead of using singular value decomposition. We designed scenarios 5 and 6 to compare their performance in tea plantation mapping.
First, both SubFus and LPPSubFus effectively solved the noise phenomenon that existed in PLIF-based classifications, and the distribution of land cover obtained was similar. Moreover, the LPPSubFus-based classifications better coincided with our knowledge of the study area than SubFus-based classifications. Second, the classification accuracy results showed that LPPSubFus yielded an improvement in OA of 3.09% and 4.17% compared to SubFus in both SVM-and NN-based classifications (Table 5), respectively. Furthermore, LPPSubFus provided better classifications for tea plantation and forest than SubFus, with both PA and UA more than 90% (Table 6). Third, compared with HS data, LPPSubFus improved the OA in SVM and NN-based classifications, but SubFus did not. The above results indicate that LPPSubFus had outstanding and more stable performance than SubFus in tea plantation mapping. However, LPPSubFus also has shortcomings, such as not identifying waterbodies well in NN-based classifications.

Benefits of Image Fusion for Tea Plantation Mapping
Due to the development of satellite sensing technology, multi-source remote sensing data applications have become more prevalent [58][59][60][61]. Data from different sensors can provide extra information that makes up for the deficiency of single remote sensing data [25]. In addition, multi-source data are more easily accessible than high-quality time-series data because weather easily influences the optical sensors [35,62]. Among the many different sensors, HS and SAR data have the advantage of low cost and easy access.
The results show that HS data fused SAR data by FLIF improved OA by up to 3.09% compared to using HS data alone (Table 5). However, the fusion of HS and SAR data by PLIF led to decreased classification accuracy generally. SAR data can provide backscattering information in surface roughness, water content, and shape [21]. However, SAR data are easily affected by speckle noise due to the entirely different optical data imaging mechanism, which negatively impacts information extraction [24]. Only by choosing an appropriate image fusion method can useful information be effectively extracted. We found that not all image fusion methods improved classification accuracy for tea plantation mapping. Compared with optical data alone, it is difficult to judge whether the fusion of optical and SAR data can achieve better results in this study. Fortunately, LPPSubFus had stable performance in fusing HS and SAR data in this study, and the classification accuracy was improved under both SVM and NN algorithms. Therefore, it is beneficial to fuse HS and SAR data via LPPSubFus for tea plantation mapping, in terms of both data acquisition and classification accuracy.

Classification Algorithm for Tea Plantation Mapping
This study used two classification algorithms to compare the performance of different fusion methods in tea plantation mapping. The results show that both SVM and NN algorithms achieved good results with the highest OA of more than 90%. The OA in NN-based classifications using scenario 1 was the same as the highest OA in SVM-based classifications. Moreover, scenario 6 with the NN algorithms achieved the highest OA of 95%, and UA and PA above 90% for tea plantations and forests. These results show that NN algorithms had better performance than SVM in this study.
For scenarios 1, 5, and 6, the OA of the NN-based classifications was higher than that of SVM-based classifications, whereas the OA of scenarios 2-4 was not. We speculated that this might be due to too many features in the classifications of scenarios 2-4, and the NN algorithms had insufficient ability to process high-dimensional data. In addition, NN algorithms were sensitive to the feature number, and the OA volatility in NN-based classifications was more evident than that in SVM-based classifications (Figure 7). SVM algorithms are not susceptible to the "dimensionality disaster" and have a more robust ability to process a large number of features [63]. Therefore, the OA did not change significantly with the feature number. The results indicate the importance of data dimensionality reduction in classification using HS and SAR data. For some classification algorithms, too many features will reduce the data processing speed and classification accuracy. LPPSubFus achieved the highest accuracy in both SVM-and NN-based classifications, thus, it has advantages compared to other fusion methods in terms of classification accuracy and compatibility with different algorithms.

Conclusions
The accurate mapping of tea plantations ensures appropriate management of tea-producing region for environmental protection. We proposed a LPPSubFus method to fuse spectral indexes, backscattering, and texture information of HS and SAR data. This study then explored the performance of SubFus, LPPSubFus, and PLIF in tea plantation mapping with SVM and NN algorithms. Through visual comparison and accuracy assessment, the following conclusions can be drawn: (1) FLIF can reduce data dimensionality and significantly improve noise phenomenon and classification accuracy. However, PLIF did not contribute to classification accuracy and even led to decreased accuracy in some cases. (2) LPPSubFus improved OA by 3% compared to SubFus. In particular, LPPSubFus with NN algorithms achieved the highest OA of 95%, and the PA and UA for tea plantations and forests exceeded 90%, which was concluded to be the best method for tea plantation mapping in this study. (3) LPPSubFus is compatible with different classification algorithms and had more stable and superior tea plantation mapping performance than PLIF and SubFus.
In summary, two advantages exist for LPPSubFus for tea plantation mapping. First, LPPSubFus provides a higher classification accuracy than PLIF and has the ability to reduce data dimensionality. Second, compared to SubFus, LPPSubFus has better performance in both SVM and NN algorithms. This study provides a new approach for future tea plantation mapping. The LPPSubFus approach used in this study could also be extended to the use of HS and SAR data in other forest mapping applications.