一種強(qiáng)非均質(zhì)性油藏自動(dòng)歷史擬合混合求解方法

張 凱, 馬小鵬, 王增飛, 劉 凡, 馬 瑋, 姚 軍

(1.中國(guó)石油大學(xué)(華東)石油工程學(xué)院,山東青島 266580; 2.北京中油瑞飛信息技術(shù)有限公司,北京 100007;3.中海油研究總院，北京 100010; 4.海洋石油高效開(kāi)發(fā)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100028;5.勝利油田油氣井下作業(yè)中心，山東東營(yíng) 257001)

目前中國(guó)大多數(shù)油田經(jīng)過(guò)長(zhǎng)期的勘探開(kāi)發(fā),已經(jīng)進(jìn)入注水開(kāi)發(fā)的中后期階段,尤其是砂巖油氣藏儲(chǔ)層之間的非均質(zhì)特征非常強(qiáng)[1-3],嚴(yán)重影響了油藏注水開(kāi)發(fā)效果。開(kāi)展強(qiáng)非均質(zhì)性油藏中高滲通道識(shí)別技術(shù)的研究,對(duì)于油氣資源開(kāi)發(fā)具有重要意義。Brigham等[4-9]利用油藏動(dòng)靜態(tài)資料等識(shí)別復(fù)雜油氣藏中的高滲通道,方法成本高,適用性局限。自動(dòng)歷史擬合技術(shù)依靠油藏?cái)?shù)值模擬,根據(jù)油藏動(dòng)態(tài)數(shù)據(jù)反演調(diào)整地質(zhì)參數(shù),具有很強(qiáng)的可靠性。目前油藏自動(dòng)歷史擬合方法眾多[10-15],但是缺乏專門針對(duì)強(qiáng)非均質(zhì)性油藏的方法。筆者對(duì)強(qiáng)非均質(zhì)性油藏自動(dòng)歷史擬合方法展開(kāi)研究,提出混合求解方法:PCA方法[16-19]與DCT方法[20-23]混合;SPSA算法[24-27]與ABC算法[28-29]混合。

1 強(qiáng)非均質(zhì)性油藏自動(dòng)歷史擬合方法

1.1 歷史擬合數(shù)學(xué)模型的建立

自動(dòng)歷史擬合問(wèn)題是一個(gè)典型的反問(wèn)題,通常結(jié)合先驗(yàn)地質(zhì)信息進(jìn)行求解能夠得到接近準(zhǔn)確解的結(jié)果[11]。為了得到準(zhǔn)確的油藏地質(zhì)特征,需要利用已知的觀測(cè)數(shù)據(jù)與油藏模型參數(shù)m之間的關(guān)系建立數(shù)學(xué)模型,表達(dá)式為

dobs=g(m)+εr.

(1)

式中,dobs為生產(chǎn)動(dòng)態(tài)數(shù)據(jù);g(·)為油藏?cái)?shù)值模擬器;εr為觀測(cè)數(shù)據(jù)誤差。

為了求解反問(wèn)題,通常對(duì)原函數(shù)進(jìn)行正則化處理,利用先驗(yàn)信息作為約束,對(duì)所求問(wèn)題作合適的形式轉(zhuǎn)換,使反問(wèn)題的解更加準(zhǔn)確。其中貝葉斯理論就是正則化方法之一,基于貝葉斯理論,觀測(cè)數(shù)據(jù)與油藏模型參數(shù)之間關(guān)系式[30]為

P(m|dobs)∝P(dobs|m)P(m).

(2)

式中,m為油藏模型參數(shù);P(m|dobs)為貝葉斯理論中的后驗(yàn)概率;P(dobs|m)為貝葉斯理論中的似然函數(shù);P(m)為貝葉斯理論中的先驗(yàn)概率。

在地質(zhì)參數(shù)m給定的條件下,獲得的觀測(cè)數(shù)據(jù)存在觀測(cè)誤差εr,一般認(rèn)為誤差符合均值為0,協(xié)方差矩陣為CD的高斯型概率分布,即εr～(0,CD)。因此可得到似然函數(shù)P(dobs|m)為

(3)

