二維演算法

❶ 二維實序列的快速傅里葉變換（FFT）

在地球物理數據處理中，經常遇到處理二維實數據的情況。例如在地震勘探中，對面波勘探數據作頻散分析解釋時，要將時間－空間域的信息轉換為頻率－波數域頻譜；在重磁異常的濾波或轉換中，要將空間域的異常f（x，y）轉換為波數域F（ω，υ）等。這些分析都需要進行二維的傅里葉變換（FFT）。

根據傅里葉變換的定義，對於連續二維函數f（x，y），其傅里葉變換對為

地球物理數據處理基礎

對於離散的二維序列f_jk（j＝0，1，…，M－1；k＝0，1，…，N－1），其傅里葉變換為

地球物理數據處理基礎

1.二維復序列的FFT演算法

對於M條測線，每條測線N個測點，構成復序列y_jk（j＝0，1，…，M－1；k＝0，1，…，N－1），根據離散傅里葉公式（8－41），其傅里葉變換為

地球物理數據處理基礎

於是，可以分兩步套用一維復FFT完成二維復FFT的計算。

（1）沿測線方向計算

對於j＝0，1，…，M－1逐測線套用一維復FFT，執行式（8－43）。定義復數組 則演算法為

1）對於j＝0，1，…，M－1，作2）～7）；

2）將y_jk輸入A₁（k），即A₁（k）＝y_jk（k＝0，1，…，N－1）；

3）計算W^r，存入W（r），即

4）q＝1，2，…，p（p＝log₂N），若q為偶數執行6），否則執行5）；

5）k＝0，1，2，…，（2^p－q－1）和n＝0，1，2，…，（2^q－1－1）循環，作

A₂（k2^q＋n）＝A₁（k2^q－1＋n）＋A₁（k2^q－1＋n＋2^p－1）

A₂（k2^q＋n＋2^q－1）＝［A₁（k2^q－1＋n）－A₁（k2^q－1＋n＋2^p－1）］·W（k2^q－1）

至k，n循環結束；

6）k＝0，1，2，…，（2^p－q－1）和n＝0，1，2，…，（2^q－1－1）循環，作

A₁（k2^q＋n）＝A₂（k2^q－1＋n）＋A₂（k2^q－1＋n＋2^p－1）

A₁（k2^q＋n＋2^q－1）＝［A₂（k2^q－1＋n）－A₂（k2^q－1＋n＋2^p－1）］·W（k2^q－1）

至k，n循環結束；

7）q循環結束，若p為偶數，將A₁（n）輸入到Y_jn，否則將A₂（n）輸入到Y_jn（n＝0，1，…，N－1）；

8）j循環結束，得到Y_jn（j＝0，1，…，M－1；n＝0，1，…，N－1）。

（2）垂直測線方向計算

對於n＝0，1，…，N－1逐一套用一維復FFT，執行式（8－44）。即

1）對於n＝0，1，…，N－1，作2）～7）；

2）將Y_jn輸入A₁（j），即A₁（j）＝Y_jn（j＝0，1，…，M－1）；

3）計算W^r存入W（r），即

4）q＝1，2，…，p（p＝log₂M），若q為偶數執行6），否則執行5）；

5）j＝0，1，2，…，（2^p－q－1）和m＝0，1，2，…，（2^q－1－1）循環，作

A₂（j2^q＋m）＝A₁（j2^q－1＋m）＋A₁（j2^q－1＋m＋2^p－1）

A₂（j2^q＋m＋2^q－1）＝［A₁（j2^q－1＋m）－A₁（j2^q－1＋m＋2^p－1）］·W（j2^q－1）

至j，m循環結束；

6）j＝0，1，2，…，（2^p－q－1）和m＝0，1，2，…，（2^q－1－1）循環，作

A₁（j2^q＋m）＝A₂（j2^q－1＋m）＋A₂（j2^q－1＋m＋2^p－1）

A₁（j2^q＋m＋2^q－1）＝［A₂（j2^q－1＋m）－A₂（j2^q－1＋m＋2^p－1）］·W（j2^q－1）

至j，m循環結束；

7）q循環結束，若p為偶數，將A₁（m）輸入到Y_mn，否則將A₂（m）輸入到Y_mn（m＝0，1，…，M－1）；

