文章编号：1004-0609(2014)03-0793-08

超高密度电法有限单元法正演与广义最小二乘反演

冯德山^{1, 2}，王鹏飞^{1, 2}，杨炳坤^{1, 2}

(1. 中南大学 地球科学与信息物理学院，长沙 410083；

2. 中南大学 有色金属成矿预测教育部重点实验室，长沙 410083)

摘  要：针对超高密度电法的全四极装置，采用三角形网格剖分和线性插值的有限单元法进行正演模拟，然后对该模拟结果开展广义最小二乘正则化反演，反演过程中对模型参数同时施加光滑模型和背景模型约束以提高反演稳定性。以一个低阻双层结构模型为例，与5种常规电法进行对比，发现超高密度电法具有数据信息量大、反演精度高的优点。对一个山谷地形下正方体低阻模型的超高密度电法数据进行正演和反演模拟，反演结果表明：该方法在较复杂地形条件下也能较好地刻画出异常体的形态及地下电性的分布，提高物探解释精度。

关键词：超高密度电法；全四极装置；有限单元法；正演模拟；广义最小二乘反演

中图分类号：P631　　     文献标志码：A

Finite element forward simulation and least square generalized inversion of ultra-high density resistivity method

FENG De-shan^{1, 2}, WANG Peng-fei^{1, 2}, YANG Bing-kun^{1, 2}

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,

Central South University, Changsha 410083, China)

Abstract: According to the full-four-pole array of ultra-high density resistivity method, the triangle mesh subdivision and the finite element method of linear interpolation were used for forward simulation, and then generalized least square regularization for inversion was used for the simulation results. At the same time, the smooth model and background model constraints were imposed in model parameters to improve the stability of inversion. Taking a low resistivity double-layer structure model as an example, ultra-high density resistivity method possesses the advantages of large volume data and high inversion accuracy compared with 5 kinds of conventional resistivity methods. The forward and inversion simulations were carried out on data of low resistance cube under valley terrain model using ultra-high density electrical method. The inversion results show that even under complex terrain conditions this method can portray the shape of abnormal body and underground electrical distribution well, thus, effectively improving the precision of geophysical interpretation.

Key words: ultra-high density resistivity method; full-four-pole array; finite element method; forward simulation; generalized least square inversion

高密度电法具有快速、便捷、准确等优点，几十年来，被广泛地应用于岩溶探测、断层破碎带追索和管线探测等领域，解决了很多工程实际问题，成为浅部勘探的主要地球物理方法^[1-3]。但伴随着工程勘探复杂化的情况，传统高密度电法的一些自身缺陷也逐渐显露出来，如同一地点因采用不同数据采集方式产生不同反演结果^[4-5]，未能充分利用电极阵造成有限数据信息量，不利于后期反演等^[6-7]。而超高密度电法是澳大利亚ZZ电阻率成像研发中心提出的一种全新电阻率法勘探方法^[8-9]。与传统高密度电法相比，它不受传统电极装置和数据采集方式的限制，在相同时间内，可采集该剖面电极间所有可能组合的电位信息。一次性布极所获得的信息量是传统高密度电法的数十倍以上，保证了后期数据处理及反演结果的精度。

鉴于超高密度电法诸多优点，它刚出现就受到工程界的广泛关注，如雷旭友等^[10]将超高密度电法应用于土洞、煤窑采空区和岩溶的勘探中，有效地提高了勘探精度；李红立等^[11]采用跨孔超高密度电法对花岗岩进行球状风化勘探实验，结果表明，此方法能准确地定位孤石发育位置；孙林^[12]将超高密度电法应用于煤矿砂体探测，可有效地指导露天矿体机械化开采；黄杰^[13]等将超高密度电法应用于破碎带追踪。

目前超高密度电法的相关研究大多数都局限于工程应用，采用该方法进行正演、反演计算的研究有关报道不多。为此，本文作者针对超高密度电法的全四极装置开展了正演模拟和反演效果分析，可为超高密度电法的进一步研究、应用与推广奠定理论基础。

1  基本原理

与常规电法类似，超高密度电法也是一种以岩石和矿石电性差异为物性基础的阵列加密型勘探方法。通过用供电电极AB向地下输入电流为I的供电电流，再测量所有电极MN上电流I在介质中产生的电位差ΔU。根据ρ=KΔU/I，求出相应的电阻率ρ。其中，K是与A、B、M和N的排列及间距有关的电极系系数。

