李瑞軍， 楊志文， 吳士文， 郝仕龍,3， 鄧榮鑫， 尹卿芳

(1.陽泉市大陽泉煤炭有限責任公司，山西 陽泉 045000； 2.華北水利水電大學 測繪與地理信息學院，河南 鄭州 450046； 3.南陽師范學院 南水北調(diào)中線水源區(qū)水安全河南省協(xié)同創(chuàng)新中心，河南 南陽 473061； 4.北京京水建設集團有限公司，北京 100193)

植被是陸地生態(tài)系統(tǒng)的重要組成部分，是連接土壤圈、大氣圈、水圈、生物圈的自然紐帶[1]。植被覆蓋度作為表征地表生態(tài)環(huán)境和植被分布特征的主要參數(shù)，在水土保持、調(diào)節(jié)區(qū)域小氣候、涵養(yǎng)水源等方面發(fā)揮著重要作用，具有明顯的季候和年際變化特征[2-3]。傳統(tǒng)植被覆蓋度測量，多數(shù)是在研究區(qū)選擇有限的幾個或幾十個典型樣點進行實地測量，雖能夠保證數(shù)據(jù)觀測的精度，但費時費力、成本高、難度大，并且難以滿足區(qū)域地表植被覆蓋度的需求[4]；而遙感技術具有快速實時、監(jiān)測范圍廣、人力物力耗費較少、非破壞性等特點，為估算區(qū)域尺度的陸表植被覆蓋度提供了強有力的手段[5]。植被指數(shù)是遙感領域中用來表征植被生長、長勢的有效指標，在眾多植被指數(shù)中，歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)與植被覆蓋度相關性最高[6-8]，已有不少研究將NDVI的像元二分模型應用于植被覆蓋動態(tài)變化監(jiān)測分析[3-4,9]。

煤炭作為我國的主要能源，是我國經(jīng)濟和社會發(fā)展的重要基礎和基本保障[10]。及時、高效地監(jiān)測礦區(qū)生態(tài)環(huán)境狀況，為礦區(qū)資源的合理開發(fā)提供決策參考，對于實現(xiàn)經(jīng)濟發(fā)展和環(huán)境效益的有機統(tǒng)一尤為重要。國內(nèi)外研究發(fā)現(xiàn)，煤炭資源的大規(guī)模、高強度開采雖然促進了當?shù)亟?jīng)濟發(fā)展，也不可避免地造成礦區(qū)內(nèi)地表塌陷、水土流失、植被破壞、地下水系紊亂，進而導致生態(tài)環(huán)境結(jié)構(gòu)、功能被破壞等問題，對地區(qū)生態(tài)環(huán)境和經(jīng)濟社會可持續(xù)發(fā)展造成了嚴重影響[11-15]，在生態(tài)環(huán)境較為脆弱的地區(qū)，該問題表現(xiàn)得尤為突出。另有研究結(jié)果表明，煤炭資源的合理開發(fā)能使礦區(qū)的土地沙漠化逆轉(zhuǎn)、植被覆蓋率提高，整個生態(tài)環(huán)境向良性方向發(fā)展[16-18]。但以上研究均是針對其它地區(qū)，對陽泉煤礦的研究大多集中在煤炭開采、景觀生態(tài)效應和生態(tài)修復等領域[19-21]，在煤礦開采對地表植被覆蓋度的時空動態(tài)變化及其影響因素研究方面還未見報道。本文以陽泉市大陽泉煤礦不同開采年份工作面為研究對象，基于Landsat TM/OLI數(shù)據(jù)并結(jié)合像元二分模型，對研究區(qū)2009年和2019年地表植被覆蓋度的分布格局及時空動態(tài)變化進行分析；將植被覆蓋圖與開采年份、各地形因子疊加統(tǒng)計，探究影響地表植被覆蓋度變化的主要因子，以期為陽泉市煤炭資源的開發(fā)和地表植被生態(tài)管理作出科學、準確的評價。

1 材料與方法

1.1 研究區(qū)概況

