田惠文，張欣欣，畢如田，朱洪芬，席 茜

(1.山西農(nóng)業(yè)大學(xué) 資源環(huán)境學(xué)院,山西 晉中 030801; 2.中國地質(zhì)大學(xué) 地球物理與空間信息學(xué)院,湖北 武漢 430074)

陸地生態(tài)系統(tǒng)碳循環(huán)在全球碳收支中起著主導(dǎo)作用[1]，而農(nóng)田生態(tài)系統(tǒng)是構(gòu)成陸地生態(tài)系統(tǒng)的三大系統(tǒng)之一，且農(nóng)田生態(tài)系統(tǒng)碳庫是全球碳庫和陸地生態(tài)系統(tǒng)碳庫的重要組成部分，在全球碳平衡中扮演著重要角色[2]。農(nóng)田生態(tài)系統(tǒng)受人類活動(dòng)影響最為顯著，其中不同的管理措施、農(nóng)田的棄耕是學(xué)者關(guān)注的重點(diǎn)，如石小霞等對華北平原不同農(nóng)田管理措施與土壤碳庫的變化關(guān)系進(jìn)行了研究[3];王海候等研究了不同秸稈還田模式與土壤碳庫特性之間的關(guān)系[4];趙威等對豫西丘陵坡地棄耕農(nóng)田植被演替對土壤碳庫的影響開展了研究[5]。此外，水土流失、土地鹽堿化及土地利用變化等也會(huì)對農(nóng)田碳庫產(chǎn)生擾動(dòng)，且這些因素都不可避免地耦合著人類活動(dòng)的影響。

工業(yè)革命以來，人類生產(chǎn)方式發(fā)生了翻天覆地的變化，空氣中CO2等氣體含量不斷上升，導(dǎo)致了溫室效應(yīng)和一系列生態(tài)環(huán)境問題[6]，其中煤炭開采是重要的原因之一。據(jù)估算，在1989—2010年由煤礦開采造成的土地?fù)p毀面積達(dá)100.01×104hm2[7]，礦區(qū)生態(tài)環(huán)境面臨著嚴(yán)峻考驗(yàn)。農(nóng)田生態(tài)系統(tǒng)碳庫主要包括植被碳庫和土壤碳庫[8-9]，近年來，有學(xué)者對煤礦區(qū)土壤碳匯進(jìn)行了研究，如余建等研究了徐州市柳新鎮(zhèn)高潛水位礦區(qū)不同復(fù)墾及利用方式的土壤有機(jī)碳、無機(jī)碳的構(gòu)成與變化特征[10];程靜霞等以焦作礦區(qū)為例，研究了煤炭開采沉陷區(qū)土壤有機(jī)碳的空間變化[11]。也有學(xué)者對礦區(qū)植被碳匯展開研究，如侯湖平等以徐州九里礦區(qū)為例，研究發(fā)現(xiàn)煤炭開采對區(qū)域農(nóng)田植被碳儲量的影響屬于失碳效應(yīng)[12];侯湖平等研究了煤礦區(qū)土地利用變化對生態(tài)系統(tǒng)植被碳儲量的影響，發(fā)現(xiàn)礦區(qū)碳總量降低[13];黃翌等以山西大同礦區(qū)為例，研究發(fā)現(xiàn)多年來礦區(qū)植被-土壤系統(tǒng)碳匯量減少[14];徐占軍等采用DNDC模型研究了煤田開采沉陷對農(nóng)田土壤和植被碳庫的擾動(dòng)[15]。煤礦區(qū)作為一個(gè)特殊的復(fù)雜地理區(qū)域，長時(shí)間、大規(guī)模和高強(qiáng)度的煤炭開采不可避免地破壞原生的礦床地質(zhì)條件，引發(fā)植被退化、水土流失及物種減少等[16-17]，尤其是地下開采造成了大面積塌陷地、土地沉陷和耕地破壞等生態(tài)環(huán)境問題，嚴(yán)重制約了農(nóng)田生態(tài)系統(tǒng)碳庫的積累[18-20]。煤礦開采、運(yùn)輸過程中也會(huì)產(chǎn)生大量的煤煙、粉塵以及有害氣體，對周邊的植被生長、土壤質(zhì)量及居民生活造成負(fù)面效應(yīng)[21]。目前，對于礦區(qū)農(nóng)田生態(tài)系統(tǒng)的碳儲存研究主要集中在礦區(qū)沉陷、采空區(qū)、損毀土地等直接受損區(qū)域上的植被固碳或土壤碳儲存，同時(shí)考慮植被-土壤碳匯變化的研究較少且忽略了煤炭開采對直接受損區(qū)以外的周邊地區(qū)。

因此，完善煤炭開采導(dǎo)致的農(nóng)田生態(tài)系統(tǒng)固碳損失的評估方法不僅對礦區(qū)土地復(fù)墾、生態(tài)修復(fù)有著極其重要的意義，還可為煤礦區(qū)生態(tài)補(bǔ)償測算提供決策與意見。筆者以山西省晉城市長河流域?yàn)槔?，對煤炭開采導(dǎo)致的農(nóng)田生態(tài)系統(tǒng)固碳損失(植被固碳和土壤碳儲存)進(jìn)行評估。本研究的主要目的:建立一種適用于煤炭開采導(dǎo)致農(nóng)田生態(tài)系統(tǒng)固碳損失的方法;量化煤炭開采對周邊農(nóng)田的間接受損范圍;估算礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失，為礦區(qū)土地復(fù)墾與生態(tài)修復(fù)提供依據(jù)。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

