李 軍 張藝藂 張彩月 謝慧真 張成業(yè) 杜夢豪 王雅穎

(1.中國礦業(yè)大學(北京)地球科學與測繪工程學院,北京 100083;2.煤炭資源與安全開采國家重點實驗室,北京 100083)

煤炭基地長時間的采礦活動對土地造成的挖損、壓占、用地類型轉變等,給生態(tài)環(huán)境帶來了不容忽視的負面影響[1-3],其中最直接的影響之一就是對地表原有植被造成損毀。遙感具有大面積快速觀測的優(yōu)勢,利用長時序遙感技術準確識別煤炭基地植被損毀區(qū)域和損毀時間,對當地生態(tài)環(huán)境質量監(jiān)管具有重要的現實意義[4-5]。

為實現植被損毀遙感識別,國內外學者提出了多種算法。主要包括Sen+Mann-Kendall算法、趨勢分析法、LandTrendr (LT)算法、Continuous Change Detection and Classification (CCDC)算法等。鐘琪等[6]、王科雯等[7]分別利用Sen +Mann-Kendall算法擬合了大寧礦區(qū)和神東礦區(qū)植被的損毀趨勢,并對變化速率進行了分級。徐佳等[8]、劉英等[9]利用趨勢分析法,通過最小二乘原理進行逐像元的一元線性擬合,計算像元回歸直線的斜率,實現了對神東礦區(qū)植被損毀的趨勢識別,并分別預測了礦區(qū)植被未來變化趨勢,分析了影響植被的驅動因子。但Sen+Mann-Kendall算法和趨勢分析法只體現了長時序植被損毀的變化趨勢,無法準確判別植被損毀的位置與時間。與之相比,LT算法和CCDC算法能夠獲取像元尺度的植被變化軌跡,進而確定植被損毀發(fā)生的位置、時間等信息,因此被廣泛應用于植被損毀識別中。在大尺度森林背景下,殷崎棟等[10]、劉珊珊等[11]利用LT算法分別實現了陜西省和浙江省森林區(qū)域的植被損毀識別。肖真等[12]、李賀等[13]利用CCDC算法,分別實現了江西省森林區(qū)域和緬甸南部橡膠林區(qū)域的植被損毀識別。在小尺度礦區(qū)背景下,YANG等[14]、楊張瑜等[15]利用LT算法分別對澳大利亞的庫拉格煤礦以及中國的紫金山金銅礦進行了植被損毀的動態(tài)識別;ZHANG等[16]利用CCDC算法實現了對江西德興銅礦的植被損毀動態(tài)識別。但需要注意的是,即使在同一區(qū)域,不同算法進行植被損毀識別的效果不盡相同[17-18]。因此,需要對LT算法和CCDC算法在同一區(qū)域的植被損毀識別結果進行進一步的定量分析與對比,確定其適用性。

近年來,Google Earth Engine (GEE)憑借其海量的數據資源和強大的超算能力成為了快速遙感云計算的利器[19]。利用GEE云計算的優(yōu)勢,LT和CCDC算法計算速度得到了極大的提升,為區(qū)域植被損毀識別提供了大規(guī)模運用的機會。但是,以往礦業(yè)領域的研究主要是針對單個礦區(qū),缺少對LT和CCDC算法以煤炭基地尺度為場景進行植被損毀的研究與應用,未曾定量對比分析這兩種算法的適用性,導致相關工作者在選擇煤炭基地尺度的植被損毀識別方法時主觀性較強,缺少科學的數據支撐。

為了解決上述問題,本研究基于GEE云計算平臺,以Landsat長時序衛(wèi)星影像數據為基礎,分別采用LT和CCDC算法對神東煤炭基地進行植被損毀識別,獲取了1990—2020年植被損毀識別結果,從植被損毀識別、不同地表區(qū)域適用性、植被損毀時間3個角度,對比分析了兩種算法識別神東煤炭基地植被損毀的適用性,為兩種算法在煤炭基地尺度的進一步應用提供了科學數據參考。

