低滲透氣藏型儲氣庫儲層物性參數(shù)的反演分析

張金冬，王雪梅，譚羽非，于克成，張?zhí)鹛?/p>

(1.哈爾濱工業(yè)大學 建筑學院，哈爾濱 150006；2.寒地城鄉(xiāng)人居環(huán)境科學與技術工業(yè)和信息化部重點實驗室(哈爾濱工業(yè)大學)，哈爾濱 150090)

常規(guī)氣藏改建的地下儲氣庫進行數(shù)值模擬時，一般不考慮儲層物性參數(shù)隨地層壓力的變化情況[1]。低滲透氣藏改建的地下儲氣庫，儲層非均質性強、物性條件復雜，在地下儲氣庫的注氣過程中，隨著氣體的不斷注入，儲層的地層壓力從原始地層壓力開始不斷增加[2-3]。對于低滲透儲層來說，地層壓力的微小變化就會引起儲層孔隙度和滲透率的變化，進而影響地下儲層的滲流能力，最終影響地下儲氣庫的注入量[4]。文獻[5]認為滲透率越低的儲層，其滲透率隨地層壓力變化的越劇烈。文獻[6]認為氣藏開采時孔隙度的變化范圍較滲透率來說小得多，對氣藏的開采結果沒有影響。但是在儲氣庫的注采過程中，不能照搬氣藏開采的理論，在擴容建庫時，強注可能引起孔隙結構的變化，進而影響巖石的壓縮性，在注采過程中也需考慮其變化規(guī)律。因此，在分析低滲透氣藏改建的地下儲氣庫的注采運行過程時，需考慮滲透率和孔隙度隨地層壓力的變化。

為了獲得低滲透氣藏型儲氣庫注氣過程中儲層物性參數(shù)的變化規(guī)律，本文首先基于地質統(tǒng)計學中的變差函數(shù)理論確定儲層滲透率和孔隙度等參數(shù)的初始分布情況。然后依據(jù)反問題理論，利用已知井點處地層壓力的實測值和計算值之差構建目標函數(shù)，實現(xiàn)了對儲層滲透率和孔隙度的反演求解。本文通過案例證明了模型的正確性并利用最小二乘法擬合得到了滲透率和孔隙度與地層壓力之間的函數(shù)關系式。

1 模型的建立及求解

1.1 反問題理論基礎

反問題從數(shù)學模型的角度來看就是模型識別問題。反問題一般通過系統(tǒng)辨識或模型辨識來完成，借助數(shù)學物理方法，通過對微分方程中未知參數(shù)的確定，來完成對源的辨識[7-8]。在滲流力學領域，反問題一般是指從某些模型參數(shù)或者模型動態(tài)推斷或者識別整個模型。在求解正問題時，觀察數(shù)據(jù)的數(shù)目一般等于或大于待求參數(shù)的數(shù)目，這時正問題的解是唯一且穩(wěn)定的。而在求解反問題時，觀察數(shù)據(jù)的數(shù)目則小于待求參數(shù)的數(shù)目，這時模型的解不唯一，需要附加一定的條件保證解的存在性和唯一性[9]。這往往通過構造目標函數(shù)，使目標函數(shù)最小化來實現(xiàn)[10]。

考慮單相滲流試井壓力的數(shù)據(jù)時，待確定的模型參數(shù)是儲層孔隙度和滲透率。觀察的數(shù)據(jù)值為地層壓力，觀察的數(shù)據(jù)個數(shù)小于待測定的參數(shù)個數(shù)，為了獲得與觀察數(shù)據(jù)的個數(shù)相吻合的模型參數(shù)的實現(xiàn)，需利用地層壓力的實測值和計算值構造目標函數(shù)并在目標函數(shù)中引入正則化參數(shù)和光滑泛函來保證模型解的唯一性和確定性。

1.2 反演的目標函數(shù)

(1)

式中：Pi(X)為通過正問題求解得到的第i個測點的計算壓力值,X為待反演的物性參數(shù)，本文中指孔隙度和滲透率，本文待求解的參數(shù)是低滲透儲層的孔隙度和滲透率，已知參數(shù)是測點地層壓力，待求解的參數(shù)大于已知參數(shù)，這樣的反問題是不適定的。為了數(shù)值求解的穩(wěn)定性，本文利用正則化方法在泛函J(X)中引入光滑泛函，用如下泛函代替J(X)：

