殷宗敏，韓守富，楊 環(huán)，馬金輝

(蘭州大學(xué)資源環(huán)境學(xué)院，甘肅 蘭州 730000)

綏德縣位于陜北黃土高原的東北部，區(qū)內(nèi)總的地勢為：北高南低，中東部無定河與黃河的分水嶺高于兩側(cè)，無定河兩側(cè)梁峁地勢高于河谷谷地，黃河谷地最低?？h域位于鄂爾多斯地臺東部，構(gòu)造活動比較弱，區(qū)內(nèi)僅發(fā)育少量平緩開闊的褶皺。隨著人類改造自然能力的增強(qiáng)，人類活動對自然環(huán)境的改變越來越大[1]。綏德縣與地質(zhì)災(zāi)害關(guān)系密切的人類工程活動主要有交通建設(shè)、水利設(shè)施建設(shè)、群眾修建房屋時(shí)削坡過陡等，這些建設(shè)生產(chǎn)過程導(dǎo)致滑坡情況時(shí)有發(fā)生，造成人員傷亡、財(cái)產(chǎn)損失。

圖1 綏德縣城周邊發(fā)育滑坡

滑坡是區(qū)內(nèi)最為嚴(yán)重的地質(zhì)災(zāi)害之一。境內(nèi)滑坡主要發(fā)育在河谷兩側(cè)邊坡上和黃土梁峁區(qū)的陡坡地段，穩(wěn)定性差，大多為黃土滑坡，如圖1所示。其發(fā)生一般預(yù)兆不明顯，運(yùn)動速度快，發(fā)生過程短促，一旦發(fā)生，往往造成人員財(cái)產(chǎn)損失或者影響交通，因此對該區(qū)域進(jìn)行地表形變監(jiān)測尤為重要。利用時(shí)序InSAR技術(shù)進(jìn)行地表形變監(jiān)測，提取研究區(qū)長時(shí)間序列地表變形數(shù)據(jù)[2]，基于地表變形速率圈定形變速率較大的不穩(wěn)定區(qū)域，該區(qū)域可能為滑坡區(qū)。通過實(shí)地調(diào)查等方式，對這些區(qū)域進(jìn)行甄別，以便進(jìn)行相應(yīng)的防災(zāi)減災(zāi)工作。由于城區(qū)人口聚集且存在發(fā)生滑坡的條件，因此對于該區(qū)域的監(jiān)測尤為重要，故選擇的研究范圍如圖2所示。

圖2 研究區(qū)范圍

1 數(shù)據(jù)來源及研究范圍

1.1 數(shù)據(jù)基本信息

本文選擇Sentinel-1A IW SLC數(shù)據(jù)，軌道為11，下載網(wǎng)址為：https://scihub.copernicus.eu/。影像截止日期為2017年9月20日，覆蓋研究區(qū)的影像共38景，數(shù)據(jù)基本特征如表1所示。

Sentinel-1A數(shù)據(jù)大的空間基線、多普勒中心和數(shù)據(jù)的時(shí)間分布情況如圖3所示，其中空間基線以0為中心，正態(tài)分布于兩側(cè)，值在150 m以內(nèi)，這為圖像高精度配準(zhǔn)提供了保障，多普勒中心分布于0～-0.1之間，影像的獲取時(shí)間以2016年12月30日為中心，且中心時(shí)間點(diǎn)以后的影像更為密集。

1.2 數(shù)據(jù)時(shí)間分布

Sentinel-1A的11軌道為升軌道數(shù)據(jù)，以2016年12月30日影像為主影像，影像的時(shí)間序列分布情況相對均勻，其時(shí)間序列分布以及相應(yīng)的空間基線、時(shí)間基線、多普勒中心的分布情況如表2所示。

表1 所用數(shù)據(jù)基本信息

表2 所用Sentinel-1A IW SLC數(shù)據(jù)時(shí)間分布表

2 數(shù)據(jù)處理方法

2.1 時(shí)序Insar參數(shù)優(yōu)化和處理

研究區(qū)地表覆蓋變化大，失相干現(xiàn)象嚴(yán)重，難以捕獲到永久散射體[3]。滑坡變形既包括緩慢蠕動，又包括快速變形，不僅僅存在線性變形，更多的是非線性變形。針對InSAR處理過程中的主要環(huán)節(jié)和目前所面臨的問題，為提高滑坡區(qū)InSAR測量精度，選擇合適的圖像連接拓?fù)洹V波、解纏方案及參數(shù)；進(jìn)行大氣延遲效應(yīng)估算，獲得適合研究區(qū)地形和大氣特征的大氣延遲效應(yīng)參數(shù)[4]。

