水塔流量估計講解

1、水塔流量估計摘要 1關(guān)鍵字 1一 .問題重述 1二 .問題分析 2三 .模型假設(shè) 3四問題求解 3（ 1 ）數(shù)據(jù)差值法 . 5（ 2 ）數(shù)值擬合法 . 8五 .結(jié)果簡單分析 12六模型檢驗 12七模型改進(jìn) 12八 .實(shí)驗小結(jié) 12九參考文獻(xiàn) 12附錄 12摘要針對水塔流量估計問題， 由于水塔流量與水位變化密切相關(guān)， 所以通過水位 變化來確定其流量； 用流速和時間的關(guān)系來確定水位的變化， 即水位差； 由物理 學(xué)中Torricelli定律確定流速，繼而解決具體問題。 我們對某一天水塔水位統(tǒng)計數(shù) 據(jù)進(jìn)行分析， 將其分為水泵供水時段和水泵不工作時段， 然后通過插值與擬合模 型分別進(jìn)行計算。計算結(jié)果如下

2、：圖 1.1 插值法算得各時段句全天的用水量（用高度計）第一未供水時 段第二未供水時 段第一供水時段混合時段全天拉格朗日插值 法145.6231258.866454.0689127.8244586.3828分段線性插值 法147.1430258.969749.6051115.0578570.7756三次樣條插值 法145.6870258.654753.3334120.0178577.6929圖 1.2 數(shù)據(jù)擬合法算得各個時段及全天的用水量（用高度計）時段第一未供水 時段第二未供水 時段第一供水時 段混合時段全天用水高度146.5150257.760546.131776.3076526.7148

3、關(guān)鍵字 ：流量 插值與擬合 torricelli 定律一 . 問題重述 :某一地區(qū)的用水管理機(jī)構(gòu)要求各社區(qū)提供以每小時多少加侖計的用水量以 及每天所用水的總量。 許多社區(qū)沒有測量流入或流出水塔的水量裝置， 他們只能 通過測量水塔每小時的水位高度來代替，其誤差不超過 0.05.但是，當(dāng)水塔中的 水位下降到最低水位 L 時水泵就自動啟動向水塔輸水直到最高水位 H，在此期間 無法測量水泵的供水量。 這樣，當(dāng)水泵正在輸水時就不容易建立水塔中水位和用 水量之間的關(guān)系，水泵每天輸水一次或兩次，每次約兩小時。當(dāng)水位降至約 27.00ft（英尺） 時，水泵開始啟動向水塔供水， 當(dāng)水位升到 35.50ft 時，

4、水泵自動停止工作。已知某一小鎮(zhèn)的水塔高度為 40ft，直徑為 57ft 的正圓柱， 下表記錄的是某一 天的水塔水位的真實(shí)數(shù)據(jù)（注： 1ft=0.3024m）。表 1.1 某一天水塔水位統(tǒng)計表時間 /s水位 /0.01ft時間 /s 水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵開動3

5、5932水泵開動82649水泵開動39332水泵開動859683475394353550899533397433183445932703340試估計任何時刻 （包括水泵正在輸水時間） 從水塔流出的水流量和一天的用 水總量。二. 問題分析：流量是單位時間內(nèi)流出水的體積，由于水塔是正圓柱形，橫截面積時常數(shù)， 所以在水泵不工作時段， 流量很容易根據(jù)水位相對時間的變化率算出。 問題的難 點(diǎn)在于如何估計水泵供水時段的流量。水泵供水時段的流量只能依據(jù)供水時段前后的流量經(jīng)差值或擬合得到。 作為 用于差值或擬合的原始數(shù)據(jù)， 我們希望水泵不工作時段的流量越準(zhǔn)確越好。 事實(shí) 上，水泵不工作時段的用水量可以由測量記

6、錄直接得到嗎， 由表中記錄的下降水 位乘以水塔的橫截面積就該時段的用水量。 這個數(shù)值可以用來檢驗數(shù)據(jù)差值或擬 合的結(jié)果。為了表示方便，我們將表中的數(shù)據(jù)全部化為國際標(biāo)準(zhǔn)單位（如下表） ，時間 用小時（ h），高度用米（ m）。表 2.1 一天內(nèi)水塔水位記錄時間（ h）水位（ m）時間（ h）水位（ m）09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.848.

