王澤龍 劉先貴 李志勇 張 輝 張傳進

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2 EnKF算法

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

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

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

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

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

下角標k—— 時間序列中的第k項，即第k時間步（k為整數(shù)）。

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

上角標f——油藏模型預測的數(shù)據(jù)；

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

至此，完成了在k時刻的后驗模型集合的估算。通過不斷重復此過程，對比不同時刻的實際測量數(shù)據(jù)，并不斷更新油藏地質(zhì)模型集合。EnKF的算法流程如圖1所示。

圖1 EnKF算法流程圖

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

1.3 ES-MDA算法

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

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

圖2 ES-MDA算法流程圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 ES-MDA算法改進型的流程框圖

2.4 新方法流程

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

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

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

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

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

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

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

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

3 應用實例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

同樣，式中預測方差可以通過σ2的下角標來表示，下角標的標準準則與預測相對誤差e相同。

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

表1 兩個方法的預測相對誤差指標統(tǒng)計表 （單位：%）

表2 兩個方法的預測方差指標統(tǒng)計表

4 結(jié)論

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

本文通過理論推導和應用研究，得到以下認識：

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

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

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