第7章 数值模型设计、应用及水文地质勘查

7.3.8.3模型识别方法及注意问题

有了上述工作的准备，可以试运行程序。

（1）模型识别方法及识别的目标
模型识别是利用某时间段已知的或试探性给定的边界条件、源汇项及其地下水动态——通常是地下水头/水位和流量（泉和自流井等），通过某种方法确定含水系统的全部参数（K，u，α，…），以及控制不足的边界条件和调查不详的源汇量的可能范围内的修正（后者要特别慎重），甚至要修正原设计的地下水流控制方程：模型维数、稳定/不稳定流、达西/非达西流等，这个过程也称为解逆问题。
解逆问题的间接解法又分为模拟-校正法（人工调整法）和自动优选法两类。前者是模拟者先按照种种方法给定一组参数的初始估计值，以运行程序，获得模拟的地下水动态，然后将其与实测值对比，根据它们的差异来修正参数，再运行程序……如此反复执行，直到满意为止。因此这种方法能够充分发挥模拟者对地区水文地质条件的理解，把握调整在上、下限参数值范围之内;优选法能摆脱模拟者个人认识的局限性，具有相对客观的优点。然而如果完全使用优选法，它是多解的（我们经历过），且其结果往往不可接受。因此，目前大多采用两者相结合的方法，前期采用模拟-校正法大幅度的调整参数，到后期精细调整期，该法已无能为力，而采用优选法（详见第6章）。
模型识别的目标是使诸模拟值趋于实测值，通常以误差的平方和为其目标函数，即（陈崇希等，2003b，2005b）

或以均方根误差（RMS)即标准差为其目标函数，即

式中:，，……，分别为待求的各个水文地质参数；、分别为对比用的观测孔数和对比时段数;为水头权因子；为模型模拟的号观测孔时段的水头；为号观测孔时段的实测水头；、、分别为对比用的泉数、自流井数和对比时段数；为流量权因子；为模型模拟的号泉或自流井，时段的流量(m3/d)；为i号泉或自流井，j时段的实测流量(m3/d)；为观测孔水头拟合所占的比重；为泉和自流井的流量拟合所占的比重。

​

 有些文献所建模型识别的目标函数没有设权因子，实质上是将所有的权因子取其为1。设权因子是由于:
1）选用的拟合点一般分布不均，大范围内仅有一个拟合点的权因子应给大值，密布在同一区的拟合点的权因子应给小值。
2）如果拟合点分布比较均匀，由于它们控制的部位/元素上的差别（是否为核心模拟区等）权因子应赋以不同值；重要的泉、自流井、河流、“岩性天窗“等附近及其在水均衡中占重要地位的地下水补给/排泄区等，权因子应赋以大值。
3）在不同的调参阶段和对象，可以突出某些拟合点的权因子，以凸显其反应，例如拟合水头或流量等。
4）水头与流量的量纲不同，因此其权重入与（1一1）的取值要考虑。
在各参数的上下限约束条件下，求取目标函数G最小，若此时的水位拟合差未达到足够小，则需分析其原因，调整相关的分区及参数值，再次求取目标函数最小。若拟合曲线第一对比时段的水头/流量误差过大，这可能是初始水头不合理所致，要调整初始值。如此反复模拟、分析，直至求得满意的结果为止。如果始终不能达到合理的结果，则要分析所建的基本模型、诸元素的模拟方法、边界条件及源汇项等是否有问题及数据输入是否有误（如运城），甚至所提供的水文地质条件是否有误（如渭北东部岩溶水模型）。整个识别的流程如图7.3.1所示。
（2）模型识别中的几个问题
1）对于法，模型运转初期若干个△t 的模拟结果不宜拟合求参。对初始条件采用法的模型，模型运转初期若干个时阶，抽水井附近结点的模拟水位通常具有明显的波动。这种波动显然违背了抽水过程的水力实质，它是数值解法本身带来的误差（陈崇希等，1985）。从这个角度来说，抽水初期的数值模拟结果不宜用来与实测水位拟合求参数。

