植被覆盖和降雨因子变化及对东北黑土区土壤侵蚀的影响

王 硕, 方海燕, 和继军

(1.首都师范大学 城市环境过程和数字模拟国家重点试验室培育基地, 北京100048; 2.中国科学院地理科学与资源研究所 陆地水循环及地表过程重点试验室, 北京 100101; 3.中国科学院大学 资源与环境学院， 北京 100049)

土壤侵蚀是世界性的环境问题之一[1]，其不仅破坏土地资源，降低土壤肥力，还会诱发洪涝灾害，严重影响生态环境和可持续发展[2]。因此，探索土壤侵蚀变化的原因，成为控制水土流失的关键。20世纪Wischmeier[3-4]建立了通用土壤流失方程(USLE)，该方程将诸多影响土壤侵蚀速率的物理和管理因素归纳为6个主要因子，分别为降雨侵蚀力因子R，土壤可蚀性因子K，坡度坡长因子LS，植被覆盖与管理因子C和水土保持措施因子P。近些年来，全球气候变化开始对土壤侵蚀产生影响，同时大规模的人为活动在加重土壤侵蚀的产生[5]。国内外都开展了这方面的工作，例如Teng[6]和Chakrabortty[7]曾在青藏高原和印度东部研究气候变化影响下的土壤侵蚀，吴昌广[8]和余新晓[9]也曾结合降雨侵蚀力与植被覆盖因子在三峡与黄土区进行土壤侵蚀研究，最终表明植被覆盖和降雨侵蚀力因子已成为影响土壤侵蚀变化的两个关键因素。

东北黑土区是中国重要的商品粮基地[10]，素有“北大荒”之称。然而由于黑土区独特的地理环境[11]，不合理的开垦耕作[12]以及近些年气候变暖的影响，使得黑土区内植被覆盖率和降雨发生了巨大变化[13]，最终导致其水土流失现象日益突出。许多研究者[14-16]曾采用RUSLE模型，在东北黑土区的诸多区域开展过工作。然而以往开展的工作中，单因子变化对侵蚀造成影响的研究相对较少，特别是基于植被覆盖和降雨因子变化对整个黑土区侵蚀影响的研究更是未见报道。近20 a，该区域植被覆盖率和降雨量变化较大，因此有必要针对这些变化的影响程度进行探讨，有利于揭示东北黑土区土壤侵蚀的空间变化特征，对于深化该区的土壤侵蚀研究更具深远意义。因此，本文利用气象、遥感等资料，采用RUSLE模型估算土壤侵蚀，并将东北黑土区划分为6个区域，从时空角度结合降雨侵蚀力因子、植被覆盖与管理因子变化来揭示土壤侵蚀空间差异及其对降雨侵蚀力和植被覆盖变化的敏感性，进而为水土保持规划提供数据支持。

1 研究区与研究方法

1.1 研究区概况

东北黑土区(38°42′—53°36′，115°24′—135°12′)包括黑龙江、吉林、辽宁及内蒙古东北部地区，总面积约为1.24×106km2，海拔在0～2 667 m之间。大兴安岭、小兴安岭和长白山分别坐落于研究区的西部、东北和东南，辽河平原、松嫩平原、三江平原三大平原位于三大山脉之间并且由南至北排开，辽河、嫩江、松花江流淌其中。研究区属于温带季风气候，四季分明。近40 a年均降雨量由东南部的900 mm衰减至西北部的300 mm左右，近70%降水来自夏季。研究区主要有黑土、黑钙土、潮土、白浆土等土壤类型。山区以森林为主，平原地区以耕地为主。为了方便研究区进行描述，本文依据松嫩平原区域、嫩江、松花江及辽河等地势与水系为边界，将研究区分为6个区域，分别为辽东半岛及长白山南部(Ⅰ区)、三江平原及长白山北部(Ⅱ区)、小兴安岭(Ⅲ区)、蒙古北部和大兴安岭北部(Ⅳ区)、辽河平原及大兴安岭中南部(Ⅴ区)及松嫩平原(Ⅵ区)。

1.2 数据来源

本文2000，2005，2010，2015和2018年0.5°×0.5°的日降雨数据来自国家气象科学中心。土壤数据来自联合国粮食与农业组织The Harmonized World Soil Database (HWSD)。数字高程影像(DEM； 90 m分辨率)数据来源于中国科学院地理空间数据云平台。2000，2005，2010，2015和2018年归一化植被指数(NDVI)与土地利用数据(100 m分辨率)均来自中国科学院资源环境科学数据中心。

1.3 研究方法

