彭藝，吳云龍，楊玉，袁華清

(1.中國地震局地震研究所，中國地震局地震大地測量重點實驗室，湖北 武漢 430071；2.防災(zāi)科技學(xué)院，河北省地震動力學(xué)重點實驗室，河北 三河 065201)

1 引 言

我國疆土遼闊、物質(zhì)資源豐富，全國各地存在著各種各樣的地質(zhì)條件，高原、山地、盆地、林地等各種地形遍布全國。這些各不相同的地質(zhì)背景，加上各地不同的氣候影響，致使各種地質(zhì)災(zāi)害的發(fā)生，這大大影響著我國的民生建設(shè)和經(jīng)濟發(fā)展。而其中，地面沉降作為一個典型的環(huán)境地質(zhì)災(zāi)害，不僅反映了沉降區(qū)域的地表變化趨勢，也表征了相應(yīng)區(qū)域的地質(zhì)變化規(guī)律，同時也可以在此基礎(chǔ)上探究其他災(zāi)害，如地震、泥石流等的成因和發(fā)展趨勢。因此，合理利用現(xiàn)有的科學(xué)監(jiān)測方法，選定較為典型的研究區(qū)域進行監(jiān)察和預(yù)測，同時對監(jiān)測結(jié)果進行分析并制訂相應(yīng)的應(yīng)對措施，是一項十分重要的研究。傳統(tǒng)滑坡沉降的監(jiān)測主要依賴于人工勘測例如三角高程測量，即使是使用了GNSS等手段獲取數(shù)據(jù)，仍然存在著小范圍、低精度、高成本以及分辨率低和數(shù)據(jù)處理周期長等各種缺陷[1]。隨著衛(wèi)星遙感領(lǐng)域的進步和發(fā)展，基于遙感手段的滑坡沉降監(jiān)測越來越成為一種重要的方法。

合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)[2]技術(shù)是一種星載測量技術(shù)，相比于傳統(tǒng)的測量方式，將Sentinel-1衛(wèi)星影像數(shù)據(jù)經(jīng)相干處理得到沉降區(qū)域的地表形變信息，從而直觀地表征研究測區(qū)的地表沉降情況，具有全天候、全天時、覆蓋廣泛、分辨率高和精度高的優(yōu)點，目前在山體滑坡、三角洲地區(qū)沉降[2]、地震火山、礦區(qū)沉降以及城市地下水開采引發(fā)的地面沉降監(jiān)測等各方面已經(jīng)得到廣泛應(yīng)用。長安大學(xué)的趙延嶺利用干涉圖堆疊技術(shù)和SBAS技術(shù)對長江三峽區(qū)域的滑坡進行了探測與識別，獲取形變信息，分析滑坡成因，做出合理預(yù)測，整體效果比較理想[3]。InSAR技術(shù)主要分為兩種：傳統(tǒng)的差分干涉測量方法(D-InSAR)和時間序列差分干涉測量，其中D-InSAR技術(shù)的構(gòu)建原理是對覆蓋同一地區(qū)的兩幅干涉圖像(一幅形變前的干涉圖像，另一幅形變后的干涉圖像)進行差分處理以獲取目標區(qū)域地表形變[4，5]。在早期的地表形變監(jiān)測中，D-InSAR技術(shù)發(fā)揮著非常重要的作用，但是隨著科學(xué)技術(shù)的發(fā)展，其對時間和基線去相干過程中產(chǎn)生的誤差、干涉圖噪音、大氣效應(yīng)和外部DEM誤差的消除越來越無法達到精度[6]，為了更好地提高影像的相干性，在此基礎(chǔ)上提出SBAS-InSAR小基線技術(shù)，選擇最短基線對影像對進行配準，提高了數(shù)據(jù)處理的精度，同時小基線技術(shù)對于DEM的數(shù)據(jù)精度要求相較于其他InSAR手段并不高，是目前應(yīng)用較為廣泛的一種地表形變監(jiān)測手段[7，8]。

