基于升降軌Sentinel-1A SBAS InSAR的西南山區(qū)滑坡隱患早期識別研究

劉排英 賀少帥 王鵬生 孫 亮

1 石家莊鐵路職業(yè)技術(shù)學(xué)院測繪工程系，石家莊市四水廠路18號，050041 2 中國航天空氣動(dòng)力技術(shù)研究院，北京市云崗西路17號，100074

滑坡具有隱蔽性強(qiáng)、破壞程度大、發(fā)生頻率高和分布廣泛等特點(diǎn)。貴州省位于中國西南地區(qū)，地貌類型主要為高原、山地、丘陵和盆地[1]，氣候溫暖濕潤、多云多雨、植被茂密，這會(huì)增大滑坡的隱蔽性和多發(fā)性，也會(huì)增加光學(xué)遙感成像和解譯滑坡的難度。《2020年度貴州省地質(zhì)災(zāi)害防治方案》提出“空天地”一體化的地災(zāi)隱患調(diào)查，使用InSAR開展廣域地表形變監(jiān)測，識別、提取出更多地質(zhì)災(zāi)害隱患點(diǎn)。

衛(wèi)星InSAR利用微波頻段的相位信息探測雷達(dá)視線向cm級甚至mm級微小地表形變，具有廣覆蓋、全天候、全天時(shí)、高重訪和主動(dòng)遙感的優(yōu)勢。由于傳統(tǒng)D-InSAR受時(shí)空去相干和大氣擾動(dòng)等限制，以PS InSAR[2]和SBAS InSAR[3]為代表的時(shí)序InSAR應(yīng)運(yùn)而生，可以在很大程度上解決D-InSAR的應(yīng)用限制。在山區(qū)環(huán)境中，分布式散射體(distributed scatterers，DS)數(shù)量遠(yuǎn)大于永久散射體數(shù)量，因此SBAS InSAR比PS InSAR更適用于山區(qū)場景的滑坡形變監(jiān)測，能夠顯著減小SAR影像的數(shù)量需求。

2014年以來，中高分辨率Sentinel-1A數(shù)據(jù)免費(fèi)向用戶開放，極大推進(jìn)了InSAR在我國的普及和應(yīng)用。針對重大地質(zhì)災(zāi)害隱患的早期識別，我國學(xué)者相繼提出空-天-地“三查”技術(shù)體系和形態(tài)-形變-形式“三形觀測”的技術(shù)思路，強(qiáng)調(diào)InSAR在地質(zhì)災(zāi)害監(jiān)測中的普查作用，使得InSAR在地質(zhì)災(zāi)害調(diào)查中逐漸得到認(rèn)可[4-5]。本文以貴州省最近一次滑塌體量大且致災(zāi)嚴(yán)重的雞場鎮(zhèn)滑坡為研究背景，獲取2018-07～2019-07覆蓋貴州省水城縣雞場鎮(zhèn)滑坡發(fā)生前的31期升軌和30期降軌Sentinel-1A數(shù)據(jù)，結(jié)合30 m分辨率的SRTM DEM數(shù)據(jù)，使用SBAS InSAR開展滑坡發(fā)生前的地表形變監(jiān)測和隱患早期識別研究。該研究可分析和評估SBAS InSAR在中國西南山區(qū)地表形變監(jiān)測中的應(yīng)用效果，為貴州省以及中國西南山區(qū)滑坡隱患調(diào)查與早期識別提供技術(shù)參考。

1 研究區(qū)域及數(shù)據(jù)

貴州省位于中國西南地區(qū)，地貌類型以高原、山地、丘陵和盆地為主，其中山地和丘陵占92.5%[1]。貴州省年平均降雨量超過1 000 mm，主要集中在夏季。凹凸不平的地貌類型和時(shí)域集中的強(qiáng)降雨特征易引發(fā)滑坡、泥石流等地質(zhì)災(zāi)害。2019-07-23水城縣雞場鎮(zhèn)發(fā)生大型滑坡(圖1白色輪廓區(qū))，該滑坡體屬于構(gòu)造侵蝕、剝蝕型中山地貌，滑坡在高海拔地帶發(fā)生，演變?yōu)檠貎蓷l溝道分布的長滑移距離型泥石流，整個(gè)滑坡區(qū)域長約1 300 m，前后緣高差約460 m，體積約為191.2×104m3，屬于典型的高位遠(yuǎn)程滑坡，滑坡體和碎屑流的主滑方向約為26°[6]。