由于本研究区较大，RUSLE模型相比于CSLE模型而言其参数更易获取，且该模型已经在其他大尺度地区成功开展过工作，C，P因子能有所借鉴。因此本文采用RUSLE模型计算土壤侵蚀强度，表达式为：

A=R·K·L·S·C·P

(1)

式中：A为年土壤侵蚀模数(t/hm2·a)；R为降雨侵蚀力因子〔MJ·mm/(hm2·h)〕；K为土壤可蚀性因子〔(t·h)/(MJ·mm)〕；L为坡长因子(无单位)；S为坡度因子(无单位)；C为植被覆盖与管理因子(无单位)；P为水土保持措施因子。

1.3.1 降雨侵蚀力R因子R值反映了降雨对土壤的潜在侵蚀能力。章文波[17]指出我国东北地区利用逐日降雨数据计算R值精度较高。因而，本研究采用日降雨数据计算R值：

(2)

式中：Mi为第i个半月时段的侵蚀力值〔(MJ·mm)/(hm2·h·a)〕;K为该半月时段内的天数；Dj为半月时段内第j天的侵蚀性日雨量，要求日雨量大于12 mm，否则以0计算;α，β为模型的待定系数。

(3)

α=21 568×β-7.189

(4)

式中：Pd12为≥12 mm的日平均雨量(mm)；Py12为≥12 mm的年平均雨量(mm)。通过以上方法计算R并对其进行插值[18]，最终得到R值分布。

1.3.2 土壤可蚀性K因子 张科利[19]表明修正后的EPIC模型在我国东北地区适用性较好。因此本研究通过Sharpley[20]和Williams[21]提出的EPIC模型方法，再结合张科利修正公式进行计算：

(5)

K=0.515 75×Kepic-0.013 83

(6)

式中：Kepic为修正前的K因子；CLA，SIL，SAN分别为土壤黏粒,粉粒和砂粒含量(%);C为有机碳含量(%); SN1=1-SAN/100。

1.3.3 坡度坡长LS因子LS反映地貌特征对土壤侵蚀的影响。考虑到我国东北地势起伏，而RUSLE是美国农业部以美国缓坡农用地为对象建立的。因此，本研究采用Liu[22-23]提出基于我国陡坡的LS因子修正方法：

(7)

(8)

式中：λ为坡长(m);θ为坡度(°)。

1.3.4 植被覆盖与管理C因子 由于研究区较大，植被覆盖度变化于0%～76%，为避免出现斑块化结果和准确获取C值，经参考黑土区相关文献[24-25]后采用蔡崇法[26]建立的C值模型计算C因子值。计算方法为：

(9)

(10)

式中：c为植被覆盖度; NDVI为归一化植被指数; NDVImax为研究区全植被覆盖的NDVI值; NDVImin为全裸地或无植被覆盖的NDVI值。

1.3.5 水土保持措施P因子 根据《黑土区水土流失综合防治技术标准》，规定<3°坡耕地实施等高改垄；3°～5°坡耕地采取地埂植物带；>5°修筑梯田；研究区内的坡耕地进行了大规模治理[28]。根据相关文献[29]赋值方法，最终将水域区P值取0，梯田P值为0.03，等高改垄与地埂植物带的旱田P值取0.352，顺耕耕作的农田及其他自然植被区P值取1。

1.4 敏感系数

本文采用敏感系数法来衡量侵蚀模数对因子变化的敏感程度。计算公式为：

(11)

式中：Q为敏感度,表示侵蚀模数在影响因子变动下的变化快慢,数值越大敏感性越强; ΔA为侵蚀模数年均变化量; ΔX为影响因子的年均变化量。

(12)

(13)

式中：A1为研究期第一年的侵蚀模数;A2为研究期最后一年的侵蚀模数;X1为研究期第一年的平均因子数值;X2为研究期最后一年的平均因子数值;T为研究期时间跨度(a)。

为了消除数量级的影响，将敏感度进行标准化处理：

(14)

式中：M为标准化敏感度，即敏感系数；Qmax为区域内敏感度最大值；Qmin区域敏感度最小值。敏感系数分级[30]详见表1。

表1 土壤侵蚀敏感系数的分级

1.5 数据处理

本文借助ArcGIS 10.2软件，利用2000，2005，2010，2015和2018年5期C因子与其他各因子的多年平均值进行运算，得到C因子变化下土壤侵蚀模数；同理，利用2000，2005，2010，2015和2018年5期R因子数据得到R因子变化下的土壤侵蚀模数；再结合5期对应年份的各个因子数据，得到研究区内各个年份的土壤侵蚀模数。最后根据水利部《土壤侵蚀分类分级标准(SL190-2007)》，将土壤侵蚀模数进行重分类，得到C因子变化下土壤侵蚀模数时空分布、降雨侵蚀力因子变化下的侵蚀模数分布和各个年份的土壤侵蚀模数分布。

