起伏地形下井-地联合电阻率观测含水陷落柱模型异常特征

孟 麟， 王 智

(长江大学电子信息学院， 荆州 434023)

在实际煤炭的探勘和开采过程中，陷落柱是不可避免会遇到的一种灾害性的地质异常体。它是由于地下水不断对上方可溶性岩石进行溶蚀和上方覆盖岩层重力的共同作用下产生的溶岩塌陷体，给煤田的开采工作带来了极大的安全隐患。为了预防安全隐患，准确勘探出陷落柱位置是最有效的方法，许多学者利用陷落柱和周围围岩之间的密度、磁性和电性差异，在地表采用地震勘探[1]、瞬变电磁法[2]、直流电法[3]和高密度电法[4]等方法提前进行陷落柱的勘探。受地表条件限制和陷落柱空间形态特征的影响，目前准确高效地探测陷落柱位置信息还存在一定的困难。特别是对于直流电法勘探来说，地表条件会对勘探的深度产生比较大的影响，且直流电法勘探的精度会随着地质体埋深的增加而降低。为了解决勘探深度的问题，井-地联合电阻率观测装置成为了较好的选择。相比与井-井(或地-地)观测装置，井-地联合观测装置对于地下深部的勘探有着异常幅值和勘探精度都更高的优点，在勘探过程中有着很好的应用前景。近年来，井-地联合观测装置已在勘察领域得到广泛应用。汤井田等[5]利用井-地电阻率法得到歧离率来确定高阻率油藏边界；赵广茂等[6]用井-地观测装置来勘探地下陷落柱的存在和含水性；王智等[7]利用井-地观测方法和改进的异常电位法，确定了适用于起伏地表的自然边界条件。

虽然井-地联合电阻率观测装置已经被广泛应用在实际工程勘探中，理论研究方法也存在许多，但是对于在起伏地表的情况下不同形态和位置的典型异常体的异常形态特征的系统研究不足。使用井-地联合电阻率观测装置针对起伏地表情况下，含水陷落柱处于不同形态和位置所产生异常情况进行系统的分析，以期为实际勘探含水陷落柱提供理论依据。

1 2.5D有限单元法数值模拟

在三维无限半空间Ω中，建立直角坐标系。假设Z轴垂直向下，那么点电源O的坐标为(xO,yO,zO),表示为O，则总电位满足微分方程为[8]

(1)

(2)

式(2)中：u(x,y,z)简写为u；ΓS和Γ分别为三维无线半空间Ω的地表边界和无穷远边界；n为地表边界ΓS的外法向向量；r表示该坐标点到点电源之间的距离。

假设电性的分布沿Y轴方向不变，对式(1)、式(2)进行傅氏变换，即将空间域的三维边界问题变成波数域的二维边界问题[8]：

(3)

式(3)中：U(x,ki,z)表示波数域的电位，简写为U，其中ki为波数(i=1,2,…,n)；K1和K0分别为一阶和零阶贝塞尔函数。式(3)对应的变分问题[8]为

δF(U)=0

(4)

采用不规则三角形单元e对求解区域进行剖分后，Fe(U′)表示每个三角形单元上波数域电位U′的函数，并将其进行总体合成后得

(5)

KU=P

(6)

解该方程组，便可求得节点m在对应波数k中的电位U(m,ki)。然后按式(7)进行傅里叶反变换便可得到空间中的电位u(m)：

(7)

式(7)中：波数(ki)与傅里叶反变换系数(gi)选自徐世浙[8]采用最优化方法得到的5个波数和对应的傅里叶反变化系数。

2 观测方式与观测装置

井中和地表电极的联合观测方法可以分为：井-地联合观测方式和地-井联合观测方式。井-地联合观测方式是指将供电电极置于井中，测量电极置于地表的观测方式；地-井联合观测方式是指将供电电极置于地表，测量电极置于井中的观测方式。两种观测方式均可以通过分别移动供电电极和测量电极和相对位置来研究不同水平位置和不同深度地下介质的视电阻率分布情况。再根据在实际勘探工作中所使用的不同观测装置，可以将其细分为井-地二极装置、地-井二极装置、井-地三极装置、地-井三极装置、井-地四极装置、地-井四极装置等。因篇幅有限，研究选择井-地二极装置为观测装置。

井-地二极测深装置，指供电电极A置于井中，另一供电电极B置于无穷远处；观测电极M置于地表，另一观测电极N置于无穷远处。再通过上下移动供电电极A和水平移动观测电极M，便可以探测不同深度和水平位置的地下介质分布。以M的横坐标为观测记录点的横坐标，A的纵坐标为观测记录点的纵坐标画出测深断面图，观测示意图如图1所示。

图1 井-地二极观测装置示意[9]Fig.1 Sketch map of pole-pole on borehole-ground joint resistivity surveying system[9]

