地球物理学报  2017, Vol. 60 Issue (5): 1918-1936   PDF    
基于递归分析和聚类的大地电磁信噪辨识及分离
李晋1,2, 汤井田1 , 燕欢2, 彭代鑫2, 徐志敏1,3     
1. 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 地球科学与信息物理学院, 长沙 410083;
2. 湖南师范大学物理与信息科学学院, 长沙 410081;
3. 西北综合勘察设计研究院, 西安 710003
摘要: 为了剖析大地电磁信号和强干扰的本质特征,进一步精细分离出微弱的大地电磁有用信号,提出基于递归分析和聚类的大地电磁信噪辨识及分离方法.首先,运用递归分析法扩展大地电磁一维时间序列的维数,分析了嵌入维数、延迟时间和判别阈值对递归图的性能,并研究了不同长度的序列对递归定量分析参数的影响情况.然后,构建典型的大地电磁强干扰类型和微弱的大地电磁有用信号样本库,针对样本库讨论了强干扰和微弱大地电磁信号之间的递归定量分析参数,分析了K均值聚类和模糊C均值聚类的信噪辨识效果.最后,对实测大地电磁数据进行信噪辨识处理,并仅对辨识为强干扰的时间段采用数学形态滤波进行噪声压制.实验结果表明,递归分析能定性及定量地描述大地电磁信号时间序列的非线性特征和原动力系统的本质规律,与聚类算法相结合能对矿集区实测大地电磁信号进行信噪辨识;处理后的卡尼亚电阻率-相位曲线更为光滑、连续,其结果更为精细地保留了大地电磁信号低频段的缓变化信息,整个低频段的大地电磁数据质量得到了明显改善.
关键词: 大地电磁      信噪辨识      信噪分离      递归分析      聚类     
Identification and spearation of magnetotelluric signal and noise based on recurrence analysis and clustering
LI Jin1,2, TANG Jing-Tian1, YAN Huan2, PENG Dai-Xin2, XU Zhi-Min1,3     
1. School of Geosciences and Info-Physics, Key Laboratory of Metallogenic Prediction of Non-Ferrous Metals and Geological Environment Monitor, Ministry of Education, Central South University, Changsha 410083, China;
2. Institute of Physics and Information Science, Hunan Normal University, Changsha 410081, China;
3. Northwest Research Institute of Engineering Investigations and Design, Xi'an 710003, China
Abstract: In order to analyze the essential characteristics of magnetotelluric (MT) sounding data and strong interferences, and further to separate weak useful MT signal, we propose a new method for identification and separation of MT signal and noise based on recurrence analysis and clustering. First, we use the recurrence analysis method to extend the dimension of MT time series and analyze the embedding dimension, delay time and determination threshold, and study the Recurrence Quantification Analysis (RQA) parameters for the sequence of different lengths. Then, we build a sample dataset for typical strong interferences and weak useful MT signal. According to the sample dataset, we discuss the RQA parameters between strong interferences and weak MT signal, and analyze the effects of signal and noise identification for the K-means clustering and fuzzy C-means clustering. Finally, the real MT data is processed through signal and noise identification, and only the time series containing strong interferences are suppressed by mathematical morphology filtering. Experimental results show that the recurrence analysis can qualitatively and quantitatively describe the nonlinear characteristics of the time series and the essential rule of the prime power system for MT. Combined with clustering algorithm permits to identify the measured MT data in ore concentration area. The Cagniard resistivity-phase curve is more smooth and continuous after using the proposed method. Moreover, the slow change information of low frequency for MT is more retained finely, and the quality of MT data for the overall low frequency band is improved significantly.
Key words: Magnetotelluric sounding data      Signal and noise identification      Signal and noise separation      Recurrence analysis      Clustering     
1 引言

大地电磁测深法 (Magnetotelluric,MT) 是基于电磁感应原理,利用天然交变电磁场来研究地下岩层的电学性质及其分布特征 (Tikhonov,1950).由于天然大地电磁场频带范围宽、信号微弱,实际观测到的大地电磁信号极易受到各种噪声的污染 (魏文博,2002).在我国“深部探测技术与实验研究”(SinoProbe) 专项实施过程中,不可避免地需要在长江中下游成矿带及典型矿集区开展大地电磁探测工作 (吕庆田等, 2011, 2015董树文等,2012).这些地区人烟稠密、经济发达,随处可见的高压电网、广播电台、通讯电缆、信号发射塔、各种金属管网,以及用于矿山开采的大功率直流电机车等众多不利因数制约了大地电磁数据本身的可靠性.大地电磁数据质量急剧下降,阻抗估算偏差严重及测量获得的视电阻率值过度失真等状况导致不能客观反映地下电性分布,严重影响了后续的电磁法反演解释水平,这些都给地球物理勘探工作带来了巨大困扰.尤其是随着人类文明的不断进步和经济的高速发展,无论是研究深部地质构造还是寻找深部隐伏盲矿体,都不可能完全避开上述复杂的电磁干扰环境 (肖晓等,2014).因此,运用现代信号处理技术剖析矿集区大地电磁信号与强干扰的本质特征,从各种强干扰中精细分离出微弱的大地电磁有用信号成为在矿集区开展大地电磁测深无法回避的实际问题,同时也是一项极具挑战性的任务.

大地电磁测深理论提出至今,噪声问题一直困扰着广大大地电磁研究者,如何消除大地电磁信号中的噪声干扰,提高大地电磁测深数据质量,是国内外长期瞩目并不断取得进展的研究课题.诸多现代信号处理方法,如自功率谱和互功率谱 (Vozoff,1972Kao and Rankin, 1977Goubau et al., 1978)、远参考道 (Gamble et al., 1979)、Robust (Egbert and Booker, 1986Sutamo and Vozoff, 1991)、高阶统计量 (王书明和王家映,2004)、Hilbert-Huang变换 (汤井田等,2008Cai et al., 2009Chen et al., 2012)、广义S变换 (景建恩等,2012)、方差比维纳滤波 (Kapple,2012)、数学形态滤波 (Li et al., 2011汤井田等, 2012a, 2012b李晋等,2014b)、子空间增强 (李晋等,2014a)、同步时间序列依赖 (王辉等,2014)、小波去噪及相关综合算法 (Trad and Travassos, 2000范翠松等,2008凌振宝等,2016) 等均被应用到该领域,并在一定程度上对大地电磁测深数据质量的改善起到了积极作用,推动了大地电磁测深法在与各种噪声竞争中不断成长.

