劉希琛，闞光遠，丁留謙，何曉燕，梁 珂

(1.流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；2.中國水利水電科學(xué)研究院，北京 100038；3.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038；4.水利部京津冀水安全保障重點實驗室，北京 100038；5.北京中水科工程集團有限公司，北京 100048)

1 研究背景

概念性水文模型用于模擬天然流域中的復(fù)雜水文現(xiàn)象，流域特征通過模型參數(shù)進行描述。由于并非所有的模型參數(shù)都可由野外觀測、實驗室測量等途徑獲取，因此大多數(shù)模型參數(shù)通常由實測降雨—徑流資料反推獲得?；谥悄軆?yōu)化算法的參數(shù)優(yōu)化方法成為率定水文模型參數(shù)的主流方法[1]。

自20世紀(jì)二三十年代開始，隨著產(chǎn)匯流理論的日趨成熟和計算機技術(shù)的逐漸進步，流域水文模型迅速發(fā)展[2-4]。迄今為止，新安江模型、TOPMODEL、SWAT、MIKE SHE、HEC-HMS、WEP等模型已廣泛應(yīng)用于流域水文模擬[5-9]。新安江模型是我國應(yīng)用最廣的水文模型[10]。進行新安江模型參數(shù)優(yōu)化時，需要反復(fù)運行模型[11]，利用優(yōu)化算法在可行域內(nèi)尋找較優(yōu)的模型參數(shù)組合，進而獲取最優(yōu)模型參數(shù)。傳統(tǒng)的優(yōu)化方法在尋優(yōu)時，僅考慮模型參數(shù)的取值范圍，并未考慮模型參數(shù)的其他約束條件[12]。以新安江模型為例，無約束條件下的參數(shù)優(yōu)化存在兩個問題：一是長歷時洪水模擬中，當(dāng)產(chǎn)流和蒸散發(fā)參數(shù)設(shè)置不當(dāng)時，易于發(fā)生土壤含水量出現(xiàn)負值的現(xiàn)象。以往研究者針對這個問題的常見解決方案是將模擬中出現(xiàn)的負值土壤含水量重置為零，這種處理相當(dāng)于引入了多余的水量，會導(dǎo)致水量不平衡并降低模擬精度；二是新安江模型的匯流參數(shù)CG、CI、CS具有明確的物理意義，他們的大小關(guān)系反映了不同水源的匯流速度，但通常的參數(shù)優(yōu)化中并未考慮匯流參數(shù)取值的大小關(guān)系，可能得出不合理的徑流模擬結(jié)果。綜上所述，在水文模型參數(shù)優(yōu)化中施加約束是十分必要的，需要開展相關(guān)研究。

本研究嘗試解決新安江模型參數(shù)優(yōu)化過程中的約束條件問題，基于新安江模型原理、參數(shù)物理意義和SCE-UA算法，設(shè)計了數(shù)值模擬實驗，探索模型參數(shù)組合與土壤含水量正負值之間的關(guān)系，在目標(biāo)函數(shù)中構(gòu)建了罰函數(shù)機制，用于處理土壤含水量出現(xiàn)負值的情況，通過罰函數(shù)引導(dǎo)參數(shù)優(yōu)化算法尋找非負土壤含水量對應(yīng)的模型參數(shù)組合。匯流參數(shù)方面，通過設(shè)計罰函數(shù)確保匯流參數(shù)滿足CG≥CI≥CS的大小關(guān)系，引導(dǎo)優(yōu)化算法得到合理的匯流參數(shù)組合。數(shù)值實驗結(jié)果表明，提出的約束SCE-UA方法能夠使優(yōu)化算法獲得的模型參數(shù)更為合理，避免土壤含水量出現(xiàn)負值和匯流參數(shù)大小關(guān)系錯誤的不合理現(xiàn)象，能夠取得更好的水量平衡效果和更高的模擬精度。

2 考慮約束條件的新安江模型參數(shù)優(yōu)化方法

