全局約束反演多道吸收補償方法

劉金濤 王 孝 王小衛 蘇 勤 劉夢麗 袁 煥

(中國石油勘探開發研究院西北分院，甘肅蘭州 730030)

0 引言

眾所周知，地震波在地下介質中傳播時，由于大地濾波作用，會導致地震波能量衰減、速度頻散，降低分辨薄層的能力。目前已提出了許多地層吸收衰減模型(Q模型)描述這一現象，其中常Q模型是地震勘探領域使用最多的地層吸收衰減模型[1]。地震波的衰減效應破壞了地震子波在時間方向的穩態特征，導致傳統的反褶積方法很難從非穩態地震記錄中準確地恢復反射系數序列。因此，基于吸收補償的提高分辨率方法具有明確的物理意義和較大應用潛力[2-8]。

反Q濾波是吸收補償的主要方法，有多種實現方式，其中，波場延拓和稀疏反演是兩類應用最為廣泛的方法。由于波場延拓方法在外推過程中存在放大函數，極易放大高頻噪聲，結果不穩定。因此，稀疏反演類反Q濾波方法應運而生[9-10]。Zhang等[11]將衰減補償看作基于反演的最小二乘問題，并用貝葉斯理論進行正則化約束；Braga等[12]在小波域實現了L2范數約束的反Q濾波；Oliveira等[13]在頻率域實現了L1范數約束的衰減補償，Wang等[14]對其進行改進，在時間域實現了L1范數約束的衰減補償。以上方法都只考慮了地震道縱向上的信息，而沒有考慮相鄰道之間的相關性，即都是單道補償方法，因而結果往往會出現橫向不連續問題。為了解決這一問題，將相鄰道之間的相關性作為先驗信息加入反演系統進行多道反演補償的方法逐漸發展起來。

Zhang等[15]提出了一種多道反演方法，該算法將橫向約束融入反演算法中，反演的結果較常規單道算法具有更好的橫向連續性。Kazemi等[16]假設反射系數序列稀疏分布，提出了一種非線性多道盲源反褶積方法，主要優點是抗噪性較強且不需要子波的先驗信息。Yuan等[17]先將數據轉化到變換域，然后在變換域數據上加一個約束項，實現了多道波阻抗反演，在一定程度上緩解了單道反演算法的橫向不連續性問題；在此基礎上，Yuan等[18]提出了一種在頻率—波數域的多道衰減補償方法，增強了地震同相軸橫向連續性，提高了吸收補償的穩定性和抗噪性。

f-x域預測濾波是描述線性事件空間預測特性的有用工具，主要被用于數據插值和去噪領域[19-22]。為提高地層吸收補償的穩定性和抗噪性，本文提出了一種基于反演的全局約束吸收補償方法。該方法首先在f-x域計算信號預測算子[23-25]，用于表征地震信號在空間上的結構關系和能量變化；然后，將該算子作為橫向正則化約束引入到縱向吸收補償系統，建立具有雙向(縱向和橫向)即全局約束的反演目標函數；最后應用模型數據和實際數據驗證了方法的有效性。

1 地層吸收模型

1.1 Q的定義

當平面波通過均勻的黏彈性介質時，其振幅衰減和速度頻散可以用一個無量綱的參數Q描述

(1)

式中：ω為圓頻率；α(ω)和v(ω)分別為衰減系數和相速度。由于品質因子Q必須為正數，所以

(2)

進一步假設

(3)

則式(1)可以近似為

(4)

1.2 修正的Kolsky-Futterman模型

Kolsky[26]假設在測量頻帶范圍內衰減系數α(ω)與頻率呈嚴格的線性關系，即

(5)

式中ω0為參考圓頻率。同時，Kolsky定義其相速度為

(6)

將式(5)和式(6)代入式(4)，可得

(7)

對于較大的品質因子，即Q(ω0)?1，有

(8)

式中γ=1/[πQ(ω0)]。則相速度和品質因子可近似為

(9)

(10)

其中自由選擇的參考頻率ω0小于最低的測量頻率。

