基于Hampel三截尾函數(shù)的儲層彈性和物性參數(shù)同步反演

李志勇,張家樹,蔡涵鵬,胡光岷

(1.電子科技大學(xué)通信與信息工程學(xué)院,四川成都611731;2.西南交通大學(xué)信號與信息處理四川省重點(diǎn)實(shí)驗(yàn)室,四川成都610031)

?

基于Hampel三截尾函數(shù)的儲層彈性和物性參數(shù)同步反演

李志勇1,張家樹2,蔡涵鵬1,胡光岷1

(1.電子科技大學(xué)通信與信息工程學(xué)院,四川成都611731;2.西南交通大學(xué)信號與信息處理四川省重點(diǎn)實(shí)驗(yàn)室,四川成都610031)

常規(guī)的地震反演方法通常假設(shè)噪聲信號服從高斯分布，以此構(gòu)建出的目標(biāo)函數(shù)在特定地區(qū)反演效果好，但其普適性較小。假設(shè)地震噪聲信號和巖石物理模型誤差服從高斯分布且?guī)в屑有悦}沖噪聲，引入Hampel三截尾函數(shù)構(gòu)造能夠同時壓制高斯噪聲和脈沖噪聲的反演目標(biāo)函數(shù)，利用牛頓法求解目標(biāo)函數(shù)得到儲層彈性和物性參數(shù)反演的迭代公式，并通過閾值的選擇適應(yīng)噪聲分布不同的地區(qū)。該算法利用Hampel三截尾函數(shù)構(gòu)建加權(quán)矩陣,能在一定程度上自適應(yīng)調(diào)節(jié)地震數(shù)據(jù)、巖石物理約束和先驗(yàn)信息之間的權(quán)重，消除地震和巖石物理噪聲的統(tǒng)計性誤差造成的反演不穩(wěn)定性，最終得到穩(wěn)定的儲層彈性和物性參數(shù)反演結(jié)果。通過理論模型測試和實(shí)際資料應(yīng)用顯示了算法的穩(wěn)定性和可靠性。

AVO反演;巖石物理反演;聯(lián)合反演;高斯噪聲;脈沖噪聲;Hampel三截尾函數(shù)

20世紀(jì)80年代,勘探工作者發(fā)現(xiàn)地震數(shù)據(jù)的反射振幅隨炮檢距的變化可以用于油氣檢測,從而發(fā)展了利用縱波AVO屬性進(jìn)行構(gòu)造解釋和巖性識別的AVO分析技術(shù)[1]。在AVO分析技術(shù)的基礎(chǔ)上發(fā)展起來的AVO反演方法,可從地震疊前資料中獲得反映地下儲層信息的縱、橫波速度和密度等信息,為儲層預(yù)測與烴類檢測提供量化依據(jù)。隨著AVO反演方法逐漸趨于成熟,人們開始利用地震、測井和巖心等數(shù)據(jù)的互補(bǔ)特性,從井間地震數(shù)據(jù)中反演地層巖石的彈性參數(shù)和物性參數(shù)[2-3]。與單純的AVO反演相比,這種多源信息聯(lián)合反演的方法不僅可以減少反演多解性的問題,還能夠同時獲得多種儲層參數(shù)，因此受到了學(xué)術(shù)界和工業(yè)界的廣泛關(guān)注。

