地球外源场约占地球总磁场的1%^[1]，它起源于近地空间电流体系，这些电流体系主要分布在电离层-磁层，行星际空间电流变化虽有影响但贡献甚微^[2]。外源场最明显的特征是具有快速的时空变化分布，且空间分布具有全球性^[1]。由于外源场携带着丰富而复杂的地球空间电磁环境信息，因此，研究如何从地球各圈层磁场中提取外源场信号，分析外源场的时空变化特征，对于认知空间电磁环境的变化规律，理解日地能量耦合机制，实时准确预报空间天气变化具有重要的科学意义和潜在的应用价值。

地磁场模型参数化的关键问题之一是对地球各圈层磁场的贡献实施有效分离。对于外源磁场的分离，早期主要根据电离层-磁层磁场的时空特征，借助于各类筛选条件(如地方时条件、纬度条件、Dst/ Kp指数条件)对外源场扰动变化较大的数据进行剔除，然后采用多项式拟合^[3]、球谐分析^[4]或沿轨滤波技术对外源场的贡献进行分离^[5]，如基于沿轨快速傅里叶变换的Kaiser滤波^[6]或Butterworth滤波^[7]，但滤波难点在于截止波长的确定^{[5, 7]}，容易引入信号泄露误差，应用时需分析未模型化的残差信号强度。尽管早期国内外学者基于卫星磁测数据，依据外源场的时空特征对内外场源分离和模型参数化进行了相关理论分析^[8-12]，但由于单颗卫星完成覆盖全球的观测周期通常为几天，难以实现全球范围内电离层-磁层磁场的准同步观测^[13]，因此，很难建立实时高精度的外源场模型。2013-11-22，Swarm卫星成功发射，不仅为地磁场提供了迄今为止的最高精度测量，而且首次实现了对地磁场进行1 h到数天时间尺度上变化磁场的全球观测，Swarm星群的高时空覆盖采样、优化的数据处理流程^{[5, 14]}，有利于分析空间电磁环境的时空变化规律，为外源场的完整分离和模型化提供了强有力的数据和技术支撑。2016年，欧空局Swarm星群应用与研究中心分别发布了最新的电离层磁场MIO_SHA_2D^[15]和磁层磁场MMA_SHA_2C^[10]系列模型，为外源场的改正提供了重要的参考和依据。

内源场和外源场的信号频率存在一定重叠^{[1, 16]}，很难通过滤波方法对其进行准确分离。若对全球或局部岩石圈磁场进行建模，在实施主磁场模型改正之后，剩余磁场主要包括岩石圈磁场、外源场及其感应磁场，主磁场残留的长波信号比较微弱。相对于外源场的全球尺度，起源于磁性岩石圈的岩石圈磁场，在卫星高度上的空间尺度通常为几百公里，甚至更小，因此在建立全球或局部岩石圈磁场模型过程中，可以尝试采用沿轨高通滤波对外源场进行分离。

高通滤波器从实现方法上可分为无限冲激响应(infinite impulse response，IIR)和有限冲激响应(finite impulse response，FIR)^[17-18]。由于外源场时变特征明显，为使不同频率分量的信号经滤波后具有精确统一的时标，本文采用具有精确线性相位特征的FIR滤波器。FIR滤波器的常用设计方法主要包括窗函数法和频率采样法^[13]。窗函数法相对简单，物理意义清晰，通过选择合适的窗函数可以有效抑制信号泄露。

为了分离剩余磁场中外源场的长波信号，本文研究基于FIR滤波器的Hamming窗函数法，通过设计合理的滤波参数，对Swarm卫星磁静日中低磁纬度(±50°)的观测数据进行沿轨高通滤波，并引入最新的MIO_SHA_2D、MMA_SHA_2C电离层-磁层磁场改正模型，对滤波方法的可行性和精度进行评估。

1 数据源与数据预处理 1.1 数据筛选

