内容摘要：Geoist的magmod模块，支持了8种常见的地磁场模型解算功能。地磁场无时不刻不在变化，一般可以分为：主磁场、岩石圈场和外源场。其中主磁场最强可以达到50000nT以上，外源场次之约1000nT，岩石圈场最弱。由于地球的地磁北和地理北不重合，地磁北总是在变化之中，历史上还曾经发生过多次的倒转。但总体来说变化不是很快，因此，依据现有观测数据可以用模型还描述地球主磁场的这种变化，通常5年更新一个模型。常说的IGRF和WMM就是最有名的两款，但是如果你要计算某一个时间地点的地磁场参考值，下载模型文件后还必须进行一系列解算才能实现，今天我们就说说Geoist软件包中，对地磁模型解算的函数调用方法。
1、支持的地磁模型
目前，支持IGRF、WMM、CHAOS、MF、EMM、LCS、SIFM、DIFI共8个模型，这些模型都是以高斯球谐分析方法为基础，结合大量地磁场实测资料建立的。下面是这些模型的详细信息：
（1） IGRF模型
国际地磁参考场模型（International Geomagnetic Reference Field，简称IGRF）是国际上通用的标准地磁模型，是一种根据高斯理论建立的用于描述地球主磁场及其长期变化在全球分布的数学模型。它是由国际地磁和高空大气物理协会（IAGA）的V-MOD工作小组建立和维护的。模型每隔5年更新一次。1968年，IAGA发布了第1代IGRF模型（IGRF-1）。2014年，IAGA发布了第12代IGRF模型（IGRF-12），至此，IGRF包括了1900-2015年（间隔5年）共24个主磁场模型，其适用的时间范围为1900.01.01-2019.12.31。IGRF-12融合了卫星（Swarm：2013-2014、Ørsted：1999-2013、CHAMP：2000-2010、SAC-C:2001-2013）数据、地面台站观测数据和地面流动磁测数据，球谐系数13阶，对应的空间分辨率为3000 km。该模型地磁场的全球估计精度约为50-300 nT。
（2） WMM 模型
世界地磁场模型 (World Magnetic Mode l，简称WMM)是一种主要用于描述地球主磁场，同时也兼顾到岩石圈磁场和海洋感应磁场长波成分的数学模型。WMM是IGRF的候选模型之一。该模型主要为美国、英国国防部、北大西洋公约组织( NATO)和国际海道测量组织(WHO)提供导航及定向服务，同时在民用导航定位系统和航向姿态测量系统中也有着广泛应用。WMM是由美国国家地理空间情报局（NGA）和英国国防地理中心（DGC）提供资助，并由美国国家地球物理数据中心（NGDC）联合英国地质调查局（BGS）共同研制的世界地磁模型。模型每隔5年更新一次。最新的WMM2015模型在2014年12月发布，有效使用期为2015.01.01-2019.12.31。WMM2015模型所使用的数据主要包括卫星磁测（Swarm：2013-2014、Ørsted：1999-2013、CHAMP：2000-2010）和地面台站时均值两种类型。该模型的球谐系数是12阶，对应的空间分辨率为3200 km。该模型地磁强度的全球估计精度约为90-170 nT。
（3） CHAOS 模型
CHAOS系列模型( CHAMP/Ørsted/SAC-C model)是描述全球地磁场（包括内源场和外源场）及其长期变化的高精度数学模型。该模型在构建过程中, 采用了一些新的改进技术，如: 重新确定资料筛选标准、矢量资料的坐标转化、外源磁场的拟合等，使模型的可靠性得以提高。它是由丹麦国家空间中心（DTU Space）建立和维护的。2006 年，DTU Space提出了第一代CHAOS模型。第一代CHAOS模型是利用1999年3月至2005年9月的CHAMP、Ørsted和SAC-C三颗卫星的高精度磁测数据计算得到的。DTU Space 在2008-2016年先后提出了xCHAOS、CHAOS－2、CHAOS－3、CHAOS－4、CHAOS－5和CHAOS-6模型。最新的CHAOS-7模型在2019年发布，有效使用期为1999-2020年。CHAOS-7 使用的数据包括卫星磁测数据（来自CHAMP, Ørsted、SAC-C、Cryosat2和Swarm卫星）和182个地面台站数据。CHAOS-7模型球谐级数展开至90阶，其中，地核磁场为1－20阶，岩石圈磁场为21－90阶。
（4） MF 模型
MF模型（Magnetic Field，简称MF）是基于 CHAMP卫星数据构建的描述岩石圈磁场的数学模型。它在推测岩石圈的组成和结构上有独特的优势，其中的长波长成分是编制大陆异常图、全球尺度海洋异常图和航磁异常图的重要依据，而且具有较高的精度。MF模型是由德国国家地球科学研究中心(GFZ)建立和维护的。第1代的MF1模型于2002年发布，用15-80阶球谐函数描述岩石圈磁场的可见部分。最新的MF7模型在2010年8月发布。该模型主要使用了2007年5月至2010年4月的CHAMP卫星数据，球谐级数展开至133阶，对应的空间分辨率为300 km。
（5） EMM模型
增强地磁场模型（Enhanced Magnetic Model，简称EMM）是描述地球主磁场和岩石圈磁场的数学模型。EMM模型是利用球谐分析法描述岩石圈磁场的模型中精度和空间分辨率最高的模型之一，其包含的地球岩石圏磁场信息更加全面精细。该模型被广泛应用于民用定位导航系统中，具有很高的实用价值。EMM模型是由美国国家地球物理数据中心（NGDC）和英国地质调查局（BGS）联合研制和维护的。该模型的前身是NGDC-720模型。最新的EMM2017模型在2017年7月发布，有效使用期为2000-2022年。EMM2017使用的数据包括卫星、 航空、海洋和地面磁测数据，球谐级数展开至790阶（地核磁场为1-15阶，岩石圈磁场为16-790阶），空间分辨率达到51 km。
（6）LCS模型
基于CHAMP和Swarm卫星数据的岩石圈磁场模型（Lithospheric model derived from CHAMP and Swarm satellite data，简称LCS）是描述全球岩石圈磁场的高精度数学模型。LCS模型是由丹麦国家空间中心（DTU Space）的Nils Olsen教授领导的科研小组于2017年首次提出的，目前只有一个版本LCS-1。该模型使用了CHAMP卫星在2006年9月至2010年9月期间以及Swarm 的Alpha 和 Charlie两颗卫星在2014年4月至2016年12月期间的磁测数据。该模型是由磁梯度数据计算得到的。LCS-1用16-185阶球谐函数表示岩石圈磁场，空间分辨率达到220 km。需要注意的是，当球谐阶数大于133时，模型会出现严重的能量泄露现象。
（7） SIFM 模型
Swarm初始磁场模型（Swarm Initial Field Model，简称SIFM）是描述地球磁场（包括主磁场和岩石圈磁场）及其长期变化的数学模型。该模型是由丹麦国家空间中心（DTU Space）的Nils Olsen教授领导的科研小组于2014年提出的。SIFM模型在构建过程中，仅使用了由欧洲航天局 (ESA)发射的Swarm系列卫星（Alpha、 Bravo和Charlie）在2013年11月至2015年1月期间的的磁测数据。该模型首次利用由低轨道卫星数据计算得到的磁梯度信息来提高模型精度，并取得了良好效果（模型估计值与卫星实测值的平均残差值仅为0.12 nT）。SIFM模型球谐级数展开至70阶，对应的空间分辨率为550 km。对静态岩石圈磁场值（15-70阶）的描述，SIFM和MF7模型基本接近。
（8） DIFI 模型
专用电离层磁场反演模型（Dedicated Ionospheric Field Inversion Model，简称 DIFI）是描述地球中低纬度地区（+/- 55度之间）在地磁静日期间的太阳宁静区（Solar-quiet，Sq）和赤道电急流(Equatorial electrojet，EEJ) 磁场（此类磁场属于地磁场中占比例极小的感应磁场）及其变化的数学模型。DIFI模型是由欧洲航天局 (ESA)提供资助，并由科罗拉多大学波尔得分校的环境科学合作研究所（CIRES）和巴黎地球物理学院（IPGP）共同研制的。该模型的计算方法于2013年首次提出。在2015-2016年期间，CIRES发布了DIFI-2015a、 DIFI-2015b和DIFI-2版本。最新的DIFI-3模型在2017年发布。该模型在构建过程中，使用了Swarm系列卫星和地面观测台站在2013年12月至2017年6月期间的磁测数据。DIFI-3模型球谐级数展开至60阶。
2、magmod的函数接口
首先，GEOIST的geoist\magmod\data目录下，已经存放了全部支持模型的球谐函数文件，打开init.py文件，可以查阅到如下信息：
WMM_2015 = join(_DIRNAME, 'WMM2015v2.COF')
WMM_2020 = join(_DIRNAME, 'WMM2020.COF')
EMM_2010_STATIC = join(_DIRNAME, 'EMM-720_V3p0_static.cof')
EMM_2010_SECVAR = join(_DIRNAME, 'EMM-720_V3p0_secvar.cof')
CHAOS6_CORE_X8 = join(_DIRNAME, 'CHAOS-6-x8_core.shc')
CHAOS6_CORE_LATEST = CHAOS6_CORE_X8
CHAOS6_STATIC = join(_DIRNAME, 'CHAOS-6_static.shc')
IGRF11 = join(_DIRNAME, 'igrf11coeffs.txt')
IGRF12 = join(_DIRNAME, 'IGRF12.shc')
IGRF13 = join(_DIRNAME, 'igrf13.COF')
SIFM = join(_DIRNAME, 'SIFM.shc')
LCS1 = join(_DIRNAME, 'LCS-1.shc')
MF7 = join(_DIRNAME, 'MF7.shc')
要想使用这些模型，在python程序中，直接import即可，如下：
from geoist.magmod.data import (
EMM_2010_STATIC, EMM_2010_SECVAR, WMM_2015, WMM_2020, IGRF13,
CHAOS6_CORE_LATEST, CHAOS6_STATIC,
IGRF11, IGRF12, SIFM,)
其次，模型有了就是要加载解算函数了，在GEOIST的geoist\magmod\magnetic_model目录下，有多种类型的解算函数接口，import方法如下：
from geoist.magmod.magnetic_model.loader_igrf import load_model_igrf
from geoist.magmod.magnetic_model.loader_wmm import load_model_wmm
from geoist.magmod.magnetic_model.loader_emm import load_model_emm
接下来，要注意的就是时间定义，表示时间除了公历日期外，还有常用到儒略日，有时候需要相互转换。直接import函数即可，方法如下：
import datetime as dt
from geoist.magmod.time_util import (
decimal_year_to_mjd2000, decimal_year_to_mjd2000_simple,mjd2000_to_decimal_year,mjd2000_to_year_fraction
)
准备好上面的步骤后，我们就可以计算了，不论什么模型，最后计算都要调用eval函数，示例如下：
wmm2020 = load_model_wmm(WMM_2020)
d2 = 2022.5
loc1 = (80.0, 0.0, 100.0)
wmm2020.eval(decimal_year_to_mjd2000(d2), loc1, 0, 0, **options)
简单来讲，就是先import所有要用的接口，然后load模型，再eval一个时间地点的模型，来计算磁场参考值即可。
3、一个完整的简单例子
通过上面的介绍，我们来看看如果要完成的计算一个模型，代码大概长什么样子。下面我们测试一个电离层模型（DIFI4）计算一个地面点在1天24小时的地磁场变化结果。
完整代码如下：
import datetime as dt
from geoist.magmod.data import DIFI4
from geoist.magmod.magnetic_model.loader_mio import (
load_model_swarm_mio_internal,
load_model_swarm_mio_external)
from geoist.magmod.time_util import decimal_year_to_mjd2000
from geoist.magmod.util import datetime_to_decimal_year
loc = (45.0, 105.0, 1.0)
options = {"f107": 70, "scale": [1, 1, -1]} #-1 is Z direction
# load DIFI4 model
difi4 = load_model_swarm_mio_external(DIFI4)
difi42 = load_model_swarm_mio_internal(DIFI4)
#get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import numpy as np
magdifi = np.zeros(3*24).reshape(24,3)
magdifi2 = np.zeros(3*24).reshape(24,3)
for i in range(24):
t1 = decimal_year_to_mjd2000(datetime_to_decimal_year(dt.datetime(2019,1,1,i,0,30)))
magdifi[i] = difi4.eval(t1, loc, 0, 0, **options)
magdifi2[i] = difi42.eval(t1, loc, 0, 0, **options)
plt.title("DIFI-4 IONOSPHERE MAGNETIC FIELD MODEL")
plt.xlabel("UTC time/h")
plt.ylabel("Sq intensity/nT")
plt.plot(magdifi[:,0],'bo', label = "North-X")
plt.plot(magdifi[:,1],'ro', label = "East-Y")
plt.plot(magdifi[:,2],'go', label = "Vertical-Z")
plt.plot(magdifi[:,0] + magdifi2[:,0],'b', label = "North-X2")
plt.plot(magdifi[:,1] + magdifi2[:,1],'r', label = "East-Y2")
plt.plot(magdifi[:,2] + magdifi2[:,2],'g', label = "Vertical-Z2")
plt.legend()
plt.show()
如果正确运行，结果可视化如图1所示：
加载模型，确定接口，计算再加上可视化，计算地磁模型，是不是很简单？
结语：
计算地磁模型，获得指定时空地点的地磁参考场，在地磁勘探、导航领域都有重要的应用。特别是对于高精度磁场参数的预报更加关键，基于magmod程序包提供程序接口可以作为开发多种地磁应用软件的基础，甚至可以进一步开发为在线的API服务接口。
下一节将继续介绍magmod针对IGRF13和WMM2000最新模型的更新支持情况。
