神東礦區(qū)植被覆蓋度時序變化與驅動因素分析及引導恢復策略

劉 英，雷少剛，陳孝楊，陳 敏，楊英明，李心慧，張旭陽，龍林麗，卞正富

(1.安徽理工大學 地球與環(huán)境學院，安徽 淮南 232001；2.中國礦業(yè)大學 礦山生態(tài)恢復教育部工程研究中心，江蘇 徐州 221116；3.深部煤礦采動響應與災害防控國家重點實驗室，安徽 淮南 232001；4.煤炭水資源保護與利用國家重點實驗室，北京 102209)

我國重點建設的14個大型煤炭基地主要集中在生態(tài)環(huán)境脆弱、水土流失嚴重的黃土高原以及新疆等地區(qū)[1]。根據煤炭工業(yè)“十三五”發(fā)展規(guī)劃，預計截止到2020年煤炭占我國一次能源消費中的比例仍占近60%，立足于國內是我國能源安全戰(zhàn)略要求，煤炭作為能源主體地位和壓艙石作用在很長一段時間內不會變化[2]。由于地理位置與氣候原因，西部干旱半干旱地區(qū)水資源占有量僅占我國的8.3%，水資源十分匱乏，生態(tài)環(huán)境十分脆弱[3-4]。加之國家煤炭開采戰(zhàn)略重心西移及集群化、高強度的開采方式，使得西部礦區(qū)的地表生態(tài)環(huán)境呈現出強烈擾動的態(tài)勢，使本就脆弱的生態(tài)環(huán)境日趨退化[5]。西部煤炭資源開發(fā)引發(fā)的生態(tài)環(huán)境約束依舊是社會各界關注焦點[6]，大規(guī)模的煤炭開采引起地表沉陷，影響塌陷區(qū)地形、水文、土壤等植物生境因子，從而導致植物生境改變、影響植被生長[7-8]。

當前我國西部礦區(qū)已開展了諸多且富有成效的植被恢復重建研究，形成了以表土改良、地表裂縫治理、植物根際環(huán)境改善以及菌根共生等為主的植被恢復與重建技術[9-11]。但是這些技術成本往往過高，大面積種植新進的喬、灌類植物對本地植物物種形成沖擊，大量開挖魚磷坑亦會引起采后生態(tài)系統(tǒng)重復擾動，且部分重建植物難以適應惡劣的氣候地理條件，管護周期長，成活率低，生長緩慢甚至衰亡，現有西部礦區(qū)重建植被可持續(xù)性不容樂觀，其生態(tài)恢復效益以及對生態(tài)環(huán)境影響還有待商榷，形成了植被恢復資金高投入與恢復低效率的矛盾局面[12]。另一方面，煤炭開采地表塌陷對植被的影響并非全是負面的，部分區(qū)域的植被覆蓋度、生物多樣性呈增長態(tài)勢[13]。所以當務之急，應嚴格遵循當地自然生態(tài)系統(tǒng)發(fā)展規(guī)律，避免對采后生態(tài)系統(tǒng)的再次擾動，科學配置，優(yōu)化布局、引導恢復生態(tài)系統(tǒng)的自修復能力[14]。弄清煤炭開采擾動區(qū)植被覆蓋度時序變化，是探索半干旱礦區(qū)植被引導型恢復有效方法重要基礎與前提。

對于半干旱礦區(qū)受損植被引導型恢復來說，需要面臨：植被在哪種破壞程度下可以實現自恢復?當需要人工引導干預時，在什么地方干預?怎么干預?干預到何種程度?解決上述基本問題，首先需要對礦區(qū)植被覆蓋度時序變化及驅動因素進行分析，這是采煤擾動礦區(qū)植被受脅迫以及采后植被恢復狀況空間感知、采后恢復區(qū)植被覆蓋度變化趨勢和年際波動程度評估的基礎，可以為礦區(qū)有序、精準進行植被引導恢復工作提供依據。筆者將基于Landsat-NDVI時序數據，以神東中心礦區(qū)為研究尺度，分析采區(qū)和非采區(qū)植被覆蓋度的時序演變規(guī)律，并依據不同時序植被覆蓋度格點回歸斜率、相關系數及標準差，對研究區(qū)不同采后恢復區(qū)植被覆蓋度變化趨勢以及年際波動程度進行評估，摸清研究區(qū)植被生長狀況，為礦區(qū)有序、有針對地進行植被恢復工作提供基礎數據。

1 研究區(qū)域概況與數據來源

1.1 研究區(qū)域概況

