張文龍，張童康

(1.自然資源陜西衛(wèi)星技術(shù)應(yīng)用中心，陜西 西安 710100； 2.中國(guó)煤炭地質(zhì)總局航測(cè)遙感局，陜西 西安 710100)

我國(guó)山地丘陵區(qū)約占國(guó)土面積的65%，地質(zhì)條件復(fù)雜，是世界上地質(zhì)災(zāi)害最嚴(yán)重、受威脅人口最多的國(guó)家之一[1-2]。近年來(lái)，我國(guó)地質(zhì)災(zāi)害處于多發(fā)態(tài)勢(shì)，主要的地質(zhì)災(zāi)害有崩塌[3]、滑坡[4]、泥石流和地面塌陷[5]。地質(zhì)災(zāi)害多發(fā)地區(qū)主要有陜北黃土高原地區(qū)[6]、山西采煤沉陷區(qū)[7]、三峽庫(kù)區(qū)、滇西北和藏東南等地區(qū)。地質(zhì)災(zāi)害防治形勢(shì)嚴(yán)峻，給人民生命安全和經(jīng)濟(jì)帶來(lái)嚴(yán)重威脅[8]。目前，利用InSAR技術(shù)對(duì)地質(zhì)災(zāi)害早期識(shí)別的研究取得了一些成功范例。張路等[9]采用相干散射體時(shí)序InSAR方法，識(shí)別出了17處持續(xù)變形中的不穩(wěn)定坡體；戴可人等[10]利用InSAR成功對(duì)潛在的滑坡進(jìn)行了早期識(shí)別；Grabriel等[11]利用D-InSAR技術(shù)探測(cè)厘米級(jí)的地表形變，測(cè)量了美國(guó)加利福尼亞州東南部的英佩瑞爾河谷灌溉區(qū)的地表形變；黃潔慧等[12]采用二軌D-InSAR技術(shù)對(duì)米林滑坡進(jìn)行形變監(jiān)測(cè)，得到米林滑坡災(zāi)前形變趨勢(shì)；張曉倫等[13]將D-InSAR技術(shù)和光學(xué)遙感技術(shù)結(jié)合對(duì)東川城區(qū)周邊滑坡災(zāi)害進(jìn)行了有效識(shí)別。傳統(tǒng)的D-InSAR技術(shù)能夠有效探測(cè)地表微小變形，但是容易受空間和時(shí)間失相干的影響，大多數(shù)干涉圖中都會(huì)存在大區(qū)域信號(hào)不相關(guān)或沒(méi)有變形信息的現(xiàn)象，而且在最終的變形相位中往往包含明顯的大氣延遲相位，嚴(yán)重影響地表變形測(cè)量的精度和效率[14]。為了剔除這些誤差，提高監(jiān)測(cè)精度，國(guó)外學(xué)者Ferreti等[15-16]提出了PS-InSAR和SBAS-InSAR技術(shù)為代表的時(shí)間序列InSAR技術(shù)。

代志宏等[14]采用PS-InSAR技術(shù)對(duì)南寧地表沉降進(jìn)行了監(jiān)測(cè)與分析，該研究可為地質(zhì)災(zāi)害預(yù)測(cè)和防治工作提供新思路。牛全福等[17]將PS-InSAR技術(shù)成功應(yīng)用于蘭州市地面形變的監(jiān)測(cè)。馮文凱等[18]利用SBAS-InSAR技術(shù)對(duì)金沙江流域沃達(dá)村滑坡進(jìn)行地表形變監(jiān)測(cè)，為類似老滑坡監(jiān)測(cè)預(yù)警提供了新的思路與借鑒。Liu Peng等[19]利用SBAS-InSAR技術(shù)處理了三峽庫(kù)區(qū)2003—2010年間EnvisatSAR圖像，識(shí)別出巴東縣2個(gè)山體滑坡，認(rèn)為季節(jié)滑坡運(yùn)動(dòng)與水位變化之間存在明顯的相關(guān)性。張?jiān)娗训萚20]嘗試采用小基線方法對(duì)一定區(qū)域范圍內(nèi)的歷史滑坡進(jìn)行監(jiān)測(cè)分析，為區(qū)域監(jiān)測(cè)提供了有益借鑒。可見(jiàn)，SBAS-InSAR方法不僅可以在區(qū)域尺度上獲取長(zhǎng)時(shí)序的緩慢地表形變信息，而且在潛在滑坡識(shí)別、已知滑坡監(jiān)測(cè)等方面具有較大的應(yīng)用潛力[21]，但是其在地質(zhì)災(zāi)害隱患大范圍早期識(shí)別中應(yīng)用較少。本文采用SBAS-InSAR技術(shù)結(jié)合高分影像對(duì)大范圍地質(zhì)災(zāi)害易發(fā)區(qū)進(jìn)行早期識(shí)別，通過(guò)SBAS-InSAR技術(shù)提取形變區(qū)，然后利用高分影像進(jìn)行遙感解譯，剔除不可靠的形變區(qū)，有效提高了地質(zhì)災(zāi)害早期隱患識(shí)別的精度。

