0 引言

作为国家地震前兆台网的重要观测手段之一，井水位观测一直在地震前兆异常的提取、地震预测研究、震情跟踪中发挥着重要作用(Roeloffs，1988；Cicerone et al，2009)。大地震发生之前，在井水位观测中经常记录到明显的异常现象，如2008年汶川8.0级(车用太等，2008；邱永平等，2009；程万正等，2013；邱永平，2013；王小娟等，2014；王军等，2018；晏锐等，2018)、2010年玉树7.1级(孔令昌等，2012)等地震。地下水位观测数据具有较宽的记录频带，不同频段的数据具有不同的应用价值。长周期信号可以用来研究震前的趋势性变化(陆明勇等，2010)以及大地震造成的地球自由振荡(Yan et al，2014、2016)等；高频成分可以用来研究同震响应(王永刚等，2016；向阳等，2017)。

井水位观测不同程度地受到观测站点的场地、自然环境、人为活动等因素的制约与干扰，不同尺度的干扰源产生的干扰信号的频段存在差异(李继业等，2015；孔庆敏等，2018)。因此，将记录分为不同的频带，然后进行进一步的研究更有意义。

频段分离的方法有多种，其中，滤波的方法最为常用。滤波利用有效信号和噪声不同的频率特性将其分离，从中提取出有效成分。但是，滤波会使信号的波形发生畸变，进而对后续研究造成影响。另外，数字信号处理方法大多直接来源于信号分析的专业工具，它们具有普适性，但是未必是处理地震前兆记录的最佳选择，因此在信号处理过程中需要考虑信号自身的特征。

本文选取张道口-1井和新10井的井水位观测数据，考虑到观测数据既包含高频的同震响应以及阶变信号、具有稳定周期的固体潮信号，又包含低频的趋势性信息的特征，故尝试使用集合经验模态分解方法(ensemble empirical mode decomposition，EEMD)将信号分为不同的频段，以期为后续的噪声压制、干扰信号去除以及前兆异常事件的识别和研究提供高质量的基础数据，并为前兆观测资料在地震预测中的应用提供科学依据。

1 方法简介

希尔伯特-黄变换(Hilbert-Huang Transform)是分析非线性、非平稳信号的有力工具(Huang et al，1998)，它包括经验模态分解(EMD)和希尔伯特变换2部分。经验模态分解基于信号自身特征将其分解为物理意义明确的几个本征模态函数(IMF)和1个低频趋势，具有自适应的特点，在地球科学、大气科学以及生物科学等诸多领域都有较好的应用效果(Huang et al，2008；Song et al，2012)。对于给定的信号X(t)，其经验模态分解可以表示为

$ X(t)=\sum\limits_{i=1}^{n} c_{i}+r_{n} $

(1)

其中，c_i为IMF；r_n为趋势项。

IMF是具有变化周期和振幅的振荡函数，满足以下特征：①对于整个IMF函数，极值点的个数与跨越零点的个数之差不超过1；②IMF的任意一点，极大值和极小值包络线的平均值为0(现实中一般不会严格等于0，但是满足给定的误差范围)。

对于给定的数字信号序列X_i(i=0，1，…，N-1)，EMD的计算流程为：

(1) 首先找出数据的极大值、极小值点，使用3次样条函数拟合出上下包络线。

(2) 计算上下包络线的均值函数m₁₀。从原信号中减去m₁₀得到残差h₁₀，为1次拟合

$ X(t)-m_{10}=h_{10} $

(2)

理想情况下，h₁₀即为1个IMF；实际中，该过程要重复多次以满足IMF的特征要求

$ h_{10}-m_{11}=h_{11}, h_{11}-m_{12}=h_{12}, h_{1(k-1)}-m_{1 k}=h_{1 k} $

(3)

(3) 当h_1k满足IMF的特征要求时，第(2)步的循环停止。记h_1k为C₁，得到第1个IMF函数。从原信号中减去c₁记为r₁，重复上述过程，依次找出c₂、c₃、…和r_n。

