Google Earth Engine平臺支持下的鐵礦區開采及植被變化遙感動態監測*

劉家興，隋翔宇，包妮沙，毛亞純，趙占國，石玉君

(1.東北大學 資源與土木工程學院，遼寧 沈陽 110819;2.中國黃金集團有限公司，北京 100011；3.中國黃金集團內蒙古礦業有限公司，內蒙古 滿洲里 021400)

中國鐵礦石90%以上以露天形式開采，具有開采范圍大、服務年限長、對周邊環境影響大等特點。遙感對地觀測技術能快速、動態、大范圍地獲取礦區地表覆蓋信息，具有周期短、實時性、動態性和信息豐富等特點，可實現露天開采過程中對礦區以及周邊環境陸面演變信息的動態提取[1]，Landsat遙感衛星數據由于數據量豐富、空間分辨率適中，在礦區環境監測中被廣泛應用。利用目視解譯[2]、植被歸一化指數[3]、像元二分模型[4]、圖像分類方法[5-8]等可以獲取礦區土地利用/覆蓋變化信息，進一步反應礦區開采對周邊生態環境的影響以及復墾效果。其中，遙感數據的時間跨度和間隔是否能夠覆蓋礦區開采的全生命周期，并且能夠全面反映礦區邊開采、邊復墾的地表變化信息，是礦區環境動態監測的關鍵。

Google Earth Engine(GEE)是Google提供的對大量全球尺度地球科學資料(尤其是衛星數據)進行在線可視化計算分析處理的平臺，可容納數PB的衛星影像和地理空間數據集目錄[9]。相比傳統的遙感處理軟件它采取云端計算的方法在谷歌云上運算，處理能力不受空間、時間的限制，可快速、批量的處理大量的衛星影像，具有簡單且功能強大的API，通過Web的交互式開發平臺使用Google提供的函數庫即可訪問云端數據并進行云計算。GEE平臺提供的最大合成植被歸一化指數(NDVI)、地表反射率數據可以為不同尺度的植被覆蓋度監測提供長時間、高精度的數據,包括亞熱帶地區林地年際變化監測[10]、黃土高原區植被覆蓋度時空變化監測[11]以及城市綠地動態監測[12]等。

本文在國家綠色礦山建設及東北老工業國土空間修復背景下，針對礦區開采以及恢復對土地利用的影響，基于GEE云平臺，選擇東鞍山鐵礦區和大孤山鐵礦區為研究區，以Landsat衛星1999～2019年數據為數據源，利用時間序列法監測采礦區周邊3 km緩沖區內的植被變化，通過NDVI時間變化以及變化速率空間分布，分析采礦活動下植被時空分布特征，利用梯度分水嶺的分割算法提取礦區地表開采破壞信息，獲得礦區近20年跨度的地表開采破壞土地面積變化，實現實時準確掌握一定時間內的礦區開采破壞土地、影響范圍內的植被時空變化特征，可以為東北老工業基地的生態文明建設、礦山土地復墾及生態恢復提供支持。

1 研究區概況與數據源

1.1 研究區概況

本文所選的研究區為位于鞍山市境內的東鞍山鐵礦區和大孤山鐵礦區，如圖 1所示。東鞍山礦區位于鞍山市南郊距市中心7 km處，至今已有五十余年的開采歷史，礦產資源儲量豐富，主要礦產品為赤鐵礦，平均品位32.46%，已累計輸出鐵礦石2.59億t，礦區作業實施土地復墾，遵循“邊開采、邊復墾”策略。大孤山礦區位于鞍山市中心東南方向12 km處，為亞洲第一大露天鐵礦，至今已有一百多年開采歷史，海拔由原地貌260 m采掘至-306 m，成為大型深坑露天鐵礦。礦石礦物種類主要為磁鐵礦、赤鐵礦，礦石品位為33.63%，儲量2.65億t。研究區植被類型以速生落葉喬木為主，包括刺槐、棉槐和火炬樹等。

圖1 研究區位置圖

1.2 GEE云平臺遙感數據

本次實驗所用數據如表 1所示，其中,Landsat 7系列的植被歸一化(NDVI)年度平均數據和年度最大合成數據用于本次研究礦區的植被變化，經過大氣校正的Landsat5、Landsat7、Landsat8的表面反射率數據用于本次研究礦區的開采變化。

表1 實驗所用數據

Landsat7年度平均NDVI合成數據集(LANDSAT/LE07/C01/T1_ANNUAL_NDVI)中的“NDVI”波段數據，為一年的第一天到最后一天，將每個像元每年的所有影像中取NDVI平均值作為該年的數據而合成。Landsat7年度最大NDVI合成數據集(LANDSAT/LE07/C01/T1_ANNUAL_GREENEST_TOA)中的“greenness”波段數據，為一年的第一天到最后一天，將每年的所有影像中最綠的像元作為合成值合成該年的數據。

