基于時序InSAR的寧夏西吉縣滑坡災害隱患識別

弓永峰,王 輝,吳學華,張 佳,2*,劉 君

(1.寧夏回族自治區國土資源調查監測院，寧夏 銀川 750004；2.中國地質大學(武漢)環境學院，湖北 武漢 430078；3.長安大學水利與環境學院，陜西 西安 710054)

早期識別是防治地質災害的重要技術手段之一[1]。對于人工調查難以到達的陡峻山區以及植被覆蓋度較高的區域，滑坡災害早期識別可以解決隱患點底數不清的問題，具有良好的地質災害預警和風險管理作用[2]。近年來，合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar,InSAR)技術以其全天候、大范圍、高精度的優勢不斷得到推廣，對于滑坡災害風險評估的發展起著重要的推進作用[3-5]。傳統的時序InSAR技術通過差分合成孔徑雷達干涉測量(Differential Interferometric Synthetic Aperture Radar，D-InSAR)、短基線合成孔徑雷達干涉測量(Small Baseline Subset Interferometric Synthetic Aperture Radar，SBAS-InSAR)、永久散射體合成孔徑雷達干涉測量(Permanent Scatterer Interferometric Synthetic Aperture Radar，PS-InSAR)等對滑坡災害進行識別。其中，D-InSAR技術是對兩景或三景影像數據進行處理，可快速獲取地表形變結果，但無法得到長時序的地表形變信息；SBAS-InSAR技術適用于分布式目標分析，可得到大范圍的面狀地表形變信息；PS-InSAR適用于點目標分析，可提取可靠性高的PS點[6-7]。利用InSAR技術可以高效地進行區域地質災害的早期識別[8]。

InSAR、光學遙感等綜合遙感測量技術的應用對大區域潛在滑坡災害的早期識別及滑坡體風險的評估具有重大的作用[9-12]。如趙寶強[13]通過InSAR技術圈定出白龍江流域部分潛在的滑坡隱患，并驗證了InSAR技術的可靠性；王戰衛[14]基于InSAR和光學遙感技術，識別出青海省大通縣144個潛在的滑坡災害隱患點，表明InSAR技術對大范圍滑坡監測的可行性；黎光艷[15]利用InSAR技術對白龍江中游滑坡進行監測，分析了該地區滑坡分布和形變規律；趙超英等[16]利用SAR數據對甘肅黑方臺滑坡形變進行了監測，研究了滑坡失穩模式；許強[17]提出構建地質災害隱患早期識別的“三查”體系，通過光學遙感、InSAR、無人機航拍等技術在滑坡監測中的綜合應用可取得較好的效果；長安大學張勤教授研究團隊和成都理工大學許強教授研究團隊結合InSAR技術和GNSS監測系統提前6 h成功預警了甘肅永靖縣突發性黃土滑坡，此次成功預警為后續滑坡監測奠定了基礎[18]。盡管InSAR技術在滑坡災害識別和監測方面的應用越來越廣泛，但由于時空失相干的影響，在地形起伏較大且植被覆蓋度高的山區，滑坡災害的識別難度依然很大[19]。因此，InSAR技術針對山區潛在滑坡災害的識別效果和適用性還需要進一步深入分析。本研究的創新在于綜合了SBAS-InSAR技術和PS-InSAR技術在提取地表形變信息方面的優勢，針對寧夏南部山區潛在滑坡災害進行早期識別，有效克服了由于地形地貌所帶來的時間和空間上失相干的局限性，更好地消除了平地效應，得到了精確度更高的地表形變結果，表明了InSAR技術在寧夏南部山區潛在滑坡災害高精度識別和監測的可行性。

寧夏南部山區屬西北黃土高原的一部分，區域內山地與平原錯落分布，高差起伏大，地表溝壑縱橫，是寧夏地質災害的高易發區。西吉縣位于寧夏南部山區，地貌以黃土丘陵為主，次為紅土丘陵，地質環境復雜，地質災害呈易發多發態勢。滑坡占西吉縣地質災害總數的86%，是西吉縣最發育的地質災害。本文以寧夏西吉縣為研究區，首先利用PS-InSAR和SBAS-InSAR融合技術對研究區2018—2021年Sentinel-1A升軌數據進行干涉處理，得到了該區域高精度的地表形變結果；然后通過目視判讀光學影像在地表形變結果中圈定出滑坡災害隱患區域，并通過野外核查，檢驗滑坡災害隱患識別結果的可靠性和有效性；最后分析了典型老滑坡形變時空特征，對典型焦灣滑坡群進行實地勘查，以為寧夏南部山區地質災害隱患早期識別、監測預警以及治理搬遷避讓等地質災害防治工作提供充分有效的數據支撐和方法參考。