多源信息聯(lián)合反演研究多于一個的目標(biāo)函數(shù)在給定區(qū)域上的最優(yōu)化,求取目標(biāo)函數(shù)最小或最大值所對應(yīng)的一組模型參數(shù),并把這組模型參數(shù)看成是滿足這些條件的“最佳解”。顯然,目標(biāo)函數(shù)不同,最佳解也就不一樣。為了求解方便,絕大多數(shù)方法都假設(shè)噪聲信號服從高斯分布[4-7],利用統(tǒng)計性方法和確定性方法求解目標(biāo)函數(shù)。前者將觀測數(shù)據(jù)和模型參數(shù)都視為隨機(jī)變量,用統(tǒng)計的方法來確定解所服從的概率分布,得到的是滿足觀測數(shù)據(jù)的解具有多大概率;而后者認(rèn)為觀測數(shù)據(jù)和模型參數(shù)都是確定量,因而用數(shù)理方程或代數(shù)方法求解反演問題,根據(jù)觀測數(shù)據(jù)求得模型參數(shù)。BOSCH等[8]在貝葉斯框架下將地震似然函數(shù)、巖石物理似然函數(shù)和先驗(yàn)信息結(jié)合,構(gòu)建服從高斯分布的后驗(yàn)概率密度函數(shù),利用統(tǒng)計性方法(蒙特卡羅采樣)獲得儲層的孔隙度和波阻抗。統(tǒng)計性方法存在不完全依賴于初始猜測、反演過程理論上不會陷入局部極值的優(yōu)點(diǎn)。但其缺點(diǎn)也十分明顯,就是計算工作量大、效率低,難以完全滿足大規(guī)模地球物理反演的要求。BOSCH等[9-10]在貝葉斯框架下利用確定性方法(最小二乘法)對信噪比相對較高的地震疊后資料進(jìn)行了孔隙度和波阻抗同步反演。李志勇等[11]在BOSCH等研究成果的基礎(chǔ)上嘗試?yán)肁VO正演模擬和Gasmann方程建立儲層物性參數(shù)與地震疊前觀測數(shù)據(jù)之間的聯(lián)系,利用最小二乘法同步反演獲得縱橫波速度、密度、孔隙度、含水飽和度和泥質(zhì)含量。研究發(fā)現(xiàn),利用最小二乘準(zhǔn)則可以方便地計算最速下降方向和修正步長,并且具有不錯的收斂速度。但對于信噪比相對較低的疊前地震數(shù)據(jù),使用最小二乘準(zhǔn)則會使誤差數(shù)據(jù)得到相同的權(quán)重,嚴(yán)重影響反演結(jié)果的穩(wěn)定性。另外,由于速度模型與埋藏深度通常存在正比關(guān)系,而儲層物性參數(shù)服從0～1的截斷分布,因此，建立儲層物性參數(shù)與彈性參數(shù)之間的巖石物理模型時,誤差同樣不完全服從高斯分布[12]。LI等[13]采用自適應(yīng)加權(quán)系數(shù)來調(diào)整地震似然函數(shù)、巖石物理似然函數(shù)和先驗(yàn)信息之間的權(quán)重,提高了多參數(shù)同步反演方法的穩(wěn)定性,但無法從根本上消除脈沖噪聲信號的影響。

對于脈沖噪聲信號,TAYLOR等[14]使用L1范數(shù)建立目標(biāo)函數(shù),得到了比傳統(tǒng)的L2范數(shù)更魯棒的解。然而,L1范數(shù)在殘差接近于0的時候存在奇異點(diǎn),容易導(dǎo)致反演算法的不穩(wěn)定。在TAYLOR等研究成果的基礎(chǔ)上,很多學(xué)者研究了基于L1范數(shù)的變形算法。如BUBE等[15]采用迭代加權(quán)算法來最小化L1/L2混合范數(shù)。GUITTON等[16]、呂曉春等[17]采用Huber范數(shù)建立目標(biāo)函數(shù),在殘差較小時采用L2范數(shù),在殘差較大時采用L1范數(shù),避免了L1范數(shù)的奇異性帶來的優(yōu)化困難。Huber范數(shù)在一定程度上減小了脈沖噪聲對反演算法的影響,但由于閾值點(diǎn)為無窮大,也沒有從根本上消除脈沖噪聲的影響。隨著脈沖噪聲比例的增加,Huber范數(shù)將會失效。考慮到Biweight范數(shù)能夠?qū)Υ笥谠O(shè)定閾值的噪聲信號予以忽略,JI[18]將Biweight范數(shù)引入到地震反演中,提高了算法對脈沖噪聲的抑制能力。劉洋等[19]提出了針對非高斯噪聲的地震疊前非高斯反演概念和思想,構(gòu)造了能同時壓制高斯噪聲和脈沖噪聲的混合范數(shù)并將其作為目標(biāo)函數(shù)。岳碧波等[20]假設(shè)地震噪聲信號服從非高斯α穩(wěn)定分布,提出了基于非高斯α穩(wěn)定分布的最小Lp范數(shù)地震反演方法。另外,考慮到廣義極值概率密度分布函數(shù)在一定情況下具有逼近任意概率密度分布函數(shù)的能力,ZHANG等[21]提出了基于廣義極值分布的AVO反演方法。

誤差協(xié)方差矩陣是L2范數(shù)的一個重要統(tǒng)計性質(zhì),其正確估計是進(jìn)行數(shù)據(jù)完美匹配的前提。由于協(xié)方差矩陣對脈沖噪聲的存在非常敏感,因此脈沖噪聲可能導(dǎo)致對協(xié)方差矩陣估計的錯誤,解決這一問題最直接的辦法就是采用魯棒性估計方法修正協(xié)方差矩陣。本文假設(shè)地震噪聲信號和巖石物理模型誤差服從高斯分布且?guī)в屑有悦}沖噪聲,引入Hampel三截尾函數(shù)構(gòu)建重加權(quán)算子來對協(xié)方差矩陣進(jìn)行修正。Hampel三截尾函數(shù)結(jié)合了Huber范數(shù)和Biweight范數(shù)的優(yōu)點(diǎn),不僅能夠減小脈沖噪聲對協(xié)方差矩陣估計的影響,而且還能將超過閾值的噪聲信號完全排除,提高算法對高斯噪聲和脈沖噪聲的抑制能力。

