曾 明， 柳建新， 陳 波, 趙廣東, 陳龍偉

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長沙 410083；3 桂林理工大學(xué) 地球科學(xué)學(xué)院，桂林 541006)

0 引言

磁異常正演模擬對于磁異常解釋及反演相當(dāng)重要，為了定性地解釋觀測磁異常，通常根據(jù)已知先驗(yàn)信息來假定物理參數(shù)來模擬地下地質(zhì)體產(chǎn)生的磁異常，再與觀測磁異常進(jìn)行比較，從而推算地下復(fù)雜地質(zhì)體的地球物理特征。磁異常正演的效率和精度對磁異常反演結(jié)果準(zhǔn)確性以及反演效率有直接的影響。磁異常ΔT正演方法可分為空間域和頻率域兩大類。

為了計算地下復(fù)雜地質(zhì)體磁異常，空間域方法一般是將地下地質(zhì)體剖分成多個便于計算的規(guī)則單元體，分別計算每個單元體在觀測點(diǎn)處產(chǎn)生的磁異常，最后累加求和得到整個地質(zhì)體產(chǎn)生的總磁異常，如球體(磁偶極子)，圓柱體，多邊形，長方體，多面體，薄板等簡單幾何體，以及任意二度體，二度半體以及三度體等磁異常表達(dá)式逐漸被推導(dǎo)出來[1-8]。由于長方體能簡單有效地近似模擬復(fù)雜異常體，被廣泛應(yīng)用于正演模擬及位場反演領(lǐng)域。郭志宏等[9]推導(dǎo)了長方體無解析奇點(diǎn)磁異常表達(dá)式。由于空間域方法的解析式較為復(fù)雜，當(dāng)研究區(qū)域較大，或磁化率分布不均勻時，需要對場源區(qū)域進(jìn)行更多、更細(xì)的剖分，空間域算法的計算耗時將會大大提高，限制了其在正演模擬及反演方面的應(yīng)用。

頻率域方法因其計算速度快，頻譜表達(dá)式簡單、緊湊也迅速發(fā)展起來。Bhattacharyya[10]基于傅里葉變換推導(dǎo)了任意磁化方向長方體二維磁異常連續(xù)頻譜表達(dá)式；Schouten[11-12]基于FFT的方法模擬水平層狀源的磁異常；Bhattacharyya[13]推導(dǎo)了任意二度體磁位頻譜表達(dá)式，并用于頻率域反演研究；Pedersen[14]推導(dǎo)了任意二度半體磁異常頻譜表達(dá)式，并改進(jìn)了多面體頻率域表達(dá)式；吳宣志[15]給出任意指向的均質(zhì)直線段、多邊形面和多面體的磁異常譜表達(dá)式，并利用它們進(jìn)行不規(guī)則三度體磁異常的正演計算，且導(dǎo)出了任意磁化方向的斜平行六面體磁異常頻譜表達(dá)式；Tontini[16]推導(dǎo)了物性呈三維高斯分布的場源的磁場頻率域表達(dá)式；Tontini[17]對磁位三維全空間頻率域正演進(jìn)行了研究，給進(jìn)一步推導(dǎo)了直立長方體、球體的磁異常頻譜表達(dá)式。

利用快速傅里葉變換(FFT)對磁異常進(jìn)行正演，相比于空間域算法計算效率有很大提高，但由于離散傅里葉變換一些固有的缺陷使其正演精度較低(如混疊失真、頻譜泄露以及邊緣效應(yīng)等的影響)。為了提高計算精度，Tontini[17]通過擴(kuò)邊的方法提高計算精度；Chai[18]通過偏移采樣的方法減少了傅里葉正演模型的計算誤差，Wu等[19]指出快速傅里葉變換(FFT)采用的“矩形積分”不能很好地逼近傅里葉震蕩積分，在偏移采樣技術(shù)的基礎(chǔ)上引入了Gauss-FFT技術(shù)提高了二維傅里葉磁異常正演的精度，并應(yīng)用到二維位場正演計算。

由于三維傅里葉正演的精度高以及在磁異常反演領(lǐng)域的潛在應(yīng)用，筆者重新推導(dǎo)了長方體三維磁異常頻譜，引入Gauss-FFT技術(shù)用于三維磁異常頻率域正演，并通過數(shù)值模型驗(yàn)證其相對于空間域在計算用時上的優(yōu)越性，以及相對于快速三維傅里葉正演在計算精度上的優(yōu)勢。

1 方法原理

1.1 3D傅里葉正演理論

