新安江模型研究的回顧和展望

陸旻皎

（1.長岡技術科學大學，日本新潟縣長岡市 940-2188；2.重慶交通大學，重慶 400060）

1 研究背景

新安江模型是由河海大學（原華東水利學院）趙人俊教授所領導的團隊研發的降雨徑流模型，于1963年提出雛形，1973年正式發布［1-2］。該模型在濕潤、半濕潤地區得到廣泛應用。其應用范圍不僅僅局限于降雨徑流模擬及洪水預報，也包括水資源規劃管理、防災、農業等與水有關的領域。為國民經濟發展發揮了重大貢獻，被列入中華人民共和國重大科技成果選集［3］。與很多世界同類知名模型，如薩克拉門托模型［4］、TOPMODEL［5］、水箱模型［6］和HBV 模型［7］等相比，在精度上也不遜色［8-9］。

新安江模型在中國及世界各國的廣泛應用，使其結構設計在不同氣候和地理條件下得到考驗。圖1和圖2為2018年8月初在百度學術和谷歌學術上搜索相關中英文關鍵詞得到的期刊論文數。可見新安江模型相關研究是相當活躍的。學術研究如此，全國水文機構也廣泛使用新安江模型。從1999年舉行的全國水文預報技術競賽［9］參賽隊伍可見一斑。此次競賽共有25 隊10 個模型參賽，其中14 個隊選擇了新安江模型。新安江模型的應用流域數應該是其他模型很難比肩的。這為其參數的區域規律研究創造了良好的條件。

新安江模型在國內水文單位取得了壓倒性的利用，但在學術層面的應用似乎有停滯或減少的趨勢。與國內期刊論文數近十年基本保持不變或微幅增加相比（圖1），英文論文數有比較明顯的下降趨勢（圖2）。是新安江模型已經趨于成熟，已經變成諸多模型中的一個，無需進一步研發呢。還是隨著社會需求的多樣化，不能滿足各類用戶的需求，不能吸引新用戶。抑或是傳統的集總式已經到了被分布式模型取代的時代。

圖1 百度學術新安江模型相關期刊論文數

圖2 谷歌學術新安江模型相關論文數

就分布式模型而言，在實際應用中的優勢并不明顯。日本所有主要流域都引進了分布式模型進行洪水預報，但精度等方面還有待改善。發揮分布式模型理論上的優勢需要高精度分布式降水資料的穩定獲取以及分布式水文參數的提取等方面的支撐。筆者認為在這些方面還有很多工作要做。因此，在今后較長的一段時間內分布式模型不可能完全取代傳統的集總式模型。模型用戶將會根據自己的需求做出選擇。作為實質標準水文模型，新安江模型將繼續發揮其作用。

本文試圖回顧有關新安江模型研發的相關研究，冀望從中得到啟發，展望模型發展的方向。如有不周之處，還請指正。

2 新安江模型研究的回顧

期待新安江模型研發更注重可用性和易用性，本文主要對模型核心概念、參數空間特性、參數率定及參數可移植性進行了回顧。

2.1 蓄滿產流和現代產流機制研究愛爾蘭工程師Mulvaney 在1851年發表的經驗公式［10］是第一個從雨量計算（洪峰）流量的公式。推理公式［11］雖然只是把式中的經驗常數替換成徑流系數，卻跨出了把產流和匯流區分處理的一大步。而1930年代初由Sherman 提出的單位線理論［12］和Horton 提出的滲透公式［13］分別提供了匯流和產流計算手段，自然形成了預報流量過程線的強有力組合，對現代水文科學，水文模擬產生了深遠的影響。全流域超滲產流的概念被廣泛接受。雨強大于下滲能力，即超滲成為產流的唯一機制，產生的地表流也被冠名為超滲地表流或Horton 地表流。Betson 于1964年證明即使是超滲產流也僅在流域的一部分發生，提出部分產流區（partial area）的概念［14］，才開始改變這種狀況。而山坡水文學研究得出的變動產流區（variable source area）的概念［15-17］對產流機制的認知產生了更大的影響。在此基礎上，Dunne 和Black 提出了飽和產流機制［18-19］，由此產生的地表流也被命名為飽和地表流。在濕潤地區，土壤發育較好，下滲能力大，超滲產流難以發生。另一方面，由于雨量豐沛，地下水面一般比較淺，加上地形作用地下水向河道附近匯集，在其附近形成飽和或接近飽和的區域，一旦降雨，這個區域就會成為或迅速成為產流區，而且會在降雨過程中擴大或縮小。TOP?MODEL［5］正是由Beven and Kirkby 提出的基于飽和產流機制的水文模型。

