馮國(guó)慶 何玉俊* 劉紅林 陳 燕 張 萍④ 薛芳芳

(①西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川成都 610500； ②中國(guó)石油新疆油田分公司采油二廠，新疆克拉瑪依 834008； ③中國(guó)石化河南油田分公司采油二廠，河南南陽(yáng) 473100； ④中國(guó)石化重慶涪陵頁(yè)巖氣勘探開發(fā)有限公司，重慶 408014； ⑤中國(guó)石油長(zhǎng)慶油田分公司油氣工藝研究院，陜西西安 710018)

0 引言

隨機(jī)模擬方法通常用于儲(chǔ)層預(yù)測(cè)[1-6]。影響隨機(jī)建模精度的關(guān)鍵因素是基礎(chǔ)資料的質(zhì)量、數(shù)量及方法的選擇[7]。只有選擇合適的建模方法，利用資料之間的相互約束關(guān)系，才能獲得更為精確的地質(zhì)模型[8]。

試井是油氣勘探開發(fā)過(guò)程中認(rèn)識(shí)儲(chǔ)層的一個(gè)重要手段。試井資料可以反映地層的橫向變化，對(duì)儲(chǔ)層平面展布特征的認(rèn)識(shí)有著重要的作用。但是，目前結(jié)合試井資料的儲(chǔ)層建模方法[9-11]較少，并且大多采用的流程為：先利用壓力恢復(fù)試井?dāng)?shù)據(jù)，反演井周圍的滲透率；再利用滲透率數(shù)據(jù)求出滲透率的變差函數(shù)；然后在變差函數(shù)的約束下模擬生成滲透率模型。該方法存在一定局限性，主要是因?yàn)橛脡毫謴?fù)試井?dāng)?shù)據(jù)求得的滲透率數(shù)據(jù)有限，根據(jù)其計(jì)算的變差函數(shù)難以保證準(zhǔn)確性和連續(xù)性。

本文將試井測(cè)試資料用于儲(chǔ)層建模方法中，首先使用傳統(tǒng)隨機(jī)建模方法結(jié)合靜態(tài)資料建立儲(chǔ)層初始地質(zhì)模型，再利用模擬退火算法結(jié)合試井動(dòng)態(tài)資料對(duì)初始模型進(jìn)行修正，建立同時(shí)符合靜態(tài)資料和試井動(dòng)態(tài)資料的地質(zhì)模型。該方法可保證模型有更好的連續(xù)性，同時(shí)修正孔隙度和滲透率模型，應(yīng)用范圍更廣。

1 方法原理

首先利用靜態(tài)資料和傳統(tǒng)的隨機(jī)地質(zhì)建模方法[12-16](如序貫高斯模擬、序貫指示模擬等)建立儲(chǔ)層的初始屬性模型，包括孔隙度模型和滲透率模型等。初始模型建立后，結(jié)合試井動(dòng)態(tài)資料(主要為試井壓力數(shù)據(jù))，利用模擬退火算法[17-21]對(duì)屬性模型進(jìn)行修正。

模型修正方法為：在初始模型中產(chǎn)生隨機(jī)擾動(dòng)，在擾動(dòng)后的參數(shù)場(chǎng)中，利用數(shù)值模擬計(jì)算單井井底壓力；然后比較模型計(jì)算的壓力與實(shí)際測(cè)試壓力之間的誤差，根據(jù)誤差判斷是否進(jìn)行下一次擾動(dòng)；如果計(jì)算的各單井的壓力與實(shí)際測(cè)試壓力之間誤差滿足要求，則結(jié)束擾動(dòng)，從而得到修正后的模型。

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

模擬退火算法源于固體退火原理。固體退火是指先將固體加熱至較高溫度后，再讓其緩慢冷卻至常溫。該過(guò)程中固體中粒子不斷轉(zhuǎn)移，由無(wú)序到有序，內(nèi)能不斷降低，故固體退火趨勢(shì)可以由內(nèi)能函數(shù)判斷。同樣，在模擬退火過(guò)程中，需要反復(fù)擾動(dòng)參數(shù)場(chǎng)，并且需要目標(biāo)函數(shù)[22]判斷是否繼續(xù)擾動(dòng)。本文采用試井壓力數(shù)據(jù)作為評(píng)判指標(biāo)，目標(biāo)函數(shù)為

(1)

