海域天然氣水合物水平井儲層關(guān)鍵參數(shù)反演方法

張 鑫, 黃芳飛, 譚玉陽, 梁前勇, 董一飛, 何 川

(1. 北京大學(xué)地球與空間科學(xué)學(xué)院,北京 100871; 2.中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局,廣東廣州 510760;3.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島 266100; )

天然氣水合物是甲烷等小型氣體分子在低溫高壓條件下與水分子結(jié)合形成的似冰狀固體,絕大部分賦存于大陸邊緣陸坡的海底沉積層中,被認(rèn)為是21世紀(jì)理想的替代能源[1-3]。2002年以來,加拿大[4]、美國[5]、日本[6]和中國進(jìn)行了一系列天然氣水合物開采試驗(yàn)。相比陸上,海域天然氣水合物直井開采試驗(yàn)取得了重大突破[7-8],但距離海域水合物經(jīng)濟(jì)可行的開采目標(biāo)仍有很大差距[9-10],而利用水平井進(jìn)行海域天然氣水合物開發(fā)是提高單井產(chǎn)氣量的有效手段。水平井開采時(shí),需要對不同井段的產(chǎn)氣量進(jìn)行動態(tài)評價(jià)。常規(guī)油氣藏開發(fā)所使用的產(chǎn)能評價(jià)方法(如試井分析法、示蹤劑分析法等)均假定目標(biāo)井段在分段處理時(shí)被完全隔離[11-13],無法將其直接應(yīng)用于海域天然氣水合物的產(chǎn)量評價(jià)。測斜儀是一種測量地層傾斜場(位移場的梯度)變化的儀器[14-17],由于海域天然氣水合物儲層通常埋藏較淺[18],開采過程中的儲層體積變化會導(dǎo)致較大幅度的地層變形,繼而引發(fā)海底傾斜場的變化,這就使得測斜儀監(jiān)測技術(shù)成為解決海域水合物開采產(chǎn)氣量評價(jià)問題的可行手段。應(yīng)用測斜儀進(jìn)行儲層參數(shù)反演時(shí),水平井不同井段體積參數(shù)間的互補(bǔ)性可能導(dǎo)致目標(biāo)函數(shù)陷入較大的收斂區(qū)域,此時(shí)反演方法所能達(dá)到的精度將直接影響著反演結(jié)果的好壞。常用的反演算法(如直接網(wǎng)格搜索算法、梯度法和模擬退火法等)[19-21]對此無法取得良好效果。筆者提出一種基于鄰域算法[22-23]的儲層參數(shù)反演方法,以南海神狐海域?yàn)槔?研究不同噪聲程度和模型間距對反演結(jié)果的影響,并與直接網(wǎng)格搜索算法、模擬退火法的反演結(jié)果進(jìn)行對比,以驗(yàn)證提出方法的有效性。

1 方法原理

1.1 海底傾斜場正演模擬方法

圖1為水平井開采單段模型與測斜儀觀測系統(tǒng)示意圖。以該水平井射孔段中心點(diǎn)在海底面的投影點(diǎn)為笛卡爾坐標(biāo)系的原點(diǎn),定義X軸、Y軸、Z軸的正方向分別為正東、正北和豎直向上方向。海域天然氣水合物水平井開采時(shí),賦存于儲層孔隙或充填于巖石骨架中的固態(tài)水合物率先在射孔段附近發(fā)生分解,并逐漸在一定范圍內(nèi)形成以射孔段為中心的水合物分解區(qū)域?？紤]到邊際效應(yīng),可不失一般性地假定該水合物分解區(qū)域以射孔段為中心呈圓柱狀水平展布。在該區(qū)域內(nèi),由固態(tài)水合物分解的部分天然氣和水將與儲層中的自由水和可能存在的游離氣一起沿著氣水通道產(chǎn)出。開采過程中,實(shí)際產(chǎn)出的天然氣包括固態(tài)水合物分解的天然氣產(chǎn)出的部分和儲層中可能賦存的游離氣,產(chǎn)出的水包括固態(tài)水合物分解水產(chǎn)出的部分和儲層中的自由水。除此之外,開采過程中還可能伴隨一定的出砂(骨架顆粒),而氣、水、砂的產(chǎn)出將導(dǎo)致該分解區(qū)域的儲層體積變化,進(jìn)而引起海底面的傾斜場變化。

圖1 水平井開采單段模型與測斜儀觀測系統(tǒng)示意圖Fig.1 Single section model of horizontal well exploitationand schematic diagram of tiltmeter observation system

