武翼飛,陳國光

(1.中國地質(zhì)科學院，北京 100037；2.中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心,江蘇 南京 210016;3.中國地質(zhì)大學(北京),北京 100083)

土壤保持功能是生態(tài)系統(tǒng)服務功能的一個重要方面,為土壤形成、植被固著、水源涵養(yǎng)等提供了重要基礎[1]。生態(tài)系統(tǒng)服務和權(quán)衡綜合評估模型(InVEST模型),是目前應用最廣泛的生態(tài)系統(tǒng)服務功能評估模型[2]。InVEST模型中的土壤保持模塊(SDR)在通用土壤流失方程(RUSLE)的基礎上進行改進,考慮到地塊自身攔截上游沉積物的能力,彌補了RUSLE方程中忽略這一水文過程的缺陷,使得土壤保持功能評估結(jié)果更加直觀、合理[3]。近些年來,InVEST模型在我國不同地區(qū)的生態(tài)系統(tǒng)服務功能評估和土壤保持研究領(lǐng)域得到了應用。周彬等[4]應用InVEST模型對北京山區(qū)土壤侵蝕進行模擬,探究不同森林類型土壤保持功能差異;劉英等[5]基于InVEST模型構(gòu)建綜合生態(tài)系統(tǒng)服務評估模型,對神東礦區(qū)生態(tài)系統(tǒng)服務功能進行分析。隨著InVEST模型的廣泛應用,唐堯等[6]、肖微[7]等發(fā)現(xiàn)模型需根據(jù)地區(qū)具體情況不斷調(diào)整參數(shù)等不足之處。

福建省大田縣銀頂格礦區(qū)曾經(jīng)是福建省主要鐵礦產(chǎn)區(qū),礦區(qū)開采歷史久,尤其是近30年運用大型機械深度開采,導致礦山開采區(qū)域山體植被破壞嚴重、生物多樣性降低;礦山邊坡土層結(jié)構(gòu)松散、透水強、土壤肥力差,且土壤和水體污染嚴重,水源涵養(yǎng)能力不足[8]。近些年，隨著福建省閩江流域山水林田湖草生態(tài)保護項目的推進,銀頂格廢棄礦山成為重點生態(tài)修復對象。本文以銀頂格礦區(qū)所處小流域為研究對象,以土壤保持功能作為評估指標,基于InVEST模型應用銀頂格礦區(qū)的降雨侵蝕力因子、土壤可蝕性因子、植被覆蓋度、DEM和生物物理系數(shù)(如土壤保持措施因子)等數(shù)據(jù),對比分析生態(tài)修復前后土壤保持功能的變化,探討土壤保持功能的影響因素,為礦區(qū)生態(tài)修復效果評估和礦區(qū)生態(tài)系統(tǒng)服務功能評估提供借鑒。

1 研究區(qū)概況

研究區(qū)位于福建省大田縣前坪鄉(xiāng)、太華鎮(zhèn)、均溪鎮(zhèn)交界處(117°48′38″～117°50′34″E、25°50′0″～25°51′28″N)的銀頂格礦區(qū)及所處的上華村流域,總面積為31.12 km2(圖1)。區(qū)內(nèi)地質(zhì)構(gòu)造復雜,巖石以鈣質(zhì)砂巖、石英砂巖為主,其次為花崗巖、凝灰?guī)r和變粒巖[9]。大田縣銀頂格礦區(qū)生態(tài)修復項目于2019年10月啟動,2021年4月完成修復。生態(tài)修復工程主要采取工程措施與生物措施相結(jié)合的方式。工程措施以水土保持、土地恢復為主,采取削坡減載措施,將采礦區(qū)挖成多級平臺,共挖土方約729 062 m3;對原來剝離挖損的土地采取回填推平措施,共回填推平土方約686 155 m3;各平臺設置鋼筋混凝土截排洪溝合計長約39 552 m; 對強酸性土壤采取“播撒石灰石+覆蓋專用防水毯”治理措施。生物措施以改良土壤、景觀恢復為主,采取挖穴種植、播種草籽、栽植草坪等措施,恢復植被綠化面積約2.29 km2。在修復工程完成后,礦山景觀得到改善,邊坡穩(wěn)定性提高,坡面流水通道暢通,存在的地質(zhì)災害隱患基本消除,酸性土壤得到有效治理,共恢復耕地面積0.21 km2,增加建設用地0.32 km2和農(nóng)業(yè)設施用地0.94 km2。