本文的主要目的是采用沿轨高通滤波对外源场信号进行改正，并通过最新的国际高精度外源场改正模型验证该方法的可行性。为此，本文选取了2015年1月Swarm A星磁静日中低磁纬度(±50°)的低频观测数据(采样频率为1 Hz)，在进行数据筛选之前，首先进行半圆轨道分离和10 s间隔的二次重采样，轨道分离的目的是为了对观测数据进行沿轨滤波。在卫星轨道高度，10 s间隔的二次重采样大约对应于75 km的空间分辨率，可减弱小尺度磁异常信号的混频效应。

考虑到MIO_SHA_2D电离层磁场模型主要描述了中低纬度地区磁静期Sq和赤道电集流磁场及其感应磁场^[15]，本文参考国际上大多数地磁场建模中所用到的地磁指数条件^{[16, 19-20]}，利用式(1)选取地磁扰动相对平静的观测时段数据：

$ \left\{ {\begin{array}{*{20}{l}} {|{\rm{Dst}}| \le 20\;\;{\rm{nT}}\;\;{\rm{ and }}\;\;{\rm{d}}|{\rm{Dst}}| \le 10\;\;{\rm{nT}}}\\ {Kp \le 2\;\;{\rm{ and }}\;\;{\rm{d}}|Kp| \le 2} \end{array}} \right. $

(1)

式中，d| Dst|和d| Kp |分别表示Dst、Kp在前3 h内的变化值。为了充分保证数据的观测质量，通过flags值剔除了Star-camera运行状态欠佳的观测数据，并对X、Y、Z矢量和标量F互差大于10 nT的观测数据予以剔除，最后，剔除了连续超过240个历元没有数据的短弧轨道。图 1给出了基于卫星观测值提取磁静期剩余磁场的预处理流程。

1.2 坐标与时间系统转换

由于MIO_SHA_2D和MMA_SHA_2C是基于地磁坐标系所建立的电离层磁场模型，因此在进行模型计算之前，需先将各测点的位置坐标从球坐标系(λ, θ)转换到地磁坐标系(Λ, Θ)，然后再将模型所计算的地磁分量从地磁坐标系(N, E, U)_dp转换到球坐标系(N, E, U)_geo。此外，MIO_SHA_2D模型还涉及到协调世界时(coordinated universal time, UTC)与地磁世界时t_m之间的转换。

1.2.1 坐标系统转换

1) 球坐标系与地磁坐标系之间转换，可通过球面三角公式实现^{[1, 3]}：

$ \left\{\begin{array}{l}{\cos \mathit{\Theta}=\cos \theta \cos \theta_{0}+\sin \theta \sin \theta_{0} \cos \left(\lambda-\lambda_{0}\right)} \\ {\sin \mathit{\mathit{\Lambda}}=\sin \theta \sin \left(\lambda-\lambda_{0}\right) / \sin \mathit{\Theta}}\end{array}\right. $

(2)

式中，(Λ, Θ)表示点P的地磁经度和地磁余纬；(λ, θ)表示点P的地心经度和地心余纬；(λ₀, θ₀)表示地磁北极N_m的地心经度和地心余纬。

2) J2000惯性坐标系与地固空间直角坐标系之间的转换。

计算地磁世界时t_m需要确定太阳在地磁坐标系中的日下点经度。对于给定的UTC时间，通常先根据JPL提供的星历文件计算太阳在J2000惯性坐标系下的位置和速度，然后按照式(3)将其转换到地心地固直角坐标系下。

$ {(x, y, z)^{\rm{T}}} = (\mathit{\boldsymbol{PR}})(\mathit{\boldsymbol{NR}})(\mathit{\boldsymbol{B}}1)(\mathit{\boldsymbol{B}}2){(X, Y, Z)^{\rm{T}}} $

(3)

式中，(PR)是岁差矩阵；(NR)是章动矩阵；(B1)是极移矩阵；(B2)是地球自转矩阵^[20-22]，其定义可参考文献[23]。

