MIC-PCA耦合算法 在径流预报因子筛选中的应用

王丽萍，李宁宁，马皓宇，纪昌明，李贵博

(华北电力大学可再生能源学院，北京 102206)

0 引 言

径流预报是根据前期水文气象要素，通过成因分析与数理统计等方法，对未来一段时期的径流进行科学的预测。准确及时的径流预报，对于争取防汛、抗旱的主动权，制定科学的水资源调度方案，确保水利设施的安全和发挥其经济效益具有重要意义[1]。传统径流预报方法一般包括物理成因法和数理统计法等[2,3]。近年来国内外学者从各个方向对径流预报的理论与方法进行了深入的研究与探索，提出了许多基于现代智能方法和数值天气预报的综合预报模型，主要包括模糊分析[4]、灰色系统理论[5]、混沌理论[6]、小波分析理论[7]、人工神经网络等[8]，这些方法对于提高径流预报结果的精度和可靠性有着重要意义。然而，相比于预报模型，目前对水文预报因子选取的研究相对较少，预报因子作为预报模型的输入侧，直接影响预报结果的准确性和短期预报的实效性，如何从众多预报因子中筛选出关于预报对象信息含量高、冗余信息少的预报因子集是当前预报工作中的重点。

预报因子(输入变量)的选取是预报模型应用于实际工程的关键问题之一。若模型的输入变量过多，会导致输入信息冗余，样本观测误差增加，模型复杂度随之提升，增大预报误差和计算时间。输入变量过少，则模型输入侧提供的预报信息不足，无法很好地解释输出变量的变化机理，难以得到准确的预报结果。预报因子的筛选主要利用因子与预报对象间的相关关系，剔除含有不相关和重复信息的因子，确保筛选出高信息量、强相关性的预报因子。现阶段，传统的预报因子选取方法有先验判断法、逐步回归法、相关系数法等。先验判断法容易受判断者主观意识的影响，缺乏客观性；逐步回归法与相关系数法只能处理线性问题，在处理非线性问题时存在较大偏差。朱双等[9]采用灰色关联分析来量化预报因子与预报对象的关联程度，并按关联度大小从众多的相关因子中挑选出对径流过程影响显著的预报因子；赵铜铁钢[10]，刘蕊鑫[11]等将互信息运用于神经网络径流预报模型输入变量的选取，并对结果进行了分析，证明该方法具有一定的实用性；纪昌明等[12]建立最大联合互信息模型进行预报因子筛选，新方法能够为预报模型提供更加科学的输入，提高模型的预报精度；闪丽洁等[13]分别以相关系数法、逐步回归法以及相关系数法-逐步回归法筛选出的因子作为预报的输入项，结果表明综合方法筛选出的预报因子组合可以取得较好的模拟效果；周育琳等[14]验证了相关系数法-主成分分析法结合的综合方法优选预报因子的效果优于单一方法。

在现有研究的基础上，本文引入一种新的衡量相关关系的度量指标——最大信息系数，并将其与主成分分析法结合，提出最大信息系数-主成分分析耦合方法(MIC-PCA)，应用于径流预报因子筛选。并选取BP人工神经网络作为预报模型，该模型广泛应用于复杂的非线性系统建模。将耦合算法筛选出的因子集输入到预报模型中以验证因子筛选的效果。研究结果表明，相比于现行方法，MIC-PCA法能够为水文预报模型提供更加准确科学的输入，提高模型的预报精度。

1 MIC-PCA耦合算法

1.1 最大信息系数

最大信息系数(The maximal information coefficient MIC)是一种基于互信息的度量二维变量间相关关系的指标，由麻省理工学院的David N.Reshef[15]等人于2011年提出。该方法是在互信息的基础上经不等间隔寻优与矫正处理[16]，相比于传统的相关关系度量指标(Pearson相关系数、Spearman相关关系、互信息等)，主要具有以下3个优点。

(1)普适性。MIC法总能找到一种网格划分，搜索最优的分割点计算其互信息值，使其能有效反应变量之间的任意函数关系(包括线性或非线性关系)，在函数的叠加等非函数关系的度量上也有优异的表现。

(2)稳健性。MIC不易受到观测样本中异常值的影响，水文资料一般序列较长，往往存在异常值，因此，MIC适合分析水文资料的相关关系。

(3)公平性。对于2组信息量相同的变量，其对应的MIC值也相同，可以放在同一等级上公平对待。当2个随机变量满足函数或函数叠加关系时，MIC依概率收敛到1；当2个随机变量相互独立时，MIC依概率收敛到0。MIC值理论上在[0,1]，相关程度一目了然，使其在不同关系中具有可比性。

