馮 斌, 張學伍, 徐 敏, 胡海峰

〔中煤科工重慶設計研究院(集團)有限公司, 重慶 400016〕

土壤侵蝕與泥沙輸移是氣候、土壤、植被、水文、地形等多種因素的綜合結果[1]，代表著陸地生態(tài)系統(tǒng)中重要的物質循環(huán)與輸移過程[2]。流域內土地利用景觀類型、數量通過改變土壤侵蝕的發(fā)生機制，改變水文結構和侵蝕系統(tǒng)，引起土壤輸移、攔截能力的變化，進而影響產沙量的增加或減少[3-4]。隨著經濟活動的發(fā)展，由于人類活動導致的土地利用變化是影響土壤侵蝕的直接驅動力，也是導致水土流失時空分異的重要原因。除了土地利用變化的影響外，流域泥沙輸移還受泥沙輸移路徑影響，泥沙輸移和連通性關系密切，水文研究領域借用連通性的概念開展了坡面—溝道水文連通、河—湖和流域水文連通的研究[5-7]。與泥沙輸移直接相關的泥沙連通性，是指流域內泥沙通過分離和輸移過程中從源到匯的傳輸程度[8]。泥沙連通性是可以很好地表示流域內泥沙輸移潛力的表征量。Borselli等人[9]開發(fā)了分布式的泥沙連通性指數(IC)，量化泥沙連通的潛在大小。Cavalli等人[10]對IC指數的坡度、泥沙貢獻面積及權重因子的方法進行了改進。以連通性為基礎，可以開發(fā)出不同空間位置的泥沙輸移比[11]，被應用于評價泥沙輸移量。流域泥沙連通性主要由流域異質性的空間格局所決定[12]，對景觀空間格局的改變均會影響泥沙連通性，進而影響泥沙輸移過程，如修建梯田、淤地壩、道路及建設用地增加等[13-14]。總體來說，在考慮流域泥沙輸移通量時，現有研究主要關注的是土地利用變化對土壤侵蝕的影響，而忽視了因土地利用變化導致的泥沙連通性的影響。在國外，隨著過去20 a間頁巖氣開發(fā)的快速發(fā)展，頁巖氣開發(fā)平臺、道路和管線占用土地和干擾環(huán)境狀況逐漸受到關注[15]，頁巖區(qū)開采需要修建大量的通行道路，進一步增加建設用地，可能會導致區(qū)域土地利用變化明顯。隨著中國頁巖氣將大規(guī)模開發(fā)利用，相關學者[16]建議頁巖氣開發(fā)設計應考慮占地、景觀破碎化的影響。本研究在考慮頁巖氣開采區(qū)建設對于土地利用的改變的影響基礎上，探尋這種景觀變化對對泥沙流失影響。本文以重慶市涪陵區(qū)中國首批頁巖氣勘查開發(fā)示范區(qū)為研究區(qū)域，分析開發(fā)前后景觀格局的變化(2010—2019年)，并利用RUSLE通用土壤流失方程計算土地利用變化導致的土壤侵蝕變化，分析土地利用對泥沙連通性和泥沙輸移的影響，研究結果可為類似區(qū)域控制土壤流失為目標的流域土地利用布局和優(yōu)化提供重要的參考價值。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

研究區(qū)位于中國重慶市涪陵區(qū)焦石鎮(zhèn)，地理坐標為：東經117°26′07′—117°38′13″，北緯29°33′50″—29°47′32″，處于四川盆地及山地過渡地帶，面積約為284 km2，地形以低山丘陵為主，研究區(qū)范圍屬于頁巖氣開采區(qū)的范圍，來源于重慶市規(guī)劃與自然資源局。研究區(qū)受地質構造和長江、烏江水系的切割，相對高差較大，為141～1 422 m。該區(qū)為亞熱帶濕潤季風氣候區(qū)，雨量充沛，區(qū)域年平均降雨1 100 mm，但分布不均，年內降水主要集中于5—9月，占全年降水量的68%左右，年均蒸發(fā)量1 076 mm，相對濕度81%；多年平均氣溫18.0 ℃，≥10 ℃積溫5 719 ℃，無霜期多年平均317 d，日照時數1 087 h；區(qū)內主導風向為東北向，年平均風速0.7 m/s。研究區(qū)土壤以紫色土和水稻土為主；地帶性植被類型屬亞熱帶常綠闊葉林，由于人為原因區(qū)域內的天然林已大量消失，目前植被主要為人工林、天然次生林及農作物。

