本發(fā)明屬于河道水流及泥沙沖淤計算領(lǐng)域，尤其涉及一種基于動態(tài)規(guī)劃逐次逼近法的河道糙率反演方法。

背景技術(shù)：

糙率是水力計算中一個相當(dāng)重要又敏感的關(guān)鍵參數(shù)，對河道水流及其沖淤變化的計算成果影響很大。同時糙率的影響因素多且復(fù)雜，不僅與周界壁面粗糙程度有關(guān)，還受河流形態(tài)及水流特性的影響，目前比較常用的推求糙率的方法主要分為傳統(tǒng)方法和智能反演算法兩大類。

傳統(tǒng)方法主要包括經(jīng)驗法、半經(jīng)驗法以及試錯法。經(jīng)驗法，即根據(jù)河床材料及其粗糙程度直接查表得到。此法過于粗糙，糙率表本身也有局限性。同時該方法也無法考慮河槽、灘面等因素的影響。半經(jīng)驗法，利用實測水文資料及地形資料先求出河段的某些平均量，然后代入謝才公式并求出謝才系數(shù)，再利用巴甫洛夫斯基公式或曼寧公式估算出該河段的平均糙率值；該方法一般只適用于地形較規(guī)則的均勻流態(tài)，對沿程變化較大的天然河流來說不太適用。試錯法，首先利用經(jīng)驗公式給定初步糙率，利用數(shù)值模擬技術(shù)進(jìn)行模擬，通過比較計算水位與實際水位之間的差距，反過來調(diào)整初始糙率，重復(fù)這一過程直至計算水位與實際水位之間的差距滿足精度要求；此方法具有較大的經(jīng)驗性和任意性，尤其在河網(wǎng)計算中，手工調(diào)試費時費力、工作量大，精度難以保證，且當(dāng)計算條件改變時，糙率也會隨之變化，需要經(jīng)常更新。

反演算法將沿程糙率看作優(yōu)化向量，將水位、流量計算值與實測值之間的誤差作為目標(biāo)函數(shù)，通過優(yōu)化算法，尋求使得目標(biāo)函數(shù)最小化的最優(yōu)向量。20世紀(jì)60年代以來，國內(nèi)外學(xué)者提出一些反演糙率的算法，一般是將各河道糙率視為定值，都沒有考慮糙率隨流量的變化，實際上各河道不同流量下的糙率是不同的。

技術(shù)實現(xiàn)要素：

為了解決上述問題，本發(fā)明提出了一種基于動態(tài)規(guī)劃逐次逼近法的河道糙率反演方法，包括

步驟1、采用一維水動力學(xué)模型模擬明渠、河道水流變化，模型中的糙率n采用優(yōu)化方法進(jìn)行率定，以計算河段在計算周期內(nèi)的水位誤差平方和最小作為目標(biāo)，建立河道糙率反演模型，對不同流量下的糙率進(jìn)行反演；

步驟2、采用動態(tài)規(guī)劃逐次逼近法進(jìn)行河道糙率反演，將不同流量級下的糙率反演分解成多階段優(yōu)化問題，每次計算只對當(dāng)前階段的糙率進(jìn)行尋優(yōu)，其它階段的糙率值保持不變。

所述步驟1中一維水動力學(xué)模型為

采用saint-venant方程組描述河道一維水流運動，其水流連續(xù)方程和運動方程分別為

式中：t為時間；x為流程；q為流量；z為水位；a為過水?dāng)嗝婷娣e；b為河寬；r為水力半徑；n為糙率；v為斷面平均流速；ql和ul為河段單位長度側(cè)向入流量和側(cè)向入流在x方向的分量；

動量修正系數(shù)α1＝(∫au²da)/(q²/a)；g為重力加速度；

采用線性化的preissmann四點隱格式對方程(1)和(2)進(jìn)行離散，并采用追趕法求解離散后的方程組，方程(2)包含待率定的參數(shù)-糙率n，以水位誤差平方和最小為目標(biāo)，建立河道糙率反演模型。

所述建立河道糙率反演模型具體為

根據(jù)水位觀測站的分布，將河道劃分若干個子河段，根據(jù)實測流量變化范圍，將流量劃分為若干等級，各子河段在不同流量級下分別采用不同的代表糙率，即各子河段在第j流量級下的糙率表示為xmj表示第m個河段在第j流量級下的糙率值，m＝1,2,…,m，其中m代表子河段總數(shù)；選擇水位過程作為目標(biāo)函數(shù)的變量，根據(jù)最小二乘判定標(biāo)準(zhǔn)構(gòu)造如下目標(biāo)函數(shù)：

式中：為第j個流量級下各斷面的糙率,j＝1,2,…,n；zi,t為t時刻第i個水位站的計算水位值；為t時刻第i個水位站的實測水位值；t、m和n分別為計算時段數(shù)、水位觀測站個數(shù)和流量級個數(shù)；