中國的水文規律及產流機制的研究比國外大約要遲20年。在大規模水利建設的推動下從1950年代開始迅速發展，至1960年代中期就獲得了一批重要的研究成果［20］。由于客觀條件的限制，中國學者主要對降雨徑流關系及流量過程線進行了非常深入的研究［2，21］，并提出蓄滿產流機制，也就是土壤水分滿足土壤蓄水容量以后才開始產流。更進一步把蓄水容量的空間分布統計歸納后投影到一條分布曲線——即蓄水容量曲線上，得以非常簡潔地表達產流量和產流區的時變性和非線性特征。Beven［22］將新安江模型恰如其分地分類為概率分布式模型（probability distributed model）［23］。蓄水容量曲線可以說是新安江模型的核心，其統計投影手法也被很多水文模型，如VIC 模型［24］、ARNO 模型［25］所借鑒，為考慮流域水文特性的異質性開創了一個新的途徑。

需要指出的是，蓄滿產流雖然與飽和產流非常相似，但還是有一些區別。超滲產流和飽和產流都是表述地表徑流的產流機制，而蓄滿產流則是描述總徑流量的產流機制。所謂蓄滿是指土層全層達到田間持水量，蓄滿以后進入土層的水分為自由水，最后在重力作用下以壤中流或地下水的形式流出。在新安江模型里，產流和分水源是分別通過引進張力水蓄水容量曲線和自由水蓄水容量曲線來處理的。

另一個潛在的問題是，蓄滿產流沒有考慮雨水進入土壤后的垂向運動，即下滲后的再分布［26］。濕潤地區一般植被良好，下滲能力很大，以日本為例，森林流域的穩定下滲率平均在200 mm/h以上［27］，很少出現超滲產流，而且下滲的水分也在很短的時間內完成再分布。加之當時模型運用以日模型為主，這樣的處理有其合理性。但是隨著數據采集技術的進步，水文資料的時間分辨率會提高。同時社會也會對洪水預報提出更高的要求，時模型的利用會逐漸增加。有必要對再分布的影響進行評估。

2.2 蓄水容量曲線如前節所述，蓄水容量曲線是新安江模型的核心。迄今為止已經有一些研究嘗試用地形、土壤和植被等信息推求蓄水容量［28-31］。這些研究的共同之處是，假設流域某一點或某一柵格的張力水蓄水容量（WM）和自由水蓄水容量（SM）為包氣帶內的所有張力水和腐殖質土層內的所有自由水，即有：

式中：θs、θf、θwp分別為飽和含水量、田間持水量和凋萎含水量；Lv、Lh分別為包氣帶厚度和腐殖質土層厚度。然后通過地形指數或土壤地形指數［5］等計算包氣帶和腐殖質土層厚度得到蓄水容量。也有部分研究嘗試從地形指數累積曲線直接推導蓄水容量曲線［32-33］。孔凡哲和宋曉猛提出用地表坡度推導蓄水容量曲線的方法［34］。值得特別提出的是，向小華等提出了一種考慮毛細管邊緣（capillary fringe）的蓄水容量［35］，改善了式（1）不考慮水分毛細管上升的缺陷，可望提高對細顆粒土壤的適用性。

上述研究不但可以幫助確定有關模型參數，更重要的是可以使模型的核心概念有更明確的物理意義。但筆者認為包氣帶厚度和蓄水容量不但取決于地形、土壤和植被，也應該受氣候條件的影響。在此整理一下基本物理概念，考慮如何引進氣候因素。

筆者定義張力水蓄水容量為最大包氣帶張力水缺水量，為：

