CN113869710A - 一种基于大规模气候异常的极端降水趋势分析方法 - Google Patents

Abstract

本发明公开了属于水文水资源技术领域的一种基于大规模气候异常的极端降水趋势分析方法。包括步骤1：收集气象水文观测数据和大尺度气候异常指标数据，再整理成完整的气象水文信息资料；步骤2：对年最大径流序列方差进行测试，识别年最大径流序列数据中的变异点，确定变异点位置；步骤3：确定年最大径流时间序列是否具有统计上的显著趋势；步骤4：对最大径流序列的统计参数进行线性以及非线性趋势变化分析；步骤5：利用赫斯特指数H估计径流时间序列的长期平稳性特征。本发明对流域历史极端降水的变化进行了全面调查，并对基于流域年最大径流序列的所有统计参数进行线性以及非线性趋势变化分析，解决了水文资料短缺地区资料移用的问题。

Description

一种基于大规模气候异常的极端降水趋势分析方法

技术领域

本发明涉及水文水资源技术领域，尤其涉及一种基于大规模气候异常的极端降水趋势分析方法。

背景技术

气候变化和人类活动是影响流域水文循环过程和水资源演变规律的两大驱动因素，所引起的水文效应已成为当前全球变化研究领域的焦点问题。气候变化对降水、蒸发、径流、土壤湿度等造成直接影响，并引起水资源在时间和空间上的重新分配以及水资源总量的改变，进而改变关键基础设施的设计风险。对人类社会水资源的开发、利用以及规划和管理等诸多环节造成严重影响，并进一步影响生态环境与社会经济的可持续发展。因此，考虑到气候变化的可能影响，传统水文频率分析基础的一致性假设将不再成立，公共安全标准可能会受到影响。世界各地发生的大洪水凸显了解决非平稳性和未来水文极端事件可能发生的变化的必要性，这些极端事件会影响水安全、水资源管理和大型水坝的运行。以往有大量研究针对长时间水文序列进行趋势性检验(如Mann-Kendall秩次相关检验)，周期性检验(如周期图法，连续谱分析法，小波分析法)，突变性检验(如Hurst系数，R/S分析，双累积曲线法)等，但往往在对某地区极端降水的趋势分析中没有考虑大尺度气候异常所导致的可能影响，即极端降水可能与大尺度气候异常有关，以及无资料地区的资料移用问题。

因此有必要发明一种基于大规模气候异常的极端降水趋势分析方法，以克服现有技术存在的缺陷，并在此基础上重点检测在大规模气候异常情况下极端降水的变化趋势和变异点，以及时间序列持续性对降水趋势的可能影响，同时解决水文资料短缺地区资料移用问题。这对水资源的合理利用和规划管理具有重要的现实意义，对进一步研究气候变化对流域水资源与防洪安全的影响及区域可持续发展具有重要的指导意义。

发明内容

本发明的目的是提出一种基于大规模气候异常的极端降水趋势分析方法，其特征在于，包括以下步骤：

步骤1：气候水文数据收集及预处理；收集气象水文观测数据和大尺度气候异常指标数据，再整理成完整的气象水文信息资料；

步骤2：变异点分析；通过统计检验方法对年最大径流序列方差进行测试，识别年最大径流序列数据中的变异点，并通过调查实地气候水文现象，确定变异点位置；

步骤3：时间趋势分析；利用非参数Mann-Kendall检验、Spearman秩相关检验以及Pearson检验确定年最大径流时间序列是否具有统计上的显著趋势；

步骤4：水文频率分析；利用GAMLSS模型对最大径流序列的统计参数进行线性以及非线性趋势变化分析；

步骤5：长期平稳性分析；利用赫斯特指数H估计径流时间序列的长期平稳性特征。

所述步骤1中的大尺度气候异常指标数据包括南方涛动指数、北大西洋涛动指数、太平洋十年涛动指数和北太平洋指数。

所述步骤2中的统计检验方法包括Pettitt检验、累积总和检验、W秩和检验、线性回归中的结构变化检验、贝叶斯变异点检验和小波分析。

