基于稀疏表示的增維疊前地震反演方法

楊 森 吳國忱*② 張明振 杜澤源 單俊臻 梁展源

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266580； ③中國石化勝利油田物探研究院，山東東營 257022)

0 引言

河流相薄儲層是目前油氣勘探的重點目標之一，該類儲層彈性參數空間變化劇烈，常規地震反演方法難以刻畫。為此，學者們在AVO理論的基礎上發展了FAVO反演技術，以減少薄層反演的不確定性，提高薄層反演精度。

Chapman等[1-3]通過動態等效介質模型得到不同頻率地震數據體所包含的儲層信息； Wilson等[4]結合瞬時譜分解與Smith-Gidlow反射方程，提出了頻變AVO反演理論； Wu等[5]利用偽Wigner-Ville分布算法實現了依賴頻率的AVO反演；張世鑫等[6]在疊后反演的過程中基于Shuey近似公式提出了縱波速度頻散梯度； 吳國忱等[7]通過正演模擬分析發現，疊前道集的地震屬性隨炮檢距和頻率變化同樣具有很大的差異性； 楊千里等[8]提出了基于貝葉斯理論的多尺度地震反演方法； 吳建魯等[9]構建了一維虛巖石物理模型，證明了中觀尺度不同流體間的相對流動是誘導地震波在地震頻帶衰減的主控因素； 李坤等[10]推導了時頻域FAVO直接反演目標方程，并提出了基于匹配追蹤譜分解的時頻域FAVO流體識別方法； 張繁昌等[11]在頻率域利用地震數據的正、余弦分量，研究了稀疏反射系數頻率域正、余弦分量協同反演方法。

FAVO反演技術不斷發展，頻率信息的充分發掘、利用是關鍵，即如何獲得高分辨率的頻譜是保障FAVO反演精度的基礎。傅里葉變換是頻譜分解的基礎，Gabor[12]利用可滑動的高斯窗提出了窗口傅里葉變換(STFT)；Morlet等[13]提出了連續小波變換(CWT)使時頻窗口可調；Stockwell等[14]省去了窗函數的選擇，提出了S變換，改善了窗寬固定的缺陷。然而這些時頻分解方法均屬于線性時頻分析技術，都會受到海森堡測不準原理的影響。為克服這種影響，學者們逐步建立了非線性的時頻分析技術。Mallat等[15]提出了匹配追蹤技術；Liu等[16]通過引入地震信號的瞬時信息，提出了動態匹配追蹤的快速算法；張繁昌等[17]提出了雙參數動態子波庫匹配追蹤快速算法，由于具自適應特征，在頻譜分解時具有較好的時頻分辨率，更適用于復雜的地震信號，可以更準確地刻畫地震信號的時頻特征，真實反映地下沉積體信息。

然而，盡管FAVO反演理論的不斷發展促進了薄層反演、流體識別等問題的解決，但因道集疊加而導致的干涉效應往往并未消除，稀疏表示技術的應用也主要偏向于時頻分析和地震屬性分析，這些因素制約了薄層刻畫的精度。因此，本文提出一種不疊加的思想，建立基于稀疏表示的增維疊前地震反演模式，有效削弱干涉調諧效應的影響，實現維度增加的高精度疊前地震反演。

需要說明的是，本文中增維既指反演輸入數據是多維的信息，又指對疊前道集的信息進一步挖掘的過程。利用譜分解方法將地震信號從炮檢距—時間剖面拓展為炮檢距—時間—頻率三維數據體，結合測網所組合的觀測系統實現了信息量的增維，增強了反演解釋描述能力。

1 方法原理