長河流域位于山西省晉城市澤州縣西北部，在112.63°～112.77°N，35.50°～35.63°E之間，下轄3個(gè)鄉(xiāng)鎮(zhèn)，土地利用類型主要有林地、草地、農(nóng)田、采礦用地、建設(shè)用地、沉陷積水區(qū)等，總面積約113.16 km2，其中農(nóng)田面積為64.55 km2(圖1)。研究區(qū)屬溫帶大陸性季風(fēng)氣候，年平均氣溫10.6 ℃，年降雨量550～600 mm，地上農(nóng)作物類型主要是玉米、冬小麥等，農(nóng)田耕作方式為小麥-玉米輪作，施肥方式為傳統(tǒng)的土壤施肥，施肥量基本相同，約為N肥1.5 t/km2，K肥0.9 t/km2，P肥0.45 t km2。研究區(qū)海拔約857 m，地形以山地、丘陵為主，境內(nèi)土壤類型為褐土和紅黏土，其中褐土面積占91.14%。長河流域位于山西沁水煤田界內(nèi)，流域內(nèi)共有天泰坤達(dá)等11個(gè)煤礦，煤炭資源豐富。2000—2015年因煤炭開采造成的采空區(qū)、沉陷區(qū)、損毀土地等面積約10.29 km2，造成了植被覆蓋度降低、水土流失和環(huán)境污染等一系列嚴(yán)重的生態(tài)問題。

圖1 研究區(qū)Fig.1 Study area

1.2 數(shù)據(jù)及其預(yù)處理

(1)礦區(qū)數(shù)據(jù)。長河流域各礦業(yè)提供的采空區(qū)、已損毀和已塌陷等土地現(xiàn)狀分布數(shù)據(jù)。

(2)土地利用數(shù)據(jù)。來源于2000年和2015年澤州縣土地利用規(guī)劃空間數(shù)據(jù)庫的1∶1萬土地覆蓋數(shù)據(jù)。

(3)農(nóng)田數(shù)據(jù)。提取自澤州縣土地利用規(guī)劃空間數(shù)據(jù)庫的1∶1萬土地覆蓋數(shù)據(jù)。

(4)施肥制度。來源于澤州縣耕地質(zhì)量監(jiān)測評價(jià)數(shù)據(jù)，共16個(gè)樣點(diǎn)。

(5)耕作方式。提取自山西省測土配方施肥數(shù)據(jù)，共528個(gè)樣點(diǎn)。

(6)氣象數(shù)據(jù)。來源于中國氣象數(shù)據(jù)網(wǎng)，包括氣溫、降水量等，時(shí)間序列為2000—2015年。

(7)遙感數(shù)據(jù)。來源于美國地質(zhì)勘探局(United States Geological Survey,USGS)網(wǎng)站的MODIS和Landsat影像，時(shí)間序列為2000—2015年。數(shù)字高程模型(DEM)由各礦區(qū)提供，精度為1 m，并采用DEM數(shù)據(jù)計(jì)算坡度、坡向。

(8)土壤類型數(shù)據(jù)。來源于山西省第2次土壤普查數(shù)據(jù)庫。

(9)地貌類型數(shù)據(jù)。來源于澤州縣國土資源局。

(10)植被數(shù)據(jù)。2013年實(shí)測的植被凈初級生產(chǎn)力(Net Primary Production,NPP)數(shù)據(jù)，實(shí)測方法見4.3。

(11)土壤數(shù)據(jù)。2015年7月在長河流域采樣，并在實(shí)驗(yàn)室測試了土壤有機(jī)碳、土壤容重等數(shù)據(jù)。

本文研究對象為農(nóng)田，所有數(shù)據(jù)均采用農(nóng)田矢量數(shù)據(jù)進(jìn)行了裁剪。以上數(shù)據(jù)均經(jīng)空間數(shù)據(jù)格式轉(zhuǎn)換、校正并統(tǒng)一投影轉(zhuǎn)化為CGCS2000坐標(biāo)系(圖1，2)。

1.3 研究方法

首先利用礦區(qū)數(shù)據(jù)，提取2000年以來因煤炭開采造成的采空、損毀及塌陷范圍，作為礦區(qū)直接受損區(qū);然后綜合考慮海拔、坡度、土壤類型、到礦區(qū)距離等因素，利用空間疊加分析，選擇出與礦區(qū)生境相似且未受煤炭開采干擾的區(qū)域作為參照區(qū);基于Landsat影像計(jì)算歸一化差值植被指數(shù)(Normalized Difference Vegetation Index,NDVI)，分析直接受損區(qū)邊緣向外的NDVI指數(shù)變化特征，確定間接受損區(qū);模擬長河流域農(nóng)田植被NPP，測定農(nóng)田土壤有機(jī)碳，并與受損區(qū)、參照區(qū)范圍疊加，計(jì)算得受損區(qū)農(nóng)田植被固碳損失量和土壤碳儲存損失量;結(jié)合農(nóng)田的土地利用變化情況，估算礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量，技術(shù)路線如圖3所示。

