神東礦區生態系統服務功能評價

劉 英，魏嘉莉，畢銀麗，岳 輝，何 雪

(1.西安科技大學 測繪科學與技術學院，陜西 西安 710054; 2.西安科技大學 地質與環境學院，陜西 西安 710054; 3.中國礦業大學(北京) 地球科學與測繪工程學院，北京 100083)

神東礦區是中國已探明儲量最大的煤田，然而該礦區生態環境極其脆弱，抗擾動能力差，煤炭開采引起的一系列生態問題如地表沉陷、植被退化以及水土流失等，進一步影響礦區生態系統服務功能，威脅礦區的可持續發展。面對大規模開采和生態保護之間的矛盾，采取科學的方法評估神東礦區的生態系統服務功能(Ecosystem Service,ES)，掌握其時空格局的演變規律及驅動因素，探索煤炭開采對其造成的影響，對兼顧神東礦區煤炭開采和生態環境的協調發展至關重要，可為礦區走資源保護型開采與生態環境治理相協調的綠色礦業之路提供科學支撐。

近年來，隨著人類生態保護意識的覺醒，生態系統服務日益被重視，一些生態系統服務評價模型應運而生，如ARIES(Artificial Intelligence for Ecosystem Services)、SoLVES(Social Values for Ecosystem Services)及InVEST(Integrate Valuation of Ecosystem Services and tradeoffs Tool)模型等，其中InVEST模型以其強大的空間表達功能和運行結果可視化功能，迅速成為生態系統服務功能定量評估和動態分析的有效手段，且該模型具有多個模塊，對應多種生態系統服務功能。國外學者利用InVEST模型對夏威夷、委內瑞拉、哥倫比亞以及亞洲印度尼西亞等區域[1-3]的ES進行了評估。國內學者集中于InVEST模型的水源供給、土壤保持、碳固定、水質的凈化及生物的多樣性等模塊，實現了多元化區域的生態系統服務功能評估，如橫斷山區[4]、韶關市[5]、白龍江流域[6]及黃河流域[7-8]等。目前國內對ES的研究日益成熟，大多將單項生態系統服務作為評價標準，缺乏對多種生態系統服務進行綜合評價的研究。

因此，筆者基于InVEST模型計算神東礦區1990—2018年的水源供給、土壤保持、碳儲存生態系統服務，對3者進行耦合建立綜合生態系統服務功能(Comprehensive Ecosystem Service,CES)評價模型，首先基于礦區尺度對神東礦區CES的空間分布特征、變化特征及聚類特征進行分析，并利用地理探測器探究CES變化的驅動因素;然后基于礦井尺度分析神東礦區不同程度開采區、礦井采區和非采區及復墾區CES的差異，探求煤炭開采對礦區CES的影響。本研究旨在從不同尺度綜合測度神東礦區CES的變化規律，為促進神東礦區綠色開采和生態可持續的協調發展提供科學依據。

1 研究區概況

神東礦區隸屬于神東煤炭基地，地處內蒙古鄂爾多斯和陜西北部交錯地帶，覆蓋黃土高原北緣與毛烏素沙漠過渡地帶東部(圖1)。礦區屬溫帶干旱、半干旱大陸性季風氣候，春季多風，夏季多暴雨，冬季長且寒冷;年降雨量為400 mm左右，年季差異較大，年均蒸散發達1 319 mm;平均風速為2.5～3.0 m/s，地勢西北高東南低。礦區以風沙和黃綿土為主，北部、東南部為黃土丘陵溝壑區，梁峁起伏，溝壑縱橫，地表支離破碎。地表原生植被種類單調，以耐旱、耐寒的沙生植物、旱生植物為主，呈稀疏灌叢景觀。礦區煤層賦存穩定，構造簡單，適宜機械化開采，主要以井工礦為主，多采用走向長壁全部垮落法開采，是全國最大的煤炭生產基地。

2 研究方法及數據

2.1 研究方法

2.1.1綜合生態系統服務功能評價模型的構建

(1)水源供給功能。InVEST模型中的產水量模塊基于水量平衡原理，利用降水量減去實際蒸散量計算每個柵格的水源供給量。計算公式為

(1)

式中，Yxj為第j類土地利用類型柵格x的年產水量，mm;Txj為第j類土地利用類型柵格x的實際年均蒸散發量，mm;Px為柵格x的年均降水量，mm。

(2)土壤保持功能。InVEST模型中的土壤保持模塊即泥沙輸移比模塊以通用土壤流失方程USLE為基礎，產生基于柵格單元的土壤保持量。本研究利用InVEST模型中土壤保持模塊計算礦區水蝕區的土壤保持量，計算公式為

