王瑣琛　戚浩　夏仕安

【摘 要】地脈動資料的研究相對于地震觀測資料的研究尚有很大開發(fā)空間，程序開發(fā)人員試圖開發(fā)出一款能夠?qū)Φ孛}動研究起到幫助作用的實用程序，因而開發(fā)了地脈動參數(shù)計算入庫程序。本文介紹了地脈動參數(shù)計算入庫程序的開發(fā)，詳細(xì)介紹了該程序的架構(gòu)以及各部分功能的實現(xiàn)，包括數(shù)據(jù)流的接收、地脈動參數(shù)計算的算法實現(xiàn)、地脈動參數(shù)的入庫以及程序界面的設(shè)計，其中著重介紹了地脈動參數(shù)計算在編程時的算法的實現(xiàn)。在本文的結(jié)尾闡述了該程序的使用、維護(hù)情況以及該程序的開發(fā)心得和積極意義。

【關(guān)鍵詞】地脈動；數(shù)據(jù)庫；算法

The Development of Calculating and Storage of Microtremor Program

WANG Suo-chen QI Hao XIA Shi-an

（Seismological Bureau of Anhui Province， Hefei Anhui 230031， China）

【Abstract】Compared with the study of seismic observations，the study of microtremor data has considerable development space now.Program developers trying to develop a practical program which is able to play a helpful role in the research of microtremor，so we developed calculating and storage of microtremor program.This article describes the development of calculating and storage of microtremor program，details the architecture of this program and the realization of the function of each part of the program，including receiving a data stream，the implementation of microtremor parameter calculation algorithm，the storage of parameter，and the design of the program interface，which focuses on the implementation of microtremor parameter calculation algorithm.At the end of this article it describes the use of the program，as well as to maintain the program and the development experience and the positive significance.

【Key words】Microtremor； Database； Algorithm

0 引言

隨著“九五”、“十五”測震臺網(wǎng)的建設(shè)，各臺網(wǎng)已經(jīng)積累了大量的波形資源，這些波形包括地震事件波形，但更多的是無震的地脈動數(shù)據(jù)。地震事件波形包含了豐富的震源、地球介質(zhì)以及臺站場地響應(yīng)等信息，已經(jīng)有許多科技人員利用它做震源參數(shù)計算、地下速度結(jié)構(gòu)反演等大量的研究工作，發(fā)揮了很大的作用。而地脈動資料一直被塵封在柜子里，很少被人問津。近幾年也有一些專家有針對性的選取了幾個地脈動指標(biāo)，通過編寫程序進(jìn)行計算，并對計算結(jié)果進(jìn)行了初步分析，研究結(jié)果顯示該項工作非常有意義。但由于沒有一個系統(tǒng)的數(shù)據(jù)處理系統(tǒng)，使得該項工作處于停滯階段。為了更大地發(fā)揮臺網(wǎng)記錄效益，最大限度利用數(shù)字記錄信號，程序開發(fā)人員準(zhǔn)備開發(fā)一套地脈動參數(shù)計算入庫程序，為地震預(yù)報研究人員提供更多的數(shù)據(jù)產(chǎn)品。

相當(dāng)多的地震學(xué)家認(rèn)為，在地震發(fā)生前有一個應(yīng)力在未來震源區(qū)集中的過程，即孕震過程。當(dāng)這一過程發(fā)展到一定階段時，孕震區(qū)內(nèi)的巖石可能會出現(xiàn)微破裂和塑性化、塑性硬化、相變等一系列現(xiàn)象，從而導(dǎo)致地震射線傳播路徑的變化，引起地震波出射角和方位角的變化，出射角、方位角這些異常變化有可能作為地震預(yù)報的一種指標(biāo)。此外，孕震區(qū)內(nèi)的震源動力學(xué)參數(shù)的變化也可能引起地脈動記錄的某些變化，根據(jù)地脈動的頻譜變化異常也可以探索預(yù)測中強(qiáng)震發(fā)生。因此實現(xiàn)上述地脈動參數(shù)的計算入庫就是該軟件的目標(biāo)。

