1.本发明涉及一种趋势检测方法，尤其涉及一种水文气象极值序列的趋势检测方法。
背景技术：2.全球变暖是对地球的主要威胁，气候变暖是人为温室气体排放增加的结果，正在逐步改变地球气候模式。根据政府间气候变化专门委员会ipcc的说法，温度升高导致的水文水循环特征急剧变化。因此，全球变暖将影响水文过程，并导致世界不同地区出现水文气象极端事件的风险增加，为了有效预测和防范水文气象极端事件发生风险，最为关键一步就是对水文气象极端事件进行有效的趋势性检测和识别。目前国内外使用最多的是非参数法的趋势性检测方法：
3.mann-kendall趋势性检测方法目前应用于水文气象极值序列趋势性检测最为广泛，由于其仅仅适用于单调性趋势的检测，对于非线性，如先增加后降低的非单调性趋势趋势的检测效果存在明显不足，且容易出现假设检验中的两类错误：“弃真”和“纳伪”错误。
技术实现要素：4.发明目的：本发明旨在提供一种水文气象极值序列的趋势检测方法，解决水文气象极值序列趋势性检测只能单调性趋势监测的问题。。
5.技术方案：本发明的一种参数化的水文气象极值序列的趋势检测方法,包括以下步骤：
6.(1)采用广义极值分布gev函数建立水文气象极值序列的趋势模型；
7.(2)基于贝叶斯推断框架估计水文气象极值序列的趋势模型的参数；
8.(3)根据趋势模型的参数估计值获得贝叶斯对数似然比，然后得出水文气象极值序列的趋势。
9.步骤(1)的趋势模型：
[0010][0011]
式中，x
t
为水文气象数据，gev分布的位置μ
t
，尺度σ
t
和形状参数k
t
是随着时间t逐年变化的。
[0012]
水文气象极值序列的趋势模型gev分布的位置μ
t
的趋势结构包括：
[0013][0014]
式中，μ0、μ1、μ2、μ3为未知待定参数。
[0015]
水文气象极值序列的趋势模型gev分布的尺度σ
t
的趋势结构包括：
[0016][0017]
式中，σ0、σ1、σ2、σ3为未知待定参数。
[0018]
水文气象极值序列的趋势模型gev分布的形状参数κ
t
的趋势结构包括：
[0019][0020]
式中，κ0、κ1、κ2、κ3为未知待定参数。
[0021]
步骤(2)包括以下步骤：
[0022]
(21)趋势模型的参数先验分布的假定；
[0023]
(22)基于gibbs的蒙特卡洛随机采样法获得潜在的后验分布参数；
[0024]
(23)获得水文气象极值序列的趋势模型参数估计值。
[0025]
步骤(2)中假设先验分布参数集为β,后验分布按照公式(5)求得：
[0026][0027]
式中，f(x
t
|β,t)各时间步时变gev参数的后验分布概率密度值。
[0028]
步骤(3)中获得贝叶斯对数似然比
[0029]
bay_lr＝2(logl
1-logl0)
ꢀꢀ6[0030]
式中，logl0为平稳模型的对数似然值，而logl1代表了其余趋势性模型的对数似然值。
[0031]
步骤(3)中bay_lr的值大于的α分位数，则水文气象极值序列存在趋势性；当α＝0.95时，获得95％显著水平的水文气象极值序列趋势性。
[0032]
有益效果：与现有技术相比，本发明具有如下显著优点：
[0033]
(1)有效克服了非参数趋势检测法对于水文气象极值序列非单调性趋势检测识别效果不佳的问题，从参数化角度有效保证水文气象极值序列的趋势检测精度。
[0034]
(2)基于贝叶斯统计推断框架的使用一方面充分保证复杂趋势模式情形下对于水文气象极值序列的分布参数估计精度不会因为引入过多的参数而受到影响，另一方面采用gibbs采样算法简化了水文气象极值序列分布参数的估计过程。
附图说明
[0035]
图1为本发明中流程示意图；
[0036]
图2为本发明各站点极值序列的自相关结果；
[0037]
图3为本发明非参数化mann-kendall趋势检测结果；
[0038]
图4为本发明非参数化mann-kendall趋势检测统计变量z值得空间分布经验密度图；
[0039]
图5为本发明特定站点的后验分布参数计算过程；
[0040]
图6为本发明特定站点的后验分布的趋势模型拟合效果图；
[0041]
图7为本发明基于bay_lr趋势检测法进行参数化的极值序列检测结果；
[0042]
图8为本发明参数非单调性结果；
[0043]
图9为本发明参数化和非参数化对于样本长度的敏感性分析结果。
具体实施方式
[0044]
下面结合附图对本发明的技术方案作进一步说明。
[0045]
为便于理解本发明，首先对水文气象极值序列非参数趋势检验法—mann-kendall趋势检验法做以下说明：
[0046]
对于双边检验，零假设是数据遵循不存在单调趋势而为了能够判断极值序列是否服从还是服从假设，存在单调趋势，mann-kendall检验统计量根据如下公式求得：
[0047][0048]
式中，xj和xk分别为极值序列中第j和第k个观测值；sgn(.)为依据xj与xk相对大小的符号判断函数；同时统计量s的方差可以按照如下公式求得：
[0049][0050]
式中，p是数据集中观测值相同的组数，tj是第j个数值相同组。那么标准化后的z值统计变量可以按如下公式计算：
[0051][0052]
在双边趋势检验中，对于给定的置信水平(显著性水平)α，若|z|≥z
1-α/
，则原假设是不可接受的，即在置信水平α(显著性检验水平)上，时间序列数据存在明显的上升或下降趋势。z为正值表示上升趋势，负值表示减少趋势，z的绝对值在大于等于1.645，1.96，2.576时表示分别通过了置信度90％，95％，99％的显著性检验。
[0053]
如图1所示，本发明的一种参数化的水文气象极值序列的趋势检测方法，包括以下步骤：
[0054]
步骤1.采用广义极值分布gev函数建立水文气象极值序列的趋势模型
[0055]
首先采用极值分析领域常用的广义极值分布函数gev进行水文气象极值序列的拟合，
[0056][0057]
式中，x
t
为水文气象数据，gev分布的位置μ
t
，尺度σ
t
和形状参数κ
t
是随着时间t逐年变化的；
[0058]
位置μ
t
的趋势结构包括：
[0059][0060]
式中，μ0、μ1、μ2、μ3为未知待定参数。
[0061]
尺度σ
t
的趋势结构包括：
[0062][0063]
式中，σ0、σ1、σ2、σ3为未知待定参数。
[0064]
形状参数κ
t
的趋势结构包括：
[0065][0066]
式中，k0、k1、k2、κ3为未知待定参数。
[0067]
水文气象极值序列的趋势模型一共包括63种趋势模式结构和1种平稳即无趋势模式。
[0068]
步骤2：基于贝叶斯推断框架估计水文气象极值序列的趋势模型的参数；
[0069]
(1)首先对水文气象极值序列平稳模式的参数进行估计。通过参数先验分布模式假定，基于gibbs采样算法的后验分布推求进行。
[0070]
(2)接着对水文气象极值序列63种趋势模式下的参数进行估计。包括如下步骤：
[0071]
(21)参数先验分布的结构如下：
[0072][0073]
式中，位置和形状参数的先验分布模式服从均值为0方差为1000的正太分布，而尺度参数服从参数为25的柯西分布。
[0074]
(22)基于gibbs采样算法进行后验分布的求解:
[0075]
(221)为参数集φ
t
设定初始化值φ0。
[0076]
(222)将初始化值代入到先验分布模型中求得先验分布的似然值lr
prio
,进而基于观测值样本序列获得gev分布似然值lr
gev
,进而得到后验密度的对数似然值为前两者之和。
[0077]
(223)不断地更新迭代，重复m步，并结合gibbs采样算法当后验密度对数似然值逐渐趋近定值时，采用算法终止。
[0078]
(23)基于m个后验分布参数集求得φ
t
的中值和参数上下界。获得水文气象极值序列的趋势模型参数估计值。
[0079]
步骤3：基于bay_lr趋势检测法进行参数化的水文气象极值序列检测
[0080]
获得了27种情形下的参数估计结果，接下来bay_lr趋势检测法可以按照如下公式求解贝叶斯对数似然比bay_lr
[0081]
bay_lr＝2(logl
1-logl0)
[0082]
其中,logl0为平稳模型的对数似然值而logl1代表了其余26种趋势性模型的对数似然值；由于以上的bay_lr统计量收敛于(卡方)分布，其中卡方分布参数δ为平稳模型和非平稳模型间的差异，为此可以用于检验零趋势假设(分布参数没有趋势)是否存在的重要依据：bay_lr的值大于的α分位数，上述零假设可在α显著性水平上被拒绝，也就代表所研究的水文气象极值序列检测序列存在趋势性。当α＝0.95时，就可以获得95％显著水平的水文气象极值序列趋势性，bay_le会返回固定的p-value值，p-vlaue值小于5％表明具有95％显著水平的趋势性。
[0083]
为了有效验证本发明提出的参数化趋势检测方法的有效性，给出以下实施例。以中国范围内839个气象站点的年极大日降雨序列作为研究对象。以中国范围内的一般气象站点观测到的日降雨量为基础数据，按照年极大值提取法获得年极大日降雨序列，样本长度为1951-2018年一共68年的数据长度。
[0084]
(1)非参数的mann-kendall趋势检测法的识别结果分析
[0085]
为了对比分析本发明的参数化的水文气象极值序列趋势检测法的效果，计算出了非参数的mann-kendall趋势法的检测结果。但由于mann-kendall的检验结果很大程度受到序列的自相关的影响，为此需要先分析降雨极值序列的自相关情况，计算结果如图2所示。839个站点中acf的p-value值都是大于0.05的表明大部分降雨序列自相关性较弱可以采用mann-kendall趋势检验法对原始极值序列进行计算分析。
[0086]
由图3可知，839个气象站点的年极大日雨量序列只有69个站点具有95％显著水平的上升或者下降趋势性，只占到了69/839＝8.2％。图4中z值的空间分布密度函数图也很好地体现了这一点。
[0087]
(2)基于贝叶斯的参数估计结果
[0088]
由于站点众多，且趋势模型种类众多，本节只是显示了站点编号为52203(42.5
°
n,93.3
°
e)的最优趋势模型的分析过程，如图5和图6所示。图5说明了由于有效的先验分布假定使得后验分布在特定的采样步骤之后能够收敛于定值，图6也很好地体现了最优趋势模型估计出的理论分布模型和经验分布最为接近的。
[0089]
(3)基于参数化的水文气象极值序列趋势性检验结果
[0090]
图7显示了最终的水文气象极值序列bay_lr趋势检测法的优选结果。由图7中结果可知，一共有330个站点体现出单调性或者非单调性趋势，总体来看参数化的趋势性检测方
法可以有效检测出330/839＝39.5％。且bay_lr趋势检测法检测出的站点编号大部分都包括了mann-kendall趋势检测的站点编号，换句话说，本发明提出bay_lr趋势检测法不但能够检测出mann-kendall检测出的单调性结果，同时也可以有效识别出多余非单调趋势和参数间存在反向趋势的情况，相比mann-kendall检测法效果提升了39.5％-8.2％＝31.3％。同时由于参数非单调性趋势和参数间反向趋势的存在图8中所示使得mann-kendall检测效果不佳，无法识别此类的极值序列趋势性，也是参数化的bay_lr趋势检测法能够比mann-kendall检测法多出31.3％的原因所在。
[0091]
(4)非参数法和参数法对于样本长度的敏感性分析
[0092]
图9中显示出了两种水文气象极值序列趋势检测方法对于样本长度的敏感性响应关系。图9中上面4幅图是参数化的趋势检测方法分别对选定的4个站点的降雨极值序列按照样本长度由41变化到60时，步长为1，bay_lr趋势检测法检测出的趋势性模型，图中lr等于0代表了序列没有趋势性，而1-26分别代表了26种的趋势性模式结构，可见总的来看非线性趋势1比较常见，而三次型的非线性趋势2模式基本上没有符合的站点。而图9下面4幅图分别是mann-kendall趋势检测法对于样本长度的敏感性分析结果，可以发现随着样本长度增加，mann-kendall的关键统计量z值是逐步增大的，并且由没有趋势逐步过渡到趋势性的情形，可见其对于样本长度过于敏感，而参数化的bay_lr趋势检测法则表现出较弱的敏感性，不管是样本长度增大还是减小，序列的趋势性都能够有效得检测出来。
技术特征：1.一种参数化的水文气象极值序列的趋势检测方法,其特征在于：包括以下步骤：(1)采用广义极值分布gev函数建立水文气象极值序列的趋势模型；(2)基于贝叶斯推断框架估计水文气象极值序列的趋势模型的参数；(3)根据趋势模型的参数估计值获得贝叶斯对数似然比，然后得出水文气象极值序列的趋势。2.根据权利要求1所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：步骤(1)的趋势模型：式中，x
t
为水文气象数据，gev分布的位置μ
t
，尺度σ
t
和形状参数κ
t
是随着时间t逐年变化的。3.根据权利要求2所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：水文气象极值序列的趋势模型gev分布的位置μ
t
的趋势结构包括：式中，μ0、μ1、μ2、μ3为未知待定参数。4.根据权利要求2所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：水文气象极值序列的趋势模型gev分布的尺度σ
t
的趋势结构包括：式中，σ0、σ1、σ2、σ3为未知待定参数。5.根据权利要求2所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：水文气象极值序列的趋势模型gev分布的形状参数κ
t
的趋势结构包括：式中，κ0、κ1、κ2、κ3为未知待定参数。6.根据权利要求1所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：步骤(2)包括以下步骤：(21)趋势模型的参数先验分布的假定；(22)基于gibbs的蒙特卡洛随机采样法获得潜在的后验分布参数；
(23)获得水文气象极值序列的趋势模型参数估计值。7.根据权利要求6所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：步骤(2)中假设先验分布参数集为β,后验分布按照公式(5)获得：式中，f(x
t
|β,t)各时间步时变gev参数的后验分布概率密度值。8.根据权利要求1所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：步骤(3)中获得贝叶斯对数似然比bay_lr＝2(logl
1-logl0)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)其中,logl0为平稳模型的对数似然值，而logl1代表了其余趋势性模型的对数似然值。9.根据权利要求8所述的参数化的水文气象极值序列的趋势检测方法,其特征在于：步骤(3)中bay_lr的值大于的α分位数，则水文气象极值序列存在趋势性；当α＝0.95时，获得95％显著水平的水文气象极值序列趋势性。
技术总结本发明公开了一种参数化的水文气象极值序列的趋势检测方法,包括以下步骤：(1)采用广义极值分布GEV函数建立水文气象极值序列的趋势模型；(2)基于贝叶斯推断框架估计水文气象极值序列的趋势模型的参数；(3)根据趋势模型的参数估计值获得贝叶斯对数似然比，然后得出水文气象极值序列的趋势。本发明从参数化角度有效保证水文气象极值序列的趋势检测精度；基于贝叶斯统计推断框架的使用一方面充分保证复杂趋势模式情形下对于水文气象极值序列的分布参数估计精度不会因为引入过多的参数而受到影响，另一方面采用Gibbs采样算法简化了水文气象极值序列分布参数的估计过程。水文气象极值序列分布参数的估计过程。水文气象极值序列分布参数的估计过程。
技术研发人员：徐鹏程 仇建春 李帆 刘赛艳 蒋新跃
受保护的技术使用者：扬州大学
技术研发日：2022.01.17
技术公布日：2022/7/5
