, 在多数磁异常处理过程中, 往往不考虑剩余磁化强度的影响。但在某些条件下, 由于剩余磁化强度的存在, 总场磁化强度的方向与感应磁化强度可能完全不一样。强剩磁和强磁性体的退磁作用改变了磁化强度的大小和方向, 导致磁异常的幅值和形态发生畸变, 影响磁测资料的反演解释[1-3]。若无总场磁化强度方向的先验信息, 常规反演方法计算出来的结果都无效。如:2005年, SHEARER[4]模拟结果表明LI等[5]的反演算法在磁化方向估计误差超过15°时会产生错误的反演结果。
剩磁条件下的磁异常反演方法主要有三大类。
第一类方法:首先估计磁异常数据的磁化方向, 然后直接反演总强度磁异常ΔT或垂直分量Za。1986年, 唐俊德[6]提出利用磁异常三分量确定磁化方向的圆曲线法; 1993年, ROEST等[7]在二维情况下, 提出基于磁场总梯度模和磁源重力异常总水平梯度互相关的估计方法; 1995年, MEDEIROS等[8]提出利用等效源磁矩反演估计磁化方向; 2001年, 甘西[9]利用磁场的垂直分量确定磁性体磁化方向; 2002年, HANEY[10]在二维情况下, 提出一种利用小波变换估计磁化方向的方法; 2004年, BILIM等[11]通过寻找磁源重力异常与重力异常相关系数的极大值来确定磁化方向; 2005年, PHILLIPS[12]给出估计磁化方向的Helbig积分法的直接算法和间接算法; 2006年, NICOLOSI等[13]运用等效源算法估计地壳结构的磁化方向, 同年, DANNEMILLER等[14]提出基于化极磁异常垂直梯度与总梯度模的互相关确定总磁化方向的方法; 2009年, GEROVSKA等[15]通过化极磁异常与磁异常模量的互相关估计磁化方向; 2010年, LI等[16]对常规的估计磁化方向方法进行对比研究, 证明了这些方法不适用于磁性体磁化方向变化较大的情况; 2015年, LIU等[17]在二维条件下, 采用连续反演的方式获得磁性体的物性分布与磁化方向。这类方法一般适用于孤立、剩磁单一的场源体, 只能估计出研究区的平均磁化方向, 往往不适用于研究区磁化方向变化较大的情况。
第二类方法:将原始磁异常转化为受剩余磁化方向影响小、和场源位置对应关系较好的转换量, 然后对转换量进行反演。1972年, NABIGHIAN[18]提出解析信号法(AS), 研究表明二维情况下总梯度模完全独立于磁化方向; 2000年, STAVREV等[19]基于磁异常模量Ta提出其它4种转换量:R模量, E模量, L模量和Q模量; 并证明这些转换量在三维情况弱敏感于磁化方向, 在二维情况下完全独立于磁化方向; 2001年, PAINE等[20]提出两种类似于解析信号(AS)的转换量, 即垂直积分的总梯度模和总梯度模的垂直积分, 并将两种转换量运用于三维条带磁异常反演, 取得比AS更好的反演结果; 2005年, SHEARER[4]进一步提出在三维情况下弱敏感于磁化方向的转换量:磁异常模量和总梯度模, 并提出直接反演这些转换量的三维反演算法; 2006年, STAVREV[21]利用上述5种模量的比率来估计简单二维源的位置及形态; 2010年, LI等[16]进一步补充磁异常模量反演的三维算法, 并提出一套剩磁条件下的三维反演策略, 将其应用于高精度航测磁数据反演并取得较好效果; 2012年, BEIKI等[22]运用同样弱敏感于磁化方向的NSS数据来确定场源体的埋深; 2013年, PILKINGTON等[23]反演了受多个磁化方向干扰地区的NSS数据, 比直接用总场磁异常反演的结果更加可靠; 2014年, LI等[24]反演了受复杂剩磁影响的磁异常模量数据。
第三类方法:直接反演磁化强度矢量。2004年, WANG等[25]首次运用反演总磁化强度矢量的三分量来获得磁化强度的方向, 这种方式更适用于估计孤立、均匀磁化地质体的磁化方向; 2009年, LELIÈVRE等[26]改进了WANG等[25]的方法, 证明这种反演方法能够降低剩余磁化强度的影响。这类方法增强了反演参数的不确定性, 加大了反演难度。
基于上述分析, 特别是转换量NSS具有弱敏感于磁化方向的特性, 本文首先介绍NSS的计算方法并与其它转换量进行了对比, 然后针对磁异常反演的特点引入模型灵敏度矩阵[27], 研究了基于加权模型参数的正则化共轭梯度反演算法, 并应用于强剩磁条件下的NSS反演, 最后利用模型和实测数据验证本文方法的有效性和适用性。
1 归一化磁源强度归一化磁源强度(NSS)是由磁偶极子的磁异常梯度张量矩阵推导而来的转换量[28]。若一个体积为v, 磁化强度分布为M的场源体产生的磁位为:
|
(1) |
式中:|r-r0|为观测点与场源之间的距离。那么, 其磁异常梯度张量Γ为:
|
(2) |
式中:Bx, By和Bz为观测磁异常B(r)的三分量。
矩阵Γ为对称矩阵并且只含有5个独立分量, 因此磁异常梯度张量矩阵Γ可以表示为:
|
(3) |
式中:l1, l2和l3为Γ的特征向量; λ1, λ2和λ3则为其对应的特征值。WILSON[28]给出磁偶极子的归一化磁源强度μ的定义:
|
(4) |
式中:Cm =10-7H/m; m为磁偶极矩。并且推导出NSS可由磁异常梯度张量Γ矩阵的特征值计算出来[22-23, 28]:
|
(5) |
由公式(4) 可以看出, NSS与距离的四次方呈反比, 与偶极子的磁矩呈正比, 与磁化方向无关。NSS具有弱敏感于磁化方向的特性, 当研究区域剩磁不均一时, 该转换类数据相比其它剩磁处理方法更具有优势。
下面通过一个简单块体模型来说明NSS在剩磁条件下的优势。块体模型是一个中心埋深为225 m, 边长为250 m的正方体, 模型的磁化率κ为0.05, 地磁场强度为50 000 nT, 地磁场倾角为90°, 地磁场偏角为0, 磁化倾角为30°, 磁化偏角为60°。图 1a为块体模型的空间位置分布图, 图 1b为块体模型产生的磁异常, 图 1c为磁异常(图 1b)转换出来的总梯度模(AS)数据, 图 1d为磁异常(图 1b)转换出来的磁异常模量(Ta)数据, 图 1e为磁异常(图 1b)转换出来的NSS数据。从图 1c至图 1e可以看出, AS, Ta, NSS都是弱敏感于磁化方向的转换量, 并且NSS与场源体的中心对应性要比同类型的转换量要好。
|
图 1 块体模型的空间位置分布(a)、块体模型产生的磁异常(b)、由磁异常转换得到的总梯度模(c)、由磁异常转换得到的磁异常模量(d)和由磁异常转换得到的归一化磁源强度(e) |
由公式(4) 可知:NSS随距离的四次方衰减, 即NSS对深部场源信息极不敏感, 所以本文引入模型灵敏度矩阵Wm[27]:
|
(6) |
NSS正演核函数A会随着深度的增加而迅速减小(即越深的网格单元的灵敏度越小), 从而导致其反演结果贴近地表而不符合实际情况。通过引入Wm模型灵敏度矩阵可消除核函数随着深度下降造成反演结果贴近地表的影响。因此, 本文采用的模型参数加权为:
|
(7) |
式中:Aw=AWm-1为加权正演算子; mw=Wmm为加权模型参数。模型参数加权以后, 所有网格单元的加权模型参数mw的综合灵敏度一致, 并且对数据d的贡献基本一致。因此反演的目标函数变为:
|
(8) |
式中:φ(mw)为由观测数据与预测数据之差所决定的数据拟合泛函; S(mw)为模型稳定泛函; α为正则化因子;Wd为数据加权矩阵; m0w为加权先验参考模型。最后, 得到共轭梯度的迭代流程[27]如下:
|
(9) |
式中:rnw为观测数据与拟合数据之间的残差; lnα为梯度; βnα为共轭系数; 
实测数据当中往往会含有一定量的噪声, LI等[5]提出用数据的标准偏离差的倒数作为数据加权因子并取得良好效果, 本文也采用这种手段, 故数据加权因子为:
|
(10) |
式中:σi为数据的标准偏离差, 若数据含噪声较少, 则Wd为1。
正则化因子α的选取至关重要, 其选取方法为:首先取α0=0进行第一次迭代试算, 我们选取α1=(‖WdAwm1w-Wdd‖2)/[S(m1w)]作为第一次迭代的正则化因子, 然后按照ZHDANOV[27]提出的自适应正则化算法, 在每次迭代求解出新的模型参数之后, 判断
|
(11) |
如果(11) 式成立, 则表示当前的正则化因子偏大, 无法使数据拟合误差明显下降, 此时, 按照(12) 式来选取下次迭代的正则化因子。
|
(12) |
式中:q可根据实际情况选择。
此外, 我们引入PORTNIAGUINE等[29]的物性约束方式, 即反演结果的物性变化范围满足mbg≤m≤mbg+mup, 其中, mbg为背景物性值, mup为物性上限。
对于该正则化共轭梯度反演流程, 以卡方误差作为迭代停止条件。
1) 当卡方误差达到1时, 停止反演迭代, 其误差的计算方式为:
|
(13) |
式中:dcal为拟合数据; dobs为观测数据; dnoise为观测数据的标准偏离差; N为数据总个数, Erms为卡方误差。
2) 迭代次数达到最大时停止反演。
3 模型试验理论模型由4个大小不同、物性不同、含有不同剩磁的块体组合而成, 表 1列出了各块体的几何参数、物性参数及磁化方向。地磁场强度为50 000 nT, 地磁场倾角为65°, 地磁场偏角为-25°。组合块体模型的空间位置分布如图 2所示, 模型观测磁异常数据(加入数据量级2%的高斯噪声再加1 nT)如图 3a所示, 由磁异常换算得到的NSS如图 3b所示。磁异常数据网格大小为1 m×1 m, 反演网格剖分为1 m×1 m×1 m, 剖分网格个数为21×21×10, 共4 410个; 地面测点在网格中心的正上方, 初始模型m0=0。
|
表 1 组合块体模型各块体的几何参数、物性参数及磁化方向 |
|
图 2 组合块体模型的空间位置分布 |
|
图 3 组合块体模型产生的磁异常(a)和由磁异常转换得到的NSS(b) |
由图 3b可以看出, NSS与场源体有良好的中心对应性。在强剩磁的干扰下, 为了体现NSS反演比常规磁异常反演更具有优势, 我们首先假设磁化方向与地磁场方向一致, 利用本文提出的加权模型参数的三维反演算法对图 3a所示磁总场异常进行反演。
图 4a为磁异常三维反演结果, 可以看出, 磁异常反演结果几乎完全错误。然后, 我们利用基于深度加权的正则化反演算法[5]、基于模型灵敏度的正则化反演算法[27]和本文的反演算法对图 3b所示的NSS进行反演。在深度加权的正则化反演算法中, 由于NSS随距离的四次方衰减, 因此选取的深度加权因子衰减系数为4。为了便于比较, 这3种反演方式都选取最小模型稳定器。图 4b为NSS基于深度加权的正则化方法[5]反演结果, 图 4c为NSS基于模型灵敏度的正则化方法[27]反演结果, 图 4d为NSS加权模型参数反演结果。从图 4中可以看出, 常规反演方法的NSS反演结果比较模糊, 只能大概反映出各个块体模型的位置、埋深及大致轮廓, 异常体整体形态与模型基本相似, 磁化率比真实值偏小; 而本文提出的加权模型参数的反演结果较好, 更清晰地刻画出地下场源体的几何形态与物性分布。由该模型的试验结果可知, 本文的模型参数加权方式比常规的深度加权因子和模型灵敏度矩阵更能降低NSS反演中核函数随距离的四次方衰减对反演结果的影响, 能更清晰地刻画出地下场源体的几何形态与物性分布。
4 实际资料试验
图 5a为加拿大西北部地区让玛丽里弗航测磁异常, 飞机飞行高度为125 m。地磁场倾角为79°, 地磁场偏角为21°。由图 5a可以看出, 图中存在两个磁化方向明显不同的局部磁异常, 采用单一磁化方向对该数据进行常规磁异常反演显然不合适。因此, 我们将磁异常进行换算得到NSS, 如图 5b所示。NSS数据与地下场源体有着良好的中心对应性。磁异常数据网格大小为0.2 km×0.2 km, 反演网格剖分为0.2 km×0.2 km×0.2 km, 剖分网格个数为26×26×10, 共6 760个; 地面测点在网格中心的正上方, 初始模型m0=0。
|
图 5 加拿大西北部地区让玛丽里弗航测磁异常(a)和由磁异常转换得到的NSS数据(b) |
我们分别采用基于深度加权的正则化反演算法和本文反演算法对图 5b所示NSS进行反演, 结果如图 6所示。对比图 6a和图 6b可以看出, NSS加权模型参数的反演结果较好, 能更清晰地刻画出场源体的位置、埋深、大致轮廓及物性分布。从图 6可以看出, 南部场源体中心埋深位于0.5 km处, 从0.2 km延伸至1.0 km; 北部场源体中心埋深位于0.8 km处, 从0.2 km延伸至1.2 km; 两个场源体的整体形态比较相似, 北部场源体磁化率较高。这两个场源体被推断为沿着古河道展布的碎屑磁铁矿[28]。
5 结论
强剩磁的存在改变了磁化方向, 进而影响了磁异常的反演和解释。归一化磁源强度(NSS)是一种由磁异常梯度张量换算而来的弱敏感于磁化方向的转换量, 其与场源的中心对应性比同类型的转换量要好。我们将加权模型参数的正则化共轭梯度反演算法用于NSS反演。模型试验和实际资料反演证明了在强剩磁的条件下, NSS的反演结果比直接反演总场磁异常更加可靠和稳定, 与常规的正则化反演方法相比, 这种模型参数加权的方式能更好地降低NSS反演中核函数随距离的四次方衰减对反演结果的影响, 其反演结果能更好地刻画场源体的几何形态与物性分布。在强剩磁的条件下, 反演NSS数据能够基本确定场源体的位置、埋深、大致轮廓及物性分布, 为磁异常的解释提供了新的思路。
| [1] |
李江霞, 李琴, 李竹强. 青格里底山区块断裂与火成岩重磁异常特征研究[J].
石油物探, 2012, 51(3): 312-318 LI J X, LI Q, LI Z Q. Study on gravity & magnetic anomaly characteristics of fractures and igneous rock in Qinggelidi mountainous exploration block[J]. Geophysical Prospecting for Petroleum, 2012, 51(3): 312-318 |
| [2] |
徐世浙, 张研, 文百红, 等. 切割法在陆东地区磁异常解释中的应用[J].
石油物探, 2006, 45(3): 316-318 XU S Z, ZHANG Y, WEN B H, et al. The application of cutting method to interpretation of magnetic anomaly in region of Ludong[J]. Geophysical Prospecting for Petroleum, 2006, 45(3): 316-318 |
| [3] |
张锡林, 姚长利. 低磁纬度地区磁异常化极方法及应用[J].
石油物探, 2003, 42(4): 541-544 ZHANG X L, YAO C L. Reduction to the pole and its applications in area of low magnetic latitude[J]. Geophysical Prospecting for Petroleum, 2003, 42(4): 541-544 |
| [4] | SHEARER S E.Three-dimensional inversion of magnetic data in the presence of remanent magnetization[D]. Colorado:Colorado School of Mines, 2005 http://geophysics.mines.edu/cgem/pdf%20files/Shearer%20Thesis%202005.pdf |
| [5] | LI Y, OLDENBURG D W. 3-D inversion of magnetic data[J]. Geophysics, 1996, 61(2): 394-408 DOI:10.1190/1.1443968 |
| [6] |
唐俊德. 任意磁性体磁化强度方向的确定[J].
桂林冶金地质学院学报, 1986(2): 75-81 TANG J D. Determination of the direction of magnetization for an arbitrary body[J]. Journal of GuiLin College of Geology, 1986(2): 75-81 |
| [7] | ROEST W R, PILKINGTON M. Identifying remanent magnetization effects in magnetic data[J]. Geophysics, 1993, 58(5): 653-659 DOI:10.1190/1.1443449 |
| [8] | MEDEIROS W E, SILVA J B C. Gravity source moment inversion:a versatile approach to characterize position and 3-D orientation of anomalous bodies[J]. Geophysics, 1995, 60(5): 1342-1353 DOI:10.1190/1.1443870 |
| [9] |
甘西. 利用微机确定磁性体磁化方向的可行性探讨[J].
湖北地矿, 2001, 15(1): 33-38 GAN X. Feasibility of the application of microcomputer method to the magnetic direction determination of the magnetic body[J]. Hubei Geology & Mineral Resources, 2001, 15(1): 33-38 |
| [10] | HANEY M. Total magnetization direction and dip from multi scale edges[J]. SEG Technical Program Expanded Abstracts, 2002, 21(1): 238-244 |
| [11] | BILIM F, ATES A. An enhanced method for estimation of body magnetization direction from pseudogravity and gravity data[J]. Computers & Geosciences, 2004, 30(2): 161-171 |
| [12] | PHILLIPS J D. Can we estimate total magnetization directions from aeromagnetic data using Helbig's integrals[J]. Earth Planets & Space, 2005, 57(8): 681-689 |
| [13] | NICOLOSI I, BLANCO-MONTENEGRO I, PIGNATELLI A, et al. Estimating the magnetization direction of crustal structures by means of an equivalent source algorithm[J]. Physics of the Earth & Planetary Interiors, 2006, 155(1/2): 163-169 |
| [14] | DANNEMILLER N, LI Y. A new method for determination of magnetization direction[J]. Geophysics, 2006, 71(6): L69-L73 DOI:10.1190/1.2356116 |
| [15] | GEROVSKA D, ARA ÚZO-BRAVO M J, STAVREV P. Estimating the magnetization direction of sources from southeast Bulgaria through correlation between reduced-to-the-pole and total magnitude anomalies[J]. Geophysical Prospecting, 2009, 57(4): 491-505 DOI:10.1111/gpr.2009.57.issue-4 |
| [16] | LI Y G, SHEARER S E, HANEY M M, et al. Comprehensive approaches to 3D inversion of magnetic data affected by remanent magnetization[J]. Geophysics, 2010, 75(1): L1-L11 DOI:10.1190/1.3294766 |
| [17] | LIU S, HU X, XI Y, et al. 2D sequential inversion of total magnitude and total magnetic anomaly data affected by remanent magnetization[J]. Geophysics, 2015, 80(3): K1-K12 DOI:10.1190/geo2014-0019.1 |
| [18] | NABIGHIAN M N. The analytic signal of two-dimensional magnetic bodies with polygonal cross-section its properties and use for automated anomaly interpretation[J]. Geophysics, 1972, 37(3): 507-517 DOI:10.1190/1.1440276 |
| [19] | STAVREV P, GEROVSKA D. Magnetic field transforms with low sensitivity to the direction of source magnetization and high centricity[J]. Geophysical Prospecting, 2000, 48(2): 317-340 DOI:10.1046/j.1365-2478.2000.00188.x |
| [20] | PAINE J, HAEDERLE M, FLIS M. Using transformed TMI data to invert for remanently magnetised bodies[J]. Exploration Geophysics, 2001, 32(4): 238-242 DOI:10.1071/EG01238 |
| [21] | STAVREV P. Inversion of elongated magnetic anomalies using magnitude transforms[J]. Geophysical Prospecting, 2006, 54(2): 153-166 DOI:10.1111/gpr.2006.54.issue-2 |
| [22] | BEIKI M, CLARK D A, AUSTIN J R, et al. Estimating source location using normalized magnetic source strength calculated from magnetic gradient tensor data[J]. Geophysics, 2012, 77(6): J23-J37 DOI:10.1190/geo2011-0437.1 |
| [23] | PILKINGTON M, BEIKI M. Mitigating remanent magnetization effects in magnetic data using the normalized source strength[J]. Geophysics, 2013, 78(3): J25-J32 DOI:10.1190/geo2012-0225.1 |
| [24] | LI S L, LI Y G. Inversion of magnetic anomaly on rugged observation surface in the presence of strong remanent magnetization[J]. Geophysics, 2014, 79(2): J11-J19 DOI:10.1190/geo2013-0126.1 |
| [25] | WANG M Y, DI Q Y, KUN X U, et al. magnetization vector inversion equations and forword and inversed 2-d model study[J]. Chinese Journal of Geophysics, 2004, 47(3): 601-609 DOI:10.1002/cjg2.v47.3 |
| [26] | LELIÈVRE P G, OLDENBURG D W. A 3D total magnetization inversion applicable when significant, complicated remanence is present[J]. Geophysics, 2009, 74(3): L21-L30 DOI:10.1190/1.3103249 |
| [27] | ZHDANOV M S. Geophysical inverse theory and regularization problems[M]. The Netherlands: Elsevier Science, 2002: 1-50. |
| [28] | WILSON H S. Analysis of the magnetic gradient tensor[J]. Canada:Defence Research Establishment Pacific, Technical Memorandum, 1985: 5-13 |
| [29] | PORTNIAGUINE O, ZHDANOV M S. Focusing geophysical inversion images[J]. Geophysics, 1999, 64(3): 874-887 DOI:10.1190/1.1444596 |