水土流失方程RUSLE

目录

1 计算降雨侵蚀力 R

2 土壤可蚀性因子 K

3 坡度因子 S 与坡长因子 L

4 地表植被覆盖管理因子 C

5 水土保持措施因子 P

6 计算RUSLE

---

RUSLE方程如下 

A = R * K * LS * C * P            单位

A   年土壤侵蚀模数                  t/（hm2·a）
R   降雨侵蚀力                     （MJ·mm）/（hm2·h·a）
K   土壤可蚀性因子                 （t·h）/（MJ·mm)
LS  坡度因子（S）与坡长因子（L）     无量纲
C   地表植被覆盖管理因子             无量纲
P   水土保持措施因子                无量纲

1 计算降雨侵蚀力 R

采用 Wischmemier 年降雨侵蚀力经验公式

1.1 将 ERA5 2021年的降雨量导出成为12个栅格

将 era5 月尺度降雨数据导出为12个栅格

1.1.1 首先看如何输出一个月的

工具栏中的多维工具中选择创建NetCDF栅格图层。

 导入文件后其它会自己出来，这里需要设置两个地方，如果想输出多维度的栅格图层，那么在波段维度里选择你要叠加的属性即可，例如：可将高度设置为波段维度，从而使创建的多波段栅格中的每个波段都表示该高度上的温度。

在像图里这样，将 波段维度 和 维度值 都留白时，输出的就是按时间叠加的单波段栅格，两个一起选会报错弹出。

当你在维度值里选择了 time ，你就可以如图选择要输出的时间点，输出后就是显示当前时间的栅格，但它仍旧是一个按时间叠加的单波段栅格

 我们需要将它导出，你可以选中，导出，也可以选择 工具栏——数据管理工具 中的复制栅格

 直接输出就可以，输出的就是按时间叠加的单波段栅格当时显示的样子

1.1.2 在模型构建器中批量操作

打开“模型”，先把需要的东西都从右边拖进去，在迭代器中添加 for 循环

 点击模型中的 “创建 NetCDF 栅格图层”，在如下图“值的选择方法”中选择BY INDEX，可以按位置选择所要输出的%值%，这样就不要求For循环的值和NetCDF文件中的值一样了

点击 for ，这些设置很好理解，如果要循环1到12月，应设置 0到11 （python运算规则）但我们现在只能输出 2到12月，原因后边说

拉出一条线指到 “创建 NetCDF 栅格图层”，在此处进行循环。拉线后会弹出选择框，选择前提条件

 从tmp_Layer1拉线到Copy Raster，选择Input Raster，双击Copy Raster，这里需要注意的是，我们要改变Output Raster Dataset 名字，如下图所示，在当中加上“%Value%”（与之前对应），文件后缀名为tif

 完成，运行就可以了，以上步骤参考：Sugarlolly的批量导出各月栅格文件-CSDN博客批量导出各月栅格文件

​

但是值得一提的是，虽然按照值的位置输出，如上图，第一个值的位置应为0，也就是从0开始迭代，但是在arcgis官方文档中说，从0开始的话，后边就没法接别的工具（我理解是这个意思）第一个只能手动导出，所以还是考虑如何用 python 吧

官方文档：在迭代中使用反馈—帮助 | ArcGIS for Desktop

​

1.2 降雨侵蚀的栅格计算

由于 ERA5 的月降水数据实际上是每日的 m (per day)，所以需要将每个数据乘以当月的日数      整个月的降雨总量（mm）：total  (mm per month) =1000* 这个月的月均降雨量*这个月的总天数

太感谢这位了E.die的博客

据里边的评论所说，蒸散发也是的

打开 Spatial Analyst Tools—地图代数—栅格计算器工具，使用下边的公式计算了黄河流域年降水量，得到的平均值为691mm，和实测值647mm较为吻合

​

​

 然后使用下边的公式逐个计算其月降水侵蚀力，最后再加起来

1.735 * Power(10,(1.51 * Log10(Power(("2021monthly01.tif" * 1000*31),2) / "2021monthly_total.tif") - 0.8188))

按照上边的公式输入栅格计算器一个个计算，公式要在word里写，不然会因为换行符问题进行奇怪的报错

​

 最后用这个公式得到国际单位的 R 

2 土壤可蚀性因子 K​

 土壤可蚀因子直接采用 “USLE_K20” 数据集

​

先把范围裁剪出来，栅格需要使用 “按掩膜提取” 裁剪

Spatial analysis——提取分析——按掩膜提取

Spatial analysis——数学分析——转为浮点型

Spatial Analyst Tools—地图代数—栅格计算器工具

​

 做出来值看着很不错，但是极值太大，参考区间在 0.0172 – 0.0869 使用

    （Spatial analysis——提取分析——按属性提取）

