黄河小浪底调水调沙问题.doc

黄河小浪底调水调沙问题内容摘要:为了确定排沙量与时间、排沙量与水流量的函数关系，我们可以用SAS软件做线性回归得到排沙量与时间的函数关系式，再利用所求函数在区间0,24上进行积分得到总排沙量1.93962亿吨。对于排沙量与水流量之间的关系，按时间分为两段进行拟合，最终用MATHLAB软件来画出图像，确定排沙量与排水量之间的函数关系式。关键词：调水调沙实验，sas,排沙量，排水量，matlab，拟合，线性回归问题的提出：在小浪底水库蓄水后，黄河水利委员会进行了多次试验，特别是2004年6月到7月进行的黄河第三次调水调沙试验具有典型的意义。这次试验首次由小浪底、三门峡和万家寨三大水库联合调度，进行接力式防洪预泄放水，形成人造洪峰进行调水调沙试验成功。这次试验的一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底的库区的沉积泥沙。在小浪底开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2720m3/s，使小浪底水库的排沙量也不断的增加。表一是由小浪底观测站从6月29日到7月10日检测到的试验数据。问题分析：1、对于问题一，所给数据中水流量x和含沙量h的乘积即为该时刻的排沙量y即：y=hx。2、对于问题二，研究排沙量与排水量的关系，从实验数据中可以看出，开始排沙量随水量增加而增加，而后随水流量的增加而减少，显然变化关系并非线性的关系，为此，把问题分为两部分，从水流量增加到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间的函数关系。模型假设：1、水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。2、时间是连续变化的，所取时间点依次为1,2,3，,24,单位时间为12h.模型的建立与求解：对于问题一，因为排沙量与时间的散点图基本符合正态曲线，如图二所示。所以，排沙量的对数与时间的函数关系就应该符合二次函数关系，因而排沙量取对数后，再与时间t进行二次回归，排沙量取自然后的数据见表2.假设排沙量与时间函数关系的数学模型是两边取对数得Lny=at2+bt+c先由表二做出排沙量的自然对数lny与时间t的散点图见图一，并利用SAS软件进行拟合，得到排沙量的自然对数与时间的回归方程为：Lny=-0.0209t2+0.4298t+10.6321由回归拟合参数表可知回归方程是显著的，因为相关系数人R2=0.9629,误差均方S2=0.0543，说明回归曲线拟合效果很好。所以排沙量与时间之间的函数关系式为图二：排沙量对时间的曲线图最后对所求出的函数关系在区间0,24之间进行积分结果为总排沙量1.93962亿吨，此与媒体报道的排沙量几乎一样。对于第二个问题，两个阶段的数据如表三表四所示表三：第一阶段试验数据序号1234567891011水流量x18001900210022002300240025002600265027002720含沙量h326075859098100102108112115表四：第二阶段的试验观测数据序号12345678910111213水流量x265026002500230022002000185018201800175015001000900含沙量h11611812011810580605040322085对于第一阶段，有表四用MATLAB作图（如图三）可以看出其变化趋势，我们用多项式做最小二乘拟合。设三次拟合函数关系h=a0+a1x+a2x2+a3x3其中a0,a1,a2,a3,为待定系数。四次拟合函数关系h= a0+a1x+a2x2+a3x3+a4x4其中a0,a1,a2,a3,a4为待定系数。图三：第一阶段水流量与排沙量之间的关系图三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数h=0.0032x2-2.4929x3,拟合效果如图四所示图四：三次多项式拟合效果，红线为拟合曲线类似的四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x3-7.4347x4,拟合效果如图五所示图五：四次多项式拟合效果，蓝线线为拟合曲线对于第二阶段，有表五用MATLAB作图可以看出其变化趋势，我们用多项式做最小二乘拟合。设三次拟合函数关系h=a0+a1x+a2x2+a3x3其中a0,a1,a2,a3,为待定系数。四次拟合函数关系h= a0+a1x+a2x2+a3x3+a4x4其中a0,a1,a2,a3,a4为待定系数。三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601.则拟合函数h=-0.9475x2+464.9601x3,拟合效果如图图六所示图六：三次拟合函数拟合效果类似的四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数h=-0.0013x2+1.1219x3-354.5952x4,拟合效果如图七所示图七：四次拟合函数拟合效果结论以及分析检验：用SAS软件做线性回归得到排沙量与时间的函数关系式为再利用所求函数在区间0,24上进行积分得到总排沙量1.93962亿吨，这与现实情况基本相符。对于第一阶段三次多项式拟合由MATLAB拟合函数求解a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数h=0.0032x2-2.4929x3对于第一阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x3-7.4347x4对于第二阶段三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601.则拟合函数h=-0.9475x2+464.9601x3,对于第二阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数h=-0.0013x2+1.1219x3-354.5952x4讨论与推广：1、对于第一个问题排沙量与时间不是严格的正态函数关系可能与实际有些偏差，此外还可以用SAS软件进行高次的多向式回归2、对于第二个问题，由于MATLAB软件的计算可能有些偏差导致拟合的函数关系可能与实际有稍微偏差，此外，还可以进行高次的拟合。附录：1、排沙量与时间的关系图像的MATLAB程序：t=1:1:24;y=57600,114000,157500,187000,207000,235200,250000,265200,2862000,2400,312800,307400,306800,300000,271400,231000,160000,111000,91000,54000,45500,30000,8000,4500; plot(t,y,r)2、对排沙量求自然对数的MATLAB程序与结果：y3=log(y)y3 = Columns 1 through 17 10.9613 11.6440 11.9672 12.1389 12.2405 12.3682 12.4292 12.4882 12.5644 12.6195 12.6533 12.6359 12.6340 12.6115 12.5113 12.3502 11.9829 Columns 18 through 24 11.6173 11.4186 10.8967 10.7255 10.3090 8.9872 8.41183、第一阶段的排沙量与水流量之间的关系MATLAB程序：x=1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720; h=32,60,75,85,90,98,100,102,108,112,115; x1=2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1500,1000,900; h1=116,118,120,118,105,80,60,50,40,32,20,8,5; plot(x,h,r:)4、第一阶段三次多项式拟合函数以及拟合效果程序与结果： A1=polyfit(x,h,3)Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. In polyfit at 80A1 = 0.0000 -0.0000 0.0032 -2.4929 z1=polyval(A1,x);plot(x,h,k+,x,z1,r)5、第一阶段四次多项式拟合函数以及拟合效果程序与结果：A2=polyfit(x,h,4)Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. In polyfit at 80A2 = -0.0000 0.0000 -0.0000 0.0121 -7.4347 z2=polyval(A2,x); plot(x,h,*,x,z2,r)6、第二阶段三次多项式拟合函数以及拟合效果程序与结果：A3=polyfit(x1,h1,3)Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. In polyfit at 80A3 = -0.0000 0.0006 -0.9475 464.9601 z3=polyval(A3,x1); plot(x,h,*,x1,z3,b)7、第二阶段四次多项式拟合函数以及拟合效果程序与结果：A4=polyfit(x1,h1,4)Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. In polyfit at 8