舒夢珵， 王彥飛

1 中海油研究總院， 北京 100028 2 中國科學(xué)院地質(zhì)與地球物理研究所油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029

?

儲層重力密度反演后驗(yàn)約束正則化方法

舒夢珵1， 王彥飛2*

1 中海油研究總院， 北京 100028 2 中國科學(xué)院地質(zhì)與地球物理研究所油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029

本文針對蒸汽輔助重力泄油(SAGD)生產(chǎn)中開發(fā)監(jiān)測問題，發(fā)展了綜合應(yīng)用地震及重力數(shù)據(jù)反演儲層密度的聯(lián)合反演算法.通過測井?dāng)?shù)據(jù)建立縱波阻抗與密度的直接關(guān)系，并推導(dǎo)出這種關(guān)系下重力與縱波阻抗數(shù)據(jù)聯(lián)合反演的計(jì)算方法，從而計(jì)算出蒸汽腔體密度分布規(guī)律.文中應(yīng)用密度反演后驗(yàn)約束正則化方法，采用Tikhonov正則化模型，通過波阻抗數(shù)據(jù)作為約束進(jìn)行聯(lián)合反演，在算法上提高了穩(wěn)定性，同時得到較高的反演精度.文中對SAGD生產(chǎn)中的理論模型進(jìn)行了方法試算，并分析了算法的誤差，最終應(yīng)用于SAGD生產(chǎn)的實(shí)際數(shù)據(jù)中，通過最終反演結(jié)果分析，該方法取得了很好的應(yīng)用效果.

SAGD生產(chǎn)； 重力； 縱波阻抗； Tikhonov正則化； 聯(lián)合反演

1 引言

地球物理數(shù)據(jù)的聯(lián)合反演方法發(fā)源于20世紀(jì)70年代中期，Vozoff和Jupp(1977)開創(chuàng)了聯(lián)合反演的先河，用迭代二階馬奎特阻尼最小二乘法實(shí)現(xiàn)了一維直流電測深(DC)和大地電磁測深(MT)資料的聯(lián)合反演，解決了層狀介質(zhì)中的各向異性問題.Savino等(1980)利用地震P波走時和重力資料聯(lián)合反演，研究華盛頓東部地區(qū)地殼上地幔密度、速度結(jié)構(gòu).Lines等(1988)使用地面地震數(shù)據(jù)、聲波測井，地面重力及井中重力等資料研究了地震、重力同步反演及順序反演方法，在反演過程中，充分利用了已有的地面地震、井下聲波測井、VSP 數(shù)據(jù)、地面及井下重力數(shù)據(jù)等資料，從而大大縮小了模型的選擇范圍，減小了反問題的多解性，強(qiáng)化了解釋過程，并分析了同步反演及順序反演的各自的應(yīng)用效果.Dobróka等(1991)對垂直地震剖面(VSP)走時數(shù)據(jù)、電法數(shù)據(jù)，采用基于最大頻率值(MFV)的加權(quán)最小二乘算法進(jìn)行聯(lián)合反演，與阻尼最小二乘算法相比，該算法具有估計(jì)誤差小以及初始模型選擇對結(jié)果影響較小兩個特點(diǎn).Sun和Schuster (1992)提出了一個在層析成像反演中多個目標(biāo)函數(shù)的極小化過程，該過程在層析成像反演中是十分有用的，特別是同時做幾種類型數(shù)據(jù)模擬，該過程將分級的優(yōu)化問題轉(zhuǎn)成為等效約束優(yōu)化的問題，從而使問題簡單化.Zeyen和Pous(1993)在具有先驗(yàn)信息的基礎(chǔ)上，如密度、磁化率、剩余磁化強(qiáng)度等，對重、磁場的聯(lián)合反演問題進(jìn)行了研究.而張貴賓等(1993)以BG理論為基礎(chǔ)，在重磁異常線性反演中將該理論與吉洪諾夫正則化方法相結(jié)合求解地下密度源(或磁源)分布及質(zhì)心(或磁質(zhì)心)位置;在重、磁非線性反演中結(jié)合應(yīng)用正則化方法和馬奎特思想給出一種確定地下密度(或磁性)界面的穩(wěn)定迭代算法——正則馬奎特法.楊輝(1998)以地震資料解釋的三維構(gòu)造圖作為先驗(yàn)信息，用重力三維正演剝離基底及基底以上界面所產(chǎn)生的重力效應(yīng)，然后對分離后的基底巖性異常用穩(wěn)健的SVD 算法來線性反演基底密度差.最后，利用重、磁、電、震、地面地質(zhì)、鉆井等資料綜合解釋了盆地的基底時代及巖性，取得了令人滿意的地質(zhì)效果.21世紀(jì)初，Parsons等(2001)應(yīng)用重力數(shù)據(jù)及地震數(shù)據(jù)開展了順序反演的工作，發(fā)展了重力數(shù)據(jù)及地震走時層析成像數(shù)據(jù)的順序綜合反演算法.Tondi等(2000,2003)綜合利用折射波、寬角度反射波走時以及重力數(shù)據(jù)進(jìn)行順序綜合反演來構(gòu)造二維速度模型.隨后又提出地震及重力數(shù)據(jù)順序聯(lián)合反演固體地球深-淺部分構(gòu)造重建的方法，取得了很好的應(yīng)用效果.Onizawa等(2002)給出了地震數(shù)據(jù)與重力數(shù)據(jù)同步反演較好的算法實(shí)例，Coutant等(2012)在其工作基礎(chǔ)上進(jìn)行改進(jìn)，并在La Soufriere地區(qū)熱液系統(tǒng)勘查中進(jìn)行了應(yīng)用.Hayashi等(2005)應(yīng)用面波及微重力觀測數(shù)據(jù)，通過建立橫波速度與密度之間的關(guān)系進(jìn)行聯(lián)合反演，并應(yīng)用阻尼最小二乘的方法進(jìn)行求解，取得了很好的效果.Moorkamp等(2011)給出MT、重力及散射波數(shù)據(jù)聯(lián)合反演的全局方法.

