黄河小浪底调水调沙问题

1黄河小浪底调水调沙问题摘要为了确定排沙量与时间、排沙量与水流量的函数关系，我们可以用SAS软件做线性回归得到排沙量与时间的函数关系式，再利用所求函数在区间[0,24]上进行积分得到总排沙量1.93962亿吨。对于排沙量与水流量之间的关系，按时间分为两段进行拟合，最终用MATHLAB软件来画出图像，确定排沙量与排水量之间的函数关系式。关键词：调水调沙实验，排沙量，matlab，拟合2目录摘要...............................1问题的提出..........................3问题分析............................3模型假设............................3模型的建立与求解....................4结论以及分析检验...................10讨论与推广.........................10参考文献...........................113一．问题的提出2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功．整个试验期为20多天，小浪底从6月19日开始预泄放水，直到7月13日恢复正常供水结束．小浪底水利工程按设计拦沙量为75.5亿立方米，在这之前，小浪底共积泥沙达14.15亿吨．这次调水调试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙．在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700立方米/每秒，使小浪底水库的排沙量也不断地增加．下面是由小浪底观测站从6月29日到7月10日检测到的试验数据：表1:试验观测数据单位：水流为立方米/秒，含沙量为公斤/立方米日期6.296.307.17.27.37.4时间8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10时间8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085现在，根据试验数据建立数学模型研究下面的问题：(1)给出估算任意时刻的排沙量及总排沙量的方法；(2)确定排沙量与水流量的变化关系。二．问题分析1、对于问题一，所给数据中水流量x和含沙量h的乘积即为该时刻的排沙量y即：y=hx。2、对于问题二，研究排沙量与排水量的关系，从实验数据中可以看出，开始排沙量随水量增加而增加，而后随水流量的增加而减少，显然变化关系并非线性的关系，为此，把问题分为两部分，从水流量增加到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间的函数关系。三．模型假设1、水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。2、时间是连续变化的，所取时间点依次为1,2,3，…,24,单位时间为12h.四．模型的建立与求解4已知给定的观测时刻是等间距的，以6月29日零时刻开始计时，则各次观测时刻分别为3600(124),1,2,...,24,itii其中：计时单位为s。第一次观察的时刻128800t，最后一次观察的时刻241022400t。记第i（1,2,...,24i）次观测时水流量为iv，含沙量为ic，则第i次观测时的排沙量为iiiycv。相关数据见表2表2插值数据对应关系单位：排沙量为kg/s日期6.296.307.17.27.37.4时间8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00节点123456789101112排沙量57600114000157500187000207000235200250000265200286200302400312800307400日期7.57.67.77.87.97.10时间8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00节点131415161718192021222324排沙量3068003000002714002310001600001110009100054000455003000080004500对于问题（1），根据所给问题的实验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，采用三次样条函数进行插值。利用MATLAB函数，求出三次样条函数，得到排沙量y=y（t）与时间的关系，然后进行积分，就可以得到总得排沙量21()ttzytdt。最后求得总得排沙量为91.84410t，计算的MATLAB程序如下：排沙量对时间的关系图像的MATLAB程序t=1:1:24;y=[57600,114000,157500,187000,207000,235200,250000,265200,286200,302400,312800,307400,306800,300000,271400,231000,160000,111000,91000,54000,45500,30000,8000,4500];plot(t,p,'b')5051015202500.511.522.533.5x105对于问题（2），研究排沙量与水流量的关系，从实验数据可以看出，开始排沙量是随着水流量的增加而增长，而后是随着水流量的减少而减少。显然，变化规律并非是线性的关系，为此，把问题分为两部分，从开始水流量增加到最大值27203m/s（即增长的过程）为第一阶段，从水流量的最大值到结束为第二阶段，分别来研究水流量与排水量的关系。画出排沙量与水流量的散点图两个阶段的数据如表3表4所示表3：第一阶段试验数据序号1234567891011水流量18001900210022002300240025002600265027002720含沙量326075859098100102108112115表4：第二阶段的试验观测数据序号12345678910111213水流量265026002500230022002000185018201800175015001000900含沙量11611812011810580605040322085排沙量与时间的关系图像的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')6051015202500.511.522.53x106第一阶段的排沙量与水流量之间的关系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:')71800190020002100220023002400250026002700280030405060708090100110120第一阶段三次多项式拟合函数以及拟合效果程序与结果：三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数h=0.0032x^2-2.4929x^3,A1=polyfit(x,h,3)z1=polyval(A1,x);plot(x,h,'k+',x,z1,'r')Warning:Polynomialisbadlyconditioned.AddpointswithdistinctXvalues,reducethedegreeofthepolynomial,ortrycenteringandscalingasdescribedinHELPPOLYFIT.Inpolyfitat76A1=1.0e+003*0.0000-0.00000.0032-2.492981800190020002100220023002400250026002700280030405060708090100110120第一阶段四次多项式拟合函数以及拟合效果程序与结果：四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4=-7.4347则拟合函数h=0.0121x^3-7.4347x^4A2=polyfit(x,h,4)z2=polyval(A2,x);plot(x,h,'*',x,z2,'r')Warning:Polynomialisbadlyconditioned.AddpointswithdistinctXvalues,reducethedegreeofthepolynomial,ortrycenteringandscalingasdescribedinHELPPOLYFIT.Inpolyfitat80A2=-0.00000.0000-0.00000.0121-7.434791800190020002100220023002400250026002700280030405060708090100110120第二阶段三次多项式拟合函数以及拟合效果程序与结果：三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3=464.9601.则拟合函数h=-0.9475x^2+464.9601x^3A3=polyfit(x1,h1,3)z3=polyval(A3,x1);plot(x,h,'*',x1,z3,'b')Warning:Polynomialisbadlyconditioned.AddpointswithdistinctXvalues,reducethedegreeofthepolynomial,ortrycenteringandscalingasdescribedinHELPPOLYFIT.Inpolyfitat80A3=-0.00000.0006-0.9475464.9601108001000120014001600180020002200240026002800020406080100120140第二阶段四次多项式拟合函数以及拟合效果程序与结果：四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3=1.1219a4=-354.5952则拟合函数h=-0.0013x^2+1.1219x^3-354.5952x^4A4=polyfit(x1,h1,4)z4=polyval(A4,x1);plot(x1,h1,'k*',x1,z4,'r:')Warning:Polynomialisbadlyconditioned.AddpointswithdistinctXvalues,reducethedegre