直流电阻率三维正演的预条件共轭梯度法研究

杨金凤，邓居智，陈 辉，匡海阳

（东华理工大学 放射性地质与勘探技术国防重点学科实验室，江西抚州 344000）

直流电阻率三维正演的预条件共轭梯度法研究

杨金凤，邓居智，陈 辉，匡海阳

（东华理工大学 放射性地质与勘探技术国防重点学科实验室，江西抚州 344000）

电阻率三维正演；SSOR－PCG；异常场法；有限差分法；E－SCAN法

0 前言

近年来，随着计算机技术的快速发展，需要巨大的计算量及计算机内存的电阻率三维正演在PC机上得以实现，并逐渐步入实用阶段。但要实现快速高效的全三维电阻率反演，仍然具有一定的难度，而正演是反演的基础，因而电阻率三维正演仍然是当前研究热点和难点［2～5］。如何快速、准确地求解离散后的线性方程组，则是直流电阻率三维正演的核心问题。

在常用的线性方程组的求解方法中，预条件共轭梯度法是较好的迭代方法，特别是不完全Cholesky共轭梯度（ICCG）法［2，6～8］。但是对于大型复杂问题，必然要增加网格剖分的数量，那么计算机内存的占用量势必要增加，且计算速度也会受到影响。这就需要在此基础上研究其它一些更适合的预处理方法，进一步优化内存占用量和计算速度［9］。基于有限差分的数值计算方法，Spitzer［10］采用对称超松弛预条件共轭梯（SSOR－CG）法求解离散的线性方程组，实现了快速、精确、实用性强的正演计算。吴小平［7］引入不完全Cholesky共轭梯度（ICCG）算法及一维稀疏存储模式，实现了电阻率三维有限差分正演的快速计算。作者在SPITZER的研究基础上，利用SSOR－CG法研究直流电阻率三维有限差分正演问题，采用第一类边界条件和混合边界条件，并引入异常场法来消除总场中由点电源导致的奇异性。

1 直流电阻率法三维正演

1.1 基本方程及其边界条件

将点电源置于水平地面上A点，I为电流强度。产生电源三维地电场电位u满足微分方程：

式中 σ（x，y，z）为电导率；f为电流源项，是源位置的函数。

当电流强度为＋I的点电流源时，源函数为：

式中 （x，y，z）、（x0，y0，z0）分别为观测点和点电流源位置；δ为狄拉克函数。

u满足如下边界条件的电位函数。

其中 Γs为区域Ω的地面边界；Γ∞为区域Ω的地下无穷远边界；n是边界外法线方向的坐标变量；r为源点到边界上点的向径；θ是边界外法线方向的单位矢量n与r的夹角。

1.2 异常场法

从上面的方程可以看出，电位函数存在电流源项，即该函数在点电流源上存在奇异点。那么当求解区域存在电流源时，如果采用总场电位法求解，往往造成电源附近解的不可信。一种常见的做法是将场源所在点从求解区内“挖除”，这样便在被“挖除”的场源点周围形成新的求解边界。这些新边界上的节点离场源很近，而场源在那里产生的正常场，一般都大大超过不均匀体产生的异常场，因而这些边界节点位置可以近似为正常场值。但是由于场源附近的电场变化梯度很大，而用差商代替微商的作法是以节点之间的电位变化为基础的，所以为提高差分格式的近似程度，必须在场源附近减小网格步长。但这样会使网格剖分和构组差分格式变得复杂、繁琐，而减小步长又势必增加总结点数，从而增加计算量和计算费用［11］。

为了提高计算精度和计算速度，避免由源奇异性引起的差分计算复杂性，作者采用将解析计算和数值计算相结合的方法，即异常场法［12、13］。假设电流源附近的电导率为σM，由供电电流源在以电导率σM为均匀半空间产生的正常场为up，与供电电流源无关的由电阻率分布异常引起的异常场为us。可得如下关系式：

将式（5）、式（6）代入基本方程，得异常场方程

其中 边界条件为：

1.3 微分方程离散化

利用有限差分法对式（7）进行数值求解时，由于方程是二阶微分方程，含有对电导率空间分布的导数，要求电导率具有一定的连续性才可以求解因此对式（7）进行体积分，由格林定理将体积分化为面积分，消除对电导率的导数要求，再利用有限差分法构成求解方程［2］，即：