2.1 新安江模型參數(shù)1980年代，趙人俊等提出了三水源新安江模型[13]，該模型廣泛應(yīng)用于我國的濕潤與半濕潤地區(qū)洪水模擬與預(yù)報工作中[14]。本研究基于三水源新安江模型，按照蓄滿產(chǎn)流概念計算降雨產(chǎn)生的總徑流量，將徑流成分劃分為飽和地面徑流、壤中水徑流和地下水徑流。河網(wǎng)匯流采用滯后演算法，河道演進采用馬斯京根法。新安江模型需要優(yōu)化的參數(shù)為16個。如表1所示。

表1 新安江模型參數(shù)及其取值范圍

本研究選擇呈村流域作為研究區(qū)域，流域地處亞熱帶季風(fēng)氣候區(qū)，徑流年內(nèi)、年際變化較大，為典型的濕潤地區(qū)，適用新安江模型。事實上，呈村流域是新安江模型誕生地，模型最初的應(yīng)用就是從這里開始的，對于本研究來說，呈村流域是較為合適的選擇。該流域總面積為298 km2，流域內(nèi)共有呈村、汪村、左龍、田里等十個雨量站。流域出口水文站為呈村站。由于流域面積較小，可設(shè)置各子流域到出口斷面的馬斯京根分段數(shù)為0，因此馬斯京根法參數(shù)X不參加優(yōu)化。為了避免參數(shù)相關(guān)性對參數(shù)優(yōu)化結(jié)果的穩(wěn)定性帶來負面影響，引入結(jié)構(gòu)性約束KG+KI= 0.7。最終需要優(yōu)化的參數(shù)個數(shù)為14個。

2.2 土壤含水量與模型參數(shù)關(guān)系分析不合理的參數(shù)組合會導(dǎo)致長系列模擬中蒸發(fā)過大，進而發(fā)生土壤含水量出現(xiàn)負值的現(xiàn)象。土壤含水量是新安江模型產(chǎn)流計算模塊的輸出，因此，在降雨資料和流域初始狀態(tài)變量固定的前提下，土壤含水量的計算結(jié)果必然與模型的產(chǎn)流參數(shù)密切相關(guān)。產(chǎn)流參數(shù)中，IM的取值相對固定，對土壤含水量影響較小。K反應(yīng)了流域蒸散發(fā)能力，K值越大蒸發(fā)越強烈，土壤含水量越??；反之，土壤含水量越大。余下的參數(shù)中，WM、B和C相互制約，對土壤含水量的影響較為復(fù)雜，故而選擇這三個參數(shù)作為研究對象。固定其他參數(shù)，對WM、B和C進行隨機采樣并運行新安江模型進行數(shù)值模擬，觀察土壤含水量是否出現(xiàn)負值。參考其他學(xué)者在呈村流域的研究成果[15]，新安江模型的參數(shù)取值范圍見表2所示。

表2 呈村流域新安江模型參數(shù)范圍

由表1展示了參數(shù)的物理意義，WM為流域張力水容量，反映了流域包氣帶總體厚度，該參數(shù)可以分解上層、下層和深層張力水容量，即WUM、WLM和WDM；C為深層蒸散發(fā)系數(shù)，該值取決于深根植物的覆蓋面積；B為張力水蓄水容量曲線指數(shù)，該值反映了流域張力水蓄水條件的不均勻性。

采用拉丁超立方體抽樣算法[15]對WM、WUM、WLM、C和B進行隨機抽樣，將抽樣結(jié)果與其他參數(shù)組成樣本參數(shù)組，將樣本參數(shù)代入新安江模型并記錄模擬結(jié)果中土壤含水量是否出現(xiàn)負值(如出現(xiàn)負值，則結(jié)果為-1，反之，結(jié)果為1)，通過觀察參數(shù)組合的分布探索土壤含水量模擬值的規(guī)律。按照以上方法，進行10萬次洪水過程隨機模擬，參數(shù)取值分布與土壤含水量的正負值見下圖(其中橙色點據(jù)為土壤含水量出現(xiàn)負值的參數(shù)組，綠色點據(jù)對應(yīng)土壤含水量始終不出現(xiàn)負值的參數(shù)組)。

