Seifi A，Kazemzadeh M B，Mohammadi H

（伊朗Amirkabir科技大學(xué)工業(yè)工程系）

0 引言

一般采用數(shù)值模擬方法預(yù)測油氣藏的累計(jì)產(chǎn)量。數(shù)值模擬屬于密集型計(jì)算過程，計(jì)算量大，采用元建模方法并建立合適的元模型，可以減小計(jì)算量、提高計(jì)算效率。元模型是一種代理數(shù)學(xué)方程，主要用于近似表示復(fù)雜體系內(nèi)的輸入-輸出關(guān)系[1]。元建?；静襟E[1-6]包括：①定義輸入變量及相應(yīng)的輸出響應(yīng)；②進(jìn)行實(shí)驗(yàn)設(shè)計(jì)（DoE），選擇樣本點(diǎn)，并確定各樣本點(diǎn)的響應(yīng)；③構(gòu)建元模型，建立輸入變量與輸出響應(yīng)之間的關(guān)系；④通過試驗(yàn)點(diǎn)驗(yàn)證元模型的準(zhǔn)確性。

本文以伊朗南部某大型裂縫性氣藏為例，基于元建模方法對該氣藏在18年規(guī)劃周期內(nèi)的累計(jì)產(chǎn)量進(jìn)行預(yù)測。建立3種不同的元模型（傳統(tǒng)二次模型、乘性模型及徑向基函數(shù)（RBF）模型），并利用基于數(shù)值模擬的氣藏累計(jì)產(chǎn)量預(yù)測結(jié)果對3種模型進(jìn)行驗(yàn)證和對比。

1 氣藏基本參數(shù)估計(jì)及數(shù)值模擬

為了對元模型的準(zhǔn)確性進(jìn)行驗(yàn)證和對比，在建立元模型之前預(yù)先估算了氣藏的基本參數(shù)（孔隙度、滲透率及含水飽和度）并對氣藏進(jìn)行了數(shù)值模擬。

本文研究對象為伊朗南部某大型裂縫性氣藏。對該氣藏進(jìn)行網(wǎng)格劃分，x、y、z方向的網(wǎng)格數(shù)分別為45、37和95，總網(wǎng)格數(shù)為158 175。對探井的FMI測井?dāng)?shù)據(jù)進(jìn)行了解釋和分析，計(jì)算了裂縫參數(shù)（方位角、傾角及密度），考慮裂縫參數(shù)的影響，采用序貫指示模擬（SIS）法和序貫高斯模擬（SGS）法[7]對各網(wǎng)格的孔隙度和滲透率進(jìn)行了估算。

為了更真實(shí)地反應(yīng)本文研究氣藏的孔隙空間分布特征，首先用裂縫面模型對孔隙度模型進(jìn)行了改進(jìn)，即根據(jù)裂縫面模型求取孔隙度分布的變差函數(shù)，然后選用相關(guān)的變差函數(shù)估算各網(wǎng)格的孔隙度。采用序貫高斯模擬（SGS）法，共獲得了孔隙度分布的30個(gè)模擬結(jié)果，并選擇了擬合效果最佳的模擬結(jié)果。

獲取了研究區(qū)13口井的巖心滲透率數(shù)據(jù)?；?個(gè)水平方向（x，y方向）和1個(gè)垂直方向（z方向）的裂縫參數(shù)計(jì)算了滲透率修正系數(shù)，并生成了 x，y和z方向的變差函數(shù)。采用序貫指示模擬（SIS）法，對裂縫參數(shù)分布進(jìn)行了模擬，進(jìn)而得到考慮裂縫影響的滲透率模型。

對滲透率和孔隙度模型進(jìn)行了三維可視化處理（見圖 1），以直觀顯示儲集層的連通性、阻流帶發(fā)育情況及高孔高滲帶的分布。

圖1 x，y和z方向的滲透率模型及孔隙度模型

基于已估算的氣藏基本參數(shù)，建立了該氣藏的Eclipse數(shù)值模型，以對元模型進(jìn)行驗(yàn)證和對比。采用Eclipse軟件在5臺并行PC機(jī)（CPU為雙核Intel 3.4 GHz）上計(jì)算了1口井生產(chǎn)18年的累計(jì)產(chǎn)量，耗時(shí)超過15 h。如果對所有158 175個(gè)網(wǎng)格進(jìn)行計(jì)算，將會耗費(fèi)很長時(shí)間，而采用元建模的方法可以提高計(jì)算效率。

2 元建模

