基于分形盒维算法的天然地震波初至拾取方法研究

张 琪， 王有学， 熊 彬， 蒋婵君， 王心宇， 曾 成

(桂林理工大学 地球科学学院，桂林 541006)

0 引言

在天然地震数据的处理中，拾取初至信号是地震波走时成像、研究地球内部结构的基础。实际工作中，根据地震波的振幅和波形特征，利用人工拾取初至信号。然而，天然地震波含有各种干扰，人工提取存在的主观偏见给提取初至信号带来了极大的不唯一性。因此，以计算机为基础的自动拾取方法毋容置疑地可以给我们提供更为准确的地震波旅行时间信息。

对于自动拾取地震波初至到时，地震科学家们提出了许多方法，诸如相关法、神经网络法、分形法等[1]。Gelchinsky等[2]提出了一种基于地震信号的特定相关属性来确定初至的方法;王继等[3]用多道互相关和最小二乘方法确定远震震相初至时间;潘树林等[4]利用优化相关法拾取地震事件的初至信号，这些方法均基于地震道的整体特征，对噪声有较好的抑制作用，但其拾取结果仍然受到续至波的影响较大。Murat等[5]提出了在低信噪比的信号中反向追踪神经网络技术能在噪声干扰背景下识别初至，但神经网络拾取初至学习速度慢，实现复杂，需要综合考虑网格结构、训练方法等问题[6];美籍数学家曼德布罗特[7]提出分形思想并于1975年创建分形理论，分形盒维数(Fractal Box)作为分形维数的一种定义，它的大小及变化可以反映信号的不规则度和复杂度。这种方法容易理解，易于数学计算且直观形象，所以分形理论在解决自动提取地震波初至到时的问题上发挥了长足的作用[8-10]。

对地震信号的抽样序列来说，将其作无量纲图形化处理可得到一条自仿射曲线，具有不规则性和近似的自相似性,可以用分数维来描述其不规则程度[11]。研究发现，在地震波初至到达之前，地震监测系统监测到的地震信号时间序列为随机噪声的响应，地震波初至到达之后，地震信号时间序列为随机噪声与地震波动的复合响应[12]。由于地震信号的分段特征，在地震波初至到达前后地震波时间序列的分形特征会发生突变，根据其特征的变化可以定量地确定初至波的到达时间[13]。

笔者基于分形盒维算法(FBM- Fractal Box Method)提出了一种天然地震波初至拾取的方法，该方法能够较好地解决自动拾取地震波初至的问题。

1 分形盒维算法的基本原理

维的定义基于δ规模下的测量法的思想。设X是Rn中任一非空有界子集，对于一段离散的数字信号x(i) (i=1，2，…，N)，x(n)∈X，将一系列尺度为δ(i) (i=1，2，…，M)的方网格对X进行覆盖(图1)，得到有效覆盖网格数Nδ(i)(X) (i=1，2，…，M)，则X的上、下盒维定义[14]分别为

(1)

(2)

如果X的上、下盒维相等，则X盒维数为

(3)

离散的数字信号无法求得式(3)中的极限值，实际计算中，通过最小二乘法得到lnNδ(i)(X)～ lnδ(i)的拟合直线，其斜率就是的分形盒维数p。

图1 数字信号分形维数的计算原理示意图Fig.1 Calculation principle of fractal dimension on digital signal

对于地震波记录来说，它是一个离散序列x(i) (i=1，2，…，N)，我们在一个相对固定的时窗Wk(k=1，2，…，N)(长度为l)内，将时窗M等分，一般M取2的幂次方。以信号零值为零线，以相同的等分间隔画平行于零线的平行线，使得最外侧的平行线包括了整个时窗内的地震信号，此时，在时窗内形成一个边长为l/(M-1)的正方形阵列。然后对时窗内的地震信号以Δt=l/(M-1)进行重采样，并对时窗内重新采样的地震信号计算与一系列尺度为δ(i)=l/(M-1)*(i-1)对应的Nδ(i)，最后求得该时窗的分形盒维数pk。按照地震记录的采样间隔移动时窗，重复上述步骤，得到每一个时窗所对应的分形盒维数，直到覆盖整个信号序列。

