赵荣春,吕玉增,韦柳椰,文达麟,刘 波

(桂林理工大学 a.地球科学学院; b.有色及贵金属隐伏矿床勘查教育部工程研究中心, 广西 桂林 541006)

随着我国经济的快速发展,对煤炭、金属矿等资源的需求也进一步加大。然而当前浅地表矿产资源面临枯竭,资源勘探已逐步向立体、纵深方向发展,大力发展深部以及外围找矿已经变得十分迫切,加强深部找矿方法的研究将具有十分重要的意义[1-2]。传统的地面激发极化法勘探深度有限,很难适应大深度找矿的需求,很大程度上影响了勘探的效果。井中激电法作为地面激电法的一种延伸,通过将供电点或测量点放置于地下,提高了激电探测的深度,同时还可以减小地表覆盖层以及人文干扰对测量的影响,充分依靠已知的井巷资源,可在探测井旁矿体、追索矿化带、井底盲矿、估算见矿深度等方面起着重要作用[3-6]。

针对地面、井中等三维激电数值模拟的问题,国内外学者进行了诸多研究并取得了长足的进展,黄俊革等[7]采用径向剖分的有限元法对井地三维电阻率法进行数值模拟,王智[2]、李长伟[8]采用非结构化网格剖分的方式对井中激电实现了数值计算,吕玉增等[9]利用有限单元法实现了地-井五方位激电观测,Wang等[10]、Avdeev等[11]分别用有限差分和积分方程法对电阻率测井进行了三维模拟,Rücker等[12-13]采用非结构化网格实现三维电阻率有限元正反演,提高了对复杂模型的处理能力。这些工作对于井中电法勘探具有很好的指导意义,但是常用的基于四面体下的结构化网格剖分方式,剖分单元在3个方向上长度不能同时放大,井眼环境与实际形状也存在差异;非结构化网格剖分可灵活地对网格进行局部加密,可处理任意形状的模型,然而非结构化网格的数据结构和算法复杂,网格节点难以控制。

本文基于柱坐标下以井眼为中心进行结构化的放射状网格剖分,使得井眼环境与实际形状一致,可减少网格方向性影响引起的误差,同时能在径向和深度方向上进行灵活的网格加密,既能有效减少网格单元,又能获得较好的网格质量,确保计算精度。针对地-井观测、单井井中观测两种方式, 讨论了不同形态异常体在井中观测时激电异常特征,探讨了异常体空间位置及钻孔环境的变化对激电观测的影响规律。

1 网格剖分及有限单元法

1.1 放射状网格剖分

放射状测网网格剖分基于柱坐标进行, 利用柱坐标和直角坐标系之间存在简单关系:x=rcosθ,y=rsinθ,z=z, 在矩形规则测网网格剖分的基础上, 把实际的观测面投影到R-θ面上进行网格剖分和单元编号,再利用直角坐标与柱坐标之间的相对关系再返回实际平面,完成对放射状测网的网格剖分,如图1、图2所示。

图1 放射状测网网格剖分平面投影示意图Fig.1 Planar projection diagrams of radial grid dissection

图2 放射状测网网格剖分立体投影示意图Fig.2 3D projection diagrams of radial grid dissection

放射状测网网格的剖分以柱坐标中轴线为中心进行放射状剖分, 中心区域为井眼所在位置,采用三棱柱单元进行剖分, 之后将每个三棱柱单元拆分为3个四面体单元(图3), 外围区域剖分为不规则的六面体单元, 每个不规则六面体单元拆分为5个四面体单元(图4), 因此放射状剖分仍然基于四面体进行, 网格剖分的间距由圆中心向外围逐步增大, 边界取足够大的区域。 以井眼为中心进行放射状的网格剖分, 使得网格大小以井眼为中心呈放射性逐步增大, 同样的剖分区域占用的网格数量大幅减少, 同时, 由于区域剖分单元以井口为中心环形分布, 单元的一个边与电场等位线基本重合, 这与井中观测电场分布的实际情况更加吻合,提高了单元节点之间的电位插值精度。

图3 三棱柱单元拆分示意图Fig.3 Mitsubishi column unit split schematic diagrams

图4 不规则六面体拆分示意图Fig.4 Irregular hexahedron split diagrams

1.2 有限单元法

针对四面体单元,可用单元上4个角点的异常电位值来表示异常电位的连续函数。对于任意四面体单元e, 其4个角点编号和对应的异常电位值可分别用1、2、3、4和u1、u2、u3、u4来表示。因此,e中任意点p的异常电位u,都可用u1、u2、u3、u4进行线性插值来近似表达,其中Ni表示形函数,它是由x、y、z构成的线性函数[3]:

(1)

为了提高计算精度,采用三维电场的异常电位法求解,三维地电断面点电源异常电位的边值问题可归纳为[14-17]

(2)

(3)

其中:K和K′为系数矩阵；u为待求异常电位向量；u0为正常电位向量。解方程组(3),可得到巷道各节点的异常电位,异常电位u和正常电位u0相加得到总电位V,进而计算视电阻率。按照Seigel理论,等效电阻率ρ*与真电阻率ρ之间的关系为[22]

ρ*=ρ/(1-η)。

(4)

视极化率的计算公式为

(5)

其中: ΔUZ为观测总电位差; ΔU1为无激电效应的一次电位差; ΔU2为二次电位差。

2 数值算例与结果分析

2.1 计算精度分析

为验证计算程序的正确性,选用低阻球体模型,分别采用三维有限元程序与异常球体解析式计算[19],并进行对比验证。使用三极测量装置(图5)在地面A处进行供电,井中观测,模型参数设置如下:ρ0=100 Ωm,ρ1=10 Ωm,r=1 m,d=2 m,h=10 m,L=10 m,H=30 m,电位观测,观测点距为1 m。计算结果如图6所示,数值计算结果与解析解基本吻合。

2.2 不同形态异常体激电观测

为探讨地-井、 单井井中观测方式下井旁不同形态异常体所引起的激电异常的特征, 在不考虑钻井环境的情况下, 设计钻井深度H为500 m, 分别设置棱柱体模型、 水平板状体模型。 棱柱体中心埋深h位于地下300 m, 与井孔的直线距离d为20 m, 且位于井孔正西方向(图7), 异常体具体尺寸如图8所示; 水平板状异常体上边界距地面h=300 m, 异常体右边界与井孔的直线距离d为15 m, 位于井孔正西方向(图9), 异常体具体尺寸如图10所示; 异常体电阻率为20 Ωm, 极化率为0.5, 围岩电阻率为200 Ωm, 极化率为0。

图5 地井模型示意图Fig.5 Surface-borehole model sketch

图6 数值解与解析解对比Fig.6 Comparison between numerical and analytical solution

图7 棱柱体模型观测示意图Fig.7 Prism model observation diagram

图8 棱柱体模型示意图Fig.8 Schematic diagram of prism model

图9 水平板状体模型观测示意图Fig.9 Horizontal plate model observation diagrams

地-井观测时,设计供电点A0位于异常体相对井口的方向上(主方位)距离L为250 m、A1位于异常体相对井口的反方向上(反方位)距离L为250 m、A2位于井口处,分别进行供电测量,无穷远B极置于井孔正北方位大于5倍AO处,测量点距MN为10 m; 单井井中观测时, 从井口往下采用类似于高密度电法中的变断面连滚A-MN三极排列进行观测,分别绘制视电阻率、视极化率异常图(图11、 图12)。地-井观测时,在不同位置供电测量,都可在异常体中心位置(300 m)处观测到明显的低阻高极化异常,相比于不同方位供电情况下的激电异常特征,井口方位供电时,低阻高极化异常现象以极值点为中心呈基本对称的形态。主方位供电时,视电阻率曲线在靠近异常体中心时先呈下降趋势并达到极小值,并在异常体下方出现视电阻率极大值点,视极化率曲线则是在靠近异常体中心时先呈上升趋势并达到极大值,在异常体下方出现视极化率极小值点;反方位供电时,视电阻率、视极化率曲线形态则与主方位供电时的曲线形态刚好相反。在实际工作中,根据井口及其他方位的激电异常情况,结合其他的已知地质信息可对判断异常体的赋存方位提供一定的参考。 单井井中观测时, 以井深为纵坐标, (1/4)AO为横坐标, 以异常体为中心, 可观测到非常明显的高阻低极化异常, 且异常形态及位置与实际模型有很好的对应关系。

图10 水平板状体模型示意图Fig.10 Schematic diagram of horizontal plate model

图11 地-井观测棱柱异常体(a, b)、水平板状异常体(c, d)激电异常曲线Fig.11 IP anomaly curves of prism anomaly(a, b)and horizontal plate anomaly(c, d)for surface-borehole observation

图12 单井井中观测棱柱异常体(a, b)、 水平板状异常体(c, d)激电异常图Fig.12 IP anomaly of prism anomaly(a, b) and horizontal plate anomaly(c, d) for single borehole observation

2.3 异常体空间位置变化对观测的影响

