基于InSAR的澀北天然氣田地表形變監(jiān)測(cè)與開(kāi)采特征反演

謝亞奇 梁衛(wèi)濤 丁 武

1 陜西省一三一煤田地質(zhì)有限公司,陜西省韓城市象山路131院,715499

天然氣作為一種相對(duì)清潔的能源,其燃燒時(shí)產(chǎn)生的二氧化碳和有害物質(zhì)少于煤炭、石油等傳統(tǒng)能源,但天然氣長(zhǎng)期開(kāi)采會(huì)引起氣田地下壓力減小,導(dǎo)致地表沉降甚至誘發(fā)地震[1]。傳統(tǒng)的氣田地表沉降監(jiān)測(cè)方法包括水準(zhǔn)測(cè)量和導(dǎo)航衛(wèi)星測(cè)量,這2種技術(shù)受到觀測(cè)條件和成本限制,只能提供離散站點(diǎn)的形變信息,無(wú)法全面揭示氣田地表沉降狀態(tài)。InSAR技術(shù)具有空間分辨率高、覆蓋范圍大、可重復(fù)觀測(cè)和無(wú)需地面設(shè)備等優(yōu)點(diǎn),被廣泛應(yīng)用于地下水、油田、礦區(qū)、火山等由地下物質(zhì)變化引起的地表形變研究[2]。

澀北氣田是青海省最大的天然氣田,也是我國(guó)西部地區(qū)重要的天然氣基地之一。經(jīng)過(guò)幾十年開(kāi)采,已經(jīng)引發(fā)嚴(yán)重的地表沉降,可能導(dǎo)致采氣設(shè)備、輸氣管道和周?chē)ㄖ锏膿p壞。因此,需要對(duì)澀北氣田地表進(jìn)行定期監(jiān)測(cè),并根據(jù)監(jiān)測(cè)結(jié)果采取必要措施,確保設(shè)備安全運(yùn)行。

為監(jiān)測(cè)澀北氣田天然氣開(kāi)采引起的地表形變,本文使用時(shí)序InSAR技術(shù),基于退役和在軌的多顆SAR衛(wèi)星數(shù)據(jù),恢復(fù)澀北氣田近22 a開(kāi)采過(guò)程的地表形變特征,并討論天然氣開(kāi)采導(dǎo)致地表發(fā)生形變的機(jī)理。

1 研究區(qū)概況

柴達(dá)木盆地蘊(yùn)含豐富的油氣資源,目前已發(fā)現(xiàn)的天然氣儲(chǔ)量主要集中在盆地東部的三湖地區(qū)[3]。區(qū)內(nèi)目前已探明澀北1號(hào)、澀北2號(hào)和臺(tái)南等3個(gè)整裝大型第四系生物氣田及駝峰山、臺(tái)吉乃爾、伊克雅烏汝等含氣構(gòu)造[4],具體分布如圖1所示。

圖中綠框?yàn)闈睔馓锏貐^(qū)SAR數(shù)據(jù)覆蓋范圍圖1 柴東三湖地區(qū)氣田及背斜構(gòu)造帶分布[5]Fig.1 Distribution of gas fields and anticline structural belts in Sanhu region

澀北氣田開(kāi)采范圍內(nèi)地層主要為第四系湖泊相沉積,儲(chǔ)層巖性主要為泥質(zhì)粉砂巖,其次為粉砂質(zhì)泥巖、含粉砂泥巖和泥巖、粉砂巖[6]。澀北氣田主要由澀北1號(hào)、澀北2號(hào)和臺(tái)南3大氣田組成,其中澀北1號(hào)氣田發(fā)現(xiàn)于1964年,于1995年試采開(kāi)發(fā),1996年柴達(dá)木盆地東部首條天然氣長(zhǎng)輸管道投入運(yùn)營(yíng),標(biāo)志著澀北天然氣田正式開(kāi)發(fā)。由于持續(xù)的大規(guī)模天然氣開(kāi)采,澀北天然氣田已產(chǎn)生顯著的地表形變,這對(duì)氣田開(kāi)采設(shè)備和天然氣管道的安全運(yùn)行形成巨大威脅。

2 InSAR數(shù)據(jù)處理與形變分析

