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

陆旻皎

（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］提出的“降维缩容”模型优化战略。通过设计不同的目标函数进行全局敏感度分析可以从不同角度分析参数空间结构，据此提出有效的参数率定方案应该是有效的途径。

另一方面，深化模型参数规律研究，减少需率定参数也是应该努力的方向。利用各种水文、气候、地理和地质等方面的数据、信息、知识、规律和法则，从模型物理概念出发理论上推估模型参数值或参数间关系是一个有效的途径。收集模型参数，形成一个覆盖全国的大数据，分析其区域规律，也是很有价值的工作。中国幅员广阔，覆盖了多个气候区，这些规律性不但适用于国内，应该可以适用于国外类似地区。可为缺资料或无资料流域的模型应用创造条件。

模型研发的最终目的就是效果良好且使用方便的模型。作为一个尝试，笔者开发了一个网页版新安江模型，用户只需注册，上传数据，就可以开始使用模型。一般在十几分钟内可以完成。不需要学习操作系统，不需要编程，大大降低了新安江模型使用的门槛。只要输入参数值，就可以执行模型。通过图示模拟结果和标示各种数值指标，使模型参数人工调试更为顺畅。作为参数率定支持功能，可估算部分参数值，推荐用户试用。某种意义上讲，我们的网页版新安江模型是一个反映和测试最新研究成果的一个平台。

在模型推广方面，目前还看不到一个官方渠道发布模型的软件、网页和各种说明书的官方版本。要在国外推广，英文说明书必不可少。除了说明书，设计一个好的新安江模型课件也会帮助新安江模型走进各个学校的课堂。一方面可以加强模型程序的中央控制，另一方面也可以提高模型公信力。同时减轻用户升级模型的负担。另外，学术界的模型应用已经越来越多样化，不仅仅局限于传统的洪水预报及水资源管理。迄今为止，已有很多水文研究人员做了各方面的工作。应该有一个平台把他们的成果整合到模型中去，形成一个能覆盖很多与水有关的应用领域的综合性模型。