1 研究區概況

西吉縣位于寧夏回族自治區南部山區(105°20′E～106°04′E，35°35′N～36°14′N)，縣內交通較便利，東西長約67 km，南北寬約74 km，總面積為3 144 km2，見圖1。西吉縣境內發育清水河、葫蘆河和祖歷河三大季節性河流，根據當地地貌特征可將該地區分為黃土丘陵溝壑區、土石山區和河谷川道3個二級地貌單元，其中黃土丘陵溝壑區占總面積的83%，地勢北高南低、東高西低，海拔在1 656～2 606 m之間。西吉縣生態環境脆弱，年平均氣溫為12.7℃，年平均降水量為570.2 mm，降水主要集中在7～9月份。西吉縣境內土地類型復雜多樣，土層深厚，土質疏松，受1920年8.5級海原大地震的影響，縣境內地震后有大量的新滑坡形成，同時因遭受地震的影響和破壞，該區域局部地區巖土體結構松散，加劇了區域內崩塌、滑坡地質災害發生的可能性[20]。

圖1 西吉縣地理位置圖

2 研究數據與方法

2.1 研究數據

本次研究所需的數據有：Sentinel-1A衛星影像數據、數字高程模型(Digital Elevation Model,DEM)數據、哨兵衛星精密軌道文件數據、光學遙感影像數據。其中，哨兵一號衛星是歐洲航天局哥白尼計劃(GMES)分別于2014和2016年升空的兩顆對地監測衛星，雙衛星以12 d為周期對地球全天候晝夜雷達成像，幅寬為250 km，主要應用于山體滑坡、城市地面沉降監測等方面，哨兵衛星影像數據下載網址為https://search.asf.alaska.edu/，數據中心免費提供同軌道雙衛星影像；DEM(SRTM)數據可提供參考地形，空間分辨率為90 m，其下載網址為https://srtm.csi.cgiar.org/；哨兵衛星精密軌道文件數據能夠對軌道誤差信息進行修正，其下載網址為http://www.gscloud.cn/；光學遙感影像數據為天地圖影像，影像拍攝時間為2021年3月27日，空間分辨率為0.48 m。考慮到數據運行的時間成本，在下載數據時篩選出成像質量較好的數據。本研究選取2018年12月至2021年6月Sentinel-1A衛星影像數據共28景，數據類型為單視復數影像(Single Look Complex，SLC)，數據模式為干涉寬帶(Interferometric Wide，IW)，軌道方向為升軌，空間分辨率為5×20 m、入射角為39.07°。Sentinel-1A衛星影像數據具體成像時間，見表1。

表1 Sentinel-1A衛星影像數據成像時間

2.2 時序InSAR方法

小基線集干涉測量(Small Baseline Subset-Interferometric Synthetic Aperture Radar，SBAS-InSAR)技術是Berardino等[6]在2002年提出的適用于分布式目標的InSAR時序分析方法，主要應用于地表形變時序監測領域的研究。其基本原理是依據同一地區時序上的N+1景SAR影像數據獲得的時空基線設置相應的合理閾值，選出公共主影像，其他影像以主影像為參照組成無數個小基線集合并生成差分干涉圖，并利用最小二乘法及矩陣的奇異值分解法求解地表形變時間序列及形變速率。

假設從研究區域獲取影像時間排列為(t0，t1,…,tn)的N+1幅SAR影像，經過處理可以得到M幅差分干涉圖，則M滿足：

(1)

假設對不同時刻ta和tb獲取的兩幅影像進行干涉處理得到差分干涉圖j(tb>ta)，則差分干涉圖j對應某一像素點(方位向x，距離向r)的相位值可以表示為

δφj(x,r)=φ(tb,x,r)-φ(ta,x,r)

(2)

式中:φ為干涉相位;j∈(1，2，…，M)；λ為雷達波長；d(tb,x,r)和d(ta,x,r)為時間ta和tb以初始時間t0為對照[d(t0,x,r)=0]沿雷達視線(LOS)方向累計的地表形變量。

為了獲取符合規律且更為直觀的地表沉降序列，生成第j幅差分干涉圖的相關時段與平均相位速度進行求積，可計算出該幅差分干涉圖的相位值，則平均相位速度與相位分別表示如下：

(3)

(4)

式中：v為平均相位速度；φ為干涉相位；t為相位獲取時間；δφj為第j幅差分干涉圖的干涉相位；k∈(1,2,…,M)。

