鲁中山区CSLE模型中植被和地形因子的优化方法

邢先双, 董明明, 郭 静, 孟 琳, 张 玉, 赵登良, 庞海威, 边 振

(1.山东省水文中心, 山东 济南 250002; 2.潍坊市水文中心, 山东 潍坊 266071; 3.济南大学 水利与环境学院, 山东 济南 250022)

土壤侵蚀对生态系统的影响日益引起人们的重视。长期以来，对土壤侵蚀产生泥沙及泥沙运移的定量研究一直是相关学者的研究重点[1]。遥感技术作为开展大范围、长时间序列土壤侵蚀研究的重要技术手段，在中国土壤侵蚀研究中发挥了巨大作用[2]。

目前，中国水土流失遥感监测方法主要有综合评判法和中国土壤流失方程(CSLE， Chinese Soil Loss Equation)模型法[3]。CSLE充分考虑了中国的土壤流失特征，对模型中的各个变量都进行了适合中国自然地理状况的系统性研究[4]，因此，相比于美国通用土壤侵蚀模型(USLE, universal soil loss equation)[5]和修正后的美国通用土壤侵蚀模型(RUSLE, revised universal soil loss equation)[6]，CSLE更适用于中国的水土流失状况的评估。近些年，随着遥感技术、地理信息系统技术、全球定位系统技术的日渐成熟，ArcGIS 10.7平台被广泛应用于土壤侵蚀因子的监测，将土壤侵蚀模型与ArcGIS 10.7平台结合的研究也越来越多。

利用多源遥感数据的优势，有助于解决长时间序列内高分辨率植被覆盖信息获取、提高土壤侵蚀因子监测的准确性[7]，解决了传统水土流失动态监测主要采用中低分辨率遥感数据，其精度难以满足当前水土流失监测与治理需要的问题。本文在前人研究的基础上，将无人机低空遥感数据融入模型计算中，对CSLE模型中的植被覆盖与生物措施因子和地形地貌因子进行参数优化及精度检验，并使用优化后的模型与国家监测的结果进行拟合分析，验证优化后模型的可行性，优化模型方法对同类或相似地区开展水土流失动态监测和水土保持工作有重要的意义。

1 研究区概况

郝峪小流域位于山东省淄博市博山区池上镇西池村(118°05′29″—118°08′06″E，36°19′52″—36°21′23″N)，属小清河支流淄河上游源头区，总面积14.4 km2，在全国水土保持区划中属于北方土石山区(北方山地丘陵区)—泰沂及胶东山地丘陵区—鲁中南低山丘陵土壤保持区。研究区海拔高度范围在350～657 m之间，相对高差307 m，地势东南高西北低，整体向西北倾斜。

研究区属暖温带半湿润大陆性季风气候，多年平均降水量715.1 mm，最大降水量1 332.7 mm(1964年)，最小降水量331.7 mm(1989年)，年平均蒸发量为1 201.4 mm，多年平均气温12.7 ℃，最高气温42 ℃，最低气温-21.8 ℃；年平均日照时数2 607 h；主导风向以南、西风为主，平均风速3.2 m/s。研究区上游顶部为片麻岩及黑云二长混合花岗岩，左半流域为各种阴影混合岩，右半流域为黑云斜长石片麻岩，流域内无断层、裂隙漏水现象。由此母质发育形成的土壤以褐土和棕壤为主。研究区植被覆盖较好，乔木有刺槐、杨树、侧柏、花椒、板栗、桃树、苹果树等，灌木黄荆、酸枣等，草本有黄草、白背草、刺猬皮等，植被覆盖率达90%。

2 研究方法

2.1 数据的获取与处理

(1) 遥感数据(2019年)。主要包括用于解译土地利用的GF-1卫星影像(分辨率2 m)；用于计算归一化植被指数(NDVI)的Landsat系列影像(分辨率30 m)；用于对土地利用解译核查、植被覆盖度精确计算的研究区内无人机航拍数据(优于0.1 m)。遥感影像坐标及投影标准：CGCS 2000国家大地坐标系，1985国家高程基准，投影方式为正轴等面积割圆锥投影(Albers投影)。