圖1 福建省大田縣(a)及研究區(qū)(b)地理位置圖Fig. 1 Geographical location map of Datian County (a) and the study area (b) in Fujian Province

2 研究方法

2.1 InVEST模型

本文利用InVEST模型中的土壤保持功能模塊(SDR)評估銀頂格礦區(qū)的土壤保持能力,通過對比修復工程前后的土壤保持狀況,分析礦區(qū)土壤保持功能恢復效果。SDR模塊又稱為泥沙輸移比模塊,計算公式為

RKLSx=R·Kx·LSx，

(1)

USLEx=R·Kx·LSx·Cx·Px，

(2)

(3)

SDRx=SEDRx+RKLSx-USLEx，

(4)

式中：SDRx是柵格x的土壤保持量，SEDRx為泥沙持留量，RKLSx為潛在土壤侵蝕量，USLEx和USLEy分別為柵格x和其上坡柵格y的實際土壤侵蝕量，指植被覆蓋下土壤流失量，以上單位均為t；SE為泥沙持留效率；R是降雨侵蝕力因子，MJ·mm·hm-2·h-1·a-1；K是土壤可蝕性因子，t·hm2·h·hm-2·MJ-1·mm-1；LS代表坡長與坡度因子；C為植被覆蓋管理因子；P為土壤保持措施因子。LS因子、C因子和P因子均是無量綱參數(shù)。

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

為評估銀頂格礦區(qū)的土壤保持功能,本文采用遙感調(diào)查和資料收集相結(jié)合的方法開展工作。土地利用數(shù)據(jù)和植被覆蓋數(shù)據(jù)是利用研究區(qū)2019年11月7日和2021年12月7日兩期空間分辨率為2 m的高分六號遙感影像數(shù)據(jù),經(jīng)過預處理、土地利用類型解譯[10-11]、實地驗證后獲得。InVEST模型的坡度與坡長因子分別選取兩期DEM數(shù)據(jù):第一期選用從NASA網(wǎng)站獲取的2019年12.5 m空間分辨率DEM數(shù)據(jù);第二期是2021年11月在研究區(qū)利用大疆無人機獲得的礦山重點修復地區(qū)數(shù)字高程模型。土壤類型數(shù)據(jù)來源于福建省主要土壤類型專題圖,由福建省1∶250 000土壤類型分布圖[12]裁剪而成;土壤粒徑、有機碳數(shù)據(jù)來源于中國土種數(shù)據(jù)庫[13]。氣象數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn),本文選取福建省19個氣象站點在2019年6月—2021年5月逐月降雨量數(shù)據(jù),通過ArcGIS克里金插值獲得研究區(qū)像元大小為12.5 m×12.5 m的月均降雨量柵格數(shù)據(jù),用于計算降雨侵蝕力因子。

2.3 各因子計算方法

2.3.1 降雨侵蝕力因子

降雨侵蝕力因子(R)是降雨引起土壤侵蝕的潛在能力,是導致土壤侵蝕的關(guān)鍵性因素之一[14],主要受區(qū)域降雨量和降雨強度的影響。周伏健等[15]提出一套適用于福建省的基于月降雨量的降雨侵蝕力估算模型,該模型在降雨較充沛地區(qū)得到了廣泛應用。本文采用此模型,利用站點月降雨量數(shù)據(jù),估算銀頂格礦區(qū)的降雨侵蝕力因子，其公式為

(5)

式中:R為降雨侵蝕力因子,MJ·mm·km-2·h-1·a-1;Pi為第i個月降雨量,mm。

經(jīng)計算，研究區(qū)內(nèi)降雨侵蝕力因子最大值為343 822 MJ·mm·km-2·h-1·a-1,最小值為331 897 MJ·mm·km-2·h-1·a-1,整體相差不大。利用ArcGIS軟件克里金空間插值生成2019年和2021年研究區(qū)降雨侵蝕力因子分布圖(圖2),與陳曉瑜等[16]關(guān)于福建省降雨侵蝕力的研究結(jié)論相符。

