CN115099453B - 一种多变量栅格化卫星遥感综合干旱风险评估方法 - Google Patents

Abstract

本发明公开了一种多变量栅格化卫星遥感综合干旱风险评估方法，包括以下步骤：获取并处理目标研究区域内各栅格的遥感数据，遥感数据包括卫星降水量、卫星增强型植被指数和卫星地表温度，进而通过标准正态化方法得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；提取综合干旱风险评估要素；计算各综合干旱风险评估要素的经验频率；构建综合干旱风险评估指标CDRA。本发明将卫星遥感与干旱多维度量结合在一起，实现了由点到面对干旱灾害进行实时精细监测，不仅可以从多角度全面的描述评估干旱而且弥补了离散站点监测空间和时间不连续的问题，对大规模干旱风险评估和生态保护具有指导意义。

Description

一种多变量栅格化卫星遥感综合干旱风险评估方法

技术领域

本发明属于环境测量领域，具体涉及一种多变量栅格化卫星遥感综合干旱风险评估方法。

背景技术

由于干旱具有形成慢、持续时间长的特点，所以经常不被人们注意，但是其造成的伤害并不亚于地震、洪涝等自然灾害。准确快速的监测到旱情，并对旱情风险进行有效准确的评估已成为亟须研究的热点科学问题。干旱指数构建是量化干旱研究的基础，通过捕捉干旱的开始、持续时间和结束等特征，其可对干旱特征进行客观的时空分析和监测，从而达到度量、对比和综合的研究目的。目前，气象干旱常用的指数有标准化降水指数(SPI)、自校正帕默尔干旱指数(sc-PDSI)等；农业干旱常用的指数有标准化土壤湿度指数(SSI)、卫星增强型植被指数(EVI)、植被健康指数(VHI)等。然而单一指数往往不能综合多种要素，从而无法综合地从气象、农业、生态等领域出发全方面地监测干旱，且干旱等级与干旱风险之间没有清晰明确的关系。已有的综合干旱风险评估多基于地面站点观测数据或陆面模型模拟数据。站点观测数据虽然可以准确表征站点周边的干旱状况，但是观测站点空间密度有限，仅靠常规站点的观测还不足以了解干旱发生发展过程的全貌；陆面模型模拟的数据虽然可以弥补站点数据的不足，但是模拟结果不确定性较高，即使基于相同的气象驱动数据，不同水文模型的模拟结果相差较大，影响了干旱监测精度。

发明内容

为了综合表征气象-生态农业干旱，并评估综合干旱事件所造成的风险，本发明提供了一种多变量栅格化卫星遥感综合干旱风险评估方法，基于卫星降水量、卫星增强型植被指数和卫星地表温度等遥感数据产品，构建了一种将卫星遥感与干旱多维度量结合在一起的栅格化卫星遥感综合干旱风险评估方法。

为解决上述技术问题，本发明采用以下技术方案：

一种多变量栅格化卫星遥感综合干旱风险评估方法，针对目标研究区域，执行以下步骤，实现对目标研究区域内各栅格的综合干旱事件进行综合干旱风险评估：

步骤1，基于目标研究区域的栅格划分，针对目标研究区域在预设历史时长段内的遥感数据，基于以该预设历史时长段内的每月为起点、预设个数月份为跨度的时间段，获得该预设历史时长段内各时间段的遥感数据，包含时间段内的卫星降水量、卫星增强型植被指数和卫星地表温度，进而得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；

步骤2，基于各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，分别针对各栅格在预设个数连续时间段的卫星遥感综合干旱指数CRSDI均小于预设卫星遥感综合干旱指数CRSDI阈值，表示该栅格发生了一次综合干旱事件，提取各栅格在预设历史时长段内分别对应各综合干旱事件的各预设综合干旱风险评估要素；

步骤3，针对各栅格对应各综合干旱事件的各预设综合干旱风险评估要素，逐栅格计算各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率；

步骤4，基于各栅格对应各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率，构建综合干旱风险评估指标CDRA，进而实现对目标研究区域内各栅格的综合干旱事件进行综合干旱风险评估。

作为本发明的一种优选技术方案，所述步骤1的具体过程如下：

步骤1.1，获取目标研究区域各栅格在预设历史时长段内的遥感数据，包括月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度，并对目标研究区域的遥感数据进行插值处理统一分辨率；

步骤1.2，针对月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度，基于以每月为起点，预设个数月份为跨度的时间段，获得预设历史时长段内各时间段的遥感数据，进而计算获得各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数、以及卫星地表温度的边际概率值；

步骤1.3，基于各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的边际概率值，计算获得各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值；

步骤1.4，基于各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值，计算获得各栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值；并基于各栅格在各时间段分别对应的累积概率值，得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI。

作为本发明的一种优选技术方案，所述步骤1.2中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的边际概率值：

步骤1.2.1，针对该栅格在预设历史时长段内的月尺度的卫星降水量，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有f项，将该栅格的卫星降水量时间段序列按照各时间段结束月份划分为12个独立的子序列，即结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星降水量均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量边际概率值F(x)：