為了求得先驗(yàn)概率,須利用已知的油藏信息建立先驗(yàn)?zāi)Ｐ?。本文中采用隨機(jī)建模[21]方法中的序貫高斯模擬方法生成多個(gè)可選的、等概率的先驗(yàn)?zāi)Ｐ汀Ｄ壳半S機(jī)建模的方法包括序貫高斯模擬、截?cái)喔咚鼓M以及分形模擬等方法[31-32]。其中序貫高斯模擬主要運(yùn)用高斯模型,根據(jù)離散的觀測(cè)點(diǎn)生成一個(gè)連續(xù)的隨機(jī)面,地質(zhì)參數(shù)符合高斯分布規(guī)律,能夠很好地反映實(shí)際地質(zhì)情況,在構(gòu)建地質(zhì)模型領(lǐng)域得到了廣泛的應(yīng)用。因此可以得到先驗(yàn)概率為

(4)

根據(jù)貝葉斯理論可以求得后驗(yàn)概率為

P(m|d)∝exp[-O(m)].

(5)

其中

式中,O(m)為歷史擬合數(shù)學(xué)模型的目標(biāo)函數(shù);mpr為先驗(yàn)?zāi)Ｐ蛥?shù)。

后驗(yàn)概率P(m|dobs)中不僅包含了觀測(cè)數(shù)據(jù),還包含了先驗(yàn)地質(zhì)信息。則歷史擬合反演問(wèn)題轉(zhuǎn)化為使后驗(yàn)概率P(m|dobs)最大化問(wèn)題,此過(guò)程稱為MAP(Maximum A Posteriori)。后驗(yàn)概率越大,自動(dòng)歷史擬合所獲得模型參數(shù)越能反映真實(shí)的地質(zhì)情況。以O(shè)(m)作為目標(biāo)函數(shù),則問(wèn)題就轉(zhuǎn)換成求解目標(biāo)函數(shù)最小值的問(wèn)題。

1.2 PCA特征提取方法

主成分分析方法(principal component analysis,PCA)是目前應(yīng)用很廣泛的一種特征提取方法,在數(shù)據(jù)降維和特征提取方面的有效性使得它在人臉識(shí)別領(lǐng)域獲得了廣泛的應(yīng)用[16-17,19]。它主要利用K-L變換從原始的數(shù)據(jù)集中提取主要特征并構(gòu)成特征向量空間,而后將原有的數(shù)據(jù)影射到特征向量空間上并得到一組投影系數(shù),能夠很好地保留原始數(shù)據(jù)的主要特征信息。利用PCA方法對(duì)先驗(yàn)?zāi)Ｐ偷姆蔷|(zhì)特征進(jìn)行提取,能夠很好地提取出先驗(yàn)?zāi)Ｐ椭懈?低)滲區(qū)域的信息,從而準(zhǔn)確地反演出強(qiáng)非均質(zhì)油藏中高(低)滲區(qū)域。

首先根據(jù)先驗(yàn)信息生成Ne個(gè)先驗(yàn)?zāi)Ｐ?記為mj(j=1,2,…,Ne),平均值為mpr。對(duì)于實(shí)際油藏的歷史擬合問(wèn)題,由于先驗(yàn)?zāi)Ｐ蜐M足序列高斯分布,則其平均模型為

(6)

(7)

(8)

直接求CM的特征值和特征向量計(jì)算量較大,因此求解矩陣V=ATA的特征值λi及對(duì)應(yīng)的正交化特征向量vi,減小計(jì)算量的同時(shí)更加準(zhǔn)確地反映出先驗(yàn)?zāi)Ｐ蛥?shù)場(chǎng)的特征,選取前p(p≤Ne?Nm)個(gè)最大特征值及對(duì)應(yīng)的特征向量;這些特征向量包含了先驗(yàn)?zāi)Ｐ椭械闹饕卣餍畔?。利用?8)求取協(xié)方差矩陣CM的正交歸一化特征向量,獲得特征向量提取矩陣Φ即為原參數(shù)場(chǎng)的特征提取場(chǎng),

Φ=[φ1,φ2,…,φe],

(9)

(10)

對(duì)圖1表示的對(duì)數(shù)滲透率場(chǎng)進(jìn)行特征提取,可見(jiàn)最大特征值λ1對(duì)應(yīng)的特征提取場(chǎng)φ1中存在1條明顯的低(高)滲通道,特征值為λ5對(duì)應(yīng)的特征提取場(chǎng)φ5中的高(低)滲通道條數(shù)增加到3條。隨著特征值的減小,對(duì)應(yīng)的特征提取場(chǎng)中的高(低)滲通道數(shù)越來(lái)越多,因此選擇對(duì)應(yīng)特征值大的特征提取場(chǎng)進(jìn)行歷史擬合能夠減少先驗(yàn)?zāi)Ｐ椭袩o(wú)關(guān)信息的干擾,更加準(zhǔn)確地反演出高(低)滲通道。

