一种基于时间域反射系数的穿煤隧道瞬变电磁超前含水夹层探测方法

本发明是一种基于时间域反射系数的穿煤隧道瞬变电磁超前含水夹层探测方法，通过此方法能反应岩层分层界面的电性分布情况，解释探测前方的水害情况，属于隧道超前地质预报。

背景技术：

0、技术背景

1、我国煤炭资源分布较广，目前国内高速公路、高速铁路的隧道建设常常会穿越煤层。煤层是一类特殊的岩土介质，具有岩石质地脆弱，结构破碎，遇水软化，工程力学性能差等特点。隧道开挖活动中，掌子面前方多存在充水、断层或裂隙等地质构造形成的含水夹层，若煤层存在渗透性突变区域，则水可以通过这些区域渗透到隧道中，支护不当的情况下，引发突水灾害。运用瞬变电磁法对穿煤隧道掌子面进行正演模拟，能够通过电磁信号的传播和响应，获取煤层电性分布信息，基于含水构造的导电性等物理性质与煤层存在差异，可判断含水体的位置。

2、在地球物理勘探中，岩石中的水或其他导电的矿物质有一定的电导率，这导致在电磁波传播过程中发生能量损耗，表现为电磁波的衰减。为了使模拟能准确，计算在有耗介质中的时域反射系数。同时，通过时间域反射系数，能更好地分析边值问题，比如岩层分界面对电磁波的影响。运用有耗介质中的直接时间域反射系数正演方法进行瞬变电磁正演模拟使模型更符合实际情况，有助于深入理解电磁场在地下岩层中的传播过程，得到岩层分层界面的电性分布情况，为掌子面前方地质结构进行解释预报提供参考价值，保障隧道施工安全。

技术实现思路

1、本发明提供了一种基于时间域反射系数的穿煤隧道瞬变电磁超前含水夹层探测方法，是针对探测穿煤隧道掌子面前方含水夹层提出的一种探测方法，以解决正演模拟计算精度及速度问题。

2、本发明的技术方案是：一种基于时间域反射系数的穿煤隧道瞬变电磁超前含水夹层探测方法，此方法的具体步骤如下：

3、step1、建立穿煤隧道掌子面煤层-含水夹层-煤层tem模型，设置激发源，描述地下结构。

4、建立，贴近隧道掌子面处布置一直立圆形线圈，其中通以阶跃电流。

5、

6、掌子面上的发射线圈在t<0时刻持续通电，在等于零的瞬间停止通电。由于电磁感应，可在掌子面前方三维空间中得到时间域阶跃电流激发下的瞬变电磁响应。

7、掌子面前方岩层为三层，中间层为含水夹层(即为异构体)，第一层和第三层为煤层。

8、step2、采用加权平均法，计算煤层-含水夹层-煤层结构的等效电磁参数，将step1的模型等效为均匀模型，然后通过直接时间域矢量有限元法求解均匀模型的瞬变电磁响应；

9、分层岩层瞬变电磁正演模拟中，需要根据几何模型，对整个模拟区域进行网格划分，分层界面的网格划分是关键的一步，它直接影响到模拟结果的准确性。同时，分层界面通常是电磁响应变化显著的地方，需要在这些区域进行精细的网格划分，但是更高分辨率的网格可能需要更多的计算资源。因此，通过加权平均法计算等效电磁参数，这种方法考虑了每个层的贡献，并根据各层的重要性进行了加权，更准确地反映分层介质的整体行为。将煤层-含水夹层-煤层等效为均匀模型可以避免分层界面进行网格划分，减少计算的复杂性，提高计算速度，尤其是在大范围的勘探和模拟中。

10、step2.1、模型中的煤层-含水夹层-煤层，每一层具有不同的电磁参数，电导率、磁导率和介电常数，采用加权平均法，计算整个岩层结构的等效电磁参数。

