0 引言

地球磁场包含内源场和外源场，其中外源场起源于空间电流体系，主要分布于电离层和磁层中(徐文耀，2003).外源变化场及其在地球内部产生的感应场虽然仅约占地磁总场的1%，但它包含了丰富的地球物理信息，如地球磁层、电离层动态结构(Yamazaki et al., 2009; 袁桂平等，2015)及地下介质电导率等信息(汤吉等，2005；柳建新等，2012；王桥和黄清华，2016)，甚至可能还含有与地震孕育相关的信息(张继红等，2010；赵国泽等，2010；丁鉴海等，2011；Xu et al., 2013; Han et al., 2016, 2017).因此，地磁场观测在空间天气监测、矿产资源勘探、自然灾害研究等领域具有重要价值.然而，随着城市建设的不断发展，地磁场观测受到的人文干扰越来越严重.这些汽车、地铁、高铁等产生的近场噪音造成了地磁脉动等变化场及其感应场难以被精确观测，严重影响了地磁数据的实际使用价值.因而，抑制这些近场噪音在电磁法勘探、空间天气及地震电磁研究中均具有重要意义，一直备受关注(Larsen et al., 1996; Egbert, 1997; Chave and Thomson et al., 2004; Takahashi et al., 2007; 汤井田等，2008；Hattori et al., 2013a, 2013b；Hattori and Han, 2018；Han et al., 2014, 2015；凌振宝等，2016).

国内外的研究学者就抑制地磁信号中的噪声干扰提出了一系列的方法，其中，自然正交分量法(Xu and kamide, 2004; 顾左文等，2009)和与其类似的主成分分析法(韩鹏等，2009; Huang, 2011)在提取地磁场主要变化特征中有着较好的应用，但由于这两种方法本质上都是观测样本的正交能量投影，并未对各观测样本进行剖分，因而当磁场信号和噪音信号在观测样本间相关性均较高时，去噪效果不理想.此外，短时傅里叶、小波变换等时频分析工具也被应用到地磁数据的去噪中(Harada et al., 2004; 谢凡，2009; 谢凡等，2011；Han et al., 2011; 汪伟明和贺巍，2018)，可在一定程度上抑制轨道交通、汽车等人文干扰，但此类方法多是基于磁场信号和噪音信号的频率特征不同，通常对频率较高的噪音信号进行滤除，因而去除噪声后的地磁信号的高频部分损失严重，容易丢失相关频段信号包含的信息.而且，当信号和噪音频段混叠时，这些时频分析方法较难奏效.汤井田等(2012)提出了一种基于信号形态特征的数学形态滤波法，可对非线性、非平稳近场噪声进行有效压制，但如何选取合适的滤波器结构元素及其尺寸有待进一步研究.也有学者尝试将奇异谱分析法应用于电车噪音去除中(Saito et al., 2011, 2013)，实现了对噪音信号的准确检测，但由于奇异谱分析方法只能反映信号时间上和空间上的线性相关关系，在剥离噪音信号中仍存在较大困难.综上，如何有效抑制/去除地磁场观测中的人文噪音仍是一个尚未彻底解决的问题.

盲源分离指仅从若干观测到的混合信号中恢复出无法直接观测的各个原始信号的过程(Jutten and Herault, 1991).这里的“盲”，是指源信号不可测，混合系统特性事先未知这两个方面.盲源分离是信号处理中一个传统而又极具挑战性的问题，通常需要假设各源信号具有某些统计特征，然后基于特定的目标函数和优化准则实现源信号和混合矩阵的估计，从而在信号源及混合机制均未知的情况下实现不同场源信号的分离(杨福生和洪波，2006).盲源分离法自提出以来，在语音信号处理(苏野平和杨荣震，2002)、地球物理信号处理(贾瑞生等，2015)、医学信号处理(Jung et al., 2000)等领域得以成功应用.由于地磁场观测数据包含空间场源、近场噪音等多个场源，因此盲源分离在地磁场数据去噪中也可能具有潜在应用价值.谢凡等(2011)将盲源分离法应用到地磁场轨道交通干扰的抑制中，取得了较好效果.但由于混合信号场源种类的不确定性、含噪声干扰的复杂性以及磁场信号和噪声干扰混合机制的不确定性，当噪音信号和空间场源信号特征相似时，分解结果存在一定的不稳定(Snoussi et al., 2007).目前，盲源分离法在地磁场数据去噪中的应用尚不成熟，仍处于探索阶段.本文基于已有研究，结合远场参考法，尝试寻找一种稳定有效的分离近场噪声的方法.首先，我们开展仿真分析，阐明此方法的原理并考察其有效性.在此基础上，我们将该方法应用于日本千叶地区实际磁场观测数据，并选用KAK国际基准台作为远场参考，进一步验证其在实际数据分析及近场噪音抑制中的效用.