為研究SBAS-InSAR技術(shù)在地面沉降監(jiān)測中的應(yīng)用，本文以四川省茂縣為例，利用SBAS-InSAR技術(shù)研究該地山體滑坡的發(fā)生以及對該地地形產(chǎn)生的影響，進而結(jié)合當(dāng)?shù)氐刭|(zhì)氣候條件推測我國西南地區(qū)的地形變化，再而分析我國全國范圍內(nèi)山體滑坡的趨勢和影響，對經(jīng)濟文化發(fā)展的影響也可以因地制宜地相應(yīng)得到預(yù)測和解決。

2 研究區(qū)域概況

四川省茂縣是四川省阿壩藏族羌族自治州的一個下轄縣，該縣處于四川省的西北部，地質(zhì)上處于青藏高原和川西平原的交界地，地勢變化大，總體呈現(xiàn)西北高東南低的趨勢。茂縣地貌以高山峽谷為主，海拔高低懸殊比較大，氣候受到西風(fēng)環(huán)境和印度洋西南季風(fēng)影響，局部氣候較為復(fù)雜，多種氣候狀況都有。常見各種氣候災(zāi)害，再加上位于岷江上游，水資源豐富，春夏季節(jié)經(jīng)常會有暴雨、冰雹、泥石流等自然災(zāi)害發(fā)生，同時也位于地震多發(fā)區(qū)，極端天氣相對較多，出現(xiàn)大型滑坡的情況也相對較多。

除當(dāng)?shù)氐牡乩須夂颦h(huán)境容易導(dǎo)致山體滑坡等的自然災(zāi)害的發(fā)生，茂縣的地質(zhì)條件特殊，這也決定著當(dāng)?shù)氐刭|(zhì)災(zāi)害的發(fā)生。根據(jù)青藏高原的整體地勢走向，自西向東將青藏高原東緣分為青藏高原地貌區(qū)、龍門山高山地貌區(qū)和四川盆地地貌區(qū)[9]。而導(dǎo)致這些地貌產(chǎn)生的原因是地下的地質(zhì)構(gòu)造，青藏高原東緣地質(zhì)構(gòu)造自西向東分別是松潘-甘孜褶皺帶、龍門山造山帶和四川盆地[10]。本文選擇的四川省茂縣位于青藏高原地貌區(qū)和龍門山高山地貌區(qū)的交界處，即松潘-甘孜褶皺帶和龍門山造山帶的交界[11]，周圍發(fā)生過多次較大的歷史地震。地質(zhì)活動活躍，地形起伏大，再加上岷江上游水源的沖刷，這里形成陡峭的峽谷地勢。此外，茂縣的地層巖石分布著砂巖、板巖、泥頁巖和碳酸鹽等，土質(zhì)易變質(zhì)變形，在地下構(gòu)造變化的過程中極其容易受到受影響，巖石易碎，容易產(chǎn)生地面沉降，如果遇到雨水沖刷，山體滑坡和泥石流無法避免[11，12]。再加上地勢崎嶇不平、交通不便，在自然災(zāi)害出現(xiàn)后經(jīng)常無法第一時間進行救援，而這里土地、礦產(chǎn)和生物資源極其豐富，歷史上也是兵家必爭的要塞之地，所以四川茂縣的山體滑坡現(xiàn)象是我們必須時刻監(jiān)測預(yù)警并制定合理措施的一項工作。

歷史上，四川茂縣發(fā)生過兩次十分嚴重的山體滑坡災(zāi)害，對這里造成了非常大的影響。1933年8月25日，四川省茂縣疊溪發(fā)生的7.5級強烈大地震誘發(fā)山體滑坡，幾乎將這個位于四川省西部的羌族聚居區(qū)全部掩埋，只剩一座破敗不堪的城隍廟，據(jù)不完全統(tǒng)計，死亡人數(shù)約有 2 500人，造成嚴重損失。2017年6月24日，茂縣疊溪鎮(zhèn)新磨村的山體突發(fā)高位垮塌，被埋100多人，還造成了 2 km的堵塞河道，1 600余米的道路掩埋，以及光纜受損有 3 km長的損害，基本中斷基礎(chǔ)通信。經(jīng)初步分析，這是一起由強降雨誘發(fā)的山體滑坡，再加上茂縣地處地震斷裂帶以及汶川地震對山體的影響，在連日降雨的情況下發(fā)生了山體滑坡。

