基于ANSYS的复杂地形条件下电法数据地形改正研究

尚建阁,李晨晖,苏丽娟,李 冰,谢 崇,邢尚鑫

(1.河南省有色金属矿产探测工程技术研究中心,河南 郑州 450016;2.河南省有色金属地质矿产局第二地质大队,河南 郑州 450016;3.河南省自然资源科技创新中心(矿山生态环境保护修复研究),河南 郑州 450016;4.河南有色地质矿产集团有限公司,河南 郑州 450016;5.河南省有色金属地质勘查总院,河南 郑州 450052;6.河南省有色金属深部找矿勘查技术研究重点实验室,河南 郑州 450052)

0 引言

随着全球能源矿产的需求量不断增加,找矿工作的程度也越来越深入。直流激发极化法是一种常规的物探方法,也是目前应用最为广泛的物探方法之一,通过该方法可以快速的获取地球内部介质的极化率、电阻率等物理参数,但地球表面并不是理论上的一个光滑的平面,而是一个地形高低不平、连绵起伏的复杂曲面,特别是在山区,地形变化更加复杂,正是由于这些复杂地形的存在,在使用直流激发极化法开展找矿过程中,发现的异常会产生畸变,可能导致淹没地质体的有用信息,呈现假异常,使研究者不能做出正确的推断解释,降低了物探综合参数利用水平[1-3]。因此,复杂地形条件下电阻率地形改正方法的研究,成为物探工作者急需解决的关键问题之一。ANSYS软件是一款大型的融电磁学、结构力学、流体力学、仿声学等于一体的大型通用有限元分析软件[4]。ANSYS具有丰富的单元库,提供了对各种物理场量的分析功能[5-6],为地球物理工作者在进行数值模拟分析时提供了强大的数学分析功能。前人运用ANSYS开展了一系列的地球物理电法数据的数值模拟,对典型的地电模型全空间介质影响进行模拟,验证了ANSYS开展电法数据模拟的正确性,均取得了较为理想的理论效果[7-9],但是对复杂地形条件下的模拟应用较少。

文章通过介绍运用ANSYS开展电法地形改正数值模拟,总结经验并运用到实际生产找矿工作中,取得了良好应用效果,为广大地球物理工作者在复杂地形条件下开展直流激发极化法数据解释推断提供技术参考。

1 方法原理

地形影响问题是从20世纪50年代伴随着直流电法在矿产资源勘查等方面的大规模应用开始受到重视的[10]。复杂的地形对直流激发极化法中电阻率参数的影响是致命的,因此,正确消除起伏地形对电阻率参数的影响,对后期地球物理成果解释推断的真实性起到了至关重要的作用[11]。

研究地形对视电阻率的影响,目前主要使用的有三种方法:物理模拟法、解析法和数值模拟法。所谓直流电法的地形改正,就是通过上述三种方法将在起伏地形上测量得到的电阻率数据转换为平地形影响下数据的过程,从而消除不同地形所引起的误差。由于物理模拟法和解析法不适用于场源和物性分布复杂的正演情况,因此,数值模拟法则成为电阻率法正演最主要的方法[12]。随着计算机技术的发展和物探解释方法的进步,数值模拟技术得到了飞跃式发展。常用的电阻率数值模拟方法有有限差分法、积分方程法、边界单元法和有限元法,其中有限元法具有较高的求解精度且理论较为成熟,能够适应任意模型形态。因此,有限元法是目前处理电法正演问题的主流方法。

ANSYS是一款能实现多场及多场耦合分析、实现前后处理、求解及多场分析统一数据库的一体化大型有限元分析软件,一般而言,ANSYS的基本分析过程可以分为三步,即创建有限元模型、施加载荷进行求解和查看分析结果。对应软件结构的三个程序模块:前处理模块(PREP7)、分析求解模块(SOLUTION)和后处理模块(POST1和POST26)。前处理模块(PREP7)的主要任务是建立结构分析的有限单元模型,它是结构分析的开始,一般步骤为分析准备、设置单元类型、设定实常数、定义材料属性、创建模型、划分网格。分析求解模块(SOLUTION)的主要任务是进行结构的计算,一般步骤:定义自由度、施加荷载、设定分析类型、给定边界条件及求解等。后处理模块(POST1和POST26)的主要任务是进行结构的分析,它包括通用后处理器(POST1)和时程后处理器(PSOT26)两种,前者主要用于静态分析的结果后处理,而后者主要用于动态分析结果的后处理。