1 技術(shù)方法

1.1 巖石物理模型

地震資料能夠提供縱、橫波速度和密度等信息,利用這些信息可以進(jìn)一步估算孔隙度、飽和度和泥質(zhì)含量等儲層物性參數(shù)。地震勘探中通常使用Gassmann方程來定量描述低頻條件下儲層物性參數(shù)發(fā)生變化時介質(zhì)有效彈性模量的相應(yīng)變化。采用基于Gassmann方程的流體替換技術(shù)可以得到不同孔隙流體條件下砂巖儲層的體積模量和剪切模量[22]:

(1)

(2)

式中:Ksat和μsat分別是流體飽和巖石的有效體積模量和有效剪切模量;Kdry和μdry分別是干巖石的有效體積模量和有效剪切模量,可以采用改進(jìn)的Hashin-Strikman彈性模量上邊界計算得到[23];K0是巖石基質(zhì)的有效體積模量;Kfl是孔隙流體的有效體積模量;φ是孔隙度。

對于任意巖石,給定各成分體積含量之后,其等效模量將處于上、下界限之間,但其精確的數(shù)值將依賴于幾何細(xì)節(jié)。由于Gassmann-Biot理論假設(shè)巖石具有一個均勻的礦物模量,因此混合礦物常常用一種“平均礦物”來描述。本文采用最簡單、但不一定是最好的界限。以砂、泥巖為例,如果巖石的礦物成分已知,那么利用Voigt-Reuss-Hill平均可以計算巖石基質(zhì)的等效體積模量K0[24]:

(3)

式中:C是泥質(zhì)含量,Ks和Kc分別是砂巖和泥巖的體積模量。

同理,當(dāng)巖石孔隙所含流體為多相流體時,如果流體的成分已知,那么利用Reuss模型[25]可以計算孔隙流體的等效體積模量Kfl:

(4)

式中:Sw是含水飽和度;Khc是油或氣的體積模量;Kbrine是地層水的體積模量。

建立巖石物理模型的另一個參數(shù)是飽和巖石的密度ρsat,該參數(shù)由巖石骨架密度ρ0和充填流體的密度ρfl按體積比例加權(quán)平均得到:

(5)

式中:ρsat,ρ0,ρfl分別為流體飽和巖石、巖石骨架和充填流體的密度。

在計算得到流體飽和巖石的體積模量Ksat,剪切模量μsat和密度ρsat之后,流體飽和巖石的縱波、橫波速度公式可以表示為:

(6)

(7)

1.2 反演目標(biāo)函數(shù)

文獻(xiàn)[11]通過假設(shè)地震噪聲信號和巖石物理模型誤差服從高斯分布,給出了基于最小二乘法的儲層彈性和物性參數(shù)同步反演目標(biāo)函數(shù):

(8)

通過研究(x,y)為何值時使S(x,y)達(dá)到最小可以得到最佳解。然而,公式(8)中的協(xié)方差矩陣Cd和Cx|y對脈沖噪聲十分敏感,容易導(dǎo)致反演算法的不穩(wěn)定。本文假設(shè)地震噪聲信號和巖石物理模型誤差服從高斯分布且?guī)в屑有悦}沖噪聲,利用Hampel三截尾函數(shù)構(gòu)建重加權(quán)算子對協(xié)方差矩陣進(jìn)行修正。新的反演目標(biāo)函數(shù)表示為:

(9)

式中:ωd和ωx|y是基于Hampel三截尾函數(shù)的重加權(quán)算子[26],可表示為:

(10)

(11)

(12)

式中:ω(e)=φ(e)/e,φ(e)=?ρ(e)/?e,ρ(e)是Hampel三截尾函數(shù);ε,Δ1和Δ2是Hampel三截尾函數(shù)的閾值,用來控制函數(shù)對脈沖噪聲的抑制能力,閾值越小,其抗脈沖噪聲的能力越強(qiáng)。圖1為Hampel函數(shù)和加權(quán)算子示意圖?？梢钥闯?當(dāng)|ei|<ε時,Hampel函數(shù)為二次函數(shù);當(dāng)ε≤|ei|<Δ2時,Hampel函數(shù)為一次函數(shù);當(dāng)誤差|ei|>Δ2時,Hampel函數(shù)為常數(shù)。由此可見,Hampel函數(shù)結(jié)合了Huber函數(shù)和Biweight函數(shù)的優(yōu)點(diǎn):當(dāng)誤差較小時,使用L2范數(shù)求解;當(dāng)存在少數(shù)脈沖噪聲時,使用L1范數(shù)求解;當(dāng)存在大量脈沖噪聲時,直接對噪聲信號進(jìn)行截斷,從根本上排除了脈沖噪聲對反演的影響。