基于Okada[24-25]和張鑫等[26]的研究,本文中用塊狀模型等效水合物分解區(qū)域內(nèi)因氣、水、砂產(chǎn)出而導(dǎo)致的儲層體積變化,模型體積即為水合物分解區(qū)域內(nèi)產(chǎn)出的氣、水、砂所占據(jù)的全部孔隙空間之和。模型長度L為該水平井段的射孔段長度,且該模型的寬W和高H可認(rèn)為相等,模型中心為射孔段中心點(diǎn)位置c(x0,y0,z0)。此時(shí),開采過程中該水合物分解區(qū)域在海底面產(chǎn)生的傾斜場變化可由塊狀模型產(chǎn)生的傾斜場變化等效表示。因此水平井多段開采導(dǎo)致的海底面傾斜場變化可通過將各段模型引起的傾斜場變化進(jìn)行矢量疊加得到。

海域天然氣水合物水平井開采時(shí),可實(shí)時(shí)獲取總產(chǎn)氣量、總產(chǎn)水量和總出砂量,而難以直接獲取不同井段的產(chǎn)氣、產(chǎn)水和出砂量,這就為多段開采的產(chǎn)氣量評價(jià)帶來困難?？紤]到開采過程中水合物分解區(qū)域內(nèi)氣、水、砂的產(chǎn)出將導(dǎo)致儲層的體積變化。當(dāng)不同井段產(chǎn)水和出砂量較小或占該段儲層變化體積的比例一致時(shí),不同井段的產(chǎn)氣量與該段水合物分解區(qū)域內(nèi)的儲層體積變化呈正相關(guān)性。此時(shí),可通過海底測斜儀監(jiān)測到的因儲層體積變化導(dǎo)致的傾斜場變化來反演不同井段的儲層體積變化量,繼而進(jìn)行海域天然氣水合物水平井開采的多段產(chǎn)氣量評價(jià)。

1.2 鄰域算法

1.2.1 算法簡介

鄰域算法[22-23]用于在多維參數(shù)空間中尋找到可較好擬合觀測數(shù)據(jù)的模型集合,其基本元素為Voronoi單元,通過將Voronoi單元中每個(gè)點(diǎn)對應(yīng)的目標(biāo)函數(shù)設(shè)為常數(shù),可得目標(biāo)函數(shù)值的一個(gè)近似,稱為NA面。NA面提供了對參數(shù)空間中不規(guī)則點(diǎn)分布的簡單非光滑插值方法。

假設(shè)通過某種方法可以產(chǎn)生新的樣本,這些新樣本能夠局部地改變Voronoi單元的大小,從而改善NA面的局部分辨率。

1.2.2 目標(biāo)函數(shù)

目標(biāo)函數(shù)的建立是地球物理反演問題的關(guān)鍵,本文中目標(biāo)函數(shù)定義為測斜儀所記錄的觀測值與模型理論值之間的相對擬合差,其計(jì)算公式為

(1)

式中,N為觀測點(diǎn)個(gè)數(shù);Mobsi和Mcali分別為第i個(gè)觀測點(diǎn)處傾斜場的實(shí)際觀測值和理論計(jì)算值;Ci為權(quán)重系數(shù), 取值范圍為0～1, 其數(shù)值取決于觀測數(shù)據(jù)的質(zhì)量。

1.2.3 算法流程

采用鄰域算法進(jìn)行體積參數(shù)反演的具體步驟可歸納如下:

(1)首先在體積參數(shù)的合理區(qū)間范圍內(nèi)隨機(jī)產(chǎn)生Ns個(gè)初始體積模型M1={m1,m2,…,mNs}。

(2)計(jì)算每個(gè)體積模型對應(yīng)的目標(biāo)函數(shù)值。

(3)選擇出目標(biāo)函數(shù)最小的Nr個(gè)體積模型M2={m1,m2,…,mNr}。

(4)在M2決定的Voronoi單元內(nèi)再次隨機(jī)生成Ns個(gè)體積模型,即每個(gè)單元內(nèi)新生成的體積模型數(shù)目為Ns/Nr。

(5)依照上述方法重復(fù)迭代,每一次都選擇出Nr個(gè)最佳體積模型,最佳體積模型的Voronoi單元也隨之不斷更新改變。

(6)當(dāng)某次迭代產(chǎn)生的新體積模型滿足收斂條件時(shí),輸出該模型作為最終結(jié)果。

2 模型試算