输入 VALUE > 0.087，提取大于0.087的部分，结果如下

 可以看到河道附近值明显偏大，我们采用

Spatial Analyst Tools—地图代数—栅格计算器工具

输入Con("raster" > 0.087, 0.087, "raster")，把栅格中大于0.087的都赋值为 0.087

3 坡度因子 S 与坡长因子 L

LS因子直接采用数据集 “LS20”

根据说明文档，根据上文方法做处理

说明文档：

4 地表植被覆盖管理因子 C

4.1 PIE 得到的 NDVI 栅格处理

基于 MOD09A1 的 2021 年黄河流域 NDVI，但是没有进行归一化处理

参考：Arcgis | 栅格数据归一化 - 知乎 (zhihu.com)

可以用栅格计算器归一化一下也可以不归一化，直接放到 envi 里计算植被覆盖度

(b1 lt 0.115)*0+(b1 gt 0.633)*1+(b1 ge 0.115 and b1 le 0.633)*((b1-(0.115))/(0.633-(0.115)))

输出后放到arcgis裁剪一下，就可以了，注意不要用arcgis处理envi默认路径的东西，处理不了

4.2 计算C

根据 蔡崇法等 的计算方法

 c为植被覆盖度，C为地表植被覆盖管理因子，使用

    Spatial Analyst Tools—地图代数—栅格计算器工具

我们注意到公式里的植被覆盖度是 0-100，而我们的是 0-1 ，先扩大100倍

 然后输入下列公式，计算出结果

Con ("veg_basin2_7" > 78.3, 0, (Con ("veg_basin2_7" > 0.1, (0.6508 - 0.3436 * (Log10("veg_basin2_7"))), 1)))

 公式参考：ArcGIS：栅格计算器的运算符和函数详解 (ppmy.cn)

 注意，其实是无法实现链式比较的，只能从左往右挨个判断。

5 水土保持措施因子 P

先裁剪了一个2020LUCC的轮廓，

然后把 slope 数据投影成 CGS WGS 1984，再重采样成和 LUCC 一样的栅格，结果如下

   （ 数据管理工具—栅格—栅格处理—重采样 ）

 使用重分类工具，按照下列要求进行重分类，先给 LUCC赋值，旱地先随便赋值一个，由于本地区没有 “苔原”，所以表里没有，先将值扩大一百倍赋值（第二列），然后再转浮点缩小100

无值 0	0	0
耕地 10	500    后期要用	0.45-0.85
林地 20	100	1
草地 30	100	1
灌木地 40	100	1
湿地 50	10	0.1
水体 60	0	0
苔原 70	——	——
人造地表 80	0	0
裸地 90	500    后期要用	0.45-0.85
冰川和永久积雪 100	0	0
海水 255	0	0

参照上述表格赋值如上，输出得到重分类的结果

Con (("Slope_Res" <= 15) & ("Reclass_Extr" == 500), 45, 0)

Con (("Slope_Res" <= 25) & ("Reclass_Extr" == 500), 65, 0)

Con (("Slope_Res" <100) & ("Reclass_Extr" == 500), 85, "Reclass_Extr")

把这些像上边一样罗列起来就可以了，输出 basin_P100.tif

然后把 basin_P100.tif 转为浮点型，再除以100

6 计算RUSLE

A = R * K * LS * C * P 

由于RULSE公式计算出的土壤侵蚀模数A的单位为 t / (hm2·a)

代入下列公式将单位转换为 t / (km2·a)

Af = 100A       

全部重采样到 LS 的分辨率

依据水利部颁布的《土壤侵蚀分类分级标准》（SL190-2007）土壤侵蚀强度分为

0 – 1000 t/(km^2 ·a)	微度 100 - 80
1000 – 2500 t/(km^2 ·a)	轻度 80 - 60
2500 – 5000 t/(km^2 ·a)	中度 60 - 40
5000 – 8000 t/(km^2 ·a)	强烈 40 - 20
8000 – 15,000 t/(km^2 ·a)	极强烈   20 - 0
> 15,000 t/(km^2 ·a)	剧 烈    0

最后再根据上边的表得分

Con ( "Raster" < 1000, (((1000 - "Raster")/1000)*20+80), 0)

Con ( "Raster" < 2500, (((2500 - "Raster")/1500)*20+60), 0)

Con ( "Raster" < 5000, (((5000 - "Raster")/2500)*20+40), 0)

Con ( "Raster" < 8000, (((8000 - "Raster")/3000)*20+20), 0)

Con ( "Raster" < 15000, (((15000 - "Raster")/7000)*20), 0)

最后分区统计计算出得分即可

收工！看完点个赞吧。

---

2023.4.22 

由于有些地区一个月都没有降水，在采用公式计算时就会因为对数函数问题导致月侵蚀度为空值

如果有些 nc 文件用arcgis打不开，可以将后缀名改为 hdf 再尝试