雷 晨,赵龙梅,朱叶琳,王 乔

（辽宁省地震局，辽宁 沈阳 110034）

0 引言

进行基于地震面波的成像，首先需要从观测台站记录到的大量观测数据中，可靠地提取出对于研究有价值的基础资料，最重要、最繁琐也最基础的工作就是面波相速度频散曲线的测量。20世纪50年代，傅里叶变换技术应用于面波频散测量，60年代初，数值滤波技术和时间变量滤波技术在面波频散测量中得到应用，此后逐渐发展建立了面波频时分析的基础。Landisman等[1]提出，可以利用双台数据互相关提取双台间的面波群速度、相速度。时至今日，用双台波形互相关法来提取面波的频散曲线已经是常用方法之一，面波频散测量工作可以借助多种较为成熟的软件开展，姚华建[2]提出了一种基于图像分析的双台面波相速度频散曲线提取方法，实现了快速、直观、准确地频散测量。近年来，基于辽宁及周边地区固定地震台站记录的资料，开展了一系列地震面波层析成像相关研究工作，在基础数据处理过程中，使用了姚华建编写的TSAnalysis（双台相速度频散曲线快速提取软件）进行资料处理，得到了较好的应用效果。

1 双台互相关方法原理

按照提取一条频散曲线所用到的资料来划分，常用的频散测量方法有单台法、双台法和双事件法等。所谓双台法，顾名思义，是基于成对台站观测数据开展相关研究的方法，且两个台站几乎位于从事件到台站的同一大圆路径上；而互相关，是两个数据之间关于时间的一个函数，反映两个数据的相匹配程度。传统的单台法在测量中需要用到关于震源的大量参数，比如震源机制、发震时刻、震中位置、初始相位等等，而测定这些参数的同时会引入大量误差项，拉低频散的测量精度。相比之下，使用双台互相关方法处理观测数据，就不用考虑震源的这些参数，可以在保留频散信息的同时去掉震源的影响，从而可以提高很多测量精度[3-6]。双台互相关法的原理如下。

双台数据互相关谱函数为：

对于n阶振型，上式(3)则可以重写为：

由式（3）、（4）可知，双台数据经过互相关可以消除震源的相位，只保留下频散信息。式（4）中第一项可以反映出面波在双台之间介质传播的各阶振型的频散信息，而第二项为交叉振型项。如果假设在统计意义下，不同振型相互独立，不同阶振型互相关后叠加会相互抵消，那么式中的第二项可以忽略不计；走时信息属于已知量，可以消去，式（4）可以变为：

当然，双台法也有一些应用限制，首先面波路径要为大圆路径，即所选双台与所选地震震中几乎位于同一大圆路径上，并且双台在震中的同侧；其次所选双台之间的距离要远小于震中距，即要使用远震事件面波观测资料；最后，所研究区域台站分布不能太稀疏，还要筛选数量足够且方位角分布较好的地震事件。

2 基于图像分析的双台相速度频散曲线快速提取软件

TSAnalysis是使用Matlab编写的一款图形界面软件[2,8]，界面如图1所示。在辽宁地区面波成像研究中使用到该软件的主要功能有：地震资料的预处理—包括地震波形读取，数据降采样，去除零漂、仪器响应，振幅归一化；窄带滤波器生成及滤波；基于多重滤波法计算双台群速度频散；基于图像分析提取双台相速度频散曲线等。

图1 TSAnalysis软件图形界面Fig.1 TSAnalysis graphic interface

在地震资料的预处理过程中，通过程序参数设置控制每一步计算，首先对大量采样率为100Hz的地震事件数据降采样，重采样频率为1Hz；然后对数据资料进行零漂校正、去仪器响应、波形归一化等一系列处理；在窄带滤波器的生成和滤波环节，选取了加Kaiser窗的有限冲击响应滤波器，为获得最好的窄带滤波效果，即尽可能抑制高阶面波和其他波的影响，很好的分离出基阶面波，Kaiser窗的特征参数β值取为9。

TSAnalysis软件的核心部分为计算能够直接反映相速度和周期之间关系的振幅矩阵并用图像的方法显示，从而基于图像清晰直观的确定实际频散曲线。为了更好的控制整个过程的正确性和精度，在使用多重滤波法处理双台记录得到每个台站群速度到时后，只选择到时连续清晰的双台资料进行窄带滤波；在采用二次样条插值方法计算双台互相关振幅矩阵时，只选取互相关最明显且相速度值符合天然地震频散特征范围的区域进行计算和插值变换。Yao et al[8,9]为提取高质量的相速度频散曲线采用了多种技巧，使这种基于图像分析技术提取瑞雷面波相速度频散的方法在保证正确性和精度的同时非常直观易于操作。

3 资料处理实例

