田惠文,張欣欣,畢如田,朱洪芬,席 茜
(1.山西農業大學 資源環境學院,山西 晉中 030801; 2.中國地質大學 地球物理與空間信息學院,湖北 武漢 430074)
陸地生態系統碳循環在全球碳收支中起著主導作用[1],而農田生態系統是構成陸地生態系統的三大系統之一,且農田生態系統碳庫是全球碳庫和陸地生態系統碳庫的重要組成部分,在全球碳平衡中扮演著重要角色[2]。農田生態系統受人類活動影響最為顯著,其中不同的管理措施、農田的棄耕是學者關注的重點,如石小霞等對華北平原不同農田管理措施與土壤碳庫的變化關系進行了研究[3];王海候等研究了不同秸稈還田模式與土壤碳庫特性之間的關系[4];趙威等對豫西丘陵坡地棄耕農田植被演替對土壤碳庫的影響開展了研究[5]。此外,水土流失、土地鹽堿化及土地利用變化等也會對農田碳庫產生擾動,且這些因素都不可避免地耦合著人類活動的影響。
工業革命以來,人類生產方式發生了翻天覆地的變化,空氣中CO2等氣體含量不斷上升,導致了溫室效應和一系列生態環境問題[6],其中煤炭開采是重要的原因之一。據估算,在1989—2010年由煤礦開采造成的土地損毀面積達100.01×104hm2[7],礦區生態環境面臨著嚴峻考驗。農田生態系統碳庫主要包括植被碳庫和土壤碳庫[8-9],近年來,有學者對煤礦區土壤碳匯進行了研究,如余建等研究了徐州市柳新鎮高潛水位礦區不同復墾及利用方式的土壤有機碳、無機碳的構成與變化特征[10];程靜霞等以焦作礦區為例,研究了煤炭開采沉陷區土壤有機碳的空間變化[11]。也有學者對礦區植被碳匯展開研究,如侯湖平等以徐州九里礦區為例,研究發現煤炭開采對區域農田植被碳儲量的影響屬于失碳效應[12];侯湖平等研究了煤礦區土地利用變化對生態系統植被碳儲量的影響,發現礦區碳總量降低[13];黃翌等以山西大同礦區為例,研究發現多年來礦區植被-土壤系統碳匯量減少[14];徐占軍等采用DNDC模型研究了煤田開采沉陷對農田土壤和植被碳庫的擾動[15]。煤礦區作為一個特殊的復雜地理區域,長時間、大規模和高強度的煤炭開采不可避免地破壞原生的礦床地質條件,引發植被退化、水土流失及物種減少等[16-17],尤其是地下開采造成了大面積塌陷地、土地沉陷和耕地破壞等生態環境問題,嚴重制約了農田生態系統碳庫的積累[18-20]。煤礦開采、運輸過程中也會產生大量的煤煙、粉塵以及有害氣體,對周邊的植被生長、土壤質量及居民生活造成負面效應[21]。目前,對于礦區農田生態系統的碳儲存研究主要集中在礦區沉陷、采空區、損毀土地等直接受損區域上的植被固碳或土壤碳儲存,同時考慮植被-土壤碳匯變化的研究較少且忽略了煤炭開采對直接受損區以外的周邊地區。
因此,完善煤炭開采導致的農田生態系統固碳損失的評估方法不僅對礦區土地復墾、生態修復有著極其重要的意義,還可為煤礦區生態補償測算提供決策與意見。筆者以山西省晉城市長河流域為例,對煤炭開采導致的農田生態系統固碳損失(植被固碳和土壤碳儲存)進行評估。本研究的主要目的:建立一種適用于煤炭開采導致農田生態系統固碳損失的方法;量化煤炭開采對周邊農田的間接受損范圍;估算礦區農田生態系統固碳損失,為礦區土地復墾與生態修復提供依據。
長河流域位于山西省晉城市澤州縣西北部,在112.63°~112.77°N,35.50°~35.63°E之間,下轄3個鄉鎮,土地利用類型主要有林地、草地、農田、采礦用地、建設用地、沉陷積水區等,總面積約113.16 km2,其中農田面積為64.55 km2(圖1)。研究區屬溫帶大陸性季風氣候,年平均氣溫10.6 ℃,年降雨量550~600 mm,地上農作物類型主要是玉米、冬小麥等,農田耕作方式為小麥-玉米輪作,施肥方式為傳統的土壤施肥,施肥量基本相同,約為N肥1.5 t/km2,K肥0.9 t/km2,P肥0.45 t km2。研究區海拔約857 m,地形以山地、丘陵為主,境內土壤類型為褐土和紅黏土,其中褐土面積占91.14%。長河流域位于山西沁水煤田界內,流域內共有天泰坤達等11個煤礦,煤炭資源豐富。2000—2015年因煤炭開采造成的采空區、沉陷區、損毀土地等面積約10.29 km2,造成了植被覆蓋度降低、水土流失和環境污染等一系列嚴重的生態問題。

