南水北調(diào)中線工程渠道糙率的系統(tǒng)辨識

楊開林，汪易森

(1.中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038;2.國務(wù)院南水北調(diào)工程建設(shè)委員會專家委員會，北京 100038)

1 前言

南水北調(diào)中線工程是緩解我國北方水資源嚴重短缺、優(yōu)化水資源配置、改善生態(tài)環(huán)境的重大戰(zhàn)略性基礎(chǔ)設(shè)施，其渠道斷面尺寸大，形狀多樣，有梯形、矩形等，是寶貴的原型試驗平臺。目前，已經(jīng)對南水北調(diào)中線京石段過水建筑物典型斷面糙率進行了原型測試，按照常規(guī)水力學(xué)方法率定的糙率在0.013 37～0.015 7之間，大部分偏高，應(yīng)用價值受到懷疑。目前，中線工程正在建設(shè)，需要準(zhǔn)確可靠的實測糙率進行設(shè)計復(fù)核，同時，隨著工程的建成，其運行水力控制也需要這一關(guān)鍵水力學(xué)參數(shù)。因此，開展中線工程渠道糙率的系統(tǒng)辨識研究，不僅具有重要的理論價值，而且具有十分重要的實用價值。

系統(tǒng)辨識是20世紀60年代從現(xiàn)代系統(tǒng)理論發(fā)展出來的科學(xué)分支，國外對河道糙率的辨識研究開展較早。Becker等(1972)[1]首先構(gòu)造了一個關(guān)于明渠非恒定流參數(shù)辨識的影響系數(shù)方法，結(jié)合改進的單純形法，最小化誤差平方或最小化最大誤差絕對值來反演糙率;Chiu等(1978)[2]以卡爾曼濾波配合觀測的水深推求曼寧糙率;Wormleaton等(1984)[3]將水深和流量的相對誤差函數(shù)作為優(yōu)化函數(shù)，利用影響系數(shù)方法和非線性最小二乘算法求解最適合的糙率以縮減水深、流量的相對誤差;Crissman等(1993)[4]以圣維南方程為基礎(chǔ)建立了河床糙率隨時間變化的預(yù)測模式;Wasantha Lal(1995)[5]采用奇異值分解方法反演糙率;Khatibi等(1997)[6]探討了準(zhǔn)則函數(shù)在不同噪聲水平中的特性和樣本量對計算結(jié)果可靠性的影響;Atanov等(1999)[7]基于伴隨方程方法，采用最小二乘原理求解目標(biāo)泛函，反演計算梯型明渠糙率;Ramesh等(2000)[8]采用連續(xù)二次規(guī)劃方法反演糙率;Yan Ding等(2004)[9]基于最優(yōu)控制理論研究了二維淺水方程中糙率的辨識問題;Wong等(2003)[10]建立了水深及糙率值之間的關(guān)系，并在進行運動波模式演算過程中，由水深變化即時修正糙率;Hsu Minghis等(2006)[11]建立了實際觀測水深與計算值的目標(biāo)函數(shù)，利用Gauss－Newton方法求解非線性的最小二乘問題。

在國內(nèi)，金忠青等(1998)[12]根據(jù)實測資料，選擇水位過程或流量過程作為目標(biāo)函數(shù)的變量，構(gòu)造各河段誤差平方和這樣一個目標(biāo)函數(shù)，采用復(fù)合形法求解目標(biāo)函數(shù)直接優(yōu)化;董文軍等(2002)[13]根據(jù)參數(shù)辨識理論對一維水流方程中的糙率進行理論推導(dǎo)，使用最小二乘法建立最優(yōu)模型的目標(biāo)函數(shù);李光熾等(2003)[14]提出利用卡爾曼濾波來求解河道的糙率;程偉平等(2005)[15]引入控制論理論，應(yīng)用帶參數(shù)的卡爾曼濾波法進行河道糙率反演分析。

上述糙率的辨識研究主要針對河道，是在假設(shè)水力學(xué)參數(shù)，如流量、水深、水位等測量值準(zhǔn)確無誤差的基礎(chǔ)上進行的。由于原型水位、流量等測量數(shù)據(jù)包括很多不確定性因素，即隨機誤差，如測量儀器的誤差;測量環(huán)境的影響，如安裝儀器測點的選擇、安裝質(zhì)量、氣溫和水溫等因素，所以在一般情況下，按照這種方法獲取的糙率受測量誤差和計算誤差的影響，不能推廣應(yīng)用到其他渠道系統(tǒng)，特別是人工渠道。