顯然,閾值的合理設(shè)置是壓制高斯和脈沖噪聲信號的關(guān)鍵。在噪聲信號服從高斯分布的假設(shè)下,如果噪聲信號大于某一給定閾值T,則其概率可表示為:

(13)

式中,σ表示高斯噪聲信號的標(biāo)準(zhǔn)差。相應(yīng)地,p{|e|>ε},p{|e|>Δ1}和p{|e|>Δ2}分別表示|e|大于ε,Δ1和Δ2的概率。ZOU等[27]提出根據(jù)數(shù)據(jù)的95.0%,97.5%和99.0%置信區(qū)間來設(shè)置ε,Δ1和Δ2,計算獲得相應(yīng)的閾值:

(14)

(15)

(16)

式中:σ(k)是迭代殘差的瞬時標(biāo)準(zhǔn)差。一種常用的方法是利用遞歸公式對瞬時標(biāo)準(zhǔn)差估計進(jìn)行濾波:

(17)

由于σ(k)容易受到脈沖噪聲的影響,這里采用一種魯棒性的估計方法:

(18)

式中:mad{·}表示平均絕對中值偏差,用來抑制脈沖噪聲對方差估計的影響。

圖1 Hampel函數(shù)(a)和由Hampel函數(shù)構(gòu)建的加權(quán)函數(shù)(b)

1.3 牛頓法迭代遞推公式

將公式(9)的梯度用泰勒公式展開,可以近似地得到:

(19)

(20)

式中:

(21)

(22)

(23)

2 效果分析

2.1 模型測試一

首先由測井曲線解釋得到地層巖石的孔隙度、含水飽和度和泥質(zhì)含量等儲層物性參數(shù)(圖2a至圖2c),然后根據(jù)巖石物理模型計算縱、橫波速度和密度等彈性參數(shù)(圖2d至圖2f),最后利用疊前正演模型合成地震角道集(圖3),測試其含脈沖噪聲時本文方法的反演效果。地震子波采用30Hz的零相位雷克子波,道集的最小炮檢距為150m,最大炮檢距為4000m,炮檢距間隔為350m。在反演過程中,初始模型和先驗(yàn)?zāi)Ｐ涂梢酝ㄟ^測井曲線濾波得到,物性參數(shù)的初始模型和先驗(yàn)?zāi)Ｐ拖嗤?/p>

圖2 利用測井曲線解釋得到的儲層參數(shù)(一)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

圖3 利用測井?dāng)?shù)據(jù)合成的地震角道集(一)

圖4 加入30%脈沖噪聲信號的地震角道集(一)

當(dāng)?shù)卣鸾堑兰?0%的隨機(jī)脈沖噪聲時(圖4),分別采用最小二乘法和本文方法進(jìn)行反演,結(jié)果如圖5和圖6所示。由圖5和圖6可以看出,兩

者均能識別反射層位,但最小二乘法的反演結(jié)果不穩(wěn)定,鋸齒現(xiàn)象明顯。這是因?yàn)槭褂米钚《朔〞箽埐畹玫较嗤臋?quán)重,脈沖噪聲的存在使地震數(shù)據(jù)得到偏大的權(quán)重(Cd偏小),巖石物理約束和先驗(yàn)信息的權(quán)重相對減小,彈性參數(shù)和物性參數(shù)的反演結(jié)果主要來自于地震數(shù)據(jù)。本文利用Hampel三截尾函數(shù)設(shè)置加權(quán)矩陣,自適應(yīng)調(diào)整Cd的大小,減小了地震脈沖噪聲在目標(biāo)函數(shù)中的權(quán)重,因而提高了算法對脈沖噪聲的抑制能力。

圖5 利用最小二乘法反演得到的儲層參數(shù)(一)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

圖6 利用本文方法反演得到的儲層參數(shù)(一)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

2.2 模型測試二

