基于奇異矩陣分解法的河道糙率反演計算方法

王玲玲,鐘 娜,成高峰

(1.河海大學(xué)水利水電學(xué)院,江蘇南京 210098;2.華東勘測設(shè)計研究院,浙江杭州 310014)

河道糙率不僅與表面粗糙程度有關(guān),還受水力要素及水流特性的影響,是一個綜合水力摩阻系數(shù),對河道水流及其沖淤變化的計算成果有很大影響.文獻[1-3]提出了幾種河道糙率的計算公式和方法,但這些計算公式和方法[1-3]大多屬于半經(jīng)驗性的,所得結(jié)果往往與實際偏差較大.工程實踐中常根據(jù)河道的實測水文資料通過試錯法推求糙率.若河網(wǎng)規(guī)模較大,用試錯法調(diào)試每一河段的糙率,工作量將十分巨大.文獻[4-5]提出的系統(tǒng)參數(shù)反演的脈沖譜、離散-優(yōu)化等方法,雖具有理論基礎(chǔ),但不能應(yīng)用于大規(guī)模欠定的河網(wǎng)糙率反分析.本文借助文獻[6]的基本思路,采用建立在對輸出誤差不斷進行迭代校正基礎(chǔ)上的反分析方法和奇異矩陣分解(singular value decomposition,SVD)法[7-8]建立了河道糙率反演計算方法,并通過編程實現(xiàn)了東江108 km河道的糙率反演,得到了較為精確的糙率分布.

1 反演計算方法的建立

1.1 理論基礎(chǔ)

假設(shè)待率定參數(shù)與有效空間內(nèi)的變量存在連續(xù)函數(shù)關(guān)系,給待率定參數(shù)一個微小的變化量,建立影響系數(shù)矩陣,若計算域中各物理量的輸出值與實測值之間誤差不斷減小,則可以求得待率定參數(shù)的空間分布.

設(shè)ykj表示k(k=1,2,…,l)時刻、j(j=1,2,…,m)位置通過數(shù)學(xué)模型計算的物理量(如測站水位),Ykj是同一變量的實際觀測值.目標(biāo)則是找出參數(shù)xi(糙率),i=1,2,…,n,使誤差εkj=ykj-Ykj在ykj的非恒定變化過程中逐漸減小.其中:n為率定糙率參數(shù)總數(shù);m為水文測站總數(shù);l為非恒定計算的時刻點總數(shù).

一維河道沿程水位ykj與待率定糙率xi之間的關(guān)系滿足圣維南方程組,其具體數(shù)值求解方法詳見文獻[9].

用矩陣表示,則式(1)可寫成

式中:A k,m×n——影響系數(shù)矩陣;k——計算的時刻點.

Ak,m×n的確定過程如下:給每個河道斷面糙率參數(shù)xi一個初始值并將其代入河道計算模型,計算出初始糙率條件下的各個測站的非恒定水位過程ykj(xi);按照1到n的順序給每一個待率定參數(shù)xi一個很小的增量(如Δxi=0.001)就得出新的參數(shù)組合,將其代入河道計算模型,計算出測站非恒定過程的水位值ykj(xi+Δx);將每次計算得出的ykj(xi+Δx)和ykj(xi)代入式(1)進行計算,就可以得到參數(shù)xi對j測站物理量的影響系數(shù)akji.當(dāng)有n個參數(shù)需要率定時,必須進行n+1次河道模型的計算.

設(shè)有m個測站,根據(jù)各個測站每一時刻計算值與實測值的誤差,可得出各個測站在整個計算時段內(nèi)的平均誤差 εsj:

寫成矩陣形式,得

其中

求解方程(5),即可得到滿足誤差 εs=O的率定參數(shù)xi的修正值矩陣Δx.對糙率初始分布進行修正,可獲得下一輪計算的新的糙率分布:

式中:r——迭代次數(shù);ρr——控制增量幅度的欠松弛參數(shù).修正過程中,對 xi采用約束條件

予以控制.其中 xub,xlb為指定的上、下邊界.

反復(fù)進行上述迭代與修正過程,直到各測站水位計算值與實測值誤差的2范數(shù)‖εs‖2在給定的誤差允許范圍內(nèi),且糙率參數(shù)修正值收斂為止.

1.2 SVD方法的應(yīng)用

上述數(shù)理問題求解結(jié)果的精度,很大程度上依賴于所形成大型代數(shù)方程組的求解方法.對于測站數(shù)m遠小于待率定斷面數(shù)n的情況,所形成的代數(shù)方程組將嚴重病態(tài),其求解成為反演計算的關(guān)鍵.SVD方法可直接用于求解m=n,m＞n,m＜n,特別是系數(shù)矩陣為奇異矩陣的問題[11].方程(5)的求解可通過將系數(shù)矩陣 A s分解成 A s=U∑VT的方法來實現(xiàn).將式(5)寫成

圖1 東江河道走勢及部分水位測點分布Fig.1 Trend of Dongjiang River and distribution of measuring points

