周 娅 ，郭 萍 ，杨 柳 ，宋培培

（1.贵州省水利水电勘测设计研究院，贵州 贵阳 550002；2.中国农业大学水利与土木工程学院，北京 100083）

岗南水库位于河北省石家庄市平山县境内，位于滹沱河干流上，是一座以防洪为主，兼顾灌溉、发电、工业及城市生活用水等综合利用的大（I）型水库。水库承担着石家庄和下游铁路、公路、华北油田及冀中平原的防洪保护任务，并保障区域农业、工业和城市用水安全[1]。小觉水文站是岗南水库入库径流控制水文站，上游来水相对天然，且受水库调节作用影响较小，长系列水文数据代表性较好。本文首先以主成分分析法提取出汛期气象因子的主成分，然后通过A Trous小波提取出各序列的随机项、趋势项和周期项，分别对各项建立BP神经网络预测模型。

1 研究区概况

滹沱河属于海河流域的子牙河水系，全长587 km，流域面积24690 km2。河流流经山西省境内的繁峙等五县市，于阳泉市盂县闫家庄进入河北省，然后由石家庄市区北穿京广铁路，向东横贯河北省，并从衡水市的安平县一直向东流至献县臧桥，在此和滏阳河相汇流入大海[1]（见图1）。流域年气温为4.0℃～8.8℃，年内降水主要集中在6月～9月，约占全年降水量的70%～80%，上游位于太行山的背风坡，雨量较小，径流量也较低，年径流深仅为20 mm～50 mm左右，而经太行山的东坡进入河北省时，降雨量增加，径流量也随之骤增，年径流深在50 mm～150 mm之间，然后经黄壁庄水库进入山麓冲积平原，年径流再次减少，径流深降至50 mm以下。

图1 滹沱河流域范围

2 方法简介

2.1 主成分分析法

1933年，Hotelling首先提出主成分分析，该方法利用降维的思想，以信息损失量很少为前提，将多个指标转换为几个综合指标，由原始变量的线性组合得到每一个综合指标，每个综合指标之间互不相关[2]。

主成分分析的数学模型为

可用矩阵表示为

经过标准化处理的矩阵就是X的相关系数矩阵R，如果主成分分析的一切计算不是从协方差矩阵出发，而是直接从样本相关系数矩阵出发，就等于先对数据进行标准化，然后再利用协方差矩阵进行主成分分析[2]。主成分分析的基本步骤是：首先结合数据，判断是否需要进行主成分分析；其次结合主成分的特征值和累计特征值贡献率来确定抽取的主成分或因子的数据；最后，进行主成分分析，将抽取出来的主成分分析存为新变量，以便继续分析[3]。

2.2 小波变换—A Trous

对水文时间序列f（t），t=1，2，…，N进行小波分解，令C0（t）=f（t），A Trous小波分解过程如下：

A Trous小波重构过程如下：

根据样本序列长度，低通滤波器选用db4小波，选择最大分解尺度P=3，采用MATLAB小波分解工具箱分解得到样本的随机项、趋势项和周期项，采用BP神经网络对各项分别进行预测，再重构得到预测模型的输出值[4]。

3 汛期月径流预测

汛期采用1969年～2009年的7月、8月、9月、10月的16个气象因子：x1五台山降水量、x2五台山风速，x3五台山气压、x4五台山气温、x5五台山相对湿度、x6五台山日照时数、x7五台山水汽压、x8原平降水量、x9原平风速、x10原平气压、x11原平气温、x12原平相对湿度、x13原平日照时数、x14原平水汽压、x15小觉站降水量、x16小觉站蒸发量，以及小觉站径流量值作为分析基础。

因子分析前，首先进行KMO检验和巴特利球体检验。KMO检验用于检查变量间的相关性和偏相关性，取值在0～1之间。KMO统计量越接近于1，变量间的相关性越强，因子分析的效果越好。实际分析中，KMO统计量在0.7以上时效果比较好；当KMO统计量在0.5以下时不适合应用因子分析法。采用主成分分析方法分析以上数据，得出汛期实测数据对应的KMO值为0.819（大于0.8），意味着变量间的相关性较强，说明变量很适合作因子分析。巴特利特球形检验法是以相关系数矩阵为基础的。它的零假设相关系数矩阵是一个单位阵，即相关系数矩阵对角线的所有元素均为1，所有非对角线上的元素均为0。巴特利特球形检验法的统计量Sig值是根据相关系数矩阵的行列式得到的。如果Sig值较大，且其对应的相伴概率值小于指定的显著水平0.05时，拒绝零假设，表明相关系数矩阵不是单位阵，原有变量之间存在相关性，适合进行主成分分析；反之，数据不适合进行主成分分析。对本次采用的数据进行Bartlett的球形度检验，其Sig值小于显著水平0.05，说明变量之间存在相关关系。所选取的影响因子的数据符合因子分析的条件，适合做主成分分析。进行主成分分析时，由表1可知前4个主成分贡献率累计百分比为88.943%，对应的初始特征值分别是 8.364、2.975、2.190、0.702，根据初始特征值计算得到主成分表达式的系数矩阵（表2），进一步计算得到汛期四个主成分 Y1、Y2、Y3、Y4作为预测输入。

