于顯亮，彭 楊，吳志毅

(華北電力大學(xué) 可再生能源學(xué)院，北京102206)

0 引 言

數(shù)值模擬是研究水流運(yùn)動(dòng)的主要手段，模型中參數(shù)取值是否合理對(duì)計(jì)算結(jié)果影響很大，其中糙率是一個(gè)相當(dāng)重要又敏感的關(guān)鍵參數(shù)，對(duì)河道水流及其沖淤變化的計(jì)算成果有很大影響。糙率的確定是一個(gè)復(fù)雜的問題，目前尚無通用的方法，一般根據(jù)經(jīng)驗(yàn)公式估算或采用試錯(cuò)法反求不同水位下的糙率值，具有較大的經(jīng)驗(yàn)性和任意性。尤其在河網(wǎng)計(jì)算中，手工調(diào)試費(fèi)時(shí)費(fèi)力、工作量大，精度難以保證，且當(dāng)計(jì)算條件改變時(shí)，糙率也會(huì)隨之變化，需要經(jīng)常更新。20世紀(jì)60年代以來隨著反演問題研究的發(fā)展，國內(nèi)外很多學(xué)者對(duì)糙率反問題也進(jìn)行了研究。Becker和Yeh[1]提出影響系數(shù)法對(duì)糙率進(jìn)行反演計(jì)算。Wasantha Lal[2]采用奇異值分解方法反演糙率。金鐘青，韓龍喜等[3]采用復(fù)合形法對(duì)糙率反問題進(jìn)行了求解。程偉平[4]采用卡爾曼濾波對(duì)河道糙率參數(shù)反問題進(jìn)行求解。張潮[5]采用BP神經(jīng)網(wǎng)絡(luò)獲取“水位觀測值-糙率”之間的映射關(guān)系，將水位實(shí)測值代入該映射關(guān)系，得到河網(wǎng)各河段的糙率反演值。劉志賢[6]采用基本遺傳算法對(duì)河網(wǎng)糙率進(jìn)行了反分析研究。李麗等[7]采用自適應(yīng)隨機(jī)搜索算法進(jìn)行糙率反演。這些研究將河道糙率視為定值，沒有考慮糙率隨流態(tài)的變化，且有些算法或收斂速度受約束條件數(shù)目影響較大，或受初值影響較大，容易早熟或陷入局部最優(yōu)，或?qū)崟r(shí)數(shù)據(jù)要求很高，需要有足夠數(shù)量的水位觀測資料。因此，有必要對(duì)河道糙率反演方法進(jìn)行進(jìn)一步研究。

對(duì)于天然河道，當(dāng)流量水位變幅較大時(shí)，糙率一般隨流量或水位變化[8]，計(jì)算時(shí)若同一個(gè)計(jì)算斷面在不同流量下取同一個(gè)糙率值無法反映水流真實(shí)的流動(dòng)過程，因此本文將糙率視為流量的函數(shù)，將河道不同流量級(jí)下的糙率反演問題轉(zhuǎn)化為一個(gè)多階段的糙率優(yōu)化問題。以河道一維恒定流計(jì)算為例，將水位誤差平方和最小作為目標(biāo)，提出了基于動(dòng)態(tài)規(guī)劃算法的河道糙率反演模型，將糙率試錯(cuò)的過程交由電腦來完成，實(shí)現(xiàn)河道糙率的自動(dòng)反演。

1 模型建立

1.1 一維恒定流模型

對(duì)于長河段長時(shí)段的非恒定水流運(yùn)動(dòng)，可將其概化為平均的洪水過程，作為恒定問題處理，不考慮流量沿程變化，此時(shí)水流運(yùn)動(dòng)可直接用運(yùn)動(dòng)方程來描述，即

(1)

式中：x為流程；Q為流量；A為過水?dāng)嗝婷娣e；Z為水位；R為水力半徑；n為糙率；g為重力加速度。

采用有限差分法對(duì)式(1)進(jìn)行求解，其離散格式為：

(2)

式中：Δxj為斷面間距，Δxj=xj-xj+1,其中xj為斷面距下游斷面的距離；ψ為數(shù)值離散權(quán)重因子，本文中ψ=0.5；角標(biāo)j和j+1表示變量為計(jì)算斷面j和j+1上的變量。

考慮到天然河道中的水位一般為緩流，且出口水流條件不同，其水面線呈壅水或降水曲線，計(jì)算時(shí)一般從下游斷面向上游斷面進(jìn)行推算。假定一系列的下游斷面水位zj+1，從中求出滿足式(2)的zj值，即為所求的上游斷面水位。式(1)中的糙率n是一待定系數(shù)，其值直接影響計(jì)算斷面水位的大小，必須慎重取值，不能簡單地選用一個(gè)固定值。本文不僅考慮糙率沿程變化，而且也考慮糙率隨流量的變化，將基于動(dòng)態(tài)規(guī)劃算法建立河道糙率反演模型，將河道不同流量級(jí)下的糙率反演問題轉(zhuǎn)化為一個(gè)多階段的糙率優(yōu)化問題，實(shí)現(xiàn)對(duì)糙率參數(shù)的自動(dòng)調(diào)整。