α1(‖M1X‖2+‖M2X‖2+‖M3X‖2)

(2)

式中：M1、M2、M3分別為x、y和z方向的二階光滑矩陣，α1為正則化參數(shù)。

數(shù)值模擬中利用隨機誤差可以得到共軛梯度法的收斂條件：

α1(‖M1X‖2+‖M2X‖2+

‖M3X‖2)

(3)

式中:N為測點個數(shù)，σ為測點壓力的殘差。

1.3 反演的方法

本文采用共軛梯度法對目標函數(shù)進行求解。共軛梯度法的循環(huán)方式為

Xn+1=Xn+αndn

(4)

式中:Xn為待反演的滲透率和孔隙度的第n次的預測值,an為迭代步長,dn為共軛梯度搜索方向，表示為

dn=-?Jα1(Xn)+βn-1dn-1

(5)

設xi為反演的滲透率與孔隙度的向量，則目標函數(shù)的梯度向量?Jα1可以表示為

(6)

其中

因此，求解目標函數(shù)的關鍵是求得儲層壓力對滲透率和孔隙度的變化率。

1.4 地層壓力對儲層滲透率和孔隙度變化率的求解

氣體在低滲透氣藏改建的地下儲氣庫中的流動屬于低速流，遵循達西滲流規(guī)律[11-12]，其控制方程可以寫成如下形式：

(7)

式中:μ為流體黏度,K為儲層絕對滲透率,P為儲層的壓力,ρ為氣體密度,φ為儲層孔隙度,c為巖石的壓縮系數(shù),q為源(匯)項，表示單位時間單位地層體積注入或采出的流量，注入井取正值，采出井取負值。

邊界條件：

P=P0,τ=0

(8)

(9)

(10)

(11)

經(jīng)分部積分可得：

(12)

[K]{P}+[D]{?P/?τ}=[E]{Q}+[R]{f}

(13)

式中：{P}為不包括給定值在內(nèi)的節(jié)點壓力向量,{f}、{Q}分別為f、Q的節(jié)點向量,[K],[D],[E],[R]是與變量相關的矩陣。

將方程中的已知量和未知量分離，有限元列式可以表示成如下形式：

(14)

將式(14)可進一步改寫為

[Ka]{P}+[Da]{?P/?τ}={W}

(15)

對式(15)進行時間上的Galerkin的差分離散有

(16)

(17)

(18)

由式(17)、(18)可知，要確定在τ時刻壓力對滲透率和孔隙度的變化率，必須已知τ-1時刻的儲層滲透率和孔隙度的分布，以此類推，要想求解各時刻地層壓力對孔隙度和滲透率的變化率，必須知道儲層滲透率和孔隙度的初始分布。

1.5 儲層參數(shù)初始分布的確定

考慮單相滲流試井壓力的數(shù)據(jù)時，待確定的模型參數(shù)是網(wǎng)格的孔隙度和滲透率，流體參數(shù)已知。觀察的數(shù)據(jù)值為地層壓力，觀察的數(shù)據(jù)個數(shù)小于待測定的參數(shù)個數(shù)，為了獲得與觀察數(shù)據(jù)的個數(shù)相吻合的模型參數(shù)的實現(xiàn)，基于地質統(tǒng)計學的變差函數(shù)理論確定儲層滲透率和孔隙度的初始分布。

在地質統(tǒng)計學的差值方法中，直接由已知數(shù)據(jù)計算計算出來的實驗變差函數(shù)可能導致奇異矩陣多解，因此一般用被證明了可以保證方程組有唯一解和非負均方差的理論模型來代替實驗變差函數(shù)，其中比較常用的理論模型是球狀模型，其表達式如下:

(19)

式中:C0為塊金值，表示在很短的距離內(nèi)變量的空間變異性,C為拱高，表示區(qū)域化變量在空間上的變異性的程度,C與C0之和表示基臺值，反映的是區(qū)域化變量在空間上的總變異性,a為變程，反映的是區(qū)域化變量的變異范圍。

對式(19)進行求解，需已知C0、C和a的值。根據(jù)文獻[13]，C0、C和a的值可由已知數(shù)據(jù)的均值和方差進行計算得到。對于球狀模型的求解，可以利用Kriging插值法實現(xiàn)。Kriging法通過引進以距離為自變量的變差函數(shù)來計算權值，由于變差函數(shù)既可以反映自變量的空間結構特性，又可以反映變量的隨機分布特性，利用Kriging方法進行空間數(shù)據(jù)差值往往可以獲得理想的效果[14]。

