区域航磁异常汇编数据的长波成分提取、评价与替换技术

张攀， 杜劲松,2,3*， 王震， 杨明豫， 陈超,2

1 中国地质大学(武汉)地球物理与空间信息学院， 武汉 430074 2 地球内部多尺度成像湖北省重点实验室， 武汉 430074 3 地质过程与矿产资源国家重点实验室， 武汉 430074

0 引言

地磁场是地球固有的天然物理场之一，其在固体地球科学、海洋科学、环境科学、资源勘探、工程勘察以及军事与国防等诸多领域发挥着重要作用(管志宁, 2005; Hinze et al., 2013; 杨华和梁月明, 2013; 熊盛青等, 2014, 2016; Li and Wang, 2018).在众多的地磁场测量方式中，航空磁测相比地面磁测具有较高的测量效率，且不受水域、森林、沼泽、沙漠和高山等特殊地貌和复杂地形的限制，同时由于飞行是在距地表一定的高度(如50 m～5 km)进行的，从而减弱了地表磁性不均匀体与环境干扰等影响(熊盛青, 2020)；此外，相对卫星磁测(Zhou et al., 2018)而言，航空磁测由于距离岩石圈较近，对于岩石圈磁场的中、短波长成分具有较高的敏感度与测量精度.因此，各个国家与地区的区域岩石圈磁场调查一般均采用航空测量方式(Ravat et al., 2009; 熊盛青等, 2013; Xiong et al., 2016; Golynsky et al., 2018; Djomani et al., 2019).

但是，各个国家和地区编制的航磁异常图及其数据集一般均是由多个数据子集汇编而成，这些子数据集之间的测量时间不同、仪器与平台各异，具有不同的飞行高度、磁场测量分辨率和精度以及测点定位精度，而且不同时期采用的数据校正方式以及相关模型(例如正常场校正所采用的主磁场模型的分辨率和精度各异)和参数也可能不相同.即使上述问题均不存在，由于地球主磁场随时间变化，一方面岩石圈感应磁化强度及感应磁场也将发生变化(Jackson, 2007; Thébault et al., 2009; Sebera et al., 2019)，另一方面即使岩石圈磁性结构不发生任何变化，总磁场强度异常即ΔT磁力异常也会由于主磁场变化而发生改变(袁晓雨等, 2015; 孙石达等, 2020).因此，区域航磁异常汇编数据往往存在低频或长波成分不可靠性(傅敖云, 1984; Ravat et al., 2003).这虽然对于岩石圈浅部的、局部的磁性体与地质构造等研究的影响较小，但是对于理解深部的、区域性岩石圈磁性构造格架及其大地构造与动力学意义具有较大影响.

相比之下，卫星磁测可以在较短的时间内比较均匀地覆盖全球，对地磁场的测量精度与分辨率也在逐步提升，随着数据处理与建模技术的快速发展，目前单纯依赖卫星数据构建的全球岩石圈磁场球谐模型可靠阶数可达90阶(Olsen et al., 2017)，对应空间分辨率为2°(在赤道地区约223 km).因此，采用这些模型的可靠波段对区域航磁异常汇编数据的长波成分进行评价与替换是提高区域磁力异常数据汇编质量的关键步骤之一.显然，构建“全波段”均可靠的区域磁力异常汇编数据，作为一个国家或地区的基础地球物理场数据，将能更好地服务于教学、科研以及各个行业和领域的广泛应用(例如辅助导航, 郑晖等, 2012).

对于全球数据汇编，可以首先对其进行球谐展开，然后将相应长波波段的球谐系数采用基于卫星磁测数据构建的全球岩石圈磁场模型的球谐系数替代(Maus et al., 2009; Lesur et al., 2016; 杜劲松等, 2017).对于区域数据汇编，尽管可以采用(超)长剖面高精度磁测数据进行调平处理，或者选择编图区内部分切割线与测线组成若干磁场调平框架从而调整各个测区的磁场水平(尹航等, 2015)，但是更为高效的方法还是在长波波段采用基于卫星磁测数据构建的全球岩石圈磁场模型预测数据进行替换.