1 研究區(qū)與數據

1.1 研究區(qū)概況

神東煤炭基地位于晉陜蒙三省交界處,地理范圍為東經109°41′～111°36′,北緯38°42′～40°06′,如圖1所示。自1985年來,神東煤炭開采行業(yè)逐漸發(fā)展,“十二五”期間,隨著國家的重點建設,形成了千萬噸礦井群的生產格局,是目前我國最大的煤炭生產基地和重要的優(yōu)質動力煤出口基地,承擔了國家能源保供的重大責任[20]。然而,該地區(qū)位于干旱半干旱地帶,屬于大陸性季風氣候,年均降雨量僅有400 mm,水資源稀少[21]。同時,研究區(qū)內晉陜蒙各地城鎮(zhèn)建設發(fā)展迅速,例如研究區(qū)西北部的鄂爾多斯市東勝區(qū),自1990年以來,建成區(qū)面積呈明顯擴張趨勢,1990年城市建成區(qū)面積為15.11 km2,2008年達到57.79 km2,2020年達到78 km2,其中2000—2008年城市擴張速率最高[22]。在自然本底條件和人類城鎮(zhèn)建設、煤炭工業(yè)活動的共同作用下,神東煤炭基地生態(tài)環(huán)境復雜多變、較為脆弱[23],履行國家能源保供任務的責任與脆弱的生態(tài)環(huán)境保護之間的矛盾突出。

1.2 數據源

本研究采用了1990—2020年的Landsat系列衛(wèi)星地表反射率數據集 (SurfaceReflectance,SR),分別為Landsat 5的Thematic Mapper (TM)、Landsat 7的Enhanced Thematic Mapper Plus (ETM+)、Landsat 8的Operational Land Imager (OLI),其中TM影像2 023景、ETM+影像2 047景、OLI影像943景。所有數據均經過輻射定標、大氣校正、去云、拼接、裁剪等處理。

2 技術方法

2.1 技術路線

本研究通過LT和CCDC算法對研究區(qū)植被損毀進行長時間序列的監(jiān)測,提取植被損毀信息。利用Google Earth的高分辨率遙感影像,通過目視解譯對神東煤炭基地的植被損毀進行精度驗證。從植被損毀識別評價、不同地表區(qū)域適用性評價、植被損毀時間提取精度評價3個方面,對兩種算法在神東煤炭基地進行植被損毀識別的適用性進行對比分析。技術路線如圖2所示。

2.2 LandTrendr算法

由KENNEDY等[24]提出的LT算法是一種時間分割算法,該算法通過對以年為間隔的時間序列數據進行分割、擬合和平滑,生成像元尺度上數據隨時間變化的軌跡,并通過軌跡獲取植被損毀發(fā)生的時間以及植被損毀程度等信息[25]。本研究軌跡構建步驟為:① 去除云或云陰影所引起的異常頂點;② 定義軌跡構造的最大段數;③ 通過去除最小角度的段對初始的時序軌跡進行簡化;④ 評估簡化的時序軌跡與原始軌跡之間的差異。因此,LT算法需要進行相關參數設置以確保識別質量,本研究LT算法的參數取值見表1。

表1 LT算法參數設置Table 1 Parameters setting of LT algorithm

2.3 CCDC算法

由ZHU等[26]提出的CCDC算法是一種基于Landsat衛(wèi)星觀測數據的像元尺度時間序列分割算法。有別于其他算法,CCDC算法會利用時間序列內所有可利用的觀測數據,即通過Fmask算法對云和云陰影進行掩膜處理,當獲得新的Landsat衛(wèi)星觀測數據時,CCDC算法通過結合每個像元所有可用的Landsat衛(wèi)星觀測結果來擬合時間序列模型。如果新的連續(xù)觀測值超過預期范圍(3倍均方根誤差),則進行斷點標記并估計新的時間序列模型,即在斷點前后生成兩個時間序列段。CCDC算法可以識別出規(guī)定時間序列內的所有斷點,并獲取植被損毀時間等信息。本研究CCDC算法參數設置見表2。

