随着城市建设的迅速发展、城市圈的快速扩张，受轻轨地铁、行车干扰、高压直流干扰的地磁台站逐年增多。为有效抑制干扰，邱颖^[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)来建立模型：

$ \mathop {{\rm{min}}}\limits_w = \sum\limits_{j = 1}^i | d(j) - {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{\varphi }}(j){|^2} + \lambda {\left\| \mathit{\boldsymbol{w}} \right\|^2} $

(1)

引入$\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{D}}(i) = {{[d(1), \cdots , d(i)]}^{\rm{T}}}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i) = [\mathit{\boldsymbol{\varphi }}(1), \cdots , \mathit{\boldsymbol{\varphi }}(i)]} \end{array}} \right.$，则最小二乘解为：

$ \mathit{\boldsymbol{w}}(i) = {[\lambda \mathit{\boldsymbol{I}} + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i)\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{(i)^{\rm{T}}}]^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i)d(i) $

(2)

式中，λ为小正数。通过矩阵求逆引理可以推出：

$ \mathit{\boldsymbol{w}}(i) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i){[\lambda \mathit{\boldsymbol{I}} + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{(i)^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i)]^{ - 1}}d(i) $

(3)

由于φ (u(i))维度较高，故构造核函数k(u(i), u(j))= φ (u(i)) φ (u(j))= φ (i) φ (j)，从而Φ ^T Φ可通过核空间来计算。通过上式可以看出，权重被明确表示为输入数据的一个线性组合：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{w}}(i) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i)\mathit{\boldsymbol{a}}(i)}\\ {\mathit{\boldsymbol{a}}(i) = {{[\lambda \mathit{\boldsymbol{I}} + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{{(i)}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(i)]}^{ - 1}}d(i)} \end{array}} \right. $

(4)

定义Q (i)=[λ I + Φ (i)^T Φ (i)]^-1, 使用上式滑动窗口结构，加入新数据后的矩阵的逆矩阵可以表示为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}(i) = \mathit{\boldsymbol{r}}{{(i)}^{ - 1}}}\\ {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}(i - 1)\mathit{\boldsymbol{r}}(i) + \mathit{\boldsymbol{z}}(i)\mathit{\boldsymbol{z}}{{(i)}^{\rm{T}}}}&{ - \mathit{\boldsymbol{z}}(i)}\\ { - \mathit{\boldsymbol{z}}{{(i)}^{\rm{T}}}}&1 \end{array}} \right]} \end{array} $

(5)

其中，$\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{z}}(i) = \mathit{\boldsymbol{Q}}(i - 1)\mathit{\boldsymbol{h}}(i)}\\ {\mathit{\boldsymbol{r}}(i) = \lambda + \mathit{\boldsymbol{\varphi }}{{(i)}^{\rm{T}}}\mathit{\boldsymbol{\varphi }}(i) - \mathit{\boldsymbol{z}}{{(i)}^{\rm{T}}}\mathit{\boldsymbol{h}}(i)}\\ {\mathit{\boldsymbol{h}}(i) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{{(i - 1)}^{\rm{T}}}\mathit{\boldsymbol{\varphi }}(i)} \end{array}} \right.$。因此，式(4)权重的扩展系数为：

$ \mathit{\boldsymbol{a}}(i) = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{a}}(i - 1) - \mathit{\boldsymbol{z}}(i)\mathit{\boldsymbol{r}}{{(i)}^{ - 1}}e(i)}\\ {\mathit{\boldsymbol{r}}{{(i)}^{ - 1}}e(i)} \end{array}} \right] $

(6)

式中，e(i)为通过期望信号与预测信号f_i-1(u(i))之差计算的预测误差：

$ {{f_{i - 1}}(u(i)) = \sum\limits_{j = 1}^{i - 1} {{\mathit{\boldsymbol{a}}_j}} (i - 1)k(u(j), u(i))} $

(7)

$ {e(i) = d(i) - {f_{i - 1}}(u(i))} $

(8)

针对地磁易受脉冲干扰信号、仪器本身和环境噪声持续干扰的特性，因此需具有逼近能力且数值计算稳定、同时能给出合理结果的核函数。研究中通常选择高斯核，其定义见文献[9]。

