2017, Vol. 32
2. 三峡大学土木与建筑学院, 宜昌 443002
2. China Three Gorges University College of civil engineering and architecture, Yichang 443002, China
电阻率法是以地壳中岩石和矿石的导电性差异为物理基础，通过观测与研究人工建立的地下介质中电流场 (稳定场或交变场) 的分布规律进行找矿和解决地质问题的一种电法勘探分支方法 (付良魁和李金铭，1980).该方法具有确定地下异常能力强，应用领域广的优点, 被广泛应用在环境、工程地质 (吴小平和徐果明，2000；Dahlin，2001)、水文 (Wilkinson et al., 2006)、考古 (Griffiths et al., 1990)、矿产资源勘察 (Bauman，2005) 等领域.电阻率层析成像是20世纪80年代末发展起来的一种新的地球物理探测方法 (白登海和于晟，1995)，成像范围从微米级至千米级 (Storz et al., 2000).除了在地表勘探之外该方法还可以应用到跨孔电阻率 (雷旭友等，2009；潘纪顺等，2010) 以及水上勘探中 (Loke and Lane Jr，2004).
随着仪器硬件以及正反演技术研究的发展，电阻率法勘探从传统的两个电极供电两个电极采集发展到现在的两个电极供电多道采集技术.电阻率反演技术相应地从1D、2D (王兴泰等，1995；阮百尧和徐世浙，1996)、2.5D (毛先进和鲍光淑，1999) 发展到3D (吴小平等，2015).采集技术和反演技术的发展使电阻率层析成像有可能解决更为复杂的地质构造问题.尽管现阶段3D电阻率成像应用还达不到2D应用水平，但在一些复杂的环境及工程问题上情况较为复杂，需要提高勘探准确度，对3D电阻率成像具有迫切的需求.不同的电阻率勘探观测系统，对地下介质的分辨率也不同.常规的电阻率成像通常采用2个电极供电2个电极采集的方式，并通过滚动4极的位置以组合成不同观测方式进行测量.近年来，多通道并行采集得到快速的发展.其中包括超高密度 (雷旭友等，2009；潘纪顺等，2010) 实现地面、井地、井井观测使电阻率法得到更广泛的应用，但其无法准确确定异常体沿垂直剖面方向的延伸长度；三维E-scan (黄俊革等，2006) 采集方式提高了采集电位的数据量，获得更好的观测角度，其对解决电阻率纵向分辨率具有一定的局限性；井地密集阵列数据采集 (农观海等，2013) 分别采集纵向及横向两个相邻测线电极之间电位差，比常规四极采集电阻率剖面更多，可提供更丰富的地电信息.
在充分利用房屋、桥梁、隧道、水利大坝、煤矿采空区等地基地质勘探阶段已有工程勘探孔基础上，本文提出井地井三维全通道电阻率成像方法.该方法同时具有超高密度电法、E-scan采集方式、井地联合等方式的优点，又充分考虑场地较小等因素提高电阻率法勘探的纵横向分辨率及成像的可靠性.
全通道采集为除去2个供电电极外其余全部电极均采集电位场数据并参与成像计算，增加了观测电位数据的采集量，例如采用32个电极采集，任选两个电极供电则选择的次数为C322，每次供电采集的电位数为30个，所以全通道采集得到的电位值数为C322×30=29760个.温纳装置采集时，32个电极最多可以采集10层，所以其数据采集个数为
本项研究从全通道采集出发，研究两个供电电极在不同位置供电时电场的分布规律以及3D静电场的差分求解算法，进而求出理论模型下的观测电位值.对于电阻率反演，最为关键问题为反演系数矩阵，因显式写出为稠密矩阵难以存储，在此我们采用伴随矩阵虚拟点源的方式求取，迭代反演采用共轭梯度方法.
1 方法理论直流电阻率层析成像主要包括正演和反演两个部分.要解决直流反演成像问题，必须首先实现电阻率正演，即求解稳定电流场.其次反演问题主要涉及到目标函数选取、搜索方向及步长、迭代方法及终止条件.我们下面主要针对电阻率正反演问题进行讨论.
1.1 3D电阻率正演因地下介质电阻率沿空间三个方向均发生变化，3D空间稳定电流场方程可表示为
|
(1) |
式中：x0、y0、z0为供电电极位置；ρ为电阻率；I为供电电流；u为电位.
对式 (1) 进行离散并用矩阵向量的形式表达，即：
|
(2) |
其中：A为五对角带状系数矩阵；u为待求解网格节点处电位值；f为-Iδ(x-x0)δ(y-y0)δ(z-z0).在自由地表处，采用镜像边界条件；在其他外边界，采用扩大模型边界区域方法，令电位值近似为零.
通过正演计算可以求取空间某一网格点处电位对地下模型网格点上所有电阻率的偏导数
|
(3) |
设uρ=∂u/∂ρ为观测电位对地下介质模型电阻率偏导数，进一步整理可得：
|
(4) |
式 (4) 和式 (1) 有着相同的方程形式，对式 (4) 进行离散化，并写成矩阵可得：
|
(5) |
其中矩阵B为稀疏的五对角矩阵，具有与矩阵A类似的结构.通过式 (5) 求逆可得：
|
(6) |
电位对电导率的偏导数计算步骤是：首先通过式 (2) 正演求得电位场向量u；然后将系数矩阵B和向量u相乘得到ψ=Bu，ψ为虚拟场源向量；最后再通过一次正演基于式 (6) 计算得到地下某点处电位对电阻率的偏导数μρ.
1.2 3D电阻率反演在利用实际采集到的电位数据值做电阻率层析成像时，需要对反演模型进行离散化.通常根据观测电极距的大小来确定模型剖分网格间距的大小，一般选择网格间距为电极距的一半.但实际的供电电极位置及观测数据位置不一定正好落在离散网格节点上，因此需要将采集到的不规则数据及供电电极位置投影到规则的离散网格上.投影可以通过双线性插值并以投影矩阵Q的方式来实现 (Pidlisecky et al., 2007).公式为
|
(7) |
其中：d为观测数据；Q为投影矩阵；u为供电在所有网格上的电位值，每一列代表一次供电模型所有网格的电位值，每一行为同一模型网格在不同次供电产生的电位值 (矩阵u列数为供电电极对数、行数为模型总网格数)；f为震源函数.
因每一次供电采集位置互不重叠，所以投影矩阵Q为正交矩阵Q-1=QT，则投影后的规则网格上分布的电位值数据为
|
(8) |
考虑到模型电阻率不能为负值以及电阻率的数量级可能变化很大，为了反演稳定，这里采取将电阻率模型参数取电导率对数的方式可以保证电阻率不为负值，如下式所示为
|
(9) |
此时，在规则网格上的观测电位u对模型m的求导为
|
(10) |
其中e-m为一主对角矩阵.
定义雅克比矩阵J为观测电位数据d对模型m的求导，有：
|
(11) |
电阻率成像目标函数ϕ(m) 分为两个部分，第一部分表示理论数据和观测数据的拟合程度，第二部分表示反演模型与参考模型的差距.通过目标函数的最小化，得到的电阻率模型既能较好的拟合观测数据又能靠近参考模型，具体表示如下：
|
(12) |
其中：dobs观测数据值；d(m) 模型正演数据值；m为迭代模型；mref为参考模型；Wd为观测数据加权矩阵；Wm为模型约束加权矩阵；β为观测数据拟合和模型约束之间的权重.
对式 (12) 所表示的目标函数，求极值，并线性化得到：
|
(13) |
其中：
δm为模型修改量；设g=JTWdTWd(d-dobs)-βWmTWm(m-mref) 为目标函数的梯度；H=JTJ+βWmTWm为近似汉森矩阵.
基于式 (13) 的最优化采用高斯-牛顿算法，初始电阻率模型一般给定一个均匀模型，采取观测视电阻率的平均值作为初始模型的电阻率值.写成线性方程组形式，有：
|
(14) |
在反演迭代时，通过采用共轭梯度迭代方法，避免显示求出雅克比、汉森矩阵，减少计算及存储难度，具体实现时通过使模型梯度与雅克比矩阵相乘的方式，降低雅克比矩阵的存储难度及避免显式计算雅克比矩阵.
当通过线性方程组求解器求出δm后，可以用式 (14) 对初始模型进行迭代，公式为
|
(15) |
其中α为模型迭代步长，可以采用线性搜索方法决定它的大小 (Armijo，1966)，相应搜索目标函数为
|
(16) |
其中
c1为一常数，通常取10-4.至此，通过共轭梯度迭代，对模型进行了更新迭代以满足反演目标函数，对模型不停迭代修正，直至满足终止条件.
2 全通道装置与温纳装置反演分辨率比较全通道采集可以实现常规高密度地表采集，同时可以实现任意观测装置方式采集，具有更好的采集方位与数据采集量，对地层电性结构具有更好的约束作用.在高密度电阻率法勘探多种装置形式中，温纳装置具有最小的装置系数，即具有最强的抗干扰能力及纵向分辨率等优点，为工程高密度勘探中采用最多的装置形式.所以本文以二维温纳装置为例说明全通道装置成像优势.图 1a所示模型背景电阻率为100 Ω·m，两个低电阻异常电阻率值为10 Ω·m.图 1b为温纳装置视电阻率剖面图.图 1c和图 1d分别为温纳反演结果和全通道装置反演结果图.通过反演结果可以看出：温纳装置与全通道装置都可以较好的分辨出异常所在的位置.但全通道电阻率成像能够更好地恢复背景电阻率值，得到的背景电阻率更平滑均匀，且低阻异常中心位置与形态收敛的更为准确.
|
图 1 全通道和温纳装置电阻率成像比较 (a) 电阻率模型 (在100 Ω·m的背景电阻率上存在两个10 Ω·m的低电阻异常); (b) 温纳装置视电阻率剖面; (c) 温纳装置反演结果; (d) 全通道装置反演结果. Figure 1 Comparison of resistivity tomography by Wenner and full channel survey (a) Resistivity model consisting of 2 low resistivity anomalies of 10 Ω·m embedded on the background resistivity of 100 Ω·m; (b) Apparent resistivity profile with Wenner survey; (c) Inverted model by the Wenner survey; (d) Inveretd resistivity model by full channel survey. |
下面从井地井装置形式、双电极供电电场分布及对典型异常合成数据反演测试说明井地井联合观测的成像效果.我们分别建立了单一低阻异常、单一高阻异常、高低阻双异常模型.
3.1 井-地-井全通道装置设计在工程地质勘探有地质勘探钻孔的情况下，分别将电极布置在地表与孔中 (图 2).这里我们根据实际工程情况设计了4个井的井-地-井三维观测方式，整个观测装置共包含103个电极，每个井中分别放置10个电极，地表放置63个电极.电极距1 m，井深12 m，井距在6~10 m之间变化, 图中矩形点为地表电极位置，圆形点为井中电极位置.
|
图 2 井-地-井电阻率观测装置 (a) 平面俯视图; (b) 立体图. Figure 2 Well-surface-well resistivity survey configuration (a) Map view; (b) 3D view. |
图 3为一次正负双电极供电时通过 (2) 式计算出的三维空间产生电场分布，其中曲面为电位等势面.通过计算出空间任意位置电位，然后根据各电极空间位置利用 (7) 式求取理论观测电位值.其中理论模型背景电阻率为100 Ω·m的均匀模型，通过比较数值计算的电位场与解析求出的电位场数值，数值计算的电位场值误差为0.01%水平，说明本文采用的有限差分的数值计算方法正确有效，满足求解精度要求.
|
图 3 正负双电极供电三维电场分布 Figure 3 Three dimensional electric field distribution for bi-electrodes power supply |
图 4为单一低阻异常模型反演结果，异常大小为2×2×2 m，距离地面4 m，异常电阻率为10 Ω·m，背景电阻率为100 Ω·m.图 5为单一高阻异常模型反演结果，异常大小与低阻模型大小相同，距离地面4 m，异常电阻率为1000 Ω·m背景电阻率为100 Ω·m.从反演结果可以看出，单一低阻和高阻异常模型的反演结果均能较好的反应出空间分布位置及形态.
|
图 4 单个低阻异常反演结果 (a) 电阻率模型; (b) 反演结果. Figure 4 3D resistivity inversion with a single low resistivity anomaly (a) Resistivity model; (b) Inversion result. |
|
图 5 单个高阻异常反演结果 (a) 电阻率模型; (b) 反演结果. Figure 5 3D resistivity inversion with a single high resistivity anomaly (a) Resistivity model; (b) Inversion result. |
图 6为高低阻双异常组合模型反演结果，异常大小均为2×2×2 m，距离地面4 m，异常距离4 m，背景电阻率为100 Ω·m，高阻异常1000 Ω·m，低阻异常10 Ω·m.从反演结果可以看出，高低阻异常中心位置、空间形态均成像较好，且没有产生多余假异常.
|
图 6 高低阻异常模型反演结果 (a) 电阻率模型; (b) 反演结果. Figure 6 3D resistivity inversion with high and low resistivity anomalies (a) Resistivity model; (b) Inversion result. |
以单一低阻异常模型三维井地井联合观测反演结果说明拟合残差及对观测数据的拟合效果.图 7为迭代100次残差衰减曲线图，从图中可以看出，随着迭代次数的增加，残差呈平稳下降，说明本文所采用的共轭梯度反演算法稳定.迭代至40次后残差变化量较小说明迭代到稳定的阶段，误差百分比从开始的100%下降到总体误差的15%左右.一般迭代终止条件包括，迭代次数、每次迭代残差下降的梯度、残差剩余量占初始残差百分比、模型更新量等.本文迭代终止及收敛条件为电阻率或地震走时残差下降量低于某一百分比时迭代终止 (如迭代残差下降量小于1%).当选择迭代残差下降1%作为判断准则时，迭代次数为44次.
|
图 7 残差衰减曲线 (迭代100次) Figure 7 Attenuation curve of residual error (iteration 100 times) |
本次三维正反演模型网格数为10×12×14共1680个网格，一次供电可计算出1680个电位值.因为本文装置电极总数为103个，除去每次两个正负供电电极，每次供电采集到101个电位值.为验证本文反演算法对观测数据的拟合程度，我们比较了反演得到的电阻率模型结果一次供电正演的电位值与理论模型对应电极供电正演得到的电位值.
图 8为理论模型一次供电的观测数据和图 7中迭代44次反演结果一次供电观测值对比图，从图中可以看出反演结果计算出的值与理论观测值基本一致，理论模型的正演电位值与反演结果正演电位值数量级均在1e2，两者之差的数量级为1e-1，即根据反演结果计算的观测电位值与根据理论模型计算得到的观测电位值误差在0.1%级别，说明反演结果与理论模型的吻合程度较高.
|
图 8 理论模型与反演得到模型正演计算的观测电位值对比 (a) 理论模型与反演得到的模型的正演得到的观测电位数值; (b) 反演得到的模型正演计算的观测电位值减去理论模型正演得到电位值之差曲线. Figure 8 Contrast electrical potential value calculated from theoretical model and inversion modeling |
全通道装置与温纳装置相比，全通道电阻率成像方法能够较好地恢复背景电阻率值，且计算得到的背景电阻率更平滑均匀，对低阻异常中心位置收敛与形态恢复更为准确.
4.2从理论模型的全通道电阻率成像反演结果可以看出：单一低阻和单一高阻及高低阻组合的异常模型的反演结果均能较好的反应出异常的空间分布形态及异常的中心位置.
4.3在采用的共轭梯度算法全通道电阻率成像反演中，随着迭代次数的增加，残差呈平稳下降，且最后达到稳定的阶段.从而说明共轭梯度算法在全通道电阻率成像反演中的收敛性及稳定性.
4.4本文所采用的三维井间及地表结合的直流电阻率装置形式，为解决工程勘探中场地较小且需要勘探一定深度需求的地质问题提供了一种有效的勘探方法，且能充分利用已有的勘探钻孔，具有较好的应用前景.
4.5通过对理论模型与全通道电阻率成像反演得到的模型所计算的电位值对比分析可以得出，本文所采用的反演方法能够取得较好的反演结果.
4.6全通道采集数据量巨大，会一定程度增加反演计算时间；同时全通道数据采集时在供电电极位置附近会产生电位变化较大，在反演计算时可适当增大电极位置处网格模型的平滑权重，以达到较为稳定的反演结果.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Bai D H, Yu S. 1995. Theory and methods of resistivity tomography[J]. Progress in Geophysics (in Chinese), 10(1): 56–75. |
| [] | Bauman P. 2005. 2-D resistivity surveying for hydrocarbons-a primer[J]. Focus, 30(4). |
| [] | Dahlin T. 2001. The development of DC resistivity imaging techniques[J]. Computers & Geosciences, 27(9): 1019–1029. |
| [] | Fu L K, Li J M. 1980. Electrical Prospecting Tutorial (in Chinese)[M]. Beijing: Geological Publishing House. |
| [] | Griffiths D H, Turnbull J, Olayinka A I. 1990. Two-dimensional resistivity mapping with a computer-controlled array[J]. First Break, 8(4): 121–129. |
| [] | Huang J G, Wang J L, Ruan B Y. 2006. A study on FEM modeling of anomalies of 3-D high-density E-SCAN resistivity survey[J]. Chinese Journal of Geophysics (in Chinese), 49(4): 1206–1214. DOI:10.3321/j.issn:0001-5733.2006.04.037 |
| [] | Lei X Y, Li Z W, Zhe J P. 2009. Applications and research of the high resolution resistivity method in explovation of caves, mined regions and Karst region[J]. Progress in Geophysics (in Chinese), 24(1): 340–347. |
| [] | Loke M H, Lane Jr J W. 2004. Inversion of data from electrical resistivity imaging surveys in water-covered areas[J]. Exploration Geophysics, 35(4): 266–271. DOI:10.1071/EG04266 |
| [] | Mao X J, Bao G S. 1999. A new technique for 2.5 dimensional resistivity imagery[J]. Geophysical and Geochemical Exploration (in Chinese), 23(2): 150–152. |
| [] | Nong G H, Huang J G, Gao W L, et al. 2013. Theoretical model in the observation method of borehole-ground dense array 3D resistivity/polarizability[J]. Chinese Journal of Engineering Geophysics (in Chinese), 10(1): 29–34. |
| [] | Pang J S, Ge W Z, Zhe J P. 2010. High-density electrical resistivity imaging technique of ground/borehole-to-surface/inter-well[J]. Journal of North China Institute of Water Conservancy and Hydroelectric Power (in Chinese), 31(2): 74–78. |
| [] | Pidlisecky A, Haber E, Knight R. 2007. RESINVM3D:A 3D resistivity inversion package[J]. Geophysics, 72(2): H1–H10. DOI:10.1190/1.2402499 |
| [] | Ruan B Y, Xu S Z. 1996. Rapid inversion of two-dimensional D[J]. C. resistivity sounding curve[J]. Geophysical and Geochemical Exploration (in Chinese), 20(6): 455–460. |
| [] | Storz H, Storz W, Jacobs F. 2000. Electrical resistivity tomography to investigate geological structures of the earth's upper crust[J]. Geophysical Prospecting, 48(3): 455–471. DOI:10.1046/j.1365-2478.2000.00196.x |
| [] | Wang X T, Wang J S, Li X Q. 1995. A new method for the reconstruction of two-dimensional resistivity image[J]. Geophysical and Geochemical Exploration (in Chinese), 19(1): 54–59. |
| [] | Wilkinson P B, Meldrum P I, Chambers J E, et al. 2006. Improved strategies for the automatic selection of optimized sets of electrical resistivity tomography measurement configurations[J]. Geophysical Journal International, 167(3): 1119–1126. DOI:10.1111/gji.2006.167.issue-3 |
| [] | Wu X P, Liu Y, Wang W. 2015. 3D resistivity inversion incorporating topography based on unstructured meshes[J]. Chinese Journal of Geophysics (in Chinese), 58(8): 2706–2717. DOI:10.6038/cjg20150808 |
| [] | Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics (in Chinese), 43(3): 420–427. DOI:10.3321/j.issn:0001-5733.2000.03.016 |
| [] | 白登海, 于晟. 1995. 电阻率层析成象理论和方法[J]. 地球物理学进展, 10(1): 56–75. |
| [] | 付良魁, 李金铭. 1980. 电法勘探教程[M]. 北京: 地质出版社. |
| [] | 黄俊革, 王家林, 阮百尧. 2006. 三维高密度电阻率E-SCAN法有限元模拟异常特征研究[J]. 地球物理学报, 49(4): 1206–1214. DOI:10.3321/j.issn:0001-5733.2006.04.037 |
| [] | 雷旭友, 李正文, 折京平. 2009. 超高密度电阻率法在土洞、煤窑采空区和岩溶勘探中应用研究[J]. 地球物理学进展, 24(1): 340–347. |
| [] | 毛先进, 鲍光淑. 1999. 2.5维电阻率成像的新方法[J]. .5维电阻率成像的新方法, 23(2): 150–152. |
| [] | 农观海, 黄俊革, 高文立, 等. 2013. 井地密集阵列三维电阻率/极化率观测方法的理论模型研究[J]. 工程地球物理学报, 10(1): 29–34. |
| [] | 潘纪顺, 葛为中, 折京平. 2010. 地面/井地/井间超高密度电阻率成像技术[J]. 华北水利水电学院学报, 31(2): 74–78. |
| [] | 阮百尧, 徐世浙. 1996. 二维直流电阻率测深曲线的快速反演[J]. 物探与化探, 20(6): 455–460. |
| [] | 王兴泰, 王劲松, 李晓芹. 1995. 二维电阻率图像重建的一种新方法[J]. 物探与化探, 19(1): 54–59. |
| [] | 吴小平, 刘洋, 王威. 2015. 基于非结构网格的电阻率三维带地形反演[J]. 地球物理学报, 58(8): 2706–2717. DOI:10.6038/cjg20150808 |
| [] | 吴小平, 徐果明. 2000. 利用共轭梯度法的电阻率三维反演研究[J]. 地球物理学报, 43(3): 420–427. DOI:10.3321/j.issn:0001-5733.2000.03.016 |