为了准确且唯一地提取区域航磁异常数据的长波成分，需要一种基函数完备、正则、调和且正交的区域岩石圈磁场建模方法，而修正球冠谐分析方法(Thébault et al., 2004, 2006a,b)正是满足这些条件且是目前较为成熟的一种区域地磁场建模方法.该方法在地球外部一个球面锥体之内对地磁场进行建模，相比于传统球冠谐分析方法(Haines, 1985a,b)，修正球冠谐分析方法考虑到两组边界条件而非一组，这两组边界条件引入了外源Legendre函数序列与Mehler函数序列，前者改善了卫星高度数据的拟合与预测效果，后者改善了数据空缺区的上、下延拓效果.相比其它方法，R-SCHA方法具有计算量小、计算精度高的明显优势，不仅可以融合地面、船载、航空与卫星等多来源的地磁场不同分量(例如：总磁场强度、磁场三分量、磁场梯度等)的观测数据，而且在分量转换以及空间延拓等方面具有较高的计算精度(Ou et al., 2013; Vervelidou et al., 2018).因此，该方法已被广泛应用于区域性主磁场、岩石圈磁场、电离层与磁层磁场等的建模之中(徐文耀等, 2011; Talarn et al., 2017; Torta, 2020).

为了对航磁异常汇编数据中不可靠的低频数据进行提取、评价和替换，本文选用二维修正球冠谐分析(R-SCHA2D)方法(Thébault, 2008)，首先对航磁异常数据进行R-SCHA2D反演建模，得到修正球冠谐系数，之后根据模型系数与波长之间的关系提取长波成分，最后使用全球岩石圈磁场模型的长波成分对区域航磁异常汇编数据中的长波成分进行评价与替换.本文首先交代了R-SCHA2D的方法原理与计算算法，然后给出了区域航磁异常汇编数据长波成分提取、评价与替换技术，最后通过仿真测试与实际应用证明了该方法的有效性与实用性.

1 修正球冠谐分析方法

1.1 方法原理

修正球冠谐分析方法的建模区域Ω是由地心距为a和b的两个曲面S1与S2以及以地心为顶点、半张角为θmax的球面椎体包围而成(如图1a所示).对于近似分布在同一高度上的数据(例如航磁异常汇编数据)，建模区域上、下球冠面近似重合，此时修正球冠谐分析方法中Mehler函数的定义不再成立.针对此问题，Thébault(2008)提出了二维修正球冠谐分析方法，其建模范围如图1b所示.

当建模区域的上、下球冠重合时，通过推导可以得到磁位异常在无源空间的级数表达形式(Thébault, 2008)：

(1)

图1 三维(a)与二维(b)修正球冠谐分析的建模区域示意图Fig.1 Sketch map of 3D (a) and 2D (b) modelling regions for revised spherical cap harmonic analysis

进而，根据磁位异常可以得到磁场三分量异常：

(2)

航磁异常数据一般均为总磁场强度异常(ΔT)数据，而磁力异常矢量的模量一般远小于主磁场强度，因此可将总磁场强度异常视为磁场三分量异常在主磁场方向上的投影(Blakely, 1995; 管志宁, 2005; Hinze et al., 2013)，即

ΔT=ΔBxcosIcosD+ΔBycosIsinD+ΔBzsinI，(3)

其中，I与D分别为地球主磁场方向的倾角和偏角.

1.2 坐标与分量转换计算方法

在模型反演计算之前，首先需要将地理坐标与地磁场三分量异常从球坐标系转换至球冠坐标系(安振昌, 1992; De Santis et al., 1996)，坐标与分量异常的转换过程分为如下三步：

首先，通过式(4)将坐标从球坐标系(φ,θ,r)转换到地心直角坐标系(x,y,z)：

