礦山恢復治理區植被物候與健康狀況遙感監測

帥 爽，張 志，呂新彪，陳 思，馬梓程，謝翠容

（1. 中國地質大學（武漢）地質調查研究院，武漢 430074；2. 中國地質大學（武漢）地球物理與空間信息學院，武漢 430074；3. 湖北省國土測繪院，武漢 430010）

0 引 言

礦產資源的開發為國民經濟發展做出了巨大貢獻，但粗放式開發導致的耕地損毀、土壤侵蝕、粉塵污染、噪音污染、水污染以及對生物多樣性的影響逐步顯現。礦山恢復治理工程對降低礦產資源開發活動對生態環境的影響至關重要[1]。隨著國家生態文明建設戰略的實施，“使礦山生態保護、修復工作貫穿整個礦山生命周期”是礦山地質環境修復與治理工作的新理念。

對因采礦活動而退化的土地進行植被重建和恢復工作具有挑戰性，特別是礦山恢復治理區域土壤養分缺乏，不適當的護理或長期的擱置可能會永久性損害植被的健康[2]。恢復治理區植被健康狀況是評價恢復治理工作效果、研究礦山退化土地環境變化的重要指標之一[3]。

傳統的植被覆蓋變化調查和植被健康調繪主要采用野外調查和實驗室研究等常規方法，通常耗費大量人力物力，而且這些方法僅適合小區域的調查和評估[1]。前人利用多期次多光譜遙感數據結合各類型植被指數開展了大量礦山恢復治理區的監測和評估，取得了良好效果。Zhang等[4]驗證了礦山恢復治理區植被覆蓋度與植被指數的關系。Zenkov等[5]利用1984年至2016年多期次遙感數據研究了保加利亞TROJANOVO露天煤礦區的土地復墾動態變化過程。Xie等[6]提出了一種用于評估礦山損毀和恢復治理狀況的新型MRAIs （Termed Mining and Restoration Assessment Indicators）指數，相較傳統植被指數，MRAIs對礦山開發、修復引起的地類細微變化有更好顯示。Yang等[7]基于NDVI時間序列，根據礦產資源開發和恢復治理特征，定義了8種NDVI時間序列變化模板，再利用一系列閾值和動態匹配方法，對研究區像元進行歸類，獲取礦山開發擾動時間和恢復治理狀況。然而，這些研究僅針對礦區內土地利用類型變化和植被覆蓋度變化進行遙感監測，忽略了恢復治理區植被健康狀況的評估，另外在遙感圖像上區分恢復治理區植被與正常植被也是急需解決的問題。

礦山恢復治理區土壤質地、結構、含水量、孔隙度和重金屬含量等土壤質量指標的差異會引起植株生長稀薄、葉片短小、嚴重黃化[8-9]。由于反應植被健康程度的生理要素（如葉綠素含量、細胞結構、水分或氮含量）的變化可以通過反射光譜來監測[10]，部分研究者利用遙感技術評估恢復治理區植被健康程度和生長狀況。Lu等[11]通過對TM3與TM4波段構成的二維特征空間分析，獲得土壤線，將TM4波段坐標軸向土壤線旋轉，消除土壤影響，進而評價特定時間植被的生長狀態。胡秀娟等[12]利用WorldView-2衛星影像，基于歸一化山地植被指數、氮素反射指數和黃光波段反射率3個因子，提出了一種新型植被健康指數（VHI），以快速、大面積地監測與評價水土流失區的植被健康狀況。Karan等[1]利用2000年和2015年Landsat遙感數據和RVI、EVI和NDVI等植被指數監測露天煤礦恢復治理區植被變化，并利用遙感方法評估植被水分與植被健康程度的關系。然而，這些研究大多利用單一時相或多年相近月份的反射率光譜數據進行植被健康程度的評估，導致研究結果只能反映植被一個或幾個生長階段對土壤質量的響應，無法獲得植被整個生長階段對土壤質量的響應特征，可能導致結論的隨機性。

