預(yù)條件共軛梯度法在地震數(shù)據(jù)重建方法中的應(yīng)用

霍志周，熊 登，張劍鋒

1 中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

2 中國科學(xué)院大學(xué)，北京 100049

3 中國石油東方地球物理公司物探技術(shù)研究中心，河北涿州 072751

1 引 言

由于地表障礙物的存在（建筑物、道路、橋梁等）、地形條件的因素（禁采區(qū)和山區(qū)、森林、河網(wǎng)地區(qū)等）以及儀器硬件等問題的存在，引起地震數(shù)據(jù)在空間方向上的不規(guī)則采樣．地震數(shù)據(jù)的不規(guī)則對DMO、FK 濾波、速度分析、多次波消除、譜估計和波動方程偏移等方法的處理結(jié)果帶來了嚴重的影響．因此研究非規(guī)則采樣地震數(shù)據(jù)的規(guī)則化對提高后續(xù)速度分析的精度、改善疊加效果、提升疊前偏移的質(zhì)量等有著非常重要的作用．基于最小平方反演的Fourier變換地震數(shù)據(jù)重建方法是目前一種常用的地震重建方法［1－6］，本文通過將Fourier變換與貝葉斯參數(shù)反演方法相結(jié)合建立目標(biāo)函數(shù)［7－9］，得到最小范數(shù)解．該方法的核心是對每一個頻率成分求解線性方程組，從不規(guī)則采樣地震數(shù)據(jù)中估計出重建數(shù)據(jù)的Fourier系數(shù)．可以使用共軛梯度法（Conjugate Gradient Method，CG）求解反演生成的線性方程組，由于該線性方程組的系數(shù)矩陣是Toeplitz矩陣，因此可以使用預(yù)條件共軛梯度法（Preconditioned Conjugate Gradient Method，PCG）［8，10］加快收斂速度，提高計算效率，保證解的穩(wěn)定性．本文主要論述基于Toeplitz矩陣的循環(huán)預(yù)條件的構(gòu)造方法［11－20］，以及在最小平方Fourier地震數(shù)據(jù)重建方法中的應(yīng)用．

2 最小平方Fourier重建原理

通常地震數(shù)據(jù)在時間方向上是規(guī)則采樣的，首先將地震數(shù)據(jù)從時間域變換到頻率域．對于沿空間方向非規(guī)則采樣的地震數(shù)據(jù)，為了求取空間波數(shù)，定義非規(guī)則離散傅里葉變換為

式中，Δxn＝（xn＋1－xn－1）/2，n＝0，…，N－1，xn表示空間采樣點位置．如果直接用（1）式計算空間波數(shù)，則由于采樣的非規(guī)則而引起誤差．因此需要利用最小平方反演來消除誤差．設(shè)帶限地震數(shù)據(jù)的波數(shù)范圍為 ［－MΔk，MΔk］ ，相應(yīng)的波數(shù)域帶寬為（2M＋1）Δk，波數(shù)域規(guī)則采樣采樣間隔為Δk，則由波數(shù)域重建任意空間位置xn的離散傅里葉反變換為

把上述離散傅里葉反變換當(dāng)作正演模型．記系數(shù)矩陣為Anm＝，不規(guī)則采樣地震數(shù)據(jù)為dn＝P（xn，ω），待求的規(guī)則采樣波數(shù)為ω）．則將公式（2）寫成矩陣形式為

在實際數(shù)據(jù)處理中，由于數(shù)據(jù)可能不完全是帶限的，所以部分空間波數(shù)成分會超出定義的頻帶范圍，這些超出的成分構(gòu)成了上述正演模型的誤差或噪音，因此在（3）式中需要噪聲項

