f－x域Cadzow濾波隨機噪聲衰減及其局限性分析

彭更新 （中石油塔里木油田分公司勘探開發(fā)研究院，新疆 庫爾勒841000）

鄧曉東 （中國石油大學(xué) （北京）油氣資源與探測國家重點實驗室，北京102249）

滿益志 （中石油塔里木油田分公司勘探開發(fā)研究院，新疆 庫爾勒841000）

劉昭 （中國石油大學(xué) （北京）油氣資源與探測國家重點實驗室，北京102249）

段文勝，崔永福 （中石油塔里木油田分公司勘探開發(fā)研究院，新疆 庫爾勒841000）

頻率－空間域 （f－x域）Cadzow濾波，也稱為f－x域eigenimage濾波，是一種基于SVD分解，利用特征剖面重構(gòu)信號壓制噪聲的方法［1］。該方法最早由Cadzow引入，用于解決譜分析問題，并在核磁共振醫(yī)學(xué)圖像處理中得到廣泛應(yīng)用［2，3］。Trickett將該方法應(yīng)用于地震噪聲壓制，并推廣到f－x－y三維隨機噪聲衰減中［4，5］。為適應(yīng)彎曲地震同相軸的處理，崔樹果等［6］、袁三一等［7］提出了局部Cadzow濾波，以便能最大限度地壓制噪聲，保留有效信號。

奇異值分解 （SVD）去噪是工業(yè)界常用的隨機噪聲壓制方法，它利用信號和噪聲在奇異值分布上的差異進行信號恢復(fù)和噪聲剔除，當?shù)卣鸱瓷錇樗椒瓷渫噍S，該方法效果很好，但不適應(yīng)多組傾斜同相軸的去噪處理。f－x域Cadzow濾波基于信號在f－x域的可預(yù)測性，將地震數(shù)據(jù)變換到f－x域，取每一個頻率切片構(gòu)建Hankel矩陣，再進行SVD分解，所得到的奇異值對應(yīng)于所有不同斜率的同相軸，所以f－x域Cadzow濾波更適合于處理傾斜同相軸。筆者主要通過多組數(shù)值模型，對Cadzow濾波的去噪性能及其局限性進行了試驗分析。

1 基本原理

假設(shè)地震數(shù)據(jù)中含K條線性同相軸和隨機噪聲頻譜，將其變換到f－x域后，表示為：

式中：T（f，x）為頻率 －空間域地震數(shù)據(jù)；f為頻率，Hz；x是地震記錄的水平位置，m；Wm（f）為子波的振幅譜；τm和lm分別為第m條線性同相軸的截距和斜率；K為線性同相軸的個數(shù)；N（f，x）為噪聲的頻譜。

設(shè)頻率f為常數(shù)，抽取一個頻率切片，構(gòu)建Hankel矩陣，表示成：

式中：S、N分別為信號與噪聲的n×（u－n＋1）Hankel矩陣。所謂的Hankel矩陣是指沿反對角線方向的矩陣元素相等，通常取n＝［u／2］使該矩陣盡可能為方陣，n＝［·］表示向上舍入。Cadzow濾波就是要利用SVD分解重構(gòu)信號矩陣S。

假設(shè)地震信號有k個傾角，根據(jù)信號在f－x域的預(yù)測性特點，f－x域某一固定頻率的信號采樣點可以表示成：

式中：Pd（d＝1，2，…，k）為預(yù)測濾波器；n為地震信號的道數(shù)。式（3）表明，Sj可以由前j－1個點預(yù)測出來。那么，由Sj構(gòu)成的Hankel矩陣S經(jīng)過矩陣的元素運算后求得的有效秩數(shù)將不超過傾角的個數(shù)k，也就是說rank（S）≤k。

對矩陣S進行SVD分解，得到：

式中：U、V為酉矩陣；H為共軛轉(zhuǎn)置；Σ為實對角線矩陣，其對角線元素σi＝ ［Σ］i，i，且有排列順序σ1≥σ2≥ …≥σr≥0。矩陣的有效秩數(shù)不超過傾角的個數(shù)，因而用降秩近似和：Fk（S）＝I1＋I2＋…＋Ik，Ii＝（i ＝ 1，2，…，k）。式中：Fk（S）為矩陣S的降秩近似和；Ii為第i個加權(quán)特征圖像；ui和vi分別為矩陣U和V的第i列。由于用此方式估計的矩陣Fk（S）一般不滿足Hankel結(jié)構(gòu)，需要沿反對角線方向作平均，表 示 成 Hk（S）。Stephenson［3］證 明，對于含k個傾角的無噪地震數(shù)據(jù)，有Hk（S）＝S。