除此之外，礦山恢復治理區土壤質量退化還會導致作物物候特征的改變[13]。作物生長期長度和生長速率等物候參數是作物生長模擬的重要指標[13-15]，也能反映土壤的質量狀況。遙感植被指數的時間序列常被用于作物物候參數的提取和監測，該技術已經成功應用于作物類型識別、產量估算和生長狀態監測[16]，成為物候學與其他學科聯系的橋梁。但目前還鮮有通過研究礦山恢復治理區作物物候特征，評估恢復治理區植被健康狀況及土壤質量，進而評價礦山恢復治理工程效果的研究。另外，對于單個礦山區域監測，MODIS數據難以滿足作物識別的精度要求[17]，Landsat數據又難以滿足植被物候信息提取所需的時間分辨率。而哨兵2號遙感數據具備小區域植被物候特征研究所需的時間分辨率和空間分辨率，同時其在植被探測特征的可見光至近紅外波段具有光譜分辨率優勢[18]，與Landsat、ASTER和其他多光譜衛星數據相比，哨兵2號數據在植被紅邊區域包含更多波段，使得最初為高光譜數據開發的植被紅邊位置（REP）計算技術可用于哨兵2號數據。

植被指數與作物產量密切相關，同時能反應作物的養分和生理信息差異。植被指數作為單產估算模型因子與作物產量建立統計相關模型是一種常用的作物產量估算方法[19]。GVI、NDVI、EVI、NDVI累積值[20]作為模型輸入變量，建立玉米、小麥等作物產量統計相關模型，獲得了較高的精度。當作物受到脅迫作用時，相應的氮、色素、酶等發生變化，可應用各種植被指數監測這些生理指標變化[21]。MCARI、NDVI、NDMI等植被數據常被用于作物色素、氮含量和冠層葉片含水量的檢測[22]。通過分析恢復治理區和正常耕作區作物各類型植被指數值的差異，有利于獲取不同區域作物養分和生理信息差異。

本研究使用2020-4-12-10-31日采集的13期哨兵2號遙感數據開展玄武巖礦山恢復治理區植被健康狀況遙感診斷指標研究。主要內容包括：1）利用哨兵2號NDVI時間序列，基于NDVI時序重構、動態閾值、曲線曲率等方法，分別提取并對比恢復治理區和正常耕作區作物關鍵物候指標，分析恢復治理區作物與正常作物生長期內生長狀況的差異；2）提取研究區多類型植被指數、植被覆蓋度、植被紅邊位置等遙感指標參數，對比恢復治理區作物與正常作物各遙感指標參數差異，評價各遙感指標區分礦山恢復治理區作物與正常作物健康狀況的能力。為評價礦山恢復治理工程效果提供依據。

1 研究區及數據源

1.1 研究區概況

研究區位于黑龍江省七臺河市桃山區一處玄武巖采石場，地貌上為玄武巖丘陵，采石場于2012年開始開采，占地類型包括開采面、中轉場地、固體廢棄物和礦石堆等。玄武巖的開采活動造成了采石場及其周邊耕地的挖損和壓占。采石場北側兩處固體廢棄物區域（圖 1中A，B區域）于2019年進行了恢復治理，恢復治理類型為土地復墾，在對A、B區域固體廢棄物堆進行削坡和表面平整的基礎上，使用附近耕地土壤進行了堆墊，表層堆墊土壤厚度約為60 cm，堆墊土有機質、氮磷鉀等土壤化學特征與周圍耕地近似。但由于堆墊過程中，重型機械的碾壓，導致表層客土被壓實而土壤容重升高，進而影響植被扎根和水分與肥料的供給。A、B區域在土地復墾的基礎上，種植了玉米，玉米的播種、施肥與周邊農田同步。采石場周邊土地利用類型主要為玉米旱地和田間道路。該區域屬于溫帶季風氣候，年平均氣溫2.4～3.9 ℃，年均降雨量約490 mm。

1.2 遙感數據及處理

哨兵2號遙感數據由歐空局（https://scihub.copernicus.eu/）免費提供。由于黑龍江省玉米作物種植期主要為4月至10月，收集了研究區2020年4月至10月云覆蓋量小于5%的共13期（4月12日—10月31日）哨兵2號L1C級原始數據。利用Sen2cor軟件對L1C級哨兵2號數據進行輻射定標和大氣校正，獲取10 m空間分辨率的地表反射率數據，用于研究區植被物候參數、植被指數、紅邊位置、植被覆蓋度等信息的提取。不同時相的高空間分辨率遙感數據（高分2號、Google Earth數據）作為補充數據被用于對比分析，判斷恢復治理時間，勾繪恢復治理區域范圍。

