土壤侵蚀方程中植被覆盖与管理因子C改进方法

摘要

本发明涉及一种土壤侵蚀方程中植被覆盖与管理因子C改进方法，引入遥感数据时空融合模型和地形调节植被指数TAVI，改进植被覆盖与管理因子C的计算方法，即利用地形调节植被指数TAVI代替传统的归一化指数NDVI来估算植被覆盖与管理因子C以消除地形引起的阴坡和阳坡差异，通过逐月植被覆盖与管理因子C与降雨侵蚀力因子R的相乘累加改进传统仅利用单一时相植被覆盖与管理因子C来计算土壤侵蚀模数，有效地匹配了植被覆盖与降雨的年内变化。本发明既消除了地形引起的阴坡和阳坡差异，又顾及了植被的季节变化特征，提高了C因子的估算精度，有效提高了土壤流失强度的估算精度与合理性。

说明书

技术领域

本发明涉及一种土壤侵蚀方程中植被覆盖与管理因子C改进方法。

背景技术

美国通用土壤流失方程及修正方程（USLE/RUSLE）是世界上应用最为广泛的土壤侵蚀预报模型，其中植被与经营管理因子C的值是模型诸因子中变化幅度最大的，可相差2—3个数量级，对土壤侵蚀最敏感，其合理估算对土壤侵蚀的准确预测尤为重要。

利用遥感数据计算的植被覆盖度FVC进行C因子的估算应用较多，目前一般利用归一化植被指数NDVI来计算FVC，没有考虑由于地形引起的阴坡和阳坡植被覆盖度在遥感影像存在差异。同时，植被具有明显的年内变化周期，利用不同月份的遥感数据所获得的地表FVC会有较大的差异。这些差异使得仅利用某一时期的影像，基于常用的植被指数如NDVI，计算获得的C因子的估算存在较大不确定性，进而导致土壤侵蚀结果不可靠。针对这两个问题，引入遥感数据时空融合模型和地形调节植被指数TAVI，改进C因子的计算方法，即利用TAVI代替传统的NDVI来估算FVC进而计算C因子来消除地形引起的阴坡和阳坡差异，通过逐月C因子与降雨侵蚀力R因子的相乘累加改进传统仅利用单一时相C因子来计算土壤侵蚀模数，有效地匹配了植被覆盖与降雨的年内变化。该改进的C因子计算方法，既消除了地形引起的阴坡和阳坡差异，又顾及了植被的季节变化特征，提高了C因子的估算精度，有效提高了土壤流失强度的估算精度与合理性。

发明内容

本发明的目的在于提供一种土壤侵蚀方程中植被覆盖与管理因子C改进方法，以克服现有技术中存在的缺陷。

为实现上述目的，本发明的技术方案是：一种土壤侵蚀方程中植被覆盖与管理因子C改进方法，包括如下步骤：

步骤S1：获取研究区所需年份内逐日低分辨率地表反射率数据和无云少云条件下中等分辨率遥感数据并进行预处理，得到逐月低分辨率地表反射率数据和中等分辨率地表反射率数据；

步骤S2：对步骤S1所得到的逐月低分辨率地表反射率数据和中等分辨率地表反射率数据进行时空融合，生成融合后的逐月中等分辨率的红光波段和近红外波段地表反射率数据；

步骤S3：计算逐月地形调节植被指数TAVI；

步骤S4：计算逐月植被覆盖度FVC；

步骤S5：计算逐月植被覆盖与管理因子C；

步骤S6：计算逐月降雨侵蚀力因子R；

步骤S7：计算土壤可蚀性因子K、坡长因子L、坡度因子S和水保措施因子P；

步骤S8：计算年土壤侵蚀模数。

进一步的，所述步骤S1的具体步骤如下：

步骤S11：收集研究区所需年份内逐日低分辨率地表反射率数据，以一设定日期数据为主，缺失部分利用当月最大合成法，生成逐月低分辨率地表反射率数据，包括红光波段和近红外波段的地表反射率数据；

