王澤龍 劉先貴 李志勇 張 輝 張傳進(jìn)

（1. 中國(guó)石油國(guó)際勘探開(kāi)發(fā)有限公司；2. 中國(guó)石油勘探開(kāi)發(fā)研究院）

隨著油氣行業(yè)的智能化轉(zhuǎn)型，在油藏的精細(xì)化管理和生產(chǎn)過(guò)程中，人們?cè)絹?lái)越重視數(shù)字化信息的實(shí)時(shí)采集，并依靠計(jì)算機(jī)對(duì)各類(lèi)生產(chǎn)數(shù)據(jù)的分析，不斷修改油藏模型，及時(shí)調(diào)整油田開(kāi)發(fā)技術(shù)政策，以便實(shí)現(xiàn)油藏的經(jīng)濟(jì)、高效開(kāi)發(fā)[1]。目前，油藏歷史擬合技術(shù)一般只能利用井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)[2]，通過(guò)人工手動(dòng)或計(jì)算機(jī)輔助的方法修改各類(lèi)參數(shù)，將油藏?cái)?shù)值模擬得到的預(yù)測(cè)井口生產(chǎn)數(shù)據(jù)與井口測(cè)量的數(shù)據(jù)進(jìn)行擬合，以便得到更為可靠的油藏模型。但由于井口生產(chǎn)數(shù)據(jù)的空間覆蓋范圍十分有限，油藏歷史擬合缺少足夠的約束條件，以此為基礎(chǔ)修正的油藏模型往往不符合實(shí)際地質(zhì)情況。

油藏歷史擬合問(wèn)題的數(shù)學(xué)本質(zhì)是求解最優(yōu)化。目前，梯度類(lèi)算法和集合類(lèi)算法是解決油藏歷史擬合問(wèn)題的兩類(lèi)常見(jiàn)的優(yōu)化算法。其中，Levenberg-Marquardt（LM）算法是常用的梯度類(lèi)算法，該算法用于小型反演問(wèn)題效率很高，但對(duì)于大型油藏模型由于需要計(jì)算和存儲(chǔ)超高維度的敏感度矩陣，造成計(jì)算量過(guò)于巨大而難以廣泛應(yīng)用[3]。集合類(lèi)算法是基于貝葉斯理論的后驗(yàn)概率來(lái)進(jìn)行最優(yōu)化求解，1960年，Kalman[4]提出卡爾曼濾波，通過(guò)對(duì)平均值和均方誤差的最小化問(wèn)題求解，來(lái)實(shí)現(xiàn)對(duì)某一過(guò)程中的狀態(tài)參數(shù)高效地進(jìn)行循環(huán)貝葉斯概率推斷；1970年，Buck[5]和Sunahara[6]將卡爾曼濾波擴(kuò)展到非線性問(wèn)題，提出擴(kuò)展卡爾曼濾波（EKF）；1994年，Evensen[7]提出集合卡爾曼濾波（EnKF），將卡爾曼濾波與蒙特卡洛方法[8]相結(jié)合，形成了集合類(lèi)歷史擬合算法的雛形；2002年，N?vdal[9]等第一次將EnKF引入石油工程領(lǐng)域，用于更新油井附近的油藏模型，隨后 Lorentzen[10]、Chen[11]、Li[12]、Jafarpour[13]等學(xué)者加入研究，使得EnKF算法成為當(dāng)前油藏歷史擬合最熱門(mén)的研究領(lǐng)域；2013年，結(jié)合集合平滑算法與EnKF，Emerick[14]提出基于多次數(shù)據(jù)吸收的集合平滑算法（ES-MDA），為油藏歷史擬合問(wèn)題提供了一種新思路。

四維地震監(jiān)測(cè)是在同一個(gè)地震采集區(qū)，間隔不同的時(shí)間部署三維地震，通過(guò)分析不同時(shí)間點(diǎn)的三維地震屬性差異來(lái)監(jiān)測(cè)油藏動(dòng)態(tài)情況，可為油氣田精細(xì)開(kāi)發(fā)管理提供覆蓋面廣、信息量大的測(cè)量數(shù)據(jù)資料，據(jù)此計(jì)算出油藏壓力和溫度、儲(chǔ)層生產(chǎn)能力和流體變化，對(duì)明確剩余油氣分布、提高采收率有著重要意義。由于四維地震數(shù)據(jù)的復(fù)雜程度及與井口生產(chǎn)數(shù)據(jù)存在差異，當(dāng)前的油藏歷史擬合尚難以實(shí)現(xiàn)對(duì)四維地震數(shù)據(jù)的有效利用。本文研究提出了一種基于四維地震數(shù)據(jù)的歷史擬合新方法，通過(guò)對(duì)四維地震數(shù)據(jù)的粗化與稀疏化處理，實(shí)現(xiàn)地震數(shù)據(jù)與井口數(shù)據(jù)的匹配，將由四維地震監(jiān)測(cè)獲得的縱波波阻抗數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)相結(jié)合，并通過(guò)建立起油藏狀態(tài)與地震波響應(yīng)的數(shù)學(xué)模型，形成油藏流體—地震屬性耦合模擬程序，改進(jìn)ES-MDA算法，實(shí)現(xiàn)四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)的聯(lián)合油藏歷史擬合。

1 常規(guī)油藏歷史擬合算法