2 时窗大小对分形盒维的影响

通过对每个时间窗口下进行分形盒维数的计算，可得到地震波信号的分形盒维数轨迹。但是，时窗长度过小，计算量增大；时窗长度过大，降低了信号识别能力。针对时窗长度选取及它对分形维的影响,又分别选取时窗长度l=0.3 s、0.6 s、1 s、2 s这四个取值进行测试，其结果如图2所示。

由图2可以看出，当时窗长度取0.6 s时，分形盒维数曲线在有效信号附近下降趋势明显(图2(c))，计算速度合适。

3 分形盒维在地震波初至拾取中应用

根据分形盒维方法，我们对广西地区宽频带地震流动台网记录的地震事件 (表1、图3(b)、图4(b)、图5(b)、图6(b))进行盒维分析。地震波采样间隔10 ms，时窗长度取0.6 s。为了求取分形盒维数的变化率，对分形盒维数曲线进行了平滑[15](图3(c)、图4(c)、图5(c)、图6(c))，其计算结果如图3～图6(b)～图6(d)所示。

由图3可以看出，在10 s以前分形盒维数p在1.5左右，而在10 s以后，分形盒维数p急剧降至1.2以下，由分形盒维数p的变化率(图3(d))自动拾取初至时间为t=9.51 s；在图4中，在3 s以前分形盒维数p在1.5左右，而在3 s以后，分形盒维数p明显地降至1.4左右，由分形盒维数p的变化率(图4(d))自动拾取初至时间为t=3.08 s。在图5和图6中，分形盒维数p具有明显变化，由分形盒维数p的变化率(图5(d)与图6(d))自动拾取初至时间分别为t=7.70 s、t=18.91 s。

地震事件时间/UTC经度/°E纬度/°N震源深度/km震级震中距/°马里亚纳群岛地震2016/7/29 21:18:26145.5418.54207.627.733.47广西梧州地震2016/7/31 9:18:13111.4724.1322.554.92.05俾斯麦群岛地震2016/8/31 3:11:35152.79-3.69476.006.849.77日本本州地震2016/11/21 20:59:50141.3537.3914.696.929.57

图3 马里亚纳群岛地震记录及其分形计算结果Fig.3 Seismic data from earthquake in Mariana islands and computed results by FBM(a)地震信号；(b)地震信号的分形盒维数；(c)平滑后的分形盒维数；(d)分形盒维数的变化率

图4 广西梧州地震记录及其分形计算结果Fig.4 Record data from earthquake in Wuzhou and computed results by FBM(a)地震信号；(b)地震信号的分形盒维数；(c)平滑后的分形盒维数；(d)分形盒维数的变化率

图5 俾斯麦群岛地震记录及其分形计算结果Fig.5 Seismic data from earthquake in Bismarck islands and computed results by FBM(a)地震信号；(b)地震信号的分形盒维数；(c)平滑后的分形盒维数；(d)分形盒维数的变化率

图6 日本本州地震记录及其分形计算结果Fig.6 Record data from earthquake in Japan's Honshu and computed results by FBM(a)地震信号；(b)地震信号的分形盒维数；(c)平滑后的分形盒维数；(d)分形盒维数的变化率

4 结论

根据分形理论，对基于分形理论的地震波初至拾取的方法进行了研究，主要结论如下：

1)分形维理论可以根据地震波初至到达前后地震波时间序列的分形特征突变，定量地确定初至波的到达时间，能够较好地解决自动拾取地震波初至的问题。

2)基于分形理论的地震波初至时间的拾取方法易于数学计算，其拾取精度可以达到0.01 s；拾取有唯一点，解决了手动拾取带来的不唯一性。对于采样间隔10.0 ms，记录长度在30 s之内的地震信号，耗时为3.42 s，以后可以通过改进算法更好的缩短拾取耗时。

3)试验结果表明，当窗口大小为0.6 s时，计算速度合适，可以保证地震波波至时间的拾取精度。