ABSTRACT

The research of earth free oscillation has become one of important means to ascertain the internal structure of earth. This paper uses the marginal spectrum method based on Hilbert Huang transform for the first time, and the Tonghe station’s borehole strain data without removing solid tide. Then we extracted earth free oscillation mode that are motivated by Wenchuan M8.0 earthquake, Chile M8.8 earthquake and the east coast of Honshu of Japan M9.0 earthquake, and the results are very close to the PREM mode. Among which, the Wenchuan earthquake is tested 24 fundamental frequency modes, the Chile and Japan earthquakes are both 76, which consist of ₀S₀ and ₀S₂~₀S₇₆. The marginal spectrum method can effectively extract the earth free oscillation because of these strong earthquakes’ mode have spectrum divide phenomenon.

Keywords:Hilbert Huang Transform, The Marginal Spectrum, The Strong Earthquakes, Free Oscillation

采用边际谱方法提取地球自由振荡

李继业¹，武晓军²，马龙辰¹，任建辉³，张彦吉²，孙鹏宇⁴，孙强⁵

¹黑龙江省地震局，黑龙江 哈尔滨

²哈尔滨市通河地震台，黑龙江 哈尔滨

³哈尔滨市防震减灾技术中心，黑龙江 哈尔滨

⁴哈尔滨市延寿地震台，黑龙江 哈尔滨

⁵哈尔滨市阿城地震台，黑龙江 哈尔滨

收稿日期：2017年6月9日；录用日期：2017年6月25日；发布日期：2017年6月28日

摘 要

地球自由振荡的研究已成为探查地球内部结构的重要手段之一。本文首次采用基于希尔伯特–黄变换的边际谱方法，利用通河台钻孔应变资料，在没有对数据进行去固体潮的前提下，提取了2008年5月12日汶川M8.0级，2010年2月27日智利M8.8级和2011年3月11日日本本州东海岸附近海域M9.0级强震的地球自由振荡。通过与PREM模型相对比，结果十分接近。其中，汶川地震检测到了24个基频信号，智利地震、日本地震检测到全部76个基频信号，包括₀S₀、₀S₂~₀S₇₆。汶川地震检测到的₀S₂、智利地震检测到的₀S₆和₀S₇等和日本地震检测到的₀S₂和₀S₃等均存在频谱分裂现象。证明边际谱方法可以有效的提取地球自由振荡。

关键词 :希尔伯特–黄变换，边际谱，强烈地震，自由振荡

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

大的地震发生时，除产生涉及地球局部运动的体波和面波外，还能激发弹性球体的地球产生固有振动，这就是地球的自由振荡。这种振动通常会持续数天甚至数周时间，周期从几十秒到近一个小时。地球自由振荡理论上是长周期面波的推广，可以作为独立于体波和面波之外的另一种地球内部结构的约束条 [1] [2] 。地球自由振荡是探查地球内部结构的重要手段之一，通过实际观测和理论分析与对比，可以研究地球内部的各向异性、非弹性以及震源机制的性质等。

自由振荡理论研究可追溯至19世纪末，Lamb [3] 详细讨论了均匀球体简单模型下自由振荡模式的分布规律。Love [4] 在Lamb的研究成果基础上，讨论了重力作用下可压缩球体振荡问题，给出了极为精确的₀S₂固有周期。直到1954年由Benioff [5] 首次在他自己设计的应变地震仪上，发现1952年11月4日勘察加地震所激发57分钟的长周期地球自由振荡。20世纪80年代初，Dziewonsiki et al. [1] 又在前人研究基础上，提出了初步参考地球模型(Preliminary Reference Earth Model，简称PREM)，计算了各向同性球体和纵向层状结构下球型和环型振荡的理论本征周期，并给出实际观测结果，极大丰富了自由振荡理论的发展。