2 固碳機(jī)制

固碳主要是指生態(tài)系統(tǒng)的碳凈吸收量，其主要是生態(tài)系統(tǒng)植被固碳和土壤碳儲存，即來自于植物光合作用固定的碳及枯枝落葉等進(jìn)入土壤形成的有機(jī)碳[22]。從生態(tài)系統(tǒng)碳固定的角度來看，植被—土壤可以看作一個(gè)統(tǒng)一的系統(tǒng)，植物的光合作用是該系統(tǒng)的碳來源，植被與土壤中的呼吸作用為碳支出，剩余的轉(zhuǎn)化為土壤有機(jī)質(zhì)(圖4)。在長河流域農(nóng)田生態(tài)系統(tǒng)中，植被固碳量近似地以植被NPP表示(植物光合作用扣除自氧呼吸后的差值)，土壤有機(jī)碳通過樣品采集，實(shí)驗(yàn)室測試獲得。因此，只需要獲取煤炭開采以來農(nóng)田受損區(qū)及參照區(qū)的植被NPP及土壤有機(jī)碳，就可以估算礦區(qū)農(nóng)田生態(tài)系統(tǒng)的碳損失量。

圖2 數(shù)據(jù)源Fig.2 Data source

圖3 技術(shù)路線Fig.3 Technology flow chart

圖4 固碳機(jī)制Fig.4 Mechanism of carbon sequestration

3 受損區(qū)及參照區(qū)的確定

土地利用變化(農(nóng)田轉(zhuǎn)變?yōu)榉寝r(nóng)田)是對區(qū)域農(nóng)田生態(tài)系統(tǒng)固碳影響最顯著的因素[23-24]。由于本研究劃定的受損區(qū)為農(nóng)田，故農(nóng)田轉(zhuǎn)化為其他土地利用類型的區(qū)域不歸入受損區(qū)，但估算礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失時(shí)，仍考慮因煤炭開采導(dǎo)致的農(nóng)田面積減少造成的固碳損失。在土地利用類型未發(fā)生改變的情況下，影響農(nóng)田生態(tài)系統(tǒng)固碳量的常規(guī)影響因素包括:氣候、地貌類型、土壤類型、耕作方式、施肥制度以及農(nóng)業(yè)科技水平等[13]。本研究假定土地利用類型未發(fā)生改變，充分考慮常規(guī)因素的影響來確定參照區(qū)和受損區(qū)，最后結(jié)合農(nóng)田土地利用類型的轉(zhuǎn)化情況估算煤炭開采導(dǎo)致的農(nóng)田生態(tài)系統(tǒng)固碳損失量。

3.1 直接受損區(qū)

從礦區(qū)各個(gè)礦業(yè)公司提供的礦區(qū)數(shù)據(jù)中，提取2000年以后由于煤炭開采造成的采空、塌陷、損毀等土地?cái)?shù)據(jù)作為直接受損區(qū)，經(jīng)計(jì)算，山地直接受損區(qū)面積為3.52 km2，丘陵直接受損區(qū)面積為1.54 km2(圖5)。

圖5 受損區(qū)和參照區(qū)Fig.5 Damaged area and reference area

3.2 參照區(qū)

固碳損失量估算的關(guān)鍵在于參照區(qū)的選擇。受損區(qū)位于山地和丘陵，由圖2可知，長河流域的耕作方式、作物種類、施肥制度及農(nóng)業(yè)科技水平基本相同，但山地和丘陵農(nóng)田的灌溉條件、管理措施不同，為了準(zhǔn)確的評估與檢驗(yàn)，本研究分別選擇山地直接受損區(qū)、丘陵直接受損區(qū)對應(yīng)的參照區(qū)C1,H1及用于檢驗(yàn)的參照區(qū)C2,H2。具體選擇方法:① 在研究區(qū)東部的非礦區(qū)分別劃定氣候、地理等條件與山地直接受損區(qū)、丘陵直接受損區(qū)相近的區(qū)域;② 在劃定的區(qū)域內(nèi)選擇出與受損區(qū)未受損時(shí)農(nóng)業(yè)生產(chǎn)及生產(chǎn)方式相近的區(qū)域作為潛在參照區(qū);③ 受損區(qū)農(nóng)田97.92%分布在褐土分布區(qū)，因此剔除掉潛在參照區(qū)中的紅黏土分布區(qū);④ 將DEM及由DEM計(jì)算得到的坡度、坡向等在ArcGIS10.5中進(jìn)行疊加分析[25]，在潛在參照區(qū)中選取與直接受損區(qū)相近的生境單元。⑤ 在選取的生境單元中尋找集中連片的區(qū)域作為最終的參照區(qū)(圖3)。

3.3 間接受損區(qū)

