一種相對保幅的低頻逆時噪聲壓制方法及其應用

陳可洋

（中國石油大慶油田有限責任公司勘探開發(fā)研究院，黑龍江 大慶 163712）

在復雜構造高精度地震波成像方法中，逆時偏移是理論較為成熟且最精確的成像方法之一。該方法采用雙程地震波波動方程，對方程的近似最少，適合于任意陡傾角、速度在縱橫向變化較為劇烈的情況，可以解決單程波偏移方法受地層傾角限制的缺陷和克?；舴蚱品椒ǖ亩嘀底邥r和計算盲區(qū)等問題，同時還可實現多次波、回轉反射波等通常認為是干擾波類型波場的準確成像（其他方法均將多次波等波場當作一次反射波進行處理，這必然引入較大誤差），且能夠準確地處理焦散面的多值走時和相位變化問題[1-5]。同時，隨著計算機技術的快速發(fā)展，特別是基于CPU/GPU高性能集群并行計算技術和大容量磁盤的快速存儲技術的出現，較大程度地改善了逆時偏移技術工業(yè)化應用的現狀[6]。鑒于上述諸多優(yōu)點，地震波疊前逆時成像技術受到地球物理學界的高度重視和廣泛應用。

眾所周知，逆時偏移方法區(qū)別于其他現行地震成像方法的顯著特點：經其處理后的結果上存在較強能量的低頻背景噪聲。若不對該噪聲進行有效壓制，則將掩蓋有效的、真實的地層細節(jié)。地球物理界對此開展了多項研究工作，并取得了豐碩的研究成果，主要可歸納為以下兩個方面：一方面對逆時成像條件進行改進。例如：Yoon等[7]提出了在零延遲互相關條件中加入波印廷矢量來消除低頻成像噪聲。Liu等[8]把全波場分解成上下行單程波波場，并運用分解出的單程波波場結合相關逆時成像條件達到消除低頻成像噪聲的目的。陳可洋等[9]首次分析了低頻噪聲成因是由以震源和檢波點及它們在速度分界面的投影為橢圓焦點的三條橢圓曲線引入，上下行波場分離逆時成像條件是消除這種噪聲的有效方法，并指出速度平滑可以削弱低頻背景噪聲，但它以損失成像精度為代價?？惮|等[10]比較了拉普拉斯算子、波印廷矢量法和上下左右行波分解法消除低頻背景噪聲的效果。Antoine等[11]采用最小二乘預測誤差濾波器消除成像噪聲。另一方面是對逆時成像結果進行后續(xù)噪聲壓制處理。例如：Zhang等[12]指出拉普拉斯算子濾波相當于成像波場角度域衰減。劉紅偉等[13]指出在疊前對數據進行相位和振幅校正，在成像后運用拉普拉斯算子濾波法消除成像噪聲。陳可洋[14]研究了低頻噪聲是由上行反射波參與相關逆時成像引起，在疊加域、成像點域、共炮點域進行了拉普拉斯算子去噪分析，同時對比了帶通濾波方法。陳康等[15]采用四階拉普拉斯算子有效地去除了低頻成像噪聲，并使濾波前、后子波的振幅和相位保持相對不變。以上低頻噪聲壓制處理均在深度域實現，低頻背景噪聲均得到了有效壓制。在前人研究基礎上，開展了一種全新的低頻噪聲壓制方法探索研究，它是基于熱力學傳導方程的擴散濾波方法，通過參數試驗和實際資料處理，以達到相對保幅低頻噪聲壓制的目的，對復雜構造地震資料高精度逆時成像具有重要的指導意義。

1 基本原理

在數字圖像處理領域，Perona和Malik[16]首次提出了用于熱力學傳導的偏微分方程，即PM方程：

式中：t為擴散時間；div為散度算子；?是梯度算子；是擴散函數，為一個有界非負的遞減函數；U為t時刻的擴散濾波結果，其中U0為t=0時刻的原始數據，即擴散濾波迭代計算的初始條件。

為了將式（1）應用于低頻逆時噪聲壓制處理，對其進行簡化。令擴散函數為某一常數g0，此時有，然后對其進行有限差分離散化，得到如下數值計算公式：

