高密度电法-专著

高密度电阻率法仍然是以岩、土导电性的差异为基础, 研究人工施加稳定电流场的作用下地中传导电流分布规律的一种电探方法。因此, 它的理论基础与常规电阻率法相同,所不同的是方法技术。高密度电阻率法野外测量时只需将全部电极( 几十至上百根) 置于观测剖面的各测点上, 然后利用程控电极转换装置和微机工程电测仪便可实现数据的快速和自动采集, 当将测量结果送入微机后, 还可对数据进行处理并给出关于地电断面分布的各种图示结果。显然, 高密度电阻率勘探技术的运用与发展, 使电法勘探的智能化程度大大向前迈进了一步。由于高密度电阻率法的上述特点, 相对于常规电阻率法而言, 它具有以下特点: ( 1 ) 电极布设是一次完成的, 这不仅减少了因电极设置而引起的故障和干扰, 而且为野外数据的快速和自动测量奠定了基础。 ( 2 ) 能有效地进行多种电极排列方式的扫描测量, 因而可以获得较丰富的关于地电断面结构特征的地质信息。

(3) 野外数据采集实现了自动化或半自动化, 不仅采集速度快( 大约每一测点需2～5s) ,而且避免了由于手工操作所出现的错误。

(4)可以对资料进行预处理并显示剖面曲线形态, 脱机处理后还可自动绘制和打印各种成果图件。

(5)与传统的电阻率法相比, 成本低, 效率高, 信息丰富, 解释方便。 （一）高密度电阻率法采集系统

早先的高密度电阻率法采集系统采用集中式电极转换方式。如图4.1所示。进行现场测量时，用多芯电缆将各个电极连接到程控式电极转换箱上。电极转换箱是一种由微片机控制的电极自动转换装置，它可以根据需要自动进行电极装置形式、极距及测点的转换。电极转换箱开关由电测仪控制，电信号由电极转换箱送入电测仪，并将测量结果依次存入存储器。

图4.1 高密度电阻率法测量系统结构示意图（集中式）

随着技术的发展，高密度电法仪日趋成熟。表现在：采用嵌入式工控机，大大提高系统的稳定性与可靠性；采用笔记本硬盘存储数据，可以满足野外长时间施工的工作需求；系统采用视窗化、嵌入式实时控制与处理软件，便于野外操作；可实现多种工作模式的转换，计算机与电测仪一体化，携带方便。新一代高密度电法仪多采用分布式设计。所谓分布式是相对于集中式而言的，是指将电极转换功能放在电极上。分布式智能电极器串联在多芯电缆上，地址随机分配，在任何位置都可以测量；实现滚动测量和多道、长剖面的连续测量。

图4.1 高密度电阻率法测量系统结构示意图（分布式）

系统可以做高密度电阻率测量，又可以同时做高密度极化率测量，应用范围宽。

图

1.常用装置

高密度电阻率法在一条剖面上布置一系列电极时可组合出十多种装置。高密度电阻率法的电极排列原则上可采用二极方式，即当依次对某一电极供电时，同时利用其余全部电极依次进行电位测量，然后将测量结果按需要转换成相应的电极方式。但对于目前单通道电测仪来讲，这样测量所费时间较长。其次，当测量

电极逐渐远离供电电极时，电位测量幅值变化较大，需要不断改变电源，不利于自动测量方式的实现。高密度电阻率法常用的装置见图，包括温纳装置（Wennerα、Wennerβ、Wennerγ）、偶极-偶极装置(Dipole- Dipole)、三极装置(Pole-Dipole、Dipole- Pole)、温纳-斯伦贝谢装置(Wenner-Schlumberger)等。

图

(a)温纳α装置；(b)温纳β装置；(c)温纳γ装置；(d)偶极-偶极装置

(e)三极装置；(f)温纳-斯伦贝谢装置

2.装置特点及视参数的计算 （1）温纳装置

在高密度电阻率法中，由于温纳装置与异常对应关系好，是常用的装置之一。最早的高密度电阻率法一般使用三电位电极系。所谓三电位电极系就是将温纳装置、偶极装置和微分装置按一定方式组合后构成的一种测量系统。这是由于电极转换需要时间，因此当连接好等距的AMNB四个电极后，可以作三次组合，依次构成温纳装置、偶极装置和微分装置，或称为温纳α装置、 温纳β装置和温纳γ装置。这样在某一测点就可以获得三个电极排列的测量参数。

温纳装置对电阻率的垂向变化比较敏感，一般用来探测水平目标体。温纳装置的装置系数是2a，相比于其它装置而言是最小的。因而同样情况下，可观测到较强的信号，可以在地质噪声较大的地方使用。另一方面，由于它的装置系数小，因此在同样电极布置情况下，它的探测深度也小。另外，温纳装置的边界损失较大。

温纳α装置、温纳β装置和温纳γ装置三种排列形式（见图），视电阻率参数及计算公式为：

Usk，k2a

IUsk，k6a

IUsk，k3a

I根据三种电极排列的电场分布，三者之间的视电阻率关系：

12sss

33对高密度电阻率法而言，由于一条剖面地表电极总数是固定的，因此，当极距扩大时，反映不同勘探深度的测点数将依次减少。图显示了温纳装置测点分布。

图温纳α装置测点分布示意图

x：最小电极距；n：间隔系数

由图可见，剖面上的测点数随剖面号增加而减小，断面上测点呈倒梯形分布，任意剖面上测点数可由下式确定：

Dn =Psum-(Pa-1)·n

式中：n—间隔系数；Dn—剖面上测点数；Psum—实接电极数；Pa—装置电极数，对三电位电极系而言，Pa =4；对三极装置，Pa =3。

如对温纳装置而言，设有30路电极，则Dn303n。当n=1时，第一条剖面上的测点数D157。令Dn1，可求出最大间隔系数为nmax9。总测点数剖面数而言，总测点数N为：

N(303n)

n19（2）偶极-偶极装置

偶极-偶极装置高灵敏度区域出现在发射偶极和接收偶极下方，这意味着本装置对每对偶极下方电阻率变化的分辨能力是比较好的。同时，灵敏度等值线几乎垂直的，因此偶极-偶极装置水平分辨率比较好，一般用来探测向下有一定延伸的目标体。相对于温纳装置，偶极-偶极装置观测的信号要小一些。

skUMN，k2n(n1)(n2)a I（3）三极装置

三极装置有更高的灵敏度和分辨率（Sasaki 1992）。同时，三极装置的两个电位电极在网格内，因此受电噪声干扰也相对小一些。与偶极-偶极装置相比，三极装置所测信号要强一些。另外，三极装置可以进行“正向”（单极-偶极）和“反向”（偶极–单极）测量，因此边界损失小。

UMN，k2n(n1)a skI（4）温纳-斯伦贝谢装置

温纳-施伦贝谢是一个变种（Pazdirek and Blaha，1996），其高灵敏度值出现在测量电极之间的正下方，有适当的水平和纵向分辨率，但探测深度小，在三维电法难以单一使用。

UMN，kn(n1)a skI可以联合使用这些装置，有的程序可联合反演。 （二）资料处理与反演解释 1.统计处理

统计处理包括以下内容：

（1）利用滑动平均计算视电阻率的有效值，例如三点平均：

