(Fujian Earthquake Agency，Fuzhou 350003，China)
DOI: 10.13512/j.hndz.2021.02.06
备注
针对地震台站勘选工作中测试数据处理繁琐的问题，提出采用QT程序开发框架和GMT绘图工具，开发一套勘选测试数据自动处理软件。该软件支持自动处理多种格式的波形数据，自动完成地震噪声处理的所有过程，产出噪声计算结果，绘制PSD、PDF以及RMS等图件。该软件应用于福建地震台网的台站勘选测试数据处理工作，提高了勘选工作效率，取得了良好的应用效果。
Aiming at the problem of tedious process in the seismic site selection work，the paper proposes to develop a set of software for automatically processing seismic noise recorded during measurements for seismic stations selection by using QT program development framework and GMT drawing tools. The software supports automatic processing of waveform data in multiple formats，automatically completes all processes of seismic noise processing，produces noise calculation results，and draws maps of noise power spectral density(PSD)，noise power spectral probability density(PDF) and noise effective amplitude values(RMS).The software is applied to seismic site selection work in Fujian Seismic Network，which improved the efficiency of the selection work and achieved good application results.
引言
在地震台站勘选过程中，需要在拟新建台址上进行连续观测48 h以上的场地噪声测试。测试结束后，按照规范要求对测试数据进行处理，获得台址地震噪声计算结果，用于分析和评估拟建台址地震噪声水平是否满足台站观测要求[1]。
处理勘选测试数据包括数据解析、数据预处理、噪声功率谱密度(PSD，Power Spectral Density)计算、速度有效振幅值(RMS，Root Mean Square)计算、绘制各种图件等多个复杂步骤，人工处理不但需要处理人员掌握复杂的数据处理方法和多种处理工具的使用，而且处理过程极为繁琐。廖诗荣与陈绯雯[2]开发了一套应用地震噪声功率谱概率密度函数(PDF，power spectrum probability density function)方法自动处理勘选测试数据的软件。该软件只能处理SEED格式的波形数据，当处理数据采集器记录的波形数据时，须借助其它软件进行数据格式转换，操作复杂，工作量大。本文吸收该系统的核心算法，在QT开发框架下结合GMT(Generic Mapping Tools)[3-4]绘图工具，开发一套地震台站勘选测试数据自动处理软件，该软件具有人机交互界面，参数配置简便，可对我国常用数据采集器记录格式文件进行批量处理，并自动产出多种结果图件。
1 软件设计
1.1 功能设计地震台站勘选测试数据处理软件主要处理从数据采集器中导出的波形文件。在国内各地震观测台网中，主要配置的数据采集器有北京港震公司的EDAS型和英国GURALP公司的CMG-DM24型。因此软件不仅需要能解析SEED、MiniSEED和EVT[5]等通用的数据交换格式的波形文件，还需要能解析EDAS数据采集器记录的DAT格式和CMG-DM24数据采集器记录的GCF格式的波形文件。
根据地震台站勘选规范要求[6]，处理结束后，需要得到每个通道的RMS总平均值和每小时RMS值的结果文件，还需要自动绘制噪声功率谱概率密度图(PDF图)、噪声加速度功率谱密度图(PSD图)、PSD值在某一频点处随时间变化图(PSD时间曲线图)、PSD值在每个频点随时间变化图(PSD时频图)、RMS值每小时柱状分布图(RMS柱状分布图)和RMS值在每个频点随时间变化图(RMS时频图)等各种图件。
另外，软件还需实现一些在勘选工作中可能用到的功能，包括通过地震仪的仪器响应文件绘制仪器的幅频特性曲线图；通过输入拟建台址的经纬度计算拟建台址距最近海岸线的距离；通过输入拟建台址的经纬度识别拟建台址所在地名并绘制标注拟建台址的点位图。
1.2 结构设计地震台站勘选测试数据处理软件采用QT开发框架，基于C++语言进行开发。在软件内部通过设计自动调用GMT4绘图命令程序实现自动绘图功能，软件系统结构如图1所示。
图1 地震台站勘选数据处理软件系统结构图
Fig.1 Structure chart of data processing software for seismic stations selection在软件界面输入台站勘选测试波形文件后，根据波形格式，选择相应的波形解析模块把压缩数据解析为原始波形记录数据并存入波形缓存池，数据处理模块从波形缓存池中按设置的计算长度截取波形数据进行PSD值计算、PDF计算和RMS计算，最后在GMT绘图模块中自动绘制PDF图、PSD图和RMS图等各类图件。
1.3 界面设计在QtDesigner软件界面设计平台上，通过拖拽的方式实现地震台站勘选测试数据处理软件界面的设计。软件界面如图2所示，界面左半部分为参数设置和处理过程显示区，右半部分处理结果显示区。参数设置包括仪器的记录频带范围、RMS值的输出单位、绘制PSD值随时间变化的频点值和仪器响应文件或灵敏度值。选择勘选测试数据文件夹和波形格式，在输出框中列出所有待处理文件，开始处理后，输出框实时显示处理过程信息。在处理结果显示区中，RMS结果框输出三分向平均RMS值，底部图片显示框可通过按钮查看产出的各类图件。勘选人员通过软件界面显示的处理结果，可以初步判断该台址的噪声环境水平是否符合地震观测要求。
图2 地震台站勘选数据处理软件界面
Fig.2 The interface of data processing software for seismic stations selection2 实现过程
3 应用示例
福建大田均溪测震台因灾重新选址新建，在拟新建台址处安放GURALP公司的CMG-40TDE-40S型地震仪连续测试48 h。采用地震台站勘选测试数据处理软件处理测试波形数据，获得PSD图、5 Hz频点处PSD值随时间变化图、PSD值每个频点处随时间变化图、PDF图、RMS每小时柱状分布图、RMS值每个频点随时间变化图和三个分向RMS总平均值等处理结果及图件。
3.1 噪声加速度功率谱密度分析噪声加速度功率谱是分析台址地震噪声水平的重要指标[5]。本软件自动产出PSD图、特定频点PSD值随时间变化图、PDF图和PSD值每个频点处随时间变化图等四类图片，从而实现对台址噪声加速度功率谱密度进行全面分析。
(1)噪声加速度功率谱密度图(PSD图)如图7所示，两根青色曲线是新全球高噪声模型(NHNM)和新全球低噪声模型(NLNM)，其他三条不同颜色曲线表示三个通道最大概率的PSD值。通过该图，可以直观的看到三分向在各频段内噪声水平，以及与NLNM比较情况。
(2)某特定频点处PSD值随时间变化图如图8所示。在规范中要求5 Hz频点处PSD值随时间变化图，在宽频带或者甚宽频带地震计还要5 s频点和100 s频点变化图，本软件支持产出用户自定义频点随时间变化图。
(3)图9是垂直分向的加速度功率谱概率密度图(PDF图)，图中横坐标为频率，纵坐标是加速度功率谱，颜色块表示概率。在PDF图中，通过观测时间段内加速度功率谱的概率分布情况，可以了解非天然地震事件持续时间占记录时间的大致比例。
(4)垂直分向上，每频点PSD值随时间变化图如图 10所示，横坐标为时间，以小时为单位，纵坐标采用对数坐标表示频率，用颜色块表示该频点该时间点处的PSD值。通过该图，可以更为全面了解噪声在时间和频率上的分布情况。
图7 噪声加速度功率谱
Fig.7 Power spectral density of noise acceleration
图8 5.0 Hz频点处PSD值随时间变化
Fig.8 PSD value changes with time at frequency of 5.0 Hz
图9 加速度功率谱概率密度
Fig.9 Probability density of acceleration power spectrum
图 10 每频点PSD值随时间变化
Fig.10 The PSD value of each frequency point changes with time3.2 噪声速度振幅有效值分析图 11是RMS值每小时的柱状分布图，一个小时的RMS值取该小时内所有中心频率处的最大概率RMS值的算术平均值。图 12是垂直分向RMS值在每频点随时间变化图，横坐标为时间，以小时为单位，纵坐标频率为对数坐标，颜色块表示RMS值大小。在RMS值柱状分布图中，可以很直观地看到RMS值在每个小时的具体值，以及整个观测期间RMS是否平稳，通过RMS值在每频点随时间变化图，可以进一步了解到噪声在频率上的分布情况，为分析噪声的特征提供了依据。另外，软件自动生成每小时RMS值的文本文件，并把所有小时的RMS值的平均值显示在软件界面中。
图 11 RMS值每小时柱状分布图
Fig.11 Hourly histogram of RMS value
图 12 RMS值每频点随时间变化图
Fig.12 The RMS value of each frequency point changes with time4 结语
地震台站勘选测试数据处理软件按照规定的数据处理方法实现自动处理勘选测试数据，并能自动产出丰富的结果图件，解决了地震台站勘选人员处理勘选测试数据流程繁琐的问题，提高了台站勘选人员的工作效率，也使处理过程更为规范化。该套软件应用于福建地震台网的勘选测试数据的处理，自动产出了该台址地噪声PSD图、PDF图和RMS图等图件和数据，为评估拟建台址的地震噪声水平提供了充足且可靠的依据，取得了良好的应用效果。
地震台站勘选测试数据处理软件还可以应用于其它地震噪声分析工作。例如，通过处理地震台站记录的波形数据，绘制台站的PDF图，可以分析地震台站观测数据质量，也可以判断仪器工作状态是否正常。
- [1] 中国地震局. GB/T 19531_1-2004地震台站观测环境技术要求第一部分：测震[S]. 北京：地震出版社，2004.
- [2] 廖诗荣，陈绯雯. 应用概率密度函数方法自动处理地震台站勘选测试数据[J]. 华南地震，2008，28(4)：82-92.
- [3] Wessel，W.H.FSmith. New improved version of Generic Mapping Tools released[J]. EOS Trans. AGU，1998，79(47)：579.
- [4] Wessel，W.H.FSmith. Gridding with continuous curvature splines in tension[J]. Geophysics，1990，55(3)，293-305.
- [5] 中国地震局监测预报司. 数字地震观测技术[M]. 北京：地震出版社，2003.
- [6] 中国地震局. JSGC-01中国数字测震台网技术规程[S]. 北京：地震出版社，2005.
- [7] 北京港震机电技术有限公司. 地震数据采集器EDAS-24IP用户指南[R]. 北京：北京港震机电技术有限公司，2008.
- [8] 北京港震机电技术有限公司编. 地震数据采集器EDAS-24GN用户指南[R]. 北京：北京港震机电技术有限公司，2016.
- [9] 杨周胜，杨晶琼，姚远，等. 地震交换标准数据的压缩技术和数据解压[J]. 地震研究，2019，42(1)：144-149.
- [10] 中国地震局. DB/T2-2003地震波形数据交换格式[S]. 北京：地震出版社，2003.
- [11] Peterson Jon. Observations and Modeling of Seismic Background Noise[R]. Virginia：U.S. Geological Survey，1993.
- [12] McNamara D.E. ，Boaz R.I. Seismic Noise Analysis System Using Power Spectral Density Probability Density Functions：A Stand-Alone Software Package[R]. Virginia：U.S. Geological Survey，2005.
2.1 波形解析波形解析模块由SEED、MiniSEED、EVT和DAT等格式的波形文件解析子模块组成。由于台站勘选测试数据一般都是直接从数据采集器中获取，本文对国内最为通用的CMG-DM24型和EDAS型数据采集器记录的波形文件在软件中的解析过程进行阐述。其中，港震EDAS数据采集器有两种型号，EDAS-24GN型和EDAS-24IP型，两种型号记录的数据文件压缩方式不同。
2.1.1 EDAS-IP型EDAS数据采集器记录的波形文件需要通过读取文件头段的前16个字符判断数据采集器的型号，“SEISMIC EVENT”或 “SEISMIC SERIES”为EDAS-24IP型，“GZDAS”则为EDAS-24GN型[7-8]。
EDAS-IP连续波形的文件由文件头块、台站块、索引块和数据块组成。文件结构如图3所示[7]。
图3 EDAS-IP型数据采集器记录文件结构
Fig.3 Record file structure of EDAS-IP data acquisition device文件头块包括事件头段、台网参数、起始时间和记录长度。在读取文件头块时主要获取事件头段中判断数据采集器型号的事件标志、台网参数中的台网代码和台站个数。在读取台站块数据时，根据台站个数循环读取每个台站和通道的信息并保存。读完文件头块和台站块后，跳过2048字节大小的索引块，读取动态数据块，包括绝对秒(即数据块的时间长度)、压缩数据块长度等参数和压缩数据。动态数据块的具体读取流程如图4所示。
图4 读动态数据块读流程图
Fig.4 Flow chart of dynamic data block reading根据块长度读取一个压缩数据块，采用文献[9-10]中的解压方法把压缩的数据解压为原始数据。数据块是按台站和通道顺序存储的，一个通道的数据大小为采样率乘绝对秒数。解压数据存入缓存时，按通道顺序依次压入由台站序号、通道序号、相对时间秒和1秒采样点数组成的四维数组的缓存池。
2.1.2 EDAS-GN型EDAS-GN型记录的波形文件总体结构与EDAS-IP类似，由文件头块、台站块、索引块和数据块组成，最大的区别是数据块，它的结构如图5所示[4]。
图5 EDAS-GN型数据采集器记录文件数据块结构
Fig.5 Data block structure of record file of EDAS-GN data collector数据区的压缩方法是，对各通道的N秒数据做前后差分diff(n)=x(n)-x(n-1)，以10个差分值diff为一组统计有效字长，然后根据字长在压缩编码表中选择使用的编码，按照选用的编码表中的字长对差分值diff 值重新编码输出。图6是一个通道的压缩数据CMPDATAFRM的结构[8]。
图6 CMPDATAFRM数据结构
Fig.6 Data structure of CMPDATAFRM解压过程为：首先计算差分组组数N=(SP*LEN)/10，SP为采样率，LEN为当前数据块长度(s)；然后读取一级索引和各差分组的二级索引，根据一级索引和差分组的二级索引，查文献[8]中的编码表，获得各差分组的编码长度；最后根据编码长度读取差分值diff(n)，并计算当前原始值X(n)=diff(n)+X(n-1)。
2.1.3 CMG-DM24型CMG-DM24型数据采集器默认的记录格式为GCF格式。在该仪器配套的采集软件SCREAM中带有GCF转MiniSeed格式的子程序gcf2msd。把该子程序移植到地震台站勘选数据处理软件中，当解析GCF格式波形文件时，自动调用gcf2msd转换子程序，把GCF格式转换为MiniSeed格式，再采用MiniSeed解析模块进行解析。
MiniSeed是SEED数据格式的纯数据块，它由若干个长度为4096字节的数据记录组成，每个数据记录包含控制头段区和数据区块组成。在控制头段中，前7个字节是由6个整数和1个字符“D”组成的数据记录标志，通过匹配这个标志循环读取数据记录进行解压，在控制头段信息中，读取原始数据的压缩格式，按文献[9]、[10]中的解压方法把数据区的数据解压为原始数据。
2.2 数据处理地震台站勘选波形数据处理的一般过程是：首先把波形截为若干个记录段，然后每个记录段按Peterson[11](1993)的方法计算PSD值，最后采用McNamara等人[12](2005)提出的概率密度函数方法计算PDF值和RMS值。本软件设计能处理最长24 h的数据文件，并将数据段进行分段处理，截取长度默认为720 s，支持用户自行设置，通常长度应大于最长记录周期期的6倍[2]。为使每个记录段的处理结果更为平滑，在截取记录段时以50%的重叠率进行重叠。噪声功率谱密度(PSD)、噪声功率谱概率密度(PDF)和噪声速度有效振幅(RMS)的计算过程如下。
2.2.1 噪声功率谱密度计算一个记录段波形数据的PSD值计算流程如下[2,12]：
(1)数据预处理。为了消除对长周期(即大于仪器记录的最长周期)功率谱估计的影响，对原始波形数据去均值并采用平均斜率法消除长周期线性趋势，处理公式为：
x(t)=u(t)-umean-Tlp(1)
其中x(t)为处理后的数据，u(t)为原始数据，umean为原始数据的算术平均值，Tlp为长周期线性趋势，其处理公式为：
(2)式(2)中，Tr为原始数据u(t)的总长度。
另外，为了减少有限长度数据序列进行FFT变换时造成的频率渗漏，将锥度为10%的正弦窗函数应用于记录段数据序列，使数据段两端平滑地衰减至零。最后PSD 值乘以1.142 857以补偿应用窗函数所造成的功率值损失。
(2)速度功率谱密度值计算。采用快速傅立叶变换把原始数据变换为以频率为自变量的函数Y(f)，然后根据公式(1)计算每个频点的功率谱密度Pv(f)。
(3)式(3)中，N为采样点个数；Δt为采样时间间隔。
(3)扣除仪器响应。为了得到地动噪声的物理量值，即速度值，功率谱Pv(f)需要进行仪器响应校正。
(4)式(4)中， H(f)为仪器的传递函数。
(4)加速度功率谱密度的计算。为了与全球标准低噪声模型(NLNM)和全球标准高噪声模型(NHNM)进行比较，采用公式(5)转换为加速度功率谱密度PSDa。
PSDa(f)=4π2f2PSDv (5)
(5)平滑处理。为了得到 PSD值在频率对数坐标中呈等间隔采样，本文采用1/3 倍频程积分对每条记录的功率谱密度做平滑处理。
(6)式(6)中，fl为低频拐角频率，fh为高频拐角频率，n为介于二者之间频率f的个数，中心频率fc以1/9倍频程为增加步长。
2.2.2 功率谱概率密度计算所有记录段PSD值计算结束后，在-200 dB到-50 dB范围的功率以1 dB为间隔划分成连续的功率窗，统计每个频点的PSD值落在对应功率窗内的记录段个数，然后根据概率密度函数公式(7)计算该频点处各个功率窗的概率[2]。
Ppsd(fc)=Npfc/Nfc (7)
式(7)中，Npfc为在中心频率fc处落在某一功率窗的记录段个数，Nfc为总记录段个数，Ppsd(fc)为fc频点处某一功率窗的概率。
2.2.3 噪声速度有效振幅计算提取中心频点在1~20 Hz频带范围内所有PSD值，按公式(6)计算各中心频点RMS值[1-2]。
(8)式(8)中，PSDv为速度功率谱密度，fc为中心频率，RBW为相对带宽。
根据每个中心频点fc的概率密度函数公式(9)计算该中心频点在各个RMS值从-210至-70 dB区间以1 dB为间隔的区间的概率[2]。
Prms(fc)=NRfc/Nfc (9)
统计每个中心频率的概率后，取各中心频率最大概率处的RMS值求算术平均值，该值作为该通道数据的最后计算RMS值。计算公式为公式(10)[2]。
(10)2.3 GMT绘图在GMT绘图模块中，主要包括绘制PSD图、PSD时间曲线图、PSD时频图、PDF图、RMS柱状分布图和RMS时频图等绘图子模块。每个绘图子模块自动生成一个BAT批处理文件，并在批处理文件中写入绘制相应图件的所有GMT命令。通过QT中的QProcess类调用执行批处理文件实现自动绘图。
本文中采用的GMT命令模块为4.5版本[3-4]。在绘制PSD时间曲线图和RMS柱状分布图中，横轴为时间单位，绘底图采用psbasemap命令，其中设置-B选项为“-Bs4Hg4H”，使坐标刻度标记和网格为4 h。绘制数值曲线和柱状图均采用psxy命令。绘制柱状图符号的参数为“-Sb0.01ub-180”，其中，-Sb为绘制柱状参数，0.001u表示柱状宽度为0.001个坐标轴单位，b-180表示柱状由-180延长到y值。
在绘制PDF图子模块中，由于PDF图是由三维数据绘制的二维图片，PDF数据文件的第1列为频率，第2列为PSD值，第3列为概率，因此，绘制前需要把PDF数据进行网格化处理，首先gawk命令把第1列频率进行对数计算，然后用blockmedian 命令进行插值，使数据成为规则的等间隔数据，最后surface命令网格化生成网格文件。绘图时，先用psbasemap命令绘制出横轴对数坐标刻度标记，然后再采用grdview命令绘制概率彩色图。
在绘制PSD时频图和RMS时频图中，数据文件均为三维数据，第1列数据为时间，第2列为频率，第3列为PSD值或RMS值。在数据文件中频率取数，然后进行插值网格化生成网格文件。采用psbasemap命令绘制横轴为时间坐标刻度标记，纵轴频率对数坐标刻度标记，然后用grdview命令绘制PSD值或RMS值彩色图。
