Google Earth Engine平臺(tái)支持下的鐵礦區(qū)開(kāi)采及植被變化遙感動(dòng)態(tài)監(jiān)測(cè)*

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

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

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

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

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

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

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

圖1 研究區(qū)位置圖

1.2 GEE云平臺(tái)遙感數(shù)據(jù)

本次實(shí)驗(yàn)所用數(shù)據(jù)如表 1所示，其中,Landsat 7系列的植被歸一化(NDVI)年度平均數(shù)據(jù)和年度最大合成數(shù)據(jù)用于本次研究礦區(qū)的植被變化，經(jīng)過(guò)大氣校正的Landsat5、Landsat7、Landsat8的表面反射率數(shù)據(jù)用于本次研究礦區(qū)的開(kāi)采變化。

表1 實(shí)驗(yàn)所用數(shù)據(jù)

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

2 研究方法

2.1 基于時(shí)間序列的礦區(qū)植被變化檢測(cè)方法

時(shí)間序列法是一種統(tǒng)計(jì)分析方法，針對(duì)某種問(wèn)題或現(xiàn)象獲取一段時(shí)間內(nèi)一系列的觀測(cè)數(shù)據(jù)，這些數(shù)據(jù)間往往具有一定的變化規(guī)律，再通過(guò)一定的數(shù)學(xué)方法計(jì)算出相應(yīng)的函數(shù)值，以此來(lái)分析某一問(wèn)題或現(xiàn)象內(nèi)部存在的規(guī)律。在遙感中時(shí)間序列法利用一段時(shí)間內(nèi)一系列的影像數(shù)據(jù)，根據(jù)影像像元灰度值在這段時(shí)間內(nèi)的變化，對(duì)每個(gè)像元生成的一元線性回歸方程進(jìn)行擬合，通過(guò)分析擬合的一元線性回歸方程的斜率正負(fù)及大小來(lái)分析該像元以及影像在該時(shí)間段內(nèi)的變化以及變化速率。

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

2.2 基于梯度分水嶺的分割算法的礦區(qū)開(kāi)采范圍提取

分水嶺分割算法是一種基于拓?fù)淅碚摰臄?shù)學(xué)形態(tài)學(xué)的分割方法，它把圖像看作是一個(gè)拓?fù)涞孛?，把圖像的每一個(gè)像素的灰度值看作該點(diǎn)的海拔高度，把每一個(gè)局部極小值稱為集水盆地，分水嶺即集水盆地的邊界。

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

(1)

最簡(jiǎn)單的梯度近似表達(dá)式如下：

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

(2)

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

(3)

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

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

(4)

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

(5)

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

3 結(jié)果與分析

3.1 研究區(qū)植被變化

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

圖2 東鞍山植被變化數(shù)據(jù)

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

圖3 大孤山植被變化數(shù)據(jù)

3.2 鐵礦區(qū)開(kāi)采范圍變化

如圖4、圖5所示，分別為東鞍山礦區(qū)和大孤山礦區(qū)的梯度分割結(jié)果與Google Earth目視解譯結(jié)果的對(duì)比。參考2019年Google Earth影像目視解譯結(jié)果對(duì)梯度分割算法提取的礦區(qū)邊界進(jìn)行精度評(píng)定，利用梯度分割算法自動(dòng)提取的東鞍山礦區(qū)提取精度達(dá)到85.60%，大孤山礦區(qū)提取精度達(dá)到86.99%。

圖4 東鞍山目視解譯及其分割結(jié)果

圖5 大孤山目視解譯及其分割結(jié)果

