利用Stacking/SBAS技術在滇西北地區滑坡隱患的識別對比

董繼紅,張 肅,梁京濤,楊 磊,趙 聰

(1.四川省綜合地質調查研究所 稀有稀土戰略資源評價與利用四川省重點實驗室 四川 成都 610081;2.四川省智慧地質大數據有限公司,四川 成都610081)

滑坡是全球最常見、分布最廣泛的地質災害類型之一,不僅造成人員傷亡和基礎設施破壞,還會形成鏈式災害,產生二次破壞[1]。據統計,1995—2014年,全球共發生3 876 處災難性滑坡,造成16萬余人死亡和1.1萬余人受傷[2]。準確識別潛在滑坡并繪制滑坡隱患分布圖是防災減災工作的重點內容[3]。預先獲知滑坡隱患的分布位置、發育特征和變形趨勢,對預測潛在風險區和防災治理尤為重要。目前滑坡隱患識別手段多基于光學影像和現場調查等,工作強度大且效率低,亟需尋找一種高效、高精度技術手段對滑坡隱患進行識別和監測。

合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)作為一種非接觸獲取地表形變的技術[4],具有全天時、全天候、穿云透霧等優點,在滑坡隱患識別中被廣泛應用[5]。隨著InSAR技術的不斷發展,形變監測從起初精度較低且提取信息單一的D-InSAR技術,發展到后來的時間序列InSAR(time series interferometric synthetic aperture radar,TS-InSAR)技術[6]。TS-InSAR技術能利用較多SAR影像,最大限度克服時空去相干、大氣延遲和數字高程模型(digital elevation model,DEM)誤差等的影響,拓展了InSAR技術的應用領域[7-8]。特別是Sandwell 等[9]提出的Stacking技術,能通過較少SAR數據生成的差分干涉圖進行相位堆疊,快速獲取地表形變結果而被廣泛應用于滑坡隱患的早期識別[10]。Berardino等[11]提出的小基線集(small baselines subset,SBAS)技術通過將較短時空基線干涉對組合,在提高監測點密度的同時獲取監測點形變量。目前上述兩種技術均被廣泛應用于滑坡隱患識別[12],如Liang等[13]以川西片區為例,開展大范圍Stacking技術與SBAS技術滑坡隱患識別對比分析,何佳陽等[14]利用InSAR和SAR技術開展雅礱江沿線西昌區域滑坡隱患識別。上述文獻對比了不同InSAR技術在滑坡隱患識別中的應用效果,但未分析隱患識別的差異原因,而且在西南高山峽谷區利用不同InSAR技術進行隱患識別和對比分析尤為重要。

本研究利用Stacking和SBAS技術處理Sentinel-1數據以獲取地表形變速率,并結合SAR數據的可視性進行隱患識別。將不同技術的隱患識別效果進行對比,詳細分析兩種InSAR技術在滑坡隱患識別中的差異性和影響因素。

1 研究區域

研究區位于云南省麗江市與怒江、大理兩州的交界處(圖1),屬云貴高原與橫斷山脈結合部位,地勢西北高,東南低。金沙江、瀾滄江和怒江流經該區域,為泛三江源區域,地貌復雜多樣,屬于低緯高原季風氣候。受斷裂、構造的影響,該區滑坡頻繁發生[15],使得當地居民的生命財產和安全長期受到威脅。因此,對該區滑坡隱患識別進行系統分析,不僅對該區域具有重要意義,也能為地質環境相似的西南山區滑坡隱患識別、治理提供參考。

2 數據和采用的技術方法

2.1 數據

收集研究區2018年11月—2021年12月期間92期Sentinel-1衛星SAR影像數據,如圖1所示,由于Sentinel-1數據幅寬較大,需根據研究區大小對SAR影像進行裁剪,圖1為數據覆蓋范圍。表1為所采用SAR數據集的基本參數。收集Sentinel-1衛星對應的POD(precise orbit ephemerides)精密軌道星歷數據,用于去除因軌道誤差引起的系統誤差。引入AW3D30 DSM數據去除InSAR干涉處理中的地形相位并輔助SAR影像地理編碼,該數據也被用來計算研究區的可視性分布情況。

