E-mail Alert Rss
 

物探与化探, 2022, 46(4): 887-896 doi: 10.11720/wtyht.2022.1377

方法研究·信息处理·仪器研制

基于Fractal模型的复电阻率法2.5D有限元数值模拟

龙秀洁,1, 陈汉波,2, 莫亚军1, 区小毅1, 卢胜辉1

1.广西壮族自治区地球物理勘察院,广西 柳州 545005

2.桂林理工大学 地球科学学院,广西 桂林 541006

Fractal model-based 2.5 D finite element modeling of complex resistivity method

Long Xiu-Jie,1, Chen Han-Bo,2, Mo Ya-Jun1, Ou Xiao-Yi1, Lu Sheng-Hui1

1. Guangxi Geophysical Investigation Institute,Liuzhou 545005,China

2. College of Earth Sciences, Guilin University of Technology,Guilin 541006,China

通讯作者: 陈汉波,(1990-),男,博士后,主要从事电磁场数值模拟与反演研究工作。Email:chenhanbo@glut.edu.cn

责任编辑: 王萌

收稿日期: 2021-07-9   修回日期: 2021-10-11  

基金资助: 桂林理工大学科研启动基金(RD2100002165)
中国博士后科学面上基金(2021MD703820)

Received: 2021-07-9   Revised: 2021-10-11  

作者简介 About authors

龙秀洁,(1989-),男,工程师,从事地球物理勘查及相关研究工作。Email: 496178747@qq.com

摘要

复电阻率法在油气资源、矿产勘查中发挥了重要的作用,为了深刻认识复电阻率法异常特征变化规律,本文对复电阻率2.5D正演问题展开研究。首先直接给出复电阻率法2.5D有限元正演所满足的变分问题,并详细地推导相应的刚度矩阵的计算过程。引入Fractal模型作为等效模型研究频谱激电异常特征。对单元内的复电导率及复电位均进行线性插值,而后,采用不完全LU分解的稳定双共轭梯度算法求解有限元线性方程组,获得异常复电位值。设计3个典型的地电模型验证了本文算法的正确性及精确性,并分析了不同装置下,不同频率的2.5D复电阻率异常响应特征。数值模拟结果表明,采用Fractal模型研究激发极化异常特征是可行、有效的;不同装置、不同频率下的复电阻率法异常特征有着显著的差异。

关键词: 有限元; 复电导率; Fractal模型; 异常复电位

Abstract

This study proposed the variational problem of 2.5D finite element forward modeling of the complex resistivity method and detailed the process of solving stiff matrix of finite equations. The Fractal model was introduced as a research model for studying the equivalent induced polarization anomalies of spectra. Furthermore, the complex conductivity and complex potential of a grid unit were linearly interpolated. Then, to obtain anomalous complex potential, finite element linear equations were solved using the biconjugate gradient stabilized method with incomplete LU decomposition. The results of three typical geoelectric models validated the correctness and accuracy of the algorithm proposed in this study. Furthermore, this study analyzed the abnormal response characteristics of 2.5D complex resistivity under different frequencies.

Keywords: finite element method; complex resistivity; Fractal model; anomalous complex potential