(4)

其次，通过式(5)对笛卡尔坐标系进行旋转，使得笛卡尔坐标系Z轴指向球冠中心：

(5)

其中，θ0与φ0为球冠中心的余纬度和经度.最后，通过式(6)将坐标从地心直角坐标系转换到球冠坐标系中：

(6)

地磁场三分量异常可通过式(7)从球坐标系转换至球冠坐标系：

(7)

其中

(8)

(9)

1.3 基于ΔT磁力异常数据的R-SCHA2D模型系数求解方法

式(1)与(2)的二维修正球冠谐级数表达式可被描述为一个线性系统：

(10)

其中，Gx、Gy、Gz分别为磁场三分量异常的二维修正球冠谐核函数矩阵，m为二维修正球冠谐系数矩阵.从式(3)可以看出，总磁场强度异常可视为磁场三分量异常的线性组合，将式(7)代入式(10)中，首先将球冠坐标系中地磁场三分量异常的核函数转换至球坐标系：

(11)

其次，根据式(3)对地磁场三分量异常的核函数进行线性组合，即可得到基于总磁场强度异常数据的修正球冠谐系数解算方程：

=ΔT，(12)

上标T表示矩阵的转置.由于观测数据不可避免地包含噪声，因此将式(12)改写为

(13)

其中数据加权矩阵为

(14)

式中εi为第i个ΔT磁力异常观测数据的标准差.式(13)为一个线性方程组，可以采用共轭梯度法进行求解，但是在实际反演计算中，系数矩阵的条件数往往很大，会严重影响到共轭梯度算法的收敛性，因此本文采用预优共轭梯度法(Pilkington, 1997)以改善系数矩阵的条件数，通过预处理使系数矩阵的特征值较集中分布，从而提高了迭代收敛速度及计算稳定性.

1.4 航磁异常数据的补空与扩边方法

对于实际的航空磁力异常数据，往往存在数据空白区以及不规则的区域边界，若不对其进行处理而直接进行建模，则在数据空白区可能构建出虚假异常以及在边界区域致使振荡异常出现.压制这些现象最直接与最有效的方法即是对数据空区进行补值以及对区域边界外围进行扩边处理(安振昌, 2003).本文采用全球岩石圈磁场高阶球谐模型，例如790阶/次的EMM2017球谐模型(www.ngdc.noaa.gov/geomag/EMM/)或800阶/次的WDMAMv2球谐模型(www.wdmam.org)，通过球谐解算得到与航磁异常汇编数据相同高度的ΔT磁力异常数据(陈康等, 2021)，从而完成对区域航磁异常汇编数据的空区补值和边界区域扩边处理.

2 长波成分的提取、评价与替换技术

对于区域航磁异常汇编数据长波成分的提取、评价与替换，主要按照如下五个步骤进行：

第一步，航磁异常数据预处理，主要包含三个方面：一是，对于实际航磁异常数据需要采用全球高阶岩石圈磁场模型进行空区补值与外围扩边处理，在目标区域外围，若没有足够扩边的数据，则会由于边界效应将影响目标区域内部的建模精度；二是，实际航磁异常数据的分辨率可能较高，考虑到长波成分提取的目标任务以及计算量与计算效率问题，区域磁场模型无须拟合高频成分，因此可以预先对补值与扩边之后的磁力异常数据进行一定程度的降尺度和平滑处理，作为建模输入磁力异常数据；三是，合理评估建模输入磁力异常数据的噪声，即每个测点磁力异常数据的标准差，若无法评价，则将数据加权矩阵每个对角线元素均设置为1 nT-1.

第三步，航磁异常长波成分的提取与评价.由于式(1)所示的基函数系具备正交性与完备性，因而修正球冠谐分析所得的各个系数是唯一的且可以确定的，而且非整数nk阶与波长l具有严格的对应关系(Thébault et al., 2006b)，即

(15a)

(15b)