表2 CCDC算法參數設置Table 2 Parameters setting of CCDC algorithm

2.4 基于GEE平臺的植被損毀識別

GEE云計算平臺集成了海量地理空間數據,具有強大的分析計算和空間可視化能力,是目前最先進的地理大數據分析平臺之一[27]。在GEE云計算平臺中,利用JavaScript編程語言,能夠快速利用LT和CCDC算法實現區(qū)域植被損毀識別[28-30]。歸一化差值植被指數(Normalized Difference Vegetation Index,NDVI)可描述植被的生長狀態(tài),是評定區(qū)域植被綜合狀況的常用指標之一[31]。因此,本研究以Landsat衛(wèi)星影像為基礎,通過NDVI指數表征植被生長狀態(tài),基于GEE云計算平臺分別使用LT算法和CCDC算法對神東煤炭基地植被進行損毀識別。其中NDVI通過對衛(wèi)星影像進行如下波段運算得到:

式中,NIR為近紅外波段地表反射率;R為紅光波段地表反射率。

另外,需要注意的是:① 部分像元不止一次發(fā)生植被損毀事件,為便于后續(xù)統(tǒng)一進行分析處理,本研究采用像元第一次發(fā)生植被損毀的時間;② LT算法會將研究的起始年份(本研究為1990年)識別為植被損毀的開始點,造成錯分現象。因此在實際應用中需要將起始年份(1990年)識別到的植被損毀結果進行剔除。

2.5 適用性分析

本研究從植被損毀識別精度、不同地表區(qū)域的適用性、植被損毀時間提取精度3個方面,以總體精度、生產者精度、用戶精度、錯分誤差、遺漏誤差等作為評價指標,對比分析了LT和CCDC算法在神東煤炭基地植被損毀識別中的適用性。其中,總體精度是被正確分類的像元總和除以總像元數,即混淆矩陣中對角線上的各類別像元數總和除以影像總像元數目。生產者精度是某一類別被正確分類的像元數(混淆矩陣對角線上該類別的數量)除以該類別的地面真實像元總數。用戶精度是某一類別被正確分類的像元數(混淆矩陣對角線上該類別的數量)除以被算法識別為該類別的像元總數。錯分誤差是在未發(fā)生植被損毀地區(qū)的樣本區(qū)域中,被錯誤識別為損毀區(qū)域的面積與樣本區(qū)域總面積的比值。遺漏誤差是在發(fā)生植被損毀地區(qū)的樣本區(qū)域中,未被識別為損毀區(qū)域的面積與樣本區(qū)域總面積的比值。

(1)植被損毀識別精度評價。本研究針對兩種算法識別結果,分別在識別到植被損毀的區(qū)域和未識別到植被損毀的區(qū)域隨機選擇了250個像元作為樣本。對于每個樣本像元,通過Google Earth提供的歷史高分辨率影像進行目視判讀,判斷該像元是否發(fā)生了植被損毀,作為精度評價的真值,生成混淆矩陣,計算總體精度、生產者精度和用戶精度。

(2)不同地表區(qū)域的適用性評價。本研究以Google Earth中神東煤炭基地原始高分辨率影像為基礎,通過目視解譯勾畫了水體、林地、城市擴張用地、露天采場4類典型樣本區(qū)域(圖1)。水體和林地樣本區(qū)域未發(fā)生植被損毀事件,在此區(qū)域識別的植被損毀為算法造成的錯分所致。城市和露天采場區(qū)域地表植被完全破壞,在此區(qū)域未識別到的植被損毀為算法造成的遺漏所致。通過樣本區(qū)域內兩種算法識別的植被損毀面積來判斷算法的錯分誤差和遺漏誤差。

(3)植被損毀時間精度評價。本研究面向算法識別出的植被損毀事件,每年隨機選取50個樣本點(共1 500個),通過Google Earth中的歷史高分辨率影像進行目視解譯,確定植被損毀的實際時間,作為精度評價的真值,生成混淆矩陣,計算總體精度、生產者精度和用戶精度。