对式（10）应用格林定理，有：

式中 si，j，k为包围区域Δvi，j，k的表面。

对求解区域进行六面体剖分，然后在网格点上利用七点差分格式离散，如图1所示。

图1 直流电阻率三维正演网格剖分和离散示意图Fig.1 3－D finite difference grid discretization

从式（11）可见，对研究区域及其边界，边界条件式（8）和式（9）可直接应用到该方程中。利用中心差商近似得代替us/n，并沿体积Δvi，j，k的边界积分，则第（i，j，k）节点及其周围六个节点上的电位满足式（12）。

对于全部有限网格的所有节点，式（12）可表示成矩阵形式：

式中 A为耦合系数；u为待求节点的异常场列向量；b为已知向量。

2 基于SSOR－PCG法的大型稀疏线性方程组求解

对于电阻率三维正演问题，大型线性方程组的求解是决定计算速度和效率的关键一步。而由前面的分析可知，电阻率三维正演所形成的矩阵是含有大量零元素的稀疏矩阵。如果使用一维变带宽存储方式，其中仍包含大量的零元素，很大程度上浪费了内存空间；而使用直接法求解该方程组，其速度又不是很理想。为此，作者采用一维非零元素的压缩存储的方式对系数矩阵的存储，并利用预条件共轭梯度法求解大型线性方程组，极大地降低了内存占用量和提高了计算速度。

2.1 一维非零元素压缩存储模式

由系数矩阵A的特点可知，其每行或每列至多含七个非零元，即含有七个非零对角元。而从式（12）可看出，节点（i，j，k）的电位值u只与其相邻六个节点的电位值有关，因此在程序中定义七个一维整型数组ipos0（μ）、…、ipos6（μ）来指示系数C0（μ）、…、C6（μ）中的每一个非零元的列号，这里μ＝1、…、imjmkm，用以指示非零元的行号，也是节点编号。其中离散的节点按照从左到右，从前到后，从上到下的顺序编号，按编号i、j、k，i＝1、…、im；j＝1、…、jm；k＝1、…、km，首先沿剖分的x方向读取，然后再沿y方向，最后沿z方向依次循环读取。这样很容易得到任意位置编号（i，j，k）的非零元的行号，按式（15）计算：

那么对应的列号也随即可以确定。

从上面的存储过程可以看出，其比按行压缩存储更为简单易行。由于每个节点都有一个位置（i，j，k），这样就不需要保存每一个节点对应的连接系数在矩阵中的行和列。而根据节点编号可以很容易地反向计算出节点的具体位置（i，j，k），同时其周围六个节点的位置和它编号μ也可以求得。在计算矩阵与列向量的乘积时，只有七个节点相对应的列向量中非零元与列向量进行相乘，这样就可以消除与零元素相乘的运算，提高了解方程的速度，所需的存储空间也得以大大减少。

2.2 SSOR－PCG迭代算法

共轭梯度迭代算法，对于系数矩阵接近单位阵时其对应条件数很小，计算速度很快；而当A的条件数cond（A）＞102时，则收敛速度就很慢。所以对于电阻率三维正演计算的类似问题，由于网格大小和物性参数可能相差几个数量级，所形成的大型稀疏线性方程组的稀疏矩阵A的特征值的变化范围很大，因而对应的条件数必然很大，所以就要对式（14）的线性系统进行预处理，改善矩阵A的条件数。然后再运用CG法解该式，计算速度将会大大提高。

作者在本文采用SSOR预条件方法［14］，先将系数矩阵A分解为：

其中 E是下三角矩阵；I是单位矩阵。

其中 L＝I＋wE；w是松弛因子。

解式（14）的SSOR预条件共轭梯度迭代过程如下：

其中 r0、p0是共轭梯度迭代的初始化向量；u0是解的初值；α、β为常数。