1 SBAS-InSAR技術(shù)

SBAS-InSAR技術(shù)主要去除或減弱了時(shí)空失相干等誤差的影響，假設(shè)在時(shí)間段內(nèi)，獲取S景SAR影像，選取主影像并與從影像進(jìn)行配準(zhǔn)。然后進(jìn)行差分干涉處理，一共得到N幅差分干涉圖，且滿足：

(1)

對(duì)于tA和(tA

(2)

假定不同干涉圖間的形變速率為υk，k+1，則tA～tB間的累積形變可表示為：

(3)

對(duì)N幅干涉條紋圖進(jìn)行三維時(shí)空相位解纏，即可求出不同SAR獲取時(shí)間的形變速率。

時(shí)間序列InSAR技術(shù)流程[23]如圖1所示。

圖1 SBAS-InSAR數(shù)據(jù)處理流程Fig.1 Data processing flow of SBAS-InSAR

2 研究區(qū)概況與實(shí)驗(yàn)數(shù)據(jù)

2.1 研究區(qū)概況

子長(zhǎng)市位于黃土高原中部，總面積為2 405 km2[24]，具體位置如圖2所示。子長(zhǎng)市屬典型的黃土高原丘陵溝壑區(qū)，基巖產(chǎn)狀近于水平，無(wú)明顯褶皺和斷裂，而節(jié)理裂隙構(gòu)造較為發(fā)育。境內(nèi)峁梁起伏，溝壑縱橫，地面支離破碎，崩塌、滑坡等地質(zhì)災(zāi)害頻繁發(fā)生[25]。

圖2 研究區(qū)位置Fig.2 Location of study area

2.2 實(shí)驗(yàn)數(shù)據(jù)

研究采用Sentinel-1數(shù)據(jù)作為InSAR監(jiān)測(cè)主要數(shù)據(jù)源，研究區(qū)范圍為Path/frame：11/116，共獲取39景Sentinel-1數(shù)據(jù)，部分?jǐn)?shù)據(jù)參數(shù)如下。

數(shù)據(jù)獲取日期共有10日，激化方式VV，入射角39.02°，分辨率20 m，拍攝模式為升軌。

在形變監(jiān)測(cè)中，為了去除干涉相位中的地形相位，在處理過(guò)程中還需要高精度的DEM數(shù)據(jù)作為輔助數(shù)據(jù)進(jìn)行配準(zhǔn)、去除地形相位等操作。本次選用最新的AW3D DEM作為輔助數(shù)據(jù)。

InSAR數(shù)據(jù)處理中，在影像數(shù)據(jù)的圖像配準(zhǔn)、干涉圖生成、相位解纏、形變圖生成等步驟中均涉及到軌道信息。此次工作利用ENVI中的SARscape拓展模塊自動(dòng)讀取精密軌道數(shù)據(jù)，并添加到干涉處理過(guò)程中。

此次遙感監(jiān)測(cè)工作主要采用國(guó)產(chǎn)高分系列(GF-1、GF-2)衛(wèi)星數(shù)據(jù)，研究區(qū)涉及GF-2數(shù)據(jù)和少量GF1數(shù)據(jù)共19景，覆蓋解譯區(qū)全區(qū)。

3 結(jié)果分析

3.1 時(shí)空基線分析

將軌道Path/Frame為11/116的Sentinel-1數(shù)據(jù)源進(jìn)行格式轉(zhuǎn)化，采用時(shí)序處理方法中的小基線集技術(shù)進(jìn)行時(shí)間序列InSAR技術(shù)處理。時(shí)間基線閾值設(shè)置為90 d，空間基線閾值設(shè)置為5%，構(gòu)建干涉像對(duì)組合，用于后續(xù)的差分干涉處理。