2.1 輸入變量和輸出響應(yīng)

本文假定與氣藏累計(jì)產(chǎn)量相關(guān)的參數(shù)為 x方向的滲透率（Kx）、y方向的滲透率（Ky）、z方向的滲透率（Kz）、孔隙度（φ）和含水飽和度（Sw），則氣藏內(nèi)某井18年內(nèi)的累計(jì)產(chǎn)量可表示為：

（3）式中，Y為輸出響應(yīng)，X為輸入變量，F(xiàn)(X)為近似響應(yīng)函數(shù)，ε為誤差。本文建立了3種不同的元

因此，元建模的輸入變量為Kx，Ky，Kz，φ和Sw，輸出響應(yīng)為Qacc。

2.2 實(shí)驗(yàn)設(shè)計(jì)

實(shí)驗(yàn)設(shè)計(jì)方法可分為兩類：依賴于模型的設(shè)計(jì)和獨(dú)立于模型的設(shè)計(jì)，一般情況下獨(dú)立于模型的設(shè)計(jì)在確定性計(jì)算機(jī)實(shí)驗(yàn)中的應(yīng)用效果更佳[8]。在獨(dú)立于模型的實(shí)驗(yàn)設(shè)計(jì)中，最常用的最優(yōu)性判斷準(zhǔn)則包括最大熵、均方誤差、極大極小和極小極大 4種。本文采用最大熵設(shè)計(jì)方法[9-10]，最大熵設(shè)計(jì)是一種序貫方法，通過搜索參數(shù)空間尋找設(shè)計(jì)點(diǎn)，進(jìn)而使相應(yīng)設(shè)計(jì)點(diǎn)的實(shí)驗(yàn)提取信息（熵）最大化，可采用遺傳算法[11]求解。

進(jìn)行氣藏元建模實(shí)驗(yàn)設(shè)計(jì)時(shí)，應(yīng)考慮儲集層特征[11]。由于元建模的目的是尋找累計(jì)產(chǎn)量最高的點(diǎn)，在對元模型進(jìn)行驗(yàn)證時(shí)只考慮高潛力點(diǎn)，不考慮元模型在低潛力點(diǎn)的擬合效果。本文研究的氣藏復(fù)雜且網(wǎng)格數(shù)量大，因此選用HIP經(jīng)驗(yàn)函數(shù)（儲量計(jì)算公式）排除研究空間內(nèi)的低潛力點(diǎn)。HIP經(jīng)驗(yàn)函數(shù)的表達(dá)式為：

利用HIP經(jīng)驗(yàn)函數(shù)可以生成氣藏網(wǎng)格空間的產(chǎn)能等值線（見圖 2），以顯示不同位置的產(chǎn)氣潛力。初始設(shè)計(jì)點(diǎn)的位置選在等值線的峰值處。基于等值線圖中的信息，選定了 200個(gè)高潛力點(diǎn)作為輸入點(diǎn)?？紤]到元模型內(nèi)參數(shù)的數(shù)量及相互關(guān)系，為了滿足估算模型參數(shù)所需的自由度，至少應(yīng)開展25次實(shí)驗(yàn)，基于最大熵準(zhǔn)則從200個(gè)高潛力點(diǎn)中選取25個(gè)點(diǎn)作為設(shè)計(jì)點(diǎn)[11]。

圖2 由HIP經(jīng)驗(yàn)函數(shù)生成的儲量等值線圖

2.3 元模型建立

元建模最關(guān)鍵的步驟是建立輸入變量與輸出響應(yīng)之間的關(guān)系，即建立元模型。元模型的一般形式為：模型，包括傳統(tǒng)二次模型、乘性模型和徑向基函數(shù)（RBF）模型，并對3種模型進(jìn)行驗(yàn)證和對比。

2.3.1 傳統(tǒng)二次模型

傳統(tǒng)二次模型的優(yōu)點(diǎn)為：可以分別確認(rèn)各因素的影響程度，有助于從模型中剔除影響較小的因素，繼而生成更平滑的模型。用于計(jì)算氣藏累計(jì)產(chǎn)量的二次模型可表示為：

上述二次模型中共包含21個(gè)系數(shù)（β0～β20），部分系數(shù)不具有明顯的統(tǒng)計(jì)意義，因而可以將相應(yīng)的項(xiàng)從模型中去除，從而簡化模型。為了確定各系數(shù)的重要性，定義原假設(shè)為 βi=0（i=0，1，…，20），備擇假設(shè)為βi≠0，采用t檢驗(yàn)法進(jìn)行假設(shè)檢驗(yàn)，顯著性水平α取為 0.1。檢驗(yàn)結(jié)果表明：僅有 7個(gè)系數(shù)（β0，β1，β3，β4，β5，β6，β7）對輸出響應(yīng)有顯著影響。因此，二次模型可簡化為：