大陽泉煤炭有限責任公司為陽泉市地方國有企業(yè)，位于陽泉市西南郊(東經(jīng)113°30′54″～113°33′6″，北緯37°48′54″～37°50′12″)。1984年開始開采，研究區(qū)地理位置及不同開采年份工作面分布如圖1所示。陽泉市位于沁水煤田的東北角，其煤炭資源儲量約1.7×1010t，埋藏深度150～500 m，易于開發(fā)，品質(zhì)優(yōu)良；它是我國最大的無煙煤生產(chǎn)基地，也是山西中東部的政治、經(jīng)濟、文化中心。境內(nèi)以山地為主，地形復雜，大部分地區(qū)海拔700 m以上，地勢西高東低，屬暖溫帶半濕潤大陸性季風氣候區(qū)，年均氣溫8～12 ℃，無霜期130～180 d，年降水量450～550 mm。

圖1 研究區(qū)位置示意圖及不同開采年份工作面分布圖

1.2 數(shù)據(jù)來源與處理

本文采用的遙感數(shù)據(jù)為2009年5月31日的Landsat-5 TM和2019年5月20日的Landsat-8 OLI影像數(shù)據(jù)，下載自中國科學院遙感與數(shù)字地球研究所對地觀測數(shù)據(jù)共享計劃網(wǎng)站(http://ids.ceode.ac.cn)，影像質(zhì)量好，空間分辨率為30 m。運用ENVI 5.1軟件對Landsat影像依次進行輻射定標、大氣校正和幾何校正，將影像的像元亮度值轉(zhuǎn)換為地表真實反射率值。Landsat-5 TM影像為四級產(chǎn)品，基于該影像完成Landsat-8 OLI影像的幾何校正。最后，利用研究區(qū)的邊界矢量數(shù)據(jù)進行裁切，得到研究區(qū)的遙感影像。

DEM數(shù)據(jù)由陽泉市大陽泉煤礦有限責任公司提供，運用ArcGIS 10.2軟件的3D分析模塊，生成空間分辨率為30 m的高程圖、坡度圖及坡向圖，通過重新調(diào)整高程、坡度和坡向分級并進行重分類，再與植被覆蓋度分布圖疊加，探討地形因子及開采年份對研究區(qū)植被覆蓋度分布的影響；利用ArcGIS 10.2軟件，在研究區(qū)隨機生成359個采樣點，提取2009—2019年植被覆蓋度的變化、開采時間、高程等信息，利用地理探測器探測各影響因子的重要性。結(jié)合研究區(qū)實際情況，對各地形因子進行如下等級劃分。

1)高程：研究區(qū)最大高程1 125 m，最低高程766 m，為方便研究，將研究區(qū)高程分為8類，即[766 m,800 m]、(800 m,850 m]、(850 m,900 m]、(900 m,950 m]、(950 m,1 000 m]、(1 000 m,1 050 m]、(1 050 m,1 100 m]、(1 100 m,1 130 m]。

2)坡度：依據(jù)《水土保持綜合治理規(guī)劃通則》(GB/T 15772—1995)，結(jié)合研究區(qū)地形及其對植被長勢的影響，將坡度分級改為(0°,5°]、(5°,8°]、(8°,15°]、(15°,25°]、(25°,35°]、(35°,70°)[22]。

3)坡向：把研究區(qū)分為8個自然坡向，分別為(0°,22.5°]和(337.5°,360°)(北坡)、(22.5°,67.5°](東北坡)、(67.5°,112.5°](東坡)、(112.5°,157.5°](東南坡)、(157.5°,202.5°](南坡)、(202.5°,247.5°](西南坡)、(247.5°,292.5°](西坡)、(292.5°,337.5°](西北坡)。

1.3 研究方法

1.3.1 植被覆蓋度反演

像元二分模型[23]假定單個像元地表S由植被覆蓋地表和裸土地表構(gòu)成。若假定像元地表完全由植被覆蓋的純像元遙感信息為Sveg，完全由裸土覆蓋的純像元遙感信息為Ssoil，則像元植被覆蓋度f為：

f=(S-Ssoil)/(Sveg-Ssoil) 。

(1)

要估算像元植被覆蓋度，只需要計算土壤和植被的純像元遙感信息即可。該模型簡單實用，通過引入?yún)?shù)減弱了土壤背景、大氣以及植被類型等的影響。