煤炭開采會(huì)對礦區(qū)周邊環(huán)境造成負(fù)面影響，植被的生長情況最能反映礦區(qū)的生態(tài)環(huán)境，因長河流域的農(nóng)作物種類基本相同，故以礦區(qū)周邊農(nóng)田NDVI指數(shù)的變化來確定間接受損區(qū)。農(nóng)田NDVI不僅受到煤炭開采的影響，也受到土地利用差異、氣候變化、地形條件等其他因素的影響，為了削弱其他因素的影響:① 水域、建設(shè)用地、未利用地的NDVI指數(shù)與農(nóng)田NDVI指數(shù)差異較大，結(jié)合高分辨率的Google Earth衛(wèi)圖進(jìn)行識別，剔除農(nóng)田邊緣的水域用地、建設(shè)用地、未利用地的混合像元，利用純凈的農(nóng)田像元來確定NDVI的空間分布情況。② 研究區(qū)面積較小，降水量、氣溫等氣候條件相近，為了更進(jìn)一步的減小降水量、氣溫等差異，分區(qū)測算山地直接受損區(qū)、丘陵直接受損區(qū)周邊的農(nóng)田NDVI指數(shù)的變化。③ 地學(xué)研究中，空間交互受到距離衰減效應(yīng)的影響，地理實(shí)體距離越近，交互強(qiáng)度越大[26]，大面積直接受損區(qū)域與附近小范圍的直接受損碎斑(碎斑為幾個(gè)像元點(diǎn)大小)的空間交互作用顯著，為減小空間交互作用的影響，忽略此類小碎斑對周邊區(qū)域的影響。

確定間接受損區(qū)的具體步驟如下:① 利用2015年Landsat影像計(jì)算礦區(qū)農(nóng)田NDVI指數(shù)，并剔除混合像元。② 從山地直接受損區(qū)外緣的像元(30 m×30 m)開始向外搜索，每30 m統(tǒng)計(jì)一次NDVI值小于參照區(qū)C1的農(nóng)田NDVI均值像元的頻數(shù)，若搜索到丘陵，則統(tǒng)計(jì)NDVI值小于參照區(qū)H1的農(nóng)田NDVI均值像元的頻數(shù)，直到此頻數(shù)變化趨于穩(wěn)定為止。③ 從丘陵直接受損區(qū)外緣的像元開始向外搜索，每30 m統(tǒng)計(jì)一次NDVI值小于參照區(qū)H1的農(nóng)田NDVI均值像元的頻數(shù)，若搜索到山地，則統(tǒng)計(jì)NDVI值小于參照區(qū)C1的農(nóng)田NDVI均值像元的頻數(shù)，直到此頻數(shù)變化趨于穩(wěn)定為止。④ 將像元頻數(shù)趨于穩(wěn)定的臨界距離定義為煤炭開采對周邊環(huán)境的影響范圍并劃入間接受損區(qū)。

圖6為長河流域煤炭開采對礦區(qū)周邊環(huán)境的影響情況，隨著與直接受損區(qū)邊界距離的增加，直接受損區(qū)周邊像元NDVI值小于參照區(qū)NDVI均值的像元頻數(shù)先波動(dòng)下降后趨于穩(wěn)定。經(jīng)計(jì)算，像元頻數(shù)距山地直接受損區(qū)990 m、丘陵直接受損區(qū)660 m時(shí)趨于穩(wěn)定，且此距離附近的植被NDVI均值分別接近參照區(qū)C1和H1的NDVI均值，表明植被隨著與直接受損區(qū)距離的增大受煤炭開采的影響逐漸減小，距離受損區(qū)990 m以外的山地、660 m以外的丘陵基本不受煤炭開采的影響。因此，本研究劃定煤炭開采對山地直接受損區(qū)、丘陵直接受損區(qū)周邊環(huán)境的影響距離分別為990,660 m，劃定后的間接受損區(qū)面積為18.99 km2(圖5)。

圖6 煤炭開采對周邊環(huán)境的影響范圍Fig.6 Disturbance distance for the coal mining to the surrounding environment

4 農(nóng)田生態(tài)系統(tǒng)固碳損失評估

4.1 植被固碳

由于無法實(shí)地測量煤炭開采以來各年的植被固碳量，因此采用覆蓋整個(gè)研究區(qū)的多年遙感數(shù)據(jù)結(jié)合生產(chǎn)力模型來模擬農(nóng)田植被NPP。由于研究區(qū)面積較小，植被生長季的變化很快，需要較高時(shí)間、空間分辨率的遙感影像對其進(jìn)行監(jiān)測。采用MODIS數(shù)據(jù)盡管可以模擬逐月的植被NPP，但受空間分辨率的限制，其NPP估算結(jié)果的空間表達(dá)不夠精細(xì)。中高分辨率影像受時(shí)間序列和云雨天氣的影響，只能做單時(shí)相或單一季節(jié)的NPP模擬。目前，遙感衛(wèi)星獲取數(shù)據(jù)的時(shí)空分辨率還不足以同時(shí)達(dá)到在空間和時(shí)間上精確監(jiān)測植被動(dòng)態(tài)變化的目的[27]。因此，筆者利用遙感數(shù)據(jù)時(shí)空融合技術(shù)，對MODIS數(shù)據(jù)和Landsat數(shù)據(jù)進(jìn)行融合獲取較高時(shí)空分辨率的遙感數(shù)據(jù)來模擬農(nóng)田植被NPP[28]。

4.1.1遙感數(shù)據(jù)融合

