几种初至自动拾取算法在地震波法隧道超前探测中的对比分析

王 伟， 高 星， 刘孝卫

(中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室， 北京 100101)

0 引言

隧道超前地质预报进行地震数据处理时，地震波初至拾取精度会直接影响隧道掘进前方异常地质体预测结果的准确度[1]，这就需要选择初至拾取精度高的算法来进行拾取，而目前的地震初至波自动拾取的方法有很多。

在20世纪70年代以前，初至的拾取主要是靠人工根据波形特征及经验来完成的，过程枯燥、耗时且拾取结果容易因为拾取者的主观因素出现偏差。在20世纪80年代至90年代末，随着计算机的普及和相关技术的发展，大量的处理软件实现了在窗口显示地震波形，通过人机交互的方式来进行初至的拾取，拾取效率得到了很大的提高，但这并没有解决耗时和主观因素干扰的问题。进入21世纪后，大量的初至自动拾取算法被提出并实现应用，包括Chen 和 Stewart于2005年提出的长短时窗算法、Sabbiobe 和 Velis于2010年改进的能量比值计算公式以及神经网络、小波变换等算法，极大地促进了初至自动拾取算法的发展。

事实上，地震波初至拾取是一项繁重复杂但非常重要的工作，现阶段还没有一种算法可以完美地解决如噪声干扰大、信噪比低等环境下地震数据的初至拾取精度低的问题，至今研究人员仍在寻找精细准确的自动拾取方法[2]。

本文所使用的数据来源于TSP-SK(tunnel seismic prediction)型隧道地质超前预报仪，如图1所示。

(a) 隧道地质超前预报仪

(b) 工作示意图

Fig. 1 TSP-SK tunnel advance forecaster and its working sketch

该仪器运行稳定，数据记录清晰，能够满足本次初至拾取对比分析对数据质量的要求。目前，隧道超前探测常用的初至自动拾取的方法主要有4种，分别为最大振幅法(TST，tunnel scatter technology使用的方法为振幅比法)、能量比法(TSP203、TSP303使用该方法)、STA/LTA(short time average /long time average，短时间的均值/长时间的均值)法、AIC(akaike information criterino)法。初至拾取的方法虽多，但其检测震相初至的理论依据却相似，即根据信号与噪音以及各震相在某领域(时域、频域或时频域)或某些属性(运动学或动力学)方面的差异来区分信号和噪声，进而得到震相初至[3-4]。由于各个算法依据的地震属性有很大的差异，因此对各初至拾取算法的拾取结果进行对比是非常必要的。本文通过对常用的4种初至拾取的算法进行对比，分析各种算法的拾取精度，力求在现阶段找到一种简单、高效的初至拾取算法。

1 初至拾取方法

1.1 最大振幅法

最大振幅法(最大能量法)，即选取地震道上能量最大的点作为地震波初至。优点是算法操作简单，计算速度快；缺点是受噪声影响大，且容易产生错误[5]。

1.2 能量比法

该方法首先选取某一特定的地震属性道，如振幅绝对值、波形长度、振幅包络等，然后根据相邻前、后窗口的属性值比值特征确定地震波初至的位置。选取合适的地震属性道是该方法成功的关键[6-7]。

本文利用能量比法拾取地震波初至时，采用变换时窗统计法记录地震初至前后时窗内的能量差异特征，其中时窗为固定大小的时间量[8-9]。能量比法的基本公式为：

(1)

(2)

式(1)—(2)中：A′为时窗内后、前能量比；x(t)为地震记录振幅值；T1为时窗起点；T0为时窗中点；T2为时窗终点；A为1个地震道的相对能量；α为稳定系数，本文取值为3[10-11]。

利用该原理，提取能量比值函数A的最大值出现的位置作为初至到时拾取点[12]，如图2所示。

图2 第2道地震记录能量比法初至到时拾取图

Fig. 2 No.2 seismic record′s first arrival time pick-up by energy ratio method

1.3 AIC法

该方法的基本原理是求取地震信号AIC函数局部最小值的位置[13-14]。对于长度为N的信号x，其AIC函数可表示为：

AIC(k)=klog{var(x[1,k])}+(N-k-1)log{var(x[k+1,N])}。

(3)

式中： AIC即峰度赤池信息量准则，由日本统计学家赤池弘次创立和发展； var表示序列的方差函数。

AIC函数在初至点峰值尖锐，准确率较高；但AIC函数受时窗的影响较大，在不同的时窗下，其局部最小值出现的位置往往不同[7,15]。

