基于L-BFGS反演算法的ΔT精确计算磁异常分量Tap方法

doi:10.11720/wtyht.2019.1383

基于L-BFGS反演算法的ΔT精确计算磁异常分量T_ap方法

甄慧翔, 杨宇山^,, 李媛媛, 刘天佑

中国地质大学(武汉) 地球物理与空间信息学院,湖北武汉 430074

Method for accurately calculating magnetic anomaly component using ΔT based on L-BFGS inversion algorithm

ZHEN Hui-Xiang, YANG Yu-Shan^,, LI Yuan-Yuan, LIU Tian-You

Institute of Geophysics and Geomatics, China University of Geosciences(Wuhan), Wuhan 430074, China

通讯作者: 杨宇山(1977-),男,副教授,主要从事重磁资料处理与解释方面的教学和科研工作。Email:samyys@126.com

责任编辑: 王萌

收稿日期: 2018-10-24 修回日期: 2019-01-16 网络出版日期: 2019-06-20

基金资助:

国家重点研发计划项目“华北克拉通辽东/胶东重要成矿区带金多金属矿深部预测及勘查示范”下设课题“空—地—井立体探测技术方法构建与示范”. 2018YFC0603803

Received: 2018-10-24 Revised: 2019-01-16 Online: 2019-06-20

作者简介 About authors

甄慧翔(1995-),男,汉族,硕士研究生,主要从事重磁资料处理及反演研究工作。Email:zhenhx1995@163.com. 。

摘要

磁法勘探理论中,将ΔT磁异常看作磁异常矢量T_a在地磁场方向的分量T_ap,是ΔT异常处理与解释的物理基础,然而这种近似存在误差,理论计算及实验已经证明这种近似所产生的误差将随着T_a异常强度的增大而迅速增加。当磁异常T_a远小于地磁场T₀时,误差影响小,可忽略,在强磁异常情况下,误差大,ΔT异常的处理解释精度会受到很大的影响。对于高精度磁法勘探而言,必须将ΔT转换成磁异常分量T_ap进行处理解释。笔者提出了基于有限储存BFGS(L-BFGS)反演算法的ΔT精确计算磁异常分量方法,首先推导了T_ap计算ΔT的正演公式,利用ΔT与T_ap的差值构建反演T_ap的目标函数,采用L-BFGS算法由ΔT解算T_ap。模型实验表明该方法计算得到的T_ap十分接近理论值,即可将误差降低两个数量级。在存在噪声与背景场情况下该方法也都能得到很好的结果。将本方法应用于福建阳山铁矿ΔT磁测资料的处理,得到了与实际更加符合的处理解释结果。

关键词： 强磁异常 ; 总场异常ΔT ; 磁异常分量T_ap ; 高精度处理解释 ; L-BFGS

Abstract

In the magnetic exploration theory, total-field anomaly ΔT is regarded as the component T_ap of the magnetic anomaly vector T_a on the main field (T₀) direction and thus constitutes the theoretical basis. However, there is an error in this approximation. Theoretical calculations and experiments have proved that this approximation error will increase rapidly as the T_a increases. When the magnetic anomaly T_a is much smaller than T₀, the influence of the error is small and negligible. In the case of a strong magnetic anomaly, the error is large, and the processing interpretation accuracy of the ΔT anomaly is greatly affected. For high-precision magnetic exploration, ΔT must be converted to a magnetic anomaly component Tap for processing and interpretation. In this paper, the method of accurately calculating the magnetic anomaly component using T_ap based on the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm is proposed. Firstly, the authors derived the forward formula for ΔT from Tap, and then constructed the objective function of Tap inversion by the difference function between ΔT and T_ap. L-BFGS algorithm was used to solve the Tap from ΔT. Model experiments show that the Tap calculated by this method is very close to the real value, which can reduce the error by two orders of magnitude. This method also yields good results in the presence of noise and background fields. The method was applied to the processing of ΔT magnetic survey data of the Yangshan iron mine in Fujian Province, and the results of processing and interpretation which are more consistent with the actual results were obtained.