每次擾動(dòng)均需要重新計(jì)算目標(biāo)函數(shù)、對(duì)比擾動(dòng)前后的目標(biāo)函數(shù)值。一般模擬退火算法使用Metropolis準(zhǔn)則[17]判斷是否繼續(xù)擾動(dòng)，但本文中的初始模型是利用靜態(tài)資料根據(jù)傳統(tǒng)隨機(jī)建模方法建立的，能在一定程度上反映基本地質(zhì)特征，因此將退火溫度設(shè)置為接近0的狀態(tài)，簡(jiǎn)化模擬退火優(yōu)化過(guò)程。同時(shí)，引入MAP算法[23]判斷是否接受擾動(dòng)，當(dāng)目標(biāo)函數(shù)減少時(shí)接受擾動(dòng)；而目標(biāo)函數(shù)增大時(shí)拒絕擾動(dòng)，這樣會(huì)使算法較快達(dá)到收斂。其數(shù)學(xué)表達(dá)式為

(2)

式中：Paccept為接受擾動(dòng)的概率；Oold、Onew分別為擾動(dòng)前、后的目標(biāo)函數(shù)值。

1.2 算法優(yōu)化

目標(biāo)函數(shù)中的模擬壓力一般通過(guò)數(shù)值模擬計(jì)算得到，但由于模擬退火過(guò)程中需要對(duì)模型初始參數(shù)場(chǎng)進(jìn)行多次擾動(dòng)，如果每次擾動(dòng)都利用數(shù)值模擬計(jì)算井底壓力，則需要花費(fèi)大量時(shí)間。為減少數(shù)值模擬次數(shù)，本文引入偏心井的兩區(qū)復(fù)合油藏壓力公式計(jì)算壓力值[24-26]。

當(dāng)油井位于油藏內(nèi)區(qū)(圖1a)時(shí)，壓力計(jì)算公式為

0≤r≤r′

(3)

(4)

圖1 兩區(qū)均質(zhì)油藏偏心井模型示意圖

當(dāng)油井位于油藏外區(qū)(圖1b)時(shí)，壓力計(jì)算公式為

(5)

Kn(srD)cosn(θ-θ′)]r1≤r≤r′

(6)

(7)

由于油藏的非均質(zhì)性，計(jì)算結(jié)果會(huì)有誤差。但在擾動(dòng)次數(shù)較少時(shí)，壓力解析解與數(shù)值模擬計(jì)算結(jié)果誤差不大。隨著退火擾動(dòng)次數(shù)增加，誤差會(huì)不斷增大。因此在擾動(dòng)次數(shù)較多時(shí)，解析解不能滿足計(jì)算精度，中間需要穿插數(shù)值模擬進(jìn)行復(fù)算，修正計(jì)算誤差。本文采取每擾動(dòng)50次后利用數(shù)值模擬計(jì)算一次，即利用解析解和數(shù)值解相結(jié)合的方法計(jì)算單井的井底壓力。這樣既保證了計(jì)算的準(zhǔn)確性，又節(jié)約了運(yùn)算時(shí)間。

在目標(biāo)函數(shù)中引入計(jì)算的壓力值，優(yōu)化后的模擬退火目標(biāo)函數(shù)為

(8)

在式(8)的基礎(chǔ)上加入調(diào)節(jié)因子αi得到

(9)

式(9)可以加快計(jì)算速度[27]。

1.3 擾動(dòng)改進(jìn)

在擾動(dòng)階段，模擬退火方法通常是先在模型中隨機(jī)挑選一個(gè)網(wǎng)格，然后在數(shù)據(jù)集中隨機(jī)挑選新的孔隙度和滲透率數(shù)據(jù)替換原有的數(shù)據(jù)，擾動(dòng)過(guò)程中沒有范圍限制。這在退火溫度較高的模型中有助于跳出局部最優(yōu)而達(dá)到全局最優(yōu)，但在退火溫度很低時(shí)卻有可能增加無(wú)效擾動(dòng)次數(shù)。因此，需要對(duì)新數(shù)據(jù)的挑選做出一定限制，以減少擾動(dòng)次數(shù)，更快達(dá)到全局最優(yōu)。

本文擾動(dòng)替換網(wǎng)格采用以目標(biāo)網(wǎng)格為中心的7×7網(wǎng)格區(qū)域(根據(jù)模型的大小和精細(xì)程度可以適當(dāng)調(diào)整區(qū)域范圍)，以其中數(shù)據(jù)的最大值和最小值作為擾動(dòng)的上限和下限改進(jìn)擾動(dòng)過(guò)程。

2 計(jì)算步驟

結(jié)合試井動(dòng)態(tài)資料的隨機(jī)模擬算法主要包含擾動(dòng)、決策和更新三個(gè)過(guò)程。