x(i)(s(i1)s(i)s(i1)3

式中，i=1,2,3，……，Dn，x(i)为i点的视电阻率有效值。 （2）计算整个测区或某一断面的统计参数

1N平均值：xx(i)，N为某一测区或某一断面上的测点数

Ni1N22/n，标准差：As(i)nxAi1[s(i)x]2/n相等马？

i1N（3）计算电极调整系数

K(L)x(i)/s(L)

其中s(L)为电极距为L时全部视电阻率观测数据平均值。 （4）计算相对电阻率

y(i)K(L)x(i)xx(i)/s(L)

通过计算相对电阻率，可以在一定程度上消除地点断面由上到下水平地层的对变化。因此，相对电阻率断面图照顾要反映地电体沿剖面的横向变化。 （5）对视参数分级

为了对视参数进行分级，首先必须按平均值和标准差关系视参数的分级间隔。间隔太小，等级过密，间隔太大，等级过稀，都不利于反映地电体的分布。一般情况下，以采用五级制为宜，即根据平均值和标准差的关系划分四个界限：

D1xA；D2xA/3；D3xA/3；D4xA

利用上述视参数的分级间隔，可将断面上各点的s(i)或y(i)划分成不同的等级用不同的符号或灰阶（灰度）表示时,便得到视参数异常灰度图，如：

s(i)D1，低阻

s(i)D1~D2，较低阻

ρs(i)D2~D3，中等

s(i)D3~D4，较高阻

ρs(i)D4，高阻

视参数的等级断面图在一定条件下能比较直观和形象反映地点面的分布特征。

统计处理原则上适应于三电位电极系中各种电极排列的测量结果，只是在考虑视电阻率参数图示时，由于偶极和微分两种排列的异常和地电体之间具有复杂得对应关系，因此一般只对温纳装置的测量结果进行统计处理。当然，随着现代高密度电法仪装置的增加，温纳-斯伦贝谢装置的测量结果也可进行统计处理。 2. 比值参数

高密度电阻率法的野外观测结果出了可以绘制相应装置的视参数断面图外，根据需要还可绘制两种比值参数图。考虑到三电位电极系中三种视参数异常的分布规律，选择了温纳β装置和温纳γ装置两种装置的测量结果为基础的一类比值参数。该比值参数的计算公式为：

T(i)s(i)s(i)

由于温纳β和温纳γ这两种装置在同一地电体上锁获得的视参数总是具有相反的变化规律，因此用该参数绘制的比值断面图，在反映地电结构的分布形态方面，远比相应装置的视电阻率断面图清晰和明确的多。

图？？是对所谓地下石林模型的正演模拟结果。模型的电性分布已如图所示，其中温纳α装置的s拟断面图几乎没有反映，而T比值断面图则清楚地反映了上述模型的电性分布。

另一类比值参数是利用联合三极装置的测量结果为基础组合而成的，其表达式为：

sA(i)sB(i) (i,i1)ABs(i1)s(i1)式中s(i)和s(i1)分别表示剖面上相邻两点视电阻率值，计算结果示于i和比值参数反映了联合三极装置歧离带曲线沿剖面水二乘向的变化i1点之间。

率。图表征比值参数在反映地电结构能力方面所作的模拟实验，视电阻率s断面图只反映了基底的起伏变化，而比值断面图却同时反映了基底起伏中的低阻构造。

2 高密度电阻率法二维地形边界元数值解法

高密度电阻率法是常规电法的一个变种，就其原理而言，与常规电法完全相同。它仍然以岩、矿石的电性差异为基础，通过观测和研究人工建立的地中稳定电流分布规律，解决水文、环境与工程地质问题。高密度电阻率法的正演问题就是传导类电法的正演问题，也就是求解稳恒点电源电流场的边值问题。

对二维地形，设起伏地面下均匀各向同性介质的电阻率为1，具有电流强度为I的稳恒点电流源位于地面任一点A(x,0,z)。域Ω的边界由1和2组成（见图1）。

根据位场理论可知，在有源域内及其边界任意一点M(x,y,z)处的电位U (x,y,z)满足：

控制微分方程 2U(x,y,z)F M∈Ω (1) 自然边界条件 Q(x,y,z)UQ(x,y,z) M∈1 (2) n本质边界条件 U(x,y,z)U(x,y,z) M∈2 (3)

式中，F2I1(MA)，Q和U分别是边界1和2上已知边值函数，这里，

Q0，UI11，rAi为源点A到场点“i”间的距离。 2rAi 可以看到，地形是二维的，即沿y轴无变化，而电位U是y偶函数，所以我们也把上边值问题称为2.5维问题。对（1）、（2）、（3）式进行余弦傅立叶变换可得：

控制微分方程 2u(x,,z)2u(x,,z)f M∈Ω (4) 自然边界条件 q(x,,z)q(x,,z) M∈1 (5)

本质边界条件 u(x,,z)u(x,,z) M∈2（6） 式中，fI1(xxA,zzA)，q0，uI1K0(rA),K0(rA)为第二类零2阶修正贝塞尔函数，是余弦傅立叶变换量或波数。

这样，便将三维偏微分方程（1）变成了而二维偏微分方程（4），即将三维空间的电位U(x,y,z)变换为二维空间的变换电位u(x,,z)。为求得u，可采用边界元法求解u所满足的亥姆霍兹方程（4）式，借助格林公式及二维介质亥姆霍兹方程的基本解，即可把u满足的亥姆霍兹方程及边界条件等价地归化为如下的边界积分方程：

iIui1K0rAuq*dqu*d （7）

2411式中，i为边界点“i”对区域的张角，u为亥姆霍兹方程的基本解，

u*k1*K1(kr)cos(n,r),k为波数，r为点(xi,zi)到点uK0(r),qn22(x,z)的矢径，n为边界的外法线方法，K1(kr)为第二类一阶修正贝塞尔函数。 采用边界元离散技术，将域Ω的边界1进行剖分，分成N1个单元。根据积分的可加性，（7）式中对边界1的积分可化为对每个单元j上的积分之和：

ciuiBij1N1juqd*j1N1jqu*d （8）

式中ciiI，Bi1K0rA。 24i 方程组（8）式仅是含有N1未知量的线性方程组，解此方程组即可求得变换电位值u(x,,z)，然后按下式：

2 U(x,y,z)u(x,,z)cos(y)d

0进行傅氏逆变换，即可求得电位值U(x,y,z)。

根据所采用的高密度电阻率法装置类型，逐点计算出某记录点处的纯地形异常视电阻率值sD，然后用“比较法”进行地形改正。地形改正公式如下：

sGs/sD/1

式中sG地形改正后的视电阻率值；s该记录点实测的视电阻率值；sD纯地形影响值，它是一个无量纲的标量，1一般取1m。

利用边界单元法计算高密度电阻率法地形边界位场问题是很有效的。但是在算法引入时，必须针对高密度电阻率法的特点，作一些技术处理。高密度电阻率法电极排列密集，并且采用了差分装置，所有这些特点都要求计算精度高，运算速度快。另外，所形成的矩阵也因测量电极到供电电极的距离变化很大呈带状分布，并且当波数较大时，矩阵中的系数几乎都接近于零，造成解的不稳定。为了解决这一问题，采用了增广矩阵法求解方程组HUB的效果较为满意。

为了保证精度，同时又减少运算次数，除采用九波数傅氏反变换外，还采用

了不等分单元剖分方案。具体做法是，在测线外，越远则单元剖分长度越大，且为最小电极距的整数陪；在测线段，则以最小电极距长度划分边界单元。为了避免r等于零时贝塞尔函数无穷大的问题，剖分结点应不与电极点位置重合，最好选取相邻电极的中点为结点。

2.1二维地电构造中点电流源场的正演有限元算法的基本原理

[4]

有限单元的思想最早是由Courant于1943年提出（倪光正）。20世纪50年代初期，由于工程分析的需要，有限元法在复杂的航空结构分析中最先得到应用，而有限元法(Finite Element Method)这个名称则由Clough于1960年在其著作中首先提出。半个世纪以来，以变分原理为基础建立起来的有限元法，因其理论依据的普遍性，不仅广泛地应用于各类结构工程，而且作为一种声誉很高的数值计算方法已被普遍推广并成功地用来解决其他工程领域中的问题，如热传导、渗流、流体力学、空气动力学、岩土力学、机械零件强度分析、电气工程问题等。

20世纪70年代初，J.H.Coggon首先将有限元法用于电法勘探，后来L.R.Rijo完善了有限元法数值计算方法，使之成为正演模拟计算的有效方法（罗延钟）。80年代初，我国的周熙襄等在引进Coggon-Rijo的算法时，将Dey和Morrison用于有限差分模拟的混合边值条件引入了有限单元法，从而发展了Coggon-Rijo的有限单元算法。此后，罗延钟等在选用边值条件和反傅氏变换的算法及波数取值等方面，又有了一些新的发展，使整个算法更臻完善，能在不做任何校正的情况下，对相当大范围内变化的电极距，取得较高精度的计算结果。90年代中期，杨进提出了迭代有限元算法，它不仅模拟复杂地球物理模型的能力强、模拟精度高、且占用计算机内存小，是对有限元方法的又一大发展。

以前的二维地电构造中点电流源场的有限元算法都是建立在常规直流电法基础之上的，而高密度电阻率法因其装置的特殊设计，所以用有限元进行正演计算时必须兼顾速度和精度。本文详细地研究了二维地电构造中点电流源场的正演有限元法的理论基础，并结合高密度电阻率法的特点，在有限元算法本身及其相关技术上作了许多改进，使之成为高密度电阻率法正演计算的有效算法。

2.1.1稳定电流场微分方程边值问题及其相应的变分问题

高密度电阻率法二维正演计算问题即点源场二维地电断面的边值问题。取三维笛卡儿坐标系的X轴为垂直于地质体走向的方向，Y为轴为地质体走向的方向，Z轴垂直指向下。我们知道，二维情况下一般设地下导电性沿Y轴方向无变化，即(x,z)，点电流源

A(xA,ZA)位于XOZ平面上。这时，虽然地下介

图2.1 点源二维区域示意图 Figure2.1 Sketch map of point-element for 2-D area 质的电导率是二维的，而场源是三维的，电位函数u(x,y,z)所满足的微分方程及边界条件为：

①

u(x,y,z)u(x,y,z)u(x,y,z)((x,z))((x,z))((x,z))f xxyyzz（2.1）

式中fI(xxA)(y)(zzA)，I为A点的电流强度。

② 在地面边界1上应满足第二类边界条件，即诺依曼条件，在电法勘探中也可称为绝缘边值条件

un0 （2.2）

1③在其余边界2上，例如在离电流源很远的地方，可设

u0 (2.3)

2这种给定已知电位的边界条件称为第一类边界条件，即狄里希莱条件，或称极限边界条件，也可以称作强加边界条件。 ④第三类边值条件

实际计算时，计算网格毕竟是有限的，式（2.3）中的无穷远边界条件很难满足，因此提出了第三类边值条件。对均匀半空间点源二维问题，地中任一点电位的变化总有下面的一般形式

u(x,y,z)cx2y2z2c (2.4) r式中c为常数，并设电源点放在坐标原点，r为点源到计算点的径向距离，因此有

ucu2rncos (2.5) nrr式中为r和边界外法线n之间的夹角，上式也可以写为

uu0 (2.6a) nr式中cos,上式便是第三类边值条件，或曰混合边界条件。

但是对于非均匀地电断面，便导不出（2.6a）式。为此，需采用更一般形式的混合边界条件，即

uu0 （2.6b） nr为与地电断面不均匀性等有关的一个系数。式中，大量的试算发现取值在0～

1之间，均值约为1/2。关于系数的取值问题暂不作讨论，在下面的公式推导中，为记述的方便，不妨暂时取=1。

⑤在所论由边界1和2包围的区域内，两种具有导电率为1和2的介质的交界面处，电位和电流密度法向分量应满足以下衔接条件,或曰普遍边界条件 (a) 由于电位的连续性，在交界面处，有

u1(x,y,z)u2(x,y,z) (2.7)

(b) 由于电流法向分量的连续性，在交界面处，有

1u1u22 (2.8) nn其中n为交界面的法线方向的单位矢量。

前已谈到，对于点源二维问题，介质的电阻率是二维分布，而场源是三维的。为使问题简化，我们常把三维问题通过傅氏变换化为二维问题。由于电位是关于y轴的偶函数，故取余弦变换对

2F()0f(x)cos(x)dx (2.9) 2f(x)F()cos(x)d0对（2.1）式按（2.9）式取傅氏变换，得

v(,x,z)v(,x,z)((x,z))((x,z))2(x,z)v(,x,z)f'(2.10) xxzz1式中f'I(xxA)(zzA)，v(,x,z)称为傅氏电位、变换电位或傅氏变

2换电位。

这样，便将三维微分方程（2.1）式变成了二维微分方程（2.10）式。当分别对若干个给定的波数值求解方程（2.10）式，计算出傅氏电位v(λ，x，y)后，再按（2.9）式作傅氏反变换，即可计算出所求的电位 u(x,y,z)2v(,x,z)cos(y)dy （2.11）

0在求解偏微分方程（2.10）时，于求解边界上，宜采用如下边值条件： ①在地面边界1上，应用绝缘边值条件，即第二类边值条件：

vn0 （2.12）

1②在其余边界2上，采用前述第三类边界条件，即混合边界条件。其导出过程如下

由于求解区的左、右和底边界远离场源和电性异常体，故在那里的电场分布可近似看成与均匀大地情况相同。设场源是位于Ak(xk,o,zk)，电流强度为Ik的

N个点电流源(k=1,2,3……N)，则对于上述边界上的某点P(x,y,z)，其电位可写成如下形式

u(x,y,z)k1NQIkry2k2 (2.13)

式中Q，—均匀大地的电阻率；rk(xxk)2(zzk)2。 2对（2.13）式作傅氏变换，得

v(,x,z)u(x,y,z)cos(y)dy

0 0k1NNQIkry2k2cos(y)dy

QIkK0(rk) (2.14)

k1式中的K0(rk) 是零阶修正贝塞尔函数。对（2.14）式沿界外法线n求方向导数，并考虑到修正贝塞尔函数的求导性质（dK0(x)/dxK1(x)），可得

Nv(,x,z)Q.IkK1(rk)cosk (2.15) nk1式中,K1(rk) 为一阶修正贝塞尔函数；k为二维断面（x—z坐标面）上，场源

Ak到所论点P的矢径rk与P点边界外法线方向n之间的夹角。

从（2.14）和（2.15）式中消去常数Q，可得上述边界上的混合边值条件

vv0 (2.16) n2式中，Ik1NknK1(rk)cosk 。

kIk1K0(rk) 于是，所提出的问题可归结为对若干个给定的波数λ值，求解电位的傅氏变换电位v所满足的下列二维偏微分方程的边值问题 ①v满足二维偏微分方程

vv((x,z))((x,z))2(x,z)vf (2.17) xxzz其中

1NfIk(xxk)(zzk) （2.18）

2k1②在地面边界1上

v0 (2.19) n1③在其余边界2上

vv0 (2.20a) n2其中 Ik1NknK1(rk)cosk (1.20b)

kIk1K0(rk)与上述二维偏微分方程边值问题（2.17式、2.19式、2.20a式）等价的变分问题为

vvJ(v)()2()22v22fvdsv2dl极值 (2.21)

xz2从二维条件下的欧拉方程边值问题出发，我们可以证明（2.17、2.19、2.20a）与（2.21）式的等价性。并且地面边值条件(2.2)式在泛函极值问题（2.21）式中没有出现，这是因为这一条件在泛函极值问题中能自然得到满足，故称此边界条件为变分问题的自然边界条件。此外，在求解区（由边界12围成的区域）内存在有限个电导率突变面时的内边界条件（2.7）、（2.8）式，在变分问题中也是自然满足的，不需另作处理。至此，用有限单元法计算二维地电构造中的点电流源场的核心，就是用数值方法解变分问题（2.21）式。

求解变分方程（2.21）式，就是要找出一个傅氏电位的空间坐标函数

v(,x,z)，以使泛函J(v)最小。有限单元法就是用来求解这一问题的数值解的计

算技术，它依据泛函欧拉方程边值问题与泛函极值问题的等价性，将微分方程的边值问题转化为相应的泛函极值问题。对于稳定电流场，根据电场总能量最小化原理，泛函的极值问题（2.21式）应满足

v2v2222J(v)()()v2fvdsvdl0 （2.22）

xz2即泛函取得极值的必要条件是它的变分为零。如果将整个求解区间Ω划分为若干个单元e（参见图2.2），e足够小，以致可以认为傅氏电位v(,x,z)在各单元内呈简单函数形式（比如，在二维条件下，可假设v(,x,z)在各三角单元内呈线性变化），于是，泛函的极值问题（或变分问题）又可简化为多元函数的极值问题，而多元函数的极值问题是大家所熟知的。这就是有限单元法的基本思想。

求出傅氏电位v(,x,z)后，利用傅氏反变换（2.11式），即可求出电位函数

u(x,y,z)。对于高密度电阻率法来讲，一般y0，即要求求出主剖面上的电位

图2.2 区域剖分示意图 u(x,z)，再利用公式

sKUMN （2.23） I即可计算出各种装置的视电阻率值。

有限单元法的优点是对形状不规则的地质体和地形起伏模拟能力强，但它的计算程序编制相当复杂。要编制出一个效度好、适应能力强的高密度电阻率法正演有限元算法程序，需要做很多理论研究和实际工作才行。

2.1.2 变分问题的离散化

为了用有限单元法求解偏微分方程的边值问题，首先需要将其等价的变分问题离散化。也就是要将（2.21）式离散化，这包括求解区的离散化（网格剖分）和泛函J(v)的离散化（导出线性方程）两个组成部分。 ① 网格剖分

由于受计算机的计算量、储存量及计算速度的限制，所有的数值计算方法都需要将无限大地中的电场分布限定在有限的求解区内进行讨论，并将连续的求解区域离散化，即网格剖分。

原则上讲，网格剖分可以根据所研究的地电条件灵活地将求解区限定为任何形状，并灵活采用任意适当的网格来剖分求解区。在二维情况下，最常用的办法是将求解区剖分为一系列互不重迭的三角单元，每个单元的顶点称为节点。剖分的规则是：1)整个剖分范围应该是越大越好，剖分范围越大，计算精度也就越高。但剖分得越大，计算量增大，计算速度降低。因此，在作网格剖分时，既要考虑精度，又要照顾速度；2)各单元的节点只能与相邻单元的节点相重，而不能成为相邻单元的边内点；3)要使每个单元内介质的电导率为常数，在不同电导率介质的内分界面上或整个求解区的外边界上，使三角形单元的一个边相互衔接，以这样的折线去逼近边界线；4)三角形单元最好是接近正三角形，不要使其中一个角很小或出现很大的钝角；5)在电场变化剧烈，电位参数的二阶或高阶导数大的地方（如：在异常体边沿、场源附近等），单元宜剖分的细一些，而在电场变化平缓的地方（如：远离异常体和场源），单元面积可取得大一些。当然，单元面积的变化应是渐变的；。