前人諸多的研究成果表明(Maceira and Ammon,2009;Gallardo and Meju,2004;Tikhotsky and Achauer,2008; Basuyau and Tiberi,2011)，聯(lián)合反演在算法上可以大致分為兩類：(1)順序反演和(2)同步反演.在數(shù)據(jù)應(yīng)用上可以分為兩類：(1)基于相同物性地球物理觀測數(shù)據(jù)之間的聯(lián)合反演.如反射地震的走時和振幅聯(lián)合反演，縱波和橫波資料的聯(lián)合反演;不同電法或電磁法所取得的觀測數(shù)據(jù)的聯(lián)合反演.這種聯(lián)合反演有天然的合理性，因?yàn)樗鼈兙谙嗤膸r石物性差異，物理基礎(chǔ)相同，觀測場之間必然存在著相關(guān)性；(2)基于不同巖石物性的地球物理觀測數(shù)據(jù)之間的聯(lián)合反演，如地震和重力，地震和MT，重力和MT 等.這種聯(lián)合的基礎(chǔ)是不同物性之間存在著相關(guān)的內(nèi)在聯(lián)系，由相關(guān)的物性參數(shù)必然會誘發(fā)相關(guān)的物理異常，這是利用多種地球物理信息進(jìn)行聯(lián)合反演的地球物理基礎(chǔ)(楊輝等，2002).

本文針對稠油開采中SAGD生產(chǎn)實(shí)際問題，分析單井與水平井組合情況下蒸汽腔體分布規(guī)律(Bulter,2004).基于密度反演后驗(yàn)約束正則化方法(Wang and Xiao,2001; 王彥飛，2007;王彥飛等，2011)，求解地震數(shù)據(jù)與重力數(shù)據(jù)聯(lián)合反演的問題.通過測井?dāng)?shù)據(jù)，建立縱波阻抗與密度的直接關(guān)系，進(jìn)行聯(lián)合反演，得到儲層含汽砂巖及含稠油砂巖的密度分布規(guī)律，最終分析SAGD生產(chǎn)中蒸汽腔體分布情況.

