1、水塔流量的估计一.问题的提出某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计其流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量。通常水泵每天供水一两次，每次约 2h(小时 )。水塔是一个高为 12.2m,直径为 17.4m 是正圆柱。按照设计,水塔水位降至约 8.2m 时,水泵自动启动,水位升到约为 10.8m 时水泵停止工作。表 1 是某一天的水位测量记录(符号“/”表示水泵启动)，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量及一天的总用水量。 表 1：水位测量记录（时
2、刻：h，水位：cm）时刻 0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97水位 968 948 931 913 898 881 896 852 839 822时刻 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93水位 / / 1082 1050 1021 994 965 941 918 892时刻 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91水位 866 843 822 / / 1059 1035 1011二、问题分析流量是单位时间流
3、出的水的体积，由于水塔是正圆柱形，横截面积是常数，在水泵不工作的时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。 水泵供水时段的流量只能靠供水时段前后的流量拟合得到，作为用于拟合的原始数据，我们希望水泵不工作的时段流量越准确越好。这些流量大体可由两种方法计算: 一是直接对表 1 中的水位用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量。二是先用表中数据拟合水位时间函数，求导数即可得到连续时间的流量。12.2m17.4m10.8m8.2m一般说来数值微分的精度不高，何况测量记录还是不等距的，数值微分的计算尤其麻烦。下面我们用第二种方法处理。有了任何时刻的
4、流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，如表 1 可知从t=0 到 t=8.97(h)水位下降了 968 822=146(cm)，乘以水塔的截面积就是这一时段的用水量。这个数值可以用来检查拟合的结果。3、模型假设1. 流量只取决于水位差，与水位本身无关。按照 Torricelli (托里切利, 1608-1647, 意大利数学家、物理学家、气压计原理发现者 )定律从小孔流出的流体的流速正比于水面高度的平方根，题目给出水塔的最低和最高水位分别是8.2m 和 10.8m(设出口的水位为零),因为 ,所以可忽略水位对10.8/2.51速度的影响。 2. 根
5、据最低和最高水位分别是 8.2m 和 10.8m 及表 1 的水位测量记录,假设水泵第 1 次供水时段为 到 ,第 2 次供水时段为 到 。9t1t 20.8t3t其中前 3 个时刻取自实测数据（精确到 0.1h）,最后 1 个时刻来自每次供水约两小时的已知条件（从记录看，每 2 次供水时段应在有记录的 22.96h 之后不久结束） 。3. 水泵工作时单位时间的供水量大致是常数，此常数大于单位时间的平均流量。4. 流量是对时间的连续函数。5. 流量与水泵是否工作无关。6. 由于水塔截面积是常数, ,为简单起22 217.4/3.78Srm见，计算中将流量定义为单位时间流出的水的高度，即水位对时
6、间变化率的绝对值（水位是下降的） ，最后给出结果时再乘以 S 即可。即：水位是时间的连续函数 ()ht水位对时间的变化率(流量) dt任何时刻的流量： ()VtSA四、模型建立1.拟合水位时间函数从表 1 测量记录看，一天有两个供水时段（以下称第 1 供水时段和第 2 供水时段）和 3 个水泵不工作时段(以下称第 1 用水时段 到 ，第 2 用0t8.97t水时段 到 第 3 用水时段 以后)。0.95t2.48t 23t对第 1、2 用水时段的测量数据分别作多项式拟合，得到水位函数 和 。为使拟合曲线比较光滑，多项式次数不要太高，一般1()ht2()ht用 36 次。由于第 3 时段只有 3
7、 个测量记录，无法对这一时段的水位作出比较好的拟合，可采用外推的办法解决。 2. 确定流量时间函数对于第 1、2 用水时段，只需将水位函数 求导数即可，(),12iiht对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第 2 供水时段流量外推，将第 3 用水时段流量包含在第 2 供水时段内，需要拟合四个流量函数。3. 一天总用水量的估计总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。 00t tVdSh五、模型求解根据表一，可以对各时段的数据进行拟合。建立时间和水位向量的函数关系。1212,th 为使拟合曲线比
8、较光滑，多项式次数不要太高，一般用 36 次。先用 3 次函数进行拟合第一时段0,9的水位、流量。设 t、h 为已输入的时刻和水位测量记录(水泵启动的 4 个时刻不输入): h=968 948 931 913 898 881 869 852 839 822 0 0 1082 1050 1021 994 965 941 918 892 866 843 822 0 0 1059 1035 1018; t=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.8
9、3 17.93 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91; x1=t(1:10);x2=x1.2;x3=x1.3; n=10;m=3; x=ones(n,1),x1,x2,x3; b,bint,r,rint,stats=regress(h(1:10),x,0.05)b =967.7356-22.10791.3586-0.0785bint =966.5367 968.9345-23.3536 -20.86221.0241 1.6932-0.1030 -0.0541r =0.2644-0.4851-0.16760.67540.0260-0.63
10、280.5374-0.46760.3089-0.0590rint =-0.3216 0.8505-1.5934 0.6232-1.3634 1.0282-0.2883 1.6392-1.2427 1.2947-1.7032 0.4376-0.5339 1.6088-1.5560 0.6207-0.8696 1.4874-0.6534 0.5353stats =1.0e+004 *0.0001 2.3222 0.0000 0.0000 rcoplot(r,rint)残差分析可以看到，拟合比较符合，拟合效果较好。经过实验，选取 3 次较合适。故以下选用三次函数进行拟合。1、拟合第 1 用水时段各时
11、刻的流量，可由如下程序代码得到： c1=polyfit(t(1:10),h(1:10),3); %用 3 次多项式拟合第 1 用水时段水位 h1=h1(t)，c1 输出 3 次多项式的系数a1=polyder(c1); %a1 输出多项式（系数为 c1）导数的系数，h1=-polyval(a1,t(1:10); %给出水位变化率 h1=h1(t)在 t(1)-t(10)上的离散值，即流量t1 = 0:0.1:9; %将第一用水时段0 , 9细分h11= -polyval(a1,tp1); %h11 输出多项式 a1 在 t11 点的函数值(取负后边为正值)，即t11 时刻的流量(水位下降的速率
12、) 。用程序进行计算，求得在0，9内各时刻的流量值(水位变化率)如下表：it0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97ih22.11 19.81 17.91 16.14 15.12 14.42 14.28 14.64 15.38 16.69第一用水段流量（水位变化率）曲线图2. 拟合第 2 用水时段11,20.8的水位、流量c2=polyfit(t(13:23),h(13:23),3); %用 3 次多项式拟合第 2 用水时段水位水位 h2=h2(t) ， c2 输出 3 次多项式的系数a2 = polyder(c2); %a2 输出多项式（系
13、数为 c2）导数的系数,h2=-polyval(a2,t(13:23); %给出水位变化率 h2=h2(t)t2 = 11:0.1:20.8; %将第二用水时段11 , 20.8细分h22 = -polyval(a2,t2); %x2 输出多项式（a2）在 tp2 点的函数值 (取负后边为正 值)，即tp2 时刻的流量(水位下降的速率) 。进行计算，求得在11,20.8内各时刻的流量值(水位变化率)如下表：it10.9512.0312.9513.8814.9815.9016.8317.9319.0419.9620.84ih22.1119.8117.9116.1415.1214.4214.281
14、4.6415.3816.6967.81第二用水段流量（水位变化率）曲线图3. 拟合第 1 供水时段9，11的流量在第 1 供水时段（t = 911）之前（即第 1 用水时段）和之后（第 2 用水时段）各取几点，其流量已经得到，用它们拟合第 1 供水时段的流量。为使流量函数在 t =9 和 t =11 连续我们简单地只取 4 个点，拟合 3 次多项式（即曲线必过这 4 个点） ，实现如下：x1 = -polyval(a1,8 9) %取第 1 时段在 t=8,9 的流量x2=-polyval(a2,11 12) %取第 2 时段在 t=11,12 的流量x12 = x1 x2; %将四个点合并c
15、12 = polyfit(8 9 11 12, x12, 3); %拟合 3 次多项式t12 =9:0.1:11; %将第一供水时段11 , 20.8 细分xx12 = polyval(c12,t12); %xx12 输出第一供水时段各时刻的流量进行计算，求得在9,11 内各时刻的流量值(水位变化率)如下表：it9.98 10.92ih25.39 32.45第一供水段流量（水位变化率）曲线图4. 拟合第 2 供水时段20.8,24的流量在第 2 供水时段之前取 t =20,20.8 两点的水流量，在该时刻之后（第 3 用水段）仅有 3 个水位记录，我们用差分得到流量，然后用这 4 个数值拟合第
16、 2 供水时段的流量如下： dt3 = diff(t(26:28); %最后 3 个时刻之后的两两之差dh3 = diff(h (26:28); %最后 3 个水位的两两之差dht3 = -dh3./dt3; % t(22)和 t(23)的流量（差商代替微商）t3 = 20 20.8 t(26) t(27)； %取第二供水时段前后各两点x3=-polyval(a2,t3(1:2),dht3; %求得 t3 各时刻的流量c3 = polyfit(t3,x3,3); % 拟合 3 次多项式t3 = 20.8:0.1:24； %将第二供水时段和第三用水时段细分xx3=polyval(c3,t3);
17、% xx3 输出第 2 供水时段及第三用水时刻的流量用程序进行计算，求得在20.8,24内各时刻的流量值(水位变化率)如下表：it22.01 22.96 23.88 24.99 25.91ih25.547 23.937 21.622 18.478 16.08第二供水段流量（水位变化率）曲线图4. 一天总用水量的估计 第 1、2 用水时段和第 1、2 供水时段流量的积分之和，就是一天总用水量。虽然诸时段的流量已表示为多项式函数，积分可以解析的算出，这里仍用数值积分计算。(1). 第一用水时段的用水量 1110001tttVdShdSh其中积分值 h1 通过梯形公式计算：10111()()()2k
18、NNt t kkttht (2). 第二用水时段的用水量 3332222tttVdShdSh321 1()()()kNNt t kkttht (3). 第一供水时段的用水量 22211112tttVdShdSh21121()()()kNNt t kkktthtt (4). 第二供水和第三用水时段的用水量 4443333tttVdShdSh01 13 ()()()2kNNt t kkttht m1 = 0.1*trapz(h1); % 第 1 用水时段用水量(按高度计), 0.1 为积分步长m2 = 0.1*trapz(h2); %第 2 用水时段用水量m12 = 0.1*trapz(xx12)
19、; %第 1 供水时段用水量m3 = 0.1*)trapz(xx3); %第 2 供水时段用水量V= (m1+m2+m12+m3)*237.8*0.01; %一天总用水量,高度是 cm，换算成 m。通过程序可以计算出各时段的用水量及一天的总用水量如下：用水量第 1 用水时段第 2 用水时段第 1 供水十点第 2 供水时段一天总用水量iV12V123V吨 347.62 613.63 119.85 178.14 1259.2其中四个时段都用三次多项式进行拟合，如果第一用水时段、第一、二供水时段用三次多项式，第二用水时段用四次多项式拟合，得另一计算结果：用水量第 1 用水时段第 2 用水时段第 1
20、供水十点第 2 供水时段一天总用水量iV12V123V吨 347.62 613.76 115.33 187.02 1263.7下图为一天用水流量的曲线图：5. 流量及总用水量的检验计算出的各时刻的流量可用水位记录的数据来检验。用水量 V1 用第 1 用水时段水位测量记录中下降高度 h1=968-822=146cm 来计算并检验，在第一用水时段水的实际用量为：V1=S* h1=237.8*146*0.01=347.19(吨)类似地，第二用水时段用去的水的高度 h2=1082-822=260cm,实际用水量为： V3=S* h2=237.8*260*0.01=618.28(吨)而计算值为： V1=
21、347.62(吨) ， V2=613.76(吨) 比较接近。供水时段流量的一种检验方法如下：供水时段的用水量加上水位上升值260 是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量） ，而两个供水时段水泵的功率应大致相等。第 1、2 时段水泵的公率可以计算如下：p1 = (V1+260*0.01*237.8)/2; %第 1 供水时段水泵的功率( 水量以高度计）tp2 = 20.8:0.1:23;xp2 = polyval(c3,tp2); % xp2 输出第 2 供水时段各时刻的流量p2 = (0.1*trapz(x3)+260*0.01)*237.8/2.2; %第 2 供水时段水泵的功率（水量仍以高度计）六.模型评价在实际中，各时段的函数拟合次数不一定相同，本题中为了方便，由第一用水时段的分析而推广至全部时段，应运了 3 次拟合。因此可以对此进行改进。本文分段处理，并应运了第 2 用水时段的外推，得到供水时段的用水量，方法简便，且有较好的实际效果，比较成功。
