1、第 卷 第期 年月地球物理学报 ，熊治涛，唐新功，李丹丹 基于非结构网格有限元的三维井地电阻率法任意各向异性正演模拟地球物理学报，（）：，：，.（），（）：，：基于非结构网格有限元的三维井地电阻率法任意各向异性正演模拟熊治涛，唐新功，李丹丹油气资源与勘探技术教育部重点实验室（长江大学），武汉 南方科技大学地球与空间科学系，深圳 武汉大学中国南极测绘研究中心，武汉 摘要地层介质的电各向异性增加了井地电阻率法响应的复杂性，开展基于电各向异性介质模型的井地电阻率法响应规律研究对于正确解释各向异性显著地区的观测数据特征至关重要 针对垂直线源井地电阻率法的任意各向异性响应模拟问题，本文提出了一种基于非结
2、构网格有限元三维正演算法，通过引入的对称正定张量来表征任意各向异性的电导率，采用非结构四面体网格有限元方法来离散电位的边值问题，通过将垂直线源等效为一系列点源问题，进而实现了任意各向异性介质中井地电阻率法的高效数值计算 通过与三个地电模型的解的对比，验证了本文数值解算法的精度和有效性 针对线源远离和垂直穿过异常体的两类模型，分别考察了当围岩或异常体为各向异性介质时的井地视电阻率响应特征 结果表明，对于各向异性地层，围岩和异常体的主轴电阻率值和旋转角均会对井地视电阻率的幅值及分布产生显著的影响 研究结果对于提高井地电阻率法的认识和资料解释水平具有重要的理论和实际意义关键词井地电阻率法；各向异性；
3、有限元；正演 ：中图分类号 收稿日期 ，收修定稿基金项目国家自然科学基金项目（，）资助第一作者简介熊治涛，男，年生，博士后，主要研究方向为电磁法勘探 ：通讯作者唐新功，男，年生，博士，教授，博士生导师，主要研究方向为电磁法勘探与地球动力学 ：，（）地 球 物 理 学 报（）卷 ，；引言电磁勘探方法作为勘探地球物理领域的重要分支，是以地下介质的电性差异为物性基础，来实现对地下地电结构的推断（底青云等，；，）不同于常规的地面电阻率方法，井地电阻率法是通过套管或将激励源放入井中向地层供入大功率电流，在地表布设测网观测由地下介质的电性变化而引起的电位异常，从而达到对地电结构高精度成像的目标 该方法的场
4、源更靠近勘探目标体，使得在横向探测范围和纵向分辨率方面均较常规的地面电阻率方法更具优势，因此被广泛地应用于固体矿产、环境工程和水文勘查以及油气开发的油田注水和压裂效果动态监测等领域（，；，；，；，；，）利用井地电阻率法来精细刻画地下复杂的三维地电结构时将会产生庞大的三维数据体，为了对这些数据体进行快速准确的反演解释，开发适应于复杂地电模型的高效与高精度的井地电阻率法三维数值模拟算法势在必行 国内外一些学者针对井地电阻率法问题开展了一系列先导性研究 和 （）基于有限差分技术实现了三维井地电阻率法正演模拟；（）利用神经网络技术实现了对日本 地区三维井地数据的快速反演；李长伟等（）通过有限元方法完成
5、了三维带地形模型的井地电位模拟，考察了点源埋深对地表异常电位分布的影响；戴前伟等（）研究了复杂线源的井地电位响应，表明线源的分布形态严重影响到地表异常电位的分布特征；等（）利用有限差分法模拟了三维井地和井间电位响应，验证了井地与井间联合应用对剩余油圈定的有效性；等（）、杨沁润等（）、等（）讨论了不同类型电流源的井地电位响应特征，并利用阻尼最小二乘法实现了对井地电位数据的反演计算；等（）利用差异场技术提高了数值解的精度，并详细研究了井地电阻率法探测效果的影响因素；孟麟等（）通过有限元方法模拟了山谷地形下 维井地电阻率法的响应；王智等（）采用非结构化网格有限元方法模拟了三维井地电阻率法的响应，考察
6、了地形对低阻异常响应的影响以上关于井地电阻率法的研究均是基于电阻率各向同性理论的假设，而电各向异性现象在地球内部是普遍存在的（，）随着电磁勘探技术的不断发展，地层介质呈现出的电各向异性现象逐渐得到了国内外学者的广泛关注 因此基于电各向异性理论的直流电法（，；，；任政勇等，；朱姣等，）、频率域电磁法（，；，；殷长春等，；，；，；邱长凯等，；，）和时间域电磁法（，；，；周建美等，；刘亚军等，）均开展了大量研究，都证明了各向异性对电磁场响应存在明显影响并会显著增加数据处理与解释的难度期熊治涛等：基于非结构网格有限元的三维井地电阻率法任意各向异性正演模拟尽管基于电各向异性理论的各种电磁方法的正反演研究
7、一直是国内外学者关注的焦点，但是目前针对任意各向异性地电模型的三维井地电阻率法的电位响应的计算却还未见报道 另外，传统井地电阻率法的三维数值模拟多是采用结构化网格来离散模型，对起伏地形和复杂地质结构的刻画会带来困难，不利于解决复杂地电模型问题；而相比之下，非结构化网格不仅能够高精度地拟合复杂模型的边界变化，而且同时具有较小的自由度 因此，本文拟采用非结构网格有限元方法开展任意各向异性介质中三维井地电阻率法的响应规律研究，详细考察并分析当异常体和围岩分别为各向异性介质时井地电阻率法的视电阻率响应规律 这对于深化认识与开展电各向异性模型的井地电阻率法的研究，提高井地电阻率法的勘探效果和应用水平都具
8、有重要的意义井地电位响应模拟算法地下稳恒点电流源的三维边值问题可表示为：（）（），（）式中，为电位；为电流强度；（，）为点源位置；为 函数；为地空界面（）的外法线；表示无穷远处的侧边界和底边界；为电导率，对于任意各向异性介质，是的对称、正定张量，一般可由一个主电导率张量经过三重欧拉旋转得到（，）：（）（）（）（）（）（），（）式中，（，）表示张量电导率的各个元素；、和为主电导率张量分别围绕电性主轴、轴旋转的角度（电性主轴坐标系 与观测坐标系 对应）；、和分别为坐标旋转矩阵，其表达式为：（），（），（）.依据变分原理，式（）对应的变分问题为：（）（）（）（）.（）井地电阻率法利用套管向地层供电时
9、，套管直径远小于套管长度和测点与套管间的距离，因而可将其简化为有限长线电流源模型来处理（，）有限长线源产生的电位，可以看成是将线源划分为段小线源产生的电位的叠加，而每段小线源均可当作点源处理 故线源激励的变分问题可表达为：（）（）（）（）.（）将式（）代入到式（）中，得到了如下表达形式：（）（）（）（）（）.（）本文通过四面体单元剖分地电模型，能够保证在较小计算量的前提下，方便的对场源和测点处进行加密剖分（如图所示）在线性四面体单元内的任意点（，）的电位可表示为：，（）式中和分别为四面体单元第个节点的插值函数和电位值地 球 物 理 学 报（）卷图场源和测点处四面体网格示意图 将式（）代入到式（
10、）中，可方便地进行单元分析并求得刚度矩阵 最终求解的大型线性方程组为：，（）式中，为大型稀疏系数矩阵，为待求节点电位，为源项表达式 最后采用 求解器求解式（）视电阻率计算公式为：，（）式中，为电位差，、为测量电极，装置系数的表达式为：（）（）（）（）（）（）（）（）（），（）式中，为或，和为发射线源的顶和底埋深，为测点埋深，为测点距源的径向距离数值算例 算法精度验证模型一：均匀半空间模型的电阻率为 ，垂直线源长度为 ，电流强度为（下同），回流极位于无穷远处（下同）由图可见，均匀半空间模型的数值解与解析解基本重合，最大相对误差仅约为 图 均匀半空间模型的数值解与解析解以及相对误差曲线图 模型二：
11、围岩电阻率为 ，低阻球体的电阻率为、半径为、球心位置为（，），点源位于坐标原点测线沿着轴，测点间距为 由图可见，低阻球体模型的数值解与解析解基本重合，最大相对误差约为 图低阻球体模型的数值解与解析解以及相对误差曲线图 模型三：为了进一步验证算法的有效性，设计了围岩电阻率为 、，异常体的规模为 、中心位置为（，）、电阻率为、的三维各向异性模型 点源位于坐标原点，采用二级装置观测，个接收点均匀分布在收发距为 的圆周上 将本文的有限元法计算结果与谱元法（）计算结果（朱姣等，）进行比对 由图可见，有限元法和谱元法的计算结果基本重合，最大相对误差约为 以上检验结果表明了本文算法的有效性和高精度 井旁异常
12、体模型的井地电阻率法响应为了研究异常体和围岩分别为各向异性介质时井地电阻率法的响应规律，设计了如图所示的模型四，其中垂直线源长度为 ，异常体的中心埋深 、中心距源 、规模为 期熊治涛等：基于非结构网格有限元的三维井地电阻率法任意各向异性正演模拟图三维各向异性模型不同算法的视电阻率（）和相对误差（）曲线图 （）（）图模型剖面（）和平面（）以及测点分布（）示意图 （）（）（）（按，顺序）测线布设方式如图 所示，以发射源为圆心，按环形测网方式呈放射状布置测线，测线近源端距源，测点间距为，仅观测测点的径向电位差（下同）围岩为各向同性、异常体为各向异性时井地电阻率法的响应特征下面考察围岩为各向同性、异常
13、体为任意各向异性模型的井地电阻率法响应 其中围岩电阻率为 ，各向异性异常体的电性参数设置见表所示 计算中除了改变模型四中异常体的电性参数外，其余参数均保持不变 表中 分别代表异常体的主电导率张量围绕电性主轴、轴旋转（如图所示）时电各向异性参数设置情况图中给出了异常体电阻率分别为 和 时各向同性模型的视电阻率响应 由于低阻体对电流的汇聚作用，导致在低阻体上方观测的电位差变小，使得视电阻率在低阻体上方出现极小值，且低阻体的电阻率越小，视电阻率值就越低 图中视表异常体的电性参数 （）（）（）电阻率异常区清晰地指示了低阻体的水平位置图中横、纵坐标表示测点坐标值（单位为，下同）（）异常体主电导率张量围绕
14、各向异性主轴轴旋转（）图给出了异常体主电导率张量围绕轴旋转时的视电阻率响应 从图中可以看出，随着旋转角的增大，视电阻率值基本不变，且与图中异常体电阻率为 的各向同性模型结果基本一致 这是 由于随着角增大，异常体在方向的电阻率始地 球 物 理 学 报（）卷图主电导率张量旋转示意图（）绕轴旋转；（）绕轴旋转；（）绕轴旋转 （）；（）；（）图各向同性模型的视电阻率响应（）；（）图中黑色线框代表异常体的水平位置，下同 （）；（），终为 且保持不变，根据异常体与源的相对位置以及按放射状布设的测线，使得异常体上方测点的电位差主要由模型轴方向的电性结构决定 因此，本模型的视电阻率响应对异常体主电导率张量围绕
15、轴旋转不敏感（）异常体主电导率张量围绕各向异性主轴轴旋转（）图为异常体主电导率张量围绕轴旋转时的视电阻率响应从图中可以看出，随着旋转角的不断增大，低阻体上方的视电阻率值逐渐减小 当 时，视电阻率的幅值和分布规律与异常体电阻率为 的各向同性模型结果（图）相近；当 时，视电阻率响应与异常体电阻率为的各向同性模型结果（图）相似这是由于随着角增大，异常体在方向的电阻率值由 减小至，进而模型的视电阻率响应由接近于异常体电阻率为 的各向同性模型的结果（图）逐渐向异常体电阻率为 的各向同性模型的结果（图）靠近（）异常体主电导率张量围绕各向异性主轴轴旋转（）图 为异常体主电导率张量围绕轴旋转时的视电阻率响应由
16、图 可见，异常体的各向异性对视电阻率响应产生了显著的影响 当旋转角从 增大到 时，异常体在方向的电阻率值由 减小到了，同时在异常体上方测量的结果主要展示模型轴方向的电性结构因此，视电阻率响应特征由接近于异常体电阻率为 的各期熊治涛等：基于非结构网格有限元的三维井地电阻率法任意各向异性正演模拟图异常体主电导率张量围绕轴旋转时视电阻率响应 图异常体主电导率张量围绕轴旋转时视电阻率响应 地 球 物 理 学 报（）卷图 异常体主电导率张量围绕轴旋转时视电阻率响应 向同性模型的结果（图）逐渐向异常体电阻率为 的各向同性模型的结果（图）靠近 异常体为各向同性、围岩为各向异性时井地电阻率法的响应特征在层理发育地区，围岩也会呈现出明显的电各向异性现象，因此研究各向异性围岩的井地电阻率法响应特征对于提高各向异性显著地区实测资料处理与解释的精度和可信度具有积极的意义 下面考察模型四（图）中任意各向异性围岩的井地电阻率法响 应 规 律，其 中 各 向 同 性 异 常 体 的 电 阻 率 为，围岩电性参数设置见表所示 计算中除了改变围岩的电性参数外，其余参数均保持不变 表中 分别代表围岩的主电导率张量围绕电性
