Hartley變換在位場(重、磁)異常導(dǎo)數(shù)計算中的應(yīng)用

馬國慶，黃大年，杜曉娟，李麗麗

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130061

Hartley變換在位場(重、磁)異常導(dǎo)數(shù)計算中的應(yīng)用

馬國慶，黃大年，杜曉娟，李麗麗

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130061

導(dǎo)數(shù)計算是位場(重磁)數(shù)據(jù)處理中必不可少的技術(shù)手段，現(xiàn)今大多采用Fourier變換來進行。Hartley變換是在Fourier變換基礎(chǔ)上定義的一種實數(shù)域運算，比Fourier變換更加對稱，所需的運算量更少。筆者推導(dǎo)出基于Hartley變換的重磁異常導(dǎo)數(shù)計算公式，計算結(jié)果與理論值之間誤差小于5%，通過理論模型證明Hartley變換可替代Fourier變換進行位場異常的導(dǎo)數(shù)計算，且受噪音干擾較小。將Hartley變換用于位場邊界識別濾波器的計算，獲得了清晰的斷裂分布。

位場；導(dǎo)數(shù)； Fourier變換；Hartley變換

0 引言

導(dǎo)數(shù)計算被廣泛應(yīng)用于重磁數(shù)據(jù)的處理與解釋[1-3]，具有突出淺源異常、區(qū)分疊加異常、確定地質(zhì)體邊界以及削弱背景異常的作用?，F(xiàn)今大多采用Fourier變換[4]來完成重磁異常導(dǎo)數(shù)的計算。Fourier變換要求數(shù)據(jù)的尺寸為2的冪次，而Hartley變換無此要求。Hartley變換是在Fourier變換基礎(chǔ)上定義的一種實數(shù)域內(nèi)的運算[5]，至少比Fourier變換減少一半的空間和時間[6]，且更加對稱。Hartley變換已經(jīng)在圖像處理、模式識別、地震波場模擬等領(lǐng)域得到廣泛應(yīng)用[7-11]；在位場數(shù)據(jù)處理中最早被用來進行波譜分析[12-13]，其得到的異常波譜也能成功地區(qū)分出背景場和局部場。Sundarajan等[14]將Hartley變換與Hilbert變換相結(jié)合進行功率譜分析，并證明采用該方法所需的計算量相對Fourier變換要少。后來人們利用Hartley變換進行位場剖面數(shù)據(jù)的相關(guān)分析與分量之間的轉(zhuǎn)換[15-16]，均取得了較好的應(yīng)用效果。

筆者根據(jù)Hartley變換的基本性質(zhì)及位場基本原理推導(dǎo)出基于Hartley變換的重磁異常二維和三維導(dǎo)數(shù)計算公式。通過理論模型試驗證明Hartley變換計算得到的導(dǎo)數(shù)與Fourier變換計算結(jié)果相接近，且與理論值誤差較小，可作為進行位場數(shù)據(jù)處理與轉(zhuǎn)換的另一種方法。

1 二維重(磁)異常頻譜及導(dǎo)數(shù)計算公式推導(dǎo)

1.1 一維Hartley變換的定義

Hartley變換是酉變換的一種，變換前后的信號熵和能量不變[5]。Hartley變換的一維形式如下:

正變換為

(1)

逆變換為

(2)

由Hartley變換定義奇函數(shù)與偶函數(shù)如下:

奇函數(shù)為

(3)

偶函數(shù)為

(4)

1.2 一維Hartley變換的性質(zhì)

性質(zhì)1：Hartley變換與自身奇函數(shù)和偶函數(shù)的關(guān)系[5-6]為

(5)

(6)

性質(zhì)2：Hartley變換與Fourier變換的關(guān)系式為

(7)

(8)

(9)

式中：PH(w)、XH(w)分別為pxy(τ)、x(t)的Hartley變換；YHe(w)、YHo(w)分別為y(t)的偶函數(shù)變換和奇函數(shù)變換。

(10)

1.3 二維重磁異常Hartley譜及導(dǎo)數(shù)計算公式

