基于TV井地2．5D直流电阻率正则化反演

蓝泽鸾，张志勇，邓居智，周 峰，李 曼

（东华理工大学放射性地质与勘探技术国防重点学科实验室，南昌 330013）

基于TV井地2．5D直流电阻率正则化反演

蓝泽鸾，张志勇，邓居智，周 峰，李 曼

（东华理工大学放射性地质与勘探技术国防重点学科实验室，南昌 330013）

针对目前地面直流电阻率法所面临的勘探深度有限的问题，开展了基于TV井地2．5D直流电阻率正则化反演研究。首先从井地直流电满足的边值问题出发，采用三角单元二次插值有限元法推导了2．5D井地电阻率正演公式。为了提高计算效率，对边界采用近似处理、以及利用图论理论的矩阵重排与填入元分析方法，实现了稀疏矩阵的直接解法，大大提高了计算效率。为了能够提高反演结果的稳定以及对异常边界分辨能力，采用L曲线法求取正则化因子、TV的稳定因子以及共轭梯度方法进行反演目标函数的求解，研究表明，反演结果高效、准确。最后对不同井地观测方式的反演结果进行对比分析，得出井地直流电阻率异常特征和分布规律。

井地直流电法；共轭梯度法；有限单元法；TV稳定因子；L曲线

0 引言

直流电阻率法已被广泛用于矿产资源勘探、环境与工程物探中，已成为物探常用方法。传统电阻率测量主要采用的是固定装置（pole-pole、dipole-dipole、Wenner以及Schlumberger）进行地表观测，从而得到电阻率曲线或拟断面图，其大致能显示地下电阻率在垂直方向或水平方向的变化。但是在某些特殊情况下，由于测区的限制或电阻率在水平方向差异很大，导致传统的电法勘探受限而无法勘探到深部的目标体，因此开展井地电阻率成像技术或反演已变成目前的研究重点［1－6］。

井地直流电阻率反演研究主要是针对井－地－井、井－地、井－井等观测方式，Zhou［7］等利用解析灵敏度函数实现了快速二维／三维井间电阻率成像；Wikinson［8］等重点讨论了矿井直流电阻率成像技术问题；国内有关井地直流电阻率反演研究同样取得了相当大的成果，王志刚等［9－11］采用了Born近似实现了井地电法数据的三维反演、井地电法的准解析近似三维反演研究，以及利用积分方程法实现了井地电法的三维模拟的研究；吕玉增等［12］实现了直流电井间三维直接成像；沈平等［13］研究了井间视电阻率的几何成像方法；岳建华等［14］开展了井－地三维电阻率成像技术的研究；徐凯军等［15］实现了基于共轭梯度法的垂直有限线源的三维电阻率反演；柯敢攀等［16］开展了井地电法的三维正反演研究；安然等［17］进行了井地电阻率三维反演研究。对目前来讲，地下情况复杂，现有的电阻率成像技术在一定程度难以满足现有的需求，这里采用了全变差正则化方法，该方法具有一定的边界识别能力，是Rudin等［18］在1992提出，主要应用于图像边界的识别问题；后来Wang等［19］利用混合正则化方法研究了地震波的反演；Acar等［20］对该方法进行了改进，分析该方法的特性。因此借鉴这些理论，采用改进的全变差稳定因子应用到井地直流电阻率反演中，为了提高反演结果的稳定性，运用基于L曲线方法求取正则化因子方法，数值计算表明，反演结果准确、可靠。

1 直流电阻率正演原理

2．5D直流电阻率所满足的数学模型归结为点源二维问题，其经过傅里叶变换后，将变成波数域二维问题［21］：

式（1）中：区域Ω为二维研究区域；Γs、Γ∞是二维区域的边界，Γs为区域Ω的地表边界，Γ∞为区域Ω的地下边界；σ为电导率；r为点电源到边界的距离；n为无穷远边界的外法向方向；k为波数；K0、K1分别为零阶、一阶第二类贝塞尔函数。

