Evaluating Metal Effects on the Reflectance Spectra of Plant Leaves during Different Seasons in Post-Mining Areas, China

This study examined the relationship between the leaf reflectance of different seasons and the concentration of heavy metal elements in leaves, such as Co, Cu, Mo, and Ni in a post-mining area. The reflectance spectra and leaf samples of three typical plants were measured and collected in a whole growth cycle (June, July, August, and September). The Red Edge Position (REP), Readjustment Normalized Difference Vegetation Index (RE-NDVI), and Photochemical Reflectance Index (PRI) were extracted and used to explore its relation with the heavy metals concentrations in leaves between different seasons. The results show that all three Vegetation Indices (VIs) were insensitive indicators for monitoring the metal effects of vegetation in different seasons, which showed similar trends. Based on this, the Continuum Removal Indices (CRIs) were proposed from the continuum removed approach and extended for detecting the effects of heavy metal pollution over a full growth cycle. The relationship between the metal concentrations and CRIs of different plants was respectively analyzed by Stepwise Multiple Linear Regression (SMLR) and Partial Least Squares Regression (PLSR). It is found that a significant correlation exists between the band depth and the concentration of Cu and Ni based on the White birch data sets using the PLSR, resulting in a small deviation from the established relationships. Compared with VIs, the approach of coupling CRIs and multiple regressions was effective for improving the estimation accuracy. The presented study provides a detection model of leaf heavy metals that can be adapted to different growing cycles, even an arbitrary growing cycle.


Introduction
One of the major environmental problems resulting from mining areas is the metal pollutants [1][2][3]. Currently, heavy metal pollution monitoring mainly depends on the traditional geochemical methods, which are time consuming with a low efficiency and unsuitable for large-scale monitoring [4,5]. However, the vegetation was shown to play an important role in element migration, which could adsorb and transport ore-forming elements to the surface and result in a near-surface medium anomaly [6]. When plants uptake excessive concentrations of heavy metals in soil, their growth will lead to a certain extent of harmful changes, such as a decrease in chlorophyll content and/or changes in the leaf internal structure [7,8]. These directly or indirectly influence the spectral feature of the plant leaves [9,10]. Therefore, exploring the spectral response of the plant leaves by metal-induced stress is of the utmost importance for monitoring the phytoremediation of metal mine-polluted land or finding the concealed deposits using remote sensing techniques, with its advantages of fast dynamic, subtle recognition, and nondestructive detection.
Early spectral analysis of the interaction between vegetation and metal dates back to the period of 1970-1980, such as the works by Collins et al. [11], Milton et al. [12], and Horler et al. [13], who defined the Red Edge Position (REP) and demonstrated shifts in the REP of vegetation reflectance-based heavy metals stress. Based on this, a substantial amount of studies explored the effects of the metals on vegetation spectra. The most fundamental characteristic is the 'blue shift' of REP that occurs when vegetation is subjected to a heavy metal concentration, which refers to the spectra that tend to shift toward the left and shorter wavelengths. Stress also tends to produce an increase in reflectance at around 680 nm absorption. Some previous studies have been conducted under controlled laboratory conditions rather than field environments. In the laboratory, the specific plants were grown at the controlling levels of the certain metal over a full growth cycle [14][15][16][17]. Many of these studies tended to identify comparatively obvious spectral features by metals stress, which often manifests as a higher simulation accuracy of heavy metal content than the field studies. However, the vegetation spectral properties are a function that depends on the many geographical variations and applying the results of the laboratory to the natural environment still needs to be explored [17][18][19]. Field studies that examine the plant response to metal-contaminated soil have mostly focused on old waste deposit sites [1][2][3], river floodplains [19,20], polluted farms [21,22], etc. In these studies, the dominant species of plant in specific growth stages were selected as the study object, and then the VIs were used to explore relationships between the metal content and spectral response. The spatial distribution of metal contaminants can be indirectly inferred from these relationships at the leaf scale or canopy level [22,23]. Although many scholars have carried out relevant research work and developed some meaningful VIs, there have been no field studies of the multiple-metals effect on the spectra of varying vegetation species in the different seasons. On the other hand, the metal trans-location properties from roots to leaves for different species of plant still need to be discussed, which is a necessary first-step for exploring the spectral response of the plant by metal-induced stress [24,25].
Several VIs have been developed and used to assess plant stress. These VIs using reflectance measurements in the red and near-infrared regions are sensitive even to small changes in vegetation status produced by heavy metals stress. The most widely used VI is REP, and the shift of REP towards shorter or longer wavelengths, respectively, was used to detect the stress or growth of green vegetation [13,26]. Currently, this fundamental spectroscopic-plant principle in these previous studies is widely used in vegetation detection [9,17]. Many researches focused on exploring the change of several existing VIs around the red edge corresponding to stress, such as the Normalized Difference Vegetation Index (NDVI) [2,22,27], Photochemical Reflectance Index (PRI) [14,28], and Ratio Vegetation Index (RVI) [16,22], and wavelet-fractal analysis [29] has also been demonstrated to be a successful indicator of a variety of metal stresses. However, these VIs were applied to different species, with no comparison with other techniques and without being applied to the plant in different seasons, and the optimum stages of vegetation samples for detecting metals were not guaranteed to have been collected in the field. Meanwhile, the majority of the well-succeeded works on the VIs developed and applied have been achieved by exploring detached leaves rather than canopies, and a variety of classic VIs (e.g., REP, PRI) were developed at leaf level. Thus, it is particularly important to explore a hyperspectral detection mode of a heavy metal at the leaf level.
The continuum removal method at a specific waveband is another spectral analysis approach that has been assessed for biochemical parameter estimation [30,31]. The individual absorption features of continuum removal can be compared by normalizing them to the common baseline [32], and this method and its derived parameters have been proved useful for mineral mapping [33] and have also been extended to estimating the foliar biochemical parameters (nitrogen, lignin, and cellulose, etc.) [34][35][36]. Limited work has used the continuum removal to assess the plant response to metal-include due to geochemical stress or for old waste deposit sites. Thus, the comparison of VIs and CRIs utilizing hyperspectral vegetation data sets at different growth stages for the detection of multi-metals pollutants is necessary and forms the focus of this study.
Multivariate regression has become a promising method to process the large amount of spectral information, for example, SMLR has been successfully used to simulate the vegetation parameters [34,37]. However, it may be impacted by multicollinearity in the processing of hyperspectral data. This typically happens when the number of wavelengths studied is larger than the number of observations [38]. In contrast, PLSR was first developed to deal with the multicollinearity problem in chemical application areas, and this statistical technique has been successfully extended to simulate vegetation biochemical properties [39,40] and biophysical properties [41,42]. An optimal selection of spectral bands and an appropriate statistical method for heavy metal estimation among the vast array of those available are required.
For the above reasons, the purposes of the research are to: (i) explore the sensitivity of the most common VIs used to monitor the metal effects in different growth stages; and (ii) extend the CRIs derived from continuum removal to analyse the plant response to multiple-metals over a full growth cycle. To achieve these objectives, this study measured the spectra of a native representative plant in the waste copper mines over a full growth cycle to extract three representative VIs and CRIs. Univariate regression with VIs and two multivariate regression methods with CRIs were applied to evaluate the relation between vegetation reflectance and metal concentration. The response of the hyperspectral indices of different plant species under different metal levels and seasons by this relation was discussed.

