一种水土保持综合治理土壤侵蚀变化量实时定量监测方法

摘要

本发明公开了一种水土保持综合治理土壤侵蚀变化量实时定量监测方法，包括：获得研究区水土保持综合治理前和治理后遥感图像；提取数字高程图DEM以及各水土保持措施遥感影像图斑；以影像图斑为单元，计算各图斑坡长因子L、坡度因子S；采用归一化植被指数NDVI计算各图斑植被覆盖度B，并根据植被覆盖度B估算该图斑植被覆盖或作物因子C；依据现有研究成果对各图斑水土保持措施因子P进行赋值；计算土壤侵蚀量减少比例SEDP＝1-SE治理后/SE治理前，其中SE治理前和SE治理后分别表示研究区水土保持综合治理前和治理后土壤侵蚀量。本发明仅通过对遥感影像信息提取即可得到土壤侵蚀量变化，投入成本相对较低、方便快捷。

说明书

技术领域

本发明涉及水土保持综合治理研究领域，特别涉及一种水土保持综合治理 土壤侵蚀变化量实时定量监测方法。

背景技术

土壤侵蚀量变化是水土保持综合治理工程成效的最直接反映，也是水土保 持综合治理工程的重要监测因子，对于水土保持综合治理工程效益评价具有重 要的参考作用。

目前，水土保持综合治理工程中，土壤侵蚀量变化仍沿用“径流小区观测 试验”、“小流域对比法”等，这些方法虽然能够真实反映治理前后土壤侵蚀量 变化，但观测指标多，耗时费工，且缺乏治理工程各个阶段的实时定量监测功 能。

因此，寻求一种能够简便的、适用于一般水土保持综合治理工程且能够实 时定量反映土壤侵蚀变化量的方法就显得十分必要。

发明内容

本发明的目的在于克服现有方法的缺点与不足，提供一种水土保持综合治 理土壤侵蚀变化量实时定量监测方法，该方法不需要野外定点观测试验设施， 仅通过对研究区遥感影像信息提取即可得到土壤侵蚀变化量，简便易用且投入 成本相对较低。

本发明目的通过以下技术方案实现：一种水土保持综合治理土壤侵蚀变化 量实时定量监测方法，包括如下步骤：

(1)获得研究区在进行水土保持综合治理前的遥感影像和治理后的遥感影 像；

(2)对遥感影像进行校正和信息提取，提取数字高程图DEM和各水土保 持措施遥感影像图斑；所述各水土保持措施遥感影像图斑是指根据实际情况进 行各种水土保持措施后得到的每种措施对应的遥感影像图斑；

(3)以各水土保持措施遥感影像图斑为计算单元，基于数字高程图DEM， 计算各图斑坡长因子L、坡度因子S；

(4)采用归一化植被指数NDVI计算各图斑植被覆盖度B，并根据植被覆 盖度B推算对应图斑植被覆盖或作物因子C；

(5)对各水土保持措施遥感影像图斑对应的水土保持措施因子P进行赋值；

(6)对水土保持综合治理前的因子L、S、C和P进行相乘运算得到SE_治理前，对水土保持综合治理后的因子L、S、C和P进行相乘运算得到SE_治理后；

(7)按照公式：SEDP＝1-SE_治理后/SE_治理前，得到研究区域分水土保持措施图 斑的土壤侵蚀量减少比例SEDP；

(8)基于各水土保持措施图斑面积，对各水土保持措施图斑对应的SEDP 按照加权平均的方法进行计算，获得整个研究区域的土壤侵蚀量变化，计算步 骤如下：

假设研究区域有n个水土保持措施图斑，各水土保持措施图斑面积为AR_i， 1≤i≤n，各水土保持措施图斑对应的水土保持综合治理前后土壤侵蚀量减少比 例依次为SEDP_i，那么以图斑面积为加权的整个研究区域土壤侵蚀量减少比例即 为：

SEDP_{整个研究区域}＝(SEDP₁×AR₁+SEDP₂×AR₂+…+SEDP_i×AR_i)/