基于稀疏表示的增維疊前地震反演方法： 首先，以目的層為篩選條件對疊前道集進行角度優選，抽取有效的單角度地震資料；其次，基于稀疏表示理論，結合匹配追蹤算法與時頻分解方法對得到的單角度地震資料進行高分辨率的時頻分解，完成炮檢距—時間(O-T)域數據體向炮檢距—頻率—時間(O-F-T)域數據體的增維；再次，以每個角度的單頻資料為基礎，利用地震速度與Aki-Richard近似公式構建無井反演低頻模型[18]，基于貝葉斯理論直接建立波阻抗與地震數據之間的映射關系，利用非線性最優化算法對初始模型進行擾動，以前一個頻率反演的結果作為后一個頻率反演的約束，逐級反演，最終得到不同角度高分辨率增維反演結果；最后，通過聯合儲層的優勢響應頻帶計算截至該頻帶的增維反演結果，實現對非均質薄儲層的刻畫。

1.1 道集優選

疊加產生的干涉效應影響著疊后和疊前(部分角道集疊加)地震反演精度。如圖1所示，振幅、頻率均不相同的兩列波疊加產生的混合帯限子波常見于地震信號(圖1a)。頻率、振幅均相同的兩列波，當時差較大時，波峰與波谷相遇疊加振幅被削弱(圖1b)；當時差較小時，波峰與波峰疊加，振幅得到增強(圖1c)。因此利用疊加地震資料難以對目標準確描述。

圖1 兩列波疊加效應

(a)振幅、頻率均不相同； (b)振幅、頻率均相同，時差大； (c)頻率、振幅均相同，時差小。藍色曲線所示地震波為綠色和紅色曲線所示地震波的疊加

為了消除疊加效應，需對目的層段所包含的數據范圍進行優選，如圖2所示，對存在畸變的角度道集進行切除，逐一提取有效炮檢距或角度范圍內的地震數據組成若干單角度地震數據體，作為下一步時頻分解的資料基礎。

1.2 基于稀疏表示理論的高分辨率時頻分解

道集優選后所得的地震信號均為單角度的地震信號，避免了因為角度疊加而可能產生的干涉，但儲層厚度變化所導致的調諧效應依然存在。為得到高分辨率的頻譜分解結果，勢必要消除波形疊加影響，將子波進行還原、分解。

在地震信號處理領域，稀疏表示理論能夠有效促進反褶積、地震數據重建等線性反演問題的求解[19]，得到信號稀疏表示的過程被稱為稀疏分解。匹配追蹤算法通過構造過完備匹配子波庫、挑選與輸入信號最匹配的原子，可以實現信號稀疏表示的自適應分解[20]，所得的匹配子波相互獨立，因而可以利用該算法分離混合帯限子波(圖1a)。

圖2 道集優選

匹配追蹤算法的核心之一是對過完備匹配子波庫的構建。過完備匹配子波庫可以理解為由一系列匹配子波所構成的子波矩陣，也可以描述為不同原子的集合。匹配子波實質上是通過對基本子波ω進行時移、調制和相位變化而產生的。經典的匹配追蹤算法以Garbor字典構建過完備匹配子波庫，為了適應信號的不同特征，可以選擇Fourier字典描述頻率；構造Dirac字典刻畫位置[19]；創建wavelet字典確定位置和尺度[21]。不同字典的選擇可解決不同的問題，Ricker子波因波形簡單、收斂較快、延遲時間短而在薄互層的時頻分析中更具優勢。并且Ricker子波具有三個控制參量，延遲時間可控，子波庫波形豐富，與地震子波匹配較為靈活[22]。因此，本文以Ricker子波為母小波構建子波字典。

假設ω(t)為匹配子波庫中滿足條件的一個基本子波，令δ=(τ,f,φ)作為匹配子波的控制參數，則匹配子波表示為

ωδ(t)=ω(t-τ)ej[2πf(t-τ)+φ]

(1)

式中：τ為時移；f為頻率；φ為相位； j為虛數單位。通過控制參量τ、f、φ便可以完成對子波的匹配，其中時移因子固定匹配子波的位置；頻率因子把基本子波的能量集中在少數的時頻點附近。故過完備匹配子波庫Φ可以表征為

Φ={ωδ(t)}δ∈Γ

(2)

式中Γ表示為δ=(τ,f,φ)的集合，則一道帶寬有限的地震信號s(t)可描述為