重力場源在z≤0空間中的重力異常f(x,z)是調(diào)和函數(shù)。已知z=0(水平面上)的函數(shù)值f(ξ,0)，求z≤0空間中的重力函數(shù)值f(x,z)，這種數(shù)學(xué)問題為狄里希萊(Dirichlet)問題(第一邊值問題)。利用格林函數(shù)可求得上半空間解的積分表達式。對于二維問題，延拓積分表達式為[17]

(11)

已知Fourier變換關(guān)系式[18]：

(12)

(13)

由性質(zhì)3可得f(x,z)的Hartley變換式為

(14)

式中，F(xiàn)H(u,0)為f(x,0)的Hartley變換。

對式(14)進行反Hartley變換可得

(15)

應(yīng)用微分定理對式(15)做積分號下求微分的計算，即可求出重、磁異常導(dǎo)數(shù)的表達式。

水平導(dǎo)數(shù)：

(16)

垂直導(dǎo)數(shù)：

(17)

從式(16)和(17)中可以看出，Hartley在進行導(dǎo)數(shù)計算時，均是實數(shù)域內(nèi)的運算，不需要虛數(shù)參與計算,因此比Fourier變換減少一半的運算量。

2 三維重(磁)異常Hartley譜及導(dǎo)數(shù)計算公式推導(dǎo)

2.1 二維Hartley變換的定義

Hartley變換二維形式與二維Fourier變換不同。對于二維實函數(shù)的Hartley變換[9]，其積分核存在2種形式：cas(ux+vy)、cas(ux)cas(vy)。為了便于計算，筆者采取可分離的第二種形式，并有如下的正、逆變換表示式為

(18)

(19)

對于二維Hartley變換也可以構(gòu)造出奇、偶函數(shù)，其形式如下:

奇函數(shù)為

(20)

偶函數(shù)為

(21)

2.2 二維Hartley變換的基本性質(zhì)

性質(zhì)5：二維Hartley變換與Fourier變換的關(guān)系式為

(22)

性質(zhì)6：二維Hartley變換的褶積關(guān)系式為

(23)

(24)

2.3 三維重(磁)異常Hartley譜及導(dǎo)數(shù)計算公式

重力場源在z≤0空間中的重力異常f(x,y,z)是調(diào)和函數(shù)。對于三維問題，延拓積分表達式為

(25)

a.一階水平導(dǎo)數(shù)；b.未擴邊一階垂直導(dǎo)數(shù)；c.擴邊后一階垂直導(dǎo)數(shù)；d.二階垂直導(dǎo)數(shù)。圖1 不同方法計算的重力異常導(dǎo)數(shù)結(jié)果Fig.1 The derivative results of gravity anomaly computed by different methods

a.理論重力異常；b.理論一階垂直導(dǎo)數(shù)；c.利用Fourier變換計算得到的一階垂直導(dǎo)數(shù)；d.利用Hartley變換計算得到的一階垂直導(dǎo)數(shù)；e.Fourier變換計算導(dǎo)數(shù)與理論導(dǎo)數(shù)差；f.Hartley變換計算導(dǎo)數(shù)與理論導(dǎo)數(shù)差；g.Fourier變換計算得到的含噪異常的一階垂直導(dǎo)數(shù)；h.Hartley變換計算得到的含噪異常計算得到的一階垂直導(dǎo)數(shù)。圖2 長方體產(chǎn)生的重力異常及其導(dǎo)數(shù)計算結(jié)果Fig.2 Gravity anomaly and derivative data generated by rectangular prism

已知Fourier關(guān)系式[18]：

(26)

由性質(zhì)5可得

(27)

(28)

因此，φ(ε,η)的Hartley變換為

(29)

由性質(zhì)6可得，f(x,y,z)的Hartley變換式為

(30)

式中：FH(u,v,z)為f(x,y,z)的Hartley變換，F(xiàn)H(u,v,0)為f(x,y,0)的Hartley變換。

對式(30)進行反Hartley變換可得

(31)