本研究采用ZHU等提出的增強(qiáng)自適應(yīng)的遙感圖像時(shí)空融合方法(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model,ESTARFM)[29]，利用MODIS的多時(shí)相信息與Landsat的高空間細(xì)節(jié)信息，合成高時(shí)空分辨率影像。首先，基于tk時(shí)刻獲取的MODIS,Landsat影像和tp時(shí)刻獲取的MODIS影像，通過計(jì)算影像間相似像元的空間分布差異及波譜差異(空間距離權(quán)重、光譜距離權(quán)重)，實(shí)現(xiàn)tp時(shí)刻的Landsat影像尺度滑動(dòng)窗口內(nèi)中心像元的反射率預(yù)測[29];然后，根據(jù)tm,tn與tp之間的MODIS影像反射率的差異，分別計(jì)算tm和tn相對于tp的時(shí)間距離權(quán)重[29];最后，用tm和tn計(jì)算的tp時(shí)刻的高分辨率影像和時(shí)間距離權(quán)重來計(jì)算tp時(shí)刻最終的融合影像[29]，其中tp時(shí)刻中心像元反射率可表示為

F(xw/2,yw/2,tp,B)=F(xw/2,yw/2,tk,B)+

(1)

式中，F(xiàn)(xw/2,yw/2,tk,B)為tk時(shí)刻的高分辨率像元值;B為影像參與計(jì)算的波段；(xw/2,yw/2)為窗口的中心像元;N為相似像元的數(shù)目；C(xi,yi,tp,B)和C(xi,yi,tk,B)為tp和tk時(shí)刻的低分辨率像元值；Wi和Vi為窗口內(nèi)各相似像元的權(quán)重和轉(zhuǎn)換系數(shù)。算法中權(quán)重函數(shù)和轉(zhuǎn)換系數(shù)可以表示為

(2)

Vi=(Fin-Fim)/(Cn-Cm)

(3)

式中，Di為由窗口中心的預(yù)測像元與窗口內(nèi)相似臨 近像元的光譜距離和空間距離計(jì)算的產(chǎn)品;Cn,Cm,Fin,Fim分別為不同時(shí)間(tm和tn)相應(yīng)的低分辨率和高分辨像元值。

為使結(jié)果更準(zhǔn)確，使用兩期數(shù)據(jù)的時(shí)間距離權(quán)重，計(jì)算tp時(shí)刻的融合影像:

F(xw/2,yw/2,tp)=TmFm(xw/2,yw/2,tp)+

TnFn(xw/2,yw/2,tp)

(4)

式中，Tm,Tn分別為tm和tn相對于tp的時(shí)間距離權(quán)重。

4.1.2NDVI計(jì)算

選用Landsat影像及融合后的影像構(gòu)建2000—2015年逐月NDVI值[28]，計(jì)算公式為

NDVI=(RNI-R)/(RNI+R)

(5)

式中，RNI為近紅外波段;R為紅光波段。

4.1.3NPP估算模型

CASA(Carnegie Ames Stanford Approach)模型是NPP估算模型中的一種光能利用率模型[31]，由NDVI、溫度、降水量、太陽輻射、植被類型及土壤類型等共同驅(qū)動(dòng)[32]。該模型中，NPP由植被的光合有效輻射和光能利用率2個(gè)因子來表示，其計(jì)算公式為

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)

式中，常數(shù)0.5為植被所能利用的太陽有效輻射占太陽總輻射的比例;SOL(x,t)為t月份像元x處的太陽總輻射量，由Penman-Monteith模型計(jì)算[33];FPAR(x,t)為植被層對入射光合有效輻射的吸收比例，可采用:

FPAR(x,t)=FPARNDVI×FPARSR×0.5

(8)

式中，常數(shù)0.5為比例系數(shù);FPARNDVI為由NDVI估算的FPAR，NDVI數(shù)據(jù)由式(5)獲得;FPARSR為由比值植被指數(shù)SR估算的FPAR。公式為

(FRARmax-FPARmin)+FPARmin

(9)

(FRARmax-FPARmin)+FPARmin

(10)

SR(x,t)=[1+NDVI(x,t)]/[1-NDVI(x,t)]

(11)

在理想條件下，植被具有最大光能利用率，在現(xiàn)實(shí)條件下最大光能利用率主要受氣溫和水分的影響，實(shí)際的光能利用率ε(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)為水分脅迫影響系數(shù)。Tε1(x,t),Tε2(x,t)和Wε(x,t)計(jì)算公式參照文獻(xiàn)[34]。驅(qū)動(dòng)模型所用到的太陽總輻射、氣溫、降水量等氣象數(shù)據(jù)采用ANUSPLIN法進(jìn)行空間插值生成[35]。

4.2 土壤碳儲存