Duijndam 等［7］通過最小平方反演估計得到非規(guī)則采樣數(shù)據(jù)d（xn，ω）的波數(shù)．從非規(guī)則采樣數(shù)據(jù)向量d中計算出未知的規(guī)則采樣的傅里葉系數(shù)向量可以歸結(jié)為求解一個不適定線性反演問題．可以使用任何所需的參數(shù)估計技術(shù)，我們認為噪音和先驗信息都是高斯分布的，假設(shè)噪音的協(xié)方差矩陣為Cn，平均值為零：n＝N（0，Cn），先驗信息為＝．可以利用貝葉斯參數(shù)反演方法通過尋找后驗概率密度函數(shù)的最大值來進行反演．后驗概率密度函數(shù)為

其中p（d｜）是似然函數(shù)，p（）表示模型向量的先驗分布．

求取p（d｜）的最大后驗解轉(zhuǎn)化為求下面目標(biāo)函數(shù)的最小化的解

這里，Cn為噪音的協(xié)方差矩陣，C~p為先驗?zāi)Ｐ偷膮f(xié)方差矩陣．最小化目標(biāo)函數(shù)得

式中，為要計算得到的規(guī)則采樣波數(shù)，Wii＝Δxi為權(quán)系數(shù)對角矩陣，AH和A分別表示Fourier正反變換矩陣，λ為高斯正則化因子．方程（7）左端的系數(shù)矩陣AHWA為Toeplitz矩陣，矩陣AHWA的病態(tài)程度受非規(guī)則采樣數(shù)據(jù)之間最大間隔Δxa的大小來控制．最大間隔Δxa越大，矩陣AHWA病態(tài)程度越大，方程就越難以收斂．求解方程（7）最常用的方法為CG 算法，其中主要的運算是矩陣和向量相乘．因為系數(shù)矩陣為Toeplitz 矩陣，可以充分利用Toeplitz矩陣的特殊結(jié)構(gòu)，使用Toeplitz矩陣與向量相乘的循環(huán)褶積快速算法［14］，通過快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）計算循環(huán)褶積來大大加快運算，該方法稱為CG＋FFT 算法．很明顯共軛梯度法迭代次數(shù)越多，解方程耗費的時間就越多，計算效率越低．因此需要引入預(yù)條件來改善系數(shù)矩陣的病態(tài)程度，以減少迭代次數(shù)．但引入預(yù)條件又為迭代法增加了額外的計算量，因此施加預(yù)條件需要把握預(yù)條件減少的解方程計算量和預(yù)條件本身增加的計算量之間的平衡，形成預(yù)條件還要考察預(yù)條件本身需要的計算量和內(nèi)存需求．總體而言，預(yù)條件必須解決如下兩方面問題：

1）預(yù)條件本身應(yīng)該很容易構(gòu)造和使用，即不需要花費大量計算時間構(gòu)造預(yù)條件；

2）預(yù)條件系統(tǒng)應(yīng)該很容易求解，即迭代收斂迅速．

預(yù)條件需要針對具體問題設(shè)計，不同的預(yù)條件有不同的特征值聚集能力．常規(guī)的Jacobi預(yù)條件和SSOR 預(yù)條件對我們的問題已經(jīng)沒有加速效果．因此我們轉(zhuǎn)而尋求能滿足如下條件的預(yù)條件：適用于致密矩陣計算，并且能充分利用Toeplitz矩陣的結(jié)構(gòu)，能很好滿足上述要求的就是循環(huán)預(yù)條件．循環(huán)預(yù)條件的最大優(yōu)勢就是利用FFT 計算預(yù)條件循環(huán)矩陣的特征值，這樣在共軛梯度法計算矩陣與向量的除法時用特征值與向量FFT 簡單相除即可．

3 Toeplitz矩陣的循環(huán)預(yù)條件

3．1 矩陣的預(yù)條件

對于求解大型線性方程組Ax＝b，加上預(yù)條件矩陣M，得到原方程組同解的方程組

預(yù)條件方程組滿足關(guān)系式