1.1 油藏歷史擬合的目標(biāo)函數(shù)

油藏自動(dòng)歷史擬合問(wèn)題可以轉(zhuǎn)為在約束條件下的最優(yōu)化問(wèn)題。這個(gè)最優(yōu)化問(wèn)題的目標(biāo)函數(shù)可以用最小二乘法誤差表示[15]，定義如下：

式中m——油藏地質(zhì)模型參數(shù)矩陣；

mf——修改后的油藏地質(zhì)模型參數(shù)矩陣；

f（x）——油藏流體數(shù)值模擬程序；

duc——實(shí)際測(cè)量的數(shù)據(jù)向量；

CM——油藏地質(zhì)模型參數(shù)m的協(xié)方差矩陣；

CD—— 計(jì)算結(jié)果f（m）和實(shí)測(cè)數(shù)據(jù)duc的協(xié)方差矩陣。

目標(biāo)函數(shù)分為兩部分，每個(gè)部分的數(shù)值均大于或等于零。為使目標(biāo)函數(shù)最小化，既需要模擬計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)相吻合，即f（m）-duc趨近于零，也需要修改后的地質(zhì)模型不要偏離初始模型太遠(yuǎn)，即f（mmf）趨近于零，用來(lái)表征更新過(guò)程受到足夠的地質(zhì)資料約束。

1.2 EnKF算法

對(duì)于油藏歷史擬合，定義擴(kuò)展?fàn)顟B(tài)向量xk為k時(shí)刻估算的狀態(tài)參數(shù)向量，其中包括了靜態(tài)地質(zhì)模型的參數(shù)、主要油藏動(dòng)態(tài)參數(shù)與實(shí)測(cè)數(shù)據(jù)[16]。

式中mk—— 靜態(tài)地質(zhì)模型參數(shù)（如孔隙度、滲透率）在時(shí)間k時(shí)的估值，其維度是Nm×1；

Pk—— 動(dòng)態(tài)參數(shù)（如油藏壓力、油氣水飽和度），其維度是Np×1；

dk—— 實(shí)測(cè)數(shù)據(jù)（如油氣產(chǎn)量、產(chǎn)水量、井底壓力或四維地震數(shù)據(jù)），其維度是Nn×1 ；

Nm，Np，Nn—— 分別是單個(gè)地質(zhì)模型中的靜態(tài)參數(shù)數(shù)量、單一時(shí)刻的動(dòng)態(tài)參數(shù)數(shù)量、單一時(shí)刻的實(shí)測(cè)數(shù)據(jù)量；

下角標(biāo)k—— 時(shí)間序列中的第k項(xiàng)，即第k時(shí)間步（k為整數(shù)）。

對(duì)于油藏流體數(shù)值模擬的過(guò)程定義如下：

上角標(biāo)f——油藏模型預(yù)測(cè)的數(shù)據(jù)；

上角標(biāo)u——更新后的數(shù)據(jù)。

至此，完成了在k時(shí)刻的后驗(yàn)?zāi)Ｐ图系墓浪?。通過(guò)不斷重復(fù)此過(guò)程，對(duì)比不同時(shí)刻的實(shí)際測(cè)量數(shù)據(jù)，并不斷更新油藏地質(zhì)模型集合。EnKF的算法流程如圖1所示。

圖1 EnKF算法流程圖

EnKF算法具有如下明顯的缺點(diǎn)：（1）應(yīng)用效果十分依賴(lài)于初始模型的質(zhì)量；（2）在EnKF更新方程中，模型參數(shù)和狀態(tài)參數(shù)同時(shí)更新，難免造成更新后的狀態(tài)參數(shù)與實(shí)際的狀態(tài)參數(shù)不一致。

1.3 ES-MDA算法

在控制優(yōu)化理論中，對(duì)時(shí)間序列的預(yù)測(cè)算法除濾波類(lèi)算法外，還有平滑（smoother）算法，即利用過(guò)去、當(dāng)前和未來(lái)的信息來(lái)估算某個(gè)時(shí)間點(diǎn)的狀態(tài)參量。如集合平滑算法（ES）對(duì)全部時(shí)間序列數(shù)據(jù)整體吸收的方式，避免了EnKF算法更新后的狀態(tài)參數(shù)與實(shí)際的狀態(tài)參數(shù)不一致的問(wèn)題，且對(duì)線性問(wèn)題的計(jì)算效率遠(yuǎn)高于EnKF。2013年，Emerick[14]等受EnKF多次迭代的啟發(fā)，將EnKF與ES相結(jié)合，通過(guò)對(duì)相同數(shù)據(jù)的重復(fù)吸收來(lái)減小修正步長(zhǎng)，形成基于多次數(shù)據(jù)吸收的集合平滑算法（ES-MDA），并實(shí)現(xiàn)了對(duì)非線性問(wèn)題的應(yīng)用。

按照上述3個(gè)步驟，使用ES-MDA更新方程對(duì)整個(gè)集合中的所有子集模型都完成修正后，即完成了一次數(shù)據(jù)吸收。重復(fù)以上步驟Na次，即是ES-MDA算法的全部過(guò)程。ES-MDA算法的流程框圖如圖2所示。

圖2 ES-MDA算法流程圖

