2018, Vol. 61
2. Department of Mathematics, Duke University, Durham, NC, 27708, USA;
3. 中国地质大学(武汉)自动化学院, 武汉 430074;
4. 复杂系统先进控制与智能自动化湖北省重点实验室, 武汉 430074
2. Department of Mathematics, Duke University, Durham, NC, 27708, USA;
3. School of Automation, China University of Geosciences, Wuhan 430074, China;
4. Hubei Key Laboratory of Advanced Control and Intelligent Automation for Complex Systems, Wuhan 430074, China
地球天然脉冲电磁场(the Earth′s natural pulse electromagnetic field, ENPEMF)是指在地表能够接收到的由天然场源所产生的一次和二次综合电磁场总场.地震孕育和其他浅地表地球动力学过程中的震磁效应，可视为有效的天然场源(詹志佳等, 2000; 杨涛等, 2004; Malyshkov and Malyshkov, 2011).俄罗斯托姆斯克理工学院的Vorobyov教授(1970)提出地壳构造-电能之间的转换过程也是ENPEMF的信号场源之一.相比于大气雷暴所引发的一系列脉冲，地下的动力学脉冲源被称为“地下的雷暴”，尤其在滑裂类地质灾害来临前，其产生的脉冲波会发生剧烈异常(Sibgatulin et al., 2007; Malyshkov and Malyshkov, 2009)，可用于滑坡前兆电磁预警和震前电磁异常监测.
ENPEMF信号具有甚低频(VLF)非平稳脉冲特点，可通过设置阈值的VLF频段传感器在地表接收.上世纪80年代俄罗斯科学家Gokhberg等(1982)将该方法应用于如何准确预测地震震级、方位、深度等方面；乌克兰学者Boryssenko和Polishchuk(1999)使用该方法进行了浅层地表管道、线缆等其他隐藏物的检测研究；而俄罗斯科学院的Malyshkov和Malyshkov教授(2009)将此方法用于滑坡的风险预警和地下油气勘探.ENPEMF信号具有时变的振幅和频率，从时频分析的角度去分析和研究，可以挖掘信号时频部分隐藏的信息，有助于研究震前ENPEMF数据的时频异常特征(郝国成等, 2015; 郝国成等, 2016).时频分析是通过约束变换，把单一时域信号在“时间-频率-能量”的平面中展开，将时序信号变换为以时间和频率为变量的函数.相比于传统的傅里叶分析，时频分析更能突出信号的时间-频率的联合分布瞬态特征.
经典的时频分析方法有短时傅里叶变换(short time Fourier transform, STFT)(Portnoff, 1980)，小波变换(wavelet transform, WT)(Daubechies and Heil, 1992)，S变换(S transform, ST)(Stockwell et al., 1996)，Cohen型(Cohen, 1989)，Wigner-Ville分布(Wigner-Ville distribution, WVD)(Boashash and Black, 1987)等.这些经典的方法虽然有效，但都存在各自的缺点，如STFT的时频聚集度较差，WVD存在交叉项等.针对这些问题，一些改进算法，如STFT-WVD(王见等, 2013)，NSTFT-WVD(郝国成等, 2016)等和新兴的时频分析算法，如基于经验模态分解方法(Empirical mode decomposition, EMD)的希尔伯特黄变换(Hilbert-Huang transform, HHT)(Huang et al., 1998)、同步压缩小波变换(Synchrosqueezing wavelet transform, SWT)(Daubechies et al., 2011)、数据驱动时频分析(Data Driven Time-Frequency Analysis, DDTFA)方法(Hou and Shi, 2012)、布拉施克分解(Blaschke decomposition, BKD)(Coifman et al., 2017)、凸优化(Convex optimization, Tycoon) (Kowalski et al., 2016)等方法应运而生.时频分析方法已经广泛应用于医学、地球物理等多种技术领域，并取得了不错的成果.在地震勘探数据分析中，邓攻等(2015)将S变换应用于深反射地震弱信号的谱分解中，有效提高了信噪比和分辨率；杨锦玲等(2017)利用STFT方法有效观测到了于田MS7.3地震震前重力扰动异常信号的频率特征；吕君等(2012)在研究北京一次小级别地震时，使用WVD方法测得震前异常次声波的主要集中频率.本文则针对ENPEMF信号非平稳的特点，提出一种基于二值化(Binarized)同步压缩小波变换的数据驱动时频分析方法(BSWT-DDTFA).SWT算法可以提供聚集度较高的时间-频率-能量分布，但抗噪声性能较差；DDTFA方法基于EMD方法和压缩感知理论，通过自适应的稀疏分解提取信号的分量并可较好地重构信号，却过于依赖相位初始值的手动赋值，并且不能描述信号的能量分布.而BSWT-DDTFA是一种两者结合的优化算法，实现了初始相位较为准确的自动赋值，可以直接得到数据的时间-频率分布, 且保留了DDTFA抗噪声性能好的特性，适用于地球天然脉冲电磁场(ENPEMF)信号的分析研究.本文以2013年芦山MS7.0地震为例，利用BSWT-DDTFA算法对地震前后8天内采集的ENPEMF数据进行时频分析，研究ENPEMF信号的震前时频分布特点.
1 地球天然脉冲电磁场信号 1.1 ENPENF场源机理讨论ENPEMF是一种新兴的地球物理应用领域的技术方法，其场源机理有多种假设.俄罗斯学者Aleksandrov等(1972)将其归因于大气雷暴，俄罗斯Malyshkov和Malyshkov教授(2009)认为其与地震密切相关.地表可获得的脉冲可映射地球内部某些活动的发生和发展，并且这些数据具有日变化和年变化的规律，携带了大量有价值的电磁信息.与此同时，随着地球物理学及其他相关学科的研究和发展，震磁机理也出现了更进一步的实验和解释.比如震源区应力的逐渐积累和重新分布致使岩石构造磁性发生的变化、地下介质发生的形变与相变、应变波传播激发电磁效应等现象，会产生压磁效应、感应磁效应、电动磁效应及热磁效应、地震-电离层耦合等相关活动(Kuo et al., 2014; Pulinets et al., 2015)，影响到地表瞬变磁场的接收，地表观测设备可以捕获到宽频谱的电磁辐射波.震前电磁效应在地震前兆信息研究中一直是热点问题，各国学者从多个角度出发研究电磁孕震信息，试图建立地震与电磁之间的潜在联系，研究地震前电磁信号的异常特点.基于微破裂机械能-电能转换机制和“地壳波导”观点，各类滑裂型地质灾害前兆及其伴生现象(力学、物理、化学等各类异常)可在地表产生甚低频(VLF)的ENPEMF信号脉冲波动.ENPEMF信号携带了大量有价值的电磁异常信息，反映地质活动对较大尺度时空的影响，本文采用合适的非平稳信号时频分析方法, 研究ENPEMF信号在强震灾害发生期间的时频异常特点，为进一步研究ENPEMF方法提供技术帮助.
1.2 ENPENF接收设备ENPEMF信号接收设备(安放在武汉九峰地震台)采用的频段最高延展至25 kHz，设置的接收频率为14.5 kHz.相对正常背景场，地震事件来临前，ENPEMF信号会有明显的异常起伏变化.即在此观测频率，有可能发现一种新的与孕震信息有关的变化规律.
接收ENPEMF信号的仪器为俄罗斯科学院托木斯克分院GR-01型设备，其接收工作原理如图 1所示.
|
图 1 ENEMF设备的工作原理 Fig. 1 Working principle of the ENPEMF signal receiving device |
图 1所示仪器记录了ENPEMF信号超过设定阈值的脉冲数目和幅度，每天6个文件，每4 h利用GPRS上传一次数据，可以用互联网远程下载.如果仪器接收到了大规模异常变化的脉冲，与正常的背景脉冲数目包络轨迹不同，则有可能包含某地异常变化的孕震信息.
1.3 ENPEMF信号武汉九峰地震台安放的3台GR-01接收设备摆放的方向为西-东(W-E)和北-南(N-S)，分为3个通道：CN1、CN2和CN3.在安装调试设备时，根据波形的包络大小和曲线形状，参照指导书校对波形，设置合适的放大阈值，得到后缀为.gr1的文件，再通过专用程序转换成.txt文件. ENPEMF数据具有以下几个特点：①设备传感器设置为甚低频段：5~25 kHz，在武汉设置的接收频率为14.5 kHz；②数据经过转换处理后，采样率为1 s；③数据存储格式为：时间-幅度-脉冲数(t-AH-NH)，时间单位为秒(s). AH为超过设定阈值的量化后的脉冲幅度，单位是原始毫伏(mV)信号程控放大后的相对数值，在这里只是一个包络大小变化的参照量；NH为超过设定阈值的量化后的脉冲数目，两者均表征地表磁场的强弱；④ENPEMF信号属于为非平稳信号，频域特点明显，对其进行相应的频率分析，了解震前信号的时频分布特点，有助于识别孕震信息.
2013年4月20日，我国四川省雅安市芦山县发生7.0级大地震，其地理位置及周边地区活动断层带如图 2所示，曾祥方(2013)和Han等(2014)研究了本次芦山地震的成因机制.本文从时频分析的角度研究芦山地震前后(15日到22日)ENPEMF信号的时频分布特点.数据来源为武汉九峰地震台GR-01设备(通道CN2)，图 3a是原始数据的包络图，可以看出，ENPEMF信号的幅值密且乱，无法从中获得有价值的分析信息.如果使用FFT全局变换的频域分析方法，只能呈现原始信号中存在哪些频率分量及相应的幅值，如图 3b所示，不能提供细节信息，还需要用时频分析方法进一步研究局部的时间-频率-能量谱的分布特点.
|
图 2 2013年4月20日芦山MS7.0地震所处地震带及周边地区活动断层图分布图(王杰等, 2013) Fig. 2 The seismic zone and distribution of surrounding active faults of Lushan MS7.0 earthquake on April 20th |
|
图 3 (a) ENPEMF原始数据的包络图(通道CN2)；(b) 4月20日数据的傅里叶变换(FFT)频谱图 Fig. 3 (a) Envelope of original ENPEMF data (channel 2); (b) Fourier transform (FFT) spectrum of the data on April 20th |
DDTFA方法对信号是完全自适应的，相对于EMD和EEMD方法而言，不仅减弱了端点效应和模态混叠现象，抗噪声性能强，而且有更加坚实的数学理论支撑(Hou and Shi, 2012). DDTFA有两个比较重要的方面，以保证对数据的完全自适应性：第一是用来分解数据的基是来源于数据本身，而不是之前设定的；第二是在包含本征模态函数的字典中寻找信号的最稀疏表示.最稀疏表示是指分解结果在分量具有物理意义的基础上保证分量个数最少.所以DDTFA的处理过程可以分为两个部分：
(1) 首先建立一个高度冗余的字典D，其定义式为：
|
(1) |
其中，V(θ, λ)表示所有比cosθ(t)平滑的函数的集合，以一组过完备傅里叶基构造的V(θ, λ)可表示为
|
(2) |
其中，
(2) 然后，通过合适的稀疏分解方法寻找信号在这个字典上的最稀疏表示，该过程可以通过求解一个非线性优化问题P0来解决.
P0：最小化M，使其满足条件式(3)：
|
(3) |
其中，fk被称为第k个本征模态函数(IMF)，ak和θk分别是第k个IMF的幅度和相位.
在DDTFA算法中，稀疏分解的过程可选择基于非线性三阶全变差(TV3)最小化的分解或者基于非线性匹配追踪法的分解两种方法.但是由于非线性TV3最小化的分解对噪声比较敏感，而非线性匹配追踪法的计算量较大，所以本文选择针对周期数据的一种基于快速傅里叶变换(FFT)的快速算法(Hou and Shi, 2012; Hou et al., 2013; 张学阳, 2012)，实验证明此快速算法对于非周期数据也具有一定的适用性.
2.2 同步压缩小波变换(SWT)同步压缩方法可以看作是一种特殊的时频重排技术，将时频能量谱中的每一个点根据一定的映射关系对应到一个新的时间和尺度/频率坐标，也就是按照一定的重排算子向其局部重心进行重新分配(Auger et al., 2013; Daubechies et al., 2011).SWT是基于连续小波变换的同步压缩算法，即对连续小波变换后的小波系数进行压缩和重新分配，目的是使信号的频率在时频图中的位置更加准确，分布更加清晰(Thakur et al., 2013; Herrera et al., 2014; Wu et al., 2013)，其主要实现步骤可以简单归纳为：
(1) 对于待分析信号s(t)，选择一个适当的母小波ψ(t)，计算其连续小波变换Ws(a, b)：
|
(4) |
其中，a是尺度因子，b是时间平移因子；
(2) 由式(5)计算时频重排算子ωs(a, b)：
|
(5) |
其中，
(3) 依据重排算子ωs(a, b)对连续小波变换的结果Ws(a, b)进行同步压缩，得到压缩后的小波系数谱Ts(ω, b)：
|
(6) |
其中，A(b)={a; Ws(a, b)≠0}, a∈A(b).算法流程如图 4所示.
|
图 4 SWT算法流程图 Fig. 4 The flow chart of SWT |
对于DDTFA方法，当非线性匹配追踪分解过程中的初始相位赋值为信号分量的平均频率时，得到的IMF分量和瞬时频率曲线比较准确.目前给DDTFA初始相位赋值的方式有两种：一种是将信号进行短时傅里叶变换，然后通过肉眼识别信号的主要频率作为初始相位，这种方式过程比较繁琐，并且会受主观因素的影响(Hou and Shi, 2012)；另外一种是通过EMD得到信号的IMF分量，然后以每个IMF分量的平均频率作为DDTFA的初始赋值(张学阳, 2012)，这种EMD-DDTFA算法虽然可以实现自适应初始赋值，但是由于EMD本身存在模态混叠现象，并且容易受噪声的影响，会导致得到初始频率赋值不准确，时频结果错误.本文提出的改进算法(BSWT-DDTFA)是通过二值化SWT的处理结果提取信号的主要频率，实现DDTFA相位初值的自动赋值.BSWT-DDTFA方法既能得到聚集度高的时频分布，又保留了DDTFA抗噪声性能好的优点.
BSWT-DDTFA方法的技术路线如图 5所示，主要包含如下五个步骤：
|
图 5 BSWT-DDTFA方法的流程图 Fig. 5 The flow chart of BSWT-DDTFA |
(1) 通过SWT得到信号x(t)的时频分布；
(2) 对SWT的时频分布进行二值化，提取信号的主要频率值；
(3) 将提取出来的频率值乘以2πt作为数据驱动时频分析方法的初始相位值θ0；
(4) 对信号x(t)在参数θ0下进行DDTFA变换，得到IMF分量和瞬时频率曲线；
(5) 将得到的IMF分量分别进行同步压缩小波变换，叠加后得到x(t)的时频分布.
3 实验仿真 3.1 BSWT-DDTFA方法的有效性BSWT-DDTFA方法能有效地分解信号，得到IMF分量和频率曲线，可以直接绘制信号的时间-频率分布图，且抗噪声性能比较优越.本节设计了两组不同的仿真实验，将仿真信号分别采用BSWT-DDTFA算法和EMD-DDTFA算法分解，其中SWT选用morlet母小波.实验结果证明了BSWT-DDTFA算法得到的IMF分量和瞬时频率曲线都优于EMD-DDTFA算法.
仿真信号x1(t)：三分量调频信号，时间区间t∈[0, 4]，采样点数为2048，信噪比SNR=1.5 dB，信号函数式为
|
(7) |
仿真信号x2(t)：二分量间断调频信号，时间区间t∈[0, 2]，采样点数为2048，信噪比SNR=1 dB，信号函数式为
|
(8) |
图 6是仿真信号x1(t)和x2(t)在含噪声情况下分别通过EMD-DDTFA和BSWT-DDTFA方法分解得到的IMF分量.对于x1(t)，由EMD-DDTFA分解的3个IMF分量几乎都是无规则的噪声分量，与信号的实际分量相似度极小，几乎不能分辨；而BSWT-DDTFA分解得到的IMF分量与实际分量相比，虽然幅度略有偏差，但每个分量的频率都很接近.对于x2(t)，EMD-DDTFA同样没有正确分解出信号的分量，而BSWT-DDTFA可以将信号按频率分解，尤其是能够将IMF1在t=1 s之后才有幅度变化的特点很好地表现出来.以上结果表明EMD-DDTFA受信号中噪声的影响很大，而BSWT-DDTFA除同样能实现初始相位的自动赋值外，还可以更为准确地分解信号.
|
图 6 不同方法分解得到的信号x1(t)和x2(t)的分量 (a)信号函数式定义的实际分量；(b) EMD-DDTFA方法得到的信号的IMF分量；(c) BSWT-DDTFA方法得到的信号的IMF分量. Fig. 6 IMFs of x1(t) and x2(t) decomposed by different methods (a) The actual components defined by the signal function; (b) IMFs decomposed by EMD-DDTFA; (c) IMFs decomposed by BSWT-DDTFA. |
如图 7所示，比较两个仿真信号由BSWT-DDTFA、EMD-DDTFA方法得到的瞬时频率曲线.由图 7a可以看出，EMD-DDTFA方法绘制的瞬时频率曲线中含有大量噪声分量，基本不能分辨信号的真实频率(虚线).噪声分量的频率大致分布在50~300 Hz，远远高于信号的真实频率，这是由于EMD在信号含噪声情况下分解不准确，计算得到的相位初值错误所导致的.相比之下，BSWT-DDTFA的频率曲线(图 7b)有很明显的改善；除在x1(t)的10 Hz分量的两个端点处存在端点飞翼外，其他与信号的真实频率曲线(虚线)几乎重合，表明与EMD-DDTFA相比，BSWT-DDTFA方法可以更准确地获得信号分量的初始相位和瞬时频率曲线.
|
图 7 仿真信号x1(t)和x2(t)的真实瞬时频率(虚线)和(a)EMD-DDTFA瞬时频率(实线); (b) BSWT-DDTFA瞬时频率(实线) Fig. 7 The actual instantaneous frequency (dashed) and (a) EMD-DDTFA instantaneous frequency (solid); (b) BSWT-DDTFA instantaneous frequency (solid) of signal x1(t) and x2(t) |
再比较SWT和BSWT-DDTFA得到的时频分布，如图 8所示.在有噪声干扰下，对于仿真信号x1(t)，SWT时频分布图中有少量噪声信号的分布，每个分量都有频率扩散的现象，使得时频分布模糊；而BSWT-DDTFA时频分布图中的噪声得到了抑制，尤其是10 Hz和5 Hz的分量更为平滑，频率聚集度更高.对于仿真信号x2(t)，频率扩散的现象同样存在，且在5 Hz以下存在较为明显的干扰分量；而BSWT-DDTFA时频分布中两个分量都更为平滑，分辨效果更好.表明BSWT-DDTFA比SWT的抗噪声性能好，在含噪声情况下可以得到较高聚集度的时频分布图.
|
图 8 不同方法得到的仿真信号x1(t)和x2(t)的时频分布图 (a) SWT时频分布；(b) BSWT-DDTFA时频分布 Fig. 8 The time-frequency distribution of signal x1(t) and x2(t) (a) The time-frequency distribution from SWT; (b) The time-frequency distribution from BSWT-DDTFA |
实际工程应用中采集的信号均为含噪声数据，所以下面通过BSWT-DDTFA方法处理不同信噪比的信号，更为直观地研究BSWT-DDTFA方法的抗噪声性能.
对于信号稀疏分解的效果，为避免肉眼观察带来的主观性和局限性，本文引入了一个参数来评价稀疏分解的效果.计算BSWT-DDTFA算法的IMF分量的绝对误差平均值，并将其作为判断BSWT-DDTFA处理结果优劣的一个性能指标，用Ei表示，其定义式为(9)式：
|
(9) |
其中，IMF(:, i)表示通过BSWT-DDTFA方法得到的IMF分量，rIMF(:, i)表示信号的实际分量，N代表采样点数. Ei值越小，表示误差越小，处理结果越好.
仿真信号sig(t)的时间区间为t∈[0, 2]，采样点数为2048，函数表达式为
|
(10) |
对仿真信号sig(t)添加不同信噪比的噪声，由BSWT-DDTFA分解得到的分量如图 9所示.在信噪比为20 dB、10 dB和0 dB时，BSWT-DDTFA可以将信号的三个分量近乎理想地分解出来；信噪比为-10 dB时，IMF2的2 s和IMF3的0.5 s处幅度出现较大误差，受到噪声的一点影响，但是整体仍能比较准确地得到信号分量.为了量化分析，表 1中给出仿真信号sig(t)由EMD-DDTFA和BSWT-DDTFA方法得到的3个IMF分量在不同信噪比之下的平均误差值，并绘制成折线图，如图 10所示.
|
图 9 不同信噪比下sig(t)由BSWT-DDTFA分解的分量 (a) SNR=20 dB; (b) SNR=10 dB; (c) SNR=0 dB; (d) SNR= -10 dB. Fig. 9 The components of sig(t) decomposed by BSWT-DDTFA with different SNR |
|
表 1 通过不同方法得到的IMF分量的绝对误差平均值 Table 1 The absolute error of the IMFs obtained by the BSWT-DDTFA method |
|
图 10 EMD-DDTFA和BSWT-DDTFA方法得到的IMF分量的绝对误差平均值 Fig. 10 The average absolute error of the IMFs decomposed by EMD-DDTFA and BSWT-DDTFA |
结合表 1中的数值和图 10中的折线图，在实验的信噪比范围内，BSWT-DDTFA的误差都小于EMD-DDTFA方法的误差；尤其是在信噪比小于15 dB时，EMD-DDTFA的误差明显增大，达到0.35以上，远远大于BSWT-DDTFA的误差.并且EMD-DDTFA的误差波动范围较大，而相比之下BSWT-DDTFA方法的误差波动范围较小.所以BSWT-DDTFA方法对噪声抑制效果较好，抗噪声性能优越.
通过图 11，比较不同程度的噪声对SWT和BSWT-DDTFA时频分布的影响.可以看出在信噪比从20 dB降低到0 dB的过程中，SWT时频分布在信号主要频率附近的模糊现象变得严重，而BSWT-DDTFA时频分布一直保持较高的聚集度，尤其是20 Hz分量，相对于SWT更为平滑；当信噪比达到-10 dB时，SWT的时频分布发生了严重的畸变，但在BSWT-DDTFA的时频-能量分布中，只有最高频率的频率分量受到了较为明显的影响，而其他两个频率分量的描述仍然比较清楚.所以图 11的实验结果充分说明了BSWT-DDTFA可以有效地抑制噪声对信号分析研究的影响.
|
图 11 不同信噪比下sig(t)由SWT和BSWT-DDTFA绘制的时频分布图 (a) SNR=20 dB; (b) SNR=10 dB; (c) SNR=0 dB; (d) SNR= -10 dB. Fig. 11 The time-frequency distribution of sig(t) plotted by SWT and BSWT-DDTFA with different SNR |
从以上仿真信号验证得到的结果图可以得出，BSWT-DDTFA算法通过BSWT来提取信号的主要频率，最终实现了数据驱动时频分析方法自适应初始赋值，并且可以得到与真实信号更加接近的IMF分量与时频分布，比EMD-DDTFA算法得到的结果更加准确，比EMD-DDTFA算法和SWT算法的抗噪声性能更好.
4 基于BSWT-DDTFA的震前ENPEMF数据的时频特点分析地球天然脉冲电磁场信号是一种典型的非平稳信号，而且由于它的场源比较复杂，含有较多的噪声，所以将BSWT-DDTFA方法应用于地球天然脉冲电磁场信号是比较合适的.将本文所选取如表 2的AH数据，进行BSWT-DDTFA处理分析，可以得到如图 12所示的结果，希望能从中了解震前ENPEMF信号的时频分布特点.
|
|
表 2 通道CN2的观测数据 Table 2 Observation data of the S-N Channel (CN2) |
|
图 12 2013年芦山地震前后ENPEMF信号的时频分布的变化 (a) 4月15日; (b) 4月16日; (c) 4月17日; (d) 4月18日; (e) 4月19日; (f) 4月20日(地震发生); (g) 4月21日; (h) 4月22日. Fig. 12 The variation in time-frequency distribution of ENPEMF data before and after the Lushan earthquake in 2013 (a) April 15th; (b) April 16th; (c) April 17th; (d) April 18th; (e) April 19th; (f) April 20th (the earthquake happened); (g) April 21st; (h) April 22nd. |
图 12是处理的芦山地震期间2013年4月15日到2013年4月22日之间的ENPEMF信号，地震发生的时间是2013年4月20日，通过BSWT-DDTFA和SWT两种方法得到了信号的时间-频率-能量分布图，此部分主要观察震前时频分布的异常特征.从图中可以看出BSWT-DDTFA比SWT更有效地屏蔽噪声的影响，将把主要的频率显现出来，更为清楚地分辨出时频的变化规律.ENPEMF信号的频率主要集中在0.1 Hz以下，在震前几天，时频分布图的0.1~0.3 Hz范围内开始出现频率分量，并逐天增多，在震前两天又恢复平静，而在震前一天再次增多并达到最大值，震后逐渐恢复.另外，ENPEMF信号的能量分布在地震前后的几天内很不稳定，会出现明显的增强或者减弱.
通过分析BSWT-DDTFA方法处理得到的ENPEMF数据的时频-能量分布图，可以发现，在地震发生前后期间，频率分量会经过增多-减少-增多-平稳的过程，能量也经过增大-减少-增大-平稳的过程.在地震之后(21日)开始逐步减少，22日的时频平面则恢复平静，频率分量及其噪声分布没有明显表现.图 12的芦山地震前后时频分布图中，频率分量的幅度/能量谱也随之变化，一定程度上对应了地震孕震区域的地质活动发展趋势.
5 结论本文针对ENPEMF信号的非平稳特点，用改进后的时频分析方法研究芦山地震前的时间-频率联合分布特性.为提高频率表示的精度和抗噪性能，解决数据驱动时频分析(DDTFA)方法的初值设置问题，提出了基于二值化同步压缩小波变换的改进算法(BSWT-DDTFA). BSWT-DDTFA算法具备对输入时序信号的自适应特点，可直接得到信号的IMF分量、频率曲线及时频分布，该算法的时间和频率聚集度、抗噪性能都得到了较好的提升.加噪的数据仿真实验结果表明，在-10 dB至20 dB的白噪声干扰环境下，对比基于经验模态分解的数据驱动时频分析(EMD-DDTFA)方法，本文的BSWT-DDTFA方法可分解得到误差更小的IMF分量，抗噪声性能也得到了较为明显的改善.同时，相较于同步压缩小波变换，BSWT-DDTFA得到的时频分布更为清晰，能够更加准确地反映非平稳信号的频率随时间的分布规律和特征.本文应用BSWT-DDTFA方法对芦山地震前后8天的ENPEMF数据进行时频分析，处理结果显示：震前ENPEMF信号的频谱宽度会出现展宽-缩窄-展宽的趋势，信号能量谱会出现增强-减弱-增强-临震的特点.
致谢 感谢审稿专家提出的宝贵意见.感谢合作导师数学系Ingrid Dubeches对本论文的方法指导.感谢杜克大学Hau-Tieng Wu和清华大学史作强提供的技术指导.感谢杜克大学iiD中心(the information initiative at Duke)、中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室对本研究提供的学术支持.感谢武汉九峰地震台对设备的精心维护.
Aleksandrov M S, Baklenova Z M, Gladshtein N D. 1972. Fluktuatsii elektromagnitnogo polya Zemli v diapazone SNCh (ULF Fluctuations of the Earth's Electromagnetic Field). Moscow: Nauka.
|
Auger F, Flandrin P, Lin Y T, et al. 2013. Time-frequency reassignment and Synchrosqueezing:An overview. IEEE Signal Processing Magazine, 30(6): 32-41. DOI:10.1109/MSP.2013.2265316 |
Boashash B, Black P. 1987. An efficient real-time implementation of the Wigner-Ville distribution. IEEE Transactions on Acoustics, Speech, and Signal Processing, 35(11): 1611-1618. DOI:10.1109/TASSP.1987.1165070 |
Boryssenko A, Polishchuk V I. 1999. Earth near-surface passive probing by natural pulsed electromagnetic field.//Proceedings (ICCEA'99) 1999 International Conference on Computational Electromagnetics and Its Applications. Beijing, China: IEEE, 529-532, doi: 10.1109/ICCEA.1999.825235.
|
Cohen L. 1989. Time-frequency distribution-a review. Proceedings of the IEEE, 77(7): 941-981. DOI:10.1109/5.30749 |
Coifman R R, Steinerberger S, Wu H T. 2017. Carrier frequencies, holomorphy, and unwinding. SIAM Journal on Mathematical Analysis, 49(6): 4838-4864. DOI:10.1137/16M1081087 |
Daubechies I, Heil C. 1992. Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics.
|
Daubechies I, Lu J F, Wu H T. 2011. Synchrosqueezed wavelet transforms:An empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis, 30(2): 243-261. DOI:10.1016/j.acha.2010.08.002 |
Deng G, Liang F, Li X T, et al. 2015. S-transform spectrum decomposition technique in the application of the extraction of weak seismic signals. Chinese Journal of Geophysics (in Chinese), 58(12): 4594-4604. DOI:10.6038/cjg20151221 |
Gokhberg M B, Morgounov V A, Yoshino T, et al. 1982. Experimental measurement of electromagnetic emissions possibly related to earthquakes in Japan. Journal of Geophysical Research:Solid Earth, 87(B9): 7824-7828. DOI:10.1029/JB087iB09p07824 |
Han L B, Zeng X F, Jiang C S, et al. 2014. Focal mechanisms of the 2013 MW6.6 Lushan, China earthquake and high-resolution aftershock relocations. Seismological Research Letters, 85(1): 8-14. DOI:10.1785/0220130083 |
Hao G C, Gong T, Dong H B, et al. 2015. Time-frequency characteristics and energy analysis of the Earth's natural pulse electromagnetic fields based on ensemble empirical mode decomposition:The Lushan MS7.0 earthquake as an example. Earth Science Frontiers (in Chinese), 22(5): 231-238. DOI:10.13745/j.esf.2015.05.019 |
Hao G C, Chen Z C, Zhao J, et al. 2016. Time-frequency analysis of the Earth's natural pulse electromagnetic field signal before and after the Lushan MS7.0 earthquake based on NSTFT-WVD transform. Earth Science Frontiers (in Chinese), 23(1): 276-286. DOI:10.13745/j.esf.2016.01.025 |
Herrera R H, Han J J, Van Der Baan M. 2014. Applications of the synchrosqueezing transform in seismic time-frequency analysis. Geophysics, 79(3): V55-V64. DOI:10.1190/geo2013-0204.1 |
Hou T Y, Shi Z Q. 2012. Data-driven time-frequency analysis. Applied and Computational Harmonic Analysis, 35(2): 284-308. DOI:10.1016/j.acha.2012.10.001 |
Hou T Y, Shi Z Q, Tavallali P. 2013. Convergence of a data-driven time-frequency analysis method. Applied and Computational Harmonic Analysis, 37(2): 235-270. DOI:10.1016/j.acha.2013.12.004 |
Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. The Royal Society, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193 |
Kowalski M, Meynard A, Wu H T. 2016. Convex Optimization approach to signals with fast varying instantaneous frequency. Applied and Computational Harmonic Analysis, 44(1): 89-122. DOI:10.1016/j.acha.2016.03.008 |
Kuo C L, Lee L C, Huba J D. 2014. An improved coupling model for the lithosphere-atmosphere-ionosphere system. Journal of Geophysical Research:Space Physics, 119(4): 3189-3205. DOI:10.1002/2013JA019392 |
Lv J, Guo Q, Feng H N, et al. 2012. Anomalous infrasonic waves before an small earthquake in Beijing. Chinese Journal of Geophysics (in Chinese), 55(10): 3379-3385. DOI:10.6038/j.issn.0001-5733.2012.10.020 |
Malyshkov Y P, Malyshkov S Y. 2009. Periodicity of geophysical fields and seismicity:Possible links with core motion. Russian Geology and Geophysics, 50(2): 115-130. DOI:10.1016/j.rgg.2008.06.019 |
Malyshkov Y P, Malyshkov S Y. 2011. Eccentric rotation of the earth's core and lithosphere: Origin of deformation waves and their practical application. The Earth's Core: Structure, Properties and Dynamics, chapter 6, 1-95.
|
Portnoff M. 1980. Time-frequency representation of digital signals and systems based on short-time fourier analysis. IEEE Transactions on Acoustics, Speech, and Signal Processing, 28(1): 55-69. DOI:10.1109/TASSP.1980.1163359 |
Pulinets S A, Ouzounov D P, Karelin A V, et al. 2015. Physical bases of the generation of short-term earthquake precursors:A complex model of ionization-induced geophysical processes in the lithosphere-atmosphere-ionosphere-magnetosphere system. Geomagnetism and Aeronomy, 55(4): 521-538. DOI:10.1134/S0016793215040131 |
Sibgatulin V G, Peretokin S A, Khlebopros R G. 2007. Fundamental peculiarities of the entropy model of energy processes in seismic areas. Earth Science Frontiers, 14(6): 222-225. DOI:10.1016/S1872-5791(08)60013-5 |
Stockwell R G, Mansinha L, Lowe R P. 1996. Localization of the complex spectrum:The S transform. IEEE Transactions on Signal Processing, 44(4): 998-1001. DOI:10.1109/78.492555 |
Thakur G, Brevdo E, Fučkar N S, et al. 2013. The synchrosqueezing algorithm for time-varying spectral analysis:Robustness properties and new paleoclimate applications. Signal Processing, 93(5): 1079-1094. DOI:10.1016/j.sigpro.2012.11.029 |
Vorobyov A A. 1970. On probability of electrical discharges in the Earth's interior. Geological I Geophysics, (12): 3-13. |
Wang J, Li J T, Lu H L, et al. 2013. Using STFT-Wigner transform to suppress the cross terms in Wiger-Ville distribution. Journal of Chongqing University (in Chinese), 36(8): 15-18, 25. DOI:10.11835/j.jssn.1000-582X.2013.08.003 |
Wu H T, Chan Y H, Lin Y T, et al. 2013. Using synchrosqueezing transform to discover breathing dynamics from ECG signals. Applied and Computational Harmonic Analysis, 36(2): 354-359. DOI:10.1016/j.acha.2013.07.003 |
Yang J L, Li Z N, Guan Y M, et al. 2017. Study on gravity disturbance before the Yutian MS7.3 earthquakes. Chinese Journal of Geophysics (in Chinese), 60(10): 3844-3852. DOI:10.6038/cjg20171014 |
Zeng X F, Luo Y, Han L B, et al. 2013. The Lushan MS7. 0 earthquake on 20 April 2013:A high-angle thrust event. Chinese Journal of Geophysics (in Chinese), 56(4): 1418-1424. DOI:10.6038/cjg20130437 |
Zhan Z J, Gao J T, Zhao C L, et al. 2000. Research on tectonomagnetism and its application to earthquake prediction. Earthquake (in Chinese), 20(S1): 127-134. DOI:10.3969/j.issn.1000-3274.2000.z1.022 |
Zhang X Y. 2012. Improved data-driven time-frequency analysis and its application[Master's thesis] (in Chinese). Changsha: National University of Defense Technology.
|
邓攻, 梁锋, 李晓婷, 等. 2015. S变换谱分解技术在深反射地震弱信号提取中的应用. 地球物理学报, 58(12): 4594-4604. DOI:10.6038/cjg20151221 |
郝国成, 龚婷, 董浩斌, 等. 2015. 基于聚类经验模态分解的地球天然脉冲电磁场时频与能量谱分析:以芦山MS7.0地震为例. 地学前缘, 22(5): 231-238. DOI:10.13745/j.esf.2015.05.019 |
郝国成, 陈忠昌, 赵娟, 等. 2016. 基于NSTFT-WVD变换的芦山MS7.0级地震前后地球天然脉冲电磁场信号时频分析. 地学前缘, 23(1): 276-286. DOI:10.13745/j.esf.2016.01.025 |
吕君, 郭泉, 冯浩楠, 等. 2012. 北京地震前的异常次声波. 地球物理学报, 55(10): 3379-3385. DOI:10.6038/j.issn.0001-5733.2012.10.020 |
王见, 李金同, 卢华玲, 等. 2013. 采用STFT-Wigner变换抑制Wigner-Ville分布交叉项. 重庆大学学报, 36(8): 15-18, 25. DOI:10.11835/j.jssn.1000-582X.2013.08.003 |
王杰, 张雄, 潘黎黎, 等. 2013. 芦山地震(MS7.0)前甲烷释放与大气增温异常. 地学前缘, 20(6): 29-35. |
杨锦玲, 李祖宁, 关玉梅, 等. 2017. 于田MS7.3地震震前重力扰动信号研究. 地球物理学报, 60(10): 3844-3852. DOI:10.6038/cjg20171014 |
杨涛, 刘庆生, 付媛媛, 等. 2004. 震磁效应研究及进展. 地震地磁观测与研究, 25(6): 63-71. DOI:10.3969/j.issn.1003-3246.2004.06.010 |
曾祥方, 罗艳, 韩立波, 等. 2013. 2013年4月20日四川芦山MS7.0地震:一个高角度逆冲地震. 地球物理学报, 56(4): 1418-1424. DOI:10.6038/cjg20130437 |
詹志佳, 高金田, 赵从利, 等. 2000. 构造磁学及其预测地震研究. 地震, 20(S1): 126-134. DOI:10.3969/j.issn.1000-3274.2000.z1.022 |
张学阳. 2012.改进的数据驱动时频分析方法及其应用[硕士论文].长沙: 国防科学技术大学. http://cdmd.cnki.com.cn/Article/CDMD-90002-1014048145.htm
|