Sx-Retention=Sx-Potential-Sx+Dx

(2)

Sx-Potential=RxKxLx

(3)

Sx=RxKxLxPxCx

(4)

(5)

式中，Sx-Retention,Dx為柵格x的土壤保持量、泥沙持留量，t；Sx-Potential為僅考慮地貌和氣候條件下柵格x的潛在土壤侵蝕量，t；Sx和Sy為考慮了管理、工程措施后柵格x及其上坡柵格y的實際侵蝕量，t；Rx，Kx，Lx，Cx，Px分別為柵格x的降雨侵蝕力因子、土壤可蝕性因子、坡長坡度因子、植被覆蓋和管理因子和土壤保持措施因子，這些因子的計算參考InVEST用戶指導手冊[9]；rx為柵格x的泥沙持留效率；ry為上坡柵格y的泥沙持留效率。

神東礦區遭受風蝕的面積約占全區的70%，因而利用土壤通用風蝕方程(Wind Erosion Equation，WEQ)監測神東礦區風力侵蝕區的土壤保持量。該模型通過土壤風蝕速率與各因子之間的經驗關系發展而來，被廣泛的應用于風蝕區的研究中。計算公式為

Ex=WxKxLxBx

(6)

Ex-Potential=WxKxLx

(7)

Ex-Retention=Ex-Potential-Ex

(8)

式中，E為柵格x的風蝕侵力模數;Ex為柵格x的實際風力侵蝕模數;Ex-Potential為柵格x的潛在風力侵蝕模數;Ex-Retention為柵格x的風力侵蝕土壤保持量;Bx為植被覆蓋因子，即研究區的植被覆蓋度;Wx為風蝕氣候因子，計算參考文獻[10-11]。

(3)碳儲存功能。碳儲量模塊基于土地利用類型和其對應的地上生物量碳庫(陸地表層上所有存活的植物材料的碳)、地下生物量碳庫(植物活的根系統中的碳)、土壤碳庫(土壤中的有機碳)和死亡有機碳庫(凋落物、倒立或站立著的已死亡的樹木的碳)，將這4種碳庫的碳儲存量相加來評價每個地圖單元和整個景觀的總碳儲存，公式為

Ci=Ci,above+Ci,below+Ci,soil+Ci,dead

(9)

(10)

其中，Ci為土地利用類型i的土壤及生物量總碳密度,t/hm2;Ci,above,Ci,below,Ci,soil,Ci,dead分別為土地利用類型i的地上生物量碳密度、地下生物量碳密度、土壤有機質碳密度和枯落物有機質碳密度;Ctotal為總碳儲量，t;Si為土地利用類型i的面積，hm2;n為土地利用類型數。考慮到四大基本碳庫已滿足研究目的及數據的難獲取性，模型中的采伐木材產品或相關木材產品斑塊的第五大碳庫不予考慮[10]。查閱相關資料和文獻[12-17]并采用生物量碳密度和土壤碳密度與降水量、氣溫的關系校正模型[18-20]得到研究區各土地利用類型的碳密度(表1)。

表1 不同土地利用類型的碳密度Table 1 Carbon density of different land use types

(4)構建綜合生態系統服務功能評價模型。考慮到神東礦區的水源供給、土壤保持、碳儲存生態系統服務功能具有相同的重要性，賦予其相同的權重以構建綜合生態系統服務功能評估模型SCE。計算公式為

(11)

式中，SEw，SEs，SEc分別為標準化后研究區水源供給、土壤保持、碳儲存生態系統服務量，取值均為[0，1];SCE為綜合生態系統服務量，取值為[0，1]，值越大表示綜合生態系統服務功能越高。

2.1.2梯度結構相似度

基于梯度的結構相似度模型(Gradient-based Structural Similarity,GSSIM)[21]從圖像組成的3個要素來解釋圖像的結構信息，公式為

GSSIM(x,y)=[l(x,y)]a[c(s,y)]β[g(x,y)]γ

(12)

其中，l(x，y)，c(x，y)，g(x，y)分別為2期圖像x，y的亮度、對比度、梯度結構度比較函數,l(x，y)，c(x，y)分別由x，y的均值、標準差計算得到，g(x，y)通過梯度幅值來計算，梯度幅值則采用Sobel算子計算的梯度得到;常數α，β和γ用于調整各分量的重要程度，取α=β=γ=1。GSSIM值越高，表示圖像x，y越相似，2者的變化越小。將GSSIM值分為3個區間:[0,0.25]為突變區，表示2期圖像x，y發生了顯著變化;(0.25,0.65]為中變區，說明變化處于中等水平;大于0.65的區域定義為低變區，說明兩期圖像變化并不明顯[22]。