长期以来，地球自由振荡提取方法和检测对象也在不断完善和扩展。Ness等 [6] 采用二阶Chebyshev数字高通滤波器提取了地球自由振荡资料；Van Camp [7] 用LSQ高通滤波器和LSQ低通滤波器消除了1998年Baleny岛M8.1级地震超导重力观测资料中重力潮汐影响，提取了地球的球型自由振荡。雷湘鄂等 [8] 运用20阶多项式拟合固体潮观测值，提取了超导重力仪在2001年6月秘鲁M7.9级地震激发的全部₀S₀~₀S₃₂球型振荡振型，并观测到了₀S₀、₀S₃的谱线分裂现象。万永革等 [9] 利用中国数字地震台网(CDSN)波形资料，采用功率谱密度估计方法，在没有去除固体潮影响的情况下，提取了2001年昆仑山口西M8.1级地震激发的₀S₄-₀S₇₆地球球型自由振荡，取得较好的效果。严珍珍等 [10] 采用谱元法，结合高性能并行计算，数值模拟了特大地震激发的弹性波在地球内部传播过程。任佳等 [11] 采用垂直摆倾斜仪和水管倾斜仪数字化观测资料，利用功率谱密度估计方法，在没有对资料进行去固体潮处理情况下，准确获得了2008年汶川M8.0级大地震激发的₀S₆~₀S₃₂基频球型自由振荡。杨跃文等 [12] 利用功率谱密度估计方法，在没有对资料进行处理的情况下，准确获得了2010年2月27日智利M8.8级大地震激发的₀S₅~₀S₃₅基频球型自由振荡。吴海波等 [13] 采用功率谱密度方法，提取了JCZ-1型地震仪检测到的地球环形自由振荡。

本文首次应用基于希尔伯特–黄变换的边际谱方法，采用同一钻孔应变的数字化资料，提取3次8级以上强震激发的地球自由振荡。通过对信号进行希尔伯特–黄变换(Hilbert-Huang Transform，简称HHT)，得到一种新的信号时频描述方式，被称为Hilbert谱，再对信号进行时间域积分，得到类似于Fourier谱的Hilbert边际谱。边际谱方法相比传统的傅里叶功率谱估计方法对数据要求更短，且其频率分辨率和估算精度都要高于傅里叶功率谱估计。

2. Hilbert-Huang变换原理

希尔伯特–黄变换是一种非线性、非平稳信号的处理方法。最早由美国华裔科学家Huang等1998年提出，他对Fourier变换基本信号和频率定义做了创造性的改进，不在认为信号是由一组固定频率和振幅的正弦信号的加权和，而是由固有模态函数(intrinsic mode function, IMF)的信号组成。希尔伯特–黄变换(HHT)方法由经验模态分解(EMD)方法和Hilbert变换组成 [14] 。

2.1. 经验模态分解(Empirical Mode Decomposition，简称EMD)

EMD方法即Huang变换，它是依据信号本身时间尺度特征，将信号分解为不同时间尺度且满足以下两个定义条件的一组IMF，每个IMF可以被认为是信号中固有的一个模态函数：1) 对于一列数据极值点和过零点数目必须相等或至多相差一点；2) 在任意点，由局部极大点和极小点构成的两条包络线平均值为零。

EMD分解过程：1) 找出信号x(t)所有极大点和极小点，将其用三次样条函数分别拟合为原数据序列上、下包络线，上、下包络线的均值为平均包络线m₁，将原数据序列减去m₁可得到一个去掉低频的新数据序列h₁。一般h₁不是一个平稳数据序列，为此重复以上过程n次，使所得平均包络线趋于零，此时h_1n就是第一个IMF(c₁)，它表示信号数据中最高频成分。2) 用x(t)减去c₁得到一个去掉高频成分的新数据序列，重复步骤1)，得到一系列c_n和最后一个不可分解的序列r_n，r_n代表x(t)的均值或趋势项。那么原序列x(t)可表示为IMF分量和一个残余项的和：

(1)

信号经分解后得到多个IMF的组合，对每个IMF分量进行Hilbert变换，即可得到每个IMF分量的瞬时频率，综合所有IMF分量的瞬时频谱得到Hilbert谱。

获得IMF分量以后，对每一阶IMF作Hilbert变换，设为y(t)：

(2)

x(t)和y(t)共同组合为一解析信号z(t)，采用极坐标形式表示：

(3)

其中：

(4)

(5)

定义瞬时频率：

(6)

