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

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

0 引 言

流域水文模型在進(jìn)行水文規(guī)律研究和解決實(shí)際生產(chǎn)問題中起著重要作用。中國現(xiàn)在最常用的流域水文模型是新安江模型[1]，模型結(jié)構(gòu)考慮了濕潤半濕潤地區(qū)的蒸散發(fā)、產(chǎn)匯流特性。采用新安江模型建立實(shí)際流域的洪水預(yù)報(bào)方案時(shí)，首先要根據(jù)流域的實(shí)際水文資料對(duì)新安江模型參數(shù)進(jìn)行率定，然后根據(jù)率定的參數(shù)制定洪水預(yù)報(bào)方案，對(duì)未來的洪水進(jìn)行預(yù)報(bào)[2]。

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

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

1 新安江模型參數(shù)線性化率定方法

1.1 新安江模型介紹

趙人俊教授于1973年在編制新安江入庫洪水預(yù)報(bào)方案時(shí)提出了新安江模型。新安江模型結(jié)構(gòu)是分散性的，主要分為蒸散發(fā)計(jì)算、產(chǎn)流計(jì)算、分水源計(jì)算、和匯流計(jì)算四個(gè)層次。模型各層次結(jié)構(gòu)的功能，計(jì)算采用的方法和參數(shù)見表1。

表1 新安江模型各層次結(jié)構(gòu)功能、計(jì)算方法和相應(yīng)參數(shù)

1.2 新安江模型參數(shù)線性率定原理

參數(shù)線性化率定方法對(duì)非線性參數(shù)函數(shù)以參數(shù)作為自變量求導(dǎo)，再通過導(dǎo)函數(shù)差分線性化，并對(duì)線性化的參數(shù)用誤差平方和目標(biāo)函數(shù)進(jìn)行率定，然后逐步逼近非線性參數(shù)的最優(yōu)值。該方法收斂、率定速度快，解決了以誤差平方和為目標(biāo)函數(shù)求解非線性函數(shù)參數(shù)增加不相關(guān)的局部優(yōu)值的問題。與傳統(tǒng)的方法相比，該方法的實(shí)用性更強(qiáng)、效果更好。

新安江模型結(jié)構(gòu)完整，計(jì)算清晰，是一個(gè)典型的非線性模型。把新安江模型看作一個(gè)非線性函數(shù)，模型的計(jì)算流量當(dāng)作該函數(shù)的因變量，模型的參數(shù)當(dāng)成該函數(shù)的自變量。首先，對(duì)新安江模型參數(shù)變量求導(dǎo)，通過導(dǎo)函數(shù)差分使參數(shù)線性化；然后，把新安江模型這個(gè)非線性函數(shù)進(jìn)行一階泰勒級(jí)數(shù)展開，構(gòu)建計(jì)算流量和模型參數(shù)之間的線性函數(shù)，通過最小二乘法求解參數(shù)。把非線性參數(shù)函數(shù)轉(zhuǎn)變?yōu)榫€性參數(shù)函數(shù)求解可以避免產(chǎn)生一些不相關(guān)的局部參數(shù)解。以流量誤差平方和、參數(shù)迭代步長收斂容差為目標(biāo)函數(shù)進(jìn)行率定，逐步逼近新安江模型參數(shù)全局最優(yōu)值。

1.3 新安江模型參數(shù)線性率定步驟

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

(1)

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

具體率定步驟：

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

(2)根據(jù)給定的參數(shù)向量、降雨資料、蒸發(fā)資料用新安江模型計(jì)算防洪斷面流量q(θ，x)和參數(shù)靈敏度矩陣[5]S，S反映了參數(shù)改變引起計(jì)算流量的改變程度。把S代入函數(shù)式(1)得函數(shù)式(2)，計(jì)算S時(shí)需要計(jì)算函數(shù)對(duì)各參數(shù)的偏導(dǎo)數(shù)，通過導(dǎo)函數(shù)差分使參數(shù)線性化，計(jì)算公式見(式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)，得到新的參數(shù)向量Qj+1和參數(shù)尋找方向Δθ。即

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

(6)

Δθ=Qj+1-Qj

(7)

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

(8)

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

(9)

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

要證明步驟(3)中尋找方向的正確性，只要證明對(duì)于任意一步尋找的參數(shù)向量θj+1對(duì)應(yīng)的流量誤差平方和Fj+1比上一步參數(shù)向量θj對(duì)應(yīng)的流量誤差平方和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)

由此可見，參數(shù)率定過程中的流量誤差平方和是遞減的。隨著循環(huán)計(jì)算次數(shù)的增加，其對(duì)應(yīng)的流量誤差平方和越來越小，最終趨向于最小值，參數(shù)也逐漸逼近最優(yōu)的參數(shù)值。

