基于GIS/RS的灾区土壤保持功能恢复效应评价

刘桐恺, 陈元哲, 董天翔

(1.南京信息工程大学 应用气象学院, 南京 210044; 2.广东省郁南县气象局, 广东 郁南 527199)

土壤保持在维护区域生态安全中具有重要作用，土壤侵蚀可导致土壤和养分的流失、土地贫瘠、土层变薄、宜耕地减少等问题。将GIS和RS技术与RUSLE相结合，对实现灾区土壤保持功能效应评价具有重要意义。2008年“5·12”汶川地震，引发了严重的地质灾害和生态系统损毁，加剧了土壤侵蚀，增加了山洪暴发的风险，对灾区居民生命财产安全、生态安全带来巨大威胁，因而进行地震灾区土壤保持恢复效应评价，掌握灾区土壤保持功能的恢复效果及变化趋势，是灾区环境保护和有关决策部门迫切需要掌握的重要信息。

近年来，国内外学者对流域土壤侵蚀量的研究层出不穷，从最初的定性评价向着定量化计算推进。其中，将“3S”技术与USLE结合对流域土壤侵蚀量的定量化评估是现阶段研究的焦点。在国外，自美国农业部282号和537号农业手册于1965年和1978年先后刊出由Wischmeier和Smith[1-2]提出的通用水土流失方程(USLE)，其后的定量化研究大多基于上述方程，经后人改进，丰富了各因子的实际意义，形成了RUSLE版本[3]。在我国，研究者以通用水土流失方程为基本形式，根据气候、土壤等具体情况进行相应的修正，建立了适合于我国的土壤流失方程，如刘元宝等[4]参照美国USLE建立了中国土壤流失方程(chinese soil equation, CSLE)。目前，国内学者主要将通用水土流失方程应用在流域[5]，而对受地质灾害影响导致土壤侵蚀极为严重的灾区的研究尚为少见。

本文基于修订的通用水土流失方程(RUSLE)，结合GIS和RS技术，选择四川省极重灾区为研究区，以ArcGIS软件作为平台，遥感数据Landsat5 TM为主要数据源，DEM等其他数据为辅助数据源，建立灾区土壤保持功能评估模型，定量化评价灾区2007年、2009年、2011年地震前后3个时段的土壤侵蚀情况，揭示地震灾区土壤保持功能恢复效果和变化趋势，为制定灾区生态重建和生态恢复规划提供参考。

1 研究区概况

本文的研究区是“5·12”汶川大地震中的汶川县、北川县等10个极重受灾区域，经纬度范围为东经102°51′46″—105°37′12″，北纬30°48′8″—33°1′26″，总面积25 802 km2，其地形复杂，高低悬殊，西高东低的特点特别明显，西部为高原、山地，海拔多在3 000 m以上，东部为盆地、丘陵，海拔多在500～2 000 m。其地域辽阔，土壤类型丰富，垂直分布明显。平原，丘陵主要为水稻土、冲积土、紫色土等，是四川省农作物主要产区，高原、山地依海拔高度分别分布不同土壤，大多数都有利于作物的生长，该土壤内富含钾、钙、磷、镁、锰、铁等元素，土质风化度低，土壤发育浅，肥力高，是研究区分布面积最广的土壤之一。气候区域表现差异显著，东部盆地大部年降水量900～1 200 mm，但在地域上，盆周多于盆底，盆西缘山地是降雨最多之地，为1 300～1 800 mm，次为盆东北和东南缘山地，为1 200～1 400 mm，盆中丘陵区降雨最少，为800～1 000 mm。在季节上，冬季降水最少，占全年总雨量的3%～5%，夏季降水最多，占全年总雨量的80%，冬干夏雨，雨热同期。高原降雨少，年降水量大部为600～700 mm，6—9月为雨季，降雨占全年总雨量的70%～90%，11—4月为干季，各月降水量小于10 mm，西南山地降水地区差异大，干湿季节分明，大部分地区年降水800～1 200 mm，木里以北与川西北高原接壤，年降水小于800 mm，安宁河东侧与东部盆地相当，年降水1 000 mm左右。雨季(6—9月)降水占全年总降水量的85%～90%，气象灾害种类多，发生频率高，范围大，是汶川地震中受灾最为严重的几个区域，可以看出具有研究的代表性。