(2) 土地利用解译与核查。土地利用解译采用室内人工目视解译与野外调查相结合的方式进行。基于所获取的高分辨率(优于2 m)遥感数据，结合ArcGIS 10.7平台对研究区土地利用类型进行目视解译，结合航拍数据进行修改、核查，确保解译精度。解译研究区2019年土地利用类型总图斑数为119个，野外调查图斑数为15个，占总图斑数的12.61%。土地利用图斑的解译属性及边界吻合正确率为95.97%。

(3) 水土流失数据。淄博郝峪小流域控制站2019年产流产沙资料与国家水土流失动态监测成果中研究区2019年土壤侵蚀模数与各因子栅格数据。

(4) DEM数据。在地理空间数据云网站(http:∥www.gscloud.cn/)获取研究区30 m分辨率DEM数据，进行检查、裁剪、重采样等处理，作为坡度、坡长数据获取以及水土保持措施解译和土壤侵蚀综合分析的基础。在流域内选取两处调查样区(图1)，利用搭载多光谱相机的无人机获取了样区的DSM数据。

图1 郝峪小流域无人机拍摄位置示意图

2.2 土壤侵蚀模型CSLE计算方法

在《2019年度水土流失动态监测技术指南》的指导下，利用已获取的研究区基础数据，基于ArcGIS 10.7平台，引入中国土壤侵蚀模型CSLE对研究区土壤侵蚀进行测算，计算模型为：

A=R·K·L·S·B·E·T

(1)

