2017年8月8日21时19分46秒在四川省北部阿坝州九寨沟县发生7.0级地震，震中位于九寨沟核心景区西部5 km附近。强震往往会导致当地的土壤侵蚀量增加，造成当地水土流失和环境破坏。例如，“5·12”汶川地震发生后，在地震影响下，受灾区的水土保持功能在不同程度上受损：四川理县震后侵蚀量相比2001年增加了12万t/a左右^[1]。茶坪山区水土严重受损，震后侵蚀量相比震前1年增加了72.3%^[2]。龙门山区相比震前1年侵蚀面积增加了86.3 km^2[3]。北川县震后土壤侵蚀面积相比2000年全国土壤侵蚀遥感的调查成果增加275.13 km^2[4]。1999年台湾集集发生的M-W7.6级大地震，在地震初期年均侵蚀深度达到6.04 mm，是震前的2.26倍^[5]。由此可见，地震会导致当地的水土保持功能受损，在侵蚀严重的区域，人民的生产生活和环境会受到严重的威胁；因此，在灾区不同区域，特别是像九寨沟这样以旅游业为主的区域，对土壤侵蚀进行定量评价，划分出受地震影响严重的区域，并分析该区域的空间分布规律，能够有助于现在或将来将有限的资源最有效的投入到灾后水土保持工作中，这对灾后生态恢复具有现实的意义。国内外已有大量关于土壤侵蚀的研究。目前土壤侵蚀的模型主要分2种：1)以观测数据建立的经验模型；2)基于侵蚀机理建立的物理模型。经验模型出现的时间比较早，且模型结构也比物理模型精简，计算更加方便。在经验模型方面，美国于1985年提出了通用土壤流失方程(USLE)，被美国农业部以农业手册的方式颁布执行。USLE模型为土壤侵蚀研究提供了一个简单实用的工具，在全世界范围内得到广泛的应用。早期该模型在降雨条件下，可以预测年均土壤侵蚀量，可以计算出由于细沟侵蚀、面蚀造成的水土流失。随着空间信息技术的发展，运用USLE模型对区域尺度上的土壤侵蚀进行定量评价成为现实。如Batjes^[6]利用0.5经度乘0.5纬度网格全球数据，通过USLE模型实现全球范围内土壤侵蚀的定量评价。我国学者通过修改USLE模型参数因子的计算方法，得到适合于我国土壤状况的土壤侵蚀模型。刘宝元等^[7]在2001基于USLE模型提出了中国土壤流失方程(CSLE)。卜兆宏等^[8]在2003年改进的侵蚀模型与USLE形式一样，只是修正了各个侵蚀因子的算法。

结合前人研究，修正USLE的侵蚀因子，得到适合于研究区的侵蚀模型，定量计算出九寨沟地震前后受灾区的土壤侵蚀强度，对地震前后土壤侵蚀参数做简要的对比分析，划分出受地震影响较大的区域。本研究可初步得出震后水土保持功能易受损区域的空间分布规律，为震后水土流失防治提供指导，对震后灾区尤其是受灾景区的水土保持工作提供科学依据。

1 研究区概况

九寨沟县地处四川盆地和青藏高原2大地貌单元过渡的深切割高山峡谷地带，海拔1 844~4 545 m，地势西北、西南高，东南低，大部分区域的坡度在30°以上。九寨沟县位于中国的西南部，主要受西南季风和东南季风交替影响，多年平均气温为7 ℃左右, 年平均降水量600 mm，雨季集中在每年的5月到10月，约占全年降雨量的3/4。森林覆盖率达54.9%，植被群落随地貌、气候、土壤和海拔高度的不同而有差异。2017年8月8日九寨沟地震震中位于九寨沟县境内，这场地震造成当地山体破碎、土质疏松、植被破坏，严重加重当地水土保持负担。根据震前、震后的遥感影像差异，确定受地震影响的区域作为研究区(图 1)。