为使程序简化，L.R.Rijo采用了图2.3所示比较规则的矩形求解区和三角形剖分网格。试算结果表明，这种在矩形网格中布置交叉对称三角形剖分网格可以足够近似地模拟一般常见的不平地形和电性异常体，同时又能节省计算量，是经典的一种剖分方法。

在完成上述单元剖分的同时，于求解区内形成了若干个节点（包括内节点和边节点），待定的傅氏电位函数

v(,x,z)在这些节点上的值，

称为节点函数值，自然它们是待定的未知值。这样，便用求解区有限个待定节点函数值来近似表征待定函数在连续空间中的分布。这些待定的节点函数值v1，v2，…，vM（M—总结点数），组成一组独立变量。所谓求待定函数的数值解，就是要确定这些节点的函数值。

②泛函的离散化 1）线性插值

为了计算二维变分问题（2.21）式中积分形式的泛函J(v)，需要知道整个求解区Ω内的v函数值。前面说到，求解区经过网格剖分后，可以用求解区内有限个待定节点函数值来近似表征待定函数在连续空间中的分布。这样，对于二维情况，就可以利用各节点的函数值在各单元内作线性插值来逐个单元计算v，也就是说，将v表示成该单元三个节点之函数值的函数，具体算法如下

设第e个单元的三个节点按逆时针方向的编号分别为i、j、m，其坐标记为

