Anomaly edge enhancement and topographic correction technology of linear source 3D borehole-to-surface electrical method
-
摘要: 基于非结构网格有限元方法开展了三维复杂地电模型的线源井地电法的高效正演模拟研究,探讨了通过求取电场响应导数来刻画目标体边界范围、采用差异场地形校正技术来消除地形影响等措施对井地电法成像的效果和精度的影响。并通过对比与解析解,验证了本文数值解算法的有效性。模型计算结果表明:积水巷道的空间位置和走向均会引起视电阻率的显著变化,视电阻率变化率的极值准确且清晰地指示了巷道边界的位置;电位的归一化总水平导数极大地提高了井地电法对目标体复杂边界位置的识别能力;地形对井地电场分布的影响也很大,其视电阻率响应与地形形状近似呈对称关系,利用差异场技术能有效地削弱地形对井地电法高精度成像的影响。Abstract: Based on the finite element method of unstructured grid, the efficient forward modeling of the borehole-to-surface electrical method derived by the linear current source under the condition of the 3D complex geoelectric model was carried out. The effects on the effectiveness and accuracy of the borehole-to-surface electrical method imaging were discussed by obtaining the electric field response derivative to characterize the boundary range of the target body, and using the difference field topography correction technology to eliminate the topographic influence. And the comparison between the numerical solution and the analytical solution verifies the effectiveness of the algorithm in this paper. The model calculation results show that the spatial position and direction of the roadway with water accumulation cause significant changes in the apparent resistivity, and the extreme value of the apparent resistivity change rate accurately and clearly indicates the boundary position of the roadway. The normalized total horizontal derivative of the electric potential greatly improves the ability of the borehole-to-surface electrical method to identify the complex boundary position of the target body. Moreover, the influence of topography on the distribution of borehole-to-surface electrical field is also serious, and its apparent resistivity response is approximately symmetrical to the shape of the topography. The difference field technique can effectively weaken the influence of topography on the high-precision imaging of the borehole-to-surface electrical method. The research results have important theoretical and practical significance for improving the data interpretation level and application effect of the borehole-to-surface electrical method.
-
前言
井地电法通过井中激发、地面接收的方式对地下目标体进行成像,该方法综合了测井与地面两种电法勘探的优势,被广泛应用于油气资源、固体矿产、工程和环境等勘察领域(Xiong et al,2020)。相较于常规的地面电阻率方法,井地电法通过井眼将电流延伸至更靠近勘探目标体的位置,这样更有利于对深部或复杂构造地区目标体的成像。针对深部小规模目标体,井地电法的异常信号可能非常微弱,导致目标体的响应特征不明显或边界识别模糊不清,以至于很难准确地圈定目标体的分布范围(戴前伟等,2013;王智,潘和平,2014),从而影响到对资源量的准确评估或对安全生产事故与环境污染的及时防治。在复杂地形起伏地区开展电磁勘探时,地形是影响电磁场分布形态的重要因素,压制地形影响是地形起伏地区地质勘探取得良好效果的关键技术之一(薛国强等,2016)。另外,随着地质勘探目标复杂度的不断提高以及对实测资料精细化处理与解释的需要,继续采用一维、二维简单模型对构造复杂地区的数据进行处理和解释已无法满足工业生产的需求(Spitzer,1995)。因此,开展三维复杂模型下井地电法的边界识别与地形校正方法研究,对提高井地电法的勘探效果具有重要的理论意义和应用价值。
边界识别是位场资料处理与解释的常用手段(Li et al,2014),但是在同属于位场范畴的井地电法中的研究及应用实例却很少。在井地电法勘探中,地质体电性横向分界面对应地表处的电流密度变化剧烈,为利用边界识别技术刻画地质体的电性分界面提供了条件。戴前伟等(2013)应用异常电位的导数来圈定目标体的边界范围;王智和潘和平(2014)通过异常电位的归一化总水平导数增强了目标体的响应特征和边界位置反映;杨沁润等(2020)和Yang等(2021)应用电位的梯度参数使目标体在地表的反映更加显著。汤井田等(2007)利用视电阻率的歧离率来确定高阻油气藏边界,但该方法对井旁盲矿的探测能力很弱。
此外,我国多山区丘陵,在这些山地开展电磁勘探,地形对电磁场响应的影响往往不容忽视(Xiong et al,2019)。尽管早在二十世纪七十年代,Ku等(1973)就指出地形对电磁场的分布存在影响,然而已有文献中针对井地电法的理论研究大多仍然是基于平坦地表的假设,在地形校正方面的研究则更少见报道。张天伦和张伯林(1995)通过水槽模拟试验给出了一种消除井地电法中地形或浅层非均匀体等干扰的方法;李长伟等(2010)通过有限元方法实现了三维带地形模型的井地电位响应计算;王智等(2018)模拟了含地形的二维极化介质模型的井地电法响应。
如何有效地离散复杂地电模型对于提高正演计算的精度和效率至关重要。相较于结构化网格,非结构化网格能更好地适应不规则的复杂模型。鉴于此,本文将基于非结构网格有限元方法开展线源三维井地电法的数值计算研究。文中通过采用对井地电法响应求取导数的方法来刻画目标体边界范围,并将差异场处理技术(Kong et al,2008;Xiong et al,2020)应用到井地电法中进行地形校正,以加强对边界的刻画和对地形的压制,提升井地电法的计算效率和成像效果。
1. 方法
在笛卡尔坐标系中,设点电流源的位置为
${\boldsymbol{r}}_0= $ (x0,y0,z0),电流强度为I,则空间中任意位置的电流密度${\boldsymbol{J}} $ 满足$$ \nabla \cdot {\boldsymbol{J}} = I{\text{δ}} ( {{r_0}} ) \text{,} $$ (1) 式中,
${\boldsymbol{J}} $ 为传导电流密度,单位为A/m2;δ为Dirac函数。稳态场的电场强度
${\boldsymbol{E}} $ 可表示为电位u的负梯度,即$$ {\boldsymbol{E}} = - \nabla u {\text{.}} $$ (2) 由式(1)和式(2)并结合微分形式的欧姆定律,得到点源电位u满足的微分方程为
$$ \nabla \cdot ( {\sigma \nabla u} ) = - I{\text{δ}} ( {{r_0}} ) \text{,} $$ (3) 式中σ为电导率,单位为S/m。
由于空气的电阻率无穷大,因此没有电流穿过地空界面
$\varGamma ^{\rm{s}}$ 流入空气中,即地表处电流密度的法向分量为零,故地表处电位u满足第二类边界条件:$$ \frac{{\partial u}}{{\partial n}} = 0 {\text{.}} $$ (4) 将侧边界和下边界(Γ∞)取得足够远,使电位u在其上衰减为零,即
$$ u=0,$$ (5) 则由式(3)—(5)可得到地下稳恒点电流源的边值问题:
$$ \left\{ {\begin{array}{*{20}{l}} {\nabla \cdot ( \sigma \nabla u ) = - I{\text{δ}} ( {r_0} ) }, \\ {\begin{array}{*{20}{l}} {\dfrac{{\partial u}}{{\partial n}} = 0}&{{\varGamma ^{\rm{s}}}} , \end{array}} \\ {\begin{array}{*{20}{l}} {u = 0}&{{\varGamma ^\infty }} {\text{.}} \end{array}} \end{array}} \right. $$ (6) 依据变分原理以及第一标量格林定理,方程(6)对应的变分问题为
$$ \left\{ {\begin{array}{*{20}{l}} { F ( u ) = \dfrac{1}{2}\displaystyle\iiint_V {\nabla u \cdot ( {\sigma \nabla u} ) {\rm{d}}V - \displaystyle\iiint_V {uI{\text{δ}} ( {r_0} ) {\rm{d}}V}}}, \\ { {\text{δ}} F ( u ) = 0} ,\\ {\begin{array}{*{20}{c}} {u = 0}&{{\varGamma ^\infty }} \end{array}} {\text{.}} \end{array}} \right. $$ (7) 井地电法通过钢套管向地层供电时,套管直径远小于套管长度,因而可将钢套管简化为有限长线电流源模型来处理(谭河清等,2004)。有限长线源产生的电位可以看成是将线源划分为n段小线源产生的电位的叠加,而每段小线源均可当作点源处理。根据叠加原理,线源激励的电位变分问题为
$$ \left\{ {\begin{array}{*{20}{l}} {{F^*} ( u ) = \displaystyle\sum\limits_{m = 1}^n {\left[ {\dfrac{1}{2}\displaystyle\iiint_V {\nabla u \cdot ( {\sigma \nabla u} ) {\rm{d}}V - \iiint_V {uI{\text{δ}} ( {r_0} ) {\rm{d}}V}}} \right]} } ,\\ {\delta {F^ * } ( u ) = 0} ,\\ {u = 0}\quad{{\varGamma ^\infty }} {\text{.}} \end{array}} \right. $$ (8) 将方程(8)中第一式展开成如下的表达形式:
$$ {F^ * } ( u ) = \sum\limits_{m = 1}^n {\left\{ {\frac{1}{2}\iiint_V {\left[ {\frac{{\partial u}}{{\partial x}}\left( {\sigma \frac{{\partial u}}{{\partial x}}} \right) + \frac{{\partial u}}{{\partial y}}\left( {\sigma \frac{{\partial u}}{{\partial y}}} \right)+ \frac{{\partial u}}{{\partial {\textit{z}}}}\left( {\sigma \frac{{\partial u}}{{\partial {\textit{z}}}}} \right)} \right]{\rm{d}}V - \iiint_V {uI{\text{δ}} ( {r_0} ) {\rm{d}}V}}} \right\}} {\text{.}} $$ (9) 由于本文采用线性的四面体单元剖分地电模型,因此单元e内任意点p(
$x, y, {\textit{z}}$ )的电位u可表示为$$ {u_p} = \sum\limits_{i = 1}^4 {{N_i}{u_i}} \text{,} $$ (10) 式中Ni和ui分别为四面体单元第i个节点的插值函数和电位值。
将式(10)代入到式(9)中,并取F*对ui的偏导可得
$$ \frac{{\partial {F^ * }}}{{\partial {u_i}}} = \sum\limits_{m = 1}^n {\sum\limits_{j = 1}^4 {\iiint_V {\left[ {\sigma {u_j}\left( {\frac{{\partial {N_i}}}{{\partial x}}\frac{{\partial {N_j}}}{{\partial x}} + \frac{{\partial {N_i}}}{{\partial y}}\frac{{\partial {N_j}}}{{\partial y}} + \frac{{\partial {N_i}}}{{\partial {\textit{z}}}}\frac{{\partial {N_j}}}{{\partial {\textit{z}}}}} \right) - {N_i}I{\text{δ}} ( {r_0} ) } \right]{\rm{d}}V}} } {\text{.}} $$ (11) 对单元e进行单元分析时,单元刚度矩阵Be中元素
${{B}}^e_{ij}$ 为$$ {{B}}_{ij}^e = \iiint_V {\sigma {u_i}\left( {\frac{{\partial {N_i}}}{{\partial x}}\frac{{\partial {N_j}}}{{\partial x}} + \frac{{\partial {N_i}}}{{\partial y}}\frac{{\partial {N_j}}}{{\partial y}} + \frac{{\partial {N_i}}}{{\partial {\textit{z}}}}\frac{{\partial {N_j}}}{{\partial {\textit{z}}}}} \right){\rm{d}}V} \text{,} $$ (12) 式中三重积分采用如下的公式进行计算:
$$ \iiint_V {{{ ( {{N_1}} ) }^k}{{ ( {{N_2}} ) }^l}{{ ( {{N_3}} ) }^b}{{ ( {{N_4}} ) }^a}{\rm{d}}V = \frac{{k!l!b!a!}}{{ ( {k + l + b + a + 3} ) !}}}6{V^e} \text{,} $$ (13) 式中V e为单元e的体积。
将单元刚度矩阵合成总体刚度矩阵,即可得到最终的大型线性方程组
$$ {\boldsymbol{BU}}{\text{ = }}{\boldsymbol{f}} \text{,} $$ (14) 式中,B为大型稀疏系数矩阵,U为待求节点电位向量,f为源项向量。最后采用PARDISO求解器求解方程(14)。
视电阻率计算公式为
$$ \rho = \frac{{\Delta u}}{I}\frac{1}{{{K_M} - {K_N}}} \text{,} $$ (15) 式中,Δu为电位差,M和N为测量电极,K为装置系数,其表达式为
$$ \begin{split}& {K_p} = \frac{1}{{4\pi ( {l_2} - {l_1} ) }}\ln \frac{{\left( {{l_2} - {{\textit{z}}_p} + \sqrt {{{\left({l_2} - {{\textit{z}}_p}\right)}^2} + r_p^2} } \right)\left( {{l_2} + {{\textit{z}}_p} + \sqrt {{{ ( {l_2} + {{\textit{z}}_p} ) }^2} + r_p^2} } \right)}}{{\left( {{l_1} - {{\textit{z}}_p} + \sqrt {{{ ( {l_1} - {{\textit{z}}_p} ) }^2} + r_p^2} } \right)\left( {{l_1} + {{\textit{z}}_p} + \sqrt {{{ ( {l_1} + {{\textit{z}}_p} ) }^2} + r_p^2} } \right)}} - \\& \frac{1}{{4\pi ( {l_{2b}} - {l_{1b}} ) }}\ln \frac{{\left( {{l_{2b}} - {{\textit{z}}_p} + \sqrt {{{ ( {l_{2b}} - {{\textit{z}}_p} ) }^2} + r_{bp}^2} } \right)\left( {{l_{2b}} + {{\textit{z}}_p} + \sqrt {{{ ( {l_{2b}} + {{\textit{z}}_p} ) }^2} + r_{bp}^2} } \right)}}{{\left( {{l_{1b}} - {{\textit{z}}_p} + \sqrt {{{ ( {l_{1b}} - {{\textit{z}}_p} ) }^2} + r_{bp}^2} } \right)\left( {{l_{1b}} + {{\textit{z}}_p} + \sqrt {{{ ( {l_{1b}} + {{\textit{z}}_p} ) }^2} + r_{bp}^2} } \right)}}, \end{split} $$ (16) 式中:p代表M或N;l1,l2,l1b和l2b分别为发射、回流线源的顶、底埋深,单位为m;zp为测点埋深,单位为m;rp和rbp分别为测点距发射、回流线源的径向距离,单位为m。
2. 模型计算及结果分析
为验证本文算法的计算精度,设计垂直线源长度为300 m、电阻率为100 Ω·m的均匀半空间模型(图1),其中回流线源位于发射线源右侧4 km处且位置固定不变(下同),电流强度为20 A (下同)。从图1b中可见,均匀半空间模型的数值解与解析解基本重合,相对误差小于1%,表明了本文算法的可靠性。由图1c可知,电流主要分布于源附件一定范围内,可见以发射源为圆心呈放射状布设测量电极(如图1a所示),有利对于捕捉微弱的有效信号。
图 1 均匀半空间模型的电极布设及计算结果(a) 电极布设;(b) 数值解与解析解对比及其相对误差;(c) 电流密度lgJ (单位为A/m2)分布Figure 1. Electrode layout and calculation results of uniform half-space model(a) Electrode layout;(b) Comparison of numerical solution with analytical solution and their relative error; (c) Distribution of current density lgJ (unit in A/m2)2.1 井地电法对高、低阻体的识别能力
模型设计如图2所示:线源长度为500 m,背景电阻率为100 Ω·m;A和B两个目标体的尺寸为200 m×200 m×40 m (按x,y,z顺序),中心埋深为170 m,中心点距源为300 m,电阻率分别为1 Ω·m和1×104 Ω·m,以发射源为中心呈左右对称分布。测线布设如图2c所示,以发射源为圆心,按环形测网方式呈放射状布置24条测线,测线近源端距源10 m,测点间距为10 m (下同),发射源位于图形中心(下同)。
图3给出了地表的视电阻率等值线图,可见高阻体和低阻体在图中均较清晰,低阻体呈现得更为清晰,高阻体在形态和水平位置的确定上相较低阻体略差。这表明,相较于高阻体,井地电法对低阻目标体的探测更加灵敏。因此可以利用井地电法圈定低阻积水采空区或隐伏矿体的分布范围,但是该方法整体上对高阻体和低阻体的边界刻画得仍然不够精确。
2.2 巷道模型的响应特征与边界识别
2.2.1 边界增强技术
在重磁资料的处理中,边界识别技术是刻画地质体边界的常用手段,其中数值类边界识别方法常利用极值或零值点来圈定地质体的边界范围。在井地电法勘探中,目标体与围岩横向边界对应地表处的电流密度变化剧烈,因而在位场资料处理中将边界识别的高阶求导和归一化总水平导数引入井地电法中,以便更加精细地刻画目标体的边界位置。
导数的求解公式由五点拉格朗日插值函数表示。一阶导数的计算公式为
$$ {f'} = L_0'{f_0} + L_1'{f_1} + L_2'{f_2} + L_3'{f_3} + L_4'{f_4} \text{,} $$ (17) 二阶导数计算公式为
$$ {f{{''}}} = L''_0{f_0} + L''_1{f_1} + L''_2{f_2} + L''_3{f_3} + L''_4{f_4} \text{,} $$ (18) 式中
$L'_i $ 和$L''_i $ (i=0,1,2,3,4)分别为拉格朗日插值函数的一阶、二阶导数,fi为相应的井地电法响应。极坐标下总水平导数的计算公式为
$$ {\rm{Thd}} = \sqrt {f_r^2 + \frac{{f_\theta ^2}}{{{r^2}}}} \text{,} $$ (19) 式中fr和fθ分别为井地电法响应沿着径向方向r和对极角θ的导数。
归一化总水平导数为
$$ {\rm{Thd}}_{{\rm{nor}}} ( {i,j} ) = \frac{{{\rm{Thd}} ( {i,j} ) }}{{\max \left[ {{\rm{Thd}} ( {i - n:i + n, j - m:j + m} ) } \right]}} \text{,} $$ (20) 式中,Thd(i,j)为点(i,j)处井地电法响应的总水平导数,m和n表示窗口大小。
2.2.2 模型设计与结果分析
巷道模型是煤田中的常见模型,研究积水巷道的响应特征对认识和防治矿井水害事故具有积极的意义。为了考察巷道的井地电法响应规律以及对积水采空区复杂边界范围圈定的有效性,这里设计了两个不同的巷道模型。其中,模型一展示不同空间展布的低阻巷道的井地电法计算结果,模型二探讨井地电法对积水采空区复杂边界的识别能力的提高。
模型一:如图4所示,线源长度为150 m,背景电阻率为100 Ω·m,巷道的中心埋深为55 m,电阻率为1 Ω·m。巷道的规模与相对位置如表1中所示。分别计算表1中四个不同走向、不同长度巷道模型的井地电法响应。同时,为了更加清晰、准确地刻画巷道边界,对测线上的视电阻率分别求取一阶和二阶导数,测线位置如图4b中黑色虚线所示。图5为巷道模型的地表视电阻率平面等值线图,图6为视电阻率一阶导数曲线图,图7为模型3的视电阻率二阶导数曲线图。
表 1 巷道规模及位置Table 1. Roadway scale and location模型 巷道规模 巷道左侧距源/m 长/m 宽/m 高/m 1 200 20 10 100 2 200 20 10 0 3 20 200 10 100 4 500 20 10 −200 由图5可见,四个不同空间展布的低阻巷道从位置到走向均得到了较为明显的反映,但其响应规律却不尽相同。沿着径向方向的巷道(图6a,b,d)的两端边界可由视电阻率的假高阻异常边界大致圈定,视电阻率低值区的长轴与巷道走向基本平行;对于垂直于径向方向的巷道(图6c),其上下两端边界均能较为准确地给出。由于源电流的影响,巷道在地表的异常显示往径向外一侧弯曲。可见,视电阻率等值线结果有助于快速确定巷道平面展布的大致位置。
虽然视电阻率平面等值线结果能够大致显示出巷道的边界位置,但是对于边界的刻画却不够清晰和准确,因而不利于积水采空区的准确定位和后期治理。为了更加清晰地圈定巷道边界,这里利用位场边界识别技术中的高阶导数方法,分别对视电阻率求取一阶和二阶导数,利用视电阻率的梯度变化特征来进一步精细刻画巷道的边界位置。在一阶导数曲线(图6)中,沿着径向方向的巷道模型1,2和4 (分别对应于图6中a,b和d)的两端边界可由一阶导数极值点(极小值或极大值)准确给出,而垂直于径向方向的巷道模型3的两侧边界仍不够清晰。对巷道模型3的响应求取二阶导数(图7)之后,二阶导数的极值点与巷道模型3两侧边界位置则对应得比较准确。由此可见,当视电阻率不足以准确地圈定目标体边界时,可以考虑使用其一阶或二阶导数来提高井地电法对巷道边界的刻画能力。需要注意的是,视电阻率的导数尽管能比较准确地给出巷道在x方向的边界(如图6和图7),但是要确定细长巷道的整体边界是困难的。
模型二:如图8所示,垮塌巷道模型的线源长度为100 m,背景电阻率为100 Ω·m;巷道尺寸为20 m×200 m×10 m (按x,y,z顺序),其中心埋深为35 m,中心点距源210 m,电阻率为1 Ω·m,垮塌部分对称分布于巷道中心的两侧且尺寸均为20 m×40 m×10 m。观测参数中除了将测线增至48条外,其余参数均不变。图9为模型响应的平面等值线图。
图 9 平面等值线图(a) 电位;(b) 视电阻率;(c) 异常电位;(d) 电位总水平导数的负对数;(e) 电位的归一化总水平导数Figure 9. Plane contour maps(a) Electric potential;(b) Apparent resistivity;(c) Abnormal electric potential;(d) Negative logarithm of the total horizontal derivative of the electric potential;(e) Normalized total horizontal derivative of the electric potential在电位等值线图(图9a)中,低阻巷道产生的异常电位被淹没在背景场中,未能引起井地电位的显著畸变;在视电阻率等值线图(图9b)中巷道的位置指示准确、形态显示清晰;在异常电位(图9c)与总水平导数(图9d)结果中同样能准确地给出巷道的中心位置,但均不能准确刻画巷道边界的范围;在电位归一化总水平导数图(图9e)中,垮塌巷道的边界变化特征显示清晰,水平展布圈定准确。由此可见,将重磁位场处理边界的方法应用到井地电法中不仅是合适的,而且采用的电位归一化总水平导数方法显著增强了井地电法刻画巷道边界的效果。
2.3 地形的井地电法响应规律及其校正技术
我国约三分之二的国土面积为山区丘陵,众多资源分布于地形起伏地区,实际勘探常受到各种复杂地形环境的影响。张天伦和张伯林(1995)的研究结果显示,存在地形时可能会淹没目标体的井地电法响应,以至于无法准确地确定目标体的边界范围。因此,有必要深入研究地形对井地电场分布的影响规律以及地形校正技术。
为探究地形的响应特征,设计了如图10所示的均匀半空间表面的金字塔形状的纯地形(地垒或地堑)模型,其中:线源长度为500 m,电阻率为100 Ω·m;金字塔地形的宽边长260 m,窄边长100 m,中心点距源300 m。计算地垒高或地堑深取为20 m和40 m时模型的井地电法响应。图11为地表视电阻率等值线图,图中白色虚线框表示金字塔地形的上、下两个截面边界;图12为视电阻率曲线图,测线位置如图10b中虚线所示。
由图11和图12可见,地形对井地电场分布的影响是显著的,其影响范围基本集中在地形分布区域。图11中视电阻率的异常值能准确地刻画出金字塔地形的截面边界和各条棱边。图12中视电阻率的极值也准确地指示了地形的截面边界位置,地垒和地堑地形的响应特征正好与地形形状相反,具体表现为地垒低阻、地堑高阻的异常现象,且地垒地形越高或地堑地形越深,对井地视电阻率响应的影响则越大。
由上述计算结果可见,地形对井地电位响应的影响显著,对地形起伏变化(特别是剧烈变化)地区的电法资料进行精细化处理,欲获得准确、可靠的地下电性结构信息,就必须消除地形的影响。为了有效地削弱地形因素对井地电法高精度成像的干扰,本文将差异场处理技术用于井地电法的地形压制中。该方法的具体处理思路如下:地形产生的异常电位值Utopo由含地形的均匀半空间电位Ubackg与不含地形的均匀半空间电位Uhasp来表达,即
$$ {U_{{\rm{topo}}}} = {U_{{\rm{backg}}}} - {U_{{\rm{hasp}}}} \text{;} $$ (21) 再从含地形模型的总场电位Utotal中减去地形产生的异常电位Utopo,即可实现地形校正处理。地形校正后的新总场电位Unew为
$$ {U_{{\rm{new}}}} = {U_{{\rm{total}}}} - {U_{{\rm{topo}}}} {\text{.}} $$ (22) 为了验证上述差异场地形校正技术的应用效果,设计一个含地形的三维井地电法模型,如图13所示。模型中线源长度为500 m,背景电阻率为100 Ω·m;代表积水采空区的低阻目标体的尺寸为200 m×200 m×20 m (按x,y,z顺序),中心埋深为160 m,中心点距源300 m,电阻率为1 Ω·m。地堑和地垒地形均设置为金字塔形状,具体参数如图13中所示。图14为未进行地形校正的地表视电阻率等值线图。
从图14中能明显地看出地形与低阻体共同引起的井地电法响应,该响应不仅体现了地下目标异常体的响应,也叠加了纯地形的结果,而且由于地形在浅地表面,其响应的影响幅度甚至远大于低阻体。也就是说,由于地形的干扰,地下目标体的响应甚至可能被淹没,这使得准确圈定低阻体的分布范围变得困难。另一方面也表明地形对井地电法响应产生的影响显著。为了提高井地电法在地形起伏地区的应用效果,有必要对地形起伏地区的观测数据首先进行地形校正处理。
图15为地形校正之后的地表视电阻率等值线图,可见经过地形校正后(如图15a,b),地形对井地电法成像的影响被显著压制,可以比较准确地给出目标体的水平位置,显著地提高了井地电法的成像精度。需要注意的是,地形校正仅消除了地形影响,并未改变由于地形存在导致接收器位置上下起伏放置所造成的影响(不同于水平地表时全部接收器到异常体的垂直距离均完全一致的情况),因此与完全水平地表时的测量结果略有差异,即:地垒地形校正后的视电阻率值(图15a)较水平地表(图15c)略大,地堑地形校正后的结果(图15b)较水平地表略小。通过地形校正后,目标体的位置和范围都可以得到较为准确的恢复,可见利用这种差异场处理技术已经显著地压制了地形对井地电法响应的影响。
2.4 复杂模型
为了探讨地形校正对目标体边界识别的改善效果,设计如图16所示地垒地堑组合模型:线源长度为500 m,背景电阻率为100 Ω·m;目标体的尺寸为200 m×200 m×20 m (按x,y,z顺序),中心埋深为160 m,中心点距源为430 m,电阻率为1 Ω·m,地垒高或地堑深20 m,地形其余参数如图16中所示。图17为该模型的地表视电阻率等值线图。
由图17可知,经过地形校正后显著增强了地下目标体的呈现效果:地形校正前(图17a)视电阻率主要体现地形的轮廓,而地形校正后(图17b)视电阻率低阻异常准确地给出了目标体的水平位置;在视电阻率归一化总水平导数结果(图17c)中,异常主要体现的是地形的棱边,不能指示目标体的边界位置,而地形校正后的视电阻率归一化导数结果(图17d)较好地呈现了目标体的边界位置。可见,消除地形的影响有助于提高井地电法对目标体边界的识别能力。
3. 讨论与结论
本文基于非结构网格有限元方法实现了三维井地模型的电位响应模拟。文中分别计算了巷道模型和带地形模型的电位响应,并将位场处理方法引入到井地电法资料处理中,利用井地电法响应的导数更准确地圈定了目标体的边界范围。同时本文将差异场处理技术用于压制地形的影响,通过求取和分离地形引起的电位响应,显著增强了井地电法对复杂模型的成像效果。通过对设计的一系列三维地电模型的计算,得到了以下几点认识:
1) 井地电法响应的变化率能显著地增强井地电法对地下目标体边界位置的识别能力,将其用于刻画目标体的边界位置能够提高井地电法的应用效果和资料解释准确度。
2) 差异场地形校正技术能有效地压制地形干扰,利用其消除地形影响对于提高地形起伏地区井地电法资料处理与解释的精度和可信度是至关重要的。
-
图 1 均匀半空间模型的电极布设及计算结果
(a) 电极布设;(b) 数值解与解析解对比及其相对误差;(c) 电流密度lgJ (单位为A/m2)分布
Figure 1. Electrode layout and calculation results of uniform half-space model
(a) Electrode layout;(b) Comparison of numerical solution with analytical solution and their relative error; (c) Distribution of current density lgJ (unit in A/m2)
图 9 平面等值线图
(a) 电位;(b) 视电阻率;(c) 异常电位;(d) 电位总水平导数的负对数;(e) 电位的归一化总水平导数
Figure 9. Plane contour maps
(a) Electric potential;(b) Apparent resistivity;(c) Abnormal electric potential;(d) Negative logarithm of the total horizontal derivative of the electric potential;(e) Normalized total horizontal derivative of the electric potential
表 1 巷道规模及位置
Table 1 Roadway scale and location
模型 巷道规模 巷道左侧距源/m 长/m 宽/m 高/m 1 200 20 10 100 2 200 20 10 0 3 20 200 10 100 4 500 20 10 −200 -
戴前伟,陈德鹏,陈勇雄,侯智超. 2013. 电法勘探中异常响应特征的增强算法及其实现[J]. 煤田地质与勘探,41(3):75–78. doi: 10.3969/j.issn.1001-1986.2013.03.018 Dai Q W,Chen D P,Chen Y X,Hou Z C. 2013. The enhanced algorithms and its implementation for the abnormal response characteristics in electrical exploration[J]. Coal Geology &Exploration,41(3):75–78 (in Chinese).
李长伟,阮百尧,吕玉增,段长生,杨庭伟,王建历. 2010. 点源场井-地电位测量三维有限元模拟[J]. 地球物理学报,53(3):729–736. Li C W,Ruan B Y,Lü Y Z,Duan C S,Yang T W,Wang J L. 2010. Three-dimensional FEM modeling of point source borehole-ground electrical potential measurements[J]. Chinese Journal of Geophysics,53(3):729–736 (in Chinese).
谭河清,沈金松,周超,董辉,房锡业,张福莱. 2004. 井-地电位成像技术及其在孤东八区剩余油分布研究中的应用[J]. 石油大学学报(自然科学版),28(2):31–37. Tan H Q,Shen J S,Zhou C,Dong H,Fang X Y,Zhang F L. 2004. Borehole-to-surface electrical imaging technique and its application to residual oil distribution analysis of the eighth section in Gudong Oilfield[J]. Journal of the University of Petroleum,China (Edition of Natural Science)
,28(2):31–37 (in Chinese). 汤井田,张继锋,冯兵,林家勇,刘长生. 2007. 井地电阻率法歧离率确定高阻油气藏边界[J]. 地球物理学报,50(3):926–931. doi: 10.3321/j.issn:0001-5733.2007.03.035 Tang J T,Zhang J F,Feng B,Lin J Y,Liu C S. 2007. Determination of borders for resistive oil and gas reservoirs by deviation rate using the hole-to-surface resistivity method[J]. Chinese Journal of Geophysics,50(3):926–931 (in Chinese).
王智,潘和平. 2014. 三维井地电阻率法异常响应特征增强算法模拟研究[J]. 石油物探,53(4):491–500. doi: 10.3969/j.issn.1000-1441.2014.04.016 Wang Z,Pan H P. 2014. Research on the enhanced algorithms of the abnormal response characteristics for 3D borehole-to-surface resistivity method[J]. Geophysical Prospecting for Petroleum,53(4):491–500 (in Chinese).
王智,吴爱平,李刚. 2018. 起伏地表条件下的井中激电井地观测正演模拟研究[J]. 石油物探,57(6):927–935. doi: 10.3969/j.issn.1000-1441.2018.06.015 Wang Z,Wu A P,Li G. 2018. Forward modeling of borehole-ground induced polarization method under undulating topography[J]. Geophysical Prospecting for Petroleum,57(6):927–935 (in Chinese).
薛国强,闫述,陈卫营. 2016. 电磁测深数据地形影响的快速校正[J]. 地球物理学报,59(12):4408–4413. doi: 10.6038/cjg20161202 Xue G Q,Yan S,Chen W Y. 2016. A fast topographic correction method for electromagnetic data[J]. Chinese Journal of Geophysics,59(12):4408–4413 (in Chinese).
杨沁润,谭茂金,李桂山,张福莱,白泽. 2020. 大斜度井和水平井井地三维电阻率数值模拟和联合反演[J]. 地球物理学报,63(12):4540–4552. doi: 10.6038/cjg2020O0137 Yang Q R,Tan M J,Li G S,Zhang F L,Bai Z. 2020. Numerical simulation and joint inversion of three-dimensional borehole-to-surface resistivity of high deviated or horizontal wells[J]. Chinese Journal of Geophysics,63(12):4540–4552 (in Chinese).
张天伦,张伯林. 1995. 消除直流电阻率三极梯度法中各种干扰的实验与研究[J]. 石油地球物理勘探,30(1):100–110. doi: 10.13810/j.cnki.issn.1000-7210.1995.01.013 Zhang T L,Zhang B L. 1995. Research on removing noises in DC resistivity trielectrode gradient measurement[J]. Oil Geophysical Prospecting,30(1):100–110 (in Chinese).
Kong F N,Johnstad S E,Røsten T,Westerdahl H. 2008. A 2.5D finite-element-modeling difference method for marine CSEM modeling in stratified anisotropic media[J]. Geophysics,73(1):F9–F19. doi: 10.1190/1.2819691
Ku C C,Hsieh M S,Lim S H. 1973. The topographic effect in electromagnetic fields[J]. Can J Earth Sci,10(5):645–656. doi: 10.1139/e73-065
Li L L,Han L G,Huang D N. 2014. Normalized edge detection,and the horizontal extent and depth of geophysical anomalies[J]. Appl Geophys,11(2):149–157. doi: 10.1007/s11770-014-0436-2
Spitzer K. 1995. A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods[J]. Geophys J Int,123(3):903–914. doi: 10.1111/j.1365-246X.1995.tb06897.x
Xiong Z T,Tang X G,Li D D,Zhang L Q. 2019. Linear source CSAMT response simulation in the 2D anisotropic formation with topography[J]. J Appl Geophys,171:103861. doi: 10.1016/j.jappgeo.2019.103861
Xiong Z T,Tang X G,Qiu W Z,Zhao C Y,Zhang L Q. 2020. New algorithm for three-dimensional borehole-to-surface apparent resistivity imaging based on unstructured mesh finite-element method[J]. IEEE Access,8:184–194. doi: 10.1109/ACCESS.2019.2961806
Yang Q R,Tan M J,Zhang F L,Bai Z. 2021. Wireline logs constraint borehole-to-surface resistivity inversion method and water injection monitoring analysis[J]. Pure Appl Geophys,178(3):939–957. doi: 10.1007/s00024-021-02674-6