神東煤田中心礦區(qū)位于鄂爾多斯高原東南部及陜北黃土高原北緣和毛烏素沙漠的東南邊緣(圖1)，地理坐標為東經109.83°～110.34°，北緯39.56°～39.19°，面積約900 km2。該區(qū)年平均氣溫6.2 ℃，極端最低氣溫-31.4 ℃，極端最高氣溫36.6 ℃。年降水量為300～400 mm，主要集中在夏季，蒸發(fā)量大于降雨量4倍以上，屬于典型的干旱、半干旱的高原大陸性氣候。地處干草原與森林草原的過渡地帶，主要植被類型為干草原、落葉闊葉灌叢和沙生類型植被，主要植物品種有沙柳、楊樹、長芒草、芨芨草、油松、檸條、沙棘、紫花苜蓿、沙竹、油蒿等。地貌地形西北高，東南低，平均海拔在+1 200 m左右。礦區(qū)東部及東北部為黃土丘陵山區(qū)，區(qū)內溝壑縱橫，形成梁峁、溝壑和土塬3種地貌，侵蝕活躍，水土流失嚴重；西部及西南部以風沙地為主的流動、半固定及固定沙地上，分布著沙地植被，主要是沙地先鋒植物群落和油蒿群落；在洼地、灘地和湖泊周圍分布有濕地植被。

圖1 研究區(qū)域位置Fig.1 Study area geographic location

1.2 數據來源