将公式（1）转化为变分问题：

对式（2）求变分并令其为零，得到线性方程组（3）。

通过求解线性方程组（3），得到波数域的电位U，然后通过傅氏反变换得到电位。

式中：r为点的位置；km是波数；gm是加权系数［21］（表1）。

表1 傅氏变换系数km和gmTab．1 Fourier transform coefficients，kmand gm

鉴于各种地表的地球物理观测方法，观测数据对浅部分辨率高、向深部逐渐降低，基于这样的认识，设计了基于二叉树的收缩网格（图1）剖分算法。该算法采用三角单元剖分，具有较好地拟合地表地形与地下地质体的能力。为了进一步提高正演计算的效率和精度，作者采用了基于图论理论的矩阵重排与填入元分析算法，实现了基于TV井地2．5D直流电阻率正则化反演算法。该算法采用了改进的Cholesky［22］分解算法，其只需要对总体系数矩阵进行一次分解，然后对于不同的右端向量进行回代即可。正演计算边界对测点的影响大，通常测点距离边界越近，边界的影响就越大，为了减少边界对测量的影响，采用了边界近似，主要是将区域分成目标区和边界区，目标区为异常体存在的区域，其网格步长采用等距剖分；边界区的网格步长采用指数递增的形式来模拟无穷边界。

图1 三角单元剖分示意图Fig．1 Sketch map of triangular element subdivision

2 基于TV井地直流电阻率正则化反演

反演是资料解释的重要环节，为了能够很好地对异常体边界进行识别，研究了基于TV井地直流电阻率正则化反演技术。然而目标函数的最优化求解与正则化因子的求取是正则化反演的关键，反演采用L曲线法，实现正则化因子的求取；利用CG方法求解反演目标函数，提高反演过程的稳定性与计算效率。

2．1 基于TV正则化反演

Tikhonov［23］正则化为：

其中：P（α，m）为目标函数，α是正则化因子；φd（m）也称数据目标函数，通常采用实测数据与预测数据的距离表示：

其中：mapr为先验模型；＾Wm为模型权重矩阵。如果以式（7）为标准形式，则通过引入矩阵＾We可以得到不同稳定因子的统一表示形式［33］。

为了能够有效地识别边界，Rudin提出了Total varivation（TV）因子

其中β为一个小的数。分析知，式（9）作为稳定因子将突出模型梯度变化，能够突出边界识别。

为了利用式（9）统一表示稳定因子，由式（9）构造模型权重函数为式（10）。

其中：m（r）为模型参数分布函数；▽m（r）为模型参数梯度；e为与计算机数值精度有关的一个很小正数。然后将公式（6）、式（7）、式（10）带入到目标函数式（5）中，即可得到基于TV稳定因了的正则化反演目标函数。

2．2 共轭梯度法（CG）

采用共轭梯度法［15，22］求解公式（12）最优化问题，其基本思想跟最速下降法一样：

通过线性搜索使目标函数最小化的方式来求解下降步长：

然后通过对目标函数进行一系列的化简，得到下降步长式（17）。

2．3 L曲线求取正则化因子

对于公式（5），正则化因子α是数据目标函数与模型目标的权重系数，它决定了反演结果的好坏。α的过大与过小都会影响反演效果，如何正确地选择正则化因子是反演的关键。

作者采用由Hansen等［24］提出的L曲线法求取正则化因子，该方法是一种基于数据误差水平未知的启发式选取正则化因子的方法，它是通过对φd（m）与φm（m）在双对数下进行刻画，其曲线形状为“L”，其正则化因子为曲线曲率最大值，一般是“隅角”对应的位置，如图2所示。

图2 L曲线Fig．2 L－curve

曲线曲率的表达式为：

其中：ρ＝φd（m）；η＝φm（m）；η′为η对α的一阶导数。

3 算例分析

3．1 均匀半空间模型

为了验证算法的正确性，设计均匀半空间电阻率为100Ω·m，计算了理论电位曲线并与数值模拟结果进行对比，计算结果如图3所示，相对误差在2%以内。

