隧道地震和瞬变电磁联合探测方法与流程

[0001]
本发明属于隧道超前地质预报和工程勘察技术领域，具体涉及一种隧道地震和瞬变电磁联合探测方法。


背景技术：

[0002]
1.反射地震和瞬变电磁探测的研究意义和研究现状
[0003]
隧道施工时，隧道超前方向的断层破碎带、岩溶溶洞、突水突泥等潜在的危险源极易引发高事故率，造成重大的经济损失和人员安全。反射波地震和瞬变电磁是两种最主要的隧道超前探测方法。隧道地震探测方法因观测系统的不同，主要有水平地震剖面-hsp(inazaki et al.,1999)，真反射层析成像-trt(neil et al.,1999；董远浪.2019)，隧道地震层析成像-tst(zhao et al.,2006；xiao et al.,2012)和tsp方法(dickmann and sander,1996；li et al.,2015；li et al.,2017；esmailzadeh et al.,2018)。然而因隧道空间的限制，现有的隧道地震方法计算出的偏移剖面仅能判定不良地质体的空间位置，无法准确判定裂隙水。
[0004]
瞬变电磁法具有对低阻异常体的分辨率较高，是隧道富水岩溶，破碎带，断层探测最优的探测方法。薛国强等(2008)提出应用“浮动薄板”理论，以二次电导微分参数数据绘制成像剖面。苏茂鑫等(2010)总结和归纳视纵向电导解释方法和特点及其应用效果。孙怀凤等(2011)提出隧道瞬变电磁多点阵列式探测方法，分析含水构造的三维空间展布规律。李术才等(2014)提出隧道瞬变电磁超前预报平行磁场响应探测方法。戚志鹏等(2015)研究隧道瞬变电磁任意共面垂直磁场全空间全域视电阻率解释方法。李凯等(2017)研究隧道瞬变电磁超前探测中不同装置型式对采集数据和成像结果的影响，通过数值试验和现场试验对比了4种掌子面收发装置布设方式：中心回线装置、固源直线接收偶极装置、固源拱形接收偶极装置和固源等距接收偶极装置。张乐文等(2018)设计均匀半空间模型和典型层状模型，测试瞬变电磁法响应探测感应电动势转换成视电阻率的适用性。邢修举等(2019)设计出大断面三维多匝小回线瞬变电磁探测技术，对岩溶隧道区突水问题做出准确的预测。
[0005]
2.地震和瞬变电磁反演的物性基础
[0006]
电磁和地震数据的物性基础之间，存在直接关系和间接关系。
[0007]
2.1.直接关系
[0008]
电磁和地震数据之间至今还没有普遍适用的公式，目前普遍应用的是faust研究总结出的电阻率-速度公式，也就是著名的faust公式：
[0009]
v＝khc
drt
[0010]
推导变换为：
[0011][0012]
faust公式表征在深度约束条件下地层电阻率曲线与声波曲线之间具有统计关系，式中v为声波速度，k、c、d为待定系数，h为深度，rt为波阻抗曲线，这个公式可通过岩石
物性测量和统计相关得到，它只适用于某些特定地区，对于不同地区式中的常量是不同的。
[0013]
2.2.间接关系
[0014]
间接联系地震和电磁数据之间的关系是通过archie(1942)公式和gassmann(1951)方程将电磁和地震中参数通过孔隙度和饱和度进行联系，但是这种联系也不是普遍适用的，主要在于archie公式和gassmann方程也是经验公式，随着地区的变化而变化，但就目前而言是研究联合反演比较好的联系。
[0015]
对于这种关系，首先通过archie公式将电导率σ和孔隙度φ及含水饱和度s
w
相联系，其公式为
[0016][0017]
其中a是曲折因子，m为孔隙度指数，n为饱和度指数，σ
w
为含盐水结构电导率。之后通过gassmann方程将孔隙度φ和油、气、水的饱和度so，sg，sw。同地震波p波速度vp进行联系，其公式为
[0018][0019]
其中
[0020]
k
sat
＝(1-β)k
ma
+β2m
[0021][0022][0023]
ρ
sat
＝(1-φ)ρ
ma
+φ(s
w
ρ
w
+s
o
ρ
o
+s
g
ρ
o
)
[0024]
β是biot系数，通常情况下其是孔隙度的函数。
[0025]
以上我们介绍了电磁和地震数据之间的直接和间接物性关系，目前研究还只用间接关系，直接关系还没有得到应用。
[0026]
现有的隧道内地震和瞬变电磁方法存在两个突出问题：
[0027]
(1)反射地震方法能准确探测掌子面前方危险源的空间位置，但无法判定是否含有裂隙水，富水构造内是否含水等。
[0028]
(2)瞬变电磁探测方法能准确判定富水构造、裂隙水等异常体，但对异常体的空间位置不敏感，无法准确探明异常体的空间位置。瞬变电磁反演结果严重依赖于初始电阻率模型，隧道环境内因围岩视电阻率初始模型的误差，致使瞬变电磁的反演多解性较强，反演结果中富水构造带的视电阻率值和空间位置都有非常强的不确定性。