2 结果与分析

2.1 模拟结果的验证

2.1.1C因子结果验证 本文采用的C值模型结果与张宪奎[27]和张雪花[16]的C因子赋值标准基本一致，因此该C值模型在东北黑土区可适用(表2)。

表2 C因子结果验证对比

2.1.2 侵蚀模数结果验证 在以往的研究中，东北黑土区采用137Cs示踪、径流小区监测与模型估算等方法开展了侵蚀研究的相关工作。例如，方华军[31]、李国强[32]与张克新[33]同样采用137Cs示踪法分别在松花江镇、拜泉县与辽东湾的得到了220，245与1 739～3 892 t/(km2·a)的结果，与本研究所得的187.62，237.46与2 568.57～4 323.36 t/(km2·a)的估算值相差不大；刘宝元[34]还在鹤山农场利用径流小区监测得到坡耕地的侵蚀模数为845～1 157 t/(km2·a)，对比所用数据为鹤山县的土壤侵蚀模数，因此小于径流小区结果。还有盛美玲[35]利用WATEM/SEDAM模型在得到的2005年土壤侵蚀模数为351.2 t/(km2·a)，顾治家[36]利用CSLE模型得到的结果为677 t/(km2·a)，Fang[29]利用RUSLE模型在东北黑土区尺度下计算2010年的侵蚀模数为943.7 t/(km2·a)，本研究结果均与以上模型估算值相近(表3)。本研究中土壤侵蚀模数偏小的原因可能来自两个方面，其一是当前所能获取的NDVI与土地利用数据分辨率过低平滑了数值[38]，降低了土壤侵蚀模数；其二则是研究时间和区域、以及赋值精细度不同所造成。

表3 土壤侵蚀模数结果验证对比

2.2 RUSLE各因子分布及变化特征

2.2.1 研究区R因子变化分布及变化特征 2000—2018年R因子呈现先增加后降低再增加的趋势。空间上，R因子由西南向东北部减少，区域上呈现：Ⅰ>Ⅱ>Ⅵ>Ⅲ>Ⅴ>Ⅳ的分布特点(图1)。从局部上看，Ⅰ区是R值最大区并且以丹东市为中心向四周降低。辽东半岛东侧与西侧受到长白山脉阻隔造成降雨量差距明显，其中2018年东西侧R值相差近2 000 (MJ·mm)/(hm2·h·a)。

图1 研究区各年份降雨侵蚀力因子分布

在Ⅱ，Ⅲ和Ⅵ区R因子由765.42 (MJ·mm)/(hm2·h·a)增至1 417.79 (MJ·mm)/(hm2·h·a)。研究区内最小R值出现在Ⅳ和Ⅴ区，该区R值稳定升高但变化幅度较小，呈现大兴安岭中部大于南北部的分布特点。

2.2.2 研究区K因子与LS因子变化分布及变化特征 研究区内K因子最高值分布在大兴安岭南部(Ⅴ)和辽河平原北部(Ⅵ)内；最低值则出现在大兴安岭北部(Ⅳ)、辽河平原西南部(Ⅴ)和松嫩平原(Ⅵ)。研究区内LS因子的高值主要分布在三大山脉，研究区内的三大平原均由低值构成。

2.2.3 研究区C因子变化分布及变化特征 2000—2018年C值降低了42%，呈现中部与西南部偏高，东部与北部偏低的特征(图2)。C因子的变化主要受夏季平均气温影响[39]。数据表明[40]，2000—2011年夏季平均气温持续升高，而2011年后夏季平均气温开始稳定，这趋势与C值变化曲线相吻合。C因子降低受政策因素的作用也不容忽视。随着三北防护林、退耕还林等工程的大规模实施，东北全区植被覆盖面积出现了稳固的增加[41]。从局部看，Ⅵ区受到政策因素影响C值呈下降趋势，18 a间C值小于0.1的面积占比由27.29%增至69.62%，为植被覆盖率增加最快的区域。Ⅰ，Ⅱ，Ⅲ和Ⅳ区的C值处于持续降低的态势，该范围曾经存在大量的高值区域(C>0.1)，至2010年高值区域面积占比由31.62%下降至12.23%，这也与夏季平均气温的逐年升高息息相关。Ⅴ区在研究期间C值低于0.1的面积占比总体增加了34.12%。该区的C值与海拔高度的气温呈现负相关[42]，在高海拔地区常年气温较低，因此造成高海拔地区C值并未发生明显改变。