Wang等[27]指出，式(6)相速度僅僅是對于ω?ω0的一個漸近表達式。由于地震勘探資料的頻率范圍相對較低(一般小于500Hz)，其修正的相速度公式為

(11)

式中：h是一個頻率不相關的常數；ωh為重新定義的調諧參數。

對于修正的相速度模型，其相應的修正Q模型為

(12)

本文所用的地層吸收模型即為上述模型。

2 基于反演的縱向約束吸收補償方法

在完全彈性介質中，可以用褶積模型簡單描述地震數據采集過程，即地震記錄可由地震子波和反射系數序列褶積外加噪聲得到

d(t)=r(t)*w(t)+n(t)

(13)

式中：d(t)為地震記錄；r(t)為反射系數；w(t)為震源子波；n(t)為噪聲。

在黏彈性介質中，由于地層的吸收衰減作用，導致震源子波會隨著傳播時間而變化。因此，需要將穩態的褶積模型推廣到非穩態

(14)

式中w(τ，t)為時變子波，且有w(0，t)=w(t)。

(15)

將式(14)離散可改寫為矩陣向量形式，有

(16)

式中N為地震道的樣點數。

式(16)的加入縱向約束的L2范數最優化問題可用目標函數表示為

(17)

式中：W為時變子波矩陣；r為反射系數向量；d為地震記錄向量；μt為縱向約束因子，主要作用為調節數據殘差項(第一項)與正則化項(第二項)的權重，使得反演結果更合理。通過極小化目標函數，式(17)的解可表示為

r=(WTW+μtI)-1WTd

(18)

3 基于反演的全局約束吸收補償方法

假設地震剖面S(t，x)只有一個振幅為常數、斜率為p的線性同相軸，轉換到頻率域，可表示為

S(f，x)=A(f)exp(j2πfxp)

(19)

式中A(f)是同相軸子波的頻譜。假設道間隔Δx為常數，有x=(k-1)Δx，k=1，2，…，K，K為地震剖面的道數。對于任意頻率f，有

Sk(f)=Aexp(jθk)

(20)

式中θ=2πfpΔx。則第k道可用第k-1道表示為

Sk(f)=a1(f)Sk-1(f)

(21)

式中a1(f)=exp(jθ)。式(21)遞歸方程是一階差分方程，也稱為一階自回歸(AR)模型。通過該模型，可以沿著空間方向遞歸地預測有效信號。

類似地，假設在時空域有M條線性同相軸，則可以用一個M階的差分方程表示為

(22)

式中ai(f)表示信號預測算子，主要用于預測無噪聲的疊加復諧波。式(22)也可寫成預測誤差形式

(23)

式中：g0=1；gi=-ai。

在實際數據中，噪聲是不可避免的。假設地震數據由有效信號和噪聲疊加而成，即

Dk=Sk+Nk

(24)

將上式代入式(23)，可得

(25)

上式為自回歸滑動平均(ARMA)模型，其自回歸分量和滑動平均分量相等。

(26)

為了更加準確地預測有效信號，AR模型中的參數M一般需大于實際同相軸的條數，而且M的確定需要應用統計準則。

假設M=2，式(26)可以寫成線性方程組

(27)

其矩陣表達式為

D-Ya=N

(28)

為了最小化噪聲能量，建立目標函數

(29)

對目標函數求偏導，并令導數為零，有

YHYa=YHD

(30)

式中上標“H”表示共軛轉置。需要注意的是，矩陣YHY為Toeplitz矩陣，因此可用Levinson遞歸法快速求解。為了獲得穩定的濾波器系數，通常需要在Toeplitz矩陣的對角線上增加一個小的擾動μ1，則式(30)改寫為

(YHY+μ1Ι)a=YHD

(31)

那么待求的信號預測算子為

a=(YHY+μ1I)-1YHD

(32)

每一個頻率的預測算子是相互獨立的，將每個頻率的一維預測算子依次排列，可得一維預測算子

(33)