分析国内外相关文献可知,已有的大地电磁强干扰压制方法均是对采集到的大地电磁数据进行整体处理,缺乏对大地电磁信号和强干扰进行信噪辨识的环节.若能从大地电磁信号和强干扰的本质特征入手剖析大地电磁信号和强干扰之间的定性及定量辨识关系,将对后续的大地电磁高精细分离起到关键作用,并能有效避免已有技术在大地电磁噪声压制中的过处理、最大限度地保留大地电磁信号的低频缓变化信息.

递归分析起源于20世纪30年代,由于工业革命的开始与普及,人们除了对地球上的物质、现象进行观察与研究外,也开始对宇宙射线的递归行为进行探索.20世纪80年代,衍生了递归图 (Recurrence Plot,RP) 的概念,即将时间序列的递归特性刻画为二维相空间 (Eckmann et al., 1987);此后不断有学者提出定量分析递归图的方法,从此递归分析走向成熟并逐渐被成功应用于盲源提取 (Soriano et al., 2011)、故障诊断 (杨栋和任新伟,2012)、脑电检测 (孟庆芳等,2014)、太阳活动 (周双等,2015)、电力系统 (Bhui and Senroy, 2016) 等领域.聚类分析是研究分类问题的一种多元统计方法,作为一种非督导分类方法,聚类分析主要用于把分类对象按照一定规则分成若干类,这些类不是事先给定的,而是根据数据的特征确定的.该方法已被广泛应用于图像分割 (Porter and Canagarajah, 1996)、模式识别 (Ferrari-Trecate et al., 2003)、数据挖掘 (孙庆先等,2008) 等领域.本文结合递归分析和聚类,尝试性地对典型的大地电磁强干扰和微弱的大地电磁信号进行定性和定量剖析,并仅对辨识为强干扰的大地电磁时间序列进行噪声压制.仿真结果表明,该方法能有效避免已有技术整体处理过程中的盲目性,其结果更为精细地反映了测点本身所固有的丰富、可靠的低频信息,卡尼亚电阻率-相位曲线的整体形态更为光滑、连续,整个低频段的数据质量得到了明显改善.

2 递归分析 2.1 递归图

递归图是一种非线性动力学分析方法,主要用于分析时间序列的周期性、非平稳性和混沌性.它以相空间重构为基础,通过构造多组向量将信号某一瞬时状态与相应时延后的瞬时状态联系起来,表达信号内在的递归特性 (Marwan et al., 2007).为此,递归图具有准确反映信号演变过程的特点,该方法通过揭示隐藏在非平稳系统中的递归特性来描述信号的本质特征,从而反映恢复后的混沌吸引子所具有的某种规律,并得到相似性、信息量和预测性等相关先验知识.

N点一维时间序列o(i) 为例,重构m维相空间得到动力学吸引子Oi如下所示:

(1)

式中,m表示嵌入维数,τ表示延迟时间.当mτ选择恰当时,Oi可用来刻画相空间轨迹.

在相空间重构的基础上,递归值Rij定义为:

(2)

式中,‖·‖表示Euclidean范数,用来计算重构相空间中吸引子OiOj的欧式距离;ε表示判别阈值,H(x) 表示Heaviside函数,在递归图上用来刻画吸引子OiOj的演变规律.

通过上述步骤即把一维时间序列o(i) 重构成了m维的相空间轨迹Oi,从动力学系统上实现了在高维空间恢复吸引子.当Rij=1时,在递归图坐标 (i, j) 位置用一个黑点来表示,代表这两点之间的距离是递归的——相互逼近,否则用一个白点 (Rij=0) 或空格来表示,代表两者之间不具有递归性——相互远离.这些沿着几何空间轨迹移动的每一个递归点都对应系统的一个状态,通过对这些黑、白点的分布情况进行研究,在拓扑等价意义下可获得原系统的动力学特性.因此,递归图是状态空间内轨道i时刻的状态关于j时刻的递归现象,其本质是一种时-时信号处理方法,它通过黑、白相间的点来描绘二维图形,从而观测多维动力学系统的内部机理.该方法可以定性反映一维时间序列的轨迹特征和递归特性,并能测定动力学系统中时间序列的非平稳性和内在相关性.

2.2 递归定量分析

递归图虽然能够准确反映非线性系统演化过程中的内在特征,并能从高维相空间重构吸引子分析信号所包含的动力学信息,但这种方法只能进行定性的分析.实际研究中,我们更希望能从定量的视角剖析递归图所反映的特征.为此,Zbilut和Webber (1992)对递归图进行研究,提出了递归定量分析 (Recurrence Quantification Analysis,RQA) 方法;该方法通过统计递归图中的点和线段来定量刻画递归图中小尺度的微观结构 (Trulla et al., 1996).常用的递归定量分析参数有:递归率 (Recurrence Ratio,RR)、确定率 (Determinism,DET)、熵 (Entropy,ENTR) 和平均对角线长度 (Lmean,Lmean),分别定义如下:

递归率 (RR):

(3)

递归率表示递归图中黑点所占总点数的百分比,用来确定递归图中递归点的密度.当信号突变时,递归率会明显增加.

确定率 (DET):

(4)

式中,lmin表示对角线结构中所取的长度初值,一般取lmin=2;P(l) 表示平行于主对角线线段中长度为l的分布概率.确定率表示平行于主对角线45°递归点的分布情况,用来描述动力学系统的确定性.

熵 (ENTR):

(5)

熵表示平行于主对角线线段长度分布的Shannon熵,用来描述递归图中对角线的复杂度.

平均对角线长度 (Lmean):

(6)

平均对角线长度表示平行于主对角线线段长度的平均值,用来反映信号变化的随机性.平均对角线长度越短,信号变化则越随机.

因此,从系统演化特征出发的RQA参数能定量描述不同系统的动力学行为,不同RQA参数表达的意义也不尽相同,将这些参数有效结合可用来分析待处理信号的本质特征.

3 聚类分析 3.1 K均值聚类

K均值 (K-means) 聚类是一种非监督实时聚类算法,由于算法简洁高效,在聚类分析中使用最为广泛 (Likas et al., 2003).算法采用距离作为相似性的评价指标,其核心思想是通过迭代将数据对象划分为不同的类别,使得评价聚类性能的标准测度函数达到最优.