图2 研究区各年份植被覆盖与管理因子分布

2.2.4 研究区P因子变化分布及变化特征 2000—2018年平均P值变化率不足1%，且变化均出现于黑龙江的北部(Ⅳ)和大兴安岭南部(Ⅴ)。相比而言，C与R因子变化率都超过了50%。因此，本文仅考虑C因子与R因子对于东北黑土区土壤侵蚀模数变化的影响。

2.3 单因子变化下的土壤侵蚀时空特征

2.3.1 研究区C因子变化下土壤侵蚀特征 在C因子变化下，2000—2018年研究区土壤侵蚀模数持续降低，呈现南部侵蚀严重、北部侵蚀轻微的特征(图3和表4)。局部上看，研究区内中部和北部主要以微轻度侵蚀为主。Ⅱ，Ⅲ和Ⅳ区侵蚀模数与C因子的年均变化量都较小，敏感系数均小于0.3，属不敏感区域。Ⅵ区C因子年均变化量达0.021，高于全区的平均水平，该因子大幅度下降势必导致该区侵蚀模数随之大幅度下降。根据敏感度计算公式，Ⅵ区敏感度为0，属不敏感区域，表明该区在单位植被覆盖度增加的情况下，对侵蚀的缓解程度不明显。以上不敏感区域大部分处于地势平坦且土壤可蚀性较低的地区，所以低值的K与LS因子削弱了植被覆盖对侵蚀变化的影响。

表4 研究区不同分区土壤侵蚀模数变化量与敏感系数

图3 研究区C变化下的土壤侵蚀模数分布

土壤侵蚀最严重的区域为Ⅰ与Ⅴ区，研究期间强度以上侵蚀向中度以下侵蚀转变的面积近63 000 km2。这两个区域的C因子年均变化量仅为Ⅵ区的1/2，但其年均侵蚀变化量却为Ⅵ区的3～4倍，分别达到了-214.33与-271.31 t/(km2·a)，同时Ⅰ与Ⅴ区的敏感系数分别为0.95和1.00，是研究区内的两个强度敏感区域。Ⅰ区和Ⅴ区分别坐落于长白山中南部与大兴安岭中南部，地势陡峭，LS值最大；并且Ⅰ区濒临渤海，降雨充沛，R值最大；而Ⅴ区是K，P值最高的区域。在其他因子值较高的情况下，C因子若发生轻微变化，势必会造成侵蚀模数发生较大的改变。因此应优先提高这两个区域的植被覆盖程度，能够更为高效的减少土壤侵蚀。

2.3.2R因子变化下土壤侵蚀时空特征 2000—2018年，受R因子变化的影响，土壤侵蚀模数呈现先增加后减少再增加的趋势，空间上呈现西南部侵蚀严重，东北部侵蚀轻微的特点(图4和表5)。局部来看，Ⅰ区年均侵蚀增加量虽为全区最大的69.81 t/(km2·a)，但该区的敏感系数仅为0.08，属不敏感区域。Ⅱ，Ⅲ和Ⅵ区侵蚀模数变化小，敏感系数均在0.1以下，是不敏感区域。以上4个区域不敏感的原因在于C因子较小，较高的植被覆盖可有效拦截降雨，减少降雨击溅侵蚀，缓解强降雨对侵蚀造成的影响。在R因子偏低且浮动较小的Ⅳ和Ⅴ区，出现较严重的侵蚀和侵蚀强度变化，敏感系数分别为1.00和0.45，属中度敏感和强度敏感。Ⅳ和Ⅴ区分布着大兴安岭与蒙古高原，地势陡峭并且包含大量的K，LS，C和P因子的高值区域，便造成这两个分区敏感。因此该区域应优先减少坡耕地，实行免耕等耕作方式，降低降雨对侵蚀所造成的影响。

图4 研究区R变化下的土壤侵蚀模数分布

表5 研究区内区域变化量与敏感系数

2.4 不同年份土壤侵蚀时空特征

研究期间东北黑土区年均侵蚀模数为893.53 t/(km2·a)，远超过该区的容许土壤流失量200 t/(km2·a)。Ⅰ和Ⅴ区侵蚀较为严重，其余地区属于微度侵蚀和轻度侵蚀且变化较小。从土壤侵蚀面积来看，从2000—2015年研究区微度侵蚀面积占比增加了3.73%，而中度以上侵蚀的面积占比缩幅为3.84%，说明研究区出现中度以上侵蚀向轻微度侵蚀转移的迹象。从2015—2018年，微度侵蚀面积缩小了55 989.67 km2，轻度侵蚀面积增长了50 462.60 km2，研究期后3 a出现微度向轻度以上侵蚀转移的情况(图5和表6)。

