隧道超前探水的动电效应数值模拟研究

Numerical Simulation of Electrokinetic Effect of Tunnel Advance Water Exploration

摘要: 隧道超前预报是工程领域重要的前沿问题。前人使用过多种方法探测隧道掘进过程中的含水构造，但是效果一般。动电效应作为一种新的地球物理方法，兼具地震方法和电磁方法的优点，有望在隧道超前探水问题上提供新的有效的思路。本文基于有限元方法对隧道探水模型中动电效应进行数值模拟。结果表明，动电效应具有其他方法不可比拟的优势。由于本方法中弹性波场与电磁波场具有天然的耦合关系，因此可以进行相互验证，从而更加快速准确地提取波场信息；震电效应中可以通过接收排列记录的位移场判断含水构造相对位置，利用电场提取含水构造轮廓信息，再利用多点激发获得的界面响应信息进行目标体定位；电震效应更是由于位移场特殊的传播规律带给我们新的思路，从而实现仅需一次单点激发即可实现含水构造定位的目标。本文的数值模拟工作为使用动电效应进行隧道超前探水提供了施工方案和分析方法，对探测影响隧道掘进的含水构造具有借鉴价值。

Abstract: Advance forecast of tunnel is an important frontier problem in engineering field. Many methods have been used to detect water-bearing structures during the tunnel drum, but the results were not satisfactory. As a new geophysical method, electrokinetic effect has the advantages of both seismic method and electromagnetic method, which is expected to provide a new and effective idea for tunnel advanced exploration. In this paper, the electrokinetic effect of tunnel water ex-ploration model is numerically simulated based on finite element method. The results show that electrokinetic effect has incomparable advantages over other methods. Due to the natural coupl-ing relationship between elastic wave field and electromagnetic wave field in this method, mutual verification can be carried out, so that wave field information can be extracted more quickly and accurately. As for seismoelectric effect, the relative positions of water-bearing structures can be determined by receiving the displacement fields recorded in alignment, and the profile information of water-bearing structure can be extracted by electric field, and the target body can be located by using the interface response information obtained by multi-point excitation. The electroseismic effect brings us a new idea due to the special propagation law of displacement field, so as to realize the goal of water-bearing structure positioning with only a single point of excitation. The numerical simulation work in this paper provides a construction scheme and analysis method for tunnel advance water exploration using electrokinetic effect, which has reference value for detecting water-bearing structures affecting tunnel driving.

1. 引言

动电现象是当今地球物理领域的研究热点，指与流体与固体介质基体之间的相对运动有关的电磁信号和地震信号。动电效应具有地震的分辨率和电磁方法对流体的敏感性。前人针对动电效应做了大量的实验与数值模拟工作。

目前震电效应机理主要是Pride [1] 推导出的一套动电耦合方程组。Warden et al. [2] 将其推广到部分饱和情况。同时高黏土含量和油藏等特殊情况也得到补充和完善 [3] [4] [5] [6] [7]。与此同时，对动电效应耦合特性的认识也在不断深入。Ren et al. [8] [9] [10] 通过广义反透射系数法分析了水平层状模型下动电效应产生的电磁场特征，并提出了隐失波对同震电磁信号的贡献。Gao et al. [11] [12] 通过全局矩阵法和割线积分法分析了天然地震产生早期电磁信号的构成与特点。Gao et al. [13] 和高玉涛等 [14] 分别通过有限差分和有限元研究了不规则地质体模型中震电波场特征。

然而，由于对动电效应理论与耦合特征认识的欠缺，在工程上少有应用，特别是在隧道超前预报上，目前在国内报道仅见于罗金和刘春生 [15]。同样，对于隧道超前预报中动电效应研究，也缺乏相应的数值模拟工作。我们认为，动电效应的一些特性，比如震电勘探中的界面响应信号和电震勘探中的地震信号，在隧道超前预报中有着其他方法不可比拟的优势。

隧道工程具有埋深大、洞线长、地质条件复杂的特点 [16]，使得地质灾害在施工过程中频发，其中突水、突泥是主要的两种灾害。这两种灾害的发生条件和控制因素实质上是致灾构造，受地质条件和施工活动影响。一般划分为岩溶类、断层类和其他成因类 [17]。传统的物探方法在隧道超前探水上存在距离和分辨率的瓶颈，或者是仪器设备不适用于隧道这种工程环境。

本文就隧道掘进过程中常见的含导水地质构造，建立相应的地质模型，通过有限单元法来模拟由动电效应激发的地震电磁场。在分析不同地质模型引起震电信号特征基础上，对突水、突泥构造进行定位，总结地震电磁场的规律，为实际情况下动电效应隧道超前预报提供参考。

