在石油物探中, 电法勘探占有重要地位。特别是在区域普查中, 它可用于研究区域构造, 确定沉积盆地基底起伏情况, 圈定含油气远景区等^[1-2]。而在进行电法勘探时, 发现地下岩石的导电性常常表现出非常明显的宏观各向异性特性, 例如板岩和页岩, 它们在层理方向与垂直层理方向的电导率之比达到了1.2~5.0^[3-5]。在解释这些区域采集到的地电数据时, 如果忽略了电阻率各向异性的影响, 往往会得到不精确甚至错误的地电构造信息^[6-7], 因此, 研究电阻率各向异性地电断面的正演算法具有重要的理论和实际意义。

对于直流电阻率各向异性数值模拟算法的研究, 国内外学者已经取得了一定的研究成果。WAIT^[8], LI等^[9]和PERVAGO等^[10]推导出了一维层状任意各向异性介质点电源电位的解析解, 并分析了各向异性对层状介质直流视电阻率的影响; 徐世浙^[3]实现了点源二维垂向各向异性地电断面的直流电场有限元算法, 并分析了各向异性对电测深曲线的影响; VERNER等^[11]基于规则四边形网格离散微分方程, 采用狄利克雷边界条件, 实现了二维直流电阻率任意各向异性有限差分数值模拟算法; BIBBY^[12]针对轴对称模型实现了直流电阻率有限元算法; ZHOU等^[13]利用高斯积分网格实现了二维和三维直流电阻率倾斜各向异性正演算法; LI等^[4]实现了三维直流电阻率任意各向异性有限元正演算法。以上这些二维、三维数值解法均基于规则的结构网格, 很难准确模拟地形起伏、倾斜界面等复杂的地质结构, 为此, WANG等^[14]利用基于非结构四面体网格的有限元法, 实现了三维直流电阻率任意各向异性正演。非结构网格不仅在模拟复杂地电结构方面相对于结构网格有明显的优势, 还特别适合模拟多尺度构造, 比如在大尺度的区域构造中存在很多小尺度的不均匀体。

本文基于完全非结构三角网格, 实现了二维直流电阻率倾斜各向异性自适应有限元数值模拟算法。算法利用对偶加权后验误差估计因子来控制非结构三角网格自动细化过程^[15-17], 网格剖分过程自动完成, 无需人工干预, 不仅可以很好地模拟起伏地形和复杂的地质构造, 还可以大大减小人为网格剖分所导致的误差, 增加了模拟的灵活性, 此外, 混合边界条件的使用能进一步提高数值解的精度。对多个理论模型进行了正演计算, 在验证算法有效性的同时, 深入分析了倾斜各向异性对直流电阻率测深的影响。

1 方法理论 1.1 控制方程和边界条件

在倾斜各向异性介质中, 电导率σ为张量, 建立如图 1所示的坐标系。

在x′y′z′坐标系中(x′与x平行, y′与y轴、z′与z轴的夹角为β), 假设沿x′, y′, z′方向上的电导率分别为σ_x, σ_y, σ_z, 则电导率张量σ′可以表示为:

