李 坤 印興耀

（①中國石油大學（華東）地球科學與技術學院，山東青島266580；②海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術功能實驗室，山東青島266071）

0 引言

疊前地震反問題作為儲層地球物理學的核心內容，是解決復雜油氣儲層特征描述及油氣判識的重要理論與方法，通常包括波動方程反演（層析成像、全波形反演等）、AVO（Amplitude Variation with Offsets，振幅隨炮檢距變化）、EI（Elastic Impedance，彈性阻抗）反演［1-4］。然而，由于帶限地震數(shù)據(jù)缺失低頻和高頻信息，且受到噪聲的干擾，疊前地震反問題表現(xiàn)出較強的不適定性，因此疊前地震反問題必須通過其他途徑（如鉆、測井數(shù)據(jù)、油藏建模數(shù)據(jù)等）補充缺失的模型信息。正則化和概率化反演（以貝葉斯理論為主）是緩解疊前地震反問題不適定性最重要的兩種手段［5］。然而，隨著油氣藏類型日益復雜化，僅僅依靠單一的模型最優(yōu)解，已經(jīng)難以滿足地球物理勘探人員對復雜油氣儲層巖性和含流體性質的精細描述的要求，故模型參數(shù)最優(yōu)解、不確定度及置信區(qū)間等勘探風險評估逐漸引起儲層地球物理學家的關注。

在地震反演的貝葉斯框架中，模型參數(shù)的先驗信息可以通過隨機變量的先驗概率密度分布函數(shù)（Prior Probability Density Function，Prior PDF）定量表征。其認為地下介質模型參數(shù)是隨機的，則地震反問題是指在多類型觀測數(shù)據(jù)集的協(xié)同約束下匹配隨機變量的隨機實現(xiàn)過程［5］。目前，以貝葉斯推斷為基礎的反問題中，最大后驗概率解（Maximum A Posteriori Estimation，MAP）［5-13］、后驗PDF的顯式解［14-20］及后驗PDF的隨機解（以馬爾可夫鏈蒙特卡洛模型為主，MCMC）［21-27］是應用最廣泛的三種概率解。貝葉斯MAP解是通過貝葉斯公式計算得到模型參數(shù)的后驗概率密度分布，以最大化后驗PDF獲取模型參數(shù)的最優(yōu)解。后驗概率的顯式解是指模型參數(shù)的后驗統(tǒng)計特征（如均值和協(xié)方差等）可以通過顯式表達式定量表征，以高斯型概率模型為應用典型。隨機解是指通過隨機采樣的方式對后驗PDF進行隨機采樣，利用大量隨機樣本點逼近理論后驗PDF的求解方式。Hansen等［28］闡述了兩類觀測數(shù)據(jù)約束下高斯概率密度模型的地質統(tǒng)計模擬問題以及貝葉斯線性反問題的顯式解。Grana等［14-15］在混合先驗分布的基礎上，采用序貫模擬方法對后驗PDF的顯式解進行單點模擬，實現(xiàn)了縱波阻抗、離散巖相的同步預測，該方法由于未引入模型參數(shù)的低頻信息，需要在模擬相對縱波阻抗后補償?shù)皖l背景。Lang等［16］在高斯概率模型的基礎上，通過后驗概率模型的顯式解實現(xiàn)了疊前地震縱、橫波阻抗和密度的隨機反演方法。Li等［27］聯(lián)合馬爾科夫鏈蒙特卡洛模型、混合高斯概率先驗信息，提出了巖相驅動下的疊后地震隨機反演方法，改善了疊后地震反演分辨率的同時，實現(xiàn)了“離散相態(tài)”與“彈性參數(shù)”的協(xié)同預測。

