0 引言

近年来，国内在直流电法正演数值模拟方面，刘树才等(2004)给出了矿井直流电法三维数值模拟存在的问题.强建科(2006)解决了存在地形起伏条件下的三维电阻率模拟及成像问题.(赵斌，2007)利用ANSYS软件对起伏地形的激发极化法进行了数值模拟.张继锋等(2013)基于起伏地形的电偶源进行了2.5维数值模拟.刘申芬等(2009)利用有限元法对井地电阻率进行了数值模拟和反演.从以上研究可以看出，常规的直流电法理论及技术发展已经十分完善，当前直流电法数值模拟主要解决矿井直流电法、地形起伏、井-地电阻率等特殊领域中的特殊问题，因此，寻求一种好的数值模拟方法及模拟软件十分重要.ANSYS软件就是一款基于有限元法的专业数值分析软件，它是融结构、热、流体、电磁学和声学于一体的大型CAE通用有限元分析软件，具有多种分析功能，包括简单线性静态分析和复杂非线性动态分析，可以用来求解结构、流体、电力、电磁场及碰撞等问题的解，包含了预处理、解题程序以及后处理和优化等模块，将有限元分析、计算机图形和优化技术等相结合，已成为解决现代工程学问题必不可少的有力工具(田永军，2007)，被广泛应用于航空航天、电子、通信、石油、化工等众多行业.

ANSYS软件在直流电法的数值模拟中有广泛的应用.汤井田等(2006)以及汤井田和王飞燕(2008)从点电源场理论出发，借助ANSYS软件通过正演模拟解与理论解的对比分析，验证了ANSYS正演模拟的正确性.高炳魁(2008)利用ANSYS软件，建立几种理论边坡模型，并对模型进行了电测深数值计算，得到边坡隐患的电位分布规律，结合物理模拟进行了对比和验证.李建华(2009)使用ANSYS软件，对井间的低阻异常体进行了数值模拟，取得不错的效果.马炳镇(2011)及张军等(2012)利用ANSYS软件对全空间下的巷道影响等问题进行了研究，分析出：巷道空腔改变稳定电流场分布，尤其在巷道迎头附近，影响极为严重.孙瑛(2012)将ANSYS软件运用到工程地质隐患探测中，模拟了工程地质隐患在直流电阻率法中产生的异常特征，取得了较为理想的效果.郭清石(2013)应用ANSYS软件，在均一地层中建立含地下溶洞槽模型，分析勘探中装置、电极间距选取问题；分析了测线偏移影响及对称四极分辨率问题.

ANSYS作为数值模拟软件，从最初模拟静电场分布，验证数值模拟的正确性，到后来借助模型分析装置类型、电极间距、巷道超前探测正演响应中，说明该软件强大的模拟功能和精确快速的计算能力，为正演计算带来了许多便利.同时，ANSYS软件的APDL语言可以灵活地实现有限元分析的众多功能，是进行ANSYS二次开发强有力的工具之一，利用该软件不仅可以实现半空间静电场的模拟，还可以实现全空间条件下异常区的响应，其众多优越性表明它在地球物理中有着广阔的应用前景.本文在矿井直流电法在全空间环境下模拟顶、底板低阻异常体对电阻率曲线的变化情况，分别建立了半空间底板异常体模型、全空间底板异常体模型、全空间顶板异常体模型、全空间顶、底板都存在异常体模型，通过电阻率剖面曲线变化，分析低阻异常体对电阻率曲线的影响.

1 ANSYS思路和过程

ANSYS思路包括以下过程：

(1)给出地球物理边值问题中的偏微分方程、边界条件及初始条件.只有对地球物理方法的原理和问题有深入的理解，才能给边值问题中的偏微分方程和边界条件以正确的描述.

(2)将地球物理问题转变为有限单元法方程.即利用变分法将微分方程的求解转化为相应的泛函的极值问题.

(3)用有限元法求解该泛函的近似解.主要涉及以下几个步骤：首先，差值，即将研究区剖分成有限个小单元，在每个单元上，将函数简化成线性函数、二次函数或高次函数.其次，计算各单元的泛函，各单元间通过单元间的节点上的函数值相互联系起来，并对各单元的泛函求和，获得整个区域上的泛函.再次，根据泛函取极值的条件，得到各节点函数值应满足的线性代数方程组.最后，求解方程组，得各节点函数值.有限单元法求解流程见图 1.

