大地震会激发持续时间很长的地球自由振荡.通过比较地球自由振荡的实际观测周期和各种模型的预测周期，获得了另外一种独立于地震体波的研究地球物理参数随深度变化的方法，为改进地球密度模型和研究地球内部构造提供了重要手段(傅承义等，1985).超导重力仪(Superconducting Gravimeter, SG)作为目前观测精度最高的重力仪，具有高度的灵敏度和稳定性，以及很宽的动态线性测量范围和极低的噪声水平，能够检测到周期从几秒到几年乃至更长的地球动力学效应(Crossley et al., 1999; Goodkind, 1999; Hinderer and Crossley, 2000, 2004)，对地球自由振荡信号的观测同样具备检测能力(Rosat et al., 2004, 2005; 雷湘鄂等，2007).利用超导重力仪观测资料检测地球自由振荡是“全球地球动力学计划”(GGP)的主要研究内容之一，也是相关地学研究的国际前沿(Crossley et al., 1999; 孙和平等，2006; Crossley and Hinderer, 2008a, 2008b).相关研究表明：由于SG具有很低的噪声水平，利用SG数据检测地球自由振荡会得到比其他观测仪器更为精确的结果(Banka and Crossley, 1999; Rosat and Hinderer, 2011).尤其是在低于1 mHz的频率范围，SG观测具有很高的信噪比(Van Camp, 1999; Rosat et al., 2003a, 2004；Widmer-Schnidrig, 2003)，可以对由于地球自转和椭率引起的低频简正模式的频谱分裂现象进行高分辨率分析(Rosat et al., 2003b, 2005; Shen and Luan, 2015).精确测定低于1 mHz的地球自由振荡简正模式的分裂频率是在不与任何弹性系数发生联系的情况下改善一维密度模型的一种方法(Widmer-Schnidrig, 2003)，极低频率的简正模式分裂还对地幔和地核的三维密度结构有高灵敏度，因此对它们的观测还可提供地球横向密度结构的约束(Ritzwoller and Lavely, 1995).研究结果表明SG将在构制频率低于1 mHz的长周期地震图和研究地球深部构造方面发挥重要作用(孙和平等，2006).

全球超导重力仪的噪声水平分析给出的相关结果表明，在小于1 mHz的频率范围内，超导重力仪的噪声水平开始低于宽频带地震仪和弹簧重力仪(Banka and Crossley, 1999)，而在该频率范围内，气压变化对重力观测的影响成为继潮汐影响以外最重要的因素.相关研究指出台站局部气压对超导重力观测数据的影响主要集中在低于1 mHz的频段上(Zürn and Widmer, 1995; Freybourger et al., 1997; Rosat et al., 2003a, 2005)，因此在利用超导重力数据研究低频频段上的自由振荡信号时，必须精细地消除气压噪声，提高微弱地球物理信号的信噪比，才能更为有效地检测地球简正模的相关特征参数.

消除气压变化对重力观测的影响的方法可分为两大类：物理模型方法和试验分析方法.物理模型方法是根据Farrell(1972)提出的弹性地球的静态负荷的重力格林函数理论，利用大气密度模型、地球模型和全球性的气温和气压观测资料，从理论上分析计算大气负荷与重力变化的关系.物理模型方法的优点是可考虑全球范围内气压变化对重力的影响，但物理模型方法不能考虑局部天气变化的影响，对全球性大气热潮效应的改正效果也并不好(Boy et al., 1998).

由于局部气压变化对重力观测的影响最大，Warburton和Goodkind(1977)建议用重力观测站的同步气压观测数据对重力观测中的气压噪声进行改正.利用线性最小二乘法，对气压变化与重力残差序列作回归分析，得到局部气压导纳值.这样得到的导纳值大约为-3 nm·s^-2·hPa^-1，但不同台站的导纳值有所不同，这种气压导纳值只是重力观测在整个频段上对气压变化的平均响应.然而当局部气压变化时，气体的运动能量在不同频段的分布是不同的，相对于高频部分，大气的低频信号具有更高的能量和较大的空间活动范围，因此会造成气压变化对重力观测的影响具有一定频率依赖性(Crossley et al., 1995；胡小刚，2006；Hu et al., 2006).为了充分地消除气压对重力观测的影响，就需要人们在不同的频段上分别估算气压导纳值.Warburton和Goodkind(1977)采用互谱分析方法(Cross-Spectral Analysis)在频率域给出了复导纳值的定义，最先提出在频域计算局部气压导纳值.Crossley等(1995)和Neumeyer(1995)利用频域上的回归方程估计了具有最小二乘意义的复气压导纳值.Hu等(2006)基于连续小波变换方法分析了局部气压对重力变化影响的时频依赖性.

本文将利用经验模态分解(Empirical Mode Decomposition, EMD)的改进方法EEMD(Ensemble EMD)，提出一种新的气压改正方法.该方法将低频自由振荡信号独立分解到一些窄带频段上，同时对台站同步气压观测进行相同的处理，并在相同频段上分别进行重力—气压变化的回归分析，得到频率依赖的气压导纳值，从而对重力残差序列进行精细的气压改正.基于该方法本文对大震后的低频频谱分裂现象进行了高分辨率分析，检测结果验证了该方法的有效性.

经验模态分解(Empirical Mode Decomposition, EMD)是希尔伯特—黄变换(Hilbert—Huang Transform，HHT)方法的重要组成部分(Huang et al., 1998)，HHT方法是一种处理非线性非稳态数据的新方法，已广泛应用到地球物理的各个领域(Huang and Wu, 2008; Gu et al., 2014).

