不規則采樣數據廣義線性反演算法的改進

劉保童，李桂花

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

引言

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

1 DFT線性反演理論

定義連續正向空間傅里葉變換為

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

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

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

式中Δxn定義為

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

為已知的非均勻數據

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

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

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

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

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

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

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

2 反演解的改進

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

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

這種非迭代的直接處理方法能夠將輸入數據的單色波分解聚焦到其最大的譜成分上.類似的思想也可用于改進Radon變換的分辨率.[6]

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

3 數值計算試驗

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

圖1 理論合成記錄

4 結論

圖2 已知的不規則采樣數據

圖3 在缺失數據的位置處填充零道

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

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

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

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

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

[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]張賢達.矩陣分析與應用[M].北京:清華大學出版社,2004.

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

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