1 程序的架構(gòu)

本程序的功能是實時接收地脈動波形數(shù)據(jù)，利用這些數(shù)據(jù)實時地計算出地脈動參數(shù)并將其存入數(shù)據(jù)庫。為實現(xiàn)實時接收地脈動數(shù)據(jù)，程序需要調(diào)用winsock建立與流服務(wù)器的連接，而與流服務(wù)器之間的通信則是通過NetSeis/IP協(xié)議；為了計算地脈動參數(shù)，程序需要通過Matlab引擎調(diào)用Matlab命令進(jìn)行相關(guān)的計算；為了將地脈動程序?qū)懭霐?shù)據(jù)庫，程序需要通過C++下的MySQL數(shù)據(jù)庫接口實現(xiàn)對MySQL數(shù)據(jù)庫的操作。圖1顯示了該程序的架構(gòu)。

圖1 地脈動參數(shù)計算入庫程序的架構(gòu)

從圖1中可以看出，主程序主要由兩個進(jìn)程構(gòu)成。數(shù)據(jù)接收進(jìn)程不間斷地運(yùn)行，負(fù)責(zé)實時接收地脈動波形數(shù)據(jù)并存入本地內(nèi)存；計算入庫進(jìn)程每隔一段時間運(yùn)行一次，調(diào)取本地內(nèi)存中的地脈動波形數(shù)據(jù)，再調(diào)用Matlab程序進(jìn)行計算求出地脈動參數(shù)并保存到數(shù)據(jù)庫中，完成上述步驟后計算入庫進(jìn)程進(jìn)入休眠狀態(tài)等待一定時間間隔后再次運(yùn)行。下面按照功能的區(qū)分詳細(xì)介紹該程序各部分的開發(fā)。

2 數(shù)據(jù)的接收

要計算地脈動參數(shù)就需要地脈動波形數(shù)據(jù)，本程序要實現(xiàn)地脈動參數(shù)的實時計算，因此采取了從數(shù)據(jù)流服務(wù)器接收實時的地脈動數(shù)據(jù)的辦法。具體做法是先用winsock提供的connect函數(shù)建立程序與流服務(wù)的連接，再利用send函數(shù)向服務(wù)器發(fā)送命令，在取得接收許可后開啟新的連接并利用recv函數(shù)接收地脈動數(shù)據(jù)。具體代碼如下：

SOCKET Pclient_1；

Pclient_1=socket（AF_INET，SOCK_STREAM，IPPROTO_IP）；

SOCKADDR_IN serveraddr_1；

serveraddr_1.sin_family=AF_INET；

serveraddr_1.sin_port=htons（5000）；

serveraddr_1.sin_addr.s_addr=inet_addr（serverIp）；

connect（Pclient_1，（struct sockaddr*）&serveraddr_1，sizeof（serveraddr_1））；

CString str_input="user"+serverUser+"＼n"；

send（Pclient_1， str_input.GetBuffer（0），str_input.GetLength（），0）；

str_input ="pass"+serverPass+"＼n"；

send（Pclient_1，str_input.GetBuffer（0），str_input.GetLength（），0）；

char ch_get_data[512]；

memset（ch_get_data，0，sizeof（ch_get_data））；

recv（Pclient_1，ch_get_data，512，0）；

str_input ="pasv rt＼n"；

send（Pclient_1，str_input.GetBuffer（0），str_input.GetLength（），0）；

memset（ch_get_data，0，sizeof（ch_get_data））；

recv（Pclient_1，ch_get_data，512，0）；

str_input ="retr"+strSta+"＼n"；

send（Pclient_1，str_input.GetBuffer（0），str_input.GetLength（），0）；

int port1，port2；

sscanf（（char*）ch_get_data，"227Real Time Data Port Entering Passive Mode（%*d，%*d，%*d，%*d，%d，%d）"，&port1，&port2）；

SOCKET Pclient_2；

Pclient_2=socket（AF_INET，SOCK_STREAM， IPPROTO_IP）；

SOCKADDR_IN serveraddr_2；

serveraddr_2.sin_family=AF_INET；