2 研究方法

2.1 基于時間序列的礦區植被變化檢測方法

時間序列法是一種統計分析方法，針對某種問題或現象獲取一段時間內一系列的觀測數據，這些數據間往往具有一定的變化規律，再通過一定的數學方法計算出相應的函數值，以此來分析某一問題或現象內部存在的規律。在遙感中時間序列法利用一段時間內一系列的影像數據，根據影像像元灰度值在這段時間內的變化，對每個像元生成的一元線性回歸方程進行擬合，通過分析擬合的一元線性回歸方程的斜率正負及大小來分析該像元以及影像在該時間段內的變化以及變化速率。

本文對于研究區植被的變化檢測，采用的指數是NDVI歸一化植被指數，選擇Landsat7系列1999～2019年年度平均NDVI數據與年度最大NDVI數據，以東鞍山礦區和大孤山礦區開采區為中心，選擇3 km緩沖區分析植被時空變化特征，利用時間序列法，以年度最大NDVI數據構建代表每個像元NDVI變化速率(slope)的一元線性回歸模型，并對其進行顯著性檢驗(p<0.05)。通過顯著性檢驗后的slope可以反映該時間序列內植被的時空變化特征，slope為正的部分依據值區間平均分為兩個等級，分別為快速增長和緩慢增長；slope為負的部分依據值區間平均分為兩個等級，分別為緩慢降低和快速降低。

2.2 基于梯度分水嶺的分割算法的礦區開采范圍提取

分水嶺分割算法是一種基于拓撲理論的數學形態學的分割方法，它把圖像看作是一個拓撲地貌，把圖像的每一個像素的灰度值看作該點的海拔高度，把每一個局部極小值稱為集水盆地，分水嶺即集水盆地的邊界。

梯度分割即計算出該圖像中不同像元灰度值所代表的梯度，對于二維圖像f(x,y)來說，它的梯度為：

(1)

最簡單的梯度近似表達式如下：

Gx=f(x,y)-f(x-1,y)

(2)

Gy=f(x,y)-f(x,y-1)

(3)

式中，Gx為f(x,y)對f(x)的偏微分；Gy為f(x,y)對f(y)的偏微分。表示了圖像在(x,y)處x方向和y方向上的梯度，可以看出，圖像的梯度等于2個相鄰像素之間的差值。

本文梯度分割的圖像選擇基于鐵礦石與其他地物光譜的差異而構建出的歸一化差異鐵指數 (Normalized Difference Iron Index，NDII)[13]：

(4)

式中，Rband5和Rband6分別是landsat8衛星band5和band6處的反射率。對于landsat5數據和landsat7數據，本文同樣使用短波紅外和近紅外波段，即band5和band4波段。根據梯度變化較大的區域，找出適合礦區開采破壞土地的分割邊界，將圖像按照此邊界進行分割。 最后采用基于最優分割尺度下的分類精度的方法計算參考樣本與提取對象的相似度，來定量評價最終提取結果的精度。相似度公式如式(5)[14]所示：

(5)

式中，f為一種對象特征；C為分類對象，即本文中提取的礦區；R為參照樣本，即Google Earth中的高分辨率影像；C∩R為提取影像和高分辨率影像共同擁有的礦區；C-R為屬于提取后影像中的礦區但不屬于高分辨率影像的礦區；R-C為屬于高分辨率影像但不屬于提取后影像的礦區。

3 結果與分析

3.1 研究區植被變化

如圖2(a)所示，為東鞍山礦區年度NDVI變化折線圖，1999～2019年期間東鞍山礦區年度最大NDVI總體呈現上升趨勢，1999年年度最大NDVI值為0.334, 2019年年度最大NDVI值為0.454。其中，在2010年NDVI值出現較大的波動，降低為0.195，主要是由于當年礦區西北部的部分農用地轉化為建設用地，以及這一期間礦區開采在東南方向形成大面積排土場。隨著礦區邊開采邊復墾工程的進行，排土場在2011年之后NDVI值逐漸上升，并在2016年礦區年度最大NDVI達到最高值0.455，說明復墾效果較好。根據礦區NDVI變化速率空間分布結果，如圖2(b)所示，開采區周邊大部分植被變化區域通過顯著性檢驗，而且緩慢增長面積占總變化面積的67.4%，說明盡管礦區開采對周邊環境產生一定的影響，但是隨著生態恢復以及土地復墾工作的開展，區域植被覆蓋度逐漸增加，復墾工作初見成效；快速降低面積占總變化面積的8.8%，主要是由于西北部及西南部城市擴張占用農用地以及礦區東南排土場堆砌致使植被銳減導致的；緩慢降低面積占總變化面積的19.0%，除分布在快速降低區域周邊外還分布在開采區域北部，主要由于邊開采邊復墾工作進行的同時，復墾工作出現一定的滯后性和恢復期。