由式(6)可看出，是时间t的单值函数，即某一时间对应某一频率，每个IMF序列在每一点的频率唯一。对每一阶IMF作Hilbert变换，并求出相应解析函数的幅值谱和瞬时频率，从而原始信号可以表示为：

(7)

上式频率和幅值都是时间的函数，可构成幅值、频率、时间的三维时频图即Hilbert幅值谱，简称Hilbert谱，记为：

(8)

基于Hilbert谱公式(8)，通过对时间的积分，可以定义边际谱为：

(9)

边际谱表达了每个频率在全局上的幅度，它代表了在统计意义上的全部累加幅度。

2.2. 边际谱的物理意义

Hilbert谱中的频率还是边际谱中的频率(即瞬时频率)，其意义与Fourier分析中的频率完全不同。Fourier变换是将任何信号分解为正弦信号的加权和，而每一个正弦信号对应着一个固定的频率(Fourier频率)和固定的幅值，因此，对于频率随时间变化的非平稳信号，Fourier变换就无能为力了。而在边际谱中，由(9)式可知边际谱提供了对频率总的幅值分布，它代表了整个信号在时间跨度上的幅值累积效应。某个频率值仅仅意味着有很大的可能性存在这一频率成分的正弦或余弦波，而其发生时刻，则在Hilbert谱中给出了精确的定位。借助于信号IMF组合，幅值调制和频率调制被清晰地分离开来。由于跨越了Fourier变换中要求的常幅值和常频率限制，也就可以突破Fourier变换仅对平稳信号有效的不足之处，使得Hilbert-Huang变换能够成功地应用于非线性、非平稳信号的处理。边际谱是Hilbert谱对时间的积分。对任意固定频率ω的时间积分就是对频率ω所有时刻t所对应的幅值求和。因此，瞬时频率ω的边际谱的含义是：信号中瞬时频率ω的总幅值的大小。可以这样认为，瞬时频率表示信号交变快慢的物理量，任一瞬时频率都有一定的幅值，将不同时间点的但具有相同瞬时频率的所有幅值加起来，就是信号中该频率的总幅值，即边际谱线的高度。但该频率不一定在所有时刻都存在，也不一定只在某一时刻出现，而是可能以不同或相同幅值出现几次。事实上，这一物理意义与Fourier频谱的物理意义是一致的，只不过在Fourier频谱中，对任意频率都要求在整个信号时间跨度上具有相同的幅值 [15] [16] 。

3. 观测点概况及资料选取

通河体应变观测井(128.65˚E，46.08˚N)位于伊春–延寿地槽褶皱带北段，地质构造属兴安岭–内蒙古地槽褶皱系、小兴安岭–张广才岭槽地褶带。主要构造形迹有北东向依舒断裂，北西向尚志断裂和岔林河断裂。深部构造莫霍面嘉荫–舒兰深度陡变带与鸡西–塔溪深度陡变带在通河附近相交，呈“X”型相交叉分布，将黑龙江省大兴安岭以东广大地区，分割成互为对顶分布的两个上地幔隆起区和两个上地幔坳陷区，构成了莫霍面起伏变化的基本轮廓 [17] 。观测井濒临郯庐断裂带北段分支之一依舒断裂，处于松嫩平原与东部山地间丘陵区，属张广才岭余脉。周围地形地貌特征属花岗岩类构造低山，台基岩性为坚硬完整的燕山期花岗岩。通河体应变观测井深81 m，探头位置79 m，安装中国地震局地壳所研制的TJ-Ⅱ型体积式应变仪，于2005年10月安装，2006年6月正式观测运行，观测精度为10⁻⁹，动态范围大于2 × 10⁴，采样率为分钟值采样。

近年来全球大地震活动频繁，2008年5月12日汶川M8.0级，2010年2月27日智利M8.8级和2011年3月11日日本本州东海岸附近海域M9.0级地震，这些罕见的大地震为人们研究地球自由振荡提供了良好的机会(表1)。资料选取这三次强震前半小时和地震后3天原始观测数据，长度4350分钟。所选观测数据固体潮清晰可辨，数据连续性好，内在质量高，受环境和人为影响小，信噪比高。为了保证数据原始可靠，没有对原始数据进行任何形式滤波处理，包括固体潮信息(图1)。通过希尔伯特–黄变换求解边际谱，提取这3次罕见强震激发的地球自由振荡信息。

