基于統計分布的Shearlet域連片數據一致性處理方法

田 坤 王德營 劉立彬 王延光 王常波 張學濤

(①中國石化勝利油田分公司物探研究院, 山東東營 257022；②山東科技大學地球科學與工程學院, 山東青島 266590)

0 引言

對進入勘探開發中、后期的老油田，地震資料處理目標逐漸由“單塊三維”轉向“連片三維”，以實現對區域整體構造及各構造帶之間轉換關系的研究[1]。受地表條件等因素的影響，連片探區的不同區塊往往采用多種震源和不同檢波器進行觀測接收，相同區塊不同期次采集的資料品質參差不齊[2]，導致相鄰區塊地震資料的振幅、頻率、波形和相位等存在較大差異[3]，嚴重影響了數據的連續對比分析和地質儲層研究工作[4]。因而，連片資料的一致性校正處理尤為重要。

目前，解決不同區塊連片資料拼接一致性問題的主要方法有反褶積參數調整法、匹配濾波法以及多種技術組合應用法。

反褶積參數調整法主要解決不同區塊數據分辨率差異的問題，王西文等[5]通過在反褶積處理環節選擇不同的預測步長，盡可能地實現有效頻譜寬度、主頻和優勢頻寬能量等的一致，這種方法在實際應用中具有一定的效果，但反褶積參數實驗相對繁瑣，對于十幾或幾十個區塊進行資料拼接的大連片處理，選取合理的參數十分耗時、費力。

匹配濾波法可分為時變、時不變兩類。時不變的匹配濾波法有子波處理法和直接匹配濾波法，其中子波處理法也稱子波整形或子波匹配，子波處理法假設連片探區相同位置處反射系數相同，認為地震記錄的差異由不同子波特征導致。王元波[6]在兩個區塊重疊區各選取多個面元的地震記錄，分別進行混道和去噪處理后，求出相對穩定的子波，計算子波整形因子并應用其消除子波差異，該方法處理質量受拼接帶資料質量的影響嚴重；宋玉龍[7]利用直接匹配濾波法對兩個區塊重疊區數據進行處理，采用最小二乘法直接計算匹配濾波算子，消除地震剖面間的振幅差異。時變的匹配濾波法包括基于小波變換的子波處理法和基于小波變換的匹配濾波法。王西文等[8]提出基于小波變換的子波處理法將重疊區的兩期數據利用不同尺度的小波函數進行分頻展開，沿時間方向開小時窗計算兩者的能量標定量，并插值得到最佳重構系數，計算各尺度對應的重構系數后進行子波重構，實現一致性處理；程金星等[9]基于小波變換的匹配濾波法將重疊區的兩期數據進行小波變換，在不同尺度上利用最小平方法分別計算兩期數據的最佳匹配濾波因子，并應用其實現一致性處理的目的；為了實現振幅、時差、頻率、相位和波形的全面校正，云美厚等[10]提出了互均衡處理技術，采用最小平方算法從疊后記錄中求取多個匹配濾波算子進行一致性處理。

多種技術組合應用法利用多種振幅、頻率、相位等一致性處理方法的組合應用解決連片數據的一致性問題。Mohan等[11]采用增益恢復、地表一致性反褶積和地表一致性均衡處理解決連片數據子波不一致問題；戴軍文等[12]利用地表一致性振幅補償、子波整形和振幅均衡處理消除連片處理中的能量不均帶來的偏移畫弧問題；裴江云等[13]針對疊前連片處理中不同區塊的非一致性問題，聯合應用幾何擴散補償、地表一致性振幅補償、地表一致性反褶積和子波整形處理等技術以實現數據的一致性處理；Greer等[14]提出了一套連片一致性處理的流程，首先通過非平穩平滑濾波處理對高分辨率結果實現振幅和頻率的均衡，隨后利用局部相似掃描方法估計和消除兩期數據剖面中的局部時移量，最后利用最小平方反演方法融合兩期數據；蔣波[15]應用非剛性匹配、匹配濾波和譜整形等方法改善連片數據處理的子波一致性；曾華會等[16]綜合利用疊前平均振幅譜匹配、混合震源相位匹配及疊前剩余時差校正等處理技術，消除了不同激發方式造成的地震資料品質差異。

匹配濾波法的一致性校正因子是通過兩塊資料的重疊區域數據計算得到，該方法適用的前提是必須存在一定范圍的資料重疊區，其處理質量嚴重依賴于重疊區的資料品質，并且一般只適用于疊后數據的拼接處理。反褶積參數調整法和多種技術組合應用法可以處理疊前數據，但參數調整繁瑣，且大都不具有時變處理能力。針對連片數據一致性處理方法存在的問題，本文提出一種基于Shearlet變換的連片數據一致性處理方法，可以處理疊前數據。該方法利用了Shearlet變換的多尺度和多方向性，在優選目標數據的基礎上，利用Alpha-trim濾波方法估計待處理數據和目標數據在各尺度和方向上的時變均值和均方差，提取并消除兩者的趨勢差異，從而實現連片數據的一致性處理，并通過模型數據和實際資料對該方法進行了驗證。