圖1 研究區Fig.1 Study area
(1)礦區數據。長河流域各礦業提供的采空區、已損毀和已塌陷等土地現狀分布數據。
(2)土地利用數據。來源于2000年和2015年澤州縣土地利用規劃空間數據庫的1∶1萬土地覆蓋數據。
(3)農田數據。提取自澤州縣土地利用規劃空間數據庫的1∶1萬土地覆蓋數據。
(4)施肥制度。來源于澤州縣耕地質量監測評價數據,共16個樣點。
(5)耕作方式。提取自山西省測土配方施肥數據,共528個樣點。
(6)氣象數據。來源于中國氣象數據網,包括氣溫、降水量等,時間序列為2000—2015年。
(7)遙感數據。來源于美國地質勘探局(United States Geological Survey,USGS)網站的MODIS和Landsat影像,時間序列為2000—2015年。數字高程模型(DEM)由各礦區提供,精度為1 m,并采用DEM數據計算坡度、坡向。
(8)土壤類型數據。來源于山西省第2次土壤普查數據庫。
(9)地貌類型數據。來源于澤州縣國土資源局。
(10)植被數據。2013年實測的植被凈初級生產力(Net Primary Production,NPP)數據,實測方法見4.3。
(11)土壤數據。2015年7月在長河流域采樣,并在實驗室測試了土壤有機碳、土壤容重等數據。
本文研究對象為農田,所有數據均采用農田矢量數據進行了裁剪。以上數據均經空間數據格式轉換、校正并統一投影轉化為CGCS2000坐標系(圖1,2)。
首先利用礦區數據,提取2000年以來因煤炭開采造成的采空、損毀及塌陷范圍,作為礦區直接受損區;然后綜合考慮海拔、坡度、土壤類型、到礦區距離等因素,利用空間疊加分析,選擇出與礦區生境相似且未受煤炭開采干擾的區域作為參照區;基于Landsat影像計算歸一化差值植被指數(Normalized Difference Vegetation Index,NDVI),分析直接受損區邊緣向外的NDVI指數變化特征,確定間接受損區;模擬長河流域農田植被NPP,測定農田土壤有機碳,并與受損區、參照區范圍疊加,計算得受損區農田植被固碳損失量和土壤碳儲存損失量;結合農田的土地利用變化情況,估算礦區農田生態系統固碳損失量,技術路線如圖3所示。
固碳主要是指生態系統的碳凈吸收量,其主要是生態系統植被固碳和土壤碳儲存,即來自于植物光合作用固定的碳及枯枝落葉等進入土壤形成的有機碳[22]。從生態系統碳固定的角度來看,植被—土壤可以看作一個統一的系統,植物的光合作用是該系統的碳來源,植被與土壤中的呼吸作用為碳支出,剩余的轉化為土壤有機質(圖4)。在長河流域農田生態系統中,植被固碳量近似地以植被NPP表示(植物光合作用扣除自氧呼吸后的差值),土壤有機碳通過樣品采集,實驗室測試獲得。因此,只需要獲取煤炭開采以來農田受損區及參照區的植被NPP及土壤有機碳,就可以估算礦區農田生態系統的碳損失量。