4. 边际谱提取结果分析

利用通河体应变三次强震数据，采用基于Hilbert-Huang变换的边际谱方法，计算了200 μHz~8200 μHz (周期122 s~5000 s)边际谱(图2)，按照₀S₀~₀S₇₆基频球型自由振荡分别进行检测。理论本征周期选用PREM模型各向异性的理论本征周期(图中虚线标注部分)。

表1. 三次强震基本参数

注：以上资料来自中国地震台网中心地震数据管理与服务系统http://www.csndmc.ac.cn/newweb/index.jsp。

图1. 通河体应变对三次强震的分钟值记录

图2. 200~4100 μHz自由震荡边际谱。(a) 200~1000 μHz；(b) 1000~2100 μHz；(c) 2000~3000 μHz；(d) 3000~4100 μHz

汶川地震检测到了24个基频信号，分别在200~1000 μHz频段检测到₀S₀、₀S₂、₀S₃和₀S₄四个，但没有检测到₀S₅ (图2(a))；1000~2000μHz频段检测到₀S₆~₀S₁₂七个(图2(b))；2000~3000 μHz频段检测到₀S₁₃、₀S₁₄和₀S₁₈三个(图2(c))；3000~4000μHz频段检测到₀S₂₄和₀S₂₅两个(图2(d))；4000~5000 μHz频段仅检测到₀S₃₆ (图3(a))；5000~6000μHz频段检测到₀S₄₆、₀S₄₇、₀S₄₉和₀S₅₄四个(图3(b))；6000~7000μHz频段未检测到信号(图3(c))；7000~8200μHz频段检测到₀S₇₄、₀S₇₅和₀S₇₆三个(图3(d))。其中₀S₃误差为3.26%，₀S₄误差为1.45%，₀S₇误差为1.21%，₀S₁₁误差为1.05%，₀S₁₂误差为1.45%，₀S₁₃误差为1.72%和₀S₁₈误差为1.48%，其它基频信号误差均在±0.8%之间。

智利地震、日本地震检测到76个基频信号，包括在200~1000 μHz频段检测到₀S₀、₀S₂、₀S₃、₀S₄和₀S₅五个(图2(a))；1000~2000 μHz频段检测到₀S₆~₀S₁₂七个(图2(b))；2000~3000 μHz频段检测到₀S₁₃~₀S₂₁九个(图2(c))；3000~4000 μHz频段检测到₀S₂₂~₀S₃₂十一个(图2(d))；4000~5000 μHz频段检测到₀S₃₃~₀S₄₃十一个(图3(a))；5000~6000 μHz频段检测到₀S₄₄~₀S₅₄十一个(图3(b))；6000~7000 μHz频段检测到₀S₅₅~₀S₆₅十一个(图3(c))；7000~8200 μHz频段检测到₀S₆₆~₀S₇₆十一个(图3(d))。其中，智利地震检测到₀S₇、₀S₁₄、₀S₁₅、₀S₁₈、₀S₂₄、₀S₂₅、₀S₄₃、₀S₅₅、₀S₇₃和₀S₇₅等十个和PREM模型本征周期一样，其余与理论值均有偏差，₀S₂、₀S₄、₀S₅、₀S₈、₀S₁₆、₀S₃₁、₀S₃₈和₀S₆₁等八个偏差较大，均大于0.6%，最大偏差为₀S₁₆的1.18% (表2)。日本地震检测到的₀S₆、₀S₁₂、₀S₁₆、₀S₂₁、₀S₂₄、₀S₄₆和₀S₅₀等七个和PREM模型本征周期一样，其余与理论值均有偏差，₀S₂、₀S₇、₀S₉、₀S₁₄、₀S₁₅、₀S₂₃、₀S₃₇和₀S₄₁等八个偏差较大，均大于0.5%，最大偏差为₀S₁₅的1.86% (表2)。