圖1 研究區(qū)及Sentinel-1A覆蓋范圍Fig.1 Study area and Sentinel-1A coverage

研究區(qū)范圍為104°35′～104°49′E、26°12′～26°27′N(圖1紅框)。實(shí)驗(yàn)采用歐空局Sentinel-1A漸進(jìn)式地形觀測模式干涉寬幅成像模式VV極化SAR數(shù)據(jù)集，右視側(cè)視角約為39°。覆蓋雞場鎮(zhèn)滑坡的Sentinel-1A數(shù)據(jù)有升軌第128軌、降軌第164軌，升軌成像時(shí)段為2018-07-25～2019-07-20，共31期；降軌成像時(shí)段為2018-07-27～2019-07-22，共30期，均為災(zāi)前影像。SAR數(shù)據(jù)集空間覆蓋范圍如圖1藍(lán)框所示。輔助數(shù)據(jù)有POD精密定軌星歷數(shù)據(jù)和分辨率為30 m的SRTM DEM數(shù)據(jù)，其中POD數(shù)據(jù)用于精化軌道坐標(biāo)，SRTM DEM數(shù)據(jù)用于SAR數(shù)據(jù)精配準(zhǔn)、去地形相位和地理編碼等。

2 研究方法

本研究使用FaSHPS算法[7]統(tǒng)計(jì)、提取同質(zhì)像元作為DS候選點(diǎn)，設(shè)定干涉對集的平均相干閾值篩選同質(zhì)像元作為DS點(diǎn)。由于研究區(qū)植被較多，去相干情況較為嚴(yán)重，因此使用像對干涉相干性最優(yōu)且連通的原則生成像對干涉連接圖[8]。使用StaMPS軟件進(jìn)行SBAS InSAR處理，反演地表形變場。整體技術(shù)流程如圖2所示。

圖2 總體技術(shù)流程Fig.2 Overall technical flow chart

貴州省水城縣地形起伏較大、植被茂密，跨季影像相互配準(zhǔn)難度較大。本文使用POD精密定軌星歷數(shù)據(jù)精化軌道坐標(biāo)，計(jì)算輔影像相對于主影像的初始偏移量，使精度達(dá)到像元級。在此基礎(chǔ)上，使用強(qiáng)度互相關(guān)算法從主影像、輔影像中分區(qū)塊提取若干同名像點(diǎn)，對每個(gè)區(qū)塊同名像點(diǎn)方位向和距離向的偏移量進(jìn)行最小二乘擬合，求得配準(zhǔn)多項(xiàng)式系數(shù)，配準(zhǔn)精度達(dá)到0.1像元。配準(zhǔn)多項(xiàng)式系數(shù)確定后，使用三次卷積內(nèi)插法將輔影像重采樣到主影像空間上。

假設(shè)共有成像日期為[t0,t1,…，tN]的N+1景SAR影像，并且所有SAR影像均已精配準(zhǔn)。通過設(shè)定時(shí)間、空間基線閾值或使用干涉相干性最優(yōu)且連通的原則，得到M個(gè)干涉像對，則干涉像對數(shù)M與影像數(shù)N+1滿足關(guān)系式：

(1)