圖2 研究區(qū)2019年(a)和2021年(b)降雨侵蝕力因子分布圖Fig. 2 Rainfall erosivity factor map of the study area in 2019 (a) and 2021 (b)

2.3.2 土壤可蝕性因子

土壤可蝕性因子(Kx)反映土壤對剝蝕與搬運的敏感性,不同土壤類型導致土壤侵蝕差異,其數(shù)值越高表示土壤越容易被侵蝕[17]。目前計算土壤可蝕性因子最為普遍是運用EPIC模型,方綱等[18]根據(jù)此模型計算出福建省主要土壤類型的土壤可蝕性因子,并應用到福建省水土流失研究中。EPIC模型中Kx的計算公式為

(6)

式中:Sa代表砂粒(粒徑0.05～2 mm)含量,%;Si代表粉粒(粒徑0.002 ～0.05 mm)含量,%;Sc代表黏粒(粒徑<0.002 mm)含量,%;C代表有機碳含量,%。

研究區(qū)內(nèi)的土壤類型主要包括黃壤、紅壤、水稻土[19],其中紅壤以黃紅壤、硅鈣質(zhì)紅壤、硅鋁鐵質(zhì)紅壤為主,水稻土以滲育水稻土和潴育水稻土為主(圖3)。在礦區(qū)生態(tài)修復工程中主要采用了削坡平整的工程措施,這對研究區(qū)土壤類型的影響有限,因此，修復前后使用同一組土壤可蝕性因子。本文從中國土種數(shù)據(jù)庫中選取鄰近地區(qū)成土母質(zhì)的土壤典型剖面的土壤粒徑數(shù)據(jù)和有機碳數(shù)據(jù)計算土壤可蝕性因子,如黃紅壤選擇位于武夷山市星村的石英砂巖成土母質(zhì)的典型剖面,紅壤選擇位于明溪縣花崗巖成土母質(zhì)典型剖面。經(jīng)公式(6)計算后,得到研究區(qū)土壤類型典型剖面的土壤可蝕性因子(表1),在土壤類型圖中賦值生成像元12.5 m×12.5 m的柵格數(shù)據(jù),最終得到研究區(qū)土壤可蝕性因子分布圖(圖3)。

表1 研究區(qū)土壤類型可蝕性因子Table 1 Soil erodibility factor of main soil types in the study area

2.3.3 坡長與坡度因子

坡長與坡度因子(LSx)表示地形地貌因素對土壤侵蝕的影響程度,是影響土壤侵蝕的重要因素,InVEST模型一般通過DEM數(shù)據(jù)計算坡度與坡長因子[20]。

運用ArcGIS軟件對修復前后的DEM數(shù)據(jù)(圖4)進行填洼處理和地形反轉(zhuǎn)處理,分別對正反地形DEM提取水流方向、計算匯流累積量、設置匯流累積量閾值(Fa)、生成河網(wǎng)、生成集水區(qū)。把集水區(qū)與反向地形集水區(qū)融合,得到匯水線與分水嶺所組成的區(qū)域即為斜坡單元。應用坡度工具提取研究區(qū)坡度,研究區(qū)內(nèi)坡度范圍為0°～66°。依據(jù)《SL277—2002水土保持監(jiān)測技術(shù)規(guī)程》[21]劃分方案,將研究區(qū)坡度劃分為6個等級:微坡(0°～5°)、較緩坡(5°～8°)、緩坡(8°～15°)、較陡坡(15°～25°)、陡坡(25°～35°)、急陡坡(>35°)。

圖4 2019年(a)和2021年(b)重點修復地區(qū)高程圖Fig. 4 DEM map of key restoration areas in 2019 (a) and 2021 (b)

InVSET模型可根據(jù)不同坡度自動選用不同的公式進行計算,由于模型是以美國坡度較小的試驗區(qū)為原型,因此需調(diào)整邊坡閾值,即大于此坡度的地區(qū)應當在采取水土保持措施條件下進行農(nóng)業(yè)活動。根據(jù)王敏等[22]相關(guān)研究,坡度受土壤侵蝕影響的臨界坡度為25°,因此本文將邊坡閾值設為25°。

當坡度<邊坡閾值時,坡度坡長因子的計算方法公式為

(7)

