一种基于条件降维重构的日径流季节性随机模拟方法与流程

本发明属于随机水文领域，更具体地，涉及一种基于条件降维重构的日径流季节性随机模拟方法。

背景技术：

径流序列作为水资源系统的输入而在系统模拟中起着至关重要的作用，常用于制定水库调度策略、评估水资源系统的风险、确定各种水力参数等。由于实测水文序列常只有几十年甚至十几年，序列较短，难以满足水资源系统分析和风险评估需要。因此，径流随机模拟方法常用于生成长序列水文数据。

常用的随机模拟方法主要是一阶自回归方法(ar(1))，该模型结构简单，概念清晰，易于实现，在水文水资源随机模拟中得到了广泛应用。但自回归模型假定实测径流数据服从正态分布且具有线性相关性，实践证明水文数据存在偏态特征；因此，在模拟偏态水文序列时，首先要进行正态变换，此过程难免会造成部分信息失真。为解决此问题，有学者提出了基于copula函数的径流随机模拟方法(car(1))，该方法有效地解决了回归模型对径流序列分布形式的正态假定问题。然而，目前基于copula函数的随机模拟方法仅考虑了径流序列的1阶滞时相关性特征，不能够充分利用实测径流信息，使得模拟精度仍有提升空间。

技术实现要素：

针对现有技术的缺陷，本发明的目的在于解决现有技术直接建立高维copula函数带来的参数求解困难、模拟精度有待提高的技术问题。

为实现上述目的，第一方面，本发明实施例提供了一种基于条件降维重构的日径流季节性随机模拟方法，该方法包括以下步骤：

s1.根据历年逐日实测日流量序列数据，拟合对应的日流量边缘分布；

s2.根据日流量边缘分布和二维copula函数，分别建立t与t-1时刻日流量、t-1与t-2时刻日流量及t与t-2时刻日流量的二维联合分布；

s3.根据t-1与t-2时刻日流量的二维联合分布、t与t-2时刻日流量的二维联合分布，建立t、t-1、t-2时刻日流量的三维联合分布；

s4.根据t、t-1、t-2时刻日流量的三维联合分布，计算已知t-1时刻日流量值和t-2时刻日流量值条件下t时刻日流量的条件分布；

s5.根据t-1时刻日流量值得到t-1时刻日流量边缘分布值，根据t-2时刻日流量值得到t-2时刻日流量边缘分布值，根据t-1、t-2时刻日流量边缘分布值和已知t-2时刻日流量值条件下t-1时刻日流量的条件分布，得到已知t-2时刻日流量值条件下t-1时刻日流量的条件分布值；

s6.令已知t-1、t-2时刻日流量值条件下t时刻日流量的条件分布值为0～1之间的随机数，根据该随机数、已知t-2时刻日流量值条件下t-1时刻日流量的条件分布值和已知t-1时刻日流量值和t-2时刻日流量值条件下t时刻日流量的条件分布，得到已知t-2时刻日流量条件下t时刻日流量的分布值；

s7.根据已知t-2时刻日流量条件下t时刻日流量的条件分布值和已知t-2时刻日流量值条件下t时刻日流量的条件分布，得到t时刻日流量的边缘分布值，根据t时刻日流量的边缘分布值和t时刻日流量的边缘分布，得到t时刻日流量模拟值。

具体地，步骤s1中，t时刻对应的日流量边缘分布f(qt)计算公式如下：

其中，fp(·)为所述边缘分布的概率密度函数，qt为t时刻日流量qt的取值。

具体地，二维copula函数c计算公式如下：

c(u1，u2)＝exp(-[(-(lnu1))^θ+(-(lnu2))^θ]^1/θ)

其中，θ为二维copula函数参数。

具体地，步骤s2包括以下步骤：

s201.分别建立t与t-1、t-1与t-2及t与t-2时刻日流量的二维联合分布；

s202.根据t、t-1、t-2时刻的日流量边缘分布，通过计算kendall秩相关系数和kendalltau偏相关系数，确定各二维联合分布中copula函数参数θ。

具体地，步骤s201具体为：

f(qt-1，qt)＝c(f(qt-1)，f(qt))＝c(u1，u3)

f(qt-2，qt)＝c(f(qt-2)，f(qt))＝c(u2，u3)

