基于樣條插值的FFT及其在重磁場正演中的應用

周印明 戴世坤 李 昆 何展翔 胡曉穎 王金海

（①中南大學地球科學與信息物理學院，湖南長沙410083；②中國石油集團東方地球物理公司綜合物化探處，河北涿州072751；③中南大學有色金屬成礦預測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室，湖南長沙410083；④南方科技大學前沿與交叉科學研究院，廣東深圳518055；⑤青海省第三地質(zhì)勘查院，青海西寧810029）

0 引言

20世紀60年代Rader［1］提出了一種快速離散傅里葉變換算法，極大提高了離散傅里葉變換的計算效率，使得快速傅里葉變換（FFT）算法得到了更加廣泛的應用，被評為20世紀十大最優(yōu)秀算法之一。早期FFT算法要求數(shù)據(jù)維數(shù)具有2n形式，隨后針對數(shù)據(jù)維數(shù)為3n、2n K（n和K均為正整數(shù)）形式的FFT算法相繼被提出［2-5］。而基于質(zhì)數(shù)分解的FFT算法的出現(xiàn)，解除了對數(shù)據(jù)維數(shù)的要求。Frigo等［6］開發(fā)了FFTW軟件包，集成各類FFT算法，實現(xiàn)了對任意維數(shù)數(shù)據(jù)的離散傅里葉變換快速計算。對于某些問題，空間或者頻率域采樣數(shù)據(jù)往往是不等間隔的，采用該采樣方式的計算方法稱為非規(guī)則采樣數(shù)據(jù)離散傅里葉變換（NDFT，Non-Uniform Discrete Fourier Transform）。Dutt等［7］、Anderson［8］、Nguyen等［9］、Potts等［10］、Fessler等［11］、Greengard等［12］、S?rensen等［13］通過改進FFT算法實現(xiàn)了NDFT。

由于傅里葉變換算子的振蕩性，傅里葉變換可以歸為振蕩函數(shù)的積分。高振蕩函數(shù)積分的計算是計算科學領域的研究熱點。目前求解振蕩函數(shù)積分的方法大致可以分為三類：Filon型方法，Levin型方法和漸進展開法。對于傅里葉變換而言，由于算子e－ikx的矩已知，大多采用Filon型方法實現(xiàn)傅里葉變 換［14-18］。Clendenin［14］采用的是線性函數(shù)；Piessens等［19］采用了更為復雜的Chebyshev多項式；Arieh［20］對基于Filon型方法的傅里葉變換算法中的積分點進行了分析，并在高震蕩方程求解中取得良好的效果；Wu等［21］將高斯積分應用于核函數(shù)的積分中，取得了較高的計算精度；Li等［22］、Dai等［23］和戴世坤等［24］對高斯快速傅里葉變換與傳統(tǒng)傅里葉變換在重、磁三維正演中的應用進行了對比研究。這些方法的區(qū)別在于剖分區(qū)間采用不同形式的多項式逼近核函數(shù)。

本文結合傅里葉變換積分與三次樣條函數(shù)，提出一種快速傅里葉變換方法。該方法特點之一是能根據(jù)變換函數(shù)譜的變化趨勢，任意地選取采樣間距，提高了傅里葉變換剖分采樣的靈活性，能兼顧計算精度和計算效率；特點之二是可根據(jù)離散后單元積分得到解析表達式，克服傳統(tǒng)方法難以獲得解析表達式的問題，進一步提高計算精度。

1 理論與方法

1.1 基于樣條插值的一維傅里葉變換

函數(shù)f（m）的一維傅里葉正變換公式為

式中：k表示波數(shù)；F（k）為波數(shù)譜。對式（1）進行離散化

式中M表示總單元數(shù)。

對單元i的內(nèi)核函數(shù)進行三次樣條插值［25］

式中系數(shù)ai、bi、ci、di可采用基于固定邊界條件的樣條插值計算得到。

將式（3）代入式（2），得到

據(jù)式（5）可得解析積分表達式為

特別地，當波數(shù)k＝0時，單元積分可簡化為

對式（7）積分可得解析表達式

1.2 基于樣條插值的二維傅里葉變換

對函數(shù)f（x，y）進行二維傅里葉正變換

式中：kx為x方向的波數(shù)；ky為y方向的波數(shù)；F（kx，ky）為波數(shù)譜。

式（9）的數(shù)值計算為二維積分，本文采取的策略是依次對x、y兩個方向分別進行一維積分

2 正確性驗證

為了分析誤差，引入相對均方根誤差［26］，一維和二維的相對均方根誤差公式分別為

式中：P和M分別是x和y方向的節(jié)點數(shù)；Bi、Bij分別表示表示一維和二維數(shù)值解；、分別表示一維和二維解析解。相對均方根誤差能突出大異常值所占的比重，避免小異常和過零異常造成的誤差失真。

為了驗證該方法的正確性，對高斯函數(shù)進行樣條插值傅里葉變換，并與解析公式進行對比。其中，一維和二維高斯函數(shù)及其傅里葉變換解析公式分別為

式中α是任意非零常數(shù)，本文令α＝0.001。