表6 研究区各年份不同土壤侵蚀强度分级面积 t/(km2·a)

图5 研究区各年份的土壤侵蚀模数分布

从空间局部上看，研究期间Ⅰ区侵蚀模数呈现减弱趋势。极强度以上侵蚀面积由9 914.18 km2减少到1 793.61 km2。该区侵蚀模数对C因子呈强度敏感，研究期C因子由0.53下降至0.17，使该区土壤侵蚀强度降低。但其地处长白山南部，地势陡峭，濒临渤海，常年雨量大，也使其侵蚀模数属于全区偏大的水平。

辽河平原(Ⅴ)是极强度以上侵蚀的主要分布区，侵蚀模数呈现先减后增的趋势。该区侵蚀对C因子与R因子皆为强度敏感。在2000—2010年，该区C因子由0.55骤降至0.19，使得侵蚀得到极大改善。在2010—2015年，该区R值由2 634.83 (MJ·mm)/(hm2·h·a)降至813.98 (MJ·mm)/(hm2·h·a)，同样改善了侵蚀情况。2015至2018年间，C因子出现30.44%的回升，导致该区侵蚀模数的少量增加。

大兴安岭(Ⅴ)地区属于对R与C因子敏感地区，呈现先增后减再增的趋势。2000年强度以上侵蚀面积为27 848.86 km2，后10 a侵蚀开始由大兴安岭南部蔓延至大兴安岭西南部，强度以上侵蚀面积达到34 475.68 km2，占该年强度以上侵蚀总面积的73.72%。前10 a中R因子增幅达42.85%，导致该区侵蚀严重化；2010—2015年，该区C与R值同时降低，改善了该区的侵蚀情况。2015—2018年，该区R和C值同时增加，强度以上侵蚀面积再次出现增加，达到了24 753.39 km2；这些年份尽管大兴安岭中部的侵蚀情况有所缓解，但大兴安岭南部与内蒙古高原接壤地带仍是高侵蚀地带，这与该区陡峭的地形有关。

Ⅰ区和Ⅴ区尽管侵蚀情况得到了相当大程度的缓解，但大部分区域仍然处于中度侵蚀的水平。而该地区属于对C因子强度敏感且地势陡峭的地区，因此应当制定相应政策减少坡耕地，增加植被覆盖的程度，能够更高效的减少当地的土壤侵蚀。Ⅳ区虽一直处于轻中度侵蚀的水平，但研究期间侵蚀模数出现了明显的增加，仅此更应该尽早的制定农田免耕，梯田等水土保持措施来缓解降雨的影响。

3 结 论

(1) 整体上，东北黑土区平均土壤侵蚀量有所下降，其原因主要为C因子和R因子的变动，其表现是部分地区降雨量的增减、植被覆盖的变化。

(2) 2000—2018年C因子呈持续降低的趋势；该因子呈现Ⅵ和Ⅴ区偏高，Ⅰ，Ⅱ，Ⅲ和Ⅳ偏小的分布特征。在该因子影响下，侵蚀模数从1 551.07 t/(km2·a)降低至665.39 t/(km2·a)。研究区Ⅰ和Ⅴ的敏感系数分别为0.95和1.00，是强度敏感的地区；Ⅱ，Ⅲ，Ⅳ和Ⅵ区敏感系数最低，属于不敏感地区。

(3) 研究期间R因子呈现先增加后降低再增加的趋势，空间上呈现：区域Ⅰ>Ⅱ>Ⅵ>Ⅲ>Ⅴ>Ⅳ特征。仅在R因子影响下，侵蚀模数呈先增加后降低再增加的趋势，Ⅴ与Ⅳ区敏感系数分别为1.00和0.45，分别属于强度敏感和中度敏感地区；Ⅰ，Ⅱ，Ⅲ和Ⅵ区的敏感系数均较低，属于不敏感地区。

(4) 多年来，东北黑土区区土壤侵蚀主要以微度和轻度侵蚀为主，两类侵蚀面积占比在92.3%～96.2%之间波动。强度侵蚀以上的地区基本分布在Ⅰ和Ⅴ区，侵蚀强度从西南到东北逐渐递减。研究期侵蚀模数从1 175.20 t/(km2·a)变化到822.07 t/(km2·a)，且2015年还出现过628.73 t/(km2·a)的最低值。