f(qt-1，qt-2)＝c(f(qt-1)，f(qt-2))＝c(u1，u2)

其中，f(qt-1，qt)表示t与t-1的时刻日流量的二维联合分布，c(，)为二维copula函数，u1＝f(qt-1)表示t-1时刻日流量的边缘分布，u2＝f(qt-2)表示t-2时刻日流量的边缘分布，u3＝f(qt)表示t时刻日流量的边缘分布。

具体地，步骤s202中，日流量qi和qj之间的kendall秩相关系数τqiqj计算公式如下：

t与t-2时刻日流量之间的kendalltau偏相关系数计算公式如下：

其中，为在已知变量qt-1的条件下qt和qt-2之间的偏相关系数，为qt和qt-2之间的kendall秩相关系数，为qt和qt-1之间的kendall秩相关系数，为qt-1和qt-2之间的kendall秩相关系数。

具体地，步骤s202中，二维copula函数c(u1，u3)或c(u1，u2)的参数θ，根据θ与kendall秩相关系数τ的关系间接得出；

二维copula函数c(u2，u3)的参数θ，根据θ与kendalltau偏相关系数的关系间接得出；

具体地，步骤s3包括以下步骤：

s301.根据二维联合分布，计算已知qt-2＝qt-2条件下qt的条件分布和已知qt-2＝qt-2条件下qt-1的条件分布；

其中，c(u3|u1)为已知u1的条件copula函数，c(u1|u2)为已知u2的条件copula函数；

s302.根据t-1与t-2时刻日流量的二维联合分布、t与t-2时刻日流量的二维联合分布，建立t、t-1、t-2时刻日流量的三维联合分布；

其中，c(，，)为三维copula函数，w1＝c(u1|u2)，w2＝c(u3|u2)，u1表示t-1时刻日流量的边缘分布，u2表示t-2时刻日流量的边缘分布，u3表示t时刻日流量的边缘分布。

具体地，步骤s4包括以下步骤：

s401.根据t、t-1、t-2时刻日流量的三维联合分布，计算copula函数的密度函数；

其中，c(，，)为三维copula函数，c(，)为二维copula函数，u1表示t-1时刻日流量的边缘分布，u2表示t-2时刻日流量的边缘分布，u3表示t时刻日流量的边缘分布；

s402.基于copula函数的密度函数，计算已知qt-2＝qt-2和qt-1＝qt-1条件下qt的条件分布；

其中，c(·)为copula函数的密度函数，h为条件分布。

第二方面，本发明实施例提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现上述第一方面所述的基于条件降维重构的日径流季节性随机模拟方法。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：

1.本发明采用条件降维重构理论构造三维copula函数，将三维copula函数转化为条件分布和二维copula函数，避免了直接建立高维copula函数带来的参数求解困难，降低了构造高维copula的难度同时比较容易估计高维copula的参数，使计算变得更容易。

2.本发明在日径流随机模拟中考虑2阶滞时相关关系，解决了现有基于copula函数的日径流随机模拟方法无法考虑2阶滞时的难题，能够很好的模拟日径流序列均值、方差、cs等统计特征，且模拟效果更好。

附图说明

图1为本发明实施例提供的一种基于条件降维重构的日径流季节性随机模拟方法流程图；

图2为本发明实施例提供的采用本发明方法与采用二维copula模拟的效果对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

事实上，径流序列多个滞时间仍具有较强相关关系，基于此，本发明提出一种基于条件降维重构的日径流季节性随机模拟方法。

如图1所示，一种基于条件降维重构的日径流季节性随机模拟方法，该方法包括以下步骤：

s1.根据历年逐日实测日流量序列数据，拟合对应的日流量边缘分布；

s2.根据日流量边缘分布和二维copula函数，分别建立t与t-1时刻日流量、t-1与t-2时刻日流量及t与t-2时刻日流量的二维联合分布；

s3.根据t-1与t-2时刻日流量的二维联合分布、t与t-2时刻日流量的二维联合分布，建立t、t-1、t-2时刻日流量的三维联合分布；

s4.根据t、t-1、t-2时刻日流量的三维联合分布，计算已知t-1时刻日流量值和t-2时刻日流量值条件下t时刻日流量的条件分布；