可见nk与球冠半张角θmax近似呈反比例关系，即球冠面积越小、nk越大，代表的波长则越短.因此，通过截断球冠谐级数的展开阶数(如K阶)而仅保留k≤K的系数，再通过正演计算即可完成长波成分的提取.对于K值的确定，可以采用波长与能量之间的谱分析通过对比航磁异常与卫磁异常的能量谱曲线从而确定.但是为了简化计算，如引言所述，目前单纯依赖卫星数据构建的全球岩石圈磁场球谐模型可靠阶数可达90阶(Olsen et al., 2017)，因此可以直接对比球谐90阶及其以下的能量谱，即令l=90得到K值，从而判断航磁异常长波成分(即k≤K)的可靠性.修正球冠谐能量谱的计算采用Vervelidou和 Thébault (2015)提出的方法，其内源与外源勒让德函数序列在球冠面上的面积归一化平均能量如式(16)所示：

(16)

第四步，航磁异常长波成分的替换.通过上步对比航磁异常与卫磁异常的能量谱曲线，确定截断波长，再采用卫磁异常的长波成分替换航磁异常的长波成分.

第五步，航磁异常数据长波成分替换之后的可靠性验证.由于卫星测量高度较高(一般高于250 km)，距离岩石圈较远，卫星磁力异常数据主要为长波成分，因此可以采用替换了和未替换长波成分的航磁异常汇编数据(ΔT)，分别与卫星磁测数据联合进行修正球冠谐建模，进而根据卫星磁力异常数据的拟合差情况判断航磁异常数据长波成分替换的必要性，这是因为若两者的长波成分不匹配则难以拟合观测数据.

3 仿真测试

为验证本文所提方法的有效性，本节首先使用全球岩石圈磁场模型通过正演计算仿真模拟航磁异常观测数据与长波成分理论数据，然后对模拟的航磁异常观测数据进行二维修正球冠谐建模进而提取其长波成分，最后将长波成分的理论数据与提取数据进行对比分析.由于澳大利亚航磁异常数据具有较高的分辨率与精度，因此在仿真测试与实际应用中均采用澳大利亚作为研究区域.

3.1 仿真数据

考虑到计算量和计算效率，仅使用EMM2017全球岩石圈磁场模型(www.ngdc.noaa.gov/geomag/EMM/)的16到400阶球谐系数，模拟计算了澳大利亚地区的总磁场强度异常数据作为测试的仿真观测数据，如图2a所示，数据分布范围为东经112°至东经154°、南纬10°至南纬44°，网格间距为0.25°×0.25°，大地高度为500 m；同时，使用该模型的16～90阶的正演计算数据作为长波成分提取结果的检核数据(图2b).

图2 仿真模拟的航磁异常观测数据(a)与长波成分理论数据(b)投影方式为等面积Hammer投影，中心经线设置为东经133°；色标左端的NaN表示缺失数据区域.Fig.2 Simulated observation data (a) of aeromagnetic anomaly and theoretical data (b) of long wavelength componentsHammer equivalent projection is adopted. Central meridian is set to be 133°E. The NaN in the left side of the color bar denotes the data gap.

3.2 测试结果及其分析

将球冠的中心设置为东经133°、南纬27°，球冠的半张角设置为27°；将Legendre函数序列的最大截断阶数设为50阶，Mehler函数序列的最大截断阶数设为5阶.对图2a所示航磁异常观测数据进行二维修正球冠谐建模，所得模型预测数据与拟合残差分别如图3a与3b所示，提取出的小于与等于球谐90阶的低频成分数据及其与检核数据之间的残差分别如图3c与3d所示.

