地球物理学报  2018, Vol. 61 Issue (4): 1352-1365   PDF    
地磁场三维曲面Spline模型
冯彦1,2, 蒋勇1, 孙涵3, 安振昌4, 黄娅2     
1. 南京信息工程大学空间天气研究所, 南京 210044;
2. 中国科学院空间天气学国家重点实验室, 北京 100190;
3. 内蒙古新天元防灾减灾研究所, 呼和浩特 010051;
4. 中国科学院地质与地球物理研究所, 北京 100029
摘要:利用149个地面实测数据以及12个子午工程测点数据,50个CHAMP卫星高度实测数据,并结合高空180 km处的50个IGRF12数据点,基于这三个高度的数据首次建立了中国地区地磁场三维曲面Spline(3D Spline)模型.在境外添加了39个测点以控制边界效应.通过CM4模型将所有测点的场源进行分离,统一通过主磁场值建模分析.通过将模拟结果与实测值、曲面Spline(2D Spline)以及Taylor(2D Taylor)、三维Taylor(3D Taylor)模型及IGRF12模型相比较,结果显示3D Spline模型的空间分布与其他模型整体趋势一致,但更为曲折,随着高度上升,3D Spline模型的要素Y的强度逐渐减弱.通过比较3D Spline、2D和3D Taylor模型对于不同高度6个缺测点的模拟值,残差和均方根偏差(RMSE),3D Spline模型的模拟效果最好,要素YZ和总强度F的RMSE值要比其他模型低50%以上.3D Spline模型在不同高度处的模拟效果主要取决于该高度附近的实测值数量和精度.
关键词: 地磁场      三维模型      曲面Spline      子午工程      CHAMP卫星      IGRF12     
The three-dimensional surface Spline model of geomagnetic field
FENG Yan1,2, JIANG Yong1, SUN Han3, AN ZhenChang4, HUANG Ya2     
1. Institute of Space Weather, Nanjing University of Information Science and Technology, Nanjing 210044, China;
2. State Key Laboratory of Space Weather, Chinese Academy of Sciences, Beijing 100190, China;
3. Xintianyuan Institute of Disaster Prevention and Reduction, Inner Mongolia, Hohhot 010051, China;
4. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
Abstract: Based on 149 surface observations and 50 CHAMP measuring data, along with newest 12 Chinese Meridian Project observations, as well as the 50 points of IGRF12 that lie in the level of 180 km, we created the three-dimensional Surface Spline (3D Spline) model of Chinese mainland for the first time. The boundary effect was controlled by 39 complimentary points. All data were removed the different field sources by CM4 model, so the created 3D Spline model expresses the regional main field. We have analyzed and compared the modeling values with that of Surface Spline (2D Spline), Taylor polynomial (2D Taylor), three-dimensional Taylor polynomial (3D Taylor) models and IGRF12. Results show the overall trends of elements X, Y, Z and F among these models are basically consistent, but the isolines of 3D Spline are relatively more tortuous. Moreover, the intensity of Y decreases while altitude increases. Through comparing the modeling values, residuals and Root-Mean-Square-Errors (RMSE) of 6 intentional missing points between 3D Spline, 2D and 3D Taylor models, results imply the former is better than others. RMSEs of Y, Z and F of 3D Spline are over 50% less than that of other models. The modeling precision of 3D Spline model at different altitudes largely depends upon the number and precision of measuring points at adjacent altitudes.
Key words: Geomagnetic field    Three-dimensional model    Surface Spline    Chinese Meridian project    CHAMP satellite    IGRF12    
0 引言

地磁场同地球电场、地球重力场并称为地球的三大地球物理场.地磁场的时空分布分析对于研究地球演化,地球内部的物理化学过程非常重要(徐文耀, 2009).

对于全球尺度而言,一般通过连续测量的地面台站和复测点数据,结合航空磁测与卫星磁测数据而联合建立全球模型以研究.其中最具代表性的模型为2014年公布的第12代国际地磁参考场模型(IGRF12)(Thébault et al., 2015).另外还有最新的第五代地磁场综合模型——CM5(Sabaka et al., 2015),它结合基于最新的SWARM卫星的level-2数据的综合反演算法而建立,不仅可以分离出地磁场的内、外源磁场,还可估算微弱的洋流磁场.最新的第五代CHAOS模型——CHAOS-5(Finlay et al., 2015)是基于连续10个月测量的SWARM星座的三颗卫星以及CHAMP、Oersted及SAC-C卫星数据建立,它除了可独立对以上六颗卫星进行数值模拟外,还可模拟过去十年间地核表面处的三次脉动事件.Olsen等(2015)完全基于SWARM卫星数据建立了SWARM基本地磁场模型(SIFM),模型1~60阶的模拟情况与基于CHAMP卫星数据所建的其他全球模型高度一致.这些新模型均基于卫星及地面数据建立,都能较好地反映全球尺度的主磁场及岩石圈磁场分布,甚至可模拟地磁外源场及其感应场的时空变化.

