一种分析旱涝急转事件特征量的方法

1.本发明涉及水文水资源应用领域，具体涉及一种分析旱涝急转事件特征量的方法。


背景技术：

2.旱涝急转是指某一区域或流域前期持续偏旱，后期发生集中强降雨，致使河流猛涨、城市内涝等状况，在短时间内由干旱状态转为洪涝状态，反映了旱涝极端事件在短期内的共存。在全球气候变暖、极端降水增多、水气循环增强的背景下，气候系统的稳定性降低。受气候变化和人类活动的影响，降水时空分布不均的情况进一步加剧，全球干旱、洪涝等极端灾害发生的频率和强度持续性增加，旱涝急转事件的频率和强度也在世界范围内持续增加，近年来已引起越来越多的国内外学者的广泛关注。
3.由于旱涝急转事件及其特征变化具有时空不连续性，对旱涝急转事件特征进行频率分析有助于制定防洪防旱措施、规划，为水资源时空分布不平衡的地区提供旱涝急转灾害预防的理论基础，在区域水资源规划与管理中具有重要意义。多研究表明旱涝灾害的强度、历时、频次等特征在非一致的气候背景(如全球变暖)驱动下具有明显变化趋势，因此，在对旱涝灾害进行频率分析时，应当对单变量时间序列或多变量相依性诊断其中是否存在变异点。在一致性情况下，估计变量的恒定参数；在非一致性情况下，估计变量的时变参数。但是，目前研究缺乏考虑旱涝急转事件的非一致性变化。在变化环境背景下，对检测到旱涝事件属性中存在的非一致性需经过讨论判别，在考虑潜在的非一致性下开展旱涝急转多变量频率分析，对气象灾害研究和防治具有十分重要的意义。


技术实现要素：