式（2）為三維擴散濾波方程，其中，i和Δx分別為x方向的離散網格節(jié)點號和網格步長；j和Δy分別為y方向的離散網格節(jié)點號和網格步長；k和Δz分別為z方向的離散網格節(jié)點號和網格步長；n和Δt分別為時間t方向的離散網格節(jié)點和離散步長；N為擴散濾波迭代次數。

對式（2）還可作進一步簡化，此時，令空間步長Δx=Δy=Δz，方程系數，忽略空間網格和時間網格步長及擴散濾波系數的影響，得：

式（2）中僅有一個參數，即擴散濾波迭代次數N?，F定義經擴散濾波N次迭代后的信號變?yōu)閁N，原始信號為U0，兩者的差為Uerr=U0-UN。經擴散濾波后的信號UN是在原始信號U0基礎上進行平滑處理，主要反映原始信號的低頻分量，即低頻背景噪聲，而殘差信號Uerr主要反映原始信號的高頻分量，即掩蓋在低頻背景噪聲之下的有效地層細節(jié)信息。所提方法的相對保幅性將在后續(xù)的實際處理試驗中進行詳細探討。

2 應用實例

2.1 低頻噪聲壓制效果分析

以FZ地區(qū)三維地震資料疊前逆時偏移處理剖面為例，采用3種低頻噪聲壓制方法和2組文中擴散濾波參數（迭代次數N）進行三維數據體去噪試驗和應用效果對比分析。圖1（a）為經照明補償后的疊前逆時偏移疊加剖面，其深度域頻譜如圖2（a）所示。圖1（b1）和圖1（b2）分別為圖1（a）經鄰點加權平滑濾波處理后的結果及分離出的低頻噪聲，對應的深度域頻譜分別為圖2（b1）和圖2（b2）。圖1（c1）和圖1（c2）分別為圖1（a）經文中擴散濾波方法處理后（迭代次數為30次）的結果及分離出的低頻噪聲，對應的深度域頻譜分別為圖2（c1）和圖2（c2）。圖1（d1）和圖1（d2）分別為圖1（a）經文中擴散濾波方法處理后（迭代次數為200次）的結果及分離出的低頻噪聲，對應的深度域頻譜分別為圖2（d1）和圖2（d2）。圖1（e1）和圖1（e2）分別為圖1（a）經高通濾波處理后（斜坡為0.1～2 Hz）的結果及分離出的低頻噪聲，對應的深度域頻譜分別為圖2（e1）和圖2（e2）。

分析圖1（a）可知，逆時偏移剖面上存在較強能量的低頻背景噪聲，掩蓋了有效的地層細節(jié)信息，且主要地震能量集中于頻譜的低頻端（圖2（a）中粉紅色箭頭所示）。從圖1（b1）和圖1（c1）分析可知，這兩種方法的低頻噪聲壓制效果相當，地層細節(jié)特征得到有效刻畫，低頻能量也得到有效壓制（對比圖2（b）和圖2（c）），但在去噪后剖面里均存在一定能量的有效信號的殘留（圖1（b2）和圖1（c2）），其中，在經鄰點加權平滑濾波處理后的剖面里存在更多的高頻信息（圖1粉紅色橢圓）。對比圖1（c1）和圖1（d1）可知，形成這兩種去噪效果差異的原因僅與文中擴散濾波方法的迭代次數大小有關，迭代次數越大，低頻噪聲壓制程度越輕，這一點也可從圖1（c2）和圖1（d2）的殘差剖面對比分析得到。其中，圖1（d2）主要為壓制掉低頻背景逆時噪聲，且無明顯的有效信號殘留。從頻譜分析可知，圖2（d1）與圖2（c1）具有相似的主頻，但前者頻帶更寬，這一點也可從對應擴散濾波剖面的頻譜分析得到，此時圖2（d2）中的頻譜能量要比圖2（c2）更集中于低頻端，因此，保留的高頻信息更多。分析圖1（d1）和圖1（e1）可知，文中擴散濾波方法與高通濾波方法具有相似的噪聲壓制效果。對比頻譜可知，這兩種去噪方法的主頻相當，其中圖2（d1）的頻譜能量分布要比圖2（e1）更廣，頻帶也更寬（圖2鮮綠色橢圓）。分析圖1（d2）和圖1（e2）可知，由于高通濾波方法采用的是逐道處理，該方法處理剖面存在橫向能量不均過度的問題（圖1藍色和紅色橢圓），這正是帶通濾波方法模糊陡傾角成像波場的主要原因，而文中所提的擴散濾波方法則是全局去噪方法。因此，不存在上述問題，且去噪后的地震剖面橫向能量過渡自然，有效信號的保真性和保幅性更好。

