巫立華，戴麗金，張麗娜，張寶劍，黃鶴凌，廖詩榮

（福建省地震局，福州 350003）

0 引言

在地震臺站勘選過程中，需要在擬新建臺址上進行連續(xù)觀測48 h以上的場地噪聲測試。測試結(jié)束后，按照規(guī)范要求對測試數(shù)據(jù)進行處理，獲得臺址地震噪聲計算結(jié)果，用于分析和評估擬建臺址地震噪聲水平是否滿足臺站觀測要求[1]。

處理勘選測試數(shù)據(jù)包括數(shù)據(jù)解析、數(shù)據(jù)預(yù)處理、噪聲功率譜密度（PSD，Power Spectral Density）計算、速度有效振幅值（RMS，Root Mean Square）計算、繪制各種圖件等多個復(fù)雜步驟，人工處理不但需要處理人員掌握復(fù)雜的數(shù)據(jù)處理方法和多種處理工具的使用，而且處理過程極為繁瑣。廖詩榮與陳緋雯[2]開發(fā)了一套應(yīng)用地震噪聲功率譜概率密度函數(shù)（PDF，power spectrum probability density function）方法自動處理勘選測試數(shù)據(jù)的軟件。該軟件只能處理SEED格式的波形數(shù)據(jù)，當(dāng)處理數(shù)據(jù)采集器記錄的波形數(shù)據(jù)時，須借助其它軟件進行數(shù)據(jù)格式轉(zhuǎn)換，操作復(fù)雜，工作量大。本文吸收該系統(tǒng)的核心算法，在QT開發(fā)框架下結(jié)合GMT（Generic Mapping Tools）[3-4]繪圖工具，開發(fā)一套地震臺站勘選測試數(shù)據(jù)自動處理軟件，該軟件具有人機交互界面，參數(shù)配置簡便，可對我國常用數(shù)據(jù)采集器記錄格式文件進行批量處理，并自動產(chǎn)出多種結(jié)果圖件。

1 軟件設(shè)計

1.1 功能設(shè)計

地震臺站勘選測試數(shù)據(jù)處理軟件主要處理從數(shù)據(jù)采集器中導(dǎo)出的波形文件。在國內(nèi)各地震觀測臺網(wǎng)中，主要配置的數(shù)據(jù)采集器有北京港震公司的EDAS型和英國GURALP公司的CMG-DM24型。因此軟件不僅需要能解析SEED、MiniSEED和EVT[5]等通用的數(shù)據(jù)交換格式的波形文件，還需要能解析EDAS數(shù)據(jù)采集器記錄的DAT格式和CMG-DM24數(shù)據(jù)采集器記錄的GCF格式的波形文件。

根據(jù)地震臺站勘選規(guī)范要求[6]，處理結(jié)束后，需要得到每個通道的RMS總平均值和每小時RMS值的結(jié)果文件，還需要自動繪制噪聲功率譜概率密度圖（PDF圖）、噪聲加速度功率譜密度圖（PSD圖）、PSD值在某一頻點處隨時間變化圖（PSD時間曲線圖）、PSD值在每個頻點隨時間變化圖（PSD時頻圖）、RMS值每小時柱狀分布圖（RMS柱狀分布圖）和RMS值在每個頻點隨時間變化圖（RMS時頻圖）等各種圖件。

另外，軟件還需實現(xiàn)一些在勘選工作中可能用到的功能，包括通過地震儀的儀器響應(yīng)文件繪制儀器的幅頻特性曲線圖；通過輸入擬建臺址的經(jīng)緯度計算擬建臺址距最近海岸線的距離；通過輸入擬建臺址的經(jīng)緯度識別擬建臺址所在地名并繪制標(biāo)注擬建臺址的點位圖。

1.2 結(jié)構(gòu)設(shè)計

地震臺站勘選測試數(shù)據(jù)處理軟件采用QT開發(fā)框架，基于C++語言進行開發(fā)。在軟件內(nèi)部通過設(shè)計自動調(diào)用GMT4繪圖命令程序?qū)崿F(xiàn)自動繪圖功能，軟件系統(tǒng)結(jié)構(gòu)如圖1所示。