针对一维地磁输入，高斯核参数估计值在[h_s/10, 10h_s]范围内选取最为合适^[9-10]，其中h_s=1.06 min{σ, R/1.34}N^-1/5，σ为输入值u的标准差，R为u的四分位差。

在设计自适应滤波器的过程中，可以根据某种准则或学习规则逐步向数据集中加入样本。假设当前数据集为$C(i) = \{ {c_j}\} _{j = 1}^{{m_i}}$，c_j为第j个数据集，m_i为集的基数，当新数据对{u(i+1), d(i+1)}出现时，立即决策u(i+1)是否加入该集合，然后根据集合设置不同的阈值，通过阈值来调节系数以跟踪非平稳性。

基于该规则，系统根据上述算法可进行更新：

1) 加入新数据对{u(i+1), d(i+1)}，计算下列值：

$ \left\{ \begin{array}{l} C(i + 1) = \{ c(i), u(i + 1)\} \\ \mathit{\boldsymbol{h}}(i + 1) = {[k(u(i + 1), {c_1}), \cdots , k(u(i + 1), {c_i})]^{\rm{T}}}\\ {f_i}(u(i + 1)) = \mathit{\boldsymbol{h}}{(i + 1)^{\rm{T}}}\mathit{\boldsymbol{a}}(i)\\ e(i + 1) = d(i + 1) - {f_i}(u(i + 1))\\ \mathit{\boldsymbol{r}}(i + 1) = \lambda + k(u(i + 1), u(i + 1)) - \mathit{\boldsymbol{h}}{(i + 1)^{\rm{T}}}\mathit{\boldsymbol{Q}}(i)\mathit{\boldsymbol{h}}(i + 1)\\ \mathit{\boldsymbol{z}}(i + 1) = \mathit{\boldsymbol{Q}}(i)\mathit{\boldsymbol{h}}(i + 1)\\ \mathit{\boldsymbol{Q}}(i + 1) = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}(i) + \mathit{\boldsymbol{z}}(i + 1)\mathit{\boldsymbol{z}}{{(i + 1)}^{\rm{T}}}/\mathit{\boldsymbol{r}}(i + 1)}&{ - \mathit{\boldsymbol{z}}(i + 1)/\mathit{\boldsymbol{r}}(i + 1)}\\ { - \mathit{\boldsymbol{z}}{{(i + 1)}^{\rm{T}}}/\mathit{\boldsymbol{r}}(i + 1)}&{1/\mathit{\boldsymbol{r}}(i + 1)} \end{array}} \right] \end{array} \right. $

2) 计算并存储系数：

$ {\mathit{\boldsymbol{a}}_i}(i + 1) = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{a}}(i) - \mathit{\boldsymbol{z}}(i + 1)e(i + 1)/\mathit{\boldsymbol{r}}(i + 1)}\\ {e(i + 1)/\mathit{\boldsymbol{r}}(i + 1)} \end{array}} \right] $

更新规则后，如果预测误差小，则修正量也会很小，误差大修正量也会相应增大。在极端情况下，如果出现异常信号，系统将学习给出预测值。但在该情况下学习系统将会不稳定，因为系统一方面要学习原曲线形态，另一方面还需根据自身学习规则进行学习。尤其是针对受高压直流干扰具有台阶的地磁数据，当垂直台阶差较大时，数据可能会失真，因此如何制定学习规则尤为重要。

2 实验对象与结果分析

地磁台站每天所产出的地磁数据包含平静变化信号和叠加在平静变化上的不规则扰动变化信号，其中不规则扰动主要包含磁扰和脉冲2大类，其频率分布较广，成分较为复杂，频率分布约为10^-3~10² Hz^[4]。通过文献[11]给出的频谱分析发现，受轨道交通影响的干扰频带主要为0.001~0.020 Hz。在噪声优势频段内，距离干扰源越近，振幅谱突跳点越多，跳点幅度越大，其中地磁Z分量对干扰源最为敏感。本文以2018-06-07九峰地磁野外测试FGM01数据的Z分量为例，由于该野外测点受到的干扰源包含轻轨行车干扰和仪器自身干扰，其中噪声信息未知，并且部分时段未受干扰，因此对台站数据进行处理具有参考价值。同时，选取2019-07-21受高压直流影响的应城地磁台FGM01数据，对该算法进行验证。

