2014, Vol. 57
2. 中国地震局地震监测与减灾技术重点实验室, 广州 510070;
3. 地震预警和重大工程安全诊断广东省重点实验室, 广州 510070
2. Key Laboratory of Earthquake Monitoring and Disaster Mitigation Technology of CEA, Guangzhou 510070, China;
3. Key Laboratory of Earthquake Early Warning and Safety Diagnosis of Major Projects of Guangdong Province, Guangzhou 510070, China
地震学研究的最基础数据对象就是地震记录数据.地震记录一定程度上反映了观测点处的实际地面运动.由于实际地震观测系统受频率特性的限制，记录数据所反映的地面运动情况是不完全的、存在畸变的.为了研究真实的地面运动情况，需要对记录数据做扣除仪器响应也就是传递函数反褶积的处理.而为了研究一些特定的地动信号，又可能需要进一步限制记录数据的频带范围.这就决定了由记录数据恢复真实地面运动及仿真不同仪器系统记录的数据处理是一项既基础又具实用价值的工作.
用于扣除仪器响应的处理方法曾被Aki和Richards(1980)归纳为有四种，其中两种完全在时间域中进行，另两种需要变换到频率域中进行.从仪器运动方程直接积分的方法，由于受记录数据中存在干扰与误差的影响，必须对记录数据的分布规律、积分初始条件附加限制条件和做出一些简化的假设 (王广福,1986).频率域的方法计算效率高，但存 在易产生高频和甚低频频谱失真的问题(Aki & Richards, 1980； 李鸿吉，1992a,1992b)，李鸿吉，1992a,1992b)提出“引用周期”的方法修改频谱以减少仪器校正带来的失真，并用现代控制论解状态方程的方法仿真仪器记录和验证仪器校正结果，取得了较好的效果，确实改善了地面位移的估计.李鸿吉的方法效果虽好，但算法复杂不利于在例行常规处理工作中应用，刘瑞丰等(1997)仍然沿用直接的Fourier分析方法进行长周期与中长周期地震仪记录的仿真.
Aki和Richards(1980)列举的第四种校正仪器响应方法是递归滤波法，他们认为当仅需要少量的递归滤波系数就能达到给定目的时，该方法比FFT方法更有效.实际情况并不简单，因为真实的地震观测系统不是少量的递归滤波系数就能完整描述的.但是，当我们研究的问题可以主要考虑地震计的低频特性和忽略高频响应带来的误差时，用少量递归滤波系数描述仪器系统就成为可能.针对机械惯性地震计模型，王广福(1986)给出了这种情况下的一种二重递归滤波方法，该方法得到的是地震计质量块位移对地面运动位移的响应，但很容易证明其结果适用于速度响应平坦地震计对地面运动速度响应的情形(见下文2节)，可实现从现代速度记录恢复地面运动速度过程及反过来仿真速度记录.由于该方法包含与地面运动数据相关的递归系数有三个，其中两个涉及历史数据，因此在进行恢复地面运动计算时需要假设两个地面运动初始值，而进行仿真计算时更是需要假设三个地面运动初始值.
一些传统的模拟记录地震仪包含电子放大器及电流计记录仪等部件，整体传递函数比单一的机械惯性地震计模型复杂，但有的实际工作(例如震级计算)又需要将现代地震记录数据仿真成这些传统仪器的记录，王洪体等(2006)为此研究了用双线性变换方法构造出仿真DD-1、DK-1地震仪的数字递归滤波器.
随着强震动观测技术的迅速发展，强震记录数据处理方法也得到迅速发展，包括计算各类地震动反应的仿真方法.金星等(2003)对众多仿真方法中较有代表性的4种做了详细研究，据此他们(2004)提出自己的递归滤波公式，可实时仿真地动速度，尝试应用于测震数据处理系统.
本文用另一种递归滤波方法实现地震记录的反褶积与仿真，仅一个递归系数与地面运动数据有关，反褶积计算不需要假设地面运动初始值，仿真计算则只需要假设一个地面运动初始值，本文同时给出 该方法的数值实验结果及应用实例.本方法是对速度记录的仿真，反褶积结果的物理量与前人的工作不同.
机械惯性地震计质量块对地面运动位移的传递函数为(王广福, 1986)：
式中ωs=2π/Ts,Ts是地震计周期,Ds是阻尼.王广福(1986)通过对此式做双线性变换得到一个共6项系数的二重递归滤波器，其中需要进行计算的系数项有4个.当考虑速度响应地震计时，地震计输出正比于质量块运动速度u(t)=Edy/dt，其拉普拉斯变换U(s)=EsY(s)，因此速度响应地震计对地面运动速度v(t)的传递函数成为：
与式(1)形式相同，只是多了个灵敏度因子E，因此王广福给出的递归滤波器同样适用于速度记录与地面运动速度间的仿真及反褶积.
我们再来考虑一下地面运动加速度的时间变化率，记之为w(t)=d3x/dt3，其拉普拉斯变换与X(s)的关系是W(s)=s3X(s).因此速度地震计对w(t)的传递函数可以表示为：
对式(3)做分式展开得到：
式中
对式(4)做冲击不变法变换(奥本海姆和谢弗, 1980)得到：
Hw(z)=
式中Δt是时域离散系统采样间隔.将s1、s2代入并化简得到：
式中:
该递归滤波器总共4个系数项，需要进行计算的系数项有3个，分别为a、b和c，仅a与地面运动相关，b和c与输出相关.由于冲击不变法变换本身会带来频谱的混叠(奥本海姆和谢弗, 1980)，应用式(6)的递归滤波器时要求记录数据是充分限带的，即地震计本身限带或数采带有适当的低通滤波器.本文用到的几种地震数据采集器符合此条件、具有低通滤波特性，且阻带衰减较大(>130 dB)，所以在应用本节 递归滤波方法的数据处理过程中不会发生频率混叠.
此外，式(6)的应用还需要注意其物理意义，用它对速度记录数据进行反褶积计算得到的量对应地动加速度的时间变化率，而不是我们通常处理的地动位移、速度或加速度.
下面从数值实验、实际事件数据仿真及地震仪方位角相对测量三个方面的应用实例说明本方法的适用性和有效性.全部实验所用地震计的型号、参 数、生产厂商及仪器出厂序列号等信息统一列于表1.
|
表1 地震计信息 Table 1 Seismometer information |
根据地震计标定的原理，标定过程相当于地震计感应到正比于标定源信号的等效地动加速度运动过程.我们对标定响应记录数据应用2节描述的递归滤波器进行反褶积处理将得到对应标定信号时间变化率的数据序列.对反褶积的结果数值积分一次应该得到与标定信号对应的数据序列.当然，实际数据中还包含有环境地动噪声的成分.
伪随机二进制序列(PRBS)标定方法是现代常用的地震计电标定方法之一(Berger et al., 1979).由于PRBS具有接近白噪声的频谱特性而且易于实现，用该方法可实测出全频带完整的地震计传递函数.PRBS序列在码元0和1之间产生信号突变，而其他时间信号不变，求其时间微商应该在PRBS码元0和1之间变化处是δ函数，因此PRBS序列信号应该是我们对2节描述方法做数值实验的一种较好的选择.
这个数值实验的基本内容是利用甚宽带地震计的PRBS标定响应数据做反褶积处理，得到一个等效于地动加速度时间变化率过程的序列，再以反褶积结果作为输入信号仿真宽带地震计的标定输出，最后与真实宽带地震计的PRBS标定响应数据进行对比.
图1是作者将表1序号为1的甚宽带地震计连接到港震公司自己的EDAS-24IP3数采，实施PRBS标定并做反褶积、数值积分处理的结果，为了图示清晰，只画出数据的初始1000个样点.图中原始PRBS信号是直接对数采标定输出采样取得的.明显看出，反褶积的结果是一系列接近于δ函数的尖锐脉冲，与预期的结果相符，而对其积分一次所得恢复的PRBS与原始PRBS信号也符合得很好，只是有高频成分缺失，这是反褶积不包括地震计的高频特性及数采的低通滤波器造成的.
 | 图1 PRBS标定反褶积处理结果 Fig.1 Deconvolution results of PRBS calibration response |
