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

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

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

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

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

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

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

0 引言

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

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

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

1.1 一維Hartley變換的定義

Hartley變換是酉變換的一種，變換前后的信號(hào)熵和能量不變[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ù)計(jì)算公式

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

對(duì)式(14)進(jìn)行反Hartley變換可得

(15)

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

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

(16)

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

(17)

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

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

2.1 二維Hartley變換的定義

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

(18)

(19)

對(duì)于二維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ù)計(jì)算公式

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

(25)

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

a.理論重力異常；b.理論一階垂直導(dǎo)數(shù)；c.利用Fourier變換計(jì)算得到的一階垂直導(dǎo)數(shù)；d.利用Hartley變換計(jì)算得到的一階垂直導(dǎo)數(shù)；e.Fourier變換計(jì)算導(dǎo)數(shù)與理論導(dǎo)數(shù)差；f.Hartley變換計(jì)算導(dǎo)數(shù)與理論導(dǎo)數(shù)差；g.Fourier變換計(jì)算得到的含噪異常的一階垂直導(dǎo)數(shù)；h.Hartley變換計(jì)算得到的含噪異常計(jì)算得到的一階垂直導(dǎo)數(shù)。圖2 長方體產(chǎn)生的重力異常及其導(dǎo)數(shù)計(jì)算結(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變換。

對(duì)式(30)進(jìn)行反Hartley變換可得

(31)

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

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

(32)

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

(33)

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

(34)

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

3 理論模型試驗(yàn)

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

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

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

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

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

ETM=

(35)

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

5 結(jié)論

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

[1] 馬國慶，杜曉娟，李麗麗. 利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進(jìn)行位場(chǎng)數(shù)據(jù)的邊界識(shí)別[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(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] 馬國慶，杜曉娟，李麗麗．解釋位場(chǎng)全張量數(shù)據(jù)的張量局部波數(shù)法及其與常規(guī)局部波數(shù)法的比較[J]．地球物理學(xué)報(bào)，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].國際地震動(dòng)態(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ù)字信號(hào)處理中的快速算法研究進(jìn)展[J]．高校應(yīng)用教學(xué)學(xué)報(bào)：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ì)地震波場(chǎng)[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變換的剖面位場(chǎng)轉(zhuǎn)換[J]. 地球物理學(xué)進(jìn)展, 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é)報(bào)，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] 穆石敏，申寧華，孫運(yùn)生.區(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] 姚長利，李宏偉，鄭元滿,等. 重磁位場(chǎng)轉(zhuǎn)換計(jì)算中迭代法的綜合分析與研究[J]. 地球物理學(xué)報(bào)，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] 馬國慶，黃大年，于平，等．改進(jìn)的均衡濾波器在位場(chǎng)數(shù)據(jù)邊界識(shí)別中的應(yīng)用[J]．地球物理學(xué)報(bào)，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

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

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

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

10.13278/j.cnki.jjuese.201401301

P631.1

A

馬國慶，黃大年，杜曉娟，等.Hartley變換在位場(chǎng)(重、磁)異常導(dǎo)數(shù)計(jì)算中的應(yīng)用.吉林大學(xué)學(xué)報(bào):地球科學(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.