Multi-Temporal Landsat-8 Images for Retrieval and Broad Scale Mapping of Soil Copper Concentration Using Empirical Models

Abstract

Mapping soil heavy metal concentration using machine learning models based on readily available satellite remote sensing images is highly desirable. Accurate mapping relies on appropriate data, feature extraction, and model selection. To this end, a data processing pipeline for soil copper (Cu) concentration estimation has been designed. First, instead of using single Landsat scenes, the utilization of multiple Landsat scenes of the same location over time is considered. Second, to generate a preferred feature set as input to a regression model, a number of feature extraction methods are motivated and compared. Third, to find a preferred regression model, a variety of approaches are implemented and compared for accuracy. In this research, 11 Landsat-8 images from 2013 to 2017 of Gulin County, Sichuan China, and 138 soil samples with lab-measured Cu concentrations collected from the area in 2015 are used. A variety a metrics under cross-validation are used for comparison. The results indicate that multi-temporal images increase accuracy compared to single Landsat images. The preferred feature extraction varies based on the regression model used; however, the best results are obtained using support vector regression and the original data. The final soil Cu map generated using the recommended data processing pipeline shows a consistent spatial pattern with a ground-truth land cover classification map. These results indicate that machine learning has the ability to perform large-scale soil heavy metal mapping from widely available satellite remote sensing images.

1. Introduction

Heavy metal concentration (HMC) is a crucial property of soil quality that affects environmental conditions and human health [1,2]. Traditionally, HMC maps are obtained by first conducting field sampling and then performing spatial interpolation. Nevertheless, the accuracy of the HMC maps obtained with the geostatistical interpolation approaches are usually limited with the spatial stationary assumption [3] and the sparse, expensive and time-consuming spatial samples process [4].

Remote sensing techniques offer efficiency, especially for studies using hyperspectral imagery to deduce soil HMC [5,6,7]. Due to the significant uncertainties in the measuring system, there are no known physical models available for the direct inversion of HMC from hyperspectral data. Nevertheless, using empirical models, researchers have had more success in estimating HMC based on hyperspectral data [7]. Currently, for HMC retrieval, most studies adopt point-based spectrometer measurements with the limitations of spatial coverage. Although some researchers employed hyperspectral images with dense and wide spatial coverage capability for mineral concentration estimation [8,9], it is challenging to obtain the hyperspectral imagery due to the lack of hyperspectral satellites and the cost of airborne systems. In contrast, multispectral satellite images are more readily available to researchers. For example, Landsat-8 imagery is freely available with the advantages of adequate spatial resolution and coverage, short revisit period, and broad spectral range over the visible, near-infrared, short wave infrared and thermal infrared portion of the spectrum. Although it is of benefit to utilize Landsat imagery for broad scale mapping HMC, such investigations are not common. Therefore, to develop a fast and efficient approach for soil HMC retrieval from Landsat-8 multispectral images is important.

The objective of this research is to investigate the use of Landsat-8 images for the estimation of soil copper (Cu) concentration. This objective is challenging due to the weak correlation between soil Cu concentration and spectral observations [10]; however, there is potential for this weak spectral response to support soil Cu retrieval. For example, wavelengths around 460, 1400, 1900 and 2200 nm have been found to be sensitive to Cu concentration in a mining area [11]. Research with respect to a preferred band number and the best spectral resolution for Cu detection have also been investigated [12]. Hence, developing a framework that can effectively enhance this weak correlation to boost retrieval accuracy is expected to have the ability to accurately map Cu concentrations. In this paper, we highlight the following key items to improve Cu concentration detection using Landsat imagery and follow with a detailed discussion of each one.

• Instead of just using a single remote sensing image (SI), multi-temporal images (MTIs) of the same area are considered;
• Instead of just using the original Landsat data, various feature extraction scenarios are evaluated;
• A number of regression models are compared and contrasted to improve overall accuracy;
• The complete data processing pipeline (multi-temporal data, feature extraction and model selection) can produce maps consistent with ancillary ground-truth information.