圖1描述了利用Cadzow濾波技術(shù)進行信號恢復(fù)的過程。隨著秩數(shù)k的增加，信號恢復(fù)程度越高；當秩數(shù)等于同相軸傾角個數(shù)（k＝3）時，信號被完全恢復(fù)出來。

圖1 Cadzow濾波信號恢復(fù)過程

2 模型試驗

圖2（a）是由1個水平同相軸和1個傾斜同相軸構(gòu)成的地震記錄；圖2（b）是加入20%隨機噪聲后結(jié)果；圖2（c）是利用f－x域Cadzow濾波去噪后的結(jié)果，隨機噪聲得到了較好壓制；圖2（d）為去噪后的殘差剖面，殘差剖面中沒有剩余相干能量，說明該方法具有較好的保幅性能。

圖2 f－x域Cadzow濾波隨機噪聲衰減

3 局限性分析

圖3（a）是一個含傾斜同相軸的地震記錄，在同相軸下方有一個能量較強的突發(fā)噪聲干擾。圖 3 （b） 是 f－x 域Cadzow濾波去噪后的結(jié)果，可以看出，Cadzow濾波之后，雖然強突發(fā)噪聲的能量得到了較好的壓制，但其剩余能量沿傾斜同相軸方向擴散，形成明顯的去噪假象，這種現(xiàn)象有可能給后續(xù)地震解釋工作造成 “陷阱”。

圖4（a）是一個含斷層的水平同相軸理論模型；圖4（b）為加入50%隨機噪聲后結(jié)果；圖4（c）為f－x域Cadzow濾波去噪后的地震記錄。由圖4可以看出，雖然Cadzow濾波較好地壓制了隨機噪聲，但兩個斷點明顯向兩側(cè)延伸，模糊了斷點精度。這種現(xiàn)象在實際復(fù)雜斷塊地震資料處理中是十分有害的，需要引起地震資料處理人員的高度重視。

圖5是含雙曲同相軸的地震記錄及其Cadzow濾波后的結(jié)果。由于雙曲同相軸的兩翼接近于線性同相軸，去噪之后地震反射畸變不大，但在雙曲同相軸的頂點附近，有效信號被嚴重壓制和畸變。該試驗表明，基于線性預(yù)測理論的Cadzow濾波去噪方法不能很好地適應(yīng)曲率較大彎曲同相軸的去噪處理。

圖3 強能量突發(fā)噪聲地震剖面及f－x域Cadzow濾波的結(jié)果

圖4 含斷層的地震記錄及f－x域Cadzow濾波去噪結(jié)果

4 結(jié) 語

雖然在處理線性同相 軸 時， f－x 域Cadzow濾波具有較強的去噪能力、較好的保幅性能，但由于該方法的理論基礎(chǔ)是基于地震信號的線性可預(yù)測性，因此，在應(yīng)用中也存在一定的局限性：當?shù)卣鹩涗浿写嬖趶娔芰客话l(fā)噪聲時，f－x 域 Cadzow 濾波之后，強噪聲能量沿地震信號方向擴散，產(chǎn)生虛假的反射同相軸；f－x域Cadzow濾波模糊了斷點位置，在對復(fù)雜斷塊疊后地震資料進行去噪處理時，應(yīng)該關(guān)注和警惕該方法對斷點產(chǎn)生的模糊效應(yīng)；在地震記錄中存在曲率變化較大的彎曲同相軸時，f－x域Cadzow濾波很容易損害大曲率有效信號的反射能量。

圖5 含雙曲同相軸的地震記錄及其Cadzow濾波后的結(jié)果

［1］Trickett S R．F－xeigenimage noise suppression ［A］．2002SEG Annual Meeting ［C］．Salt Lake，2002－10－6～11．

［2］Cadzow J A．Signal enhancement——a composite property mapping algorithm ［J］．IEEE Transactions on Ascoustics Speech an Signal Processing，1988，36（1）：49～62．

［3］Stephenson D S．Linear prediction and entropy maximum methods in NMR spectroscopy ［M］．Oxford：Pergamon Press，1988．

［4］Trickett S．F－xyeigenimage noise suppression ［J］ ．Geophysics，2003，68（2）：751～759．

［5］Trickett S．F－xy Cadzow noise suppression ［A］．2008SEG Annual Meeting ［C］．Las Vegas，2008－11－9～14．

［6］崔樹果，朱凌燕，王建花 ．f－x域Cadzow技術(shù)分塊壓制隨機噪聲及其應(yīng)用 ［J］．石油物探，2012，51（1）：43～50．

［7］Yuan S Y，Wang S X．A local f－xCadzow method for noise reduction of seismic data obtained in complex formations ［J］．Petroleum Science，2011，8 （3）：269～277．