2 模拟实验

利用ANSYS软件进行地形改正,其一般思路:① 根据电法工作的实际情况,依照测点高程建立物理模型,并要根据计算精度,对模型进行网格剖分和边界条件的设定;② 根据工作精度的要求,合理选择单元类型,利用ANSYS程序进行纯地形条件下各测点的视电阻率值计算;③ 通过比值法,对实测视电阻率数据进行地形影响的压制。计算公式采用ρs改=(ρs实/ρs地形)×ρ0,通过ANSYS有限元分析得到的纯地形视电阻率数据,对实测数据进行比值压制,最终得到改正后的ρs改数据。总的来说,利用ANSYS软件进行地形改正需要正确记录实际地形高程数据,利用软件所提供的有限元分析软件计算出纯地形条件下影响的数值,再利用实际测量获得的电阻率数值与之进行比值法进行地形影响消除得到改正后的真实数值,基本处理流程见图1。

图1 ANSYS有限元数据分析流程图

--------------------------------------------------------------------------------------------------------------

下面利用ANSYS编写命令流加以说明。

(1)!进入前处理器,建模

/prep7

/PNUM,area,1

et,1,plane230 ! 定义单元类型

mp,1,resistivity ! 定义材料电阻率参数

*ask,KP_N,8,, ! 定义创建地形模型的点数

*dim,KP_NX,ARRAY,KP_N,,

*dim,KP_NY,ARRAY,KP_N,, ! 定义地形节点的XY数值

k,I,KP_NX(I),KP_NY(I) ! 定义关健节点

a,1,2,3,4,5,6,7,8…… ! 由点创建面

(2)进入求解器!将模型网格划分,定义边界条件和施加荷载

asel,s,,,1,2

aatt,1

asel,s,,,4

aatt,2

asel,s,,,1,4

esize,,20 ! 定义网络划分尺度

mshape,1,2d ! 定义剖分类型

amesh,all ! 开始网络划分

lsscale,all

nsel,s,ext ! 选择边界上的节点

nsel,u,loc,y

*get,Nnod,node,,count

*get,Nmin,node,,num,min

*do,J,1,Nnod,1

Dist1=SQRT((NX(Nmin)-NAX)**2+(NY(Nmin)-NAY(Nmin))**2)

Dist2=SQRT((NX(Nmin)-NBX)**2+(NY(Nmin)-NBY(Nmin))**2)

V0=Currt*Res1*(1/Dist1-1/Dist2)/(2*Pi)

D,Nmin,VOLT,V0 ! 给边界上的节点赋值

Nmin=ndnext(Nmin)

*enddo

fk,2,amps,1

fk,11,amps,-1 ! 施加电流荷载

finish

/sol

antype,static

! /status,solu ! 定义分析类型和求解

solve

FINISH

(3)! 进入后处理器,查看结果

/post

*SET,DJ,0.2

*SET,MN,0.2

*DIM,V1,ARRAY,20,1,

*DIM,Ra,TABLE,20,2,

*DO,I,1,20,1

V1(I,1)=VOLT(NODE(2+DJ*(I-1),6,0))-VOLT(NODE(2+DJ*I,6,0)) ! 求模拟电位差

AM=2+DJ*(I-1)

AN=2+DJ*I

BM=6-DJ*(I-1)

BN=6-DJ*I

KK=2*3.14156/(1/AM-1/AN-1/BM+1/BN)

Ra(I,1)=2+I*DJ

Ra(I,2)=V1(I,1)*kk/10

*enddo

FINISH

/GCOL,1,dis

/GCOL,2,Resistivity

*VPLOT,Ra(1,1),Ra(1,2),2

FINISH