2.1 傾斜場正演模擬

圖2為海域天然氣水合物水平井開采的簡化地質(zhì)模型。長度為300 m的水平井包含A、B兩個(gè)井段,分別用藍(lán)色和紅色表示。假定水平井沿東西向展布,以其中點(diǎn)在海底面上的投影為笛卡爾坐標(biāo)系的坐標(biāo)原點(diǎn),X軸、Y軸、Z軸分別對應(yīng)著正東、正北和豎直向上方向。忽略開采區(qū)域海底的地勢起伏,35個(gè)海底測斜儀(依次編號為1～35)以原點(diǎn)為中心呈網(wǎng)格狀均勻布設(shè),網(wǎng)格間距100 m。

圖2 海域天然氣水合物水平井開采簡化地質(zhì)模型Fig.2 Simplified geological model for gas hydrate horizontal well exploitation in sea area

以中國南海神狐海域的天然氣水合物儲層為例,該儲層深度范圍在201～278 m (海底以下深度),中心深度為240 m,儲層飽和度為33%、孔隙度為33%[8]。2017年南海天然氣水合物試采的總產(chǎn)氣量為3.09×106m3。標(biāo)準(zhǔn)狀態(tài)下,1 m3的固態(tài)水合物可分解為164 m3天然氣。若試采產(chǎn)出的天然氣完全來自于儲層中固態(tài)水合物的分解,那么儲層中分解的固態(tài)水合物體積約為1 884 m3。由于儲層中1 m3的固態(tài)水合物分解為0.2 m3天然氣和0.8 m3水,故此時(shí)造成的儲層體積變化為376.8 m3。根據(jù)理想氣體狀態(tài)方程,在神狐海域儲層原位溫度壓力狀態(tài)下,1 m3的游離氣在標(biāo)準(zhǔn)狀態(tài)下的產(chǎn)出的天然氣體積為130 m3。若試采產(chǎn)出的天然氣完全來自于儲層中的游離氣,此種情況下儲層的體積變化為2 377 m3。因此試采造成的儲層體積變化區(qū)間為376.8～2 377 m3。考慮到水平井開采能夠進(jìn)一步增大產(chǎn)能,本文中取模型的等效區(qū)間為1 000～4 000 m3。為便于對比模型體積對正演結(jié)果的影響,取A、B井段產(chǎn)氣量的等效體積分別為2 000和3 000 m3。若射孔段長度為50 m,則模型A的長、寬、高分別為50、6.3和6.3 m,模型B的長、寬、高分別為50、7.7和7.7 m。對應(yīng)的等效模型中心點(diǎn)坐標(biāo)(如圖黃色圓點(diǎn)所示)分別為(-100, 0,-240)和(100, 0,-240),即等效模型間距200 m。

圖3為A、B井段在各觀測點(diǎn)產(chǎn)生的傾斜場變化。

圖3 傾斜場正演模擬結(jié)果Fig.3 Forward modeling of tilt field

傾斜場矢量(用箭頭表示)是由東西向與南北向這兩個(gè)方向的傾斜場分量合成得到的,數(shù)值單位為微弧(μR),它指的是該位置處地層傾斜的角度值。黃色圓點(diǎn)表示水平井不同井段射孔段中心位置在海底面上的投影?？梢钥闯?不同井段在海底面產(chǎn)生的傾斜場以其中點(diǎn)在海底面的投影為中心,存在明顯的邊界區(qū)域。A、B井段產(chǎn)生的傾斜場集中分布于觀測系統(tǒng)西、東兩側(cè)。圖3(c)是A、B井段各自傾斜場的矢量疊加,黃色圓點(diǎn)表示觀測系統(tǒng)的中心位置。在其他參數(shù)一致的情況下,由于B井段相比A井段等效體積更大,對總傾斜場的貢獻(xiàn)更大,因此傾斜場表現(xiàn)為自觀測系統(tǒng)中心向東的偏移特征。

2.2 模型參數(shù)反演

為了排除體積參數(shù)選取的特殊性對反演結(jié)果的影響,在1 000～4 000 m3的體積區(qū)間內(nèi)隨機(jī)生成100組體積參數(shù)分別作為不同井段等效體積的真實(shí)值并正演模擬得到各觀測點(diǎn)處的傾斜場值作為反演的輸入值,再結(jié)合諸如射孔位置、儲層厚度、地層參數(shù)等已知信息,采用本文中提出的鄰域算法選取適當(dāng)參數(shù)(在本文中,鄰域算法的兩個(gè)調(diào)節(jié)參數(shù)Ns和Nr的取值分別為50和5)進(jìn)行地球物理反演得到模型的體積參數(shù)(反演解),如圖4所示。反演解與真實(shí)值高度重合,表明應(yīng)用鄰域算法進(jìn)行天然氣水合物儲層參數(shù)反演的可行性。

