-
Notifications
You must be signed in to change notification settings - Fork 109
模型核心算法
在进行产流计算之前，首先需要计算前期土壤蓄水容量。利用前期降雨蒸发数据，采用流域蒸发计算模型逐日计算。
- 当降雨 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 =
注意计算时，一定有$a>W_0$，$P-E>R$，从流域蓄水容量曲线上很容易看出，如果有不符合上述条件的情况出现，一定要检查是计算机计算的小误差还是因为自己计算或者设置参数有问题导致的。
关于不透水面积，书上使用计算公式Wmm'=(1+B)/(1-IMP) * Wm来处理，那得到的R应该就是已经考虑不透水面积的了，即R是针对整个流域的径流深。不过书上没有给出这块的证明，也没有这部分的推导，暂时我也没有推这部分，觉得稍微有点问题。是不是直接在透水面积上计算，然后最后再将最终计算结果转换到流域面积上更合理一些？目前先按照书本上的做了。
水源划分的原因在于后面汇流阶段不同径流成分的水滴其汇流表现不同，汇流时间不一样，因此需要区别处理。 这部分书上解释地不是太详细，这里给出个人的理解，如下所述。
首先流域产流用的标志是流域上一点的土壤含水量超过田间持水量。接下来，对于流域而言，土壤好比地下径流蓄水库，遵循蓄满产流的定义，土壤在达到田间持水量之后，会直接补给地下水，直至该“水库”蓄满，然后蓄满产流，产生地表径流。而该“水库”蓄满的含义就是指产流区域上一点的水量达到或高于自由水蓄水容量，高出的部分就是地表径流。（这块和前面蒸发分几层没有关系，它们互相之间应该是不对应的，因为这里是对自由水讨论的，前面蒸散发是对张力水讨论的）
而自由水蓄水容量在产流面积上分布也是不均匀的，考虑用和产流计算时同样的技巧，即用概率的方式表达，表达的形式和之前产流的方式是一致的————产流面积上的自由水蓄水容量曲线，那么就可以分出地表水和地下水了。
具体的表达如下：
从上式可以推出(推导过程和产流计算时候的一致): 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，就得到壤中流；自由水深×地下径流出流系数×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是模型参数。
- 当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} $$
关于对卷积的理解可参考说明。 我觉得上式的理解难点有两处，一是卷积运算里函数的取反，即翻转，以及滑动叠加，这两个运算的结合； 二是卷积运算是一个序列运算，每个时刻的运算结果不仅考虑了当前时刻的输入响应，还跟过去所有时刻的输入响应有关，而跟哪些过去有关，是利用卷操作来指定的。 放到净雨和单位线的卷积，就是说一个净雨它的作用是在单位线每个时段上都留下一个输出，从第一个时段净雨输入到单位线中，它不仅在第一个时段留下了输出，它给后续所有时段都留下了输出，那么到最后一个时段，它也给它后续多个时段留下了输出。所以才有了上述表达式。
汇流研究的其实主要是两类问题，一是刚刚提到的已知净雨过程和流域单位线推求出口断面流量过程的“预测”问题；另一个镜像的问题就是已知净雨过程和相应地出口断面流量，如何反演单位线的“识别”问题。所以问题的重点其实在于如何要从历史数据中反演出单位线。
推求单位线有从不同角度，按不同理论的很多种方法。
比较经典的是分析法求时段单位线。
- 选择若干场时空分布均匀，历时较短的降雨形成的单峰洪水来分析；
- 点绘流量过程线后分割基流及前期洪水退水，计算本次洪水的直接径流深（具体过程见下一节分析）；
- 对该场次降雨进行产流计算的结果作为净雨过程；
- 根据净雨郭晨估计对应的地面径流流量过程线，利用以下公式推算单位线： $$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；
- 修正单位线：检查总量是不是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值，也可以用优选的方法来计算。
本项目给出的代码就先以上面谈及的单位线为例。
需要注意，在计算过程中由于需要利用实测数据进行参数识别，而前面计算的产流结果是分水源的，所以需要对实测数据进行分水源分析。为了分析地面径流，需要把地下径流从实测数据中割除出去。《工程水文学》上介绍的方法有较为简便的斜线分割法，还有经验公式法。这部分的编程处理略微有些麻烦，具体的数学表达还需考虑。
应用单位线时还有几类实际问题需要补充说明：
- 单位时段不同的单位线之间的转换
- 实际中时变和非线性特性的影响
- 无水文资料地区
应用单位线计算汇流时，常常因净雨量时段长和单位线的时段长不一致而引起误差。因此需要不同各时段长的单位线转换。
如何理解？ 从S曲线开始说起。前面提到流域时段单位线实际上是单位矩形函数的系统响应。 单位矩形函数可以表达为两个阶跃函数的差（那么不同时段长对应的就是阶跃函数的相对位置不同），阶跃函数求导就是冲激函数，所以容易推求冲激响应积分就是阶跃响应，那么阶跃响应的差就是时段单位线。 定义流域S曲线：$S(t)=\int_0^t u(0,t)$，即单位阶跃响应，则时段单位线$u(\triangle t,t)=S(t)-S(t-\triangle t)$。 流域S曲线的定义就是流域上分布均匀，且一直维持在1个单位强度的净雨所形成的流域出口断面流量过程。
注意实际中，都是离散条件下的计算：
了解了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$时段对应的单位线了。公式如下：
对于瞬时单位线，由于实际中运算是离散的，所以一般都是先把瞬时单位线转换为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}$$
联立上述各式，容易求得:
书中新安江模型里，令壤中流和地下径流机制相同，因此总结上述内容。壤中流RSS和地下径流RG分别流入壤中流水库和地下径流水库，经过蓄水库的调蓄（记两水库的消退系数为KKSS、KKG），成为地下水对河网的总入流TRSS、TRG。计算公式如下：
模型的消退系数也是按日模型给定的，所以进行实时洪水预报时，计算时段小于24h时，需要作相应地变化。
设模型计算所取时段长为△t（h），R为△t内的净雨，则：
计算完毕之后，汇总：TR(t)=TRS(t)+TRSS(t)+TRG(t)
单元面积河网汇流这部分，可以考虑，也可以不考虑，看前面怎么理解单位线吧，《水文预报》上说由于单元流域的面积一般不打而且河道较短，对水流运动的调蓄作用小，可以忽略，合并在前面的单位线中。
如果考虑的话，一般采用线性水库或者滞后演算法进行。滞后演算法等价于线性渠和线性水库的串联，《水文预报》书上的公式是：
河道演算。新安江模型一般采用马斯京根法进行河道演算。河道流量演算是以圣维南方程组为理论基础，利用上断面流量过程演算下断面流量过程。马斯京根法是对圣维南方程组近似得到的水量平衡方程和槽蓄方程进行求解的一种方法。现在直接求解圣维南方程的数值解法也比较多，也较常应用。关于圣维南方程的理解，可以参考这个，和前面一样，如果想复习下相关的微积分基础知识，推荐参考这个，关于纳维尔斯托克斯方程和圣维南方程之间的关系，简单查到了这个。NS方程实在是有点复杂，需要专门学流体力学的人来解释了，我是不太会了，这里因为要用水文方法求解，就从圣维南方程开始说起，其他的就简单带过。
关于流体，一开始学的时候就讲了是按欧拉描述学的，分析的是场，而学刚体的时候是拉格朗日的对物体的描述，有较大区别。因为一块流体是会变形的，太复杂，所以想象的时候，想象空间中向量场的可视化会好一点。在一个空间中，质量守恒，能量守恒，动量守恒。圣维南方程就是个一维的平衡方程组，是相对简单的。
圣维南方程用于明渠水力学，基本假设是：流体密度是常数，满足连续性方程，且应力应变关系有牛顿流体的本构关系。（未完待续。。。）
水文法更简单，还简化了圣维南方程，得到了水量平衡方程和槽蓄方程。马斯京根法建立了马斯京根槽蓄曲线方程，将其与水量平衡方程联立，可得到马斯京根流量演算方程：
其中，C0，C1，C2都是KE，XE的函数。