3) 地磁观测分量从地磁坐标系转换到球坐标系下，可通过如下步骤实现：

① 在地磁坐标系参考框架下，利用式(4)将站心地磁计算分量(N, E, U)_dp转换到地磁空间直角坐标系(X, Y, Z)_dp；

$ (X, Y, Z)_{\mathrm{dp}}=\boldsymbol{R}_{1} \cdot(N, E, U)_{\mathrm{dp}} $

(4)

$ \boldsymbol{R}_{1}=\left[\begin{array}{ccc}{-\sin \mathit{\Phi} \cos \mathit{\Lambda}} & {-\sin \mathit{\Lambda}} & {-\cos \mathit{\Phi} \cos \mathit{\Lambda}} \\ {-\sin \mathit{\Phi} \sin \mathit{\Lambda}} & {\cos \mathit{\Lambda}} & {-\cos \mathit{\Phi} \sin \mathit{\Lambda}} \\ {\cos \mathit{\Phi}} & {0} & {-\sin \mathit{\Phi}}\end{array}\right] $

(5)

式中，R₁为旋转矩阵；Φ为点P的地磁纬度。

② 利用式(6)将(X, Y, Z)_dp转换到地心球空间直角坐标系(X, Y, Z)_geo；

$ (X, Y, Z)_{\mathrm{geo}}=\boldsymbol{R}_{2}^{-1} \cdot(X, Y, Z)_{\mathrm{dp}} $

(6)

$ \boldsymbol{R}_{2}=\left[\begin{array}{ccc}{\cos \theta_{0}} & {0} & {-\sin \theta_{0}} \\ {0} & {1} & {0} \\ {\sin \theta_{0}} & {0} & {\cos \theta_{0}}\end{array}\right]\left[\begin{array}{ccc}{\cos \lambda_{0}} & {\sin \lambda_{0}} & {0} \\ {-\sin \lambda_{0}} & {\cos \lambda_{0}} & {0} \\ {0} & {0} & {1}\end{array}\right] $

(7)

式中，R₂为旋转矩阵。

③ 在地心球坐标参考框架下，利用式(8)将(X, Y, Z)_geo转换为站心地磁观测分量(N, E, U)_geo。

$ (N, E, U)_{\mathrm{geo}}=\boldsymbol{R}_{3} \cdot(X, Y, Z)_{\mathrm{geo}} $

(8)

$ \boldsymbol{R}_{3}=\left[\begin{array}{ccc}{-\sin \varphi \cos \lambda} & {-\sin \varphi \sin \lambda} & {\cos \varphi} \\ {-\sin \lambda} & {\cos \lambda} & {0} \\ {-\cos \varphi \cos \lambda} & {-\cos \varphi \sin \lambda} & {-\sin \varphi}\end{array}\right] $

(9)

式中，R₃为旋转矩阵；φ为点P的地心纬度。

1.2.2 时间系统转换与计算

1) 地磁世界时t_m的计算^[15]：

$ t_{m}=\left(180^{\circ}-\varphi_{d, s}\right) / 15 $

(10)

式中，φ_{d, s}为地磁坐标系日下点的地磁经度。

2) 格里高利历与儒略日JD之间的转换：

$ \begin{array}{l} {\rm{JD}} = {\mathop{\rm int}} (365.25y) + {\mathop{\rm int}} (30.6001(m + 1)) + \\ D + h{\rm{ / 24 + 1720981}}{\rm{.5}} \end{array} $

(11)

式中，当M≤2时，则y=Y－1，m=M+12；当M>2时，则y=Y，m=M；h=H+Min/60+S/3 600，其中Y为年，M为月，D为日，H、Min、S分别为时、分、秒；JD为儒略日；int表示向下取整。

2 基于Hamming窗的FIR高通滤波器参数的确定