2）降雨入渗补给系数α大于1。一般地说，降雨入渗补给系数α小于等于1。就总体而言，此约束条件没有问题。但在山前倾斜平原的上部，特别是在沟谷出口附近或扇间洼地处，由于山区坡面流的作用使局部洼地的入渗补给量远远大于自身降雨的入渗量，致使某时段α大于1（如山东、山西某地模拟过程中曾遇到此情况）。
3）河水单位长度入渗补给量。对于不能用4.10节中“河流与地下水间补给、排泄的模拟“的小河，其河水单位长度入渗补给量q通常采用q=α，Q计算（Q.为河流的流量；αr为单位长度河流入渗补给系数）。然而人们发现，不仅α。随空间变化，也随时间而变。α，值有两个因素起作用:
a.如黑河中游地区的梨园河，当新改道时其渗漏量很大，随着时间延续淤积形成，其渗漏量减小。
b、随季节降雨强度形成河床流量的差异而变，在洪水期，河流量剧增，对于矩形断面的均匀流，流量与河流水深的5/3次方成正比，而河流的入渗补给量与垂向水力梯度成正比，两者随河水位上升的增速不等，前者要快得多，这一点要特别注意（山西某模拟者也发现此问题）。
4）不允许潜水位高于地面。若某时阶解得某格/结点潜水位高于地面，一定要核对地面标高及影响人渗强度的诸因素，并实地核对此事件出现的可能性。若是否定的，则计算机程序应将潜水位自动控制在地面标高处，且其地面下游可能有更多的入渗量。

5）“拟合情况良好，说明所建模型可靠“。这是模拟报告常见的提法。“拟合是否满意“目前尚难以提出合理的准则，虽然有人尝试提出过。评审数值模拟技术报告，首要的是水文地质条件概化的合理性及模拟方法的科学性，在此基础上才能说“拟合是否满意“，因为这是多解的，特别是拟合点不够多、拟合时间不够长的条件。如果模拟是为了地下水开发利用的规划，那么模拟预测应在100年的量级上，而当前的模型识别时段一般不超过10年，大量的为5年以内，甚至仅有1~2年实测时段。依笔者的意见，应在3~5年以上。
6）“每个参数分区至少应有一个观测孔“吗？对于大尺度模型，例如近10年的中国地质调查局地下水资源-地质环境保护项目的模型及含有弱透水层的地面沉降模型以及其他三维流等模型，这是不可能的。因为仅就弱透水层而言，除了专门研究外基本上无观测孔。实际上这个“至少应有“也是没有必要的，因为某分区参数的改变也会在邻区或更远处的观测孔水位动态得以反映（与模型识别时段长度及参数分布等因素有关）。同样，弱透水层参数（如垂向渗透系数K'，）的变化，会明显影响相邻含水层观测孔水头的动态。当然，如果不计成本，观测点愈多愈好。然而反过来，我们却极少见到利用泉和自流井等的流量来拟合的报道。虽然拟合流量的难度比水位要困难得多，但这很重要，以降低解的不唯一性，增加数值模型的可靠性；我们也极少见到利用混合观测孔水位来拟合的报道。另外，延长观测时间/增加拟合时间长度（虽然通常提为“至少一个水文年“，但“一个水文年“实在不足）也极其重要（比多打观测孔要节省经费），它增加了模型的可靠性。
7）拟合分析-调整参数过程，不要一个孔、一个孔地分析，特别是调参之初。而是一个地貌单元（如整个洪冲积扇），或延伸至其前缘的冲积平原的一组观测点综合分析-调整，因为它们的参数与动态互相联系。此外，有时不宜一个参数、一个参数地调整，而需要两个或三个参数联合调整，特别是模型识别的中后期。
8）拟合分析-调整参数，是系统工程，是对水文地质条件的深化认识过程。由于参数初始估计值的赋值-般难以接近实际值，因此第一次调参往往难以下手。一种可参考的方法是:选择某范围或全区，轮流微动某个参数群（如水平渗透系数K、给水度p、垂向渗透系数K.、入渗补给系数α等），即交替乘以1.1~1.2及0.8~0.9，小步调整其初始估计值。并观察、分析其模拟水头/流量的反应，以获得初步的经验。如此还可相对缩小范围进行，以获得进一步的参数一动态特征或规律。虽然目前计算机的模拟速度很快，但我们要避免“单纯运算“的做法，要重视分析模拟结果，分析模拟结果的时间通常要多于甚至远多于模拟的时间。
9）关于稳定流拟合。不少报告单独给出稳定流所得的拟合情况，而且往往将所有拟合水位置于一张小比例尺的图中，难以看出拟合的差值，此做法不妥，因为“稳定流“与实际条件不符；有的报告仅给出某一时刻的拟合情况，也欠妥。“模型识别是建立数值模型中最重要的步骤“，也是费时最多的环节，当模型设计合理，模拟方法正确的前提下，模型的可靠性主要取决于此。因此所有实测点的整个拟合过程（曲线）都应给出，指明测点的位置，并分析之。
10）把握流场特征。在拟合分析-调整分区和参数的过程中，模拟者务必把握在边界条件和源汇分布控制下区域的流场特征，特别是三维流场，在此基础上，根据Darcy定律和水均衡原理来调整分区和参数。调整拟合曲线，重在动态平均值的趋势（升降斜率）的一致，其次是波幅及峰谷出现时间的拟合。如果仅调整一个参数，一般会有多解。例如:某观测孔模拟的水头动态曲线降速明显偏大，可以加大该孔水平上游的渗透系数/导水系数，或减小该孔水平下游的渗透系数/导水系数，或加大该孔垂向上游的垂直渗透系数，或减小该孔垂向下游的垂直渗透系数，或者（如果处于诸潜水层入渗补给区）加大入渗补给系数/补给量，或者（如果处于潜水层蒸发排泄区）减小蒸发排泄系数/排泄量（对于涉及水均衡的后两者，还应进行水均衡诸项变化的分析），给水度p也能调整水头动态曲线的升降速度等，或者上述诸点的综合。仅分析至此是不足的，必须进一步分析，每一个调整方案对附近其他拟合点动态的影响，是有利于还是不利于拟合。显然，仅针对此孔此时段而言，有多种方案。
11）不存在某时刻平面等水头线的拟合。笔者认为:拟合，只是对每一观测点动态曲线的拟合，不存在某时刻平面等水头线的拟合﹔对于平面等水头线，除拟合点外，其他曲线段都是插值而得，不同的插值方法，拟合结果不同。某些实例的地下水动态拟合得不太理想，但平面等水头线却拟合得不错。“
12）地下水模拟元素的模拟方法的失真。在短期的模型识别过程中有可能通过其他参数的调整而暂时在一定程度上掩盖其模拟方法上的问题，然而对于长期的预测，其结果将会与实际的差距拉大。
 