2.2 技術方法

本研究基于Stacking和SBAS技術獲取地表形變結果,因SAR數據為側視雷達成像,在地形起伏區域存在幾何畸變現象,故首先采用改進后的R指數進行可視性計算,分析地形特征對結果的影響。

2.2.1R指數

SAR衛星采用側視雷達成像,因此在地形起伏較大的區域容易形成幾何畸變,不僅降低了SAR影像的質量,還會對隱患識別造成漏判、誤判,因此有必要準確識別幾何畸變區域,以提高InSAR識別隱患的準確率[15]。常見的幾何畸變現象包括疊掩、陰影和透視收縮。

SAR影像的幾何畸變與雷達衛星的入射角、方位角、地面坡度、坡向有關,本研究采用由Notti等[16]提出的R指數模型,以及Ren等[20]為識別較遠被動疊掩區域而提出的改進R指數公式計算SAR影像幾何畸變區域。

R=sin{θ+arctan[tan(α)×cos(φ-β)]}×Sh×La×Fa。

(1)

其中:θ為衛星的入射角;φ為衛星的方位角;α為地面坡度;β為坡向;Sh為山體陰影系數,陰影區域取0,其他區域取1.0;La為頂底倒置(Layover)系數,主動頂底倒置和被動頂底倒置區域的值為0,其他區域的值為1.0;Fa為較遠被動疊掩系數,主動疊掩區域和被動疊掩區域的值為0,其他區域的值為1.0。Sh、La、Fa三個系數由ArcGIS計算獲得,通過計算可視性分析,可直觀了解不同地形特征對結果的影響。

2.2.2 Stacking技術

Stacking技術通過對多幅差分干涉圖加權平均獲取地表形變速率結果。該技術的前提是假設地表形變趨勢為線性形變,通過相位堆疊有效抑制大氣相位和DEM誤差,從而提高形變信息的精度[17]。該技術通過處理少量SAR數據快速獲取大范圍滑坡隱患分布圖,結果可靠,被廣泛應用于滑坡隱患的早期識別[18]。

2.2.3 SBAS技術

SBAS技術基于配準后的多幅SAR影像,根據設置的時間基線閾值和空間基線閾值生成小基線集,進行差分干涉處理,經濾波降低相位噪聲增加高相干點位,再對濾波后的干涉圖進行相位解纏,利用奇異值分解(singular value decomposition,SVD)小基線集得到形變速率,最后對形變速率積分得到監測時間段內的累積形變量[19]。

3 數據處理過程

利用GAMMA軟件對所獲取的92景Sentinel-1A數據進行InSAR數據處理,具體處理流程如圖2所示。主要包括以下5個步驟:

1) 數據預處理:選擇位于中間的一期SAR影像作為主影像,根據研究區范圍對其進行裁剪,完成主影像與DEM的配準、裁剪,得到SAR坐標系下的DEM數據;同時完成主影像與其他SAR數據配準;限制時間基線不超過60 d,垂直基線不大于±250 m,生成時空基線,共獲取380對組合干涉對。

2) 干涉工作流處理:對上述組合干涉對進行主輔影像共軛相乘得到差分干涉圖;采用精密軌道數據去除平地相位;采用DEM數據模擬并剔除地形相位;采用自適應濾波的方法對差分干涉圖進行濾波處理,以消除或減弱噪聲,生成相干系數圖;根據R指數計算陰影疊掩區域,同時因水域沒有有效干涉信息,對差分干涉圖進行掩膜處理,去除陰影疊掩區域和水域。采用最小費用流(minimum costflow,MCF)法進行相位解纏,使用相干性掩模避開相干性較低、相位不可靠的區域,設定0.3為相干值的閾值,低于此相干值的區域設為空值并不參與解纏計算。

3) Stacking計算對相位解纏結果進行相位堆疊計算,獲取形變速率結果,然后將形變轉換到LOS向,借助DEM數據進行地理編碼獲取地理坐標系下形變結果。

4) 時間/空間域變形估算