(xi,zi)、(xj,zj)、(xm,zm)，对应的节点函数值为vi、vj、vm（参见图2.2）。

图2.3 L.R.Rijo网格剖分示意图 假设各单元内，函数v(,x,z)是线性变化的，即

v(x,z)abxcz （因为积分变量是一个常数，不妨把 v(,x,z)记为v(x,z)）

（2.24）

式中的的系数a、b、c可由单元上的三个节点的函数值和坐标算出。对于单元e，三个节点上都满足（2.24）式，故有

viabxiczivjabxjczj （2.25） vmabxmczm按克莱姆(Cramer)法则求解上述线性方程组，得

vixixjxmzizm1xizizj zm avjvmzj1xj1xm(aiviajvjamvm)/2 （2.26a）

同样，有

b(bivibjvjbmvm)/2 （2.26b） c(civicjvjcmvm)/2 （2.26c）

式中

aixjzmxmzjajxmzixizmamxizjxjzibizjzmbjzmzibmzizjcixmxjcjxixm （2.27） cmxjxi单元面积

11xj21xm1xizizjzm1(bicjbjci) （2.28） 2 将（2.26a、2.26b、2.26c）式代入（2.24）式，便得到单元e内函数v线性插值的近似表示式

v(x,z)1aibixcizviajbjxcjzvjambmxcmzvm2

（2.29） 或

v(x,z)Ni(x,z)viNj(x,z)vjNm(x,z)vm （2.30）

式中的Ni、Nj、Nm称为基函数，分别为

Ni(x,z)(aibixciz)/2Nj(x,z)(ajbjxcjz)/2 （2.31） Nm(x,z)(ambmxcmz)/2(x,z)e

如果在每个单元上都作出了v(x,z)的这种近似，就能得到整个求解区内

v(x,z)的总体近似函数，这个函数在各个单元内是线性的，即它是一个分片为线

性的函数。对任意两个相邻单元来说，近似函数在公共边上的值，被两个端节点的函数值唯一决定，所以，总体近似函数在整个求解区自然是连续的。

2）单元分析

至此，二维变分问题（2.21）式中的泛函——整个求解区Ω上的积分J(v)，可以分解为各个单元e上的积分之和，即

J(v)Je(v) （2.32）

e而Je(v)不难根据各单元e内函数v的线性插值近似表示式（2.29）或（2.30）算出。因此，首先应从分析一个三角单元e出发。对于当前问题的变分提法（2.21式），三个顶点分别为i，j，m的三角形单元e（参见图2.2）上的泛函

vvJe()2()22v22fvdsv2dl

xze2vv=e()2()2ds2ev2ds

xzee-vIk(xxk)(zzk)dsv2dl (2.33)

ek12N式中，e表示单元e内的电导率值，由于它在每一个单元是常数，因此e可以提到积分外。同样，积分变量也可以提到积分外。

根据（2.24）式及（2.26b）、（2.26c），有

vb(bivibjvjbmvm)/2 （2.34a） xvc(civicjvjcmvm)/2 （2.34b） z可见，它们只与节点坐标及节点函数值有关，在一个单元内是常数，故可提到积分号外。于是，可将（2.33）式中第一项积分表示为

v2v2(x)(z)ds

evv=()2()2dxdz

zex1(bivibjvjbmvm)2(civicjvjcmvm)2 (2.35) 4同时，根据狄拉克函数的性质，有

(xxA)(zzA)dxdze1Ae (2.36)

0Ae于是（2.33）式中积分的第三项为

IiviIjvjImvmvIk(xxk)(zzk)ds (2.37) ninjnmk1eN式中，Ii、Ij、Im分别表示在i、j、m点的供电电流强度。若场源（供电电极）不在某个节点上，则该点上的电流值为零；ni、nj、nm分别表示以节点i、j、

m为顶点的单元数。（2.33）式的第二项积分

2vdsNi(x,y)viNj(x,y)vjNm(x,y)vmdxdz ee2 vi2222N(x,z)dxdzv(N(x,z))dxdz jjiee22vm(Nm(x,z))dxdz2vivjNi(x,z)Nj(x,z)dxdzee2vjvmNj(x,z).Nm(x,z)dxdze(2.

2vmviNm(x,z)Ni(x,z)dxdze38 )

为了计算上式中的积分，我们首先来看基函数Ni(x,z)、

Nj(x,z)、Nm(x,z)的几何意

义。对于三角单元e（i,j,m）中的任一点P(x,y)（见图2.4），由（2.27）式可算得

图2.4 基函数几何意义示图 aibixcizxjxmzjzmzzjzm1zj1zmx1xj1xmz1x1xj1xm2pjm

式中，pjm为三角形pjm的面积。故由（2.31）式 得

Ni(x,z)aibixcizpjm （2.39）

2即Ni(x,z)表示以P(x,z) 为一顶点，jm为对边的三角形面积pjm 与单元面积

之比值。Nj(x,z) 和Nm(x,z) 的几何意义可以类推。由此可得Ni(x,z)、

Nj(x,z) 和Nm(x,z) 的下述性质

Ni(x,z)Nj(x,z)Nm(x,z)1 （2.40）

Ni(xi,zi)1Ni(xj,zj)0 (1-41) Ni(xm,zm)0以及

Nj(xi,zi)0Nj(xj,zj)1 (1-42) Nj(xm,zm)0另一方面，由（2.31）式，可得

1(aibixciz)2 (1-43)

1Nj(ajbjxcjz)2Ni可将（2.43）式看作是从（x,z）到（Ni，Nj）的坐标变换。由（2.41）和（2.42） 式可知，这种变换将（x,z） 平面上的i(xi,zi)、j(xj,zj)、m(xm,zm) 点分别变换为（Ni，Nj） 平面上的i'(1,0)、j'(0,1)、m'(0,0) 点（见图2.4），其面积单元之比等于雅可比(Jacobian)式

NidNidNjxNjdxdzxNibiz2Njbjz2ci21(bcbc)

ijjicj422考虑到（2.28）式

bicjbjci2

故有

dxdz2dNidNj （2.44）

由此可计算

(ijm)e[Ni(x,z)]2dxdz22NidNidNj ijm'e'''1Ni2NidNjdNi

00122Ni2(1Ni)dNi01 （2.45） 6同样

(ijm)eNi(x,z)Nj(x,z)dxdz2NNdNdNijiijm'e'''j

1Ni1N2NiNjdNjdNi2i(1Ni)2dNi 020012=(Ni2Ni2Ni3)dNi01 （2.46） 12可以将以上两式类推到（2.38）式中的其余积分，并写成如下一般形式

Np(x,z)Nq(x,z)dxdze(1pq) （2.47） 12式中

pq10(pq)(pq) p、q分别表示i、j 或m。

将（2.47）式代入（2.38）式得

2vdse22(viv2jvmvivjvjvmvmvi) （2.48） 6现在，只剩下(2.33)式中最后一项（线）积分了，而此项积分可仿照第二项（面）积分的方法求得。应该指出，此（线）积分是由求解区边界2上的边值条件引入的，故只需对属于2 的单元边界计算该积分。设单元边界jm 在求解区边界2 上，可以近似地将jm 上的系数视为常数（取为j或m点的值j 或

m）；同时，电导率e亦为常数，皆可提到积分号外。于是，问题归结为计算积分

2vdljmjmN(x,z)viiNj(x,z)vjNm(x,z)vmdl

2根据前述Ni、Nj、Nm的几何意义，可以写出它们在jm 边上的表示式（参见图2.5）

Ni(x,z)pjm/0 Nj(x,z)pmi/1l/ljm Nm(x,z)pij/l/ljm

故有

图2.5 含边界的三角单元上线积分计算用图 ljm2vdl jm[(10ljmlljm)vjlljmvm]2dl

