陳可洋， 陳樹民， 李來林， 吳清嶺， 范興才， 劉振寬， 王建民

（1.中國石油大慶油田有限責(zé)任公司勘探開發(fā)研究院，黑龍江大慶 163712；2.中國石油大慶油田有限責(zé)任公司勘探事業(yè)部，黑龍江大慶 163453）

基于擴散濾波方法的地震數(shù)據(jù)修復(fù)技術(shù)

陳可洋1， 陳樹民1， 李來林1， 吳清嶺1， 范興才1， 劉振寬2， 王建民1

（1.中國石油大慶油田有限責(zé)任公司勘探開發(fā)研究院，黑龍江大慶 163712；2.中國石油大慶油田有限責(zé)任公司勘探事業(yè)部，黑龍江大慶 163453）

將基于偏微分方程的圖像修復(fù)技術(shù)引入到地震資料數(shù)字處理領(lǐng)域，提出基于擴散濾波方法的地震數(shù)據(jù)修復(fù)技術(shù).將待修復(fù)數(shù)據(jù)區(qū)域附近的未受損數(shù)據(jù)按照一定的擴散函數(shù)規(guī)則，通過擴散濾波方法將其擴散至待修復(fù)區(qū)域，每次迭代計算僅更新待修復(fù)區(qū)域內(nèi)的數(shù)據(jù)，未受損區(qū)域內(nèi)的數(shù)據(jù)保持不變，根據(jù)迭代前后的殘差與閾值比較，最終確定迭代終止的條件.地震數(shù)據(jù)插值處理以及地震數(shù)據(jù)的局部修復(fù)處理2個應(yīng)用實例表明，采用擴散濾波方法的地震數(shù)據(jù)修復(fù)技術(shù)可以修復(fù)地震數(shù)據(jù)，能夠有效恢復(fù)缺失的地震波場信息，在地震資料數(shù)字處理中有一定的應(yīng)用價值.

擴散濾波；數(shù)據(jù)修復(fù)；地震插值；地震處理；圖像修復(fù)

0 引言

圖像資料是人類獲取信息的一個重要來源，它可以用信號處理的方法對其進行描述.在數(shù)字圖像處理中，一項重要的研究內(nèi)容就是圖像的修復(fù)，它是對圖像上信息缺損區(qū)域進行圖像信息填充的過程，其主要思想就是用現(xiàn)有的圖像信息來自動恢復(fù)缺失的信息，并使觀察者無法察覺到圖像曾經(jīng)缺損或已被修復(fù).數(shù)字圖像的修復(fù)技術(shù)已廣泛應(yīng)用于諸多工程研究領(lǐng)域，如文物信息的修復(fù)、醫(yī)學(xué)圖像的修復(fù)、污跡的消除、舊照片中折痕的修復(fù)、圖像中文字的去除等等，甚至可以利用圖像修復(fù)技術(shù)達到一些特殊的藝術(shù)效果.因此，對數(shù)字圖像修復(fù)方法的研究是圖像處理領(lǐng)域的熱點之一，且具有重要的應(yīng)用價值［1－6］.目前，國內(nèi)外學(xué)者開展了諸多圖像修復(fù)技術(shù)方面的研究，并形成了較多成熟的圖像修復(fù)算法，主要包括基于偏微分方程和紋理合成.其中基于偏微分方程的算法是最早出現(xiàn)的圖像修復(fù)模型，其算法種類已較為豐富，其構(gòu)建思路是利用物理學(xué)中的熱傳導(dǎo)方程，將待修復(fù)區(qū)域周圍的已知信息沿著等照度線方向延伸到待修復(fù)區(qū)域中.其中具有代表性的圖像修復(fù)模型有：用三階偏微分方程來模擬平滑傳輸過程的BSCB模型［3］，用三階偏微分方程來模擬曲率驅(qū)動擴散（CDD）模型［5］等等.