First, to overcome the limited spectral resolution of Landsat-8, we adopt the use of Landsat-8 MTI, as opposed to just single images, to enhance the spectral information. Landsat-8 imagery plays an important role in studies on land cover classification [13], land surface temperature and climate change [14,15], agriculture [16], and vegetation property retrieval [17,18,19]. More specifically, Landsat-8 has been used for studies on the retrieval of soil properties, e.g., soil moisture [20,21] and soil salinity [22], while there are few applications of Landsat-8 data for soil HMC retrieval due to their limited spectral sensitivity for soil heavy metal components. Landsat-8 data have a lower spectral resolution compared to hyperspectral data, so to improve detection of Cu concentrations, we combine MTIs for HMC retrieval.

The use of all bands of multi-temporal images has been applied to problems such as land cover classification [23] and unmixing [24] using multispectral data as well as quantitative retrieval of land surface properties from synthetic aperture radar (SAR) data [25,26]. These studies leverage MTI with repeated observations based on the assumption of the temporal correlation effect of the spectral-temporal profile. As an example, for crop classification, the crop categories are supposed to keep constant or highly correlated during the multi-temporal observation period.

We know of no published research that demonstrates the use of MTI for the retrieval of soil Cu concentration. The research reported in this paper involves the investigation of the benefits of using MTIs from 2013 to 2017 over using a single time observation in 2015 for retrieving the Cu concentration in 2015. Considering that the HMC has a strong impact on soil and vegetation [8], the assumption is that the Cu concentration in 2015 is correlated with the temporal observations and that pixels with similar Cu concentrations tend to have similar temporal observation patterns. Moreover, compared to SI, the MTI is captured under different environmental and illumination conditions, and therefore, using MTIs can statistically reduce the influence of some less relevant factors, e.g., the soil moisture content.

Second, to extract the Cu-relevant information in the high-dimensional temporal spectral cube, we explore the use of several common feature extraction techniques. The process of feature extraction projects features in the original space into a new feature space with lower dimensionality, hopefully finding more meaningful information to improve the learning performance and build better generalizable models [27]. The application of feature selection algorithms can produce better prediction accuracy [28,29,30,31]. Although the importance of feature selection for HMC retrieval from imagery has been realized, feature importance is often not evaluated in a systematic way. Furthermore, the relationship between the feature importance and different models is not considered and discussed. In our research, we systematically evaluate three commonly used feature extraction methods (principal component analysis (PCA), minimum noise fraction transform (MNF) and isometric feature mapping (ISOMAP)) and choose the preferred feature combination based on the importance of different features relative to different regression algorithms.

Third, to identify the preferred regression algorithm, we adopt and compare several common regressors: support vector regression (SVR), partial least square regression (PLSR) and artificial neural network (ANN). Due to a lack of knowledge of the underlying physics in the soil Cu retrieval problem, uncertainty from the false understanding of the mathematical relationship between features and observations can harm the inversion accuracy. Simply using any empirical model with the improper complexity or structure is not appropriate. Therefore, to select a regression model from a group of empirical models based on unbiased model evaluation is important. However, the absence of external cross-validation for unbiased accuracy measure and model selection is common [29,32,33,34,35,36] and raises uncertainties in retrieval results. Since test sets do not influence the model training [32], accuracy measures without external cross-validation is biased. Therefore, it is necessary to perform test set validation. To address this, we systematically evaluate and select regression models based on cross-validation conducted on the whole sample set [37].

Finally, we propose a data processing framework that involves the efficient data selection, feature extraction, and regression model. We use 11 multi-temporal Landsat-8 image scenes and 138 soil samples with Cu concentration over the Gulin County in Sichuan, China. Using this framework, we can explore the interaction between these three key components and identify a data processing pipeline that can increase the accuracy of Cu concentration estimation. The remainder of the paper is organized as follows. Section 2 demonstrates the theoretical background of using empirical models to retrieve soil Cu concentration from the Landsat data. Section 3 describes the study area, data and methods used for the retrieval. Section 4 shows the results, including the accuracies for different models and the pixel-based map of the Cu concentration distribution over the study area. Section 5 discusses the results and uncertainties. Section 6 summarizes the paper and draws conclusions.

