移樣離散傅里葉變換在重磁勘探中的應(yīng)用

柴玉璞 萬(wàn)海珍

(東方地球物理公司綜合物化探處，河北涿州 072751)

0 引言

重磁資料反演是解釋人員從重磁資料中獲取地質(zhì)信息的重要手段，而絕大多數(shù)反演方法都是以正演方法為基礎(chǔ)構(gòu)建的。因此重磁數(shù)據(jù)的正演一直是地球物理學(xué)家研究的重點(diǎn)和熱點(diǎn)之一。

重磁正演問(wèn)題分兩大類(lèi)： 空間域正演和波數(shù)域正演。空間域正演的優(yōu)點(diǎn)是直觀、精度高；缺點(diǎn)是公式推導(dǎo)較難，導(dǎo)出結(jié)果繁瑣，對(duì)于較復(fù)雜的形體，則難以導(dǎo)出正演解析表達(dá)式。空間域正演的另一缺點(diǎn)是計(jì)算速度慢。波數(shù)域正演的優(yōu)點(diǎn)是公式推導(dǎo)簡(jiǎn)單，導(dǎo)出結(jié)果簡(jiǎn)潔。對(duì)于很多空間域無(wú)法導(dǎo)出的正演公式，在波數(shù)域中很容易導(dǎo)出。波數(shù)域正演的另一優(yōu)點(diǎn)是計(jì)算速度快。波數(shù)域正演的主要缺點(diǎn)是精度低。為了提高波數(shù)域正演計(jì)算精度，1988年，Chai等[1]提出了波數(shù)域最佳偏移抽樣法(偏移量為0.22倍波數(shù)域采樣間隔)，將重磁波數(shù)域正演精度提高了幾十倍。2014年，Wu等[2]將柴玉璞[3-4]提出的移樣離散傅里葉變換(SFT)與高斯節(jié)點(diǎn)積分相結(jié)合，提出了一種高精度快速位場(chǎng)波數(shù)域正演方法，對(duì)重磁勘探技術(shù)的發(fā)展做出了很大貢獻(xiàn)，但尚有以下不足。

(1)作者對(duì)所提方法的論證是從位函數(shù)譜的性質(zhì)出發(fā)，未能充分展現(xiàn)高斯FFT算法廣泛的應(yīng)用領(lǐng)域。

(2)文中缺少將高斯節(jié)點(diǎn)積分引入傅里葉變換數(shù)值計(jì)算的數(shù)學(xué)演繹過(guò)程，不免讓讀者認(rèn)為文中公式(19)僅源于觀察、公式(20)僅源于猜想。而此種感覺(jué)乃數(shù)學(xué)論證欠嚴(yán)謹(jǐn)所致。

(3)文中未給出偏移量與高斯節(jié)點(diǎn)的換算關(guān)系、權(quán)系數(shù)與高斯求積系數(shù)之間的換算關(guān)系，也未直接給出區(qū)間[0，1]內(nèi)大于2的高斯抽樣點(diǎn)(即偏移量)和權(quán)系數(shù)，無(wú)法實(shí)現(xiàn)編程。

本文從移樣離散傅里葉變換理論出發(fā)，對(duì)高斯FFT算法展開(kāi)論證，將其應(yīng)用范圍擴(kuò)展到任意有界函數(shù)的正、反傅里葉變換。文中給出了將高斯節(jié)點(diǎn)積分引入傅里葉變換數(shù)值計(jì)算的嚴(yán)謹(jǐn)數(shù)學(xué)演繹過(guò)程，并在此過(guò)程中導(dǎo)出了偏移量與高斯節(jié)點(diǎn)的換算關(guān)系以及權(quán)系數(shù)與高斯求積系數(shù)之間的換算關(guān)系。這樣借助數(shù)學(xué)手冊(cè)[5]中的高斯型求積節(jié)點(diǎn)和求積系數(shù)數(shù)據(jù)表，即可實(shí)現(xiàn)編程。

1 DFT變換的推廣

目前采用的離散傅里葉變換(DFT)表達(dá)式為

(1)

(2)

式(2)稱(chēng)為DFTξη變換，其中0≤ξ,η≤1。這一推廣的正確性顯而易見(jiàn)，把式(2)中的一個(gè)等式代入另一個(gè)等式，即可證明兩式的互逆性。利用新確立的DFTξ η變換可以表述一個(gè)重要定理——傅里葉變換離散化定理。

2 移樣離散傅里葉變換理論和算法