本文的目的是根據(jù)水力學(xué)原理，考慮渠道斷面形狀、底坡、渠長變化的影響，建立渠道沿程糙率與粗糙高度ks和水力半徑R的函數(shù)關(guān)系，然后通過數(shù)學(xué)變換提出系統(tǒng)辨識的線性模型，最后以南水北調(diào)中線工程原型觀測資料為基礎(chǔ)，應(yīng)用系統(tǒng)辨識的方法消除水力測量隨機誤差的干擾，得到通用的渠道沿程糙率經(jīng)驗公式。

2 數(shù)學(xué)模型

2.1 渠道沿程糙率的計算

對于恒定非均勻流渠道，渠道沿程糙率計算方程為

式(1)中，Q為流量，m3/s;y為水深;A為斷面面積，m2;R為水力半徑，m;L為渠道長度，m;n為渠道糙率:s0為底坡;g為重力加速度，m/s2;下標(biāo)“1”為渠道進口;下標(biāo)“2”為渠道出口;ζ為渠道局部阻力系數(shù);ˉA為渠道過水?dāng)嗝嫫骄担?/p>

在工程設(shè)計中，國內(nèi)工程常常利用式(3)計算渠道沿程糙率，即

由此可以看出，糙率n不僅與反映渠道表面平整度的ks值有關(guān)，還與水力半徑R的大小有關(guān)，n將隨R的加大而加大。ks的取值在工程初期一般取0.000 61。

美國墾務(wù)局根據(jù)已建渠道的實測資料和室內(nèi)試驗分析，推薦采用下述方法確定混凝土渠道的糙率(王光謙等，2008)[16]

從式(3)和(4)可知，渠道糙率是等效粗糙高度ks和水力半徑R的函數(shù)，可以寫成統(tǒng)一形式

當(dāng)令系數(shù)c1=18和c1lgc2=19.55－18lgks時，則式(5)與國內(nèi)現(xiàn)有經(jīng)驗公式(3)相同。當(dāng)令系數(shù)c1=1/0.056 5=17.7 和 c1lgc2=17.7lg9 711=70.6時，則式(5)與美國墾務(wù)局經(jīng)驗公式(4)中的第二式相同。美國墾務(wù)局與我國經(jīng)驗公式系數(shù)的不同，主要反映了模型試驗和原型觀測的差別。南水北調(diào)中線工程干渠斷面尺寸在我國人工渠道中是最大的，根據(jù)原型實測觀測數(shù)據(jù)重新率定系數(shù)c1和c2的值，具有十分重要的實際意義。

2.2 渠道沿程糙率系統(tǒng)辨識線性模型

為了便于應(yīng)用系統(tǒng)辨識的最小二乘法確定系數(shù)c1和c2，式(5)可以改寫為

式(6)中

由式(1)和(6)可得

整理得

式(8)中，

在系統(tǒng)參數(shù)辨識過程中，線性方程式(8)的系數(shù)a1、a2和b是已知量，w1和w2是未知量，即需要辨識的參數(shù)。

對于m條渠道的系統(tǒng)，對渠道i式(8)可改寫為

其矩陣形式為

式(11)中，

2.3 系統(tǒng)辨識的最小二乘法

實測參數(shù)，如流量、水深，存在各種噪聲(誤差)，系數(shù)矩陣A和向量B是依賴于實測參數(shù)，所以按照式(11)只能得出辨識參數(shù)W的估計值，即

式(13)中，ε = [ε1，ε2，…，εM]T，其中上標(biāo)T為轉(zhuǎn)置符號。ε稱為擬合誤差，或者殘差。

式(16)中，矩陣ATA是對稱矩陣，當(dāng)它為非奇異(正則矩陣)時，有

上述結(jié)果是在認為觀測數(shù)據(jù)具有相同可信度的基礎(chǔ)上推導(dǎo)出來的，即對每個殘差εi給予相同的權(quán)，當(dāng)對每個殘差項加以不同的權(quán)，并令P為加權(quán)矩陣時，則廣義最小二乘法的準(zhǔn)則函數(shù)為

本文采用了全主元消去法直接求解式(16)和式(19)?？梢宰C明，最小二乘估計是無偏估計、一致估計和有效估計(徐枋同等，1999)[17]，可以有效消除隨機誤差的干擾。

3 渠道沿程糙率的系統(tǒng)辨識

下面將以南水北調(diào)中線京石段應(yīng)急供水工程實測資料為依據(jù)，考慮渠道斷面形狀、底坡、渠長變化的影響，應(yīng)用系統(tǒng)辨識的最小二乘法確定渠道沿程糙率與等效粗糙高度ks和水力半徑R的函數(shù)關(guān)系。

3.1 實測數(shù)據(jù)