其中，x指代任意一时间段的卫星降水量，各子序列的长度均为f_a，i指代x在对应的子序列中排序为第i个；

步骤1.2.2，对该栅格在预设历史时长段内的月尺度的卫星增强型植被指数，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有m项，将该栅格的卫星增强型植被指数时间段序列按照各时间段结束月份划分为12个独立的子序列，即结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星增强型植被指数均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星增强型植被指数边际概率值F(y)：

其中，y指代任意一时间段的卫星增强型植被指数，各子序列的长度均为m_a，j指代y在对应的子序列中排序为第j个；

步骤1.2.3，对该栅格在预设历史时长段内的月尺度的卫星地表温度，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有r项，将该栅格的卫星地表温度时间段序列按照各时间段结束月份划分为12个独立的子序列，即结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星增强型植被指数均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星地表温度边际概率值F(z)：

其中，z指代任意一时间段的卫星增强型植被指数，各子序列的长度均为r_a，k指代z在对应的子序列中排序为第k个；

上述步骤中，卫星降水量时间段序列、卫星增强型植被指数时间段序列和卫星地表温度时间段序列长度一致，即f＝m＝r，f_a＝m_a＝r_a。

作为本发明的一种优选技术方案，所述步骤1.3中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值：

步骤1.3.1，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值F(x,y)：

F(x,y)＝P(X≤x,Y≤y)＝C(F(x),F(y))＝C(u,v)

式中，x指代任一时间段的卫星降水量；y指代任一时间段的卫星增强型植被指数；F(x)＝u为卫星降水量x所在时间段的边际概率值；F(y)＝v为卫星增强型植被指数y所在时间段的边际概率值；C为Copula函数；

其中，Copula函数为对称Archimedean Copulas家族中的Gumbel Copula函数，

式中，θ₁为参数，其取值范围为[1,∞)；exp为以自然常数e为底的指数函数，ln为以常数e为底数的对数函数；

步骤1.3.2，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值F(x,y,z)：

F(x,y,z)＝P(X≤x,Y≤y,Z≤z)＝C(F(x),F(y),F(z))＝C(u,v,w)

式中，z指代任一时间段的卫星地表温度；F(z)＝w为卫星地表温度z所在时间段的边际概率值；

其中，Copula函数为非对称Archimedean Copulas中的Gumbel Copula函数，

式中，θ₂为参数，其取值范围为[1,∞)。

作为本发明的一种优选技术方案，所述步骤1.4中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值，并基于各栅格在各时间段分别对应的累积概率值，得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；

步骤1.4.1，采用以下公式计算该栅格在各时间段的二元联合概率值与三元联合概率值之差，即事件(X≤x,Y≤y,Z>z)的联合概率q：

P(X≤x,Y≤y,Z>z)＝P(X≤x,Y≤y)-P(X≤x,Y≤y,Z≤z)

＝C(F(x),F(y))-C(F(x),F(y),F(z))＝C(u,v)-C(u,v,w)＝q

步骤1.4.2，将各时间段的联合概率q按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段的二元联合概率值与三元联合概率值之差的累积概率值F(q)：

其中，q指代任意一时间段的二元联合概率值与三元联合概率值之差，I指代q的排序为第I个，N指代时间段序列长度，即时间段总数；

上述步骤中，二元联合概率值与三元联合概率值之差所得对应的时间段序列与卫星降水量时间段序列、卫星增强型植被指数时间段序列和卫星地表温度时间段序列长度一致，即N＝f＝m＝r；

步骤1.4.3，基于该栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值，通过标准正态化方法得到该栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，其计算公式为：

式中，

为标准正态分布。

作为本发明的一种优选技术方案，所述步骤1.3中针对计算卫星降水量和卫星增强型植被指数的二元联合概率值，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值，分别针对各栅格，通过以下方法，分别选取计算卫星降水量和卫星增强型植被指数的二元联合概率值对应的最优Copula函数，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值对应的最优Copula函数：

基于Archimedean Copulas中的Gumbel Copula函数、Clayton Copula函数和Frank Copula函数，以及Meta-elliptic Copulas中的Student t Copula函数和GaussianCopula函数，基于该5种Copula函数，通过以下公式，以赤池信息量准则AIC、贝叶斯信息准则BIC和方根误差准则RMSE三项中至少两项值最小作为最优的判断标准，从该5种Copula函数中分别选取计算计算卫星降水量和卫星增强型植被指数的二元联合概率值对应的最优Copula函数，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值对应的最优Copula函数；

式中，MSE为均方误差，N为样本容量，即时间段总数；b为Copula函数的参数个数；Pe_n为第n个时间段的联合分布经验频率，P_n为第n个时间段Copula函数计算的联合概率值；

第n个时间段的卫星降水量和卫星增强型植被指数的联合分布经验频率H(x_n,y_n)为