(AR₁+AR₂+…+AR_i)。

优选的，所述监测方法基于GIS软件实现，步骤(3)中在GIS软件中生成 研究区域坡长因子L和坡度因子S专题图，步骤(4)中在GIS软件中生成研究 区域植被覆盖或作物因子C专题图；步骤(5)中在GIS软件中生成研究区域水 土保持措施因子P专题图；步骤(6)中，通过GIS软件空间要素叠加分析对因 子L、S、C和P进行相乘叠加运算分别得到SE_治理前、SE_治理后；步骤(7)中， 对步骤(6)计算得到的SE_治理前、SE_治理后在GIS软件中进行叠加分析，得到研究 区域各水土保持措施图斑的土壤侵蚀量减少比例SEDP。

优选的，所述步骤(2)中，对遥感影像进行校正和信息提取具体步骤是：

(2-1)在ENVI软件下基于之前地面定标点对遥感影像进行影像配准和几 何校正，同时在ENVI软件下对遥感影像进行大气校正；

(2-2)基于校正后的遥感影像在GIS中进行信息提取，提取的信息包括两 个方面：生成影像区域的数字高程图DEM和提取各水土保持措施遥感影像图斑。

优选的，所述步骤(3)中，采用美国通用水土流失方程推荐公式，计算研 究区内各图斑的坡长因子L、坡度因子S，其计算方法采用下述方法中的任意一 种：

(3-1-1)坡长因子L的计算公式如下：

L＝(λ/22.13)^m；

m＝β/(1+β)；

β＝(sinθ/0.0896)/[3.0(sinθ)^0.8+0.56]；

坡度因子S的计算公式如下：

其中，θ为图斑平均坡度值，λ为各图斑平均坡长，不超过300m；

(3-1-2)坡长因子L的计算公式如下：L＝(λ/22.13)^0.44；

坡度因子S的计算公式如下：S＝17/(1+e^{2.3-6.1sin(θ)})-1.5；

其中λ为各图斑平均坡长；θ为图斑平均坡度值。

更进一步的，所述步骤(3)中，将各图斑的坡长因子L、坡度因子S的乘 积作为一个坡长坡度因子，记为LS，其计算方法采用下述方法中的任意一种：

(3-2-1)LS＝(λ/22.13)^m×(21.91sin(θ)-0.96)；

其中LS为坡长坡度因子；λ为各图斑平均坡长，不超过300m；θ为图斑平 均坡度值；m为坡长指数，当θ≤0.5°，m＝0.2；当0.5°<θ≤1.5°，m＝0.3； 当1.5°<θ≤3°，m＝0.4；当θ>3°，m＝0.5；

(3-2-2)LS＝(λ/22.13)^m×(65.4sin(θ)²+4.56sin(θ)+0.065)；

其中LS为坡长坡度因子；λ为各图斑平均坡长；θ为图斑平均坡度值；m 为坡长指数；该公式针对坡度大于2.86°适用，此时m＝0.5；

(3-2-3)LS＝(A/22.13)^m×(sin(θ)/0.0896)ⁿ；

其中LS为坡长坡度因子，A为图斑面积，θ为图斑平均坡度值；m和n为 常数，分别取0.4～0.6和1.2～1.3；

(3-2-4)LS＝(λ/22.0)×(θ/5.16)^1.3；

其中LS为坡长坡度因子，λ为各图斑平均坡长；θ为图斑平均坡度值，该 公式适用于坡度大于5.14°的研究区。

优选的，所述步骤(4)中，计算植被覆盖或作物因子C的方法具体为：

(4-1)在ENVI软件中计算各图斑植被指数NDVI；

(4-2)计算各图斑植被覆盖度B，参照如下公式：

B＝(NDVI_C-NDVI_MIN)/(NDVI_MAX-NDVI_MIN)；

其中，NDVI_C是计算当前图斑NDVI值，NDVI_MAX和NDVI_MIN分别为全幅 遥感影像中NDVI的最大和最小值；

(4-3)计算各图斑植被覆盖或作物因子C，公式如下：