式中：Lv，max為包氣帶最大厚度；θ為體積含水量；z為垂向坐標，以上方為正。包氣帶厚度是隨著時間變化的，其最大值代表該地點的干燥程度，也受氣候影響。由于最大厚度一般出現于很長的無降雨期的末期，含水量垂向分布比較穩定，可以近似地假設毛管張力ψ和重力處于平衡狀態，即：

如果用Brooks-Corey 模型：

表達土壤水分特征曲線，可得含水量如下：

式中：ψb，λ為Brooks-Corey 模型參數；θr為滯留含水量。田間持水量出現在：

對式（3）積分可得：

通過式（6），只要有土壤參數和地下水資料就可以估算張力水蓄水容量。

2.3 模型參數和率定新安江模型參數按照模型結構分為4 層：蒸散發計算參數；產流計算參數；分水源計算參數；匯流計算參數［36］。表1列出了三水源新安江模型相關參數及其物理意義，取值范圍和敏感度［36-37］。由于本文會引用筆者所建的網頁版新安江模型［38］，其敏感度［39］也列在表中。

筆者的網頁版新安江模型與原版三水源新安江模型的不同之處僅在于：引進雨量折算系數；省略子流域間匯流；省略壤中流和地下水的河網匯流。引進雨量折算系數是考慮多雨地區雨量誤差有可能導致蒸散發能力折算系數出現不合理數值［40］。省略子流域間匯流和省略壤中流和地下水的河網匯流是因為筆者的網頁版新安江模型主要目的是教學和研究，對象流域較小，一般不需要分子流域。而且河網匯流的時間尺度一般遠小于壤中流和地下水的時間尺度，可以忽略。另一方面，連結各子流域的河網匯流無論在信息提取和模型化都已經非常成熟，有很多水力學或水文學匯流模型可供選擇。劃分子流域在大流域的建模中帶來的不確定性相對較小。從教學和研究的觀點看，理解中小流域內的水文規律相對比較重要。

兩個模型有14 或15 個模型參數。雖然這些參數都有明確的物理意義，參數率定即使對專業的水文工作者來說也非易事。在1999年舉行的全國水文預報技術競賽中，新安江模型的14 個隊在同一個流域用同一資料得到的參數和模擬結果還是存在較大的差異［41］。專業人員尚且如此，更何況跨專業模型用戶。怎么降低參數率定的不確定性是一個亟待解決的重要課題。

理解參數空間的結構對參數率定，無論是手動還是自動率定，都是非常重要的。所謂結構是指參數的取值范圍，參數的敏感度及參數間關聯性。趙人俊等［1，42］通過大量的模型率定歸納出表1所示的參數范圍，并發現4 個層次間的參數獨立性較好［43］。在參數敏感度方面，發現敏感參數只有7 個，分別為：K（蒸散發能力折算系數）、SM（表層土自由水蓄水容量）、KI（表層自由水壤中流日出流系數）、KG（表層自由水地下水日出流系數）、CG（地下水日消退系數）、L（河網蓄水滯后時間，表中未列出）和CS（地表水日消退系數）。

表1 三水源新安江模型參數及其取值范圍和敏感度

敏感度分析可以幫助理解參數空間的結構。由于水文現象的非線性性很強，應該盡可能使用反映參數空間整體情況的全局敏感度分析方法。Zhang 等［37］用普適似然度（GLUE）方法對新安江模型參數敏感度進行分析，指出對模型輸出最為敏感的3 個參數為K、SM和CS。K控制水量極為敏感，SM和CS影響流量過程線非常敏感。Lu 和Li［39］以實測及計算日流量計算年月日過程線的Nash 效率系數為目標函數，在年月日這3 個不同的時間尺度上對新安江模型參數用Morris 法進行了全局敏感度分析。結果表明，在不同的時間尺度不同的參數組合顯示其敏感度。在所有時間尺度上，Cp和Cep具有壓倒性敏感度。在固定Cp和Cep后，分水源和匯流參數，SM、KG、CS和CI（壤中流消退系數）在日尺度上比較敏感。

