2. 国土资源部航空地球物理与遥感地质重点实验室, 北京 100083
2. Key Laboratory of Airborne Geophysics and Remote Sensing Geology, Ministry of Land and Resources, Beijing 100083, China
磁异常形态不仅反映构造格架,也能在纵向上揭示地壳的深部信息,地质构造无论出露于地表还是隐伏于地下其演化形迹多能在航磁图中清晰反映.但不同于重力异常较简单的对应关系,磁异常不仅取决于磁性场源自身,还取决于现今地磁场的强度和方向,即使不考虑剩磁,磁异常也较同源重力异常复杂许多.为便于解释航磁异常,1957年Baranov[1-2]提出化极概念,以消除倾斜磁化影响,简化异常解释.随后Bhattacharyya[3]通过傅里叶级数实现化极处理.Kanasewich和Agarwal[4]建议用快速傅里叶变换实现Bhattacharyya算法.Gunn[5]则给出了统一的频率域位场转换表达.目前,化极处理基本上均通过快速傅里叶变换在频率域中实现,是一种重要的磁异常转换处理方法.
除难以解决剩磁问题外,化极处理在低纬度地区也面临困难.频率域化极算子属于放大类型的滤波器,当磁纬度较低特别是磁赤道处,频率域化极算子在垂直于地磁偏角方向变得极不稳定甚至奇异,表现出沿主磁场方向的条带状干扰.Baranov[1]建议化极中磁倾角应大于30°,随后的算法将其限制为16.5°[2].Silva[6]认为倾角小于15°的频率域化极不具实用性.在解决低纬度甚至磁赤道处化极困难时,最直接办法是用化赤代替化极[7],但水平磁化资料解释难度大,化赤处理应用极少.Silva[6]采用等效源方式实现低纬度化极处理.Hansen和Pawlowski[8]提出用Wiener滤波解决低纬度化极问题,Keating和Zerbo[9]曾对该方法进行改进.Li和Oldenburg [10-11]采取反演依次实现了空间域和频率域的化极处理.在国内,张培琴等[12]根据补偿圆滑滤波思想提出余弦削顶方向滤波和双曲函数滤波方法处理低纬度化极.姚长利等[13-14]针对频率域化极算子中扇形放大区域分别提出了压制因子和阻尼因子,实现低纬度化极.骆遥等[15]利用概率成像进行等效反演实现低纬度化极,此外提出在空间域迭代实现磁赤道处化极[16].处理低纬度磁异常时,解析信号理论一直受到重视,例如MacLeod等[17]曾对低纬度化极与总梯度模解释方法进行对比,但需要指出的是三维解析信号理论并不完善,不同于二度体磁异常梯度模与磁化方向无关,三度体磁异常的总梯度模仍受到磁化方向影响.
在处理大面积磁测资料时,用某一固定的磁化方向进行化极显然并不合适,区域内如果主磁场方向变化较大,化极结果将出现畸变.Arkani-Hamed[18]提出用等效层形式进行变倾角化极(DRTP),Swain[19]针对低纬度问题对其进行了改进,后续Arkani-Hamed[20]又提出对自身DRTP算法的改进.Cooper和Cowan[21]采用泰勒级数展开进行DRTP处理,但一阶以上的化极场导数的含义尚不明确,此外该方法本身仍面临低纬度问题.国内针对大面积磁测资料化极时,多采用分带、分块变倾角化极或滑动窗口变倾角化极方法[12].这类处理本质上与通常的化极没有区别,只是涉及数据的拼接或编图问题,但分带、分块处理中可能损失磁异常的长波信息.
国内外学者对化极特别是低纬度化极进行过大量研究,对理论模型化极能得到较精确的结果,实际资料处理中能够获得满意的化极结果,但化极处理特别是主流的频率域化极在理论上仍存在某些问题.例如频率域中化极算子存在扇形放大区域,低纬度时尤为明显[22],但化极算子如何作用于波谱并不能像向下延拓算子等相位不变的滤波器那样易于说明,只能给出化极算子的振幅或相位[7].当然,某些问题源自傅里叶变换本身,为此Chai[23]提出偏移抽样理论以期提高傅里叶变换的数值精度,并应用于磁倾角为3°的低纬度模型进行化极试验[24].本文从Hartley变换性质出发,提出用Hartley变换实现化极,分析了频率域化极算子并给出化极过程中波谱的具体变化.同时,结合现有低纬度化极方法实现Hartley变换的低纬度化极,理论模型计算和实际资料处理表明该方法具有实用性.
2 Hartley变换化极 2.1 Hartley变换及性质1942年Hartley [25]提出的一种类似傅里叶变换的积分变换,即Hartley变换.1983年傅里叶分析专家Bracewell[26]给出了离散Hartley变换(DHT).随后,不少学者提出了Hartley变换的快速算法[27-28],即快速Hartley变换(FHT).文献表明DHT所需的内存和计算时间要比傅里叶变换节省近一半,有大量关于Hartley变换研究及应用的成果发表在IEEE出版的杂志上.二维Hartley变换则存在着两种不同的形式,其区别在于积分核,Sundararajan[29]曾对此进行了讨论.根据Brancewell[30]的建议,本文采用标准形式的二维Hartley变换.二维Hartley变换为[31]
(1) |
(2) |
其离散形式为
(3) |
(4) |
其中cas (θ)=cosθ+sinθ.可以看出,如果将DHT定义中的M-1N-1因子换成M-1/2N-1/2,其逆变换也相应乘以M-1/2N-1/2,就可以得到一种严格对称的变换对.实际上,同傅里叶变换一样,上述因子的具体位置是无关紧要的,可以证明任何信号经两次Hartley变换均能完全恢复.由于DHT变换均在实数空间中进行,一张图像就可以表示Hartley变换的频谱,而傅里叶变换后的谱则需要分别用实部和虚部表述,或是借助幅度和相位两张图像表示,正因如此相位谱在位场资料处理中几乎从未真正使用,多数情况仅是利用功率谱分析幅度.Hartley变换后其频谱的独特优势,则有望为后续谱分析提供行之有效的手段.
为了推导Hartley变换化极的表达式,首先给出Hartley变换中位场转换的某些表达.同傅里叶变换类似,Hartley变换中存在如下微分关系:
(5) |
其中HT表示Hartley变换,f(x,y)的Hartley谱为H(u,v),u、v分别是x、y方向上的波数.此外,对于位场资料特有的解析延拓和垂向导数表达,有:
(6) |
(7) |
其中,h<0,
ΔT可以近似为磁异常矢量在主磁场方向的投影,其定义为:
(8) |
根据重磁位泊松公式:
(9) |
可以将ΔT用同源重力位U的二阶梯度张量表示,有:
(10) |
其中,Hax、Hay、Za分别是磁异常矢量的三个分量,B0={Bx,By,Bz}是主磁场方向的单位矢量(不考虑剩磁情况),μ0是真空中的磁导率,M是磁化强度,G是万有引力常数,ρ是同源地质体的密度.
对式(10)进行Hartley变换,有:
(11) |
其中,HT(u,v)是ΔT(x,y)的Hartley谱,H(u,v)是同源重力位U的Hartley谱,u、v分别是x、y方向上的波数,
令u=-u、v=-v,则:
(12) |
根据式(11)和(12),可以得到:
(13) |
当垂直磁化时,有:
(14) |
其Hartley谱为
(15) |
因此,Hartley域的化极谱可表述为
(16) |
其中,
当水平磁化即Bz=0时,有:
(17) |
此时因传统的频率域化极算子中不存在复数而与Hartley域化极算子在形式上完全一致,Hartley变换化极同样面临低纬度困难.为了进一步说明低纬度问题,将式(16)的化极谱写成径向形式,有:
(18) |
其中,HPa=
上述Hartley域化极算子对谱幅度的改变分析同样适用于傅里叶变换,因为Hartley变换与傅里叶变换存在密切的联系.例如Hartley变换谱与傅里叶变换的实部和虚部之差等价,此外利用Hartley变换也可以实现快速傅里叶变换.正因如此,目前的一些低纬度化极处理如压制因子或阻尼因子等方法主要对图 1a的进行了处理,没有精细的考虑图 1b中波谱放大的情况.当然,这是由于傅里叶变换本身造成的问题,姚长利等[13-14]建议频谱压制的扇形区间大小为12°,实质上正是综合考虑上述两种情况中对功率谱综合放大的一种折衷结果.实现低纬度Hartley变换化极的精确方法在于分别考虑对HPa和HPb两种滤波算子的压制改造,而不是以往综合地考虑在一个扇形波数域内进行压制.实现低纬度Hartley域化极中对HPa和HPb两种化极算子的改造需要考虑下列三个问题:(1)对化极算子的改造不易过大,否则将失去化极谱特征;(2)改造后的化极算子仍要充分光滑;(3)对HPa和HPb两种化极算子的改造程度应尽可能相同.彭利丽等[34]进一步对姚长利提出的阻尼因子进行改进,使其不仅压制谱的实部,同时也改造虚部,实际上是在考虑问题(3).综合考虑上述三种要求,选取下面的滤波算子对Hartley变换化极算子进行改造,即
(19) |
这里称
为检验Hartley变换化极方法,选取Hansen和Pawlowski[8]提出的理论模型进行化极处理.模型长宽均为20m,高2 m,顶面埋深1 m的磁性长方体(图 3a),计算64×64大小的磁异常网格,网格横纵间距均为1m.图 3b给出了垂直磁化条件下的理论磁异常(化极异常).此处分别采用了不同斜磁化条件下ΔT进行检验,为准确说明问题,化极处理中未对模型数据进行包括数据扩充在内的任何处理. 图 4给出了Hartley变换化极的处理效果,限于篇幅仅展示磁倾角为60°、磁偏角为30°(图 4a)和磁倾角为15°、磁偏角为30°(图 4c)两组模型的化极效果,分别代表典型的中高纬度和低纬度地区.可以看出,图 4b和图 4d给出化极结果与图 3b的理论化极场无论是在形态上还是数值上都非常一致.需要指出的是,该模型埋藏极浅且数据量很小,位场转换的实际误差较大,图 4c给出的斜磁化条件属于低纬度,在未进行包括扩边、滤波等特殊处理条件下能给出图 4d这一化极效果,本身即能够说明化极处理的准确性.
检验低纬度化极效果时,则选取了磁赤道处理论模型进行化极检验,图 5a给出了水平磁化条件下的理论异常场.采用式(19)的滤波形式进行低纬度处理,其中I′=20°,图 5b给出化极后的效果.可以看出低纬度磁赤道处化极效果明显,一定程度上由于滤波原因,部分异常未能完全恢复,其效果不如专门提出的磁赤道化极方法[16]理想,但对比使用余弦削顶方向滤波或双曲函数滤波对同样数据进行的化极处理,几者间的化极效果是较一致的.对比文献[34]对类似模型化极进行的化极处理对比,也可以确认图 5b的化极效果.
为了进一步说明低纬度化极中不同滤波参数的影响,此处选取了两种较极端的情况对图 5a的磁赤道磁异常模型进行化极处理,图 6a和图 6b分别为采用I′=3°和I′=60°的化极效果.图 6a表明较接近主磁场方向的视磁化方向能够给出更为接近理论模型的化极结果,此时化极算子的放大作用较大,出现了轻微的南北向条带状干扰-这表明视磁化倾角很小时理论模型仍能取得较好的化极效果.图 6b表明当视磁化倾角较大时化极场更趋于稳定,化极异常与原始异常趋近,近似的互为倒相.磁异常倒相解释方法目前已应用于低纬度磁测解释,例如郭志宏等[35]采用通过ΔT剖面负磁异常的“反切”做法,使得计算中高磁纬度区磁异常场源深度的切线法及系数表可直接用于低磁纬度地区,一定程度上解决了低磁纬度地区用切线法计算场源深度的问题.方迎尧等[36]通过对大量实测的航磁资料研究,分析了不同磁化条件下ΔT、Za、Ha的关系,提出了低磁纬度地区ΔT异常倒相180°的解释方法,用于南海低纬度地区航磁资料解释.上述大角度视磁化倾角的磁赤道化极表明,磁异常倒相解释方法本身就是一种滤波压制较大的化极处理,这也揭示了骆遥等[16]提出磁赤道处化极迭代中采用拟合异常差值的“倒相”这一校正方案的合理性.用倒相方法处理低纬度磁异常或化赤异常具有其自身合理性,一定程度上简化了低纬度磁异常解释.
由于前面的理论模型计算已经表明Hartley变换的化极的效果,限于篇幅这里不再展示Hartley变换在中、高纬度地区实际资料化极的结果,仅给出一个低纬度化极实例进行说明.利用Hartley变换对南海某区域磁异常数据进行化极处理,该区内的磁倾角非常接近0°属磁赤道地区.该区域由跨度不大,其主磁场方向的变化很小,这里直接采用式(19)的滤波方式进行Harley变换化极,视磁化倾角取20°.图 7a和图 7b分别为实测的总场磁异常和化极结果.化极后磁异常仍保持稳定,没有出现畸变,异常的形态和幅度发生了显著的改观.该资料已采取过多种化极方法进行处理[15-16],虽然以往该资料的化极处理更侧重于反算原始磁异常的拟合程度,但就位场转换的结果而言,图 7b化极效果也较为合理.实际资料处理表明,通过结合现有的一些低纬度化极的滤波处理方式即可直接实现Hartley变换进行低纬度化极处理.
采用Hartley变换建立了Hartley变换化极方法,并结合现有低纬度化极方法实现了利用Hartley变换的低纬度化极.通过对Hartley变换化极算子分析,揭示了化极处理中具体的波谱变化,为进一步解决低纬度化极理论困难提供了条件.文献表明快速Hartley变换所需的内存和计算时间要比傅里叶变换节省近一半,建立的Hartley域化极方法在实际资料处理中将更具优势,能够实现更大规模的资料处理,同时也为逐点化极的实用化提供了更为高效的技术手段.理论模型计算和实际资料处理表明建立的Hartley域化极方法是准确可靠的,能够在磁测资料处理中发挥积极作用.
附录A 位场二维Hartley变换的解析延拓算子根据外部Dirichlet问题的公式,在整个h<0上半空间,位场有
(A1) |
(A1)可以写成褶积形式,有
(A2) |
由于-h/(x2 +y2 +h2)3/2是偶函数,式(A2)的Hartley变换满足类似傅里叶变换中的褶积定理,其Hartley变换为
(A3) |
其中HT表示Hartley变换.借助位场转换中一个最基本的积分公式
(A4) |
其中z>0,则-h/(x2+y2+h2)3/2的Hartley变换有
(A5) |
即
(A6) |
因此,不论采用傅里叶变换还是Hartley变换,其位场解析延拓的频率域算子具有相同的形式,频率域解析延拓的步骤和过程也非常类似.
致谢感谢复审专家指正公式中存在的错误.
[1] | Baranov V. A new method for interpretation of aeromagnetic maps: Pseudo-gravimetric anomalies. Geophysics , 1957, 22(2): 359-383. DOI:10.1190/1.1438369 |
[2] | Baranov V, Naudy H. Numerical calculation of the formula of reduction to the magnetic pole. Geophysics , 1964, 29(1): 67-79. DOI:10.1190/1.1439334 |
[3] | Bhattacharyya B K. Two-dimensional harmonic analysis as a tool for magnetic interpretation. Geophysics , 1965, 30(5): 829-857. DOI:10.1190/1.1439658 |
[4] | Kanasewich E R, Agarwal R G. Analysis of combined gravity and magnetic fields in wave number domain. Journal of Geophysical Research , 1970, 75(29): 5702-5712. DOI:10.1029/JB075i029p05702 |
[5] | Gunn P J. Linear transformations of gravity and magnetic fields. Geophysical Prospecting , 1975, 23(2): 300-312. DOI:10.1111/gpr.1975.23.issue-2 |
[6] | Silva J B C. Reduction to the pole as an inverse problem and its application to low-latitude anomalies. Geophysics , 1986, 51(2): 369-382. DOI:10.1190/1.1442096 |
[7] | Blakely R J. Potential Theory in Gravity and Magnetic Applications. New York: Cambridge University Press, 1996 . |
[8] | Hansen R O, Pawlowski R S. Reduction to the pole at low latitude by Wiener filtering. Geophysics , 1989, 54(12): 1607-1613. DOI:10.1190/1.1442628 |
[9] | Keating P, Zerbo L. An improved technique for reduction to the pole at low latitudes. Geophysics , 1996, 61(1): 131-137. DOI:10.1190/1.1443933 |
[10] | Li Y, Oldeuburg D W. Reduction to the pole using equivalent sources. 60th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2000: 386-389. |
[11] | Li Y G, Oldeuburg D W. Stable reduction to the pole at the magnetic equator. Geophysics , 2001, 66(2): 571-578. DOI:10.1190/1.1444948 |
[12] | 张培琴, 赵群友. 低纬度区航磁异常变倾角磁化方向转换方法. 物探化探计算技 , 1996, 18(3): 206–214. Zhang P Q, Zhao Q Y. Methods of the magnetic direction transform of aeromagnetic anomalies with differential inclinations in low magnetic latitudes. Computing Techniques for Geophys. Geochem. Expl. (in Chinese) , 1996, 18(3): 206-214. |
[13] | 姚长利, 管志宁, 高德章, 等. 低纬度磁异常化极方法-压制因子法. 地球物理学报 , 2003, 46(5): 690–696. Yao C L, Guan Z N, Gao D Z, et al. Reduction to the pole of magnetic anomalies at low latitude with suppression filter. Chinese J. Geophys. (in Chinese) , 2003, 46(5): 690-696. |
[14] | 姚长利, 黄卫宁, 张聿文, 等. 直接阻尼法低纬度磁异常化极技术. 石油地球物理勘探 , 2004, 39(5): 600–606. Yao C L, Huang W N, Zhang Y W, et al. Reduction-to-the-pole at low latitude by direct damper filtering. Oil Geophys. Prosp. (in Chinese) , 2004, 39(5): 600-606. |
[15] | 骆遥, 薛典军. 基于概率成像技术的低纬度磁异常化极方法. 地球物理学报 , 2009, 52(7): 1907–1914. Luo Y, Xue D J. Stable reduction to the pole at low magnetic latitude by probability tomography. Chinese J. Geophys. (in Chinese) , 2009, 52(7): 1907-1914. |
[16] | 骆遥, 薛典军. 磁赤道处化极方法. 地球物理学报 , 2010, 53(12): 2998–3004. Luo Y, Xue D J. Reduction to the pole at the geomagnetic equator. Chinese J. Geophys. (in Chinese) , 2010, 53(12): 2998-3004. |
[17] | MacLeod I, Jones K, Dai T F. 3-D analytic signal in the interpretation of total magnetic field data at low magnetic latitudes. Expl. Geophys. , 1993, 24(4): 679-688. DOI:10.1071/EG993679 |
[18] | Arkani-Hamed J. Differential reduction-to-the-pole of regional magnetic anomalies. Geophysics , 1988, 53(12): 1592-1600. DOI:10.1190/1.1442441 |
[19] | Swain C J. Reduction-to-the-pole of regional magnetic data with variable field direction, and its stabilisation at low inclinations. Expl. Geophys. , 2000, 31(2): 78-83. DOI:10.1071/EG00078 |
[20] | Arkani-Hamed J. Differential reduction to the pole: revisited. Geophysics , 2007, 72(1): L13-L20. DOI:10.1190/1.2399370 |
[21] | Cooper G R J, Cowan D R. Differential reduction to the pole. Computers & Geosciences , 2005, 31(8): 989-999. |
[22] | Li X. Magnetic reduction-to-the-pole at low latitudes: observations and considerations. The Leading Edge , 2008, 27(8): 990-1002. DOI:10.1190/1.2967550 |
[23] | Chai Y P. Shift sampling theory of Fourier transform computation. Sci. China (Ser E) , Sci. 1997, 40(1): 21-27. DOI:10.1007/BF02916587 |
[24] | Chai Y P. A-E equation of potential field transformations in the wavenumber domain and its application. Appl. Geophys. , 2009, 6(3): 205-216. DOI:10.1007/s11770-009-0032-z |
[25] | Hartley R V L. A more symmetrical Fourier analysis applied to transmission problems. Proc. IRE. , 1942, 30(3): 144-150. DOI:10.1109/JRPROC.1942.234333 |
[26] | Bracewell R N. Discrete Hartley transform. J. Opt. Soc. Am. , 1983, 73(12): 1832-1835. DOI:10.1364/JOSA.73.001832 |
[27] | Bracewell R N. The fast Hartley transform. Proc. IEEE. , 1984, 72(8): 1010-1018. DOI:10.1109/PROC.1984.12968 |
[28] | Hou H S. The fast Hartley transform algorithm. IEEE Trans. Computers. , 1987, 36(2): 147-156. |
[29] | Sundararajan N. 2-D Hartley transforms. Geophysics , 1995, 60(1): 262-267. DOI:10.1190/1.1443754 |
[30] | Bracewell R N. Fourier Transform and Its Applications. (3rd edition). New York: McGraw-Hill, 2000 . |
[31] | Bracewell R N. The Hartley Transform. New York: Oxford University Press, 1986 . |
[32] | 魏雅利, 骆遥. 基于Hartley变换的剖面位场转换. 地球物理学进展 , 2010, 25(6): 2102–2108. Wei Y L, Luo Y. 2D potential field transformation based on Hartley transform. Progress in Geophys. (in Chinese) , 2010, 25(6): 2102-2108. |
[33] | Mendonca C A, Silva J B C. A stable truncated series approximation of the reduction-to-the-pole operator. Geophysics , 1993, 58(8): 1084-1090. DOI:10.1190/1.1443492 |
[34] | 彭利丽, 郝天珧, 姚长利, 等. 低纬度磁异常化极方法应用效果对比. 地球物理学进展 , 2010, 25(1): 151–161. Peng L L, Hao T Y, Yao C L, et al. Comparison of the application effects of the reduction-to-the-pole methods at low magnetic latitudes. Progress in Geophys. (in Chinese) , 2010, 25(1): 151-161. |
[35] | 郭志宏, 于长春, 周坚鑫. 低磁纬度区ΔT剖面磁异常场源深度计算的切线法. 物探与化探 , 2003, 27(5): 391–394. Guo Z H, Yu C C, Zhou J X. The tangent technique of ΔT profile magnetic anomaly in the low magnetic latitude area. Geophysical & Geochemical Exploration (in Chinese) , 2003, 27(5): 391-394. |
[36] | 方迎尧, 张培琴, 刘浩军. 低磁纬度地区ΔT异常解释的途径与方法. 物探与化探 , 2006, 30(1): 48–54. Fang Y Y, Zhang P Q, Liu H J. Approaches to the interpretation of magnetic ΔT anomalies in the low magnetic latitude area. Geophysical & Geochemical Exploration (in Chinese) , 2006, 30(1): 48-54. |