11、将岩层视为平行于掌子面的平面，计算每一层的体积，即厚度hi与平面面积s的乘积；计算总体积，将所有层的体积相加，得到整个分层岩层的总体积v；计算每一层的体积分数fi，其体积分数等于该层的体积除以总体积。公式如下：

12、

13、进一步计算各电磁参数的加权平均值，得到等效电磁参数：

14、

15、其中，σeff为等效电导率；μeff为等效磁导率；εeff为等效介电常数。

16、step2.2、从时间域麦克斯韦方程组出发，推导时域电磁场控制方程，引入第一类狄利克雷边界条件，并应用伽辽金法，得到矢量有限元方程。

17、时域麦克斯韦方程为：

18、

19、麦克斯韦方程组中的物质关系方程为：

20、

21、其中，b为磁感应强度(wb/m2)；e为电场强度(v/m)；h为磁场强度(a/m)；d为电位移矢量(c/m2)；j为电流密度；σ为电导率(s/m)；μ为磁导率(h/m)；ε为介电常数(h/m)。

22、将等效电磁参数代入(5)式得到均匀模型的物质方程：

23、

24、将其运用于麦克斯韦方程中进行求解。

25、在各向同性、均匀的非色散介质中，不考虑位移电流的条件下，将中心回线源作为外加电流源，时间域电磁场满足如下麦克斯韦方程：

26、

27、其中，μ0为自由空间磁导率(h/m)；js为源电流密度。

28、为了描述电磁场在无源区介质交界面上的行为，确保电场和磁场的光滑性和连续性，满足以下边界条件，解出电磁场：

29、

30、其中，界面单位法向量为n。

31、在无穷远处，即电场或者磁场在无穷远边界上的切向分量为零，采用狄利克雷边界条件：

32、

33、将(7)式电磁方程改写为：

34、

35、

36、对(10)式两边取旋度，然后将(11)式代入可得：

37、

38、运用矢量恒等式得到电场控制方程：

39、

40、时域电磁场根据扩散规律在地层中传播方法在时域内直接求解电场扩散方程。

41、对于step2.2中的矢量有限元方程，利用galerkin加权余量法，建立电场矢量波方程，并将其与矢量恒等式和散度定律相结合，建立电场galerkin方程。

42、时域电场控制方程相应的余量为：

43、

44、通过对计算区域中的单个四面体单元的余量进行加权积分，得到有限元矢量变分方程：

45、

46、其中，f表示矢量基函数。

47、令b＝f，根据矢量恒等式将(15)式的第一项积分进行分布积分：

48、

49、然后根据高斯公式，将(16)式中的第二部分积分转化为面积分：

50、

51、通过矢量公式a·(b×c)＝(a×b)·c，引入无穷远的边界条件e(r,t)|γ＝0，可将(17)式化简为：

52、

53、最后，得出有限元法的矢量变分控制方程：

54、

55、其中，μ0为自由空间磁导率(h/m)；js为源电流密度；f为矢量基函数；t为响应时间；

56、step2.3、将计算区域剖分成有限个单元，采用非结构化的一阶四面体单元进行剖分计算区域，采用矢量whitney型插值基函数，在每个小单元内逼近解的分布。

57、对于任意给定的一个四面体单元，四个顶点的编号分别为i＝1,2,3,4，而在相应的顶点上对应的棱边编号分别为(1)、(2)、(3)、(4)、(5)、(6)。

58、剖分单元内的函数为f(x,y,z),将电场切线方向上的自由度加到各单元的棱边，求解区域内任一四面体单元ωe内的函数可近似表示为：

59、fe(x,y,z)＝ae+bex+cey+dez      (20)

60、联立以下四个顶点函数，求解出待定系数ae、be、ce、de：

61、

62、其中，fie为第i个顶点上的函数值。

63、将求解的待定系数代入(21)式得：

64、

65、其中，四面体内的函数ωe单元中第i个节点的插值基函数为：

66、

67、

68、最终求出该插值基函数表达为：

69、