从图 1(b)可以看出，仪器自身正常噪声在0.2 nT以下，干扰噪声叠加其上，噪声分布不均匀；从频率图(图 1(c))可以看出，噪声频率布满整个频率带，并未凸显优势频率；从图 1(d)可以看出，最大噪声频段集中在0.010 Hz左右，表明野外观测数据的干扰源较多，同时还相互叠加，呈现高度非线性。

为了根据台站仪器特性制定相应的学习规则，本文采用未受人为干扰的应城台数据进行地磁秒数据理论分析，其中GM3数据Z分量的原始曲线和一阶差分曲线如图 2所示。

从应城台6 d的原始曲线与一阶差分图(图 2)可以看出，未受干扰的Z分量背景噪声基本稳定在±0.2 nT以内，不排除其他干扰，少部分在±0.4 nT以内，整体不超过±0.8 nT。因此可将本次学习规则设为：

$ {c_j} = \left\{ \begin{array}{l} {\rm{正常区}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} |du| \le 0.2\\ {\rm{接受区}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.2 < |du| \le 0.4\\ {\rm{勉强接受区}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.4 < |du| \le 0.8\\ {\rm{异常区}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.8 < |du| \end{array} \right. $

其中, du=u(i+1)-u(i), 2组数据一阶差分统计与初始化参数如表 1、2所示。

从2个表中一阶差分统计结果可以看出，绝大部分数据都在±0.2以内，少量分布于其他区域。为简化计算，本文设定学习规则为：在±0.2以内为正常区，参考值设为1表示采纳；将其他区域参数设为0表示拒绝，拒绝区将学习给出估计值。

针对具有明显垂直梯度的曲线，如方波曲线(图 3)，如果只对0.2以内为正常接收的学习规则，滤波学习后在上沿和下沿之间会存在一段以0.2为斜率的折线。为了改善该类失真情况，采用一阶差分学习方法。从图 3可以看出，曲线存在1个突跳点，利用规则滤掉突跳点后一阶差分为一条直线，然后归算至曲线，原曲线的上沿将下移至下沿而形成完整的曲线，从而克服台阶问题。

采用上述规则学习相同时段恩施地磁台数据并将结果作为期望值，滤波后结果与频谱分析如图 4所示。

从图 4中滤波曲线与一阶差分曲线可以看出，噪声值降低直至正常区，随着时间的推移加载在曲线上的噪声幅度减小。通过频谱分析发现，在0.05 Hz以后，曲线表现为有序、逐渐降低的形态，说明此后数据基本不受干扰或干扰较小，表明该算法具有自主学习能力并能不断优化滤波参数。而6~12 h时段存在噪声，这是因为该算法需要一个学习适应的过程。从图 5可以看出，通过对一阶差分进行修正，然后反算至原始数据可以消除台阶干扰，这将为台站自动消除数据台阶提供新方法。目前针对高压直流输电干扰，主要做法是联合电力部门进行人工抑制。这2种方法会在台阶断点处略有差距，但走势一致，不会对地磁数据后期分析造成影响。

为进一步验证算法的实用性与可靠性，本文选取2019-09-24~09-30受轻轨、地铁与施工干扰的九峰地磁台GM4磁力仪H分量数据作为实验对象，为分辨滤波曲线，将同时间段内滤波曲线减小50 nT(图 6)。从图中可以看出，前期学习期间仍存在干扰现象，但随着时间推移，从25日(即86 400 s)后曲线形态开始恢复正常，并能正常反映地磁场变化。上述分析表明，该算法泛化学习能力强，具有其他算法不具备的优越性。

3 结语

通过对受干扰的地磁数据进行频谱分析发现，轨道交通对地磁信号的噪声优势频段为0.005~0.020 Hz，主要干扰频段为0.010~0.020 Hz。轻轨交通、仪器自身噪声、行车干扰若同时加载在地磁信号上则无优势频段，呈现出非线性。基于地磁干扰特征与频谱分析，本文提出基于最小二乘的核自适应滤波算法，并制定相应的学习规则，推导递推公式。

通过实际干扰数据分析可以看出，核自适应滤波算法作为一种非线性滤波方法，能通过实时调节滤波参数从而抑制干扰，为一种完全自适应的无需人工干预的滤波算法。相比于其他算法，该算法对于地磁数据的处理更具有适用性。本文首次将核自适应滤波算法应用于地磁干扰抑制研究，以期为地磁干扰抑制提供新方法。