7.3.8.4“拟合总是不佳“的分析
当拟合点较多，特别是不回避关键拟合点（泉和自流井的流量和混合观测孔水位等）且拟合时间较长的条件下，“拟合总是不佳“可能存在以下四方面的原因。前两点是结构性的，不改变基本模型及模拟/处理方法，无法提高拟合结果。
（1）所建模型的基本类型与实际问题不符
如7.2.1.1节所讨论的，如果一个实际为三维流的水文地质条件，采用了二维流模型刻画；一个弱透水层为各向异性的多层含水系统三维流的问题，采用了准三维流模型刻画；一个多层的三维流问题，采用了分层过少的模型刻画（与抽水井、观测孔的类型及位置有关）；一个不稳定流问题，采用了稳定流模型刻画等，这类问题若不改变基本模型，不可能拟合好。反过来，如果拟合结果可以接受，并不能说明所建模型的基本类型与实际问题相符。可能是现有的用来拟合观测孔与抽水井的类型及位置尚未明显地暴露出三维流特征。一旦增加观测孔或加入其他类型的抽水井，就可能显示出三维流特征。这就是笔者为什么强调:“拟合好或差不是检验建模合理性的唯一标准，更重要的也是最根本的是视水文地质实体的概化是否合理，模拟方法是否符合机理，正因为如此，我们十分强调对水文地质条件的正确分析和建模论证“（陈崇希等，1999c）。
（2）地下水元素的模拟方法失真或仿真性不佳
“防止模拟失真，提高仿真性“是数值模拟的核心（陈崇希，2003a）。不言而喻，一个失真的或仿真性不佳的模拟成果作为地下水资源的可持续开采量或防止地质环境恶化的依据，将会产生什么后果？模型设计者必须深入分析水文地质条件，在诸类边界条件和源汇控制下的地下水运动的机制，科学地分析现有软件对地下水诸元素模拟的合理性；对于不科学的模拟方法，必须予以纠正，一味地遵照、模仿前人的做法要不得。Andersonet al. （1992）也认为:模拟者“应判断模拟过程中使用的方法是否合理“。但遗憾的是，我国应用性的地下水模拟技术报告，通常见不到对其所采用软件有关诸地下水元素模拟方法科学性的评价，拿来就用，似乎什么问题也没有。我国的有关专著、论文也极少评论国内外软件、文献有关模拟方法的仿真性，更少涉及模拟方法的改进，特别是经常遇到的地下水流诸模拟元素（如:泉，非完整河，混合抽水井，自流井，混合观测孔，降雨入渗补给的滞后性等）。科学发展的核心是创新，跟踪不是创新。
与其说采用一个失真的模拟成果作为地下水资源可持续开采量或防止地质环境恶化的依据来指导区域地下水管理，不如说暂时模拟不了该区的有关地下水问题，而作为攻关问题有待解决。与其说做出一个错误的结论，不如暂时不做出结论，有待进一步研究。这就是笔者强调，“防止模拟失真，提高仿真性是数值模拟的核心“（陈崇希，2003a）的初衷。
（3）输入数据的系列长度及其可靠性问题
这可能属于勘查阶段（精度）问题，刻画水文地质条件的信息偏少，或者未能充分收集已有的资料。后者应避免，关于前者，例如；
1地质孔间距过大，含水系统的结构、岩性难以控制；2）水文地质点偏少或过少；
3）边界条件控制资料缺乏或过少；
4）抽水井、观测孔及自流井等必要信息缺乏；5）抽水井、泉水流量数据不足或精度不高；
6）主要河流的水位/流量缺乏实测值或者推测值误差过大；
7）地形高程（主要是潜水浅埋区）、泉口和自流井口等标高的精度不高。
上述诸点基于水文地质勘查尚属于初级阶段，拟合不佳是可以理解的。但要求:1）模拟者完整地、全面地搜集资料；
2）各类资料对所建模型可靠性的地位、重要性，模拟者应做到心中有数。3）对一些关键性数据应尽量补充搜集、测量。
（4）模拟者的学术技术基础、经验、责任心等问题
水文地质模拟者的水平（理论、方法和经验等）可高可低，但“实事求是“及“责任心“不可无。有了后两点素质，其他技术方面的水平是可以逐步提高的。例如:大量的观测孔不去搜集或随意去掉；渤海湾海水的最大深度才24m，20m海水等深线离海岸距离超过50km之遥，而所建沙岸滨海区海水入侵模型却将其刻画为直立海岸——切穿整个模型至底面，这是技术问题还是责任心问题？模型中存在混合抽水井、混合观测孔，但论文、专著中却一字不提它们的存在，更不说明是如何模拟或处理的；模型中存在自流井，但论文、专著中却不提及这些重要的水文地质信息，以回避其模拟方法；研究区存在泉的，不去区分属于上升泉还是下降泉这个最基本的水文地质特征，更不可能由此提出符合机理的模拟方法……总之，如果初次遇到某些新问题不会模拟是完全可以理解的，但是应该说明如何处理；将问题掩盖的做法，是科学的大敌。
模拟工作的“短平快“是要不得的，模型识别过程中，既要“时时把握从水文地质条件来调整有关参数“，又要“模型识别过程是对水文地质条件的再认识“，要把握好看来是矛盾的两个对立方面，然而这是对模拟者综合素质水平的评价准则。模拟过程需要花费很多时间和一定的经费，但与水文地质勘查相比，后者所投入的经费更多，工期更长，工作更艰苦。模拟是对前期所有水文地质勘查成果的综合、升华与量化提高，许多重要结论（地下水可持续开采量的评价、地面沉降、海水入侵、生态环境保护，或矿山排水、岩溶塌陷及水利工程相关问题）都来自模拟结果，模拟者应负责任地完成此项任务，这也是对前人从野外现场艰苦获得水文地质勘查成果的尊重。
如果已充分收集现有的资料，所设计的概念模型反映了实际的水文地质条件，采用符合流动机理的模拟方法，由熟悉水文地质条件、掌握模拟方法、具备模拟经验和责任心强的高素质模拟工作者精心地操作，那么拟合不佳，则反映的是勘查阶段（认识自然的阶段）问题。
 