2. Theoretical Background

The retrieval of soil Cu concentration from remote sensing images is based on the principle that Cu concentration can be represented as a function of the observed spectral bands, which, however, is very complex and unknown. The objective of this work is to establish an empirical model to retrieve the soil Cu concentration from Landsat-8 MTIs, which can be formulated as:

$$y=f_k\left(g_j\left(x_t\right)\right),$$

(1)

where $x_t$ is the Landsat-8 images observed at time point t, $g_j(·)$ is a function of data transformation for the feature extraction and selection and $f_k(·)$ is the regression model that maps the feature representation of $x_t$ to y, which is associated with the soil Cu concentration. The subscripts t, j, k of x, $g(·)$, $f(·)$ indicate, respectively, the variations of choices over data preparation, feature configuration and regressor selection, which are key issues in improving the performance of the empirical model. Therefore, to achieve an enhanced empirical model for Cu concentration estimation, this paper tries to investigate these three key issues by (1) using multi-temporal satellite observations temporally close to 2015 to enhance the spectral information of images, (2) extracting meaningful features using different approaches, i.e., PCA, MNF and ISOMAP, and selecting the most effective feature subset as model input and (3) employing multiple regressors, i.e., SVR, PLSR and ANN, to estimate the Cu concentration and choose the best regressor for mapping the Cu distribution over the study area. This research framework allows investigating the interactions among the three key components in the data inversion system and thereby enables the identification of the preferred combination of algorithms to boost the retrieval accuracy of Cu concentration.

3. Materials and Methods

We design an experiment to retrieve soil Cu concentration accurately and effectively by exploring key issues in empirical data inversion, i.e., data preparation, feature extraction and model selection. We compare three regressors that estimate the soil Cu concentration based on four features generated from eleven multi-temporal Landsat-8 images. We employ the cross-validation technique for bias-reduced accuracy measures. The Cu concentration map is generated using the preferred model with its preferred input features. Figure 1 shows the workflow.

3.1. Study Area and Soil Samples

The study area is located in Shiping village, Gulin County, Sichuan Province, China. The area has a subtropical climate, with annual average rainfall 748.4–1184.2 mm [38]. Due to many mining and industrial areas, the area is known to contain heavy metal contamination. The location of the study area and the 3D representation of the remote sensing image are shown in Figure 2. The soil in this area is mainly comprised of mollic gleysols, gleyic luvisols, planosols and eutric fluvisols according to the soil map databased provided by the Food and Agriculture Organization of the United Nations. A total of 138 soil samples were collected in this study area in 2015, and the Cu concentrations of samples were analyzed and recorded.

3.2. Multi-Temporal Images of Landsat-8

In this paper, 11 single date images of Landsat-8 ranging from 2013 to 2017 (path 128, row 40 and 41) are obtained from Geospatial Data Cloud website [39] and selected as the multi-temporal images with the cloudiness coverage lower than $10\%$. The imagery acquisition dates are summarized in

Table 1. Landsat-8 image acquisition dates.

2013	2014	2015	2016	2017
16 June 2013	6 August 2014 a	3 April 2015 a	8 June 2016	19 February 2017
	6 August 2014 b	3 April 2015 b	26 July 2015
	9 October 2014	8 July 2015
	28 December 2014

“a” represents for the image from row 40 of the satellite; “b” represents for the image from row 41 of the satellite.

. The downloaded images are level 1T products after topographic correction and have a 30 m spatial resolution except for the panchromatic band (band 8) with a spatial resolution of 15 m, which is not used in this study. Bands 10 and 11 of the 28 December 2014 imagery were destroyed due to a broken file and had to be excluded. As a result, the combined MTI comprises 108 bands in total.

3.3. Feature Selection

3.3.1. Feature Extraction

Feature extraction projects variables in the original space into a new feature space with a lower dimensionality, finding more meaningful information to improve the model learning performance and build better generalizable models [27]. In this study, three representative feature extraction methods are employed for extracting features from multi-temporal images of Landsat-8. Principal component analysis (PCA), minimum noise fraction (MNF) and isometric feature mapping (ISOMAP) are applied for feature extraction from multi-temporal images. The obtained features extracted by these three methods have 5, 5 and 10 dimensions separately.