該步驟是在上述第2)步完成之后進行,①相鄰點間參數估計:將點目標連接構成不規則三角網,依據點間連接關系求解相鄰點差分相位差;②殘余高程計算和線性變形:依據基線組合,估算相鄰點間的線性變形速率和DEM誤差;③殘余相位低通濾波:從差分干涉相位中減去步驟①中差分相位得到殘余相位,對殘余相位進行空間域低通濾波得到濾波后的殘余相位;④奇異值分解處理:根據短基線像對組合關系,對上一步得到的濾波后殘余相位進行奇異值分解(Singular Value Decomposition,SVD)處理,求解每個影像對應時刻的大氣相位和非線性變形相位;⑤大氣相位和非線性變形相位計算:對奇異值分解得到的大氣相位和非線性變形相位進行空間域高通濾波,得到大氣相位,并對濾波后的相位序列進行時域低通濾波,得到非線性變形相位。

5) 形變量計算

將上一步獲取的非線性形變結果和線性形變結果相加,根據時間基線參數獲得形變量結果,將結果轉換到視線向形變,利用DEM數據進行基準修正和地理編碼,獲取地理坐標系下形變速率和時間序列結果。

4 結果分析

4.1 可視性分析

假設SAR衛星為太陽所在位置,利用表1中Sentinel-1數據參數,通過ArcGis軟件對衛星高度角和方位角計算Sh、Fa和La系數值,然后通過式(1)計算得到最終的SAR地形可視性分布圖(圖3)。當R>sinθ,為好可視性區域;當0

(a)地形可視性分布;(b)幾何畸變面積統計

4.2 地表形變結果分析

利用Stacking和SBAS技術分別對2018年11月—2021年12月的Sentinel-1數據進行處理,得到研究區雷達視線向形變速率(圖4),正值表示地表形變靠近衛星飛行方向;負值表示地表形變背向衛星飛行方向。分別利用Stacking和SBAS技術獲取的LOS向形變速率為-83～51和-79～29 mm/yr。對形變參考區域的平均值和標準差進行統計分析,確定以-10～10 mm/yr作為相對穩區,該閾值外的速率對應不同程度的形變,值越大表示形變越強烈,其中形變最大區域主要集中在河谷兩側。通過對比監測結果發現,SBAS技術獲取的結果點位較稀疏,究其原因一方面受到幾何畸變的影響,另一方面由某些區域相干性不連續導致。

圖4 沿LOS向平均形變速率圖

4.3 相關性分析

對兩種不同InSAR技術獲取的Stacking和SBAS結果進行相關性分析(圖5),發現兩者線性相關性系數達0.72,表明兩種方法得到的結果具有較高的一致性。通過相關性分析,一方面驗證兩種技術獲取結果的一致性,另一方面可驗證內符合精度。

圖5 Stacking與SBAS技術年平均形變速率相關性圖

4.4 地質災害隱患識別對比分析

為了對比Stacking和SBAS兩種技術對地質災害隱患的識別能力,借助該區域的DEM數據和光學影像數據,分別對結果進行解譯和統計。Stacking和SBAS技術分別識別出隱患點45處和36處,結合野外驗證得到滑坡隱患分布圖(圖6),圖中綠色圓點為Stacking和SBAS技術共同識別的滑坡隱患,紅色圓點為僅被Stacking技術識別出來的隱患點。表2為不同InSAR技術識別的滑坡隱患數目,可以看出,Stacking技術識別的滑坡隱患數目較多,但SBAS技術識別隱患的準確率更高。

圖6 不同技術識別滑坡隱患分布圖

為分析兩種技術識別滑坡隱患差異的原因,對識別的32處隱患點的最大LOS向形變速率和R指數進行統計(表3)。由表3可以看出,能被兩種技術共同識別的隱患點中,SBAS和Stacking技術能夠識別的最小形變速率分別為14.8、15.6 mm/yr,且R指數均大于0.8,屬于好可視性區域;僅被Stacking技術識別出的隱患點最小形變速率為14.2 mm/yr,R指數為0.52～0.63,屬于透視收縮區域。

表3 不同InSAR技術識別滑坡特征統計表

為進一步分析兩種技術識別的差異性,選擇兩處典型滑坡點進行分析:第一處滑坡點是被兩種技術共同識別的中村滑坡;第二處滑坡點為僅被Stacking技術識別的營盤鎮滑坡。