即矩陣M－1A的條件數(shù)遠小于矩陣A的條件數(shù)．如果M對A的近似程度非常高，則通常M－1A的特征值接近1，譜半徑為：ρ（M－1A）≤1，迭代計算時方程（8）就會迅速收斂．一般而言，實際計算中M－1A的特征值在某點（一般為1．0）附近集中程度越高，則求解方程（8）收斂速度越快．在共軛梯度迭代中加預(yù)條件時，并不需要計算矩陣乘積M－1A，而是在每次迭代時求解另外一個線性方程

其中rk＋1和zk＋1分別為共軛梯度第k＋1次迭代時的殘差向量．

公式（10）計算時，若預(yù)條件M為循環(huán)矩陣，則可以利用循環(huán)褶積快速算法．當(dāng)系數(shù)矩陣A為Toeplitz矩陣時，其循環(huán)預(yù)條件通常采用一定的規(guī)則從Toeplitz矩陣A中提取元素，并組合為循環(huán)矩陣，然后利用循環(huán)矩陣與向量相除的循環(huán)褶積快速算法，因此計算效率比大部分其它預(yù)條件都高．

3．2 常規(guī)循環(huán)預(yù)條件

考慮到通常Toeplitz矩陣的主對角線和其旁邊的元素為強占優(yōu)，可以簡單的把前1/2列拷貝過來，再卷繞回去，作為循環(huán)預(yù)條件，稱為Strang循環(huán)預(yù)條件［15］．循環(huán)預(yù)條件還可以通過使

達到最小化來定義，這里‖·‖F(xiàn)表示Frobenius范數(shù)，稱為T Chan循環(huán)預(yù)條件［16－17］．

通常Strang 預(yù)條件矩陣本身條件數(shù)大于T Chan預(yù)條件．另外Strang預(yù)條件不能保持原系數(shù)矩陣的正定性，而T Chan 則能保證正定．雖然如此，但Strang預(yù)條件收斂速度并不一定比T Chan預(yù)條件慢．T Chan循環(huán)預(yù)條件由于能很好保持原系數(shù)矩陣的正定性，所以很常用．但由于T Chan循環(huán)預(yù)條件在實際應(yīng)用中仍不能取得滿意效果，所以通常不直接使用T Chan預(yù)條件，而是用來構(gòu)造其它類型的循環(huán)預(yù)條件．

4 生成函數(shù)與循環(huán)預(yù)條件

4．1 生成函數(shù)

從數(shù)學(xué)角度可用生成函數(shù)來刻畫Toeplitz矩陣的特征值分布情況．n×n階的Toeplitz矩陣A的第j個特征值可以用生成函數(shù)f（x）在 ［－π，π］ 的取值f（2πj/n）來表示．令ak為Toeplitz矩陣A的第k個元素，aj，k＝aj－k，0≤j，k＜n，n≥1，則定義f為Toeplitz矩陣A的生成函數(shù)，有

根據(jù)矩陣生成函數(shù)理論，良態(tài)Toeplitz矩陣的生成函數(shù)f（x）＞0，即為非負函數(shù)，由于特征值f（2πj/n）＞0恒成立，所以Toeplitz矩陣正定．病態(tài)Toeplitz矩陣是由于生成函數(shù)有零值．生成函數(shù)的零值越多，則Toeplitz矩陣病態(tài)越嚴重．相應(yīng)地，使預(yù)條件后的系數(shù)矩陣具有更良態(tài)的譜，也等同于使該系數(shù)矩陣對應(yīng)生成函數(shù)的零值減少或者完全沒有．

病態(tài)Toeplitz矩陣的循環(huán)預(yù)條件C，從生成函數(shù)角度可以看作是用另一個正函數(shù)g（x）與A的生成函數(shù)f（x）進行褶積，以此對A矩陣生成函數(shù)中的零值進行匹配，同時保證褶積后的函數(shù)能收斂到f（x），即滿足

g（x）＊（x）作為C的生成函數(shù)，為f（x）的平滑近似，并且在xj＝2πj/n，（0≤j＜n）的取值應(yīng)非常接近，這樣能充分保證預(yù)條件矩陣C與原矩陣A非常接近．即有