因而MIC在相关关系辨识、特征选择等方面的应用范围更广，将MIC引入到径流预报因子筛选是可行的。

1.2 主成分分析

主成分分析(Principal Component Analysis，PCA)由美国运筹学家Salty在1977年提出，该方法[17]利用降维思想，在保证数据信息损失最小的前提下，经线性变换和舍弃一小部分信息，把多变量转化为少数几个综合变量(即主成分)。PCA的基本原理是借助于一个正交变换，将其分量相关的原随机向量转化成其分量不相关的新随机向量，在代数上表现为将原随机向量的协方差矩阵变换成对角形阵，在几何上表现为将原坐标系变换成新的正交坐标系，使之指向样本点散布最开的若干个正交方向，然后对多维变量系统进行降维处理，使之能以一个较高的精度转换成低维变量系统。对于一个矩阵来说，将其对角化即产生特征根及特征向量的过程，也是将其在标准正交基上投影的过程，而特征值对应的即为该特征向量方向上的投影长度，因此该方向上携带的原有数据的信息最多。主成分空间内每一个主成分表示转换后有效的新特征。张辉[18]等选用主成分分析法进行经济学指标筛选。迟国泰[19]等通过主成分分析法删除了因子负载小的指标。说明了主成分分析法在因子筛选及降维中的实用性。

传统主成分分析法的计算步骤简述如下。

(1)形成样本矩阵，样本标准化处理。

(2)计算样本矩阵的协方差矩阵。

(3)对协方差矩阵进行特征值分解，计算主成分贡献率及累计贡献率。

(4)按照一定的累计贡献率选取最大的n个特征值对应的特征向量组成投影矩阵， 得到降维后的新样本矩阵。

1.3 MIC-PCA

筛选预报因子主要目的是剔除含有较多冗余信息的因子，冗余信息既包括有效信息含量低，对预报对象影响不显著的信息，又包括重复信息，无论因子的信息含量高或低都可能含有重复信息。所以预报因子的筛选过程分为2个步骤：一是要筛选出信息含量高，对预报结果和预报精度影响显著的因子。这一步剔除了冗余信息中的无关信息，信息含量低的因子与预报变量相关程度很低甚至无关，输入到预报模型中不仅会增加模型复杂度，还会影响预报结果的准确性，可能造成对训练样本的拟合度不够，即“欠学习”。二是要剔除高信息量因子中重复信息含量较高的因子。强相关关系的高信息量因子之间存在着较大的信息重叠部分，即因子之间存在较多的信息冗余。这样的因子输入到预报方法中会使预报模型过为复杂，影响预报结果的准确性，可能产生“过学习”，增加预报模型的复杂程度。

现有的预报因子筛选方法，一部分以相关系数法、互信息法等为代表，其只考虑最大化因子集的信息含量，而无法保证筛选出的因子集中重复信息少；另一部分以主成分分析法为代表，只考虑对变量空间进行降维，获得代表原始变量的主要成分，而无法保证筛选出的因子集信息量最大化。目前，很少有方法可以达到2者兼顾，而且单一的方法也很难同时满足这2点要求。因此，考虑对这2种方法进行有机结合，为预报因子的选取工作提供了新的思路。

基于此，本文提出最大信息系数-主成分分析耦合方法(MIC-PCA)，应用于2阶段径流预报因子筛选中。第1阶段，通过最大信息系数有效衡量预报因子与实际径流系列的相关关系，筛选出与实际径流系列相关性强的因子，这些因子通常为高信息量因子，对于预报变量影响显著；第2阶段，从高信息量因子中剔除掉含有较多重叠信息的因子。

主成分分析的目的就是 “去冗余”。“去冗余”的目的就是使保留下来的维度间的相关性尽可能小，传统方法中以协方差来描述相关性，本文提出用MIC特征矩阵替代协方差矩阵来衡量“信息”的多少。MIC矩阵的主对角线上的元素与协方差矩阵主对角线上的元素一样具有“方差”的含义，其他元素是两两维度间的相关性度量。因此，用MIC特征矩阵替代协方差矩阵是合理的。通过计算预报因子之间的最大信息系数以建立预报因子MIC特征矩阵，采用主成分分析法剔除信息重叠的因子，得到最终的预报因子集，可保证因子集含有最大信息量的同时重叠信息较少。MIC特征矩阵替代协方差矩阵具有如下优点。

