王保麗 藺 營 張廣智② 印興耀②

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580；②海洋國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測技術(shù)功能實(shí)驗(yàn)室，山東青島266071)

0 引言

地震反演是獲得地下介質(zhì)彈性參數(shù)、物性參數(shù)及含流體性質(zhì)參數(shù)的重要方法，是定量描述含油氣儲(chǔ)層的核心地震技術(shù)。根據(jù)不適定地球物理反問題求解方式的差異，地震反演分為確定性反演和隨機(jī)反演兩大類。基于最小化原理求解反問題最優(yōu)解或近似解的方法，通常稱為確定性反演，其估計(jì)結(jié)果是一個(gè)相對(duì)平滑的唯一解。依據(jù)統(tǒng)計(jì)學(xué)理論的隨機(jī)反演方法充分考慮了觀測數(shù)據(jù)及反演過程的不確定性，消除了確定性反演中的平滑性，能夠得到地下介質(zhì)隨機(jī)模型的多個(gè)實(shí)現(xiàn)，從而獲得多個(gè)反演結(jié)果[1-6]，在獲取反演最優(yōu)解的同時(shí)有效評(píng)價(jià)反演結(jié)果的不確定性。

地下介質(zhì)普遍存在非均質(zhì)特征，含油氣儲(chǔ)層的非均質(zhì)特征主要表現(xiàn)為儲(chǔ)層空間分布范圍變化大，因此儲(chǔ)層的彈性、物性及含流體性質(zhì)在空間上是變化而不均勻的[7-8]。通常使用統(tǒng)計(jì)學(xué)的隨機(jī)介質(zhì)模型描述地下復(fù)雜非均勻介質(zhì)，利用均值、方差、自相關(guān)函數(shù)等參數(shù)表征地下介質(zhì)空間變化的統(tǒng)計(jì)特性[9-12]。為精細(xì)研究地下介質(zhì)特性及其參數(shù)變化，需要借助隨機(jī)介質(zhì)模型描述地下介質(zhì)的小尺度非均勻特性，進(jìn)而有效刻畫地下油氣藏的細(xì)節(jié)信息。相比確定性反演，地震隨機(jī)反演方法融合了地層格架、測井、地震等多尺度信息，可以獲得更豐富的地下地層信息，使反演結(jié)果的分辨率更高[13-15]，可以描述非均質(zhì)復(fù)雜儲(chǔ)層。

近幾年，隨機(jī)反演方法在儲(chǔ)層預(yù)測、流體檢測及頁巖甜點(diǎn)預(yù)測等方面獲得了較好的效果。前人[16-19]利用隨機(jī)反演方法求取儲(chǔ)層彈性參數(shù)、物性參數(shù)以及流體因子。郭同翠等[20]利用疊前地質(zhì)統(tǒng)計(jì)學(xué)反演方法預(yù)測頁巖甜點(diǎn)的空間展布。Azevedo等[21]利用疊前地震隨機(jī)反演獲得縱橫波速度、密度及巖相數(shù)據(jù)。De Figueiredo等[22-23]基于巖石物理先驗(yàn)?zāi)Ｐ偷穆?lián)合貝葉斯反演獲得了彈性參數(shù)、物性參數(shù)等儲(chǔ)層特征參數(shù)。Mohamed等[24]結(jié)合地震屬性與疊后地震隨機(jī)反演獲得了砂巖相預(yù)測模型。

以貝葉斯推斷為基礎(chǔ)的地震隨機(jī)反演方法整合了未知模型參數(shù)的先驗(yàn)概率密度分布信息，通過似然函數(shù)建立模型參數(shù)與觀測地震數(shù)據(jù)間的關(guān)系，獲得待反演參數(shù)的后驗(yàn)概率密度函數(shù)[22,24-29]。在反演過程中，先驗(yàn)信息模型用于描述地下地質(zhì)體的空間變化特征，對(duì)于獲得準(zhǔn)確的反演結(jié)果非常重要，有利于縮小模型參數(shù)的求解空間。但常規(guī)的地震隨機(jī)反演方法以測井?dāng)?shù)據(jù)為硬數(shù)據(jù)，利用變差函數(shù)表征地下地質(zhì)體的空間展布特征[13-16,24-26,30-31]。然而，基于測井?dāng)?shù)據(jù)的隨機(jī)模擬不能描述復(fù)雜地下地質(zhì)體的空間變化特征，從而降低了非均質(zhì)儲(chǔ)層地震定量表征的穩(wěn)定性和可靠性。