圖2 數據源Fig.2 Data source

圖3 技術路線Fig.3 Technology flow chart

圖4 固碳機制Fig.4 Mechanism of carbon sequestration
土地利用變化(農田轉變為非農田)是對區域農田生態系統固碳影響最顯著的因素[23-24]。由于本研究劃定的受損區為農田,故農田轉化為其他土地利用類型的區域不歸入受損區,但估算礦區農田生態系統固碳損失時,仍考慮因煤炭開采導致的農田面積減少造成的固碳損失。在土地利用類型未發生改變的情況下,影響農田生態系統固碳量的常規影響因素包括:氣候、地貌類型、土壤類型、耕作方式、施肥制度以及農業科技水平等[13]。本研究假定土地利用類型未發生改變,充分考慮常規因素的影響來確定參照區和受損區,最后結合農田土地利用類型的轉化情況估算煤炭開采導致的農田生態系統固碳損失量。
從礦區各個礦業公司提供的礦區數據中,提取2000年以后由于煤炭開采造成的采空、塌陷、損毀等土地數據作為直接受損區,經計算,山地直接受損區面積為3.52 km2,丘陵直接受損區面積為1.54 km2(圖5)。

圖5 受損區和參照區Fig.5 Damaged area and reference area
固碳損失量估算的關鍵在于參照區的選擇。受損區位于山地和丘陵,由圖2可知,長河流域的耕作方式、作物種類、施肥制度及農業科技水平基本相同,但山地和丘陵農田的灌溉條件、管理措施不同,為了準確的評估與檢驗,本研究分別選擇山地直接受損區、丘陵直接受損區對應的參照區C1,H1及用于檢驗的參照區C2,H2。具體選擇方法:① 在研究區東部的非礦區分別劃定氣候、地理等條件與山地直接受損區、丘陵直接受損區相近的區域;② 在劃定的區域內選擇出與受損區未受損時農業生產及生產方式相近的區域作為潛在參照區;③ 受損區農田97.92%分布在褐土分布區,因此剔除掉潛在參照區中的紅黏土分布區;④ 將DEM及由DEM計算得到的坡度、坡向等在ArcGIS10.5中進行疊加分析[25],在潛在參照區中選取與直接受損區相近的生境單元。⑤ 在選取的生境單元中尋找集中連片的區域作為最終的參照區(圖3)。
煤炭開采會對礦區周邊環境造成負面影響,植被的生長情況最能反映礦區的生態環境,因長河流域的農作物種類基本相同,故以礦區周邊農田NDVI指數的變化來確定間接受損區。農田NDVI不僅受到煤炭開采的影響,也受到土地利用差異、氣候變化、地形條件等其他因素的影響,為了削弱其他因素的影響:① 水域、建設用地、未利用地的NDVI指數與農田NDVI指數差異較大,結合高分辨率的Google Earth衛圖進行識別,剔除農田邊緣的水域用地、建設用地、未利用地的混合像元,利用純凈的農田像元來確定NDVI的空間分布情況。② 研究區面積較小,降水量、氣溫等氣候條件相近,為了更進一步的減小降水量、氣溫等差異,分區測算山地直接受損區、丘陵直接受損區周邊的農田NDVI指數的變化。③ 地學研究中,空間交互受到距離衰減效應的影響,地理實體距離越近,交互強度越大[26],大面積直接受損區域與附近小范圍的直接受損碎斑(碎斑為幾個像元點大小)的空間交互作用顯著,為減小空間交互作用的影響,忽略此類小碎斑對周邊區域的影響。
確定間接受損區的具體步驟如下:① 利用2015年Landsat影像計算礦區農田NDVI指數,并剔除混合像元。② 從山地直接受損區外緣的像元(30 m×30 m)開始向外搜索,每30 m統計一次NDVI值小于參照區C1的農田NDVI均值像元的頻數,若搜索到丘陵,則統計NDVI值小于參照區H1的農田NDVI均值像元的頻數,直到此頻數變化趨于穩定為止。③ 從丘陵直接受損區外緣的像元開始向外搜索,每30 m統計一次NDVI值小于參照區H1的農田NDVI均值像元的頻數,若搜索到山地,則統計NDVI值小于參照區C1的農田NDVI均值像元的頻數,直到此頻數變化趨于穩定為止。④ 將像元頻數趨于穩定的臨界距離定義為煤炭開采對周邊環境的影響范圍并劃入間接受損區。
圖6為長河流域煤炭開采對礦區周邊環境的影響情況,隨著與直接受損區邊界距離的增加,直接受損區周邊像元NDVI值小于參照區NDVI均值的像元頻數先波動下降后趨于穩定。經計算,像元頻數距山地直接受損區990 m、丘陵直接受損區660 m時趨于穩定,且此距離附近的植被NDVI均值分別接近參照區C1和H1的NDVI均值,表明植被隨著與直接受損區距離的增大受煤炭開采的影響逐漸減小,距離受損區990 m以外的山地、660 m以外的丘陵基本不受煤炭開采的影響。因此,本研究劃定煤炭開采對山地直接受損區、丘陵直接受損區周邊環境的影響距離分別為990,660 m,劃定后的間接受損區面積為18.99 km2(圖5)。