相对于岩石圈磁场，外源场的空间分布具有全球性，其信号波长较长，频率较低，因此可借助高通滤波器对外源场信号进行分离。高通滤波器从实现方法上可分为无限冲激响应IIR和有限冲激响应FIR。考虑到外源场时变特征明显，不同频率分量的信号经滤波后需要精确统一的时标，本文选用具有精确线性相位特性的FIR滤波器。FIR滤波器的常用设计方法主要包括窗函数法和频率采样法^[24]。窗函数法是从时域进行设计，频率采样法是从频域进行设计。窗函数法相对简单、物理意义清晰，通过选择合适的窗函数可以有效抑制信号泄露，本节将重点研究利用Hamming窗函数法分离外源场信号的FIR高通滤波器参数设置问题。

2.1 归一化的截止频率

Hamming窗函数的基本参数包括通带、阻带截止频率、过渡带宽及滤波器长度。卫星磁测数据中，信号频率f与卫星速度v、波长分辨率λ之间的关系为f=v/λ^[25]。当卫星速度一定时，通带截止频率f_p决定了输出信号的波长分辨率。若已知岩石圈磁场和外源场的临界频率，则可以明确确定滤波器的通带截止频率，事实上岩石圈磁场和外源场信号频率存在一定的重叠，因此需要设计合理的滤波参数，进而有效抑制吉布斯效应。

阻带截止频率f_s取决于过渡带宽Δω，Δω=ω_p－ω_s，对于Hamming窗，过渡带宽与滤波器的窗口长度N之间的关系为：

$ \Delta \omega=\omega_{p}-\omega_{s} \geqslant 6.6 \pi / N $

(12)

式中，ω_p、ω_s分别为滤波器的通带和阻带数字截止频率，