2.1 SAR數(shù)據(jù)與InSAR預(yù)處理

為研究澀北氣田天然氣開(kāi)采引起的地表形變,本文收集5個(gè)軌道的ERS-1/2、ENVISAT、ALOS和Sentinel-1A/B等多顆衛(wèi)星SAR數(shù)據(jù),時(shí)間跨度為1996～2018年,共168景數(shù)據(jù),詳細(xì)信息見(jiàn)表1。

表1 SAR數(shù)據(jù)信息Tab.1 SAR data Information

由于研究區(qū)氣候干燥、植被覆蓋非常少,2幅相隔數(shù)年的SAR影像仍有可能保持相干性?；诖?本文在選取InSAR干涉對(duì)時(shí),基于短空間基線和長(zhǎng)時(shí)間基線原則生成1 030幅干涉圖,從中選取597幅高質(zhì)量干涉圖進(jìn)行后續(xù)時(shí)間序列分析(圖2)。使用DEOS/DORIS、OPOD等精密軌道數(shù)據(jù)精化SAR數(shù)據(jù)的軌道參數(shù),使用30 m分辨率的SRTM DEM數(shù)據(jù)模擬并去除地形相位。為提高干涉圖的相干性,對(duì)各軌道的干涉圖進(jìn)行多視處理,使干涉圖分辨率與DEM分辨率相同。為進(jìn)一步提高干涉圖的信噪比,對(duì)干涉圖進(jìn)行自適應(yīng)空間域?yàn)V波處理,相位解纏使用基于Delaunay三角網(wǎng)的最小費(fèi)用流算法。

2.2 時(shí)序InSAR處理

為恢復(fù)澀北氣田自正式開(kāi)發(fā)以來(lái)地表形變隨時(shí)間的演化特征,本文采用小基線集技術(shù)(small baseline subset, SBAS)求解每個(gè)軌道數(shù)據(jù)的形變時(shí)間序列[7]。SBAS算法通過(guò)設(shè)置較小的時(shí)空基線閾值獲取高質(zhì)量干涉圖,以降低時(shí)空失相干和DEM誤差的影響。此外,通過(guò)時(shí)間域高通濾波和空間域低通濾波來(lái)估計(jì)軌道誤差和大氣延遲等空間相關(guān)誤差,從而獲取高精度的形變速率圖(圖3)和形變時(shí)間序列(圖4)。

圖4 點(diǎn)A在雷達(dá)視線向形變時(shí)間序列Fig.4 Time series of LOS deformation of point A

SAR側(cè)視成像僅能獲取地表形變?cè)诟鬈壍览走_(dá)視線向的投影信息,而聯(lián)合同一時(shí)期的升軌和降軌SAR數(shù)據(jù),不僅能提高監(jiān)測(cè)數(shù)據(jù)的時(shí)間分辨率,也可為獲取監(jiān)測(cè)區(qū)的多維形變信息提供基礎(chǔ)數(shù)據(jù)?？紤]到InSAR技術(shù)對(duì)南北方向形變不敏感,本文采用多維小基線集技術(shù)(multidimensional small baseline subset, MSBAS)[8]對(duì)2014年以來(lái)收集的升降軌Sentinel-1數(shù)據(jù)進(jìn)行東西方向和垂直方向的二維形變分解。MSBAS技術(shù)是在SBAS-InSAR技術(shù)基礎(chǔ)上發(fā)展而來(lái),利用在時(shí)間和空間上重疊的不同軌道升降軌SAR影像數(shù)據(jù),實(shí)現(xiàn)多視角影像時(shí)間分辨率的融合,解算研究區(qū)在垂直向和東西向的二維形變速率及累積形變時(shí)間序列。二維時(shí)間序列分解公式如下:

(1)

式中,A為SAR影像的時(shí)間間隔矩陣;S={SE,SU}={-cosθsinφ,cosφ}為投影矩陣,θ為方位向角度,φ為視線入射角;λ為正則化參數(shù);VE和VU為待求的EW向和UD向分量;D為D-InSAR技術(shù)獲取的視線向觀測(cè)值;L為正則化系數(shù)矩陣。本文利用MSBAS技術(shù)成功獲取澀北氣田地區(qū)二維形變速率(圖5)和形變時(shí)間序列(圖6),可為后續(xù)的氣田開(kāi)采時(shí)間序列參數(shù)反演提供形變數(shù)據(jù)。