PDF (7643KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

龙秀洁, 陈汉波, 莫亚军, 区小毅, 卢胜辉. 基于Fractal模型的复电阻率法2.5D有限元数值模拟[J]. 物探与化探, 2022, 46(4): 887-896 doi:10.11720/wtyht.2022.1377

Long Xiu-Jie, Chen Han-Bo, Mo Ya-Jun, Ou Xiao-Yi, Lu Sheng-Hui. Fractal model-based 2.5 D finite element modeling of complex resistivity method[J]. Geophysical and Geochemical Exploration, 2022, 46(4): 887-896 doi:10.11720/wtyht.2022.1377

0 引言

复电阻率法是一种重要的地球物理勘探方法,具有电性参数多的特点,在金属矿勘探领域起着重要的作用,较其他物探方法具有抗干扰能力强,多参数对比解释可提供更丰富的异常信息的优点。近年来,该方法被广泛应用于石油地震勘探、水文地质调查等领域[1]

国内外学者对于复电阻率法数值模拟做了大量的研究,取得众多成果。如国内的罗延钟等[2-3]、张桂清等[4]、刘菘等[5]。国外的研究者,如G.W. Hohmann等[6]、A. Weller等[7]、D.W. Oldenburg等[8]等。然而,大多对于频谱激电法的研究均基于Core-Core模型,且没有考虑复电导率呈线性变化的情况。因此,本文在前人的基础上,基于异常复电位法,采用Fractal模型[9]研究复电导率呈线性变化的复电阻率法2.5D有限元数值模拟,并对其相应的异常响应特征进行分析,以期能为复电阻率法野外实际工作提供理论依据。

1 复电阻率法边值问题

1.1 Fractal 模型

Fractal 模型是由B.R.P. Rocha提出用来描述岩、矿石激发极化现象的数学模型,B.R.P. Rocha等[9]对大量岩、矿石标本及露头进行测量,结果表明,岩、矿石的激发极化响应所引起的复电阻率可用以下公式进行表述:

ρ(ω)=ρ01-m1-11+1+uδr(1+v)rh,

其中:ρ0为岩、矿石未极化的直流电阻率;m为介质的极化率(无量纲);δr为矿物颗粒电阻率百分比;rh=1/(1+iωτ0);u=iωτ(1+v); v=(iωτf)-η; τ为双层振荡相关的松弛时间常数;τ0为采样松弛时间常数;τf为分形松弛时间;η为介质的分形几何直接相关的参数(无量纲)。

式(1)对于定量描述岩矿石的激发极化现象提供了一种新的手段,相对于常用的core-core模型来说,其提供了更多的岩石物性参数,从而可以更详细地解释岩矿石激发极化效应所产生的原因及影响因素。图1展示了不同η的复电阻率的振幅和相位随着频率的变化规律。计算参数为:ρ0=100 Ω·m、m=0.5、δr=1.0、τ=10-6 s、τ0=10-12 s、τf=100、η分别取0.05、0.20、0.35、0.5、0.65、0.80、0.95。由图1可以看得到复电祖率ρ(ω)幅值随着频率的增大而递减,在低频时,且η取较小值时,复电祖率幅值趋向于ρ0。而相位φ的绝对值与η成正比,在低频及高频时趋向于0。

图1

图1   不同η时的复电阻率振幅(a)及相位频谱(b)

Fig.1   The amplitude (a) and phase (b) of complex resistivity for different frequency exponent


1.2 基于Fractal模型的复电阻率变分问题

点源条件下,复电阻率法2.5D变分问题如下所示:

Fu-=Ω12σu-2+12σk2u-2+σ'u-0u-+σ'k2u-0u-dΩ+ΓkK1krAcosrA,n-K1krBcosrB,nK0krA-K0krB12σu-2+σ'u-0u-dΓδFu-=0

其中:Ω为研究区域;Γ为无穷远边界;σ为介质的复电导率;σ0为均匀介质的复电导率(电源点处的电导率);σ'=σ-σ0 为介质的异常复电导率,k为波数;K0K1分别为零阶和1阶贝塞尔函数;rArB分别为观测点距源点的距离;u-0为正常复电位。可通过以下公式进行计算:

u-0=IK0(krA)-IK0(krB)2πσ0,

式中:I为电流大小。

2 有限单元法分析

2.1 有限单元剖分

采用矩形网格进一步细分为4个三角形网格剖分方式对研究区域进行剖分,如图2所示:将方程(2)对研究区域Ω和无穷边界Γ的积分分为各个单元的积分之和:

F(u-)=ΩΩ12σ(u-)2+12σk2u-2dΩ+ΩΩ[σ'u-0u-+σ'k2u-0u-]dΩ+ ΩΓk[K1(krA)cos(rA,n)-K1(krB)cos(rB,n)]K0(krA)-K0(krB)12σu-2dΓ+ ΩΓk[K1(krA)cos(rA,n)-K1(krB)cos(rB,n)]K0(krA)-K0(krB)σ'u-0u-dΓ δFu-=0 

图2

图2   有限单元网格剖分示意

Fig.2   Sketch of finite element mesh


在单元内对复电位及复电导率进行线性插值。线性插值的母单元与子单元如图3所示。则单元内的复电位及复电导率则可表示为

u-=N1u-1+N2u-2+N3u-3,

σ=N1σ1+N2σ2+N3σ3,

图3

图3   线性插值的母单元(a)及子单元(b)

Fig.3   Parent element (a) and sub-element(b) of linear interpolation


其中:N1N2N3为插值形函数;u-1u-2u-3为三角单元节点复电位;σ1σ2σ3为单元节点复电导率。形函数的具体形式如下所示:

N1= 12Δ(a1x+b1y+c1);N2= 12Δ(a2x+b2y+c2);N3=12Δa3x+b3y+c3),

其中:a1=y2-y3 ;a2=y3-y1 ;a3=y1-y2 ;b1=x3-x2 ;b2=x1-x3 ;b3=x2-x1 ;c1=x2y3-x3y2 ;c2=x3y1-x1y3 ;c3=x1y2-x2y1 ;Δ=12(a1b2-a2b1) 为三角单元面积。