与传统的信号分解方法相比，经验模态分解(EMD)采用的分解方法是直观的、自适应的，不包括任何先验信息，以数据本身各种内在的特征时间尺度(Characteristic Time Scale)特性为基础直接在时间域上进行.EMD分解并不涉及到任何频域的内容，完全由自身性质出发，从而克服了傅氏变换在处理非线性、非平稳信号的不足，同时回避了小波方法需要人为选取小波基的局限性.该分解方法建立在一个简单的假设基础之上，那就是任何信号，包括非稳态信号，在任意给定的时间段内，都是由一些简单的固有振动模式组成，它们相互叠加在一起从而构成信号本身.EMD分解的目的就是希望从信号中提取出这些固有的振动模式，我们将这些固有振动模式命名为“本征模态函数”(Intrinsic Mode Function，IMF)，并且要求它们满足下面2个基本条件：

(1) 整个时间序列的极值点的数量(即极大值和极小值数目之和)与过零点的数量相等，或最多相差1个.

(2) 在时间序列的任意一点上，由极大值插值确定的包络线与由极小值插值确定的包络线的均值始终为零.

根据上面给出的本征模态函数(IMF)的定义，通过EMD分解中所谓的一个“筛选” (Sifting Process)过程，就可以将任意信号分解成一组IMF和一个残余量，对该方法的详细介绍请参考文献Huang等(1998)，本文在这里不再赘述.

Wu和Huang(2009)基于噪声辅助数据分析(NADA)理论(通过增加白噪声避免模态混叠，有助于提取数据中的真实信号)，对EMD方法作出了进一步改进，提出EEMD(Ensemble EMD)方法.EEMD方法在EMD的分解过程基础上增加了如下过程：

(1) 为目标数据增加了一个白噪声序列；

(2) 将添加了加白噪声的数据利用EMD分解成IMF分量；

(3) 反复重复步骤(1)和(2)，但是每次采用不同的白噪声序列；

(4) 计算相应的IMF分量的集合平均值(Ensemble Means)，并将其作为最终的分解结果.

EEMD方法成功地解决了EMD分解中的模态混叠问题(Mode Mixing)，使分解得到的IMF具有更为集中的频率信息.该方法为下一步进行具有频率依赖特性的精细化气压改正提供了有效的数学工具.

当地震等干扰因素存在时，我们并不能通过直接对观测数据进行回归分析得到正确的气压导纳值.要正确评估气压对重力观测的影响，必须在没有其他噪声的基础之上才能完成.为了得到各个频段上气压变化对重力的不同影响，就要求我们选用平静时期的连续重力观测信号，并利用与其同步观测的气压数据，将它们分解到不同的频段，再通过线性回归分析得到正确的气压导纳值.

基于噪声水平分析结果(Rosat and Hinderer, 2011; 江颖等，2016)，首先选择每个台站连续数天的平静期观测，在去掉潮汐信号后，对重力残差和同步气压观测分别进行EEMD分解，将两种类型的观测信号分解到相同的窄带频段上，然后计算相应的气压导纳值.在对平静期观测去潮汐信号时，可以利用调和分析得到精确的潮汐参数，通过计算合成潮的方式完成潮汐改正.

频率依赖的大气导纳值可由线性回归方程解算得到，公式为

(1)

式中g₀(Δf)和P(Δf)是在子频段Δf中的重力变化和气压变化.利用最小二乘方法求解线性回归方程，计算|g_c(Δf)|²为最小值时的回归系数c(Δf)，即气压导纳值.通过这种方法计算的大气导纳值是频率依赖的，因而可更有效地消除气压变化对重力观测的影响(Crossley et al., 1995).

由于EEMD分解得到的IMF的个数与数据长度有直接关系(Wu and Huang, 2009)，这里我们选用为期10天的分钟值记录，结果证明这样的方式可将两类观测数据分解到合适的窄带频段上.图 1给出了法国Strasbourg台站在一段为期10天的平静期内的重力残差及其经EEMD分解得到的IMF结果.图 2是该台站同步观测的气压变化及其EEMD分解结果.结果显示：对经过潮汐改正后得到的重力残差和同步气压观测进行EEMD分解后，分别一共可以得到12个IMF和一个残余趋势项.比较两类数据的时间序列，可以发现重力残差和气压变化的原始信号相关性很大，说明在去掉潮汐的影响后，气压是导致重力变化的主要来源.为了精密确定地球自由振荡的相关参数，对气压影响进行改正是很有必要的.

图 1

Fig. 1

图 1 法国Strasbourg台站在平静期的重力残差及其EEMD分解结果 数据长度为10天，纵坐标单位均为nm·s-2. Fig. 1 Gravity residuals and its decomposing results based on EEMD during quiet period at Strasbourg in France Data length is 10 days, and the unit of vertical axis is nm·s-2.

图 2

Fig. 2

图 2 法国Strasbourg台站的同步气压观测及其EEMD分解结果 数据长度为10天，其中气压观测信号减去了自身的平均值，纵坐标单位均为百帕(hPa). Fig. 2 Simultaneous atmospheric pressure and its decomposing results based on EEMD during quiet period at Strasbourg in France Data length is 10 days, and the unit of vertical axis is nm·s-2.