7、228.978.2222.02水泵開動29.98水泵開動22.96水泵開動10.93水泵開動23.8810.5910.9510.8224.9910.3512.0310.5025.9110.18三. 模型假設(shè)（1）假設(shè)流量只取決于水位差，與水位本身無關(guān)。由水塔的最低和最高水位分 別是 8.1648m（27x0.3024 ）和 10.7352m（35.50x0.3024 ），設(shè)出口的水位為零， 從小孔流出的流速正比于水面高度的平方根，可忽略水位對流速的影響。（2）假設(shè)流量是時間的連續(xù)函數(shù)。為計算簡單，不妨將流量定義成單位時間流 出水的高度，即水位對時間變化率的絕對值（水位時下降的） ，得到結(jié)果后在

8、乘 以水塔橫截面積即可。（3）已知水塔橫截面積為 S=（57*0.3048）2* /4=237.8（平方米）。四問題求解首先根據(jù)表 2.1的數(shù)據(jù)，用 MATLAB軟件作出水位 時間散點(diǎn)圖 4.1。程序見 附錄 1圖 4.1 水位 時間散點(diǎn)圖下面計算流量與時間的關(guān)系。 根據(jù)數(shù)據(jù)散點(diǎn)圖 1.0,一種簡單的處理方法是先 將表 2.0中的數(shù)據(jù)分為三段， 然后對每一段數(shù)據(jù)做如下處理： 設(shè)某段數(shù)句為 （x0， y0），（x1,y1）,（xn，yn），相鄰數(shù)據(jù)中點(diǎn)的平均流速采用公式：流速 =（左端點(diǎn)的水位 -右端點(diǎn)的水位） /區(qū)間長度，即： x（i 1） xiv（）=（yi-y（i+1）/（x （i+1）-