7.3.8.5灵敏度分析
模型识别所得的参数，是在一定边界条件及源汇（还包含模型分层、参数分区以及模拟方法）等条件下所得的。然而它们不可能如此准确，例如地下水开采量改变5%，那么对模拟响应量（模拟的地下水头或流量等）将带来多大影响？又如，整个模拟区的某主干河流长180km，而只有3个水文站有实测河水位资料，其他沿河格/结点的河水位都是插值而得，它不可能是准确的，那么若水位改变0.1m，对模拟响应量带来多大影响？若影响大，是灵敏的，则识别模型的可靠性差；反之，是不灵敏的，识别模型的可靠性好。
就模型识别所得的参数而言，这些参数也是不可能完全准确的，若某个/些参数改变10%，对模拟响应量带来多大影响？若影响大，是灵敏的，则识别的模型可靠性差；反之，是不灵敏的，识别的模型可靠性好。
灵敏度分析是对所识别模型不确定性/可靠性的一种量化评估。如何量化？定义灵敏度（系数）为

式中:S，为改变第j个参数p，（如渗透系数、给水度、入渗补给系数等）对第i个观测点模拟量H；（如水头、流量等）的灵敏度。N。为观测点数（也可为结点数）；N为参数总数。具体做法是将上式差分化。
上式的灵敏度具有不同量纲，为此也可改写成统一的量纲为一的形式，即

