基于Landsat8 與高分數據的礦山植被動態監測研究

吳鳳敏，鄭稚棚，余 靜，梁 星，張靈犀，楊光譜，余 洋，趙珍妮

（1. 重慶市地理信息和遙感應用中心，重慶 401147）

目前，對于植被遙感動態監測有大量研究[1-9]，其中Landsat 數據應用最為廣泛，該類數據可以免費下載，且適用于大尺度范圍的長時間序列監測[10-22]。基于此，研究以重慶市中梁山片區為實驗區，利用Landsat8 數據對礦山植被長時間序列情況進行分析，同時結合高分數據對礦山植被短時間序列進行動態監測，并對于不同分辨率數據的植被監測結果及關聯性進行分析，希望能夠為后續礦山修復治理過程中植被狀況動態監管提供技術支撐。

1 研究區概況及研究方法

1.1 研究區概況

中梁山片區地理坐標為106°23′20″～106°24′25″E，29°26′10″～29°31′50″N，總面積為28.63 km2，最高點位于中部山頂，最高高程為695.8 m，最低點位于山腳，最低高程為270 m。研究區內礦山有70個，總面積為191.36 hm2，位于中梁山背斜之南傾伏端，以低山地貌為主。研究選擇中梁山片區為示范區域，主要由于該區域露天礦山數量較多，存在已完成修復治理、正在修復治理和關閉未治理等多種類型，在礦山植被狀況動態監測研究中具有一定代表性（圖1）。

圖1 中梁山片區研究區域

研究利用landsat8 數據對2012—2021 年中梁山片區露天礦山植被指數進行變化分析，利用高分數據進行2017—2021 年短時間序列變化分析。Landsat8 數據來源于地理空間數據云，高分數據包括高分1號和高分6 號數據，來源于重慶市地理信息和遙感應用中心。為盡量保證植被監測年度生長一致，Landsat8數據主要采用夏天時期影像， 數據獲取時間分別為2013-08-19、2014-07-21、2015-07-08、2016-07-11、2017-08-14、2018-05-13、2019-03-29、2020-09-07、2021-06-06，根據1∶5 000DEM數據計算該區域平均高程為0.37 km。高分數據數據獲取從2017—2021年，部分數據因為不完全覆蓋實驗區需要進行拼接，每景數據平均高程同樣根據1∶5 000DEM進行計算（表1）。

表1 高分衛星數據使用總體情況表

1.2 主要技術流程

研究采用ENVI5.3 軟件首先對Landsat8 數據和高分數據進行預處理，分別進行植被指數計算，并對礦山植被狀況進行現狀及變化分析，最后對兩者的關聯性進行詳細對比分析（圖2）。研究選擇NDVI 作為礦山植被監測指數，由于礦山生態修復過程中植被生長變化較大，NDVI 能夠較好地反映植被指數從低值到高值區間的動態變化，與植被分布密度呈線性相關。數據預處理主要包括輻射定標、大氣校正、正射校正等流程。輻射定標處理主要指將影像的灰度（D）值轉化成輻射亮度值或表觀反射率；大氣校正主要是將輻射定標后的表觀反射率或輻射亮度值轉化為地表反射率的過程；正射校正需基于地形高程模型（DEM）來對影像中所有的像元執行地形變形的校正，讓圖像滿足正射投影的要求，可根據影像自帶的RPC參數和ENVI中的DEM信息，對影像數據進行正射校正。

圖2 Landsat8數據和高分數據處理及關聯性分析技術流程

采用ENVI軟件進行輻射定標時，Landsat8數據可以直接采用自帶模型Radiometric Calibration進行輻射定標，高分數據需要網上下載并安裝China Satellites Support插件，同時根據中國資源衛星應用中心發布的輻射定標參數（http://www.cresda.com/CN/index.shtml）。由于NDVI計算取值范圍為[-1.0，1.0]，研究對該范圍外的作為異常值進行剔除。根據已有學者研究成果，裸土NDVI 理論上接近于0，但由于受眾多因素影響，其取值范圍在-0.1～0.2之間。

2 Landsat植被動態監測

2.1 現狀情況分析