1.3 野外調查

2020年8月24日開展野外調查工作，主要調查恢復治理區作物與正常作物的生長狀態，調查中發現礦山恢復治理區與周邊正常耕地玉米作物長勢存在明顯差異。礦山恢復治理區玉米植株生長稀薄，植株高度明顯低于周邊正常玉米作物，葉片短小且黃化現象明顯，生長狀態明顯劣于周邊正常耕作區域玉米作物（圖2）。

2 研究方法

研究方法包含兩部分，首先是建立研究區作物NDVI時序重構曲線并提取物候參數，在對哨兵2號影像進行大氣校正獲取反射率數據的基礎上，計算研究區NDVI時間序列，對比、選取適當濾波方法，重構NDVI時序曲線，提取關鍵物候參數，對比分析整個生長期內恢復治理區作物與正常作物生長狀況差異。其次是分別提取并對比恢復治理區作物與正常作物植被指數、植被覆蓋度、紅邊位置等遙感指標值差異，評估各遙感指標區分恢復治理區作物與正常作物的能力，和診斷作物健康程度的可行性和敏感性。具體研究方法和流程如圖3。

2.1 NDVI時序重構

在對研究區哨兵2號數據集進行大氣校正獲取反射率數據的基礎上，計算獲得NDVI時間序列。由于云層和陰影的影響，NDVI時序數據仍然存在較大噪聲，因此在應用之前必須對其進行重構[23]。以往研究中對NDVI時序重構方法的優劣沒有統一觀點[10,23]，本研究中，對比常用濾波方法（Savitzky-Golay（S-G）濾波、Double Logistic（D-L）擬合與Asymmetric Gaussian（A-G）濾波）的去除噪聲、保留植被生長特征的效果，對NDVI時序重構效果進行定性評價，同時使用統計分析對各濾波方法NDVI時序重構效果進行定量評價，選取最適合的去噪方法。S-G濾波器是一種基于滑動窗口的加權平均算法，它通過n階擬合多項式來計算某一點附近固定數量的點的平滑值，達到降噪效果。在S-G濾波器中，多項式次數和窗口大小的設置至關重要。D-L擬合和A-G濾波方法是采用特定函數對時間序列中最大值和最小值附近的數據進行擬合，并將局部擬合函數組合成整體擬合函數，增強擬合函數的靈活性，使其符合NDVI時間序列數據的復雜行為。

玉米營養生長期生長速率表現出“慢—快—慢”的規律，其中拔節期生長速率最快[24]。從圖4可以看出，這3種擬合方法均去除了大部分噪聲，在NDVI達到峰值前，3種擬合曲線均呈“S”型，符合玉米的生長趨勢，其中S-G擬合曲線與原始NDVI數據更接近，尤其在生長起始階段更能反映原始NDVI曲線的變化趨勢，而D-L和A-G的擬合曲線更平滑，這可能是因為S-G濾波器強調細節，而D-L擬合和A-G濾波方法重建了整體變化趨勢[10]。3種濾波方法NDVI時序重構效果的統計結果如表1所示，D-L擬合方法的對原始數據的增大比其他兩種方法更為明顯。標準差、協方差和相關系數反映了重構曲線與原始曲線的接近程度（即保真能力），標準差、協方差越小，相關系數越大，曲線的保真能力越好[10]。3種方法中，S-G濾波結果的標準差和均方根誤差均最小、相關系數最大，分別為0.191、0.024 2和0.993。與其他方法相比，S-G濾波器不僅能去除噪聲，而且能更好地擬合NDVI原始數據。

表1 三種濾波方法NDVI時序重構效果的統計參數Table 1 Statistical parameters of NDVI time series reconstruction effect of three filtering methods

2.2 物候參數提取

玉米的生長大致分為營養生長（種子發芽到抽雄期結束）和生殖生長（抽穗開始到成熟）兩個階段[24]。在對NDVI時序重構的基礎上，利用曲線曲率法和動態閾值法[25]分別提取恢復治理區和正常耕作區玉米的關鍵物候期（出苗期、拔節期、抽雄期、成熟期和物候期間隔時間）。

