重、磁方法具有成本低、周期短、覆盖面广等特点，已成为解决空间探测、大地构造、区域构造，固体矿产勘探、能源勘探、地震及地质灾害预测、环境及工程检测等方面问题不可或缺的重要手段.目前，位场数据处理和转换分别在空间域和频率域中进行，频率域中的位场数据处理和转换主要是基于Fourier变换(侯重初和刘秀芳，1978；Nagendra R, 1984；吴宣志等，1987；朱文孝等，1989；穆石敏等，1990)、余弦变换(张凤旭等, 2006, 2007)、Hilbert变换(Nabighian，1972；Stanley and Green, 1976；Mohan et al., 1982；张凤旭等，2005；骆遥等，2011)、小波变换(李宗杰等，1997；Li and Oldenburg, 1998)以及Hartley变换(Rao et al., 1996；魏雅利和骆遥，2010；骆遥，2013；马国庆等，2014)来实现的.由于频率域位场数据处理和转换有计算公式简单、计算效率高的特点，是位场数据处理和转换的研究重点.Hartley变换是电气工程师Hartley(1942)提出的一种类似于Fourier变换的实数域积分变换方法.著名的Fourier分析专家Bracewell(1983)首次提出了离散Hartley变换(DHT)的定义.随后不少学者(Bracewell，1984；Hou，1987)提出了Hartley变换的快速算法(FHT).理论上，DHT所需的内存和计算时间要比复数域DFT(离散Fourier变换)节省近一半(Bracewell，1986)，已有部分关于Hartley变换研究及应用的文章发表在IEEE所属期刊上，其广泛应用及显示出的重要作用也为地球物理学家所重视.自1990年开始，Saatcilar等(1990)最早将Hartley变换引入到地球物理领域.同年，Marobhe I M应用Hartley变换解释二维板状体磁异常，这是国外最早将Hartley变换方法应用于重、磁数据处理和解释中.Saatcilar和Ergintav(1991)利用Hartley变换求解弹性波方程，使得Hartley变换在地震数据处理领域得到快速发展，与传统的复数域Fourier变换方法相比，其计算量大大减少，计算时间也比复数域Fourier变换减少近一半.随后，国外学者纷纷将Hartley变换应用于重磁位场数据处理和转换中(Sundararajan and Brahmam, 1998；Kadirov，2000；Al-Garni Mansour and Sundararajan, 2012)，而我国学者自20世纪90年代开始研究Hartley变换，并且最早将其应用于地震波场模拟以及偏移成像等方面工作(刘迎曦等，1993；缪林昌，1993；尧德中等，1994；周辉和何樵登，1995；张文生，2003).直到魏雅利和骆遥(2010)第一次将Hartley变换引入到剖面位场数据处理和转换中，推导了Hartley变换位场解析延拓及垂向n阶导数因子，从而填补了国内在利用Hartley变换研究位场处理和转换方面的空白.随后骆遥(2013)将Hartley变换应用于磁异常低纬度化极，为低纬度化极提供了新的方法技术.马国庆等人(2014)研究了Hartley变换在位场导数计算方面的应用，通过理论模型试验，验证了Hartley变换导数计算结果与Fourier变换导数计算结果接近.而在位场数据处理和转换中，Hartley变换和Fourier变换的优劣研究却很少.

本文通过整理Hartley变换和Fourier变换的定义、性质、计算量以及在位场数据处理和转换中的频率响应，并通过理论模型测试和实际资料处理试验来研究Hartley变换和Fourier变换在位场数据处理和转换中的优劣.

(1) Hartley变换和Fourier变换的定义

设一维实函数f(x)的Hartley变换及其逆变换定义为(Hartley，1942)

(1)

其中cas(2πux)=cos(2πux)+sin(2πux).

对于二维Hartley变换的定义，存在两种不同的形式，其区别在于积分核不同.Sundararajan(1995)曾对此进行讨论，根据Bracewell(2000)的建议，本文采用标准形式的二维Hartley变换，则二维连续实函数f(x, y)的Hartley变换及其逆变换定义为

(2)

其中cas[2π(ux+vy)]=cos[2π(ux+vy)]+sin2π(ux+vy)].

可以看出，Hartley变换及其逆变换具有完全相同的形式，构成了严格对称的变换，并且Hartley变换在实数域空间进行.

而对于一维函数f(x)的Fourier变换及其逆变换定义为(布隆什坦和谢缅佳也夫，1965)