3 結果與分析

3.1 植被損毀識別結果

基于GEE云計算平臺,使用LT算法和CCDC算法對神東煤炭基地進行了植被損毀識別,結果分別如圖3、圖4所示。圖中A、B區(qū)域分別為圈劃的鄂爾多斯市東勝區(qū)的城市建成區(qū)和研究區(qū)露天采場集中分布區(qū)域。

圖3 LT算法植被損毀識別結果Fig.3 Identification results of vegetation damage of LT algorithm

圖4 CCDC算法植被損毀識別結果Fig.4 Identification results of vegetation damage of CCDC algorithm

對比圖3和圖4可知:① 總體而言,LT和CCDC算法在研究區(qū)西部識別到的植被損毀多于東部地區(qū)。LT算法能夠識別到的植被損毀明顯少于CCDC算法。② 在鄂爾多斯市東勝區(qū)的城市建成區(qū)(A區(qū)域)中,LT算法基本無法識別植被損毀,而CCDC算法的植被損毀識別效果較好,且算法所識別的植被損毀多發(fā)生于2000—2008年,這一現象與東勝區(qū)城市擴張歷史相吻合。③ 在B區(qū)域中,LT算法能夠識別到的植被損毀比A區(qū)域顯著增多,但仍有部分露天采場區(qū)域的植被損毀無法識別(圖3(e)中圈定區(qū)域)。相較之下,CCDC算法相同區(qū)域(圖4(e)中圈定區(qū)域)識別到的植被損毀能夠更好地吻合露天采場邊界。

3.2 植被損毀識別精度評價

利用選擇的樣本像元,對兩種算法的植被損毀識別結果進行了驗證,精度分別如表3和表4所示。

表3 LT算法識別結果的混淆矩陣Table 3 Confusion matrix of the identifying results with LT algorithm

由表3、表4可知:LT算法的總體精度為73.6%,CCDC算法的總體精度為84.4%,LT算法比CCDC算法低10.8個百分點。LT算法“損毀像元”的生產者精度和“非損毀像元”識別的用戶精度均比較低,分別為66.9%和54%。CCDC算法“非損毀像元”的生產者精度以及“損毀像元”的生產者精度、用戶精度都達到了80%以上,“非損毀像元”的用戶精度為78.8%,也顯著高于LT算法對應的54.0%。通過與已有研究結果進行對比發(fā)現,本研究所得總體精度略低于以往的研究結果(LT算法74.4%～93%,CCDC算法85.1%～99.4%)[10-16]。

3.3 不同地表區(qū)域適用性評價

兩種算法對于研究區(qū)植被損毀識別的遺漏誤差、錯分誤差分別見表5、表6。

表5 植被損毀遺漏誤差Table 5 Omission errors of vegetation damage

表6 植被損毀錯分誤差Table 6 Commission errors of vegetation damage

由表5、表6可知:在城市擴張區(qū)域,CCDC算法的遺漏誤差為21.6%,LT算法的遺漏誤差為91.0%,比CCDC算法高69.4個百分點。在露天采場區(qū)域,LT算法的遺漏誤差為63.0%,CCDC算法的遺漏誤差為8.6%,比LT算法低54.4個百分點。由此可見:LT算法僅能識別部分露天采場區(qū)域的植被損毀,對于城市擴張所造成的植被損毀基本無法識別。在林地區(qū)域,LT和CCDC算法的錯分誤差分別為0.7%和11.6%,兩種算法都能有效地避免錯分現象。但在水體區(qū)域,LT和CCDC算法的錯分誤差分別為3.7%和76.0%,CCDC算法難以避免水體區(qū)域的錯分現象。