步骤S12：收集研究区所需年份内无云少云条件下中等分辨率多光谱遥感数据，对遥感数据进行几何校正、大气校正预处理，生成中等分辨率地表反射率数据，如果是不同传感器的中等分辨率遥感数据，需要进行光谱归一化，让不同传感器相似波段的地表反射率数据具有可比性；

步骤S13：以中等分辨率地表反射率数据为基准，对逐月低分辨率红光波段和近红外波段的地表反射率数据进行几何配准。

进一步的，所述步骤S2的具体步骤如下：

步骤S21：根据所得到的对应月份低分辨率和中等分辨率地表反射率数据对的数量和待计算月份情况，选择不同的时空融合模型；

步骤S22：当仅得到一对低分辨率和中等分辨率地表反射率数据时，基于时空自适应反射率融合模型STARFM，利用该对地表反射率数据和预处理获得的逐月低分辨率地表反射率数据进行时空融合，生成除已有的中等分辨率数据外的其他月份中等分辨率地表反射率数据：

（1）

其中，M为得到的中等分辨率地表反射率数据，M_p为拟融合生成的中等分辨率地表反射率数据，L_o为低分辨率地表反射率数据，为窗口大小；为参与融合的中心像元，为归一化权重系数，由光谱、时间与空间距离三个维度构成；x、y分别代表空间位置横坐标和纵坐标，其下标分别代表不同的位置；t代表时间，t₁为已获取的一个时相数据的获取时间，t_p为拟融合预测的数据时间；

步骤S23：当得到两对或以上低分辨率和中等分辨率地表反射率数据时，选择月份最接近待融合时相的两对数据，和，和，基于增强型时空自适应反射率模型ESTARFM，利用这两对地表反射率数据和预处理获得的逐月低分辨率地表反射率数据进行时空融合，生成除已有的中等分辨率数据外的其他月份中等分辨率地表反射率数据；

先利用公式（1）分别生成拟融合预测数据M_p1和M_p2，

（2）

（3）

再利用如下公式加权计算拟融合预测数据M_p：

（4）

其中，T₁、T₂为时间权重，由如下公式计算：

（5）

步骤S24：根据步骤S22和步骤S23，生成全年1-12月逐月中等分辨率的红光波段和近红外波段地表反射率数据。

进一步的，所述步骤S3中，逐月地形调节植被指数TAVI的计算方法如下：

（6）

其中，表示近红外波段地表反射率，表示红光波段地表反射率，表示研究区红光波段地表反射率的最大值，表示地形调节因子；

利用“极值优化”方法进行的计算，具体步骤如下：

（1）选择样区，检查遥感影像质量，在复杂地形山区选择一定大小的面状区域，确保样区影像“噪声”干扰最小、具有强烈的地形影响，若研究区较小，以整个研究区为样区；

（2）影像分类，应用非监督分类或监督分类方法把遥感影像中的植被划分成阴坡与阳坡两大类；

（3）优化匹配，设计循环程序，令从0开始，以0.001为间隔，依次递增，同时考察TAVI阴坡部分的最大值M_TAVI阴与阳坡部分的最大值M_TAVI阳，当M_TAVI阴与M_TAVI阳相等或近似相等时，退出循环，得到优化结果；若当累增至5时，M_TAVI阴与M_TAVI阳还不满足条件，则返回步骤（1），重新选择样区，或者返回步骤（2），调整分类参数，再重新计算，直至M_TAVI阴与M_TAVI阳满足相等或相近的条件。

进一步的，所述步骤S4中，逐月植被覆盖度FVC的计算方法如下：

（7）

其中，TAVI_S和TAVI_V分别代表全裸土和全植被覆盖的TAVI值。

进一步的，所述步骤S5中，逐月植被覆盖与管理因子C计算方法如下：

（8）。

进一步的，所述步骤S6中，逐月降雨侵蚀力因子R是利用收集的研究区的逐月或逐日的降雨量数据，根据区域特征计算而得。