Study Area
The study area is located in the Duobaoshan copper mine of the Lesser Khingan, Heihe City, in the Midwest of Heilongjiang Province, China. According to the retrieval data from the website of Weather China, Heihe experiences a monsoon-influenced humid continental climate, and the average annual rainfall is 490~540 mm, close to two-thirds of that in the months of June to August. It is rich in polymetallic mineral resources of Co, Cu, Mo, and Ni. The two sampling points were distributed in the two main ore deposits located roughly 3 km apart, which are large Cu-Mo deposits of the porphyry type. In addition, a control group was collected from Woduhe forest farm situated to the north of the mine. In sum, three sampling points were selected for the experiment, including one control point (S 0 ) in the forest farm and two points (S 1 and S 2 ) in the copper mine (as shown in Figure 1).
In the two main ore deposits, due to artificial exploitation, two vertical abandoned pits of accumulated waste rock have formed, which was easy to use for producing the ore-forming elements' migration and resulted in a near-surface medium anomaly. The herbaceous and woody specimens have always coexisted at sampling points. The major plant species include White birch, Mongolian oak, and Cotton grass (as shown in Figure 2). As the vegetation coverage is relatively high in the two sampling points, this provides an opportunity to assess the effects of the metal on leaf reflectance. Meanwhile, the repair of abandoned pits using the vegetation is being conducted by Duobaoshan Copper Ltd. (Heihe, Heilongjiang, China). Thus, this is an ideal area for the study of the spectral response of metals-stressed leaves from the perspectives of the mineral environment and resources. In the two main ore deposits, due to artificial exploitation, two vertical abandoned pits of accumulated waste rock have formed, which was easy to use for producing the ore-forming elements' migration and resulted in a near-surface medium anomaly. The herbaceous and woody specimens have always coexisted at sampling points. The major plant species include White birch, Mongolian oak, and Cotton grass (as shown in Figure 2). As the vegetation coverage is relatively high in the two sampling points, this provides an opportunity to assess the effects of the metal on leaf reflectance. Meanwhile, the repair of abandoned pits using the vegetation is being conducted by Duobaoshan Copper Ltd. (Heihe, Heilongjiang, China). Thus, this is an ideal area for the study of the spectral response of metals-stressed leaves from the perspectives of the mineral environment and resources.

Field Ground Data Collection
ASD FieldSpec Pro (Analytical Spectral Devices, Boulder, CO, USA) is a common portable spectroradiometer, which has been widely used in the detection and analysis of stress caused by environmental pollution [16,17,19,43]. In this study, leaf reflectance spectra were measured in situ using an ASD FieldSpec 3 High-Resolution Pro spectrometer, with a 25° field of view (FOV) through a permanent fiber optic probe. This instrument detects electromagnetic radiation in the wavelength range of 350-2500 nm; a sampling interval of 1.4 nm between 350-1000nm and 2 nm between 1000-2500 nm; and a spectral resolution of 3 nm (@ 350-1000 nm), 8.5 nm (@ 1000-1900 nm), and 6.5 nm (@ 1900-2500 nm).
For measuring leaf spectra, in order to reduce the effects of the sun's position and ensure the minimum amount of samples required for element testing, two to three leaves were randomly selected from each plant and cut at the leaf stalk using steel scissors. Each leaf was placed on a black paperboard in the open ground, to minimize the background spectral noise or radiation transmitted through the leaf. The foreoptics were aligned vertically and the height of the foreoptics was kept within 5 cm above the leaf surface at a nadir position. The leaf area observer by the ASD sensor was less than 1 cm in diameter, which ensured that the FOV filled the leaf without being influenced by the surroundings. Prior to each leaf reflectance measurement, a reflectance reference measurement with a Spectralon panel (Labsphere Inc., North Sutton NH, USA) was measured under the same conditions immediately. The final conversion to spectral reflectance was done by dividing the measured spectrum of the leaf samples by that of the Spectralon panel. The ASD was configured to average 15 consecutively acquired spectra for each record. For each plant, the mean of all spectra from two to three leaves in different positions was calculated for subsequent analysis, in order to average out the differences in reflectance as the result of different angles. All spectral measurements were taken under cloudless or near-cloudless conditions between 11:30 a.m. and 2:00 p.m. In order