柴玉璞在其發(fā)表在《中國(guó)科學(xué)》的論文[3]中證明了一個(gè)重要定理——傅里葉變換離散化定理。該定理揭示了函數(shù)離散值與其譜離散值之間的關(guān)系，并據(jù)此給出了嚴(yán)謹(jǐn)?shù)摹⑼黄菩缘母道锶~變換數(shù)值計(jì)算理論——移樣離散傅里葉變換理論。

2.1 傅里葉變換離散化定理

設(shè)f(x)與F(u)為一傅里葉變換對(duì)，且均無(wú)無(wú)窮型間斷點(diǎn)。依據(jù)下式

(3)

將f(x)和F(u)離散化，并加權(quán)折疊成兩個(gè)長(zhǎng)度為N的復(fù)序列(dΔ=1/N)，則兩序列恰好構(gòu)成一DFTξη變換對(duì)

(4a)

或

(4b)

2.2 移樣離散傅里葉變換理論和算法

將式(4a) 中k=0和l=0的兩項(xiàng)置于方程左邊，其余項(xiàng)放在右邊，可得

DFTξ η[f(n+ξ)Δ]-F[(m+η)d]

(5a)

上式即為DFTξ η算法的誤差方程。同理可由式(4b) 導(dǎo)出IDFTξ η算法的誤差方程

IDFTξ η{F[(m+η)d]}-f[(n+ξ)Δ]

(5b)

式(5a)、式(5b)中，左端第一項(xiàng)分別為DFTξ η和IDFTξ η算法表達(dá)式，第二項(xiàng)分別為離散真譜和函數(shù)真值表達(dá)式，右端為誤差表達(dá)式。

實(shí)際上，式(5)以數(shù)學(xué)語(yǔ)言表述了一個(gè)嚴(yán)謹(jǐn)且具突破性的傅里葉變換數(shù)值計(jì)算理論——移樣離散傅里葉變換。其嚴(yán)謹(jǐn)性在于它不僅給出了算法表達(dá)式，而且給出了誤差表達(dá)式；其突破性在于它給出了更廣義、更普遍的新算法—— SFT算法(式(2))。利用這一新算法，既可由函數(shù)的任意一組等間隔抽樣值計(jì)算任意一組等間隔頻點(diǎn)上的譜值，也可以由任意一組等間隔頻點(diǎn)上的譜值計(jì)算任意一組函數(shù)等間隔抽樣值。這就為將高斯節(jié)點(diǎn)積分引入傅里葉變換數(shù)值計(jì)算創(chuàng)造了條件。因此，將高斯節(jié)點(diǎn)積分引入傅里葉變換數(shù)值計(jì)算，是移樣離散傅里葉變換理論的一個(gè)重要應(yīng)用。

3 將高斯節(jié)點(diǎn)積分引入傅里葉變換數(shù)值計(jì)算

高斯節(jié)點(diǎn)積分問(wèn)題有嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)證明[6]，比較麻煩，也很深?yuàn)W，涉及拉格朗日插值、多項(xiàng)式正交性、勒讓德函數(shù)等問(wèn)題。作為應(yīng)用者，僅將其結(jié)論深入淺出地表述如下： 對(duì)于區(qū)間[-1，1]上有界函數(shù)f(x)的積分，可用其M個(gè)點(diǎn)的函數(shù)值加權(quán)求和高精度逼近，其條件是這M個(gè)點(diǎn)需是M階勒讓德多項(xiàng)式pM(x)的M個(gè)根。高斯節(jié)點(diǎn)積分的數(shù)學(xué)表達(dá)為

(6)

(7)

高斯節(jié)點(diǎn)積分是針對(duì)實(shí)函數(shù)的。由于復(fù)數(shù)譜的實(shí)部和虛部都是實(shí)函數(shù)，所以正、反傅里葉變換都可以引入高斯節(jié)點(diǎn)積分。所謂將高斯節(jié)點(diǎn)積分引入傅里葉變換數(shù)值計(jì)算，就是將傅里葉積分的N個(gè)等間隔子區(qū)間[nΔ,(n+1)Δ]上的積分，轉(zhuǎn)換為區(qū)間[-1，1]上的高斯積分。眾所周知，對(duì)于變量t在區(qū)間[a，b]上的積分可通過(guò)變量替換

(8)

轉(zhuǎn)換為區(qū)間[-1，1]上的積分，即

(9)