2.1.3局部莫蘭指數

探究生態系統服務的空間聚類特征可有效把控其強弱分布的走向。空間聚類的高-高聚類要求要素具有高值且被其他具有高值的要素所包圍，低-低聚類則要求要素具有低值且被其他具有低值的要素所包圍。局部莫蘭指數(local Moran’s I)可將空間關聯結構模式可視化以衡量局域空間自相關性，對ES的空間聚集性進行測度[23]，其中空間聯系局域指標ILSA(Local Indicators of Spatial Assoction)的公式為

(13)

2.1.4地理探測器

地理探測器[24]中的因子探測器可探測自變量對因變量空間分異性的解釋程度，用q值度量，q值越大表示自變量對因變量的解釋力越強，反之則越弱;交互作用探測器可以識別不同影響因子的交互作用，評估2個自變量協同作用對因變量的解釋力。本研究利用地理探測器探究神東礦區生態系統服務空間分異的驅動因子。基于自然和人為因素2個方面，選取降水、坡度、高程、植被覆蓋度及土地利用類型5個影響因子。地理探測器要求輸入變量為類別數據，參考王勁峰提出的數據離散化方法并根據相關文獻的先驗知識[6,24-25]，將坡度按照<5°，5°～10°，10°～15°，15°～20°，20°～25°，25°～30°，30°～35°，>35°分為8級，植被覆蓋度按照<0.3，0.3～0.4，0.4～0.5，0.5～0.6，0.6～0.7，0.7～0.8，0.8～0.9，0.9～1.0分為8級，將高程、多年平均降雨量按等距分類法分為9類，土地利用類型依據類別進行分類。王勁峰指出提取柵格格網點密度越大計算精度越高，考慮到軟件可容納的最大行為32 767，因此基于ArcGIS隨機提取神東礦區30 000個點，作為地理探測器的運行數據。

2.2 數據來源

InVEST模型與WEQ模型的運行以多源數據為依托，本研究的數據來源與說明見表2。縱觀神東礦區的開采歷程，經歷了4個階段：① 1985—1990年，處于開發起步階段;② 1991—2000年，礦區確立了“高起點、高技術、高質量、高效率、高效益”的礦井建造方針，處于探索提高階段;③ 2001—2008年，礦區對小煤礦進行整合，進入集成創新開發階段;④ 2009—2018年，礦區進入快速發展跨越階段，高產、高效、高清潔的開采方式成為主流。因此，選取1990,2000,2005,2010,2015,2018年探究神東礦區各開采周期的綜合生態系統服務功能的演變規律。另外，本研究用到的遙感影像有1990年的Landsat 5和2018年的Landsat 8影像，對其進行輻射定標、大氣校正等預處理。

表2 數據來源與說明Table 2 Data source and description

3 結果與分析

3.1 基于礦區尺度的生態系統服務功能評價

3.1.1神東礦區CES時空分布特征

圖2為神東礦區各年份的SCE，神東礦區SCE由0.448 6(1990年)下降到0.382 5(2000年)后持續增加至0.471 6(2015年)，2018年又下降到0.453 2，整體以每年0.009的速度呈增加趨勢。1990年為神東礦區未開采期，自然生態系統受人為干擾較少，生態系統服務功能為較高水平;開采時期2000年生態系統服務功能減弱，但又在2005—2015年不斷增強，且2010年超出1990年水平，這與神東礦區實行煤炭開采與礦區生態環境治理工作并重的方針政策密切相關，1990年雖為未開采時期，但是其植被覆蓋極其稀少，神東礦區開采時期實行植被修復，植被覆蓋率由開發初的3%～11%提高到60%以上。岳輝等[34]的研究也表明神東礦區生態環境質量在2005—2015年總體趨于穩定并向生態環境質量好的方向發展。

圖2 神東礦區1990—2018年SCEFig.2 SCE of Shendong mining area from 1990 to 2018