Field Ground Data Collection
ASD FieldSpec Pro (Analytical Spectral Devices, Boulder, CO, USA) is a common portable spectroradiometer, which has been widely used in the detection and analysis of stress caused by environmental pollution [16,17,19,43]. In this study, leaf reflectance spectra were measured in situ using an ASD FieldSpec 3 High-Resolution Pro spectrometer, with a 25 • field of view (FOV) through a permanent fiber optic probe. This instrument detects electromagnetic radiation in the wavelength range of 350-2500 nm; a sampling interval of 1.4 nm between 350-1000 nm and 2 nm between 1000-2500 nm; and a spectral resolution of 3 nm (@ 350-1000 nm), 8.5 nm (@ 1000-1900 nm), and 6.5 nm (@ 1900-2500 nm).
For measuring leaf spectra, in order to reduce the effects of the sun's position and ensure the minimum amount of samples required for element testing, two to three leaves were randomly selected from each plant and cut at the leaf stalk using steel scissors. Each leaf was placed on a black paperboard in the open ground, to minimize the background spectral noise or radiation transmitted through the leaf. The foreoptics were aligned vertically and the height of the foreoptics was kept within 5 cm above the leaf surface at a nadir position. The leaf area observer by the ASD sensor was less than 1 cm in diameter, which ensured that the FOV filled the leaf without being influenced by the surroundings. Prior to each leaf reflectance measurement, a reflectance reference measurement with a Spectralon panel (Labsphere Inc., North Sutton NH, USA) was measured under the same conditions immediately. The final conversion to spectral reflectance was done by dividing the measured spectrum of the leaf samples by that of the Spectralon panel. The ASD was configured to average 15 consecutively acquired spectra for each record. For each plant, the mean of all spectra from two to three leaves in different positions was calculated for subsequent analysis, in order to average out the differences in reflectance as the result of different angles. All spectral measurements were taken under cloudless or near-cloudless conditions between 11:30 a.m. and 2:00 p.m. In order to reduce the effects of the illumination changes, the Spectralon panel was measured for 20 min of use or the light changed significantly in twenty minutes.
Leaf samples were collected almost synchronously with the leaf spectral reflectance measurements in the field. In addition to the collection of the measured leaves, the plant roots were immediately collected for evaluating the translocation characteristics of heavy metal from the root to leaf by the calculation of the Translocation Factor (TF). The TF for metals within a plant was expressed as the ratio of [metal] Leaf /[metal] Root , which has been widely used to evaluate the capability of plants to accumulate the metal from roots to above ground parts [24,25]. All samples were put into sample bags with labels, which were stored in an insulated cabinet filled with ice. The ice was changed every day for keeping the samples fresh. The samples were sent to a commercial laboratory and analyzed for total metal concentrations of Co, Cu, Mo, and Ni by flame atomic absorption spectrometry (AAS).
Field work was initially conducted during four months of a typical plant growing season from June to September 2013, which corresponded to the tiller, booting, mature, and fall growth stages of vegetation. We considered that only sampling points from 2013 were used for building and validating the model, which might cause instability of the model, and we were prepared to collect additional

