地球主磁场约占地球总磁场强度的95%，主要起源于地球液态外核磁流体发电机过程，地球主磁场控制着地磁场的全球分布，并伴随着千年尺度的长期变化和百万年尺度的极性倒转。研究地球主磁场的空间分布特征，对认知三维地幔电导率结构，理解核幔边界(core-mantle boundary, CMB)动量、能量耦合机制，预测地磁强度变化和极性倒转具有重要意义^[1]。

欧洲太空局(European Space Agency, ESA)地球磁场探测卫星任务SWARM于2013-11-22成功实施，主要研究并加深对地球磁场及其变化规律的认识。SWARM任务包含3颗不同高度的极轨卫星，其中两颗极轨卫星轨道高度460 km，第3颗轨道高度530 km，3颗极轨卫星组成星座对地磁场进行探测。SWARM任务的主要目标是提供地磁场的最佳测量及首个从一小时到数年时间尺度上变化的全球表示，更具挑战性的工作是分离各种磁场来源，SWARM将同时获得地球内部场源和电离层-磁层电流系统的时空特性^[2]。相较于单卫星任务，SWARM数据将显著提高现有磁场模型的精度，并提供高分辨率、高精度的近地磁场的新模型，为调查迄今为止未被发现的地球内部特征提供了新的前景^[3]。SWARM任务以更好的空间分辨率进行长期多卫星同步观测，将进一步提高地核磁场及其长期变化模型的精度。

本文给出了地球主磁场模型反演中SWARM卫星磁测数据预处理流程和基本步骤，分析了主磁场信号分离的各种相关模型的改正方法，并利用精化后的SWARM卫星磁测数据，采用高斯球谐分析方法建立了全球主磁场模型WHU-MM01，联合EMM2015、IGRF12和POMME10模型并对所建模型的精度进行了评估和检验。

1 地球主磁场模型球谐理论

全球地磁场建模的经典方法是高斯球谐分析(spherical harmonic analysis, SHA)，该理论严格依据位场理论建立。球谐分析法使用截断球谐级数的和代表全球地磁场的时空分布，为现今地磁学的蓬勃发展奠定了坚实基础。应用研究表明，球谐分析依然是反演全球地磁场模型的主要方法^[4]。

地磁位场的球谐函数表达式如下：

$V\left( {r, \theta , \lambda } \right) = a\mathop \sum \limits_{n = 1}^\infty \mathop \sum \limits_{m = 0}^n {\left( {\frac{a}{r}} \right)^{n + 1}}\left[ {g_n^m{\rm{cos}}m\lambda + h_n^m{\rm{sin}}m\lambda } \right]P_n^m\left( {{\rm{cos}}\theta } \right)$

(1)

式中，V是主磁场标量位；a为地球半径；r为地心距；θ、λ分别为地理余纬度和地理经度；g_n^m、h_n^m为通常所称的地磁场球谐系数；P_n^m(cosθ)是n阶m次施密特准归一化的缔合勒让德函数。对磁标量所在测点的直角坐标系的3个方向求负梯度，可得到相应的地磁分量^[5]，即：

$ \begin{array}{*{20}{c}} {X = - {B_\theta } = \frac{{\partial V}}{{r\partial \theta }} = }\\ {\mathop \sum \limits_{n = 1}^\infty \mathop \sum \limits_{m = 0}^n {{\left( {\frac{a}{r}} \right)}^{n + 2}}\left[ {g_n^m{\rm{cos}}m\lambda + h_n^m{\rm{sin}}m\lambda } \right]\frac{{\partial P_n^m\left( \theta \right)}}{{\partial \theta }}}\\ {Y = {B_\lambda } = - \frac{{\partial V}}{{r{\rm{sin}}\theta \partial \lambda }} = }\\ {\mathop \sum \limits_{n = 1}^\infty \mathop \sum \limits_{m = 0}^n {{\left( {\frac{a}{r}} \right)}^{n + 2}}\left[ {g_n^m{\rm{cos}}m\lambda - h_n^m{\rm{sin}}m\lambda } \right]\frac{{mP_n^m\left( \theta \right)}}{{{\rm{sin}}\theta }}}\\ {Z = - {B_r} = \frac{{\partial V}}{{\partial r}} = - \mathop \sum \limits_{n = 1}^\infty \mathop \sum \limits_{m = 0}^n \left( {n + 1} \right){\rm{ \times }}}\\ {{{\left( {\frac{a}{r}} \right)}^{n + 2}}\left[ {g_n^m{\rm{cos}}m\lambda + h_n^m{\rm{sin}}m\lambda } \right]P_n^m\left( \theta \right)} \end{array} $

(2)

式(2)即为地磁场三分量的高斯球谐表达式，N为阶次的截断阶数；球谐系数g_n^m、h_n^m的个数分别为N(N+3)/2和N(N+1)/2，球谐系数共有N(N+2)个。若已知地球表面及其外部空间任意一点的地磁场3分量，就可以利用球谐分析方法计算得到地磁场的球谐展开系数。

2 SWARM卫星磁测数据的预处理与精化

本文选取轨道高度在450~500 km的SWARM A星和C星的地磁观测数据反演地磁场模型，时间跨度为2014.0~2016.0。作为卫星磁测数据和模型反演之间的桥梁与纽带，卫星磁测数据的预处理在建模中有着不可或缺的重要地位，磁测数据预处理的方法与结果直接影响到后期反演模型的精度和可靠性^[6]。

在对数据进行处理之前，引入各向异性误差可以减弱小尺度磁场信号的混频效应，进而有效控制最终模型^[7]。首先以10 s为间隔对数据进行二次重采样，对应的空间分辨率约为70 km。