影像tA和影像tB生成的第i幅干涉圖(tA

為消除Sentinel-1A數(shù)據(jù)在山地、丘陵等地區(qū)的陰影、疊掩現(xiàn)象，本文對研究區(qū)成像時(shí)段相對一致的升軌和降軌數(shù)據(jù)分別進(jìn)行獨(dú)立處理，獲取視線向形變速率(正值表示地表向靠近衛(wèi)星方向發(fā)生形變，負(fù)值表示向遠(yuǎn)離衛(wèi)星方向發(fā)生形變)。根據(jù)InSAR成像幾何可知，視線向形變量在水平方向只能分解出東西方向分量，無南北方向分量，這也是衛(wèi)星InSAR成像幾何的固有缺陷。水平形變量VE-W和垂直形變量VV表達(dá)式為：

VE-W=VLOS·sinα，VV=VLOS·cosα

(2)

式中，VLOS為雷達(dá)視線向形變速率，α為降軌和升軌入射角。

為了更好地分析坡體形變特征，將VLOS投影至斜坡方向：

(3)

式中，γ為雷達(dá)衛(wèi)星視向線方向和斜坡坡度方向間夾角，詳細(xì)公式見文獻(xiàn)[9]。

3 數(shù)據(jù)處理與結(jié)果分析

3.1 SBAS InSAR處理

使用POD精密定軌星歷數(shù)據(jù)精化Sentinel-1A軌道坐標(biāo)，生成單視復(fù)數(shù)(single look complex, SLC)影像，并按照研究區(qū)范圍裁剪SLC數(shù)據(jù)。使用干涉相干性最優(yōu)且像對連通的原則分別對31期(ID:0～30)升軌影像和30期(ID:0～29)降軌影像進(jìn)行干涉像對連接，分別生成64個(gè)干涉對和68個(gè)干涉對(圖3)。

圖3 干涉像對Fig.3 Interference image pairs

對研究區(qū)SLC數(shù)據(jù)進(jìn)行SBAS InSAR處理：1)選擇升軌、降軌第1期(ID:0)影像作為超級主影像，多視視數(shù)設(shè)為方位向×距離向=1×4，使用SRTM DEM進(jìn)行主輔影像精配準(zhǔn)、差分干涉處理，計(jì)算相干圖、幅度圖，分辨率分別為13.9 m×14.5 m(升軌)和14.0 m×14.5 m(降軌)；2)使用Goldstein算法對差分干涉圖進(jìn)行濾波處理；3)設(shè)定相干系數(shù)閾值為0.2進(jìn)行掩膜，使用MCF算法對差分干涉圖進(jìn)行相位解纏；4)設(shè)定平均相干閾值為0.41，升軌共提取26 038個(gè)DS點(diǎn)，降軌共提取50 049個(gè)DS點(diǎn)；5)在DS點(diǎn)上建立和解算觀測方程，使用時(shí)空濾波器去除大氣相位和高程殘差相位，使用SVD算法解算DS點(diǎn)的時(shí)間序列形變；6)使用SRTM DEM對反演的地表形變和幅度圖進(jìn)行地理編碼。

3.2 結(jié)果分析與討論

升軌、降軌Sentinel-1A影像SBAS InSAR反演的2018-07～2019-07研究區(qū)地表雷達(dá)視線向平均形變速率結(jié)果如圖4所示。DS點(diǎn)空間分布存在3種情況：1)西向陽坡DS點(diǎn)多數(shù)存在于升軌Sentinel-1A影像上，少量存在于降軌影像上；2)東向陽坡DS點(diǎn)多數(shù)存在于降軌Sentinel-1A影像上，少量存在于升軌影像上；3)地勢較為平坦的DS點(diǎn)同時(shí)存在于升軌和降軌Sentinel-1A影像上。上述現(xiàn)象是由衛(wèi)星雷達(dá)右視側(cè)視成像造成，可通過融合升降軌時(shí)序InSAR形變場提升地勢起伏較大山區(qū)、丘陵的形變監(jiān)測能力。

雞場鎮(zhèn)滑坡體(圖4紫色區(qū)域)在發(fā)生滑移前植被覆蓋度高，DS點(diǎn)較稀疏，上沿及周邊DS點(diǎn)形變速率約為0，說明雞場鎮(zhèn)滑坡發(fā)生前滑坡體上沿及周邊區(qū)域并未出現(xiàn)大型地表形變，表現(xiàn)為穩(wěn)定狀態(tài)。查詢雞場鎮(zhèn)滑坡發(fā)生前兩周天氣狀況可知(表1)，雞場鎮(zhèn)滑塌體在滑坡發(fā)生當(dāng)天有強(qiáng)降雨，在滑坡發(fā)生前兩周內(nèi)有11 d出現(xiàn)不同程度降雨。因此初步判斷雞場鎮(zhèn)滑坡是連續(xù)降雨、強(qiáng)降雨誘發(fā)的突發(fā)性滑坡，已超出Sentinel-1A的12 d重訪周期監(jiān)測能力，這與Wang等[1]的研究結(jié)論一致。

圖4 研究區(qū)Sentinel-1A視線向形變結(jié)果Fig.4 Sentinel-1A line-of-sight deformation results in the study area

表1 水城縣歷史天氣