(3)

式中：α=[α1,…,αk]稱為s(t)在完備匹配子波庫Φ中的表示系數；ωδk表示第k次的匹配子波；αk、τk、fk、和φk分別為控制第k次匹配子波ω的振幅、中心時間、波峰頻率和相位。

控制參量是通過掃描過完備匹配子波庫、實現信號全局尋優自適應分解而確定的。構建過完備匹配子波庫Φ后，從中選出與給定信號s最為匹配的匹配子波ω1，滿足內積〈s,ωi〉為匹配子波庫中所有原子與s內積中最大的一個

〈s,ω1〉>〈s,ωi〉i≠1ωi∈Φ

(4)

式中ωi為子波字典內的第i個子波原子。一次迭代后，信號s(t)被分解為

s(t)=〈s,ω1〉ω1+Rs1

(5)

式中Rs1是第一次迭代后的殘差信號。〈s,ω1〉越大，代表子波匹配度越高。

對殘差信號Rs1重復相同的匹配過程，即從Φ中篩選出與Rs1最為匹配的第二個原子。不斷重復上述過程，則迭代第k次(k≥0)有

Rsk=〈s,ωk-1〉ωk-1+Rsk-1

(6)

(7)

式中：Rsk為第k次迭代后的殘差信號；M為匹配子波的個數。

經過k次迭代分解后，原始信號s(t)可以用k個匹配子波的合成近似表示。當殘差信號Rsk足夠小，即信號的稀疏表示跟信號的近似程度足夠好或殘差信號滿足設定的閾值時，迭代停止。

如圖3所示，設計11個Ricker子波疊加合成的單道信號，并在0.8s附近以兩個相距小于λ/4的子波模擬出一個混合帯限子波。可以看出，合成信號被重新分解為11個主要的子波，而位于0.8s附近的混合帶限子波也被重新分解為2個互相獨立的子波，有效地還原了其所蘊含的真實信息。匹配追蹤算法具有良好的稀疏分解效果，能夠有效地保存信號的原始信息(圖4)，保障后續時頻分解的準確性。

在傳統的信號分析和處理中，傅里葉變換得到的振幅譜或相位譜反映地震信號的平均信息，而不能得到局部信息[23]。短時傅里葉變換受時窗長度限制，不同的時窗參數對結果產生較大的影響，無法同時兼顧時間分辨率和頻率分辨率。得益于匹配追蹤算法的優勢，在利用匹配追蹤算法將信號s(t)稀疏分解后，由于匹配子波集中的每一個匹配子波ωδn(t)都具有獨立性，因此可以對每一個匹配子波分別進行短時傅里葉變換。對低頻子波選取大的時窗、高頻子波選用較小的時窗，如此便可以克服常規短時傅里葉變換時窗固定不可調、分辨率單一的缺點，得到具有較高分辨率的時頻分布。基于匹配追蹤算法的STFT時頻分布為

圖3 基于匹配追蹤算法的稀疏分解結果

圖4 基于匹配追蹤算法的分解重構結果

(8)

式中gn(t)是與子波ωδn(t)對應的時窗。雖然在整體上短時傅里葉變換的時窗是不斷改變的，但是僅對某一確定的主頻時窗是固定的，因此可以極大地提升時頻分辨率。進而利用匹配追蹤算法得到的子波集，選擇不同的時頻分解方法就可以有效地提升分解精度，并且有效地避免干涉影響。

對圖3的合成信號選取不同方法進行時頻表征，結果如圖5所示。可以看出，短時傅里葉變換過程中不會產生交叉項，但因為時窗不可調，時頻表征的分辨率受到限制(圖5a)；經過匹配追蹤算法的改良，時頻分辨率得到了有效提升(圖5c)；相對于短時傅里葉變換，利用Wigner-Ville分布雖然分辨率有所提升，但產生了交叉干擾項，出現偽影(圖5b)；利用匹配追蹤算法的Wigner-Ville分布在保證時頻分辨率的同時，消除了交叉項的干擾(圖5d)。因此，相比于對信號直接進行時頻表征，基于匹配追蹤稀疏分解子波集解析得到的時頻分析結果能夠保證更高的時頻分辨率。