相对于大尺度的全球模型,一般采用区域地磁场模型以研究小尺度的区域地磁场.常用的区域地磁模型有Taylor多项式(2D Taylor)模型(Le Mouel,1969)、球冠谐模型(SCH)(Haines,1985)、曲面Spline模型(2D Spline)(安振昌和徐元芳,1981)、矩谐模型(RH)(Alldredge,1981)、Legendre多项式模型(LP)等.

2D Spline模型首先由安振昌等(1982)引入地磁学领域,该模型具有计算简单、使用方便、过点插值等优点(冯彦等,2010a),较适于模拟磁场分布较为复杂的区域(冯彦等,2010b).基于2D Spline模型,陈斌等(2016)建立了2010.0年代的“中国地磁参考场曲面样条模型”,并分析了地磁各要素的分布情况.Feng等(2016)基于三维Taylor多项式(3D Taylor)模型和2D Spline模型研究了1960.0、1970.0、1980.0、1990.0和2000.0年代中国地区的地磁异常场,结果显示曲面Spline模型相较其他模型可反映更多的磁异常信息.冯彦等(2013)还基于墨西哥境内的49个地面测点结合2D Spline和CM4模型,研究了该区域的磁异常.高金田等(2006)利用2D Spline模型研究了1900—1936年中国部分地区的地磁场及其长期变化.尽管以上研究都取得了理想的结果,但建模时都假设地面测点位于同一平面,并未考虑高度因素,而地磁场的内源场源于地球液态外核(徐文耀,2009),并随距离的三次方衰减.通过实际测量发现,总强度F的衰减率约为20 nT/km.而中国大陆地区的地表海拔高度范围为-0.15~8.8 km,因此普通的2D建模会产生约-10~180 nT的精度损失.

考虑到2D Spline模型的优点,本研究拟基于该模型首次建立精度更高(可模拟空间磁场分布)的三维曲面Spline(3D Spline)模型,建模数据包括地面台站及CHAMP卫星的实测矢量数据,并采用了集我国最先进观测技术的子午工程磁测数据,并考虑了地磁异常场、地磁外源场(电离层、磁层磁场及其感应场)的影响.另外为了提高模型精度,在地面和卫星高度的中间高度(180 km)位置添加了补充点,从而联合三种高度的数据进行三维建模.柳士俊等(2011)曾首次将2D Taylor模型进行了三维化,该模型也被用于分析卫星高度的地磁场(蒋勇等,2015),并认为其利用较低的截断阶数即可反映更多的地磁场信息.现将3D Spline建模结果分别与常用的2D Spline模型、2D和3D Taylor模型等相比较.文章第1节阐述建模所涉及的数据与方法,第2节展示建模结果,最后给出相应的结论和讨论.

1 数据与方法 1.1 数据

拟建3D Spline模型的数据包括地面台站数据、子午工程数据、CHAMP卫星实测矢量数据以及IGRF12的模型值.

1.1.1 地面数据

1998—2001年进行了全国范围的地面台站及复测点地磁测量.共进行了39个台站数据的采集,另外117个复测点的总强度F、磁偏角D以及磁倾角I数据都通过质子旋进磁力仪CZM-2和磁通门磁力仪DIM-100获得,其精度分别达到了1.5 nT和±0.1′,而FDI的误差分别为±1.5 nT、±0.3′和±0.3′.复测点数据的日变化以及其他干扰已通过距离最近的地磁台站进行了消除,且测量数据通过参考台站数据的时均值而统一归算至统一时间点.其中7个台站为子午工程台站(暂时不考虑在内),最终适用的2000.0年地面测点为149个.