超高密度电法采用多条多芯电缆，当采用四电极组合时，对所有测量电极MN进行数据采集，不再区分观测装置方式，一次采集可兼含偶极排列、温纳排列和施伦贝谢尔排列等，故称为“全四极装置观测方式”。在全四极装置中，以64根电极为例^[14]，将每个排列的64根电极分为奇数组32个(1，3，5，…，61，63)和偶数组32个(2，4，6，…，62，64)的2组电缆，则电极变换步骤如下：1) 确定供电电极AB，选定A为奇数组，B为偶数组；2) 采用供电电极奇偶配对的采集方式，先固定A电极为1号电极，然后依次变换偶数组电极B，即此时供电电极AB变化规律分别为1-2，1-4，…，1-64，当B电极选定后，在AB附近(最好在AB间的中部)选定一接地条件良好电极作为M极，剩余电极作为N极，组成多个MN测量电极同时进行数据采集；3) 依次变换奇数组电极A，直到奇数组最大，每次移动完电极A后，供电电极B和测量电极MN选定方式与步骤2)相同。在一次通电过程中同时测量其他电极所能组成的MN的电位差，这样可以得到最多61个电位(MN1，MN2，MN3，…，MN60，MN61)数据信息。而奇数组32个电极和偶数组32个电极互相配对作为供电电极，即有32×32=1024 (次)供断电过程，若在保证所有电极接触良好的条件下，每次供电可最多同时采集61个电位差数据，故总的数据量大约最多为32×32×61=62464，采集数据量超过常规方法40倍以上。

2  超高密度电法有限单元正演模拟

超高密度电法也是在基于人工直流电场的作用下，以岩石和矿石的电性差异为基础观测地下电场分布规律的勘探方法。假设地下介质的导电性沿走向无变化，则点电流源场各节点电位的计算所满足的变分问题为^[15]

     (1)

式中：Ω为模拟区域；Γ_∞为边界条件；δ为狄拉克函数；σ为电导率；I为供电电极的电流；r_A是电源点(x_A，z_A)至积分点的距离；k是傅氏变换波数；n边界处法线向量；和分别为零阶和一阶修正贝塞尔函数。

2.1  网格剖分

为了更精确地拟合地形和复杂地电体，并注意到岩矿石的物性参数变化特性，有限元网格剖分方式采用三角形单元，如图1所示。

2.2  三角单元线性插值

图1  有限元网格三角剖分及三角单元示意图

Fig. 1  Schematic diagrams of FEM triangle gridding (a) and triangle element (b)

在三角形单元上，如图1所示记三角形3个顶点按逆时针方向分别为1、2和3。每个顶点坐标分别为(x₁，z₁)、(x₂，z₂)和(x₃，z₃)，顶点的电位为U₁、U₂和U₃，电导率为σ₁、σ₂和σ₃。则三角单元中电位和电导率函数为

                     (2)

式中：

                         (3)

                                (4)

                                (5)

                             (6)

                           (7)

2.3  单元分析

令

                            (8)

                             (9)

                     (10)

                     (11)

则

               (12)

式中：

       (13)

                   (14)

                        (15)

       (16)

其中：c₁₁=(6，2，2)σ，c₁₂=(2，2，1)σ，c₁₃=(1，1，2)σ，c₂₁=(1，2，1)σ，c₂₃=(1，2，2)σ，c₃₁=(2，1，2)σ，c₃₂=(1，2，2)σ，c₃₃=(2，2，6)σ。

假设三角单元Δ的12边落在Γ_Δ上，由于无穷远边界离电源比较远，即cos(r_A，n)≈1，可将看作常数C，提至积分符号外，则有

                 (17)

式中：，，l₁₂是边界单元三角形的12边边长，，，，，其余。

2.4  总体合成

根据式(1)~(4)得到单元积分：

                   (18)

式中：。将单元列向量扩展成全体节点的列向量，其中N为节点总数，将3×3的系数矩阵拓展成N×N的矩阵K，将列向量拓展成N维列向量S。再将泛函F(U)对U求变分，并令其为0，则超高密度电法有限元正演模拟公式为

KU=S                                     (19)

由上可知，经单元总体合成的K矩阵为大型稀疏矩阵，为加快求解速度，避免占据大量内存，在集成系数矩阵K时，本文作者采用变带宽、一维数组方式只存储系数矩阵中非零元素，同时记录非零元素在矩阵中的位置。

