基于壓縮感知的反射系數域沿層L2范數約束去強屏蔽方法

張軍華 王 靜 王延光 劉立彬 李紅梅 王喜安

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②中國石油大學(華東)深層油氣重點實驗室，山東青島 266580； ③中國石化勝利油田分公司勘探開發研究院，山東東營 257015； ④中國石化勝利油田分公司物探研究院，山東東營 257022)

0 引言

利用地震資料研究薄層一直是難點，其中強屏蔽使下伏儲層信號變弱，導致難以精細描述與預測儲層。去除強屏蔽的方法包括四類：

(1)去相似背景的方法，通過分離地震記錄的背景信息和有用信息，從而去除與地層特性無關的信息，由汪恩華等[1]首先提出。趙銘海[2]、李軍亮[3]用于提取薄層反射特征。

(2)多子波變換方法，將地震道分解成多個不同頻率、振幅、相位的子波成分，然后根據實際需要重構子波，舍棄強屏蔽層對應的子波，從而達到去強屏蔽的效果。趙爽等[4]利用地震信號不同尺度的子波分解與重構去除煤層強屏蔽。張軍華等[5]應用多子波分解重構方法剝離強屏蔽層，提高了儲層預測精度。向龍斌等[6]利用多子波分解和地震旋回分析方法剝離火成巖強反射。金成志等[7]引進子波主成分分解思想，基于長、短旋回分解在Wheeler域提取并剝離了高頻短旋回強屏蔽信息。江馀等[8]提出變分模態分解方法消除強反射層，分辨率高、分解完備并避免了模態混疊效應。

(3)基于匹配追蹤的方法，基于投影追蹤、逐步遞推的小波算法，通過將每一次迭代得到的殘差信號不斷地進行最佳匹配實現。Wang[9]提出了基于自由多尺度的匹配追蹤去除強屏蔽方法，利用地震時頻譜分解剝離強屏蔽； 隨后，采用多道匹配追蹤算法解決強屏蔽的橫向連續性問題[10]。李海山等[11]成功地應用匹配追蹤算法分離煤系地層強屏蔽。張在金等[12]引入沿層信息，利用匹配追蹤技術在橫向定位強屏蔽層并剝離，再利用低頻伴影現象驗證剝離效果。

(4)基于反演的方法，基于已知先驗信息，將地震數據按照某種基函數分解而獲得反射系數。Levy等[13]提出了稀疏脈沖地震反演方法，可提高強反射之下薄層弱反射的分辨率。Chopra等[14]、Puryear等[15]提出了基于反射系數的薄層反演理論，可較好地分辨薄儲層。張銳等[16]引進壓縮感知理論在深度域反演薄層反射系數，提高了儲層識別精度。陳祖慶等[17]基于壓縮感知理論結合譜反演，利用地震資料部分譜信息恢復反射系數(稀疏反演)，改善了資料的分辨率。顧雯等[18]利用地震波形指示反演技術精確地預測了強屏蔽之下的薄儲層。宋維琪等[19]基于較大反射信號約束的觀測矩陣，聯合應用局部壓縮感知和區域壓縮感知的處理方法提高地震資料分辨率，從而增強弱信號。夏紅敏等[20]基于壓縮感知的譜反演，通過梯度投影稀疏重構實現地震譜反演，反演結果的信噪比高和橫向連續性好。

第四類方法是目前的研究熱點，尤其是壓縮感知[21-24]是一種稀疏采樣重構方法，具有很好的稀疏表征能力，根據地震資料可得到高分辨率的反射系數信息。常規的壓縮感知方法均基于反射系數稀疏反演提高地震資料分辨率，而對于強屏蔽層信息和儲層弱信息耦合疊加情況而言，效果并不理想，原因在于在剝離強反射層的同時往往會損失弱反射層的信息，從而不能有效預測儲層。

本文利用反射系數域分辨率高和不存在子波重疊的優勢，提出了基于壓縮感知的反射系數域沿層L2范數約束去強屏蔽方法，該方法基于壓縮感知理論，先根據時域反射系數域稀疏特性利用沿層信息分離強屏蔽與弱信號，再進行稀疏反演，最后將去強屏蔽后的反射系數與原子波褶積，從而獲得去強屏蔽層的高分辨率結果。理論模型測試和實際資料應用取得了較好效果，證明了方法的可行性。

1 算法

1.1 基本原理

設N×1維反射系數r(t)(t為時間)可以分解為n個采樣點處反射系數幅值αn和延時τn的線性組合[17]

(1)

式中：N為地震數據的采樣點個數； δ(·)為脈沖函數。引入強屏蔽層信息(一般通過人工拾取同相軸得到)，r(t)可以分解為

(2)

式中k表示強屏蔽層位置信息。

設子波為w(t)，噪聲為n(t)，則合成記錄為

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

(3)

對式(3)進行傅里葉變換，得

W(f)e-i2πfτk+N(f)

(4)

式中αk和τk分別為強屏蔽層位置處的反射系數幅值和時間位置。針對不同的頻率分量fm(m=1，…，M，M≥k)，將上式離散化，即