為了更好地驗證頻譜分解效果，進一步設計如圖6所示的復雜信號模型。對圖6b的復雜信號利用匹配追蹤算法進行稀疏分解，并以此為基礎在分解子波集內進行STFT變換。為更直觀地對比分解效果，改用譜圖形式進行對比，得到了如圖7b所示的頻譜圖，它與圖7a展示的原譜圖基本一致，證明了該方法的有效性。

譜分解技術拓寬了信號的維度，增加了所含的信息量。而基于稀疏表示理論(本文采用匹配追蹤算法)的時頻分解方法極大地提高了信號的時頻分辨率，保障了頻率信息的精確度，為后續道集多維分解與增維反演提供了基礎。

1.3 疊前道集的多維分解

疊前道集的地震屬性隨炮檢距和頻率變化，并且具有很大的差異性。因此，有別于傳統的AVO分析，疊前反演除了要考慮炮檢距，還應該考慮頻率的選擇。將上述的稀疏表示求解方法應用于疊前地震資料多維分解，能夠有效克服干涉效應，提高地震數據的時頻分辨率，且在分解過程中不會引入假頻，為增維疊前反演提供可靠的數據支持。

圖5 不同算法下的時頻表征

圖6 復雜信號模型

在疊前地震反演過程中，對于地球物理反演問題，模型參數和數據是以某種方式聯系起來的，數學上常構建為

d=Px+n

(9)

式中：d∈RN為含噪聲地震數據；P(RT→RN)為有界線性算子，T和N表示數據維度；x∈RT為待求解數據；n為噪聲。

觀測得到的地震數據d涵蓋了不同尺度頻率成分，可以表征為

(10)

式中：fi表示地震數據分解的頻率為iHz；fc為設定的分解截止頻率；dfi是頻率為fi的地震數據。

結合稀疏表示地震信號x=Φα以及稀疏正則化算法，不考慮噪聲干擾，基于L1范數的稀疏約束反演問題可表示為以下泛函問題

(11)

式中：χ為正則參數；d=Px=PΦα,α=[α1,…,αi]，因此不同頻率下的地震數據分量可描述為

(12)

式中：αi為字典的表示系數；Φδ(fi)={ωδ(fi)(t)}d(fi)∈Γ,0

圖8為A研究區a測線28°角道集地震數據按不同頻率尺度分解圖。可以看出：對于淺層(1.0～1.4s)河流相薄儲層(紅框所示)，分解所得頻率越高的單頻剖面，信息越豐富(圖8b～圖8f)；對于深部潛山(2.2～3.0s)，分解所得頻率越低的單頻剖面，潛山內部的地震響應越清晰。

圖7 模型分解效果對比

圖8 A研究區a測線稀疏時頻域地震數據多頻率分解圖

因此，不同深度的目標優勢響應頻段并不相同。對道集展開多維分解，在優勢響應頻段內可有效提升目標反演分辨率，進而提升儲層識別能力。

1.4 增維疊前地震反演方法

地震反演可分為疊后與疊前反演。疊后反演通過全角度疊加，雖然提高了資料品質，但弱化了儲層巖性及流體性質信息[24]。疊前反演以佐普里茲方程及其近似式為基礎[25]，采用部分角度疊加的策略，考慮子波在反演過程中的影響，得到了分辨率更高的反演結果，因此被廣泛應用。為消除干涉影響，提高對薄儲層的刻畫能力，發揮所拓展的頻率分量優勢，本文以貝葉斯數理統計理論為框架，開展增維疊前地震反演。

假定目標待反演的L組參數為m=[m1,m2,m3,…,mL]，為避免干涉效應所處理得到的增維道集數據集合

(13)

