1 引 言

直流电阻率勘探以其成像效果好、适应性强、方便快捷等诸多优点被广泛应用于岩土工程、矿产资源、环境等领域的勘查工作中.近年来该方法在观测空间、应用功能、观测尺度等方面均取得了新的拓展.在观测空间方面，直流电阻率法勘探由地表无限半空间发展到矿井、隧道、钻孔等三维全空间中，面对的勘探对象更加复杂；在应用功能方面，直流电阻率法被引入到突水、滑坡等地质灾害实时监测、地震前兆监测等新领域，与传统的勘探功能不同，实时监测功能要求电阻率法勘探必须具备快速的实时反演成像方法；而在观测尺度方面，由野外大尺度勘探拓展到实验室尺度下的精细探查(如岩石试块的破裂监测等)，对反演精度要求较高.综上所述，随着直流电阻率法应用范围的拓展，对其反演精度、反演速度等的要求越来越高，在反演精度和反演速度等方面对三维电阻率反演方法进行深入的研究具有重要的实际意义和理论价值.

目前三维电阻率反演方法面临着多解性、不适定性及求解速度慢等问题和困难.国内外许多学者对此问题进行了深入的研究^[1-5]，用于三维反演的方法有α中心方法、Born 近似理论、Zohdy反演法、最大似然反演理论和最小二乘反演理论等.其中，Sasaki将有限单元法和最小二乘法相结合，在最小二乘准则中加入光滑约束，系统地研究了敏感度矩阵、光滑度矩阵的计算方法，为使用有限单元法进行电阻率反演奠定了基础^[2].此后基于光滑约束的最小二乘线性反演方法成为电阻率反演的主要方法.

阮百尧、黄俊革等利用电极互换定律法通过一次正演得到敏感度矩阵，将体积因子作为光滑约束，使得三维电阻率反演的精度和深部网格的分辨率得到显著提高，但是其反演速度较低^{[6, 7]}.底青云等从微分方程的积分解出发推导出偏导数矩阵的解析表达式，实现了敏感度矩阵的快速计算^[8].JieZhang，吴小平等将共轭梯度方法应用于E-SCAN 测量方式的正演模拟与反演成像，其突出之处是避免了敏感度矩阵的直接计算，但对于E-SCAN 方式每次反演中实际上需要多次正演才能求得敏感度矩阵，耗时较长，且由于没有进行敏感度矩阵的直接求解，使得无法检验敏感度矩阵的正确性^{[9, 10]}.宛新林等通过对光滑度矩阵元素进行适当改进，使之适用于各种情况下光滑度矩阵的求取，并用最小二乘正交分解法求解反演方程，有效地提高了计算速度^[11].AdamPidlisecky等开发了基于Matlab 语言的三维电阻反演程序，使用预条件双共轭梯度法求解正演中的大型线性方程组，使用不精确高斯牛顿法作为反演算法，取得了较好的效果^[12].

上述研究工作为三维电阻率反演奠定了基础，但是三维电阻率反演还存在一些未能解决的问题，主要有以下两个方面：(1)非惟一性问题是三维电阻率反演的固有问题，有时表现得较为严重，除了光滑约束之外，尚未提出其他约束方法，急需提出合理的约束方法以便有效地改善反演问题的非惟一性和反演效果；(2)反演耗时较长，成为制约三维电阻率反演在实际应用中推广的瓶颈，有待对反演计算效率进行优化，以便提高反演速度.

2 基于不等式约束与光滑约束的最小二乘反演方法

本文采用以三维有限元数值正演为基础的反演方法，设三维有限元模型的网格数量为n_x×n_y×n_z= m个(其中n_x，n_y，n_z分别为x，y，z方向上的网格数)，设模型参数为m= (ρ₁，ρ₂，…，ρ_m)^T，视电阻率观测数据为d_obs = (ρ_s1，ρ_s2，…，ρsn)^T，其中n为观测数据的数量.

三维电阻率反演是非线性问题，将其线性化并正则化，得到如下方程：

(1)

式中，A为敏感度矩阵，Δ_m为模型参数增量向量，Δd为观测数据d_obs与正演理论值d_m 的残差向量，d_m 是根据给定的模型参数由数值正演得到的理论观测数据.

2.1 光滑约束的施加

三维电阻率反演问题往往表现为混定问题，往往导致方程(1)为病态方程.为解决该问题，将光滑约束引入反演方程.所谓光滑约束，就是指相邻网格的电阻率光滑过渡，也就是使相邻网格电阻率差异极小.对于第i个网格而言，光滑约束可表示如下：