由图3a与3b可以看出，模型较好地拟合了观测数据，仅少部分高频异常数据受限于较低的展开阶数而没有得到拟合，一方面展开阶数太高将大幅增加计算量，另一方面部分高频信号的未拟合不会影响长波成分的提取质量.由图3c与3d可以看出，从仿真数据中提取出的低频成分在分布形态上与检核数据一致，残差的最大值为77.72 nT、最小值为-47.68 nT、平均值为0.26 nT、均方差为±7.74 nT，低频成分的能量在分离中仅具有轻微的泄露，这说明本文提出的长波成分提取技术具有有效性.

图3 模型预测数据(a)及其拟合差(b)、提取的长波成分数据(c)及其误差(d)Fig.3 Predicted data (a) and misfit (b) of the model, extracted data of long wavelength components (c) and its errors (d)

图4为模拟航磁异常(EMM2017模型16～400阶/次)与理论长波成分(EMM2017模型16～90阶/次)的球冠谐频谱曲线.由于建模区域的半张角为27°，因此能谱曲线球冠谐阶数的区间大小为180/27、初始阶数为90/27.由图4可以看出，理论长波成分与航磁异常低于球谐90阶的长波成分的频谱非常一致，这进一步说明本文的长波长成分提取方法是可靠的，此外还说明对航磁异常数据进行降尺度处理不会影响其长波成分的提取.

图4 模拟航磁异常与理论长波成分的频谱曲线Fig.4 Spectrum curves of simulated aeromagnetic and theoretical long-wavelength magnetic anomaly fields

4 实际应用

为了进一步说明本文所提方法在实际应用中的有效性以及必要的数据处理步骤，本节将长波长成分提取、评价与替换技术应用于澳大利亚地区的航磁异常汇编数据.

4.1 澳大利亚航磁异常汇编数据

如图5所示，本文所使用的航磁异常数据为第七代澳大利亚总磁场强度异常网格数据(Djomani et al., 2019).该网格数据的网格间距约为3弧秒(约80 m)，相较于第六代澳大利亚总磁场强度异常网格数据，新加入了234组航磁测量数据，且该网格数据已被调平到地形之上80 m的高度，网格数据所在的坐标系为GDA94坐标系.

图5 澳大利亚航磁异常数据(a)与陆地高程(b)Fig.5 Aeromagnetic anomaly data (a) and land elevation (b) of theAustralia

4.2 数据处理

由于原始数据的数据量较大，在局部区域具有上万nT的异常值存在，且存在数据空缺区(图5a)，因此在使用修正球冠谐分析方法建模之前，需要对原始数据进行以下四步的处理：

(1)考虑到计算量，首先将原始数据降尺度为0.05°×0.05°间隔的网格数据；

(2)使用EMM2017全球岩石圈磁场模型的16～790阶，对数据空区进行补值，如图6a所示；

(3)经过降采样处理后，部分高幅值异常变成了奇异点，故对降采样之后的数据再进行基于滑动球冠的高斯滤波，高斯滤波器的滤波半径设置为100 km，滑动球冠的半张角为0.25°，对高斯滤波后的0.05°×0.05°网格间距数据(图6b)进行均值降采样到0.25°×0.25°网格间距(图6c)，并计算了各个网格点的标准差(图6d)，作为反演建模的观测数据误差；

(4)将数据从GDA94大地坐标系转换到球冠坐标系.

图6 澳大利亚航磁异常数据处理结果(a) 降采样为0.5°×0.5°网格间距的数据扩边结果； (b) 高斯平滑结果； (c) 高斯平滑滤除的高频异常；(d) 0.25°×0.25°滑动窗口统计计算得到的高频异常的标准差分布.Fig.6 Processed results of aeromagnetic anomaly data of the Australia(a) Extended result based on the downscaled grid data with resolution of 0.5°×0.5°; (b) The result by Gaussian smooth; (c) The filtered high-frequency anomalies by Gaussian smooth; (d) Standard deviation map of the filtered high-frequency anomalies by statistically calculation in a moving window with a size of 0.25°×0.25°.

4.3 长波成分的提取、评价与替换