对于经EEMD分解得到的长周期IMF(从IMF7到IMF12)，直接从时间域里就能看到重力残差和气压变化具有相似的形态，相关性(负相关)比较高，并且IMF的频率依照分解顺序依次降低.而处于相对高频段的IMF单从时间域很难比较它们的频率大小以及两类数据的相互关系.因此为了验证这两种类型的信号经过EEMD分解后是否被分解到相同的频段上，本文对重力残差和气压变化的前6个IMF进行傅氏变换，将它们转化到频率域内进行相互比较.图 3给出了气压变化和重力残差经过EEMD分解以后得到的前6个IMF的振幅谱之间的比较.从图中可以看到，一方面每种信号的IMF严格按照分解顺序从先到后频率依次降低，而且没有出现模态混叠的情况；另一方面两种信号对应的IMF集中在相同的频段上，除IMF-1稍有不同外，其他每组对应的IMF的振幅谱形态具有很高的相似性.

图 3

Fig. 3

图 3 同步气压和重力残差经过EEMD分解以后得到的前6个IMF的振幅谱 气压的振幅单位为帕斯卡(Pa)，重力残差的振幅单位为nm·s-2，所有振幅谱经过37个点的Parzen谱窗的平滑处理. Fig. 3 Amplitude spectra of first 6 IMFs of simultaneous atmospheric pressure and gravity residuals The unit of pressure′s amplitude is Pa and the gravity residual′s is nm·s-2. All amplitude spectra are smoothed by a 37-point Parzen spectral window.

图 3的结果表明：经过EEMD分解以后，重力残差和同步气压观测被分解到了相同的频段上，完全可以满足信号分解的目的.下面我们利用公式(1)，针对重力残差和气压观测每组对应的IMF，计算频率依赖的大气导纳值，结果列于表 1.另外为了分析气压变化对重力的影响程度，表 1同时也给出了每组对应的IMF的相关系数.

表 1 Table 1

表 1 重力观测和气压观测相互对应的IMF之间的导纳值和相关系数 Table 1 Admittances and correlation coefficients of the same IMF between gravity and atmospheric pressure 子频段 导纳值 相关系数 子频段 导纳值 相关系数 IMF-1 0.3312 0.0199 IMF-7 4.1479 0.9286 IMF-2 0.4165 0.0540 IMF-8 2.4778 0.9118 IMF-3 1.7939 0.5352 IMF-9 2.4952 0.9583 IMF-4 3.5496 0.8569 IMF-10 2.8154 0.9798 IMF-5 3.9893 0.9213 IMF-11 2.5760 0.9888 IMF-6 3.5674 0.9484 IMF-12 3.7962 1.0000 注：导纳值和相关系数均取了绝对值, 导纳值的单位为nm·s-2·hPa-1.

表 1 重力观测和气压观测相互对应的IMF之间的导纳值和相关系数 Table 1 Admittances and correlation coefficients of the same IMF between gravity and atmospheric pressure

表 1的结果显示，气压变化和重力残差的相关系数从IMF-1到IMF-12基本呈现出一种单调上升的趋势，也就是说频率越低，这两种信号的相关性越强；尤其从IMF-4开始，气压变化和重力的相关性全部保持在0.85以上，由图 3可以看到恰恰是以IMF-4为起点，两种信号的优势频率开始低于1 mHz，这再次说明了气压对重力的影响主要集中在1 mHz以内的频率段上.每组IMF计算得到的气压导纳值的大小都有所不同，但并没有明显的单调趋势.在高频段，从IMF-3到IMF-1，气压变化和重力的相关系数越来越低，而导纳值也越来越小.

为了验证基于EEMD分解得到的两种信号IMF间的相关性是否正确，本文利用信号相干分析方法直接从频率域比较气压变化和重力残差序列在各频率的相关性.信号x和y的相干函数可定义为(Kay, 1988)：

(2)

式中，S_xx(f)、S_yy(f)分别为信号x和y的自功率谱密度函数，S_xy(f)为互功率谱密度函数.相干函数直接从频域角度揭示了两种信号的相关关系，就统计学角度而言，相干函数值介于0.1~0.3为弱相关，0.3~0.5为中等相关，0.5~1.0为强相关.相干值越大，表明在相应频段重力观测受气压变化的影响越大.

利用公式(2)，对上述法国Strasbourg台站为期10天的重力残差序列和同步观测的气压变化做相干分析，结果由图 4给出.为了便于比较，图 4同时给出了两种信号处于不同频段的各IMF间的相关系数的绝对值.在不同频段做相关分析可以看作是相干分析在特定频率点的表述，从图 4中可以看到，相干分析的结果和各频段IMF间相关系数，在趋势上是一致的：在小于4 mHz范围内，两条曲线基本重合；特别在小于1.5 mHz的频段里，两种信号的相干值和相关系数的绝对值均大于0.8，说明在该频段重力残差受到了气压变化的强烈影响.通过比较两种相关分析方法的结果，验证了本文基于EEMD分解得到的两种信号的IMF间相关特性的正确性.

图 4

Fig. 4

图 4 法国Strasbourg台站为期10天的重力残差序列和同步观测的气压变化的相干分析结果和这两种信号处于不同频段的各IMF间的相关系数的绝对值的比较(为了更好的比较，相干分析结果经过101个点的帕曾谱窗的平滑处理) Fig. 4 Coherence analysis of gravity residuals and simultaneous atmospheric pressure at Strasbourg in France and comparison of absolute values of correlation coefficients for their corresponding IMFs located at different frequency bands (For better comparison, the results of coherence analysis are smoothed by a 101-point Parzen spectral window)