2.2 单元分析

方程(4)中的第1项

Ω12σ(u-)2dΩ=Ω12σu-x2+u-y2dxdy=12u-eTK1eu-e,

其中:

K1e=(σi+σj+σm)12Δaiai+bibiaiaj+bibjaiam+bibmaiaj+bibjajaj+bjbjajam+bjbmaiam+bibmajam+bjbmamam+bmbm,

方程(4)积分第2项

Ω12σk2u-2dΩ=12σk2u-eTΩNNTdΩu-e=12u-eTK2eu-e,

其中

K2e=k260Δ6σi+2σj+2σm2σi+2σj+σm2σi+σj+2σm2σi+2σj+σm2σi+6σj+2σmσi+2σj+2σm2σi+σj+2σmσi+2σj+2σm2σi+2σj+6σm,

方程(4)积分第3项

Γk[K1(krA)cos(rA,n)-K1(krB)cos(rB,n)]K0(krA)-K0(krB)12σu-2dΓ=12u-eTK3eu-e,

其中:

K3e=k[K1(krA)cos(rA,n)-K1(krB)cos(rB,n)]K0(krA)-K0(krB)3σi+σjσi+σj0σi+σjσi+3σj0000,

方程(4)积分第4项

Ωσ'u-0u-dΩ=u-eTK1e'u-e0,

其中:

K'1e=σ'i+σ'j+σ'm12Δaiai+bibiaiaj+bibjaiam+bibmaiaj+bibjajaj+bjbjajam+bjbmaiam+bibmajam+bjbmamam+bmbm,

方程(4)积分第5项

Ωσ'k2u-0u-dΩ=u-eTK'2eu-e0,

其中:

K'2e=k260Δ6σ'i+2σ'j+2σ'm2σ'i+2σ'j+σ'm2σ'i+σ'j+2σ'm2σ'i+2σ'j+22σ'i+6σ'j+2σ'mσ'i+2σ'j+2σ'm2σ'i+σ'j+2σ'mσ'i+2σ'j+2σ'm2σ'i+2σ'j+6σ'm,

方程(4)积分第6项

Γk[K1(krA)cos(rA,n)-K1(krB)cos(rB,n)]K0(krA)-K0(krB)σ'u-0u-dΓ=u-eTK'3eu-e0,

式中:

K'3e=k[K1(krA)cos(rA,n)-K1(krB)cos(rB,n)]K0(krA)-K0(krB)3σ'i+σ'jσ'i+σ'j0σ'i+σ'jσ'i+3σ'j0000,

2.3 总体合成

在单元内,将式(9)、式(11)、式(13)、式(15)、式(17)、式(19)的积分结果进行相加后,再扩展成由全体节点组成的矩阵或列阵:

F(u-)= 12u-T(K1e+K2e+K3e) u-+ u-T(K1e'+ K2e'+ K3e') u-0= 12u-TKeu-+ u-TKe'u-0,

式中:u-Tu-0为所有场值的列向量;KeK'e为扩展矩阵。

2.4 求解线性方程组

对方程(20) 进行求变分:

δF(u-)=δ u-TKeu-u-TK'eu-0,

δF(u-)=0,有

δ u-TKeu-=-δ u-TKe'u-0,

又因δu-T≠0,有

Keu-=K'eu-0

求解方程(22)便可得波数域中的异常复电位值。本文采用不完全LU分解的稳定双共轭梯度算法来求解线性方程组,具体求解过程请详见文献[10],在此不再赘述。求解出异常复电位后,再通过傅里叶反变换便可得到三维空间中的异常复电位。

2.5 复电阻率的计算

求解出三维空间中的异常复电位,便可按照以下公式计算复电阻率及相位值:

ρs=k u-(M)-u-(N)I0,
φs=arctan Re(ρs)Im(ρs)· 180π,

其中:k为装置系数;u-Mu-N为接收电极MN处异常复电位;ρ0为按式(1)计算出来的均匀介质的复电阻率值。

3 模型计算

3.1 模型1

为了检验文中算法及所编制的程序的正确性,设置一个两层水平层状极化介质模型,厚度为6 m,电阻率为100 Ω。Fractal模型参数如表1所示。供电频率f=0.125 Hz,将2.5D有限元数值解与一维数字滤波算法所计算出来的解析解进行对比,如图4所示。由图4可知,复电阻率的实分量和虚分量数值解与解析解基本一致,从而验证了本文所述算法及所编制的程序的正确性。

表1   模型1参数

Table 1  parameters of the first model

nh/mρ0/(Ω·m)mδrτ/μsτf/msητ0/ps
1
2
5
100
10
0.906
0.885
4.959
4.877
23.43
18.72
10
100
0.20
0.44
1.0
1.0

新窗口打开| 下载CSV


图4

图4   2.5维有限元数值解与一维数字滤波解析解对比

Fig.4   Comparison of 2.5D FEM result with numerical filter wave numerical solution


3.2 模型2

模型2 为极化半空间中有一异常体。具体几何参数如图5所示,Fractal模型参数如表2所示。图6图7为频率为0.125、0.25、0.5、1.0 Hz的复电阻率法2.5D有限元数值模拟结果。由图6图7可看到,复电阻率的各分量都很好地反映出异常体的存在,且异常响应范围与异常体的埋藏位置基本吻合。在计算的频率范围内,复电阻率的相位及虚分量均为负值,偶极装置下的复电阻率振幅和相位均呈“八”字型异常特征,特别值得注意的是复电阻率的虚分量,其幅值随着频率的增大而减小,且远小于实分量。

图5

图5   模型2断面示意

Fig.5   Cross section of second model


表2   模型2的Fractal 模型参数

Table 2  Fractal parameters of the second model

Kindρ0/(Ω·m)mδrτ/μsτf/msητ0/ps
半空间
异常体
5200
1470
0.906
0.885
4.959
4.877
23.43
18.72
10
100
0.20
0.44
1.0
1.0

新窗口打开| 下载CSV


图6

图6   给定频率的偶极—偶极测量方式复电阻率振幅(左列)和相位(右列)拟断面

a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz

Fig.6   The pseudo section of amplitude and phase of apparent complex resistivity in dipole-dipole array for giving the frequencies


图7

图7   给定频率的对称四极测量方式复电阻率实部(左列)和虚部(右列)拟断面

a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz

Fig.7   The pseudo section of real and imaginary component of apparent complex resistivity in four-pole array for giving the frequencies


3.3 模型3

为了进一步验证本文算法的有效性,设计一个极化半空间中有双异常体模型。异常体大小均为12 m×6 m,埋深维4 m,两者相距12 m。异常体电导率为0.01 s/m,背景介质的电导率为0.001 s/m。具体几何参数如图8所示,Fractal模型参数如表3所示。图9~12为频率为0.125、0.25、0.5、1.0 Hz的不同测量装置下的复电阻率法2.5D有限元数值模拟结果。由图9~12可看到,复电阻率的各分量都很好地反映出异常体的存在,且异常响应范围与异常体的埋藏位置基本吻合。在计算的频率范围内,4种测量装置下的视复电阻率的虚分量均为负值;对称四极装置下,复电阻率的实部和虚部呈现双“柱体”型异常特征,如图9所示;如图10所示,偶极装置下的复电阻率实部和虚部均呈双“八”字型异常特征;从图11可知,三极测量装置,复电阻率的实部和虚部的异常特征与对称四极类似,不过虚部的右侧部分的值大于左侧;如图12可知,二极装置下,可以看到复电阻率的异常位置与模型非常吻合,其纵向分辨率较于三极和对称四极测量装置高。同样,我们可以观察到,在4种测量装置下,复电阻率的实部随着频率的增大,产生的变化不大;不同于实部分量,复电阻率的虚分量,其幅值随着频率的增大而减小,且远小于实分量。

图8

图8   模型3断面示意

Fig.8   Cross section of second model


表3   模型3的Fractal 模型参数

Table 3  Fractal parameters of the second model

Kindρ0/
(Ω·m)
mδrτ/μsτf/msητ0/ps
半空间10000.9064.95923.43100.201.0
异常体
(左)
1000.8854.87718.721000.441.0
异常体
(右)
1000.8854.87718.721000.441.0

新窗口打开| 下载CSV


图9

图9   给定频率的对称四极测量方式复电阻率实部(左列)和虚部(右列)拟断面

a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz

Fig.9   The pseudo section of real and imaginary component of apparent complex resistivity in four-pole array for giving the frequencies


图10

图10   给定频率的偶极—偶极测量方式复电阻率实部(左列)和虚部(右列)拟断面

a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz

Fig.10   The pseudo section of real and imaginary component of apparent complex resistivity in dipole-dipole array for giving the frequencies


图11

图11   给定频率的三极剖面测量方式复电阻率实部(左列)和虚部(右列)拟断面

a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz

Fig.11   The pseudo section of real and imaginary component of apparent complex resistivity in pole-dipole array for giving the frequencies


图12

图12   给定频率的二极测量方式复电阻率实部(左列)和虚部(右列)拟断面

a、b—0.125 Hz;c、d—0.25 Hz;e、f—0.5 Hz;g、h—1.0 Hz

Fig.12   The pseudo section of real and imaginary component of apparent complex resistivity in pole-pole array for giving the frequencies


4 结论

本文实现了复电阻率法2.5D有限元数值模拟,采用异常复电法,避免了源的奇异性的影响,提高了数值解的精度;引入Fractal 模型对频率域激发极化响应进行研究,可以为更好解释激发极化效应提供一种新的手段;对复电导率及复电位进行线性插值,可以模拟更为真实的野外地质情况;分析了所设计的地电模型的复电祖率异常响应特征,可为野外实际勘探提供理论依据。数值模拟结果验证了本文基于Fractal 模型的CR2.5D有限元数值模拟算法的正确性和精确性,表明引入Fractal 模型研究复电导率呈线性变化的CR响应特征是可行的、有效的。

参考文献

杨振威, 郑伟, 李晓斌, .

频谱激电法的发展与展望

[J]. 物探与化探, 2015, 39(1):22-28.

[本文引用: 1]

Yang Z W, Zhen W, Li X B, et al.

The development and prospect of the spectral induced polarization method

[J]. Geophysical and Geochemical Exploration, 2015, 39(1):22-28.

[本文引用: 1]

罗延钟, 张桂青. 频率域激电法原理[M]. 北京: 地质出版社, 1988.

[本文引用: 1]

Luo Y Z, Zhang G Q. Principle of frequency domain IP method[M]. Beijing: Geological Publishing House, 1988.

[本文引用: 1]

罗延钟, 方胜.

视复电阻率频谱的一种近似反演方法

[J]. 地球科学, 1986, 11(1) : 9-102.

[本文引用: 1]

Luo Y Z, Fang S.

An approximate invrtsion of the apparent complex resistivity spectrum

[J]. Earth Science, 1986, 11(1) : 9-102.

[本文引用: 1]

张桂青, 崔先文, 罗延钟.

一种反演频谱激电法视频谱求取真参数的方法

[J]. 地质与勘探, 1987, 23(4) : 48-54.

[本文引用: 1]

Zhang G Q, Cui X W, Luo Y Z.

The Determination of intrinsic parameters by inversing IP apparent spectrum

[J]. Geology and Exploration, 1987, 23(4) : 48-54.

[本文引用: 1]

刘崧, 官善友, 高鹏飞.

求极化椭球体真 Cole-Cole 参数的联合谱激电反演

[J]. 地球物理学报, 1994, 37(S1) : 542-551.

[本文引用: 1]

Liu S, Guan S Y, Gao P F.

Joint sip inversion for estimation of intrinsiccole-cole parameters of apolarizable ellipsoid

[J]. Chinese Journal of Geophysice, 1994, 37(S1) : 542-551.

[本文引用: 1]

Hohmann G W.

Three-dimensional induced polarization and electromagnetic modelling

[J]. Geophysics, 1975, 40 (2) :309-324.

DOI:10.1190/1.1440527      URL     [本文引用: 1]

Weller A, Seichter M, Kampke A.

Induced-polarization modelling using complex electrical conductivities

[J]. Geophys. J. Int., 1996, 127(1):387-398.

DOI:10.1111/j.1365-246X.1996.tb04728.x      URL     [本文引用: 1]

Oldenburg D W, Li Y.

Inversion of induced polarization data

[J]. Geophysics, 1994, 59 (9) :1327-1341.

DOI:10.1190/1.1443692      URL     [本文引用: 1]

Rocha B R P, Habashy T M. Fractal geometry, porosity and complex resistivity. I: From rough pore interfaces to hand specimens, in: Developments in Petrophysics[M], London: Geological Society Publishing House, London, 1995.

[本文引用: 2]

柳建新, 蒋鹏飞, 童孝忠, .

不完全LU分解预处理的BICGSTAB算法在大地电磁二维正演模拟中的应用

[J]. 中南大学学报:自然科学版, 2009, 40(2):484 -491.

[本文引用: 1]

Liu J X, Jiang P F, Tong X Z, et al.

Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling

[J]. Journal of Central South University: Science and Technology, 2009, 40(2):484 -491.

[本文引用: 1]

/

京ICP备05055290号-3
版权所有 © 2021《物探与化探》编辑部
通讯地址:北京市学院路29号航遥中心 邮编:100083
电话:010-62060192;62060193 E-mail:whtbjb@sina.com