3 数值算例及分析

3.1 算法验证

在完成有限单元法中网格剖分的部分时，使用Gmsh软件[10]来建立典型的地质模型并完成局部加密的非结构化网格剖分。在存储大型稀疏矩阵时使用CSR(compressed sparse row)格式[11]，减少程序对内存的占用。对于求解正演过程中的线性方程组，选用Eigen库中的BiCGSTAB(biconjugate gradient stabilized method)求解器。在计算多个点电源时，使用MPI(message passing interface)进行并行计算，提高计算效率。为了验证算法的正确性，首先选择具有解析解的简单模型进行计算。

二层介质的地电断面模型如图2所示。图2中，ρ1、ρ2分别表示第一层与第二层介质的电阻率，第一层断面的高度为15 m。表2给出了单个点电源(O)时，地表不同测点(V)的电位。

由表1数据表明，计算解的最大相对误差发生在距离点源最近的测量点上达到5%，其余测点的相对误差均在1%左右，证明了编译算法的正确性。

3.2 算例

使用井-地二级装置进行观测。将钻井井口与地质体重心连线在地表上的投影，记为X测线。井中供电电极A′和地表观测电极M′的点距均设置为1 m。测量电极在地表X测线0～50 m移动。钻井井口在X测线14 m处垂直向下，记为Z测线，供电电极在Z测线范围5～20 m内移动。

图2 二层介质地电断面Fig.2 Two-layers geo-electric cross section

陷落柱是煤层下伏可溶性岩石被地下水溶蚀后引起上覆岩层冒落而形成的古岩溶塌陷体。含水陷落柱在电阻率特性上显示低阻特性，为了系统地对比各种装置和不同位置模型的观测视电阻率异常，围岩电阻率ρ1为100 Ω·m，含水陷落柱电阻率ρ2为5 Ω·m。模拟的含水陷落柱的模型类似为“子弹形”，顶部埋深为8 m，底部深度16 m。通过改变模型中心距离点源所在的钻井的垂直距离来分析其对观测视电阻率异常的影响，设定的距离分别为11、21 m。若含水陷落柱模型顶点的位置确定时，改变模型的倾斜角度，分别为左倾30°、直立90°和右倾30°，对比倾斜角度的不同对观测视电阻率异常的影响。地质模型的网格剖分采用局部加密的方式，平坦地形下直立含水陷落柱模型及网格剖分如图3所示。

3.2.1 平坦地形

由图4可知，当在平坦地形的情况下含水陷落柱模型的视电阻率断面图中异常只存在独立的椭圆形的低阻异常，存在这样的低阻异常主要是因为含水陷落柱模型相对围岩电阻率为低阻，在电场中吸引电流。在X测线方向，不论距离钻井的距离和倾斜角度如何变化，低阻异常中心始终位于含水陷落柱模型的右侧；在Z测线方向，随着测点深度的增加，在含水陷落柱模型的上部分的深度区间中，在水平方向上的等值线相对于模型下部分的深度区间中的等值线更加密集，在垂直方向等值线的变化较为平缓。

表1 算法所得的电位计算解与解析解

图3 平坦地形下直立含水陷落柱模型示意Fig.3 The upright water-bearing collapse column model in flat terrain

图4 含水陷落柱模型视电阻率断面(平坦地形)Fig.4 The resistivity section of water-bearing collapse column model(flat terrain)

对比图4(a)、图4(c)、图4(e)和图4(b)、图4(d)、图4(f)可以看出，低阻异常的最大值和距离钻井的距离有关，在含水陷落柱模型距离钻井11 m的位置时，低阻异常的最大值为13%左右；在含水陷落柱模型距离钻井21 m位置时，低阻异常的最大值减小到10%左右。同时模型在距离钻井11 m的位置时，低阻异常中心的深度与含水陷落柱模型顶端的深度一致，水平方向上差距2～4 m，当距离钻井的距离为21 m时，低阻异常的中心发生上移，但在水平方向上的差距基本不变化。

在含水陷落柱模型距离钻井的位置发生改变时，改变模型倾斜角度所产生影响的趋势基本不发生变化。当模型左倾30°时，低阻异常被向右上方挤压；当模型右倾30°时，低阻异常被向左下方拉伸。但无论如何改变倾斜角度，在含水陷落柱模型底部向下方向的视电阻率等值线是最为稀疏的。

图5 纯山谷地形视电阻率断面Fig.5 The resistivity section of valley terrain

3.2.2 山谷地形