(3)

对于二维函数f(x, y)的Fourier变换及其逆变换定义为

(4)

从Fourier变换和Hartley变换的定义形式对比可以看出，二者定义形式基本相似.Hartley正、反变换定义式具有完全相同的核函数，且为实数；而Fourier正、反变换定义式中的核函数互为共轭，且为复数.

(2) Hartley变换和Fourier变换的性质

Hartley变换具有与Fourier变换相似的特征，也拥有Fourier变换相似的大部分性质(甘利灯，1992).下面给出一维和二维Hartley变换的性质.

设一维连续实函数f(x)、f₁(x)和f₂(x)，其Hartley变换为H(u)=HT[f(x)]，H₁(u)=HT[f₁(x)]，H₂(u)=HT[f₂(x)]，其Fourier变换为F(u)=FT[f(x)]，F₁(u)=FT[f₁(x)]，F₂(u)=FT[f₂(x)].其中u是x方向的频率，“HT”代表Hartley变换算子，“FT”代表Fourier变换算子，表 1给出了一维Hartley变换和Fourier变换的性质(甘利灯，1992).

表 1 (Table 1)

表 1 一维Hartley变换与Fourier变换性质对比表 Table 1 The contrast table of properties of 1D Hartley transform and Fourier transform

表 1 一维Hartley变换与Fourier变换性质对比表 Table 1 The contrast table of properties of 1D Hartley transform and Fourier transform

设二维连续实函数f(x, y)、f₁(x, y)和f₂(x, y)，其Hartley变换谱分别为H(u, v)=HT[f(x, y)]，H₁(u, v)=HT[f₁(x, y)]，H₂(u, v)=HT[f₂(x, y)]，其Fourier变换谱分别为F(u, v)=FT[f(x, y)]，F₁(u, v)=FT[f₁(x, y)]，F₂(u, v)=FT[f₂(x, y)]，其中，u、v分别是x和y方向的频率.“HT”作为Hartley变换算子，“FT”作为Fourier变换算子，表 2给出了二维Hartley变换与Fourier变换的性质(甘利灯，1992).

表 2 (Table 2)

表 2 二维Hartley变换与Fourier变换性质对比表 Table 2 The contrast table of properties of 2D Hartley transform and Fourier transform

表 2 二维Hartley变换与Fourier变换性质对比表 Table 2 The contrast table of properties of 2D Hartley transform and Fourier transform

从表 1和表 2可以看出，Hartley变换的时移定理、频移定理、时域微分定理和褶积定理均需将原来的频谱做关于零频的对称谱，然后对其频谱和对称谱均乘以相应的因子才可得到相应变换后的结果；而基于Fourier变换的性质不需要求取关于零频的对称谱.但对于褶积定理而言，当实函数f₁(x, y)和f₂(x, y)中有一个函数为偶函数时，那么这两个函数的Hartley变换褶积频谱与Fourier变换褶积频谱形式相同(Saatcilar et al., 1990)，即HT[f₁(x, y)×f₂(x, y)]=H₁(u, v)·H₂(u, v).因此，通过Fourier变换和Hartley变换性质对比可以发现，Fourier变换的大部分性质比Hartley变换的性质更加简洁.

(3) Hartley变换和Fourier变换快速算法的计算量

设M_e表示实乘运算次数，A_e表示实加运算次数，M=2^s和N=2^t为点数.则一维快速Fourier变换的计算量分别为(蒋增荣等，1993)

二维快速Fourier变换的计算量分别为(蒋增荣等，1993)

一维、二维快速Hartley变换的计算量分别为(成礼智，1988)

通过快速Hartley变换计算量与快速Fourier变换计算量对比可得，实数域快速Fourier变换计算量约为复数域快速Fourier变换计算量的一半，近似与快速Hartley变换计算量相当.

这里的位场转换包括位场延拓、各阶方向导数转换、重磁异常分量转换和磁异常磁化方向转换.

设U_t1(x, y, z)，U_t2(x, y, z)分别为位函数U(x, y, z)在单位矢量t₁=(l₁, m₁, n₁)和t₂=(l₂, m₂, n₂)方向的位场分量.U_s1(x, y, z)，U_s2(x, y, z)分别表示同一分量在单位矢量s₁=(α₁, β₁, γ₁)和s₂=(α₂, β₂, γ₂)磁化方向的磁异常.U_{t, s}(x, y, z)表示地磁场方向为t =(l, m, n)和磁化方向为s =(α, β, γ)的磁异常分量.则有关Fourier变换和Hartley变换在位场数据处理和转换中的频率响应对比见表 3和表 4(据侯重初和刘秀芳，1978；魏雅利和骆遥，2010；骆遥，2013)所示.