对图6b所示航磁异常数据进行二维修正球冠谐分析，参数设置与仿真测试的相同，提取出的小于和等于球谐90阶的长波成分与短波成分数据分别如图7a与7b所示.前已述及，由于没有观测数据的约束，在区域磁场建模时，数据空区会产生虚假异常，进而在处理转换中引入较大误差，故对数据空区进行补值是有必要的，但是由于不同数据集之间的差异，补值同样也会在数据中引入错误的低频成分，这可以明显地从图6a中看出，在数据缝合线两边的数据存在明显的系统差异，从分离出的长波长成分(图7a)中也可以看出，对数据进行补值也引入了错误的长波长成分.除了在补值区域与航磁异常数据之间的缝合区域之外，对比图2b与图7a可知，由航磁数据提取的长波长成分与EMM2017全球岩石圈磁场模型中的低频成分(16～90阶)的空间分布形态基本一致，但是磁力异常的幅值变化存在较大差异.

图7 提取的澳大利亚航磁异常的长波成分(a)与短波成分(b)Fig.7 Extracted data sets of long (a) and short (b) wavelength components from aeromagnetic anomaly data of the Australia

由航磁数据提取的长波长成分(图7a)与EMM2017全球岩石圈磁场模型中的低频成分(16～90阶)(图2b)在幅值方面存在较大差异.但是，该差异主要体现在航磁异常数据区域，而在补值与扩边区域差异微弱，这从侧面说明计算方法是正确的、计算结果是可靠的.EMM2017全球岩石圈磁场模型是在EMAG2v3全球网格数据(Meyer et al., 2017)的基础上进行球谐分析构建的，其16～133阶采用了MF7模型(MF6模型的升级版, Maus et al., 2008)进行了替换.为了说明MF7模型16～90阶的可靠性，本文又采用最新的全球岩石圈磁场模型GRIMM_L120(Lesur et al., 2013)的16～90阶、LCS-1(Olsen et al., 2017)与CM6(Sabaka et al., 2020)的16～90阶、CHAOS7.6(Finlay et al., 2020)的21～90阶球谐系数解算了澳大利亚地区航空高度的磁力异常，如图8所示，其统计参数见表1.通过比较可以看出，五个模型解算的长波磁力异常无论是在空间分布形态还是幅值变化方面均不存在较大的差异，这说明图7a所示的航磁异常长波成分自身是不可靠的.

图8 不同全球岩石圈磁场模型解算的澳大利亚磁力异常长波成分(a) GRIMM_L120模型的16～90阶； (b) LCS-1模型的16～90阶； (c) CM6模型的16～90阶； (d) CHAOS7.6模型的21～90阶.Fig.8 Magnetic anomaly data sets with long wavelength components over the Australia from different global lithospheric magnetic field modelsGRIMM_L120 (a), LCS-1 (b) and CM6 (c) models with spherical harmonic (SH) degrees 16～90,and CHAOS7.6 (d) model with SH degrees 21～90.

图9为卫磁异常、原始航磁异常与替换过低频的航磁异常的频谱曲线，可见未替换过低频的航磁异常与卫磁异常在球冠谐30阶以下的低频成分存在较大差异，鉴于卫磁异常的可靠阶数可达90阶，因此证明原始航磁异常中的低频成分是不可靠的，尤其是30阶以下的波长成分.在90阶以上，单纯基于卫星磁测数据构建的全球岩石圈磁场模型(如GRIMM_L120、LCS-1、CM6和CHAOS7.6)由于空间分辨率有限导致其能量低于EMM2017和航磁异常的能量，EMM2017模型由于联合了卫磁和航磁异常数据因而其能谱曲线与航磁异常的能谱曲线比较一致，而在200阶以上，EMM2017的能量略高于航磁异常的能量，这是由于在提取航磁异常长波成分时对原始数据进行了降尺度和高斯平滑处理，但是这并不影响航磁异常长波成分提取的可靠性.此外，通过图9还可以发现，在所有单纯基于卫星磁测数据构建的全球岩石圈磁场模型中，CHAOS7.6与LCS-1模型与实际航磁异常数据的一致性最好，尤其是在球冠谐110阶以上的波长成分，但是CHAOS7.6模型未给出16～20阶的岩石圈磁场球谐系数，因而其能量在球冠谐30阶以下偏低.

