劉保童，李桂花

（1.天水師范學(xué)院 物理與信息科學(xué)學(xué)院，甘肅 天水 741001；2.山東科技大學(xué) 地質(zhì)科學(xué)與工程學(xué)院，山東 青島 266590）

引言

譜分析是一個非常活躍的研究領(lǐng)域，可能的應(yīng)用相當廣泛，詳細的介紹可從參考文獻[1]中看到.從本質(zhì)上講，譜分析是一個欠定的線性反問題，它可用于從非均勻采樣點重建信號，這時，算法從非均勻采樣信號估計離散傅里葉變換（DFT），最后用反DFT得 到 均 勻 采 樣 的 信 號.Oldenburg（1976）把Backus-Gilbert線性反演理論應(yīng)用于具有噪聲和空缺數(shù)據(jù)的傅里葉變換，Sacchi et al.（1998）提出了一種利用貝葉斯（Bayesian）方法的傅里葉變換反演方法，針對稀疏分布傅里葉系數(shù)的估計，他們在迭代反演中使用了先驗柯西（Cauchy）概率密度函數(shù).在DFT反演地震道重建方面，Duijndam et al.（1999）利用了一個與空間采樣間隔相對應(yīng)的對角規(guī)則化矩陣使不規(guī)則采樣數(shù)據(jù)的DFT反演穩(wěn)定化，這是在數(shù)據(jù)空間中使用約束.[2]Wang（2003）利用由初始的DFT估計所計算的模型的協(xié)方差矩陣作為對角加權(quán)矩陣，給出了一種稀疏約束最小平方反演方法，所不同的是在模型空間中使用約束.[3]本文從Moore-Penrose廣義逆的角度討論DFT譜估計這一線性反演問題，在算法中引入了一個實對角正定約束矩陣，用低頻分量的傅里葉圖像對高頻傅里葉圖像的估計進行約束，提高了模型空間的分辨率.作者在微型計算機上編寫程序?qū)崿F(xiàn)了這一算法，用有空缺道的不規(guī)則采樣多道記錄模型進行數(shù)值計算試驗，證明該方法是有效的.

1 DFT線性反演理論

定義連續(xù)正向空間傅里葉變換為

式中x是空間變量，kx是波數(shù)，而ω=2πf，f是時間頻率.反變換由下式給出

時間傅里葉變換對用類似的方式定義，只不過指數(shù)所用符號的約定不同.對于沿x規(guī)則采樣的帶限數(shù)據(jù)，眾所周知，與（1）式有關(guān)的離散傅里葉變換（DFT）為

其中Δx的選擇應(yīng)足夠小以避免空間傅里葉域假頻.在不規(guī)則采樣的情況下，一種獲得傅里葉域數(shù)據(jù)的直接方法是使用黎曼和（the Riemann sum），（1）式中的積分用一個與實際采樣位置（x0，x1，…，xN-1）對應(yīng)的和來代替：

式中Δxn定義為

ln=(xn+xn-1)2是兩個樣點之間的中心點.我們稱變換（4）為非均勻離散傅里葉變換（NDFT），NDFT不能精確地還原重現(xiàn)傅里葉譜.

為已知的非均勻數(shù)據(jù)

現(xiàn)在考慮采樣間隔為Δkx，數(shù)據(jù)帶限區(qū)間在[-MΔkx，MΔkx]總共有Mp=2M+1個樣點的一個規(guī)則采樣傅里葉域，對任一空間位置xn，與NDFT對應(yīng)的離散反傅里葉變換由下式給出

且是精確的，Δkx的選擇應(yīng)足夠小以避免空間x上的假頻.將不規(guī)則空間采樣位置 (x0，x1，…，xN-1)所對應(yīng)的方程（6）的N個式子結(jié)合起來可寫成矩陣形式

利用（7）式求P～是一個線性反演問題，存在解的非唯一性.

由于A是一個N×Mp的矩陣，A的逆應(yīng)該用Moore-Penrose廣義逆來求解.[4]當N＞Mp，即方程（7）為超定方程時，（7）式的最小方差解

當N＜Mp，即方程(7)為欠定方程時，(7)式的最小范數(shù)解

式中λ為阻尼因子，一般取0.1～1即可.(AAH+λ2I)-1等價于一個反褶積算子，它提高了變換域的分辨率.但（8）、（9）式給出的這個反演解還不夠理想，分辨率仍不能滿足信噪分離和數(shù)據(jù)重建的要求，因此，尋求反問題（7)式的具有更高分辨率的解顯得非常重要.