由圖1可知，在降雨資料和流域初始狀態(tài)變量固定的前提下，WM大于171.16 mm(圖中黑線位置)時，模擬結(jié)果的土壤含水量就不會出現(xiàn)負值；C越大，則土壤含水量出現(xiàn)負值的可能性越大(即橙色點據(jù)越密集)，這與C增大導(dǎo)致深層蒸散發(fā)量變大加劇土壤蒸發(fā)的物理情景相一致；采樣結(jié)果表明參數(shù)B并未體現(xiàn)出明顯規(guī)律性。三個參數(shù)中，WM的規(guī)律最為明顯：WM越大，土壤含水量出現(xiàn)負值的可能性越小。這一規(guī)律可用于指導(dǎo)約束條件罰函數(shù)的構(gòu)建。

圖1 WM、B、C分布與土壤含水量關(guān)系Fig.1 Relationship map between WM，B，C distribution and soil water content

2.3 匯流參數(shù)物理意義分析CS、CI和CG分別反映了坡面匯流、壤中水徑流和地下水徑流的匯流速度，且消退系數(shù)與匯流速度呈負相關(guān)關(guān)系[14，16]。同一流域內(nèi)，坡面匯流速度≥壤中流匯流速度≥地下水匯流速度，因此，我們可以推斷，如果匯流參數(shù)的取值符合其物理意義，則必須滿足約束條件CG≥CI≥CS。這一關(guān)系用于設(shè)計參數(shù)優(yōu)化過程中模型匯流參數(shù)的罰函數(shù)。

2.4 約束條件罰函數(shù)的構(gòu)建為了將約束條件應(yīng)用于模型參數(shù)優(yōu)化，本研究將約束條件以罰函數(shù)的形式引入SCE-UA參數(shù)優(yōu)化算法的目標(biāo)函數(shù)中。在參數(shù)優(yōu)化的目標(biāo)函數(shù)計算過程中引入兩個罰函數(shù)，罰函數(shù)實現(xiàn)兩種功能：第一，令不符合約束條件的參數(shù)樣本在種群進化中失去競爭優(yōu)勢；第二，要反映出不合格的參數(shù)樣本違反約束條件的程度。根據(jù)這兩項要求，我們設(shè)計如下兩個罰函數(shù)：

(1)保證土壤含水量始終非負的罰函數(shù)

如果參數(shù)樣本計算的土壤含水量出現(xiàn)負值，則目標(biāo)函數(shù)值調(diào)整為：

f=λ+λ×(WMmax-WM)

(1)

式中：λ為懲罰因子，取值根據(jù)優(yōu)化目標(biāo)確定；WMmax為WM的上限，其余參數(shù)意義同前。

(2)保證約束條件CG≥CI≥CS的罰函數(shù)

如果參數(shù)樣本不滿足CG≥CI≥CS，則目標(biāo)函數(shù)值調(diào)整為：

f=λ+λ×(|min(0，CG-CI)|+|min(0，CI-CS)|)

(2)

3 結(jié)果及討論

3.1 徑流模擬結(jié)果分析本研究通過數(shù)值模擬實驗，驗證參數(shù)優(yōu)化的效果[17]。兩組實驗均基于研究流域1986—1999年共計14年的日降雨—徑流資料，率定新安江模型參數(shù)[18]。其中，一組實驗采用未引入約束條件的SCE-UA算法；另一組實驗采用引入雙約束條件的SCE-UA算法(簡稱約束SCE-UA算法)。每組實驗均進行20次參數(shù)優(yōu)化并對結(jié)果進行分析。需要注意的是：未引入約束條件的情況下，當(dāng)遭遇土壤含水量出現(xiàn)負值時將其置零，以確保優(yōu)化能夠連續(xù)進行。表3是未引入約束條件的20次優(yōu)化的統(tǒng)計結(jié)果，表4是引入約束條件之后的20次優(yōu)化的統(tǒng)計結(jié)果。

表3 未引入約束條件的SCE-UA算法優(yōu)化新安江模型參數(shù)統(tǒng)計結(jié)果

表4 約束SCE-UA算法優(yōu)化新安江模型參數(shù)統(tǒng)計結(jié)果

由表3與表4可知，引入約束條件之前，所有的模擬結(jié)果中都出現(xiàn)了土壤含水量小于0的情況，且第8、9、11、12、20組參數(shù)，出現(xiàn)CI≥CG的現(xiàn)象。引入約束條件之后，模擬結(jié)果和匯流參數(shù)組合均未出現(xiàn)不滿足物理意義的情況。