2. 理论方法

2.1. 动电耦合控制方程组

本文基于Pride动电耦合控制方程组，对二维SHTE模式进行数值模拟。考虑系数型偏微分方程形式，在频域求解，最后进行傅里叶变换将动电波场由频率域转换为时间域。

$\frac{\partial }{\partial x}\left(G\frac{\partial u_y}{\partial x}\right)+\frac{\partial }{\partial z}\left(G\frac{\partial u_y}{\partial z}\right)+\left(w^2\rho -w^2\frac{{\rho }_f^2}{\tilde{\rho }}\right)u_y+iw{\rho }_{f}L{E}_y=-F_y+\frac{{\rho }_f}{\tilde{\rho }}{f}_y$ (1)

$\frac{{\partial }^2{E}_y}{\partial x^2}+\frac{{\partial }^2{E}_y}{\partial z^2}+i{w}^3\mu {\rho }_{f}L{u}_y+\mu w^2\left(\tilde{\epsilon }+\tilde{\rho }{L}^2\right)E_y=-iw\mu L{f}_y-iw\mu C_y$ (2)

其中，u和E分别为固体介质位移和电场强度；F、f和C分别是作用于固体介质、流体介质的力载荷和电流载荷；下标y表示因变量的y向分量，w是频率； $\rho$ 、 ${\rho }_f$ 和 $\tilde{\rho }$ 分别是固体介质、流体介质密度和等效密度，L是动电耦合系数，下标为0表示其静态形式； $\mu$ 是磁导率； $\epsilon$ 和 $\tilde{\epsilon }$ 分别是介电常数和等效介电常数。相关系数详见Pride & Haartsen [18] 和Haartsen & Pride [19]。

以亥姆霍兹方程形式表示式(1)和(2)，记作

$-\frac{\partial }{\partial x}\left(d\frac{\partial B}{\partial x}\right)-\frac{\partial }{\partial z}\left(d\frac{\partial B}{\partial z}\right)+eB-f=0$ (3)

其中，d是扩散系数，e是吸收系数，f是源项，B是因变量。

SHTE边界条件分为两类，其一是空气/孔隙介质界面处，切应力等于0，电磁场水平分量连续。记作：

$n\cdot {\left(c\nabla u_y\right)|}_{z=z_{\text{air}}}=0$ (4)

$n\cdot {\left(c\nabla E_y\right)|}_{z=z_{\text{air}}}=n\cdot {\left(c\nabla E_y\right)|}_{z=z_{\text{f}}}$ (5)

在地下介质分界面，位移–应力–电磁分量连续，

$n\cdot {\left(c\nabla B_y\right)|}_{z=z_{\text{f(n)}}}=n\cdot {\left(c\nabla E_y\right)|}_{z=z_{\text{f(n}+1)}}$ (6)

式(4)~(6)中，n是界面法线方向单位矢量，z表示两种不同介质共同界面的位置，下标air表示空气，f表示地下与空气接触的介质，f(n)表示地下第n类介质，f(n + 1)同理。

以v为权函数对计算域和边界积分，通过分部积分得到式(3)和边界条件的弱解形式

${\displaystyle \underset{\Omega }{\iint }\left(\frac{\partial v}{\partial x}d\frac{\partial B_y}{\partial x}+\frac{\partial v}{\partial z}d\frac{\partial B_y}{\partial z}\right)}+{\displaystyle \underset{\Omega }{\iint }veB\text{d}\Omega }-{\displaystyle \underset{\Omega }{\iint }vf\text{d}\Omega }=0$ (7)

2.2. 有限单元分析

求解域采用三角形单元剖分，剖分单元如图1所示。

Figure 1. Schematic diagram of triangular element

图1. 三角形单元示意图

节点坐标按逆时针依次是： $\left(x_i,y_i\right)$ 、 $\left(x_j,y_j\right)$ 和 $\left(x_k,y_k\right)$。面积用 $\Delta$ 表示。 $N_1^c\left(x,y\right),N_2^c\left(x,y\right),N_3^c\left(x,y\right)$ 分别定义分别为小三角形 ${\Delta }_{pjk},{\Delta }_{pki},{\Delta }_{pij}$ 与 $\Delta$ 的比值，易得