2.2.1 動態閾值法

利用J?nsson等[25]提出的動態閾值法，提取研究區玉米作物的出苗期和抽雄期。將作物生長曲線上升階段，距離NDVI最小值為NDVI最大值與最小值差值10%的位置定義成作物營養生長的開始期，即出苗期；在作物生長曲線下降階段，將距離NDVI最大值為NDVI最大值與最小值差值10%的位置定義為生殖生長的開始期，即抽雄期，如圖5a所示。

2.2.2 曲線曲率法

利用李艷等[24]、J?nsson等[25]使用的曲線曲率法，通過計算NDVI時序重構曲線曲率極值確定玉米的拔節期和成熟期，NDVI時序重構曲線的曲率通過計算重構曲線的一階導數獲得，曲率絕對值最大值對應時刻為農作物生長（衰老）速率變化最大的時期，重構曲線上升曲線曲率最大值點為玉米營養生長的拔節期；重構曲線下降曲線曲率最小值點為玉米生殖生長的成熟期，如圖5b。

2.2.3 物候期間隔時間

此次研究中，還選取了不同物候期之間的間隔天數作為物候指標，分析整個生長期內恢復治理區作物與正常作物生長期整體生長狀況差異。具體指標計算如下

式中DOY為關鍵物候期的年積日天數，L1、L2和L3分別指4個關鍵物候期（出苗期、拔節期、抽雄期和成熟期）之間的間隔天數。

2.3 植被指數

野外調查中發現，相比正常玉米作物，恢復治理區玉米作物明顯黃化。為對比各植被指數區分恢復治理區作物與正常作物、診斷礦山恢復治理區植被健康狀況的能力，選取6種常用于綠色植被色素檢測的經典植被指數，對比其對于區分恢復治區作物與正常耕作區作物和反應作物生長狀況的能力。6種植被指數分別為歸一化植被指數（NDVI）[26]，綠色歸一化植被指數（GNDVI）[27]，用以估算葉綠素含量的多時相中分辨率成像光譜儀（MERIS）的陸地葉綠素指數（MTCI）[28]，與葉綠素a、葉綠素b和類胡蘿卜素的單位面積濃度具有較強相關性的特征色素簡單比值指數（PSSRA）[29]，用以估算植被冠層葉綠素含量的倒紅邊葉綠素指數（IRECI）[30]和反應葉綠素含量變化的改進的葉綠素吸收指數（MCARI）[31]，如表2。

表2 植被指數及其計算公式Table 2 Vegetable indices and their formulas

2.4 植被紅邊

綠色植被在680至780 nm之間反射光譜中最大斜率對應的波長位置被定義為植被的紅邊位置（REP，Red-edge Position）[32]。REP對于植被葉片與冠層綠葉素含量差異、環境條件差異對葉綠素含量估算的影響不敏感，可用于評估植被的葉綠素含量、物候期生長狀況、健康狀況以及重金屬污染狀況等[32]。通過計算REP可以判斷植被冠層的生長期、病蟲害狀況以及重金屬含量等。為了從不同類型的光譜數據中提取REP參數，前人提出了多種REP提取方法，其中最常用的方法包括最大一階導數法、多項式擬合技術、逆高斯方法和線性外推技術等。

根據哨兵2號數據在紅邊區域的波段設置，本次研究中使用四點線性插補法[32]提取目標像元紅邊位置。該方法假設紅邊的反射率ρ曲線可以用一條位于780 nm和670 nm反射率中點附近的直線來近似，定義紅邊位置為

對于哨兵2號遙感數據，紅邊位置則定義為

2.5 植被覆蓋度

前人研究證實了植被覆蓋度和NDVI之間的顯著線性相關關系[4]，所以通過建立二者之間的轉換關系可以提取植覆蓋度信息。利用像元二分法模型對植被覆蓋度進行估算

式中FVC為植被覆蓋度，NDVIsoil為全裸土（無植被覆蓋）區域的歸一化植被指數值，NDVIveg為全植被覆蓋區域的歸一化植被指數值。研究中通常選取累積像元數占比5%和95%處的NDVI值作為純裸土（NDVIsoil）和純植被（NDVIveg）的NDVI值[33]。

