基于耦合协调度的区域风水复合侵蚀风险评价方法

摘要

本发明公开了一种基于耦合协调度的区域风水复合侵蚀风险评价方法，具体包括如下步骤：首先确定研究区域，收集研究区边界shp文件、气象数据、NDVI数据；计算研究区逐月的风蚀气候侵蚀力、降雨侵蚀力；接着计算风蚀气候侵蚀力、降雨侵蚀力、NDVI三维耦合协调度；根据概率统计理论，计算风蚀气候侵蚀力、降雨侵蚀力、NDVI三维耦合协调度的理论累积概率反映风水复合侵蚀风险；根据Mann‑Kendall趋势检验法，即可识别出研究区风水复合侵蚀风险变化趋势的空间分布结果。

说明书

技术领域

本发明属于土壤侵蚀风险评价技术领域，涉及一种基于耦合协调度的区域风水复合侵蚀风险评价方法。

背景技术

土壤侵蚀是一种在全球广泛发生的自然灾害，严重影响人类社会的生产和生活。风力和水力过程一直被认为是土地退化与荒漠化的重要驱动力。水蚀和风蚀是两种最普遍发生的土壤侵蚀类型。特别是在干旱或者半干旱地区，两种侵蚀交互作用形成的复合侵蚀过程和效果不同于单一的风蚀和水蚀。风水复合侵蚀是指风力、水力对同一对象(区域)的共同作用和交替作用，从而产生有别于单一营力的侵蚀、搬运及沉积过程。全球陆面大约17.5％的面积受到风水复合侵蚀的影响。剧烈的风水复合侵蚀导致的土地退化造成了全球性的环境问题，是当前研究的热点问题之一。

风水复合侵蚀主要发生在半干旱地区、海岸、河岸和湖滨地带。在过去的几十年里，关于风水复合侵蚀的研究取得了丰富的研究成果，主要涉及试验小区的风水复合侵蚀观测、不同时空尺度风水复合侵蚀特征分析以及风水复合侵蚀过程中的泥沙搬运、堆积、再分配和对环境的影响。虽然目前的研究取得了许多进展，然而这些成果主要集中在小区尺度的侵蚀泥沙监测和特征分析，以及从风蚀或者水蚀角度对风蚀复合侵蚀地区的侵蚀过程进行研究，然后将两者结合起来进行分析，缺少从区域尺度对风水复合侵蚀力进行耦合评估的理论和方法。因此提出基于耦合协调度的区域风水复合侵蚀风险评价方法是十分重要的。

发明内容

本发明的目的是提供一种基于耦合协调度的区域风水复合侵蚀风险评价方法，采用该方法能够快速准确的识别区域风水复合侵蚀风险的空间分布和变化趋势。

本发明所采用的技术方案是，基于耦合协调度的区域风水复合侵蚀风险评价方法，具体按照以下步骤实施：

步骤1，确定研究区域，收集研究区域边界shp文件、气象站的逐日气象数据和逐月的NDVI栅格数据；

步骤2，利用步骤1收集到的研究区气象站的逐日气象数据，计算每个气象站逐月的风蚀气候侵蚀力，并通过插值和掩膜提取计算获得研究区域的逐月风蚀气候侵蚀力栅格数据；

步骤3，利用步骤1收集到的研究区气象站的逐日气象数据，计算每个气象站逐月的降雨侵蚀力，并通过插值和掩膜提取计算获得研究区域的逐月降雨侵蚀力栅格数据；

步骤4，基于耦合协调度理论，利用步骤1收集到的研究区逐月NDVI数据、步骤2计算得到的研究区逐月风蚀气候侵蚀力栅格数据、步骤3计算得到的研究区逐月降雨侵蚀力栅格数据，计算风蚀气候侵蚀力、降雨侵蚀力、NDVI三维耦合协调度；

步骤5，基于步骤4的计算结果，根据概率统计理论，选择几种概率分布函数作为备用选择，采用线性矩法分别估计几种概率分布函数的参数，再分别利用每种备用概率分布函数计算相应的理论累积概率；同时计算经验累积概率，最后根据经验累积概率和理论累积概率RMSE最小的原则，选出风蚀气候侵蚀力、降雨侵蚀力、NDVI三维耦合协调度的最优理论累积概率；

步骤6，在步骤5选取的最优理论累计概率的基础上，对逐年的理论累积概率取平均值即可获得研究区若干年平均风水复合侵蚀风险空间分布结果；

步骤7，在步骤5选取的最优理论累计概率的基础上，根据Mann-Kendall趋势检验法，即可识别出研究区风水复合侵蚀风险变化趋势的空间分布结果。

本发明的特点还在于：

步骤2的具体过程为：

步骤2.1，根据步骤1收集到的气象数据，利用FAO Penman-Monteith方法计算每个气象站的逐月潜在蒸散发，FAO Penman-Monteith方法的具体计算公式如下：

