楊開(kāi)林，汪易森

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

1 前言

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

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

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

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

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

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

2.1 渠道沿程糙率的計(jì)算

對(duì)于恒定非均勻流渠道，渠道沿程糙率計(jì)算方程為

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

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

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

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

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

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

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

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

式(6)中

由式(1)和(6)可得

整理得

式(8)中，

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

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

其矩陣形式為

式(11)中，

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

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

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

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

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

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

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

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

3.1 實(shí)測(cè)數(shù)據(jù)

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

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

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

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

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

圖1 實(shí)測(cè)糙率n與水力半徑RFig.1 Calibrated channel roughness n，versus hydraulic radius R

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

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

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

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

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

或者

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)語(yǔ)

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

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

研究也表明，現(xiàn)有國(guó)內(nèi)常用沿程糙率經(jīng)驗(yàn)公式計(jì)算結(jié)果除一段渠道與實(shí)測(cè)值吻合外，其余渠道均遠(yuǎn)遠(yuǎn)小于實(shí)測(cè)值。本文沿程糙率經(jīng)驗(yàn)公式與美國(guó)墾務(wù)局經(jīng)驗(yàn)公式計(jì)算結(jié)果非常接近，可以作為其他人工混凝土渠道工程設(shè)計(jì)的依據(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)的水力計(jì)算及參數(shù)反問(wèn)題[J].水動(dòng)力學(xué)研究與進(jìn)展，1998，13(3):280－285.

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

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

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

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

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