2.3.2 乘性模型

以乘性模型為基礎(chǔ)，可構(gòu)建一種近似表示輸入變量與輸出響應(yīng)間關(guān)系的替代模型，該近似模型強(qiáng)調(diào)氣藏靜態(tài)參數(shù)間相互作用對氣藏產(chǎn)量的影響。乘性模型中各因素的影響具有倍增性，例如，在任何給定點(diǎn)，如果含水飽和度高，即使?jié)B透率和孔隙度很高，氣藏產(chǎn)量仍然很低。因此，為了確保每口井均有較高的產(chǎn)量，泄氣半徑內(nèi)的所有氣藏參數(shù)均需有利于增產(chǎn)。構(gòu)建乘性模型的思路與HIP經(jīng)驗(yàn)函數(shù)一致。用于計(jì)算氣藏累計(jì)產(chǎn)量的乘性模型可表示為：

上述乘性模型中包含各氣藏參數(shù)的方差（用 V(·)表示），以表征氣藏物性垂向變化以及因?qū)Ω鲄?shù)屬性模型進(jìn)行網(wǎng)格粗化所導(dǎo)致的不確定性。模型中包含 11個(gè)系數(shù)（γ0～γ10），基于25個(gè)設(shè)計(jì)點(diǎn)的氣藏?cái)?shù)值模擬數(shù)據(jù)，采用非線性回歸分析法，對各系數(shù)進(jìn)行了估算。擬合后的乘性模型為：

2.3.3 RBF模型

以徑向基函數(shù)（RBF）模型表示輸入變量與輸出響應(yīng)間的關(guān)系必須具備 2個(gè)主要條件：目標(biāo)函數(shù)為噪聲函數(shù)，且噪聲因子的影響隨著徑向點(diǎn)與參考點(diǎn)間距離的增加而嚴(yán)格遞減或遞增。RBF模型中包含1個(gè)連接界面，以銜接輸入變量與目標(biāo)函數(shù)。用于計(jì)算氣藏累計(jì)產(chǎn)量的RBF模型可表示為：

其中

（8）式中，hj(X)為徑向基函數(shù)，表示徑向點(diǎn)（Aj）處屬性向量Xj的非線性界面函數(shù)估計(jì)值。確定模型形式的關(guān)鍵在于選定合適的徑向基函數(shù)，常用的標(biāo)準(zhǔn)徑向基函數(shù)包括高斯RBF、薄板樣條RBF、多二次RBF、逆多二次RBF等[12]。這些徑向基函數(shù)各有優(yōu)點(diǎn)，但高斯RBF的使用頻率更高[12-15]。對本文的研究對象而言，隨著與參考點(diǎn)（Ar）間距離的增大，徑向點(diǎn)對目標(biāo)函數(shù)的影響減弱，因此，高斯RBF是較好的選擇，其表達(dá)式為：

（9）式中，c是一個(gè)標(biāo)量，可以通過獲取更為精確的近似值的方式進(jìn)行優(yōu)化，且通常與模型中的一些方差相關(guān)?！琗r-Xj‖2為Xr-Xj的歐氏范數(shù)，Xr為參考點(diǎn)處屬性微量，Xr與Xj均為包含475個(gè)元素的長向量（5個(gè)氣藏參數(shù)、95個(gè)層位）。

在建立本文所研究氣藏的RBF模型時(shí)，某個(gè)設(shè)計(jì)點(diǎn)被視為參考點(diǎn) Ar，鄰近網(wǎng)格的位置可以被視為參考點(diǎn)周圍的徑向點(diǎn)Aj，圖 3顯示了給定的參考點(diǎn)周圍的徑向點(diǎn)分布情況。用于估算 RBF模型中各系數(shù)（λ0，λ1，…，λk）的數(shù)據(jù)集由 25個(gè)設(shè)計(jì)點(diǎn)及每個(gè)設(shè)計(jì)點(diǎn)周圍的24個(gè)徑向點(diǎn)的氣藏參數(shù)組成，基于這些數(shù)據(jù)，采用線性回歸分析方法即可估算模型中各系數(shù)。

2.4 元模型驗(yàn)證