歸一化植被指數(shù)(NDVI)是遙感圖像中近紅外波段反射值與紅光波段反射值之差與二者之和的比值，是反映植被生長狀態(tài)和植物空間分布密度的重要參數(shù)。其比值形式也能夠消除大部分與地形、云陰影、太陽角和儀器定標等引起的輻照度變化，以加強對植被的響應。把NDVI代入像元二分模型，其公式為：

fc=(NDVI-NDVIsoil)/ (NDVIveg-NDVIsoil)。

(2)

式中：fc為植被覆蓋度；NDVI為歸一化植被指數(shù)；NDVIsoil為裸土或無植被覆蓋區(qū)域像元的NDVI值；NDVIveg為完全植被覆蓋區(qū)域像元的NDVI值，即純植被像元的NDVI值。

理論上，NDVIsoil的值應接近于0，NDVIveg代表全植被覆蓋像元的最大值。但是，由于受光照條件、影像質(zhì)量、植被類型等眾多因素影響，影像不可避免地存在噪聲，兩者的值均會發(fā)生變化。文中，NDVIsoil取影像中NDVI累計概率在5%附近的值；NDVIveg取影像中NDVI累計概率在95%附近的值[3]，以此可反演得到研究區(qū)兩個時期的植被覆蓋度。

1.3.2 地理探測器

地理探測器是王勁峰等提出的一種基于空間分異理論，揭示其背后驅(qū)動因子的新的統(tǒng)計學方法[24-25]，主要包括因子探測器、交互作用探測器、風險區(qū)探測器等。傳統(tǒng)植被變化驅(qū)動因子分析方法對模型假設條件和數(shù)據(jù)要求較多，當因子過多時，計算過程較為繁瑣[26]；而地理探測器模型假設制約條件較少，已有學者將其用于荒漠化、NDVI、植被覆蓋度等時空演變的驅(qū)動因子研究[27-30]。因子探測器的具體公式為：

(3)

式中：q為度量空間分異性的指標，表示某因子解釋了q×100%的植被覆蓋度；h為影響因子的分層數(shù)，h=1、2、…、L；Nh和N分別為影響因子的層h和全區(qū)的樣本數(shù)；δh和δ分別為層h和研究區(qū)的植被覆蓋度的方差。

q的值域為[0,1]，q值越大，表明該因子對研究區(qū)植被覆蓋度變化的解釋力越強，反之則越弱。

交互探測器用于評估因子共同作用是否會增加或減弱對植被覆蓋度的解釋力。首先，計算兩種因子X1和X2對因變量的q值；其次，計算因子交互的q值，并對q(X1)、q(X2)與q(X1∩X2)進行比較。

2 結(jié)果與分析

2.1 植被覆蓋度空間分布特征

為更好地研究開采時間對地表植被生長的影響情況，選取遠離研究區(qū)且沒開采過的區(qū)域作為對照。因植被覆蓋度為[0,1]，為便于進行植被覆蓋動態(tài)變化分析，結(jié)合植被覆蓋實際情況，將其分為5級：低植被覆蓋度[0.0,0.2]、中低植被覆蓋度(0.2,0.4]、中植被覆蓋度(0.4,0.6]、中高植被覆蓋度(0.6,0.8]、高植被覆蓋度(0.8,1.0]，通過密度分割處理，生成植被覆蓋度分布圖。

研究區(qū)2009年和2019年植被覆蓋度的空間分布如圖2所示，研究區(qū)內(nèi)的空值區(qū)域?qū)儆诜侵脖桓采w區(qū)，已掩膜去除。

圖2 研究區(qū)2009年和2019年植被覆蓋度的空間分布

由圖2可知：2009年研究區(qū)整體覆蓋度較好，多數(shù)區(qū)域達到了中高植被覆蓋度等級，但東北方向區(qū)域(2000年以前開采)的植被覆蓋度明顯低于其它區(qū)域，個別區(qū)域植被覆蓋度甚至低于0.2；到2019年，研究區(qū)整體植被覆蓋度略有改善，東北方向(2000年以前開采)植被覆蓋度也明顯提高，植被生長狀況明顯改善。東北方向地表植被覆蓋度較差，可能是由于該區(qū)距離城區(qū)較近，受人類活動影響較大；近年來，隨著生態(tài)保護和地表植被恢復工程的開展，到2019年，該區(qū)域植被覆蓋度有一定提高，但仍略低于其它區(qū)域。

