孫偉 劉峻峰 高海英 劉新 宋盼盼
1.國家管網集團華中分公司 2.國家管網集團科技數字本部 3.紫光軟件系統有限公司
我國地勢多山地、高原及丘陵,而油氣管道線路較長,在鋪設過程中不可避免地會經過地形復雜區域,尤其是山區、丘陵等地形起伏較大的區域。油氣管道在運行過程中常常受到地質災害的威脅與破壞,一旦發生,不僅會造成管道變形、斷裂,導致油氣泄漏、管線停輸,帶來巨大的經濟損失,還有可能引發火災、爆炸等事故[1-4]。因此,為了盡可能地避免、減少地質災害對輸油管道的威脅與破壞,確保能源運輸的長效安全,應采取科學、有效的方法與手段對管道沿線發育或潛在的地質災害進行全面識別、監測與預警。
地質災害的發生發展具有一定的隱蔽性、潛伏性和突發性,GPS測量、北斗定位技術等監測手段雖然精度高,但只能獲取離散的監測點,且監測周期較長,不易發現隱蔽性較強的隱患災害點[5-7]。星載合成孔徑雷達干涉測量(InSAR)技術具有覆蓋范圍廣、監測精度高,可進行全天時、全天候監測地表形變的特點[8]。作為雷達遙感中的先進技術,采用時間序列雷達遙感圖像的干涉雷達和差分干涉雷達技術,可以精確地測定地面的微小位移變化。該技術已經在地面沉降監測、滑坡、地震、火山活動及其他地質災害評估與監測等方面得到了廣泛的應用[9]。2005年,吳立新等采用5景ERS1/2數據,對唐山市和開灤礦區使用差分干涉測量技術(D-InSAR)進行了監測,得到了間隔時間超過半年的地下采礦及采水導致的形變[10]。2008年,葛大慶等以廊坊市主城區為研究對象,介紹了短基線差分干涉紋圖集的地表形變場監測方法[11]。2018年,張路教授以大渡河上游丹巴縣為研究區域,采用時序InSAR技術與歷史累積的先進對地觀測衛星(ALOS PALSAR)和先進的合成孔徑雷達(ASAR)數據成功地識別出17處不穩定斜坡,通過與實測數據對比,驗證了InSAR技術的有效性和優勢[12]。2019年,何朝陽等通過永久散射體合成孔徑雷達干涉測量(PS-InSAR)技術,利用四川省理縣47景SentineL-1A數據,成功獲取該區域3處形變體,形變結果與GPS數據顯示的變形趨勢一致,驗證了InSAR結果的合理性[13]。2020年,戴可人教授使用Sentinel-1A數據,利用時序InSAR技術對雅礱江流域雅江縣-木里縣段的高山峽谷區域進行了滑坡災害隱患廣域早期識別,成功探測到8處隱患區域[14]。同年,張磊以中緬油氣管道貴州省晴隆段為主要研究對象,使用InSAR時序分析技術全面識別地質災害早期隱患點位置,成功地識別到地表明顯變化圖斑41個,并對形變區危險等級進行劃分[15]。由此可以證明,InSAR技術可為地質災害防治工作提供重要的數據支撐。
目前,國內使用InSAR技術對輸油管道地表地質災害監測并結合水土流失問題導致地質災害發生的研究較少,因此本研究根據研究區內地形地貌特征,基于Sentinel-1A數據利用小基線集(SBAS)干涉疊加技術,對成品油管道進行地質災害識別研究,討論InSAR技術在輸油管道地質災害監測方面的適用性、監測成果的可靠性,并結合現場核查資料進一步分析地質災害對輸油管道安全運作的危害程度,以便更準確地為地質災害防治和水土保持措施提供數據支持。
Berardino和Lanari等在2002年首次提出SBAS-InSAR技術,基于多主影像的組合方式,通過設置時間基線和空間基線閾值將同一地區多期合成孔徑雷達(SAR)數據組合成短時間短基線的干涉對[16],從而通過干涉疊加算法來獲得研究區的形變速率和形變量。干涉對是利用D-InSAR技術得到的攜帶研究區域形變信息的相位差數據,該數據則是通過衛星在同一區域兩次成像所得到的SAR數據做相位差分而形成。本研究以重軌干涉測量的成像原理來介紹D-InSAR監測地表形變的原理。D-InSAR的幾何原理如圖1所示。