$ \mathit{\boldsymbol{\sigma '}} = \left( {\begin{array}{*{20}{c}} {{\sigma _x}}&0&0\\ 0&{{\sigma _y}}&0\\ 0&0&{{\sigma _z}} \end{array}} \right) $

(1)

σ′绕x′轴旋转到与计算坐标系xyz重合, 得到新坐标系xyz中的电导率张量$ \underset{=}{\mathop{\mathit{\boldsymbol{ \sigma}}\rm{ }}}\, $=R_xσ′R_x^T, R_x= $ \left( \begin{matrix} 1 & 0 & 0 \\ 0 & \rm{cos}\mathit{\beta } & \rm{-sin}\mathit{\beta } \\ 0 & \rm{sin}\mathit{\beta } & \rm{cos}\mathit{\beta } \\ \end{matrix} \right)$。经过坐标旋转后, xyz坐标系中的电导率$\underset{=}{\mathop{\mathit{\boldsymbol{ \sigma}}\rm{ }}}\, $为:

$ \underline{\mathit{\boldsymbol{\underline \sigma }}} = \left( {\begin{array}{*{20}{c}} {{\sigma _x}}&0&0\\ 0&{{\sigma _y}{{\cos }^2}\beta + {\sigma _z}{{\sin }^2}\beta }&{\frac{1}{2}\left( {{\sigma _y} - {\sigma _z}} \right)\sin 2\beta }\\ 0&{\frac{1}{2}\left( {{\sigma _y} - {\sigma _z}} \right)\sin 2\beta }&{{\sigma _y}{{\sin }^2}\beta + {\sigma _z}{{\cos }^2}\beta } \end{array}} \right) $

(2)

当点电源I位于地表, 地下介质为电阻率各向异性时, 电位u满足的微分方程为^[3]:

$ \nabla \cdot \left( {\underline{\mathit{\boldsymbol{\underline \sigma }}} \nabla u} \right) = - 2I\delta \left( A \right) $

(3)

式中:δ(A)表示以A为中心的狄拉克函数。

将(2)式代入(3)式, 并考虑到$ \underset{=}{\mathop{\mathit{\boldsymbol{ \sigma}}\rm{ }}}\, $与x无关, 可得到:

$ \begin{array}{l} \frac{\partial }{{\partial y}}\left[ {\left( {{\sigma _y}{{\cos }^2}\beta + {\sigma _z}{{\sin }^2}\beta } \right)\frac{{\partial u}}{{\partial y}} + \frac{1}{2}\left( {{\sigma _y} - {\sigma _z}} \right) \cdot } \right.\\ \left. {\sin 2\beta \frac{{\partial u}}{{\partial z}}} \right] + \frac{\partial }{{\partial z}}\left[ {\frac{1}{2}\left( {{\sigma _y} - {\sigma _z}} \right)\sin 2\beta \frac{{\partial u}}{{\partial y}} + \left( {{\sigma _y}{{\sin }^2}\beta + } \right.} \right.\\ \left. {\left. {{\sigma _z}{{\cos }^2}\beta } \right)\frac{{\partial u}}{{\partial z}}} \right] + {\sigma _x}\frac{{{\partial ^2}u}}{{\partial {x^2}}} = - 2I\delta \left( A \right) \end{array} $

(4)

(4) 式是一个三维微分方程, 这里选择用傅里叶变换将它转换为二维问题。将(4)式关于x方向进行傅里叶变换, 得到:

$ \begin{array}{l} \frac{\partial }{{\partial y}}\left[ {\left( {{\sigma _y}{{\cos }^2}\beta + {\sigma _z}{{\sin }^2}\beta } \right)\frac{{\partial U}}{{\partial y}} + \frac{1}{2}\left( {{\sigma _y} - {\sigma _z}} \right) \cdot } \right.\\ \left. {\sin 2\beta \frac{{\partial U}}{{\partial z}}} \right] + \frac{\partial }{{\partial z}}\left[ {\frac{1}{2}\left( {{\sigma _y} - {\sigma _z}} \right)\sin 2\beta \frac{{\partial U}}{{\partial y}} + \left( {{\sigma _y}{{\sin }^2}\beta + } \right.} \right.\\ \left. {\left. {{\sigma _z}{{\cos }^2}\beta } \right)\frac{{\partial U}}{{\partial z}}} \right] - {\sigma _x}{k^2}U = - I\delta \left( A \right) \end{array} $

(5)

式中:U为波数域的电位值。

或者:

$ \nabla \cdot \left( {{\mathit{\boldsymbol{ \tau }}} \nabla U} \right) - {k^2}{\sigma _x}U = - I\delta \left( A \right) $

(6)

式中: $ \nabla = \frac{\partial }{{\partial \mathit{y}}}{\mathit{e}_\mathit{y}}{\rm{ + }}\frac{\partial }{{\partial \mathit{z}}}{\mathit{e}_\mathit{z}};\mathit{\boldsymbol{\tau = }}\left[{\begin{array}{*{20}{c}} {{\mathit{\tau }_{{\rm{11}}}}}&{{\mathit{\tau }_{{\rm{12}}}}}\\ {{\mathit{\tau }_{{\rm{21}}}}}&{{\mathit{\tau }_{{\rm{22}}}}} \end{array}} \right]$, τ₁₁=σ_ycos²β+σ_zsin²β, τ₂₂=σ_ysin²β+σ_zcos²β, τ₁₂=τ₂₁=1/2(σ_y-σ_z)sin2β; k为x方向的波数。

假定边界Γ_∞位于无穷远处, 使得电阻率各向异性不均匀体在该边界上的影响近似为零, 这时Γ_∞上的电位可近似按点电源的电位确定。假设Γ_∞周围介质是均匀各向异性的, 且各向异性主轴与坐标轴重合, 这时点电源的电位可表示为:

$ u = \frac{{c\sqrt {{\sigma _x}} }}{{\sqrt {{x^2} + {\lambda _y}{y^2} + {\lambda _z}{z^2}} }} $

(7)

其中, λ_y=σ_x/σ_y, λ_z=σ_x/σ_z, c为比例系数。对(7)式进行傅里叶变换:

$ \begin{array}{l} U = \int_0^\infty {\frac{{c\sqrt {{\sigma _x}} }}{{\sqrt {{x^2} + {\lambda _y}{y^2} + {\lambda _z}{z^2}} }}\cos kx{\rm{d}}x} \\ \;\;\;\; = c\sqrt {{\sigma _x}} {K_0}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right) \end{array} $

(8)

式中:K₀为第二类零阶修正贝塞尔函数。基于(8)式对y求偏导数, 得到:

$ \frac{{\partial U}}{{\partial y}} = - c\sqrt {{\sigma _x}} {K_1}k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} \frac{{k{\lambda _y}y}}{{\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} }} $