圖1 地震臺站勘選數(shù)據(jù)處理軟件系統(tǒng)結(jié)構(gòu)圖Fig.1 Structure chart of data processing software for seismic stations selection

在軟件界面輸入臺站勘選測試波形文件后，根據(jù)波形格式，選擇相應(yīng)的波形解析模塊把壓縮數(shù)據(jù)解析為原始波形記錄數(shù)據(jù)并存入波形緩存池，數(shù)據(jù)處理模塊從波形緩存池中按設(shè)置的計算長度截取波形數(shù)據(jù)進行PSD值計算、PDF計算和RMS計算，最后在GMT繪圖模塊中自動繪制PDF圖、PSD圖和RMS圖等各類圖件。

1.3 界面設(shè)計

在QtDesigner軟件界面設(shè)計平臺上，通過拖拽的方式實現(xiàn)地震臺站勘選測試數(shù)據(jù)處理軟件界面的設(shè)計。軟件界面如圖2所示，界面左半部分為參數(shù)設(shè)置和處理過程顯示區(qū)，右半部分處理結(jié)果顯示區(qū)。參數(shù)設(shè)置包括儀器的記錄頻帶范圍、RMS值的輸出單位、繪制PSD值隨時間變化的頻點值和儀器響應(yīng)文件或靈敏度值。選擇勘選測試數(shù)據(jù)文件夾和波形格式，在輸出框中列出所有待處理文件，開始處理后，輸出框?qū)崟r顯示處理過程信息。在處理結(jié)果顯示區(qū)中，RMS結(jié)果框輸出三分向平均RMS值，底部圖片顯示框可通過按鈕查看產(chǎn)出的各類圖件?？边x人員通過軟件界面顯示的處理結(jié)果，可以初步判斷該臺址的噪聲環(huán)境水平是否符合地震觀測要求。

圖2 地震臺站勘選數(shù)據(jù)處理軟件界面Fig.2 The interface of data processing software for seismic stations selection

2 實現(xiàn)過程

2.1 波形解析

波形解析模塊由SEED、MiniSEED、EVT和DAT等格式的波形文件解析子模塊組成。由于臺站勘選測試數(shù)據(jù)一般都是直接從數(shù)據(jù)采集器中獲取，本文對國內(nèi)最為通用的CMG-DM24型和EDAS型數(shù)據(jù)采集器記錄的波形文件在軟件中的解析過程進行闡述。其中，港震EDAS數(shù)據(jù)采集器有兩種型號，EDAS-24GN型和EDAS-24IP型，兩種型號記錄的數(shù)據(jù)文件壓縮方式不同。

2.1.1 EDAS-IP型

EDAS數(shù)據(jù)采集器記錄的波形文件需要通過讀取文件頭段的前16個字符判斷數(shù)據(jù)采集器的型號，“SEISMIC EVENT”或“SEISMIC SERIES”為EDAS-24IP型，“GZDAS”則為EDAS-24GN型[7-8]。

EDAS-IP連續(xù)波形的文件由文件頭塊、臺站塊、索引塊和數(shù)據(jù)塊組成。文件結(jié)構(gòu)如圖3所示[7]。

圖3 EDAS-IP型數(shù)據(jù)采集器記錄文件結(jié)構(gòu)Fig.3 Record file structure of EDAS-IP data acquisition device

文件頭塊包括事件頭段、臺網(wǎng)參數(shù)、起始時間和記錄長度。在讀取文件頭塊時主要獲取事件頭段中判斷數(shù)據(jù)采集器型號的事件標(biāo)志、臺網(wǎng)參數(shù)中的臺網(wǎng)代碼和臺站個數(shù)。在讀取臺站塊數(shù)據(jù)時，根據(jù)臺站個數(shù)循環(huán)讀取每個臺站和通道的信息并保存。讀完文件頭塊和臺站塊后，跳過2048字節(jié)大小的索引塊，讀取動態(tài)數(shù)據(jù)塊，包括絕對秒（即數(shù)據(jù)塊的時間長度）、壓縮數(shù)據(jù)塊長度等參數(shù)和壓縮數(shù)據(jù)。動態(tài)數(shù)據(jù)塊的具體讀取流程如圖4所示。