表 3 (Table 3)

表 3 一维Hartley变换和Fourier变换在位场数据处理和转换中的频率响应对比表 Table 3 The contrast table of frequency response based on 1D Hartley transform and Fourier transform in potential field data processing and transformation

表 3 一维Hartley变换和Fourier变换在位场数据处理和转换中的频率响应对比表 Table 3 The contrast table of frequency response based on 1D Hartley transform and Fourier transform in potential field data processing and transformation

表 4 (Table 4)

表 4 二维Hartley变换和Fourier变换在位场处理和转换中的频率响应对比表 Table 4 The contrast table of frequency response based on 2D Hartley transform and Fourier transform in potential field data processing and transformation

表 4 二维Hartley变换和Fourier变换在位场处理和转换中的频率响应对比表 Table 4 The contrast table of frequency response based on 2D Hartley transform and Fourier transform in potential field data processing and transformation

综上，Hartley变换和Fourier变换的位场延拓表达式相同，Hartley变换的导数因子和方向转换因子要比Fourier变换的复杂.总的来说，在位场数据处理和转换中，Fourier变换的频谱响应表达式要比Hartley变换的频谱响应表达式更加简洁.

为了验证Hartley变换在位场数据处理和转换中的计算精度和计算时间，我们设计了单一直立六面体模型(图 1)进行测试，模型长、宽范围分别为-400~400 m(图中白色边框为模型的具体位置)，顶深50 m，厚750 m，模型体密度为0.7×10³ kg/m³，磁化强度为1.0 A/m，地磁场偏角为5°，倾角为50°.计算面为z=0平面规则网，计算点的变化范围在x方向上从-1000~1000 m，点距为10 m；y方向从-1000~1000 m，点距为10 m.图 1给出了在z=0 m处的理论重力异常平面图，图 2给出了在z=0 m处的理论磁力异常平面图，所有模型测试的扩边方法均用最小曲率扩边方法(王万银等，2009).

图 1

Figure 1

图 1 理论模型重力异常平面图 Figure 1 The plane map of the theoretical model's gravity anomaly

图 2

Figure 2

图 2 理论模型ΔT磁异常平面图 Figure 2 The plane map of the theoretical model's ΔT magnetic anomaly

(1) 位场解析延拓测试

为了验证Fourier变换和Hartley变换在位场解析延拓中的计算精度和计算时间.我们选取了观测面为z=0 m处的理论重力异常，分别向上延拓20 m和向下延拓20 m的重力异常，计算结果如图 3所示.

图 3

Figure 3

图 3 Fourier变换和Hartley变换重力异常延拓平面图 (a)和(d)分别为高度-20 m和20 m的理论重力异常；(b)和(c)分别为基于Hartley变换和Fourier变换的重力异常向上延拓20 m的结果；(e)和(f)分别为基于Hartley变换和Fourier变换的重力异常向下延拓20 m的结果. Figure 3 The plane map of the continued gravity anomaly based on Fourier transform and Hartley transform (a) and (d)Are the theoretical gravity anomaly in -20 m and 20 m; (b) and (c) Are the continued gravity anomaly based on Fourier transform and Hartley transform in -20 m, respectively; (e) and (f) Are the continued gravity anomaly based on Fourier transform and Hartley transform in 20 m, respectively.

图 3a是向上延拓20 m的理论重力异常平面图，图 3b和图 3c分别是对图 1所示重力异常利用Fourier变换和Hartley变换向上延拓20 m的重力异常平面图.图 3d是向下延拓20 m的理论重力异常图，图 3e和图 3f分别是对图 1所示重力异常利用Fourier变换和Hartley变换向下延拓20 m的重力异常平面图.从图 3对比发现，Fourier变换与Hartley变换的延拓结果基本相同.通过计算Fourier变换和Hartley变换向上延拓20 m的重力异常与该高度的理论重力异常的均方误差分别0.066×10^-6 m/s²和0.065×10^-6 m/s²，Fourier变换与Hartley变换向下延拓20 m的重力异常与该高度的理论重力异常的均方误差分别为0.152×10^-6 m/s²和0.151×10^-6 m/s².从均方误差上可以看出其计算精度相近.就计算量而言，实数域Fourier变换与Hartley变换的计算时间大约分别为0.802 s和0.796 s(计算机配置：内存4 GB，处理器Intel(R)Core(TM)i5 CPU 760 @ 2.80 GHz 2.79 GHz，显卡ATI Radeon HD 5750，下同).因此，Hartley变换与实数域Fourier变换位场延拓具有相同的精度，且计算时间也基本相同.