圖1 實際地震資料疊前逆時偏移低頻噪聲壓制前后及其殘差剖面Fig.1 Practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression sections and their residual sections

圖2 實際地震資料疊前逆時偏移低頻噪聲壓制前后及其殘差的深度域歸一化頻譜Fig.2 Depth normalized amplitude spectrum of practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression and their residual section spectrum

2.2 相對保幅性分析

為了進一步驗證文中所提方法的相對保幅性，從第2.1小節(jié)3種方法的三維處理數據體滿覆蓋位置提取出一道逆時偏移記錄（圖3）和1 km深度處的水平切片，同時從水平切片中提取出一條滿覆蓋位置的主測線幅值（圖4）進行相對保幅性分析。

圖3 實際地震資料疊前逆時偏移低頻噪聲壓制前后深度方向波形曲線及其殘差對比Fig.3 Wave in the depth direction of practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression and their residual curves comparison

圖3（a）和圖4（a）分別為原始逆時偏移深度域信號波形和水平切片振幅曲線（com，黑色）和經不同去噪方法處理后的有效信號，其中，smo（粉紅色）代表鄰點加權平滑濾波處理結果，dif（藍色）代表擴散濾波方法處理后（迭代次數為30次）的結果，dif1（紅色）代表擴散濾波方法處理后（迭代次數為200次）的結果，bp（深黃色）代表高通濾波處理后（斜坡為0.1～2 Hz）的結果。圖3（b）和圖4（b）分別為原始逆時偏移深度域信號波形和與水平切片振幅曲線不同去噪方法處理后的低頻噪聲，其中com-smo（粉紅色）代表鄰點加權平滑濾波法去除的低頻噪聲，com-dif（藍色）代表擴散濾波方法（迭代次數為30次）去除的低頻噪聲，com-dif1（紅色）代表擴散濾波方法（迭代次數為200次）去除的低頻噪聲，com-bp（深黃色）代表高通濾波法（斜坡為0.1～2 Hz）去除的低頻噪聲。

圖4 實際地震資料疊前逆時偏移低頻噪聲壓制前后水平切片中某inline方向曲線及其殘差對比Fig.4 Wave along an inline direction in horizontal slice of practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression and their residual curves comparison

分析圖3（a）可知：由于原始信號振幅較大，原始逆時偏移結果的波形曲線與經不同去噪方法處理后的結果差異較小，但從去除的低頻噪聲曲線中可以明顯地看出不同去噪方法的壓制效果（圖3（b）），其中文中擴散濾波方法（迭代次數為30次）去噪效果與鄰點加權平滑濾波法效果相當，前者更加平滑，與文中擴散濾波方法（迭代次數為200次）效果相比，迭代次數越大，地層細節(jié)特征將得到越多的保留，低頻波形曲線過度自然光滑，但如果該參數取得過大，則達不到低頻噪聲壓制的目的，此時的去噪結果將保留更多的低頻能量。因此，迭代次數的選擇需根據試驗和處理要求來實現優(yōu)化選擇。而高通濾波法（斜坡為0.1～2 Hz）去除結果存在一定的波形抖動（圖3（b）天藍色橢圓），這與頻域斜坡的選取存在一定的影響，若斜坡選取過陡，則吉普斯效應就會出現在處理結果中。雖然高通濾波方法在橫向相對保幅性方面存在缺陷，但仍是目前低頻噪聲壓制中的一種高效的方法，因其采用逐道處理方式，需求內存較小，不需要多次迭代計算。