1.2 研究方法

1.2.1 數據來源 研究區(qū)的頁巖氣開發(fā)始于2012年，到2018年研究區(qū)的頁巖氣開發(fā)的附屬設施相關建設已基本完成，因此選擇2010—2019年為研究時間段。土壤侵蝕模擬數據主要包括日降雨、地形、土地利用、土壤理化參數等。日降雨數據來源于涪陵區(qū)氣象局，2010—2019年期間日降雨數據，主要用于計算降雨侵蝕力。地形數據是空間柵格大小12.5 m的DEM數據，該數據由ALOS(advanced land observing satellite)衛(wèi)星相控陣型L波段合成孔徑雷達(PALSAR)采集，數據下載網址https:∥urs.earthdata.nasa.gov/users/new，主要用于計算侵蝕模擬中的坡長和坡度因子。土地利用來源于Google earth下載的遙感影像(2010和2019年)，根據中國的《土地利用現狀分類》標準和區(qū)域的土地利用特點，將土地利用劃分為水田、旱耕地、林地、草地、建筑用地(頁巖氣開采地、道路)。土壤可蝕性的計算依靠有機質和砂、粉、黏3種粒徑顆粒的百分比含量，相關數據來源于國家地球系統(tǒng)科學數據中心。本研究利用徑流小區(qū)的泥沙流失量數據來驗證本文的計算結果，徑流小區(qū)的數據來源于2019年重慶市水土保持公報，從重慶市水利局網站公開下載，本研究共采用了公報中的16個徑流小區(qū)數據，這些小區(qū)泥沙流失量數據觀測年份是2018年部分降雨事件結果(表1)。

表1 研究區(qū)徑流小區(qū)數據

1.2.2 土壤流失的模擬 土壤流失模擬主要包括土壤侵蝕和泥沙輸移比的模擬。土壤侵蝕的模擬采用修正的通用土壤流失方程(RUSLE)計算空間尺度的土壤侵蝕強度，泥沙輸移比基于土壤連通性指數來反映泥沙的輸移。具體計算公式[17]為：

SSYi=Ei·SDRi

(1)

式中：SSYi指柵格i的產沙模數〔t/(hm2·a)〕;Ei柵格i的土壤侵蝕量〔t/(hm2·a)〕; SDRi是柵格i的泥沙輸移比。其中土壤侵蝕量的估算公式為:

Ei=R·Ki·LSi·Ci·Pi

(2)

式中：R是降雨侵蝕力〔(MJ·mm)/(hm2·h·a)〕；K是土壤侵蝕能力〔(t·h)/(MJ·mm)〕；LSi是地形因素；C是植被因素；P是水土保持工程措施因素。

(1) 降雨侵蝕力(R)。降雨侵蝕力是指一次降雨總動能(E)與該次降雨的最大30 min雨強(I30)的乘積，亦稱EI30指數。但是由于這類詳細的降雨數據獲取困難，在研究中通?；谌战涤瓿叨葦祿?，把日降雨量大于12 mm的降雨稱為侵蝕性降雨[18]，本文利用章文波等[22]提出日降雨侵蝕力計算方法，該方法應用在中國的第一次水土流失調查中，在中國有較為應用的普適性。

(3)

式中：k=1,2,…,m為一年中侵蝕性降雨的天數；Pk為日降雨大于12 mm的侵蝕性降雨量；Pd12,Py12分別為一年中平均侵蝕降雨量和整個研究時段的平均侵蝕降雨量。