2 材料与方法 2.1 土壤侵蚀量预测模型

本文的侵蚀模型采用通用土壤侵蚀模型(USLE)来定量计算震前、震后的土壤流失量。目前，该模型在我国很多地方都被广泛应用，如黄炎和等^[9]和周伏建等^[10]运用该模型分别在福建省和闽东南地区进行了实践。耿豪鹏等^[11]应用USLE模型在黄土高原西部进行了实践。王娇等^[12]在USLE模型的理论指导下，研究了太行山区水土保持的空间分布特征。赵磊等^[13]将USLE模型的各个因子进行相乘得出滇池某流域的土壤侵蚀量。USLE模型方程为

$ A=RKLSPC。$

(1)

式中：A是土壤流失量，t/(hm²·a)；R是降雨可蚀性的因子，MJ·mm/(hm²·h·a)；K是土壤可蚀性的因子，t·hm²·h/(MJ·mm·hm²)；L、S分别为坡长和坡度因子；P是水土保持措施因子；C是作物和覆盖因子。

土壤侵蚀强度预测技术路线如图 2所示。

2.2 数据源

基础数据及图件主要包括2015年和2018年5 m分辨率DEM，四川行政区划矢量图，2015年8月21日和2017年8月9日遥感影像，四川1:100万土壤类型图，2015年1:10万土地利用栅格图，九寨沟县及周边县城(1981—2010)20年的降雨数据。DEM由四川测绘局提供，遥感影像来源于Landsat 4~5和高分2号，土壤类型图来源于中国第2次普查数据库，相关参数参考《四川土种志》。降雨数据来自中国气象数据网。土地利用分布图来自中国科学院地理科学与资源研究所。在ArcGIS中采用样条插值方法将降雨数据生成5 m×5 m的栅格数据。所有成果均采用WGS1984地理坐标系和高斯-克吕格投影，该坐标投影系统被认为是中国应用最广泛的大地坐标投影系统，栅格大小设为5 m×5 m。

3 土壤侵蚀因子提取 3.1 降雨侵蚀因子

降雨侵蚀因子与降雨历时、降雨量、降雨类型、雨强等因素有关，它反映了降雨特征对土壤侵蚀的影响。由于降雨侵蚀力很难测得，一般依据降雨参数来估算。本研究采用我国2011年全国水土流失普查使用的降雨可蚀性计算方法^[14]，其模型如下：

$ R=\frac{1}{N}\sum\limits_{i=1}^{N}{\left( \alpha \sum\limits_{j=1}^{N}{P_{{{\text{d}}_{ikj}}}^{\beta }} \right)}。$

(2)

其中：

$ \alpha =\text{ }21.239{{\beta }^{-\text{ }7.396\text{ }7}}; $

(3)

$ \beta =0.624\ 3+\frac{27.346}{{{P}_{{{\text{d}}_{0}}}}}; $

(4)

$ {{P}_{{{d}_{0}}}}=\frac{1}{N}\sum\limits_{i=1}^{N}{\sum\limits_{k=1}^{12}{\sum\limits_{J=1}^{M}{P_{{{\text{d}}_{ikj}}}^{{}}}}}; $

(5)

$ R=\sum\limits_{k=1}^{12}{{{R}_{k}}}。$

(6)

式中：R为降雨侵蚀性，MJ·mm/(hm²·h·a)；N为计算数据序列长度，次；P_dijk为第i年第k月第j次侵蚀性降雨量，mm；α和β为模型参数，通过式(3)和(4)估算；P_d0为侵蚀性降雨的多年平均值，mm；M为第i年第k月侵蚀性降雨的回数(日降雨量≥12 mm作为侵蚀性降雨)；R_k为第k月的降雨侵蚀性，MJ·mm/(hm²·h·a)。感谢四川省九寨沟管理局提供的诺日郎雨量站和九寨天堂雨量站2016至2018的降雨数据。通过以上公式计算出诺日郎雨量站和九寨天堂雨量站的降雨可蚀力，再通过ArcGIS中的插值方法得到研究区R值栅格图。由于地震对土壤可蚀性、植被覆盖度、地形和水土保持措施有影响，对降雨条件没有影响。笔者在于探讨地震对研究区土壤侵蚀的影响，为了排除降雨侵蚀因子对研究结果的影响，假设震前震后降雨侵蚀力相同。

