2. 中国地质科学院地球物理地球化学勘查研究所, 河北廊坊 065000
2. Institute of Geophysical and Geochemical Exploration, Chinese Academy of Geological Sciences, Hebei Langfang 065000, China
磁异常化极是磁法数据处理中的重要方法.处理过程是将磁化强度和磁异常投影的倾角换算成90°,换算结果的意义相当于将异常体放置于磁北极的环境中产生的异常场的垂向分量,目的是使异常场的形态简化,从而增强异常场与异常体的对应关系,有利于推断解释.在定性解释中,甚至可以直接利用磁异常的特征确定构造分布,如断裂位置及其走向、岩体的分布范围、沉积盆地的大致规模等.
地球物理工作者很早就开始重视磁异常化极工作并进行了大量的研究工作.Baranov (1957)提出了将磁异常化到磁北极的概念(简称化极),并采用空间域褶积的方式实现了化极处理.Bhattacharyya (1965)首次在频率域将磁异常频谱与化极转换的频率域滤波因子相乘,得到简捷的频率域化极方法.快速傅里叶变换的出现使得此方法迅速进入真正实用化.由于原理清晰、实现简单、计算快速,频率域化极方法成为化极处理的常用方法.但是随着纬度的降低,磁异常的形态以及与异常体的对应关系会变得非常复杂,化极运算也会变得不稳定,常规化极方法无法有效实现化极转换.磁异常化极因子属于放大性的转换因子,是不稳定算子,在低纬度地区当磁倾角的绝对值趋近于零时,化极因子的放大作用会急剧增强.由于化极因子为复数因子,所以很难直接判断化极算子对磁场波谱的变化过程.为此,人们进行了大量的研究,也取得了不少实用的旨在改善化极效果的频率域低纬度化极方法技术(Macleod et al., 1993;Silva, 1986;Hansen and Pawlowski, 1989;吴健生和王家林,1992;Mendonça and Silva, 1993;Keating and Zerbo, 1996;张培琴和赵群友,1996;Li and Oldenburg, 2001;郭志宏等,2003;Li,2008;骆遥和薛典军, 2009, 2010;方迎尧等,2006;姚长利等, 2003, 2004, 2012;石磊等,2012;林晓星和王平,2012;柴玉璞,2012;骆遥,2013).其中有一类方法是通过对化极因子进行压制从而可以明显改善化极结果,该类化极方法主要包括伪倾角压制法(石磊等,2012)、平方余弦压制法(姚长利等,2003)、维纳滤波法(Hansen and Pawlowski, 1989)、双曲正弦压制法(姚长利等,2004)、直接阻尼法(林晓星和王平,2012)、变频双向阻尼法(柴玉璞,2012).通过对化极算子的分析总结得到:伪倾角压制法和平方余弦压制法是在原化极算子上乘一个函数,实施直接压制;维纳滤波法、直接阻尼法和变频双向阻尼法是在原化极算子的分母中加一个函数(大于零),实施间接压制(或称阻尼);双曲正弦压制法则是在原化极算子的分母中加一个函数(大于零),并从其分子中减一个函数(大于零),同时实施直接压制和间接压制.尽管三种方法压制的具体方式不同,但上述方法均会对磁场频谱中相位谱进行改动.骆遥(2013)在实数空间利用Hartley变换进行磁异常化极,对改造后的化极算子进行谱分析,结果表明滤波作用在压制化极算子幅度的同时也在一定程度上改变了相位谱.
根据频谱理论可知,频谱是由振幅谱和相位谱组成,振幅谱具有“因子相乘”特性,相位谱具有“因子相加”的特性.因此如果将化极因子和磁场频谱均分解为振幅谱和相位谱,在滤波过程中由于相位谱只进行加减运算而不会存在有极值的问题,只需要考虑化极因子的振幅谱在极端情况的应对措施,并可以参照向下延拓算子和垂向导数算子一类的滤波器进行压制或阻尼处理.本文在原有的低纬度化极因子处理方法的基础上,提出了对原有化极过程及化极因子的改进方法,从而改善低纬度磁场的化极效果,主要讨论对阻尼因子方法的改进,并通过模型计算说明改进阻尼法的计算效果.
2 化极原理及分析磁异常化极计算在空间域为复杂的褶积计算,而在频率域为简单的乘积运算.所以在频率域中的转换计算只需要乘上相应的转换因子,然后反变换到空间域就可得到转换后的空间域结果.当有多个转换因子的时候,可以在频率域内一次性完成频谱转换,这是频率域处理转换的优势.频率域中的频谱是由振幅谱和相位谱组成,滤波过程可以分解为振幅谱相乘和相位谱相加,公式表述为:
(1) |
式中:Q为磁场频谱,AQ为磁场振幅谱函数,φQ为磁场相位谱函数,H为化极因子频谱,AH为化极因子振幅谱函数,φH为化极因子相位谱函数.
频率域中的化极因子H(u, v)的表达式为:
(2) |
式中:qk=
将式(2)中的qk分解为振幅谱和相位谱可得:
(3) |
式中:
当磁场方向和磁化方向一致时I2=I3,I0=I1,因为化磁极时I2=I3=
(4) |
式中:
式(4)化极因子中的振幅谱函数和相位谱函数是I0和θ的函数.相位谱函数φH在I0和θ-D0为任何值的情况下都是有限且连续,振幅谱AH在I0和θ-D0同时趋近于0和±
(5) |
当θ→D0±
图 1表明化极因子中的相位谱值位于[-π, π]区间内,地磁倾角I0的值越小,曲线变化越剧烈.当I0=0的极端情况时曲线变为正负均分的方波形状,φH的函数值为:
(6) |
式(6)中φH的绝对值除个别点为零值外其他值都是π,即化极因子的相位谱基本上是将原始异常的相位谱改变180°.化极过程中化极因子中的相位谱与磁场的相位谱只进行加减运算,所以这部分计算是稳定的.图 2表明化极因子振幅谱的滤波作用主要集中于垂直主磁场方向的扇形区域内,化极的纬度越低对波谱的调制作用越大,扇形区变得越尖锐,会导致相应区域的计算过程不稳定,出现扇形死亡地带.因此化极因子在低纬度化极过程只需要对振幅谱进行处理即可.
3 阻尼化极方法
姚长利等(2004)提出了一种直接阻尼的低纬度化极方法,能在低纬度地区取得良好的效果.通过对化极因子进行改造可以防止化极数据在“死亡带”内出现不稳定的现象.其原理是在化极因子的分母上直接加上一个很小的值,这样就可以使化极因子在D0±
(7) |
式中:α=θ-D0,d是一个很小的量,α0为一较小的角度.Iα为低纬度特征角,表示I0小于该角度就采取针对性的低纬度处理措施.该参量在变倾角地区化极进行分块处理,会使得结果进行拼接很自然.带阻尼因子的化极因子如下:
(8) |
通过上述对化极因子的改造,可以使磁异常化极在低纬度地区变得稳定.该方法在理论研究和实际应用中均有重要意义,但是该阻尼因子同时会改变化极因子的相位谱.通过第2节分析,化极因子中的相位谱不应该进行改动,仅改动振幅谱.因此本文提出只需要将阻尼因子加入到化极因子的振幅谱中即可,最终的化极因子公式由式(4)变为:
(9) |
则化极因子中的振幅谱加入阻尼因子后的特征如图 3,与图 2对比可以看出,加入阻尼因子后化极因子的振幅谱函数变为有限且连续.
模型为一长方体(图 5a),长、宽分别为20 m,厚度为2 m,上顶埋深为1 m,计算64×64的网格总场异常,网格纵、横间距都是1 m.图 5b为垂直磁化条件下的理论磁异常,即化极异常.模型参数设置为磁倾角I0=0°,磁偏角D0=30°.为准确说明问题,化极处理过程中未进行扩边处理.参量d的大小直接关系到阻尼因子的阻尼程度,d越大阻尼作用越强.为了对比阻尼化极方法和改进的阻尼化极方法的效果,同时也为了选取最佳的参量d,有必要分析化极结果与参量d的关系.本文选取参量d为自变量,化极结果与理论化极结果差值的离差值为变量进行分析(图 4),在磁倾角I0=0°的极端情况下,该对比结果可以在一定程度上反映两种方法的稳定性和准确性.图 4中的阻尼法的误差曲线随着阻尼作用增加整体趋势变小,但是曲线呈锯齿状,说明阻尼法相对不稳定;改进的阻尼法误差曲线也是随着阻尼作用的增加整体趋势变小,曲线整体光滑无锯齿出现,曲线值一般不大于阻尼法的曲线值,并且在d=0.04附近出现极小值.
根据图 4曲线的结果选取了参量d的两个数值0.04和0.4进行对比分析.图 5c和图 5d分别为d=0.04时阻尼法和改进阻尼法的化极结果,通过对比化极结果不难看出,改进阻尼法的化极效果有明显的改善,异常形态和异常大小均与理论值更为接近;图 5e和图 5f分别为d=0.4时阻尼法和改进阻尼法的化极结果,两者的结果基本一致,改进的阻尼法的效果略好,但等值线形态均不如d=0.04的等值线效果好.试验结果表明改进的阻尼方法更加稳定且准确度更高,而且可以通过上述误差曲线图为参量d的选取提供一定的参考依据.
4.2 多体组合模型计算效果实际地质情况一般为多个异常体的组合,理论上多体组合模型计算与单体模型计算的性质是相同的,但是组合异常体的磁异常相互叠加,异常变化比单体模型复杂.本文设计了多个异常体的组合模型以模拟磁赤道地区的情况,来检验方法的化极效果.基本参数选取磁倾角I0=0°,磁偏角D0=0°,产生的水平磁化ΔT异常等值线见图 6a,理论化极结果见图 6b,数据大小为128×64.采用余弦阻尼法和改进的阻尼法分别进行化极计算,结果见图 7a和图 7b所示.通过对比化极结果不难发现,改进的阻尼法在多体模型中的化极应用效果也有明显的改善.
理论公式表明改进的阻尼方法能更好地保留化极因子中的信息,并可以减少阻尼因子的干扰.数值模拟实验表明,改进阻尼算法对低纬度的磁场数据能取得较好的化极效果.控制参量与阻尼方法一致,意义明确.本文为其他压制类方法化极效果的改进提供了可行的思路.
致谢感谢论文评审过程中专家提出的宝贵意见.
Baranov V. 1957. A new method for interpretation of aeromagnetic maps:pseudo-gravimetric anomalies. Geophysics, 22(2): 359-383. DOI:10.1190/1.1438369 | |
Bhattacharyya B K. 1965. Two-dimensional harmonic analysis as a tool for magnetic interpretation. Geophysics, 30(5): 829-857. DOI:10.1190/1.1439658 | |
Cai Y P. 2012. Discussion of wavenumber domain RTP methods at low latitudes. Oil Geophysical Prospecting (in Chinese), 47(3): 496-505. | |
Fang Y Y, Zhang P Q, Liu H J. 2006. Approaches to the interpretation of magnetic ΔT anomalies in the low magnetic latitude area. Geophysical & Geochemical Exploration (in Chinese), 30(1): 48-54. | |
Guo Z H, Yu C C, Zhou J X. 2003. The tangent technique of ΔT profile magnetic anomaly in the low magnetic latitude area. Geophysical & Geochemical Exploration (in Chinese), 27(5): 391-394. | |
Hansen R O, Pawlowski R S. 1989. Reduction to the pole at low latitudes by Wiener filtering. Geophysics, 54(12): 1607-1613. DOI:10.1190/1.1442628 | |
Keating P, Zerbo L. 1996. An improved technique for reduction to the pole at low latitudes. Geophysics, 61(1): 131-137. DOI:10.1190/1.1443933 | |
Li X. 2008. Magnetic reduction-to-the-pole at low latitudes:observations and considerations. The Leading Edge, 27(8): 990-1002. DOI:10.1190/1.2967550 | |
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 | |
Lin X X, Wang P. 2012. An improved method for reduction to the pole of magnetic field at low latitude-the method of frequency conversion bidirectional damping factor. Chinese Journal of Geophysics (in Chinese), 55(10): 3477-3484. DOI:10.6038/j.issn.0001-5733.2012.10.031 | |
Luo Y, Xue D J. 2009. Stable reduction to the pole at low magnetic latitude by probability tomography. Chinese Journal of Geophysics (in Chinese), 52(7): 1907-1914. DOI:10.3969/j.issn.0001-5733.2009.07.026 | |
Luo Y, Xue D J. 2010. Reduction to the pole at the geomagnetic equator. Chinese Journal of Geophysics (in Chinese), 53(12): 2998-3004. DOI:10.3969/j.issn.0001-5733.2010.12.024 | |
Luo Y. 2013. Hartley transform for reduction to the pole. Chinese Journal of Geophysics (in Chinese), 56(9): 3163-3172. DOI:10.6038/cjg20130929 | |
MacLeod I N, Jones K, Dai T F. 1993. 3-D analytic signal in the interpretation of total magnetic field data at low magnetic latitudes. Exploration Geophysics, 24(3-4): 679-688. | |
Mendonça C A, Silva J B C. 1993. A stable truncated series approximation of the reduction-to-the-pole operator. Geophysics, 58(8): 1084-1090. DOI:10.1190/1.1443492 | |
Shi L, Guo L H, Meng X H, et al. 2012. The modified pseudo inclination method for magnetic reduction to the pole at low latitudes. Chinese Journal of Geophysics (in Chinese), 55(5): 1775-1783. DOI:10.6038/j.issn.0001-5733.2012.05.035 | |
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 | |
Wu J S, Wang J L. 1992. Improving the effect of reducing magnetic anomaly to pole in low magnetic latitude area by using directional high cut filter. Oil Geophysical Prospecting (in Chinese), 27(5): 670-677. | |
Yao C L, Guan Z N, Gao D Z, et al. 2003. Reduction-to-the-pole of magnetic anomalies at low latitude with suppression filter. Chinese Journal of Geophysics (in Chinese), 46(5): 690-696. | |
Yao C L, Huang W N, Zhang Y W, et al. 2004. Reduction-to-the-pole at low latitude by direct damper filtering. Oil Geophysical Prospecting (in Chinese), 39(5): 600-606. | |
Yao C L, Li H W, Zheng Y M, et al. 2012. Research on iteration method using in potential field transformations. Chinese Journal of Geophysics (in Chinese), 55(6): 2062-2078. DOI:10.6038/j.issn.0001-5733.2012.06.028 | |
Zhang P Q, Zhao Q Y. 1996. Methods of the magnetic direction transform of aeromagnetic anomalies with differential inclinations in low magnetic latitudes. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 18(3): 206-214. | |
柴玉璞. 2012. 从化极算法误差方程看各种波数域低纬度化极方法. 石油地球物理勘探, 47(3): 496–505. | |
方迎尧, 张培琴, 刘浩军. 2006. 低磁纬度地区ΔT异常解释的途径与方法. 物探与化探, 30(1): 48–54. | |
郭志宏, 于长春, 周坚鑫. 2003. 低磁纬度区ΔT剖面磁异常场源深度计算的切线法. 物探与化探, 27(5): 391–394. | |
林晓星, 王平. 2012. 一种改进的低纬度磁场化极方法--变频双向阻尼因子法. 地球物理学报, 55(10): 3477–3484. DOI:10.6038/j.issn.0001-5733.2012.10.031 | |
骆遥, 薛典军. 2009. 基于概率成像技术的低纬度磁异常化极方法. 地球物理学报, 52(7): 1907–1914. DOI:10.3969/j.issn.0001-5733.2009.07.026 | |
骆遥, 薛典军. 2010. 磁赤道处化极方法. 地球物理学报, 53(12): 2998–3004. DOI:10.3969/j.issn.0001-5733.2010.12.024 | |
骆遥. 2013. Hartley变换化极. 地球物理学报, 56(9): 3163–3172. DOI:10.6038/cjg20130929 | |
石磊, 郭良辉, 孟小红, 等. 2012. 低纬度磁异常化极的伪倾角方法改进. 地球物理学报, 55(5): 1775–1783. DOI:10.6038/j.issn.0001-5733.2012.05.035 | |
吴健生, 王家林. 1992. 用高阻方向滤波器提高低磁纬度地区磁异常化极效果. 石油地球物理勘探, 27(5): 670–677. | |
姚长利, 管志宁, 高德章, 等. 2003. 低纬度磁异常化极方法--压制因子法. 地球物理学报, 46(5): 690–696. | |
姚长利, 黄卫宁, 张聿文, 等. 2004. 直接阻尼法低纬度磁异常化极技术. 石油地球物理勘探, 39(5): 600–606. | |
姚长利, 李宏伟, 郑元满, 等. 2012. 重磁位场转换计算中迭代法的综合分析与研究. 地球物理学报, 55(6): 2062–2078. DOI:10.6038/j.issn.0001-5733.2012.06.028 | |
张培琴, 赵群友. 1996. 低磁纬度区航磁异常变倾角磁方向转换方法. 物探化探计算技术, 18(3): 206–214. | |