(2) 土壤可蝕性因子(K)。土壤在雨滴打擊、徑流沖刷等外力作用下被分撒和搬運的難易程度被定義為土壤可蝕性。EPIC(erosion-productivity impact calculator)模型[19]應用土壤質地及有機質兩個較為易獲得的因子變量來計算土壤可蝕性，其表達式為：

(4)

(5)

式中：K為土壤侵蝕因子；SAN為沙粒含量；SIL為黏粒含量；CLA為粉粒含量；C為土壤有機碳含量。

(3)LS地形因子。LS地形因子的計算利用Desmet和Govers[20]在考慮水流匯集和溢流后，采用匯水面積來估算某段坡的LS值：

(6)

式中：(LS)i,j為柵格(i,j)處的地形因子;Si,j為坐標(i,j)處的坡度因子;ASi,j-out為坐標(i,j)出口處單位匯水面積;ASi,j-in為坐標(i,j)入口單位匯水面積;m為坡長指數。

(4) 植被因子(C)和水土保持措施因子(P)。植被因子和土地利用有關，我們根據土地利用進行賦值，水土保持工程措施因子主要是一些工程措施和耕作措施，這兩個因子，借鑒相關研究結果獲取[21-23]，結果如圖1所示。

圖1 研究區(qū)主要的土壤侵蝕因子

1.2.3 泥沙連通性指數 本文采用Borselli等[9]提出的徑流泥沙連通性指數IC表征泥沙連通性。圖2給出了泥沙連通性的計算過程，以參考點R的泥沙連通性為例，該方法考慮了參考點R的匯流區(qū)上游部分的基本情況和下游泥沙輸移路徑的基本特征。其上坡模塊Dup考慮了侵蝕源區(qū)坡度、面積和土地利用計算潛在徑流和泥沙向下游流動的可能性，下坡模塊Ddn表示泥沙沿著徑流路線到達指定溝道或者流域出口的可能性，屬于“流動力”式方法?？傮w的計算公式為：

圖2 泥沙連通性的計算過程[14]

(7)

1.2.4 泥沙輸移比計算及驗證 基于泥沙連通性，柵格單元尺度的泥沙輸移比可以按以下公式[14]計算。

(8)

式中：SDRi指柵格i的泥沙輸移比; SDRmax指最大的泥沙輸移比,本研究中為1, IC0及k為待定參數。

利用徑流小區(qū)的實測泥沙流失數據確定IC0和k的值，利用其中的10個徑流小區(qū)確定IC0和k的值，利用剩下的6個徑流小區(qū)驗證計算結果。IC0和k的值確定過程如下(表2)：首先根據徑流小區(qū)的信息計算每個徑流小區(qū)的土壤侵蝕量，然后根據土壤侵蝕量與實測泥沙流失量計算每個徑流小區(qū)的泥沙輸移比，再根據徑流小區(qū)的信息計算泥沙連通性，最后通過泥沙輸移比的計算公式應用最小的誤差試錯法得到推薦的IC0及k組合系數，分別是0，1。

表2 徑流小區(qū)的土壤侵蝕、泥沙連通性、泥沙輸移比及流失量計算

本研究選取Nash-Suttcliffe模型效率系數(NS)作為評價模擬結果的指標。Nash-Suttcliffe模型效率系數(NS)，體現實測與模擬量過程的擬合程度的好壞，其表達式為：

(9)

注：NS為模型效率系數。

2 結果與分析

2.1 土地利用變化

根據2010，2019年遙感影像解譯結果獲得了相應年份的土地利用(表3)，主要表現為耕地面積的減少及其他主要土地利用類型的增加。減少最大的是旱耕地，由2010年的42.73%減少到2019年的35.86%，水田減少幅度稍小從10.03%減少到8.37%。林地增加幅度最大，由2010年的27.06%增加到2019年的31.1%，建筑用地由2010年的280 hm2增加到2019年的980 hm2。

表3 研究區(qū)2010-2019年土地利用變化

2.2 土壤侵蝕變化