其中，x_n指代第n个时间段的卫星降水量，y_n指代第n个时间段的卫星增强型植被指数，N_b指代各时间段的样本中同时满足X≤x_n,Y≤y_n的样本个数；

第n个时间段的变量卫星降水量、卫星增强型植被指数和卫星地表温度的联合分布经验频率H(x_n,y_n,z_n)为：

其中，z_n指代第n个时间段的卫星地表温度，N_t指代各时间段的样本中同时满足X≤x_n,Y≤y_n,Z≤z_n的样本个数。

作为本发明的一种优选技术方案，所述步骤2中基于预设卫星遥感综合干旱指数CRSDI＝-0.5为阈值，分别针对各栅格在预设个数连续时间段的卫星遥感综合干旱指数CRSDI均小于预设卫星遥感综合干旱指数CRSDI阈值，表示该栅格发生了一次综合干旱事件，提取各栅格在预设历史时长段内分别对应各综合干旱事件的各预设综合干旱风险评估要素；各预设综合干旱风险评估要素包括干旱历时、干旱烈度、干旱烈度峰值及干旱后生态恢复期；

干旱历时是指连续小于CRSDI阈值的时间段个数；

干旱烈度为连续小于CRSDI阈值的各时间段CRSDI累加和的绝对值；

干旱烈度峰值指在连续小于CRSDI阈值的各时间段中最小CRSDI的绝对值；

干旱后生态恢复期获取过程如下：

从综合干旱事件结束后的第一个月开始，逐月计算植被生长变化E，计算公式为：

式中，E_h为第h个月的植被生长变化，y_h为综合干旱事件结束后第h个月的卫星增强型植被指数，

为y_h所对应月份卫星增强型植被指数在预设历史时长段内的年平均值；E_h<0表明生态未恢复，干旱后生态恢复期为综合干旱事件结束后首次E_h≥0时h所对应的月份个数。/>

作为本发明的一种优选技术方案，所述步骤3中分别针对目标研究区域的各栅格，执行以下步骤，逐栅格计算各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率；

将该栅格在预设历史时长段内经历的各综合干旱事件的干旱历时由小到大的次序排列、干旱烈度由小到大的次序排列、干旱烈度峰值由小到大的次序排列、干旱后生态恢复期由小到大的次序排列；

该栅格在任一次综合干旱事件中干旱历时的经验频率D为：

其中，g为该栅格在预设历史时长段内共经历的综合干旱事件总次数，该次综合干旱事件的干旱历时排序为第J个；

该栅格在任一次综合干旱事件中干旱烈度的经验频率S为：

其中，该次综合干旱事件的干旱烈度排序为第K个；

该栅格在任一次综合干旱事件中干旱烈度峰值的经验频率P_E为：

其中，该次综合干旱事件的干旱烈度峰值排序为第H个；

该栅格在任一次综合干旱事件中干旱后生态恢复期的经验频率T为：

其中，该次综合干旱事件的干旱后生态恢复期排序为第L个。

作为本发明的一种优选技术方案，所述步骤4中，

基于各栅格对应各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率，基于权重组合法得到各次综合干旱事件的风险评估指标CDRA，其计算公式为：

CDRA＝α×D+β×S+γ×P_E+δ×T

式中，α、β、γ和δ分别为D、S、P_E和T所对应的权重，α+β+γ+δ＝1。

本发明的有益效果：本发明提供了一种多变量栅格化卫星遥感综合干旱风险评估方法，首先获取并处理研究区内的遥感数据，包括月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度；分别逐栅格计算卫星降水量、卫星增强型植被指数和卫星地表温度序列的边际概率值；分别逐栅格计算卫星降水量和卫星增强型植被指数的二元联合概率值及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值；逐栅格计算二元联合概率值与三元联合概率值之差所得序列的累积概率值；进而通过标准正态化方法得到卫星遥感综合干旱指数CRSDI；然后提取综合干旱风险评估要素；计算各综合干旱风险评估要素的经验频率；最终构建综合干旱风险评估指标CDRA。本发明将卫星遥感与干旱多维度量结合在一起，而遥感数据能够提供覆盖范围大、时空连续的气象、水文和植被信息，具有重访周期短、空间分辨率高、数据获取方便、资料客观等优点，实现了由点到面对干旱灾害进行实时精细监测，不仅可以从多角度全面的描述评估干旱而且弥补了离散站点监测空间和时间不连续的问题，对大规模干旱风险评估和生态保护具有指导意义，进而减少干旱灾害产生的巨大社会效益。

附图说明

图1为本发明一种多变量栅格化卫星遥感综合干旱风险评估方法的流程图。

具体实施方式

下面的实施例可使本专业技术人员更全面地理解本发明，但不以任何方式限制本发明。