圖1 PCA方法特征提取場(chǎng)Fig.1 PCA method to extract feature field

根據(jù)PCA方法的K-L變換原理,目標(biāo)函數(shù)中協(xié)方差矩陣的逆矩陣可以被近似表示為

(11)

實(shí)現(xiàn)特征提取場(chǎng)Φ代替先驗(yàn)?zāi)Ｐ蛥?shù)場(chǎng)M進(jìn)行自動(dòng)歷史擬合。此時(shí),定義新的參數(shù)s(s∈Rp)對(duì)油藏模型參數(shù)進(jìn)行參數(shù)化變換,

(12)

則目標(biāo)函數(shù)變?yōu)?/p>

(13)

其中

使用PCA特征提取方法可將油藏模型的歷史擬合問(wèn)題由Nm維降低到p維,并且有效避免了協(xié)方差矩陣求逆的過(guò)程,大規(guī)模油藏進(jìn)行自動(dòng)歷史擬合時(shí)由于擬合參數(shù)的維數(shù)過(guò)于龐大,CM的逆矩陣計(jì)算十分耗時(shí),避免目標(biāo)函數(shù)中協(xié)方差矩陣CM的逆矩陣計(jì)算能夠大幅度提高計(jì)算速度。結(jié)合優(yōu)化求解方法,通過(guò)更新參數(shù)s降低目標(biāo)函數(shù)值,并基于式(13)反求油藏參數(shù),實(shí)現(xiàn)強(qiáng)非均質(zhì)性油藏歷史擬合問(wèn)題的求解。

1.3 基于離散余弦變換的PCA降維方法

PCA方法對(duì)于高滲通道區(qū)分明顯的油藏模型歷史擬合,能夠很好地反演出接近真實(shí)地質(zhì)情況的油藏模型。但對(duì)于地質(zhì)參數(shù)如孔隙度、滲透率等分布復(fù)雜的油氣藏,PCA方法并不能夠很好地反演出準(zhǔn)確的儲(chǔ)層參數(shù),因此在PCA方法的基礎(chǔ)上提出PCA方法與DCT方法混合進(jìn)行特征提取。離散余弦變換方法(discrete cosine transform,DCT)是將一系列的離散點(diǎn)表示為頻率不同的余弦函數(shù)的疊加。在圖像處理領(lǐng)域,利用DCT方法良好的能量壓縮性能,通過(guò)保留大系數(shù)來(lái)保存大部分的圖像信息。在進(jìn)行自動(dòng)歷史擬合時(shí)選取低頻系數(shù)對(duì)先驗(yàn)?zāi)Ｐ瓦M(jìn)行降維,能夠保留先驗(yàn)?zāi)Ｐ偷拇蟛糠中畔?從而更加精準(zhǔn)地描述油藏地質(zhì)信息。

對(duì)油藏先驗(yàn)?zāi)Ｐ蛥?shù)場(chǎng)M進(jìn)行離散余弦變換,D是參數(shù)場(chǎng)離散余弦變換后的DCT系數(shù)矩陣。定義DCT變換矩陣W1和W2,其元素分別定義為

(14)

(15)

(16)

對(duì)圖2表示的滲透率場(chǎng)進(jìn)行DCT變換得到降維后的特征提取場(chǎng),該變換能夠更好的保留先驗(yàn)數(shù)據(jù)的信息,增強(qiáng)了歷史擬合的魯棒性。

本文中利用PCA和DCT結(jié)合進(jìn)行先驗(yàn)?zāi)Ｐ偷奶卣魈崛?PCA特征提取的參數(shù)場(chǎng)ΦP為(p/2)×Ne維,DCT變換特征提取場(chǎng)ΦD為(p/2)×Ne維。將兩者組合得到特征提取場(chǎng)Φ代入原來(lái)的PCA方法過(guò)程即可實(shí)現(xiàn)PCA方法與DCT方法的混合,

Φ=[ΦP;ΦD].

(17)

圖2 DCT方法特征提取Fig.2 Discrete cosine transform to extract feature

1.4 歷史擬合目標(biāo)函數(shù)的優(yōu)化求解