應(yīng)用微分定理對式(31)做積分號下求微分運算，即可求出重、磁異常導(dǎo)數(shù)的表達式。

x方向?qū)?shù)：

(32)

y方向?qū)?shù)：

(33)

z方向?qū)?shù)：

(34)

從式(32)、(33)和(34)中可以看出，三維導(dǎo)數(shù)計算公式與二維導(dǎo)數(shù)計算公式的推導(dǎo)過程不同，但最終的導(dǎo)數(shù)形式是一致的。

3 理論模型試驗

以二維情況為例：水平無限延伸圓柱體半徑R=50 m，中心點埋深h=100 m，與圍巖密度差ρ=1 g/cm3。分別利用離散Fourier和離散Hartley變換計算z=0觀測面上重力異常的導(dǎo)數(shù)(圖1)，并統(tǒng)計2種方法計算得到的導(dǎo)數(shù)與理論導(dǎo)數(shù)之間的均方差(表1)。

從表1中可以看出，采用Hartley變換計算得到的垂直導(dǎo)數(shù)的精度要略高于Fourier變換計算結(jié)果的精度，可作為位場導(dǎo)數(shù)計算的另一種選擇。

表1 不同方法計算導(dǎo)數(shù)與理論導(dǎo)數(shù)的均方差

Table 1 The error between theoretical derivative and computed derivative by different methods

Fourier變換Hartley變換一階水平導(dǎo)數(shù)2.04×10-63.41×10-6一階垂直導(dǎo)數(shù)2.86×10-52.68×10-5二階垂直導(dǎo)數(shù)1.33×10-32.72×10-4

噪聲是實際數(shù)據(jù)處理中不可避免的影響因素，試驗Hartley變換和Fourier變換在存在噪聲情況下的導(dǎo)數(shù)計算效果。圖2為埋深分別為15 m和20 m的長方體產(chǎn)生的重力異常。從圖2中可以看出，Hartley變換計算結(jié)果的幅值與理論值更為接近。從圖2g、h可以看出，Hartley變換受噪音干擾較小，依舊能很好地完成導(dǎo)數(shù)的計算，而Fourier變換的導(dǎo)數(shù)結(jié)果受噪音影響相對較大。這主要是由于在頻率域中轉(zhuǎn)換因子具有明顯的噪聲放大作用[19]，而Hartley變換為實數(shù)域內(nèi)的運算，不會明顯地增大噪聲的干擾。

4 實際數(shù)據(jù)應(yīng)用

為了進一步試驗本文方法的計算效果，分別利用Hartley變換與Fourier變換完成邊界識別濾波器-增強型?圖(enhanced theta map,ETM)[20]的計算，其表達式為

ETM=

(35)

式中，mean代表均值。利用不同方法計算得到的增強型?圖對中國西北部朱日和地區(qū)磁異常進行邊界識別處理，計算結(jié)果(圖3)均可以清晰地反映出地層之間的界限。

5 結(jié)論

Hartley變換是在Fourier變換基礎(chǔ)上定義的一種實數(shù)域計算，筆者推導(dǎo)出了基于Hartley變換的位場(重、磁)異常導(dǎo)數(shù)計算公式。通過理論模型試驗證明，Hartley變換計算結(jié)果與Fourier變換計算結(jié)果接近，與理論值之間誤差較小。將Hartley變換和Fourier變換應(yīng)用于位場數(shù)據(jù)邊界識別濾波器的計算，均能清晰地獲得地層之間的界限特征，因此,Hartley變換可作為位場數(shù)據(jù)處理與轉(zhuǎn)換的另一種選擇。

[1] 馬國慶，杜曉娟，李麗麗. 利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進行位場數(shù)據(jù)的邊界識別[J]. 吉林大學(xué)學(xué)報：地球科學(xué)版， 2011, 41(增刊1): 345-348. Ma Guoqing, Du Xiaojuan, Li Lili. Edge Detection of Potential Field Data Using Correlation Coefficients of Horizontal and Vertical Derivatives[J]. Journal of Jilin University: Earth Science Edition, 2011, 41(Sup.1): 345-348.