所述步骤3中的非参数Mann-Kendall检验具体如下：

假设有n个样本量(x₁,…,x_n)的时间序列，对于所有k,j≤n，且k≠j，x_k和x_j的分布是不同的，计算检验统计量s为：

其中

当n>10时，标准的正态统计变量通过下式计算：

其中，S为正态分布，均值为0；方差var(s)＝n(n-1)(2n+5)/18；

对于统计值Z，Z大于0时，是增加趋势；Z小于0时，是减少趋势；Z的绝对值在大于1.28、1.64和2.32时，分别表示通过了置信度为90％、95％和99％的显著性检验。

所述步骤4中的GAMLSS模型具体如下：

假设某一时刻t(t＝1,2,…,n)的随机变量观测值y_t服从概率密度函数f(y_t|θ^t)，θ＝(θ_t1,θ_t2,…,θ_tp)是时刻t对应的分布/统计参数向量，p是分布参数的个数，n是观测值的个数；记y＝(y₁,y₂,…,y_n)^T为独立观测值组成的向量，θ_k为所有时刻的第k个分布参数组成的向量，θ_k＝(θ_1k,θ_2k,…,θ_nk)^T，k＝1,2,…,p；记gk(·)表示θ_k与相应的解释变量X_k和随机效应项之间的单调函数关系，表示为：

式中：η_k为长度为n的向量，β_k＝(β_1k,β_2k,…,β_Ikk)^T为长度为I_k的回归参数向量，X_k为n×I_k的解释变量矩阵，Z_jk为一个已知的n×q_jk固定设计矩阵，γ_jk是一个q_jk维的正态分布随机变量向量，Z_jkγ_jk表示第j项随机效应，q_jk表示第j项随机效应中的随机影响因子维数；前两个参数θ₁和θ₂定义为位置参数向量和尺度参数向量；用v和τ表示随机变量序列的偏度向量和峰度向量；

如果不考虑随机效应对分布参数的影响，即对于k＝1,2,…,p时，令J_k＝0，那么GAMLSS模型为全参数模型：

g_k(θ_k)＝η_k＝X_kβ_k

如果假定随机变量Y服从三参数概率分布，那么GAMLSS模型变成：

g₁(μ)＝X₁β₁

g₂(σ)＝X₂β₂

g₃(ν)＝X₃β₃

当研究随机变量分布参数的变化与时间t的关系时，解释变量矩阵为：

得到分布参数与时间变量的函数关系：

GAMLSS模型关于回归参数β的似然函数为：

以似然函数取值最大为目标函数，采用RS算法来估计回归参数β的最优值；GAMLSS模型的全局拟合偏差GD定义为：

式中

为回归参数估计值所对应的对数似然函数；同时引入广义AIC准则进行判断，防止模型的过度拟合，其定义为：

GAIC＝GD+#df

式中：df为模型中的整体自由度，#是惩罚因子；GAIC值最小的模型作为最优模型。

所述步骤5中利用赫斯特指数H估计径流时间序列的具体方法如下：

若H＝0.5，径流时间序列用随机游走来描述；若0.5<H<1，径流时间序列存在长期记忆性；若0≤H<0.5，径流时间序列为反持续性。

本发明的有益效果在于：

本发明弥补了传统长期径流序列变异点分析、趋势分析、水文频率分析中缺乏考虑大规模气候异常现象的不足；对流域历史极端降水的变化进行了全面调查，并对基于流域年最大径流序列的所有统计参数进行线性以及非线性趋势变化分析，以弥补常规方法在水文序列趋势分析中的不足。此外，通过区域洪水频率分析解决了水文资料短缺地区资料移用的问题。

附图说明

图1为本发明一种基于大规模气候异常的极端降水趋势分析方法的流程图。

具体实施方式

本发明提出一种基于大规模气候异常的极端降水趋势分析方法，下面结合附图和具体实施例对本发明做进一步说明。