对于包含n个数据对象的数据集X,采用K-means聚类过程如下:

首先,任意选取K个对象作为初始聚类中心C.然后,根据最小距离原则将数据对象分配至最接近的聚类中心形成一类.最后,更新聚类中心、重新分配数据对象,并不断重复上述步骤直至标准测度函数开始收敛为止,即每个聚类不再发生变化.

标准测度函数J(C, X) 定义如下:

(7)

式中,d(Ci, Xj)=‖CiXj‖表示聚类中心Ci与数据对象Xj之间的欧式距离.

由于K-means聚类中K值需要事先给定,若选到噪声或孤立点,则往往导致算法迭代次数增加、聚类性能下降.

3.2 模糊C均值聚类

传统聚类分析方法把每个待处理对象严格划分到某一类别,具有非此即彼的特性.然而,在许多实际分类问题中,分类对象之间并没有明确的界限,反映出亦此亦彼的特征.因此,随着模糊集理论的提出,人们开始用模糊的方法来处理聚类问题;由于模糊聚类能够更为客观地反映现实世界,为此该方法逐渐引领了聚类分析的发展趋势.

模糊C均值聚类 (Fuzzy C-Means Clustering,FCM) 是一种最为典型的无监督模糊聚类方法 (Bezdek,1973).与K-means聚类不同,FCM算法在实现过程中不需要人为干预,是一种柔性的模糊划分.该算法把N个数据对象X划分为c个模糊类,通过计算每个类的聚类中心vi,使得目标函数值最小.在模糊分类中,每一个数据不是严格地划分为某一类,而是以一定的隶属度归于某一类 (Rezaee,2010).

目标函数定义如下:

(8)

式中,uik表示第k个待分类对象属于第i个聚类中心的隶属度,dik=‖xkvi‖表示数据对象xk与聚类中心vi之间的欧式距离,V={v1, v2, …, vc}表示c个类的聚类中心,J(U, V) 表示各类对象与聚类中心加权平均距离之和,m∈[1, ∞) 是一个加权指数,U=(uik)c×n表示由每个对象与相应聚类中心隶属度构成的隶属矩阵,约束条件如下:

(9)

(10)

模糊C均值聚类即在式 (9) 和式 (10) 的聚类准则下求UV,使得J(U, V) 最小.这是关于 (U, V) 的约束优化问题,通过拉格朗日乘法对目标函数求导,可得到viuik分别如下所示:

(11)

(12)

最后,根据求解的隶属度矩阵U中元素的取值,可以确定待分类对象的归属.

4 仿真实验分析 4.1 递归图参数选取

为了分析递归图的递归性能,构造一组测试信号进行仿真实验.测试信号包含高斯随机信号、正弦信号以及两者叠加的正弦带噪信号,其时域波形图和递归图如图 1所示.

图 1 测试信号时域波形图和递归图 Fig. 1 Time domain waveform figures and recurrence plots of signal test

分析图 1可知,高斯随机信号递归图中的递归点均匀分布在主对角线附近,空白区域呈随机、离散状态,较好地体现了该信号所具有的随机分布的递归特征;正弦信号递归图中的递归点平行于主对角线,且呈现出周期性的棋盘状结构,充分表明了该信号本身所具有的周期性的递归特性;正弦带噪信号的递归图大体上仍呈棋盘结构,但递归点有离散趋势,较好地反映了高斯随机信号和正弦信号两者之间的递归特征.

表 1所示为测试信号的递归定量分析参数.

表 1 测试信号递归定量分析参数 Table 1 RQA parameters of signal test

分析表 1可知,高斯随机信号的各项参数指标均为最低,从RQA参数的定义可以解释为该信号包含的随机性应最强;正弦信号的各项参数指标均为最高,表明此时信号具有非常明显的确定率和周期性;正弦带噪信号的各项递归参数介于两者之间,说明该信号包含的周期性不是特别明显,但比随机状态强烈.

接下来,以测试信号中的高斯随机噪声和正弦信号为例,分别讨论嵌入维数、延迟时间和判别阈值对递归图中递归特性的影响情况.

图 2所示为延迟时间 (τ=1) 和判别阈值 (ε=1) 一定时,不同嵌入维数 (m=3, 4, 5, 6) 递归图.

图 2 不同嵌入维数递归图 Fig. 2 Recurrence plots of different embedding dimensions

图 2可知,选择合适的嵌入维数对反映待分析信号的递归特性起着至关重要的作用.当嵌入维数m减小时,黑点逐渐增加,直至覆盖整个相空间轨迹;当嵌入维数m增加时,黑点逐渐减小,反映相空间轨迹特性的递归特征逐渐减弱、直至消失.

图 3所示为嵌入维数 (m=3) 和判别阈值 (ε=1) 一定时,不同延迟时间 (τ=1, 2, 3, 4) 递归图.

图 3 不同延迟时间递归图 Fig. 3 Recurrence plots of different delay times

图 3可知,不同信号的延迟时间所起的作用不尽相同.对于周期性的正弦信号,延迟时间影响比较明显;当延迟时间增加时,黑点减少,白色区域增大,体现相空间轨迹特性的棋盘结构随之减弱.对于高斯随机信号,延迟时间并不敏感,递归点仍然均匀分布在主对角线两侧.

图 4所示为嵌入维数 (m=3) 和延迟时间 (τ=1) 一定时,不同判别阈值 (ε=1, 2, 3, 4) 递归图.

图 4 不同判别阈值递归图 Fig. 4 Recurrence plots of different determination thresholds

图 4可知,判别阈值的选取对待分析信号的递归特性起到关键作用.当判别阈值过大时,黑色区域变大,递归点将湮没相空间轨迹,此时递归图的递归特性难以分析与观测.

通过分析大量仿真实验发现,当m=3、τ=1、ε=1时,反映递归图相空间的吸引子能够完全展开,相空间轨迹的递归特性可以得到较好的刻画.

图 5所示为一段实测MT数据在不同的嵌入维数 (m)、延迟时间 (τ) 和判别阈值 (ε) 时的递归图.

图 5 实测MT数据递归图 Fig. 5 Recurrence plots of measured MT data

分析图 5可知,针对测试信号讨论的不同嵌入维数、延迟时间和判别阈值下递归特性的结论适合于分析一般的实测MT数据.当m=3、τ=1、ε=1时,实测MT数据的递归特性能得到较好的展示,这为后续选取实测MT数据中的嵌入维数、延迟时间和判别阈值提供了依据.