serveraddr_2.sin_port=htons（port1*256+port2）；

serveraddr_2.sin_addr.s_addr=inet_addr（serverIp）；

connect（Pclient_2，（struct sockaddr*）&serveraddr_2，sizeof（serveraddr_2））；

while（1）

{

memset（ch_get_data，0，sizeof（ch_get_data））；

recv（Pclient_2，ch_get_data，512，0）；

}

上述代碼實現(xiàn)了程序不斷從流服務(wù)器接收地脈動數(shù)據(jù)的功能，這些地脈動數(shù)據(jù)是以512字節(jié)為單位的MiniSEED格式存儲的。要進(jìn)行地脈動參數(shù)的計算，程序要先解析MiniSEED格式的數(shù)據(jù)，提取其中的臺站信息以及地脈動波形數(shù)據(jù)。程序?qū)⑦@些地脈動數(shù)據(jù)轉(zhuǎn)換成int型變量，臨時存儲在一個2維數(shù)組中以備調(diào)用。該二維數(shù)組的一個維度表示不同的時刻，另外一個維度表示不同的臺站以及不同的通道。當(dāng)內(nèi)存中保存有用于計算地脈動參數(shù)的波形數(shù)據(jù)時，程序就將使用這些數(shù)據(jù)計算出地脈動參數(shù)，下面將介紹程序是如何用地脈動波形數(shù)據(jù)計算出地脈動參數(shù)的。

3 算法介紹

地脈動參數(shù)的計算方法是本程序的核心，該軟件能夠計算的地脈動參數(shù)包括：方位角、視出射角、三個分量的時間線性度以及對應(yīng)的平均半周期、空間線性度（？琢1和？琢2）以及頻譜參數(shù)，其中頻譜參數(shù)包含三個分量的卓越頻率、最大譜值、最低頻率、最小譜值、高頻衰減系數(shù)以及頻譜衰減段斜線擬合程度。

方位角的計算使用了偏振分析法，該方法的原理是構(gòu)造地脈動振動形成的偏振橢球，并將其用一個協(xié)方差矩陣來表示，求出該矩陣的特征值和特征向量進(jìn)而求出地脈動的方位角。下面介紹該方法在程序中的具體實現(xiàn)步驟。先截取需要進(jìn)行地脈動參數(shù)計算的地脈動波形數(shù)據(jù)。為方便計算將得到的數(shù)據(jù)的三分量進(jìn)行去偏移量，代碼示例如下：

//start為數(shù)據(jù)起始點，end為數(shù)據(jù)結(jié)束點。Data.EW[i]、Data.NS[i]、Data.UD[i]為三個分量的地脈動數(shù)據(jù)。

double sum.EW=0，sum.NS = 0，sum.UD=0；//這三個變量先用來存儲三分量數(shù)據(jù)的求和的值再轉(zhuǎn)化為偏移量。

int i=start；

while（i

{

//下面求出各分量數(shù)據(jù)總和。

sum.EW=sum.EW+Data.EW[i]；

sum.NS=sum.NS+Data.NS[i]；

sum.UD=sum.UD+Data.UD[i]；

i=i+1；

}

//下面求出各分量數(shù)據(jù)的偏移量，其中（end-start）為每個通道數(shù)據(jù)的個數(shù)。

sum.EW=sum.EW/（end-start）；

sum.NS=sum.NS/（end-start）；

sum.UD=sum.UD/（end-start）；

//最后對三分量的數(shù)據(jù)進(jìn)行去偏移量。

i=start；

while（i

{

Data.EW[i]=Data.EW[i]-sum.EW；

Data.NS[i]=Data.NS[i]-sum.NS；

Data.UD[i]=Data.UD[i]-sum.UD；

i= i+1；

}

數(shù)據(jù)經(jīng)過程序上述步驟的處理，被消去了偏移量。下面用經(jīng)過去偏移量處理的數(shù)據(jù)構(gòu)造一個協(xié)方差矩陣，代碼如下：

//矩陣A為所求的協(xié)方差矩陣。

i=start；