$C=\left(\begin{array}{cc}1,& B=0;\\ 0.6508-0.3436\lg \left(B\right),& 0<B<78.3\%;\\ 0,& B=1;\end{array}\right).$

更进一步的，所述步骤(4)中，植被覆盖或作物因子C通过查表赋值的方 法得到，具体为：

表1植被覆盖或作物因子C查询表1

将所述植被分为林草地和农地两种类型，对于林草地，根据其植被或作物 类型，结合步骤(4-2)计算得到的植被覆盖度B，通过查询表1，得到对应的植被 覆盖或作物因子C；对于农地，根据其植被或作物类型，结合耕作方式，得到对 应的植被覆盖或作物因子C。

优选的，所述步骤(4)中，植被覆盖或作物因子C通过查表赋值的方法得 到，具体为：

表2植被覆盖或作物因子C查询表2

土地类型 植被覆盖度B 植被覆盖或作物因子C 坡地 10 0.99 乔木林地 30 0.83 混交林 30 0.83 经济林 20 0.89 草地 50 0.75 梯田 20 0.89 荒地 80 0.68 居民地 20 0.89

将所述植被根据土地类型分为坡地、乔木林地、混交林、经济林、草地、 梯田、荒地、居民地，通过查询表2，得到对应的植被覆盖度B、植被覆盖或作 物因子C。

优选的，所述步骤(5)中，对各水土保持措施遥感影像图斑对应的水土保 持措施因子P进行赋值，其赋值方法采用下面中的任意一种：

(5-1)将水土保持措施分为以下几类：耕作措施、工程措施、林草措施、 综合措施，其中耕作措施包括横坡耕作、横坡植物篱、垄作、平作、免耕，工 程措施包括水平梯田、截流沟、水平沟，林草措施包括水保林、经果林、林草 复合模式，综合措施包括封禁治理、人工林，根据上述分类，采用表3对各水 土保持措施遥感影像图斑对应的水土保持措施因子P进行赋值：

表3水土保持措施因子P查找表1

(5-2)将水土保持措施分为等高带状种植、跨坡等高条植、等高耕作、条 植、等高梯田、水平沟、等高边界、植物篱、生物覆盖、造林，根据上述分类， 采用表4对其对应的水土保持措施因子P进行赋值：

表4水土保持措施因子P查找表2

水土保持措施 P值 等高带状种植 0.4 跨坡等高条植(非在等高线上) 0.5 等高耕作 0.55 条植(有草或作物带) 0.15 等高梯田 0.1 水平沟 0.1 等高边界 0.2 植物篱 0.35 生物覆盖 0.35 造林 0.4–0.6

(5-3)将水土保持措施根据坡度分为顺坡耕作、等高沟垄种植、等高带状 耕作、水平梯田根据上述分类，采用表5对其对应的水土保持措施因子P进行 赋值：

表5水土保持措施因子P查找表3

坡度/度 顺坡耕作 等高沟垄种植 等高带状耕作 水平梯田 2-5 1.0 0.3 0.5 0.1 5.1-9 1.0 0.4 0.7 0.16 9.1-16 1.0 0.5 0.9 0.3 16.1-20 1.0 0.7 1.0 0.4 20.1-25 1.0 0.9 1.0 0.5

(5-4)将水土保持措施根据用地类型分为坡地、乔木林地、混交林、经济 林、草地、梯田、荒地、居民地，见表6：

表6水土保持措施因子P查找表4

用地类型 坡地 乔木林地 混交林 经济林 草地 梯田 荒地 居民地 P值 1.0 1.0 0.2 0.69 0.2 0.35 1.0 1.0

根据上述分类，采用表6对其对应的水土保持措施因子P进行赋值。

优选的，所述步骤(6)中，基于美国通用水土流失方程，所述SE′_治理前和SE′_治理后的计算公式如下：

SE′_治理前＝R_治理前×K_治理前×L_治理前×S_治理前×C_治理前×P_治理前；

SE′_治理后＝R_治理后×K_治理后×L_治理后×S_治理后×C_治理后×P_治理后；