3 研究方法

針對四川茂縣的山體滑坡監(jiān)測，本文采用合成孔徑雷達(Interferometric Synthetic Aperture Radar，InSAR)的技術(shù)進行測量，相對于一些需要人工進行操作的三角高程測量、GNSS測量等，InSAR技術(shù)具有全天候、全天時、覆蓋廣泛、分辨率高和精度高的優(yōu)點，是目前對山體滑坡進行監(jiān)測精度最高也最合理的方法。在2002年，Berardino等人在InSAR技術(shù)的基礎(chǔ)上提出了小基線技術(shù)，獲取地表形變以研究大尺度上的形變[8]，后經(jīng)Hopper等人逐步優(yōu)化，逐步提高了傳統(tǒng)InSAR技術(shù)時間采樣率，減弱空間失相干造成的影響和地形對差分的影響，增強了形變監(jiān)測的可靠性，在地表形變監(jiān)測中發(fā)揮著非常重要的作用。

本文選取了2017年3月15日～6月19日期間四川茂縣以及周邊區(qū)域的8景Sentinel-1A數(shù)據(jù)，利用Hooper[13]提出的SBAS-InSAR數(shù)據(jù)處理模型進行處理，對影像數(shù)據(jù)進行匹配，然后進行干涉處理、濾波處理和解纏，得到最小基線構(gòu)成的三角網(wǎng)后，再進行反演得到地表形變的速率，最后提取時間序列曲線，并進行時序分析，找到滑坡和疑似滑坡的區(qū)域。

3.1 算法模型

在進行反演之前，要先對影像的相位進行計算。根據(jù)衛(wèi)星影像獲取的N+1幅SAR影像，形成M幅干涉圖，生成干涉相位，再利用外部數(shù)據(jù)去除地形相位影響，根據(jù)最小基線原則，利用三維相位解纏的方法進行解纏，之后進行濾波處理以消除各種影響。

在相應(yīng)時間間隔內(nèi)，使用同一傳感器獲取覆蓋同一地區(qū)的N+1幅影像，根據(jù)集合內(nèi)基線距最小、集合間基線距最大的原則，通過給定閾值，可生成干涉圖M幅(M和N滿足(N+1)/2≤M≤[N(N+1)]/2)，第i幅干涉圖任意像元(x，y)處得到的干涉相位可表示為：

фobs，i(x，y)=фtopo，i(x，y)+фdefo，i(x，y)+фatm，i(x，y)+фnoise，i(x，y)

(1)

其中фobs表示纏繞相位，фtopo表示地形相位，фdefo表示DEM誤差引起的相位，фatm表示大氣噪聲引起的相位，фnoise表示一起熱噪聲等其他相位，這些相位都可以通過濾波來去除。

之后再經(jīng)過圖像匹配、干涉處理、去除地形相位、解纏計算、濾波處理后進行最小二乘計算，得到地表形變速率。其中最小二乘計算的公式如下：

AVlos=φobs

(2)

Vlos=A+φobs

(3)

(4)

提取形變時序后，在各相位點上聯(lián)立方程組，利用最小二乘進行解算。利用式(2)求矩陣A+，使用的是矩陣的奇異值分解法(SVD)求其廣義逆，即A=USVT。再利用形變量的差和時間差得到各個時間點的形變速率矩陣即式(4)，帶入式(3)可得到形變速率的矩陣[1，2]，即通過A+=VS+UT求得A+，最終求得Vlos[6，8，13，14，15]。