圖5 MSBAS二維平均速率Fig.5 2D average velocity from MSBAS

圖6 點(diǎn)A在水平方向(東西向)和垂直方向形變時(shí)間序列Fig.6 Deformation time series of point A in horizontal and vertical directions

2.3 InSAR形變結(jié)果分析

利用1996～2018年5個(gè)軌道的SAR數(shù)據(jù)獲取澀北3大氣田的形變速率圖(圖3)和形變時(shí)間序列結(jié)果。由圖3可見(jiàn),天然氣開(kāi)采產(chǎn)生的地表形變主要分布在3大氣田開(kāi)采區(qū)內(nèi),其空間分布與已探明的澀北氣田氣藏空間分布(圖1)相似。T362、T269和T490軌道形變結(jié)果顯示,1996～2010年澀北1號(hào)氣田是3大氣田中形變量和形變范圍最大的地區(qū)。而T70和T77軌道形變結(jié)果顯示,2014～2018年臺(tái)南氣田已經(jīng)成為該地區(qū)形變量和形變范圍最大的地區(qū),說(shuō)明從2014年開(kāi)始,臺(tái)南氣田開(kāi)采速度明顯增加。從5個(gè)軌道的形變速率圖可以看出,澀北3大氣田形變場(chǎng)位置分布存在一定差異,即形變場(chǎng)在不同SAR成像幾何中的形變位置不同,說(shuō)明天然氣開(kāi)采產(chǎn)生的地表形變存在明顯的水平分量,而不單是垂直方向形變。

為分析澀北1號(hào)氣田從正式開(kāi)發(fā)至2018年約22 a的地表形變信息,提取澀北1號(hào)氣田沉降中心A點(diǎn)的形變時(shí)間序列結(jié)果(圖4)。InSAR監(jiān)測(cè)結(jié)果顯示,2003年之前A點(diǎn)基本保持穩(wěn)定,僅有微弱隆升;2003～2005年A點(diǎn)隆升近6 cm,由于這2 a多時(shí)間缺少觀測(cè)數(shù)據(jù),尚不清楚其具體形變過(guò)程;2005年A點(diǎn)開(kāi)始快速下沉,至2018年仍保持下沉趨勢(shì),且下沉速率明顯增加。

由圖5可見(jiàn),垂直向形變量約占雷達(dá)視線向形變量的80%。由圖6可見(jiàn),澀北1號(hào)氣田地表形變中心在垂直方向和東西方向基本表現(xiàn)為線性形變,且垂直方向形變量遠(yuǎn)大于水平方向形變量。

3 大地測(cè)量數(shù)據(jù)建模

為分析澀北1號(hào)氣田因天然氣開(kāi)采導(dǎo)致地表產(chǎn)生形變的機(jī)理,本文基于均勻彈性半空間位錯(cuò)模型,采用GBIS軟件[9]分析澀北1號(hào)氣田在2015～2018年間形變?cè)醋兓^(guò)程。該軟件采用馬爾科夫鏈蒙特卡洛(MCMC)方法,結(jié)合Metropolis-Hasting算法對(duì)收斂于貝葉斯后驗(yàn)分布的模型參數(shù)進(jìn)行隨機(jī)抽樣,計(jì)算每個(gè)參數(shù)概率密度函數(shù)的估計(jì)值,并提供參數(shù)在95%置信區(qū)間的最優(yōu)解。

本文選擇GBIS軟件中Sill模型,利用垂直方向和東西方向形變時(shí)間序列結(jié)果,計(jì)算澀北1號(hào)氣田在2015～2018年間形變?cè)吹拈L(zhǎng)度、寬度、深度、走向、參考點(diǎn)位置(X,Y)和張量等參數(shù)。圖7為利用升降軌Sentinel-1衛(wèi)星垂直方向和東西方向累積形變反演得到的2018-03-11模擬形變和殘差形變,圖8為反演得到的澀北1號(hào)氣田形變?cè)瓷疃群腕w積隨時(shí)間變化曲線。

