黄河小浪底调水调沙问题
黄河小浪底调水调沙问题
黄河小浪底调水调沙问题内容摘要:为了确定排沙量与时间、排沙量与水流量的函数关系，我们可以用SAS软件做线性回归得到排沙量与时间的函数关系式，再利用所求函数在区间[0,24]上进行积分得到总排沙量1.93962亿吨。对于排沙量与水流量之间的关系，按时间分为两段进行拟合，最终用MATHLAB软件来画出图像，确定排沙量与排水量之间的函数关系式。
关键词:调水调沙实验，sas,排沙量，排水量，matlab，拟合，线性回归问题的提出:在小浪底水库蓄水后，黄河水利委员会进行了多次试验，特别是2004年6月到7月进行的黄河第三次调水调沙试验具有典型的意义。这次试验首次由小浪底、三门峡和万家寨三大水库联合调度，进行接力式防洪预泄放水，形成人造洪峰进行调水调沙试验成功。这次试验的一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底的库区的沉积泥沙。在小浪底开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2720m^3/s，使小浪底水库的排沙量也不断的增加。表一是由小浪底观测站从6月29日到7月10日检测到的试验数据。
问题分析:1、对于问题一，所给数据中水流量x和含沙量h的乘积即为该时刻的排沙量y即:y=hx。
2、对于问题二，研究排沙量与排水量的关系，从实验数据中可以看出，开始排沙量随水量增加而增加，而后随水流量的增加而减少，显然变化关系并非线性的关系，为此，把问题分为两部分，从水流量增加到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间的函数关系。
模型假设:1、水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。
2、时间是连续变化的，所取时间点依次为1,2,3，…,24,单位时间为12h. 模型的建立与求解:<一>对于问题一，因为排沙量与时间的散点图基本符合正态曲线，如图二所示。所以，排沙量的对数与时间的函数关系就应该符合二次函数关系，因而排沙量取对数后，再与时间t进行二次回归，排沙量取自然后的数据见表2. 假设排沙量与时间函数关系的数学模型是
at^2，bt，c y,e
两边取对数得
Lny=at^2+bt+c
先由表二做出排沙量的自然对数lny与时间t的散点图见图一，并利用SAS软件进行拟合，得到排沙量的自然对数与时间的回归方程为:
Lny=-0.0209t^2+0.4298t+10.6321
由回归拟合参数表可知回归方程是显著的，因为相关系数人R^2=0.9629,误差均方S^2=0.0543，说明回归曲线拟合效果很好。
所以排沙量与时间之间的函数关系式为
,0.0209t^2，0.4289t，10.6312 y,e
图二:排沙量对时间的曲线图
最后对所求出的函数关系在区间[0,24]之间进行积分
结果为总排沙量1.93962亿吨，此与媒体报道的排沙量几乎一样。
24,0.0209t^2，0.4289t，10.631212*60*60*dte,
<二>对于第二个问题，两个阶段的数据如表三表四所示 0
表三:第一阶段试验数据
序号 1 2 3 4 5 6 7 8 9 10 11
水流量1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 x 含沙量32 60 75 85 90 98 100 102 108 112 115 h
表四:第二阶段的试验观测数据
序号 1 2 3 4 5 6 7 8 9 10 11 12 13
水流量2650 2600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900 x
含沙量116 118 120 118 105 80 60 50 40 32 20 8 5 h
对于第一阶段，有表四用MATLAB作图(如图三)可以看出其变化趋势，我们用多项式做最小二乘拟合。
设三次拟合函数关系h=a0+a1x+a2x^2+a3x^3
其中a0,a1,a2,a3,为待定系数。
四次拟合函数关系h= a0+a1x+a2x^2+a3x^3+a4x^4 其中a0,a1,a2,a3,a4为待定系数。
图三:第一阶段水流量与排沙量之间的关系图
三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数 h=0.0032x^2-2.4929x^3,拟合效果如图四所示
图四:三次多项式拟合效果，红线为拟合曲线
类似的四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合
函数
h=0.0121x^3-7.4347x^4,拟合效果如图五所示
图五:四次多项式拟合效果，蓝线线为拟合曲线
对于第二阶段，有表五用MATLAB作图可以看出其变化趋势，我们用多项式做
最小二乘拟合。设三次拟合函数关系h=a0+a1x+a2x^2+a3x^3 其中a0,a1,a2,a3,为待定系数。
四次拟合函数关系h= a0+a1x+a2x^2+a3x^3+a4x^4 其中a0,a1,a2,a3,a4为待
定系数。
三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601.则拟合函数 h=-0.9475x^2+464.9601x^3,拟合效果如图图六所示
图六:三次拟合函数拟合效果
类似的四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数
h=-0.0013x^2+1.1219x^3-354.5952x^4,拟合效果如图七所示
图七:四次拟合函数拟合效果
结论以及分析检验:
<一>用SAS软件做线性回归得到排沙量与时间的函数关系式为
,0.0209t^2，0.4289t，10.6312 y,e
再利用所求函数在区间[0,24]上进行积分得到总排沙量1.93962亿吨，这与现实情况基本相
符。
<二>对于第一阶段三次多项式拟合由MATLAB拟合函数求解
a0=a1=0,a2=0.0032,a3=-2.4929.
则拟合函数h=0.0032x^2-2.4929x^3
对于第一阶段四次多项式拟合由MATLAB拟合函数求解出
a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x^3-7.4347x^4 对于第二阶段三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-
0.9475,a3= 464.9601.则拟合函数h=-0.9475x^2+464.9601x^3,
对于第二阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-
0.0013，a3= 1.1219
a4=-354.5952则拟合函数h=-0.0013x^2+1.1219x^3-354.5952x^4
讨论与推广:1、对于第一个问题排沙量与时间不是严格的正态函数关系可能与实际有
些偏差，此外还可以用SAS软件进行高次的多向式回归
2、对于第二个问题，由于MATLAB软件的计算可能有些偏差导致拟合的函数关系可能与实
际有稍微偏差，此外，还可以进行高次的拟合。
附录:1、排沙量与时间的关系图像的MATLAB程序:
t=1:1:24;
y=[57600,114000,157500,187000,207000,235200,250000,265200,2862000,24 00,312800,307400,
306800,300000,271400,231000,160000,111000,91000,54000,45500,30000,80 00,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.4118 3、第一阶段的
排沙量与水流量之间的关系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 80
A1 =
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 80
A2 =
-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 80
A3 =
-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 80
A4 =
-0.0000 0.0000 -0.0013 1.1219 -354.5952 >> z4=polyval(A4,x1);
>> plot(x1,h1,'k*',x1,z4,'r:')
【参考文献】
【1】韩中庚，数学建模方法及其应用，北京，解放军信息工程程大学，2005
年6月，第一
版。
【2】盛祥耀，《高等数学》，北京，高等教育出版社，1992年8月，第一
版。
【3】张卓.SAS软件的应用，基于ARMA模型的商品销售额的预测分析，2005 【4】刘娜，在SAS中拟合ARCH/GARCH模型，2005