1 盲源分离法 1.1 基本原理

盲源分离在20世纪80年代由法国学者JeannyHerault和Christian Jutten提出，他们首先基于递归神经网络模型和学习算法，实现了2个独立源信号混合的分离(Herault and Jutten, 1986).盲源分离研究的核心问题是如何在源信号和混合系统特性无法精确获知的情况下，从混迭信号(观测信号)中分离出各源信号.它的基本瞬时线性模型可用下式表示：

(1)

式中X(t)为观测信号，假设已知m个观测样本，则X(t)=[x₁(t), x₂(t), …x_m(t)]^T，S(t)为源信号，可表示为：S(t)=[s₁(t), s₂(t), …s_n(t)]^T.一般情况下观测样本个数应不小于源个数，即m≥n.A为m×n混合矩阵，表征观测系统的混合特征.本研究使用最常用的线性混合模型，各源信号S(t)经过不同线性组合叠加得到观测样本信号X(t)，线性组合系数矩阵即为A.

当观测样本数和源的个数相同时，通过对(1)式中的混合矩阵A求逆，再将求得的逆矩阵左乘观测信号X(t)，便可得源信号S(t).通常这个过程用一个盲源分离的求解模型来表示：

(2)

其中，Y(t)为分离信号，是源信号的估计Y(t)=Ŝ(t)，W是分离矩阵.理论上讲，分离矩阵具有任意性，式(2)有无穷多个解.但当我们假设源信号具有某些统计特征时，比如，各源信号彼此统计独立时，理论上可求得唯一最优解.这种基于信号间统计独立性求解的方法被称为独立成分分析法，是一种最常见的盲源分离方法.独立成分分析法原理简单明了，且源信号独立性假设合乎绝大多数观测系统的客观规律，其在信号处理领域得到了广泛的应用(Comon and Jutten, 2010).因而本研究也使用独立成分分析法进行磁场数据的盲源分离.

1.2 算法实现

在众多基于独立成分的判定准则中，负熵以其稳健性强、计算效率较高等优点在盲源分离中得到了越来越多的应用.负熵与信号的非高斯性等价，非高斯性越强则负熵越大.根据概率统计论中的中心极限定理，由多个独立源信号组成的混合信号的高斯性强于任何一个独立信号，换句话说，单个源信号的非高斯性强于混合源信号的非高斯性.因此，在基于负熵的独立成分算法中，只要判定出混合信号中某个信号拥有最大的负熵值，即可认定该信号为一个源信号.目前比较成熟的FastICA算法(Hyvarinen, 1999)，即是通过判定信号的负熵来实现混合信号分离.

在上述算法中，负熵可定义为

(3)

(4)

式中，H(·)为信号的熵，J(·)为负熵，p(·)是随机变量的概率密度函数，Y_gauss表示同Y具有相同方差的高斯随机变量.在FastICA算法中，由于概率密度函数难以估计，进而采用了非线性函数代替概率密度函数的近似方式来进行优化判定.近似负熵公式为

(5)

上式中，E(·)是均值运算，G(·)代表非线性函数，一般地，可根据信号的高斯性选择不同的非线性函数代入到式(5)中进行运算(Hyvarinen, 1999)，常用非线性函数有

(6)

(7)

(8)

本文选用反映信号高阶统计特征及非高斯性的式(8)作为非线性函数(Hyvarinen, 1999).将式(2)及(8)代入到式(5)中得到式(9)：

(9)

在上式中，只需要求得一个分离矩阵使得J(Y)最大便可得到分离信号.一般采用牛顿迭代法搜索J(Y)极值，进而求得分离矩阵.

需要特别指出的是，从式(1)可以看出，由于混合矩阵和各个相互独立的源都是未知的，如果对独立源乘上某个标量，同时对混合矩阵除以一个相同的标量，混合信号仍保持不变.因此，在盲源分离中，源信号的方差、能量是不能确定的.为了避免解的不唯一性，提高优化搜索的计算效率，在牛顿迭代搜索求解中对源信号进行方差归一化处理，因而得到的源通常存在幅度失真.另一方面，从公式(1)还可以看出，如果把混合矩阵的列向量重新排序，源信号的行向量做相对应的顺序调整，混合信号不变.因此，在盲源分离中，各独立源信号的顺序不能确定.综上，基于信号独立性的盲源分离给出的是“波形保持解”，只能提供源信号的波形特征，不能提供信号绝对强度信息.

