基于二维数据的加权非线性共轭梯度三维反演

韩 雪，强建科，鲁 凯，李俊营，满开峰，毛先成



基于二维数据的加权非线性共轭梯度三维反演

韩 雪1，强建科2,3，鲁 凯2，李俊营2，满开峰2，毛先成2,3

( 1.广东省地质物探工程勘察院，广东 广州 510800；2.中南大学 地球科学与信息物理学院，湖南 长沙 410083；3.中南大学 有色金属成矿预测教育部重点实验室，湖南 长沙 410083)

针对当前电阻率数据的解释仍以二维反演为主，而地下为三维结构的情况，研究起伏地形条件下三维源二维采集方式获取的多条测线测深数据的三维反演问题。正演算法采用三棱柱剖分的有限单元法，可以模拟起伏地形下的复杂模型；在求解大型线性方程组时采用了不完全乔利斯基分解回代技术，大大提高了三维多点电源的正演计算速度；在计算偏导数矩阵过程中，依据互换原理，只需有效组合正演时每个点电源所对应的节点电位便可得到偏导数矩阵，节约大量计算时间。其次三维反演算法采用加权正则化共轭梯度法，抑制了严重病态的非线性问题。模型计算结果表明：加权非线性共轭梯度反演算法既保证了反演稳定收敛，又能适合起伏地形反演，达到了满意的效果。

电阻率三维反演；起伏地形；共轭梯度反演；二维数据三维反演

1 引 言

直流电阻率法在探测岩溶、采空区等方面应用广泛[1,2]，然而在实际勘查中经常会遇到地形起伏情况，实测数据会发生严重的畸变，这对后期数据的反演解释工作带来了非常大的困难。对于三维反演，偏导数矩阵的计算时间和计算精度，直接影响反演计算的时间和精度，因此，偏导数矩阵求取是反演问题的关键之一。在电阻率一维反演中，偏导数矩阵常常采用差分法近似获得，由于模型参数少，即使多次调用正演程序，总体速度仍然很快，但对于二、三维反演来说，模型参数很多，差分法显然是不合适的。为了快速方便地求出精度较高的偏导数矩阵，国内外的地球物理学家进行了不懈的努力。Tripp等(1984)[3]提出了一种利用电位函数与模型参数间的简单关系来计算偏导数的近似方法，当时是用在二维地电断面的反演问题中，取得了很好的效果。以后，很多地球物理学家的偏导数计算大都采用这个办法。阮百尧等(1999，2001)[4,5]对二维偏导数矩阵的求取做了详细论述；黄俊革(2003)[6]利用互换定律在三维方面的应用作了研究，提高了计算速度；Papadopoulos等(2011)[7]改进了雅可比矩阵的计算量，去除对反演没有重要贡献的参数，在保证精度的基础上，减少内存，提高反演速度。本文采用互换定律来计算偏导数矩阵。

Park等(1991)[8]发表了有限差分的三维反演算法，开始了比较系统的反演方法研究。Sasaki(1994)[9]阐述了基于有限单元法的三维电阻率反演，但对深部单元的分辨能力不强，反演得到的电阻率值与实际模型有较大差别。Zhang等(1995)[10]使用共轭梯度算法对三维模型进行了正反演研究。吴小平等(1999，2000，2005)[11-13]、刘海飞等(2011)[14]研究了E-SCAN测量方式在非平坦地形上的三维正反演问题，取得了一些好的效果，但在实际工程勘探中，受成本和周期的影响，获取三维勘探原始数据存在较大困难，而且E-SCAN测量方式是基于二极电位测量装置的，分辨率低于梯度装置，目前野外勘探仍然以多剖面的二维测量为主。黄俊革等(2004)[15]对三维电阻率反演作了改进，提出“体积因子”、将全局反演改为局部反演等措施，提高了反演的速度和精度。强建科等(2007，2013)[16,17]提出三棱柱单元剖分的有限元方法可以模拟大倾角的复杂地形，同时改进了偏导数矩阵，解决了反演深度偏浅的弊病。Gunther等(2006)[18]实现了任意地形的三维电阻率反演，正演采用的是非结构化网格技术的有限单元法。Oldenborger等(2009)[19]利用点扩散函数实现了三维电阻率成像反演。