按照表3與表4設(shè)置新安江模型參數(shù)值，得到的模擬流量值與實測流量值的散點圖如圖2與圖3所示，(散點圖NO.1—NO.20對應(yīng)著表中模擬洪水NO.1—NO.20)，用皮爾森相關(guān)系數(shù)R來衡量未引入約束條件的SCE-UA算法與約束SCE-UA算法的性能，相關(guān)性越高，模擬徑流過程精度越高，算法性能越強。

圖2 引入約束條件前的水文模擬結(jié)果散點圖Fig.2 Scatter plot of hydrological simulation results without introducing constraint

圖3 引入約束條件后的水文模擬結(jié)果散點圖Fig.3 Scatter plot of hydrological simulation results with constraints

圖2與圖3記錄了兩種算法的皮爾遜相關(guān)系數(shù)R，表5統(tǒng)計了模擬徑流過程的NSE與R，計算了兩者的均值、方差和極差(最大值與最小值之差)。從表5中可知，兩種優(yōu)化算法的NSE均值分別為0.8406和 0.8400(接近《水文情報預(yù)報規(guī)范》甲級標(biāo)準(zhǔn))，R均值分別為0.9177與0.9173(極強相關(guān))，這表明兩組模擬徑流過程的精度較高且非常接近；兩種優(yōu)化算法的NSE方差與極差、R方差與極差都接近于0，這表明兩種算法模擬徑流的質(zhì)量非常穩(wěn)定，引入約束條件之后，優(yōu)化結(jié)果的可靠性不會受到影響。

綜上所述，約束條件不僅可以提供符合物理意義的模型參數(shù)集，而且不會降低優(yōu)化算法的性能，是一種可靠的優(yōu)化算法。

3.2 算法穩(wěn)定性分析歐氏距離衡量多維空間中兩個點之間的絕對距離，該值越接近于0，說明多維空間中兩個點越接近。以某個實驗組內(nèi)所有樣本的歐氏距離的均值與方差代表算法的穩(wěn)定性，其中歐氏距離的均值越小，代表樣本間的差異越??；歐氏距離的方差越小，代表樣本之間歐氏距離的差異越小。兩種算法的歐氏距離的統(tǒng)計值見表6。

表6 兩種算法歐氏距離統(tǒng)計值

如表6所示，約束SCE-UA算法歐氏距離的均值和方差都優(yōu)于SCE-UA算法，這說明約束SCE-UA算法的優(yōu)化結(jié)果之間更加接近，該算法有更優(yōu)秀的穩(wěn)定性。

為了進一步論證約束SCE-UA算法在穩(wěn)定性方面的優(yōu)勢，本文計算了各實驗組的參數(shù)方差，方差越小，表明算法越穩(wěn)定。兩種優(yōu)化算法對應(yīng)的優(yōu)化結(jié)果參數(shù)方差見表7。

表7 優(yōu)化獲得的模型參數(shù)的方差

表8 水文模擬結(jié)果徑流量相對誤差統(tǒng)計

由表7可知，除WUM、WDM與C以外，約束SCE-UA 算法的優(yōu)化結(jié)果參數(shù)方差均不大于傳統(tǒng)SCE-UA算法的優(yōu)化結(jié)果。其中，兩種優(yōu)化算法之間WM方差的差異尤其顯著，約束SCE-UA的WM方差僅為傳統(tǒng)SCE-UA的1.1%。我們對表 3、表4與表7的結(jié)果開展了分析，認為約束 SCE-UA 算法獲取的WUM、WDM與C方差偏大現(xiàn)象出現(xiàn)的可能原因如下：新安江模型參數(shù)中，WM反映流域包氣帶總體厚度，該參數(shù)在優(yōu)化過程中可被拆分為WUM、WLM和WDM，分別代表上層、下層和深層的張力水容量。由于模型結(jié)構(gòu)設(shè)計的原因，WUM、WLM、WDM和C存在明顯的“異參同效”現(xiàn)象，導(dǎo)致了其分布規(guī)律較為復(fù)雜，這可能是導(dǎo)致方差偏大的誘因之一。此外，本研究統(tǒng)計方差時，樣本點個數(shù)為 20，這一樣本容量相對較小，可能不足以反映所有模型參數(shù)分布的規(guī)律，這可能是導(dǎo)致方差偏大的誘因之二。