3.3.2. Feature Importance Measure

The permutation feature importance measure used to evaluate different groups of features. The rationale is that when a variable is permuted, the connection between this variable and the response is lost. If this variable is associated with the response, using the permuted variable and the remaining unpermuted variables for prediction would lead to the prediction accuracy decreasing. Therefore, by calculating the difference in the prediction accuracy before and after permuting a variable, the importance of this variable can be quantified [40,41].

We explore the importance of four groups of features, i.e., original bands (OB) of MTIs, PCA, MNF and ISOMAP, by removing one of them each time and calculating the difference value of prediction accuracy before and after removing this group of features. In this manner, we can obtain the importance of each group of features related to different models used for predicting different kinds of HMC. Then, we select the “best” feature combination using forward feature selection [42], where the feature groups are added one by one from the most important to the least important until the prediction accuracy not increasing any more.

3.4. Regression Methods

Three regression methods are compared for the retrieval soil Cu concentration from Landsat-8 images. PLSR, ANN and SVM are all popular and efficient regression methods that are commonly used in soil property studies [8,18,43]. PLSR is based on the linear model, ANN represents the data-driven machine learning method and SVR implements regression by building a decision boundary.

3.4.1. Partial Least Square Regression, PLSR

PLSR is a specific form of multivariate linear regression [10], which has been widely used for soil property estimation [22,44,45]. PLSR assumes that the dependent variable can be estimated via a linear combination of explanatory variables [10]. By using the partial least-squares approach, PLSR obtains the projections of both dependent variables and explanatory variables, maximizing the covariance between projections, and then establishes the regression model based on these projections [46]. PLSR has the advantages of effectively dealing with collinear and noisy independent variables, modelling multiple response variables simultaneously [44,47] and having the inferential capability for possible linear relationships [44].

3.4.2. Artificial Neural Network, ANN

ANNs have become a popular machine learning regression algorithm over the last decade [46]. Compared to the linear regression models, ANNs can fit the data and model the complicated underlying relationships more flexibly [48]. Though ANN is good at processing nonlinear and complex problems, even with the inaccurate and noisy data [49], its performance highly relies on the network structure with a proper number of layers and neurons [46].

We use a multiple-layer feed-forward backpropagation network with three layers, including an input layer, a hidden layer and an output layer, with a learning rule based on the Levenberg–Marquardt algorithm. This common structure was adopted in related studies [44,50,51] and turned out to be effective for soil property retrieval. Tan-sigmoid and linear activation functions are applied to the hidden layer and the output layer separately in this study. A common hyperparameter optimization method, grid search [52], is adopted to find the optimal neuron number of the hidden layer in the set of $\{10, 15, 20, 25\}$.

3.4.3. Support Vector Regression, SVR

SVR is based on the same principles as the support vector machine for classification by minimizing error and individualizing the hyperplane, which maximizes the margin. As a kernel-based regression method, SVR solves nonlinear regression problems by transferring the data from the original space to a higher dimensional space via a kernel function [53]. In $ϵ$-SVR, the goal is to find an estimation function as flat as possible that has a deviation under $ϵ$ from the targets of the training data [54]. In this study, $ϵ$-SVR with a radial basis kernel is adopted to build the HMC inverse model from Landsat-8 imagery. The performance of SVR deeply depends on the hyperparameters [55]. There are two main hyperparameters to set, i.e., cost and $ϵ$. The cost is the penalty associated with errors larger than $ϵ$, and $ϵ$ represents the minimally required precision. Optimal hyperparameters are set by performing the grid search in a discretized two-dimensional parameter space.

3.5. Accuracy Measure

To quantify model performances based on different regression methods and different input features, four parameters, including the adjusted coefficient of determination (R${}^2$), root mean square error (RMSE), mean absolute error (MAE) and standard error (SE), are calculated. R${}^2$ indicates the goodness of fit of the model, ranging from 0 to 1. RMSE indicates the absolute estimation errors [44]. The total MAE is attributed by each error in proportion to its magnitude rather than its square [56]. SE is calculated by taking the standard deviation and dividing it by the square root of the sample size.