圖7 澀北1號(hào)氣田Sentinel-1衛(wèi)星垂直方向和東西方向累積形變模擬結(jié)果Fig.7 Cumulative displacement in horizontal and vertical directions from Sentinel-1 satellites of Sebei No.1 gas field

圖8 基于二維形變時(shí)間序列反演Sill模型地下開(kāi)采體積和深度變化Fig.8 Estimated underground gas exploitation volume and depth change of Sill model using the time series of 2D deformation

由圖7可見(jiàn),Sentinel-1衛(wèi)星垂直方向和東西方向累積形變反演結(jié)果的模擬殘差均很小,說(shuō)明該反演方法較為可靠。由圖8可見(jiàn),形變?cè)吹奶烊粴獾叵麦w積在2015～2018年逐漸增大,即天然氣地下開(kāi)采體積呈近似線性增加趨勢(shì)。形變?cè)瓷疃茸兓瘎t分為兩個(gè)階段:2015-01～2016-08為臺(tái)階式變化,開(kāi)采深度增加至2 100 m;2016-08以后開(kāi)采深度逐步穩(wěn)定在2 300m左右。

4 討 論

4.1 2015～2017年澀北1號(hào)天然氣開(kāi)采量估計(jì)

本文通過(guò)InSAR形變結(jié)果反演得到2015-01～2018-03澀北1號(hào)天然氣田地下氣體體積減少約6 000萬(wàn)m3,2015～2017年地下氣體體積分別減少約1 646萬(wàn)m3、1 635萬(wàn)m3和19.63萬(wàn)m3。根據(jù)2016～2018年《中國(guó)石油青海油田公司年鑒》可知,2015～2017年澀北1號(hào)氣田實(shí)際氣體產(chǎn)量分別為13.72億m3、13.07億m3和15.98億m3[10-12]。在相同時(shí)間段內(nèi),InSAR反演得到的地下氣體體積與實(shí)際氣體產(chǎn)量之間存在巨大差異,這歸因于不同環(huán)境條件(如溫度和壓力)下氣體體積的變化。本文估計(jì)地下氣體體積變化的條件為澀北1號(hào)氣田2015～2017年平均壓力為920 kPa、平均溫度為319.26 K[10-12];而《中國(guó)石油青海油田公司年鑒》統(tǒng)計(jì)結(jié)果則是基于中國(guó)標(biāo)準(zhǔn)狀態(tài),即壓力為101.325 kPa、溫度為293.15 K。為增加評(píng)估結(jié)果的可靠性,根據(jù)理想氣體狀態(tài)方程,將地下條件下的氣體體積變化轉(zhuǎn)換到中國(guó)標(biāo)準(zhǔn)狀態(tài)。轉(zhuǎn)換后的氣體體積變化結(jié)果為13.72億m3、13.07億m3和15.98億m3,由此可見(jiàn),InSAR反演得到的天然氣體積與《中國(guó)石油青海油田公司年鑒》統(tǒng)計(jì)結(jié)果的差異極小,可以定量評(píng)估InSAR反演方法在監(jiān)測(cè)地下氣體體積變化方面的可靠性。

此外,以下因素也會(huì)導(dǎo)致InSAR估計(jì)的天然氣開(kāi)采量與《中國(guó)石油青海油田公司年鑒》統(tǒng)計(jì)結(jié)果之間存在差異:1)本文使用均勻反演模型,反演過(guò)程中僅采用單一形變?cè)吹钠骄疃群推骄鶑埩坑?jì)算地下氣體體積,由于氣藏區(qū)層數(shù)多,層深和層厚不等,因此會(huì)產(chǎn)生一定誤差;2)天然氣開(kāi)采通常伴隨有少量水和砂,會(huì)造成一定差異;3)根據(jù)孔隙彈性理論,地下氣體壓力降低可能會(huì)抵消一定的地下體積變化[13]。

4.2 天然氣開(kāi)采過(guò)程中的地表形變模式