2 重力正演模型

我們首先給出密度模型和觀測系統(tǒng)：采用地表觀測和接收數(shù)據(jù)，g(xi)為在第i個觀測點(diǎn)的重力接收數(shù)據(jù)，dj為地下介質(zhì)密度模型網(wǎng)格化后第j個網(wǎng)格的密度值(通常是未知的)，如圖1所示.

圖1 密度模型及觀測系統(tǒng)示意圖Fig.1 Map of density model and acquisition geometry

根據(jù)圖1給出的重力觀測模式，可以得到重力計(jì)算的物理模型如下：

(1)

其中，fk，j為給定重力場觀測網(wǎng)格后的距離系數(shù)，定義作

(2)

其中，dx,dz分別為地下測量區(qū)域網(wǎng)格的長度和寬度；(xk,zk)為地下觀測點(diǎn)坐標(biāo)，γ為牛頓萬有引力常數(shù)(γ=6.67×10-11N·m2/kg2).

上式可簡寫為如下形式：

(3)

其中，F(xiàn)為位置距離算子，d為地下每個網(wǎng)格單元的密度值，g為地面觀測到的重力異常值.

3 反演模型正則化

在實(shí)際問題中，由于測量誤差的存在，觀測到的重力異常數(shù)據(jù)通常是帶噪音的.假設(shè)g為重力異常理論值，實(shí)際數(shù)據(jù)記為gδ，并且假定

(4)

于是(3)式并不能精確求得,一般需要求解一個最小二乘問題：

(5)

然而直接求解由(5)式導(dǎo)出的正規(guī)化方程組是嚴(yán)重的病態(tài)問題，這是因?yàn)镕*F的條件數(shù)要遠(yuǎn)遠(yuǎn)大于F的條件數(shù).因此我們需要研究合適的正則化技巧.

3.1 Tikhonov正則化模型

由于觀測不足(n≤m)，重力反演問題通常是不適定的，因此直接求解(5)式相關(guān)的最小二乘問題通常會帶來不穩(wěn)定的計(jì)算效果(王彥飛等，2011).為克服問題的不適定性，本文建立如下基于2-范數(shù)的密度反演約束最優(yōu)化模型：

(6)

(F*F+αI)(d-d0)=F*(gδ-Fd0).

(7)

于是密度模型可由公式

d(α)=d0+(F*F+αI)-1F*(gδ-Fd0)

(8)

獲得.現(xiàn)在的關(guān)鍵問題是求解合適的正則參數(shù)α.

3.2 計(jì)算正則參數(shù)與正則解

利用偏差原理，公式(7)中α滿足如下非線性方程(WangandXiao,2001)：

(9)

對φ(α)=0求根，可以通過牛頓求根公式得到，即

(10)

注意到

(11)

于是公式(10)中的φ′(α)可表達(dá)為

(12)

由以上推導(dǎo)，根據(jù)正則化反演原理，本文給出如下可以估計(jì)正則參數(shù)同時可以反演重力密度的正則化反演算法：

(1) 取初始迭代值α0=0.01，δ=0.01，ε=0.0001，kmax=30，其中ε為算法終止迭代誤差，kmax為終止迭代最大步數(shù).給定位置距離算子F，重力異常值gδ，初始先驗(yàn)密度模型d0，并令k∶=0.

(2) 利用高斯消去法解如下所示方程組:

(13)

(3) 計(jì)算φ(αk)及φ′(αk):

(14)

(4) 求取參數(shù)αk+1:

(5) 終止準(zhǔn)則：確定求解是否小于解的合理誤差范圍，若

則終止反演算法，輸出反演結(jié)果d(αk).

4 數(shù)值試驗(yàn)

4.1 簡單模型正反演模擬