2 仿真分析

为了进一步阐明盲源信号分离的工作原理，并考察其有效性，我们首先进行仿真分析.根据实际观测，地磁场数据中通常包含两种场源：第一种为远场空间信号，源于电离层及磁层，在较大空间范围内信号具有很好一致性；第二种为近场人文干扰，其中最常见的是电车、高铁、地铁等干扰，是一种局部性信号.基于此，我们构造了如图 1所示的两个源信号：(a)周期为20 s的调幅正弦波信号，模拟地磁连续性脉动Pc3；(b)周期为19 s的方波信号，模拟电车干扰信号.经如式(10)所示的矩阵混合后，加入1%的随机白噪音使其更贴近实际观测信号，得到如图 2所示的混合信号X₁.

(10)

将混合信号X₁按照1.2节中的算法进行盲源分离，得到了图 3所示的分离信号Y₁，及分离矩阵：

(11)

对比图 1中源信号可知，分离信号Y₁较好的刻画了源信号的波形特征，但随机噪音相对幅度似乎有所放大.图 3b中正弦波分离信号仍包含着部分方波信号，与图 1中的正弦波信号存在一定差异.为了进一步提升分离效果，我们尝试使用参考台站进行盲源分离.

在实际地磁数据中，远场参考台通常认为是不含有近场噪音.因此我们对混合矩阵A₁稍作修改，将其中一个混合信号的方波输入置零，得混合矩阵A₂：

(12)

然后对图 1中的两个信号源进行混合，加入1%的随机白噪音使其更贴近实际观测信号，得到如图 4所示的混合信号X₂.图 4a中的信号不含有方波噪音，模拟参考台信号；图 4b中的信号与图 2b中的信号相同，模拟目标台站信号.

考虑到参考台站不含有近场噪音，其本身可视为一个源信号，因此参考台站信号应和分离所得源信号之一具有较高相似度.基于此，对混合信号X₂进行盲源分离时，我们可进行进一步的监督学习，即：在考察各个信号源与参考台信号相似度的基础上，搜索负熵极值，获得同时满足独立性和与参考台相关性的最优解.基于上述思想，对混合信号X₂进行盲源分离，得到了图 5所示的分离信号Y₂，及分离矩阵：

(13)

对比图 1及图 3可发现，当使用参考台时，分离效果明显提升，可较准确地得到源信号的波形信息.但无论是否使用参考台，分离所得信号与源信号均存在幅度差异，且还可能存在180°相位差(图 5)，这在其对应的分离矩阵W₂中亦有所体现，即对应的源信号混合系数为负值.

将式(2)两侧乘以分离矩阵W，可重构混叠信号(观测信号)：

(14)

基于此式可以进一步发现，虽然独立源Y乘上某个标量，同时对混合矩阵W除以一个相同的标量，重构信号保持不变，故源的强度不能确定，但各独立源与其对应的混合系数的乘积是确定的，即对于每一个重构的混合信号(观测样本)，各独立源的贡献是唯一的：

(15)

式中为重构信号的行向量，对应各观测样本，i=1, 2, …, m，Y_j(t)为盲源分离所得独立源，n为独立源个数.通常，源的个数与观测样本数相同，即m=n.基于参考台数据信息，将噪音源的贡献置零，即可实现目标台站的近场噪音剥离.

图 6展示了X₂中的目标台站重构信号与原始混合信号的对比.其中(a)为正弦波信号在原始混合信号(黑色)及在重构信号中的贡献(红色)；(b)为方波信号在原始混合信号中(黑色)及在重构信号中的贡献(红色)；(c)为原始混合信号(黑色)，重构信号(红色)，及将重构信号中的方波信号置零所得的去噪信号(绿色).从中可以发现，基于参考台的盲源信号分离法可完全实现信号重构，并能很好的将混合信号(观测样本)中包含的不同源的成分分开，从而准确剥离目标台站干扰信号，提取观测数据中的空间场源信号.

基于以上仿真分析可知，盲源信号分离法可有效分离不同种类的信号，可能在地磁场数据处理及去噪中存在应用价值.当使用参考台数据时，可提升信号分离质量，因此在接下来实测数据分析中，我们将使用参考台站数据.

3 实测数据处理