表1列出了從唐河倒虹吸出口至漕河渡槽進口渠道特征參數(shù)一覽表，渠道數(shù)為8，其中渠道1～5是同形梯形渠道，它們的斷面形狀、尺寸、底坡完全相同;渠道6和8也是同形梯形渠道，它們的斷面形狀、尺寸、底坡完全相同;渠道7是隧洞，過水?dāng)嗝媸蔷匦巍Ｐ枰f明的是，渠道5和渠道6不是連續(xù)渠道，并且所有渠道中的流動都是非均勻流，水面線是雍水曲線。

表2列出各渠段平均水力半徑R和計算曼寧糙率n一覽表，其渠道編號與表1是一一對應(yīng)的，根據(jù)該表可以畫出如圖1所示n－R曲線。

從圖1可見，受糙率噪聲影響，糙率n隨水力半徑R變化波動較大，下面將應(yīng)用系統(tǒng)辨識的最小二乘法減少噪聲干擾，確定渠道沿程糙率與等效粗糙高度ks和水力半徑R的函數(shù)關(guān)系。

表1 渠道特征參數(shù)Table 1 Characteristic parameters of channels

表2 渠道參數(shù)Table 2 Parameters of channels

圖1 實測糙率n與水力半徑RFig.1 Calibrated channel roughness n，versus hydraulic radius R

3.2 渠道沿程糙率參數(shù)的最優(yōu)估計

當(dāng)取每座阻水橋局部阻力系數(shù)ζ為0.12時，根據(jù)式(7)和式(9)可得

把上述系數(shù)矩陣和向量代人式(16)，可得辨識參數(shù)W的最優(yōu)估計為W^LS=[ ]67.6 28.0T，因為w1=c1lgc2，w2=c1，代人式(5)可得

當(dāng)考慮每段渠道長度不同對每個殘差εi的影響時，可取加權(quán)矩陣P為對角線矩陣且對角線元素pii為

采用廣義最小二乘法參數(shù)估計，有

或者

需要說明，雖然式(20)和式(22)中沒有顯示出現(xiàn)等效粗糙高度ks，但是公式中已經(jīng)包含了ks的影響。

4 幾個重要渠道沿程糙率經(jīng)驗公式比較分析

目前，國內(nèi)外都有一些渠道沿程糙率計算經(jīng)驗公式，比較典型的是國內(nèi)常用經(jīng)驗公式(3)和美國墾務(wù)局經(jīng)驗公式(4)，下面根據(jù)南水北調(diào)中線京石段應(yīng)急供水工程實測資料把它們與本文兩個經(jīng)驗公式(20)和公式(22)進行比較分析。應(yīng)用上述4個經(jīng)驗公式，可以計算得到如表3所示結(jié)果，其中渠道編號與表1和表2一一對應(yīng)，實測沿程糙率是南水北調(diào)中線京石段應(yīng)急供水工程實測率定資料。

表3 典型經(jīng)驗公式渠道糙率計算比較表(按照水力半徑R的大小排列)Table 3 Comparison of channel roughness under four formulas(according to the hydraulic radius R)

根據(jù)表3數(shù)據(jù)可以畫出如圖2所示曲線，其中曲線1和2分別對應(yīng)本文經(jīng)驗公式(20)和(22);曲線3對應(yīng)國內(nèi)常用經(jīng)驗公式(3)，ks=0.000 61;曲線4對應(yīng)美國墾務(wù)局經(jīng)驗公式(4);點5是南水北調(diào)中線京石段應(yīng)急供水工程實測率定。

從圖2可得如下結(jié)論。

1)國內(nèi)常用沿程糙率經(jīng)驗公式計算結(jié)果，除一段渠道與實測值吻合外，其余渠道均遠遠小于實測值。國內(nèi)經(jīng)驗公式n的平均值約為0.013 6，實測平均值約為0.014 9，兩者平均相對偏差約為9.5%。造成誤差大的原因主要是該公式是根據(jù)實驗室資料得到的，并且ks=0.000 61與實際工程可能差別較大。

2)本文沿程糙率經(jīng)驗公式(20)和公式(22)與美國墾務(wù)局經(jīng)驗公式計算結(jié)果非常接近。從表3和圖2可見，本文沿程糙率經(jīng)驗公式(22)和美國墾務(wù)局經(jīng)驗公式對應(yīng)水力半徑R的n值偏差小于0.000 1，相對偏差小于0.7%。說明南水北調(diào)中線京石段應(yīng)急供水工程渠道建設(shè)質(zhì)量達到國外同等水平。

3)本文沿程糙率經(jīng)驗公式(22)計算南水北調(diào)中線京石段應(yīng)急供水工程n的平均值約為0.014 8，實測n的平均值約為0.014 9，兩者相對偏差小于0.7%。