(9)

式中:K₁为第二类一阶修正贝塞尔函数。由(8)式可得$ \mathit{c = U/}\left[{\sqrt {{\mathit{\sigma }_\mathit{x}}} {\mathit{K}_{\rm{0}}}\left( {\mathit{k}\sqrt {{\mathit{\lambda }_\mathit{y}}{\mathit{y}^{\rm{2}}} + {\mathit{\lambda }_\mathit{z}}{\mathit{z}^{\rm{2}}}} } \right)} \right]$, 代入(9)式, 得到:

$ \frac{{\partial U}}{{\partial y}} + \frac{{{K_1}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right)}}{{{K_0}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right)}}\frac{{k{\lambda _y}y}}{{\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} }}U = 0 $

(10)

同理, 基于(8)式对z求偏导数, 并代入比例系数c, 有:

$ \frac{{\partial U}}{{\partial z}} + \frac{{{K_1}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right)}}{{{K_0}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right)}}\frac{{k{\lambda _z}z}}{{\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} }}U = 0 $

(11)

空气不导电, 地表Γ_s电流密度J的法向分量为零, 故有:

$ \mathit{\boldsymbol{J}}\left| {_n} \right. = {\left( { - \mathit{\boldsymbol{\sigma }}\nabla u} \right)_n}\left| {_{{\mathit{\Gamma }_s}}} \right. = 0 $

(12)

对(12)式进行傅里叶变换, 得到:

$ \mathit{\boldsymbol{\tau }}\frac{{\partial U}}{{\partial n}} = 0 $

(13)

(6) 式、(10)式、(11)式和(13)式组成了电阻率倾斜各向异性二维点电源直流电阻率波数域的边值问题。

1.2 有限元近似

将边值问题转化为有限元方程的方法主要有两种, 即加权余量法和变分法^[18], 前者利用虚功原理, 后者利用最小位能原理, 两种方法是等价的。本文利用加权余量法推导电阻率倾斜各向异性二维点电源直流电阻率边值问题的有限元方程, 将等式(6)两边乘以转换电位的任意小量δU, 并在模型区域Ω内积分, 可得:

$ \begin{array}{l} \int\limits_\mathit{\Omega } {\nabla \cdot \left( {\mathit{\boldsymbol{\tau }}\nabla U} \right)\delta U{\rm{d}}\mathit{\Omega }} - \int\limits_\mathit{\Omega } {{k^2}{\sigma _x}U\delta U{\rm{d}}\mathit{\Omega }} + \\ \;\;\;\;\;\;\int\limits_\mathit{\Omega } {I\delta \left( A \right)\delta U{\rm{d}}\mathit{\Omega }} = 0 \end{array} $

