彭國(guó)民 徐凱軍* 杜潤(rùn)林 劉 展

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

1 引言

儲(chǔ)層物性參數(shù)是儲(chǔ)層預(yù)測(cè)和開(kāi)發(fā)的重要參數(shù)。地震勘探一直是油氣勘探和開(kāi)發(fā)的主要工具，一方面能夠提供地下精細(xì)構(gòu)造信息，另一方面利用巖石物理模型和地震數(shù)據(jù)反演能夠得到儲(chǔ)層的孔隙度、含油(氣)飽和度、滲透率、泥質(zhì)含量等物性參數(shù)，進(jìn)而對(duì)儲(chǔ)層的含油氣性進(jìn)行評(píng)價(jià)[1,2]。但是，利用地震資料估算儲(chǔ)層物性參數(shù)時(shí)，僅孔隙度和泥質(zhì)含量的估算相對(duì)可靠，對(duì)實(shí)際生產(chǎn)有一定的指導(dǎo)作用[3]。僅僅利用地震資料估算含油(氣)飽和度存在多解性，有可能導(dǎo)致錯(cuò)誤的油氣勘探開(kāi)發(fā)決策。例如，由于含氣飽和度的變化對(duì)地震波速度的影響不明顯，利用地震數(shù)據(jù)估算儲(chǔ)層的含氣飽和度有可能產(chǎn)生假象[4]。

海洋可控源電磁(MCSEM)勘探作為一種新的海上油氣勘探技術(shù)，能夠更好地識(shí)別含油氣薄層等微小地質(zhì)目標(biāo)體，被認(rèn)為是地震勘探的有效輔助手段[5,6]。在MCSEM勘探中，發(fā)射的電磁信號(hào)對(duì)高阻的含油氣儲(chǔ)層比較敏感，從儲(chǔ)層折射回來(lái)的電磁信號(hào)含有儲(chǔ)層的電阻率信息，利用儲(chǔ)層巖石物理關(guān)系可以將儲(chǔ)層電阻率與孔隙度、含流體飽和度聯(lián)系起來(lái)。因此，在根據(jù)地震資料確定儲(chǔ)層構(gòu)造格架的前提下，利用MCSEM數(shù)據(jù)能夠預(yù)測(cè)儲(chǔ)層內(nèi)的流體類(lèi)型和含油氣飽和度。

通過(guò)上述分析，聯(lián)合反演MCSEM和地震數(shù)據(jù)能夠得到更加可靠的儲(chǔ)層物性參數(shù)預(yù)測(cè)結(jié)果，降低單一數(shù)據(jù)反演的多解性。利用MCSEM數(shù)據(jù)和地震AVA數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)時(shí)，為了進(jìn)行這兩類(lèi)數(shù)據(jù)的正演，需要利用巖石物理關(guān)系將待估算的儲(chǔ)層物性參數(shù)與電磁屬性(電導(dǎo)率)、地震屬性(速度、密度)聯(lián)系起來(lái)。Hoversten等[7]、Abubakar等[8]和Gao等[9]在確定性框架下基于經(jīng)驗(yàn)性或統(tǒng)計(jì)性的巖石物理關(guān)系聯(lián)合反演MCSEM數(shù)據(jù)和地震AVA數(shù)據(jù)，直接估算儲(chǔ)層物性參數(shù)。而Hou等[10]和Chen等[11]則在貝葉斯框架下構(gòu)建聯(lián)合反演目標(biāo)函數(shù)，實(shí)現(xiàn)MCSEM數(shù)據(jù)和地震AVA數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)。與確定性聯(lián)合反演相比，基于貝葉斯框架的聯(lián)合反演能夠給出待反演參數(shù)的概率分布信息(均值、模式及置信區(qū)間)，可定量分析反演結(jié)果的不確定性。在中國(guó)，由于MCSEM勘探技術(shù)研究起步較晚，關(guān)于利用MCSEM和地震數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)的研究鮮有報(bào)道。杜潤(rùn)林[12]和徐凱軍等[13]利用地震數(shù)據(jù)和MCSEM數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)，取得了較好的反演結(jié)果。

模擬退火(SA)算法是一種常見(jiàn)地球物理反演全局尋優(yōu)算法，在求解重力與地震[14]、大地電磁與重力[15]、大地電磁與地震[16-18]數(shù)據(jù)聯(lián)合反演問(wèn)題中已取得較好的效果。