隨機(jī)擾動(dòng)梯度近似算法[24-26](simultaneous perturbation stochastic approximation,SPSA)是通過(guò)添加擾動(dòng)和重復(fù)迭代方法尋找線性或非線性目標(biāo)函數(shù)的局部最優(yōu)值,適用于多維變量系統(tǒng),同其他隨機(jī)優(yōu)化方法(如遺傳算法)比較,可以與任何數(shù)值模擬器耦合。在實(shí)際應(yīng)用中算法的近似梯度與實(shí)際梯度方向十分接近,具有很強(qiáng)的收斂性。但在迭代的過(guò)程中,由于步長(zhǎng)的選取不當(dāng),可能會(huì)導(dǎo)致目標(biāo)函數(shù)陷入局部最優(yōu)。本文中提出人工蜂群算法(ABC)與SPSA算法混合對(duì)目標(biāo)函數(shù)進(jìn)行求解。ABC算法[28-29]是一種較為系統(tǒng)的群集智能隨機(jī)優(yōu)化算法,具有操作簡(jiǎn)單、控制參數(shù)少、搜索精度較高和魯棒性較強(qiáng)的特點(diǎn),全局搜索能力較好,但局部搜索能力較差。將ABC算法的全局搜索模式引入SPSA算法中,可以提高SPSA算法的全局搜索能力。

(18)

(19)

(20)

在迭代的過(guò)程中,由于步長(zhǎng)的選取具有局限性,可能會(huì)導(dǎo)致目標(biāo)函數(shù)陷入局部最優(yōu),此時(shí)需要引入ABC算法的搜索模式尋求全局最優(yōu)值。在ABC算法中如果在一個(gè)食物源(問(wèn)題的解)的位置處不能在預(yù)先設(shè)定的循環(huán)次數(shù)limit內(nèi)找到更優(yōu)的解,則該食物源被放棄并且被偵察蜂找到的新食物源所代替,偵察蜂搜索新的食物源依據(jù)為

sk+1=smin+rand[0,1](smax-smin).

(21)

將該搜索模式引入SPSA算法中,改進(jìn)SPSA算法流程如圖3所示。

圖3 SPSA與ABC算法混合優(yōu)化流程Fig.3 SPSA and ABC algorithm hybrid optimization flow chart

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

2.1 PCA方法自動(dòng)歷史擬合實(shí)例

對(duì)一個(gè)油氣水三相的油藏模型進(jìn)行PCA方法歷史擬合研究,采用Eclipse軟件作為油藏?cái)?shù)值模擬器。該油藏存在三條明顯的高滲通道,需要反演調(diào)整的地質(zhì)參數(shù)為滲透率。用序貫高斯模擬隨機(jī)生成100個(gè)先驗(yàn)?zāi)Ｐ?取部分如圖4所示。利用PCA方法對(duì)目標(biāo)函數(shù)中的先驗(yàn)?zāi)Ｐ蛥?shù)矩陣進(jìn)行特征提取,利用ABC算法與SPSA算法混合進(jìn)行優(yōu)化求解,滲透率場(chǎng)的反演結(jié)果如圖5所示?？梢钥闯?當(dāng)?shù)郊s第5步時(shí),注水井Inj-1、生產(chǎn)井Pro-2、Pro-5連線上出現(xiàn)了明顯的高滲條帶,當(dāng)?shù)?5步時(shí),生產(chǎn)井Pro-4、Pro-3連線上同樣出現(xiàn)了一條明顯的滲透率條帶。由于生產(chǎn)井Pro-6是單口井生產(chǎn),與其他井之間不存在連通性,所以滲透率值變化很慢,直到迭代到25步時(shí),才獲得一條模糊的滲透率條帶,至此反演所獲得的滲透率場(chǎng)已經(jīng)非常接近真實(shí)模型的滲透率場(chǎng)。結(jié)果表明,利用PCA方法提取地質(zhì)參數(shù)中的非均質(zhì)特征能夠很好地反演出油藏中的高(低)滲通道;同時(shí)將高維地質(zhì)參數(shù)數(shù)據(jù)轉(zhuǎn)化為低維參數(shù)使收斂速度更快。

2.2 PCA結(jié)合DCT方法歷史擬合實(shí)例