3  超高密度电法反演

超高密度电法反演建立在正演的基础上，其原理是根据野外采集的数据或正演模拟得到与之对应的预测数据，将其与观测值在最小二乘法下构造一个误差函数进行比较，并通过不断修正模型参数使误差函数取得极小值，这样修正后的模型参数就是地下地质体的真实参数^[16]。考虑到在电法反演过程中电阻率变化较大，故采用对数值表示模型电阻率参数^[17]。超高密度电法的线性反演方程为

                                 (20)

将该方程进行正则化，其方程为

                            (21)

式中：Δd为数据残差矢量，其值等于实测视电阻率与模拟的视电阻率之间的对数值差值( ，i=1，2，…，n)；Δm为模型参数的改正向量(，j=1，2，…，m)；A为偏导数矩阵()，ρ_j为第j个网格节点的电阻率。

在最小二乘法反演过程中，不可避免地存在信息不足或观测误差，导致矩阵A的某些奇异值很小，使矩阵A^TA严重病态。为抑制方程的病态，对模型参数同时施加光滑模型^[18]和背景模型约束^[19]：

       (22)

式中：m为模型参数向量；m_b为背景模型；C为光滑度矩阵，可用一阶或者二阶平滑度方式进行计算。如果式(22)中C为单位矩阵，则表示仅施加了背景模型约束。

这样可以构造相应的目标函数Φ为

        (23)

式中：右边第一项为通常数据拟合差的最小二乘项，第二项为同时引入光滑模型和背景模型约束项。其中λ为正则化因子或拉格朗日乘子。将式(23)两端对Δm^T进行求导并令其为0，可得到最小二乘线性反演方程为

     (24)

式(24)也可等效为下面的线性方程组

                   (25)

从方程组(25)中求得的模型修改量Δm代入下式：

                          (26)

如此便可得到新的预测模型参数矢量m^(k)。重复上述操作直到实测数据和模拟数据之间的平均均方根误差满足要求为止。其中，平均均方根误差R_MS定义为

                          (27)

4  模型反演实例

在实际的超高密度电法勘探中，如果61道数据同时采集，则n根电极将获得61n²/4组数据，可见电极数与数据量呈平方关系。为了减少计算量和节约反演时间并能说明超高密度电法的准确性、实效性及强抗干扰能力，本文作者选用48根电极进行模拟。设计一个低阻正方体上覆低阻地层的结构模型，采用超高密度电法全四极装置采集方式进行正演模拟及反演成像，并将该反演结果与5种常规排列装置进行效果对比分析。

4.1  低阻混合结构模型

图2所示为低阻混合结构模型示意图。如图2所示，模型背景电阻率ρ₀=100 Ω·m；模型电阻率分别为ρ₁=5 Ω·m和ρ₂=10 Ω·m。图3所示为超高密度电法全四极装置数据采集方式电阻率反演图；图4~8分别为偶极排列、α排列、β排列、γ排列和施伦贝谢尔排列电阻率反演图，其反演图蓝色方框(见网络电子版，下同)对应模型投影区域。

图3所示为超高密度电法全四极装置电阻率反演图。分析图3可知，红色区域为低阻地带，电阻率为5~10 Ω·m，绿色区域为背景电阻率分布，电阻率约为100 Ω·m，该反演结果在表达浅部低阻薄层、低阻异常体大小、空间位置、电阻率大小等参数上与模型非常一致，说明超高密度电法反演数据精度极高；而在图4所示偶极排列电阻率反演图中，浅部的低阻薄层在图中反映并不明显，矩状低阻异常体位置虽有反映，但异常体的位置较真实位置明显上移；再分析图5~8可知：α排列、β排列、γ排列和施伦贝谢尔排列与偶极排列的电阻率反演结果类似，浅部低阻薄层在图中反映并不明显，尽管它们都能有效指示矩状低阻异常体的存在，但不能准确指明异常体的大小和位置。对比图3和4~8可知，相对其他常规排列装置超高密度电法，超高密度电法能够更有效地反映地下的电性分布状况，达到准确探测的效果。

图2  低阻混合结构模型示意图

Fig. 2  Schematic diagram of low resistivity mixed model

图3  超高密度电法全四极装置电阻率反演图

Fig. 3  Map of electronic resistivity inversion with full-four-pole array of ultra-high density resistivity method

图4  偶极排列电阻率反演图

Fig. 4  Map of electronic resistivity inversion with dipole-dipole array