为了综合表征气象-生态农业干旱，并评估综合干旱事件所造成的风险，本发明基于卫星降水量、卫星增强型植被指数和卫星地表温度等遥感数据产品，构建了一种将卫星遥感与干旱多维度量结合在一起的栅格化卫星遥感综合干旱风险评估方法。本发明提供一种多变量栅格化卫星遥感综合干旱风险评估方法，本实施例以青藏高原为例，卫星降水量采用IMERG V06产品数据，卫星增强型植被指数和卫星地表温度采用MODIS V6产品数据，针对目标研究区域，执行以下步骤，如图1所示，实现对目标研究区域内各栅格的综合干旱事件进行综合干旱风险评估：

步骤1，基于目标研究区域的栅格划分，针对目标研究区域在预设历史时长段内的遥感数据，该预设时长段越长，所得数据越精确，最好是30年左右，基于以该预设历史时长段内的每月为起点、预设个数月份为跨度的时间段，如1、3和6个月等，以3个月为例，则各时间段分别为1-3月、2-4月、3-5月等，获得该预设历史时长段内各时间段的遥感数据，包含时间段内的卫星降水量、卫星增强型植被指数和卫星地表温度，进而得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，其中，每个时间段的数据反映的是每个时间段结束月份的情况；

所述步骤1的具体过程如下：

步骤1.1，获取目标研究区域各栅格在预设历史时长段内的遥感数据，包括月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度，并对目标研究区域的遥感数据进行插值处理统一分辨率，采用双线性插值法将分辨率统一为0.1°；

步骤1.2，针对月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度，基于以每月为起点，预设个数月份为跨度的时间段，获得预设历史时长段内各时间段的遥感数据，各时间段的遥感数据即为每月遥感数据之和，进而计算获得各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数、以及卫星地表温度的边际概率值；

所述步骤1.2中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的边际概率值：

步骤1.2.1，针对该栅格在预设历史时长段内的月尺度的卫星降水量，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有f项，将该栅格的卫星降水量时间段序列按照各时间段结束月份划分为12个独立的子序列，即结束月份相同的时间段在同一个子序列中，即每年的1-3月归集在一个子序列中，每年的2-4月归集在一个子序列中，以此类推；并分别针对各子序列，将子序列中各时间段的卫星降水量均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量边际概率值F(x)：

其中，x指代任意一时间段的卫星降水量，各子序列的长度均为f_a，i指代x在对应的子序列中排序为第i个；

步骤1.2.2，对该栅格在预设历史时长段内的月尺度的卫星增强型植被指数，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有m项，将该栅格的卫星增强型植被指数时间段序列按照各时间段结束月份划分为12个独立的子序列，即结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星增强型植被指数均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星增强型植被指数边际概率值F(y)：

其中，y指代任意一时间段的卫星增强型植被指数，各子序列的长度均为m_a，j指代y在对应的子序列中排序为第j个；

步骤1.2.3，对该栅格在预设历史时长段内的月尺度的卫星地表温度，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有r项，将该栅格的卫星地表温度时间段序列按照各时间段结束月份划分为12个独立的子序列，即结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星增强型植被指数均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星地表温度边际概率值F(z)：

其中，z指代任意一时间段的卫星增强型植被指数，各子序列的长度均为r_a，k指代z在对应的子序列中排序为第k个；

上述步骤中，卫星降水量时间段序列、卫星增强型植被指数时间段序列和卫星地表温度时间段序列长度一致，即f＝m＝r，f_a＝m_a＝r_a。

步骤1.3，基于各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的边际概率值，计算获得各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值；

所述步骤1.3中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值：

步骤1.3.1，对于所有x,y∈R,R∈(-∞,∞)，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值F(x,y)：

F(x,y)＝P(X≤x,Y≤y)＝C(F(x),F(y))＝C(u,v)

式中，x指代任一时间段的卫星降水量；y指代任一时间段的卫星增强型植被指数；F(x)＝u为卫星降水量x所在时间段的边际概率值；F(y)＝v为卫星增强型植被指数y所在时间段的边际概率值；C为Copula函数；

其中，Copula函数为对称Archimedean Copulas家族中的Gumbel Copula函数，

式中，θ₁为参数，其取值范围为[1,∞)；exp为以自然常数e为底的指数函数，ln为以常数e为底数的对数函数，采用极大似然估计法进行参数估计，计算出θ₁；

步骤1.3.2，对于所有x,y,z∈R,R∈(-∞,∞)，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值F(x,y,z)：

F(x,y,z)＝P(X≤x,Y≤y,Z≤z)＝C(F(x),F(y),F(z))＝C(u,v,w)

式中，z指代任一时间段的卫星地表温度；F(z)＝w为卫星地表温度z所在时间段的边际概率值；

其中，Copula函数为非对称Archimedean Copulas中的Gumbel Copula函数，

式中，θ₂为参数，其取值范围为[1,∞)，采用极大似然估计法进行参数估计，计算出θ₂。

另一方面，所述步骤1.3中针对计算卫星降水量和卫星增强型植被指数的二元联合概率值，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值，分别针对各栅格，通过以下方法，分别选取计算卫星降水量和卫星增强型植被指数的二元联合概率值对应的最优Copula函数，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值对应的最优Copula函数：