3 結果與分析

3.1 不同區域作物物候特征提取與對比

利用曲線曲率法和動態閾值法在恢復治理區（A區和B區）和正常耕作區（C區）內各像元重構的NDVI曲線上分別提取了出苗期、拔節期、抽雄期和成熟期等4個玉米作物關鍵物候指標。A、B、C三區域各像元出苗期、拔節期、抽雄期和成熟期均值如圖6。可以看出，A區和B區玉米作物各關鍵物候期均晚于C區，C區平均出苗期為第147天，A區和B區平均出苗期分別為第159天和第152天；A區和B區平均拔節期分別為第195天和第192天，晚于C區的平均第183天；A、B、C三區域平均抽雄期比較接近，分別為第249天、第249天和第248天；A區和B區平均成熟期分別為第296天和第297天，晚于C區的平均第275天。顯示出恢復治理區與正常耕作區玉米作物生長狀況的差異。從3個區域各物候指標的統計參數（表3）來看，正常玉米種植區內部作物各物候指標標準差較小，物候指標較穩定，而恢復治理區內部作物各物候指標標準差較大，各像元間各物候指標變化明顯，顯示出復墾土地區域內作物生長狀況存在差異。另外，三區域各物候指標值的振幅區間存在較大重疊，不利于區分恢復治理區和正常耕作區作物。

表3 訓練區像元關鍵物候期的統計參數Table 3 The statistical parameters of the key phenological periods o the training areas

同時對比了訓練區像元各關鍵物候期間間隔的時間長度，如圖7。L1、L2和L3的統計參數如表4。整體而言，C區域拔節期與抽雄期間時間長度（L2）比A、B區域長，而A、B區域出苗期與拔節期間時間長度（L1）和抽雄期與成熟期間時間長度（L3）大于C區域。其中A、B區出苗期與拔節期間時間長度（L1）與C區相差較小，A、B區域平均出苗期與拔節期間時間長度（L1）分別比C區域分別長約1～2 d和3～4 d，這可能是恢復治理區土壤質地的差異對玉米作物生長前期影響較弱。隨著作物逐漸生長，土壤質地的差異也逐步顯現，表現為A、B區域平均L2比C區域分別短約10～11 d和7～8 d，抽雄期與成熟期間時間長度（L3）分別比C區域長20～21 d和21～22 d。相比之下，C區域L3參數與A、B區域相差更大。同時，統計參數（表4）也反映，3個區域內部像元L3標準差最小，且其L3取值振幅重疊度也最小，利于區分恢復治理區 和正常耕作區植被及其生長狀況。

表4 訓練區像元關鍵物候期間隔時間長度的統計參數Table 4 The standard deviation, maximum and minimum of the four phenological periods from the three training areas

3.2 不同區域作物植被指數對比

利用A、B、C三區域NDVI均值差異最大的第235天（2020年8月22日）哨兵2號數據計算A、B、C三區各像元NDVI、GNDVI、MTCI、PSSRA、IRECI、MCARI等植被指數，并對比恢復治理區（A區和B區）和正常耕作區（C區）像元各植被指數值的分布情況，如圖8。

結果顯示正常作物區像元NDVI、GNDVI、PSSRA和IRECI值振幅區間與恢復治理區無重疊，能較好區分恢復治理區與正常耕作區作物及其生長狀態差異，而恢復治理區與正常耕作區像元MTCI值和MCARI值分布區間有較大重疊，難以區分三區域作物及其生長狀態的差異。C區像元NDVI、GNDVI、PSSRA和IRECI值整體大于A區和B區，其中C區像元PSSRA值振幅區間與A區和B區差異最大，C區PSSRA值振幅下限（累積像元數占比5%處值）比A區和B區PSSRA值振幅上線（累積像元數占比95%處值）分別高12.61%和44.74%，C區NDVI、GNDVI和IRECI值振幅下限比A區振幅上線分別高1.94%，3.82%和5.81%，比B區振幅上線分別高7.20%、5.73%和42.59%。

3.3 不同區域作物紅邊位置對比