利用RUSLE通用土壤流失方程模擬了2010及2019年的土壤侵蝕狀況。根據《土壤侵蝕強度分類分級標準(SL190-2007)》，將研究區(qū)的土壤侵蝕強度劃分為微度、輕度、中度、強度以及極強度5個級別(圖4)。2010和2019年間侵蝕格局幾乎沒有發(fā)生較大的改變(表4)，研究區(qū)域的2019年平均土壤侵蝕模數是2.78 t/(hm2·a)，2010年的侵蝕模數比2019年大，是3.12 t/(hm2·a)。整個研究區(qū)大部分區(qū)域2019年土壤侵蝕級別屬于微度，具體來說，52.64%的區(qū)域土壤侵蝕小于2 t/(hm2·a)。級別屬于輕度侵蝕及以下的區(qū)域占了整個研究區(qū)域的95.62%。

表4 不同年份土壤侵蝕分級統(tǒng)計

圖4 研究區(qū)2010和2019年土壤侵蝕分級

2010和2019年間的土壤侵蝕格局變化不大(圖4)，并統(tǒng)計了2019年不同土地利用侵蝕結果(表5)。從統(tǒng)計結果可知，旱耕地是土壤侵蝕的主要策源地，其土壤侵蝕模數是6.39 t/(hm2·a)，是水田約5倍，整個旱耕地面積占研究區(qū)域35.86%，卻貢獻了侵蝕量的82.59%。林地和草地的土壤侵蝕模數均很低，分別是0.58，0.87，因其很低的土壤侵蝕率，這兩種土地利用面積占了研究區(qū)的52.6%，其侵蝕量只貢獻了13.24%。本研究支持了相關研究的結果[24-25]，不管是在三峽庫區(qū)或是相鄰的川中丘陵區(qū)旱耕地是主要的侵蝕源地，雖然其平均侵蝕速率屬于輕度，但由于土壤侵蝕有著巨大的時空異質性，需要重點關注。

表5 2019年研究區(qū)不同土地利用侵蝕狀況

2.3 泥沙連通性及泥沙輸移比的變化

以2010—2019年不同的土地利用為基礎計算了這兩年的連通性指數，并基于這些指數獲得柵格尺度的泥沙輸移比(圖5—6)。在ArcGIS中統(tǒng)計出不同土地利用類型的連通性指數IC和泥沙輸移比SDR(表6)。整體來說，整個流域的連通性指數呈現減小變化趨勢，具體來說從2010年的-0.46減小到2019年的-0.65。耕地的連通性指數減少幅度最大，如水田由2010年的-0.31減小到2019年的-0.52，旱耕地由2010年的-0.38減小到2019年的-0.66。在2010和2019年，連通性指數最大的是水域，連通指數最小的是林地，主要原因是林地植被覆蓋高，且根系的固土作用不僅是其本身侵蝕指數小，而且具有攔截泥沙的功能，因此使其連通指數小。

表6 研究區(qū)2010-2019年不同土地利用的IC及SDR

圖5 研究區(qū)2010-2019年泥沙連通性指數(IC)

由于連通性指數減小使泥沙輸移比也隨之減小，其中旱耕地減小最大，由2010年的0.29減小到2019年的0.22，建筑用地的泥沙輸移比減小幅度最小。根據2019年各土地利用的泥沙輸移比比較發(fā)現，各主要的土地利用的泥沙輸移比在0.2～0.27之間，差異不是非常大，但是由于不同土地利用土壤侵蝕模數差異較大導致不同土地利用的土壤流失量差異較大。

圖6 研究區(qū)2010-2019年泥沙輸移比(SDR)

根據研究區(qū)的土壤侵蝕模數和泥沙輸移比疊加，統(tǒng)計得到不同土地利用的泥沙流失量(表7)。旱耕地由于土壤侵蝕模數高，因此依然土壤流失模數最高的一類土地利用。旱耕地從2010—2019年泥沙流失量變化幅度最大，由1.97 t/(hm2·a)減小到2019年的1.41 t/(hm2·a)。整個研究區(qū)產沙模數由0.83 t/(hm2·a)減小到0.62 t/(hm2·a)，泥沙流失量從2.38×104t減少到1.76×104t。

表7 研究區(qū)2010-2019年不同土地利用土壤流失量