$ \left\{\begin{array}{l}{\omega_{p}=f_{p} \times 2 \pi / f_{s t}} \\ {\omega_{s}=f_{s} \times 2 \pi / f_{s t}}\end{array}\right. $

(13)

相应的归一化的模拟截止频率分别为：

$ \left\{\begin{array}{l}{f_{p}^{{N}}=f_{p} \times 2 / f_{s t}} \\ {f_{s}^{{N}}=f_{s} \times 2 / f_{s t}}\end{array}\right. $

(14)

2.2 滤波器长度

在给定通带、阻带截止频率的情况下，可根据式(12)近似计算滤波器窗口长度N。滤波器窗口长度并不具有唯一性，在确定窗口长度时，应着重考虑两个原则：①滤波器长度N的奇偶性，对于长度相当的滤波器，N为偶数时性能比N为奇数时好，但显著的缺点是滤波时产生的时延(N-1)/2不具有整数特性，需进行内插^[24]；②窗口长度越短，滤波器群时延越小；窗口长度越长，过渡带减小，幅频响应更加精准，但当信号频率与采样频率之比很小时，可能引起类似码间串扰效应，无法得到稳定输出。

3 算例分析

本文分析了磁静日期间2015年1月Q1~Q10(平静日)10 d的Swarm A星观测数据，采用图 1的数据筛选方法选取Q1~Q10观测时段相对平静的半圆轨道，然后采用EMM2015模型(1~720阶)对内源场进行改正，提取剩余磁场。

以2015-01-20某条半圆轨道为例，图 2和图 3分别给出了地磁纬度为±50°范围内沿轨方向上X、Y、Z分量的剩余磁场强度及径向剩余磁场Z分量的幅频响应特性。从图 2可以看出，在中低纬度地区，剩余磁场强度为±10~±20 nT；随着纬度的升高，X、Y分量的剩余磁场表现较为平稳，而Z分量则表现为增强趋势，但其强度无明显突变，说明在中低纬度地区的磁静期，相对于水平方向，外源场对径向分量Z的影响相对较大。从图 3可以看出，剩余磁场主要集中在低频段，表现为长波信号，这些长波信号主要源于大尺度磁层磁场，其次源于夜间电离层磁场及EMM2015模型的残差项。

下面分别采用Hamming窗函数法和MMA_SHA_2C、MIO_SHA_2D模型对剩余磁场的外源场成分进行改正和精度分析。根据式(12)，设定滤波窗口长度为87阶，通带和阻带波长分别为1 200 km和3 000 km。图 4给出了两种方法计算的径向分量(Z)强度，图 5和图 6分别给出了两种方法改正后的径向残差(ΔZ)大小及相应的幅频响应特性，其中黑色曲线代表模型改正结果，红色曲线代表滤波结果。

从图 4可以看出，滤波结果与模型改正结果的一致性吻合较好。图 5表明：①采用两种方法改正后的径向残差均优于±3 nT，对比图 2可以看出，两种方法对大尺度的剩余磁场均能作出有效改正；②两种方法改正后的残差互差约为2 nT，相对于模型改正，滤波方法对大尺度的磁场信号改正更为彻底，但滤去的信号并不能完全界定为外源场，因为这些大尺度的磁场含有未模型化的主磁场信号，甚至可能含有大尺度的磁异常信号；③在沿轨方向上，经滤波和模型改正后的残差分布具有相似的变化特征，表明两种方法对大尺度剩余磁场的改正具有一致性，图 5所描述的径向残差Z分量具有相似的幅频响应特性，验证了该滤波方法的可行性。由图 6可知，滤波后的低频信号强度大大减弱，这表明基于Hamming窗函数法能有效滤去大尺度的剩余磁场信号。

为了进一步分析滤波器长度对剩余磁场长波信号改正的精度影响，在给定通带波长1 200 km、阻带波长3 000 km的情况下，以外源场模型改正后的残差信号强度为参考，采用Hamming窗函数法，对2015年1月Q1~Q10共10 d经筛选后的14条半圆轨道的剩余磁场进行沿轨高通滤波，通过引入均方根误差(root mean square error，RMSE)，对比分析不同窗口长度滤波后的残差与模型改正后的残差，表 1给出了统计分析的结果。

从表 1可以看出，随着滤波器长度的增加，在水平方向上，RMSE呈增大趋势，在径向上，RMSE的最小值在增大，最大值在减小，其均值呈现减小趋势，表明剩余磁场水平分量和径向分量对滤波器长度具有不同的响应特性。因此，为了均衡滤波器在水平方向和径向上的滤波精度，滤波器的长度不宜选得过小或过大。

4 结语

本文利用Swarm A星2015年1月的磁静日观测数据，采用Hamming窗函数法沿轨高通滤波对大尺度的外源场信号进行分离，并借助于最新的MIO_SHA_2D和MMA_SHA_2C电离层-磁层磁场模型验证了该方法的可行性，同时给出了滤波器窗口长度的确定方法，为滤波器的参数设计提供参考。研究表明：

1) 经EMM2015(1~720阶)模型改正后的剩余磁场主要以外源场信号为主，其次源于EMM2015模型的残差项，利用Hamming窗函数法的沿轨高通滤波可对中低纬度地区大尺度的外源场信号实施有效分离，其结果与外源场模型的计算值一致性吻合较好；

2) 采用Hamming窗函数法和外源场模型改正后的径向残差均优于3 nT，相对于模型改正，滤波方法对大尺度的磁场信号分离更为彻底，若以外源场模型为参考，则滤波后的残差约为2 nT，表明该方法对外源场的分离具有较高的精度，但滤波存在信号泄漏问题，因此，该方法只能为外源场信号的分离提供借鉴和参考，其可靠性尚需作进一步的论证；

3) 剩余磁场水平分量和径向分量对滤波器的窗口长度具有不同的响应特性，在给定通带波长1 200 km、阻带波长3 000 km的情况下，将滤波器长度设置为87~129阶，可均衡滤波器在水平方向和径向上的滤波精度。