为了验证本文所提方法的实用性，我们选用日本千叶地区布设的FDG、KYS、UCU三个地磁观测台于2005年1月11日观测的一段长900 s的数据进行分析(台站详情及磁力仪参数见Hattori et al., 2004).同时选用距离此三个台站约120 km的国际基准台Kakioka(KAK，36.232°N, 140.186°E)作为远场参考.图 7给出了上述四个台站的空间分布.从图中可以看到，FDG、KYS、UCU三个台站间距离电车轨道较近，易受电车干扰.图 8展示了4个台站的南北分量的秒采样磁场实测数据(去均值后)，从中可以看到，四个台站的磁场变化比较一致，主要呈现的是地磁脉动信号特征，但位于千叶的三个台站包含明显的电车干扰信号(类方波信号和脉冲信号).基于已有的研究可知(Hattori et al., 2013a)，观测数据中的信号源主要有两类：空间场源和近场人文干扰.但每一类源又可能包含不同形态的信号，比如空间场源中可能含有连续性脉动(Pc)和不规则脉动(Pi), 人文噪音可能包含类方波信号和短时脉冲等.因此，在实际数据分析中，很难事先确定源的个数.作为一种尝试，我们将四个台站的观测数据作为观测样本矩阵，进行盲源信号分离.

分离结果如图 9所示，共得到四个相互独立的信号源：IC₁，IC₂，IC₃，IC₄，及分离矩阵W：

(16)

从图 9中可以发现，分离信号IC₃在形态与远参考台站KAK的观测信号高度相似，因此可判定其是空间场源信号.分离信号IC₁呈现出了原始观测信号中电车干扰的波形形态，且出现的类方波的时间点与原始观测信号中干扰出现的时间点基本一致，因此判定其为近场电车噪音.分离信号IC₂在第680~850 s间也存在较明显脉冲/类方波信号，因此也可判定其为近场噪音.分离信号IC₄较难判定属于何种信号，鉴于其包含的短时脉冲信号与IC₁中电车信号具有同步性，因此也将其归为近场噪音.

根据仿真分析可知，基于式(15)，通过分离矩阵和估计的源信号可重构观测信号(理论上重构的观测信号与原始观测信号一致)：

(17)

由于我们在盲源分离时使用的FastICA算法对源信号进行了方差归一化处理，因此分离矩阵W的系数可反映对应源在各观测样本中的贡献.考察式(16)，发现其第三列的数值较其他三列大，根据式(17)可知，W的第三列是IC₃的系数，说明空间场源信号在各观测样本中占据绝对能量贡献，即各台站磁场数据中主要成分为空间场源信号，此结果与图 8所展示的结果完全一致，从侧面印证了IC₃为空间场源信号.进一步考察分离矩阵W第三列数据可以发现，各系数差异很小，说明空间场源信号在各台站贡献较为一致，与已有的研究相符(Hattori et al., 2013b).

为了更直观阐明各观测样本包含的独立源成分，我们基于式(17)计算各重构观测信号中每一种源的能量贡献，结果如图 10所示.在该图中，两个横轴分别为独立源和台站观测样本，纵轴为源在台站观测样本中的能量贡献.首先，分析KAK台站各源能量贡献可知IC₃占据绝对多数，其他三个源在KAK台站的贡献均不足其1%.这与KAK台站为远场参考台，几乎不含有电车干扰的事实相符，同时也说明IC₁，IC₂，IC₄为近场噪音，他们对远场参考台的能量贡献非常小.此外，IC₃在各台站的能量贡献均占据绝对主导，证明其为全局性信号，这为判定IC₃是空间场源信号提供了有力依据.考察IC₁在各台站的能量贡献可知，其在参考台站KAK的能量贡献几乎为零，但在其他三个台站KYS、FDG、UCU的能量贡献显著，因此可确认IC₁为近场噪音的主要成分.综上，加入参考台的盲源信号分离可为源的性质确定提供重要参考信息.

基于上述分析，我们可以判定分离信号IC₁、IC₂、IC₄为近场噪声干扰，按照上节仿真分析中的方法，将它们的贡献置零，可得各台站去除电车干扰之后的信号.图 11给出了剔除电车干扰前后的各台站观测信号及两者之差.从(a)，(c)，(e)中可以看到，位于千叶的三个观测台站中的电车干扰信号得到了明显的抑制/剔除；而在(g)图中，参考台KAK原始信号与去噪信号几乎重合，且(h)图显示，他们间的差值小于0.03 nT，几乎接近仪器的分辨率，说明KAK台站包含近场噪音极小，从侧面证明了所使用分离方法的合理性和有效性.此外，(b)，(d)，(f)图所示的去噪前后观测数据差，实际上反映了各台站总体的近场噪音，其幅度在各台站间存在差异，KYS和UCU台站噪音信号幅度大致相当，FDG台站明显较弱.观察图 7台站分布可知，FDG台离主要噪音源电车轨道的距离最远，而KYS和UCU台离电车轨道距离大致相同，因此，当地下电导率大致均匀时，理论上FDG接受到的来源于电车轨道的信号应该较弱，KYS和UCU接收到的信号则应相对较强.这一推论与图(b)，(d)，(f)的观测事实吻合，表明主要的近场噪音源是电车，同时也说明了分离结果较为精确.

