孫 鵬

（新疆水利水電勘測(cè)設(shè)計(jì)研究院有限責(zé)任公司，新疆烏魯木齊 830000）

0 引言

合成孔徑雷達(dá)干涉測(cè)量技術(shù)（Interferometric Synthetic Aperture Radar，InSAR）因其具有高精度、高分辨率、全天候等優(yōu)點(diǎn)，是監(jiān)測(cè)地表形變中一種方法。該技術(shù)應(yīng)用在城市地表沉降監(jiān)測(cè)、礦區(qū)地表沉降監(jiān)測(cè)、山體滑坡監(jiān)測(cè)等項(xiàng)目中，發(fā)揮了極大的作用。InSAR 技術(shù)演化出Stacking-InSAR 技術(shù)與SBAS-InSAR 技術(shù)，并可將二者結(jié)合起來作為一種組合方法應(yīng)用在地表沉降、滑坡監(jiān)測(cè)中［1-2］。

本文利用InSAR 技術(shù)對(duì)新疆阿勒泰地區(qū)布爾津滑坡進(jìn)行監(jiān)測(cè)與分析。Sentinel-1 系統(tǒng)在周邊區(qū)域擁有豐富的InSAR 數(shù)據(jù)積累、高重訪周期與大面積的影像覆蓋范圍，收集基于2014 年10 月20 日至2021 年2 月9 日采集的Sentinel-1 升軌影像，數(shù)據(jù)處理以增強(qiáng)型SBAS-InSAR 集成時(shí)序分析為主要技術(shù)手段，獲取布爾津滑坡區(qū)域形變的趨勢(shì)、范圍和梯度等信息，并從中篩選出地表形變變化量達(dá)到預(yù)警臨界值，和發(fā)生形變突變的可疑區(qū)域，作為地表形變監(jiān)測(cè)預(yù)警的重點(diǎn)工作區(qū)域。

1 基本原理與技術(shù)優(yōu)勢(shì)

1.1 Stacking-InSAR技術(shù)

干涉圖疊加（Stacking-InSAR）技術(shù)的數(shù)據(jù)處理基礎(chǔ)是進(jìn)行經(jīng)典二軌差分InSAR 技術(shù)處理。Stacking-InSAR 技術(shù)是在獲取M 幅相位解纏的差分干涉圖后，將解纏相位加權(quán)疊加，估計(jì)出平均相位變化速率。由于大氣相位分量在空間域具有低頻特性、時(shí)間域具有高頻特性，故Stacking 技術(shù)最大程度地減少了大氣誤差，提高了形變速率解算精度［3-4］。該方法假設(shè)每一幅獨(dú)立的解纏相位圖中，在不考慮其他噪聲的分布情況下，區(qū)域的形變速率為總形變量與時(shí)間的比值，即線性形變速率，大氣延遲誤差隨機(jī)且相等。通常情況下，為保證結(jié)果的準(zhǔn)確性，選擇時(shí)間和空間基線都比較短的數(shù)據(jù)進(jìn)行組合。

1.2 SBAS-InSAR技術(shù)

時(shí)序分析采用小基線集干涉測(cè)量技術(shù)（Small Baseline Subset，SBAS）充分考慮大量的干涉圖信息，能有效地消除和削弱解纏粗差、大氣誤差以及DEM誤差等因素的影響，從而獲取更高精度的地表形變成果［5］。基于SBAS 技術(shù)的數(shù)據(jù)處理主要解算步驟如下：

1）對(duì)N 幅SAR 影像數(shù)據(jù)按一定時(shí)空基線條件進(jìn)行干涉組合處理，形成M 幅干涉圖，利用已有DEM（本次采用SRTM）作為外部高程數(shù)據(jù)，進(jìn)行差分處理，生成差分干涉圖，以去除地形及平地效應(yīng)影響；

2）對(duì)差分干涉圖進(jìn)行自適應(yīng)濾波處理以去除相位噪聲影響，對(duì)濾波后的差分干涉圖進(jìn)行相位解纏，生成解纏相位圖；

3）根據(jù)基線條件和干涉相位信息估算高程誤差及線性形變相位，在原始干涉相位中減去估算的線性形變相位得到殘余相位，此時(shí)的殘余相位中主要包括非線性形變及大氣相位；

4）解纏此殘余相位，并補(bǔ)償線性形變相位部分得到完整的形變相位，此時(shí)相位中還包含有大氣相位的影響；

5）對(duì)去除線性形變的殘余相位進(jìn)行空域低通和時(shí)域高通濾波處理以分離出大氣相位；

6）在形變相位中減去大氣相位影響，得到形變相位值；

7）基于SVD 的形變求解；