在地震勘探中，由于野外采集常受到施工條件的綜合影響（如地形、障礙物等），在采集數(shù)據(jù)中常常在某些位置出現(xiàn)地震記錄的道間距過大、異常道、數(shù)據(jù)不完整或數(shù)據(jù)缺失等問題，這意味著一些重要地層的地震響應(yīng)信息的丟失，并可能在各種處理過程中形成許多不必要的噪聲.從目前的地震插值技術(shù)的應(yīng)用狀況來看，地震道插值方法可以在一定程度上恢復(fù)這些缺失的地震波場信息，可以提高偏移成像的精度，并切實有效地壓制噪音和抑制空間假頻出現(xiàn).因此，地震道插值是提高地震資料處理和成像精度的一種有效方法，也是地震數(shù)據(jù)處理中的一個難點與熱點問題［7－9］.

由于地震數(shù)據(jù)和圖像數(shù)據(jù)在數(shù)值計算方面的相似性，因此近幾十年來，圖像處理技術(shù)已逐漸引入到地震資料處理和解釋中，并取得了較好的應(yīng)用效果［10－14］，例如中值濾波技術(shù)、圖像邊緣增強技術(shù)、圖像細(xì)化算法、圖像噪聲壓制技術(shù)等等，這些技術(shù)的廣泛應(yīng)用對于實際地震資料的預(yù)處理和修飾性處理等方面發(fā)揮了非常重要的作用.在前人研究的基礎(chǔ)上，我們將圖像修復(fù)技術(shù)引入到地震資料數(shù)字處理領(lǐng)域，以基于偏微分方程的擴散濾波方法為例，詳細(xì)探討實際地震數(shù)據(jù)的修復(fù)方法，并在地震數(shù)據(jù)的插值處理和地震數(shù)據(jù)的局部修復(fù)等方面進行了實例應(yīng)用，該技術(shù)的引入以期為實際地震資料處理工作提供方法指導(dǎo).

1 基本理論

Perona和Malik［15］首次提出了用于擴散濾波處理的偏微分方程

式（1）中，t為擴散時間（可不考慮時間步長dt，常用迭代次數(shù)N代替），div為散度算子，?是梯度算子，U為t時刻的擴散濾波結(jié)果，其中U0為t＝0時刻的原始地震數(shù)據(jù)，即擴散濾波迭代計算的初始條件，g（·）是擴散函數(shù)，為一個有界非負(fù)的遞減函數(shù)，Perona和Malik［15］建議其形式為

與數(shù)字圖像修復(fù)技術(shù)相類似，地震數(shù)據(jù)修復(fù)技術(shù)的主要思想可以表述為：將待修復(fù)區(qū)域附近的信息沿著等照明度方向擴散到待修復(fù)區(qū)域，并按照一定的規(guī)則控制擴散方向，從而達到修復(fù)地震數(shù)據(jù)的目的，即用偏微分方程中已知的初始條件和邊界條件來求得未知的待修復(fù)區(qū)域中的地震信息.分析地震數(shù)據(jù)可知，與地震同相軸相垂直的方向是地震數(shù)值變化最快的方向，而順著地震同相軸的方向是地震數(shù)值變化最慢的方向，即圖像修復(fù)技術(shù)中的等照明度方向.根據(jù)待修復(fù)區(qū)域周圍未受損區(qū)域中的地震數(shù)據(jù)信息來修正待修復(fù)區(qū)域中的每一個地震采樣點數(shù)值，且修復(fù)過程只改變待修復(fù)區(qū)域內(nèi)的地震數(shù)值，未受損區(qū)域內(nèi)的地震數(shù)值保持不變.基于擴散濾波方法的地震數(shù)據(jù)修復(fù)技術(shù)與熱傳導(dǎo)方程中能量的傳遞過程類似，地震數(shù)據(jù)的信息是在一定的約束條件下以振幅能量的形式向待修復(fù)區(qū)域進行擴散，最終使整個地震數(shù)據(jù)達到穩(wěn)定狀態(tài).這種約束條件就是根據(jù)地震數(shù)據(jù)本身的結(jié)構(gòu)特性（即擴散函數(shù)）來控制各個空間方向的擴散強度.然后根據(jù)最新迭代得到后的地震數(shù)據(jù)和前一次迭代得到的地震數(shù)據(jù)之間的差值與給定的閾值進行比較，如果該差值小于給定的閾值，則可以結(jié)束迭代，并得到修復(fù)處理后的地震數(shù)據(jù).若該差值大于給定的閾值，則繼續(xù)迭代計算，直到小于給定的閾值后結(jié)束.