为了进一步考察去噪效果，我们对去噪前后的数据进行谱分析.图 12为各台站去噪前后的谱对比，其中参考KAK台的频谱几乎不变，而位于千叶的KYS、FDG、UCU台站的信号谱在0.02~0.2 Hz频段去噪后出现了明显的下降.而基于已有研究(Saito et al., 2011)，电车干扰多出现在这一频段，因此谱分析的结果也可以说明近场电车噪音得到了有效抑制，获得了较好的去噪效果.图 13将三个目标台站与参考台的谱进行了对比，(a)为去噪前数据，(b)为去噪后数据.从(a)中可以看出，位于千叶的三个目标台在0.02 Hz以上的频段相对参考台信号明显较强，且三个台的谱较为一致，说明这三个台站受到了较为一致的近场噪音干扰.经过盲源分离去噪后，三个台站的近场噪音均得到了较好的抑制，去噪后的谱与参考台的谱吻合较好，如(b)所示.

4 讨论

如何有效抑制近场人文噪音一直是地球电磁学领域的重要研究课题.传统的方法多基于时频分析，如短时傅里叶分析、小波分析等，因受海森堡测不准原理制约，这些基于母函数的分析方法难以同时在时间和频率域中达到高精度，对于突变信号如电车产生的方波干扰及脉冲噪音处理效果较差.此外，当信号和噪音频段混叠时，这些时频分析方法效果往往不理想.在本文的仿真分析中，方波信号的周期是19 s，正弦波信号周期是20 s，两者极为接近，但基于参考台的盲源信号分离法仍可很好的将两者分离，展示了此方法的有效性和潜在优势.

在盲源分离的实际应用中，如何确定源的种类和性质是一个难题，直接影响分离结果的进一步使用.本文结合磁场的空间场源信号在超低频(ULF)段较大范围内近似相同这一特性，将远场参考植入盲源分离法中，提升了分离效果.因参考台不含有近场噪音，其本身可视为一个已知源，可对分离结果进行进一步约束优化，因此引入参考台数据有助于确定分离结果中源的种类和性质.实测数据的处理结果表明，基于参考台的盲源分离法可有效分离不同场源信号，将分离所得信号与参考台信号对比，并考察各源能量贡献，可较明确判定源的种类，从而识别近场噪音源和空间场源，实现近场噪音的抑制/剥离.需要注意的是，参考台的数据质量及空间场源信号在台站间的一致性，会对盲源分离的结果产生影响，因此在实际数据处理中，找到合适的参考台站是使用本方法的关键.

图 11中，基于盲源分离结果，我们对信号进行了重构进而较精确的得到了各个台站的近场噪音.这些近场噪音不仅包含电车干扰源的特征，也包含源到观测台站地下电导率信息，因为近场噪音在从源(电车)到观测台站的传播过程中，其衰减特征与传播路径上的电导率密切相关.因此准确提取的近场噪音有可能在浅层电性结构的研究中具备使用价值.在未来的研究中，可尝试开发一套类似地震学中背景噪音的方法，从精确提取的电磁近场噪音中获取地下电性结构信息.此外，在大地电磁勘探中，近场噪音普遍存在，本方法可实现近场信号和空间场源信号的分离，能较准确的对振幅谱进行校正，在主动源及被动源勘探中均具备潜在的应用价值.

5 结论

本研究将远场参考法和盲源信号分离法相结合，提出了一种能有效抑制地磁场观测数据中近场噪音的方法.通过仿真分析和实测数据检验，得到如下结论：

(1) 盲源分离法可以实现不同种类信号的分离，当两种信号的周期接近时，此方法依然有效.当使用参考台数据时，可对分离结果进一步约束优化，提升分离效果.

(2) 基于参考台的盲源分离法可有效分离地磁场观测数据中的近场噪音和空间场源信号，可实现近场噪音的抑制/剥离.参考台站数据在确定场源种类和性质时具有重要作用，是评价分离结果的一个重要参考.

上述结论表明基于参考台的盲源分离法是一种较好的时间域去噪方法，在地磁场数据处理、电磁法勘探等领域具备潜在使用价值.应用其提取到的近场噪音可能包含地下电导率信息，值得进一步研究.

致谢  本文所使用参考台数据由日本气象厅提供.