根據對參數空間結構的理解，就可以設計一套參數率定方案。趙人俊等的方案［42］是：以水量平衡為目標函數率定蒸發及產流層的參數，以過程線絕對誤差為目標函數率定匯流層的參數，以過程線對數絕對誤差為目標函數率定分水源層的參數。而Li 和Lu［44］則把參數分成3 組：第一組包括兩個折算系數，第二組包括分水源及匯流的參數，第三組包括產流和蒸發的參數，然后依次優化這些參數。第一組和第三組擬合年過程線，第2 組擬合日過程線。其和趙人俊等的方案有不謀而合之處，具體見下節。

2.4 模型參數自動率定自從概念性水文模型問世以來，參數自動率定一直是模型開發者和用戶的愿望。也被認為是增加模型可用性和易用性的重要途徑。大多數自動率定以優化算法為基礎。隨著計算機性能提高，涌現出一批出色的優化算法，如遺傳算法、粒子群算法和模擬退火算法等等。正如Duan 等［45］指出的那樣，即使是一僅有6 個參數的簡化水文模型，其目標函數的響應曲面已非常復雜，有很多局部極值。這種復雜性的主要原因包括參數間相互作用， 響應曲面的非凸性和微分的不連續性。因此，在早期的模型優化研究中，較多使用Rosenbrock 法和單純形法等直接尋優的方法，而不太使用基于梯度的方法。從Duan 等的結果不難知道這些局部優化方法很難搜到全局最優。也就是說這些方法可以或多或少改善模型參數，卻很難找到最好的參數組合。最近的參數優化研究較多使用全局優化算法［46-49］，也顯示了全局優化算法的優越性。特別是SCE-UA（或SCEM-UA）是一個考慮水文模型特性的全局優化算法，在國內外得到廣泛利用［50-51］。

毫無疑問，相關研究對新安江模型自動率定是非常有意義的。但必須指出，即使是全局優化算法，也只是提高探索得到全局最優的概率。不應過分仰賴優化算法，而應該對優化結果進行充分的論證。筆者曾經用SCEM-UA 法對表1中15 個參數進行優化，雖然每次優化都能得到一組不錯的參數組合，但每次都不同。而且每次調用模型50 萬次都達不到收斂標準。Li 和Lu［44］根據敏感度分析結果，提出了多步優化方案。把參數分成3 組：第一組包括兩個折算系數；第二組包括分水源及匯流的參數；第三組包括產流和蒸發的參數。從一組常用參數值出發，依次用SCEM-UA 法優化第一組，第二組和第三組參數。前面的優化結果在后續優化中使用。第一組和第三組最大化年過程線效率系數，第2 組擬合最大化日過程線效率系數。為了確認方案的收斂性，筆者也曾經循環這三組參數的優化，第二次、第三次結果變化不大，也就是說不用循環就可以得到較好的結果。在理想條件下，即使用人工生成流量資料時，方案基本能夠得到生成流量時使用的參數。在使用實際資料時，也能得到比較滿意的結果，包括精度和參數的合理性。收斂速度第一組最快，第二組次之，第三組最慢。模型調用次數一般在10 萬到15 萬之間，遠遠低于設定的優化程序模型調用上限50 萬次。

根據Li 和Lu［44］的結果，Lu 和Li［52］提出了“降維縮容”的模型優化戰略。“降維”是通過理解參數空間的結構把高維空間分割成幾個低維子空間。在這些子空間上進行優化，可以提高優化的收斂性和得到全局最優的可能性。“縮容”分兩個部分。一個是在不失一般性的情況下盡可能縮小參數取值范圍，使參數空間容積縮小。如利用退水曲線推估地下水消退系數，然后據此設定比較小的取值范圍。另一個是考慮參數間的相互關系限制取值范圍，如設定CI ＜CG。雖然這會限制理論上可能出現而實際上很少出現的壤中流出流慢于地下水的現象，但卻可以縮小參數空間，并避免新安江模型結構上潛在的異參同效［53］。從模型結構不難發現，兩組參數值和互換不會對模型結果產生任何影響。因此在CI和CG的取值范圍有重疊時，參數空間內會出現完全一模一樣的兩個子空間，使絕大多數優化算法面臨困難。