圖4 讀動態(tài)數(shù)據(jù)塊讀流程圖Fig.4 Flow chart of dynamic data block reading

根據(jù)塊長度讀取一個壓縮數(shù)據(jù)塊，采用文獻(xiàn)[9-10]中的解壓方法把壓縮的數(shù)據(jù)解壓為原始數(shù)據(jù)。數(shù)據(jù)塊是按臺站和通道順序存儲的，一個通道的數(shù)據(jù)大小為采樣率乘絕對秒數(shù)。解壓數(shù)據(jù)存入緩存時，按通道順序依次壓入由臺站序號、通道序號、相對時間秒和1秒采樣點數(shù)組成的四維數(shù)組的緩存池。

2.1.2 EDAS-GN型

EDAS-GN型記錄的波形文件總體結(jié)構(gòu)與EDAS-IP類似，由文件頭塊、臺站塊、索引塊和數(shù)據(jù)塊組成，最大的區(qū)別是數(shù)據(jù)塊，它的結(jié)構(gòu)如圖5所示[4]。

圖5 EDAS-GN型數(shù)據(jù)采集器記錄文件數(shù)據(jù)塊結(jié)構(gòu)Fig.5 Data block structure of record file of EDAS-GN data collector

數(shù)據(jù)區(qū)的壓縮方法是，對各通道的N秒數(shù)據(jù)做前后差分diff（n）=x（n）-x（n-1），以10個差分值diff為一組統(tǒng)計有效字長，然后根據(jù)字長在壓縮編碼表中選擇使用的編碼，按照選用的編碼表中的字長對差分值diff值重新編碼輸出。圖6是一個通道的壓縮數(shù)據(jù)CMPDATAFRM的結(jié)構(gòu)[8]。

圖6 CMPDATAFRM數(shù)據(jù)結(jié)構(gòu)Fig.6 Data structure of CMPDATAFRM

解壓過程為：首先計算差分組組數(shù)N=（SP*LEN）/10，SP為采樣率，LEN為當(dāng)前數(shù)據(jù)塊長度（s）；然后讀取一級索引和各差分組的二級索引，根據(jù)一級索引和差分組的二級索引，查文獻(xiàn)[8]中的編碼表，獲得各差分組的編碼長度；最后根據(jù)編碼長度讀取差分值diff（n），并計算當(dāng)前原始值X（n）=diff（n）+X（n-1）。

2.1.3 CMG-DM24型

CMG-DM24型數(shù)據(jù)采集器默認(rèn)的記錄格式為GCF格式。在該儀器配套的采集軟件SCREAM中帶有GCF轉(zhuǎn)MiniSeed格式的子程序gcf2msd。把該子程序移植到地震臺站勘選數(shù)據(jù)處理軟件中，當(dāng)解析GCF格式波形文件時，自動調(diào)用gcf2msd轉(zhuǎn)換子程序，把GCF格式轉(zhuǎn)換為MiniSeed格式，再采用MiniSeed解析模塊進行解析。

MiniSeed是SEED數(shù)據(jù)格式的純數(shù)據(jù)塊，它由若干個長度為4096字節(jié)的數(shù)據(jù)記錄組成，每個數(shù)據(jù)記錄包含控制頭段區(qū)和數(shù)據(jù)區(qū)塊組成。在控制頭段中，前7個字節(jié)是由6個整數(shù)和1個字符“D”組成的數(shù)據(jù)記錄標(biāo)志，通過匹配這個標(biāo)志循環(huán)讀取數(shù)據(jù)記錄進行解壓，在控制頭段信息中，讀取原始數(shù)據(jù)的壓縮格式，按文獻(xiàn)[9]、[10]中的解壓方法把數(shù)據(jù)區(qū)的數(shù)據(jù)解壓為原始數(shù)據(jù)。

2.2 數(shù)據(jù)處理