(2)

式中，Δm_i为第i个网格的模型参数修正量，Δm_i^F，Δm_i^B，Δm_i^L，Δm_i^R，Δm_i^U，Δm_i^D分别表示第i个网格的前后、左右和上下的网格的电阻率修正量.可将整个模型的光滑约束用矩阵形式表示：

(3)

式中C为光滑度矩阵.施加光滑约束后的三维电阻率反演方程如下：

(4)

式中λ 为拉格朗日常数，λ 的大小决定了光滑约束的权重.

2.2 不等式约束的施加

多解性问题是地球物理反演的固有问题^[13-14]，为进一步减小反演问题的多解性，将模型参数的变化范围作为先验信息，将表征模型参数变化范围的不等式约束

(5)

施加到反演方程中.式中，m_i为第i个网格的电阻率，ρmin_i和ρmax_i分别为第i个网格的电阻率的下限和上限.需要指出的是，电阻率的变化范围可以是根据一般物理常识获得的较为宽泛的一个范围，也可以是根据钻孔等其他方式获得的较为精确的一个范围.

综合考虑光滑约束和不等式约束，提出如下三维电阻率反演的目标函数：

(6)

在式(6)中,光滑约束可按照式(4)的形式施加到反演方程中,而如何将不等式约束施加到反演方程中是本文的一个关键性难题.为此采用H.J.Kim的做法^[15],定义一个新的参数向量X来表征模型参数,X与m的关系表示为

(7)

式中X_i为向量X的元素.

根据式(7)可得到X的扰动量ΔX与Δm之间的关系：

(8)

得到ΔX_i后便可得到下一代模型参数m^{(k+ 1)},如式(9).由式(9)可以看出第i个网格模型参数被严格限制在ρ_mini~ρ_maxi之间,通过这种方式将不等式约束施加到了反演方程中.式中,m^_i(k+1)为第k+1代的第i个模型参数,

(8)

2.3 三维电阻率反演流程

本文所采用的最小二乘法三维电阻率反演的流程见图 1，具体如下：

(1) 首先设定网格电阻率的初值，一般设定为与观测数据背景值接近的均一初始模型；

(2) 通过数值正演得到相应的理论观测数据d_m；

(3) 进行反演收敛判断，若理论观测数据与实际观测数据之间的误差满足收敛判据，将此时得到的模型参数作为反演结果输出.反之进行下一步计算；

(4) 计算敏感度矩阵A和光滑度矩阵C，求解反演方程，得到ΔX；(5)计算得到新一代模型参数，执行第(2)步，进入下一循环.

需要指出的是，反演收敛的判据为rus<ε_inv，其中rus为观测数据d_obs 与正演理论值d_m 之间的均方误差，ε_inv 为反演收敛的容许值.

(10)

3 三维电阻率反演计算效率优化

在整个反演流程中，最为耗时的是敏感度矩阵A的计算和反演方程的求解，提高二者的求解速度对于优化三维电阻率反演的计算效率至关重要.因此提出了三维电阻率反演计算效率优化方案，在该方案中，Cholesky分解法被用来求解敏感度矩阵计算中的多个点源场的正演问题，该方法只需对总体系数矩阵进行一次分解，然后对不同的右端向量进行回代即可；预条件共轭梯度法被引入到反演方程的求解中，将雅可比迭代中的对角阵作为预处理矩阵，具有求逆方便、无需内存空间的特点，有效地加快了收敛速度.需要说明的是，因为将雅可比迭代的块对角矩阵作为预条件矩阵，所以将其称为雅可比预条件共轭梯度法(Jacobi Preconditioning Conjugate Gradient，JPCG).

3.1 基于CholesKy分解算法的敏感度矩阵求取

敏感度矩阵A表示的是观测数据对模型参数的偏导数矩阵，即观测参数ρsi随网格电阻率ρ_j的一次变化率表示为

(11)

本文采用电极互换定律法求解敏感度矩阵，其基本思路见文献[7].该方法需要分别将每个测量点作为点电源，分别计算每个点源供电时的电场分布，也就是说每求解一次敏感度矩阵需要进行n次点源电场的正演.方程

(12)

为三维点源电场有限元正演中的大型线性方程组，式中：u和u₀ 分别为含有各节点异常和正常电位值的向量，K′ 和K分别为正常和异常电位向量的总体系数矩阵.对于一次敏感度矩阵求取过程中的n次点源三维场正演而言，K保持不变，而右端向量随点源的变化而变化.