圖3 參考點(diǎn)周圍的徑向點(diǎn)分布情況

表1為數(shù)值模擬得到的25個(gè)設(shè)計(jì)點(diǎn)累計(jì)產(chǎn)量（用Eclipse軟件計(jì)算）以及采用3種元模型建模后求得的累計(jì)產(chǎn)量，表中的誤差為各模型的輸出響應(yīng)（即累計(jì)產(chǎn)量）相對數(shù)值模擬計(jì)算結(jié)果的誤差。為了檢驗(yàn)各元模型的有效性，選擇 7個(gè)試驗(yàn)點(diǎn)（在最小泄氣半徑范圍內(nèi)），分別采用數(shù)值模擬方法和3種元模型計(jì)算累計(jì)產(chǎn)量（見表2）。

表1 設(shè)計(jì)點(diǎn)數(shù)值模擬和元建模計(jì)算的累計(jì)產(chǎn)量

表2 試驗(yàn)點(diǎn)數(shù)值模擬和元建模計(jì)算的累計(jì)產(chǎn)量

由表1可知：①在25個(gè)設(shè)計(jì)點(diǎn)，基于二次模型、乘性模型、RBF模型預(yù)測的累計(jì)產(chǎn)量均與數(shù)值模擬預(yù)測結(jié)果有較好的一致性，回歸分析的復(fù)相關(guān)系數(shù)值分別為94.6%、96.6%和96.6%；②以25個(gè)設(shè)計(jì)點(diǎn)數(shù)值模擬預(yù)測累計(jì)產(chǎn)量為參照，RBF模型的預(yù)測準(zhǔn)確性優(yōu)于另外2種模型。

由表2可知：①在26號點(diǎn)，二次模型、乘性模型和 RBF模型的累計(jì)產(chǎn)量預(yù)測準(zhǔn)確性都較差，說明 26號點(diǎn)為低潛力點(diǎn)，這與實(shí)際情況一致；②在 7個(gè)試驗(yàn)點(diǎn)，二次模型、乘性模型、RBF模型的平均誤差分別為15.5%、14.1%和6.6%；③若剔除26號點(diǎn)，二次模型、乘性模型、RBF模型的平均誤差分別降為7.6%、8.2%和1.6%。

綜上所述，二次模型、乘性模型和RBF模型均有較好的累計(jì)產(chǎn)量預(yù)測準(zhǔn)確性，RBF模型的累計(jì)產(chǎn)量預(yù)測準(zhǔn)確性優(yōu)于另外兩種模型。主要原因可能是RBF模型的輸入變量更多，因此，如果能獲取各網(wǎng)格的氣藏參數(shù)數(shù)據(jù)，RBF模型是累計(jì)產(chǎn)量預(yù)測的最佳選擇。

2.5 元模型應(yīng)用

利用上述 3種元模型，可在幾秒鐘內(nèi)計(jì)算出整個(gè)氣藏空間中任何位置的累計(jì)產(chǎn)量，如果采用數(shù)值模擬方法，則需耗時(shí)數(shù)月。圖 4為基于二次模型求得的氣藏累計(jì)產(chǎn)量等值線圖，累計(jì)產(chǎn)量在12.8×106m3以上的區(qū)域?yàn)楦邼摿^(qū)，可作為建設(shè)生產(chǎn)井的備選區(qū)。

圖4 基于二次模型求得的氣藏累計(jì)產(chǎn)量等值線圖

3 結(jié)論

在建立元模型之前預(yù)先估算了氣藏的基本參數(shù)（孔隙度、滲透率及含水飽和度）并對氣藏進(jìn)行了數(shù)值模擬。

將氣藏內(nèi)某點(diǎn)18年規(guī)劃周期內(nèi)的累計(jì)產(chǎn)量表示為x方向的滲透率、y方向的滲透率、z方向的滲透率、孔隙度和含水飽和度的函數(shù)。采用HIP經(jīng)驗(yàn)函數(shù)排除了研究空間內(nèi)的低潛力點(diǎn)，并基于最大熵準(zhǔn)則選取了25個(gè)高潛力點(diǎn)作為設(shè)計(jì)點(diǎn)。

建立了3種元模型（二次模型、乘性模型和RBF模型），基于各模型預(yù)測了25個(gè)設(shè)計(jì)點(diǎn)和7個(gè)試驗(yàn)點(diǎn)處的累計(jì)產(chǎn)量，并結(jié)合數(shù)值模擬預(yù)測結(jié)果進(jìn)行了驗(yàn)證和對比。結(jié)果表明：與數(shù)值模擬方法相比，元建模方法的計(jì)算時(shí)間大幅縮短；各模型均有較好的累計(jì)產(chǎn)量預(yù)測準(zhǔn)確性；RBF模型的累計(jì)產(chǎn)量預(yù)測準(zhǔn)確性優(yōu)于另外2種模型。