3.6. Cross-Validation Estimation

To guarantee an unbiased evaluation of models, the test dataset should be separated from the training dataset and follow the same distribution as the training set [57]. The accuracy measures evaluated should be based on both the training set and the test set to avoid overestimating the generalization ability of complex models due to the overfitting effect [41]. In the k-fold cross-validation, samples are randomly partitioned to k equal-sized subsets, and only one subset is used as the test set. The process is repeated k times so that each of the subsets is used as a test set once. In this study, 6-fold cross-validation is used for model evaluation. Therefore, given 138 soil samples, 115 samples are used to train the model and the remaining 23 samples are used to evaluate the model performance on the test dataset. This is repeated 6 times. In order to obtain results independent of a particular partitioning, this procedure is repeated 20 times, and then $6\times 20=120$ models are obtained. Both training and test accuracies of each model are calculated.

For training each model, 5-fold cross-validation is adopted to set hyperparamters that guarantee the best performance of each method, including the number of principal components for PLSR, cost and $ϵ$ value for SVR, as well as neuron numbers and learning rate for ANN. Hyperparameters that generate the highest R${}^2$ on test data would be selected for model training.

4. Results

4.1. Comparison between Single and Multi-Temporal Images of Landsat-8

Both the combined MTIs and SIs of Landsat-8 are used to train the regressors (i.e., PLSR, ANN and SVR) for the prediction of the soil Cu concentration. The 20 repeated 6-fold cross-validation is adopted to evaluate the performance of each model using all indicators: R${}^2$, RMSE, MAE and SE. The mean value of R${}^2$, as well as the standard deviation of R${}^2$, are shown in

Table 2. The mean and standard deviation values of regression indicators on three regressors over 20 repeated 6-fold cross-validation using MTL and SI on the training dataset and the test dataset.

		PLSR		SVR		ANN
		Mean	Std.dev	Mean	Std.dev	Mean	Std.dev
Adjusted R${}^2$	MTI	0.568	0.131	0.641	0.160	0.476	0.197
	SI	0.368	0.148	0.433	0.237	0.249	0.185
RMSE	MTI	16.997	3.178	15.515	3.426	20.651	5.664
	SI	16.257	3.014	14.636	3.508	23.305	6.213
MAE	MTI	15.215	3.021	12.498	3.914	17.673	5.444
	SI	19.493	2.752	16.711	4.857	22.994	5.310
SE	MTI	15.215	3.021	9.099	2.691	12.910	4.270
	SI	19.493	2.752	11.874	3.915	14.504	4.461

The highest accuracies are in bold format.

. The mean values of R${}^2$ on MTIs are higher than that on SI for all regressors, and the standard deviation of R${}^2$ decreases when using MTIs rather than SIs. For MAE and SE, the mean values of MTI are lower than those on SI, and the standard deviation is slightly higher than SI. The mean values of RMSE of PLSR and SVR do not show an obvious advantage of MTIs, given that the mean values increase when using MTIs. For ANN, RMSE reaches a lower RMSE mean and a lower standard deviation.

4.2. Feature Evaluation and Selection

Following the method of feature importance measurement introduced in Section 3.3.2, we first combine all the features to C1 and drop PCA, MNF and ISOMAP one by one to form C2, C3, C4 and C5 separately. The prediction accuracy of these five feature combinations on PLSR, SVR and ANN are summarized in

Table 3. Mean and standard deviation of prediction accuracies of Cu concentration estimated by PLSR, SVR and ANN using 20 repeated 6-fold cross-validations.

		PLSR	SVR	ANN	PLSR	SVR	ANN
		Mean			Std.dev.
Adjusted R${}^2$	C1	0.608	0.600	0.507	0.134	0.169	0.183
	C2	0.103	0.254	0.144	0.119	0.233	0.154
	C3	0.556	0.626	0.498	0.117	0.143	0.206
	C4	0.614	0.593	0.528	0.109	0.164	0.186
	C5	0.601	0.577	0.498	0.132	0.158	0.186
