地下水数值模拟07

第六章

反求水文地质参数的数值方法一、基本概念正演问题与反演问题正演问题（正问题）：

——在已知地下水流动的微分方程及水文地质参数K、M、W、μ、μ*和边界条件的前提下，求解渗流区域内的水头分布规律和流量。反演问题（逆问题）：——根据地下水的天然动态或抽水试验的观测资料研究所选用的方程类型是否适当、确定方程的参数和检验定解条件。一、基本概念反求参数问题的适定性解的存在性——解是否存在？解的唯一性——解是否唯一？解的稳定性——解是否稳定？根据实际资料反求渗流区的水文地质参数(有时包括边界流量)，这样的解是否存在？根据实际资料反求的水文地质参数是否唯一？当实测资料有微小误差时，反求的水文地质参数的误差是否也微小，即水文地质参数是否连续依赖于实测资料。可以满足很难满足很难满足解的唯一性例1均质等厚二维承压含水层稳定流动问题在同一区域同一边界条件下，不管导水系数T取何值，水头分布相同。例2◆正问题：H为未知量，T为已知量可以求得：◆逆问题：T为未知量可以求得：需要补充相应条件以确定参数若已知若单从观测数据反求参数，可能存在多种解！！解的稳定性例水头较小的误差，可能会引起所求参数较大的误差!!若实测水头存在误差ε，即：从数学上讲：虽然ε很小，但其导数可能很大。若数值模型反演问题解的适定性以非均质承压二维非稳定流为例(i,j)(i-1,j)(i,j+1)(i,j-1)(i+1,j)如图所示的有限差分网格，其差分方程为：◆正问题：H为未知量◆逆问题：T、μ*、W为未知量数值模型反演问题解的适定性根据方程个数N和未知参数个数M（m1+m2+m3）之间的关系：（1）N<M:方程组可能存在无穷多组解——不满足唯一性（3）N>M:一般不存在任何一组参数使所有方程同时满足。然而，根据问题本身的物理特性，应该存在一组参数使这些方程基本得到满足（2）N=M:当系数矩阵行列式不等于0，能求出唯一解（需要足够多资料）当系数矩阵行列式的值很小，则水头观测的微小变化可能造成参数的很大变化——不满足稳定性反求参数问题的适定性反求参数问题本身不一定是唯一的，也不一定是稳定的。但是对于实际问题，我们可以根据对水文地质条件的初步认识以及通过其它手段得到一些辅助的资料和参数的约束条件从而使得反求参数问题在一定程度上是唯一的和稳定的。二、反求参数的直接解法求解思路：指在地下水流动微分方程(或描述地下水流动的数值模型)中，将水头值作为已知量，将待求的参数(往往包括源汇项及边界流量)等作为未知量，直接求解未知参数的方法。常用方法：局部直接求逆法数学规划法应用现状：对数据误差十分敏感，对观测资料有过高的要求，因而目前还难以应用。局部直接求逆法例已知不同时刻局部区域水头H和源汇项W的实测值，需要求解T和μ*设：优点：不需要初始条件，也不需要边界条件。对均质各向同性、等厚的渗流区来说，只需要知道某个局部区域在两个不同时刻的水头值和垂直方向的水量交换的实测资料，因此，这个方法称为“局部直接求逆法”。难点：如何求出方程中水头对时间和空间的导数。实际计算发现系数矩阵对应的行列式之值通常很小，因此求参数住往会产生很大的误差数学规划法当方程个数N>未知参数个数M时，方程为超定方程组，一般不存在任何一组参数使所有方程同时满足。