式中：A表示土壤侵蚀模数〔t/(hm2·a)〕;R表示降雨侵蚀力因子〔(MJ·mm)/(hm2·h·a〕;K表示土壤可蚀性因子〔(t·h)/(MJ·mm)〕;LS表示坡长坡度因子,无量纲;B表示植被覆盖于生物措施因子,无量纲;E表示水土保持工程措施因子,无量纲;T表示耕作措施因子,无量纲。

2.2.1 地形地貌因子LS地形地貌因子包括坡长因子和坡度因子，表征的是地形对坡面土壤侵蚀速率的控制程度。采用符素华等[8]提出的公式计算坡度、坡长因子，计算公式分别为：

坡长因子计算公式：

(2)

式中：λ为坡长(m);m为坡长指数,随坡度而变。

(3)

坡度因子计算公式：

(4)

式中：S为坡度因子(无量纲);θ为坡度(°)。

2.2.2 植被覆盖与生物措施因子B植被盖率指植物群落总体或各个体的地上部分的垂直投影面积与样方面积之比的百分数。一切形式的植被覆盖，均可不同程度地抑制水土流失[9]。基于Landsat系列中分辨率多光谱影像(包括蓝、绿、红和近红外4个波段)及无人机低空遥感数据，根据土地利用类型，结合24个半月降雨侵蚀力因子比例计算B因子。园地、林地和草地采用公式计算，其余土地利用类型直接查表1进行赋值。

(1) 园地、林地B因子计算公式：

(5)

(6)

式中：WRi为第i个半月降雨侵蚀力占全年侵蚀力比例,取值范围为0～1; SLRi为第i个半月园地、林地和草地的土壤流失比例,无量纲,取值范围为0～1。

(2) 灌木林地SLRi计算公式：

(7)

(8)

(9)

式中：NDVI为归一化植被指数; NIR为近红外波段的反射率;R为可见光红波波段的反射率; FVC为基于NDVI计算的植被覆盖度,取值范围为0～1。

(3) 果园、有林地和其他林地SLRi计算公式：

SLRi=0.444 68×e-3.200 96×GD-

0.040 99×eFVC(1-GD)+0.025

(10)

2.2.3 其他因子(R，K，E，T因子) 为研究L，S，B因子优化后对模型计算结果的影响，本研究的R，K，E，T因子与国家水土流失动态监测成果中的栅格数据保持一致，利用研究区边界裁剪后代入公式进行计算。

2.3 模型评价方法

为使CSLE模型更适合鲁中山区实际情况，本研究从提高原始数据空间分辨率角度，利用无人机正射影像与倾斜摄影技术，对模型中的地形地貌因子和植被覆盖与生物措施因子进行了参数优化，并进行优化精度验证。采用的因子优化及验证方法为： ①拟合优度统计分析。基于SPSS分析平台，使用相关性分析工具，确定两组数据拟合曲线，并得出拟合优度的统计量：确定系数R2。R2衡量的是回归方程整体的拟合程度，表达的是变量之间的总体关系其取值范围为：0～1，越接近1说明模型的效率越高[10]。本研究使用该方法对CSLE模型中的植被覆盖度进行了拟合分析，优化了模型中的植被覆盖与生物措施因子。 ②直方图匹配规则。在本研究所选取数据的时间内，流域的地形地貌可以认为是不变的，因此可以使用直方图匹配方法对CSLE模型中参数地形地貌因子(LS)优化结果进行评价。直方图指的是LS因子分布中每一因子值与其出现频数的统计关系，用横坐标表示研究区的LS因子值，纵坐标表示频数。

对于像地形地貌因子这样客观不变的参数，其因子值的统计分布差异是由于研究人员所选取的研究方法导致。相关研究表明，基于不同分辨率的地形数据所提取的LS因子分布差异较大，地形数据的空间分辨率越高，LS因子值越精确[11]。在此研究基础上，通过无人机航摄获取的高程数据提取的高分辨率地形地貌因子LShigh可对基于30 m分辨率DEM数据提取的LSmedium进行修正。

表1 非园地、林地、草地的B因子赋值

基于ArcGIS 10.7平台，提取LShigh和LSmedium的因子值累计频率，根据直方图匹配规则，即从累计频率曲线中指定任意因子值LSh，查看其对应的LShigh的累计频率，然后查看LSmedium中对应频率值的因子值LSm，当无法得到对应的频率值时，使用差分法计算得到LSm。用LSh替换修正LSm，即可得到修正后的，LSmedium修正，并得到参数修正的拟合曲线。通过对比分析尺度转换前后的LS因子值频率曲线的直方图相似度(HS)，间接验证尺度转换模型的精度。用两组数据的累计频率直方图面积差相对指数表示直方图相似度，其计算公式为：

(11)

式中：Shigh表示LShigh的累计频率分布曲线面积；Smedium表示LSmedium的累计频率分布曲线的面积,可以通过微积分算得。HS值域为[0,1],HS值越接近1,两组数据的相似度越高。

3 结果与分析

3.1 土壤侵蚀模型因子的计算与优化结果

3.1.1 地形地貌因子 为提高山东省水土保持动态监测项目中的模型计算精度，本研究以具有完整水文资料的鲁中山区淄博市郝峪小流域为研究区，利用ASTER GDEM数据(空间分辨率30 m)，结合ArcGIS 10.7平台中的水文分析工具提取并计算得出基于中分辨率数据的研究区地形地貌因子并重采样为10 m分辨率，记为LSmedium。利用无人机对研究区中两个有代表性地形的区域进行倾斜摄影，提取DSM数据(空间分辨率优于0.1 m)(图2)。通过点云数据滤波处理，获取地面点后编辑得到的DSM消除了树木和建筑对地形的影响，可作为高分辨率的DEM数据。相同方法计算LS因子并重采样为10 m分辨率，得出高分辨率数据的地形地貌因子LShigh。

在无人机倾斜摄影区域选取100个像元，利用ArcGIS 10.7工具将整个研究区及选取的100个像元的坡度进行分级，结合耕地和园林草的坡度分级将坡度分为8组(表2)。从表2可以看出，研究区和所选的100个像元坡度都以 8°～15°和15°～25°为主，坡度在0°～2°和>35°时最少。同样，将整个研究区和选取的100个像元的LS因子分成8组(表3)。从表3可以看出，研究区和所选的100个像元中LS因子值都以10～20，20～30为主，LS值在70～80时最少。由此可以说明所选取的100个像元值具有代表性，可进行尺度转换。

表2 研究区及所选数据坡度的像元数目分布

图2 地理空间数据云提取的DEM以及无人机遥感影像提取的DSM数据

表3 研究区及所选数据LS因子的像元数目分布

将提取LSmedium和LShigh值进行相关曲线绘制，建立中分辨率与高分辨率地形地貌因子的尺度转换公式，使用该公式对整个项目区LS因子进行修正，得到修正的研究区LSmedium修正。

(12)

利用直方图匹配规则，设置修正步长选取为0.1，绘制LSmedium和LSmedium修正曲线(图3)。

尺度转换前后，通过LS的分布频率和累计分布频率图(图3)分析表明，LSmedium经过尺度转换后生成的LSmedium修正频率曲线与LShigh比较接近，这表明两组数据的统计分布非常相似，研究区在LS因子值小于30的区域，LS因子值呈现错位分布，研究区大于30的区域，二者分布几乎重合，这说明LS因子尺度转换模型有一定效果。

图3 研究区LS因子尺度转换前后频率及累计频率分布

经计算，尺度转换后的LS因子值频率曲线相似度为0.91，这表明两组因子值的相似度较高，尺度转换后的LS因子可以用来计算研究区的土壤侵蚀。选取无人机倾斜摄影范围内不同高程梯度共30个像元的LSmedium修正值与无人机倾斜摄影提取的DSM数据进行拟合分析，线性相关曲线R2为0.697 1。进一步运用统计分布和直方图相似度进行精度验证。统计分布是指统计分析LS因子转换前后的平均值、标准差、最大值和最小值以及频率和累计频率曲线图。LShigh和LSmedium以及LSmedium修正的频率曲线图(图3)表明，LSmedium频率曲线峰值较小，且峰值附近频率高，数据较集中于12～25之间。LShigh和LSmedium修正频率曲线特征相似，峰值较LSmedium有所增大，且数值分布更加均匀。统计结果(表4)表明，修正后的LSmedium修正平均值、标准差相较LSmedium与LShigh更加接近，LSmedium修正整体较中分辨率的LSmedium增大4.6%，标准差较LSmedium增大8.5%。以上频率特征的差异主要是由于中分辨率DEM数据由于空间分辨率较低，局部高程差异被消除，高程变化趋于平滑，高程变幅也有所减小，导致坡度、坡长因子值和空间变幅整体小于实际数据。而无人机DSM数据由于空间分辨率有较大提高，能较精确的反映细微的地形差异，因此基于无人机的DSM计算的LS因子LShigh和修正后的LSmedium修正能够更准确地反映实际地形地貌信息，其精度较LSmedium有所提高。

表4 尺度转换前后LS因子值统计结果

3.1.2 植被覆盖因子与生物措施因子 基于不同分辨率的遥感数据提取的植被覆盖度会有较大的差异，相较于Landsat系列卫星影像，搭载多光谱相机的无人机平台可以提取更为精确的植被覆盖数据。但是后者的耗时较大，成本较高，所以本研究于植物充分生长且云量相对较少的9月，在研究区选取两个样地进行无人机影像采集，经处理计算后生成分辨率为0.1 m的无人机正射影像。

基于Landsat影像提取的24期植被覆盖度栅格数据中有两期与无人机正射影像的时相接近。使用栅格计算器对两期植被覆盖度栅格数据求均值并重采样后得到分辨率为10 m的中分辨率覆盖度栅格数据FVCmedium。使用Creat Fishnet工具创建与FVCmedium像元边界完全重合的解译网格(大小为10 m×10 m)。使用Zonal工具分区统计各10 m×10 m网格内植被覆盖度，得到10 m分辨率高精度植被覆盖度栅格数据 FVCUVA。基于采样空间分布均匀和各土地利用类型栅格数占栅格总数比例一致的原则，选取1 000个采样点，使用ArcGIS 10.7中Sample工具提取FVCUVA及对应的FVCmedium像元值，共采集1 000个数值对。各土地利用类型采样点选取情况详见表5。

表5 各土地利用类型栅格数据采样比例

对两组栅格数据进行拟合分析，不同函数关系下的相关系数(表6)表明，FVCUVA和 FVCmedium呈现较好的相关性，各类函数的R2范围在0.565 2～0.686 8之间，选择R2最大的幂函数作为植被覆盖度的尺度转换模型，其表达式为：

(13)

式中：FVChigh为FVCUVA与FVCmedium最优拟合方程计算后的优化后植被覆盖度栅格数据,FVCmedium为同时相中分辨率影像计算的植被覆盖度。

FVChigh和FVCmedium的幂函数关系分布散点图如图4所示。

表6 FVChigh和FVCmedium拟合分析

图4 FVChigh和FVCmedium栅格数值分布

利用上文公式对基于中分辨率的FVC进行修正，得到研究区高分辨率植被覆盖度FVChigh。在无人机拍摄区域有林地和果园根据覆盖度梯度抽取100个像元进行修正后植被覆盖度与正射影像覆盖度FVCUVA拟合分析(图5)，R2为0.718 7，表明优化后的植被覆盖度FVChigh更接近实际植被覆盖度。

3.2 基于多源遥感的水土流失监测结果与分析

利用ArcGIS 10.7平台将郝峪小流域的土壤侵蚀数据从山东省土壤侵蚀数据提取出，再通过优化后的模型计算出郝峪小流域的土壤侵蚀强度，结合《2019年度水土流失动态监测技术指南》中对土壤侵蚀量的分级标准，利用重分类工具进行侵蚀强度分级(图6)。

图5 郝峪小流域FVChigh和FVCUVA拟合分析

图6 郝峪小流域原侵蚀强度和优化后侵蚀强度等级分布

利用ArcGIS 10.7平台的分区统计工具对国家监测和优化后的的土壤侵蚀强度进行面积统计，结果详见表7。从表7可以看出，优化LS因子和B因子后的模型计算所得水土流失面积较国家监测的水土流失面积增加0.85 km2，且全部为轻度侵蚀。为更好地验证优化后模型的可行性，随机从国家监测的土壤侵蚀数据和优化模型后的土壤侵蚀数据中提取100组数据进行拟合分析(图7)，R2为0.700 3，说明优化后的土壤侵蚀数据与国家监测的土壤侵蚀数据具有良好的相关性，同时又存在一定差异。由线性拟合方程可知，优化模型因子后，侵蚀模数变化范围有所减小，优化模型因子后的土壤侵蚀模数计算结果，在较小值区域大于国家监测数据，即国家监测成果中侵蚀模数小于947.53 t/(hm2·a)的区域经优化计算后，数值有所增大；反之在侵蚀模数值较大区域，优化计算结果小于国家监测成果。

表7 国家监测数据与优化模型后土壤侵蚀强度面积分析

鲁中山区中度及以上侵蚀主要发生在灌木林地和旱地。原始数据空间分辨率的提高，使植被覆盖度较中分辨率数据计算结果有所增大，减小B因子值；而地形因子计算中能够更好地识别梯田区域从而减小LS因子值。轻度和微度侵蚀主要发生在有林地，为植被覆盖度较高区域，优化后的B因子有所增大，导致侵蚀模数和轻度侵蚀强度面积的增加。同时由于微度和轻度侵蚀面积(合计14.41 km2)占研究区总面积(14.44 km2)的99.79%，因此微度和轻度侵蚀间相互转移的概率最大，且侵蚀模数在某一阈值范围内，对应的侵蚀强度不会发生变化，由实际计算结果统计可知，中度及以上侵蚀区域内的土壤侵蚀模数有所变化，但侵蚀强度未发生变化。

经过现场调查发现，研究区近年来旅游业发展和生态建设效果良好，不同坡面、坡度与土地利用类型区水土流失状况差异不大，与优化后变化特征一致，但仍缺乏定量资料，验证优化后模型的计算精度。

图7 优化模型后的数据与国家监测数据的相关分析

4 结 论

(1) 本研究基于多元遥感数据，通过无人机高分辨率影像与中分辨率数据的拟合分析与精度验证，从空间分辨率的提高方面，提高了DEM和FVC解译精度。利用直方图匹配规则，实现了LS因子的尺度转换，验证结果表明，LS因子频率曲线相似度为0.91，可以进行研究区的土壤侵蚀模型计算。对多源遥感数据植被覆盖度提取的尺度转换进行了探讨，得出R2为0.686 8的幂函数转换模型。

(2) 对CSLE模型中因子进行优化的基础上计算了研究区土壤侵蚀强度，并与该区域国家监测土壤侵蚀数据进行了对比分析。结果表明，优化模型因子后的研究区侵蚀模数空间异质性有所减小，轻度侵蚀面积有所增加，中度及以上侵蚀面积未变化，但侵蚀模数总体减小，这与LS因子和B因子优化后数值变化特征一致，说明模型因子的优化对监测结果具有显著影响，但对侵蚀模数计算精度的提高是否具有显著作用，仍需从土壤侵蚀机理上结合进一步试验研究加以验证。

(3) 不同空间分辨率的遥感数据源对土壤侵蚀的各因子计算均会造成一定的影响，优化后的LS因子和B因子更加接近无人机DSM数据处理的LS值和正射影像下的植被覆盖情况，即以上因子精度均有所提高，继而有助于提高模型计算结果精度，使得基于优化后因子计算的土壤侵蚀模数更符合当地实际情况。但本研究目前仅对地形地貌因子以及植被覆盖与生物措施因子进行了尺度转换和精度提升，对其他因子如何利用高分辨率数据源进行优化，并提高土壤侵蚀模型计算精度仍需进一步探讨。