2.2 圖像連接拓?fù)?、濾波

通過不同方案和參數(shù)的設(shè)定，提高滑坡區(qū)InSAR測量精度[5]，擴(kuò)大InSAR測量范圍，針對特定的SAR數(shù)據(jù)集，建立適當(dāng)?shù)膱D像拓?fù)洌纾夯趩沃鲌D像的星型結(jié)構(gòu)、基于空間基線最佳的多主圖像拓?fù)浣Y(jié)構(gòu)、基于時(shí)間基線的多主圖像拓?fù)浣Y(jié)構(gòu)，分析其棄除大氣效應(yīng)的不同效果從而選擇最佳圖像拓?fù)洌疚淖罱K選擇星型結(jié)構(gòu)；選擇適合研究區(qū)地表覆蓋、變形規(guī)律的濾波方式，如：Boxcar Goodstein和Adaptive Filter等，并進(jìn)行參數(shù)調(diào)整，分析其對干涉效果和最終變形結(jié)果的影響，本文最終選擇Adaptive Filter濾波方法。

2.3 剔除或減輕地形影響

選擇了不同來源和分辨率的DEM數(shù)據(jù)[6]，測試其分離地形效應(yīng)的效果，獲得能有效分離地形效應(yīng)的DEM精度大小和分辨率等指標(biāo)，這里選擇90米分辨率dem；采用多角度觀測來減緩因地形產(chǎn)生的SAR圖像幾何畸變(透視收縮、陰影、疊掩)，增加SAR測量點(diǎn)密度。

圖3 軌道數(shù)據(jù)集基本參數(shù)圖

2.4 大氣延遲效應(yīng)估算和剔除

Processing Parameters分區(qū)是SARProZ軟件反演大氣相位的主要功能模塊，其通過干涉計(jì)算獲得各像素點(diǎn)的干涉相位，并依據(jù)研究區(qū)域的地表形變情況以及干涉相位中各分量的物理及統(tǒng)計(jì)特性建立其與高程誤差、線性形變、熱形變等因子參數(shù)的纏繞相位函數(shù)模型，根據(jù)大氣相位分析像素點(diǎn)網(wǎng)鄰近永久散射體目標(biāo)(像素點(diǎn))間的相位差分關(guān)系，采用相關(guān)算法考察相鄰像素點(diǎn)差分干涉相位關(guān)系對此函數(shù)模型進(jìn)行求解，解算出各像素點(diǎn)的高程誤差、線性形變等貢獻(xiàn)因子相位，并對殘余相位中的非線性形變與失相干噪聲相位進(jìn)行分解，最終計(jì)算獲得精確的大氣延遲分量纏繞相位值。大氣延遲相位是InSAR測量的主要干擾誤差源[7]，需要對該誤差進(jìn)行剔除，其處理步驟為：通過調(diào)整參數(shù)，有效的選取PS候選點(diǎn)；對選取的PS候選點(diǎn)構(gòu)建合理的分析像素點(diǎn)網(wǎng)；結(jié)合研究區(qū)域的現(xiàn)狀及數(shù)據(jù)集的情況，給出合理的參數(shù)，分析大氣延遲相位[8]；選擇合理的參考點(diǎn)并檢驗(yàn)各參數(shù)相位分量的合理性。

圖4 研究區(qū)形變速度圖

3 結(jié)果與分析

3.1 實(shí)驗(yàn)結(jié)果

將覆蓋研究區(qū)的不同時(shí)段的多期SAR影像按照成像時(shí)間先后順序，以選取的公共主影像為基準(zhǔn)，將其余所有影像都配準(zhǔn)并取樣到主影像空間[9]，總共可形成K個(gè)干涉對。通過相位差分處理，可得到K幅干涉相位圖。然后借助提取影像中永久散射體(PS)目標(biāo)，通過分離地形數(shù)據(jù)誤差、大氣延遲誤差等一系列的誤差項(xiàng)，從而提高PS的地表形變估計(jì)精度，最終得到可靠的地表形變結(jié)果。將獲取的可靠的時(shí)間序列形變結(jié)果導(dǎo)出為csv和kml文件，將kml結(jié)果在谷歌地球中顯示，其結(jié)果如圖4所示，其中形變速度值范圍為-15mm/year至15mm/year,顏色越深表示形變速度越大，接近-15mm/year為紅色，接近15mm/year為藍(lán)色。

