在上述方程中,只有2个独立综合参数,这决定了3个参数无法同时反演确定。补给强度参数包括降水和渠系渗漏入渗补给系数、抽水流量有关参数和泉流量有关参数等。
5)对于非稳定流模型,由于把观测孔地下水位动态的时间序列作为识别依据,模型反演不一定需要观测孔的数目与待求参数的数目相同。以有限元模型反求渗透系数为例,从单元内节点之间关联系数的计算式可以看出,关联系数和渗透系数有成正比的关系,可以表示为
eeGPL??PLKe (4.5.20)
式中:Ke为单元渗透系数(假定为均质);?pL为取决于单元几何特征的形状引子。图4.19给出了一个三角形四边形混合有限元网格,边界节点均为定水头节点,含有0、1、2共3个内节点,含水层在A—B线左侧渗透系数为KL,右侧渗透系数为KR,则有限元模型是由3个方程构成的方程组
?KL?L?KR?R??t0?H0k?1?KL?01H1k?1?KR?02H2k?1?F0k (4.5.21)
k?1KL?10H0??KL?11??t1?H1k?1?F1k (4.5.22) k?1k?1KR?20H0??KR?22??t2?H2?F2k (4.5.23)
式中:?L和?R是节点0两侧的几何因子;?pL表示节点p和L之间关联系数中的几何因子;
?tp表示与时间和其他因素有关的系数。其中右端向量?F?可一直追溯到初始条件。在上述
方程组中,原待求变量加上需要反演确定的两个渗透参数,则共有5个待求变量,而独立方程只有3个。从定解条件上来看,待求2个参数KL、KR必须用至少2个节点的水头作为已知变量才能反演。
kA01KL2KRB
图4.19 含有5个内节点的有限元模型
这说明从固定时刻的定解条件上讲,待求参数的数目应与观测孔的数目一致。然而,以下情况仍然需要考虑:①已知节点1的观测水位动态,能否同时反演2个参数KL和KR?
即在已知H1的情况下,通过上述方程组能否求解出KL、KR、H0、H2。显然,由于待求变量有4个,而固定时刻的数值模型方程组中只有3个方程,因此通过第一个时步的运算是解决不了这个问题的。但是,在进入下一个时间步长之后,又会产生3个方程,2个待求变量(在第二个时刻:节点1水头已知,而节点0和节点2的水头未知)。加上第一个时刻的3个方程,现在共有6个方程,而待求变量也只有6个,因此应用唯一解。对于任何一次时步的推进,都会使方程的增加数大于待求变量的增加数,这种方程数和待求变量数的变化过程可用图4.20来表示。对于具有n个待求参数和m个观测孔的反演问题,由于没增加1个时步,方程数与待求变量数之间的差距就缩小m,因为只要时步数k满足km?n,模型反演问题就是唯一解或超定的。这导致一个重要的结论:少量的观测孔也可以反演大量的参数,只要时间剖分的时步数足够大。当时步大于临界值时,反演问题往往是超定的,可寻求一个使拟合误差最小的最优解。②用节点1的观测水位动态反演渗透参数KL和KR是否可靠?这个问题等价为:渗透系数KL和KR得取值很不敏感性分析进行判断。
6)模型识别需要紧密结合对水文地质条件和地下水流动机理得分析。对于数值模型,起求解的结果并不仅仅由参数决定,还取决于模型对物理现实的刻画程度。模型的简化性质,如用分区参数代替连续分布的参数、用平面二维流模型处理三维流问题等,
k?0:初始水头分布; n个待反求参数; m个已知观测节点 方程数:Ek=0 待求变量数:Uk=n k=k+1 新增Np个节点方程 新增Np?m个待求节点水待求变量数:方程数:Ek?Ek?1?Np Uk?Uk?1?Np?m Ek?Uk??Ek?1?Uk?1??m Ek?Uk?0 定解问题或超定问题
图 4.20 非稳定流数值模型那个反演的适定性随时间的变化趋势
往往对模型参数识别的可靠性和物理意义有影响。在某些工程领域,由于数值模型识别出来的参数往往和试验值相差甚远,有人提出“模型参数”的概念。造成这种现象的一个重要原因是物理意义的参数在简化模型中可能是和一些综合参数混合在一起的,综合参数中包含了不可测的一些简化系数或因子,并没有独立性。当这些系数或因子与物理参数同时被识别时,它们各自都是不在准确的,但它们形成的综合参数或无法用参数表示的综合效应可能是唯一的。
三、反求参数的优化方法
数值模型反演的方法包括直接法和间接法。直接法是吧已知的水头分布和水头变化率代人数值模型的方程组,直接计算出其中的参数。这种方法操作复杂而且存在很大的不稳定性,实际应用很少。目前使用的反演方法绝大多数为间接法,即先假定一组参数的取值进行模拟,然后比较计算结果和观测结果,如果误差较大则调整参数数值重复计算,直到满足拟合标准。
当参数的数目少于3个时,可以采用试算法,即目估计算的动态曲线与实测动态曲线的差距,然后根据地下水动力学常识或经验增大或减小数值,协调不同参数之间的配比关系,反复进行直到满意为止。这种方法需要经验,而且在判断拟合效果上带有主观性,当参数较多时很难操作。相对而言,自动优化方法比较合理。
最优化方法存在一个判别函数,用来评价模型拟合的效果,可表示为
J?g??1Np,Nk2p?1,k?12???Hp?tk??Hp(tk)?????Hs??1Ni,Nk2i?1,k?12???Qi?tk??Qi?tk???(4.5.24)
??Qs??式中:g为参数组的一个取值状态;H和H?分别为模型计算水头和实测水头;Q和Q?分别为模型计算和实测的地下水径流量;?Hs和?Qs为两个参考值用于数据的归一化。最优化方法的目标是判别函数达到最小值,并把达到目标的参数取值状态作为模型反演的结果
gop??g1,g2,g3,...,gM?:J?gop??min?J?g1?,J?g2?,J?g3?,...,J?gM?? (4.5.25)
式中:gop是优化得到的参数状态,任何一个参数状态gi可表示为
?p1i??p1?min,p1?max???p2i??p2?min,p2?max? (4.5.26) gi??...??pji??pj?min,pj?max????式中:pji为第j个参数在状态gi中的取值;pj?min和pj?max分别为第j个参数的最小和最大允许值;n为待求参数的数目。
如果只有两个待求参数,可以采用穷举法进行优化,即在参数空间中取优选窗口的点阵,每个状态点都运行模型得到一个判别函数的结果。最后,根据判别函数的等值线图确定最优参数状态(图4.21)。穷举法计算工作量较大,如每个参数取值范围都剖分成19个状态点,则需要重复运行模型100次,当模型节点数很大时,计算工作量相当巨大。
为了避免穷举法的缺陷,可以采用其他的选法,这些优选方法能够从一个参数状态快速地搜寻下一个使判别函数更小的参数状态。常用的优化法有逐个修正法、单纯形法和共轭梯度法等,最近几年发展起来的方法有神经网络算法、遗传算法和蚁群算法等。最优化算法的改进方向一般是提高搜索速度、增强稳健性和复杂条件的适应性。然而,当需要反演确定的参数很多时,采用新的最优化算法并不一定能够减少计算工作量。
J(p1,p2) P2 P2-max 3 4 P2-op 1 2 4 P2-min O P1-min P1-op P2-max
图4.21 含两个待求参数的穷取优化法