基于Archimedean Copulas中的Gumbel Copula函数、Clayton Copula函数和Frank Copula函数，以及Meta-elliptic Copulas中的Student t Copula函数和GaussianCopula函数，基于该5种Copula函数，通过以下公式，以赤池信息量准则AIC、贝叶斯信息准则BIC和方根误差准则RMSE三项中至少两项值最小作为最优的判断标准，从该5种Copula函数中分别选取计算计算卫星降水量和卫星增强型植被指数的二元联合概率值对应的最优Copula函数，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值对应的最优Copula函数；

式中，MSE为均方误差，N为样本容量，即时间段总数；b为Copula函数的参数个数；Pe_n为第n个时间段的联合分布经验频率，P_n为第n个时间段Copula函数计算的联合概率值；

第n个时间段的卫星降水量和卫星增强型植被指数的联合分布经验频率H(x_n,y_n)为

其中，x_n指代第n个时间段的卫星降水量，y_n指代第n个时间段的卫星增强型植被指数，N_b指代各时间段的样本中同时满足X≤x_n,Y≤y_n的样本个数；

第n个时间段的变量卫星降水量、卫星增强型植被指数和卫星地表温度的联合分布经验频率H(x_n,y_n,z_n)为：

其中，z_n指代第n个时间段的卫星地表温度，N_t指代各时间段的样本中同时满足X≤x_n,Y≤y_n,Z≤z_n的样本个数。

理论上每一个时间段均可遴选出各自最优偏态分布曲线及参数，但算法复杂且工作量巨大，同时通用性显著降低。考虑到同一变量具有相对稳定的统计特征，选择对不同时间尺度都拟合较好的同一种分布曲线但不同参数组合来处理时间段和季节性差异的影响。

步骤1.4，基于各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值，计算获得各栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值；并基于各栅格在各时间段分别对应的累积概率值，得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI。

所述步骤1.4中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值，并基于各栅格在各时间段分别对应的累积概率值，得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；

步骤1.4.1，采用以下公式计算该栅格在各时间段的二元联合概率值与三元联合概率值之差，即事件(X≤x,Y≤y,Z>z)的联合概率q：

P(X≤x,Y≤y,Z>z)＝P(X≤x,Y≤y)-P(X≤x,Y≤y,Z≤z)

＝C(F(x),F(y))-C(F(x),F(y),F(z))＝C(u,v)-C(u,v,w)＝q

步骤1.4.2，将各时间段的联合概率q按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段的二元联合概率值与三元联合概率值之差的累积概率值F(q)：

/>

其中，q指代任意一时间段的二元联合概率值与三元联合概率值之差，I指代q的排序为第I个，N指代时间段序列长度，即时间段总数；

上述步骤中，二元联合概率值与三元联合概率值之差所得对应的时间段序列与卫星降水量时间段序列、卫星增强型植被指数时间段序列和卫星地表温度时间段序列长度一致，即N＝f＝m＝r；

步骤1.4.3，基于该栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值，通过标准正态化方法得到该栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，其计算公式为：

式中，

为标准正态分布。

基于卫星遥感综合干旱指数CRSDI，将干旱划分为无旱、轻旱、中旱、重旱和特旱等不同的等级。-1<CRSDI≤-0.5，划分为轻旱；-1.5<CRSDI≤-1，划分为中旱；-2<CRSDI≤-1.5，划分为重旱；CRSDI≤-2，划分为特旱。

步骤2，基于各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，分别针对各栅格在预设个数连续时间段的卫星遥感综合干旱指数CRSDI均小于预设卫星遥感综合干旱指数CRSDI阈值，表示该栅格发生了一次综合干旱事件，提取各栅格在预设历史时长段内分别对应各综合干旱事件的各预设综合干旱风险评估要素；

所述步骤2中基于游程理论，预设卫星遥感综合干旱指数CRSDI＝-0.5为阈值，分别针对各栅格在预设个数连续时间段的卫星遥感综合干旱指数CRSDI均小于预设卫星遥感综合干旱指数CRSDI阈值，表示该栅格发生了一次综合干旱事件，即当CRSDI序列在一定时段内连续小于阈值时则出现负游程，表明发生了综合干旱事件，提取各栅格在预设历史时长段内分别对应各综合干旱事件的各预设综合干旱风险评估要素；各预设综合干旱风险评估要素包括干旱历时、干旱烈度、干旱烈度峰值及干旱后生态恢复期；

干旱历时是指连续小于CRSDI阈值的时间段个数；

干旱烈度为连续小于CRSDI阈值的各时间段CRSDI累加和的绝对值；

干旱烈度峰值指在连续小于CRSDI阈值的各时间段中最小CRSDI的绝对值；

干旱后生态恢复期获取过程如下：

从综合干旱事件结束后的第一个月开始，逐月计算植被生长变化E，计算公式为：

