基于Landsat時(shí)序數(shù)據(jù)的高潛水位煤礦區(qū)植被擾動(dòng)分析

方良成，陳永春，安士凱，徐燕飛，殷夢(mèng)杰，李志輝，趙 萍，陳業(yè)禹

(1.淮南礦業(yè)(集團(tuán))有限責(zé)任公司，安徽 淮南 232001；2.平安煤炭開采工程技術(shù)研究院有限責(zé)任公司，安徽 淮南 232001；3.合肥工業(yè)大學(xué) 資源與環(huán)境工程學(xué)院 空間信息智能分析與應(yīng)用研究所，安徽 合肥 230009；4.安徽大學(xué) 資源與環(huán)境工程學(xué)院 安徽省礦山生態(tài)修復(fù)工程實(shí)驗(yàn)室，安徽 合肥 230601)

我國是全球最大的煤炭開采和消費(fèi)國[1]，煤炭在我國的一次性能源結(jié)構(gòu)中長(zhǎng)期保持主體能源地位。然而，煤炭資源的開發(fā)利用也不可避免地造成了地表塌陷、大氣污染、水污染、土地侵蝕、農(nóng)作物減產(chǎn)等一系列生態(tài)環(huán)境問題，阻礙礦區(qū)社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展[2-4]。

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

1 研究區(qū)概況

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

2 數(shù)據(jù)獲取與處理

2.1 Landsat時(shí)間序列NDVI數(shù)據(jù)的構(gòu)建

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

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

表1 Landsat 系列數(shù)據(jù)基本信息Table 1 Basic information of Landsat series data

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

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

2.2 沉陷積水區(qū)提取

水體對(duì)太陽光具有強(qiáng)吸收性，在近紅外、短波紅外的波長(zhǎng)范圍幾乎可吸收全部的入射能量?？梢岳盟w在可見光波段和近紅外或短波紅外波段的反差來構(gòu)建水體指數(shù)[18]。MNDWI(Modified Normalized Difference Water Index)除能抑制植被的背景噪聲外，還可有效抑制土壤及建筑物對(duì)水體提取的干擾[19]。MNDWI自提出至今已在國內(nèi)外得到廣泛應(yīng)用與認(rèn)可，因此，本次研究采用MNDWI水體指數(shù)來進(jìn)行沉陷積水區(qū)的提取，MNDWI的計(jì)算公式為：

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

3 NDVI時(shí)序變化特征分析

將各年份NDVI影像按如下標(biāo)準(zhǔn)(表2)進(jìn)行分級(jí)，結(jié)果如圖2所示。

表2 NDVI等級(jí)劃分標(biāo)準(zhǔn)Table 2 NDVI classification standard

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

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

圖2 2007—2018年NDVI分級(jí)Fig.2 NDVI classification from 2007 to 2018

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

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

4 NDVI空間分異特征

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

4.1 熱點(diǎn)分析

熱點(diǎn)分析(Getis-Ord Gi*)用來識(shí)別具有空間統(tǒng)計(jì)性的高值聚集與低值聚集，即熱點(diǎn)(Hot Spots)與冷點(diǎn)(Cold Spots)的分布。熱點(diǎn)的分布表示NDVI高值集中出現(xiàn)的位置，冷點(diǎn)的分布表示NDVI低值集中出現(xiàn)的位置，從而，可根據(jù)NDVI熱點(diǎn)和冷點(diǎn)區(qū)域的空間演變探測(cè)人類活動(dòng)對(duì)礦區(qū)植被生長(zhǎng)的影響。表達(dá)式和原理如下：

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

從圖5可知，2007—2018年間熱點(diǎn)聚集和冷點(diǎn)聚集區(qū)域總體上均在增大。具體來看，2007—2008年，農(nóng)作物為成熟季節(jié)，大部分農(nóng)作物被收割，熱點(diǎn)聚集和冷點(diǎn)聚集的區(qū)域較?。?009—2012年熱點(diǎn)聚集和冷點(diǎn)聚集的區(qū)域明顯增大；2013—2014年間熱點(diǎn)區(qū)和較熱點(diǎn)區(qū)的數(shù)量明顯減少，次熱點(diǎn)區(qū)和冷點(diǎn)區(qū)的數(shù)量有所增加。2015年以后，熱點(diǎn)和冷點(diǎn)聚集區(qū)域又開始增加，且逐漸趨于穩(wěn)定。