本文中采用高斯白噪聲探究隨機(jī)噪聲對反演結(jié)果的影響,它是一類概率密度函數(shù)服從高斯分布(即正態(tài)分布)的隨機(jī)噪聲,均值μ為0,標(biāo)準(zhǔn)差σ為1。在反演過程中,分別將5%、10%、15%這3組不同程度的隨機(jī)噪聲加入到真實(shí)值中作為實(shí)際觀測值進(jìn)行模型參數(shù)反演。圖5和6分別為不同噪聲水平下采用鄰域算法(NE)和直接網(wǎng)格搜索算法(GS)對100組體積參數(shù)進(jìn)行反演得到的反演解與真實(shí)值的對比及其誤差統(tǒng)計(jì)。其中圖5為從100組體積參數(shù)中以每5組為間隔選取20組反演結(jié)果得到的。圖6的縱坐標(biāo)及誤差棒分別代表了誤差的平均值及其標(biāo)準(zhǔn)差?？梢钥闯?反演效果隨著噪聲程度的增大而變差。噪聲程度達(dá)到15%時(shí),大部分反演解仍與真實(shí)值保持較高的吻合度,表明鄰域算法具有良好的抗噪性。兩種反演方法的誤差基本相等,而直接網(wǎng)格搜索算法的反演誤差稍大(圖6(d))的原因在于該算法采用等網(wǎng)格剖分,其搜索結(jié)果僅位于有限劃分的網(wǎng)格點(diǎn)上。

上述結(jié)果是利用處理器為i5-8250U、主頻為1.6 GHz的筆記本電腦計(jì)算得到的,鄰域算法反演的平均計(jì)算耗時(shí)約37 s,而直接網(wǎng)格搜索算法的平均計(jì)算耗時(shí)約293 s。雖然理論上當(dāng)網(wǎng)格搜索算法的網(wǎng)格剖分足夠密集時(shí),可以搜索到全局最優(yōu)解,但其計(jì)算耗時(shí)巨大,實(shí)時(shí)性差。鄰域算法利用的是模型目標(biāo)函數(shù)值排序,因而只需將目標(biāo)函數(shù)的相對值計(jì)算出來即可,這一做法大大降低了計(jì)算難度,提高了計(jì)算效率,節(jié)省計(jì)算耗時(shí)。

圖4 無噪聲條件下A、B井段的體積參數(shù)反演結(jié)果Fig.4 Volume parameters inversion results of A and B sections without noise

圖5 不同噪聲水平、不同算法下A、B井段的體積參數(shù)反演結(jié)果Fig.5 Volume parameters inversion results of A and B sections by different algorithms at different noise levels

為定量探究不同噪聲水平下模型間距對反演結(jié)果的影響,當(dāng)水平井劃分為兩個(gè)井段時(shí),在20～280 m的模型間距范圍內(nèi)取步長為20 m,在隨機(jī)生成的100組模型體積參數(shù)(真實(shí)值)正演得到的傾斜場值中加入不同程度的隨機(jī)噪聲,接著通過不同反演算法進(jìn)行體積參數(shù)反演并統(tǒng)計(jì)分析反演結(jié)果。圖7為無噪聲時(shí)分別運(yùn)用鄰域算法和模擬退火算法(SA)進(jìn)行體積參數(shù)反演的結(jié)果誤差統(tǒng)計(jì)。可以看出,當(dāng)?shù)刃Ｐ烷g距約大于80 m時(shí),兩種方法的反演誤差近乎相等,但當(dāng)?shù)刃Ｐ烷g距較小(小于80 m)時(shí),利用模擬退火法進(jìn)行參數(shù)反演的結(jié)果誤差較大,反演結(jié)果與真實(shí)值相偏離,而此時(shí),利用鄰域算法進(jìn)行模型參數(shù)反演的結(jié)果與真實(shí)值更為接近。

圖6 不同噪聲水平、不同算法下A、B井段的體積參數(shù)反演誤差統(tǒng)計(jì)Fig.6 Error statistics of volume parameters inversion results of A and B sections by different algorithms at different noise levels

