Fusion of Ground and Airborne Magnetic Data Using Multi-Layer Equivalent Source Method
-
摘要: 随着地磁场观测数据的不断积累，高效利用观测数据成为亟待解决的问题.以往研究表明，仅对单一观测手段获得的实测数据进行分析与解释，往往很难达到目前解决相关地质问题的精度要求.因为各观测方法获得的地磁场数据集通常在分辨率、精度、高程及覆盖范围方面存在局限性和差异性，造成单一数据集仅能有效表征地磁场某一频段信息的问题.而解决该问题的一种有效途径是数据间的融合.为此，基于等效源方法，提出一种多层等效源技术方案，应用于航空和地面磁测数据融合，提高地面数据插值补空、扩边及航空数据下延的精度.该方法针对观测信息的频谱特征，采用3个位于不同深度的等效源层模拟实测数据；较传统的单层等效源方法，减少了等效源设置的盲目性，增强了观测信息在等效源模型中分配的有序性和结构性.理论实验表明，多层等效源模型设置具有更高的计算精度，航空与地面磁测数据融合可以起到显著的相互丰富及改善的作用.最后，将该方法应用于湖北金牛火山岩盆地航空与地面磁测数据融合，获得了丰富的、平面上规则分布的地磁场数据.Abstract: With the accumulation of measured magnetic data, it is becoming urgent to use these data efficiently. Previous studies have shown that it is difficult to meet the accuracy requirement for solving geological problems by using the measured data from a single observation method. Because of the limitations and differences in resolution, accuracy, elevation and range of magnetic data obtained by various methods, single dataset can only effectively reflect the information over a certain range of wavelength of magnetic field. An effective way to solve this problem is the fusion of data. Therefore, based on the equivalent source method, a multi-layer equivalent source technology is proposed in this paper, which can be applied to the fusion of ground and airborne magnetic data to improve the accuracy of interpolation, continuation, extension, transformation, et al. According to the spectral characteristics of observation data, three-layer equivalent sources at different depths are used to fit the measured data. Compared with the traditional single-layer equivalent source method, it can reduce the blindness for setting equivalent sources, and improve the ordering and structural performance for allotting observation information into equivalent sources. Synthetic experiment shows that the three-layer model has higher computational accuracy, and data fusion can significantly improve each dataset. Finally, the method is applied to the fusion of ground and airborne magnetic data in Jinniu basin, Hubei, and abundant magnetic data with regular distribution on a plane are obtained.
-
Key words:
- data fusion /
- multi-layer equivalent source /
- interpolation /
- extension /
- data transformation /
- geophysics
-
0. 引言
地磁场是地球的天然物理场之一，其中蕴含着地球内部磁性与温度结构等丰富信息，能够为地球构造运动及地质演化过程提供重要指示.因此，地磁场数据在矿产资源勘查、环境监测和地质调查等领域具有不可替代的作用.随着测量设备与技术的发展，利用星载、机载、船载、钻井及人力等观测手段获得了大量不同高度、不同分辨率和不同分布范围的地磁场测量数据（Stolz et al., 2006；Clark，2013；杜劲松和陈超，2015；Zhou et al., 2015；张婉等，2018），为研究地磁场和相关地球科学问题提供了宝贵的数据基础.但是，各观测手段获得的数据没有统一标准，而且分辨率单一、观测高度与观测范围有限，仅能表征地磁场某一频段范围的信息（Maus et al., 2009）.随着对资料解释准确性要求的不断提高，利用单一技术手段获得的位于孤立平面或曲面上的实测数据，无论是所包含的地磁场信息量还是覆盖范围，已逐渐不能满足需求.因此，针对大量的实测数据，如何有效加以利用是一个具有实际意义的研究课题.当前，采用数学物理方法将各数据集进行融合不失为一种高效且经济的策略；通过融合能够扩充数据的频段范围，各数据集可在观测范围和观测高度方面相互弥补.为此，国内外学者提出了许多融合技术方法，也取得了较好的研究成果.
地磁场数据融合的方式通常是基于各个独立的观测数据集建立地磁场模型，将各数据集中包含的地磁场信息反映到所建立的模型中，再通过模型正演过程恢复地磁场以实现数据的融合.一种简洁的方法是将磁场展开成函数的形式，通过实测数据解算函数中相应的系数以实现对磁场的描述，如球谐函数、球冠谐函数和矩谐函数等（Li et al., 1995；Ou et al., 2013）.另一种统计类方法也是地磁场建模的传统方法，包括最小二乘、最小二乘配置和输入-输出方法（Barzaghi et al., 2008；高新兵等，2013）.再者，如Stockes、Molodensky、Hotine等积分类方法（Featherstone，2003；Sjöberg，2005）、基函数方法（吴怿昊等，2016）等，也被用于多源观测数据的融合以构建重磁场模型（Andrews et al., 2015）.以上方法在一定程度上满足了数据融合的需求，但也存在一定的局限和不足，如仅能解算单一类型的数据、要求数据水平规则网格分布、需要获知数据与噪声之间的互功率谱和协方差函数或计算精度较低等，难免影响其实用性.
除此之外，等效源方法也是一种有效且常用的数据融合技术手段.其利用一组简单的虚拟源（即等效源）代替真实场源进行重磁场建模，通过等效源正演对重磁场信息进行恢复.该方法可以直接处理原始点位上的观测数据，能够联合解算不同类型、高度及分辨率的磁测数据，具有很好的适用性和计算精度.等效源方法最早由Dampney（1969）用于重力场的空间延拓. Bhattacharyya and Chen（1977）提出磁偶层等效源进行磁异常转换.之后该方法被广泛应用于重磁场勘探领域数据的插值与网格化（Cordell and Grauch, 1982）、空间延拓（Oliveira Jr et al., 2013）、曲化平（安玉林等，2013）、梯度和分量转换（Barnes and Lumley, 2011）以及化极（Li et al., 2014）等处理.此外，等效源方法还被应用于球坐标系中卫星磁测数据的处理与解释（Asgharzadeh et al., 2008），包括岩石圈磁场建模（Purucker et al., 2000）、卫星磁测数据曲化平（或高度归一化）（Ravat et al., 1995）、梯度与分量转换（Purucker，1990）、多源数据融合（Kim et al., 2007）、向下延拓和化极（Whaler，1994）等方面.
在磁场测量工作中，实际地面磁测往往存在观测面起伏、空白区、测点非均匀分布以及观测范围有限等问题，对数据进一步的处理分析工作造成不便.本文针对以上问题，采用等效源方法建立物理模型，融合航空磁测数据与地面磁测数据，实现地面数据的曲化平、填补空白区、扩边和上延，以及航空数据的下延；同时利用观测数据进行磁场分量及梯度的转换，为后续分析解释工作提供可靠的、丰富的基础数据.传统的等效源方法一般在地表之下某一固定深度设置单层等效源以简化场源模型，提高计算速度.然而实际场源分布在地表之下较广的深度范围，所产生的磁场响应对应着观测数据中不同频段的信息.虽然布设在近地表的单层等效源能够很好地恢复观测数据，但在长距离空间延拓以及数据转换中会存在较大误差（庞旭林，2012），很难较完整地获取观测数据中所含的信息.因此，本文提出多层等效源模型，用于航空和地面磁测数据的融合，实现了将不同频段、高度和分辨率的数据有序且结构化地分配到等效源模型中，提高了数据融合的有效性，为获得高精度的观测数据扩边、插值补空、延拓及转换结果提供了充分保证.
1. 技术方法
基于位场的等效性，等效源方法通过设置一组简单的虚拟源（即等效源）代替真实场源，拟合观测数据建立等效源模型，进而利用等效源模型正演计算进行重磁场数据的处理与转换.根据位场理论中的基本原理，位场作为空间位置的函数F（x，y，z）可以与场源m（ξ，η，ζ）之间建立如下关系：
$$ \begin{aligned} &F(x, y, z)= \\ &\iiint\limits_{\Omega} A(x, y, z, \xi, \eta, \zeta) m(\xi, \eta, \zeta) \mathrm{d} \xi \mathrm{d} \eta \mathrm{d} \zeta, \end{aligned} $$ (1) 其中，A为积分核函数，表示场与场源之间的物理几何关系.如果F（x，y，z）与m（ξ，η，ζ）均为离散分布，并假设离散位场数据为M个，离散场源单元为N个，则式（1）可以表示为
$$ \begin{aligned} &F\left(x_{i}, y_{i}, z_{i}\right)=\sum\limits_{j=1}^{N} G\left(x_{i}, y_{i}, z_{i}, \xi_{j}, \eta_{j}, \zeta_{j}\right) m\left(\xi_{j}, \eta_{j}, \zeta_{j}\right) \\ &i=1, 2, 3, \cdots M, \end{aligned} $$ (2) 式中，
$$ \begin{aligned} &G\left(x_{i}, y_{i}, z_{i}, \xi_{j}, \eta_{j}, \zeta_{j}\right)= \\ &\iiint\limits_{\Omega_{j}} A\left(x_{i}, y_{i}, z_{i}, \xi, \eta, \zeta\right) \mathrm{d} \xi \mathrm{d} \eta \mathrm{d} \zeta, \end{aligned} $$ (3) Ωj表示离散场源单元的积分区域.式（2）采用矩阵形式可表示为
$$ \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{Gm}}{\rm{, }} $$ (4) 式中d=（dg1，dg2，dg3，…，dai，dai+1，…，daM）为数据向量，其中包括地面数据{dgi}和航空数据{dai}，G={Ggij；Gaij}是M×N阶核矩阵，m=（m1，m2，m3，…，mN）是场源物性向量.
对于磁场数据，常采用的等效源单元类型包括磁偶极子（或磁荷）、磁偶层和均匀磁化长方体（Bhattacharyya and Chan, 1977；安玉林等，2013；Li et al., 2014）.结合本文的研究目的，希望应用等效源模型进行磁场数据延拓和转换时具有更好的保幅性，根据以往研究中对3种等效源单元磁场响应随计算高度衰减特性的对比（李端等，2018），认为采用在深度方向具有一定延伸的长方体作为等效源单元，有助于“保持”更多深部场源信息（谢汝宽等，2015；李端等，2018），能够满足保幅性要求.均匀磁化长方体等效源单元磁场正演计算采用郭志宏等（2004）给出的解析公式.
多层等效源模型是指在传统单层等效源之下增加等效源层.理论上，等效源层数越多越接近真实场源分布，获得的处理结果也越接近理论值，但相应的计算成本会大幅度增加.因此，依据数据中所含信息的频段特点，本文提出三层等效源技术方案（图 1），具体设计思路如下：
图 1 三层等效源模型示意图RLPS表示径向对数功率谱估算；zf, zs, zt分别表示3个等效源层中心埋深Fig. 1. Three-layer equivalent source model（1）第一层等效源设置针对地面观测数据，用于模拟浅层场源产生的磁场信息（以中-短波长信息为主）.该层等效源采用传统单层等效源方法进行设置，等效源位于观测面之下2~6倍点距深度（Dampney，1969；Xia and Sprowl, 1991），其在“吸收”数据短波长信息的同时不易产生“虚假”高频信号.等效源单元水平尺寸与地面测点间距（或平均间距）一致（MacLennan and Li, 2013），垂向厚度等于或略大于水平尺寸.等效源层可以采用平面展布，但当地形起伏剧烈时应设置为平行于地形的曲面.
（2）第二层等效源的设置目的是用于模拟中-深层场源产生的磁场信息（可对应航空数据的中-短波长和地面数据中-长波长信息）.可对地面磁测数据和航空磁测数据进行频谱分析（Syberg，1972），进而估计该层等效源的深度设置参数.本文综合航空数据与地面数据的对数功率谱曲线的中-低频段，获得场源平均深度作为第二层等效源的参考深度.该层等效源单元的水平尺寸与航空磁测数据点距一致，厚度可取不大于6倍水平尺寸的厚度（谢汝宽等，2015）.由于距地表较远，第二层等效源采用水平层展布，以简化构建等效源模型的难度.
（3）第三层等效源则针对航空数据进行设置，以存储深层场源磁场信息；其深度依据航空磁测数据频谱分析估计的深部场源平均深度进行取值.等效源单元水平尺寸可参考观测数据覆盖范围选取第二层等效源单元水平尺寸的3~10倍.如此设置有两个目的，其一是通过增加几何尺寸提高深部等效源的敏感性，避免反演过程中物性集中于浅层等效源；其二是减少等效源单元的总数，从而减少计算量.该层等效源单元的垂向厚度及空间展布形式与第二层相同.
本文所提出的三层等效源模型结构有利于高效且有序地模拟磁测数据的主要频段信息.而如何通过等效源物性分布将磁测数据中的信息反映到各层等效源上是该方法的另一关键问题.在等效源的空间位置和几何形态确定之后，根据式（4）确立的线性关系采用反演算法获得等效源物性分布，可将式（4）改写为：
$$ {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{Gm}} = {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{d}}, $$ (5) 基于本文等效源设置规则，上式中观测数据数量一般会少于等效源单元个数，且核矩阵规模较大，使得矩阵的条件数变大，造成反演过程不稳定.因此，在反演过程中加入阻尼因子λ，如式（6），能够有效稳定反演过程（Tikhonov and Arsenin, 1977）.
$$ \left({{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}} + \lambda \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{m}} = {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{d}}, $$ (6) 式中，I为单位阵.由于GTG是对称正定矩阵，因此可以被分解为RΛRT，Λ为GTG特征值组成的对角矩阵，R则是对应的特征向量矩阵，加入阻尼因子后，Λ被改变为
$$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\tau }}_1} + \lambda }&{}&{}&{}\\ {}&{{\mathit{\boldsymbol{\tau }}_2} + \lambda }&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{\mathit{\boldsymbol{\tau }}_i} + \lambda } \end{array}} \right], $$ (7) τi为特征值，令τmax与τmin表示特征值的最大值和最小值.计算（GTG + λI）的条件数κ可得
$$ \boldsymbol{\kappa}=\left|\frac{\boldsymbol{\tau}^{\max }+\lambda}{\boldsymbol{\tau}^{\min }+\lambda}\right|. $$ (8) 显然，加入阻尼因子后能够有效抑制条件数的增加，但当阻尼因子取值较大时会影响反演计算的收敛速度.可以利用理论模型，采用不同阻尼因子进行试验，根据收敛速度的变化确定最佳阻尼因子取值范围（Li et al., 2020）.
三层等效源分布在不同的深度，在反演时会出现物性值集中在浅层的现象（Silva et al., 2014），深层等效源会失去作用，从而不能实现涵盖数据中各频段信息的目的.即使如前述增加了深层等效源单元的几何尺寸，在改善深层等效源物性分布方面仍有一定不足（李端等，2018）.产生以上现象的原因主要是由于正演核函数元素随距离的衰减，Li and Oldenburg（1996）在反演过程中通过在阻尼项加入深度加权函数来平衡衰减，显著改善了反演结果.此外，Pilkington（1997）采用预条件共轭梯度方法，通过预条件矩阵的方式加入深度加权函数取得了相同的效果.出于简化计算的目的，本文采用第二种方式，即利用深度加权函数构造预条件矩阵P，并加入到式（6）中，得出式（9）.预条件矩阵在此处除了深度加权的作用，还可以对式（6）中左侧的矩阵进行预处理以提高计算效率.
$$ \boldsymbol{P}\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{G}+\lambda \boldsymbol{I}\right) \boldsymbol{m}=\boldsymbol{P} \boldsymbol{G}^{\mathrm{T}} \boldsymbol{d}, $$ (9) 由于本文加入深度加权的目的是针对深层等效源，因此采用变加权的方式构造预条件矩阵，如式（10）所示.
式中，z是等效源深度（图 1），β是加权指数，与磁场数据的衰减有关.当等效源单元的深度与几何尺寸确定后，可以通过测试等效源单元磁场响应随埋深变化的衰减来确定β的取值（Li and Oldenburg, 1996）.
$$ \boldsymbol{P}=\left[\begin{array}{ccccccc} z_{f}^{0} & & & & & & & \\ & \ddots & & & & & & \\ & & z_{f}^{0} & & & & & \\ & & & z_{s}^{\beta} & & & & \\ & & & & \ddots & & & \\ & & & & & z_{s}^{\beta} & & & \\ & & & & & & z_{t}^{\beta} & & \\ & & & & & & &\ddots & \\ & & & & & & & &z_{t}^{\beta} \end{array}\right]. $$ (10) 2. 理论模型实验
为检验方法的应用效果，设计理论模型模拟地形表面和500 m高度平面上（航空测量高度）的实测ΔT磁异常.理论模型如图 2a所示，由6个不同大小和埋深的均匀磁化长方体组成，由于除化极以外，利用等效源方法进行磁异常转换时受场源磁化方向影响较小（李端，2018），因此本文不考虑剩磁与退磁，令模型磁化倾角与偏角为（45°，5°），模型参数见表 1.图 2b显示模拟的起伏地表，其中红框表示地面磁测数据范围，航空磁测数据分布与地形范围一致.在地面磁测数据中设置3个空白区（图 2c），根据观测数据分别建立单层和三层等效源模型，利用等效源方法和两种常规方法（傅里叶变换和克里金插值方法）对观测数据进行插值补空、扩边、曲化平（平面高度为350 m）、延拓、梯度与分量转换处理，对比3种方法的计算精度.通过理论模型正演的ΔT磁异常、磁异常梯度（ΔTx，ΔTy，ΔTz）以及磁场三分量（Hax，Hay，Za）如图 3a，3b与3c所示，其中x表示北方向，y表示东方向，z表示垂向且向下为正.
图 2 理论模型及其地面上观测ΔT磁异常a.组合模型；b.地形高程（红框为地面观测范围）；c.模拟地面观测ΔT磁异常Fig. 2. Synthetic models and their theoretical ΔT anomaly on ground表 1 理论模型参数Table Supplementary Table Parameters of synthetic models模型体 中心埋深
(m)磁化强度
A·m-1）磁化方向 1 1 000 1.5 倾角I=45°
偏角A=5°2 900 2 3 1 600 2.5 4 4 000 4 5 5 500 3 6 6 500 5 
图 3 不同高度正演ΔT磁异常(a~c)和磁异常梯度正演结果(d~i)以及磁场分量正演结果(j~o)a，b，c分别为地面、350 m和500 m高度平面上ΔT计算结果；d~f为地面上磁异常梯度计算结果；g~i为350 m高度平面上相应计算结果；j~l为地面上磁场三分量计算结果；m~o为350 m高度平面上相应计算结果Fig. 3. Theoretical ΔT anomaly at different elevations (a-c), theoretical magnetic gradients (d-i) and theoretical components (j-o) on ground and plane at 350 m三层等效源模型结构设置如下：第一层等效源位于地表之下平行于地表起伏的曲面之上，并作为单层等效源方法的等效源模型使用.第二层与第三层等效源呈水平面展布，深度分别由地面与航空磁测数据对数功率谱进行估计（图 4）.三层等效源的平面分布范围均与航空磁测数据范围一致，各层等效源的设置参数见表 2.
图 4 理论模拟地面与航空磁异常对数功率谱曲线图中红色与蓝色直线为对数功率谱曲线分段拟合线；a，b分别为地面与航空磁测数据计算曲线Fig. 4. Radial logarithmic power spectrum of observed magnetic data表 2 理论模型试验等效源层参数Table Supplementary Table Parameters of equivalent source layers of the theoretical model层号 单元体几何参数
(长×宽×高)(m)等效层顶面深度 展布形态 1 200×200×200 距地表600 m 随地形起伏 2 200×200×600 距0 m平面2 500 m 水平 3 2 000×2 000×6 000 距0 m平面5 500 m 水平 基于深层等效源单元的几何形态，根据衰减曲线拟合程度（图 5）选择预条件矩阵中β的取值为（0，2，2）分别对应由浅至深各层等效源单元.进而利用预条件共轭梯度算法获得三层等效源模型物性分布.
图 5 等效源单元磁场响应随埋深变化的衰减曲线a.对应200 m×200 m×600 m等效源单元；b.对应2 000 m×2 000 m×6 000 m等效源单元Fig. 5. Magnetic attenuation with depths of equivalent source cells用单层和三层等效源模型正演获得地面上规则网格分布的ΔT磁异常、磁异常梯度以及磁场三分量，实现地面磁测数据的插值补空、扩边、航空数据的下延以及梯度与分量转换；计算350 m高度平面上的相应结果，对地面数据进行上延曲化平处理.将观测面近似为一个平面（考虑到高差仅为一倍点距），基于傅里叶变换将航空数据下延至地表平均高程，利用克里金插值方法对地面及下延数据进行统一的网格化，获得地面ΔT磁异常的插值补空与扩边结果；此外，利用克里金方法对地面数据进行插值网格化，并基于傅里叶方法将地面数据和航空数据分别延拓至350 m高度平面，再次采用克里金方法对两个数据集进行统一的插值，进而与等效源方法的曲化平结果进行对比.
本文将计算值与理论值作差，获得各方法计算误差的分布.图 6中显示为ΔT磁异常对应的计算误差，其中本文三层等效源方法的计算误差最小，绝大部分区域的误差在异常平均幅值的1%以内（图 6b）；单层等效源方法的计算结果在边界处存在明显误差（图 6a）；相对于等效源方法，傅里叶-克里金方法的计算误差较大，整体的计算效果较差，认为将起伏观测面近似为水平面进行数据处理是造成误差的主要原因（李端，2018）；350 m高度平面上的误差统计显示了类似的结果，其中图 6f所示在地面观测范围边界处存在的高误差值，原因在于傅里叶变换和克里金插值计算过程中的误差积累，使得350 m高度平面上地面数据的上延结果与航空数据的下延结果不匹配.
图 6 三种方法计算ΔT磁异常与理论值绝对误差a~c分别为单层、三层等效源方法与傅里叶‒克里金方法在地面上的计算结果与理论值差值；d~f分别为350 m高度平面上对应差值Fig. 6. Distribution of absolute errors between the computational ΔT and theoretical values磁异常梯度的计算误差如图 7所示，从图中可以看到三层等效源方法的转换结果整体上优于单层等效源方法.单层等效源方法的计算结果仍然在边界处存在较大误差，且随着计算高度的增加误差也随之增加.
图 7 单层与三层等效源方法计算磁异常梯度与理论值绝对误差a~c与d~f分别为单层与三层等效源方法在地面上计算（ΔTx, ΔTy, ΔTz）与理论值差值；g~i与j~l分别是单层和三层等效源方法在350 m高度平面上对应差值Fig. 7. Distribution of absolute errors between the computational magnetic gradients and theoretical values磁场三分量转换结果对应的计算误差（图 8）说明了相同的问题.而且，相对于磁异常梯度的转换，在三分量转换中，三层等效源方法的计算误差要远低于单层等效源方法的计算误差.
图 8 单层与三层等效源方法计算磁场分量与理论值绝对误差a~c与d~f分别为单层与三层等效源方法在地面上计算（Hax, Hay, Za）与理论值差值；g~i与j~l分别是单层和三层等效源方法在350 m高度平面上对应差值Fig. 8. Distribution of absolute errors between the computational components and theoretical values通过模型实验，能够证明本文所提出的方法有效提高了计算精度，相对于传统单层等效源方法，显著改善了等效源技术在磁场数据融合方面的应用效果.而且，等效源方法相对于常规方法（如傅里叶变换和克里金插值方法），其最大优势在于能够直接处理原始点位上的数据，确保分量及梯度转换结果的同源属性，以及同时实现插值、延拓和多种参量间的转换.
3. 实测数据应用
实测数据来自于湖北金牛火山岩盆地的地面和航空观测ΔT磁异常.由于自然（地形）和人为（房屋、高压线）因素的影响，地面磁测工作极易出现空白区，且测点分布不规则，研究区测点的平均间距为200 m（图 9a）；地形最高点288 m，最低点6 m，东侧地势较高，西侧地势相对较低，高程观测数据分布与磁测点位一致（图 9b）.航空数据的规则性和连续性较地面观测要好，且覆盖范围较大，航测数据观测高度为500 m，点距为200 m.本文利用航空（图 9c）和地面磁测数据的融合计算，对地面磁测数据进行规则网格化、补空、曲化平、扩边和参量转换处理.采用IGRF12模型估算研究区磁化倾角与偏角为（46.42°，-4.4°），不考虑剩磁及退磁影响（李端，2018）.通过磁测数据对数功率谱估算场源深度（图 10），从而建立三层等效源模型（表 3），计算得到过地形最高点平面上的ΔT磁异常（图 11）、磁异常梯度及磁场三分量（图 12）.
图 9 实测地面及航空ΔT磁异常a.地形高程；b.地面实测ΔT磁异常；c.航空测量ΔT磁异常Fig. 9. Measured ΔT anomaly and elevation
图 10 实测地面与航空磁异常对数功率谱曲线图中红色与蓝色直线为对数功率谱曲线分段拟合线；a，b分别为地面与航空磁测数据计算曲线Fig. 10. Radial logarithmic power spectrum of observed magnetic data表 3 实测数据应用等效源层参数Table Supplementary Table Parameters of equivalent source layers of the actual application层号 单元体几何参数
(长×宽×高)(m)等效层顶面深度 展布形态 1 200×200×200 距地表200 m 随地形起伏 2 200×200×600 距0 m平面2 600 m 水平 3 2 000×2 000×6 000 距0 m平面6 000 m 水平 
图 11 计算结果与观测数据拟合差以及288 m高度平面上ΔT计算结果a.恢复地面磁测数据拟合差；b.恢复航空磁测数据拟合差；c.计算ΔTFig. 11. Difference between calculation ΔT and observation values and computational ΔT on the plane at 288 m
图 12 288 m高度平面上磁异常梯度及磁场分量转换结果a~c分别为转换ΔTx，ΔTy与ΔTz；d~f分别为转换Hax，Hay与ZaFig. 12. Transformations on the plane at 288 m计算所获得288 m高度平面上丰富的磁场信息，细化了航空磁测数据，填补了地面数据中存在的空白区并使其分布范围得到扩充，减弱了数据中的地形效应，较好地反映了研究区的地质构造特征，多个高磁异常可能对应该区域的隐伏岩体.磁异常梯度转换结果，为确定场源的空间位置提供了更加细节的信息，有助于提高最终解释成果的准确性.磁场分量转换成功从观测获得的标量数据中提取了地磁场矢量信号，利用矢量数据将极大地提高研究区内地质问题解释推断成果的可靠性.
4. 结论与建议
本文在前人研究工作的基础上，针对地面与航空磁测数据融合问题，根据观测数据中信号的波长特点，提出三层等效源技术方案以期获得更完整的地磁场信息，并通过融合实现地面磁测数据的插值补空、上延、扩边、航空数据下延以及磁异常梯度和磁场分量间的转换.根据理论模型试验和实测数据应用，得到以下4点认识：
（1）鉴于单层等效源方法计算结果中存在较大误差，尤其是分量转换结果误差最为显著.本文认为仅利用一层等效源很难使观测数据中包含的短-中-长波长信息优化或有序的分配到等效源模型中；由于反映到单层等效源模型中的信息极易产生混叠与丢失，在进一步的转换处理中必然会对计算结果产生严重影响.
（2）三层等效源是针对观测数据的频段特点提出的一种等效源模型，其设置参数根据观测数据所反映的场源平均深度、数据分辨率以及数据空间分布范围及形态进行取值，能够保证观测信息在等效源模型中的合理分配，增强了等效源模型的针对性和结构性，减少了等效源设置的盲目性，显著改善了计算精度，充分保证了数据融合的有效性.
（3）根据等效源模型多层的特点，基于预条件共轭梯度算法在反演计算过程中加入深度加权函数，避免等效源物性过度集中于浅层，确保深层等效源能够发挥其应有的作用.通过构建预条件矩阵引入深度加权的方式简洁且有效，不增加计算难度，并可以预处理核矩阵以提高计算稳定性.
（4）相对于磁场数据处理常规方法，等效源方法具有3点优势，其一是可直接对原始点位上的观测数据进行计算，即不要求观测数据分布在水平面或规则网格之上；其二是可进行多种参量间的转换；最后，有助于确保各转换结果的同源性.因此，该方法具有更好的适用性和经济性.
综上，本文提出的三层等效源技术方案较好地实现了地面与航空磁测数据间的融合，提高了磁场数据处理的计算精度.但需要强调的是，由于深部增加了等效源层，计算量也会随之增加；当面对大数据量时，需要较高性能的计算设备和较长的计算时间.因此，在后续的工作中，建议优化反演算法以减少迭代次数，或结合矩阵压缩技术减少计算量，从而提高计算速度.
致谢: 感谢《湖北大冶金牛火山岩盆地1∶2.5万重磁调查项目》、《湖北大冶金牛火山岩盆地重磁异常研究及深部找矿预测》项目以及《地质过程与矿产资源国家重点实验室自主研究课题》对本研究的支持.中国地质大学（武汉）地球物理与空间信息学院梁青博士、刘营博士以及柳一鸣、涂晨明等同学对本研究提供了宝贵的建议，刘双博士在反演工作方面给予了指导性的帮助，在此一并感谢！ -
图 1 三层等效源模型示意图
RLPS表示径向对数功率谱估算；zf, zs, zt分别表示3个等效源层中心埋深
Fig. 1. Three-layer equivalent source model
图 2 理论模型及其地面上观测ΔT磁异常
a.组合模型；b.地形高程（红框为地面观测范围）；c.模拟地面观测ΔT磁异常
Fig. 2. Synthetic models and their theoretical ΔT anomaly on ground
图 3 不同高度正演ΔT磁异常(a~c)和磁异常梯度正演结果(d~i)以及磁场分量正演结果(j~o)
a，b，c分别为地面、350 m和500 m高度平面上ΔT计算结果；d~f为地面上磁异常梯度计算结果；g~i为350 m高度平面上相应计算结果；j~l为地面上磁场三分量计算结果；m~o为350 m高度平面上相应计算结果
Fig. 3. Theoretical ΔT anomaly at different elevations (a-c), theoretical magnetic gradients (d-i) and theoretical components (j-o) on ground and plane at 350 m
图 4 理论模拟地面与航空磁异常对数功率谱曲线
图中红色与蓝色直线为对数功率谱曲线分段拟合线；a，b分别为地面与航空磁测数据计算曲线
Fig. 4. Radial logarithmic power spectrum of observed magnetic data
图 5 等效源单元磁场响应随埋深变化的衰减曲线
a.对应200 m×200 m×600 m等效源单元；b.对应2 000 m×2 000 m×6 000 m等效源单元
Fig. 5. Magnetic attenuation with depths of equivalent source cells
图 6 三种方法计算ΔT磁异常与理论值绝对误差
a~c分别为单层、三层等效源方法与傅里叶‒克里金方法在地面上的计算结果与理论值差值；d~f分别为350 m高度平面上对应差值
Fig. 6. Distribution of absolute errors between the computational ΔT and theoretical values
图 7 单层与三层等效源方法计算磁异常梯度与理论值绝对误差
a~c与d~f分别为单层与三层等效源方法在地面上计算（ΔTx, ΔTy, ΔTz）与理论值差值；g~i与j~l分别是单层和三层等效源方法在350 m高度平面上对应差值
Fig. 7. Distribution of absolute errors between the computational magnetic gradients and theoretical values
图 8 单层与三层等效源方法计算磁场分量与理论值绝对误差
a~c与d~f分别为单层与三层等效源方法在地面上计算（Hax, Hay, Za）与理论值差值；g~i与j~l分别是单层和三层等效源方法在350 m高度平面上对应差值
Fig. 8. Distribution of absolute errors between the computational components and theoretical values
图 9 实测地面及航空ΔT磁异常
a.地形高程；b.地面实测ΔT磁异常；c.航空测量ΔT磁异常
Fig. 9. Measured ΔT anomaly and elevation
图 10 实测地面与航空磁异常对数功率谱曲线
图中红色与蓝色直线为对数功率谱曲线分段拟合线；a，b分别为地面与航空磁测数据计算曲线
Fig. 10. Radial logarithmic power spectrum of observed magnetic data
图 11 计算结果与观测数据拟合差以及288 m高度平面上ΔT计算结果
a.恢复地面磁测数据拟合差；b.恢复航空磁测数据拟合差；c.计算ΔT
Fig. 11. Difference between calculation ΔT and observation values and computational ΔT on the plane at 288 m
图 12 288 m高度平面上磁异常梯度及磁场分量转换结果
a~c分别为转换ΔTx，ΔTy与ΔTz；d~f分别为转换Hax，Hay与Za
Fig. 12. Transformations on the plane at 288 m
表 1 理论模型参数
Table 1. Parameters of synthetic models
模型体 中心埋深
(m)磁化强度
A·m-1）磁化方向 1 1 000 1.5 倾角I=45°
偏角A=5°2 900 2 3 1 600 2.5 4 4 000 4 5 5 500 3 6 6 500 5 
下载: 导出CSV
表 2 理论模型试验等效源层参数
Table 2. Parameters of equivalent source layers of the theoretical model
层号 单元体几何参数
(长×宽×高)(m)等效层顶面深度 展布形态 1 200×200×200 距地表600 m 随地形起伏 2 200×200×600 距0 m平面2 500 m 水平 3 2 000×2 000×6 000 距0 m平面5 500 m 水平
下载: 导出CSV
表 3 实测数据应用等效源层参数
Table 3. Parameters of equivalent source layers of the actual application
层号 单元体几何参数
(长×宽×高)(m)等效层顶面深度 展布形态 1 200×200×200 距地表200 m 随地形起伏 2 200×200×600 距0 m平面2 600 m 水平 3 2 000×2 000×6 000 距0 m平面6 000 m 水平
下载: 导出CSV
-
[1] An, Y.L., Chai, Y.P., Zhang, M.H., et al., 2013. An Optimal Model of the Equivalent Source for Reduction-to-Plane of Potential Field on Uneven Surface and the New Method to Deduce Unit Potential Field Expression of the Optimal Model. Chinese Journal of Geophysics, 56(7): 2473-2483 (in Chinese with English abstract). doi: 10.1002/cjg2.20045/full [2] Andrews, S. B., Moore, P., King, M. A., 2015. Mass Change from GRACE: A Simulated Comparison of Level-1B Analysis Techniques. Geophysical Journal International, 200(1): 503-518. https://doi.org/10.1093/gji/ggu402 [3] Asgharzadeh, M. F., von Frese, R. R. B., Kim, H. R., 2008. Spherical Prism Magnetic Effects by Gauss-Legendre Quadrature Integration. Geophysical Journal International, 173(1): 315-333. https://doi.org/10.1111/j.1365-246X.2007.03692.x [4] Barnes, G., Lumley, J., 2011. Processing Gravity Gradient Data. Geophysics, 76(2): I33-I47. https://doi.org/10.1190/1.3548548 [5] Barzaghi, R., Tselfes, N., Tziavos, I. N., et al., 2008. Geoid and High Resolution Sea Surface Topography Modelling in the Mediterranean from Gravimetry, Altimetry and GOCE Data: Evaluation by Simulation. Journal of Geodesy, 83(8): 751-772. https://doi.org/10.1007/s00190-008-0292-z [6] Bhattacharyya, B.K., Chan, K.C., 1977. Reduction of Magnetic and Gravity Data on an Arbitrary Surface Acquired a Region of High Topographic. Geophysics, 42(42): 1411-1430. https://doi.org/10.1190/1.1440802 http://adsabs.harvard.edu/abs/1977Geop...42.1411B [7] Clark, D.A., 2013. New Method for Interpretation of Magnetic Vector and Gradient Tensor Data Ⅱ: Application to the Mount Leyshon Anomaly, Queensland, Australia. Exploration Geophysics, 44(2): 114-127. https://doi.org/10.1071/EG12066 [8] Cordell, L., Grauch, V.J.S., 1982. Reconciliation of the Discrete and Integral Fourier Transforms. Geophysics, 47(2): 237-243. https://doi.org/10.1190/1.1441330 [9] Dampney, C. N. G., 1969. The Equivalent Source Technique. Geophysics, 34(1): 39-53. doi: 10.1190/1.1439996 [10] Du, J.S., Chen, C., 2015. Progress and Outlook in Global Lithospheric Magnetic Field Modelling by Satellite Magnetic Measurements. Progress in Geophysics, 30(3): 1017-1033 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTOTAL-DQWJ201503005.htm [11] Featherstone, W.E., 2003. Software for Computing Five Existing Types of Deterministically Modified Integration Kernel for Gravimetric Geoid Determination. Computer and Geosciences, 29(2): 183-193. https://doi.org/10.1016/S0098-3004(02)0074-2 doi: 10.1016/S0098-3004(02)00074-2 [12] Gao, X.B., Li, S.S., Li, H., et al., 2013. Application of Point Mass Model and Least Square Collocation in Multi-Source Gravity Data Fusion. Geodesy and Geodynamics, 33(1): 145-149 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTOTAL-DKXB201301034.htm [13] Guo, Z.H., Guan, Z.N., Xiong, S.Q., 2004. Cuboid ΔT and Its Gradient Forward Theoretical Expressions without Analytic Odd Points. Chinese Journal of Geophysics, 47(6): 1131-1138 (in Chinese with English abstract) [14] Kim, H. R., von Frese, R. R. B., Taylor, P. T., et al., 2007. Improved Magnetic Anomalies of the Antarctic Lithosphere from Satellite and Near-Surface Data. Geophysical Journal International, 171(1): 119-126. https://doi.org/10.1111/j.1365-246X.2007.03516.x [15] Li, D., 2018. Reconstruction Method of Gravity and Magnetic Fields by Equivalent Sources (Dissertation). China University of Geosciences, Wuhan (in Chinese with English abstract). [16] Li, D., Chen, C., Liang, Q., et al., 2018. Reconstruction of Discrete Data Using Three-Tier Equivalent Source with Variable Size. Earth Science, 43(3): 873-886 (in Chinese with English abstract). [17] Li, D., Liang, Q., Du, J., et al., 2020. Transforming Total-Field Magnetic Anomalies into Three Components Using Dual-Layer Equivalent Source. Geophysical Research Letters, 47(3): e2019GL084607. https://doi.org/10.1029/2019GL084607 doi: 10.1029/2019GL084607 [18] Li, J.C., Chao, D.B., Ning, J.S., 1995. Spherical Cap Harmonic Expansion for Local Gravity Field Representation. Manuscripta Geodaetica, 20: 265-277. http://www.ingentaconnect.com/content/ssam/03408825/1995/00000020/00000004/art00004 [19] Li, Y., Nabighian, M., Oldenburg, D.W., 2014. Using an Equivalent Source with Positivity for Low-Latitude Reduction to the Pole without Striation. Geophysics, 79(6): J81-J90. https://doi.org/10.1190/GEO2014-0134.1 doi: 10.1190/geo2014-0134.1 [20] Li, Y., Oldenburg, D.W., 1996. 3-D Inversion of Magnetic Data. Geophysics, 61(2): 394-408. https://doi.org/10.1190/1.1443968 [21] MacLennan, C.A., Li, Y.G., 2013. Denoising Multicomponent CSEM Data with Equivalent Source Processing Techniques. Geophysics, 78(3): 125-135. https://doi.org/10.1190/GEO2012-0226.1 http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SEGEAB000029000001000774000001&idtype=cvips&gifs=Yes [22] Maus, S., Barckhausen, U., Berkenbosch, H., et al., 2009. EMAG2: A 2-Arc Min Resolution Earth Magnetic Anomaly Grid Compiled from Satellite, Airborne, and Marine Magnetic Measurements. Geochemistry, Geophysics, Geosystems, 10(8): Q08005. https://doi.org/10.1029/2009GC002471 doi: 10.1029/2009GC002471/abstract [23] Oliveira Jr, V.C., Barbosa, V.C.F., Uieda, L., 2013. Polynomial Equivalent Layer. Geophysics, 78(1): 1-13. https://doi.org/10.1190/geo2012-0196.1 http://adsabs.harvard.edu/abs/2013Geop...78G...1O [24] Ou, J. M., Du, A. M., Thébault, E., et al., 2013. A High Resolution Lithospheric Magnetic Field Model over China. Science China Earth Sciences, 56(10): 1759-1768. https://doi.org/10.1007/s11430-013-4580-y [25] Pang, X.L., 2012. Research on Reduction of Aeromagnetic Anomalies by Means of Equivalent Source Technology (Dissertation). China University of Geosciences, Beijing (in Chinese with English abstract), . [26] Pilkington, M., 1997. 3-D Magnetic Imaging Using Conjugate Gradients. Geophysics, 62(4): 1132-1142. https://doi.org/10.1190/1.1444214 [27] Purucker, M.E., 1990. The Computation of Vector Magnetic Anomalies: A Comparison of Techniques and Errors. Physics of the Earth and Planetary Interiors, 62: 231-245. https://doi.org/10.1016/0031-9201(90)90168-W [28] Purucker, M.E., Ravat, D., Prey, H., et al., 2000. An Altitude-Normalized Magnetic Map of Mars and Its Interpretation. Geophysical Research Letters, 27(16): 2449-2452. doi: 10.1029/2000GL000072 [29] Ravat, D., Langel, R.A., Purucker, M., et al., 1995. Global Vector and Scalar Magsat Magnetic Anomaly Maps. Journal of Geophysical Research, 100: 20111-20136. doi: 10.1029/95JB01237 [30] Silva, J.B.C., Santos, D.F., Garabito, G., 2014. Harmonic and Biharmonic Biases in Potential Field Inversion. Geophysics, 79(1): G15-G25. https://doi.org/10.1190/GEO2013-0137.1 doi: 10.1190/geo2013-0137.1 [31] Sjöberg, L. E., 2005. A Local Least-Squares Modification of Stokes' Formula. Studia Geophysica et Geodaetica, 49(1): 23-30. https://doi.org/10.1007/s11200-005-1623-7 [32] Stolz, R., Zakosarenko, V., Schulz, M., et al., 2006. Magnetic Full-Tensor SQUID Gradiometer System for Geophysical Application. The Leading Edge, 25(2): 178-180. https://doi.org/10.1190/1.2172308 [33] Syberg, F.J.R., 1972. A Fourier Method for the Regional-Residual Problem of Potential Fields. Geophysical Prospecting, (20): 47-75. https://doi.org/10.1111/j.1365-2478.1972.tb00619.x [34] Tikhonov, A.N., Arsenin, V.Y., 1977. Solution of Ill-Posed Problem. Mathematics of Computation, 32(144): 491-491. [35] Whaler, K.A., 1994. Downward Continuation of Magsat Lithospheric Anomalies to the Earth's Surface. Geophysical Journal International, 116: 267-278. https://doi.org/10.1111/j.1365-246X.1994.tb01797.x [36] Wu, Y.H., Luo, Z.C., Zhou, B.Y., 2016. Regional Gravity Modeling Based on Heterogeneous Data Sets by Using Poisson Wavelets Radial Basis Functions. Chinese Journal of Geophysics, 59(3): 852-864 (in Chinese with English abstract). doi: 10.6038/cjg20160308 [37] Xia, J., Sprowl, D.R., 1991. Correction of Topographic Distortion in Gravity Data. Geophysics, 56(4): 537-541. https://doi.org/10.1190/1.1443070 [38] Xie, R.K., Wang, P., Duan, S.L., et al., 2015. Analysis of the Reduction of Aeromagnetic Gradients Data to a Horizontal Plane. Progress in Geophysics, 30(6): 2836-2840 (in Chinese with English abstract). http://search.cnki.net/down/default.aspx?filename=DQWJ201506051&dbcode=CJFD&year=2015&dflag=pdfdown [39] Zhang, W., Zhang, X.J., Tong, J., et al., 2018. Gravity and Magnetic Anomaly Characteristics and Its Geological Interpretation in Rizhao and Lianyungang Areas. Earth Science, 43(12): 4490-4497 (in Chinese with English abstract). http://en.cnki.com.cn/Article_en/CJFDTotal-DQKX201812017.htm [40] Zhou, J., Meng, X., Guo, L., et al., 2015. Three-Dimensional Cross-Gradient Joint Inversion of Gravity and Normalized Magnetic Source Strength Data in the Presence of Remanent Magnetization. Journal of Applied Geophysics, 199: 51-60. https://doi.org/10.1016/j.jappgeo.2015.05.001 http://www.sciencedirect.com/science/article/pii/S0926985115001512 [41] 安玉林, 柴玉普, 张明华, 等, 2013. 曲化平用最佳等效源模型及其单位位场表达式推导的新方法. 地球物理学报, 56(7): 2473-2483. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201307032.htm [42] 杜劲松, 陈超, 2015. 基于卫星磁测数据的全球岩石圈磁场建模进展与展望. 地球物理学进展, 30(3): 1017-1033. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201503005.htm [43] 高新兵, 李珊珊, 李海, 等, 2013. 点质量模型与最小二乘配置在多源重力数据融合中的应用. 大地测量与地球动力学, 33(1): 145-149. https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201301034.htm [44] 郭志宏, 管志宁, 熊盛青, 2004. 长方体ΔT场及其梯度场无解析奇异点理论表达式. 地球物理学报, 47(6): 1131-1138. doi: 10.3321/j.issn:0001-5733.2004.06.029 [45] 李端, 2018. 基于等效源技术的重磁场重构方法(博士学位论文). 武汉: 中国地质大学. [46] 李端, 陈超, 梁青, 等, 2018. 基于三层变尺度等效源的离散重力数据重构. 地球科学, 43(3): 873-886. doi: 10.3799/dqkx.2017.513 [47] 庞旭林, 2012. 航磁异常数据曲面延拓等效源法技术研究(硕士学位论文). 北京: 中国地质大学. [48] 吴怿昊, 罗志才, 周波阳, 2016. 基于泊松小波径向基函数融合多源数据的局部重力场建模. 地球物理学报, 59(3): 852-864. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201603008.htm [49] 谢汝宽, 王平, 段树岭, 等, 2015. 航磁梯度数据曲化平分析. 地球物理学进展, 30(6): 2836-2840. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201506051.htm [50] 张婉, 张玄杰, 佟晶, 等, 2018. 日照-连云港地区重磁异常特征及其构造意义. 地球科学, 43(12): 4490-4497. doi: 10.3799/dqkx.2018.518 -
下载:
点击查看大图
计量
- 文章访问数: 554
- HTML全文浏览量: 219
- PDF下载量: 36
- 被引次数: 0