在研究區(qū)內(nèi)劃定出A、B、C、D、E五個(gè)形變區(qū)，主要分布于北盤江兩側(cè)斜坡體上(圖4黑框)。將各形變區(qū)升降軌視線向形變沿豎直方向、東西水平方向分解，沿坡向合成，結(jié)果如圖5所示。豎直形變量向上(抬升)為正值，向下(沉降)為負(fù)值；水平形變量向西為正值，向東為負(fù)值。各形變區(qū)沿坡向形變量及方向示意圖中箭頭線段長度與坡向形變量成正比。由圖5等高線分布情況可知：1)山坡坡度相對較陡處坡向形變量相對較大；2)地勢平坦地區(qū)DS點(diǎn)密度明顯大于地勢陡峭處DS點(diǎn)密度；3)地勢平坦地區(qū)主要表現(xiàn)為豎直形變，地勢陡峭處同時(shí)存在豎直形變和水平形變。5個(gè)形變區(qū)的潛在威脅目標(biāo)如圖5所示：1)形變區(qū)A。形變分布在斜坡體、采礦區(qū)或加工廠，形變可能與地下采礦和礦物加工的抽排水有關(guān)，主要威脅對象有工廠、居民住房、學(xué)校和道路；2)形變區(qū)B。形變分布在斜坡體，形變可能由斜坡體失穩(wěn)造成，主要威脅對象有居民住房和道路；3)形變區(qū)C。形變分布在斜坡體、加工廠，部分斜坡體具有明顯的老滑坡痕跡，形變可能由礦物加工的抽排水或斜坡失穩(wěn)引起，主要威脅對象有居民住房、工廠、道路和河流；4)形變區(qū)D。形變分布在斜坡體，形變可能由斜坡體失穩(wěn)造成，主要威脅對象有居民住房和道路；5)形變區(qū)E為露天采礦區(qū)，附近有煤炭倉儲(chǔ)中心、煤炭加工廠和發(fā)電廠等，形變可能與采礦和礦物加工的抽排水有關(guān)，主要威脅對象有居民住房、道路和附近工廠。

4 結(jié) 語

使用SBAS InSAR對升軌、降軌Sentinel-1A數(shù)據(jù)進(jìn)行處理，提取貴州省水城縣研究區(qū)5個(gè)明顯形變區(qū)。結(jié)合Google Earth推斷引起地表形變的主要原因有斜坡體失穩(wěn)、地下/露天采礦和礦物加工的抽排水等，潛在威脅對象包括居民住房、學(xué)校、工廠和道路等。雞場鎮(zhèn)滑坡發(fā)生前升軌、降軌Sentinel-1A的SBAS InSAR形變場并未出現(xiàn)明顯形變信息，且在滑坡發(fā)生當(dāng)天及前兩周內(nèi)有11 d出現(xiàn)不同程度降雨，因此初步判斷雞場鎮(zhèn)滑坡是由外部因素(如降雨)觸發(fā)的突發(fā)性滑坡。

SBAS InSAR能有效應(yīng)用于山區(qū)滑坡隱患早期識別和形變監(jiān)測，但對于外部因素(如降雨)觸發(fā)的突發(fā)性滑坡，衛(wèi)星雷達(dá)存在重訪周期較長的局限性，如Sentinel-1A為12 d，Sentinel-1A/B雙星為6 d。因此，隨著衛(wèi)星數(shù)量增加，可通過多星組網(wǎng)構(gòu)建衛(wèi)星星座，縮短衛(wèi)星雷達(dá)的重訪周期，提升突發(fā)性滑坡隱患的早期識別和形變監(jiān)測能力。本研究在貴州水城縣滑坡隱患調(diào)查中所用技術(shù)、方法及實(shí)驗(yàn)結(jié)論可為貴州省以及中國西南山區(qū)滑坡隱患調(diào)查與早期識別提供參考。

致謝：感謝英國利茲大學(xué)Andy Hooper教授提供StaMPS軟件，感謝中山大學(xué)蔣彌教授提供FaSHPS軟件，感謝歐空局提供Sentinel-1A數(shù)據(jù)，感謝美國航空航天局提供SRTM DEM數(shù)據(jù)。

圖5 各形變區(qū)豎直形變、水平形變、沿坡向合成形變量與方向Fig.5 Vertical and horizontal deformation, deformation and direction along the slope in each deformation area