4.2 不同长度序列对递归定量分析参数的影响

由于实测大地电磁数据量庞大,为了分段处理采集到的大地电磁信号,有必要讨论不同长度序列对递归定量分析参数的性能影响情况.

图 6所示为具有相同频率、不同长度 (L=300、400、500、600、700) 的正弦信号时域波形图和递归图.

图 6 不同长度序列的时域波形图和递归图 Fig. 6 Time domain waveform figures and recurrence plots for the sequence of different lengths

表 2所示为上述不同长度序列的递归定量分析参数.

表 2 不同长度序列的递归定量分析参数 Table 2 RQA parameters for the sequence of different lengths

分析对比图 6表 2可知,对于具有类似特征的待分析信号,即使信号序列的长度发生变化,递归图的形态特征以及递归定量分析参数值均相对稳定,其递归特性并没有发生突变,这为后续处理实测大地电磁数据的分段奠定了基础.

4.3 样本库构建

为了剖析矿集区典型强干扰与天然微弱大地电磁信号之间的定量辨识关系,首先必须获取矿集区含典型强干扰类型的MT信号和几乎未受到电磁干扰的MT信号.

根据前期的研究结果表明,在庐枞和铜陵矿集区采集到的海量大地电磁测深数据中,通常充斥大量类方波、类充放电三角波和类脉冲干扰.这些大尺度强干扰的能量和幅值远超过正常大地电磁信号几个数量级,微弱的大地电磁有用信号几乎被完全湮没 (汤井田等,2012c).庐枞矿集区某测点电道和磁道的大地电磁测深数据时间片段如图 7所示.

图 7 庐枞矿集区某测点大地电磁测深数据时间片段 Fig. 7 Time segment of magnetotelluric sounding data in the Luzong ore concentration area

图 7可知,由于测点所处位置面临复杂的电磁干扰环境,导致电道 (Ex/Ey) 和磁道 (Hy/Hx) 中伴随出现大量上述强干扰类型.这些强干扰的存在往往引起卡尼亚电阻率曲线呈45°渐近线快速上升、相位趋近于0°或180°左右的近源效应,测点本身所固有的电性结构信息则完全无法解释 (汤井田等,2015).为此,文中将这3类最常见的强干扰类型作为矿集区典型强干扰进行分析和讨论.

图 8所示为青海省海西蒙古族藏族自治州某测点电道和磁道大地电磁测深数据时间片段.

图 8 青海地区某测点大地电磁测深数据时间片段 Fig. 8 Time segment of magnetotelluric sounding data in the Qinghai area

图 8可知,该测点电道和磁道的时间域波形中均未观测到如图 7所示矿集区中常见的大尺度强干扰异常波形;测点的卡尼亚电阻率曲线光滑、连续,没有出现45°上升的近源效应,且视电阻率值相对平稳,相位在50°左右,这些均符合天然大地电磁场的正常逻辑.

由于青海地区人烟稀少、基本无明显电磁干扰,且采集到的大地电磁数据具备天然大地电磁场信号的本质特征,为此文中将在该地区采集到的大地电磁数据作为几乎未受到电磁干扰的大地电磁信号进行分析和讨论.

接下来,从上述矿集区任选含类方波干扰、类充放电三角波干扰和类脉冲干扰的MT时间序列各50个样本,以及从电磁干扰相对不明显的青海地区任选50个几乎未受到电磁干扰的MT信号组成实验样本库,进一步探讨矿集区典型强干扰类型和大地电磁微弱信号之间的信噪辨识关系.样本库中每个样本的长度均为240个采样点,共计200个样本.

4.4 信噪辨识

为了甄辨大地电磁信号和典型强干扰,从样本库中随机抽取一组几乎未受到电磁干扰的MT信号,以及含类方波干扰、类充放电三角波干扰和类脉冲干扰的MT信号进行仿真实验,时域波形图和递归图如图 9所示,相应的递归定量分析参数如表 3所示.

图 9 各类实测大地电磁信号时域波形图和递归图 Fig. 9 Time domain waveform figures and recurrence plots of all kinds of measured MT signal
表 3 各类实测大地电磁信号的递归定量分析参数 Table 3 RQA parameters of all kinds of measured MT signal

分析图 9表 3可知,几乎未受到电磁干扰的MT信号其递归图中的递归点均匀分布,时间序列变得无法预测;反映在相空间轨迹的聚集度明显分散,递归图几乎毫无规律,各项RQA参数值均低于其他三类含强干扰的MT信号;根据RQA参数的物理意义可知,此时信号应最无规律可循,并呈随机、不确定状态,符合天然MT信号是随机分布的特征.含类方波干扰和类充放电三角波干扰的MT信号其递归点密集分布在与主对角线平行的直线两侧,递归图中出现大片聚集的黑色或白色区域,且呈现出很有规则的图形、层次感分明;其各项RQA参数指标明显增大,从动力学机理可以解释为,此时信号所包含的确定性和可预测性最为强烈.含类脉冲干扰的MT信号由于只在很小的时间段内出现突变,与几乎未受到电磁干扰的MT信号在波形形态上非常相似,为此递归图中并未出现特别规则的图案,但相空间轨迹的聚集度有缓慢集中的趋势;各项RQA参数指标仅比几乎未受到电磁干扰的MT信号略高,说明该段信号不具备含类方波干扰和类充放电三角波干扰时明显强烈的类周期性.

图 10所示为样本库中四类大地电磁信号递归定量分析参数分布图,其中m=3、τ=1、ε=1.

图 10 递归定量分析参数分布 Fig. 10 Distribution of RQA parameters

图 10可知,几乎未受到电磁干扰的MT信号其4组RQA参数值均具有很小的波动性,含强干扰的3类MT信号特征曲线虽然出现跳变和混叠现象,但与几乎未受到电磁干扰的MT信号RQA参数相比,仍具有良好的区分度.尤其是确定率和平均对角线长度,含典型强干扰的MT信号其确定率几乎保持在同一水平,说明这3类强干扰中包含的规则度明显增加;几乎未受到电磁干扰的MT信号的平均对角线长度则远小于其他3类含强干扰的MT信号,从RQA参数的物理意义可解释为,此时信号包含的随机状态信息应最多、复杂度最高.

