基于Landsat時序數據的高潛水位煤礦區植被擾動分析

方良成，陳永春，安士凱，徐燕飛，殷夢杰，李志輝，趙 萍，陳業禹

(1.淮南礦業(集團)有限責任公司，安徽 淮南 232001；2.平安煤炭開采工程技術研究院有限責任公司，安徽 淮南 232001；3.合肥工業大學 資源與環境工程學院 空間信息智能分析與應用研究所，安徽 合肥 230009；4.安徽大學 資源與環境工程學院 安徽省礦山生態修復工程實驗室，安徽 合肥 230601)

我國是全球最大的煤炭開采和消費國[1]，煤炭在我國的一次性能源結構中長期保持主體能源地位。然而，煤炭資源的開發利用也不可避免地造成了地表塌陷、大氣污染、水污染、土地侵蝕、農作物減產等一系列生態環境問題，阻礙礦區社會經濟的可持續發展[2-4]。

植被作為生態環境的重要組成部分，其長勢是衡量礦區生態環境的重要指標之一[5-8]。近年來，Landsat系列遙感數據的免費開放為植被長時序變化分析提供了數據源。歸一化差值植被指數NDVI(Normalized Difference Vegetation Index)是近紅外波段與紅光波段反射率的差與和之比，可以靈敏地對植被進行檢測，反映地表植被生長狀況、植被覆蓋度、生物量及環境要素影響[6-9]。眾多學者基于NDVI對煤礦區生態環境監測做了大量研究，蔡宗磊等[10]基于GF-1與SPOT6數據，以NDVI等多個植被指數為自變量，對北方露天煤礦區植被覆蓋度進行估測，發現基于NDVI的模型精度最高；黃海[11]利用MODIS和Landsat遙感數據反演礦區NDVI等多個植被指數，分析礦區地表植被和土壤濕度的空間時序變化特征和演化機理，有效地揭示西部黃土礦區煤炭開采對地表植被與土壤濕度的影響；黃翌等[12]基于像元二分法，以NDVI值為參數計算半干旱煤礦區植被覆蓋度，探討礦區植被覆蓋度的空間相關性和空間異質性的變化；李晶等[13]結合NDVI和綜合森林(IFZ)重建阿巴拉契亞地區1987年以來27 a間LUCC動態變化特征，探究采礦擾動隨開采規模變化的關系。上述研究對象多為露天煤礦或干旱半干旱礦區，針對我國東部水熱條件較好的高潛水位煤礦區研究尚不多見[14]。筆者以安徽淮南顧橋礦采煤沉陷區為研究對象，利用2007—2018年Landsat影像構建NDVI時間序列數據，通過NDVI年際變化分析、熱點分析、聚類與異常值分析、剖面線分析，研究煤炭開采沉陷對高潛水位平原礦區植被的擾動效應，以期為礦區生態環境評價提供參考。

1 研究區概況

淮南礦區位于安徽省中北部的淮南市，地理坐標為東經116o20'～117o14'，北緯32o22'～33o22'，處于暖溫帶和亞熱帶的過渡地區[15]。礦區地理位置優越，氣候適宜，煤炭儲量豐富，植被以農作物為主，成片林地與草地較少，是我國華東地區重要的煤糧生產基地。由于地下水位埋藏淺，可采煤層厚度大、層數多、埋深大，重復開采擾動多，淮南礦區煤礦地表沉陷具有規模大、積水深、影響時間長等特點，導致了突出的生態環境問題。顧橋礦位于淮南礦區中西部，正式投產于2007年4月28日，是亞洲地下開采規模最大的礦井之一，煤炭開采形成的采煤沉陷區具有典型性，因此，選取顧橋礦采煤沉陷區為研究對象，探討煤炭開采對周圍植被造成的擾動效應，研究區域地理位置及范圍如圖1所示。

2 數據獲取與處理

2.1 Landsat時間序列NDVI數據的構建

圖1 研究區地理位置和范圍Fig.1 Geographical location and scope of the research area

本文基于2007—2018年Landsat遙感影像進行長時間序列NDVI指數的構建。為盡可能保證時相一致，選用的12期影像數據采集時間均為每年的4、5月份，影像分辨率為30 m，質量較好，基本無云覆蓋，具體情況見表1。