對各個時間段地表沉降速度進行積分，即可獲取各個時段的地表沉降量，具體表示為矩陣積分形式：

Bv=δφ

(5)

式中當系數矩陣B(M×N)滿秩或秩虧時，研究區域地表形變速率可分別利用最小二乘法和奇異值分解法求解，再根據地表形變速率求出相應的累積地表形變量。

永久散射體干涉測量(Persistent Scatterers-InSAR，PS-InSAR)技術是由Ferretti等[21]提出的針對具有高相干性的點目標進行差分干涉測量，其基本原理是選取同一地區時間連續的多景影像數據，通過設定適當的基線閾值從影像中選取1景為主影像，其余為從影像，并將所有影像與主影像配準后生成時間序列干涉相對，從研究區域中獲取穩定的具有高相干性的永久散射體點，再使用DEM對干涉相對進行差分干涉處理，得到每個永久散射體點的差分干涉相位，其又包含相位的4個分量：

φdint=φdem+φdef+φatm+φnoise

(6)

式中:φdint為每個永久散射體點的差分干涉相位;φdem為DEM誤差造成的地形相位;φdef為雷達視線方向(LOS)地表形變引起的形變相位;φatm為大氣延遲引起的大氣相位;φnoise為噪聲引起的噪聲相位。

通過對永久散射體進行時間序列分析，根據各相位分量特性將差分干涉相位中大氣、噪聲、DEM等相位分離出來，獲取到每個永久散射體點的地表形變相位，從而準確地監測出研究區域地表位移變化值。通過該方法識別出SAR影像中后向散射特征穩定的點，可以避免人為選取過程中存在不穩定GCP點帶來的誤差。

2.3 山區地表形變監測

山區地形復雜，溝壑縱橫，高差起伏較大，植被生長繁茂，導致傳統調查監測手段無法獲取區域性的地表損傷，從而影響對地質災害危險性和風險性的判別。雖然PS技術局限于線性形變，但可以獲取相干性高的GCP點；而SBAS技術在分析分布式目標方面具有很好的應用效果。因此，本文基于兩種方法的優勢互補，對PS-InSAR、SBAS-InSAR技術進行融合用于獲取地表形變監測結果，相比分別利用PS-InSAR或SBAS-InSAR技術可有效減小SAR數據在處理過程中造成的誤差，增加結果的準確性。基于SBAS-InSAR技術處理數據過程中，在軌道精煉時導入了利用PS-InSAR技術獲取的高相干性的GCP點，以此來去除殘余的相位信息，更好地消除了平地效應，提高了SBAS-InSAR技術處理結果的精確性。其具體操作步驟如下：

(1) 數據預處理：由于原始數據量較大，通過對研究區進行裁剪，可減少數據處理時間，提高處理效率。

(2) 利用PS-InSAR技術提取高相干性的GCP點：設置形變采樣頻率、殘余高度采樣頻率、時空相關濾波、相干系數閾值等參數，減小大氣相位誤差，獲取高質量、穩定的GCP點。

(3) 利用SBAS-InSAR技術進行差分干涉處理：通過設置一定的時空基線閾值，生成干涉相對，再進行配準、去平、濾波、相位解纏及相干計算，得到去平和濾波后的相位圖及解纏圖。

(4) 軌道精煉和重去平：導入PS-InSAR技術獲取的GCP點，消除平地效應。

(5) 時間/空間域形變計算：通過兩次SBAS反演有效去除大氣相位，得到形變速率及形變量。

融合PS技術的SBAS-InSAR技術流程如圖2所示。

圖2 融合PS技術的SBAS-InSAR技術流程圖

3 研究結果與分析

3.1 潛在滑坡災害隱患識別

利用SARScape 5.2軟件對研究數據進行處理，將結果導入Arcgis 10.3軟件中進行密度分割，并添加光學遙感影像底圖，最終得到西吉縣雷達視線(LOS)向年平均地表形變速率在空間上的變化圖。根據時序InSAR技術監測結果顯示，相干點目標(CTs)為824 557個，平均相干點密度約為262 CTs/km2，能夠滿足與滑坡有關方面研究的要求[11]。根據斜坡的巖性特征、數據反演精度和實地調查結果，以及在前人研究基礎上[12]，設置±10 mm/a為雷達視線向地表形變速率閾值，借助光學影像對潛在的滑坡進行識別，得到西吉縣雷達視線向年平均地表形變速率圖和識別的潛在滑坡災害隱患點，如圖3所示。圖中地表形變速率正值說明地表運動相對衛星方向靠近；地表形變速率負值說明地表運動相對衛星方向遠離。

