数学建模-DHUMCM11016-A.doc

2011年东华大学数学建模竞赛论文赛题编号（A）（河流径流量预测与相互关系分析）参赛队号：DHUMCM11016参赛队员：王鹤翔（信息科学与技术学院，班级 自动化0902，学号090900417）罗世操（信息科学与技术学院，班级 自动化0902，学号090900424）俞承昊（信息科学与技术学院，班级 电气0902，学号090900425）2011年6月1日摘要:分别运用人工神经网络BP和RBF方法对各水文站的日径流量进行了拟合分析和预测,并对两者的输入变量的选择和个数的确定以及隐藏层、输出层单元数的确定等关键技术问题进行了探讨。由预测结果可以看到在输入变量中加入前几天的日径流量可以较大的提高预测的精度。最后对两者的预测效果进行了比较，无论拟合的过程还是预测检验的结果都与实测值基本符合，但同时也看到BP神经网络在汛期时由于水量变化过大，预测效果不及RBF神经网络。说明RBF预测模型更为合理。对于各水文站径流之间关系的确定，我们采用多元线性回归的逐步回归的方法，逐步剔除相关系数较小的水文站，最总确定它们之间的制约关系。对于大型水库的位置的选取，先从各流域的年径流量出发，找到年径流量在10亿立方米以上的流域作为参考点。然后观察到各水文站径流随季节的变化一致，剔除时间因素对选址影响，主要考虑地域的平衡。我们建立了一个理想化的模型来估计水库对水资源地域的重分配。理想地认为相邻的各水文距离为“1”，由于水是自高向低流，所以也理想地认为上游的海拔成一个固定倾角线性上升，这样水库拦截下的水资源可按等差数分配给每个水文站。最后根据水资源在空间上的匀称程度决定水库的位置。关键字: 径流量预测；神经网络；制约关系；水库选址问题的重述2011年，中国山东、山西、河北、湖北等地遭受百年一遇的特大旱灾，3月旱灾蔓延至江西、湖南等地，截至3月，全国耕地受旱面积1.16亿亩农田受旱，272万人饮水困难。旱情的预测和水资源的调节与利用已成为我们共同关注的话题。分析河流动态变化，快速而准确的进行水量预报，更好地利用河流水量，是我们刻不容缓的任务。根据河流动态的变化，快速准确的进行水量水情预报，为可控调节用水，旱情缓解提供一定的技术支持。图1是某流域一河流水文站网络概化图。为了能够更好地分析其动态水流变化，在河流各段分别设置了径流检测点，又称水文站，用于检测实时水量的变化。在此河流流域共设置距离不等的18个水文站，在附件的Excel文件中，记录了各个控制断面在一年内各天的径流变化的均值。这里假设径流的变化只与平均水位、降雨量之间密切相关，同时，各个水文站之间径流量也会相互作用。1、试用各水文站的平均水位、降雨量数据预测未来3天内各水文站的径流量（产汇流与径流存在区别，径流量的大小还会受到水库调蓄等因素的影响）。2、请分析各个控制断面间径流量相互制约关系，表1给定了流域多年平均水量消耗，实际可看作用水比例。3、如果让你在此流域设置二个大型水库，用于平衡调节各地区径流量及抗旱蓄水，请问如何设置此水库位置，并分析原因。一、水文站的未来3天径流量预测1 问题的分析各个水文台的径流量因素涉及到各个综合复杂的因素，一个水文台在预测径流量时必须得考虑当日的降水量和当日的平均水位值。如何给出一个水文台的径流量预测的参考模型尤为重要。由于径流量的值受到各种环境复杂因素的影响，应在主要的环境因素（降水量和平均水位）建立模型。模型的建立就是根据现有数据从降水量和平均水位出发在基本约束条件下误差最小的径流量标准。再以此为依据对现在各水文台的径流量进行量化分析，提出合理的参照标准。2 模型假设针对本问题，提出以下合理假设：(1)文中所搜数据能较真实反映各个水文台的事实情况。(2)水文台的周围温度，用水量以及预期短期内不考虑降水量暴增应保持稳定。(3)对所搜数据能较好适用于对现阶段的水文台径流量分析,径流量的短期微小变化对总体预测没有太大影响。3 模型的建立与求解3.1 预测径流量模型的建立基于BP神经网络的水文台径流量预测模型径流量的模型BP 神经网络的数学模型的建立，首先是基于下面最基本的数学公式。沿信息传播的方向，给出网络的状态方程，用和 表示第层第 个神经元的输入和输出，网络的各层输入输出关系可描述为：第一层（输入层）：将输入引入神经网络。 第二层（隐层）： 其中是激发函数，可以取不同形式，如：S 函数：高斯基函数：以及径向基函数、样条基函数、小波函数等。第三层（输出层）： 基于本文模型的条件：i=2,l=5,k=1，以上就完成了神经网络的基本网络构建。下一步是要确定网络的学习算法。学习的基本思想是：通过一定的算法调整网络的权值，使网络的实际输出尽可能接近期望的输出。在本网络中，采用误差反传（BP）算法来调整权值（即 BP 神经网络）。其基本原理是：当网络的输入（即相应影响因素）为 时，网络的实际输出为，网络的期望输出（即实际的径流量值）为 。定义学习的目标函数为：（均方差法）同时 BP 算法通过下列公式来调整权值，使目标函数最小：（最速下降法）其中为学习率。其中具体分析每一层神经网络，可以得到：其中基于以下公式，偏差逐步反传：得出各因素的处理后权值： 3.2 模型的求解1.准备和数据预处理神经网络是一种适宜处理具有残缺结构和能够分析含有错误成分的算法,它能够在信息含糊、不确定、不完整、存在矛盾及假象等复杂环境中处理分析数据；并且神经网络所具有的自学能力，使得传统数学算法应用最为困难的有效数据获取工作，转换为网络的变结构调节过程,从而大大方便了各种不同应用对象的建模与分析，进而可以对一些复杂问题做出合理的判断决策,做出有效的预测和估计。 训练数据如下： 表 1 HS2站1,2,3月前20天逐日降水量表和平均水位表1月2月3月30.0830.2630.74082175071956082031.17072875072630072740073189073632.93114352083731.58113442163317173030.93132786122387112031.43101929092034103330.951月2月3月 0.6 5.9 3.2 6.8 24.0 2.6 4.6 24.1 4.6 1.2 16.0 6.1 0.6 25.0 1.4 2.7 18.1 6.1 11.7 2.3 3.3 2.6 0.2 3.7 6.0 2.1 20.6 2.8 8.6 3.8 0.7 17.8 0.5 5.2 2.2 采用 MATLAB神经网络工具箱进行回水河径流模拟包括以下步骤:( 1) 从数据文件中导入训练样本和检验样本。P（一月份前十天的降水量和平均水位）为训练样本输入矢量；T（一月份前十天的径流量实测值）为训练样本目标矢量，PP（后十天）为检验样本输入矢量，TT为检验样本目标矢量。(2)神经网络初始化net=newff(30.07 30.11;0 4.6,5,1,tansig,purelin,traingd);。以上分别设置了 BP 神经网络的隐含层神经元数、 输出层神经元数、 神经元采用的激活函数以及网络训练算法。利用函数 initff( ) 对网络各层权重和阈值进行初始化。此外还需设置一些网络训练参数:net. trainParam. epochs = 50000; % 设 置网 络 训练 次数 为50000次。net. trainParam. show = 50; % 设置每隔 50 次显示当前网络训练误差。net. trainParam. goal= 0.0001; %设置网络训练目标误差为 0. 1。（3）对BP神经网络进行训练和仿真。net,tr=train(net,P,T); %利用函数 train( ) 对网络进行训练。A= sim( net, P) ; AA= sim( net, PP) ; %利用函数 sim( ) 对训练年份和检验年份日流量进行仿真。( 4) 将上述代码保存为 M 文件, 运行该文件即可进行日流量的模拟。模拟效果：表 2 是对HS2站1月份的前20天的径流量的实测值和BP网络模型下的预测值表2 1月 径流量实测值一月径流量预测值 0.60.61320.60.60210.530.52990.60.57280.530.52660.530.53290.530.53530.530.53250.530.53530.870.86950.60.71270.870.86901.521.40461.71.65041.11.38970.990.74700.870.87180.750.72910.680.67010.750.7291然而BP神经网络在数值偏差过大的时候拟合效果欠佳，所以以下引入RBF网络对数据进行拟合（以下是HS1站的数据）：3.3 RBF神经网络以HS1站为例，HS1站数据取于3月1日到4月3日，预测3月31日后3天的径流量。其中，RBf的输入变量为前3天的径流量，当日降水量、水位、径流量，输出目标为后3天的径流量。训练量本为前27组数据，测试样本为31日的一组数据。运用MATLAB进行数据拟合，预测：拟合预测效果见下图clear; p=p1;T=T1;p_train=p(:,1:27);T_train=T(:,1:27);p_test=p(:,28); T_test=T(:,28);PN1,minp,maxp,TN1,mint,maxt = premnmx(p_train,T_train);%归一化处理PN2 = tramnmx(p_test,minp,maxp);TN1,mint,maxt = premnmx(T_train);TN2 = tramnmx(T_test,mint,maxt); goal=0.01; spread=1.5; mn=28;df=1;net=newrb(PN1,TN1,goal,spread,mn,df) y=sim(net,PN2)%预测 y = -0.7567 -0.6850 -0.5582y1=postmnmx(y,mint,maxt)y1 = 518.8337 571.5238 664.7098由右下图可知，对于前两天的预测效果较好，对于第三天有一定偏差，但不大。基本符合预期。 HS1 站径流仿真与预测 注：红线为实际径流值 绿色加号为仿真结果 蓝圈为未来三天径流值同理，在相同时间内取值可以得到其他水文站的后三天径流值。HS2 站径流仿真与预测 HS3站径流仿真与预测HS4 站径流仿真与预测 HS5 站径流仿真与预测HS6 站径流仿真与预测 HS7站径流仿真与预测 HS8站径流仿真与预测 HS9站径流仿真与预测 HS10站径流仿真与预测 HS11站径流仿真与预测 HS12站径流仿真与预测 HS13站径流仿真与预测HS14站径流仿真与预测 HS15站径流仿真与预测HS16站径流仿真与预测 HS17站径流仿真与预测 HS18站径流仿真预测4 径流量的求解由上述结果可得降水量和平均水位对径流量求解如下：选取HS11站的降水量和平均水位为样本进行讨论预测6月5日到6月7日的径流量：原始数据选取为HS11站在6月份的样本数据；最终由MatlabM编程运行可得径流量为：T5= 173.4210，T6= 173.4210，T7=173.4230结果的分析：HS11站在6月5日到6月7日的径流量实际值为：TT5=174;TT6=170;TT7=175 二:各水文站之间的相互制约1. 问题分析假设各水文站之间只有直接相关的水文站才会相互影响，不直接相关的水文站影响微小；因此将18个水文站分为，1.HS1，HS2，HS3，HS8;2. HS2,HS1,HS8;3. HS3,HS1,HS5，HS4，HS6;4. HS4,HS3,HS5;5. HS5,HS3,HS4; 6. HS6，HS7，HS3; 7. HS7,HS6;8. HS8，HS12，HS1，HS9，HS11;9. HS9,HS11,HS12,HS8;10. HS10,HS11;11. HS11,HS10，HS9，HS8,HS12; 12.HS12,HS13,HS14,HS15,HS8,HS11,HS9;13.HS13,HS12,HS14,HS15;14.HS14,HS15,HS13,HS12;15.HS15,HS16,HS12,HS14,HS13;16.HS16,HS15,HS17;17.HS17,HS15,HS16,HS18;18.HS18,HS172.确立模型以上18个组第一个水文台作为因变量Y，后几个水文台作为自变量X，以下径流量均为每年各月的平均径流量。运用MATLAB工具箱进行回归分析。以下采用逐步回归的方法1. x1=2.1 11.0 36.6 50.2 58.5 31.0 7.8 37.3 16.7 8.5 3.3 2.7;x2=32.9 100.0 231.0 451.0 584.0 242.0 144.0 294.0 172.0 88.2 56.5 44.7;x3=71.3 145.0 414.0 687.0 990.0 423.0 204.0 354.0 242.0 128.0 114.0 83.1;y=80.3 203.0 702.0 1120.0 1830.0 642.0 266.0 678.0 439.0 147.0 130.0 99.4;x=x1 x2 x3;stepwise(x,y)因为x3显著性最差，所以移去x3Variables have been created in the current workspace. b=regress(y,x)b = 0.5345 0.0834 1.64383.分析结果所以HS1与HS2，HS3，HS8之间的关系为Y= 0.5345+0.0834*x1+1.6438*x2,即HS1径流量=HS2径流量*0.0834+HS3径流量*1.6438+0.5345同1可得个水文站的径流量为：1.x1=80.3 203.0 702.0 1120.0 1830.0 642.0 266.0 678.0 439.0 147.0 130.0 99.4;x2=71.3 145.0 414.0 687.0 990.0 423.0 204.0 354.0 242.0 128.0 114.0 83.1;y=2.1 11.0 36.6 50.2 58.5 31.0 7.8 37.3 16.7 8.5 3.3 2.7;x=x1 x2;stepwise(x,y)b=regress(y,x)HS2（径流量以下省略）=0.0057+0.0585*HS12.y=32.9 100.0 231.0 451.0 584.0 242.0 144.0 294.0 172.0 88.2 56.5 44.7;x=x1 x2 x3 x4;（xi为各组作为自变量HSj的径流量，下同）stepwise(x,y)b=regress(y,x)HS3=0.0337+0.7103*HS1-1.6908*HS4+2.8705*HS53.y=2.9 13.3 30.7 72.9 77.7 19.6 8.6 41.2 14.8 8.2 4.7 3.4;x=x1 x2;stepwise(x,y)b=regress(y,x)HS4=1.4261*HS5-0.00354.y=2.9 8.8 16.6 50.1 54.6 16.1 11.2 31.8 16.8 7.2 4.4 4.1;x=x1 x2;stepwise(x,y)b=regress(y,x)HS5=0.4087*HS4+0.04135.y=11.9 31.7 66.4 156.0 200.0 80.3 48.2 101.0 56.6 32.4 23.3 19.9;x=x1 x2;rstool(x,y,purequadratic)beta,rmseHS6=5.6643+0.3474*HS3-3.7393*HS7-0.0001*HS3+0.5527*HS76.x=11.9 31.7 66.4 156.0 200.0 80.3 48.2 101.0 56.6 32.4 23.3 19.9;X=ones(12,1) x;Y=0.8 2.2 3.8 10.1 11.5 3.5 2.1 6.6 3.1 1.3 1.0 0.9; b,bint,r,rint,stata=regress(Y,X);b,bint,statsHS7=0.0617*HS6-0.34587.y=71.3 145.0 414.0 687.0 990.0 423.0 204.0 354.0 242.0 128.0 114.0 83.1;x=x1 x2 x3 x4;stepwise(x,y)b=regress(y,x)HS8=0.2936+101.0359*HS9+1.3466*HS11-0.2328*HS128.y=0.2 0.8 1.1 0.9 0.6 0.8 0.3 0.3 0.4 0.2 0.2 0.2;x=x1 x2 x3;rstool(x,y,purequadratic)beta,rmse HS9=0.1614-0.0004*HS8+0.0015*HS11+0.0001*HS129.x=35.2 36.0 124.0 192.0 315.0 140.0 79.4 125.0 69.0 30.6 45.4 26.5;X=ones(12,1) x;Y=5.2 16.6 31.4 61.2 76.5 33.2 22.7 31.0 17.9 13.5 9.8 8.3;b,bint,r,rint,stata=regress(Y,X);b,bintHS10=2.0859+0.2481*HS1110.y=35.2 36.0 124.0 192.0 315.0 140.0 79.4 125.0 69.0 30.6 45.4 26.5;x=x1 x2 x3 x4;rstool(x,y,purequadratic)beta,rmseHS11=9.8692+0.2520HS8+117.0822*HS9+2.1422*HS10-1.0614*HS12+0.0002*HS82-101.5969*HS92-0.0419*HS102+0.0042*HS12211.y=37.0 155.0 187.0 140.0 101.0 169.0 72.4 55.9 98.0 48.2 38.0 35.8;x=x1 x2 x3 x4 x5 x6;stepwise(x,y)Variables have been created in the current workspace. regress(y,x)HS12=-0.6781+212.2589*HS9+0.7556*HS11-2.2428*HS14+1.0797*HS13+0.1808*HS1512.y=5.2 16.6 31.4 61.2 76.5 33.2 22.7 31.0 17.9 13.5 9.8 8.3;x=x4 x5 x6;stepwise(x,y)Variables have been created in the current workspace. regress(y,x)HS13=0.02503+0.1505*HS12+0.0067*HS1513.y=26.7 74.1 187.0 373.0 490.0 176.0 111.0 168.0 124.0 73.0 48.4 43.4;x=x4 x5 x6;stepwise(x,y)Variables have been created in the current workspace. regress(y,x)HS14=-0.1517+6.5471*HS13-0.0358*HS1514.y=26.1 83.0 52.2 53.4 296.0 372.0 105.0 76.3 74.2 67.2 35.0 28.0;x=x4 x5 x6 x7;stepwise(x,y)Variables have been created in the current workspace. regress(y,x)HS15=0.3366*HS12+1.6417*HS13-1.0202*HS14+17.9060*HS1615.y=3.1 6.7 11.7 22.1 31.4 12.3 9.5 13.5 7.5 6.4 5.0 4.2;x=x6 x7;stepwise(x,y)Variables have been created in the current workspace. regress(y,x)HS16=0.015+0.1973*HS1716.y=9.3 20.4 57.6 112.0 137.0 39.1 30.6 57.6 37.7 24.5 15.9 12.0;x=x6 x7 x8;stepwise(x,y)Variables have been created in the current workspace. regress(y,x)HS17=-0.0496+3.0835*HS16+3.0132*HS1817. x=9.3 20.4 57.6 112.0 137.0 39.1 30.6 57.6 37.7 24.5 15.9 12.0; X=ones(12,1) x; Y=1.4 2.4 5.1 15.3 18.1 5.8 4.8 7.0 4.1 3.5 2.5 2.0; b,bint,r,rint,stata=regress(Y,X); b,bintHS18=0.013+