張騰，謝帥，黃波，范景輝，陳建平，童立強(qiáng)

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

0 引言

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

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

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

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)源介紹

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

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

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

表2 Sentinel-1 數(shù)據(jù)信息Tab.2 Sentinel-1 data information

2 技術(shù)原理及數(shù)據(jù)處理

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

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

(1)

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

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

圖2 數(shù)據(jù)處理流程圖Fig.2 Data processing flow chart

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

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

3 結(jié)果與分析

3.1 滑坡探測(cè)空間格局分析

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

(a)ALOS-2數(shù)據(jù)集探測(cè)的活動(dòng)滑坡結(jié)果 (b)Sentinel-1數(shù)據(jù)集探測(cè)的活動(dòng)滑坡結(jié)果

3.2 滑坡探測(cè)結(jié)果驗(yàn)證及分析

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

(a)a處 (b)e處

4 結(jié)論

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

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

2)結(jié)果表明，數(shù)據(jù)集在時(shí)間基線不同的情況下，其探測(cè)能力和優(yōu)勢(shì)是不同的。在保持一定相干性的基礎(chǔ)上，時(shí)間基線越長(zhǎng)，越有利于識(shí)別出具有較小形變速率的活動(dòng)滑坡；時(shí)間基線越短，則更有利于識(shí)別具有較大形變速率的活動(dòng)滑坡。長(zhǎng)、短時(shí)間基線數(shù)據(jù)結(jié)合，能對(duì)滑坡的慢速運(yùn)動(dòng)和快速運(yùn)動(dòng)實(shí)現(xiàn)組合識(shí)別。

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

志謝：衷心感謝JAXA EO-RA2項(xiàng)目(編號(hào)：P3073002)提供ALOS-2數(shù)據(jù)，ESA提供Sentinel-1數(shù)據(jù)。