RMSE	C1	16.384	14.741	20.338	3.543	3.806	5.407
	C2	14.809	20.707	27.151	4.858	4.910	8.275
	C3	17.920	15.314	20.522	4.569	3.519	6.533
	C4	16.365	15.335	19.598	3.571	3.385	5.628
	C5	16.532	15.515	21.511	2.952	3.426	5.741
MAE	C1	14.475	13.698	16.717	2.955	4.140	4.722
	C2	24.426	20.835	26.282	3.854	5.491	5.900
	C3	15.872	13.193	17.028	2.961	3.675	5.576
	C4	14.246	13.702	16.423	2.786	3.865	4.732
	C5	14.691	14.340	17.378	2.738	3.636	4.957
SE	C1	14.475	9.189	12.595	2.955	2.831	3.931
	C2	24.426	12.995	16.786	3.854	3.920	5.942
	C3	15.872	9.552	12.931	2.961	2.884	5.040
	C4	14.246	9.481	12.125	2.786	2.624	3.852
	C5	14.691	9.579	13.118	2.738	2.594	3.826

The highest accuracies are in bold format.

. By subtracting the adjusted R${}^2$ mean values of C2, C3, C4 and C5 from that of C1 separately, we obtain the feature importance of each feature group, which is plotted in Figure 3. The importance score of these features based on RMSE, MAE and SE are all calculated in addition to the importance score based on R${}^2$. Except for the RMSE-based importance order of features, other indicators demonstrate the same importance rank. For all three regressors, OB is always the most important feature group for the retrieval of Cu concentration. PCA is the second important feature group for both PLSR and ANN, while SVR takes PCA as the worst input feature. The importance orders of three regressors are listed in

Table 4. Features importance in descreasing order for PLSR, SVR and ANN.

Importance Order	PLSR	SVR	ANN
1	OB	OB	OB
2	PCA	ISOMAP	PCA
3	ISOMAP	MNF	ISOMAP
4	MNF	PCA	MNF

.

Once the importance order of feature groups for each regressor was obtained, we use the forward feature selection strategy to determine the best feature combination by adding a feature once from the most important to the least important to the model input. For example, for PLSR, we use the OB as the input to train a regressor and test its accuracy. Next, OB together with the second important feature, PCA, is used, then ISOMAP and MNF are added one by one. For all three regressors, we tested the Cu prediction accuracy on different numbers of the feature. The mean values of test R${}^2$ are plotted in Figure 4.

For PLSR and ANN, the overall trend is that R${}^2$ increases first and then decreases as the feature number increases. For SVR, the R${}^2$ has reached the highest value when using OB only and then goes down with the increasing number of features. By identifying the peak positions of three regressors, the best feature combination for each regressor is identified.

4.3. Model Selection

To fairly compare three regressors—PLSR, SVR and ANN—we choose the most preferred feature combination for each regressor to ensure that we are comparing the best performances of these three regressors. We use OB as the input for SVR training, the combination of OB and PCA for PLSR and the combination of OB, PCA and ISOMAP for ANN. Moreover, the hyperparameters of the three models are determined by the grid search approach. Then, the prediction accuracies of three regressors are plotted in Figure 5, including the mean and standard deviation values of RMSE, MAE and ANN, since the R${}^2$ statistics have been displayed in Figure 4. SVR obtains the lowest mean values of RMSE, MAE and SE. PLSR performs worse than SVR and better than ANN.

4.4. Soil Cu Concentration Mapping

Given the operations of data preparation, feature extraction and model comparison, the OB feature is selected to input into the trained SVR model with the highest prediction R${}^2$, which is the “best” model among 120 models obtained during the 20 repeated 6-fold cross-validations. The adjusted R${}^2$ of this model reaches 0.87 on the training set and 0.84 on the test set (see Figure 6). By inputting the OB of the images over the entire study area to the preferred SVR model, the Cu concentration distribution map is generated, which is shown in Figure 7a.

