張騰，謝帥，黃波，范景輝，陳建平，童立強

(1.中國地質大學(北京)地球科學與資源學院，北京 100083；2.北京市交通委員會，北京 100073；3.河北省水文工程地質勘查院，石家莊 050021；4.中國自然資源航空物探遙感中心，北京 100083)

0 引言

在中國，山地和丘陵地區(qū)占據國土面積的65%，而滑坡是山地丘陵地區(qū)最常見的自然地質災害?；轮饕怯苫顒訕嬙旌蜌v史地震影響的內部動力地質作用和降雨等外部自然因素兩方面控制產生的[1]，滑坡一旦形成，就會對居住在山區(qū)的人們有著致命的影響，因此提前預警滑坡的形成對防止危險事故的發(fā)生至關重要，而斜坡的運動特征是滑坡最重要的指標之一，是滑坡發(fā)生的先兆。常規(guī)滑坡災害點監(jiān)測(如水準和全球定位系統(tǒng)(global positioning system，GPS)測量)要耗費大量的人力物力，十幾個到上百個觀測點的布置需要耗費大量的時間，地理條件偏遠和自然環(huán)境惡劣的山區(qū)等地往往無法開展專門的地面觀測。合成孔徑雷達差分干涉測量技術(differential interferometric synthetic aperture Radar,DInSAR)的出現有效改善了這一難題。在滑坡災害的相關研究中，范景輝[2]利用合成孔徑雷達(synthetic aperture Radar,SAR)數據、反射體及GPS測量對范家坪滑坡變形進行監(jiān)測，綜合分析滑坡影響因素；Liu等[3]和王振林等[4]用ALOS-2數據對中國西南部的滑坡變形進行監(jiān)測，取得了不錯的成果；王振林等[4]和Intrieri等[5]利用Sentinel-1數據對滑坡變形進行時序性分析監(jiān)測。

四川省茂縣地區(qū)位于西南山區(qū)，地質災害頻發(fā)，滑坡災害造成巨大人員傷亡。本文綜合過往研究經驗，獲取了2015—2017年長波長L波段的ALOS-2和短周期觀測的Sentinel-1兩種不同分辨率的遙感衛(wèi)星數據，并結合資源三號(ZY-3)衛(wèi)星制作的高分辨率5 m數字高程模型(digital elevation model，DEM)數據和航天飛機雷達地形測繪(shuttle Radar topography mission，SRTM)30 m分辨率DEM數據，使用兩軌法和Stacking方法探測茂縣中部地區(qū)的滑坡形變，分析茂縣滑坡的形變特征，為地區(qū)的防災減災提供相關參考。

1 研究區(qū)概況及數據源

1.1 研究區(qū)概況

四川省茂縣位于四川省阿壩藏族羌族自治州，是中國西南地區(qū)典型的山區(qū)縣(圖1)。縣域內跨越岷江和涪江上游高山河谷地帶，發(fā)育著茂汶斷層和龍門山斷層，斷裂順著岷江河谷發(fā)育，存在著明顯的破碎帶和脆弱的地質條件。當地的年降雨量集中于5—10月，分布不均，最大日降雨量達到75.2 mm，瞬時降雨量大，很容易導致地質災害的發(fā)生[6]。地區(qū)內以千枚巖、板巖的變質巖為主，而在岷江流域的河谷沖刷沿岸坡上堆積了大量的松散土塊，基礎設施的防治條件較差，極易造成土塊崩塌。沿岸半山腰一旦崩塌，地形較陡，落差較大，很容易發(fā)展成巨型、巖質或基巖滑坡[1]。地區(qū)處于南北方向地震帶即松潘和龍門山地震帶的中部地區(qū)，同時九頂山華夏系構造帶從南向北斜貫該地區(qū)，受構造因素運動影響較劇烈。地區(qū)內地震現象常有發(fā)生，1933年疊溪7.5級地震和2008年的汶川8.0級大地震更是減弱了茂縣地區(qū)內陡坡的穩(wěn)定性，在暴雨等外部誘發(fā)因素下，區(qū)內的滑坡運動極有可能加劇，需要特別防范。2017年6月，茂縣疊溪鎮(zhèn)新磨村突發(fā)山體高位垮塌，造成河道堵塞，100余人被掩埋。當地山高坡陡，位于流域上游而建設的公路、水電、水利設施等人工建設項目，更是有著巨大的隱患，而對山區(qū)游客、當地居民和通行車輛存在著潛在的人身安全威脅[7-9]。在自然和人為的綜合作用下，該區(qū)域為重要地質災害隱患地區(qū)，區(qū)內的白布村滑坡群便是典型案例。研究區(qū)地形和數據覆蓋情況如圖1所示，圖中紅色框表示具體研究區(qū)域，水藍色框表示Sentinel-1數據挑選的burst覆蓋范圍，黃色和深藍色框分別代表ALOS-2 SM2級別和SM3級別數據覆蓋范圍。

