低滲透油藏CO2驅(qū)試井解釋方法

蘇玉亮，蔡明玉，孟凡坤，范理堯，李 蕾

（1.中國石油大學(xué)（華東）石油工程學(xué)院，山東青島 266580；2.中國石油大學(xué)勝利學(xué)院，山東東營 257061）

受制于低滲透油藏水驅(qū)注入能力，CO2驅(qū)已逐漸成為提高低滲透油藏采收率的主要方法之一，具有廣泛的應(yīng)用前景［1-4］。在CO2注入過程中開展試井測試，可了解油藏性質(zhì)及井的狀況，對指導(dǎo)后續(xù)生產(chǎn)具有重要意義［5-6］。目前，針對CO2驅(qū)，提出的多為復(fù)合油藏試井模型，TANG等率先揭示了CO2驅(qū)流體分布特征，將油藏劃分為CO2區(qū)、CO2-原油過渡區(qū)及未波及原油區(qū)［7］；隨后，AMBASTHA 等建立了三區(qū)復(fù)合試井解釋模型，分析試井曲線變化規(guī)律［8］，ISSAKA 等假定過渡區(qū)流體性質(zhì)為指數(shù)變化，給出了變性質(zhì)過渡區(qū)復(fù)合模型半解析解［9］。在中國，王敬瑤運(yùn)用傳統(tǒng)三區(qū)復(fù)合模型，分析某CO2驅(qū)區(qū)塊試井?dāng)?shù)據(jù)，但未考慮CO2驅(qū)流體分布特征［10］；朱建偉等基于兩區(qū)復(fù)合模型，考慮流體性質(zhì)變化所產(chǎn)生壓降，但試井曲線無明顯變化規(guī)律［11］；蘇玉亮等充分考慮了CO2驅(qū)替過程中流體冪律型、非連續(xù)型等變化方式，建立了三區(qū)復(fù)合模型，分析其對試井曲線的影響［12-14］。近年來，眾多學(xué)者又考慮應(yīng)力敏感、水平井、垂直裂縫直井等復(fù)雜條件，建立多區(qū)復(fù)合模型，研究試井曲線變化特征［15-17］。

針對目前研究集中于試井曲線分析，而缺乏對CO2驅(qū)試井解釋方法研究的現(xiàn)狀，考慮啟動壓力梯度、應(yīng)力敏感等低滲透油藏特性以及流體性質(zhì)變化，建立改進(jìn)復(fù)合區(qū)試井模型，并綜合SPSA 自動搜索優(yōu)化算法，最終形成CO2驅(qū)試井解釋方法流程。

1 CO2驅(qū)試井模型建立與求解

1.1 物理模型

由封閉邊界低滲透油藏CO2驅(qū)三區(qū)復(fù)合模型剖面（圖1）可見，一區(qū)為驅(qū)替波及區(qū)（CO2區(qū)）；二區(qū)為CO2-原油過渡區(qū)（CO2+原油輕組分），該區(qū)流體性質(zhì)與流體組分及油藏壓力、溫度等緊密相關(guān)，難以用特定函數(shù)統(tǒng)一表征；三區(qū)代表驅(qū)替未波及區(qū)（原油區(qū)）。溫度（304.14K）和臨界壓力（7.382MPa），故注入CO2實際處于超臨界狀態(tài)，其性質(zhì)類似于液體，因此可簡化模型為單相滲流。同時考慮低滲透油藏非達(dá)西滲流的特點(diǎn)，在模型中加入表征低滲透特點(diǎn)的

一般地層條件下的溫度和壓力均超過CO2臨界啟動壓力梯度及地層壓力敏感性。模型的基本假設(shè)條件如下：①地層為水平、等厚且恒溫的封閉地層，油藏原始地層壓力在各處相等。②注入井CO2注入量恒定。③忽略重力和毛管力的影響，考慮表皮污染和井筒存儲效應(yīng)。④一區(qū)、三區(qū)中流體均微可壓縮，且流體性質(zhì)均一。⑤二區(qū)中流體黏度及壓縮系數(shù)隨半徑呈連續(xù)（冪律、指數(shù)）、非連續(xù)等多種變化形式。⑥整個流動過程為單相滲流，受流體性質(zhì)差別影響，各區(qū)啟動壓力梯度不同，均服從非達(dá)西滲流規(guī)律，滲流速度表達(dá)式為：