​

这样，式（7.3.10）似乎合理了，但是该式含义是当某参数p，改变某百分数时，H的改变量占原量的百分数。然而，同一模型在低海拔的华北平原与高海拔的青藏高原的灵敏度就不相同，因为水头的“基数”不同，前者小，相对误差大﹔后者反之。为避免此缺点，也可将模拟响应量H的“基数“改为模拟区水头/流量的最高值与最低点值之差（H.maxH.min），有时参数p也可如此处理。即

还可采取其他形式灵敏度系数的定义。
为了综合性分析，也可采用上述模型识别的目标函数G的形式（7.3.7）代替单一i点的模拟响应量H；。
上述是逐个改变一个参数对模拟响应量H；的影响（单因素灵敏度分析），这并非灵敏度分析的全部，如果要获得更加全面的灵敏度分析信息，还可同时改变两个、三个……参数（多因素灵敏度分析），特别是某些相关的组合参数，如:潜水层的渗透系数与重力给水度；承压层的渗透系数与弹性给水度﹔潜水层的重力给水度与入渗补给系数﹔潜水层的渗透系数、重力给水度与入渗补给系数﹔潜水层的渗透系数、重力给水度、相邻弱透水层的垂向渗透系数与入渗补给系数等以获得模拟响应量H；的反应。
另外，计算灵敏度系数时，其中的导数要差分化，其参数改变量△p的取值大小也会影响灵敏度系数的大小，因为作为实际复杂含水系统而言，灵敏度系数一般是非线性的。如此，若要获得更多的灵敏度信息，就要取不同的△p，其工作量可想而知。
对于不稳定流的不同模拟时间长度，一般具有不同的灵敏度系数，当然模型识别的末时刻是应该计算的，至于其他时间的信息视需要而定。模型识别末时刻的灵敏度系数不一定是最大的。
实际上，模型识别自动优选参数过程包含了相当的灵敏度系数的计算。
模型识别所得的参数，是在一定边界条件及源汇项（还包含模型分层、参数分区以及水文地质诸元素的模拟方法）等条件下所得的。
尽管做了上述如此多的灵敏度系数的计算（且尚未包含模型分层、参数分区以及模拟方法等），但如何评价模型的可靠性，仍然难以给出客观的标准，有待进一步研究。特别是，模型更重要的用途是预测-管理，然而对预测结果的可靠性/误差的分析，目前尚无合适的方法。另外，灵敏度分析目标是评估模型识别所得参数的可靠性，而不能提高模拟预测的可靠性，且工作量又大。更重要的是，地下水诸元素模拟失真或仿真性不佳带来的模拟响应量的误差将如何反映到参数的灵敏度系数上？基于此（当前国内对于实际复杂含水系统，例如近几年的地质调查项目，参数接近或超过100个），笔者建议:将主要精力放在水文地质条件分析，合理建模，提高仿真技术及拟合上。
 

7.3.8.6模型检验
如果研究区的实测动态资料较长，可将时段分为模型识别和模型检验两个阶段，将模型识别获得的参数用来预测，并与其后的另一时段实测动态资料进行对比，以说明其可靠性，此即所谓模型检验。另一种做法是，将全部实测动态时段作为模型识别时段，即无模型检验，特别是实测动态时段的前后期具有不同的源汇/边界作用时（例如前期无降水，后期有降水﹔前期某区无抽水，后期有抽水﹔前期某区无灌溉，后期有灌溉﹔前期某区为承压状态，后期为无压状态等），将模型识别和模型检验合二为一更为有利。由于延长模型识别时段会提高模型识别的可靠性，而模型检验只是判断模型识别结果的好与差，而不能改善模型识别的结果，况且模型识别的拟合情况也能判断其结果的好与差。基于此，并考虑到当前我国的观测数据系列一般都不够长，因此笔者更倾向于将全部实测动态时段作为模型识别时段。
 

​

/d