式中：θ={θ1,θ2,θ3,…,θn}表示目的層段所涵蓋的角度范圍，θn為篩選的最大角度；f={f1,f2,…,fc}為目標層的優勢響應頻帶。

根據貝葉斯理論，將數據—模型參數空間中的后驗概率密度(PPDF)表示為

(14)

式中dθf是入射角為θ、頻率為f的地震數據。

假設地震數據中的噪聲完全服從高斯分布，則似然函數可約定為

(15)

(16)

積分矩陣

(17)

(18)

對式(18)兩邊同取對數，得到目標泛函

J(m)=(Gm-dθf)T(Gm-dθf)+

(19)

式中μ和λ別為柯西約束項權值與趨勢背景約束項權值

(20)

由貝葉斯理論中模型參數的后驗概率分布構建反演目標函數，采用最大后驗概率解求取模型參數的方法，對目標函數(式(18))求導，整理得到反演方程

(21)

其中

(22)

表1 求解m流程

增維疊前地震反演流程如圖9所示。首先，從角道集中篩選、抽取n個角度的單角度地震數據，并對每一個單角度地震數據結合匹配追蹤算法按照fc個頻率進行頻率分解； 其次，依托速度場與Aki-Richard近似公式構建不同角度的低頻模型，以單角度地震數據所截取頻帶中的最低頻率數據d(θ1,f1)作為初始輸入，按照式(19)求解，得到該角度頻率的彈性阻抗結果EI(θ1,f1)； 再次，以該角度第二個頻率d(θ1,f2)的地震數據作為輸入，以第一個頻率d(θ1,f1)數據的彈性阻抗結果EI(θ1,f1)作為低頻模型，繼續按照式(19)求解，得到d(θ1,f2)的彈性阻抗結果EI(θ1,f2)； 依次循環，直到得到該角度截止頻率處數據d(θ1,fc)的反演結果EI(θ1,fc)，則該結果為此角度的分頻逐級反演結果； 最后，由此轉入下一個角度數據體，重復上述步驟，直到所有截取角度反演完成，得到n個增維彈性阻抗反演結果。

圖9 增維疊前反演流程圖

2 模型測試

為了驗證本文提出的增維疊前地震反演方法的有效性，采用經典的EAGE推覆體模型中的一條二維剖面(圖10)進行試算，以驗證增維反演結果的精度。

首先，基于稀疏表示理論對相應的褶積模型進行分解，得到不同頻率的地震數據；其次，依托模型速度場構建低頻模型，以分解得到的5Hz地震數據為初始輸入，得到如圖11a所示的彈性阻抗反演結果；再次，以圖11a所示的阻抗數據為模型輸入，以10Hz的地震數據作為輸入，得到10Hz的彈性阻抗反演結果；以此類推逐級反演，最終得到了各頻率的逐級反演結果。因數據量過多，本文只截取了5、25、35、45Hz的反演結果進行展示。由模型增維疊前反演結果(圖11)可以看出，可由低頻信息恢復模型的背景部分，高頻信息恢復出薄層小目標細節，逐級進行更新刻畫，可以有效地恢復局部薄層(200ms～300ms，黃色所示)細節。

圖10 彈性阻抗模型

圖11 不同頻率增維疊前地震反演結果

為了更直觀地分析反演效果，抽取單道反演結果(圖12)。可見本文方法反演所得的彈性阻抗曲線與實際彈性阻抗曲線相比具有較高的吻合度。

圖12 抽取的單道反演結果

3 實際資料分析

為了測試增維疊前地震反演方法的實際應用效果，選取A研究區進行實驗，目的層段為館陶組河流相薄儲層。目前已有鉆井63口，其中失利井30口。失利的主要原因之一為儲層薄、反演精度不高。該區Well-1井鉆遇厚度為7.18m與11.9m的兩套油層。