圖1 笛卡爾坐標(biāo)系下磁場源離散化示意圖Fig.1 Schematic diagram of the discretization of magnetic field sources

任何磁性分布地質(zhì)體都可以用許多磁化率均勻分布的長方體模型單元來模擬，通過計算各模型單元體在觀測點(diǎn)產(chǎn)生的磁異常，再對所有模型單元磁異常進(jìn)行簡單的疊加求和，從而得到整個地質(zhì)體產(chǎn)生的磁異常。如圖1所示為笛卡爾坐標(biāo)系下，給定磁化率大小及方向的場源離散化示意圖，z軸向下為正，整個場源區(qū)域按照x、y、z方向排列順序被均分為M×N×L棱柱體。各單元棱柱體的中心位置坐標(biāo)為(ξa,ηb,ζc)，x、y、z方向各單元體的間距分別為Δx、Δy、Δz；a、b、c分別對應(yīng)x、y、z方向上的模型單元體的位置。假定各個直角棱柱體的磁化率大小為j(ξa,ηb,ζc)，磁化強(qiáng)度方向?yàn)?l,m,n)，正常場方向?yàn)?l0,m0,n0)，大小為T0。正演觀測點(diǎn)P(x,y,z)為正演區(qū)域內(nèi)任意單元體中心點(diǎn)。

在空間域中，由任一單元體在觀測點(diǎn)P(x,y,z)產(chǎn)生的磁位U(x,y,z)可以表示為[20]：

(1)

其中：u0為真空磁導(dǎo)率，

(2)

則磁異常在空間域的積分表達(dá)式可以表示為：

(3)

其中:?/?/t、?//?/s分別表示正常磁場和磁化的方向余弦:

(4)

將式(4)代入式(3)并進(jìn)行三維傅里葉變換：

(5)

其中kx、ky、kz分別為x、y、z方向的波數(shù)。

(6)

其中k為波數(shù)，

根據(jù)傅里葉變換微分定理：

(7)

將式(6)、式(7)代入式(5)得到：

(lkx+mky+nkz)·

e-i(kxξ+kyη+kzζ)dξdηdζ

(8)

對整個場源區(qū)域模型單元體磁異常頻譜進(jìn)行累加求和到：

n0kz)(lkx+mky+nkz)

F[j(ξa,ηb,ζc)]

(9)

e-i(kxξ+kyη+kzζ)ΔxΔyΔz

(10)

再對式(9)中F[ΔT(kx,ky,kz)]進(jìn)行三維傅里葉逆變換從而得到整個場源區(qū)域磁異常。

1.2 3D Gauss-FFT技術(shù)

在實(shí)際應(yīng)用中，通常是采用離散數(shù)據(jù)計算磁異常，這就意味著式(10)應(yīng)該用相應(yīng)的離散傅里葉變換，但離散傅里葉變換因?yàn)榻財嗪碗x散化會不可避免地產(chǎn)生混疊失真、頻率泄漏、邊緣效應(yīng)等一些誤差[20]。Wu等[19]研究了二維Gauss-FFT方法來提高二維離散傅里葉逆變換的精度。

為了便于計算，我們用g(xa,yb,zc)表示各單元體中心的坐標(biāo)位置，x、y、z方向的觀測間隔為Δx、Δy、Δz，點(diǎn)S為第一個單元體的中心坐標(biāo)為(x0,y0,z0)，則觀測數(shù)據(jù)的離散坐標(biāo)g(xa,yb,zc)為：

xa=x0+aΔx;a=1,2,…,M-1

yb=y0+bΔy;b=1,2,…,N-1

zc=z0+cΔz;c=1,2,…,L-1

(11)

對于場源區(qū)域被均分M×N×L棱柱體，在尼奎斯特頻率約束條件下的離散波數(shù)可以表示為[22]:

(12)

則其三維離散傅里葉正變換(DFT)和三維離散傅里葉逆變換(IDFT)分別表示為：

e-i(kxpxa+kyqyb+kzτzc)ΔxΔyΔz

(13)

ei(kxpxa+kyqyb+kzτzc)ΔkxΔkyΔkz

(14)

式(13)、式(14)是采用矩形積分來近似替代連續(xù)傅里葉變換，式(13)具有與式(10)相同的形式，為了提高三維傅里葉變換的精度，采用Gauss積分的方法代替矩形積分以提高三維離散傅里葉變換的精度。一維的Gauss積分公式如下：

(15)

其中:h表示Gauss積分節(jié)點(diǎn)數(shù)，F(xiàn)h和th表示在[-1,1]區(qū)間上的高斯系數(shù)和高斯節(jié)點(diǎn)值。

