直流电透视三维反演技术在探测含水陷落柱中的应用研究

张 硕，赵文龙

(中铁二院工程集团有限责任公司，成都 610031)

0 引言

我国工业和国民经济的持续高速发展，都离不开能源产业和矿产资源强有力的支撑。而煤炭资源的开采已是国家获取能源矿产的主要途径，随着开采速度的加快，其开采深度和开采规模也不断扩大，矿井突水、淹井事故已经成为煤矿的重大灾害之一，给各大井田的安全生产带来威胁。据统计国内80%的矿井突水事故是由于构造引起的[1-2]。陷落柱一旦成为地下水的通道，将会对矿井的安全生产造成巨大威胁，特别是陷落柱可以导通多层含水层，使整体的水力联系增强，可将煤层底板下承压的奥灰水导入矿井，造成重大灾害[3-4]。

近年来，直流电法技术蓬勃发展，在探测巷道顶底板隔水层厚度，含水、导水构造，断裂破碎带以及工作面内部小构造等都取得了显著的效果，但在井下探测采煤工作面底板围岩中隐伏的含水，导水构造方面仍是一个难题[5]。由于直流电的穿透特性及高阻煤层和巷道空间对电流的屏蔽，矿井直流电透视技术对于煤层底板内的隐伏导水构造具有良好的响应规律。目前，已有几种典型的建场方法，二维和三维条件的正演计算也颇有成果，资料处理的方法主要为曲线对比法和地电成像法，相应的反演成像研究正在开展[6-7]，但相对还很欠缺。随着计算机硬件的提升，直流电阻率法的三维反演技术不断发展，在起伏地形条件，约束性反演等方面取得了较大进步[8-10]。为了进一步丰富直流电透视资料解释方法，笔者以直流电法三维正反演程序与直流电透视方法相结合，采用巷道间AM法，直接测量电位，对隐伏含水陷落柱的探测进行了详细研究。

1 正演边界条件及数值模拟

1.1 正演边界条件

在稳定电流场中，由欧姆定律的微分形式出发，根据电荷守恒定理，高斯定理等可以得出点电源三维电场分布方程：

▽·(σ(x,y,z)▽u(x,y,z))=f

(1)

即：

(2)

因此正演要解决的问题，即是求解式(2)在满足一定边界条件下的解，笔者采用混合边界条件，Dey等提出了混合边界条件[14]：

(3)

其中：n为边界Γ∞外法线方向上的坐标变量；r为电流源到边界上的距离；θ为n与r的夹角。该方法可以保持电位在边界上的物理特性，可以减少边界附近的网格数量，减少计算量。

1.2 模型设置

为了研究矿井直流电透视法对煤层底板陷落柱的响应特征，采用等比原则，设置如图1所示地电模型，采用有限差分法、采用A-M观测方式进行正演计算：如图1(a)所示为数值模型的三维结构示意图，从上到下分别为顶板围岩、煤层、底板围岩，图中蓝色方形体为拟设定的模拟导水陷落柱。图1(b)为煤层底板的水平切片，图1(c)为图1(a)沿红色线框切下的断面图。整个计算区域为190 m×130 m×124 m的长方体区域(网格为38×26×25)，整个区域存在两种网格剖分形式，在煤层内网格剖分为5 m×5 m×4 m，煤层厚度为4 m；其他区域网格剖分大小为5 m×5 m×5 m。煤层电阻率值设为1 000 Ω·m，顶底板电阻率值设为100 Ω·m ，巷道的电阻率设为3 000 Ω·m ，回采工作面两侧巷道间距设为100 m。陷落柱的大小为20 m×20 m×30 m，电阻率值为10 Ω·m ，异常体模型的顶面位于煤层底板下深0 m处，在x方向上位于30个电极的中心，在y方向上位于两巷道的中间位置。

图1 数值模型参数设置说明图Fig.1 Numerical model parameter setting instructions(a)装置模型示意图;(b)煤层底板切片;(c)YOZ断面

模型中电极分别布设于两侧巷道的底板上，供电电极A布设30个，分别编号为A1、A2…、A30，点距为5 m，测量电极同样布设30个，分别编号为M1、M2…、M30，测量点距为5 m。当A1点供电时，M1到M30号测点分别进行观测，同样当A2点供电时，M1到M30号测点分别进行观测，这样直到A30点供电完成，一次正演共测得900个数据。