表1給出了研究區(qū)2009年和2019年植被覆蓋度的統(tǒng)計結(jié)果。由表1可知：①近10年來，研究區(qū)中植被覆蓋度及其以上的面積占比均超過了88.86%，植被覆蓋總體情況良好；中高植被覆蓋度和高植被覆蓋度的面積均呈增加趨勢。②2009年中植被覆蓋度及其以下的面積占比為41.37%，所占比例遠高于2019年的20.60%，多出了4.70×105m2；而中高植被覆蓋度的面積在2009年占55.05%，到2019年，占比增加到了64.76%；高植被覆蓋度面積增加比例最為顯著，2019年的面積占比比2009年增加了308.89%，達到了14.64%。

表1 2009年和2019年各等級植被覆蓋度的面積及占比

2.2 地形因子對植被覆蓋度的影響

山區(qū)地形特征是一個多維變量，地表植被除了受降水、氣溫等氣象因子影響外，在小范圍區(qū)域，地形因子對植被覆蓋度空間分布的影響尤為重要。地形因子主要影響不同區(qū)域的光、熱、水、土等資源，進而影響植物生長。為探究礦區(qū)植被覆蓋度隨地形因子變化的空間分布規(guī)律，本文提取了高程、坡度、坡向等地形因子，并劃分等級，分析地形因子對植被覆蓋度空間分布的影響。

2.2.1 植被覆蓋度隨高程的變化特征

在一定高程范圍內(nèi)，山地存在氣溫垂直遞減和相對濕度隨高程升高而增加的規(guī)律；同時，隨著高程的變化，不同高度上接受到的太陽輻射、熱量，氣溫以及土壤性質(zhì)等也會發(fā)生變化，進而影響地表植被覆蓋度的分布[31]。表2為研究區(qū)高程分級統(tǒng)計表，在[766 m,1 130 m]范圍內(nèi)，不同高程所占面積比例隨高程的增加而呈先上升后降低的趨勢；研究區(qū)高程主要分布在(950 m,1 100 m]，占總面積的67.99%。

表2 研究區(qū)高程分級統(tǒng)計表

圖3給出了植被覆蓋度與高程的疊加分析結(jié)果。由圖3可以看出：2019年所有高程范圍內(nèi)的植被覆蓋度均高于2009年的；隨著高程的增加，植被覆蓋度逐漸升高，達到900 m時逐漸平穩(wěn)，1 050 m之后又逐漸降低，(900 m,1 050 m]高程范圍內(nèi)的區(qū)域是植被生長狀況良好的區(qū)域。這是因為低高程區(qū)，雖然水分、熱量和氣溫都適合植被的生長需要，但該區(qū)域受人類活動干擾大；隨著高程升高，人類活動影響減少，熱量和水分仍適合植被生長，植被覆蓋度開始上升，在(900 m,1 050 m]高程范圍內(nèi)，植被生長最好；當高程大于1 050 m后，隨著高程升高，氣溫明顯下降，熱量成為限制植被生長的主要因素，則植被覆蓋度下降[22]。

圖3 2009年和2019年植被覆蓋度隨高程的變化

2.2.2 植被覆蓋度隨坡度的變化特征

坡度表示局部地表的傾斜程度，直接影響物質(zhì)能量的轉(zhuǎn)換方式與程度，與熱量、水分條件、土壤厚度等因素關系密切，可對土壤水分和養(yǎng)分的分布進行調(diào)控，改變土壤的基本屬性，在很大程度上影響地表植被的分布[32]。研究區(qū)坡度分級統(tǒng)計情況見表3。由表3可以看出：研究區(qū)地表坡度主要在(8°,35°]范圍內(nèi)，該坡度范圍的面積占研究區(qū)總面積的比例為83.14%；而坡度為(35°,70°)的陡坡區(qū)域面積只占研究區(qū)總面積的7.12%，說明整個研究區(qū)的地勢較為平坦，起伏不大。