擾動(dòng)過(guò)程就是擾動(dòng)儲(chǔ)層的參數(shù)場(chǎng)，即隨機(jī)改變儲(chǔ)層中某位置孔隙度和滲透率，從而導(dǎo)致模型計(jì)算的井底壓力產(chǎn)生變化，再利用解析方法計(jì)算擾動(dòng)后參數(shù)場(chǎng)的目標(biāo)函數(shù)。然后持續(xù)擾動(dòng)，當(dāng)達(dá)到退火溫度下的最大擾動(dòng)次數(shù)時(shí)，停止擾動(dòng)。

決策過(guò)程就是判斷擾動(dòng)過(guò)程中參數(shù)場(chǎng)擾動(dòng)前后目標(biāo)函數(shù)值的大小，從而判斷是否接受擾動(dòng)；判斷是否達(dá)到極限擾動(dòng)值；決定是否利用數(shù)模模擬值修正目標(biāo)函數(shù)；判斷是否滿足循環(huán)條件和停止條件。

更新過(guò)程就是更新網(wǎng)格擾動(dòng)次數(shù)和接受網(wǎng)格擾動(dòng)次數(shù)，更新儲(chǔ)層的滲透率、孔隙度的網(wǎng)格分布。如果接受擾動(dòng)，那么Oold將被Onew代替，而新的目標(biāo)函數(shù)與網(wǎng)格擾動(dòng)后的滲透率和孔隙度有關(guān)，目標(biāo)函數(shù)的準(zhǔn)確性將會(huì)影響到?jīng)Q策過(guò)程的可靠性。如果參數(shù)場(chǎng)改變了N次，那么利用解析解計(jì)算的目標(biāo)函數(shù)就有N次誤差(由于解析解采用儲(chǔ)層參數(shù)平均值導(dǎo)致)；因此，當(dāng)達(dá)到一定的擾動(dòng)次數(shù)后，引入數(shù)值模擬更新擾動(dòng)后參數(shù)場(chǎng)的壓力，重新計(jì)算目標(biāo)函數(shù)，使目標(biāo)函數(shù)更準(zhǔn)確。

結(jié)合試井動(dòng)態(tài)資料的隨機(jī)模擬方法的具體步驟如下。

(1)生成初始參數(shù)場(chǎng)并計(jì)算初始目標(biāo)函數(shù)Oi，令Ncal=0，Nswap=0， ΔOmax=0，Ntry=0。

(2)在模型中隨機(jī)挑選一個(gè)網(wǎng)格，使用Gibbs采樣[28-29]對(duì)參數(shù)場(chǎng)進(jìn)行擾動(dòng)，在擾動(dòng)過(guò)程中，根據(jù)鄰近網(wǎng)格屬性值設(shè)定擾動(dòng)范圍，加快收斂速度，重新計(jì)算目標(biāo)函數(shù)O，Ntry=Ntry+1。

(3)內(nèi)循環(huán)。判斷ΔO與ΔOmax的大小，若ΔO≥ΔOmax，則拒絕擾動(dòng)；再判斷Ntry是否到達(dá)Nmax，若未達(dá)到則返回步驟(2)，若達(dá)到則終止循環(huán)并重新開始步驟(1)；若ΔO<ΔOmax，則接受擾動(dòng)，轉(zhuǎn)到步驟(4)。

(4)根據(jù)新的參數(shù)場(chǎng)，利用壓力公式計(jì)算目標(biāo)函數(shù)O，則Ncal=Ncal+1，Nswap=Nswap+1。

(5)判斷Ncal是否等于Nsmax，若相等，則利用模擬器計(jì)算目標(biāo)函數(shù)O，Ncal=0。判斷O與Omin的大小，若O

(6)外循環(huán)。判斷Nswap與Ntot的大小，若Nswap>Ntot，則輸出結(jié)果；若Nswap

上述步驟中：Ncal為用解析解計(jì)算模擬壓力次數(shù)；Nswap為互換網(wǎng)格次數(shù)；Ntry為當(dāng)前退火溫度下擾動(dòng)次數(shù)；Nmax為當(dāng)前退火溫度下最大擾動(dòng)次數(shù)；Nsmax為用解析解計(jì)算模擬壓力的極限次數(shù)；Ntot為總的接受擾動(dòng)的最大次數(shù)； ΔO為目標(biāo)函數(shù)差值；Omin為最小目標(biāo)函數(shù)值。

3 實(shí)例分析

選擇某注水開發(fā)油田的一反五點(diǎn)注采井網(wǎng)模型為實(shí)例，研究區(qū)范圍為800m×800m，共四口采油井(A、B、C和D井)和一口注水井(I井)，相關(guān)參數(shù)見表1。