式中，E_h为第h个月的植被生长变化，y_h为综合干旱事件结束后第h个月的卫星增强型植被指数，

为y_h所对应月份卫星增强型植被指数在预设历史时长段内的年平均值；E_h<0表明生态未恢复，干旱后生态恢复期为综合干旱事件结束后首次E_h≥0时h所对应的月份个数。

步骤3，针对各栅格对应各综合干旱事件的各预设综合干旱风险评估要素，逐栅格计算各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率；

所述步骤3中分别针对目标研究区域的各栅格，执行以下步骤，逐栅格计算各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率；

将该栅格在预设历史时长段内经历的各综合干旱事件的干旱历时由小到大的次序排列、干旱烈度由小到大的次序排列、干旱烈度峰值由小到大的次序排列、干旱后生态恢复期由小到大的次序排列；

该栅格在任一次综合干旱事件中干旱历时的经验频率D为：

其中，g为该栅格在预设历史时长段内共经历的综合干旱事件总次数，该次综合干旱事件的干旱历时排序为第J个；

该栅格在任一次综合干旱事件中干旱烈度的经验频率S为：

其中，该次综合干旱事件的干旱烈度排序为第K个；

该栅格在任一次综合干旱事件中干旱烈度峰值的经验频率P_E为：

其中，该次综合干旱事件的干旱烈度峰值排序为第H个；

该栅格在任一次综合干旱事件中干旱后生态恢复期的经验频率T为：

其中，该次综合干旱事件的干旱后生态恢复期排序为第L个。

步骤4，基于各栅格对应各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率，构建综合干旱风险评估指标CDRA，进而实现对目标研究区域内各栅格的综合干旱事件进行综合干旱风险评估，可以根据CDRA指标的大小来评估，越大说明干旱风险越大。

所述步骤4中，基于各栅格对应各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率，基于权重组合法得到各次综合干旱事件的风险评估指标CDRA，其计算公式为：

CDRA＝α×D+β×S+γ×P_E+δ×T

式中，α、β、γ和δ分别为D、S、P_E和T所对应的权重，α+β+γ+δ＝1，令α＝β＝γ＝δ＝0.25。

已有的综合干旱风险评估多基于地面站点观测数据或陆面模型模拟数据。站点观测数据虽然可以准确表征站点周边的干旱状况，但是观测站点空间密度有限，仅靠常规站点的观测还不足以了解干旱发生发展过程的全貌；陆面模型模拟的数据虽然可以弥补站点数据的不足，但是模拟结果不确定性较高，即使基于相同的气象驱动数据，不同水文模型的模拟结果相差较大，影响了干旱监测精度。随着遥感技术的发展，基于可见光、红外和微波波段等遥感数据的干旱监测方法逐渐发展起来，其能够提供覆盖范围大、时空连续的气象、水文和植被信息，具有重访周期短、空间分辨率高、数据获取方便、资料客观等优点。因此，迫切需要构建一种基于卫星遥感的综合干旱指数，并基于该指数提出多变量栅格化卫星遥感综合干旱风险评估方法。

综上分析可知，构建的多变量栅格化卫星遥感综合干旱风险评估方法在干旱监测、干旱风险评估方面具有一定的可靠性和综合性等优势，为地面观测站点少的区域提供了更加综合准确的评估方法，实现了由点到面对干旱灾害进行精细监测与风险评估，这对大规模干旱风险评估和生态保护具有指导意义，进而减少干旱灾害产生的巨大社会效益。

以上所述仅为本申请的部分优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (5)

1.一种多变量栅格化卫星遥感综合干旱风险评估方法，其特征在于：针对目标研究区域，执行以下步骤，实现对目标研究区域内各栅格的综合干旱事件进行综合干旱风险评估：

步骤1，基于目标研究区域的栅格划分，针对目标研究区域在预设历史时长段内的遥感数据，基于以该预设历史时长段内的每月为起点、预设个数月份为跨度的时间段，获得该预设历史时长段内各时间段的遥感数据，包含时间段内的卫星降水量、卫星增强型植被指数和卫星地表温度，进而得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；

所述步骤1的具体过程如下：

步骤1.1，获取目标研究区域各栅格在预设历史时长段内的遥感数据，包括月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度，并对目标研究区域的遥感数据进行插值处理统一分辨率；

步骤1.2，针对月尺度的卫星降水量、卫星增强型植被指数和卫星地表温度，基于以每月为起点，预设个数月份为跨度的时间段，获得预设历史时长段内各时间段的遥感数据，进而计算获得各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数、以及卫星地表温度的边际概率值；

步骤1.3，基于各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的边际概率值，计算获得各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值；

步骤1.4，基于各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值，计算获得各栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值；并基于各栅格在各时间段分别对应的累积概率值，得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；

步骤2，基于各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，分别针对各栅格在预设个数连续时间段的卫星遥感综合干旱指数CRSDI均小于预设卫星遥感综合干旱指数CRSDI阈值，表示该栅格发生了一次综合干旱事件，提取各栅格在预设历史时长段内分别对应各综合干旱事件的各预设综合干旱风险评估要素；