本文在應(yīng)用Archie公式和Gassmann方程分析電磁屬性和地震屬性對(duì)孔隙度和含水飽和度變化的靈敏度基礎(chǔ)上，提出了一種基于模擬退火優(yōu)化算法的MCSEM數(shù)據(jù)和地震AVA數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)方法，通過(guò)理論模型驗(yàn)證了本文方法預(yù)測(cè)儲(chǔ)層物性參數(shù)的優(yōu)越性。

2 靈敏度分析

Archie[19]提出了經(jīng)典的Archie公式，該公式把地層電導(dǎo)率σ與孔隙度φ和含水飽和度Sw聯(lián)系起來(lái)，奠定了利用測(cè)井資料定量評(píng)價(jià)油氣儲(chǔ)層飽和度的理論基礎(chǔ)。Gassmann[20]提出了流體替換方程，該方程建立了縱波速度vP、橫波速度vS、密度ρ與孔隙度φ、含油飽和度So、含水飽和度Sw、含氣飽和度Sg之間的關(guān)系，是儲(chǔ)層定量研究的重要成就。為了更好地理解反演結(jié)果，本文首先利用經(jīng)驗(yàn)性的Archie公式和Gassmann方程分析電磁屬性和地震屬性對(duì)孔隙度和含水飽和度變化的靈敏度。Archie公式的具體形式為

(1)

式中：a為與巖性有關(guān)的系數(shù)；m為孔隙度(膠結(jié))指數(shù)；n為飽和度指數(shù)。Gassmann方程的具體形式為

(2)

式中

(3)

其中：Ksat、μsat、ρsat分別表示流體飽和巖石的體積模量、剪切模量和密度；Kma、μma、ρma分別為巖石骨架的體積模量、剪切模量和密度；Ko、Kw、Kg分別為油、水、氣的體積模量；ρo、ρw、ρg分別為油、水、氣的密度；Co、Cw、Cg分別為油、水、氣的校正項(xiàng)；β為Biot系數(shù),一般是孔隙度的函數(shù)，本文采用Nur[21]提出的臨界孔隙度模型,定義為

(4)

式中φc為臨界孔隙度，超過(guò)該孔隙度固體變?yōu)閼腋☆w粒。

本文僅探討雙相介質(zhì)儲(chǔ)層(含油和水)的電磁屬性和地震屬性對(duì)孔隙度和含水飽和度變化的靈敏度。Archie公式中的參數(shù)設(shè)置為：a=1.0、m=1.2、n=2.4、σw=3.0S/m，Gassmann方程中的參數(shù)設(shè)置為ρma=2.56g/cm3、ρw=1.05g/cm3、ρo=0.75g/cm3、φc=0.4、μma=44GPa、Kma=37GPa、Kw=2.81GPa、Ko=0.75GPa、Co=1、Cw=1。

圖1a為由Archie公式得到的電導(dǎo)率隨孔隙度和含水飽和度的變化趨勢(shì)，從圖中可看出，地層電導(dǎo)率受含水飽和度的變化影響更明顯。圖1b～圖1d分別為由Gassmann方程得到的縱波速度、橫波速度、密度隨孔隙度和含水飽和度的變化趨勢(shì)，從圖中可看出，這三個(gè)地震屬性均對(duì)孔隙度的變化非常敏感，而對(duì)含水飽和度變化的敏感性較差，這表明利用地震數(shù)據(jù)反演能夠提供更可靠的孔隙度估計(jì)，而不能得到可靠的含水飽和度估計(jì)。

圖1 電導(dǎo)率(a)、縱波速度(b)、橫波速度(c)、密度(d)隨孔隙度φ和含水飽和度Sw的變化趨勢(shì)

由以上分析可知，單獨(dú)利用電磁或地震資料難以同時(shí)得到可靠的孔隙度和含水飽和度估計(jì)，這也說(shuō)明了利用MCSEM和地震數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)的必要性。

3 聯(lián)合反演方法

3.1 模型參數(shù)

本文僅考慮雙相介質(zhì)儲(chǔ)層(含油和水)，采用圖2所示的方式進(jìn)行模型參數(shù)化。儲(chǔ)層參數(shù)包括孔隙度φ和含油飽和度So，儲(chǔ)層上覆地層和下伏地層的參數(shù)包括電導(dǎo)率σ和縱波速度vP、橫波速度vS、密度ρ。由于待反演的參數(shù)為儲(chǔ)層的孔隙度φ和含油飽和度So，對(duì)于儲(chǔ)層的上覆和下伏地層的電導(dǎo)率σ和縱波速度vP、橫波速度vS、密度ρ可以事先通過(guò)單獨(dú)的電磁反演和地震反演獲得。

圖2 模型參數(shù)化