可以采用SAGD稠油熱采的簡化模型進(jìn)行試算(楊立強(qiáng)等,2007;孟巍等,2006;劉尚奇等，2007)，如圖2中的模型所示.正、反演模擬的數(shù)據(jù)參數(shù)為:

(1)網(wǎng)格分布：1)水平方向70個；2)垂直方向30個；

(2)網(wǎng)格大小為：10 m×10 m；

(3)觀測間隔為10 m.這樣的網(wǎng)格分布與實(shí)際測量中的觀測系統(tǒng)一致.

4.2 理論模型正演計(jì)算

圖2中的模型為SAGD稠油生產(chǎn)中雙直井與水平井組合方式下蒸汽腔體分布模型.模型給出了水平井垂直于紙面的情況.

模型中異常體(圖2藍(lán)色部分)表示由于注入蒸汽、稠油被驅(qū)替的情況.產(chǎn)生的蒸汽代替稠油的空間位置，并造成儲層區(qū)域重力值的變化.

模型的密度分別按照工區(qū)內(nèi)含稠油礫巖巖石密度及含蒸汽礫巖巖石密度的平均值給出：(1) 含稠油礫巖巖石密度為2.43 g·cm-3(圖2中褐色區(qū)域)；(2) 含蒸汽礫巖巖石密度為1.98 g·cm-3(圖2中藍(lán)色區(qū)域).

由于實(shí)際情況中稠油與蒸汽腔體有過渡帶存在，為了使模型更為符合實(shí)際情況，文中對模型異常體邊界進(jìn)行線性平滑處理，以模擬油-水-汽的過渡變化關(guān)系.

對模型進(jìn)行重力異常的正演計(jì)算，可以得到圖3a的重力異常正演結(jié)果.

本文提出應(yīng)用縱波阻抗數(shù)據(jù)作為約束進(jìn)行重力聯(lián)合反演的方法，模型的縱波阻抗數(shù)據(jù)如圖4所示.

4.3 理論模型反演計(jì)算

下面根據(jù)第2節(jié)給出的重力模型并利用我們提出的正則化方法進(jìn)行反演計(jì)算.位置距離算子矩陣如圖5所示.計(jì)算表明，該位置距離矩陣是極度壞條件的，條件數(shù)達(dá)到3.1140×1021，因而必須應(yīng)用正則化反演技巧.

應(yīng)用圖3a中理論正演重力數(shù)據(jù)進(jìn)行反演，可得到圖6所示密度反演結(jié)果.為了驗(yàn)證反演算法的精度，給出反演結(jié)果與真實(shí)模型的誤差分布如圖7所示.文中在正演的觀測數(shù)據(jù)中加入隨機(jī)噪音，以模擬真實(shí)數(shù)據(jù)情況并說明提出算法的抗噪效果.信噪比計(jì)算公式為r=10 log(S/N)，其中S和N分別是信號和噪聲的平均功率，單位為dB.圖8給出信噪比為90時的重力正演數(shù)據(jù).把初始密度模型(如圖3b)帶入反演模型，應(yīng)用上述正則化反演算法，經(jīng)過迭代計(jì)算最終可以得到反演的模型密度如圖9所示.為了驗(yàn)證反演算法的精度，給出反演結(jié)果與真實(shí)模型的誤差分布如圖10所示.

圖2 SAGD生產(chǎn)二維密度模型:雙直井與水平井組合二維模型(水平井垂直于紙面)(密度單位為g·cm-3)Fig.2 2D forward density model in SAGD production: double straight wells with horizontal well joint 2D model (horizontal well perpendicular to the paper) (density unit: g·cm-3)

圖3 模型正演結(jié)果: (a) 模型重力正演結(jié)果； (b)蒸汽腔體密度分布模型(密度單位為：g·cm-3)Fig.3 Forward simulation result: (a) gravity forward result; (b) density distribution model of steam chambers (density unit: g·cm-3)

圖4 相對波阻抗模型(色標(biāo)中數(shù)字無量綱)Fig.4 Relative P-wave impedance model