进一步的，所述步骤S7中，土壤可蚀性因子K、坡长因子L、坡度因子S和水保措施因子P的计算方法如下：

利用土壤属性数据，基于侵蚀—生产力影响计算模型EPIC进行土壤可蚀性因子K的计算：

（9）

其中：SAN为砂粒含量（%）；SIL为粉砂含量（%）；CLA为粘粒含量（%）；c为有机碳含量(%)；SN₁=1-SAN/100；

利用数字高程模型DEM数据，计算坡长因子L和坡度因子S：

（10）

（11）

（12）

其中：θ是坡度，λ 为坡长，m 为坡长指数；

利用中等分辨率遥感数据，分类获得土地利用/覆盖图层，根据土地利用/覆盖类别对水保措施因子P进行赋值。

进一步的，所述步骤S8中，年土壤侵蚀模数计算方法如下：

（13）

其中：R_i为逐月的降雨侵蚀力因子，C_i为逐月的植被覆盖与管理因子。

相较于现有技术，本发明具有以下有益效果：本发明所提出的一种土壤侵蚀方程中植被覆盖与管理因子C改进方法，克服了现有C因子估算方法中没有考虑地形引起的阴坡和阳坡差异，以及不考虑植被季节变化而引起的植被覆盖差异问题，提高了C因子的估算精度与稳定性，进而提高利用土壤流失方程估算土壤侵蚀模数的精度和合理性。

附图说明

图1为本发明一种土壤侵蚀方程中植被覆盖与管理因子C改进方法的流程图；

图2为本发明实施例中逐月C因子空间分布图；

图3为本发明实施例中不同计算方法下的水土流失强度分布图；

图4为本发明实施例中野外实地考察点分布图；

图5为本发明实施例中不同计算方法下的水土流失率曲线图；

图6为本发明实施例中不同计算方法下的单月对比图。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

如图1所示，本发明的一种土壤侵蚀方程中植被覆盖与管理因子C改进方法，包括如下步骤：

步骤S1：获取研究区所需年份内逐日低分辨率地表反射率数据，进行预处理后得到逐月低分辨率的红光波段和近红外波段的地表反射率数据，获取研究区所需年份内无云少云条件下一个时期或若干时期的中等分辨率遥感数据，进行预处理后得到一个时期或若干时期的中等分辨率地表反射率数据；

所述步骤S1的具体步骤如下：

步骤S11：收集研究区所需年份内逐日低空间分辨率（空间分辨率250米）高重访周期的MODIS（中分辨率成像光谱仪）第1-2波段（即红光波段和近红外波段）地表反射率数据，以每月15日左右日期数据为主，缺失部分利用当月最大合成法，生成逐月低分辨率地表反射率数据，包括红光波段和近红外波段的地表反射率数据；

步骤S12：收集研究区所需年份内无云少云条件下一个时期或若干时期的中等分辨率（空间分辨率10-30米）多光谱遥感数据，对遥感数据进行几何校正、大气校正预处理，生成一个时期或若干时期的中等分辨率地表反射率数据，如果是不同传感器的中等分辨率遥感数据，需要进行光谱归一化，让不同传感器相似波段的地表反射率数据具有可比性；

步骤S13：以中等分辨率地表反射率数据为基准，对MODIS逐月低分辨率红光波段和近红外波段的地表反射率数据进行几何配准。

步骤S2：对步骤S1所得到的逐月低分辨率地表反射率数据和一个时期或若干时期中等分辨率地表反射率数据进行时空融合，生成融合后的逐月中等分辨率的红光波段和近红外波段地表反射率数据；

所述步骤S2的具体步骤如下：

步骤S21：根据所得到的对应月份低分辨率和中等分辨率地表反射率数据对的数量和待计算月份情况，选择不同的时空融合模型；

