新安江模型參數線性化率定研究

2022-06-21 09:53:08王紅艷王新龍

水力發電 2022年1期

王紅艷，王新龍，周倩，楊宇

(1.寧波弘泰水利信息科技有限公司，浙江寧波 315000；2.寧波市水利水電規劃設計研究院有限公司，浙江寧波 315000；3.江蘇省水文水資源勘測局徐州分局，江蘇徐州 221000)

0 引言

流域水文模型在進行水文規律研究和解決實際生產問題中起著重要作用。中國現在最常用的流域水文模型是新安江模型[1]，模型結構考慮了濕潤半濕潤地區的蒸散發、產匯流特性。采用新安江模型建立實際流域的洪水預報方案時，首先要根據流域的實際水文資料對新安江模型參數進行率定，然后根據率定的參數制定洪水預報方案，對未來的洪水進行預報[2]。

流域水文模型除了模型結構要合理外，模型參數的率定也是一個十分重要的環節[3]，模型參數的好壞直接關系著洪水預報方案的優劣。傳統的參數率定方法為人機交互率定，其主要根據人們的先驗知識進行判斷、估計。這種方法雖然簡單易行，但調試參數需要大量時間且得到的參數不一定是最優的[4]。目前，國內外研究較多的是模型參數自動率定方法，參數的自動率定就是人們預先編制好一個參數自動尋優的計算程序，然后輸入水文資料和模型參數初始值，通過自動尋優計算得到最優的模型參數值[5-7]。

傳統上的參數優化的方法是數學尋優法，幾乎都是以誤差平方和最小為目標函數，在目標函數響應曲面上尋找參數最優值，通過一階函數求導為零得到參數的最優值。這種方法對于線性參數具有較好的效果，但是用于優選非線性參數時，會給非線性參數增加不相關的局部優值。為了避免在非線性參數優選時出現這些問題，包為民教授提出了非線性函數參數的線性化率定方法[8]，在函數響應曲面上尋找參數的最優值。本文把參數線性化方法應用于新安江模型，并通過理想流域和實際流域驗證該方法的可行性和高效性。

1 新安江模型參數線性化率定方法

1.1 新安江模型介紹

趙人俊教授于1973年在編制新安江入庫洪水預報方案時提出了新安江模型。新安江模型結構是分散性的，主要分為蒸散發計算、產流計算、分水源計算、和匯流計算四個層次。模型各層次結構的功能，計算采用的方法和參數見表1。

表1 新安江模型各層次結構功能、計算方法和相應參數

1.2 新安江模型參數線性率定原理

參數線性化率定方法對非線性參數函數以參數作為自變量求導，再通過導函數差分線性化，并對線性化的參數用誤差平方和目標函數進行率定，然后逐步逼近非線性參數的最優值。該方法收斂、率定速度快，解決了以誤差平方和為目標函數求解非線性函數參數增加不相關的局部優值的問題。與傳統的方法相比，該方法的實用性更強、效果更好。

新安江模型結構完整，計算清晰，是一個典型的非線性模型。把新安江模型看作一個非線性函數，模型的計算流量當作該函數的因變量，模型的參數當成該函數的自變量。首先，對新安江模型參數變量求導，通過導函數差分使參數線性化；然后，把新安江模型這個非線性函數進行一階泰勒級數展開，構建計算流量和模型參數之間的線性函數，通過最小二乘法求解參數。把非線性參數函數轉變為線性參數函數求解可以避免產生一些不相關的局部參數解。以流量誤差平方和、參數迭代步長收斂容差為目標函數進行率定，逐步逼近新安江模型參數全局最優值。

1.3 新安江模型參數線性率定步驟

把新安江模型這個非線性函數進行多元一階泰勒級數展開，構建計算流量q(θ，x)和模型參數θ之間的線性函數

(1)