⑦考慮CO2注入過程中地層壓敏效應(yīng)，其滲透率變化可表示為：

1.2 數(shù)學(xué)模型

根據(jù)物理模型的假設(shè)，推導(dǎo)滲流數(shù)學(xué)模型并進(jìn)行無因次化后可得：

模型初始條件為：

（3）式中，在R1D＜rD＜R2D的條件下，V，W，Q的定義式分別為：

當(dāng)黏度及壓縮系數(shù)隨半徑成解析函數(shù)形式變化時，對于給定的rD可以直接計算出該半徑處的V，W，Q值；對于黏度及壓縮系數(shù)隨半徑成不規(guī)則變化時，對于給定的rD可以采用插值法（如樣條插值）計算出該半徑處的V，W，Q值。

（3）—（8）式中其他各無因次參數(shù)定義式為：

1.3 模型求解

借鑒文獻(xiàn)［14］中所采用的求解方法，首先對滲流方程及邊界條件（（1）—（4）式）進(jìn)行對數(shù)坐標(biāo)的轉(zhuǎn)化，并對主控方程（（3）式）進(jìn)行離散可得：

設(shè)一區(qū)、二區(qū)及三區(qū)網(wǎng)格數(shù)分別為N1，N2，N3，則對于一區(qū)1＜i＜N1：

對于三區(qū)N2＜i＜N3：

內(nèi)邊界（i=1）離散為：

其中：

一區(qū)-二區(qū)（i=N1）、二區(qū)-三區(qū)（i=N2）交界面及外邊界（i=N3）離散方程分別為：

對（27）—（30），（32）—（34）式進(jìn)行統(tǒng)一整理，可得到線性代數(shù)方程組：

利用Newton-Raphson方法構(gòu)造牛頓迭代矩陣：

即可求解關(guān)于壓力的非線性方程組：

2 SPSA搜索優(yōu)化算法

SPSA（Simultaneous Perturbation Stochastic Approximation）方法最先由SPALL［18］于1992年提出，與傳統(tǒng)的有限差分近似法（梯度法）不同，該方法通過設(shè)置一定的隨機(jī)變量，對所有控制變量進(jìn)行同步擾動，以此來獲得擾動梯度，進(jìn)行損失函數(shù)等性能指標(biāo)函數(shù)的計算。盡管擾動梯度是隨機(jī)的，但通過選擇適當(dāng)?shù)碾S機(jī)變量，可以保證對于最小化問題來講恒為下山方向，且期望值為真實梯度。算法的計算流程包括以下6 個步驟［19］：①模型參數(shù)初始化及選擇。設(shè)定計數(shù)標(biāo)記k=0。選擇初始估計值θ0，根據(jù)迭代系數(shù)αk和ck計算表達(dá)式（αk=a/（A+k+1），ck=c/（k+1）γ），確定非負(fù)系數(shù)a，c，A，α和γ。②隨機(jī)擾動向量生成。運(yùn)用Bernoulli±1 分布生成p維擾動向量Δk，為提高早期迭代穩(wěn)定性，采用多次計算后的平均梯度。③目標(biāo)函數(shù)估計?；谀壳爸蹈浇耐綌_動，根據(jù)步驟①和②中設(shè)定的ck與Δk，計算兩者目標(biāo)函數(shù)值。④梯度的近似計算。生成同步擾動變化梯度，其表達(dá)式為：

gk()中p個變量反映中所有元素同步擾動；⑤參數(shù)迭代更新。對進(jìn)行更新的表達(dá)式為：