表3 研究區(qū)坡度分級統(tǒng)計表

圖4給出了研究區(qū)植被覆蓋度與坡度的疊加分析結(jié)果。由圖4可以看出：2019年研究區(qū)的植被生長狀況在所有坡度區(qū)均優(yōu)于2009年的；坡度較低區(qū)域的植被覆蓋度變化趨勢較為平緩，坡度大于8°后逐漸升高，在坡度為(25°,35°]時植被覆蓋度達到最大值，2019年的植被覆蓋度最高達到了0.73；坡度大于35°之后，植被覆蓋度開始降低。

圖4 2009年和2019年植被覆蓋度隨坡度的變化

總體來講，除坡度為(35°,70°)的陡坡區(qū)之外，坡度越大，植被覆蓋度越高。這是由于地勢平坦區(qū)域是人類生產(chǎn)活動影響較大的區(qū)域，植被覆蓋度較?。浑S著坡度增加，坡度較大區(qū)域受人類活動影響較小，該區(qū)域土壤持水能力和營養(yǎng)成分都適合植被生長，尤其有利于草地、灌木等覆蓋度較高的植被生長；而坡度為(35°,70°)的地區(qū)，地勢過于陡峭，水土流失嚴重，土壤養(yǎng)分也不易保存，植物生長受限。

2.2.3 植被覆蓋度隨坡向的變化特征

坡向表示每一個柵格高程值變化量的改變方向，主要影響太陽輻射和土壤含水量，從而改變植被的生長和分布特征[22]。表4為研究區(qū)坡向分級統(tǒng)計表。由表4可以看出：除南坡的面積占比僅6.32%以外，其余坡向的面積占比為10.02%～16.71%，各個坡向上的分布面積比較均勻。圖5給出了研究區(qū)植被覆蓋度與坡向的疊加分析結(jié)果。

表4 研究區(qū)坡向分級統(tǒng)計表

圖5 2009年和2019年植被覆蓋度隨坡向變化

由圖5可以看出：2019年植被生長狀況在各個坡向上均優(yōu)于2009年的；2009年和2019年各坡向植被覆蓋度最大值與最小值差值均為0.06，各坡向間差異不大，總體表現(xiàn)出東坡、南坡和東南坡植被生長狀況最差，西北坡生長狀況最好的特點。這是由于研究區(qū)屬半濕潤大陸性季風氣候，年降雨量較少，僅450～550 mm，東坡、南坡和東南坡方向雖光熱條件充足，但蒸發(fā)量過大，易使土壤水分含量不足；西北坡接受太陽輻射較少，土壤水分易于保持，光熱條件適宜，更有利于植被生長。因此，西北坡植被生長狀況優(yōu)于東坡、南坡和東南坡的[29]。

2.3 植被覆蓋度的時空演變及其驅(qū)動因素

2.3.1 植被覆蓋度的時空變化趨勢

為揭示近10年來植被生長狀態(tài)的空間變化特征，利用ArcGIS軟件對2009年和2019年的植被覆蓋度圖進行柵格疊加運算，得到研究區(qū)近10年來植被覆蓋度的時空變化分布圖，如圖6所示。從圖中可以看出，2009—2019年研究區(qū)內(nèi)植被覆蓋度的變化整體上呈增加趨勢。表5為植被覆蓋度變化等級的面積統(tǒng)計，近10年來，植被覆蓋度增加的區(qū)域面積占研究區(qū)總面積的81.70%，植被覆蓋度降低的區(qū)域面積僅占總面積的18.30%，植被覆蓋度增加的區(qū)域面積遠高于降低區(qū)域的面積，且以[0.0,0.2]的植被覆蓋度變化等級所占面積比例最大，達到了71.64%。結(jié)果表明，近10年來，研究區(qū)植被覆蓋度整體呈上升趨勢，植被生長狀況明顯改善。

圖6 2009—2019年研究區(qū)植被覆蓋度的時空變化

表5 2009—2019年植被覆蓋度變化的面積統(tǒng)計