Vegetation Spectral Analysis
(1) REP: The red edge is defined as the slope of the steep rise in the reflectance spectra of the green plant between 690 and 740 nm [9]. A variety of algorithms for quantitative analysis of the REP, such as: (i) the first difference transformation [44], which needs narrow and wider bandwidth spectra and can enhance the absorption features that might be impacted by background absorptions [44,45]; (ii) the inverted Gaussian method, which smoothens the spectral signature in the red-edge region using the inverted Gaussian function with complicated calculation procedures [46]; (iii) and the linear methods, which only use three or four fixed wavelengths to calculate the REP [47,48], and which has the advantage of rapidity; however, it overestimates the REP when compared to the other two methods [48]. Meanwhile, the latter two methods did not make full use of the very fine spectral resolution of the spectroradiometer that was used [49]. Thus, we used the maximum of the first derivative to calculate REP that can be expressed as Equation (1), and REP is the wavelength of the maximum of FDT.
where FDT is the first-difference transformation at a wavelength i, at the midpoint between wavebands j and j + 1; R λ(j) is the reflectance at the j waveband; R λ(j+1) is the reflectance at the j + 1 waveband; and ∆ λ is the difference in wavelengths between j and j+1.
(2) VIs: VIs are a linear or nonlinear combination of values of two or more bands. Two types of two-band VIs were calculated, including RE-NDVI and PRI, which were confirmed for detection of the metal stress through vegetation reflectance in the literature [17,22,50].
(I) RE-NDVI: The RE-NDVI was based on the chlorophyll index developed by Gitelson and Merzlyak [51] and showed a high sensitivity to pigment changes, exhibiting a better stability than NDVI when applied across a wide range of species and plant functional types [52]. The RE-NDVI was derived from hyperspectral data and calculated by the following formula: (II) PRI: The PRI was originally developed by Gamon, Peñuelas, and Field to assess the rapid changes in the pigment concentrations due to a variety of stress conditions [28], and thus serves as an estimate of metal concentration. PRI was defined as follows: In Equations (2) and (3), R xxx refers to leaf reflectance at wavelength xxx in nanometers.
(3) Continuum Removal Indices (CRIs): The continuum is a convex hull (continuum line) fitted over the top of a spectrum to connect local spectral maxima [32]. The parameters of the continuum removal methodology are derived by calculating the continuum removed reflectance (CRR). The CRR is defined as the original reflectance (OR) value divided by the value of the continuum line (CL) for each waveband i of the absorption feature [47].
The start and end spectral values of a specific wavelength interval are on the hull, which are given the number 1 by the continuum removed. As a result, the relative value of the output curves will be between 0 and 1, where the absorption troughs are enhanced.
Previous studies show that the damage of plants caused by heavy metal decreases photosynthesis in vivo replacement of the central Mg 2+ ion in the chlorophyll molecule by a heavy metal ion [7]. Although the replacement amount depends on the metal variety, the plant photosynthesis is more or less inhibited by the replacement, which could affect the major absorption pit ranging from 550 nm to 750 nm [9]. In the short-wavelength infrared bands (SWIR, 800 nm~2500 nm), the vegetation spectrum is dominated by the internal leaf cellular structure and water absorption, which is in contrast to the near-infrared bands where chemical absorption is largely masked by water [26,30,53]. However, determining the metal stress from reflectance on a fresh leaf is extremely complex due, among other things, to the strong effect of water [54]. The visible bands (VIS) have a stronger capacity in the detection and analysis of heavy metals stress than the SWIR bands, which was confirmed by a previous study [14,22]. Therefore, we were interested in isolating the specific absorption feature in this study, and linear continua were fitted between the start (550 nm) and endpoint (750 nm) of the absorption feature and then continuum removal was applied. Three variables were derived by continuum removal as follows: (I) The band depth (BD) at each wavelength (i) within the absorption feature was calculated by subtracting the continuum removed reflectance (CRR) from 1 [55]: (II) A normalization method using band indices was developed to reduce the effects of the soil background and water [55]. The normalized band depth ratio (BDR) was calculated by dividing the BD of each wavelength i by the band depth center (BD c ), which is the maximum BD [32].
(III) The normalized band depth index (NBDI) was inspired by the form of the NDVI, which was calculated using the following formula [32,55]: where i is the particular wavelength between 550 nm and 750 nm.

Statistical Analysis
Analysis of variance (ANOVA) was used to ascertain the level of the statistically significant differences between the sampling points, and the normality and homogeneity of variance were evaluated prior to ANOVA, which was performed in SPSS 13 (SPSS, Inc., Chicago, IL, USA). The significance level (p) for ANOVA was classified as weakly significant (0.1 ≤ p < 0.05), significant (0.05 ≤ p < 0.01), or highly significant (p ≤ 0.01) [56]. Pearson's correlation coefficient (r) was used to identify the strength of a linear association between two variables in this study, and the strength of the correlation was set as weak (|r| ≤ 0.39), moderate (0.40 ≤ |r| ≤ 0.59), or strong (|r| ≥ 0.60) [57][58][59][60]. Two regression analysis methods (SMLR and PLSR) were performed to evaluate the relationship between metal concentrations and the leaf reflectance between 550 nm~750 nm, which were tested for the data sets of all vegetation samples and three divided vegetation samples.
Since the number of samples is limited, independent testing samples were not used to measure the validity of the models, but we used a cross-validation procedure. The cross-validation method has the ability to detect outliers and provide nearly impartial estimations of the simulation error [61,62]. It was conducted by assigning each sample population into groups, and each sample was evaluated by the rest of the samples. The statistical results were assessed according to coefficient of determination (R 2 ) and Root Mean Square Error of Cross-Validation (RMSE CV ). Since there is no corresponding module in SPSS, two regression analyses (i.e., SMLR and PLSR) were performed using the Statistics Toolbox within Matlab TM.
Stepwise regression fits a dependent dataset (i.e., metal concentration) using a linear combination of independent variables (i.e., spectral indices). SMLR computes the F-statistic for each wavelength and metal concentrations significance levels (p-values) for wavelengths to enter and to stay in the equation, which ends when none of the wavelengths outside the equation have a significance (p-values) at or below the enter level and all wavelengths in the equation are significant above the stay level. p-values to enter and move wavelengths were set at 0.01 and 0.02 [37,63], respectively, in order to have a highly significant level between the selection of the wavelengths and the metal concentrations. In comparison, PLSR is a method used to reduce the large amount of highly collinear spectral variables into a limited number of uncorrelated latent variables or factors while combining the most useful information, which is similar to the principal components approach [37]. Whereas principal component regression is based on the decomposition of spectral data alone, PLSR conducts the decomposition on both the spectral data and the response variable simultaneously. This process greatly avoids the latent over fitting problem while using as much information as possible [9,64].