根据ANSYS软件进行地形改正的思路和分析流程,分别开展了不同坡度角的山谷地形和山脊地形对直流激发极化法电阻率影响的数值模拟(图2)。图2a为山脊模型,θ为坡度角,材料电阻率为100 Ω·m,为了研究山脊纯地形对视电阻率的影响程度,分别进行了坡度为0°、10°、20°、30°、40°的模型计算。由图2a可见,在山脊纯地形上得到的视电阻率曲线呈现为低阻特征,并且随着坡度角的增大,其低阻特征更加明显。图2b为山谷模型,θ为坡度角,材料电阻率为100 Ω·m,同样也进行了坡度为0°、10°、20°、30°、40°的模型计算。由图2b可见,在山谷纯地形上得到的视电阻率曲线呈现为高阻特征,并且随着坡度角的增大,其高阻特征更加明显。由图2可见,纯地形异常主要发生在角域顶点附近,也就是说主要发生在坡度变化的地方。异常的大小与坡度角的大小有关,即与地形起伏的陡缓程度关系密切。形态相同的角域,异常幅度随坡度的增大而增大。

图2 不同地形条件影响下模拟电阻率曲线图

为了更好地研究地形对异常体电阻率曲线的影响,在数值模拟器中加入异常体进行数值模拟,图2c为山脊地形下高阻直立薄板电阻率曲线图,θ为坡度角,背景电阻率为100 Ω·m,直立高阻体电阻率为1000 Ω·m。由图2c可见,在水平地面时,直立薄板视电阻率异常表现为高阻特征,随着地形坡度的变大,异常曲线表现为低阻特征,尤其当坡度角为30°、40°时,基本为纯地形引起的异常曲线特征。图2d为不同山谷地形下低阻水平薄板电阻率曲线,θ为坡度角,背景电阻率为100 Ω·m,水平薄板的电阻率为10 Ω·m。由图2d可见,在水平地面时,水平薄板电阻率曲线表现为低阻特征,随着地形坡度的变大,异常曲线表现为高阻特征,尤其当坡度角为30°、40°时,基本为纯地形引起的异常曲线特征。由此可以看出,地形起伏对异常体引起的电阻率曲线影响很大,特别是当地形切割严重时,对电阻率资料的影响更大,使异常体电阻率异常发现较大畸变,甚至出现相反的结果。图2e为高阻直立薄板山脊地形改正前后的电阻率曲线。由图2e可见,蓝色曲线为地改前的电阻率曲线,当存在地形影响时,高阻体也表现为低阻特征;红色曲线为地改后的电阻率曲线,通过比值法改正后,电阻率曲线表现出高阻特征,与实际地下介质情况相符。图2f为高阻直立薄板山脊地形改正前后的电阻率曲线。由图2f可见,蓝色曲线为地改前的电阻率曲线,当存在地形影响时,由于地形的影响,低阻体在电阻率曲线上也表现为高阻特征;红色曲线为地改后的电阻率曲线,通过比值法改正后,电阻率曲线表现出低阻特征,与实际地断面介质情况相符。通过电流源有限元分析计算,其结果充分表明,地形对电阻率的影响很大,特别是地形起伏较大地区,影响程度更大。通过地形模拟实验可以总结得出地形对电阻率能产生较大影响,具有正地形展现出低阻特性,负地形展现出高阻特性,这样的影响随着地形的复杂程度不断提高显现的越发明显,甚至能抵消地下目标体所产生的真实电阻率影响,使得目标体的真实特征不能被探获。利用ANSYS软件进行快速地形改正有效地解决了地形起伏对电阻率的影响,使反映的介质电阻率回归为真实。

在使用ANSYS软件进行数值分析的过程中,网格剖分的大小直接影响着计算结果的精度,所以网格剖分越细,计算精度也就有所提高。但是,网格划分的大小也影响着计算工作量的大小,直接影响着计算速度。因此,在应用中应权衡这两种因素,做到综合考虑。图3为不同网格剖分密度下的山脊地形电阻率曲线。由图3可见,随着剖分尺度的大小变化,曲线的光滑程度也发生变化,且剖分尺度越大,曲线振荡程度越小,说明计算精度越高。在ANSYS稳态电流有限元分析中,程序没有提供远边界单元类型,在实际应用时,须构建一个远边界模型,对研究区和边界区进行不同的网格划分。

