井中激发极化(井中激电)时将供电装置/测量装置/供电装置和测量装置放入井中使其靠近异常体, 以提高激发极化强度或激发极化响应强度。相较井井观测和地井观测, 井地观测具有观测数据量大, 可精确反演的优点。目前主要采用井地电阻率法来圈定油气藏的边界^[1], 张天伦^[2-3]分别在井中无套管与井中有套管的情况下进行物理水槽实验以确定油气藏边界, 并证明了利用井地电阻率法确定异常体边界的可行性; 张天伦等^[4-6]利用非无穷远处三极剖面法确定油气藏边界, 并在新疆等地成功发现了小块油气藏; 汤井田等^[7]利用井地电阻率法得到的歧离率确定高电阻率油藏边界。

激发极化的正、反演以直流电阻率法正、反演为基础。直流电阻率正演主要采用有限差分法^[8-15]、有限单元法^[16-22]和边界单元法^[23-25]。有限差分法受限于曲面边界的处理而有限单元法的优势在于处理曲面边界; 边界单元法适用于纯地形因素的起伏地表模拟, 具有速度快、精度高与节约内存等优点, 但无法实现复杂模型的模拟; 有限单元法作为一种主要的直流电阻率正演方法, 处理复杂的模型和边界时有较强的灵活性和适应性, 可以自动满足其内部边界条件。在实际勘探中, 地形的起伏会造成高电阻率或低电阻率干扰, 并对反演结果造成不可预计的影响, 有限单元法适合模拟复杂物性和起伏地表模型。直流电阻率正演存在电源奇异性的问题, 即源点对应节点处的电位无穷大, 引起模拟误差增大, 目前解决该问题的方法主要有异常电位法^{[14, 26]}和源点附近局部加密网格^{[19, 20]}两种方法。在源点附近加密网格虽然提高了局部的模拟精度, 但增加了计算量, 此外网格结构的任意性也破坏了有限单元法固有的收敛结构, 在提高局部精度同时降低了整体的收敛速度^[27]。LOWRY等^[26]提出了异常去除技术, 该技术将总电位分为一次场电位(含奇异值部分)和二次场电位(非奇异值部分), 一次场电位如平坦地形下的均匀半空间、水平层状介质或者垂直接触面等电位, 可用简单模型的背景场电位方程解析解来替代, 将该解析解代入原总电位满足的方程, 可消除方程右端电源项狄拉克函数引起的奇异性。起伏地表条件下, 简单模型的背景场电位方程无解析解, 因此无法直接使用异常电位法。BLOME^[25]利用边界单元法对背景场电位方程单独求解; PENZ^[28]重新定义了异常电位法, 背景场电位仍然为平坦地形下均匀半空间电位方程的解析解, 将起伏地表与地下不均匀体电位方程作为异常场电位方程, 将平坦地形下均匀半空间电位方程的解析解作为背景场方程的解, 采用广义有限差分法对异常电位的新边界条件进行求解。

目前对井中激发极化法的应用以及正、反演研究非常少^[29-31], 大多数的电阻率法与激发极化法的正、反演研究都是基于平坦地形的。但实际应用中地形的影响不可忽略, 因此本文基于PENZ^[28]的研究, 将改进的异常电位法引入起伏地表条件下井中激电的正演模拟。基于二维线源场, 采用加权余量法得到有限单元积分方程, 并引入直接解法求解器求解最终的线性方程组。将几个典型的起伏地表模型试算结果与采用传统的总电位法得到的结果进行比较, 验证了本方法的正确性与有效性。

1 井中激电正演模拟算法 1.1 直流电阻率法二维线源边值问题

在直角坐标系中, 假设Z轴垂直向下, 源点A的坐标为(x_A, z_A), 则在二维线源场中总电位满足以下方程:

$ \begin{array}{l} \frac{\partial }{{\partial x}}\left[ {\sigma \left( {x,z} \right)\frac{{\partial u\left( {x,z} \right)}}{{\partial x}}} \right] + \frac{\partial }{{\partial z}}\left[ {\sigma \left( {x,z} \right)\frac{{\partial u\left( {x,z} \right)}}{{\partial z}}} \right]\\ \;\;\; = - I\delta \left( {x - {x_A}} \right)\delta \left( {z - {z_A}} \right)\;\;\;\;\;{\rm{in}}\;\mathit{\Omega } \end{array} $

(1)

式中:δ是狄拉克函数; σ为地层模型的电导率; Ω为整个计算区域; in表示在区域内。求解式(1)还需要相应的边界条件, 边界分为地表边界Γ_s与无穷远边界Γ_∞。因地表无电流流向空气中, 故满足纽曼边界条件:

$ \frac{{\partial u}}{{\partial \mathit{\boldsymbol{n}}}} = 0\;\;\;\;{\rm{on}}\;{\mathit{\Gamma }_{\rm{s}}} $

(2)

式中: n为边界Γ_s的外法向方向; Γ_∞为无穷远边界; on表示在边界上。COGGON^[32]采用了以下两种边界条件:

$ u = 0 $

(3)

$ \frac{{\partial u}}{{\partial \mathit{\boldsymbol{n}}}} = 0 $

(4)

试算结果表明, 采用(3)式计算得到的电位偏小, 采用(4)式计算得到的电位偏大, 因此DEY等^[10]引入了以下混合边界条件:

$ \frac{{\partial u}}{{\partial \mathit{\boldsymbol{n}}}} + \frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}u = 0\;\;\;{\rm{on}}\;{\mathit{\Gamma }_\infty } $

(5)

式中: n为边界Γ_∞的外法向方向; r为源点到无穷远边界Γ_∞的距离; r为源点到无穷远边界Γ_∞的有向距离; (r, n)为r, n之间的夹角。(1)式, (2)式和(5)式共同称为总电位的边值问题。

1.2 改进的异常电位法

(1) 式右端项包含的狄拉克函数会造成源点处的奇异值, 从而引起源点附近的模拟误差增大。目前有两种解决该问题的方法, 一种为加密源点附近的网格^[20]; 另一种为文献[26]提出的异常消除技术。该技术将总电位u分离为一次场电位u_p(奇异值部分)与二次场电位u_s(非奇异值部分), 一次场电位为平坦地形下的均匀半空间、层状空间或者垂直接触面的电位等, 可使用简单模型的背景场电位方程解析解来替代, 总电位定义为:

$ u = {u_{\rm{p}}} + {u_{\rm{s}}} $

(6)

异常电位的边值问题如下^{[14, 18, 33]}:

$ \left\{ {\begin{array}{*{20}{c}} {\nabla \cdot \left( {\sigma \nabla {u_{\rm{s}}}} \right) = - \nabla \cdot \left[ {\left( {\sigma - {\sigma _0}} \right)\nabla {u_{\rm{p}}}} \right]}&{{\rm{in}}\;\mathit{\Omega }}\\ {\frac{{\partial {u_{\rm{s}}}}}{{\partial \mathit{\boldsymbol{n}}}} = 0}&{{\rm{in}}\;{\mathit{\Gamma }_{\rm{s}}}}\\ {\frac{{\partial {u_{\rm{s}}}}}{{\partial \mathit{\boldsymbol{n}}}} + \frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}{u_{\rm{s}}} = 0}&{{\rm{in}}\;{\mathit{\Gamma }_\infty }} \end{array}} \right. $

(7)

