土壤侵蚀是指地表土壤及母质在自然或人为因素影响下, 被破坏、剥蚀、搬运和沉积的过程^[1]。目前已成为全球主要的生态环境问题, 不仅会引起土地退化、土壤肥力下降, 诱发泥石流、滑坡等自然灾害, 而且还会造成河流泥沙淤积, 导致河床抬高, 严重影响河流的自然泄洪能力及水利水电工程的建设和使用^[2-3]。据统计, 世界水库每年因沉积物造成的库容损失约0.5%—1%, 蓄水量将会在50年内减少到目前的一半左右^[2]。

湄公河长度和沉积物载荷量位居东南亚首位, 年径流量高达4.75×10¹¹ km^3[3]。流域拥有丰富的自然生态系统, 为沿岸居民提供了食物、交通和经济等众多方面支持^[4]。流域地跨多个纬度, 海拔差异较大, 地形复杂多样, 河段多狭窄陡峭, 具有丰富的水电开发潜能, 但部分国家面临着水电开发技术较为落后和资金短缺等问题。中国电力企业具有先进的开发技术和丰富的建设经验, 且可提供充足的资金支持, 在促成中国与中国-中南半岛经济走廊内各国开展水电贸易合作方面发挥了重要作用^[5]。该经济走廊是中国参与水电项目投资建设的主要区域。统计发现中国在海外参与投资建设的水电项目83个, 近40%位于该经济走廊。水利水电项目施工阶段会改造地形, 破坏地表覆盖, 增加土壤侵蚀风险^[6]。此外, 大湄公河次区域经济合作机制和中国-中南半岛经济走廊项目的开展, 加快了流域内大坝建设速度, 目前已完成和正在建设133座, 一旦全部完工, 96%的向三角洲输送的河流沉积物将会被困在河流中, 严重威胁流域居民的生存和生态系统的平衡^[7]。因此, 对湄公河流域土壤侵蚀状况进行分析具有重要意义。

优先级理论在水土流失治理中起至关重要的作用^[8]。研究表明, 多数流域的大部分沉积物是由关键区域产生的, 仅选定这些关键区域实施水土流失治理, 就会大大降低全流域侵蚀量^[9]。因此, 确定湄公河各子流域水土流失治理的优先顺序, 对实施具有针对性和成本效益最大化的水土保持措施具有直接指导意义。本文在气候、土壤、遥感等数据的基础上, 运用GIS技术和USLE模型, 得到湄公河流域土壤侵蚀强度及其空间分布特征。通过联合信息熵方法计算土壤侵蚀与影响因子的关系, 并基于优先级理论确定子流域的水土流失优先治理级别。

1 研究区概况

湄公河发源于中国, 流经中国、缅甸、泰国、老挝、柬埔寨和越南6个国家(图 1)。河流全长4.9×10³ km, 流域面积达7.95×10⁵ km², 年均径流量为1.5×10⁴ m³/s, 年均沉积物负荷量为1.6×10⁸ t^[10]。其中永珍以上为上游, 永珍到巴色为中游, 巴色以下为下游。流域主要受两种季风影响：5—10月初受西南季风影响, 潮湿多雨, 为雨季；11月—次年3月中旬受东北季风影响, 干燥少雨, 为旱季^[11]。流域地处8°—34°N, 94°—109°E, 南北横跨二十多个纬度, 海拔差约6 km。气温和降水由南到北递增, 年均温最低为-15℃、最高为29℃, 年累降水量最低为440 mm、最高为3500 mm。

2 研究方法 2.1 子流域提取

针对湄公河干流和主要一级支流进行子流域提取。基于ArcGIS中水文模块, 以一级支流于干流的汇流点为泻流点, 提取干流和各一级支流的不嵌套子流域。由于湄公河流经6个国家, 干流及各支流在各国命名存在差异, 为便于研究, 文章对湄公河流域各支流及干流各段所对应子流域进行重新编号。研究区共33个子流域, 编号依照从上至下, 从左至右的顺序：干流13个(M1—M13), 一级支流20个(T1—T20)。

2.2 通用土壤流失方程模型