圖5 模型的位置距離算子矩陣(縱橫坐標(biāo)分別為距離矩陣的維數(shù)大小，色標(biāo)表示元素的大小)Fig.5 Matrix of relative position operator

圖6 通過理論正演重力數(shù)據(jù)反演密度結(jié)果 (密度單位為g·cm-3)Fig.6 Density inversion result map via gravity data without noise (density unit: g·cm-3)

圖7 反演密度誤差結(jié)果(密度單位為g·cm-3)Fig.7 Error map of density inversion (density unit: g·cm-3)

圖8 信噪比為90的正演重力數(shù)據(jù)Fig.8 Forward gravity data curve with S/N=90

圖9 反演密度結(jié)果(密度單位為g·cm-3)Fig.9 Density inversion result (density unit: g·cm-3)

圖10 反演密度誤差結(jié)果(密度單位為g·cm-3)Fig.10 Error propagation of density inversion (density unit: g·cm-3)

圖11 迭代得到的正則參數(shù)αFig.11 Regularization parameter α

值得指出的是，正則參數(shù)α在密度反演中起著至關(guān)重要的作用.應(yīng)用3.2節(jié)中的算法迭代得到的正則參數(shù)序列如圖 11所示.正則參數(shù)分布曲線說明，隨著迭代的進(jìn)行，問題的條件數(shù)逐漸得到改善，正則解越來越逼近問題的真實(shí)解(王彥飛，2007).

4.4 實(shí)際數(shù)據(jù)試驗(yàn)

針對SAGD開發(fā)區(qū)，觀測重力結(jié)果如圖12所示.利用阻抗Z和密度ρ對該區(qū)內(nèi)兩口觀察井的測井曲線進(jìn)行分析，建立二者之間的二次多項(xiàng)式函數(shù)關(guān)系表達(dá)式如下：

(15)

其中，ρ為密度測井結(jié)果，Z為測井曲線計(jì)算得到的縱波阻抗結(jié)果.利用最小二乘法，可以得到公式(15)中的參數(shù)c0=20719000.0，c1=-20580.7，c2=5.77248.由公式(15)所得阻抗Z和密度ρ的二次回歸曲線如圖13.

利用阻抗Z和密度ρ的關(guān)系式，得到縱波阻抗數(shù)據(jù)如圖14所示.利用縱波阻抗數(shù)據(jù)作為約束，把觀測數(shù)據(jù)帶入反演算法進(jìn)行實(shí)際重力數(shù)據(jù)的反演計(jì)算，密度反演結(jié)果如圖15所示.我們從反演結(jié)果中可以看到，水平井注、采汽位置及垂直井注汽位置附近有低密度分布區(qū)，文中初步判斷為SAGD生產(chǎn)中蒸汽替換部分.同時在低密度區(qū)周圍有相對高密度區(qū)分布，文中推測這些相對高密度區(qū)是由含水砂巖及含稠油砂巖混合相態(tài)介質(zhì)組成.

文中的密度反演結(jié)果直接反映了儲層密度的變化，該結(jié)果可以用來分析含汽砂巖、含稠油砂巖、油汽水混合砂巖的空間分布區(qū)域.這說明我們提出的聯(lián)合反演方法可以對調(diào)整SAGD生產(chǎn)中的注采方案起到輔助作用.

5 結(jié)論

本文針對SAGD生產(chǎn)中開發(fā)監(jiān)測問題，發(fā)展了綜合應(yīng)用地震及重力數(shù)據(jù)反演儲層密度的聯(lián)合反演算法.本文取得了如下研究成果：

(1) 提出應(yīng)用縱波阻抗數(shù)據(jù)作為約束進(jìn)行重力聯(lián)合反演的方法；

圖12 觀測線重力異常分布Fig.12 Distribution curve of observed abnormal gravity data

圖13 測井曲線中密度與波阻抗曲線二次回歸關(guān)系曲線Fig.13 The square regression curve of density and P-wave impedance logs