其中 Δz=VTΔx,d=U-1εs.由此可以得到

式中λi為∑對角線上的各元素.如果至少存在一個 λi的值為0或接近于0,則矩陣為奇異矩陣.此時,對應(yīng)于該 λ值的1/λ均賦值為0.通過這種方式可以得到該奇異矩陣的唯一解.另外,控制式(8)中最小的奇異值 λmin的舍入誤差對于減小結(jié)果誤差也很有效.因為 λmin很小,則1/λ很大,結(jié)果的誤差會被無限放大.對于非奇異矩陣,用λmax/λmin代替奇異矩陣中的控制條件 λmin,可防止該矩陣成為病態(tài)矩陣[12].

2 反演計算方法的應(yīng)用

將上述反演計算方法應(yīng)用于廣東東江108km長河段的糙率計算.

2.1 水文資料及計算工況

東江是珠江水系三大河流之一,干流全長520 km,集水面積27040 km2.本次計算從三方塘水文站至崗頭村水文站,全長108km,沿途設(shè)有18個水位測點,河道走勢及部分測站位置如圖1所示.

沿河劃分64個計算斷面,即63個河段,空間步長Δx=108.331m,時間步長Δt=20min;河床糙率初值為0.03,欠松弛參數(shù)為0.8,上、下邊界取 x ub=0.1,x lb=0.001.根據(jù)2005年6月2日15:00—4日12:00約50h時間域內(nèi)16個測點的實測水位資料反演63個河段的糙率值.本算例是個欠定問題,所形成的代數(shù)方程組將嚴重病態(tài).

采用上游三方塘水位站和下游崗頭村水位站的實測水位變化過程作為上、下游邊界條件(圖2).以50h恒定流的計算結(jié)果作為非恒定流計算的初始條件.

2.2 糙率反演結(jié)果及分析

2.2.1 糙率系數(shù)計算收斂過程

采用上述糙率反演計算方法,經(jīng)過8次迭代計算,得到64個斷面糙率系數(shù)的穩(wěn)定結(jié)果.圖3和圖4所示為部分斷面糙率系數(shù)計算收斂過程和相應(yīng)各水位測站計算值與實測值誤差的降低過程.

圖2 上、下游邊界水位過程線Fig.2 Hydrographs of upstream and downstream water levels

圖3 沿程若干斷面糙率系數(shù)的收斂過程Fig.3 Convergence process of roughness coefficients of various sections

圖4 水位測點計算值與實測值誤差的降低過程Fig.4 Reduction process of errors between calculated and measured values

由圖3和圖4可見,本文計算方法在得到穩(wěn)定的糙率系數(shù)時,計算水位與實測水位誤差也穩(wěn)定減小.

2.2.2 糙率反演結(jié)果

圖5所示為8次迭代計算后的糙率沿程分布.

圖5表明,東江河床的糙率范圍在0.02～0.045之間,在平原河流糙率的經(jīng)驗值范圍內(nèi).結(jié)合地形資料分析可知:糙率在0.02～0.03之間的河段主要表現(xiàn)為單一河槽的河段,且河槽中植被較少;糙率大于0.03的河段則主要表現(xiàn)為灘槽共存的復(fù)式河段,灘地有較為稠密的雜草,阻水能力較強.因此,由糙率反演結(jié)果也可以大致判斷河段的斷面形態(tài)以及植被的稠密程度.

圖5 8次迭代計算后的糙率沿程分布Fig.5 Distribution of roughness coefficients after 8 iterations

圖6 糙率率定過程中纜西鎮(zhèn)水位非恒定過程Fig.6 Unsteady p rocess of water levels in roughness calibration at Lanxi town

2.2.3 非恒定水位逼近過程

雖然沿程各水位測點的實測水位與計算水位的差值在整個非恒定過程中是穩(wěn)定收斂的(圖4),但不足以說明非恒定過程各個時刻計算值與實測值吻合程度的提高,因此需要對糙率率定過程各測站水位計算值與實測值的非恒定變化過程進行比較.圖6為糙率率定過程纜西鎮(zhèn)測站各迭代步水位計算值與實測值的比較結(jié)果.

由圖6可以看出,隨著糙率的迭代收斂,測站水位的計算值趨向于實測值的過程:由糙率初值(0.03)計算的水位非恒定過程,與實測值相差很大;采用糙率率定模型對初值進行修正后獲得的糙率新值,計算的水位過程與實測值的吻合程度,明顯好于修正前;隨著率定過程的進行,水位結(jié)果逐漸趨向?qū)崪y值;當(dāng)糙率系數(shù)經(jīng)過8次迭代時,該糙率計算的非恒定過程與實測資料基本吻合.可見,糙率的尋優(yōu)過程,就是各個測站各個時刻的水位計算值逐漸趨向?qū)崪y值的過程.同時說明,在糙率率定時,采用水位計算值與實測值的誤差在整個非恒定過程中逐漸減小作為目標(biāo)函數(shù)是可行且易于實現(xiàn)的.