所述步骤2中基于预设卫星遥感综合干旱指数CRSDI＝-0.5为阈值，分别针对各栅格在预设个数连续时间段的卫星遥感综合干旱指数CRSDI均小于预设卫星遥感综合干旱指数CRSDI阈值，表示该栅格发生了一次综合干旱事件，提取各栅格在预设历史时长段内分别对应各综合干旱事件的各预设综合干旱风险评估要素；各预设综合干旱风险评估要素包括干旱历时、干旱烈度、干旱烈度峰值及干旱后生态恢复期；

干旱历时是指连续小于CRSDI阈值的时间段个数；

干旱烈度为连续小于CRSDI阈值的各时间段CRSDI累加和的绝对值；

干旱烈度峰值指在连续小于CRSDI阈值的各时间段中最小CRSDI的绝对值；

干旱后生态恢复期获取过程如下：

从综合干旱事件结束后的第一个月开始，逐月计算植被生长变化E，计算公式为：

式中，E_h为第h个月的植被生长变化，y_h为综合干旱事件结束后第h个月的卫星增强型植被指数，

为y_h所对应月份卫星增强型植被指数在预设历史时长段内的年平均值；E_h<0表明生态未恢复，干旱后生态恢复期为综合干旱事件结束后首次E_h≥0时h所对应的月份个数；

步骤3，针对各栅格对应各综合干旱事件的各预设综合干旱风险评估要素，逐栅格计算各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率；

所述步骤3中分别针对目标研究区域的各栅格，执行以下步骤，逐栅格计算各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率；

将该栅格在预设历史时长段内经历的各综合干旱事件的干旱历时由小到大的次序排列、干旱烈度由小到大的次序排列、干旱烈度峰值由小到大的次序排列、干旱后生态恢复期由小到大的次序排列；

该栅格在任一次综合干旱事件中干旱历时的经验频率D为：

其中，g为该栅格在预设历史时长段内共经历的综合干旱事件总次数，该次综合干旱事件的干旱历时排序为第J个；

该栅格在任一次综合干旱事件中干旱烈度的经验频率S为：

其中，该次综合干旱事件的干旱烈度排序为第K个；

该栅格在任一次综合干旱事件中干旱烈度峰值的经验频率P_E为：

其中，该次综合干旱事件的干旱烈度峰值排序为第H个；

该栅格在任一次综合干旱事件中干旱后生态恢复期的经验频率T为：

其中，该次综合干旱事件的干旱后生态恢复期排序为第L个；

步骤4，基于各栅格对应各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率，构建综合干旱风险评估指标CDRA，进而实现对目标研究区域内各栅格的综合干旱事件进行综合干旱风险评估；

所述步骤4中，基于各栅格对应各综合干旱事件中各预设综合干旱风险评估要素分别对应的经验频率，基于权重组合法得到各次综合干旱事件的风险评估指标CDRA，其计算公式为：

CDRA＝α×D+β×S+γ×P_E+δ×T

式中，α、β、γ和δ分别为D、S、P_E和T所对应的权重，α+β+γ+δ＝1。

2.根据权利要求1所述的一种多变量栅格化卫星遥感综合干旱风险评估方法，其特征在于：所述步骤1.2中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的边际概率值：

步骤1.2.1，针对该栅格在预设历史时长段内的月尺度的卫星降水量，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有f项，将该栅格的卫星降水量时间段序列按照各时间段结束月份划分为12个独立的子序列，结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星降水量均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量边际概率值F(x)：

其中，x指代任意一时间段的卫星降水量，各子序列的长度均为f_a，i指代x在对应的子序列中排序为第i个；

步骤1.2.2，对该栅格在预设历史时长段内的月尺度的卫星增强型植被指数，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有m项，将该栅格的卫星增强型植被指数时间段序列按照各时间段结束月份划分为12个独立的子序列，结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星增强型植被指数均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星增强型植被指数边际概率值F(y)：

其中，y指代任意一时间段的卫星增强型植被指数，各子序列的长度均为m_a，j指代y在对应的子序列中排序为第j个；

步骤1.2.3，对该栅格在预设历史时长段内的月尺度的卫星地表温度，基于每月为起点，预设个数月份为跨度的时间段，该栅格卫星降水量的时间段序列共有r项，将该栅格的卫星地表温度时间段序列按照各时间段结束月份划分为12个独立的子序列，结束月份相同的时间段在同一个子序列中，并分别针对各子序列，将子序列中各时间段的卫星增强型植被指数均按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段分别对应的卫星地表温度边际概率值F(z)：

其中，z指代任意一时间段的卫星增强型植被指数，各子序列的长度均为r_a，k指代z在对应的子序列中排序为第k个；

上述步骤中，卫星降水量时间段序列、卫星增强型植被指数时间段序列和卫星地表温度时间段序列长度一致，f＝m＝r，f_a＝m_a＝r_a。