約束條件：

糙率取值范圍約束

式中：分別為糙率取值的下限和上限；

最大水位誤差約束

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

abs(zi,t-z′i,t)≤εi(5)

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

所述步驟2中采用動態(tài)規(guī)劃逐次逼近法進(jìn)行河道糙率反演的基本步驟為

201、根據(jù)流量級將糙率反演劃分為n個階段，并將每個階段各斷面糙率在內(nèi)進(jìn)行離散；

202、給出每個階段各子河段滿足約束條件的初始糙率值j＝1,2,…,n；

203、固定第2個至第n個階段下各子河段糙率值，利用動態(tài)規(guī)劃法確定第1個階段各子河段糙率的最優(yōu)值表示第1個流量階段下第m個河段第一次迭代計算得到的糙率值；m＝1,2,…,m，其中m代表子河段總數(shù)；在對第1個階段糙率優(yōu)化計算的過程中，目標(biāo)函數(shù)值要包含第2至第n個階段初始糙率下的計算結(jié)果；

204、將第1階段優(yōu)化后的糙率值以及第3至第n個階段下的糙率值固定，對第2階段進(jìn)行同樣的優(yōu)化計算，得到第2個階段的最優(yōu)糙率值表示第2個流量階段下第m個河段第一次迭代計算得到的糙率值；

205、按上述方式依次對階段j＝4,…,n進(jìn)行優(yōu)化計算，其余階段的糙率固定不變，直到計算完所有階段；

206、重復(fù)第203～205步驟，進(jìn)行第2次、第3次……迭代計算，直到前后兩次優(yōu)化后的目標(biāo)函數(shù)滿足精度要求或達(dá)到迭代次數(shù)為止，計算結(jié)束，則最后一次的優(yōu)化計算結(jié)果即為所求的最優(yōu)解。

有益效果

本發(fā)明的目的是在糙率反演中通過采用計算機智能調(diào)參，給出不同流量級下的糙率變化，各流量下對應(yīng)的糙率可通過插值得出。本發(fā)明不需考慮后效性的影響，考慮糙率隨流量的變化情況得到不同流量下的糙率值。本發(fā)明可避免其他方法存在的受初值影響較大、容易早熟或陷入局部最優(yōu)的問題，具有降低維數(shù)、減少計算量、計算收斂速度較快等特點。將計算機智能調(diào)參代替人工調(diào)參，節(jié)省人力物力，避免手工調(diào)試的不確定性，為提高水力仿真精度及實現(xiàn)糙率精確率定奠定了良好的基礎(chǔ)。

附圖說明

圖1采用動態(tài)規(guī)劃逐次逼近法進(jìn)行河道糙率反演流程圖；

圖2檜溪水位觀測站各迭代步水位計算值與實測值比較圖；

圖3各典型斷面第4次迭代計算水位誤差圖；其中圖a為大新號斷面，圖b為屏山斷面，圖c為新市斷面，圖d為冒水?dāng)嗝?，圖e為下河壩斷面，圖f為檜溪斷面，圖g為新灘斷面，圖h為溪洛渡斷面。

具體實施方式

下面結(jié)合附圖，對本發(fā)明作詳細(xì)說明。本發(fā)明提出了一種基于動態(tài)規(guī)劃逐次逼近法的河道糙率反演方法，包括

步驟1、采用一維水動力學(xué)模型模擬明渠、河道水流變化，模型中的糙率n采用優(yōu)化方法進(jìn)行率定，以計算河段在計算周期內(nèi)的水位誤差平方和最小作為目標(biāo)，建立河道糙率反演模型，對不同流量下的糙率進(jìn)行反演；

步驟2、采用動態(tài)規(guī)劃逐次逼近法進(jìn)行河道糙率反演，將不同流量級下的糙率反演分解成多階段優(yōu)化問題，每次計算只對當(dāng)前階段的糙率進(jìn)行尋優(yōu)，其它階段的糙率值保持不變。

所述步驟1中一維水動力學(xué)模型為

采用saint-venant方程組描述河道一維水流運動，其水流連續(xù)方程和運動方程分別為

式中：t為時間；x為流程；q為流量；z為水位；a為過水?dāng)嗝婷娣e；b為河寬；r為水力半徑；n為糙率；v為斷面平均流速；ql和ul為河段單位長度側(cè)向入流量和側(cè)向入流在x方向的分量；

動量修正系數(shù)α1＝(∫au²da)/(q²/a)；g為重力加速度；

采用線性化的preissmann四點隱格式對方程(1)和(2)進(jìn)行離散，并采用追趕法求解離散后的方程組，方程(2)包含待率定的參數(shù)-糙率n，以水位誤差平方和最小為目標(biāo)，建立河道糙率反演模型。