滿足上述關(guān)系后預(yù)條件系統(tǒng)C－1A的特征值就接近1，使預(yù)條件后的方程組迭代時能迅速收斂．滿足該要求的g（x）均與δ函數(shù)很接近，如后文的高階B樣條核函數(shù)．用生成函數(shù)和一些核函數(shù)（如Dirichlet核、Fejér核）褶積，可以生成多種循環(huán)預(yù)條件，包括常見的Strang預(yù)條件、T Chan預(yù)條件等．其中Strang預(yù)條件可以看作是Dirichlet核與生成函數(shù)褶積后生成的，而T Chan預(yù)條件則能用Fejér核和生成函數(shù)褶積后生成．生成函數(shù)的優(yōu)點是能夠比較直觀和真實的反映矩陣的病態(tài)程度．但根據(jù)生成函數(shù)與核函數(shù)褶積構(gòu)造的預(yù)條件，在處理特別病態(tài)情況也即是原始矩陣生成函數(shù)有多個零值時，效果仍然不理想，因此還需要深入研究．另一個很關(guān)鍵的問題是，在構(gòu)造生成函數(shù)預(yù)條件時，如果每次都要求推導(dǎo)Toeplitz矩陣對應(yīng)的生成函數(shù)，則不但構(gòu)成極大的計算負擔(dān)，而且操作復(fù)雜，效率很低．因此，利用生成函數(shù)構(gòu)造預(yù)條件的方法主要集中于不需要顯式寫出Toeplitz矩陣的生成函數(shù)，而是通過Toeplitz矩陣的元素與核函數(shù)作用就能直接構(gòu)造預(yù)條件的方法．

4．2 廣義Jackson核函數(shù)循環(huán)預(yù)條件

廣義Jackson 核循環(huán)預(yù)條件［18］Kn，2r是用Toeplitz矩陣的生成函數(shù)與廣義Jackson 核函數(shù)褶積，然后用褶積后的新函數(shù)作為生成函數(shù)，求取對應(yīng)的循環(huán)矩陣作為原Toeplitz矩陣A的循環(huán)預(yù)條件．但是廣義Jackson 核預(yù)條件方法不需要顯式求出Toeplitz矩陣的生成函數(shù)．令x∈ ［－π，π］，定義Kn，2r（x）核函數(shù)為

其中r為褶積核的階數(shù)．Kn，2r為歸一化常數(shù)，滿足Kn，2r（x）dx＝1．廣義Jackson 核函數(shù)作為循環(huán)預(yù)條件的優(yōu)點是對所有正的r，n，都能保證生成的預(yù)條件是正定的．另外，高階預(yù)條件可以由低階預(yù)條件與自身褶積得到，這樣構(gòu)造預(yù)條件相對簡單．

給定m×m階Toeplitz矩陣A，其對應(yīng)的廣義Jackson循環(huán)預(yù)條件為m×m循環(huán)矩陣則有，即r（n－1）＜m≤rn．

核函數(shù)與生成函數(shù)褶積

其中

由上式可見，Kn，2r＊f僅僅依賴于元素ak．為了構(gòu)造循環(huán)預(yù)條件，僅需要Kn，2r＊f在，0≤j＜m處的值，而不需要知道生成函數(shù)信息．

4．3 B樣條循環(huán)預(yù)條件

B樣條循環(huán)預(yù)條件［19］是將Toeplitz矩陣A的生成函數(shù)與B 樣條函數(shù)褶積，將褶積結(jié)果作為生成函數(shù)，求取對應(yīng)的循環(huán)矩陣作為預(yù)條件．B樣條預(yù)條件方法也不需知道Toeplitz矩陣A的生成函數(shù)，而是直接提取A的元素aij與核函數(shù)進行組合．

為了構(gòu)造循環(huán)預(yù)條件，定義B樣條函數(shù)為