同樣利用2020年8月22日哨兵2號數據計算A、B、C三區各像元植被紅邊位置指數（REIP），A、B、C三區各像元REIP值分布如圖9。可見，相較正常健康作物（C區），恢復治理區作物（A、B）紅邊明顯向短波方向移動，作物光譜的“紅邊藍移”現象反映了恢復治理區作物與正常耕作區作物健康狀況的差異。從REIP值的振幅區間上看，正常耕作區（C區）像元REIP值較集中，分布在728.58 nm至729.93 nm之間，而恢復治理區像元REIP值較分散，A區REIP值分布于714.98 nm至726.22 nm之間，B區REIP值分布于716.90 nm至728.38 nm之間，同時，恢復治理區像元REIP值振幅區間與正常作物區無重疊，表明植被紅邊位置指數（REIP）能較好區分恢復治理區、正常耕作區作物及其生長狀態。C區REIP值振幅下限比A區和B區振幅上線分別高0.03%和0.32%。

3.4 不同區域作物植被覆蓋度對比

通過對研究區NDVI圖像像元值進行統計，確定累積像元數為5%和95%處的NDVI值，即純裸土（NDVIsoil）和純植被（NDVIveg）的NDVI值，分別為0.221和0.879，并依據公式6計算了A、B、C三區域的植被覆蓋度（FVC），三區域像元FVC值分布如圖10。結果顯示正常耕作區（C區）像元FVC分布十分集中（0.937至0.968），而恢復治理區（A區、B區）像元FVC分布比較分散，A區像元FVC值分布在0.376至0.914，B區像元FVC值分布在0.537至0.851。C區像元FVC值整體大于A區和B區，同時C區FVC值振幅區間與A區和B區沒有重疊，表明FVC能較好區分恢復治理區、正常耕作區作物及其生長狀態。C區FVC值振幅下限比A區、B區FVC值振幅上限分別高2.53%和10.17%。

4 討 論

4.1 不同區域作物物候參數差異分析

本研究利用哨兵2號影像時間序列和S-G濾波方法重構了研究區NDVI時間序列，并基于動態閾值、曲線曲率方法等方法提取并對比了礦山恢復治理區和正常耕作區玉米作物關鍵物候期和各物候期間時間間隔等關鍵物候參數，驗證了礦山恢復治理區作物與正常耕作區作物物候特征的差異。相較正常作物，恢復治理區作物出苗期、拔節期、抽雄期、成熟期不同程度推遲，這可能是礦山恢復治理區土壤質地、含水量差異，甚至重金屬脅迫減緩了玉米的生長速率，抑制了胚根的生長[34]。另外，恢復治理區與正常耕作區作物生長狀態的差異是逐步顯現的，表現為A、B、C區域像元各關鍵物候期間的時間間隔差異逐漸增大，這與前人對于重金屬脅迫下水稻物候特征變化的研究結論類似[10,15]。各物候參數中，抽雄期與成熟期間時間長度（L3）對于診斷恢復治理區與正常耕作區作物生長狀態效果最好，訓練區內部像元L3取值較穩定（標準差最小，表4），且A區、B區和C區的L3取值振幅重疊度也最小。同時，此次研究驗證了哨兵2號數據用于提取、分析作物物候信息的可行性，哨兵2號數據在可見光、近紅外波段空間分辨率為10 m和20 m，重返周期5 d，一定程度上能夠彌補Landsat數據（重返周期16 d）用于植被物候信息研究中時間分辨率不足以捕捉到快速發生的變化（如農田綠化）的問題[18]，以及MODIS、AVHRR和MERIS等數據用于非均質區域植被物候特征研究時的混合像元問題[17]。研究區周邊佳木斯、集賢農業氣象監測站點監測數據顯示，該區域春玉米平均拔節期與成熟期分別為第180天至190天和第265至275天[35]，與此次提取的正常耕作區玉米作物拔節期（平均第183天）和成熟期（平均第275天）數據大致相符，驗證了本文數據源和方法應用于玉米物候信息提取的可靠性。仍需要注意的是，由于部分時間段內（如第250天至275天，第275天至300天）影像云量過大無法使用，導致此次研究中影像時間序列間隔不一致，可能會影響物候信息提取的精度。下一步工作中，將嘗試采用時空融合算法，利用MODIS等數據彌補部分期次影像云量過大的問題。