while（i

{

A[0][0]=A[0][0]+Data.EW[i]*Data.EW[i]/（end-start）；

A[0][1]=A[0][1]+Data.EW[i]*Data.NS[i]/（end-start）；

A[0][2]=A[0][2]+Data.EW[i]*Data.UD[i]/（end-start）；

A[1][0]=A[0][1]；

A[1][1]=A[1][1]+Data.NS[i]*Data.NS[i]/（end-start）；

A[1][2]=A[1][2]+Data.NS[i]*Data.UD[i]/（end-start）；

A[2][0]=A[0][2]；

A[2][1]=A[1][2]；

A[2][2]=A[2][2]+Data.UD[i]*Data.UD[i]/（end-start）；

i=i+1；

}

協(xié)方差矩陣A構(gòu)造完畢后，求出該矩陣的特征值和特征向量，進(jìn)而求出所選地脈動數(shù)據(jù)的方位角?？梢允褂肕atlab軟件直接求出矩陣A的特征值和特征向量，將矩陣A傳入Matlab中，調(diào)用Matlab命令“[V，D]=eig（A）”，得到矩陣A的全部特征值，構(gòu)成對角陣D，而A的特征向量構(gòu)成V的列向量。用矩陣V求方位角的代碼如下：

//求方位角angle。

angle=atan2（V[0]，V[1]）*180/PI；//PI為圓周率的取值。

anlge=90-angle；

while（angle<0）

angle=angle+360；

至此，地脈動參數(shù)之一的方位角得以求出。

視出射角的計算方法也是采用偏振分析法，并且與計算方位角一樣，都需要解出特征向量矩陣V。在求取方位角并計算出特征矩陣V后，用下面的代碼可以直接求出視出射角。

//求視出射角iangle。

iangle=acos（-V[2]）*180/PI；//PI為圓周率的取值。

下面介紹時間線性度和平均半周期的計算方法。對于某一分量上一段時間長度的地脈動波形數(shù)據(jù)，取其過零線（平衡點）的一系列時間值：

將每個過零點的數(shù)據(jù)的時刻存儲在一個數(shù)組X中，將這些點的序號存儲在另一個數(shù)組Y中。將數(shù)組X，Y傳入Matlab，調(diào)用“R=corrcoef（X，Y）”命令，可以求出兩組數(shù)據(jù)的線性相關(guān)系數(shù)并保存在R中，調(diào)用“A=polyfit（X，Y，1）”命令可以求出兩組數(shù)據(jù)線性擬合的斜率并保存在A中。R與A即時間線性度？酌與平均半周期Th。

空間線性度表示地脈動實際射線軌跡和初動射線的符合程度?？臻g線性度的計算需要用到之前求方位角時用到的協(xié)方差矩陣A的特征值矩陣D。

對于完全均勻的介質(zhì)有？琢1=？琢2=1；介質(zhì)非均勻性越強(qiáng)，？琢1和？琢2的值越小，程序中計算空間線性度的代碼如下：

//求空間線性度alpha。

alpha[0]=1-D[1]/D[0]；

alpha[1]=1-D[2]/D[0]；

頻譜參數(shù)的計算較為復(fù)雜，在編程實現(xiàn)頻譜參數(shù)的計算時，Matlab軟件發(fā)揮了重要的作用。首先為了去掉地脈動波形中的直流分量先進(jìn)行歸一化處理，具體做法是先對數(shù)據(jù)進(jìn)行去偏移量，再把每一個數(shù)據(jù)除以該段數(shù)據(jù)的均方差。

由于在求解方位角時已經(jīng)對數(shù)據(jù)進(jìn)行過去偏移量的處理，在這里只需要再將數(shù)據(jù)除以均方差即可完成歸一化處理。為了最大限度地保持原始數(shù)據(jù)的信息量，沒有使用更多的處理方法。選4096個數(shù)據(jù)點采用FFT方法直接求功率譜，即先求數(shù)據(jù)點的FFT交換，得到波形的傅里葉譜，再從傅里葉譜求的波形的功率譜。將處理好的數(shù)據(jù)Yi傳入Matlab，保存在數(shù)組Y中，依次調(diào)用如下Matlab命令：