3.根据权利要求2所述的一种多变量栅格化卫星遥感综合干旱风险评估方法，其特征在于：所述步骤1.3中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值，以及卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值：

步骤1.3.1，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量和卫星增强型植被指数的二元联合概率值F(x,y)：

F(x,y)＝P(X≤x,Y≤y)＝C(F(x),F(y))＝C(u,v)

式中，x指代任一时间段的卫星降水量；y指代任一时间段的卫星增强型植被指数；F(x)＝u为卫星降水量x所在时间段的边际概率值；F(y)＝v为卫星增强型植被指数y所在时间段的边际概率值；C为Copula函数；

其中，Copula函数为对称Archimedean Copulas家族中的Gumbel Copula函数，

式中，θ₁为参数，其取值范围为[1,∞)；exp为以自然常数e为底的指数函数，ln为以常数e为底数的对数函数；

步骤1.3.2，采用以下公式计算获得该栅格在各时间段分别对应的卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值F(x,y,z)：

F(x,y,z)＝P(X≤x,Y≤y,Z≤z)＝C(F(x),F(y),F(z))＝C(u,v,w)

式中，z指代任一时间段的卫星地表温度；F(z)＝w为卫星地表温度z所在时间段的边际概率值；

其中，Copula函数为非对称Archimedean Copulas中的Gumbel Copula函数，

式中，θ₂为参数，其取值范围为[1,∞)。

4.根据权利要求3所述的一种多变量栅格化卫星遥感综合干旱风险评估方法，其特征在于：所述步骤1.4中分别针对目标研究区域内的各栅格，执行以下步骤，获得各栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值，并基于各栅格在各时间段分别对应的累积概率值，得到各栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI；

步骤1.4.1，采用以下公式计算该栅格在各时间段的二元联合概率值与三元联合概率值之差，获得事件(X≤x,Y≤y,Z>z)的联合概率q：

P(X≤x,Y≤y,Z>z)＝P(X≤x,Y≤y)-P(X≤x,Y≤y,Z≤z)

＝C(F(x),F(y))-C(F(x),F(y),F(z))＝C(u,v)-C(u,v,w)＝q

步骤1.4.2，将各时间段的联合概率q按由小到大的次序排列，采用以下公式计算获得该栅格在各时间段的二元联合概率值与三元联合概率值之差的累积概率值F(q)：

其中，q指代任意一时间段的二元联合概率值与三元联合概率值之差，I指代q的排序为第I个，N指代时间段序列长度；

上述步骤中，二元联合概率值与三元联合概率值之差所得对应的时间段序列与卫星降水量时间段序列、卫星增强型植被指数时间段序列和卫星地表温度时间段序列长度一致，N＝f＝m＝r；

步骤1.4.3，基于该栅格在各时间段分别对应的二元联合概率值与三元联合概率值之差所得的累积概率值，通过标准正态化方法得到该栅格在各时间段分别对应的卫星遥感综合干旱指数CRSDI，其计算公式为：

式中，

为标准正态分布。

5.根据权利要求1所述的一种多变量栅格化卫星遥感综合干旱风险评估方法，其特征在于：所述步骤1.3中针对计算卫星降水量和卫星增强型植被指数的二元联合概率值，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值，分别针对各栅格，通过以下方法，分别选取计算卫星降水量和卫星增强型植被指数的二元联合概率值对应的最优Copula函数，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值对应的最优Copula函数：

基于Archimedean Copulas中的Gumbel Copula函数、Clayton Copula函数和FrankCopula函数，以及Meta-elliptic Copulas中的Student t Copula函数和Gaussian Copula函数，基于该5种Copula函数，通过以下公式，以赤池信息量准则AIC、贝叶斯信息准则BIC和方根误差准则RMSE三项中至少两项值最小作为最优的判断标准，从该5种Copula函数中分别选取计算计算卫星降水量和卫星增强型植被指数的二元联合概率值对应的最优Copula函数，以及计算卫星降水量、卫星增强型植被指数和卫星地表温度的三元联合概率值对应的最优Copula函数；

式中，MSE为均方误差，N为样本容量；b为Copula函数的参数个数；Pe_n为第n个时间段的联合分布经验频率，P_n为第n个时间段Copula函数计算的联合概率值；

第n个时间段的卫星降水量和卫星增强型植被指数的联合分布经验频率H(x_n,y_n)为

其中，x_n指代第n个时间段的卫星降水量，y_n指代第n个时间段的卫星增强型植被指数，N_b指代各时间段的样本中同时满足X≤x_n,Y≤y_n的样本个数；

第n个时间段的变量卫星降水量、卫星增强型植被指数和卫星地表温度的联合分布经验频率H(x_n,y_n,z_n)为：

其中，z_n指代第n个时间段的卫星地表温度，N_t指代各时间段的样本中同时满足X≤x_n,Y≤y_n,Z≤z_n的样本个数。

Images