(1)协方差矩阵本身含有单位，需要变量进行归一化处理以消除单位带来的影响。

(2)协方差只能表征变量间的线性关系，而MIC可以衡量变量间的非线性关系甚至非函数关系[20]。

(3)MIC基于信息熵，本身就含有“信息”的属性，可以比协方差更好地度量信息量。

设初始预报因子集为X={X1,X2,…,Xn}，预报对象为Y，MIC[21,22]-PCA的具体计算步骤如下。

(1)在2个因子的散点图(集合D)上进行x×y划分，将其元素按x值划分到x个格子中，按y值划分到y个格子中。

(2)计算集合D的点落在给定的网格G上所得到的频率分布D|G，(允许某几个网格内没有落入数据集D中的点)。为了保证既不会因为网格过为细密造成每个样本点都有自己的小格，而导致即便对随机数据也有MIC≠0，又不会因为网格过为稀疏而不能精确地反映数据集的信息，需要合理选择网格G的划分上界B(n)，ω(1)

(3)计算不同的网格G确定不同的概率分布，即I(D|G) ，表示点集基于分布D|G的互信息。

(4)在基于x×y网格G的所有可能分布D|G的互信息中找到最大的互信息maxI(D|G)记为I*(D,x,y)。然后将不同大小网格G上的最大互信息标准化，使得各最大互信息均在(0,1)内，保证了公平性，得到二维数据集D的特征矩阵M(D)。

(1)

(5)从步骤(4)得到的结果中找出最大者，即为最大信息系数。

MIC(D)=maxx y

(2)

(6)分别计算各预报因子与实测径流间的最大互信息系数MIC(Xi,Y)，并按照从大到小的顺序依次排列，筛选出排名靠前的若干因子组成新的因子集X′={X1,X2,…,Xm}，m≤n。

(7)分别计算因子集X′中各因子之间的MIC值MIC(Xi,Xj)(i,j≤n)，组成MIC特征矩阵。

(8)计算MIC特征矩阵的特征值并使其按从大到小顺序排列λ1≥λ2≥…λm≥0。分别求出对应于特征值λi的特征向量ei(i=1,2,…,m)，要求‖ei‖=1。

(9)计算主成分贡献率及累计贡献率。一般取累计贡献率达85%～95%的p个特征值所对应的p个主成分，即得到最终的预报因子集X″={X1,X2,…,Xp}。

贡献率：

(3)

累计贡献率：

(4)

运用MIC-PCA进行预报因子筛选的流程图如图1所示。

2 实例分析

2.1 初始预报因子集选取

本文选取雅砻江流域打罗水文站日径流预报作为研究对象。打罗水文站位于官地水电站下游，官地水电站是雅砻江水电基地下游卡拉至江口河段规划的5个梯级电站之一。由于官地水电站上游未设有入库水文站，官地的入库径流一般通过“反推”得到。径流变化对于水电站水库调度、发电及农业灌溉等有重要影响，准确预报打罗水文站的径流量可以提前推算出官地的入流情况，对官地水电站提前制定调度方案，指导电站科学合理运行具有重要意义。打罗站位于雅砻江干流上，其径流主要由支流九龙河(乌拉溪站控制)汇入雅砻江干流(三滩站控制)形成，上游有泸宁站汇入。径流的传递具有时延性，所以上游站点已发生的当日流量直接影响着下游站点未来的流量。另外，由于流量序列具有自相关性，打罗站前期流量对当日流量的预报有着重要作用。

为验证本文提出的因子筛选方法的可行性与优越性，选取打罗站2009年汛期5-10月日流量作为预报对象，即Y={Qdl,t}，选取1998-2008年共10 a汛期(5-10月)打罗水文站前1 d至前5 d天然日径流和上游水文站三滩站、乌拉溪站、泸宁站的当日径流及前1 d至前5 d流量实测数据作为候选预报因子集：X={Qdl,t-1,…,Qdl,t-5,Qln,t,…,Qln,t-5,Qst,t,…,Qst,t-5,Qwlx,t,…,Qwlx,t-5}，采用MIC-PCA耦合法进行预报因子筛选。

图1 MIC-PCA法预报因子筛选流程Fig.1 Flow chart of prediction factor screening by MIC-PCA method

其中Qdl,t-1,…,Qdl,t-5代表打罗站前1 d至前5 d天然日径流，其他符号含义以此类推。所选各站径流资料由雅砻江公司提供，该公司之前也利用此数据资料作过相关的分析，可以充分反映官地上游来水情势以及水文变化情况，选取的10 a资料，时间跨度大，能够体现该流域径流在时间上的变化，该数据具有代表性好，可靠性高等特点。文中提及的水库、水文站点的布设示意图如图2所示。