地震臺站勘選波形數(shù)據(jù)處理的一般過程是：首先把波形截為若干個記錄段，然后每個記錄段按Peterson[11]（1993）的方法計算PSD值，最后采用McNamara等人[12]（2005）提出的概率密度函數(shù)方法計算PDF值和RMS值。本軟件設(shè)計能處理最長24 h的數(shù)據(jù)文件，并將數(shù)據(jù)段進行分段處理，截取長度默認(rèn)為720 s，支持用戶自行設(shè)置，通常長度應(yīng)大于最長記錄周期期的6倍[2]。為使每個記錄段的處理結(jié)果更為平滑，在截取記錄段時以50%的重疊率進行重疊。噪聲功率譜密度（PSD）、噪聲功率譜概率密度（PDF）和噪聲速度有效振幅（RMS）的計算過程如下。

2.2.1 噪聲功率譜密度計算

一個記錄段波形數(shù)據(jù)的PSD值計算流程如下[2，12]：

（1）數(shù)據(jù)預(yù)處理。為了消除對長周期（即大于儀器記錄的最長周期）功率譜估計的影響，對原始波形數(shù)據(jù)去均值并采用平均斜率法消除長周期線性趨勢，處理公式為：

其中x（t）為處理后的數(shù)據(jù)，u（t）為原始數(shù)據(jù)，umean為原始數(shù)據(jù)的算術(shù)平均值，Tlp為長周期線性趨勢，其處理公式為：

式（2）中，Tr為原始數(shù)據(jù)u（t）的總長度。

另外，為了減少有限長度數(shù)據(jù)序列進行FFT變換時造成的頻率滲漏，將錐度為10%的正弦窗函數(shù)應(yīng)用于記錄段數(shù)據(jù)序列，使數(shù)據(jù)段兩端平滑地衰減至零。最后PSD值乘以1.142 857以補償應(yīng)用窗函數(shù)所造成的功率值損失。

（2）速度功率譜密度值計算。采用快速傅立葉變換把原始數(shù)據(jù)變換為以頻率為自變量的函數(shù)Y（f），然后根據(jù)公式（1）計算每個頻點的功率譜密度Pv（f）。

式（3）中，N為采樣點個數(shù)；Δt為采樣時間間隔。

（3）扣除儀器響應(yīng)。為了得到地動噪聲的物理量值，即速度值，功率譜Pv（f）需要進行儀器響應(yīng)校正。

式（4）中，H（f）為儀器的傳遞函數(shù)。

（4）加速度功率譜密度的計算。為了與全球標(biāo)準(zhǔn)低噪聲模型（NLNM）和全球標(biāo)準(zhǔn)高噪聲模型（NHNM）進行比較，采用公式（5）轉(zhuǎn)換為加速度功率譜密度PSDa。

（5）平滑處理。為了得到PSD值在頻率對數(shù)坐標(biāo)中呈等間隔采樣，本文采用1/3倍頻程積分對每條記錄的功率譜密度做平滑處理。

式（6）中，fl為低頻拐角頻率，fh為高頻拐角頻率，n為介于二者之間頻率f的個數(shù)，中心頻率fc以1/9倍頻程為增加步長。

2.2.2 功率譜概率密度計算

所有記錄段PSD值計算結(jié)束后，在-200 dB到-50 dB范圍的功率以1 dB為間隔劃分成連續(xù)的功率窗，統(tǒng)計每個頻點的PSD值落在對應(yīng)功率窗內(nèi)的記錄段個數(shù)，然后根據(jù)概率密度函數(shù)公式（7）計算該頻點處各個功率窗的概率[2]。

式（7）中，Npfc為在中心頻率fc處落在某一功率窗的記錄段個數(shù)，Nfc為總記錄段個數(shù)，Ppsd（fc）為fc頻點處某一功率窗的概率。

2.2.3 噪聲速度有效振幅計算

提取中心頻點在1~20 Hz頻帶范圍內(nèi)所有PSD值，按公式（6）計算各中心頻點RMS值[1-2]。

式（8）中，PSDv為速度功率譜密度，fc為中心頻率，RBW為相對帶寬。

