起伏地形条件下电阻率法三维反演
韩雪, 朱伟国
广东省地质物探工程勘察院,广东 广州 510800

作者简介: 韩雪(1988-),女,助理工程师,从事地球物理勘查工作。Email:kahanxue@126.com

摘要

地形起伏下的三维反演问题是当前电阻率法研究的一个难题。为了更好地解决上述问题,采用加权正则化共轭梯度法实现起伏地形三维电阻率反演算法。该方法引入了加权正则化思想,显著降低了迭代时目标函数发散的问题,提高了反演稳定性。对比了两种消除反演中地形影响的方法,结果表明,直接带地形的电阻率法三维反演具有更好的分辨率,能有效地消除地形所造成的误差,但在起伏角度偏大,如河流、堤坝等接近垂直角度时,使用此方法会使得反演发散得不到满意的结果;此时采用基于地形校正的方法有一定的效果。

关键词: 电阻率法; 三维反演; 加权正则化共轭梯度法; 起伏地形
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2017)01-0136-05
A discussion on 3D inversion result of the resistivity method under the condition of undulate topography
HAN Xue, ZHU Wei-Guo
Geological and Geophysical Engineering Exploration Institute of Guangdong Guangzhou 510800, China
Abstract

3D resistivity inversion under the condition of undulate topography is a difficult problem in the current resistivity method. In order to solve this problem, the authors use a weighted regularized conjugate gradient method to accomplish three-dimensional resistivity inversion algorithm of undulating terrain. The method introduces weighted regularization thought which can significantly reduce the problem of the objective function iteration divergence and improve the stability of the inversion. A comparison of the two methods in eliminating topographic influence in the inversion shows that the direct use of the terrain resistivity method of 3D inversion can yield better resolution, which can effectively eliminate the error caused by terrain. Nevertheless, under the condition that the rolling angle is almost vertical such as the cases of rivers and dams. The continual utilization of this method can hardly get a satisfactory result because it causes the inverse divergence. In such cases, the use of topographic correction method has some advantageous effects.

Keyword: resistivity method; 3D inversion; weighted regularized conjugate gradient method; undulate topography
0 引言

电阻率法作为电法勘探的一个重要分支, 在普查金属、非金属矿产、水文、工程、环境地质调查以及勘查能源等方面都发挥着重要作用[1, 2, 3, 4]。在进行电法勘探的实际工作中, 时常会遇到各种起伏地形的情况, 如山脊、山谷、陡坎等地形。地形起伏时测得的视电阻率数据包含了有用异常和地形异常, 当地形起伏较大且复杂时, 依据电阻率资料进行解释就可能会得出错误的结果。

对如何消除起伏地形的影响, 国内外很多学者进行了研究, 常用的办法是通过数值模拟用纯地形响应对观测数据进行地形校正, 正演模拟采用边界单元法[5, 6]、有限差分法[7]或有限单元法[8]。地形校正的方法对于削弱地形的影响有一定的效果, 但只有在地下电性均匀时才会完全消除地形的影响, 为了克服这一问题, 国内外学者研究了包含地形的二维或三维反演[9, 10, 11, 12, 13, 14]。由于实际的地下电性结构多为三维地质体, 在野外勘探仍以二维测量方式为主的前提下, 就需要充分利用由高密度电阻率仪采集的多条二维剖面数据来展开起伏地电条件电阻率法三维反演。本文采用加权正则化共轭梯度法来求解最优化问题, 由于引入了加权正则化思想, 显著降低了迭代时目标函数发散的问题, 提高了反演的稳定性。通过对理论模型的分析给出了地形校正和直接带地形反演两种消除地形影响方法的优劣, 在地形不平条件下为电阻率法资料解释提供参考。

1 三维反演算法
1.1 目标函数

针对电阻率法反演的不适定性[15], 引入Tikhonov正则化方法[16, 17, 18, 19, 20], 得到目标函数

fα(m, d)=ϕ(m)+αS(m), (1)

其中, m表示模型参数电阻率, 单位为Ω · m; d为电阻率法数据观测值, 即视电阻率值, 单位为Ω · m; ϕ (m)=A(m)-d2是数据目标函数, 它是预测数据和实测数据之差的平方和, A为非线性正演算子, 文中采用有限单元法实现; α 是正则化因子; S(m)是模型目标函数, 文中采用基于模型与先验信息mapr的最小二范数:S(m)=‖ (m-mapr)‖ 2

为了提高解的稳定性, 引入观测数据权值和模型参数权值, 使可靠的数据在反演中发挥更大的作用, 由此目标函数可表示为

fα(m, d)=WdA(m)-Wdd2+αWm(m-mapr)2, (2)

式中, WdWm分别表示数据和模型参数的权系数矩阵。

1.2 正则化因子的选择

正则化因子α 是在最佳拟合和最合理稳定之间起平衡作用的纽带。当α 很小时, 求得的模型非常粗糙, 会导致反演的不稳定性和不正确性; 当α 很大时, 求得的模型会与先验信息非常接近, 数据得不到很好的拟合。因此, 正则化因子的选择直接关系到反演结果的好坏。