为了解决信号的局部微小扰动给EMD计算带来不稳定的问题以及EMD计算存在的混频现象，Wu等(2009)改进了EMD算法，通过多次叠加白噪声以增加计算结果的稳定性，即为集合经验模式分解(EEMD)。

2 张道口-1井水位数据处理 2.1 井水位基本情况

张道口-1井深约1150m，位于沧县隆起上双窑凸起与白塘口凹陷交界的白塘口西断裂带上升盘(王建国等，2007)。选取2017年4~10月的水位分钟值数据，波形如图 1(a)所示，图 1(b)给出了傅里叶变换振幅谱。由图 1(b)可见，在12.00、12.44、24.00h时有3个峰值分别对应固体潮的半日波和全日波。

2.2 EEMD处理

将EEMD方法应用到上述水位观测数据中，得到18个IMF，其中，第1~10个如图 2所示，第11~18个主要包含长周期趋势性信息。由图 2可见，第1~4个IMF表征信号的高频信息至少包含2个高频事件；第5~7个IMF表征具有固定振动周期的信号；而第8~10个IMF则主要表征信号的长周期成分。我们将第1~4个IMF叠加在一起，视为信号的高频分量；将第8~18个IMF叠加在一起，视为信号的低频趋势(图 3)。由图 3可见，叠加之后，对于高频事件的识别和低频趋势的判定就变得较为容易。

2.3 水位稳定周期成分分析

IMF5~7的波形图、小时-天谱图和傅里叶振幅谱如图 4所示。由图 4可见，IMF5的小时-天谱图信号每天有2个峰值和谷值，说明信号具有约12h的周期；同时，峰谷值出现的时间随天数的增加而向后推移，说明信号的周期略大于12h；IMF5的傅里叶振幅谱也表明这段信号的主周期约为12h；这些特征与M2固体潮的周期特征非常吻合。不仅如此，IMF5的小时-天谱图在横向上还表现出14个彩色条带，说明这组信号存在周期为半个月的旋回。由图 4还可见，IMF7的小时-天谱图和傅里叶振幅谱同样表现出信号略大于24h的周期特征和周期为半个月的旋回特征；IMF6的傅里叶振幅谱信号兼有12、24h周期的成分，说明EEMD方法在分解的时候存在混叠现象。

IMF5~7信号都具有明显的固体潮变化特征，对固体潮的研究具有较重要的科研价值，例如利用固体潮观测数据可以探究地球深部结构(Lau et al，2017)等，其中关键的一步是将固体潮记录提取出来。有人曾尝试使用小波分析方法提取地下水位观测的潮汐变化(杨从杰等，2005)；然而，小波分析是基于给定的小波函数，不具备自适应性(Huang et al，1998；Huang et al，2008)。本文使用EEMD方法得到的IMF5具有清晰的M2波特征，IMF7具有全日波特征，可以为后续研究提供基础数据。虽然IMF6混叠了半日波和全日波的信号，可能不适合直接应用；但是，可以将IMF5~7的信号叠加起来，联合研究半日波、全日波信号特征。

2.4 水位高频成分分析