9、xi）每段數(shù)據(jù)首尾點(diǎn)的柳樹采用下面的公式計算：v(x0)=(3y0-4y1+y2)/(x2-x0),v(xn)=(-3yn+4y(n-1)-y(n-2).用以上公式算得流速與時間之間的數(shù)據(jù)表 4.0 如下：表 4.1 流速與時間關(guān)系數(shù)據(jù)表時間( h)流速( cm/h )時間( h)流速( cm/h )029.8913.4229.030.4621.7414.4326.361.3818.4815.4426.092.39516.2216.3724.733.4116.3017.3823.644.42515.3218.4923.425.4413.0419.5025.006.4515.4520.4023.

10、867.46513.9820.8422.178.4516.3522.02水泵開動8.9719.2922.96水泵開動9.98水泵開動23.8827.0910.93水泵開動24.4321.6210.9533.5025.4518.4811.4929.6325.9113.3012.4931.52由表 4.1作出流速 時間散點(diǎn)圖 4.2，程序見附錄 2：圖 4.2 流速 時間關(guān)系散點(diǎn)圖 下面用數(shù)據(jù)差值和數(shù)據(jù)擬合兩種方法來估計水塔水流量：（1）數(shù)據(jù)差值法由表 4.1，對水泵不工作時段 1和水泵不工作時段 2 采用差值方法。 可以得 到任意時刻的流速， 從而可以知道任意時刻的流量。 對于水泵工作時段數(shù)據(jù)太

11、少， 我們將它與水泵工作時段 2 合并一同進(jìn)行差值。由于最后一段水泵不工作時段數(shù) 據(jù)太少，我們將它與水泵工作時段 2 合并一同進(jìn)行差值處理 （該時段一下簡稱混 合時段）。這樣，總共需要對四段數(shù)據(jù)（第 1，2 未供水時段，第 1供水時段，混 合時段）進(jìn)行差值處理。A. 第 1 未供水時段數(shù)據(jù)差值求解， 分別用拉格朗日插值法，分段線性差值及三次 樣條差值三種方法算出流量函數(shù)和用水量（用水高度） 。圖 4.3 中曲線 lglr， fdxx，scyt 分別表示用拉格朗日差值法，分段線性差值法 及是三次樣條差值法對第 1 未供水時段數(shù)據(jù)差值得到的流速曲線， 程序見附錄 3圖 4.3 第 1 未供水時段流

12、速曲線B. 第 2 未供水時段數(shù)據(jù)差值求解，分別用拉格朗日插值法，分段線性差值及三 次樣條差值三種方法算出流量函數(shù)和用水量（用水高度） 。圖 4.4 中曲線 lglr， fdxx，scyt 分別表示用拉格朗日差值法，分段線性差值法 及是那次樣條差值法對第 2 未供水時段數(shù)據(jù)差值得到的流速曲線， 程序代碼見 附錄 4：圖 4.4 第 2 未供水時段流速曲線C. 第 1 供水時段數(shù)據(jù)差值求解， 分別用拉格朗日插值法， 分段線性差值及三次樣 條差值三種方法算出流量函數(shù)和用水量（用水高度） 。圖 4.5 中曲線 lglr ， fdxx ，scyt 分別表示用拉格朗日差值法，分段線性差 值法及是那次樣條

13、差值法對第 2 未供水時段數(shù)據(jù)差值得到的流速曲線， 程序代碼 見附錄 5：圖 4.5 第 2 未供水時段流速曲線D. 混合時段數(shù)據(jù)差值求解，分別用拉格朗日插值法，分段線性差值及三次樣條差 值三種方法算出流量函數(shù)和用水量（用水高度） 。圖 4.2 中曲線 lglr， fdxx，scyt 分別表示用拉格朗日差值法，分段線性差值法 及是那次樣條差值法對第 2 未供水時段數(shù)據(jù)差值得到的流速曲線 : 程序執(zhí)行后顯示結(jié)果，程序代碼見附錄 6：圖 4.6 混合時段流速曲線綜合得到：表 4.2 差值法算得個時段及全天的用水總量（以高度計）時段及方法第 1 未供水 時段第 2 未供水 時段第 1 供水時 段混合

14、時段全天拉格朗日插 值法145.6231258.866454.0689127.8244586.3828分段線性差 值法147.1430258.969749.6051115.0578570.7756三次樣條差 值法145.6870258.654753.3334120.0178577.6929圖4.6是用分段線性及三次樣條差值方法算得全天水塔水流速與時間的關(guān)系 曲線，圖中曲線 scyt 表示三次樣條差值曲線，曲線 fdxx 表示分段線性差值的曲 線，程序代碼見附錄 7：圖 4.7 全天水塔流速與時間的關(guān)系曲線(2)數(shù)值擬合法1) 擬合水位 時間函數(shù)根據(jù)表 4.1 中的數(shù)據(jù)，分別對第 1、 2 未供

15、水時段的測量數(shù)據(jù)直接做三次多 項式擬合，得到水位與實(shí)踐的關(guān)系函數(shù)。程序中，變量 x1 存放了以 0.1 為時間步長第 1 未供水時段各個時刻的水位 高度。圖 4.7 繪出的是第 1 未供水時段水位與時間的三次多項式擬合曲線， 程序 代碼見附錄 8：圖 4.8 第 1 未供水時段水位與時間的三次多項式擬合曲線下圖4.9 ，Matlab程序?qū)Φ?2未供水時段的數(shù)據(jù)進(jìn)行擬合 ,得到的第2未供水時 段水位與時間三次多項式擬合曲線， Matlab程序見附錄 9.圖 4.9 第 2 未供水時段水位與時間的三次多項式擬合曲線2）確定流量 時間函數(shù)對第 1,2未供水時段的水位求導(dǎo)可得流量， 用三次多項式擬合第

16、 1,2 未供水 時段的流速與時間關(guān)系曲線，程序代碼見附錄 10：圖 4.10 第 1,2 未供水時段的流速與時間關(guān)系曲線（三次多項式擬合）與圖 4.6 中相應(yīng)部分段對比發(fā)現(xiàn)，用三次多項式擬合效果不是很好，若改用5 次多項式擬合則得效果較好的第 1,2 未供水時段的流速與時間關(guān)系曲線，見圖 4.9（程序代碼見附錄 11）:圖 4.11 第 1,2 未供水時段的流速與時間關(guān)系曲線（五次多項式擬合）第 1 供水時段的流速則用前后時刻的流速擬合得到。 為使流速函數(shù)在 t=9 和 t=11 處連續(xù)，只取 4 個點(diǎn)，用 3 次多項式擬合得第 1 供水時段流速與時間關(guān)系 曲線圖 4.12，同圖 4.1 中