图2 水文站点布设示意图Fig.2 Layout of hydrological stations

2.2 基于MIC-PCA法的预报因子筛选

第1阶段，计算因子集中各因子与打罗站当日径流间的最大信息系数MIC(Xi,Y)，计算结果如图3所示。分析图3中的计算结果可以得到以下结论。

(1)打罗站当日径流序列与其自身的MIC值为1，是完全相关关系，以该点作为参照点，那么打罗站当日径流序列与打罗站t-1～t-5、泸宁站和三滩站t～t-5的径流序列的MIC值均在0.6以上，说明这些预报因子与预报对象相关关系强，而与乌拉溪站的最大信息系数均位于0.6以下，且明显远小于其他因子，说明乌拉溪站径流与打罗站径流的相关性较弱，所含信息量相比于其他站点较少，可以剔除。

(2)打罗站、泸宁站和三滩站的最大信息系数曲线差距很小，基本重叠，表明这3个站点在同一时刻的因子所含信息量接近，之间可能存在重复信息，需要考虑进一步剔除简化；再次，每一个站点的曲线都大致呈下降趋势，说明打罗站当日径流序列与各站径流序列的相关程度随着时间间隔增加而减弱，即预报因子与预报对象的时间差越小，其包含的预报信息越多，该因子在预报中所起的作用越大。

(3)泸宁站与三滩站的曲线在t时刻到达最高点，且与t-1时刻接近，可以推测泸宁站与三滩站的径流只需不到1 d即可流经打罗站，乌拉溪站在t-1时刻到达最高点，且与t、t-2时刻的最大信息系数值非常接近，所以乌拉溪站的径流要经过1～2 d才能抵达打罗站。这与各站点间的布设位置和距离的实际情况相一致。综上所述，第1阶段筛选后得到的预报因子集为X′={Qdl,t-1,…,Qdl,t-5,Qln,t,…,Qln,t-5,Qst,t,…,Qst,t-5}，因子集由23个因子初步简化为17个。

图3 MIC(Xi, Y)计算结果Fig.3 Calculation results of MIC(Xi,Y)

第2阶段，在MIC法筛选出的因子集的基础上，采用改进的主成分分析法剔除冗余的预报因子，精简预报集。计算因子集X′中各因子之间的最大信息系数MIC(Xi,Xj)，得到17×17的最大信息系数矩阵；计算MIC矩阵的特征值向量及各因子的主成分贡献率向量。计算结果见表1。

表1 MIC-PCA法第2阶段计算结果Tab.1 Second phase calculation results of MIC-PCA method

显然，前3个因子的累计贡献率达到99%以上，说明前3个因子含有的预报信息量占17个因子的99%以上，其余因子间的相关关系很强，可以互相取代，存在较多冗余信息。说明泸宁站、三滩站与打罗站的径流十分接近，可能3个站在地理位置上距离邻近。因此完全可以用这3个因子代替17个因子，则最终筛选出的预报因子集为X″={Qdl,t-1,…,Qdl,t-3}，因子集由17个因子简化为3个。

2.3 方法对比

为了验证MIC-PCA法的有效性，本文选取BP人工神经网络作为预报模型，将筛选得到的预报因子集的观测样本作为模型的输入样本，对预报模型进行训练和测试，通过径流系列的预测结果检验预报因子集的筛选效果。采用1998-2008年共2 024组汛期数据进行模型训练，神经元传递函数选取双曲正切函数，隐含层神经元数量设为200，训练次数10 000次。预测打罗站2009年5月1日至10月31日的日平均流量过程，并采用因子筛选方法中应用较多的互信息法(MI)，与单一最大信息系数法(MIC)、传统主成分分析法(PCA)和本文提出的MIC-PCA法的因子筛选结果和预报结果进行对比。预报结果的评定采用计算时间、平均绝对误差MAE、均方根误差RMSE、确定性系数DC和预报合格率QR这5项指标，对各模型的拟合结果进行检验，评定指标的计算公式见表2。

表2 预报准确度计算指标和计算公式Tab.2 The calculation index of forecast accuracy

表2式中：MAE代表预测值与实际值的偏差绝对值求平均；RMSE代表预测值与实际值之差的平方根的平均值；DC代表预报过程与实测过程之间的吻合程度，其值越接近1，偏离程度越小；QR代表预报合格点所占比重，短期径流预报误差小于15%为合格预报；N为序列长度；Q(i)为实际流量值；Qf(i)为预测流量值；n为预报合格次数。