圖1 研究區(qū)地形及ALOS-2、Sentinel-1數據覆蓋范圍Fig.1 Terrain of the study area and ALOS-2 and Sentinel-1 data coverage

1.2 數據源介紹

ALOS-2衛(wèi)星是由日本航空航天局(Japan Aerospace Exploratin Agency,JAXA)于2014年5月發(fā)射的高級陸地觀測衛(wèi)星。其上搭載PALSAR-2傳感器,工作波段為L波段(1.2 GHz段),最高分辨率為：3 m(距離向)×1 m(方位向)，能在復雜大氣條件下全天候工作，最大觀測范圍為490 km×350 km。相對于短波波段的SAR，ALOS-2能更好地去除植被在斜坡上的去相關效應，更加準確地對地表信息進行觀察，及時監(jiān)測到一些自然災害，例如滑坡。本文使用的ALOS-2數據參數如表1所示。因為本次實驗采集到的ALOS-2數據量較少，無法使用多時相InSAR技術，便采用常規(guī)DInSAR技術中的兩軌法探測地表形變。

表1 ALOS-2 PALSAR-2數據信息統(tǒng)計Tab.1 ALOS-2 PALSAR-2 data information

Sentinel-1衛(wèi)星由歐空局(European Space Agency,ESA)運營，搭載了C波段(5.6 cm)傳感器，獲取數據時間上具有穩(wěn)定的連續(xù)性，其IW成像模式，幅寬250 km，地面空間分辨率5 m(距離向)×20 m(方位向)，采用新型循序掃描地形觀測(terrain observation with progression scans，TOPS)成像技術，數據量豐富，下載方便。研究采集了45景Sentinel-1升軌數據，日期為2015年11月26日—2017年12月21日，采用Stacking技術探測地表形變。Sentinel-1數據參數如表2所示。在本次研究中，DEM采用了兩種不同分辨率的數據。一個為基于我國的ZY-3的兩景立體像對通過ENVI和PCI軟件生成的5 m DEM數據。ZY-3衛(wèi)星配備三線陣測繪相機和多光譜相機，全色波段分辨率2.1 m，幅寬52 km[10]，具有廣泛的覆蓋范圍，受大氣和地形的影響較小，適合進行偏遠地區(qū)的測圖工作。在本次研究中用ZY-3衛(wèi)星制作的5 m分辨率的DEM輔助ALOS-2數據實施圖像配準和地理校正[11]，在InSAR差分干涉的過程中，高精度的高程數據能夠提高配準精度、削弱高程殘余相位，5 m的數字表面模型能更好地滿足研究的需要。另一個的DEM地形數據是由美國國家地理空間情報局(National Geospatial-Intelligence Agency，NGA)和美國國家航空航天局(National Aeronautics and Space Administration，NASA)免費提供的30 m分辨率的SRTM-DEM數據，被用于輔助更大范圍內Sentinel-1 數據的圖像配準和地理校正。

表2 Sentinel-1 數據信息Tab.2 Sentinel-1 data information

2 技術原理及數據處理

DInSAR兩軌法是對同一地區(qū)的兩景SAR影像進行干涉差分。它需要一個外部DEM來模擬和消除地形的影響，外部DEM被轉為地形相位后，地形因素被消除，從而得到形變相位。兩軌法首先需要對外部DEM 和主影像進行精確匹配，把DEM 的高程模擬成地形干涉條紋形式(即模擬地形相位)，然后就可以從干涉圖中移除DEM 模擬的地形相位，經過濾波及相位解纏后，從而得到形變信息。