W(fm)e-i2πfmτk+N(fm)

(5)

記Amn=W(fm)e-i2πfmτn，由于e-i2πfmτn表示部分傅里葉矩陣，因此Amn為子波對角矩陣W(fm)與e-i2πfmτn的乘積組合。由于e-i2πfmτn和W(fm)均具有兩兩不相關特性，故Amn滿足有限約束條件(Restricted Isometry Property，RIP)，可以利用壓縮感知進行信號稀疏表示再重構。

將式(5)寫成矩陣形式

(6)

y=Ax′+ALx+N

(7)

根據壓縮感知理論，通過L1范數最小化策略，可得

(8)

式中λ為正則項，可控制L1范數所占比重。本文創新性地引入沿層信息ALx，在反演的同時利用高分辨率的反射系數分離耦合疊加的強屏蔽層信號與儲層弱信號，有利于強屏蔽層的拾取和剝離。由式(8)可見：第一項為沿層L2范數約束項或誤差估計項，表示利用沿層控制去強屏蔽層，也使迭代結果收斂； 第二項是稀疏約束項，促進結果的稀疏程度。

為了更好地約束層信息，引入窗函數，其寬度為T，則改進后的沿層對角矩陣為

La=diag(li)N×N

(9)

1.2 方法實現流程

由于式(8)中的x未知，通過引進迭代思路求解，構成以下目標函數

(10)

依次迭代求解x和x′。本文采用快速迭代閾值算法[25]求解目標函數，正則化初始值設為最大振幅的10%。

根據褶積公式，去強屏蔽后的反射系數x′與原始子波重構可以得到去強屏蔽的處理結果

s*=x′*w(t)

(11)

具體實現過程如下：

(1)輸入信號y、矩陣A、層位信息L，設定參數λ和正則化初始值；

(2)計算x和x′，利用快速迭代閾值算法求解x，再求解x′；

(3)將x′與子波褶積進行重構，s*=x′*w(t)；

(4)輸出最優化解s*。

本文方法沿層按列處理地震數據剝離強屏蔽層，即利用強屏蔽層信息(對應向量的一個位置)與反射系數的乘積剝離該列的強屏蔽層信息，特別在處理多道地震剖面時，采用三列加權平均方法提取子波，保證列與列之間的橫向連續性。

2 理論模型測試

2.1 一維模型測試

為了證明本文方法的去強屏蔽效果，首先建立一維含強屏蔽蓋層的薄互層正演模型(圖1a)，共201個采樣點，采樣率為2ms，選用30Hz零相位雷克子波作為震源進行正演模擬得到一維模型處理結果(圖1)。可見：在褶積模型中儲層弱反射被強屏蔽蓋層干擾，無法識別儲層(圖1a)； 去強屏蔽反射系數將強屏蔽蓋層與儲層分離(圖1c)，而且沒有損失儲層信息； 去強屏蔽地震記錄的波形后移，波谷能量減弱，增強和凸顯了儲層信息(圖1d黑色箭頭處)。圖2為一維模型去強屏蔽前、后時頻譜。由圖可見：去強屏蔽前，強能量團凸顯，掩蓋了下部儲層弱信息(圖2a)； 去強屏蔽后，強屏蔽層下部110ms處的殘留弱能量團對應儲層弱信息，說明本文方法在去除強屏蔽的同時保留了儲層弱反射信息，而且方法的抗噪性較強。

圖1 一維模型處理結果(a)褶積模型； (b)含30%噪聲地震記錄； (c)去強屏蔽反射系數； (d)去強屏蔽地震記錄本例是單道地震資料，沿層信息在100ms處，窗函數寬度為4ms

2.2 二維模型測試

為驗證本文方法的橫向連續性，設計了一個上部為強屏蔽蓋層的砂體反射系數模型(圖3a)，共51道，101個采樣點，采樣率為2ms，選用30Hz雷克子波作為震源進行正演模擬得到二維模型處理結果(圖3)。在實際處理過程中，采用逐列處理方式，并利用三列加權平均方法提取子波，以保證列與列之間的橫向連續性。可見，含30%高斯隨機噪聲地震記錄中強屏蔽蓋層反射和噪聲影響儲層弱反射的識別(圖3c)，去強屏蔽地震記錄中儲層弱反射清晰可見，也增強了其他層的反射(圖3e)。

圖2 含30%噪聲地震記錄(a)、去強屏蔽地震記錄(b)時頻譜

圖3 二維模型處理結果

為了考察本文方法的有效性，采用匹配追蹤算法對圖3c進行去強屏蔽層處理(圖4)。對比匹配追蹤方法去除的強屏蔽信息(圖4b)和圖3f發現，匹配追蹤算法去除的強屏蔽層位置出現偏差(圖4b紅色虛線框)，而本文方法剝離的強屏蔽層位置準確(圖3f)。分別對兩種方法的處理結果進行時頻分析(圖5)。可見：原始地震記錄30Hz時頻切片(圖5a)中強屏蔽層下方低頻能量團被覆蓋； 經本文方法去強屏蔽層后，強屏蔽層下方60ms處的弱能量團顯現，其余位置的能量團也相對增強(圖5b的黑色虛線對應弱反射位置)； 匹配追蹤方法去強屏蔽后，殘留強反射能量，沒有顯現強屏蔽層下方的弱能量團(圖5c紅色虛線框)。時頻分析結果說明本文方法較匹配追蹤方法去除強屏蔽層的效果更好，且不會損失弱反射信息。