Metal Elements in Leaf and Root
The variations of concentrations for Co, Cu, Mo, and Ni within all plant samples of the different sampling points are shown in Table 1. Among all the samples that were analyzed, the accumulation of Co, Cu, Mo, and Ni was significantly (p < 0.05) higher in the leaves of S 1 and S 2 , which was approximately three to five times that in S 0 . Meanwhile, the element content of S 2 was generally higher than that in S 1 . Comparing different seasons, the contents of four elements in August and September were higher than in the other two months. The average contents of four elements in S 1 and S 2 were relatively lower after repair by Duobaoshan Copper Ltd. (Heihe, Heilongjiang, China) (i.e., after 2016 data are in), which were still above those in S 0 . There was a ten-fold difference in the metal concentrations of leaves between Cu/Ni and Co/Mo. Co, Cu, and Mo concentrations in White birch and Cotton grass were higher than in Mongolian oak. Cotton grass had the lowest Ni concentration when comparing all the kinds of plant. The correlation coefficients (r) between the metal concentration in the leaf and root of different plant species at S 1 and S 2 sampling points were calculated and the results are shown in Table 2. White birch, Mongolian oak, and Cotton grass showed significant positive correlations and the significance level was p < 0.01. Among them, the average value of r for White birch and all four elements was the highest (r = 0.75), followed by Mongolian oak (r = 0.71), and Cotton grass (r = 0.65). This implies that the metal concentrations in the leaf increased linearly as the root concentration of the respective plant increased. Meanwhile, the TF of different plants at S 1 and S 2 sampling points is shown in Table 2. At S 1 and S 2 sampling points, for all elements, the TF of White birch was greater than 1.5. For Co, Cu, and Ni, the TF of Mongolian oak was greater than 1.0. For Cu and Mo, the TF of Cotton grass was greater than 1.0.

Analysis of REP, PRI and RE-NDVI
The averaged leaf reflectance spectra of different plant species at S 0 , S 1 , and S 2 sampling points are given in Figure 3a. In order to explain the consistency of spectral measurements, the vertical error bars were added in the form of standard deviations to Figure 3a (see Figure 3b). However, as space is limited, different plant species of all sampling points are not provided, so take S 1 sampling point as an example.
The differences in the reflectance of three plants between three sampling points were mainly located after 550 nm or 750 nm. With the influence of multiple factors such as leaf internal factors (water, leaf structure, etc.) and the external environment, the formation mechanism of vegetation reflectance after 750 nm was complicated, which resulted in the differences in reflectance and the standard deviation after 750 nm was relatively larger (as shown in Figure 3b). In comparison, the differences in reflectance located in the range of 550 nm~750 nm were relatively smaller, and there are some rules in terms of the change of green peak and red edge. Compared with the S 0 point, the green peak becomes shallow and the slope of the red edge becomes slow in the reflectance of three plants at the S 2 point that has a relatively higher heavy metals concentration. However, cotton grass displayed the opposite change, which is probably because the moderate consumption of heavy metal can promote the growth of plants.
Three VIs (REP, PRI, and RE-NDVI) results of different sampling points and plant species obtained during the different seasons are given in Figure 4. The plants from S 1 and S 2 showed no significant difference in VIs compared to plants of S 0 , and there were no more unified changing rules between the three points. In June, REP, PRI, and RE-NDVI were lower in S 0 than in S 2 ; however, the REP and PRI exhibited the opposite feature in July and August, which were higher in S 0 compared to S 1 and S 2 . In September, for the three VIs, there are no uniform rules between S 0 , S 1 , and S 2 , whereas the RE-NDVI of the S 1 and S 2 appears greatly different to S 0 . Differences in VIs were relatively small and were not in order between the three dominant plants. However, three VIs of the Mongolian oak in September were relatively lower, which was a relatively special result, and this was due to the senescence rate of the Mongolian oak in September being relatively quicker than in the other two plants. Three VIs (REP, PRI, and RE-NDVI) results of different sampling points and plant species obtained during the different seasons are given in Figure 4. The plants from S1 and S2 showed no significant difference in VIs compared to plants of S0, and there were no more unified changing rules between the three points. In June, REP, PRI, and RE-NDVI were lower in S0 than in S2; however, the REP and PRI exhibited the opposite feature in July and August, which were higher in S0 compared to S1 and S2. In September, for the three VIs, there are no uniform rules between S0, S1, and S2, whereas the RE-NDVI of the S1 and S2 appears greatly different to S0. Differences in VIs were relatively small and were not in order between the three dominant plants. However, three VIs of the Mongolian oak in September were relatively lower, which was a relatively special result, and this was due to the senescence rate of the Mongolian oak in September being relatively quicker than in the other two plants.
For different sampling points and different plant species, both VIs appeared to have significant differences (p < 0.05) during different seasons. Compared with July and August, the REP in June and September tended to shift toward the left to shorter wavelengths. The general trend for the changes in REP was that June was less than July and August and greater than September. The range of the 'red shift' from June to July and August increased by 13 nm~15 nm, and the range of the 'blue shift' form July and August to September descended by 18 nm~20 nm. Compared with PRI and REP, it was shown that both exhibited similar change trends from June to September. The PRI was negative in June and September and significantly lower than July and August. In contrast with REP and PRI, the RE-NDVI displayed no significant difference between June, July, and August. Nonetheless, the RE-NDVI in September was significantly lower than the other three months. For different sampling points and different plant species, both VIs appeared to have significant differences (p < 0.05) during different seasons. Compared with July and August, the REP in June and September tended to shift toward the left to shorter wavelengths. The general trend for the changes in REP was that June was less than July and August and greater than September. The range of the 'red shift' from June to July and August increased by 13 nm~15 nm, and the range of the 'blue shift' form July and August to September descended by 18 nm~20 nm. Compared with PRI and REP, it was shown that both exhibited similar change trends from June to September. The PRI was negative in June and September and significantly lower than July and August. In contrast with REP and PRI, the RE-NDVI displayed no significant difference between June, July, and August. Nonetheless, the RE-NDVI in September was significantly lower than the other three months.
In the next study, the data in S 1 and S 2 points were grouped into four data sets (all vegetation samples, White birch, Mongolian oak, Cotton grass) based on the plant species, to explore the relationship between VIs and the concentration of metal elements. The correlation coefficient (r) between the concentration of Co, Cu, Mo, and Ni and three VIs for four different data sets (n = 90) were calculated. The results are shown in Figure 5, where the correlativity between the concentration of four heavy metal and three VIs was relatively weak. The |r| value range between 0.122 and 0.505, and the best correlation, appears for the Cu concentration and REP of the White birch data set (r = −0.505). Compared with the r value that was obtained from the three VIs, the REP exhibited a better performance, and there was no significant difference between the other two VI S . Of the four metal elements, the correlations of Cu and Ni were stronger than Co and Mo. The VIs of most data sets show a negative correlation with the contents of heavy metal, except for the correlation between the VIs of the partial data set and Mo. Comparing the four different data sets, the results showed no significant differences (p > 0.05) between the r value of each metal element.  In the next study, the data in S1 and S2 points were grouped into four data sets (all vegetation samples, White birch, Mongolian oak, Cotton grass) based on the plant species, to explore the relationship between VIs and the concentration of metal elements. The correlation coefficient (r) between the concentration of Co, Cu, Mo, and Ni and three VIs for four different data sets (n = 90) were calculated. The results are shown in Figure 5, where the correlativity between the exhibited a better performance, and there was no significant difference between the other two VIS. Of the four metal elements, the correlations of Cu and Ni were stronger than Co and Mo. The VIs of most data sets show a negative correlation with the contents of heavy metal, except for the correlation between the VIs of the partial data set and Mo. Comparing the four different data sets, the results showed no significant differences (p > 0.05) between the r value of each metal element.

