2010, Vol. 53
2. 桂林理工大学地球科学学院, 桂林 541004
2. College of Earth Sciences, Guilin University of Technology, Guilin 541004, China
近些年来，随着固体矿产资源的供需矛盾加剧，人们急需开拓“第二找矿空间”，即寻找深矿、盲矿，这必然加大地球物理找矿的难度.电法是非常有效的找矿方法，但是传统的电法勘探随着目标体的埋深增加，分辨率随之降低.基于井中供电地表测量电位分布的井地电法，提高了对深部地质体的勘探能力，在深部探矿中具有很好的应用前景.对井地电法进行三维正演数值模拟以及开发相关的反演技术，具有很好的理论和实际应用价值，目前一些国内外学者开始关注这一工作并进行了一些相关的研究，如Wang等[1]，Avdeev等[2]分别用有限差分和积分方程法对电阻率测井进行了三维模拟，Nam等[3]针对层状介质的情况，用傅里叶展开和非正交坐标系实现了侧向测井的快速三维有限元模拟，岳建华等[4]利用有限差分法实现了井地电法的三维数值模拟，王志刚等[5]，徐凯军等[6]，柯敢攀等[7]利用有限差分法实现了垂直有限线源井地电法的三维正演模拟，取得了较好的研究成果.
由于实际矿藏勘探过程中需要面对复杂的地形地貌和复杂的深部结构特征，在正演模拟中如何更好地选择剖分区域，对起伏地形和倾斜井的情形做适当处理等均需要进一步地展开.此外，井中电法在深度方向上要求大范围的剖分，从而加大了数值模拟的计算量，这使得现有的正演计算效率不足以满足野外工作的需求.有限元法的单元划分灵活，适用于复杂地电模型的数值模拟，在地球物理的数值模拟中得到了广泛的应用，如阮百尧等[8, 9]用有限元法实现了电阻率测深的三维有限元模拟.本文基于这些问题，采用有限单元法，考虑到地形起伏和倾斜井的情形，以井为中心选择与井轴一致倾斜的计算区域，将其离散成三个侧面为平行四边形，上下面随地形起伏的平行三棱柱单元，利用仿射变换技术，定义了从子单元到母单元的坐标映射，给出了精确的单元积分，用压缩存储方式存储刚度矩阵，用预处理-共轭梯度法来求解线性方程组，提高了计算效率，为反演效率的提高奠定了基础.
2 方法原理 2.1 地下点源场的边值问题和变分问题地下三维点源场边值问题可表述为[10];
|
(1) |
等价的变分问题是[10];
|
(2) |
其中，u代表电位，I是供电电流强度，σ为介质的电导率，r为电源点A到测点的距离，n是边界的外法向向量，Ω是研究区域，Γs为研究区域的地面边界，Γ∞为研究区域的地下边界.
2.2 区域离散与单元分析为了处理复杂地形，我们将地形起伏转化为z坐标的变化，使向下方向的网格剖分随地形起伏，对网格节点的z坐标作出修正; z=z-p，p是相应的地面点的地形高程.正地形时p＞0，水平时p=0，负地形时p＜0.
为了能处理地形起伏和井倾斜的情况，同时考虑到用尽可能少的网格达到既满足计算精度又提高计算效率的目的，我们以井为中心采用矩形区域剖分，远离中心区域网格间距逐渐增大，整个剖分区域与井保持相同的角度倾斜，将区域离散成平行五面体单元，见图 1.网格单元的上下底为相互平行的三角形，三个侧面为平行四边形，是平行五面体单元，因此可以采用仿射变换将子单元映射到母单元.记直角坐标系为Ⅰ[O；e1，e2，e3]，母单元所在的局部坐标系为Ⅱ[O′；e′1，e′2，e′3]，子单元到母单元的变换如图 2所示.
|
图 1 网格剖分示意图 Fig. 1 Schematic diagram for finite element mesh |
|
图 2 平行五面体单元的变换 Fig. 2 Coordinate transformation for parallel tri-prism elements |
五面体的六个顶点记为; pi（xi，yi，zi）（i=1，2，…，6），定义仿射坐标变换;
|
则p=Ap′+p1，从而
|
其中，det（A）是矩阵A的行列式，A-1代表A的逆矩阵，bij是A中元素aji的代数余子式.
五面体单元可以看作是六面体单元的退化，单元中电位采用三线性插值，在母单元上选取插值函数为[11];
|
则
式中ui为五面体第i个顶点的电位值.
记
|
则有u=NTue，其中，上标T表示取矩阵的转置.
将（2）式中右端各项离散成单元积分之和，单元的电位列向量ue扩展成全体节点的电位列向量u，将6×6的单元刚度矩阵扩展成单元系数矩阵，全部单元系数矩阵相加，得到
|
其中，b=（0，…，0，I，0，…，0）T，即b是一个与电源点相应的位置取电流值，其他位置值为0的向量.由δF（u）=0得到线性方程组;
|
在单元集成过程中，用一个二维数组ijk存储单元的节点信息，用一个一维数组存储边界信息，在单元积分时，在含有边界面的单元上可同时加上相应的边界积分.由于总刚度矩阵K为大型的稀疏正定矩阵，为了节约存储空间和提高计算效率，我们采用CSR（行压缩存储）方式对其进行存储.具体而言，将总刚度矩阵存储于三个一维数组（a，ja，ia）中，通过读取ijk数组的信息，可判断每个节点属于哪几个单元，从而决定三个一维数组（a，ja，ia）的结构，实现CSR存储.
2.3 线性方程组的求解由于经过有限元离散而生成的刚度矩阵是对称正定阵，我们的正演问题也就转化为求解一个大型稀疏对称正定线性系统.共轭梯度法是求解此类问题的一种高效迭代类算法，该迭代算法收敛的快慢依赖于系数矩阵的谱分布情况; 当矩阵的特征值分布于一个很长的区间上时，矩阵的条件数变大，收敛速度可能会变得很慢.人们通常采用预处理技术来提高迭代的收敛速度，通常的预条件子有不完全LU分解预处理阵，多项式预处理矩阵等.经过验证，本文涉及的问题所生成的刚度矩阵的不完全Cholesky分解是不稳定，为此本文采用了SGS（SymmetricGauss-Seidel）预条件子，将K分裂为K=D-E-ET，其中对角阵D=diag（K），-E是K的严格下三角部分，从而
|
这里，
|
由于在生成M过程中及用前代和回代求解Mzi-1=ri-1时，需要更频繁地读取对角元素，且对角元均为非零元，为便于存取和提高计算效率，L，U按照MSR（改进的行压缩存储）格式存储，将对角元素取逆，存储在一维数组的最前面.这样本文使用的预处理-共轭梯度法（PCG）[12~14]伪码如下;
将K分裂为L，U；M=LU
给定初值x0，计算r0=b-Kx0
对i=1，2，…循环
求解Mzi-1=ri-1
ρi-1=ri-1Tzi-1
判断，若i=1执行
ρ1=z0
否则，执行
βi-1=ρi-1/ρi-2
pi=zi-1+βi-1pi-1
结束判断
qi=Kpi
αi=ρi-1/pi-1Tqi
xi=xi-1+αipi
ri=ri-1-αiqi
检查收敛性，不满足收敛要求则继续循环结束循环
3 算例依据上述原理和算法，我们编制了可以处理复杂地形和倾斜井情形的井-地电法三维正演模拟程序.程序在Intel（R）Core（TM）2DuoCPU，2.66GHz，3.00GB内存，WindowsXP下运行.以下面的模型4为例，网格剖分60×50×50，单元数为300000，节点数为158661，做仿射坐标变换，将子单元映射到母单元，可直接求出单元积分，则单元刚度矩阵计算及总刚度矩阵集成只需2.45s，而采用等参变换，做高斯数值积分来求解单元刚度矩阵，这个过程却高达2024s，在精度范围内二者得到了完全相同的刚度矩阵，用预处理-共轭梯度法求解158661阶压缩存储后的线性方程组，需要16.3s，做一次完整的正演计算，包括读取模型等，总耗时为32.26s.
模型1; 三层层状大地，第一层电阻率50Ωm，厚度为5m；第二层电阻率100Ωm，厚度为10m；第三层电阻率20Ωm.供电电极在地面，测量装置采用对称四极，解析解采用数字滤波法计算，数值解与解析解对比最大相对误差为4.86%，见图 3.
|
图 3 数值解与解析解对比图 Fig. 3 Comparison between analytical and numerical solutions |
模型2; 如图 4所示，为一90°角域二维山脊地形，均匀介质，电阻率为10Ωm，电源点在地下，位于x=-3m处，在图中用A（+I）表示.表 1为地表电位数值解与解析解的对比，最大相对误差为2.9%，平均相对误差为0.66%.
|
图 4 90°角域山脊地形 Fig. 4 90° ridge topography mode |
|
表 1 90°角域地形数值解与解析解对比 Table 1 Comparison between analytical and numerical solutions of a 90° ridge topography model |
模型3; 模型如图 5所示，异常体为一长10m、宽4m、高2m的长方体，走向为y方向，顶部埋深6m，电阻率为10Ωm，背景电阻率为200Ωm，井直立，穿过异常体中心，图中用A1，A2，A3，A4来标示供电点.图 7为供电点置于井内深度h=3m、h=6m、h=8m、h=12m时地表电位异常等值线图，虚线框标出异常体在地面的投影.由图中可以看出，当供电点在异常体上方时出现负异常；当供电点在异常体下方，在低阻体边界处出现两个正异常极值，可以清楚地看出低阻异常体的位置.
|
图 5 模型3示意图 Fig. 5 Schematic diagram of model 3 |
|
图 6 模型4示意图 Fig. 6 Schematic diagram of model 4 |
|
图 7 直立井地面电位异常等值线图（等值线中数字单位为V） Fig. 7 Ground potential anomaly contour section for the vertical well model (unit: V) (a)h=3m; (b)h=6m; (c)h=8m; (d)h=12m. |
模型4; 模型如图 6所示，异常体同模型3，顶部埋深6m，异常体中心坐标为（-3.5，0，7），电阻率为10Ωm，围岩电阻率200Ωm.设井轴位于xoz平面上，与z轴夹角为30°，斜穿过异常体.图 8分别为供电点置于井中深度h=3.5m、h=6m、h=8m、h=11m时地面电位异常等值线图，这里h指的是从地面到供电点的垂直距离.虚线框标出异常体在地面的投影.从图中可以看出，电位异常等值线在x方向不再对称.由于低阻体吸引电流的作用，当供电点位于异常体下方靠近异常体时，靠近点源方向出现负异常，在离开点源方向，靠近异常体方向出现正异常极值，两个正异常极值出现在异常体的两端.随着供电点埋深增加，供电点远离异常体，正异常极值由两个变为一个，沿低阻体走向拉长.
|
图 8 倾斜井地面电位异常等值线图（等值线上数字单位为V） Fig. 8 Ground potential anomaly contour section for the deviated well model (unit: V) (a)h=3.5m; (b)h=6m; (c)h=8m; (d)h=11m. |
模型5; 模型在x方向有地形起伏.如图 9所示，图中用粗线和十字标记分别标识井和供电点的位置，井的方向同模型4，井轴与z轴夹角为30°，背景电阻率为20Ωm.图 9a为均匀大地情形下xz断面的电位等值线图；图 9b为含有异常体情形下xz断面的电位等值线图，其中左上的虚线方框所在位置是一4m×4m×4m低阻体，电阻率为1Ωm，右下的虚线方框所在位置是4m×5m×6m的高阻体，电阻率500Ωm.由图中可以看出，与均匀大地的情形相比，异常体的存在使电位等值线的分布发生改变，低阻体所在位置与背景区域相比是电位分布低梯度区，低阻体上的电位梯度小，等值线稀疏；反之，高阻体上的电位梯度大，等值线密集.
|
图 9 起伏地形xz断面电位等值线图（等值线上数字单位为V） (a)均匀大地; (b)含低阻和高阻异常体. Fig. 9 xzpotential contour section for the undulating surface topography model (unit: V) (a) Homogeneous; (b) Embedded low and high resistivity blocks. |
为了实现起伏地形、倾斜井情况下点电源井地电法的正演，本文基于仿射坐标变换、压缩存储和预条件krylov子空间技术等研究了三维有限元数值模拟算法，并编制了相应的程序，利用五个数值模型对程序的可靠性和计算效率进行测试分析.数值结果表明，针对平行多面体单元，相比于采用高斯数值积分，本文采用的单元分析方法可以大大提高计算效率.在数值实验部分，我们通过具有解析解的模型1和模型2模拟结果与解析解的对比，验证了本文的算法程序是可靠的，满足了正演计算的精度要求，且计算效率较高.进一步地，模型3和模型4对应穿过异常体的垂直井和倾斜井的情况，模拟结果表明，当地下存在异常体时，在地表可以观测到电位异常且电位异常的分布与异常体在地面的投影位置一致，当供电点电源位于异常体的上方和下方，地表电位异常出现相反的极值；且井的倾斜和供电点位置的变化等都会影响到地面电位异常的分布，在实际工作中必须对这些影响因素综合考虑而作出判断.模型5给出了一个起伏地形的情形，电位等值线断面图的畸变与高低阻异常体对应很好，进一步证明了本文算法是可靠的，可以实现起伏地形条件下的模拟.由于边界效应的存在，靠近边界的部分模拟值误差较大，因此如何能有效地减小剖分区域，而又能消除由于人为地截断计算区域对计算精度的影响，需要进一步地予以解决.
| [1] | Wang T, Signorelli J. Finite-difference modeling of electromagnetic tool response for logging while drilling. Geophysics , 2004, 69(1): 152-160. DOI:10.1190/1.1649383 |
| [2] | Avdeev D B, Kuvshinov A V, Pankratov, et al. Three dimensional induction logging problems, part 1:An integral equation solution and model comparisons. Geophysics , 2002, 67(2): 413-426. DOI:10.1190/1.1468601 |
| [3] | Nam M J, Pardo D, Toores-Verdín C. Simulation of DC dual-laterolog measurements in complex formations:A Fourier-series approach with nonorthogonal coordinates and self-adapting finite elements. Geophysics , 2009, 74(1): E31-E43. DOI:10.1190/1.3000681 |
| [4] | 岳建华, 刘志新. 井-地三维电阻率成像技术. 地物理学进展 , 2005, 20(2): 407–411. Yue J H, Liu Z X. Three dimension resistivity tomography of mine-ground. Progress in Geophysics (in Chinese) , 2005, 20(2): 407-411. |
| [5] | 王志刚, 何展翔, 刘昱. 井地直流电法三维数值模拟及异常规律研究. 工程地球物理学报 , 2006, 3(2): 87–92. Wang Z G, He Z X, Liu Y. Research of three-dimensional modeling and anomalous rule on borehole-ground DC method. Chinese Journal of Engineering Geophysics (in Chinese) , 2006, 3(2): 87-92. |
| [6] | 徐凯军, 李桐林. 垂直有限线源三维地电场有限差分正演研究. 吉林大学学报(地球科学版) , 2006, 36(1): 137–147. Xu K J, Li T L. The forward modeling of three-dimensional geoelectric field of vertical finite line source by finite-difference method. Journal of Jilin University (Earth Science Edition) (in Chinese) , 2006, 36(1): 137-147. |
| [7] | 柯敢攀, 黄清华. 井地电法的三维正反演研究. 北京大学学报(自然科学版) , 2009, 45(2): 264–272. Ke G P, Huang Q H. 3D forward and inversion problems of borehole-to-surface electrical method. Acta Scientiarum Naturalium Universitatis Pekinensis (in Chinese) , 2009, 45(2): 264-272. |
| [8] | 阮百尧, 熊彬, 徐世浙. 三维地电断面电阻率测深有限元数值模拟. 地球科学 , 2001, 26(1): 73–77. Ruan B Y, Xiong B, Xu S Z. Finite element method for modeling resistivity sounding 3-D geo-electric section. Earth Science-Journal of China University of Geosciences (in Chinese) , 2001, 26(1): 73-77. |
| [9] | 阮百尧, 熊彬. 电导率连续变化的三维电阻率测深有限元模拟. 地球物理学报 , 2002, 45(1): 131–138. Ruan B Y, Xiong B. A finite element modeling of three-dimensional resistivity sounding with continuous conductivity. Chinese J.Geophys. (in Chinese) , 2002, 45(1): 131-138. |
| [10] | 徐世浙. 地球物理中的有限单元法. 北京: 科学出版社, 1994 : 178 -188. Xu S Z. The Finite Element Method in Geophysics (in Chinese). Beijing: Science Press, 1994 : 178 -188. |
| [11] | 《数学手册》编写组. 数学手册. 北京: 人民教育出版社, 1979 : 990 -998. 《Math Handbook》Compile Group. Math Handbook (in Chinese). Beijing: Education Publishing House, 1979 : 990 -998. |
| [12] | Saad Y.Iterative Methods for Sparse Linear Systems, PWS, Bostion, 1996 http://www.oalib.com/references/16139760 |
| [13] | Barrett R, Berry M, Chan T, et al.Templates for the Solution of Linear Systems:Building Blocks for Iterative Methods.SIAM, Philadelphia, 1994 http://www.oalib.com/references/18991244 |
| [14] | 吴小平, 汪彤彤. 利用共轭梯度算法的电阻率三维有限元正演. 地球物理学报 , 2003, 46(3): 428–432. Wu X P, Wang T T. A 3D finite element resistivity forward modeling using conjugate gradient algorithm. Chinese J.Geophys. (in Chinese) , 2003, 46(3): 428-432. |