8）對(duì)形變進(jìn)行地理編碼，獲取WGS-84 坐標(biāo)系下的形變成果圖。

1.3 總體技術(shù)路線

本文根據(jù)監(jiān)測(cè)區(qū)域基本情況，計(jì)劃出合理可行的技術(shù)路線。應(yīng)用雷達(dá)監(jiān)測(cè)數(shù)據(jù)，采用Stacking-InSAR 技術(shù)與SBAS-InSAR 技術(shù)相結(jié)合方式實(shí)現(xiàn)了變形數(shù)據(jù)提取。根據(jù)地質(zhì)變形情況，及時(shí)建立災(zāi)害預(yù)警與防治措施。本文總體技術(shù)路線如圖1 所示。

圖1 Stacking-InSAR 和SBAS-InSAR 時(shí)序分析提取地表形變技術(shù)路線圖

2 研究區(qū)域與數(shù)據(jù)源

2.1 研究區(qū)域

布爾津滑坡位于新疆阿勒泰地區(qū)布爾津縣境內(nèi)，滑坡位于水電站庫(kù)區(qū)的上游區(qū)域。該區(qū)域山體受雨水的沖刷導(dǎo)致地質(zhì)情況不穩(wěn)定，很容易出現(xiàn)山體滑坡現(xiàn)象，造成河流堵塞帶來的發(fā)電站隱患問題。研究區(qū)域范圍內(nèi)無居住區(qū)，地勢(shì)起伏較大；最低點(diǎn)海拔+203.4 m，最高點(diǎn)海拔+1125.9 m，平均海拔+841.9 m。該滑坡的研究范圍如圖2 所示。

圖2 研究區(qū)域范圍

為比較本文提及的Stacking-InSAR 和SBASInSAR 技術(shù)在滑坡體監(jiān)測(cè)中的有效性，在滑坡區(qū)域布設(shè)監(jiān)測(cè)點(diǎn)，并采用傳統(tǒng)極坐標(biāo)法獲取監(jiān)測(cè)點(diǎn)位移。為便于體現(xiàn)各監(jiān)測(cè)點(diǎn)位移方向，建立監(jiān)測(cè)坐標(biāo)系，X方向?yàn)槠叫泻拥?，Y方向?yàn)榇怪焙拥?，通過坐標(biāo)變化量直觀表示各點(diǎn)位移趨勢(shì)，監(jiān)測(cè)基準(zhǔn)網(wǎng)及監(jiān)測(cè)點(diǎn)布置如圖3 所示。

圖3 監(jiān)測(cè)基準(zhǔn)網(wǎng)及監(jiān)測(cè)點(diǎn)布置圖

2.2 數(shù)據(jù)源

2.2.1 影像數(shù)據(jù)

為避免SAR 影像幾何畸變的影響，所用SAR數(shù)據(jù)為升軌Sentinel-1 數(shù)據(jù)，采用TOPS（Terrain Observation With Progressive Scans）掃描方式，共計(jì)153 景，覆蓋時(shí)間段為2014 年10 月20 日至2021年2 月9 日。數(shù)據(jù)基本參數(shù)如表1 所示，影像采集范圍如圖4 所示。

表1 升軌SAR 數(shù)據(jù)的基本參數(shù)

圖4 研究范圍數(shù)據(jù)覆蓋圖

2.2.2 DEM獲取

形變監(jiān)測(cè)過程中，數(shù)字高程模型（DEM）是非常重要的外部數(shù)據(jù)，在生成差分干涉圖時(shí)需要DEM 來模擬地形相位，此外在雷達(dá)坐標(biāo)系和地理坐標(biāo)系相互轉(zhuǎn)化（地理編碼）的過程中也需要高精度的DEM。

美國(guó)太空總署（NASA）和國(guó)防部國(guó)家測(cè)繪局（NIMA）于2000 年聯(lián)合發(fā)射的“奮進(jìn)號(hào)”航天飛機(jī)完成了為期11 天的航天飛機(jī)雷達(dá)制圖任務(wù)（Shuttle Radar Topography Mission，SRTM）。經(jīng)過兩年的數(shù)據(jù)處理，最終獲得了平面精度±20 m，高程精度±16 m 的全球數(shù)字高程模型，分辨率包含1 弧秒和3 弧秒兩種，分別對(duì)應(yīng)30 m 和90 m。本項(xiàng)目選用30 m分辨率的DEM 數(shù)據(jù)用于地形相位的消除。

2.2.3 極坐標(biāo)差分法獲取數(shù)據(jù)