Results of Continuum Removal Analysis
(1) SMLR: All four CRIs (OR, BD, BDR, and NBDI) used the SMLR on the four different data sets to develop relationships between the four metal elements' (Co, Cu, Mo, and Ni) concentrations and the vegetation reflectance. The stepwise regression was run with a narrow range of multiple F to enter and move on OR, BD, BDR, and NBDI data and metal concentrations. First, the stepwise regression was applied to select the optimized bands between 550 nm and 750 nm to be included in the model. Then, the optimal bands were used for calculating the cross-validated statistics.
The results from the analyses of the four variables by the continuum-removed absorption features are shown in Table 3. Considering the R 2 and RMSECV values, most of the correlations established are weak to moderate, R 2 varies between 0.21 and 0.58 and the average R 2 of all four data sets is 0.37, and the RMSECV produced the same result as the R 2 . However, comparing the simulative ability of the four CRIs, BD, BDR and NBDI were better than OR, BD was relatively good, and there were no significant differences between BDR and NBDI. The simulating ability of Cu and Ni was superior to that of Co and Mo, and the White birch data set to Cu (R 2 from 0.47 to 0.58) and the Cotton grass data set to Ni (R 2 from 0.42 to 0.56) were expected to have better performances. Compared with the other three data sets, the simulating ability of all vegetation samples to all elements was relatively steady.
Between one and six optimized wavebands were selected, and two-three optimized bands were in the majority. The numbers of optimized bands from the BD data and NBDI data, respectively, were the highest (3~6) and the lowest (1~3). As shown in Figure 6, the optimized bands for simulating Co element content were changes in a wide range and ranged from 640 nm to 690 nm. The bands for Cu and Mo were more centrally distributed, which were all located around 570 nm and 665 nm, whilst Ni was more concentrated close to 650 nm and 670 nm.