圖6 煤炭開采對周邊環境的影響范圍Fig.6 Disturbance distance for the coal mining to the surrounding environment
由于無法實地測量煤炭開采以來各年的植被固碳量,因此采用覆蓋整個研究區的多年遙感數據結合生產力模型來模擬農田植被NPP。由于研究區面積較小,植被生長季的變化很快,需要較高時間、空間分辨率的遙感影像對其進行監測。采用MODIS數據盡管可以模擬逐月的植被NPP,但受空間分辨率的限制,其NPP估算結果的空間表達不夠精細。中高分辨率影像受時間序列和云雨天氣的影響,只能做單時相或單一季節的NPP模擬。目前,遙感衛星獲取數據的時空分辨率還不足以同時達到在空間和時間上精確監測植被動態變化的目的[27]。因此,筆者利用遙感數據時空融合技術,對MODIS數據和Landsat數據進行融合獲取較高時空分辨率的遙感數據來模擬農田植被NPP[28]。
4.1.1遙感數據融合
本研究采用ZHU等提出的增強自適應的遙感圖像時空融合方法(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model,ESTARFM)[29],利用MODIS的多時相信息與Landsat的高空間細節信息,合成高時空分辨率影像。首先,基于tk時刻獲取的MODIS,Landsat影像和tp時刻獲取的MODIS影像,通過計算影像間相似像元的空間分布差異及波譜差異(空間距離權重、光譜距離權重),實現tp時刻的Landsat影像尺度滑動窗口內中心像元的反射率預測[29];然后,根據tm,tn與tp之間的MODIS影像反射率的差異,分別計算tm和tn相對于tp的時間距離權重[29];最后,用tm和tn計算的tp時刻的高分辨率影像和時間距離權重來計算tp時刻最終的融合影像[29],其中tp時刻中心像元反射率可表示為
F(xw/2,yw/2,tp,B)=F(xw/2,yw/2,tk,B)+

(1)
式中,F(xw/2,yw/2,tk,B)為tk時刻的高分辨率像元值;B為影像參與計算的波段;(xw/2,yw/2)為窗口的中心像元;N為相似像元的數目;C(xi,yi,tp,B)和C(xi,yi,tk,B)為tp和tk時刻的低分辨率像元值;Wi和Vi為窗口內各相似像元的權重和轉換系數。算法中權重函數和轉換系數可以表示為