事實上，方程（7）的求解是一個不適定的反問題，因而只能通過加先驗約束來實現(xiàn)，且一般都要迭代.[1]本文下面給出一種非迭代的反演解，可有效地提高變換分辨率，是對（8）、（9）的改進.

2 反演解的改進

由前面的分析可知，反問題（7）是一個約束的欠定最小平方問題，因此，按照廣義最小平方理論[5]，本文給出（7）的一個解

式中，G=(gil)Mp×Mp是一個實對角正定約束矩陣，其對角元素為

這種非迭代的直接處理方法能夠?qū)⑤斎霐?shù)據(jù)的單色波分解聚焦到其最大的譜成分上.類似的思想也可用于改進Radon變換的分辨率.[6]

一旦使用改進的反演方法獲得了尖銳的DFT譜，我們就可以對它實現(xiàn)反傅里葉變換產(chǎn)生數(shù)據(jù)空間中期望的規(guī)則采樣輸出.

3 數(shù)值計算試驗

圖1是一個理論合成二維數(shù)據(jù)體，共60道，道間距20米，時間采樣率2毫秒，有三個同相軸.為了檢驗本文提出的方法的有效性，將圖1中的13道記錄數(shù)據(jù)去除，即假設(shè)已知圖2所示的缺失不規(guī)則采樣數(shù)據(jù)，我們要由它來重建恢復(fù)規(guī)則采樣的完全數(shù)據(jù).為看清楚缺失道所在的位置，在缺失數(shù)據(jù)的位置處填充零道，見圖3.如用普通的最小平方反演解（9）式來重建，所得到的結(jié)果顯示在圖4中，與圖3中的記錄比較很容易看出，這一結(jié)果實際上等價于在缺失樣點位置處插入了零道.這和理想結(jié)果（圖1中的記錄）相差甚遠.現(xiàn)在用本文提出的算法，由圖2中的記錄來重建非均勻數(shù)據(jù)，利用（10）式估計的F-K譜顯示在圖6中，而利用非均勻離散傅里葉變換（NDFT）所得到的圖2中記錄的F-K譜如圖5所示，顯然，改進后的算法使估計出的DFT譜具有很高的分辨率.利用圖6中尖銳的F-K譜，通過反傅里葉變換就可以產(chǎn)生數(shù)據(jù)空間中期望的規(guī)則采樣輸出，重建結(jié)果顯示在圖7中，它與圖1中顯示的期望輸出符合的很好.

圖1 理論合成記錄

4 結(jié)論

圖2 已知的不規(guī)則采樣數(shù)據(jù)

圖3 在缺失數(shù)據(jù)的位置處填充零道

圖4 普通的最小平方反演解重建結(jié)果

圖5 用NDFT所得到的圖2中記錄的F-K譜

不規(guī)則采樣或缺失數(shù)據(jù)的傅里葉重建是數(shù)據(jù)正則化處理的主要方法之一，本文從Moore-Penrose廣義逆出發(fā)，采用最小平方線性反演方法估計DFT譜，在解中引入了一個實對角正定約束矩陣對普通的最小平方算法做了改進，提高了模型空間的分辨率.

圖6 由改進的線性反演解所得到的圖2中記錄的F-K譜

圖7 用本文提出的改進反演解重建所得結(jié)果

[1]SACCHI M D,ULRYCH T J,and WALKER C J.Interpolation and extrapolation using a high-resolution discrete Fourier transform[J].IEEE Transactions on Signal Processing,1998,46(1):31-38.

[2]DUIJNDAM A J W,SCHONEWILLE M A,and HINDRIKS C O H.Reconstruction of band-limited signals,irregularly sampled along one spatial direction[J].Geophysics,1999,64(2):524-538.

[3]YANGHUA Wang.Sparseness-constrained least-squares inversion:Application to seismic wave reconstruction[J].Geophysics,2003,68(5):1633-1638.

[4]張賢達.矩陣分析與應(yīng)用[M].北京:清華大學(xué)出版社,2004.

[5]TARANTOLA A.Inverse problem theory:Methods for data fitting and model parameter estimation[M].Elsevier Science Publ.,1987.

[6]劉保童,朱光明.一種頻率域提高Radon變換分辨率的方法[J].西安科技大學(xué)學(xué)報,2006,26(1):112-116.