Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 

Magnetic source parameter estimation based on variable depth DEXP method of the local wavenumber

  • XiaoMei LIU , 1, 2 ,
  • YanGuo WANG , 1, 2, *
Expand
  • 1 National Key Laboratory of Uranium Resources Exploration-Mining and Nuclear Remote Sensing, Nanchang 330013, China
  • 2 School of Geophysics and Measurement Control Technology, East China University of Technology, Nanchang 330013, China

Received date: 2024-01-06

  Online published: 2025-01-14

Copyright

Copyright ©2024 Progress in Geophysics. All rights reserved.

Abstract

Local wavenumber is a common method for magnetic data processing and interpretation. DEXP method can quickly image the spatial position of magnetic source, but the spatial imaging resolution is poor, easy to be affected by superposition ies, and easy to produce false information. In this paper, we propose to improve the local wave number DEXP method, and construct the variable depth imaging function by introducing the depth scale factor. This method can use the maximum value to reflect the field source space position, and estimate the magnetic source structure index through the extreme value and the inversion depth value. At the same time, combining different depth scale factors and different orders of imaging results can effectively improve the reliability and applicability of the method. Model tests show that the improved method has stronger spatial imaging resolution and more accurate estimates of field source parameters relative to the conventional local wavenumber DEXP method. Finally, the new method was applied to the ground magnetic anomaly in the youth point area of Nenbei Farm in Heilongjiang Province, and obtained good application results.

Cite this article

XiaoMei LIU , YanGuo WANG . Magnetic source parameter estimation based on variable depth DEXP method of the local wavenumber[J]. Progress in Geophysics, 2024 , 39(6) : 2337 -2344 . DOI: 10.6038/pg2024HH0501

0 引言

磁法勘探是矿产勘查最为常用的地球物理方法之一,磁法数据处理是地质解释的重要依据.由于局部波数在二维空间不受磁化方向影响,因此无需化极处理,而备受关注与使用,基于局部波数的磁源参数估计方法可以快速完成场源参数的计算(Thurston and Smith, 1997; Smith et al., 1998),是一类磁法常用的自动解释工具.Smith和Salem(2005)采用局部波数的各种组合来估计磁源深度和形状(构造指数);Salem等(2005)利用原始局部波数和相位转换后局部波数进行异常体的反演,并介绍了三维应用效果;Keating(2009)利用局部波数及其垂直导数对异常体的深度及构造指数进行估计,但是该方法会增大噪声的干扰;崔莉和王万银(2011)认为向上延拓某一高度来进行局部波数反演,可以解决高频干扰,从而提高局部波数的反演精度;李丽丽等(2012)采用不同位置和不同高度磁异常的局部波数直接计算场源体的深度和构造指数,不需要计算局部波数的导数,且降低了噪声干扰;马国庆等(2012)提出了采用张量局部波数法来进行位场全张量数据的解释,通过理论模型证明张量局部波数法可以很好地完成位场全张量数据的反演工作;Ma和Li(2013)利用不同水平位置与不同高度局部波数组合进行异常的解释,获得了良好的效果;马国庆等(2014)提出归一化局部波数法进行重磁异常的解释,并采用了多种归一化函数试算.然而,这类方法由于使用了不同阶次的导数,因此对噪声较为敏感.
Abbas和Fedi(2014)提出了基于局部波数DEXP法(Depth From Extreme Points),该方法利用成像极大值来识别场源空间位置,采用极值的数值及深度值来估算构造指数,由于该方法使用了稳定的向上延拓计算,因此结果具有一定的稳定性,但该方法成像分辨率较低而更易受叠加异常的影响.针对于此,本文借鉴田野和王彦国(2023)提出的变深度成像磁源估计方法,通过引入一个深度尺度因子,提出了磁源变深度局部波数DEXP成像方法.并通过模型试验和实例应用验证了本文方法的可行性、有效性与实用性,同时与常规局部波数DEXP法进行了对比分析.

1 理论基础