圖4 匹配追蹤方法對圖3c去強屏蔽處理結果(a)匹配追蹤方法去強屏蔽結果； (b)匹配追蹤方法去除的強屏蔽信息

圖5 本文方法和匹配追蹤方法去強屏蔽層前、后的時頻切片(a)圖3c的30Hz時頻切片； (b)圖3e的30Hz時頻切片； (c)圖4a的30Hz時頻切片

3 實際地震資料應用

本文方法是基于壓縮感知的稀疏重構方法，在處理實際地震數據過程中，子波提取是非常重要的環節，因為子波的作用體現在以下兩方面：①基于壓縮感知對地震記錄進行稀疏反演，就是要消除子波信息干擾，從窄頻帶信息的地震記錄重構寬頻帶反射系數信息，子波會影響反演精度，而實際資料的原始子波是未知的，會造成反演誤差； ②利用反射系數進行拓頻處理，尤其是基于褶積模型拓展高頻信息，要以原始子波為基礎，通過對比頻譜信息范圍選取寬頻帶子波進行拓頻，因此子波會影響拓頻效果。本文針對實際資料采用多道加權復賽譜提取子波技術，該方法基于同態反褶積技術，對多道復賽譜序列疊加后求平均，反射系數的平均復賽譜會趨于零，而子波的復賽譜保持不變。根據該原理，通過濾波分離反射系數和子波，從而提取子波。

將本文方法應用于A區，圖6為過井原始地震剖面和本文方法去除強屏蔽層處理結果。由圖可見，去強屏蔽層后，強屏蔽層位置振幅減弱，強屏蔽層下方出現弱反射信息(圖6b紅色箭頭處)，并且與井曲線對應較好。圖7和圖8分別為去強屏蔽層前、后的20Hz和34Hz時頻切片。由圖可見：在去強屏蔽層后的20Hz時頻切片上，強反射同相軸下方能量團增強，出現低頻伴影現象(圖7b紅色箭頭處)； 在去強屏蔽層后的34Hz時頻切片上，強反射同相軸下方也出現儲層弱能量團信息(圖8b紅色箭頭處)。圖9為去強屏蔽層前、后沿層能量半時屬性切片。由圖可見，去強屏蔽層前沿層能量半時屬性值基本不變，無法展示下伏灘壩砂體弱反射儲層(圖9a)； 經去強屏蔽層處理，井點及整個儲層的相關性變好，可以有效識別有利儲層區域(圖9b)。

圖6 過井原始地震剖面(a)和本文方法去除強屏蔽層處理結果(b)選取的地震數據時間范圍為2000～2300ms，采樣率為2ms。紅線表示強屏蔽層位置，綠線為另一個標志層位置，下同。沿層信息通過人工解釋層位獲得

圖7 去強屏蔽層前(a)、后(b)的20Hz時頻切片

圖8 去強屏蔽層前(a)、后(b)的34Hz時頻切片

圖9 去強屏蔽層前(a)、后(b)沿層能量半時屬性切片(上、下開窗10ms)

4 結論

當上覆地層存在厚層烴源巖等特殊層位時，在地震剖面上會呈現強振幅特征，此類強屏蔽層往往會掩蓋下伏儲層的信息，不利于準確預測儲層。為了有效去除強屏蔽層的影響，本文提出了基于壓縮感知的反射系數域沿層L2范數約束去強屏蔽方法。該方法基于壓縮感知理論，首先根據時域反射系數域稀疏特性，利用沿層信息分離強屏蔽層與儲層，再進行稀疏反演； 最后將去強屏蔽后的反射子波與原子波褶積，從而獲得去強屏蔽層的高分辨率結果，其優勢在于高分辨率的反射系數能夠分離強屏蔽層與儲層的信息，有利于拾取和剝離強屏蔽層。理論模型測試和實際地震資料應用的結果表明，文中方法的計算步驟以及參數選取較簡單，處理后的地震數據可以有效凸顯強屏蔽層以下的儲層弱反射信息，證明了本文方法的有效性和和可行性。

在實際應用中，以下兩點需要關注：

(1)為了更好地進行層位約束，引進了一個窗函數，其寬度會影響強屏蔽層的剝離效果，在實際應用中，要具體分析地震數據，一般根據強屏蔽層的厚度確定窗函數寬度；

(2)本文方法的難點在于沿層準確拾取強屏蔽層位置，目前方法大多基于人工追蹤層位，存在一定的不確定性，而且難以保證精度。因此在應用文中方法時，要盡量結合井位信息，確保層位拾取的準確性。