基于CNN网络和遥感数据的土壤侵蚀评估方法

周海宏，雷 磊，姚普及，吴 健，王 良，刘 涛，景 龑，潘晓彤

(1.国网陕西省电力公司，陕西 西安 710048；2.国网陕西省电力公司电力科学研究院，陕西 西安 710100；3.陕西省东庄水利枢纽工程建设有限责任公司，陕西 咸阳 713208)

土壤侵蚀是一个复杂的4阶段动态过程，包括土壤的剥离、破坏、迁移和随后的沉积物沉积[1].侵蚀是一种自然发生的过程，它通过水和风的物理力量或耕作等与农业活动有关的力量，侵蚀农田表层土，从而影响所有地貌.土壤侵蚀既可以是一个渐进的过程，持续相对不被注意，也可能以惊人的速度发生，造成表层土壤的严重流失.土壤侵蚀会去除有机质和重要营养物质，阻碍植被生长，从而对整个生物多样性产生负面影响[2-3].这种特殊现象被认为是对土壤肥力和生产力的最大威胁.它改变了土壤的物理、化学和生物特性，导致潜在农业生产力下降，特别是如今世界人口不断增长的情况下会引起人们对粮食安全的担忧.此外，土壤侵蚀的影响超出了肥沃土地的损失，因为它还可能导致溪流和河流的污染和沉积加剧，关闭水道，导致鱼类和其他海洋物种减少.退化的土地也降低了蓄水能力，有时导致洪水加剧.

目前，国内外学者在土壤侵蚀测量方面采用了许多不同的研究方法.最初，用于估计土壤可蚀性风险的最常见模型是通用土壤流失方程(USLE)，该方程后来被修订为修订的通用土壤流失方程(RUSLE)[4].综合或单独使用许多创新技术，如卫星遥感[5]、地质统计学[6]、野外光谱学、机器学习[7]以及现场和实验室土壤反射率联合测量等等.在试图调查土壤侵蚀状况时，需要收集大量空间分布良好的样本进行进一步分析，这通常是一个昂贵和耗时的过程.因此，将统计学和遥感相结合是一个有效的土壤侵蚀检测方法.特别是地球观测(EO)卫星技术的日趋成熟，为地质监测提供丰富的数据基础.因此，有学者已经综合利用卫星遥感数据、地理信息系统方法和RUSLE方法监测和估计土壤侵蚀率[8-9].然而，上述方法采用的USLE和RUSLE方程都是经验模型，其参数包含不确定性.

1 土壤流失估计

RUSLE代表了气候、土壤、地形和土地利用如何影响由雨滴冲击和地表径流引起的细沟和沟间土壤侵蚀.它已被广泛应用于估算土壤侵蚀损失和评估土壤侵蚀风险.此外，RUSLE可以指出更好的开发和保护计划，以控制不同土地覆盖条件下的侵蚀，如农田、牧场和受干扰的林地.RUSLE方程是：

A=R×K×LS×C×P.

(1)

式中，A为年平均土壤流失量，单位为t ha-1ya-1年；R为降雨侵蚀因子，单位为MJ mm ha-1h-1ya-1；K为土壤可蚀性因子，单位为t h MJ-1mm-1；LS为地形因子(无量纲)，C为覆盖管理因子(无量纲)，P为保护措施因子(无量纲).

1.1 降雨因子R

在USLE和RUSLE模型中，降雨侵蚀力R因子是降雨侵蚀力的一个指标，它是降雨引起侵蚀的潜在能力.R因子计算公式是：

(2)

其中，R为降雨侵蚀指标，Pi为月降雨量.

1.2 土壤因子K

土壤侵蚀K值计算公式是：

K=10-3×(160.8-2.31x1+0.38x2+2.26x3+1.31x4+14.67x5).

(3)

其中，K为土壤侵蚀指标，x1为细砾(1～3 mm)%，x2为细砂(0.05～0.25 mm)%，x3为粗粉砾(0.01～0.05 mm)%，x4为细粉砾(0.005～0.01 mm)%，x5为有机质(10 g/kg).

1.3 地形因子LS

本文采用数字高程模型(DEM)计算LS因子.LS因子在每个DEM像素中计算公式如下，

(4)

式中，Lij-in是网格单元(i,j)的斜率长度；Aij-in是坐标(i,j)的网格单元面积，单位为m2；D是网格单元尺寸，单位为m；m是通用土壤流失方程(USLE)L因子的长度指数；xij为(sinαij+cosαij).

此外，在RUSLE中，m根据细沟和层间侵蚀的比率β而变化，m和β计算公式如下:

(5)

式中，β随坡度变化而变化.

β值由以下公式得出：

(6)

此外，坡度由以下公式计算：

(7)

其中θ是坡度，单位为度.

1.4 地表覆盖管理因子C