3.2 計算流程

本文使用SBAS-InSAR模型，對衛(wèi)星影像圖完成干涉處理、相位解纏和濾波處理，最后得到形變速率并分析其形變規(guī)律[16]。即先選擇的主影像，將副影像與其進行配準，并進行干涉處理，參照DEM高程數(shù)據(jù)去除地形相位影響，根據(jù)小基線原則，生成小基線集干涉圖，設(shè)定閾值后篩選SDFP候選點，通過估算其穩(wěn)定性，選取穩(wěn)定性高的點以去除估算誤差，在將此時得到的形變差利用三維相位解纏的方法進行解纏；最后進行濾波處理去除噪聲，提取形變時序，利用最小二乘法，同時利用矩陣的奇異值分解法建立方程組運算。具體流程如圖1所示：

圖1 SBAS處理數(shù)據(jù)流程圖

4 實驗結(jié)果

本文選取2017年3月15日～6月19日期間茂縣8幅Sentinel-1A影像數(shù)據(jù)，根據(jù)茂縣地理位置103°39′46″E，32°4′47″N，裁剪至合適大小，基于SBAS-InSAR模型，將各影像兩兩連接匹配，生成連接圖和時間基線，并自動定義超級主影像為20170315，再將各像對之間的匹配之后會對相距較遠的兩個數(shù)據(jù)進行解纏，得到最小基線的影像對配準三角網(wǎng)(圖2)。

圖2 影像對連接圖、時間基線分布圖、解纏后三角網(wǎng)

將解纏后的影像對進行干涉處理、濾波處理和3D解纏。參照拼接好的茂縣區(qū)域的DEM高程數(shù)據(jù)，選取閾值為0.2進行DEM法干涉處理，同時選擇高斯濾波去除噪聲提高影像對相干性，并使用最小費用法進行3D解纏。本文選擇相干性較好的影像對20170315-20170327的相干系數(shù)圖、干涉圖、解纏后影像圖進行展示(圖3)，并剔除相干性不好的影像對。

圖3 20170408-20170514相干系數(shù)圖、干涉圖、解纏后影像圖

選取相干性較好的影像對20170408-20170514，選取45個GCP點進行軌道精化和重去平以提高影像對的配準精度，選取GCP點時遵循盡量均勻分布，在未發(fā)生形變區(qū)域處人工選取的原則。經(jīng)兩次反演，分別進行二次解纏，以達到優(yōu)化干涉像對、去除DEM殘余誤差的目的，以此反演出地表形變速率和殘余地形，以及使用大氣濾波以去除大氣信息、估算形變組分。將兩次反演的SAR坐標系下地表信息轉(zhuǎn)換至WGS-84坐標系下，通過內(nèi)插和平滑的方法進行地理編碼，得出地表形變平均速率，經(jīng)密度分割后選取測區(qū)特征點繪制時間序列曲線分析其走勢，并結(jié)合地質(zhì)氣候條件進行分析。

5 結(jié)果分析

將測區(qū)的地表形變速率進行密度分割(圖4)，該區(qū)域最大抬升速率大約為 50 mm/a，即圖中紅色的部分，而且分布在這里的點非常少，說明出現(xiàn)抬升的地方相對來說比較少。同時也可以看出，最大沉降速率可達到 -100 mm/a，如圖4所示藍色部分分布的點較多，說明這段時間內(nèi)出現(xiàn)了比較嚴重的滑坡。

圖4 形變速率分布圖