图3 数值解与解析解电位对比曲线图Fig．3 Potential contrast curve of numerical solution and analytical solution

3．2 模型I

设置如图4所示的地电模型，即均匀半空间存在一个倾斜低阻异常体，异常体的电阻率大小为50Ω·m，背景电阻率的大小为100Ω·m。设计了井－地－井的观测方式，井所在的位置如图5虚线所示，并且反演数据由正演数据加入1%、10%、20%以及30%随机噪声得到，对比不同随机噪声的反演结果，如图5所示。

图5为不同噪声水平的反演结果，随着噪声逐渐增强，反演结果不精细，反演得到的电阻率断面逐渐变得模糊，很难分辨出异常体的倾斜方向以及大致范围。由图5（a）可以看出，低阻异常体清晰可见，其倾斜的方向也能够反应出异常体倾向，异常范围集中，且背景电阻率与真实电阻率相近。随着噪声水平的增加，异常范围相应地反映出异常体所在的位置，但是异常体的倾向难以判别，且反演断面图存在一些假异常，背景电阻率与真实电阻率值相差较大，井所在的位置反演结果随着噪声水平的增加，影响较大。

3．3 组合模型

在均匀半空间存在两个异常体（图6），M1为高阻异常体，电阻率值为200Ω·m；M2为低阻异常体，电阻率值为50Ω·m，背景电阻率值为100 Ω·m。设计两种井地的观测方式，分别为井－地－井与地－井－地，井所在的位置如图7虚线所示，并且反演数据由正演数据加入1%、10%、20%以及30%随机噪声得到，通过对比两种观测方式不同随机噪声的反演结果，其反演结果如图7所示。

图7为不同噪声水平的反演结果，其反演结果都能大致反映出异常所在的位置，但是随着噪声增加，反演结果不精细，且井所在的位置的反演结果误差较大。从图7可以看出，当随机噪声为10%、20%和30%时，反演断面图上存在小区域的假异常，围岩电阻率与真实电阻率误差较大；随机噪声为1%时，背景电阻率与真实电阻率相近，异常体的反演相对集中，其值与真实值相近。

图8为不同噪声水平的反演结果，反演结果都能大致反应出异常所在的位置，但是随着噪声增加，反演结果不精细，且井所在的位置的反演结果误差较大。从反演断面图可以看出，当随机噪声为10%、20%和30%时，反演断面图上存在小区域的假异常，围岩电阻率与真实电阻率误差较大；随机噪声为1%时，背景电阻率与真实电阻率相近，异常体的反演相对集中，其值与真实值相近。

对比图7与图8反演结果发现，井－地－井测量装置相对单井反演效果更明显，异常范围更集中，能够很好地反映出异常的边界，易于圈定异常体的范围。随着噪声水平的增加，井－地－井测量方式相对单井反演结果影响较小，易于圈定异常的范围。

图4 地电模型断面图Fig．4 Sketch of the model

4 结论

作者在前人的基础上，开发了基于TV井地2．5D直流电阻率正则化反演算法，并进行模型试算，试算表明，程序计算正确，算法稳定高效。研究工作得到以下成果：

1）井－地－井测量方式相对单井反演效果更明显，异常范围更集中，能够很好地反映出异常的边界，易于圈定异常体的范围。

2）基于L曲线求取正则化因子的方法以及TV稳定因子的选择，提高了反演结果的稳定性。

3）近似的边界处理以及基于图论理论的矩阵重排，提高了正演计算效率及精度。

图5 不同噪声水平的井－地－井电阻率反演断面图Fig．5 Different noise level of the well－to－well resistivity inversion section

图6 地电模型断面图Fig．6 Sketch of the model

图7 不同噪声水平的地－井－地反演断面图Fig．7 Different noise level of the Surface－Borehole－Surface resistivity inversion section

图8 不同噪声水平的井－地－井反演断面图Fig．8 Different noise level of the well－to－well resistivity inversion section