作者选择了实际运行中的广东区域地震台网两个台站的设备进行仿真实验，地震计情况见表1序号2组合，数采型号同为港震公司EDAS-24IP3.对这两个台站实施了参数相同的PRBS标定：PRBS信号代码bin10r，码元宽度200 ms，信号幅度900 counts，序列输出2次.然后对120 s甚宽带摆的标定响应反褶积并仿真成60 s宽带摆的响应.
图2a显示了实测及仿真计算的波形结果，仿真的标定响应与实测的BBVS-60响应符合得很好，标定信号结束后的暂态响应不同是两台被测仪器标定电路的偏置电流不一致所造成的.图2b则是实测及仿真数据的频谱结果，仿真的频谱与实测的BBVS-60响应频谱也基本重合.反褶积与仿真的数值实验结果验证了2节给出的递归滤波方法是适用的，物理意义明确，可与实测信号对应.
 | 图2 PRBS标定实测与仿真结果 (a) PRBS标定BBVS仿真与实测响应；(b) BBVS仿真与实测响应频谱. Fig.2 Real and simulated results of PRBS calibration (a) BBVS real & simulated response of PRBS calibration; (b) Spectra of BBVS real & simulated response. |
广州测震台安装有JCZ-1T超宽带地震计(见表1序号3的组合)，数采是港震公司的EDAS-24IP6.DS-4A短周期地震计是临时安装的，接珠海泰德公司的TDE-324CI数采.两套仪器同时记录到 了2012年8月31日菲律宾附近海域的M7.3地 震，作者从JCZ-1T的BB通道记录数据用本文第2节的递归滤波方法仿真DS-4A记录，UD向结果见图3，除数采输入量程不同外，仿真结果与实际短周期记录一致.
 | 图3 地震事件实测与仿真比较 (a) JCZ-1T UD向超宽带记录；(b)仿真的UD向短周期记录；(c) DS-4A UD向短周期实测记录. Fig.3 Real and simulated earthquake records (a) JCZ-1T UD ultra broadband record; (b) Simulated SP record; (c) DS-4A UD real SP record. |