图3 二维纯山脊地形电阻率曲线图(不同剖分尺度)

3 应用实例

河南省西部山区是我国重要的有色金属成矿区带之一,该区域分布有多个大大小小的有色贵金属矿床,为河南省矿产资源需求提供了保障[13-14]。研究区位于河南省西部洛宁县,熊耳山北坡,其大地构造位置处于中朝准地台南缘、华熊台缘凹陷、崤山—鲁山拱褶断束中部(图4)。地层出露主要为太古宇太华群的角闪斜长片麻岩、熊耳群的安山岩和第四系。区内岩浆活动频繁,岩浆岩种类较多,其中以燕山期酸性岩最为发育,并与本区内生金属矿产的形成关系极为密切。断裂构造十分发育,并以NE—NNE向断裂为主,且规模较大。研究区围岩蚀变主要为硅化、绢云母化、绿泥石化、黄铁矿化、泥化等。硅化表现为石英沿破碎带充填,形成巨厚层的硅化碎裂岩带出现。构造蚀变带内矿化与蚀变主要有黄铁矿化、方铅矿化、碳酸盐化、高岭土化、褐铁矿化、硅化、绢云母化等,金属矿物表现形式为星点状、浸染状、细脉状、细脉侵染状、团块状。

图4 区域地质矿产图(a)与研究区地质简图(b)

表1为研究区电性参数统计结果,由表1可见,研究区矿石(含矿化体)与围岩存在着较大的物性差异,目标体电性表现特征为低阻高极化特征,为在该研究区开展物探工作提供物性前提。

表1 研究区电性参数统计结果

本次物探研究工作是在试验区K9号矿脉上进行,采用时间域激电偶极-偶极测深装置,点距20 m,隔离系数为5,MN极距40 m。仪器为V8,最大供电电流为14 A,6个电通道方式。图5为试验剖面测量结果,由图5a可见,极化率断面图,极化率高值异常与K9号矿脉基本对应,反映出地形对极化率测量结果影响较小,而图5b中所反映电阻率异常与矿脉不对应,并且展现出了高电阻率的特征,这一高电阻率特征是地下矿脉真实电阻率与地形影响产生的电阻率叠加影响产生的综合反映,其低阻异常向左侧地形相对低缓的一侧进行偏离。这是由于地形切割较大,山谷地形纯电阻率异常反映为高阻,致使矿体引起的低电阻异常被掩盖,发生偏离现象。图5c为采用ANSYS软件计算出该地形所产生的响应值与实际探测的电阻率值的比值进行地形改正计算,其结果可见低电阻异常位置向右侧偏移,并与K9号矿脉基本对应,彻底消除了复杂地形对探测结果的影响效应。结合图5a和图5c结果显示与岩矿石物性结果一致,表现为低阻高极率特征,证明采用ANSYS软件进行电阻率地形改正使测量结果回归真实反映,效果显著。

4 结论

(1)激发极化法是地质找矿中应用最为广泛的地球物理方法之一,但是其电阻率数据收复杂地形影响也最为严重,负地形将形成高电阻率反映,而正地形将产生低电阻率反映,从而实际物探测量所获得电阻率带有地形加载的影响,使得大地介质的真实电阻率得不到有效体现。

(2)ANSYS软件是一款能实现多场及多场耦合分析、实现前后处理、求解及多场分析统一数据库的一体化大型有限元分析软件,能快速准确进行多种场景数值模拟计算,能够对复杂地形条件下不同装置类型的电阻率开展二维、三维地形改正,利用该软件有限元分析功能模拟计算出纯地形下产生的电阻率响应,通过比值法将野外所获取的电阻率值中地形影响消除,获得大地介质真实电阻率,实现复杂地形对电阻率值的零影响效应,为物探解释的可靠性提供保障。

(3)应用ANSYS软件在研究区复杂地形条件下开展电阻率地形改正实际应用,对比结果显示,运用该软件进行地形改正方法可行,效果显著。