為提高解纏的精確性，選取可用于3D解纏的連接圖。選取39景Sentinel-1數(shù)據(jù)共計(jì)134個(gè)干涉像對(duì)組合，干涉圖組合圖如圖3—圖5所示。

圖3 時(shí)空基線連接Fig.3 Space-time baseline connection

圖4 時(shí)間和位置的連接Fig.4 Connection of time and position

圖5 3D解纏連接Fig.5 3D untangled connection

3.2 差分干涉圖

對(duì)差分干涉圖進(jìn)行空間濾波和相位解纏，獲取差分干涉圖(圖6)。研究中，為避免基線誤差的影響，采用GCP點(diǎn)結(jié)合地面高程值和對(duì)應(yīng)的解纏干涉相位進(jìn)行基線優(yōu)化，優(yōu)化后再次進(jìn)行差分干涉、濾波和相位解纏，獲取解纏后的差分干涉相位。本次選取其中一景延安市子長(zhǎng)市瓦窯堡鎮(zhèn)(Path/frame：11/116)差分干涉圖和相干系數(shù)圖，如圖6所示(雷達(dá)坐標(biāo)系下)。

圖6 SAR干涉對(duì)查分結(jié)果Fig.6 SAR interference matching results

從差分干涉圖上能夠直接發(fā)現(xiàn)由形變引起的干涉條紋，見(jiàn)圖6中的虛線框位置，同時(shí)在相干系數(shù)圖(圖6 (c))可以看出這些形變區(qū)域明顯。為進(jìn)一步分析潛在的形變區(qū)域，將所有建立連接的影像對(duì)進(jìn)行差分干涉處理，選取干涉質(zhì)量較優(yōu)的差分干涉像對(duì)，高相干點(diǎn)識(shí)別方法選擇時(shí)序相干均值大于0.35相干目標(biāo)點(diǎn)進(jìn)行SBAS-InSAR反演，對(duì)去除差分干涉相位進(jìn)行時(shí)間序列分析，對(duì)殘差相位采用空間域?yàn)V波方法進(jìn)行大氣延遲相位提取，獲取大氣延遲相位后，再次進(jìn)行差分干涉得到差分干涉圖。此時(shí)仍然存在恒定相位和越變相位，進(jìn)行軌道精煉和重去平，消除上述相位引起的誤差，保證干涉圖質(zhì)量，提高SBAS-InSAR反演結(jié)果的精度。對(duì)比時(shí)序SAR像對(duì)獲取的相干系數(shù)和解纏結(jié)果，找到多數(shù)干涉圖中相干性高的區(qū)域并在此選擇GCP點(diǎn)，進(jìn)行重去平處理。選擇線性模型，通過(guò)Goldstein濾波消除大氣效應(yīng)引起的相位誤差，得到最終的時(shí)序形變量，并將得到的平均形變速率圖及時(shí)序形變量進(jìn)行地理編碼，轉(zhuǎn)換到地理坐標(biāo)系下。

3.3 SBAS-InSAR監(jiān)測(cè)結(jié)果分析

基于SARscape軟件，進(jìn)行了基線估算、配準(zhǔn)，Sentinel-1數(shù)據(jù)配準(zhǔn)精度均優(yōu)于0.001個(gè)像元，滿足干涉的要求；通過(guò)差分干涉、濾波、解纏等，利用39景Sentinel-1數(shù)據(jù)獲取了2018年1月—2020年2月陜北子長(zhǎng)地區(qū)的地表形變信息(沿雷達(dá)視線向，即LOS向)。監(jiān)測(cè)結(jié)果包括平均形變速率、形變區(qū)分布、累計(jì)形變量等信息(圖7、圖8、圖9)。利用視坡夾角將沿雷達(dá)視線方向的形變速率轉(zhuǎn)化成沿斜坡坡度方向的形變速率，以表達(dá)沿斜坡面的地表形變信息，形變速率圖和累計(jì)形變量圖中紅色代表遠(yuǎn)離衛(wèi)星視線方向，即代表沉降；藍(lán)色代表靠近衛(wèi)星視線方向，即代表抬升。

圖7 地表形變速率Fig.7 Surface deformation rate

圖8 地表累計(jì)形變Fig.8 Cumulative deformation diagramof surface deformation

圖9 形變區(qū)分布Fig.9 Distribution of deformation areas