从以上的计算流程可以看出，预条件矩阵很容易求得，且不增加计算内存需求。这样在整个计算过程中，主要就是系数矩阵A与任一列向量的乘积，以及预条件矩阵与任一列向量的乘积。而A中每一行的非零元素不参加计算，再加上用一维稀疏存储，所以可以快速求得。关键是如何快速的求得预条件矩阵M的逆与任一向量的乘积，直接求解很困难。注意到L是下三角阵，而LΤ是上三角阵，可先将ti＝（LLΤ）－1ri转化为求方程组LLΤti＝ri，然后用高斯消去法先顺带、后回带解得ti。由于矩阵L中仅有少量的非零元素，并以稀疏存储方式存储，因此顺带和回带过程仅涉及少量非零元素的操作，均可以非常简单、快速地完成。

松弛因子w的选择是控制迭代收敛的重要参数，对模型给予不同的w值试算结果如下页图2所示。从图2可看出，w在1.3～1.7之间迭代次数变化较小，对迭代过程影响不是很大；而w取“0”时，该迭代是发散的；最终取w值为“1.5”，此时的收敛速度最快。SSOR－PCG迭代的收敛标准是是第一次迭代后残差的L2范数。

由此可见，单独的SSOR－PCG迭代技术或系数矩阵稀疏存储技术均不能有效的解决问题，但将两者有机地结合，再加上合适的松弛因子的选择就可形成快速高效的迭代算法。

3 模型算例

3.1 垂直接触带模型

垂直接触带模型如下页图3（b）所示，左边介质的电阻率是1Ωm，右边的电阻率为10Ωm。单

其预条件矩阵定义为：位点电源位于左边介质上，距界面6m，采用39× 39×20网格剖分，计算机配置为Intel Core 3.0GHz双核处理器，内存3.0GB，32为XP系统。SSOR－PCG算法迭代546次收敛，耗时3.69s。图3（a）是单位点电源的总电场法，异常场法计算结果和理论电位曲线的比较。异常电场法计算的电位曲线和理论电位曲线几乎完全吻合，平均误差为0.69%；而总场法计算的电位曲线相差很大，其平均误差为9.03%，尤其是在源点附近误差很大，这是由源的奇异性引起的。

图2 对不同的松弛因子w值SSOR－PCG算法所需的迭代次数Fig.2 The number of iteration required by SSOR－PCG for different relaxation factors

3.2 层状介质模型

模型2为二层介质，第一层电阻率为10Ωm，第二层的电阻率为100Ωm，厚度为3m。采用39 ×39×20的网格剖分，SSOR－PCG迭代742次，耗时5.23s。

表1是单位点电源的总场法、异常场位法计算结果ut、us和理论电位值u的比较。异常场法计算的平均误差△uS为0.16%，而总场法的平均误差△ut为6.48%，精度提高一个多量级。特别是在电源附近，总场法计算误差较大，异常场法计算精度较高。从表1还可看出，利用源点移除法，去除源点后总场法的平均误差为3.45%，误差减小了近一半，但还是比异常场法的平均误差高，其平均误差只有0.12%。可见异常场法对于提高计算精度，有显著的效果。

3.3 低阻异常模型

为了说明SSOR－PCG迭代算法在直流电阻率三维正演模拟中计算的可靠性，作者设计了一个典型的低阻异常模型，设地下有一个40m×40m× 20m，电阻率为1Ωm的低阻长方体，如图4所示其顶部埋深为20m，围岩的电阻率为100Ωm；图4（b）是电性分界线在地面（xoy平面）的投影。将计算区域剖分为37×37×21个单元。

作者采用的是高密度三维观测方法ESCAN［15］，将所布设的电极网格化为21×21，其中x方向和y方向均布设21个电极，点距、线距均为5m。图4（b）标出地面上E－SCAN法测量的441个电极位置。

图5为y＝0低阻异常体中心E－SCAN装置视电阻率拟断面图，其中横坐标x为算术坐标，纵坐标AM/2为对数坐标。由图5可见，低阻模型会产生一个低阻异常区和高阻异常区，低阻异常区范围为AM/2从2.5m至20m，x方向从－40m～40m变化。高阻异常区位于低阻异常体下方，横向范围与模型基本一致。