圖3 西吉縣雷達視線(LOS)向年平均地表形變速率圖和識別的潛在滑坡災害隱患點

由圖3可知：研究區內共識別出11個特大潛在的滑坡災害隱患點，其中有3個歷史崩塌，均屬于斜坡在自然和人工條件下遭受破壞變形，工程地質學上統稱其為不穩定斜坡類地質災害[22-23]。大氣降水的頻次以及地殼活動能夠在一定程度上誘發滑坡、泥石流等斜坡類地質災害[24-26]。據調查，1996—1998年間研究區強降雨的頻次高，發生滑坡、崩塌的點多，完全符合此時段內降水頻次高且強度大的特點。西吉縣境內地層以第四系黃土為主，由于黃土的垂直節理非常發育，斜坡邊緣常發育平行于斜坡走向的張裂縫帶，寬數米，加之黃土結構松散、孔隙率高以及黃土特有的濕陷性等，尤其在強降雨時，地表水入滲速度很快，土體含水量較大，其抗拉、抗剪強度大大降低，連續性降水導致地表水持續對斜坡的不斷滲透，最終造成滑坡災害發生。

3.2 典型潛在滑坡形變的時空特征分析

3.2.1 向家村老滑坡

向家村老滑坡位于西吉縣城西南方向約26 km處，其光學影像圖和InSAR結果圖如圖4所示。圖中紅色界線區域為滑坡范圍；黃色界線區域代表滑坡早期識別加劇范圍。

圖4 向家村老滑坡光學影像圖和InSAR結果圖

由圖4可以看出，向家村老滑坡滑坡體呈不規則長舌狀，周界清晰可見，后緣發生明顯形變，滑坡體中下部有耕地分布，滑坡下部有道路橫穿。

圖5為向家村老滑坡的地表形變速率分布圖。

圖5 向家村老滑坡的地表形變速率分布圖

由圖5可知，向家村老滑坡最大地表形變速率位于滑坡體形變區后緣，其值達到-17 mm/a，滑坡范圍超出原有滑坡邊界。

為了對向家村老滑坡滑坡體的形變趨勢進行深入研究，選取1#～5#5個樣點(具體位置如圖5標記所示)形變進行了時序分析，得到向家村老滑坡形變的時間序列圖，見圖6。其中，1#樣點位于該滑坡體后緣，該處發生明顯形變，累計地表形變量達-50 mm；2#樣點位于該滑坡體后側邊緣，自2018年12月起該樣點形變持續呈下降趨勢，累計地表形變量達-30 mm；3#樣點位于滑坡體中間位置，該樣點形變趨于加速狀態，累計地表形變量為10.21 mm；4#樣點位于該滑坡治理區左側，監測期間該樣點累計地表形變量僅為-2 mm；5#樣點位于該滑坡治理區右側，該樣點形變在監測期間內基本趨于穩定狀態。這是因為：1#樣點、2#樣點位于該滑坡后緣，由于連續降雨造成松散物源下移，導致滑坡后緣形變明顯；3#樣點處緩慢抬升，主要是由于該滑坡后緣部分松散巖土下滑在該處堆積所致；4#樣點、5#樣點位于滑坡治理區附近，由于人為活動的影響，使該處地表形變趨于穩定狀態。通過有效整治后該滑坡體前緣處地表形變基本保持穩定狀態，而滑坡體中后緣地表形變趨勢仍呈加速狀態且累計地表形變量已較大，一旦平衡被打破，滑坡體下方居民的生命和財產安全將無法得到保證，且還可能會影響道路交通安全。

圖6 向家村老滑坡的形變時間序列圖

3.2.2 焦灣滑坡

焦灣滑坡為大型滑坡群，位于西吉縣平峰鎮焦灣村焦灣組，該處滑坡主要由1920年、1970年和2008年發生在濫泥河左岸黃土斜坡上的幾次大地震誘發所致，其滑坡體較陡，上覆黃土，土體破碎，落水洞極其發育，滑坡前緣面臨沖溝溝底有少量的新近系紅泥巖出露。該地區為典型的黃土丘陵溝壑區，黃土梁整體走向西北，地勢由東南向西北傾斜，峁梁起伏，溝壑縱橫，滑坡發育于黃土丘陵斜坡上。根據InSAR結果顯示，該大型滑坡群在監測時間內又產生了新的形變，具有老滑坡復活的可能性。結合光學遙感解譯和InSAR技術對焦灣滑坡災害隱患點進行識別，選取焦灣滑坡群中有形變的區域進行分析， 其結果如圖7所示。圖中紅色邊界內代表滑坡范圍。通過對InSAR得到的量化結果進行分級并用不同的顏色顯示，圖7(b)中藍色部分表示滑坡形變區域，6#樣點處地表形變速率為16 mm/a，7#樣點處地表形變速率為13 mm/a。7#樣點位于焦灣滑坡邊界外部，隨時有滑塌的危險。

