地脈動參數計算入庫程序的開發

王瑣琛　戚浩　夏仕安

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

【關鍵詞】地脈動；數據庫；算法

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 引言

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

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

1 程序的架構

本程序的功能是實時接收地脈動波形數據，利用這些數據實時地計算出地脈動參數并將其存入數據庫。為實現實時接收地脈動數據，程序需要調用winsock建立與流服務器的連接，而與流服務器之間的通信則是通過NetSeis/IP協議；為了計算地脈動參數，程序需要通過Matlab引擎調用Matlab命令進行相關的計算；為了將地脈動程序寫入數據庫，程序需要通過C++下的MySQL數據庫接口實現對MySQL數據庫的操作。圖1顯示了該程序的架構。

圖1 地脈動參數計算入庫程序的架構

從圖1中可以看出，主程序主要由兩個進程構成。數據接收進程不間斷地運行，負責實時接收地脈動波形數據并存入本地內存；計算入庫進程每隔一段時間運行一次，調取本地內存中的地脈動波形數據，再調用Matlab程序進行計算求出地脈動參數并保存到數據庫中，完成上述步驟后計算入庫進程進入休眠狀態等待一定時間間隔后再次運行。下面按照功能的區分詳細介紹該程序各部分的開發。

2 數據的接收

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

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）；

}

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

3 算法介紹

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

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

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

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

int i=start；

while（i

{

//下面求出各分量數據總和。

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

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

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

i=i+1；

}

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

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

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

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

//最后對三分量的數據進行去偏移量。

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；

}

數據經過程序上述步驟的處理，被消去了偏移量。下面用經過去偏移量處理的數據構造一個協方差矩陣，代碼如下：

//矩陣A為所求的協方差矩陣。

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；

}

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

//求方位角angle。

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

anlge=90-angle；

while（angle<0）

angle=angle+360；

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

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

//求視出射角iangle。

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

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

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

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

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

//求空間線性度alpha。

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

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

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

由于在求解方位角時已經對數據進行過去偏移量的處理，在這里只需要再將數據除以均方差即可完成歸一化處理。為了最大限度地保持原始數據的信息量，沒有使用更多的處理方法。選4096個數據點采用FFT方法直接求功率譜，即先求數據點的FFT交換，得到波形的傅里葉譜，再從傅里葉譜求的波形的功率譜。將處理好的數據Yi傳入Matlab，保存在數組Y中，依次調用如下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為高頻衰減系數。

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

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

4 地脈動參數的入庫

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

表1 存儲地脈動參數的表中的所有字段

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

表2 地脈動參數類型代碼的含義

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

5 程序的界面以及使用

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

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

6 結束語

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

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

【參考文獻】

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

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

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

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

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

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

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

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

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

[10]楊迪雄，王偉.近斷層地震動的頻譜周期參數和非平穩特征分析[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.