图1为本发明基于大规模气候异常的极端降水趋势分析方法的流程图。共分为如下五个步骤，分别为：气候水文数据收集及预处理、变异点分析、时间趋势分析、水文频率分析和长期平稳性分析。

步骤一：气候水文数据收集及预处理

经过良好质量控制的长期气象水文数据集对实现可靠的年径流极值变异点检测和非一致性分析至关重要。应在能力范围内尽可能选取较多的气象水文观测站、站点最短记录大于50年且流域覆盖多种地理类型的数据，以确保研究结果的统计有效性，获得具有统计意义的变异点和趋势性的代表性历史数据。

1.气象水文观测数据收集

1.1地面观测数据：

首先通过国家气象局等相关政府职能部门获取某流域国家级地面站自建站以来的气压、气温、降水量、蒸发量、相对湿度、风向风速、日照时数和地表温度等气象要素的日观测数据；其次通过水利部以及对应各省市、自治区水利部门获取某流域水文站记录的水文及河道要素(如日水位和流量等)、雨量站的降水要素以及水利工程的关键信息(如水库下泄流量等)，采集场次洪水高精度的气象水文观测数据。

1.2实地调查数据：

气象水文观测受到时间、空间的限制，收集的资料往往不能满足实际研究需要。因此将通过气象水文调查(流域调查、水量调查、历史洪水与暴雨调查、枯水和旱情调查以及其他专项调查等)来补充观测的不足，使气象水文资料更加系统、完整。并拟对该流域干支流上的河道形态数据(如河道糙率、河道比降、河床岩性和地下水埋深等)以及基础的下垫面信息(如植被覆盖、土壤类型和高程等)进行调查；此外，还将充分收集该流域水库调度规则及关键信息(如水库历史信息和水库设计标准等)，对所收集的水库群调度格局等关键信息进行调查核实。

2.大尺度气候异常指标数据收集

通过收集国际气候学界广泛研究的南方涛动指数、北大西洋涛动指数、太平洋十年涛动指数和北太平洋指数四种气候异常信息，构建大尺度气候异常数据集。

2.1南方涛动指数Southern Oscillation Index(SOI)：气象学上用南太平洋大溪地与达尔文两地的气压差来衡量，反映了厄尔尼诺现象的活跃程度。具体关系是：南方涛动指数出现持续性的负值，该年有厄尔尼诺现象。相反地，如果南方涛动指数出现持续性的正值,该年有反厄尔尼诺现象(拉尼娜现象)。

2.2北大西洋涛动North Atlantic oscillation(NAO)：用亚速尔群岛的蓬塔代尔加达斯和冰岛的阿库雷里之间的气压基，作为北大西洋涛动指数。北大西洋涛动强，表明两个活动中心之间的气压差大，北大西洋中纬度的西风强，为高指数环流。这时墨西哥湾暖流及拉布拉多寒流均增强，西北欧和美国东南部因受强暖洋流影响，出现暖冬；同时为寒流控制的加拿大东岸及格陵兰西岸却非常寒冷。反之北大西洋涛动弱，表明两个活动中心之间的气压差小，北大西洋上西风减弱，为低指数环流。这时西北欧及美国东南部将出现冷冬，而加拿大东岸及格陵兰西岸则相对温暖。

2.3太平洋十年涛动(Pacific Decadal Oscillation，PDO)：PDO是一种以10年周期尺度变化的以中纬度太平洋盆地为中心的强、周期性的海洋大气气候变化模式。变换周期通常为20～30年。PDO的特征为太平洋北纬20度以北区域表层海水温度异常偏暖或偏冷。在太平洋十年涛动“暖相位”(或“正相位”)期间西太平洋偏冷而东太平洋偏暖，在“冷相位”(或“负相位”)期间西太平洋偏暖而东太平洋偏冷。

2.4北太平洋指数North Pacific Index(NPI)：NPI是Alaska湾附近Aleutian低压区域上平均的月海平面气压异常。北太平洋冬季大气环流的主要模态表现为阿留申低压(AL)的强度变化。NPI是对这一模态的敏感指数，可以用来揭示冷季的北太平洋海平面气压(SLP)变率。

