激电二维数值模拟研究

陈永凌，蒋首进，谢 丹

（1．武警黄金十二支队，成都 610059；2．武警黄金十支队，昆明 65000）

激电二维数值模拟研究

陈永凌1，蒋首进2，谢 丹1

（1．武警黄金十二支队，成都 610059；2．武警黄金十支队，昆明 65000）

从点源二维地电问题出发，采用有限单元法进行了地电场进行数值模拟，采用自适应三角剖分来实现起伏地表的模拟 ，针对双边三极装置，实现了多种模型的正演研究；通过多种模型的正反演，总结异常产生的规律，为激电法的分析提供了有效的信息。

激发极化法；有限元法；双边三极

0 引言

有限元方程是由Coggon在20世纪70年代由电磁场能量最小原理导出的，通过计算点源二维地点问题的变分方程，得到其相应的解。20世纪80年代，颇瑞德摩A．B研究了三维地电有限元数值模拟。二十一世纪初，Matesihat计算了二维模型使用为函数场源的三维场源的响应。我国从1975年开始在电法勘探中引入有限元方法，取得了很好的效果。上世纪80年代钟本善等［7］发表了基于点源二维地电的有限元正演数值模拟以及数值模拟过程中的若干问题，进行了有限单元正演模拟，推导了地电场模拟的有关方程；徐世浙［6］通过网格内电导率连续变化，采用矩形网格剖分对地电场进行了数值模拟；阮百尧［8］提出电导率和电位双线性插值，在此过程中采用线性变化的电导率分块，有效地提高了有限元正演数值模拟的精度；阮百尧［8］针对结构单元电导率等参数与实际不相符合的问题，相继提出双线性插值的模拟方法和连续线性变化数值模拟，并且优化了系数矩阵，十分有效地提高了计算精度和速度。

作者在前人的基础上，运用有限单元法进行了数值模拟，采用自适应三角剖分来实现起伏地表的模拟，通过稳定电流场的基本方程的导出、点源二维地电问题的讨论、边值问题的讨论等。在这过程中采用双二次插值基函数和Galerkin法将微分方程离散化，利用代数多重网格法求解方程组，波数域—空间域的傅立叶反变换，采用分段波数一维连分法直接计算反变换中的定积分［5］。

1 基本原理

1．1 点源二维地电问题

对于二维地电问题，我们通常假设平行于地质体走向为y轴，其倾向方为x轴，向下方向为z轴，由于地质体走向方向物性条件不变，则地电参数σ只随x、z的方向发生变化，在y方向上不发生变化，即

这时候电位的基本微分式方程变为

由于在二维情况下地电参数σ只随x、z的方向发生变化，在y方向上不发生变化，为简化计算条件，我们使用傅氏变换

将坐标系空间中的电位U变换为二维波数域中的电位φ，式中k为空间波数。

对式（3）进行傅氏变换可以得到

式（4）就是在点源二维地电断面时变换后的电位φ必须满足的基本微分方程式。通过上述变换将地电三维问题转换为二维来描述，其中k值只是一个参数，实际工作中k值为一常数，由不同的k值来求得相应的φ值，再通过傅氏反变换来求取未知电位U。

1．2 稳定电流场的边值问题

流场的分布是在整个地下介质空间分布的，而我们需要计算的只是其中的一部分空间，所以需要在区域边界上给傅氏电位φ进行边界条件的计算。

边界条件同常分为三类：

其中：第一类边界条件又叫“狄里希莱条件”；第二类边界条件又叫“诺依曼条件”；第三类边界条件又叫“混合边界条件”。

上面第三类边界条件为修正后的第三类边界条件，式中λ为边界修正参数，η由式（5）确定。

式中：lλ为λ点的供电电流强度；rλ为测点到场源距离为λ点的点矢径为边界处法线向量。

由于狄里希莱条件会产生计算值比解析解的值小；诺依曼条件会产生计算值比解析解的值大。混合边界条件具有能在远处边界面上保持电位的物理特性，且又不需要首先段设一个“正常”的电场分布的优点，这样可以减少边界附近放延展网格的数量，从而减少了计算的工作量。同时因为在边界上保持电位的连续，没有电性的突变，对地下边界的反射作用也消除。由于混合边界条件具有较高的模拟精度，较好的反演效果，所以一般情况下在计算过程中采用混合边界条件。

1．3 地电二维有限元法基本原理

从上面的推导中得到了点源二维目标函数

的极小值相对应。

其中：φ为目标函数；D为研究区域；Γ为研究区域D的边界。

1．4 网格剖分示意图

