Seifi A，Kazemzadeh M B，Mohammadi H

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

0 引言

一般采用數(shù)值模擬方法預測油氣藏的累計產量。數(shù)值模擬屬于密集型計算過程，計算量大，采用元建模方法并建立合適的元模型，可以減小計算量、提高計算效率。元模型是一種代理數(shù)學方程，主要用于近似表示復雜體系內的輸入-輸出關系[1]。元建模基本步驟[1-6]包括：①定義輸入變量及相應的輸出響應；②進行實驗設計（DoE），選擇樣本點，并確定各樣本點的響應；③構建元模型，建立輸入變量與輸出響應之間的關系；④通過試驗點驗證元模型的準確性。

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

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

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

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

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

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

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

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

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

2 元建模

2.1 輸入變量和輸出響應

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

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

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

2.2 實驗設計

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

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

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

圖2 由HIP經驗函數(shù)生成的儲量等值線圖

2.3 元模型建立

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

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

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

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

2.3.2 乘性模型

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

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

2.3.3 RBF模型

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

其中

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

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

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

2.4 元模型驗證

圖3 參考點周圍的徑向點分布情況

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

表1 設計點數(shù)值模擬和元建模計算的累計產量

表2 試驗點數(shù)值模擬和元建模計算的累計產量

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

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

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

2.5 元模型應用

利用上述 3種元模型，可在幾秒鐘內計算出整個氣藏空間中任何位置的累計產量，如果采用數(shù)值模擬方法，則需耗時數(shù)月。圖 4為基于二次模型求得的氣藏累計產量等值線圖，累計產量在12.8×106m3以上的區(qū)域為高潛力區(qū)，可作為建設生產井的備選區(qū)。

圖4 基于二次模型求得的氣藏累計產量等值線圖

3 結論

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

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

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

符號注釋：

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

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

[2] 龐偉, 陳德春, 張仲平, 等. 非均質油藏水平井分段變密度射孔優(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.