利用该原理，在初步选定的初至到时点附近建立1个新的时窗，计算出这个时窗内的AIC函数值，并提取能量函数AIC的最小值出现的位置作为初至到时拾取点[3,16]，如图3所示。

图3 第2道地震记录AIC法初至到时拾取图

Fig. 3 No.2 seismic record′s first arrival time pick-up by AIC method

1.4 STA/LTA法

Allen于1978年提出地震波STA/LTA 法[17]，其主要原理是根据地震波形特征函数的长短时均比值等特征拾取初至[18-19]。当信号到达时，STA要比LTA变化更快，相应的STA/LTA值会出现明显的增加，当比值大于设定的阈值时，就可判定初至波到达[20-21]。该方法容易受到时窗长度和滑动步长的影响[22-23]。

(4)

(5)

(6)

式(4)—(6)中： STA、LTA为短、长时窗的振幅均值；B为信号的振幅;X为阈值，经过不断地对比拾取分析，本文取值为11[24-25]。

利用该原理，设计滑动时窗并依次计算STA/LTA的值，短时窗的值为固定值，而长时窗的值不断依次累加，这样可以减少突变点对初至提取的干扰。计算的同时不断判断STA/LTA的值是否超过阈值，判断流程如图4所示。若该值大于阈值X(X=1，取该阈值是通过对数据多次运算分析得到)，则将该点作为初至到时拾取点，如图5所示[1,26]。

图4 STA/LTA算法流程图

图5 第2道地震记录STA/LTA法初至到时拾取图

Fig. 5 No.2 seismic record′s first arrival time pick-up by STA/LTA method

2 对比分析

2.1 数据采集

数据来源于云南某一隧道超前探测实例，采用TSP-SK三分量加速度检波器单点接收、多点激发的方式，检波点远离隧道掌子面，采用1—19号有效激发数据，单炮药量70 g，炮间距1.5 m，最小炮检距20 m，0.062 5 ms采样，采样点数为7 218[27]。现场布置如图1(b)所示，有效激发19炮记录的x方向分量地震数据如图6所示。由图6可以看出： 各个道地震记录88 ms之前的记录占主要成分的为地震波能量，88 ms后的记录参杂了较多成分的高频声波，这部分为干扰波，必须在后续的处理中予以消除；第11道的地震记录受到的干扰比较大，因此后续分析应适当减少各个初至拾取方法对第11道数据初至拾取结果的参考权重；初至到时从第1道到第19道是逐渐减小的，大体上呈线性减小，因此在确定第1道的初至后，可以确定在首道初至到时的时窗长度内基本可以拾取到所有地震记录道的初至到时，这样可以大大减少算法的数据计算量，提高拾取效率。

图6 TSP-SKA地震记录图

2.2 数据预处理

2.2.1 地震记录去噪处理

由图6可以看出，由于爆炸声波或者环境噪音等干扰因素使得波形记录中参杂了许多噪声。本文利用巴特沃兹带通滤波器，滤去主频带范围外的噪声，即将80～1 000 Hz内的数据保留，有效消除由于噪声引起的高低频干扰。滤波后的部分道数据记录如图7和图8所示，可以看出，经过滤波去除了高频干扰信号。

(a) 地震记录原始图

(b) 地震记录滤波图

2.2.2 提取特征函数

提取特征函数是描述地震波某类特征差异的函数，可以用CFi表示。不同学者提出了不同形式的特征函数，目的是为了反映地震波到达时信号的变化规律[28-29]。本文使用的特征函数为：

(7)

(a) 地震记录原始图

(b) 地震记录滤波图

选用此特征函数的原因是，其既可以反映波形信息的振幅变化，又可以反映频率变化。推导过程如下。

设

x(n)=Bcos[ω(n)Δt+Φ]，n=1,2,3,…。

(8)

式中：ω为圆频率；Φ为相位角。

则

x(n)=Bcos[ω(n)Δt+Φ]；

(9)

x(n-1)=Bcos[ω(n-1)Δt+Φ]；

(10)

x(n+1)=Bcos[ω(n+1)Δt+Φ]。

(11)

根据三角恒等式：

cos(2α)=2cos2α-1=1-2sin2α。

(12)

将式(10)和式(11)相乘得：

x(n-1)·x(n+1)=B2cos2[ω(n)Δt+Φ]-B2sin2ω。

(13)

将式(9)带入式(13)可得：