Keywords： magnetic anomaly ; total-field anomaly ΔT ; magnetic anomaly component T_ap ; high-precision processing and interpretation ; L-BFGS

PDF (6013KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

甄慧翔, 杨宇山, 李媛媛, 刘天佑. 基于L-BFGS反演算法的ΔT精确计算磁异常分量T_ap方法. 物探与化探[J], 2019, 43(3): 598-607 doi:10.11720/wtyht.2019.1383

ZHEN Hui-Xiang, YANG Yu-Shan, LI Yuan-Yuan, LIU Tian-You. Method for accurately calculating magnetic anomaly component using ΔT based on L-BFGS inversion algorithm. Geophysical and Geochemical Exploration[J], 2019, 43(3): 598-607 doi:10.11720/wtyht.2019.1383

0 引言

20世纪70年代,ΔT测量因其具备无需调平和定向等便捷性,逐渐取代了磁异常分量Z_a测量,成为磁法勘探中的主要测量方式^[1,2]。质子磁力仪和光泵磁力仪都是测量地磁总场,然后减去地磁正常场得到ΔT。在磁法勘探理论中,总场异常ΔT,即地磁总场与地磁正常场的模量差,常被近似看作磁异常矢量T_a在地磁正常场方向上的投影分量T_ap,因此磁法处理解释中存在ΔT近似误差。磁法勘探理论^[3,4,5,6,7]指出,当磁异常远小于地磁场时,近似误差对后续处理的影响可以忽略。这种近似的优势在于,ΔT相当于是磁异常矢量T_a在地磁场方向上的投影分量,ΔT与磁性体磁化强度M之间成了简单的线性关系,而与地磁正常场T₀无直接关系。后续分量转换等计算及反演解释都将得到简化,所以在弱磁异常条件下,将ΔT近似看作T_ap便捷有效。

Blakely^[3]指出可以对ΔT近似处理的缘由是由于航空磁测测得的磁异常较小,近似误差E也较小(图1),所以在处理中忽略ΔT近似误差。磁法勘探理论教材中普遍提到,当T₀≫T_a时,误差才可以被忽略。Hinze^[6]通过实验说明当T₀为30 000 nT时,如果T_a为500 nT,误差只有4 nT,而当T_a为 10 000 nT时误差可以达到1 667 nT。实际工作条件并非都是航空磁测和卫星磁测。在具有强磁异常的条件,如埋深较浅、规模较大的矿体或井中ΔT磁测探头靠近矿体,磁异常非常大甚至超过地磁场,ΔT近似误差E不容忽略。如今磁法勘探中,工作方法往往不管磁异常T_a强度大小,普遍将总场异常ΔT看作磁异常T_a在地磁场方向上的投影分量T_ap。后续的化极、延拓、分量转换、磁异常模量计算等处理解释方法都是在这种近似下进行的,究其缘由,一方面是仍未意识到ΔT近似误差的影响,另一方面是因为缺乏代替近似处理的解决方案。袁晓雨等^[8,9]通过实验发现,随着磁异常的增大,E按近似指数关系迅速增大,在强磁异常条件下不容忽视,否则对后续的定量计算和分析解释会造成较大的影响。误差最为直观的影响就是磁异常分量转换过程,因为ΔT并不是磁异常T_a的分量,而处理解释中,将它直接看作磁异常T_a分量T_ap,误差是必然存在的。其他处理解释方法如化极、磁异常模量求取等都与分量转换密切相关。垂直磁化条件下,磁异常的形态以及磁异常与磁性体的关系都比较简单,便于进行地质解释,因此处理中常常对ΔT磁异常作化到地磁极处理。化极理论上也应是对地磁场方向上的投影分量T_ap进行处理^[10,11],所以ΔT化极存在近似误差影响。图2所示,对一个单球体模型强磁异常ΔT和T_ap化极结果对比。模型球体中心坐标(0,0),埋深150 m,半径为130 m,磁化强度50 A/m,地磁场强度T₀为50 000 nT,地磁倾角45°,偏角0°。图2可以看出,强磁异常情况下,T_ap化极结果良好。ΔT化极后,不仅在磁性体处有较大误差,在磁性体周围也存在虚假异常,异常体的位置形状都受到影响,异常整体向磁偏角方向偏移,对后续的精细化处理解释有很大影响。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 总场异常ΔT及其他相关物理量关系示意

Fig.1 Schematic diagram of total field anomaly ΔT and other related physical quantities

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 强磁异常情况下ΔT与T_ap化极结果对比

a—球体模型ΔT;b—球体模型T_ap;c—球体模型ΔT化极结果;d—球体模型T_ap化极结果;红圈表示球体模型投影,黑点表示球体模型的中点

Fig.2 Reduction to the pole (RTP) results comparison of ΔT and T_ap under strong magnetic anomalies

a—ΔT of the sphere model; b—T_ap of the sphere model; c—RTP of the sphere model ΔT; d—RTP of the sphere model T_ap;the red circle in the figure represents the projection of the sphere model and the black dot represents the midpoint of the sphere model

为此笔者提出一种方法,不对ΔT进行近似处理,直接计算得到T_ap,提高处理解释精度。笔者首先推导了ΔT与T_ap的差值公式,由此得到根据T_ap计算ΔT的正演公式。然后利用正演公式构建了反演磁异常分量T_ap的目标函数,最后利用L-BFGS优化算法解该非线性优化问题得到磁异常分量T_ap。

1 方法原理

1.1 误差原理

下面根据磁法勘探理论,推导误差公式。地磁场总场T是地磁正常场(后面简称正常场)T₀与磁异常T_a的矢量叠加

(1)

T = T_{0} + T_{a},

总场异常ΔT与磁异常T_a不同,ΔT是地磁场总场T与正常场T₀的模量差

(2)

Δ T = | T | - | T_{0} | = T - T_{0},

总场异常ΔT与T_ap的差值E,也就是近似处理导致的误差,表示为

(3)

E = Δ T - T_{ap} 。

图1中,由三角余弦定理可知

(4)

T = \sqrt[]{T_{0}^{2} + T_{a}^{2} + 2 T_{0} T_{a} \cos θ},

由式(2)、式(4)可以得到

(5)

Δ T^{2} + T_{0}^{2} + 2 Δ T T_{0} = T_{0}^{2} + T_{a}^{2} + 2 T_{0} T_{ap},

(6)

2 T_{0} (Δ T - T_{ap}) = T_{a}^{2} - Δ T^{2},

然后通过式(3)和式(6)得到误差公式

(7)

E = \frac{T_{a}^{2} - Δ T^{2}}{2 T_{0}} 。

根据误差公式,结合三角关系可知

(8)

T_{a} \geq Δ T,

由此可以得到两个性质

(9)

E \geq 0,

(10)

Δ T \geq T_{ap} .

1.2 构建目标函数

目标函数的构建是结合已知的误差公式(7),得到一个可以由磁异常分量T_ap计算ΔT的正演公式。根据正演公式来进行反演计算。

由误差公式可知,对于实测数据,任意一个测点都满足下式

(11)

Δ T = T_{ap} + \frac{T_{a}^{2} - Δ T^{2}}{2 T_{0}} 。

式(11)中,在不大的区域内正常场T₀可看作固定方向的常数,除了ΔT、T_ap,还存在变量T_a。根据位场理论磁异常T_a与T_ap可以相互转换,假设有m×n二维平面网格上总测点数为N,通过磁异常分量换算,利用磁异常分量T_ap平面数据在频率域乘上转换因子,得到磁异常T_a三个方向分量H_ax、H_ay、Z_a。然后通过矢量合成可求得T_a平面数据。过程如下:

(12)

F [T_{ap}] = \int_{-}^{\infty} \int_{-}^{\infty} T_{ap (x, y, 0)} e^{- uix} e^{- viy} d x d y,

(13)

\begin{array}{l} H_{ay} = f (\frac{2 π iv}{q_{t_{0}}} F [T_{ap}]), \\ Z_{a} = f (\frac{2 π i (u^{2} + v^{2})^{1 / 2}}{q_{t_{0}}} F [T_{ap}]), \\ T_{a}^{2} = H_{ax}^{2} + H_{ay}^{2} + Z_{a}^{2}, \end{array}

(14)

T_{a}^{2} = f {(\frac{2 π iu}{q_{t_{0}}} F [T_{ap}])}^{2} + f {(\frac{2 π iv}{q_{t_{0}}} F [T_{ap}])}^{2} + f {(\frac{2 π i (u^{2} + v^{2})^{\frac{1}{2}}}{q_{t_{0}}} F [T_{ap}])}^{2}

, q_{t_{0}} = 2 π [i (L_{0} u + M_{0} u) + N_{0} (u^{2} + v^{2})^{\frac{1}{2}}] 。

其中:F[]表示二维傅里叶正变换,f()表示二维傅里叶反变换,L₀、M₀、N₀为地磁场单位矢量t₀的方向余弦。

上面介绍了T_ap计算T_a的过程,式(11)中 $T_{a}^{2}$ 可由式(14)计算得到,由此得到T_ap计算ΔT的正演公式

(15)

\begin{array}{l} Δ T = T_{ap} + \frac{f {(\frac{2 π iu}{q_{t_{0}}} F [T_{ap}])}^{2}}{2 T_{0}} + \frac{f {(\frac{2 π iv}{q_{t_{0}}} F [T_{ap}])}^{2}}{2 T_{0}} + \\ \frac{f {(\frac{2 π i (u^{2} + v^{2})^{1 / 2}}{q_{t_{0}}} F [T_{ap}])}^{2}}{2 T_{0}} - \frac{Δ T^{2}}{2 T_{0}} 。 \end{array}

为了计算便捷,避免解方程,式(15)右侧的ΔT可直接采用观测得到的总场异常ΔT_obs。现在根据正演公式可以构建一个求解T_ap的非线性优化问题,目标函数如下

(16)

ϕ = \min_{T_{ap}} ‖ Δ T_{pred} - Δ T_{obs} ‖_{2}^{2} 。

其中:ΔT_pred和ΔT_obs分别表示T_ap正演得到的总场异常和观测得到的总场异常。在不考虑傅里叶变换误差时,该目标函数存在唯一的全局极小值0,因此研究问题可归结为无约束非线性大规模优化问题^{[12,13,14,15]}。求解式(16)最优化问题,通常采用迭代算法,通过迭代更新T_ap最小化ϕ(T_ap)。下文用m表示T_ap,其迭代格式如下

(17)

m_{k + 1} = m_{k} + α_{k} d_{k},

其中:步长α_k由线性搜索策略确定,d_k为搜索方向。对于该优化问题,共轭梯度法与有限内存的BFGS(L-BFGS)^{[14,15,16,17,18]}方法是常见的选择。与流行的共轭梯度法相比,L-BFGS方法作为拟牛顿方法的变体,具有更快的收敛速度,这有助于减少迭代次数^{[19,20,21,22]}。采用L-BFGS方法,搜索方向d_k在L-BFGS方法中可以用Hessian矩阵的逆近似 $H_{k}^{- 1}$ 和梯度g_k表示:

(18)

d_{k} = - H_{k}^{- 1} g_{k} 。

由于文中目标函数的梯度g_k难以直接计算,利用差分代替微分来计算梯度。L-BFGS算法^{[23,24,25,26,27]}无需存储矩阵 $H_{k}^{- 1}$ ,仅存储最近的n步的向量组{s_i,y_i $}_{i = 0}^{k}$ 就能计算 $H_{k}^{- 1}$ ,计算过程如下:

(19)

\begin{array}{l} H_{k + 1}^{- 1} = ({V^{T}}_{k} {V^{T}}_{k - 1} \dots {V^{T}}_{k - n + 1}) H_{k 0}^{- 1} (V_{k - n + 1} \dots {V^{T}}_{k - 1} V_{k}) + \\ ({V^{T}}_{k} {V^{T}}_{k - 1} \dots {V^{T}}_{k - n + 2}) (ρ_{0} s_{0} {s^{T}}_{0}) (V_{k - n + 2} \dots {V^{T}}_{k - 1} V_{k}) + \\ ⋮ \\ {V^{T}}_{k} (ρ_{k - 1} s_{k - 1} {s^{T}}_{k - 1}) V_{k} + ρ_{k} s_{k} {s^{T}}_{k}, \end{array}

其中:V_k=I-ρ_ky_k ${s^{T}}_{k}$ ,ρ_k=1/( ${y^{T}}_{k}$ s_k),y_k=g_k₊₁-g_k,s_k=m_k₊₁-m_k。

2 模型实验

为验证该方法的可靠性,在理论模型上进行验证。采用单球体模型和组合模型,并在此基础上做了弱磁性组合模型、加噪声组合模型、叠加背景场组合模型,模型参数见表1。分别验证该方法在弱磁性、有噪声以及有背景场的数据环境下的应用效果。下文模型实验均采用32×32网格数据,设地磁场强度为50 000 nT,地磁倾角45°偏角30°,采样点线距为100 m(如图3所示)。

表1 模型几何参数和磁化强度

Table 1 Model geometric parameters of and magnetization

模型	中心X坐标 /m	中心Y坐标 /m	半径 /m	中心埋深 /m	磁化强度J /(A·m^-1)
单球体模型	50	50	250	300	120
组合模型A	-300	300	250	350	120
组合模型B	300	-300	400	500	120

新窗口打开| 下载CSV

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 模型水平面与剖面投影

a—单球体模型;b—组合球体模型

Fig.3 Model horizontal and profile projection

a—single sphere model; b—combined sphere model

2.1 单球体模型

如图4所示,磁异常分量T_ap、总场异常ΔT和通过反演得到的T_ap在异常形态上非常相似。但是它们的幅值存在很大差异。计算它们的差异如图5所示。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 单球体模型T_ap、ΔT和计算得到的磁异常分量T_ap对比

a—磁异常在地磁场方向上的投影分量T_ap;b—总场异常ΔT;c—反演得到的磁异常分量T_ap;

Fig.4 Comparison of T_ap, ΔT and calculated T_ap

a—the projection component T_ap of the magnetic anomaly in the direction of the normal geomagnetic field; b—the total field anomaly ΔT; c—the magnetic anomaly component T_ap obtained by the inversion

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 单球体模型反演误差

a—ΔT与T_ap的差值(即初始偏差);b—反演得到的T_ap与理论T_ap存在的偏差

Fig.5 Single sphere model inversion error

a—difference between ΔT and T_ap (ie, initial deviation); b—error between the inverted T_ap and the theoretical T_ap

图4中可以看到ΔT近似误差达到8 000 nT,对后续处理解释会有严重影响。于是使用L-BFGS算法解该问题。利用ΔT初始化T_ap,最终最大误差收敛到20 nT左右。最后反演得到的T_ap与模型T_ap非常接近,表明该方法是有效的。由于目标函数的计算中的包含傅立叶变换,导致结果存在一些误差。但它与原来的近似误差相比已经降低了两个数量级。图6对反演得到的T_ap磁异常分量化极,化极结果与理论模型几乎没有差异,位置和形态都对应非常好。

2.2 组合模型

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 反演得到的磁异常分量T_ap的化极结果

图中红圈表示球体模型投影,黑点表示球体模型的中点

Fig.6 The result of the inversion of the magnetic anomaly component T_ap

The red circle in the figure represents the projection of the sphere model, and the black dot represents the midpoint of the sphere model

如图7所示,T_ap、ΔT和反演计算得到的T_ap也在异常形状上非常相似。但是它们的幅值存在很大差异。计算它们的差异如图8所示。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 组合球体模型T_ap、ΔT和反演得到的磁异常分量T_ap对比

a—磁异常在地磁场方向上的投影分量T_ap;b—总场异常ΔT;c—反演得到的磁异常分量T_ap;

Fig.7 Comparison of T_ap, ΔT and calculated T_ap

a—the projection component T_ap of the magnetic anomaly in the direction of the earth’s magnetic field; b—the total field anomaly ΔT; c—the magnetic anomaly component T_ap obtained by the inversion

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 组合球体模型反演误差

a—ΔT与T_ap的差值(即初始偏差);b—反演得到的T_ap与理论T_ap存在的偏差

Fig.8 Inversion error of combined sphere model

a—difference between ΔT and T_ap (ie, initial deviation); b—error between the inverted T_ap and the theoretical T_ap

由图7可见,组合球体模型正演的ΔT幅值在-12 000 nT到28 000 nT之间,而ΔT与T_ap的最大差值达到8 000 nT,由此可见强磁条件ΔT计算误差很大,不容忽视。与单球体模型类似,利用ΔT初始化T_ap,反演得到的T_ap与模型T_ap的初始最大偏差从 8 000 nT下降到40 nT。由图8可以看出该方法在叠加异常上,迭代结果的误差相较于单球体模型有所增加,说明模型的复杂程度对反演结果存在一定影响。不过原始计算误差还是被降低了两个数量级,可以满足计算要求。图9对反演得到的磁异常分量T_ap化极,化极结果与理论模型几乎没有差异,位置和形态都对应非常好。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 反演得到的磁异常分量T_ap的化极结果

红圈表示球体模型投影,黑点表示球体模型的中点

Fig.9 The RTP result of the magnetic anomaly component T_ap obtained by inversion

The red circle in the figure represents the projection of the sphere model, and the black dot represents the midpoint of the sphere model

2.3 弱磁情况组合模型

前面的两个实验中,强磁条件下,在计算后都还保留一定的误差,如果在弱磁条件下,计算结果的误差会相应降低,还是会固定在数十nT。为此笔者做了弱磁条件下的实验。在弱磁性模型实验中,保持其他参数不变,仅将磁化强度设为12 A/m。图10可见反演计算得到的磁异常分量T_ap与理论值的误差也相应较小,说明傅里叶变换产生的剩余误差同磁异常大小相关,而且该方法在磁异常不同时,都可以将ΔT与T_ap的初始误差降低两个数量级。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 弱磁异常条件下组合球体模型反演得到的T_ap误差

a—ΔT与T_ap的差值(即初始偏差);b—反演得到的T_ap与理论T_ap存在的偏差

Fig.10 error of T_ap obtained by inversion of combined sphere model under weak magnetic anomaly

a—difference between ΔT and T_ap (ie, initial deviation); b—error between the inverted T_ap and the theoretical T_ap

2.4 叠加噪声模型

实际数据往往存在随机噪声,于是本文做了叠加噪声模型,对模型ΔT叠加0~100 nT的随机误差,然后本文方法迭代计算反演T_ap,反演结果图11c仍带有随机噪声,去除随机噪声后得到结果图11d,可以看出该方法基本不受随机噪声影响。

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 叠加噪声的模型实验

a—0~100 nT随机噪声;b—叠加了随机噪声的ΔT的近似误差; c—叠加噪声的ΔT反演得到的T_ap存在的误差; d—反演得到的T_ap消出噪声后的误差结果

Fig.11 Model experiment for superimposed noise

a—0~100 nT random noise; b—approximation error of ΔT superimposed with random noise; c—error of T_ap obtained by ΔT inversion of superimposed noise; d—T_ap error result obtained by inversion (after noise elimination)

2.5 叠加背景场模型

叠加背景场模型,是对模型正演得到的ΔT叠加一个背景场。然后对该叠加场采用本文方法反演计算T_ap,计算结果是仍受背景场影响的T_ap,当去除背景场后得到最终结果图12d,可以看出背景场对迭代计算结果存在一定的影响,不过影响不大。

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 叠加背景场的模型实验

a—背景场;b—叠加了背景场的ΔT的近似误差;c—反演得到的T_ap存在的误差;d— 反演得到的T_ap消除背景场后的误差结果

Fig.12 The model experiment of the superimposed background field

a—the background field; b—approximation error of ΔT (overlay background field); c—error of T_ap obtained by ΔT inversion of superimposed background field; d—T_ap error result obtained by inversion (after eliminating the background field)

3 实测资料计算

将上述方法应用于福建阳山铁矿的实际磁测资料处理中(图13)。研究区内矿体磁性较强,产生的磁异常普遍达到数千nT,最高可以达到10 000 nT以上。模型试验结果显示,由于ΔT处理解释时带有近似误差,化极处理后,不仅会影响强磁性体化极后的幅值位置和形状,还会在强磁异常体附近区域产生较大的虚假异常,如图13c所示。

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 阳山铁矿实测ΔT化极结果与计算得到的T_ap化极结果对比

a—实测资料ΔT化极结果;b—计算得到的T_ap化极结果;c—二者化极结果的差值

Fig.13 Comparison of ΔT reduction to the pole result of measured data and T_ap RTP result obtained by inversion

a—ΔT RTP result of measured data; b—T_ap RTP result obtained by inversion; c—difference between two polarization results

将该方法应用到福建阳山实测资料中,ΔT与计算得到的T_ap化极结果(图13a、b)对比可见,不仅在强磁性体位置处幅值存在较大差异,而且在异常体附近,特别是强磁性体的北侧有很大范围的计算误差,这与模型实验中得到的结论相吻合,该误差在精细化处理解释中不可忽略。剖面化极曲线在曲线形状与极值点位置上存在差异,通过对二者采用组合直立长方体自动反演,得到不同的结果(图14),对比发现见矿钻孔ck40、未见矿钻孔zk4298与ΔT反演结果相互矛盾,与计算得到的T_ap反演结果对应良好。

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 阳山铁矿ΔT与计算得到的T_ap剖面反演对比

a—沿剖面AB,ΔT化极和T_ap化极曲线及其反演拟合结果;b—观测得到的ΔT数据反演结果;c—计算得到的T_ap数据反演结果;见矿钻孔ck40、ck42,未见矿钻孔ck39、zk4298

Fig.14 Profile inversion comparison of ΔT and calculated T_ap

a—ΔT reduction to the pole curve, T_apreduction to the pole curve and their inversion fitting results along the section AB; b—inversion result of the observed ΔT data; c—inversion result of the calculated T_ap data; See ore through dring bore ck40, ck42, no mine ore is found in borehole ck39, zk4298

4 结论

实验证明磁法勘探中将ΔT看作磁异常T_a在地磁场方向上的投影T_ap,并在这种近似下进行处理解释,在具有强磁异常的情况下,如地面磁测时埋深较浅、规模较大的矿体或井中ΔT磁测探头靠近矿体,磁异常非常大甚至超过地磁场,ΔT近似误差E不可忽略。本文提出的利用L-BFGS优化算法反演计算T_ap的方法可以得到十分接近真实T_ap的结果,存在的反演结果的误差主要来源于目标函数中的傅里叶变换。总体上该方法可以将ΔT最大近似误差降低两个数量级,在此基础上进行处理解释可以得到更加精确的结果,实测资料也证明了该结论。随着磁法勘探的发展和计算能力的提高,使用反演得到的T_ap数据进行处理解释,可以很大程度的提高处理解释的精度,符合高精度处理解释的需求。不过该方法的目标函数计算复杂度较高,因此进行大范围T_ap反演时计算量很大,未来在如何提高方法的计算效率方面仍需进一步研究。

The authors have declared that no competing interests exist.

作者已声明无竞争性利益关系。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Parasnis D

Principles of Applied Geop

[M]. Berlin: Springer Netherlands, 1979.