2 三维建模

利用ANSYS软件的APDLP平台进行正演模拟的流程如(图 2)所示，首先要进行三维建模.利用APDLP平台编写程序语句建立模型，程序语句简练、易懂能实现各种复杂的模型.本文以建立一个长为300 m宽为200 m高为80的半空间模型(图 3)为例子说明一下建模过程.设长轴坐标为为x轴，宽为z轴，高为y轴.源代码如下所示：

/clear！清理命令窗内容

/prep7!建立模型名称

k, 1, 0, 0, 100

k, 2, 300, 0, 100

k, 3, 300, -80, 100

k, 4, 0, -80, 100

k, 5, 0, 0, -100

k, 6, 300, 0, -100

k, 7, 300, -80, -100

k, 8, 0, -80, -100！建立8个点

v, 1, 2, 3, 4, 5, 6, 7, 8！将8个点连成体

vovlap, all！将所有的体连成一块

numcmp, all！对体进行重新编号

save！保存文件

模型建立完成后，首先对模型赋属性，赋属性的目的是：使模型具有大地电阻率特性；其次，对模型进行网格划分，网格划分可分为自动划分和人为划分，一般情况下选择四面体对模型划分，但在划分的过程中，可能出现划分不规则点的情况，需要人为调节网格单元大小；再次对模型赋边界条件，注意：不要在采集数据位置赋上边界条件，否则程序无法运行；然后对模型赋供电电流，为了简化运算，一般电流设为1 A；最后，求解电阻率值.数据采集完成后将数据导出，根据需求进行成图分析.

在半空间条件下，建立长300 m，宽为200 m，高度为80 m，围岩电阻率为200的模型，将电阻率为1半径为5 m的球体放在x=75，y=0，z=-25 m处，该模型称为半空间异常模型，记作模型1.在模型1基础上，加一长度为300 m，宽为200 m，高为80 m的顶板围岩，异常体的位置保持不变，该模型称为全空间底板异常体模型，记作模型2.在模型2的基础上，将底板异常体位置移到x=75，y=0，z=25 m坐标处，大小形态保持不变，该模型称为全空间下顶板异常体模型，记作模型3.在模型3的基础上，在x=75 m，y=0，z=-25位置加入一个跟顶相同的异常体，形成顶、底板对称异常体模型，该模型称为全空间条件下顶、底板异常体模型，记作模型4.

以模型1为例，说明在ANSYS平台下模型电阻率的求解过程.本模拟数据采集点位于y=0，z=0的x轴上采集，共使用30个电极，电极间距为10 m，采用MNB三极装置，剖面数为14，滚动数为14，每个模型共采集196个数据点.数据采集完成后，将四个模型相同地层的电阻率剖面变化图(图 4~图 9)，根据电阻率剖面的变化，分析各种模型低阻异常体响应.

3 电阻率曲线分析 3.1 全空间与半空间变化分析

首先分析全空间与半空间电阻率变化特征.通过图 4~图 9可以看出半空间和全空间的相同点：同一水平的电阻率都表现出“V”形态，在异常体处出现电阻率最低值，在异常体两侧出现极大值.原因：低阻异常体对电流的“吸引”，导致异常体两侧电流密度增大，两侧出现电阻率极大值.在模型地层浅部，电阻率曲线变化明显，深部变化不明显.

不同点：①半空间环境下相邻电极间电阻率曲线波动小，变化平缓，而全空间环境下电阻率曲线变化波动大.②异常体处最低点的变化不同.在浅层形态情况：异常体上方电阻率最低值从低到高排列：半空间异常体，全空间底板异常体，顶板异常体，顶、低板异常体.深层情况：全空间环境下异常体电阻率变化曲线趋向同一形态，但顶、底板异常体，电阻率曲线“V”行开口更宽.全空间顶板异常体和底板异常体的电阻率曲线变化形态相同.

从上述可知全空间环境使异常体对电流密度的吸引减弱，最小值点减小，同时引起相邻电极电阻率的波动和震荡，电阻率曲线变化更剧烈.

3.2 异常体位置对电阻率曲线的影响分析