為了充分利用已知地震、測井?dāng)?shù)據(jù)中蘊(yùn)含的地下地層空間結(jié)構(gòu)信息，本文基于隨機(jī)介質(zhì)理論，依據(jù)地震和測井?dāng)?shù)據(jù)得到精確表征地層空間變化特征的非均勻介質(zhì)特征參數(shù)，構(gòu)建非均勻介質(zhì)特征參數(shù)模型作為后續(xù)隨機(jī)反演的先驗(yàn)信息模型，通過優(yōu)化算法求解后驗(yàn)概率密度函數(shù)，形成了基于非均勻介質(zhì)特征參數(shù)的地震隨機(jī)反演方法。

1 方法原理

本文引入隨機(jī)介質(zhì)構(gòu)建貝葉斯隨機(jī)反演的先驗(yàn)?zāi)Ｐ?，利用混合型自相關(guān)函數(shù)表征非均勻介質(zhì)的空間分布特征，從已知地震、測井?dāng)?shù)據(jù)中估算縱向、橫向自相關(guān)長度及自相關(guān)角度等非均勻介質(zhì)特征參數(shù)，獲得反映地下復(fù)雜儲(chǔ)層小尺度空間變化特征的信息，為后續(xù)反演提供先驗(yàn)?zāi)Ｐ汀?/p>

1.1 隨機(jī)介質(zhì)基本原理

隨機(jī)介質(zhì)是利用統(tǒng)計(jì)方法描述地下介質(zhì)不同尺度非均勻性的介質(zhì)模型，一般可以用均值、標(biāo)準(zhǔn)差和自相關(guān)函數(shù)等一階、二階統(tǒng)計(jì)量表征地下介質(zhì)的空間隨機(jī)變化特征。隨機(jī)介質(zhì)模型可以表征地下復(fù)雜介質(zhì)的非均勻特性，有效反映含油氣儲(chǔ)層的小尺度變化特征。在二階平穩(wěn)假設(shè)條件下，隨機(jī)介質(zhì)模型定義為[9-11]

m(t,x)=m0(t,x)+σm(t,x)

(1)

其中

σm(t,x)=δm(t,x)f(t,x)

(2)

式中：m0(t,x)為模型的均值，可描述介質(zhì)的平均特性，表征介質(zhì)的大尺度非均勻性，x為橫坐標(biāo)，t為時(shí)間；σm(t,x)描述平均特性上的隨機(jī)擾動(dòng)，表征介質(zhì)的小尺度非均勻性；δm(t,x)一般為測井?dāng)?shù)據(jù)的標(biāo)準(zhǔn)差；f(t,x)為空間分布特征，是滿足二維自相關(guān)函數(shù)R(t,x)的均值為0、標(biāo)準(zhǔn)差為1的二維隨機(jī)序列。因此，可以使用均值為零的平穩(wěn)空間隨機(jī)過程表示地下介質(zhì)的小尺度非均勻性。

1.2 非均勻介質(zhì)特征參數(shù)

混合型自相關(guān)函數(shù)表征地下介質(zhì)的小尺度空間分布特征，通過引入粗糙度因子，混合型自相關(guān)函數(shù)綜合了高斯型和指數(shù)型的特點(diǎn)。高斯型自相關(guān)函數(shù)較好地描述了單尺度平滑的非均勻介質(zhì)。指數(shù)型自相關(guān)函數(shù)描述了隨機(jī)介質(zhì)的多尺度、自相似的特性，能更好地表征實(shí)際介質(zhì)的多尺度特性?；旌闲妥韵嚓P(guān)函數(shù)則更靈活地描述地下介質(zhì)，適用性更強(qiáng)[9-10]。因此，本文采用混合型自相關(guān)函數(shù)R(t,x)描述非均勻介質(zhì)的空間變化特征，即