Hastie等［29］基于貝葉斯推斷理論框架，系統(tǒng)闡述了混合概率模型在“離散變量”反演與分類領域的意義?；旌细怕誓Ｐ驮谪惾~斯推斷中是實現(xiàn)疊前地震“連續(xù)彈性參數(shù)”和“離散相態(tài)”協(xié)同預測的基礎?；诖耍疚囊罁?jù)巖性驅動的模型參數(shù)統(tǒng)計先驗信息，把地層巖性（如砂、泥巖）與模型參數(shù)和觀測地震數(shù)據(jù)聯(lián)系在一起，提出了一種基于混合概率模型的時頻聯(lián)合域疊前地震概率化彈性阻抗反演方法；綜合利用時域、頻域地震、低頻整合先驗信息及已知模型點四類數(shù)據(jù)集實現(xiàn)對彈性阻抗、離散巖相的序貫協(xié)同預測；此外，將非線性邊界修正算法引入疊前地震流體因子的直接提取過程，旨在提高疊前地震反演的穩(wěn)定性，克服彈性參數(shù)提取過程中易出現(xiàn)“超界解”的問題。通過理論模型測試和實際勘探實例證實疊前地震概率化彈性阻抗反演方法的穩(wěn)定性和可行性。

1 方法原理

1.1 時頻聯(lián)合域地震正演模型構建

不同變換域地震反演是通過考慮不同域內地震響應與合成地震數(shù)據(jù)間的匹配程度以搜索最優(yōu)模型參數(shù)的過程［5］。頻率域反演相比時間域反演具有頻率分量自動解耦、頻率選擇自由、多分量逐級尋優(yōu)及分辨率高的特點。因此，將混合域卷積模型應用于地震數(shù)據(jù)的正演

式中：S（ω，θi）、W（ω，θi）分別為入射角為θi的地震頻譜和子波頻譜，ω為角頻率；rPP（t，θi）為基于Biot-Gassmann孔隙彈性理論的Russell線性AVO近似方程［30-33］；j為虛數(shù)單位；t為時間；N（ω）為噪聲頻譜。式（1）可以重組為矩陣形式

假設η為利用實際測井數(shù)據(jù)構建的彈性阻抗低頻先驗數(shù)據(jù)，λL為低頻背景的正則化因子，D為正則化對角矩陣。聯(lián)合式（3）與η，可將正演模型改寫為

式中：H為輸入的三種類型已知數(shù)據(jù)集；P為正演核矩陣；Cωω、Ctt和Cηη分別是和η的協(xié)方差矩陣；Cωt、Cωη分別是、η之間的協(xié)方差矩陣；Ctη是與η之間的協(xié)方差矩陣；CH為三種數(shù)據(jù)集H的聯(lián)合協(xié)方差矩陣。

1.2 混合概率模型驅動的彈性阻抗反演方法

混合概率模型是多個“單峰”PDF的線性疊加，旨在量化待反演模型參數(shù)的“多峰”分布特征，是“單峰”PDF的推廣形式［5］。待反演模型參數(shù)的先驗概率密度分布往往受到儲層巖性的影響，儲層巖性不同將導致模型參數(shù)的統(tǒng)計特征差異［5，15，19，27］。具有“多峰”分布特征的混合概率模型是實現(xiàn)“離散型”—“連續(xù)型”模型參數(shù)協(xié)同預測的關鍵。高斯型PDF是數(shù)學、物理及工程等領域中最普遍的統(tǒng)計分布，具有線性變換性質，易于求解地球物理反問題的顯式解［16-18，28-29］。因此，本文將高斯混合PDF作為彈性阻抗反演模型的先驗PDF。其中，彈性阻抗的先驗均值、方差及高斯混合PDF中的峰態(tài)數(shù)量由實際測井數(shù)據(jù)的解釋結果和統(tǒng)計特征獲得。

針對貝葉斯線性地震反問題，若目標工區(qū)中儲層巖性類別的數(shù)量為M，將式（6）代入式（A-8）～式（A-9）（附錄A），則模型最優(yōu)解可用后驗概率分布的顯式解進行量化，推導得到三類數(shù)據(jù)集（即和η）協(xié)同約束下的第k個高斯概率分量pk（m|H）的顯式解，其后驗均值為

式中：λk＝p（zk）為第k種離散巖性zk的先驗概率；表示均值為、協(xié)方差為的高斯PDF，分別為已知數(shù)據(jù)集H的訓練均值和訓練協(xié)方差。針對實際工區(qū)可以通過測井解釋數(shù)據(jù)統(tǒng)計獲取，其中Fk是第k種巖性出現(xiàn)的頻數(shù)。