由于地球扁率和自转影响，地球某些本征振荡不是简并的，而是存在谱线分裂现象。汶川地震检测的₀S₂ (图2(a))、智利地震检测的₀S₆和₀S₇ (图2(b))等和日本地震检测的₀S₂和₀S₃ (图2(a))等均存在频谱分裂现象，由边际谱物理意义可知，振荡频谱分裂都是真实存在的。智利地震检测的₀S₉和₀S₁₀ (图2(b))等和日本地震检测的₀S₄和₀S₅ (图2(a))等也存在频谱分裂的可能。

5. 讨论与结论

本文首次采用基于希尔伯特–黄变换理论的边际谱方法，利用国内目前纬度最高的钻孔应变观测资料，提取了通河台钻孔应变2008年5月12日汶川M8.0级，2010年2月27日智利M8.8级和2011年3月11日日本本州东海岸附近海域M9.0级强震激发的地球自由振荡信息，认识如下：

1) 通河体应变准确检测到了2008年5月12日汶川M8.0级、2010年2月27日智利M8.8级和2011年3月11日的日本本州东海岸附近海域M9.0级地震激发的地球自由振荡信号，整个提取过程是在没有对数据进行去固体潮的前提下进行的。汶川地震检测到了24个基频信号，智利地震、日本地震检测到全部76个基频信号，包括₀S₀、₀S₂~₀S₇₆。汶川地震检测的两个信号、智利地震检测的十个信号、日本地震检测的七个信号和PREM模型本征周期一样，其余与理论值均有偏差。造成这种偏差的原因是不同震源机制激发的地球自由振荡具有不同的波谱结构和振荡特征。PREM模型仅考虑地球非均匀性均是随半径变化的，对横向不均匀性对振荡的影响研究还比较少有关。

2) 汶川地震检测的₀S₂、智利地震检测的₀S₆和₀S₇等和日本地震检测的₀S₂和₀S₃等信号均存在频谱分裂现象。究其原因，地球内部很多因素都会影响到地球自由振荡谱形态，如地球结构横向不均匀性及径向不均匀性、地球旋转、地球椭率、地球内部各向异性、物质衰减(与品质因子倒数有关)以及震源机制性质等因素。由于地球自转和椭率影响，地球自由振荡本征频率不是简并的，会出现地球自由振荡谱线(或谱峰)分裂现象。总之，这些因素会导致地球自由振荡谱分裂、耦合及谱峰变化等很多效应。这些振型谱值分裂为进一步研究地球内核的结构提供了较可靠的信息。

3) 本文应用边际谱方法提取地球自由振荡，与传统的傅里叶功率谱估计方法不同。边际谱从统计意

图3. 4000~8200 μHz自由震荡边际谱。(a) 4000~5000 μHz；(b) 5000~6000 μHz；(c) 6000~7000 μHz；(d) 7000~8200 μHz

表2. 球型基频自由振荡周期的观测值与PREM理论值对比

Continued

义上表征了整组数据每个频率点的累积幅值分布，而傅里叶频谱某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组份。边际谱可以处理非平稳信号，如果信号中存在某一频率能量出现，就表示一定有该频率的振动波出现，也就是说，边际谱能比较准确地反映信号的实际频率成分，而傅里叶变换只能处理平稳信号。

4) 希尔伯特变换只是单纯求解信号瞬时振幅，频率和相位，有可能出现没有意义的负频率；希尔伯特–黄变换先将信号进行EMD分解，得到各个不同尺度分量，然后对每一个分量进行希尔伯特变换得到有实际意义的瞬时频率。但是边际谱中还存在其他谱峰，这些谱峰是否是谐波振型和环形振型信息，还需进一步研究。

基金项目

中国地震局地震科技星火计划项目(XH16012)、中国地震局监测预报司震情跟踪青年课题(2016010132)共同资助。

文章引用

李继业,武晓军,马龙辰,任建辉,张彦吉,孙鹏宇,孙 强. 采用边际谱方法提取地球自由振荡
Using the Marginal Spectrum Method to Extract Earth Free Oscillation of Strong Earthquakes[J]. 应用物理, 2017, 07(06): 173-183. http://dx.doi.org/10.12677/APP.2017.76024

参考文献 (References)