借鑒圖像修復(fù)模型及其基本原理［1］，我們構(gòu)建如下針對地震數(shù)據(jù)的修復(fù)模型為

式（2）中，t代表迭代的時間，Δt代表迭代的時間步長，Utn代表對t時刻地震數(shù)據(jù)Ut的修正量，可表示為Utn＝δLt·Nt，代表Lt沿著擴散濾波Nt方向的變化量，其中δLt＝div（g（‖?U‖）·?U）．當(dāng)該修正量Utn變化不明顯時，δLt·Nt≈0此時Ut＋Δt≈Ut，則可以結(jié)束迭代計算，從而得到地震數(shù)據(jù)修復(fù)后的處理結(jié)果.

具體的步驟

1）確定待插值樣點在數(shù)據(jù)體中的空間位置，并建立待插值樣點的數(shù)據(jù)索引；

2）將原始數(shù)據(jù)作為擴散濾波方程的邊界條件，采用式（1）進行數(shù)據(jù)的迭代計算；

3）根據(jù)地震數(shù)據(jù)的修復(fù)模型（式（2））對待修復(fù)數(shù)據(jù)區(qū)域進行局部迭代更新；

4）待修復(fù)區(qū)迭代前后的殘差絕對值為0或小于門檻值作為迭代收斂和終止的條件，若大于門檻值，則繼續(xù)2）至4）步驟，直至達到迭代收斂條件；

5）對原始數(shù)據(jù)不做任何修改，最終輸出修復(fù)后的地震數(shù)據(jù).

2 應(yīng)用實例

2.1 地震數(shù)據(jù)的插值處理

地震數(shù)據(jù)處理中常見的就是由于采集因素導(dǎo)致的道間距過大的問題，若對這種資料直接進行后續(xù)的地震常規(guī)處理和成像處理均將引入空間假頻問題，直接影響后續(xù)的精細(xì)地震解釋的質(zhì)量，這增加了實際地震勘探的風(fēng)險.這里我們以兩種最常見的地震數(shù)據(jù)類型為例開展基于擴散濾波方法地震數(shù)據(jù)修復(fù)技術(shù)的應(yīng)用研究，主要針對疊前和疊后地震數(shù)據(jù)的插值處理.

圖1（a1）和圖1（a2）分別為某地區(qū)地震疊前時間偏移剖面和理論合成的單炮模擬記錄，圖1（b1）和圖1（b1）分別為圖1（a1）和圖1（a2）經(jīng)地震數(shù)據(jù)抽稀后的結(jié)果，圖1（c1）和圖1（c2）分別為圖1（b1）和圖1（b1）經(jīng)擴散濾波方法地震數(shù)據(jù)修復(fù)后的結(jié)果.分析圖1（a）和圖1（c）可知，原始數(shù)據(jù)和經(jīng)擴散濾波方法修復(fù)后的地震同相軸連續(xù)性等特征基本保持一致，同時，原始疊后剖面（圖2（a））和炮記錄（圖2（b））中被抽掉的地震道波形曲線與文中方法插值恢復(fù)后的地震波波形曲線吻合度較高.由此可見，基于擴散濾波方法的地震修復(fù)技術(shù)在疊前和疊后地震數(shù)據(jù)插值處理中取得了較好的應(yīng)用效果.

圖1 擴散濾波方法在疊前和疊后地震數(shù)據(jù)插值中的應(yīng)用Fig.1 Application of diffusion filtermethod in pre-stack and post-stack seismic data interpolation

2.2 地震數(shù)據(jù)的局部修復(fù)處理