通過混合后驗均值和后驗協(xié)方差的顯式表達式（式（7）、式（8）），即可計算三類已知數(shù)據(jù)集H約束下的待反演模型參數(shù)最優(yōu)解。為了更好地分析待反演模型的不確定性，采用基于單點實現(xiàn)的序貫模擬算法（Sequential Simulation）（圖1）對后驗PDF進行隨機采樣。此外，將先前的模擬點集ms作為貝葉斯地震反演中下一個模擬點mi的約束條件。因此，針對第k個高斯概率分量pk（m|H），在待模擬點位置的后驗均值可改寫為

圖1 利用序貫模擬算法的時頻聯(lián)合域地震反演示意圖［16-17］

在貝葉斯推斷中，混合先驗概率模型搭建了儲層彈性阻抗與離散巖性之間的關系，在彈性阻抗、觀測數(shù)據(jù)等特征屬性的協(xié)同約束下，即可獲得儲層巖性的后驗概率密度分布；然后，依據(jù)最大后驗概率解獲取儲層巖性的最優(yōu)解。基于此，從被優(yōu)選的目標分量中隨機采樣得到巖相約束下的彈性阻抗數(shù)據(jù)。Hansen等［28］和Lang等［16］提出了2類條件數(shù)據(jù)約束下的地震反問題的線性解。基于貝葉斯分類，可以得到多類數(shù)據(jù)集協(xié)同約束的第k個高斯概率分量的后驗權重

為序貫模擬過程中條件數(shù)據(jù)的先驗均值。待模擬點mi的可以用來判斷、識別離散的儲層巖性。當離散巖性分類（優(yōu)選混合概率模型中的高斯分量）之后，待反演模型參數(shù)（即彈性阻抗）將從被選擇的高斯分量中隨機采樣獲得。值得注意的是，當訪問下一個模擬點時，先前的模擬點被視為已知數(shù)據(jù)，且已知數(shù)據(jù)的鄰域（矩形或圓形鄰域）選擇對于提高地震序貫模擬的計算效率具有重要意義，鄰域半徑通常設置為一個波長至兩個波長之間。

1.3 疊前地震Gassmann流體參數(shù)反演方法

復雜油氣儲層的飽和巖石為干巖石骨架與多種流體填充物構成的多孔雙相介質，宏觀巖石模量受干巖石骨架和孔隙流體的綜合影響。Russell等［30-31］基于Biot-Gassmann孔隙彈性理論推導了巖石飽含流體情況下的線性地震AVO反射系數(shù)方程

式中：θ為入射角；f是Gassmann流體項，反映巖石孔隙流體的作用；μ為巖石剪切模量，反映巖石骨架的作用；ρ為飽和巖石密度；為干巖石骨架的縱橫波速度平方比；是飽和巖石的縱橫波速度平方比；Δlnf、Δlnμ和Δlnρ分別為Gassmann流體項、剪切模量和密度的地層反射率。

由式（15）可知，彈性阻抗反演比疊前地震AVO反演具有更高的計算效率和穩(wěn)定性［1，4］，因此提出了基于模型正則化和邊界約束的改進彈性阻抗反演方法，以預測表征油氣儲層含流體性質的Gassmann流體項等彈性參數(shù)。與式（15）對應的流體彈性阻抗方程可以改寫為

式中：F、μ和ρ分別表示井旁道實測Gassmann流體項、剪切模量和密度構成的向量；n為重采樣后的測井數(shù)據(jù)樣點數(shù)。

采用阻尼最小二乘法（或共軛梯度法）求解式（18），獲得a′（θi）、b′（θi）和c′（θi）。令彈性參數(shù)的相對值表示為Rf、Rμ、Rρ和REI，即

式中N為單道地震數(shù)據(jù)的樣點數(shù)。為了實現(xiàn)非線性方程的線性化，對式（16）等號兩邊取自然對數(shù)，并將式（19）代入，得

式中：Aθi＝diag［a（θi）］N；Bθi＝diag［b（θi）］N；Cθi＝diag［c（θi）］N；diag［·］表示對角矩陣運算；Q為入射角的數(shù)量。

