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

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

0 引言

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

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

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

1 理論與方法

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

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

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

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

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

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

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

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

特別地，當(dāng)波數(shù)k＝0時(shí)，單元積分可簡(jiǎn)化為

對(duì)式（7）積分可得解析表達(dá)式

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

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

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

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

2 正確性驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 磁場(chǎng)數(shù)值模擬

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

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

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

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

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

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

圖6 不同K時(shí)本文方法計(jì)算的磁異常場(chǎng)（左）及梯度張量（右）相對(duì)均方根誤差

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

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

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

4 應(yīng)用

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

圖10 本文方法計(jì)算的重力異常場(chǎng)圖

5 結(jié)論

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

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

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

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