ES-MDA算法因提出時(shí)間較晚，僅有少量學(xué)者在地下水和油藏模擬中進(jìn)行了少量理論研究，均取得了較為理想的理論性數(shù)值試驗(yàn)結(jié)果[17]，該算法在油藏歷史擬合的適用性還需要進(jìn)一步研究。對(duì)于四維地震數(shù)據(jù)結(jié)合井口生產(chǎn)數(shù)據(jù)的油藏歷史擬合，理論上ESMDA算法契合了多維度時(shí)間序列數(shù)據(jù)的整體同化要求，但由于原始算法中的預(yù)測(cè)步驟中僅進(jìn)行了一次模擬預(yù)測(cè)，無(wú)法滿(mǎn)足地震屬性—油藏流體耦合模擬的要求，并且四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)在時(shí)間和空間兩個(gè)維度上都存在很大差異，故仍需要對(duì)ES-MDA算法進(jìn)行優(yōu)化改進(jìn)。

2 基于四維地震數(shù)據(jù)的歷史擬合

2.1 四維地震監(jiān)測(cè)對(duì)于油藏歷史擬合的作用

油藏歷史擬合現(xiàn)有的方法僅通過(guò)擬合井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)來(lái)修正油藏的模型參數(shù)。誠(chéng)然，油藏狀態(tài)的動(dòng)態(tài)變化過(guò)程，可以反映在井口的生產(chǎn)動(dòng)態(tài)上，但由于井口的覆蓋范圍有限，無(wú)法準(zhǔn)確實(shí)現(xiàn)對(duì)于整個(gè)油藏狀態(tài)的表征。從物理意義上講，由于井口動(dòng)態(tài)稀疏有限的井口動(dòng)態(tài)生產(chǎn)數(shù)據(jù)必然會(huì)造成反向模型求解出現(xiàn)嚴(yán)重的欠定和多解性，所求結(jié)果往往與實(shí)際情況不符。

隨著油田生產(chǎn)，不斷有油氣從油藏中采出，伴隨著油藏內(nèi)流體交換及狀態(tài)的變化，還會(huì)造成油藏壓力、溫度的變化。其中，油藏壓力和流體飽和度等因素的變化，在四維地震監(jiān)測(cè)過(guò)程中，也會(huì)造成地震波響應(yīng)的改變。其中，縱波波阻抗是常見(jiàn)表征油藏狀態(tài)和流體變化的地震屬性指標(biāo)之一。

地震波阻抗反演的地震解釋技術(shù)早在20世紀(jì)70年代就開(kāi)始形成，經(jīng)過(guò)幾十年發(fā)展，已經(jīng)能比較準(zhǔn)確地獲得波阻抗數(shù)據(jù)，在歷史擬合中若成功應(yīng)用，可以增強(qiáng)地質(zhì)模型中流體變化特征的可信度。

伴隨著油氣開(kāi)采的過(guò)程，一方面，溶解氣在儲(chǔ)層內(nèi)脫氣、氣相推動(dòng)油水液相，以及流體的變化造成油藏儲(chǔ)層內(nèi)孔隙中流體密度減?。涣硪环矫?，油藏內(nèi)地層壓力增加，則會(huì)造成孔隙膨脹，油藏巖石密度減小[19]。這兩方面情形都會(huì)降低地震波反射傳播速度，從而減小地震縱波波阻抗，且可以通過(guò)有色反演技術(shù)表現(xiàn)出來(lái)。相反，若在油氣開(kāi)采過(guò)程中，壓力無(wú)法得到及時(shí)補(bǔ)充，造成油藏壓力下降，導(dǎo)致孔隙被壓縮，提高了油藏的巖石和流體整體密度，地震縱波波阻抗數(shù)值就會(huì)上升[20]。

由此可見(jiàn)，在油氣田生產(chǎn)過(guò)程中發(fā)生的油藏狀態(tài)變化，包括油藏壓力、流體狀態(tài)、飽和度等指標(biāo)，都可以通過(guò)四維地震監(jiān)測(cè)技術(shù)來(lái)反映。四維地震監(jiān)測(cè)技術(shù)主要包括四維地震數(shù)據(jù)的采集、處理和解釋?zhuān)⑼ㄟ^(guò)巖石物理學(xué)理論將反射地震波屬性與油藏狀態(tài)屬性聯(lián)系起來(lái)，從而實(shí)現(xiàn)對(duì)整個(gè)油藏范圍全覆蓋的準(zhǔn)動(dòng)態(tài)監(jiān)測(cè)。

2.2 四維地震數(shù)據(jù)的稀疏化處理

首先，按照地質(zhì)—油藏模型耦合的通常處理方法，將處理后的地震縱波波阻抗數(shù)據(jù)體進(jìn)行網(wǎng)格粗化。經(jīng)過(guò)網(wǎng)格粗化以后，地震屬性數(shù)據(jù)體與油藏流體模擬的模型尺度相一致，在獲得油藏流體狀態(tài)（經(jīng)油藏流體數(shù)值模擬后）的基礎(chǔ)上，通過(guò)巖石物理學(xué)模型正演可直接得到相應(yīng)地震波阻抗數(shù)據(jù)。但由于油藏模型網(wǎng)格通常高達(dá)上百萬(wàn)網(wǎng)格量級(jí)，鑒于當(dāng)前電腦硬件計(jì)算能力限制，仍需對(duì)地震波阻抗數(shù)據(jù)進(jìn)一步進(jìn)行稀疏化處理。