文章选用通用土壤流失模型(USLE)进行土壤侵蚀评估, 模型是Wischmeier和Smith于1965年根据落基山以东地区收集到的近万个径流小区的试验资料, 进行系统统计分析提出的经验性土壤侵蚀模型, 全面考虑了降雨、地形和土壤可蚀性等因素^[12-14]。方程如下：

(1)

式中, A为土壤侵蚀量, 单位t hm^-2 a^-1, R为降雨侵蚀力因子, 单位MJ mm hm^-2 a^-1, K为土壤可侵蚀性因子, 单位t hm² h MJ^-1 hm^-2 mm^-1, L和S为坡长、坡度因子, C为植被覆盖管理措施因子, P为土壤保持措施因子, L、S、C、P皆为无量纲数据。

2.2.1 降雨侵蚀力因子(R)

R因子反映了由降雨引起的潜在土壤侵蚀能力, 是评价一个地区土壤侵蚀状况的主要动力指标。本文采用Wischmeier月尺度计算公式^[12]：

(2)

式中, P_i为月降雨量(mm), P为年降雨量(mm)。降水数据来自WorldClim, 是Robert等人基于ANUSPILIN插值法, 整合了不同尺度气象数据建立的数据集^[15], 分辨率为30″。

2.2.2 土壤可蚀性因子(K)

K因子反映了土壤被侵蚀营力分离、搬运的难易程度, 不同的土壤类型受影响程度不同。本文采用Sharply等在EPIC模型中的提取K因子计算公式^[16]：

(3)

式中, SAND为土壤砂粒含量, SILT为土壤粉粒含量, CLAY为土壤粘粒含量, C为土壤有机碳含量, SN1=1-SAND/100。土壤属性数据来自全球连续土壤数据集SoilGrids, 该数据集是Hengl等在HWSD数据的基础上建立的, 包含土壤有机碳含量和pH值等, 分辨率为250 m^[17]。

2.2.3 坡度-坡长因子(S-L)

S和L因子反映了地形对土壤侵蚀的影响。本文借助ArcGIS 10.3的表面分析和水文分析模块, 基于GTOPO30 DEM数据, 计算得到坡度数据。坡长λ为从坡面径流的起点开始, 垂直于等高线沿坡向下直到坡度减缓到足以发生沉积的水平距离。根据坡长和坡度间的负相关关系^[18], 坡长取值如下表所示：

结合USLE模型中地形因子方程和Liu等的坡度计算公式^[19], 得出地形因子计算方程：

(4)

(5)

式中, θ为坡度值, m为坡长指数。

2.2.4 植被覆盖管理措施与土壤保持措施因子(C-P)

C因子是指在相同的降雨、土壤和地形条件下, 有植被覆盖的坡面与连续休闲裸露坡面的土壤流失量之比, 值介于0—1之间^{[12, 20]}。由于不同植被类型和管理方式对土壤侵蚀的作用不同, 本文将土地利用类型分为林地、水浇地等10大类, C因子值参考亚洲和越南土壤侵蚀资料制定^[22-23]。P因子是指特定保持措施下与标准条件下的土壤流失量比值, 值域范围为0—1。基于前人的研究, 结合研究区的土地利用类型特点进行P因子赋值^[21-24]。

土地利用数据为GlobCover 2009, 该数据是欧洲太空局和联合国环境规划署等共同参与完成, 数据运用LCCS分类系统, 包含22种土地利用类型, 分辨率为300 m。

2.3 信息熵计算方法

空间上连续而状态属性离散的地理数据称为离散空间场, 离散空间场之间往往隐含着不同程度的相关, 需对其进行定量计算和表达, 从而得到相关规律^[28-29]。本文利用信息熵计算方法, 计算土壤侵蚀与各因子的空间相关性。信息熵是指对一个随机变量的信息量或不确定性的度量, 可用来计算空间离散面状分布变量之间的相关性^[25], 公式如下：

(10)

式中, p(x_i)为发生事件x_i的概率, n为可能发生的事件总数。

当x, y构成二维随机变量时, 其联合分布概率为p(x_i, y_i), 则x, y的联合熵为^[26]：

(11)

为便于计算, 常用指标K表示x, y之间的相关程度^[27], 公式为：

(12)

K的范围为[0, 1], K值越大表示x和y的相关程度越高, K=0则表示不相关。

2.4 优先级理论