其中，R为降雨侵蚀力，单位为MJ mm h^-1hm^-2a^-1；K为土壤可蚀性因子， 单位为t hm²h MJ^-1hm^-2mm^-1；L为坡长因子；S为坡度因子；C为植被覆盖或 作物因子；P为水土保持措施因子；

对于同一特定研究区域，治理前后R、K基本相等，因此，土壤侵蚀量减少 比例SEDP的计算公式如下：

SEDP＝1-SE′_治理前/SE′_治理后

＝1-SE_治理后/SE_治理前

＝1-(L_治理后×S_治理后×C_治理后×P_治理后)÷(L_治理前×S_治理前×C_治理前×P_治理前)；

若SEDP＝0，则表示治理前后土壤侵蚀量未发生改变；若SEDP<0，则表示 治理后土壤侵蚀量增加；若SEDP>0，则表示治理后土壤侵蚀量降低。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明不需要野外建设定点观测试验设施，仅通过对遥感影像信息的提 取即可得到土壤侵蚀量变化情况，从而对水土保持综合治理工程土壤流失变化 量进行实时定量监测，投入成本低、方便快捷。

2、本发明对土壤侵蚀变化量实时定量监测是对研究区进行全尺度覆盖分析， 与以往抽样调查方法有明显不同。

3、本发明设立了一个指标土壤侵蚀量减少比例SEDP，计算中忽略降雨侵 蚀力(R)和土壤可蚀性(K)变化，在不计算R、K参数的情况下，仍可以实 时定量获得研究区土壤侵蚀量变化情况，具有计算简单的优点。

附图说明

图1是本实施例的算法流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方 式不限于此。

实施例1

本实施例一种水土保持综合治理土壤侵蚀变化量实时定量监测方法，基于 以下工作原理：

基于美国通用水土流失方程(USLE)SE＝K×R×L×S×C×P，其中SE为土壤 侵蚀量，t hm^-2a^-1；K为土壤可蚀性因子，t hm²h MJ^-1hm^-2mm^-1；R为降雨侵蚀 力，MJ mm h^-1hm^-2a^-1；L为坡长因子；S为坡度因子；C为植被或作物因子；P 为水土保持措施因子。那么，土壤侵蚀量减少比例SEDP即为：

SEDP＝1-SE′_治理前/SE′_治理后

＝1-(R_治理后×K_治理后×L_治理后×S_治理后×C_治理后×P_治理后)÷(R_治理前×K_治理前×L_治理前×

S_治理前×C_治理前×P_治理前)         (1)

针对某一特定研究区域，其降雨侵蚀力(R)与土壤可蚀性(K)与水土保 持工程实施与否并无物理性内在因果关系，可认为降雨侵蚀力(R)与土壤可蚀 性(K)基本未发生改变。即为R_治理后≈R_治理前和K_治理后≈K_治理前，那么公式(1)可 化简为：

SEDP＝＝1-SE′_治理前/SE′_治理后

＝1-SE_治理后/SE_治理前

＝1-(L_治理后×S_治理后×C_治理后×P_治理后)÷(L_治理前×S_治理前×C_治理前×P_治理前)    (2)

式中SE′_治理前和SE′_治理后分别表示水土保持工程实施后和前研究区土壤侵蚀量； L_治理后、S_治理后、C_治理后、P_治理后和L_治理前、S_治理前、C_治理前、P_治理前分别表示水土保持 工程实施后和前研究区域坡长因子、坡度因子、植被或作物因子、水土保持措 施因子。

可以看出，土壤侵蚀量减少比例SEDP存在以下几个可能值：若SEDP＝0， 则表示治理前后土壤侵蚀量未发生改变；若SEDP<0，则表示治理后土壤侵蚀量 增加了；若SEDP>0，则表示治理后土壤侵蚀量降低了。从而实现对水土保持综 合治理工程土壤侵蚀变化量动态定量监测评价。

在实际应用中，还可以针对具体项目研究问题，设定不同等级SEDP值， 进而对水土保持治理工程的土壤侵蚀变化量进行定量评价。