1 方法原理

1.1 Shearlet變換的基本原理

Guo等[17-18]基于合成小波理論，利用仿射系統把幾何分析和多尺度分析結合起來構造了Shearlet變換，其后又將該變換推廣到三維空間。Shearlet變換的構造方法如下。

(2)定義Shearlet波原子ψj,s,m為

(1)

式中：MAS=SsAa；m代表位置平移量;x為位置變量。令ψ∈L2(R2)，其中，L2(R2)表示二維實數域上平方可積函數空間，并滿足下列條件：

函數f的Shearlet變換可定義為

Sf(j,s,m)=〈f,ψj,s,m〉

(2)

式中〈·，·〉表示內積。

圖1是取不同a、s值所對應的水平方向和垂直方向的Shearlet波原子在頻域的支撐[19]，其分割方式與Curvelet變換相似。

圖1 不同a、s值對應的水平方向(左)和垂直方向(右)的Shearlet波原子在頻域的支撐

1.2 一致性校正因子估算

對連片探區中各區塊數據分別進行地表一致性振幅補償和地表一致性反褶積處理，受資料品質、地層吸收差異等因素的影響，連片數據的時間、頻率在空間上仍存在一定的差異。

選取一個處理的目標模型消除這種差異，凌云等[20]利用時頻空間域球面發散與吸收補償方法選取模型炮作為目標模型；熊艷梅等[21]優選部分炮檢距數據作為目標模型，進行非剛性匹配處理實現數據的一致性處理。本文借鑒后者做法，在連片探區數據品質相對較好的區塊內選取一定數量品質較好的數據作為目標模型，連片探區中的其他數據都以該目標數據為校正標準，估計Shearlet域的一致性校正因子，消除目標模型數據與待處理數據的統計差異量，實現連片探區的一致性處理。

一致性校正因子的估計包括兩個方面，一個是反映目標模型數據時、頻、空間變化的統計趨勢估計，另一個是待處理數據的時、頻、空間變化趨勢的提取。

1.2.1 目標模型數據的時空變趨勢估計

在工區A中選取N炮質量相對較好的炮集數據作為模型數據，記為Xi，其中i=1,2，…,N為模型炮的炮號。對Xi進行Shearlet變換，為便于表述，此處僅給出二維的情況，式(2)中位置參數m用二維變量mk1和mk2表示。變換后的Shearlet系數計算公式為

SX(i,j,s,mk1,mk2)=〈xi,ψj,s,mk1,mk2〉

(3)

Alpha-trim均值濾波又稱非線性平滑濾波[22]，可以消除絕大部分極值因素的影響，同時保留數據非線性變化趨勢。利用Alpha-trim均值濾波確定Shearlet系數SX(i,j,s,mk1,mk2)的時變均值，計算方法如下：

(1)由小到大排序系數；

(2)拋棄總數中α比例的極大值和極小值部分，其中0≤α≤0.5；

(3)計算剩余部分的均值。

當α=0時該方法退化為均值濾波；當α=0.5時該方法退化為中值濾波。

為了提高Shearlet系數時空變趨勢估計結果的穩定性，對多個模型炮集數據的估計結果進行統計平均處理，得到尺度j、方向s下的Shearlet系數多炮統計平均時空變均值和時空變均方差

(4)

(5)

針對不同j和s循環上述步驟，直至全部遍歷，完成模型炮集數據的時空變趨勢估計。

1.2.2 待處理數據的時空變趨勢估計

由于Shearlet變換是基于頻域中的二進尺度分割進行，因而在相同頻帶寬度下，高頻端分解的尺度個數比低頻端少，頻帶分割不精細，嚴重影響了頻率一致性校正的精度。針對這一問題，在頻域中分別計算模型數據和待處理數據的統計頻譜，以模型數據的統計頻譜為期望，計算頻率一致性校正算子，對待處理數據進行統計頻率一致性校正處理，處理結果記為Y。

(6)

式中：c、d表征二維高斯函數兩個維度；λ為均方差。λ值越大，高斯窗函數越平坦，濾波效應越明顯；λ值越小，高斯窗函數越陡峭，濾波效應降低。

(7)

(8)

經過二維Alpha-trim滑動濾波處理和二維高斯濾波處理后得到的時空變均值和均方差結果，基本消除了局部擾動異常，并且能夠穩定地反映Shearlet系數隨時間和空間的變化趨勢。

