2018, Vol. 61
等效源方法是利用虚拟场源模拟实测异常的一种常用位场数据处理方法，可以用于位场数据的空间延拓(包括曲面延拓)、梯度计算以及分量转换等.
根据高斯通量定理，场源的位函数沿包围场源的封闭曲面法向积分与曲面所包含的场源总量有关，而与场源分布无关.依据此原理，理论上在满足位场不变的前提下可以任意地构造“等效”场源，但是所能观测的场分布是有限且离散的，因此合理地构造“等效源”是等效源方法在应用中取得较好效果的一个关键问题.
在等效源方法的应用中，许多学者基于Dampney(1969)的思想将等效源设计成单层模型，例如：在曲面磁测数据处理中，Bhattacharyya和Chan(1977)应用磁偶层等效源进行磁异常转换；Silva(1986)在观测面之下构建单层离散磁偶极子进行磁异常化极；侯重初等(1985)、Hansen和Miyazaki(1984)、王万银等(王万银等, 1991; 王万银和潘作枢, 1996)基于磁偶层等效源实现曲线和曲面上的磁异常转换；安玉林等(2013)提出了“模型变换法”极大地简化了磁偶层正演公式，为该方法的实际应用提供很好的借鉴和参考；Li和Oldenburg(2000)在有限的深度范围设置单层连续直立棱柱体等效源实现低纬度磁异常化极；徐世浙(2004)将单层等效源布置于地表实现航磁数据向下延拓.而关于单层等效源深度的选择，大多采用Dampney(1969)以及Xia和Sprowl(1991)提出的设置准则.
选择单层等效源且将其布置于近地表是目前等效源方法的主要特点.然而，由于观测数据的局限性，近地表的单层等效源模型难以准确地模拟实测数据，主要的问题是位场信号中的长波长成分不易被保全.对此，国内外学者提出了多层等效源模型的方法，如Needham(1970)、吴晓平(1984)等利用多层点质量模型构建等效源保证重力场向上延拓的精度，于德武(2004)采用两层水平圆柱体模拟无限下延的二维直立板体并且进行磁异常换算，庞旭林(2012)利用多层磁偶极子等效源模拟不同深度场源以实现航磁数据的高度调平(或“曲化平”)和向上延拓.
为了保证高精度、高分辨率地重构位场，通常需要设置密集的等效源单元体，因而需要求解大型线性方程组.因此，计算效率及稳定性也是等效源方法实际应用的另一个关键问题.Leão和Silva(1989)通过滑动窗口分块反演将大型核矩阵分解成稳定性较强的小维度矩阵，Mendonça和Silva(1994)对数据进行稀疏压缩，Li和Oldenburg(2010)以及刘天佑等(2007)基于小波压缩技术降低核矩阵维数，Oliveira等(2013)利用多项式描述窗口内等效源分布以减少反演参数等.在稳定性和噪声压制方面，Li和Oldenburg(2001)、Barnes和Lumley(2011)则利用了正则化的思想.
前人的研究为合理地构建等效源模型以及实现快速而稳定的计算提供了很好借鉴.然而，前述的等效源方法依然存在不足，多数等效源方法对应用效果的评价只注重与实测数据的拟合程度，而忽略了等效源分布在物理上的合理性(杜劲松和陈超, 2015).尤其当深部场源产生的位场信号波长大于观测数据分布的空间范围时，在浅部设置的单层等效源往往会因观测数据分布空间有限而无法完整模拟深部场源信号，以致利用等效源重构的位场在向上延拓时会产生畸变.物理上，等效源分布于整个地下空间是最合理的，有可能得到接近真实场源的等效源模型，但是如此无异于地下空间三维物性反演，计算效率和难度可想而知.而多层等效源方法事实上是单层等效源与地下全空间等效源之间的一种折中方法，是一个相对合理的选择，既能保证等效源模型充分表达磁异常不同频段成分，又不失方法简洁与高效的优势.然而目前的多层等效源方法在布设方案中缺少等效源深度设置方面的讨论.此外，等效源单元体形状的选择也会影响方法的应用效果.本文针对以上问题进行了相应研究，提出一种基于多层等效源的技术方案应用于磁异常数据的处理与转换，在保证计算效率的同时，可有效增强模型的合理性进而改善等效源方法的应用效果.
1 方法原理 1.1 基本理论与方程根据位场原理，位场函数F(x, y, z)与场源的分布函数m(ξ, η, ζ)在场源以外空间Ω的关系可用积分方程表示，即
|
(1) |
式中A为积分核函数，是表征场与场源之间物理几何关系的函数.若位场数据与场源是离散的，且每个离散单元内物性是均匀的，则公式(1)可写成
|
(2) |
其中
|
(3) |
式中Γj为第j个离散场源单元的积分区域.假设离散位场数据有M个，用di表示第i个观测点的位场数据，则公式(2)可以简化成
|
(4) |
将公式(4)写成矩阵形式，即
|
(5) |
式中d为M阶向量，m为N阶向量，G为M×N阶矩阵.
1.2 等效源单元体选择等效源单元体的选择对方法的应用效果具有一定影响.离散的等效源单元可以采用不同的形式，包括磁偶极子、磁偶层及均匀磁化规则形体等.关于等效源单元体的选择，如谢汝宽等(2015)认为单元体在深度方向延伸过小将影响深部信息的获取，以致于向上延拓时会产生较大误差.磁偶极子、磁偶面元和长方体的对比试验(见附录A)表明，采用长方体有助于“吸收”更多的深部场源信息，且随着等效源深度增加，可增大等效源单元体尺寸，以保持应有的敏感度并减少等效源单元体数量.因此本文采用长方体作为等效源单元体构建等效源，其正演计算采用郭志宏等(2004)给出的解析式.
假设等效源单元体受到均匀磁化，则公式(4)可以写成
|
(6) |
式中，κj为第j个等效源单元体磁化强度，di表示第i个计算点处等效源单元体的磁总场异常强度.
因此，可以将等效源单元体的物性参数作为未知量，其空间位置与几何形态信息作为已知量给定，并假设磁测数据总数为M，等效单元总数为N，根据位场叠加原理建立等效源物性与磁异常之间的线性关系，利用矩阵形式表示如下
|
(7) |
式中，dM=(d1, d2, …, dM)表示磁异常数据向量，κN=(κ1, κ2, …, κN)表示等效源物性向量，GM×N表示正演核函数矩阵.
1.3 等效源设置如前所述，利用单层等效源模拟真实场源会导致信号“泄漏”.模型试验表明，“泄漏”的信号分布在不同的频段(见附录B)，在地下深部设置更多的等效源层是必要的.显然，多层等效源的层数及其深度设置是方法的关键，层数越多、计算成本越高(见附录C).此外，深部等效源层需要适当向测区外扩展，尽可能拟合测区内长波长(或低频)信号，以减小边界效应和深部信号采样不足造成的失真.通过理论模型试验，并借助等效源单元体尺寸可变性和场源深度估计技术，本文提出的等效源设置方案是利用三层等效源的合理配置模拟不同深度场源.具体思路如下：
(1) 第一层等效源设置
第一层等效源用于拟合磁异常的主体部分，尤其是高频成分.等效源单元体尺寸可依据观测点平均间距而定(MacLennan and Li, 2013).该等效源层可置于观测面之下的水平面上，也可置于起伏观测面之下的等深曲面上.综合前人的建议(Dampney, 1969; Xia and Sprowl, 1991)，该等效源层深度设置在3~6倍于磁测数据平均点距为宜，这样既能较好地“吸收”观测数据中的高频信息，又不易产生虚假噪声.
(2) 第二层等效源设置
第二层等效源设置的目的是模拟磁测异常的中-低频成分，弥补第一层等效源对场源信息吸收的不足.根据模型试验，第二层等效源单元体几何尺寸可选择第一层的1.5~3倍，横向展布方式可参考第一层等效源，并将该层等效源设置范围横向拓展至略大于实测数据范围.第二层等效源的深度是重要参数之一，需要对场源深度进行分析.可利用磁异常数据的谱分析(如Ravat et al., 2007; Li, 2013)获得深部场源的等效平均深度；也可利用欧拉反褶积的方法(如Thompson, 1982; Reid and Thurston, 2014)进行估计.当观测面起伏较大时，可能会对功率谱估算造成影响，但将其作为该层等效源深度的参考仍是可行的.
(3) 第三层等效源设置
理论模型试验表明，当研究区域磁异常存在一个低频背景场时，上述两层等效源也很难得到好的拟合效果.这种低频背景场在实际资料中经常存在，如深部磁性界面或是居里等温面深度变化引起的磁异常.因此，第三层等效源的设置旨在拟合区域低频背景场，同时减弱边界效应.该等效源层单元体可选用厚度较大的长方体或磁荷面，尺寸为第二层的2~5倍，设置在第二层等效源之下，呈水平展布，且分布范围为2~5倍的数据范围.等效源层的深度可依据当地居里等温面或磁性界面深度确定.在没有先验信息的情况下，根据理论模型试验，该等效源层深度可以选择3~10倍的第二层等效源(平均)深度.
1.4 等效源物性参数反演多层等效源的反演方式可以分为两类，一类是由深至浅，分步、逐级反演，这种方案在一定程度上使计算过程复杂化，不易控制拟合误差；且易导致深层等效源人为分配更多的场源信息，而浅部等效源得不到有效分配；另一类是多层等效源同时反演，可使各层等效源依据异常场信息自动分配.本文采用三层等效源同时反演的策略.
根据本文等效源的设置规则，(7)式中观测数据数量一般会少于等效源单元体数量(M < N)，即为欠定方程，因此利用正则化方法(Tikhonov and Arsenin, 1977)建立反演目标函数：
|
(8) |
式中Wd为加权对角矩阵，对角元素为1/σj，其中σj为第j个数据的误差标准差，旨在平衡不同测点对数据拟合差的贡献.考虑到反演问题是不稳定的，加入正则化项能够有效提高解的稳定性，λ为正则化参数，根据观测数据受干扰程度确定λ取值，随干扰信号由弱变强，λ取值相应由小增大.Wm为对角矩阵，其对角元素为
|
(9) |
进而利用预优共轭梯度法求解.
预优共轭梯度法是求解线性方程的有效方法之一(Pilkington, 1997; 刘双等, 2013)，通过预处理使核矩阵的特征值比较集中分布，从而提高收敛速度.具体方法是在方程(9)两边同乘以预优矩阵P，使P(GTWdTWdG+λWmTWm)≈I，本文借鉴Pilkington(1997)的方法，取预优矩阵为
|
(10) |
式中，z表示单元体中心到观测面的垂直距离(或平均垂直距离)，β为衰减系数，与场的衰减速率有关, 可以参考Li和Oldenburg(1996)的“深度加权因子”.本文根据模型试验，对磁总场强度数据建议取值3~6效果较好.加入预优矩阵，方程(9)变成
|
(11) |
这里预优因子zβ在改善矩阵特征值的同时，还起到了深度加权的作用，使反演得到的等效源物性能够更好地分布在不同深度的等效层上，得到更为合理的等效源模型.
1.5 磁异常数据转换将反演得到的等效物性κinv带入相关待重构或转换磁异常参量的正演公式，可得
|
(12) |
式中，drec表示重构或转换磁异常参量向量，可以是磁总场强度异常、或其不同方向上的梯度，或磁异常矢量的分量，G*代表计算参量相应的正演核矩阵.
2 理论模型试验针对前面提出的问题，此处设计一个二维组合模型，分别采用单层等效源模型和三层等效源模型进行模拟与对比试验.
二维组合模型及其理论磁异常如图 1所示，其中模型参数见表 1，模型包括不同尺寸的长方水平柱体和深部磁性场源(图 1b)，分布于起伏地形(高程变化0~180 m)之下；地磁场感应强度为50000 nT，地磁倾角I=40°，剖面取磁北方向；并将深部场源分布范围设计为观测剖面长度的2.5倍.
|
图 1 二维模型及其磁异常 (a)组合模型理论磁异常；(b)组合模型及三层等效源设置示意图. Fig. 1 2D synthetic models and theoretical magnetic anomaly (a) Theoretical magnetic anomaly of models; (b) 2D models and equivalent source (ES) layers. |
|
表 1 二维模型参数 Table 1 Parameters of 2D models |
本试验设计了三层等效源模型，其中第一层等效源模型同时被用作单层等效源方法进行测试.考虑到理论磁异常(图 1a)剖面的离散点距取10 m，第一层等效源层顶距地表深度为60 m；根据异常频谱分析结果(图 2)设置第二层等效源层顶深度为400 m；第三层等效源层顶距地面平均深度约为1400 m.三层等效源参数设置见表 2.
|
图 2 二维模型理论磁异常对数功率谱 Fig. 2 The logarithmic power spectrum of 2D theoretical magnetic anomaly |
|
|
表 2 等效源层参数 Table 2 Parameters of equivalent source layers |
分别采用单层等效源方法和三层等效源方法对二维模型(图 1b)进行了试验.模拟计算了剖面地表、高程为180 m平面(曲化平)以及高程为300 m平面(向上延拓)上的磁异常ΔT(图 3)，同时也计算了ΔT水平和垂直方向导数即ΔTx和ΔTz(图 4)以及水平磁异常Hax和垂直磁异常Za(图 5)，并将它们与模型理论值进行对比.
|
图 3 单层与三层等效源计算的ΔT磁异常与模型理论异常比较 (a), (b), (c)分别为起伏观测面和不同高度(h=180 m和h=300 m)平面上ΔT计算值与模型理论值；(d), (e)和(f)分别为计算值与模型理论值之差值. Fig. 3 Comparison of magnetic anomaly ΔT calculated by ES with the model′s theoretical anomaly (a), (b), (c) show the theoretical ΔT and those calculated by single layer and three layers ES on observation surface, planes at h=180 m and h=300 m respectively. (d), (e), (f) show their differences respectively. |
|
图 4 单层与三层等效源计算的ΔT梯度与模型理论异常梯度比较 (a), (b), (c)和(g), (h), (i)分别为起伏观测面和不同高度(h=180 m和h=300 m)平面上ΔTx和ΔTz计算值与模型理论值；(d), (e), (f)和(j), (k), (l)分别为ΔTx和ΔTz计算值与理论值之差值. Fig. 4 Comparison of ΔT derivatives calculated by ES with the model′s theoretical derivatives (a), (b), (c) and (g), (h), (i) show the theoretical ΔTx and ΔTz and those calculated by single layer and three layers ES on observation surface, planes at h=180 m and h=300 m, respectively. (d), (e), (f) and (j), (k), (l) show their differences respectively. |
|
图 5 单层与三层等效源计算磁异常分量与模型理论值比较 (a), (b), (c)和(d), (e), (f)分别为起伏观测面和不同高度(h=180 m和h=300 m)平面上Hax和Za计算值与模型理论值. Fig. 5 Comparison of components of magnetic anomaly calculated by ES with the model′s theoretical components (a), (b), (c) and (d), (e), (f) show the theoretical Hax and Za and results calculated by single layer and three layers ES on observation surface, planes at h=180 m and h=300 m respectively. |
从试验结果来看，在起伏观测面上，单层等效源方法和三层等效源方法均能很好地拟合理论异常(见图 3a和图 3d)，但是随着距观测面高度的增加，单层等效源拟合误差显著增大(见图 3e和图 3f)，尤其在剖面两端，而三层等效源模型的计算值较好地拟合了理论异常.
磁异常水平梯度和垂直梯度是反映异常变化细节特征的参量.因此，本试验引入了磁异常水平梯度和垂直梯度来评价两种方法的效果.一般而言，依据起伏地形上的实测离散数据直接计算磁异常梯度会产生较大的误差，尤其在地形起伏大且磁异常变化剧烈的地带.等效源方法模拟实测异常也会产生这样的误差.图 4显示，单层等效源方法和三层等效源方法获得的地面磁异常梯度与理论值均存在误差，但三层等效源方法的误差明显较小，根据本试验结果统计，起伏地面上单层等效源模型计算结果的均方差约为三层等效源模型的10倍，尽管上延后这种误差会显著减小，但是单层等效源模型计算结果仍存在区域性的差异.可见多层等效源方法可有效提高磁异常导数换算的精度.
当实测数据未能包含地下磁性体异常场完整信息时，利用实测磁异常进行异常分量转换存在着多解性.本试验中加入的背景磁异常有可能造成异常分量转换误差.图 5给出的结果表明，本试验中无论在地表还是上延不同高度，三层等效源模型转换异常分量与理论值均十分接近，而单层等效源模型则因多解性影响存在显著的系统误差.
为进一步验证多层等效源方法的计算效果，本文设计了起伏地表情况下的三维理论模型.三维模型由一组不同尺度和深度的均匀磁化模型体组成(图 6b)，模型参数见表 3.地磁场感应总强度为50000 nT，地磁倾角I=45°，地磁偏角D=10°.起伏地表最高点100 m，最低点0 m，如图 6a所示.观测范围1 km×1 km，测点间距10 m.在模型异常中加入了由起伏磁荷面正演的理论区域异常，模拟实际资料中的低频背景场信息.
|
图 6 三维模型及其观测面上理论总磁异常 (a)起伏地形；(b)三维模型组合；(c)模型磁异常；(d)背景场；(e)总磁异常. Fig. 6 3D synthetic models and their theoretical magnetic anomaly on observation surface (a) Uneven topography; (b) 3D models; (c) Models′ magnetic anomaly; (d) Background magnetic field; (e) Total magnetic anomaly. |
|
|
表 3 三维模型参数 Table 3 Parameters of 3D models |
等效源模型与二维模型试验中的类似，第一层等效源顶面深设置为距观测面30 m；根据径向对数功率谱估计(图 7)，第二层等效源顶面埋深距观测面110 m；第三层等效源单元选用磁荷面，其深度设置在海平面(h=0 m)之下800 m，且分布范围3倍于数据区域；三层等效源参数设置见表 4.
|
图 7 三维模型理论总磁异常径向对数功率谱 Fig. 7 The radial logarithmic power spectrum of total theoretical magnetic anomaly |
|
|
表 4 等效源层参数 |
基于三层等效源模型，可以得到地面、h=100 m(曲化平)与h=200 m(向上延拓)平面上磁异常(ΔT)、磁异常梯度(ΔTx, ΔTy, ΔTz)及分量(Hax, Hay, Za)的模拟结果.简洁起见，本文仅以磁异常(ΔT)、磁异常垂向梯度(ΔTz)及垂向分量(Za)计算结果为例进行分析.
从图 8中可以看出，起伏地面上的磁异常值得到了很好的重构，且曲面延拓结果(h=100 m与h=200 m平面上)也具有很高的拟合精度.通过与理论值对比统计发现，在各结果中绝大部分区域相对误差均在1%以内，与理论值均方差未超过0.09 nT；转换的磁异常垂向梯度(ΔTz)结果与理论值形态特征吻合较好(图 9)，绝大部分区域的相对误差控制在1%左右.在起伏地面、h=100 m与h=200 m平面上，均方差统计分别为0.07 nT/m、0.001 nT/m和0.001 nT/m；磁异常及其梯度的完好拟合说明本文方法能较好地“还原”原始磁异常.从磁异常的垂向分量转换来看(图 10)，在上述三种情况下，等效源计算的磁异常垂向分量与理论值均方差均未超过3 nT，说明本文方法在背景场干扰下也能有效压制多解性.
|
图 8 三层等效源计算ΔT磁异常与理论值的比较 (a), (b), (c)与(d), (e), (f)分别为起伏观测面和不同高度(h=100 m和h=200 m)平面上ΔT模型理论值与计算值；(g), (h), (i)分别为模型理论值与计算值之差值. Fig. 8 Comparison of ΔT calculated by three layers ES with the model′s theoretical anomaly (a), (b), (c) and (d), (e), (f) show the theoretical ΔT and those calculated by three layers ES on observation surface, planes at h=100 m and h=200 m respectively. (g), (h), (i) show their differences respectively. |
|
图 9 三层等效源计算垂向梯度ΔTz与理论值的比较 (a), (b), (c)与(d), (e), (f)分别为起伏观测面和不同高度(h=100 m和h=200 m)平面上ΔTz模型理论值与计算值；(g), (h), (i)分别为模型理论值与计算值之差值. Fig. 9 Comparison of ΔTz calculated by three layers ES with the model′s theoretical ΔTz (a), (b), (c) and (d), (e), (f) show the theoretical ΔTz and those calculated by three layers ES on observation surface, planes at h=100 m and h=200 m respectively. (g), (h), (i) show their differences respectively. |
|
图 10 三层等效源计算垂向分量Za结果与理论值的比较 (a), (b), (c)与(d), (e), (f)分别为起伏观测面和不同高度(h=100 m和h=200 m)平面上Za模型理论值与计算值；(g), (h), (i)分别为模型理论值与计算值之差值. Fig. 10 Comparison of Za calculated by three layers ES with the model′s theoretical Za. (a), (b), (c) and (d), (e), (f) show the theoretical Za and those calculated by three layers ES on observation surface, planes at h=100 m and h=200 m respectively. (g), (h), (i) show their differences respectively. |
实例所用数据来自广西某地区地面实测ΔT磁异常.研究区地形复杂，呈北高南低的趋势，地形最高点776 m，最低点17 m，山谷沟壑遍布整个测区(见图 11a).观测数据的平均测点间距200 m.由图 11b可见，研究区内磁异常信息较为丰富，其中高幅值磁异常主要分布在西部和南部区域，区内磁异常主要反映了隐伏岩体分布及磁性基底起伏.采用三层等效源方法对实测数据进行处理，根据该区已有的地质与地球物理工作成果，三层等效源模型由浅至深设置为：第一层、二层等效源均随地形起伏展布，采用长方体作为等效单元.第一层等效源顶面深距观测面3倍平均点距(600 m)，单元体水平尺寸200 m×200 m，厚度200 m，第二层等效源顶面深度距观测面1200 m(图 11e)，等效单元水平尺寸400 m×400 m，厚度400 m，第三层等效源采用磁荷面，呈水平展布，深度取海平面以下9000 m，单元体水平尺寸1.2 km×1.2 km，分布面积约为研究区面积2.5倍.
|
图 11 研究区地形、实测磁异常及由等效源计算的磁异常 (a)研究区地形；(b)地表实测磁异常；(c)等效源计算的地面磁异常；(d)上延至h=900 m平面上的磁异常；(e)实测磁异常径向对数功率谱. Fig. 11 Topography, observed magnetic anomaly and reconstructive results in study area (a) Topography (observation surface); (b) Observed magnetic anomaly; (c) Reconstructive magnetic anomaly on observation surface; (d) Reconstructive magnetic anomaly on the plane at the level of 900 m; (e) The radial logarithmic power spectrum of observed data. |
根据拟合差统计，绝大部分区域的相对误差未超过3%，符合对该地区干扰信号的估计.因此，利用反演得到的三层等效源模型所计算的磁异常(图 11c)及其上延结果(图 11d)能够充分反映实测异常的区域与局部特征.上延至h=900 m平面上的磁异常，较好地反映了近南北向分布的深部地质构造特征，多个高磁异常区与已确定的岩体也有较好的对应(冯建辉，2016).磁异常梯度计算结果(图 12)很好地突出了场源细节信息，并有效控制了误差和干扰信号的放大；分量转换结果(图 13)提供了磁异常空间矢量信息，可为进一步研究区内场源空间特征及性质提供可靠的基础数据.
|
图 12 地面上磁异常梯度及分量转换结果 Fig. 12 Calculated results of derivatives and components on observation surface |
|
图 13 高程为900 m平面上磁异常梯度及分量转换结果 Fig. 13 Calculated results of derivatives and components on the plane at the elevation of 900 m |
本文在前人的研究基础上，尝试通过设计三层等效源获取更完整的场源信息，并应用于磁异常数据的处理与转换.通过模型试验和实际应用，得到了以下四点认识：
(1) 当场源分布在不同深度范围内，尤其当观测数据中含有背景场信息时，单层等效源方法无法构建出合理的等效源模型；布设浅、中、深三个不同深度的等效源层，可以更完整地获得磁异常信息，进而可以得到更高精度的磁异常曲化平及上延结果.
(2) 多层等效源模型相对单层等效源模型，其计算量与计算难度有所增加，但是本文通过等效源单元参数的合理设置，结合场源深度估计，可减少计算量和等效源深度设置的盲目性；预优共轭梯度方法的应用，可解决反演过程中的“趋肤效应”问题，并提高了计算效率.
(3) 对等效源方法应用效果的评价，传统的方法采用等效源模拟磁异常与实测数据拟合程度来评判.根据本文模型试验，仅凭磁异常拟合差检验是不够的，还应考虑磁异常梯度甚至磁异常转换结果进行联合评判.
(4) 由于三层等效源模型能够更完整地获取场源信息，因此利用所获得的等效源模型进行磁异常参量转换(如分量换算、化极等)更为可靠.
综上，本文提出的三层等效源技术方案能有效增强等效源模型模拟实测异常完整信息的能力，提高磁异常转换的计算精度，可以为后续磁异常数据的处理和解释奠定良好的基础.
致谢中国地质大学(武汉)地球物理与空间信息学院刘双副教授在反演工作中提供了指导和帮助，华中科技大学引力中心孙石达博士对本研究提供了宝贵的意见，在此一并感谢.
附录A为研究三种等效单元的磁异常随高度增加的衰减情况，以h=0 m平面为基准面；在基准面之下设置三种等效源单元，具有相同的磁距及中心埋深(h=-150 m)，磁化倾角和偏角均与地磁倾角I=45°和偏角D=10°相同；计算三种等效源单元在基准面上的磁异常作为基准值，并随高度每增加100 m计算一个磁异常值Fi，最后计算相对于基准值的磁异常衰减率Ri，即
|
从衰减率曲线上可以明显看出(图 A1)，长方体具有更缓慢的衰减特性.
|
图 A1 (a) 三种等效源单元体的磁异常随高度衰减；(b)三种等效源单元体示意图 图中dipole为磁偶极子、plate-element为磁偶面元、rectangular prism为均匀磁化长方体. Fig. A1 (a) Magnetic anomaly attenuation characteristics of three equivalent sources; (b) Three types of equivalent sources The three equivalent sources are magnetic dipole (Dipole), magnetic plate element (Plate-element) and magnetic rectangular prism (Rectangular prism) |
利用正文中的二维模型，在300 m平面上，对比模型理论磁异常与单层等效源和三层等效源计算磁异常的对数功率谱及相位谱，通过误差统计，可以明显地反映出单层等效源的计算结果在各频段上均存在信号“泄露”.
|
图 B1 对数功率谱与相位谱分析 (a1), (a2)为磁异常理论值、单层与三层等效源计算值的对数功率谱及其拟合差；(b1), (b2)为相应的相位谱结果及拟合差. Fig. B1 (a1) The logarithmic power spectrums of theoretical ΔT and those calculated by ES for single layer and three layers; (b1) The results of phase spectrums; (a2), (b2) The differences between theoretical and calculated values respectively |
设计二维组合模型如图 C1所示，观测剖面水平，剖面长度7500 m，高程h=0 m，地磁倾角I=60°，剖面取磁北方向.
|
图 C1 2D磁化模型空间位置及其在观测面上的磁异常响应 Fig. C1 2D synthetic models and theoretical magnetic anomaly |
利用观测磁异常分别建立单层、双层、三层、四层及五层等效源模型(见图 C2)，其中双层等效源模型又根据深层等效源埋深不同分为两种(如图 C2(b、c)所示).另需说明，图 C2中显示的等效源位置为其各层顶面埋深位置.
|
图 C2 等效源模型示意图(各层顶面位置) (a)单层；(b)双层_1; (c)双层_2; (d)三层; (e)四层; (f)五层. Fig. C2 Equivalent source models for different number of layer (showing the top surface in the figure) (a) Single layer; (b) Double layers_1; (c) Double layers_2; (d) Three layers; (e) Four layers; (f) Five layers |
对比不同等效源模型在观测面及750 m高度平面上的磁异常、磁异常梯度及分量(以垂向梯度及分量为例)的计算结果，并统计与理论值的拟合差.
对比分析图 C3及图 C4中的计算结果，单层及双层等效源模型的计算结果中存在较大的误差，当等效源层数增加至三层及以上时，计算精度得到了显著提高.但是在三层的基础上继续增加等效源层数，计算精度并没有明显的提升，即三、四、五层等效源模型的计算精度基本处于同一水平，却增加了计算难度以及降低了计算速度.所以，综合考虑计算效率及等效源模型的合理性，文中提出了三层等效源模型的设置方案，确保在浅、中、深三个深度范围内均有等效源分布，以应对实际场源分布情况，而又不过多增加方法的计算成本.
|
图 C3 观测面上及750 m高度上计算结果及其与理论值的拟合差分布 图中左列为各等效源模型在观测面上的计算磁异常、磁异常垂向梯度和垂向分量及其与理论值的拟合差；右列为各等效源模型在750 m高度上对应计算结果及其与理论值的拟合差. Fig. C3 Calculated results and their difference relative to theoretical values on observation surface and the altitude of 750 m Calculated magnetic anomaly, gradient, comparison and their difference relative to theoretical values on observation surface are listed in the left column, the corresponding results on the altitude of 750 m are listed in the right column |
|
图 C4 对数功率谱分析 (a)磁异常理论值、单层与双层等效源计算值的对数功率谱; (b)磁异常理论值、三层、四层与五层等效源计算值的对数功率谱; (c)计算值与理论值拟合差统计. Fig. C4 Logarithmic power spectrum analysis (a), (b) Logarithmic power spectrums of theoretical ΔT and those calculated by ES for single layer, double, three, four and five layers; (c) show the differences between theoretical and calculated values |
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 (in Chinese), 56(7): 2473-2483. DOI:10.6038/cjg20130733 |
Barnes G, Lumley J. 2011. Processing gravity gradient data. Geophysics, 76(2): 133-147. |
Dampney C N G. 1969. The Equivalent Source Technique. Geophysics, 34(1): 39-53. DOI:10.1190/1.1439996 |
Du J S, Chen C. 2015. Progress and outlook in global lithospheric magnetic field modelling by satellite magnetic measurements. Progress in Geophysics (in Chinese), 30(3): 1017-1033. DOI:10.6038/pg20150305 |
Feng J F. 2016. Extraction and explanation of gravity and magnetic anomaly features in Shedong Area of Guangxi[Master's thesis] (in Chinese). Wuhan: China University of Geosciences (Wuhan).
|
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 (in Chinese), 47(6): 1131-1138. |
Hansen R O, Miyazaki Y. 1984. Continuation of potential fields between arbitrary surfaces. Geophysics, 49(6): 787-795. DOI:10.1190/1.1441707 |
Hou Z C, Cai Z X, Liu K F. 1985. Interpretation system of potential field transformation on an uneven surface from the potential of double layer. Acta Geophysica Sinica (in Chinese), 28(4): 410-418. |
Leão J B C, Silva J B C. 1989. Discrete linear transformations of potential field data. Geophysics, 54(4): 479-507. |
Li C F, Wang J, Lin J, et al. 2013. Thermal evolution of the North Atlantic lithosphere:New constraints frommagnetic anomaly inversion with a fractal magnetization model. Geochemistry, Geophysics, Geosystems, 14(12): 5078-5105. DOI:10.1002/ggge.v14.12 |
Li Y G, Oldenburg D W. 1996. 3-D inversion of magnetic data. Geophysics, 61(2): 394-408. DOI:10.1190/1.1443968 |
Li Y G, Oldenburg D W. 2000. Reduction to the pole using equivalent sources. //2000 SEG Annual Meeting, 6-11 August 2000, Calgary, Alberta. Calgary: Society of Exploration Geophysicists.
|
Li Y G, Oldenburg D W. 2001. Stable reduction to the pole at the magnetic equator. Geophysics, 66(2): 571-578. DOI:10.1190/1.1444948 |
Li Y G, Oldenburg D W. 2010. Rapid construction of equivalent sources using wavelets. Geophysics, 75(3): L51-L59. DOI:10.1190/1.3378764 |
Liu S, Feng J, Gao W L, et al. 2013. 2D inversion for borehole magnetic data in the presence of significant remanence and demagnetization. Chinese Journal of Geophysics (in Chinese), 56(12): 4297-4309. DOI:10.6038/cjg20131232 |
Liu T Y, Yang Y S, Li Y Y, et al. 2007. The order_depression solution for large_scale integral equation and its application in the reduction of gravity data to a horizontal plane. Chinese Journal of Geophysics (in Chinese), 50(1): 290-296. |
MacLennan C A, Li Y G. 2013. Denoising multicomponent CSEM data with equivalent source processing techniques. Geophysics, 78(3): E125-E135. DOI:10.1190/geo2012-0226.1 |
Mendonça C A, Silva J B C. 1994. The equivalent data concept applied to the interpolation of potential field data. Geophysics, 59(5): 722-732. DOI:10.1190/1.1443630 |
Needham P E. 1970. The formation and evaluation of detailed geopotential models based on point masses. Department of Geodetic Science and Surveying, Ohio State University.
|
Oliveira V C Jr, Barbosa V C F, Uieda L. 2013. Polynomial equivalent layer. Geophysics, 78(1): G1-G13. DOI:10.1190/geo2012-0196.1 |
Pang X L. 2012. Research on reduction of aeromagnetic anomalies by means of equivalent source technology[Master's thesis] (in Chinese). Beijing: China University of Geosciences (Beijing).
|
Pilkington M. 1997. 3-D magnetic imaging using conjugate gradients. Geophysics, 62(4): 1132-1142. DOI:10.1190/1.1444214 |
Ravat D, Pignatelli A, Nicolosi I, et al. 2007. A study of spectral methods of estimating the depth to the bottom of magnetic sources from near-surface magnetic anomaly data. Geophysical Journal International, 169(2): 421-434. DOI:10.1111/gji.2007.169.issue-2 |
Reid A B, Thurston J B. 2014. The structural index in gravity and magnetic interpretation:Errors, uses, and abuses. Geophysics, 79(4): J61-J66. DOI:10.1190/geo2013-0235.1 |
Silva J B C. 1986. Reduction to the pole as an inverse problem and its application to low-latitude anomalies. Geophysics, 51(2): 369-382. DOI:10.1190/1.1442096 |
Thompson D T. 1982. EULDPH:A new technique for making computer-assisted depth estimates from magnetic data. Geophysics, 47(1): 31-37. DOI:10.1190/1.1441278 |
Tikhonov A N, Arsenin V Y. 1977. Solution of ill-posed problems. Mathematics of Computation, 32(144): 491-491. |
Wang W Y, Pan Z S, Li J K. 1991. Continuation methods for curved surface of the three-dimensional high-precision gravity and magnetic potential field. Geophysical & Geochemical Exploration (in Chinese), 15(6): 415-422. |
Wang W Y, Pan Z S. 1996. Data processing and transformation methods of potential fields on an arbitray surface with dipole layer potential technique. Journal of Xi'an Engineering University (in Chinese), 18(3): 69-76. |
Wu X P. 1984. Point-mass model of local gravity field. Acta Geodaetica et Cartographic Sinica (in Chinese), 13(4): 250-258. |
Xia J H, Sprowl D R. 1991. Correction of topographic distortion in gravity data. Geophysics, 56(4): 537-541. DOI:10.1190/1.1443070 |
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 (in Chinese), 30(6): 2836-2840. DOI:10.6038/pg20150650 |
Xu S Z, Shen X H, Zou L J, et al. 2004. Downward continuation of aeromagnetic anomaly from flying altitude to terrain. Chinese Journal of Geophysics (in Chinese), 47(6): 1127-1130. |
Yu D W. 2004. Conversions of magnetic anomaly with the equivalent sourse method. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 26(2): 133-135. |
安玉林, 柴玉璞, 张明华, 等. 2013. 曲化平用最佳等效源模型及其单位位场表达式推导的新方法. 地球物理学报, 56(7): 2473-2483. DOI:10.6038/cjg20130733 |
杜劲松, 陈超. 2015. 基于卫星磁测数据的全球岩石圈磁场建模进展与展望. 地球物理学进展, 30(3): 1017-1033. DOI:10.6038/pg20150305 |
冯建辉. 2016. 广西社垌地区重磁异常特征提取与解释[硕士论文]. 武汉: 中国地质大学(武汉).
|
郭志宏, 管志宁, 熊盛青. 2004. 长方体裣T场及其梯度场无解析奇点理论表达式. 地球物理学报, 47(6): 1131-1138. |
侯重初, 蔡宗熹, 刘奎俊. 1985. 从偶层位出发建立曲面上的位场转换解释系统. 地球物理学报, 28(4): 410-418. |
刘双, 冯杰, 高文利, 等. 2013. 强剩磁强退磁条件下的二维井中磁测反演. 地球物理学报, 56(12): 4297-4309. DOI:10.6038/cjg20131232 |
刘天佑, 杨宇山, 李媛媛, 等. 2007. 大型积分方程降阶解法与重力资料曲面延拓. 地球物理学报, 50(1): 290-296. |
庞旭林. 2012. 航磁异常数据曲面延拓等效源法技术研究[硕士论文]. 北京: 中国地质大学(北京). http://cdmd.cnki.com.cn/Article/CDMD-11415-1012364846.htm
|
王万银, 潘作枢, 李家康. 1991. 三维高精度重磁位场曲面延拓方法. 物探与化探, 15(6): 415-422. |
王万银, 潘作枢. 1996. 偶层位曲面位场数据处理及转换方法. 西安地质学院学报, 18(3): 69-76. |
吴晓平. 1984. 局部重力场的点质量模型. 测绘学报, 13(4): 250-258. |
谢汝宽, 王平, 段树岭, 等. 2015. 航磁梯度数据曲化平分析. 地球物理学进展, 30(6): 2836-2840. DOI:10.6038/pg20150650 |
徐世浙, 沈晓华, 邹乐君, 等. 2004. 将航磁异常从飞行高度向下延拓至地形线. 地球物理学报, 47(6): 1127-1130. |
于德武. 2004. 用等效磁源法进行磁异常转换. 物探化探计算技术, 26(2): 133-135. |