由圖3可知，神東礦區CES在空間上具有一致性，均以西北部最弱，整體呈由西北向東南逐漸增強的格局。1990年礦區CES在西北部高低值混合分布，東南部CES則以高值為主;2000年的CES較1990年有所減弱，尤以西北部減弱明顯，CES低值區連片性有所增強，東南部CES也呈現退化趨勢;2005，2010，2015年礦區CES持續好轉，2005,2010年西北部礦井密集處的CES輕微增強，2010年礦區西南部CES有所增強，2015年礦區西北部的CES變化并不明顯，集中在東南部呈現增強態勢;2018年CES全區較2015年有所減弱，與1990年相比，礦區西北部CES呈減弱趨勢，南部CES得到增強。整體來看，神東礦區CES往好的方向發展。神東礦區CES呈現西北部最弱，且由西北向東南逐漸增強的格局，可能是由于神東礦區西北部自1990年后工礦用地不斷擴張，加之一部分未利用地穩定分布于此，兩者的各項CES均較弱，導致西北部CES處于低水平，而東南部有部分草地和耕地覆蓋，可在一定程度上抑制土壤流失、增強碳匯;另外，礦區降雨自西北向東南逐漸增多，水源供給受降雨分布的影響較為明顯，且東南區域土壤侵蝕以水蝕為主，高降雨量區域的降雨侵蝕力較高會導致土壤潛在侵蝕大，土壤保持量因此表現為高值，因此水源供給和土壤保持均表現為西北弱東南強的特征。

圖3 神東礦區1990—2018年生態系統服務空間分布Fig.3 Spatial distribution of ecosystem services in Shendong mining area from 1990 to 2018

3.1.2神東礦區CES空間變化特征

為定量分析神東礦區CES的空間變化特征，基于MATLAB平臺利用GSSIM模型計算神東礦區1990—2000，2000—2005，2005—2010，2010—2015，2015—2018，1990—2018年6個時期的GSSIM并生成空間分布圖(圖4)。統計各變化區所占比例并結合圖4可知，神東礦區CES空間變化在1990—2000年以低變區為主(占比為45.69%)，突變區占比16.12%，在礦區西北部分布連片，該區域在1990—2000年之間建成投產了寸草塔、大柳塔等礦井，這可能是導致此區域CES發生突變的原因;神東礦區CES空間變化在2000—2005，2005—2010年以低變區為主(占比為64.14%，51.37%)，在全區均有分布，中變區與突變區呈細碎狀分布在礦區西部以及礦井密集的西北部，這是由于2000—2010年，礦區大部分的植被有所改善，植被退化的區域主要集中在礦區的西部和西北部地區;2010—2015年礦區東南部突變區明顯增多，礦區CES發生明顯變化的區域由西北部轉移到東南部，這是由于2010年之后，礦區植被覆蓋逐漸呈自西向東逐漸增加的規律，加之2013年礦區內降雨量驟然增加(降雨量達661 mm，高出平均51.2%)，東南部的植被狀況長勢越發良好，水源供給與碳儲量均有所增強，與2010年相比形成突變;2015—2018年全區以低變區為主，突變有所減少仍分布在礦區東南部;1990—2018年仍以低變區為主(40.04%)，全區均有分布，突變區次之(30.72%)，呈細條斑塊狀分布，西北部較其余區域相對密集。

圖4 1990—2018年各時段GSSIM空間分布Fig.4 Spatial distribution of GSSIM in various periods from 1990 to 2018

為進一步分析神東礦區突變區和中變區CES的變化狀態，以1990—2018年GSSIM影像為例在圖上任意選取突變區3個樣點區(A，B，C)與中變區3個樣點區(E，F，G);并利用1990，2018年Landsat影像和1990—2018年的生態系統服務影像差值圖進行驗證(圖5)。由圖5可知，1990—2018年，突變樣點區A的植被覆蓋得到了明顯改善，結合1990和2018年的CES差值圖可發現，此樣點區的CES有所增強;突變樣點區B由低植被覆蓋區變成了工礦用地，差值圖顯示此樣點區的CES減弱;突變樣點區C由低值被覆蓋區變為露天開采區的排土場，差值圖顯示其CES有所減弱;中變樣點區D，E的地表覆蓋發生了明顯改變，均由低植被覆蓋區變為露天開采區或建筑區，差值圖顯示中變樣點區D，E的CES均呈減弱趨勢;中變樣點區F的植被覆蓋有所增加，CES為增加狀態。由此可知，神東礦區CES減弱的區域往往伴隨著工礦用地擴張對植被覆蓋區的侵占，而神東礦區CES增強的區域植被覆蓋度呈增加趨勢，因此神東礦區CES的空間變化與植被覆蓋和用地類型的變化密切相關。

圖5 突變樣點區(A，B，C)和中變樣點區(D，E，F)Fig.5 Mutation sample point area (A,B,C) and intermediate variation sample point area (D,E,F)

