一种基于梯度参数寻优的降水正态化分析方法及系统

1.本发明涉及水文数据处理技术领域，更具体地，涉及一种基于梯度参数寻优的降水正态化分析方法及系统。


背景技术：

2.降水数据是重要的水文气象观测资料，围绕降水数据开展建模分析，是开发降水数据产品、分析流域干旱事件和开展水文预报的有效途径。受降水的自然属性影响，降水数据通常呈现非正态分布：一方面，降水通常呈现正偏分布，具有偏度与峰度较高的特点；另一方面，降水具有天然的下边界，即降水量的最小值为零，导致降水数据呈现离散-连续的混合分布。而目前许多统计分析方法是在正态分布的前提下开展理论推导，因此降水数据的非正态特征将导致其建模分析过程更加复杂，并对统计分析的结果产生一定的影响。
3.面向降水的非正态特征，目前常用的方法是采用如log变换、box-cox变换与log-sinh变换等正态变换方法将非正态分布的降水数据转化为服从正态分布的数据，再进一步进行建模分析。其中，不同的变换方法具有不同的变换参数，并且各个参数对于正态变换的影响不同，常用的方法是根据经验设置变换参数，但经验性的参数设置难以适应不同气候条件下的降水分布特征，由此得到的降水数据分析结果的准确性有待提高。


技术实现要素：