To compare the retrieval and mapping of the Cu concentration to the traditional mapping approach based on interpolation, an interpolation map, shown in Figure 7b, is produced using exponential-kriging interpolation [58] on the ArcGIS platform. The exponential semi-variogram model is used as the semi-variance function, where the parameters, including the range, sill and nugget, are default values calculated internally. We set the number of points to 60 within the default search radius. Figure 7a shows more details about the Cu concentration distribution and has a larger spatial coverage compared with Figure 7b.

5. Discussion

5.1. Benefits of Using Multi-Temporal Landsat-8 Images

The experiment results of the comparison between MTIs and SIs (see Section 4.1) show that using MTIs to establish the inversion model of Cu concentration achieves better prediction precision, which is consistent with our assumption that MTIs can enhance the spectral information for a better estimation of Cu concentration. MTI techniques have been adopted to many other applications, e.g., crop classification [59,60] and spectral unmixing [24]. MTIs turned out to be more efficient for classification [59,60], unmixing [24] and land surface property retrieval [25,26] by many studies. Previouse research took advantage of the temporal correlation effect among the MTI observations. For example, for crop type identification, the assumption is that the crop type dictates the temporal patterns in the MTI data, i.e., different crop types tend to have different temporal behaviours. Similarly, we assume the Cu concentration in 2015 is correlated with the MTI observations from the year 2013 to the year 2017, which also generates better inversion accuracy than just using SI.

5.2. Feature Selection

Feature selection is conducted by first ranking four feature groups according to the permutation feature importance and then choosing the useful feature groups via the forward feature selection. As shown in Table 4, all three regressors take OB as the most valuable feature, given that without OB ranking first for all regressors, the estimation accuracy generated by the combination of other remaining features largely decreases. We attribute this importance of OB to its rich spectral information, which can be regarded as multiple repeated observations of the study area. Except for OB, the ranking the feature importance varies with different regressors, which is ascribed to different feature preferences of regressors. For example, PCA is the second most important feature group for both PLSR and ANN but is the least important for SVR (see Figure 3). A related work [41] reports the difference in feature preference between different classifiers on classification tasks.

This paper evaluates one feature importance under the existence of all other features. The importance of a certain feature to a regressor cannot be defined independently without taking into account other input features because its importance can vary with different feature combinations that input into the regressor with this feature. For example, the PCA feature group per se is very helpful for SVR but is less valuable with the existence of OB, MNF and ISOMAP. On the contrary, a feature can be more useful when with other features than itself [42].

5.3. Regressor Comparison

Based on four indicators obtained with the 20 repeated 6-fold cross-validations, SVR turns out to be the model with the best capability of generalization, given that the mean R${}^2$ of prediction is 0.641 and the maximum R${}^2$ of prediction is 0.844, as well as the smallest RMSE, MAE and SE among three regressors. PLSR also performs well with the mean prediction R${}^2$ 0.618 and much better than ANN with the mean prediction R${}^2$ 0.528. Generally speaking, SVR outperforms PLSR, and PLSR outperforms ANN at any feature group number (see Figure 4). In a related study [51], Cu estimation by ANN is also inferior to multiple linear regression. SVR has been demonstrated to have a stronger generalization ability compared to ANN [61,62]. We are not indicating that ANN has no potential capability for soil HMC because not many different architectures with dense grid search for hyperparameters are tested exhaustively. Some research demonstrated that ANN performed well on soil property estimation [18,44]. Particularly for the HMC retrieval from Landsat images, ANN outperforms PLSR, and the R${}^2$ exceeds 0.8 for the concentration estimation of As, Cd, Ni and Pb [50]. This inconsistency in ANN performance between this work and other works is attributed to the different heavy metal categories, study area, training samples, network structures and hyperparameters, which need further study.

5.4. Cu Concentration Mapping

To evaluate the Cu concentration map, the classification map of this study area is introduced, which is provided by China Land Consolidation and Rehabilitation. The red regions on the classification map (Figure 7c) represent the mining area, the positions of which are consistent with the distribution of red areas representing a high Cu concentration in the Cu distribution map (Figure 7a). In the areas of farmland and forest (yellow and green on the classification map), the Cu concentration is much lower (blue areas of Figure 7a). The Cu concentration map can be explained by the classification map by assuming that Cu contamination is caused by the mining activities of humans.