當(dāng)磁盤上存儲的地震數(shù)據(jù)受到外界復(fù)雜環(huán)境因素影響的情況下，該數(shù)據(jù)中將引入一些異常數(shù)值，直接影響后續(xù)的利用.與正常數(shù)值相比，若這些異常值特征差異較為明顯，則可根據(jù)一些算法將這些數(shù)值篩選出來并做好待修復(fù)的標(biāo)記，再采用地震數(shù)據(jù)修復(fù)技術(shù)來實現(xiàn)待修復(fù)數(shù)據(jù)區(qū)域的修復(fù)處理.為此，我們以理論合成的單炮數(shù)據(jù)為例，模擬被破壞的待修復(fù)地震數(shù)據(jù).

圖3（a）為理論合成的單炮記錄，圖3（b）為基于能流密度矢量篩選出來的左行波數(shù)據(jù)來描述待修復(fù)數(shù)據(jù)，圖3（c）基于能流密度矢量篩選出來的右行波數(shù)據(jù)來描述待修復(fù)數(shù)據(jù)，圖3（d）對圖3（b）的地震數(shù)據(jù)修復(fù)結(jié)果，圖3（e）為對圖3（c）的地震數(shù)據(jù)修復(fù)結(jié)果.分析圖3可知，待修復(fù)的數(shù)據(jù)由于采用篩選的方法，數(shù)據(jù)中存在大量的空值區(qū)域（圖3（b）和圖3（c）），數(shù)據(jù)缺失信息較多，且同相軸的連續(xù)性較差.而采用基于擴散濾波方法的地震數(shù)據(jù)修復(fù)技術(shù)后（圖3（d）和圖3（e）），同相軸的連續(xù)性增強，空值區(qū)域內(nèi)的地震數(shù)據(jù)也得到了有效的恢復(fù).由此可見，地震數(shù)據(jù)修復(fù)技術(shù)在實際地震數(shù)據(jù)修復(fù)中可以取得較好的應(yīng)用效果.

圖2 原始和插值恢復(fù)后的地震道波形曲線Fig.2 Seismic traces wave forms of original trace and interpolated recovered trace

圖3 擴散濾波方法在地震數(shù)據(jù)修復(fù)處理中的應(yīng)用Fig.3 Application of diffusion filtermethod in seismic data repair processing

3 結(jié)束語

將圖像修復(fù)技術(shù)引入地震資料的數(shù)字處理領(lǐng)域，構(gòu)建了地震數(shù)據(jù)的修復(fù)技術(shù)，并給出具體的步驟.該技術(shù)采用擴散濾波方法，將擴散函數(shù)作為每次迭代修復(fù)的約束條件，從而確定擴散的方向和強度，同時將未受損區(qū)域內(nèi)的地震數(shù)據(jù)擴散至待修復(fù)區(qū)域內(nèi)，并使未受損區(qū)域內(nèi)的地震數(shù)據(jù)數(shù)值始終保持不變，僅迭代更新待修復(fù)區(qū)域內(nèi)的數(shù)據(jù).將迭代前后的殘余與給定閾值進行比較，可最終確定迭代終止的條件.通過地震數(shù)據(jù)的插值處理和地震數(shù)據(jù)的局部修復(fù)處理應(yīng)用表明，本文提出的地震數(shù)據(jù)修復(fù)技術(shù)可以較好地實現(xiàn)地震數(shù)據(jù)缺失信息的修復(fù)，有效恢復(fù)了部分缺失的地震信號，在實際地震資料處理中具有一定的應(yīng)用潛力.該地震數(shù)據(jù)修復(fù)技術(shù)避免了常規(guī)地震插值方法在F-K域或T-K域進行變換處理并受插值角度的限制等缺點，該方法在時空域來進行處理，因此具有計算過程簡單易實現(xiàn)的優(yōu)點，適合于實際地震資料處理中異常道（或野值道）被剔除后波場信息的有效恢復(fù)和三維地震數(shù)據(jù)體加密插值處理等.若待修復(fù)樣點過密且已知地震信息較少時，文中方法計算結(jié)果的誤差將較大.

［1］劉庚農(nóng).圖像修復(fù)技術(shù)研究［D］.合肥：合肥理工大學(xué)，碩士學(xué)位論文，2011.