1.4 新安江模型參數(shù)線性率定指標(biāo)

在模型參數(shù)的優(yōu)選中，目標(biāo)函數(shù)的選取直接影響參數(shù)率定的效率和結(jié)果。本文根據(jù)對(duì)模型成果的要求，選擇的子目標(biāo)函數(shù)主要有5個(gè)。

(1)計(jì)算流量和實(shí)測(cè)流量的誤差平方和

(13)

式中，Qx為計(jì)算的流量值，m3/s；Q0為實(shí)測(cè)流量值，m3/s。

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

(3)徑流深的相對(duì)誤差

(14)

式中，Rc(i)為計(jì)算徑流深，mm；Ro(i)為實(shí)測(cè)徑流深，mm。

(4)洪峰的相對(duì)誤差

(15)

式中，Qmax，c為計(jì)算洪峰，m3/s；Qmax，o為實(shí)測(cè)洪峰，m3/s。

(5)確定性系數(shù)DC，主要用來評(píng)價(jià)洪水預(yù)報(bào)過程和實(shí)測(cè)流量過程之間的擬合程度，確定性系數(shù)越大，表明計(jì)算流量和實(shí)測(cè)流量過程擬合得越好。即

(16)

2 理想流域驗(yàn)證

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

表2 參數(shù)初值取值范圍

針對(duì)新安江模型SM、KI、KG、CS、CI、CG、KE、XE統(tǒng)一進(jìn)行新安江模型線性率定，檢驗(yàn)該方法在高維數(shù)的參數(shù)空間的應(yīng)用效果。具體的率定步驟如下：

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

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

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

表3 理想模型率定結(jié)果

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

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

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

3.1 烏溪溝流域次洪參數(shù)線性化率定

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

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

表4 參數(shù)初值范圍

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

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

表5 烏溪溝流域參數(shù)率定成果

表6 烏溪溝流域次洪參數(shù)

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

圖1 31830703次洪水模擬過程

圖2 31800530次洪水模擬過程

表7

3.2 烏溪溝流域次洪模型檢驗(yàn)

烏溪溝流域次洪參數(shù)率定選取2008年～2011年之間12場(chǎng)洪水作為檢驗(yàn)洪水，計(jì)算結(jié)果見表8，12場(chǎng)洪水全部合格，合格率為100%。其中，2次洪水模擬過程見圖3、4。

表8 次洪檢驗(yàn)結(jié)果

圖3 31080827次洪水模擬過程

綜合評(píng)定表明，31場(chǎng)洪水參與模擬，全部合格，合格率為100%；12場(chǎng)洪水參與檢驗(yàn)，全部合格，合格率為100%，平均確定性系數(shù)為0.923。根據(jù)GB/T 22482—2008《水文情報(bào)預(yù)報(bào)規(guī)范》，預(yù)報(bào)等級(jí)為甲等。

圖4 31090627次洪水模擬

4 總 結(jié)

本文把包為民教授提出的參數(shù)線性化率定方法應(yīng)用于新安江模型，建立了新安江模型參數(shù)線性化率定方法，并通過理想流域應(yīng)用和實(shí)際流域應(yīng)用驗(yàn)證了該方法的的可行性和高效性。該方法具有以下優(yōu)點(diǎn)：

(1)新安江模型參數(shù)線性化率定方法在模型函數(shù)曲面上求解參數(shù)，解決了在流量誤差平方和目標(biāo)函數(shù)曲面上求解參數(shù)時(shí)增加不相關(guān)的局部優(yōu)值的問題。

(2)新安江模型參數(shù)線性化率定方法用一階泰勒級(jí)數(shù)函數(shù)以參數(shù)作為自變量展開，用最小二乘法確定參數(shù)的尋找方向，用參數(shù)迭代步長收斂容差和流量誤差平方作為目標(biāo)函數(shù)，尋找路徑始終在目標(biāo)函數(shù)曲面的低值槽附近，尋找過程的路徑短，最終找到的參數(shù)為參數(shù)的全局最優(yōu)值。

(3)新安江模型參數(shù)線性化率定中參數(shù)的收斂速度快，不受參數(shù)空間維數(shù)和參數(shù)初始值的影響，參數(shù)率定結(jié)果穩(wěn)定。取不同個(gè)數(shù)、不同的初始值的參數(shù)組合，經(jīng)過線性化率定均可以較快地收斂到穩(wěn)定的參數(shù)最優(yōu)值附近。

綜上所述，新安江模型參數(shù)線性化率定方法合理有效，計(jì)算思路簡單明確，具有較好的實(shí)際應(yīng)用效果。