互信息和MIC按照80%的单侧置信度选取因子，主成分分析选取99%累积贡献率的因子，各种方法筛选出的因子集及预报效果见表3。

表3 各方法筛选结果与预报准确度指标值Tab.3 The screening results of various methods and the index value of prediction accuracy

结果显示：互信息和最大信息系数都是将与预报对象相关关系强的因子筛选出来，保证筛选出的因子集信息量最大。最大信息系数作为MIC-PCA的第1阶段，目的是筛选出含有高信息量的因子，与互信息相比，因其衡量非线性、非函数关系的能力更强，更容易发现互信息无法检测到的因子间的关联信息，故筛选出的因子更多。但所筛选出的因子集作为BP人工神经网络的输入表现不佳，计算时间最长，预报精度和合格率都最低，说明这些因子间信息重叠部分较大，影响了预报的准确度。

互信息与主成分分析的计算时间相差不大，但主成分分析的预报效果明显较好，主成分分析法基于各因子的线性组合，将系数大的识别为主成分，并通过降维剔除了大多冗余信息。但主成分分析法无法识别与预报对象非线性关系较强的预报因子，故遗漏了该部分因子所包含的信息，拟合效果欠佳。互信息法筛选出的因子集中保留了过多的因子，可能存在信息冗余，降低了BP人工神经网络的学习效果，拟合效果不佳。

MIC-PCA在所有方法中表现最佳，预报因子集最为精简，计算时间最短。虽然因子数量少，但衡量预报准误差的MAE和RMSE均最小，DC最接近1，说明预报径流与实际径流的偏差小，拟合精度高，且预报合格的点数多，预报效果最佳。采用MIC-PCA筛选出的因子既保证了信息量充足，又剔除了较多重复信息，使因子集合理精简保证了高质量、精准的模型输入，因此预报模拟效果好。这说明MIC-PCA可以为官地水电站的入库站径流预报提供准确的预报因子集。

采用的互信息法、主成分分析法和MIC-PCA法用于径流预报，3种方法的预报结果及实际径流过程线如图4所示。

图4 以各方法的因子集进行径流预报的结果Fig.4 The result of runoff forecast based on the factor set of each method

从图4可以看出，无论是在低流量还是高流量的情况下，MIC-PCA的预报径流曲线与实际径流曲线十分接近，预报准确度较高；互信息法在低流量附近和峰值附近拟合偏差较大，前期5月1日至6月10日预报值偏大，之后普遍偏低。根据互信息法筛选出的因子集进行预报偏差较大，平均每个预测点与真实值的相差350 m3/s，不利于水库汛期的安全运行；MIC-PCA法与主成分分析法的预测径流比较贴近，但可以看出MIC-PCA法的的预报径流曲线与实际径流曲线更加吻合，尤其是在极端值附近。实测径流序列最大峰值出现于2009年8月13日，为8 090 m3/s，MIC-PCA法预测为7 913 m3/s，而主成分分析法预测为7 273 m3/s，较实际值明显偏小，若使用该预报值作为调度决策的参考，使官地水电站推算出的入库径流量偏小，会发生预留的防洪库容不足的后果，对下游及水库安全造成损害。而MIC-PCA法的预报值与实际值始终较为接近，预报合格率高达93%，可以为官地水电站的调度决策提供更加可靠的参考。

3 结 论

为解决现有径流预报因子筛选方法较少能够兼顾信息量高、冗余信息少等要求，本文将最大信息系数衡量变量间复杂相关关系的能力，与主成分分析法降维去除冗余信息的能力有机结合，提出MIC-PCA耦合算法，并应用于径流预报因子筛选。第1阶段，计算各预报因子与实测径流间的最大互信息系数从而筛选信息含量高的因子；第2阶段，采用改进的主成分分析法，即以预报因子之间的最大信息系数矩阵代替传统方法中的协方差矩阵，剔除信息重叠的因子，最终得到筛选后的预报因子集。本文将该方法应用于打罗水文站的日径流预报中，并将该方法筛选出的预报因子集与其他多种因子筛选方法的筛选结果分别作为BP人工神经网络的输入以验证该方法的有效性，综合各项结果可以看出，该方法能够从众多预报因子中筛选出高信息量低冗余度的因子，为预报模型提供精确合理的输入，有助于提高预报模型的预报效果。综上，本文提出的MIC-PCA耦合算法在预报因子筛选方面具有实用价值。

□