月尺度的SOI，NAO，PDO和NPI指数数据可以从全球气候观测系统Global ClimateObserving System(GCOS)获得。年度和季节的时间序列值由通过计算各个时期的平均值得到。

3.气象水文信息预处理

各种气象水文观测的原始站点数据和大尺度气候异常数据都将按科学的方法和统一的格式整理、分析、统计、提炼成系统、完整的气象水文信息资料，以支持研究变化环境下某流域年径流极值序列的变异点检验、趋势性分析、水文频率分析等问题。

步骤二：变异点分析

通过目前国际流行的多种统计检验方法，包括Pettitt检验，累积总和检验，W秩和检验，线性回归中的结构变化检验，贝叶斯变异点检验和小波分析对年最大径流序列方差进行测试。借助多种检验方式，重点在于识别年最大径流序列数据中主要的变异点。并通过调查实地气候水文现象，确定变异点位置。

步骤三：时间趋势分析

要确定年最大径流时间序列是否具有统计上的显著趋势，可以使用非参数Mann-Kendall检验和Spearman秩相关检验，以及参数最小二乘线性回归(Pearson)检验。在检测到一个变异点后，这些检测被用于分析两个子系列的趋势：一个在原始记录中检测到的变异点之前和一个在之后。如果没有检测到变异点，这些测试将应用于整个序列。

在时间序列趋势分析中，Mann-Kendall检验法是世界气象组织推荐并已广泛使用的非参数检验方法，最初由Mann和Kendall提出，许多学者不断应用Mann-Kendall方法来分析降水、径流、气温和水质等气候水文要素时间序列的趋势变化。Mann-Kendall检验不需要样本遵从一定的分布，也不受少数异常值的干扰，适用于水文、气象等非正态分布的数据，计算简便。

假设有n个样本量(x₁,…,x_n)的时间序列，对于所有k,j≤n,且k≠j，x_k和x_j的分布是不同的，计算检验统计量s，公式如下：

其中：

当n>10时，标准的正态统计变量通过下式计算：

S为正态分布，均值为0；方差var(s)＝n(n-1)(2n+5)/18。对于统计值Z来说，大于0时，是增加趋势；小于0时，则是减少趋势。Z的绝对值在大于1.28、1.64、2.32时，分别表示通过了置信度90％、95％和99％的显著性检验。

Pearson相关系数是最常用的相关系数，又称积差相关系数，取值-1到1，绝对值越大，说明相关性越强。该系数的计算和检验为参数方法，适用条件如下：(1)两变量呈直线相关关系，如果是曲线相关可能不准确；(2)极端值会对结果造成较大的影响；(3)两变量符合双变量联合正态分布。Spearman秩相关系数对原始变量的分布不做要求，适用范围较Pearson相关系数广，即使是等级资料，也可适用。但其属于非参数方法，检验效能较Pearson系数低。

步骤四：水文频率分析

步骤三提到的常用趋势分析方法一般只考虑了水文序列均值的线性趋势变化，这与实际情况是不相符的。事实上，水文序列的趋势会存在于水文序列的任何统计参数，不仅仅是均值，还包括方差、偏态系数等。另外，水文序列的趋势可能是线性的，也可能是非线性的。因此步骤四将引入位置、尺度、形状的广义可加模型(Generalized Additive Modelsfor Location,Scale and Shape，简称GAMLSS)，基于该流域年最大径流序列的所有统计参数进行线性以及非线性趋势变化进行分析，以弥补常规方法在水文序列趋势分析中的不足。

GAMLSS模型是(半)参数回归模型，可以描述随机变量序列的任何统计参数与解释变量之间的线性或非线性关系。同时，GAMLSS模型可以描述的随机变量分布函数类型范围比较广泛，包括一系列的高偏度和高峰度的离散和连续分布，尤其适合拟合具有超峰度和平顶峰度、高度正偏/负偏而不服从传统指数分布的随机变量序列。该模型在经济学、医学、生态学等一些领域中已经得到了较为广泛的应用。

