为检验地震计是否处于正常状态，需要自动定时地给地震计发出一组阶跃信号，通过分析地震计相应的原始输出信号确定仪器的传递函数，并判断仪器是否处于正常状况^[1-2]，以便及时进行相应的校正和维修处理。《测震台网运行管理细则(2015修订版)》规定，短周期和宽频带仪器每月进行一次阶跃标定，甚宽频带和超宽频带仪器每半年进行1次阶跃标定。阶跃标定可检验仪器的运行状况，但也会影响仪器的正常观测^[3-4]。若标定过程中突发地震，则无法利用记录的数据进行各种分析。因此，如何有效去除阶跃标定信号，对恢复地震信号及后续分析处理具有十分重要的意义。目前关于从原始信号中去除阶跃标定信号方面的研究较少，本文以湖北地震台网中麻城地震台记录的数据为例，提出一种去除地震计阶跃标定信号的方法，将理论阶跃标定信号与多项式拟合方法相结合，从阶跃标定信号中提取真实的地震信号，并利用PSD检验对本文方法进行验证。

1 方法概述

地震计的传递函数T(s)可近似表示为^[5]：

$ T(s) = \frac{{Y(s)}}{{\dot X(s)}} = \frac{{{S_v}{s^2}}}{{{s^2} + 2\xi \omega s + {\omega ^2}}} $

(1)

式中，s为复频率，即s=jω，Y(s)为地震计输出的拉普拉斯变换，(s)为地震动速度信号的拉普拉斯变换，S_v为地震计的速度灵敏度，ω、ξ分别为低频端截止频率和阻尼系数。

阶跃响应的时域响应公式为^[6]：

$ y(t) = \frac{{{e^{ - \xi \omega t}}}}{{\sqrt {1 - {\xi ^2}} \omega }}\sin (\sqrt {1 - {\xi ^2}} \omega t) $

(2)

甚宽频带地震计的低频端周期T=120 s、阻尼系数ξ=0.707，360 s后幅度变化未达到最大值的3/1 000 000，故阶跃信号宽度取T的3倍时长可满足要求，实际标定时间取1 200 s。

由傅里叶级数可知，阶跃信号由不同正弦波合成。例如宽度为1 200 s、幅度为1的阶跃信号的基频ω₀=2×π/2 400，阶跃信号可分解为^[7]：

$ P(t) = \frac{4}{{\rm{ \mathsf{ π} }}}\sum\limits_{n = 1}^\infty {\sin (n{\omega _0}t} )/n $

(3)

式中，P(t)为阶跃信号，n为正奇数。由式(3)可知，阶跃标定是以ω₀为基频、3次及以上奇数次谐波频率之和。阶跃标定信号为全频带信号，因此无法采用常规滤波方法进行去除。

式(1)为简化后的地震计传递函数，实际的传递函数则更为复杂。用简化函数计算得到的地震计周期、阻尼与理论值略有差异，因此不能简单地只去除理论阶跃标定信号，还需进行后续的相关处理和计算，真实有效地还原地震信号。本文研究思路为：1)根据《测震台网运行管理细则(2015年修订版)》规定，阶跃标定幅度信噪比≥40 dB，即响应曲线的幅值为地脉动信号幅值的100倍以上，这样才能推算出比较准确可靠的周期和阻尼系数。本文选取安静时段的标定信号，根据计算的结果得到理论的时域响应值，然后从原始信号中去除此理论响应值。2)通过合适的方法去除1)中未完全去除的干扰，使得标定时段及标定前后各1 h的功率谱大致相同。本文采用2次多项式拟合的方法，第1次多项式为18阶，第2次则只需3阶。3)利用上述2个步骤得出的经验公式对地震时段的标定信号进行处理。

本文以麻城台CTS-1E甚宽频带地震计垂直向数据为例。CTS-1E甚宽频带地震计为CTS-1甚宽频带地震计的升级产品，具有甚宽频带、大动态范围、高灵敏度的特点，观测频带为50 Hz~120 s，灵敏度为2 000 V ·s/m^[8]。该地震计配备的数据采集器为EDAS-IP24，量程为20 V，转换因子为3 192 nV/count，采样率为100 sps。

2 方法实现

选取麻城台CTS-1E甚宽频带地震计安静时的标定数据及标定前后各1 h的数据(图 1(a))，图 1(b)为完整标定数据。图 1(b)中标定数据的单向最大幅度为1 049 365 count，360 s后理论阶跃幅度值不足3 count。图 1(a)中噪声信号约为标定最大幅度的1/5 000，满足阶跃标定要求。