由于天然氣的大量開(kāi)采,地下物質(zhì)逐漸減少,導(dǎo)致地表產(chǎn)生明顯形變。圖9為天然氣開(kāi)采剖面示意圖,從中可看出天然氣開(kāi)采與地表形變之間的關(guān)系,紅線表示由于地下天然氣開(kāi)采產(chǎn)生的正斷層,黃色箭頭表示地表水平形變方向。由于地下物質(zhì)開(kāi)采導(dǎo)致地面下沉,形變中心周?chē)鷧^(qū)域地表也朝著形變中心發(fā)生水平位移。從圖5和圖6可以看出,澀北3大氣田形變區(qū)在快速下沉的同時(shí)還伴隨有明顯的水平形變。

圖9 天然氣開(kāi)采剖面示意圖[14]Fig.9 Section diagram of gas production

由于SAR為側(cè)視成像,這會(huì)導(dǎo)致單一軌道InSAR觀測(cè)對(duì)水平形變的敏感度較低。本文利用同一時(shí)期的升降軌SAR數(shù)據(jù)進(jìn)行聯(lián)合解算,獲得研究區(qū)二維(垂直方向和東西方向)形變時(shí)間序列(圖6),可更全面地展示天然氣開(kāi)采過(guò)程中地表形變變化過(guò)程。本文選用InSAR二維形變場(chǎng)作為氣田形變?cè)捶囱莸挠^測(cè)數(shù)據(jù),其反演精度優(yōu)于基于單一軌道InSAR獲得的形變場(chǎng)。

為分析澀北1號(hào)氣田的開(kāi)采過(guò)程,繪制該氣田1997～2013年日產(chǎn)氣量曲線圖(圖10)。由圖可知,2003年以前生產(chǎn)規(guī)模較小,2004年開(kāi)始大規(guī)模開(kāi)采,2007年以后開(kāi)采速度明顯提高。

圖10 澀北1號(hào)氣田開(kāi)發(fā)曲線圖[15]Fig.10 Development curve of Sebei No.1 gas field

結(jié)合InSAR時(shí)間序列形變結(jié)果(圖4)可以看出,1996～2003年澀北1號(hào)氣田地表形變不明顯,該時(shí)間段內(nèi)澀北1號(hào)氣田處于探索試采階段,開(kāi)采井?dāng)?shù)量和產(chǎn)氣量較少,因此未產(chǎn)生明顯的地表形變。隨著澀北1號(hào)氣田的不斷開(kāi)發(fā),部分老井地層壓力逐步下降,澀北1號(hào)氣田在2003～2006年平均單井日產(chǎn)氣量遞減明顯[16]。為提高老井天然氣產(chǎn)量,通過(guò)注氣增壓方式將底層水排出,這也導(dǎo)致2003～2005年間澀北1號(hào)氣田地表快速隆升。2005～2007年,澀北氣田加快產(chǎn)能建設(shè),引入水平井技術(shù),極大提高了天然氣產(chǎn)量。自2005年起,澀北1號(hào)氣田地表開(kāi)始出現(xiàn)快速下沉,最大形變區(qū)(圖3中A點(diǎn))在2005～2011年形變速率達(dá)到30 mm/a,2014～2018年形變速率更是達(dá)到70 mm/a。

5 結(jié) 語(yǔ)

1)二維InSAR形變時(shí)間序列表明,在天然氣開(kāi)采區(qū)域,地表以垂直形變?yōu)橹?、水平形變?yōu)檩o,而在東臺(tái)吉乃爾湖地區(qū)則表現(xiàn)為單純的垂直方向地面沉降。

2)基于彈性位錯(cuò)理論反演的澀北1號(hào)氣田開(kāi)采過(guò)程顯示,開(kāi)采初期開(kāi)采量和開(kāi)采深度成正比。該氣田開(kāi)采深度從2015年的1 500 m逐漸加深至2017年的2 300 m。

3)根據(jù)理想氣體狀態(tài)方程計(jì)算InSAR反演的地下天然氣儲(chǔ)層體積變化在標(biāo)準(zhǔn)狀態(tài)下的體積,結(jié)果顯示,InSAR反演結(jié)果與實(shí)際統(tǒng)計(jì)數(shù)據(jù)差異極小,表明InSAR監(jiān)測(cè)和反演結(jié)果具有可靠性。