2 数据获取与遥感专题信息提取

根据研究的实际情况和需要，本次建模所用数据源分为两类：遥感数据源与非遥感数据源。遥感数据源选用2007年、2009年和2011年美国陆地卫星Landsat5 TM影像，非遥感数据源主要包括研究区不同年份土壤类型分布数据、DEM(数字高程模型)数据、四川省1∶50 000行政区划图以及气象站降雨信息数据，其中DEM为ASTER-GDEM 30 m数据。

2.1 遥感数据源

Landsat5 TM卫星是美国于1984年3月发射的光学对地观测卫星。它所获得的图像是迄今为止在全球应用最为广泛，成效最为显著的地球资源卫星遥感信息源，同时它也是目前在轨运行时间最长的光学遥感卫星。它的轨道高度为705 km，轨道倾角为98.2°，运行周期为98.9 min，刈幅宽度为185 km，重访周期为16 d。本文主要以Landsat5 TM数据为遥感影像数据源，由于研究区范围比较大，故每个年份都需要6景影像才能完全覆盖整个研究区域，各时段具体数据信息见表1。

表1 遥感影像统计详情

注：数据申请来源于地球系统科学数据共享平台。

2.2 非遥感数据源

为了与Landsat5-TM保证具有一致的空间分辨率，本文选用的DEM数据为利用ASTER生产的30 m GDEM数据第二版本产品，该版本产品相比第一版本精度得到了显著提高[6]。ASTER GDEM数据是世界上迄今为止可为用户提供的最完整的全球数字高程数据，它填补了航天飞机测绘数据中的许多空白，为各学科、领域的用户和研究人员提供了需要的高程和地形信息，这些信息对整个地球科学都具有很大价值。该数据可免费在地理空间数据云(http:∥www.gscloud.cn/)下载，本文中对GDEM数据的处理主要包括影像拼接和数据裁剪，2007年灾区DEM数据如图1所示。本文中另一主要非遥感数据源为中国气象数据共享服务网(http:∥cdc.cma.gov.cn)提供的降雨量月值数据和各气象站的经纬度数据，时间范围为2007年、2009年和2011年各个月份，为了获取与遥感数据相匹配的栅格化数据，本文采用GIS空间插值方法，根据气象站点的经纬度信息，通过对离散的气息数据进行克里金插值，获得空间上连续并且像元大小与遥感影像一致、投影相同的气象栅格图。

图1 灾区DEM数据

2.3 遥感数据处理与信息提取

文中对3个年份遥感数据的处理过程主要包括遥感影像辐射校正、几何校正、裁切、镶嵌、专题信息提取等[7-8]，其中专题信息提取过程中，鉴于本文的研究目的，需要建立研究区土地利用分类体系，参考国家颁布的“土地利用分类标准”[9]，考虑到研究区实际情况，制定乔木、园地、建设用地、未利用地、水体、灌丛、耕地和草地8个类别的土地利用分类体系。

上述土地利用分类体系，主要借助Google Earth上的高分辨率卫星影像，在TM影像上建立训练样区，综合考虑分类时间和分类精度等因素，每类地物均在影像上均匀选取50个感兴趣区，运用最大似然分析算法进行监督分类，最后对分类后的结果进行后处理和精度评定，详细的专题信息提取方法和精度验证可参考文献[10]，本文得到3个年份的土地利用分类总体精度分别为77.34%，82.59%，83.61%，其精度符合模型设计要求，图2分别为不同年份分类结果图。

图2 研究区不同年份土地利用分类

3 空间数据建模与分析

3.1 土壤保持功能评估总体技术路线

土壤保持功能效应恢复评估建模由降雨侵蚀力因子R值、土壤可蚀性因子K值、植被覆盖因子C值、土壤保持措施因子P值、坡度坡长因子LS值5大影响因素构成。借鉴国内外土壤侵蚀定量评估理论，建立各个因素的定量评估模型，将所有空间数据统一到WGS_1984_UTM_Zone_48N空间参考坐标下，统一栅格像元大小为100 m。在此基础上，结合修订的通用水土流失方程(RUSLE)，利用ArcGIS空间分析工具，得到地震前后3个时段灾区土壤侵蚀情况，进而进行对比分析土壤保持功能恢复效果和变化趋势。