3.2 土壤可蚀性因子

研究区共发育了14种土壤，其中：山地石灰土所占面积最大为25.48%；山地暗棕壤和山地棕壤共占总面积的34%左右；亚高山草甸土和高山草甸土共占24%左右^[15]。一般土壤流失方程中的土壤可蚀性因子是指标准地块中每单位降雨可蚀性系数的土壤流失量。目前国际上估计K值通用的方法是Wischmeier等^[16]的诺谟图。国内有很多学者曾直接引用诺谟方程计算K值，评价我国不同地区的土壤侵蚀情况；然而，直接引用该方法得出的值与实验实测值相差较大。卜兆宏等^[17]研究了运用诺谟方程与查图表法得出的K值其之间的差异，给出了有机质土壤的修正系数(表 1)。

笔者根据文献[15]和全国2次普查结果和四川土种志，我们得出研究区各类土壤的相关属性。根据不同类型土壤的有机质质量分数和土壤质地，获得不同土壤类型的K值。由于资料难获取，只有震前土壤图。在本研究中，震前、震后统一使用震前土壤可蚀因子。

3.3 坡长坡度因子提取

坡越长，坡度越大，坡上水流速度和冲刷力越强，侵蚀力越强；但由于空间数据不确定性的缘故，通过DEM直接提取坡度务必会带来误差。汤国安等^[18]研究得出5 m精度的DEM是保证地形描述精度的误差较小分辨率尺度。关于提取坡度的方法，陈楠等^[19]对比了6种不同的提取坡度的方法，最后得出三阶反距离平方差分算法和三阶反距离权差分法误差较小。为了确保提取的地形数据误差影响最小，基于5 m分辨率的地形数据，采用三阶反距离平方权差分算法, 提取研究区区域的地面坡度。笔者运用的土壤侵蚀方程中，把径流起点到中断点或截留点的水平距离称之为坡长。坡长因子是指在其他条件相同的情况下，单位面积任意坡度的土壤流失量比上单位面积标准坡长的土壤流失量的值。其计算公式^[20]为

$ L=\text{ }{{(\lambda /22.1)}^{m}}。$

(7)

式中：L为根据22.1 m标准坡长上转化得到的坡长因子；λ为坡长，m，λ取值如下：

$ λ=af。$

(8)

式中：f为上坡来水流入流像元的总数，个：a为像元边长，m。