首先由測井曲線解釋得到巖石的孔隙度、含水飽和度和泥質(zhì)含量等儲層物性參數(shù)(圖7a至圖7c);然后根據(jù)巖石物理模型計算縱、橫波速度和密度等彈性參數(shù)(圖7d至圖7f,紅色曲線),并在彈性參數(shù)中加入脈沖噪聲(圖7d至圖7f,黑色曲線);最后利用含脈沖噪聲的彈性參數(shù)合成地震角道集(圖8),并在地震角道集中加入30%的隨機(jī)脈沖噪聲(圖9),分別測試巖石物理模型含脈沖噪聲與巖石物理模型和地震角道集都含脈沖噪聲兩種情況下本文方法的反演效果。地震子波采用30Hz的零相位雷克子波,道集的最小炮檢距為150m,最大炮檢距為4000m,炮檢距間隔為350m。在反演過程中,初始模型和先驗(yàn)?zāi)Ｐ涂梢酝ㄟ^測井曲線濾波得到,物性參數(shù)的初始模型和先驗(yàn)?zāi)Ｐ拖嗤?/p>

當(dāng)巖石物理模型含脈沖噪聲、地震角道集不含脈沖噪聲時(圖8),分別采用最小二乘法和本文方法進(jìn)行反演,結(jié)果如圖10和圖11所示。由圖10和圖11可以看出,兩者均能識別反射層位,但最小二乘法的彈性參數(shù)反演結(jié)果與真實(shí)值相差較大。這是因?yàn)槊}沖噪聲的存在使巖石物理模型得到偏大的權(quán)重(Cx|y偏小),地震角道集和先驗(yàn)信息的權(quán)重相對減小,彈性參數(shù)的反演結(jié)果主要來自于巖石物理約束。本文利用Hampel三截尾函數(shù)設(shè)置加權(quán)矩陣,自適應(yīng)調(diào)整Cx|y的大小,減小了巖石物理脈沖噪聲在目標(biāo)函數(shù)中的權(quán)重,提高了算法對脈沖噪聲的抑制能力。

圖7 利用測井曲線解釋得到的儲層參數(shù)(二)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

當(dāng)巖石物理模型和地震角道集都含脈沖噪聲時(圖9),分別采用最小二乘法和本文方法進(jìn)行反演,結(jié)果如圖12和圖13所示。由圖12和圖13可以看出,兩者均能識別反射層位,但使用最小二乘法得到的彈性參數(shù)反演結(jié)果不穩(wěn)定且偏離真實(shí)值。這是因?yàn)樽钚《朔〞箽埐顢?shù)據(jù)得到相同的權(quán)重,增加了脈沖噪聲在目標(biāo)函數(shù)中的權(quán)重。同時,地震數(shù)據(jù)和巖石物理約束條件以及先驗(yàn)信息之間的相對關(guān)系變得不穩(wěn)定,反演結(jié)果既有可能來自于地震數(shù)據(jù),也有可能來自于巖石物理約束。本文方法能夠自適應(yīng)調(diào)整地震數(shù)據(jù)、巖石物理約束和先驗(yàn)信息之間的相對關(guān)系,提高算法的穩(wěn)定性,因此彈性參數(shù)和物性參數(shù)的反演結(jié)果都與真實(shí)值吻合較好。

圖8 利用測井?dāng)?shù)據(jù)合成的地震角道集(二)

圖9 加入30%隨機(jī)脈沖噪聲后的地震角道集(二)

圖10 利用最小二乘法反演得到的儲層參數(shù)(二)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

圖11 利用本文方法反演得到的儲層參數(shù)(二)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

圖12 利用最小二乘法反演得到的儲層參數(shù)(二)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

圖13 利用本文方法反演得到的儲層參數(shù)(二)a 孔隙度; b 含水飽和度; c 泥質(zhì)含量; d 縱波速度; e 橫波速度; f 密度

2.3 應(yīng)用實(shí)例

利用蘇里格氣田實(shí)際數(shù)據(jù)對本文方法的應(yīng)用效果進(jìn)行了研究。研究區(qū)目的層是上古生界二疊系底部H8—S1砂層組,屬沼澤背景下的辮狀河沉積,天然氣成藏條件復(fù)雜。同時,構(gòu)造對天然氣聚集不起主要控制作用,儲層經(jīng)過強(qiáng)烈的成巖作用改造,有效孔隙以次生孔隙為主,具有低孔、低滲和非均質(zhì)性嚴(yán)重的特點(diǎn)。圖14為蘇里格氣田A測線地震剖面,該測線共有601道,道間距為25m,時間采樣率為2ms。該剖面經(jīng)過6口井,各井所在CDP號分別為1665,1691,1792,1856,2007和2044,圖中黑色曲線為H8層位。

圖14 蘇里格氣田疊后地震剖面