图5  α排列装置电阻率反演图

Fig. 5  Map of electronic resistivity inversion with α array

图6  β排列装置电阻率反演图

Fig. 6  Map of electronic resistivity inversion with β array

图7  γ排列装置电阻率反演图

Fig. 7  Map of electronic resistivity inversion with γ array

图8  施伦贝谢尔排列电阻率反演图

Fig. 8  Map of electronic resistivity inversion with γ Schlumberger array

4.2  山谷结构模型

如图9所示，模型背景电阻率ρ₀=100 Ω·m，极化率η₀=1.0%；模型边长为4 m×4 m，中心位置为(24 m，-7 m)，模型顶部与地表相距3 m，模型电阻率ρ₁=5 Ω·m，极化率η₁=5.0%。图10所示为超高密度电法全四极装置电阻率反演图。其中，反演图蓝色方框对应模型投影区域。

图9  山谷结构模型示意图

Fig. 9  Schematic diagram of valley structure model

图10  山谷结构模型全四极装置电阻率反演图

Fig. 10  Map of electronic resistivity inversion with full-four-pole array for valley structure model

在起伏地形条件下分析图10，绿色区域为背景电阻率分布情况，电阻率约为10 Ω·m，红色区域与蓝色方框吻合得较好，为正方形低阻体，电阻率为5~15 Ω·m，虽然在模拟的背景和低阻正方体接触带有一约为35 Ω·m的过度带，但总体上全四极装置的电阻率反演图能很好地反映正方体的中心埋深，对模型异常体的大小和深度也能很好地拟合。可见，在山谷地形条件下，全四极装置也能较好地对地质体进行勘探，清晰地反映地下电场特征。

5  结论

1) 超高密度电法改变常规电法勘探中数据采集方式，采用多电极、多通道、全自动采集技术，野外采集一次性布极，大大提高了工作效率，且可采集到几十倍于常规电法的数据，克服了在同一地点因不同常规数据采集方法导致不同反演结果的缺点，提高了反演精度。

2) 以基于三角形剖分的有限单元法进行高密度电法正演模拟，能够有效地拟合地形和复杂地电体，正演结果正确可靠；以广义最小二乘法对超高密度电法进行反演成像，结果表明超高密度激电法比单一排列装置反演精度高，反演效果好，其探测效果优于单一排列装置的。

REFERENCES

[1] 李金铭. 地电场与电法勘探[M]. 北京: 地质出版社, 2005: 208-215.

LI Jin-ming. Electric field and electrical prospecting[M]. Beijing: Geological Publishing Heuse, 2005: 208-215.

[2] 钟 韬, 邓艳平. 高密度电法在西部岩溶地区勘探中的应用[J]. 工程地球物理学报, 2009, 6(S1): 80-85.

ZHONG Tao, DENG Yan-ping. The applied research in the western karst region by high density resistivity method[J]. Chinese Journal of Engineering Geophysics, 2009, 6(S1): 80-85.

[3] 张延凯, 王 章, 吴 迪, 蔡嗣经. 高密度电法在姑山矿防渗帷幕工程中的应用[J]. 有色金属: 矿山部分, 2010, 62(2): 67-70.

ZHANG Yan-kai, WANG Zhang, WU Di, CAI Si-jing. Application of high density resistivity method in the watertight screen project in Gushan mine[J]. Nonferrous Metals: Mining Section, 2010, 62(2): 67-70.

[4] 刘 飞. 高密度电阻率法正反演数值模拟研究[D]. 北京: 中国地质大学, 2012.

LIU Fei. The research of high density resistivity method of positive and negative numerical modeling[D]. Beijing: China University of Geosciences, 2012.

[5] PAPADOPOULOS N G, TSOKAS G N, DABAS M, YI M J, KIM J H, TSOURLOS P. Three-dimensional inversion of automatic resistivity profiling data[J]. Archaeological Prospection, 2009, 16(4): 267-278.

[6] 李美梅. 高密度电阻率法正反演研究及应用[D]. 北京: 中国地质大学, 2010.

LI Mei-mei. Research on forward modeling and inversion of multi-electrode resistivity method[D]. Beijing: China University of Geosciences, 2010.

[7] 陈 宁, 雷 宛, 肖宏跃, 吕志强, 胡承林. 边界系数法在高密度电阻率法正反演解释中的应用效果研究[J]. 成都理工大学学报: 自然科学版, 2011, 38(3): 284-290.