在USLE/RUSLE模型中，覆盖系数(C)是一个反映种植方式对土壤侵蚀率影响的指数，它以土地利用为基础.本文以NDVI、土壤调节植被指数(SAVI)和改良土壤调节植被指数(MSAVI)(高分遥感影像不同波段组合确定[10]),3种植被指数作为神经网络的输入层，C因子为输出层.

1.5 土壤保持措施因子P

遥感影像基础上P值需要根据土壤的类型和不同流域管理治理等资料来确定.一般情况下无保持措施地区取1.0，梯田地区取0.15.

2 基于CNN的C因子评估过程

2.1 数据集

数据集为2016年至2021年5景高分辨率地理参考正射影像(分辨率为0.25～0.5 m).所有图像的特性(例如空间分辨率、颜色分布和光照条件)略有不同，但却是在7月底至9月初的植被生长季节进行记录.

此外，根据国家测绘局1980年代测绘的1∶5万地形图，建立了空间分辨率为30 m的数字高程模型(digital elevation model，DEM)，利用研究区的比例尺对土地覆盖等级进行了评价.

2.2 训练过程

用于训练CNN模型的数据包括NDVI、SAVI、MSAVI数据和训练标签.训练标签是将土地利用图和遥感影像相叠加，提取各个土地利用类型内植被指数值，并根据不同土地利用类型的C值作为标签(表1所示).

表1 不同土地利用类型的C值

CNN网络结构如图1所示.首先，网络由卷积层1和ReLU激活构成，后面是最大层化层.对于每个最大层化层，生成的特征映射的大小将减半.在后续卷积层中，应用一系列具有ReLU激活的转置卷积层，以恢复原始图像大小.最后，一个1×1卷积层和一个像素级的softmax激活函数提供最终的C因子输出.softmax函数将每个数据的激活重新调整为[0，1]间隔.

图1 CNN网络结构

文中利用均方误差损失对CNN进行训练，则预测值f(x)与样本真实值c的损失函数计算如下，

(8)

其中n为样本的个数.ci和fxi分别为第i个样本的真实值及其对应的预测值.

3 案例分析

选取密云水库流域数字高程图及遥感图像，并对其进行土壤侵蚀分析.图2所示为该区域覆盖范围示意图.接下来根据第1节相关知识对用降雨因子R、地形因子LS、地表覆盖管理因子C、土壤因子K、土壤保持措施因子P进行计算.

图2 密云水库流域区域覆盖范围示意图

首先，通过GIS对整个流域的R因子输入图进行样条插值.图3所示为案例区域R值的空间分布.可以看出，R值随降水特征由西北向东南逐渐增大.该区域的最小和最大R值分别为0.351和 206.214 MJ mm ha-1h-1ya-1.其次，根据所研究区域地貌特征及公式(3)，K值范围为0.117至 0.397 5 t h MJ-1mm-1.图4所示为案例区域K值的空间分布.

图3 R因子空间分布示意图

图4 K因子空间分布示意图

接着利用第2节方法对C因子进行评估.根据前述可知，数据集为2016年至2021年5景高分辨率地理参考正射影像，将其转化为NDVI、SAVI、MSAVI数据和训练标签，共包含 1 292 个训练样本.仿真环境使用TensorFlow完成模型搭建，训练过程参数如表2所示.图5所示为使用CNN网络计算的C因子空间分布图，值在0.004 1～0.108 9之间，且受人为干扰影响，在低洼处较高.此外，由于包含区域有限，假定该区域的P因子值为1.

表2 网络训练部分参数

图5 C因子空间分布示意图

图6所示为30个实测采样点与CNN评估C值之间的散点图.可以看出C因子与CNN反演的相关系数r为0.929，标准差为0.048.虽然有个别采样点与评估结果有偏差，但大部分采样点能够正确反映C因子结果.

图6 利用CNN计算C因子

接着采用空间分辨率为30 m的地形图和式(4)、(5)绘制了图7所示的坡度小于或等于21%的坡度因子(LS)空间分布图.可以看出，LS值与地表起伏有直接关系.该区域储层的最高LS值计算为20.63.

图7 LS因子空间分布示意图 图8 土壤侵蚀空间分布图

在完成数据输入程序并编制R、K、C、P和LS图作为数据层后，在GIS环境下对其进行倍增，绘制出显示研究区土壤侵蚀空间分布图，以确定每个区域的侵蚀风险，如图8所示.可以看出在研究区内，黄土丘陵区由于土壤可蚀性而具有较大的侵蚀风险.密云水库流域一半以上属于极低侵蚀风险等级，土壤流失量小于10 t ha-1ya-1.土壤流失量从流域东南向西北增加，西北部土壤流失量最大，超过 150 t ha-1ya-1.

4 结语

本文对基于遥感图像的土壤侵蚀评估进行了研究，在此基础上，创新性的提出了利用CNN模型求解RUSLE中C因子.通过仿真分析，结果表明本文所提方法能够利用遥感图像进行土壤侵蚀风险评估.