式中，cm為歸一化常數(shù)，使Bm（0）＝1．則Bm（x）關(guān)于中心對稱，．對任意m＝1，2，…，m階函數(shù)Qm（x）在x＝0，1，…，m結(jié)點上定義為

基于Bm（x）的離散函數(shù)定義為

因此基于Bm（x）的B 樣條循環(huán)預(yù)條件Cm定義為（只寫出循環(huán)矩陣Cm的首列）

其中ak和－k為Toeplitz矩陣A中的元素．由此可見B樣條預(yù)條件并不需要生成函數(shù)，而是直接與A中的元素作用．高階B 樣條預(yù)條件通常能更好適應(yīng)病態(tài)Toeplitz方程組，但階數(shù)高意味著預(yù)條件計算量的增加，可能會比較明顯的增加求解方程組的時間，實際使用也一般就取到4階．B樣條函數(shù)具有很好的對稱性，而且階數(shù)越高，系數(shù)的函數(shù)圖形越接近δ函數(shù)，對應(yīng)B樣條預(yù)條件矩陣與原矩陣A更接近，從而有更好的迭代加速效果．

4．4 其它核函數(shù)構(gòu)造的循環(huán)預(yù)條件

沿用生成函數(shù)與非負核函數(shù)褶積構(gòu)造新的循環(huán)預(yù)條件的思路，Chan［20］給出了幾種核函數(shù)并推導(dǎo)了相應(yīng)的循環(huán)預(yù)條件．我們選用了效果相對較好的兩種循環(huán)預(yù)條件，對應(yīng)的預(yù)條件矩陣C首列為

1）Hamming 循環(huán)預(yù)條件：

2）von Hann 循環(huán)預(yù)條件：

5 數(shù)值算例

本文提供兩類數(shù)值算例，以驗證所述循環(huán)預(yù)條件的收斂速度和處理病態(tài)Toeplitz矩陣的能力．所有的算例均采用預(yù)條件共軛梯度迭代法求解方程．迭代中止均采用誤差限‖r‖2/‖b‖2≤1．0×10－7．通過觀察需要的迭代次數(shù)多少，可以看出該預(yù)條件對病態(tài)矩陣的適應(yīng)能力，以及不同循環(huán)預(yù)條件對同一問題收斂能力的相互對比．

5．1 生成函數(shù)構(gòu)造病態(tài)Toeplitz矩陣

第一類數(shù)值算例從純數(shù)學(xué)角度出發(fā)，根據(jù)給定的生成函數(shù)構(gòu)造對應(yīng)的病態(tài)Toeplitz 矩陣，將Toeplitz矩陣作為方程Ax＝b的A矩陣，b向量全設(shè)為1．令x∈ ［－ π，π］ ，則給定生成函數(shù)x2（π4－x4），對應(yīng)的病態(tài)Toeplitz矩陣為T1，對應(yīng)的元素為

表1為不同預(yù)條件求解T1x＝b需要的迭代次數(shù)．T1的行列大小見左邊第一列．從第二列開始依次為以單位矩陣I為預(yù)條件（無預(yù)條件）、1 階到4階廣義Jackson 核函數(shù)循環(huán)預(yù)條件（Kn，2，Kn，4，Kn，6，Kn，8）、3階到5階B樣條循環(huán)預(yù)條件（B3，B4，B5），以及Hamming和von Hann核函數(shù)循環(huán)預(yù)條件對應(yīng)的迭代次數(shù)．其中Kn，2對應(yīng)T Chan循環(huán)預(yù)條件．如果迭代次數(shù)達到n次，共軛梯度還未收斂，則以“－”替代．由表1可見，預(yù)條件能大幅減少迭代次數(shù)．并且當(dāng)n從32增加到512時，矩陣T1條件數(shù)也相應(yīng)從4．1433×102增加到1．0154×105，但預(yù)條件迭代次數(shù)基本不隨n的增加而增加．所以預(yù)條件對求解系數(shù)矩陣為病態(tài)Toeplitz 矩陣的線性方程組具有非常重要的作用．