根据土壤侵蚀严重程度对子流域的治理次序进行排序, 以期达到成本效益最大化^{[8, 28]}。本文将33个子流域水土流失治理次序进行排序, 并分为四个等级：第一级, 子流域平均土壤侵蚀模数为0—500 t km^-2 a^-1；第二级, 子流域平均土壤侵蚀模数为500—2500 t km^-2 a^-1；第三级, 子流域平均土壤侵蚀模数为2500—5000 t km^-2 a^-1；第四级, 子流域平均土壤侵蚀模数为5000—800 t km^-2 a^-1。

3 结果与分析 3.1 湄公河土壤侵蚀因子空间分布

计算得到流域各土壤侵蚀因子(图 2)。R因子值域为1129—40828 MJ mm hm^-2 a^-1, 均值为6923 MJ mm hm^-2 a^-1。流域上游R值普遍较低, 中、下游R值相对较高, 尤其是干流M7—M9和支流T6、T8、T20都高于均值。K因子值域为0.02—0.05 t a MJ^-1 mm^-1, 均值为0.03 t a MJ^-1 mm^-1。K因子分布差异较小, 呈中、下、上游依次降低的分布特征。LS因子值域为0.04—11.00, 均值为1.62。下游LS因子值较低, 中、上游值较高。流域中游和下游东部的C因子值相对较低, 而上游和下游西部的C因子值较高, 尤其是干流M1、M3、M9和支流T1、T2、T7、T9、T11—T17。

3.2 湄公河流域土壤侵蚀分布特征

据中国水利部土壤侵蚀分级标准进行分级统计^[1], 流域内年土壤侵蚀总量为1.06×10⁹ t a^-1, 平均土壤侵蚀模数为1.98×10³ t km^-2 a^-1, 属轻度侵蚀。流域内60.09%的区域为微度侵蚀, 侵蚀总量占比3.88%, 主要包括上游中部(M3—M6、T3—T6)、中游(M9、M10、T7—T17)及下游的T19、M12和M13。39.91%的区域存在不同程度侵蚀, 达到强烈侵蚀的面积占比为12.68%, 侵蚀总量占比却高达71.84%, 包括上游部分地区(M4—M6、T4—T6)、中游的M7、M9、T8和T10、下游的T20。

3.3 湄公河流域土壤侵蚀与各影响因子的关系

土壤侵蚀强度是多因子共同作用的结果, 分析其与不同影响因子之间的关系, 对科学开展水土流失治理工作有重要指导作用。本文基于联合信息熵计算方法, 分别计算了土地利用类型、坡度、海拔、降雨侵蚀力和土壤可蚀性与土壤侵蚀强度空间相关性。由下表可知, 湄公河流域土壤侵蚀中, 土地利用类型、坡度和海拔起主导作用。

3.3.1 土地利用类型与土壤侵蚀强度的关系

湄公河流域以旱地、林地和灌丛为主, 面积占比达87.42%。从侵蚀强度分布看, 裸地/稀疏植被各侵蚀强度分布差异不大, 近53%的面积达到强烈侵蚀；灌丛以微度和轻度侵蚀为主, 超过该类型总面积的一半；草地以轻度、微度和中度侵蚀为主, 面积占比达86.66%；其余各类型皆以微度和轻度侵蚀为主, 面积占比都在90%左右。

由平均侵蚀模数可知, 水体、永久冰雪和城市建设用地基本不发生侵蚀。湿地、水浇地和林地为微度侵蚀, 面积占比为36.44%, 侵蚀量占比为1.69%；旱地和草地为轻度侵蚀, 面积占比为38.69%, 侵蚀量占比为24.36%；灌丛为强烈侵蚀, 面积占比为23.75%, 侵蚀量占比高达70.58%；裸地/稀疏植被属于极强烈侵蚀, 面积占比为0.52%, 侵蚀量占比高达3.41%。

3.3.2 坡度与土壤侵蚀强度的关系

据水利部土壤侵蚀分级标准将坡度分级, 流域坡度以0—5°为主, 面积占比高达75.69%, 主要分在湄公河上游的南部地区、中游和下游, 8°以上主要分布在上游北部地区、中游的M7、T8和下游的T20。