Y=fft（Y）；

p=Y.*conj（Y）/N；

F=100*（0：（N/2-1））/ N；

P=p（1：N/2）；

[a1，b1]=max（P）；//其中a1為最大譜值。

c1=F（b1）；//c1為卓越頻率。

a=P>（a1*0.707）；

[a2，b2]=max（a）；

c2=P（b2）；//c2為最小譜值。

d2=F（b2）；//d2為最低頻率。

z=F>0.61 & F<2.00；

lgs=log10（P（z>0））；

lgf=log10（F（z>0））；

R=corrcoef（lgf，lgs）；//R為高頻衰減系數(shù)。

A=polyfit（lgf，lgs，1）；//A為頻譜衰減段斜線擬合程度。

至此，所有地脈動參數(shù)全部計算得出，接下來程序要做的就是將計算得到的地脈動數(shù)據(jù)儲存起來以用于分析研究。

4 地脈動參數(shù)的入庫

本程序使用MySQL數(shù)據(jù)庫存儲地脈動參數(shù)。對于某一臺站記錄到的一段三分量的地脈動數(shù)據(jù)，共有28個地脈動參數(shù)，分別為：方位角與視出射角（共2個）、每個分向的時間線性度和平均半周期（共6個）、空間線性度？琢1和？琢2（共2個）、頻譜參數(shù)每個分向各6個（共18個）。對于單個地脈動參數(shù)，需要包含以下信息：所屬臺網(wǎng)、臺站、通道名、參數(shù)類型、參數(shù)的值、參數(shù)的單位、參數(shù)被計算出的時刻。綜合上述考慮，可將所有地脈動參數(shù)存儲在同一張表內(nèi)，用不同的字段存儲地脈動數(shù)據(jù)包含的不同信息。表1列出了存儲地脈動參數(shù)的表中的所有字段。

表1 存儲地脈動參數(shù)的表中的所有字段

其中用于存儲地脈動參數(shù)的類型的Type字段中的信息用字母代碼表示，這是為了方便數(shù)據(jù)的查詢。表2顯示了地脈動參數(shù)類型代碼的含義。

表2 地脈動參數(shù)類型代碼的含義

程序在計算出地脈動參數(shù)后會將地脈動數(shù)據(jù)信息按照表的格式存入數(shù)據(jù)庫。

5 程序的界面以及使用

圖2展示了本程序的界面。

程序界面的左邊是臺站管理區(qū)域，用于添加或者刪除要接收數(shù)據(jù)并進(jìn)行地脈動參數(shù)計算入庫的臺站。程序界面的中間為數(shù)據(jù)庫配置區(qū)域，在這個區(qū)域輸入連接數(shù)據(jù)庫的參數(shù)，以及數(shù)據(jù)庫表名，計算得出的地脈動參數(shù)將保存于該數(shù)據(jù)庫中。程序界面的右側(cè)上部分是流服務(wù)配置區(qū)域，在該區(qū)域輸入用于接收數(shù)據(jù)流的服務(wù)器的參數(shù)，程序?qū)⑦B接至該服務(wù)器并從該服務(wù)器接收地脈動數(shù)據(jù)。程序界面的右側(cè)下部分是地脈動參數(shù)計算的相關(guān)參數(shù)的控制區(qū)域，在該區(qū)域中可以設(shè)置計算窗長的長度以及處理間隔的時長。窗長的含義是用于計算地脈動參數(shù)時所選取的地脈動數(shù)據(jù)的長度，以秒為單位計算。處理間隔是指程序每兩次計算地脈動參數(shù)之間的時間間隔，以分鐘為單位計算。上述參數(shù)中的流服務(wù)器參數(shù)、臺站信息、數(shù)據(jù)庫參數(shù)都同步保存在程序根目錄下config文件夾中的配置文件config.txt中。在配置好各項參數(shù)后，點擊右下方的“開始處理”按鈕，程序就會從流服務(wù)器接收地脈動數(shù)據(jù)并計算出地脈動參數(shù)，再將計算出的地脈動參數(shù)存入數(shù)據(jù)庫中。