为了进一步分析不同频段气压变化对重力观测的不同影响，下面针对处于不同频段的IMF分别研究气压变化和重力残差的相互关系.由表 1可以看到处于低频段的IMF7到IMF12，重力和气压观测之间的相关性很强.从后文的结果可知长周期IMF(从IMF7到IMF12)的优势频率都集中在小于0.1 mHz的频段上，因此我们将重力和气压变化的长周期IMF分别相加并研究它们各自的和信号.图 5a和b给出了气压变化和重力残差的原始信号与它们各自的IMF7到IMF12的和信号的比较.从图中我们可以看到，IMF7到IMF12的和信号与各自本身的原始信号之间在大的变化趋势上很相像，其差别主要来自两方面：一个是IMF7到IMF12的和信号中没有包含趋势项，另外就是将小振幅的高频信息排除在和信号之外.因此气压和重力信号在趋势上的相关性主要来自于0.1 mHz以下的贡献.图 5c给出了两类和信号的归一化振幅谱，两条谱峰基本重叠在了一起，同样说明了气压变化和重力在0.1 mHz以下存在着极高的相关性.另外我们计算了两类和信号对应的气压导纳值和相关系数，分别为-3.094 nm·s^-2·hPa^-1和0.991，这些数值反映了气压变化和重力在低于0.1 mHz的频段上的相互关系.

图 5

Fig. 5

图 5 (a) 和(b)中黑色实线分别为气压观测和重力残差，蓝色实线和红色实线分别为气压变化和重力残差的IMF7到IMF12的和信号；(c)蓝色实线和红色实线分别为气压和重力残差的IMF-7到IMF-12和信号的归一化振幅谱 Fig. 5 (a) Atmospheric pressure during 10 days. The black line denotes original observations and the blue line denotes the sum signals from IMF-7 to IMF-12 of original pressure observations based on the EEMD; (b) Gravity residuals during 10 days. The black line denotes original observations and the blue line denotes the sum signals from IMF-7 to IMF-12 of gravity residual based on the EEMD; (c) Normalized amplitude spectrum of sum signals from IMF7 to IMF12. The blue line denotes atmospheric pressure and the red line denotes gravity residuals

那么从0.1 mHz到1 mHz的频段上，气压变化对重力观测又会带来多大的影响呢？这就需要分析处于更高频段上的IMF之间的相互关系了.这里我们以IMF-5为例，给出了气压变化和重力之间在时域(时间序列)和频域(振幅谱)里的比较结果，见图 6.图 6a展示了法国Strasbourg台站在上述平静期里IMF-5的局部信号.在为期2天半的时间序列里，重力变化与气压变化走势完全一致，曲线基本重合，其中气压变化已乘上IMF-5的气压导纳值-3.989 nm·s^-2·hPa^-1(见表 1)转化为了重力值.图 6b是气压变化和重力残差的IMF-5的归一化振幅谱，计算振幅谱的数据长度仍为十天.振幅谱表明两种信号的IMF-5的频段主要集中在0.2~0.4 mHz，同时两类信号的谱峰也较为吻合.

图 6

Fig. 6

图 6 重力残差和气压观测的IMF-5之间的比较 (a)局部时间序列：红色虚线为气压变化，已乘上气压导纳值-3.989 nm·s-2·hPa-1转化为重力值，黑色实线为重力变化，两者的相关系数为-0.9213；(b)归一化振幅谱：红色和黑色实线分别表示气压和重力的归一化振幅谱，并均经过了11个点的Parzen谱窗的平滑处理. Fig. 6 Comparison of IMF-5 between gravity residuals and atmospheric pressure (a) Part of time series. The red dotted line denotes the atmospheric pressure multiplied by admittance -3.989 nm·s-2·hPa-1, the black solid line denotes the gravity residual and their correlation coefficient is -0.9213;(b) Normalized amplitude spectrum. The red solid line and the black solid line denote atmospheric pressure and gravity residual respectively and they are smoothed by an 11-point Parzen spectral window.

从上述分析中可以看到，在不同频段上的气压导纳值是有一定区别的，长周期IMF(从IMF7到IMF12)和IMF-5的气压导纳值相差接近1.因此如果使用单一的气压导纳值进行气压改正可能会在不同频段引入较大的气压噪声.分析结果表明：通过EEMD分解我们可以充分详细地了解台站局部气压在不同频段对重力观测造成的影响，以此为基础对地震后的重力观测进行气压改正是我们降低气压噪声，提高信噪比，高精度检测微弱低频自由振荡信号的重要保证.

利用本文引入的EEMD方法改正台站气压变化对重力观测的影响，选取2004年苏门答腊大地震后GGP台网中6个超导重力台站的8组高质量分钟值观测数据，检测低频自由振荡的频谱分裂现象，同时检验利用EEMD方法进行气压改正的实际效果.台站和仪器信息如下：法国的Strasbourg台站(代码：ST，仪器编号：C026，下同)，澳大利亚的Canberra台站(CB，C031)，意大利的Medicina台站(MC，C023)，比利时的Membach台站(MB，C021)，以及两个双球型SG台站，德国的Bad Homburg台站(H1, CD030_L; H2, CD030_U)和南非的Southland台站(S1, CD037_L; S2, CD037_U).

首先选取地震发生之前的一段平静期(为期10天)来评估各台站气压变化在不同频段对重力观测的影响，各台仪器的计算结果列于表 2.表 2的第一列给出了利用傅氏振幅谱估计的每个IMF对应的优势频率范围.作为比较，由重力残差和气压观测的原始信号计算得到的气压导纳值和它们的相关系数一并列于表 2中.