假設(shè)三維Gauss積分在x、y、z方向的高斯節(jié)點(diǎn)數(shù)為Ix、Iy、Iz，則高斯系數(shù)及節(jié)點(diǎn)值在[0,1]區(qū)間的值分別為(λix, αix)、(λiy, αiy)、(λiz, αiz)，其中ix=1、2…、Ix，iy=1、2…、Iy，iz=1、2…、Iz。首先對波數(shù)進(jìn)行Gauss偏移得到偏移頻譜：

(16)

(17)

由式(16)可以看出,Gauss偏移頻譜是在f(xa,yb,zc)標(biāo)準(zhǔn)三維離散傅里葉變換(DFT)的基礎(chǔ)上乘以ei(αix Δkxxa+αiyΔkyyb+αizΔkz)Gauss偏移因子。由(17)式可以看出，只需要在標(biāo)準(zhǔn)三維離散傅里葉逆變換(IDFT)的基礎(chǔ)乘以ei(αix Δkxxa+αiyΔkyyb+αizΔkz)Gauss偏移因子再進(jìn)行Gauss積分求和。

2 模型算例

通過設(shè)計簡單模型，驗(yàn)證3D Gauss-FFT磁異常正演算法的可靠性，并與空間域方法以及標(biāo)準(zhǔn)3D FFT算法進(jìn)行比較。

2.1 單個長方體正演計算

設(shè)計地下三維空間大小為x:0 m～2 000 m，y：0 m～2 000 m，z：0 m～1 000 m，將地下三維模型均分成200×200×100個直立長方體模型單元，每個模型單元的大小為10 m×10 m×10 m，場源所在位置為：x：800 m～1 200 m，y：700 m～1 300 m，z：700 m～800 m，如圖2(a)所示為y=995剖面的磁化率分布圖。場源所在區(qū)域的磁化率大小k=0.01 SI，剩余的模型空間看成背景場，其模型單元的磁化率大小設(shè)置為k=0 SI，正常場大小T0為50 000 nT，磁化方向與正常場方向相同且垂直向下。

圖2 單個長方體模型y=595 m縱切面磁異常及其誤差Fig.2 The magnetic anomalies and errors of a cuboid model computed by 3D Fourier forward modeling along the section of y=595 m(a)磁化率分布模型；(b)模型理論磁異常值；(c)、(e)分別為標(biāo)準(zhǔn)3D FFT正演縱切面磁異常值及其誤差；(d)、(f)分別為4點(diǎn)3D Gauss-FFT正演縱切面磁異常值及其誤差

圖3 單個長方體模型 y=595 m與x=995 m切面交線處磁異常值及其誤差Fig.3 The magnetic anomalies and errors along the profile x=595m and y=995 m of cuboid model(a)空間域及標(biāo)準(zhǔn)3D FFT以及3D Gauss-FFT磁異常值；(b)標(biāo)準(zhǔn)3D FFT以及3D Gauss-FFT磁異常值絕對誤差

圖4 階梯、長方體混合模型y=695 m縱切面磁異常及其誤差Fig.4 The magnetic anomalies and errors of a cuboid and dyke model computed by 3D Fourier forward modeling along the section of y=695 m(a)磁化率分布模型；(b)模型理論磁異常值；(c)、(e)分別為標(biāo)準(zhǔn)3D FFT正演縱切面磁異常值及其誤差；(d)、(f)分別為4點(diǎn)3D Gauss-FFT正演縱切面磁異常值及其誤差

圖2(b)為空間域算法y=595 m縱切面的理論磁異常值，可見其在異常體位置有一個明顯的磁異常。圖2(c)、圖2(e)分別為標(biāo)準(zhǔn)3D FFT算法計算的磁異常值和標(biāo)準(zhǔn)3D FFT算法與空間域理論值的絕對誤差，其絕對誤差最大值為7 nT，離場源越遠(yuǎn)誤差越大，圖2(d)、圖2(f)分別為4點(diǎn)3D Gauss-FFT正演縱切面磁異常值及其與空間域理論值的絕對誤差，其在異常體邊界部分誤差稍有增大，絕對誤差最大值為0.04 nT。通過對比可知，4點(diǎn)3D Gauss-FFT正演算法極大地提高了模型正演精度，相對于標(biāo)準(zhǔn)3D FFT算法提高了幾個數(shù)量級。

