-
Notifications
You must be signed in to change notification settings - Fork 109
模型参数率定
新安江模型是概念性模型，在进行计算与率定之前，需要对模型参数的物理过程归纳，然后初定参数范围。关于水文模型参数的率定需要人工经验，很多文献都有提及，一直以来都是想要去尝试改进的点，能力所限就不展开了。
- 蒸散发能力折算系数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.15
0.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作为指标，其数学表达式为：
根据我国SL250-200《水文情报预报规范》，调试时通常以年径流绝对误差最小及确定性系数DC最大为目标函数。本项目代码部分就先以此为示例进行参数率定。不过不影响接下来的说明，另外，后期有需求再补充其他寻优目标。
许可误差是依据预报精度的使用要求和实际预报技术水平等综合确定的误差允许范围。由于洪水预报方法和预报要素的不同，对许可误差规定亦不同。包括洪峰预报许可误差，峰现时间预报许可误差，径流深预报许可误差和过程预报许可误差。这些根据书本上的说明，还是挺繁琐的。
预报项目精度评定规定:
定义合格预报为合格预报次数与预报总次数之比，表示多次预报总体的精度水平。一次预报的误差小于许可误差时称为合格预报。
预报项目精度按合格率QR或确定性系数DC的大小分为甲乙丙三个等级。
当进行实时作业预报时，还需要对洪峰预报时效等级进行分析。用洪峰预报时效性系数描述:
形式上包括人机交互率定（假设参数，机器运算，比较实测值与计算值，人工经验调整参数，直到最优）和优化算法自动率定。
新安江模型参数率定一般分为“日模型”和“次洪模型”两大部分。模型可靠性取决于模型研制中使用的水文资料的质量和代表性。洪水预报方案要求使用样本数量不少于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}$，可解得
次洪水模型在调试时通常以洪水总量、洪峰值、峰现时间按许可误差统计合格率最高和确定性系数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算法等都有精英保留的运算，因此使用框架即可，保留这些新个体即形成了下一代种群，然后循环计算，直至达到前面设置好的终止条件，即优化结束。
前面提到，在实际调参过程中，辅助一些人工经验会更好，可是我并没有什么经验，只能依靠书本上的知识，或者说如果在一个很少有人做的流域开展水文预报，那这个时候从问题本身去进行分析可能也是必不可少的，现在很多可视化的工具可以辅助我们进行分析。因为可能涉及高维，所以在可视化过程中可能需要一些降维操作。