(2)
Vi=(Fin-Fim)/(Cn-Cm)
(3)
式中,Di為由窗口中心的預測像元與窗口內相似臨 近像元的光譜距離和空間距離計算的產品;Cn,Cm,Fin,Fim分別為不同時間(tm和tn)相應的低分辨率和高分辨像元值。
為使結果更準確,使用兩期數據的時間距離權重,計算tp時刻的融合影像:
F(xw/2,yw/2,tp)=TmFm(xw/2,yw/2,tp)+
TnFn(xw/2,yw/2,tp)
(4)
式中,Tm,Tn分別為tm和tn相對于tp的時間距離權重。
4.1.2NDVI計算
選用Landsat影像及融合后的影像構建2000—2015年逐月NDVI值[28],計算公式為
NDVI=(RNI-R)/(RNI+R)
(5)
式中,RNI為近紅外波段;R為紅光波段。
4.1.3NPP估算模型
CASA(Carnegie Ames Stanford Approach)模型是NPP估算模型中的一種光能利用率模型[31],由NDVI、溫度、降水量、太陽輻射、植被類型及土壤類型等共同驅動[32]。該模型中,NPP由植被的光合有效輻射和光能利用率2個因子來表示,其計算公式為
NPP(x,t)=APAR(x,t)ε(x,t)
(6)
式中,APAR(x,t)為像元x在t月吸收的光合有效輻射;ε(x,t)為像元x在t月的光能利用率。
植被的光合有效輻射取決于太陽總輻射和植被對光合有效輻射的吸收比例,公式為
APAR(x,t)=SOL(x,t)FPAR(x,t)×0.5
(7)
式中,常數0.5為植被所能利用的太陽有效輻射占太陽總輻射的比例;SOL(x,t)為t月份像元x處的太陽總輻射量,由Penman-Monteith模型計算[33];FPAR(x,t)為植被層對入射光合有效輻射的吸收比例,可采用:
FPAR(x,t)=FPARNDVI×FPARSR×0.5
(8)
式中,常數0.5為比例系數;FPARNDVI為由NDVI估算的FPAR,NDVI數據由式(5)獲得;FPARSR為由比值植被指數SR估算的FPAR。公式為
(FRARmax-FPARmin)+FPARmin
(9)
(FRARmax-FPARmin)+FPARmin
(10)
SR(x,t)=[1+NDVI(x,t)]/[1-NDVI(x,t)]
(11)
在理想條件下,植被具有最大光能利用率,在現實條件下最大光能利用率主要受氣溫和水分的影響,實際的光能利用率ε(x,t)可估算為
ε(x,t)=εmaxTε1(x,t)Tε2(x,t)Wε(x,t)
(12)
其中,εmax為理想條件下的最大光能利用率;Tε1(x,t)和Tε2(x,t)分別為低溫和高溫對光能利用率的脅迫作用;Wε(x,t)為水分脅迫影響系數。Tε1(x,t),Tε2(x,t)和Wε(x,t)計算公式參照文獻[34]。驅動模型所用到的太陽總輻射、氣溫、降水量等氣象數據采用ANUSPLIN法進行空間插值生成[35]。
土壤有機碳與植被密切相關,植被凋落物經微生物分解,可以轉化為土壤中的有機質。煤炭開采導致植被破壞,勢必引起土壤有機碳含量發生變化。本研究于2015年7月,利用GPS定位技術,本著采樣點集中不分散的原則,按照區域農田土壤1 km×1 km采樣網格布局法確定采樣點,采樣時以網格中心為原點、5 m為半徑的范圍內,結合礦區的實際地形,在研究區內采樣116個(圖2)。土壤采樣時,以蛇形法設置3個點,采用多點混合取樣[36],每個點采集耕層(0~20 cm)、中層(20~40 cm)、底層(40~60 cm)共3個深度的土樣,密封帶回實驗室,采用環刀法測量土壤容重[37],采用重鉻酸鉀容量法測定土壤有機碳[38]。
土壤有機碳密度被用于衡量單位面積一定深度土體中有機碳儲量[39],計算公式為
SOCj=tjbjOj(1-Vj)
(13)
其中,SOCj為第j層土壤碳密度,j為土壤層;tj為土層厚度;bj為土壤容重;Oj為土壤有機碳含量;Vj指第j層土壤中直徑>2 mm的礫石所占的體積比,由于研究區農田土壤含礫石量基本為0,可忽略不計。
本研究于2013年在研究區農田選擇了2塊等面積的樣地,分別位于直接受損區和間接受損區,樣地均為小麥-玉米輪作的旱地,土地理化性質和農田管理措施一致,植被各部分比重采用重鉻酸鉀容量法—外加熱法測定,農田植被NPP采用“收割法”測定[40]。為了驗證植被NPP模擬值的準確性,將2013年模型模擬結果與2013年實測值進行對比。由表1可知,模擬的結果與實測值差異較小,且植被NPP值均為直接受損區小于間接受損區,表明本研究估算的NPP基本可靠。
表1 NPP模擬值與實測值
Table 1 Simulated and measured NPP values gC/m2