利用式(9)將子區(qū)間[nΔ,(n+1)Δ]上的積分轉(zhuǎn)換為區(qū)間[-1，1]上的高斯節(jié)點(diǎn)積分

(10)

進(jìn)而可得

(11a)

(11b)

由此可見(jiàn)，傅里葉積分

(12a)

可用下式高精度逼近

(12b)

傅里葉積分

(13a)

可用下式高精度逼近

(13b)

由勒讓德多項(xiàng)式的性質(zhì)可知，高斯節(jié)點(diǎn)具有關(guān)于區(qū)間中心的對(duì)稱(chēng)性。式(4b)中的f為實(shí)函數(shù)，顯然，利用該式(其左端只保留l=0項(xiàng)，即忽略截?cái)嗾`差)可證明，移樣離散傅里葉反變換中，對(duì)稱(chēng)高斯節(jié)點(diǎn)序列的譜具有共軛性。因此，式(13b)可簡(jiǎn)化為

(13c)

式(12b)和式(13c)即是一維傅里葉正、反變換數(shù)值計(jì)算的實(shí)用公式。

柴玉璞在其專(zhuān)著[4]中將一維DFTξ η和IDFTξ η算法誤差方程推廣到二維，其形式與一維完全一致(式(5))。對(duì)于二維IDFTξ η算法誤差方程，略去有限效應(yīng)項(xiàng)，令ξ=0，并將左端k1=k2=0項(xiàng)移到右端，得到

f[(n1+k1N)Δ,(n2+k2N)Δ]

(14)

二維傅里葉變換相當(dāng)于進(jìn)行兩次一維傅里葉變換，因此，若將高斯節(jié)點(diǎn)積分引入二維傅里葉變換(方法同一維，只是分別對(duì)兩重積分實(shí)施)

(15a)

則其積分可用下式高精度逼近

(15b)

式中ξ1,i、ξ2,i分別表示第i個(gè)移樣離散傅里葉變換中x軸和y軸上的偏移量。同理，將高斯節(jié)點(diǎn)積分引入二維傅里葉反變換

(16a)

則上式可高精度逼近為

(16b)

式中η1,i、η2,i分別表示第i個(gè)移樣離散傅里葉變換中u軸和v軸上的偏移量。由式(14)可判斷，二維傅里葉反變換的第一重變換中，其對(duì)稱(chēng)高斯節(jié)點(diǎn)序列的譜無(wú)共軛性，但第二重變換其對(duì)稱(chēng)高斯節(jié)點(diǎn)序列的譜具有共軛性。因此式(16b)可簡(jiǎn)化為

(16c)

式(15b)和式(16c)即是二維傅里葉正、反變換數(shù)值計(jì)算的實(shí)用公式。式中兩重一維移樣離散傅里葉變換的高斯節(jié)點(diǎn)數(shù)理論上可以不同，但一般情況下無(wú)此必要，故本文未作區(qū)別。

柴玉璞在其專(zhuān)著[4]中詳細(xì)論述了SFT與DFT的關(guān)系，并指出SFT仍可借用FFT實(shí)現(xiàn)快速計(jì)算。其實(shí)，只需將DFTξη的定義式(式(2))展開(kāi)，上述關(guān)系和結(jié)論便一目了然。

4 位場(chǎng)波數(shù)域正演中的應(yīng)用

位場(chǎng)是有界函數(shù)，而且其譜呈負(fù)指數(shù)規(guī)律快速衰減，所以高斯FFT算法可用于位場(chǎng)反傅里葉變換，且效果特別好。位場(chǎng)反傅里葉變換，即位場(chǎng)波數(shù)域正演，也是重磁勘探的重要課題之一，也特別需要借助這一數(shù)學(xué)方法提高計(jì)算精度。所以本文選擇位場(chǎng)波數(shù)域正演作為這一數(shù)學(xué)方法的應(yīng)用領(lǐng)域，并沿用Wu等[2]對(duì)這一快速、高精度位場(chǎng)波數(shù)域正演方法的命名——高斯傅里葉正演。

設(shè)計(jì)一個(gè)由5個(gè)同尺寸的直立方柱體組成的模型。5個(gè)直立方柱體的高度均為2km，頂面埋深為1km，頂面面積為10km×10km，剩余密度均取2.0g/cm3。這5個(gè)方柱體的平面位置：其中一個(gè)的中心與圖1的中心點(diǎn)重合，其余4個(gè)的中心分別與圖1的四條邊的中點(diǎn)重合。正演計(jì)算面積為64km×64km，網(wǎng)格距為0.5km，網(wǎng)格數(shù)據(jù)點(diǎn)為128×128。圖1是采用空間域方法計(jì)算的模型重力異常。圖2為高斯抽樣點(diǎn)數(shù)分別為2、4、6時(shí)的高斯傅里葉正演結(jié)果及其對(duì)應(yīng)的正演誤差。高斯抽樣點(diǎn)數(shù)為2、4、6時(shí)的均方根誤差分別為13.300、0.060、0.001mGal，說(shuō)明隨著高斯抽樣點(diǎn)數(shù)的增加，誤差急劇減小，當(dāng)高斯抽樣點(diǎn)增加到4時(shí)，精度已可滿(mǎn)足常規(guī)解釋要求。

圖1 方柱體組合模型重力異常圖(空間域正演)

圖2 方柱體組合模型高斯抽樣點(diǎn)數(shù)為2(a)、4(b)、6(c)時(shí)的剩余重力異常圖(左)及誤差分布(右)

高斯傅里葉正演的一般規(guī)律是，場(chǎng)源離計(jì)算窗中心越遠(yuǎn)，其譜的波動(dòng)頻率越高，則需要更多的高斯抽樣點(diǎn)才能達(dá)到期望精度，特別是在源跨越計(jì)算窗邊、延伸到計(jì)算窗外的情況。當(dāng)高斯抽樣點(diǎn)為2時(shí)，中心方柱體的正演誤差已經(jīng)很小，以至于可以忽略，但周?chē)膫€(gè)方柱體的正演誤差仍相當(dāng)大(圖2a右)；當(dāng)高斯抽樣點(diǎn)為4時(shí)，周?chē)膫€(gè)方柱體的正演誤差已微乎其微(圖2b右)；當(dāng)高斯抽樣點(diǎn)為6時(shí)，整個(gè)計(jì)算區(qū)域內(nèi)的誤差已經(jīng)降至非常低(約10-3mGal)(圖2c右)。實(shí)際上此時(shí)的誤差是譜的截?cái)嗾`差，它不再隨高斯抽樣點(diǎn)的增加而減小。