式中，q(θj，x)=[q(θj，x1)，q(θj，x2)，…，q(θj，xL)]T為L組函數計算值，即根據新安江模型參數θj和水文資料系列x計算的防洪斷面流量過程q(θj+1，x)=[q(θj+1，x1)，q(θj+1,x2)，…，q(θj+1，xL)]T，L組真參數相應的函數值，如果模型沒有任何誤差，其值為實測的斷面流量系列，后面計算時用實測流量系列代入；θj=[(θ1，j，θ2，j，…，θn，j)]T為函數的參數變量，即新安江模型的參數變量；x=(x1，x2，…，xL)T為函數的輸入變量，即實測的雨量、蒸發等水文資料；e=(e1，e2，…，eL)T為函數的誤差，即模型的誤差。

具體率定步驟：

(1)根據參數的物理意義和流域的實際情況給定一組新安江模型參數初始值θ0。

(2)根據給定的參數向量、降雨資料、蒸發資料用新安江模型計算防洪斷面流量q(θ，x)和參數靈敏度矩陣[5]S，S反映了參數改變引起計算流量的改變程度。把S代入函數式(1)得函數式(2)，計算S時需要計算函數對各參數的偏導數，通過導函數差分使參數線性化，計算公式見(式3)、(式4)。如下

q(θj+1，x)=q(θj，x)+S(θj+1-θj)+e

(2)

(3)

Δθi，j=0.001θi，j

(4)

(5)

(3)用最小二乘法[6]求解式(2)，得到新的參數向量Qj+1和參數尋找方向Δθ。即

θj+1=θj+(STS)-1ST(q(θj+1，x)-q(θj，x))

(6)

Δθ=Qj+1-Qj

(7)

(4)以流量誤差平方和最小為目標函數確定參數尋找方向上最優步長比例系數b(1≥b>0)，得出尋找方向上最優的參數向量θj+1。即

(8)

θj+1=θj+b(STS)-1ST(q(θj+1，x)-q(θj，b，x))

(9)

(5)判斷新的參數θj+1與上一步參數θj之差是否小于給定的迭代步長收斂容差ε，即|θj+1-θj|<ε，如果小于則尋優結束，得到最優參數θj+1；否則，令θj=θj+1轉步驟(2)繼續循環計算，直到得到滿足條件的最優參數。

要證明步驟(3)中尋找方向的正確性，只要證明對于任意一步尋找的參數向量θj+1對應的流量誤差平方和Fj+1比上一步參數向量θj對應的流量誤差平方和Fj小即可。證明如下[7]

Fj+1=(q-qj+1)T(q-qj+1)

(10)