3.2 土壤保持功能指标分析

3.2.1 降雨侵蚀力因子 降雨侵蚀力是降雨引起土壤侵蚀的潜在表现，它与降雨量、降雨时长、降雨强度、降雨动能、雨滴大小、瞬时雨率等有关，年降雨侵蚀因子(R)的计算方法很多[11-13]，但尚无通用的计算方法，一般将计算EI得到的降雨侵蚀力称为经典方法，经典算法需降雨过程资料，工作难度大，且资料难以获取。因而，多采用其他资料代替降雨动能计算降雨侵蚀力，这种方法称为简易算法。本文根据灾区实际情况选择刘秉政[12]的简易计算公式：

(1)

式中：R为降雨侵蚀力因子[(MJ·mm)/(hm2·h·a)]；P6-9为年6—9月降雨量之和(mm)；Pre为年雨量(mm)。

获得不同年份气象站点降雨数据后，按照图4建立降雨侵蚀力评估模型，图3为3个时段灾区降雨侵蚀力因子空间分布图。

图3 研究区不同年份降雨侵蚀力因子分布

3.2.2 土壤可蚀性因子 土壤可蚀性因子(K)值的大小表示土壤对侵蚀的敏感性及降水所产生的径流量与径流速率的大小，反映土壤对侵蚀外营力(降雨、径流、地形、地表植被和人为活动)剥蚀和搬运的敏感性，是影响土壤流失量的内在因素[14]，土壤结构、有机质含量和土壤剖面渗透性的影响[15]，尤其与土壤机械组成和土壤有机质含量的相关性较高[16-17]，影响K因子的有多方面，但一般说来，质地越粗或越细的土壤有较低K值，而质地适中的反而有较高的K值。估算K值的方法很多，通常根据实测的E值，应用通用水土流失方程反推获取K值,但获取大面积的实测E值是不可能的，本文K值采用Shirazi等[18]在EPIC(erosion productivityimpact calculator)模型中的计算公式：

(2)

式中：K为土壤可蚀性因子[(t·hm2·h)/(MJ·mm·hm2)]；SAN为砂粒含量(%)；SIL为粉粒含量(%)；CLA为黏粒含量(%)；C′为有机碳含量(%)；SN1=1-SAN/100。图4为3个时段灾区土壤可蚀性因子空间分布图。

3.2.3 植被覆盖因子 植被覆盖因子(C)值反映的是所有覆盖和管理变量对土壤侵蚀的综合作用，研究表明，C因子要受到诸如植被、作物种植顺序、生产力水平、生长季长短、栽培措施、作物残余管理、降雨时间分布等众多因素控制。C因子作为度量土壤侵蚀量的一个重要参数，其值介于0～1之间。基本没有土壤侵蚀危险的地区被赋予0，1值被赋予给那些很容易受到侵蚀的地区，归一化植被指数(NDVI)是估算C因子最普遍的数据，国内外不少研究者建立了植被指数和C的关系式[18-19]，本文采用Jong[20]在1994年得出的由NDVI计算USLE-C的公式(不包括水体和建筑用地等非植被区域)：

(3)

式中：C为地表植被覆盖因子(无量纲)；α和β是决定NDVI-C关系曲线的参数，由于建筑用地等的土壤侵蚀强度很小，因而将其C值赋0。

图5为3个时段灾区植被覆盖因子空间分布图。

图4 研究区不同年份土壤侵蚀力因子分布

图5 研究区不同年份植被覆盖因子分布

3.2.4 土壤保持措施因子 土壤保持措施因子(P)值是指特定水土保持措施的土壤流失量与相应未实施保持措施的顺坡耕作地块的土壤流失量之比值，通常的侵蚀控制措施有：等高耕作、修梯田等。措施好的P值小，反之P值大，P值介于0～1之间。0值代表根本不发生侵蚀的地区，而1值代表了未采取任何控制措施的地区,参照有关研究结果[21-23]，本文中P值按照不同土地利用类型赋值，其中非耕地区域认为是未采取任何土壤保持措施赋值为1，对于耕地，通常坡度越大，土壤保持作用越突出，因此耕地依据表2按照坡度范围赋值。