如圖6所示， 東鞍山礦區(qū)面積整體呈現(xiàn)減小趨勢(shì)。1999～2007年間，東鞍山礦區(qū)開(kāi)采破壞面積不斷減小，在這幾年間，由于鋼鐵工業(yè)的快速發(fā)展和政府政策的導(dǎo)向造成了鐵礦石需求量增大，鐵礦石價(jià)格升高，中國(guó)鐵礦資源的開(kāi)采在這幾年明顯加快了速度，鐵礦石產(chǎn)量和年消耗成品礦數(shù)量大幅增加，以滿足鋼鐵工業(yè)增長(zhǎng)的需求。而在2007～2008年，這一年間的開(kāi)采面積從1.104 3 km2增加到1.443 6 km2，原因是在該年鞍鋼集團(tuán)擴(kuò)大露天礦開(kāi)采面積，在東鞍山處主要擴(kuò)大區(qū)域?yàn)榈V山的東北部。自2008年起，礦區(qū)的面積不斷減小至0.997 4 km2，這是由于此時(shí)東鞍山鐵礦加強(qiáng)“邊開(kāi)采，邊綠化”的措施以及鞍山市發(fā)布青山工程規(guī)劃[15]。自2008年起，政府對(duì)東鞍山附近露天采礦場(chǎng)以及多個(gè)排土場(chǎng)進(jìn)行了復(fù)墾、綠化，因此礦區(qū)面積在近幾年間不斷減小。

圖6 東鞍山礦區(qū)變化監(jiān)測(cè)

如圖7所示，大孤山礦區(qū)面積整體呈現(xiàn)出比東鞍山礦區(qū)面積更大的減小幅度。1999～2007年間，大孤山礦山面積整體呈減小趨勢(shì)。2002年相比2001年礦區(qū)面積有大幅度的增長(zhǎng)，這是因?yàn)樵?002年，大孤山西南部排土場(chǎng)的增加。2007～2009年間，為了恢復(fù)大孤山生產(chǎn)能力，延長(zhǎng)礦山壽命，大孤山鐵礦北部擴(kuò)大工程正式開(kāi)工，這一擴(kuò)建工程使這幾年間的礦區(qū)面積增加了1.176 3 km2。在2009年后，政府實(shí)施“邊開(kāi)采，邊復(fù)墾”的政策，恢復(fù)大孤山地區(qū)綠化程度以及礦區(qū)生態(tài)環(huán)境，礦區(qū)開(kāi)采面積減少了2.883 6 km2。

圖7 大孤山礦區(qū)變化監(jiān)測(cè)

4 結(jié)論與討論

本文通過(guò)利用Google Earth Engine(GEE)云平臺(tái)，使用Landsat系列數(shù)據(jù)對(duì)位于鞍山市鐵礦群的東鞍山鐵礦區(qū)和大孤山鐵礦區(qū)進(jìn)行了近20年的植被覆蓋變化監(jiān)測(cè)以及礦區(qū)開(kāi)采面積變化監(jiān)測(cè)，通過(guò)分析數(shù)據(jù)及結(jié)果得出以下結(jié)論：

(1)基于NDVI時(shí)間序列能夠從時(shí)空反映露天開(kāi)采植被的變化，1999～2019年間，東鞍山、大孤山礦區(qū)及影響范圍內(nèi)NDVI指數(shù)整體呈增加趨勢(shì)，尤其礦區(qū)排土場(chǎng)由于進(jìn)行土地復(fù)墾和植被恢復(fù)，NDVI值呈現(xiàn)顯著增加，說(shuō)明礦區(qū)復(fù)墾工作初見(jiàn)成效。

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

(3)利用GEE云平臺(tái)強(qiáng)大的數(shù)據(jù)庫(kù)及其數(shù)據(jù)處理功能，實(shí)現(xiàn)了對(duì)鞍山市東鞍山鐵礦區(qū)和大孤山鐵礦區(qū)的植被變化以及礦區(qū)開(kāi)采面積變化的動(dòng)態(tài)監(jiān)測(cè)，突破了本地硬件條件的限制，實(shí)現(xiàn)在線數(shù)據(jù)處理及分析，且大量節(jié)省了查找數(shù)據(jù)與下載數(shù)據(jù)的時(shí)間。