上述仿真结果表明,递归图具有从高维相空间重构吸引子的特征,其相空间轨迹能对大地电磁时间序列确定性成分的存在和周期性成分的嵌入进行刻画;该方法较好地描述了非线性系统的内在特征,且能从某种角度辨识信号所包含的动力学行为,适合定性判断大地电磁时间序列的非稳态动态变化.反映递归图中微观结构的RQA参数在区分微弱大地电磁信号和典型强干扰时具有较好的噪声鲁棒性,可以更为直观地体现大地电磁时间序列的动力学信息及递归特性,适合定量区分和甄辨大地电磁信号和典型强干扰.

图 11所示为针对样本库中200个样本提取RQA参数进行K-means聚类的轮廓图.其中,聚类中心选取为[15,100].

图 11 样本库K-means聚类轮廓图 Fig. 11 Silhouette figure of K-means clustering for the sample dataset

分析图 11可知,K-means聚类将样本库中200个样本分为两类,50个几乎未受到电磁干扰的MT信号归为第1类,其余150个含强干扰的MT信号归为第2类,聚类效果良好.

图 12所示为针对样本库中200个样本提取RQA参数进行模糊C均值聚类的效果图.

图 12 样本库模糊C均值聚类效果图 Fig. 12 Fuzzy C-means clustering of the sample dataset

图 12可知,模糊C均值聚类自动将样本库中几乎未受到电磁干扰MT信号和含强干扰MT信号辨识为两类,聚类效果良好.由于K-means聚类的关键在于选取合适的聚类中心,当样本数据凌乱时,只能随机选取聚类中心的初值,导致K-means聚类的效果不够理想;模糊C均值聚类通过优化目标函数得到每个样本点对所有类中心的隶属度,从而达到对样本数据自动分类的目的,因此该方法更适合于海量实测大地电磁数据的甄辨与分类,后续实测大地电磁数据均采用模糊C均值聚类进行处理.

5 实测资料处理 5.1 时间域处理效果

文中实测MT时间序列来源于V5-2000大地电磁测深仪采集到的MT数据,以TSL、TSH文件为例,由于低频段的MT数据为全时段采集 (采样率为24 Hz),中、高频交替采集 (采样率分别为320 Hz和2560 Hz),因此反映测点电性结构信息的数据主要集中在低频段,一般而言只需要对TSL文件中的MT数据进行处理.为此,基于递归分析和聚类的大地电磁信噪辨识及分离方法的具体步骤如下:首先,将低频TSL文件中的MT时间序列 (ExEyHxHy) 进行分帧,每240个采样点 (10 s) 为1帧,帧与帧之间不重叠.然后,求解电道 (ExEy) 中每1帧信号的RQA参数,针对RQA参数采用模糊C均值聚类将MT时间序列分为两类.接着,通过分析对比两类RQA参数值的大小,判断几乎未受到电磁干扰MT信号与含强干扰MT信号所在的类,并仅对辨识为强干扰的MT数据运用数学形态滤波进行信噪分离处理.最后,将形态滤波处理后的MT信号与辨识为低频缓变化的MT信号相结合,重构MT有用信号.处理过程中,磁道 (HxHy) MT信号则遵循相关性的原则,即当电道Ex/Ey中某1帧辨识为典型强干扰时,磁道Hy/Hx相同时刻对应的帧同样做相应的形态滤波处理.

图 13所示为一段含大尺度类方波干扰和类充放电三角波干扰的MT实测数据信噪辨识效果图.

图 13 MT实测数据信噪辨识效果 Fig. 13 Signal and noise identification of measured MT data

分析图 13可知,含大尺度类方波干扰和类充放电三角波干扰的MT信号,经本文方法首先通过对MT序列提取鲁棒性的RQA参数,然后根据RQA参数进行模糊C均值聚类后,可以准确地辨识出MT序列中典型强干扰和低频缓变化信息出现的时间段,为后续有针对性地压制MT强干扰指明了方向.

表 4所示为图 13中所标识位置帧对应的递归定量分析参数.分析表 4可知,第1帧和第3帧所对应的RQA参数明显低于第2帧和第4帧对应的RQA参数,尤其是Lmean参数差异显著,从RQA参数的物理意义可以判断第2帧和第4帧所包含的确定性和规律性应明显强于第1帧和第3帧.对比观测图 13中MT实测序列发现,第1帧和第3帧的MT序列中并未出现类似于典型强干扰的异常波形,而第2帧和第4帧的MT序列中包含典型的大尺度类方波干扰和类充放电三角波干扰.由此可知,通过提取MT序列中鲁棒性的RQA参数,可以从非线性动力学上定量辨识大尺度强干扰和微弱的MT有用信号.

表 4 MT实测数据所标识位置对应的递归定量分析参数 Table 4 RQA parameters of measured MT data for the identified location

图 14所示为图 13中对应的MT实测数据信噪分离效果图.

图 14 MT实测数据信噪分离效果 Fig. 14 Signal and noise separation of measured MT data

分析图 14可知,形态滤波整体处理虽然能够压制大尺度类方波和类充放电三角波干扰,但由于方法缺少信噪辨识环节,原本微弱的MT信号的缓变化信息也被随之过处理;处理后的MT信号均在基线附近,导致卡尼亚电阻率曲线低频段的MT数据质量堪忧.本文方法在信噪辨识的基础上,可以有针对性地对辨识为大尺度强干扰的MT信号进行信噪分离,较好地避免了已有技术整体处理时的盲目性;处理后的信号在压制大尺度强干扰的同时,有效地保留了几乎未受到电磁干扰或受干扰程度不明显时间段的缓变化信息,重构的MT信号中包含了更多丰富的反映MT低频信息的细节成分.

5.2 实测点处理效果

图 15图 16所示分别为庐枞矿集区测点B30242和B30743的卡尼亚电阻率-相位曲线.从图 15图 16可知,两测点原始数据的视电阻率曲线大约在50~0.05 Hz频段以近45°渐进线快速上升;在5~0.05 Hz频段,对应的相位呈0°或±180°,近14个频点的数据完全失真;在0.05 Hz左右,xy方向和yx方向的视电阻率值竟然超过100000 Ωm左右;在0.05~0.01 Hz频段,视电阻率值突然下降至1000 Ωm左右.该现象属于典型的近源效应,其结果已不能真实反映地下电性结构.

图 15 测点B30242原始数据的卡尼亚电阻率-相位曲线 Fig. 15 Cagniard resistivity-phase curve of original data for site B30242
图 16 测点B30743原始数据的卡尼亚电阻率-相位曲线 Fig. 16 Cagniard resistivity-phase curve of original data for site B30743