式中:LSx為坡度坡長因子;Fa為匯流累積量閾值;Cs為柵格大小,本文DEM數(shù)據(jù)柵格大小為12.5 m;s為百分比坡度,%;n為坡長指數(shù),取值參考公式(8),

(8)

當坡度>邊坡閾值時,坡度坡長因子的計算方法為

LS= 0.08β0.35s0.6，

(9)

(10)

式中:β為根據(jù)流向確定的柵格大小,m;s為百分比坡度,%。

2.3.4 植被覆蓋管理因子

植被覆蓋管理因子(C)是植被覆蓋對土壤侵蝕的抑制能力,與植被覆蓋度、土地利用類型相關(guān)。C值介于0～1,數(shù)值越大說明該區(qū)域表面的植被覆蓋度越低,土壤侵蝕越嚴重[23]。應用EVNI軟件提取出修復前后的歸一化植被指數(shù)(NDVI)[24],再通過NDVI估算出研究區(qū)植被覆蓋度。植被覆蓋管理因子的計算采用目前適用較廣的蔡崇法模型[25],公式為

(11)

式中:C為植被覆蓋管理因子,取值范圍為[0,1];FVC為植被覆蓋度,%。

在ArcGIS軟件運用上述公式計算得到研究區(qū)修復前后的植被覆蓋管理因子柵格圖。為符合InVEST模型數(shù)據(jù)輸入格式,需要以土地利用類型為

單位確定植被覆蓋管理因子的平均值。利用ArcGIS軟件的空間分析工具統(tǒng)計出各個土地利用類型的植被覆蓋管理因子平均值,其中將水面和交通運輸用地類型分別修正為0和1,最后將因子賦值到研究區(qū)土地利用類型圖中,得到研究區(qū)2019年和2021年植被覆蓋管理因子分布圖(圖5)。

圖5 研究區(qū)2019年(a)和2021年(b)植被管理覆蓋因子分布圖Fig. 5 Vegetation cover-management factor maps of the study area in 2019 (a) and 2021 (b)

2.3.5 土壤保持措施因子

土壤保持措施因子(Px)是反映采取水土保持措施后削弱土壤侵蝕能力的因子,取值范圍為[0,1]。取值越小則代表該地區(qū)采取的水土保持措施對土壤侵蝕的抑制越有效,數(shù)值為0表示該地區(qū)無侵蝕;取值越高則代表水土保持措施越差,數(shù)值為1表示該地區(qū)未采取任何土壤保持措施?，F(xiàn)今常采用土地利用類型賦值的方法[26],同一種土地利用類型采用的土壤保持措施大致相同。因此，本文生態(tài)修復前后同一種土地利用類型的P值也相同。借鑒周來等[27]、楊冉冉等[28]相關(guān)研究成果,確定研究區(qū)的土壤保持措施因子(表2):河流、坑塘、溝渠等水面賦值為0,該土地利用類型不會發(fā)生土壤侵蝕;交通用地、農(nóng)村宅基地、工業(yè)用地等建設類用地,由于地面硬化幾乎不發(fā)生土壤侵蝕,因此，賦予有效位數(shù)的最小值0.001;生態(tài)修復前的采礦用地P值賦值為1;其余植被覆蓋的用地類型根據(jù)植物的抗侵蝕能力賦予0～1的值。

表2 不同土地利用類型的土壤保持措施因子表Table 2 The support practice factors of different land use types

2.3.6 模型參數(shù)及運行

植被覆蓋管理因子和土壤保持措施因子需要制作成生物物理參數(shù)表格(csv格式),表中對相應的土地利用類型進行編碼。矯正參數(shù)K、IC0和SDRmax值分別采用InVEST模型系統(tǒng)[22]默認值2、0.5、0.8。打開InVEST模型中的土壤保持模塊(SDR),將所需數(shù)據(jù)導入運行模型,得到運算結(jié)果、柵格數(shù)據(jù)和矢量數(shù)據(jù)。

3 結(jié)果與分析

3.1 土壤保持功能分析

運用InVEST模型的土壤保持模塊獲得研究區(qū)修復前后的潛在土壤侵蝕量、實際土壤侵蝕量、泥沙持留量和土壤保持量等數(shù)據(jù)。

3.1.1 土壤侵蝕量空間分布及時空變化