表 2 Table 2

表 2 8组超导重力观测在不同频段的气压导纳值以及气压变化与重力的相关系数 Table 2 Admittances of 8 SGs in different frequency bands and correlation coefficients between atmospheric pressure and gravity 仪器 ST H1 H2 CB MC MB S1 S2 原始重力残差 c 2.763 2.798 2.824 3.125 2.760 2.806 3.385 3.326 Corr. 0.995 0.995 0.996 0.991 0.985 0.990 0.949 0.937 IMF-2 c 0.417 0.694 0.695 0.513 0.606 0.441 0.526 0.365 (1.5~3.5 mHz) Corr. 0.054 0.144 0.139 0.078 0.069 0.032 0.102 0.071 IMF-3 c 1.794 2.304 2.227 2.385 1.656 2.475 2.953 2.676 (0.8~1.5 mHz) Corr. 0.535 0.557 0.542 0.603 0.285 0.438 0.694 0.663 IMF-4 c 3.550 3.483 3.468 3.415 3.215 3.355 2.983 2.848 (0.4~0.8 mHz) Corr. 0.857 0.794 0.780 0.842 0.640 0.773 0.805 0.810 IMF-5 c 3.989 3.970 3.909 4.143 3.907 3.982 3.176 2.904 (0.2~0.4 mHz) Corr. 0.921 0.903 0.910 0.889 0.859 0.897 0.848 0.830 IMF-6 c 3.567 3.301 3.482 3.013 3.588 3.447 3.123 2.486 (0.1~0.2 mHz) Corr. 0.948 0.891 0.915 0.801 0.825 0.938 0.773 0.769 IMF-7到IMF-12 c 3.094 3.105 3.222 3.163 2.469 2.946 3.505 3.415 (小于0.1 mHz) Corr. 0.991 0.994 0.992 0.982 0.940 0.976 0.932 0.931 注：导纳值和相关系数均取了绝对值，c为气压导纳值，单位：nm·s-2·hPa-1，Corr.为相关系数.

表 2 8组超导重力观测在不同频段的气压导纳值以及气压变化与重力的相关系数 Table 2 Admittances of 8 SGs in different frequency bands and correlation coefficients between atmospheric pressure and gravity

EEMD方法是在原信号的基础上通过加白噪声再求集合平均(ensemble average)的方式作为最终的分解结果，当具体应用到SG数据的分解时，会对最先分解得到的处于高频段的IMF-1带来较大的扰动.因此每次计算得到的IMF-1的气压导纳值会有较大出入，结果不稳定.幸运的是我们主要关注的是低频自由振荡信号，频率段处于1.5 mHz以下，IMF-1不稳定的结果不会给我们带来不好的影响，因此这里将IMF-1的相关结果忽略.实际上，在IMF-1(优势频率集中在3.5~8.3 mHz，见图 3)和IMF-2所在的频段上，重力与局部气压波动的相关性很低，气压对重力变化的影响较小，因而不必在此频段进行气压改正.

从表 2可看出，在频率范围低于1.5 mHz的低频自由振荡频段，局部气压导纳值具有很强的频率依赖性，在不同频段的区别较大.如在IMF-5所在的低频段(0.2~0.4 mHz)，气压导纳值的绝对值较大，接近4 nm·s^-2·hPa^-1，气压变化与重力变化相关性高，相关度基本在85%以上.而在IMF-3所在的高频段(0.8~1.5 mHz)气压导纳值的绝对值相对较小，气压变化与重力变化相关性相对较低.IMF7到IMF12的和信号主要包括了原始信号去趋势后的长周期变化，与原始信号在变化趋势上较为相似，但它们各自计算得到的气压导纳值仍然有所不同.因此如果利用由原始信号得到的时域气压导纳值，进行不分频段的局部气压改正，则不仅不能完全消除低频段的气压影响，而且还会在较高频段注入气压噪声.

表 2中的每一行数据基本处于相同的区间范围，例如IMF-2的相关系数在0.1上下浮动，浮动范围大部分不超过0.05，其导纳值均在0.5 nm·s^-2·hPa^-1左右；IMF-5的相关系数在0.9上下浮动，浮动范围大部分也不超过0.05，其导纳值多数在4 nm·s^-2·hPa^-1左右，并且都是各台仪器所有子频段上的最大值.在其他的IMF中，不论是气压导纳值还是相关系数也都有类似的表现.每个台站内部的各个IMF相对变化趋势也大致相同.参数变化的这种一致性可能反应了气压对重力观测影响的某种内在本质.

对震后去掉潮汐信号以后的重力残差记录和同步气压观测分别做EEMD分解，将信号分解到以上几个频段，然后利用各个频段上的气压导纳值分别进行气压改正.我们利用傅氏变换将去掉气压影响的每个IMF转化到频率域得到他们各自的振幅谱，图 7给出了IMF3到IMF5的振幅谱结果.结果显示将处于不同频率段上的简正模分解到了不同的IMF中，使我们可以单独利用一个IMF或几个IMF的组合来评估某一个简正模式.图 8给出了这三个IMF的和信号对应的振幅谱.可以看到，在低于1.5 mHz的频段上，除了检测到了所有的球型振荡以外，还清楚地检测到了多个环型振荡的耦合现象.

图 7

Fig. 7

图 7 经气压改正后IMF-3到IMF-5的振幅谱 Fig. 7 Amplitude spectra of IMFs from IMF-3 to IMF-5 after atmospheric pressure correction based on the EEMD

