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

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

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

0 引言

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

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

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

1 理論與方法

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

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

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

式中M表示總單元數。

對單元i的內核函數進行三次樣條插值［25］

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

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

據式（5）可得解析積分表達式為

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

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

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

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

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

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

2 正確性驗證

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

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

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

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

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

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

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

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

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

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

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

3.1 重力場數值模擬

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

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

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

表1所示為不同K時利用本文方法與高斯FFT方法計算地面重力場及其張量梯度耗時統計。可以看出，隨著K的增大，樣條插值FFT計算時間呈線性增長。在計算精度相近的情況下，即波數K＝71時，樣條插值FFT的計算時間為0.74s，高斯FFT的計算時間為4.85s，即前者約為后者的1/6。基于NG＝4的高斯FFT計算時間相當于傳統FFT計算量的16倍，而本文算法僅需計算一次積分，因而在相同的計算精度下，基于樣條插值的傅里葉變換計算速度高于高斯FFT算法。

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

表1 重力模型樣條FFT與高斯FFT計算時間統計

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

3.2 磁場數值模擬

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

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

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

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

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

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

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

表2 磁力數據樣條FFT與高斯FFT計算時間對比

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

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

4 應用

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

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

5 結論

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

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

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

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