陳可洋

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

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

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

1 基本原理

在數(shù)字圖像處理領(lǐng)域，Perona和Malik[16]首次提出了用于熱力學(xué)傳導(dǎo)的偏微分方程，即PM方程：

式中：t為擴(kuò)散時(shí)間；div為散度算子；?是梯度算子；是擴(kuò)散函數(shù)，為一個(gè)有界非負(fù)的遞減函數(shù)；U為t時(shí)刻的擴(kuò)散濾波結(jié)果，其中U0為t=0時(shí)刻的原始數(shù)據(jù)，即擴(kuò)散濾波迭代計(jì)算的初始條件。

為了將式（1）應(yīng)用于低頻逆時(shí)噪聲壓制處理，對(duì)其進(jìn)行簡(jiǎn)化。令擴(kuò)散函數(shù)為某一常數(shù)g0，此時(shí)有，然后對(duì)其進(jìn)行有限差分離散化，得到如下數(shù)值計(jì)算公式：

式（2）為三維擴(kuò)散濾波方程，其中，i和Δx分別為x方向的離散網(wǎng)格節(jié)點(diǎn)號(hào)和網(wǎng)格步長(zhǎng)；j和Δy分別為y方向的離散網(wǎng)格節(jié)點(diǎn)號(hào)和網(wǎng)格步長(zhǎng)；k和Δz分別為z方向的離散網(wǎng)格節(jié)點(diǎn)號(hào)和網(wǎng)格步長(zhǎng)；n和Δt分別為時(shí)間t方向的離散網(wǎng)格節(jié)點(diǎn)和離散步長(zhǎng)；N為擴(kuò)散濾波迭代次數(shù)。

對(duì)式（2）還可作進(jìn)一步簡(jiǎn)化，此時(shí)，令空間步長(zhǎng)Δx=Δy=Δz，方程系數(shù)，忽略空間網(wǎng)格和時(shí)間網(wǎng)格步長(zhǎng)及擴(kuò)散濾波系數(shù)的影響，得：

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

2 應(yīng)用實(shí)例

2.1 低頻噪聲壓制效果分析

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

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

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

圖2 實(shí)際地震資料疊前逆時(shí)偏移低頻噪聲壓制前后及其殘差的深度域歸一化頻譜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 相對(duì)保幅性分析

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

圖3 實(shí)際地震資料疊前逆時(shí)偏移低頻噪聲壓制前后深度方向波形曲線及其殘差對(duì)比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）分別為原始逆時(shí)偏移深度域信號(hào)波形和水平切片振幅曲線（com，黑色）和經(jīng)不同去噪方法處理后的有效信號(hào)，其中，smo（粉紅色）代表鄰點(diǎn)加權(quán)平滑濾波處理結(jié)果，dif（藍(lán)色）代表擴(kuò)散濾波方法處理后（迭代次數(shù)為30次）的結(jié)果，dif1（紅色）代表擴(kuò)散濾波方法處理后（迭代次數(shù)為200次）的結(jié)果，bp（深黃色）代表高通濾波處理后（斜坡為0.1～2 Hz）的結(jié)果。圖3（b）和圖4（b）分別為原始逆時(shí)偏移深度域信號(hào)波形和與水平切片振幅曲線不同去噪方法處理后的低頻噪聲，其中com-smo（粉紅色）代表鄰點(diǎn)加權(quán)平滑濾波法去除的低頻噪聲，com-dif（藍(lán)色）代表擴(kuò)散濾波方法（迭代次數(shù)為30次）去除的低頻噪聲，com-dif1（紅色）代表擴(kuò)散濾波方法（迭代次數(shù)為200次）去除的低頻噪聲，com-bp（深黃色）代表高通濾波法（斜坡為0.1～2 Hz）去除的低頻噪聲。

圖4 實(shí)際地震資料疊前逆時(shí)偏移低頻噪聲壓制前后水平切片中某inline方向曲線及其殘差對(duì)比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）可知：由于原始信號(hào)振幅較大，原始逆時(shí)偏移結(jié)果的波形曲線與經(jīng)不同去噪方法處理后的結(jié)果差異較小，但從去除的低頻噪聲曲線中可以明顯地看出不同去噪方法的壓制效果（圖3（b）），其中文中擴(kuò)散濾波方法（迭代次數(shù)為30次）去噪效果與鄰點(diǎn)加權(quán)平滑濾波法效果相當(dāng)，前者更加平滑，與文中擴(kuò)散濾波方法（迭代次數(shù)為200次）效果相比，迭代次數(shù)越大，地層細(xì)節(jié)特征將得到越多的保留，低頻波形曲線過(guò)度自然光滑，但如果該參數(shù)取得過(guò)大，則達(dá)不到低頻噪聲壓制的目的，此時(shí)的去噪結(jié)果將保留更多的低頻能量。因此，迭代次數(shù)的選擇需根據(jù)試驗(yàn)和處理要求來(lái)實(shí)現(xiàn)優(yōu)化選擇。而高通濾波法（斜坡為0.1～2 Hz）去除結(jié)果存在一定的波形抖動(dòng)（圖3（b）天藍(lán)色橢圓），這與頻域斜坡的選取存在一定的影響，若斜坡選取過(guò)陡，則吉普斯效應(yīng)就會(huì)出現(xiàn)在處理結(jié)果中。雖然高通濾波方法在橫向相對(duì)保幅性方面存在缺陷，但仍是目前低頻噪聲壓制中的一種高效的方法，因其采用逐道處理方式，需求內(nèi)存較小，不需要多次迭代計(jì)算。

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

3 結(jié)束語(yǔ)

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

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

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

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

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

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

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

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

[6]劉紅偉，李博，劉洪，等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn)[J].地球物理學(xué)報(bào)，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]陳可洋，吳清嶺，范興才，等.地震波疊前逆時(shí)偏移脈沖響應(yīng)研究與應(yīng)用[J].石油物探，2013，52（2）：163-170.

[10]康瑋，程玖兵.疊前逆時(shí)偏移假象去除方法[J].地球物理學(xué)進(jìn)展，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]劉紅偉，劉洪，鄒振，等.地震疊前逆時(shí)偏移中的去噪與存儲(chǔ)[J].地球物理學(xué)報(bào)，2010，63（9）：2171-2180.

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

[15]陳康，吳國(guó)忱.逆時(shí)偏移拉普拉斯算子濾波改進(jìn)算法[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.