樣地丘陵直接受損區丘陵間接受損區本研究258.57268.42實測值260.21269.08
圖7(a)為模型模擬得出的農田植被年均NPP分布圖,長河流域西部礦區的農田植被NPP小于東部非礦區。假設礦區農田土地利用類型未發生變化,則2000—2015年因煤礦開采導致的受損區植被固碳量為3 936.16 t/km2,直接受損區植被固碳量為3 867.99 t/km2,間接受損區植被固碳量為3 954.33 t/km2(表2)。圖7(b)~(d)為通過克里金插值獲得研究區農田土壤有機碳的空間分布圖,長河流域西北部土壤有機碳含量較低,中部、東部土壤有機碳含量較高,受秸稈還田等農田管理措施的影響,農田土壤有機碳含量為耕層>中層>底層。將土壤有機碳單位與植被固碳量單位統一,假設礦區農田土地利用類型未發生變化,則2000—2015年因煤碳開采導致的受損區土壤碳儲存量為6 737.05 t/km2,直接受損區土壤碳儲存量為6 438.93 t/km2,間接受損區土壤碳儲存量為6 816.49 t/km2(表2)。礦區農田生態系統植被固碳量小于土壤碳儲存量,主要是由于農作物在生長期中固定了大量的碳,其中一部分以農產品的形式輸出,剩余的基本進行秸稈還田處理進入大氣或者土壤,所以農田土壤中儲存了更多的有機碳。

圖7 農田植被年均NPP與土壤有機碳空間分布Fig.7 Spatial distribution map of mean annual NPP and soil organic carbon in farmland
受損區內農田植被受煤炭開采、氣溫、降水量和地表人類活動的共同影響;參照區受氣溫、降水量和地表人類活動的影響。朱阿興等提出兩個區域的地理環境越相似,這兩個區域的目標變量值越相近[41],由于整個研究區面積較小,經計算受損區、參照區的平均氣溫差值為0.2 ℃,平均降水量差值為3.1 mm,氣候條件相近;地表人類活動主要是種植模式、農業科技水平、施肥方式和施肥量,多年來研究區種植模式為小麥-玉米輪作,農業科技水平相同,施肥方式和施肥量相近。因此,參照區與受損區未受損時的農田生態系統固碳能力是相近的,礦區農田固碳損失量主要受煤炭開采活動的影響。這樣,本研究確定的參照區可很好的模擬受損區未受損時的情景。同時,煤炭開采導致研究區土地利用類型發生變化,而土地利用變化(農田轉變為非農田)是對區域農田生態系統固碳影響最顯著的因素,評估固碳損失量必須考慮到土地利用變化的影響。煤炭開采后部分農田轉變為采礦用地和沉陷積水區,農田面積減少。因此,煤炭開采導致的農田固碳損失量應包括農田轉化的采礦用地、沉陷積水區的固碳量。煤炭開采導致的礦區農田生態系統固碳損失量計算流程:
表2 農田生態系統固碳量
Table 2 Carbon sequestration in farmland ecosystems