為定量分析突變區、中變區及低變區CES的變化規律，將各時間段CES變化差值圖與對應時段突變區、中變區及低變區的矢量疊加，統計不同變化區生態系統服務量增加和減少區域所占面積百分比(表3)。由表3可知，1990—2000年CES的下降主要來源于低變區(占比為32.71%)，中變區次之(28.64%);2000—2005年、2005—2010年，盡管突變區和中變區增加區域面積在持續增加，但是CES的增加均以低變區為主(占比分別為53.14%，42.43%);2010—2015年，突變區、中變區和低變區CES增加區域比例依次為19.19%，21.15%，20.85%，此階段突變和中變增加的區域對SCE增加的貢獻增大;2015—2018年，低變區再次成為CES的下降主要貢獻者。1990—2018年，突變區和中變區增加和減少的區域并無明顯差異，低變區增加的區域大于減少的區域，且低變區占據全區的40.04%，說明1990—2018年CES的變化仍以低變區變化為主。總體上，神東礦區CES的變化主要是由于礦區內低變區的變化導致，中變區和突變區對神東礦區CES變化的貢獻力次之。

表3 生態系統服務量增加和減少區域所占面積百分比Table 3 Increase and decrease percentage of CES in study area

3.1.3神東礦區CES空間聚類分析

本研究基于GeoDa空間分析軟件，利用局部莫蘭指數對神東礦區CES的空間聚類進行測度，以有效把控礦區CES強弱分布走向。由圖6可知，神東礦區CES的空間聚類特征以不顯著為主，低-低聚集次之，高-低聚集最少。不顯著區域分布在礦區中部及東南部區域，低-低聚集區由1990年的西北部、南部及河流兩岸破碎分布轉變為2000—2018年的西北部連片分布，在1990年礦區開采初期，西北部CES持續處于最弱水平;高-高聚集區的面積較小，零落分布在河流兩岸和礦區南部的部分區域，由土地利用類型圖可知河流沿岸分布著CES較強的耕地，所以神東礦區河流沿岸的CES水平較高;1990年低-高聚集區的分布處于礦區西北部，2000年以后主要分布在礦區東南部，該區域1990年被草地和耕地覆蓋，2000年以后被局部的工礦用地和建設用地破碎分割，因此呈低值區被高值區包圍的格局。整體上看，神東礦區的空間聚類特征明顯，呈“西北部最弱，河流沿岸最強，東南部不顯著”的格局。此格局可能是氣候、植被、人類活動等綜合因素的影響所致，礦區降雨自西北向東南逐漸增多，西北部因采礦造成的植被覆蓋破壞均為西北部低-低聚類區提供了動力。因此，應將礦區西北部作為實行“煤炭開采與礦區生態環境治理工作并重”政策的傾斜區，同時對河流沿岸的高-高聚類區實施保護政策，穩定此格局，并采用提高植被覆蓋、及時復墾廢棄礦區及改進綠色開采技術等策略，提高神東礦區的生態系統服務功能。

3.1.4神東礦區CES變化影響因素分析

由地理探測器結果可知，對神東礦區CES解釋力最大的因子為土地利用(q=0.173)，其次為坡度(q=0.144)，然后依次為高程(q=0.106)，降雨量(q=0.102)、植被覆蓋度(0.057)，且各影響因子均通過了置信水平99%的顯著性檢驗。土地利用的解釋力最高是由于不同地類的固碳能力存在差異，其蒸散能力與泥沙持留效果的不同導致水源供給和土壤保持功能也存在差異;坡度和高程對生態系統服務功能的影響主要體現在坡度小高程低的區域人為活動較為頻繁，對生態系統的擾動較大;降雨量對水源供給的作用較為明顯，土壤保持和碳儲量受到間接影響，所以解釋能力次于坡度與高程;植被覆蓋度高的區域碳匯較強，可抵御降雨侵蝕及其他物理性侵蝕，產水能力高而蒸散系數低，對生態系統功能均有正向促進作用，然而盡管神東礦區的植被覆蓋在逐漸改善，但沙地植被占絕對優勢，弱化了對生態系統服務的促進作用。因此，神東礦區應著重進行以提高植被覆蓋為導向的土地利用結構優化方式，人工調控植被演替方向與速度，優化草本為主、草灌結合的林分結構，采用微生物復墾等科學方法復墾礦區沉陷地，對未利用地加大綠化力度等。

交互探測器的統計結果(表4)表明，神東礦區生態系統服務量的影響因子交互作用均高于單因子。土地利用與其他因子交互作用居于首位，且與高程的交互作用最為明顯(q=0.276);植被覆蓋度與其余因子交互作用提升明顯，與高程、坡度和土地利用的交互均為植被覆蓋度單因子的3倍以上，與土地利用的交互作用提升最多;高程與坡度的交互作用為單因子的1.5～2倍;降雨量與其他因子的交互作用提升并不明顯，與土地利用的交互作用相對較高。分析可知，神東礦區的生態治理應以優化土地利用格局為先導，根據礦區高程的分布特征實施草地修復、沉陷地復墾等政策;神東礦區在坡度30°以下的區域隨坡度的遞加CES不斷增強，30°以上發生減弱，有研究指出[6]隨著坡度的遞增，土壤保持功能整體呈先增強后減弱的趨勢，因此神東礦區應對坡度與植被覆蓋、人類活動進行動態監測，掌握CES強度分布的大致走向，將工作重點放在植被覆蓋度較低的陡坡(30°以上)區域，也可在緩坡區域(15°以下)實施生態管理措施以平衡采礦及其他人類活動的地表擾動，以期提高神東礦區生態系統服務能力。