符號注釋：

Kx——x方向的滲透率，10-3μm2；Ky——y 方向的滲透率，10-3μm2；Kz——z方向的滲透率，10-3μm2；φ——孔隙度，%；Sw——含水飽和度；Qacc——?dú)獠啬尘?8年內(nèi)的累計(jì)產(chǎn)量，106m3；Vb——網(wǎng)格體積，m3；Y——輸出響應(yīng)；X——輸入變量；F(X) ——近似響應(yīng)函數(shù)；ε——誤差；β0—β20——二次模型中的系數(shù)；V(·)——方差；γ0—γ10——乘性模型中的系數(shù)；λ0—λk——RBF模型中的系數(shù)；hj(X)——徑向基函數(shù)；n——變量個(gè)數(shù)；k——參考點(diǎn)周圍徑向點(diǎn)個(gè)數(shù)；c——RBF模型中的標(biāo)量；Xr——參考點(diǎn) Ar處的屬性向量；Xj——徑向點(diǎn) Aj處的屬性向量。

[1] Güyagüler B. Optimization of well placement and assessment of uncertainty[D]. Stanford： Stanford University, 2002.

[2] 龐偉, 陳德春, 張仲平, 等. 非均質(zhì)油藏水平井分段變密度射孔優(yōu)化模型[J]. 石油勘探與開發(fā), 2012, 39(2)： 214-221.Pang Wei, Chen Dechun, Zhang Zhongping, et al. Segmentally variable density perforation optimization model for horizontal wells in heterogeneous reservoirs[J]. Petroleum Exploration and Development,2012, 39(2)： 214-221.

[3] Yeten B, Durlofsky L J, Aziz K. Optimization of nonconventional well type, location and trajectory[R]. SPE 77565, 2002.

[4] Foroud T, Seifi A, Hassani H. Surrogate-based optimization of horizontal well placement in a mature oil reservoir[J]. Petroleum Science and Technology, 2012, 30(11)： 1091-1101.

[5] Mohammadi H, Seifi A, Foroud T. A robust Kriging model for predicting accumulative outflow from a mature reservoir considering a new horizontal well[J]. Journal of Petroleum Science and Engineering, 2012, 82/83： 113-119.

[6] 何東博, 王麗娟, 冀光, 等. 蘇里格致密砂巖氣田開發(fā)井距優(yōu)化[J]. 石油勘探與開發(fā), 2012, 39(4)： 458-464.He Dongbo, Wang Lijuan, Ji Guang, et al. Well spacing optimization for Sulige tight sand gas field, NW China[J]. Petroleum Exploration and Development, 2012, 39(4)： 458-464.

[7] Etminan A H, Seifi A. An improved model for geostatistical simulation of fracture parameters and their effect on static and dynamic models[J]. The Open Petroleum Engineering Journal, 2008, 1： 47-57.

[8] Mckay M D, Beckman R J, Conover W J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics, 2000, 42(1)： 55-61.

[9] Sacks J, Welch W, Mitchell T, et al. Design and analysis of computer experiments (with discussion)[J]. Statist Sci, 1989, 4(4)： 409-423.

[10] Ko C W, Lee J, Queyranne M. An exact algorithm for maximum entropy sampling[J]. Operations Research, 1995, 43(4)： 684-691.

[11] Seifi A, Karimifar M J. Entropy based spatial design： A genetic algorithm approach (case study)[J]. World Academy of Science,Engineering and Technology, 2008, 45： 637-646.

[12] Chen V C P, Tsui K L, Barton R R, et al. A review on design,modeling and applications of computer experiments[J]. IIE Transactions,2006, 38(4)： 273-291.

[13] Ball K, Sivakumar N, Ward J D. On the sensitivity of radial basis interpolation to minimal data separation distance[J]. Constructive Approximation, 1992, 8(4)： 401-426.

[14] Venkataraman R. Application of the method of experimental design to quantify uncertainty in production profiles[R]. SPE 59422, 2000.

[15] Jakobsson S, Patriksson M, Rudholm J, et al. A method for simulation based optimization using radial basis functions[J].Optimization and Engineering, 2010, 11(4)： 501-532.