Skip to content

Home

Jump to bottom
OuyangWenyu edited this page May 22, 2019 · 4 revisions

Welcome to the hydro-model-xaj wiki!

hydro-model-xaj

新安江水文模型 下面各个公式的推导过程打出来还是有些麻烦，这里就主要以文字为主进行介绍了，如果有人 star 或者 fork 这个项目的话，我就把下面算法里的公式的推导过程打出来。

本项目有了第一个相对完整的可测试版本v0.1.1，仍在开发中。。。

项目介绍

根据中国水利水电出版社出版的，由武汉大学叶守泽老师和河海大学詹道江老师等合编的《工程水文学》第三版教材中新安江模型的原理，结合河海大学芮孝芳老师《水文学原理》中的相关知识，并重点参考了河海大学包为民老师的《水文预报》第4版，编写的新安江三水源模型的 python 版本。

使用说明

运行test.py中的测试函数即可，目前测试函数包括直接调用模型的测试函数，以及率定的测试函数，简单修改程序即可测试任意一个。测试使用的数据较少，并不符合实际预报规范，目前项目只是为了梳理新安江模型原理及其运算过程，实际应用还在本项目代码基础上进一步开发。

目前文档使用markdown编写，浏览器上显示格式不完全，建议使用vscode编辑器+相应插件后浏览。

如果觉得项目中文字和代码对原理的理解有误，或者出现代码运行错误，请在issues中留言。

主要内容

数据处理：

首选txt 格式，适合 nc 格式的地方也尝试使用 nc 数据格式，比如有经纬度的时间序列数据——雷达测雨数据等。原因可以参考：Python programming guide for Earth Scientists。nc 格式数据 python 接口官方资料：netcdf4-python，补充一些资料如下：Generating NetCDF files with Python；Python-NetCDF reading and writing example with plotting。

新安江模型核心算法

前期土壤含水量计算

在进行产流计算之前，首先需要计算前期土壤蓄水容量。利用前期降雨蒸发数据，采用流域蒸发计算模型逐日计算。

  • 当降雨 P 大于蒸发 E，即 P-E=PE>0。 说明流域蓄水量增加，这部分增加的水量在流域蒸发模型中，首先被用来补充上层土壤水，在流域上层蓄水量达到最大值后，再补充下层，同理最后补充深层。在判断表层水是否补充充足时，需要计算产流量 R，产流量计算遵循蓄满产流原理（见下一小节）。在整个过程中，只有表层有蒸散发，且按蒸散发能力蒸发，K 是折算系数，EM 是实测的水面蒸发值，下层和深层为 0。
  • 当 PE<0 时。 说明流域失水，那么首先蒸发上层水，按上层蒸发能力蒸发，上层蒸发完之后再消耗下层的含水量，下层的蒸发量与下层含水量成正比，与蓄水容量成反比，在两层蒸发模型中根据剩余蒸发量，按比例直接计算可得。但是，当下层含水量与下层蓄水容量之比小于深层蒸散发系数 C 时，需要利用三层蒸发模型。

上面一段蒸发为什么这么计算，书上没有解释，个人理解： 首先蒸散发的对象是张力水W，张力水在模型里指的应该是田间持水量以下的水分，即毛管水。产流主要是自由水，模型里的自由水主要是指重力水和径流，自由水被概化为流域地下自由水蓄水水库，有蓄水和下泄。在蒸散发中，对张力水分析，下层含水量与下层蓄水容量之比小于深层蒸散发系数 C，深层蒸散发系数可以理解为深层蒸发量和深层蒸发能力之比，同样假设其等于深层含水量与深层蓄水容量之比，那么下层含水量与下层蓄水容量之比理解为下层蒸散发系数，也就是说下层蒸散发系数小于深层时，深层水应该承担一部分的蒸发量。但是如果下层也按深层的蒸散发系数蒸发，可以承担剩余的蒸发量，那就不需要深层水了。这里推理一下，应该是认为下层水的蒸散发系数不是定值，看下层含水量情况以及和深层蒸散发系数比较的情况，因为C$\times$WLM标识了毛管断裂含水量，因此当下层含水/下层蓄水容量 大于深层蒸散发系数时，可以认为下层土壤的水分输移是稳定的，下层含水量是可以一直有深层补充的，即下层和深层是关联的，因此采用二层蒸发模型即可；或者就算小，但是因为C$\times$(EP-EU)标识了土壤水的扩散能力，按照深层蒸散发系数蒸发也能满足除去表层蒸发之后的蒸发量时，同样可以认为深层水和下层水是联系的；如果不满足这两种情况，那就说明深层水到下层水的输移已经打断，需要深层水独立分担一部分了。

最后汇总三层蒸发和三层蓄水。

整体计算流程图如下所示：

  • PE>=0 时
st=>start: 执行
sub=>subroutine: 计算R
cond2=>condition: WU+PE-R>WUM
op2=>operation: WU=WU+PE-R
cond3=>condition: WU+PE-R-WUM+WL>WLM
op3=>operation: WU=WUM,WL=WLM,WD=W+PE-R-WU-WL
op4=>operation: WU=WUM,WL=WU+WL+PE-R-WUM
op5=>operation: EU=K·EM,ED=EL=0
io=>inputoutput: E=EU+EL+ED,W=WU+WL+WD
e=>end: 结束

st->sub->cond2
cond2(yes)->cond3
cond2(no)->op2->op5
cond3(yes)->op3->op5
cond3(no)->op4->op5
op5->io->e
  • PE<0 时 以下PE均是PE的绝对值，即E-P
st=>start: 执行
cond4=>condition: WU>PE
op7=>operation: EU=K·EM,EL=ED=0,WU=WU+PE
op12=>operation: .
op8=>operation: EU=WU+P,WU=0
cond5=>condition: WL>C·WLM
op9=>operation: EL=(K·EM-EU)·WL/WLM,WL=WL-EL,ED=0
cond6=>condition: WL>C·(K·EM-EU)
op10=>operation: EL=C·(K·EM-EU),WL=WL-EL,ED=0
op11=>operation: EL=WL,WL=0,ED=C·(K·EM-EU)-EL,WD=WD-ED
io=>inputoutput: E=EU+EL+ED,W=WU+WL+WD
e=>end: 结束

st->cond4
cond4(yes)->op7->op12->io
cond4(no)->op8->cond5
cond5(yes)->op9->io
cond5(no)->cond6
cond6(yes)->op10->io
cond6(no)->op11->io

流域产流计算