图9 航磁异常与卫磁异常的频谱曲线Fig.9 Spectrum curves of aeromagnetic and satellite magnetic anomaly fields

此外，由表1可以看出，航磁异常提取的长波成分的平均值与全球岩石圈磁场模型的平均值差异较大，这是由于这些差异主要来源于长波成分，即球冠谐30阶以下的波长部分.

由于CM6模型(Sabaka et al., 2020)采用了Ørsted、SAC-C、CHAMP与Swarm大量的卫星磁测数据，并且采用了综合法进行地磁场建模，各种起源的地磁场之间的耦合性更强，其主磁场与变化磁场也被广泛用于磁力异常数据的归算处理，因此综合考虑之后选取CM6模型的16～90阶球谐系数解算的磁力异常数据对提取的澳大利亚航磁异常长波数据进行替换.具体计算步骤为：

(1)采用构建的二维修正球冠谐模型解算澳大利亚原始航磁异常数据网格节点上的长波成分磁力异常，将其从原始航磁异常数据之中扣除，得到剩余的短波成分磁力异常；

表1 澳大利亚总磁场强度异常长波成分统计参数表Table 1 Statistic parameters of magnetic anomaly data sets with long wavelength components over the Australia

(2)采用CM6模型的16～90阶球谐系数解算澳大利亚原始航磁异常数据网格节点上的长波成分磁力异常，将其加入步骤(1)长波成分提取之后的剩余短波成分磁力异常，即可得到最终的航磁异常数据(图10a)；

(3)若需要对航磁异常空区进行补值以及对外围进行扩边，则可以将EMM2017模型的16～90阶球谐系数采用CM6模型的16～90阶球谐系数进行替换，得到新的16～790阶的全球岩石圈磁场球谐系数之后，再根据空区待补值与待扩边数据的坐标进行磁力异常解算，将其与第(2)步处理结果组合，即可得到空区补值与扩边之后的航磁异常格网数据(图10b).

进一步地，为了验证航磁异常汇编数据长波成分的替换效果，分别采用替换了和未替换长波成分的航磁异常汇编数据(ΔT)，与CHAMP卫星磁测数据(ΔBx、ΔBy与ΔBz)联合进行修正球冠谐建模.其中，所采用的CHAMP卫星磁测数据如图11所示，其处理流程和相关参数与全球岩石圈磁场模型GRIMM_L120(Lesur et al., 2013)的一致.采用相同的建模参数，基于替换了长波成分的航磁异常汇编数据与CHAMP卫星磁测数据联合建模预测的卫星测点处磁力异常及其拟合差如图12所示，而基于未替换长波成分的航磁异常汇编数据与CHAMP卫星磁测数据联合建模预测的卫星测点处磁力异常及其拟合差如图13所示，两种情况之下的相关统计参数见表2与表3.由此可以看出，替换了长波成分的航磁异常汇编数据与卫星磁测数据融合效果更好，而由航磁异常汇编数据提取的长波成分与卫星磁测数据存在不相容性或非一致性.

但是，从图12和表3可以看出，即使替换了航磁异常长波成分，在卫星磁测数据拟合残差中还存在一些趋势信号，笔者多次调试航磁异常与卫磁异常的拟合权重，但是依然无法拟合这些残存信号.Vervelidou等(2018)认为这些残余信号来源于卫磁数据处理中未完全消除的外源场(例如电离层和磁层)信号，且分量之间在物理上存在不耦合性，因此在建模时无法完全拟合.为了验证这种观点，笔者又计算了GRIMM_L120、LCS-1、CM6和CHAOS7.6模型对卫磁数据的拟合情况，发现这些残余信号的确存在，如图14所示.对比图12与图14可以发现，本文所构建模型的拟合效果更好，这是由于航磁异常90阶以上的中短波长成分的可靠性更高.除此之外，笔者认为这些无法拟合的趋势信号还有可能来源于主磁场与岩石圈磁场之间的混叠效应.

