地下水数值模拟基础：计算矩阵[A]的方法

修正的Gauss-Newton法需要计算偏导数=1,2,…,L;i=1,2,…,n)。根据f l的表达式,上述偏导数的计算最终必然归结为计算,此处的H是用参数k模拟计算得到的水头在选定的观测点和观测时刻的值,它除了是坐标、时间的函数外,还是所选参数k的函数,即对二维问题有H=H(x,y,t,k)。这个偏导数在解逆问题中起着重要作用,有时把它称为灵敏度系数(W.Yeh,1986)。

计算它的简单方法是选用差分格式,如中心式差分格式来近似,为此有

k i的增量Δk i的取值既不能太大,以减少由式(4-42)近似计算引起的截断误差;也不能太小,以便能观测到H值有一定的变化。通常和k i成一定比例,如Bar d(1974)建议的那样,Δk i=ak i,式中比例系数a在下列范围内取值,10-5≤a≤10-2。对于实际问题一般通过试算来确定。这种方法有时称为影响系数法(W.Yeh,1986)。

另一种称为灵敏度方程法(W.Yeh,1986),它是和模拟计算结合起来一起进行的。兹以下列二维水流方程为例加以说明:

式(4-43)对k i求偏导数,得

相应的初始条件和边界条件为

式中式(4-43)和式(4-44)左端在形式上是完全相同的,只是未知变量不同而已;右端两者也是相似的,式(4-43)中的w相当于式(4-44)中方括号内的项。因此,两个方程可以用相同的程序来求解。(www.xing528.com)

基于这样的认识,如用有限元法求解,与式(4-43)对应、离散得到的代数方程组为

式中:矩阵[D]和[P]的元素分别为

将式(4-45)对k i求导可以得到

上式左端圆括号内的项显然与式(4-45)中{F}相当。其中矩阵[D ki]和[P ki]的元素分别为

为了计算被积函数中的,首先要将T和S表示为k 1,k 2,…,k n的函数,事实上已经这样做了,所以能很方便地利用数值积分求出它对k i的导数来,于是矩阵[D ki]和[P ki]也就形成了。具体运算时,先解式(4-45)求出{H}。有了{H},利用初始条件和差分公式可以方便地求得。这样式(4-46)左端圆括号内各项就确定了。据此运用相同的解法,不难从式(4-46)解得{H ki},即我们需要的偏导数为了求得和参数组k 1,k 2,…,k n对应的水头H及其偏导数需要同时解n+1个有相同系数矩阵的方程组。

需要指出的是由组成的矩阵在利用修正的Gauss-Newton法反求水文地质参数的过程中有重要意义,如果计算不精确,则识别阶段所得结果也是不可靠的。