根據(jù)每個中心頻點fc的概率密度函數(shù)公式（9）計算該中心頻點在各個RMS值從-210至-70 dB區(qū)間以1 dB為間隔的區(qū)間的概率[2]。

統(tǒng)計每個中心頻率的概率后，取各中心頻率最大概率處的RMS值求算術(shù)平均值，該值作為該通道數(shù)據(jù)的最后計算RMS值。計算公式為公式（10）[2]。

2.3 GMT繪圖

在GMT繪圖模塊中，主要包括繪制PSD圖、PSD時間曲線圖、PSD時頻圖、PDF圖、RMS柱狀分布圖和RMS時頻圖等繪圖子模塊。每個繪圖子模塊自動生成一個BAT批處理文件，并在批處理文件中寫入繪制相應(yīng)圖件的所有GMT命令。通過QT中的QProcess類調(diào)用執(zhí)行批處理文件實現(xiàn)自動繪圖。

本文中采用的GMT命令模塊為4.5版本[3-4]。在繪制PSD時間曲線圖和RMS柱狀分布圖中，橫軸為時間單位，繪底圖采用psbasemap命令，其中設(shè)置-B選項為“-Bs4Hg4H”，使坐標(biāo)刻度標(biāo)記和網(wǎng)格為4 h。繪制數(shù)值曲線和柱狀圖均采用psxy命令。繪制柱狀圖符號的參數(shù)為“-Sb0.01ub-180”，其中，-Sb為繪制柱狀參數(shù)，0.001u表示柱狀寬度為0.001個坐標(biāo)軸單位，b-180表示柱狀由-180延長到y(tǒng)值。

在繪制PDF圖子模塊中，由于PDF圖是由三維數(shù)據(jù)繪制的二維圖片，PDF數(shù)據(jù)文件的第1列為頻率，第2列為PSD值，第3列為概率，因此，繪制前需要把PDF數(shù)據(jù)進行網(wǎng)格化處理，首先gawk命令把第1列頻率進行對數(shù)計算，然后用blockmedian命令進行插值，使數(shù)據(jù)成為規(guī)則的等間隔數(shù)據(jù)，最后surface命令網(wǎng)格化生成網(wǎng)格文件。繪圖時，先用psbasemap命令繪制出橫軸對數(shù)坐標(biāo)刻度標(biāo)記，然后再采用grdview命令繪制概率彩色圖。

在繪制PSD時頻圖和RMS時頻圖中，數(shù)據(jù)文件均為三維數(shù)據(jù)，第1列數(shù)據(jù)為時間，第2列為頻率，第3列為PSD值或RMS值。在數(shù)據(jù)文件中頻率取數(shù)，然后進行插值網(wǎng)格化生成網(wǎng)格文件。采用psbasemap命令繪制橫軸為時間坐標(biāo)刻度標(biāo)記，縱軸頻率對數(shù)坐標(biāo)刻度標(biāo)記，然后用grdview命令繪制PSD值或RMS值彩色圖。

3 應(yīng)用示例

福建大田均溪測震臺因災(zāi)重新選址新建，在擬新建臺址處安放GURALP公司的CMG-40TDE-40S型地震儀連續(xù)測試48 h。采用地震臺站勘選測試數(shù)據(jù)處理軟件處理測試波形數(shù)據(jù)，獲得PSD圖、5 Hz頻點處PSD值隨時間變化圖、PSD值每個頻點處隨時間變化圖、PDF圖、RMS每小時柱狀分布圖、RMS值每個頻點隨時間變化圖和三個分向RMS總平均值等處理結(jié)果及圖件。

3.1 噪聲加速度功率譜密度分析

噪聲加速度功率譜是分析臺址地震噪聲水平的重要指標(biāo)[5]。本軟件自動產(chǎn)出PSD圖、特定頻點PSD值隨時間變化圖、PDF圖和PSD值每個頻點處隨時間變化圖等四類圖片，從而實現(xiàn)對臺址噪聲加速度功率譜密度進行全面分析。