(3)

式中：a、b分別為橫向(x方向)、縱向(t方向)自相關(guān)長度；θ為自相關(guān)角度；η為粗糙度因子，η=0為高斯型自相關(guān)函數(shù)，η=1為指數(shù)型自相關(guān)函數(shù)。

綜合式(1)～式(3)可知，描述地下非均勻介質(zhì)的特征參數(shù)包括均值、標(biāo)準(zhǔn)差、縱向和橫向自相關(guān)長度及自相關(guān)角度等，通過這些參數(shù)可以模擬產(chǎn)生滿足相應(yīng)的自相關(guān)函數(shù)且具有指定均值和標(biāo)準(zhǔn)差的隨機(jī)介質(zhì)模型。

1.3 非均勻介質(zhì)特征參數(shù)的估計(jì)及建模

為從實(shí)際地震數(shù)據(jù)中獲得地下地層空間結(jié)構(gòu)特征參數(shù)，首先建立地震記錄、地震子波與統(tǒng)計(jì)特征參數(shù)模型之間的函數(shù)關(guān)系[11]。根據(jù)隨機(jī)介質(zhì)理論可知，在連續(xù)介質(zhì)中，波阻抗Z(t,x)與反射系數(shù)r(t,x)的關(guān)系為

(4)

式中：Z0(t,x)為表征低頻分量的背景波阻抗；δz(t,x)為表征高頻分量的擾動(dòng)量。

由式(4)、線性褶積模型以及卷積的微分性質(zhì)可知，地震數(shù)據(jù)功率譜Ss(ω,kx)滿足

(5)

對(duì)SδZ(ω,kx)進(jìn)行傅里葉反變換得到隨機(jī)介質(zhì)的自相關(guān)函數(shù)R(t,x)，進(jìn)而估算出縱向、橫向自相關(guān)長度和自相關(guān)角度等參數(shù)。

利用估算的特征參數(shù)求取二維自相關(guān)函數(shù)及其功率譜，得到非均勻介質(zhì)的振幅譜，在區(qū)間[0，2π)產(chǎn)生均勻分布的隨機(jī)相位信息。綜合振幅譜和相位譜得到非均勻隨機(jī)介質(zhì)的頻譜函數(shù)，由傅里葉反變換得到二維隨機(jī)序列，對(duì)其進(jìn)行標(biāo)準(zhǔn)化，使其均值為零、標(biāo)準(zhǔn)差為1，再結(jié)合由測井?dāng)?shù)據(jù)得到的標(biāo)準(zhǔn)差和均值，從而構(gòu)建非均勻隨機(jī)介質(zhì)模型[12]。

2 非均勻介質(zhì)特征參數(shù)分析

由于均值表征介質(zhì)的大尺度平均特性，標(biāo)準(zhǔn)差反映偏離平均值的程度，為更好地突顯小尺度非均勻特性，本文給定具體的均值和方差，討論不同縱向、橫向自相關(guān)長度以及不同自相關(guān)角度對(duì)非均勻波阻抗介質(zhì)模型的影響。

圖1為不同縱向、橫向自相關(guān)長度及不同自相關(guān)角度的非均勻介質(zhì)模型。由圖可見：①自相關(guān)函數(shù)反映介質(zhì)在空間的相關(guān)程度，a展示了介質(zhì)在水平方向的相關(guān)范圍，反映介質(zhì)的橫向尺度；b展示了介質(zhì)在垂向的相關(guān)范圍，反映介質(zhì)的縱向尺度(圖1a、圖1b)。②θ反映隨機(jī)介質(zhì)模型的擾動(dòng)方向(圖1c、圖1d)，通過θ模擬地下介質(zhì)層位的方向變化。