1.2 河道糙率反演模型

(3)

式中：Zi,t為t時(shí)刻通過上述模型計(jì)算得到的第i個(gè)水位站的水位值；Z′i,t為t時(shí)刻第i個(gè)水位站的實(shí)測水位值；T和M分別為計(jì)算時(shí)段數(shù)和水位觀測站個(gè)數(shù)。

糙率取值應(yīng)在一個(gè)符合實(shí)際的范圍內(nèi)，即將糙率取值范圍作為約束條件，可表示為：

(4)

此模型的目標(biāo)函數(shù)是要獲得使計(jì)算河段在計(jì)算周期內(nèi)的水位誤差平方和最小的糙率值，是對(duì)河道糙率取值的綜合考慮，不排除個(gè)別時(shí)段某個(gè)斷面水位誤差過大的情況，因此應(yīng)給出最大水位誤差約束，即：

abs(Zi,t-Z′i,t)≤εi

(5)

式中：εi為第i個(gè)河段的最大水位誤差。

2 模型求解方法

恒定流計(jì)算時(shí)，河段沿程水位僅與當(dāng)前時(shí)段的進(jìn)出口邊界條件有關(guān)，與前一時(shí)刻的沿程水位無關(guān)，因此不同流量級(jí)之間并沒有時(shí)間上的聯(lián)系，所以在計(jì)算過程中，只需計(jì)算不同流量級(jí)不同狀態(tài)(即糙率)下的水位值，不需要對(duì)其進(jìn)行決策。在計(jì)算完成后，從第一階段開始直到最后一個(gè)階段，逐階段進(jìn)行決策，從而得到整個(gè)階段的最優(yōu)解。對(duì)于第i個(gè)子河段第k個(gè)階段，其遞推方程可表示為：

Fik(xik)=min[f(xik)+Fi,k-1(xi,k-1)]

(6)

式中：f(xik)是子河段i在第k階段的糙率取值為xik時(shí)的水位計(jì)算值與實(shí)測值之差的平方；Fik(xik)表示子河段i從第一階段到第k階段的最小水位計(jì)算值與實(shí)測值之間的誤差平方之和。

圖1給出了基于動(dòng)態(tài)規(guī)劃算法的河道糙率反演模型計(jì)算流程。

圖1 基于動(dòng)態(tài)規(guī)劃算法的河道糙率反演模型計(jì)算流程圖Fig.1 Flowchart of roughness inverse based on dynamic programming

3 實(shí)例應(yīng)用

將上述模型和方法應(yīng)用到向家壩庫區(qū)河道恒定流計(jì)算的糙率反演中。向家壩壩址位于金沙江下游，上距溪洛渡水電站壩址157 km，下游距離宜賓市河道里程32 km，回水長度156.6 km。計(jì)算范圍為溪洛渡水文站J99至向家壩壩前J18(距向家壩壩址約1.3 km)，全長約149.446 km，沿程共設(shè)有82個(gè)斷面，斷面間距一般為500～2 000 m。庫區(qū)設(shè)有大新號(hào)、屏山、新市、冒水、下河壩、檜溪、新灘、溪洛渡等水位(文)站，具體如圖2所示。

計(jì)算時(shí)將河道從下游到上游依次劃分為壩前～大新號(hào)段、大新號(hào)～屏山段、屏山～新市段、新市～冒水段、冒水～下河壩段、下河壩～檜溪段、檜溪～新灘段和新灘～溪洛渡段，共8個(gè)子河段，并假定每個(gè)子河段內(nèi)的斷面糙率值相同。

圖2 向家壩庫區(qū)典型斷面分布示意圖Fig.2 Distribution of typical cross sections in Xiangjiaba reservoir

采用2007年6月19日2時(shí)至2007年8月15日20時(shí)實(shí)測水位流量資料對(duì)向家壩庫區(qū)河道進(jìn)行糙率反演，入庫流量變化范圍為2 361.6～18 221.06 m3/s。按流量大小將糙率反演計(jì)算分為2 000～3 000、3 000～4 000、4 000～5 000、5 000～7 000、7 000～8 000、8 000～9 000、9 000～10 000、10 000～15 000和15 000～20 000 m3/s等9個(gè)階段，并假定每個(gè)階段的糙率離散區(qū)間為[0.01, 0.12]。