圖14 縱波阻抗剖面(圖中虛線矩形為直井注汽位置；空心圓圈為注、采汽水平井)Fig.14 P-wave impedance section (dotted rectangle: steam injection range of straight well; hollow circles:produced horizontal well)

圖15 密度反演結(jié)果(圖中虛線矩形為直井注汽位置；空心圓圈為注、采汽水平井)(密度單位為g·cm-3)Fig.15 Density inversion section (dotted rectangle: steam injection range of straight well; hollow circles: produced horizontal well) (density unit: g·cm-3)

(2) 提出密度反演后驗(yàn)約束正則化方法，特別是給出了后驗(yàn)選取正則參數(shù)的計(jì)算方法；

(3) 建立了阻抗和密度的二次回歸模型，針對理論模型和實(shí)際數(shù)據(jù)，通過波阻抗數(shù)據(jù)作為約束進(jìn)行聯(lián)合反演，在算法上提高了穩(wěn)定性，同時得到較高的反演精度.

(4) 通過聯(lián)合反演在實(shí)際數(shù)據(jù)中的應(yīng)用可以看到，反演得到的低密度分布區(qū)域與已知井注汽位置吻合很好，且反演得到的低密度區(qū)連通性特征與正演模型形態(tài)相似，這些特征均驗(yàn)證了聯(lián)合反演算法在實(shí)際生產(chǎn)中的有效性.

同時，文中的聯(lián)合反演方法對信噪比及初始模型的要求較高，在實(shí)際應(yīng)用中需要較高的數(shù)據(jù)信噪比及較為精確的初始模型作為約束.

致謝 十分感謝審稿人提出的寶貴意見，使得論文內(nèi)容更加充實(shí).感謝遼河油田提供相關(guān)數(shù)據(jù)并允許此文的發(fā)表,同時感謝對于此工作給予重大幫助的同事們.

Basuyau C, Tiberi C. 2011. Imaging lithospheric interfaces and 3D structures using receiver functions, gravity, and tomography in a common inversion scheme.Computers&Geosciences, 37(9): 1381-1390.

Bulter R M. 2004. The behavior of non-condensable gas in SAGD: A rationalization.JournalofCanadianPetroleumTechnology, 43(1): 28-34.

Coutant O, Bernard M L, Beauducel F, et al. 2012. Joint inversion of P-wave velocity and density, application to La Soufrière of Guadeloupe hydrothermal system.GeophysicalJournalInternational, 191(2): 723-742.Dobróka M, Gyulai á, Ormos T, et al. 1991. Joint inversion of seismic and geoelectric data recorded in an underground coal mine.Geophys.Prosp., 39(5): 643-665.Gallardo L A, Meju M A. 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints.J.Geophys.Res., 109(B3), doi: 10.1029/2003JB002716.Hayashi K, Matsuoka T, Hatakeyama H. 2005. Joint analysis of a surface-wave method and micro-gravity survey.JournalofEnvironmentalandEngineeringGeophysics, 10(2): 175-184.

Lines L R, Schultz A K, Treitel S. 1988. Cooperative inversion of geophysical data.Geophysics, 53(1): 8-20.

Liu S Q, Wang X C, Gao Y R, et al. 2007. SAGD process with the combination of vertical and horizontal wells in super-heavy oil reservoir.PetroleumExplorationandDevelopment(in Chinese), 34(2): 234-238.Maceira M, Ammon C J. 2009. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure.J.Geophys.Res., 114(B2), doi: 10.1029/2007JB005157.

Meng W, Jia D, Xie J N, et al. 2006. Optimization of geological design for SAGD process in super heavy oil reservoirs using a combination of vertical and horizontal wells.JournalofDaqingPetroleumInstitute(in Chinese), 30(2): 44-47.

Moorkamp M, Heincke B, Jegen M, et al. 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data.Geophys.J.Int., 184(1): 477-493.

Onizawa S, Mikada H, Watanabe H, et al. 2002. A method for simultaneous velocity and density inversion and its application to exploration of subsurface structure beneath Izu-Oshima volcano, Japan.Earth,PlanetsandSpace, 54(8): 803-817.