在信號(hào)處理領(lǐng)域，稀疏表示（sparse representation）可以對(duì)連續(xù)性的或密集的信號(hào)進(jìn)行稀疏化處理，從而以較小的數(shù)據(jù)體量來(lái)表達(dá)大體量矩陣的關(guān)鍵性信息。對(duì)于信號(hào)s如果只有m（m＜N）個(gè)元素為非零值，其余元素均為零，那么該信號(hào)s被稱(chēng)為m-稀疏信號(hào)[21]。求解稀疏表示的過(guò)程被稱(chēng)為稀疏分解。在通常情況下，一個(gè)信號(hào)可以用K個(gè)信號(hào)原子φk的線性疊加來(lái)表示：

其中，α一般被稱(chēng)為信號(hào)在數(shù)據(jù)字典φ中的表示稀疏系數(shù)（representation coefficients）。字典φ是由信號(hào)原子φk為縱向組成的矩陣，通常被歸一化為方差，即[22]。

對(duì)于四維地震中的地震波阻抗數(shù)據(jù)體，表示的是空間上不同油藏網(wǎng)格中的不同油藏狀態(tài)下的縱波波阻抗數(shù)據(jù)值，雖然不具有信號(hào)那樣時(shí)間上的連續(xù)性，但考慮到一定范圍內(nèi)地質(zhì)和油藏屬性的連續(xù)性，可以仿照信號(hào)數(shù)組，受信號(hào)處理領(lǐng)域中的稀疏表示啟發(fā)，通過(guò)設(shè)計(jì)一個(gè)空間上的單位稀疏矩陣作為稀疏系數(shù)，來(lái)實(shí)現(xiàn)對(duì)超大型四維地震縱波波阻抗數(shù)據(jù)的稀疏化處理。四維地震監(jiān)測(cè)中，分別在t1，t2，t3，...，tN時(shí)刻進(jìn)行三維地震采集，經(jīng)地震數(shù)據(jù)處理與解釋?zhuān)约熬W(wǎng)格粗化后，每個(gè)時(shí)間點(diǎn)的縱波地震波阻抗數(shù)據(jù)為Ip_large（t），其維度與油藏網(wǎng)格的維度一致。對(duì)每個(gè)時(shí)間點(diǎn)的可采用相同的稀疏系數(shù)矩陣。Ip_4D為稀疏處理后的四維地震波阻抗矩陣，其過(guò)程可通過(guò)式（21）表示：

2.3 歷史擬合ES-MDA算法優(yōu)化

原始ES-MDA算法在油藏歷史擬合中由于僅吸收了井口動(dòng)態(tài)生產(chǎn)數(shù)據(jù)，其預(yù)測(cè)步驟為運(yùn)行全時(shí)段的油藏流體模擬程序，計(jì)算出油田生產(chǎn)歷史的完整井口生產(chǎn)動(dòng)態(tài)預(yù)測(cè)數(shù)據(jù)。由于需要耦合油藏流體數(shù)值模擬與地震屬性正向模擬，原始ES-MDA算法的預(yù)測(cè)步驟僅做一次油藏流體模擬，無(wú)法滿(mǎn)足要求，故需要對(duì)ES-MDA算法進(jìn)行優(yōu)化改進(jìn)。

首先，建立起油藏流體和地震波阻抗的耦合模擬與油藏地質(zhì)參數(shù)的數(shù)學(xué)關(guān)系。油藏流體數(shù)值模擬不僅需要輸出井口動(dòng)態(tài)生產(chǎn)數(shù)據(jù)，為了進(jìn)行地震屬性的正演模擬，也需要輸出油藏狀態(tài)函數(shù)。油藏模擬器的計(jì)算過(guò)程可以通過(guò)式（22）表示。

地震屬性的正演模擬可根據(jù)前面介紹的地震縱波波阻抗的計(jì)算方法完成，其計(jì)算可以通過(guò)式（23）簡(jiǎn)要表示。

式中I（x）——地震縱波波阻抗正演模擬程序；

Ip——模擬計(jì)算得到的地震縱波波阻抗數(shù)據(jù)。

將油藏流體模擬器輸出的油藏狀態(tài)參數(shù)作為地震縱波波阻抗正演模擬的輸入數(shù)據(jù)，再將計(jì)算得到的地震縱波波阻抗數(shù)據(jù)與油藏流體模擬器輸出的井口動(dòng)態(tài)生產(chǎn)數(shù)據(jù)相結(jié)合，即可實(shí)現(xiàn)油藏流體—地震屬性的耦合模擬，其數(shù)學(xué)模型可以通過(guò)式（24）表示。

GI(m)——油藏流體—地震屬性耦合模擬程序，其程序架構(gòu)圖如圖3所示。

圖3 油藏流體—地震屬性耦合模擬程序架構(gòu)圖

使用地震縱波波阻抗數(shù)據(jù)和井口生產(chǎn)數(shù)據(jù)進(jìn)行歷史擬合，還需要對(duì)這兩類(lèi)動(dòng)態(tài)數(shù)據(jù)在空間和時(shí)間兩個(gè)維度上進(jìn)行匹配。一方面，由于四維地震數(shù)據(jù)采集的時(shí)間間隔往往是1年以上，遠(yuǎn)遠(yuǎn)高于井口生產(chǎn)數(shù)據(jù)采集的時(shí)間間隔；另一方面，地震縱波波阻抗數(shù)據(jù)的空間密度要遠(yuǎn)遠(yuǎn)高于井口的生產(chǎn)數(shù)據(jù)。因此，需要對(duì)四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)進(jìn)行數(shù)據(jù)匹配。