[2] Ma G, Li L. Edge Detection in Potential Fields with the Normalized Total Horizontal Derivative[J]. Computers & Geosciences, 2012, 41: 83-87.

[3] 馬國慶，杜曉娟，李麗麗．解釋位場全張量數(shù)據(jù)的張量局部波數(shù)法及其與常規(guī)局部波數(shù)法的比較[J]．地球物理學(xué)報，2012，55(7)：2450-2461. Ma Guoqing, Du Xiaojuan, Li Lili. Comparison of the Tensor Local Wavenumber Method with the Conventional Local Wavenumber Method for Interpretation of Total Tensor Data of Potential Fields[J]. Chinese Journal of Geophysics, 2012，55(7)：2450-2461.

[4] Blakely R J. Potential Theory in Gravity and Magnetic Applications[M].Cambridge: Cambridge University Press, 1996: 256.

[5] Harteley R V L. A More Symmetrical Fourier Anal-ysis Applied to Transmission Problems[J]. Proc IRE, 1942, 30(2): 144-150.

[6] 余品能，路凌云.離散Hartley變換的一種快速遞歸算法[J].石油地球物理勘探，1998，33(5)：591-596. Yu Pinneng, Lu Lingyun. A Fast Recursive Algorithm of Discrete Hartley Transform[J]. Oil Geophysical Prospecting，1998，33(5)：591-596.

[7] 彭軍，羅奇峰. Hartley變換在地震波研究中的應(yīng)用[J].國際地震動態(tài)，2008(12)：6-8. Peng Jun, Luo Qifeng. The Application of Hartley Transform to the Study on Seismic Wave[J]. Recent Developments in World Seismolog, 2008(12)：6-8.

[8] 余品能，蔣增榮．圖像及數(shù)字信號處理中的快速算法研究進展[J]．高校應(yīng)用教學(xué)學(xué)報：A輯，1991， 6(2)：302-316. Yu Pinneng, Jiang Zengrong. Advances in the Study of Fast Algorithms in Image and Digital Signal Processing[J]. Applied Mathematics: A Journal of Chinese Universities, 1991，6(2)：302-316.

[9] 甘利燈．哈特萊變換及其性質(zhì)[J]．石油地球物理勘探，1992，27(5)：605-616. Gan Lideng. Hartley Transform and Its Prosperities[J]. Oil Geophysics Prospecting, 1992，27(5)：605-616.

[10] 周輝，何樵登．利用Hartley變換模擬各向異性介質(zhì)地震波場[J]．石油地球物理勘探，1995，30(5)：593-601. Zhou Hui, He Qiaodeng. Seismic Wave Modeling in Anisotropic Medium by Hartley Transform[J]. Oil Geophysics Prospecting, 1995，30(5)：593-601.

[11] Saatcilar R. The Use of Hartley Transform in Geo-physical Applications[J]. Geophysics,1990, 55(11):1488-1495.

[12] Bravo R, Escalona O, Mora F, et al. Vector Spectral Combination of Orthogonal Leads Using the Hartley Transform for Late Potential Analysis in Chagas Disease[J]. Computers in Cardiology, 1997, 10:415-417.

[13] Mansour A, Garni A, Narasimman S . Hartley Spectral Analysis of Self-Potential Anomalies Caused by a 2-D Horizontal Circular Cylinder[J]. Arabian Journal of Geosciences, 2012, 5(6): 1341-1346.

[14] Sundararajan N, Srinivas Y. Fourier-Hilbert Versus Hartley-Hilbert Transform with Some Geophysical Applications[J]. Journal of Applied Geophysics, 2010,71(4):157-161.

[15] 魏雅利, 駱遙.基于Hartle變換的剖面位場轉(zhuǎn)換[J]. 地球物理學(xué)進展, 2010,25(6): 2102-2108. Wei Yali, Luo Yao. 2D Potential Field Transformation Based on Hartley Transform[J]. Progress in Geophysics, 2010, 25(6): 2102-2108.