Parsons T, Blakely R, Brocher T. 2001. A simple algorithm for sequentially incorporating gravity observations in seismic travel time tomography.Int.Geol.Rev., 43(12): 1073-1086.

Savino J M, Rodi W L, Masso J F. 1980. Simultaneous inversion of multiple geophysical data sets for earth structure. 45th Annual International Meeting, SEG.

Sun Y H, Schuster G T. 1992. Hierarchic optimizations for smoothing and cooperative inversion. Annual Meeting Abstracts, Society of Exploration Geophysicists, 745-748.Tikhotsky S, Achauer U. 2008. Inversion of controlled-source seismic tomography and gravity data with the self-adaptive wavelet parametrization of velocities and interfaces.Geophys.J.Int., 172(2): 619-630.

Tondi R, De Franco R, Barzaghi R. 2000. Sequential integrated inversion of refraction and wide-angle reflection travel times and gravity data for two-dimensional velocity structures.Geophys.J.Int., 141(3): 679-698.

Tondi R, De Franco R, Biella G. 2003. Deep and shallow solid-earth structures reconstructed with sequential integrated inversion (SII) of seismic and gravity data.SpaceScienceReviews, 108(1-2): 115-130.

Vozoff K, Jupp D L B. 1977. Effective search for a buried layer: An approach to experimental design in geophysics.SocietyofExplorationGeophysicists, 8(1): 6-15.

Wang Y F. 2007. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press.

Wang Y F, Stepanova I E, Strakhov V N, et al. 2011. Inverse Problems in Geophysics and Solution Methods (in Chinese). Beijing: Higher Education Press.

Wang Y F, Xiao T Y. 2001. Fast realization algorithms for determining regularization parameters in linear inverse problems.InverseProblems, 17(2): 281-291.

Yang H. 1998. Basement density inversion using gravimetric and seismic data and the integrative interpretation.OGP(in Chinese), 33(4): 496-502.

Yang H, Dai S K, Song H B, et al. 2002. Overview of joint inversion of integrated geophysics.ProgressinGeophysics(in Chinese), 17(2): 262-271.

Yang L Q, Chen Y M, Wang H Y, et al. 2007. Physical and numerical simulation of steam assisted gravity drainage with vertical and horizontal well combination in extra heavy oil reservoir.JournalofChinaUniversityofPetroleum(in Chinese), 31(4): 64-69.Zeyen H, Pous J. 1993. 3-D joint inversion of magnetic and gravimetric data with a priori information.GeophysicalJournalInternational, 112(2): 244-256.

Zhang G B, Shen N H, Wang X C, et al. 1993. The program system on generalized linear cooperative inversion of the potential anomalies.JournalofChangchunUniversityofEarthSciences(in Chinese), 23(2): 197-204.

附中文參考文獻(xiàn)

劉尚奇, 王曉春, 高永榮等. 2007. 超稠油油藏直井與水平井組合SAGD技術(shù)研究. 石油勘探與開發(fā), 34(2): 234-238.

孟巍, 賈東, 謝錦男等. 2006. 超稠油油藏中直井與水平井組合SAGD技術(shù)優(yōu)化地質(zhì)設(shè)計(jì). 大慶石油學(xué)院學(xué)報(bào), 30(2): 44-47.

王彥飛. 2007. 反演問題的計(jì)算方法及其應(yīng)用. 北京: 高等教育出版社.

王彥飛, 斯捷潘諾娃I E, 提塔連科V N等. 2011. 地球物理數(shù)值反演問題. 北京: 高等教育出版社.

楊輝. 1998. 重力、地震聯(lián)合反演基巖密度及綜合解釋. 石油地球物理勘探, 33(4): 496-502.

楊輝, 戴世坤, 宋海濱等. 2002. 綜合地球物理聯(lián)合反演綜述. 地球物理學(xué)進(jìn)展, 17(2): 262-271.