选取2013年11月25日5:56:53在151.2oE，45.57oN发生的矩震级6级地震事件作为计算实例，地震震源深度35km。

首先加载全部台站记录的地震波形文件和所有设备的仪器响应文件，设置重采样频率为1Hz。设置加Kaiser窗的窄带滤波器的测试中心周期为40s，Kaiser窗的特征参数β值取为9，滤波器带宽设置为1s，于是生成的窄带滤波器如图2所示。

图2 窄带滤波器实例(a.时域；b.频域)Fig.2 Narrowband filter example(a.Time domain, b.Frequency domain)

在程序筛选符合位于同一大圆路径上的两个台站时，限定两个最大偏移角α和β分别为3o和5o（α表示震中到两个台站的方位角差，β表示震中到较近台站和两台站连线之间的方位角差）。此实例中所选取地震记录是位于双台大圆路径上的DHT(敦化台)和ZUH(遵化台)两个台站的波形记录，经过程序降采样和去仪器响应并窄带滤波后提取的波形如图3所示，DHT台站距离震中较近，可见图中上部为DHT台站记录波形，下部为ZUH台站记录波形。

图3 经程序初步处理的大圆路径上双台波形记录Fig.3 Waveform of dual stations on great circle path

图4 群到时图像Fig.4 Image of group arrival time

根据此前试算辽宁地区数据时的经验，利用现有数据提取周期15s至120s的相速度频散是可行的，所以使用软件提取频散过程中对周期范围略放大，起始周期10s，终止周期130s，周期间隔设置为1s。首先由提取的双台记录得到双台群到时图像，如图4所示，颜色标识归一化之后每个周期波包能量情况，红色代表高能量，蓝色相反，绿色曲线即为程序追踪的群到时曲线。图中可见红色区域比较清晰且连续，说明这一组数据的质量合格，可以继续进行下一步计算。加移动窗过滤可以进一步滤除高阶振型和噪声的影响，于是得到更加清晰的群速度到时图像，如图5所示。

图5 加移动窗处理后的群到时图像Fig.5 Image of group arrival time by add moving time window

至此，可以通过程序计算得到互相关振幅矩阵并直接用图像法显示，如图6所示，可以很清晰的看到互相关最为明显的区域(即黄颜色最深的区域)，并操作程序提取相速度频散曲线。如图中黑色实线为程序自动搜索得到的频散曲线，红色点为依照周期范围设置所提取的满足信噪比条件的各周期频散点，这些数据可进入下一步分析。

图6 基于图像分析的双台相速度频散曲线提取实例Fig.6 Example for the determination of inter-station Rayleigh wave phase velocity dispersion curves based on the image analysis technique

4 资料及处理结果

依照本次研究区域范围，向辽宁省地震局监测中心和国家测震台网数据备份中心[10-11]申请了辽宁省内及周边邻省共77个固定台站宽频带地震计记录的2012年1月至2014年6月间，矩震级满足大于5.5级且小于7.5级，震源深度小于100km，震中距满足大于10o小于100o的全部700多个地震事件垂直向波形数据。经过计算和挑选，最终从数百个地震事件中筛选出瑞利面波发育较好，具有比较完备的震中距和方位角覆盖，且频散特征明显的105个地震事件用于本次研究。

图7 双台路径分布图Fig.7 Distribution of paths

使用上述基于图像分析的双台相速度频散曲线提取方法，借助TSAnalysis软件，共提取超过3000条周期范围在15s到120s内的相速度频散曲线(包括重复路径)。经过仔细筛选和计算对比，仍保留了超过2600条的高质量重复路径相速度频散曲线，可见这种方法的质量控制效果还是较为高效及可信的。通过对重复路径频散资料取算数平均数，最后获得超过1100条独立双台路径相速度频散资料，对于整个辽宁地区及周边的覆盖情况如图7所示，图中亮黄色标识辽宁省边界。图8展示了从15s到120s之间不同中心周期的独立路径分布情况，整体覆盖情况均满足下一步工作开展的需求。

图8 各周期独立路径数Fig.8 Number of paths at each interest period

5 结论

本文利用辽宁及邻省固定地震台观测资料，使用TSAnalysis软件，应用基于图像分析的双台相速度频散曲线提取方法处理了105个地震事件的波形记录，成功提取了1100条独立双台路径相速度频散资料，对于整个辽宁地区及周边区域的覆盖情况良好，高质量的相速度频散资料为后续开展的地震面波层析成像工作奠定了良好基础。相速度频散的识别和分析是很多项地球物理研究的基础，基于图像分析的双台相速度频散曲线提取方法除了本文展示的天然地震面波数据处理，还可以应用于噪声成像及其他研究的基础数据处理。相速度频散的有效提取对于更好的研究浅部构造及深部地球动力学过程有着积极的意义。