基于減小數(shù)據(jù)體量、節(jié)省計(jì)算時(shí)間、保證運(yùn)算數(shù)據(jù)穩(wěn)定性這3個(gè)原則，數(shù)據(jù)匹配工作需要實(shí)現(xiàn)以下3個(gè)目標(biāo)：（1）將地震波阻抗的空間維度降低。通過(guò)在空間上設(shè)計(jì)合理的稀疏采樣的方法，將地震縱波波阻抗數(shù)據(jù)的空間體量降低2～3個(gè)數(shù)量級(jí)，從而達(dá)到減小數(shù)據(jù)體量的效果。（2）保持井口生產(chǎn)數(shù)據(jù)的時(shí)間維度。將地震縱波波阻抗數(shù)據(jù)在時(shí)間維度上插入到井口生產(chǎn)數(shù)據(jù)上，以保證ES-MDA算法所吸收的數(shù)據(jù)有足夠反映油藏動(dòng)態(tài)變化的時(shí)間密度，以及保證后驗(yàn)概率更新的可靠性。（3）不同類(lèi)型實(shí)際測(cè)量數(shù)據(jù)的數(shù)值大小不同，數(shù)據(jù)吸收算法會(huì)把數(shù)值大的數(shù)據(jù)給予更高的影響權(quán)重，且如果數(shù)值的數(shù)量級(jí)差距過(guò)大，也會(huì)造成算法的不穩(wěn)定。為此，需要設(shè)計(jì)實(shí)測(cè)數(shù)據(jù)矩陣，使得不同類(lèi)型的實(shí)際測(cè)量數(shù)據(jù)的數(shù)量級(jí)一致，如式（26）所示。

式中Dobs—— 匹配后的包含井口生產(chǎn)動(dòng)態(tài)和四維地震縱波波阻抗的實(shí)際測(cè)量數(shù)據(jù)；

標(biāo)量fd，fI—— 分別是井口生產(chǎn)數(shù)據(jù)和稀疏化處理后的四維地震縱波波阻抗數(shù)據(jù)的匹配系數(shù)；

dobs—— 包含原油產(chǎn)量、產(chǎn)水量、注水井井底流壓等井口生產(chǎn)動(dòng)態(tài)的實(shí)際測(cè)量數(shù)據(jù)。

根據(jù)式（26），數(shù)據(jù)匹配的關(guān)鍵在于匹配系數(shù)設(shè)計(jì)。良好的匹配系數(shù)設(shè)計(jì)會(huì)使得需要擬合的實(shí)際測(cè)量數(shù)據(jù)dobs中的每個(gè)非零元素?cái)?shù)值的數(shù)量級(jí)相當(dāng)，從而實(shí)現(xiàn)在歷史擬合過(guò)程中，每種類(lèi)型的實(shí)際測(cè)量數(shù)據(jù)都能得到充分利用和擬合，也保證了程序的數(shù)值穩(wěn)定性。之后，對(duì)ES-MDA算法中的調(diào)整實(shí)際測(cè)量數(shù)據(jù)步驟進(jìn)行優(yōu)化改進(jìn)，結(jié)合四維地震與井口動(dòng)態(tài)的數(shù)據(jù)匹配方法，根據(jù)式（27）對(duì)匹配后的實(shí)際測(cè)量數(shù)據(jù)加上白色噪聲。

式中Duc—— 加誤差后的包含井口生產(chǎn)動(dòng)態(tài)和四維地震波阻抗的數(shù)據(jù)。

通過(guò)上述對(duì)預(yù)測(cè)步驟和調(diào)整實(shí)測(cè)數(shù)據(jù)步驟的改進(jìn)，形成了ES-MDA算法的改進(jìn)型，其算法流程圖如圖4所示。

圖4 ES-MDA算法改進(jìn)型的流程框圖

2.4 新方法流程

綜合上文所述，筆者提出了一套基于四維地震波阻抗數(shù)據(jù)與生產(chǎn)數(shù)據(jù)的油藏自動(dòng)歷史擬合新方法體系，主要步驟如下：

（1）儲(chǔ)層隨機(jī)建模形成初始油藏地質(zhì)模型集合。利用測(cè)井?dāng)?shù)據(jù)、巖心數(shù)據(jù)、地震數(shù)據(jù)等所有可利用的資料，根據(jù)儲(chǔ)層類(lèi)型特點(diǎn)，選取合理的儲(chǔ)層隨機(jī)建模方法，建立一整套初始油藏地質(zhì)模型集合；

（2）井口生產(chǎn)數(shù)據(jù)預(yù)處理。分析所有的油井生產(chǎn)歷史數(shù)據(jù)，去除掉明顯異常值，盡量減少產(chǎn)量的高頻波動(dòng)，選取適當(dāng)?shù)臅r(shí)間間隔。

（3）地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)匹配。原始地震數(shù)據(jù)體的空間密度往往非常高，而時(shí)間密度卻遠(yuǎn)低于井口生產(chǎn)數(shù)據(jù)。按照前面介紹的方法，在空間維度上降低地震數(shù)據(jù)體量，時(shí)間維度上保持井口生產(chǎn)數(shù)據(jù)的密度，并形成實(shí)際測(cè)量數(shù)據(jù)向量dobs。