依據(jù)153 景的SAR 影像的成像時(shí)間，選取2014 年10 月20 日至2021 年2 月9 日之間的SAR數(shù)據(jù)，每次監(jiān)測(cè)方法保持一致，采用相同觀測(cè)人員、儀器使用極坐標(biāo)差分法觀測(cè)。結(jié)合每次觀測(cè)成果，分析得出各監(jiān)測(cè)點(diǎn)觀測(cè)成果相對(duì)首期成果水平位移量和垂直位移量。位移量的X方向?yàn)檎硎狙睾酉蛏嫌畏较颍粗睾酉蛳掠畏较?，Y方向?yàn)檎硎敬怪焙拥老蚝拥婪较?，反之垂直河道向山頂方向；H方向?yàn)? 表示無變化，為負(fù)表示下沉。為了驗(yàn)證Stacking-InSAR 與SBAS-InSAR 對(duì)沉降邊緣與沉降中心的監(jiān)測(cè)能力，選取了圖2 中的6 個(gè)特征點(diǎn)，如表2 所示。其中，1-1、1-3、2-1 和2-2 屬于沉降邊緣區(qū)域，分布于滑坡區(qū)域西北側(cè)；1-3 和1-4屬于沉降中心區(qū)域，2-3 和2-4 位于滑坡區(qū)東南側(cè)。

表2 極坐標(biāo)法獲取的變形點(diǎn)坐標(biāo)值

3 數(shù)據(jù)處理與結(jié)果分析

對(duì)獲取的153 景Sentinel-1 數(shù)據(jù)進(jìn)行雷達(dá)差分干涉測(cè)量處理，按照時(shí)間序列構(gòu)建152 個(gè)干涉對(duì)，其中2020.08.15-2020.08.29 的垂直基線最長(zhǎng)，達(dá)到了200.89 m，2016.04.15-2016.04.20 的垂直基線最短，數(shù)值為43.18 m。按照?qǐng)D1 所示的Stacking-InSAR 和SBAS-InSAR 技術(shù)處理流程對(duì)152 個(gè)干涉對(duì)處理，獲取到相鄰影像之間的形變信息，通過時(shí)序累加，得到各時(shí)間節(jié)點(diǎn)上的累計(jì)沉降情況。

依據(jù)干涉圖的相干性最大原則，選取2020 年8月20 日的影像為主影像，其余影像與其配準(zhǔn)，分別設(shè)定時(shí)間基線閾值為100 d，空間基線閾值為220 m。將所有滿足閾值的限定的兩景SAR 影像兩兩組合，形成了152 個(gè)時(shí)序差分干涉對(duì)進(jìn)行處理。按照本文的技術(shù)路線圖，采用Stacking-InSAR 與SBAS-InSAR 技術(shù)處理152 個(gè)干涉對(duì)，得到研究區(qū)域形變信息。布爾津滑坡153 景SAR 數(shù)據(jù)計(jì)算獲得的2014 年10 月至2021 年2 月地表形變累計(jì)沉降量如圖5 所示。布爾津滑坡2014 年10 月至2021年2 月地表形變速率如圖6 所示。從這兩幅圖可以看出：布爾津滑坡區(qū)域這段時(shí)期內(nèi)最大累計(jì)形變量達(dá)-550 mm，最大形變速率達(dá)-76.9 mm/a。圖6中，負(fù)值表示地表正在遠(yuǎn)離衛(wèi)星，正值表示地表正在靠近衛(wèi)星。從圖6 可以看出：新疆阿勒泰布爾津滑坡區(qū)域分布有三處明顯變形體，變形體1 平均形變-6.0 mm/a，最大形變-25.4 mm/a；變形體2 平均形變-4.8 mm/a，最大形變-17.1 mm/a；變形體3 平均形變-19.0 mm/a，最大形變-76.9 mm/a。

圖5 地表累計(jì)形變圖

圖6 2014-2021 平均形變速率圖

為了更直觀感受研究區(qū)域內(nèi)的時(shí)空演變規(guī)律，這里展示SBAS-InSAR 的時(shí)序結(jié)果。2014 年10月、12 月與2021 年1 月、2 月時(shí)空演變情況如圖7-圖10 所示。

圖7 2014 年10 月時(shí)空演變結(jié)果

圖8 2014 年12 月時(shí)空演變結(jié)果

圖9 2021 年1 月時(shí)空演變結(jié)果

圖10 2021 年2 月時(shí)空演變結(jié)果