$\{\begin{array}{l}{N}_1^c\left(x,y\right)=\frac{1}{2\Delta }\left(a_{i}x+b_{i}y+c_i\right)\\ N_2^c\left(x,y\right)=\frac{1}{2\Delta }\left(a_{j}x+b_{j}y+c_j\right)\\ N_3^c\left(x,y\right)=\frac{1}{2\Delta }\left(a_{k}x+b_{k}y+c_k\right)\end{array}$ (8)

式(4)中，

$a_i=y_j-y_k,b_i=x_k-x_j,c_i=x_j{y}_k-x_k{y}_j;$

$a_j=y_k-y_i,b_j=x_i-x_k,c_j=x_k{y}_i-x_i{y}_k;$

$a_k=y_i-y_j,b_k=x_j-x_i,c_k=x_i{y}_j-x_j{y}_i;$

$\Delta =\frac{1}{2}\left(a_i{b}_j-a_j{b}_i\right)$

用Galerkin法的插值函数对式(7)进行展开，权函数 $v=N_i^e,i=1,2,3$。则三角单元c有积分如下

${\displaystyle \underset{j=1}{\overset{3}{\sum }}{B}_j}{\displaystyle \underset{{\Omega }_c}{\iint }d\left(\frac{\partial N_i^c}{\partial x}\frac{\partial N_j^c}{\partial x}+\frac{\partial N_i^c}{\partial z}\frac{\partial N_j^c}{\partial z}\right)\text{d}\Omega }+{\displaystyle \underset{j=1}{\overset{3}{\sum }}{B}_j}{\displaystyle \underset{\Omega }{\iint }e{N}_i^c{N}_j^c\text{d}\Omega }-{\displaystyle \underset{\Omega }{\iint }f{N}_i^c{N}_j^c\text{d}\Omega }=0$ (9)

假设单元内介质参数与积分无关，则式(9)中积分可记为

$\{\begin{array}{l}{D}_c=d_c{\displaystyle \underset{{\Omega }_c}{\iint }\left(\frac{\partial N_i^c}{\partial x}\frac{\partial N_j^c}{\partial x}+\frac{\partial N_i^c}{\partial z}\frac{\partial N_j^c}{\partial z}\right)\text{d}\Omega }\\ E_c=e_c{\displaystyle \underset{{\Omega }_c}{\iint }{N}_i^c{N}_j^c\text{d}\Omega }\\ F_c=f_c{\displaystyle \underset{{\Omega }_c}{\iint }{N}_i^c{N}_j^c\text{d}\Omega }\end{array}$ (10)

将式(8)代入式(10)中，单元矩阵分别为

$\{\begin{array}{l}{D}_c=\frac{d_c}{4\Delta }\left[\begin{array}{ccc}{a}_i^2+b_i^2& a_i{a}_j+b_i{b}_j& a_i{a}_k+b_i{b}_k\\ a_j{a}_i+b_j{b}_i& a_j^2+b_j^2& a_j{a}_k+b_j{b}_k\\ a_k{a}_i+b_k{b}_i& a_k{a}_j+b_k{b}_j& a_k^2+b_k^2\end{array}\right]\\ E_c=\frac{e_c\Delta }{12}\left[\begin{array}{ccc}2& 1& 1\\ 1& 2& 1\\ 1& 1& 2\end{array}\right]\\ F_c=\frac{f_c\Delta }{12}\left[\begin{array}{ccc}2& 1& 1\\ 1& 2& 1\\ 1& 1& 2\end{array}\right]\end{array}$ (11)

式(3)可以表示成单元矩阵形式

$D_c{B}_c+E_c{B}_c-F_c=0$ (12)

其中，D_c是单元扩散系数矩阵，E_c是单元吸收系数矩阵，F_c是单元源向量，B_c是单元节点场值向量。

2.3. 计算效率

我们在计算域外加上20 m厚的扩散域，计算发现，当扩散域与计算域网格单元尺寸相同时，在一定的时间内不会对计算域动电信号产生干扰。计算域为80 m × 80 m。计算域与吸收层网格大小一样，计算域采用自由三角形网格，吸收层采用规则四边形网格，最大单元尺寸为0.5 m。本文采用并行稀疏直接求解器(MUMPS)求解系数型偏微分方程。程序在两台机器上运行，台式机处理器为6核的Intel (R) Core (TM) i5-10500 CPU @ 3.10 GHz，机带RAM8.00 GB；笔记本电脑处理器是8核的Intel (R) Core (TM) i7-10870 H CPU @ 2.20 GHz，机带RAM16.00 GB。计算效率如表1所示：

Table 1. The computational efficiency of finite element method

表1. 有限元法的计算效率

3. 数值模拟

3.1. 正确性验证

我们在全空间模型下通过互易性检验和与半解析解对比验证计算方法的正确性。激励源为1 N的点单力或1 A∙m⁻²的点电流，坐标点位于r₁ (x = 30 m, y = 50 m)，接收点位于r₂ (x = 70 m, y = 50 m)，源的时间分布由Ricker子波给出。中心频率1000 Hz，时间延迟4 ms，取最大时间t_max为40 ms。解析最小波长λ_min的最大网格单元大小不得超过λ_min/1.5。这样一方面占了源能量的最大部分，另一方面不会产生降低计算速度。我们在域级别上使用二维高斯函数G来近似点源处的狄拉克δ。

$G\left(x,z\right)=\frac{1}{2\pi }{\text{e}}^{\left[-\frac{1}{2}\frac{{\left(x-x_0\right)}^2+{\left(z-z_0\right)}^2}{dS}\right]}$ (13)

互易性检验是验证方法正确性的有效手段，根据

$G_v^C\left(x_r,t;x_s,t_0\right)=G_E^F\left(x_s,t;x_r,t_0\right)$ (14)

式(12)是对Pride & Haartsen [18] 关于互易性描述的变形， $G_v^C$ 代表以电流源为激励的速度场，是位移场对时间的一阶导数， $G_E^F$ 代表以单力点源为激励的电场，x的下标s代表激励源，r代表接收点，t下标0代表源开始激发的时间，无下标表示接收时间。互易性原理认为地震源和电流源及对应的接收点互换位置，前者产生的电场和后者生成的速度场数值相等。

如图2所示，电流源下的电震波场与单力点源下的震电波场，在互换激励源与接收点位置后，得到的电场和速度场波形是一致的，方法的正确性得到验证。

Figure 2. Reciprocity verification of full-space models

图2. 全空间模型的互易性验证

同样，我们与高玉涛等 [14] 做了对比(图3)，说明本文有限元方法的正确性。

Figure 3. Comparison between numerical and semi-analytical solution in a homogeneous medium

图3. 均匀介质全空间模型下数值解和半解析解对比

3.2. 隧道探水数值模拟

隧道掘进过程中遇到的含导水构造主要有断层破碎带、岩溶含水地质构造(溶洞、岩溶管道和地下河)和较复杂的洞穴系统。我们将这三种类型盖化为平板、球体及简单形体的组合体。本文中只涉及二维形式的动电效应数值模拟，所以我们对模型做了二维的简化，如图4所示。数值模拟介质参数见表2，介质2表示含水构造，介质1表示含水构造外其他介质。

Figure 4. Tunnel advance detection-schematic diagram of 2D model of water-bearing structure. (a) Aquifer model; (b) Karst model; (c) Fault fracture zone model

图4. 隧道超前探测–含水构造二维模型示意图。(a) 含水层模型；(b) 岩溶模型；(c) 断层破碎带模型

Table 2. Properties of tunnel water exploration model

表2. 隧道探水模型介质参数

3.2.1. 含水层模型

含水层是隧道掘进过程中可能遇到的地质构造，一般充填泥沙和块石，一旦揭穿，会给隧道和施工人员造成危害。含水层构造一般规模较大，结构简单，所以我们将其作为第一个模拟对象。

图5中(a)、(b)和(c)给出了震电效应中三个不同时刻的波场快照，在17.5 ms时，地震波恰好传播到含水层。可以看到，在含水层–围岩界面处，有蓝色亮斑，此处接收到的电信号E_y要强于计算域其他地方。主要原因是含水层激发的动电效应更强，同时电信号是同震电场与非同震电场的叠加。在蓝色亮斑外还有一圈红色的“晕”，即界面响应的表现，是由S波激发出磁场，再经过磁场感应出电场，在地震剖面图上显示为“平”的事件(图6(a))；在20 ms时，随着S波继续向含水层传播，电场相对于17.5 ms时更强(图5(b))，但是在含水层与围岩下界面，颜色相对于周围更“暗”，主要原因是下层围岩激发的动电效应相对于含水层而言更弱。在22.5 ms时，含水层和围岩中同时激发的动电效应强度对比更加鲜明(图5(c))。

图5中(d)、(e)和(f)显示了三个不同时刻电震效应的波场快照。(d)图是4.0625 ms时(Ricker子波延时4 ms)波场的分布，与震电效应有很大的不同。由于一开始激励源处激发电磁场，“瞬时”(相对于地震波波速大三个数量级)传到了含水层-围岩界面处，并于上下界面形成平行于界面的地震波向反方向传播，在20 ms到达隧道掌子面。在20 ms附近(图5(e))，不仅有地震波反射回隧道前方，同时还有在隧道掌子面边缘处拐角形成的球面波。这种现象在(f)图显示更加明显。

图6显示了x轴坐标位于45 m至55 m，间隔0.5 m，z轴坐标位于60 m接收排列动电剖面图。(a)记录了震电效应的E_y，(b)是记录了电震效应的u_y。电场E_y中同时存在辐射电场信号、界面响应信号和同震电场信号。前两者同向轴斜率接近0，后者有一定斜率。电震剖面图中记载了多处“平”的事件，这是因为地震波前面平行含水层界面传播。在4 ms附近为辐射电磁波信号，17.5 ms时上界面处激发的地震波到达接收排列，20 ms时下界面激发的地震波传播至接收排列。

图7展示了坐标为(45 m, 65 m)的接收点的动电信号。(a)是震电效应的电场记录，(b)是电震效应的位移场记录。(a)中可以看出，18 ms左右对应含水层上界面处界面响应，20 ms是含水层下界面处界面响应；(b)中，20 ms左右对应含水层上界面处电震效应转换处的S波到达接收点，22 ms对应含水层下界面处记录的S波。这样的对比也体现出动电效应的优势。虽然震电和电震在一定条件下可以相互转换，但是在面对不同地质结构时，由于传播规律不一样，可以提供不同但是却紧密相关的信息。之所以震电效应中界面效应会快2 s记录到，是因为这个时间对应的是SH波从震源位置传播到相对于含水层垂点，而电震记录中这个时间对应的是SH波从含水层上界面传播到接收点。

Figure 5. The wave field snapshots of seismoelectric effect and electroseismic effect at t = 17.5 ms, t = 20 ms and t = 22.5 ms, where (a), (b) and (c) are electric fields under seismoelectric effect, and (e), (f) and (g) are displacement fields under electroseismic effect

图5. 震电效应和电震效应在t = 17.5 ms、t = 20 ms和t = 22.5 ms时刻下的波场快照，其中(a)、(b)和(c)是震电效应下的电场，(e)、(f)和(g)是电震效应下的位移场

Figure 6. Electrokinetic effect record of receiver array (parallel to the tunnel plane, z = 60 m, x = (45:5:55) m). (a) Electric field record under seismoelectric effect; (b) Displacement field record under electroseismic effect

图6. 接收排列(平行于隧道掌子面，z = 60 m，x = (45:5:55) m)的动电效应记录。(a) 震电效应下的电场记录；(b) 电震效应下的位移场纪录

Figure 7. electrokinetic effect record of receiving point at (45 m, 65 m). (a) Electric field record under seismoelectric effect; (b) Displacement field record under electroseismic effect

图7. 接收点位于(45 m，65 m)的动电效应记录。(a) 震电效应下的电场记录；(b) 电震效应下的位移场纪录

3.2.2. 岩溶模型

岩溶形体一般有溶洞、地下河、岩溶管道三类。溶洞一般沿地质结构面发育，有层面溶洞和沿构造裂隙溶洞两类，断面一般为近圆形或椭圆形。地下河是接近水平的洞穴系统，二维情况下一般可用平板模型代替，其他情况我们使用圆代替(二维模型)。岩溶通道是由洞穴和管道联系起来的洞穴系统，直径0.01~0.1 m腔体为溶管，直径0.1~1 m左右腔体称为管道，二维情况下我们当作圆处理。隧道施工过程中，如果揭露岩溶管道或导致地下水流方向改变，会造成大型突水突泥事故。

图8中子图(a)、(b)和(c)显示了岩溶模型不同时刻震电效应的波场快照。(a)图中，当地震波到达岩溶构造区域时，在地震波前面与岩溶构造相切的地方，有界面响应产生，激发的电磁波“瞬时”传到了接收点。(b)图中地震波传到了岩溶构造远离激励源的另一端，且波前面与岩溶构造界面相切，同样产生界面响应。在(a)图所代表的12.5 ms和(b)图所代表的16.25 ms时间内，地震波在岩溶构造-围岩界面上引起同震电场，同时反射回围岩，故波场快照中岩溶构造外部有环状的地震波场，在(b)图及20 ms的(c)图有明显表现。

Figure 8. The wave field snapshots of seismoelectric effect and electroseismic effect, where (a), (b) and (c) are electric fields under seismoelectric effect at t = 12.5 ms, t = 16.25 ms and t = 20 ms respectively, and (e), (f) and (g) are displacement fields under electroseismic effect at t = 4.0625 ms, t = 10 ms and t = 20 ms respectively

图8. 震电效应和电震效应的波场快照，其中(a)、(b)和(c)是震电效应在t = 12.5 ms、t = 16.25 ms和t = 20 ms时刻下的电场，(e)、(f)和(g)是电震效应在t = 4.0625 ms、t = 10 ms和t = 20 ms时刻下的位移场

图8中后三幅子图分别显示了4.0625 ms、10 ms和20 ms时电震效应引起的位移场的波场快照。与震电效应形成的电场传播规律不同的是，电震效应中电磁波瞬时传到了岩溶构造，然后以岩溶构造为“源”，向外传播地震波。我们在(d)图中可以看到两处波场振幅较大处，分别是激励源的位置和岩溶构造的位置。(e)图中显示了两道由岩溶构造区域向外传播的地震波，且波峰之间距离恰好为岩溶构造的直径，原因是外围的地震波是激励源处电磁波“瞬时”到达岩溶构造后转换为地震波向外传播，类似于震电效应中的界面响应，而后至的地震波是界面处激发的地震波向岩溶构造内传播，然后一部分能量继续透射处岩溶构造，向最外围的波场反方向传播，由于圆的对称性，显示为仅相差一个岩溶构造直径距离的多道地震波。(f)图对不同道地震波场的强弱显示更为明显，同时也显示了隧道拐角处可以激发球面波。

图9显示了x轴坐标位于45 m至55 m，间隔0.5 m，z轴坐标位于60 m接收排列的动电记录。(a)是震电效应下的电场E_y，(b)是电震效应下的位移场u_y。(a)图中我们可以看到四个“平”的事件，分别是位于4 ms、11 ms、15 ms和17.5 ms附近。第一个为辐射电场；第二个对应地震波波前面与岩溶构造相切时激发的界面响应，由此产生的电场“瞬时”传到接收排列；第三个为地震波前离开岩溶构造、与界面相切时，产生第二次界面响应，电场“同时”被接收排列接收；第四个是岩溶构造内部的地震波在遇到界面后反射回岩溶构造内部，再次回到第一次激发界面响应的地方，激发出第三次界面响应。(b)图显示了一次“平”的事件和3次不水平的同向轴。“平”的时间对应于电磁波“瞬时”到达接收排列，并激发出地震波；3次非水平的同向轴解释同图8后三幅子图显示的电震效应的波场快照。从两幅动电记录图中波至时间和相对振幅的对比可以论证上述结论，同时也说明了震电和电震的互易性。

Figure 9. Electrokinetic effect record of receiver array (parallel to the tunnel plane, z = 60 m, x = (45:5:55) m). (a) Electric field record under seismoelectric effect; (b) Displacement field record under electroseismic effect

图9. 接收排列(平行于隧道掌子面，z = 60 m，x = (45:5:55) m)的动电效应记录。(a) 震电效应下的电场记录；(b) 电震效应下的位移场纪录

图10显示了位于(45 m, 65 m)处接收点的动电信号，(a)为震电效应下的电场E_y，(b)为电震效应下的位移场u_y。波到时与成因与上述其他两种表现形式类似。从这组图中我们主要是能够更清晰地看出波幅值的差异。以震电效应为例，相对于PSVTM模式，SHTE模式辐射电磁场大于界面响应和同震电磁场信号，且界面响应信号大于同震电磁信号。电震效应同理。

Figure 10. Electrokinetic effect record of receiving point at (45 m, 65 m). (a) Electric field record under seismoelectric effect; (b) Displacement field record under electroseismic effect

图10. 接收点位于(45 m, 65 m)的动电效应记录。(a) 震电效应下的电场记录；(b) 电震效应下的位移场纪录

3.2.3. 断层破碎带模型

断层是地壳受力发生断裂，沿断裂面两侧岩块发生的显著相对位移的构造，从小于1 m至数百千米规模不等。李术才等 [17] 将断层类致灾构造分为富水断层型、导水断层型和阻水断层型。其中，阻水型断层会使得含水层排泄通道不畅，地下水位富集，水位上升。一旦隧道揭穿阻水结构，断层阻水作用失效，隧道即成为地下水排泄空间，发生突水突泥。研究断层破碎带动电效应对快速准确避免突水突泥灾害具有重大意义。我们在设计二维断层破碎带模型时，将破碎带和断层均设计为含水构造，此外，隧道按当前方向掘进可揭穿断层结构。

图11中前三幅子图是断层模型震电效应不同时刻的波场快照，在7.0313 ms时，地震波到达破碎带边缘，引起界面响应，在破碎带内部可以看到红色亮斑，是同震电场和界面响应共同贡献的效果。此外，计算域空间显示为蓝色，是因为界面处电磁波以电磁波速度传遍计算域，故此时在任一点放置接收器均可接收到电场信号。在8.4375 ms时，我们看到能量主要聚集在断层破碎带中，主要是因为断层破碎带能激发更强的动电响应，此外，上断层与围岩界面处界面响应明显。在10 ms时，地震波到达下断层下界面，激发出界面响应。

图11(d)~(f)显示了断层模型电震效应不同时刻的结果。(d)图显示激励源发出电流，“瞬时”就传播到断层破碎带，并引起电震效应。(e)图显示了6.25 ms时界面处的电震效应，含水层-围岩处电磁波转换而来的地震波向相反方向运动，波前面平行于界面。(f)图显示了10 ms时位移场，与(e)图相比，位移场波场分布显示得更加清晰。

图12显示了x轴坐标位于45 m至55 m，间隔0.5 m，z轴坐标位于60 m接收排列的动电记录。(a)是震电效应下的电场E_y，(b)是电震效应下的位移场u_y。(a)图中我们可以看到三处较为明显的斜率为0的同向轴，分别是辐射电场、破碎带上界面激发的界面响应和上断层上界面与下断层上界面的界面响应的叠加。(b)图与(a)图存在较多的不同，除了与电流源直接相关的S波有关外，位移场记录了两道斜率不为零且不相等的同向轴。在5 ms至9 ms时间段的同向轴为破碎带激发出的S波，另外一条斜率较缓的位于10 ms附近的同向轴是下断层的下界面处激发的S波。

Figure 11. The wave field snapshots of seismoelectric effect and electroseismic effect, where (a), (b) and (c) are electric fields under seismoelectric effect at t = 7.0313 ms, t = 10 ms and t = 8.4375 ms respectively, and (e), (f) and (g) are displacement fields under electroseismic effect at t = 4.0625 ms, t = 6.25 ms and t = 10 ms respectively

图11. 震电效应和电震效应的波场快照，其中(a)、(b)和(c)是震电效应在t = 7.0313 ms、t = 10 ms和t = 8.4375 ms时刻下的电场，(e)、(f)和(g)是电震效应在t = 4.0625 ms、t = 6.25 ms和t = 10 ms时刻下的位移场

Figure 12. Elctrokinetic effect record of receiver array (parallel to the tunnel plane, z = 60 m, x = (45:5:55) m). (a) Electric field record under seismoelectric effect; (b) Displacement field record under electroseismic effect

图12. 接收排列(平行于隧道掌子面，z = 60 m，x = (45:5:55) m)的动电效应记录。(a) 震电效应下的电场记录；(b) 电震效应下的位移场纪录

图13展示了位于(45 m, 65 m)接收点的动电信号，(a)为震电效应下的电场E_y，(b)为电震效应下的位移场u_y。(a)图中，在4 ms处记录了与震源直接相关的辐射电磁波，7 ms左右记录到破碎带处的界面响应记录，由于破碎带的厚度没有达到震电效应电磁波的分辨率，故以电场振幅加强形式存在。10 ms附近记录的是下断层下界面处界面响应。对于电震效应来说，相对于震电效应，分析复杂结构要更加困难，因为电流源激发的电磁场在含水构造中进行动电转换后以界面轮廓的形式向外部传播。(b)图中记录了三种较为明确的波形，4 ms与电流源相关的位移场，和震电效应中的辐射电场对应，7 ms左右由破碎带上界面激发出的S波，11 ms左右由下断层下界面激发的S波。从波场快照中可以对单点接收的动电记录进行验证。

Figure 13. Electrokinetic effect record of receiving point at (45 m, 65 m). (a) Electric field record under seismoelectric effect; (b) Displacement field record under electroseismic effect

图13. 接收点位于(45 m, 65 m)的动电效应记录。(a) 震电效应下的电场记录；(b) 电震效应下的位移场纪录

3.3. 含水构造定位

在隧道施工过程中，伴随着各种震动，比如爆破、冲击钻等，我们可以利用这些作为震源进行震电勘探，与施工过程中产生的其他震动相比，具有震源位置明确，不占用施工场地、不影响工程进度和强度较大的优点。或者我们采用锤击震源与电流源，在人工选取合适震源位置后，可以实现多点激发和振幅叠加。隧道动电探测激励源和接收排列布置如图14所示。在施工过程中，可按如下步骤判断含水构造大致位置及粗略的几何特征：

Figure 14. Location diagram of water-bearing structure

图14. 含水构造定位示意图

1) 设置两组接收排列，一组接收排列平行于隧道掌子面，另一组垂直于掌子面，分别记为Array₁和Array₂；

2) 以S₁作为激励源位置，通过平行掌子面接收排列的地震记录判断含水构造与震源的相对位置关系;

3) 将垂直于隧道掌子面的接收排列布置靠近在含水构造一方；