［2］王麗坤.三維各向異性擴散濾波方法研究與應(yīng)用 ［D］.成都：成都理工大學(xué)，碩士學(xué)位論文，2011.

［3］Bertalmio M，Sapiro G，Caselles V，et al.Image inpainting［C］∥Proceedings of International Conference on Computer Graphics and Interactive Techniques，New Orleans，Louisiana，USA，2000：417－424.

［4］Chan T F，Shen JH.Mathematicalmodels for local non-texture inpainting［J］.SIAMJournal of Applied Mathematics，2002，62（3）：1019－1043.

［5］Chan T F，Shen JH.Non-texture inpainting by curvature-driven diffusions（CDD）［J］.Journal of Visual Communication and Image Representation，2001，12（4）：436－449.

［6］Chan T F，Kang SH，Shen JH.Euler's elastica and curvature based inpainting［J］.SIAMJournal of Applied Mathematics，2002，63（2）：564－592.

［7］Zhang WY.Amplitude-preserved prestack interpolation of seismic data with 3-D Fourier transform ［D］.Beijing：China University of Petroleum.Engineering Master's Thesis，2011.

［8］陳可洋，范興才，吳清嶺，等.提高逆時偏移成像精度的疊前插值處理研究與應(yīng)用 ［J］.石油物探，2013，52（4）：409－416.

［9］陳可洋，吳清嶺，李來林，等.松遼盆地三維地震資料連片處理關(guān)鍵技術(shù)及其應(yīng)用效果分析 ［J］.巖性油氣藏，2012，24（2）：87－91.

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

［11］陳可洋，楊微，吳清嶺，等.地震反射波與散射波波場分離方法初探［J］.巖性油氣藏，2013，25（2）：76－81.

［12］Lu Y H，Lu WK.Edge-preserving polynomial fittingmethod to suppress random seismic noise［J］.Geophysics，2009，74（4）：69－73.

［13］Fehmers G C，Hocker C FW.Fast structural interpretation with structure-oriented filtering［J］.Geophysics，2003，68（4）：1286－1293.

［14］Lavialle O，Pop S，Germain C，et al.Seismic fault preserving diffusion［J］.Journal of Applied Geophysics，2007，61（2）：132－141.

［15］Perona P，Malik J.Scale-space and edge detection using anisotropic diffusion［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12（7）：629－639.

Seism ic Data Repair Technology Based on Diffusion Filtering Method

CHEN Keyang1，CHEN Shumin1，LILailin1，WU Qingling1，F(xiàn)AN Xingcai1，LIU Zhenkuan2，WANG Jianmin1
（1．Exploration and Development Research Institute ofDaqing Oilfield Company Limited ofCNPC，Daqing，Heilongjiang 163712，China；2．Exploration Company ofDaqing Oil Field Co．Ltd ofCNPC，Daqing，Heilongjiang 163453，China）

We introduce figure repair technology based on partial differential equation into seismic data digital processing，and put forward seismic data repair technology based on diffusion filteringmethod.The technology diffuses undamaged area's data to data to be repaired using diffusion filteringmethod with a certain diffusion function.In each iteration，only the data to be repaired is updated，and undamaged area's data is remained unchanged.Final iteration termination condition can be ensured by comparing difference with threshold before and after iterations.Two application examples of seismic data interpolation processing and seismic data local repair processing show that the technology can achieve the purpose of repairing seismic data.It recovers effectively lost seismic wave field information.Therefore，it can be used in practical seismic data digital processing.

diffusion filter；data repair；seismic interpolation；seismic processing；figure repair

date：2013－08－17；Revised date：2013－11－26

P315.3

A

1001-246X（2014）04-0465-06

2013－08－17；

2013－11－26

國家重點基礎(chǔ)研究發(fā)展計劃（973）（2009CB219307）資助項目

陳可洋（1983－），男，碩士，工程師，主要從事高精度地震波傳播模擬與逆時成像、多線程并行計算與模塊開發(fā)、實際地震資料數(shù)字處理方法研究，E-mail：keyangchen＠163.com