選取2000-01—2016-12月的共337景Landsat 5 TM和8 OLI遙感影像數據，時間分辨率為16 d，空間分辨率為30 m，影像下載列/行號為127/33。所有Landsat影像均從美國地質調查局(USGS)網站(https://glovis.usgs.gov/)下載。由于遙感數據不可避免的而受到云層覆蓋的影響，數據選取時剔除云量大于5%的Landsat影像。為了提高精度和減少地形和大氣的影響，利用軟件ENVI 5.1(Exelis Visual Information Solutions USA)對遙感影像進行預處理。所有獲取的Landsat TM和OLI遙感影像數據經輻射定標、大氣校正、配準、裁剪等預處理，并根據研究區(qū)域的矢量邊界進行裁剪。所有影像統(tǒng)一采用UTM投影，WGS84坐標系。所有影像以2000年1月份的影像為基準配準，以保證影像之間誤差在1個像元以內。

1.3 礦區(qū)開采工藝及開采情況

該礦區(qū)煤炭資源賦存條件好，目前開采的礦井主要包括大柳塔礦、補連塔礦、哈拉溝礦、上灣礦等，以大采高一次采全高綜合機械化采煤為主要采煤方法，300 m寬工作面較為普遍，走向長度1 500～5 220 m不等，采高通常在2.0～7.0 m，由于綜合機械化程度高，平均采速每日可推進10.0 m以上，有的甚至達到20 m/d，地表最大下沉值在2 320～5 887 mm，最大水平位移可達3 284 mm，最大下沉速度可達837 mm/d，地表綜合邊界角、超前影響角與綜合邊界角一般分布在50°～70°，綜合移動角一般為60°～70°，裂縫角一般為65°～85°，研究區(qū)部分工作面屬于初采，大量工作面屬于重復采動。煤炭資源大規(guī)模、高強度開采引起的地表沉陷具有下沉速度快、下沉量集中，造成地表破損嚴重，地裂縫、塌陷坑等非連續(xù)形變發(fā)育等特點。

2 研究方法

2.1 植被指數與覆蓋度的計算方法

歸一化植被指數(Normalized Difference Vegetation Index,NDVI)可描述植被的生長狀況，是定量評價植被蓋度常用指標之一[15]。其計算時由于植物葉片組織對綠光，尤其是近紅外光有很強的反射，對藍光(470 nm)和紅光(650 nm)有很強的吸收。因此，在Landsat影像中提取植被信息時多采用紅外與紅波段圖像反射值比值，這樣可消除部分輻射誤差。其計算公式為

(1)

式中，INDV為歸一化植被指數；NNIR,NR分別為近紅外和紅光帶的實際表面反射率。

植被覆蓋度(Vegetation Coverage Fraction，VCF)指某一地域植物垂直投影面積與該地域面積的比值。目前應用最廣泛的是利用歸一化植被指數(如INDV)，根據像元二分模型計算得到植被覆蓋度指數[16]，其計算公式為

(2)

(3)

(4)

2.2 植被覆蓋度時序變化趨勢計算方法

依據神東礦區(qū)2000—2016年時序植被覆蓋度數據，以格網中包含像素植被覆蓋度的平均值，通過對每一個像元進行一元線性回歸方程擬合，得到每個像元回歸直線的斜率θslope，當θslope<0時，表示植被覆蓋度呈減少趨勢；當θslope>0，表示植被覆蓋度呈升高趨勢。采用擬合相關系數R來判定變化趨勢顯著性，當R<0時，表示植被覆蓋度呈減少的趨勢；當R>0，表示植被覆蓋度呈升高的趨勢[18]。用t分布檢驗2個變量是否真正相關及其相關性的顯著水平[19]。計算公式為

(5)

(6)

其中，θslope為一元線性回歸方程斜率；ti為年份序號，i=1，2，…，15；xi為第i年生長季植被覆蓋度均值；R為相關系數；n為總年份數，在本研究中n=15。依據式(5)計算得到每個像元回歸直線的斜率θslope，根據t檢驗結果對應的擬合相關系數R，將變化趨勢及顯著性分成10個等級，即：①θslope<0，0.010，0.10，0.050，0.0250，0.010，0.01

2.3 植被覆蓋度時序波動程度計算方法

利用時序變異系數(Standard Deviation，SD)來表示植被時序波動程度，時序變異系數值越大，說明植被覆蓋度年際波動越大，受外界干擾程度越嚴重；反之，則植被覆蓋度年際波動越小，受干擾程度低。其計算公式為

(7)

其中，DS為植被時序變異系數。結合式(5)得到的θslope，將植被覆蓋度波動程度按正負分成10個等級，即：①θslope<0，DS≥0.4，劃分為“極高幅度降低”；②θslope<0，DS∈[0.30,0.40)，劃分為“高幅度降低”；③θslope<0，DS∈[0.15,0.30)，劃分為“中幅度降低”；④θslope<0，DS∈[0.10,0.15)，劃分為“低幅度降低”；⑤θslope<0，DS∈[0,0.10)，劃分為“極低幅度降低”；⑥θslope>0，DS∈[0,0.10)，劃分為“極低幅度升高”；⑦θslope>0，DS∈[0.10,0.15)，劃分為“低幅度升高”；⑧θslope>0，DS∈[0.15,0.30)，劃分為“中幅度升高”；⑨θslope>0，DS∈[0.30,0.40)，劃分為“高幅度升高”；⑩θslope>0，DS≥0.4，劃分為“極高幅度升高”。

3 礦區(qū)植被覆蓋度變化及驅動因素分析

3.1 基于Landsat的神東中心礦區(qū)植被NDVI時序變化分析

由于各工作面開采時間不同，主要選取各礦井2014,2011,2009,2007,2004和2001年開采的典型工作面區(qū)域為采區(qū)，以神東中心礦區(qū)為對照區(qū)，對2000—2016年間植被NDVI時序變化進行統(tǒng)計(圖2)，圖中紅色箭頭為不同煤炭開采年份。從整體時序變化來看，2000—2016年神東中心礦區(qū)植被INDV整體呈現出明顯物候性周期變化。植被INDV的最大值出現在每年7月下旬至8月中旬，最小值出現在每年12月份到1月份。從植被INDV年際變化來看，2000—2005年植被INDV變化幅度變化較小，總體在0.28左右；2005—2008年植被INDV小幅上升，然后2008—2010植被INDV呈下降趨勢，相比于2008年峰值降低了近19.15%；2010—2014年植被INDV出現明顯上升，相較于2010年峰值增加了近31%，且以2014年的峰值最大；2014—2016年植被INDV有小幅降低，整體在0.40左右波動。但是從整個時序區(qū)間來看，神東中心礦區(qū)植被INDV呈明顯的增大趨勢，這說明植被覆蓋情況總體上有所好轉。原因在于：地方政府和神華集團進行了大面積礦區(qū)植被重建，完成水土保持以及綠化、大面積防風固沙和多項小流域水土保持綜合治理措施，對整個礦區(qū)植被覆蓋度有明顯提升，由圖2可以看出，截止到2016年，神東中心礦區(qū)植被INDV從2000年0.20左右提高到0.45～0.60。

圖2 2000—2016年神東中心礦區(qū)與不同開采年份工作面INDV(30 m)時序變化Fig.2 Time-series changes of INDV (30 m) of Shendong central mining area and working face in different mining years from 2000 to 2016

由于研究區(qū)不同開采年限不同，由此引起的植被INDV變化具有差異性，為此將不同開采年限工作面時序INDV求平均值，減去中心礦區(qū)時序INDV，得到不同開采年限工作面INDV差異值時序變化情況。統(tǒng)計發(fā)現，INDV差異值時序變化與該工作面煤炭開采年份密切相關。煤炭開采前工作面地表INDV約等于甚至大于中心礦區(qū)INDV，INDV差異值均在0附近上下波動，變化幅度為±0.03；時序INDV差異值大于0，說明該工作面采前INDV高于作為對照的中心礦區(qū)INDV，這主要由于這些被統(tǒng)計的工作面主要分布在大柳塔礦井硬梁地覆蓋區(qū)，特別是5-2煤層開采區(qū)，植被長勢相對較好，分布著大面積的楊樹、檸條、沙柳、油蒿、長芒草、百里香、沙棘等，而中心礦區(qū)還包括大面積風沙覆蓋區(qū)，這些地區(qū)植被多以低矮沙生草本、灌木為主，其植被覆蓋度相對較低。隨著煤炭開采活動的推進，INDV差異值逐漸增大，且小于0，說明該煤炭開采工作面INDV小于中心礦區(qū)INDV，從工作面采前到采后INDV差異值從“大于0”或者“0附近波動”轉變“小于0”，對應的采區(qū)非采區(qū)INDV表現為“采區(qū)INDV不小于等于非采區(qū)”轉變?yōu)椤安蓞^(qū)INDV小于非采區(qū)”，這體現了隨著采煤擾動影響地表植被INDV相對于對照區(qū)降低的變化趨勢。通過對采后INDV差異值時序變化趨勢進行擬合發(fā)現：采后5 a，工作面地表INDV差異值呈增大趨勢；采后7 a，工作面地表INDV差異值開始降低；采后12 a工作面地表INDV差異值逐漸趨向于0。說明采后5 a內，相對于對照區(qū)，工作面采區(qū)植被INDV的變化表現為持續(xù)降低的過程；采后7 a，受擾動區(qū)植被開始恢復，工作面地表INDV與對照區(qū)INDV差異值開始降低；至采后12 a，工作面采區(qū)植被INDV基本能夠恢復至對照區(qū)水平。