技术实现要素：

[0029]
针对上述问题，本发明提供了一种隧道地震和瞬变电磁联合探测方法，以解决隧道瞬变电磁反演成像时初始模型问题，提高掌子面前方富水异常体的反演精度。结合隧道
地震和瞬变电磁方法，提高可预报精度。
[0030]
为达到上述目的，本发明提供的具体技术方案如下：
[0031]
一种隧道地震和瞬变电磁联合探测方法，其关键技术在于，包括以下步骤：
[0032]
步骤1)：从反射地震的偏移成像数据中提取波阻抗信息，确定不同距离段的波阻抗系数值；
[0033]
步骤2)：依据不同距离段的波阻抗系数值反推计算出掌子面前方围岩的速度值；
[0034]
步骤3)：利用步骤2)计算出掌子面前方围岩的围岩速度值，依据地震速度和电阻率的关系得到的波阻抗曲线，将波阻抗曲线转换成围岩的电阻率值初始模型；
[0035]
步骤4)：根据步骤3)得到的围岩的电阻率值初始模型实现瞬变电磁数据的反演；
[0036]
所述步骤4)中，围岩的电阻率值初始模型实现瞬变电磁数据的反演的过程包括以下粒子群优化算法的基本步骤：
[0037]
步骤5)：结合隧道地震和瞬变电磁成像结果，提供合理的预测结论，提高掌子面前方异常体的可预报精度。
[0038]
通过上述设计，应用地震偏移成像数据转化提取的电阻率值初始模型，减小的电阻率值初始模型对反演结果的误差，降低瞬变电磁反演的多解性，提高瞬变电磁反演结果中对富水构造带的空间位置和含水量的反演精度。
[0039]
进一步描述，所述步骤1)中不同距离段的波阻抗系数值为：
[0040][0041]
z
k
是不同距离段的波阻抗系数值，r
k+1
为前一段的振幅值，r
k
为后一段的振幅值，h
max
为最大距离。
[0042]
再做进一步描述，步骤2)中，所述掌子面前方围岩的速度计算内容为：
[0043]
假设检波器的安装位置至掌子面距离段的波阻抗系数值为：
[0044]
zk
j j＝h
geophone
,1,2,....,h
face
；
[0045]
该段距离内平均波阻抗系数值zk
average
为：
[0046][0047]
其中，h
geophone
是检波器的安装位置，h
face
是掌子面的位置；
[0048]
所述掌子面前方围岩的速度值vrock
j
为：
[0049]
vrock
j
＝vdirect*z
k
/zk
average j＝h
face
,1,2,....,h
max
；k＝0,1,2,....,h
max
[0050]
zk
average
是平均波阻抗系数，z
i
是不同距离段的波阻抗系数，vdirect是直达波速度。
[0051]
再做进一步描述，所述步骤3)中，依据地震速度和电阻率的关系得到的波阻抗曲线，将波阻抗曲线转换成围岩的电阻率值初始模型的过程包括以下步骤：
[0052]
所述地震速度和电阻率的直接关系为：
[0053]
v＝khc
drt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)；
[0054]
根据公式(1)得到波阻抗曲线rt的公式为：
[0055][0056]
式中v为声波速度；k、c、d为待定系数，h为围岩的深度；
[0057]
将公式(2)中的所有待定系数均设置1，忽略围岩的深度h的变化，转换得到围岩的电阻率值初始模型rt'：
[0058][0059]
再做进一步描述，围岩的电阻率值初始模型实现瞬变电磁数据的反演的过程采用粒子群优化算法，算法包括以下的基本步骤：
[0060]
(1)将波阻抗曲线转换成围岩的电阻率和深度参数，建立瞬变电磁法反演初始模型进行计算；
[0061]
(2)计算反演目标函数方程：
[0062]
min φ(v,t)＝||d(v,t)-x(v,t)||2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0063]
其中d(v,t)表示野外实测感应电位数据；x(v,t)表示根据正演模型每次迭代计算的感应电位数据，定义误差的二范数为目标函数，求解最优解的过程就是当φ(v,t)取最小值时，对应输入的视电阻率值和深度参数，然后进行迭代计算，当迭代后目标函数优于迭代前时，将当前的反演参数作为最优参数，否则保留原反演参数；
[0064]
(3)判断：是否达到最高迭代次数，如否，迭代次数增加1后进入步骤(4)，如果达到误差要求范围内或达到最高迭代次数，停止迭代，输出视电阻率和深度参数；
[0065]
(4)持续更新模型的反演参数视电阻率值和深度参数，计算新的目标函数，返回步骤(2)，直到满足终止条件，停止迭代，输出视电阻率和深度参数；
[0066]
(5)根据反演得到了视电阻率值和深度参数进行成像分析。
[0067]
通过上述设计，充分结合隧道地震和瞬变电磁方法的优势，避免单一方法的局限性，提高掌子面前方异常体的可预报精度。
[0068]
本发明的有益效果：应用地震偏移成像数据转化提取的围岩的电阻率值初始模型，减小围岩的电阻率值初始模型对反演结果的误差，降低瞬变电磁反演的多解性，提高瞬变电磁反演结果中对富水构造带的空间位置和含水量的反演精度，可以很好地满足实际应用的需要。结合隧道地震和瞬变电磁成像结果，提供合理的预测结论，提高掌子面前方异常体的可预报精度。成本低、易于实现且传输速率高、可靠性好等优点，可以实现较好的灵活度和空间位置偏移冗余度。
附图说明
[0069]
图1为本发明的主流程图；
[0070]
图2为云南某隧道地震数据采集和瞬变电磁数据采集工作示意图；
[0071]
图3为云南某隧道纵波偏移剖面图；
[0072]
图4为云南某隧道电阻率值初始模型图；
[0073]
图5为云南某隧道应用平均电阻率模型和初始电阻率模型分别反演得到的电阻率模型图；
[0074]
图6为四川某隧道地震数据采集和瞬变电磁数据采集工作示意图；
[0075]
图7为四川某隧道纵波偏移剖面图；
[0076]
图8为四川某隧道电阻率值初始模型图；
[0077]
图9为四川某隧道分别应用平均电阻率模型和初始电阻率模型反演得到的电阻率模型图；
[0078]
图10为本发明的粒子群优化算法流程图。
具体实施方式
[0079]
下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。
[0080]
从图1可以看出，一种隧道地震和瞬变电磁联合探测方法，包括以下步骤：
[0081]
步骤1)：从反射地震的偏移成像数据中提取波阻抗信息，确定不同距离段的波阻抗系数值；
[0082]
其中，所述步骤1)中不同距离段的波阻抗系数值为：
[0083][0084]
z
k
是不同距离段的波阻抗系数值，r
k+1
为前一段的振幅值，r
k
为后一段的振幅值，h
max
为最大距离。
[0085]
步骤2)：依据不同距离段的波阻抗系数值反推计算出掌子面前方围岩的速度值；
[0086]
其中，步骤2)假设检波器的安装位置至掌子面距离段的波阻抗系数值为：
[0087]
zk
j
,j＝h
geophone
,1,2,....,h
face
；
[0088]
该段距离内平均波阻抗系数值zk
average
为：
[0089][0090]
其中，h
geophone
是检波器的安装位置，h
face
是掌子面的位置；
[0091]
所述掌子面前方围岩的速度值vrock
j
为：
[0092]
vrock
j
＝vdirect*z
k
/zk
average j＝h
face
,1,2,....,h
max
；k＝0,1,2,....,h
max
[0093]
zk
average
是平均波阻抗系数，z
i
是不同距离段的波阻抗系数，vdirect是直达波速度。
[0094]
步骤3)：利用步骤2)计算出掌子面前方围岩的围岩速度值，依据地震速度和电阻率的关系得到的波阻抗曲线，将波阻抗曲线转换成围岩的电阻率值初始模型；
[0095]
其中，所述步骤3)中，依据地震速度和电阻率的关系得到的波阻抗曲线，将波阻抗曲线转换成围岩的电阻率值初始模型的过程包括以下步骤：
[0096]
所述地震速度和电阻率的直接关系为：
[0097]
v＝khc
drt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)；
[0098]
根据公式(1)得到波阻抗曲线rt的公式为：
[0099][0100]
式中v为声波速度；k、c、d为待定系数，h为围岩的深度；
[0101]
将公式(2)中的所有待定系数均设置1，忽略围岩的深度h的变化，转换得到围岩的电阻率值初始模型rt'：
[0102][0103]
步骤4)：根据步骤3)得到的围岩的电阻率值初始模型实现瞬变电磁数据的反演；
[0104]
从图10可以看出，围岩的电阻率值初始模型实现瞬变电磁数据的反演的过程采用粒子群优化算法，算法包括以下的基本步骤：
[0105]
(1)将波阻抗曲线转换成围岩的电阻率和深度参数，建立瞬变电磁法反演初始模型进行计算；
[0106]
(2)计算反演目标函数方程：
[0107]
min φ(v,t)＝||d(v,t)-x(v,t)||2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0108]
其中d(v,t)表示野外实测感应电位数据；x(v,t)表示根据正演模型每次迭代计算的感应电位数据，定义误差的二范数为目标函数，求解最优解的过程就是当φ(v,t)取最小值时，对应输入的视电阻率值和深度参数，然后进行迭代计算，当迭代后目标函数优于迭代前时，将当前的反演参数作为最优参数，否则保留原反演参数；
[0109]
(3)判断：是否达到最高迭代次数，如否，迭代次数增加1后进入步骤(4)，如果达到误差要求范围内或达到最高迭代次数，停止迭代，输出视电阻率和深度参数；
[0110]
(4)持续更新模型的反演参数视电阻率值和深度参数，计算新的目标函数，返回步骤(2)，直到满足终止条件，停止迭代，输出视电阻率和深度参数；
[0111]
(5)根据反演得到了视电阻率值和深度参数进行成像分析。
[0112]
在实施例1中，课题组研究成员于2019年12月份在云南昭通市某公路隧道开展地震(tsp)见图1中a和瞬变电磁超前地质预报见图1中b。超前钻孔显示掌子面前方19.76米见涌水，涌水量约2-3m3/min，总量约5万方。掌子面岩性以中风化白云质灰岩、泥质灰岩为主。隧道内不良地质体主要为岩溶构造区、暗河以及下伏劣质煤层。岩石呈深灰色、灰白色。锤击回弹难碎，属较坚硬岩。
[0113]
结合图3、4和5可以看出，根据一种隧道地震和瞬变电磁联合探测方法的步骤1)，tsp数据的纵波偏移剖面，分别在掌子面前方19m图3中的a处、35m图3中的b处、75m图3中的c处、102m图3中的d处、处存在强阻抗界面。钻孔显示掌子面前方19.76米见涌水，与tsp偏移剖面中19m图3中的a处的强阻抗界面相一致。
[0114]
根据一种隧道地震和瞬变电磁联合探测方法的步骤3)依据地震波速与瞬变电磁的视电阻率的直接关系，提取出掌子面前方60m的电阻率初始模型图4。图3显示掌子面前方约20m和32m为相对低视电阻率，约为500ω
·
m，42m和52m为相对较高的视电阻率，约为600ω
·
m。结合步骤4)瞬变电磁数据反演时，视电阻率的反演结果非常依赖于初始模型。数据反演时，分别给定平均视电阻率值为600ω
·
m和初始电阻率模型实现瞬变电磁反演。两个剖面的反演结果差异较大，尤其是在掌子面前方40m～55m。给定平均视电阻率模型时，反演剖面中异常区域较多，发现三个低视电阻率区域，详见图5中的a中的a、b和c处。而给定初始电阻率模型图4时，通过约束反演减少异常体区域数量，在掌子面前方16m～24m图5中的b中的a处、29m～31m(图5中的b中的b处)显示低视电阻率区域，而38m～54m图5中的b中的c处显示高电阻率区域。已开挖的地质资料显示，掌子面前方19m为富水溶洞，与瞬变电磁反演结
果中a区域相吻合图5中的b。
[0115]
实施例2，课题组成员于2019年12月份在四川某隧道开展tsp和瞬变电磁超前地质预报工作。折多山隧道为青藏高原区域内典型的极高风险隧道，海拔高度在3600m～4600m。隧道区域内不良地质主要表现为活动断裂、崩塌与危岩、岩爆、大变形、涌突水、瓦斯、强震区、积雪冻土等。正在开挖的隧道区段为山前洪积物，碎石和破碎带比较发育，掌子面未见明显的漏水。tsp图6中的a和瞬变电磁图6中的b探测的现场工作照。
[0116]
图7为tsp数据的纵波偏移剖面，分别在掌子面前方12m图7中的a处、25m图7中的b处、49m图7中的c处和70m图7中的d处发现强阻抗界面。后期开挖验证掌子面前方52m处发现大量的破碎岩石从掌子面上方倾泻而下，与49m图7中的a处的强阻抗界面相吻合。
[0117]
从图7所示的纵波偏移剖面，依据地震波速与瞬变电磁的视电阻率的直接关系，提取出掌子面前方100m的初始电阻率模型图8。图8显示掌子面前方约0m～5m段为较高的视电阻率，约为500ω
·
m，而在50m～100m段为相对较低的视电阻率，约为350ω
·
m。
[0118]
瞬变电磁数据反演时，视电阻率的反演结果非常依赖于初始模型。数据反演时，分别给定平均视电阻率值为400ω
·
m和初始电阻率模型图8实现瞬变电磁反演图9。两个剖面的反演结果差异较大，给定平均视电阻率模型时，反演结果的分辨率较低，并且整个视电阻率值都偏低，约为140ω
·
m，而在掌子面前方85m～95m为较高的视电阻率值，约为220ω
·
m图9中的a。给定初始电阻率模型图8时，反演结果的异常区域较为集中，在掌子面前方48m～60m图9中的b中的a处为相对较低视电阻率，约为300ω
·
m。已开挖的地质资料显示，掌子面前方52m出现极为破碎的岩石，与瞬变电磁反演结果中a区域相吻合图9中的b。
[0119]
本发明提供的隧道地震波阻抗约束的瞬变电磁反演方法，具有以下优点：应用地震偏移成像数据转化提取的电阻率初始模型，减小初始模型对反演结果的误差，降低瞬变电磁反演的多解性，提高瞬变电磁反演结果中对富水构造带的空间位置和含水量的反演精度，可以很好地满足实际应用的需要。
[0120]
以上所述实施例仅表达了本发明的实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。