本发明涉及一种径流时间序列模拟方法，具体涉及一种月径流模拟方法。
背景技术：
对径流时间序列的模拟是水资源管理方案制定的基础，不仅反映水文系统的一些内在规律，更能用于洪水及干旱分析、水文及水资源模拟等，对水资源管理方案制定至关重要的支撑工作。同时作为水资源系统的输入数据，尤其作为水库群的来水数据，径流模拟的随机性也是水资源系统风险的主要来源。更加精确的径流模拟及预测对水资源的规划及管理，特别对于水库群的联合调度具有巨大的现实意义。
对于时间序列模拟的关键在于对历史数据相依结构及统计参数的模拟，早期模型包括ARIMA有参模型及基于自举法的无参模型。无参模型虽然避免了繁杂的参数求解过程，但只能模拟现有径流数据包含的信息，对数据的完整性要求较高。ARIMA模型使用更加广泛，但其基于正态边际分布及线性等假设限制了其精确度。Copula是一种新型的求解多维联合分布的函数，其灵活的结构可避免传统模型的缺点，在径流模拟中有着巨大的潜力。而现有研究多局限在二维copula，结构更加灵活功能更加强大的多元copula少有研究。
技术实现要素：
发明目的：本发明的目的在于针对现有技术的缺点，提供一种基于D藤copula函数的月径流模拟方法，适用于当前气候变化条件下适用于偏态分布的季节性月径流模拟模型，该模型为流域水资源统一调度提供更为有力的依据。
技术方案：本发明提供了一种基于D藤copula函数的月径流模拟方法，包括以下步骤：
(1)采用图形方法构建月径流密度分布函数拟合模型；
(2)采用混合D藤copula函数构建各月月径流相依结构模拟模型；
(3)采用基于Kendall系数的相关度检验法确定各月径流混合D藤copula维度；
(4)采用基于极大似然法和AIC法的有序法推求D藤copula参数；
(5)采用条件函数构建基于条件混合D藤copula联合分布函数的月径流模拟模型。
进一步，步骤(1)选取三种或以上接近历史数据分布特点的单变量分布函数作为备选函数，分别将历史径流序列拟合各备选函数，并绘制累计分布函数图检验拟合效果，选出最佳分布函数。采用极大似然法求解边际分布函数参数如下：
式中，Θ为边际函数的参数向量，L为分布函数的似然值。
进一步，气候变化条件下月径流时间序列具有季节性、偏态分布等特点，步骤(2)不同月份的径流采用不同结构的混合D藤copula模拟其与之相关月径流间的相依结构。D藤copula中各树由一系列二维copula以及相应条件copula组成，混合D藤copula中这些基础结构copula由不同种类的二维copula组成，最大程度上提供更加灵活的结构，以模拟不同的相依类型。模型将从Elliptical copula、Archimedean copula、二参数Archimedean copula以及Clayton copula、Gumbel copula的旋转类型中选择需确定的参数，有混合D藤copula的维度、变量以及组成其树结构的二维copula类型。
假设t月径流与其前d-1个连续月的径流相关，其联合分布用混合D藤copula函数表示为：
式中，t＝1,2,...,12，ft代表联合分布密度函数，Xt代表t月径流，xt代表相应取值，i代表连线，j代表树，k代表月径流变量，代表相应下角标的D藤copula的密度函数，可表示为：
式中，Fd(xd)为边际分布函数，fd(xd)为密度函数；
以四维模型为例，四个月径流的联合分布可表示为：
式中，F(·)为月径流边际分布，F·|·(·)为月径流条件分布。
进一步，维度，即与当前月径流相关的月径流变量数目，是构建各月联合分布模型的重要参数。步骤(3)采用Kendall系数衡量各月径流间的相依关系，从而挑选相关径流作为预测因子，D藤copula的变量即为预测因子及当月径流，其维度为变量个数；
两变量间的Kendall系数计算如下：
式中，τ为Kendall’s tau值，n为月径流观测个数，Pn及Qn分别表示一致的径流数对和不一致的径流数对数量；
对于不同两月径流变量X1与X2无相关关系的假定，τn服从期望为0，方差为的正态分布；则，当取置信度为α＝5％时，
假定不成立，即两变量相关。
进一步，步骤(4)采用顺序法从上至下顺次确定各树T1,...,Td-1中结构copula类型及其参数，依据误差函数挑选copula类型后估计参数；
采用赤池信息量准则作为误差函数进行误差分析，选出组成D藤copula的最佳二维Copula函数系列：
AIC＝-2log(L(θ|data))+2V (6)
式中，L为密度函数的极大似然值，θ为参数向量，V为参数数目，data为历史径流数据，log(L(θ|data))表示为：
式中，m、n代表月径流观测个数和总数，c(·,·)为二维copula密度函数，u1,m为月径流变量X1第m个观测值的边际分布，u2,m为月径流变量X2第m个观测值的边际分布；
确定某基础copula类型后采用极大似然法计算参数值，如下：
式中，为似然函数最大时参数值。
进一步，步骤(5)已知当月径流及与之相关的月径流联合分布函数，当月径流由条件联合分布函数的反函数及与当月径流相关的月径流值获得：
假设当月径流为xt，其与另外d-1个月的径流相关，则xt可由条件函数的反函数求得：
xt＝F-1(wt|xt-d+1,...,xt-1) (9)
其中
式中，wt为服从(0，1)均匀分布随机数，F-1(·|·)是F(·|·)的反函数，Ct,t-d+1|t-d+2,...,t-1代表相应下角标的D藤copula的分布函数。
有益效果：本发明构建了基于混合D藤copula的月径流模拟模型，该模型依据Kendall’s tau值构建各月的相关结构，置信值在95％以上即为相关，采用混合D藤copula计算相关结构，克服了传统模型线性和稳态等问题，模拟结果更准确可靠。
附图说明
图1为本发明的流程图；
图2为案例流域图；
图3为唐乃亥站点月径流密度分布函数估计图；
图4为唐乃亥站点各月径流|Z|值图；
图5为黄河唐乃亥站点11月月径流模型数结构图；
图6为黄河唐乃亥站点月径流模拟结果图。
具体实施方式
下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
一种月径流模拟方法，如图1所示，本案例选取黄河干流唐乃亥站点的月径流进行时间序列模拟，黄河流域图及该站地理信息如图2所示，利用该站站点1956年至2010年的还原月径流资料：
(1)构建各月径流边际分布拟合模型
根据唐乃亥月径流历史数据分布特点，选取Weibull函数、lognormal函数以及gamma函数分布作为备选函数，分别对各月径流序列进行函数拟合，绘制累计分布函数图如图3所示，经检验，gamma函数为最佳分布函数，并按式(1)进行参数计算，结果如下表：
(2)采用基于Kendall系数的相关度检验法确定各月径流模拟模型维度
根据式(4)及式(5)依次计算各月径流间的Kendall系数及相应|Z|值，探究相关关系，选取各月径流相关的月径流，如图4所示，图中横坐标Lag为滞后月，纵坐标为Z值。
(3)采用基于极大似然法和AIC法的有序法推求模型参数
根据式(6)及式(8)顺次推求树T1中各组成copula的种类及其相应参数，在根据树T1的结构计算其条件copula作为树T2的变量，依次按上述步骤计算至树Td-1各基础copula的种类及其相应参数。11月径流的混合D藤copula结构如图5所示。统计各月混合D藤copula中各基础copula种类如下表所示：
(4)构建基于条件混合D藤copula的月径流模拟模型
随机生成服从(0，1)均匀分布随机数wj，并已知当月径流xt相关的在其前发生的月径流值xt-d+1,…,xt-1，可根据式式(9)及式(10)推求xt，重复以上步骤顺次推求其余月径流，获得月径流时间序列。生成500组与历史数据等长的月径流序列如图6所示，模拟序列完全涵盖历史序列。