70、当观察点位于节点i，基函数当观察点位于与节点i相对的三角形内，基函数

71、令体积坐标可以得到：

72、

73、四面体单元内的任一棱边i，whitney矢量基函数为：

74、

75、其中，第i棱边的两个端点表示为i1和i2；表示连接点的棱边长。在第i棱上，其矢量基函数在其余5个棱边上，

76、有限元数值解被表示为空间和时间函数的展开，可以用上述矢量基函数将四面体内任一点的电场强度进行展开：

77、

78、其中，为位置r处的基函数，用于描述电场在空间中的分布；为不同时间的电场是四面体单元中棱边i上的投影。

79、step2.4、采用后向差分格式将偏微分方程转化为差分形式，得到矩阵方程，计算区域单元分解后，将单元矩阵组装成整个系统的全局矩阵，即生成整体刚度矩阵。

80、在时域电磁场矢量有限元控制方程中，存在电场随时间变化的偏导，必须先利用后向差分方法对其进行时间离散，将代入(19)式得到有限元矩阵方程：

81、

82、转化为矩阵形式：

83、

84、

85、

86、其中，为四面体单元每条棱边上第n时刻的电场投射值的列向量；jn为

87、第n个时刻的电流密度；ae为系数矩阵。

88、将ae分解为a1e和a2e：

89、

90、

91、a1e转换为矩阵形式,分解化简后得到:

92、

93、其中，k1、k2、k3分别为：

94、

95、a2e转换为矩阵形式,分解化简后得到:

96、

97、其中，k4为：

98、

99、同理，可以将矩阵be分解为b1e和b2e：

100、

101、

102、b1e转换为矩阵形式,分解化简后得到:

103、

104、其中，为已知量是上一时刻求得的电场值，表示为：

105、

106、b2e转换为矩阵形式,分解化简后得到:

107、

108、将每个单元矩阵汇总整合，形成整个计算区域的刚度矩阵，整体刚度矩阵与节点位移关联起来，形成大型线性方程组。

109、ae＝b或

110、其中，ne为剖分四面体单元数量。

111、step2.5、采取直接法mumps求解器来求解由有限元法离散化产生的线性方程组，得到每个时间步的场解。

112、瞬变电磁探测一般采集到的是感应电动势，需要用法拉第定律将每个时间步的电场强度解转换为磁感应强度。各个方向的磁感应强度变化率为：

113、

114、基于以上步骤得到均匀半空间等效模型的瞬变电磁响应，磁感应强度变化率

115、step3、考虑电磁场在介质中传播时的响应，通过时间域传播模式法，来得到电、磁在分层有耗岩层中的时域传播模式矩阵，构建直接时间域反射系数不同极化模式方程，利用数值积分算法求解，得到不同极化模式的煤层-含水夹层-煤层界面时域反射系数。

116、所述step3中包括，推导有耗分层岩层的时域传播模式矩阵来得出直接时间域反射系数不同极化模式方程，用时间域反射系数来描述平面波在有耗煤层-含水夹层-煤层中的时间域电磁响应，具体包括：

117、通过时间域反射系数来描述煤层-含水夹层-煤层的瞬变电磁响应；

118、在tem探测中，一次场是由发射线圈激发的电磁场产生，反射系数通常涉及一次场的产生。因此可通过一次场求取分层有耗岩层时间域反射系数。

119、电磁波垂直入射到分层有耗岩层表面。为表示瞬变电磁中心回线源激发的平面波的传播，用如下麦克斯韦方程描述电磁场在岩层中的行为，其中岩层的电磁性质由电导率σ和介电常数ε决定，考虑了电磁场在岩层中传播时的响应，进而推导电、磁在分层有耗岩层的时域传播模式矩阵。

120、

121、为了描述不同介质中波的传播特性，确保波在传播过程中保持相位一致性。电场和磁场可以表示为如下形式：

122、

123、其中，电磁场是时域t和空间z的函数，且是一个沿z向传播的平面波；c为波速。