选择一种指数函数:

αk=α1qk-1; k=1, 2, 3, ; 0< q< 1(3)

在初始迭代时取α 0=0, 正则化参数α 1由初始迭代计算所得, 表达式为

α1=WdA(m1)-Wdd2Wm(m1-mapr)2(4)

该方法简单易行, 保证了目标函数的收敛性且不需要任何初始的经验值, 减少了计算量。

1.3 反演算法流程

电阻率反演算法的流程如图1所示, 针对目标函数 fα(m, d)的最小值问题, 利用加权正则化共轭梯度法求解, 使得反演算法收敛、稳定。

图1 加权正则化共轭梯度三维电阻率反演算法流程

2 结果分析
2.1 偶极装置在起伏地形上的异常

研究所用装置类型为偶极装置, 该装置起伏地形的视电阻率异常不同于三极装置, 由于布极方式和记录点规定有所不同使得该装置异常较为复杂。起伏地形用不同角度的二维陡坎地形来模拟, 地下为均匀介质, 电阻率为1 Ω · m。图2是该模型下偶极装置的视电阻率等值线, 起伏角度分别为30° 、45° 、60° 和85° , 以斜坡为中心, 呈不对称“ 八字形” , 左侧为“ 低阻” 异常, 即极小值在斜坡低高程一侧, 右侧为“ 高阻” 异常, 即极大值在斜坡抬升一侧, 受影响的深度随极距变大。α 等于30° 时极小值约为介质电阻率的0.6倍, 极大值约为介质电阻率的1.4倍; α 等于45° 时极小值约为介质电阻率的0.4倍, 极大值约为介质电阻率的1.6倍; α 等于60° 时极小值约为介质电阻率的0.2倍, 极大值约为介质电阻率的1.8倍; α 等于85° 时极小值约为介质电阻率的0.1倍, 极大值约为介质电阻率的3.5倍; 故异常幅度随着起伏角度的增加而增大, 且α 从60° 到85° 之间的变化最为剧烈, α 等于85° 时的极大值约为α 等于60° 时的2倍。

2.2 三维反演算例分析

采用了2种消除电阻率三维反演中地形影响的方法, 一是常用的“ 比较法” , 将观测数据的视电阻率值与单纯地形的数值模拟值相比来进行地形改正然后对改正后的数据进行反演; 二是抛开地形校正这一概念, 直接将地形带入反演算法中即带地形电阻率法三维反演。

模型一:如图3所示为坡度30° 二维陡坎地形模型, 陡坎相对高度为4 m, 均匀半空间中有一个 4 m× 4 m× 4 m的低阻立方体模型, 电阻率为1 Ω · m, 异常体的中心坐标为(0 m, 0 m, -4 m), 围岩电阻率为100 Ω · m, 测量点距为1 m, 测线沿x方向布置线距为2 m, 每条剖面上有25个供电点, 共有5条剖面。

通过带地形的正演数值模拟获得“ 实测数据” , 再由该数据得到背景电阻率, 以此作为反演初始模型电阻率。不考虑地形因素, 直接进行水平地面电阻率三维反演(图4), 可见在剖面的左边为8 m× 6 m大面积的低阻异常, 视电阻率小于30 Ω · m, 在剖面右边为条带状高阻异常, 视电阻率大于110 Ω · m, 这与真实模型相比显而易见产生了很大的偏差, 对解释结果造成错误结论, 因此起伏地形条件下的电阻率三维反演必须消除地形影响。

图2 偶极装置纯地形异常的视电阻率等值线断面(单位:Ω · m)

图3 坡度30° 二维陡坎地形低阻体(模型一)示意

图4 不考虑地形测线y=0 m的垂直断面

采用两种消除地形影响的方法分别进行反演。图5为基于地形改正在y=0 m处的垂直断面和在z=-4 m处的水平断面, 相比于不考虑地形的反演取得了一定的效果, 基本反映出了低阻异常体的存在, 但中心位置有所偏离且视电阻率值相对真实电阻率偏高, 可见地形改正难以完全消除地形影响。 图6为直接将地形带入反演在y=0 m处的垂直断面和在z=-4 m处的水平断面, 三维反演结果能很好地还原低阻体的位置、形状和大小, 与地下真实模型接近, 没有出现多余的构造, 相比于基于地形改正的反演结果, 三维反演显示了更好的分辨率, 得到较好的反演效果。

在模型一的基础上设计坡度为60° 的二维陡坎地形模型, 陡坎相对高度为6.5 m, 异常体的中心坐标为(0 m, 0 m, -2 m), 其他参数同模型一。将地形直接带入反演, 图7为带入地形时该模型在y=0 m处的垂直断面。该图反映了低阻体的位置、形状和大小, 与真实模型接近, 没有出现多余的构造。继续增加起伏坡度, 设计坡度约90° 的二维陡坎地形模型, 陡坎相对高度为5 m, 异常体的中心坐标为(0 m, 0 m, -3 m), 其他参数同模型一。当起伏坡度增加, 如果仍然采用带地形的三维反演, 由于拐点处的值急剧增大或减少使得反演发散, 这时只能用比较法消弱地形影响之后做平地形的反演。图8为该模型y=0 m处的垂直断面, 可以看出异常中心向左上偏, 但该图可以反映出异常体的大致位置和大小。

