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

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

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

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

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

1 研究區(qū)概況

柴達(dá)木盆地蘊(yùn)含豐富的油氣資源,目前已發(fā)現(xiàn)的天然氣儲量主要集中在盆地東部的三湖地區(qū)[3]。區(qū)內(nèi)目前已探明澀北1號、澀北2號和臺南等3個整裝大型第四系生物氣田及駝峰山、臺吉乃爾、伊克雅烏汝等含氣構(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

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

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

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

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

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

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

2.2 時序InSAR處理

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

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

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

(1)

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

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

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

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

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

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

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

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

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

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

圖7 澀北1號氣田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 基于二維形變時間序列反演Sill模型地下開采體積和深度變化Fig.8 Estimated underground gas exploitation volume and depth change of Sill model using the time series of 2D deformation

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

4 討 論

4.1 2015～2017年澀北1號天然氣開采量估計(jì)

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

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

4.2 天然氣開采過程中的地表形變模式

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

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

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

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

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

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

5 結(jié) 語

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

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

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