表1 Landsat 系列數據基本信息Table 1 Basic information of Landsat series data

在對遙感影像進行輻射定標、大氣校正、交叉定標[16]、裁剪等圖像預處理后，根據如下NDVI計算公式，計算得到每一期影像的NDVI值。

式中：NIR、R分別代表近紅外、紅光波段的像元反射率值，分別對應TM或ETM+影像4、3波段，或OLI影像5、4波段。NDVI值的范圍為[–1,1]，反映了植被的生長狀態；負值與云、雪、水等相關，0一般代表裸土或巖石，在有植被生長的區域，正值越接近1表示植被覆蓋度越高[17]。

2.2 沉陷積水區提取

水體對太陽光具有強吸收性，在近紅外、短波紅外的波長范圍幾乎可吸收全部的入射能量。可以利用水體在可見光波段和近紅外或短波紅外波段的反差來構建水體指數[18]。MNDWI(Modified Normalized Difference Water Index)除能抑制植被的背景噪聲外，還可有效抑制土壤及建筑物對水體提取的干擾[19]。MNDWI自提出至今已在國內外得到廣泛應用與認可，因此，本次研究采用MNDWI水體指數來進行沉陷積水區的提取，MNDWI的計算公式為：

式中：GREEN、MIR分別代表綠光、短波紅外波段的像元反射率值，分別對應TM或ETM+影像2、5波段，或OLI影像3、6波段。

3 NDVI時序變化特征分析

將各年份NDVI影像按如下標準(表2)進行分級，結果如圖2所示。

表2 NDVI等級劃分標準Table 2 NDVI classification standard

對各年份不同等級NDVI進行統計，如圖3所示，可知研究區NDVI值主要分布在0.4～1，所占百分比均在70%以上。隨著煤炭的開采，地表發生沉陷積水，NDVI值小于0.1的區域所占百分比逐年增加。研究區內各年份沉陷積水區的擴張(圖2)與開采工作面一致，表現為南北向延伸和東向擴展。西北部雖也有少量開采工作面，但由于地表為生產和城鎮生活區，開采過程中根據相關規范留設大量保護煤柱，因此，沉陷區未進一步向西擴展[20]。

利用MNDWI提取的沉陷積水區范圍制作掩膜，統計掩膜后研究區年際NDVI均值與變異系數，如圖4所示。年際NDVI均值在不同階段有所波動，變異系數總體呈上升趨勢。變異系數最小的時間為2007年與2008年，表明煤炭開采初期，研究區內植被的總體生長狀況未發生明顯改變。隨著煤炭的開采，整體的NDVI離散程度逐漸增大。年均值最小的也是2007年和2008年，由文獻[5]可知，礦區存在大量耕地，植被NDVI值變化主要受研究區內農作物耕種和收割活動影響，5—6月為冬小麥成熟及收獲季節，NDVI呈下降趨勢。結合當地歷史氣溫數據及降水數據發現，2007、2008年由于水熱條件好，農作物成熟期提前，農作物大多被收割，導致其NDVI與其他年份相比較低。尤其是2007年，年平均氣溫為17.4℃，高出常年1.8℃，全年除7月份平均氣溫低于常年0.5℃外，其余11個月份平均氣溫均比往年偏高且降雨量豐富，雨水偏多60%以上[21]。

圖2 2007—2018年NDVI分級Fig.2 NDVI classification from 2007 to 2018

圖3 2007—2018年NDVI占比Fig.3 The proportion of NDVI from 2007 to 2018

圖4 研究區2007—2018年NDVI均值與變異系數變化Fig.4 The variety of mean value and variation coefficient from 2007 to 2018

4 NDVI空間分異特征