为研究地质异常体空间位置变化对井中激电观测的影响规律, 本节采用地-井观测方式, 分别改变低阻棱柱异常体的埋深及井旁距离进行研究, 其中模型参数及供电点位置与上节一致,如图7、 图8所示。 设置以下两种试验条件： ①保持异常体、 围岩的电性参数不变, 只改变异常体埋深, 分别在距离地表h=120、 180、 240、 300 、 360 m设置棱柱异常体; ②保持异常体、 围岩的电性参数不变, 只改变异常体在水平方向上的位置, 分别在距离井孔d=15、 25、 35、 45、 55 m设置棱柱异常体。 绘制视电阻、 视极化率异常曲线如图13所示。

通过不同埋深条件下异常曲线的对比,可见不同埋深条件下所获得的激电异常曲线的形态基本一致,激电异常振幅随着异常体埋深的增加逐渐变小,这是由于测点逐渐远离供电点,导致测量时所获得的电信号逐渐减弱所致;异常体与井孔的水平距离对激电观测效果具有很大的影响,当异常体距离井孔的水平距离较近时,激电异常较为明显,随着异常体与井孔距离的增加,激电异常幅值快速减小。 综合分析可知, 异常体埋深与水平位置的变化都会对激电异常观测产生影响,在实际的野外观测中,当测点深度过大时,可适当增加供电电流,以提供更强的电信号,由于异常体井旁距离对观测效果有很大影响,这对井中激电观测的开展带来了很大的制约,在实际的生产实践中,应尽量使观测井孔与异常体靠近,以获得更好的探测效果。

图13 地-井观测异常体变埋深(a, b)、 变水平位置(c, d)激电异常曲线Fig.13 IP anomaly curves for anomaly bodies with varying buried depth(a, b) and horizontal position(c, d) for surface-borehole

2.4 钻孔环境对地-井观测的影响

钻孔是一个特殊的工作环境,井中观测需要在井中布设电极,同时需要保证电极与井壁接触良好,一般以水作为导电介质,保证电极与井壁充分接触,或者采用刷子电极,让电极附着在井壁上,但是两种情况下的钻孔环境与围岩都有着较大的电性差异,因此,了解钻孔环境对井中观测的影响是很有必要的。下面结合当前矿山验证孔、普查孔的实际孔径大小并进行适当放大,当钻孔含水或具有较高的电阻率时,分析钻孔环境对井中观测的影响规律。

设计井深500 m, 供电点A位置距离井孔中心正西方位距离L为200 m, 围岩电阻率为1 000 Ωm, 极化率为0。 ①井液电阻率为60 Ωm, 极化率为0, 分别设计钻孔直径(D)分别为0、 100、 200、 400、 800 mm; ②钻孔直径为200 mm, 钻孔(井液)电阻率分别为60、 100、 200、 2 000、 6 000 Ωm, 极化率为0。 绘制视电阻率异常曲线图14。

可见改变钻孔直径大小,对地-井观测带来的影响主要体现在井的浅部位置以及井的末端,钻孔直径越大,对观测数据的影响也越大,影响范围也逐渐变宽,尤其是对浅部数据的影响应该引起重视,井底部分的视电阻率会出现较大的跳跃,这是由于井底与围岩存在较大的电阻差异,使得靠近井底与围岩分界面时,测点处电流被高阻围岩排斥使得电流密度降低,导致视电阻率减小。当钻孔电阻率低于围岩时,观测结果受到的影响较大；当钻孔电阻率较大时,对观测的结果则相对较小,影响的范围同样集中在浅部以及井的底部；当钻孔电阻率高于围岩电阻率时,靠近井底部的测点其电流被低阻围岩吸引,使得电流密度升高导致电阻率曲线向右弯曲。

图14 钻孔环境对地-井观测的影响Fig.14 Effect of borehole environment for surface-borehole observation

3 结 论

采用基于柱坐标下的放射状网格剖分方式,对井中激电观测进行了有限元数值模拟,通过分析得出以下结论:

(1)通过对不同形态异常体进行计算,结果表明低阻高极化异常形态与实际模型具有很好的对应性,验证了本文数值计算结果的正确性,满足异常分析的要求。

(2)异常体埋深与水平位置的变化都会对激电异常观测产生明显的影响,其中异常体井旁距离的大小决定了观测的效果,给井中激电观测的开展带来了很大的制约,在实际的生产实践中,应尽量使观测井孔与异常体靠近,以获得更好的探测效果。

(3)钻孔环境带来的影响主要集中体现在井的浅部位置以及井的末端,钻孔直径越大,对观测数据的影响越明显,钻孔在低阻环境下也会对观测产生较大的影响。因此，在实际工作中,对井中观测数据的处理与解释必须考虑钻孔环境带来的影响。