3.高密度数据反演
在地球物理学中，地球物理反演是利用在地球表面观测到的物理现象推测地球内部介质物理状态的空间变化及物性结构。如果把地球物理问题分为资料采集、数据处理和反演解释三个阶段的话，那么，资料采集是基础，数据处理是手段，反演解释才是地球物理工作的目的。
在高密度电阻率法中，仅根据高密度电阻率法的视参数等值线断面图或分级灰度图来进行解释显然是很不够的。为了获得地下介质电性分布更为精确的图案，必须进行二维电阻率反演。
下面介绍光滑约束最小二乘反演方法
假如我们假定二维电阻率反演所使用的模型由大量电阻率为常数的矩形单元组成（见图4.2）。传统的方法是用非线性最优化迭代方法去确定每一矩形元的电阻率。而利用光滑约束最小二乘反演方法（deGroot-Hedlin and Constable，1990）可以计算每一体元的电阻率，并使得计算的理论值与实测视电阻率值之残差达到最小。最小二乘公式是
?JTJ??CTCp?JTg (4.3)
?式中，J是偏导数雅可比矩阵，?是阻尼因子，g是计算的理论视电阻率和实测视电阻率之差矢量矩阵，p是模型参数的校正矢量矩阵。二维光滑滤波系数c是去保持模型参数在其值连续变化时有一定的光滑度（Sasaki,1992）。要计算的是模型参数的校正矢量p。
反演过程主要分为三步：
（1）计算当前模型下的视电阻率值，一般用有限元法或有限差差分法； （2）计算偏导数雅可比矩阵J； （3）求解线性方程组（4.3）。
假如能避免计算第二步，则可以节约大量的机时。下面我们就讨论这个问题。
4.4.1 均匀介质模型
均匀介质模型是最为理想化的初始模型。为计算模型的偏导数值，1990年McGillivray和Oldenburg详细地讨论了确定其偏导数的数学分析方法，尤其是对地电问题的偏导数，讨论尤为详细。对于均匀地电模型，有可能使用电位的解析解和格林公式去确定其偏导数。对电阻率为?的均匀半空间，其泊松公式
?2???Is??Xs? (4.4)
式中，?是位于点Xs的电流源Is产生的电位。假如给上述公式一个扰动，地下电阻率产生??的变化，则电位将产生??的变化，Park和Van（1991）推导出以下公式
???????2???????d? (4.5)
v这里，假定在每一体元中电阻率的变化为一常数，而其余的地方都为零。参数?'是虚电流源在测量电极处产生的电位。对均匀半空间情况，位于原点（0，0，0）处的电流源在空间任意一点产生的电位为
??和
?Is2?x?y?z?222?1 (4.6)
2????Is2?[(x?a)?y?z]22212 (4.7)
这里，测量电极位于点（a，0，0）。在计算出?、?'的散度后，式（4.5）可写为
??Isx(x?a)?y2?x2??dxdydz (4.8) 2332222222???4?v?x?y?z?[(x?a)?y?z]2当???0时，上式的左边简化为偏导数。积分中的右端项是均匀情况下的Frechet导数。
公式（4.8）对均匀半空间上单极―单极装置在地下一个小体元（x，y，z）上所测得的电位也是适用的。
由公式（4.8）可获得二维矩形元的偏导数，不过要在x、y有限积分后面乘以y方向的无限积分。Baker（1992）建议矩形元的排列方式与拟断面图上的数据点的排列方式一致。当然，靠近地表采用较薄的矩形元，而靠近底部采用较厚的矩形元将获得更好的结果。对于一个矩形元的偏导数可由下式给出（见图4.3）
??Is???4?2
z2x2??z1x1?????x(x?a)?y2?x2?x2?y?z22?32[(x?a)?y?z]2223dxdydz (4.9)
2 图4.3 影响二维面元偏导数计算的矩形元的参数 C1和P1分别是电流电极和电位电极 Figure4.3 The rectangle parameter working on 2-D unit partial derivative 为了简化，令
??Fy????x(x?a)?y2?x2?x2?y2?z2?32[(x?a)2?y2?z2]3dxdydz (4.10)
2这样，公式（4.9）可写为
??Is???4?2z2x2z1x1??Fydxdz (4.11)
公式（4.10）中的积分可用解析方法解出。方法如下：
把（4.10）改写为
Fy???0?x?02?y?z222212??x?x?a??y2xa23222?z2232?dy (4.12)
??令
?x?y?z??x?x?a??y?z232?dy?2?x2?z2，?2?(x?a)2?z2
则（4.12）式变为
Fy???0?x2?y?z?02212??x?x?a??y23222?z232?dy (4.13)
?2xa???2?y???12?y232?dy上述积分结果可写为（Gradshteyn and Ryzhik，1965）
2??2E(k)??2K(k)xa[?2??2E(k)?2?2K?k??Fy??? (4.14) 2?22222??????????????????式中，K(k)和E(k)分别是第一类和第二类完全椭圆积分（Press et al.，1988）。由于公式（4.13）的解要求?>?，因此上式仅当x?0.5a时才有意义。 当x?0.5a时，需要按下式重写公式（4.12）
Fy???0?x?02?y?z222232??x?x?a??y2a?x?a?23222?z2212?dydy
(4.15)
??令
?x?y?z??x?x?a??y?z232? ?2?(x?a)2?z2，?2?x2?z2
同样可得x<0.5a时的表达式
2??2E(k)??2K(k)?x?a?a[?2??2E(k)?2?2K?k??Fy??? (4.16) 2?22222??????????????????特殊情况下，当x=0.5a时，令
?2?(0.5a)2?z2
则有
Fy???20??2?y22?dy???a20??2?y23?dy (4.17)
上述积分结果可写为（Gradshteyn and Ryzhik，1965）
?13a2?Fy???3? (4.18) 5?16???2?总结如下：
（1）当x>0.5a时，Fy为
2??2E(k)??2K(k)xa[?2??2E(k)?2?2K?k???? (4.19) Fy??222222??????????????????式中
?2?x2?z2,?2?(x?a)2?z2k?(???)220.5a
????0其中，K(k)和E(k)分别是第一类和第二类完全椭圆积分（Press等，1988）。 （2）当x<0.5a时，Fy为
2??2E(k)??2K(k)xa[?2??2E(k)?2?2K?k???? (4.20) Fy??222222??????????????????这里
?2??x?a?2?z2,?2?x2?z2
（3）当x=0.5a时，Fy为
?13a2?Fy???3? (4.21) 5?16a??2a