(2) 位场导数转换测试

为了对比Fourier变换和Hartley变换位场导数计算精度和计算时间，我们选用图 1理论重力异常，分别利用Hartley变换和Fourier变换计算x、y和z方向一阶导数(图 4).从图 4可以看出，Hartley变换与Fourier变换位场导数计算结果基本相同.Fourier变换在计算x、y和z方向一阶导数结果与理论异常导数之间的均方误差分别为1.51×10^-4 s^-2，1.51×10^-4 s^-2，3.39×10^-3 s^-2，Hartley变换在计算x、y和z方向一阶导数结果与理论异常导数之间的均方误差分别为1.51×10^-4 s^-2，1.50×10^-4 s^-2，3.39×10^-3 s^-2.实数域Fourier变换和Hartley变换位场导数转换所用时间(x，y，z方向导数计算时间基本相同)分别为0.792 s和0.780 s.因此，Hartley变换与Fourier变换在位场导数计算方面具有相同的计算精度，且与实数域Fourier变换计算时间基本相同.

图 4

Figure 4

图 4 Fourier变换和Hartley变换重力异常导数转换结果平面图 (a)，(b)和(c)分别是水平x方向一阶导数理论值，Fourier变换和Hartley变换计算结果；(d)，(e)和(f)分别是水平y方向一阶导数理论值，Fourier变换和Hartley变换计算结果；(g)，(h)和(i)分别是z方向一阶导数理论值，Fourier变换和Hartley变换计算结果. Figure 4 The plane map of the derivative transformation of gravity anomaly based on Fourier transform and Hartley transform (a), (b) and (c) Are theoretical, Fourier transform and Hartley transform of gravity anomalies about horizontal x direction first order derivative, respectively; (d), (e) and (f) Are theoretical, Fourier transform and Hartley transform of gravity anomalies about horizontal y direction first order derivative, respectively; (g), (h) and (i) Are theoretical, Fourier transform and Hartley transform of gravity anomalies about vertical z direction first order derivative, respectively.

(3) 磁异常分量转换和磁化方向转换测试

根据图 2所示磁力异常进行化极处理，理论化极磁异常如图 5a所示.图 5b和图 5c分别是基于Fourier变换和Hartley变换得到的化极磁异常，从图 5可以看出，Fourier变换和Hartley变换的化极磁异常基本相同.我们也计算了频率域化极磁异常与理论化极磁异常之间的均方差，实数域Fourier变换化极与理论化极磁异常之间的均方差为4.395 nT，Hartley变换化极与理论化极磁异常之间的均方差为4.394nT，从均方差可以看出，Hartley变换化极与Fourier变换化极具有相同的精度.在计算时间上，实数域Fourier变换用时0.802 s，Hartley变换化极用时0.794 s，可以看出Hartley变换与实数域Fourier变换计算时间基本相同.

图 5

Figure 5

图 5 Fourier变换和Hartley变换磁异常化极结果平面图 (a)，(b)和(c)分别是理论，Fourier变换和Hartley变换化极磁异常. Figure 5 The plane map of the reduction-to-the-pole magnetic anomaly based on Fourier transform and Hartley transform (a), (b) and (c) Are the Theoretical, Fourier transform and Hartley transform reduction-to-the-pole magnetic anomaly, respectively.

本次实际资料处理选取我国某地调平后的航磁资料进行处理，现给出实际资料的总场磁异常如图 6所示，然后利用Hartley变换和Fourier变换方法对该区的实际资料进行处理，从而比较Hartley变换和Fourier变换在位场处理和转换方面的效果，实际资料处理扩边方法均用最小曲率扩边方法.(王万银等，2009)

图 6

Figure 6

图 6 ΔT磁力异常平面图 Figure 6 The plane map of ΔT magnetic anomaly

(1) 位场解析延拓