在GAMLSS模型中，假设某一时刻t(t＝1,2,…,n)的随机变量观测值y_t服从概率密度函数f(y_t|θ^t)，θ＝(θ_t1,θ_t2,…,θ_tp)是时刻t对应的分布/统计参数向量，p是分布参数的个数，n是观测值的个数。记y＝(y₁,y₂,…,y_n)^T为独立观测值组成的向量，θ_k为所有时刻的第k个分布参数组成的向量，θ_k＝(θ_1k,θ_2k,…,θ_nk)^T，k＝1,2,…,p。记gk(·)表示θ_k与相应的解释变量X_k和随机效应项之间的单调函数关系，一般表示为：

式中：η_k为长度为n的向量，β_k＝(β_1k,β_2k,…,β_Ikk)^T为长度为I_k的回归参数向量，X_k为n×I_k的解释变量矩阵，Z_jk为一个已知的n×q_jk固定设计矩阵，γ_jk是一个q_jk维的正态分布随机变量向量，Z_jkγ_jk表示第j项随机效应，q_jk表示第j项随机效应中的随机影响因子维数。以上模型的前两个参数θ₁和θ₂通常被定义为位置参数向量和尺度参数向量；当具体采用μ和σ符号时，则分别表示随机变量的均值向量和均方差(或变差系数)向量。如果分布中还有其他参数，那么就被定义为形状参数。对于大多数分布，最多只有两个形状参数，通常用v和τ表示随机变量序列的偏度向量和峰度向量。如果不考虑随机效应对分布参数的影响，即对于k＝1,2,…,p时，令J_k＝0，那么GAMLSS模型就成为一个全参数模型：

g_k(θ_k)＝η_k＝X_kβ_k

如果假定随机变量Y服从三参数概率分布，那么GAMLSS模型的通用表达式就变成：

g₁(μ)＝X₁β₁

g₂(σ)＝X₂β₂

g₃(ν)＝X₃β₃

当主要研究随机变量分布参数的变化与时间t的关系时，解释变量矩阵就可以表示为：

则可以得到分布参数与时间变量的函数关系：

GAMLSS模型关于回归参数β的似然函数为：

以似然函数取值最大为目标函数，可以采用RS算法(Rigby and Stasinopoulosalgorithm)来估计回归参数β的最优值。GAMLSS模型的全局拟合偏差GD定义为：

式中

为回归参数估计值所对应的对数似然函数。同时引入广义AIC准则(Generalized Akaike Information Criterion，GAIC)进行判断，防止模型的过度拟合，其定义为：

GAIC＝GD+#df

式中：df为模型中的整体自由度，#是惩罚因子。GAIC值最小的模型就作为最优模型。

本发明中的计算将由基于R软件平台的GAMLSS程序包完成。GAMLSS程序包中提供了大量的概率分布函数类型，由于水文计算中常用的样本统计参数主要有均值、均方差(或变差系数)和偏态系数，通过这三种统计参数可以比较确切地描述水文系列的分布特征。因此考虑了四种广泛使用的双参数概率分布(Gumbel、Gamma、对数正态和Weibull分布)来对选定的年最大径流序列数据进行水文频率分析。其中时间t作为协变量以解释时间序列的非平稳性。

设计洪水可采用流量或雨量资料进行推算，但目前水文站网密度比较小，很多河流上没有足够数量和足够观测时间的水文站。中小型水利工程数量繁多，分布广泛，在计算设计洪水时所在流域往往缺乏足够的、甚至全无暴雨洪水资料。因此，研究无资料地区设计洪水计算方法具有重要意义。目前国际上广泛应用的区域洪水频率分析的主要目的是运用区域分析技术解决水文资料短缺地区资料移用问题。其基本思想是，在推求研究测站的某种估计量时，通过考虑某种区域总体的平均特征，改善估计量的性能。

