動態(tài)規(guī)劃算法在河道糙率反演中的應(yīng)用

于顯亮，彭 楊，吳志毅

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

0 引 言

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

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

1 模型建立

1.1 一維恒定流模型

對于長河段長時段的非恒定水流運動，可將其概化為平均的洪水過程，作為恒定問題處理，不考慮流量沿程變化，此時水流運動可直接用運動方程來描述，即

(1)

式中：x為流程；Q為流量；A為過水斷面面積；Z為水位；R為水力半徑；n為糙率；g為重力加速度。

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

(2)

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

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

1.2 河道糙率反演模型

(3)

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

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

(4)

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

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

(5)

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

2 模型求解方法

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

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

(6)

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

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

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

3 實例應(yīng)用

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

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

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

采用2007年6月19日2時至2007年8月15日20時實測水位流量資料對向家壩庫區(qū)河道進行糙率反演，入庫流量變化范圍為2 361.6～18 221.06 m3/s。按流量大小將糙率反演計算分為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個階段，并假定每個階段的糙率離散區(qū)間為[0.01, 0.12]。

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

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

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

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

4 結(jié) 語

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

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

圖4屏山水文站和溪洛渡水文站實測流量過程比較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)的水力計算及參數(shù)反問題[J]. 水動力學(xué)研究與進展(A輯)，1998，13(3)：280-285.

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

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

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

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

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