3 討 論

頁巖氣開采區(qū)的土地利用變化主要表現為耕地減少，采礦用地和草地增加，這種變化格局和整個三峽庫區(qū)趨勢一致[26]。由于我們的研究時間是2010—2019年，這段時間的退耕還林工作已完成即沒有新增加的退耕用地，在此期間的耕地減少主要是因為耕地撂荒導致的，因為耕地撂荒在沒有人為干擾的情況下很快會發(fā)展為草地，因此本研究可以發(fā)現草地的增加幅度最大。相關研究已經證實了重慶是耕地撂荒率較高的地區(qū)，史鐵丑[27]在重慶通過調研比對發(fā)現，耕地撂荒率為18.0%，撂荒耕地以旱地為主，比例為82.4%。

本研究通過RUSLE對研究區(qū)域進行土壤侵蝕的模擬，本研究區(qū)域高差相對較大，但是RUSLE模型已成功對相似地形條件的三峽庫區(qū)成功應用[28]。通過模擬發(fā)現不同的土地利用的土壤侵蝕差異巨大，主要的侵蝕源是旱耕地，達6.39 t/hm2，這與川中丘陵區(qū)不同土地利用侵蝕模數排序相似，然而與川中丘陵區(qū)的耕地相比，該區(qū)域的侵蝕模數偏小[25]。本研究沒有考慮頁巖氣開發(fā)中修建道路及鉆井對土壤擾動造成的土壤流失，因此可能會稍微低估區(qū)域的土壤侵蝕量。

從小區(qū)、坡面到流域尺度，泥沙源區(qū)的連通及響應單元間的空間連通主導著侵蝕的發(fā)生和發(fā)展，進而影響了區(qū)域的產沙[29-30]。本研究發(fā)現在研究區(qū)域近10 a的土地利用變化使泥沙連通性減弱，正如前文所述由于耕地撂荒和頁巖氣的開采使建設用地增加，景觀結構上表現為土地斑塊增加更加破碎化。相關研究發(fā)現，建筑設施等起到阻礙作用的景觀元素削弱泥沙輸移路徑連通，形成較低的連通度[31]。

泥沙連通指數(IC)是分布式的即可以分布整個研究區(qū)，該指數將侵蝕產沙與泥沙輸出結合起來，具有較強的過程基礎，實質上是一個結合侵蝕與泥沙輸移的簡化模型[30]。因此，Vigiak[11]通過改進提出了基于泥沙連通指數的泥沙輸移比。在本研究區(qū)域近10 a的土地利用使區(qū)域的連通指數減少，同樣地泥沙輸移比也相應減少。整個研究區(qū)的泥沙輸移比從0.27減少到0.22。高旭彪[32]根據同位素示蹤法測定了流域土壤侵蝕模數和調查的平均淤積模數總結了紫色土丘陵區(qū)的泥沙輸移比公式，根據此公式本研究區(qū)的泥沙輸移比在0.21～0.25之間，說明基于連通指數建立的分布式泥沙輸移模型具有較好的適用性。

4 結 論

(1) 2010—2019年研究區(qū)域隨著頁巖氣的開發(fā)及耕地撂荒，加速了土地利用結構的變化，建設用地、林草地增加，而相應的耕地減少。

(2) 在2010—2019年，因土地利用格局變化導致的土壤侵蝕變化微弱，整個研究區(qū)內的土壤侵蝕量皆小，從2010年的3.12 t/(hm2·a)變?yōu)?019年為2.78 t/(hm2·a)，整個區(qū)域的侵蝕的策源地是旱耕地，其侵蝕量占整個研究區(qū)的82%。

(3) 在2010—2019年，整個研究區(qū)的連通性指數IC呈現減小變化趨勢，從2010年的-0.46減小到2019年的-0.65，整個區(qū)域泥沙輸移比減小，其中旱耕地減小幅度最大。侵蝕量變化雖不大，但因泥沙輸移比的減少，使整個區(qū)域產沙模數由0.83 t/(hm2·a)減少到2019年的0.62 t/(hm2·a)。