17、相應(yīng)部分段較為吻合，程序代碼見附錄 12：圖 4.12 第一供水時段流速與時間的三次多項式似合曲線10對于混合時段，在第 2 供水時段之前取 t=19.96,20.84 兩點(diǎn)的流速，用第三 末供水時段的 3 個記錄做差分得到兩個流速數(shù)據(jù) 21.62,18.48，然后用這 4 個數(shù)據(jù) 做 3 次多項式似合得到混合時段的流速與時間的關(guān)系曲線如圖 4.13 所示，同數(shù) 據(jù)插值法得到的曲線圖 4.1 中相應(yīng)曲線段較為吻合，程序代碼見附錄 13：圖 4. 11 混合時段流速與時間的三次多項式擬合曲線3）估計一天的用水總量分別對供水的兩個時段和不供水的兩個時段積分 （流量對時間積分） 并求和 得到一天的用

18、水總量約為 526.8935（此數(shù)據(jù)是用水總高度 cm） .表 4.3.1 列出各 段用水量，同數(shù)據(jù)插值法算得的表 4.1 數(shù)據(jù)相比，較為吻合。表 4.3 數(shù)據(jù)擬合算得各個時段的用水高度（一）時段第 1 未供水 時段第 2 未供水 時段第 1 供水時 段混合時段全天用水高度145.67260.6646.6073.9635526.89354）流量及用水總量的檢驗 計算出各時刻的流量可以用水位記錄數(shù)據(jù)的數(shù)值微分來檢驗， 各時段的用水 高度可以用實(shí)際記錄的水位下降高度來檢驗。 例如，算得第 1 未供水時段的用水 高度是 145.67，而實(shí)際記錄的水位下降高度為 968-822=146，兩者是吻合的同

19、樣 地，算得第 2 未供水時段的用水高度是 260.66，而實(shí)際記錄的水位下降高度為 1082-822=260，兩者也是吻合的。從算法設(shè)計和分析可知，計算結(jié)果與各時段所用的擬合多項式的次數(shù)有關(guān)。表 4.4 給出的是對第 1 ，2 未供水時段分別用 5 ，6 次多項式擬合得到的用水高度。11表 4.4 數(shù)據(jù)擬合法算得各個時段的用水高度（二）時段第 1 未供水 時段第2 未供水 時段第 1 供水時 段混合時段全天用水高度146.5150257.760546.131776.3076526.7148五. 結(jié)果簡單分析對于數(shù)據(jù)差值法，由表 4.1 可以看出，使用三次樣條差值法得到的第 1,2 未 供水時

20、段的用水高度結(jié)果 145.6870和258.6547與表中記錄的下降高度 146和 260 相差不大， 說明差值結(jié)果與原始數(shù)據(jù)比較吻合。 用三次樣條差值法估計全天的用 水量約為 539.4450x237.8x10=1282800.21（升）。對 于 數(shù) 據(jù) 擬 合 法 ， 由 表 4.3.2 可 得 全天 的 用 水 總 量約 為526.7148x237.8x10=1252527.7944（升），同數(shù)據(jù)差值法得到的結(jié)果也很接近。六模型檢驗七模型改進(jìn)八. 實(shí)驗小結(jié)本次實(shí)驗主要進(jìn)行水塔水量的估算。 第一種估計為數(shù)據(jù)差值方法， 我們用了 三種不同的差值方法進(jìn)行估計。 在求解的過程中， 可以熟悉數(shù)據(jù)差