3.2 神東中心礦區(qū)植被覆蓋度時序變化趨勢分析

采用式(5)，(6)計算得到不同格點植被覆蓋度一元線性回歸斜率θslope(圖3(a))和對應的相關系數R(圖3(b))，利用ArcGIS疊置分析功能對計算結果進行疊加得到神東礦區(qū)植被覆蓋度空間變化趨勢圖(圖3(c))。

對圖3結果進行統(tǒng)計得到：隨著時間的推移，神東中心礦區(qū)植被覆蓋度均值整體呈現小幅增加趨勢，植被覆被度呈升高與降低的區(qū)域面積分別占中心礦區(qū)總面積的72.35%和27.65%，其中θslope>0.01的面積占研究區(qū)總面積的52.03%，研究區(qū)域植被覆蓋度以正向變化為主，由此可見，煤炭開采引起的地表塌陷對植被覆蓋度的影響并不全是負面的。θslope<0，且R∈[-1,-0.89)內，通過了0.01顯著性檢驗，植被受損最為嚴重，植被覆蓋度極顯著性降低，面積約為4.81 km2，占研究區(qū)域面積的0.62%。θslope<0，且R在[-0.89,-0.71)內，通過了0.025顯著性檢驗，植被受損較為嚴重，植被覆蓋度較極顯著降低，其面積約為9.39 km2，占整個研究區(qū)域面積的1.21%。θslope<0，且R在[-0.71,-0.58)內，通過了0.025顯著性檢驗，植被覆蓋度顯著降低，其面積約為20.89 km2，占整個研究區(qū)域面積的2.69%。θslope<0，且R在[-0.58,-0.32)內，通過了0.05顯著性檢驗，植被覆蓋度較顯著降低，其面積約為108.26 km2，占整個研究區(qū)域面積的13.92%。值得注意的是，還有大約占整個中心礦區(qū)面積的20.78%，161.69 km2的區(qū)域地表植被覆蓋度變化θslope在[-0.01,0.01]內，盡管整體上神東中心礦區(qū)植被覆蓋度呈上升趨勢，但是在煤炭開采活動最為集中、擾動最為頻繁的中心礦區(qū)植被覆蓋度升高幅度較低，甚至有大面積的植被退化區(qū)域。上述植被覆蓋度降低的區(qū)域主要在神東中心礦區(qū)各礦井工作面上覆地表小區(qū)域成片分布。因此，在礦山植被修復重建時，有針對性的對植被退化區(qū)域進行識別，然后對識別出的退化區(qū)域植被進行引導修復，可以降低礦山植被修復的成本，提高植被重建的效率，這也是引導型恢復提出的立足點。