圖15是A測線地震角道集的方差分布圖,可以看出,沿著該測線方向的方差變化較大。在這種情況下,地震數(shù)據(jù)、巖石物理約束和先驗(yàn)信息之間的相對權(quán)重關(guān)系十分復(fù)雜,基于噪聲信號高斯分布的假設(shè)往往會導(dǎo)致反演結(jié)果不穩(wěn)定。圖16a是位于CDP1856處的實(shí)際地震角道集,可以看出,1970ms處角道集的振幅遠(yuǎn)遠(yuǎn)大于上覆地層的振幅。圖16b是利用最小二乘法反演結(jié)果生成的地震角道集,圖16c是實(shí)際角道集與合成角道集的殘差分布圖,可以看出,當(dāng)上覆地層的地震殘差很小時,1970ms處的地震殘差仍然很大。顯然,圖16c中的角道集殘差已經(jīng)不再滿足噪聲信號服從高斯分布的假設(shè),假設(shè)地震殘差服從高斯分布且?guī)в屑有悦}沖噪聲更為合理。

圖15 地震角道集的方差分布

圖17和圖18分別是采用最小二乘法和本文方法進(jìn)行反演得到的結(jié)果。可以看出,當(dāng)巖石物理模型和地震角道集都含有脈沖噪聲時,最小二乘法的反演結(jié)果并不理想。一方面,巖石物理模型的脈沖噪聲使反演得到的彈性參數(shù)偏離初始模型的趨勢;另一方面,地震角道集的脈沖噪聲使反演得到的物性參數(shù)變得不可靠。本文方法通過設(shè)置加權(quán)矩陣自適應(yīng)調(diào)整地震角道集、巖石物理約束和先驗(yàn)信息之間的相對權(quán)重,在一定程度上能夠壓制脈沖噪聲的影響,提高算法的穩(wěn)定性。

圖16 CDP1856處地震角道集a 實(shí)際角道集; b 合成角道集; c 實(shí)際角道集與合成角道集的殘差

沿著H8層位曲線,彈性參數(shù)與物性參數(shù)的初始模型和先驗(yàn)?zāi)Ｐ涂梢酝ㄟ^6口井的內(nèi)插、外推模擬得到。圖19是利用本文方法反演得到的縱波速度、橫波速度剖面,圖20是利用本文方法反演得到的密度、孔隙度、含水飽和度和泥質(zhì)含量剖面?？梢钥闯?彈性參數(shù)和物性參數(shù)的反演結(jié)果穩(wěn)定性較好,且反演結(jié)果與測井曲線較為吻合。

圖17 CDP1856處最小二乘法反演結(jié)果a 縱波速度; b 橫波速度; c 密度; d 孔隙度; e 含水飽和度; f 泥質(zhì)含量

圖18 CDP1856處本文方法反演結(jié)果a 縱波速度; b 橫波速度; c 密度; d 孔隙度; e 含水飽和度; f 泥質(zhì)含量

圖19 縱波速度(a)和橫波速度(b)反演結(jié)果

圖20 密度(a)、孔隙度(b)、含水飽和度(c)和泥質(zhì)含量(d)反演結(jié)果

3 結(jié)束語

儲層彈性和物性參數(shù)同步反演是提取巖性參數(shù)的一條重要途徑。隨著勘探目標(biāo)由構(gòu)造油氣藏轉(zhuǎn)向巖性油氣藏,同步反演的應(yīng)用越來越廣泛。本文假設(shè)地震噪聲信號和巖石物理模型誤差服從高斯分布且?guī)в屑有悦}沖噪聲,利用Hampel三截尾函數(shù)構(gòu)建加權(quán)矩陣對反演目標(biāo)函數(shù)進(jìn)行約束,通過反演殘差自適應(yīng)調(diào)整目標(biāo)函數(shù)中的加權(quán)矩陣,提高了同步反演算法對高斯噪聲和脈沖噪聲的適應(yīng)能力。

當(dāng)觀測數(shù)據(jù)不完全服從高斯分布時,多源信息約束反演本身是嚴(yán)重不穩(wěn)定的。利用Hampel三截尾函數(shù)構(gòu)建加權(quán)矩陣能在一定程度上自適應(yīng)調(diào)節(jié)各個目標(biāo)函數(shù)之間的權(quán)重,消除地震和巖石物理噪聲的統(tǒng)計性誤差造成的反演不穩(wěn)定性,最終得到穩(wěn)定的儲層彈性和物性參數(shù)反演結(jié)果。

[1] OSTRANDER W J.Plane-wave reflection coefficients for gas sands at nonnormal angles of incidence[J].Geophysics,1984,49(10):1637-1648