ljm2ljm2lllll2222)dlvm2dl2vjvm(2)dl =vj(12llljmljmjmjm00ljm0=

ljm32(v2 jvmvjvm) jm2 （2.49）

如果ij或 mi 位于边界2上，则可类推得

2vdlijlijv32iv2jvivj ij2 （2.50）

或

mi2vdvlmi2(vmvi2vmvi) mi2 （2.51） 3将（2.35）、(2.37)、 (2.48)、（2.49）、（2.50） 和（2.51）代入（2.33）式，得单元e上的泛函

Jee(bv4iibjvjbmvm)2(civicjvjcmvm)2

e222(viv2 +jvmvivjvjvmvmvi) 6IiviIjvjImvm -()

ninjnmeijlij(vi2v2jvivj)/32ejmljm(v2jvmvivj)/322l(vvvmvi)/3emimimi(ij2)(jm2) （2.52） (mi2)这样，便将单元泛函表示成该单元各节点函数值vi、vj和vm的多元函数，将所有单元的Je(v)相加，便得到整个求解区的J(v)

J(v)Je(v) （2.53）

e于是，前述变分问题（泛函的极值问题）就变为多元函数的极值问题。我们知道，多元函数J(v)取得极值的必要条件是其偏导数为零。即J(v)对每一变量v1、

v2、…、vM的偏导数为零

J(v)vrvrJee(v)eJe(v)0 （2.54） vr(r1,2,,M)

我们仍就从一个单元分析入手，由（2.52）式可计算包含i、j、m三个节点之单元e的泛函Je 对各节点的函数值的一阶偏导数：

Je0 v1…

2ilij/321Jee(bibicici)0vi32l/32imi(ij2)(jm2)vi

(mi2)jlij/3210+e(bibjcicj)6020210+e(bibmcicm)6l/32mimieeeKiiviKijvjKimvmIi/ni

(ij2)(jm2)vj

(mi2)(ij2)I(jm2)vmi

ni(mi2)…

JeeeKevKvKjiijjjjmvmIj/nj vj…

JeeeeKmiviKmjvjKmmvmIm/nm vm…

Je0 vM可以将以上各式写成矩阵形式：

i列 j列 m 列

Jev1JeviJevjJevmJevMeKiiKejiKemieKijKejjeKmjeKimKejmeKmmv1viIi/nii行vjIj/njj行 vmIm/nmm行vM 矩阵[Ke]及列矢量{Ie}中的虚点均为零元素或记为

{Je}[Ke]{V}{Ie} (2.55)

式中：M维列矢量{Je}的元素为

Je （r=1,2,…，M） （2.56） vrJreM维列矢量{V}的元素为各节点的函数值vr（r=1,2,…，M） M维列矢量{I}的元素为

Ir/nrI0ere

(ri,j,m)(ri,j,m) (2.57)

MM阶单元矩阵[Ke]的非零元素为

Krse1G2Gslsr'e(brbscrcs) (2.58)

632s,r,r'i,j,m;sr';2G1(rs)(rs)

eKrs表示式中最后一项仅当单元e的边sr'在求解区边界2上才出现；其中的s是

按（2.20b）式对节点s算出的系数值。

3）总体合成

(2.55)式是整个求解区中某一个单元e的泛函Je对各节点的函数值的一阶偏导数所形成的线性方程组的矩阵形式。将所有单元的Je相加，便得到整个求解区的泛函的变分

JJee=0

（2.59） 写成M维列矢量

{J}{Je}0 （2.60）

e({J}的元素为JrJ，r1,2,,M) vr将（2.55）代入（2.60）式得

[K]VI0

eeee或写成

KVI （2.61）

式中，刚度矩阵K[Ke]，其元素为所有单元矩阵[Ke]的相应元素之和

eKrsKrs （2.62）

ee(r,s=1,2，…，M)

场源列矢量IIe，由（2.55）可知，其元素IrIre为节点r上的供电

ee电流强度（r=1,2,…，M），如果供电电极不在节点r上，则Ir0。

（2.61）式便是由变分问题（2.22）式离散化后导出的线性方程组。解此方程组便可对给定的波数确定各节点电位的傅氏变换电位v，并可通过反傅氏变换计算各节点的电位值u。然后，根据公式（2.23）即可计算出视电阻率值。

2.2 二维地电构造中点电流源场的正演有限元算法计算的相关技术

2.2.1方程组的简化

采用图2.3所示的三角形网格能对求解区作比较精细的剖分，因而模拟地形或电性异常体及单元内作线性插值的近似程度都较好。但其缺点是节点数目较多，所形成的线性方程组阶数较高，因而计算量较大。下面介绍一种算法，把图2.3中所示网格中位于各矩形中心之节点函数值从方程（2.61）中消除，以降低待求解的线性方程组的阶数，从而大大减少计算量。

我们考察图2.3所示网格中的一个小矩

图2.6 网格中的一个小矩形单元 Figure2.6 Little rectangle unit in mesh 形单元E（图2.6），其四个角的节点编号设为i、j、m、k ，显然有

kMi (2.63)

mMj式中,M为矩形网格的纵向节点数（参见图2.3），该矩形单元的中心节点设为a，其与四个角节点相连，将此矩形单元分为四个三角形单元，依次表示为e1,e2,e3和e4，它们的面积都相同，且为

1hxhz (2.64) 4其中，hx和hz分别为矩形单元沿x和z方向的步长。

设各三角形单元内电导率是均匀的，但各三角形单元之间，电导率可以不相同，分别以1、2、3、4表示e1、e2、e3、e4内介质的电导率。按(2.58)式可写出各三角形单元(e1、e2、e3、e4)的矩阵[Ke1],、[Ke2]、[Ke3]和[Ke4]，

它们之和便给出矩形单元E的矩阵KE

[KE][Ke1][Ke2][Ke3][Ke4]

（2.65）

在计算这些矩阵的元素时应注意，在当前情况下，br和cr(r=i,j,m,k或a)都可表示为十分简单的形式。比如，在e1中，按（2.26）式有：

bizjzahz/2， cixaxihx/2 bjzazihz/2，cjxiaxahx/2 bazizjhx ，caxjxi0

于是，按（2.58）式可算出[Ke1]的各非零元素为

12212 K(bici)23e1iihxhz12hz2hx2()= 1434hxhz421hx12=()hxhz 2hxhz121hzKe1jj12(bc)2j2j123

2hz2hxhxhz12() = 14434hxhz21hx12 =()hxhz

2hxhz121hzKe1aa12(bc)2a2a123

hxhz12(h0) = 134hxhz212zhz12 =21.hxhz

hx12Ke1ijKe1ji12(bibjcicj)126

hz2hx2hxhz12() = 1464hxhz421hx12 =()hxhz

2hxhz241hzKe1iaKe1ai12(bibacica)126

hz2hxhz12(0) = 1264hxhz21hz12 =1hxhz

hx24Ke1jaKe1aj12(bjbacjca)126

hxhz12hz2(0) = 164hxhz221hz12 =1hxhz

hx24同样可计算出[Ke2]的各非零元素

Ke2jjhx22()hxhz 2hxhz122hzKe2mmhx22()hxhz 2hxhz122hzKe2aahz2222hxhz

hx12Ke2jmKe2mjhzhx22()hxhz 2hxhz242Ke2jaKe2ajhz222hxhz

hx24hx222hxhz

hz24Ke2maKe2am[Ke3]的各非零元素

Ke3mmhx32()hxhz 2hxhz123hzKe3kkhx32()hxhz 2hxhz123hzKe3aahz3223hxhz

hx12Ke3kmKe3mkhx32()hxhz 2hxhz243hzKe3maKe3amhz323hxhz

hx24hz223hxhz

hx24Ke3kaKake3[Ke4]的各非零元素

Ke4kkhx42()hxhz 2hxhz124hz4hzKe4iihx42()hxhz 2hxhz12hx4224hxhz

hz12Ke4ikKe4aaKe4kihzhx42()hxhz 2hxhz244Ke4kaKe4akhx424hxhz

hz24hx424hxhz

hz24Ke4iaKe4ai将以上各式代入（2.65）式，便可算出矩形单元Ｅ的矩阵[KE]之各非零元素，比如

同样可得 =KEKe1e2e3e4iiiiKiiKiiKii

=ke1ii00ke4ii

1hzh2x14hzhx2hh)h()2=(4xhzhxhz xz122hxhz12=

14(hzhx22h6hxhz) (2.66a) xhzKEKeKe12hzhxjjjj1jj22hh2(hxhz) (2.66b) xz6KEehzh2mmKmm2Ke323mm2(hxhxhz) (2.66c) xhz6KEe4kkKekk3K(hzhkk432hx2hxhz) (2.66d) xhz6KEe1e2e3e4aaKaaKaaKaaKaa

2(12)hzh2(h234224)x1hxhz (2.66e) xhz12EKEKe1hzh2Kxijjiij12(h12hxhz) (2.66f)

xhzKE2jmKEmjKe2hzh2jm2(xh12hxhz) (2.66g) xhzKEmkKEkmKe33(hz2hhxh2mkhxhz) (2.66h) xz12

KKKEkiEike4kihzhx2(hxhz) (2.66i) 2hxhz124KKKEiaEaie4iahxhz(14)2K41hxhz (2.66j)

hzhx24e1iae2jiaKEjaKKKEaje1jahx(12)2hz12hxhz (2.66k)

hxhz24hxhz(23)223hxhz (2.66l)

hzhx24KEmaKEamKe2maKe3maKEkaKEakKKe3kae4kahx(34)2hz34hxhz (2.66m)

hxhz24除上列各元素之外，矩阵[KE]的其余元素皆为零。

由于节点a位于矩阵单元E内部，其相邻的节点只有E的四个顶角（i、j、m、k）,而与其它节点无直接联系，故除所论单元矩阵[KE ]外，其它单元矩阵内部都不包含与a有关的元素，即刚度矩阵[K]= KE 中第a行和第a列上

E的元素，都与单元矩阵[KE ]中之相应元素相同。其中的非零元素只有

EKaaKaa

EKaiKiaKia

KajKjaKEja

EKamKmaKma EKakKkaKka

由此，可将线性方程组（2.61 ）中第a个方程及其它（共四个）包含有节点函数值va的方程写出如下

EEEEEKaavaKaiviKajvjKamvmKakvk0EEEEEKiavaKiiviKijvjKimvmKikvkIiEEEE KEvKvKvKvKvI （2.67）jaajiijjjjmmjkkjEEEEEKmavaKmiviKmjvjKmmvmKmkvkImEEEEEKkavaKkiviKkjvjKkmvmKkkvkIk应该指出，以上第一个方程的右端项Ia0，这是因为矩形单元中心节点a是一个虚设的节点，供电极不会置于该处，故该点的供电电流强度Ia为零；其次，除第一个方程之外，其余四个方程都包括有别的项，这里没有全写出来；此外，

EEEKsa除与a有关的系数（Kas ,s=i、j、m、k）外，其余系数一般KrsKrs(r,s=i、

j、m、k),即合成矩阵[K]与单元矩阵[Ke ]的对应元素一般并不相同，而上述方程中的系数乃是合成矩阵[K]的元素。

通过消元一次，将va从（2.68 ）式的后四个方程中消去，于是得到

EEEEEKaavaKaiviKajvjKamvmKakvk0''''KiiviKijvjKimvmKikvkIi K'jiviK'jjvjK'jmvmK'jkvkIj （2.68）

''''KmiviKmjvjKmmvmKmkvkIm''''KkiviKkjvjKkmvmKkkvkIk式中

'EEEKiiKiiKiaKai/Kaa'EEEKmmKmmKmaKam/Kaa'EEEKijK'jiKijKiaKaj/Kaa''EEEKmkKkmKmkKmaKam/Kaa''EEEKimKmiKimKiaKam/KaaEEK'jjKjjKEK/Kjaajaa'EEEKkkKkkkkaKak/kaa'EEK'jmKmjKjmKEjakam/Kaa ''EEEKkiKikKkiKkaKai/Kaa'EEK'jkKkjKjkKEK/Kjaakaa从（2.68 ）式可看出，经过一次消元处理后，便只有第一个方程中包含有节点函数值va 。在作数值模拟时，通常并不需要确定矩形中心节点的函数值va。这样，就可去掉该节点函数值和包含它的方程。对每一个矩形单元都照此处理，便可去掉大量（将近一半）待定节点函数值和方程数，从而使线性方程组的阶数大为降低，使计算量成倍地减少。

从（2.68）式还可以看出，消掉节点函数值va将使刚度矩阵[K]的有关元素

发生变化，改变量只与单元矩阵[Ke]中a行或a列上的元素有关。为计算方便，只需把这些改变量引入单元矩阵[Ke]中的相应元素即可。这样，对各个矩形单元便组成新的单元矩阵[Ke]，其非零元素为

(E)E2EE2EKiiKiiE(Kia)/Kaa K(jjE)KEjj(Kja)/Kaa

(E)EE2E(E)EE2EKmmKmm(Kma)/KaaKkk(Kka)/Kaa Kkk

(E)EE(E)(E)EEEEKijK(jiE)KijEKiaKEja/Kaa KjmKmjKjmKjaKma/Kaa (E)(E)EEEE(E)(E)EEEEKmkKkmKmkKmaKka/KaaKikKkiKkaKia/Kaa Kki (E)(E)EEEE)(E)EEKimKmiKiaKma/Kaa K(jkKkjKkaKE/Kjaaa

如果矩形单元E紧靠着求解区边界，则上式中的某些元素还应加上与混合边值条件有关的项。由（2.58）式可知，若E紧靠左边界（节点i和j在边界2上），则下列元素改为

2(E)E2EKiiKiiE(Kia)/Kaaihz1 (2.69a)

32E2EK(jjE)KE(K)/Kjhz1 (2.69b) jjjaaa31(E)EEKijK(jiE)KijEKiaKE/Kjhz1 (2.69c) jaaa3若E紧靠右边界（节点m和k在边界2上），则需改变的元素为

2(E)EE2EKmmKmm(Kma)/Kaamhz3 (2.70a)

32(E)EE2EKkkKkk(Kka)/Kaakhz3 (2.70b)

31(E)(E)EEEEKmkKkmKmkKmaKka/Kaakhz3 (2.70c)

3若E紧靠底边界（节点j 和m在边界2上），则需改变下列元素

2E2EK(jjE)KE(K)/Kjhx2 (2.71a) jjjaaa32(E)EE2EKmmKmm(Kma)/Kaamhx2 (2.71b)

31E)(E)EEEK(jmKmjKEKK/Kmhx2 (2.71c) jmjamaaa3将(2.66)式代入(2.69)、(2.70)或(2.71)式，便可具体计算出新的(或最终

的)单元矩阵[K(E)]的各非零元素。而各矩形单元矩阵[K(E)] 之和便构成新的(或最终的)刚度矩阵[K]

KK(E) （2.72）

E此新的刚度矩阵[K]不再包括代表各矩形中心节点的行和列,它是一个阶数等于求解区内矩形网格节点数的方阵，其规模比先前包括矩形中心节点时大为减小。同时,按(2.61)式构成的线性方程组的节点函数值列矢量{V}和场源列矢量{I},也都不包括各矩形中心节点的对应量，其维数也等于求解区内矩形网格上的节点数。自然，解此线性方程组最终算出的结果,只包括这些矩形网格节点的函数值。

应该指出,上述消掉矩形中心节点函数值，使线性方程组降阶的方法，不同于直接对求解区进行矩形网格剖分的作法。它本质上仍然是作较细的三角网格剖分，因而对起伏地形及电性异常体的模拟和线性插值处理的近似程度都比直接作矩形网格剖分好。不过，它形成刚度矩阵[K]的计算量较后者大,并且存放各三角形单元电导率所占用的存贮量也较多。为了克服此缺点，可对内部电性均匀(1234)和不均匀（1、2、3和4不完全相同）的矩形单元分别进行处理，对导电性均匀的矩形单元（在整个求解区内占大多数），只存贮一个电导率值，并且可用简化的公式计算[K(E)]的元素。

2.2.2 线性方程组的解法

有限元法导出的线性方程组（2.61）式系数矩阵（刚度矩阵[K]）具有以下特点：

1）由于Krs=Ksr，即单元矩阵[K]是对称的，而总体的刚度矩阵[K]乃是单元矩阵[Ke]之和，所以刚度矩阵[K]也是对称的。此外，数学上还可以证明[K]是正定的； 2）刚度矩阵的阶数等于求解区的总节点数。通常，在求解高密度电阻率法二维正演问题时，节点数可达几千甚至更大，所以[K]又是大型的。不过，在刚度矩阵[K]中，只有对角元素及其附近的元素不为零，其余皆为零元素，故一般地说，[K]是稀疏矩阵。并且，如果对节点作恰当的编号，使中心节点与其相邻节点之

e

e

e

编号相差尽量小，则刚度矩阵[K]的非零元素将比较集中地分布在对角线附近的带状区域内，因此，[K]又是带状的。

总之，有限单元法在求解电法勘探正演问题时，所得到线性方程组的系数矩阵[K]是一个大型、稀疏和带状的对称、正定矩阵。有限单元法大部分时间要花在解线性方程组上，特别是在进行反演时，解方程组的耗时所占比例更大。因此，选用合适的方程组的解法，也是十分重要的。在有限单元法中，求解线性方程组最为常用的方法有迭代法（如SOR法和CG法）和直接法。实践表明，对这种大型、稀疏和带状的对称、正定系数矩阵，最好采用乔勒斯基（Cholesky）分解法求解。

乔勒斯基分解是将对称、正定的系数矩阵[K]分解为一个同阶的下三角矩阵[L]与其转置矩阵[LT]之积（故也称LLT分解），即

[K][L][L]T （2.73） 于是方程（2.61）变成

[L][L]T{U}{I} （2.74） 它等价于以下两个方程

[L]T{U}{W} （2.75）

[L]{W}{I} （2.76）

先解第二个方程组，按如下“回代”公式计算中间变量{W}的各元素

Wi(IiLipWp)/Liip1i1(i1,2,,N) （2.77）

然后，以算出的中间变量{W}作为已知右端项，求解第一个方程组，按如下“回代”公式计算{U}的各元素

U(WLU)/Liipipiipi1NN（2.78） (iN,N1,,1)

以上两式中约定(....)0和(....)0

1N10下三角矩阵[L]的非零元素Lij按下列递推公式算出

j1(KijLipLjp)p11i12Lij(KiiLip)2p10jiji （2.79）

ji(i1,2,,N)在按（2.79）式作乔勒斯基分解，计算下三角矩阵[L]的对角元素时，包含几次开方运算，它的计算量较大，计算精度稍差。为避开开方运算，可采用乔勒斯基分解的变形━LDLT分解，即将系数矩阵[K]分解为三个同阶矩阵之积 [K][L][D][L]T （2.80）

其中，[L]为单位下三角矩阵，[D]为一对角矩阵，它们的元素可按下列递推公式算出

j1~~LaijKijaipjpp1 ~/dLijaiji~LdiKiiaipipp1i1ji （2.81） ji(i1,2,,N)在对[K]作了LDLT分解之后，方程（2.61）变成

[L][D][L]T{U}{I} （2.82） 它等价于以下两个方程

[L]{W}{I} （2.83） [L]T{U}[D]1{W} （2.84） 它的解可按下列“回代”公式算出

WiIiLipWpp1i1(i1,2,,N) （2.85）

和

UiWi/diLpiUppi1N(iN,N1,,1) （2.86）

上述解线性方程组的直接法，相对于迭代法的缺点是要占用较多计算机储存单元。所以，在电子计算机上实现上述计算时，还应注意以下技巧，以尽量减少计算量和储存量：

（1）因为刚度矩阵[K]是对称的，所以只需计算和储存其一半元素，通常只需形成和储存下半个矩阵；

（2）尽量重复利用刚度矩阵[K]、场源列矢量{L}以及其中间矩阵来储存新生成的矩阵（如[L]等），而不是再开辟新矩阵；

（3）因[K]是稀疏的，因此，甚至并不需要储存其下半个矩阵的全部元素，而只储存其中的非零元素，因而可进一步大量节省计算机的储存量，这种储存方法称为“紧缩存贮方法”。数学上已证明，带状的对称、正定矩阵[K]作LLT或LDLT分解后，下三角矩阵[L]也是带状的，且其半带宽S不超过[K]的半带宽，故在采用紧缩存贮的条件下，[K]和[L]仍可先后占用相同的存贮单元。

如前所述，在采用图2.3所示矩形求解区和网格剖分时，实用的只是矩形网格上的节点，由于求解区通常沿水平（X）方向较长，节点数（N）较多，而沿垂直方向（Z）较短，节点数（M）较少，故将节点从左边开始，逐列由上至下进行编号。这样，刚度矩阵[K]的半带宽S最小（S=M+2）。

由于对称性，实际上只需要存贮[K]或[L]下三角中半带宽内的元素，故采用二维数组的形式存放这些元素。二维数组的列数等于半带宽S，而其行数与总节点数（N×M）相同。采用这种紧缩存贮方式，对用LLT或LDLT分解法解有限单元法的线性方程组是经济和方便的。并且，由于采用上述紧缩存贮，形成刚度矩阵[K]和分解、计算下三角矩阵[L]所涉及的元素数目都大为减少，所以不仅节省了存贮单元，而且也节省了计算量。

在计算高密度电阻率法正演问题时，需要计算不同供电电极的位置的电场分布。改变供电电极的位置不仅影响线性方程组的右端(场源)项{I}，而且使刚度矩阵[K]中与2边界节点有关的元素发生变化，这是因为（2.69a2.71c）各式中的[IKK1(rK)cosk]k2n[Ik1nKK0(rk)] 与供电点A位置有关。为了减

少计算量，在电极距不太大或求解区范围足够大时，可选用装置（主要是供电极）处于某一适中位置时的η值，固定不变。当装置移动或变化，供电电极位置改变时，只改变线性方程组的（2.61）的右端项{I}，而其系数（刚度）矩阵[K]保持不变。这样，在计算不同供电位置的节点函数值时，只需形成和分解一次刚度矩

阵[K]；而对每一供电位置，改变其右端项后由已分解好的[L]和L矩阵作回代，

T便可算出相应的节点函数值。因为与形成和分解矩阵[K]相比，回代的计算量甚小，所以采用上述处理方法计算多个供电位置的电场，增加的计算量并不大。这是提高电剖面法计算速度的一条重要措施。

2.2.3 视电阻率的计算

对于给定的供电电极位置经反傅氏变换（见附录B）计算出给定的各地面节点之电位后，便可进而计算给定装置给定位置上的视电阻率s。 比如，对于温纳装置（在高密度电阻率法中有人称之为

装置，

AM=MN=NB=na）。若给定的供电极A（I）和B（-I）之节点编号为nA和nB，则方程组（2.61）的右端列矢量{I}，除第nA和nB个元素分别为InAI和InBI之外，其它元素皆为零。由此算出A、B间各节点的电位值后，便可按下式计算测线上不同位置的视电阻率 sUnMUnNKI （2.87)

式中，UnM和UnN分别是测量电极M和N所在节点（编号为nM和nN）的电位值；装置系数K按下式计算

K(xAxM)(xAxN)2na (2.88)

xMxN式中的xA,xB,xM,xN分别是节点nA,nB,nM,nN的横坐标值，a为最小电极距，n为电极隔离系数。

若使用偶极装置（在高密度电阻率法中有人称之为装置，AB=MN=a，AM=BN=（n+1）•a）则

Kn(n1)(n2)a (2.89) 若采用三极装置（AM=na ，B->∞，MN=a），则

K4na (2.90)

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z2221 (4.6)

2Is2[(xa)yz]22212 (4.7)

这里，测量电极位于点（a，0，0）。在计算出、'的散度后，式（4.5）可写为

Isx(xa)y2x22dxdydz (4.8) 33222222v4xyz2[(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2E(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2E(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 (4.18) Fy3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2a16a

4.4.2 计算偏导数矩阵

为了去获得有限范围内矩形元的偏导数值，需要计算一个关于Fy函数的双重积分（见公式4.11）。很明显，公式（4.19）、（4.20）没有简单的解析解，这样，我们有必要利用高斯求积方法去计算其数值解。通常情况下，高斯求积方法比其他方法（梯形求积方法、龙贝格求积方法）更能得到一个比较精确的解。其近似公式为

Is42z1z2x2x1Fyx,zdxdz

AIs11Fu,vdudv (4.22) 211y4AIs42WL1K1lnnkKWEFyu,vdudv

式中

u=(2x-x1-x2)/(x2-x1), v=(2z-z1-z2)/(z2-z1), A=0.25(x2-x1)(z2-z1)

注意转换后横坐标u、v的积分范围是[-1，1]。参数nx、nz是x和z方向的求值数。wk和we是权系数，以上各参数的值可参见Churchhouse（1981）。

计算数的大小应依据所计算的矩形元与电极间的距离而定。当矩形元靠近电极时，Fy随x、z的变化很快。当矩形元最接近某一个电极时，x、y方向的计算数分别是10和8。当矩形元跨越两个电极间距以上时，x、y方向的计算数分别取4和3。随着矩形元与电极之间的距离的增加，x、y方向的计算数逐渐减少。

为了计算雅可比矩阵中的每一个元素，需要计算二维模型中所有可能的两根电极（电流电极和电位电极）的组合条件下所有矩形元的偏导数值（见图4.2）。图4.2中有21根电极，可划分为63个矩形元，总的组合数是26460（21×20×63）。实际上，由于所讨论的问题具有对称性，计算量可以减半。此外，对于均

匀介质模型，偏导数矩阵中的很多元素是一样的。例如，对于矩形元2对电极2、3的偏导数与矩形元3对电极3、4的偏导数是一样的；并且与矩形元4对电极4、5的偏导数是一样的等等。这样，偏导数矩阵的值可减少到7200个。

一个矩形元的偏导数值仅仅依赖于边角处的矩形元对于电极跨度x与z的比。假如电极之间的间距衡为常数，反演所用的模型中矩形元的排列也保持同样的水平间距，那么这个比值将保持不变，并与电极间距无关。这样，偏导数矩阵的值仅仅只需计算一次，这样就大大地节约了计算偏导数矩阵的时间。需要说明的是，仅仅需要计算单极-单极排列的偏导数矩阵，因为高密度电阻率法中任意一个电极序列都可以从单极-单极排列中引出。

当然，偏导数矩阵即可由有限元或有限差分方法得到（Sasaki，1992）。Loker和Barker（1994）用有限差分方法与上述解析方法算出的偏导数矩阵进行了对比，其均方误差小于5%。

4.4.3 方法步骤

（1）均匀介质初始模型的电阻率q0可设置为测量的视电阻率值的平均值，即

1mq0fi (4.23)

mi1（2）雅可比矩阵计算出来后，再选择一个合适的阻尼因子（通常为0.05），光滑约束最小二乘公式（4.1）就建立起来了。阻尼因子的值依赖于数据中的随机噪声水平（Sasaki，1992）。噪声水平高，阻尼因子的值要大一些；噪声水平低，阻尼因子的值可适当小一些。对同样大小的矩形元，电极排列的响应随着深度的增加而降低。在较深的位置（隔离系数较大）上，为了保持反演过程的稳定，平滑矩阵中的元素的值应增大（Sasaki，1989）。每增加一个隔离系数，其值大约提高10%。

（3）解最小二乘公式（4.3）就可以得到模型参数的改正量。其迭代公式是

q1q0p (4.24)

因为受所选的阻尼因子的影响，因此，最好是用不同的阻尼因子重复计算，这样将获得较好的结果。

数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数数|数数数数

加一个：视参数等值线断面色谱图图

（三）高密度电阻率法的应用 1.野外工作技术 （1）测网布置

地球物理工作的测区一般是由地质任务确定的。对主要应用于工程及环境地质调查中的高密度电法而言, 按工程地质任务所给出的测区往往是非常有限的，我们只能在需要解决工程问题的有限范围内布设测线、测网，可供选择的余地往往很少，这是一般工程物探经常遇到的情况。测网布设除了建立测区的坐标系统外，还包含了技术人员试图以多大的网度和怎样的工作模式去解决所给出的工程地质问题，在这里，经验和技巧非常重要。特殊情况下，高密度电阻率法可布设不规则的测线和测网，尽可能在有限的测区内获得更多的测量数据。 （2）装置选择

如Wenner，Dipole-Dipole、Pole-Dipoler和Wenner-Schlumberger装置等。通常使用的装置还如上述四种类型。不同的测量系统基本上以这几种装置为主, 但也各有特点, 有的高密度电阻率仪提供了十多种装置以供选择。不同装置可联合使用, 也可根据需要单独使用。选择一个合适的工作装置应考虑：（a）探测目标的特性、(b)仪器灵敏度以及(c)场地噪声本底水平，更要考虑：d）装置对地下电阻率水平或垂向变化分辨能力、(e)探测深度、（f）有效探测范围以及（g）信号强度。

（3）最小电极距和排列长度的选择

最小电极距和排列长度的选择取决于地质对象的大小和埋藏深度。要保证有足够的横向分辨率，探测目标体横向上至少要有2～3根电极通过。同时，由于高密度电阻率法实际上是一种二维探测方法,所以在保证最大极距能够探测到主要地质对象的前提下，还要考虑围岩背景也能在二维断面图中得到充分的反映。如对小而深的探测目标体，要求较小的电极间距和较多的电极数。

对于长剖面，可以通过电极的移动来获得连续的断面数据。图？是温纳-斯

伦贝谢装置通过两次移动来获得18x剖面长度的例子。一般地，在剖面对接时要重叠三个点，重叠点的数据取两次测量的平均值。

图 温纳-斯伦贝谢装置移动测量示意图

2. 高密度电阻率法在工程与环境地质中的应用

近年来, 高密度电阻率法在场地勘察、公路及铁路隧道选线、坝基及桥墩选址、采空区及地裂缝调查以及水库渗漏研究等领域得到广泛应用, 取得了明显的地质效果和显著的经济效益。下面用几个高密度电阻率法的实际例子来说明该方法的应用。

（1）在煤气管道探测中的应用

场地地形有微小的起伏，测线左侧是水田，右侧是沙石公路。水田一侧地势较低，公路一侧地势较高。实测时，最小电极距为0.3m，电极数为30，N=9。

图7a 是温纳装置视电阻率等值线断面图。从图中可以看出，地下介质视电阻率在30～80m范围内，反映了地下耕织土及其下部亚粘土的电阻率；在16～24号电极之间，浅部出现100～120m相对高阻电阻率值，是由于沙石公路引起的。大约在16～17号电极之间，见明显的视电阻率等值线封闭圈，指示了煤气管道的位置。其上部视电阻率等值线局部密集，这是由于采用明挖埋设煤气管后填埋性质不同的渣料所致。

图b，是Ts视参数等值线断面图。在煤气管道位置，出现了Ts高值等值线圈，其异常较s异常明显。

煤气管道中心位置埋深对应N=4，此时装置极距为1.8m，估算埋深大约为0.6m左右，这一结果也为雷达探测和实地调查结果所证实。

1x23456789101112131415161718192021222324252627282930x0.0n=11.02.03.04.05.06.07.08.0(m)n=2n=3n=4n=5n=6n=7n=8n=9

(a)

1x234567891011121314151617181920212223242526272829x0.0n=11.02.03.04.05.06.07.08.0(m)n=2n=3n=4n=5n=6n=7n=8n=9 (b)

1x23456789101112131415161718192021222324252627282930x0.0n=11.02.03.04.05.06.07.08.0(m)n=2n=3n=4n=5n=6n=7n=8n=9

(c)

图煤气管道上温纳装置视参数等值线断面图

(a)温纳装置s等值线断面图；(b)Ts等值线断面图；(c)Ts等值线断面色谱图 （2）在城市自来水管探测中的应用

测量场地为城市水泥街道，宽18M。探测目标为路面下顺着街道延伸的城市污水排水管道，目的是确定管道的位置。测量场地水泥路面表层混凝土厚20cm，其下为水稳层，厚约25～40cm。再往下为压实层，厚约1.0m，主要有砂、石和

土组成，电阻率较高，约为300~400m。压实层以下过度到自然状态粘土层，厚6.0 m左右，电阻率约为40～60m。基底是强风化泥质粉沙岩，电阻率约为15～25m。

测试路段水泥路面混凝土配合比为四车（建筑工地常用的手推翻斗车）石、三车沙、两包水泥。按此比例，并用同一类型沙、石料和同一型号水泥制作了多块混凝土标本。28天后，标本测量结果为49m，混凝土路面实际上是导电的。水泥主要由硅酸盐组成，其中夹有少量的煤灰和粘土。分析混凝土导电的原因，含有少量的煤灰是一个因素，其次可能是在混凝土凝固后，夹在其中的水份象小气泡一样成蜂窝状永久地保留在里面，从而使得混凝土导电。

我们试验了多种电极，其中以由半截酸奶瓶、铜钉和掺有食盐的粘土所组成的接地电极效果最好，且制作方便，可重复使用，如图1所示。

铜钉粘土酸奶瓶

实际工作中，为保证电极与水泥地面接触良好，必须保持瓶中粘土有一定的湿度。一般在跑极时，电极应放在盛有泥浆的塑料盆里，当路面上有较厚的灰尘时，应先用湿抹布把灰尘擦净。

该场地水泥路面下埋设有横穿街道的自来水水管，在电性也表现为低阻。高密度电阻率法测量时使用了30个自制电极，垂直于水管对称布置。最小电极距

a=0.3m，剖面长度=8.7m。实测时，获得了s、s、s三种参数，室内进行了比值处理，计算了Ts参数。s、s、s等值线断面图均很好地反映了水管位

置，特别是Ts等值线断面图对水泥路面及其以下各介质分布特征反映相当清晰，见图6。

1x23456789101112131415161718192021222324252627282930x1.02.03.04.05.06.07.08.0(m)n=1n=2n=3n=4n=5n=6n=7n=8n=9

自来水管上温纳装置视电阻率等值线断面图

从视电阻率等值线断面图可以看出，异常的分布特征与地下介质的电性特征

是一致的。

此外，在水泥路面上，我们还进行了二极、三极剖面法测量，均获得了良好的效果。

m360340320300sAsB2802602402202001801601401.02.03.04.05.06.07.08.0(m)

自来水管上联合剖面装置视电阻率等值线断面图

（3）在溶洞探测中的应用

沪-昆高速公路通过浙江省常山县某灰岩区，该区地下溶洞发育。测量时采用三极滚动装置，布设总电极60根，分为两串，每串30根。最小相邻电极距

x=2m，N=1～30，最大电极距AO=59m。扫描N=1～30正好一串，这样就可以将前一串电极移至第二串电极之后继续进行测量，直至整个测线测完为止。图5是其中某路段高密度电阻率法探测图像，正演采用有限元法，反演采用光滑约束最小二乘法反演法，在电阻率图像上，溶洞表现为低阻，围岩相对为高阻。本结果已为续后的钻探资料证实。

水平距离(m)00.45.19.616.620.040.060.080.0100.0120.0140.0160.0180.0200.0深度(m)11.223.549.4103溶洞217455953溶洞1998

图5 沪-昆高速公路浙江常山境内某灰岩段电阻率二维反演解释成果图 （4）在隧道工程地质勘查中的应用

浦城－南平高速公路从桩号K80+300.000～K164+260.000路段地貌及地质条件较差，路线穿过最高山海拔1091米，跨越崇阳溪、麻阳溪两条较大河流，在84公里范围内隧道总长约10公里。因此隧道位工程地质条件的探测是路线勘查及设计最重要的内容之一。其中祝源隧道最长，左线4270m（起点桩号：Zk80+105.000，终点桩号：Zk84+375.000），右线4340m（起点桩号：Yk80+090+000，终点桩号：Yk84+430.000）。

共布设电极120根，最小相邻电极距x=10m，测量时采用三极滚动装置，N=1～90。为了增加探测深度，除了加大供电电流外，还采用了一些特殊的技术措施。

祝源隧道电阻率二维反演解释断面色谱图

祝源隧道电阻率二维反演解释断面彩色分级断面图

结合遥感、地调与钻孔资料，推断了F1～F10共10条构造，见表。

隧道 构造编号 名称 F1 F2 F3 F4 祝源 隧道 F5 F6 F7 F8 F9 F10 位置 YK80+340 YK80+850 YK81+468 YK82+300 YK82+800 YK83+100 YK83+500 YK83+770 YK83+900 YK84+460 倾角 ∠85° ∠80° ∠70° ∠80° ∠85° ∠65° ∠70° ∠70° ∠65° ∠85° 规模较大，向下延伸至隧道洞身 规模较大，向下延伸至隧道洞身 规 模

图？？是六（安）—黄（尾）高速yk48+950.000～yk49+250.000段转步园隧道电阻率二维反演解释断面色谱图。采用温纳装置，总电极数60，最小相邻电极距x=10m。

图？？是根据该隧道地震资料绘制的转步园隧道地质解释剖面图，与地调与钻孔资料吻合。

转步园隧道地质解释剖面图

从以上两个实例可以看到，一般情况下，高密度电法解释必须结合其它物探资料以及钻探、地质、水文和物性资料等才能得到合理的结论。 （5）在防空洞探测中的应用

地下人防工事探测是城市物探中经常遇到的难题之一。防空洞虽然埋深不大，但空间位置复杂，填充物不明，且地表大多是杂填土，并堆积了厚厚的建筑垃圾，受场地范围限制，探测条件很不理想。图6是经过高密度电阻率法二维成像技术处理后在某防空洞上得到的电阻率图像，测量时采用温纳装置，单位电极距x=0.5米，使用了60根电极。图中清晰可见防空洞的位置（高阻），与雷达探测结果以及实地调查十分吻合。

图6 某防空洞上的高密度电阻率法图像