Stacking技術對生成的多景差分干涉圖進行解纏后估算出其線性相位速率[12]，可以在較少的時序數據基礎上獲取年平均形變速率，一定程度上減少大氣誤差、提高形變精度，從而完成滑坡的探測。其原理是基于最小二乘法對N組觀測值線性回歸，估算公式為：

(1)

式中：Δtj為第j組干涉圖的時間基線；φj為第j組解纏的差分干涉相位；ph_rate為線性相位速率；N為生成的差分干涉圖的數量。

InSAR數據處理的具體步驟如圖2所示。第一步，ALOS-2數據配準是將其與高分辨率的DEM進行地理編碼后生成初始查找表，計算軌道偏移量生成多項式精化配準查找表，多次重采樣后進行亞像元級配準；Sentinel-1數據配準是將采集的多個SAR數據選出覆蓋研究區(qū)范圍的條帶，按照研究區(qū)范圍進行拼接、裁剪，隨后選擇一張圖像作為主圖像，利用DEM地形數據和精密的軌道參數對SAR數據進行初步配準，再建立主輔影像配準的查找表進行精確配準，最后通過強度互相關最大化的方法估算殘余偏移量，精確配準后再基于頻譜多樣性進行重疊區(qū)配準。第二步，分別將配準后的兩種數據進行影像多視處理，ALOS-2的多視比為4∶2，Sentinel-1的多視比為5∶1。第三步，將Sentinel-1數據通過限定時間基線和空間基線的原則生成67對干涉對，進行像對差分干涉處理(空間基線限定150 m，時間基線限定24 d)；選擇相同級別的ALOS-2數據兩兩組合生成干涉對。第四步，將ALOS-2數據干涉對進行自適應濾波，以及去除大氣延遲相位，選擇相干圖中相關系數較高、位于穩(wěn)定區(qū)域的像元作為參考點以進行相位解纏，便得到了研究區(qū)的形變相位。第五步，針對Sentinel-1數據將多個差分條紋圖疊加求平均后進行濾波，并利用最小費用流算法進行相位解纏，然后對多景解纏的差分干涉圖根據式(1)計算像元的線性形變速率。相比于DInSAR技術，Stacking方法可以有效地抑制大氣效應與DEM誤差的影響。

圖2 數據處理流程圖Fig.2 Data processing flow chart

ALOS-2數據和Sentinel-1數據的干涉對生成的差分干涉圖如圖3所示。從圖3可以看出，ALOS-2的L波段比Sentinel-1的C波段，干涉相位受地表植被影響更小，克服時間去相干的能力更佳，干涉相位信息更加豐富。

(a)ALOS-220160306-20160320干涉對 (b)Sentinel-120160301-20160325干涉對

3 結果與分析

3.1 滑坡探測空間格局分析

InSAR形變結果通常采用水準儀或GPS觀測數據進行驗證。研究區(qū)位于陡峭的山區(qū)，地理條件惡劣，山高險峻，無法提供精準的地面測量數據來定量驗證InSAR形變結果。本研究將采用光學影像特征和已有歷史調查資料作為輔助驗證數據，對ALOS-2與Sentinel-1兩種數據的InSAR探測結果進行比較和交叉驗證，以間接驗證InSAR探測滑坡形變的可靠性。兩種數據集探測的形變結果如圖4所示。圖4(a)顯示了基于ALOS-2數據集獲取的雷達視線向(LOS)年均地表形變信息，形變速率結果集中于-80～-10 mm/a。如圖4(b)所示，基于Sentinel-1A數據集提取出的點目標分布較均勻密集，形變速率結果分布集中于-90～-10 mm/a。對比兩種數據集探測的形變結果可以得到初步結論，研究區(qū)出現了多個滑坡變形中心，主要分布于研究區(qū)的中部，山谷的兩側。形變方向主要以兩岸向河流滑動為主，形變速率空間分布不均一，從邊緣到中心逐漸增加，呈現圓形或橢圓形。由于本次實驗采用的是升軌SAR數據在河谷西側形變信息比東側更為豐富，主要原因是河谷東側斜坡受壓縮、疊掩等幾何畸變影響導致可靠結果難以獲取[13]。綜合分析兩種數據集的形變速率結果，圈定了8個活動滑坡，分別為黃草坪滑坡a，白布村滑坡b和c，窩窩店滑坡d和e，蘇家坪滑坡f，水草溝滑坡g，立角村滑坡h。兩種數據的探測結果在其中7個活動滑坡處的空間分布能互相對應，分別是a，c，d，e，f，g和h區(qū)域。從圖4(b)可以看出，a和b區(qū)域的形變量較大，形變速率最大絕對量達到了200 mm/a。綜合以上可以得出，采用少量長波段ALOS-2影像和DInSAR兩軌法，能夠在具有一定植被覆蓋度的山區(qū)探測到較為明顯的滑坡地表形變；而基于長時序Sentinel-1影像和Stacking時序分析方法，探測效果比采用ALOS-2影像的兩軌法更佳。