平均土壤侵蚀模数随坡度的增加呈先增加后下降趋势, 于8—15°坡度带达到峰值。其中0—5°为轻度侵蚀；5—8°和25°以上为中度侵蚀；8—25°为强烈侵蚀。流域内侵蚀主要发生在0—15°坡度带, 占总侵蚀量的94.82%, 尤其是5—15°, 面积占比为22.75%, 侵蚀量占比高达62.25%。各侵蚀强度在各坡度带分布差异较大(图 4), 0—5°呈递减趋势, 以微度和轻度侵蚀为主, 面积占比为89.03%；5—8°和15—25°皆呈先下降后增加再下降的趋势, 前者以微度和极强烈侵蚀为主, 后者以微度和强烈侵蚀为主；8—15°各等级差异相对较小, 以微度、中度和剧烈侵蚀为主；25°以上区域以微度、强烈和极强烈侵蚀为主。除0—5°坡度带达到强烈侵蚀的面积占比较少(仅为4.55%), 其余坡度带到达强烈侵蚀的面积占比皆高于35%。

3.3.3 海拔与土壤侵蚀强度的关系

据湄公河流域海拔分布状况, 将流域海拔分级(表 7)。流域内海拔以0—2000 m为主, 海拔为0—500 m主要分布在中、下游, 4000 m以上区域主要包括干流M1—M3和支流T1、T2。

平均土壤侵蚀模数随海拔的上升呈先增加后下降的趋势, 于500—1000 m范围达到峰值。其中0—500 m和2000—4000 m为轻度侵蚀, 500—2000 m和4000 m以上为中度侵蚀。土壤侵蚀主要发生在500—2000 m的区域, 面积占比为30.82%, 侵蚀量占比高达61.75%。从各高程带侵蚀强度分布看(图 5), 0—500 m呈递减趋势, 以微度和轻度侵蚀为主；500—1000 m呈先减少后增加再减少的趋势, 以微度和轻度侵蚀为主；1000—2000 m除微度侵蚀占46%的面积占比外, 其余各侵蚀强度差异不大；2000—4000 m范围内的2个高程带皆呈递减趋势, 以微度和轻度侵蚀为主；4000 m以上则呈现先增加后减少的趋势, 以微度和轻度侵蚀为主。

3.4 子流域优先治理次序及空间分布特征

统计各子流域年土壤侵蚀量, 进行子流域土壤侵蚀强度分级(图 6)。同时对33个子流域优先治理次序进行排序及优先治理分级划分(表 8)。

湄公河33个子流域中微度侵蚀9个, 轻度侵蚀12个, 面积占比为26.98%和32.92%, 侵蚀量占比为3.88%和27.10%, 优先治理等级分别定为第四级和第三级。9个子流域为中度侵蚀, 面积占比为33.32%, 侵蚀总量占比高达48.90%, 优先治理等级定为第二级。干流M7和支流T5、T6达到强烈侵蚀强度, 面积占比为6.78%, 侵蚀量占比高达20.12%, 优先治理等级定为第一级。

基于联合信息熵计算第一级的3个子流域的土地利用类型、坡度、海拔、降雨侵蚀力、土壤可蚀性与土壤侵蚀强度的空间相关性。结果表明3个子流域土壤侵蚀的主要影响因素皆为土地利用类型、坡度和海拔。

4 讨论 4.1 湄公河流域土壤侵蚀情况

基于USLE模型, 分析湄公河流域土壤侵蚀空间分布特征, 结果表明湄公河流域土壤侵蚀模数范围为0—5.26×10⁴ t km^-2 a^-1, 平均土壤侵蚀模数为1.98×10³ t km^-2 a^-1, 为轻度侵蚀。Hoang等研究指出该流域土壤侵蚀模数范围为0—4.16×10⁴ t km^-2 a^-1, 与本研究结果相似, 但Suif等研究显示该流域土壤侵蚀模数范围为0—3.22×10³ t km^-2 a^-1, 其最高值仅为本研究结果的1/10, 是由于USLE各因子计算过程中赋值标准的差异以及海拔等数据分辨率差异造成的, 说明土壤侵蚀模数绝对值有效性需进一步规范。我们发现尽管各研究中土壤侵蚀模数绝对值有差异, 但其空间分布规律是相似的, 流域内大约10%的区域为强烈侵蚀, 侵蚀量可达全流域总量的70%, 主要包括上游部分地区(尤其是M4—M6、T4—T6)、中游的M7、T8、M9和T10四个子流域和下游的T20^[33-34]说明湄公河流域土壤侵蚀强烈区域较为集中, 应引起其流经区域的高度重视, 做好水土保持管理, 减少强烈侵蚀的面积和侵蚀量。

