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

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

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

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

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

1 研究區(qū)概況

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

2 研究方法及數(shù)據(jù)

2.1 研究方法

2.1.1綜合生態(tài)系統(tǒng)服務功能評價模型的構建

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

(1)

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

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

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的泥沙持留效率。

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

Ex=WxKxLxBx

(6)

Ex-Potential=WxKxLx

(7)

Ex-Retention=Ex-Potential-Ex

(8)

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

(3)碳儲存功能。碳儲量模塊基于土地利用類型和其對應的地上生物量碳庫(陸地表層上所有存活的植物材料的碳)、地下生物量碳庫(植物活的根系統(tǒng)中的碳)、土壤碳庫(土壤中的有機碳)和死亡有機碳庫(凋落物、倒立或站立著的已死亡的樹木的碳)，將這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為土地利用類型數(shù)。考慮到四大基本碳庫已滿足研究目的及數(shù)據(jù)的難獲取性，模型中的采伐木材產(chǎn)品或相關木材產(chǎn)品斑塊的第五大碳庫不予考慮[10]。查閱相關資料和文獻[12-17]并采用生物量碳密度和土壤碳密度與降水量、氣溫的關系校正模型[18-20]得到研究區(qū)各土地利用類型的碳密度(表1)。

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

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

(11)

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

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

2.1.3局部莫蘭指數(shù)

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

(13)

2.1.4地理探測器

地理探測器[24]中的因子探測器可探測自變量對因變量空間分異性的解釋程度，用q值度量，q值越大表示自變量對因變量的解釋力越強，反之則越弱;交互作用探測器可以識別不同影響因子的交互作用，評估2個自變量協(xié)同作用對因變量的解釋力。本研究利用地理探測器探究神東礦區(qū)生態(tài)系統(tǒng)服務空間分異的驅動因子?；谧匀缓腿藶橐蛩?個方面，選取降水、坡度、高程、植被覆蓋度及土地利用類型5個影響因子。地理探測器要求輸入變量為類別數(shù)據(jù)，參考王勁峰提出的數(shù)據(jù)離散化方法并根據(jù)相關文獻的先驗知識[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類，土地利用類型依據(jù)類別進行分類。王勁峰指出提取柵格格網(wǎng)點密度越大計算精度越高，考慮到軟件可容納的最大行為32 767，因此基于ArcGIS隨機提取神東礦區(qū)30 000個點，作為地理探測器的運行數(shù)據(jù)。

2.2 數(shù)據(jù)來源

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

表2 數(shù)據(jù)來源與說明Table 2 Data source and description

3 結果與分析

3.1 基于礦區(qū)尺度的生態(tài)系統(tǒng)服務功能評價

3.1.1神東礦區(qū)CES時空分布特征

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

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

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

圖3 神東礦區(qū)1990—2018年生態(tài)系統(tǒng)服務空間分布Fig.3 Spatial distribution of ecosystem services in Shendong mining area from 1990 to 2018

3.1.2神東礦區(qū)CES空間變化特征

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

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

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

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

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

表3 生態(tài)系統(tǒng)服務量增加和減少區(qū)域所占面積百分比Table 3 Increase and decrease percentage of CES in study area

3.1.3神東礦區(qū)CES空間聚類分析

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

3.1.4神東礦區(qū)CES變化影響因素分析

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

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

表4 各項生態(tài)系統(tǒng)服務影響因子交互作用q值Table 4 q values of interaction of impact factors

3.2 礦井尺度CES差異分析

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

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

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

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

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

3.2.2采區(qū)和非采區(qū)的CES差異分析

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

表5 各礦井采區(qū)與非采區(qū)生態(tài)系統(tǒng)服務量及變化比Table 5 Ecosystem service quantity and change ratio between mining area and non-mining area

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

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

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

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

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

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

4 討 論

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

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

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

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

5 結 論

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

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