4.本发明考虑水文序列的时变统计特征，构造非一致背景下的旱涝急转事件特征量分析方法，对旱涝急转事件特征量进行分析，以此获取更加可靠的风险防范信息，及时对旱涝急转事件特征值风险作出评估，对旱涝急转灾害的分析与防治具有十分重要的意义。
5.本发明的技术方案如下：
6.一种分析旱涝急转事件特征量的方法，包括如下步骤：
7.s1：确定研究区域，收集流域内的日尺度降水量数据；
8.s2：基于流域的降水量数据，计算标准化加权平均降水swap指数，识别流域的旱涝急转过程，提取旱涝急转事件的烈度和历时特征；
9.s3：采用adf(单位根检验)分别检验烈度特征和历时特征的时间序列是否具有一致性；
10.当时间序列一致性成立时，采用参数恒定的概率分布模型分别对烈度特征和历时特征进行拟合，并基于sbc(schwarz bayesian information criterion)准则和aic(akaike information criterion)准则的评价，得到烈度边缘分布函数和历时边缘分布函数；
11.当时间序列一致性不成立时，采用广义累加模型(gamlss)分别对烈度特征和历时特征进行拟合，得到烈度边缘分布函数和历时边缘分布函数；
12.s4：采用copula函数似然比方法(copula-based likelihood-ratio test，clr test)检验烈度边缘分布函数和历时边缘分布函数之间的一致性假设是否成立；
13.当一致性假设成立时，采用参数恒定的copula联合分布函数模拟得到历时和烈度的联合分布函数；
14.当一致性假设不成立时，采用参数时变的copula联合分布函数模拟得到历时和烈度的联合分布函数；
15.s5：基于历时和烈度的联合分布函数，计算历时和烈度“或”与“且”的联合重现期，并采用最可能组合法(most-likely weight function)推求得到历时和烈度的联合设计值，将该设计值作为旱涝急转事件灾害防治的风险阈值。
16.进一步，步骤s2中标准化加权平均降水swap指数的计算过程如下：
17.swap指数是以当日旱涝状态受前期旱涝状态和当日降水的影响为前提的指数，其计算方法是假定同一日加权平均降水(wap)序列服从gamma分布，通过γ分布对加权平均降水序列进行拟合，最后对累积概率分布进行标准正态化处理，加权平均降水序列的计算公式如下所示：
[0018][0019]
wn＝(1-a)a
n-1
[0020]
式中，wapi表示加权平均降水序列，pn表示第n天的降水，wn表示pn的权重，a是权重随时间衰减的参数，an是第n天权重随时间的衰减参数，n是前期降水影响天数，a与n的选取根据研究区当地土壤-水系统而定，取a＝0.9，n＝45。
[0021]
进一步，步骤s2中，基于计算的swap指数，采用游程理论识别提取流域的旱涝急转过程，具体步骤如下：
[0022]
对某一旱涝指数时间序列，设定指数临界值x0、x1和时间长度临界值d0、d1和d2，当指数持续低于x0的时间不小于d0时，认为出现干旱；当指数持续高于x1的时间不小于d1时，认为出现洪涝；当干旱与洪涝之间相隔的时间不大于d2时，认为构成了旱涝急转事件；根据旱涝指数划分标准，选用其中旱等级上界与中涝等级下界对应的数值作为x0和x1，即x0＝-1和x1＝1；参考各类急转事件识别研究的游程理论识别过程，确定基于日值指数的识别时长阈值d0＝d1＝10天，干旱状态与洪涝状态之间的急转期时长阈值d2＝8天，根据提取的旱涝急转事件，得到每次旱涝急转事件发生的历时和烈度。
[0023]
进一步，步骤s3中，当时间序列一致性成立时，拟合得到烈度边缘分布函数和历时边缘分布函数的具体过程如下；
[0024]
选取在气象水文研究中较常用的概率分布模型作为备选模型，具体包括指数(exponential,exp)、伽马(gamma)、三参数的广义伽马(generalized gamma,gamma3)、正态(normal)、对数正态(log-normal,lnorm)、逻辑斯蒂(logistic,logis)、韦布尔(weibull)一共7个概率分布模型，通过这7个备选模型分别对烈度特征和历时特征进行拟合，并通过sbc准则对上述7个备选模型的拟合优度进行评价；
[0025]
sbc的计算公式如下：
[0026]
sbc＝gd+log(ξ)
×
df＝-2 log l(θ)+log(ξ)
×
df
[0027]
式中，gd为全局拟合偏差，df是自由度数目，θ和ξ分别代表概率分布参数的估计值和参数个数，l(
·
)代表对数似然函数值；
[0028]
为避免造成分布的过拟合，同时也考虑通过aic准则对上述7种备选模型的拟合优度进行评价；
[0029]
aic的计算公式如下：
[0030]
aic＝2ξ-2 log l(θ)
[0031]
最后，在所有通过检验的备选模型中，选出使得aic和sbc的计算值同时达到最小的分布模型为相对最优模型，并将该相对最优模型作为烈度边缘分布函数或历时边缘分布函数。
[0032]
进一步，对于备选的正态和对数正态概率分布模型，其对变量的定义区间为[-∞,+∞]，而待检验的烈度特征和历时特征均大于0，因此需将上述2个备选模型的定义区间修正为正值区间，用以下公式表示修正后的累积概率：
[0033][0034]
式中，f(x|θ)为概率密度函数，θ为概率分布参数。
[0035]
进一步，步骤s3中，当时间序列一致性不成立时，拟合得到烈度边缘分布函数和历时边缘分布函数的具体过程如下；
[0036]
采用可综合考虑位置、尺度和形状参数的广义累加模型(generalized additive models for location,scale and shape，gamlss)进行边缘分布拟合，在广义累加模型中，假设响应变量y的个独立观测值服从概率分布f(yi|θi)；
[0037]
其中，θi＝(μi,σi,vi,τi)是概率分布参数组成的向量，μ和σ为位置和尺度参数，v和τ为形状参数，分别表征分布的歪度和峭度；
[0038]
通过连接函数建立概率分布参数与其解释变量的联系，表达式为：
[0039][0040]
即：
[0041][0042]
式中，gk(
·
)代表连接函数；θk、μ、σ、v、τ是长度为的列向量；
[0043]
xk是已知矩阵，对应个时段的解释变量值；k表示待分析序列的序号；是待求解的参数向量，维度为jk；z
jk
是已知设计矩阵；γ
jk
是由随机变量组成的列向量。
[0044]
在本发明中，将识别出的旱涝急转发生时间(drought onset)作为解释变量，建立急转历时与烈度概率分布参数与旱涝急转发生时间的函数关系，并采用gamlss模型的半参数累加形式，将概率分布参数表达为旱涝急转起始时间的三次样条函数，有利于捕捉可能存在的非线性关系。本发明以最大化似然函数为目标，计算备选概率分布的时变参数，并依据sbc准则和aic准则从7种备选分布模型中筛选最优的烈度边缘概率函数和历时边缘概率函数。
[0045]
进一步，步骤s4中，采用copula函数似然比方法检验烈度边缘分布函数和历时边缘分布函数之间的一致性假设是否成立的过程如下；
[0046]
copula函数似然比方法的基本假设是，固定copula函数类型不变，通过检测其参数θc的变化来分析多变量水文序列相依结构强度的变点；
[0047]
已知y1,y2,...,yd是d维随机变量观测值组成的时间序列，其中yi＝(y
1,i
,y
2,i
,...,y
n,i
),i＝1,2,...,d，观测值所在总体服从概率分布f(y1|θ
c,1
),f(y2|θ
c,2
),...,f(yn|θ
c,n
)，θ
c,i
是不同时段对应的copula函数参数，在copula函数似然比检验中，原假设是多变量的相依结构不存在变点(即所有观测值服从同一概率分布)，表达公式如下：
[0048]
h0:θ
c,1
＝θ
c,2
＝...＝θ
c,n
＝η0[0049]
备择假设h1为：存在k
*
且使得θ
c,1
＝θ
c,2
＝...＝θ
c,k*
＝η1，并且θ
c,k*
＝...＝θ
c,n
＝η2，同时η1≠η2；构造基于copula函数的似然比λk，当λk小于给定阈值时，拒绝原假设h0，并认为多变量相依结构存在变点k
*
；
[0050][0051]
式中，ln(
·
)代表整个样本序列的似然函数；lk(
·
)和分别是检验点k之前和之后的似然函数；ui表示双变量序列的边缘分布概率；c(
·
)表示copula概率密度函数；和分别是copula函数参数η0、η1和η2的极大似然估计值，改写λk为对数形式，得到：
[0052][0053]
进一步，构造对数形式的统计检验量如下：
[0054][0055]
当zn大于给定显著性水平对应的阈值时，拒绝原假设h0，似然比统计量服从渐近分布：
[0056][0057]
式中，z表示给定的阈值，当z趋近于无穷大时，h＝l＝[ln(n)]
3/2
；
[0058]
h、l、o表示渐近分布的参数；p是存在一个变点的copula函数的参数数目；将代入上式，如果计算得到p值小于阈值，则在该阈值显著性水平下拒绝原假设，一致性假设不成立，反之则成立。
[0059]
进一步，步骤s4中，根据copula函数建立的历时和烈度的联合分布函数表达为：
[0060]
p(d,s)＝c(ud,us|θc)
[0061]
式中，d、s分别表示烈度和历时；ud,us表示烈度和历时的边缘分布；
[0062]
与gamlss模型类似，为了估算copula的时变参数，引入解释变量，用于描述copula时变参数与解释变量之间的函数关系；
[0063]
g(θc)＝β0+xβ1[0064]
式中，g(θc)表示当copula参数是θc的连接函数；
[0065]
θc＝(θ
c,1
,θ
c,2
,...,θ
c,n
)
t
是copula时变参数组成的列向量；
[0066]
x＝(x1,x2,...,xn)
t
代表解释变量的n个观测值；β0和β表示待求解的参数向量；
[0067]
下面根据历时和烈度关系的一致性假设成立与否，分情况给出连接函数的表达式：
[0068]
(1)当两变量关系保持一致性时，copula参数是不随时间变化的恒定值，连接函数表示成：
[0069]
g(θc)＝constant
[0070]
(2)当一致性假设破坏时，copula参数是随时间变化的函数，因而引入干旱事件的发生时间time，作为解释变量，连接函数改写为：
[0071]
g(θc)＝β0+time
·
β1[0072]
进一步，采用边际函数推断法，通过最大化对数似然函数，求解copula参数θ
c,t
，公式如下：
[0073]
[0074]
式中，u
d,i
表示第i时段的烈度边缘分布；u
s,i
表示第i时段的历时边缘分布；di表示第i时段的烈度；si表示第i时段的历时；θd表示烈度边缘分布的参数；θs表示历时边缘分布的参数；θ
d,i
表示第i时段的烈度边缘分布的参数；θ
s,i
表示第i时段的历时边缘分布的参数；c(
·
)是copula的概率密度函数；fd(
·
)和fs(
·
)分别代表历时、烈度边缘分布的概率密度函数；在边际函数推断法(ifm)中，由于先前已经通过拟合历时和烈度的边缘分布得到上述公式等号右侧第二、三项的最大值，所以只需最大化等号右边第一项，即可求得最优的β0和β1，并将β0和β1的估计值代入计算copula函数的估计值；
[0075][0076]
式中，timei表示第i时段的发生时间。
[0077]
本发明根据描述极限条件下变量之间相依结构的不同，常用的copula函数可以分为椭圆形(elliptical)copula、阿基米德(archimedean)copula、极值(extreme value)copula以及混合型copula等，目前在水文领域应用较为广泛的为椭圆形copula和阿基米德copula。本发明共采用了一种椭圆形copula(gaussian copula)以及四种阿基米德copula函数(clayton copula、frank copula、gumbel copula和joe copula)，共五种copula函数作为备选单参数的备选copula函数。拟合后需对选取的copula函数进行拟合优度检验。依照贝叶斯信息准则(bayesian information criterion,bic)，选取bic值最小的函数作为最优联合分布copula函数，建立历时和烈度的联合分布函数。
[0078]
进一步，步骤s5中，计算历时和烈度“或”与“且”的联合重现期的过程如下：
[0079]
旱涝急转的单变量特征p(x≥x
*
)的重现期计算如下：
[0080][0081]
其中，t为旱涝急转事件的重现期，f(x
*
)为x
*
事件概率函数；interval为所有旱涝急转事件平均发生的时间间隔，单位为年，即某一次旱涝急转事件开始到下一次旱涝急转事件开始之间的时间长度；若在时长为n
year
年内取样得到n
alter
场急转事件，interval计算公式为：
[0082][0083]
多变量的旱涝急转频率分析具有多种组合的联合重现期，包括“或”联合重现期和“且”联合重新期，“或”联合重现期是指旱涝急转事件随机变量有一个或以上超过给定的阈值计算得到的重现期，用or操作符“∨”表示；“且”联合重现期是指旱涝急转事件随机变量均超过给定的阈值计算得到的重现期，因此“且”重现期又被称为“联合超越重现期”或“同现重现期”，用and操作符“∧”表示；
[0084]
不同视角下的重现期计算公式如下；
[0085]
在单变量视角下：
[0086]
[0087][0088]
式中，ts、td分别为历时、烈度的单变量视角下重现期，fs(s)和fd(d)分别表示历时和烈度发生的概率；
[0089]
在双变量视角下：
[0090][0091][0092]
式中，表示在历时和烈度两个变量视角下的“或”联合重现期，表示在历时和烈度两个变量视角下的“且”联合重现期，f
s,d
(s,d)为历时与烈度的联合发生概率；
[0093]
基于上述各式计算得到历时和烈度“或”与“且”的联合重现期。
[0094]
进一步，步骤s5中，采用最可能组合法推求得到历时和烈度的联合设计值的过程如下：
[0095]
旱涝急转事件的联合设计值是指在某一设计重现期水平下，由急转特征组成的急转联合密度f(u,v)达到最大值时的多变量组合(u
ml
,v
ml
)，即该组合出现的可能性最大，因此，通过构造以下方程求解：
[0096]
[u
ml
,v
ml
]＝argmaxf(u,v)
[0097]
f(u,v)＝c(u,v)
·
f(u)
·
f(v)
[0098]
其中，f(u,v)为急转双变量的联合概率密度函数，c(u,v)为对应联合copula分布的概率密度函数，f(u)、f(v)为变量的边缘密度函数，[u
ml
,v
ml
]为所求的多变量的联合设计值组合；
[0099]
采用最可能组合法推求得到联合设计值的计算步骤如下：
[0100]
(1)、基于copula联合分布函数，采用monte carlo法模拟多组变量组合(u
ml
,v
ml
)，此取组数取值为100000；
[0101]
(2)、根据“或”、“且”重现期的计算公式，求出上一步的双变量组合结合对应的联合重现期；
[0102]
(3)、根据确定的重现期水平，筛选满足精度要求(e＝10-4
)的变量组合，其中联合概率密度值最大的变量组合即为所推求的联合设计值。
[0103]
相比现有技术，本发明具有的有益效果如下：
[0104]
本发明提供一种基于非一致性理论分析旱涝急转事件特征量的方法，通过运用时变参数copula函数建立非一致性背景下旱涝急转特征相依性较强的历时-烈度之间的联合分布，并计算各重现期下的旱涝急转特征的联合设计值，以此联合设计值作为对旱涝急转事件的风险阈值，可对旱涝急转事件特征值风险作出更进一步地评估，对气象灾害研究和防治具有十分重要的意义。
附图说明
[0105]
图1为本发明的流程示意图。
具体实施方式
[0106]
附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本专利的限制。
[0107]
实施例1：
[0108]
如图1所示，一种分析旱涝急转事件特征量的方法，包括如下步骤：
[0109]
s1：确定研究区域，收集流域内的日尺度降水量数据；
[0110]
s2：基于流域的降水量数据，计算标准化加权平均降水swap指数，识别流域的旱涝急转过程，提取旱涝急转事件的烈度和历时特征；
[0111]
s3：采用adf(单位根检验)分别检验烈度特征和历时特征的时间序列是否具有一致性；
[0112]
当时间序列一致性成立时，采用参数恒定的概率分布模型分别对烈度特征和历时特征进行拟合，并基于sbc(schwarz bayesian information criterion)准则和aic(akaike information criterion)准则的评价，得到烈度边缘分布函数和历时边缘分布函数；
[0113]
当时间序列一致性不成立时，采用广义累加模型(gamlss)分别对烈度特征和历时特征进行拟合，得到烈度边缘分布函数和历时边缘分布函数；
[0114]
s4：采用copula函数似然比方法(copula-based likelihood-ratio test，clr test)检验烈度边缘分布函数和历时边缘分布函数之间的一致性假设是否成立；
[0115]
当一致性假设成立时，采用参数恒定的copula联合分布函数模拟得到历时和烈度的联合分布函数；
[0116]
当一致性假设不成立时，采用参数时变的copula联合分布函数模拟得到历时和烈度的联合分布函数；
[0117]
s5：基于历时和烈度的联合分布函数，计算历时和烈度“或”与“且”的联合重现期，并采用最可能组合法(most-likely weight function)推求得到历时和烈度的联合设计值，将该设计值作为旱涝急转事件灾害防治的风险阈值。
[0118]
本实施例所使用的气象资料可来源于国家气象科学数据中心产品“中国地面气候资料日值数据集(v3.0)”，所选用的流域可为珠江流域，收集1961-2020年珠江流域内75个气象站点的降水量的逐日数据。
[0119]
在本实施中，步骤s2中标准化加权平均降水swap指数的计算过程如下：
[0120]
swap指数是以当日旱涝状态受前期旱涝状态和当日降水的影响为前提的指数，其计算方法是假定同一日加权平均降水(wap)序列服从gamma分布，通过γ分布对加权平均降水序列进行拟合，最后对累积概率分布进行标准正态化处理，加权平均降水序列的计算公式如下所示：
[0121][0122]
wn＝(1-a)a
n-1
[0123]
式中，wapi表示加权平均降水序列，pn表示第n天的降水，wn表示pn的权重，a是权重随时间衰减的参数，an是第n天权重随时间的衰减参数，n是前期降水影响天数，a与n的选取根据研究区当地土壤-水系统而定，取a＝0.9，n＝45。
[0124]
在本实施中，步骤s2中，基于计算的swap指数，采用游程理论识别提取流域的旱涝急转过程，具体步骤如下：
[0125]
对某一旱涝指数时间序列，设定指数临界值x0、x1和时间长度临界值d0、d1和d2，当指数持续低于x0的时间不小于d0时，认为出现干旱；当指数持续高于x1的时间不小于d1时，认为出现洪涝；当干旱与洪涝之间相隔的时间不大于d2时，认为构成了旱涝急转事件；根据旱涝指数划分标准，选用其中旱等级上界与中涝等级下界对应的数值作为x0和x1，即x0＝-1和x1＝1；参考各类急转事件识别研究的游程理论识别过程，确定基于日值指数的识别时长阈值d0＝d1＝10天，干旱状态与洪涝状态之间的急转期时长阈值d2＝8天，根据提取的旱涝急转事件，得到每次旱涝急转事件发生的历时和烈度。
[0126]
在本实施中，步骤s3中，当时间序列一致性成立时，拟合得到烈度边缘分布函数和历时边缘分布函数的具体过程如下；
[0127]
选取在气象水文研究中较常用的概率分布模型作为备选模型，具体包括指数(exponential,exp)、伽马(gamma)、三参数的广义伽马(generalized gamma,gamma3)、正态(normal)、对数正态(log-normal,lnorm)、逻辑斯蒂(logistic,logis)、韦布尔(weibull)一共7个概率分布模型，通过这7个备选模型分别对烈度特征和历时特征进行拟合，并通过sbc准则对上述7个备选模型的拟合优度进行评价；
[0128]
sbc的计算公式如下：
[0129]
sbc＝gd+log(ξ)
×
df＝-2 log l(θ)+log(ξ)
×
df
[0130]
式中，gd为全局拟合偏差，df是自由度数目，θ和ξ分别代表概率分布参数的估计值和参数个数，l(
·
)代表对数似然函数值；
[0131]
为避免造成分布的过拟合，同时也考虑通过aic准则对上述7种备选模型的拟合优度进行评价；
[0132]
aic的计算公式如下：
[0133]
aic＝2ξ-2logl(θ)
[0134]
最后，在所有通过检验的备选模型中，选出使得aic和sbc的计算值同时达到最小的分布模型为相对最优模型，并将该相对最优模型作为烈度边缘分布函数或历时边缘分布函数。
[0135]
在本实施中，对于备选的正态和对数正态概率分布模型，其对变量的定义区间为[-∞,+∞]，而待检验的烈度特征和历时特征均大于0，因此需将上述2个备选模型的定义区间修正为正值区间，用以下公式表示修正后的累积概率：
[0136][0137]
式中，f(x|θ)为概率密度函数，θ为概率分布参数。
[0138]
在本实施中，步骤s3中，当时间序列一致性不成立时，拟合得到烈度边缘分布函数和历时边缘分布函数的具体过程如下；
[0139]
采用可综合考虑位置、尺度和形状参数的广义累加模型(generalized additive models for location,scale and shape，gamlss)进行边缘分布拟合，在广义累加模型中，假设响应变量y的个独立观测值服从概率分布f(yi|θi)；
[0140]
其中，θi＝(μi,σi,vi,τi)是概率分布参数组成的向量，μ和σ为位置和尺度参数，ν和τ为形状参数，分别表征分布的歪度和峭度；
[0141]
通过连接函数建立概率分布参数与其解释变量的联系，表达式为：
[0142][0143]
即：
[0144][0145]
式中，gk(
·
)代表连接函数；θk、μ、σ、ν、τ是长度为的列向量；
[0146]
xk是已知矩阵，对应个时段的解释变量值；k表示待分析序列的序号；是待求解的参数向量，维度为jk；z
jk
是已知设计矩阵；γ
jk
是由随机变量组成的列向量。
[0147]
在本实施中中，将识别出的旱涝急转发生时间(drought onset)作为解释变量，建立急转历时与烈度概率分布参数与旱涝急转发生时间的函数关系，并采用gamlss模型的半参数累加形式，将概率分布参数表达为旱涝急转起始时间的三次样条函数，有利于捕捉可能存在的非线性关系。本发明以最大化似然函数为目标，计算备选概率分布的时变参数，并依据sbc准则和aic准则从7种备选分布模型中筛选最优的烈度边缘概率函数和历时边缘概率函数。
[0148]
在本实施中，步骤s4中，采用copula函数似然比方法检验烈度边缘分布函数和历时边缘分布函数之间的一致性假设是否成立的过程如下；
[0149]
copula函数似然比方法的基本假设是，固定copula函数类型不变，通过检测其参数θc的变化来分析多变量水文序列相依结构强度的变点；
[0150]
已知y1,y2,...,yd是d维随机变量观测值组成的时间序列，其中yi＝(y
1,i
,y
2,i
,...,y
n,i
),i＝1,2,...,d，观测值所在总体服从概率分布f(y1|θ
c,1
),f(y2|θ
c,2
),...,f(yn|θ
c,n
)，θ
c,i
是不同时段对应的copula函数参数，在copula函数似然比检验中，原假设是多变量的相依结构不存在变点(即所有观测值服从同一概率分布)，表达公式如下：
[0151]
h0:θ
c,1
＝θ
c,2
＝...＝θ
c,n
＝η0[0152]
备择假设h1为：存在k
*
且使得θ
c,1
＝θ
c,2
＝...＝θ
c,k*
＝η1，并且θ
c,k*
＝...＝θ
c,n
＝η2，同时η1≠η2；构造基于copula函数的似然比λk，当λk小于给定阈值时，拒绝原假设h0，并认为多变量相依结构存在变点k
*
；
[0153][0154]
式中，ln(
·
)代表整个样本序列的似然函数；lk(
·
)和分别是检验点k之前和之后的似然函数；ui表示双变量序列的边缘分布概率；c(
·
)表示copula概率密度函数；和分别是copula函数参数η0、η1和η2的极大似然估计值，改写λk为对数形式，得到：
[0155][0156]
进一步，构造对数形式的统计检验量如下：
[0157][0158]
当zn大于给定显著性水平对应的阈值时，拒绝原假设h0，似然比统计量服从渐近分布：
[0159][0160]
式中，z表示给定的阈值，当z趋近于无穷大时，h＝l＝[ln(n)]
3/2
；
[0161]
h、l、o表示渐近分布的参数；p是存在一个变点的copula函数的参数数目；将代入上式，如果计算得到p值小于阈值，则在该阈值显著性水平下拒绝原假设，一致性假设不成立，反之则成立。
[0162]
在本实施中，步骤s4中，根据copula函数建立的历时和烈度的联合分布函数表达为：
[0163]
p(d,s)＝c(ud,us|θc)
[0164]
式中，d、s分别表示烈度和历时；ud,us表示烈度和历时的边缘分布；
[0165]
与gamlss模型类似，为了估算copula的时变参数，引入解释变量，用于描述copula时变参数与解释变量之间的函数关系；
[0166]
g(θc)＝β0+xβ1[0167]
式中，g(θc)表示当copula参数是θc的连接函数；
[0168]
θc＝(θ
c,1
,θ
c,2
,...,θ
c,n
)
t
是copula时变参数组成的列向量；
[0169]
x＝(x1,x2,...,xn)
t
代表解释变量的n个观测值；β0和β表示待求解的参数向量；
[0170]
下面根据历时和烈度关系的一致性假设成立与否，分情况给出连接函数的表达式：
[0171]
(1)当两变量关系保持一致性时，copula参数是不随时间变化的恒定值，连接函数表示成：
[0172]
g(θc)＝constant
[0173]
(2)当一致性假设破坏时，copula参数是随时间变化的函数，因而引入干旱事件的发生时间time，作为解释变量，连接函数改写为：
[0174]
g(θc)＝β0+time
·
β1[0175]
进一步，采用边际函数推断法，通过最大化对数似然函数，求解copula参数θ
c,t
，公式如下：
[0176][0177]
式中，u
d,i
表示第i时段的烈度边缘分布；u
s,i
表示第i时段的历时边缘分布；di表示第i时段的烈度；si表示第i时段的历时；θd表示烈度边缘分布的参数；θs表示历时边缘分布的参数；θ
d,i
表示第i时段的烈度边缘分布的参数；θ
s,i
表示第i时段的历时边缘分布的参数；c(
·
)是copula的概率密度函数；fd(
·
)和fs(
·
)分别代表历时、烈度边缘分布的概率密度函数；在边际函数推断法(ifm)中，由于先前已经通过拟合历时和烈度的边缘分布得到上述公式等号右侧第二、三项的最大值，所以只需最大化等号右边第一项，即可求得最优的β0和β1，并将β0和β1的估计值代入计算copula函数的估计值；
[0178][0179]
式中，timei表示第i时段的发生时间。
[0180]
本发明根据描述极限条件下变量之间相依结构的不同，常用的copula函数可以分为椭圆形(elliptical)copula、阿基米德(archimedean)copula、极值(extreme value)copula以及混合型copula等，目前在水文领域应用较为广泛的为椭圆形copula和阿基米德copula。本发明共采用了一种椭圆形copula(gaussian copula)以及四种阿基米德copula函数(clayton copula、frank copula、gumbel copula和joe copula)，共五种copula函数作为备选单参数的备选copula函数。拟合后需对选取的copula函数进行拟合优度检验。依照贝叶斯信息准则(bayesian information criterion,bic)，选取bic值最小的函数作为最优联合分布copula函数，建立历时和烈度的联合分布函数。
[0181]
在本实施中，步骤s5中，计算历时和烈度“或”与“且”的联合重现期的过程如下：
[0182]
旱涝急转的单变量特征p(x≥x
*
)的重现期计算如下：
[0183][0184]
其中，t为旱涝急转事件的重现期，f(x
*
)为x
*
事件概率函数；interval为所有旱涝急转事件平均发生的时间间隔，单位为年，即某一次旱涝急转事件开始到下一次旱涝急转事件开始之间的时间长度；若在时长为n
year
年内取样得到n
alter
场急转事件，interval计算公式为：
[0185][0186]
多变量的旱涝急转频率分析具有多种组合的联合重现期，包括“或”联合重现期和“且”联合重新期，“或”联合重现期是指旱涝急转事件随机变量有一个或以上超过给定的阈值计算得到的重现期，用or操作符“∨”表示；“且”联合重现期是指旱涝急转事件随机变量均超过给定的阈值计算得到的重现期，因此“且”重现期又被称为“联合超越重现期”或“同现重现期”，用and操作符“∧”表示；
[0187]
不同视角下的重现期计算公式如下；
[0188]
在单变量视角下：
[0189][0190][0191]
式中，ts、td分别为历时、烈度的单变量视角下重现期，fs(s)和fd(d)分别表示历时和烈度发生的概率；
[0192]
在双变量视角下：
[0193][0194][0195]
式中，表示在历时和烈度两个变量视角下的“或”联合重现期，表示在历时和烈度两个变量视角下的“且”联合重现期，f
s,d
(s,d)为历时与烈度的联合发生概率；
[0196]
基于上述各式计算得到历时和烈度“或”与“且”的联合重现期。
[0197]
在本实施中，步骤s5中，采用最可能组合法推求得到历时和烈度的联合设计值的过程如下：
[0198]
旱涝急转事件的联合设计值是指在某一设计重现期水平下，由急转特征组成的急转联合密度f(u,v)达到最大值时的多变量组合(u
ml
,v
ml
)，即该组合出现的可能性最大，因此，通过构造以下方程求解：
[0199]
[u
ml
,v
ml
]＝argmax(u,v)
[0200]
f(u,v)＝c(u,v)
·
f(u)
·
f(v)
[0201]
其中，f(u,v)为急转双变量的联合概率密度函数，c(u,v)为对应联合copula分布的概率密度函数，f(u)、f(v)为变量的边缘密度函数，[u
ml
,v
ml
]为所求的多变量的联合设计值组合；
[0202]
采用最可能组合法推求得到联合设计值的计算步骤如下：
[0203]
(1)、基于copula联合分布函数，采用monte carlo法模拟多组变量组合(u
ml
,v
ml
)，此取组数取值为100000；
[0204]
(2)、根据“或”、“且”重现期的计算公式，求出上一步的双变量组合结合对应的联合重现期；
[0205]
(3)、根据确定的重现期水平，筛选满足精度要求(e＝10-4
)的变量组合，其中联合概率密度值最大的变量组合即为所推求的联合设计值。
[0206]
本发明提供一种基于非一致性理论分析旱涝急转事件特征量的方法，通过运用时变参数copula函数建立非一致性背景下旱涝急转特征相依性较强的历时-烈度之间的联合分布，并计算各重现期下的旱涝急转特征的联合设计值，以此联合设计值作为对旱涝急转事件的风险阈值，可对旱涝急转事件特征值风险作出更进一步地评估，对气象灾害研究和防治具有十分重要的意义。
[0207]
实施例2：
[0208]
本实施例提供一种分析旱涝急转事件特征量的装置，包括依次通讯连接的采集模块、提取特征模块、筛选模块、检验模块及输出模块；
[0209]
其中，采集模块：用于收集流域内的日尺度降水量数据；
[0210]
提取特征模块：基于采集模块的数据，计算标准化加权平均降水swap指数，识别流域的旱涝急转过程，提取旱涝急转事件的烈度和历时特征；
[0211]
筛选模块：用于根据提取特征模块得到的烈度特征和历时特征，采用adf(单位根检验)分别检验烈度特征和历时特征的时间序列是否具有一致性；
[0212]
当时间序列一致性成立时，采用参数恒定的概率分布模型分别对烈度特征和历时特征进行拟合，并基于sbc(schwarz bayesian information criterion)准则和aic(akaike information criterion)准则的评价，得到烈度边缘分布函数和历时边缘分布函数；
[0213]
当时间序列一致性不成立时，采用广义累加模型(gamlss)分别对烈度特征和历时特征进行拟合，得到烈度边缘分布函数和历时边缘分布函数；
[0214]
检验模块：用于根据筛选模块得到的烈度边缘分布函数和历时边缘分布函数，采用copula函数似然比方法(copula-based likelihood-ratio test，clr test)检验烈度边缘分布函数和历时边缘分布函数之间的一致性假设是否成立；
[0215]
当一致性假设成立时，采用参数恒定的copula联合分布函数模拟得到历时和烈度的联合分布函数；
[0216]
当一致性假设不成立时，采用参数时变的copula联合分布函数模拟得到历时和烈度的联合分布函数；
[0217]
输出模块：基于检验模块得到的历时和烈度的联合分布函数，计算历时和烈度“或”与“且”的联合重现期，并采用最可能组合法(most-likely weight function)推求得到历时和烈度的联合设计值，最后将该设计值作为旱涝急转事件灾害防治的风险阈值进行输出。
[0218]
其中，上述的模块可集成在一个系统中进行实现，最后通过输出口将结果输出到显示屏上进行显示。
[0219]
实施例3：
[0220]
本实施例提供一种电子设备，包括处理器、通信接口、存储器和通信总线，其中，处理器，通信接口，存储器通过通信总线完成相互间的通信；存储器，用于存放计算机程序；处理器，用于执行存储器上所存放的程序，实现上述实施例1中的分析旱涝急转事件特征量的方法。
[0221]
本实施例还提供一种计算机可读存储介质，该计算机可读存储介质内存储有计算机程序，计算机程序被处理器执行时实现上述的分析旱涝急转事件特征量的方法。
[0222]
显然，本发明的上述实施例仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。