以一個(gè)油氣水三相的復(fù)雜油藏模型,采用Eclipse軟件作為油藏?cái)?shù)值模擬器,進(jìn)行PCA與DCT混合特征提取方法歷史擬合研究,對(duì)滲透率場(chǎng)進(jìn)行反演。采用序貫高斯模擬隨機(jī)生成100個(gè)先驗(yàn)?zāi)Ｐ?選取部分如圖6所示。然后利用PCA與DCT方法混合對(duì)先驗(yàn)?zāi)Ｐ瓦M(jìn)行特征提取,最后利用SPSA與ABC算法混合進(jìn)行求解,結(jié)果如圖7所示,可以看出求得的滲透率場(chǎng)在生產(chǎn)井Pro-1、Pro-4和注水井Inj-5、Inj-9的連線上出現(xiàn)了高滲條帶;在生產(chǎn)井Pro-3和注水井Inj-4、Inj-7以及生產(chǎn)井Pro-2和注水井Inj-2、Inj-6連線上也出現(xiàn)了高滲條帶。在交叉井之間,利用該方法反演出的滲透率值也十分接近真實(shí)滲透率場(chǎng)。對(duì)比PCA與DCT方法混合與傳統(tǒng)的SVD方法進(jìn)行歷史擬合所獲得的結(jié)果,可以看出相比PCA結(jié)合DCT方法,SVD方法反演的滲透率場(chǎng)中高滲區(qū)域并不明顯,與真實(shí)滲透率場(chǎng)的差距較大。與單一的PCA方法進(jìn)行歷史擬合相比,PCA與DCT混合的方法能夠獲得更加精細(xì)的擬合結(jié)果,滲透率場(chǎng)更加接近真實(shí)地質(zhì)情況。

圖4 初始隨機(jī)滲透率場(chǎng)實(shí)現(xiàn)Fig.4 Initial random permeability field

圖5 滲透率場(chǎng)反演過(guò)程Fig.5 Permeability field inversion process

圖6 隨機(jī)滲透率場(chǎng)實(shí)現(xiàn)Fig.6 Random permeability field

圖7 滲透率場(chǎng)反演結(jié)果Fig.7 Results of history matching

圖8和9分別為累積指標(biāo)和單井指標(biāo)擬合結(jié)果。由圖8可以看出,PCA與DCT方法混合擬合得到的累積產(chǎn)水量與地層壓力隨時(shí)間的變化曲線能夠很好地反映真實(shí)生產(chǎn)數(shù)據(jù)的變化;由圖9可以看出,擬合得到的生產(chǎn)井Pro-1的日產(chǎn)油量以及生產(chǎn)井Pro-3的日產(chǎn)水量隨時(shí)間的變化曲線很好地反映了真實(shí)生產(chǎn)數(shù)據(jù)的變化。滲透率場(chǎng)與生產(chǎn)數(shù)據(jù)的擬合結(jié)果基本滿足實(shí)際需要,因此提出的PCA結(jié)合DCT方法可以作為一種有效的方法進(jìn)行復(fù)雜強(qiáng)非均質(zhì)性油藏的自動(dòng)歷史擬合。

圖8 累積指標(biāo)的擬合結(jié)果Fig.8 History matching results of cumulative index

圖9 單井指標(biāo)的擬合結(jié)果Fig.9 History matching results of single well index

3 結(jié) 論

(1)將SPSA算法與ABC算法混合對(duì)目標(biāo)函數(shù)進(jìn)行求解,SPSA算法能夠很好地尋找目標(biāo)函數(shù)的局部最優(yōu)值,近似梯度與實(shí)際梯度方向十分接近使SPSA算法具有很強(qiáng)的收斂性,引入ABC算法的全局搜索模式提高SPSA算法的全局搜索能力。SPSA算法與ABC算法混合提高了模型求解速度與精度。

(2)利用PCA特征提取的方法基本能夠識(shí)別出油藏中的高滲通道,同時(shí)將高維地質(zhì)參數(shù)數(shù)據(jù)轉(zhuǎn)化為低維參數(shù)使優(yōu)化算法收斂速度更快;但是由于單一的PCA方法不能提取出較為全面的儲(chǔ)層信息,滲透率場(chǎng)反演結(jié)果與真實(shí)滲透率場(chǎng)存在較大偏差。

(3)PCA和DCT方法混合能夠提取更為全面的地質(zhì)特征,并且提高歷史擬合的魯棒性,很好地反演出油藏中的高滲通道區(qū)域并且反演結(jié)果更加接近真實(shí)場(chǎng)。提出的混合求解方法可以作為一種有效的方法進(jìn)行該類油藏的自動(dòng)歷史擬合。