21、值的理論和方 法。第二種估計方法為數(shù)據(jù)擬合法， 用多項式進(jìn)行擬合，得到水塔水流量的估 計。九參考文獻(xiàn)附錄：下面程序是對上面問題的補(bǔ)充，均由 matlab 軟件編寫：附錄 1：x=xlsread( stgj.xls, sheet1 , A2:A29 );y=xlsread( stgj.xls, sheet1 , B2:B29 );plot(x,y, b* )xlabel( x 軸 );ylabel(y 軸 );附錄 2：x=xlsread(stgjls.xls, sheet1, A2:A32 );y=xlsread(stgjls.xls, sheet1, B2:B32 );plot(x,y,b*

22、 )xlabel( x軸 );ylabel(y 軸 );12附錄補(bǔ)充 ： 拉格朗日函數(shù)文件f unction y=lglrcz(x0,y0,x)n=length(x0); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j=kp=p*(z-x0(j)/(x0(k)-x0(j); endends=p*y0(k)+s;end=s;附錄 3：t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97;v=29.89,21.74,18.48,16.22,16.30,15.32

23、,13.04,15.45,13.98,16.35,19.27; t0=0:0.1:8.97;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,spline);sancytjf=0.1*trapz(scyt)plot(t,v,* ,t0,lglr,r,t0,fdxx,g,t0,scyt,b) gtext(lglr)gtext(fdxx)gtext(scyt)13附錄 4t=10.95,11.49,12.49,13.42,14.43,15.

24、44,16.37,17.38,18.49,19.50,20.4 0,20.84;v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.8 6,22.17;t0=10.95:0.1:20.84; lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr) fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx) scyt=interp1(t,v,t0, spline );sancytjf=0.1*trapz(scyt)plot(t,v, * ,t0,lglr, r

25、,t0,fdxx, g ,t0,scyt, b) gtext( lglr )gtext( fdxx ) gtext( scyt )附錄 5：t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84;v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,33.50,29.63,31.52,29.03,26.36,26.

26、09,24.73,23.64,23.42,25.00,23.86,22.17;t0=8.97:0.1:10.95;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0, spline );sancytjf=0.1*trapz(scyt)plot(t,v, * ,t0,lglr, r ,t0,fdxx, g ,t0,scyt, b ) gtext( lglr ) gtext( fdxx ) gtext( scyt )14附錄 6：t=10.9

27、5,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84,23.88,24.43,25.45,25.91;v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17,27.09,21.62,18.48,13.30;t0=20.84:0.1:22.96;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=i

28、nterp1(t,v,t0, spline );sancytjf=0.1*trapz(scyt)plot(t,v, *,t0,lglr, r,t0,fdxx, g,t0,scyt, b ) gtext( lglr ) gtext( fdxx )附錄 7：t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49 ,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84,23.88,24 .43,25.45,25.91;v=29.89,21.74,18.4

29、8,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29, 33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22. 17,27.09,21.62,18.48,13.30;t0=0:0.1:25.91; fdxx=interp1(t,v,t0);fdxxjf=trapz(fdxx) scyt=interp1(t,v,t0, spline );sancytjf=trapz(scyt)plot(t,v, * ,t0,fdxx, g ,t0,scyt, b ) gtext( f

30、dxx ) gtext( scyt )附錄 8：t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:1

31、0),h(1:10),3); % 輸入項 3 表示擬合多項式的次數(shù)tp1=0:0.1:8.9;x1=polyval(c1,tp1);plot(tp1,x1)15附錄 9：t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.4

32、1,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),5); % 輸入項 5 表示擬合多項式的次數(shù)tp1=0:0.1:8.9;x1=polyval(c1,tp1);plot(tp1,x1)附錄 10：c1=polyfit(t(1:10),h(1:10),3); % c2=polyfit(t(11:21),h(11:21),3); a1=polyder(c1);a2=polyder(c2); tp1=0:0.01:8.97;tp2=10.95:0.01:20.84; x13=-polyval(a1,tp1)

33、;x113=-polyval(a1,0:0.01:8.97); wgsysl1=100*trapz(tp1,x113); % x14=-polyval(a1,7.93,8.97); % x23=-polyval(a2,tp2);x114=-polyval(a2,10.95:0.01:20.84); wgsysl2=100*trapz(tp2,x114); % x24=-polyval(a2,10.95,12.03); % x25=-polyval(a2,19.96,20.84); % subplot(1,2,1)plot(tp1,x13*100) subplot(1,2,2) plot(tp2,x23*