1.2.3 一致性校正因子

按照上述方法計算得到模型數據和待處理數據相對穩定的時空變均值和均方差后，通過下式計算一致性校正因子

(9)

(10)

1.3 Shearlet域一致性校正處理

計算出Shearlet域一致性校正因子后，按照下式調整待處理數據Shearlet系數的時變均值

(11)

(12)

在實際數據處理時，應注意通過定義時窗避開直達波和淺層折射波對一致性校正因子估計質量的影響。

2 模型試算與實際資料處理

為了驗證本文方法的可靠性和有效性，分別對模型數據和實際數據進行了一致性處理和分析。

2.1 模型試算

設計一個速度模型，在模型地表的不同位置利用主頻為35Hz、極大值為2.5的Riker子波和主頻為30Hz、極大值為1的Riker子波模擬激發，分別正演得到炮集數據A和數據B(圖2)。由圖可見，數據A(圖2a)的能量比數據B(圖2b)強；由頻譜統計圖可以看出兩者存在頻率差異，數據A的主頻略高并且頻帶略寬(圖2c)；由數據A、數據B的振幅統計直方圖可以看出(圖2d)，二者的振幅幅值在零值附近分布居多，但數據A的幅值分布相對分散，較小值明顯多于數據B，其統計直方圖較“胖”，表明數據A的均方差較大。

圖2 正演模型數據及其統計頻譜、振幅分析

分別計算數據A和數據B中所有采樣點值的均值E

(13)

式中：M表示炮集數據中所有采樣點的個數；Al表示第l個采樣點的幅值。

再分別計算數據A、數據B的均方差σ

(14)

由式(13)計算出數據A、數據B的均值分別為EA=5.06×10-5和EB=1.20×10-4，兩者均值差異很小，接近零值(均值反映了數據中直流分量的強弱)。利用式(14)計算出數據A、數據B的均方差分別為σA=2.7213和σB=0.6275，說明炮集數據的能量越強，其對應的均方差越大。下面以數據A作為目標模型，對數據B進行一致性處理。

利用圖2c中數據A、B的統計頻譜，設計頻率一致性校正因子，對數據B進行頻域統計一致性校正處理。圖3為處理結果，對比分析可以看出，處理前(圖2b)、后(圖3a)炮集數據的能量變化不明顯。由圖3b可知，經過頻率統計一致性處理后，數據B的統計頻譜與目標模型數據A的統計頻譜基本相吻合。對比處理前(圖2d)、后(圖3c)炮集數據的振幅統計直方圖可知，處理后數據B中接近零值的數據個數急劇下降，朝著遠離零值的趨勢變化，但整體來看，數據B的振幅統計直方圖依然較“瘦”，數據的均方差略有增加，但與數據A的振幅統計直方圖相比，依然存在較大的差異。

對數據A(圖2a)和頻率一致性處理后的數據B(圖3a)進行Shearlet變換。圖4顯示f-k域的Shearlet尺度分割及其變換系數的對比。圖4a為圖2a的f-k譜，圖4b為某個尺度和方向上的Shearlet的f-k域分割窗，利用該分割窗分別對數據A和頻域一致性處理后的數據B進行分割，得到對應的Shearlet系數(圖4c、圖4d)。可以看出，由于模型數據和待處理數據的反射波能量、傾角和時空域分布都存在較大差異，經過Shearlet變換后，在相同尺度和方向上的系數也存在較大差異。

圖4c、圖4d中Shearlet系數的幅度反映了局部數據在該尺度和方向的能量強弱，Shearlet局部系數依然是波形的形態特征。對于這種形態而言，均方差能夠反映數據在特定尺度和方向的能量強弱，一般來說，數據能量越強，其Shearlet系數的均方差越大。

利用本文方法對圖4c、圖4d的Shearlet系數進行一致性校正因子估算。圖5a給出了二維滑動窗內數據由小到大排列及Alpha-trim濾波參數α=0.2時剔除了極值后的數據對比。由圖可見，數據較大、較小值的數量均少但數值幅度較大，這會嚴重影響均值計算結果的質量，尤其是當數據中存在局部異常噪聲的情況下。經Alpha-trim均值濾波處理剔除一定比例的較大、較小值后，能夠消除局部極值對時變均值估計的影響，提高估計的穩定性。圖5b為圖4d的Shearlet系數經過一致性校正處理后的結果。綜合對比可見，處理后數據B系數的幾何特征(傾角、各反射間幾何關系等)都與處理前保持一致，但能量分布的相對關系有較大改變，整體上與目標模型數據的Shearlet系數(圖4c)一致。

圖3 頻域統計一致性校正處理結果及其統計分析