数据预处理包括数据格式、时间和坐标系统的转换，数据筛选和非主磁场信号的模型改正^[8]。

本文选取了地方时为21:00至次日03:00的数据，最大限度削弱地磁静日变化Sq对主磁场信号的干扰。同时为了弱化磁暴、磁层大尺度对流等地磁活动对卫星磁测数据的影响，引入Dst、Kp指数条件，选取对应地磁活动平静期间的数据^[7]：

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

(3)

式中，d|Dst|和d|Kp|分别表示前3 h内的地磁扰动指数的变化率。

为消除粗差对模型反演的影响，本文将各地磁要素实际观测值与现有先验模型(如IGRF12)相应值进行综合对比，利用两者差值剔除可能存在的粗差^{[9, 10]}。具体的判断条件如下：

$ \left\{ \begin{array}{l} {\rm{|}}{X_{{\rm{obv}}}} - {X_m}{\rm{|}} > = 100\;{\rm{nT}}\\ {\rm{|}}{Y_{{\rm{obv}}}} - {Y_m}{\rm{|}} > = 50\;{\rm{nT}}\\ {\rm{|}}{Z_{{\rm{obv}}}} - {Z_m}{\rm{|}} > = 100\;{\rm{nT}}\\ {\rm{|}}{F_{{\rm{obv}}}} - {F_m}{\rm{|}} > = 100\;{\rm{nT}} \end{array} \right. $

(4)

式中，X_obv和X_m分别为地磁要素观测值和IGRF12模型值(Y、Z分量类似)，每个地磁分量粗差的阈值通过已有模型相应分量的误差大小经验来确定。一般情况下，地磁场水平分量X和垂直Z分量受到外源的干扰相对较大，其阈值稍大于Y分量。

除此之外，本文结合观测数据的Quality值，对用于确定卫星姿态的Star-Camera及其平台的运行状态进行判断^[11]。最后，为确保观测数据的质量，本文剔除了连续超过240个历元没有观测数据的短弧轨道。同时，为了避免模型在水平和垂直方向上因数据分布不均而产生震荡效应，本文对观测数据进行了格网化处理，格网的空间分辨率为2°×2°。

为了从卫星磁测数据中提取主磁场信号，需对上述挑选的数据进行非主磁场改正，消除地壳磁场及外源场的影响^[12]。本文选用了EMM2015模型(16~720阶)对地壳磁场信号进行剔除，对于电离层-磁层磁场的改正，则分别选用了ESA发布的DIFI2电离层磁场模型和CIRES发布的POMME10的模型。图 1给出了剔除非主磁场后SWARM卫星观测值的径向磁场分量。为了分析非主磁场信号改正后的残差大小，从改正后的信号中剔除了IGRF12主磁场模型值^[13]。从图 1可以看出，在包含中低纬范围内的大部分区域，SWARM观测值外源场残差分量小于5 nT，在两极地区，残差分量略有增大，但绝大部分残差小于25 nT。

3 反演结果与精度评估

基于高斯球谐分析方法，利用精化后的SWARM卫星磁测数据计算得到了地磁主磁场模型WHU-MM01球谐系数，最高阶数15，这些系数可用于计算2014.0~2016.0地球表面及近地空间任何一点处的地磁分量及其长期变化。

本文基于WHU-MM01模型反演系数，计算了同一观测历元同一点上模型反演值与卫星观测值的差值和均方根误差。比较分析发现，各地磁要素的模型拟合误差绝对值均优于29 nT，X、Y、Z分量的均方根误差分别优于6 nT、7 nT和3 nT。图 2给出了WHU-MM01模型相对于卫星观测值的拟合结果。

从图 2可以看出，在中低纬度区域，模型拟合误差较小，但在南北极区域，模型拟合误差相对较大。在两极地区，非主磁场信号改正后的径向残差分量依然较大，因此，造成拟合误差较大的原因可能与南北极区域地磁扰动较大密切相关。

为了进一步验证模型的可靠性和稳定性，本文基于WHU-MM01模型的反演系数，以2°×2°的格网分辨率，计算了地球表面上格网点的主磁场地磁要素X、Y、Z、F的强度，将其与EMM2015(1~15阶)模型的计算结果进行对比，并给出了其与EMM2015(1~15阶)、IGRF12(1~13阶)和POMME10(1~15阶)模型的互差比较。图 3给出了WHU-MM01模型与EMM2015(1~15阶)模型地磁要素差值ΔX、ΔY、ΔZ、ΔF差值的等值线分布图。

由图 3可以看出，WHU-MM01模型与EMM2015模型吻合较好，在中低纬度区域，地磁要素X、Y、Z、F的互差较小，在两极地区，互差相对较大，但整体上差异均优于3 nT。

表 1给出了WHU-MM01模型与EMM2015、IGRF12和POMME10模型对应的互差全球分布统计结果。WHU-MM01与EMM2015、IGRF12和POMME10等3个模型的比较结果表明，WHU-MM01与EMM2015一致性最好，F分量标准差为2.2 nT，与IGRF12互差较大，而Z分量标准差达到了10.8 nT。综上所述，WHU-MM01模型与EMM2015模型的精度最为接近。

4 结束语

通过计算WHU-MM01模型值和SWARM卫星观测值的差值，从内符合精度对模型进行检验，结果表明模型拟合误差较小，X、Y、Z地磁分量的均方根误差分别优于6 nT、7 nT和3 nT；通过与国际主流的EMM2015(1~15阶)、IGRF12(1~13阶)和POMME10(1~15阶)地球主磁场模型进行对比，发现其与EMM2015、IGRF12和POMME10模型精度相当，量级互差3~10 nT，接近国际同类模型的精度。