本文研究的是目前野外勘探常采用的三维源二维采集方式获取的三维数据，依此数据进行三维起伏地形反演算法。

2 三维反演思路

研究反演的第一步是正演问题。本文采用三维起伏地形有限元正演算法[16]，地下空间采用三棱柱(五面体)单元剖分，这种单元可以模拟起伏地形，插值函数选为三线性插值函数，外边界采用混合边界[20]。

在实现三维反演中，采用固定网格大小，只考虑视电阻率对地下各个网格电导率的响应关系，可以减少未知参数数量，便于求解。

在通常情况下，直流电阻率法的正演问题可以表示成如下形式，即

A(m)=d

(1)

式(1)中：A是非线性正演算子，这里用有限元法实现；m是模型参数电导率，单位为S/m；d是观测数据，即视电阻率值，单位为Ω·m。

而反演过程就是根据非线性正演算子A和观测数据矢量d求出模型参数m。

m=A-1d

(2)

由于式(2)是病态的，即解是非惟一和不稳定的，因此要根据正则化原理求解病态问题，很多学者将正则化反演方法应用于直流电的三维反演[21-27]，目标函数方程如下：

Pα(m,d)=φ(m)+αS(m)=min

(3)

rn=A(mn)-d

(4)

(5)

(6)

(7)

(8)

3 偏导数矩阵

在共轭梯度求解中，需要知道偏导数矩阵，即视电阻率对模型单元电导率的灵敏度矩阵

F=[aij]

其中i=1,2,…,Ns,

j=1,2,…,Mm

(9)

式(9)中，ρsi是实测得到的视电阻率数据，共有Ns个；σj是第j个网格电导率，共有Mm个模型参数。

以上偏导数矩阵计算，可归结到网格节点电位对模型中各个网格电导率的偏导数矩阵计算问题。

例如，对于偶极-偶极装置来说，视电阻率计算公式为

(10)

对上式两端求电导率偏导数，有

(11)

由此看出，任何排列视电阻率对模型电导率的偏导数都可用点电流源的电位对电导率的偏导数来表达。

在正演计算时，有限元的单元集成方程为

KU=P

(12)

式(12)中，K为系数矩阵，U为单元节点电位，P为电源项和边界条件。

在方程(12)两边，对地下节点K标识的网格单元电导率σK求偏导数，由于场源项P与地下电导率无关，故求导结果为

(13)

比较方程(12)和(13)，可将方程(13)理解为，地面节点电位U对任意网格单元电导率σK的偏导数，可以由方程(13)右端项(新场源项)求得，而该新场源项仅与该单元(三棱柱)矩阵系数和6个节点电位有关，因此，把新场源项可以看作位于这6个节点上的6个微小点电流源：

(14)

式(14)中Kij(i=1，2，…，6；j=1，2，…，6)为三棱柱单元的刚度矩阵之诸元素。当对求解区的网格进行剖分时，所用步长是相等的，即求解区采用均匀剖分，这时，式(14)的系数矩阵对于所有单元都是相同的。这给计算带来了方便。

(15)

将(14)式代入(15)式，有：

(16)

有了地面节点的电位对模型电导率的一阶偏导数，按照(11)式可以组合出偶极装置的视电阻率对电导率的偏导数矩阵F，即雅可比矩阵。获得雅可比矩阵F后，就可进行反演成像计算了。