1.3 数值模拟结果分析

选取A1、A10、A15供电时测量的数据绘制成电阻率曲线图(图2)。当A1电极供电时，从M1号电极到M30号电极所测得视电阻率值逐渐增大，最大值达到105 Ω·m 以上；当A10电极供电时，从M1号电极到M30号电极所测得视电阻率值逐渐增大，但整体值较A1号电极供电时要小；当A15号电极供电时，测得视电阻率曲线关于M15号测点对称。整体计算视电阻率值均大于底板围岩电阻率，即底板存在低阻异常体时二极装置测得视电阻率值出现反异常现象[15-16]。

图2 A1、A5、A15供电视电阻率曲线图Fig.2 Apparent resistivity curve of A1、A5、A15 Power Supply(a)观测模型;(b)视电阻率曲线

2 反演成像理论

笔者采用拟高斯-牛顿法进行反演，运用共轭梯度算法求取模型修改量，提高计算速度[8]。同时，为了减少内存的消耗，程序采用不直接求取雅克比矩阵J的计算方法[9]。

2.1 反演方程

设参数模型为m=(ρ1,ρ2,...,ρk)T，其中k为网格数量，观测数据dobs=(ρs1,ρs2,...,ρsn)T，其中n为观测数据量，模型正演计算数据为dcal。首先构造反演目标函数：

(4)

引入拉格朗日常数λ和光滑度矩阵C，得到：

λ(CΔm)T(CΔm)

(5)

如果函数φ(m)的值趋于极小，则函数φ(m)的导数等于零：

(6)

因为dcal(mk+1)=dcal(mk)+J·(Δmk)，所以式(6)可写为：

(JTJ+λCTC)Δm=JTΔd

(7)

当完成非线性问题的线性化，求解该方程组是三维电阻率反演成像中的重要的问题。解上述反演方程的困难在于三维问题的雅克比偏导数矩阵维数很大，无法直接求取，因此采用不直接求取J的方法。

2.2 雅克比矩阵J的处理

采用不直接求取J的方法，利用一次拟正演求取雅克比偏导数矩阵J与一个向量的乘积，从而大大降低了计算机的内存占用，提高了计算速度[11-14]。

以一个3×3的矩阵进行运算，假设观测数据个数为3，模型网格数量为3，以此为例，雅克比偏导数矩阵表达式为：

(8)

所以

(9)

在正演过程中由于

(10)

因此

(11)

(12)

(13)

所以

(14)

其中

(15)

所以两边求导得

(16)

由于A为对称矩阵，所以

(17)

最终JTΔd可写成如下形式：

(18)

这样JTΔd的计算量就相当于以as1、as2、…、asK为场源的一次正演。按照这种方法，同样可以利用一次拟正演计算出J·x的值，大大降低了反演过程中计算机内存开销，缩短计算时间。

2.3 共轭梯度算法

3 数值模型实验

3.1 实验模型

建立三维空间的含水陷落柱模型，将观测区域剖分为14 m×14 m×12 m的网格，设陷落柱模型规格为2 m×2 m×2 m，电阻率值为10 Ω·m ，底板围岩电阻率为100 Ω·m ，覆盖煤层电阻率为2 000 Ω·m 。陷落柱的中心位置坐标为(0 m，0 m，-1 m)，两平行巷道间距离为9 m，左侧为一连通巷道，空间位置如图3所示，图3(b)和图3(c)分别为导水陷落柱在xoy平面和zoy平面上的投影。

图3 含水陷落柱及电极布设模型Fig.3 The numerical model of collapsed column and electrode location(a)三维地电模型;(b)沿xoy切片；(c)沿zoy切片

3.2 反演结果分析

反演中取模型光滑因子为0.05，经过五次迭代得到结果进行切片，如图4所示。蓝色低阻区域电阻率小于60 Ω·m ，呈椭球状，反演结果可以很好地反映含水体的位置。

图4 反演电阻率结果Fig.4 Inversion resistivity results

4 物理实验及结果分析