綜上所述，約束SCE-UA算法的穩(wěn)定性要優(yōu)于SCE-UA算法。造成這種現(xiàn)象的原因是約束條件引導(dǎo)算法在符合物理意義的參數(shù)空間進行優(yōu)化，得到的參數(shù)更為合理和穩(wěn)定。

3.3 徑流過程水量平衡分析為進一步驗證約束SCE-UA算法在水量平衡方面的優(yōu)越性，兩組實驗中均統(tǒng)計了徑流量相對誤差，引入約束條件前后水文模擬徑流量相對誤差統(tǒng)計結(jié)果如圖4所示。

圖4 引入約束條件前后水文模擬結(jié)果徑流量相對誤差箱型圖Fig.4 Boxplot of relative error of runoff before and after the introduction of constraints for hydrological simulation

由圖4可知，約束SCE-UA對應(yīng)樣本的REV中，存在兩個異常值，在進行統(tǒng)計時，已將異常值除去。通過圖表結(jié)果可知，SCE-UA與約束SCE-UA的REV均值分別為10.121和8.788，REV中位數(shù)分別為10.452與8.521，約束SCE-UA算法的水文模擬結(jié)果徑流量相對誤差的中位數(shù)和平均值明顯小于SCE-UA算法的。SCE-UA與約束SCE-UA的REV上下四分位差分別為1.785和0.518，極差分別為2.741和1.310，約束SCE-UA算法的水文模擬結(jié)果徑流量相對誤差的取值更穩(wěn)定。結(jié)果表明，與傳統(tǒng)SCE-UA算法

相比，約束SCE-UA算法沒有引入額外的水量，能保證水文過程的水量平衡且具有更高的模擬精度。

綜上所述，約束SCE-UA算法在保證優(yōu)化算法性能的同時，不但可以使模型參數(shù)和狀態(tài)變量更符合物理意義，還能保證模擬徑流過程的水量平衡，提升徑流模擬的精度，是一種更為可靠的參數(shù)優(yōu)化方法。

4 結(jié)論

基于以上研究，我們可以得出以下結(jié)論：

(1)由數(shù)值實驗可知，在降雨資料和流域初始狀態(tài)變量固定的前提下，WM對土壤含水量的正負影響很大，WM越大，土壤含水量出現(xiàn)負值的可能性越小，反之，則土壤含水量出現(xiàn)負值的可能性越大。

(2)本文提出了參數(shù)優(yōu)化過程中的兩個約束條件：①土壤含水量始終非負，針對“土壤含水量出現(xiàn)負值”問題，該約束條件提供了一種更好的解決方式，即要求產(chǎn)流參數(shù)WM足夠大；②匯流參數(shù)大小關(guān)系始終滿足CG≥CI≥CS，該約束條件反映出同一流域內(nèi)，坡面匯流速度≥壤中流匯流速度≥地下水匯流速度的關(guān)系。兩種約束條件分別反映了參數(shù)優(yōu)化過程中模擬徑流和匯流參數(shù)的物理意義。

(3)引入約束條件之后的SCE-UA算法較原SCE-UA算法有以下優(yōu)點：①水文模擬結(jié)果更合理，水量平衡效果更好，兩組數(shù)值實驗的REV統(tǒng)計結(jié)果表明，與傳統(tǒng)SCE-UA算法相比，約束SCE-UA算法沒有引入額外的水量，能保證水文過程的水量平衡且具有更高的模擬精度；②水文模擬過程和匯流參數(shù)的物理意義均得到滿足，模型參數(shù)優(yōu)化結(jié)果不會出現(xiàn)違反物理意義的不合理組合；③優(yōu)化結(jié)果的穩(wěn)定性更高，兩種算法歐氏距離統(tǒng)計結(jié)果表明，約束SCE-UA得出的樣本之間差異較小，穩(wěn)定性更強。