2020, Vol. 63
2. 中国地质大学(武汉)地球物理与空间信息学院 地球内部多尺度成像湖北省重点实验室, 武汉 430074;
3. 中国地质大学(武汉)地质过程与矿产资源国家重点实验室, 武汉 430074
2. Hubei Subsurface Multi-scale Imaging Key Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
3. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan 430074, China
岩石圈磁异常是由岩石圈内目标地质体与周围环境间存在的磁学性质差异所引起, 同样利用对应的磁异常可以获取并研究该磁性差异, 进而分析目标体的组成、结构及其演化历史.当前磁法观测量主要包括总磁场强度及其三方向梯度、地磁场三分量和磁场梯度张量四种数据类型.其中, 应用最早且最广泛的观测量是总磁场强度, 其测量值与地磁正常场强度幅值的差值就是总强度磁异常.地磁场三分量观测由于受搭载平台稳定性影响较大, 目前主要见于卫星平台和海洋拖曳平台观测(Cohen and Achache, 1990; Engels et al., 2008; 吴招才等, 2011).相对于总磁场强度及其三分量观测来说, 总磁场梯度测量和磁场梯度张量测量具有受地磁场日变和搭载平台影响较小、包含信息量更丰富等方面的优势, 同时梯度数据本身对于浅部场源也具有更强的分辨能力, 因而应用前景更为广阔, 目前国内外也已有多个投入实际应用的案例(Christensen and Rajagopalan, 2000; 管志宁等, 2002; Schmidt and Clark, 2006; 张青杉等, 2010; 郭华等, 2013; Luo et al., 2015).
对于仅有总强度磁异常数据的研究对象来说, 为了从已知异常数据中获取更多关于场源体的信息, 前人提出了多种方法对总强度磁异常进行转换计算来获取其相应的化极异常、磁异常三分量、三方向梯度、磁异常梯度张量以及其他多种参量的非线性组合如磁异常模量、总梯度模、归一化磁源强度等.这些方法主要分为两类, 其一是基于傅里叶变换或希尔伯特变换来进行磁场分量转换和导数计算(Bhattacharyya, 1965; Lourenco and Morrison, 1973; Nabighian, 1984).该类方法要求待转换的数据分布在水平规则格网之上, 因而适用范围相对较窄, 但优势在于转换过程无需大量的存储空间和计算时间.其二则是等效源方法(Equivalent Source Method, ESM), 该类方法最初由Dampney(1969)提出, 现已被应用于大量实际工作中, 不仅包括上面所述的分量转换和导数计算(Emilia, 1973; Li S L and Li Y G, 2014), 还包括数据去噪(Salem et al., 2010)、数据曲化平(Hansen and Miyazaki, 1984; Nakatsuka and Okuma, 2006)等多种数据处理工作.等效源方法的原理是通过构建等效场源, 利用磁场与场源的物理关系, 以等效场源正演得到的预测数据与当前实测数据的拟合差大小作为评价标准, 如果两者拟合差小到可接受范围, 那么此时的等效场源就可以用于正演计算其他目标参量.由于等效源方法对于数据的观测高度是否统一、数据分布是否均匀没有特殊要求, 因而具有更广泛的应用范围, 但对于存储空间和计算时间的需求则相对较大.
目前上述两类方法的正常应用均有一个假设的前提, 即认为总强度磁异常近似等于磁异常矢量在地磁正常场方向上的分量.而该假设的成立及其后续应用则需依托于两个条件(Blakely, 1996)：磁异常矢量的幅值远远小于地磁正常场的幅值, 和地磁正常场的方向在研究区域内可近似认为是恒定值.将总强度磁异常近似当作磁异常矢量的分量具有许多优势.第一, 其从无物理意义的模量差值成为了具备物理意义的磁异常矢量在特定方向的分量, 为其与磁异常矢量另外三个固定方向分量之间的转换和相同场源不同磁化强度方向之间的转换提供了理论基础.第二, 由于正常场方向通常已知, 该方向上的分量则是三个固定方向上分量的线性组合, 因而构建场源体和磁异常之间的物理关系就比较简单, 尤其体现在(等效)场源物性反演过程中.第三, 结合研究区域正常场方向恒定, 可得总强度磁异常作为正常场方向的分量满足拉普拉斯方程, 因此其是调和势场, 为计算其深度方向上的梯度提供了理论基础.一般来说, 对于磁法勘探的研究尺度以上两个条件都是成立的(Blakely, 1996; 管志宁, 2005), 因而采用上述常规方法就可以满足相应数据处理工作的需求.
但是, 对于具备一定规模高磁性场源体的特殊地区, 其对应总强度磁异常幅值可能高达10000 nT.如此高幅值相比于地磁正常场全球平均值(约50000 nT)来说, 会使得上述第一个条件无法满足, 进而导致总强度磁异常近似等同于磁异常矢量在正常场方向上的投影的假设不再成立.袁晓雨等(2015)对该高幅值磁异常计算过程中存在的误差进行了较为详细的分析研究, 其研究结果明确表明在强磁性体引起的高幅值总强度磁异常的数据处理、反演和解释过程中, 现有方法会产生较大的误差.袁晓雨(2016)增加了计算误差对于磁场转换和反演结果的影响分析, 同时指出基于严格无误差的方式处理总强度磁异常是尚需解决的问题.
本文提出一种基于等效源的方法来尝试解决这一问题.本方法不通过数学上的近似来赋予总强度磁异常物理含义, 而是直接利用其数学计算过程(即地磁场总强度模量与正常场强度模量的差值), 来构建等效场源与总强度磁异常之间的数学物理关系, 进而利用最优化理论求取等效源并最终正演计算其他转换参量.本文将首先介绍该方法的基本原理和求解过程, 然后通过模型试验证实该方法的正确性和有效性, 最后将其应用于实测高值磁异常数据并给出相关讨论和最终结论.
1 方法理论独立于等效源方法的具体计算过程, 等效源所选取的等效场源类型和分布同样是该方法极为重要的一个方面.针对不同的研究范围、尺度和目标, 合适的等效源类型和分布选择也不尽相同.针对本文主要研究对象为磁法勘探, 通常研究区域范围和深度较小, 本文选择与三维物性反演相同的场源剖分方式, 将稍大于数据覆盖范围及一定深度范围内的整个空间划分为紧邻且尺寸相同的长方体.Li S L和Li Y G(2014)采用这种等效源设置方案用于起伏地形下由总强度磁异常转换计算磁异常总模量, 并给出以此方式划分场源的优势在于能够同时恢复观测数据中的低频和高频信息, 而且相对于层状等效源来说也省去了如何确定最优等效层深度的难题.此外, 本文选择三维空间等效源也是为了便于后续将该方法应用于高幅值总强度磁异常的三维磁化率成像计算过程.在确定场源构建方式之后, 下面将介绍新方法的改进之处和计算流程.
1.1 总强度磁异常及其反演雅可比矩阵为了区分两个数值相近的磁异常, 本文将总强度磁异常记作ΔT1, 磁异常矢量在正常场方向上的投影记作ΔT2.假设观测点处磁异常三分量为(Hax, Hay, Za), 研究区域正常场方向为固定值(I0, D0), 且正常场三个分量分别为(H0x, H0y, Z0), 此时有
|
(1) |
以及
|
(2) |
由于磁异常三分量与场源之间的关系是线性关系, 即第i个观测点的三分量数据, 是所有场源模型体在该观测点处正演所得对应三分量的线性相加之和, 此过程可以用式(3)来表示：
|
(3) |
其中, M指空间剖分模型单元总个数; mj指第j个模型单元的磁化率; Gx、Gy、Gz分别是指磁异常三分量的正演敏感度核矩阵.在式(3)基础上, 结合式(2)可知ΔT2与场源模型之间同样是线性正演关系, 因此常规的等效源方法利用该线性关系根据异常数据反算等效源分布.然而根据式(1)则可知ΔT1是磁异常三分量的非线性组合, 因而其与场源模型之间的关系也就是非线性的, 利用其反演计算等效源分布的过程也就比常规方法稍显复杂.
将式(3)代入式(1), 令正常场强度幅值T0i=
|
(4) |
其对第j个模型单元磁化率mj的一阶偏微分Jij即为
|
(5) |
由式(5)可见, 该偏微分结果是与所有物性参数相关的非独立量, 也就是说非线性关系中观测数据对于特定单元物性参数的敏感度还与其他单元的物性参数相关.敏感度核矩阵J(雅可比矩阵)即由此一阶偏微分Jij构成.
1.2 总目标函数构建及其求解过程基于等效源方法的总目标函数中, 用于量化及评价等效磁化率正演的预测数据dpre与观测数据dobs之间差异的数据拟合差函数Φd(m)是最重要的一项.虽然对于等效源方法来说, 位场的等效性和结果的多解性对于方法本身具有正向作用, 但是实际观测数据中的噪声仍是难以避免的.为了压制数据中噪声的影响, 根据正则化思想加入模型稳定项Φm(m)来增强反演过程的稳定性同样非常必要(Aster et al., 2005).本文参照Li和Oldenburg(1996, 1998)构建的总目标函数Φ(m)见式(6).
|
(6) |
其中, Wd是对角矩阵, 其中元素对应观测数据的误差值的倒数; F(m)表示由式(4)确定的模型空间物性参数m与预测数据dpre之间的数学关系; mref为参考模型; Wm由所选择的模型目标函数类型所确定, 详细形式及其中具体参数的含义可参考Li和Oldenburg(1996, 1998); λ是平衡两者相互权重的正则化参数.
由于式(6)中F(m)与m的关系是非线性的, 因此求取总目标函数Φ(m)最小的最优化过程也是非线性问题.针对该非线性最优化问题, 本文采用加权迭代最小二乘方法(Iteratively Reweighted Least Squares, IRLS)来求解, 其基本思路就是将非线性问题转化为多次线性问题的迭代累加.
具体求解思路为采用第(n-1)次迭代的观测数据拟合差计算第n次迭代计算的物性模型修正量.假设初始磁化率模型为m0; 第(n-1)(n≥1)次迭代的物性模型为m(n-1), 其正演预测结果为d(n-1)=F(m(n-1)), 则第n次迭代计算过程可视为用物性模型变化量Δm拟合观测值dobs与d(n-1)的差值, 因而第n次迭代反演目标函数由式(6)演化为
|
(7) |
式(7)对Δm求偏导并令等式为零, 即得到求取第n次Δm的最终表达式(8), 进而第n次最终结果可由mn=m(n－1)+Δm计算得到.如此迭代计算, 直至满足收敛条件(如Δm中绝对值最大项不超过0.001以及总目标函数变化量小于当前总目标函数的1%)则终止计算, 当前物性结果即可作为最终等效源物性结果输出, 之后便可用于正演计算其他转换异常.
|
(8) |
由于正则化参数决定了数据拟合差函数和模型目标函数两者的相对重要性, 因而该参数的确定是计算过程中非常重要的一个步骤.目前常见的确定方法有两种, L曲线法和广义交叉验证法(Li and Oldenburg, 2003; Farquharson and Oldenburg, 2004).这里需要特别说明的是, 等效源方法更加侧重于恢复数据, 因而该方法所使用的正则化参数相对于三维物性反演过程中所确定的最优正则化参数来说, 建议设置的更小一些, 以增加式(6)中数据拟合差函数的相对权重.
约束条件主要是指物性范围约束.对于等效源方法来说, 等效场源对于其磁化率的正负没有特殊要求.但若将该方法的应用延伸至磁化率反演成像, 则需要考虑加入正值约束.加入正值约束的目的不仅在于这样符合真实的物理现象, 还在于正值约束对于磁化率成像结果的聚焦程度具有正向的促进作用.因此, 如果满足加入正值约束的条件, 即使是等效源方法对等效源磁化率正负没要求, 本文也建议加入正值约束.位场反演计算所用的物性范围约束方法较多, 例如强制约束法(Boulanger and Chouteau, 2001; 姚长利等, 2007)、对数变换法(Barbosa et al., 1999; Kim et al., 1999)、对数障碍法或称内点法(Li and Oldenburg, 2003; Lelièvre and Oldenburg, 2009a)以及乘子法(Zhang et al., 2015)等.本文通过对上述方法的应用效果进行试验对比, 认为对数障碍法中额外加入的对数障碍项会使整个计算过程的收敛流程更为稳定, 所以推荐使用该方法用于物性范围约束.方法具体实现思路可参考Li和Oldenburg(2003)获取详细信息.
2 模型试验模型试验选用两个长方体作为场源, 场源形态和位置如图 1所示.对于磁异常幅值较高的情况, 以图 1中所示的场源分布, 需要两个异常体本身处于高磁状态, 即具有较高的磁化率或总磁化强度幅值, 并且通常会伴随有较强的剩余磁化强度和自退磁现象.模型试验过程中均假设两个异常体存在强剩磁或自退磁现象, 且分别考虑了两种情况：两个异常体总磁化强度方向一致并假设该方向已知, 和两个异常体总磁化强度方向不一致且数据处理过程中均假定该方向未知.由于两种情况所得数据处理和分析结果均非常相似, 所以本文仅选取相对更复杂的第二种情况进行说明和结果展示.
|
图 1 模型设置三维展示图(两个异常体等效磁化率均为3 SI) Fig. 1 A synthetic example containing two cubes both with an effective susceptibility of 3 SI |
模型试验将首先对比ΔT1和ΔT2的差异; 然后采用两种常规处理方法(频率域转换和等效源方法)对异常ΔT1进行转换和计算, 选用转换和计算得到的磁异常垂直分量(Za)、磁异常总模量(Ta)和归一化磁源强度(Normalized Source Strength, NSS)三个参量, 与理论值作差来评价常规方法计算结果; 再然后利用本文所提方法对ΔT1进行处理, 同样计算上述三个参量及其与理论值的差异; 最后通过对以上差异值进行统计和对比分析, 对新方法的应用效果进行评价.
此处选择用于对比的三个参量的转换和计算过程包括了分量转换、模量计算和梯度计算, 以尽量保证对比的全面性.另外, 需要明确的是虽然化极异常是一种非常重要的转换异常, 但是对于模型试验第二种情况, 即异常体存在剩磁但磁化方向未知且与正常场方向不一致的条件下, 常规处理方法如基于频率域转换的方法是难以计算得到正确化极异常的, 因而采用磁异常垂直分量Za来代替化极异常进行对比分析工作.
本文常规处理方法中频率域转换算法采用Gerovska和Araúzo-Bravo(2006)提供的Matlab程序MMTrans, 并在程序原有基础上对其进行拓展加入了计算NSS的相关内容, 有关程序已开源可供下载(https://github.com/KICIOLLO/MMTrans-update); 常规等效源方法由于依托于总强度磁异常与场源磁化率分布的正演关系是线性的近似假设, 因此常规等效源方法是线性反演过程, 相当于是本文新方法中多次迭代线性过程中的一次, 实现过程相对简单.
首先, 令两个异常体的等效磁化率均为3 SI, 西侧异常体的总磁化强度方向为(Im, Dm)=(45°, 10°), 东侧异常体的总磁化强度方向为(Im, Dm)=(-45°, 80°); 假设当地的地磁正常场强度幅值为57000 nT, 倾角和偏角分别为45°和10°, 相当于地磁正常场三分量(H0x, H0y, Z0)=(39693 nT, 6999 nT, 40305 nT).假设观测高度是地面0 m以上25 m, 正演计算得到的相关磁异常及参量如图 2所示, 包括ΔT1(图 2a)、ΔT2(图 2b)及其两者的差异值(图 2c)、磁异常垂直分量Za(图 2d)、磁异常总模量Ta(图 2e)和归一化磁源强度NSS(图 2f).从图 2中可以看出, 由于强剩磁的存在, 东侧异常体对应的ΔT1、ΔT2和Za的形态以负值为主, 且其形态分布特征明显与仅存在感应磁化强度的情况有较大差异; 相对而言, Ta和NSS则与场源体水平位置对应良好, 且该两个参量所受场源体磁化强度方向影响较小的特征, 也是其被应用于强剩磁和强自退磁效应下定性和定量解释的重要原因(Li et al., 2010; Pilkington and Beiki, 2013; Liu et al., 2013; Krahenbuhl and Li, 2017).还需说明的一点是, 图 2c中ΔT1与ΔT2的差值是非负的, 袁晓雨等(2015)从数学上证实了该差值的非负性, 并结合模型试验认为如将ΔT1当作ΔT2进行后续处理如磁化率三维成像会得到偏大于真实值的结果.
|
图 2 图 1中两个异常体总磁化强度方向不同时模型正演磁异常 (a)总强度磁异常(ΔT1); (b)磁异常矢量在地磁正常场方向的投影(ΔT2); (c) ΔT1与ΔT2的差异; (d)磁异常矢量垂向分量(Za); (e)磁异常总模量(Ta); (f)归一化磁源强度(NSS).各图中黑色实线框标示了异常体水平方向位置和范围. Fig. 2 The forward modeling of the synthetic model shown in Fig. 1 when the magnetization direction of two cubes are different (a) The total-field magnetic anomaly (ΔT1); (b) The projection of the anomalous field vector onto the direction of the geomagnetic reference field (ΔT2); (c) The difference between ΔT1 and ΔT2; (d) The vertical component of the anomalous field vector (Za); (e) The total amplitude of the anomalous field vector (Ta); (f) The normalize source strength (NSS).The black solid lines in each figure delight the horizontal locations of the two cubes. |
由于本试验中两个异常体具有不同的磁化方向, 而多个异常体的不同磁化强度方向通常难以利用已有特定方法(Dannemiller and Li, 2006; Gerovska et al., 2009; Rao et al., 2016; Li et al., 2017)估算得到, 因而此处假设等效源磁化强度方向是当地地磁正常场的方向.而在此假设条件下, 等效源计算过程正值约束条件是否加入则需要慎重考虑.比如该试验中的总强度磁异常ΔT1(图 2a)中东侧的高值负异常, 由于该高值负异常南方向没有伴生相应的高值正异常, 如采用正常场方向的正值场源正演结果则难以有效地拟合该数据.因此, 针对该模型试验, 在等效源磁化方向未知的情况下, 无论常规等效源方法还是本文新方法, 均未对等效场源添加正值约束条件.
然后, 根据实际资料处理流程得到的异常ΔT1按照常规处理思路, 将其当作是ΔT2异常的近似值, 分别采用两种常规方法, 即基于频率域的转换计算方法和常规的等效源方法对ΔT1异常进行处理.此处需要说明的是, 本文不探讨不同方法应用过程中观测噪声对于方法处理结果的影响, 在后续处理过程中所用ΔT1异常均未包含噪声, 在常规等效源方法和新方法反演计算等效源时假设所有观测点数据观测误差为1 nT; 并且后续的转换结果不再展示异常图, 仅展示转换和计算结果与图 2中对应理论异常值的差异.
利用频率域转换方法得到的Za、Ta和NSS异常与理论值的差异分别如图 3a、3b和3c所示; 常规等效源方法所采用的正则化参数为102.0, 得最终数据拟合差函数值为4.6×10-6, 之后利用所得等效场源正演计算得到的Za、Ta和NSS异常与理论值的差异分别如图 3d、3e和3f所示.各个差异结果对应的异常范围、均值和标准差分列于表 1.两种常规方法对于磁异常分量Za和总模量Ta的转换计算结果, 与真实值的差异在形态上较为相似, 两者幅值范围均约为600 nT, 且其标准差均在70 nT左右; NSS由于在计算磁异常分量的基础上增加了梯度运算, 而梯度运算过程对于转换计算过程中的差异具有增强效果, 导致两种常规方法得到的NSS结果与理论值的差异在形态特征上具有较大的差异.
|
图 3 ΔT1常规方法处理结果与图 2中对应真实值的差异 (a), (b)和(c)分别是采用频率域方法转换计算得到的Za, Ta和NSS结果与图 2中对应真实异常的差值; (d), (e)和(f)则分别是采用常规等效源方法获取等效源后再次正演计算得到的Za, Ta和NSS结果与图 2中对应真实异常的差值. Fig. 3 The differences between the true anomalies shown in Fig. 2 and their corresponding transformed or calculated anomalies of ΔT1 using general methods (a), (b), and (c) are the differences between Za, Ta, and NSS, which are transformed using method based on FFT, and their corresponding true values shown in Fig. 2, respectively; (d), (e), and (f) are the differences between Za, Ta, and NSS, which are calculated using general equivalent source method, and their corresponding true values shown in Fig. 2, respectively. |
|
表 1 模型试验不同方法处理结果与对应真实值差异统计表 Table 1 Statistics of the differences between the true anomalies and their corresponding calculated anomalies from ΔT1 using different methods for synthetic model |
最后, 对ΔT1异常采用新方法计算等效源磁化率分布, 同样采用正则化参数为102.0, 得最终数据拟合差函数值为4.1×10-3; 之后利用等效源正演计算三个参量值, 并与真实值作差, 如图 4所示.为了与图 3中结果更清晰地对照, 成图时选用了与图 3中相应参量相同的色标.图 4非常明显地展示了Za(图 4a)和Ta(图 4b)与对应真实值差异的幅值更集中于零值附近; 表 1中数据统计结果表明, 两个参量结果的幅值范围和标准差均仅有常规方法处理结果的1/5~1/10.以上试验结果充分说明了相对于常规转换和计算方法, 新方法在计算磁异常分量和总模量方面能够获取与真实值更为接近的结果.但是, 对于参量NSS来说, 新方法获得的结果在异常形态和幅值上与常规方法处理的结果相比均没有明显的提升.这也是需要注意的一点, 本文新方法用于转换计算与原始观测异常同等阶次的异常时, 能够提高最终结果的质量; 但是在计算原始观测数据的高阶梯度异常时, 可能由于离散点观测数据信息量不充足的原因, 相较于常规方法则并未明显提高计算结果的质量.
|
图 4 ΔT1新方法处理结果与图 2中对应真实值的差异 (a), (b)和(c)分别是采用新方法计算得到的Za, Ta和NSS结果与图 2中对应真实异常的差值. Fig. 4 The differences between the true anomalies shown in Fig. 2 and their corresponding calculated anomalies of ΔT1 using the new method (a), (b), and (c) are the differences between Za, Ta, and NSS, which are calculated using the new equivalent source method, and their corresponding true values shown in Fig. 2, respectively. |
本文研究区对应的磁异常最大幅值超过10000 nT, 是非常良好的新方法试验数据.该部分将选择磁异常总模量Ta用于数据处理结果之间的对比工作; 处理方法和过程分别选择频率域方法转换计算和新方法求取等效源并正演计算.由于实际研究区并无Ta真值可供参考, 因此仅根据上述两种方法处理结果的差值来体现和反映两种方法处理结果的差异程度和新方法应用的必要性.
该铁矿为大型隐伏条带状铁建造(Banded Iron Formation, BIF)矿体, 铁矿储量估计达到10亿吨(Wu et al., 2015).图 5所示为该铁矿部分地面实测总强度磁异常, 可见其磁异常幅值范围为-3000 nT至12000 nT; 当地地磁正常场强度幅值为54000 nT, 倾角和偏角分别为57.9°和-7.1°.
|
图 5 研究区实测总强度磁异常 Fig. 5 The total-field magnetic anomaly over the study region |
根据图 5所示磁异常和正常场方向, 采用频率域方法计算得到的磁异常总模量Ta如图 6a所示.继而利用新等效源方法, 采用正值约束, 三维等效源模型空间剖分为65×91×34=201110个大小均一、边长为50 m的正方体单元.令数据误差均为1 nT, 采用正则化参数为106.0; 经新方法等效源反演计算, 得最终数据拟合差函数值为723, 将反演得到等效源进行正演计算得到的Ta异常如图 6b所示.两种方法所得Ta结果的差异见图 6c, 图中显示两者差异值最大幅值可高达800 nT, 较好地体现了两种方法处理结果的区别.
|
图 6 频率域方法转换Ta结果(a)与新方法计算Ta结果(b)及两者差异(c) Fig. 6 The transformed or calculated Ta anomalies using method based on FFT (a) and new equivalent source method (b) and the difference (c) between the above two Ta anomalies |
在方法理论部分讨论等效源设置时, 已说明本文中选用的设置方案是三维物性反演所用的三维等效源设置方案; 在确定新方法总目标函数的过程中, 同样是按照物性反演的思路来构建的模型稳定项(或称模型目标函数).因此, 本文中提出的新方法, 理论上可直接用于三维磁化率反演成像.为了保证计算得到的三维磁化率分布结果符合真实物理规律, 能够提供关于场源体在形态和磁化率相对强弱等方面的信息, 新方法应用于三维磁化率反演还需要满足如下特定条件：整个反演空间场源体的总磁化强度方向已知, 反演空间深度范围能够包含所有目标体, 所有观测点位置对应的正常场强度和方向已知, 以及三维反演计算过程中对磁化率施加正值约束.
针对场源体总磁化强度方向较为复杂的情况, 当前已有学者提出了相应的磁化强度矢量反演方法(如Lelièvre and Oldenburg, 2009b; Li and Sun, 2016; Liu et al., 2017).由于总强度磁异常受该磁化强度矢量方向的影响比较大, 因而是用于反演计算磁化强度矢量较好的源数据之一.并且前述所列方法同样均未考虑高幅值异常的情况.因此与本文所提新方法采用相同的处理方式, 将数学上直接获取总强度磁异常的过程应用于已有磁化强度矢量反演方法中也是后续将开展的内容.
由于本文主要研究对象是磁法勘探, 研究区域通常较小, 整个区域内地磁正常场方向可认为是恒定值, 因此所介绍的主要研究内容是磁异常总模量幅值较高时对于数据处理的影响及其解决方案.然而, 当研究对象尺度较大时, 整个研究区域内地磁正常场方向可能存在较大的变化, 其直接影响就是导致总强度磁异常不再满足拉普拉斯方程(Blakely, 1996), 而此时对于磁异常分量转换、化极、梯度和模量计算等处理工作均应采用等效源方法进行计算.本文所提新方法同样适用于该种情况.但是通常尺度较大的研究对象也对应了较大的研究区域范围和深度, 如需获取空间分辨率较高的物性分布结果, 对于计算硬件需求也就更高; 同时由于新方法相对于常规等效源方法来说需要的内存量更大(前者所需内存是后者的4倍), 并且前者的非线性反演过程所花费时间也更多.再者, 尺度较大的研究对象对应的磁异常由于观测高度较高因而幅值相对较低, 如大地水准面以上4 km高度处EMAG2(Maus et al., 2009)磁异常幅值通常介于±400 nT之间.这种条件下常规等效源方法和新方法处理结果的差异尚未有学者分析.因此, 在此条件下是否一定要采用新方法也将是后续的研究内容之一.
5 结论本文提出了一种基于等效源的总强度磁异常非线性处理方法.该方法不依赖于常规条件和假设, 而是根据总强度磁异常的获取过程, 即实测总磁场强度幅值和对应正常场强度幅值之差, 来构建目标函数用于求取等效源, 尤其适用于具有高幅值总强度磁异常的情况.在对新方法的总目标函数构建思路、求解过程和关键技术细节进行介绍的基础上, 开展模型试验对新方法的计算结果进行验证.模型试验结果表明, 对于Za、Ta转换异常来说, 采用新方法处理结果与真实值差异, 为常规方法所得结果与真实值差异的1/5甚至更小; 结合实际观测资料的处理结果共同体现了采用新方法用于高幅值磁异常数据处理的有效性和必要性.同时, 本文也对该方法相应的总强度磁异常处理思路应用于三维磁化率反演成像、磁化强度矢量反演、大区域磁化率成像等工作的可行性进行了讨论, 与之相关问题将作为后续内容开展进一步研究.
致谢 两位匿名审稿人对本文提供了非常详细且有效的建议; 本文所用实际磁测数据来源于中国冶金地质总局矿产资源研究院; 制图所用工具包括Gmsh(Geuzaine and Remacle, 2009)和GMT(Wessel et al., 2013)两款开源软件; 张壹、张双喜对于本文开源程序MMTrans-update贡献了部分代码; 牛向龙博士、梁青博士和胥鸿睿博士对于本文内容和结构提出了宝贵意见和建议.特此致谢.
Aster R C, Borchers B, Thurber C H. 2005. Parameter Estimation and Inverse Problems. London, UK: Academic Press.
|
Barbosa V C F, Silva J B C, Medeiros W E. 1999. Gravity inversion of a discontinuous relief stabilized by weighted smoothness constraints on depth. Geophysics, 64(5): 1429-1437. DOI:10.1190/1.1444647 |
Bhattacharyya B K. 1965. Two-dimensional harmonic analysis as a tool for magnetic interpretation. Geophysics, 30(5): 829-857. DOI:10.1190/1.1439658 |
Blakely R J. 1996. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press.
|
Boulanger O, Chouteau M. 2001. Constraints in 3D gravity inversion. Geophysical Prospecting, 49(2): 265-280. DOI:10.1046/j.1365-2478.2001.00254.x |
Christensen A, Rajagopalan S. 2000. The magnetic vector and gradient tensor in mineral and oil exploration. Preview, 84: 77. |
Cohen Y, Achache J. 1990. New global vector magnetic anomaly maps derived from Magsat data. Journal of Geophysical Research:Solid Earth, 95(B7): 10783-10800. DOI:10.1029/JB095iB07p10783 |
Dampney C N G. 1969. The equivalent source technique. Geophysics, 34(1): 39-53. DOI:10.1190/1.1439996 |
Dannemiller N, Li Y G. 2006. A new method for determination of magnetization direction. Geophysics, 71(6): L69-L73. DOI:10.1190/1.2356116 |
Emilia D A. 1973. Equivalent sources used as an analytic base for processing total magnetic field profiles. Geophysics, 38(2): 339-348. DOI:10.1190/1.1440344 |
Engels M, Barckhausen U, Gee J S. 2008. A new towed marine vector magnetometer:methods and results from a Central Pacific cruise. Geophysical Journal International, 172(1): 115-129. DOI:10.1111/j.1365-246X.2007.03601.x |
Farquharson C G, Oldenburg D W. 2004. A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems. Geophysical Journal International, 156(3): 411-425. DOI:10.1111/j.1365-246X.2004.02190.x |
Gerovska D, Araúzo-Bravo M J. 2006. Calculation of magnitude magnetic transforms with high centricity and low dependence on the magnetization vector direction. Geophysics, 71(5): I21-I30. DOI:10.1190/1.2335516 |
Gerovska D, Araúzo-Bravo M J, Stavrev P. 2009. Estimating the magnetization direction of sources from southeast Bulgaria through correlation between reduced-to-the-pole and total magnitude anomalies. Geophysical Prospecting, 57(4): 491-505. DOI:10.1111/j.1365-2478.2008.00761.x |
Geuzaine C, Remacle J F. 2009. Gmsh:A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11): 1309-1331. DOI:10.1002/nme.2579 |
Guan Z N, Hao T Y, Yao C L. 2002. Prospect of gravity and magnetic exploration in the 21st century. Progress in Geophysics (in Chinese), 17(2): 237-244. DOI:10.3969/j.issn.1004-2903.2002.02.008 |
Guan Z N. 2005. Geomagnetic Field and Magnetic Exploration (in Chinese). Beijing: Geological Publishing House.
|
Guo H, Wang P, Xie R K. 2013. A study of geological interpretation with the tri-axial aeromagnetic gradients. Progress in Geophysics (in Chinese), 28(5): 2688-2692. DOI:10.6038/pg20130551 |
Hansen R O, Miyazaki Y. 1984. Continuation of potential fields between arbitrary surfaces. Geophysics, 49(6): 787-795. DOI:10.1190/1.1441707 |
Kim H J, Song Y, Lee K H. 1999. Inequality constraint in least-squares inversion of geophysical data. Earth, Planets and Space, 51(4): 255-259. DOI:10.1186/BF03352229 |
Krahenbuhl R A, Li Y G. 2017. Investigation of magnetic inversion methods in highly magnetic environments under strong self-demagnetization effect. Geophysics, 82(6): J83-J97. DOI:10.1190/geo2016-0676.1 |
Lelièvre P G, Oldenburg D W. 2009a. A comprehensive study of including structural orientation information in geophysical inversions. Geophysical Journal International, 178(2): 623-637. DOI:10.1111/j.1365-246X.2009.04188.x |
Lelièvre P G, Oldenburg D W. 2009b. A 3D total magnetization inversion applicable when significant, complicated remanence is present. Geophysics, 74(3): L21-L30. DOI:10.1190/1.3103249 |
Li J P, Zhang Y T, Yin G, et al. 2017. An approach for estimating the magnetization direction of magnetic anomalies. Journal of Applied Geophysics, 137: 1-7. DOI:10.1016/j.jappgeo.2016.12.009 |
Li S L, Li Y G. 2014. Inversion of magnetic anomaly on rugged observation surface in the presence of strong remanent magnetization. Geophysics, 79(2): J11-J19. DOI:10.1190/geo2013-0126.1 |
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. 1998. 3-D inversion of gravity data. Geophysics, 63(1): 109-119. DOI:10.1190/1.1444302 |
Li Y G, Oldenburg D W. 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method. Geophysical Journal International, 152(2): 251-265. DOI:10.1046/j.1365-246X.2003.01766.x |
Li Y G, Shearer S E, Haney M M, et al. 2010. Comprehensive approaches to 3D inversion of magnetic data affected by remanent magnetization. Geophysics, 75(1): L1-L11. |
Li Y G, Sun J J. 2016. 3D magnetization inversion using fuzzy c-means clustering with application to geology differentiation. Geophysics, 81(5): J61-J78. DOI:10.1190/geo2015-0636.1 |
Liu S, Hu X Y, Liu T Y, et al. 2013. Magnetization vector imaging for borehole magnetic data based on magnitude magnetic anomaly. Geophysics, 78(6): D429-D444. DOI:10.1190/geo2012-0454.1 |
Liu S, Hu X Y, Zhang H L, et al. 2017. 3D magnetization vector inversion of magnetic data:Improving and comparing methods. Pure and Applied Geophysics, 174(12): 4421-4444. DOI:10.1007/s00024-017-1654-3 |
Lourenco J S, Morrison H F. 1973. Vector magnetic anomalies derived from measurements of a single component of the field. Geophysics, 38(2): 359-368. DOI:10.1190/1.1440346 |
Luo Y, Wu M P, Wang P, et al. 2015. Full magnetic gradient tensor from triaxial aeromagnetic gradient measurements:Calculation and application. Applied Geophysics, 12(3): 283-291. DOI:10.1007/s11770-015-0508-y |
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. DOI:10.1029/2009GC002471 |
Nabighian M N. 1984. Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms:fundamental relations. Geophysics, 49(6): 780-786. DOI:10.1190/1.1441706 |
Nakatsuka T, Okuma S. 2006. Reduction of magnetic anomaly observations from helicopter surveys at varying elevations. Exploration Geophysics, 37(1): 121-128. |
Pilkington M, Beiki M. 2013. Mitigating remanent magnetization effects in magnetic data using the normalized source strength. Geophysics, 78(3): J25-J32. DOI:10.1190/geo2012-0225.1 |
Rao C F, Yu P, Zhang L L, et al. 2016. Estimating the magnetization direction of magnetic bodies through correlation between reduced-to-the-pole anomaly and the normalized source strength.//2016 SEG International Exposition and Annual Meeting. Dallas, Texas: Society of Exploration Geophysicists, 1627-1630.
|
Salem A, Lei K X, Green C, et al. 2010. Removal of cultural noise from high-resolution aeromagnetic data using a two stage equivalent source approach. Exploration Geophysics, 41(2): 163-169. |
Schmidt P W, Clark D A. 2006. The magnetic gradient tensor:Its properties and uses in source characterization. The Leading Edge, 25(1): 75-78. DOI:10.1190/1.2164759 |
Wessel P, Smith W H F, Scharroo R, et al. 2013. Generic mapping tools:Improved version released. Eos, Transactions American Geophysical Union, 94(45): 409-410. |
Wu H Y, Niu X L, Zhang L C, et al. 2015. Geology and geochemistry of the Macheng Algoma-type banded iron-formation, North China Craton:Constraints on mineralization events and genesis of high-grade iron ores. Journal of Asian Earth Sciences, 113: 1179-1196. DOI:10.1016/j.jseaes.2015.05.024 |
Wu Z C, Gao J Y, Luo X W, et al. 2011. Marine measurement of the three-component geomagnetic field. Progress in Geophysics (in Chinese), 26(3): 902-907. DOI:10.3969/j.issn.1004-2903.2011.03.015 |
Yao C L, Zheng Y M, Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces. Chinese J. Geophys. (in Chinese), 50(5): 1576-1583. |
Yuan X Y, Yao C L, Zheng Y M, et al. 2015. Error analysis of calculation of total field anomaly due to highly magnetic bodies. Chinese J. Geophys. (in Chinese), 58(12): 4756-4765. DOI:10.6038/cjg20151235 |
Yuan X Y. 2016. Analysis of calculation error of total field anomaly and high precision processing in strongly magnetic environments[Master's thesis] (in Chinese). Beijing: China University of Geosciences (Beijing).
|
Zhang Q S, Ma F L, Xu L Y. 2010. Study on the 3D gradient aeromagnetic survey. Geology and Exploration (in Chinese), 46(6): 1087-1091. |
Zhang Y, Yan J G, Li F, et al. 2015. A new bound constraints method for 3-D potential field data inversion using Lagrangian multipliers. Geophysical Journal International, 201(1): 267-275. DOI:10.1093/gji/ggv016 |
管志宁, 郝天珧, 姚长利. 2002. 21世纪重力与磁法勘探的展望. 地球物理学进展, 17(2): 237-244. DOI:10.3969/j.issn.1004-2903.2002.02.008 |
管志宁. 2005. 地磁场与磁力勘探. 北京: 地质出版社.
|
郭华, 王平, 谢汝宽. 2013. 航磁全轴梯度数据地质解释优势研究. 地球物理学进展, 28(5): 2688-2692. DOI:10.6038/pg20130551 |
吴招才, 高金耀, 罗孝文, 等. 2011. 海洋地磁三分量测量技术. 地球物理学进展, 26(3): 902-907. DOI:10.3969/j.issn.1004-2903.2011.03.015 |
姚长利, 郑元满, 张聿文. 2007. 重磁异常三维物性反演随机子域法方法技术. 地球物理学报, 50(5): 1576-1583. DOI:10.3321/j.issn:0001-5733.2007.05.035 |
袁晓雨, 姚长利, 郑元满, 等. 2015. 强磁性体ΔT异常计算的误差分析研究. 地球物理学报, 58(12): 4756-4765. DOI:10.6038/cjg20151235 |
袁晓雨. 2016.强磁异常ΔT的计算误差及高精度处理转换分析研究[硕士论文].北京: 中国地质大学(北京). http://cdmd.cnki.com.cn/Article/CDMD-11415-1016183927.htm
|
张青杉, 麻丰林, 许丽云. 2010. 航空三维磁梯度测量方案研究. 地质与勘探, 46(6): 1087-1091. |