為分析LT和CCDC算法在不同地表類型區(qū)域造成遺漏的原因,本研究通過GEE云計算平臺提取了兩種算法擬合的像元NDVI時序曲線,典型示例分別如圖5和圖6所示。圖5所示的像元在2007年成為露天采場,NDVI值出現明顯降低,但LT算法擬合的曲線并未在此處進行分段,沒有成功識別出該損毀事件。圖6所示的像元于2000年左右經歷了城市快速擴張,2000年后未再出現圖中虛線框中每年7、8月的周期性高值現象,NDVI數值的時序分布明顯改變,但該變化未被CCDC算法識別出?？傮w而言,在干旱半干旱地區(qū)發(fā)生植被損毀事件后,盡管NDVI值的分布明顯改變,但是變化幅度仍然較小,與植被茂盛地區(qū)(如以往研究中的森林區(qū)域)NDVI的“斷崖式”變化存在較大差異。這導致算法在植被稀疏的干旱半干旱地區(qū)僅能識別出NDVI整體呈降低的變化趨勢,并未實現對斷點變化的識別。受此影響,LT算法僅能識別部分露天采場區(qū)域的植被損毀,對于城市擴張所造成的植被損毀基本無法識別。CCDC算法通過均方根誤差對時間序列模型進行分段,因此該算法在露天采場和城市擴張用地區(qū)域的識別效果較好,相較于LT算法更適用于神東煤炭基地的植被損毀識別。

圖5 LT算法損毀識別的遺漏誤差Fig.5 Omission errors of the damage identification of LT algorithm

圖6 CCDC算法損毀識別的遺漏誤差Fig.6 Omission errors of the damage identification of CCDC algorithm

總體上,兩種算法在神東煤炭基地的識別精度略低于以往研究的原因在于,以往研究均針對森林區(qū)域或小區(qū)域的單個礦區(qū),而神東煤炭基地覆蓋面積大,且位于干旱半干旱地區(qū),植被覆蓋度較低、生態(tài)環(huán)境的地理空間異質性強。一方面,植被覆蓋度低容易導致植被損毀時NDVI時序曲線的下降斷點不明顯,算法識別不出植被損毀現象。另一方面,在煤炭基地尺度的大面積植被損毀自動識別中,LT和CCDC算法自適應性弱,一個算法在整個區(qū)域采用一套相同的參數設置,無法針對不同位置異質性強的特點自適應優(yōu)化調整。總之,兩種算法在神東煤炭基地進行大面積植被損毀識別中效果仍然較好,只是精度稍有降低,根據相關原因分析,兩種算法仍存在改進的空間。

3.4 植被損毀時間精度評價

對于能夠識別到植被損毀的像元,本研究進一步評估了LT和CCDC算法在識別損毀時間上的精度,結果分別如圖7和圖8所示。

圖7 LT算法植被損毀時間的精度驗證結果Fig.7 Precision verification results of vegetation damage time based on LT algorithm

圖8 CCDC算法植被損毀時間的精度驗證結果Fig.8 Precision verification results of vegetation damage time based on CCDC algorithm

兩種算法識別植被損毀時間的總體精度分別為89.2%和79%,LT算法總體精度比CCDC算法提高了10.2個百分點。以“十二五”規(guī)劃的起始時間(2011年)對植被損毀時間序列進行劃分,2011—2020年LT算法總體精度由2011之前的87.1%提高至92%,而CCDC算法總體精度由2011之前的74.1%提高至88.8%,總體精度提升幅度較大。

由圖7和圖8可知:在生產者精度方面,LT算法的生產者精度主要分布于75%～98%,均值為88.7%。CCDC算法的生產者精度變化幅度較大,均值為81.6%。在某些年份中,兩種算法的生產者精度較低,例如1999年,LT算法和CCDC算法的生產者精度分別為68.1%和41.2%,均為歷史最低值。在用戶精度方面,LT算法的用戶精度主要分布于80%～96%,均值為89.4%,用戶精度低值出現在2010年,為78%。CCDC算法的用戶精度均值為79%,用戶精度變化幅度較大,低值出現于1992年和2000年,不足60%。

為分析總體精度提高原因,本研究統(tǒng)計了兩種算法在通過目視解譯不同地表區(qū)域(水體、林地、城市擴張用地、露天采場)典型樣本區(qū)中識別的植被損毀面積與時間信息,如圖9所示。