同時選取4個地勢較高地特征點和1個地勢較低地特征點繪制時間序列曲線(圖5，圖6)，分析監(jiān)測區(qū)域的地面沉降狀況。5個特征點的地表在3月～4月初整體形變十分微小，基本沒有變化，甚至?xí)形⑿〉牡乇硖?，但是?月～5月初，出現(xiàn)了明顯的高程值變化，特征點處呈現(xiàn)較陡的曲線走向，說此處發(fā)生了一定程度的山體滑坡，不過在5月中旬，5個特征點的高程變化值又集中在0的上下，說明5月中旬這里沒有發(fā)生大的地勢形變事件，即沒有影響較大的自然災(zāi)害發(fā)生，而在5月末～6月中旬這段時間，時間序列曲線卻出現(xiàn)了驟變，如圖前4個特征點的時間序列曲線急速下降為負值，第5個點的時間序列曲線迅速上升，并且按照曲線的走向，在6月19日之后仍會按照這個趨勢下降和上升，說明此段時間內(nèi)此地發(fā)生了較為嚴重的山體滑坡，并且隨著時間的增加，山體滑坡會繼續(xù)發(fā)生，而且甚至可能會有更加嚴重的自然災(zāi)害發(fā)生，可能會對附近居民的生活以及道路產(chǎn)生十分嚴重的影響。

圖5 特征點1-4時間序列曲線

圖6 特征點5時間序列曲線

四川茂縣的山體形變發(fā)生上述地勢變化，與其所處的地勢條件和該處的氣候變化有著非常大的關(guān)系。這里地處四川盆地西北部，山地較多，四周分布著海拔為 1 000 m～3 000 m的青藏高原、云貴高原、巫山、大巴山和大婁山等幾座高山，周圍群山環(huán)繞，西北邊緣是很長的龍門山脈，中間地勢低矮，可明顯分為邊緣山地和中間盆地兩部分。同時，四川盆地也是我國著名的紅層盆地，其地表巖石主要是由紫紅色砂巖和頁巖組成，它們極易風(fēng)化發(fā)育成紫色土，土壤中含有豐富的鈣、磷、鉀等營養(yǎng)元素，土壤極其肥沃，但是土層淺薄，母巖疏松，非常易于崩解，很容易在大雨天氣隨著雨水被沖走造成泥石流，嚴重的甚至可以造成山體滑坡。除此之外，這里屬于亞熱帶季風(fēng)性濕潤氣候，夏季平均溫度24℃～28℃，極端高溫達到36℃～42℃，年平均降水量在 1 000 mm～1 300 mm之間，但是年內(nèi)降水分布不均，整體氣候呈現(xiàn)冬干、春旱、夏澇、秋綿雨的分布趨勢，70%～75%的雨量都集中在6月～10月的夏季，高溫多雨的夏季很容易發(fā)生連續(xù)暴雨天氣，由此造成山體滑坡的大型自然災(zāi)害現(xiàn)象。

所以正是由于四川茂縣處于這樣的地形和氣候條件中，才會出現(xiàn)如圖所示的時間序列的變化。圖中3月份5個特征點的高程變化值基本都集中在0的上下，會有微小的形變，因為四川盆地的地表巖石易于風(fēng)化，而此時正值春天，相對來說，風(fēng)比較大，再加上紫色土土質(zhì)疏松，容易崩解，巖石的風(fēng)化以及土壤的流失都是這里地勢發(fā)生變化的原因。根據(jù)2017年全年茂縣月降水量(圖7)分析，4月份開始進入夏天，雨水天氣逐漸增多，處在山脈地區(qū)再加上土壤疏松，泥土隨著雨水就從地勢高的地方流失到地勢較低的地方，也就產(chǎn)生了圖中所示的前4個特征點高程變化值是負值，而第5個特征點的高程變化值則為正的情況，且明顯從圖中分析得到，第5個特征點的地勢相較于前4個都比較低。5月中旬，根據(jù)圖中茂縣平均氣溫分布圖以及降水量分布圖，可以發(fā)現(xiàn)5月份茂縣的平均溫度并不算高，而且降水量也相對較少，因為5月17日～18日這段時間，茂縣遭受了大風(fēng)、冰雹自然災(zāi)害的襲擊，所以地勢沒有大規(guī)模的變化。

圖7 2017年茂縣月降水量和平均溫度圖