3.2 巖石物理關(guān)系

構(gòu)建巖石物理關(guān)系通常有三種方法：巖石物理學(xué)基本理論、經(jīng)驗(yàn)性的巖石物理模型和實(shí)際測(cè)井資料擬合，目前工業(yè)界廣泛應(yīng)用的是經(jīng)驗(yàn)巖石物理模型(例如Archie公式和Gassmann方程)和實(shí)際測(cè)井資料擬合[22-24]。經(jīng)驗(yàn)巖石物理模型的參數(shù)一般是利用測(cè)井資料通過(guò)線性回歸方法求得。對(duì)于實(shí)際測(cè)井資料擬合，通常先對(duì)測(cè)井資料進(jìn)行統(tǒng)計(jì)分析，然后假定儲(chǔ)層物性參數(shù)與電性參數(shù)、彈性參數(shù)滿足某種線性關(guān)系或多項(xiàng)式關(guān)系，同樣是利用測(cè)井資料通過(guò)多元回歸方法求得假定關(guān)系式中的系數(shù)。

本文利用前文靈敏度分析中給出的經(jīng)驗(yàn)性巖石物理模型(Archie公式和Gassmann方程)將孔隙度和含油飽和度與電性參數(shù)和彈性參數(shù)聯(lián)系起來(lái)。關(guān)于巖石物理模型參數(shù)的選取對(duì)聯(lián)合反演結(jié)果的影響，Gao等[9]通過(guò)僅其中變化一個(gè)或兩個(gè)參數(shù)研究巖石物理模型參數(shù)的誤差對(duì)反演結(jié)果的影響。

3.3 聯(lián)合反演算法

本文構(gòu)建的聯(lián)合反演目標(biāo)函數(shù)為

(5)

模擬退火(SA)是一種非線性的全局優(yōu)化算法，該方法不需要給定初始值，理論上一般可以收斂到全局極小值，且不需要求解靈敏度矩陣。該方法計(jì)算效率可能會(huì)稍低，但是對(duì)于一維反演，其計(jì)算效率能夠滿足實(shí)際應(yīng)用需求。本文采用SA算法求解聯(lián)合反演目標(biāo)函數(shù)，其包括三個(gè)關(guān)鍵步驟：接受概率、模型擾動(dòng)及退火方案[25]。在SA算法中，接受概率采用由Boltzmann概率分布給出

(6)

式中：k為Boltzmann常數(shù)，本文取值為1；T為當(dāng)前溫度； ΔE=E2-E1，E為能量值(即目標(biāo)函數(shù)值)。模型擾動(dòng)采用如下方式

(7)

(8)

Ti=T0exp(-i1/N)

(9)

式中：T0為初始溫度；N為一常數(shù)，本文取值為2[25]。

對(duì)于電磁場(chǎng)正演，采用擬解析解計(jì)算電磁場(chǎng)的分量[26]。地震AVA正演采用精確的Zoeppritz方程模擬地震AVA數(shù)據(jù)[27]

(10)

式中：RPP、RPS分別為縱、橫波反射系數(shù)；TPP、TPS分別表示縱、橫波透射系數(shù)；下標(biāo)“1”和“2”分別代表界面上、下變量；α1、α2分別為縱波的反射角和透射角；β1、β2分別為橫波的反射角和透射角。由Zoeppritz方程計(jì)算得到縱波反射系數(shù)RPP，再與給定的地震子波(通常為Ricker子波)進(jìn)行褶積得到合成地震數(shù)據(jù)。

4 模型試驗(yàn)

設(shè)計(jì)一個(gè)層狀模型驗(yàn)證MCSEM數(shù)據(jù)和地震AVA數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)的優(yōu)越性(圖3)。設(shè)定海水層厚度為1000m，海底下地層分為三大層，中間層為待反演的含油儲(chǔ)層，分為五小層，每層的厚度為20m，每層的孔隙度和含油飽和度見(jiàn)圖3(以海底為z方向的零點(diǎn))，其中第①層和第③層的電導(dǎo)率為0.147S/m、縱波速度為5508m/s、橫波速度為3703m/s、密度為2.4g/cm3。

MCSEM的發(fā)射源位于海底上方50m處，發(fā)射頻率分別為0.10、0.25、0.50、1.00、2.00Hz，海底布設(shè)10個(gè)電磁采集站，采集站間距為500m，且第一個(gè)采集站到發(fā)射源的水平距離為1000m。地震AVA道集的入射角分別為5°、10°、15°、20°、25°、30°、35°、40°，雷克子波的頻率為45Hz。對(duì)于MCSEM本文采用的是電場(chǎng)分量Ex與磁場(chǎng)分量Hy振幅的比值：Ex/Hy。