(14)

利用矢量恒等式∇·(τ∇U)δU=∇·(τ∇U·δU)-τ∇U·∇δU和高斯公式$ \oint\limits_{\mathit{L}}{\mathit{A}\cdot \mathit{n}\rm{d}\mathit{\Gamma =}\iint\limits_{\mathit{S}}{ \nabla \cdot \mathit{A}\rm{d}\mathit{x}\rm{d}\mathit{y}}}$, 由(14)式可得:

$ \begin{array}{l} \int\limits_\mathit{\Gamma } {\mathit{\boldsymbol{\tau }}\frac{{\partial U}}{{\partial n}}\delta U{\rm{d}}\mathit{\Gamma }} - \int\limits_\mathit{\Omega } {\mathit{\boldsymbol{\tau }} \nabla U \cdot \nabla \delta U{\rm{d}}\mathit{\Omega }} - \int\limits_\mathit{\Omega } {{k^2}{\sigma _x}U\delta U{\rm{d}}\mathit{\Omega }} + \\ \;\;\;\;\;\;\int\limits_\mathit{\Omega } {I\delta \left( A \right)\delta U{\rm{d}}\mathit{\Omega }} = 0 \end{array} $

(15)

将边界条件(10)式、(11)式、(13)式代入(15)式, 可得到:

$ \begin{array}{l} \int_\mathit{\Omega } {\left[ {\mathit{\boldsymbol{\tau }}\nabla U \cdot \left( {\nabla \delta U} \right) + {k^2}{\sigma _x}U\delta U} \right]{\rm{d}}\mathit{\Omega }} + \\ \int_{{\mathit{\Gamma }_\infty }} {k{\sigma _x}\frac{{{K_1}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right)}}{{{K_0}\left( {k\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} } \right)}}\frac{{y\cos \theta + z\sin \theta }}{{\sqrt {{\lambda _y}{y^2} + {\lambda _z}{z^2}} }} \cdot } \\ \;\;\;\;\;U\delta U{\rm{d}}\mathit{\Gamma } = \int_\mathit{\Omega } {I\delta \left( A \right)\delta U{\rm{d}}\mathit{\Omega }} \end{array} $

(16)

其中, θ为Γ_∞的外法向方向与y轴的夹角。将模型区域Ω离散成若干个三角单元, 则(16)式的积分等于所有三角单元的积分总和, 形成如下线性方程组:

$ \mathit{\boldsymbol{KU}} = \mathit{\boldsymbol{P}} $

(17)

式中:U是所有网格节点处电位值构成的列向量; P=(0, …, 0, I/2, 0, …, 0)^T, 只有电源点对应的节点处其值为I/2;系数矩阵K对称且高度稀疏。为了节省内存空间, 采用按行压缩存储方式, 并利用直接求解器SuperLU^[19]求解线性方程组(17)。对于若干个波数k, 解方程(17)即可求得波数域的电位U, 再对U进行反傅里叶变换^[20], 即可求得空间域所有网格节点处的电位u, 再利用三角单元中3个节点处的电位值进行线性插值可以求得接收点处的电位值。根据不同的测量装置^[21], 由接收点处的电位值就可以计算出视电阻率值。

2 基于后验误差估计的网格自动细化

本文中用到的基于后验误差估计的网格自动细化的基本理论在文献[15]中已经有详细论述, 自适应网格细化有限元正演流程如图 2所示。首先, 输入初始网格和模型参数数据, 并计算初始网格的有限元近似解; 然后, 利用有限元近似解计算初始网格中每个三角单元的局部误差估计值和相应对偶问题的有限元解; 最后, 计算每个三角单元的后验误差估计因子值, 对后验误差估计因子值比较大的三角单元进行网格细化; 再次计算网格细化后的有限元解, 并且重复迭代上述过程, 直到有限元数值解收敛为止。

3 模型计算

