王琐琛　戚浩　夏仕安

【摘 要】地脉动资料的研究相对于地震观测资料的研究尚有很大开发空间，程序开发人员试图开发出一款能够对地脉动研究起到帮助作用的实用程序，因而开发了地脉动参数计算入库程序。本文介绍了地脉动参数计算入库程序的开发，详细介绍了该程序的架构以及各部分功能的实现，包括数据流的接收、地脉动参数计算的算法实现、地脉动参数的入库以及程序界面的设计，其中着重介绍了地脉动参数计算在编程时的算法的实现。在本文的结尾阐述了该程序的使用、维护情况以及该程序的开发心得和积极意义。

【关键词】地脉动；数据库；算法

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构造完毕后，求出该矩阵的特征值和特征向量，进而求出所选地脉动数据的方位角。可以使用Matlab软件直接求出矩阵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。

空间线性度表示地脉动实际射线轨迹和初动射线的符合程度。空间线性度的计算需要用到之前求方位角时用到的协方差矩阵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.