受以煤炭開發為主的多種因素共同作用，礦區植被覆蓋度表現出明顯的局部依賴性和異質性[22]。熱點分析(Getis-Ord Gi*)能夠運用Getis-Ord Gi*統計對數據集中的每個要素進行計算，并在空間上獲取發生聚類的位置，尋找具有顯著統計學意義的熱點，即要素在具有高值特征的同時，鄰近要素也同樣具有高值特征[21]，可以有效揭示植被擾動中的局部效應。聚類與異常值分析(Anselin Loacal Moran I)根據要素位置和要素值對空間的自相關性進行判斷，表達的是空間鄰近位置上數值間的相互趨同或背離的傾向。因此，空間關聯指數Getis-Ord Gi*、Anselin Loacal Moran I可根據空間位置的高低值簇分布情況，有效揭示研究區空間聚簇特征，對礦區植被擾動研究具有重要意義。

4.1 熱點分析

熱點分析(Getis-Ord Gi*)用來識別具有空間統計性的高值聚集與低值聚集，即熱點(Hot Spots)與冷點(Cold Spots)的分布。熱點的分布表示NDVI高值集中出現的位置，冷點的分布表示NDVI低值集中出現的位置，從而，可根據NDVI熱點和冷點區域的空間演變探測人類活動對礦區植被生長的影響。表達式和原理如下：

式中：Gi*為Getis-Ord Gi*統計；wi,j為以距離規則定義的空間權重，同樣空間范圍相鄰為1，不相鄰為0；xj為要素j的屬性值；n為要素總和。Gi*統計出的結果即為Z得分。對于具有顯著統計學意義的正值Z得分，Z得分表示標準差的倍數，Z得分越高，熱點的聚類就越緊密，即NDVI高值聚集的情況越明顯。反之，負值的Z得分越高，冷點的聚類就越緊密，即NDVI低值聚集的情況越明顯。按照顯著性水平分別為“0.01水平(雙側)上顯著相關”、“0.05水平(雙側)上顯著相關”、“0.1水平(雙側)上顯著相關”與“不相關”，將結果顯示為冷點區(Cold Spot～99% Confidence)、較冷點區(Cold Spot～95%Confidence)、次冷點區(Cold Spot～90% Confidence)、無顯著性(Not Significant)、熱點區(Hot Spot～99%Confidence)、較熱點區(Hot Spot～95% Confidence)、次熱點區(Hot Spot～90% Confidence)。

從圖5可知，2007—2018年間熱點聚集和冷點聚集區域總體上均在增大。具體來看，2007—2008年，農作物為成熟季節，大部分農作物被收割，熱點聚集和冷點聚集的區域較小；2009—2012年熱點聚集和冷點聚集的區域明顯增大；2013—2014年間熱點區和較熱點區的數量明顯減少，次熱點區和冷點區的數量有所增加。2015年以后，熱點和冷點聚集區域又開始增加，且逐漸趨于穩定。

從空間分布上看，隨著煤炭的開采，積水沉陷面積逐年變大，冷點圍繞沉陷區逐年增加。同時，在永幸河與德上高速沿線冷點區增加的現象也較為明顯。德上高速與永幸河是礦區煤炭運輸的主要道路和居民生活的軸線，植被覆蓋度變化十分劇烈，人類活動對地表植被擾動程度顯著，大面積的建筑修建使該處的植被覆蓋度急劇降低。2007—2008年間熱點區分布相對較少，主要分布于研究區北部，2009年后熱點分布增多，幾乎覆蓋整個研究區；2013—2014年間熱點分布區減少，大多轉化為較熱點和次熱點；2015年后熱點區又開始在北部聚集。

4.2 聚類與異常值分析

聚類與異常值分析(Anselin Loacal Moran I)能有效識別具有統計顯著性的空間異常值，發現研究區內NDVI高值與低值聚集的邊界以及高低值分布的異常模式，揭示植被的干擾效應。其工作原理如下：

式中：Ii為Local Moran’s I指數；xi為要素i的屬性；為對應屬性的平均值。

Anselin Loacal Moran I可區分具有顯著統計性0.05水平(雙側)上顯著相關的“高–高”聚類(High-High Cluster)、“低值被高值包圍”類(Low-High Outlier)、無顯著性類(Not Significant)、“高值被低值包圍”類(High-Low Outlier)及“低–低”聚類(Low-Low Cluster)。

圖5 2007—2018年NDVI熱點分析Fig.5 Hot spots analysis of NDVI from 2007 to 2018