对于大型线性方程组求解而言，PCG 法的求解速度一般要远高于Cholesky分解法.但对于上述情况，Cholesky分解法往往优于PCG 算法.因为针对n次正演，Cholesky 法只需对K进行一次分解，然后对不同的右端向量进行回代计算即可，而PCG 法需要对每次正演的方程组进行单独求解，当n取值较大时，PCG 法的求解速度明显慢于Cholesky 分解法.

为了对比Cholesky 分解法和JPCG 法在求解敏感度矩阵中的计算效率，假设模型单元数为43600个，观测数据n=80，进行80次三维点源场正演的耗时如表 1.可见，对于单个线性方程组(一次正演)JPCG 法具显著的速度优势，而对于较多点电源情况下的敏感度矩阵求解，Cholesky分解法具有较为明显的优势.因此本文采用Cholesky分解法来求解敏感度矩阵.

3.2 基于jJPCG算法的反演方程求解^{[16, 17]}

由表 1可看出，对于只有单个右端向量的线性方程组，JPCG 法的求解速度要远高于Cholesky分解法，因此本文采用JPCG 法求解反演方程(4).

PCG 法通过引入预条件矩阵M达到了降低系数矩阵条件数的目的，有效地加速了收敛.对于反演方程，具体的PCG法求解流程如图 2，其中B为方程的右端向量，r是方程的右端与方程左端的差向量，ε 为算法收敛的容许误差，z，p，α，β 为迭代过程的中间参数.在PCG 法迭代过程中，矩阵乘A^TA，C^TC和方程组Nz(i+1)=r(i+1) 的求解是决定收敛速度的关键.

3.2.1 矩阵乘的处理

若直接计算A^TA和C^TC，耗时较长，内存空间消耗较大.因为PCG 法只需要计算A^TA或C^TC与某个列向量的乘积，可利用矩阵乘的交换定律来计算：

(13)

式中b为一个列向量，y₁ =Ab，y₂ =Cb.其计算量比直接执行矩阵乘的运算大大减小.

3.2.2 预条件矩阵的选取

在PCG 算法中，解方程组Mz⁽ⁱ⁺¹⁾=r⁽ⁱ⁺¹⁾ 是较为关键的步骤，本文采用雅可比迭代的块对角矩阵作为预处理矩阵.对于反演方程(4)，预条件矩阵的具体形式为

(14)

由于该预条件矩阵是对角矩阵，故求逆非常容易，大大提高了计算速度.

4 反演实例 4.1 合成数据的反演算例

为了评价最小二乘法三维电阻率反演的效果，利用合成观测数据进行三维反演.原模型如图 3，在均匀半空间中存在两个电阻率为10Ωm 的低阻体，围岩的电阻率为300Ωm，其中低阻体一的顶部埋深为2m，低阻体二的顶部埋深为6m.采用施伦贝谢尔装置测量，共30根电极，电极间距为2m，共布置了5 条测线(其中测线一位于x=0位置，其余测线之间的间距为2m).

利用5条测线的合成观测数据进行反演，经过对视电阻率剖面图的分析对不等式约束进行了定义，因为观测数据中的视电阻率基本在300Ωm 以内，因此将整个模型中网格的电阻率值上限设为350Ωm(有一定的冗余)，下限设为0Ωm.而在5.0m<z<10.0m 与22.0m<z<28.0m 范围内视电阻率数值稍低，基本在240Ωm 以内，因此将该范围内网格电阻率上限值设为300Ωm，下限值设为0Ωm.

建立三维模型，模型中分为目标区和非目标区，目标区是参与反演的区域，也是数据观测区域，网格的步长均匀，用均匀网格剖分.非目标区是边界区域，网格的步长倍增，以模拟无限远边界.该算例中目标区域的网格数量为15×15×40=9000个，模型的反演迭代次数为10 次，耗时约160 min(计算机配置为主频1.83GHz，内存1.0GB).图 4为施加不等式约束后的反演结果，可见两个低阻体在位置、规模、形状、电阻率值等特征方面均与原模型基本一致，将反演结果中电阻率低于60Ωm 的网格提取出来，发现提取出的常体与原模型中的低阻体完全一致.

为了对比验证，同时进行了传统的最小二乘反演，未施加不等式约束，反演结果如图 5，可见反演效果明显降低，尤其对于低阻体二成像精度较低，成像效果较差.

可见，虽然本算例中在施加不等式约束时限定的电阻率变化范围较宽泛，但这种不等式约束对反演过程中搜索方向的引导作用非常明显，借助不等式约束和反演效率优化方案，最小二乘反演方法可得到较为精确的反演结果，且反演计算速度也满足勘探工作的要求.

