黄河小浪底调水调沙问题

黄河小浪底调水调沙问题摘要：本文利用插值拟合的方法通过Matlab工具模拟出了排沙量与时间、排沙量与水流量的函数关系，并且求出了总排沙量为1.704亿吨。整个模型简单且方便计算，其中排沙量与水流量的函数关系为分段函数。关键词：调水调沙Matlab插值拟合1一、问题重述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.假设所给数据客观准确的反应了现实情况2.假设所给数据遵循一定规律变化，即是连续的3.假设模型中不需要考虑一些外在因素4.假设可将时间化为等分的时间点进行计算三、符号说明t:时间或时间点v:水流量S:含沙量V:排沙量2四、问题分析假设水流量和含沙量都是连续的，那么某一时刻的排沙量V=v(t)S(t)，其中v(t)为t时刻的水流量，而S(t)为t时刻的含沙量。通过观察数据，这些数据是每个12小时采集一次，所以我们可以将时间设为时间点t，依次为1，2，3，……，24，单位时间为12h。为了找到排沙量与时间的关系，我们就要先找到水流量和含沙量与时间的关系，一但找到水流量和含沙量与时间的关系，那么所要求的问题也就不难解决了。五、模型的建立与求解通过分析，我们假设水流量和含沙量都是连续的，那么我们开始对问题“(1)给出估算任意时刻的排沙量及总排沙量的方法”进行求解。我们通过Matlab工具将所知道的数据显示为直观的图像，如下所示，具体程序见附录的%tuxing.m。通过观察图像，我们可以看出其变化并不光滑，而且也没有特定的表现出服从某种分布的趋势。但是为了得到具体的计算函数，我们就必须对数据进行拟合，所以通过Matlab先利用spline方法对数据进行插值，从而提高精确度，使图像变得光滑，然后利用多项式进行拟合，当多项式次数越高拟合也越准确，但是由于数据受到的影响较多，所以这里的数据也不是准确值，因此我们可以只取三次进行拟合，也方便了后续的计算。于是我们分别对含沙量和水流量进行插值拟合，便可以得到下面图像和结果，具体程序见附录%hansha.m和%liuliang.m。3所得到的拟合函数为：y=0.014*x^{3}-1.3*x^{2}+21*x+16即含沙量与时间的关系式为：S=0.014*t^3-1.3*t^2+21*t+16所得到的拟合函数为：y=0.13*x^{3}-14*x^{2}+2.4e+002*x+1.5e+0034即水流量与时间的关系式为：v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003因为某一时刻的排沙量V=v(t)S(t)，所以我们可以将所拟合出来的多项式带入上式，通过Matlab进行计算可以得到下面答案，程序见附录%jisuan.m。ans=91/50000*t^6-73/200*t^5+2429/100*t^4-14573/25*t^3+2866*t^2+35340*t+24000即排沙量与时间的关系为：V=0.0018*t^6-0.365*t^5+24.29*t^4-582.92*t^3+2866*t^2+35340*t+24000由于这里的多项式次数过高，不便于计算和传播，所以我们可以对其再进行一次拟合，有下面结果，程序见附录%paisha.m。所以拟合后的函数为V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004，通过图像可以看出排沙量与时间服从正态分布，所以也可以化成的形式e的指数形式进行拟合，这里就不再重复计算。我们得到了拟合函数，下面就可以计算出这几天的总排沙量，通过Matlab编程可以计算出定积分，结果如下，程序详见附录%jisuan.m。ans=170366976000即总含沙量为1.704亿吨。下面我们对问题“(2)确定排沙量与水流量的变化关系。”进行分析计算。以下所有相关程序见附录%paishui.m，下面就不重复说明。我们先利用Matlab将排沙量和水流量的相关数据反映到图像中。5通过观察可以看出，其关系是分段的，所以我们按时间进行分段拟合，拟合原理同问题（1）相同，于是可以得到分段前后的拟合多项式。y=-7.5e-005*x^{3}+0.43*x^{2}-5.2e+002*x+3.6e+0046y=2.3e-005*x^{3}-0.066*x^{2}+1.9e+002*x-1.9e+005综上，我们可以得到排沙量与水流量的关系式为-7.5e-5*v^3+0.43*v^2-5.2e+2*v+3.6e+40=t9V=2.3e-5*v^3-0.066*v^2+1.9e+2*v-1.9e+59=t=24六、模型评估本模型的优点是：建模简单，方便计算，适用度广。但也有最大的缺点为：精确度较低。为了减少误差，我们可以通过增大模型中拟合多项式的次数。当然在日后的模型改进中可以加入误差评估系统，来对模型进行完善。附录%tuxing.mT=1:24;S=[32607585909810010210811211511611812011810580605030262085];W=[18001900210022002300240025002600265027002720265026002500230022002000185018201800175015001000900];7subplot(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.mT=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.mT=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)8%paisha.mt=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.mt=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('后半段图')