由形變速率圖和累計(jì)形變圖可知，本次監(jiān)測(cè)子長(zhǎng)地表最大形變速率為-76 mm/a，累計(jì)最大形變量為-152 mm。形變區(qū)域主要分布在子長(zhǎng)市瓦窯堡鎮(zhèn)和余家坪鎮(zhèn)，其中子長(zhǎng)市南部的形變區(qū)域相較于北部較多，主要集中于南部的瓦窯堡鎮(zhèn)、余家坪鎮(zhèn)和楊家園子鎮(zhèn)；北部的澗峪岔鎮(zhèn)和南溝岔鎮(zhèn)的形變區(qū)域分布較少。

子長(zhǎng)市地表形變主要為煤礦開(kāi)采導(dǎo)致地表塌陷變形，其變形區(qū)形態(tài)在InSAR結(jié)果上較為規(guī)則，形態(tài)多呈長(zhǎng)條狀或類似矩形，與煤礦開(kāi)采區(qū)形態(tài)相關(guān)。

3.4 典型點(diǎn)分析

子長(zhǎng)市新磘上滑坡(ZC1),累積形變?nèi)鐖D10所示。

圖10 滑坡點(diǎn)ZC1的累積形變(沿LOS向形變)Fig.10 Cumulative deformation of ZC1 at landslidepoint (deformation along LOS direction)

從圖10中可以看出，滑坡點(diǎn)ZC1有明顯的形變信號(hào)，結(jié)合光學(xué)影像和地形分析，此處滑坡在雷達(dá)視線的背坡向，其從高向低滑動(dòng)，在InSAR監(jiān)測(cè)結(jié)果中表現(xiàn)出負(fù)的形變(遠(yuǎn)離衛(wèi)星)。滑坡區(qū)域中某特征點(diǎn)的形變曲線顯示(圖11)，此處滑坡在2018年1月—2019年11月間不同時(shí)刻有明顯的形變趨勢(shì)。

圖11 ZC1滑坡區(qū)中某特征點(diǎn)的時(shí)序形變曲線Fig.1 Time-series deformation curve of a characteristic point in ZC1 landslide area

4 結(jié)論

基于SBAS-InSAR技術(shù)開(kāi)展陜北子長(zhǎng)地區(qū)地表形變InSAR監(jiān)測(cè)，獲取了2018年1月—2020年2月子長(zhǎng)地區(qū)的地表沉降信息(沿雷達(dá)視線向，即LOS向)。子長(zhǎng)地表最大形變速率為-76 mm/a，累計(jì)最大形變量為-152 mm。形變區(qū)域主要分布在子長(zhǎng)市瓦窯堡鎮(zhèn)和余家坪鎮(zhèn)。

通過(guò)基于InSAR地表形變監(jiān)測(cè)、綜合遙感技術(shù)監(jiān)測(cè)、野外調(diào)查等工作，排除偽隱患點(diǎn)，形成了對(duì)區(qū)域地質(zhì)災(zāi)害發(fā)育及嚴(yán)重程度的正確認(rèn)識(shí)。初步總結(jié)了一套InSAR形變分類方法，由于地表形變的形成因素較多以及InSAR數(shù)據(jù)獲取、處理過(guò)程中的誤差影響，InSAR監(jiān)測(cè)獲取的形變到地質(zhì)災(zāi)害隱患確定過(guò)程尤為重要，對(duì)形變點(diǎn)(區(qū))進(jìn)行分類，有益于地質(zhì)災(zāi)害隱患點(diǎn)的獲取，為隨后數(shù)據(jù)的處理優(yōu)化也提供了借鑒。初步總結(jié)了一套適用于子長(zhǎng)地區(qū)地質(zhì)災(zāi)害InSAR監(jiān)測(cè)技術(shù)體系，本次基于InSAR技術(shù)開(kāi)展區(qū)域性地質(zhì)災(zāi)害監(jiān)測(cè)預(yù)警，綜合應(yīng)用高分辨率多光譜光學(xué)影像遙感、合成孔徑雷達(dá)干涉測(cè)量、三維傾斜攝影測(cè)量等對(duì)地觀測(cè)技術(shù)，獲取區(qū)域性地質(zhì)災(zāi)害的分布特征、變形狀況，結(jié)合實(shí)地調(diào)查核實(shí)工作，圈定并監(jiān)測(cè)重點(diǎn)地表形變區(qū)，確定地質(zhì)災(zāi)害隱患點(diǎn)，編制地質(zhì)災(zāi)害隱患分布圖。