4.2 试验数据的反演算例

为了进一步评价最小二乘法三维电阻率反演的效果，将电阻率层析成像系统用来监测煤矿突水模型试验中岩层断裂以及突水通道形成的过程，并对试验数据进行了反演计算.

4.2.1 试验基本情况

物理模型试验的示意图(如图 6)给出了岩层的组成以及含水构造、监测层的位置.其中防突层上方的含水层中预置了含水构造，以模拟作用在防突层上的水压力.防突层自上至下由软岩层(厚8cm)、监测层(厚1cm)、硬岩层(厚5cm)和软岩层(厚2cm)组成^[18].试验材料为新型流固耦合试验材料，由于该试验材料不导电，故在防突层中设置了电阻率层析成像监测层.监测层中的材料采用黏土等导电材料.由于监测层位于防突层之中，故监测层中电阻率信息的变化可直观地反映采动应力和水压力作用下防突层的变化信息.监测层的位置和测线布置见图 6，电极总数为30 根，采用施伦贝谢尔装置形式，每组成像数据包含196个视电阻率值.

4.2.2 试验数据反演

笔者在文献[19]中对试验数据进行了细致的分析，但存在一个问题，那就是反演得到的裂隙分布与实际试验情况有所偏差.为了解决上述问题，本文采用基于不等式约束和光滑约束的最小二乘反演方法对试验数据进行反演计算，得到了突水灾害发生之前02:35:47时刻监测层的电阻率结构，如图 7.

经过对视电阻率剖面图的分析，对不等式约束进行了设定，整个模型中网格的电阻率值上限为400Ωm，下限为0Ωm，而在0.25m<z<0.50m 与0.60m<z<0.90 m 范围的网格电阻率上限值为800Ωm，下限值为0Ωm.

三维反演模型在x，y，z方向上的网格数分别为70×2×20=2800 个，达到收敛所需的迭代次数为7次，总耗时为520s(计算机主频为1.83GHz，内存为1.0GB)，可见上述反演方法具有迭代次数少，收敛速度快的优点.

4.2.3 反演结果与实际情况对比

02:35:47 时刻监测层电阻率结构图表明，在之前裂隙发展的基础上，中部的高阻裂隙几乎贯通，而左侧的高阻裂隙的扩展速度亦加快.根据裂隙的状态可以推断防突层将在短时间内于中部发生断裂坍塌.而试验过程中，防突层(含监测层)中部于02:38:24时刻发生断裂，最后防突层的左侧发生断裂，图 8是防突层发生坍塌突水后的照片，可见中部和左部裂隙的位置与反演结果中的高阻裂隙的位置基本一致.与文献[19]中相应反演结果对比，可发现施加不等式约束后，反演结果更接近最优解，与实际情况基本一致，反演的多解性显著降低.

上述分析表明以本文提出的三维电阻率反演方法为核心技术的电阻率层析成像监测系统可将裂隙扩展、岩层破断等重要前兆过程以成像的方式直观形象表达出来，有效地捕捉到突水前兆信息，为突涌水灾害的预测提供了重要的参考和指导.

5 结 论

(1) 本文将介质电阻率取值范围作为先验信息和约束条件，将表征介质电阻率取值范围的不等式约束施加到三维电阻率反演方程中，该约束从理论上讲有利于降低反演问题的多解性，可有效地改善反演效果.

(2) 本文基于Cholesky 分解算法和预条件共轭梯度算法(PCG)的特点提出了反演计算效率优化方案，实现了三维电阻率观测数据的快速稳定反演.

(3) 反演算例表明，借助不等式约束和反演效率优化方案，最小二乘反演方法可得到较为精确的反演结果，有效地提高了反演计算效率，具有良好的推广前景.

(4) 对煤矿采区模型试验中电阻率层析成像监测试验数据的反演实例表明，本文的三维电阻率反演方法可以满足实时监测工作中电阻率勘探数据的快速反演的要求，可实现实验室尺度下的介质电阻率结构的精细探查.本文三维电阻率反演方法的上述特点有利于直流电阻率勘探方法在地质灾害实时监测这一新领域的推广.

(5) 本文中所提出的三维电阻率反演方法中不等式约束的施加是关键，针对某一具体的观测数据而言，可以根据一般物理常识或者钻孔等方式获得不等式约束范围，但通常情况下可以根据视电阻率数据的特点来确定不等式约束，当然这需要建立在对多种地电结构进行深入正演研究的基础上，这是本文的难点.