圖3 神東中心礦區(qū)植被覆蓋度時序變化趨勢Fig.3 Time-series variation trend of vegetation coverage in Shendong central mining area

3.3 神東中心礦區(qū)植被覆蓋度時序波動程度分析

采用式(7)計算得到不同格點植被覆蓋度時序波動程度(圖4(a))，并利用ArcGIS疊置分析功能對計算結果與一元線性回歸斜率θslope的正負情況(圖3(a))進行疊加，得到神東礦區(qū)植被覆蓋度正向、負向波動情況(圖4(b))。

中心礦區(qū)年際間植被覆蓋度變異系數(DS)處于0.04～0.46 (圖4(a))，根據劃分的植被覆蓋度波動等級，極低幅度、低幅度、中幅度、高幅度和極高幅度波動區(qū)的面積分別占中心礦區(qū)總面積的1.11%，34.71%，56.71%，6.09%和1.38%，以中、低幅度變化為主。結合θslope的正負，研究區(qū)域東北部、中部、南部大面積地區(qū)植被波動等級為低幅度降低，該地區(qū)為荒漠化嚴重地區(qū)，植被覆蓋度本來就很低，加之煤炭井工開采地下水會被抽空，影響地下水對地表水的補給，植被生長的立地條件會受到極大程度的破壞，進而影響地表受擾動區(qū)植物的生長與恢復，導致該區(qū)域植被生長與恢復狀況較差。中心礦區(qū)東南部、西南部、西北部年際間植被覆蓋度變化主要以中幅度變化為主，結合區(qū)域θslope的正負，植被覆蓋度變化表現為中幅度升高，這與神東礦區(qū)大力開展植被恢復活動，先后建立完善了專門水保綠化機構，制定了礦區(qū)防風固沙與水土保持規(guī)劃，以林草植被建設為主，完成礦區(qū)植被重建以及大面積常規(guī)防風固沙和多相小流域水土保持綜合治理工程，致力于提高該地區(qū)植被覆蓋度與多樣性有關[3]，如，在大柳塔、補連塔、哈拉溝等采煤塌陷區(qū)建立了植被修復重建應用的國家級示范區(qū)。此外，居民區(qū)、城市道路與河道兩邊以及煤礦周邊輔助設施建設區(qū)域，植被受人為活動影響較大，植被覆蓋度變化幅度較大。

圖4 神東中心礦區(qū)植被覆蓋度時序波動程度Fig.4 Time series fluctuation degree of vegetation coverage in Shendong central mining area,coefficient of variation

3.4 植被覆蓋度變化驅動因素分析

3.4.1植物覆蓋度與土壤理化性質關系分析

以大柳塔礦5-2和2-2部分工作面作為實驗樣地，選取4個采煤工作面，共建立18個植被調查樣方，在每個樣方內隨機建立10 m×10 m植被群落學調查大樣方，并利用5點取樣法在樣地區(qū)域內設置5個2 m×2 m小樣方，并在每個小樣方利用對角線法設立1個1 m×1 m基本調查單元，同時也作為土壤理化性質數據取樣的網格單元，然后利用Canoco 4.5軟件大柳塔礦井18個樣方植被與土壤理化性質典范對應關系進行分析。主要考慮各樣方中與植被生長密切相關的土壤含水量(SWC)、有機質(OM)、容重(BD)、pH值、全氮(TN)、總磷(TP)、全鉀(TK)、總有機碳量(TOC含量)、土壤電導率(EC)及各樣方中24種主要建群種植物、個體數量與蓋度，分析結果如圖5所示，圖5中紅色、藍色、青色分別為土壤理化性質、植物名稱和樣方編號，紅色剪頭線段的長短代表各土壤理化性質因子對不同植物及樣方的影響大小，紅色與藍色箭頭線段之間的夾角代表相關性大小，銳角表示呈正相關，鈍角則表示呈負相關，夾角越小相關性越大。

圖5 植被-土壤理化性質典范對應關系Fig.5 Canonical correspondence between vegetation and soil properties

從植被與土壤理化性質典范對應關系分析結果來看，第3象限上的因子可以解釋67.2%的植被分布狀況。土壤含水量SWC與第1排序軸(橫軸)第3象限相關系數達到0.727 2，與第2排序軸(橫軸)第3象限相關系數達到0.746 7，說明SWC為影響植被生長的第1重要土壤理化性質因子。從植被樣方在排序圖上的分布，可以看出，有12個樣方主要分布在第3象限，并且24種主要植物建群種有13種分布在第3象限，如：紫花苜蓿、楊樹、檸條、華北白前、沙棘、斜莖黃耆、狗娃花、胡枝子、地角兒苗、芨芨草等。此外，土壤有機質(OM)與第1排序軸相關系數為0.649 4，表明其也是影響植被分布較為重要的因子，其他重要影響因子還包括容重(BD)、全氮(TN)等。諸多研究表明采煤沉陷裂縫破壞了土體結構，增加了土壤水分的蒸發(fā)面，是裂縫區(qū)土壤水分降低的主要原因[12,20-21]。