4) 在S₁震源位置施加电性源，分析两组接收排列的位移场信号，若Array₁或Array₂剖面图同向轴水平，可以说明含水构造为平行或垂直于隧道掌子面；若同相轴有斜率，且斜率相同，可以判断含水构造有规则平面，但是平面与排列不成特殊角度关系，通过斜率可以计算出交角；若斜率不同，可以判定含水构造为不规则体；

5) 对电震和震电进行多次操作，将剖面记录叠加。如果是震电，更换源的位置，拾取接收排列界面响应信息重构含水构造轮廓；如果是电震，在S₁位置激发，通过同样在S₁激发的震电界面响应信息确定对应同向轴意义，然后对各接收点位移场做时深转换，进而描述含水构造轮廓。具体操作及分析以断层破碎带定位为例。

对于含水构造定位，我们提供了震电和电震两种方式。在隧道中，我们将接收点记作R₁~R₁₂号，分别对应(45 m, 90 m)，(45 m, 85 m)，(45 m, 80 m)，(45 m, 75 m)，(45 m, 70 m)，(45 m, 65 m)，(45 m, 60 m)，(47 m, 60 m)，(49 m, 60 m)，(51 m, 60 m)，(53 m, 60 m)，(55 m, 60 m)，震源位于(50 m, 60 m)。以电震为例，我们可以提取接收点界面响应时间，进行时深转换，以接收点所在位置为圆心，界面响应转换后得到的深度为半径画圆弧，标记为A₁~A₁₁。利用A₁~A₁₁可以基本确定靠近震源一方的破碎带的轮廓L₁。在电震效应中，如果利用上述的方法，我们并不能得到上断层靠近震源的一侧的轮廓。因为当电流源激发时，在含水层引起电震效应，转换后的地震波波前面会按照L₁与L₂的轮廓向电流源方向传播，所以位于隧道掌子面处的接收点R₇~R₁₂不会接收到来自L₂的信号。这时，我们需要震电信号作为补充，假设在(50 m, 60 m)处沿Z轴向下锤击，即F_z = 1 N的单力点源，因为自震源处向L₂可以引垂线，所以在上断层与围岩的上界面处会产生界面响应，假设震源位置分别在R₁₀~R₁₂处，以震源为圆心，向隧道掌子面前方画弧，即可定位L₂。