由图5可见，纯山谷地形的视电阻率断面图会产生高阻、低阻异常并存的情况，均呈椭圆形分布，与前人研究得出的地形与视电阻率呈左右对称的规律相符，这是由于在山谷地形中存在空气，阻碍了电流的通过。低阻异常位于山谷地形的右下方，当山谷地形临近钻井时，低阻异常中心距离山谷地形更近，位于X=28～30 m，Z=5 m处，异常最大值为22%；若山谷地形远离钻井时，低阻异常中心向下偏移，位于X=38～40 m，Z=6 m处，异常最大值为14%。高阻异常的中心位置在山谷地形的左下方，临近钻井时，其中心位于X=14～15 m，Z=5～6 m处，异常最大值为14%；远离钻井时，其中心略微向右偏移，位于X=25～26 m，Z=5 m处，异常最大值为11%。在图5中等值线最密集的地方都位于山谷地形底端的中心处。

根据图6(a)、图6(c)、图6(e)和图5(a)对比，即在临近钻井的纯山谷地形下方存在含水陷落柱模型。视电阻率断面图中依旧是高阻和低阻异常并存的情况，改变倾斜角度对高阻和低阻异常的最大值并无影响。对比纯山谷地形，此时低阻异常的最大值为36%，高阻异常的最大值改变为10%。当含水陷落柱模型左倾30°时，模型位于山谷右部低阻异常下方，视电阻率异常基本无变化；当含水陷落柱模型直立时，由于模型的低阻特性吸引电流，是视电阻率等值线向高阻异常方向挤压。模型的上部分面积小，相对应地对电流的吸引力小，等值线变化小，下部分面积相对大，对等值线影响增大；当含水陷落柱模型右倾30°时，断面图视电阻率异常不再显示为高阻异常和低阻异常左右对称，在模型的影响下阻断了高阻异常对下方的影响，产生了一个新的低阻异常，异常中心与模型的最低点平行，最大值为6%。

根据图6(b)、图6(d)、图6(f)和图5(a)对比，低阻异常的最大值变为30%，高阻异常的最大值变为12%，断面图的视电阻率等值线趋势基本无变化，是因为山谷地形对下方含水陷落柱模型产生的异常产生了遮挡效果。

根据图7(a)、图7(c)、图7(e)和图5(b)对比，即在远离钻井的纯山谷地形下方存在一个含水陷落柱模型。与纯山谷地形产生的一个低阻异常和一个高阻异常不同，此时视电阻率断面上显示两个低阻异常和一个高阻异常，从左到右依次为低阻-高阻-低阻，分别是由含水陷落柱模型阻断高阻异常产生和山谷地形产生。纯山谷地形产生的高阻和低阻异常中心位置不变动，低阻异常最大值改变为22%，高阻异常的最大值改变为2%，基本和围岩视电阻率相等。由含水陷落柱模型阻断山谷高阻异常而产生的低阻异常的中心位于模型的左侧，在Z测线方向和模型的中心平行，在X测线方向和模型中心距离5 m。由于含水陷落柱在有低阻特性，对电流有吸引作用使山谷地形中部垂直X测线的视电阻率等值线向模型方向凸起，而偏移角度基本不影响此趋势。当存在偏移角度时对由模型产生的低阻异常的形态有所影响，当右倾30°时，对低阻异常由向左下方的拉伸趋势；当左倾30°时，对低阻异常向右方拉伸趋势。

根据图7(b)、图7(d)、图7(f)和图6(b)、图6(d)、图6(f)对比可得，当远离钻井即远离测点电源的时候，含水陷落柱模型对纯山谷地形产生的影响减弱，但基本的趋势仍然存在。

在Z测线方向上，5～20 m区域内为观测区域，此区域使用Surfer软件根据数据自动填充图6 含水陷落柱模型视电阻率断面(山谷地形临近钻井)Fig.6 The resistivity section of water-bearing collapse column model(valley terrain near drilling)

在Z测线方向上，5～20 m区域内为观测区域，此区域使用Surfer软件根据数据自动填充图7 含水陷落柱模型视电阻率断面(山谷地形远离钻井)Fig.7 The resistivity section of water-bearing collapse column model(valley terrain away from drilling)

3.2.3 山脊地形

由图8可知，在纯山脊地形的视电阻率断面图中产生了高阻、低阻异常并存的情况，地形右边存在高阻异常，左边存在低阻异常与纯山谷地形产生的视电阻率异常正好相反，这也符合前人的研究结果。当山脊地形临近钻井时，低阻异常中心位于X=15 m，Z=5 m处，异常值最大为10%，高阻异常中心位于X=26～27 m，Z=5 m处，异常最大值为17%；当山脊地形远离钻井时，低阻异常中心位于X=26～27 m，Z=5 m处，异常的最大值为10%左右，高阻异常中心位于X=37～38 m，Z=10 m处，异常值最大为10%。在改变山脊地形与钻井的距离的过程中，对低阻异常造成的影响不大，但是使高阻异常发生了明显的改变，使其向下拉伸的趋势严重且产生的高阻异常最大值也下降近1/2。