圖1 不同縱向、橫向自相關(guān)長度及不同自相關(guān)角度的非均勻介質(zhì)模型(a)a=10,b=30,θ=0°；(b)a=30,b=10,θ=0°；(c)a=b=15，θ=0°；(d)a=b=15，θ=20°構(gòu)建的隨機(jī)介質(zhì)模型的網(wǎng)格數(shù)為400×400，網(wǎng)格間距為dx=dt=1，波阻抗均值為6×106kg·m-2·s-1，標(biāo)準(zhǔn)差為6×105kg·m-2·s-1

總之，通過分析不同尺度的非均勻介質(zhì)參數(shù)可知，縱向、橫向自相關(guān)長度和自相關(guān)角度等統(tǒng)計(jì)量描述了介質(zhì)彈性參數(shù)的空間擾動(dòng)，可以刻畫介質(zhì)的空間變化特性。因此，借助非均質(zhì)介質(zhì)特征參數(shù)構(gòu)建的模型反映了地下介質(zhì)的非均質(zhì)特性，為后續(xù)反演提供了更可靠的先驗(yàn)信息。

3 基于非均勻介質(zhì)參數(shù)的隨機(jī)反演

基于上述方法原理構(gòu)建表征地下介質(zhì)非均質(zhì)特性的先驗(yàn)信息模型，在貝葉斯理論框架下，聯(lián)合似然函數(shù)得到反演目標(biāo)函數(shù)，進(jìn)而利用優(yōu)化算法求解目標(biāo)函數(shù)，最終形成基于非均勻介質(zhì)特征參數(shù)的隨機(jī)反演方法。

3.1 貝葉斯理論

貝葉斯理論融合了待反演參數(shù)模型和已知數(shù)據(jù)的先驗(yàn)信息，通過似然函數(shù)轉(zhuǎn)化為后驗(yàn)信息，得到模型參數(shù)的概率估計(jì)，其后驗(yàn)概率密度函數(shù)為[32]

(6)

式中：p(m|d)為模型參數(shù)m的后驗(yàn)概率密度函數(shù)，表示觀測數(shù)據(jù)d已知時(shí)m的分布規(guī)律；p(m)為m的先驗(yàn)概率密度函數(shù)，表示d未知時(shí)m的分布規(guī)律；p(d|m)為似然函數(shù)，表示m與d的擬合程度；p(d)為d的分布情況。

3.2 目標(biāo)函數(shù)

本文以構(gòu)建的非均勻介質(zhì)模型為地質(zhì)統(tǒng)計(jì)先驗(yàn)?zāi)Ｐ?，并結(jié)合待反演參數(shù)與地震數(shù)據(jù)之間的正演關(guān)系構(gòu)建似然函數(shù)，進(jìn)而以貝葉斯理論為基礎(chǔ)[32]，得到表征最終反演結(jié)果的后驗(yàn)概率密度函數(shù)。為減少反演結(jié)果的不確定性，增加了非均勻介質(zhì)特征參數(shù)的約束項(xiàng)。因此，最終建立的目標(biāo)函數(shù)為

(7)

3.3 非?？焖倭孔油嘶饍?yōu)化算法

非?？焖倭孔油嘶?Very Fast Quantum Annealing，VFQA)算法是量子退火算法(Quantum Annealing，QA)的改進(jìn)算法，采用依賴于溫度的似Cauchy分布產(chǎn)生新的擾動(dòng)模型，可以加快算法的收斂速度。

為提高隨機(jī)反演的計(jì)算效率，本文采用VFQA算法擾動(dòng)更新模型數(shù)據(jù)，擾動(dòng)方式為

(8)

式中：mj∈[Ai,Bi]為修改后模型參數(shù)，[Ai,Bi]∈[Nx,Nt]為模型的搜索空間，j∈[Nx,Nt]，Nx×Nt為模型的大小；mi為修改前模型參數(shù)，i∈[Nx,Nt]；yi為由u產(chǎn)生的隨機(jī)變量，u∈[0,1]為均勻分布的隨機(jī)數(shù)；T為當(dāng)前溫度。

(9)

式中k、L均為可調(diào)參數(shù)。

3.4 基于非均勻介質(zhì)特征參數(shù)的隨機(jī)反演

