地震波视出射角求取软件的C++实现

李 稳,吴 晨,窦海岳,吕子强

(1.中国地震局地球物理勘探中心,郑州 450002;

2.山东省地震局,济南 250014)

地震波视出射角求取软件的C++实现

李 稳1,吴 晨2,窦海岳2,吕子强2

(1.中国地震局地球物理勘探中心,郑州 450002;

2.山东省地震局,济南 250014)

在分析数字地震波形记录数据格式和信号处理流程的基础上,选用C++作为编程语言开发了专用的地震波视出射角求取软件。该软件利用了C++语言的优势,具有高性能数值计算能力和高级交互式界面设计。软件中使用的波形资料处理流程对地震信号处理研究有参考价值。具体算法均在处理大量实际资料的过程中得到了严格检验,高效可靠。

地震波视出射角;联合定位;地震波理论与应用;C++

0 引言

通过充分挖掘地震波所携带的丰富信息,增加独立的观测参数以减轻地震定位反问题的多解性是地震科学研究中的重要课题。地震波视出射角参数是一个独立的观测参数。理论计算和实际地震资料处理结果表明,利用地震波视出射角参数与到时数据进行联合定位可有效提高地震定位对震源深度的分辨,从而提高地震定位精度[1]。

本文在分析数字地震波形数据格式和处理流程的基础上,选用C++作为编程语言,在Visual C++6.0和MFC 6.0环境下编写了ProcessDSR可视化交互分析软件。该软件专用于求取数字地震记录中的地震波视出射角参数。

1 地震波视出射角的概念

人们在地球表面观测和接收地震波。地震波入射射线与地面的夹角是真出射角,常用e表示。地面位移矢量与地面的夹角称为视出射角,常用¯e表示[2]。视出射角的计算公式为公式(1)。视出射角与真出射角的关系式为公式(2)。图1是视出射角与真出射角示意图。

图1 视出射角与真出射角示意图

式中,AZ为地动位移的垂直分量,ANS为N S分量,AEW为EW分量。VP为纵波速度,VS为横波速度,其比值称为介质的泊松比。

2 数字地震波形记录数据格式

数字地震波形记录以特定的数据格式存储,常见的有 SEED、SAC、EVT等。甘肃省地震台网提供的波形记录采用EVT格式。

EVT格式包含事件头段结构、台网参数块结构、台站参数结构、通道参数结构、索引块结构和数据块结构几个部分。地震波原始数据样本经STEIM2压缩算法压缩后保存在数据块结构中[3]。

3 数字地震波处理流程

台网记录到的地震波原始信号包含着一些干扰因素,在求视出射角参数之前需对其进行扣除仪器响应、滤波等常规处理。本文参考前人的研究成果[4],设计了规范有效的数字波形资料处理流程,见图2。

图2 求取地震波视出射角参数的流程图

4 软件实现

C++是静态类型语言的典型代表。在高级系统程序设计(如用户界面设计、图像和声音的编排等)和科学计算等领域有明显的优势[5]。地震波视出射角参数求取软件既需实现复杂、高效的数字信号处理算法,又要提供高性能的用户界面设计。经过分析讨论,我们选用C++作为编程语言进行软件开发。在Visual C++6.0和MFC 6.0环境下编写、编译,在win32平台上进行系统测试。

4.1 数据读取

在EVT数据格式中存在由某些字段值决定内存分配的情况,因此在定义类时需用到一级和二级指针。为了安全,必须禁止类对象之间的按值传递[6]。另外,为避免内存泄露,在使用指针前将其初始化,在不再使用时释放其所指向的内存并将其复位为0,同样非常重要。

设计了CDSR类用来存储EVT波形数据。该类具有基本的数据读入、提取功能,并通过声明私有拷贝构造函数的方法禁止按值传递。

ReadDSR()函数的实现利用了 EVT文件将二进制数据封装在文件头、台站参数、波形数据块3个部分的特点。首先建立好嵌套结构;利用fread语句读取文件头字段信息;根据其中提供的台站总数计算台站参数部分所占内存,然后将台站参数部分以数据块的模式读入;再根据台站参数部分提供的采样率、记录长度、数采字长等信息分配内存,读取波形数据块。

将数据读入CDSR类的实例之后,可调用DistillData()函数在地震波初至附近提取4096个采样点(通常为81.92s)的计算数据。该取值既能保证截取的信号具有足够高的频率分辨力,又适宜于进行快速傅氏变换和小波变换。

4.2 扣除仪器响应