6 結(jié)束語

經(jīng)過了2個月的試運(yùn)行觀察，地脈動參數(shù)計算入庫程序運(yùn)行穩(wěn)定，能夠保證地脈動參數(shù)實時地計算入庫。地脈動數(shù)據(jù)的存儲方式經(jīng)檢驗也是可靠而有效的，用于存儲地脈動數(shù)據(jù)的MySQL數(shù)據(jù)庫按照每分鐘計算一次的地脈動參數(shù)數(shù)據(jù)量，在二個月內(nèi)積累了2.42GB。在對該數(shù)據(jù)庫進(jìn)行查詢操作時，數(shù)據(jù)庫運(yùn)行流暢，能夠按照查詢命令給出的查詢范圍快速地給出所需要的地脈動參數(shù)以便用于分析研究。綜上所述該程序的開發(fā)實現(xiàn)了預(yù)期的功能。

地脈動參數(shù)計算入庫程序的開發(fā)，使得安徽測震臺網(wǎng)的大量地脈動數(shù)據(jù)得到了更充分的利用，為地震預(yù)測預(yù)報研究工作提供了更加豐富的資料。同時在開發(fā)該程序時運(yùn)用到的一些編程技術(shù)和算法——例如數(shù)據(jù)流的實時接收技術(shù)、SEED格式的解析、C++程序與Matlab程序的聯(lián)合開發(fā)、C++程序與MySQL數(shù)據(jù)庫系統(tǒng)的聯(lián)合開發(fā)等技術(shù)以及各類地脈動參數(shù)的算法也為安徽測震臺網(wǎng)今后的科研起到了推進(jìn)作用。

【參考文獻(xiàn)】

[1]和躍時.利用數(shù)字地震記錄計算震中方位角方法簡介[J].東北地震研究，2004，20（1）：48-53.

[2]馬亮，盧建旗，朱敏，郭曉云.偏振分析法計算震中方位角的精度分析[J].防災(zāi)科技學(xué)院學(xué)報，2014，16（1）：36-47.

[3]戚浩，黃顯良，汪貴章，夏仕安，張炳，程鑫.安徽地震臺網(wǎng)地脈動處理系統(tǒng)的建設(shè)[J].地震地磁觀測與研究，2011，32（3）：161-163.

[4]李志雄，袁錫文，朱航，丁軍，陳金燕，明穗花.海南數(shù)字地脈動參數(shù)處理系統(tǒng)[J].地震，2008，28（2）：87-92.

[5]王梅德，韓艷杰，郭祥云，于仁寶，賈炯，趙祖虎，鞠勇.地脈動在大震前的異常變化研究[J].地震研究，2014，37（1）：74-78.

[6]夏英杰，倪四道，曾祥方.汶川地震前地脈動信號的單臺法研究[J].地球物理學(xué)報，2011，54（10）：2591-2596.

[7]劉東旺，凌學(xué)書，何小偉，童國林.地震波時間線性度在安徽及鄰區(qū)中強(qiáng)震中短期預(yù)報中的應(yīng)用研究[J].地震地磁觀測與研究，2000，21（6）：49-57.

[8]陳化然，馮德益.大震前地磁場時空間線性度變化的初步研究[J].地震地磁觀測與研究，1994，15（5）：22-27.

[9]鄧亞虹，徐平，李喜安，范文.層狀場地脈動卓越頻率對應(yīng)的理論計算深度研究[J].振動工程學(xué)報，2015，28（3）：366-373.

[10]楊迪雄，王偉.近斷層地震動的頻譜周期參數(shù)和非平穩(wěn)特征分析[J].地震工程與工程振動，2009，29（5）：27-35.

[11]Kelli Hardesty，Lorraine W.Wolf，Paul Bodin.Noise to signal：A microtremor study at liquefaction sites in the NewMadrid Seismic Zone[J].Geophysics：Journal of the Society of Exploration Geophysicists，2010，75（3）：B83.