实践证明使用井下地震仪可以有效降低地震记录的环境噪声水平从而提高信噪比，但其水平向拾震器的方位不容易确定.如果在井底预先安装带锁定装置的底座，可使用井下陀螺仪定向，但井下陀螺仪很昂贵而且易损坏，因此很多井下地震仪的安装没做陀螺定向.一种折中的方法是同时在井口附近安装一台地面地震仪用相关分析的方法间接定向.而地表地震计有时也需要相对定向.吕永清等(2007)对相关分析的间接定向方法做了详细介绍.英国Güralp公司在其数采控制软件Scream(4.2以上版本)中集成有该功能，美国Geotech公司也有相应功能的SMARTAngle软件，这两款软件有一个共同点是在井下地震仪与地面地震仪中选择其一的两个水平向数据进行坐标旋转与另一个的NS向做相关分析给出结果.由于实际水平向地震计的正交性可能存在微小误差、更重要的是记录数据中存在有噪声，对两个水平向数据均做相关分析所得计算角度可能不尽相同，吕永清等(2007)据此建议使用两个水平向计算角度的均值.
吕永清等介绍的相关分析间接定向方法与Holcomb(2002)的方法有区别，区别不在算法而仅在于结果表示方法.由于Holcomb使用了极坐标图示法显示其计算结果，负的相关值与正的相关值重叠，使得Holcomb没能从图示结果分辨出地震计的正确方位.而吕永清等使用了直角坐标图示法，可以正确区分出正与负的相关系数角度关系.
吕永清等的文章中没有对分析的样本数据长度做说明，而Holcomb对井下观测数据分析的样本数据长度仅为256点.吕永清等与Holcomb的文章均说明记录数据中地脉动以外的随机噪声对进行方位角相关分析非常不利，为此作者选择较大的样本数据长度(例如100 sps采样，1小时的数据)以减少随机噪声带来的影响，这个选择令本文的计算结果有较好的稳定性和一致性.
当应用于相关分析间接定向的参考地震计与待测地震计具有不同频带时，吕永清等(2007)也只是建议选用两套仪器响应都平坦的部分进行同频带数字滤波，但根据作者的实践经验，在此情况下只做同频带数字滤波不一定能取得好的结果.为此作者尝试了使用本文第2节的递归滤波方法先将两套仪器中周期较长者的记录仿真成周期较短者的记录，然后再做同频带滤波.对几种不同频带的地震计进行多次实验都取得了较好的效果，下面以作者所做两个有仪器测量了地震计方位角的实验说明情况；两个实验均是用一台60 s的宽带地震计与一台2 s的短周期地震计做相关分析，但宽带地震计为不同生产厂家的产品；计算分别用有仿真和无仿真的方法对相同的记录数据进行；考虑到短周期地震仪的工作频带范围，作者选择滤波频带为0.5~5.0 Hz.用于测量地震计方位角的仪器是北京耐威时代科技有限公司的NV-NF301寻北仪，该仪器内置动力调谐陀螺，方位角测量不受环境磁场影响，角度测量范围0°~360°，测量精度0.3°，分辨率0.01°；每个角度用寻北仪重复测多次，取平均值以减少测量误差.
第一个实验用了表1序号4的宽带和短周期地震计组合，同时接在一台港震公司的EDAS-24GN6数采上，数采设置采样率为100sps，最小相位滤波器，±10 V输入量程.两台地震计安装在同一个基墩上，水平间距50 cm，寻北仪测得方位角分别 T34294为317.7°,09178为350.4°，角度差为32.7°. 对于无仿真方法，计算地震计EW向旋转34.0°时达到最高相关系数0.9761，NS向旋转33.5°时达到最高相关系数0.9744，平均值与寻北仪测量结果相差1.05°.图4a是两个水平向地震计以4°步长全方位旋转计算的相关系数曲线，图4b是在4°步长计算最佳角度基础上再以0.1°步长细化旋转计算的相关系数曲线，图4c、d、e、f是两台地震计两个水平向记录数据以最佳旋转角度旋转前后的同频带滤波结果局部波形图.对于有仿真方法，计算地震计EW向旋转33.1°时达到最高相关系数0.9993，NS向旋转32.7°时达到最高相关系数0.9989，平均值与寻北仪测量结果相差0.2°.图5是与图4对应的计算结果.显然经仿真计算的相关系数更高、旋转滤波后的波形图符合得更好而旋转角度更接近寻北仪测量值，结果优于未仿真计算的结果.
 | 图4 T34294对09178无仿真计算结果 (a)全方位相关系数图；(b)最佳旋转角度附近相关系数；(c)旋转前EW向数据；(d)旋转前NS向数据；(e)旋转后EW向数据；(f)旋转后NS向数据. Fig.4 Correlation results of T34294 vs 09178 without simulation (a) Raw correlation coefficients; (b) Fine correlation coefficients; (c) Original EW data;(d) Original NS data; (e) Rotated EW data; (f) Rotated NS data. |
 | 图5 同图4，但为T34294对09178仿真计算结果 Fig.5 Same as Fig.4, but for correlation results of T34294 vs 09178 with simulation |