表1 不同預(yù)條件求解T1x＝b需要的迭代次數(shù)Table 1 Iteration numbers for T1x＝bfor different preconditioned systems

圖1為方程組T1x＝b當(dāng)n＝32時的矩陣特征值分布，圖中橫坐標(biāo)表示特征值模的大小，縱坐標(biāo)沒有特殊含義，當(dāng)圖中有多個預(yù)條件系統(tǒng)時，表示預(yù)條件系統(tǒng)的編號．對比可見，各種核函數(shù)循環(huán)預(yù)條件預(yù)后系統(tǒng)的特征值（圖1（b—d））均有效地壓縮了原病態(tài)Toeplitz矩陣T1特征值分布空間（圖1a），預(yù)條件后的特征值系統(tǒng)均聚集在1．0左右．而且隨著階數(shù)增加，特征值在1．0左右聚集越多，從理論上講需要的迭代次數(shù)應(yīng)該更少，但實際求解需要的迭代次數(shù)并沒有顯著減少，這可能和相互間特征值分布差異仍然不夠明顯有關(guān)．

5．2 地震數(shù)據(jù)模型算例

為模擬地震資料處理中面臨的病態(tài)反演情況，設(shè)計了一組隨機分布的一維坐標(biāo)來模擬個N非規(guī)則采樣地震道的偏移距x．一維非規(guī)則采樣地震數(shù)據(jù)的加權(quán)Fourier重建需要計算如下Toeplitz矩陣T構(gòu)成的線性方程組：

其中Δk＝2π/X為波數(shù)域采樣間隔，X為排列長度，λ為Gauss正則化因子，取λ＝0．01，w為權(quán)系數(shù)．T的病態(tài)程度與排列的非規(guī)則采樣程度相關(guān)．實驗中對規(guī)則排列的xl（間距Δx）加上隨機擾動量，計算中選取xl＝xl＋3Δx（0．5－rand（））．同樣令b＝（1，1，…，1），求解方程組Tx＝b．取N＝32，64，128，256，512，計算迭代次數(shù)如表2所示．

表2 模擬Fourier重建非規(guī)則采樣的反演迭代次數(shù)表Table 2 Iteration number by simulating Fourier reconstruction for irregular sampling

由表2可見，在處理該問題時，預(yù)條件還是有一定的效果．實際計算時發(fā)現(xiàn)即使Toeplitz矩陣大小達到了512×512，Toeplitz矩陣的特征值也很好的被壓縮在2．0 以內(nèi)．但是迭代次數(shù)上的優(yōu)化則沒有理論數(shù)據(jù)效果明顯．而且與理論數(shù)據(jù)正好相反的是，高階核函數(shù)反而沒有低階核函數(shù)的收斂效果好．因此對該類問題的預(yù)條件還需要進一步的探索．

圖2為n＝32時方程（30）的特征值分布，圖中橫縱坐標(biāo)和含義與圖1中相同．可見預(yù)條件確實將病態(tài)Toeplitz矩陣的特征值進行了大幅壓縮．但也可以看出高階核函數(shù)在實際問題中沒有起到更好的預(yù)條件效果，高階和低階核函數(shù)的預(yù)條件系統(tǒng)特征值分布差異很小．經(jīng)計算發(fā)現(xiàn)，B 樣條由3 階到5階，預(yù)條件系統(tǒng)條件數(shù)由66．1 上升至95．9；廣義Jackson 預(yù)條件系統(tǒng)階數(shù)由2 階到8 階，條件數(shù)由67．6上升至118．9，接近原始病態(tài)條件數(shù)159．4．因此高階預(yù)條件不適用該問題的解決．

6 地震數(shù)據(jù)重建算例