考慮到Gassmann流體項、剪切模量及密度的大尺度先驗背景可以從實際測井數(shù)據(jù)中獲得，故將彈性參數(shù)的已知先驗信息用L2范數(shù)表示，并引入到Gassmann流體項等參數(shù)的反演中，提高彈性參數(shù)反演的可靠性，則彈性參數(shù)反演的目標泛函和更新量分別為

式中：Δξf、Δξμ和Δξρ均為先驗數(shù)據(jù)的迭代殘差；先驗信息的正則化算子Df、Dμ、Dρ是單位矩陣；ξf、ξμ和ξρ是彈性參數(shù)的先驗背景；λf、λμ和λρ分別是Gassmann流體項、剪切模量和密度的先驗約束權重。式（21）～式（22）未考慮待反演模型的下限和上限信息，然而實際地震數(shù)據(jù)受到噪聲的干擾，反演結果往往會超出模型邊界的解（固有不穩(wěn)定性引起）。為了改善該問題，引入非線性邊界轉化算法。假定當前迭代次數(shù)為y，則中間變量的更新量為

2 理論模型測試

為了驗證混合概率模型驅動的疊前地震概率化反演可行性，利用重采樣后的測井數(shù)據(jù)設計了一維理論模型（圖2）。圖2展示了Gassmann流體項、剪切模量、密度、小角度彈性阻抗的先驗概率密度分布及理論離散巖相。采用30Hz零相位Ricker小波、2ms采樣間隔和精確Zoeppritz方程，合成了信噪比S/N＝1情況下的AVO地震數(shù)據(jù)，AVO地震數(shù)據(jù)的入射角分別為10°、20°和30°（圖3）。

第一步是基于貝葉斯混合概率模型的地震彈性阻抗反演處理。不同角度（10°、20°和30°）的彈性阻抗與離散巖相的估計結果如圖4～圖6所示。從圖4～圖6可以看出，50次模擬實現(xiàn)（灰線）的均值解（紅線）與理論模型（綠線）的擬合度較高。不同角度情況下，50次離散巖相判識的最大后驗概率解（圖4c、圖5c和圖6c的左半部分）均與理論巖相（圖4c、圖5c和圖6c的右半部分）具有高度一致性。與之對應的，50次模擬實現(xiàn)的巖相誤差主要出現(xiàn)在地層界面上及先驗均值之間的模型參數(shù)過渡區(qū)，該現(xiàn)象產(chǎn)生的根本原因是這些區(qū)域的模型參數(shù)在不同后驗PDF中的概率是相近的。因此，基于貝葉斯混合概率模型的地震彈性阻抗反演方法可以有效地實現(xiàn)彈性阻抗和離散巖相的概率化預測，為后續(xù)的邊界約束疊前地震概率化彈性阻抗反演奠定了數(shù)據(jù)基礎。

圖2 實際測井數(shù)據(jù)

圖3 S/N＝1情況下入射角為10°（a）、20°（b）和30°（c）的合成地震數(shù)據(jù)

圖4 S/N＝1情況下小角度（10°）彈性阻抗與離散巖相的估計結果

圖5 S/N＝1情況下中角度（20°）彈性阻抗與離散巖相的估計結果

圖6 S/N＝1情況下大角度（30°）彈性阻抗與離散巖相的估計結果

第二步是基于彈性阻抗的Gassmann流體項、剪切模量等彈性參數(shù)的同步反演。圖7展示了S/N＝1情況下反演得到的f、μ和ρ。由圖可見，本方法的50次模擬實現(xiàn)（灰線）的彈性參數(shù)預測結果被有效地限制在上、下限（棕色虛線）范圍內，與理論模型高度一致；f預測結果的相對誤差為12.8%，μ的相對誤差為17.9%，證明了該邊界約束彈性阻抗反演算法在參數(shù)提取方面的可行性；存在噪聲干擾的情況下，f的估計結果受隨機噪聲的影響程度較小，μ的反演精度比f低，與理論分析相一致，主要是由于μ對地震AVO反射系數(shù)的貢獻度低導致的。