土壤有機(jī)碳與植被密切相關(guān)，植被凋落物經(jīng)微生物分解，可以轉(zhuǎn)化為土壤中的有機(jī)質(zhì)。煤炭開采導(dǎo)致植被破壞，勢必引起土壤有機(jī)碳含量發(fā)生變化。本研究于2015年7月，利用GPS定位技術(shù)，本著采樣點(diǎn)集中不分散的原則，按照區(qū)域農(nóng)田土壤1 km×1 km采樣網(wǎng)格布局法確定采樣點(diǎn)，采樣時(shí)以網(wǎng)格中心為原點(diǎn)、5 m為半徑的范圍內(nèi)，結(jié)合礦區(qū)的實(shí)際地形，在研究區(qū)內(nèi)采樣116個(gè)(圖2)。土壤采樣時(shí)，以蛇形法設(shè)置3個(gè)點(diǎn)，采用多點(diǎn)混合取樣[36]，每個(gè)點(diǎn)采集耕層(0～20 cm)、中層(20～40 cm)、底層(40～60 cm)共3個(gè)深度的土樣，密封帶回實(shí)驗(yàn)室，采用環(huán)刀法測量土壤容重[37]，采用重鉻酸鉀容量法測定土壤有機(jī)碳[38]。

土壤有機(jī)碳密度被用于衡量單位面積一定深度土體中有機(jī)碳儲量[39]，計(jì)算公式為

SOCj=tjbjOj(1-Vj)

(13)

其中，SOCj為第j層土壤碳密度,j為土壤層;tj為土層厚度;bj為土壤容重;Oj為土壤有機(jī)碳含量;Vj指第j層土壤中直徑>2 mm的礫石所占的體積比，由于研究區(qū)農(nóng)田土壤含礫石量基本為0，可忽略不計(jì)。

4.3 固碳損失分析

本研究于2013年在研究區(qū)農(nóng)田選擇了2塊等面積的樣地，分別位于直接受損區(qū)和間接受損區(qū)，樣地均為小麥-玉米輪作的旱地，土地理化性質(zhì)和農(nóng)田管理措施一致，植被各部分比重采用重鉻酸鉀容量法—外加熱法測定，農(nóng)田植被NPP采用“收割法”測定[40]。為了驗(yàn)證植被NPP模擬值的準(zhǔn)確性，將2013年模型模擬結(jié)果與2013年實(shí)測值進(jìn)行對比。由表1可知，模擬的結(jié)果與實(shí)測值差異較小，且植被NPP值均為直接受損區(qū)小于間接受損區(qū)，表明本研究估算的NPP基本可靠。

表1 NPP模擬值與實(shí)測值
Table 1 Simulated and measured NPP values gC/m2

樣地丘陵直接受損區(qū)丘陵間接受損區(qū)本研究258.57268.42實(shí)測值260.21269.08

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

圖7 農(nóng)田植被年均NPP與土壤有機(jī)碳空間分布Fig.7 Spatial distribution map of mean annual NPP and soil organic carbon in farmland

受損區(qū)內(nèi)農(nóng)田植被受煤炭開采、氣溫、降水量和地表人類活動(dòng)的共同影響;參照區(qū)受氣溫、降水量和地表人類活動(dòng)的影響。朱阿興等提出兩個(gè)區(qū)域的地理環(huán)境越相似，這兩個(gè)區(qū)域的目標(biāo)變量值越相近[41]，由于整個(gè)研究區(qū)面積較小，經(jīng)計(jì)算受損區(qū)、參照區(qū)的平均氣溫差值為0.2 ℃，平均降水量差值為3.1 mm，氣候條件相近;地表人類活動(dòng)主要是種植模式、農(nóng)業(yè)科技水平、施肥方式和施肥量，多年來研究區(qū)種植模式為小麥-玉米輪作，農(nóng)業(yè)科技水平相同，施肥方式和施肥量相近。因此，參照區(qū)與受損區(qū)未受損時(shí)的農(nóng)田生態(tài)系統(tǒng)固碳能力是相近的，礦區(qū)農(nóng)田固碳損失量主要受煤炭開采活動(dòng)的影響。這樣，本研究確定的參照區(qū)可很好的模擬受損區(qū)未受損時(shí)的情景。同時(shí)，煤炭開采導(dǎo)致研究區(qū)土地利用類型發(fā)生變化，而土地利用變化(農(nóng)田轉(zhuǎn)變?yōu)榉寝r(nóng)田)是對區(qū)域農(nóng)田生態(tài)系統(tǒng)固碳影響最顯著的因素，評估固碳損失量必須考慮到土地利用變化的影響。煤炭開采后部分農(nóng)田轉(zhuǎn)變?yōu)椴傻V用地和沉陷積水區(qū)，農(nóng)田面積減少。因此，煤炭開采導(dǎo)致的農(nóng)田固碳損失量應(yīng)包括農(nóng)田轉(zhuǎn)化的采礦用地、沉陷積水區(qū)的固碳量。煤炭開采導(dǎo)致的礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量計(jì)算流程:

表2 農(nóng)田生態(tài)系統(tǒng)固碳量
Table 2 Carbon sequestration in farmland ecosystems

t/km2

假設(shè)礦區(qū)農(nóng)田土地利用類型未發(fā)生變化，則土地利用現(xiàn)狀為農(nóng)田的植被固碳損失量計(jì)算公式為

Lv=(NPPr-NPPd)Sd+(NPPr-NPPi)Si

(14)

式中，NPPr,NPPd和NPPi分別為參照區(qū)、直接受損區(qū)和間接受損區(qū)的NPP均值;Sd,Si分別為直接受損區(qū)和間接受損區(qū)面積。