（4）運(yùn)行ES-MDA自動(dòng)歷史擬合程序。設(shè)定好數(shù)據(jù)吸收次數(shù)Na和加權(quán)吸收α，輸入初始地質(zhì)模型集合及數(shù)據(jù)配合后的實(shí)際測(cè)量數(shù)據(jù)dobs，運(yùn)行ESMDA自動(dòng)歷史擬合程序，不斷調(diào)用油藏流體模擬程序，修正得到后驗(yàn)地質(zhì)模型集合。

（5）地質(zhì)模型驗(yàn)證與預(yù)測(cè)數(shù)據(jù)分析。ES-MDA油藏自動(dòng)歷史擬合得到后驗(yàn)地質(zhì)模型，需要仔細(xì)對(duì)照地質(zhì)模型是否符合所在儲(chǔ)層類(lèi)型成藏過(guò)程的地質(zhì)形態(tài)，再對(duì)后驗(yàn)地質(zhì)模型運(yùn)行油藏流體模擬程序，查看預(yù)測(cè)的生產(chǎn)數(shù)據(jù)是否與實(shí)際測(cè)量的動(dòng)態(tài)數(shù)據(jù)相符合。

如果歷史擬合結(jié)果不符合預(yù)期，則需要重回第二步和第三步，重新進(jìn)行數(shù)據(jù)匹配和ES-MDA自動(dòng)歷史擬合程序的參數(shù)設(shè)計(jì)，直到得到較好的結(jié)果?；谒木S地震數(shù)據(jù)的油藏自動(dòng)歷史擬合程序架構(gòu)圖見(jiàn)圖5。

圖5 基于四維地震數(shù)據(jù)的油藏自動(dòng)歷史擬合程序架構(gòu)圖

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

3.1 歷史擬合案例設(shè)置

本文研究區(qū)以布倫特油田某海相砂巖油藏為原型，建立油藏模型作為真實(shí)油藏地質(zhì)情況的參照。該油藏為簡(jiǎn)單背斜構(gòu)造，邊界由斷層和地質(zhì)層位構(gòu)成，內(nèi)部無(wú)明顯斷層，油藏邊界為密閉斷層，不含活躍邊底水，缺少外界能量供給。對(duì)于海相沉積環(huán)境，通常在海平面上升時(shí)期進(jìn)行沉積，在海平面下降時(shí)則沒(méi)有沉積物沉降，但海的作用影響到海岸，從而會(huì)造成沉積邊界形成[23]。油藏具有較強(qiáng)的非均質(zhì)性，存在未知的流體高滲通道。

參照油藏模型由40×120×20個(gè)正交立方網(wǎng)格組成，每個(gè)網(wǎng)格尺寸是75m×75m×4m，其中78720個(gè)網(wǎng)格為有效儲(chǔ)層網(wǎng)格，其余為非儲(chǔ)層死網(wǎng)格。其中，孔隙呈高斯分布；海相砂巖儲(chǔ)層滲透率基本呈對(duì)數(shù)高斯分布（圖6）。

圖6 對(duì)數(shù)滲透率參照模型與頻率分布直方圖

通過(guò)序貫高斯模擬建立100套孔隙度和對(duì)數(shù)滲透率的初始地質(zhì)模型。對(duì)數(shù)滲透率初始模型的平均值和其中4個(gè)初始模型如圖7所示。

圖7 對(duì)數(shù)滲透率初始地質(zhì)模型

由于序貫高斯模擬使用的是隨機(jī)路徑條件下根據(jù)變差函數(shù)的克里金的估算，故初始地質(zhì)模型的不確定性主要表現(xiàn)在估算點(diǎn)的空間位置上。初始地質(zhì)模型的平均地質(zhì)參數(shù)在空間分布上較為均勻，難以判斷油藏的地質(zhì)特征；而對(duì)單獨(dú)的初始地質(zhì)模型則有明顯的地質(zhì)屬性空間分布規(guī)律。這說(shuō)明未經(jīng)歷史擬合的初始地質(zhì)模型在空間上存在很大的不確定性。為減小地質(zhì)模型的不確定性，需要納入油藏動(dòng)態(tài)數(shù)據(jù)，在此基礎(chǔ)上進(jìn)行貝葉斯后驗(yàn)概率估計(jì)，從而提高地質(zhì)模型的置信度。

為保持油藏壓力，在油田生產(chǎn)初始階段即采用注水采油方式生產(chǎn)，有5口注水井和5口采油井，以反五點(diǎn)法方式注水，P1、P2、P3、P4和P5為采油井，Ip1、Ip2、Ip3、Ip4和Ip5為注水井。已知該油藏18年的井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)，如生產(chǎn)井產(chǎn)油量、產(chǎn)水量，以及注水井的井底流壓。對(duì)該油田海底設(shè)置永久地震檢波器，每2年進(jìn)行一次四維地震采集，經(jīng)過(guò)處理與解釋?zhuān)@得縱波波阻抗數(shù)據(jù)體。

下面通過(guò)只使用井口生產(chǎn)數(shù)據(jù)和綜合使用四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)進(jìn)行油藏自動(dòng)歷史的研究來(lái)檢驗(yàn)基于四維地震數(shù)據(jù)的油藏歷史擬合方法的可靠性，探討四維地震數(shù)據(jù)對(duì)于修正油藏模型的重要意義。