為了闡述基于彈性阻抗的Gassmann流體項、剪切模量等彈性參數(shù)的同步反演方法與常規(guī)反演方法之間的差異，圖8展示了不同反演方法得到的Gassmann流體項反射率的頻譜計算結果。由圖可見，最小二乘確定性反演方法僅能獲取模型參數(shù)的“最優(yōu)解”，反演結果的頻帶范圍主要與輸入地震數(shù)據(jù)的頻寬有關，高頻信息不能得到有效恢復（圖8b）；本文方法反演結果的頻譜（圖8c、圖8d）比常規(guī)最小二乘反演方法具有更寬的頻帶，尤其是在模型參數(shù)的高頻信息恢復方面，該概率化反演方法可以更有效地恢復地下介質高頻信息（圖8d）。

圖7 S/N＝1情況下50次隨機模擬實現(xiàn)的彈性參數(shù)及其相對誤差估計結果

圖8 S/N＝1情況下不同反演方法得到的Gassmann流體項反射率的頻譜分析對比

3 實際資料處理

以中國華南地區(qū)Well-1井含油儲層的地震勘探為例，驗證本文混合概率模型驅動的疊前地震彈性阻抗反演方法的實用性。

研究區(qū)砂巖儲層主要為含油層與油水同層。由圖9可以看出，小角度彈性阻抗可以較好地區(qū)分優(yōu)質砂巖儲層（含油砂巖、含油水砂巖，圖9a橢圓框）與泥巖（干層）（圖9a）；相比中角度、大角度彈性阻抗，小角度彈性阻抗可以更好地區(qū)分砂巖儲層（圖9b）。不同巖性情況下，小角度彈性阻抗的統(tǒng)計特征（先驗均值、協(xié)方差）差異更大。然而，中角度與大角度彈性阻抗在砂巖儲層和泥巖干層位置的先驗均值卻比較接近，故在概率化推斷過程中將存在較強的不確定性。由于小、中、大角度彈性阻抗具有相似的統(tǒng)計特征“拐點”，因此本文將彈性阻抗巖相判識結果的“交集”作為目標工區(qū)巖相預測的依據(jù)。

圖9 彈性阻抗的巖石物理參數(shù)交會與統(tǒng)計

圖10 儲層含流體性的巖石物理參數(shù)交會與統(tǒng)計

圖11 小角度（5°～16°）疊加的地震剖面、彈性阻抗及離散巖相的反演結果

由圖10可見，Gassmann流體項較好地區(qū)分了油層、油水同層及泥巖干層，同時Gassmann流體項、密度等參數(shù)服從混合高斯概率密度分布，表現(xiàn)為“多峰”的概率統(tǒng)計特征?；诖?，可將Russell近似公式作為“流體相”驅動疊前地震反演的地球物理映射關系，即模型參數(shù)為表征孔隙流體的Gassmann流體項、干巖石骨架的剪切模量、巖石密度及孔隙含流體類型（流體相）。

由圖11～圖13可以看出：①由于概率化地震反演的不確定性，單次模擬結果中仍存在序貫隨機采樣的干擾（圖11c～圖11d、圖12c～圖12d、圖13c～圖13d）。因此，地震概率化反演往往需要多次隨機實現(xiàn)的均解或最大條件概率解（圖11d、圖12d、圖13d）。對比圖11c和圖11d，10次實現(xiàn)的均值解可以有效地抑制單次實現(xiàn)中的隨機干擾。②井旁道的彈性阻抗估計結果（紅線）和砂泥巖模型的巖性判識結果與實際測井彈性阻抗計算結果、砂巖儲層解釋結果高度一致，可以較準確地識別出砂巖油藏的發(fā)育位置（圖11e～圖11h、圖12e～圖12h、圖13e～圖13h）。