由于開采前土壤有機(jī)碳含量已無法獲取，因此以受損區(qū)與參照區(qū)內(nèi)土壤有機(jī)碳的差值來計(jì)算土壤有機(jī)碳的損失量。假設(shè)礦區(qū)農(nóng)田土地利用類型未發(fā)生變化，則土地利用現(xiàn)狀為農(nóng)田的土壤碳儲存損失量計(jì)算為

Lp=(SOCr-SOCd)Sd+(SOCr-SOCi)Si

(15)

式中，SOCr,SOCd和SOCi分別為參照區(qū)、直接受損區(qū)和間接受損區(qū)的農(nóng)田土壤碳密度(式(13))。

考慮到農(nóng)田土地利用變化對農(nóng)田生態(tài)系統(tǒng)固碳損失量的影響，礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量：

L=Lv+Lp+UA

(16)

式中，U為農(nóng)田土地利用類型未發(fā)生變化時(shí)受損區(qū)農(nóng)田生態(tài)系統(tǒng)單位面積固碳量(表2);A為2000—2015年農(nóng)田面積減少量，即農(nóng)田轉(zhuǎn)化為采礦用地、沉陷積水區(qū)的面積(由2000,2015年土地利用現(xiàn)狀疊加估算，2000—2015年直接受損區(qū)內(nèi)農(nóng)田面積減少約0.62 km2，間接受損區(qū)農(nóng)田面積減少約0.78 km2)。

在相近的氣候、地貌、土壤、耕作方式、農(nóng)業(yè)科技水平和施肥制度下，煤炭開采導(dǎo)致直接受損區(qū)地表沉陷，農(nóng)田植被的生長環(huán)境發(fā)生了嚴(yán)重破壞，甚至在雨季出現(xiàn)長期的積水，植被生命活動(dòng)受到抑制，植被生產(chǎn)力下降甚至徹底消失，進(jìn)而導(dǎo)致土壤有機(jī)碳含量下降，單位面積的農(nóng)田生態(tài)系統(tǒng)固碳量降低，由式(13)～(16)得，直接受損區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量為12 947.55 t，其中因農(nóng)田面積減少造成的固碳損失為6 390.29 t，土地利用現(xiàn)狀為農(nóng)田的固碳損失量為6 557.26 t;間接受損區(qū)受煤炭開采導(dǎo)致的煤煙、粉塵以及有害氣體的影響，植被生長、土壤質(zhì)量受到的負(fù)面影響相較直接受損區(qū)要小，故直接受損區(qū)固碳降低幅度大于間接受損區(qū)，由式(13)～(16)得，間接受損區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量為17 790.13 t，其中因農(nóng)田面積減少造成的固碳損失為8 401.24 t，現(xiàn)有農(nóng)田固碳損失量為9 388.89 t。綜上所述，煤炭開采造成的生態(tài)環(huán)境破壞及農(nóng)田面積減少是礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失的重要原因。

總之，與未受煤炭開采得參照區(qū)相比，在整個(gè)煤炭開采的受損區(qū)范圍內(nèi)，農(nóng)田生態(tài)系統(tǒng)固碳損失總量為30 737.68 t，其中直接受損區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量為12 947.55 t，占總損失的42.12%，間接受損區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量為17 790.13 t，占總損失的57.88%。本研究估算的煤炭開采導(dǎo)致的礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失精度相對較高，這主要是因?yàn)?以往研究忽略了間接受損區(qū)，本研究在減弱了其他因素帶來影響的同時(shí)，量化了間接受損區(qū)，并評估了直接受損區(qū)和間接受損區(qū)的固碳損失量，使礦區(qū)生態(tài)系統(tǒng)固碳損失評估更為精確;煤炭開采會(huì)導(dǎo)致原始農(nóng)田發(fā)生土地利用類型上的改變，本研究考慮了土地利用類型變化對農(nóng)田植被碳庫的影響，保證了煤炭開采對農(nóng)田生態(tài)系統(tǒng)固碳損失量評估的精度。煤炭開采對區(qū)域農(nóng)田生態(tài)系統(tǒng)固碳量的影響屬于失碳效應(yīng)。

4.4 方法驗(yàn)證

為了驗(yàn)證本研究方法的適用性:① 采用參照區(qū)M2,H2來估算間接受損范圍及農(nóng)田生態(tài)系統(tǒng)固碳損失量，經(jīng)計(jì)算，煤炭開采對直接受損區(qū)周邊環(huán)境的影響范圍為山地990 m、丘陵660 m，影響范圍與上文結(jié)果一致;煤炭開采導(dǎo)致的農(nóng)田生態(tài)系統(tǒng)固碳損失量為29 274.57 t，與上文結(jié)果僅相差4.76%，即使同一區(qū)域的植物自身生產(chǎn)力也會(huì)存在一定的差異，故不同參照區(qū)計(jì)算農(nóng)田生態(tài)系統(tǒng)固碳損失量也會(huì)帶來小范圍的誤差。② 直接受損區(qū)位于長河流域的8個(gè)煤礦區(qū)中，將8個(gè)煤礦的直接受損區(qū)、間接受損區(qū)農(nóng)田生態(tài)系統(tǒng)單位面積固碳損失量與研究區(qū)進(jìn)行對比，相對偏差分布如圖8所示。根據(jù)礦區(qū)資料顯示天安岳圣山、天泰坤達(dá)煤礦進(jìn)行了復(fù)墾，單位面積固碳損失量較小;天安盈盛開采規(guī)模擴(kuò)大，生態(tài)破壞加劇，單位面積固碳損失較大，其余煤礦區(qū)與研究區(qū)整體單位面積固碳損失量的相對偏差較小，表明該方法評估煤礦開采導(dǎo)致的礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失量是適用的。

