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

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

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

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

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

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

1.1 研究區(qū)概況

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

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

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

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

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

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

2 研究方法

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

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

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

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

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

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

(1)

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

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個相鄰像素之間的差值。

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

(4)

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

(5)

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論與討論

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

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

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

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