3.2 方法1：井口生產(chǎn)數(shù)據(jù)的油藏歷史擬合

這里僅對(duì)井口生產(chǎn)數(shù)據(jù)進(jìn)行擬合，用以驗(yàn)證ESMDA歷史擬合算法程序的可靠性。擬合的井口生產(chǎn)數(shù)據(jù)有3種：（1）油井產(chǎn)油量；（2）油井產(chǎn)水量；（3）注水井井底流壓。經(jīng)大量試驗(yàn)確定[24]，ES-MDA歷史擬合中，設(shè)置數(shù)據(jù)吸收次數(shù)Na=4，每次數(shù)據(jù)吸收的α分別為α1=α2=α3=α4=4。

井口生產(chǎn)數(shù)據(jù)歷史擬合結(jié)果如圖8所示。圖中綠色細(xì)實(shí)線為初始地質(zhì)模型通過(guò)油藏流體模擬預(yù)測(cè)的油田生產(chǎn)數(shù)據(jù)，包括油田日產(chǎn)油量、油井P1的日產(chǎn)油量和日產(chǎn)水量、注水井I1的井底流壓。圖中，紅色粗虛線為實(shí)際測(cè)量的生產(chǎn)數(shù)據(jù)，黑色細(xì)實(shí)線為經(jīng)過(guò)歷史擬合修正地質(zhì)模型后，在油藏流體模擬器上預(yù)測(cè)得到的生產(chǎn)數(shù)據(jù)結(jié)果。

圖8 方法1：井口生產(chǎn)數(shù)據(jù)歷史擬合結(jié)果

圖9為修正后的對(duì)數(shù)滲透率模型?？梢钥闯?，對(duì)于生產(chǎn)井的產(chǎn)油量和產(chǎn)水量擬合結(jié)果良好，黑色細(xì)實(shí)線基本圍繞在紅色粗虛線附近，且分布遠(yuǎn)比綠色細(xì)實(shí)線更集中，說(shuō)明經(jīng)ES-MDA算法修正地質(zhì)模型后，對(duì)后驗(yàn)地質(zhì)模型進(jìn)行油藏流體模擬預(yù)測(cè)，基本上擬合了實(shí)際測(cè)量的井口生產(chǎn)數(shù)據(jù)。

圖9 方法1：修正后的對(duì)數(shù)滲透率模型

3.3 方法2：四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)的油藏歷史擬合

在前面針對(duì)井口生產(chǎn)數(shù)據(jù)的油藏歷史擬合基礎(chǔ)上，將通過(guò)智能地震處理解釋技術(shù)[25]獲得的四維地震的縱波波阻抗數(shù)據(jù)也納入歷史擬合，真正實(shí)現(xiàn)了四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)的油藏自動(dòng)歷史擬合研究[16]。

首先，分析各類(lèi)型動(dòng)態(tài)數(shù)據(jù)進(jìn)行數(shù)據(jù)匹配，設(shè)置實(shí)際測(cè)量數(shù)據(jù)匹配矩陣，如式（28）所示。將油井產(chǎn)油量fQ、油井產(chǎn)水量fW、注入井井底流壓fB、地震縱波波阻抗的匹配系數(shù)fI分別設(shè)為1、10、0.005和0.1，從而得到數(shù)據(jù)匹配后的實(shí)際測(cè)量數(shù)據(jù)向量Dobs。

運(yùn)行基于四維地震數(shù)據(jù)的油藏自動(dòng)歷史擬合程序，井口動(dòng)態(tài)生產(chǎn)數(shù)據(jù)的擬合結(jié)果如圖10所示。對(duì)比圖8，圖10中的各項(xiàng)生產(chǎn)數(shù)據(jù)擬合結(jié)果要明顯更好。圖中黑色細(xì)實(shí)線不僅都圍繞在紅色粗虛線附近，而且其分布更集中于實(shí)際測(cè)量數(shù)值。說(shuō)明同時(shí)擬合四維地震數(shù)據(jù)與井口生產(chǎn)數(shù)據(jù)，要比僅對(duì)井口生產(chǎn)數(shù)據(jù)進(jìn)行歷史擬合數(shù)值收斂性增強(qiáng)很多，對(duì)后驗(yàn)地質(zhì)模型進(jìn)行油藏流體數(shù)值模擬的預(yù)測(cè)結(jié)果更接近于真實(shí)測(cè)量數(shù)據(jù)，間接表明了后驗(yàn)地質(zhì)模型的置信度的提高。修正后的對(duì)數(shù)滲透率模型如圖11所示。

圖10 方法2：井口生產(chǎn)數(shù)據(jù)歷史擬合結(jié)果

圖11 方法2：修正后的對(duì)數(shù)滲透率模型

將后驗(yàn)地質(zhì)模型與方法1的結(jié)果對(duì)比，不難發(fā)現(xiàn)方法2的后驗(yàn)地質(zhì)模型集合的平均地質(zhì)參數(shù)分布擁有更明顯的特征，更接近于真實(shí)的油藏情況。且單個(gè)集合中的子集地質(zhì)模型相互之間的差異也減小更多。說(shuō)明加入地震數(shù)據(jù)進(jìn)行歷史擬合的后驗(yàn)地質(zhì)模型的置信度更高的同時(shí)，也證實(shí)4次數(shù)據(jù)吸收的次數(shù)已足夠充分。增加地震數(shù)據(jù)作為修正條件之后，ES-MDA算法的數(shù)據(jù)吸收次數(shù)明顯減少。