第二个实验用了表1序号5的宽带和短周期地震计组合，同时接在一台港震公司的EDAS-24IP6数采上，数采参数设置与第一个实验一致.两台地震计也安装在同一个基墩上，水平间距50 cm，寻北仪 测得方位角分别08200为251.4°,11M46为206.3°， 角度差为45.1°.本次实验连续记录了27 h，对每小时的记录数据均进行了计算，角度差及相关系数结果见表2，图6、7是第一小时数据与图4、5对应的计算结果.全部计算均显示经仿真的相关系数高于未仿真的相关系数.夹角计算结果则显示未仿真的 计算离散性较大，从44.7°到50.6°都有，平均值46.9°； 经仿真计算的夹角则离散性较小，仅分布在45.4°到46.3°之间，平均值45.8°，比未仿真计算结果更接近寻北仪测量值.
|
|
表2 基墩实验二相关分析结果 Table 2 Correlation results of second experiment on pier |
 | 图6 同图4，但为08200对11M46无仿真计算结果 Fig.6 Same as Fig.4, but for correlation results of 08200 vs 11M46 without simulation |
 | 图7 同图4，但为08200对11M46仿真计算结果 Fig.7 Same as Fig.4, but for correlation results of 08200 vs 11M46 with simulation |
综合以上实验分析结果，在作者的实验条件下，两台地震计频带不相同，经仿真计算的结果显得更可靠、误差更小.
作者还在一个安装了英国Güralp公司井下观测系统的台站做实验，与系统安装时厂家给出的方位角做比较.该台站为广东区域地震台网的石榴岗台，地震计型号、参数见表1序号6、7的组合，地震计安装深度约220 m，安装时用地面地震计比测方法求出井下地震计方位角并设入CMG-DM24数采做了旋转矫正，地面地震计用罗盘定向，即系统输出 的NS对应地磁北南向.该台站在距离井口约100 m 外地面以下约15 m深的地下室中有浇注于基岩上的摆墩.
作者先是在距井口约1.5 m处地面安放了表1序号6组合中的BBVS-60地震计，用寻北仪定向，并用EDAS-24IP3数采接GPS天线对时记录；井下与地面记录数据采样率均为100sps.取了每段半小时共5段记录数据进行相关计算，对井下CMG-3TB的数据分别采取仿真和无仿真处理，根据广东地区约为2~3 s的地脉动周期选择0.2~0.6 Hz的带通滤波器进行同频带滤波，计算结果角度差平均值无仿真时-0.43°，有仿真时-0.17°(详情见表3，角度以井下CMG-3TB摆对地理坐标的偏离计算，顺时针方向为正).本实验有无仿真处理角度计算结果较接近，但相关系数是有仿真处理的结果稍高.
|
|
表3 井下仪器地脉动记录相关分析结果 Table 3 Microseism correlation results of downhole records |
作者之后在地下室中的摆墩上安放了一台同是Güralp产品的CMG-40T地震计(见表1序号7组合)，用寻北仪定向，用EDAS-24GN3数采接GPS天线对时记录，数据采样率同为100sps.实验系统连续记录了近1个月，与井下系统同时记录到了一些发生在不同方位的远震，作者挑选了6个记录较 清晰的事件数据进行处理(见表4)，地震参数取自中国地震台网中心发布的地震目录.对井下CMG-3TB的数据也分别采取仿真和无仿真处理，滤波频带则根据事件数据频谱分2种情况：云南地震震中距较小、频率成分稍高，采用10~2 s的通带，其他5个地震采用15~5 s的通带.计算结果角度差平均值无仿真时2.23°，有仿真时1.33°，详情见表5.本实验无仿真处理结果明显离散性大，相关系数更是远低于有仿真处理的结果.
|
|
表4 地震事件信息 Table 4 Earthquake events information |
|
|
表5 井下仪器事件记录相关分析结果 Table 5 Earthquake correlation results of downhole records |
根据数值实验及实际记录数据的仿真情况，对本文介绍的递归滤波方法的应用条件与使用效果归纳出以下几点：
(1)该方法有一定的适用条件，即要求记录数据是充分限带的或者忽略高频响应带来的误差是可以接受的；
(2)由于递归滤波系数比双线性变换方法得到的系数少，使用起来会更方便；
(3)该方法对速度记录反褶积得到的数据对应地动加速度的时间变化率；
(4)该方法在时域中处理，须假设少数几个初始值，但这较易从实际观测环境做出合理的假设，未发现在频率域方法中对长周期频谱修正时易出现的情况：不做修正或修正不合适会带来明显的长周期畸变李鸿吉，1992a,1992b).
综合来说，在充分注意适用条件及处理结果物理意义的情况下，应用该方法很容易取得令人满意的反褶积与仿真结果.
致 谢 作者的同事林伟高级工程师及劳谦、陈建涛、王国望、黄勇、贾雨等人参加了在石榴岗台井边地面及地下室的实验工作；中国地震局地球物理研究所杨大克研究员及两位匿名审稿专家对本文提出了一些宝贵、有益的建议，给本文增色许多；作者在此一并表示感谢！| [1] | Aki K, Richards P G. 1980. Quantitative Seismology Theory and Methods. Vol. II. San Francisco: W H Freeman and Company: 584-585. |
| [2] | Berger J, D C Agnew, R L Parker, W E Farrell. 1979. Seismic System Calibration: 2. Cross-spectral Calibration Using Random Binary Signals. Bull. Seism. Soc. Am. 69(1) : 271-288. |
| [3] | Holcomb L G. 2002. Experiments in Seismometer Azimuth Determination by Comparing the Sensor Signal Outputs with the Signal Output of an Oriented Sensor. U. S. Geol. Surv. Open-File Rept. 02-183. |
| [4] | Jin X, Ma Q, Li S Y. 2003. Comparison of four numerical methods for calculating seismic response of SDOF system. Earthquake Engineering and Engineering Vibration (in Chinese), 23(1): 18-30. |
| [5] | Jin X, Ma Q, Li S Y. 2004. Real-time simulation of ground velocity using digital accelerograph record. Earthquake Engineering and Engineering Vibration (in Chinese), 24(1): 49-54. |
| [6] | Li H J. 1992a. Examination of broadband seismic records recovery of CDSN in the frequency domain with the principle of modern control theory. Acta Seismologica Sinica (in Chinese), 14(1): 90-99. |
| [7] | Li H J. 1992b. Determination of the ground displacement using FFT and modern control engineering. Chinese J. Geophys. (in Chinese), 35 (1): 37-43. |
| [8] | Liu R F, Chen P S, Dang J P, et al. 1997. The application of imitation for broad band digital seismic recording. Seismological and Geomagnetic Observation and Research (in Chinese), 18(3): 7-12. |
| [9] | Lü Y Q, Cai Y X, Cheng J L. 2007. Orientation for seismometer with coherence analysing method. Journal of Geodesy and Geodynamics (in Chinese), 27(4): 124-127. |
| [10] | Oppenheim A V, Schafer R W. 1980. Digital Signal Processing (in Chinese). Beijing: Science Press, 146-150. |
| [11] |
Wang G F. 1986. Calibration and Checking of Seismograph Systems (in Chinese). Beijing: Seismological Press, 255-271. |
| [12] | Wang H T, Chen Y, Zhuang C T. 2006. A method on construction of real time seismic waveform emulation filter with bilinear transformation. Seismological and Geomagnetic Observation and Research (in Chinese), 27(1): 68-73. |
| [13] | 金星 ,马强,李山有. 2003. 四种计算地震反应数值方法的比较研究. 地震工程与工程振动, 23(1): 18-30. |
| [14] | 金星,马强, 李山有. 2004. 利用数字强震仪记录实时仿真地动速度. 地震工程与工程振动, 24(1): 49-54. |
| [15] | 李鸿吉.1992a. 用现代控制论原理检验中国数字地震台网宽频带地震记录的频率域复原. 地震学报, 14(1): 90-99. |
| [16] | 李鸿吉.1992b. 用FFT和现代控制论方法恢复地面位移. 地球物理学报, 35(1): 37-43. |
| [17] | 刘瑞丰,陈培善, 党京平等. 1997. 宽频带数字地震记录仿真的应用. 地震地磁观测与研究, 18(3): 7-12. |
| [18] | 吕永清,蔡亚先, 程骏玲. 2007. 确定地震计安装方位的相干性分析方法. 大地测量与地球动力学, 27(4): 124-127. |
| [19] | 奥本海姆 A V, 谢弗 R W. 1980. 数字信号处理. 北京: 科学出版社, 146-150. |
| [20] | 王广福.1986. 地震观测系统的标定与检查. 北京: 地震出版社, 255-271. |
| [21] | 王洪体,陈阳, 庄灿涛. 2006. 使用双线性变换构造实时地震波形仿真滤波器的方法研究. 地震地磁观测与研究, 27(1): 68-73. |