楊立強(qiáng), 陳月明, 王宏遠(yuǎn)等. 2007. 超稠油直井-水平井組合蒸汽輔助重力泄油物理和數(shù)值模擬. 中國石油大學(xué)學(xué)報(bào)(自然科學(xué)版), 31(4): 64-69.

張貴賓, 申寧華, 王喜臣等. 1993. 位場廣義線性綜合反演系統(tǒng)的建立. 長春地質(zhì)學(xué)院學(xué)報(bào), 23(2): 197-204.

(本文編輯 胡素芳)

The posterior constrained regularization method for reservoir density inversion

SHU Meng-Cheng1, WANG Yan-Fei2*

1CNOOCResearchInstitute，Beijing100028,China2KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China

We present the result of a seismic and gravity joint reservoir density inversion algorithm applied to monitoring the process of steam assisted gravity drainage (SAGD). The method can be used to analyze the density distribution of steam chambers in SAGD in the case of combination of the single well with horizontal wells.Gravimetric density and seismic impedance are combined as a posteriori constraint to formulate the Tikhonov regularization model. Minimization of the Tikhonov regularization model can yield a reasonable inversion result of the density distribution of steam chambers in SAGD. First, a regression model of the seismic impedance and the density through well logging data is built, which can be used to calculate the density constraint data through the seismic P wave data and gives the density constraint. Then, an optimization model based on L2 norm is established. In solving the minimization model, Euler equations along with proper choice of the regularization parameter are simultaneously solved. Utilizing the micro-gravimetric Bouguer gravity anomaly data, the joint inversion scheme based on seismic impedance data and gravimetric data is performed, which yields the density distribution of the reservoir. By analyzing properties of the reservoir density, the density distribution of steam chambers in SAGD can be drawn.Faced with the problem of development and monitoring during SAGD production, this paper investigates the joint inversion method based on the seismic and gravimetric data. The results obtained are as follows: (1) applying the P wave data as a constraint to invert the gravimetric density is proposed; (2) a posteriori regularization method is utilized to calculate the density, in particular, a posteriori choice of the regularization parameter is performed; (3) a quadratic regression model is built for the seismic impedance and the density, theoretical and field data applications reveal that the joint inversion algorithm is stable and can yield high precision of inversion results; (4) it is shown from the field data applications using our joint inversion method that the inverted low density distribution area coincides well with the known well gas injection position, and the connectivity of the low density area is similar to the forward model, hence proves the feasibility of the joint inversion algorithm in practice. We also remark that the joint inversion method relies on the proper initial model and high SNR data; this may be a requirement for practical usage of the joint inversion method.For the development and monitoring problem of SAGD production, we develop a joint inversion method based on the seismic data and gravimetric data. Using a posteriori choice of the regularization parameter and a quadratic regression model for the seismic impedance and the density, an inversion procedure is performed by solving an Euler equation. The inversion results can well explain and draw the density distribution of steam chambers in SAGD.

SAGD production; Gravity; P-wave impedance; Tikhonov regularization; Joint inversion

10.6038/cjg20150622.

國家杰出青年科學(xué)基金項(xiàng)目(41325016)與國家自然科學(xué)基金面上項(xiàng)目(11271349)資助.

舒夢珵，男，中國科學(xué)院地質(zhì)與地球物理研究所博士，研究方向?yàn)榭碧降厍蛭锢?

*通訊作者 王彥飛，男，研究員，2002年于中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院獲博士學(xué)位，主要從事計(jì)算及勘探地球物理領(lǐng)域的研究工作. E-mail:yfwang@mail.iggcas.ac.cn

10.6038/cjg20150622

P631

2014-06-24，2015-05-23收修定稿

舒夢珵，王彥飛. 2015. 儲層重力密度反演后驗(yàn)約束正則化方法.地球物理學(xué)報(bào)，58(6):2079-2086,

Shu M C, Wang Y F. 2015. The posterior constrained regularization method for reservoir density inversion.ChineseJ.Geophys. (in Chinese),58(6):2079-2086,doi:10.6038/cjg20150622.