图 17图 20所示分别为上述两测点经形态滤波整体处理和本文方法处理后的卡尼亚电阻率-相位曲线.

图 17 测点B30242经形态滤波整体处理的卡尼亚电阻率-相位曲线 Fig. 17 Cagniard resistivity-phase curve of overall processing by morphological filtering for site B30242
图 18 测点B30743经形态滤波整体处理的卡尼亚电阻率-相位曲线 Fig. 18 Cagniard resistivity-phase curve of overall processing by morphological filtering for site B30743
图 19 测点B30242经本文方法处理的卡尼亚电阻率-相位曲线 Fig. 19 Cagniard resistivity-phase curve of the proposed method for site B30242
图 20 测点B30743经本文方法处理的卡尼亚电阻率-相位曲线 Fig. 20 Cagniard resistivity-phase curve of the proposed method for site B30743

图 17图 18可知,两测点原始数据经形态滤波整体处理后,视电阻率曲线近45°上升的现象有所缓解,相位曲线有所抬升.然而,在0.1~0.0005 Hz超低频段,视电阻率曲线均呈下掉趋势,尤其是xy方向的视电阻率曲线不光滑、视电阻率值突跳明显,对应的相位曲线变化紊乱、误差棒增大.究其原因,形态滤波整体处理虽能有效压制大尺度强干扰,但一些反映测点深部构造信息的低频缓变化成分也被随之剔除;处理后的数据均在基线附近,其结果俨然已不能精确反映地下电性结构.

分析图 19图 20可知,经本文方法先对大地电磁信号进行信噪辨识,再有针对性地去除强干扰后,两测点的视电阻率曲线整体形态更加平稳,视电阻率值相对稳定、最大值下降了1个数量级;相位曲线显著抬升,相位达到50°左右的正常水平,近源干扰得到了明显压制.在5~0.05 Hz频段,近14个频点的数据质量均得到了有效改善;相比形态滤波整体处理,在0.1~0.0005 Hz超低频段,视电阻率曲线下掉的现象有明显缓解,且卡尼亚电阻率-相位曲线光滑、连续,误差棒减小.

为进一步说明本文方法的处理效果,图 21给出了测点B30242经本文方法处理前后0.3 Hz电磁场极化方向对比图.

图 21 测点B30242处理前后0.3 Hz电磁场极化方向对比 Fig. 21 Comparison of polarization directions at 0.3 Hz for site B30242 before and after processing

分析图 21可知,测点原始数据的电磁场极化方向在某些极化角度上高度集中,表现出一定的规则性,从电磁场极化方向的基本原理可以说明,该测点采集过程中受到了强烈的电磁噪声干扰 (Weckmann et al., 2005).经本文方法处理后,电磁场极化方向较原始数据有明显改善,极化角度呈离散趋势,其随机无序性符合天然大地电磁场极化方向的基本特征,这在一定程度上可以说明该测点所包含的强噪声干扰得到了有效压制.

综上可知,本文方法在压制典型强干扰的同时,尽可能地保留了测点本身所固有的丰富、真实的低频信息,整个低频段的大地电磁数据质量得到了明显改善,其结果为后续提升电磁法反演解释水平提供了更为可靠、置信的大地电磁数据.

6 结论

针对现有技术在大地电磁信噪分离过程中缺乏信噪辨识环节、导致处理后的数据往往丢失低频缓变化这一问题,本文从非线性行为学的角度出发,剖析了微弱大地电磁信号和典型强干扰时间序列的递归特性,讨论了两者之间的定性及定量辨识关系,并结合聚类分析对实测大地电磁数据的递归定量分析参数进行信噪辨识,研究了基于递归分析和聚类的大地电磁信噪辨识及分离方法.通过实验仿真和实测数据处理,结果表明:递归图的相空间轨迹能对大地电磁时间序列确定性成分的存在和周期性成分的嵌入进行刻画,适合定性判断大地电磁时间序列的非稳态动态变化;递归定量分析能直观地表征大地电磁时间序列的动力学信息与递归特性,适合定量区分和甄辨微弱大地电磁信号和典型强干扰类型,且具有较好的噪声鲁棒性.结合聚类分析对大地电磁时间序列做信噪辨识处理,低频段的大地电磁缓变化信息得到了更为精细的保留,有效地避免了噪声压制过程中的盲目性;仅对辨识为强干扰的大地电磁数据运用数学形态滤波处理后的大地电磁时间序列中尽可能多地反映了天然大地电磁场所固有的本质特征,卡尼亚电阻率-相位曲线更为光滑、连续,整个低频段的大地电磁数据质量得到了明显改善.本文方法为今后在矿集区开展高精细的大地电磁噪声压制提供了一条新的研究思路,具有重要的理论意义与实际应用价值.

由于方法通过获取大地电磁时间序列的递归定量分析 (RQA) 参数来定量描述大地电磁时间序列的非线性特征,并结合模糊C均值聚类和数学形态滤波辨识、分离含强干扰的大地电磁数据,当测点包含的噪声干扰其时域波形具有大尺度形态特征或形态特征较为规则时,方法能取得较为满意的效果.然而,矿集区面临复杂的噪声干扰环境,当采集的大地电磁原始数据中的噪声干扰其形态特征不规则时,信号与噪声的分界逐渐模糊,信噪辨识的精度将有所下降.为此,在递归分析的基础上,如何寻求更为鲁棒、稳健的大地电磁信噪辨识特征参数,以及结合智能算法进一步提升聚类精度,将是今后的研究重点.