Miller和Singh(1994)给出了一种可以均衡不同异常的倾斜角法,其二维空间表达式为:
$Tilt=\arctan \frac{T_z}{T_x}, $
其中TxTz是磁异常Txz方向上的导数.
一阶局部波数k1被定义为解析信号局部相位的变化率(Thurston and Smith, 1997):
$\begin{aligned}k_1 & =\frac{\partial}{\partial x} \arctan \left(\frac{T_z}{T_x}\right) \\& =\frac{1}{|A(x)|^2}\left(\frac{\partial^2 T}{\partial x \partial z} \cdot \frac{\partial T}{\partial x}-\frac{\partial^2 T}{\partial x^2} \cdot \frac{\partial T}{\partial x}\right), \end{aligned}$
其中解析信号表示为:
$A=\sqrt{\left(\frac{\partial T}{\partial x}\right)^2+\left(\frac{\partial T}{\partial z}\right)^2}.$
Smith等(1998)通过研究,发现局部波数与场源的水平位置x0、深度z0及构造指数N有关,即:
$k_1(x, z)=\frac{(N+1)\left(z_0-z\right)}{\left[\left(x-x_0\right)^2+\left(z-z_0\right)^2\right]} \text {, }$
式中(x, z)是观测位置,(x0, z0)是场源位置,其中构造指数N与场源几何形状有关(Thompson, 1982),N=0,1,2分别对应台阶,岩脉和水平圆柱体.
二阶局部波数可以写成(Smith et al., 1998):
$k_2=\frac{\partial}{\partial x} \arctan \left(\frac{T_{z z}}{T_{x z}}\right)=\frac{(N+2)\left(z_0-z\right)}{\left[\left(x-x_0\right)^2+\left(z-z_0\right)^2\right]} \text {, }$
Abbas等(2014)给出了n阶局部波数表达式:
$\begin{aligned}k_n & =\frac{\partial}{\partial x} \arctan \left(\frac{\partial^n T}{\partial z^n} / \frac{\partial^n T}{\partial x \partial z^{n-1}}\right) \\& =\frac{(N+n)\left(z_0-z\right)}{\left[\left(x-x_0\right)^2+\left(z-z_0\right)^2\right]}, \end{aligned}$
从公式(6)可以看出,局部波数表达式与场源位置(x0, z0)和构造指数N有关,因此,可以构建一种函数用于估计场源参数.这里,定义一个基于局部波数的变深度成像函数LWDInβ
$\operatorname{LWDI}_n^\beta=(-z)^\beta \cdot k_n=\frac{(N+n)\left(z-z_0\right)(-z)^\beta}{\left[\left(x-x_0\right)^2+\left(z-z_0\right)^2\right]}, $
其中,β是深度尺度因子,z≤0,为向上延拓高度,函数LWDInβ单位为mβ-1.
公式(7)对xz方向取导数,令其等于0,则:
$\left\{\begin{array}{l}\frac{\partial \operatorname{LWDI}_n^\beta}{\partial x}=\frac{2(N+n)\left(z-z_0\right)\left(x-x_0\right)(-z)^\beta}{\left[\left(x-x_0\right)^2+\left(z-z_0\right)^2\right]^2}=0 \\\frac{\partial \operatorname{LWDI}_n^\beta}{\partial z}=\frac{(N+n)\left\{\left[(\beta+1) z-\beta z_0\right]\left(x-x_0\right)^2-\left[\beta z_0+(1-\beta) z\right]\left(z-z_0\right)^2\right\}}{\left[\left(x-x_0\right)^2+\left(z-z_0\right)^2\right]^2}=0\end{array}, \right.$
求解公式(8)得:
$x=x_0, z=\frac{\beta}{\beta-1} z_0, $
由于z≤0,z0≥0,因此,0 < β < 1,公式(9)表明,函数LWDInβ在坐标$ \left(x_0, \frac{\beta}{\beta-1} z_0\right)$处取得极大值.因此,可以利用该函数极值识别出场源的空间位置.
当场源位置(x0, z0)确定后,还可以利用公式(7)计算场源构造指数N
$N=\frac{\operatorname{LWDI}_n^\beta\left(x_0, z_0\right) \cdot z_0^{1-\beta}}{\beta^\beta(1-\beta)^{1-\beta}}-n, $
虽然该方法使用了高阶导数,但同时也运用了稳定的向上延拓,因此结果具有一定稳定性.另外,选取较小的深度尺度因子β时,成像分辨率较高,但易受噪声干扰影响;而选取较大的β时,结果稳定性较强,但空间成像分辨率较低.还需指出的是,当β=0.5时,该方法与常规局部波数DEXP法(Abbas et al, 2014)完全相同.
由于函数LWDInβ是非线性函数,因此成像结果中可能存在虚假极值,不过这些虚假信息往往表现为偶极特征(同一水平位置,相近高度上存在极大值与极小值),与真实场源呈现的单极特征明显不一致,因此易于识别.
为了方便使用本文方法,给出了方法的计算流程:
(1) 选择深度尺度因子β和延拓高度间隔Δz,则缩放后的高度间隔为$ \frac{\beta}{1-\beta} \Delta z$
(2)计算出不同延拓高度磁异常导数:$ \frac{\partial^n T}{\partial z^n}, $$ \frac{\partial^n T}{\partial x \partial z^{n-1}}, \frac{\partial^{n+1} T}{\partial x^n \partial z}, \frac{\partial^{n+1} T}{\partial z^{n+1}}$,获取kn
(3) 计算函数LWDInβ,并对高度轴按$ \frac{\partial^n T}{\partial z^n}$比例回放,根据极大值拾取场源位置;
(4) 利用成像极大值数值和深度值,采用公式(10)计算场源构造指数.

2 模型试验

2.1 单一模型

为了测试本文方法的可行性与有效性,建立了一个宽为1 m,上顶埋深为15 m的岩脉,磁化倾角为45°.图 1给出了该岩脉在地面产生的理论磁异常及含40%随机噪声的磁异常.
图1 单一岩脉磁异常及40%噪声磁异常

Fig 1 Magnetic anomaly and 40% noise-added magnetic anomaly

图 2是无噪声时深度尺度因子β=0.3、0.5时的1阶局部波数DEXP成像结果,可以看出,两个成像结果图中均存在一个明显的极大值,位于(100 m, 15 m)处,极大值位置与岩脉上顶埋深位置完全一致.但β=0.3时成像结果更加聚焦,β=0.5时(常规方法)的成像结果较为发散,说明改进方法对模型体位置反映更加明显.
图2 无噪声时不同深度尺度因子β下的1阶局部波数DEXP成像结果

(a)β=0.3;(b)β=0.5.

Fig 2 The first-order local wavenumber DEXP image applied to noise-free magnetic anomaly

(a)β=0.3;(b)β=0.5.

图 3是含40%随机噪声时深度尺度因子β=0.5、0.6时的1阶局部波数DEXP成像结果,可以看出,β=0.5时(常规方法)成像结果图已不能识别出场源位置,但β=0.6时的1阶成像结果图可以清晰准确地识别出场源位置.
图3 含噪声时不同深度尺度因子β下的1阶局部波数DEXP成像结果

(a)β=0.5;(b)β=0.6.

Fig 3 The first-order local wavenumber DEXP image applied to noise-corrupted magnetic anomaly

(a)β=0.5;(b)β=0.6.

图 2ab图 3b的极值分别为0.1634 m-0.7,0.2591 m-0.5和0.3422 m-0.4,利用公式(10)计算得到的构造指数分别为1.00、1.00和0.93,与岩脉的理论构造指数1基本完全一致.
该试验表明,选择较小的深度尺度因子时,成像分辨率更高,而选择较大的深度尺度因子时,抗干扰能力更强,因此灵活选择深度尺度因子可以提高方法的有效性及适用性.

2.2 叠加模型

为了验证新方法在复杂情况下的应用效果,建立了由1个台阶、2个岩脉及2个水平圆柱体组成的叠加模型,模型参数见表 1图 4是叠加模型的磁异常.
表1 叠加模型参数

Table 1 Overlays the model parameters

模型体编号 模型体类型 边长或半径/m 上顶或质心埋深/m 下底/m 有效磁化倾角/(°) 磁化强度/(A·m-1)
台阶(N=0) 30 60 0.1
岩脉(N=1) 0.5 20 45 10
岩脉(N=1) 0.5 30 60 10
水平圆柱体(N=2) 5 25 45 10
水平圆柱体(N=2) 5 50 60 10
图4 叠加模型磁异常

Fig 4 Magnetic anomalies generated by multi-source model

图 5是不同β的1阶局部波数DEXP成像结果,当β=0.01时(图 5a),DEXP成像图能够较好地反映出所有场源的空间位置,极值分别位于(101 m, 30 m)、(300 m, 20 m)、(500 m, 30 m)、(699 m, 25 m)和(900 m, 61 m),显然除了埋深最大的第二个圆柱体深度存在明显偏差外,其他场源上的识别结果与理论值完全一致,利用公式(10)计算得到的场源构造指数分别是0、1.1、1.1、2.1和3.2,同样只有埋深最大的场源构造指数误差最大;当β=0.2时(图 5b),DEXP成像图包含了4个极大值,分别位于(101 m, 33 m)、(300 m, 22 m)、(500 m, 35 m)及(698 m, 28 m),依次对应着台阶、两个岩脉及第一个圆柱体的空间位置,根据深度估计值和极值数值,利用公式(10)计算得到的构造指数分别为0.2、1.2、1.4和2.4,此时不仅识别的场源深度值和构造指数估计值与理论值存在一定偏差,而且无法识别出第二圆柱体;当β=0.5时(图 5c),DEXP成像图已无法反映出任一场源.
图5 叠加模型磁异常不同β的1阶局部波数DEXP成像结果

(a)β=0.01;(b)β=0.2;(c)β=0.5.

Fig 5 The first-order local wavenumber DEXP image applied to the synthetic magnetic anomaly

(a)β=0.01;(b)β=0.2;(c)β=0.5.

图 6是不同β的2阶局部波数DEXP成像结果,可以看出,β=0.1和0.2时的局部波数DEXP成像图均可以较好地反映出所有场源的空间位置,β=0.1时的5个场源深度估计值分别是30、20、30、25、49 m,与理论深度值基本完全吻合,对应的构造指数估计值依次是0、1、1、2和1.9,同样与理论值较吻合;β=0.2时成像图反映的场源深度分别位于30、20、29、25、49 m,对应的构造指数估计值分别是0、1、0.9、2和1.9,显然场源参数估计值与理论值非常接近;当β=0.5时的成像图仅能有效地识别出第一个岩脉的位置(300 m, 21 m)和估计出构造指数N=1.1,无法反映出其他四个模型体的空间位置.
图6 叠加模型磁异常不同β的2阶局部波数DEXP成像结果

(a)β=0.1;(b)β=0.2;(c)β=0.5.

Fig 6 The second-order local wavenumber DEXP image applied to the synthetic magnetic anomaly

(a)β=0.1;(b)β=0.2;(c)β=0.5.

该叠加模型试验表明,常规1阶局部波数DEXP法(β=0.5)在该叠加模型试验中已失效,而2阶局部波数DEXP法也仅能识别出1个场源位置,再次显示了本文方法具有更强的空间成像分辨率和更高的场源空间位置识别精度.

3 实例应用

为了验证本文方法在实际数据中的应用效果,选取了黑龙江嫩北农场青年点地区的地面磁异常数据进行测试.研究区整体位于弱磁性的二长花岗岩分布,但目前发现该区域及紧邻区域存在褐铁矿化蚀变带和闪长岩脉.图 7是研究区1:1万磁异常,其中点距为20 m,可以看出,磁异常场整体呈现西低东高的异常特征,工区内存在数条明显的条带磁异常,其中中部存在2条近NS向条带磁异常已被证实为褐铁矿化蚀变带引起的,其他NE向的条带磁异常推测是由闪长岩脉引起的.这里选取了图 7中的AA′剖面进行测试,该剖面基本垂直于所有构造走向,对应的磁异常及一、二阶解析信号分别见图 8图 9图 10, 从解析信号图中可以清晰地看出该剖面存在5个极大值.图 11图 12分别是β=0.2、0.3、0.5的1、2阶局部波数DEXP成像结果,表 2给出了不同情况下的磁源参数估计结果.从图 11图 12可以看出,当β=0.2和0.3时的1阶局部波数DEXP可以获得5个有效磁源位置,而β=0.5时的1阶局部波数DEXP和任意β时2阶局部波数DEXP均只能获得4个磁源位置. 从表 2可以看出,不同β的1、2阶局部波数DEXP法获得磁源空间位置和构造指数基本一致,5个磁源分别位于(100 m, 40 m)、(340 m, 50 m)、(760 m, 40 m)、(1320 m, 40 m)和(2580 m, 40 m)处,其中第1、4、5个磁源的构造指数接近1,表明这三个磁源接近岩脉,而第2、3个磁源的构造指数却明显大于1,说明了这2个场源具有一定的宽度,恰好这2个场源对应的是NS走向的褐铁矿化蚀变带,而其他3个近似脉状的磁源对应磁异常是NE向,推测为3条闪长岩脉.
图7 黑龙江嫩北农场青年点地区磁异常图

Fig 7 Magnetic anomaly of Qingniandian area in Nenbei Farm, Heilongjiang Province

图8 AA′剖面磁异常

Fig 8 Magnetic anomaly of the profile AA′

图9 AA′剖面磁异常1阶解析信号

Fig 9 The first-order analytic signals of the profile AA′

图10 AA′剖面磁异常2阶解析信号

Fig 10 The second-order analytic signals of the profile AA′

图11 AA′剖面不同β的1阶局部波数DEXP成像结果

(a)β=0.2;(b)β=0.3;(c)β=0.5.

Fig 11 The first-order DEXP image with different β applied to the magnetic anomaly of the profile AA′

(a)β=0.2;(b)β=0.3;(c)β=0.5.

图12 AA′剖面不同β的2阶局部波数DEXP成像结果

(a)β=0.2;(b)β=0.3;(c)β=0.5.

Fig 12 The second-order DEXP image with different β applied to the magnetic anomaly of the profile AA′

(a)β=0.2;(b)β=0.3;(c)β=0.5.

表2 AA′剖面不同β不同阶次局部波数DEXP法的磁源参数估计表

Table 2 Magnetic source parameters of AA′ profile estimated by local wavenumber DEXP method with different β

场源 阶次 β 场源位置/m EV N β 场源位置/m EV N β 场源位置/m EV N
1 1 0.2 (100, 40) 0.06 1.01 0.5 (100, 40) 0.08 0.99 0.3 (100, 30) 0.17 0.88
2 (100, 40) 0.09 1.04 (100, 40) 0.13 1.10 (100, 40) 0.24 1.10
2 1 (340, 50) 0.07 1.89 (340, 50) 0.10 1.89
2
3 1 (760, 40) 0.07 1.51 (760, 40) 0.10 1.56 (760, 30) 0.19 1.51
2 (760, 50) 0.10 2.01 (760, 50) 0.15 2.10 (760, 40) 0.29 1.69
4 1 (1320, 40) 0.06 0.98 (1320, 40) 0.08 0.99 (1320, 40) 0.16 1.02
2 (1320, 50) 0.09 1.48 (1320, 40) 0.12 1.45 (1320, 40) 0.23 0.97
5 1 (2580, 40) 0.06 1.02 (2560, 50) 0.08 0.98 (2560, 50) 0.16 1.29
2 (2580, 40) 0.09 1.04 (2580, 40) 0.12 1.06 (2580, 40) 0.23 0.99
该应用实例表明,改进的局部波数DEXP法能够识别出更多的磁源,而且利用不同阶次、不同β的成像结果可以相互印证,增加结果的可靠性,提高方法的实用性.

4 结论

本文在分析不同阶次局部波数表达式与磁源参数关系的基础上,借鉴局部波数DEXP法,提出了变深度局部波数DEXP法,该方法通过引入一个深度尺度因子β,使得方法具有更强的灵活性,还可以结合不同β值的成像结果来提高方法的可靠性.模型试验结果表明,选用较小的β可以获得更高的空间成像分辨率和场源位置识别精度,而选用较大的β可以获得更稳定的计算结果,通过灵活选择深度尺度因子β,可以显著提高局部波数DEXP方法的场源识别能力.在实例应用中,提供了6种不同情况的局部波数DEXP成像结果,使得本文方法在磁源参数估计方面更加可靠.

感谢审稿专家提出的修改意见和编辑部的大力支持!

Abbas M A Fedi M. Automatic DEXP imaging of potential fields independent of the structural index Geophysical Journal International 2014 199 3 1625 1632

DOI

Abbas M A Fedi M Florio G. Improving the local wavenumber method by automatic DEXP transformation Journal of Applied Geophysics 2014 111 250 255

DOI

Cui L Wang W Y. The technique for application of local wavenumber to magnetic interpretation Geophysical & Geochemical Exploration 2011 35 6 779 784

Keating P. Improved use of the local wavenumber in potential-field interpretation Geophysics 2009 74 6 L75 L85

DOI

Li L L Du X J Ma G Q. Improved local wavenumber methods in the interpretation of magnetic fields Journal of Jilin University (Earth Science Edition) 2012 42 4 1179 1185

DOI

Ma G Q Du X J Li L L. Comparison of the tensor local wavenumber method with the conventional local wavenumber method for interpretation of total tensor data of potential fields Chinese Journal of Geophysics 2012 55 7 2450 2461

DOI

Ma G Q Huang D N Li L L . A normalized local wavenumber method for interpretation of gravity and magnetic anomalies Chinese Journal of Geophysics 2014 57 4 1300 1309

DOI

Ma G Q Li L L. Alternative local wavenumber methods to estimate magnetic source parameters Exploration Geophysics 2013 44 4 264 271

DOI

Miller H G Singh V. Potential field tilt—A new concept for location of potential field sources Journal of Applied Geophysics 1994 32 2-3 213 217

DOI

Salem A Ravat D Smith R . Interpretation of magnetic data using an enhanced local wavenumber (ELW) method Geophysics 2005 70 2 L7 L12

DOI

Smith R S Salem A. Imaging depth, structure, and susceptibility from magnetic data: The advanced source-parameter imaging method Geophysics 2005 70 4 L31 L38

DOI

Smith R S Thurston J B Dai T F . iSPITM—the improved source parameter imaging method Geophysical Prospecting 1998 46 2 141 151

DOI

Thompson D T. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data Geophysics 1982 47 1 31 37

DOI

Thurston J B Smith R S. Automatic conversion of magnetic data to depth, dip, and susceptibility contrast using the SPI (TM) method Geophysics 1997 62 3 807 813

DOI

Tian Y Wang Y G. Estimating magnetic source parameter using variable depth imaging method based on analytic signal amplitude Progress in Geophysics 2023 38 5 2135 2146

DOI

万银. 局部波数在磁异常解释应用中的方法技术 物探与化探 2011 35 6 779 784

丽丽 晓娟 国庆. 改进的局部波数法及其在磁场数据解释中的应用 吉林大学学报(地球科学版) 2012 42 4 1179 1185

DOI

国庆 晓娟 丽丽. 解释位场全张量数据的张量局部波数法及其与常规局部波数法的比较 地球物理学报 2012 55 7 2450 2461

DOI

国庆 大年 丽丽 . 重磁异常解释的归一化局部波数法 地球物理学报 2014 57 4 1300 1309

DOI

彦国. 基于解析信号变深度成像的磁源参数估计方法 地球物理学进展 2023 38 5 2135 2146

DOI

Outlines

/