对上述的数值模拟结果进行物理实验的验证，采用拟高斯-牛顿法对物理实验数据反演。由于高阻煤层和巷道空间对电流的屏蔽作用，在巷道底板布设电极的电场分布特征与地面半空间中布设电极的电场分布特征极为一致，所以在物理模拟中采用了地面半空间模型，采用一定浓度的含食盐泥浆模拟实际的含水陷落柱，如图5所示,其中观测电极为A1、A2、…、A30号，相邻供电点距离为0.1 m，测量电极命名为M1、M2、…、M30号，相邻测点距离为0.1 m，陷落柱大小为40 cm×40 cm×40 cm，电阻率值约5 Ω·m 。

图5 物理模拟实验设计与现场图Fig.5 Design and site of physical simulation experiment(a)物理模拟实验示意图;(b)实验现场

4.1 实测数据分析

在拟定异常区域开挖之前进行一次观测，得出不存在导水陷落柱时的观测值，然后根据实验设计，按照设定的尺寸开挖、充水，再次进行测量，得到两组视电阻率曲线的对比图。

实验数据提取M点测得的电位值，采用二极(AM)装置系数计算电阻率值。为了较为直观地对实验数据进行分析，选取A15号供电电极供电时由30个测量电极测得的数据绘制电阻率曲线图，如图6所示，从M1到M30点观测时，整体视电阻率表现出先增大后减小的趋势，在M14号电极观测的视电阻率值最高。物理实验结果与数值模拟结果具有一致性。

图6 A15号电极供电视电阻率曲线图Fig.6 Apparent resistivity curve of A15 power supply(a)A15号电极供电观测;(b)视电阻率曲线

结合数值模拟结果分析，测区AM之间存在低阻异常时，中间电极A供电时所测得电阻率值反而增大，出现了反异常现象。该现象是由于低阻体在地表，对电流产生吸引，地表AM间整体电阻下降，导致测点M处测得电位升高。沈平等[15]在井间电阻率成像方法的研究中也提出了这种反异常的现象，并对此进行了解释。汤井田等[16]在研究不同装置下点源球体的近似解与精确解对比时，对A-M法出现正反异常的条件做了详细的论述，认为极距AM应小于异常体的埋深与半径之比，否则会出现反异常的现象。

4.2 实测数据反演

我们用以上理论对上述物理实验数据进行反演成像处理，反演初始模型电阻率值设为22 Ω·m ，模型网格为35×35×30，经过5次迭代用时42 s，拟合误差下降到5%。同一模型参数，采用直接求取雅克比偏导数矩阵J的方式进行计算迭代一次的时间为765 s。利用一次拟正演计算雅克比矩阵和一个向量乘积的方法，可以大大缩短计算时间。

物理模型实验反演结果如图7所示，实际模拟陷落柱在x方向和y方向坐标为(0 m，0 m)，中心深度为0.2 m。由反演结果可以看出，在探测区域的中心位置(实际模拟含水陷落柱所在位置)，出现了一个明显的低阻异常区域。通过z方向上切片可以看出，低阻异常区的中心深度为0.6 m，模拟低阻异常区域实际挖掘的深度为0.2 m，考虑到在实验进行中泥浆水不断下渗，所以实际低阻异常的中心位置下移。该区域在水平位置垂向位置上都和模拟低阻体所在位置具有一致性，反演结果能较好地反映模拟低阻异常体所在三维空间中的位置。从反演结果中可以看出存在一处明显的高阻异常区域，由实际实验区域电阻率分布不均造成。

图7 实测数据反演结果切片图Fig.7 The map of 3D inversion used measured data

5 结论与建议

通过对陷落柱模型的数值模拟和物理模拟，以及实验结果的三维反演成像研究得出以下结论：

1)通过对导水陷落柱模型的正演计算结果可以看出，采用常规视电阻率计算公式进行计算时，在低阻陷落柱区域计算结果为高阻，出现反异常现象，测点距离陷落柱越远，视电阻率越接近围岩电阻率值。

2)在进行三维反演时，采用不直接求取雅克比矩阵J的方法，利用一次拟正演求取雅克比矩阵J与一个向量的乘积，从而大大降低了计算机的内存占用，加快了计算速度。

3)对物理实验数据进行三维反演，结果很好地反映出了导水陷落柱的实际位置，数值模拟的视电阻率计算结果和实测数据的反演结果，为实际电透视探测陷落柱的解释提供了借鉴。

4)以后将继续针对不同的观测系统，不同空间位置的异常体进行研究，进一步总结探测规律。