圖1中S1和S2分別代表不同時刻在不同空間位置的衛星傳感器,將S1設定為主影像成像傳感器,S2設定為輔影像成像傳感器;R1和R2分別代表在雷達視線向(LOS)上S1與S2與地表之間的距離;B代表S1與S2之間的距離,稱為空間基線,它在LOS向與S1垂直視線向的投影;α為空間基線B與水平方向之間的夾角,θ為主影像雷達入射角;h為傳感器與地表之間的距離,P為地表監測點,P’為P點發生形變后的位置。在不考慮大氣、噪聲等因素的影響時,S1和S2獲取的SAR影像進行差分干涉得到的干涉相位如式(1)。
(1)
式中:Δφ為相位差,(°);λ為雷達入射波波長,m;B⊥為垂直基線長度,m;h為傳感器與地表之間的距離,m;R1為S1傳感器與地面之間的斜距,m;θ為主影像雷達入射角,(°);Δd為地表形變量,mm。
式(1)等式右側前式為地形相位,后式為形變相位,那么去除地形相位后就可得到形變相位φdef,如式(2)。
(2)
獲取N幅SAR影像,生成M幅差分干涉圖(見圖2)。SBAS-InSAR技術利用這些干涉圖,估計出平均形變速率之后,通過時域的高通濾波和空間域的低通濾波去除大氣的影響,使用奇異值分解(SVD)方法估計形變時間序列。該方法可以處理多個子集,聯合同一傳感器獲取的不同時間段的干涉結果進行求解,提高結果的精度,降低數字高程模型(DEM)誤差和大氣相位延遲的干擾。SVD方法在很大程度上解決了由于時空基線過長而導致的失相干和大氣效應問題,且滿足兩個重要要求:①使用所有小基線子集組合的方式增加了時間采樣率;②得到了更為連續的形變圖,提高了形變圖空間采樣率[17-18]。

該技術數據處理流程如圖3所示。利用加入相應的精密軌道數據的若干景SAR數據進行干涉處理,剔除失相干較為嚴重的干涉對,使用干涉效果好的干涉結果進行時序處理,減小大氣效應和DEM數據帶來的誤差,得到研究區形變速率信息。

本研究選取的成品油管道區段位于江西省南部,該區段共包含318個樁,管道長度16.56 km,起于江西省吉安市萬安縣,止于贛州市的贛縣。研究區段范圍內海拔高度最低38 m,最高590 m,相對高差552 m,總體呈現出東南、西南、西北部地勢偏高,主要地形以中低山、丘陵為主,北部和中部為低丘崗地。管線所經區域,地表主要為第四系砂土、亞砂土等組成的黃色土壤,基巖由白堊系粉砂巖、礫巖等組成。管線經過山體一側形成順向坡,或者順坡上行,形成橫向坡。管道鋪設和施工道路等形成人工切坡,地表覆蓋層薄,植被恢復較差,巖土體裸露現象較多,附著能力差,在雨水的沖刷下,水土流失現象明顯,導致易發生坡體崩塌和滑坡等地質災害。
共計使用46景C波段Sentinel-1A升軌SAR數據,時間覆蓋范圍為2019年1月8日-2020年7月13日,具體信息見表1。SentineL-1A衛星數據的成像模式是干涉寬幅模式,在該模式下衛星采用TOPS技術,該技術的一大特點是以80 km寬,20 km長的脈沖帶為成像基礎,在方位向相鄰脈沖帶之間存在1.5 km的重疊區域,在3個子條帶之間存在2 km的重疊區域,從而保證數據對地面覆蓋的連續性[19-20]。表1所列為Sentinel-1 A影像數據主要參數。