标度洪水方法是一种常用的区域频率分析方法。假定有N个水文测站，第i个站点的年最大洪峰流量资料长度记为n_i，观测的洪峰系列记为Q_ij,j＝1,…,n_i；i＝1,…,N。若这些站点除了洪水尺度系数不同之外洪水频率分布的线型和参数是完全一样的，这些站点就组成了一个所谓的“均匀水文分区”。记Q_i(P)为第i个站点处超过概率为P(0<P<1)的设计洪水,根据标度洪水法，有：

Q_i(p)＝μ_iq(p)

式中μ_i为第i个站点标度洪水；q(P)为区域频率曲线。标度洪水μ_i取值一般为相应站点的洪水频率分布曲线期望值E[Q_i]。区域频率曲线的线型一般可以事先确定，但是其分布参数待定。这些分布参数可以是根据常规矩法或者概率权重矩法求得的变差系数CV和偏态系数CS，也可以是根据线性矩法求得的线性变差系数τ和偏态系数τ₃。

首先对每个站点的洪峰系列无因次化，生成系列：

q_ij＝Q_ij/Q_i(j＝1,...,n_i；i＝1,...,N)

对所有站点的参数估计值进行加权平均,得到区域频率曲线的参数估计值：

即得到估计的区域频率曲线：

对任何一个站点，超过概率P所对应的设计洪峰流量为：

步骤五：长期平稳性分析

根据赫斯特指数估计径流数据的长期平稳性特征。如果H＝0.5，表明时间序列可以用随机游走来描述；如果0.5<H<1，表明时间序列存在长期记忆性；如果0≤H<0.5，表明粉红噪声(反持续性)即均值回复过程。也就是说，只要H≠0.5，就可以用分形布朗运动来描述该时间序列数据。

目前国际上已经开发了几种概念算法来检测径流序列的长期平稳特征。英国科学家赫斯特(赫斯特指数)对尼罗河进行长期的水文观测，采用的数据分析方法，称为变标度极差分析法(简称R/S分析法)。通过分析认为各年的流量存在着一定的时间相关性。赫斯特指数显示出洪涝干旱具有变化的长程效应。其中被广泛应用的是赫斯特开发的重新调整范围(R/S)统计算法、聚合方差方法和去趋势波动分析。

聚合方差法如下所示。对于具有标准偏差σ的离散时间步长为i的平稳过程y_i：

Y_i ^(k)＝(Y_i+...+Y_i-k+1)/k

表示时间尺度k处的聚合过程，标准差为σ^(k)。对于足够大的k，Y_i ^(k)代表综合的水文气候过程；k＝30通常用于根据基本缩放属性估计H：

以上方法弥补了传统长期径流序列变异点分析、趋势分析、水文频率分析中缺乏考虑大规模气候异常现象的不足；对流域历史极端降水的变化进行了全面调查，并对基于流域年最大径流序列的所有统计参数进行线性以及非线性趋势变化分析，以弥补常规方法在水文序列趋势分析中的不足。此外，通过区域洪水频率分析解决了水文资料短缺地区资料移用的问题。

此实施例仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (6)

1.一种基于大规模气候异常的极端降水趋势分析方法，其特征在于，包括以下步骤：

步骤1：气候水文数据收集及预处理；收集气象水文观测数据和大尺度气候异常指标数据，再整理成完整的气象水文信息资料；

步骤2：变异点分析；通过统计检验方法对年最大径流序列方差进行测试，识别年最大径流序列数据中的变异点，并通过调查实地气候水文现象，确定变异点位置；

步骤3：时间趋势分析；利用非参数Mann-Kendall检验、Spearman秩相关检验以及Pearson检验确定年最大径流时间序列是否具有统计上的显著趋势；

步骤4：水文频率分析；利用GAMLSS模型对最大径流序列的统计参数进行线性以及非线性趋势变化分析；

步骤5：长期平稳性分析；利用赫斯特指数H估计径流时间序列的长期平稳性特征。