圖9 1991—2020年不同地表區(qū)域損毀面積統(tǒng)計Fig.9 Statistics of damaged area of different surface areas from 1991 to 2020

由圖9可知:2011年以后,露天采場區(qū)域植被損毀面積顯著提升,這一現象與“十二五”規(guī)劃實施后神東煤炭基地的蓬勃發(fā)展有著密切聯系。LT和CCDC算法在露天采場典型樣本區(qū)域的植被損毀識別遺漏誤差優(yōu)于城市擴張典型樣本區(qū)域,因此露天采場面積的增加帶動了總體精度得到不同程度的提高。

本研究進一步統(tǒng)計了CCDC算法和LT算法識別損毀時間的誤差分布,以“誤差=實際損毀時間-算法識別植被損毀時間”方法計算了誤差并取絕對值,統(tǒng)計結果如圖10所示。根據統(tǒng)計信息,在LT算法中,準確識別損毀時間的比例為88.7%,誤差絕對值為1 a的比例為7%;對于CCDC算法,準確識別損毀時間的比例為79.1%,誤差絕對值為1 a的比例為10.9%。換言之,LT和CCDC算法識別植被損毀時間的誤差≤1 a的比例分別為95.7%和90%。上述結果說明,對于兩種算法已經識別出的植被損毀區(qū)域,其損毀時間的識別精度均較高。

圖10 LT和CCDC算法誤差分布Fig.10 Error distribution of LT and CCDC algorithms

4 結論與討論

本研究基于GEE云計算平臺,從植被損毀識別精度評價、不同地表區(qū)域的適用性評價、植被損毀時間精度評價等3個方面定量對比分析了LT和CCDC算法在神東煤炭基地進行植被損毀識別的適用性。主要得到以下結論:

(1)在植被損毀識別上,LT和CCDC算法總體精度分別為73.6%和84.4%。從總體精度角度上看,識別效果仍然較好,但是在神東煤炭基地尺度的精度略低于以往研究在其他區(qū)域的精度。

(2)在不同地表區(qū)域適用性分析中,LT和CCDC算法都能較好地避免林地區(qū)域的錯分誤差。LT算法僅能識別露天采場區(qū)域的部分植被損毀,遺漏誤差較大,且基本無法識別到城市擴張所造成的植被損毀,但是CCDC算法在這兩類區(qū)域的識別效果較好。在水體區(qū)域,LT算法顯著優(yōu)于CCDC算法。

(3)對于植被損毀時間,LT和CCDC算法識別的總體精度分別為89.2%和79%,誤差在1 a內的結果分別占95.7%和90%,對于損毀時間識別效果很好。

總體而言,相較于LT算法,CCDC算法更適用于城市擴張明顯、水體面積很少的神東煤炭基地植被損毀識別。本研究成果為神東煤炭基地生態(tài)環(huán)境質量監(jiān)管提供了數據參考,更為兩種算法在其他煤炭基地尺度的進一步應用提供了方法優(yōu)選借鑒。

LT和CCDC算法不依賴于影像分類、植被物候模型等條件,僅使用長時序Landsat影像即可快速、自動化地實現煤炭基地尺度上大范圍長時序的植被損毀識別,獲取植被損毀發(fā)生的時間、位置等信息。但兩種算法都依賴于植被時序軌跡的變化幅度,本身都存在一些局限性,具體來說:① 在干旱半干旱地區(qū),原始植被覆蓋度低,部分像元發(fā)生植被損毀時兩種算法擬合的時序軌跡變化幅度不明顯,導致損毀植被識別的遺漏;② 兩種算法只能識別植被損毀事件,無法對自然因素引起的波動和人類活動導致的突變進行區(qū)分。針對上述問題,在未來研究中,需要面向采礦活動對植被損毀的時序場景特點,研究一種能夠準確自適應識別煤炭基地內大范圍露天礦群植被損毀事件的新方法。