《1、 引言》
1、 引言
随着全球经济的发展和现代化建设需求的增加，隧道已成为复杂地质条件下交通基础设施建设的最佳选择之一。隧道施工环境地质条件复杂，各种地质灾害，如塌方或突水，通常是由岩溶和断裂带等不良地质体引起的；由于隧道施工发生在地质体中，如果不能准确预测隧道周围的不良地质体（带），隧道施工可能会面临大量人员伤亡和财产损失[1]。因此开展隧道超前探测工作，准确、及时地探明掌子面前方的隐伏地质灾害，确保施工现场安全，避免突发隧道地质灾害至关重要。
隧道超前探测是一种在全空间条件下在掌子面或迎头后方布置观测系统，超前探测前方隐伏地质灾害的地球物理技术。目前用于隧道超前探测的地球物理方法主要包括地震类和电磁法类[2‒5]，其中地震波超前预报方法探测距离长，对大型地质体预报较准确，且具有较好的地质结构面识别效果，已成为隧道内常用的超前探测方法。20世纪70年代末，德国和英国工程师利用槽波的埃里震相（airy seismic phase）进行巷道工作面前方的地质构造探测；20世纪90年代初，瑞士勘测技术有限公司（Swiss Surveying Technology Co., Ltd.）开发了一套完整的隧道地震预测系统——隧道地震超前预报（TSP）系统；20世纪90年代末，一家美国工程公司开发了真反射层析成像（TRT）技术；21世纪初与Inazaki等[7]提出了隧道垂直地震剖面法，然后由Alimoradi等[9]和Zhao等[8]在21世纪初开发应用。基于隧道建设需求的增加，到目前地震类超前地质预报方法已发展了包含隧道地震超前预报、水平地震剖面（HSP）、真反射层析成像、隧道地震层析成像（TST）、随钻隧道地震（TSWD）、隧道地质预测（TGP）和地震波负视速度法等多个分支[10‒15]。其中隧道地震超前预报基于异常地质体和围岩之间地震波波速的差异，激发地震波传播到隧道围岩中，对采集的地震数据进行处理，获得隧道掌子面前方异常地质体的分布和岩石力学参数，进而在隧道施工过程中进行预警和指导。但由于隧道地震超前预报观测系统受隧道环境限制，采集的数据量无法满足对隧道地质体波速的高计算精度要求，采用超前预报数据处理方法得到的成像结果无法达到对不良地质体精细探测的需求，只适用于简单地质情况[16]，因此亟需采用适用于隧道超前地质预报的精细探测手段来改变这一现状。
20世纪80年代，Poirier和Lailly等[17‒19]提出了基于最小二乘法的时域全波形反演（FWI），并将这一概念引入地震勘探领域。相比常规利用单一反射波或初至波数据处理方法进行物性参数反演成像，全波形反演充分使用了全波场的信息，达到了更高的分辨率，因此该方法在地震波场反演重建中备受推崇，并在地震速度建模研究与应用中受到了越来越广泛的关注[20‒22]。由于全波形反演的非线性和周期跳跃，目标函数具有多个局部极小值，这使得反演明显依赖于初始模型[23‒25]。为了减少反演结果对初始模型的依赖性，早期研究人员提出了一种时域多尺度全波形反演方法，对地震数据进行滤波，将低频反演结果作为高频反演的初始模型进行不同尺度的频率反演[26]。20世纪90年代，Pratt等[27]将全波形反演的理论扩展到频率域，称为频率域全波形反演方法。频率域全波形反演中低频分量的反演结果可直接作为高频分量的初始模型，直接实现多尺度反演的效果，减少对初始模型的依赖[27‒31]。由于这一优点，频率域全波形反演在地震勘探中得到了广泛应用。
在隧道超前探测方面，Musayev等[32‒33]首先将频率域中的全波形反演方法应用于隧道，并讨论了全波形反演方法能否进行隧道速度场的反演。Nguyen等[34‒35]提出了二维（2D）隧道地震波全波形反演的全局优化程序，实现了对隧道前方扰动区的高精度识别。Bharadwaj等[36]开发了一种地震预测系统，可以在隧道掘进机前方直接成像。Lamert等[37]提出了适用于隧道空间的弹性波时域全波形反演方法，预测隧道掌子面前方的扰动区域。国内研究中，Zhang等[38]使用基于全波形反演方法的探地雷达（GPR）来识别隧道掌子面前方20 m范围内的不良地质体；使用声波全波形反演方法预测隧道掌子面前方大型地质体的速度界面；Feng等[40]利用全波形反演探地雷达剖面改进了隧道衬砌缺陷的重建。
本文基于隧道空间狭小、观测系统受限的特点，针对目前隧道超前探测中探测精度不足的问题，将频率域声波全波形反演方法应用于隧道超前探测。以隧道低速异常体为例，构建了隧道低速异常体模型和基于隧道空间的观测系统；采用频率域声波全波形反演进行隧道速度波场重建；通过对比频率域多尺度全波形反演的不同频率组选择策略，讨论适用于隧道环境地震方法的频率组选择策略；同时讨论了隧道超前探测中观测系统侧边长度对反演效果的影响；最后采用复杂模型验证了该方法以及参数选择策略的有效性。
《2、 频率域多尺度全波形反演原理》
2、 频率域多尺度全波形反演原理
《2.1 频率域声波全波形反演原理》
2.1 频率域声波全波形反演原理
各向同性介质中的二维频率域声波波动方程公式如下[41]：
（1）
式中，k(x, z)为体变模量；ρ(x, z)为密度；d(x, z, ω)和s(x, z, ω)分别为波场压力函数与震源函数；(x, z)为空间位置向量；ω表示频率。
波场压力函数d(x, z, ω)与震源函数s(x, z, ω)的关系是线性的，因此离散化的二维声波方程等式可以简化为如下的大型稀疏线性方程组：
（2）
式中，是关于频率和介质属性的阻抗矩阵。考虑模型空间网格数为(
选择观测波场与计算波场差值的最小二范数作为频率域声波全波形反演的目标函数，具体表示为：
（3）
式中，
（4）
式中，T是雅可比矩阵的转置，该矩阵是波场对模型参数的偏导；∆d*是波场残差的复共轭；Re是复数的实部。在迭代中加入不同的模型参数扰动，通过计算E(m)的梯度函数从而得到最小化观测和计算波场之间差值的模型参数[44]。将反演误差阈值和迭代次数作为模型更新的终止条件，以保证以实际代价收敛。当反演终止时，得到满足条件的模型介质参数。
《2.2 频率域多尺度反演频率选择策略》
2.2 频率域多尺度反演频率选择策略
《2.2.1. 多尺度频率选择策略》
2.2.1. 多尺度频率选择策略
频率域多尺度反演方法可以相对提髙全波形反演对初值的稳健性，总体上是先从低频出发开始反演，然后再逐级提高频率进行反演。由于因低频数据反演陷入局部极值而导致不收敛的概率较小，因此可以估计相对较好的初始模型，提高反演过程的稳定收敛来获取更合理的解。频率域全波形反演通过低频反演得到大尺度的信息，中高频反演信息刻画细节特征，采用低频反演结果信息作为后续中高频率反演的初始模型，直接实现多尺度反演[45]。
因此，在频率域全波形反演中对频率组的选择值得研究，Sirgue和Pratt [46]基于二维频率-波数域的关系，提出若要满足垂直波数覆盖则需满足当前频率值对应的最大波数与下一频率对应的最小波数相等：
（5）
式中，
（6）
式中，α即为最大入射角的余弦函数值；hmax是观测系统最大半偏移距；z为反射层深度。
这一策略在大偏移距观测系统全波形反演中可较好实现。然而除偏移距外，这一策略还要求高密度的检波器分布和横向速度变化较小的模型条件，因此具有一定的局限性。隧道地震超前探测观测系统偏移距小，而且地质情况较复杂，频率范围大，采用该公式得到的频率组不仅限制成像深度更会影响成像效果。Wang等[47‒48]基于观测系统采样率不足的情况，依据抗混叠条件下垂直波数Δ
（7）
式中，c0为背景速度；Δf表示相邻两频率值间隔；zmax为最大成像深度。
这一公式给出符合频率组间隔要求的等分条件，这种频率采样策略已被证明具有更广泛的应用可能性；尽管由于低频区域采样点的减少，反演的收敛速度相对较慢。
《2.2.2. 隧道空间的频率域多尺度频率选择策略》
2.2.2. 隧道空间的频率域多尺度频率选择策略
隧道地震超前探测方法观测系统照明角度小、主频高、频率分布范围大，常规的频率组选择策略在隧道全波形反演中不再适用。针对这些情况，本文考虑综合两种频率选择策略进行组合的频率组采样策略。
（8）
式中，α是最大入射角的余弦函数值；c0是背景速度；Δf表示相邻两频率值间隔。此外，对于弹性波全波形反演，建议使用横波速度为背景速度。
公式（8）采用公式（6）得到覆盖垂直波数的低频加密频率组（图1），保证初始模型反演大尺度上的收敛方向。当由公式（6）得到的频率组无法满足公式（7）的抗混叠波数采样原则时，采用公式（7）得到抗混叠的线性增加高频频率组，在反演优化过程中得到局部最优值，保证对反演模型的细节刻画。
《图1》
图1 组合频率组选择策略曲线示意图。在低频范围内，频率间隔Δf随频率组数线性增加，当Δf无法满足抗混叠条件时，频率间隔选择为满足抗混叠的等值Δf。
《3、 参数对比和结果讨论》
3、 参数对比和结果讨论
《3.1 模型与观测系统设计》
3.1 模型与观测系统设计
地震隧道超前探测方法是在隧道掌子面及边墙内布置深度在1 m左右的激发孔，在孔内通过人工激发地震波向隧道围岩中传播，当围岩波阻抗发生变化时，一部分地震波会被反射回来，另一部分地震波将会继续向前传播。反射回的地震波由接收器接收并传回电脑主机形成地震波记录，处理后的记录可预报隧道掌子面前方的地质变化，为隧道施工提供可靠的地质资料[49‒50]。基于实际隧道参数，本文设计图2（a）所示的隧道低速异常体模型，模型大小为200 m（X轴方向）× 30 m（Z轴方向），隧道长100 m，掌子面宽12 m，网格步长1 m；隧道内部填充空气，隧道围岩速度为4000 m·s-1，在掌子面前方14 m、46 m、84 m深度处分别包含三个半径为3 m的圆形低速异常体，速度为3000 m·s-1。图2（b）为观测系统示意图，包含掌子面和垂直于掌子面的已开挖侧边，范围为51~101 m，震源点间距2 m，检波点间距1 m，激发深度为1 m。采用200 Hz Ricker子波作为频率域正演模拟的主频。图2（c）~（e）表示50 Hz、200 Hz和450 Hz的单频压力场切片。如图所示，低频对应长波长反应模型大尺度的信息，当速度保持不变且频率增加时，单频波场的波长变短，低速异常体与围岩波长之间的差异增大，关于异常体位置的细节变得更加明显。
《图2》
图2 隧道模型与观测系统设计以及频率域单频压力场。（a）隧道低速异常体模型，其中X是隧道模型长度，Z代表宽度；（b）隧道U型观测系统示意图；（c）~（e）震源位于掌子面中部的50 Hz、200 Hz、450 Hz的单频波场。
采样间隔选择为0.5 ms，记录长度为160 ms，通过傅里叶逆变换得到时间域声波地震记录如图3所示。图3（a）中震源点位于坐标(51, 9)处，是隧道侧边观测系统的起点。图3（b）中的炮点位于隧道掌子面中部，坐标为(102, 15)。对时间域地震记录进行分析，这些记录包含直达波和三个低速异常体界面引起的绕射波，记录在隧道拐角位置出现了波形和能量的变化。
《图3》
图3 隧道异常体模型的声波地震记录。（a）震源点位于隧道侧边的时间域声波地震记录；（b）震源点位于掌子面中部的时间域声波记录。
《3.2 不同频率选择策略全波形反演讨论》
3.2 不同频率选择策略全波形反演讨论
频率域多尺度反演原理是自低频到高频将前一个频率反演结果作为下一个频率正演的初始模型依次进行反演。由于低频反演提供长波信息，而高频反演提供详细的波传播信息，因此多尺度全波形反演收敛稳定并刻画细节[51]。采用图2（b）所示的观测系统，对图2（a）所示的隧道异常体模型进行频率域声波全波形反演。反演过程中采用速度为3000 m·s-1的均匀速度模型作为初始模型，正演模拟同理论模型一致选择主频为200 Hz的Ricker子波。根据Ricker子波的频谱，选择由低到高三个频率（50 Hz、200 Hz、450 Hz），依次反演得到均方根误差（RMSE）收敛曲线和速度模型，如图4所示。
《图4》
图4 隧道异常体模型的声波地震记录。（a）50 Hz、200 Hz、450 Hz三个频率依次反演的迭代曲线；（b）~（d）50 Hz、50 Hz和200 Hz以及50 Hz、200 Hz和450 Hz频率的反演结果。（a）中均方根误差曲线展示了迭代中三个频率的收敛过程，（b）、（d）描述了从低频到高频的不同尺度反演结果，展示了反演结果从大尺度信息到异常体位置的精细刻画。
基于公式（6）~（8）并考虑数据采集中最大炮检距与勘探深度关系，依据子波频谱得到介于10 Hz到500 Hz的三组频率组，三组频率分别对应公式（6）的S频率组、公式（7）的W频率组以及公式（8）的C频率组。子波频谱及频率组参数曲线如图5（a）、（b）。
《图5》
图5 声波全波形反演参数。（a）震源子波频谱；（b）基于公式（6）~（8）得到的三组频率组参数曲线分别为W频率组、S频率组、C频率组。
以图5（b）得到的三组频率作为多尺度反演迭代的频率组，分别对图6（a）所示背景速度初始模型进行除频率组外其他参数一致的频率域声波全波形反演，得到图6（b）~（d）不同频率选择策略的频率域声波全波形反演速度成像结果。与图4（d）中三个频率反演的速度结果相比，图6（b）~（d）所示的反演结果速度分辨率明显提高，而三组不同频率组反演结果的具体速度剖面也存在差异。
《图6》
图6 初始模型和不同频率组的反演结果模型。（a）用于全波形反演的初始模型；（b）~（d）基于公式（6）~（8）得到的W、S、C频率组反演结果速度模型；（e）过速度模型坐标轴Z = 15的速度剖面曲线。反演结果揭示了不同频率组选择策略在不同深度异常体反演中的反演分辨率差异。
从两个方面分析反演成像的效果：一是反演结果在大尺度上对异常体位置的敏感性；二是反演结果对异常体速度刻画的准确性。首先，比较图6（b）~（d）的反演结果对异常体位置反演的敏感性：对于浅部地层异常体位置反演效果，图6（b）所示W频率组的反演相对准确，W频率组第一个异常体位置存在明显的速度变化，但中、深部异常体位置相对较差，尤其位于隧道掌子面前方84 m深度处的异常体，可见该区速度变化范围明显大于异常体所在范围；图6（c）所示的S频率组反演结果对异常体前缘位置的反演效果更准确，速度曲线在异常体位置与相邻位置之间有较明显的速度变化，说明反演过程准确收敛于该位置，这也符合S频率组在低频范围内覆盖垂直波数，因而保证了大尺度上收敛方向的准确性；图6（d）所示C频率组对异常体位置的反演效应在浅、中、深部地层位置都非常稳定，中剖面速度曲线在异常体位置有明显的速度变化，反演收敛位置准确。对比可见图6（c）的S频率组和图6（d）的C频率组反演结果对异常体位置的敏感度较高，较能准确判断最优化反演的收敛方向。
其次通过图6（e）过异常体纵剖面和图7（a）~（c）过不同深度异常体剖面的速度曲线对比，讨论频率组选择对异常体速度刻画的准确性。对比可见三组频率的反演结果对浅部地层异常体的速度反演精度都较高，对中、深部地层异常体的速度反演精度都比浅部地层的精度差, 总体规律表现为反演位置的深度越深，成像结果精度越差。在三组频率单独对比中，S频率组对速度的反演精确度相对较差；W频率组以及C频率组的速度反演精度较接近，尤其是浅部地层的反演结果与初始模型差异极小，在中部和深部的反演结果中C频率组体现了更高的分辨率。
《图7》
图7 三组反演结果不同异常体速度对比曲线。（a）过速度模型坐标轴X = 114的第一个异常体的速度曲线；（b）过速度模型坐标轴X = 147的第二个异常体的速度曲线；（c）过速度模型坐标轴X = 186的第三个异常体的速度曲线。
综上所述，可知基于隧道空间的频率域声波多尺度全波形反演可以得到掌子面前方异常体成像，并且本研究提出的组合频率组选择策略在频率域反演中表现更好，能得到深部分辨率更高的速度反演结果。
《3.3 观测系统设计影响》
3.3 观测系统设计影响
在基于隧道超前探测的全波形反演理论模拟中，为增加对隧道前方地质体有效信息的接收，隧道地震超前预报观测系统采用了U型观测系统。参考地震波在隧道空间和介质的传播路径，以及地震波传播过程中的能量损失，本研究认为隧道两侧观测系统的扩展范围在达到一定长度后将不再影响反演结果。因此设计5组线性增加的不同侧边长度（10 m、30 m、50 m、70 m、90 m）U形隧道观测系统，对采用公式（8）得到的频率组进行频率域多尺度声波全波形反演。得到5组反演结果如图8（b）~（f），具体的速度剖面结果在图9中显示。
《图8》
图8 初始模型和不同侧边长度反演结果速度模型。（a）用于反演的初始速度模型；（b）~（f）分别为10 m、30 m、50 m、70 m和90 m的隧道观测系统侧边长度反演结果速度模型。可知不同观测系统下得到的反演速度结果及不同深度异常体在反演结果上的差异。
《图9》
图9 不同反演结果速度对比曲线。（a）过观测系统侧边长度为10 m、30 m、50 m反演结果速度模型坐标轴Z = 15的速度剖面对比曲线；（b）过观测系统侧边长度为10 m、30 m、50 m反演结果速度模型坐标轴Z = 15的速度剖面对比曲线。
对比图8（b）~（d）和图9（a）可知，随着双侧观测系统排列长度的增加，反演成像结果精度增加明显，其中侧边长度为10 m时对深部无法成像，侧边长度为30 m时深度异常体位置开始出现速度变化，侧边长度为50 m时对深部异常体能有分辨率较高的成像效果。对比图8（d）~（f），侧边观测系统的长度变化超过一定范围后对反演的成像精度影响趋于减小，到达一定范围后，增加侧边观测系统长度对反演成像结果无明显影响。通过图9（b）异常体横轴中线速度曲线对比，可知70 m侧边长度的观测系统的速度反演结果与90 m侧边长度的反演结果差异极小，相对误差小于0.4%。对于该模型隧道的掌子面宽度，考虑内置检波器及震源深度，选择侧边长度为70 m的观测系统即可得到相对最优的成像结果。
本文通过多次对比和摸拟，综合得出一定范围内，侧边观测系统排列长度越长，反演结果的分辨率越高。当侧边观测系统长度大于掌子面宽度5倍后，继续增加侧面观测系统对反演结果影响可忽略不计。这一结论可为隧道超前探测的实际施工以及后续隧道全波形反演的理论研究提供建议。
《3.4 复杂模型试算及应用》
3.4 复杂模型试算及应用
为验证前文讨论的参数在隧道超前探测频率域声波全波形反演的效果，本文设计建立了如图10（a）所示包含异常体、低速带以及复杂构造形态的隧道地层模型，该模型以波速为4500~5500 m·s-1的灰岩岩层为背景，包含速度为3000 ~4000 m·s-1的低速薄层异常体地层等复杂地质条件。用于反演中的初始模型采用真实模型的模糊速度，如图10（b）所示。采用三种频率选择策略作为频率组参数选择依据，使用侧边长为70 m的隧道观测系统作为隧道观测系统，最终得到如图10（c）~（e）所示速度反演结果。
《图10》
图10 复杂模型全波形反演结果对比。（a）复杂模型示意图；（b）用于复杂模型全波形反演的初始模型；（c）~（e）基于公式（6）得到的W频率组、基于公式（7）得到的S频率组以及基于公式（8）得到的C频率组反演结果速度模型。浅部地层三组反演结果与真实模型之间差异较小，在中、深部地层反演中图（e）的反演结果速度表现更好。
图10（c）~（e）反演速度成像结果显示，频率域的多尺度全波形反演能对隧道空间复杂地层中地层所在的位置和速度进行反演，其中对于越接近掌子面的地层信息刻画越准确，整体成像较准确，对深部的地层细节刻画相对浅部地层较差。本文所提出的基于公式（8）的组合频率组反演结果速度模型在中、深部地层反演中表现更好，分辨率优于其他两组频率组的反演结果。进一步对组合型频率组反演结果进行速度对比，如图11所示。比较图中的剖面速度[图11（d）]，可见由反演结果得到的整体速度归位比较准确，对地层位置的识别精准，在深部地层速度相对起伏明显的位置存在一些误差。反演结果与实际模型在不同位置的相对误差为0.3%~8%，从浅到深误差呈增加趋势，距离掌子面较近位置误差最小，最远端误差为6%左右。
《图11》
图11 C频率组反演结果速度对比。（a）~（c）初始模型、真实模型及反演结果速度模型；（d）坐标轴Z = 15的反演结果速度模型速度剖面对比曲线。
针对隧道超前探测的复杂模型多尺度频率域声波全波形反演结果证明，在合适参数条件下频率域全波形的方法可以成功反演隧道掌子面前方地层信息，并得到关于掌子面前方复杂地层信息的分辨率较高的成像结果。
所有的理论研究和讨论都是为了向实际应用发展。基于此，使用如图12所示隧道超前探测观测系统和本研究中讨论的参数进行简单的现场全波形反演验证。
《图12》
图12 野外隧道超前探测观测系统示意图。
对采集的时间域地震数据X分量（隧道长轴方向）进行滤波以及波场分离等操作后得到声波记录的分量，然后进行傅里叶变换，得到频率域后用于频率域的全波形反演。首先对隧道超前探测的X分量记录[图13（a）]进行滤波和波场分离，获得P波分量记录[图13（b）] [52]，采用如图13（c）所示直达波速度初始模型进行组合频率组的频率域声波全波形反演，得到反演速度结果模型如图13（d）所示。
《图13》
图13 隧道地震记录和速度模型。（a）接收器R1采集的X分量隧道记录；（b）处理后的对应X分量记录；（c）用于反演的初始速度模型；（d）采用组合的频率组反演结果速度模型。
根据图13（d）中的反演结果速度模型可知，低速区位于距离隧道掌子面35~45 m处，速度为3600 m·s-1左右。在随后的开挖过程中，在隧道掌子面前方约30 m处遇到速度为3700 m·s-1的软弱层。这一结果证明了本文提出的参数在进行隧道超前勘探的全波形反演在实际应用中是有发展潜力的，可以为今后的隧道建设工程提供指导，进而有助于安全高效地进行隧道工程建设。
《4、 结论》
4、 结论
本文将二维频率域声波全波形反演方法应用于隧道地震超前探测中，完成隧道超前探测频率域声波全波形反演，并通过对比讨论频率组选择策略以及观测系统对反演效果的影响，具体分析总结如下：
（1）由于隧道空间狭小，观测系统受限，而隧道地震记录具有主频较高，频率分布较广的特点，因此常规频率组选择策略不再适用。本文提出符合隧道地震方法超前探测的组合型频率组选择策略，该策略将覆盖垂直波数的低频选择策略和抗混叠的高频选择策略相结合。通过与传统频率选择策略进行对比，这一组合频率组选择策略在隧道空间的全波形反演中可以获得更好的成像分辨率。
（2）隧道探测中观测系统严重受限，因此向隧道掌子面两侧边界延伸的U型观测系统在理论模拟和实际应用中被广泛使用，但侧边长度的选择策略对隧道全波形反演结果的影响一直未被阐明。因此本文通过线性增加隧道观测系统侧边长度，进行多次对比和摸拟，得出一定范围内，侧边观测系统长度越长，反演结果分辨率越高；当侧边观测系统长度大于掌子面宽度5倍后，继续增加侧面观测系统长度对反演结果影响可忽略不计。
（3）对隧道复杂模型的频率域声波全波形反演表明，在选取合适参数条件下，该方法可以成功地反演隧道前方的地层信息，获得复杂地质信息的高分辨率结果。隧道超前探测野外采集数据的全波形反演结果进一步验证了隧道超前探测的全波形反演方法参数选择策略的有效性，以及该方法向实际应用方向发展的潜力。鉴于以上结论，本研究可为今后该方向的理论研究和实际应用提供思路和参考。
京公网安备 11010502051620号