3.4.2植物覆蓋度與地下水位埋深關系分析

植被生長過程中，水資源扮演著舉足輕重的角色，而對于生長在半干旱區(qū)的植被而言，區(qū)域內降水量低，因此地下水位埋深對其生長狀態(tài)起著至關重要的作用，是控制半干旱礦區(qū)地表植被生長狀況的又一關鍵決定因素[7]。有諸多學者對我國半干旱地區(qū)植物覆蓋度與地下水位埋深關系進行了研究。范立民團隊較早關注到神東礦區(qū)煤炭資源大規(guī)模開采與地下水資源破壞這一矛盾，指出該礦區(qū)煤層淺埋區(qū)采煤會引起大面積地下水漏失、減少河流基流量，從而破壞植被生態(tài)系統(tǒng)[22]。徐海量等[23]對我國西部干旱地區(qū)不同地下水位對植被的影響研究后認為3.5 m為該地區(qū)草本植被生態(tài)恢復的最低水位閾值。李丹[24]利用Gamma分布函數對位于鄂爾多斯高原中部的海流兔河流域生長季多年平均INDV隨地下水埋深的變化，發(fā)現地下水埋深在小于10.0 m時，INDV隨地下水埋深增大而減??；地下水埋深大于10.0 m時，INDV與地下水埋深基本沒有關系。金曉媚等[25]對海流兔河流域植被分布進行研究時認為，地下水位埋深對植被的影響范圍是1.0～5.0 m，當地下水位大于5.0 m 時，決定植被的主要因素就變成了降雨量和土壤含水量。由此可見，不同研究區(qū)地下水埋深對植被的影響閾值不同，但是地下水埋深對植被的影響閾值總體處于0～10.0 m。

對于神東中心礦區(qū)，根據卞正富[12]、雷少剛[3]等對神東礦區(qū)地下水埋深與植物群落關系長期現場調查結果，認為地下水埋深3.2 m和8.0 m為影響植被INDV的2個重要閾值(圖6)。

圖6 地下水位埋深與植被群落類型的關系 (依文獻[26]補充)Fig.6 Relationship between groundwater depth and vegetation community types (supplemented by reference[26])

地下水埋深在0～3.2 m內，對根系較淺的濕生植被分布影響最大，隨著地下水位埋深的增加，與濕生植被INDV相關性降低，此時土壤水與地下水之間作用強烈，土壤含水量飽和與非飽和狀態(tài)交替頻繁，濕生植被的生長與地下水的變化密切相關；當地下水埋深在3.2～8.0 m時，這些根系較淺的濕生植被與地下水位埋深之間不再有顯著的聯系，取而代之的是一些根系較長的旱生植被，此時地下水和土壤水之間雙向聯系，地下水仍可以通過蒸發(fā)與毛細管作用對土壤水進行補給；當地下水埋深大于8.0 m時，植被生長主要依靠的是少量存儲在土壤中的大氣降水和大氣凝結水，而與地下水埋深之間的關系不再顯著，地表植被覆蓋格局將進一步由旱生植被將演替為沙生植被[3,26]。

3.4.3植被NDVI時序變化與降雨量、溫度關系分析

氣候條件(降雨量、溫度)是決定地表植被生長和分布的重要因素[27]。利用位于礦區(qū)的大柳塔氣象站收集到的研究區(qū)域2001—2016年的氣象數據，分析月平均溫度和降雨量與神東中心礦區(qū)月植被INDV的相互關系(圖7)。神東中心礦區(qū)植被INDV與月平均氣溫和降雨量具有相似的年周期變化，如圖7所示，這也使得植被INDV與月平均氣溫(R2=0.640 1，P<0.05)和降雨量(R2=0.651 1，P<0.05)具有較高的相關性。從2000—2016年，神東中心礦區(qū)降雨量和INDV均呈增加趨勢(線性擬合斜率分別為0.004 8,0.000 38)，而溫度相對較穩(wěn)定(線性擬合斜率為-0.009 2)，且INDV和平均降雨量在每年6—8月份會有較大幅度的升高，溫度在此期間平均氣溫卻無明顯差異，根據圖7，2007—2016年間中心礦區(qū)植被INDV的升高、INDV差異值的變化與降雨量的增加也極具同步性，說明降雨量是植被INDV波峰值的最敏感氣候因子，對中心礦區(qū)植被覆蓋度的變化有極其重要的影響。