Fj+1=
[q-qj-bS(θj+1-θj)]T[q-qj-bS(θj+1-θj)]=
Fj-b(q-qj)TS(θj+1-θj)-b(θj+1-θj)TST(q-qj)+
b2(θj+1-θjTSTS(θj+1-θj)

(11)

把式(6)代入式(11)得

Fj+1=Fj-2(q-qj)TS(STS)-1ST(q-qj)+
(q-qj)TS(STS)-1ST(q-qj)=
Fj-(q-qj)TS(STS)-1ST(q-qj)

(12)

由此可見，參數率定過程中的流量誤差平方和是遞減的。隨著循環計算次數的增加，其對應的流量誤差平方和越來越小，最終趨向于最小值，參數也逐漸逼近最優的參數值。

1.4 新安江模型參數線性率定指標

在模型參數的優選中，目標函數的選取直接影響參數率定的效率和結果。本文根據對模型成果的要求，選擇的子目標函數主要有5個。

(1)計算流量和實測流量的誤差平方和

(13)

式中，Qx為計算的流量值，m3/s；Q0為實測流量值，m3/s。

(2)迭代步長收斂容差ε，反映參數率定過程中同一個參數相鄰兩個計算值之間的距離，用來判別參數率定過程是否結束的條件，若同一個參數相鄰的兩步計算值之間的距離小于ε，則認為率定的參數值已收斂到最優值附近。

(3)徑流深的相對誤差

(14)

式中，Rc(i)為計算徑流深，mm；Ro(i)為實測徑流深，mm。

(4)洪峰的相對誤差

(15)

式中，Qmax，c為計算洪峰，m3/s；Qmax，o為實測洪峰，m3/s。

(5)確定性系數DC，主要用來評價洪水預報過程和實測流量過程之間的擬合程度，確定性系數越大，表明計算流量和實測流量過程擬合得越好。即

(16)

2 理想流域驗證

選取福建東張流域作為理想流域，用給定的一組新安江模型參數值，即參數真值(表2)、實際蒸發資料、實際雨量資料計算流域出口斷面的流量過程，把該流量過程作為理想實測流量。然后根據這些水文資料進行新安江模型參數線性化率定，如果率定的參數結果能夠回到參數真值，則說明新安江模型參數線性化率定方法合理有效，可以把該方法用于實際流域，進一步研究該方法的實際應用效果。

表2 參數初值取值范圍

針對新安江模型SM、KI、KG、CS、CI、CG、KE、XE統一進行新安江模型線性率定，檢驗該方法在高維數的參數空間的應用效果。具體的率定步驟如下：

(1)在參數的物理意義范圍內隨機給SM、KI、KG、CS、CI、CG、KE、XE賦初值，各參數取值范圍如表2。

(2)根據不同的參數初值分別進行新安江模型參數線性化率定，SM、KI、KG、CS、CI、CG、KE、XE的迭代步長容差均取0.000 1，即當率定過程中SM、KI、KG、CS、CI、CG、KE、XE同時滿足迭代容差要求時率定結束。

(3)計算結果見表3。從表3可以發現，SM、KI、KG、CS、CI、CG、KE、XE取不同初始值，經過新安江模型參數線性率定方法均可以回到真值附近。參數率定過程中平均率定次數為83次，平均最終流量誤差平方和從2 786 117.38(m3/s)2下降到112(m3/s)2，參數和誤差平方和的收斂速度非?？?，率定效果很好。

表3 理想模型率定結果

3 實際案例應用

洈水流域地跨湖北省、湖南省，流域面積954.2 km2，屬亞熱帶過渡性季風氣候區。烏溪溝水文站位于洈水水庫中上游，由于該水文站以上的子流域受人類活動影響較小，流域資料的可靠性性、一致性、代表性比較好。

本次選取烏溪溝1974年～2011年43場洪水過程作為次洪資料。其中，1974年～2007年的31場洪水資料作為率定期資料，2008年～2011年12場洪水資料作為檢驗期資料。

3.1 烏溪溝流域次洪參數線性化率定

用新安江模型參數線性化率定方法對烏溪溝流域的次洪模型分水源參數(SM、KI、KG)、坡面匯流參數(CS、CI、CG)、河道匯流參數(KE、XE)進行統一率定，率定步驟和結果如下：

(1)在給定范圍內隨機給上述賦初值(見表4)，分別隨機取20組參數初值進行計算，且保持KI+KG=0.7。

表4 參數初值范圍

(2)根據不同的參數初值分別進行參數線性化率定，各參數的迭代步長收斂容差均取0.001，即當以上各參數同時滿足迭代步長容差時，率定結束。

(3)參數率定結果見表5。從表中可以發現，各參數取不同初始值，經過新安江模型參數線性率定方法均可以收斂到固定的優值附近，這些參數優值均在其物理范圍內。參數平均率定次數為163次，收斂速度很快，效果很好。次洪率定最終參數見表6。

表5 烏溪溝流域參數率定成果

表6 烏溪溝流域次洪參數

(4)次洪模擬結果見表7，31場洪水參與模擬全部合格，合格率達到了100%。平均實測徑流深為98.9 mm，平均計算徑流深為103.3 mm，徑流深平均相對誤差為2.8%。平均實測洪峰為996.3 m3/s，平均計算洪峰為927.8 m3/s，平均洪峰相對誤差為-5.9%，平均峰現時差提前0.2 h，平均確定性系數為0.938。其中，2次洪水模擬過程見圖1、2。