由地震數(shù)據(jù)和測井?dāng)?shù)據(jù)估計(jì)非均勻介質(zhì)參數(shù)，基于估計(jì)的參數(shù)構(gòu)建隨機(jī)介質(zhì)模型作為反演的先驗(yàn)信息；結(jié)合地震數(shù)據(jù)與待反演參數(shù)之間的關(guān)系構(gòu)建似然函數(shù)；然后在貝葉斯理論框架下得到待反演參數(shù)的后驗(yàn)概率分布函數(shù)，進(jìn)而得到反演目標(biāo)函數(shù)；最后利用VFQA優(yōu)化算法優(yōu)化目標(biāo)函數(shù)得到最終的反演結(jié)果(圖2)。

圖2 基于非均勻介質(zhì)特征參數(shù)的隨機(jī)反演流程

4 模型試算與實(shí)例分析

4.1 模型測試

為驗(yàn)證上述方法的有效性，選用二維波阻抗模型(圖3)測試、分析，該模型具有很明顯的非均勻特征。在具體反演過程中，分別在第10、20、30和40道處抽取四口偽井，并給定主頻為30Hz的雷克子波，將計(jì)算得到的模型的合成地震記錄作為實(shí)際地震數(shù)據(jù)進(jìn)行反演。圖4為估計(jì)的b、a及θ。由圖可見，估計(jì)的三個(gè)參數(shù)與模型數(shù)據(jù)吻合較好，如在第250個(gè)采樣點(diǎn)處，θ(圖4c)與模型數(shù)據(jù)的角度變化趨勢一致，b(圖4a)和a(圖4b)也能反映兩個(gè)薄層的空間變化特征。

圖3 波阻抗模型數(shù)據(jù)在第250個(gè)采樣點(diǎn)附近有兩個(gè)薄層

圖4 估計(jì)的b(a)、a(b)及θ(c)

另外，考慮到模型的非均質(zhì)性，在參數(shù)估計(jì)時(shí)采用加窗處理，由地震數(shù)據(jù)目標(biāo)層段的非均質(zhì)特性確定滑動(dòng)窗口的尺度，得到的特征參數(shù)相當(dāng)于在窗口范圍內(nèi)對(duì)模型進(jìn)行了平滑處理。因此，估計(jì)的參數(shù)不能與模型數(shù)據(jù)完全一致，但已足夠反映模型數(shù)據(jù)的空間趨勢變化特征。

基于上述非均勻介質(zhì)模型的構(gòu)建方法構(gòu)建非均勻介質(zhì)先驗(yàn)信息模型(圖5)，其反映了模型數(shù)據(jù)(圖3)的空間展布特征。圖6為模型數(shù)據(jù)的隨機(jī)反演結(jié)果與常規(guī)稀疏約束脈沖確定性反演結(jié)果。由圖可見，隨機(jī)反演結(jié)果的分辨率更高(圖6a)。圖7、圖8分別為第20道和第40道的隨機(jī)反演結(jié)果、常規(guī)確定性反演結(jié)果與原始模型數(shù)據(jù)對(duì)比?？梢?，相對(duì)于確定性反演結(jié)果(圖8)，隨機(jī)反演結(jié)果在幅值和形狀上均與模型數(shù)據(jù)匹配較好(圖7)。圖9為QA算法、VFQA算法接受概率隨迭代次數(shù)的變化曲線。

圖5 非均勻介質(zhì)先驗(yàn)信息模型

圖6 模型數(shù)據(jù)的隨機(jī)反演結(jié)果(a)與常規(guī)稀疏約束脈沖確定性反演結(jié)果(b)

圖7 第20道(左)和第40道(右)隨機(jī)反演結(jié)果和原始模型數(shù)據(jù)對(duì)比

圖8 第20道(左)和第40道(右)確定性反演結(jié)果和原始模型數(shù)據(jù)對(duì)比

由圖可見，VFQA算法在迭代約100次時(shí)即可收斂(圖9b)，而QA算法則需要迭代近400次才可收斂(圖9a)。

圖9 QA算法(a)、VFQA算法(b)接受概率隨迭代次數(shù)的變化曲線