下面的模型计算中, 设定所有模型的走向方向均为x方向, 采用二维网格剖分器triangle^[22]剖分网格, 每次网格细化的三角单元数占三角单元总数的20%, 网格细化的终止条件是相邻两次有限元解的相对误差小于1%。所有计算均在3.20 GHz CPU, 4G RAM的个人计算机上完成。

3.1 模型一

图 3为两层电阻率各向同性模型, 其参数为:ρ_x¹=ρ_y¹=ρ_z¹=5 Ω·m, h₁=10 m, β₁=0;ρ_x²=ρ_y²=ρ_z²=50 Ω·m, β₂=0。图 4为两层电阻率各向异性模型, 其参数为:ρ_x¹=ρ_y¹=10 Ω·m, ρ_z¹=2.5 Ω·m, h₁=10 m, β₁=0;ρ_x²=ρ_y²=ρ_z²=50 Ω·m, β₂=0。

采用二极装置, 设电源点位于坐标原点。接收电极沿着y轴正方向布置, 极距从1m逐渐增加到50m, 共有20个接收电极。为了选择适合本模型的离散化波数, 我们采用阮百尧等^[23]提出的适合单个点电源情形的波数和加权系数进行计算, 如表 1所示。模拟网格大小为300m×800m, 经过11次迭代后, 网格剖分细化终止, 有限元数值解收敛。具体的网格剖分信息和计算结果如图 5和图 6所示。

从图 5b可以看出, 电源点处和接收点处的网格加密程度要明显高于其它区域, 实现了网格局部加密的目的, 这不仅可以使得有限元数值解快速收敛于精确解, 还能减少计算资源要求, 提高计算效率。从图 6可以看出, 两层电阻率各向同性和各向异性模型的二维直流电阻率数值解都与解析解^[8]吻合得很好, 除电源点附近数值解的相对误差较大外, 其它接收点处的相对误差都在1%以内。

从电阻率各向同性模型和电阻率倾斜各向异性模型的视电阻率曲线可以看出(图 6a), 对于电阻率倾斜各向异性模型, 虽然第一层的平均电阻率$ \mathit{\rho =}\sqrt{{{\mathit{\rho }}_{\mathit{x}}}\cdot {{\mathit{\rho }}_{\mathit{z}}}}$=5 Ω·m与电阻率各向同性模型的相同, 但它们的电测深视电阻率曲线却有很大的不同, 倾斜各向异性层的影响使得电阻率倾斜各向异性模型的视电阻率值要明显偏大。这说明, 在解释电阻率各向异性介质的电测深曲线时, 采用平均电阻率近似的方法不合适, 必须考虑电阻率各向异性介质的实际情况, 才能得到正确的解释结果。

3.2 模型二

假定在ρ₀=5 Ω·m的低阻均匀介质中埋藏有一个电阻率倾斜各向异性正方体, 并假定其在走向x方向上无限延伸, y方向和z方向上的边长都为5m, 其电阻率为ρ_xx=ρ_zz=100 Ω·m, ρ_yy=5 Ω·m, 不均匀体距离水平地面的距离为0.5m, 如图 7所示。

采用单极装置, 设电源点位于坐标原点, 接收电极沿着y轴正负两个方向布设, 其距坐标原点的距离为:0.45, 0.90, 1.35, 1.80, 2.25, 2.70, 3.15, 3.60, 4.05, 4.50, 4.95, 5.40m。我们依然采用阮百尧等^[23]提出的波数和加权系数(表 1), 自适应有限元模拟网格大小为300m×600m。分别计算电阻率各向异性倾角β=30°, 45°, 60°, 75°情形下的视电阻率测深曲线, 结果如图 8所示。从图 8可以看出, 当电阻率各向异性倾角不同时, 视电阻率曲线异常的峰值会发生变化, 随着电阻率各向异性倾角的增大, 视电阻率曲线的峰值会变小, 可以用地下电流密度的分布特征来解释这种现象。在二维倾斜各向异性介质中, 电流密度是电场的横向分量和纵向分量的线性组合^[24]。这时, σ_yx=σ_yz=0, yz平面内y, z方向上电流密度分别为:

$ {j_y} = \left( {0,y.z} \right) = - {\sigma _{yy}}\frac{{\partial u}}{{\partial y}} - {\sigma _{yz}}\frac{{\partial u}}{{\partial z}} $

(18)

$ {j_z} = \left( {0,y.z} \right) = - {\sigma _{yz}}\frac{{\partial u}}{{\partial y}} - {\sigma _{zz}}\frac{{\partial u}}{{\partial z}} $

(19)

式中:σ_yz, σ_yy, σ_zz分别表示测点电导率σ在不同方向上的分量, $ \mathit{\boldsymbol{ \sigma}}\rm{ }\rm{=}\left( \begin{matrix} {{\mathit{\sigma }}_{\mathit{xx}}} & {{\mathit{\sigma }}_{\mathit{xy}}} & {{\mathit{\sigma }}_{\mathit{xz}}} \\ {{\mathit{\sigma }}_{\mathit{yx}}} & {{\mathit{\sigma }}_{\mathit{yy}}} & {{\mathit{\sigma }}_{\mathit{yz}}} \\ {{\mathit{\sigma }}_{\mathit{zx}}} & {{\mathit{\sigma }}_{\mathit{zy}}} & {{\mathit{\sigma }}_{\mathit{zz}}} \\ \end{matrix} \right)$。

将接收电极置于地面以下, 并按照(18)式和(19)式分别计算电流密度y方向和z方向上的分量。

图 9显示了当电阻率各向异性倾角β=30°, 45°, 60°, 75°时, 电流密度矢量在yz平面内的分布情况。因为电流密度在远离电源点的地方衰减很快, 电流密度的分布与r²(r表示测点与点电源I之间的距离)成反比, 所以, 在远离电源点的区域电流密度值比较小。电流密度主要沿着低阻介质流动, 其大致流动方向和低阻地层倾斜的方向一致, 角度约等于倾角, 所以, 随着低阻地层倾角的增大, 地面测得的电位差就会变小, 计算出来的视电阻率值也就会随之变小。

3.3 模型三

假设在ρ₀=100 Ω·m的地下均匀介质中埋藏有一个低阻倾斜各向异性长方体, 它在x方向上无限延伸, y方向和z方向上的边长分别为6m和2m, 其电阻率为ρ_xx=ρ_yy=4 Ω·m, ρ_zz=25 Ω·m, 长方体距离水平地面的距离为2m, 如图 10所示。利用PIDLISECKY等^[25]的最优化波数计算程序, 获得了适用于极距范围0.25~220.00m, 双点电源情形下的7个波数和加权系数(表 2)。

图 11显示了当间距因子n=5时偶极装置不同各向异性倾角情形的视电阻率曲线, 各向异性长方体的主轴电阻率ρ_xx, ρ_yy, ρ_zz为4, 4, 25 Ω·m, 各向异性倾角β=0, 30°, 60°, 90°, 电源电极和测量电极距为AB=MN=2m, 测量范围为-20~20m。从图 11可以看出, 各向异性倾角β对视电阻率的影响很大。尽管二维模型是对称的, 但当β=30°, 60°时, 视电阻率曲线却不对称, 视电阻率曲线的极小值偏离了中心位置并偏向了左边, 这时的视电阻率曲线与向右倾斜的电阻率各向同性长方形板体的视电阻率曲线非常相似。

4 结论

本文基于完全非结构三角网格, 实现了二维直流电阻率倾斜各向异性自适应有限元数值模拟算法。利用对偶加权后验误差估计因子控制非结构三角形网格自动细化过程, 网格剖分过程自动完成, 无需人工干预, 不仅可以很好地模拟起伏地形和复杂的地质构造, 还可以大大减小人为网格剖分所导致的误差。

电阻率各向异性倾角不仅会影响视电阻率曲线的峰值大小, 还会影响视电阻率曲线的对称性, 从而对异常体的形态和分布位置的解释都产生影响, 所以, 在地质解释时考虑地下介质电阻率各向异性的影响很有必要。利用电流密度矢量图可以方便地解释二维倾斜各向异性模型视电阻率曲线的特征。