表1 Sentinel-1A影像數據主要參數波段波長/cm軌道方向空間分辨率/m幅寬/km入射角/(°)C5.6升軌5×2025039.386
根據研究區特點及數據優勢,時間基線閾值設置為48天,空間基線閾值設置為300 m,去除相干性較差的干涉對后,參與時序InSAR處理的干涉對共計164對。軌道數據采用成像21天之后發布的POD精密軌道數據,DEM數據使用的是30 m的SRTM DEM數據。根據InSAR監測結果進行現場實地調查,獲取到隱患災害點的現場形變信息。
使用架構于ENVI遙感圖像處理軟件之上的SARScape軟件進行SBAS-InSAR處理,主要包括數據讀取、影像裁剪、干涉處理、去平地效應、自適應濾波、相位解纏,以及第1次估算平均位移速率和DEM校正系數、第2次估算平均位移速率和DEM校正系數,最后獲得研究區時間序列形變和平均形變速率,對形變速率異常值進行了剔除,得到了更加準確的研究區形變速率圖。
采用SBAS-InSAR技術對輸油管道重點區段地表形變進行監測,得到研究區平均形變速率,需要說明的是,監測到的平均形變速率均為LOS方向的形變值。圖4所示為采用InSAR技術獲取的輸油管道沿線區域平均形變速率及管線地表的坡度、坡向數據。從圖4(a)可看出,形變速率圖中存在多處空洞,致使空洞區域無法提取形變信息。從3個方面分析其原因:①由于InSAR技術本身是以不同時期相同區域SAR數據之間的相干性為基礎來監測地表形變,在其執行過程中,會設置相干性閾值,本研究相干性閾值為0.3,即兩景SAR數據同一位置的像素相干性大于0.3時,參與SBAS-InSAR干涉疊加分析,小于0.3的像素則會被剔除,而被剔除的區域就會在結果中以空值形式顯示;②當SAR入射角小于地表坡度,SAR脈沖照射到位于傳感器掃描方向一致的斜坡時,脈沖到達斜坡頂部和底部的斜距差小于實際地面之間的距離,致使SAR影像上的斜面長度被縮短了,出現了“透視收縮”[21],通過疊加分析形變速率和坡度數據可以知道,該區域最大坡度為37.65°,而本研究使用的SAR衛星影像入射角為39.386°,因此一定不是“透視收縮”引起的形變速率空洞;③使用的Sentinel-1A衛星運行軌道方向為升軌,即衛星繞地球由南偏東向北偏西飛行,傳感器掃描方向與飛行方向一致,那么Sentinel-1A對與其相垂直的山體的敏感程度較弱,且在坡度較大時易在SAR影像形成“陰影”。分析存在空洞的區域,由于坡度偏小,南北和東西方向山體均勻分布,不足以將形變空洞歸結于此。

通過上述分析可知,形變速率存在空洞的主要原因是由于時間序列SAR影像的不斷變化造成影像失相干,在數據運算過程中失相干區域被剔除,從而無形變結果。
根據試驗區形變速率結果對輸油管道全線形變情況以及影響范圍進行討論分析。受InSAR技術自身的限制,形變監測結果存在一定誤差,在實際分析形變結果時會將這部分誤差剔除掉,即形變速率在誤差絕對值范圍內浮動,將其視為未發生地表形變。利用全線路空間域形變速率直方圖統計結果(見圖5),設置誤差允許值(ε=-7.3 mm),將該值作為判斷是否發生地表形變的固定閾值。為了準確地判斷地表是否發生形變,將形變速率在一定長度范圍內的一階導數不等于零作為限定條件,限定條件為:v≥-7.3 mm,v/(p·20 m)≠0,其中v為形變速率,p為像素個數。當形變速率小于-7.3 mm,且形變速率在一定長度范圍內的一階導數不等于零時,確定該區域發生形變。

利用上述限定條件提取全線地表形變情況,共獲取到4處形變區域,標記為A1、A2、A3、A4。其中,A2形變區域發生在斜坡上(見圖6),其他均發生在平原區域,因此重點對A2區域進行形變分析。為了能更加清晰地反應形變對管線的影像程度和范圍,根據時序InSAR技術獲取的平均形變速率圖,繪制了該不穩定斜坡東西方向和南北方向的平均形變速率剖面曲線。東西方向沿斜坡蠕動方向進行繪制(見圖6(a)中的PP’線段),南北方向沿管道上行方向進行繪制(見圖6(a)中的NN’線段),東西向和南北向平均形變速率剖面曲線如圖6(b)、圖6(c)所示。該形變區主要覆蓋管道東側斜坡,該斜坡最大高程298 m、最小174 m、相對高差124 m、坡面長度376 m、坡度約32°。
由形變速率圖和剖面曲線圖可知,該形變區涉及范圍0.02 km2,東西向寬度320 m,南北向長度144 m,也就是說輸油管道在該位置144 m內存在安全隱患,從東西向剖面圖中可以看出(見圖6(b)),輸油管道兩側存在兩處沉降中心,東側斜坡平均形變速率最大達到-43 mm/a,西側平坦地表最大形變速率達到-35 mm/a,輸油管道穿越位置約-30 mm/a。依據Colesanti C教授的經驗總結[22],對于C波段SAR數據,當斜坡的LOS向形變速率絕對值大于2 mm/a時,可認為斜坡處于不穩定狀態。顯然,該斜坡已處于不穩定狀態,存在滑坡的可能,該監測結果與現場監測成果進行核驗,發現形變區位置及范圍與實地調查結果基本一致。