⑥迭代終止條件判斷。返回到第②步，運(yùn)用k+1 代替k。當(dāng)相鄰2 次迭代步驟間計算結(jié)果差值較小，滿足設(shè)定閾值，或達(dá)到最大迭代步數(shù)時終止運(yùn)算。

3 試井解釋與自動擬合方法

對于低滲透油藏CO2驅(qū)的試井解釋，商業(yè)試井軟件不易考慮其低滲透及CO2驅(qū)的特性。因此，可構(gòu)建低滲透油藏CO2驅(qū)試井解釋自動擬合方法，其基本思路為：先通過常規(guī)試井解釋方法得到初始參數(shù)值，再通過自動搜索算法結(jié)合改進(jìn)的模型進(jìn)行擬合矯正。其具體流程如下：①擬定初始解釋參數(shù)，可采用常規(guī)試井解釋軟件，得到基于均質(zhì)復(fù)合區(qū)模型的解釋值作為初始值。②根據(jù)礦場、實驗或者數(shù)值模擬數(shù)據(jù)，給定啟動壓力梯度、壓敏系數(shù)以及過渡區(qū)參數(shù)的變化形式。③根據(jù)上述參數(shù)重新計算壓力及壓力導(dǎo)數(shù)曲線。④計算目標(biāo)函數(shù)。⑤利用SPSA 算法在變量定義域(xmin,xmax)中進(jìn)行搜索，使得目標(biāo)函數(shù)值最小。⑥調(diào)整搜索步長和搜索次數(shù)，使所求的x滿足精度要求，并符合實際情況。⑦繪制優(yōu)化后的壓力及壓力曲線圖。

目標(biāo)函數(shù)為：

目標(biāo)函數(shù)采用的是對壓力曲線和壓力導(dǎo)數(shù)曲線同時進(jìn)行擬合。由于目前已有的低滲透油藏CO2實際試井?dāng)?shù)據(jù)極為有限，因此這里采用文獻(xiàn)［3］中CO2試井?dāng)?shù)據(jù)進(jìn)行算例示范操作。在該實例中，因采用壓力恢復(fù)試井，測試壓力隨時間的增長而增大，初始增大速率較小，后期逐步穩(wěn)定，由于3 個區(qū)間存在界面，使得壓力增大速率出現(xiàn)局部突變。通過商業(yè)軟件采用二區(qū)復(fù)合油藏模型（表1）初步擬合參數(shù)，初步試井參數(shù)解釋結(jié)果為：半徑為31.3 m，表皮系數(shù)為2.72，井儲系數(shù)為0.004 4，流度比為1.09，導(dǎo)壓系數(shù)為0.37，滲透率為31.3 mD，指數(shù)I為0.7，所得初始擬合曲線見圖2。

將基礎(chǔ)數(shù)據(jù)及初始值導(dǎo)入自動擬合程序。假設(shè)該油藏為低滲透油藏，且注二氧化碳時啟動壓力梯度為0.01 MPa/m，地層壓敏系數(shù)為0.1 mD/MPa。同時，假定在上述二區(qū)模型中存在過渡帶。由初步試井參數(shù)解釋結(jié)果可以看出，兩區(qū)流度比值相差較小，因此黏度的變化在此算例中忽略不計，仍設(shè)為恒定值。在過渡區(qū)域內(nèi)壓縮系數(shù)的變化可簡化為呈指數(shù)形式變化（由組分模型計算歸納得到），即：此時所得優(yōu)化前曲線見圖2。

表1 油藏參數(shù)Table1 Reservoir parameters

圖2 優(yōu)化曲線前后對比Fig.2 Comparison of curves before and after optimization

由于所需搜索的各個參數(shù)在數(shù)量級以及數(shù)值上具有較大差異，因此，為了減小搜索范圍，提高搜索速度，采用對各參數(shù)的變化倍數(shù)進(jìn)行搜索的方法。即初始搜索值為全為1的向量。根據(jù)油藏實際，可給定參數(shù)變化倍數(shù)的上下限。最終優(yōu)化所得曲線見圖2，其最小二乘值變化如圖3所示。