3 反演計算方法的驗證

為了驗證糙率反演計算方法的精度,采用上述糙率率定結(jié)果,對2005年6月10日17:00—12日18:00該河段各測點的水位和流量進行了計算,并與M IKE11軟件計算結(jié)果以及實測結(jié)果進行了對比.圖7所示為檳榔潭整個計算時段水位本文方法計算值與M IKE11軟件計算值以及實測值的比較,圖8所示為6月11日16:00該河段流量本文方法計算值與M IKE11軟件計算值的比較以及本文方法計算值與觀音閣1號、2號和3號3個測站實測值的比較.

圖7 檳榔潭水位本文方法計算值與M IKE11軟件計算值以及實測值的比較Fig.7 Comparison among values calculated by present method and software M IKE11 and measured ones of water levels at Binglangtan

圖8 流量本文方法計算值與MIKE11軟件計算值以及觀音閣實測值的比較Fig.8 Comparison among values calculated by present method and software M IKE11 and measured ones of discharges at Guanyinge

圖7表明,水位本文方法計算值與實測值的誤差較小,絕對誤差在10 cm以內(nèi),而與MIKE11軟件計算值的誤差較大.圖8表明,2種方法的流量計算結(jié)果均與實測結(jié)果接近,而本文方法的流量計算結(jié)果更接近于實測結(jié)果.由此可以看出,本文糙率反演方法得到的結(jié)果滿足精度要求,具有一定的可靠性.

4 結(jié) 語

為了獲得精確的糙率系數(shù),采用奇異矩陣分解法對河道糙率進行率定,形成了一套包括模型率定、水動力場預(yù)測的計算模型,并通過編制程序來實現(xiàn)算法.東江108km天然河道糙率率定結(jié)果表明,該糙率率定方法精度較高,可用于河道過流能力及復(fù)雜河網(wǎng)水動力模擬分析,可解決欠定的糙率率定問題并提高模型的實用性和洪水預(yù)報精度.

[1]何建京,王惠民.流動形態(tài)對曼寧糙率系數(shù)的影響研究[J].水利水電技術(shù),2002,22(6):22-24.(HE Jian-jing,WANG Huim in.The effects of types of flow on manning's roughness coefficient[J].Hydrology,2002,22(6):22-24.(in Chinese))

[2]程進豪,安連華,王華,等.黃河山東段河床糙率分析[J].水利學(xué)報,1997,28(1):39-43.(CHENG Jin-hao,AN Lian-hua,WANG Hua,et al.An analysis of the bed roughness in Shandong reaches of the Yellow River[J].Journal of Hyd raulic Engineering,1997,28(1):39-43.(in Chinese))

[3]李榕.關(guān)于影響曼寧粗糙系數(shù)n值的水力因素探討[J].水利學(xué)報,1989,20(12):62-66.(LI Rong.The study on hydraulic factors effect on manning's roughness coefficient n[J].Journal of Hydraulic Engineering,1989,20(12):62-66.(in Chinese))

[4]WANG Ling-ling.The algorithms and applications for determining permeability coefficients by pu lse spectrum technique[J].Journal of Hydrodynamics:Ser B,2002,14(2):112-115.

[5]戴會超,王玲玲.參數(shù)控制反問題的求解及其在壩工滲流中的應(yīng)用[J].水力發(fā)電學(xué)報,2005,24(6):72-77.(DAI Hui-chao,WANG Ling-ling.The algorithms of parameter-control inverse problems and their applications in seepage of dam engineering[J].Journal of Hydroelectric Engineering,2005,24(6):72-77.(in Chinese))

[6]WASANTHA LAL A M.Calibration of riverbed roughness[J].Journal of Hydraulic Engineering,1995,121(9):664-671.

[7]楊克君,曹叔尤,劉興年.復(fù)式河槽綜合糙率計算方法比較與分析[J].水利學(xué)報,2005,36(7):780-786.(YANG Ke-jun,CAO Shu-you,LIU Xing-nian.Analysis on methods for predicting composite roughness of river channel with compound cross section[J].Journal of Hydraulic Engineering,2005,36(7):780-786.(in Chinese))

[8]潘光斌,陳光裕.基于SVD的測量系統(tǒng)建模方法研究[J].電子測量與儀器學(xué)報,2006,20(5):60-62.(PAN Guang-bin,CHEN Guang-yu.Design of indirectly measuring model based on SVD[J].Journal of Electronic Measurement and Instrument,2006,20(5):60-62.(in Chinese))

[9]李光熾,王船海.實用河網(wǎng)水流計算[M].南京:河海大學(xué)出版社,2002.

[10]BECKER L,YEH W W-G.Identification of parameters in unsteady open channel flows[J].Water Resour Res,1972,8(4):956-965.

[11]MENEKE W.Geophysical data analysis:discrete inverse theory[M].New York:Academic Press,Inc.,1984.

[12]鐘娜.復(fù)式河道一維數(shù)值模擬[D].南京:河海大學(xué),2007.