3.4 不同方法后驗(yàn)地質(zhì)模型置信度定量分析

為了定量分析對(duì)比不同油藏歷史擬合方法的后驗(yàn)地質(zhì)模型置信度，設(shè)計(jì)了兩個(gè)統(tǒng)計(jì)學(xué)指標(biāo)，用于對(duì)比這兩種方法的后驗(yàn)地質(zhì)模型進(jìn)行油藏流體數(shù)值模擬預(yù)測(cè)結(jié)果的準(zhǔn)確性[26]。

后驗(yàn)地質(zhì)模型集合的平均地質(zhì)參數(shù)模型的油藏流體數(shù)值模擬的預(yù)測(cè)相對(duì)誤差e，定義如下：

式中，預(yù)測(cè)相對(duì)誤差可以通過(guò)e的下角標(biāo)來(lái)表示。油田生產(chǎn)的最后一天，對(duì)于油井產(chǎn)油量、產(chǎn)水量的預(yù)測(cè)相對(duì)誤差，以eoPi和ewPi表示；對(duì)于注水井井底流壓的預(yù)測(cè)相對(duì)誤差以ebhpIi來(lái)表示。其中下角標(biāo)中的i表示模型編號(hào)，P和I則分別表示采油井與注水井。

（2）后驗(yàn)地質(zhì)模型集合的平均地質(zhì)參數(shù)模型的油藏流體數(shù)值模擬的預(yù)測(cè)方差，定義如下：

同樣，式中預(yù)測(cè)方差可以通過(guò)σ2的下角標(biāo)來(lái)表示，下角標(biāo)的標(biāo)準(zhǔn)準(zhǔn)則與預(yù)測(cè)相對(duì)誤差e相同。

表1和表2分表列出對(duì)于預(yù)測(cè)相對(duì)誤差e和預(yù)測(cè)方差σ2的值。對(duì)比表中兩種方法的預(yù)測(cè)誤差和方差結(jié)果可知，方法2得到的后驗(yàn)地質(zhì)模型運(yùn)行油藏流體模擬程序之后的各項(xiàng)預(yù)測(cè)指標(biāo)，明顯好于方法1的預(yù)測(cè)，這也驗(yàn)證了四維地震數(shù)據(jù)對(duì)于油藏歷史擬合結(jié)果的積極意義。

表1 兩個(gè)方法的預(yù)測(cè)相對(duì)誤差指標(biāo)統(tǒng)計(jì)表 （單位：%）

表2 兩個(gè)方法的預(yù)測(cè)方差指標(biāo)統(tǒng)計(jì)表

4 結(jié)論

為了充分利用四維地震監(jiān)測(cè)技術(shù)提供的覆蓋全油藏的狀態(tài)信息數(shù)據(jù)，解決現(xiàn)有油藏自動(dòng)歷史擬合方法僅能處理井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)造成的數(shù)據(jù)利用程度低、擬合效果差、油藏模型更新可靠性低等問(wèn)題，本文基于貝葉斯理論方法，優(yōu)化改進(jìn)了ES-MDA算法，提出一種基于四維地質(zhì)歷史擬合的新方法，大幅提高了油藏自動(dòng)歷史擬合的質(zhì)量，并在某海相砂巖油藏的歷史擬合應(yīng)用中取得了良好效果。

本文通過(guò)理論推導(dǎo)和應(yīng)用研究，得到以下認(rèn)識(shí)：

（1）針對(duì)四維地震數(shù)據(jù)的空間密度要遠(yuǎn)高于井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)，而時(shí)間密度遠(yuǎn)低于井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)的問(wèn)題，提出了四維地震數(shù)據(jù)的粗化與稀疏化處理方法，并通過(guò)設(shè)計(jì)匹配系數(shù)形成實(shí)際測(cè)試數(shù)據(jù)矩陣，實(shí)現(xiàn)了四維地震數(shù)據(jù)和井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)的匹配。

（2）通過(guò)建立起油藏流體和地震波阻抗的耦合模擬與油藏地質(zhì)參數(shù)的數(shù)學(xué)關(guān)系，形成油藏流體—地震屬性耦合模擬程序，并結(jié)合數(shù)據(jù)匹配方法，改進(jìn)優(yōu)化ES-MDA算法，建立了基于四維地震數(shù)據(jù)的歷史擬合新方法。

（3）同時(shí)擬合四維地震數(shù)據(jù)與井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)，比僅對(duì)井口生產(chǎn)數(shù)據(jù)進(jìn)行歷史擬合所得到的數(shù)據(jù)擬合結(jié)果更好，后驗(yàn)地質(zhì)模型的置信度更高，且增加了四維地震數(shù)據(jù)作為修正條件之后，也能明顯減少ES-MDA算法的數(shù)據(jù)吸收次數(shù)，體現(xiàn)了四維地震監(jiān)測(cè)的價(jià)值。四維地震與井口生產(chǎn)動(dòng)態(tài)數(shù)據(jù)的油藏歷史擬合具有實(shí)際意義和應(yīng)用前景。