［1］SMTH，N．C，VOZOFF，K．Two-dimensional dc resistivity inversing for dipole－dipole data［J］．IEEETrans．Geosci．Remote Sensing，1984，GE－22 （02）：21－28．

［2］ZHOU B，GREENHALGH，S A．A synthetic study on cross-hole resistivity imaging with different electrode arrsys［J］．Geophys，1997，28（2）：1－5．

［3］ZHANG，J，MACKIE，R，MADDEN，T．3－D resistivity forward modeling and inversion using conjugate gra-dients［J］．Geophysics，1995，60（6）：1313－1325．

［4］OLDENBURG，D W，MCGILLICRAY，P R，ELLIS，R G．Generalised subspace methods for largescale inverse problem［J］．Geophys．1993，114（4）：12 －20．

［5］LI Y G，OLDENBURG，D W．Approximate inverse mapping in dc resistivity problems［J］．Geophys．1992，109（5）：343－362．

［6］MAURIELLOP，PATELLA，D．Resistivity anomaly imaging by probability tomography［J］．Geophys．1999，47（2）：411－429．

［7］BING ZHOU STEWART A，Greenhalgh．Rapid 2-D／3-D crosshole resistivity imaging using the analytic sensitivity function［J］．GEOPHYSICS，2002，67（3）：755－765

［8］WILKINSON P B，CHAMBERS J E，et al．Extreme sensitivity of cross-hole electrical resistivity tomography measurements to geometric errors［J］．Geophys．2008，173（3）：49－62．

［9］王志刚，何展翔，魏文博．Born近似快速三维反演井地电法数据［J］．球物理学进展，2007，22（02）：508－513．

WANG ZH G，HE ZH X，WEI W B．Fast 3Dinversion borehole ground electrical method data based on born on approximation［J］．progress in geophsycis，2007，22（02）：508－513．（In Chinese）

［10］王志刚，何展翔，魏文博，等．井地电法的准解析近似三维反演研究［J］．石油地球物理勘探，2007，42（2）：220－225．

WANG ZH G，HE ZH X，WEI W B，et al．3－D quasi－analytic approximate inversion of borehole－to－surface electric data［J］．OGP，2007，42（2）：220－225．（In Chinese）

［11］王志刚，何展翔，魏文博．积分方程法三维模拟井地电法并行算法研究［J］．物探化探计算技术，2007，29（5）：425－436．

WANG ZH G，HE ZH X，WEI W B．Parallel algorithm research of integral equation method to 3Dbore －hole surface electromagnetic modeling［J］．Computing techniques for geophysical and geochemical exploration，2007，29（5）：425－436．（In Chinese）

［12］吕玉增，阮百尧，黄俊革．直流电井间三维直接成像［J］．物探化探计算技术，2003，25（01）：60－64．

LU Y Z，RUAN B Y，HUANG J G．The 3－D immediate crosshole tomography with direct current［J］．2003，25（01）：60－64．（In Chinese）

［13］沈平，强建科，李永军，等．井间视电阻率的几何成像方法［J］．中南大学学报：自然科学版，2010，41（3）：1079 －1084．

SHEN P，QIANG J K，LI Y J，et al．Geometry image method of crosshole apparent resistivity［J］．Journal of Central South University（Sicence and Technology），2010，41（3）：1079－1084．（In Chinese）

［14］岳建华，刘志新．井－地三维电阻率成像技术［J］．地球物理学进展，2005，20（2）：407－411．

YUE J H，LIU ZH X．Three dimension resistivity tomography of mine ground［J］．progress in geophusics，2005，20（2）：407－411．（In Chinese）

［15］徐凯军，李桐林，张辉，等．基于共轭梯度法的垂直有限源三维电阻率反演［J］．煤田地质与勘探，2006，34（3）：69－71．

XU K J，LI T L，ZHANG H，et al．3Dresistivity inversion of vertical finite line source using conjugate gradients［J］．Coal geology ＆exploration，2006，34（3）：69 －71．（In Chinese）