在使用有限元法求解研究区域的时候，通常需要将研究区域剖分成许多连续的小单元。网格剖分种类繁多，通常有以下几种剖分类型：①非结构化三角网格剖分；②矩形三角剖分；③四边形剖分；④三角形剖分。四边形剖分相对其他剖分较简单，且收敛速度较快，但四边形剖分对地形起伏及复杂模型模拟较为粗超，不能较好地模拟其真实的形态，不适应于起伏地形条件及复杂条件的地电模型，且在不同方向上电性结构变化也不相同；非结构化三角网格剖分具有三角单元的特征，对各种约束条件满足，但由于网格数量非常大，导致计算速度较慢；三角网格剖分具有满足地形起伏条件，可以很好地模拟异常体形态，满足地电模型的约束条件和三角单元的特性，且网格节点数量适合，计算速度较快，因而成为一种常用的网格剖分方法。

结合以上结论，采用了如图1所示的方式进行网格剖分，解决了网格大小与精度以及计算量之间存在的一系列矛盾，也较好地适应了地形的起伏。

网格剖分的原则一般来讲有以下几点：

1）各个单元的节点只能与相邻的单元节点相重合，但是不能成为相邻的单元的边上或者边内的点。

2）每个单元内介质的电导率为常数，在不同电导率介质的分界面或者求解区的边界上，使用三角形的一个边互相连接，以模拟边界线的近似折线取代边界线。

3）三角形的三条边或者三个顶点的差别不能太大。

4）在某些区域物性比变化剧烈，点位参数的二阶或高阶导数数值较大的区域，单元剖分应剖分细密一些，而在电场变化平缓的地方，网格可以分稀疏些［5－6］。

图1 网格剖分示意图Fig．1 The grid subdivision schemes

2 模型计算

本节主要对固定点源测深采用了有限元进行了正演模拟，就多个模型实例进行数值模拟的结果分析总结。由于考虑到数据量的问题，因此在本章模型计算后成正演曲线时，取了其中的测深供电点中的7个点成图，其位置分别为：A 1＝－50 m、A 2＝－30 m、A 3＝－10 m、A 4＝0 m、A 5＝10 m、A 6＝30 m、A 7＝50 m。

2．1 二层模型

模型为分为两层（图2），第一层的电阻率为50 Ω·m，极化率为η＝5．0%，厚度为15 m；第二层的电阻率为100Ω·m，极化率为η＝1．0%。

图2 二层层状模型示意图Fig．2 Diagram of two layered model

从双层模型正演测深点的曲线（图3及图4），我们可以看出数值模拟较好地体现了地下有的异常体。

2．2 单正方柱体模型

模型为10 m×10 m低阻、高极化正方柱体（图5），中心位置为（0 m，－10 m），电阻率为50Ω·m，极化率为η＝5．0%；背景电阻率为100Ω·m，极化率为η＝1．0%。

图3 二层层状模型固定点源测深正演模拟各点源测深视电阻率曲线Fig．3 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in two layered model

图4 二层层状模型固定点源测深正演模拟各点源测深视极化率曲线Fig．4 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in two layered model

图5 单正方柱体模型示意图Fig．5 Diagram of single square cylinder model

从单正方柱体模型正演测深点的曲线（图6及图7），我们可以得到数值模拟较好地体现了地下有的异常体。

图6 单正方柱体模型固定点源测深正演模拟各点源测深视电阻率曲线Fig．6 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in single square cylinder model

图7 单正方柱体模型固定点源测深正演模拟各点源测深视极化率曲线Fig．7 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in single square cylinder model

2．3 正地形模型

2．3．1 含异常体

模型为一个正地形模型（图8），起伏高度为4 m，背景电阻率为100Ω·m，极化率为1%，正方柱体电阻率为50Ω·m，极化率为5%，大小为10 m× 10 m中心位置为（0 m，－10 m）。

从模型正演测深点的曲线（图9及图10）看出，电阻率明显受到了地形的影响，表现为在地形起伏段的两端出现了较大的电阻率值，因为①当供电点靠近山脊时，山脊部分，电流密度相对小，测量得到电位值低，根据视电阻率公式可得，此时视电阻率偏低，呈低阻；②当供电点接近山峰时（或者说供电点进入山脊部分），此时电流密度在山脊部分相对较大，测量得到的电位值偏高，视电阻率值变高，呈高阻。但是正地形模型正演测点的极化率曲线未受到影响，这是因为极化率未受到地形影响。

图8 正地形模型示意图Fig．8 Diagram of positive landform model

图9 正地形模型固定点源测深正演模拟各点源测深视电阻率曲线Fig．9 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

图10 正地形模型固定点源测深正演模拟各点源测深视极化率曲线Fig．10 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

2．3．2 纯地形

依据图8中的凸地形，进行纯地形响应研究。

模型正演测深点曲线得出（图11及图12），电阻率受到了地形的影响，山脊会造成一个低阻的假异常，然而极化率却不会受到影响。