[2] 胡華鋒,印興耀,吳國忱.基于貝葉斯分類的儲層物性參數(shù)聯(lián)合反演方法[J].石油物探,2012,51(3):225-232 HU H F,YIN X Y,WU G C.Joint inversion of petrophysical parameters based on Bayesian classification[J].Geophysical Prospecting for Petroleum,2012,51(3):225-232

[3] 印興耀,孫瑞瑩,張廣智,等.基于分形高頻初始模型和低頻先驗(yàn)信息的物性參數(shù)隨機(jī)反演[J].石油物探,2014,53(5):537-544 YIN X Y,SUN R Y,ZHANG G Z,et al.Stochastic inversion of reservoir physical property parameters based on high-frequency initial model from fractal and low-frequency prior information[J].Geophysical Prospecting for Petroleum,2014,53(5):537-544

[4] 潘昱潔,李大衛(wèi),楊鍇.確定性反演和隨機(jī)反演對井約束條件的需求分析[J].石油物探,2011,50(4):345-349 PAN Y J,LI D W,YANG K.A comparison between the requirements of multi-well constrained conditions in deterministic inversion and stochastic inversion[J].Geophysical Prospecting for Petroleum,2011,50(4):345-349

[5] 張繁昌,代榮獲,劉漢卿,等.疊前地震道集AVA多頻信息同時反演[J].石油物探,2014,53(4):453-460 ZHANG F C,DAI R H,LIU H Q,et al.Multi-frequency AVA simultaneous inversion for pre-stack seismic gathers[J].Geophysical Prospecting for Petroleum,2014,53(4):453-460

[6] 桂金詠,高建虎,雍學(xué)善,等.致密儲層敏感彈性參數(shù)疊前同步反演方法[J].石油物探,2015,54(5):541-550 GUI J Y,GAO J H,YONG X S,et al.A prestack simultaneous inversion method for sensitive elastic parameters of tight reservoir[J].Geophysical Prospecting for Petroleum,2015,54(5):541-550

[7] 王官超,杜啟振.基于包絡(luò)目標(biāo)函數(shù)的彈性波波形反演[J].石油物探,2016,55(1):133-141 WANG G C,DU Q Z.Elastic full waveform inversion based on envelop objective function[J].Geophysical Prospecting for Petroleum,2016,55(1):133-141

[8] BOSCH M,CARA L,RODRIGUES J,et al.A Monte Carlo approach to the joint estimation of reservoir and elastic parameters from seismic amplitudes[J].Geophysics,2007,72(6):O29-O39

[9] BOSCH M.The optimization approach to lithological tomography:combining seismic data and petrophysics for porosity predictio[J].Geophysics,2004,69(5):1272-1282

[10] BOSCH M,CARVAJAL C,RODRIGUES J,et al.Petrophysical seismic inversion conditioned to well-log data:methods and application to a gas reservoir[J].Geophysics,2009,74(2):O1-O15

[11] 李志勇,錢峰,胡光岷,等.儲層彈性與物性參數(shù)地震疊前同步反演的確定性優(yōu)化方法[J].地球物理學(xué)報,2015,58(5):1706-1716 LI Z Y,QIAN F,HU G M,et al.Prestack seismic joint inversion of reservoir elastic and petrophysical parameters using deterministic optimization method[J].Chinese Journal of Geophysics,2015,58(5):1706-1716