為探究開采時間對近10年來植被覆蓋度變化的影響，從圖6中提取了2009—2019年不同開采年份工作面的植被覆蓋度變化量；同時，為降低不同時期影像的誤差影響，提取的不同開采年份工作面的植被覆蓋度變化量需同時減去對照區(qū)的變化量。圖7為不同開采年份工作面2009—2019年植被覆蓋度變化量。由圖7可知：植被覆蓋度的變化以2007年為界限分為兩段，開采年份早于2007年的煤礦，除1984年、1990年、1998—2000年開采工作面的植被覆蓋度變化量為負值外，其余開采年份植被覆蓋度變化量均大于0，而2007年后開采的所有煤礦開采工作面的植被覆蓋度變化量均小于0；與其它年份相比，1984年開采工作面的植被覆蓋度變化明顯異常。這可能是由于其工作面距城區(qū)較近，受人類活動影響較大，且工作面面積較小，也增大了分析誤差。綜上可知，煤礦開采后，雖然短期內(nèi)會對地表植被生長狀況造成影響，但自然生態(tài)系統(tǒng)的自我調(diào)節(jié)能力輔以人工保護措施，會使生態(tài)系統(tǒng)逐漸向良性循環(huán)方向發(fā)展，植被生長狀況也會不斷得到改善[16-18]。

圖7 2009—2019年不同開采年份工作面植被覆蓋度變化量

2.3.2 植被覆蓋度影響因子的地理探測

由影響因子的探測結(jié)果發(fā)現(xiàn)，各因子q值排序為：開采時間>高程>坡向>坡度。開采時間是影響植被覆蓋度時空變化的最主要因素，其q值為0.14；高程為次要影響因素，其q值為0.12；坡度和坡向雖然是影響植物生長的重要因素，但這2個因子的q值均低于0.05，且坡度的q值最低，說明2009—2019年植被覆蓋度變化主要受開采時間和高程的影響，坡度和坡向的直接影響較小，這與前人的研究結(jié)論一致[30,33]。

表6為不同影響因子的交互作用的探測結(jié)果，從表6中可以看出，不同影響因子的兩兩交互作用均大大增加了單因子對植被覆蓋度變化的解釋力。交互作用q值排序為：開采時間∩高程>高程∩坡向>開采時間∩坡向>高程∩坡度>開采時間∩坡度>坡度∩坡向。開采時間、高程和坡向的兩兩交互作用較強，q值均在0.33以上，坡度與坡向的交互作用q值最低，僅0.18。

表6 不同影響因子的交互作用的探測結(jié)果

3 結(jié)語

1)2009年和2019年，研究區(qū)植被生長總體情況良好，中植被覆蓋度及其以上面積占比均超過了88.86%，中高植被覆蓋度和高植被覆蓋度面積均呈增加趨勢；空間分布上，東北方向2000年前開采工作面的地表植被覆蓋度明顯低于其它區(qū)域。

2)從地形因子對植被覆蓋度的影響分析發(fā)現(xiàn)，植被覆蓋度總體隨高程和坡度增加呈現(xiàn)先增加后減小的趨勢，在高程(900 m,1 050 m]和坡度(25°, 35°]的地帶達到最大值；在坡向上呈現(xiàn)出西北坡最高，而東坡、南坡和東南坡相對較低的趨勢，且高程對植被覆蓋度分布的影響大于坡度和坡向。

3)對比2009—2019年植被覆蓋度的時空變化發(fā)現(xiàn)，2007年以前開采的大多區(qū)域植被覆蓋度有明顯改善，而2007年后開采的區(qū)域植被覆蓋狀況呈退化趨勢，表明雖然短期內(nèi)地下采煤會對植被覆蓋度造成影響，但生態(tài)系統(tǒng)的自我調(diào)節(jié)能力再輔以人工保護措施，會使生態(tài)系統(tǒng)逐漸向良性循環(huán)方向發(fā)展，植被生長狀況逐步改善。

4)因子探測結(jié)果同樣表明，開采時間是影響研究區(qū)植被覆蓋度時空變化的最主要因素，高程次之，坡度的q值最低；交互作用探測中，以高程和開采時間的共同作用最強，說明人為因素與自然因素的協(xié)同能夠增強對植被覆蓋度變化的影響力。