利用Fourier变换和Hartley变换对图 6所示磁异常进行解析延拓.图 7a、图 7b是利用Fourier变换向上延拓300 m和向下延拓300 m的磁异常结果，图 7c和图 7d是利用Hartley变换向上延拓300 m和向下延拓300 m的磁异常结果.从图中可以发现，Hartley变换与Fourier变换位场延拓结果基本相同，在向上延拓中，Fourier变换与Hartley变换之间偏差的均方差为0.0956 nT，在向下延拓中，Fourier变换与Hartley变换之间偏差的均方差为0.158nT，这个偏差是由于数值计算误差造成，而不是方法本身造成.

图 7

Figure 7

图 7 Fourier变换和Hartley变换磁异常延拓平面图 (a)和(c)分别是基于Fourier变换和Hartley变换的磁异常向上延拓300 m结果；(b)和(d)分别是基于Fourier变换和Hartley变换的磁异常向下延拓300 m结果. Figure 7 The plane map of continuation magnetic anomaly based on Fourier transform and Hartley transform (a) and (c) Are continued magnetic anomaly based on Fourier transform and Hartley transform in -300 m, respectively; (b) and (d) Are continued magnetic anomaly based on Hartley transform and Hartley transform in 300 m, respectively.

在计算时间上，实数域Fourier变换所用时间为0.754 s，Hartley变换所用时间0.733 s.可以看出Hartley变换所用时间与实数域Fourier变换基本相同，其精度也相同.

(2) 位场导数计算

分别利用Fourier变换和Hartley变换计算图 6所示磁异常x、y、z方向一阶导数(图 8).从计算结果可以看出，Fourier变换和Hartley变换各方向导数计算结果基本相同.Hartley变换和Fourier变换位场导数在x、y和z三个方向一阶导数之间偏差的均方差分别为7.23×10^-6 nT/m、5.61×10^-6 nT/m和1.45×10^-4 nT/m.就计算量而言，Hartley变换和实数域Fourier变换的计算时间分别为0.7488 s和0.7686 s.因此，从计算结果可以看出，Hartley变换位场导数与Fourier变换位场导数处理结果基本相同，两者之间的偏差均方误差较小，表明这两种方法得到的计算结果基本相同；从计算时间来看，Hartley变换与实数域Fourier变换计算时间基本相同.

图 8

Figure 8

图 8 Fourier变换和Hartley变换磁异常导数转换结果平面图 (a)，(c)和(e)分别是基于Fourier变换的x、y和z方向的一阶导数计算结果；(b)，(d)和(f)分别是基于Hartley变换的x、y和z方向的一阶导数计算结果. Figure 8 The plane map of the derivative transformation of magnetic anomaly based on Fourier transform and Hartley transform (a), (c) and (e) Are the derivative of magnetic anomaly based on Fourier transform about x, y and z, respectively; (b), (d) and (f) Are the derivative of magnetic anomaly based on Hartley transform about x, y and z, respectively.

(3) 磁异常化极处理

对于磁异常分量转换和磁化方向转换，利用Fourier变换和Hartley变换对图 6所示磁力异常进行化极处理(图 9).通过对比发现，Fourier变换与Hartley变换化极结果基本相同，Hartley变换与Fourier变换之间的偏差均方差为0.0401 nT.在计算时间上，实数域Fourier变换用时0.8654 s，Hartley变换用时0.8634 s.因此，Fourier变换与Hartley变换所需计算时间基本相同，其精度也基本相同.

图 9

Figure 9

图 9 Fourier变换和Hartley变换化极磁异常平面图 (a)和(b)分别为Fourier变换和Hartley变换化极磁异常结果. Figure 9 The plane map of reduction-to-the-pole magnetic anomaly based on Fourier transform and Hartley transform (a) and (b) Are reduction-to-the-pole magnetic anomaly base on Fourier transform and Hartley transform respectively.

通过整理Hartley变换与Fourier的定义、性质、快速算法的计算量以及在位场数据处理和转换上的频率域响应，对Hartley变换和Fourier变换进行对比研究.通过对比认为，Hartley变换的性质、频率响应比Fourier变换复杂，Hartley变换的快速算法与实数域Fourier变换的快速算法计算量基本相同；通过理论模型测试得出其计算精度基本相同；通过理论模型测试和实际资料处理试验验证了其计算量基本相同.

因此，在位场数据处理和转换中，由于Hartley变换比Fourier变换的频率响应复杂，而且其快速算法的计算量和计算精度基本相当，从使用的角度来看，Fourier变换更具有优势.