分析图 3(a)中的2个大振幅高频事件，它们的发生时刻分别为2017年7月18日7时35分和9月8日12时50分左右。以第2个事件为例，该波形显示出同震响应事件的特征。通过查询地震目录(http://www.ceic.ac.cn/)，对比理论到时发现，该事件与2017年9月8日12：49：15发生在墨西哥沿岸近海的8.2级地震在时间上非常吻合，可能是该地震的同震响应。使用同样的方法可以确认，第1个大振幅高频事件是2017年7月18日07：34：13科曼多尔群岛地区7.8级地震的同震响应。

3 新10井水位数据处理 3.1 井水位基本情况

新10井位于乌鲁木齐市南部(43.7°N，87.62°E)，海拔高度为1056m，构造上属于柳树沟-红雁池逆冲断裂。新10井自1980年9月成井以来，井水位清晰地记录到2017年8月8日九寨沟M_S7.0、8月9日精河M_S6.6地震的同震响应(向阳等，2017)。选择2016~2018年的记录用于研究，波形图如图 5 (a)所示。

3.2 基于EEMD的井水位同震分析

将EEMD方法应用到2016~2018年新10井的静水位分钟值观测记录中，共得到20个IMF。根据前文的分析，将IMF1~4叠加在一起得到信号的高频成分(图 5 (b))用于后续同震响应的研究；将IMF11~20叠加在一起得到信号的低频成分(图 5 (c))。图 5中的IMF5~7为具有固有周期的信号分量。通过查询地震目录(http://www.ceic.ac.cn/)得知，新10井周围500km范围内2016~2018年共发生了4次5级以上地震，地震事件的详细信息如表 1所示。

放大信号的高频成分，研究上述4次地震的同震响应，波形如图 6所示。由图 6可见，新10井记录到了除事件1之外其他事件的同震响应。事件2、3的同震响应表现为振荡型，振幅先增大后减小，最大幅度分别为2.36×10^-3、2.97×10^-3m。事件4的同震响应表现为下降型，对应原始记录的阶降，最大幅度为2.13×10^-3m。事件3、4的震中距相当；但是事件3的震级比事件4大0.9，这可能是造成二者同震响应幅度存在差异的原因；事件2与事件3间同震响应幅度的差别为震级差和震中距差共同作用的结果。分析地震事件之前的异常记录发现，IMF10的分量中，事件2~4之前都存在明显的大幅度振荡(图 5 (b))。

查询地震目录得知，2016~2018年全球共发生了40次7级以上地震，将这些事件的发震时刻标在用EEMD方法得到的高频分量信号上发现(图 7)，新10井的高频信号记录到了大部分地震事件的同震响应，它们在记录上表现为一系列的脉冲信号。放大发震时刻附近的波形，可以进一步确认新10井记录到了40次地震中的31次事件的同震响应。有些脉冲信号没有对应的地震事件，可能来自记录中的台阶变化或发生在周边的震级较小的地震。

新10井记录的地震同震响应的最大幅度随震中距和震级的分布情况如图 8所示。由图 8可见，总体上最大幅度随震中距增大而减小，随震级增大而增大。我们构建最大振幅的对数随震中距和振幅线性变化的关系式，使用上述31次事件结果拟合得到的最佳公式为

$ \lg A=-2.31 \times 10^{-3} D+0.75 M-7.70 $

(4)

其中，A为最大振幅；D为震中距，单位为°；M为震级。图 8还给出了震级为7.0~8.2的理论计算结果。由图 8可见，它们在整体上与观测结果吻合。图 9给出了理论值与观测值误差的绝对值。由图 9可见，除了3次7.8级地震的误差较大以外，其他的误差约为0.4，表明式(4)可以为新10井静水位观测的同震响应的定量研究提供帮助。

4 讨论与结论

本文首先使用EEMD方法将张道口-1井地下水位观测数据分为18个IMF；通过对它们进行线性组合得到高频分量(IMF1~4)、具有固定周期的分量(IMF5~7)和低频趋势分量(IMF8~18)。划分之后，从高频分量识别包含同震响应在内的孤立事件变得较为容易；同时，识别出来的同震响应事件与测震数据吻合很好，这为定量地研究地下水位观测机理提供了思路。具有稳定周期的分量表现出固体潮半日波和全日波的信号特征，为进一步研究提供了关键的数据。而低频分量则可以为地下水位观测长趋势性变化研究提供条件和基础。在此基础上，将EEMD方法应用到2016~2018年新10井静水位观测分钟值数据，从高频信号中对比研究了31次7级以上地震的同震响应振幅，得到了同震响应随距离和震级的变化关系为lgA=-2.31×10^-3D+0.75M-7.70。

致谢: 地下流体观测数据和地震数据来自中国地震台网中心；研究中用到的EEMD程序包来自台湾“中央大学”数据分析方法研究中心(http://rcada.ncu.edu.tw/research1.html)。两位审稿人的宝贵意见大大提升了本文的质量。在此一并表示感谢。