步骤S22：当仅得到一对低分辨率和中等分辨率地表反射率数据时，基于时空自适应反射率融合模型STARFM，利用该对地表反射率数据和预处理获得的逐月低分辨率地表反射率数据进行时空融合，生成除已有的中等分辨率数据外的其他月份中等分辨率地表反射率数据：

（1）

其中，M为得到的中等分辨率地表反射率数据，M_p为拟融合生成的中等分辨率地表反射率数据，L_o为低分辨率地表反射率数据，为窗口大小；为参与融合的中心像元，为归一化权重系数，由光谱、时间与空间距离三个维度构成；x、y分别代表空间位置横坐标和纵坐标，其下标分别代表不同的位置；t代表时间，t₁为已获取的一个时相数据的获取时间，t_p为拟融合预测的数据时间；

步骤S23：当得到两对或以上低分辨率和中等分辨率地表反射率数据时，选择月份最接近待融合时相的两对数据，和，和，基于增强型时空自适应反射率模型ESTARFM，利用这两对地表反射率数据和预处理获得的逐月低分辨率地表反射率数据进行时空融合，生成除已有的中等分辨率数据外的其他月份中等分辨率地表反射率数据；

先利用公式（1）分别生成拟融合预测数据M_p1和M_p2，

（2）

（3）

再利用如下公式加权计算拟融合预测数据M_p：

（4）

其中，T₁、T₂为时间权重，由如下公式计算：

（5）

步骤S24：根据步骤S22和步骤S23，生成全年1-12月逐月中等分辨率的红光波段和近红外波段地表反射率数据。

步骤S3：计算逐月地形调节植被指数TAVI；

逐月地形调节植被指数TAVI的计算方法如下：

（6）

其中，表示近红外波段地表反射率，表示红光波段地表反射率，表示研究区红光波段地表反射率的最大值，表示地形调节因子；

利用“极值优化”方法进行的计算，具体步骤如下：

（1）选择样区，检查遥感影像质量，在复杂地形山区选择一定大小的面状区域，确保样区影像“噪声”干扰最小、具有强烈的地形影响，若研究区较小，以整个研究区为样区；

（2）影像分类，应用非监督分类或监督分类方法把遥感影像中的植被划分成阴坡与阳坡两大类；

（3）优化匹配，设计循环程序，令从0开始，以0.001为间隔，依次递增，同时考察TAVI阴坡部分的最大值M_TAVI阴与阳坡部分的最大值M_TAVI阳，当M_TAVI阴与M_TAVI阳相等或近似相等时，退出循环，得到优化结果；若当累增至5时，M_TAVI阴与M_TAVI阳还不满足条件，则返回步骤（1），重新选择样区，或者返回步骤（2），调整分类参数，再重新计算，直至M_TAVI阴与M_TAVI阳满足相等或相近的条件。

由于光学遥感影像的波段信息是地物反射太阳辐射强度的度量，一年四季受太阳高度角变化（方位角变化较小，可忽略）的影响，相同研究区同一波段在光学遥感影像阴坡与阳坡的信息差异相应变化，总体表现为夏季随太阳高度角增大，差异减小（夏至日最小），冬季随太阳高度角降低，差异增大（冬至日最大），因此也相应夏季减小，冬季增大。

步骤S4：计算逐月植被覆盖度FVC；

利用计算的逐月TAVI，采用像元二分模型计算逐月植被覆盖度FVC(用百分数表示，值为0-100)：

（7）

式中，TAVI_S和TAVI_V分别代表全裸土和全植被覆盖的TAVI值。

步骤S5：计算逐月植被覆盖与管理因子C；

利用计算的逐月FVC，计算逐月植被覆盖与管理因子C：

（8）

步骤S6：计算逐月降雨侵蚀力因子R；

利用收集的逐月或逐日降雨数据，根据不同区域的特点计算逐月降雨侵蚀力因子R。

在本实施例中，以福建省长汀县为研究区，每月的R因子计算如下：

R_i=-2.6398+0.3046P_i，i=1,2…12（9）