[12] LI Z,LIU Z,SONG C,et al.Generalized Gaussian distribution based adaptive mixed-norm inversion for non-Gaussian noise[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:3926-3930

[13] LI Z Y,SONG B B,ZHANG J S,et al.Joint elastic and petrophysical inversion using prestack seismic and well log data[J].Exploration Geophysics,2015,47(4):331-340

[14] TAYLOR H L,BANKS S C,MCCOY J F.Deconvolution with the1 norm[J].Geophysics,1979,44(1):39-52

[15] BUBE K P,LANGAN R T.Hybrid1/2 minimization with applications to tomography[J].Geophysics,1997,62(4):1183-1195

[16] GUITTON A,SYMES W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1319

[17] 呂曉春,顧漢明,成景旺.基于 Huber函數(shù)的頻率域全波形反演[J].石油物探,2013,52(5):544-552 LV X C,GU H M,CHEN J W.Huber norm based full waveform inversion in frequency domain[J].Geophysical Prospecting for Petroleum,2013,52(5):544-552

[18] JI J.Robust inversion using biweight norm and its application to seismic inversion[J].Exploration Geophysics,2012,43(2):70-76

[19] 劉洋,張家樹,胡光岷,等.疊前三參數(shù)非高斯反演方法研究[J].地球物理學(xué)報,2012,55(1):269-276 LIU Y,ZHANG J S,HU G M,et al.A study of the three-term non-Gaussian pre-stack inversion method[J].Chinese Journal of Geophysics,2012,55(1):55-64

[20] 岳碧波,彭真明,張啟衡.基于α穩(wěn)定分布的地震反演方法[J].地球物理學(xué)報,2012,55(4):1307-1317 YUE B B,PENG Z M,ZHANG Q H.Seismic inversion method based onαdistribution[J].Chinese Journal of Geophysics,2012,55(4):1307-1317

[21] ZHANG J,LV S,LIU Y,et al.AVO inversion based on generalized extreme value distribution with adaptive parameter estimation[J].Journal of Applied Geophysics,2013,98(1):11-20

[22] MAVKO G,MUKERJI T,DVORKIN J.The rock physics handbook:tools for seismic analysis of porous media[M].Cambridge:Cambridge University Press,2009:168-169

[23] GRANA D,DELLA R E.Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion[J].Geophysics,2010,75(3):O21-O37

[24] Hill R.The elastic behaviour of a crystalline aggregate[J].Proceedings of the Physical Society(Section A),1952,65(5):349-354

[25] REUSS A.Berechnung der flie?grenze von mischkristallen auf grund der plastizit?tsbedingung für einkristalle[J].Journal of Applied Mathematics and Mechanics,1929,9(1):49-58

[26] HAMPEL F R,RONCHETTI E M,ROUSSEEUW P J,et al.Robust statistics:the approach based on influence functions[M].New York:Wiley,1986:379-380

[27] ZOU Y,CHAN S C,NG T S.A robust m-estimate adaptive filter for impulse noise suppression[J].Proceedings of IEEE International Conference on Acoustics,Speech,and Signal Processing,1999:1765-1768

(編輯：戴春秋)

Simultaneous inversion on elastic and physical properties of reservoirbased on Hampel’s three-part redescending function

LI Zhiyong1,ZHANG Jiashu2,CAI Hanpeng1,HU Guangmin1

(1.SchoolofCommunicationandInformationEngineering,UniversityofElectronicScienceandTechnologyofChina,Chengdu611731,China;2.SichuanKeyLabofSignalInformationProcessing,SouthwestJiaotongUniversity,Chengdu610031,China)

The discrepancies between geophysical measurements and forward modeling data are commonly modeled as Gaussian errors by conventional seismic inversion methods.However,misfit functions constructed by this distribution always have good inversion result in specified regions but have weak adaptability.We assume Gaussian distribution with additional impulse distribution for both seismic noise data and elastic property deviations from the rock-physics model.In this case,Hampel’s three-part redescending estimate function is included in the misfit function,which could perform better to suppress both Gaussian and impulse errors as a robust measure.Minimizing the misfit function classically leads to the Newton’s optimization method,with which elastic and physical properties of reservoir are obtained.This approach can be applied to the areas with different distribution of errors by choosing the degree of threshold.The relative contributions of the rock-physics constraints and prior information are adaptively adjusted,which can increase the stability of the algorithm and ensure a better inversion.Both synthetic and real seismic data examples show that the method can obtain stable and reliable inversion results.

AVO inversion,rock-physics inversion,joint inversion,Gaussian noise,impulse noise,Hampel three-part redescending function

2016-03-03;改回日期：2016-04-26。

李志勇(1985—),男,博士在讀,主要從事地震儲層反演方面的研究。

胡光岷(1966—),男,教授,現(xiàn)主要從事網(wǎng)絡(luò)行為學(xué)與網(wǎng)絡(luò)安全、網(wǎng)絡(luò)體系結(jié)構(gòu)與協(xié)議、信號與信息處理技術(shù)的應(yīng)用等方面的研究與教學(xué)工作。

國家自然科學(xué)基金聯(lián)合基金項(xiàng)目(U1562218)資助。

P631

A

1000-1441(2017)02-0261-12

10.3969/j.issn.1000-1441.2017.02.013

This research is financially supported by the Joint Funds of the National Natural Science Foundation of China (Grant No.U1562218).

李志勇,張家樹,蔡涵鵬,等.基于Hampel三截尾函數(shù)的儲層彈性和物性參數(shù)同步反演[J].石油物探,2017,56(2):-272

LI Zhiyong,ZHANG Jiashu,CAI Hanpeng,et al.Simultaneous inversion on elastic and physical properties of reservoir based on Hampel’s three-part redescending function[J].Geophysical Prospecting for Petroleum,2017,56(2):-272