表4 各項生態系統服務影響因子交互作用q值Table 4 q values of interaction of impact factors

3.2 礦井尺度CES差異分析

3.2.1不同開采強度下的CES差異分析

神東礦區煤層構造簡單，主要以井工礦開采為主，是我國最大的井工煤礦開采基地，為了實現分區探究礦區的生態系統服務功能，通過不同區域的開采強度進行進一步分析。范立民等[35]指出煤礦的開采強度可由平面開采強度和立面開采強度決定，即開采強度主要取決于工作面長度和采厚大小2個因素。本研究嘗試探究礦區內不同開采強度下CES的差異及變化規律。神東礦區開采強度的劃分依據文獻[35-36]的研究成果，將神東礦區劃分為開采強度極高、強度高、開采強度中等和低強度4級，得到神東礦區開采強度分區圖(圖7)。極高強度開采區位于烏蘭木倫河和窟野河之間，此處分布著大柳塔、活雞兔等千萬噸礦井;高強度采區的特點是工作面較為連續，采煤高度大(1.3～4.5 m)，單位面積采出煤炭量高;中低強度采區的煤炭作業方式并未全部實現機械化開采，少部分為整合煤礦，產能低，工作面較為破碎。

圖7 神東礦區不同開采強度分布Fig.7 Distribution map of different mining intensity in Shendong mining area

將不同程度開采區的矢量與礦區生態系統服務功能圖疊加，對不同程度開采區的SCE進行統計(圖8)。從各程度開采區的SCE差異來看，6個年份SCE大小順序:中強度開采區>低強度開采區>未開采區>高強度開采區>極高強度開采區，中強度開采區的綜合生態系統服務最強，多年SCE平均值達到0.502 8;極高強度開采區的綜合生態系統服務最弱，SCE平均值為0.430 8。各個程度開采區的SCE均在1990—2000年和2015—2018年呈下降趨勢，1990—2000年未開采區SCE下降量為各區下降總量的22.34%，該值在2015—2018年減少為15.09%，而高強度采區和極高強度采區CES下降量在1990—2000年為各區下降總量的7.52%，18.19%，2015—2018年上升為21.67%，26.94%;2015—2018年與1990—2000年相比，導致礦區SCE減少的未開采區的貢獻率變小，而高強度、極高強度的貢獻率變大了，說明神東礦區在各個開采區實行的生態保護措施并未完全抵消掉采礦活動帶來的負面影響。而且，在2010—2015年，極高強度的SCE發生下降，其他程度采區則處于上升趨勢，說明極高強度開采區對SCE的增長具有抑制作用。張文凱等[36]的研究也表明，神東礦區中強度開采區的NPP最高，極強度開采區的NPP最低，與本研究結果一致。

圖8 不同開采強度區SCEFig.8 SCE quantities in different mining intensity zones

3.2.2采區和非采區的CES差異分析

為進一步探究采礦活動對生態系統服務的影響，選取神東礦區的7個主要礦井對其采區和非采區的生態系統服務量及變化比進行統計(表5)。從采區與非采區的差異來看，補連塔礦井的CES均值由采區大于非采區，其變化比由2005年的34.48%大幅下降為2010年的0.92%之后保持低值，說明2010年之后補連塔礦井采區和非采區的差距大幅縮小;大柳塔礦井采區均值均小于非采區，變化比介于-12.32%～-5.30%，說明采區與非采區差異較為穩定;活雞兔礦井采區均值大于非采區，2015年之后小于非采區，且2015年之后采區非采區差異增大;哈拉溝礦井采區均值小于非采區，變化比持續下降至2015年的-10.67%后上升，說明哈拉溝礦井采區、非采區的差異先變小后變大;烏蘭木倫礦井采區均值小于非采區，采區非采區差異十分顯著(變化比為-35.41%～-26.88%)，且該礦井的CES最弱，可能是由于該礦井位于干旱少雨、植被覆蓋低的毛烏素沙漠邊緣;榆家梁礦井由最初的采區大于非采區在2010年后逆轉為采區小于非采區，變化比介于-3.85%～1.96%;石圪臺礦井采區均值小于非采區，變化比由2010年的-14.03%上升為-5.49%，說明非采區和采區的差異在變小。