124、半空间分界面取z＝0，和z＝0-和z＝0+分别表示掌子面前半空间和后半空间，z＝0-为真空区域。半空间分界面取z＝0，和z＝0-和z＝0+分别表示掌子面前半空间和后半空间，z＝0-为真空区域。e+表示半空间界面的透射波，e-表示反射波，即：

125、

126、其中，和分别代表反射系数中与时间有关的部分和与时间无关的部分；“*”是卷积积分运算符。

127、在tem探测中，反射系数通常由一次场产生，一次场的电场e和磁场h由发射线圈激发产生，根据切向电场和磁场在界面上的连续性，描述界面处电场的关系，表达电场在界面处的变化和整个电磁场之间的关系，得出：

128、

129、其中，t1、t2是对应于两种不同损耗介质的传输矩阵，其参数分别为ε1、μ1、σ1和ε2、μ2、σ2；元素α和β是用阻抗来构建的传输矩阵；由于在半空间中没有上行波(z＞0)，e-(0+,t)＝0。即(49)式转化为:

130、

131、由(50)式即可得到电、磁在三层有耗岩层中的时域传播模式矩阵，可得：

132、

133、(50)式中

134、

135、其中，i为单位矩阵；和为阻抗矩阵，和为和的逆矩阵。

136、用阻抗来构建传输矩阵的元素α和β分别为：

137、

138、

139、其中，是由za和zb构成的总体阻抗矩阵，描述了整个岩层中传播波和反射波的行为；和可由(55)式确定：

140、

141、za和zb分别为：

142、

143、

144、其中，ε、σ为三层岩层总的介电常数和电导率，总的介电常数和电导率来捕捉多层介质对电磁波传播的影响，可以将其视为step2.1中所求的等效介电常数εeff和等效电磁率σeff；za为层内传播的波阻抗矩阵，zb为一层传播到下一层的边界条件；ηm为磁场的波阻抗，ηe电场的波阻抗，均与时间t相关。

145、

146、

147、其中，i0和i1为修正的贝塞尔函数；s为拉普拉斯域中的复变量；t'为描述积分范围的虚拟变量。

148、是的逆矩阵，通过定义来进行求解：

149、

150、(60)式中ξm和ξe满足第二类伏尔泰方程，表达了电场和磁场的传统阻抗与考虑介电常数变化修正后的波阻抗之间的关系。

151、

152、其中，η′e为修正后的电场波阻抗，η′m为修正后的电磁波阻抗，考虑了介电常数的影响。

153、记解出上式方程可得：

154、

155、(62)式两边拉普拉斯变换可得：

156、

157、应用拉普拉斯变换

158、

159、将(58)式和(59)式分别代入(63)式可得：

160、

161、

162、将(64)式和(65)式代入(60)式可以得到将(53)式和(54)式代入(51)式，可得：

163、

164、因为e+(0-,t)具有任意性，将(66)式化简为以下两个方程：

165、

166、

167、求解(67)式可得：

168、

169、于是在各向同性岩层中：

170、

171、(68)式整理后可得：

172、

173、其中

174、

175、且

176、

177、

178、(73)和(74)式中ε1'＝ε1-ε0，ε2'＝ε1-ε0，并且

179、

180、

181、将(75)和(76)式分别代入(73)和(74)式可得：

182、

183、

184、上式中χe(t)和χm(t)逆拉普拉斯变换可得：

185、

186、类似的式(71)式可以分为te模式和tm模式

187、

188、其中，rtm(t)、rte(t)分别为水平极化和垂直极化的时间域反射系数；rtm、rte分别为水平极化和垂直极化与时间无关部分的反射系数；ηtm、ηte分别为为磁场的波阻抗和电场的波阻抗。

189、通过数值积分算法求解(80)式得到不同极化模式下的有耗分层岩层界面时间域反射系数。