式中:σ为介质电导率; σ₀为背景场电导率^[14]。当源点在地面时, 平坦地形条件下均匀半空间的解析解为u_p=I/(πσ₀)ln(1/r₀), 其中r₀为源点到观测点的距离。(7)式不包含点源项, 且在起伏地表条件下背景场电位方程无解析解, 通常采用边界单元法^[23-25]来进行数值求解, 增加了计算量。在起伏地表条件下, 本文采用PENZ等^[28]提出的方法, 将总电位分为一次场电位和二次场电位, 一次场电位仍然按平坦地形条件下的正常电位计算u_p=1/(2πσ₀r₀), 修改自然边界条件, 定义新的自然边界条件为:

$ \frac{{\partial {u_{\rm{s}}}}}{{\partial \mathit{\boldsymbol{n}}}} = - \frac{{\partial {u_{\rm{p}}}}}{{\partial \mathit{\boldsymbol{n}}}} = \frac{I}{{{\rm{ \mathsf{ π} }}{\sigma _0}r}}\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)\;\;\;\;{\rm{in}}\;{\mathit{\Gamma }_{\rm{s}}} $

(8)

由(8)式可知, 一次场电位满足非齐次的边界条件, 理论上说明存在流向空气的电流, 将这部分的电位定义为“虚拟电位”, 因为实际中不存在流向空气的电流^[34]。我们认为二次场电位u_s是由起伏地表与地下不均匀体产生的异常电位。由(5)式, (6)式和(7)式可得无穷远边界条件为^[34]:

$ \frac{{\partial {u_{\rm{s}}}}}{{\partial \mathit{\boldsymbol{n}}}} + \frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}{u_{\rm{s}}} = - \frac{{\partial {u_{\rm{p}}}}}{{\partial \mathit{\boldsymbol{n}}}} - \frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}{u_{\rm{p}}} = 0\;\;\;\;{\rm{in}}\;{\mathit{\Gamma }_\infty } $

(9)

(7) 式、(8)式和(9)式组成了起伏地表条件下基于线源的异常电位法的边值问题, 通常采用有限差分法、有限单元法等方法求解该类边值问题, 本文采用有限单元法并使用非结构化网格对地电模型及起伏地表进行剖分近似。

1.3 有限单元-加权余量法

有限单元法求解变分问题, 必须先得到对应的有限单元积分方程, 将满足边界条件的微分方程转换为对应的积分方程, 通常采用两种方法:里兹(Ritz)法和加权余量法。里兹法利用变分原理将边界条件转化为等价泛函的极值问题; 加权余量法引入权函数将微分方程转化为弱形式的积分方程, 当基函数为权函数时称为伽辽金(Galerkin)法^[33]。变分问题的推导较复杂, 依据标准的变分原理只能处理自伴正定, 且边界条件为第一类或第三类齐次的微分算子方程, 大量的偏微分方程均无对应的泛函。通常采用较为简单的Galerkin方法求解变分问题, 前人的研究表明, 如果边值问题存在相应的泛函, 采用Galerkin法会得到与采用Ritz法相同弱形式的积分方程。

将区域Ω剖分成一系列的三角形单元e, 节点i, j, m为任意3个节点, 如图 1所示。

假设单元内的电位是线性分布的, 则单元中的电位u可以表示为:

$ u = {N_i}{u_i} + {N_j}{u_j} + {N_m}{u_m} $

(10)

式中:u_i, u_j, u_m分别为i, j, m这3个节点上的电位; N_i, N_j, N_m分别为3个节点对应的形函数:

$ {N_h} = \frac{1}{{2\Delta }}\left( {{a_n}x + {b_n}y + {c_n}} \right),\;\;\;\;n = i,j,m $

其中, 中间变量a_i, a_j, a_m, Δ分别为:

$ \left\{ \begin{array}{l} {a_i} = {y_i} - {y_m},{b_i} = {x_m} - {x_n}\\ {a_j} = {y_m} - {y_i},{b_j} = {x_i} - {x_m}\\ {a_m} = {y_i} - {y_j},{b_m} = {x_j} - {x_i}\\ \Delta = \frac{1}{2}\left( {{a_i}{b_j} - {a_j}{b_i}} \right) \end{array} \right. $

根据加权余量法原理, (7)式的余量为:

$ R = \nabla \cdot \left( {\sigma \nabla {u_{\rm{s}}}} \right) + \nabla \cdot \left[ {\left( {\sigma - {\sigma _0}} \right)\nabla {u_{\rm{p}}}} \right] $

(11)

将形函数N_i作为权函数, 得到三角形单元e的加权余量积分表示并令其等于0:

$ \begin{array}{*{20}{c}} {\int_{{\mathit{\Omega }^{\rm{e}}}} {N_i^{\rm{e}}R{\rm{d}}x{\rm{d}}y} = 0,}&{i = 1,2,3} \end{array} $

(12)

式中:N_i^e为三角形单元e在节点i的形函数。

将(11)式代入(12)式可得:

$ \begin{array}{l} \int_{{\mathit{\Omega }^{\rm{e}}}} {N_i^{\rm{e}}\left[ {\frac{\partial }{{\partial x}}\left( {\sigma \frac{{\partial {u_{\rm{s}}}}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}\left( {\sigma \frac{{\partial {u_{\rm{s}}}}}{{\partial z}}} \right) + \frac{\partial }{{\partial x}}\left( {\sigma '\frac{{\partial {u_{\rm{p}}}}}{{\partial x}}} \right) + } \right.} \\ \;\;\;\;\;\;\left. {\frac{\partial }{{\partial z}}\left( {\sigma '\frac{{\partial {u_{\rm{p}}}}}{{\partial z}}} \right)} \right]{\rm{d}}x{\rm{d}}y = 0,\;\;\;\;i = 1,2,3 \end{array} $

(13)

式中:σ′为异常电导率, σ′=σ-σ₀。(13)式可表示为:

$ \begin{array}{l} \begin{array}{*{20}{c}} {\frac{\partial }{{\partial x}}\left( {\sigma \frac{{\partial u}}{{\partial x}}} \right) \cdot N_i^{\rm{e}} = \frac{\partial }{{\partial x}}\left( {\sigma \frac{{\partial u}}{{\partial x}} \cdot N_i^{\rm{e}}} \right) - \sigma \frac{{\partial u}}{{\partial x}}\frac{{\partial N_i^{\rm{e}}}}{{\partial x}}}\\ {\frac{\partial }{{\partial z}}\left( {\sigma \frac{{\partial u}}{{\partial z}}} \right) \cdot N_i^{\rm{e}} = \frac{\partial }{{\partial z}}\left( {\sigma \frac{{\partial u}}{{\partial z}} \cdot N_i^{\rm{e}}} \right) - \sigma \frac{{\partial u}}{{\partial z}}\frac{{\partial N_i^{\rm{e}}}}{{\partial z}}} \end{array}\\ \;\;\;\;\;\;\;\;\;\;\;i = 1,2,3 \end{array} $

(14)

根据散度定理, 将(10)式代入(14)式可得:

$ \begin{array}{l} \int_{{\mathit{\Omega }^{\rm{e}}}} { - \sigma \left( {\frac{{\partial N_i^{\rm{e}}}}{{\partial x}}\frac{{\partial N_j^{\rm{e}}}}{{\partial x}} + \frac{{\partial N_i^{\rm{e}}}}{{\partial z}}\frac{{\partial N_j^{\rm{e}}}}{{\partial z}}} \right)u_{{\rm{si}}}^{\rm{e}}{\rm{d}}x{\rm{d}}z} - \\ \;\;\;\int_{{\mathit{\Omega }^{\rm{e}}}} {\sigma '\left( {\frac{{\partial N_i^{\rm{e}}}}{{\partial x}}\frac{{\partial N_j^{\rm{e}}}}{{\partial x}} + \frac{{\partial N_i^{\rm{e}}}}{{\partial z}}\frac{{\partial N_j^{\rm{e}}}}{{\partial z}}} \right)u_{{\rm{pi}}}^{\rm{e}}{\rm{d}}x{\rm{d}}z} + \\ \;\;\;\int_\mathit{\Gamma } {\sigma \frac{{\partial {u_{\rm{s}}}}}{{\partial n}}N_i^{\rm{e}}{\rm{d}}\mathit{\Gamma }} + \int_\mathit{\Gamma } {\sigma '\frac{{\partial {u_{\rm{p}}}}}{{\partial n}}N_i^{\rm{e}}{\rm{d}}\mathit{\Gamma }} = 0\\ \;\;\;i,j = 1,2,3 \end{array} $

(15)

式中:N_j^e为三角形单元e在节点j的形函数; u_pi^e, u_si^e分别为i节点的一次场和二次场电位; Γ为区域Ω的边界, Γ=Γ_s+Γ_∞。结合边界条件(8)和(9)式可得到以下方程:

$ \begin{array}{l} \int_{{\mathit{\Omega }^{\rm{e}}}} {\sigma \left( {\frac{{\partial N_i^{\rm{e}}}}{{\partial x}}\frac{{\partial N_j^{\rm{e}}}}{{\partial x}} + \frac{{\partial N_i^{\rm{e}}}}{{\partial z}}\frac{{\partial N_j^{\rm{e}}}}{{\partial z}}} \right)u_{{\rm{si}}}^{\rm{e}}{\rm{d}}x{\rm{d}}z} + \int_{{\mathit{\Omega }^{\rm{e}}}} {\sigma '\left( {\frac{{\partial N_i^{\rm{e}}}}{{\partial x}} \cdot } \right.} \\ \;\;\;\;\left. {\frac{{\partial N_j^{\rm{e}}}}{{\partial x}} + \frac{{\partial N_i^{\rm{e}}}}{{\partial z}}\frac{{\partial N_j^{\rm{e}}}}{{\partial z}}} \right)u_{{\rm{pi}}}^{\rm{e}} - \int_{{\mathit{\Gamma }_{\rm{s}}}} {\sigma \frac{{I\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{{\rm{ \mathsf{ π} }}r{\sigma _0}}}N_i^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_{\rm{s}}}} + \\ \;\;\;\;\int_{{\mathit{\Gamma }_\infty }} {\sigma \frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}N_i^{\rm{e}} \cdot N_j^{\rm{e}} \cdot u_{{\rm{si}}}^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_\infty }} + \\ \;\;\;\;\int_{{\mathit{\Gamma }_{\rm{s}}}} {\sigma '\frac{{I\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{{\rm{ \mathsf{ π} }}r{\sigma _0}}}N_i^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_{\rm{s}}}} + \int_{{\mathit{\Gamma }_\infty }} {\sigma '\frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}N_i^{\rm{e}} \cdot } \\ \;\;\;\;N_j^{\rm{e}} \cdot u_{{\rm{pi}}}^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_\infty } = 0,\;\;\;\;i,j = 1,2,3 \end{array} $

(16)

(16) 式可表示为以下矩阵形式:

$ \sigma \left( {\mathit{\boldsymbol{K}}_1^{\rm{e}} + \mathit{\boldsymbol{K}}_2^{\rm{e}}} \right)\mathit{\boldsymbol{u}}_{\rm{s}}^{\rm{e}} = - \sigma '\left( {\mathit{\boldsymbol{K}}_1^{\rm{e}} + \mathit{\boldsymbol{K}}_2^{\rm{e}}} \right)\mathit{\boldsymbol{u}}_{\rm{p}}^{\rm{e}} + {\mathit{\boldsymbol{g}}^{\rm{e}}} $

(17)

式中: u_p^e, u_s^e分别为三角形单元e的一次场和二次场电位向量; K₁^e, K₂^e, g^e均为中间变量矩阵; 对应的矩阵元素K₁^e, K₂^e, g^e分别为:

$ \left\{ \begin{array}{l} K_1^{\rm{e}} = \int_{{\mathit{\Omega }^{\rm{e}}}} {\left( {\frac{{\partial N_i^{\rm{e}}}}{{\partial x}}\frac{{\partial N_j^{\rm{e}}}}{{\partial x}} + \frac{{\partial N_i^{\rm{e}}}}{{\partial z}}\frac{{\partial N_j^{\rm{e}}}}{{\partial z}}} \right){\rm{d}}x{\rm{d}}z} \\ K_2^{\rm{e}} = \int_{{\mathit{\Gamma }_\infty }} {\frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}N_i^{\rm{e}} \cdot N_j^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_\infty }} \\ {g^{\rm{e}}} = \int_{{\mathit{\Gamma }_{\rm{s}}}} {\left( {\sigma - \sigma '} \right)\frac{{I\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{{\rm{ \mathsf{ π} }}r{\sigma _0}}}N_i^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_{\rm{s}}}} \\ \;\;\;\; = \int_{{\mathit{\Gamma }_{\rm{s}}}} {\frac{{I\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{{\rm{ \mathsf{ π} }}r}}N_i^{\rm{e}}{\rm{d}}{\mathit{\Gamma }_{\rm{s}}}} \end{array} \right. $

(18)

由文献[33]可得:

$ \mathit{\boldsymbol{K}}_1^{\rm{e}} = \frac{\sigma }{{4\Delta }}\left[ {\begin{array}{*{20}{c}} {{a_i}}&{{b_i}}\\ {{a_j}}&{{b_j}}\\ {{a_m}}&{{b_m}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{a_i}}&{{a_j}}\\ {{b_i}}&{{b_j}} \end{array}} \right] $

假设边界落在ij上, 有线积分∫_lN_i^aN_j^bdΓ=[(a!b!)/(a+b+1)!]l, 其中a, b为非负整数, l为ij边的长度, 故有:

$ \mathit{\boldsymbol{K}}_2^{\rm{e}} = \frac{l}{6}\frac{{\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{r\ln \left( {1/r} \right)}}\left[ {\begin{array}{*{20}{c}} 2&1&0\\ 1&2&0\\ 0&0&0 \end{array}} \right] $

$ {\mathit{\boldsymbol{g}}^{\rm{e}}} = \frac{l}{2}\frac{{I\cos \left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{n}}} \right)}}{{{\rm{ \mathsf{ π} }}r}}\left[ {\begin{array}{*{20}{c}} 1&1&0 \end{array}} \right] $

令K_e=σ(Ke₁+Ke₂), K′ _e=σ′(Ke₁+Ke₂), Ke₁, Ke₂, g^e为单元节点组成的矩阵并相加, (17)式可表示为以下的线性方程组:

$ \mathit{\boldsymbol{K}}{\mathit{\boldsymbol{u}}_{\rm{s}}} = - \mathit{\boldsymbol{K'}}{\mathit{\boldsymbol{u}}_{\rm{p}}} + \mathit{\boldsymbol{g}} $

(19)

式中: u_p, u_s分别为全体单元节点的一次场电位和二次场电位矩阵; K=∑K_e, K′=∑K′ _e, g=∑g^e, 求解此线性方程组得到二次场电位u_s, 将其代入(6)式, 得到最终的总电位。

当地下介质无极化特性时, 正演模拟得到的一次场电位为u₁; 当地下电阻率为ρ的介质存在极化特性时, 正演模拟得到极化场总电位为u, 两次正演模拟结束之后, 得到视极化率如下(基于二级装置):

$ {\eta _{\rm{s}}} = \frac{{u - {u_1}}}{u} \times 100\% = \frac{{{u_2}}}{u} \times 100\% $

(20)

式中:u₂为二次场电位。

2 数值算例及分析

网格剖分是有限单元法重要的环节, 因为非结构化网格可以局部加密网格, 因而其在精确模拟各种复杂模型时有着天然的优势。对于非结构化网格, 我们利用IntelMKL库提供的Pardiso求解器求解线性方程组。为了检验改进的异常电位法的效果, 首先利用二维线源条件下特殊地形存在解析解这一优势, 来验证该方法的正确性; 然后选择两个理论模型:山谷与山脊模型检验该方法对起伏地表的计算效果; 最后分析地形对二次场电位的影响规律, 并与总电位法得到的结果进行对比分析。

2.1 算法验证

在二维线源的情况下, 经保角变换可得到特殊地形下的解析解, 保角变换的原理如下:

取变换函数:

$ \omega = \sqrt {{z^2} + {d^2}} $

(21)

式中:z为Z平面里的任意一点的坐标, ω为W平面里任意一点的坐标, d为如图 2a所示的边界DE长度。

将Z平面(图 2a)的边界CDEFG及其下部平面变换成W平面的实轴及其下部平面(图 2b), 由(21)式可得:

$ z = \sqrt {{\omega ^2} - {d^2}} $

(22)

将z=x+iy, w=u+iv代入式(22)并展开得:

$ \begin{array}{l} x = \\ \pm \sqrt {\frac{{\sqrt {\left( {{u^2} - {v^2} - {d^2}} \right) + {{\left( {2uv} \right)}^2}} + \left( {{u^2} - {v^2} - {d^2}} \right)}}{2}} \\ y = \\ \pm \sqrt {\frac{{\sqrt {\left( {{u^2} - {v^2} - {d^2}} \right) + {{\left( {2uv} \right)}^2}} - \left( {{u^2} - {v^2} - {d^2}} \right)}}{2}} \end{array} $

(23)

式中:x的正负与u相同; y的正负与v相同。对于给定的d和w(u, iv), 可根据(23)式计算w(u, iv)的对应点z在Z平面内的坐标(x, iy); 对于W平面内的直线R′S′, 可根据(23)式得到Z平面内对应的曲线RS。W平面内(相当于平坦地形下)的电位为:

$ u = \frac{I}{{{\rm{ \mathsf{ π} }}\sigma }}\ln \left( {\frac{1}{r}} \right) $

(24)

式中:I为供电电流。

图 3为二维纯山谷地形, 厚度为20 m, 电流强度为1 A, 电导率为1.0 S/m, 五角星为源位置。表 1为改进的异常电位法求得的总电位与解析解、总电位法的结果对比。该异常电位法得到的结果和解析解的误差小, 同时提高了源点附近的总电位精度。

2.2 计算实例

模型1:图 4, 图 5分别为山谷与山脊模型, 五角星均为源位置, 均匀介质, 各项参数如下:电阻率为1 Ω·m, 三极测量装置中测量电极M, N变化范围为0~200 m, 测量电极距2 m, 山谷模型网格剖分的计算区域为200 m×100 m, 扩展边界区域为10 000 m, 三角形网格单元数为21 151, 节点数为10 864;山脊模型网格剖分的计算区域为200 m×100 m, 扩展边界长度为10 000 m, 三角形网格单元数为21 261, 节点数为10 919。

从图 6与图 7可以看出采用总电位法得到的山谷模型视电阻率曲线出现了高阻异常, 山脊模型的视电阻率曲线出现了低阻异常, 这与前人研究得出的地形与视电阻率呈镜像对称的规律相符^[35]; 在源点附近, 采用总电位法得到的视电阻率出现了畸变, 这也验证了改进的异常电位法的正确性, 与总电位法得到的结果相比, 改进的异常电位法得到的视电阻率精度更高。

模型2:该山谷模型的各项参数如图 8所示, 围岩电阻率ρ₀=100 Ω·m, 极化率η₀=0.01, 低阻异常体电阻率ρ₁=5 Ω·m, 极化率η₁=0.2, 异常体宽度20 m, 厚度10 m, 埋深10 m, 网格剖分的计算区域为200 m×100 m, 扩展边界长度为10 000 m, 生成三角形网格单元数为16 169, 节点数为8 253。五角星为源在井中的位置, 采用井地二极观测方式, 视电阻率计算公式为:ρ_s=πu/ln(1/r)。

采用改进的异常电位法对该模型进行正演模拟, 得到了线源在不同深度时的视电阻率曲线, 为了便于比较起伏地表条件下异常体的视电阻率响应特征, 还给出了纯山谷地形的视电阻率曲线和包含低阻异常体的平坦地形视电阻率曲线(图 9a, 图 9b)。从图 9a可以看出, 山谷地形造成了高阻异常的畸变, 随着源深度的增加, 山谷地形造成的高阻异常逐渐变小; 图 9b为包含低阻异常体的平坦地形(异常体的各项参数与模型2一致)随着源深度的增加视电阻率的变化规律, 源在异常体上方时, 视电阻率曲线表现出高阻特征, 源在异常体顶端时(源在地下20 m处), 视电阻率最大; 源在异常体下方时, 视电阻率曲线表现出低阻特征, 当源在异常体底部时(源在地下30 m处), 视电阻率的幅值达到最小, 随着源深度增加, 视电阻率曲线趋于平缓; 从图 9c可以看出, 地形因素造成的高阻异常与低阻异常体产生的异常叠加, 掩盖了原有低阻异常体的视电阻率变化规律, 造成了高阻异常的假象, 这会产生错误的地震解释成果; 当源位于地下20 m处时, 分别采用总电位法和改进的异常电位法对模型2进行正演模拟, 为了便于对比分析, 又采用总电位法对纯山谷地形与含低阻异常体的平坦地形(异常体的各项参数与模型2一致)进行正演模拟, 共获得了4条视电阻率曲线, 如图 9d所示, 可以看出起伏地表条件下, 总电位法与异常电位法的结果几乎一致, 再次验证了异常电位法的正确性与有效性, 纯山谷地形导致的视电阻率曲线异常与含有低阻异常体的平坦地形产生的异常相似, 这也意味着起伏地表对视电阻率影响非常大。

从图 10a的视极化率曲线可以看出, 源在低阻高极化异常体(异常体的各项参数与模型2一致)上方时, 视极化率低于围岩的极化率, 为低值异常; 源在该异常体顶端时, 视极化率变化辐值达到最大; 源在异常体下方时, 视极化率曲线变化趋势相反。图 10b为模型2的视极化率曲线, 视极化率的变化规律几乎和图 10a一致, 仅值不同。图 10c与图 10d分别为含低阻高极化异常体的平坦地形与模型2的视极化率断面图, 从图 10c中可以看到, 两个异常的边界几乎对应了低阻高极化异常体的上、下边界(白色虚线为异常体位置); 图 10d中两个异常边界与图 10c中的异常边界几乎一样, 这也说明了山谷地形未导致激电假异常, 仅对视极化率的值产生影响。

3 结论

针对起伏地表条件下传统的异常电位法中背景场电位方程无解析解这一问题, 从二维线源出发, 采用有限单元法, 对起伏地表条件下的二维地电模型进行正演模拟, 在前人的研究基础之上, 重新定义异常电位法, 将起伏地表与地下不均匀体视作异常场, 采用加权余量法对新的边界条件进行推导, 得到了对应的有限单元积分方程, 引入直接法求解器, 求解最终的线性方程组得到总电位。将二维线源条件下山谷模型的解析解以及总电位法得到的结果与改进的异常电位法得到的结果对比, 验证了该方法的正确性与有效性。含低阻高极化异常体的山谷模型正演模拟结果表明, 山谷地形改变了原有的低阻高极化异常体变化规律, 造成了高阻假异常, 但未导致激电假异常, 仅影响视极化率的值。