1. 引言
中长期径流预报作为水文预报的重要组成部分，对洪水调控，兴利调度和水资源系统规划管理有着重要意义。近些年来，随着科学理论水平的不断提高，诸如人工神经网络(ANN)，支持向量机(SVM)，遗传规划(GP)等智能算法也广泛运用到中长期径流预报领域 [1] [2] 。朱双等提出基于GM-SVM耦合的月径流预测模型，发现对径流变化剧烈的汛期表现出更优越的预测性能 [3] ，巴欢欢等将SSA-LSSVM等方法用于水布垭水库，发现经过数据前处理，径流预报精度得到提高 [4] 。然而，中长期径流由于本身水文过程复杂，受到诸多不确定因素影响，单一模型虽然有着各自的独立优势，然而其预报效果始终有限。概率预报方面，刘章君等提出贝叶斯概率洪水预报模型，所得预报结果略优于现有的HUP模型和BP-BFS模型，但该模型假定流量过程服从一阶马尔科夫过程 [5] ；李彩林等基于相似过程衍生的月径流概率预报模型进行研究，针对汛期，过渡期，非汛期分别进行误差统计 [6] 。本文根据SSA-ANN和SSA-SVM确定性模型的预报值，通过copula函数建立多变量联合分布，给出实测流量的条件概率分布，确立可靠的概率预报区间，不仅能发布确定的预报结果，还能给出一定置信系数下未来径流变化范围，从而更利于水库月调度计划的制定。
2. 研究方法
2.1. 奇异谱分析
奇异谱分析(Singular Spectrum Analysis, SSA)是一种广义功率谱。它根据所观测到的时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解、重构，从而提取出代表原时间序列不同成分的信号，如长期趋势信号、周期信号、噪声信号等，从而对时间序列的结构进行分析，并可进一步预测 [7] 。
2.1.1. 建立相空间
给定一个非零的时间序列，采用动力系统分维数估计的处理方法，按时迟滞排列建立相空间，将一维时间序列转化为多维时间序列：
(1)
式中：X称为相空间中的轨迹矩阵，
称为窗口长度。
2.1.2. 奇异值变换
计算XXT并求得其L个特征值及其相应特征向量。通过分组和重构成分，提取
个具有贡献作用的有用成分进行构建。该重建序列是被SSA滤去噪声后的有用序列，蕴含原序列的周期成分；原始序列与重建序列的差值即为噪声 [8] 。
2.2. 人工神经网络模型
人工神经网络(Artificial Neural Network, ANN)广泛运用于科学，工程等诸多领域 [9] 。人工神经网络作为由许多非线性且紧密联系的神经元构成的信息处理系统，表现出并行性，非线性映射能力，鲁棒性和容错性，自学习和自适应等特点。通常，一个神经网络可以看作三层模型，即输入层，隐含层和输出层，层与层之间通过权重连接。其中，输入层作为数据输入层，隐含层作为数据处理层，输出层给出数据处理后结果。
三层前向神经网络由于结构简单，又能解决较多实际问题，使得这种模型如今仍受到诸多使用。以BP网络为例，网络学习可以看作两部分：第一步是前向学习，输入训练样本，通过构架好的网络结构和最近一次迭代得到的权值和阈值，依次从前向后计算输出层各神经元结果。第二步反馈修正，利用输出数据对层间权值和神经元阈值进行修改，从后往前计算它们对总误差的影响梯度，比较预先设定精度，据此判断模型学习是否达到标准。神经网络模型较传统模型的优点在于，无需用显式数学表达复杂的水文物理成因过程。
2.3. 支持向量机模型
支持向量机(Support Vector Machines, SVM)是基于结构风险最小化原则，将最优分类问题转化为求解凸二次规划问题，得到全局最优解，较好地解决局部极小值的问题，同时在一定程度上克服了“维数灾”和“过学习”等传统困难，因此在诸多领域广泛应用 [10] 。给定一组训练集
(
作为输入向量，
是期望值)，SVM的回归方程可由下式给出：
(2)
式中：
为超平面的权值向量，b为偏置项。
根据支持向量机理论，可将线性回归方程转化为优化问题求解，该约束优化问题可用拉格朗日形式求解。方程满足Mercer条件 [11] ，故引入核函数
。核函数反映的是高维特征空间的内积，实现了从低维空间不可分到高维空间可分的功能。几种常用的核函数包括线性核函数、多项式核函数、Gauss径向基核函数、Fourier核函数等。采用以σ为参数的Gauss径向基函数(RBF)，其表达形式为：
(3)
2.4. 基于Copula函数的径流概率预报模型
2.4.1. Copula函数
Sklar提出可以将一个m维联合分布函数分解为m个边缘分布函数和一个Copula函数。Nelsen于1999年给出了copula函数的严格定义，即是把随机变量
的m维联合分布函数
分布函数与各自的边缘分布
相连接的函数。即存在一个m-Copula函数C，使得对任意
[12] ：
(4)
式中：
为Copula函数的相关性参数。
为评价模拟效果，引用综合评价指标均方根误差RMSE，定义为
(5)
式中：
和
分别为经验频率与理论频率；
为资料系列长度。RMS值越小，表明经验频率与理论频率之间拟合越好。
2.4.2. 条件概率
令
，
分别代表ANN和SVM确定性预报流量，
表示待预报的实际流量；
、
、
分别为
、
、
的取值。令
的边缘分布函数为
，相应的概率密度函数为
。根据Copula函数得到
和
的联合分布
，
为相应概率密度函数；
、
和
的联合分布可写为
(6)
给定确定性模型预报结果
、
后，
的条件分布函数可表示为：
(7)
对h求导得条件密度函数：
(8)
其中，
，
分别为三维Copula，实测流量的概率密度函数。对式(8)进行复化梯形公式计算，得到中值确定性预报及相应预报区间。
2.5. 模型性能指标
为评判模型的模拟预测能力，采用《水文情报预报规范》推荐的四个指标：即1) 纳什效率系数(NS)，2) 水量平衡系数(WB)，3) 年均最大径流的相对误差(Remax)，4) 相关性系数®。
3. 实例研究
丹江口水库位于湖北省丹江口市，汉江与其支流丹江汇合口下游800 m处，多年平均入库水量为394.8亿m3。开展丹江口水库的径流中长期预报，不仅是确定南水北调中线可调水量的重要依据，也是制定丹江口水库防洪、供水、发电等优化调度方案的重要参考；同时中长期径流预报的预见期较长，能为水库管理者提供更多统筹规划的时间，使水库尽可能获得最大的经济社会效益。但也正是因为较长的预见期，导致预报精度不高和预测结果的不确定性较大，因此开展丹江口水库中长期径流预报的研究对丹江口水库各项功能的充分实现具有重要意义。本文采用1954~2013年丹江口水库的月径流数据，开展概率预报研究。
3.1. 数据预处理
3.1.1. SSA参数的选择
在进行ANN和SVM确定性预报之前，先对丹江口1954~2013年的月径流序列进行SSA降噪处理。由章节1.1的描述可知，SSA中需要确定两个重要参数：窗口长度L和贡献成分p。目前研究中，窗口长度的选择没有具体规则，但可采用基于均方误差最小原则的窗口长度选取方法。由于月径流本身具有季节性特点，故L值最大不超过12，同时，窗口长度必须大于1，则窗口长度L取值范围为[2,12]之间的整数值。通过计算ANN模型和SVM模型在窗口长度L取[2,12]时，在检验期内的预测值与实测值的均方误差，结果表明两者都是在L = 11时均方误差达到最小。因此运用奇异谱分析对数据进行预处理时，ANN模型和SVM模型的最优窗口长度均取11。
确定窗口长度L后，原始序列被分解为L个子序列。通过互相关函数法(CCF)确定原序列有用成分。表1
Table 1. Cross-correlation function values between subseries and original series for different L
表1. 不同L值子序列与原序列的互相关系数统计表
为各个窗口长度下子序列与原始序列的互相关系数统计值。从表中可知，子序列与原序列互相关系数为正，即为贡献成分；反之，互相关系数为负，即为噪声序列；依此类推，即可得到各个窗口长度下的贡献成分的个数及相应的重组序列和噪声，将贡献成分叠加构成重建序列RC。
3.1.2. 预报因子的选取
分别选取ANN和SVM模型模拟预测丹江口水库的月径流。在利用数据驱动模型预报之前，选择合适的模型输入因子尤为重要。在筛选ANN模型和SVM模型的过程中，引入两种统计性指标(自相关系数ACF和偏相关系数PACF)，比较不同滞时(1~24)的径流序列的相关函数值，选取对应较大自相关系数值的滞时参数。图1给出了不同滞时情况下，窗口长度为L = 11时重组序列的自相关函数和偏相关函数值。从图中可知，当滞时为1时，序列的自相关系数值最大，当滞时为12时，自相关系数值次之。然而，对于智能算法预测模型而言，预测因子越多，模型所能捕捉的序列信息也就越多，同时由偏相关函数值看出，滞时为12时，其达到95%的置信区间。故为了尽可能地保证模型精度，选取前12个月的径流数据作为ANN和SVM模型重构序列的预报因子。
3.2.实测流量的条件概率分布
3.2.1. 边缘分布的确定
边缘分布函数可以采用任意分布，选用的标准是使得假定理论分布与经验分布的标准差最小。我国常采用P-III分布作为边缘分布，其概率密度函数为：
(9)
式中：
、
和
分别为P-III型分布的形状、尺度和位置参数；
为伽马函数。
，
，
三个变量的P-III分布参数采用适线法估计，并采用Kolmogorov-Smirnov (K-S)检验法对这些分布拟合进行检验，结果见表2。在5%的显著性水平下(临界值为0.0511)，各变量均通过了检验。
3.2.2. 联合分布的建立
利用ANN与SVM预报值
，
，建立二维copula联合分布。通过计算计算二维Copula函数的RMSE值，
Figure 1. The autocorrelation and partial autocorrelation function for reconstructed series under L = 11
图1. L = 11重构序列的自相关和偏相关函数值
Table 2. Estimated parameters of marginal distribution and K-S test
表2. 边缘分布参数估计结果和K-S检验
且由于
、
与
之间存在正相关关系，最终选用Archimedean Copula函数族中的GumbelCopula 函数构造联合分布，其秩相关系数
为0.85，根据
与
的关系，得到copula参数为6.72。利用
、
与
构建三维联合分布，对其对称参数，非对称参数进行极大似然法估参，比较各RMSE值，最终采用非对称GumbelCopula函数作为变量的联合分布。
3.2.3. 概率预报
考虑P-III分布样本空间
，同时考虑到丹江口水库月径流系列的特点，取
，
。将得到的H取值区间[70,15000] m3/s以
为间距离散为149300等分。
则H的第
个取值为
(10)
将获得的边缘分布和联合分布的参数结果代入式(6)~(8)，当已知
时，可得到关于流量的条件概率密度曲线
。
3.3. 结果分析
选取丹江口水库1954~2013年月径流数据进行概率预报研究，在选用ANN和SVM智能模型预报前，先对资料系列数据进行处理。采用1954~2002年资料进行模型训练，对模型训练数据与实测数据建立联合分布关系，采用2003~2013年数据进行确定性预报和概率性预报检验。其中，概率性预报以中值形式给出确定性预报结果，其模拟预测结果见表3。
Table 3. Performance indices of models during training and testing periods
表3. 模型率定期和检验期各模型的模拟结果
Figure 2. Comparison of observed and probability forecasted monthly inflow and 90% confidence intervals during testing period
图2. 检验期月径流实测值与概率预报与90%置信区间对比图
由表3得出，SSA-ANN和 SSA-SVM月径流预报模型率定期和检验期纳什效率系数在80%左右。相较确定性预报而言，概率性预报的效率系数在一定范围内有所提高。水量平衡系数也接近于1，表明模型模拟预测的总水量小于实际总水量，可能原因在于洪峰的预测偏小；年均最大径流相对误差在15%左右波动；相关性系数达到或接近90%。
为了直观分析确定性模型，概率模型预报效果，图2给出了检验期确定性预报与概率预报中值及实测流量过程对比。由图2可以直观地看出，预测径流与实测径流的变化规律模拟的基本一致，但预测值在年最大径流处的预测效果较差，可能原因是SSA将径流极大值处的波动信号视为噪声信号，将其过滤，从而导致极值变化规律在模型中丢失。另外，根据数理统计原理，给定显著性水平α = 0.1，计算条件概率5%和95%的分位数，它们分别代表着90%的流量预报区间的置信下限和上限值。检验期内，实测流量出现峰值时，由于确定性预报模型SSA-ANN和SSA-SVM值偏小，导致概率预报无法包含实测流量，有待进一步研究。但实测流量绝大部分位于预报区间之内，表明概率区间预报基本可靠，可以为防洪决策提供更多的信息，使得预报人员在决策中能定量考虑各种不确定性，实现预报与决策的有机结合。
4. 结论
水库来水量的预报，牵涉到地质、气候条件等因素影响，对于现有方法和模型，还是有不完善之处预报结果不完全一致给决策工作带来困难。对多个模型预报结果分析，通过选取合适的copula函数，构造多变量联合分布，进一步给出概率预报区间以及中值预报结果。通过模型分析比较，得到以下结论：
1) 概率预报可以在一定程度上整合不同模型的预测优势，提高预报精度。与单个模型方法相比，概率预报模型具有更大的优势。
2) 基于Copula方法的概率预报可以给出置信水平下的预报区间，有利于决策人员定量考虑预报的不确定性。
基金项目
国家自然科学基金(51539009)和十三五国家重点研发(2016YFC0402206)资助项目。