圖7 焦灣滑坡光學影像圖和InSAR結果圖

為了更加直觀深入地了解焦灣滑坡群形變情況及滑坡成因，對該滑坡展開了野外實地勘查。野外實地勘查結果顯示：該滑坡位于焦灣村西北部，地理位置為東經105°33′34″、北緯35°47′53″，滑坡全貌如圖8(a)所示；該滑坡周界清晰，呈圓弧狀，后壁較陡，高程為2 060 m，高約25 m，滑坡體較陡，坡度近65°，滑坡體上土體較為破碎，現已被人工改造為耕地[見圖8(b)]；滑坡體遭流水侵蝕又出現新的侵蝕溝，使滑坡體穩定性下降，坡壁發育有落水洞[見圖8(c)]；該滑坡周邊路面塌陷現象嚴重[見圖8(d)]；該滑坡滑體土主要為粉土，顏色為黃褐—黑褐色，均質結構，土層松散，抗剪強度低，垂直節理發育，含粉砂、暗色礦物，發育蟲孔、根孔，巖土體含水率較低，處于可塑狀態；該滑坡滑床主要巖性為晚更新世黃土，黑褐色，巖土體含水率較高，垂直節理發育。焦灣滑坡屬于以地震為主要觸發因素的滑坡，該滑坡前緣受河流沖刷作用，在雨季易發生蠕滑；滑坡區地層主要為更新統風積黃土和第四系滑坡堆積物，巖土體裂隙發育，物理力學性質較差，土體較破碎，降水易入滲，可能會引起滑坡整體失穩[27]。因此，對焦灣滑坡的防治應以監測為主，如有條件可在該滑坡體前緣輔以防護工程，以減少河流的沖蝕作用，一旦發現有形變跡象，立即采取防災措施。

圖8 焦灣滑坡的野外實地勘查結果

經分析，焦灣滑坡滑坡體裂縫、落水洞發育，應該及時填壓，以防止雨水入滲；裸露的坡面易于地表水的下滲，應合理植樹造林，改善坡面植被條件；滑坡體前緣存在臨空面，應盡量避免可能擾動滑坡的各種人類工程活動，并嚴禁開挖坡腳；合理利用滑坡土地，不宜在滑坡上興建有利于降水入滲的各種水利工程；建立滑坡災害監測網絡，加強雨季老滑坡動態的監測，發現形變跡象應立即上報主管部門，以便制定防災應急對策。

4 結 論

本文基于哨兵Sentinel-1A衛星影像數據，以寧夏南部山區的西吉縣滑坡為研究區域，實施了全面監測和深入分析，對篩選出的部分滑坡災害隱患監測點形變進行時間序列分析，同時對滑坡災害隱患識別監測過程中時序InSAR技術的應用方法進行了研究，得到結論如下：

(1) 依據PS和SBAS融合技術，對2018—2021年期間寧夏西吉縣累計地表形變量及形變速率進行監測，結果顯示研究區域內最大累計地表形變量為-55 mm，其中西南側的滑坡體為西吉縣境內地表形變較大的區域。

(2) 根據西吉縣雷達視線向地表形變速率圖和光學影像進行潛在的滑坡災害早期識別，在西吉縣周邊覆蓋范圍內識別出11個特大潛在的滑坡災害隱患點，其與實地勘查結果高度吻合。

(3) 對典型地質災害點向家村老滑坡形變進行時序分析表明，在2019年7月至2021年3月監測期內，向家村老滑坡下部治理區地表形變保持平穩，滑坡體上部累計地表形變量增加明顯，需要重視監測和防災減災工作。

(4) 焦灣滑坡目前處于基本穩定狀態，但在雨季尤其是大雨季節常出現蠕滑現象，坡后容易積水，應修建排水溝渠，排截地表水。

致謝:感謝中國自然資源航空物探遙感中心提供基礎數據；ASF DAAC(https://asf.alaska.edu/)為本文提供了歐空局哥白尼Sentinel-1 數據；本文采用了SARscape、ArcGIS等軟件進行數據處理，在此一并表示感謝！