从上述计算步骤可以看出，在形成雅可比矩阵的过程中，其计算量是很小的，因为上述步骤中用到的不同供电点对应的所有网格节点电位，都来自于正演计算中保留下来的；单元刚度矩阵数据也是正演时保存的。计算雅可比矩阵只是利用了正演的计算结果，通过换算得到的。当然，如果只需地面节点电位的正演结果，不需要保存整个断面空间的电位函数值。因此，为了计算雅可比矩阵需要更多的计算机内存空间以保存每一个地面电极供电时的断面空间的电位值，而这对于现在的计算机条件来说，是可以做到的。也就是说，现在只需作一次正演就可以得到雅可比矩阵，使得反演所需的时间大大地减少了，这种求雅可比矩阵的方法使得高密度电法的三维反演成像成为可能，为保障高密度电法的地质效果提供了条件。

4 三维反演算例

下面的算例中均采用偶极—偶极装置。

模型1：简单模型如图1(a)所示，均匀半空间中有一个4 m×4 m×4 m的低阻立方体模型，电阻率为10 Ω·m， 顶埋深为3 m， 中心埋深5 m，围岩电阻率为100 Ω·m，测量点距1 m，线距1 m，每条剖面上25个供电点，共设计5条剖面。通过正演获得“实测数据”。在反演中，初始模型电阻率取“实测视电阻率数据”的平均值作为背景电阻率。

图1(b)是偶极-偶极装置三维正演垂直切片图，可以看出5条视电阻率等值线异常呈标准的“八字形”低阻特点，且前后左右对称。在低阻体正上方三条测线(Y=-1 m、Y=0 m、Y=1 m)视电阻率差异较大，两边测线(Y=-2 m、Y=2 m)视电阻率差异变小。

图1(c)是三维反演结果，经过约1小时8次迭代，拟合误差小于5 %，这5条断面图很好地还原了低阻体的大小、形状和埋深，说明该算法的收敛性很好。

模型2：复杂模型如图2(a)所示，均匀半空间中有一个4 m×4 m×4 m的高阻立方体模型和一个长宽高为6 m×4 m×6 m低阻长方体模型，高阻立方体电阻率为500 Ω·m，顶埋深为2 m，低阻立方体电阻率为10 Ω·m，顶埋深为4 m，两个异常相隔4 m，围岩电阻率为100 Ω·m，测量点距1 m，线距1 m，每条剖面上31个供电点，共设计5条剖面。模拟“实测数据”由三维正演得出。初始模型电阻率取“实测视电阻率数据”的平均值作为背景电阻率。

图2 复杂模型三维反演结果Fig.2 3D inversion result of complex model

图3 起伏地形低阻体三维正反演结果Fig.3 The result of forward calculation and inversion under undulating terrain 3D low resistance body model

对上述复杂模型模拟的视电阻率数据进行三维反演，评估该算法对高低阻组合模型的分辨能力。图2(b)是三维反演结果水平切片，沿纵深方向取4个切片，从图中可看出，在Z=3 m切片上，左侧有一个高电阻率异常，其异常中心与真实模型一致；在Z=5 m、7 m切片中，右侧有一个低电阻率异常，其异常中心与真实模型一致。由此得出，该反演算法能够从复杂的组合视电阻率异常中反演出不同深度的高、低阻体，成像效果明显。

图2(c)是三维反演结果纵向切片，从5条断面图上看出，在X方向上，该算法能够清楚地识别间隔4 m的两个高低电阻率异常，且异常中心与真实模型一致，成像效果很好。图2(d)对Y=0 m断面和Z=5 m切片上实测数据与预测模型数据进行对比，二者视电阻率等值线形态一致，说明预测模型与真实模型非常接近。

模型3：在模型1的基础上增加一个三维山脊地形，如图3(a)所示，坡度约30°，斜坡长8 m，山脊高4 m，其他参数同模型1。