綜上所述，本文經(jīng)驗公式(22)考慮了渠道長度的影響，可以作為人工混凝土渠道工程的計算依據(jù)。

圖2 渠道沿程糙率與水力半徑關(guān)系曲線Fig.2 Curves of channel roughnessversus hydraulic radius

5 結(jié)語

本文的創(chuàng)新點是提出了渠道沿程糙率的系統(tǒng)辨識模型，該模型考慮了渠道幾何參數(shù)，如斷面形態(tài)、長度、底坡等的影響，將渠道沿程糙率與粗糙高度ks和水力半徑R的關(guān)系用對數(shù)函數(shù)描述，依據(jù)水力學(xué)原理，然后通過數(shù)學(xué)變換提出了適合最小二乘法進行系統(tǒng)辨識的線性模型。

同時，以南水北調(diào)中線京石段應(yīng)急供水工程實測資料為依據(jù)，假設(shè)每座阻水橋局部阻力系數(shù)均為0.12，考慮渠道斷面形狀、底坡、渠長變化的影響，應(yīng)用最小二乘法得到的渠道沿程糙率計算公式為

研究也表明，現(xiàn)有國內(nèi)常用沿程糙率經(jīng)驗公式計算結(jié)果除一段渠道與實測值吻合外，其余渠道均遠遠小于實測值。本文沿程糙率經(jīng)驗公式與美國墾務(wù)局經(jīng)驗公式計算結(jié)果非常接近，可以作為其他人工混凝土渠道工程設(shè)計的依據(jù)。

[1]Becker L，Yeh W W G.Identification of parameters in unsteady open channel flows[J].Water Resources Research，1972，8(4):956－965.

[2]Chiu C L，Emmanuel O I.Kalman filtering in open channel flow estimation[J].Journal of the Hydraulics Division，1978，104(8):1137－1151.

[3]Wormleaton P R， Karmegam M.Parameter optimization in flood routing[J].Journal of Hydraulic Engineering，1983，110(12):1799－1814.

[4]Crissman R D，Chiu C L.Uncertainties in flow modeling and forecasting for Niagara River[J].Journal of Hydraulic Engineering，1993，119(11):1231－1250.

[5]Wasantha Lal A M.Calibration of riverbed roughness[J].Journal of Hydraulic Engineering，1995，121(9):664 －671.

[6]Khatibi R H，Williams J J R，Wormleaton P R.Identification problem of open － channel friction parameters[J].Journal of Hydraulic Engineering，1997，123(12):1078 －1088.

[7]AtanovG A， Evseeva E G，Meselhe E A.Estimation of roughness profile in trapezoidal open channel[J].Journal of Hydraulic Engineering，1999，125(3):309 －312.

[8]Ramesh R，Datta B，Bhallamudi S M，et al.Optimal estimation of roughness in open－channel flows[J].Journal of Hydraulic Engineering，2000，126(4):299 －303.

[9]Yan Ding， Jia Yafei，Sam S Y Wang，et al..Identification of manning’s roughness in shallow water flows[J].Journal of Hydraulic Engineering，2004，130(6):501 －510.

[10]WongT S W， Zhou M C.Kinematic wave parameters and time of travel in circular channel revisited[J].Advances in Water Resources，2003，26:417 －425.

[11]Hsu Minghis， Fu Jincheng，Liu Wencheng.Dynamic routing model with real－time roughness updating for flood forecasting[J].Journal of Hydraulic Engineering，2006，132(6):605 －619.

[12]金忠青，韓龍喜.復(fù)雜河網(wǎng)的水力計算及參數(shù)反問題[J].水動力學(xué)研究與進展，1998，13(3):280－285.

[13]董文軍， 楊則燊.一維圣維南方程的反問題研究與計算方法[J].水利學(xué)報，2002(9):61－65.

[14]李光熾，周晶晏，張貴壽.用卡爾曼濾波求解河道糙率參數(shù)反問題[J].河海大學(xué)學(xué)報，2003，31(5):490－493.

[15]程偉平，劉國華.基于廣義逆理論的河網(wǎng)糙率反演研究[J].浙江大學(xué)學(xué)報:工學(xué)版，2005，39(10):1603－1608.

[16]王光謙，黃躍飛，魏加華，等.南水北調(diào)中線工程總干渠糙率綜合論證[J].南水北調(diào)與水利科技，2006，4(1):8－14.

[17]徐枋同，李永華.系統(tǒng)辨識理論與實踐:在水電控制工程中的應(yīng)用[M].北京:中國電力出版社，1999.