土地利用变化对泥沙连通性及泥沙输移的影响——以重庆市涪陵页岩气开采区为例

冯 斌, 张学伍, 徐 敏, 胡海峰

〔中煤科工重庆设计研究院(集团)有限公司, 重庆 400016〕

土壤侵蚀与泥沙输移是气候、土壤、植被、水文、地形等多种因素的综合结果[1]，代表着陆地生态系统中重要的物质循环与输移过程[2]。流域内土地利用景观类型、数量通过改变土壤侵蚀的发生机制，改变水文结构和侵蚀系统，引起土壤输移、拦截能力的变化，进而影响产沙量的增加或减少[3-4]。随着经济活动的发展，由于人类活动导致的土地利用变化是影响土壤侵蚀的直接驱动力，也是导致水土流失时空分异的重要原因。除了土地利用变化的影响外，流域泥沙输移还受泥沙输移路径影响，泥沙输移和连通性关系密切，水文研究领域借用连通性的概念开展了坡面—沟道水文连通、河—湖和流域水文连通的研究[5-7]。与泥沙输移直接相关的泥沙连通性，是指流域内泥沙通过分离和输移过程中从源到汇的传输程度[8]。泥沙连通性是可以很好地表示流域内泥沙输移潜力的表征量。Borselli等人[9]开发了分布式的泥沙连通性指数(IC)，量化泥沙连通的潜在大小。Cavalli等人[10]对IC指数的坡度、泥沙贡献面积及权重因子的方法进行了改进。以连通性为基础，可以开发出不同空间位置的泥沙输移比[11]，被应用于评价泥沙输移量。流域泥沙连通性主要由流域异质性的空间格局所决定[12]，对景观空间格局的改变均会影响泥沙连通性，进而影响泥沙输移过程，如修建梯田、淤地坝、道路及建设用地增加等[13-14]。总体来说，在考虑流域泥沙输移通量时，现有研究主要关注的是土地利用变化对土壤侵蚀的影响，而忽视了因土地利用变化导致的泥沙连通性的影响。在国外，随着过去20 a间页岩气开发的快速发展，页岩气开发平台、道路和管线占用土地和干扰环境状况逐渐受到关注[15]，页岩区开采需要修建大量的通行道路，进一步增加建设用地，可能会导致区域土地利用变化明显。随着中国页岩气将大规模开发利用，相关学者[16]建议页岩气开发设计应考虑占地、景观破碎化的影响。本研究在考虑页岩气开采区建设对于土地利用的改变的影响基础上，探寻这种景观变化对对泥沙流失影响。本文以重庆市涪陵区中国首批页岩气勘查开发示范区为研究区域，分析开发前后景观格局的变化(2010—2019年)，并利用RUSLE通用土壤流失方程计算土地利用变化导致的土壤侵蚀变化，分析土地利用对泥沙连通性和泥沙输移的影响，研究结果可为类似区域控制土壤流失为目标的流域土地利用布局和优化提供重要的参考价值。

1 研究区概况与研究方法

1.1 研究区概况

研究区位于中国重庆市涪陵区焦石镇，地理坐标为：东经117°26′07′—117°38′13″，北纬29°33′50″—29°47′32″，处于四川盆地及山地过渡地带，面积约为284 km2，地形以低山丘陵为主，研究区范围属于页岩气开采区的范围，来源于重庆市规划与自然资源局。研究区受地质构造和长江、乌江水系的切割，相对高差较大，为141～1 422 m。该区为亚热带湿润季风气候区，雨量充沛，区域年平均降雨1 100 mm，但分布不均，年内降水主要集中于5—9月，占全年降水量的68%左右，年均蒸发量1 076 mm，相对湿度81%；多年平均气温18.0 ℃，≥10 ℃积温5 719 ℃，无霜期多年平均317 d，日照时数1 087 h；区内主导风向为东北向，年平均风速0.7 m/s。研究区土壤以紫色土和水稻土为主；地带性植被类型属亚热带常绿阔叶林，由于人为原因区域内的天然林已大量消失，目前植被主要为人工林、天然次生林及农作物。

1.2 研究方法

1.2.1 数据来源 研究区的页岩气开发始于2012年，到2018年研究区的页岩气开发的附属设施相关建设已基本完成，因此选择2010—2019年为研究时间段。土壤侵蚀模拟数据主要包括日降雨、地形、土地利用、土壤理化参数等。日降雨数据来源于涪陵区气象局，2010—2019年期间日降雨数据，主要用于计算降雨侵蚀力。地形数据是空间栅格大小12.5 m的DEM数据，该数据由ALOS(advanced land observing satellite)卫星相控阵型L波段合成孔径雷达(PALSAR)采集，数据下载网址https:∥urs.earthdata.nasa.gov/users/new，主要用于计算侵蚀模拟中的坡长和坡度因子。土地利用来源于Google earth下载的遥感影像(2010和2019年)，根据中国的《土地利用现状分类》标准和区域的土地利用特点，将土地利用划分为水田、旱耕地、林地、草地、建筑用地(页岩气开采地、道路)。土壤可蚀性的计算依靠有机质和砂、粉、黏3种粒径颗粒的百分比含量，相关数据来源于国家地球系统科学数据中心。本研究利用径流小区的泥沙流失量数据来验证本文的计算结果，径流小区的数据来源于2019年重庆市水土保持公报，从重庆市水利局网站公开下载，本研究共采用了公报中的16个径流小区数据，这些小区泥沙流失量数据观测年份是2018年部分降雨事件结果(表1)。