表1 儲(chǔ)層參數(shù)表

利用五口井測(cè)井解釋的孔隙度、滲透率建立該井組的屬性模型。為了簡(jiǎn)化研究，模型垂向上假設(shè)為1層，平面網(wǎng)格步長(zhǎng)為20m，總網(wǎng)格數(shù)為1600(40×40×1)個(gè)，采用序貫高斯模擬方法建立初始孔隙度、滲透率模型如圖2所示。

圖2 初始孔隙度(a)和滲透率(b)模型

建立初始模型參數(shù)場(chǎng)后，需要利用試井測(cè)試的井底壓力對(duì)模型的參數(shù)場(chǎng)進(jìn)行修正。研究區(qū)四口采油井均有壓力降落試井資料。利用A、B、C三口井的資料構(gòu)建目標(biāo)函數(shù)，采用前述模擬退火算法擬合壓力數(shù)據(jù)，將D井作為后驗(yàn)井進(jìn)行符合率分析、驗(yàn)證。

由圖3可以看出，初始模型計(jì)算的各井的井底壓力與測(cè)試壓力差別較大，因此有必要進(jìn)行井底壓力擬合。在井底壓力擬合過(guò)程中，共經(jīng)歷675次擾動(dòng)，接受擾動(dòng)次數(shù)為352次。若擾動(dòng)未經(jīng)優(yōu)化，則需擾動(dòng)1000次以上，因此擾動(dòng)接受效率約提高48%。

用Eclipse數(shù)模軟件進(jìn)行7次數(shù)值模擬計(jì)算，計(jì)算時(shí)間約為15.0min，其中每次擾動(dòng)及判定需要1s，共用時(shí)11.5min，每次數(shù)值模擬需要30s，共用時(shí)3.5min；如不利用解析解對(duì)壓力進(jìn)行近似計(jì)算，數(shù)值模擬的計(jì)算次數(shù)將達(dá)到352次，增加345次的數(shù)模計(jì)算時(shí)間，共用時(shí)約為173min。

最終使模型計(jì)算的壓力與試井測(cè)試的壓力基本一致，其中A、B、C三口井平均誤差在0.7%以內(nèi)，后驗(yàn)井D井誤差也在1%以內(nèi)(表2)，滿足精度要求。

表2 優(yōu)化前后實(shí)測(cè)值與計(jì)算值平均誤差%

井名優(yōu)化前優(yōu)化后先驗(yàn)井A1.980.47先驗(yàn)井B2.230.60先驗(yàn)井C13.770.97后驗(yàn)井D1.400.87

圖4為經(jīng)過(guò)試井壓力擬合后的滲透率和孔隙度模型，與擬合前的模型(圖2)相比，經(jīng)過(guò)模擬退火的擾動(dòng)后，儲(chǔ)層參數(shù)的整體分布趨勢(shì)基本一致。修正后的模型不僅能反映儲(chǔ)層靜態(tài)參數(shù)的分布特征，同時(shí)新模型很好地符合了試井的壓力數(shù)據(jù)，比原模型更接近真實(shí)的地質(zhì)模型。

圖3 A井(a)、B井(b)、C井(c)和D井(d)的壓力變化曲線

圖4 經(jīng)過(guò)試井壓力擬合(本文方法)后的孔隙度(a)和滲透率(b)模型

4 結(jié)論

在傳統(tǒng)隨機(jī)建模方法基礎(chǔ)上，結(jié)合試井壓力數(shù)據(jù)對(duì)初始模型加以修正，使模型符合動(dòng)態(tài)資料，可以提高模型精度。

(1)以試井壓力數(shù)據(jù)作為評(píng)判標(biāo)準(zhǔn)建立目標(biāo)函數(shù)可以使優(yōu)化后的模型反映壓力動(dòng)態(tài)。

(2)引入偏心井的兩區(qū)復(fù)合油藏壓力公式，交替使用解析公式和數(shù)值模擬，既可保證計(jì)算精度，又可減少數(shù)模次數(shù)，節(jié)約計(jì)算時(shí)間，提高了實(shí)用性。

(3)對(duì)模擬退火中的擾動(dòng)過(guò)程進(jìn)行了改進(jìn)，在每次擾動(dòng)過(guò)程中，以鄰近的網(wǎng)格作為參考，對(duì)擾動(dòng)的范圍進(jìn)行了限制，擾動(dòng)后被接受的效率提高了約48%，有效減少了擾動(dòng)次數(shù)。