所述建立河道糙率反演模型具體為

根據(jù)水位觀測站的分布，將河道劃分若干個子河段，根據(jù)實測流量變化范圍，將流量劃分為若干等級，各子河段在不同流量級下分別采用不同的代表糙率，即各子河段在第j流量級下的糙率表示為xmj表示第m個河段在第j流量級下的糙率值，m＝1,2,…,m，其中m代表子河段總數(shù)；選擇水位過程作為目標(biāo)函數(shù)的變量，根據(jù)最小二乘判定標(biāo)準(zhǔn)構(gòu)造如下目標(biāo)函數(shù)：

式中：為第j個流量級下各斷面的糙率,j＝1,2,…,n；zi,t為t時刻第i個水位站的計算水位值；z′i,t為t時刻第i個水位站的實測水位值；t、m和n分別為計算時段數(shù)、水位觀測站個數(shù)和流量級個數(shù)；

約束條件：

糙率取值范圍約束

式中：分別為糙率取值的下限和上限；

最大水位誤差約束

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

abs(zi,t-z′i，t)≤εi(5)

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

所述步驟2中采用動態(tài)規(guī)劃逐次逼近法進(jìn)行河道糙率反演的基本步驟為

201、根據(jù)流量級將糙率反演劃分為n個階段，并將每個階段各斷面糙率在內(nèi)進(jìn)行離散；

202、給出每個階段各子河段滿足約束條件的初始糙率值j＝1,2,…,n；

203、固定第2個至第n個階段下各子河段糙率值，利用動態(tài)規(guī)劃法確定第1個階段各子河段糙率的最優(yōu)值表示第1個流量階段下第m個河段第一次迭代計算得到的糙率值；m＝1,2,…,m，其中m代表子河段總數(shù)；在對第1個階段糙率優(yōu)化計算的過程中，目標(biāo)函數(shù)值要包含第2至第n個階段初始糙率下的計算結(jié)果；

204、將第1階段優(yōu)化后的糙率值以及第3至第n個階段下的糙率值固定，對第2階段進(jìn)行同樣的優(yōu)化計算，得到第2個階段的最優(yōu)糙率值表示第2個流量階段下第m個河段第一次迭代計算得到的糙率值；

205、按上述方式依次對階段j＝4,…,n進(jìn)行優(yōu)化計算，其余階段的糙率固定不變，直到計算完所有階段；

206、重復(fù)第203～205步驟，進(jìn)行第2次、第3次……迭代計算，直到前后兩次優(yōu)化后的目標(biāo)函數(shù)滿足精度要求或達(dá)到迭代次數(shù)為止，計算結(jié)束，則最后一次的優(yōu)化計算結(jié)果即為所求的最優(yōu)解。

將上述模型和方法應(yīng)用到向家壩庫區(qū)河道非恒定流計算的糙率反演中。計算范圍為溪洛渡水文站至向家壩壩前距向家壩壩址約1.3km，全長約149.446km，沿程共設(shè)有82個斷面，斷面間距一般為500～2000m。

經(jīng)過4次迭代后，河段糙率基本達(dá)到穩(wěn)定。表1給出了典型斷面各流量級下不同迭代次數(shù)的糙率系數(shù)計算結(jié)果。由表1可見，本發(fā)明提出的計算方法的收斂性較好，經(jīng)過4次迭代后，各典型斷面在不同流量級下的糙率基本可達(dá)到一穩(wěn)定值。

表1典型斷面不同流量級下不同迭代次數(shù)的糙率系數(shù)

表2給出了典型斷面不同迭代次數(shù)下的水位計算平均絕對誤差。由表2可見，隨著迭代次數(shù)的增加，水位計算值與實測水位誤差也逐漸減小，經(jīng)4次迭代后各典型斷面的最大平均誤差不超過15cm。

表2典型斷面不同迭代次數(shù)下水位平均絕對誤差

圖2為檜溪水位觀測站不同迭代次數(shù)時水位計算值與實測值變化過程的比較結(jié)果。由圖2可見，隨著糙率的迭代收斂，檜溪斷面的水位計算值逐漸趨向?qū)崪y值：按糙率初始值0.03計算得到的水位值與實測值相差很大；采用本發(fā)明建立的糙率反演模型經(jīng)過1次迭代計算后得到的水位值已經(jīng)開始接近實測水位值，明顯優(yōu)于根據(jù)糙率初始值的計算結(jié)果；隨著迭代次數(shù)的增加，水位計算值逐漸逼近實測值，經(jīng)過4次迭代后得到的水位計算值已非常接近實測值。

圖3給出了第4次迭代計算后，各典型斷面計算得到的水位與實測水位的誤差值。由圖3可見，各典型斷面的水位誤差值大部分都在-0.2～0.2m之間，說明本發(fā)明建立的河道糙率反演模型具有很好的反演效果。