表1 研究区径流小区数据

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 研究区主要的土壤侵蚀因子

1.2.3 泥沙连通性指数 本文采用Borselli等[9]提出的径流泥沙连通性指数IC表征泥沙连通性。图2给出了泥沙连通性的计算过程，以参考点R的泥沙连通性为例，该方法考虑了参考点R的汇流区上游部分的基本情况和下游泥沙输移路径的基本特征。其上坡模块Dup考虑了侵蚀源区坡度、面积和土地利用计算潜在径流和泥沙向下游流动的可能性，下坡模块Ddn表示泥沙沿着径流路线到达指定沟道或者流域出口的可能性，属于“流动力”式方法。总体的计算公式为：

图2 泥沙连通性的计算过程[14]

(7)

1.2.4 泥沙输移比计算及验证 基于泥沙连通性，栅格单元尺度的泥沙输移比可以按以下公式[14]计算。

(8)

式中：SDRi指栅格i的泥沙输移比; SDRmax指最大的泥沙输移比,本研究中为1, IC0及k为待定参数。

利用径流小区的实测泥沙流失数据确定IC0和k的值，利用其中的10个径流小区确定IC0和k的值，利用剩下的6个径流小区验证计算结果。IC0和k的值确定过程如下(表2)：首先根据径流小区的信息计算每个径流小区的土壤侵蚀量，然后根据土壤侵蚀量与实测泥沙流失量计算每个径流小区的泥沙输移比，再根据径流小区的信息计算泥沙连通性，最后通过泥沙输移比的计算公式应用最小的误差试错法得到推荐的IC0及k组合系数，分别是0，1。

表2 径流小区的土壤侵蚀、泥沙连通性、泥沙输移比及流失量计算

本研究选取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 研究区2010-2019年土地利用变化

2.2 土壤侵蚀变化

利用RUSLE通用土壤流失方程模拟了2010及2019年的土壤侵蚀状况。根据《土壤侵蚀强度分类分级标准(SL190-2007)》，将研究区的土壤侵蚀强度划分为微度、轻度、中度、强度以及极强度5个级别(图4)。2010和2019年间侵蚀格局几乎没有发生较大的改变(表4)，研究区域的2019年平均土壤侵蚀模数是2.78 t/(hm2·a)，2010年的侵蚀模数比2019年大，是3.12 t/(hm2·a)。整个研究区大部分区域2019年土壤侵蚀级别属于微度，具体来说，52.64%的区域土壤侵蚀小于2 t/(hm2·a)。级别属于轻度侵蚀及以下的区域占了整个研究区域的95.62%。

表4 不同年份土壤侵蚀分级统计

图4 研究区2010和2019年土壤侵蚀分级

2010和2019年间的土壤侵蚀格局变化不大(图4)，并统计了2019年不同土地利用侵蚀结果(表5)。从统计结果可知，旱耕地是土壤侵蚀的主要策源地，其土壤侵蚀模数是6.39 t/(hm2·a)，是水田约5倍，整个旱耕地面积占研究区域35.86%，却贡献了侵蚀量的82.59%。林地和草地的土壤侵蚀模数均很低，分别是0.58，0.87，因其很低的土壤侵蚀率，这两种土地利用面积占了研究区的52.6%，其侵蚀量只贡献了13.24%。本研究支持了相关研究的结果[24-25]，不管是在三峡库区或是相邻的川中丘陵区旱耕地是主要的侵蚀源地，虽然其平均侵蚀速率属于轻度，但由于土壤侵蚀有着巨大的时空异质性，需要重点关注。

表5 2019年研究区不同土地利用侵蚀状况

2.3 泥沙连通性及泥沙输移比的变化

以2010—2019年不同的土地利用为基础计算了这两年的连通性指数，并基于这些指数获得栅格尺度的泥沙输移比(图5—6)。在ArcGIS中统计出不同土地利用类型的连通性指数IC和泥沙输移比SDR(表6)。整体来说，整个流域的连通性指数呈现减小变化趋势，具体来说从2010年的-0.46减小到2019年的-0.65。耕地的连通性指数减少幅度最大，如水田由2010年的-0.31减小到2019年的-0.52，旱耕地由2010年的-0.38减小到2019年的-0.66。在2010和2019年，连通性指数最大的是水域，连通指数最小的是林地，主要原因是林地植被覆盖高，且根系的固土作用不仅是其本身侵蚀指数小，而且具有拦截泥沙的功能，因此使其连通指数小。