結(jié)合地形指標(biāo)、ps點(diǎn)的形變速度值和聚集情況，勾畫出可能的滑坡區(qū)。發(fā)現(xiàn)其多分布于山坡部分且一側(cè)為道路或者河流谷底(見圖4)，符合自然條件下滑坡的發(fā)生區(qū)域范疇，說明基于時(shí)間序列的PS點(diǎn)能夠識別滑坡。

3.2 結(jié)果分析

視線方向的形變結(jié)果在空間呈現(xiàn)聚集程度不一的點(diǎn)要素，其中空間分布呈現(xiàn)聚集狀態(tài)的點(diǎn)要素，其聚集區(qū)域往往代表地表變化大的區(qū)域，在人工活動區(qū)域多表現(xiàn)為建筑開挖，而非人工影響區(qū)域則多表現(xiàn)為滑坡，在谷歌地圖中對點(diǎn)聚集區(qū)域沿著山坡繪制可能滑坡區(qū)，并進(jìn)行滑坡的野外驗(yàn)證。對滑坡前緣進(jìn)行調(diào)察時(shí)，對比滑動面、滑坡前緣邊界，需要注意的是，滑坡前緣一般存在剪出口(剪切破壞明顯凸起甚至反翹的部位)，如果滑坡是由于坡腳開挖引起，前緣邊界應(yīng)在坡腳填土前緣。

野外驗(yàn)證發(fā)現(xiàn)：先前圈定的滑坡區(qū)域，在現(xiàn)場基本都找到了滑動痕跡(見圖5)。其中，圖片1現(xiàn)場為施工區(qū)域，坡面呈現(xiàn)滑動跡象，且前方道路塌陷，表示該區(qū)域地表發(fā)生了較大形變，與Insar計(jì)算的結(jié)果吻合；圖片2為一個(gè)滑坡體前的窯洞內(nèi)部情況，由于山體滑動，向前產(chǎn)生推力，導(dǎo)致墻體夯土向窯洞內(nèi)側(cè)脫落；圖片3為無定河一側(cè)山谷中山體滑動后形成的階梯狀痕跡，其下窯洞存在較大的安全風(fēng)險(xiǎn)；圖片4為公路一側(cè)的山體，其存在明顯的滑坡后壁，滑動痕跡明顯。通過實(shí)地的考察對比分析，證明利用時(shí)序Insar技術(shù)進(jìn)行地表形變監(jiān)測，對形變速度較大且形變點(diǎn)聚集的區(qū)域進(jìn)行圈定，能夠找出滑坡區(qū)域，指導(dǎo)防災(zāi)減災(zāi)工作。

圖5 野外滑坡驗(yàn)證對比圖

在眾多滑坡體中，選定3個(gè)明顯的單體滑坡(單體滑坡輪廓見圖6)，獲取該區(qū)域的形變曲線，進(jìn)行滑坡單體變形研究。定位滑坡體上的形變點(diǎn)，并在時(shí)序Insar處理過程中獲得的csv文件里面找出該點(diǎn)對應(yīng)的各時(shí)間段形變值，將其繪制成折線圖(見圖7)。

圖6 滑坡單體輪廓圖

圖7 坡體變形曲線圖

從A坡變形曲線可以看到，坡體A在2016年8月份左右發(fā)生了較大規(guī)模的滑動，之后呈現(xiàn)較緩的滑動，由于該折線整體呈現(xiàn)下降趨勢，說明滑動在持續(xù)。

B坡變形曲線中，發(fā)生較大規(guī)?；瑒拥臅r(shí)間點(diǎn)為2016年9月份前后，與A滑坡體類似，形變曲線呈現(xiàn)下降趨勢，形變曲線兩區(qū)間段的形變量差值大于A滑坡體，說明其不穩(wěn)定性比A滑坡體強(qiáng)。

C坡變形曲線與A、B有所不同，其在2016年10月形變量開始趨于穩(wěn)定，在2017年7月之前，形變速度基本為正值，在7月份以后呈現(xiàn)較大變化，說明較大滑動發(fā)生在2017年7月，C滑坡體形變量相對小于A、B，說明C滑坡體較A、B穩(wěn)定。