图 8

Fig. 8

图 8 由IMF-3到IMF-5的和信号求得的振幅谱 竖直虚线给出了低频段所有的球型振荡和环型振荡的PREM模型理论值.结果显示除了得到了所有球型振荡的谱峰以外，还清楚地检测到了多个环型振荡的耦合现象. Fig. 8 Amplitude spectrum of the sum signals from IMF-3 to IMF-5 The predicted frequencies for PREM of all spheroidal and toroidal modes below 1.5 mHz are indicated by vertical dotted lines. In addition that all spheroidal modes are obtained, the plot also shows that some spheroidal-toroidal mode couplings are clearly detected.

本节选取了两个周期最长的基频球型振荡₀S₂和₀S₃，利用频率域自回归分析方法(Chao and Gilbert, 1980)估计谱峰的本征频率，以此检测它们的频谱分裂现象.实际上我们只需要选取相应频段上单独的一个IMF即可得到相关的简正模，例如₀S₂的理论频率在0.3 mHz左右，只需选取IMF-5就可将其完整检测到，因此本文选择IMF-5和IMF-4分别研究两个低阶基频球型振荡₀S₂和₀S₃的频谱分裂现象.

表 3给出了8组数据对₀S₂分裂谱峰频率的检测结果，同时计算了它们的加权平均值，并且列出了两篇文献中基于2004年Sumatra地震后的超导重力观测数据的相关结果(Rosat et al., 2005; Hu et al., 2006)及两种模型(1066A和PREM)的理论预测值.结果显示，除澳大利亚Canberra台站m=0的谱峰未被激发以外，其他各台站的每个谱峰均得到了相应的检测结果.另外南非Southland台站的两组数据检测到的m=0的谱峰具有较大的误差限，从其谱图中可以看到，在该台激发的m=0的谱峰并不明显.对比来看，本文结果和前人的观测值及两种模型的理论预测值都保持较好的一致性.

表 3 Table 3

表 3 0S2分裂谱峰的频率观测值和两个模型的理论预测值 Table 3 Observed frequencies of the splitting of 0S2 and their theoretical predictions based on two Earth models 台站 0S2的频率(μHz) m=-2 m=-1 m=0 m=1 m=2 平均值 ST 299.802±0.019 304.468±0.059 309.240±0.019 313.786±0.059 318.417±0.013 309.142±0.038 H1 299.796±0.020 304.391±0.050 309.354±0.033 313.907±0.037 318.301±0.036 309.150±0.039 H2 299.823±0.025 304.530±0.049 309.352±0.029 313.892±0.039 318.279±0.031 309.176±0.039 MC 299.748±0.025 304.419±0.057 309.128±0.043 313.775±0.039 318.250±0.026 309.064±0.040 MB 299.726±0.058 304.089±0.055 309.387±0.014 313.707±0.046 318.031±0.032 308.988±0.039 S1 299.843±0.029 305.062±0.206 309.293±0.621 313.755±0.047 318.329±0.019 309.256±0.080 S2 299.799±0.025 305.081±0.078 309.404±0.319 313.975±0.029 318.309±0.017 309.314±0.047 CB 299.983±0.032 304.494±0.055 — 313.673±0.057 318.470±0.018 309.155±0.038 加权平均 299.815±0.009 304.493±0.022 309.313±0.012 313.827±0.015 318.326±0.008 309.145±0.015 Rosat等(2005) 299.970±0.066 304.620±0.110 309.290±0.061 313.890±0.096 318.450±0.033 309.240±0.097 Hu等(2006) 299.970±0.025 304.645±0.056 309.341±0.047 313.892±0.064 318.454±0.034 309.260±0.034 1066A* 299.9310 304.7944 309.4921 314.0210 318.3800 309.3237 PREM* 299.8479 304.7817 309.3972 314.0127 318.3099 309.2699 *:基于1066A和PREM地球模型计算得到的理论预测值.

表 3 0S2分裂谱峰的频率观测值和两个模型的理论预测值 Table 3 Observed frequencies of the splitting of 0S2 and their theoretical predictions based on two Earth models

₀S₂是周期最长的自由振荡信号，其振荡频率相对最为接近地球的自转频率，因而₀S₂的分裂最为明显.图 9给出了利用4个台站的IMF-5得到的振幅谱，结果显示单独一个IMF就能得到全部的分裂谱峰.

图 9

Fig. 9

图 9 4个台站的IMF-5得到的0S2的谱峰分裂频谱图 观测数据为苏门答腊地震5小时以后，240个小时的分钟值记录，频谱计算前均加了hanning窗，4个台站的5个谱峰都完全地分裂开来. Fig. 9 Spectra of IMF-5 of gravity residuals showing frequency splitting of 0S2 into 5 completely resolved singlets at 4 individual stations The observations are based on 240 hours minute-records after 5 hours since the 2014 Sumatra earthquake. The Hanning window is applied before calculating frequency spectra.

利用法国Strasbourg台在2004年苏门答腊大地震5小时后为期546小时的分钟值记录检测了₀S₃的频谱分裂现象，并利用重力残差的IMF-4给出了该台检测得到的分裂谱峰(见图 10)，其中m=0的谱峰没有被激发.表 4给出了₀S₃分裂谱峰的频率检测结果，同时表中还列出了两个前人的检测结果和Roult等(2010)利用PREM+re地球模型计算得到的₀S₃分裂谱峰本征频率的理论值.