（1）噪聲加速度功率譜密度圖（PSD圖）如圖7所示，兩根青色曲線是新全球高噪聲模型（NHNM）和新全球低噪聲模型（NLNM），其他三條不同顏色曲線表示三個通道最大概率的PSD值。通過該圖，可以直觀的看到三分向在各頻段內(nèi)噪聲水平，以及與NLNM比較情況。

圖7 噪聲加速度功率譜Fig.7 Power spectral density of noise acceleration

（2）某特定頻點處PSD值隨時間變化圖如圖8所示。在規(guī)范中要求5 Hz頻點處PSD值隨時間變化圖，在寬頻帶或者甚寬頻帶地震計還要5 s頻點和100 s頻點變化圖，本軟件支持產(chǎn)出用戶自定義頻點隨時間變化圖。

圖8 5.0 Hz頻點處PSD值隨時間變化Fig.8 PSD value changes with time at frequency of 5.0 Hz

（3）圖9是垂直分向的加速度功率譜概率密度圖（PDF圖），圖中橫坐標(biāo)為頻率，縱坐標(biāo)是加速度功率譜，顏色塊表示概率。在PDF圖中，通過觀測時間段內(nèi)加速度功率譜的概率分布情況，可以了解非天然地震事件持續(xù)時間占記錄時間的大致比例。

圖9 加速度功率譜概率密度Fig.9 Probability density of acceleration power spectrum

（4）垂直分向上，每頻點PSD值隨時間變化圖如圖10所示，橫坐標(biāo)為時間，以小時為單位，縱坐標(biāo)采用對數(shù)坐標(biāo)表示頻率，用顏色塊表示該頻點該時間點處的PSD值。通過該圖，可以更為全面了解噪聲在時間和頻率上的分布情況。

圖10 每頻點PSD值隨時間變化Fig.10 The PSD value of each frequency point changes with time

3.2 噪聲速度振幅有效值分析

圖11是RMS值每小時的柱狀分布圖，一個小時的RMS值取該小時內(nèi)所有中心頻率處的最大概率RMS值的算術(shù)平均值。圖12是垂直分向RMS值在每頻點隨時間變化圖，橫坐標(biāo)為時間，以小時為單位，縱坐標(biāo)頻率為對數(shù)坐標(biāo)，顏色塊表示RMS值大小。在RMS值柱狀分布圖中，可以很直觀地看到RMS值在每個小時的具體值，以及整個觀測期間RMS是否平穩(wěn)，通過RMS值在每頻點隨時間變化圖，可以進一步了解到噪聲在頻率上的分布情況，為分析噪聲的特征提供了依據(jù)。另外，軟件自動生成每小時RMS值的文本文件，并把所有小時的RMS值的平均值顯示在軟件界面中。

圖11 RMS值每小時柱狀分布圖Fig.11 Hourly histogram of RMS value

圖12 RMS值每頻點隨時間變化圖Fig.12 The RMS value of each frequency point changes with time

4 結(jié)語

地震臺站勘選測試數(shù)據(jù)處理軟件按照規(guī)定的數(shù)據(jù)處理方法實現(xiàn)自動處理勘選測試數(shù)據(jù)，并能自動產(chǎn)出豐富的結(jié)果圖件，解決了地震臺站勘選人員處理勘選測試數(shù)據(jù)流程繁瑣的問題，提高了臺站勘選人員的工作效率，也使處理過程更為規(guī)范化。該套軟件應(yīng)用于福建地震臺網(wǎng)的勘選測試數(shù)據(jù)的處理，自動產(chǎn)出了該臺址地噪聲PSD圖、PDF圖和RMS圖等圖件和數(shù)據(jù)，為評估擬建臺址的地震噪聲水平提供了充足且可靠的依據(jù)，取得了良好的應(yīng)用效果。

地震臺站勘選測試數(shù)據(jù)處理軟件還可以應(yīng)用于其它地震噪聲分析工作。例如，通過處理地震臺站記錄的波形數(shù)據(jù)，繪制臺站的PDF圖，可以分析地震臺站觀測數(shù)據(jù)質(zhì)量，也可以判斷儀器工作狀態(tài)是否正常。