綜上可以看出，通過時(shí)間序列Insar獲取的該區(qū)域形變點(diǎn)屬性值能夠反映出滑坡體的變形過程，說明Insar技術(shù)能夠很好的識別出滑坡體。

4 討論與結(jié)論

4.1 討論

基于時(shí)間序列Insar技術(shù)，能夠獲得緩慢的地表形變過程，其在平原地區(qū)具有較好的效果[10]，而在山區(qū)由于植被的覆蓋以及山體對雷達(dá)信號的遮擋，降低了Ps-Insar技術(shù)的可靠性，但是可以通過參數(shù)的調(diào)整，逐步剔除部分不準(zhǔn)確相位，提高結(jié)果的可靠性。此外，參數(shù)的選擇具有極強(qiáng)的地域性以及技術(shù)性，需要通過不斷的試驗(yàn)來獲取最佳的參數(shù)，而試驗(yàn)又存在主觀性，因此只能不斷接近最佳的參數(shù)集。

由于衛(wèi)星獲取的形變方向?yàn)橐暰€方向[11]，衛(wèi)星的軌道固定，同時(shí)山體不同部分的坡度坡向不一致，導(dǎo)致獲取的形變結(jié)果與坡體形變的真實(shí)情況存在一定的誤差，而目前還沒有很好的方法解決這些偏差，有待進(jìn)一步考究。

4.2 結(jié)論

基于PS-InSAR技術(shù)，利用Sentinel影像數(shù)據(jù)，通過影像裁剪、影像配準(zhǔn)、地形相位估計(jì)、差分干涉、高相干點(diǎn)選取、相位解纏、變形提取和誤差估計(jì)等一系列數(shù)據(jù)處理手段計(jì)算的黃土高原重要城鎮(zhèn)時(shí)間序列地表變形速率具有很高的可靠性，能夠很好的識別黃土高原區(qū)滑坡區(qū)，滑坡體上形變點(diǎn)的形變曲線能較好的反映出坡體的變形情況。

[1]黃潤秋，許強(qiáng). 中國典型災(zāi)難性滑坡[M]. 北京：科學(xué)出版社.2008.

[2]Michele Crosetto , Oriol Monserrat , María Cuevas-González , Núria Devanthéry, Bruno Crippa ．Persistent Scatterer Interferometry[J]．ISPRS Journal of Photogrammetry and Remote Sensing，2015.

[3]Dong S, Samsonov S, Yin H, et al. Time-series analysis of subsidence associated with rapid urbanization in Shanghai, China measured with SBAS InSAR method[J]. Environmental Earth Sciences, 2014.72(3):677-691.

[4]Schl gel R, Doubre C, Malet J P, et al. Landslide deformation monitoring with ALOS/PALSAR imagery: A D- InSAR geomorphological interpretation method.[J]. Geomorphology.2014.231: 314-330.

[5]陳強(qiáng)．基于永久散射體雷達(dá)差分干涉探測區(qū)域地表形變的研究[D]．成都：西南交通大學(xué)博士論文.2006.

[6]傅文學(xué)，田慶久，郭小方，等．技術(shù)及其在地表形變監(jiān)測中的應(yīng)用現(xiàn)狀與發(fā)展[J]．地球科學(xué)進(jìn)展.2006.

[7]陳強(qiáng)，李永樹，劉國祥．干涉雷達(dá)永久散射體識別方法的對比分析[J].遙感信息.2006.

[8]劉國祥，陳強(qiáng)，丁曉利．基于雷達(dá)干涉永久散射體網(wǎng)絡(luò)探測地表形變的算法與實(shí)驗(yàn)結(jié)果[J]．測繪學(xué)報(bào).2007.

[9]汪魯才．星載合成孔徑雷達(dá)干涉成像的信息處理方法研究[D]．湖南大學(xué)博士論文.2005.11.

[10]李德仁，縻明生，王艷．永久散射體雷達(dá)干涉測量技術(shù)[J]．武漢大學(xué)學(xué)報(bào)(信息科學(xué)版).2004.

[11]彭曙蓉．高分辨率合成孔徑雷達(dá)干涉測量技術(shù)及其應(yīng)用研究[D]．湖南大學(xué).2009.