2.根据权利要求1所述的基于大规模气候异常的极端降水趋势分析方法，其特征在于，所述步骤1中的大尺度气候异常指标数据包括南方涛动指数、北大西洋涛动指数、太平洋十年涛动指数和北太平洋指数。

3.根据权利要求1所述的基于大规模气候异常的极端降水趋势分析方法，其特征在于，所述步骤2中的统计检验方法包括Pettitt检验、累积总和检验、W秩和检验、线性回归中的结构变化检验、贝叶斯变异点检验和小波分析。

4.根据权利要求1所述的基于大规模气候异常的极端降水趋势分析方法，其特征在于，所述步骤3中的非参数Mann-Kendall检验具体如下：

假设有n个样本量(x₁,…,x_n)的时间序列，对于所有k,j≤n，且k≠j，x_k和x_j的分布是不同的，计算检验统计量s为：

其中

当n>10时，标准的正态统计变量通过下式计算：

其中，S为正态分布，均值为0；方差var(s)＝n(n-1)(2n+5)/18；

对于统计值Z，Z大于0时，是增加趋势；Z小于0时，是减少趋势；Z的绝对值在大于1.28、1.64和2.32时，分别表示通过了置信度为90％、95％和99％的显著性检验。

5.根据权利要求1所述的基于大规模气候异常的极端降水趋势分析方法，其特征在于，所述步骤4中的GAMLSS模型具体如下：

假设某一时刻t(t＝1,2,…,n)的随机变量观测值y_t服从概率密度函数f(y_t|θ^t)，θ＝(θ_t1,θ_t2,…,θ_tp)是时刻t对应的分布/统计参数向量，p是分布参数的个数，n是观测值的个数；记y＝(y₁,y₂,…,y_n)^T为独立观测值组成的向量，θ_k为所有时刻的第k个分布参数组成的向量，θ_k＝(θ_1k,θ_2k,…,θ_nk)^T，k＝1,2,…,p；记gk(·)表示θ_k与相应的解释变量X_k和随机效应项之间的单调函数关系，表示为：

式中：η_k为长度为n的向量，β_k＝(β_1k,β_2k,…,β_Ikk)^T为长度为I_k的回归参数向量，X_k为n×I_k的解释变量矩阵，Z_jk为一个已知的n×q_jk固定设计矩阵，γ_jk是一个q_jk维的正态分布随机变量向量，Z_jkγ_jk表示第j项随机效应，q_jk表示第j项随机效应中的随机影响因子维数；前两个参数θ₁和θ₂定义为位置参数向量和尺度参数向量；用v和τ表示随机变量序列的偏度向量和峰度向量；

如果不考虑随机效应对分布参数的影响，即对于k＝1,2,…,p时，令J_k＝0，那么GAMLSS模型为全参数模型：

g_k(θ_k)＝η_k＝X_kβ_k

如果假定随机变量Y服从三参数概率分布，那么GAMLSS模型变成：

g₁(μ)＝X₁β₁

g₂(σ)＝X₂β₂

g₃(ν)＝X₃β₃

当研究随机变量分布参数的变化与时间t的关系时，解释变量矩阵为：

得到分布参数与时间变量的函数关系：

GAMLSS模型关于回归参数β的似然函数为：

以似然函数取值最大为目标函数，采用RS算法来估计回归参数β的最优值；GAMLSS模型的全局拟合偏差GD定义为：

式中

为回归参数估计值所对应的对数似然函数；同时引入广义AIC准则进行判断，防止模型的过度拟合，其定义为：

GAIC＝GD+#df

式中：df为模型中的整体自由度，#是惩罚因子；GAIC值最小的模型作为最优模型。

6.根据权利要求1所述的基于大规模气候异常的极端降水趋势分析方法，其特征在于，所述步骤5中利用赫斯特指数H估计径流时间序列的具体方法如下：

若H＝0.5，径流时间序列用随机游走来描述；若0.5<H<1，径流时间序列存在长期记忆性；若0≤H<0.5，径流时间序列为反持续性。

Images