含水构造定位要充分利用震电和电震的信息。震电效应中的位移场信息与电场信息间存在一定的数量关系，比如，同震电场的到时与位移场是一致的；如果接收排列与发生界面响应平面平行，那么界面响应到时仅为同震电场到时一半。同时，根据互易性原理，震电效应中电场与电震效应中的速度场到时相同，可以互相验证，确定波类型，明确对应到时与波速，从而更准确地提取距离信息。

4. 结论

本文采用有限单元法模拟了单力点源下激发的震电和电震波场。我们首先在全空间模型下验证了方法的可行性，然后针对隧道掘进过程中可能遇到的三种突水情形做数值模拟，通过正演结果对异常体定位，并针对震电效应和电震效应分别做了阐述，结论如下：

1) 在用于隧道超前探水的震电勘探中，我们可以利用同震电磁信号和界面响应信号。波场快照上，可以将前者看做以地震波的速度传播，后者以电磁波速度传播。在二维SHTE模式中，有辐射电磁波、水平偏振的横波和由水平偏振横波转换的电磁波。在不了解含水构造情况下，我们要充分利用震电效应和电震效应这两种勘探方法，通过动电效应为弹性波与电磁波提供的天然耦合关系寻找波场之间的数量关系，从而确定含水构造位置和轮廓。

2) 对于震电信号而言，利用界面响应信号，多次更换震源位置，利用电场界面响应信号中震源与含水构造–围岩界面特殊的几何关系，可以实现含水构造定位。

3) 对于电震信号，一次激发，多点接收，利用位移场信号中接收点与含水构造–围围岩界面特殊的几何关系，同样可以实现定位。

参考文献