图3(b)为5条测线上的视电阻率等值线断面图，从图上看出，三维山脊对视电阻率异常产生了非常严重的干扰，致使山脊下视电阻率值变得很大，其视电阻率为围岩电阻率的3倍多，两侧为低阻异常。与图1(b)比较看出，山脊产生的虚假高视电阻率异常完全破坏了低阻体产生的“八字形”视电阻率异常特征。

图3(c)为起伏地形三维反演结果，从图中看出，反演算法正确反演出了低阻体的位置、形状和大小，反演结果与地下真实模型相符，没有得出多余的构造。由此得出，直接带入地形进行三维反演才能够有效地解决由地形产生的虚假异常，从而得到地下真实的电性分布特征，同时证明该算法反演成像稳定可靠。

5 结 论

偏导数矩阵求解是电阻率三维反演的一个关键问题，根据稳定电流场的互换原理，推导出含有起伏地形信息的雅可比矩阵公式，通过组合各节点的电位，能够快速得到雅可比矩阵，节约大量计算时间。在前人研究的基础上，提出加权正则化非线性共扼梯度方法，优化了权重系数和正则化因子，从而实现起伏地形电阻率三维反演算法，算例表明该方法成像效果良好。在求解大型线性方程组时采用了不完全乔利斯基分解回代技术，大大提高了三维多点电源的正演计算速度问题。

[1]肖敏,陈昌彦,白朝旭,等.北京地区浅层采空区高密度电法探测应用分析[J].工程地球物理学报,2014,11(1):29-35.

[2]袁忠明,韦乙杰.高密度电阻率法偶极装置在某道路工程岩溶探测中的应用[J].工程地球物理学报,2015,12(1):72-76.

[3] Tripp A C, Hohmann G W, Swift C M. Two-dimensional resistivity inversion[J].Geophysics,1984,49 (10):1 708-1 717.

[4]阮百尧,村上裕,徐世浙.电阻率激发极化法数据的二维反演程序[J].物探化探计算技术,1999,21(2): 116-125.

[5]阮百尧.视电阻率对模型电阻率的偏导数矩阵计算方法[J].地质与勘探,2001,37(6):39-41.

[6]黄俊革.三维电阻率/激化率有限单元正演模拟与反演成像[D].湖南:中南大学,2003.

[7]Papadopoulos N G, Tsourlos P, Papazachos C, et al.An algorithm for fast 3D inversion of sur-face electrical resistivity tomography data: application on imaging buried antiquities[J]. Geophysical Prospecting,2011,59(3):557-575.

[8]Park S K, Van G P.Inversions of pole-pole data for 3-D resistivity structure beneath arrays of electrodes[J].Geophysics,1991,56(7):951-960.

[9]Sasaki Y.3-D resistivity inversion using the finite-element method[J]. Geophysics,1994,59(12): 1 839-1 848.

[10]Zhang J, Mackie R L, Madden T R. 3-D resistivity forward modeling and inversion using conjugate gradients[J]. Geophysics,1995,60(5):1 313-1 325.

[11]吴小平,徐果明.电阻率三维反演中偏导数矩阵的求取与分析[J].石油地球物理勘探,1999,34(4): 363-373.

[12]吴小平,徐果明.利用共轭梯度法的电阻率三维反演研究[J].地球物理学报,2000,43(3):420-427.

[13]吴小平.非平坦地形条件下电阻率三维反演[J].地球物理学报,2005,48(4):932-936.

[14]刘海飞,柳建新,郭荣文,等.起伏地形三维激电连续介质模型快速反演[J].吉林大学学报(地球科学版),2011,41(4):1 212-1 218.

[15]黄俊革,阮百尧,鲍光淑.基于有限单元法的三维地电断面电阻率反演[J].中南大学学报(自然科学版), 2004,35(2):295-299.

[16]强建科,罗延钟.三维地形直流电阻率有限元法模拟[J].地球物理学报,2007,50(5):1 606-1 613.