4.2 湄公河流域土壤侵蚀的主要因素

基于联合信息熵发现土地利用类型、坡度和海拔是影响湄公河流域土壤侵蚀的主要因素。陈龙等研究指出湄公河流域上游澜沧江段土壤侵蚀与坡度和海拔密切相关, 姚华荣等对湄公河流域上游研究表明土地利用变化对土壤侵蚀变化具有一定影响^[35-36]。研究显示裸地/稀疏植被是研究区侵蚀程度最严重的地区, 较低的植被覆盖使得植被根系固结土壤能力较弱, 较薄的枯枝落叶层也使得降水可以直接冲击地面, 形成较强的雨滴溅蚀和地表径流, 加快该区域的侵蚀。研究还发现在海拔500—2000 m, 坡度8—15°的灌丛地区是土壤侵蚀强烈的区域, 由于该区域的降水较高, 坡度较大, 国内外研究证明, 坡度与侵蚀强度之间有显著正相关关系。相比之下, 低海拔, 低坡度, 具有较大植被盖度(林地、草地和农田)或受人类活动影响地面被大面积固化的建设用地土壤侵蚀较弱^[33-38]。因此, 减少流域内裸地和稀疏植被面积, 适当进行植树造林等水土保持工程能有效缓解流域土壤侵蚀状况。

4.3 湄公河流域土壤侵蚀优先治理区

基于优先级理论对湄公河子流域分析, 33个子流域中M7、T5、T6达到强烈侵蚀, 表明侵蚀量一般主要来自流域的一小部分区域^{[3, 43-44]}。干流M7和支流T5、T6是湄公河水土流失治理中优先重点治理区, 优先治理等级为第一级。通过联合信息熵计算发现土地利用类型和坡度是这些区域土壤侵蚀治理应重点关注因素。因此, 在其水土流失治理过程中, 除了推广新能源, 减少自然林地和灌丛的开发利用外, 应实施适地适树措施, 大力开展防护林和封山育林, 增加林地和灌丛的覆盖率, 同时对旱地实施坡耕地退耕、改梯田平种, 耕作、轮作和培肥为一体化的水保措施, 以降低土壤侵蚀强度和土壤侵蚀量^[41-43]。此外, 政府应加大对条件比较恶劣的居民区的扶持力度, 增强基础设施建设, 减轻土壤侵蚀诱发的滑坡、泥石流等自然灾害造成的损失。

5 结论

湄公河流域平均土壤侵蚀模数为1.98×10³ t km^-2 a^-1, 为轻度侵蚀, 40%区域存在不同程度侵蚀, 部分地区达到强烈侵蚀强度, 主要包括上游部分地区(干流M4—M6和支流T4—T6)、中游的M7、M9、T8、T10和下游的T20子流域, 是未来土壤侵蚀重点治理区域。

湄公河流域土壤侵蚀主要影响因素为土地利用类型、坡度和海拔。其中, 土地利用类型土壤侵蚀模数呈现, 裸地/稀疏植被>灌丛>草地>旱地>林地>水浇地>湿地>水体、永久冰雪和城市建设用地；土壤侵蚀模数随坡度的增加呈先增加后下降趋势, 于8—15°坡度带达到峰值；土壤侵蚀模数随海拔的上升呈先增加后下降趋势, 于500—1000 m范围达到峰值。整体来看, 剧烈程度土壤侵蚀主要发生在海拔为500—2000m, 坡度为8—25°, 土地利用类型为裸地/稀疏植被和灌木的地区。

基于优先级理论分析, 33个子流域可以划分为四个优先治理等级, 其中21个子流域为微度和轻度侵蚀, 优先治理等级为第四级和第三级。干流M7和支流T5、T6达到强烈侵蚀, 优先治理等级为第一级, 是湄公河流域水土流失治理时应优先考虑的区域。