為了研究在輸油管道附近監測到的4個形變區域在時間序列上的形變過程,利用SBAS-InSAR技術獲取的時間序列累積形變數據繪制了4個形變區形變中心點在2019年1月-2020年7月的累積形變曲線圖(見圖7)。考慮到斜坡失穩及水土流失等地質災害受雨水影響嚴重[23-26],故收集并統計了監測時間段內研究區降雨情況。結合降雨情況,對監測到的4個區域形變過程及成因進行分析說明。
從圖7中可以看出,4個區域從2019年2月開始,均發生緩慢形變。其中,A1、A2、A3在2019年3月開始發生形變,A4在2019年2月開始發生沉降,至2019年8月趨于平緩,累積形變值分別達到-36 mm、-44 mm、-59 mm、-53 mm,在之后近4個月內形變較為穩定。根據降雨情況也可以看出,研究區在2019年1月-2019年7月,降雨天數均在10天以上,降雨量也隨之增大,極有可能是增加斜坡和易積水區域的形變量和形變速率的主要因素。2020年1月,A2斜坡、A3和A4平原區沉降速率加劇,分別在2020年3月中旬和2020年5月初達到沉降谷值,其中A2斜坡的形變值達到-61 mm,相比去年8月增加25 mm。根據降雨情況可以推斷,由于從2020年1月開始到2020年6月結束,研究區降雨天數增加后減少,期間存在連續降雨情況,斜坡以及斜坡下水溝經過長時間雨水浸泡,坡面和溝壁土質軟化,從而形變增加;A3和A4平原區形變值分別達到-79 mm、-68 mm,該區域為平原,因此可判斷地質災害類型主要是沉降和水土流失。2020年2月,A1緩坡區形變速率增加,至2020年5月達到谷值-74 mm,之后形變區域平緩,與降雨天數成正相關。監測時間段內形變速率的變化與月降雨天數之間的相關性達到0.69。由此可以推斷,地表水土流失、斜坡失穩的主要形成原因是降雨,降雨對地表和斜坡的沖擊將不斷使之發生不均勻形變。因此,在每年雨季應及時對輸油管道周邊易積水、斜坡不穩定的區域加強防范,避免發生重大地質災害。

將SBAS-InSAR技術應用于輸油管道重點區段的地質災害監測,使用46景覆蓋輸油管道的Sentinel-1A數據獲取該區域時間序列形變信息。結合坡度、坡向數據對形變速率結果存在空洞進行分析,結果表明:
(1) 不同時間段的SAR影像存在失相干是造成最終空洞的主要原因。
(2) 當山體坡度大于SAR傳感器入射角且坡向與衛星運行軌道方向一致時,易出現透視收縮和陰影,而使最終的形變速率結果出現空洞。
因此,在后續的應用當中,應結合穿透能力更強的SAR數據進行輸油管道地表形變監測。通過分析由SBAS技術獲取的形變速率值,獲取全線路發生形變區域有4處,并對其中一處發生在斜坡的形變區域進行了形變速率和影響范圍的分析,發現該斜坡最大年平均形變速率達到-43 mm/a,對管道及周邊0.02 km2造成影響。最后引入降雨天數數據,分析這4處形變區域與降雨之間的關系,結果表明降雨增多是致使地表發生形變的原因之一。
InSAR技術作為地質災害監測的新興技術手段,能夠在廣域、短周期內對地表進行形變監測,及時發現地表形變異常,能夠為輸油管道地質災害防治提供重要數據支撐。不可忽略的是,InSAR技術是以同一區域不同時間點的SAR之間的相干性作為算法實現的基本條件。因此,在地表相干性較差時會出現形變速率空值,無法獲取形變的情況,這是InSAR存在的不足,也是目前眾多學者不斷探討的問題。后續的實驗過程中,將嘗試降低相干性閾值,降低空洞概率,或使用波長更長的SAR數據對失相干區域進行二次監測。