另一個在文獻中很少被提到的是模型初始值問題。優化過程中，每次調用模型都必須設定初始值。他們對模型輸出的影響會持續一段時間。地表水短一些，幾小時到幾天，主要因流域大小而異。壤中流較長，從幾天到幾十天。地下水最長，為數月到數年［36］。土壤水分則在數月到一年左右，氣候越干燥影響時間越長［54］。計算目標函數時應盡可能剔除這段時間的模型結果，以避免優化參數值受影響。優化次洪模型時更應充分考慮初始值問題。

2.5 模型參數可移植性參數可移植性一方面是模型運用上的需求，也是對模型普適性的檢驗。日模型參數移植到小時進行次洪預報需要時域上的可移植性，資料短缺流域或無資料流域的洪水預報需要空域上的可移植性。可從時間，空間兩個方面討論［55］。

趙人俊等［42］指出新安江模型分水源和匯流相關參數受模型時段長影響。分水源參數受雨量時段內均化影響比較復雜，而匯流相關參數可以進行時段轉換。Jie 等［56］對這些參數和模型時段長作了深入的分析，并建立了參數值和時段長的經驗關系。陸旻皎［57］根據線性水庫理論得出：

式中：CΔt為時段長Δt的消退系數；T為線性水庫的蓄泄系數。并據此得到消退系數時段轉換公式：

舉例來說日消退系數是時消退系數的24 次方。因此在討論消退系數時必須明確說明時段長度。

理論上講，式（8）可以用來進行任何時段之間的轉換。但在長時段向短時段轉換時，長時段消退系數往往有很大誤差，不實用。陸旻皎［57］考慮雨量的時段內不均勻性，簡單討論消退系數所具有的不確定性，建議率定參數或進行時段轉換時確保轉換前時段長，也就是說CΔt＞0.82 以保證退水段有大約80%洪量。趙人俊和王佩蘭［36］指出地下水日消退系數一般為0.98 ～0.998，壤中流日消退系數經常達到0.9，可以用式（8）從日消退系數轉換得到時消退系數進行次洪計算。而河網蓄水消退系數決定于河網地貌等，匯流時間在中小流域一般為數十小時，日消退系數會有很大的不確定性。通過轉換得到時消退系數很困難。往往需要用次洪模型進行率定。總而言之，除了不同時段長度參數間的關系外，還必須考慮轉換前參數值所帶有的不確定性。

空域上的可移植性主要指模型參數的跨流域使用，與參數區域規律的內涵有很大重疊。對解決缺資料或無資料地區水問題具有重要意義。

歷代水文學家都非常重視這個問題，做了大量的工作。趙人俊等［42］，陳志明［58］詳盡分析了SM與地質條件，巖性的關系，為無資料地區確定SM推求提供了一個途徑。趙人俊［59］在1991年提出了CS隨時間、流量變化的方程。徐倩等［60］對黃山地區13 個流域進行參數率定，得出了CS與流域面積的經驗關系。陸旻皎［57］從線性水庫理論出發得到式（7）所示消退系數和蓄泄系數T的關系，并提出蓄泄系數與流域面積的經驗關系。在2.2 小節提及的研究也為從數值地理信息提取蓄水容量曲線提供了思路。另一方面夏自強［61］提出流域平均蓄水容量與干旱指數有關，顯示了氣候條件的影響。

Li 和Lu［62］結合年水量平衡方程式和Budyko 假設，提出了用年雨量、徑流量及蒸發能力推算Cep的關系如下：

式中：R為年徑流量；Pg為年雨量；ζ為資料干旱指數，即資料年平均蒸發量與年平均雨量之比；α為常數。式（9）在優化時可以用來減少優化參數。如果雨量不需要折算，即Cp =1 時，

就可以直接估算網頁版新安江模型的Cep，即原版新安江模型的K。這是一個很重要的參數。