4.2 實(shí)例分析

本文選取中國M油田疊后地震數(shù)據(jù)(圖10)測試、分析，其中間部分反射同相軸變細(xì)，局部分辨率較高。圖11為對(duì)圖10估計(jì)的b、a及θ。由圖可見：在圖10的局部分辨率較高位置處，b的中間部分?jǐn)?shù)值偏小(圖11a)，這是由于非均質(zhì)體的縱向尺度減小所致；a反映了同相軸的橫向變化(圖11b)；θ反映了同相軸的傾斜變化，顯示地震同相軸先趨于平緩，然后隨著道數(shù)增加傾角逐漸增大(圖11c)。

圖10 疊后地震數(shù)據(jù)縱向采樣間隔為4ms，時(shí)間范圍為1.8～2.5s，共有200道，井位于第100道處

圖11 對(duì)圖10估計(jì)的b(a)、a(b)及θ(c)

首先計(jì)算測井?dāng)?shù)據(jù)的標(biāo)準(zhǔn)差等參數(shù)，然后由克里金插值建立測井?dāng)?shù)據(jù)插值模型作為均值，進(jìn)而利用上述建模方法估計(jì)特征參數(shù)構(gòu)建地質(zhì)統(tǒng)計(jì)隨機(jī)介質(zhì)模型(圖12)。該模型充分融合了已知測井和地震數(shù)據(jù)蘊(yùn)含的地下地質(zhì)信息，可以作為后續(xù)反演的先驗(yàn)信息模型。圖13為實(shí)際數(shù)據(jù)隨機(jī)反演與常規(guī)確定性反演結(jié)果。由圖可見，相對(duì)于常規(guī)確定性反演結(jié)果(圖13b)，基于非均勻介質(zhì)統(tǒng)計(jì)特征參數(shù)的隨機(jī)反演結(jié)果的分辨率較高，且橫向連續(xù)性也較好，更準(zhǔn)確地指示了含氣儲(chǔ)層的位置(圖13a黑色橢圓處)。

圖12 構(gòu)建的非均勻介質(zhì)先驗(yàn)信息模型

圖13 實(shí)際數(shù)據(jù)隨機(jī)反演(a)與常規(guī)確定性反演(b)結(jié)果

5 結(jié)論

地震隨機(jī)反演是描述地下復(fù)雜儲(chǔ)層的一種非常有效的方法，可以獲得高分辨率的儲(chǔ)層彈性參數(shù)數(shù)據(jù)體。本文提出了一種基于非均勻介質(zhì)特征參數(shù)的隨機(jī)反演方法，借助隨機(jī)介質(zhì)理論從已知地震和測井?dāng)?shù)據(jù)中提取描述儲(chǔ)層非均質(zhì)特性的縱向、橫向自相關(guān)長度和自相關(guān)角度等統(tǒng)計(jì)特征參數(shù)，為后續(xù)隨機(jī)反演提供可靠的先驗(yàn)信息模型，實(shí)現(xiàn)了波阻抗高分辨率反演。模型試算和實(shí)例分析表明，由于充分融合了已知地震和測井?dāng)?shù)據(jù)中蘊(yùn)含的地下信息，估計(jì)的非均勻介質(zhì)特征參數(shù)描述了非均質(zhì)儲(chǔ)層彈性參數(shù)的空間擾動(dòng)特性，反映了儲(chǔ)層的空間結(jié)構(gòu)特征，為隨機(jī)反演提供了可靠的地質(zhì)統(tǒng)計(jì)先驗(yàn)信息，高分辨率隨機(jī)反演的精度較高。但是，由于需要從地震數(shù)據(jù)中估計(jì)特征參數(shù)，地震數(shù)據(jù)的品質(zhì)影響反演效果。另外，由于加窗處理使估計(jì)的特征參數(shù)與地震數(shù)據(jù)不完全一致，也會(huì)影響后續(xù)反演精度。因此，需要進(jìn)一步研究利用地震數(shù)據(jù)估計(jì)特征參數(shù)的有效方法。