在地磁观测中有各种不同类型的干扰,其中,磁暴干扰是一种常见的干扰。如何快速准确识别磁暴干扰,对地磁数据分析工作具有重要意义。研究表明,大磁暴与7级以上地震活动存在相关性(张秀玲等,2017),采用图像处理方法可实现对地磁记录图的数字化识别(董星宏等,2009);支持向量机方法建模可对变化磁场进行预测(易世华等,2013);信号指纹技术结合短时傅里叶变换、小波变换等方法可实现对地磁干扰信号类型的识别(徐鹏深等,2018)。本文利用离散Gabor变换方法对地磁数据进行分析,使用支持向量机实现对正常数据与磁暴数据的自动分类识别,此外,针对地磁观测中其他类型的干扰,建立大数据样本库,也可使用该方法对未知类型数据进行自动分类识别,以减少人工分析判断工作。
1 离散Gabor变换为了对信号的局部特征进行分析,得到局部时间范围内频谱的变化,Dennis Gabor(1946)提出了一种时频变换方法——Gabor变换。Gabor变换是真正意义上的时频分析,它不仅能在整体上提供信号的所有频谱信息,还能提供局部时间内信号变化程度的信息(陈红等,2011)。
根据傅里叶变换可知,对时间变量的采样会导致频域的周期性,而对频率的采样又会导致时域的周期性。由于同时需要对时间和频率进行离散化,因此,离散Gabor展开只适用于离散时间的周期信号,设信号XL(k)周期为L,则离散Gabor展开定义为
$ {X_L}\left(k \right) = \sum\limits_{m = 0}^{M - 1} {} \sum\limits_{n = 0}^{N - 1} {{a_{mn}}} g\left({k - m\Delta M} \right){\exp ^{j2\pi nk\Delta N}} $ | (1) |
式中,ΔM与ΔN分别为时间和频率的采样间隔,而M、N分别为时间和频率采样的样本个数;核函数g(k)为高斯函数,amn为Gabor系数。
离散Gabor变换定义如下
$ {a_{mn}} = \sum\limits_{k = 0}^{L - 1} {{X_L}\left(k \right)} \gamma *\left({k - m\Delta M} \right){\exp ^{ - j2\pi nk\Delta N}} $ | (2) |
其中,分析窗函数γ(k)为核函数g(k)的对偶窗函数(上标*表示复共轭)。过采样率定义为
$ \alpha = \frac{L}{{\Delta M\Delta N}} $ | (3) |
若将MΔM = NΔN = L代入式(3),则
$ \alpha = \frac{{MN}}{L} $ | (4) |
当α = 1时,离散Gabor变换是临界采样的;当α>1时,离散Gabor变换是过采样的,即Gabor展开系数个数多于信号样本个数,此时,Gabor展开式有冗余。
2 支持向量机支持向量机是基于统计学理论的机器学习算法(蒋一然等,2019)。它以结构风险最小化为原则,能兼顾训练误差与测试误差的最小化。支持向量机的处理过程中,将所有样本分为训练样本和测试样本2类,首先,对训练样本进行训练,根据样本特征和样本类别训练形成分类器;然后,将未知类别的测试样本特征送入分类器测试,调用分类器对测试样本类别作出判断,实现对未知类别测试样本的分类识别。支持向量机算法的基本思想是,在n维空间中,找出可以将2类点划分开的最优分类平面, 使得该分类面不仅能对2类样本进行正确划分,而且分类间隔最大。支持向量机是一个二分器,它通过分类函数对测试样本集合进行预测,以实现对2类样本的分类。
3 计算结果与分析 3.1 秒数据Gabor变换整理收集陕西榆林台、河北红山台及昌黎台、陕西乾陵台和甘肃天水台共200组GM4型磁通门磁力仪的地磁秒数据,每组秒数据有D、H、Z三分量,故每组数据大小为86 400×3。数据分为2大类,一类是100组训练数据,另一类是100组测试数据,训练数据和测试数据各包括5个地磁台正常数据和磁暴干扰数据各10组。
2019年1月12日,榆林台GM4-1观测数据正常,无任何干扰。图 1为正常数据的H分量。由图 1可见,磁场变化平缓,曲线平滑,具有明显的日变形态。对H分量采用144个Gabor系数、过采样度为144、高斯窗长度为145的Gabor变换,得到图 2所示的Gabor变换谱图,该谱图显示了时频面内能量的分布,能量值就是Gabor系数模的平方,主要集中在低频处。
2018年4月20日,榆林台GM4-1型磁通门磁力仪的地磁秒数据受地磁暴干扰,磁暴发生于4月20日00:22至21:00时,急始变幅D:1.30′,H:25.8 nT,Z:5.7 nT,最大K指数为5。图 3为磁暴干扰日的H分量。由图 3可见,曲线波动剧烈,H分量呈现正脉冲变化,数据变化幅度较大。图 4为磁暴干扰日H分量的Gabor变换谱图。
对5个地磁台的200组地磁秒数据进行Gabor变换后,计算Gabor谱图的均值与方差,将其作为特征值得到数据特征值矩阵,将该矩阵做为样本数据。图 5为榆林台10组训练数据Gabor谱图的均值和方差。由图 5可见,正常数据及磁暴干扰数据经Gabor变换后,三分量Gabor谱图的均值与方差的数量级与分布均不同。
对100组训练样本的特征值进行训练,根据样本特征值和样本类别训练形成分类器。测试时,将100组测试样本的特征值送入分类器,调用分类器作出判决,该方法对100组测试样本的识别率为94%。图 6为5个地磁台测试样本识别率对比。由图 6可见,5个地磁台的正常数据测试样本均能100%识别;对于磁暴干扰数据测试样本,榆林台识别率90%,红山台90%,昌黎台90%,乾陵台70%,天水台100%。
5个地磁台站的磁暴错误识别数据日期如表 1所示。由表 1可见,2016年10月12日22时—14日23时发生磁暴,榆林台在2016年10月12日处于磁暴初相,磁暴刚开始时磁场值高于平静值,但扰动变化不大且时间较短,故识别错误。2016年7月19日23时—21日17时发生磁暴,红山台、昌黎台、乾陵台在2016年7月21日处于磁暴恢复相,磁场总扰动强度逐渐减弱至恢复平静状态,被错误识别。2016年3月11日5时—12日13时与2015年1月7日6时—8日18时发生磁暴,乾陵台在2016年3月12日与2015年1月8日均处于磁暴恢复相,也被错误识别。
本文对榆林台、红山台、昌黎台、乾陵台和天水台的正常数据与磁暴干扰数据进行分类识别分析,实验结果表明:①对无干扰的地磁正常数据均可正确识别,说明正常数据变化平稳,特征相对明显;②对磁暴干扰数据基本可以识别,50组测试数据中仅有6组错误识别,主要集中于磁暴的初相与恢复相,此时磁暴数据变化较小且持续时间较短,特征不明显,磁暴主相数据变化剧烈,均可正确识别;③磁暴在全球同时开始,同步变化,但是具有明显的经纬度差异。其中,红山台和昌黎台相距较近,所以磁暴变化基本相同,识别结果一致。其他台站由于经纬度不同,受磁暴干扰数据变化程度也不同;④对于本文的识别方法,增加样本数量能提高识别率,实现地磁数据磁暴干扰的自动判定,减少人工分析过程,避免人为识别错误。
陈红, 彭真明, 王峻, 等. 2011. 地震信号分数阶Gabor变换谱分解方法及应用[J]. 地球物理学报, 54(3): 867-873. |
董星宏, 李西京, 张国强, 等. 2009. 地磁记录图数字化识别的研究[J]. 地震地磁观测与研究, 30(6): 49-55. |
蒋一然, 宁杰远. 2019. 基于支持向量机的地震体波震相自动识别及到时自动拾取[J]. 地球物理学报, 62(1): 361-373. |
徐鹏深, 滕云田, 于子叶, 等. 2018. 基于信号指纹的地磁异常识别算法[J]. 地震学报, 40(1): 79-88. |
易世华, 刘代志, 何元磊, 等. 2013. 变化地磁场预测的支持向量机建模[J]. 地球物理学报, 56(1): 127-135. |
张秀玲, 柳正. 2017. 太阳活动周磁暴与亚洲MS7.0以上地震活动关系[J]. 地震地磁观测与研究, 38(5): 57-61. |
Dennis Gabor. 1946. Theory of communication[J]. Journal Institute of Electrical Engineers, 93: 429-457. |