式中：ET

G、γ、u

G＝0.14(T

式中：T

R

R

R

式中：a为反照率，取值为0.23；R

其中，R

R

式中：n为实际日照时数，h；N为最大可能日照时数，h；n/N为相对日照；a

其中，R

式中：G

式中：e

步骤2.2，根据步骤2.1计算的每个气象站逐月潜在蒸散发结果，计算每个气象站逐月的风蚀气候侵蚀力，具体公式如下：

式中：C为风蚀气候侵蚀力，

步骤2.3，采用ArcGIS软件对气象站逐月的风蚀气候侵蚀力进行克里金插值，得到每月的风蚀气候侵蚀力栅格数据，使用黄河流域边界shp文件在ArcGIS软件中进行裁剪得到逐月黄河流域风蚀气候侵蚀力栅格文件。

步骤3的具体过程为：

步骤3.1，首先根据步骤1收集到的逐日气象数据，利用Zhang的公式计算每个气象站逐月的降雨侵蚀力，具体公式如下：

式中：R

α和β是模型中的两个重要参数，可以通过如下公式进行计算：

β＝0.8363+18.144/P

α＝21.586β

步骤3.2，根据步骤3.1计算的每个气象站每半月的降雨侵蚀力值，求和便可以得到每月的降雨侵蚀力值；

步骤3.3，采用ArcGIS软件对气象站逐月的降雨侵蚀力进行克里金插值，得到每月的降雨侵蚀力栅格数据，使用黄河流域边界shp文件在ArcGIS软件中进行裁剪得到逐月黄河流域降雨侵蚀力栅格文件。

步骤4具体如下：

步骤4.1，根据步骤1收集到的黄河流域逐月NDVI栅格数据、步骤2计算得到的黄河流域逐月风蚀气候侵蚀力栅格数据、步骤3得到的黄河流域逐月降雨侵蚀力数据，分别构建NDVI子系统、风蚀气候侵蚀力子系统、降雨侵蚀力子系统，其中NDVI子系统由每年12个月的月NDVI指标构成，风蚀气候侵蚀力子系统由每年12个月的月风蚀气候侵蚀力指标构成，降雨侵蚀力子系统由每年12个月的月降雨侵蚀力指标构成，上述每个指标需要进行归一化，指标归一化处理方法具体如下：

当u

当u

式中，u

计算每个子系统每项指标在研究期每年的贡献度：

式中，G

计算信息熵：

式中，e

计算权重：

式中，w

计算NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统对总系统有序度的贡献的具体方法如下：

式中，u

步骤4.2，根据多元系统的耦合度模型，,得到黄河流域NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统之间的耦合度模型，如下公式(31)～(33)：

D＝(C·T)

T＝au

式中，D为NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统三维系统的耦合协调度；C为耦合度；u

步骤5的具体过程为：

步骤5.1，根据步骤4中计算得到的NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统三维系统的耦合协调度D，利用如下公式，计算耦合协调度的经验累积概率：

其中，P

步骤5.2，采用线性矩法分别估计几种概率分布函数的参数，具体为：

D为耦合协调度，分布函数为F(D)，且存在反函数G(F)，G(F)又称为分位函数，对于耦合协调度D

D

其中：

步骤5.3，利用4种备用理论概率分布函数指数分布、耿贝尔分布、皮尔逊三型分布、广义极值分布，分别计算相应的理论累积概率：

P

式中，P

步骤5.4，分别计算经验累积概率P

式中，RMSE为均方根误差；n为样本个数；P

步骤7中具体实施如下：

利用Mann-Kendall趋势检验方法对步骤5计算得到的研究区每个栅格点耦合协调度逐年最优理论累积概率P

式中：S为Mann-Kendall趋势检验统计量；Var(S)为S的方差；Z为进行正太标准化处理后的Mann-Kendall趋势检验统计量；P

本发明的有益效果是，本发明基于耦合协调度的区域风水复合侵蚀风险评价方法，计算结果可靠、计算方法简单，具有明确的理论意义，可以充分反映区域风水复合侵蚀风险大小的空间分布，并识别风险可能增加和减少的区域，为区域水土保持规划的制定和实施提供科学的支撑。

附图说明

图1是本发明基于耦合协调度的区域风水复合侵蚀风险评价方法中实现技术流程图；

图2是本发明基于耦合协调度的区域风水复合侵蚀风险评价方法中黄河流域若干年平均耦合协调度空间分布结果示意图；

图3是本发明基于耦合协调度的区域风水复合侵蚀风险评价方法中黄河流域若干年平均风险空间分布结果示意图；

图4是本发明基于耦合协调度的区域风水复合侵蚀风险评价方法中黄河流域若干年平均风险等级空间分布结果示意图；

图5是本发明基于耦合协调度的区域风水复合侵蚀风险评价方法中黄河流域风险变化趋势结果示意图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明基于耦合协调度的区域风水复合侵蚀风险评价方法，具体包括如下步骤：

步骤1，确定研究区域，收集研究区域边界shp文件、气象站的逐日气象数据和逐月的NDVI栅格数据；

步骤2，利用步骤1收集到的研究区气象站的逐日气象数据，计算每个气象站逐月的风蚀气候侵蚀力，并通过插值和掩膜提取计算获得研究区域的逐月风蚀气候侵蚀力栅格数据；步骤2的具体过程为：

步骤2.1，首先根据步骤1收集到的研究区逐日气象数据，利用FAO Penman-Monteith方法计算每个气象站的逐月潜在蒸散发，FAO Penman-Monteith方法的具体计算公式如下：

式中：ET

G、γ、u

G＝0.14(T

式中：T

R

R

R

式中：a为反照率，取值为0.23；R

其中，R

R

式中：n为实际日照时数，h；N为最大可能日照时数，h；n/N为相对日照；a

其中，R

式中：G

式中：e

步骤2.2，根据步骤2.1计算的每个气象站逐月潜在蒸散发结果计算每个气象站逐月的风蚀气候侵蚀力，具体公式如下：

式中：C为风蚀气候侵蚀力，

采用ArcGIS软件对气象站逐月的风蚀气候侵蚀力进行克里金插值，得到每月的风蚀气候侵蚀力栅格数据，使用研究区边界shp文件在ArcGIS软件中进行裁剪得到逐月研究区风蚀气候侵蚀力栅格文件；

步骤3，利用步骤1收集到的研究区气象站的逐日气象数据，计算每个气象站逐月的降雨侵蚀力，并通过插值和掩膜提取计算获得研究区域的逐月降雨侵蚀力栅格数据；步骤3的具体步骤为：

步骤3.1，首先根据步骤1收集到的研究区逐日气象数据，利用Zhang的公式计算每个气象站逐月的降雨侵蚀力，具体公式如下：

式中：R

α和β是模型中的两个重要参数，可以通过如下公式进行计算：

β＝0.8363+18.144/P

α＝21.586β

步骤3.2，根据步骤3.1计算的每个气象站每半月的降雨侵蚀力值，求和便可以得到每月的降雨侵蚀力值。采用ArcGIS软件对气象站逐月的降雨侵蚀力进行克里金插值，得到每月的降雨侵蚀力栅格数据，使用研究区边界shp文件在ArcGIS软件中进行裁剪得到逐月研究区降雨侵蚀力栅格文件；

步骤4，基于耦合协调度理论，利用步骤1收集到的研究区逐月NDVI数据、步骤2计算得到的研究区逐月风蚀气候侵蚀力栅格数据、步骤3计算得到的研究区逐月降雨侵蚀力栅格数据，计算风蚀气候侵蚀力、降雨侵蚀力、NDVI三维耦合协调度；

步骤4的具体步骤为：步骤4.1，首先根据步骤1收集到的研究区逐月NDVI栅格数据、步骤2计算得到的研究区逐月风蚀气候侵蚀力栅格数据、步骤3得到的研究区逐月降雨侵蚀力数据，分别构建NDVI子系统、风蚀气候侵蚀力子系统、降雨侵蚀力子系统，其中NDVI子系统由每年12个月的月NDVI指标构成，风蚀气候侵蚀力子系统由每年12个月的月风蚀气候侵蚀力指标构成，降雨侵蚀力子系统由每年12个月的月降雨侵蚀力指标构成。上述每个指标需要进行归一化，指标归一化处理方法具体如下：

当u

式中u

计算每个子系统每项指标在研究期每年的贡献度：

式中，G

式中，e

计算权重：

式中，w

计算NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统对总系统有序度的贡献的具体方法如下：

式中u

步骤4.2，根据多元系统的耦合度模型，很容易得到NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统之间的耦合度模型，可以表示成以下形式：

式中，D为NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统三维系统的耦合协调度；C为耦合度；u

步骤5，基于步骤4的计算结果，根据概率统计理论，选择几种概率分布函数作为备用选择，采用线性矩法分别估计几种概率分布函数的参数，再分别利用每种备用概率分布函数计算相应的理论累积概率；同时计算经验累积概率，最后根据经验累积概率和理论累积概率RMSE最小的原则，选出风蚀气候侵蚀力、降雨侵蚀力、NDVI三维耦合协调度的最优理论分布函数及其理论累积概率；

步骤5的具体步骤为：

其中，P

步骤5.2，采用线性矩法分别估计几种概率分布函数的参数，具体为：D为耦合协调度，分布函数为F(D)，且存在反函数G(F)，G(F)又称为分位函数，对于耦合协调度D

其中：

步骤5.3，利用4种备用理论概率分布函数指数分布、耿贝尔分布、皮尔逊三型分布、广义极值分布，分别计算相应的理论累积概率：

P

式中，P

步骤5.4，分别计算经验累积概率P

式中，RMSE为均方根误差；n为样本个数；P

步骤7中具体实施如下：利用Mann-Kendall趋势检验方法对步骤5计算得到的研究区每个栅格点耦合协调度逐年最优理论累积概率P

式中：S为Mann-Kendall趋势检验统计量；Var(S)为S的方差；Z为进行正太标准化处理后的Mann-Kendall趋势检验统计量；P

实施例

如图1所示，以黄河流域为例，技术流程如图1所示，具体按照以下步骤实施：步骤1，首先收集基础数据：收集黄河流域及其附近区域1982-2015年289个气象站逐日降水量、风速、最高气温、最低气温、平均相对湿度、日照时数；1982-2015年逐月NDVI数据；黄河流域边界shp文件；

步骤2，实施如下；步骤2.1，首先根据步骤1收集到的黄河流域289个气象站1982-2015年逐日气象数据，利用FAO Penman-Monteith方法计算每个气象站的逐月潜在蒸散发，FAO Penman-Monteith方法的具体计算公式如下：

式中：ET

G、γ、u

G＝0.14(T

式中：T

R

R

式中：a为反照率，取值为0.23；R

式中：n为实际日照时数，h；N为最大可能日照时数，h；n/N为相对日照；a

式中：G

式中：e

步骤2.2，根据步骤2.1计算的每个气象站逐月潜在蒸散发结果，计算每个气象站逐月的风蚀气候侵蚀力，具体公式如下：

式中：C为风蚀气候侵蚀力，

步骤3，具体实施如下；步骤3.1，首先根据步骤1收集到的黄河流域289个气象站1982-2015年逐日气象数据，利用Zhang的公式计算每个气象站逐月的降雨侵蚀力，具体公式如下：

式中：R

β＝0.8363+18.144/P

步骤3.2，根据步骤3.1计算的每个气象站每半月的降雨侵蚀力值，求和便可以得到每月的降雨侵蚀力值。采用ArcGIS软件对气象站逐月的降雨侵蚀力进行克里金插值，得到每月的降雨侵蚀力栅格数据，使用黄河流域边界shp文件在ArcGIS软件中进行裁剪得到逐月黄河流域降雨侵蚀力栅格文件；

步骤4，具体实施如下；步骤4.1，首先根据步骤1收集到的黄河流域1982-2015年逐月NDVI栅格数据、步骤2计算得到的黄河流域逐月风蚀气候侵蚀力栅格数据、步骤3得到的黄河流域逐月降雨侵蚀力数据，分别构建NDVI子系统、风蚀气候侵蚀力子系统、降雨侵蚀力子系统，其中NDVI子系统由每年12个月的月NDVI指标构成，风蚀气候侵蚀力子系统由每年12个月的月风蚀气候侵蚀力指标构成，降雨侵蚀力子系统由每年12个月的月降雨侵蚀力指标构成。上述每个指标需要进行归一化，指标归一化处理方法具体如下：

当u

式中，e

步骤4.2，根据多元系统的耦合度模型，很容易得到黄河流域NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统之间的耦合度模型，可以表示成以下形式：

步骤5，具体实施如下；步骤5.1，根据步骤4中计算得到的NDVI子系统、风蚀气候侵蚀力子系统和降雨侵蚀力子系统三维系统的耦合协调度D，利用如下公式，计算耦合协调度的经验累积概率：

步骤5.2，采用线性矩法分别估计几种概率分布函数的参数，具体为：D为耦合协调度，分布函数为F(D)，且存在反函数G(F)，G(F)又称为分位函数，对于耦合协调度D

步骤5.3，利用4种备用理论概率分布函数指数分布、耿贝尔分布、皮尔逊三型分布、广义极值分布，分别计算相应的理论累积概率：

P

式中，P

步骤5.4，分别计算经验累积概率P

式中，RMSE为均方根误差；n为样本个数；P

步骤6，在步骤5选取的最优理论分布函数的理论累计概率的基础上，利用ArcGIS软件对逐年的理论累积概率取平均值即可获得黄河流域1982-2015年34年平均风水复合侵蚀风险空间分布结果；黄河流域1982-2015年风水复合侵蚀风险空间分布如图3所示。通过ArcGIS软件自然断点法分为3类即可获得黄河流域风水复合侵蚀风险等级空间分布，如图4所示。

步骤7，利用Mann-Kendall趋势检验方法对步骤5计算得到的黄河流域每个栅格点耦合协调度逐年最优理论累积概率P

式中：S为Mann-Kendall趋势检验统计量；Var(S)为S的方差；Z为进行正太标准化处理后的Mann-Kendall趋势检验统计量；P