利用各種水文、氣候、地理和地質等方面的數據，信息、知識、規律和法則推求模型參數是比較理想的路子。更多是用有資料流域的參數值和各種水文，氣候，地理，地質等方面的數據建立統計關系。如前所述，這方面新安江模型在應用流域數上有絕對優勢。每個流域的模型參數受率定者的主觀因素影響也許有很大的差異，但覆蓋全國的大數據一定會揭示中國境內模型參數的規律性。中國幅員廣闊，覆蓋了很多氣候區，這些規律性應該不僅僅適用于中國，對模型在全世界范圍的拓展也會有舉足輕重的影響。

3 新安江模型的未來展望

通過幾個重要方面回顧新安江模型相關研究，不難看出經過廣泛的應用，已經積累了大量的實例，經驗和知識，為新安江模型的新的發展和突破奠定了良好的基礎。

水文模型從根本上來說是水文規律認知的數理表現。沒有正確的認知，就沒有正確的水文模型。毫無疑問，水文規律研究是正確認知的根本。從水文現象的觀察、觀測、分析到抽象化，模型化，再到模型的應用和檢驗是一個周而復始，螺旋上升的過程。新資訊，如數字地理信息、遙測遙感資料，新技術都會提供新的視角和視野，允許我們更深入理解水文規律，最后反映到認知和模型中去。

水文規律研究當然也應該包括相關領域。舉例來說，流域產沙［63-64］及化學物質輸移［65］就和地表徑流的產流機制有密切關系；流域內生態環境［66］也與濕潤狀況、植被和土壤等息息相關。在這些領域的應用可以給水文模型提出新的需求，而來自相關領域的信息可以間接或直接地檢驗水文模型的合理性。受條件限制，水文學家經常僅能用流量的擬合程度來判斷模型的優劣。流量以外的檢驗手段因此具有特殊的重要意義。

從模型運用的角度來看，與其他現代水文模型相同，新安江模型有十多個參數而且相互影響，即使是水文工作者參數率定也非易事。依靠高速計算機和全局優化算法盲目優選參數也很難保證結果的最優性。如何正確高效率定參數就成為模型應用的一個門檻。無論是參數人工調試還是自動率定，理解參數空間的結構都是必不可少的。趙人俊等的方案［42］和Li 和Lu 提出的多步優化方案［44］都符合Lu 和Li［52］提出的“降維縮容”模型優化戰略。通過設計不同的目標函數進行全局敏感度分析可以從不同角度分析參數空間結構，據此提出有效的參數率定方案應該是有效的途徑。

另一方面，深化模型參數規律研究，減少需率定參數也是應該努力的方向。利用各種水文、氣候、地理和地質等方面的數據、信息、知識、規律和法則，從模型物理概念出發理論上推估模型參數值或參數間關系是一個有效的途徑。收集模型參數，形成一個覆蓋全國的大數據，分析其區域規律，也是很有價值的工作。中國幅員廣闊，覆蓋了多個氣候區，這些規律性不但適用于國內，應該可以適用于國外類似地區。可為缺資料或無資料流域的模型應用創造條件。

模型研發的最終目的就是效果良好且使用方便的模型。作為一個嘗試，筆者開發了一個網頁版新安江模型，用戶只需注冊，上傳數據，就可以開始使用模型。一般在十幾分鐘內可以完成。不需要學習操作系統，不需要編程，大大降低了新安江模型使用的門檻。只要輸入參數值，就可以執行模型。通過圖示模擬結果和標示各種數值指標，使模型參數人工調試更為順暢。作為參數率定支持功能，可估算部分參數值，推薦用戶試用。某種意義上講，我們的網頁版新安江模型是一個反映和測試最新研究成果的一個平臺。

在模型推廣方面，目前還看不到一個官方渠道發布模型的軟件、網頁和各種說明書的官方版本。要在國外推廣，英文說明書必不可少。除了說明書，設計一個好的新安江模型課件也會幫助新安江模型走進各個學校的課堂。一方面可以加強模型程序的中央控制，另一方面也可以提高模型公信力。同時減輕用戶升級模型的負擔。另外，學術界的模型應用已經越來越多樣化，不僅僅局限于傳統的洪水預報及水資源管理。迄今為止，已有很多水文研究人員做了各方面的工作。應該有一個平臺把他們的成果整合到模型中去，形成一個能覆蓋很多與水有關的應用領域的綜合性模型。