表2 不同坡度范围耕地的P值

按照技术流程，建模过程将耕地区域与非耕地区域区分开，分别按照不同的形式赋值，最后再将二者合并得出整个研究区域的土壤保持措施因子空间分布(图6)。

图6 研究区不同年份土壤保持措施因子分布

3.2.5 坡度坡长因子 坡度因子S是指在其他条件相同的情况下，任意坡度下的单位面积土壤流失量与标准小区坡度下单位面积土壤流失量之比。美国通用土壤流失方程中的坡度因子S计算关系式，是以3%～18%缓坡条件下天然径流小观测资料而建立的[24]。通用土壤流失方程允许计算的最大坡度为18%(约10°)，考虑到研究区的坡度情况，本文S值采取Liu等[25]的坡度因子计算公式：

(4)

式中：S为坡度因子；θ表示坡度角。

坡长因子L是指在其他条件相同的情况下，任意坡长的单位面积土壤流失量与标准坡长单位面积土壤流失量之比[26]，通过下式从DEM中提取。

L=(λ/22.13)m

(5)

式中：L为坡长因子；λ为水平坡长(m)，22.13是标准小区坡长(m)。由于地震对于研究区的整体DEM影响并不大再加上数据更新的限制，本文中坡度坡长因子均使用震前坡度坡长因子，计算结果如图7所示。

图7 研究区不同年份坡度坡长因子分布

3.3 灾区土壤保持功能评估模型总体构建

在建立以上各个因子的子模型后，通过评价地震前后灾区平均侵蚀模数的变化反映土壤保持生态功能恢复效应。根据修订的通用土壤流失方程(RUSLE)[3]，土壤平均侵蚀模数的计算公式如下：

A=R×K×LS×C×P

(6)

式中：R为降水侵蚀指标[(MJ·mm)/(hm2·h·a)]；K为土壤可蚀性因子[(t·h)/(MJ·mm)]；LS为坡度坡长因子(无量纲)；C为地表植被覆盖因子(无量纲)；P为土壤保持措施因子(无量纲)。

按照上述原理与公式，在ArcGIS中运用空间建模工具模块，构建灾区土壤保持功能评估总体模型，并将各年份数据分别置于模型中，获取研究区土壤平均侵蚀模数，对结果进行统计分析可知，总面积2.58万km2的研究区内，2007年区域土壤侵蚀模数平均为3 389.57 t/(hm2·a)，年侵蚀总量为8.75×107t，2009年土壤侵蚀模数平均为5 580.60 t/(hm2·a)，年侵蚀总量为1.44亿t,增加了64.64%，2011年土壤侵蚀模数平均为4 236.97 t/(hm2·a)，年侵蚀总量为1.09亿t，相比2008年减少了24.08%，土壤保持功能有所恢复。参考中国土壤侵蚀强度分类分级标准[27],按以下规则进行土壤侵蚀强度等级划分：Ⅰ级微度侵蚀平均侵蚀模数小于2 t/(hm2·a)；Ⅱ级轻度侵蚀平均侵蚀模数2～25 t/(hm2·a)；Ⅲ级中度侵蚀平均侵蚀模数25～50 t/(hm2·a)；Ⅳ级强度侵蚀平均侵蚀模数50～80 t/(hm2·a)；Ⅴ级极强度侵蚀平均侵蚀模数80～150 t/(hm2·a)；Ⅵ级剧烈侵蚀平均侵蚀模数大于150 t/(hm2·a)。对灾区不同年份平均土壤侵蚀模数进行分级(图8)并统计各等级侵蚀面积与侵蚀模数(表3)。