式中I為頻點數。將每個頻率的一維預測算子作用于含隨機噪聲的數據，可預測出有效信號

(34)

(35)

可見，該信號預測誤差算子除了在輸出位置的系數為-1以外，其他系數與信號預測算子的系數完全相同。地震信號的空間可預測性主要來源于地震同相軸的空間展布特征(即構造特征)，信號預測誤差算子直接反映了地震反射結構特征。需要注意的是，預測誤差算子作用于地震數據的結果不是有效信號，而是預測誤差。

現有的縱向約束吸收補償并未考慮到地震道橫向之間的關系。為了彌補這一缺陷，將地震信號的反射結構特征(即信號預測誤差算子)也作為一種約束引入到吸收補償算法。基于反演的全局約束吸收補償方法具體計算過程如下。

(1)通過解最小二乘的問題

(36)

(2)由預測算子可得到預測誤差算子

P=(a-M，…，a-2，a-1，-1，a1，a2，…，aM)

(37)

(3)通過反演

(38)

計算去噪地震數據。式中μ2為權衡參數。通過多次迭代求解式(38)，可得頻率域去噪數據。

(39)

式中μx為橫向約束因子。新的目標函數的解可表示為

(40)

本文提出的基于反演的全局約束多道吸收補償方法的橫向約束算子由原始的地震剖面求得，包含了原始數據的反射結構信息，能夠較明顯地提高吸收補償結果的穩定性和抗噪性，同時還保持了同相軸的橫向連續性。

4 應用

4.1 Marmousi模型試算

為驗證本文方法相對于縱向約束吸收補償方法的優越性，利用Marmousi模型進行測試。

圖1a是原始反射系數與主頻為30Hz的Ri-cker子波褶積得到的Marmousi模型模擬數據，可以看出，Marmousi模型構造復雜，既有平滑的同相軸，也有大型的背斜構造和一些小型斷裂。圖1b是Marmousi模型的衰減數據，其中Q設為50。圖1c和圖1d是信噪比(SNR)分別為5和20的加噪衰減數據，可以看出，由于地層吸收的影響，地震資料的分辨率降低，原有的構造反射特征變模糊。另外，隨機噪聲污染了地震反射信號。

分別采用傳統的縱向約束方法和本文的全局約束方法對圖1b所示的無噪衰減數據進行補償，結果如圖2所示。其中縱向約束吸收補償過程中μt=0.1，本文方法中μt=0.1，μx=1.0。比較圖2a與圖2b可以看出，在無噪情況下，縱向約束吸收補償方法和本文方法均能很好地恢復地震波能量，提高地震資料的分辨率。對比圖中紅色箭頭處可知，本文方法的結果橫向更加連續、平滑，構造結構更清晰，弱反射得到了較好的恢復，證明了本文方法的效果更好。圖3為無噪衰減數據兩種方法補償前、后單道頻譜對比，可以看出，兩種方法都有效補償了中高頻的振幅衰減，并且拓寬了原始數據的地震頻帶，但本文方法補償結果的頻帶略寬，包含更豐富的高頻有效信息。

圖1 Marmousi模型數據(a)原始模擬數據； (b)衰減模擬數據； (c)SNR=5的加噪衰減數據； (d)SNR=20的加噪衰減數據

圖2 無噪衰減數據兩種方法補償結果(a)縱向約束補償； (b)本文方法

對于含噪衰減數據，縱向約束吸收補償結果雖然恢復了地震波能量，提高了地震資料的分辨率，但是吸收補償結果的信噪比較低，地震同相軸的橫向連續性較差(圖4a)。本文方法的補償結果不僅能夠提高地震資料的分辨率，使得構造反射結構更清晰，同時也能保證補償結果的信噪比，橫向上更加連續(圖4b)，如圖中紅色箭頭所示，本文方法恢復出了更多的細節，表明本文方法比縱向約束吸收補償方法具有更強的穩定性和抗噪性。