B2sin2ω=x(n)2-x(n-1)·x(n+1)。

(14)

部分道数据的特征函数如图9和图10所示，通过与图7和图8对比可知，特征函数可以更加突出、准确地反映出地震记录波形的各个属性的变化，对提高初至拾取精度效果非常明显，论证如下。

图9 第8道特征函数图

图10 第16道特征函数图

为了方便，本文选STA/LTA法对特征函数的作用予以对比说明，采用10道地震记录数据。首先使用滤波后的波形数据作为能量比法的基础分析数据进行初至拾取，再利用特征函数作为基础分析数据进行初至拾取，可以得到表1中的结果。经计算可以得到，以原始数据为基础数据的初至拾取平均差值为0.547 ms，而以特征函数为基础的平均差值只有0.292 ms，因此使用特征函数对于提高初至拾取算法的精度有非常明显的作用。

表1以原始记录为基础和以特征函数为基础的初至拾取到时对比表

Table 1 First arrival time pick-up comparison based on original record and characteristic function ms

2.3 4种初至波拾取方法对比

利用MATLAB数据处理平台，结合4种算法的原理编写其自动拾取函数，并利用4种不同的地震波初至自动拾取算法对同一组隧道超前探测地震记录事件进行初至拾取，可以得到19道地震数据的初至到时，见表2。

由表2可以得知，最大振幅法拾取对信噪比的要求很高，容易受到环境杂波的干扰，在下文的对比中将不再对最大振幅法进行过多的分析。为了保障手动拾取初至到时结果的科学性，通过不同人进行10多次的拾取，然后取其平均值，这样可以较大程度地减少人为因素的干扰。

表2 地震波初至到时拾取表

不同算法拾取结果的折线统计图见图11。

图11 不同算法拾取结果的折线统计图

Fig. 11 Polyline statistical graph of different pick-up consequences

由图11可以看出，AIC法拾取的初至到时大部分较手动拾取的到时提前，能量比法相对于手动拾取结果也大部分提前，而STA/LTA法的拾取结果与手动拾取结果的波形较为一致。 2—15道STA/LTA法与手动拾取的结果极为接近；第1道AIC法和能量比法比较接近手动拾取结果；16—19道几种方法的初至拾取点均比手动拾取的点要延后，其主要原因之一是炮检距缩小使得算法空间较小，从而影响拾取精度。从折线图可以初步判断，STA/LTA法对于超前探测地震波初至自动拾取的精度较高。下文通过增加地震事件数，继续利用这几种方法进行拾取精度对比分析，验证初步结论的准确性。

选取不同隧道超前探测项目中记录到的数据，以1次隧道超前探测产生的所有地震数据合称为1次地震事件。事件1至事件4表示4次超前探测产生的地震数据，其原始数据如图12所示。利用几种算法继续进行初至拾取，将拾取的各个道到时与手动拾取的结果进行差值，并计算所有地震事件拾取结果的均方差值，结果见表3。

(a) 事件1

(b) 事件2

(c) 事件3

(d) 事件4

Table 3 Deviation of four first arrival time pick-up and manual pick-up ms

由表3可知，4种初至拾取精度最高的是STA/LTA法，最低的是最大振幅法，精度由高到低依次为STA/LTA法、能量比法、AIC法、最大振幅法。

3 结论与讨论

在地震数据质量及信噪比一定的条件下，地震波初至时间拾取算法影响初至时间拾取精度，通过对比分析人工拾取法与4种自动拾取算法提取的初至时间，得出如下结论：

1)特征函数可以很大程度上提高初至拾取所依赖的数据的质量，突出初至波到来时的各种震相变化幅度。

2)初至自动拾取方法减小了人工拾取P波到时的不稳定性及个体因素的差异性，自动拾取初至算法拾取的初至到时稳定、耗时小，效率高。通过对地震事件中19道地震数据测试，人工拾取法平均耗时为12 min，自动拾取法平均耗时为5.68 s。

3)在自动拾取算法中，STA/LTA算法拾取的初至时间精度最高。通过对STA/LTA法、AIC法、能量比法、最大振幅法4种地震波初至时间对比分析，STA/LTA法与人工拾取法的结果平均差值为0.46 ms，能量比法为0.89 ms，AIC法为1.03 ms，最大振幅法为4.89 ms。

几种初至自动拾取方法对于干扰大、信噪比低的数据，其精度均会降低。有待针对低信噪比数据，引入抗干扰能力更强的算法获得更加高效、准确的拾取结果。