圖4 Shearlet變換的FK域分割及其變換系數對比

圖5 Shearlet系數一致性校正處理

在Shearlet域完成各個尺度和方向的一致性校正處理，將處理后的系數反變換到時空域得到最終結果(圖6)。對比圖2a、圖3b和圖6a可以看出，經過Shearlet域一致性處理后，數據A和數據B的振幅強度在時空域的變化特征基本一致，并且反射、繞射波的特征也與處理前的數據特征保持一致，表明處理前、后炮集數據的幾何特征和波動特征都得到較好的保持。進一步分析統計頻譜(圖6b)可以看出，經過本文方法處理后，數據B的統計頻譜與模型數據A基本一致，在80～120Hz頻段內，處理后數據B的統計頻譜能量有一定的提升，略高于數據A，側面反映了隨著頻率的增加Shearlet變換對頻率一致性的處理能力有所減弱。處理后數據B和數據A的振幅統計直方圖(圖6c)顯示，兩者分布形態基本相似，除零值附近的分布特征略有不同外，其他幅度上的分布基本一致，表明處理后數據和模型數據在振幅統計上具有一致性。

圖6 Shearlet域一致性校正結果及其統計分析

2.2 實際資料處理

對實際連片區塊C、區塊D的炮集數據進行一致性校正處理，進一步驗證本文方法的有效性。區塊C(圖7a上)和區塊D(圖7下)的數據都進行了地表一致性振幅補償和地表一致性反褶積處理。對比可見，區塊C數據的信噪比高于區塊D數據。圖7c的統計頻譜顯示兩塊數據的頻帶寬度大致相當，區塊D數據頻帶略寬，并且40～60Hz范圍內的能量明顯強于區塊C。圖7d振幅統計直方圖顯示兩者的振幅分布存在一定的差異，區塊D數據的均方差略大，區塊C數據振幅落在-1000～1000范圍內的個數明顯多于區塊D。

圖8為區塊D數據進行頻域一致性校正處理的結果，對比圖8a和圖7b炮集數據可以看出，兩者無明顯差異。圖8b統計頻譜對比顯示，經過頻域一致性校正處理后，區塊D數據的統計頻譜與區塊C數據的統計頻譜基本一致。圖8d中振幅統計直方圖對比可以看出，頻域一致性校正處理后區塊D數據振幅統計直方圖的均方差略有降低，與區塊C數據的振幅統計直方圖吻合度有所提高。

圖7 實際觀測數據及其統計分析

圖8 頻域統計一致性校正處理

圖9為區塊D數據采用本文方法一致性處理的最終結果，對比圖9a和圖7b的炮集數據可以看出，一致性處理后剖面的能量分布更加均衡，局部噪聲也有所壓制。對比處理后的統計頻譜(圖9b)可以看出，本文方法能夠較好地校正兩區塊數據間的統計頻帶差異，處理后兩者的主頻、頻帶寬度等都基本達到一致。對比處理后兩區塊數據的振幅統計直方圖(圖9c)可見，本文方法極大地消除了兩區塊數據間的振幅統計分布差異，實現了兩塊數據振幅、頻率的時空一致性，從而驗證了本文方法的有效性。

圖9 Shearlet域一致性校正結果及其統計分析

圖10進一步對比了區塊D數據處理前、后地震道的波形變化，從圖中可以看出，處理前區塊C和區塊D的地震記錄的波形存在較大差異，區塊D數據的振幅明顯強于區塊C數據。經過本文方法一致性處理后，區塊D數據的振幅與區塊C數據的振幅范圍基本一致，進一步驗證了本文方法一致性處理效果。

圖10 區塊C和處理前、后區塊D數據地震道波形對比

需指出的是，受采集因素的影響，不同批次可控震源與炸藥震源采集的數據往往存在較大的相位差異，在實際連片數據處理時，應進一步消除不同數據之間相位差異。上述方法僅消除了振幅、頻率的時空上的差異，還需要采用小相位轉換或相位旋轉等處理手段消除不同區塊間的相位差，最終實現連片數據振幅、頻率、波形和相位的一致性處理。

3 結束語

基于Shearlet變換的連片數據一致性處理方法是在Shearlet域通過二維Alpha-trim均值滑動濾波處理提取時變的均值，并進一步計算均方差，通過調整待處理數據和目標數據間的統計分布達到連片數據的振幅、頻率在時空域相對一致的目的。不同頻帶的 Shearlet變換尺度分割精細程度存在差異，高頻段一致性校正能力較弱，為了解決這一問題，在Shearlet域一致性處理前，先在頻域進行統計頻譜的一致性處理。用本文方法進行模型試算和實際資料處理及對比分析，結果證實了方法的可靠性和有效性，具有一定的實際應用價值。