t/km2
假設礦區農田土地利用類型未發生變化,則土地利用現狀為農田的植被固碳損失量計算公式為
Lv=(NPPr-NPPd)Sd+(NPPr-NPPi)Si
(14)
式中,NPPr,NPPd和NPPi分別為參照區、直接受損區和間接受損區的NPP均值;Sd,Si分別為直接受損區和間接受損區面積。
由于開采前土壤有機碳含量已無法獲取,因此以受損區與參照區內土壤有機碳的差值來計算土壤有機碳的損失量。假設礦區農田土地利用類型未發生變化,則土地利用現狀為農田的土壤碳儲存損失量計算為
Lp=(SOCr-SOCd)Sd+(SOCr-SOCi)Si
(15)
式中,SOCr,SOCd和SOCi分別為參照區、直接受損區和間接受損區的農田土壤碳密度(式(13))。
考慮到農田土地利用變化對農田生態系統固碳損失量的影響,礦區農田生態系統固碳損失量:
L=Lv+Lp+UA
(16)
式中,U為農田土地利用類型未發生變化時受損區農田生態系統單位面積固碳量(表2);A為2000—2015年農田面積減少量,即農田轉化為采礦用地、沉陷積水區的面積(由2000,2015年土地利用現狀疊加估算,2000—2015年直接受損區內農田面積減少約0.62 km2,間接受損區農田面積減少約0.78 km2)。
在相近的氣候、地貌、土壤、耕作方式、農業科技水平和施肥制度下,煤炭開采導致直接受損區地表沉陷,農田植被的生長環境發生了嚴重破壞,甚至在雨季出現長期的積水,植被生命活動受到抑制,植被生產力下降甚至徹底消失,進而導致土壤有機碳含量下降,單位面積的農田生態系統固碳量降低,由式(13)~(16)得,直接受損區農田生態系統固碳損失量為12 947.55 t,其中因農田面積減少造成的固碳損失為6 390.29 t,土地利用現狀為農田的固碳損失量為6 557.26 t;間接受損區受煤炭開采導致的煤煙、粉塵以及有害氣體的影響,植被生長、土壤質量受到的負面影響相較直接受損區要小,故直接受損區固碳降低幅度大于間接受損區,由式(13)~(16)得,間接受損區農田生態系統固碳損失量為17 790.13 t,其中因農田面積減少造成的固碳損失為8 401.24 t,現有農田固碳損失量為9 388.89 t。綜上所述,煤炭開采造成的生態環境破壞及農田面積減少是礦區農田生態系統固碳損失的重要原因。
總之,與未受煤炭開采得參照區相比,在整個煤炭開采的受損區范圍內,農田生態系統固碳損失總量為30 737.68 t,其中直接受損區農田生態系統固碳損失量為12 947.55 t,占總損失的42.12%,間接受損區農田生態系統固碳損失量為17 790.13 t,占總損失的57.88%。本研究估算的煤炭開采導致的礦區農田生態系統固碳損失精度相對較高,這主要是因為:以往研究忽略了間接受損區,本研究在減弱了其他因素帶來影響的同時,量化了間接受損區,并評估了直接受損區和間接受損區的固碳損失量,使礦區生態系統固碳損失評估更為精確;煤炭開采會導致原始農田發生土地利用類型上的改變,本研究考慮了土地利用類型變化對農田植被碳庫的影響,保證了煤炭開采對農田生態系統固碳損失量評估的精度。煤炭開采對區域農田生態系統固碳量的影響屬于失碳效應。
為了驗證本研究方法的適用性:① 采用參照區M2,H2來估算間接受損范圍及農田生態系統固碳損失量,經計算,煤炭開采對直接受損區周邊環境的影響范圍為山地990 m、丘陵660 m,影響范圍與上文結果一致;煤炭開采導致的農田生態系統固碳損失量為29 274.57 t,與上文結果僅相差4.76%,即使同一區域的植物自身生產力也會存在一定的差異,故不同參照區計算農田生態系統固碳損失量也會帶來小范圍的誤差。② 直接受損區位于長河流域的8個煤礦區中,將8個煤礦的直接受損區、間接受損區農田生態系統單位面積固碳損失量與研究區進行對比,相對偏差分布如圖8所示。根據礦區資料顯示天安岳圣山、天泰坤達煤礦進行了復墾,單位面積固碳損失量較小;天安盈盛開采規模擴大,生態破壞加劇,單位面積固碳損失較大,其余煤礦區與研究區整體單位面積固碳損失量的相對偏差較小,表明該方法評估煤礦開采導致的礦區農田生態系統固碳損失量是適用的。