Results of Continuum Removal Analysis
(1) SMLR: All four CRIs (OR, BD, BDR, and NBDI) used the SMLR on the four different data sets to develop relationships between the four metal elements' (Co, Cu, Mo, and Ni) concentrations and the vegetation reflectance. The stepwise regression was run with a narrow range of multiple F to enter and move on OR, BD, BDR, and NBDI data and metal concentrations. First, the stepwise regression was applied to select the optimized bands between 550 nm and 750 nm to be included in the model. Then, the optimal bands were used for calculating the cross-validated statistics.
The results from the analyses of the four variables by the continuum-removed absorption features are shown in Table 3. Considering the R 2 and RMSE CV values, most of the correlations established are weak to moderate, R 2 varies between 0.21 and 0.58 and the average R 2 of all four data sets is 0.37, and the RMSE CV produced the same result as the R 2 . However, comparing the simulative ability of the four CRIs, BD, BDR and NBDI were better than OR, BD was relatively good, and there were no significant differences between BDR and NBDI. The simulating ability of Cu and Ni was superior to that of Co and Mo, and the White birch data set to Cu (R 2 from 0.47 to 0.58) and the Cotton grass data set to Ni (R 2 from 0.42 to 0.56) were expected to have better performances. Compared with the other three data sets, the simulating ability of all vegetation samples to all elements was relatively steady.
Between one and six optimized wavebands were selected, and two-three optimized bands were in the majority. The numbers of optimized bands from the BD data and NBDI data, respectively, were the highest (3~6) and the lowest (1~3). As shown in Figure 6, the optimized bands for simulating Co element content were changes in a wide range and ranged from 640 nm to 690 nm. The bands for Cu and Mo were more centrally distributed, which were all located around 570 nm and 665 nm, whilst Ni was more concentrated close to 650 nm and 670 nm.
(2) PLSR: The relationships between elements' concentration and the four indices (OR, BD, BDR, and NBDI) were modeled using PLSR. The optimum number of PLS factors preventing over-fitting was determined by setting the criterion and adding an additional factor that must reduce the RMSE CV by >2%. The number of factors in the final model ranged from two to six. The derived assessing statistical indicators are given in Table 4. The simulated ability of PLSR was significantly superior to SMLR, R 2 varies between 0.32 and 0.77, and the average R 2 of all four data sets is 0.49. In addition to BD, BDR can also obtain a preferable simulation accuracy, and the lowest R 2 and the highest RMSE CV were still produced by using OR. Similar to SMLR, the simulating ability of the Cu and Ni was superior to that of Co and Mo, and a good accuracy was obtained for the White birch data set to Cu (R 2 from 0.43 to 0.77). However, the Mongolian oak data set to Ni (R 2 from 0.45 to 0.74) was expected to have a moderate performance, which was different from SMLR.   (2) PLSR: The relationships between elements' concentration and the four indices (OR, BD, BDR, and NBDI) were modeled using PLSR. The optimum number of PLS factors preventing over-fitting was determined by setting the criterion and adding an additional factor that must reduce the RMSECV by >2%. The number of factors in the final model ranged from two to six. The derived assessing statistical indicators are given in Table 4. The simulated ability of PLSR was significantly superior to SMLR, R 2 varies between 0.32 and 0.77, and the average R 2 of all four data sets is 0.49. In addition to BD, BDR can also obtain a preferable simulation accuracy, and the lowest R 2 and the highest RMSECV were still produced by using OR. Similar to SMLR, the simulating

Discussion
The vegetation analysis of this study revealed a latent effect of the formerly disposed ore deposits at S1 and S2. Since S1 and S2 are probably three to five times the metal concentration of the plant leaves at S0, it is clear that two sites are affected by heavy metal pollution. The wasted time of the pit in S2 was relatively short, according to the characteristic of most rainfall in summer, thus we