4.2 遙感植被指標敏感性分析

研究結果表明，REIP、NDVI、GNDVI、PSSRA、IRECI和FVC等遙感指標能較好區分礦山恢復治理區與正常耕作區作物及其生長狀態差異。整體表現為，正常耕作區像元遙感指標值振幅區間相對收斂，而恢復治理區作物像元各遙感指標值振幅區間則比較發散（圖8-10；表5），表明正常耕作區作物生長狀態比較統一，而恢復治理區作物生長狀態有一定差異，反應了恢復治理區內作物養分生理要素狀況的差異。同時，正常耕作區像元REIP、NDVI、GNDVI、PSSRA、IRECI和FVC值整體大于恢復治理區。研究結果顯示NDVI能較好地識別恢復治理區植被及其生長狀況，與前人得出的NDVI能可靠評估植被健康的結論一致[36]。利用正常耕作區像元各遙感指標值振幅區間下限（累積像元數占比5%處值）高于恢復治理區像元各遙感指標值振幅區間上限（累積像元數占比95%處值）的百分比來定量對比各遙感指標對于區分恢復治理區作物及其生長狀況的敏感程度（如圖11），同時計算研究區像元各遙感指標值的標準差（表5）。結果顯示C區像元PSSRA值振幅區間下限高于A區和B區PSSRA值振幅區間上限的百分比最高，分別為12.61%和44.72%，表明PSSRA指數對于區分恢復治理區和正常耕作區作物及其生長狀況差異最敏感。各遙感指標中，恢復治理區像元的GNDVI值標準差最小（表5），表明恢復治理區像元的GNDVI值范圍更集中，對于反應恢復治理區作物的整體生長狀況更穩定；而正常耕作區像元的NDVI值標準差最小（表 5），表明NDVI對于反應正常耕作區作物的整體生長狀況更穩定。

表5 訓練區像元各遙感指標統計Table 5 Statistics of the remote sensing indexes from the three training areas

4.3 不同區域作物養分生理信息差異分析

不同類型區作物遙感物候參數、植被指數等指標差異反應了作物光譜特征的差異，而這種差異實質上是作物養分生理要素的差異造成的。研究結果表明，正常耕作區像元REIP、NDVI、GNDVI、PSSRA、IRECI和FVC值整體大于恢復治理區。前人研究顯示NDVI、GNDVI、PSSRA等遙感指標與植被健康狀況、葉片葉綠素含量呈顯著正相關[27,29,36]，表明恢復治理區作物葉綠素含量小于正常耕作區作物，這與野外驗證的情況相符（恢復治理區玉米作物葉片明顯黃化）。

另外，研究表明GNDVI對于反應恢復治理區作物的整體生長狀況更穩定，而PSSRA指數對于區分恢復治理區和正常耕作區作物及其生長狀態差異最敏感。前人證實了PSSRA、GNDVI與玉米作物綠色葉面積指數的強相關性[37-38]，Nguy-Robertson等通過對多個玉米種植區實驗研究，發現當玉米作物綠色葉面積指數大于2m2/m2時，PSSRA與玉米作物綠色葉面積指數的相關性強于DNVI、EVI、SAVI（Soil-Adjusted Vegetation Index）等11種植被指數[38]。王來剛等利用無人機多光譜數據進行玉米葉面積指數反演，發現GNDVI與玉米作物綠色葉面積指數的相關性強于DNVI、EVI、SAVI、OSAVI（Optimized Soil-Adjusted Vegetation Index）4種植被指數[37]。所以恢復治理區與正常耕作區玉米作物光譜特征、PSSRA、GNDVI值的差異可能與兩區域玉米作物葉面積指數的差異有關。野外調查也證實了這一點，恢復治理區玉米作物株高和葉片大小明顯小于正常耕作區，且葉片黃化明顯，導致恢復治理區玉米葉面積指數與綠色葉面積指數小于正常耕作區玉米。