參考已有研究，按照無植被覆蓋區[-1.0，0.2），低植被覆蓋區[0.2，0.3），中等植被覆蓋區[0.3，0.4），高植被覆蓋區[0.4，1.0）進行分級。結果顯示，2021年中梁山片區礦山高植被覆蓋區面積為81.09 hm2，占中梁山片區礦山總面積的42.38%；中等植被覆蓋區面積為29.35 hm2，占比15.34%；低植被覆蓋區面積為30.72%，占比16.05%，表明中梁山片區礦山植被整體情況較好，但存在部分礦山還在修復治理過程中，植被指數較小或者為推土情況無植被覆蓋（表2）。高植被覆蓋區中，NDVI 以0.4～0.5 區間為主，該區間礦山面積占高植被覆蓋區面積的50%以上，＞0.6的面積占比為46.88%。以單個礦山來看，中梁山片區70 個礦山中，平均NDVI 超過0.4 的礦山數量有36 個，平均NDVI 為0.3～0.4 的礦山為19 個，平均NDVI 為中等的礦山共有55 個，說明隨著礦山修復治理工作的不斷開展，礦山的植被指數狀況總體較好，大部分礦山植被恢復較好，部分礦山由于正處于修復治理過程中，可能存在部分區域推填土或者剛種植植被，植被指數較低。

表2 2013—2021年中梁山片區礦山植被指數NDVI分級表

2.2 年際變化分析

研究以每隔1 a NDVI值對礦山范圍內數據進行對比分析，2013—2021年中梁山片區植被指數NDVI呈現先減小后增加趨勢（圖3、4），主要波谷在2017年左右，由于礦山開采時期植被破壞較大，植被指數較低，在修復治理過程中，也可能存在部分區域推填土破壞植被的情況；隨著礦山不斷修復，前期的植被指數從較低的覆蓋指數逐漸增加，最終修復完成后植被不斷恢復，植被指數不斷增加，表明礦山生態修復治理取得了較好的成績。其中，礦山高植被覆蓋區（NVDI＞0.4）從2013 年的104.81 hm2下降至2017 年的42.17 hm2，減少面積為62.64 hm2；從2017—2021 年面積又逐漸增加，2021年的高植被覆蓋區面積達到81.09 hm2。低植被覆蓋區（0.2＜NDVI≤0.3） 具有同樣的規律，2013—2017 年面積減少12.64 hm2，2017 年到2021 年面積增加17.54 hm2。2017—2021年間高植被覆蓋區增加的趨勢較低植被覆蓋區更大，說明這個期間為礦山植被恢復的主要時期，最開始以低植被覆蓋的治理初期隨著時間的推移，植被逐漸變好，并且變為高植被覆蓋區較多。

圖3 2013年和2021年中梁山片區植被指數NDVI分布圖

研究發現，存在個別礦山植被指數到2021年總體情況較低，比如中梁山片區萌特礦山，經實地調查，該礦山于2019年6月開始工程治理，目前正處于修復治理過程中，項目周期為66 個月，到2025 年為止，存在部分區域因修復治理過程中對周邊土地的踩踏以及修復治理最開始種植小樹苗在影像上不明顯，因此總體上植被指數較低。

圖4 中梁山片區礦山2013—2021年植被指數變化情況

2.3 變化斜率分析

研究采用一元線性回歸方法模擬每個柵格的NDVI變化趨勢，重點分析礦山范圍內不同修復治理區域植被變化情況。對于NDVI 時間序列數據，每個像元對應有若干年的時間序列數值，對這些數值進行線性擬合，獲取直線的斜率，像元斜率表明了在該時間序列中該像素所代表的植被指數的演化趨勢。主要計算公式如下：

式中，k為斜率；xi為具體年份；yi為每年的NDVI值。

2017—2021 年中梁山片區礦山范圍內斜率＞0 的像元總面積為85.15 hm2，占礦山總面積的44.49%，表明這些區域內植被逐漸變好；斜率＜0 的像元面積為106.21 hm2，表明這些區域植被狀況有變差的趨勢。參考已有研究[13]，將NDVI 變化程度分為7 級（表3），圖5 為以30 m×30 m 像元大小，統計中梁山片區礦山范圍內所有像元的斜率值分布情況。中梁山片區礦山范圍內像元明顯改善面積有55.22 hm2，占像元總面積的28.85%；輕微改善及中度改善的面積合計為26.14 hm2，占比13.66%，存在一般以上斜率為負值，植被情況較差的情況，主要原因是中梁山片區面積最大的礦山重慶萌特礦山處于修復治理過程中。雖然有一部分區域在種植樹木，但由于數據還在恢復中，植被指數變化較小，且治理過程中存在部分推填土情況，因此植被指數總體變小。

表3 2013—2021年中梁山片區礦山像元NDVI斜率分級（Landsat8數據）

圖5 中梁山片區礦山范圍內單個像元斜率分布情況

研究通過ArcGIS 工具將計算斜率值與中梁山片區礦山像元進行空間關聯，可以看出中梁山片區中部及東部礦山修復情況較好，植被斜率大部分為正數，右下角萌特礦山和中部部分礦山目前還處于修復治理過程中或者關閉未治理情況，植被覆蓋較差，如圖6所示。

圖6 2013—2021年中梁山片區礦山范圍內像元植被指數斜率分布圖

3 高分數據植被動態監測

3.1 現狀情況分析

研究按照與Landsat8 數據同樣分析方法對高分數據計算的植被指數進行分析，由于分辨率較高（8 m），礦山植被指數監測結果較landsat8 數據更精細。中梁山片區整體植被情況較好，大部分區域均為高植被覆蓋區域（NDVI ＞0.4），70個礦山總體上位于中高植被覆蓋區域范圍內，2021 年有植被覆蓋區域面積為150.49 hm2，其中高植被覆蓋區域94.84 hm2，占礦山總面積的49.56%。與landsat8 計算結果相比，高植被覆蓋區面積較大，說明高分數據對于植被指數的響應更敏感。中等植被覆蓋區面積為29.48 hm2，占比為15.41%（圖7）。

圖7 2021年中梁山片區高分數據NDVI監測分布圖

3.2 變化斜率分析

由于高分數據只有2017—2021 年5 a 的數據，因此進行短時間序列斜率變化分析（表4，圖8），圖8為以8 m×8 m像元大小，統計中梁山片區礦山范圍內所有像元的斜率值分布情況。結果表明中梁山片區礦山范圍內NDVI斜率＞0的像元面積為78.94 hm2，占礦山總面積的41.25%；斜率＜0的像元面積為112.42 hm2，總體上超過礦山總面積的一半，NDVI 斜率的分布情況與Landsat8 數據分析結論基本一致，僅部分分級數值略有差異。

圖8 中梁山片區礦山范圍內單個像元斜率分布情況（高分數據）

表4 2017—2021年中梁山片區礦山像元NDVI斜率分級（GF數據）

從分級情況看，礦山范圍內嚴重退化區域面積最大，為108.51 hm2，主要是由于中梁山片區萌特礦山正在修復治理過程中造成。明顯改善區域面積居第二，占礦山總面積的39.45%，表明其他大部分礦山區域的植被都是變好的，說明礦山修復治理取得了良好的效果。

4 植被指數關聯性分析

研究選擇2019 年數據進行關聯性分析，利用Landsat8 （2019-03-29） 和高分數據（2019-07-28）分別計算NDVI，得到兩者空間分布圖。結果顯示中梁山片區范圍內高分數據在高植被覆蓋區域分級粒度更細，NDVI 大于0.8 在圖上分布較為明顯，小于0.2區域面積總體較少，主要由于高分數據分辨率為8 m，而Landsat8分辨率為30 m，在像元尺度上部分高值和低值區域進行了一定綜合，因此在中間部分值的像元數量較多（圖9，表5）。

表5 中梁山片區礦山2019年Landsat8和高分數據計算NDVI分級情況對比

圖9 2019年Landsat8和高分數據計算NDVI空間分布圖對比

通過對中梁山片區礦山范圍內的兩種影像計算的DNVI 進行計算得出，高分數據有植被覆蓋區域（NDVI ＞0.2）為82.68 hm2，明顯小于Landsat8計算出的有植被覆蓋區域面積（132.93 hm2），分析造成趨勢情況不一致的原因有以下幾種：一是分辨率原因導致部分Landsat8 像元綜合了周圍的NDVI 值較高的植被區域導致像元內植被指數整體較高；二是由于兩類NDVI數據在進行正射糾正時均采用ENVI自帶的DEM數據，可能存在部分像元位置有一定偏移；三是兩類NDVI 數據的時相不一致，對于植被生長期情況的計算結果本身也有一定差異。

研究利用中梁山片區范圍對2 種數據進行植被指數關聯性分析，通過ArcGIS 中Spatial Analyst Tools 工具進行計算。計算結果顯示，通過Landsat8 和高分數據計算NDVI 相關系數為0.739 3，表明2 種遙感影像在計算結果上呈現較高的正相關性。同時，研究以Landsat8 像元為網格（30 m × 30 m），計算高分數據NDVI 平均值，分析空間上同一像元NDVI 值差異情況，結果表明，高分數據與Landsat8 數據計算NDVI值呈現花瓣狀相關特征，兩者在高值區和低值區呈聚攏狀態，在中間值區呈離散狀態，總體上表現為中間寬兩頭窄的花瓣形。以線性關系式表達兩者特征，決定系數R2=0.692，表明兩者雖然對于同一像元的值存在差異，但總體上還是具有較好的線性相關性（圖10）。

圖10 2019年Landsat8和高分數據計算NDVI相關性分析

5 結論與討論

通過對Landsat8 數據和高分數據進行預處理，對礦山修復治理過程中植被狀況進行時間序列的動態分析，并對2 種數據計算處理的植被指數進行關聯性分析，希望通過礦山植被遙感監測能夠對未來礦山修復治理監管提供基礎支撐。本研究仍然存在一定不足之處，需要進一步探討。