图 10

Fig. 10

图 10 法国Strasbourg台站2004年苏门答腊大地震5小时后为期546个小时的分钟值记录对0S3的检测结果 竖直虚线给出了PREM+re地球模型的理论预测值，频谱图显示m=0的谱峰在该台未被激发. Fig. 10 Observed frequencies of the splitting of 0S3 based on 546 hours minute-records after 5 hours since the 2014 Sumatra earthquake at Strasbourg The predicted frequencies for PREM+re are indicated by vertical dotted lines. The singlet m=0 was not excited at Strasbourg.

表 4 Table 4

表 4 0S3的分裂谱峰的频率观测值和PREM+re模型的理论预测值 Table 4 Observed frequencies of the splitting of 0S3 and their theoretical predictions based on the PREM+re model 台站 0S3的频率(μHz) m=-3 m=-2 m=-1 m=0 m=1 m=2 m=3 法国Strasbourg 461.667±0.027 464.205±0.023 466.411±0.012 - 470.782±0.013 472.699±0.033 474.719±0.025 Rosat等(2005) 461.67±0.030 464.24±0.400 466.39±0.018 - 470.84±0.011 472.66±0.048 474.74±0.021 Buland等(1979) 461.67 464.08 466.17 468.83 470.91 473.17 474.50 PREM+re 461.986 464.123 466.272 468.439 470.606 472.791 474.989

表 4 0S3的分裂谱峰的频率观测值和PREM+re模型的理论预测值 Table 4 Observed frequencies of the splitting of 0S3 and their theoretical predictions based on the PREM+re model

这里需要指出，虽然普通的气压改正方法也能够清楚地观测到₀S₂的5个分裂谱峰，但是更为精细的气压改正可以提高分裂谱峰的检测精度.表 3中给出的Rosat等(2005)和Hu等(2006)₀S₂分裂谱峰的频率观测值是利用2004年Sumatra地震后同一个超导重力台站(法国的Strasbourg)的观测记录检测得到.然而Rosat等(2005)是基于一个名义气压导纳值进行气压改正后得到的检测结果，其给出的误差水平相对较大；而Hu等(2006)是利用小波分解的方法对重力观测记录进行了气压改正，观测值的误差水平就大大降低，并与本文给出的结果相当.另外，在数据长度的选择上，普通的气压改正方法需要更长的数据长度才能完全将分裂谱峰完全分开，例如Rosat等(2005)选择了528个小时的数据检测得到了Bad Homburg台站观测的₀S₂的5个分裂谱峰，利用888个小时的数据检测得到了Strasbourg台站观测的₀S₃的分裂谱峰.而本文仅分别利用了240个小时和546个小时的数据就完成了对₀S₂和₀S₃谱峰分裂的检测工作.

一阶球型振荡₂S₁的周期仅次于₀S₂，地震激发的₂S₁的信号强度往往非常微弱，地震仪记录的数据由于受到噪声的影响，从未检测到₂S₁的分裂现象.因此利用在低频段具有更低噪声水平的超导重力数据检测₂S₁引起了国内外学者的广泛兴趣，科学家通过不同方法尝试检测这一简正模式及其频谱分裂现象(Courtier et al., 2000; Rosat et al., 2005, 2003b; Hu et al., 2006).相关结果显示₂S_l的分裂谱峰十分微弱，尤其是m=0的谱峰几乎完全被气压噪声掩盖，因此在检测₂S_l的谱峰分裂现象时，准确的气压改正是必不可少的检测手段.

由于受到频率分辨率和品质因子的影响，过长或过短的数据均不利于检测₂S₁的分裂谱峰的频率，必须选择合适的数据长度才能更好地分离出₂S₁的3个分裂谱峰.经过多次尝试，本文最终采用了苏门答腊大地震5小时后168个小时的重力观测数据.由于₂S₁的理论频率在0.4 mHz左右，因此，选择分别经过气压改正后的IMF-5与IMF-4的和信号来研究₂S₁的频谱分裂现象.

首先选取ST台站单台超导重力仪的观测记录，并利用它在两个频段上的气压导纳值，IMF-4 (0.4~0.8 mHz)为-3.550 nm·s^-2·hPa^-1和IMF-5 (0.2~0.4 mHz)为-3.989 nm·s^-2·hPa^-1，对重力残差信号进行气压改正.另外，作为比较我们还选择了单一的时域气压导纳值-2.763 nm·s^-2·hPa^-1 (即由原始重力残差在时域中计算得到的气压导纳值)直接对原始重力残差进行气压改正.图 11a给出了这两种方法在相应频段的检测结果.从图中我们可以看到基于EEMD方法进行的气压改正效果与普通方法相比，显著提高了₂S₁谱峰的分辨率，尤其m=0谱峰的信噪比明显得到提高，使我们可清楚地观测到₂S₁的3个分裂谱峰.未经气压改正的检测结果也在图中一并给出，该结果中m=0的谱峰基本看不到，再次说明了精细化气压改正的必要性.

图 11

Fig. 11

