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