通過統(tǒng)計得出,2019年銀頂格礦區(qū)的實際土壤侵蝕總量為149 599.7 t,平均土壤侵蝕模數(shù)為4 806.4 t·km-2·a-1;2021年銀頂格礦區(qū)的實際土壤侵蝕總量為81 373.3 t,平均土壤侵蝕模數(shù)為2 614.4 t·km-2·a-1?？偭可蟻砜?2021年的實際土壤侵蝕量比2019年減少了68 226 t,平均土壤侵蝕模數(shù)降低了2 192.1 t·km-2·a-1。根據(jù)中華人民共和國水利部頒布的《SL190—2007土壤侵蝕分類標準》[29]劃分方案,土壤侵蝕強度劃分為微度侵蝕、輕度侵蝕、中度侵蝕、強烈侵蝕、極強烈侵蝕、劇烈侵蝕6種類型(表3),依據(jù)此方案得到2019年和2021年銀頂格礦區(qū)土壤侵蝕強度分級圖(圖6)。對比發(fā)現(xiàn)，2021年中度及以上侵蝕強度地區(qū)面積明顯減少,強烈及以上侵蝕強度地區(qū)面積占比從2019年的21.0%降低至2021年的8.7%,表明生態(tài)修復工程對土壤侵蝕治理效果明顯。

表3 不同土壤侵蝕強度統(tǒng)計結(jié)果Table 3 Statistical results of different soil erosion intensity

圖6 研究區(qū)2019(a)和2021年(b)土壤侵蝕強度圖Fig. 6 Soil erosion intensity map of the study area in 2019 (a) and 2021 (b)

經(jīng)計算,修復后仍有部分地區(qū)土壤侵蝕較嚴重。實地調(diào)查后發(fā)現(xiàn),目前仍有部分修復地區(qū)尚未被利用或裸露,以及部分強酸土壤治理區(qū)的植被恢復進度較慢,這些地區(qū)植被覆蓋情況尚未改善,加之2021年降雨侵蝕力因子相對降低,都對修復結(jié)果產(chǎn)生一定影響。因此，生態(tài)修復后還應繼續(xù)加強植物管護,推進礦區(qū)土地資源高效、合理利用。

3.1.2 土壤保持量空間分布及時空變化

應用ArcGIS軟件得到2019年和2021年銀頂格礦區(qū)土壤保持強度圖(圖7)。對運算結(jié)果進行統(tǒng)計分析表明,2019年銀頂格礦區(qū)的土壤保持總量為17 546 658.6 t,土壤保持強度為563 748.2 t·km-2·a-1;2021年銀頂格礦區(qū)的土壤保持總量為17 822 085.7 t,土壤保持強度為572 597.1 t·km-2·a-1。相較于2019年,2021年土壤保持總量和土壤保持強度都明顯增加,土壤保持總量增加了275 427.1 t,土壤保持強度增加了8 848.9 t·km-2·a-1。尤其是研究區(qū)北部的礦山修復重點區(qū)在礦山修復工程的作用下,植被覆蓋增加、土壤泥沙截留能力提升、土壤保持強度明顯增強,表明生態(tài)修復工程對提升整體土壤保持能力效果較好。

圖7 研究區(qū)2019年(a)和2021年(b)土壤保持強度圖Fig. 7 Soil conservation intensity map of the study area in 2019 (a) and 2021 (b)

3.2 土壤保持功能的影響因素

3.2.1 土壤保持與土地利用類型的關(guān)系

由于不同土地利用類型的土壤保持能力存在差異,本文將2021年研究區(qū)的土壤保持量依據(jù)不同土地利用類型進行分區(qū)統(tǒng)計(表4)。統(tǒng)計結(jié)果顯示:土壤保持總量從大到小依次是喬木林地、其他林地、水田、其他草地、采礦用地、旱地、竹林地、農(nóng)村宅基地、交通運輸用地、園地、設施用地;土壤保持強度從高到低依次是水田、其他林地、喬木林地、其他草地、采礦用地、竹林地、旱地、交通運輸用地、農(nóng)村宅基地、園地、設施用地。其中,喬木林地、其他林地的土壤保持總量最大,共占研究區(qū)土壤保持總量的77.1%,表明植被因子對土壤保持具有關(guān)鍵性作用。在今后的土壤保持防治工作中,應加大對林草地的保護力度,采取適當?shù)乃帘３执胧﹣碓黾悠渫寥辣３帜芰?。由于水田采取了較好的土壤保持措施,且在地形較平緩地區(qū),有效減弱了土壤侵蝕現(xiàn)象。研究區(qū)園地的植被覆蓋相對較低,抗侵蝕能力較低,故土壤保持強度差。設施用地的土壤保持強度最差,主要因為研究區(qū)大部分的設施用地是由原來的采礦用地平整而來,相關(guān)設施還未完成建設,大片土地裸露,水土保持能力差。