考虑到缺乏地面数据的高度值,现假设所有测点高度不随时间变化,通过全球陆地1 km基准高程数值(http://www.ngdc.noaa.gov/mgg/topo/)获取了以1弧分间隔的高度数据值.通过选取与实测点最接近的网格值而近似获取所有实测数据的高程数据,其误差不超过0.0083°.

拟采用最新的子午工程数据结合现有地面数据和卫星数据联合建模.该工程沿着东经120°和北纬30°共有12个地磁观测台站进行数据采集.对于地磁场测量,共分为两部分.第一部分为绝对观测,所使用的仪器包括质子旋进磁力仪、Overhauser磁力仪以及D/I的磁通门磁力仪.前两种仪器提供了总强度的绝对测量,后一种仪器提供了其他要素的基线值.第二部分为相对观测,该测量采用了磁通门磁力仪以及感应式磁变仪.相对观测为记录磁场变化的测量(Wang, 2010).这12个子午台站于2010.0年的实测矢量数据也参与三维建模.本研究在地面共采用了44个台站(2000.0年的台站与子午台站之和)数据以及117个复测点数据,共161个地面数据进行研究.

1.1.2 CHAMP卫星数据

本研究卫星数据来源于德国地球科学研究中心,具体包括Overhauser磁力计(OVM)测量的level-1标量数据以及磁通门磁力计(FGM)测量的level-1矢量数据.由于矢量数据在卫星测量过程中时刻发生着位置和强度的变化,因此需要首先对level-1矢量数据通过时间校正、数据转换、温度校正等操作才能转换为FGM level-2数据(Flechtner et al., 2010),而该阶段的数据还不能直接建模,还需将FGM参考系的数据转换为NEC(北向、东向、垂直)参考系数据,具体须通过(1)式实现:

(1)

式中R(FGM→ASC)为将FGM参考系数据转换至卫星光具座(ASC)参考系数据;R(ASC→ICRF)为将基于level-2的卫星姿态数据由ASC参考系数据转换至国际天球参考系(ICRF)数据;R(ICRF→ITRF)为将ICRF参考系数据转换至国际地球参考系(ITRF)数据;R(ITRF→NEC)为将ITRF参考系数据转换至NEC参考系数据.另外还需将卫星数据的全球时换算成地方时.参考现广泛使用的CHAOS-3等模型的选取标准,对卫星数据按以下方法进一步筛选:

(1) 选择18:00—23:00LT,00:00—05:00LT时段数据,即夜侧数据;

(2) 选择Kp≤20时间段的数据;

(3) 选择|Dst|≤20 nT时间段的数据.

本研究范围为18°N—54°N,73°E—136°E.实验表明对于某一时间点前后60天内的卫星数据经归算后可获取较为理想的测点分布.由于测点的分布并不均匀,另外在局部地区某一时间点每个卫星测点的权重应一致,因此对归算后的所有数据进行网格化处理,即对每个格点的值取其±1°内所有值的平均值.在境内地面实测点相对稀疏的区域均匀选取CHAMP卫星数据,从而尽量使所有建模数据均匀分布.另外为了减少相邻轨道的差异,只选择卫星上升轨道的数据以建模.现基于2005.0年的卫星网格数据点,考虑到地面实测数据的重要性,故只选取50个卫星网格数据,卫星数据与地面数据之比约为1:3.

由于建模所用数据的不一致,为了便于分析,拟通过移除场源,对所有测点的主磁场进行建模分析.另外为了提高三维模型的精度,考虑到地面与卫星数据存在约360 km的数据真空,我们利用IGRF12模型在180 km处(地面与卫星高度的中间距离)添加50个均匀分布的测点进行建模.假设岩石圈磁场不随时间变化,且主磁场为线性变化,利用IGRF12统一将2010.0年的12个地面子午工程数据、2005.0年的50个CHAMP卫星数据归算至2000.0年.为了约束由数值拟合而产生的边界异常,在境外基于均匀选取法(UA)(Feng et al., 2015)在卫星高度选取了39个IGRF12点作为补充,结合2000.0年地面149个数据点和180 km处50个IGRF12点,实际共有300个测点进行建模.所有测点的分布见图 1.

图 1 中国及邻近地区的所有测点的二维(a)和三维(b)分布. Albers Equal-Area Conic投影 ●地面台站;★子午工程测点;× 180 km处IGRF12补充点;+ CHAMP卫星测点;▲卫星高度的IGRF12补充点. Fig. 1 The 2D (a) and 3D (b) distribution of all measured points over Chinese and its adjacent areas. Albers Equal-Area Conic projection ● Surface points; ★ Meridian points; × IGRF12 points at 180km; + CHAMP points; ▲ IGRF12 points at satellite level.
1.2 建模方法

本研究的重点为基于2D Spline模型而建立的3D Spline模型.

1.2.1 2D和3D Spline模型

2D Spline模型(安振昌等,1982)的表达式为:

(2)

(3)

其中,W为任意地磁分量,x为地理纬度,y为地理经度;ri2=(xix)2+(yiy)2N为测点数;a0a1a2Fi为需求系数,ε为控制表面曲率变化的小量,取值为1×10-7,共有N+3个系数.

在2D Spline模型的基础上,我们效仿3D Taylor多项式模型的方法,首先对(2)式添加高度项a3z,从而加入高度项的系数a3,而由于增加了一项,为了获得基于恰定方程组的确定解,对(3)式增加了,从而构造了含有N+4个系数的N+4个方程.具体的表达式如下:

(4)

(5)

(4) 式的各参数与2D Spline模型一致,z为所增加的高度数据,而ri2现在添加了所有测点与某一高度值差的平方项,即ri2=(xix)2+(yiy)2+(ziz)2a0a1a2a3Fi为需求系数.对于(4)、(5)式分别代入已知测点的纬度、经度、高度及强度数值,通过高斯消元法即可求得对应不同地磁要素的3D Spline模型系数,基于系数回代模型,即可模拟出研究区域内不同空间位置的模型值.

1.2.2 2D和3D Taylor模型

2D Taylor多项式(Le Mouel,1969)的表达式如下:

(6)

式中W表示任意磁场要素;N为截断阶数;φ, λ分别为各地磁测点的纬度和经度;φ0, λ0为多项式的展开原点,Aij为相应各分量的模型系数,每个模型有N(N+1)/2个系数.

通过添加高度项,并通过系数的完整展开而得到3D Taylor模型(柳士俊等,2011),其表达式如下:

(7)

h为各地磁测点的高度;φ0, λ0, h0为多项式的展开原点,φ0=36°N,λ0=104.5°E,h0=180 km,Aijk为相应各分量的模型系数,每个模型有N3个系数,所有系数通过最小二乘法求取.

1.2.3 CM4模型

本研究基于所有测点的主磁场进行建模.地面和子午工程测点已经过了外源场校正,因此需要对于地面和子午工程测点移除岩石圈磁场值和外源感应场,对于CHAMP卫星数据移除岩石圈磁场、电离层、磁层磁场及其外源感应场.现通过CM4模型移除上述场源,它们的表达式可参阅CM4的相关文章(Sabaka et al., 2002, 2004).

1.2.4 IGRF12模型

为了更好地建立三维模型,在地面和卫星测点之间,即180 km高度处,基于最新的IGRF12模型在境内添加了50个均匀分布的主磁场测点值.IGRF12模型于2014年发布(Thébault et al., 2015),它基于Gauss提出的球谐函数而建立.包含了1900—2015年的24个主磁场模型以及2015—2020年的长期变化模型,IGRF12的主磁场截断阶数为13,长期变化的截断阶数为8.

2 结果

通过场源移除及时间归算的所有测点数据,结合(2)—(7)式建立了中国及邻近地区的地磁场2D、3D Spline模型及2D、3D Taylor模型.现对这四种区域模型及IGRF12进行各地磁要素的比较分析.

2.1 地磁场形态分布比较

在建模过程中,由于2D Spline模型未考虑高度因素,因此其模拟结果产生了较大的偏差,具体在2.2节中阐述.而通过比较均方根偏差(RMSE)以及各分量分布形态,认为9阶的2D Taylor模型与3阶的3D Taylor模型可较好地模拟地磁场.

为了验证所建3D Spline模型,还将其直接与基于实测点和Kriging方法绘制的图形比较.Kriging方法是求最优、线性、无偏的空间内插方法,能将预测误差最小化而最优地表示空间分布.现在固定高度——1.11 km(所有地面实测点的平均高度)处将所建模型与直接基于实测点以及Kriging方法的插值分布图比较,另外IGRF12模型也参与比较.为了节省篇幅,选取北向要素X、东向要素Y、垂直要素Z和总强度F进行比较.具体见图 2-5.

图 2 1.11 km处的要素X分布图(单位:nT) Fig. 2 The distribution of X at 1.11 km (unit: nT) (a) 3D Spline; (b) Kriging; (c) IGRF12; (d) 2D Taylor; (e) 3D Taylor.
图 3 1.11 km处的要素Y分布图(单位:nT) Fig. 3 The distribution of Y at 1.11 km (unit: nT) (a) 3D Spline; (b) Kriging; (c) IGRF12; (d) 2D Taylor; (e) 3D Taylor.
图 4 1.11 km处的要素Z分布图(单位:nT) Fig. 4 The distribution of Z at 1.11 km (unit: nT) (a) 3D Spline; (b) Kriging; (c) IGRF12; (d) 2D Taylor; (e) 3D Taylor.
图 5 1.11 km处的总强度F分布图(单位:nT) Fig. 5 The distribution of F at 1.11 km (unit: nT) (a) 3D Spline; (b) Kriging; (c) IGRF12; (d) 2D Taylor; (e) 3D Taylor.

在1.11 km处,要素X都为正值分布,强度随纬度的增加而降低.3D Spline模型的分布与IGRF12和3D Taylor较为接近,但是等值线略微曲折.Kriging方法假设所有测点都在同一高度,将所有测点进行空间内插,由于境内测点较为稠密,以此绘制的分布图会产生较多的磁异常区.9阶的2D Taylor模型同样假设所有测点位于统一平面,由于未考虑高度参数而会产生系数的估算偏差,进而导致了在模型正演时会产生部分虚假信息.根据图 2,Kriging法和2D Taylor模型可基本反映要素X的分布趋势.

对于要素Y,其分布基本为π形,Kriging法显示出了较为曲折的变化,尤其在等值线弯曲处,而3D Spline模型也同样有此特点,IGRF12与Taylor模型的等值线则较为平滑,主要由于它们都为数值拟合模型.相对而言,3D Taylor模型的分布则略有差异,且在西北和东部的正负区域的极值略大.

要素ZF的分布较为接近,而比较结果与要素X类似.总体而言,5种模型的基本分布趋势一致,而3D Spline模型的等值线分布较为曲折,能反映更多的地磁信息.

为了考察3D Spline模型的空间模拟情况,现基于该模型绘制了要素Y在地面(1.11 km)、高空(180 km)以及卫星高度(367.99 km,所有卫星测点的平均高度)处的网格图,并与相应高度的IGRF12相比较.另外基于三种高度的实测数据,结合相应高度的39个IGRF12补充点,利用2D Spline模型绘制的分布图也加入比较.具体如图 6-8.

图 6 1.11 km处的要素Y分布图(单位:nT) Fig. 6 The distribution of Y at 1.11 km (unit: nT) (a) 3D Spline; (b) IGRF12; (c) 2D Spline.
图 7 180 km处的要素Y分布图(单位:nT) Fig. 7 The distribution of Y at 180 km (unit: nT) (a) 3D Spline; (b) IGRF12; (c) 2D Spline.
图 8 367.99 km处的要素Y分布图(单位:nT) Fig. 8 The distribution of Y at 367.99 km (unit: nT) (a) 3D Spline; (b) IGRF12; (c) 2D Spline.

根据图 6-8,在地面1.11 km处,2D和3D模型的等值线都较为弯曲,从分布而言,通过观察0等值线的位置,2D模型和IGRF12较为相似,而在180 km处,2D和3D模型的等值线开始逐渐光滑,但在形态分布上,3D模型与另两种模型仍存在区别,而在卫星高度,无论是形态还是强度,三种模型均高度相似.上述结果表明,3D模型在不同高度处基本可以表示出要素Y的分布,其与2D及IGRF12的差异主要源于建模数据不同,3D模型是基于所有300个实测点而建成,比仅基于200个地面测点而建的2D模型以及仅十余个国内测点的IGRF12更为稠密,因此等值线也更为弯曲,理论而言,其形态也应更为准确.而在更高高度,如180 km和卫星高度处,由于测点数量都为89个,相对于地面大幅减少,导致了等值线的光滑以及分布形态与IGRF12更为接近的趋势.另外比较三种高度的分布,发现要素Y的强度随高度升高而逐步减弱,而形态随高度升高则变化不大.

另外,为了考察3D Spline模型在垂直方向的模拟情况,我们绘制了要素Y分量在8个高度处的模拟分布图.

根据图 9,在除了地面(1.11 km)、高空(180 km)和卫星高度(367.99 km)的其他高度处,模型模拟结果仍存在问题,例如在50 km处,Y都为负值,在100 km处进一步衰减,而在180 km处则又恢复了π形分布,然后随着高度上升而分布逐步简化,当在360 km处,其分布与367.99 km处非常接近.上述结果反映了在缺少测点的高度处,3D Spline模型的模拟结果仍不准确,而在靠近建模测点的高度处,则可得到较为接近的分布,该特性也反映了不同高度实测数据对于3D Spline模型建模的重要性.

图 9 基于3D Spline模型的不同高度处的Y分量分布(单位:nT) Fig. 9 Distributions of Y based on 3D Spline model at different altitudes (unit: nT) (a) 1.11 km; (b) 50 km; (c) 100 km; (d) 180 km; (e) 200 km; (f) 300 km; (g) 360 km; (h) 367.99 km.
2.2 缺测点的模拟比较

为了进一步验证所建立的3D Spline模型,还需要验证模型在已知测点外其他测点处的拟合程度.由于Spline模型是基于过点插值建模,因此模型值会经过所有的实测点,即满足插值条件:F(xi)=P(xi)(i=1, 2, …, N),其中F(xi)为测量值,P(xi)为插值多项式的节点值,N为所有测点数.

为了验证模型的空间模拟的准确性,我们在不同位置、不同高度处任意提取6个实测点不参与建模,再基于剩下294个测点建立了3D Spline模型,基于该模型计算出6个缺测点的模型值,另外同样通过2D Spline、9阶2D和3阶3D Taylor模型估算这6个测点的模型值,并进行比较.具体结果见表 1-4.

表 1 6个缺测点的要素X的模拟比较(单位:nT) Table 1 Modeling and comparison of X of six missing points (unit: nT)
表 2 6个缺测点的要素Y的模拟比较(单位:nT) Table 2 Modeling and comparison of Y of six missed points (unit: nT)
表 3 6个缺测点的要素Z的模拟比较(单位:nT) Table 3 Modeling and comparison of Z of six missed points (unit: nT)
表 4 6个缺测点的总强度F的模拟比较(单位:nT) Table 4 Modeling and comparison of F of six missed points (unit: nT)

根据表 1-4,2D Spline模型的模拟绝对值高于其他模型模拟值约2个数量级.这主要由于2DSpline模型建模时,假设所有不同高度的测点均位于同一平面,因此在计算系数时会产生误差.从数值模拟角度而言,2D Spline模型的系数矩阵的条件数为1.15×1022,远大于1,故系数矩阵为严重病态矩阵,为奇异矩阵,其所得到的系数解有较大偏差.从空间拟合的角度而言,由于相邻测点的高度相差太大,很难进行最佳拟合,故系数的估算会存在误差,从而影响到回代计算的网格值.通过观察,相比实测值,3D Spline和3D Taylor模型的偏差较小,而2D Taylor模型则次之.

为了更好地观察相对于实测值的模拟值,我们绘制了6个缺测点的实测值以及3D Spline、9阶2D和3阶3D Taylor模型的模拟值图(图 10).

图 10 6个缺测点的实测值与拟合值 Fig. 10 Measuring points and their modeling values of six missed points

根据图 10,6个缺测点处的实测值与3D Spline、3D Taylor模型几乎看不出偏差,而与2D Taylor模型还是存在误差,尤其是要素ZF.

由于数值基数太大,图 10无法观察模型相对于实测值的偏差,为了进一步比较,我们绘制了3D Spline、3阶3D与9阶2D Taylor模型相对实测值的残差图图 11.

图 11 6个缺测点的残差值 Fig. 11 The residual values of six missed points

根据图 11,基于2D Taylor模型的残差都要明显大于其他两种模型.相比3D Spline和3D Taylor模型,在点(41.84°N,123.42°E)和点(33°N,95°E)处前者要低于后者,而其他测点则不明显.

除了比较残差,RMSE也是考察模型模拟精度的一个重要指标,我们计算了3D Spline、2D和3D Taylor模型对于6个缺测点的RMSE值.

根据表 5,9阶2D Taylor模型的RMSE值要大于其他两模型1个数量级以上,而除了X的数值比3D Taylor模型略大,3D Spline模型对于要素YZF的RMSE值均小于3阶3D Taylor模型50%以上.这说明了3D Spline模型对于不同高度测点的模拟要好于其他3种区域模型.

表 5 3D Spline、2D和3D Taylor模型的RMSE值(单位:nT) Table 5 RMSEs of 3D Spline, 2D and 3D Taylor models (unit: nT)
2.3 系数分析

表 6列出了3D Spline模型与分别基于地面、高空和卫星高度实测值的2D Spline模型的部分(前5行)系数,从系数角度分析其差异.

表 6 3D Spline与2D Spline模型的系数 Table 6 Coefficients between 3D Spline and 2D Spline models

观察表 6,3D模型与不同高度的2D模型系数还是存在较为明显差异,这点从图 6-8的差异即可发现.模型的首项系数要远大于之后的系数值,从模型系数权重角度而言,首项系数基本决定了模拟的分布趋势和强度.对比发现在地面处的2D模型各要素的首项系数要更为接近3D模型,这是因为地面测点为200个,要远多于高空及卫星高度测点,这也解释了在地面部分,两模型的等值线都更为曲折,且强度接近.而在高空及卫星高度,由于测点数相等(89个)且相对较少,因此这两个高度处所计算的首相系数较为接近,所绘制的图形也较为光滑.图 7图 8显示3D模型仍可以较好地模拟高空和卫星高度的地磁分布,这也反映了3D模型的系数可较好地控制空间拟合.

3 结论与讨论

通过结合较新的子午工程数据,以及高精度的地面和CHAMP卫星实测值,并添加了180 km处的补充点,首次建立了中国地区的主磁场3D曲面Spline模型,通过地磁要素的形态和缺测点的比较分析,我们认为所建模型可较好地模拟不同高度的地磁场分布.基于已有结果,得到如下结论:

(1) 由于考虑了高度因素及建模特点,3D Spline模型有以下优点:因过点插值而无需讨论截断阶数的选取,不存在龙格效应(Feng et al., 2015),根据残差和RMSE的结果,3D Spline模型优于2D Spline、2D和3D Taylor模型.

(2) 3D Spline模型相较Taylor等其他基于数值拟合的区域模型,可反映更多的地磁场信息,该特点可用于反映区域地磁场精细分布研究,甚至磁异常分布.

(3) 根据建模结果,3D Spline模型在不同高度处的分布主要依赖接近该高度处的实测点数量和精度,如地面的网格分布由于测点多而更为曲折,而高空及卫星高度处由于测点少,分布较为平滑.而在缺少测点的高度处,模拟结果仍不准确.

本研究在高空和卫星高度选用的测点相对较少,这样做的目的是为了突出地面测点的重要性.今后若可获取2000.0年以后的中国地区地面复测点数据,以及航磁和SWARM卫星数据,有望进一步提高3D Spline模型的模拟精度.

目前有多种进行区域地磁场三维建模的方法,如球冠谐模型、矩谐模型等,本研究提供的方法具有使用简便、模拟精度高等特点,可作为基于海量卫星数据而开展地磁场区域磁场研究的新选择.

致谢

感谢空间天气学国家重点实验室开放课题,中国地震局地球物理研究所以及德国地球科学研究中心对研究提供的帮助.

参考文献
Alldredge L R. 1981. Rectangular harmonic analysis applied to the geomagnetic field. J. Geophys. Res., 86(B4): 3021-3026. DOI:10.1029/JB086iB04p03021
An Z C, Xu Y F. 1981. Methods of computation of geomagnetic field at greater altitude in a local region. Chin. J. Space Sci. (in Chinese), 1(1): 68-73.
An Z C, Xu Y F, Xai G H, et al. 1982. Mathematical method to express the distribution of geomagnetic field and its secular variation in a local area. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese), 25(Suppl.): 711-717.
Chen B, Ni Z, Xu R G, et al. 2016. The geomagnetic field in China and neighboring regions for the 2010.0 epoch. Chinese J. Geophys. (in Chinese), 59(4): 1446-1456. DOI:10.6038/cjg20160425
Feng Y, An Z C, Sun H, et al. 2010a. A study on model of geomagnetic normal field of China Region. Advances in Earth Science (in Chinese), 25(7): 723-729.
Feng Y, Pan J J, An Z C, et al. 2010b. Calculation and analysis of geomagnetic field horizontal gradients in China. Chinese J. Geophys. (in Chinese), 53(12): 2899-2906.
Feng Y, Sun H, Mao F, et al. 2013. Study on crustal magnetic anomaly over Mexico and its adjacent areas. Progress in Geophys. (in Chinese), 28(3): 1355-1364. DOI:10.6038/pg20130329
Feng Y, Sun H, Jiang Y. 2015. Data fitting and modeling of regional geomagnetic field. Applied Geophysics, 12(3): 303-316. DOI:10.1007/s11770-015-0500-6
Feng Y, Jiang Y, Jiang Y, et al. 2016. Regional magnetic anomaly fields:3D Taylor polynomial and surface spline models. Applied Geophysics, 13(1): 59-68. DOI:10.1007/s11770-016-0533-5
Finlay C C, Olsen N, Tøffner-Clausen L. 2015. DTU candidate field models for IGRF-12 and the CHAOS-5 geomagnetic field model. Earth, Planets and Space, 67: 114. DOI:10.1186/s40623-015-0274-3
Flechtner F M, Gruber T, Güntner A, et al. 2010. System Earth via Geodetic-Geophysical Space Techniques. Heidelberg, Berlin: Springer-Verlag.
Gao J T, An Z C, Gu Z W, et al. 2006. Distribution of the geomagnetic field and its secular variations expressed by the surface Spline method in China (a part) for 1900-1936. Chinese J. Geophys. (in Chinese), 49(2): 398-407.
Haines G V. 1985. Spherical cap harmonic analysis. J. Geophys. Res., 90(B3): 2583-2592. DOI:10.1029/JB090iB03p02583
Jiang Y, Jiang Y, Feng Y, et al. 2015. Regional geomagnetic modelling based on the magnetic data of CHAMP satellite and 3D Taylor polynomial. Chinese J. Geophys. (in Chinese), 58(9): 3121-3132. DOI:10.6038/cjg20150909
Le Mouel J. 1969. Sur la distribution des elements magnetiques en France[Ph. D. thesis]. University de Paris.
Liu S J, Zhou X G, Sun H, et al. 2011. The three dimension Taylor polynomial method for the calculation of regional geomagnetic field model. Progress in Geophys. (in Chinese), 26(4): 1165-1174. DOI:10.3969/j.issn.1004-2903.2011.04.005
Olsen N, Hulot G, Lesur V, et al. 2015. The Swarm Initial Field Model for the 2014 geomagnetic field. Geophys. Res. Lett., 42(4): 1092-1098. DOI:10.1002/2014GL062659
Sabaka T J, Olsen N, Langel R A. 2002. A comprehensive model of the quiet-time, near-Earth magnetic field:phase 3. Geophys. J. Int., 151(1): 32-68. DOI:10.1046/j.1365-246X.2002.01774.x
Sabaka T J, Olsen N, Purucker M E. 2004. Extending comprehensive models of the Earth's magnetic field with Ørsted and CHAMP data. Geophys. J. Int., 159(2): 521-547. DOI:10.1111/gji.2004.159.issue-2
Sabaka T J, Olsen N, Tyler R H, et al. 2015. CM5, a pre-Swarm comprehensive geomagnetic field model derived from over 12 yr of CHAMP, Ørsted, SAC-C and observatory data. Geophys. J. Int., 200(3): 1596-1626. DOI:10.1093/gji/ggu493
Thébault E, Finlay C C, Beggan C D, et al. 2015. International Geomagnetic Reference Field:the 12th generation. Earth, Planets and Space, 67: 79. DOI:10.1186/s40623-015-0228-9
Wang C. 2010. New chains of space weather monitoring stations in China. Space Weather, 8(8). DOI:10.1029/2010SW000603
Xu W Y. 2009. Physics of Electromagnetic Phenomena of the Earth (in Chinese). Hefei: University of Science and Technology of China Press.
安振昌, 徐元芳. 1981. 局部地区高空磁场的计算方法. 空间科学学报, 1(1): 68–73.
安振昌, 徐元芳, 夏国辉, 等. 1982. 表示局部地区地磁场及其长期变化分布的数学方法. 地球物理学报, 12(增): 711–717.
陈斌, 倪喆, 徐如刚, 等. 2016. 2010.0年中国及邻近地区地磁场. 地球物理学报, 59(4): 1446–1456. DOI:10.6038/cjg20160425
冯彦, 安振昌, 孙涵, 等. 2010a. 中国地区地磁正常场模型研究. 地球科学进展, 25(7): 723–729.
冯彦, 潘剑君, 安振昌, 等. 2010b. 中国地区地磁场水平梯度的计算与分析. 地球物理学报, 53(12): 2899–2906.
冯彦, 孙涵, 毛飞, 等. 2013. 墨西哥及邻近地区地壳磁异常场研究. 地球物理学进展, 28(3): 1355–1364. DOI:10.6038/pg20130329
高金田, 安振昌, 顾左文, 等. 2006. 用曲面Spline方法表示1900-1936年中国(部分地区)地磁场及其长期变化的分布. 地球物理学报, 49(2): 398–407.
蒋勇, 姜乙, 冯彦, 等. 2015. 基于CHAMP卫星与三维Taylor多项式模型的区域地磁建模研究. 地球物理学报, 58(9): 3121–3132. DOI:10.6038/cjg20150909
柳士俊, 周小刚, 孙涵, 等. 2011. 用于区域地磁场模型计算的三维Taylor多项式方法. 地球物理学进展, 26(4): 1165–1174. DOI:10.3969/j.issn.1004-2903.2011.04.005
徐文耀. 2009. 地球电磁现象物理学. 合肥: 中国科学技术大学出版社.