基于Gauss-FFT技術高精度三維磁異常ΔT正演

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

(1.中南大學 地球科學與信息物理學院，長沙 410083；2.有色資源與地質災害探查湖南省重點實驗室，長沙 410083；3 桂林理工大學 地球科學學院，桂林 541006)

0 引言

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

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

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

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

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

1 方法原理

1.1 3D傅里葉正演理論

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

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

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

(1)

其中：u0為真空磁導率，

(2)

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

(3)

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

(4)

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

(5)

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

(6)

其中k為波數，

根據傅里葉變換微分定理：

(7)

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

(lkx+mky+nkz)·

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

(8)

對整個場源區域模型單元體磁異常頻譜進行累加求和到：

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)]進行三維傅里葉逆變換從而得到整個場源區域磁異常。

1.2 3D Gauss-FFT技術

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

為了便于計算，我們用g(xa,yb,zc)表示各單元體中心的坐標位置，x、y、z方向的觀測間隔為Δx、Δy、Δz，點S為第一個單元體的中心坐標為(x0,y0,z0)，則觀測數據的離散坐標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)

對于場源區域被均分M×N×L棱柱體，在尼奎斯特頻率約束條件下的離散波數可以表示為[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)是采用矩形積分來近似替代連續傅里葉變換，式(13)具有與式(10)相同的形式，為了提高三維傅里葉變換的精度，采用Gauss積分的方法代替矩形積分以提高三維離散傅里葉變換的精度。一維的Gauss積分公式如下：

(15)

其中:h表示Gauss積分節點數，Fh和th表示在[-1,1]區間上的高斯系數和高斯節點值。

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

(16)

(17)

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

2 模型算例

通過設計簡單模型，驗證3D Gauss-FFT磁異常正演算法的可靠性，并與空間域方法以及標準3D FFT算法進行比較。

2.1 單個長方體正演計算

設計地下三維空間大小為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剖面的磁化率分布圖。場源所在區域的磁化率大小k=0.01 SI，剩余的模型空間看成背景場，其模型單元的磁化率大小設置為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)分別為標準3D FFT正演縱切面磁異常值及其誤差；(d)、(f)分別為4點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)空間域及標準3D FFT以及3D Gauss-FFT磁異常值；(b)標準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)分別為標準3D FFT正演縱切面磁異常值及其誤差；(d)、(f)分別為4點3D Gauss-FFT正演縱切面磁異常值及其誤差

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

圖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)空間域及標準3D FFT以及3D Gauss-FFT磁異常值；(b)標準3D FFT以及3D Gauss-FFT磁異常值絕對誤差

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

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

為了進一步驗證3D Gauss-FFT的在較復雜模型計算精度及計算用時上的優越性，以及其在斜磁化條件下，且磁化強度與正常場方向不同時的適用性，設計一個階梯狀、長方體混合復雜模型。設計地下三維空間大小為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，剩余的模型空間看成背景場，其模型單元的磁化率大小設置為k=0 SI，正常場大小T0為50 000 nT，磁化強度方向磁傾角為45°，磁偏角為0°，正常場方向磁傾角為30°，磁偏角為0°。

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

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

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

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

3 結論

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