表4 不同土地利用類型的土壤保持量統(tǒng)計結(jié)果Table 4 Statistical results of soil conservation in different land use types

3.2.2 土壤保持與坡度的關(guān)系

將研究區(qū)土壤保持量分布圖與DEM圖進行對比,發(fā)現(xiàn)土壤保持能力與坡度具有相關(guān)關(guān)系。將2021年坡度分級數(shù)據(jù)與土壤保持數(shù)據(jù)、植被覆蓋度進行分區(qū)統(tǒng)計,得出不同坡度下的植被覆蓋度、土壤保持強度和土壤保持總量(表5)。結(jié)果表明,坡度為0°～5°的地區(qū)土壤保持強度最弱,植被覆蓋度也最低;從較緩坡至陡坡,植被覆蓋度提高5%,土壤保持強度增加59 269 t·km-2·a-1。進一步統(tǒng)計發(fā)現(xiàn),隨著坡度的升高,植被覆蓋度也相應提高,土壤保持強度也相應提高。經(jīng)過分析和實地驗證發(fā)現(xiàn),礦區(qū)內(nèi)此前受到嚴重破壞的高坡度地區(qū)已基本得到治理,現(xiàn)有的高坡度地區(qū)植被未受破壞,從而保持了較高的土壤保持強度。這也反映出,在使用土壤保持模塊(SDR)時考慮到地塊自身攔截上游沉積物能力的情況下,可以計算出更為真實的土壤保持強度和保持總量。

表5 不同坡度等級的土壤保持總量Table 5 Soil conservation for different slope gradients

不同坡度等級的土壤保持強度結(jié)果表明,在植被覆蓋保持較高的情況下,泥沙的截留能力和植物根系對土壤的固結(jié)能力也會逐漸增強,從而抵消坡度對土壤保持的不利影響。因此，在銀頂格地區(qū)水土保持工作中,要加強植被保護,提高土壤泥沙截留能力,以增強土壤保持功能。

4 結(jié)論

(1)銀頂格礦區(qū)修復后,研究區(qū)土壤保持強度增加了8 848.9 t·km-2·a-1,土壤保持總量增加了275 427.1 t,表明研究區(qū)土壤保持強度明顯提高，土壤保持總量有所提升。土壤侵蝕強烈及以上強度地區(qū)面積明顯減少,占比從21.0%降低至8.7%,表明生態(tài)修復工程提升了研究區(qū)土壤保持能力。

(2)不同的土地利用類型中,水田、林地的土壤保持強度較強,有植被覆蓋的土地類型的土壤保持能力遠大于人類開發(fā)地區(qū)。從較緩坡至陡坡,植被覆蓋度提高5%,土壤保持強度增加59 269 t·km-2·a-1,隨著坡度的升高及植被覆蓋度提高,泥沙截留能力和植物根系對土壤的固結(jié)能力也會逐漸增強,表明在地表植被不受破壞、保持較高植被覆蓋的情況下,坡度對土壤保持強度的影響逐漸減弱。在土壤侵蝕防治工作中可以通過改善土地利用類型、加大對林草地的保護力度、提高植被覆蓋度、加強陡坡生態(tài)植被保護等措施來提升土壤保持能力。

(3)InVEST模型的土壤保持模塊(SDR)考慮到地塊自身攔截上游沉積物的能力,彌補了RUSLE方程的缺陷,使得計算結(jié)果更加真實、準確。模型評估結(jié)果為柵格格式,實現(xiàn)了結(jié)果的可視化,可以更直觀、便利地反映研究區(qū)土壤侵蝕和土壤保持情況。