图 11 (a) 2004年苏门答腊大地震5小时后，法国Strasbourg超导重力观测记录对2S1的检测结果，数据长度为168小时：黑色实线为EEMD方法进行气压改正后的结果，红色实线为普通方法(即利用单一的时域气压导纳值-2.763 nm·s-2·hPa-1)进行气压改正后的结果，蓝色虚线为未经气压改正得到结果；(b)三组数据的积谱结果，除法国Strasbourg外，另两组数据均来自于德国的Bad-Homburg台站 Fig. 11 (a) Observation of 2S1 at Strasbourg using 168 hours minute-records after 5 hours since the 2014 Sumatra earthquake. The black solid line denotes the observed results after atmospheric pressure correction based on the EEMD, the red solid line denoted the observed results after traditional atmospheric pressure correction, i.e. using a single nominal admittance of -2.763 nm·s-2·hPa-1 and the blue dotted line denotes observed results without atmospheric pressure correction; (b) Product spectrum of 3 sets of gravity data recorded by SGs at Strasbourg and Bad-Homburg

除此之外利用EEMD方法对Bad-Homburg的两组观测数据(H1和H2)进行气压改正，在它们各自的IMF的和信号中也检测到了₂S₁的3个分裂谱峰.多台站迭积的方法能够压制台站噪声，放大公共信号，提高信噪比，可以更精确地确定简正模的本征频率.因此这里采用一种简单的频率域迭积方法—积谱密度估计，迭积三组观测的频谱结果，以此减少台站环境噪声对₂S₁检测的影响，从而得到更准确的频率检测结果.积谱密度估计作为一种频率域的迭积方法，将来自不同台站观测数据的频谱相乘后再求其几何平均，在信噪比较高的情况下，多个台站的积谱只显示这些台站记录到的公共信号，在一个或几个台站没有记录到的信号将被消除，而且影响单个台站的系统误差将被压制(Hu et al., 2006).图 11b给出了三组观测数据的积谱密度估计结果.从图中可以看到，通过迭积的方法分裂谱峰所在频段附近的噪声水平明显降低，信噪比得到提高.

表 5给出了各台站及其积谱密度估计的结果，表中同时列出了Rosat等(2003a)基于2001年秘鲁地震后5组SG观测数据得到的迭积结果，Hu等(2006)采用小波方法进行气压改正后给出的2004年苏门答腊地震后三组观测的积谱结果，以及基于1006A和PREM模型计算得到的₂S₁分裂频率的理论值.从表 5本文的计算结果中可以看到，由于₂S₁的信号强度较弱，其信噪比较低，单台站得到的分裂谱峰的误差评定结果虽然没有₀S₂好(尤其是m＝0的情况)，但仍优于Rosat等(2003a)利用5组SG观测数据得到的迭积结果；本文三组观测的积谱结果，通过放大谐信号提高了信噪比，其误差评定结果有较大提高，与Hu等(2006)利用小波方法进行气压改正后给出积谱结果相当.另外由积谱结果得到的3个分裂谱峰的检测频率与单台站的检测频率相比，与理论值更为相符.

表 5 Table 5

表 5 2S1的分裂谱峰的频率观测值和两个模型的理论预测值 Table 5 Observed frequencies of the splitting of 2S1 and their theoretical predictions based on two Earth models 2S1的频率(μHz) m= -1 m=0 m=1 ST 397.412±0.145 404.193±0.150 410.687±0.097 H1 398.175±0.078 404.078±0.190 410.712±0.040 H2 398.404±0.106 403.610±0.290 410.707±0.071 三组数据的积谱结果 398.507±0.083 403.677±0.093 410.979±0.049 Rosat等(2003a)5组数据迭积 398.600±0.190 404.900±0.210 411.100±0.180 Hu等(2006)3组数据迭积 398.854±0.048 404.331±0.075 410.835±0.041 1066A 398.708 404.690 410.880 PREM 398.750 404.727 410.948

表 5 2S1的分裂谱峰的频率观测值和两个模型的理论预测值 Table 5 Observed frequencies of the splitting of 2S1 and their theoretical predictions based on two Earth models

在低于1 mHz的频段里，气压变化对重力观测的影响不可忽略，且具有频率依赖特性.不同频段上的气压导纳值存在一定差异，用单一的气压导纳值进行气压改正过于粗略，不仅不能完全消除气压变化对重力观测带来的影响，还可能在部分频段引入较大的气压改正噪声.

本文基于EEMD方法，将重力和气压变化分解成处于不同频段的本征模态函数(IMF)，通过计算具有频率依赖特性的气压导纳值，更精细地消除气压变化对重力观测的影响.利用2004年苏门答腊地震后GGP台网中6台超导重力仪的8组观测数据检测了频率小于1.5 mHz的低频地球自由振荡及其频谱分裂现象.结果证明：利用频率依赖的气压导纳值进行气压改正后得到的重力残差，能更有效的用于检测低频自由振荡信号，得到精度更高的低阶球型振荡的分裂谱峰，同时还检测到了部分环型振荡在重力观测中的耦合现象.对低频地球振荡的高分辨率检测结果验证了基于EEMD分解提出的气压改正方法的有效性，同时再次证明了超导重力仪在低频地球自由振荡检测中的优越性.

不论是检测₀S₂还是₀S₃的分裂频率，本文得到的观测结果虽然都与基于地球模型得到的理论预测值接近，但仍然存在一定的差异.即使考虑到误差范围，大部分台站₀S₂分裂频率的平均值均小于由模型给出的理论预测平均值.在如此高信噪比的情况下得到的检测结果，仍然与地球模型的理论预测值存在一定的出入，说明现有一维地球模型还有进一步改进的空间.在低于1 mHz的频段里还存在大量微弱的的核模信号，例如地球固态内核平动振荡(即Slichter模)，这些核模信号的高精度检测同样依赖于精细的气压改正及有效的检测方法的应用，相关研究将在今后的工作中进一步开展.