为了消除波形记录中各频段灵敏度的影响,必须扣除仪器响应[7-8]。首先根据频率响应参数进行频点灵敏度插值(利用友元函数MyChaZhi()进行计算);然后对提取的波形数据做傅氏变换,将所得复数序列的幅值除以频点灵敏度插值结果;最后进行反傅氏变换,所得复数序列的实部即为地震波数据扣除仪器响应后的结果。

实现傅氏变换时我们选用了时间抽选奇偶分解傅氏变换算法,代码如下:

4.3 滤波

长周期低频信号和高频干扰信号对地震波有效信号有着严重的影响,可设计数字滤波器将它们滤除。

Ⅰ 采用窗函数法设计线性相位的 FIR带通滤波器[9]。

(1)根据实际问题确定技术指标

选取初始的滤波器通带下边缘截止频率为2.5Hz,通带上边缘截止频率为20.0Hz,过渡带宽为1.25Hz,这些参数以后可在交互分析界面中进行调整。要求阻带衰减不小于-60dB。

(2)用因果离散时间系统逼近能满足上述技术指标的滤波器

用有限长的单位冲击响应 h(n)逼近无限长的理想单位冲击响应hd(n)的有效方法是用一个有限长的窗函数序列 w(n)来截取 hd(n)。考虑到幅度要求指标(阻带衰减不小于-60dB),我们选用了Blackman窗作为窗函数,其阻带衰减为-74dB。最终计算出的 h(n)的表达式为:

式中,ωc1=0.075π,ωc2=0.825π,它们是两个理想低通滤波器的截止频率;N=221,是单位冲击响应 h(n)的长度;α=110,是为满足线性相位响应所必需的移位值。

(3)系统实现

得到了以单位冲击响应表示的系统表达式以后,滤波器既可以用硬件来实现,还可以在计算机上通过编程用软件实现。

实现滤波器的算法很多,本软件采用了由Grant R.Griffin 发布的fir—split算法(文件fir—algs—1-0.c[10]),该算法运算速度快,且避免了使用环形缓存。

Ⅱ 小波分析滤波

考虑到地震波信号频率成分随时间变化的特点,进一步采用了小波分析滤波。小波分析滤波可通过在二进小波分析(常用Mallat二进小波算法)中适当地选择尺度因子来实现[11-13]。

4.4 倾斜校正与水平基线校正

进行倾斜校正和水平基线校正的目的是去掉记录中的直流分量。倾斜校正可通过最小二乘法线性回归实现。最小二乘法线性回归的计算公式和实现代码如下:

设回归直线为:y=a+bx用最小二乘法拟合回归直线的解为:

水平基线校正即先求得地震记录序列的平均值,再将该值从整个序列中减去。水平基线校正应在倾斜校正之后进行。其实现代码如下:

4.5 数值积分

数字地震记录是离散的数据集,不适合用牛顿-莱布尼茨公式积分,可用梯形公式(1阶的牛顿-柯特斯公式)等数值积分法积分[14]。其计算公式和实现代码如下:

4.6 二次滤波

一次积分运算在频率域中相当于频谱除以圆频率[4],所以在积分运算中,原始信号的低频成分被放大了,甚至将有用的信号淹没。因此,在积分运算之后还应再次进行滤波处理。

4.7 处理结果屏幕显示与交互分析

在地震波视出射角求取软件中除保证数据处理功能的稳定高效之外,提供一个高性能的交互分析界面也非常重要。

本论文参考“EDSP-IAS”软件界面设置,应用文档/视图结构和非模态对话框实现交互分析功能。在软件设计过程中严格遵循面向对象编程的思想,软件的修改、维护非常方便。该软件提供的具体功能包括:交互选择计算数据、通过下拉组合框设置滤波器参数、通过单选按钮选择显示速度记录还是位移记录、通过复选框选择是否在屏幕上叠加原始记录、通过拖动滑动条设置放大倍数等。计算结果均在屏幕上即时显示。

5 实际应用

利用该软件分析处理了2002年12月玉门地震序列的数字地震波形记录,共求得247条地震波的视出射角参数。应用地震波到时与视出射角联合定位方法得到74个地震事件的定位结果。将联合定位结果与目前工作中常用的双差法定位结果进行了对比研究。

参照走时残差的定义,定义视出射角残差为理论值与观测值的差值,见公式(3)。则平均视出射角残差的计算公式为公式(4),视出射角残差的均方差的计算公式为公式(5)。

表1显示了不同方法定位结果的残差分析结果。可见,增加对视出射角参数利用的联合定位结果,具有更小的平均走时残差和走时残差均方差。这说明,联合定位结果不但整体上更加准确,而且更加稳定。

表1 玉门地震序列地震定位残差分析结果表