8）n循環結束，得到二維復序列的傅氏變換Y_mn（m＝0，1，…，M－1；n＝0，1，…，N－1），

所求得的Y_mn是復數值，可以寫為

Y_mn＝R_mn＋iI_mn （m＝0，1，…，M－1；n＝0，1，…，N－1）

其中，R_mn，I_mn的值也是已知的。

2.二維實序列的FFT演算法

對於二維的實序列，我們把其看作是虛部為零的復序列，套用上述的二維復序列FFT方法來求其頻譜演算法上也是可行的，但勢必會增加大量的無功運算。因此，有必要研究二維實序列FFT的實用演算法，同一維實序列FFT的實現思路一樣，同樣把二維實序列按一定的規律構造成二維復序列，調用二維復序列FFT，然後通過分離和加工得到原實序列的頻譜。

例如采樣區域有2 M條測線，每條測線有N個點，並且M，N都是2的整數冪，需要計算實樣本序列x_jk（j＝0，1，2，…，2 M－1；k＝0，1，2，…，N－1）的傅氏變換：

地球物理數據處理基礎

類似於一維實序列FFT的思想，直接建立下面的二維實序列FFT演算法：

（1）將一個二維實序列按偶、奇線號分為兩個二維子實序列，分別作為實部和虛部組合為一個二維復序列。即令

地球物理數據處理基礎

（2）調用二維復FFT過程，求出y_jk的二維傅氏變換Y_mn的復數值：

地球物理數據處理基礎

式中：R_mn，I_mn是Y_mn的實部和虛部。

（3）利用R_mn，I_mn換算X_mn的值。

前兩步容易實現，下面分析第（3）步的實現。

記h_jk，g_jk的傅氏變換為H_mn，G_mn。根據傅里葉變換的定義，我們導出X_mn與H_mn，G_mn的關系式：

地球物理數據處理基礎

式中，H_mn，G_mn為復數，我們用上標r和i表示其實部和虛部，將上式右端實部、虛部分離

地球物理數據處理基礎

其中：

地球物理數據處理基礎

下面的任務是將H_mn，G_mn各分量與通過二維復FFT求出的R_mn，I_mn值聯系起來。為此先給出奇、偶分解性質和類似於一維情況的三個二維傅氏變換性質：

（1）奇偶分解性

任何一個正負對稱區間定義的函數，均可唯一地分解為如下偶（even）、奇（odd）函數之和：

地球物理數據處理基礎

（2）周期性

地球物理數據處理基礎

（3）復共軛性

地球物理數據處理基礎

現在我們來建立R_mn，I_mn與H_mn，G_mn的關系。對Y_mn作奇偶分解：

地球物理數據處理基礎

根據線性性質

地球物理數據處理基礎

對照式（8－54）和式（8－55），得

地球物理數據處理基礎

由於h_jk，g_jk是實函數，根據復共軛性質，上面兩式對應的奇偶函數相等。即

地球物理數據處理基礎

再由奇偶分解性和周期性，得

地球物理數據處理基礎

將式（8－57）代入式（8－50），得

地球物理數據處理基礎

再利用H_mn，G_mn周期性及復共軛性，可以得到m＝M/2＋1，…，M－1；n＝0，1，…，N－1的傅氏變換，即

地球物理數據處理基礎

將式（8－50）中M，N改為M－m，N－n，並將上式代入，得

地球物理數據處理基礎

由式（8－58）、式（8－59）和式（8－61）即可得到原始序列x_jk（j＝0，1，…，2M－1；n＝0，1，…，N－1）在m＝0，1，…，M－1；n＝0，1，…，N－1區間的傅氏變換X_mn。

具體二維實序列的FFT演算法如下：

（1）令h_jk＝x_2j，k，g_jk＝x_2j＋1，k，形成

y_jk＝h_jk＋ig_jk （j＝0，1，…，2 M－1；n＝0，1，…，N－1）

（2）調用二維復序列FFT過程，即從兩個方向先後調用一維復FFT演算法式（8－43）和式（8－44），求得y_jk的二維傅氏變換Y_mn的復數值：

Y_mn＝R_mn＋iI_mn （m＝0，1，…，M－1；n＝0，1，…，N－1）

（3）用下列公式由R_mn，I_mn的值換算X_mn的值：

地球物理數據處理基礎

地球物理數據處理基礎