图6为低阻异常体E－SCAN装置在AM/2分别取3.54m、10.31m、15.2m、27.5m时的视电阻率等值线图。图6（a）中由于其探测深度还未到低阻异常体，视电阻率基本无变化；图6（b）及图6（c），由于受低阻异常体影响，在中心位置有一个明显的低阻异常区，范围比低阻异常体大；下页图6（d）其探测深度穿过异常体时，在中心位置出现一个明显半径大约为20m近似圆形高阻假异常区由此可见，E－SCAN装置视电阻率受低阻异常体影响，随着AM的增大由低阻异常转向高阻异常异常范围也随之变小，并且E－SCAN装置对异常体横向分辨较好，纵向分辨率差。

4 结论

作者在本文详细推导了异常场法所形成的方程，并对方程进行差分离散得到大型稀疏线性方程组，将求解该大型稀疏线性方程组的SSOR－PCG迭代算法、系数矩阵的稀疏存储技术，以及消除点电源奇异性的异常场法三者有机的结合在一起。通过对二个简单模型的模拟结果与解析解对比，结果表明该方案的数值解与理论解基本吻合，说明将解析计算和数值计算相结合，并利用该算法求该正演计算中的异常电位，能消除点电源奇异性影响，在保持快速计算的基础上，显著提高了计算精度，而且内存需求大大减少。另外对低阻地电模型的正演模拟结果表明，SSOR－PCG是电阻率三维有限差分稳定高效的计算方法，实现了快速精确的电阻率三维正演计算，为其三维反演奠定了基础。

［1］ 傅良魁.电法勘探教程［M］.北京：地质出版社，1983.

［2］ DEY A，MORRISON H F.Resistivity modeling for arbitrarily shaped three－dimensional structures［J］.Geophysics，1979，44（4）：615.

［3］ MA Q Z.The boundary elemental method for 3－D DC resistivity modeling in layered earth［J］.Geophysics，1998，63（2）：399.

［4］ 阮百尧，熊彬，徐世浙.三维地电断面电阻率测深有限元数值模拟［J］.地球科学，2001，26（1）：73.

［5］ 鲁晶津，吴小平，KLAUS SPITZER.直流电阻率三维正演的代数多重网格方法［J］.地球物理学报，2010，53（3）：700.

［6］ LI Y G，SPITZER K.Three－dimensional DC resistivity forward modeling using finite elemental in comparision with difference solutions［J］.Geophys.J.Int.，2002，151（3）：924.

［7］ 吴小平，徐果明，李时灿.利用不完全Cholesky共轭梯度法求解电源三维地电场［J］.地球物理学报1998，41（6）：838.

［8］ WU X P.A 3－D finite－element algorithm for DC resistivity modeling using the shifted incomplet Cholesky conjugate gradient method［J］.Geophys.J Int.，2003，154（3）：947.

［9］ 刘斌，李术才，李树忱，等.基于预条件共轭梯度法的直流电阻率三维有限元正演研究［J］.岩土工程学报2010，32（12）：1846.

［10］SPITZER K.A 3－D finite－difference algorithm fo DC resistivity modeling using conjugate gradien methods［J］.Geophys.J.Int.，1995（123）：903.

［11］罗延钟，张桂青.电子计算机在电法勘探中的应用［M］.武汉：武汉地质学院出版社，1987.

［12］LOWRY T，ALLEN M B，SHIVER P N.Singularity removal：a refinement of resistivity modeling techniques［J］.Geophysics，1989，54（6）：766.

［13］ZHAO S，YEDLIN M J.Some refinements on the finite－difference method for 3－d dc resistivity modeling［J］.Geophysics，1996，61（5）：1301.

［14］吴小平，汪彤彤.利用共轭梯度算法的电阻率三维有限元正演［J］.地球物理学报，2003，46（3）：428.

［15］黄俊革，王家林，阮百尧.三维高密度电阻率ESCAN法有限元模拟异常特征研究［J］.地球物理学报，2006，49（4）：1206.

book=147,ebook=147

1001—1749（2012）03—0303—07

P 631.3＋22

A

10.3969/j.issn.1001－1749.2012.03.11

杨金凤（1984－），女，硕士，主要从事电法勘探正反演研究。

国家科技支撑计划（2009BAB43B03；2011BAB04B03）

2012－01－05改回日期：2012－01－09