孫偉 劉峻峰 高海英 劉新 宋盼盼
1.國家管網(wǎng)集團(tuán)華中分公司 2.國家管網(wǎng)集團(tuán)科技數(shù)字本部 3.紫光軟件系統(tǒng)有限公司
我國地勢多山地、高原及丘陵,而油氣管道線路較長,在鋪設(shè)過程中不可避免地會經(jīng)過地形復(fù)雜區(qū)域,尤其是山區(qū)、丘陵等地形起伏較大的區(qū)域。油氣管道在運(yùn)行過程中常常受到地質(zhì)災(zāi)害的威脅與破壞,一旦發(fā)生,不僅會造成管道變形、斷裂,導(dǎo)致油氣泄漏、管線停輸,帶來巨大的經(jīng)濟(jì)損失,還有可能引發(fā)火災(zāi)、爆炸等事故[1-4]。因此,為了盡可能地避免、減少地質(zhì)災(zāi)害對輸油管道的威脅與破壞,確保能源運(yùn)輸?shù)拈L效安全,應(yīng)采取科學(xué)、有效的方法與手段對管道沿線發(fā)育或潛在的地質(zhì)災(zāi)害進(jìn)行全面識別、監(jiān)測與預(yù)警。
地質(zhì)災(zāi)害的發(fā)生發(fā)展具有一定的隱蔽性、潛伏性和突發(fā)性,GPS測量、北斗定位技術(shù)等監(jiān)測手段雖然精度高,但只能獲取離散的監(jiān)測點(diǎn),且監(jiān)測周期較長,不易發(fā)現(xiàn)隱蔽性較強(qiáng)的隱患災(zāi)害點(diǎn)[5-7]。星載合成孔徑雷達(dá)干涉測量(InSAR)技術(shù)具有覆蓋范圍廣、監(jiān)測精度高,可進(jìn)行全天時(shí)、全天候監(jiān)測地表形變的特點(diǎn)[8]。作為雷達(dá)遙感中的先進(jìn)技術(shù),采用時(shí)間序列雷達(dá)遙感圖像的干涉雷達(dá)和差分干涉雷達(dá)技術(shù),可以精確地測定地面的微小位移變化。該技術(shù)已經(jīng)在地面沉降監(jiān)測、滑坡、地震、火山活動及其他地質(zhì)災(zāi)害評估與監(jiān)測等方面得到了廣泛的應(yīng)用[9]。2005年,吳立新等采用5景ERS1/2數(shù)據(jù),對唐山市和開灤礦區(qū)使用差分干涉測量技術(shù)(D-InSAR)進(jìn)行了監(jiān)測,得到了間隔時(shí)間超過半年的地下采礦及采水導(dǎo)致的形變[10]。2008年,葛大慶等以廊坊市主城區(qū)為研究對象,介紹了短基線差分干涉紋圖集的地表形變場監(jiān)測方法[11]。2018年,張路教授以大渡河上游丹巴縣為研究區(qū)域,采用時(shí)序InSAR技術(shù)與歷史累積的先進(jìn)對地觀測衛(wèi)星(ALOS PALSAR)和先進(jìn)的合成孔徑雷達(dá)(ASAR)數(shù)據(jù)成功地識別出17處不穩(wěn)定斜坡,通過與實(shí)測數(shù)據(jù)對比,驗(yàn)證了InSAR技術(shù)的有效性和優(yōu)勢[12]。2019年,何朝陽等通過永久散射體合成孔徑雷達(dá)干涉測量(PS-InSAR)技術(shù),利用四川省理縣47景SentineL-1A數(shù)據(jù),成功獲取該區(qū)域3處形變體,形變結(jié)果與GPS數(shù)據(jù)顯示的變形趨勢一致,驗(yàn)證了InSAR結(jié)果的合理性[13]。2020年,戴可人教授使用Sentinel-1A數(shù)據(jù),利用時(shí)序InSAR技術(shù)對雅礱江流域雅江縣-木里縣段的高山峽谷區(qū)域進(jìn)行了滑坡災(zāi)害隱患廣域早期識別,成功探測到8處隱患區(qū)域[14]。同年,張磊以中緬油氣管道貴州省晴隆段為主要研究對象,使用InSAR時(shí)序分析技術(shù)全面識別地質(zhì)災(zāi)害早期隱患點(diǎn)位置,成功地識別到地表明顯變化圖斑41個,并對形變區(qū)危險(xiǎn)等級進(jìn)行劃分[15]。由此可以證明,InSAR技術(shù)可為地質(zhì)災(zāi)害防治工作提供重要的數(shù)據(jù)支撐。
目前,國內(nèi)使用InSAR技術(shù)對輸油管道地表地質(zhì)災(zāi)害監(jiān)測并結(jié)合水土流失問題導(dǎo)致地質(zhì)災(zāi)害發(fā)生的研究較少,因此本研究根據(jù)研究區(qū)內(nèi)地形地貌特征,基于Sentinel-1A數(shù)據(jù)利用小基線集(SBAS)干涉疊加技術(shù),對成品油管道進(jìn)行地質(zhì)災(zāi)害識別研究,討論InSAR技術(shù)在輸油管道地質(zhì)災(zāi)害監(jiān)測方面的適用性、監(jiān)測成果的可靠性,并結(jié)合現(xiàn)場核查資料進(jìn)一步分析地質(zhì)災(zāi)害對輸油管道安全運(yùn)作的危害程度,以便更準(zhǔn)確地為地質(zhì)災(zāi)害防治和水土保持措施提供數(shù)據(jù)支持。
Berardino和Lanari等在2002年首次提出SBAS-InSAR技術(shù),基于多主影像的組合方式,通過設(shè)置時(shí)間基線和空間基線閾值將同一地區(qū)多期合成孔徑雷達(dá)(SAR)數(shù)據(jù)組合成短時(shí)間短基線的干涉對[16],從而通過干涉疊加算法來獲得研究區(qū)的形變速率和形變量。干涉對是利用D-InSAR技術(shù)得到的攜帶研究區(qū)域形變信息的相位差數(shù)據(jù),該數(shù)據(jù)則是通過衛(wèi)星在同一區(qū)域兩次成像所得到的SAR數(shù)據(jù)做相位差分而形成。本研究以重軌干涉測量的成像原理來介紹D-InSAR監(jiān)測地表形變的原理。D-InSAR的幾何原理如圖1所示。
圖1中S1和S2分別代表不同時(shí)刻在不同空間位置的衛(wèi)星傳感器,將S1設(shè)定為主影像成像傳感器,S2設(shè)定為輔影像成像傳感器;R1和R2分別代表在雷達(dá)視線向(LOS)上S1與S2與地表之間的距離;B代表S1與S2之間的距離,稱為空間基線,它在LOS向與S1垂直視線向的投影;α為空間基線B與水平方向之間的夾角,θ為主影像雷達(dá)入射角;h為傳感器與地表之間的距離,P為地表監(jiān)測點(diǎn),P’為P點(diǎn)發(fā)生形變后的位置。在不考慮大氣、噪聲等因素的影響時(shí),S1和S2獲取的SAR影像進(jìn)行差分干涉得到的干涉相位如式(1)。
(1)
式中:Δφ為相位差,(°);λ為雷達(dá)入射波波長,m;B⊥為垂直基線長度,m;h為傳感器與地表之間的距離,m;R1為S1傳感器與地面之間的斜距,m;θ為主影像雷達(dá)入射角,(°);Δd為地表形變量,mm。
式(1)等式右側(cè)前式為地形相位,后式為形變相位,那么去除地形相位后就可得到形變相位φdef,如式(2)。
(2)
獲取N幅SAR影像,生成M幅差分干涉圖(見圖2)。SBAS-InSAR技術(shù)利用這些干涉圖,估計(jì)出平均形變速率之后,通過時(shí)域的高通濾波和空間域的低通濾波去除大氣的影響,使用奇異值分解(SVD)方法估計(jì)形變時(shí)間序列。該方法可以處理多個子集,聯(lián)合同一傳感器獲取的不同時(shí)間段的干涉結(jié)果進(jìn)行求解,提高結(jié)果的精度,降低數(shù)字高程模型(DEM)誤差和大氣相位延遲的干擾。SVD方法在很大程度上解決了由于時(shí)空基線過長而導(dǎo)致的失相干和大氣效應(yīng)問題,且滿足兩個重要要求:①使用所有小基線子集組合的方式增加了時(shí)間采樣率;②得到了更為連續(xù)的形變圖,提高了形變圖空間采樣率[17-18]。
該技術(shù)數(shù)據(jù)處理流程如圖3所示。利用加入相應(yīng)的精密軌道數(shù)據(jù)的若干景SAR數(shù)據(jù)進(jìn)行干涉處理,剔除失相干較為嚴(yán)重的干涉對,使用干涉效果好的干涉結(jié)果進(jìn)行時(shí)序處理,減小大氣效應(yīng)和DEM數(shù)據(jù)帶來的誤差,得到研究區(qū)形變速率信息。
本研究選取的成品油管道區(qū)段位于江西省南部,該區(qū)段共包含318個樁,管道長度16.56 km,起于江西省吉安市萬安縣,止于贛州市的贛縣。研究區(qū)段范圍內(nèi)海拔高度最低38 m,最高590 m,相對高差552 m,總體呈現(xiàn)出東南、西南、西北部地勢偏高,主要地形以中低山、丘陵為主,北部和中部為低丘崗地。管線所經(jīng)區(qū)域,地表主要為第四系砂土、亞砂土等組成的黃色土壤,基巖由白堊系粉砂巖、礫巖等組成。管線經(jīng)過山體一側(cè)形成順向坡,或者順坡上行,形成橫向坡。管道鋪設(shè)和施工道路等形成人工切坡,地表覆蓋層薄,植被恢復(fù)較差,巖土體裸露現(xiàn)象較多,附著能力差,在雨水的沖刷下,水土流失現(xiàn)象明顯,導(dǎo)致易發(fā)生坡體崩塌和滑坡等地質(zhì)災(zāi)害。
共計(jì)使用46景C波段Sentinel-1A升軌SAR數(shù)據(jù),時(shí)間覆蓋范圍為2019年1月8日-2020年7月13日,具體信息見表1。SentineL-1A衛(wèi)星數(shù)據(jù)的成像模式是干涉寬幅模式,在該模式下衛(wèi)星采用TOPS技術(shù),該技術(shù)的一大特點(diǎn)是以80 km寬,20 km長的脈沖帶為成像基礎(chǔ),在方位向相鄰脈沖帶之間存在1.5 km的重疊區(qū)域,在3個子條帶之間存在2 km的重疊區(qū)域,從而保證數(shù)據(jù)對地面覆蓋的連續(xù)性[19-20]。表1所列為Sentinel-1 A影像數(shù)據(jù)主要參數(shù)。
表1 Sentinel-1A影像數(shù)據(jù)主要參數(shù)波段波長/cm軌道方向空間分辨率/m幅寬/km入射角/(°)C5.6升軌5×2025039.386
根據(jù)研究區(qū)特點(diǎn)及數(shù)據(jù)優(yōu)勢,時(shí)間基線閾值設(shè)置為48天,空間基線閾值設(shè)置為300 m,去除相干性較差的干涉對后,參與時(shí)序InSAR處理的干涉對共計(jì)164對。軌道數(shù)據(jù)采用成像21天之后發(fā)布的POD精密軌道數(shù)據(jù),DEM數(shù)據(jù)使用的是30 m的SRTM DEM數(shù)據(jù)。根據(jù)InSAR監(jiān)測結(jié)果進(jìn)行現(xiàn)場實(shí)地調(diào)查,獲取到隱患災(zāi)害點(diǎn)的現(xiàn)場形變信息。
使用架構(gòu)于ENVI遙感圖像處理軟件之上的SARScape軟件進(jìn)行SBAS-InSAR處理,主要包括數(shù)據(jù)讀取、影像裁剪、干涉處理、去平地效應(yīng)、自適應(yīng)濾波、相位解纏,以及第1次估算平均位移速率和DEM校正系數(shù)、第2次估算平均位移速率和DEM校正系數(shù),最后獲得研究區(qū)時(shí)間序列形變和平均形變速率,對形變速率異常值進(jìn)行了剔除,得到了更加準(zhǔn)確的研究區(qū)形變速率圖。
采用SBAS-InSAR技術(shù)對輸油管道重點(diǎn)區(qū)段地表形變進(jìn)行監(jiān)測,得到研究區(qū)平均形變速率,需要說明的是,監(jiān)測到的平均形變速率均為LOS方向的形變值。圖4所示為采用InSAR技術(shù)獲取的輸油管道沿線區(qū)域平均形變速率及管線地表的坡度、坡向數(shù)據(jù)。從圖4(a)可看出,形變速率圖中存在多處空洞,致使空洞區(qū)域無法提取形變信息。從3個方面分析其原因:①由于InSAR技術(shù)本身是以不同時(shí)期相同區(qū)域SAR數(shù)據(jù)之間的相干性為基礎(chǔ)來監(jiān)測地表形變,在其執(zhí)行過程中,會設(shè)置相干性閾值,本研究相干性閾值為0.3,即兩景SAR數(shù)據(jù)同一位置的像素相干性大于0.3時(shí),參與SBAS-InSAR干涉疊加分析,小于0.3的像素則會被剔除,而被剔除的區(qū)域就會在結(jié)果中以空值形式顯示;②當(dāng)SAR入射角小于地表坡度,SAR脈沖照射到位于傳感器掃描方向一致的斜坡時(shí),脈沖到達(dá)斜坡頂部和底部的斜距差小于實(shí)際地面之間的距離,致使SAR影像上的斜面長度被縮短了,出現(xiàn)了“透視收縮”[21],通過疊加分析形變速率和坡度數(shù)據(jù)可以知道,該區(qū)域最大坡度為37.65°,而本研究使用的SAR衛(wèi)星影像入射角為39.386°,因此一定不是“透視收縮”引起的形變速率空洞;③使用的Sentinel-1A衛(wèi)星運(yùn)行軌道方向?yàn)樯?,即衛(wèi)星繞地球由南偏東向北偏西飛行,傳感器掃描方向與飛行方向一致,那么Sentinel-1A對與其相垂直的山體的敏感程度較弱,且在坡度較大時(shí)易在SAR影像形成“陰影”。分析存在空洞的區(qū)域,由于坡度偏小,南北和東西方向山體均勻分布,不足以將形變空洞歸結(jié)于此。
通過上述分析可知,形變速率存在空洞的主要原因是由于時(shí)間序列SAR影像的不斷變化造成影像失相干,在數(shù)據(jù)運(yùn)算過程中失相干區(qū)域被剔除,從而無形變結(jié)果。
根據(jù)試驗(yàn)區(qū)形變速率結(jié)果對輸油管道全線形變情況以及影響范圍進(jìn)行討論分析。受InSAR技術(shù)自身的限制,形變監(jiān)測結(jié)果存在一定誤差,在實(shí)際分析形變結(jié)果時(shí)會將這部分誤差剔除掉,即形變速率在誤差絕對值范圍內(nèi)浮動,將其視為未發(fā)生地表形變。利用全線路空間域形變速率直方圖統(tǒng)計(jì)結(jié)果(見圖5),設(shè)置誤差允許值(ε=-7.3 mm),將該值作為判斷是否發(fā)生地表形變的固定閾值。為了準(zhǔn)確地判斷地表是否發(fā)生形變,將形變速率在一定長度范圍內(nèi)的一階導(dǎo)數(shù)不等于零作為限定條件,限定條件為:v≥-7.3 mm,v/(p·20 m)≠0,其中v為形變速率,p為像素個數(shù)。當(dāng)形變速率小于-7.3 mm,且形變速率在一定長度范圍內(nèi)的一階導(dǎo)數(shù)不等于零時(shí),確定該區(qū)域發(fā)生形變。
利用上述限定條件提取全線地表形變情況,共獲取到4處形變區(qū)域,標(biāo)記為A1、A2、A3、A4。其中,A2形變區(qū)域發(fā)生在斜坡上(見圖6),其他均發(fā)生在平原區(qū)域,因此重點(diǎn)對A2區(qū)域進(jìn)行形變分析。為了能更加清晰地反應(yīng)形變對管線的影像程度和范圍,根據(jù)時(shí)序InSAR技術(shù)獲取的平均形變速率圖,繪制了該不穩(wěn)定斜坡東西方向和南北方向的平均形變速率剖面曲線。東西方向沿斜坡蠕動方向進(jìn)行繪制(見圖6(a)中的PP’線段),南北方向沿管道上行方向進(jìn)行繪制(見圖6(a)中的NN’線段),東西向和南北向平均形變速率剖面曲線如圖6(b)、圖6(c)所示。該形變區(qū)主要覆蓋管道東側(cè)斜坡,該斜坡最大高程298 m、最小174 m、相對高差124 m、坡面長度376 m、坡度約32°。
由形變速率圖和剖面曲線圖可知,該形變區(qū)涉及范圍0.02 km2,東西向?qū)挾?20 m,南北向長度144 m,也就是說輸油管道在該位置144 m內(nèi)存在安全隱患,從東西向剖面圖中可以看出(見圖6(b)),輸油管道兩側(cè)存在兩處沉降中心,東側(cè)斜坡平均形變速率最大達(dá)到-43 mm/a,西側(cè)平坦地表最大形變速率達(dá)到-35 mm/a,輸油管道穿越位置約-30 mm/a。依據(jù)Colesanti C教授的經(jīng)驗(yàn)總結(jié)[22],對于C波段SAR數(shù)據(jù),當(dāng)斜坡的LOS向形變速率絕對值大于2 mm/a時(shí),可認(rèn)為斜坡處于不穩(wěn)定狀態(tài)。顯然,該斜坡已處于不穩(wěn)定狀態(tài),存在滑坡的可能,該監(jiān)測結(jié)果與現(xiàn)場監(jiān)測成果進(jìn)行核驗(yàn),發(fā)現(xiàn)形變區(qū)位置及范圍與實(shí)地調(diào)查結(jié)果基本一致。
為了研究在輸油管道附近監(jiān)測到的4個形變區(qū)域在時(shí)間序列上的形變過程,利用SBAS-InSAR技術(shù)獲取的時(shí)間序列累積形變數(shù)據(jù)繪制了4個形變區(qū)形變中心點(diǎn)在2019年1月-2020年7月的累積形變曲線圖(見圖7)??紤]到斜坡失穩(wěn)及水土流失等地質(zhì)災(zāi)害受雨水影響嚴(yán)重[23-26],故收集并統(tǒng)計(jì)了監(jiān)測時(shí)間段內(nèi)研究區(qū)降雨情況。結(jié)合降雨情況,對監(jiān)測到的4個區(qū)域形變過程及成因進(jìn)行分析說明。
從圖7中可以看出,4個區(qū)域從2019年2月開始,均發(fā)生緩慢形變。其中,A1、A2、A3在2019年3月開始發(fā)生形變,A4在2019年2月開始發(fā)生沉降,至2019年8月趨于平緩,累積形變值分別達(dá)到-36 mm、-44 mm、-59 mm、-53 mm,在之后近4個月內(nèi)形變較為穩(wěn)定。根據(jù)降雨情況也可以看出,研究區(qū)在2019年1月-2019年7月,降雨天數(shù)均在10天以上,降雨量也隨之增大,極有可能是增加斜坡和易積水區(qū)域的形變量和形變速率的主要因素。2020年1月,A2斜坡、A3和A4平原區(qū)沉降速率加劇,分別在2020年3月中旬和2020年5月初達(dá)到沉降谷值,其中A2斜坡的形變值達(dá)到-61 mm,相比去年8月增加25 mm。根據(jù)降雨情況可以推斷,由于從2020年1月開始到2020年6月結(jié)束,研究區(qū)降雨天數(shù)增加后減少,期間存在連續(xù)降雨情況,斜坡以及斜坡下水溝經(jīng)過長時(shí)間雨水浸泡,坡面和溝壁土質(zhì)軟化,從而形變增加;A3和A4平原區(qū)形變值分別達(dá)到-79 mm、-68 mm,該區(qū)域?yàn)槠皆?,因此可判斷地質(zhì)災(zāi)害類型主要是沉降和水土流失。2020年2月,A1緩坡區(qū)形變速率增加,至2020年5月達(dá)到谷值-74 mm,之后形變區(qū)域平緩,與降雨天數(shù)成正相關(guān)。監(jiān)測時(shí)間段內(nèi)形變速率的變化與月降雨天數(shù)之間的相關(guān)性達(dá)到0.69。由此可以推斷,地表水土流失、斜坡失穩(wěn)的主要形成原因是降雨,降雨對地表和斜坡的沖擊將不斷使之發(fā)生不均勻形變。因此,在每年雨季應(yīng)及時(shí)對輸油管道周邊易積水、斜坡不穩(wěn)定的區(qū)域加強(qiáng)防范,避免發(fā)生重大地質(zhì)災(zāi)害。
將SBAS-InSAR技術(shù)應(yīng)用于輸油管道重點(diǎn)區(qū)段的地質(zhì)災(zāi)害監(jiān)測,使用46景覆蓋輸油管道的Sentinel-1A數(shù)據(jù)獲取該區(qū)域時(shí)間序列形變信息。結(jié)合坡度、坡向數(shù)據(jù)對形變速率結(jié)果存在空洞進(jìn)行分析,結(jié)果表明:
(1) 不同時(shí)間段的SAR影像存在失相干是造成最終空洞的主要原因。
(2) 當(dāng)山體坡度大于SAR傳感器入射角且坡向與衛(wèi)星運(yùn)行軌道方向一致時(shí),易出現(xiàn)透視收縮和陰影,而使最終的形變速率結(jié)果出現(xiàn)空洞。
因此,在后續(xù)的應(yīng)用當(dāng)中,應(yīng)結(jié)合穿透能力更強(qiáng)的SAR數(shù)據(jù)進(jìn)行輸油管道地表形變監(jiān)測。通過分析由SBAS技術(shù)獲取的形變速率值,獲取全線路發(fā)生形變區(qū)域有4處,并對其中一處發(fā)生在斜坡的形變區(qū)域進(jìn)行了形變速率和影響范圍的分析,發(fā)現(xiàn)該斜坡最大年平均形變速率達(dá)到-43 mm/a,對管道及周邊0.02 km2造成影響。最后引入降雨天數(shù)數(shù)據(jù),分析這4處形變區(qū)域與降雨之間的關(guān)系,結(jié)果表明降雨增多是致使地表發(fā)生形變的原因之一。
InSAR技術(shù)作為地質(zhì)災(zāi)害監(jiān)測的新興技術(shù)手段,能夠在廣域、短周期內(nèi)對地表進(jìn)行形變監(jiān)測,及時(shí)發(fā)現(xiàn)地表形變異常,能夠?yàn)檩斢凸艿赖刭|(zhì)災(zāi)害防治提供重要數(shù)據(jù)支撐。不可忽略的是,InSAR技術(shù)是以同一區(qū)域不同時(shí)間點(diǎn)的SAR之間的相干性作為算法實(shí)現(xiàn)的基本條件。因此,在地表相干性較差時(shí)會出現(xiàn)形變速率空值,無法獲取形變的情況,這是InSAR存在的不足,也是目前眾多學(xué)者不斷探討的問題。后續(xù)的實(shí)驗(yàn)過程中,將嘗試降低相干性閾值,降低空洞概率,或使用波長更長的SAR數(shù)據(jù)對失相干區(qū)域進(jìn)行二次監(jiān)測。