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

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

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

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

1 模型的建立及求解

1.1 反問題理論基礎(chǔ)

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

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

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

(1)

式中：Pi(X)為通過正問題求解得到的第i個(gè)測(cè)點(diǎn)的計(jì)算壓力值,X為待反演的物性參數(shù)，本文中指孔隙度和滲透率，本文待求解的參數(shù)是低滲透儲(chǔ)層的孔隙度和滲透率，已知參數(shù)是測(cè)點(diǎn)地層壓力，待求解的參數(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ù)值模擬中利用隨機(jī)誤差可以得到共軛梯度法的收斂條件：

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

‖M3X‖2)

(3)

式中:N為測(cè)點(diǎn)個(gè)數(shù)，σ為測(cè)點(diǎn)壓力的殘差。

1.3 反演的方法

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

Xn+1=Xn+αndn

(4)

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

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

(5)

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

(6)

其中

因此，求解目標(biāo)函數(shù)的關(guān)鍵是求得儲(chǔ)層壓力對(duì)滲透率和孔隙度的變化率。

1.4 地層壓力對(duì)儲(chǔ)層滲透率和孔隙度變化率的求解

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

(7)

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

邊界條件：

P=P0,τ=0

(8)

(9)

(10)

(11)

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

(12)

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

(13)

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

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

(14)

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

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

(15)

對(duì)式(15)進(jìn)行時(shí)間上的Galerkin的差分離散有

(16)

(17)

(18)

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

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

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

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

(19)

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

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

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

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

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

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

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

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

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

2 算例分析

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

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

圖1 儲(chǔ)氣庫(kù)儲(chǔ)層井位分布Fig.1 Well location distribution of underground gas storage

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

表2 觀測(cè)井點(diǎn)處的孔隙度和滲透率Tab.2 Porosity and permeability at observation wells

2.1 模型驗(yàn)證

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

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

2.2 儲(chǔ)層滲透率和孔隙度變化的計(jì)算

根據(jù)已知井點(diǎn)處滲透率和孔隙度的初始值，可以計(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插值法可以得到儲(chǔ)氣庫(kù)初始滲透率和孔隙度分布等值線和云圖，如圖3、4和圖5、6所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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