分析圖4（a）可知，原始水平切片的振幅曲線與不同去噪方法處理后的結果差異較明顯，去噪后的橫向能量均比去噪前減弱，其中擴散濾波方法（迭代次數為30次）和鄰點加權平滑濾波法的低頻噪聲壓制效果較強，擴散濾波方法（迭代次數為200次）和高通濾波法（斜坡為0.1～2 Hz）較弱，但彼此間的地層細節(jié)差異不明顯。從去除的低頻噪聲水平切片振幅曲線中可以明顯地看出不同去噪方法的壓制效果（圖4（b）），其中文中擴散濾波方法的（迭代次數為30次）去噪效果與鄰點加權平滑濾波法效果相當，但與文中擴散濾波方法（迭代次數為200次）效果相比，后者的地層細節(jié)特征將得到更大程度的保留，其低頻波形曲線過度更光滑且自然，而在高通濾波法（斜坡為0.1～2 Hz）的去除結果里存在一定的波形抖動（圖4（b）綠色橢圓），這與其逐道處理方式有關，導致了處理結果在橫向能量不保幅的問題。

3 結束語

1）引入一種基于熱力學傳導的擴散濾波方法用于有效削弱或消除低頻逆時噪聲，在離散方程系數一定情況下，其簡化計算公式僅受迭代次數一個參數影響，不同迭代次數實現不同低頻噪聲壓制程度的去噪效果。因此，根據參數試驗和處理要求可實現迭代次數的優(yōu)化選擇，同時該方法簡單易實現，且可以推廣應用于二維或一維數據的低頻噪聲壓制中。

2）分析三維地震資料逆時偏移結果的低頻噪聲壓制效果表明：三種方法均可以有效壓制低頻噪聲。文中所提的擴散濾波法具有最佳的相對保幅性，去噪結果的頻帶也更寬，主頻段譜能量更均勻。鄰點加權平滑濾波法的低頻噪聲壓制程度較強，其去噪效果與較小擴散濾波迭代次數情況等效。

3）高通濾波法具有逐道快速低頻噪聲壓制的特點，仍是一種有效的低頻噪聲壓制方法，其低頻噪聲壓制程度可控（由濾波參數控制），但在三維數據體的橫向和縱向上的相對保幅性較差，易受吉普斯效應的影響，同時對陡傾角成像波場存在一定的模糊性。

[1]陸基孟.地震勘探原理（上、下冊）[M].北京：中國石油大學出版社，2004.

[2]凌云.地震數據采集、處理解釋一體化實踐與探索[M].北京：石油工業(yè)出版社，2007.

[3]渥·伊爾馬滋.地震資料分析——地震資料處理、反演和解釋（上、下冊）[M].劉懷山，王克斌，童思友，等譯.北京：石油工業(yè)出版社，2006.

[4]陳可洋.高階彈性波波動方程正演模擬及逆時偏移成像研究[D].大慶：大慶石油學院,2009.

[5]陳可洋.地震波逆時偏移方法研究綜述[J].勘探地球物理進展，2010，33（3）：153-159.

[6]劉紅偉，李博，劉洪，等.地震疊前逆時偏移高階有限差分算法及GPU實現[J].地球物理學報，2010，63（7）：1725-1733.

[7]Yoon K,Marfurt K J.Reverse-time migration using the poynting vector[J].Exploration Geophysics,2006,37(1):102-107.

[8]Liu F Q,Zhang G,Morton S A,et al.Reverse-time migration using one way wavefield imaging condition[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:2170-2174.

[9]陳可洋，吳清嶺，范興才，等.地震波疊前逆時偏移脈沖響應研究與應用[J].石油物探，2013，52（2）：163-170.

[10]康瑋，程玖兵.疊前逆時偏移假象去除方法[J].地球物理學進展，2012，27（3）：1163-1172.

[11]Antoine G,Bruno K,Biondo B.Least square attenuation of reverse time migration artifacts[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006:2348-2352.

[12]Zhang Y,James S.Practical issues of reverse time migration:true amplitude gathers,noise removal and harmonic source encoding[C].CPS/SEG Beijing 2009 International Geophysical Cofference&Exposition,2009:5.

[13]劉紅偉，劉洪，鄒振，等.地震疊前逆時偏移中的去噪與存儲[J].地球物理學報，2010，63（9）：2171-2180.

[14]陳可洋.基于拉普拉斯算子的疊前逆時噪聲壓制方法[J].巖性油氣藏，2011，23（5）：87-95.

[15]陳康，吳國忱.逆時偏移拉普拉斯算子濾波改進算法[J].石油地球物理勘探，2012，47（2）：249-255.

[16]Perona P,Malik J.Scale-space and edge detection using anisotropic diffusion[J].IEEE Tran Pattern Analysis and Machine Intelligence,1990,12(7):629-639.