圖7 神東中心礦區(qū)植被INDV與月平均溫度、降雨量的變化Fig.7 Changes of INDV and monthly mean temperature and rainfall in Shendong mining area

3.4.4采礦與土地損毀植被覆蓋的影響傳遞過程分析

(1)采空區(qū)上方巖層破斷與運動過程。井工煤炭開采在地底形成大面積采空區(qū)，破壞巖層應力的平衡狀態(tài)，引起巖體內部應力重新分布，導致巖層破斷和運動，但是不同類型巖層的破壞變形又具有差異性、非同步的特點。隨著工作面推進軟巖層破斷和運動形式表現為同步斷裂，而對于硬巖層，其巖層破斷和運動形式表現為間斷式斷裂。硬巖層破斷移動后形成的塊體重新分布，勢必會改變采空區(qū)整個上覆巖層初始狀態(tài)，因此，硬巖層(關鍵層)破斷形成的塊體主導著巖層采空區(qū)上覆巖層運動。隨著采空區(qū)硬巖層破斷與關鍵層破斷為巖塊，則上覆巖層運動最終形成如圖8所示的關鍵層主導的采動覆巖運動概貌，通常在采空區(qū)上方形成垮落帶、斷裂帶和彎曲下沉帶。當關鍵層巖塊破斷后巖層運動重新趨于平穩(wěn)時，巖層裂隙將重新閉合，形成重新壓實區(qū)；對于巖層發(fā)生劇烈運動，礦壓、裂隙與塌陷劇烈變化區(qū)域稱為裂縫發(fā)育區(qū)[28]。

圖8 采空區(qū)上覆巖層運動分區(qū)[28]Fig.8 Movement zoning of overlying strata in goaf[28]

(2)巖層破斷與運動到土體沉降變形過程。隨著巖層破斷和運動過程之后是地表土體沉降變形過程，通常在采空區(qū)上方地表形成大范圍的塌陷盆地。根據開采沉陷學理論，塌陷盆地中不同區(qū)域的下沉程度不同，一般將下沉盆地劃分為拉伸變形區(qū)、壓縮變形區(qū)和中性區(qū)[29]?！叭齾^(qū)”的形成原因以及對土體結構的影響有所差異，對于下沉盆地中間區(qū)被稱為中性區(qū)，該區(qū)地表下沉值最大，一般先會出現臨時性裂縫，煤炭開采過程中隨工作面的推進同時發(fā)育，當工作面推過裂縫后，該區(qū)破斷后巖層運動重新趨于平穩(wěn)，大部分裂縫將逐步閉合，其相對移動和變形值較小；塌陷盆地內邊緣區(qū)被稱為壓縮變形區(qū)，位于采空區(qū)邊界附近到最大沉降點之間，地表土體移動向盆地中心傾斜，使得土體產生壓縮變形而形成擠壓型裂縫；塌陷盆地外緣區(qū)被稱為拉伸變形區(qū)，位于采空區(qū)邊界到盆地邊界之間，地表沉降程度不均勻使得地表土體移動向盆地中心方向傾斜，土體產生拉伸變形，當拉伸變形超過極限抗拉強度后，地表將產生拉伸型裂縫，其特點為發(fā)育寬度、深度大、難以自恢復，破壞原始土地結構[5]，對表土水分的保持以及植物的生長的影響較大，是塌陷區(qū)植被恢復需要重點采取工程措施加以整治的區(qū)域。

(3)土體沉降變形到植物生長立地條件變化過程。隨著采空區(qū)土體塌陷變形過程后，地表植物生境條件隨之發(fā)生改變，主要表現為地形地貌條件、地下水位埋深、SWC、土壤理化性質等的改變[5,12]。地下煤炭開采對SWC的影響在不同塌陷區(qū)位存在空間異質性。中性區(qū)位于地表下沉盆地的中心，沿開采工作面傾向上其地形往往呈“兩邊高，中間低”狀，該地形有利于地表徑流的匯集，導致中性區(qū)SWC高于壓縮區(qū)、拉伸區(qū)。壓縮區(qū)土體被擠壓，土壤孔隙度降低、硬度升高、入滲能力降低，不利于土壤水分保持，所以相較于中性區(qū)，壓縮區(qū)的SWC通常更低。拉伸區(qū)地表存在永久性裂縫，在一定程度上改變地表徑流方向和匯水條件，當地表裂縫與采空區(qū)直接連通時，地表徑流將直接流進采空區(qū)(圖9)。此外，永久性裂縫還會導致土壤水分的蒸發(fā)面增加，土壤水的蒸發(fā)散失速度加快，而裂縫周邊的土體受拉伸作用導致結構破壞，也會增大土壤的孔隙比，不利于土壤水分保持，地表水分流失進一步加重[20]。