比較不同計(jì)算精度對(duì)計(jì)算結(jié)果的影響，將各階段糙率離散點(diǎn)數(shù)分別設(shè)置為12，23和56，則對(duì)應(yīng)的計(jì)算精度分別為0.01，0.005和0.002。表1給出了典型斷面在不同離散點(diǎn)數(shù)下的水位計(jì)算平均誤差。由表1可見，隨著離散點(diǎn)數(shù)的增加，典型斷面水位的平均絕對(duì)誤差都有所降低，當(dāng)離散點(diǎn)數(shù)為56時(shí)，各典型斷面水位平均絕對(duì)誤差小于15 cm。

表1 不同離散點(diǎn)數(shù)下的典型斷面水位平均絕對(duì)誤差Tab.1 The mean absolute error of water stage of typical cross sections at different discrete points

圖3給出了典型斷面在計(jì)算精度為0.002時(shí)水位計(jì)算值與實(shí)測值的差值。從圖3可見，采用上述模型和方法進(jìn)行反演糙率計(jì)算得到的水位過程線與實(shí)測資料吻合良好，各典型斷面水位計(jì)算值與實(shí)測值的差值大部分集中在-0.2～0.2 m以內(nèi)，說明本文建立的河道糙率反演模型具有良好的模擬效果，計(jì)算精度較高。

本文采用恒定流模型模擬水流運(yùn)動(dòng)，流量沿程不變，實(shí)際水流流量沿程是有變化的，如圖4給出了溪洛渡水文站和屏山水文站的流量過程。由圖4可見，兩者還是存在一定的差別，這會(huì)影響某些流量下糙率計(jì)算精度，從而增加這些時(shí)段水位計(jì)算誤差。此外，計(jì)算中沒有考慮河床沖淤變化也會(huì)使水位計(jì)算產(chǎn)生誤差。

4 結(jié) 語

本文基于動(dòng)態(tài)規(guī)劃算法，以水位誤差平方和最小為目標(biāo)，建立了河道恒定流糙率反演模型，對(duì)河道不同流量級(jí)下的糙率進(jìn)行優(yōu)化。將所建模型應(yīng)用于向家壩庫區(qū)恒定流計(jì)算的糙率反演中，給出向家壩庫區(qū)河道不同流量級(jí)下的糙率值。計(jì)算結(jié)果表明，本文建立的模型具有良好的模擬效果，且糙率區(qū)間離散點(diǎn)數(shù)越多，模擬的效果越理想，當(dāng)離散點(diǎn)數(shù)為56、計(jì)算精度為0.002時(shí)得到的典型斷面水位平均絕對(duì)誤差不超過15 cm。本研究將動(dòng)態(tài)規(guī)劃算法用于河道糙率反演計(jì)算，用智能調(diào)參代替人工調(diào)參，避免了人工調(diào)試的不確定性，對(duì)于大型河網(wǎng)糙率率定和水流模擬，該方法具有較大的實(shí)用價(jià)值和推廣意義。

圖3 典型斷面在計(jì)算精度為0.002時(shí)水位計(jì)算值與實(shí)測值的差值Fig.3 Difference of calculated and measured water stage of typical cross sections, calculation precision is 0.002

圖4屏山水文站和溪洛渡水文站實(shí)測流量過程比較Fig.4 Comparison of Pingshan hydrological station and Xiluodu hydrological station flow processes

□

[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] Wasantha Lal A M. Calibration of riverbed roughness [J]. Journal of Hydraulic Engineering，1995，121(9)：664-671.

[3] 金忠青，韓龍喜，張 健. 復(fù)雜河網(wǎng)的水力計(jì)算及參數(shù)反問題[J]. 水動(dòng)力學(xué)研究與進(jìn)展(A輯)，1998，13(3)：280-285.

[4] 程偉平，毛根海. 基于帶參數(shù)的卡爾曼濾波的河道糙率動(dòng)態(tài)反演研究[J]. 水力發(fā)電學(xué)報(bào)，2005，24(2)：123-127.

[5] 張 潮. 基于機(jī)器學(xué)習(xí)的河網(wǎng)糙率反演[D]. 杭州：浙江大學(xué)，2008.

[6] 劉志賢. 基于遺傳算法的河網(wǎng)糙率反分析研究[D]. 杭州：浙江大學(xué)，2008.

[7] 李 麗，王加虎，王建群，等. 自適應(yīng)隨機(jī)搜索算法在河網(wǎng)數(shù)學(xué)模型糙率反演中的應(yīng)用[J]. 水利水電科技進(jìn)展，2011，31(5)：64-67.

[8] 齊鄂榮，羅 昌. 庫區(qū)河道非恒定流糙率的選取及特性[J]. 武漢大學(xué)學(xué)報(bào)(工學(xué)版)，2003，36(2)：1-5.