根據(jù)井旁道彈性阻抗的巖石物理參數(shù)交會與統(tǒng)計特征分析（圖9），角度彈性阻抗可以有效地識別砂巖儲層，且將小、中、大3個角度彈性阻抗的巖相判識結果的“交集”作為砂泥巖中儲層的最終判識結果（圖14a），可以降低巖相判識的不確定性，與實際砂巖油層解釋結果吻合較好。圖14b展示了邊界約束疊前彈性阻抗反演預測的Gassmann流體項f，從圖中可以看出，Gassmann流體項在油砂儲層位置存在明顯低值異常，且砂泥巖－巖相與流體識別結果在預測位置上保持了較高的一致性，證明了本文混合概率模型驅動的疊前地震概率化反演方法在巖性預測和儲層流體識別中具較好的實用性，應用前景廣泛。

圖12 中角度（16°～26°）疊加的地震剖面、彈性阻抗及離散巖相的反演結果

圖13 大角度（26°～30°）疊加的地震剖面、彈性阻抗及離散巖相的反演結果

圖14 疊前地震概率化彈性阻抗反演的巖相和Gassmann流體項f的估算結果

4 結論

本文依據(jù)地下介質待反演模型參數(shù)的先驗概率服從混合型概率密度模型，提出了基于時域、頻域地震、低頻整合先驗信息及已知模型數(shù)據(jù)點四類條件數(shù)據(jù)集協(xié)同約束下的混合概率模型驅動的疊前地震時頻聯(lián)合域彈性阻抗反演方法。主要結論如下。

（1）混合概率模型是實現(xiàn)“離散型”—“連續(xù)型”模型參數(shù)協(xié)同預測的關鍵，只有當不同離散相的待反演參數(shù)具有不同的統(tǒng)計特征時，才能夠實現(xiàn)地震反演的離散相態(tài)判識。高斯混合PDF與高斯PDF的乘積仍可以表征為標準的混合高斯PDF，此時，“連續(xù)彈性參數(shù)”、“離散相態(tài)”的后驗均值與后驗協(xié)方差具有顯式表達式，易于采用序貫模擬隨機采樣算法求解。

（2）相比常規(guī)確定性反演方法，序貫模擬算法是通過逐點遍歷的方式實現(xiàn)模型參數(shù)的隨機采樣，先前模擬點將作為下一個模擬點的約束條件，有助于實現(xiàn)“連續(xù)彈性參數(shù)”與“離散巖相”的概率化協(xié)同判識，便于實現(xiàn)和理解，使多數(shù)據(jù)協(xié)同約束的地球物理反問題變得清楚易懂。

（3）理論測試和實際應用驗證了該方法的可行性及穩(wěn)定性。應用實例表明，混合概率模型驅動的疊前地震彈性阻抗反演的巖相和概率多重實現(xiàn)，與實際測井數(shù)據(jù)和油藏解釋結果具有較高的一致性，有助于復雜巖性油氣藏的儲層刻畫和流體直接檢測，應用前景廣泛。

附錄A

貝葉斯理論是在當前觀測樣本數(shù)據(jù)的情況下，結合對模型參數(shù)的先驗認知，計算得到模型參數(shù)后驗概率密度分布（PDF）的統(tǒng)計學習方法。

式中：p（m|H）表示后驗概率密度分布；p（H|m）為觀測數(shù)據(jù)的似然概率密度分布；m為表征該概率模型的未知參數(shù)；p（m）被稱為未知參數(shù)的先驗信息；p（H）為已知觀測數(shù)據(jù)的PDF，與θ的取值無關［34］。若待反演模型參數(shù)服從高斯混合PDF

式中：M為高斯混合概率模型中所包含概率分量的數(shù)量；；n為單道地震數(shù)據(jù)的樣點數(shù)；為模型均值；為模型參數(shù)協(xié)方差。在高斯似然函數(shù)的前提下

式中：m為被選擇頻率分量的個數(shù)；H為觀測地震數(shù)據(jù)；P為正演核矩陣；CH為協(xié)方差矩陣。利用貝葉斯公式（式（A-1））推導后驗PDF（貝葉斯公式的分母與待估計模型參數(shù)無關，可忽略分母項）

將式（A-6）改寫為標準的高斯型PDF形式

則式（A-7）所示的第k個標準型后驗PDF分量的后驗均值的顯式解為

式（A-8）、式（A-9）即描述了高斯混合概率模型驅動下的標準后驗PDF的顯式解。