灰色数学黄河2.doc

灰色数学作业黄河小浪底调水调沙问题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) 确定排沙量与水流量的变化关系2、数据处理将数据处理成每12个小时为一个计算单位。任意时刻的排砂速度由该时刻的水流量与含砂量相乘得出。取每个时间单位的排砂速度平均值与单位时间乘积来计算每单位时间的排砂量，累加得到总排砂量，为便于计算， 6.29日以前排砂量约化为0。计算结果见表2。表2 排砂速度与累计排砂量计算值 日期6.296.307.17.2时间8:0020:008:0020:008:0020:008:0020:00序号12345678排砂速度57.6114157.5187207235.2250265.2总排砂量00.370660.957101.701222.552263.507414.555445.66827日期7.37.47.57.6时间8:0020:008:0020:008:0020:008:0020:00序号910111213141516排砂速度286.2302.4312.8307.4306.8300271.4231总排砂量6.859308.130679.4595010.7991412.1258113.436514.6707215.7559日期7.77.87.97.10时间8:0020:008:0020:008:0020:008:0020:00序号1718192021222324排砂速度160111915445.53084.5总排砂量16.600517.1858217.6221417.9353418.1502618.3133418.3954218.42242单位：排砂速度为吨每秒，总沙量为千万吨3、问题一、估算任意时刻的排沙量及总排沙量的方法3.1、任意时刻总排砂量与时间的关系表2显示出数据与时间之间有较好的规例性，通过建立DMG模型，直接对总排砂量进行模拟，结果显示平均相对误差为11%。通过对残差进行正弦曲线模拟（图1），修正预测结果后，结果较为满意（图2），见表3。起点1的值约化为0，DMG模型的起点不能为0，故以时间点2为起点计算。模拟的DMG参数a = 0.03421 b = 1.21364DMG 1-AGO模拟方程为：x1pa(K)=1098.8834* e(-0.03421*K) + 35.4768*K 1096.7951残差的模拟方程为：(K)=1.66*(SIN(3.14159*(K-11.78)/10.3)+0.73162总的还原方程为x0pa(K)= x1pa(K)-x1pa(K-1)+ (K), K=3, K24时为模拟值。表3 总排砂量的模拟结果序号K总排砂量模拟值残差(k)相对误差(k)20.3706560.3706630.9570960.953640.003453 0.36%41.7012161.70271-0.001492 0.09%52.5522562.519260.032995 1.29%63.5074083.432780.074631 2.13%74.555444.460260.095183 2.09%85.6682725.604780.063500 1.12%96.8592966.855300.004000 0.06%108.1306728.18784-0.057164 0.70%119.4595049.56775-0.108249 1.14%1210.79913610.95306-0.153925 1.43%1312.12580812.29840-0.172596 1.42%1413.43649613.55934-0.122847 0.91%1514.6707214.69654-0.025825 0.18%1615.75590415.679510.076394 0.48%1716.60046416.489470.110995 0.67%1817.18582417.121200.064620 0.38%1917.62214417.583610.038533 0.22%2017.93534417.8989560.036389 0.20%2118.15026418.100890.049370 0.27%2218.31334418.2314410.081903 0.45%2318.39542418.337160.058263 0.32%2418.42242418.46494-0.042517 0.23%精度检验为：平均相对误差(n-1) = 0.73% 一 级均方差比值 C= 0.0131 一级小误差概率 p=1.0 一级残差平方和 =0.1435这个结果是令人满意的（图2）。3.2、任意时刻排砂速度与时间的关系3.2.1、直接还原法由3.1总排砂量与时间的关系的计算公式及数据，直接对每一时刻的排砂速度进行还原，由于误差的累积效应，模拟精度不太好，结果下表。表4 由累计排砂总量对每时刻排砂速度的直接还原结果序号K排砂速度模拟值残差(k)相对误差(k)157.657.6002114114.000003157.5156.1041.396 0.887%4187190.673-3.673 1.964%5207187.34919.651 9.494%6235.2235.564-0.364 0.155%7250240.1109.890 3.956%8265.2289.748-24.548 9.256%9286.2289.189-2.989 1.044%10302.4327.719-25.319 8.373%11312.8311.1231.677 0.536%12307.4330.215-22.815 7.422%13306.8292.62213.778 4.497%14300291.1398.861 2.954%15271.4235.33736.063 13.288%16231219.73311.267 4.877%17160155.2424.758 2.974%18111137.223-26.223 23.624%199176.84914.151 15.550%205469.138-15.138 28.034%2145.524.34721.153 46.491%223036.087-6.087 20.290%23812.853-4.853 60.667%244.546.300-41.800 928.885%精度检验为： 扣除最后五个特别小而误差又特大的数：平均相对误差(n-1) 平均相对误差(n-1)= 52% 越界 = 6.16% 二 级均方差比值 均方差比值C= 0.1656 一级 C= 0.2000 一级小误差概率 小误差概率p=1.0 一级 p=1.0 一级残差平方和 =7366 残差平方和 =4882扣除最后几个数后（图3），误差仍然较大，可考虑其它方法。3.2.2、非线性还原法如果要进一步提高模拟的精度，可以考虑直接对排砂速度进行模拟。对表2计算的排砂速度分析，其1AGO有较好的规律性，用与上面同样的方法( DMG模型)，对排砂速度的1AGO模拟，结果显示相对误差最大时间点7（K=7）为2.41%，平均相对误差为0.8%，结果相当好。还原为各个计算点的速度后，计算其残差及相对误差，得到结果很差，个别相对误差达到320%，平均相对误差为24%，已超过模型得合理范围。原因在于原始数据的离散性太大，而且最后得两个极小值属于系统衰竭时的情况，由于数值太小，稍有误差就产生较大的还原变异。若扣除这两个特殊点，平均相对误差可以达到6.17%，与表-4直接还原的结果差不多。若对残差再次用正弦曲线分段模拟修正结果，扣除突变值，可得平均相对误差为1.78%。DMG模型排砂速度2-AGO模拟公式：x2pa(K+1)=260789.264* e(-0.03205*K) + 8417.155*K - 260731.663(K=1,2,23)排砂速度1-AGO残差修正公式：（K）=430*SIN(3.1416*(K-11.2)/10.9) + 183.851 排砂速度1-AGO还原公式x1pa(K+1)= x2pa(K+1)-x2pa(K)+（K）排砂速还原公式:排砂速度1-AGO还原公式x0pa(K+1)= x1pa(K+1)-x1pa(K)，结果见表5。表5-未经经残差修正排砂速度的模拟结果序号K排砂速度模拟值残差(k)相对误差(k)157.657.62114114.05391 -0.05391 0.05%3157.5160.03474 -2.53474 1.61%4187176.74970 10.25030 5.48%5207199.87688 7.12312 3.44%6235.2226.84608 8.35392 3.55%7250254.79105 -4.79105 1.92%8265.2280.78526 -15.58526 5.88%9286.2302.08272 -15.88272 5.55%10302.4316.34387 -13.94387 4.61%11312.8321.82801 -9.02801 2.89%12307.4317.53632 -10.13632 3.30%13306.8303.29369 3.50631 1.14%14300279.76263 20.23738 6.75%15271.4248.38806 23.01194 8.48%16231211.27767 19.72233 8.54%17160171.02746 -11.02746 6.89%18111130.50715 -19.50715 17.57%199192.62297 -1.62297 1.78%205460.07746 -6.07746 11.25%2145.535.14604 10.35396 22.76%223019.48894 10.51106 （35%）23814.01379 -6.01379 （75.2%）244.518.80034 -14.30034 （317.8%）精度检验为：（扣除最后三个误差极大值）平均相对误差(n-1) =6.17% 三 级均方差比值 C= 0.117 一级小误差概率 p0=1 一级残差平方和 =3118.093如表-5的残差项用正弦曲线分段拟合，对残差采用还原。如图4拟合曲线为分三段：C1(K)= 14*SIN(K-12)*3.1415926*2/10)-3 ,4K11C(K)= C2(K)= 25*SIN(K-12.5)*3.1415926*2/9.5-0.04)-3 ,11K16C3(K)= 15*SIN(K-19.5)*3.1415926*2/7+0.01)-2.5 ,K16公式为：x0pa(K+1)x0pa(K+1)+C(K)表6-经残差修正排砂速度的模拟结果序号K排砂速度模拟值残差(k)相对误差(k)157.657.62114114.05000 -0.053910.05%3157.5160.03000 -2.534741.61%4187187.06449 -0.064490.03%5207210.19167 -3.191671.54%6235.2232.07507 3.124931.33%7250251.79105 -1.791050.72%8265.2269.55627 -4.356271.64%9286.2285.76793 0.432070.15%10302.4300.02908 2.370920.78%11312.8310.59902 2.200980.70%12307.4305.47976 1.920240.62%13306.8307.45912 -0.659120.21%14300297.12824 2.871760.96%15271.4270.36530 1.034700.38%16231227.33316 3.666841.59%17160156.70705 3.292952.06%18111113.41734 -2.417342.18%199183.75018 7.249827.97%205450.93438 3.065625.68%2145.547.30261 -1.802613.96%223028.62231 1.377694.59%23811.36379 -3.36379(42%)244.54.47993 0.020070.45%精度检验为：平均相对误差(n-1) = 1.77% 二 级均方差比值 C= 0.0296 一级小误差概率 p=1 一级残差平方和=226.93模拟图见图5，可见，经修正后精度有所提高，与未修正的及与图-3相比结果有很大提高。由于整个系统在越过最大值后急速趋近于零，从数据对比上看，仍然显出最后几个数据偏差较大的问题。若将系统分成两端计算（K=18为界），建立GM(1,1)或DMG模型，经验算结果都差不多。3.2.3、收尾小数据处理最后一段系统趋近于衰竭，上面几种方法都显示模拟结果偏差较大，不宜做预测使用。如果对收尾段的小数据感兴趣，可对其单独建立模型。试算了几种模型(经修正后)，预测值都为负值，只有Verhulst模型预测结果为正值，而水的含砂量决不可能为负值，所以Verhulst模型较符合实际。模拟图见图-6。数据见表-7。1-AGO还原公式x1pa(K+1)= 99.98278/( 0.2895 + 0.611240*e(-0.90070*(K-18)原始数据累减还原为：x0pa（K+1）= x1pa(K+1)- x1pa(K)平均相对误差(n-1) =18.1% 三级均方差比 C= 0.265 一级小误差概率 p=1.0 一级对原始数据的1AGO一次GM(1,1)残差模型修正后再还原，平均相对误差(n-1) =8.0% 二级均方差比值 C= 0.119 一级小误差概率 p=1.0 一级对1AGO的残差修正方程为 （K+1）=16.347445*EXP(-0.449923*(K-18)-2.05124总的累减还原方程为：x0pa（K+1）= x1pa(K+1)- x1pa(K)+ （K+1）-（K）具体结果见下表：表7 Verhulst模型预测收尾数据序号181920212223242526原始数据111915445.53084.5Verhulst修正前74.899470.208446.417824.046510.89774.6381.9210.787修正后89.19647.53950.19530.8309.3624.3142.5451.1854、问题二、排砂速度与水流量的关系4.1、灰色关联度分析对表1的排砂速度及水流量进行排序，得表8得数据。表8-排砂速度及水流量排序序号123456789101112水流量90010001500175018001800182018501900200021002200排砂速度4.583045.55457.691111114160157.5187序号131415161718192021222324水流量220023002300240025002500260026002650265027002720排砂速度231207271.4235.2250300265.2306.8286.2307.4302.4312.8首先进行灰色关联度检验，先均值化，再零化，计算灰色绝对关联度为0.80；若先区间化，再零化，计算灰色绝对关联度为0.92；可见排砂速度与水流速度之间有较好的因果关系。由于曲线（图-7a）显出明显的洪峰前后同等流量排砂速度的差异，应此以洪峰为界将数据分成两段分别模拟，更符合实际情况。4.2、洪峰前流水量与排砂量关系4.2.1、线性回归拟合法图7a显示洪峰前流水量与排砂量有较好的线性关系，用表2中的洪峰前原始数据，按最小二乘法，进行一阶线性回归拟合：方程为x0(w (k)=0.25057w (k)373.3845，w (k)为水流量表9 线性拟合数据表序号K水流量排砂速度模拟值残差(k)相对误差(k)1180057.677.642-20.04234.79%21900114102.69911.3019.91%32100157.5152.8134.6872.98%42200187177.8709.1304.88%52300207202.9274.0731.97%62400235.2227.9847.2163.07%72500250253.041-3.0411.22%82600265.2278.098-12.8984.86%92700302.4303.155-0.7550.25%102720312.8308.1664.6341.48%模拟的平均相对误差为6%，去掉第一点平均相对误差为3%，有较好的精度。但考虑实际情况，流水速度与含砂量之间受各种环境因素的影响，关系复杂，局部简化为简单线性关系可以，推广的总体情况时，就显得有点过于简单化了。为此，根据已有洪峰前的数据特点，采用线性内插法，将数据化等（流量）时距简化处理，结果见表10。4.2.2、本例采用非线性拟合方法，具体如下：表10 洪峰前水流量与排砂量关系表序号12345678910水流量1800190020002100220023002400250026002700排砂速度57.6114135.75157.5187207235.2250265.2302.4分别用GM(1,1)模型、和DMG模型对数据进行拟合，结果都差不多，GM(1,1)模型的累减还原和差分还原，不修正时平均相对误差在4.5%（图7），一次残差还原平均相对误差在3.5%左右；用DMG模型模拟不修正时平均相对误差在8%，一次GM(1,1)残差还原平均相对误差下降到在2.16%。考虑简化处理，本例取GM(1,1)模型模拟。GM(1,1)模型还原方程为：1AGO： x1pa(k+1)=1113.768* e(0.1091*k) -1056.168累减还原方程：x0pa(k+1)=x1pa(k+1)-x1pa(k)具体见表-11 GM(1,1)模型累减还原数据（未经修正）。表-11 GM(1,1)模型累减还原数据序号K水流量排砂速度模拟值残差(k)相对误差(k)1180057.621900114128.37709 -14.37709 12.61%32000135.75143.17432 -7.42432 5.47%42100157.5159.67713 -2.17713 1.38%52200187178.08212 8.91788 4.77%62300207198.60854 8.39146 4.05%72400235.2221.50091 13.69909 5.82%82500250247.03194 2.96806 1.19%92600265.2275.50577 -10.30577 3.89%102700302307.26160 -5.26160 1.74%精度检验.平均相对误差(n-1) =4.5% 二级均方差比值C= 0.120 一级小误差概率 p=1.0 一级残差平方和=746.8744.3、洪峰后流水量与排砂量关系对与洪峰后的数据采用相同的数据处理方法，生成表12的数据。表12洪峰后水流量与排砂量关系表序号k12345678910水流量2700260025002400230022002100200019001800排砂速度310.83306.8300285.7271.4231195.5160127.354序号k111213141516171819水流量17001600150014001300120011001000900排砂速度42.3233.13025.521.216.812.484.5对表12的数据建立Verhulst模型，模拟效果见表13数据、图-8。Verhulst模型的方程为：x0pa(k+1)= 165.6027 / ( 0.51851 + 0.01427* exp(0.53277*k)总的还原方程为：x0pa(k+1)x0pa(k+1)+43.0533*exp(0.00684*(k-1)+29.388*exp(0. 1291e-4*(k-1)+67.514未经残差修正时，均方差比值最大C= 0.071，小误差概率 p=1.0，精度为一级，只有平均相对误差在32.3%，超出了模型有效范围。对残差建立GM(1,1) 残差二次修正后平均相对误差下降到12%，C 0.047 ，小误差概率 p=1.0，精度为一级。本数据序列的收尾段实测数据偏少，模拟的偏差较大，若扣除K=14以后几个非实测数据，则模拟数据未经修正的平均相对误差为11%，一次修正下降到6%，二次下降到4.8%。精度检验：平均相对误差(n-1) =4.8% 二级均方差比值 C= 0.041 一级小误差概率 p=1.0 一级Verhulst模型还原数据表：表13 二次修正Verhulst模型模拟