為進(jìn)一步研究該滑坡形變的時(shí)間演化，在3 個(gè)形變體上分別提取3 個(gè)特征點(diǎn)，獲得滑坡上這些點(diǎn)2014 年10 月至2021 年2 月的時(shí)間序列形變，如圖11-圖16 所示。變形體1 最大的累積形變出現(xiàn)在其P3 點(diǎn)上，這個(gè)時(shí)間段內(nèi)累積形變達(dá)到了-174mm；變形體2 最大的累積形變出現(xiàn)在其P2 點(diǎn)上，這個(gè)時(shí)間段內(nèi)，累積形變達(dá)到了-223 mm；變形體3 最大的累積形變出現(xiàn)在其P3 點(diǎn)上，這個(gè)時(shí)間段內(nèi)，累積形變達(dá)到了-499 mm。3 個(gè)變形體在2014 年10 月至2021 年2 月存在明顯的非線性形變趨勢(shì)?？梢钥吹剑?018 年1 月至2021 年2月，3 個(gè)變形體均出現(xiàn)快速形變現(xiàn)象，形變量比較大，其中變形體3 從2016 年年初就開始出現(xiàn)加速的情況。

圖11 變形體1 特征點(diǎn)示意圖

圖12 變形體1 地表變形速率及特征點(diǎn)累計(jì)變形量

圖13 變形體2 特征點(diǎn)示意圖

圖14 變形體2 地表形變速率及特征點(diǎn)累計(jì)變形量

圖15 變形體3 特征點(diǎn)示意圖

圖16 變形體3 地表形變速率及特征點(diǎn)累計(jì)變形量

4 精度驗(yàn)證

4.1 傳統(tǒng)測(cè)量?jī)x器結(jié)果比對(duì)

為了驗(yàn)證Stacking-InSAR 和SBAS-InSAR 技術(shù)在滑坡監(jiān)測(cè)方面的精度，采用傳統(tǒng)測(cè)繪方法在滑坡體上布置一定的測(cè)點(diǎn)?？紤]到滑坡體不同位置移動(dòng)情況不同，在沉降邊緣區(qū)域選擇1-1 號(hào)點(diǎn)位，該點(diǎn)位于滑坡體西北側(cè)；沉降中心區(qū)域選擇1-3、1-4 號(hào)點(diǎn)位；2-4 號(hào)點(diǎn)位位于滑坡區(qū)東南側(cè)。各點(diǎn)位分別比較了水平方向與垂直方向移動(dòng)量。以傳統(tǒng)測(cè)繪儀器獲取的數(shù)值為真值，為了減小位置差錯(cuò)產(chǎn)生的誤差，以監(jiān)測(cè)點(diǎn)為中心，將一定范圍內(nèi)的SBAS點(diǎn)求取平均值作為SBAS 監(jiān)測(cè)值。表3、表4 統(tǒng)計(jì)了在2018 年1 月至2021 年2 月時(shí)間段內(nèi)二者比較結(jié)果。

表3 水平方向移動(dòng)量比較

表4 垂直方向移動(dòng)量比較

由表3、表4 可知：不管是用InSAR 技術(shù)還是傳統(tǒng)測(cè)繪儀器技術(shù)，變形體位移下降趨勢(shì)相同，位移量相近、數(shù)據(jù)相差較小?？梢哉J(rèn)為采用Stacking-InSAR 和SBAS-InSAR 技術(shù)的方式監(jiān)測(cè)地表形變能達(dá)到亞厘米級(jí)的精度。產(chǎn)生誤差的主要原因是兩種方法監(jiān)測(cè)點(diǎn)位置存在偏差，不能完全對(duì)應(yīng)。采用InSAR 技術(shù)監(jiān)測(cè)地面沉降，精度能夠達(dá)到厘米級(jí)甚至毫米級(jí)別，能夠滿足監(jiān)測(cè)項(xiàng)目需求。

5 結(jié)語

本文利用2014 年10 月至2021 年2 月期間覆蓋新疆阿勒泰地區(qū)的多景Sentinel-1 數(shù)據(jù)，應(yīng)用Stacking-InSAR 與SBAS-InSAR 技術(shù)相結(jié)合的方式獲取了該地區(qū)該時(shí)間段的地表累計(jì)形變量圖，分析得知該地區(qū)最大累計(jì)形變量達(dá)-550 mm，最大形變速率達(dá)-76.9 mm/a，可見最大變形體一直處于持續(xù)的形變狀態(tài)。提出以下建議：

1）實(shí)地獲取下游水電站相關(guān)數(shù)據(jù)，對(duì)該滑坡地質(zhì)情況進(jìn)行調(diào)查，做好防護(hù)工作與監(jiān)測(cè)風(fēng)險(xiǎn)的情況的工作。

2）為保證變形體能夠得到實(shí)時(shí)監(jiān)測(cè)，可結(jié)合自動(dòng)化監(jiān)測(cè)技術(shù)，保證變形體全方位、多角度、可視化監(jiān)測(cè)。