表6 研究区2010-2019年不同土地利用的IC及SDR

图5 研究区2010-2019年泥沙连通性指数(IC)

由于连通性指数减小使泥沙输移比也随之减小，其中旱耕地减小最大，由2010年的0.29减小到2019年的0.22，建筑用地的泥沙输移比减小幅度最小。根据2019年各土地利用的泥沙输移比比较发现，各主要的土地利用的泥沙输移比在0.2～0.27之间，差异不是非常大，但是由于不同土地利用土壤侵蚀模数差异较大导致不同土地利用的土壤流失量差异较大。

图6 研究区2010-2019年泥沙输移比(SDR)

根据研究区的土壤侵蚀模数和泥沙输移比叠加，统计得到不同土地利用的泥沙流失量(表7)。旱耕地由于土壤侵蚀模数高，因此依然土壤流失模数最高的一类土地利用。旱耕地从2010—2019年泥沙流失量变化幅度最大，由1.97 t/(hm2·a)减小到2019年的1.41 t/(hm2·a)。整个研究区产沙模数由0.83 t/(hm2·a)减小到0.62 t/(hm2·a)，泥沙流失量从2.38×104t减少到1.76×104t。

表7 研究区2010-2019年不同土地利用土壤流失量

3 讨 论

页岩气开采区的土地利用变化主要表现为耕地减少，采矿用地和草地增加，这种变化格局和整个三峡库区趋势一致[26]。由于我们的研究时间是2010—2019年，这段时间的退耕还林工作已完成即没有新增加的退耕用地，在此期间的耕地减少主要是因为耕地撂荒导致的，因为耕地撂荒在没有人为干扰的情况下很快会发展为草地，因此本研究可以发现草地的增加幅度最大。相关研究已经证实了重庆是耕地撂荒率较高的地区，史铁丑[27]在重庆通过调研比对发现，耕地撂荒率为18.0%，撂荒耕地以旱地为主，比例为82.4%。

本研究通过RUSLE对研究区域进行土壤侵蚀的模拟，本研究区域高差相对较大，但是RUSLE模型已成功对相似地形条件的三峡库区成功应用[28]。通过模拟发现不同的土地利用的土壤侵蚀差异巨大，主要的侵蚀源是旱耕地，达6.39 t/hm2，这与川中丘陵区不同土地利用侵蚀模数排序相似，然而与川中丘陵区的耕地相比，该区域的侵蚀模数偏小[25]。本研究没有考虑页岩气开发中修建道路及钻井对土壤扰动造成的土壤流失，因此可能会稍微低估区域的土壤侵蚀量。

从小区、坡面到流域尺度，泥沙源区的连通及响应单元间的空间连通主导着侵蚀的发生和发展，进而影响了区域的产沙[29-30]。本研究发现在研究区域近10 a的土地利用变化使泥沙连通性减弱，正如前文所述由于耕地撂荒和页岩气的开采使建设用地增加，景观结构上表现为土地斑块增加更加破碎化。相关研究发现，建筑设施等起到阻碍作用的景观元素削弱泥沙输移路径连通，形成较低的连通度[31]。

泥沙连通指数(IC)是分布式的即可以分布整个研究区，该指数将侵蚀产沙与泥沙输出结合起来，具有较强的过程基础，实质上是一个结合侵蚀与泥沙输移的简化模型[30]。因此，Vigiak[11]通过改进提出了基于泥沙连通指数的泥沙输移比。在本研究区域近10 a的土地利用使区域的连通指数减少，同样地泥沙输移比也相应减少。整个研究区的泥沙输移比从0.27减少到0.22。高旭彪[32]根据同位素示踪法测定了流域土壤侵蚀模数和调查的平均淤积模数总结了紫色土丘陵区的泥沙输移比公式，根据此公式本研究区的泥沙输移比在0.21～0.25之间，说明基于连通指数建立的分布式泥沙输移模型具有较好的适用性。

4 结 论

(1) 2010—2019年研究区域随着页岩气的开发及耕地撂荒，加速了土地利用结构的变化，建设用地、林草地增加，而相应的耕地减少。

(2) 在2010—2019年，因土地利用格局变化导致的土壤侵蚀变化微弱，整个研究区内的土壤侵蚀量皆小，从2010年的3.12 t/(hm2·a)变为2019年为2.78 t/(hm2·a)，整个区域的侵蚀的策源地是旱耕地，其侵蚀量占整个研究区的82%。

(3) 在2010—2019年，整个研究区的连通性指数IC呈现减小变化趋势，从2010年的-0.46减小到2019年的-0.65，整个区域泥沙输移比减小，其中旱耕地减小幅度最大。侵蚀量变化虽不大，但因泥沙输移比的减少，使整个区域产沙模数由0.83 t/(hm2·a)减少到2019年的0.62 t/(hm2·a)。