[17]Qiang Jianke, Han Xue, Dai Shikun.3D DC Resistivity Inversion with Topography Based on Regularized Conjugate Gradient Method[J]. International Journal of Geophysics,2013.Article ID 931876, 9 pages, 2013. doi:10.1155/2013/931876.

[18]Günther T, Rücker C, Spitzer K. Three-dimensional modelling and inversion of dc resistivity data incorporating topography-II. Inversion[J]. Geophys.J.Int.,2006,166(2):506-517.

[19]Oldenborger G A, Routh P S. The point-spread function measure of resolution for the 3-D electrical resistivity experiment[J]. Geophys.J.Int.,2009,176(2):405-414.

[20]Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped three-dimensional structures [J]. Geophysics,1979,44(4):753-780.

[21]Ellis R G, Oldenburg D W.T he pole-pole 3-D DC resistivity inverse problem: a conjugate gradient approach[J]. Geophys.J.Int.,1994,119(1):187-194.

[22]LaBrecque D J, Miletto M, Daily W, et al. The effects of noise on Occam's inversion of resistivity tomography data[J]. Geophysics,1996,61(2):538-548.

[23]Yi M J, Kim J H, Song Y, et al. Three-dimensional imaging of subsurface structures using resistivity data[J]. Geophysical Prospecting,2001,49(4):483-497.

[24]Pain C C, Herwanger J V, Worthington M H, et al. Effective multidimensional resistivity inversion using finite-element techniques[J]. Geophys.J.Int.,2002,151(3):710-728.

[25]Pidlisecky A, Haber E, Knight R. RESINVM3D: A 3D resistivity inversion package[J]. Geophysics, 2007,72(2):H1-H10.

[26]Marescot L, Palma Lopes S, Green A G. Nonlinear inversion of geoelectric data acquired across 3D objects using a finite-element approach[J]. Geophysics,2008,73(3):F121-F133.

[27]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.

[28]Zhdanov M S. Geophysical inverse theory and regularization problems[M].New York:Elsevier, 2002.

The 3D Inversion with Non-linear Conjugate Gradient Based on 2D DC Resistivity Data

Han Xue1，Qiang Jianke2,3，Lu Kai2，Li Junying2，Man Kaifeng2，Mao Xiancheng2,3

( 1.Geologic&GeophysicalEngineeringExplorationInstituteofGuangdong,GuangzhouGuangdong510800,China;2.SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,ChangshaHunan410083,China;3.KeyLaboratoryofMetallogenicPredictionofNonferrousMetals,MinistryofEducation,CentralSouthUniversity,ChangshaHunan410083,China)

In view of the current resistivity data interpretation, the 2D inversion is still taken as the main point, but actually the 3D structure is used more underground. This paper presents a research on the 3D inversion for the multiple line sounding data over rugged terrain area. Triangular-prism mesh finite element method is employed in forward algorithm, which can be used to simulate the complex models with the topography. In solving large linear equations, incomplete Cholesky decomposition and back substitution technique are used, which can significantly improve the speed of 3D multi-source modeling. According to the interchangeable principle, the partial derivative matrix could be acquired by combining the potential of forward simulation, which saves large amount of time. The weighted regularized conjugate gradient is applied to 3D inversion, which solves the severely ill-posed non-linear problem. The model calculation results show that this method could both guarantee the stable inversion convergence and be suitable for rugged topography inversion to reach satisfactory results.

3D resistivity inversion; rugged terrain; conjugate gradient inversion; 3D inversion with 2D data

1672—7940(2016)05—0561—09

10.3969/j.issn.1672-7940.2016.05.001

国家自然科学基金(编号：41472301；41174104))

韩 雪(1988-)，女，助理工程师，硕士，主要从事电法勘探研究。E-mail:kahanxue@126.com)

强建科(1967-)，男，副教授，博士，主要从事电磁法正反演研究。E-mail:qiangjianke@163.com

P631.3 文献标码： A

2016-05-28