表1 解释的总方差（汛期）

表2 主成分表达式的系数矩阵（汛期）

由于小波分析具有将时间序列分解为多层的能力，可分离季节性趋势时间序列中的趋势项、周期项和随机项，分离后各项所含频率成分单一，有利于建模预测。本研究将主成分分析提取得到的四个主成分Y1、Y2、Y3、Y4和小觉站月径流值Z进行A Trous小波变换分解，提取得到各数据序列的随机项（Y1D、Y2D、Y3D、Y4D、ZD），趋势项（Y1A、Y2A、Y3A、Y4A、ZA），周期项（Y1P、Y2P、Y3P、Y4P、ZP）。在此基础上，采用MATLAB自带BP神经网络训练工具箱（ntstool）建立各项预测模型，将以上数据作为输入信号从输入层进入网络，经隐含层逐层处理后到达输出层。每层的神经元状态只影响下一层的神经元状态。如果输出层得不到期望输出，则误差反向传播，根据误差调整网络权值和阈值，从而使BP神经网络预测输出不断逼近目标值，从而建立小觉站汛期月径流预测模型，模型采用的是训练较快的trainscg函数。

3.1 随机项

随机项预测模型网络输入层、隐含层和输出层神经元个数分别为5个、2个和1个。

图2 汛期随机项各数据集观测值和预测值回归分析图

图3 汛期随机项观测值和预测值对照图

图2 训练结果显示，训练数据集、测试数据集和检验数据集中预测值与观测值的相关系数R分别为0.9297、0.9553、0.9300，观测值与预测值的相关系数为0.9332，实测值和预测值吻合精度很高，同时模型对应数据集的均方误差（MSE）分别为 2.6×10-2、1.3×10-2、1.9×10-2，说明此模型也可以较好的模拟随机项数据集，从预测值和实测值对比图（图3）可以看出，模型对随机项的模拟精度很高，所选取的气象因子能很好的反映小觉站汛期径流量变化的随机性。

3.2 趋势项

趋势项预测模型网络结构输入层、隐含层和输出层神经元个数分别为5个、2个和1个。

图4 汛期趋势项各数据集观测值和预测值回归分析图

图5 汛期趋势项观测值和预测值对照图

图4训练结果显示，训练数据集、测试数据集和检验数据集中预测值与观测值的相关系数R分别为0.9970、0.9965、0.9986，观测值与预测值的相关系数为0.9972，实测值和预测值吻合精度很高，同时模型对应数据集的均方误差（MSE）分别为 1.4×10-3、1.2×10-3、1.0×10-3，说明此模型也可以较好的模拟趋势项数据集，从预测值和实测值对比图（图5）可以看出，模型对趋势项的模拟精度很高，可见选取的气象因子能很好的反映小觉站汛期径流量变化的趋势。

3.3 周期项

周期项预测模型网络结构输入层、隐含层和输出层神经元个数分别为5个、4个和1个。

图6 汛期周期项各数据集观测值和预测值回归分析图

图7 汛期周期项观测值和预测值对照图

图6 训练结果显示，训练数据集、测试数据集和检验数据集中预测值与观测值的相关系数R分别为0.8614、0.8631、0.9289，观测值与预测值的相关系数为0.8737，实测值与预测值吻合精度较高，同时模型对应数据集的均方误差（MSE）分别为 3.4×10-2、2.2×10-2、2.5×10-2，说明此模型也可以较好的模拟周期项数据集，从预测值和实测值对比图（图7）可以看出选取的气象因子能较好好的反映小觉站汛期径流量变化的周期性。

4 小结

图8 汛期观测值和预测值对照图

本文采用主成分分析法提取出汛期气象因子的主成分，然后采用A Trous小波提取出序列的随机项、趋势项和周期项，分别对各项建立神经网络预测模型，得到较好的预测结果，模型能很好的反映随机项、趋势项的变化，对周期项变化的响应较随机项和趋势项的差，但也具有较高的精度，说明所建立的模型能很好的反映径流的随机性和趋势性。采用A Trous小波重构法得到径流预测值，得到汛期重构的预测值和小觉水文站实测值的对比图（图8），从预测值与实测值对比来看，模型能较好的模拟小觉站月径流数据序列中的极值。