分析結果如圖6所示，2007—2010年間，顧橋礦附近的“高–高”聚類逐漸消失，“低–低”聚類現象有明顯增加。2011—2012年“高–高”聚類又再次集中出現；2013—2016年“高–高”聚類效應明顯減弱，基本只有“低–低”聚類現象出現且集中出現在沉陷積水區、崗河、永幸河和德上高速附近；2017年以后，“高–高”聚類的區域面積又明顯增大。在整個研究時間段內，研究區出現的空間聚集情況主要為“高–高”聚類及“低–低”聚類，無異常出現，表明礦區植被局部空間差異程度小，且聚集現象的增加和減少都是整體性的。

4.3 剖面線分析

顧橋礦正式投產于2007年4月28日，開采方式為地下開采，2008年地表開始出現沉陷。為排除其他人類活動干擾，分析煤礦開采沉陷對周邊植被長勢的影響，以2008年沉陷積水區幾何中心為起點,選擇沉陷積水區正北方向剖面線(圖1中AB線)對NDVI進行分析。該方向上地表覆蓋類型以耕地為主，且變化較小。根據NDVI變化趨勢，分別分析煤炭開采不同階段對沉陷積水區周邊植被的影響，將NDVI增減趨勢與上一年相比無明顯變化的點定義為NDVI穩定值(圖7和表3)，圖中的虛線表示沉陷起點到影響邊界(NDVI值穩定的位置)的距離。

圖6 2007—2018年NDVI聚類與異常值分析Fig.6 Clustering and outliner analysis of NDVI from 2007 to 2018

圖7 2007—2018年NDVI變化剖面線分析Fig.7 NDVI variation profile analysis from 2007 to 2018

由圖7和表3可知，2007—2009年，采煤造成的沉陷積水區顯著增大，NDVI值急劇下降的范圍由沉陷積水邊緣外的360 m擴大到420 m。2010—2013年，沉陷積水區緩慢增大，NDVI受影響的范圍由150 m緩慢增大至210 m，且NDVI逐年下降。2013年，因近岸水生植物生長，810 m范圍內NDVI呈上升趨勢，之后較2012年仍呈下降趨勢。2014—2018年，沉陷積水區幾乎沒有擴張，NDVI影響范圍保持為510 m，但在720～1 230 m距離區間內，NDVI仍為逐年下降趨勢。其中，2014年，1 650～1 800 m出現了較大的NDVI谷值，經比對高分辨率衛星影像，此處為蔬菜大棚。2017年，650～800 m的NDVI值較2016年也因近岸水生植物的生長出現增長；2018年，800～900 m 范圍內NDVI發生急劇下降，經比對該處為魚塘。

表3 2007—2018年沉陷積水影響分析數據Table 3 The impact of subsidence water from 2007 to 2018

綜合上述分析可知，煤炭開采當年對植被擾動影響較小，開采一年后，隨著沉陷積水的出現和面積的增長，擾動影響范圍不斷增加；當沉陷積水區穩定后，其影響范圍仍隨時間推移持續增加，最后趨于穩定，說明采煤沉陷對植被的擾動存在時序滯后性和時空累積性。

5 結論

a.通過分析不同等級NDVI所占比例及礦區各年份NDVI的變異系數發現，NDVI值大于0.4的植被在研究區占主導地位，所占百分比均在70%以上，但NDVI值的變異系數總體呈上升趨勢，表明植被生長狀況整體良好，但植被覆蓋的離散程度在增大。

b.NDVI熱點分析、聚類和異常值分析結果顯示，2007—2018年間，研究區內NDVI具有明顯的聚集特征，局部空間差異程度小，均為“高–高”聚類或“低–低”聚類，且聚集現象的增加和減少都是整體性的。總體表現為熱點在減少，冷點在增加，二者間的轉化主要發生在沉陷積水區、德上高速和永幸河附近，主要由采煤沉陷、工礦和居民區建設引起。

c.受煤炭開采影響，沉陷積水區周圍一定范圍內存在明顯的植被擾動效應，隨著沉陷積水區范圍的增長和時間的推移，影響范圍逐漸增大，最后趨于穩定，具有時序滯后性和時空累積性特征。