下面结合图1，对本实施例一种水土保持综合治理土壤侵蚀变化量实时定量 监测方法的步骤说明如下：

步骤1：获得研究区在进行水土保持综合治理前和治理后的遥感影像。

步骤2：在ENVI软件下基于之前地面定标点对遥感影像进行影像配准和几 何校正，同时在ENVI软件下对遥感影像进行大气校正；基于校正后的遥感影像 在GIS中进行信息提取，提取的信息包括两个方面：生成影像区域的数字高程 图DEM和提取各水土保持措施遥感影像图斑。本实施例所述水土保持措施分为 以下几类：耕作措施、工程措施、林草措施、综合措施，其中耕作措施包括横 坡耕作、横坡植物篱、垄作、平作、免耕，工程措施包括水平梯田、截流沟、 水平沟，林草措施包括水保林、经果林、林草复合模式，综合措施包括封禁治 理、人工林。上述的各水土保持措施类型并非实际生产中的全部类型。因此， 本方法在实际应用中可根据实际情况进行“水土保持措施”分类。

步骤3：以各图斑为计算单元，基于数字高程图DEM，计算各图斑的坡长 因子L、坡度因子S。

坡长因子L的计算公式如下：

L＝(λ/22.13)^m；

m＝β/(1+β)；

β＝(sinθ/0.0896)/[3.0(sinθ)⁰_.⁸+0.56]；

坡度因子S的计算公式如下：

其中，θ为计算图斑平均坡度值，λ为各图斑平均坡长，不超过300m。

本实施例中计算坡长因子L和坡度因子S的方法来自于已有研究成果中的 一种，在实际使用时也可因地制宜采用其他适用性理想的计算方法。例如：

(1)将坡长与坡度作为一个因子(LS)，计算方法如下：

①LS＝(λ/22.13)^m×(21.91sin(θ)-0.96)

其中LS为坡长坡度因子；λ为各图斑平均坡长，不超过300m；θ为计算图 斑平均坡度值；m为坡长指数，当θ≤0.5°，m＝0.2；当0.5°<θ≤1.5°，m＝0.3； 当1.5°<θ≤3°，m＝0.4；当θ>3°，m＝0.5。

②LS＝(λ/22.13)^m×(65.4sin(θ)²+4.56sin(θ)+0.065)

其中LS为坡长坡度因子；λ为各图斑平均坡长；θ为计算图斑平均坡度值； m为坡长指数；这里公式针对坡度大于2.86°适用，此时m＝0.5。

③LS＝(A/22.13)^m×(sin(θ)/0.0896)ⁿ

其中LS为坡长坡度因子，A为计算图斑面积，θ为计算图斑平均坡度值； m和n为常数，分别取0.4～0.6和1.2～1.3。

④LS＝(λ/22.0)×(θ/5.16)^1.3

其中LS为坡长坡度因子，λ为各图斑平均坡长；θ为计算图斑平均坡度值。 该公式适用于坡度大于5.14°的研究区。

(2)将坡长L与坡度S分别计算，方法如下：

L＝(λ/22.13)^0.44

S＝17/(1+e^{2.3-6.1sin(θ)})-1.5

其中λ为各图斑平均坡长；θ为计算图斑平均坡度值。

步骤4：在ENVI软件中计算各图斑植被指数NDVI，计算各图斑植被覆盖 度B，参照如下公式：

B＝(NDVI_C-NDVI_MIN)/(NDVI_MAX-NDVI_MIN)；

其中，NDVI_C是计算当前图斑NDVI值，NDVI_MAX和NDVI_MIN分别为全幅 遥感影像中NDVI的最大和最小值。

参照蔡崇法公式计算各图斑植被覆盖或作物因子C，公式如下：

$C=\left(\begin{array}{cc}1,& B=0;\\ 0.6508-0.3436\lg \left(B\right),& 0<B<78.3\%;\\ 0,& B=1;\end{array}\right).$

这里C因子计算方法采用的已有研究成果的一种，在实际使用中也可因地 制宜采用其他适用性理想的计算方法。例如：