图5 模型一地形改正反演结果

图6 模型一带地形反演结果

图7 带地形测线y=0 m的反演断面

图8 地形改正测线y=0 m的反演断面

3 结论

利用加权正则化共扼梯度方法实现起伏地形三维电阻率反演成像算法, 模型计算结果表明, 该反演算法较稳定、快速, 反演结果良好; 在进行有地形的电阻率法三维反演时必须要考虑到地形的影响, 否则会导致不正确的解释。文中采用两种方法来消除地形影响, 通过对理论模型的研究分析得出:在地形起伏角度不大于60° 时, 将地形直接带入电阻率法三维反演能取得较好的结果; 在地形起伏角度接近垂直时, 拐点处的异常变得很复杂, 此时不能再将地形直接带入反演算法, 只能依靠地形校正的方法。

The authors have declared that no competing interests exist.

参考文献
[1] 李金铭. 地电场与电法勘探[M]. 北京: 地质出版社, 2005. [本文引用:1]
[2] Ramirez A, Daily W, Binley A, et al. Detection of leaks in underground storage tanks using electrical resistance methods[J]. Journal of Environmental and Engineering Geophysics, 1996, 1(3): 189-203. [本文引用:1]
[3] Oldenburg D W, Li Y G, Ellis R G. Inversion of geophysical data over a copper gold porphyry deposit: A case history for Mt. Milligan[J]. Geophysics, 1997, 62(5): 1419-1431. [本文引用:1]
[4] Kemna A, Kulessa B, Vereecken H. Imaging and characterisation of subsurface solute transport using electrical resistivity tomography(ERT) and equivalent transport models[J]. Journal of Hydrology, 2002, 267(3): 125-146. [本文引用:1]
[5] 徐世浙. 电阻率二维地形改正的地质效果[J]. 地质与勘探, 1989, 25(5): 43-45. [本文引用:1]
[6] 汤洪志, 刘庆成, 龚育龄. 边界单元法在高密度电阻率法二维地形改正中的应用效果[J]. 物探与化探, 2001, 25(6): 457-459. [本文引用:1]
[7] Fox R C, Hohmann G W, Killpack T J, et al. Topographic effects in resistivity and induced-polarization surveys[J]. Geophysics, 1980, 45(1): 75-93. [本文引用:1]
[8] Holcombe H T, Jiracek G R. Three-dimensional terrain corrections in resistivity surveys[J]. Geophysics, 1984, 49(4): 439-452. [本文引用:1]
[9] Tong L T, Yang C H. Incorporation of topography into two-dimensional resistivity inversion[J]. Geophysics, 1990, 55(3): 354-361. [本文引用:1]
[10] Yi M J, Kim J H, Song Y, et al. Three-dimensional imaging of subsurface structures using resistivity data[J]. Geophysical Prospectig, 2001, 49(4): 483-497. [本文引用:1]
[11] 吴小平. 非平坦地形条件下电阻率三维反演[J]. 地球物理学报, 2005, 48(4): 932-936. [本文引用:1]
[12] Günther T, Rücker C, Spitzer K. Three-dimensional modelling and inversion of dc resistivity data incorporating topography — II. Inversion[J]. Geophysical Journal International, 2006, 166(2): 506-517. [本文引用:1]
[13] Papadopoulos N G, Yi M J, Kim J H, et al. Geophysical investigation of tumuli by means of surface 3D Electrical Resistivity Tomography[J]. Journal of Applied Geophysics, 2010, 70(3): 192-205. [本文引用:1]
[14] Plattner A, Maurer H R, Vorloeper J, et al. 3-D electrical resistivity tomography using adaptive wavelet parameter grids[J]. Geophysical Journal International, 2012, 189(1): 317-330. [本文引用:1]
[15] Tikhonov A N, Leonov A S, Yagola A G. Non-linear Illposed Problems[M]. London: Chapman and Hall, 1998. [本文引用:1]
[16] Ellis R G, Oldenburg D W. The pole-pole 3-D DC-resistivity inverse problem: a conjugate gradient approach[J]. Geophysical Journal International, 1994, 119(1): 187-194. [本文引用:1]
[17] Li Y, Oldenburg D W. Inversion of 3-D DC resistivity data using an approximate inverse mapping[J]. Geophysical Journal International, 2007, 116(3): 527-537. [本文引用:1]
[18] Pain C C, Herwanger J V, Worthington M H, et al. Effective multidimensional resistivity inversion using finite-element techniques[J]. Geophysical Journal International, 2002, 151(3): 710-728. [本文引用:1]
[19] Pidlisecky A, Haber E, Knight R. RESINVM3D: A 3D resistivity inversion package[J]. Geophysics, 2007, 72(2): H1-H10. [本文引用:1]
[20] Marescot L, Lopes S P, Rigobert S, et al. Nonlinear inversion of geoelectric data acquired across 3D objects using a finite-element approach[J]. Geophysics, 2008, 73(3): F121-F133. [本文引用:1]