首先我們選取SMAART 理論模型數(shù)據(jù)來計算，如圖3a所示的一個單炮地震記錄，總共有348道，其中有57 道是空道，圖3b 為用最小平方Fourier方法重建的結(jié)果，可以看到地震數(shù)據(jù)被很好重建．圖4中我們抽取其中的一道數(shù)據(jù)進行了重建前后的比較，可以看出重建幾乎完全恢復(fù)了原始的數(shù)據(jù)．我們對某油田的一個實際二維共偏移距地震數(shù)據(jù)進行了重建，該數(shù)據(jù)總共有2278 道，其中有684是空道．為了便于清楚的顯示地震剖面，圖5a所示的是其中選取的300 道數(shù)據(jù)，時間軸從1．6s到3．6s．圖5b為用最小平方Fourier方法重建的結(jié)果，可以看到重建道與原始道具有很好的連續(xù)性．圖6a為重建前時間偏移的結(jié)果，圖6b為重建后時間偏移的結(jié)果，對比兩幅圖可以看出，重建后的偏移剖面信噪比更高，構(gòu)造更加清晰，同相軸的連續(xù)性也更好．表3給出了預(yù)條件共軛梯度法與幾種常用的求解方法所需計算時間的比較，可以看出預(yù)條件共軛梯度法在求解速度上具有一定的優(yōu)勢．

圖5 （a）含有空道的共偏移距地震數(shù)據(jù)和（b）最小平方Fourier反演的重建結(jié)果Fig．5 （a）Common offset seismic data containing null traces（a）and its reconstruction by the least squares Fourier inversion（b）

圖6 （a）重建前數(shù)據(jù)的時間偏移結(jié)果和（b）重建后數(shù)據(jù)的時間偏移結(jié)果Fig．6 Result of time migration before data reconstruction（a）and result of time migration after data reconstruction（b）

表3 不同求解線性方程組方法的計算時間比較Table 3 Comparison of computing time for different methods of solving linear equations

預(yù)條件采用3階B 樣條循環(huán)預(yù)條件，以上算例運算結(jié)果均在Interl（R），Core（TM）2，Quad CPU Q9550 2．83GHz，內(nèi)存2GB的微機上計算得到，操作系統(tǒng)為Centos4．7，編程語言為C．

7 結(jié) 論

本文針對最小平方Fourier地震數(shù)據(jù)重建方法所形成的線性方程組的求解問題，論述了Toeplitz矩陣循環(huán)預(yù)條件的基本理論，深入探討了根據(jù)生成函數(shù)構(gòu)造針對病態(tài)Toeplitz矩陣的循環(huán)預(yù)條件的方法，給出了2種核函數(shù)的循環(huán)預(yù)條件構(gòu)造方式．用預(yù)條件共軛梯度法求解方程組，保證解的穩(wěn)定性和提高了收斂速度．在最后我們給出了幾組數(shù)值算例，驗證了我們構(gòu)造的預(yù)條件對理論數(shù)據(jù)具有明顯效果，對實際數(shù)據(jù)也具有一定的效果．

（References）

［1］ Schonewille M A．Fourier reconstruction of irregularly sampled seismic data［Ph．D．thesis］．Delft：Delft University of Technology，2000．

［2］ 劉喜武，劉洪，劉彬．反假頻非均勻地震數(shù)據(jù)重建方法研究．地球物理學(xué)報，2004，47（2）：299－305．

Liu X W，Liu H，Liu B．A study on algorithm for reconstruction of de－alias uneven seismic data．ChineseJ．Geophys．（in Chinese），2004，47（2）：299－305．

［3］ 孟小紅，劉國峰，周建軍．大間距地震數(shù)據(jù)重建方法研究．地球物理學(xué)進展，2006，21（3）：687－691．

Meng X H，Liu G F，Zhou J J．The study of reconstruction of large gap seismic data．ProgressinGeophys．（in Chinese），2006，21（3）：687－691．

［4］ 孟小紅，郭良輝，張致付等．基于非均勻快速傅里葉變換的最小二乘反演地震數(shù)據(jù)重建．地球物理學(xué)報，2008，51（1）：235－241．