s5.根据t-1时刻日流量值得到t-1时刻日流量边缘分布值，根据t-2时刻日流量值得到t-2时刻日流量边缘分布值，根据t-1、t-2时刻日流量边缘分布值和已知t-2时刻日流量值条件下t-1时刻日流量的条件分布，得到已知t-2时刻日流量值条件下t-1时刻日流量的条件分布值；

s6.令已知t-1、t-2时刻日流量值条件下t时刻日流量的条件分布值为0～1之间的随机数，根据该随机数、已知t-2时刻日流量值条件下t-1时刻日流量的条件分布值和已知t-1时刻日流量值和t-2时刻日流量值条件下t时刻日流量的条件分布，得到已知t-2时刻日流量条件下t时刻日流量的分布值；

s7.根据已知t-2时刻日流量条件下t时刻日流量的条件分布值和已知t-2时刻日流量值条件下t时刻日流量的条件分布，得到t时刻日流量的边缘分布值，根据t时刻日流量的边缘分布值和t时刻日流量的边缘分布，得到t时刻日流量模拟值。

步骤s1.根据历年逐日实测日流量序列数据，拟合对应的日流量边缘分布。

日流量序列是指多个日流量组成的序列，经过步骤s1之后，得到365个日流量边缘分布。本实施例中，边缘分布为中国水利水电工程水文计算规范推荐的采用的p-iii分布，也可以是其他分布。

p-iii分布的概率密度定义为：

其中，q表示日流量，α、β和δ分别为p-iii分布的形状、位置和尺度参数，α＞0，β＞0，δ≤q＜∞，γ(·)为gamma函数。本发明采用线性矩法估计p-iii分布的参数。

t时刻日流量边缘分布f(qt)计算公式如下：

假设t、t-1、t-2时刻日流量分别表示为qt、qt-1和qt-2，其取值分别为qt、qt-1和qt-2，本发明的目的是在已知日流量值qt-1和qt-2的前提下，计算日流量值qt。

步骤s2.根据日流量边缘分布和二维copula函数，分别建立t与t-1时刻日流量、t-1与t-2时刻日流量及t与t-2时刻日流量的二维联合分布。

本发明针对的是日径流季节性随机模拟，因此，本发明中“时刻”代表日尺度。步骤s2包括以下步骤：

s201.分别建立t与t-1、t-1与t-2及t与t-2时刻日流量的二维联合分布，各二维联合分布中copula函数参数θ未知。

u1＝f(qt-1)

u2＝f(qt-2)

u3＝f(qt)

其中，u1表示t-1时刻日流量的边缘分布，u2表示t-2时刻日流量的边缘分布，u3表示t时刻日流量的边缘分布。已知日流量值qt-1和qt-2，则日流量的边缘分布值u1和u2也已知。

c(u1，u2)＝exp(-[(-(lnu1))^θ+(-(lnu2))^θ]^1/θ)

t与t-1时刻日流量、t-1与t-2时刻日流量及t与t-2时刻日流量的二维联合分布分别表示如下：

f(qt-1，qt)＝c(f(qt-1)，f(qt))＝c(u1，u3)

f(qt-2，qt)＝c(f(qt-2)，f(qt))＝c(u2，u3)

f(qt-1，qt-2)＝c(f(qt-1)，f(qt-2))＝c(u1，u2)

其中，f(qt-1，qt)表示t与t-1的时刻日流量的二维联合分布，c为copula函数。

s202.根据t、t-1、t-2时刻的日流量边缘分布，通过计算kendall秩相关系数和kendalltau偏相关系数，确定各二维联合分布中copula函数参数。

qi和qi之间的kendall秩相关系数可定义为：

此方法能较容易地估计二维copula函数的参数，因此，c(u1，u3)和c(u1，u2)的参数可通过以上关系式估计。但计算t与t-2时刻日流量序列的相关关系时，应移除t-1时刻日流量的影响，因此需要计算偏相关系数，kendalltau偏相关系数的计算表达式为：

其中，为在已知变量qt-1的条件下qt和qt-2之间的偏相关系数，为qt和qt-2之间的kendall秩相关系数，为qt和qt-1之间的kendall秩相关系数，为qt-1和qt-2之间的kendall秩相关系数。

二维copula函数c(u1，u3)或c(u1，u2)的参数θ，可根据θ与kendall秩相关系数τ的关系间接得出。