圖2 東鞍山植被變化數據

如圖3(a)所示，為大孤山礦區年度NDVI變化折線圖，1999～2019年期間大孤山礦區年度最大NDVI總體呈現上升趨勢，1999年年度最大NDVI值為0.360, 2019年年度最大NDVI值為0.447。數據方差為0.112 8，相比于東鞍山NDVI方差0.040 5，變化波動較大，2017年由于年均降水量較低，出現NDVI較低值0.349，之后NDVI指數一直在緩慢上升。如圖3(b)所示，為大孤山礦區經過顯著性檢驗的植被變化空間分布結果，植被顯著減少區域占總變化面積的7.3%，主要分布在因城市擴張而占用農用地的東北部、尾礦庫周邊以及由于礦區開采而形成的未復墾的排土場和剝離區等。植被增加面積占總變化面積的74.9%，其中，礦區南部的排土場NDVI一直呈現緩慢增長趨勢，且有部分區域為快速增長，同樣說明1999～2019年間排土場的復墾工作已經初見成效。

圖3 大孤山植被變化數據

3.2 鐵礦區開采范圍變化

如圖4、圖5所示，分別為東鞍山礦區和大孤山礦區的梯度分割結果與Google Earth目視解譯結果的對比。參考2019年Google Earth影像目視解譯結果對梯度分割算法提取的礦區邊界進行精度評定，利用梯度分割算法自動提取的東鞍山礦區提取精度達到85.60%，大孤山礦區提取精度達到86.99%。

圖4 東鞍山目視解譯及其分割結果

圖5 大孤山目視解譯及其分割結果

如圖6所示， 東鞍山礦區面積整體呈現減小趨勢。1999～2007年間，東鞍山礦區開采破壞面積不斷減小，在這幾年間，由于鋼鐵工業的快速發展和政府政策的導向造成了鐵礦石需求量增大，鐵礦石價格升高，中國鐵礦資源的開采在這幾年明顯加快了速度，鐵礦石產量和年消耗成品礦數量大幅增加，以滿足鋼鐵工業增長的需求。而在2007～2008年，這一年間的開采面積從1.104 3 km2增加到1.443 6 km2，原因是在該年鞍鋼集團擴大露天礦開采面積，在東鞍山處主要擴大區域為礦山的東北部。自2008年起，礦區的面積不斷減小至0.997 4 km2，這是由于此時東鞍山鐵礦加強“邊開采，邊綠化”的措施以及鞍山市發布青山工程規劃[15]。自2008年起，政府對東鞍山附近露天采礦場以及多個排土場進行了復墾、綠化，因此礦區面積在近幾年間不斷減小。

圖6 東鞍山礦區變化監測

如圖7所示，大孤山礦區面積整體呈現出比東鞍山礦區面積更大的減小幅度。1999～2007年間，大孤山礦山面積整體呈減小趨勢。2002年相比2001年礦區面積有大幅度的增長，這是因為在2002年，大孤山西南部排土場的增加。2007～2009年間，為了恢復大孤山生產能力，延長礦山壽命，大孤山鐵礦北部擴大工程正式開工，這一擴建工程使這幾年間的礦區面積增加了1.176 3 km2。在2009年后，政府實施“邊開采，邊復墾”的政策，恢復大孤山地區綠化程度以及礦區生態環境，礦區開采面積減少了2.883 6 km2。

圖7 大孤山礦區變化監測

4 結論與討論

本文通過利用Google Earth Engine(GEE)云平臺，使用Landsat系列數據對位于鞍山市鐵礦群的東鞍山鐵礦區和大孤山鐵礦區進行了近20年的植被覆蓋變化監測以及礦區開采面積變化監測，通過分析數據及結果得出以下結論：

(1)基于NDVI時間序列能夠從時空反映露天開采植被的變化，1999～2019年間，東鞍山、大孤山礦區及影響范圍內NDVI指數整體呈增加趨勢，尤其礦區排土場由于進行土地復墾和植被恢復，NDVI值呈現顯著增加，說明礦區復墾工作初見成效。

(2)基于歸一化差值鐵指數的梯度分割算法，能夠準確、快速、自動提取礦區開采破壞面積，1999～2019年間，在綠色礦山建設以及開采邊界優化的基礎上，礦區開采破壞面積整體不斷縮減，同時受到鐵礦石價格和需求影響，開采破壞面積在這一期間有輕微波動和起伏。

(3)利用GEE云平臺強大的數據庫及其數據處理功能，實現了對鞍山市東鞍山鐵礦區和大孤山鐵礦區的植被變化以及礦區開采面積變化的動態監測，突破了本地硬件條件的限制，實現在線數據處理及分析，且大量節省了查找數據與下載數據的時間。