圖8 煤礦區農田固碳損失量Fig.8 Carbon sequestration loss of farmland in coal mine area
煤炭開采活動損壞了原生地貌和地表植被,造成土地利用類型變化,引起植被-土壤碳匯量減少,導致了農田生態系碳儲量的大量流失[12,17,42]。當前,有學者采用多年礦區植被NPP的差值來估算煤炭開采導致的生態損失[13],但由于各年氣溫、降水等環境因素的變化較大,不能較好的說明礦區生態損失由煤炭開采所引起;也有學者采用沉陷區周邊區域作為參照區來評估沉陷區植被固碳或土壤碳儲存的損失[14,43-44],雖然相同年份內沉陷區與參照區氣溫、降水等環境因素基本一致,然而參照區本身已受到煤炭開采的間接影響,使得評估的結果偏低。本研究在基于傳統的開采沉陷、采空及損毀土地確定礦區受損范圍的基礎上,提出了測算周邊NDVI指數變化趨勢以量化間接受損區的方法。該方法可適用于大范圍、長時間序列的地表植被覆蓋區,可更準確的確定煤炭開采對地表生態的影響邊界。同時,本研究考慮了土地利用類型變化對農田生態系統碳庫的影響,提高了礦區農田生態系統固碳損失的評估精度,也為研究礦區時空生態系統固碳累積變化及后續生態補償機制等問題提供借鑒。
已有研究表明,合理的土地復墾及差別整治有利于礦區生態系統碳匯量增加,一定程度上可抑制煤炭開采造成的生態環境破壞[10,45]。利用遙感影像來量化礦區生態受損范圍,可精確估算礦區固碳損失量,動態監測礦區生態環境及未來土地復墾的生態效益與碳匯量變化。因此,結合礦區農田生態系統固碳損失程度,科學合理的進行地貌重塑、土壤重構、植被重建及生物多樣性恢復等,對于提升農田生產能力、重建農田生態系統、緩解氣候變暖具有重要意義。本研究對長河流域2000—2015年的農田生態系統直接受損區和間接受損區的固碳損失量進行了估算,奠定了礦區碳增匯經濟潛力評價的基礎,也為礦區土地復墾、生態修復與重建提供更加準確和實際的指導。然而,本研究提出的評估方法僅適用于中國華北丘陵山區,針對其他區域的實施需依據實際情況進行調整。同時,在今后的研究中應建立長期的土壤養分監測站點,以確定煤炭開采情況下農田土壤碳儲存的動態變化及礦區農田生態系統固碳損失的機理機制。
(1)研究區內由于煤炭開采導致的農田直接受損面積為5.06 km2,間接受損面積為18.99 km2。在間接受損區內,生態環境的損毀程度隨著與直接受損區距離的增大而逐漸減弱。
(2)礦區農田生態系統植被固碳量為3 936.16 t/km2,小于土壤碳儲存量的6 737.05 t/km2。
(3)研究區因煤炭開采導致的農田生態系統固碳損失量為30 737.68 t。其中,直接受損區固碳損失量為12 947.55 t,間接受損區固碳損失量為17 790.13 t。煤炭開采對礦區農田生態系統固碳量的影響屬于失碳效應。
(4)煤炭開采導致的生態系統固碳損失評估及生態環境損毀的研究中,間接受損區也應該納入評估范圍。