表5 各礦井采區與非采區生態系統服務量及變化比Table 5 Ecosystem service quantity and change ratio between mining area and non-mining area

采區與非采區生態系統服務量的絕對差異僅能表明2者的CES是否存在明顯差異，差異的來源可能是采區與非采區生態環境的差異性或礦區生態環境治理的影響。因此為深入分析地下采礦活動的影響，需進一步對比兩者的相對差異。神東礦區大規模建設開始于20世紀90年代，以1990年為采礦前的初始狀態(此時采礦活動并未進行，采區、非采區未受采礦活動的影響)，以此時礦井內采區與非采區初始差異變化比作為參考基準，采后狀態選取2018年，若采礦后的差異變化比小于參考基準，說明礦井采區受到了地下開采活動的負面影響，相反，則可能受到了生態保護措施的正面影響[37-38]。

將神東礦區的7個主要礦井采區與非采區范圍與1990年、2018年的生態系統服務量疊加即可得到采區與非采區的初始差異和采后差異變化。由表6可知，① 補連塔礦井各個年份的采后變化比均高于初始變化比(差值為0.51%～7.74%)，說明采礦并未對采區CES產生顯著的負面影響，可能受到了生態治理的正面影響;② 大柳塔礦井采后變化比在2010年之前大于采前變化比，2010年之后小于采前變化比，因此該礦井在2010年以后受到了采礦活動的負面影響;③ 活雞兔礦井采后變化比在2010年之前大于采前變化比，2010年之后小于采前變化比，差值由1.54%(2005年)下降為-15.60%(2010)年，該礦井在2010年以后受到了采礦活動的負面影響;④ 烏蘭木倫礦井的SCE最低，采前采區SCE大于非采區，采后采區SCE小于非采區，且初始與采后狀態的差異變化最顯著(采后變化比較采前下降23.10%～36.28%)，說明采區與非采區差異在采后大幅增加，采礦活動對采區CES的負面影響較為顯著;⑤ 哈拉溝礦井采前采區SCE大于非采區，采后采區SCE小于非采區，且采后變化比小于采前變化比(差值為-10.98%～-8.90%)，說明哈拉溝礦井的CES受到采礦活動的負面影響;⑥ 榆家梁礦井采后變化比小于采前變化比(差值為-4.08%～-0.57%)，且該礦井未采取植被改善措施，說明地下采礦活動對采區CES產生了負面影響;⑦ 石圪臺礦井采區CES低于非采區，且采后變化比小于采前變化比(差值為-6.51%～-4.58%)，說明石圪臺采區受到了采礦活動的負面影響。總體上，補連塔礦井采區的CES受到了生態治理措施的正面影響，大柳塔、活雞兔礦井采區的CES在2010年以后受到了采礦的負面影響，烏蘭木倫礦井采區CES受到采礦的顯著負面影響，石圪臺、哈拉溝及榆家梁礦井采區的CES自采礦以來持續受到采礦的負面影響，或采礦的負面影響高于生態保護措施的正面影響。

表6 各礦井采區和非采區相對變化比Table 6 Relative change ratio of mining area and non-mining area in each mine

為緩解礦區干旱、貧瘠、地表沉降植被受損等生態問題，神東礦區于2008年開始和中國礦業大學(北京)合作開展微生物復墾關鍵技術研究與試驗，主要是通過從當地土壤中篩選出適宜的叢枝菌根真菌，培養后接種于植物根系，擴大植物根系對土壤水分和養分吸收能力與吸收量，利用菌絲修復了斷根，以提高植物成活率和植株生長量。

本研究為探究礦區微生物復墾工作對生態系統服務功能的影響，選取大柳塔礦井采煤沉陷地微生物復墾試驗區，探究其CES對微生物復墾的響應。將大柳塔采區分為復墾試驗區和其他區域，統計SCE的值和變化比(表7)。由表7可知，2000—2005年復墾試驗區SCE與采區其他區域的變化比呈下降趨勢(變化比分別為2.64%，2.03%)，2010，2015，2018年變化比呈持續增加態勢(變化比分別為5.10%，6.36%，9.45%)，說明復墾實驗區的SCE與采區其他區域的差距越來越大，生態系統服務功能往變強的方向發展。礦區微生物修復是一個生態持續演變過程[39]，在2008年之后礦區微生物復墾對CES持續起正向促進作用。畢銀麗等[39]指出微生物復墾后，土壤質地發生改變、土壤保水能力增強，生態效應也得到了改善。

表7 大柳塔采區與復墾試驗區的變化比Table 7 Change ratio between Daliuta mining area and reclamation experimental area