［16］柯敢攀，黄清华．井地电法的三维正反演研究［J］．北京大学学报：自然科学版，2009，45（03）：264－272．

KE G P，HUANG Q H．3DForward and inversion problems of borehole－to－surface electrical method ［J］．Acta scientiarum Naturalium Universitatis Pekinensis2009，45（03）：264－272．（In Chinese）

［17］安然，李桐林，徐凯军．井地三维电阻率反演研究［J］．地球物理学进展，2007，22（1）：247－249．

AN R，LI T L，XU K J．Wellsurface 3Dresistivity inversion［J］．Progress in geophsycis，2007，22（1）：247 －249．（In Chinese）

［18］RUDIN L I，OSHER S，FATEMI E．Nonliner total variation based noise removal algorithms［J］．Phsica D：Nonliner Phenomena，1992，60（1－4）：259－268．

［19］WANG Y F，CUI Y，YANG C C．Hybrid regularization methods for seismic reflectiveity inversion［J］．Int，J．Geomath，2011，2（1）：87－112．

［20］ACAR R，VOGEL C R．Analysis of total variation penalty methods［J］．Inverse Problems，1994，（10）：1217－1229．

［21］徐世浙．地球物理中的有限单元法［M］．北京：科学出版社，1994．

XU SH ZH．The finite element method in geophysics ［M］．Beijing：Science Press，1994．（In Chinese）

［22］吴小平，徐果明，李时灿，等．利用不完全Cholesky共轭梯度法求解点源三维地电场［J］．地球物理学报，1998，41（06）：848－855．

WU X P，XU G M，LI SH C，et al．The calculation of three－dimensional geoelectric field of point source by incomplete cholesky conjugate gradient method［J］．Chinese Journal Geophysics．1998，41（06）：848－855．（In Chinese）

［23］TIKHONOV A N，ARSENIN V Y．Solution of ill－Posed Problems［M］．Washington．DC：V H Winston and Sons，1977．

［24］HANSEN，C．Rank－deficient and discrete ill－posed problems［M］．Numerical aspects of liner inversion：Department of mathematical modeling，Technical University of Denmark Lyngby，1998．

2．5D inversion of borehole-surface DC resistivity based on the total variation function

LAN Ze-luan，ZHANG Zhi-yong，DENG Ju-zhi，ZHOU Feng，LI Man
（Key Laboratory of Radioactive Geology and Exploration Technology Fundamental Science for National Defense，East China Institute of Technology，Nanchang 330013，China）

The 2．5Dregularization inversion for borehole－to－ground DC resistivity has been implemented in this paper，which is based on total variations（TV）to facing the problem that the ground DC resistivity exploration depth was limited．To start from the well satisfied boundary value problem，using the triangular unit finite element quadratic interpolation derived to 2．5Dforward modeling formulation．While the approximate boundary treatment，as well as the use of graph theory matrix rearrangement and fill element analysis method to achieve a direct solution of sparse matrices，greatly improving the computational efficiency．In order to improve the stability of the inversion results and the ability to distinguish the exception of boundaries，using L－curve to achieve regularization factor，stability factor of the TV and the conjugate gradient method for solving the objective function，which have shown that efficiency and accuracy．Finally，the inversion results in different borehole－toground observation way were on the basis of contrastive analysis to reach the anomaly characteristics and distribution rule of borehole－to－ground dc resistivity．

borehole-surface electrical method；conjugate gradient method；finite element method；total variation；L-curve

P 631．3

A

10．3969／j．issn．1001-1749．2015．03．02

1001-1749（2015）03-0273-07

2014-10-19 改回日期：2014-12-28

国家自然科学基金（41304055，41304056）；东华理工大学研究生创新基金（DYCA13026）；江西省研究生创新基金（YC2013－S177）

蓝泽鸾（1988－），女，硕士，主要从事电法数值模拟研究，E-mail：ZFECIT＠163．com。