4.本发明为克服上述现有技术所述的降水正态化分析方法难以适应不同气候条件下的降水分布特征，数据分析准确性有待提高的缺陷，提供一种基于梯度参数寻优的降水正态化分析方法及系统。
5.为解决上述技术问题，本发明的技术方案如下：
6.一种基于梯度参数寻优的降水正态化分析方法，包括以下步骤：
7.s1、获取待分析的降水数据；
8.s2、构建正态变换模型，对所述降水数据进行正态变换，得到正态变量z；其中，所述正态变换模型中包括相应的正态变换参数；
9.s3、令所述正态变量z服从正态分布，构建所述正态变量z的联合概率密度函数；
10.s4、基于所述正态变换模型和所述联合概率密度函数，构建用于参数寻优的似然函数；其中，待寻优的参数包括正态分布参数和所述正态变换参数；
11.s5、通过推导所述似然函数的解析梯度向量对所述似然函数进行寻优，至满足预设的终止条件，得到满足使所述似然函数取最大值的最优参数；
12.s6、基于所述最优参数对所述正态变换模型进行更新，并对所述降水数据进行正态变换及建模分析，得到降水正态化分析结果。
13.本技术方案中，通过推导似然函数的解析梯度向量对所述似然函数进行寻优，简化正态变换的参数寻优过程，同时完成不同正态变换的参数估计，得到适应不同气候条件下的降水分布特征的降水正态化分析结果。
14.进一步地，本发明还提出一种基于梯度参数寻优的降水正态化分析系统，应用于上述降水正态化分析方法。其中，降水正态化分析系统包括：
15.数据采集模块，用于获取待分析的降水数据；
16.正态变换模块，用于根据预设的正态变换模型，对所述降水数据进行正态变换，得到正态变量z；
17.正态分布模块，用于使正态变量z服从正态分布，构建所述正态变量z的联合概率密度函数；
18.寻优模块，用于基于所述正态变换模型和所述联合概率密度函数，构建用于参数寻优的似然函数，通过推导似然函数的梯度向量解析表达式对所述似然函数进行寻优，至满足预设的终止条件，得到满足使所述似然函数取最大值的最优参数后，对所述正态变换模块中的正态变换模型进行更新；
19.分析模块，用于根据经过寻优更新后的正态变换模块输出的正态变量z进行建模分析，输出降水正态化分析结果。
20.与现有技术相比，本发明技术方案的有益效果是：本发明通过推导似然函数的梯度向量解析表达式并采用极大似然估计方法进行寻优，可根据降水的不同分布特性自适应的进行参数寻优，以适应不同气候条件下的降水分布特征，进而降低了水文气象工作者开展降水正态化工作的难度。
附图说明
21.图1为实施例1的基于梯度参数寻优的降水正态化分析方法的流程图。
22.图2为实施例2的正态变换前后标准化降水数据的频率分布直方图。
23.图3为实施例2的正态变换前后标准化降水的正态分布分位图。
24.图4为实施例3的降水正态化分析系统的架构图。
具体实施方式
25.附图仅用于示例性说明，不能理解为对本专利的限制；
26.为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；
27.对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
28.下面结合附图和实施例对本发明的技术方案做进一步的说明。
29.实施例1
30.本实施例提出一种基于梯度参数寻优的降水正态化分析方法，如图1所示，为本实施例的降水正态化分析方法的流程图。
31.本实施例提出的基于梯度参数寻优的降水正态化分析方法中，包括以下步骤：
32.s1、获取待分析的降水数据。
33.s2、构建正态变换模型，对所述降水数据进行正态变换，得到正态变量z；其中，所述正态变换模型中包括相应的正态变换参数。
34.s3、令所述正态变量z服从正态分布，构建所述正态变量z的联合概率密度函数。
35.s4、基于所述正态变换模型和所述联合概率密度函数，构建用于参数寻优的似然函数。
36.其中，待寻优的参数包括正态分布参数和所述正态变换参数。
37.s5、通过推导所述似然函数的解析梯度向量对所述似然函数进行寻优，至满足预设的终止条件，得到满足使所述似然函数取最大值的最优参数。
38.s6、基于所述最优参数对所述正态变换模型进行更新，并对所述降水数据进行正态变换及建模分析，得到降水正态化分析结果。
39.本实施例中通过构建似然函数并通过推导所述似然函数的解析梯度向量进行寻优，可根据降水的不同分布特性自适应的进行参数寻优，以适应不同气候条件下的降水分布特征，进而降低了水文气象工作者开展降水正态化工作的难度。
40.在一可选实施例中，所构建的正态变换模型基于log变换、box-cox变换或log-sinh变换中的一种或多种。
41.其中，记x＝[x1,x2,..,xn]为降水数据的n个样本，而z＝[z1,z2,..,zn]表示经过正态变换后相应的正态变量。
[0042]
则基于log变换的正态变换模型的表达式如下：
[0043]zlog
(x；c)＝log(x+c)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0044]
式中，log(
·
)表示自然对数函数，c为log变换的参数，常取非负数，用于处理x＝0时log变换无意义的情况。
[0045]
log变换对x的一阶导数为：
[0046]z′
log
(x；c)＝(x+c)-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0047]
由此可见，基于log变换的正态变换模型中仅有1个待寻优参数，为参数c。
[0048]
基于box-cox变换的正态变换模型的表达式如下：
[0049][0050]
式中，z
box-cox
(
·
)表示经过box-cox变换的正态变量集合，λ1、λ2为box-cox变换的正态变换参数。
[0051]
其中，λ1的取值范围为[-2,2]。且从式(3)可知，当λ1＝0时，box-cox变换等同于log变换，此时参数λ2的作用与参数c相同，用于处理x＝0时box-cox变换无意义的情况。λ2常取非负数，同时也可以固定λ2为0或其他正数。
[0052]
box-cox变换对x的一阶导数为：
[0053][0054]
由此可见，基于box-cox变换的正态变换模型中的待寻优参数为λ1、λ2。
[0055]
基于log-sinh变换的正态变换模型的表达式如下：
[0056][0057]
式中，z
log-sinh
(
·
)表示经过log-sinh变换的正态变量集合，α、β为log-sinh变换的正态变换参数。其中，当β参数趋近于无穷大时，log-sinh变换与log变换的作用相似。
[0058]
log-sinh变换对x的一阶导数为：
[0059][0060]
coth(
·
)表示双曲余切函数。由此可见，基于log-sinh变换的正态变换模型中的待寻优参数为α、β。
[0061]
在一可选实施例中，所述s3步骤中，包括以下步骤：
[0062]
s3.1、将所述降水数据中小于或等于删失阈值x0的数值视为删失值；
[0063]
s3.2、假设经过删失处理后的正态变量z服从正态分布，用于构建所述联合概率密度函数。
[0064]
进一步地，所述删失阈值x0为等于0或略大于0的实数。
[0065]
本实施例考虑到由于降水的下边界为零，因此降水数据呈现离散-连续混合分布。而在传统的降水正态化过程中，通常采用添加偏移系数的方式来处理降水的零值，但未考虑降水的混合分布对参数估计的影响。本实施例通过基于删失阈值对降水数据进行处理，将降水数据转为连续分布，相比传统的处理方式能够完整的考虑降水的混合分布对参数估计的影响，估计结果更合理。
[0066]
进一步假设经过删失处理后的正态变量z服从正态分布，构建联合概率密度函数：
[0067][0068]
式中，zi∈z表示正态变量z中的第i个经过正态变换的降水数据样本；μz和σz表示正态变量z服从正态分布的均值与标准差；pn(
·
)表示正态变量z服从正态分布的概率密度函数，φn(
·
)表示正态变量z的累积分布函数；ω1表示降水数据大于删失阈值x0的样本索引集合，数量记为n1，ω0表示降水数据小于或等于删失阈值x0的样本索引集合，数量记为n0，且n＝n0+n1。
[0069]
基于式(7)，得到用于参数寻优的似然函数的表达式如下：
[0070][0071]
式中，θ表示似然函数p(x|θ)中的参数集合，包括正态分布参数μz和σz以及正态变换参数；j表示正态变换的雅可比矩阵。
[0072]
对于似然函数p(x|θ)，常取其对数形式，即可得到：
[0073][0074]
式中，z
′
(xi)表示相应正态变换的一阶导数，对于log、box-cox与log-sinh变换，分别采用式(2)、式(4)与式(6)中的形式；erf(
·
)表示误差函数。
[0075]
进一步地，本实施例对所述似然函数及其梯度信息进行寻优，至满足预设的终止条件，得到满足使所述似然函数取最大值的最优参数。其中，待寻优的参数包括正态分布参数μz和σz，以及正态变换参数c，λ1、λ2，α和β。
[0076]
本实施例采用极大似然估计方法进行寻优，即找到一组参数使得式(9)中的logp(x|θ)取得最大值。
[0077]
进一步地，在一可选实施例中，s5步骤中，其具体步骤包括：
[0078]
s5.1、设置待寻优的参数的起始点θ0；
[0079]
s5.2、采用拟牛顿法，基于梯度向量对所述似然函数进行迭代寻优，至满足预设的终止条件，得到满足使所述似然函数取最大值的最优参数。
[0080]
其迭代求解公式如下：
[0081][0082]
式中，θ
k+1
与θk表示在第k+1与k次迭代过程中待寻优的参数的取值；gk表示参数集合θ在似然函数中组成的梯度向量在第k次迭代中的取值；表示第k次迭代中海森矩阵的逆矩阵。
[0083]
本实施例考虑到传统的全局寻优算法计算量较大，寻优所需时间较长，并且受到局部最优值的影响，采用拟牛顿法作为寻优算法，同时，基于式(9)的对数似然函数logp(x|θ)，推导出关于不同参数的梯度的解析解，以此作为梯度信息，目的在于为算法的搜索提供方向，提高算法的搜索效率，从而较快的寻找到最优参数值。常用的方法包括dfp算法(davodpm-fletcher-powell)和bfgs算法(broyden-fletcher-goldfard-shano)等。
[0084]
进一步地，对数似然函数的梯度由对数似然函数logp(x|θ)关于参数的一阶偏导构成。对于log、box-cox与log-sinh变换，均需要估计正态化变量z的均值μz和标准差σz。
[0085]
其中，均值μz的一阶偏导表示为：
[0086][0087]
标准差σz的一阶偏导表示为：
[0088][0089]
而对于不同的正态变换方法具有不同的参数。对于log变换，logp(x|θ)关于参数c的一阶偏导为：
[0090][0091]
对于box-cox变换，logp(x|θ)关于参数λ1的一阶偏导为：
[0092][0093]
当λ1＝0时，logp(x|θ)关于参数λ1的一阶偏导为：
[0094][0095]
对于box-cox变换，logp(x|θ)关于参数λ2的一阶偏导为：
[0096][0097]
对于log-sinh变换，logp(x|θ)关于参数α的一阶偏导为：
[0098][0099]
对于log-sinh变换，logp(x|θ)关于参数β的一阶偏导为：
[0100][0101]
基于各参数的一阶偏导数，可得到对数似然函数的梯度向量，具体为：
[0102]
对于log变换，对数似然函数的梯度向量为：
[0103][0104]
对于box-cox变换，对数似然函数的梯度向量为：
[0105][0106]
对于log-sinh变换，对数似然函数的梯度向量为：
[0107][0108]
进一步地，在一可选实施例中，s5.1步骤中，设置待寻优的参数的起始点θ0的步骤包括：
[0109]
设置参数θ的初始估计值θ0，并假设参数θ0在其范围内服从均匀分布：
[0110]
θ0～u(b
l
,bu)
[0111]
其中b
l
与bu分别表示参数θ的下边界与上边界；
[0112]
从参数θ0的均匀分布中随机抽取若干随机点作为拟牛顿法的初始点进行求解。
[0113]
本实施例中，考虑似然函数可能存在多个局部最优解，从参数θ0的均匀分布中随机抽取若干随机点作为拟牛顿法的初始点进行求解，最后从多个求解结果中，选择使得-logp(x|θ)取得最小值(即使得logp(x|θ)取得最大值)的参数组合，作为最终得到的参数优化结果θ
opt
：
[0114][0115]
式中，参数优化结果θ
opt
包括正态分布参数和正态变换参数的优化取值。
[0116]
进一步地，在一可选实施例中，s5.2步骤中，所述终止条件包括以下至少一种条件：
[0117]
(1)当前迭代中梯度向量的取值gk小于预设阈值εg；
[0118]
(2)在两次迭代过程中似然函数取值变化小于预设阈值ε
p
。
[0119]
条件(1)中，当梯度向量||gk||小于阈值εg时，认为似然函数已经收敛，则当前的参数集合θk为所求解结果，即为最优参数。
[0120]
条件(2)中，表示为：此时认为似然函数已经收敛，当前的参数集合θk为所求解结果，即为最优参数。
[0121]
本实施例适用于log、box-cox以及log-sinh变换的参数寻优，通过python编程语言实现其数学建模过程，有利于实现自动化降水正态变换。其中，本实施例通过构建似然函数并采用极大似然估计方法进行寻优，推导出关于不同参数的梯度的解析解，以适应不同气候条件下的降水分布特征，从而有效提高降水数据分析结果的准确性。
[0122]
实施例2
[0123]
本实施例应用实施例1提出的一种基于梯度参数寻优的降水正态化分析方法，提出一具体实施过程。
[0124]
本实施例以全球降水数据产品global precipitation climatology centre的月尺度降水为输入数据，对输入的降水数据进行正态变换，包括log、box-cox与log-sinh变换。其具体步骤如下：
[0125]
s1、采用netcdf格式的global precipitation climatology centre数据集，通过
python第三方库xarray中的open_dataset函数读取，进而提取出7月份的全球降水观测数据，得到待分析的降水数据，并存储在名为xr_gpcc_precip的变量中。
[0126]
s2、对s1中采集得到的降水数据进行数学建模分析。
[0127]
其中，设置删失阈值x0＝0.01，存储在名为threshold的变量中，降水数据中小于等于x0的数值均替换为x0，并将其位置索引记录在名为mask的变量中，同时将大于x0与小于等于x0的样本个数分别存储在名为n1和n0的变量中。
[0128]
然后通过numpy与scipy中的log函数完成自然对数计算，sinh函数完成双曲正弦计算，power函数完成指数计算，erf函数完成误差函数计算，从而完成正态变换与似然函数的构建。
[0129]
s3、在s2的基础上，对对数似然函数进行梯度分析，并将数学过程封装成函数。
[0130]
s3.1、根据s2中的对数似然函数形式，推导出其关于不同参数的梯度，包括正态分布参数与正态变换参数，其中关于参数μz、σz、c，λ1、λ2，α和β的梯度计算分别定义为函数grad_mu、grad_sigma、grad_c、grad_l1、grad_l2、grad_alpha和grad_beta。
[0131]
s3.2、将似然函数中联合概率密度计算、log变换雅可比矩阵计算、box-co x变换雅可比矩阵计算和log-sinh变换雅可比矩阵计算分别定义为函数norm、j ac_log、jac_boxcox、jac_logsinh，并结合s3.1中函数封装为类likelihood funct ion。
[0132]
s3.3、分别定义log、box-cox和log-sinh变换为函数log、boxcox和logsinh函数，并封装成类powertrans。
[0133]
s3.4、将以上定义的函数与类保存成为.py格式的文件。
[0134]
s4.通过python中import语句调用s3中封装的函数，对global precipitation climatology centre中逐个网格进行降水数据正态变换。
[0135]
s4.1、设置多个初始点，即各个参数不同的初始估计值，基于s2与s3中对数似然函数及其梯度，通过scipy中的optimize.minimize函数采用拟牛顿法进行参数寻优。
[0136]
s4.2从不同初始点得到的结果中，选择对数似然函数取值最大的一组，其相应的参数寻优结果作为降水正态变换的参数。
[0137]
s4.3、依次对不同的网格与不同的正态变换(log、box-cox和log-sinh变换)进行应用，完成全部降水的正态变换。
[0138]
s4.4、以其中三个网格为例，通过matplotlib中的pyplot.hist函数绘制正态变换前后标准化降水数据的频率分布直方图，如图2所示。
[0139]
图2中，第1行表示的原始降水数据整体呈现正偏的特点，而第2～4行的结果呈现出左右对称结果，表明本发明能够有效估计正态变换参数，较好的实现降水数据正态化。
[0140]
s4.5、通过pyplot.scatter函数绘制变换前后标准化降水的正态分布分位图，如图3所示。
[0141]
图3中，第一行表示降水原始数据的结果，可以看到散点图偏离于1：1线，表明原始数据不符合正态分布，而2～4行表示的不同变换的结果，整体沿1：1线分布，表明基于本发明得到的正态化变量服从正态分布。
[0142]
s4.6、采用scipy中的skew、kurtosis、shapiro和pearsonr函数以及numpy对正态变换后的降水进行正态性检验，并通过basemap绘制计算得到的偏度系数、峰态系数、shapiro-wilk检验p值和filliben r statistic值。
[0143]
从不同的正态性检验指标可以看出，相比于原始降水数据，正态变换后的降水的正态性整体上有了显著的提高，表明本发明能够有效估计log、box-cox和log-sinh不同正态变换的参数，以适应全球不同的气候特点，具有良好的应用效果和稳定性。
[0144]
实施例3
[0145]
本实施例提出一种基于梯度参数寻优的降水正态化分析系统，应用于实施例1提出的基于梯度参数寻优的降水正态化分析方法。如图4所示，为本实施例的降水正态化分析系统的架构图。
[0146]
本实施例提出的基于梯度参数寻优的降水正态化分析系统中，包括：
[0147]
数据采集模块，用于获取待分析的降水数据。
[0148]
正态变换模块，用于根据预设的正态变换模型，对所述降水数据进行正态变换，得到正态变量z。
[0149]
正态分布模块，用于使正态变量z服从正态分布，构建所述正态变量z的联合概率密度函数。
[0150]
寻优模块，用于基于所述正态变换模型和所述联合概率密度函数，构建用于参数寻优的似然函数，通过推导似然函数的梯度向量解析表达式对所述似然函数进行寻优，至满足预设的终止条件，得到满足使所述似然函数取最大值的最优参数后，对所述正态变换模块中的正态变换模型进行更新。
[0151]
分析模块，用于根据经过寻优更新后的正态变换模块输出的正态变量z进行建模分析，输出降水正态化分析结果。
[0152]
本实施例中，通过寻优模块对log、box-cox以及log-sinh变换的参数寻优，推导出关于不同参数的梯度的解析解，以适应不同气候条件下的降水分布特征，从而有效提高降水数据分析结果的准确性。
[0153]
进一步地，在一可选实施例中，所述正态变换模块包括log变换单元、box-cox变换单元和log-sinh变换单元中的至少一个。
[0154]
进一步地，在一可选实施例中，所述数据采集模块还用于根据预设的删失阈值x0，将所述降水数据中小于或等于删失阈值x0的数值视为删失值，然后发送至正态变换模块中得到相应的正态变量z，再输入正态分布模块使经过删失处理后的正态变量z服从正态分布，用于进一步构建正态分布概率密度函数。
[0155]
进一步地，在一可选实施例中，所述寻优模块采用拟牛顿法，基于梯度向量对所述似然函数进行迭代寻优。其中，从参数θ0的均匀分布中随机抽取若干随机点作为拟牛顿法的初始点进行求解，最后从多个求解结果中，选择使得-logp(x|θ)取得最小值(即使得logp(x|θ)取得最大值)的参数组合，作为最终得到的参数优化结果θ
opt
。
[0156]
相同或相似的标号对应相同或相似的部件；
[0157]
附图中描述位置关系的用语仅用于示例性说明，不能理解为对本专利的限制；
[0158]
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。