In recent years, heavy metal contamination in soil has emerged as a critical environmental concern worldwide. Copper (Cu), as an essential trace element for living organisms, becomes toxic when its concentration exceeds a certain threshold. Traditional laboratory‑based methods for quantifying Cu concentrations in soil are accurate but suffer from high costs, long turnaround times, and limited spatial coverage. To overcome these limitations, remote sensing technologies, especially hyperspectral imaging, have been widely explored for rapid and non‑destructive estimation of soil heavy metals. The development of low‑altitude UAV drone platforms equipped with hyperspectral sensors has further enabled flexible, high‑resolution data acquisition at the field scale, making it possible to capture subtle spectral variations induced by heavy metal stress. This study aims to develop a robust inversion framework for soil Cu content by integrating multiple spectral transformation techniques, novel three‑band spectral indices, and an advanced machine learning optimizer. The effectiveness of the proposed method is demonstrated through a case study conducted in a historical non‑ferrous metal processing area in northeastern China.
The primary objectives are: (1) to compare the performance of two‑dimensional (2D) spectral indices, three‑dimensional vegetation indices (3D‑SI), and newly constructed three‑dimensional baseline indices (3D‑BI) in terms of their correlation with soil Cu content; (2) to evaluate the impact of spectral preprocessing methods, namely Savitzky‑Golay (SG) smoothing, first derivative (FD) transformation, and continuum removal (CR), on model accuracy; (3) to optimize the hyperparameters of extreme gradient boosting (XGBoost) using the Newton‑Raphson‑based optimizer (NRBO) and examine the prediction improvements.
Materials and Methods
Study Area and Data Acquisition
The study area is an abandoned non‑ferrous metal processing plant located in a city of northeastern China. Due to decades of copper smelting and surface treatment activities, the topsoil has been heavily contaminated with copper. A total of 63 sampling points were established following a 40 m × 40 m grid, with some high‑pollution sub‑areas refined to a 20 m × 20 m grid. At each point, five sub‑samples were collected from the 0–0.5 m depth and mixed to form one composite sample. The samples were air‑dried, sieved through a 10‑mesh nylon sieve, and then subjected to chemical analysis using inductively coupled plasma mass spectrometry (ICP‑MS) to obtain the reference Cu concentrations.
Simultaneously, a UAV drone (DJI M600) was deployed to acquire hyperspectral imagery over the site. The UAV drone was equipped with a Resonon PikaL hyperspectral sensor covering 389–1023 nm with a spectral resolution of 4.2 nm (150 bands) and a spatial resolution of 15 cm. The flight was conducted on clear days between 11:00 and 14:00 local time at an altitude of 110 m, with 55% forward and side overlap. Preprocessing included radiometric calibration, atmospheric correction, and geometric mosaicing using dedicated software (Sbgcenter, Megacube, and ENVI). The final hyperspectral cube was used for all subsequent analyses.
Spectral Preprocessing
To reduce noise and enhance spectral features, three preprocessing strategies were applied to the raw reflectance obtained from the UAV drone data: Savitzky‑Golay smoothing (SG), first derivative transformation (FD), and continuum removal (CR). The SG smoothing used a window size of 11 and polynomial order of 2. The FD transformation was computed as the first derivative of the smoothed reflectance, which helps to eliminate baseline drift and resolve overlapping peaks. The CR technique normalizes the spectrum by dividing the reflectance by the convex hull envelope, emphasizing absorption features. These preprocessed spectra are referred to as OS (original SG‑smoothed), FDS (first derivative), and CRS (continuum removal).
Variable Iterative Space Shrinkage Approach (VISSA) Feature Selection
After preprocessing, the high‑dimensional spectral data still contain many redundant bands. The Variable Iterative Space Shrinkage Approach (VISSA) was employed to select the most informative bands for Cu estimation. VISSA iteratively shrinks the search space by evaluating the contribution of each variable through model population analysis. The optimal number of selected variables was determined when the root mean square error of validation (RMSEV) reached a minimum. For our dataset, 42 bands were retained as the optimal feature subset. These bands were then used to construct spectral indices.
Construction of Spectral Indices
Three categories of spectral indices were constructed using the VISSA‑selected bands:
- Two‑dimensional spectral indices (2D): Normalized Difference Index (NDVI), Difference Index (DSI), Ratio Index (RSI), and Soil‑Adjusted Vegetation Index (SAVI). Their formulations are:
$$ \alpha = \frac{\rho_i – \rho_j}{\rho_i + \rho_j} \tag{1} $$
$$ \beta = \rho_i – \rho_j \tag{2} $$
$$ \gamma = \frac{\rho_i}{\rho_j} \tag{3} $$
$$ \varepsilon = \frac{(1+L)(\rho_i – \rho_j)}{\rho_i + \rho_j + L} \tag{4} $$
where \(L\) is the soil adjustment factor (set to 0.5 in this study).
- Three‑dimensional vegetation indices (3D‑SI): Enhanced Vegetation Index (EVI), Normalized Spectral Index (NSI), Structure Insensitive Pigment Index (SIPI), and Triangular Spectral Index (TSI):
$$ \zeta = 2.5 \times \frac{\rho_k – \rho_j}{\rho_k + 6\rho_j – 7.5\rho_i + 1} \tag{5} $$
$$ \eta = \frac{\rho_k – \rho_j}{\rho_i} \tag{6} $$
$$ \theta = \frac{\rho_k – \rho_i}{\rho_k – \rho_j} \tag{7} $$
$$ \mu = 0.5 \times [120(\rho_k – \rho_i) – 200(\rho_j – \rho_i)] \tag{8} $$
- Newly constructed three‑dimensional baseline indices (3D‑BI): ABI1, ABI2, ABI3, and ABI4, designed to capture subtle interactions between bands that are sensitive to heavy metal stress:
$$ \xi = \frac{\rho_j + \rho_k}{\rho_i} \tag{9} $$
$$ \tau = (\rho_i – \rho_j) – (\rho_j – \rho_k) \tag{10} $$
$$ \psi = \frac{\rho_j – \rho_k}{(\rho_i – \rho_j) – (\rho_j – \rho_k)} \tag{11} $$
$$ \phi = \frac{\rho_j – \rho_k}{\rho_i + \rho_j + \rho_k} \tag{12} $$
For each category, all possible band combinations were enumerated using the VISSA‑selected bands (42 bands) or the full spectrum depending on the index structure. The top performing indices were identified through Spearman correlation analysis with measured Cu content.
Correlation Analysis
Spearman rank correlation coefficient (\(|r|\)) was used to evaluate the relationship between each spectral index and soil Cu concentration. The index with the highest absolute correlation coefficient under each preprocessing strategy was selected for subsequent modeling. The analysis was performed separately for 2D, 3D‑SI, and 3D‑BI categories.
Machine Learning Models
Two modeling approaches were compared: XGBoost and NRBO‑XGBoost. XGBoost is a powerful gradient boosting algorithm that builds an ensemble of decision trees. Its performance heavily depends on hyperparameter tuning. The NRBO (Newton‑Raphson‑Based Optimizer) is a metaheuristic algorithm inspired by the Newton‑Raphson search rule (NRSR) and includes a trap‑avoidance operator (TAO) to escape local optima. NRBO automatically searches for the optimal hyperparameters of XGBoost, namely the learning rate (eta), the number of estimators (n_estimators), and the maximum tree depth (max_depth). The optimization bounds were set as: eta ∈ [0.03, 0.1], n_estimators ∈ [100, 1000], max_depth ∈ [3, 10]. The population size was 10, and the stopping criterion was based on a fixed number of iterations (100). The best configuration identified was eta = 0.03, n_estimators = 135, and max_depth = 3.
The entire dataset (63 samples) was split into training (80%, 50 samples) and testing (20%, 13 samples) sets. Statistical checks confirmed that the two subsets had similar Cu content distributions (t‑test p = 0.810). Model performance was evaluated using the coefficient of determination (\(R^2\)) and root mean square error (\(\sigma_{\text{RMSE}}\)):
$$ R^2 = 1 – \frac{\sum_{i=1}^n (y_i – \hat{y}_i)^2}{\sum_{i=1}^n (y_i – \bar{y})^2} \tag{13} $$
$$ \sigma_{\text{RMSE}} = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i – \hat{y}_i)^2} \tag{14} $$
Results and Discussion
Feature Selection and Correlation Analysis
VISSA selected 42 optimal bands from each preprocessing dataset. Using these bands, the Spearman correlation coefficients for the best index in each category were computed. Tables 1–3 summarize the optimal band combinations and the corresponding correlation coefficients for OS, FD, and CR preprocessing, respectively. For the two‑dimensional indices, the highest correlation was obtained by FD‑NDVI (|r| = 0.484) at bands 642 nm and 802 nm. For the three‑dimensional vegetation indices, FD‑TSI achieved |r| = 0.544. Among the new three‑dimensional baseline indices, FD‑ABI3 exhibited the strongest correlation with |r| = 0.570 using bands at 850 nm, 885 nm, and 520 nm. These results indicate that first derivative transformation significantly enhances the sensitivity of spectral indices to Cu content, and that three‑band indices generally outperform two‑band indices in capturing the complex spectral response caused by Cu contamination.
| Index category | Best index | Band combination (nm) | |r| |
|---|---|---|---|
| 2D | NDVI | 466, 470 | 0.322 |
| 2D | RSI | 894, 988 | 0.328 |
| 2D | DSI | 466, 470 | 0.323 |
| 2D | SAVI | 466, 470 | −0.322 |
| 3D‑SI | EVI | 524, 458, 561 | 0.524 |
| 3D‑SI | NSI | 385, 466, 475 | 0.361 |
| 3D‑SI | SIPI | 549, 385, 558 | −0.491 |
| 3D‑SI | TSI | 495, 483, 466 | −0.451 |
| 3D‑BI | ABI1 | 458, 574, 512 | −0.506 |
| 3D‑BI | ABI2 | 466, 483, 491 | −0.432 |
| 3D‑BI | ABI3 | 483, 495, 458 | −0.496 |
| 3D‑BI | ABI4 | 466, 475, 898 | −0.341 |
| Index category | Best index | Band combination (nm) | |r| |
|---|---|---|---|
| 2D | NDVI | 642, 802 | 0.484 |
| 2D | RSI | 454, 466 | 0.438 |
| 2D | DSI | 537, 470 | 0.425 |
| 2D | SAVI | 553, 470 | 0.484 |
| 3D‑SI | EVI | 802, 397, 553 | −0.544 |
| 3D‑SI | NSI | 446, 553, 483 | 0.455 |
| 3D‑SI | SIPI | 993, 454, 881 | 0.535 |
| 3D‑SI | TSI | 454, 470, 553 | −0.544 |
| 3D‑BI | ABI1 | 553, 650, 446 | −0.307 |
| 3D‑BI | ABI2 | 553, 470, 591 | −0.528 |
| 3D‑BI | ABI3 | 850, 885, 520 | −0.570 |
| 3D‑BI | ABI4 | 780, 890, 997 | 0.560 |
| Index category | Best index | Band combination (nm) | |r| |
|---|---|---|---|
| 2D | NDVI | 625, 570 | 0.342 |
| 2D | RSI | 723, 719 | 0.337 |
| 2D | DSI | 625, 570 | 0.333 |
| 2D | SAVI | 723, 710 | 0.354 |
| 3D‑SI | EVI | 854, 642, 650 | 0.494 |
| 3D‑SI | NSI | 625, 1020, 389 | 0.395 |
| 3D‑SI | SIPI | 655, 574, 961 | 0.446 |
| 3D‑SI | TSI | 702, 710, 723 | −0.370 |
| 3D‑BI | ABI1 | 625, 659, 579 | −0.356 |
| 3D‑BI | ABI2 | 541, 579, 625 | −0.367 |
| 3D‑BI | ABI3 | 566, 562, 1002 | 0.453 |
| 3D‑BI | ABI4 | 710, 723, 591 | 0.376 |
The results consistently show that the first derivative transformation (FD) yields the highest correlations across all index categories. The best overall Spearman correlation was |r| = 0.570 achieved by FD‑ABI3. This demonstrates that the combination of derivative preprocessing and a three‑band baseline index effectively extracts subtle spectral information related to Cu contamination. In contrast, continuum removal produced lower correlations, possibly due to its sensitivity to baseline curvature that may obscure Cu‑related features.
Model Performance Comparison
Using the best indices from each category and preprocessing, we built XGBoost and NRBO‑XGBoost models for Cu content inversion. Table 4 presents the training and testing results for all combinations. The dataset notation follows the pattern: “preprocessing + index category” (e.g., OS‑2D). The best performance was obtained using FD‑3D‑SI (first derivative + three‑dimensional vegetation indices) with NRBO‑XGBoost, yielding a training \(R^2\) of 0.933, training \(\sigma_{\text{RMSE}}\) of 75 mg·kg⁻¹, testing \(R^2\) of 0.725, and testing \(\sigma_{\text{RMSE}}\) of 163 mg·kg⁻¹. Compared to the unoptimized XGBoost (testing \(R^2\) = 0.570), NRBO‑XGBoost improved the testing \(R^2\) by 27%.
Other notable combinations include FD‑3D‑BI with NRBO‑XGBoost (testing \(R^2\) = 0.655) and OS‑3D‑SI with NRBO‑XGBoost (testing \(R^2\) = 0.633). The baseline XGBoost models generally produced lower testing accuracies, indicating that hyperparameter optimization using the Newton‑Raphson‑based scheme is essential for achieving robust predictions.
| Dataset | Model | Training \(R^2\) | Training \(\sigma_{\text{RMSE}}\) (mg·kg⁻¹) | Testing \(R^2\) | Testing \(\sigma_{\text{RMSE}}\) (mg·kg⁻¹) |
|---|---|---|---|---|---|
| OS‑2D | XGBoost | 0.600 | 196 | 0.305 | 236 |
| OS‑2D | NRBO‑XGBoost | 0.642 | 178 | 0.412 | 221 |
| FD‑2D | XGBoost | 0.785 | 135 | 0.435 | 211 |
| FD‑2D | NRBO‑XGBoost | 0.816 | 123 | 0.515 | 201 |
| CR‑2D | XGBoost | 0.806 | 240 | 0.443 | 236 |
| CR‑2D | NRBO‑XGBoost | 0.844 | 222 | 0.352 | 188 |
| OS‑3D‑SI | XGBoost | 0.769 | 138 | 0.618 | 187 |
| OS‑3D‑SI | NRBO‑XGBoost | 0.877 | 104 | 0.633 | 176 |
| FD‑3D‑SI | XGBoost | 0.833 | 116 | 0.570 | 213 |
| FD‑3D‑SI | NRBO‑XGBoost | 0.933 | 75 | 0.725 | 163 |
| CR‑3D‑SI | XGBoost | 0.577 | 241 | 0.340 | 281 |
| CR‑3D‑SI | NRBO‑XGBoost | 0.788 | 133 | 0.357 | 256 |
| OS‑3D‑BI | XGBoost | 0.562 | 154 | 0.433 | 253 |
| OS‑3D‑BI | NRBO‑XGBoost | 0.757 | 147 | 0.454 | 212 |
| FD‑3D‑BI | XGBoost | 0.805 | 121 | 0.404 | 266 |
| FD‑3D‑BI | NRBO‑XGBoost | 0.898 | 95 | 0.655 | 148 |
| CR‑3D‑BI | XGBoost | 0.683 | 179 | 0.423 | 223 |
| CR‑3D‑BI | NRBO‑XGBoost | 0.735 | 157 | 0.650 | 171 |
The superiority of the NRBO‑XGBoost model can be attributed to its ability to explore the hyperparameter space more effectively using the Newton‑Raphson search rule and trap‑avoidance operator. This avoids suboptimal local minima that plague manual or simple grid search. Furthermore, the use of first derivative preprocessing and three‑dimensional vegetation indices (FD‑3D‑SI) provided the most informative spectral features for Cu estimation, as the derivative enhances local absorption features related to metal–organic complexes and clay minerals, while the three‑band indices capture their synergistic interactions.
The predicted Cu concentrations from the best model (NRBO‑XGBoost, FD‑3D‑SI) were spatially mapped across the study area. The predicted values ranged from 497 mg·kg⁻¹ to 1301 mg·kg⁻¹. Although these values are below the Chinese standard risk screening level for construction land (18000 mg·kg⁻¹), the spatial variability is pronounced: the ratio of high to low predictions reaches a factor of 2.6. This indicates that the contamination is highly heterogeneous, which conventional sparse sampling would likely fail to capture. The model successfully identified localized hotspots, providing a valuable tool for targeted remediation efforts. The validation using 13 independent samples confirmed a testing \(R^2\) of 0.725 with a scatter around the 1:1 line, demonstrating the model’s reliability for field applications.
Conclusion
This study presented a comprehensive framework for retrieving soil Cu content from UAV drone hyperspectral imagery by coupling spectral transformation, three‑band spectral indices, and Newton‑Raphson‑based optimization of XGBoost. The main findings are:
- First derivative transformation (FD) significantly improved the correlation between spectral indices and Cu content compared to SG smoothing and continuum removal. The best FD‑ABI3 index achieved a Spearman correlation of |r| = 0.570.
- Three‑dimensional spectral indices (both vegetation and baseline types) consistently outperformed two‑dimensional indices in correlation magnitude, demonstrating the advantage of using multi‑band interactions for heavy metal estimation.
- The NRBO‑XGBoost model, which automatically searches for optimal hyperparameters using a Newton‑Raphson metaheuristic, outperformed the default XGBoost model by a large margin. The best combination (FD‑3D‑SI + NRBO‑XGBoost) yielded a testing \(R^2\) of 0.725 and a testing \(\sigma_{\text{RMSE}}\) of 163 mg·kg⁻¹, compared to \(R^2\) = 0.570 for standard XGBoost.
- The spatial mapping capability of the UAV drone platform, combined with the optimized model, allowed us to visualize the heterogeneous distribution of Cu across the site, highlighting potential hotspots that require priority attention.
Despite the promising results, this study has limitations. The sample size (63) is relatively small, and the model’s generalization ability should be further validated across diverse environments and seasons. The influence of factors such as soil moisture, organic matter, and vegetation cover at the time of UAV drone flight could not be fully decoupled. Future work should incorporate multitemporal UAV drone acquisitions and include ancillary environmental covariates to enhance model robustness. Additionally, transferability to other heavy metals (e.g., Pb, Zn, As) and different soil types should be explored. Overall, the proposed methodology demonstrates the potential of UAV drone hyperspectral remote sensing as a rapid and accurate approach for soil heavy metal mapping, supporting environmental monitoring and risk assessment.