表2 航磁数据与卫星磁测数据联合建模的模型预测数据统计参数表Table 2 Statistic parameters of model predictions by joint modelling of aeromagnetic data and satellite magnetic measurements

表3 航磁数据与卫星磁测数据联合建模的数据拟合残差统计参数表Table 3 Statistic parameters of data misfits by joint modelling of aeromagnetic data and satellite magnetic measurements

5 结论

针对区域航磁异常汇编数据可能存在长波成分的不可靠性，本文基于二维修正球冠谐分析(R-SCHA2D)以及单纯基于卫星磁测数据构建的全球岩石圈磁场模型提出区域航磁异常汇编数据的长波成分提取、评价与替换技术.首先，使用EMM2017全球岩石圈磁场模型计算了仿真数据与检核数据，仿真测试的结果证明了该方法可以正确地分离出数据中的低频成分；然后，将本文提出的方法应用于第七代澳大利亚总磁场强度异常网格数据，采用单纯由卫星磁测数据构建的可靠长波成分替换了由数据补值与扩边以及原始航磁异常数据中自带的一些不可靠的长波成分，得到了各个频段均更可靠的澳大利亚及邻区的航磁异常数据；最后，采用替换了和未替换长波成分的航磁异常汇编数据，与CHAMP卫星磁测三分量数据联合进行修正球冠谐建模，结果显示，替换了长波成分的航磁异常汇编数据与卫星磁测数据融合效果更好，而由航磁异常汇编数据提取的长波成分与卫星磁测数据存在不相容性或非一致性.仿真测试与实际应用均证明了本文所提的区域航磁异常数据长波成分提取、评价与替换技术具有有效性和实用性，因此可以在区域航磁异常数据汇编中发挥重要作用.

图10 最终合成的澳大利亚航磁异常数据(a) 空区未补值与未扩边； (b) 空区补值与扩边.Fig.10 Finally compiled aeromagnetic anomaly data of the Australia(a) The result of no data supplement and no extension; (b) The result after data supplement and extension.

图13 使用未替换过低频成分的航磁数据与卫星数据联合建模的模型预测数据与拟合残差(a)与(b)分别为航磁异常(ΔT)预测数据与拟合残差； (c)与(d)分别为卫星磁测北向分量(ΔBx)的预测数据与拟合残差； (e)与(f)分别为卫星磁测东向分量(ΔBy)的预测数据与拟合残差； (g)与(h)分别为卫星磁测径向分量(ΔBz)的预测数据与拟合残差.Fig.13 Predictions and data misfits of the model by joint modelling of aeromagnetic anomaly data whose long wavelength components are not replaced and satellite magnetic measurements(a) and (b) are predicted data and data misfit of aeromagnetic anomaly, respectively; (c) & (d), (e) & (f), and (g) & (h) are predicted data and data misfit of northern component (ΔBx), eastern component (ΔBy) and radial component (ΔBz) of satellite magnetic anomaly, respectively.

图14 CM6模型对卫星磁测数据的拟合差分布(a) 北向分量(ΔBx)； (b) 东向分量(ΔBy)； (c) 径向分量(ΔBz).Fig.14 Data misfit distributions of satellite magnetic data by using CM6 model(a) Northern component (ΔBx), eastern component (ΔBy) and radial component (ΔBz).

致谢感谢Geoscience Australia提供澳大利亚航空磁力异常汇编数据！球面投影图件采用了Generic Mapping Tools(GMT)软件进行绘制，两位匿名审稿专家也为本文提出了宝贵的修改建议，在此一并表示衷心感谢！