[16] 孫鶴泉,沈永明,王永學(xué),等. Hartley變換在互相關(guān)分析中應(yīng)用研究[J].大連理工大學(xué)學(xué)報，2004，44(2)：285-286. Sun Hequan, Shen Yongming, Wang Yongxue, et al. Application of Hartley Transform to Cross-Correlation Analysis[J]. Journal of Dalian University of Technology, 2004，44(2)：285-286.

[17] Bracewell R N．Discrete Hartley Transform[J]．Journal of the Optional Society of America，1983，73：1832-1835.

[18] 穆石敏，申寧華，孫運生.區(qū)域地球物理數(shù)據(jù)處理方法及其應(yīng)用[M].長春：吉林科學(xué)技術(shù)出版社，1990：7-10. Mu Shimin, Shen Ninghua, Sun Yunsheng. Regional Geophysical Data Processing Method and Its Application[M]. Changchun: Jilin Science & Technology Press, 1990: 7-10.

[19] 姚長利，李宏偉，鄭元滿,等. 重磁位場轉(zhuǎn)換計算中迭代法的綜合分析與研究[J]. 地球物理學(xué)報，2012,55(6)：2062-2078. Yao Changli, Li Hongwei, Zheng Yuanman, et al. Research on Iteration Method Using in Potential Field Transformation[J]. Chinese Journal of Geophysics, 2012,55(6)：2062-2078.

[20] 馬國慶，黃大年，于平，等．改進的均衡濾波器在位場數(shù)據(jù)邊界識別中的應(yīng)用[J]．地球物理學(xué)報，2012,55(12)：4288-4295. Ma Guoqing, Huang Danian, Yu Ping, et al. Application of Improved Balancing Filters to Edge Identification of Potential Field Data[J]. Chinese Journal of Geophysics, 2012, 55(12): 4288-4295.

Hartley Transform in the Application of the Derivatives of Potential Field (Gravity and Magnetic) Data

Ma Guoqing, Huang Danian, Du Xiaojuan, Li Lili

CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130021,China

Derivative is an indispensable tool in the processing of potential field data. Now we usually use the Fourier transform to complete the computation, but this method is sensitive to noise, so cannot compute higher-order derivatives. Hartley transform (HT) is a real-valued function that is defined on the basis of Fourier transform (FT), however, symmetry character of HT is better and computational complexity of HT is smaller compared with these properties of FT. We derive the derivative computation equations of gravity and magnetic anomaly based on the basic property of HT. We demonstrate the HT on theoretical anomalies, and errors between the results computed by HT, which is insensitive to noise, and the theoretical values are less than 5%, so the HT can substitute the FT to compute the derivatives of potential field data. We also apply the HT to finish the computation of edge detection filters and obtain clearer distribution of faults.

potential field; derivative; Fourier transform; Hartley transform

10.13278/j.cnki.jjuese.201401301.

2013-06-29

國家科技專項項目(SinoProbe-09-01 201011078)；中國地質(zhì)調(diào)查局地質(zhì)礦產(chǎn)調(diào)查評價專項項目(GZH003-07-03)

馬國慶(1984-)，男，講師，博士，主要從事位場數(shù)據(jù)處理與解釋方面的研究，E-mail:magq08@mails.jlu.edu.cn

杜曉娟(1957-)，女，教授，主要從事位場數(shù)據(jù)解釋方面的研究，E-mail:dtdxj@jlu.edu.cn。

10.13278/j.cnki.jjuese.201401301

P631.1

A

馬國慶，黃大年，杜曉娟，等.Hartley變換在位場(重、磁)異常導(dǎo)數(shù)計算中的應(yīng)用.吉林大學(xué)學(xué)報:地球科學(xué)版,2014,44(1):328-335.

Ma Guoqing, Huang Danian, Du Xiaojuan,et al.Hartley Transform in the Application of the Derivatives of Potential Field (Gravity and Magnetic) Data.Journal of Jilin University:Earth Science Edition,2014,44(1):328-335.doi:10.13278/j.cnki.jjuese.201401301.