圖5為SNR=5含噪衰減數據利用兩種方法補償結果的單道對比，可以看出，本文方法補償后與原始數據更相似，無論振幅還是相位，縱向補償效果均較好(紅色箭頭所示)。圖6為圖5結果的頻譜對比，可以看出，在含噪情況下，兩種方法均能補償中高頻的振幅衰減，并且拓寬了補償前數據的地震頻帶。此外，縱向約束補償結果的高頻部分略微強于本文方法補償結果，有可能是前者對高頻噪聲的過度放大所致。

圖3 無噪數據兩種方法補償前、后單道頻譜對比(a)第100道； (b)第103道； (c)第123道

圖4 SNR=5(左)和SNR=20(右)含噪數據兩種方法補償結果對比(a)縱向約束方法； (b)本文方法

以均方根振幅為衡量地震剖面橫向連續性，圖7為SNR=5含噪衰減數據兩種方法補償結果的橫向連續性對比，可見本文方法補償結果橫向上與補償前數據更加接近，橫向上抖動較小(紅色箭頭所示)，體現了本文方法在保持剖面橫向連續性和穩定性方面的優勢。

圖5 SNR=5含噪數據兩種方法補償結果單道對比(a)第100道； (b)第103道； (c)第123道

圖6 SNR=5含噪數據兩種方法補償結果單道頻譜對比(a)第100道； (b)第103道； (c)第123道

圖7 SNR=5含噪數據兩種方法補償結果的橫向連續性對比

4.2 實際數據應用

圖8a為二維實際地震疊加剖面，采用縱向約束吸收補償方法和本文方法的處理結果分別如圖8b、圖8c所示。可以看出，縱向約束補償方法雖然可以提高地震資料的分辨率，但補償后地震剖面中的噪聲被放大，信噪比較低；經本文方法處理后，地震資料的能量衰減得到有效補償，地震分辨率得到一定提升，而且補償后地震記錄信噪比較高，噪聲得到有效壓制，地震剖面的同相軸更連續。由此可見，本文提出的基于反演的全局約束吸收補償方法比傳統的縱向約束吸收補償方法具有更強的穩定性和抗噪性，其結果更利于后續的解釋和反演。

圖9為圖8數據的局部放大對比，可以看出，相較于原始剖面(圖9a)，本文方法補償剖面(圖9c)同相軸增多、更加清晰，提高了分辨率和信噪比，構造更加清楚可辨；相較于縱向約束方法補償結果(圖9b)，本文方法補償結果信噪比較高，橫向更連續。

圖10為實際地震剖面兩種方法補償后的單道頻譜對比，可以看出，兩種方法均能有效補償地震數據中高頻振幅衰減，一定程度上拓寬了數據的地震頻帶，其中縱向約束方法處理結果高頻部分相對較強，是高頻噪聲放大效應的體現。

圖8 實際地震剖面兩種方法補償結果對比(a)原始； (b)縱向約束方法； (c)本文方法

圖9 實際地震剖面兩種方法補償結果放大對比(a)原始剖面； (b)縱向約束補償剖面； (c)全局約束補償剖面

圖10 實際地震剖面兩種方法補償結果的單道頻譜對比(a)第100道； (b)第200道； (c)第300道

5 結論

(1)本文將傳統的單道吸收補償方法拓展到了多道吸收補償，由地震信號的可預測表達計算信號預測誤差算子，并將其作為橫向約束算子引入多道反演框架中，提出了一種基于反演的全局約束多道吸收補償方法。

(2)本文方法的目標函數表達式包含了縱向和橫向兩個方向上的約束信息，既考慮了縱向上的地震道稀疏分布信息，又考慮了橫向上的地震道間的連續性。

(3)模型數據和實際數據補償結果表明，本文提出的基于反演的全局約束多道吸收補償方法具有更強的抗噪性和更高的穩定性，可有效減小補償對于高頻噪聲的放大效應，提高了地震記錄分辨率和信噪比，同相軸橫向上更連續。

(4)由于基于反演的全局約束多道吸收補償方法屬于多道反演算法，運算量大，計算效率較低，以后可從數值解法和矩陣運算等方面進行優化以提高計算效率。