圖4a為單獨(dú)利用海洋可控源電磁數(shù)據(jù)反演的結(jié)果，從圖中可看出反演得到的含油飽和度比孔隙度更接近于真實(shí)值。圖4b為單獨(dú)利用地震AVA數(shù)據(jù)反演的結(jié)果，從圖中可看出地震數(shù)據(jù)單獨(dú)反演的孔隙度相較于飽和度稍好，這與前文靈敏度分析中得出的結(jié)論是一致的。圖4c為MCSEM數(shù)據(jù)與地震數(shù)據(jù)聯(lián)合反演的結(jié)果，從中可看出相比于圖4a、圖4b，聯(lián)合反演得到的含油飽和度和孔隙度估算值均更接近于真實(shí)值。

圖5為對(duì)聯(lián)合反演結(jié)果進(jìn)行正演計(jì)算得到的MCSEM響應(yīng)值與理論值的對(duì)比，圖6為對(duì)聯(lián)合反演結(jié)果進(jìn)行正演計(jì)算得到的地震AVA響應(yīng)值與理論值的對(duì)比。從圖5可看出，不同發(fā)射頻率的電磁數(shù)據(jù)的擬合誤差均比較小；從圖6可見(jiàn)，對(duì)于地震數(shù)據(jù)，小入射角的地震道的擬合誤差比較小，而大入射角的擬合效果較差。

圖3 理論模型示意圖

圖4 MCSEM數(shù)據(jù)反演(a)、地震數(shù)據(jù)反演(b)及聯(lián)合反演(c)結(jié)果

圖5 不同發(fā)射頻率的電磁數(shù)據(jù)擬合

為了說(shuō)明采用SA優(yōu)化算法的優(yōu)勢(shì)，通過(guò)設(shè)置兩組不同的初始值、利用基于梯度下降的Occam反演方法對(duì)該理論模型進(jìn)行反演，并進(jìn)行對(duì)比分析。首先給出接近于真實(shí)值的初始值，即孔隙度初始值分別設(shè)為0.15、0.25、0.10、0.20、0.20，含油飽和度初始值分別設(shè)為0.5、0.4、0.6、0.8、0.4，反演結(jié)果如圖7a所示，可看出兩種方法的反演結(jié)果都與真實(shí)值比較接近；設(shè)定與真實(shí)值相差較大的初始值，即孔隙度初始值分別為0.4、0.4、0.4、0.4、0.4，含油飽和度初始值分別為0.4、0.4、0.4、0.4、0.4，反演結(jié)果如圖7b所示，可看出Occam方法反演結(jié)果與真實(shí)值偏差較大，而SA算法仍能得到較好的反演結(jié)果。這是因?yàn)镾A算法不依賴于初始模型，因此能得到較好的反演結(jié)果，而Occam反演方法基于梯度下降，對(duì)初始模型的選取有較強(qiáng)的依賴性。

圖6 地震AVA數(shù)據(jù)擬合

圖7 設(shè)定目標(biāo)層不同孔隙度和含油飽和度初值時(shí)Occam和SA反演結(jié)果

5 結(jié)論和認(rèn)識(shí)

本文實(shí)現(xiàn)了MCSEM數(shù)據(jù)與地震AVA數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)。該方法利用經(jīng)驗(yàn)性的Archie公式和Gassmann方程建立儲(chǔ)層物性參數(shù)與電性參數(shù)和彈性參數(shù)之間的聯(lián)系，構(gòu)建歸一化的聯(lián)合反演目標(biāo)函數(shù)，并采用模擬退火算法進(jìn)行求解。模型試驗(yàn)表明基于模擬退火優(yōu)化算法的MCSEM與地震AVA數(shù)據(jù)聯(lián)合反演儲(chǔ)層物性參數(shù)具有優(yōu)越性。

由于需要利用巖石物理關(guān)系建立儲(chǔ)層物性參數(shù)與電性參數(shù)和彈性參數(shù)之間的聯(lián)系，不正確的巖石物理關(guān)系有可能導(dǎo)致錯(cuò)誤的反演結(jié)果，因此，在實(shí)際應(yīng)用中需謹(jǐn)慎處理巖石物理關(guān)系。本文中僅考慮雙相介質(zhì)儲(chǔ)層(含油和水)，實(shí)際情況是儲(chǔ)層內(nèi)含油、氣、水的情況遠(yuǎn)比雙相介質(zhì)復(fù)雜，需進(jìn)一步探討。