圖13為Well-1井井旁疊前角道集。首先對1200ms處目的層段按優勢響應角度進行篩選，選擇品質良好的角度范圍；其次對目的層的優勢頻段進行評估，選取儲層優勢響應頻帶。如圖13b所示，抽選道集中的28°地震數據具體分析，可見本區目的層1200ms處薄儲層段為高頻響應。進一步通過時頻分析(圖13c)處理，得出本區目的層段的優勢響應頻帶為45～80Hz。

為驗證增維疊前反演方法對于薄儲層的識別精度，針對研究區過Well-1、Well-2井的b測線抽選6°～28°的地震記錄作為觀測數據，利用稀疏表示結合短時傅里葉變換對該范圍內的單角度剖面以2Hz為間隔逐一進行高分辨率的譜分解。以此為基礎，在每一個角度內對單頻波從低到高重排，由低頻向高頻移動，在貝葉斯理論框架下，以前一個頻率的反演結果作為后一個頻率反演的約束，逐級反演，最終得到了增維疊前地震反演結果。

圖13 Well-1井井旁疊前道集及頻譜

圖14為b測線利用不同反演方法所得結果。為了展示反演效果，選取測線所過Well-1、Well-2兩井的測井解釋結果進行比對。可以看出，由于受到疊加影響，常規疊后反演剖面上黑色圓圈處的兩套油層(圖14a)混疊在一起難以區分。圖14b為16°～24°的疊加數據輸入所得到的疊前反演結果，與疊后反演相比，疊前地震反演結果信噪比雖有所降低，但大套層組間的界面變得更加清晰(黑色圓圈所示)，但黑色圓圈處的兩套油層(圖14b)仍難以分辨。單角度地震數據直接反演(圖14c)避免了地震資料疊加時產生的干涉效應，反演結果分辨率更高，黑色圓圈處的兩套油層已經隱約可見。而考慮了頻率因素的增維疊前反演結果 (圖14d)與測井解釋結果吻合更好，有效地區分出了1250ms(圖中黑色圓圈所示)處的兩套薄層，對于儲層的整體刻畫也更加精細。以上結果證明在不疊加思想下考慮頻率增維反演結果具有更高的描述精度，能夠更準確地對薄儲層進行識別和刻畫。

為了進一步驗證增維疊前地震反演方法的應用效果，截取目的層1250ms館陶組沿層切片(圖15)，分別采用不同的反演方法對有利儲層進行刻畫。疊后反演方法得到的彈性阻抗(圖15a)能夠大致刻畫出主河道的輪廓，但無法分辨分支河道；以部分角度疊加的疊前彈性阻抗反演(圖15b)在一定程度上加強了對分支河道的識別能力，但主河道上游(主河道自北向南流動)與切片左上部分(藍框內)分支河道的連續性仍較差； 28°單角度地震數據反演結果(圖15c)不僅能夠刻畫主河道的輪廓，并且主河道上游輪廓更加清晰，能很好地識別分支河道； 利用28°單角度地震數據按照增維策略分頻逐級反演切片(圖15d)上，主河道、分支河道可清晰識別，還揭示出在主河道下游因干涉效應而被掩蓋的一分支河道砂體(圖15d黑框處)，這在勘探中得到證實。

圖14 A研究區b測線不同波阻抗反演方法剖面對比

圖15 A研究區館陶組不同反演方法的沿層波阻抗切片

4 結論

本文引入不疊加的疊前反演理念，既不對角度進行疊加，同時又兼顧頻率因素的影響進行頻譜分解，有效地避免了干涉效應造成的信息畸變與缺失。基于貝葉斯反演理論，結合FAVO的反演思想，提出基于稀疏表示的增維疊前地震反演方法，得到了具有更高分辨率的反演結果，有效地提升了薄層識別能力。得到了以下認識。

(1)疊加會使部分區域的信息產生畸變，弱化反演精度。不疊加思想能夠有效地避免這種影響，并對真實信息進行還原。

(2)考慮頻率因素的反演可以有效地提升儲層描述能力。

(3)本文提出的增維高分辨疊前反演方法較傳統反演方法具有更高的分辨能力，能夠對河道砂體進行有效識別。