5月中旬的冰雹等極端天氣結(jié)束之后，正式進入夏季，在亞熱帶季風(fēng)性濕潤氣候的影響下，氣溫持續(xù)上升，降水量也開始持續(xù)增加。此后降水集中，強降雨天氣持續(xù)，在雨水的沖刷下，疏松的土壤開始發(fā)生滑坡，所以5月中旬以后，地勢高的4個特征點的高程變化值發(fā)生驟變且為負值，而地勢較高的第5個特征點則發(fā)生了抬升。因為當(dāng)?shù)氐貏葑兓?，歷史地震致使茂縣巖體結(jié)構(gòu)發(fā)生破壞和裂縫，在持續(xù)的強降雨的沖刷下，具有單面山斜坡結(jié)構(gòu)的后緣巖層被拉斷[12]，發(fā)生了順層滑動，誘發(fā)持續(xù)的泥土流失。截止到本次試驗研究的6月19日，特征點持續(xù)發(fā)生沉降和抬升，并且按照時間序列曲線，在未來的時間，該地區(qū)還會按照這個趨勢繼續(xù)沉降和抬升，而且曲線變化很陡，很有可能發(fā)生更為嚴重的自然災(zāi)害，如大型山體滑坡。而事實也確實如此，2017年6月24日6時左右，茂縣疊溪鎮(zhèn)新磨村山體突發(fā)了高位垮塌，道路和村民被掩埋，是一起十分嚴重的自然災(zāi)害事件，對當(dāng)?shù)氐木用裨斐闪耸謬乐氐娜松砗拓敭a(chǎn)傷害。

6 結(jié) 論

經(jīng)過上述的實驗過程，對2017年3月～6月之間四川茂縣的衛(wèi)星影像數(shù)據(jù)，基于小基線技術(shù)，對衛(wèi)星影像進行了干涉，濾波和解纏等，再進行兩次反演得到其形變速率，最后對其進行地理編碼，找到特征點的時間序列，通過不同時間段的高程值變化來研究該處的形變情況，并同時結(jié)合當(dāng)?shù)厝暝缕骄邓蜌鉁財?shù)據(jù)分析其成因，最終得到的結(jié)論有：

(1)四川省茂縣在2017年3月15日～6月19日期間，地表發(fā)生了一定程度的沉降，最大沉降速率和抬升速率分別達到 -100 mm/a和 50 mm/a，且地面沉降隨著降水和溫度的變化而改變。

(2)3月茂縣附近整體地表形變不大，直至5月中旬，出現(xiàn)一定幅度的沉降，到5月中旬出現(xiàn)停止沉降的現(xiàn)象，之后到6月，沉降急速增加，并且有按照這個速度持續(xù)沉降的趨勢。

(3)結(jié)合2017年茂縣月平均降水和溫度數(shù)據(jù)，以及當(dāng)?shù)氐刭|(zhì)條件和構(gòu)造，分析得出：處在地質(zhì)構(gòu)造帶交界處的茂縣地質(zhì)條件特殊，巖體結(jié)構(gòu)已遭受破壞，持續(xù)的降水容易誘發(fā)山體滑坡的發(fā)生。

以上就是基于小基線技術(shù)對3月～6月四川茂縣衛(wèi)星影像數(shù)據(jù)的形變速率相干性分析，并結(jié)合當(dāng)年月平均溫度和降水?dāng)?shù)據(jù)以及當(dāng)?shù)氐牡刭|(zhì)構(gòu)造背景分析，得出了此處可能會在將來發(fā)生大型山體滑坡的結(jié)論，整體吻合當(dāng)?shù)鼐唧w滑坡變化。依據(jù)SBAS-InSAR技術(shù)對地質(zhì)活動活躍的地區(qū)進行定期監(jiān)測，并結(jié)合當(dāng)?shù)氐刭|(zhì)條件和降水?dāng)?shù)據(jù)做出分析和預(yù)測，方法合理且結(jié)果可信，對城市的發(fā)展和規(guī)避一些自然災(zāi)害帶來的損失具有重要意義。