圖7 不同模型間距、不同算法下A、B井段的體積參數(shù)反演結(jié)果Fig.7 Volume parameters inversion results of A and B sections by different inversion algorithms at different model spacing

圖8為不同模型間距時(shí)反演的目標(biāo)函數(shù)分布。可以看到,反演的目標(biāo)函數(shù)呈“單峰”狀。當(dāng)模型間距為100 m時(shí),全局最優(yōu)解落在相對較小的收斂區(qū)域內(nèi)。此時(shí),采用模擬退火法和鄰域算法均可得到較好的反演效果。而當(dāng)模型間距為20 m時(shí),全局最優(yōu)解所在的收斂區(qū)域相對較大。此時(shí),隨著迭代次數(shù)的增加,鄰域算法可實(shí)現(xiàn)對目標(biāo)區(qū)域更高精度的搜索。因此相比模擬退火法,鄰域算法的反演結(jié)果更接近真實(shí)值。

圖9為不同噪聲水平、不同模型間距下,利用鄰域算法反演的A、B井段體積參數(shù)與真實(shí)值的誤差統(tǒng)計(jì)。可以發(fā)現(xiàn),反演誤差隨著噪聲程度的增大而增大,隨著模型間距的增大而減小,即反演誤差與噪聲程度呈正相關(guān)性,與模型間距呈負(fù)相關(guān)性。

圖8 不同模型間距下反演目標(biāo)函數(shù)對比Fig.8 Contrast diagram of inversion objective function at different model spacing

圖9 不同噪聲水平、不同模型間距下A、B井段的體積參數(shù)反演誤差統(tǒng)計(jì)Fig.9 Error statistics of volume parameters inversion results of A and B sections at different noise levels and different model spacing

當(dāng)長度為300 m的水平井包含自西向東展布的A、B、C 3個(gè)井段時(shí),等效模型中心坐標(biāo)分別為(-100,0,-240)、(0,0,-240)和(100,0,-240),等效模型間距100 m。此時(shí)對隨機(jī)生成的100組體積參數(shù)在不同噪聲水平下運(yùn)用鄰域算法進(jìn)行反演,并統(tǒng)計(jì)分析反演結(jié)果(圖10)。同樣可以發(fā)現(xiàn),反演的體積誤差隨著噪聲程度的增大而增大。相比A和C井段,中間井段B的誤差均值與標(biāo)準(zhǔn)差都更大,這是由于B井段受到A和C井段的體積互補(bǔ)效應(yīng)。

當(dāng)水平井劃分為3個(gè)井段時(shí),在10～140 m的模型間距范圍內(nèi)取步長為10 m。圖11為無噪聲條件下運(yùn)用鄰域算法對不同模型間距下的體積參數(shù)進(jìn)行反演的結(jié)果誤差統(tǒng)計(jì)。圖12為不同噪聲水平不同模型間距下,采用鄰域算法所得體積參數(shù)的反演結(jié)果誤差統(tǒng)計(jì)。同樣可以發(fā)現(xiàn),反演誤差與噪聲程度呈正相關(guān)性,與模型間距呈負(fù)相關(guān)性。由于體積的互補(bǔ)特性,中間井段的反演誤差偏大。

圖10 不同噪聲水平下A、B、C井段的體積參數(shù)反演誤差統(tǒng)計(jì)Fig.10 Error statistics of volume parameters inversion results in A, B and C sections at different noise levels

圖11 不同模型間距下A、B、C井段的體積參數(shù)反演誤差統(tǒng)計(jì)Fig.11 Error statistics of volume parameters inversion results in A, B and C sections at different model spacing

圖12 不同噪聲水平不同模型間距下A、B、C井段反演的結(jié)果誤差統(tǒng)計(jì)Fig.12 Error statistics of inversion results of A, B and C sections at different noise levels and different model spacing

3 結(jié) 論

(1)與直接網(wǎng)格搜索法相比,提出的方法在保證收斂效果的前提下計(jì)算效率更高,實(shí)時(shí)性更強(qiáng);與模擬退火法相比,提出的方法在模型間距較小時(shí),反演精度更高。

(2)通過提出的方法進(jìn)行參數(shù)反演時(shí),反演誤差與噪聲程度呈正相關(guān)性,與模型間距呈負(fù)相關(guān)性。

(3)當(dāng)不同井段的產(chǎn)水和出砂量較小或占該段儲層變化體積的比例一致時(shí),提出的方法可為海域天然氣水合物水平井產(chǎn)氣量的動態(tài)評價(jià)提供算法基礎(chǔ)。