式中，P_i（i=1,2…12）是每个月的降雨量。

步骤S7：计算土壤可蚀性因子K、坡长因子L、坡度因子S和水保措施因子P；

利用土壤属性数据，基于侵蚀—生产力影响计算模型EPIC进行土壤可蚀性因子K的计算：

（10）

其中：SAN为砂粒含量（%）；SIL为粉砂含量（%）；CLA为粘粒含量（%）；c为有机碳含量(%)；SN₁=1-SAN/100；

利用数字高程模型DEM数据，计算坡长因子L和坡度因子S：

（11）

（12）

（13）

其中：θ是坡度，λ 为坡长，m 为坡长指数；

利用中等分辨率遥感数据，利用监督分类或非监督分类方法获得土地利用/覆盖图层，根据水土流失特点一般分为林地、农田、水体、不透水面和裸地等几大类。根据土地利用/覆盖类别对水保措施因子P进行赋值。

在本实施例中，林地取1.0，农田0.3，水体0，不透水面0，裸地1.0。

步骤S8：计算年土壤侵蚀模数（逐月累加法）：

（14）

其中：R_i为逐月的降雨侵蚀力因子，C_i为逐月的植被覆盖与管理因子。

目前常用的月年尺度法年土壤侵蚀模数计算方法如下：

（15）

其中C_k为根据所获得影像选用的某一月份中等分辨率遥感数据计算，一般采用归一化植被指数NDVI，而不是采用TAVI；A_k为利用相应月份获得的年土壤侵蚀模数。

在本实施例中，为了下文对改进方法评估的方便，分别利用1-12月的C因子数据计算了12种情况下的年土壤侵蚀模数。如图2所示，为本发明计算生成的长汀县2014年逐月C因子空间分布图，图3为本发明逐月累加法与最差植被覆盖的2月和最好植被覆盖的8月的水土流失强度分布图，是由计算生成土壤侵蚀模数，并根据《土壤侵蚀分类分级标准》（SL190-2007）进行侵蚀强度分级进行分级分类获得的。表1为分别利用1-12月C因子根据公式（15）计算获得的结果统计和利用本发明改进方法公式（14）生成的结果统计。从表1、图2和图3可以看出，在相同年总降雨侵蚀力的条件下（坡长坡度、土壤和水保措施等不变），用不同月份的植被覆盖度估算的C因子代表年尺度的C因子，结果差异非常大。2月年尺度法计算的水土流失总面积比逐月累积法高33.7%，8月年尺度法计算的水土流失总面积比逐月累积法低28.4%。

本发明的逐月累加法年土壤侵蚀模数A_Re的公式中每个月的C因子采用TAVI计算，弥补了常用NDVI无法消除地形影响的不足；由于Ci和Ri都有年内变化，常用公式（15）没有考虑降雨量和植被的年内变化特征，本发明的方法充分考虑了两者年内变化特征及相互匹配问题，理论上更为合理。

在本实施例中，利用野外调查、统计数据及结果时空分布等方面对改进后方法的精度、合理性等方面进行评估。

如图4所示，为长汀县2013-2014年野外实地考察点分布，利用2013-2014年野外实地考察结果验证逐月累加法水土流失估算结果，考察中水土流失强度评级由长汀县水土保持试验站人员现场认定，正确率80.43%。因此以逐月累加法作为参考，分析不同月份C因子估算的月年尺度结果的变化是可行的。

如图5、图6所示，利用统计分析也可以表明，仅利用单一月份影像估算C因子对结果影响较大，从2月到8月，水土流失率逐月下降，8月份以后逐月上升；利用植被覆盖最差月份2月估算的不同强度的水土流失面积与植被覆盖最好月份8月的结果有较大差别，而逐月累加法不同水土流失强度的面积处于两者之间，不会出现极端的情况。本方法较好地克服了植被的季相变化引起的结果差异，与传统以某个月代表全年植被覆盖估算C因子相对，估算结果比较稳定合理。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。