(a)ALOS-2數據集探測的活動滑坡結果 (b)Sentinel-1數據集探測的活動滑坡結果

3.2 滑坡探測結果驗證及分析

對比兩種數據得到的探測結果，可以看出圈定的活動滑坡空間分布能互相對應，側面驗證了本實驗的結果準確性。不過因為不同傳感器、入射角、以及時間、空間失相關等影響，兩個數據集得到的形變結果在形變值有略微差異。ALOS-2數據結果未探測到a和b活動滑坡區(qū)的原因可能是因為數據集具有較長的空間基線，區(qū)域內年變化較大，數據相干性差，導致區(qū)域內形變點太少甚至沒有點。c，d，e，g，h等區(qū)域Sentinel-1A的形變結果值分布較小，這可能是由于采集45景2015年11月26日—2017年12月21日之間的線性年平均形變速率小于ALOS-2數據采集時段的線性形變速率。本研究對圈定的a，b，c，d，e，g等6處進行了野外地質調查。調查發(fā)現，a，b，c，d，e等遙感識別的活動滑坡存在明顯的地表形變，滑坡堆積體、山坡小范圍坍塌、公路裂縫、陡坎等滑坡形變特征在多處出現，對于一些地勢起伏較大的區(qū)域g，野外調查時只能遠眺，發(fā)現河溝對岸山坡小范圍滑塌，有大塊板狀碎石。比如，從圖5(a)可以看出在圈定活動滑坡區(qū)域黃草坪a的右邊緣存在拉裂縫、堆積物，從圖5(b)可以看出在圈定的活動滑坡區(qū)域窩窩店滑坡e前緣，存在不斷滑動的跡象，右側拉裂縫明顯，有潛在危險，整體判斷為古滑坡復活。

(a)a處 (b)e處

4 結論

滑坡運動是滑坡災害發(fā)生的先兆，為防治災害，探測早期的滑坡形變分布至關重要。而不同的滑坡發(fā)育階段也有著不同的形變特征，形變速度也各有不同，所以就需要運用多種數據進行組合探測。

1)本文利用高分辨率的ALOS-2數據和5 m分辨率的DEM數據探測茂縣白布村的滑坡形變，高分辨率的DEM能夠輔助實現InSAR影像的高精度配準，并削弱高程殘余誤差，為應用DInSAR方法獲得更為可靠的形變速率圖提供了有利條件。而重訪周期較短的Sentinel-1數據基于Stacking技術生成的形變速率圖與ALOS-2數據表現了空間分布上的一致性。野外調查顯示，在依據兩種InSAR結果圈定的典型活動滑坡上，存在實地變形特征。

2)結果表明，數據集在時間基線不同的情況下，其探測能力和優(yōu)勢是不同的。在保持一定相干性的基礎上，時間基線越長，越有利于識別出具有較小形變速率的活動滑坡；時間基線越短，則更有利于識別具有較大形變速率的活動滑坡。長、短時間基線數據結合，能對滑坡的慢速運動和快速運動實現組合識別。

3)結果顯示，采用少量長波段ALOS-2影像和DInSAR兩軌法，能夠在具有一定植被覆蓋度的山區(qū)探測到滑坡地表形變；當應用C波段的Sentinel-1影像時，通過積累連續(xù)的多景數據，并應用時序分析方法，探測效果比采用ALOS-2影像的兩軌法效果更佳。

志謝：衷心感謝JAXA EO-RA2項目(編號：P3073002)提供ALOS-2數據，ESA提供Sentinel-1數據。