方法(1)：

方法(2)：

土地类型 植被覆盖度B 植被覆盖或作物因子C 坡地 10 0.99 乔木林地 30 0.83 混交林 30 0.83 经济林 20 0.89 草地 50 0.75 梯田 20 0.89 荒地 80 0.68 居民地 20 0.89

步骤5：对各图斑水土保持措施因子(P)进行赋值。

本实施例参见陈正发、张有全等研究结果，建立一P因子查找表，表如下：

表1 P因子查找表

上述P因子查找表来源于已有研究成果的一种，在实际使用中也可因地制 宜采用其他适用性理想的计算方法。例如：

方法(1)：

等高带状种植 0.4 跨坡等高条植(非在等高线上) 0.5 等高耕作 0.55 条植(有草或作物带) 0.15 等高梯田 0.1 水平沟 0.1 等高边界 0.2 植物篱 0.35

生物覆盖 0.35 造林 0.4–0.6

方法(2)：

坡度/度 顺坡耕作 等高沟垄种植 等高带状耕作 水平梯田 2-5 1.0 0.3 0.5 0.1 5.1-9 1.0 0.4 0.7 0.16 9.1-16 1.0 0.5 0.9 0.3 16.1-20 1.0 0.7 1.0 0.4 20.1-25 1.0 0.9 1.0 0.5

方法(3)：

用地类型 坡地 乔木林地 混交林 经济林 草地 梯田 荒地 居民地 P值 1.0 1.0 0.2 0.69 0.2 0.35 1.0 1.0

步骤6：根据上述步骤计算得到的坡长因子L、坡度因子S、植被覆盖或作 物因子C、水土保持措施因子P，在GIS软件中生成研究区水土保持综合治理工 程实施前后L、S、C和P值专题图。这里的专题图为栅格格式，即为一般的TIFF、 JPEG等格式即可。

步骤7：基于美国通用水土流失方程，所述SE′_治理前和SE′_治理后的计算公式如下：

SE′_治理前＝R_治理前×K_治理前×L_治理前×S_治理前×C_治理前×P_治理前；

SE′_治理后＝R_治理后×K_治理后×L_治理后×S_治理后×C_治理后×P_治理后；

其中，R为降雨侵蚀力，单位为MJ mm h^-1hm^-2a^-1；K为土壤可蚀性因子， 单位为t hm²h MJ^-1hm^-2mm^-1；L为坡长因子；S为坡度因子；C为植被覆盖或 作物因子；P为水土保持措施因子；对于同一特定研究区域，治理前后R、K基 本相等，因此可通过GIS软件空间要素叠加分析，对因子L、S、C和P进行相 乘运算得到水土保持综合治理前和治理后的4因子乘积值，分别记为SE_治理前、 SE_治理后。

步骤8：计算某一水土保持措施图斑对应的土壤侵蚀量减少比例SEDP，公 式如下：

SEDP＝1-SE_治理后/SE_治理前

＝1-(L_治理后×S_治理后×C_治理后×P_治理后)÷(L_治理前×S_治理前×C_治理前×P_治理前)；

若SEDP＝0，则表示治理前后土壤侵蚀量未发生改变；若SEDP<0，则表示 治理后土壤侵蚀量增加；若SEDP>0，则表示治理后土壤侵蚀量降低。

步骤9：通过加权平均等方法获得整个研究区土壤侵蚀变化量及生成相关图 件。计算步骤如下：

假设研究区域有n个水土保持措施图斑，各水土保持措施图斑面积为AR_i， 1≤i≤n，各水土保持措施图斑对应的水土保持综合治理前后土壤侵蚀量减少比 例依次为SEDP_i，那么以图斑面积为加权的整个研究区域土壤侵蚀量减少比例即 为：

SEDP_{整个研究区域}＝(SEDP₁×AR₁+SEDP₂×AR₂+…+SEDP_i×AR_i)/

(AR₁+AR₂+…+AR_i)。

SEDP_{整个研究区域}即为整个研究区域土壤侵蚀量减少比例。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实 施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、 替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。