設定一維傅里葉變換的取值范圍為［－100m，100m］，采樣點數(shù)為101，波數(shù)采樣范圍為［－0.3，0.3］，采樣點數(shù)為101；二維變換x、y方向點坐標范圍均為［－100m，100m］，采樣點數(shù)均為101個，兩個方向波數(shù)范圍均為［－0.3，0.3］，采樣點數(shù)均為101。一維和二維傅里葉變換的波數(shù)域和空間域采樣均采用等間隔剖分。

圖1為一維正、反傅里葉變換數(shù)值解與解析解曲線。可以看出，正、反傅里葉變換的數(shù)值解與解析解吻合度很高，正、反傅里葉變換的相對均方根誤差分別約為4.5×10－6和4.2×10－7，驗證了本文一維傅里葉正、反變換理論的正確性。

圖2為二維正、反傅里葉變換數(shù)值解與解析解等值線?？梢钥闯?，正、反傅里葉變換的數(shù)值解與解析解吻合度很高，正、反傅里葉變換的相對均方根誤差分別約為6.0×10－6和1.4×10－5，驗證了本文二維傅里葉正、反傅里葉變換理論的正確性。從計算結果可以看出，一維正反變換與二維正反變換的計算精度變化規(guī)律不同，這是由于不同的核函數(shù)譜的變化規(guī)律不同，因此其正、反傅里葉變換的計算精度不存在特定的規(guī)律性。

圖1 一維傅里葉正（左）、反（右）變換數(shù)值解與解析解對比

圖2 二維傅里葉正（左）、反（右）變換數(shù)值解與解析解

3 基于樣條插值傅里葉變換的重磁場三維數(shù)值模擬

在頻率域重磁場數(shù)值模擬中，傅里葉變換是關鍵步驟之一。設計連續(xù)介質(zhì)模型，對基于樣條插值傅里葉變換的重磁場三維數(shù)值模擬的計算效率與精度進行了研究，其中重磁場三維數(shù)值模擬采用空間波數(shù)混合域三維數(shù)值模擬方法［22-24］。文獻［25］表明基于高斯—FFT的重磁場數(shù)值模擬具有較高的計算精度，本文算例將高斯點數(shù)NG＝4的高斯FFT法重磁場數(shù)值模擬結果作為連續(xù)介質(zhì)的解析解。本文算例中，令kx＝ky＝K。測試計算機配置為CPU-Inter Core i5-4590，主頻為3.3GHz，內(nèi)存為16GB。

3.1 重力場數(shù)值模擬

設計連續(xù)介質(zhì)模型，其垂直方向密度不變，水平方向剩余密度ρ的空間變化可表示為

模擬區(qū)域大小為20km（x）×20km（y）×10km（z），坐標原點位于水平范圍中心；異常體區(qū)域大小為20km（x）×20km（y）×3km（z），頂面埋深為3km。模擬區(qū)域水平采樣間隔均為200m，垂向采樣間隔為100m，網(wǎng)格剖分數(shù)為101×101×101；異常體區(qū)域波數(shù)采樣范圍為－1.5×1.0－2～1.5×1.0－2。

圖3為不同波數(shù)下模擬的地面重力異常場及其梯度張量的相對均方根誤差?？梢钥闯觯S著波數(shù)的增大，重力異常場及梯度張量的相對均方根誤差逐漸減小，當波數(shù)K為71時，其計算精度與高斯FFT計算結果接近。當K繼續(xù)增大時，梯度張量的計算誤差有增大趨勢，原因在于隨著波數(shù)的增大，基于樣條插值傅里葉變換的計算精度逐漸高于4點高斯FFT的計算精度，若把高斯FFT作為解析解，誤差緩慢增大。

表1所示為不同K時利用本文方法與高斯FFT方法計算地面重力場及其張量梯度耗時統(tǒng)計?？梢钥闯觯S著K的增大，樣條插值FFT計算時間呈線性增長。在計算精度相近的情況下，即波數(shù)K＝71時，樣條插值FFT的計算時間為0.74s，高斯FFT的計算時間為4.85s，即前者約為后者的1/6?；贜G＝4的高斯FFT計算時間相當于傳統(tǒng)FFT計算量的16倍，而本文算法僅需計算一次積分，因而在相同的計算精度下，基于樣條插值的傅里葉變換計算速度高于高斯FFT算法。

圖3 模擬的重力異常場（左）及梯度張量（右）的相對均方根誤差曲線

表1 重力模型樣條FFT與高斯FFT計算時間統(tǒng)計

圖4為K＝71時地面重力異常場及梯度張量的數(shù)值解，圖5為圖4數(shù)據(jù)的絕對誤差。從圖4和圖5可以看出，重力異常場的絕對誤差最大約為0.05mGal，重力梯度張量的絕對誤差最大約為0.06E，均比數(shù)值解（圖4）低3個數(shù)量級，表明其計算精度高，K選取71是合適的。

3.2 磁場數(shù)值模擬

設計一個連續(xù)介質(zhì)模型，垂直方向磁化率不變，水平方向磁化率κ的空間變化為