参考文献
Bezdek J C. 1973. Cluster validity with fuzzy sets. Journal of Cybernetics, 3(3): 58-73. DOI:10.1080/01969727308546047
Bhui P, Senroy N. 2016. Application of recurrence quantification analysis to power system dynamic studies. IEEE Transactions on Power Systems, 31(1): 581-591. DOI:10.1109/TPWRS.2015.2407894
Cai J H, Tang J T, Hua X R, et al. 2009. An analysis method for magnetotelluric data based on the Hilbert-Huang transform. Exploration Geophysics, 40(2): 197-205. DOI:10.1071/EG08124
Chen J, Heincke B, Jegen M, et al. 2012. Using empirical mode decomposition to process marine magnetotelluric data. Geophysical Journal International, 190(1): 293-309. DOI:10.1111/gji.2012.190.issue-1
Dong S W, Li T D, Cheng X H, et al. 2012. Progress of deep exploration in mainland China: A review. Chinese J. Geophys., 55(12): 3884-3901. DOI:10.6038/j.issn.0001-5733.2012.12.002
Eckmann J P, Kamphorst S O, Ruelle D. 1987. Recurrence plots of dynamical systems. Europhysics Letters, 4(9): 973-977. DOI:10.1209/0295-5075/4/9/004
Egbert G D, Booker J R. 1986. Robust estimation of geomagnetic transfer functions. Geophysical Journal International, 87(1): 173-194. DOI:10.1111/gji.1986.87.issue-1
Fan C S, Li T L, Wang D Y. 2008. Treatment of wavelet transform for square wave noise in MT data. Journal of Jilin University (Earth Science Edition), 38(S1): 61-63.
Ferrari-Trecate G, Muselli M, Liberati D, et al. 2003. A clustering technique for the identification of piecewise affine systems. Automatica, 39(2): 205-217. DOI:10.1016/S0005-1098(02)00224-8
Gamble T D, Goubau W M, Clarke J. 1979. Magnetotellurics with a remote magnetic reference. Geophysics, 44(1): 53-68. DOI:10.1190/1.1440923
Goubau W M, Gamble T D, Clarke J. 1978. Magnetotelluric data analysis: removal of Bias. Geophysics, 43(6): 1157-1166. DOI:10.1190/1.1440885
Jing J E, Wei W B, Chen H Y, et al. 2012. Magnetotelluric sounding data processing based on generalized S transformation. Chinese J. Geophys., 55(12): 4015-4022. DOI:10.6038/j.issn.0001-5733.2012.12.013
Kao D W, Rankin D. 1977. Enhancement of signal-to-noise ratio in magnetotelluric data. Geophysics, 42(1): 103-110. DOI:10.1190/1.1440703
Kapple K N. 2012. A data variance technique for automated despiking of magnetotelluric data with a remote reference. Geophysical Prospecting, 60(1): 179-191. DOI:10.1111/gpr.2012.60.issue-1
Li J, Tang J T, Xiao X. 2011. De-noising algorithm for magnetotelluric signal based on mathematical morphology filtering. Noise & Vibration Worldwide, 42(11): 65-72.
Li J, Tang J T, Wang L, et al. 2014a. Noise suppression for magnetotelluric sounding data based on signal subspace enhancement and endpoint detection. Acta Phys. Sin., 63(1): 019101.
Li J, Tang J T, Xiao X, et al. 2014b. Magnetotelluric data processing based on combined generalized morphological filter. Journal of Central South University (Science and Technology), 45(1): 173-185.
Likas A, Vlassis N, Verbeek J J. 2003. The global K-means clustering algorithm. Pattern Recognition, 36(2): 451-461. DOI:10.1016/S0031-3203(02)00060-2
Ling Z B, Wang P Y, Wan Y X, et al. 2016. A combined wavelet transform algorithm used for de-noising magnetotellurics data in the strong human noise. Chinese J. Geophys., 59(9): 3436-3447. DOI:10.6038/cjg20160926
Lü Q T, Dong S W, Tang J T, et al. 2015. Multi-scale and integrated geophysical data revealing mineral systems and exploring for mineral deposits at depth: A synthesis from SinoProbe-03. Chinese J. Geophys., 58(12): 4319-4343. DOI:10.6038/cjg20151201
Lü Q T, Shi D N, Tang J T, et al. 2011. Probing on deep structure of middle and lower reaches of the Yangtze Metallogenic Belt and typical ore concentration area: a review of annual progress of SinoProbe-03. Acta Geoscientica Sinica, 32(3): 257-268.
Marwan N, Romano M C, Thiel M, et al. 2007. Recurrence plots for the analysis of complex systems. Physics Reports, 438(5-6): 237-329. DOI:10.1016/j.physrep.2006.11.001
Meng Q F, Chen S S, Chen Y H, et al. 2014. Automatic detection of epileptic EEG based on recurrence quantification analysis and SVM. Acta Phys. Sin., 63(5): 050506.
Porter R, Canagarajah N. 1996. A robust automatic clustering scheme for image segmentation using wavelets. IEEE Transactions on Image Processing, 5(4): 662-665. DOI:10.1109/83.491343
Rezaee B. 2010. A cluster validity index for fuzzy clustering. Fuzzy Sets and Systems, 161(23): 3014-3025. DOI:10.1016/j.fss.2010.07.005
Soriano D C, Suyama R, Attux R. 2011. Blind extraction of chaotic sources from mixtures with stochastic signals based on recurrence quantification analysis. Digital Signal Processing, 21(3): 417-426. DOI:10.1016/j.dsp.2010.12.003
Sun Q X, Chen Q P, Fang T, et al. 2008. Multi-scale spatial data mining model based on fuzzy clustering and its application in mine. Journal of Shanghai Jiaotong University, 42(2): 194-197.
Sutamo D, Vozoff K. 1991. Phase-smoothed robust M-estimation of magnetotelluric impedance functions. Geophysics, 56(12): 1999-2007. DOI:10.1190/1.1443012
Tang J T, Hua X R, Cao Z M, et al. 2008. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data. Chinese J. Geophys., 51(2): 603-610. DOI:10.3321/j.issn:0001-5733.2008.02.034
Tang J T, Li J, Xiao X, et al. 2012a. Mathematical morphology filtering and noise suppression of magnetotelluric sounding data. Chinese J. Geophys., 55(5): 1784-1793. DOI:10.6038/j.issn.0001-5733.2012.05.036
Tang J T, Li J, Xiao X, et al. 2012b. Magnetotelluric sounding data strong interference separation method based on mathematical morphology filtering. Journal of Central South University (Science and Technology), 43(6): 2215-2221.
Tang J T, Xu Z M, Xiao X, et al. 2012c. Effect rules of strong noise on Magnetotelluric (MT) sounding in the Luzong ore cluster area. Chinese J. Geophys., 55(12): 4147-4159. DOI:10.6038/j.issn.0001-5733.2012.12.027
Tang J T, Liu Z J, Liu F Y, et al. 2015. The denoising of the audio magnetotelluric data set with strong interferences. Chinese J. Geophys., 58(12): 4636-4647. DOI:10.6038/cjg20151225
Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of the Earth's crust. Dokl. Acad. Nauk. SSSR, 73(2): 295-297.
Trad D O, Travassos J M. 2000. Wavelet filtering of magnetotelluric data. Geophysics, 65(2): 482-491. DOI:10.1190/1.1444742
Trulla L L, Giuliani A, Zbilut J P, et al. 1996. Recurrence quantification analysis of the logistic equation with transients. Physics Letter A, 223(4): 255-260. DOI:10.1016/S0375-9601(96)00741-4
Vozoff K. 1972. The magnetotelluric method in the exploration of sedimentary basins. Geophysics, 37(1): 98-141. DOI:10.1190/1.1440255
Wang H, Wei W B, Jin S, et al. 2014. Removal of magnetotelluric noise based on synchronous time series relationship. Chinese J. Geophys., 57(2): 531-545. DOI:10.6038/cjg20140218
Wang S M, Wang J Y. 2004. Application of higher-order statistics in magnetotelluric data processing. Chinese J. Geophys., 47(5): 928-934. DOI:10.3321/j.issn:0001-5733.2004.05.027
Weckmann U, Magunia A, Ritter O. 2005. Effective noise separation for magnetotelluric single site data processing using a frequency domain selection scheme. Geophysical Journal International, 161(3): 635-652. DOI:10.1111/gji.2005.161.issue-3
Wei W B. 2002. New advance and prospect of magnetotelluric sounding (MT) in China. Progress in Geophysics, 17(2): 245-254.
Xiao X, Wang X Y, Tang J T, et al. 2014. Conductivity structure of the Lujiang-Zongyang ore concentrated area, Anhui Province: constraints from magnetotelluric data. Acta Geologica Sinica, 88(4): 478-495.
Yang D, Ren X W. 2012. Structure damage detecting using singular entropy of recurrence matrix. Journal of Vibration and Shock, 31(3): 60-63.
Zbilut J P, Webber C L. 1992. Embeddings and delays as derived from quantification of recurrence plots. Physics Letters A, 171(3-4): 199-203. DOI:10.1016/0375-9601(92)90426-M
Zhou S, Feng Y, Wu W Y. 2015. Chaos and fractal properties of solar activity phenomena at the high and low latitudes. Acta Phys. Sin., 64(24): 249601.
董树文, 李廷栋, 陈宣华, 等. 2012. 我国深部探测技术与实验研究进展综述. 地球物理学报, 55(12): 3884–3901. DOI:10.6038/j.issn.0001-5733.2012.12.002
范翠松, 李桐林, 王大勇. 2008. 小波变换对MT数据中方波噪声的处理. 吉林大学学报 (地球科学版), 38(S1): 61–63.
景建恩, 魏文博, 陈海燕, 等. 2012. 基于广义S变换的大地电磁测深数据处理. 地球物理学报, 55(12): 4015–4022. DOI:10.6038/j.issn.0001-5733.2012.12.013
李晋, 汤井田, 王玲, 等. 2014a. 基于信号子空间增强和端点检测的大地电磁噪声压制. 物理学报, 63(1): 019101.
李晋, 汤井田, 肖晓, 等. 2014b. 基于组合广义形态滤波的大地电磁资料处理. 中南大学学报 (自然科学版), 45(1): 173–185.
凌振宝, 王沛元, 万云霞, 等. 2016. 强人文干扰环境的电磁数据小波去噪方法研究. 地球物理学报, 59(9): 3436–3447. DOI:10.6038/cjg20160926
吕庆田, 董树文, 汤井田, 等. 2015. 多尺度综合地球物理探测:揭示成矿系统、助力深部找矿——长江中下游深部探测 (SinoProbe-03) 进展. 地球物理学报, 58(12): 4319–4343. DOI:10.6038/cjg20151201
吕庆田, 史大年, 汤井田, 等. 2011. 长江中下游成矿带及典型矿集区深部结构探测——SinoProbe-03年度进展综述. 地球学报, 32(3): 257–268.
孟庆芳, 陈珊珊, 陈月辉, 等. 2014. 基于递归量化分析与支持向量机的癫痫脑电自动检测方法. 物理学报, 63(5): 050506.
孙庆先, 陈秋平, 方涛, 等. 2008. 基于模糊聚类的多尺度空间数据挖掘模型及其矿山应用. 上海交通大学学报, 42(2): 194–197.
汤井田, 化希瑞, 曹哲民, 等. 2008. Hilbert-Huang变换与大地电磁噪声压制. 地球物理学报, 51(2): 603–610. DOI:10.3321/j.issn:0001-5733.2008.02.034
汤井田, 李晋, 肖晓, 等. 2012a. 数学形态滤波与大地电磁噪声压制. 地球物理学报, 55(5): 1784–1793. DOI:10.6038/j.issn.0001-5733.2012.05.036
汤井田, 李晋, 肖晓, 等. 2012b. 基于数学形态滤波的大地电磁强干扰分离方法. 中南大学学报 (自然科学版), 43(6): 2215–2221.
汤井田, 徐志敏, 肖晓, 等. 2012c. 庐枞矿集区大地电磁测深强噪声的影响规律. 地球物理学报, 55(12): 4147–4159. DOI:10.6038/j.issn.0001-5733.2012.12.027
汤井田, 刘子杰, 刘峰屹, 等. 2015. 音频大地电磁法强干扰压制试验研究. 地球物理学报, 58(12): 4636–4647. DOI:10.6038/cjg20151225
王辉, 魏文博, 金胜, 等. 2014. 基于同步大地电磁时间序列依赖关系的噪声处理. 地球物理学报, 57(2): 531–545. DOI:10.6038/cjg20140218
王书明, 王家映. 2004. 高阶统计量在大地电磁测深数据处理中的应用研究. 地球物理学报, 47(5): 928–934. DOI:10.3321/j.issn:0001-5733.2004.05.027
魏文博. 2002. 我国大地电磁测深新进展及瞻望. 地球物理学进展, 17(2): 245–254.
肖晓, 王显莹, 汤井田, 等. 2014. 安徽庐枞矿集区大地电磁探测与电性结构分析. 地质学报, 88(4): 478–495.
杨栋, 任新伟. 2012. 基于递归矩阵奇异熵的损伤识别方法. 振动与冲击, 31(3): 60–63.
周双, 冯勇, 吴文渊. 2015. 太阳高纬和低纬活动现象的混沌与分形特征. 物理学报, 64(24): 249601. DOI:10.7498/aps.64.249601