根据图8(a)、图9可得，在临近钻井的纯山脊地形下方存在含水陷落柱模型时无论模型离钻井的距离和倾斜角度如何变化，在断面中都显示两个低阻异常和一个高阻异常，从左到右依次为低阻-高阻-低阻。

由图9(a)、图9(c)、图9(e)可以看出，在含水陷落柱模型距离钻进11 m时，左部低阻异常最大值为13%，同时低阻异常中心向下拉伸，在Z测线方向基本与异常体模型的中心平行，含水陷落柱模型倾斜角度为右倾30°时，低阻异常有向左拉伸的趋势。中部高阻异常最大值下降为7%，此时含水陷落柱模型将高阻异常向下的趋势阻断。右部低阻异常最大值为7%，低阻异常的中心随含水陷落柱模型的角度倾斜而改变，左倾30°时低阻异常中心向下偏移；直立状态时低阻异常中心和含水陷落柱模型的中心平行；右倾30°时低阻异常中心向上偏移。由图9(b)、图9(d)、图9(f)可以看出，当含水陷落柱模型距离钻井21 m时，左部低阻异常的最大值和形态与纯山脊地形对比都无明显区别。中部高阻异常最大值下降到12%，右部低阻异常的最大值为7%。右部低阻异常中心随含水陷落柱模型倾斜角度的变化与图9(a)、图9(c)、图9(e)显示的变化趋势相同。

根据图8(b)、图10可得，在远离钻井的纯山脊地形下方存在含水陷落柱模型的断面图与远离钻井的纯山脊地形的断面图相比，都存在一个左部的低阻异常和一个右部的高阻异常，不同之处在于存在右部的高阻异常由于受到异常体模型的影响其异常中心向下偏移严重。

由图10(a)、图10(c)、图10(e)可以看出，在含水陷落柱模型距离钻井11 m时，左部的低阻异常的最大值上升到19%，右部高阻异常的最大值下降到2%。当含水陷落柱模型右倾30°和直立时，视电阻率等值线变化不明显，但是含水陷落柱模型左倾30°时，低阻异常左部视电阻率等值线向左拉伸趋势明显。由图10(b)、图10(d)、图10(f)可以看出，在含水陷落柱模型距离钻井21 m时，左部的低阻异常的最大值上升到14%，右部高阻异常的最大值下降到3%。相比左倾30°和右倾30°，在直立状态时高阻异常中心下降的趋势最大。

通过对井-地二极装置各种地形情况和不同位置的含水陷落柱模型的正演结果分析，可以得到如下结果。

(1)在平坦地形情况下，含水陷落柱模型会产生一个低阻异常。随着模型的位置和倾斜角的改变，视电阻率断面图中的低阻异常的形态和最大值也跟着发生改变。

(2)在纯山谷或纯山脊地形的情况下，在视电阻率断面图中从左到右依次产生高阻-低阻或低阻-高阻异常。随着地形距离钻井位置改变，视电阻率断面图中显示的趋势不变，但是产生异常的形态和最大值发生改变。

在Z测线方向上，5～20 m区域内为观测区域，此区域使用Surfer软件根据数据自动填充图8 纯山脊地形视电阻率断面Fig.8 The resistivity section of ridge terrain

在Z测线方向上，5～20 m区域内为观测区域，此区域使用Surfer软件根据数据自动填充图9 含水陷落柱视电阻率断面(山脊地形临近钻井)Fig.9 The resistivity section of water-bearing collapse column model(ridge terrain near drilling)

在Z测线方向上，5～20 m区域内为观测区域，此区域使用Surfer软件根据数据自动填充图10 含水陷落柱视电阻率断面(山脊地形远离钻井)Fig.10 The resistivity section of water-bearing collapse column model(ridge terrain away from drilling)

(3)当含水陷落柱模型处于山谷地形或山脊地形的产生的低阻异常区域时，模型在视电阻率断面图中显示的异常不明显。起伏地形对其产生了遮挡作用。

(4)当含水陷落柱模型处于山谷地形或山脊地形的产生的高阻异常区域时，视电阻率断面图中产生低阻-高阻-低阻的异常情况，根据除地形产生的两个异常外的另一个低阻异常可以判断异常体的形态和位置信息。

4 结论

针对起伏地形条件下，存在含水陷落柱对地下产生的异常无系统分析的问题，采用有限单元法对井-地电极联合观测装置电阻率测深的二维地质体视电阻率断面进行数值计算。共包含平坦地形、山谷地形和山脊地形三种地形条件，在两种起伏地形中又讨论了其距离钻井不同位置的影响。每种地形条件下均设置了距离钻井11 m或21 m，倾斜角度为左倾30°、直立90°和右倾30°共6种情况的含水陷落柱模型。详细地分析了在这些不同情况下，含水陷落柱模型在视电阻率断面中产生的影响，研究结果为实际含水陷落柱勘察过程提供了一定的理论依据。