中村滑坡位于蘭坪縣,坡向朝東。從光學影像(圖7(a))可以看出植被覆蓋較差,滑坡前緣以瀾滄江為界,由Stacking和SBAS結果可以看出滑坡體存在兩處明顯變形區域,即滑坡的中前部右側和中后部左側。對比兩種技術可知,Stacking結果的有效監測點更密集,SBAS結果形變速率更大。為進一步分析不同InSAR技術的監測結果和R指數值,沿圖7(b)中的剖面線提取形變曲線,如圖8(a)所示,其中圖7(c)中紅色箭頭標注區域對應圖8(a)中灰色陰影區域,兩種技術識別的變形區域一致,形變速率均超過-15 mm/yr,R指數均大于0.8。因此,能被兩種技術共同識別的條件是形變速率大于-10 mm/yr且斜坡處于好可視性區域。另外,造成SBAS結果點位稀疏的主要原因是部分區域的相干性在時間域不連續。而SBAS技術獲取的形變速率更大是因為Stacking技術是對形變相位的加權平均,適用于線性形變,因此監測過程會丟失非線性形變信息。

營盤鎮滑坡位于蘭坪縣,坡向朝西,光學影像如圖9(a)所示。SBAS結果無明顯形變信息且監測點位稀疏,在R指數圖中絕大部分區域為透視收縮區域。為深入分析不同InSAR技術和R指數對結果的影響,沿圖9(b)中黑色剖線提取形變曲線,圖9(c)中紅色箭頭區域對應圖8(b)灰色陰影區域,兩種技術在變形區域均出現明顯的變形趨勢,但SBAS結果中僅有4個像素值超過10 mm/yr,且最大值不超過12 mm/yr;在Stacking結果中陰影區域速率均超過10 mm/yr,部分區域變形速率超過12 mm/yr,R指數小于0.6,屬于透視收縮區域。因此,該區域未被SBAS技術識別的原因是坡向朝西,在升軌數據中處于透視收縮區域,造成在地理編碼之后監測點位稀疏,形變信息喪失;在Stacking結果中有明顯形變信息的原因是兩者求解形變速率的原理不同。但根據SAR數據成像幾何關系可知該隱患點在降軌數據中應處于好可視性區域,可以彌補SBAS結果在該區域識別的劣勢,增加SAR數據的可觀測區域。

為證明降軌數據可以提高升軌數據的好可視性范圍,對營盤鎮滑坡的降軌Sentinel-1數據進行處理,并提取該滑坡體不同位置的3個時序監測點進行變形特點分析(見圖10中的P1、P2、P3),對比圖9(d)可以發現,降軌SBAS結果中變形特征明顯,最大形變速率為-38 mm/yr,表明降軌數據可以提升升軌數據在透視收縮區域的監測效果。圖11為營盤鎮滑坡時間序列曲線圖,發現該滑坡中前部(點P2、P3)變形強烈,且變形趨勢接近,最大形變量為-83.3 mm,目前仍處于持續變形狀態,存在較大隱患。

圖10 營盤鎮形變速率圖

圖11 營盤鎮滑坡時序曲線圖

5 結論

以云南大理為例,利用Stacking和SBAS技術對Sentinel-1數據進行處理,獲取研究區地表形變速率,并結合SAR數據可視性開展滑坡隱患識別差異性及影響因素分析,結果表明:

1) 兩種InSAR技術獲取的結果線性相關性系數為0.72,具有較高的一致性。

2) 透視收縮區域占研究區總面積的50%以上,在透視收縮區域內,Stacking技術識別效果優于SBAS技術,因此采用Stacking技術識別的滑坡隱患數量多于SBAS技術,結果表明Stacking技術可以提高隱患識別能力。

3) 結合野外驗證結果顯示Stacking技術識別準確率為71.1%,SBAS技術識別準確率較高,為76.5%,表明SBAS技術識別結果可靠性更強。

4) 研究表明在西南山區開展滑坡隱患識別,應結合升降軌數據,以減少透視收縮區和疊掩陰影區面積占比,數據處理方法首選Stacking技術。此外,為提高滑坡隱患識別準確率并獲取時序變形結果應結合使用SBAS技術。