根据问题本身的物理特性，应该存在一组参数使这些方程基本得到满足。设对应于参数组（k1,k2,k3,…,km）,第i个方程存在剩余求一组参数，使得剩余Ri达到最小。可表示为剩余的平方加权和最小故构建最优化问题（规划问题）如下：目标函数约束条件优点：所求参数较好的符合实测资料难点：需要较多的观测资料，对数据误差敏感三、反求参数的间接解法求解思路：先给待定的水文地质参数假设一组初值，通过解正演问题计算相应的水头分布，然后将计算水头值与实测水头值进行对比，看二者拟合程度如何。常用方法：试估——校正法数学规划法应用现状：较为稳定，因而应用较多，但被广泛采用的有限。间接解法的常用公式故解逆问题转化为最优化问题（规划问题）如下：先给待定的水文地质参数｛k｝假设一组初值，通过解正演问题计算出ti时刻的水头分布Hj(ti)，设该点相应的水头观测值为Hjobj(ti)

，则计算值与实测值之间的拟合程度的衡量标准为：非线性约束优化问题误差绝对值误差平方和加权平方和试估——校正法优点：除用正演问题的程序外，不需要其它计算程序。充分发挥解题人员的能动性。难点：当待求参数很多时，反复调整的过程可能延续很长。缺乏一个收敛准则，很难求得最优参数。根据研究区水文地质条件和已有的抽水试验资料初步拟定一组参数值，通过解正演问题计算出各结点各时刻的水头值，然后将计算水头值与实测水头值进行拟合对比，如拟合不好，则对给出的参数初值进行调整，再按正演问题计算。重复这一过程，直到计算水头值与实测水头值之差足够小为止。逐个修正法基本思想：——初步选定一组参数｛k(0)

｝，逐个修正其中每一个分量ki0，全部修正完后便得到一组参数的改进值｛k(1)

｝。特点：以单因素优选法为基础，对参数进行逐个修正，能在满足约束条件下逐步减小目标函数值。但收敛速度不快，只有参数个数不多（n<10）且初值选取得比较好时才能体现出优越性。步骤：怎么找到最优值呢？——单因素优选法例：0.618法、二次插值法给定参数初值：保持其余n-1个参数值不变，对第一个参数按单因素优选法在变化范围内选出参数的改进值，得到参数：保持其余n-1个参数值不变，对第二个参数按单因素优选法在变化范围内选出参数的改进值，得到参数：重复上述步骤，直至全部参数修改一遍，得到参数的改进值检查收敛条件。若满足则停止运算，否则以改进值代替初值，重复第一步0.618法设函数f(x)在区间[a,b]内具有单峰性，即函数在区间上有唯一极小点。若在此区间之内任取两点a1和b1，且a1

＜b1，然后计算这两点的函数值，则可能出现以下三种情况：

怎样取a1、b1？二次插值法假设评价函数E对单个参数而言，可近似的看作抛物线关系，即二次函数。将此抛物线最低点相对应的参数值作为该参数的最优值。取三个不同的k1值，确定抛物线方程中a1，a2单纯形搜索法多维直接搜索法，可同时修正所有的待求参数。所谓单纯形是指在n维空间中具有n+1个顶点的多面体。利用单纯形的顶点，计算其函数值并加以比较，从中确定有利的搜索方向和步长，找到一个较好的点取代单纯形中较差的点，组成新的单纯形来代替原来的单纯形，

例需要确定两个水文地质参数k1，k2取3组不同的（k1,k2），在二维平面上对应三个点pHpLpG单纯形单纯形的基本操作pHpLpG评价函数最大评价函数最小评价函数次大过pH点并穿过其余二点的中点pC的方向为合适的寻找方向pC（1）反射：沿pHpC方向取点pR，使：pHpC=pCpR（2）压缩：若ER≥EH，说明pR前进的太远，需要适当后退则在pH和pR之间另取新点pSpR（3）扩张：若ER<EH，说明沿pHpC方向还可以前进的更远则在pHpR延长线上另取新点pE若EE≤ER，扩张成功，pS=pE若EE>ER，扩张不成功，pS=pR反射压缩扩张（4）收缩：若ES≥EG，说明用pS代替pH不过有多大改善，则将原单纯形缩小，组成新的单纯形n个水文地质参数，给出n+1组。α——扩张因