探讨地磁秒数据干扰抑制中核自适应滤波算法的运用随着城市建设的迅速发展、城市圈的快速扩张,受轻轨地铁、行车干扰、高压直流干扰的地磁台站逐年增多为有效抑制干扰,邱颖[1]采用小波技术对其进行抑制,万永革等[2]利用傅里叶和相关分析2种方法对干扰进行抑制,郑江蓉等[3]提出基于极值特征尺度的干扰辨识方法,谢凡等[4,5]提出小波阈值去噪模型与基于图像处理技术的数学形态学滤波算法,李季、向阳等[6,7]采用基于HHT变换域的局部自适应方法进行研究然而这些算法只对具有明显特征性的干扰源具有良好的抑制效果,当面对多种干扰源同时出现的情况,如何不失真地滤除干扰噪声将尤为重要目前,智能算法越来越成熟,本文将探索在自动获取数据的同时可自主调节滤波参数以使滤波算法最优化,该类算法可更好地应用在台站数据处理中1、核自适应滤波模型核自适应滤波器的主要思想为:用一个线性自适应滤波算法来实现非线性信号处理,同时通过选择一个固定的映射,将线性输入数据映射到高维的特征空间中,然后在该空间中使用线性滤波算法核自适应可提供一个通用的自适应滤波器,在滤波器权重中嵌入一个持续增长的存储结构,采用径向基函数网,通过数据学习网络拓扑并调整滤波参数,实现随环境变化而实时调整参数的功能[8,9]。
假设目标是在输入-输出样本序列{u(1),d(1)},{u(2),d(2)},…,{u(N),d(N)}的基础上学习一个连续的输入-输出映射,考虑到地磁数据特点,设输出为一维输出为能够学习任意的非线性映射关系,将输入u(i)变换到高维特征空间F中φ(u(i)),在高维特征空间中通过线性运算来求解低维空间的非线性滤波为简化计算,令φ(i)=φ(u(i)),在新的样本序列{φ(i),d(i)}中通过最小化代价函数估计权重w(i)来建立模型:式中,λ为小正数通过矩阵求逆引理可以推出:w(i)=Φ(i)[λI+Φ(i)TΦ(i)]−1d(i) (3)由于φ(u(i))维度较高,故构造核函数k(u(i),u(j))=φ(u(i))φ(u(j))=φ(i)φ(j),从而ΦTΦ可通过核空间来计算通过上式可以看出,权重被明确表示为输入数据的一个线性组合:{w(i)=Φ(i)a(i)a(i)=[λI+Φ(i)TΦ(i)]−1d(i) (4)定义Q(i)=[λI+Φ(i)TΦ(i)]-1,使用上式滑动窗口结构,加入新数据后的矩阵的逆矩阵可以表示为:Q(i)=r(i)−1[Q(i−1)r(i)+z(i)z(i)T−z(i)T−z(i)1] (5)其中,。
因此,式(4)权重的扩展系数为:a(i)=[a(i−1)−z(i)r(i)−1e(i)r(i)−1e(i)] (6)式中,e(i)为通过期望信号与预测信号fi-1(u(i))之差计算的预测误差:fi−1(u(i))=∑j=1i−1aj(i−1)k(u(j),u(i)) (7)e(i)=d(i)−fi−1(u(i)) (8)针对地磁易受脉冲干扰信号、仪器本身和环境噪声持续干扰的特性,因此需具有逼近能力且数值计算稳定、同时能给出合理结果的核函数研究中通常选择高斯核,其定义见文献[9]针对一维地磁输入,高斯核参数估计值在[hs/10,10hs]范围内选取最为合适[9,10],其中hs=1.06min{σ,R/1.34}N-1/5,σ为输入值u的标准差,R为u的四分位差在设计自适应滤波器的过程中,可以根据某种准则或学习规则逐步向数据集中加入样本假设当前数据集为C(i)={cj}mij=1,cj为第j个数据集,mi为集的基数,当新数据对{u(i+1),d(i+1)}出现时,立即决策u(i+1)是否加入该集合,然后根据集合设置不同的阈值,通过阈值来调节系数以跟踪非平稳性。
基于该规则,系统根据上述算法可进行更新:1)加入新数据对{u(i+1),d(i+1)},计算下列值:2)计算并存储系数:ai(i+1)=[a(i)−z(i+1)e(i+1)/r(i+1)e(i+1)/r(i+1)]更新规则后,如果预测误差小,则修正量也会很小,误差大修正量也会相应增大在极端情况下,如果出现异常信号,系统将学习给出预测值但在该情况下学习系统将会不稳定,因为系统一方面要学习原曲线形态,另一方面还需根据自身学习规则进行学习尤其是针对受高压直流干扰具有台阶的地磁数据,当垂直台阶差较大时,数据可能会失真,因此如何制定学习规则尤为重要2、实验对象与结果分析地磁台站每天所产出的地磁数据包含平静变化信号和叠加在平静变化上的不规则扰动变化信号,其中不规则扰动主要包含磁扰和脉冲2大类,其频率分布较广,成分较为复杂,频率分布约为10-3～102Hz[4]通过文献[11]给出的频谱分析发现,受轨道交通影响的干扰频带主要为0.001～0.020Hz在噪声优势频段内,距离干扰源越近,振幅谱突跳点越多,跳点幅度越大,其中地磁Z分量对干扰源最为敏感本文以2018-06-07九峰地磁野外测试FGM01数据的Z分量为例,由于该野外测点受到的干扰源包含轻轨行车干扰和仪器自身干扰,其中噪声信息未知,并且部分时段未受干扰,因此对台站数据进行处理具有参考价值。
同时,选取2019-07-21受高压直流影响的应城地磁台FGM01数据,对该算法进行验证从图1(b)可以看出,仪器自身正常噪声在0.2nT以下,干扰噪声叠加其上,噪声分布不均匀;从频率图(图1(c))可以看出,噪声频率布满整个频率带,并未凸显优势频率;从图1(d)可以看出,最大噪声频段集中在0.010Hz左右,表明野外观测数据的干扰源较多,同时还相互叠加,呈现高度非线性为了根据台站仪器特性制定相应的学习规则,本文采用未受人为干扰的应城台数据进行地磁秒数据理论分析,其中GM3数据Z分量的原始曲线和一阶差分曲线如图2所示图12018-06-07FGM01Z分量原始曲线及频谱分析图22019-07-12～07-18应城台GM3Z分量原始曲线与一阶差分曲线从应城台6d的原始曲线与一阶差分图(图2)可以看出,未受干扰的Z分量背景噪声基本稳定在±0.2nT以内,不排除其他干扰,少部分在±0.4nT以内,整体不超过±0.8nT因此可将本次学习规则设为:cj正常区 |du|≤0.2接受区 0.2<|du|≤0.4勉强接受区 0.4<|du|≤0.8异常区 0.8<|du|其中,du=u(i+1)-u(i),2组数据一阶差分统计与初始化参数如表1、2所示。
表1野外勘测数据Z分量一阶差分范围统计与参数选定表2受高压直流干扰的应城FGM01Z分量一阶差分统计与参数选定从2个表中一阶差分统计结果可以看出,绝大部分数据都在±0.2以内,少量分布于其他区域为简化计算,本文设定学习规则为:在±0.2以内为正常区,参考值设为1表示采纳;将其他区域参数设为0表示拒绝,拒绝区将学习给出估计值针对具有明显垂直梯度的曲线,如方波曲线(图3),如果只对0.2以内为正常接收的学习规则,滤波学习后在上沿和下沿之间会存在一段以0.2为斜率的折线为了改善该类失真情况,采用一阶差分学习方法从图3可以看出,曲线存在1个突跳点,利用规则滤掉突跳点后一阶差分为一条直线,然后归算至曲线,原曲线的上沿将下移至下沿而形成完整的曲线,从而克服台阶问题图3方波曲线滤波前后对比采用上述规则学习相同时段恩施地磁台数据并将结果作为期望值,滤波后结果与频谱分析如图4所示从图4中滤波曲线与一阶差分曲线可以看出,噪声值降低直至正常区,随着时间的推移加载在曲线上的噪声幅度减小通过频谱分析发现,在0.05Hz以后,曲线表现为有序、逐渐降低的形态,说明此后数据基本不受干扰或干扰较小,表明该算法具有自主学习能力并能不断优化滤波参数。
而6～12h时段存在噪声,这是因为该算法需要一个学习适应的过程从图5可以看出,通过对一阶差分进行修正,然后反算至原始数据可以消除台阶干扰,这将为台站自动消除数据台阶提供新方法目前针对高压直流输电干扰,主要做法是联合电力部门进行人工抑制这2种方法会在台阶断点处略有差距,但走势一致,不会对地磁数据后期分析造成影响图4野外勘测数据滤波与频谱分析图5受高压直流干扰的应城台滤波前后Z分量对比与滤波后一阶差分曲线为进一步验证算法的实用性与可靠性,本文选取2019-09-24～09-30受轻轨、地铁与施工干扰的九峰地磁台GM4磁力仪H分量数据作为实验对象,为分辨滤波曲线,将同时间段内滤波曲线减小50nT(图6)从图中可以看出,前期学习期间仍存在干扰现象,但随着时间推移,从25日(即86400s)后曲线形态开始恢复正常,并能正常反映地磁场变化上述分析表明,该算法泛化学习能力强,具有其他算法不具备的优越性图62019-09-24～09-30九峰地磁台GM4磁力仪原始曲线与滤波后曲线对比3、结语通过对受干扰的地磁数据进行频谱分析发现,轨道交通对地磁信号的噪声优势频段为0.005～0.020Hz,主要干扰频段为0.010～0.020Hz。
轻轨交通、仪器自身噪声、行车干扰若同时加载在地磁信号上则无优势频段,呈现出非线性基于地磁干扰特征与频谱分析,本文提出基于最小二乘的核自适应滤波算法,并制定相应的学习规则,推导递推公式通过实际干扰数据分析可以看出,核自适应滤波算法作为一种非线性滤波方法,能通过实时调节滤波参数从而抑制干扰,为一种完全自适应的无需人工干预的滤波算法相比于其他算法,该算法对于地磁数据的处理更具有适用性本文首次将核自适应滤波算法应用于地磁干扰抑制研究,以期为地磁干扰抑制提供新方法参考文献：[1]邱颖.地电场观测中已知源干扰抑制研究[D].北京:中国地震局地震预测研究所,2008.[2]万永革,盛书中,曾献军,等.地磁与地电混合观测台地磁的地电测数干扰处理方法[J].地震,2009,29(2):81-87.[3]郑江蓉,谢凡,张秀霞,等.基于极值特征尺度的干扰辨识技术在地磁观测中的应用研究[J].地球物理学进展,2012,27(2):528-536.[4]谢凡,滕云田,胡星星.数学形态滤波在地磁数据干扰抑制中的应用[J].地球物理学进展,2011,26(1):147-156.[5]谢凡.地磁观测数据中人工电磁干扰抑制技术研究[D].北京:中国地震局地球物理研究所,2011.[6]李季,潘孟春,唐莺,等.基于形态滤波和HHT的地磁信号分析与预处理[J].仪器仪表学报,2012,33(10):2175-2080.[7]向阳,孙吉泽,赖爱京,等.HHT变换在地磁数据处理的应用[J].地震地磁观测与研究,2015,36(6):48-52.[9]陈乾.核自适应滤波算法研究[D].武汉:华中师范大学,2014.[11]张明东,马骥,尚先旗,等.天津GM4磁通门磁力仪受轨道交通干扰的谱分析[J].内陆地震,2015,29(2):154-161.罗棋,周洋,吴艳霞,李查玮.核自适应滤波算法在地磁秒数据干扰抑制中的研究[J].大地测量与地球动力学,2020,40(07):708-712.基金:中国地震局地震研究所基本科研业务费专项(IS201866292-5).。