黄河小浪底调水调沙问题

黄河小浪底调水调沙问题主讲人：杜琴，陈晓梅，陈玲华，靳文进，方长箭内容摘要•为了确定排沙量与时间、排沙量与水流量的函数关系，我们可以用Matlab软件做线性回归得到排沙量与时间的函数关系式，再利用所求函数在区间[0,24]上进行积分得到总排沙量1.93962亿吨。对于排沙量与水流量之间的关系，按时间分为两段进行拟合，最终用Matlab软件来画出图像，确定排沙量与排水量之间的函数关系式。问题的提出•在小浪底水库蓄水后，黄河水利委员会进行了多次试验，特别是2004年6月到7月进行的黄河第三次调水调沙试验具有典型的意义。这次试验首次由小浪底、三门峡和万家寨三大水库联合调度，进行接力式防洪预泄放水，形成人造洪峰进行调水调沙试验成功。这次试验的一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底的库区的沉积泥沙。在小浪底开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2720m^3/s，使小浪底水库的排沙量也不断的增加。表一是由小浪底观测站从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含沙量326075909098100102108112115116日期7.57.67.77.87.97.10时间8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085问题分析：•1、对于问题一，所给数据中水流量x和含沙量h的乘积即为该时刻的排沙量y即：y=hx。•2对于问题二，研究排沙量与排水量的关系，从实验数据中可以看出，开始排沙量随水量增加而增加，而后随水流量的增加而减少，显然变化关系并非线性的关系，为此，把问题分为两部分，从水流量增加到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间的函数关系。模型假设•1、水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。•2、时间是连续变化的，所取时间点依次为1,2,3，…,24,单位时间为12h.模型的建立与求解•通过分析，我们假设水流量和含沙量都是连续的，那么我们开始对问题“(1)给出估算任意时刻的排沙量及总排沙量的方法”进行求解。•我们通过Matlab工具将所知道的数据显示为直观的图像，如下所示通过观察图像，我们可以看出其变化并不光滑，而且也没有特定的表现出服从某种分布的趋势。•图1•但是为了得到具体的计算函数，我们就必须对数据进行拟合，所以通过Matlab先利用spline方法对数据进行插值，从而提高精确度，使图像变得光滑，然后利用多项式进行拟合，当多项式次数越高拟合也越准确，但是由于数据受到的影响较多，所以这里的数据也不是准确值，因此我们可以只取三次进行拟合，也方便了后续的计算。于是我们分别对含沙量和水流量进行插值拟合，便可以得到下面图像和结果所得到的拟合函数为：y=0.014*x^{3}-1.3*x^{2}+21*x+16即含沙量与时间的关系式为：S=0.014*t^3-1.3*t^2+21*t+16•图2所得到的拟合函数为：y=0.13*x^{3}-14*x^{2}+2.4e+002*x+1.5e+003即水流量与时间的关系式为：v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003•图3因为某一时刻的排沙量V=v(t)S(t)，所以我们可以将所拟合出来的多项式带入上式，通过Matlab进行计算可以得到下面答案即排沙量与时间的关系为：V=0.0018*t^6-0.365*t^5+24.29*t^4-582.92*t^3+2866*t^2+35340*t+24000由于这里的多项式次数过高，不便于计算和传播，所以我们可以对其再进行一次拟合，有下面结果所以拟合后的函数为V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004，通过图像可以看出排沙量与时间服从正态分布，所以也可以化成的形式e的指数形式进行拟合•图4我们得到了拟合函数，下面就可以计算出这几天的总排沙量，通过Matlab编程可以计算出定积分，结果如下%jisuan.msymst;S=0.014*t^3-1.3*t^2+21*t+16;v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003;V=v*S;simple(V);symst;V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004;int(12*60*60*V,t,0,24)ans=170366976000即总含沙量为1.704亿吨•下面我们对问题“(2)确定排沙量与水流量的变化关系。”进行分析计算。•我们先利用Matlab将排沙量和水流量的相关数据反映到图像中。通过观察可以看出，其关系是分段的，所以我们按时间进行分段拟合，拟合原理同问题（1）相同，于是可以得到分段前后的拟合多项式•图5y=-7.5e-005*x^{3}+0.43*x^{2}-5.2e+002*x+3.6e+004•图6y=2.3e-005*x^{3}-0.066*x^{2}+1.9e+002*x-1.9e+005图7••综上，我们可以得到排沙量与水流量的关系式为•-7.5e-5*v^3+0.43*v^2-5.2e+2*v+3.6e+40=t9•V=•2.3e-5*v^3-0.066*v^2+1.9e+2*v-1.9e+59=t=24附件一•%tuxing.m•T=1:24;•S=[12011810580603260758590981001021081121151161185030262085];•W=[18001900210022002300240025002600265027002720265026002500230022002000185018201800175015001000900];•subplot(2,1,1);•plot(T,S);•holdon;•plot(T,S,'.');•title('时间与含沙量关系');•xlabel('时间t/12h');ylabel('含沙量/公斤每立方米');•subplot(2,1,2);•plot(T,W);•holdon;•plot(T,W,'.');•title('时间与水流量关系');•xlabel('时间t/12h');ylabel('水流量/立方米每秒');附件二•%hansha.m•T=1:24;•S=[32607585909810010210811211511611812011810580605030262085];•x=1:0.1:24;•y=interp1(T,S,x,'spline');•plot(T,S,'.',x,y);•title('时间与含沙量关系拟合图');•xlabel('时间t/12h');ylabel('含沙量/公斤每立方米');附件三•%liuliang.m•T=1:24;•W=[18001900210022002300240025002600265027002720265026002500230022002000185018201800175015001000900];x=1:0.1:24;y=interp1(T,W,x,'spline');plot(T,W,'.',x,y);title('时间与水流量关系拟合图');xlabel('时间t/12h');ylabel('水流量/立方米每秒');附件四%jisuan.msymst;S=0.014*t^3-1.3*t^2+21*t+16;v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003;V=v*S;simple(V);symst;V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004;int(12*60*60*V,t,0,24)ans=91/50000*t^6-73/200*t^5+2429/100*t^4-14573/25*t^3+2866*t^2+35340*t+24000附件五•%paisha.m•t=1:24;•V=0.0018*t.^6-0.365*t.^5+24.29*t.^4-582.92*t.^3+2866*t.^2+35340*t+24000;•plot(t,V);•title('时间与排沙量关系图')附件六•%paishui.m•t=1:24;•v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;•V=95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;•plot(v,V,'.');•title('整理图')附件七•figure;•t=1:9;•v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;•V=95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;•plot(v,V,'.');•title('前半段图')•figure;•t=10:24;•v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;•V=95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;•plot(v,V,'.');•title('后半段图')六、模型评估•本模型的优点是：建模简单，方便计算，适用度广。但也有最大的缺点为：精确度较低。为了减少误差，我们可以通过增大模型中拟合多项式的次数。当然在日后的模型改进中可以加入误差评估系统，来对模型进行完善。••谢谢大家！