Comparing the obtained Cu concentration map using the proposed pipeline with the interpolation Cu concentration map, it is observed that the distributions of Cu in these two maps share a similar distribution pattern with different spatial coverages and resolutions. Figure 7a is mapped based on pixels, and Figure 7b is produced based on sample points, so Figure 7a contains more spatial details. Another advantage of the proposed method is that the trained model has better extrapolation abilities, while the prediction area of Figure 7b is limited due to the poor extrapolation capability of interpolation methods.

5.5. Uncertainty Analysis

Although the proposed data processing pipeline leads to fairly accurate soil Cu concentration estimation, there is still a discrepancy between the estimated values and the true values of the Cu concentration, which could be caused by uncertainties in both (i) feature extraction and regression models and (ii) data.

The model uncertainty is fundamentally caused by the lack of knowledge about the true relationship between Cu concentration and remote sensing observations. In this paper, an empirical approach consisting of feature extraction and regression models is used to approximate this true relationship. Nevertheless, this approximation will inevitably lead to model bias, especially if the feature extraction and regression models are not optimized. Previous research indicates that feature selection and model selection can produce better prediction accuracy [28,29,38], where cross-validation and bootstrapping approaches are used to better utilize limited training samples [30,37,63,64]. Although this paper tries to reduce this model bias by comparing and identifying the “best” feature extraction and regression approaches and the “best” data processing pipeline, there might still be remaining model bias due to limitations of our research. For example, this study is still limited due to the limited number of feature extraction approaches and regression methods, as well as the coarse grid search density. Furthermore, for model evaluation matrices, only some classical matrices are used in this paper. However, using more matrices, e.g., residual prediction deviation or the ratio of performance to interquartile distance, might help to evaluate models more comprehensively, which could be considered for future work.

The data uncertainty is essentially caused by the inexactness, incompleteness, and insufficiency of the training samples. For example, the spatial distribution and quantity of soil samples might not be able to well-characterize the heterogeneity of the study area. The samples are not uniformly distributed in all land use types due to access limitation. The samples from the forestland are few, which could easily lead to bias for the Cu estimation over the forestland. Furthermore, although the Cu concentration of soil samples are collected and measured by a professional institution, there might be inaccuracy and even errors in the provided results.

6. Conclusions

This paper has proposed a pipeline for retrieving and mapping soil Cu concentration from Landsat MTIs. First, repeated satellite observations were utilized to enhance the spectral information of images. Second, the preferred feature is extracted and selected by systematically evaluating the feature importance. Third, the preferred model is selected by the bias-reduced model evaluation.

MTI was demonstrated to be able to make a more reliable and accurate estimation of Cu concentration using empirical regression methods. By utilizing 11 images of Landsat-8 captured in different time phases, the mean adjusted R${}^2$ obtained by SVR using 20 repeated 6-fold cross-validations on 138 soil samples increases from 0.433 to 0.641. Based on the OB of the MTIs, feature extraction and selection further facilitate the improvement of Cu concentration estimation. The mean R${}^2$ of PLSR and ANN increase from 0.568 to 0.618 and from 0.476 to 0.528 separately, indicating the necessity and benefit of feature extraction and selection. Using cross-validation estimation, PLSR, SVR and ANN are compared unbiasedly based on their best performances using their favourite feature combinations. Although ANN is a popular regression method, in our work, SVR outperforms ANN by achieving a mean R${}^2$ of 0.641, which is 21.4% higher than ANN. RMSE, MAE and SE also support the highest generalization capability of SVR.

The preferred model with the highest R${}^2$ obtained by SVR is selected to estimate the Cu concentration in soil over the study area. Compared to the interpolation map, the Cu concentration distribution map generated by the recommended pipeline gives the pixel-based Cu estimation with more spatial detail and wider spatial coverage. It also shows a consistent spatial pattern with the ground-truth land cover classification map. The results show this model’s ability to perform large-scale soil HMC mapping from widely avaliable satellite imagery.