$ m = \left\{ {\begin{array}{*{20}{l}} {0.2,}&{\theta < 1\% ;}\\ {0.3,}&{1\% \le \theta < 3\% ;}\\ {0.4,}&{3\% \le \theta < 5\% ;}\\ {0.5,}&{\theta \ge 5\% }。\end{array}} \right. $

(9)

式中m为坡长因子指数。

由于美国通用土壤侵蚀方程中的坡度因子S计算关系式是在缓坡条件下通过径流小区实验测得，USLE模型所允许的最大坡度为18°，研究区陡坡处不适用该方程，考虑到研究区的坡度情况，运用分段的方式计算坡度因子，缓坡采用通用土壤侵蚀方程中的坡度因子计算关系式，陡坡采用Liu Baoyuan等^[21]的计算公式：

$ S=\left\{ \begin{array}{*{35}{l}} 10.80\sin \theta +0.03, & \theta <{{5}^{{}^\circ }}; \\ 16.80\sin \theta -0.5, & {{5}^{{}^\circ }}\le \theta <{{10}^{{}^\circ }}; \\ 21.91\sin \theta -0.96, & \theta \ge {{10}^{{}^\circ }} 。\\ \end{array} \right. $

(10)

式中：S为坡度因子，量纲一；θ为坡度，(°)。

求出坡度后带入式(11)计算出S因子；通过水文分析板块提取像元总数和像元边长，使用栅格数据叠加分析并结合式(7)得出L因子。

3.4 水土保持因子

水土保持因子是采取特殊措施后的水土流失量比上沿坡种植时的水土流失量的值。再无任何水土保持措施的水土保持因子为1，将地震造成的崩塌、泥石流、滑坡、碎屑流等视为无任何水土保持措施的特殊侵蚀类型，将其赋值为1，其他情况P值在0~1之间。其他情况的P值笔者参考周正朝等^[22]的土壤侵蚀模型研究进行赋值。由于研究区的地形地貌较为复杂，况且地震诱发的地质灾害造成地物更加支离破碎加大识别难度，采用计算机自动分类的方法造成的误差较大；故震后信息提取通过震后高分2号2.5 m分辨率的影像采用人工目视解译的方法提取数据。对震前、震后不同水土保持因子赋值(表 3)。

3.5 植被覆盖因子

植被覆盖因子是指在同一土壤、同一坡度、同一降雨条件下，特定作物或植被的土壤流失量与未耕作过的土地的土壤流失量的比，比值为0~1。在实践中，土壤流失与植被覆盖度的关系逐渐被大多数学者所接受。由蔡崇法等^[23]研究可知，植被覆盖度与植被覆盖因子之间相关性很大，植被覆盖因子可以按照如下公式计算：

$ \begin{array}{l} C = \\ \left\{ {\begin{array}{*{20}{c}} {1,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}&{c = 0;}\\ {0.6508 - 0.3436{{\log }_{10}}\left( {c \times 100} \right),}&{0 < c < 78.3\% ;}\\ {0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}&{c \ge 78.3\% }。\end{array}} \right. \end{array} $

(11)

式中：C为植被覆盖因子；c为植被覆盖度，%。

利用研究区的I计算c值，其公式为

$ c=\left( I-{{I}_{\min }} \right)/\left( {{I}_{\max }}-{{I}_{\min }} \right)。$

(12)

式中：I为归一化植被覆盖指数；I_max为研究区的最大值，I_min为研究区的最小值。

利用震前震后的TM遥感影像在ERDAS中提取I，计算研究区的归一化植被覆盖指数，利用式(12)计算研究区植被覆盖度。当植被覆盖度>78.3%时，地表的侵蚀量以可忽略不计，当植被覆盖度为0时，C赋值为1，当0 < 植被覆盖度 < 78.3%时，C值根据式(11)计算。

4 结果与分析 4.1 土壤侵蚀空间分布

笔者参考SL 190—2007《土壤侵蚀分类分级标准》，按照水利侵蚀强度进行土壤侵蚀强度分级，其分类标准如表 4。

各因子在统一的坐标系统下生成栅格大小为5 m×5 m的栅格图。在ArcGIS栅格计算器中将九寨沟震前、震后各个因子进行叠加相乘，得到研究区震前震后土壤侵蚀强度空间分布图(图 3)。九寨沟地震对九寨沟县产生了很大的影响，以震前的资料做为参照，震后共解译出同震地质灾害1 883处，面积9.7万~23.1万m²，植被破坏面积达到8.11 km²。地震对当地的水土保持造成了很大的破坏，和震前相比震后有大量的土壤和岩石物质释放到与陡峭的斜坡相连的河床中，这些斜坡是剧烈土壤侵蚀的能源，在这次地震后的几年里，这些区域将持续保持剧烈的土壤侵蚀。运用相同的方法定量计算出震前、震后的土壤侵蚀情况(表 5)。

为了验证震后精度误差，采用常规分析方法进行分析，误差公式为

$ E={{A}_{\text{g}}}-{{A}_{\text{v}}}。$

(13)

式中：A_g为栅格求算的面积，m²；A_v为用矢量数据求出的面积，m²；E为面积误差，m²：为正时表示比实际面积大，为负时情况相反。

面积误差率计算公式为

$ Y=E\text{ }/\text{ }{{A}_{\text{v}}}。$

(14)

式中：Y为误差率，%；Y为正时表示比实际面积大，为负时情况相反。

震后各侵蚀等级的面积误差和误差率如表 6所示。

各个侵蚀等级的误差率在2%左右，在误差允许范围之内。

从表 5可以得出：研究区震前震后的侵蚀强度跨越了轻微度到剧烈6个等级，震前研究区的侵蚀面积为752 km²，土壤侵蚀总量为122万9 262 t/a，平均侵蚀模数为1 634.65 t/(km²·a)；震后侵蚀面积为807 km², 土壤侵蚀总量为193万5 785 t/a，侵蚀强度为2 378.11 t/(km²·a)。侵蚀参数变化来看，地震后侵蚀面积增加55 km²，侵蚀总量提高57%，侵蚀强度增强743.46 t/(km²·a)。由表 5可知震前、震后的侵蚀等级均为轻度侵蚀，研究区震前、震后土壤侵蚀以强度侵蚀及以下侵蚀等级为主，占到全区域土壤侵蚀的90%以上。在侵蚀强度类别中轻微度、轻度和中度的面积占比在震后全部下降，强度侵蚀及以上侵蚀的面积占比都有所上升，说明在地震的作用下，区域的土壤侵蚀强度等级有向强度侵蚀等级及以上侵蚀等级转移的趋势。为了更直观地看出震后土壤侵蚀变化的空间分布规律，将震前和震后的土壤侵蚀强度空间分布图进行相减，得到震后土壤侵蚀强度变化图。根据上文得出的侵蚀规律，将震后土壤侵蚀强度变化图以强度增量平均值为界划分为易受地震影响区(以下简称易影响区)和不易受地震影响区(不易影响区)，如图 4所示。

4.2 震后易影响区的空间分布规律 4.2.1 易影响区的空间分布规律与距发震断层距离的关系

九寨沟地震地区的地质构造较为复杂，根据李渝生等^[24]在野外发现多条北西—南东走向的断层，断层分布见图 5。利用ArcGIS的空间分析功能，以见图 5中的发震断层为基础，每1 km生成一个缓冲区，然后算出每个缓冲区内易影响区的面积和侵蚀强度，以得出易影响区的面积和侵蚀强度与距断层距离之间的关系。如图 6所示，在距发震断层2 km范围内，震后易影响区的面积较大，随着距断层距离逐渐增大，易影响区的面积也逐渐减小，在距发震断层6 km以外，几乎没有易影响区。且发现易影响区的侵蚀强度相比震前同区域的侵蚀强度也呈现一定的规律，在易影响区范围内，随着距断层距离逐渐增大，易影响区的侵蚀强度增量急剧减小。

4.2.2 易影响区的空间分布规律与坡度的关系

坡度作为通用水土保持公式中的因子，可想而知，坡度在很大程度影响易影响区的空间分布，笔者将每隔10°生成一个坡度带，最终将坡度划分为8个坡度带。对易影响区的坡度带分布规律进行分析，计算结果如表 7所示。

可以看出，易影响区的土壤侵蚀模数与坡度有很好的相关性，随着坡度的增加，土壤侵蚀模数也随着增加。土壤侵蚀量最大的坡度在40°~50°范围内，侵蚀总量和侵蚀面积分别占到易影响区总量的42.78%和41.50%，在30°~60°坡度带内，侵蚀量和侵蚀面积分别占到易影响区总量的90.60%和89.25%。由此可知，震后易影响区的分布范围主要集中在30°~60°坡度带，尤其是40°~50°坡度带。因此震后加大对该坡度带的水土保持工作是很有必要的。

4.2.3 易影响区的空间分布规律与高程的关系

九寨沟地震发生在四川盆地与青藏高原过渡的高山峡谷地带，地形起伏较大，如表 8所示，土壤侵蚀与高程有很大的相关性。易受地震影响的区域面积虽海拔的升高而增加，达到一定高度后，随海拔的升高而减小，在2 000~3 000 m易受地震影响的区域面积最大为5 km²，在海拔2 000~4 000 m易影响区的面积占到总易影响区的99%。其主要原因在于小于2 000 m的地区，地形较平，坡度较小，并且植被覆盖较好，水土保持措施较好，受地震影响也较小。地震后戴岚欣等^[25]共解译出地质灾害1 883处，地质灾害点主要分布在2 600~3 200 m。地质灾害发生的地方导致山体破碎，土壤松动，增大了土壤侵蚀的可能性。

4.2.4 易影响区的空间分布规律与坡向的关系

一般情况下，坡向可以通过对日照时长和土壤水分的影响从而对土壤侵蚀造成影响。易受影响区每一坡向级的面积和侵蚀总量如表 9所示，北坡和东北坡的易受影响区的面积最大，分别占到易受影响区总面积的24.59%和16.78%，其主要原因在于北坡和东北坡受日照时间短，水资源较少，导致该区域的植被覆盖度少，土壤较松软，且植被对土壤的加固能力弱，在地震来临时，易受到影响。东南坡的占比面积少于东北坡占比面积，但在侵蚀总量上，东南坡的侵蚀总量大于东北坡的侵蚀总量。其原因在于，研究区位于中国的西南部主要受东南季风和西南季风的影响，导致东南坡区域的降雨偏多，最终使土壤侵蚀总量增加。

5 结论

笔者结合GIS和USLE定量计算了震前、震后九寨沟地震灾区的土壤侵蚀量，划分出易受地震影响的侵蚀区域和不易受地震影响的侵蚀区域，并定量分析了易受地震影响的侵蚀区域的空间分布规律与距发震断层、坡度、高程和坡向的关系。本研究可为灾后重建工作提供理论依据与指导，尤其是对震后灾区的水土保持工作提供很好的参考价值。

1) 从侵蚀参数变化来看，相比震前，震后侵蚀面积增加8.1%，侵蚀总量提高28%，侵蚀强度增强18.39%。在侵蚀强度类别中轻微度、轻度和中度的面积占比在震后全部下降，强度侵蚀及以上侵蚀的面积占比都有所上升，因此震后水土保持工作应重视对中度侵蚀及以上的侵蚀区域的治理。

2) 对震后易受地震影响区域的空间分布规律和距断层距离进行分析。结果表明距发震断层距离越近，易受地震影响区域的面积越大，在距发震断层2 km范围内，震后易影响区的面积最大，随着距断层距离逐渐增大，易影响区的面积也逐渐减小，在距发震断层6 km以外，几乎没有易影响区。且发现易影响区的侵蚀强度相比震前同区域的侵蚀强度也呈现一定的规律，在易影响区范围内，随着距断层距离逐渐增大，易受影响区的侵蚀强度增量急剧减小，因此应该加强对距发震断层2 km区域范围的水土保持工作。

3) 研究区易受地震影响的区域主要发生在高程2 000~4 000 m, 占到总面积的99%，在30°~60°坡度带内，侵蚀面积分别占到易影响区总面积的89.25%，易受地震影响的区域在北坡、东北坡和东南坡的侵蚀面积达到61.07%；因此高程在2~4 km、坡度在30°~60°和坡向在北坡、东北坡和东南坡的区域为最易受地震影响的区域，加强该区域的水土保持工作是非常有意义的。