圖5 階梯、長方體混合模型y=695 m與x=995 m 切面交線處磁異常值及其誤差Fig.5 The magnetic anomalies and errors along the profile x=695m and y=995 m of dyke and cuboids hybrid model model(a)空間域及標(biāo)準(zhǔn)3D FFT以及3D Gauss-FFT磁異常值；(b)標(biāo)準(zhǔn)3D FFT以及3D Gauss-FFT磁異常值絕對誤差

圖3中通過y=595 m縱切面與x=995 m縱切面交線處空間域正演磁異常，標(biāo)準(zhǔn)3D FFT正演磁異常以及4點(diǎn)3D Gauss-FFT正演磁異常的對比，以及其絕對誤差對比，更直觀地看出4點(diǎn)3D Gauss-FFT正演算法在計算精度上的優(yōu)越性。

2.2 階梯、長方體混合模型正演

為了進(jìn)一步驗(yàn)證3D Gauss-FFT的在較復(fù)雜模型計算精度及計算用時上的優(yōu)越性，以及其在斜磁化條件下，且磁化強(qiáng)度與正常場方向不同時的適用性，設(shè)計一個階梯狀、長方體混合復(fù)雜模型。設(shè)計地下三維空間大小為x：0 m～2 000 m，y：0 m～2 000 m，z：0 m～1 000 m，將地下三維模型均分成200×200×100個直立長方體模型單元，每個模型單元的大小為10 m×10 m×10 m，源體所在位置如圖4(a)所示，表示為y=995 m剖面的磁化率分布圖。階梯狀異常體磁化率大小k1=0.01 SI，長方體磁化率大小k2=0.02 SI，剩余的模型空間看成背景場，其模型單元的磁化率大小設(shè)置為k=0 SI，正常場大小T0為50 000 nT，磁化強(qiáng)度方向磁傾角為45°，磁偏角為0°，正常場方向磁傾角為30°，磁偏角為0°。

圖4、圖5為在斜磁化條件下，且磁化強(qiáng)度方向和正常場方向不同時，階梯、長方體混合模型y=695 m縱切面的磁異常及其誤差，可見3D標(biāo)準(zhǔn)FFT最大絕對誤差17.33 nT，而3D Gauss-FFT正演最大絕對誤差為0.08 nT，相比3D標(biāo)準(zhǔn)FFT依然在計算精度上占有的優(yōu)勢，表明3D Gauss-FFT正演算法的有很好的適用性。

根據(jù)該階梯、長方體混合模型，表1從計算用時，絕對誤差最大值以及內(nèi)存需求三個方面對比空間域算法，標(biāo)準(zhǔn)3D FFT以及3D Gauss-FFT算法。結(jié)果表明，空間域算法計算結(jié)果為理論值，但計算用時巨大；標(biāo)準(zhǔn)3D FFT用時最少，但精度最低，只有3D Gauss-FFT算法既具有較高的精度又節(jié)約了計算時間。2點(diǎn)Gauss-FFT算法耗時少，但精度不算太高，4點(diǎn)和8點(diǎn)Gauss-FFT算法在精度上差不多，但4點(diǎn)Gauss-FFT算法計算用時少，因此4點(diǎn)Gauss-FFT算法一般足以達(dá)到正演的精度要求。

表1 各種正演算法的計算性能比較

測試計算機(jī)的配置為CPU-Intel(R) Core(TM)i5-4590，主頻為3.30 GHz，內(nèi)存為16 GB。

3 結(jié)論

針對空間域和頻率域兩種算法各自的優(yōu)點(diǎn)與不足，引入了高精度3D Gauss-FFT技術(shù)用于磁異常正演，并通過模型驗(yàn)證表明了3D Gauss-FFT磁異常正演在計算時間和計算精度上的絕對優(yōu)勢，4點(diǎn)Gauss-FFT正演相比于空間域算法在計算時間上減少了三個數(shù)量級，但相比于標(biāo)準(zhǔn)3D FFT算法計算時間有所增加；在計算精度上顯著降低了標(biāo)準(zhǔn)3D離散傅里葉變換(DFT)引起的誤差，提高了多個數(shù)量級；但內(nèi)存需求相比于空間域算法稍有增加。并通過與2點(diǎn)及8點(diǎn)Gauss-FFT正演結(jié)果進(jìn)行比較，認(rèn)為4點(diǎn)Gauss-FFT在計算精度及計算效率上已達(dá)到相應(yīng)的正演要求。筆者只研究3D Gauss-FFT磁異常正演，今后可以對磁異常梯度張量、或者其他領(lǐng)域進(jìn)行正演計算，在計算時間上，可以加入并行化計算進(jìn)一步提高計算效率，并可以將該算法應(yīng)用到磁異常三維磁化率反演。