190、根据瞬变电磁掌子面探测需求，选取合适的极化模式，计算这一模式下的时域反射系数。如遇到复杂情况下，可选取两种模式进行计算，进而提供更全面的解释。

191、具体的，通过(80)、(70)、(77)、(78)式可得到分层有耗介质的界面时域反射系数。(80)式为第二类沃尔泰拉方程，求解该方程时，可以利用数值积分方法，运用梯形法则求解该方程：

192、

193、

194、step4、根据掌子面探测需求，选取反射系数极化模式，将求得的平面波分层界面时间域反射系数与均匀等效模型的瞬变电磁响应进行卷积，得到掌子面前方煤层-含水夹层-煤层的瞬变电磁响应。本发明可以通过平均相对误差和计算时间来评价该模型及方法是否有益，提高了正演模拟精确度和速度。

195、瞬变电磁中心回线探测通常可以考虑水平极化和垂直极化两种情况，取决于电流源的激励方式。水平极化表示电流源方向与地面水平，垂直极化表示电流源方向与地面垂直。这两种极化方式在瞬变电磁法中的选择涉及到研究对象的性质、勘测的深度、地下介质的电导率分布等因素。通过变化激发线圈的电流方向，可以获得不同方向上的电磁响应，进而提供更多信息以推断地下介质的性质。

196、通过step3求得垂直、水平极化模式下的平面波分层界面时间域反射系数step2求得均匀半空间等效模型的瞬变电磁响应掌子面前方煤层-含水夹层-煤层垂直、水平极化的瞬变电磁响应通过卷积计算得到和

197、

198、

199、由于穿煤隧道掌子面含水夹层探测，是以掌子面为基准面，探测垂直方向上的岩层界面的变化，因此选择垂直极化进行探测，更能提供详细的地质信息。

200、最后，穿煤隧道掌子面前方煤层-含水夹层-煤层瞬变电磁响应为：

201、

202、传统的瞬变电磁正演方法需要在正演模拟后处理阶段将结果频率域转换为时间域响应，而该方法及模型中均匀半空间等效模型的瞬变电磁响应和时间域反射系数均在时间域上进行数值求解，可避免频率-时间转换误差。最后将两者进行卷积得到掌子面前方煤层-含水夹层-煤层的瞬变电磁响应，这一过程通过数值方法进行近似计算。

203、通过计算时间t0来评价正演模拟的计算速度，模拟值与实际观测值之间的平均相对误差mae来评价正演模拟的精确度：

204、

205、其中，n为数据点的数量；xi为模拟值；yi为实际观测值。

206、若正演模拟结果经评价达到规定量化指标，则表明使用穿煤隧道瞬变电磁超前含水夹层探测时间域反射系数法进行正演模拟具有有益效果。

207、本发明的有益效果是：首先隧道掌子面前方煤层-含水夹层-煤层等效模型的提出，可以避免正演模拟准备阶段分层界面网格的划分，进而减少计算阶段的复杂性，提高正演模拟的计算速度；然后，时间域反射系数应用于隧道超前地质预报中,计算平面波在有耗煤层-含水夹层-煤层界面的时间域反射系数，根据地层中的传播和反射进行建模可以提供更详细、更准确的数据，提高正演模拟精确度，同时，根据掌子面探测需求，选取合适的极化模式，计算这一模式下的时域反射系数，通过这种分解可以减少计算的复杂性，若遇到复杂情况，te和tm模式使得在求解电磁场的边值问题时，更容易地考虑和理解两种模式的相互影响，进而提供更全面的解释；最后，将平面波分层界面时间域反射系数与均匀等效模型的瞬变电磁响应进行卷积得到掌子面前方煤层-含水夹层-煤层含水夹层的瞬变电磁响应，因为两者结果均通过数值算法时间域上求解而得，能够避免正演模拟后处理阶段的频率-时间转化误差，得到岩层界面的电性分布情况，为成像解释预报提供瞬变电磁响应正演模拟数据，更加准确地判断含水夹层位置及形态。