本节析全空间条件下异常体位置对电阻率曲线变化的影响，因此模型1不做分析.从图 4~图 6中可以看出，电阻率曲线最低点从大到小排列依次：模型4、模型3、模型2.原因分析：模型4中顶板异常体向上“吸引”作用和底板向下“吸引”作用，与模型2相比，模型4顶板异常体对电流向上吸引作用，导致向下吸引作用“减弱”，因此，最低点大于模型2；模型3中顶板异常体向上“吸引”电流，因此，最低值大于其他情况.图 7~图 9可以看出，电阻率最低点虽然相等，但低阻体两侧电阻率值不同，从大到小排列依次是模型3、模型2、模型4.模型4中顶、底板异常体对电流向下“吸引”表现“加强”作用，导致异常体周围电流密度减少，电阻率值变小.而模型2和3随着距离的增加，顶板异常体向上“吸引”电流，减弱了向下“吸引”作用，因此模型3的值大于模型2的值，但随着深度的增加，电阻率值趋向围岩电阻率.

由上可知，模型3异常体对电阻率曲线影响主要集中在底板地层的浅层，对于地层的深层影响不大；模型3异常体使地层电阻率值减小，电阻率整体减少幅度不大，但相邻电极间电阻率值波动大.模型4异常体以底板异常体位置为标准，异常体上方地层，对电流吸引表现为“减弱作用”，使电阻率值变大；底板异常体下方地层表现为“加强作用”，使电阻值减小.模型2异常体对电阻率曲线影响深度大、范围广.

综上所述，根据异常体对电阻率曲线影响变化，可以确定异常体的位置，并根据探测目标的需求，剔除异常体的干扰情况.否则异常体将引起目标区的假异常，影响探测精度.

4 工程应用

某矿31506工作面位为15煤层3采区运输上山西侧，15煤层位于石炭~二叠系太原组下部，下伏矿井主要充水含水层为徐灰、草灰、奥灰.三采区范围内15煤层距离奥灰59.5~67.2 m，平均61.85 m.井田范围内-195 m水平以上动水补给循环条件较好、溶洞裂隙发育，属富水性强的非均质岩溶裂隙承压含水层，强富水部位主要位于断层破碎带附近，强富水层段自顶界面向下70 m以内.随着埋藏深度加大，其动水补给条件逐渐变差，岩溶裂隙发育状况及地下水循环条件逐渐减弱，但由于奥灰岩溶发育在横向上具有不均匀性，15煤层开采受奥灰水威胁.为查清31506工作面运煤巷、轨道巷底板徐灰、奥灰富水性情况，指导31506工作面开采防治水工作.

本次井下高密度二维电法勘探数据采集使用的是WJDJ-3型高密度电阻率系统，采用三极装置的采集方法，其采集参数如下：开设道数60道，测量层数29，测量滚动道距10 m采用三极装置进行数据采集.数据采集完成后，没有进行顶板异常干扰处理，经过反演得到图 10，剔除顶板异常体后得到图 11.

对比两图可以看出，两图共同点就是在地层深处的徐、奥灰对顶板地层的影响明显，不同的是图 10在地层浅部存在很多低阻异常区，而图 11浅部地层很少存在低阻异常区.最明显的区别在于图 10中存在x=65 m、x=165 m、x=295 m三处异常区，图 11在这三处存在异常区.原因：图 10没有对顶板异常区干扰造成进行剔除，导致假异常出现.在这三处“假异常”区中，最重要的是x=295 m处，如果按照图 10解释，那么x=295 m就是危险地带，容易将徐、奥灰含水层导通，容易发生突水事故，但经过剔除干扰后，“假异常”区消除.为了验证解释的准确性，对x=105 m和x=215 m和x=295 m三处进行打钻验证，经验证在x=105、x=215 m处钻孔深20 m钻孔涌水量分别是43 m³/h，35 m³/h，说明以上两处富水强，应该进行注浆加固工作面底板.而在x=295处钻孔深39 m时涌水量为3 m³/h，说明该处地层不富水，从而验证了图 11正确性.

5 结论

(1)利用ANSYS软件实现了不同环境下的直流电法数值模拟，尤其实现了全空间条件下的直流电法模拟.

(2)半空间环境下异常体对电阻率曲线影响更明显，全空间环境下顶板异常体主要影响底板浅层电阻率曲线的变化.

(3)经工程验证，基于电阻率曲线变化特征分析，成功识别出旁侧干扰，为准确探测目标区富水性打下基础.

致谢 由衷地感谢审稿专家和编辑部老师的宝贵建议和意见！