圖9 采煤塌陷對植被生長立地條件變化過程示意Fig.9 Change process diagram of plant growth site conditions caused by mining

通過沉陷過程不同變形區(qū)位土壤含水量進行時序監(jiān)測后發(fā)現：土體塌陷變形對“三區(qū)”土壤含水量的影響總體表現為，中性區(qū)SWC>壓縮區(qū)SWC>拉伸區(qū)SWC(圖10)。采后1個月中心區(qū)、壓縮區(qū)、拉伸區(qū)土壤水開始出現了較為明顯的降低，尤其是永久性裂縫大量發(fā)育的拉伸區(qū)土壤水下降最為明顯，從采前的11.8%最低降至5.5%；壓縮區(qū)土壤含水量從采前的11.9%最低將至8.4%；中性區(qū)隨著臨時性裂縫產生與自恢復，土壤含水量呈波動狀，在8.9%～12.1%波動。SWC在煤炭開采地表不同沉陷區(qū)位的空間異質性也是沉陷區(qū)植被覆蓋度時空變化特征不一致的主要原因。

圖10 開采前后不同塌陷區(qū)位土壤水分時序變化Fig.10 Temporal variation of SWC at different subsidence locations before and after mining

4 主要結論與建議

4.1 主要結論

(1)從2000—2016年神東中心礦區(qū)植被NDVI時序上整體呈物候性周期變化，2005年后中心礦區(qū)植被INDV明顯提升，在2014年的INDV峰值最大。INDV差異值時序變化與該工作面煤炭開采年份密切相關。

(2)植被覆蓋度呈升高與降低的區(qū)域面積分別占中心礦區(qū)總面積的72.35%和27.65%。中心礦區(qū)年際間植被覆蓋度標準差處于0.04～0.46，植被覆蓋度波動以中、低幅度變化為主。

(3)土壤含水量是影響植被生長的第1重要土壤理化性質因子，影響植被生長較為重要的因子包括有機質含量、容重和全氮含量等。地下水埋深3.2 m和8.0 m為影響神東礦區(qū)植被INDV的2個重要閾值，地下水埋深通過影響土壤含水量影響地表植被類型，當地下水埋深大于3.2 m后，根系較淺的濕生植被演替為旱生植被；當地下水埋深大于8.0 m后，旱生植被演替為沙生植被。降雨量是中心礦區(qū)植被覆蓋度最敏感氣候因子。

4.2 神東礦區(qū)植被引導型恢復建議

半干旱礦區(qū)植被引導型恢復首要任務需要對工作面不同沉陷區(qū)位進行劃分，在此基礎上重點對拉伸區(qū)地表永久裂縫進行識別和治理，這是采后沉陷區(qū)土壤水分條件恢復重要基礎與前提。需要加強對沉陷區(qū)植物生長與土壤水分之間的耦合關系研究，分析植物生長受到干旱脅迫及植物死亡的相關土壤含水量閾值，可以為需要引導恢復時，土壤含水量提升到何種程度提供依據。

對于地下水位埋深在3.2 m以內的區(qū)域，植被恢復引導主要以降低植被對潛水的依賴為主，比如：種植抗旱性較強的本地植物、改善植被空間結構，但是這種工程類嘗試還比較少，其可行性還有待進一步研究論證。對于地下水位大于3.2 m的區(qū)域，筆者建議可以采用高空間分辨率的遙感影像對研究區(qū)地表植被進行精細分類，結合本文得到的植被INDV時序變化趨勢及波動程度，兼顧對植被生長立地條件的考慮(包括坡度、坡向、高程、土壤類型及其他理化性質)，以植被INDV呈正向波動區(qū)植被組成情況作為植被引導恢復時的植物配置、覆蓋度控制的參考依據，且植物配置時盡量選擇具有耐土壤貧瘠、抗旱性強、根系發(fā)達的本地植物物種。此外，應該加強對礦區(qū)地下水資源分布的探測，綜合考慮地下水的采前埋深和采后埋深，對采礦引起地下水變化導致的植被影響評價，加強對采空區(qū)積水的循環(huán)利用，有計劃地優(yōu)先開采無、少地下水的區(qū)域，從源頭減少或延緩對地下水的破壞。