Meng X H，Guo L H，Zhang Z F，et al．Reconstruction of seismic data with least squares inversion based on nonuniform fast Fourier transform．ChineseJ．Geophys．（in Chinese），2008，51（1）：235－241．

［5］ Zwartjes P，Gisolf A．Fourier reconstruction with sparse inversion．GeophysicalProspecting，2007，55（2）：199－221．

［6］ Zwartjes P M，Sacchi M D．Fourier reconstruction of nonuniformly sampled，aliased seismic data．Geophysics，2007，72（1）：V21－V32．

［7］ Duijndam A J W，Schonewille M A，Hindriks C O H．Reconstruction of band－limited signals，irregularly sampled along one spatial direction．Geophysics，1999，64（2）：524－538．

［8］ 熊登．疊前地震數(shù)據(jù)規(guī)則化、重建及噪音壓制［博士論文］．北京：中國科學(xué)院研究生院，2008．

Xiong D．Prestack seismic data regularization，reconstruction and noise suppression［Ph．D．thesis］（in Chinese）．Beijing：Graduate University，Chinese Academy of Sciences，2008．

［9］ 高建軍，陳小宏，李景葉等．基于非均勻Fourier變換的地震數(shù)據(jù)重建方法研究．地球物理學(xué)進展，2009，24（5）：1741－1747．

Gao J J，Chen X H，Li J Y，et al．Study on reconstruction of seismic data based on nonuniform Fourier transform．ProgressinGeophys．（in Chinese），2009，24（5）：1741－1747．

［10］ 李冰，劉洪，李幼銘．三維地震數(shù)據(jù)離散光滑插值的共軛梯度法．地球物理學(xué)報，2002，45（5）：691－699．

Li B，Liu H，Li Y M．3－D seismic data discrete smooth interpolation using conjugate gradient．ChineseJ．Geophys．（in Chinese），2002，45（5）：691－699．

［11］ 梅金順，劉洪．ω循環(huán)型邊界條件．地球物理學(xué)報，2003，46（6）：835－841．

Mei J S，Liu H．ω－circulant boundary condition．ChineseJ．Geophys．（in Chinese），2003，46（6）：835－841．

［12］ 梅金順，劉洪．Toeplitz方程組的近似計算．地球物理學(xué)進展，2003，18（1）：128－133．

Mei J S，Liu H．Approximate computation of Toeplitz systems of equations．ProgressinGeophys．（in Chinese），2003，18（1）：128－133．

［13］ 梅金順，劉洪．預(yù)條件方程組及其應(yīng)用．地球物理學(xué)報，2004，47（4）：718－722．

Mei J S，Liu H．Preconditioned equation sets and their applications．ChineseJ．Geophys．（in Chinese），2004，47（4）：718－722．

［14］ Chan R H，Ng N K．Conjugate gradient methods for Toeplitz systems．SIAMReview，1996，38（3）：427－482．

［15］ Strang G．Introduction to Applied Mathematics．Cambridge：Wellesley－Cambridge Press，1986．

［16］ Chan T F．An optimal circulant Preconditioner for Toeplitz systems．SIAMJ．Sci．Stat．Comput．，1988，9（4）：766－771．

［17］ Chan R H．Circulant preconditioners for Hermitian Toeplitz systems．SIAMJ．MatrixAnal．Appl．，1989，10（4）：542－550．

［18］ Chan R H，Yip A M，Ng M K．The best circulant preconditioners for hermitian Toeplitz systems．SIAMJ．Numer．Anal．，2001，38（3）：876－896．

［19］ Chan R H，Tso T M，Sun H W．Circulant preconditioners from B－splines．Proceedings to the SPIE Symposium on Advanced Signal Processing：Algorithms，Architectures，and Implementations VII，1997，3162：338－347．

［20］ Chan R H， Yeung M C． Circulant preconditioners constructed from kernels．SIAMJ．Numer．Anal．，1992，29（4）：1093－1103．