图 2为3段数据所对应的PSD，图中黑色曲线分别为全球高噪声模型(NHNM)和低噪声模型(NLNM)，绿色、蓝色和红色曲线分别为阶跃标定信号、标定前1 h时信号和标定后1 h信号相对1 m²/s³的1/3倍频程PSD。从图 1(a)可以看出，标定前1 h信号由于在2 250 s左右存在干扰，其PSD在0.1 Hz以下的低频端比标定后1 h信号的PSD大，最大差值约为18 dB(图 2)。

通过EDAS-ICHECK软件计算上述标定数据对应的仪器参数，得到T=119.983 8 s，阻尼D=0.713 512。在去除理论阶跃曲线时，需要先将原始标定信号输出的初始值(在本文中取标定启动前5 s的平均值，正向为110 count，反向为35 count)赋给理论阶跃起始点，再将实际最高值(在本文中取最大和最小时平坦处±0.2 s的平均值，正向为1 049 475 count，反向为-1 049 525 count)作为理论阶跃的最大幅度。

图 3(a)为原始信号与理论阶跃信号之差，上道为正向标定后600 s信号，下道为反向标定后600 s信号。从大致趋势来看，2次结果的前300 s存在一长期波动信号，且方向相反。前3个点幅度变化较大，正向为1 352 count，反向为1 254 count，此结果可能为将理论阶跃曲线按简化的二阶高通模型计算而忽略50 Hz低通所造成，但考虑低通后的信号仍与原结果基本一致。因此，本文采用简化的二阶高通模型进行数据处理，由于前3个点变化较大故用4~6点来替代。对标定前300 s的长周期波动信号进行多项式拟合得到拟合曲线，然后将该段信号与拟合曲线作差。实验结果表明，至少需要13阶才能基本消除长周期信号，18阶之后基本无变化。图 3(b)为在实际处理中采用18阶进行处理的结果。

从图 1(a)可以看出，标定数据在标定开始和结束后360 s，仪器信号在一段时间内仍继续漂移。图 4上道为该段数据去除理论标定信号和第1次多项式拟合信号的结果，数据处理结果表明，该段数据还需进行第2次多项式拟合。图 4下道为对标定开始300 s后的数据进行3阶多项式拟合的结果。图 5为2次去除拟合信号后的PSD与标定前后各1 h信号的PSD，去除第2次拟合后信号的PSD与相对安静时标定信号的PSD已基本一致。

3 地震信号去除阶跃标定

以麻城台CTS-1E甚宽频带地震计记录的一次远震事件为例，探讨从原始信号中去除阶跃标定信号的方法。选取2015-04-17 02 :07 :44希腊克里特岛(35.1°N，26.8°E)6.1级地震为例，震源深度20 km。地震计在02 :26 :15.47开始阶跃标定，宽度为1 200 s，图 6(a)为CTS-1E记录的此次地震的垂向分量信号。从图 6可以清晰地辨别直达P波震相到时，后续震相到时由于受阶跃标定信号的影响，振幅和周期发生突变，无法判断出后续震相具体到时。利用本文方法将原始信号与理论标定信号作差，然后分别去除18阶多项式拟合和3阶多项式拟合结果，得到真实的地震信号(图 6(b))。图 6(b)不仅可以清晰地辨别P波到时，后续震相到时也被清晰恢复。

4 结语

去除地震计记录中的标定信号对恢复地震信号、减小信号处理误差及提升计算精度方面具有重要意义。本文以麻城台安静时的记录数据为例，得到去除阶跃标定信号的步骤为：

1) 利用EDAS-ICHECK软件计算出阶跃标定信号的周期和阻尼系数，以此得到理论阶跃标定信号并放大至与实际信号相同的初始值及幅值；

2) 将数采记录信号f₁与理论阶跃信号f₀作差，得到去除阶跃标定后的信号f₂=f₁－f₀；

3) 对去除阶跃标定信号的前300 s进行18阶多项式拟合得到拟合信号p₁，并将f₂与p₁作差得到信号f₃= f₂－p₁；

4) 对阶跃标定300 s后的信号进行3阶多项式拟合得到拟合信号p₂，并将信号f₃与拟合信号p₂作差，得到最终的真实信号f=f₃－p₂。

分析表明，去除阶跃标定后的地脉动信号与标定前后的地脉动信号的PSD较为接近，该方法可用于去除地震计原始记录信号中的阶跃标定信号，并可有效恢复远震信号的不同震相。

本文对CTS-1E地震计信号去除阶跃标定的方法进行分析，其他地震仪的处理方法类似，不同仪器的多项式拟合可能稍有差别，应结合观测数据进行具体分析。