图11 正地形模型固定点源测深正演模拟各点源测深视电阻率曲线Fig．11 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

图12 正地形模型固定点源测深正演模拟各点源测深视极化率曲线Fig．12 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

2．4 负地形模型

2．4．1 含异常体

模型为一个负地形模型（图13），背景电阻率为100Ω·m，极化率为1%，正方柱体电阻率为50Ω ·m，极化率为5%，大小为10 m×10 m中心位置为（0 m，－10 m）。

图13 负地形模型示意图Fig．13 Diagram of negative landform model

从模型正演测深点曲线得到（图14及图15），电阻率也受到了地形的影响，表现为地形起伏段的两端出现了较低的电阻率值，因为①当供电点靠近山谷时，山谷部分，电流密度相对大，测量得到电位值高，根据视电阻率公式可得，此时视电阻率偏高，呈高阻；②供电点进入山谷部分，此时电流密度相对山谷部分较小，测量得到的电位值偏低，视电阻率值变低，呈低阻。但是负地形模型正演测点的极化率曲线未受到影响。

图14 负地形模型固定点源测深正演模拟各点源测深视电阻率曲线Fig．14 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in negative landform model

图15 负地形模型固定点源测深正演模拟各点源测深视极化率曲线Fig．15 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in negative landform model

2．4．2 纯地形

依据图13中的负地形模型，进行纯地形响应研究。

从模型正演测深点曲线得到（图16及图17），电阻率受到了地形的影响，山谷会造成造成一个高阻的假异常，然而极化率缺不会受到影响。

图16 正地形模型固定点源测深正演模拟各点源测深视电阻率曲线Fig．16 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

图17 正地形模型固定点源测深正演模拟各点源测深视极化率曲线Fig．17 Apparent resistivity curve about simulation of the numerical relation of sounding point and fixed point source in positive landform model

4 结论

本文针对固定点源测深装置，设计了不同模型进行数值模拟，通过正演数值模拟，得到新的认识：

1）从双层模型、单正方柱体模型正演测深点的曲线，我们可以看出数值模拟较好地体现了地下有的异常体。

2）从正地形的模型（含异常体）正演测深点的曲线得到，电阻率明显受到了地形的影响，表现为在地形起伏段的两端出现了较大的电阻率值；负地形的模型（含异常体）正演测深点曲线，电阻率也受到了地形的影响，表现为地形起伏段的两端出现了较低的电阻率值。但正、负地形模型正演测点的极化率曲线未受到影响，这是因为极化率未受到地形影响。

3）从纯地形的模型正演测深点的曲线得到，电阻率也明显受到了地形的影响，表现为正地形会产生一个假的低阻异常，负地形会出现假的高阻异常；极化率未受到影响。

［1］ SMITH J T，BOOKER J R．Rapid Inversion of Twoand Three－Dimensional Magnetotelluric data［J］．Geophys Res．，1991，96：3905－3922．

［2］ 毛先进．边界积分方程和以其为基础的电阻率层析成像［D］．长沙：中南工业大学，1997．

［3］ 熊彬，阮百尧．电位双二次变化二维地电断面电阻率测深有限元数值模拟IJ］．地球物理学报，2002，45：285－294．

［4］ 罗延钟，孟水良．关于用有限单元法计算二维构造点电源场的几个问题［J］，地球物理学报，1986，29（6）：37－45．

［5］ 强建科，罗延钟，熊彬．固定点源测深激电异常研究［J］．地球物理学进展．2005，20（4）：1176－1183．

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

［7］ 周熙襄，钟本善．点源二维电法正演的有限元单元法［J］．物探化探计算技术，1983，5（3）：19－40．

［8］ 阮百尧，徐世浙．电导率分块线性变化二维地电断面电阻率测深有限元数值模拟［J］．地球科学，1998，23（3）：303－307．

Research about two－dimensional IP numerical simulation

CHEN Yong-lin1，JIANG Shou-jing2，XIE Dan1
（1．The 12th corps of Gold troops of CAPF，Chengdu 610059，China；2．The 10th corps of Gold troops of CAPF，Kunming 65000，China）

In this paper，starting from the question of point source and dimensional geoelectric field，we use finite element method to simulate geoelectric field，triangle subdivision algorithm to rolling surface，and various models to complete forward simulation according to the characteristic of bilateral three-pole device．By means of forward simulation and Inversion of various models，we have summarized some features about abnormity to offer some useful information for analysis of Induced polarization．

induced polarization method；finite element method；bilateral three-pole device

P 631．3＋24

A

10．3969／j．issn．1001-1749．2014．05．05

1001-1749（2014）05-0541-06

2014-02-25 改回日期：2014-07-23

陈永凌（1988－），男，硕士，主要研究方向为矿产地球物理，E-mail：476811938＠qq．com