二维copula函数c(u2，u3)的参数θ，可根据θ与kendalltau偏相关系数的关系间接得出。

步骤s3.根据t-1与t-2时刻日流量的二维联合分布、t与t-2时刻日流量的二维联合分布，建立t、t-1、t-2时刻日流量的三维联合分布。

s301.根据二维联合分布，计算已知qt-2＝qt-2条件下qt的条件分布和已知qt-2＝qt-2条件下qt-1的条件分布。

其中，c(u3|u1)为已知u1的条件copula函数，c(u1|u2)为已知u2的条件copula函数。

s302.根据t-1与t-2时刻日流量的二维联合分布、t与t-2时刻日流量的二维联合分布，建立t、t-1、t-2时刻日流量的三维联合分布。

其中，w1＝c(u1|u2)，w2＝c(u3|u2)。

步骤s4.根据t、t-1、t-2时刻日流量的三维联合分布，计算已知t-1时刻日流量值和t-2时刻日流量值条件下t时刻日流量的条件分布。

s401.根据t、t-1、t-2时刻日流量的三维联合分布，计算copula函数的密度函数。

s402.基于copula函数的密度函数，计算已知qt-2＝qt-2和qt-1＝qt-1条件下qt的条件分布。

其中，c为copula函数的密度函数，h为条件分布。

步骤s5.根据t-1时刻日流量值得到t-1时刻日流量边缘分布值，根据t-2时刻日流量值得到t-2时刻日流量边缘分布值，根据t-1、t-2时刻日流量边缘分布值和已知t-2时刻日流量值条件下t-1时刻日流量的条件分布，得到已知t-2时刻日流量值条件下t-1时刻日流量的条件分布值。

由于t-1和t-2时刻日流量已知，则可拟合其对应的边缘分布，进而求得在已知t-2时刻日流量条件下t-1时刻日流量的条件分布w1＝c(u1|u2)的值。

步骤s6.令已知t-1、t-2时刻日流量值条件下t时刻日流量的条件分布值为0～1之间的随机数，根据该随机数、已知t-2时刻日流量值条件下t-1时刻日流量的条件分布值和已知t-1时刻日流量值和t-2时刻日流量值条件下t时刻日流量的条件分布，得到已知t-2时刻日流量条件下t时刻日流量的分布值。

随机生成服从(0，1)均匀分布随机数ε，令h(w2|w1)＝ε，通过计算反函数h^-1(w2|w1)得到条件分布w2＝c(u3|u2)的值。

步骤s7.根据已知t-2时刻日流量条件下t时刻日流量的条件分布值和已知t-2时刻日流量值条件下t时刻日流量的条件分布，得到t时刻日流量的边缘分布值，根据t时刻日流量的边缘分布值和t时刻日流量的边缘分布，得到t时刻日流量模拟值。

在w2已知的条件下，通过求反函数c^-1(u3|u2)得到t时刻日流量的边缘分布值u3，进而计算反函数f^-1(u3)获得t时刻日流量qt。

实施例：现选取长江上游和珠江流域7个站点的日径流数据进行案例研究，7个站点分别为长江宜昌站、金沙江屏山站、岷江高场站、沱江李家湾站、嘉陵江北碚站、乌江武隆站和西江梧州站，宜昌站实测径流序列长度为1876-2015年，梧州站实测径流序列长度为1962-2008年，其他站点实测径流序列长度为1952-2015年，检验一种基于条件降维重构的日径流季节性随机模拟方法的模拟效果。

表1模拟和实测的统计特征值及其相对误差

根据日模拟和实测日流量的基本统计数据，200个样本的统计特征值如图2所示，从左到右分别是均值、sd、cs三个统计特征值效果对比，可得模拟日流量的统计特征值与实测日流量序列的较为一致。另外，表1给出了基于本发明的7个日流量站模拟和实测日流量序列统计特征值的相对误差。结果表明，本发明能够保留实测径流序列的统计特性。

表2给出了基于传统自回归方法、二维copula函数随机模拟法和本发明方法的对比结果。由表2可知，本发模拟的日流量统计特征值(包括均值、方差、偏态系数)的相对误差均小于其他两种方法。结果表明，本的模拟效果优于其他两种方法。

表2基于sarima、二维copula函数和本发明的模拟日流量统计特征值相对误差

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