理論和實(shí)踐都證明，無(wú)論場(chǎng)源分布如何，高斯傅里葉正演都能通過(guò)增加高斯抽樣點(diǎn)達(dá)到預(yù)期的精度，且仍能保持位場(chǎng)波數(shù)域正演計(jì)算速度快的優(yōu)點(diǎn)。

5 結(jié)論

(1)快速傅里葉變換(FFT)解決了傅里葉變換數(shù)值計(jì)算的效率問(wèn)題，但其抽樣點(diǎn)需是一個(gè)傅里葉積分被劃分成的數(shù)個(gè)子區(qū)間的端點(diǎn)；高斯節(jié)點(diǎn)積分具有很高精度，但其抽樣點(diǎn)須是子區(qū)間內(nèi)的某些點(diǎn)。移樣離散傅里葉變換(SFT)理論的確立，客觀上為解決這一矛盾創(chuàng)造了條件。因此可以說(shuō)，在地球物理學(xué)家將FFT 和高斯節(jié)點(diǎn)積分相結(jié)合、提出位場(chǎng)高斯傅里葉正演算法的創(chuàng)新中，SFT功不可沒(méi)。

(2)本文以嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)演繹，不僅導(dǎo)出了傅里葉積分可以數(shù)個(gè)移樣離散傅里葉變換(SFT)的加權(quán)求和高精度地逼近的結(jié)論，而且導(dǎo)出了偏移量與高斯節(jié)點(diǎn)的換算關(guān)系、權(quán)系數(shù)與高斯求積系數(shù)之間的換算關(guān)系。這為位場(chǎng)高斯傅里葉正演方法的推廣提供了嚴(yán)謹(jǐn)?shù)睦碚撝魏途幊虒?shí)現(xiàn)上的可行性和方便條件。

(3)高斯FFT算法已經(jīng)在位場(chǎng)波數(shù)域正演中得到成功應(yīng)用。因?yàn)楦咚构?jié)點(diǎn)積分理論和移樣離散傅里葉變換理論的充分條件，都是任意有界函數(shù)，因此有理由相信，該方法未來(lái)有可能在更廣闊的領(lǐng)域獲得更多應(yīng)用。

感謝中國(guó)石油集團(tuán)東方地球物理公司賈繼軍、鄧玉友和孫喜明三位高級(jí)工程師在計(jì)算機(jī)操作和繪圖方面提供的幫助。