圖8 煤礦區(qū)農(nóng)田固碳損失量Fig.8 Carbon sequestration loss of farmland in coal mine area

5 討 論

煤炭開采活動(dòng)損壞了原生地貌和地表植被，造成土地利用類型變化，引起植被-土壤碳匯量減少，導(dǎo)致了農(nóng)田生態(tài)系碳儲量的大量流失[12,17,42]。當(dāng)前，有學(xué)者采用多年礦區(qū)植被NPP的差值來估算煤炭開采導(dǎo)致的生態(tài)損失[13]，但由于各年氣溫、降水等環(huán)境因素的變化較大，不能較好的說明礦區(qū)生態(tài)損失由煤炭開采所引起;也有學(xué)者采用沉陷區(qū)周邊區(qū)域作為參照區(qū)來評估沉陷區(qū)植被固碳或土壤碳儲存的損失[14,43-44]，雖然相同年份內(nèi)沉陷區(qū)與參照區(qū)氣溫、降水等環(huán)境因素基本一致，然而參照區(qū)本身已受到煤炭開采的間接影響，使得評估的結(jié)果偏低。本研究在基于傳統(tǒng)的開采沉陷、采空及損毀土地確定礦區(qū)受損范圍的基礎(chǔ)上，提出了測算周邊NDVI指數(shù)變化趨勢以量化間接受損區(qū)的方法。該方法可適用于大范圍、長時(shí)間序列的地表植被覆蓋區(qū)，可更準(zhǔn)確的確定煤炭開采對地表生態(tài)的影響邊界。同時(shí)，本研究考慮了土地利用類型變化對農(nóng)田生態(tài)系統(tǒng)碳庫的影響，提高了礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失的評估精度，也為研究礦區(qū)時(shí)空生態(tài)系統(tǒng)固碳累積變化及后續(xù)生態(tài)補(bǔ)償機(jī)制等問題提供借鑒。

已有研究表明，合理的土地復(fù)墾及差別整治有利于礦區(qū)生態(tài)系統(tǒng)碳匯量增加，一定程度上可抑制煤炭開采造成的生態(tài)環(huán)境破壞[10,45]。利用遙感影像來量化礦區(qū)生態(tài)受損范圍，可精確估算礦區(qū)固碳損失量，動(dòng)態(tài)監(jiān)測礦區(qū)生態(tài)環(huán)境及未來土地復(fù)墾的生態(tài)效益與碳匯量變化。因此，結(jié)合礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失程度，科學(xué)合理的進(jìn)行地貌重塑、土壤重構(gòu)、植被重建及生物多樣性恢復(fù)等，對于提升農(nóng)田生產(chǎn)能力、重建農(nóng)田生態(tài)系統(tǒng)、緩解氣候變暖具有重要意義。本研究對長河流域2000—2015年的農(nóng)田生態(tài)系統(tǒng)直接受損區(qū)和間接受損區(qū)的固碳損失量進(jìn)行了估算，奠定了礦區(qū)碳增匯經(jīng)濟(jì)潛力評價(jià)的基礎(chǔ)，也為礦區(qū)土地復(fù)墾、生態(tài)修復(fù)與重建提供更加準(zhǔn)確和實(shí)際的指導(dǎo)。然而，本研究提出的評估方法僅適用于中國華北丘陵山區(qū)，針對其他區(qū)域的實(shí)施需依據(jù)實(shí)際情況進(jìn)行調(diào)整。同時(shí)，在今后的研究中應(yīng)建立長期的土壤養(yǎng)分監(jiān)測站點(diǎn)，以確定煤炭開采情況下農(nóng)田土壤碳儲存的動(dòng)態(tài)變化及礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳損失的機(jī)理機(jī)制。

6 結(jié) 論

(1)研究區(qū)內(nèi)由于煤炭開采導(dǎo)致的農(nóng)田直接受損面積為5.06 km2，間接受損面積為18.99 km2。在間接受損區(qū)內(nèi)，生態(tài)環(huán)境的損毀程度隨著與直接受損區(qū)距離的增大而逐漸減弱。

(2)礦區(qū)農(nóng)田生態(tài)系統(tǒng)植被固碳量為3 936.16 t/km2，小于土壤碳儲存量的6 737.05 t/km2。

(3)研究區(qū)因煤炭開采導(dǎo)致的農(nóng)田生態(tài)系統(tǒng)固碳損失量為30 737.68 t。其中，直接受損區(qū)固碳損失量為12 947.55 t，間接受損區(qū)固碳損失量為17 790.13 t。煤炭開采對礦區(qū)農(nóng)田生態(tài)系統(tǒng)固碳量的影響屬于失碳效應(yīng)。

(4)煤炭開采導(dǎo)致的生態(tài)系統(tǒng)固碳損失評估及生態(tài)環(huán)境損毀的研究中，間接受損區(qū)也應(yīng)該納入評估范圍。