Discussion
The vegetation analysis of this study revealed a latent effect of the formerly disposed ore deposits at S 1 and S 2 . Since S 1 and S 2 are probably three to five times the metal concentration of the plant leaves at S 0 , it is clear that two sites are affected by heavy metal pollution. The wasted time of the pit in S 2 was relatively short, according to the characteristic of most rainfall in summer, thus we consider that there is less loss of mineral elements by flow erosion, which is probably the main reason that the element content of S 2 was higher than S 1 . The heavy metals readily combine with other substances in the leaves to form the chelate, which were not easy to discharge and easily accumulated [65,66]. The metal complexes tended to increase with plant growth, leading to the element contents in August or September being relatively higher. The major sources of pollution identified were open pits, as the metal ions migrated from the open pits by dissolution and infiltration into the soil, which were then absorbed by the plant root. There is good significance between the metal concentration in the root and leaf (r = 0.63~0.74), which indicated that the stem transports metal ions from the root, which are subsequently accumulated in the leaf. However, another factor that should not be neglected is the dust containing metal ions. Considering the formation mechanism of the vegetation spectrum, the reflectance after 750 nm with great difference appears to follow no discernible laws between the three sampling points. However, 550 nm~750 nm with little difference can be regarded as an indicator for heavy metals, which was in conformity with previous studies [14,22].
According to previous results, the shift of REP is strongly correlated with photosynthetic activity and shows a greatly sensitive indicator of the health of the vegetation. In addition, two types of two-band VIs (RE-NDVI and RPI) have been confirmed for detection of the metal stress on the vegetation in the literature. Thus, we attempted to correlate the metal concentrations to three VIs of the vegetation in different seasons. As shown in Figures 4 and 5, it is clearly shown that three VIs of the four months were significantly different. The correlation coefficients between REP, RE-NDVI, PRI, and the foliar concentration of heavy metals were relatively poor. Despite many studies having shown that when a plant is stressed by the heavy metals, the spectra tend to occur as the blue shift of the red edge, they only focused on a specific growth stage. However, the photosynthetic activity through various stages of senescence should be considered when investigating the sensitivity of VIs to heavy metals. The minor change of VIs is a comprehensive result for the absorbing material (chlorophyll and other pigments), which is influenced by the increasing photosynthetic activity and decreasing excessive heavy metals. The photosynthetic activity was gradually strengthened from June to August, which is usually accompanied by an increase in the indicative bands of the pigment as the plant absorbs more energy in the photosynthetic process. Thus, the shifts or changes of the three VIs well reflected a full growth cycle of vegetation, which was insensitive to the metal stress. Despite this, the VIs of the leaf spectra can still reflect the heavy metal influence on vegetation to a certain degree. The right amount of heavy metal can accelerate plant growth in the early growth stage, which can be expressed by comparing the three VIs of S 1 /S 2 points to the S 0 point in June (see Figure 4). The uptake and transport ability of heavy metals by the plant root are obviously improved and increased with plant growth. Heavy metals enriched in plant leaves caused the REP and PRI to be higher in S 0 compared to the S 1 and S 2 in July and August. Comparing the three VIs, as the most classical indicative factor of plant healthy status, REP obtained a better performance than the other two VIs; however, this result may be influenced by the illumination geometry with changes of the solar altitudinal angle of the study area. RE-NDVI was less sensitive to heavy metals and did not significantly differ between S 0 , S 1 , and S 2 points, which can also be confirmed by the correlation coefficient between the heavy metals concentrations and RE-NDVI (see Figure 5, |r| < 0.35).
The relationship between metal and continuum-removed absorption features was explained by multivariate analyses (SMLR and PLSR). According to the results in Table 4, the average regression determination coefficient for the four elements obtained from three CRIs by PLSR and SMLR are respectively increased by 26.7% and 38.9% compared to OR, and BD yielded relatively higher R 2 and lower RMSE CV , which are, respectively, 44.7% and 53.4%. CRIs benefited from the continuum removal method and intensified the difference of the absorption strength, which shows that the continuum removal methods are better than original reflectance for monitoring elements. Among the two statistical methods, the PLSR achieved the highest R 2 values, which is consistent with the researchers demonstrating that the multivariate method improved the simulation performance, such as PLSR combined with continuum-removed absorption features for foliar biochemical parameter estimation. REP, RE-NDVI, and PRI were designed based on the optical mechanism of the plants while being combined by the fixed band. Figure 5 show that the VIs had a poor response on the plants of a full growth cycle affected by heavy metals. On the contrary, the optimal bands or factors were selected by the SMLR or PLSR, which might not have comparatively clear and recognized physical significance. However, most of the wavelengths selected for the four metal elements using SMLR are positioned close to the known Chlorophyll a and b absorption wavelengths (570 nm, Penuelas et al. [67]; 640 nm, Kumar et al. [68]; 660 nm, Curran, [55]), which are strongly influenced by the photosynthesis of vegetation.
As for the different metals, simulation models of Cu and Ni produced higher accuracies. Co and Mo showed a higher deviation from the established relationships, resulting in smaller R 2 values and larger RMSE CV . The explanation of this result is that the two mines are large Cu-Mo deposits of the porphyry type, which leads to the background value of Cu being relative higher. Despite this, the study area is rich in polymetallic mineral resources of Co, Cu, Mo, and Ni, all of which have been only found in the soil. There was a ten-fold difference in the metal concentrations of leaves between Cu/Ni and Co/Mo (see Table 1), which suggests that the plants possess the properties of heavy metal elements absorbing selectively. Meanwhile, the difference between the four elements in the TF indicated that the absorption capacity and tolerance of different plants to the four elements are diverse to an extent [25]. In addition, in the full growth cycle of vegetation, the various heavy metals existed as a mixture and their interactions strongly influence actually occurring affects, such as the synergistic and antagonistic action between different elements [69]. Cu and Ni are known as phytotoxins, which may have caused damage to the chlorophyll and leaf structure of the plants [70]. In this study, Co/Mo are probably not the principal reason for the spectral stress indicators. By comparing the capabilities for simulating the metal concentration with different data sets, satisfactory results were obtained using White birch. Meanwhile, White birch had a higher TF for all elements, which can be used as a dominant species for metal uptake in the hope that phytoremediation would be an alternative solution to mitigate the hot spots in the study site, and this is consistent with the main plant varieties planted in 2015 by Duobaoshan Copper Ltd. (Heihe, Heilongjiang, China). Therefore, White birch is one of the most susceptible indicator species for monitoring heavy metal pollution or for repairing ecological environments in this region.
The simulation accuracy of heavy metal is relatively lower in this study compared with in the greenhouse experiments. The determination coefficient (R 2 ) of the four elements ranged from 0.46~0.72, and was still at a middle level, which was derived from the BD data of all vegetation samples using PLSR. Here, a number of limiting factors need to be considered. Firstly, the stress levels are controlled in the greenhouse experiments and the researchers could pay attention to rule out factors so as to minimize the environmental impact. However, in the field, the element accumulation became stronger with the vegetation growth; meanwhile, the content of chlorophyll, water, macro-nutrients, and environment-induced noise had a corresponding change. So, many factors may have influenced the plants' spectra, so further experiments are required to identify the potential causes of spectral variability shown here. Secondly, heavy metal was a relatively independent factor affecting vegetation in the greenhouse experiments. The transport of heavy metal was limited by several factors in the field, such as rain wash, illumination changes, wind blowing, etc. The heavy metal concentrations of vegetation leaves were at a lower level, which results in variable features of vegetation spectrum not being significant.

Conclusions
In summary, our study compared the simulated ability of the three VIs and CRIs of heavy metal concentrations in leaves of different seasons. The estimation models were established by CRIs within the specific spectral range (550 nm~750 nm). From the result of the SMLR method, the majority of wavelengths selected for the four metal elements are positioned close to the known chlorophyll absorption positions. Comparing the SMLR and PLSR method, it appears that PLSR provides a better option for hyperspectral data mining. BD has the highest accuracy and is most sensitive to the metal-induced stress in the vegetation of different seasons, which is validated by different plant species in different seasons. Thus, when applying the results of this field study to airborne or satellite hyperspectral image sin this region, the model established by BD of White birch using PLSR is a relatively ideal choice. However, applying the result to evaluate pollution levels by the airborne hyperspectral image will need to consider additional factors, such as the spatial resolution, which is one of the primary problems. Meanwhile, the quality of the results of such investigations is influenced not only by the indices which best assess the vegetation status, but also by the correlation with the environmental on-site conditions. In conclusion, our study provides evidence of leaf spectrum alterations in different seasons caused by heavy metal pollution. We suppose that the presented approach can be applied to hyperspectral data of vegetation during different seasons and has potential for the monitoring of the heavy metal content of plants, and thus, the assessment of reclamation quality in post-mining regions.