从土壤侵蚀面积来看，2007—2009年的破坏期，由于地震的影响导致微度侵蚀、极强度侵蚀、剧烈侵蚀所占面积比例在增加，增幅分别为33.25%，74.93%，277.87%，中度侵蚀和强度侵蚀面积所占比例分别减少了6.23%，14.76%，面积分别减少了378 607，85 855 hm2，而从2009—2011年的恢复期，微度侵蚀、极强度侵蚀、剧烈侵蚀的面积所占比例有所减少，分别减少9.96%，2.79%，7.75%，减少了255 650，71 551，198 899 hm2，相比较破坏期，中度侵蚀和强度侵蚀所占面积有所增加，分别增加了1.57倍和0.73倍。可见，2007—2009年灾区土壤侵蚀从侵蚀范围上呈现中度侵蚀、强度侵蚀向微度侵蚀、极强度侵蚀和剧烈侵蚀转移的迹象，地震导致灾区土壤侵蚀加剧，土壤保持功能遭到削弱，而恢复期则由微度侵蚀、极强度侵蚀、剧烈侵蚀向中度侵蚀和强度侵蚀转移，土壤保持功能得到一定程度的恢复。

图8 研究区不同年份土壤侵蚀强度分级

侵蚀等级2007年侵蚀面积比例侵蚀模数比例2009年侵蚀面积比例侵蚀模数比例2011年侵蚀面积比例侵蚀模数比例微度侵蚀27.190.1836.230.1826.270.13轻度侵蚀26.7910.4120.563.6021.516.49中度侵蚀23.3624.108.605.5522.1018.56强度侵蚀11.6820.858.349.4214.3920.72极强度侵蚀7.5022.7213.1225.5510.3325.09剧烈侵蚀3.4821.7413.1555.705.4029.01

从土壤侵蚀模数上来看(图9)，从2007—2009年，地震导致极强度侵蚀和剧烈侵蚀的侵蚀模数发生急剧增加，中度侵蚀和强度侵蚀侵蚀量在减少，而与破坏期相反的是从2009—2011年，中度侵蚀和强度侵蚀的侵蚀模数有所增加，极强度侵蚀和剧烈侵蚀的侵蚀量在减少，这说明地震过后两年灾区土壤保持功能得到恢复。

图9 地震前后不同年份不同侵蚀强度土壤侵蚀模数

图10为研究区不同区域的土壤保持功能恢复图，可以明显看出“5·12”大地震的第一震中汶川县和第二震中北川县二者土壤侵蚀跳跃比较厉害，并且相对来说恢复的速度也不如其他地方，这说明，离地震中心越近的地方，土壤保持功能受到的削弱程度越大，其恢复能力也越差，因而，建议在进行灾后重建的工作中，应当重点加强对震中及其周边地区的改善。

图10 研究区不同区域土壤保持功能恢复图

4 结 论

(1) 灾区2007年土壤侵蚀模数平均为3 389.57 t/(hm2·a)，年侵蚀总量为8.75×107t，2009年土壤侵蚀模数平均为5 580.60 t/(hm2·a)，年侵蚀总量为1.44亿t,增加了64.64%，2011年土壤侵蚀模数平均为4 236.97 t/(hm2·a)，年侵蚀总量为1.09亿t，相比2009年减少了24.08%，土壤保持功能有所恢复。

(2) 从土壤侵蚀面积和土壤侵蚀模数两个角度来看，2007—2009年灾区土壤侵蚀呈现中度侵蚀、强度侵蚀向微度侵蚀、极强度侵蚀和剧烈侵蚀转移的现象，地震导致灾区土壤侵蚀加剧，土壤保持功能遭到削弱，而恢复期则由微度侵蚀、极强度侵蚀、剧烈侵蚀向中度侵蚀和强度侵蚀转移，土壤保持功能得到一定程度的恢复。

(3) 离地震中心越近的地方，土壤保持功能受到的削弱程度越大，其恢复能力也越差，因而，有关部门在进行灾后重建的工作中，应当重点加强对震中及其周边地区的改善。

本文基于GIS和RS技术，结合遥感、降雨等实测数据，对由地质灾害诱发土壤侵蚀极为严重的灾区进行了创新性研究，对比地震前后3个阶段，从理论层面估算了研究区土壤保持功能恢复状况，揭示了研究区灾后土壤保持功能恢复情况和变化趋势，为相关部门制定灾区生态重建和生态恢复规划提供参考，具有重要意义。