产流计算基本思路：流域上每点含水量达到田间持水量时产流，但是各个点的田间持水量和初始水量都是不同的，如何表达呢。模型采用了一种概率方式，利用流域蓄水容量曲线表示。整个流域有初始土壤含水量 W0，其空间分布在蓄水容量曲线上进行表达；PE 的值及其分布同样可以在曲线上表达；那么产流量就可在图上表示出来。计算公式如下：

  • P-E+a < $W_{mm^{'}}$时 $$ R = PE-\int_a^{PE+a}[1-\varphi(W_{m^{'}})]dW_{m^{'}} = P - E - (W_m-W_0)+W_m(1-\frac{a+PE}{W_{mm^{'}}})^{1+b} $$
  • 当 P-E+a $\geqq W_{mm^{'}}$时 $$R=P-E-(W_m-W_0)$$

因此，需要先根据流域初始含水量求出 a 值。这里也是为什么要提前多天计算模型开始计算时间的流域土壤含水量的原因，因为直接在模型起算时间计算时，计算蒸发需要知道 R，而计算 R 需要知道初始土壤含水量，而初始含水量未知，a 无法求，产流量不可知，无法计算蒸发，会形成死循环。而从多天前开始起算，保证前期有降雨使土壤含水量达到过田间持水量，或者长期干旱，流域蓄水量基本为 0，都便于后面的计算。 在经过前面的计算，已知起算时间流域含水量的情况下，可以根据下式求得 a 值： a = ${W_{mm^{'}}}[1-(1-\frac{W_0}{W_m}^{\frac{1}{1+b}})]$

注意计算时，一定有$a>W_0$，$P-E>R$，从流域蓄水容量曲线上很容易看出，如果有不符合上述条件的情况出现，一定要检查是计算机计算的小误差还是因为自己计算或者设置参数有问题导致的。

关于不透水面积，书上使用计算公式Wmm'=(1+B)/(1-IMP) * Wm来处理，那得到的R应该就是已经考虑不透水面积的了，即R是针对整个流域的径流深。不过书上没有给出这块的证明，也没有这部分的推导，暂时我也没有推这部分，觉得稍微有点问题。是不是直接在透水面积上计算，然后最后再将最终计算结果转换到流域面积上更合理一些？目前先按照书本上的做了。

水源划分

水源划分的原因在于后面汇流阶段不同径流成分的水滴其汇流表现不同，汇流时间不一样，因此需要区别处理。 这部分书上解释地不是太详细，这里给出个人的理解，如下所述。

首先流域产流用的标志是流域上一点的土壤含水量超过田间持水量。接下来，对于流域而言，土壤好比地下径流蓄水库，遵循蓄满产流的定义，土壤在达到田间持水量之后，会直接补给地下水，直至该“水库”蓄满，然后蓄满产流，产生地表径流。而该“水库”蓄满的含义就是指产流区域上一点的水量达到或高于自由水蓄水容量，高出的部分就是地表径流。（这块和前面蒸发分几层没有关系，它们互相之间应该是不对应的，因为这里是对自由水讨论的，前面蒸散发是对张力水讨论的） 而自由水蓄水容量在产流面积上分布也是不均匀的，考虑用和产流计算时同样的技巧，即用概率的方式表达，表达的形式和之前产流的方式是一致的————产流面积上的自由水蓄水容量曲线，那么就可以分出地表水和地下水了。 具体的表达如下： $$\frac{FS}{FR}=1-(1-\frac{SMF^{'}}{SMMF})^{EX}$$ 其中，SMF'、SMMF、EX、FS/FR分别表示产流面积上某一点的自由水蓄水容量，产流面积上最大的一点的自由水蓄水容量，流域自由水蓄水容量曲线的指数，产流面积上的各点自由水蓄水容量小于等于SMF'的流域面积 占产流面积的比例。这里需要注意量纲问题，前述容量均为水深量纲。

从上式可以推出(推导过程和产流计算时候的一致): SMF=$\frac{SMMF}{1+EX}$，AU=SMMF$[1-(1-\frac{S_0}{SMF})^{\frac{1}{1+EX}}]$ 其中，SMF表示产流面积上的自由水平均蓄水容量深，S表示自由水在产流面积上的平均蓄水深。

注意$S_0$是上一时段的值，是上一时段计算完之后（计算完之后的意思是每个时段最后还需要重新进行自由水蓄水水库水量平衡计算，在这个计算之后是一个时段计算完之后）的产流面积上的平均自由水深，S是该时段计算开始时候（计算是对一个时段整体计算的，因此就是这个时段的）的产流面积上的平均自由水深，时段间的产流面积是有变化的，因此，如果不能直接用$S_0$直接表达本时段的初始产流面积上平均自由水深，根据工程水文学和水文预报书上的记录，是令S=$S_0\cdot FR_0/FR$，这意思就是上时段计算之后和本时段计算之初（也就是本时段）产流面积上平均自由水深×产流面积相等。然后对本时段后续计算来说，分水源计算开始时的产流面积上的平均自由水深就用S表示。

有了上述表达，就可以推求地面径流的表达式了。 在 FS/FR ~ SMF' 关系图中，S是平均的自由水蓄水水库水深，如果流域产流面积上对应的初始状态是AU，那么接下来当产流量即净雨增加时，在图上增加S对应的SMF'坐标上移，多出来的部分，在曲线左侧的就是地表径流RS了，在右侧的就是地下径流了。 那么现在问题的关键就在于在这样一个关系图中，纵坐标增量如何理解。 纵坐标的增量其实就是产流，即净雨在产流面积上的表达。可以认为是全流域的产流深折算到产流区域上的，即R_产流面积上 = R/FR（注意，FR指的是产流面积占流域总面积的比例），同时，在产流区域上，净雨直接可以由P-E表达，因此也可以得到R/FR=P-E=PE。自然地，在关系图上，增量都是在产流面积上表达的，所以如果用RS表示流域上的地表径流，那么在图上，代表地面径流的曲线左侧的增量部分应该也折算到产流面积上，即曲线左侧阴影面积应该是RS/FR。根据以上分析，可以有如下表达式：

  • 当R/FR+AU<SMMF时，$$RS=FR*[PE-(SMF-S)+SMF(1-\frac{PE+AU}{SMMF})^{1+EX}]$$
  • 当R/FR+AU>=SMMF时，$$RS=FR(PE-(SMF-S))$$

另外，在上述公式中，还存在几个未知量，包括产流面积FR，FR上的平均自由水深S，对应的AU，FR上的自由水平均蓄水容量深SMF，以及FR上最大一点的自由水蓄水容量深SMMF，所以现在问题就转换为如何求这些量。 因为计算产流的时候已经求得了产流的相关数据，因此产流面积易得；SMF与SMMF之间有固定关系表达（见前面公式），所以知道SMMF就可求得SMF；AU可由S求得，所以求出S即可。那么关键就是如何求SMMF和S了。S是可以设置初值；而SMMF则通过另一条曲线来进行分析。 回到流域，定义流域产流面积上最大点自由蓄水容量与产流面积占总面积的比例之间的关系，同样为抛物线型。 $$FR=1-(1-\frac{SMMF}{SMM})^{EX}$$ 则可以求得：SMMF=SMM[1-$(1-FR)^{\frac{1}{EX}}$]，且SMM=SM(1+EX)。 SM流域平均自由水蓄水容量和次幂EX，对一个流域来说，在计算中是定值，因此由流域模型参数来定义，根据资料进行率定。因此各个变量均可求得。

最后，还需要处理地下水流如何划分为壤中流和地下径流的问题。 自由水蓄水线性水库，按比例处理，自由水深×壤中流出流系数×FR，就得到壤中流；自由水深×地下径流出流系数×FR，就得到地下径流。 自由水深是产流面积上的自由水深，就是 FS/FR ~ SMF' 关系图中增量曲线右侧部分。 当PE+AU>=SMMF时，就是SMF，即产流面积上的平均自由水蓄水容量深； 当0<PE+AU<SMMF时，就是△S+S=PE-RS/FR+S。 注意它们都是产流面积上定义的，所以最后还得乘以FR，以转换到流域层面上。 最后计算S来为下一时段的计算准备，即计算时段末产流面积上的平均自由水深，那就是自由水蓄水库中水减去流失的水分，即产流面积上自由水深-（壤中流+地下径流）/FR（注：壤中流和地下径流要折算到产流面积上，所以除以FR）。

汇总以上分析，可以得到：

S设置初值，后面迭代计算。 SM和EX是模型参数。

$$FR=\frac{R}{PE}$$ $$SMM=SM(1+EX)$$ $$SMMF=SMM[1-(1-FR)^{\frac{1}{EX}}]$$ $$SMF=\frac{SMMF}{1+EX}$$ $$AU=SMMF[1-(1-\frac{S}{SMF})^{\frac{1}{1+EX}}]$$

  • 当PE+AU>=SMMF时： $$RS=(PE+S-SMF)FR$$ $$RSS=SMF·KSS·FR$$ $$RG=SMF·KG·FR$$ $$S=SMF-\frac{RSS+RG}{FR}$$
  • 当0<PE+AU<SMMF时： $$RS={PE-SMF+s+SMF[1-\frac{PE+AU}{SMMF}]^{EX+1}}FR$$ $$RSS=(PE-\frac{RS}{FR}+S)KSS·FR$$ $$RG=(PE-\frac{RS}{FR}+S)KG·FR$$ $$S=S+PE-\frac{RS+RSS+RG}{FR}$$
  • PE+AU<=0时： $$RS=RSS=RG=S=0$$

最后更新S的时候，是做自由水蓄水库的水量平衡计算，这里面存在差分计算的误差问题。《水文预报》书上解释是：因为将R作为时段初的入流量进入自由水蓄水库的，而实际上它是时段内均匀进入的，这会造成向前差分的误差。 这个理解是这样的：主要明确下误差的数学原因就行，因为R是时段均匀进入的，所以实际计算的时候，应该是每个时段微分去计算，而做差分的话也应该细致一些，如果直接用时段末和时段初做计算，就是一步向前差分，那这样误差会比较大。 所以《工程水文学》和《水文预报》两本书上都是按5mm去分段净雨，然后逐步去计算的. 事实上，我个人理解，前面的S=$S_0\cdot FR_0/FR$之所以是合理的，实际上就是因为将净雨均分到各个小时段上减少差分的步长来减少误差。因为产流面积的变化是伴随着产流过程发生的，只有当两个时段很接近的时候，才能近似认为产流面积上的自由水量不变。然后接下来的各个推导也可以很好地保持原来的推导形式。而如果直接考虑微分推导的话，RS的计算公式可能很难给出解析表达式，所以这可能也是为什么书上要么干脆没解释，要么解释得不详细的原因。（如果是这么理解的话，实际上可以尽可能缩小计算步长？）

最后特别说明：

  • 以上量纲不是1的变量中： 在流域上定义的有：R,SM,SMM,RS,RSS,RG； 在产流面积上定义的有：S,PE,SMMF,SMF,AU。

  • 模型出流系数计算中的设定： 一般情况下，出流系数是按日模型给定的，在实时洪水预报时，计算时段都小于24h；此外，计算的时候对净雨做了5mm分段，因此模型出流系数要做出对应的变化。 设模型计算所取时段长为△t（h），R为△t内的净雨，则计算步长为时段数乘以量级步长： $$N=\frac{24}{△t}[INT(\frac{R}{5})+1]$$ 注意加1，是因为通常情况下净雨是不会被5整除的，不过为了避免错误，程序中需要进行判断，这里就不赘述了。 计算步长内的壤中流出流系数KSSD和地下水出流系数KGD与其日模型出流系数KSS和KG的关系为： $$KSSD=\frac{1-[1-(KG+KSS)]^{\frac{1}{N}}}{1+\frac{KG}{KSS}}$$ $$KGD=KSSD·\frac{KG}{KSS}$$ 书上解释不太细，我试着反推了一下，由于把壤中流和地下径流系数分开去处理会非常繁琐，且壤中流和地下径流的计算思路是一致的，将每个时段的KSS+KG设为KSSD+KGD，然后可以反推出KSS+KG和KSSD+KGD的关系为$1-(KSSD+KGD)=[1-(KSS+KG)]^{\frac{1}{N}}$。 这个在《水文预报》中有给出。原因的话我觉得从物理上理解可能更容易一些。 从线性水库的角度出发，线性水库出流和蓄量之间是线性关系，N步差分最终蓄量和一步差分最终蓄量相等，可以很容易得到上面的等式。 然后再对壤中流和地下径流分析，无论如何差分，两者之间的相对关系是固定的，即KSS/KG稳定，那么有KSSD/KGD=KSS/KG，联立该式与上式可以求出KGD和KSSD的最终表达式。 关于为什么不直接用每个单独时段的KSS和KG做参数，这是因为每个计算时段的净雨量不同，划分小段后每段的KSS和KG就不同，总是需要这个KSSD、KGD和KSS、KG的关系表达式的。

  • 每个时段计算开始时S的取值： 每个时段开始计算的时候，需要根据下式 $$FR_0·S_0=FR·S$$ 来求计算开始时所需的S值，因为时段间产流面积是不同的，所以这个计算中的FR值，每个时段应该是不一样的。 而用FR=R/PE计算时，因为R值是计算时段的R值，不是按5mm净雨划分之后的，所以这里如果依然用R/PE，可能不太合理，所以给出了这个表达式。 这也是为什么需要给出FR0值得原因。

  • FR和S的初值如何设置? 前辈们都将初值设为一个较小值，比如FR0=0.001,S0=0。但是否合理其实需要分析，如果是久旱无雨，这种情况应是合理的，可是前期土壤蓄水充足，则显然不太合理。较好的方式可能是在计算前期影响雨量的时候就把这部分计算也完成，可以在较早时候令FR0=0.001,S0=0，然后一直计算到当前时段。这可能也是通常模型计算需要一些预热计算的原因之一。这里就先按照FR0=0.001,S0=0进行设置计算了。

  • 计算中单个时段出现S>SMF的情况： 如果出现这种情况说明，S的计算超出了合理范围，需要进行修正。这种情况是由$FR_0·S_0=FR·S$导致的，其计算本身属于估算，因此产生了这种情况。

  • 《水文预报》书上有个公式——$FR_{\frac{\Delta t}{N}}=1-(1-FR)^{\frac{1}{N}}$： 给出这个公式的原因据书本上介绍是由于产流面积FR是随着自由水蓄水容量的变化而变化的，当计算时段长改变以后，它也要做相应的改变。 这个意思就是按PE的净雨和分段之后的5mm净雨，产流面积显然是不同的。所以FR值会有所变化。因为5mm对应的产流值不知道，所以没法直接用5mm算，除非从前面产流计算的时候算5mm就能停下来，可是前面的计算都是连续的过程，无法对应处理。所以才要寻找两个不同净雨之间的FR的对应关系。 这个公式从形式上看，和前面地下出流系数分段推求相似。也是做了线性化的处理，公式的意思是不产流面积的比例每小段的累乘等价于总的不产流面积比例。但是首先书上并没有给出很严谨的推论，只是直接给出了结论，所以这个地方有待解释！因为我也没解释不了，加上实践上看，应该是可行的，所以就先使用这个公式了。

  • PE<0的情况： 书上没有给出这种情况的处理方法。PE<0时，流域没有产流面积，但是从地下线性水库概化角度，从物理层面分析，如果前一个时段是有产流的，那么其自由水蓄水量是没有完全泄流的，所以只要地下水库有自由水蓄水，那么自由水应能继续出流，并仍按出流系数进行。如果没有产流就令计算时段初的S=0，那么就相当于凭空消失了上一时段末计算留下的自由水蓄量，这是不符合水量平衡原则的，这一原则还是应该保证的。 但是如果这样，这里还会和前面产流计算不匹配，即前面计算产流时，当P-E<0时，产流量为0，而这里又有了产流量，前后似乎矛盾了。 这应该怎么解释呢？书本上没有讲。我认为目前暂时比较妥善的解决方法是主要依据水量平衡原则，并仍利用已有的公式，$FR_0·S_0=FR·S$，如果上一时段末的S已知，假设没有产流，地下水库没有增加，那么本时段初计算时，S应该保持不变，即$FR_0·S_0=FR_{时段初}·S_{时段初}$，所以$FR_{时段初}=FR_0$，然后剩余计算以地表径流为0，地下径流线性出流为原则进行计算，这样可以处理这一问题。一直计算到模拟时段结束。

  • 不透水面积的影响： 不透水面积对应的流域点地下水不能概化为地下自由水蓄水水库，直接产生地表径流，所以上述的分析都是在除不透水面积之外的流域上开展的，对水深量纲没有影响，不过后面算水量的时候需要留意，可以在这里先转换到全流域的水深。 如果前面产流计算的时候，用了书上的公式Wmm'=(1+B)/(1-IMP) * Wm，那得到的R是不是已经考虑不透水面积的了？那到这边分水源的时候有没有影响？ 个人认为如果IMP真是的不透水面积/全流域面积，那应该是有影响的，FS/FR ~ SMF'那张图会变，因为产流面积上有一部分不透水，那么自由水蓄水能力小于等于SMF'的流域面积应该包含一部分不透水面积，在曲线上应该是不能从原点开始的了，而是从横轴上某一点开始。 但是这里我认为应该没有真的把IMP考虑进产流模型中，只是在WMM'的处理上用了一下，实际上后面的产流模型根本就还是在透水面积上做的，不然个人认为直接套用IMP=0时公式是不太合理的。 所以计算还是在透水面积上开展的，这样也比较合理，因为如果前面产流计算是在透水面积上开展的，公式的推导都是没问题的，然后利用1-imp把结果折算到了整个流域面积上，虽然我个人认为“应该先计算，最后再折算”应该更合理一些，但是目前书本上都是这样的，并且一般产流计算问题都不大，所以产流部分关于不透水面积的计算暂且这样。 所以分水源计算开始前，如果产流值是整个流域的，需要先折算到透水面积上，然后在分水源计算完毕之后，有必要再把计算结果折算到整个流域面积上。 PS：总之前面这最后四个补充的说明总觉得里面有些瑕疵，我是无能为力了，就先这样了。如果有高人指点，请赐教！！！

汇流计算

三水源汇流计算包括地表径流、壤中流和地下径流的坡地汇流，每个单元上有单元的河网汇流，单元面积以下有河道汇流。

地面径流汇流

可以采用单位线，也可以采用线性水库。二者应该是一套理论下来的，根据实际情况，现实中基本上是用单位线做汇流计算的，所以这里就以单位线为主进行分析。虽然实际生产中貌似各地都有做好的单位线，但因为气候变化和人类活动影响，也不知道还靠不靠谱，总之世界在变，还是从理论上出发，来分析，这样不但严谨，后面还可以结合最新的进展做研究。

如果学过信号与系统之类的课程，一定不会对卷积、传递函数这些概念感到陌生。为什么说到这些和水文好像有点远的东西，因为实际上单位线的思路是比较“系统”科学的。 对一个流域单元来说，在汇流阶段，它体现出一种对净雨过程的推移和坦化作用，即观察出口断面流量，其出口流量过程相比于净雨过程，总是洪峰迟，洪峰小。那么如何研究这一过程呢，用系统科学思维，可以把单元流域看成一个接受净雨做输入而输出流量的系统，那么汇流计算的问题就转换为了对一个输入输出系统的研究。 对于动力系统的研究方法应该挺多的，比如微积分里学的微分方程、相位图等等，在经典的信号与系统教材里讲的是如何分析系线性时不变系统的特性，毕竟线性时不变是最简单的，而且非线性的也可以用线性逼近。而在分析流域的汇流特性时，经典的单位线就是基于系统是线性时不变系统假设的。 对一个线性时不变系统，分析的基础是单位冲激响应，可以参考科普视频1和视频2。

单位冲激响应就是瞬时单位线。单位矩形函数是时段净雨的形式了，单位净雨输入到系统的响应就是时段单位线。有了时段单位线，利用卷积的离散形式，和净雨过程做运算即可得到单元出口流量了，每个时段都这么算，最后将计算结果对应时段累加即可。公式如下所示： $$Q_i=\sum_{j=1}^m \frac{h_j}{10}q_{i-j+1} \begin{cases} i= 1,2,...,l\ j= 1,2,...,m\ i-j+1=1,2,...,n \end{cases} $$

关于对卷积的理解可参考说明。 我觉得上式的理解难点有两处，一是卷积运算里函数的取反，即翻转，以及滑动叠加，这两个运算的结合； 二是卷积运算是一个序列运算，每个时刻的运算结果不仅考虑了当前时刻的输入响应，还跟过去所有时刻的输入响应有关，而跟哪些过去有关，是利用卷操作来指定的。 放到净雨和单位线的卷积，就是说一个净雨它的作用是在单位线每个时段上都留下一个输出，从第一个时段净雨输入到单位线中，它不仅在第一个时段留下了输出，它给后续所有时段都留下了输出，那么到最后一个时段，它也给它后续多个时段留下了输出。所以才有了上述表达式。

汇流研究的其实主要是两类问题，一是刚刚提到的已知净雨过程和流域单位线推求出口断面流量过程的“预测”问题；另一个镜像的问题就是已知净雨过程和相应地出口断面流量，如何反演单位线的“识别”问题。所以问题的重点其实在于如何要从历史数据中反演出单位线。

推求单位线有从不同角度，按不同理论的很多种方法。

比较经典的是分析法求时段单位线。

  1. 选择若干场时空分布均匀，历时较短的降雨形成的单峰洪水来分析；
  2. 点绘流量过程线后分割基流及前期洪水退水，计算本次洪水的直接径流深（具体过程见下一节分析）；
  3. 对该场次降雨进行产流计算的结果作为净雨过程；
  4. 根据净雨郭晨估计对应的地面径流流量过程线，利用以下公式推算单位线： $$q_i=\frac{Q_i-\sum_{j=2}^m{\frac{h_j}{10}q_{i-j+1}}}{\frac{h_1}{10}} \begin{cases} i= 1,2,...,n\ j= 1,2,...,m \end{cases}$$ n为单位线的时段数，n=l-m-1；
  5. 修正单位线：检查总量是不是10mm净雨量，是否为单峰，过程是否平滑，如果不是需要修正。

上式中注意分母是$h_1$，这个结果在卷积运算方程组列出之后可以很清楚地看到。

用最小二乘法等优化求解方法也可以推算单位线。

各时段净雨和地面径流可以写作向量q和Q，单位线作用于各个时段净雨，可以写作：Q=hq。式中， $$h=\left[ \begin{matrix} h_1 & & \ h_2 & . & h_1 \ ... & . & h_2 \ h_m & . & ... \ & & h_m \ \end{matrix} \right], q=\left[ \begin{matrix} q_1 \ q_2 \ ... \ q_n \end{matrix} \right], Q=\left[ \begin{matrix} Q_1 \ Q_2 \ ... \ Q_L \end{matrix} \right] $$ 其中，n=l-m+1。h是一个l*n的m对角阵，q是n维向量，Q是l维向量。方程个数多于变量个数，无解，这时使用最小二乘法求优化解是常见方法，求解$h^Thq=h^TQ$即可。 这里涉及的线性代数可参考视频1和视频2。

最小二乘法求解时，有时也会出现锯齿状震荡或负值，因此需要加一些约束条件求解更好。

以上是一种最常见的求时段单位线的方式，个人理解是一种直接求离散结果的方法，而实际上单位线的形式、识别方法是很多的，我认为比较简单的是通过概念性元素模拟来理解分析。从流域汇流宏观考量，模拟推移和坦化的作用。比较常用的有以下几种：

  • 只传播，不坦化，没有变形的，概化为线性渠道，用来模拟推移作用；
  • 线性水库，表达坦化作用；
  • 线性面积-时间曲线，体现汇流的分布式入流状态。

它们的串并联组合，可以构成流域汇流的基本结构，以此来组成流域单位线的基本表达式，然后再根据实测数据来推求表达式的相关参数，以“识别”出流域的单位线。

比如《工程水文学》和《水文预报》上都重点讲解的Nash单位线，其把流域概化为了n个串联的线性水库。然后推求出数学表达式。 为了根据实测数据推求单位线中的参数n和k，通过统计方式进行处理，分析单位线的一阶原点矩和二阶中心矩，进而可以推出n和k值，也可以用优选的方法来计算。

本项目给出的代码就先以上面谈及的单位线为例。

需要注意，在计算过程中由于需要利用实测数据进行参数识别，而前面计算的产流结果是分水源的，所以需要对实测数据进行分水源分析。为了分析地面径流，需要把地下径流从实测数据中割除出去。《工程水文学》上介绍的方法有较为简便的斜线分割法，还有经验公式法。这部分的编程处理略微有些麻烦，具体的数学表达还需考虑。

应用单位线时还有几类实际问题需要补充说明：

  1. 单位时段不同的单位线之间的转换
  2. 实际中时变和非线性特性的影响
  3. 无水文资料地区

应用单位线计算汇流时，常常因净雨量时段长和单位线的时段长不一致而引起误差。因此需要不同各时段长的单位线转换。

如何理解？ 从S曲线开始说起。前面提到流域时段单位线实际上是单位矩形函数的系统响应。 单位矩形函数可以表达为两个阶跃函数的差（那么不同时段长对应的就是阶跃函数的相对位置不同），阶跃函数求导就是冲激函数，所以容易推求冲激响应积分就是阶跃响应，那么阶跃响应的差就是时段单位线。 定义流域S曲线：$S(t)=\int_0^t u(0,t)$，即单位阶跃响应，则时段单位线$u(\triangle t,t)=S(t)-S(t-\triangle t)$。 流域S曲线的定义就是流域上分布均匀，且一直维持在1个单位强度的净雨所形成的流域出口断面流量过程。

注意实际中，都是离散条件下的计算： $$S(t)=\sum_{j=0}^k q_j(\triangle t,t)$$ S(t)表示第k时段末S曲线的纵坐标，$m^3/s$；$q_j$表示时段为$\triangle t$的单位线第j时段末的纵坐标，$m^3/s$；$\triangle t$是单位线时段，h。其反映连续多时段单位净雨深所形成的的流量过程线。 上式其实是一个卷积运算，只不过因为是离散单位线和离散阶跃函数卷积，阶跃值为1，每个时段都是一个1，所以形式上就只留下了单位线，但是阶跃函数仍然控制着叠加的范围，其实书上的公式那个k是无穷大的，但是阶跃函数的每个离散值都相当于一个冲激函数的离散值，它只有1个值，对应到响应上也只能有单位线时段个输出，所以使得为1的阶跃函数和不为0的单位线值卷积运算的数目到后面都是固定的，这也是为什么S曲线后面变为直线的原因，所以才有了这个k是单位线时段个数的说明。 从离散的角度分析，也很容易理解为什么阶跃响应是冲激响应的积分了。所以书上这点说的不是太清楚。

了解了S曲线及其离散形式之后，接下来分析为什么单位线时段长不一致会引起误差。 由于单位线的定义是单位时段10mm净雨对应的出口断面流量，假设单位线时段变长了，按照单位线定义，是单位时段10mm净雨对应的流量，以前单位时段是$\triangle t_0$，现在是△t（△t>$\triangle t_0$），这时候如果还是直接使用单位线，意思就是认为$\triangle t_0$内10mm净雨产生的流量，在每个△t上的大小还是原来单位线的大小，举例，比如原来1个小时10mm净雨产生的流量大小共三个时段，2-3-1，现在2个小时共10mm净雨产生的流量大小不变，即2个小时10mm净雨对应的流量是2（2h）-3（2h）-1（2h），注意和1h 5mm 2-3-1有区别，不过不管怎样，都显然是不合理的。总之，时段变化的话，原来时段的单位线肯定是不能用的。

那能不能进行时段间的单位线转换？ 可以的。关键是把握不变的东西，单位线的时段在变，形式在变，但是不变的是流域本身的性质，即单位冲激响应是不变的，离散状态下，就是说单位强度的净雨输入的响应是不变的。但是仅凭这一点还推不出不同时段间的单位线之间的关系，还有线性时不变系统的线性，即n个单位强度的净雨输入的响应是单位强度的净雨输入的响应的n倍。所以利用S曲线，可以求得强度为10/$\triangle t_0$的净雨形成的$\triangle t$时段的流量过程线，这样和$\triangle t$单位线形成的流量过程线能在时段上对应起来，然后他们都是单位强度净雨形成的流量过程线的倍数，倍数分别是10/$\triangle t_0$和10/$\triangle t$，两式联立，自然求得了$\triangle t$时段对应的单位线了。公式如下： $$q(\triangle t,t)=\frac{\triangle t_0}{\triangle t}[S_0(t)-S_0(t-\triangle t)]$$ 书上用的是S(t)-S(t-△t)，虽然下面注明了是时段为$\triangle t_0$的S曲线，但是个人认为容易误导，所以还是用$S_0$表示较好。

对于瞬时单位线，由于实际中运算是离散的，所以一般都是先把瞬时单位线转换为S(t)曲线，然后用S曲线推求无因次的单位线，最后再转换为单位线。

由于时变非线性因素的影响，在建立单位线时，需要进一步处理。

对于时段单位线：

  • 如果各次洪水求得的单位线差别不大，可以取平均。
  • 如果差别较大，需要分析影响单位线变化的主要因素，如降雨强度、暴雨中心位置等，分类求单位线。

对于瞬时单位线：

常用净雨强度反映非线性的影响。 有瞬时单位线的滞时与时段净雨量之间的关系。$T_{L,0}=C_t r^{-\alpha}=nk$。利用实测资料可以建立$T_{L,0}$和r之间的经验关系，据此可以推求不同净雨强度的瞬时单位线，然后在一场降雨中，多时段净雨量各自采用不同的瞬时单位线计算出流量过程，再线性叠加，即为流域出流量。

对于无水文资料地区，建立单位线要素与自然地理特征之间的相关关系，有助于汇流计算。

地下径流汇流

根据地下水流运动的基本微分方程可以导出地下水径流的流域汇流模型，但是在应用这类模型时需要地下水位，有关的水文地质和土壤特性等数据，较难实现。新安江模型采用的是以水量平衡方程和线性水库的蓄泄关系为基础的水文学方法，线性水库演算法。 实际上也是一种最简单的单位线方法，即一个线性水库的单位线方法，可以先求解微分方程计算瞬时单位线，然后再推求时段单位线，也可以直接计算方程的差分形式，直接计算结果。

前面已经提过线性水库了，这里再写出其表达，主要是两个方程，一个是朴素的水量平衡方程，另一个就是线性水库的特点： $$\begin{cases} I_g-Q_g=\frac{dW_g}{dt}\ W_g=K_gQ_g \end{cases}$$ 差分计算书本上采用的是中心差分形式（如果对差分有些疑问，可以参考这篇博客，如果想要再看看差分是怎么和微分联系的，可以简单复习下微积分），即$\frac{dW_g}{dt}=\frac{W_g(t+0.5\triangle t)-W_g(t-0.5\triangle t)}{\triangle t}$，令$W_g(t-0.5\triangle t)=W_{g1},W_g(t+0.5\triangle t)=W_{g2}$，那么上述两式变为: $$\begin{cases} \bar{I_g}-\bar{Q_g}=\frac{W_{g2}-W_{g1}}{\triangle t}\ W_{g1}=K_gQ_{g1}\ W_{g2}=K_gQ_{g2} \end{cases}$$ 联立上述各式，容易求得: $$Q_{g2}=\frac{\triangle t}{K_g+0.5\triangle t}\bar{I_g}+\frac{K_g-0.5\triangle t}{K_g+0.5\triangle t}Q_{g1}$$ 其中，$\bar{I_g}$指的是地下净雨对应的流量，则$\bar{I_g}=\frac{1000\cdot RG\cdot F}{3600\triangle t}$，RG表示地下净雨深，mm；F表示流域面积，$km^2$。令$KKG=\frac{K_g-0.5\triangle t}{K_g+0.5\triangle t}$，则上式课简化为$$Q_{g2}=Q_{g1}·KKG+RG·(1-KKG)·U$$

书中新安江模型里，令壤中流和地下径流机制相同，因此总结上述内容。壤中流RSS和地下径流RG分别流入壤中流水库和地下径流水库，经过蓄水库的调蓄（记两水库的消退系数为KKSS、KKG），成为地下水对河网的总入流TRSS、TRG。计算公式如下： $$TRSS(t)=TRSS(t-1)·KKSS+RSS(t)·(1-KKSS)·U$$ $$TRG(t)=TRG(t-1)·KKG+RG(t)·(1-KKG)·U$$ 其中，U是单位转换系数，因为前面的计算量纲都是水深，这里转换为流量.U=$\frac{F}{3.6\Delta t}$，F是流域面积，单位是$km^2$。 从公式中，可以看出为什么取KKSS和KKG两个系数，并取名为消退系数了。

模型的消退系数也是按日模型给定的，所以进行实时洪水预报时，计算时段小于24h时，需要作相应地变化。 设模型计算所取时段长为△t（h），R为△t内的净雨，则： $M=\frac{24}{△t}$， 计算步长内的壤中流蓄水库消退系数KKSSD和地下水蓄水库的消退系数KKGD分别与其相应地日模型消退系数KKSS和KKG的关系为：$KKSSD=KKSS^{\frac{1}{M}}$，$KKGD=KKG^{\frac{1}{M}}$。 这个比之前分水源计算的那个简单，直接就是累积就行了，容易得到。

计算完毕之后，汇总：TR(t)=TRS(t)+TRSS(t)+TRG(t)

单元面积河网汇流

一般采用线性水库或者滞后演算法进行。滞后演算法等价于线性渠和线性水库的串联，《水文预报》书上的公式是： $Q_t=CR\cdot Q_{t-1}+(1-CR)\cdot QT_{t-L}$ 我个人反推了下，这个公式等价于先一个线性渠然后串联一个线性水库。

单元面积以下河道汇流

河道演算。新安江模型一般采用马斯京根法进行河道演算。河道流量演算是以圣维南方程组为理论基础，利用上断面流量过程演算下断面流量过程。马斯京根法是对圣维南方程组近似得到的水量平衡方程和槽蓄方程进行求解的一种方法。现在直接求解圣维南方程的数值解法也比较多，也较常应用。关于圣维南方程的理解，可以参考这个，和前面一样，如果想复习下相关的微积分基础知识，推荐参考这个，关于纳维尔斯托克斯方程和圣维南方程之间的关系，简单查到了这个。因为这里要用水文方法求解，并且我流体力学、水动力学都基本空白了，所以暂时就不赘述了。

水文法简化了圣维南方程，得到了水量平衡方程和槽蓄方程。马斯京根法建立了马斯京根槽蓄曲线方程，将其与水量平衡方程联立，可得到马斯京根流量演算方程： $Q_{下,2}=C_0Q_{上,2}+C_1Q_{上,1}+C_2Q_{下,1}$ C0，C1，C2都是KE，XE的函数。

实测径流分析

一次洪水流量过程除包括本次洪水所形成的地面径流、表层流径流和地下径流外，往往还包括前期洪水尚未退完的部分水量及非本次降雨补给的深层地下径流，因此如果要对模型计算结果进行验证，需要对观测径流量进行分析计算。其次，还需要对洪水流量过程中的不同水源成分进行划分，才能在计算单位线时，有实测地表径流数据。

退水曲线分析

一个常见的流域出口断面的实测流量过程是前后洪水首尾相接的。因为流域出口断面流量是由不同水源的径流成分组成的，它们汇流过程不同，尤其是地下水，退水延滞性强，因此复式洪水过程很常见。 为了分割场次洪水过程以计算场次洪水总径流深或分割地下水，需要确定流域退水曲线，所以要分析流域的退水特性和规律。分析方法有多种，为了编程计算的便利，采用退水指数方程的方式，其余方式可参考书本——《水文预报》，在此不赘述。

退水指数方程： $$Q_t=Q_0e^{-\frac{t}{K}}$$ 式中，$Q_t$，$Q_0$分别为t时刻流量和起始退水流量，K是常数。 $$W_t=KQ_t$$ $W_t$为t时刻的蓄量。此式表明，当泄流量恒定时，K是泄完蓄水量所需时间，相当于流域平均汇流时间。

给出退水指数方程的差分形式，可以得到： $$Q_{t+1}=Q_te^{-\frac{1}{K}}$$ 令$C_g=e^{-\frac{1}{K}}$，则$K=-\frac{1}{lnC_g}$，$C_g$为常系数，反映退水速率，称为流量消退系数。

消退系数可由由计算时段始末的两个实测退水流量确定：$C_g=\frac{Q_{t+1}}{Q_t}$，由于实测流量的观测误差等，选择多组观测值进行计算较好，根据最小二乘法，书上给出如下公式： $$\hat{C_g}=\frac{\sum_{i=1}^nQ_{i,1}Q_{i,2}}{\sum_{i=1}^nQ_{i,1}^2}$$ 选择的样本代表性越佳越好。 实际编程计算中，采用矩阵运算很容易，也很容易推出上述公式，详见代码。

流量过程分割

在制作了流域退水曲线之后，可用于分割场次洪水过程线，计算次洪总径流深，分割地表径流和地下径流。 这部分以往常常手工计算，因此书本上没有给出十分适于编程的方式，结合《水文预报》和《工程水文学》两本书，给出相对来说便于编程计算的方法，表述如下（参考《工程水文学》图4-6）：

因为基流不是本次降雨形成的，首先从流量过程线将其分割出去（所以地下径流一般指的是本次降雨形成的浅层地下径流，深层地下径流比较稳定，流量也小），一般取历年最枯流量的平均值或本年汛前最枯流量用水平线分割，令为直线ED；

然后根据洪水流量过程找到场次洪水起涨点A；

根据A和退水曲线计算出其与ED交点F；

用退水曲线段反算A至刚大于B停止；

将B-A-F曲线与B点开始的场次洪水过程线做差值；

倒序检查极小值点，在B-C'之间选择一个靠近C'的极小值带你作为分离点；

然后按退水曲线计算退水过程即可；

注意在离散状态下计算地表径流值；

排除掉上次洪水的退水过程，也可以得到地下径流值。

模型参数率定

模型参数概念分析

新安江模型是概念性模型，在进行计算与率定之前，需要对模型参数的物理过程归纳，然后初定参数范围。关于水文模型参数的率定需要人工经验，很多文献都有提及，一直以来都是想要去尝试改进的点，能力所限就不展开了。

  • 蒸散发能力折算系数K K是影响产流量计算的重要且敏感的参数，产流计算中，其控制着水量平衡。蒸散发模型中，普遍应用的一个物理量是流域蒸散发能力EP。流域蒸散发能力是指供水充分情况下的流域日蒸散发量。它决定于气象因素、土壤特性及植被状况等，基于物理过程建立的理论需要资料较多，所以目前广泛采用的是借助水文蒸发皿的观测值来推求流域蒸散发能力并采用蒸散发能力折算系数做转换。一般KC<1.0
  • 不透水面积比IMP 不透水面积比可用地图进行分析，也可用历史旱期小洪水资料来分析，因为干旱期降雨形成小洪水可以认为完全是由不透水面积产生的，求出该场次洪水的径流系数，其值即为IMP。天然流域IMP=0.01~0.02. 人类活动的加剧使其有明显增大的趋势。
  • 流域蓄水容量——面积分布曲线指数B 流域蓄水容量面积分布曲线指数反映划分单元流域张力水蓄水分布的不均匀程度，一般面积越大，取值越大，几平方公里的B=0.1左右，几百到千平方公里的B=0.2~0.3；几千平方公里的B=0.4左右。
  • 流域平均张力水容量WM 流域平均张力水容量反映了流域干旱程度。分为上下深三层。可以根据前期久旱，久旱之后发生的，降雨特别大的，使全流域产流的历史降雨径流资料来确定，久旱W=0，流域蓄满W=WM，利用水量平衡，根据P,E,R即可计算出WM。可以根据经验判断，南方湿润地区WM约为120-150mm，半湿润地区WM约为150-200mm。WUM约可取20mm，WLM可取60-90mm。
  • 深层蒸散发扩散系数C C值主要取决于流域内深根植物的覆盖面积。目前根据经验，在南方多林地区C=0.150.20；在北方半湿润地区C=0.090.15。
  • 自由水蓄水容量SM SM反映表土蓄水能力，其值受降雨资料时段均化的影响，当以日为时段长时，在土层薄的山区，其值为10mm，或更小一些；在土深林茂透水性强的流域，其值可取50mm或更大一些，一般流域在10~20mm。当计算时段减小时，其值要加大。
  • 自由水蓄水容量面积分区曲线指数EX EX值反映流域自由水蓄水分布的不均匀程度，答题发拧了饱和坡面流产流面积的发展。由于饱和坡面流由坡脚向坡上发展时，产流面积的增加逐渐减慢，故认为EX应大于1.0，一般取值1.0~1.5.
  • 自由水蓄水库对地下水和壤中流的日出流系数KSS、KG KSS反映表层土的渗透性，KG反映了基岩和深层土壤的渗透性。在模型中，这两个出流系数是并联的，其和KSS+KG代表出流的快慢，区别代表壤中流和地下径流的比。1-（KG+KSS）为消退系数，它表明了直接径流的退水快慢。一般，中等流域退水历时在3d左右，取KSS+KG=0.7；若退水历时为2d，取0.8；若退水历时远大于3d，表示深层壤中流在作用，应考虑用壤中流消退系数来处理。
  • 壤中流消退系数KKSS 若无壤中流，则KKSS=0；若壤中流丰富，则KKSS=0.9，相当于汇流时间为10d。
  • 地下水消退系数KKG 可以根据枯季地下径流的退水规律来推求，KKG=0.95-0.998，大致相当于20-500天。
  • 河网单位线 单位线取决于河网的地貌特征，按前述方法，根据历史数据推求。
  • 河网蓄水消退系数CR，滞后时间L 取决于河网的地貌条件，通过模型模拟检验。
  • 马斯京根法参数KE,XE 取决于河道特征和水力特性，可以根据河道的水力特性采用水力学方法或者水文学方法推求。

参数率定算法

由于数据所限，实际应用中，多用系统识别和人工经验结合的方法来率定参数。我是没什么经验的，所以这部分大部分直接从《水文预报》中截取了，在这里打出来是给自己过一遍用的，可以直接参考原书。

参数敏感性分析

敏感性分析是一种对不确定性的分析技术，分析不确定源可以帮助更好地理解水文模型。对不敏感的参数，可以粗略一些货根据一般经验固定下来，不参加优选。新安江模型参数课分蒸散发计算、产流计算、分水源计算和汇流计算四个层次，分别对应K、WUM、WLM、C；WM、B、IMP；SM、EX、KG、KSS；KKSS、KKG、CR、L、KE、XE。

参数相关性分析

模型中如果存在相关度较高的参数，解就不稳定。新安江模型有些参数之间不独立性既存在于层次之内，也存在于层次之间。这种相关性使得异参同效现象可能较显著，其解就不稳定，也不唯一，模型率定难选最优解。因此采用一些结构性的参数约束，可以减少这种现象，使率定优化更容易。

实际应用中发现，新安江模型有些参数之间的不独立性既存在于层次之内，也存在于层次之间。

  • 层次之间参数的相关性分析。

    • 第一、二层次之间参数。 当第二层次中参数B有变化时，对总径流量R的计算结果会产生一定的影响，因此就会影响总的水量平衡，也就影响了第一层次参数的调试结果。因为B只在局部蓄满产流时起作用，当全流域蓄满时就没有作用了，因此，B可根据次洪的降雨径流关系求出，因此与第三层次参数无关。随着参数B值的增大， 计算的年径流有增大的趋势， 但影响很小，对水源划分的影响更小。 WM不影响蒸散发计算， 因此与第一层次参数无关。但WM与B有关，因此，对产流产生一些间接的影响。 天然流域IMP很小，影响不大，但都市化地区较大，对产流有一定影响,也与其他参数无关。
    • 第二、三层次之间参数。 由于产流计算采用蓄满产流，在分水源以前，总径流量R已经计算好，所以第三层次参数完全不影响第二层次参数。
    • 第三、四层次之间参数。 在分水源计算结束后，所求得的是河网总入流。汇流计算只处理河网汇流问题，与水源划分无关。所以，第三、四层次之间参数在性质上是完全独立的。但在优化参数时，都是根据流域出口断面的流量过程线，因此在定量上有一定的相关性。但流量过程线与这两个层次间参数的关系，可以通过流量过程线的分段处理来解决。在高水部分流量基本上是由地面径流和壤中流组成，主要调整参数SM、EX、KG+KSS、UH(CR、L)；退水段尾部流量基本上是由壤中流和地下径流组成，主要调整参数KG/KSS、KKSS、KKG; 低水部分小流量基本上是由地下水组成，主要调整参数KG/KSS、KKG。因此，若分段采用不同的目标函数，可以克服某些参数之间相互不独立的问题。

  • 同一层次中参数的相关性分析。

    • 第一层次中参数。 在第一层次中，若加大WUM、WLM、C的值，计算的蒸散发量E值就会增加。因此，为了控制水量平衡，调试时就会减小K的值。由于WUM与WLM都有一定的变化范围，所以这种影响是有限的。深层蒸散发扩散系数C值只对干旱季节起作用，由于湿润地区很少用到深层蒸发计算，所以一般情况下它是不敏感的，但对半湿润地区，它则是重要的。
    • 第二层次中参数。 在第二层次中，如果流域蓄水容量一面积分布曲线保持不变，则WM值愈大，B值愈小，两者并不是独立的，它们共同确定了流域蓄水容量面积分布曲线。WM不影响蒸散发计算。WM有一个约束条件，不能出现负值，若出现负值，WM要酌情加大。
    • 第三层次中参数。 新安江模型中，第三层次参数是敏感和重要的，参数相互之间的关系也相对比较复杂。SM与EX之间是不独立的，它们共同确定了自由水蓄水容量面积分布曲线。存在的问题和WM与B之间存在问题大致相同。但WM与B之间关系可以根据降雨径流相关图推求出，而SM与EX之间的关系则没有类似的途径可以解决，只能采用优化检验的方法分析。 方法是先根据经验调试好全部参数,然后将其他参数固定不变，只调试SM与EX, 确定最优解的范围。EX值大体上反映了流域自由水分布的不均匀程度。据经验，其值的最优范围在1.0~1.5，变幅不大，可令EX=1.5。 SM与KG+KSS间存在相关关系。若KG+KSS和EX固定，SM就决定了地面径流的多少，KG/KSS就决定了壤中流和地下径流的比例。当SM增大时，RS要减小；若同时减小KG+KSS，则RS可以保持不变；若KG/KSS也不变，则RSS与RG也保持不变。这3个参数之间是不独立的，这种不独立性是由线性水库结构所造成的，经验上，可以采用结构性约束KG+KSS=O.7来解决。这等于把参数减为两个，SM与KG的独立性增强了，可优化得出惟一解，只有这种解才存在各流域之间的可比性，找出区域性规律。 分析方法是先根据经验调试好全部参数，然后将其他参数固定不变，只调试SM与KG/KSS, 确定最优解的范围。由上述分析可见，分水源层次参数的独立性问题特别复杂，必须要加结构性约束才能解决，如不加约束，则各种模型的分水源参数的最优解可能有多种组合。
    • 第四层次中参数。 在汇流层次中，KKSS的作用是弥补KG+KSS=O.7的不足。它取决于退水段退水流量的快慢，与其他因素无关，因此相对比较独立，它对整个过程的影响远不及SM与KG/KSS明显。 KKG取决千地下径流的退水快慢，也相对比较独立，用枯季径流的资料很容易确定。 UH(或L、CR) 决定于流量过程线的中高水部分，因此与第三层次参数之间是比较独立的。L与KS的功能不同，前者处理平移，后者处理坦化，相互间是很独立的。UH(或L、CR) 与KE、XE之间具有相关性。解决的方法是：若单元面积汇流快一些，则河网汇流就可以慢一些，反之亦然，使它们相互间有一定的补偿作用。
结果评定

率定之前，需要明确洪水预报结果是如何评定的。 洪水预报的对象是控制断面的洪水要素，包括洪峰流量（水位）、峰现时间、洪量（径流量）和洪水过程等，对应其精度评定项目。

误差指标一般采用绝对误差、相对误差、确定性系数三种。 水文要素的预报值减去实测值为预报误差，预报误差的绝对值为绝对误差，多个绝对误差值得平均值表示多次预报的平均误差水平； 相对误差即预报误差除以实测值，多个相对误差绝对值的平均值为多次预报的平均相对误差； 洪水预报过程与实测过程之间的吻合程度可用确定性系数DC作为指标，其数学表达式为： $$DC=1-\frac{\Sigma^n_{i=1}[y_{ci}-y_{oi}]^2}{\Sigma^n_{i=1}[y_{oi}-\bar{y_{o}}]^2}$$ 式中，DC一般取两位小数；$y_c$为预报值，$m^3/s$；$y_{o}$为实测值，$m^3/s$；$\bar{y_{o}}$为实测值均值，$m^3/s$；n为资料系列长度。DC就是著名的NSE。

根据我国SL250-200《水文情报预报规范》，调试时通常以年径流绝对误差最小及确定性系数DC最大为目标函数。本项目代码部分就先以此为示例进行参数率定。不过不影响接下来的说明，另外，后期有需求再补充其他寻优目标。

许可误差是依据预报精度的使用要求和实际预报技术水平等综合确定的误差允许范围。由于洪水预报方法和预报要素的不同，对许可误差规定亦不同。包括洪峰预报许可误差，峰现时间预报许可误差，径流深预报许可误差和过程预报许可误差。这些根据书本上的说明，还是挺繁琐的。

预报项目精度评定规定:

定义合格预报为合格预报次数与预报总次数之比，表示多次预报总体的精度水平。一次预报的误差小于许可误差时称为合格预报。 $$QR=\frac{n}{m}\times 100%$$

预报项目精度按合格率QR或确定性系数DC的大小分为甲乙丙三个等级。

当进行实时作业预报时，还需要对洪峰预报时效等级进行分析。用洪峰预报时效性系数描述: $$CET=\frac{EPF}{TPF}$$ 式中，CET为时效性系数；EPF为有效预见期，指发布预报时间至本站洪峰出现的时距；TPF为理论预见期，指主要降雨停止或预报依据要素出现至本站洪峰出现的时距。也分为甲乙丙三级。

率定算法

形式上包括人机交互率定（假设参数，机器运算，比较实测值与计算值，人工经验调整参数，直到最优）和优化算法自动率定。

新安江模型参数率定一般分为“日模型”和“次洪模型”两大部分。模型可靠性取决于模型研制中使用的水文资料的质量和代表性。洪水预报方案要求使用样本数量不少于10年的水文气象资料，应包括大中小水各种代表性的年份，并保证有足够场次洪水资料，湿润地区不少于50次，干旱地区不少于25次。

对于日模型： 对一般小流域，可以不划分单元面积，用流域平均面雨量计算，大流域通常划分单元面积。 通过日模型调试可以确定第一、二层次的参数，第三层次的部分参数如KG、KSS、EX，以及第四层次部分汇流参数如KKG等。 通过日模型的调试可为次洪模型调试提供初始状态变量，如流域张力水蓄水量W、WU、WL、WD，自由水蓄水量S及退水流量等。 比较多年总径流量R，这是最基本的水量平衡校核。如有误差，首先调整K值，K是影响蒸发计算最敏感的参数，它控制着水量平衡。若发现季节性变化对流域蒸散发的影响很大，则应考虑按不同季节分别率定K值。 比较每年的径流深。若干旱年份与湿润年份有明显的系统误差，可调整WUM和C值。减小WUM可使少雨季节的蒸散发量减小，而对很干旱的季节则影响不大。加大C值可使很干旱季节的蒸散发量增加，而对雨季则影响不大。 比较年内干旱季与湿润季之间的差别。在南方流域，主要分析伏旱季的蒸散发计算是否正确。如伏旱季以后的初次洪水具有明显的系统误差，则应调整WUM和C值。 比较枯季地下径流。若枯季地下径流有系统偏大或偏小，则应调整KG和KSS值，即调整地下径流所占比重；若枯季地下径流的消退系统偏快或偏慢，则应调整KKG以调整其退水的快慢。 参数B和IMP在湿润地区并不敏感，可以通过比较小洪水的径流总量进行调整。 通过日模型的调试，可率定出第一、二层次的参数，第三层次和第四层次的部分参数。其余参数可通过次洪模型进行调试。

对次洪模型： 次洪模型以△t为时段长。新安江模型参数中有一些参数与时段长无关，如K、WM、WUM、WLM、WDM、B、C、IMP和EX。这些参数在日模型与次洪模型中可以通用。换言之，这类参数在日模型中调试好后，在次洪模型中就不需要再进行调试。 有些参数则与时段长有关，如SM、KG、KSS、KKG、KKSS、L、CR。 其中KG、KSS、KKG、KKSS虽然与时段长有关，但它们在次洪模型中的值与在日模型中的值存在一定的转化关系，实际上与日模型可以相互通用。 L、CR在日模型与次洪模型中不通用，需重新进行调试。 SM由于受降雨资料在时段内被均化的影响很大，时段愈短，所得的SM值愈正确；时段过长，例如以日作为时段长，则求得的SM值将显著减小。因此，SM值在日模型与次洪水模型中不通用。 当SM改变后，KG与KSS也相应地要改变。根据经验，令日模型的参数为$SD_D、KG_D、KI_D$，次洪水模型的参数为$SM_{\triangle t}、KG_{\triangle t}、KSS_{\triangle t}$，则由$SM_D•(KG_D+KSS_D)=SM_{\triangle t}\cdot (KG_{\triangle t}+KSS_{\triangle t}）和KG_D/KSS_D = KG_{\triangle t}/KSS_{\triangle t}$，可解得 $$KG_{\triangle t}= KG_D•SM_D/SM_{\triangle t};KSS_{\triangle t}= KI_D•SM_D/SM_{\triangle t}$$

次洪水模型在调试时通常以洪水总量、洪峰值、峰现时间按许可误差统计合格率最高和确定性系数DC最大为目标函数。许可误差是依据预报精度的使用要求和实际预报技术水平等综合确定的误差允许范围，由于洪水预报方法和预报要素的不同，对许可误差规定亦不同，这部分可以直接查资料。

比较洪水径流总量。 影响计算次洪径流总量的主要因素除降雨外，另外就是流域初始蓄水量W。有时需调整产流参数，在产流参数确定的情况下，可以 通过调整水源的比重来改善次洪径流量。可调整SM和KG，这两个参数值越大，地下径流的比重也就越大，次洪径流量就会减小。

比较洪峰值。 洪峰流量基本上由地面径流和壤中流组成，它主要取决于单位线、SM、KKSS等参数。首先调整SM，当SM确定后，可调整KKSS参数，单位线。在以河道汇流为主的区间流域，特别是大河，也可调整参数XE。

比较峰现时间。 主要调整KE、河段数n和滞时L。减小KE可使计算洪峰提前，否则相反；增大河段数n和滞时L可使计算洪峰滞后。

通过日模型和次洪模型调试，若发现某些参数不一致或明显不合理，应协调或调整这些参数后重新进行日模型和次洪模型的计算，优选出一组合理的最佳的参数。

人机交互率定参数方法的主要特点是在寻找新的参数值时，根据人们的先验知识去判断、估计。显然这种判断、估计不是唯一的，有时甚至可能是错误的。人机交互率定参数方法虽简单易行，但当无先验知识时，调试需要花费大量的计算机机时。

所以现在文献中一般都是采用人工经验和自动率定相结合的方式，现在在文献中比较少再看到经典的罗森布朗克法和单纯形法了，毕竟太经典，也没啥改进空间，目前比较常采用进化算法对参数进行率定。本项目使用遗传算法进行计算。关于遗传算法的基本原理可以参考资料。讲解比较详细。网上可以使用的程序包较多，因为java下常用的多目标优化程序包有MOEAFramework，所以算是为了和它一致，选择一个由同一作者开发的python版本多目标优化程序来使用——Platypus。

结合模型率定简单介绍遗传算法：

  • 初始化种群 随机地从参数的搜寻空间中选取m个点，然后进行编码，初始化种群。 编码方式较多，模型里都是实数运算，所以采用一种实值编码的方式即可。 模型的参数都是有物理意义的，所以都有相应地合理取值范围，因此在计算时，应设置模型参数的可行域，这样也可以缩小寻优空间，使函数的寻优过程能更快地收敛。 这部分可以参考《水文预报》书本上对每个参数的描述及给的参考取值范围，归纳一下如下： 系数类的参数大部分范围的取值在[0,1]之间，K值也有大于1的，不过一般小于1； 时间类的参数有KE和L，KE按运算时段取值即可，而滞时则需要观察流域的流量过程线进行分析判断，在本项目中，L是河网汇流计算使用的参数，因此也不会太长，不知道如何判断的话，可以采用试算的方式； 剩下的是量纲为mm的参数，书上没有推荐的是SM，即表层自由水蓄水容量曲线方次，这个也可以尝试先给出一个相对宽松的范围。
  • 计算适应度 计算评价种群个体的适应度，即计算对应的目标函数值，模型里就是计算给定的场次洪水的模型计算结果的与实际观测结果对比的结果，涉及到多场次和多种对比指标，对于多场次可以求平均，对于目标函数，可以设置不同的目标函数进行计算，在设置目标函数之后，需要设置终止优化的准则，设置最大迭代次数的情况较多，也可以设置目标函数值得收敛容差，即两次优化结果之间差距很小，就可以终止运算了，也可以设置参数迭代步长收敛容差；
  • 进化操作 进化操作，包括选择，重组和变异。这些运算的具体形式都有多种，一般选择当前比较常用的方法，比如在多目标遗传算法中常用NSGAII算法，现在常见的重组算法是SBX模拟二进制交叉算法和polynomial mutation，注意不同编码对应不同的重组算法，编码和重组算法是需要正确对应起来的。因此在设置运算时，可以设置这些，并且一般还需要指定一些运算参数，比如交叉率、变异率等等，当然也可以使用默认值；
  • 形成新种群 进化操作完成之后，得到的是育种后代，为了保持种群规模，可以执行重插入操作，像NSGAII算法等都有精英保留的运算，因此使用框架即可，保留这些新个体即形成了下一代种群，然后循环计算，直至达到前面设置好的终止条件，即优化结束。

前面提到，在实际调参过程中，辅助一些人工经验会更好，可是我并没有什么经验，只能依靠书本上的知识，或者说如果在一个很少有人做的流域开展水文预报，那这个时候从问题本身去进行分析可能也是必不可少的，现在很多可视化的工具可以辅助我们进行分析。因为可能涉及高维，所以在可视化过程中可能需要一些降维操作。

参与贡献

  1. Fork 本项目
  2. 新建 Feat_xxx 分支
  3. 提交代码
  4. 新建 Pull Request
Clone this wiki locally