本文確定儲層初始參數(shù)分布時，觀測井點的物性參數(shù)已知，通過擬合變差函數(shù)曲線,采用Kriging方法確定低滲透氣藏改建儲氣庫儲層參數(shù)的分布。

1.6 目標函數(shù)的求解

首先利用Kriging插值法確定儲層滲透率和孔隙度的初始分布，然后利用參數(shù)的分布確定各時刻地層壓力對滲透率和孔隙度的變化率，將求得的變化率式(17)、(18)代入式(6)，然后使用Newton-Raphson法進行求解。求解的具體步驟如下：

1)根據(jù)低滲透氣藏型儲氣庫各已知井點處的滲透率和孔隙度，確定儲層參數(shù)的初始分布情況，計算待測井點處地層的壓力值。

2)測量已知井點處的壓力值，將計算值和測量值代入方程(1)，判斷是否滿足收斂條件，若滿足執(zhí)行步驟4)，不滿足則執(zhí)行步驟3)。

3)將求得新的壓力對滲透率和孔隙度的變化率代入方程(1)～(5)，然后重復步驟1)和2)。

4)輸出得到的滲透率和孔隙度的值，結束程序。

2 算例分析

借鑒文獻[15]儲氣庫的注采運行數(shù)據(jù)，通過注氣流量和觀察井的壓力，利用滲流反問題構建目標函數(shù)計算注氣后儲層物性參數(shù)的變化。考慮如圖1所示的儲層，假定上、下、左、右均為不滲透邊界，在其中選取面積為400 m×400 m的區(qū)域，該儲層分布有9口井，其中1、3、5、7和9號井以定注入量q=21×104Nm3/d注入，其他井為觀測井。計算過程中地層物性參數(shù)和其他參數(shù)值見表1。

表1 計算模型的物性參數(shù)與其他參數(shù)Tab.1 Physical parameters and other parameters of calculation model

圖1 儲氣庫儲層井位分布Fig.1 Well location distribution of underground gas storage

測點壓力的殘差為σd=0.006 MPa，在注入井和觀測井處的初始孔隙度和滲透率已知，見表2。其中，滲透率的單位為μm2。

表2 觀測井點處的孔隙度和滲透率Tab.2 Porosity and permeability at observation wells

2.1 模型驗證

為驗證模型的正確性，本文首先根據(jù)地質統(tǒng)計學中的變差函數(shù)理論，以已知井點處的滲透率和孔隙度為基礎，確定了儲層參數(shù)的初始分布。通過計算地層壓力對滲透率和孔隙度的變化率，利用共軛梯度法確定了注氣后儲層參數(shù)的分布。利用新的儲層參數(shù)重新計算了低滲透氣藏型儲氣5#注入井的壓力，將計算結果和5#注入井的已測壓力進行擬合，并將其與滲透率和孔隙度當作常數(shù)的傳統(tǒng)計算方法進行對比，結果如圖2所示。由圖2可知，重新計算的5#注入井壓力和已測壓力的曲線基本重合，二者誤差最大為2.56%，平均值為0.99%，說明計算結果準確。而采用傳統(tǒng)計算方法得到的5#注入井壓力和已測壓力二者誤差最高可達7.89%。比較采用兩種方法重新計算獲得的5#注入井壓力可知，在注氣初期，二者相差不大，曲線基本重合。而在注氣中后期，二者之間的差值逐漸增大，在注氣結束時達到最大，二者相差1.1 MPa。

圖2 重新計算得到的5#注入井壓力和已測壓力對比曲線Fig.2 Comparison of recalculated and measured pressure of 5# injection well

2.2 儲層滲透率和孔隙度變化的計算

根據(jù)已知井點處滲透率和孔隙度的初始值，可以計算出滲透率和孔隙度的均值分別為2.05×10-3μm2和0.18，方差分別為0.165和0.005，根據(jù)均值和方差可取變程a為162 m，C0和C的值分別為0和0.002 5，將a、C0和C的值代入式(19)，利用Kriging插值法可以得到儲氣庫初始滲透率和孔隙度分布等值線和云圖，如圖3、4和圖5、6所示。