图3、图4分别是玉门地震序列近垂直于断层走向的深度剖面图和断裂、三维震中分布图。可见联合定位结果在震源深度方面更加收敛,显示出了明显的倾角,能更有效地反映发震构造的深部信息。据此勾画出的断层性质与前人的地质构造研究成果[15-16]相吻合。

图3 玉门地震序列深度剖面图

图4 玉门地震序列断裂、三维震源位置分布图

通过残差分析和与地质研究成果的对照分析,可以看出增加对视出射角参数的利用后所得的定位结果更加准确、可靠。

6 结论与讨论

地震定位受诸多因素的影响。在现有条件下,通过增加对地震波视出射角(尤其是震中附近台站监测到的地震波视出射角)的利用,获得更准确的地震深度信息,对研究发震断层的性质及其空间展布有实用价值。我国已在全国范围内建立了数字地震观测台网,地震波到时与视出射角联合定位方法具有广泛适用性。

C++语言具有结构规范、便于调试和类型安全的特点,在数值/科学计算、通用应用程序设计等领域具有广泛的应用前景。地震波分析处理软件既需具备高性能数值计算能力又要提供高级交互式界面设计,宜选用C++语言开发实现。

在频率域中求取视出射角参数,以提高计算效率和准确度是可行的研究方向,相关研究已经取得了初步成果[17]。另外,完善该软件使其可以直接处理多种不同格式(如SEED、SAC等)的地震数据,也是一个切实的问题。

[1] 李稳,张元生,何斌.地震波到时与视出射角联合定位方法研究[J].西北地震学报,2009,31(3):207-210.

[2] 傅淑芳,刘宝诚.地震学教程[M].北京:地震出版社,1991:83-85.

[3] 马中华.台站数据共享技术研究[D].甘肃:中国地震局兰州地震研究所,2008.

[4] 陈运泰,吴忠良,王培德,等.数字地震学[M].北京:地震出版社,2000:55-65.

[5] Bjarne Stroustrup.C++语言的设计和演化[M].裘宗燕译.北京:机械工业出版社,2002:150-154.

[6] Bruce Eckel.C++编程思想(第二版)[M].刘宗田,袁兆山,潘秋菱,等译.北京:机械工业出版社,2002:256-257.

[7] 李媛媛,吴东.山西数字遥测地震台网十五勘选子台台址地动噪声分析[J].华北地震科学,2004,22(1):60-62.

[8] 杨晶琼,杨周胜.云南地震数字遥测台网子台噪声分析[J].地震研究,2005,28(1):86-89.

[9] 陈玉东.数字信号处理[M].北京:地质出版社,2005:146-148.

[10] Grant R and Griffin.Digital Signal Processing Software[CP/OL].2008-11-28.www.dspguru.com/sw/lib/fir—algs—1-0.c.

[11] 连海宁,谢礼立.三种变换在强震记录分析方面的研究[J].华北地震科学,2007,25(3):28-33.

[12] 刘泽民.小波分析在数字地震资料中的应用[J].地震地磁观测与研究,2004,1(25):88-91.

[13] 曾宪伟.利用小波包变换识别地震和爆破[D].甘肃:中国地震局兰州地震研究所,2008.

[14] Pallab Ghosh.数值方法(C++描述)[M].北京:清华大学出版社,2008:362-371.

[15] 何文贵,郑文俊,赵广堃,等.2002年12月14日甘肃玉门5.9级地震的发震构造研究[J].地震地质,2004,26(4):688-697.

[16] 陈柏林,刘建生,张永双,等.玉门断裂全新世活动特征及其与玉门地震的关系[J].地质论评,2005,51(2):138-142.

[17] 何斌,张元生,李稳.计算地震初至波视出射角方法[J].西北地震学报,2010,32(1):11-15.

The C++Realization of Seismic Apparent Emergence Angle Calculation Softw are

LI Wen1,WU Chen2,DOU Hai-yue2,LV Zi-qiang2
(1.The Geophysical Exploration Center of CEA,Zhengzhu 450002,China;2.Earthquake Administration of Shandong Province,Jinan 250014,China)

On the basis of analysis on digital seismic wave data format and signal processing,using the C++language,we developed software for seismic apparent emergence angle calculating.Utilized the advantages of C++,the software has high-performance numerical computing ability and advanced interface design.The seismic wave processing design used in the software is valuable for the study of seismic signal processing.The computerized algorithms had stood a severe proof in the process of processing arithmetic real data.They are efficient and reliable.

apparent emergence angle;multiparameter seismic location;theory and applications of seismic waves;C++

P315.69

A

1003-1375(2010)04-0009-06

2010-04-01

国家自然科学基金(40874029)资助;

李稳(1983-),男(汉族),河南开封人,助理工程师,主要从事地震波理论与应用研究.E-mail:liwen—forwork@163.com