同時，恢復治理區表層堆墊土土壤容重的升高會影響植被扎根和水分的吸收。NDMI（Normalized Difference Moisture Index）與冠層的水分高度相關，與NDVI相比，NDMI可以更緊密地跟蹤植物生物量和水分脅迫的變化，多項研究證實了NDMI在非均質環境中描繪植被葉片含水量的能力[22,39]。為驗證恢復治理區與正常耕作區玉米葉片含水量的差異，對比了A，B，C區像元NDMI值的分布情況（圖12），同時分析了各色素相關植被指數與NDMI的相關性（表 6，圖 13）。圖12顯示，正常耕作區玉米作物NDMI值明顯高于恢復治理區玉米作物，表明正常耕作區玉米作物葉片含水量大于恢復治理區玉米作物，這是由于恢復治理區表層堆墊土土壤容重的升高影響了玉米作物水分的吸收。另外表6與圖13顯示，各與色素相關的植被指數中，PSSRA與NDMI相關性最強（R2=0.778），表明相較其他植被指數，PSSRA能更好的反應作物葉片含水量的差異。

表6 植被指數與NDMI統計關系Table 6 Statistical correlation between vegetation indices and NDMI

綜上，恢復治理區與正常耕作區玉米作物光譜特征、各遙感指標值的差異可能受兩區域玉米作物葉綠素含量、葉面積指數和葉片含水量等生理信息差異的綜合影響。另外礦山開采引起的水土污染和重金屬富集也可能導致植被的光譜脅迫特征，下一步工作中進一步補充獲取土壤質量、重金屬元素數據，分析其對玉米作物光譜的影響。本次研究方法是通過檢測因作物養分生理信息差異引起的光譜特征差異，來區分恢復治理區作物與健康作物，并評估作物的健康狀況，可用于檢測其他礦山恢復治理區域內植被生理信息差異引起的光譜特征差異。但由于其他礦山恢復治理區影響植被養分生理信息差異的具體類型（葉綠素類、氮素、水分、雜草、重金屬等）可能不同，需要對各遙感診斷指標的敏感性進行再分析。

5 結 論

基于哨兵2號遙感數據，分別提取了恢復治理區和正常耕作區作物遙感物候指標、各類型植被指數、紅邊位置、植被覆蓋度等遙感指標，并評估了各遙感指標對于區分正常耕作區與礦山恢復治理區作物，診斷作物生長狀態和健康程度的可行性，主要取得了以下結論：

1）研究區礦山恢復治理區與正常耕作區作物生長狀況的差異明顯。恢復治理區玉米作物出苗期、拔節期、抽雄期、成熟期等關鍵物候期不同程度推遲。同時，恢復治理區與正常耕作區作物生長狀態的差異隨著作物生長階段的推移逐步顯現，表現為恢復治理區與正常耕作區作物像元各關鍵物候期間的間隔天數差異逐漸增大。

2）對比評估了紅邊位置指數（REIP）、植被覆蓋度（FVC）、各類型植被指數等遙感指標對于區分正常耕作區與礦山恢復治理區作物，診斷作物健康程度的能力。REIP、NDVI、綠色歸一化植被指數（GNDVI）、特征色素簡單比值指數（PSSRA）、倒紅邊葉綠素指數（IRECI）和FVC等遙感指標能較好區分礦山恢復治理區與正常耕作區作物及其健康狀況差異，其中恢復治理區像元的GNDVI值標準差最小，對于反應恢復治理區作物的整體生長狀況更穩定。恢復治理區與正常耕作區像元的PSSRA值振幅區間差異最大，PSSRA指數對于區分恢復治理區和正常耕作區作物及其健康狀況差異最敏感。

3）通過野外調查結合遙感分析，發現恢復治理區與正常耕作區玉米作物光譜特征、各遙感指標值的差異可能受兩區域玉米作物葉綠素含量、葉面積指數和葉片含水量等生理信息差異的綜合影響。

本次研究驗證、評估了各遙感指標對于診斷礦山恢復治理區植被健康狀況的可行性，為利用遙感技術進行礦山恢復治理效果評估提供了一種技術思路。下一步嘗試將本文方法應用于區域性礦山恢復治理區植被健康狀況的評估，進一步驗證遙感診斷指標對其他礦山恢復治理區域植被健康狀況評估的適用性。