從空間分布上看，隨著煤炭的開采，積水沉陷面積逐年變大，冷點(diǎn)圍繞沉陷區(qū)逐年增加。同時(shí)，在永幸河與德上高速沿線冷點(diǎn)區(qū)增加的現(xiàn)象也較為明顯。德上高速與永幸河是礦區(qū)煤炭運(yùn)輸?shù)闹饕缆泛途用裆畹妮S線，植被覆蓋度變化十分劇烈，人類活動(dòng)對(duì)地表植被擾動(dòng)程度顯著，大面積的建筑修建使該處的植被覆蓋度急劇降低。2007—2008年間熱點(diǎn)區(qū)分布相對(duì)較少，主要分布于研究區(qū)北部，2009年后熱點(diǎn)分布增多，幾乎覆蓋整個(gè)研究區(qū)；2013—2014年間熱點(diǎn)分布區(qū)減少，大多轉(zhuǎn)化為較熱點(diǎn)和次熱點(diǎn)；2015年后熱點(diǎn)區(qū)又開始在北部聚集。

4.2 聚類與異常值分析

聚類與異常值分析(Anselin Loacal Moran I)能有效識(shí)別具有統(tǒng)計(jì)顯著性的空間異常值，發(fā)現(xiàn)研究區(qū)內(nèi)NDVI高值與低值聚集的邊界以及高低值分布的異常模式，揭示植被的干擾效應(yīng)。其工作原理如下：

式中：Ii為L(zhǎng)ocal Moran’s I指數(shù)；xi為要素i的屬性；為對(duì)應(yīng)屬性的平均值。

Anselin Loacal Moran I可區(qū)分具有顯著統(tǒng)計(jì)性0.05水平(雙側(cè))上顯著相關(guān)的“高–高”聚類(High-High Cluster)、“低值被高值包圍”類(Low-High Outlier)、無顯著性類(Not Significant)、“高值被低值包圍”類(High-Low Outlier)及“低–低”聚類(Low-Low Cluster)。

圖5 2007—2018年NDVI熱點(diǎn)分析Fig.5 Hot spots analysis of NDVI from 2007 to 2018

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

4.3 剖面線分析

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

圖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年，采煤造成的沉陷積水區(qū)顯著增大，NDVI值急劇下降的范圍由沉陷積水邊緣外的360 m擴(kuò)大到420 m。2010—2013年，沉陷積水區(qū)緩慢增大，NDVI受影響的范圍由150 m緩慢增大至210 m，且NDVI逐年下降。2013年，因近岸水生植物生長(zhǎng)，810 m范圍內(nèi)NDVI呈上升趨勢(shì)，之后較2012年仍呈下降趨勢(shì)。2014—2018年，沉陷積水區(qū)幾乎沒有擴(kuò)張，NDVI影響范圍保持為510 m，但在720～1 230 m距離區(qū)間內(nèi)，NDVI仍為逐年下降趨勢(shì)。其中，2014年，1 650～1 800 m出現(xiàn)了較大的NDVI谷值，經(jīng)比對(duì)高分辨率衛(wèi)星影像，此處為蔬菜大棚。2017年，650～800 m的NDVI值較2016年也因近岸水生植物的生長(zhǎng)出現(xiàn)增長(zhǎng)；2018年，800～900 m 范圍內(nèi)NDVI發(fā)生急劇下降，經(jīng)比對(duì)該處為魚塘。

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

綜合上述分析可知，煤炭開采當(dāng)年對(duì)植被擾動(dòng)影響較小，開采一年后，隨著沉陷積水的出現(xiàn)和面積的增長(zhǎng)，擾動(dòng)影響范圍不斷增加；當(dāng)沉陷積水區(qū)穩(wěn)定后，其影響范圍仍隨時(shí)間推移持續(xù)增加，最后趨于穩(wěn)定，說明采煤沉陷對(duì)植被的擾動(dòng)存在時(shí)序滯后性和時(shí)空累積性。

5 結(jié)論

a.通過分析不同等級(jí)NDVI所占比例及礦區(qū)各年份NDVI的變異系數(shù)發(fā)現(xiàn)，NDVI值大于0.4的植被在研究區(qū)占主導(dǎo)地位，所占百分比均在70%以上，但NDVI值的變異系數(shù)總體呈上升趨勢(shì)，表明植被生長(zhǎng)狀況整體良好，但植被覆蓋的離散程度在增大。

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

c.受煤炭開采影響，沉陷積水區(qū)周圍一定范圍內(nèi)存在明顯的植被擾動(dòng)效應(yīng)，隨著沉陷積水區(qū)范圍的增長(zhǎng)和時(shí)間的推移，影響范圍逐漸增大，最后趨于穩(wěn)定，具有時(shí)序滯后性和時(shí)空累積性特征。