CHEN Ning, LEI Wan, XIAO Hong-yue, Zhi-qiang, HU Cheng-lin. Application of the boundary coefficient method to the explanation of high-density resistivity method’s forward modeling and inversion[J]. Journal of Chengdu University of Technology: Science & Technology Edition, 2011, 38(3): 284-290.

[8] ZHE J P, GREENHALGH S, MARESCOT L. Multi-channel, full waveform and flexible electrode combination resistivity imaging system[J]. Geophysics, 2007, 72(2): 57-64.

[9] GREENHALGH S A, ZHOU B, GREEN A G. Solutions, algorithms and inter-relations for local minimisation search geophysical inversion[J]. Journal of Geophysics and Engineering, 2006, 3: 101-113.

[10] 雷旭友, 李正文, 折京平. 超高密度电阻率法在土洞、煤窑采空区和岩溶勘探中应用研究[J]. 地球物理学进展, 2009, 24(1): 340-347.

LEI Xu-you, LI Zheng-wen, ZHE Jing-ping. Applications and research of the high resolution resistivity method in exploration of caves, mined regions and karst region[J]. Progress in Geophysics, 2009, 24(1): 340-347.

[11] 李红立, 张 华, 汪传斌. 跨孔超高密度电阻率法在花岗岩球状风化体勘探中的试验研究[J]. 工程勘察, 2010, 8: 88-92.

LI Hong-li, ZHANG Hua, WANG Chuan-bin. Experimental study on the cross-hole ultra-density resistivity method used in the exploration for the spheric lightly-weathered granite[J]. Geotechnical Investigation & Surveying, 2010, 8: 88-92.

[12] 孙 林. 超高密度电法在煤矿砂体探测中的应用[J]. 煤矿开采, 2012, 17(3): 25-27.

SUN Lin. Application of high density electrical method in detecting sandbody in coal mine[J]. Coal Mining Technology, 2012, 17(3): 25-27.

[13] 黄 杰, 钟 韬, 马文德. 超高密度电法在追索破碎带中的应用[J]. 物探化探计算技术, 2009, 31(6): 586-589.

HUANG Jie, ZHONG Tao, MA Wen-de. Ultra-high density resistivity method applied to exploration of fracture zone[J]. Computing Technology for Geophysical and Geochemical Exploration, 2009, 31(6): 586-589.

[14] 钟 韬. 超高密度电法在探测采空区中的应用研究[D]. 成都: 成都理工大学, 2008.

ZHONG Tao. The applied research into the exploration of the gob by multi-channel full waveform and flexible electrode combination resistivity imaging system[D]. Chengdu: Chengdu University of Technology, 2008.

[15] 戴前伟, 肖 波, 冯德山, 刘海飞, 王鹏飞. 基于二维高密度电阻率勘探数据的三维反演及应用[J]. 中南大学学报: 自然科学版, 2012, 43(1): 293-300.

DAI Qian-wei, XIAO Bo, FENG De-shan, LIU Hai-fen, WANG Peng-fei. 3-D inversion of high density resistivity method based on 2-D exploration data and its application[J]. Journal of Central South University: Science and Technology, 2012, 43(1): 293-300.

[16] MUKANOVA B, ORUNKHANOV M. Inverse resistivity problem: Geoelectric uncertainty principle and numerical reconstruction method[J]. Mathematics and Computers in Simulation, 2010, 80(10): 2091-2108.

[17] RODI W, MACKIE R L. Non-linear conjugate gradients algorithm for 2-D magnetotelluric inversion[J]. Geophysics, 2001, 66(1): 174-187.

[18] MURAT B, LEYLA S. Two-dimensional resistivity imaging in the Kestelek boron area by VLF and DC resistivity methods[J]. Journal of Applied Geophysics, 2012, 82: 1-10.

[19] LIN C H, TAN H D, TONG T. Parallel rapid relaxation inversion of 3D magnetotelluric data[J]. Applied Geophysics, 2009, 6(1): 77-83.

(编辑  陈卫萍)

基金项目：国家自然科学基金资助项目(41074085)；教育部新世纪优秀人才支持计划项目(NCET-12-0551)；湖湘青年创新创业平台培养对象；中南大学升华育英计划项目资助

收稿日期：2013-03-20；修订日期：2013-10-28

通信作者：冯德山，副教授，博士；电话：0731-88836145；E-mail：fengdeshan@126.com