圖3 優(yōu)化過程中最小二乘值變化Fig.3 Change of target least-square value during optimization

根據(jù)最終優(yōu)化的向量，可得出優(yōu)化后的數(shù)值（表2）。由表2 可見，在該算例中，滲透率、井儲系數(shù)、導(dǎo)壓系數(shù)以及指數(shù)I的變化最為明顯。除了井儲系數(shù)，其余參數(shù)都與油藏低滲透特性及非均質(zhì)流體有直接關(guān)系。當(dāng)存在啟動壓力梯度及壓敏效應(yīng)時，理論曲線會向上抬高；而滲透率增大時，會使理論曲線相對下降。因此優(yōu)化后，滲透率解釋值增大。而當(dāng)考慮非均質(zhì)壓縮系數(shù)流體時，指數(shù)I與導(dǎo)壓系數(shù)與其直接相關(guān)，因此其參數(shù)解釋值變化較大。

表2 優(yōu)化后參數(shù)變化值Table2 Parameters after optimization

4 結(jié)論

考慮低滲透油藏CO2驅(qū)替過程中存在啟動壓力梯度以及壓敏效應(yīng)的特點(diǎn)，同時對CO2驅(qū)所導(dǎo)致的流體性質(zhì)變化運(yùn)用解析或插值函數(shù)進(jìn)行描述，建立了低滲透油藏CO2驅(qū)非均質(zhì)流體三區(qū)復(fù)合模型，并進(jìn)行了數(shù)值求解。

基于改進(jìn)的復(fù)合模型，結(jié)合SPSA 搜索優(yōu)化算法，形成了一套低滲透油藏CO2驅(qū)試井解釋方法，并對實際數(shù)據(jù)進(jìn)行了算例計算。

參數(shù)修正解釋結(jié)果表明，當(dāng)考慮低滲透油藏低滲透特性時，對儲層滲透率的參數(shù)解釋影響較大；當(dāng)考慮流體非均質(zhì)分布時，與之相關(guān)的指數(shù)系數(shù)和導(dǎo)壓系數(shù)等的參數(shù)解釋值會受到較大影響。

符號解釋

rw——井筒原始半徑，cm；R——區(qū)域半徑，m；下標(biāo)1，2，3——一區(qū)、二區(qū)和三區(qū)；v——流體滲流速度，cm/s；K——地層滲透率，mD；μ——流體黏度，mPa·s；?p——壓力梯度，10-1MPa/cm；G——啟動壓力梯度，10-1MPa/cm；ε——壓敏系數(shù)，10-1MPa/cm；p——壓力，10-1MPa；下標(biāo)D——無因次；r——儲層徑向半徑，cm；S——井筒表皮系數(shù)；t——時間，s；V，W，Q——表征區(qū)域二中流體性質(zhì)變化的中間變量，無因次；η12，η13——一區(qū)和二區(qū)、一區(qū)和三區(qū)導(dǎo)壓系數(shù)比；C——井儲系數(shù)，cm3/10-1MPa；pw——井底壓力，10-1MPa；M12，M23——一區(qū)和二區(qū)、二區(qū)和三區(qū)流體流度比；Ct——地層壓縮系數(shù)，10 MPa-1；h——地層厚度，cm；q——注入量，cm3/s；Bg——注入流體的體積系數(shù)，cm3/cm3；pi——原始地層壓力，10-1MPa；?——地層孔隙度；i——網(wǎng)格數(shù)，1＜i＜N3；j——時間步；x——對數(shù)坐標(biāo)轉(zhuǎn)換結(jié)果，x=lnr；F——?dú)埐罹仃?；J——系數(shù)矩陣；P——壓力矩陣；gk()——同步擾動變化梯度——第k次迭代的估計值；Δki——向量Δk中的第i個元素（可能是隨機(jī)變量±1）；——一區(qū)外邊界處的壓縮系數(shù)，10 MPa-1；I——壓縮系數(shù)變化指數(shù)，初始值設(shè)為0.7。