4 討 論

鑒于已有生態系統服務的文獻缺乏對生態系統服務的綜合評價，本研究通過耦合單項生態系統服務綜合測度神東礦區的生態系統服務量，并從礦區和礦井尺度考量了研究區生態系統服務的變化規律，彌補了對神東礦區生態系統服務多尺度研究的空缺。

大量學者對神東礦區及主要礦井的生態環境質量、植被覆蓋度及土壤濕度進行了研究。筆者[40]指出神東礦區植被狀況2000—2015年有改善趨勢且未來會有較弱的退化趨勢，本研究中神東礦區生態系統服務在2000—2015年逐漸增強，2018年發生減弱，與前人研究結果一致并相互驗證。基于礦井尺度，筆者[37]指出2010年以前，采礦活動并未對大柳塔礦井采區的植被產生影響，這與本研究得到的大柳塔礦井2010年之后逐漸受到采礦負面影響的結論一致;岳輝等[34]指出烏蘭木倫采區生態環境質量最差，本研究烏蘭木倫CES也處于最低水平;雷少剛[38]指出榆家梁礦井在2005年前受到采礦活動的負面影響，本研究結果也表明地下采礦活動對榆家梁采區CES產生了負面影響。本研究更新了基于神東礦區礦井尺度的研究年限，進一步厘清了主要礦井是否受采礦活動的干擾，客觀評價了礦井生態保護措施的實施成果。筆者等[22]指出礦區土壤濕度變化與地表覆蓋高度相關且與采礦活動有關，這與本文得到的神東礦區CES與植被覆蓋、土地利用類型變化密切相關保持一致。

本文構建了神東礦區綜合生態系統服務評價模型，通過前人的研究進行間接驗證，證明了該模型的綜合評估結果能夠定量反映神東礦區生態系統服務的演變規律。建議神東礦區將西北部作為實行“煤炭開采與礦區生態環境治理工作并重”政策的傾斜區，并著重進行以提高植被覆蓋為導向的土地利用結構優化，利用微生物復墾等科學方法復墾礦區廢棄土地等。研究成果可為綜合評價區域生態系統服務拋磚引玉，而區域的多項生態系統服務的耦合機制仍然是一個具有挑戰性的難題，應進一步強化決策過程與機理耦合，依托多源數據獲取生態系統服務功能供給、需求、目標管理等參數信息，致力于多項生態系統服務多尺度優化集成模型的開發與應用。因而本研究下一步將選取多種生態系統服務功能，尋求多種耦合機制，以求得更加精確的實驗結果，促使生態系統服務評估為區域生態系統服務可持續管理提供更精準可靠的決策支持。

CES受自然、社會和經濟多種復雜因素的影響，探究生態系統服務的驅動因素時，本研究所選因子多基于自然因素，進一步研究可融入人口密度、經濟指標、煤炭產量等，以探究自然社會和經濟因素對神東礦區生態系統功能變化的綜合影響。另外，基于生態系統服務功能與人類福祉的緊密聯系性，也可探尋生態系統服務與社會經濟系統的耦合關系，為區域提升生態系統服務功能、實現可持續發展提供生態決策。

5 結 論

(1)礦區尺度上，① 神東礦區SCE由1990年的0.448 6下降為2000年的0.382 5，持續增加至2015年的0.471 6又下降為2018年的0.453 2，整體呈增加趨勢，速率為0.009 a-1;空間上神東礦區西北部CES最弱，整體呈由西北向東南部逐漸增強的格局;② 神東礦區CES空間變化特征明顯，礦區CES變化主要由礦區內低變區的變化導致，中變區和突變區次之，神東礦區CES空間變化與植被覆蓋、用地類型的變化密切相關;③ 神東礦區CES空間聚類表現為“西北部最弱、河流沿岸最強、東南部不顯著”的格局;④ 神東礦區CES變化的影響因子解釋力大小為:土地利用類型>坡度>高程>降雨>植被覆蓋度，且因子兩兩交互的驅動力較單因子提升明顯。

(2)礦井尺度上，① 神東礦區不同程度開采區CES大小依次為中強度開采區、低強度開采區、未開采區、高強度開采區、極高強度開采區;② 補連塔礦井采區的CES受生態治理措施的正面影響，使得采區CES大于非采區;大柳塔、活雞兔礦井采區CES在2010年以后受到了采礦的負面影響;烏蘭木倫、石圪臺、哈拉溝及榆家梁礦井采區的CES自采礦以來持續受到采礦的負面影響，或采礦帶來的負面影響高于生態保護措施帶來的正面影響;③ 大柳塔復墾實驗區CES往變強的方向發展。