模擬區(qū)域為20km（x）×20km（y）×10km（z），坐標原點位于水平范圍中心；異常區(qū)域大小為20km（x）×20km（y）×3km（z），頂面埋深為3km。模擬區(qū)域網(wǎng)格剖分數(shù)為101×101×101；異常區(qū)域空間水平采樣間隔均為200m，垂向采樣間隔為100m。背景磁異常場為35A/m，磁傾角為45°，磁偏角為5°，波數(shù)K采樣范圍為－1.5×1.0－2～1.5×1.0－2。圖6所示為不同K時，采用樣條插值FFT計算的地面磁異常場及其梯度張量的相對均方根誤差?？梢钥闯?，隨著K的增加，磁異常場及梯度張量的相對均方根誤差逐漸減小；當K＝101時，磁異常場及梯度張量的相對均方根誤差小于0.1%，其計算精度與高斯FFT接近，當波數(shù)繼續(xù)增大時，均方根誤差降低的趨勢變緩。

圖4 K＝71時重力異常Gx（a）、Gy（b）、Gz（c）及梯度張量Gxx（d）、Gxy（e）、Gyy（f）、Gzx（g）、Gzy（h）、Gzz（i）的數(shù)值解

圖5 K＝71時重力異常Gx（a）、Gy（b）、Gz（c）及梯度張量Gxx（d）、Gxy（e）、Gyy（f）、Gzx（g）、Gzy（h）、Gzz（i）絕對誤差分布

表2所示為不同K時分別用本文方法及高斯FFT計算地面磁場及其梯度張量所耗時間對比?？梢钥闯觯S著K的增加，樣條插值FFT計算時間呈線性增長。在計算精度相近的情況下，即K＝101時本文方法耗時1.52s，高斯插值FFT算法耗時6.77s，即樣條插值FFT耗時約為高斯插值FFT的1/4。

圖7為波數(shù)K＝101時磁異常場及梯度張量的數(shù)值解，圖8為波數(shù)K＝101時磁異常場及梯度張量的絕對誤差。由圖7與圖8可以看出，磁異常場的絕對誤差最大約為0.05n T，磁梯度張量的絕對誤差最大約為2×10－5n T/m，均比數(shù)值解（圖7）小2個數(shù)量級，表明本文方法計算精度高，且本實驗中K選取101是合適的。

圖6 不同K時本文方法計算的磁異常場（左）及梯度張量（右）相對均方根誤差

表2 磁力數(shù)據(jù)樣條FFT與高斯FFT計算時間對比

圖7 K＝101時本文方法計算的磁異常場Bx（a）、By（b）、Bz（c）及梯度張量Bxx（d）、Bxy（e）、Byy（f）、Bzx（g）、Bzy（h）、Bzz（i）的數(shù)值解

圖8 K＝101時本文方法計算的磁異常場Bx（a）、By（b）、Bz（c）及梯度張量Bxx（d）、Bxy（e）、Byy（f）、Bzx（g）、Bzy（h）、Bzz（i）絕對誤差

4 應用

為了進一步驗證基于樣條插值FFT方法對實測數(shù)據(jù)的應用效果，筆者利用“地理空間數(shù)據(jù)云”平臺下載了安徽省數(shù)字高程數(shù)據(jù)（經(jīng)度為117.6°～118.4°，緯度為30.5°～31.2°），對此數(shù)據(jù)進行數(shù)值模擬。該區(qū)地形如圖9所示，境內(nèi)中部丘陵、崗地起伏，呈北東向展布。假設該地區(qū)地下地層為水平層狀介質(zhì)［26］，采用基于樣條插值FFT的空間波數(shù)混合域三維數(shù)值模擬方法計算其重力異常場，結果如圖10所示?？梢钥闯?，重力高的區(qū)域與圖9所示起伏地形范圍吻合很好，說明基于樣條插值FFT方法應用于實際數(shù)據(jù)的處理是可靠的，具有較好的適應性。

圖10 本文方法計算的重力異常場圖

5 結論

針對傳統(tǒng)FFT存在截斷效應的問題，本文結合FFT與三次樣條插值算法，提出了一種基于樣條插值的FFT算法。該方法充分利用樣條插值擬合的高階連續(xù)性及采樣靈活的優(yōu)點，兼顧計算精度和計算效率，實現(xiàn)了高效、高精度的FFT。得到如下結論：

（1）利用高斯函數(shù)FFT，對比了數(shù)值解與解析解，驗證了基于樣條插值算法的一維、二維FFT的正確性。

（2）將算法應用于頻率域重、磁位場的計算，設計連續(xù)介質(zhì)模型，對比基于樣條插值FFT與高斯插值FFT的數(shù)值解，前者的計算精度更高。

（3）對于連續(xù)介質(zhì)模型，在計算精度相近的情況下，利用基于樣條插值FFT進行重、磁場模擬計算，耗時統(tǒng)計數(shù)據(jù)表明，樣條FFT的計算效率明顯高于高斯FFT，前者的計算時間約為后者的1/6（重力數(shù)據(jù)）和1/4（磁力數(shù)據(jù)）。