圖3 儲層初始滲透率分布等值線圖Fig.3 Contour map of initial reservoir permeability distribution

圖4 儲層初始孔隙度分布等值線圖Fig.4 Contour map of initial reservoir porosity distribution

圖5 儲層初始滲透率分布云圖Fig.5 Cloud map of initial reservoir permeability distribution

圖6 儲層初始孔隙度分布云圖Fig.6 Cloud map of initial reservoir porosity distribution

對比圖3、圖5和圖4、圖6可以發(fā)現(xiàn)，低滲透氣藏改建地下儲氣庫在注氣初期，滲透率和孔隙度的分布趨勢基本一致，二者呈現(xiàn)一定的相關性。

利用獲得的儲層滲透率和孔隙度初始分布作為已知條件，將其代入到式(14)、(15)中，計算得到各時刻地層壓力對滲透率和孔隙度的變化率，將計算得到的變化率代入到最優(yōu)化目標函數(shù)中，并利用1#、3#、5#、7#和9#注采井的測量壓力作為擬合條件，反演得到的低滲透氣藏改建儲氣庫儲層滲透率和孔隙度的分布情況如圖7、8和圖9、10所示。

圖7 反演儲層滲透率分布等值線圖Fig.7 Contour map of calculated permeability distribution

圖8 反演儲層孔隙度分布等值線圖Fig.8 Contour map of calculated porosity distribution

圖9 反演儲層滲透率分布云圖Fig.9 Cloud map of calculated permeability distribution

圖10 反演儲層孔隙度分布云圖Fig.10 Cloudmap of calculated porosity distribution

對比圖3、圖7和圖4、圖8可以發(fā)現(xiàn)，隨著氣體的不斷注入，在儲層各個位置的滲透率和孔隙度也隨之變化，并且各位置的滲透率和孔隙度變化幅度都不相同，二者之間的相關性不再一致。對比圖3、圖5和圖7、圖9可以發(fā)現(xiàn)，低滲透儲層滲透率經(jīng)過注氣后其分布規(guī)律發(fā)生了較大的變化，其中變化最大的位置是在5#注采井附近，其滲透率從1.66×10-3μm增加到2.81×10-3μm。由此可見，對于低滲透儲層改建的地下儲氣庫來說，要分析其注采氣過程，必須考慮滲透率隨著孔隙流體壓力的變化情況。

2.3 滲透率和孔隙度與地層壓力的擬合

為了分析低滲透儲層滲透率和孔隙度隨著孔隙流體壓力的變化規(guī)律，本文對滲透率變化幅度最大的注入井附近的各地層壓力下的滲透率和孔隙度進行了擬合，得到的滲透率和孔隙度隨孔隙流體壓力的變化結果如圖11、12所示。

圖11 滲透率隨地層壓力變化規(guī)律Fig.11 Variation of permeability with formation pressure

圖12 孔隙度隨地層壓力變化規(guī)律Fig.12 Variation of porosity with formation pressure

圖11、12反映了低滲透儲層型地下儲氣庫在注氣時，滲透率和孔隙度隨孔隙流體壓力的變化規(guī)律。從圖11、12以及擬合的表達式可以看出，低滲透氣藏型儲氣庫隨著天然氣的不斷注入，儲層地層壓力不斷增加，滲透率和孔隙度也隨之增加。但是相比較滲透率隨地層壓力的變化而言，孔隙度隨著儲層壓力的增加變化極其緩慢，在建立方程時可以忽略。但是在低滲透氣藏型儲氣庫中，滲透率隨著壓力變化極大，不可忽略，在建立儲氣庫的注采方程時，不可當作常數(shù)處理，在建立滲流方程時需將其擬合成儲層地層壓力的表達式。

3 結 論

1)本文首先利用地質統(tǒng)計學中的變差函數(shù)理論，依據(jù)已知觀測井處的孔隙度和滲透率，確定了儲層初始孔隙度和滲透率的分布，并基于反問題理論，利用儲氣庫井點處地層壓力的實測值和計算值的差值構建目標函數(shù)，反演得到了注氣后儲層孔隙度和滲透率的分布；

2)利用觀測井的壓力實測數(shù)據(jù)，證明了構建的反問題模型的正確性，并將反演的滲透率和孔隙度與壓力之間的關系進行擬合，得到了滲透率和孔隙度與地層壓力之間的關聯(lián)式。