吳東霖,葛偉鵬,2,魏聰敏,張波,2

(1.中國(guó)地震局蘭州地震研究所,甘肅 蘭州 730000;2.蘭州地球物理國(guó)家野外科學(xué)觀測(cè)研究站,甘肅 蘭州 730000)

0 引言

中國(guó)地質(zhì)災(zāi)害分布廣泛且種類繁多,其中滑坡所占比例超過70%[1]?；聻?zāi)害在我國(guó)分布廣泛、發(fā)生頻率高,每年都會(huì)造成巨大的人員傷亡、財(cái)產(chǎn)損失和環(huán)境破壞,尤其是在我國(guó)西南、西北地區(qū)。如2017年6月24日四川茂縣新磨村發(fā)生山體高位滑坡,造成83人失聯(lián)或死亡,64戶農(nóng)房被掩埋,岷江支流堵塞近2 km[2]。因此滑坡災(zāi)害的形變監(jiān)測(cè)和隱患識(shí)別意義重大。

合成孔徑雷達(dá)干涉測(cè)量技術(shù)(InSAR)由于其全天候、全天時(shí)、范圍廣、高精度、低成本等特點(diǎn),近年來逐漸被國(guó)內(nèi)外學(xué)者投入滑坡研究中并取得了較好的效果。1996年Achache等[3]將DInSAR技術(shù)引入滑坡領(lǐng)域,觀測(cè)研究Saint-Etienne-de-Tinee滑坡的形變,得到的結(jié)果與地表測(cè)量結(jié)果具有高度一致性,證明InSAR技術(shù)可用于滑坡形變監(jiān)測(cè);2000年Vietmeier等[4]通過觀測(cè)La Valette滑坡的形變速率,進(jìn)一步證明了InSAR技術(shù)在滑坡變形監(jiān)測(cè)中的可行性和準(zhǔn)確性;2001年Ferretti等[5]利用PS-InSAR技術(shù)監(jiān)測(cè)Ancona地區(qū)滑坡,證實(shí)了PS點(diǎn)形變監(jiān)測(cè)精度可達(dá)1 mm;2004年Hilley等[6]利用PS-InSAR技術(shù)監(jiān)測(cè)美國(guó)Hayward斷層附近區(qū)域,并指出降水可能引發(fā)滑坡;2007年程滔[7]在研究區(qū)布設(shè)角反射器(CR),結(jié)合PS-InSAR對(duì)陜西子長(zhǎng)縣滑坡進(jìn)行了監(jiān)測(cè)研究;2010年王騰[8]提出一種準(zhǔn)永久散射技術(shù)(QPS-InSAR),并以此監(jiān)測(cè)了巴東新城區(qū)的兩個(gè)滑坡;2016年康亞[9]采用InSAR Stacking、IPTA、SBAS-InSAR和PS-InSAR技術(shù)對(duì)西南山區(qū)進(jìn)行形變監(jiān)測(cè),成功探測(cè)出多處滑坡災(zāi)害點(diǎn);2020年張毅等[10]綜合利用InSAR和無人機(jī)測(cè)繪技術(shù)提出一種預(yù)測(cè)潛在滑坡位置及范圍的方法,并得到了很好的驗(yàn)證。

基于上述研究,本文提取2014年10月—2019年12月覆蓋黃河流域劉家峽—蘭州段的升軌Sentinel-1數(shù)據(jù),利用PS-InSAR技術(shù)結(jié)合GPS形變率對(duì)數(shù)據(jù)進(jìn)行改正,得到研究區(qū)LOS向形變速率;將研究結(jié)果與GPS結(jié)果及前人研究進(jìn)行對(duì)比,驗(yàn)證研究結(jié)果的有效性,并進(jìn)一步分析形變典型區(qū)的形變特征和潛在危險(xiǎn)性。

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

1.1 研究區(qū)概況

本文研究區(qū)是黃河流域劉家峽—蘭州段(圖1),主要位于甘肅省臨夏州北部永靖縣,西接青海省民和縣,南與積石山縣、東鄉(xiāng)縣隔黃河相望,東北與蘭州市接壤,東南與定西市接壤。該區(qū)域鄰近青藏高原和黃土高原過渡帶,屬隴西中新生界盆地構(gòu)造,是地震多發(fā)地帶。境內(nèi)平均海拔約2 000 m,地形復(fù)雜,山川交錯(cuò),河谷縱橫,屬盆地邊沿的次高山群及高原淺山丘陵區(qū),大部分地面被黃土覆蓋,地勢(shì)東西高、中部低[11]。區(qū)內(nèi)水資源較為豐富,黃河干流從積石山縣由西南進(jìn)入劉家峽水庫(kù)后轉(zhuǎn)向東北,在劉家峽鎮(zhèn)與支流洮河交匯,后又轉(zhuǎn)向西北,經(jīng)太極鎮(zhèn)、小川—大川盆地后再轉(zhuǎn)向東北,經(jīng)鹽鍋峽鎮(zhèn)與支流湟水匯合,后一路向東進(jìn)入蘭州市區(qū)。研究區(qū)屬大陸型半干旱氣候,陽(yáng)光充足,年日照2 200～2 800 h,氣候溫和,年平均氣溫8～9 ℃,年平均降水280 mm,60%～70%的降水集中在5—8月[11]。

圖1 研究區(qū)活動(dòng)斷裂與地震分布圖(圖中紅緣黑色箭頭表示本文所用7個(gè)GPS站點(diǎn)的水平速度)Fig.1 Distribution of active faults and earthquakes in the study area

1.2 數(shù)據(jù)簡(jiǎn)介

本文選取歐空局2014年發(fā)射的Sentinel-1A衛(wèi)星的干涉寬幅(interferometric wide-swath,IW)模式、Level-1 SLC (single look complex)級(jí)別、升軌的C波段SAR數(shù)據(jù)為數(shù)據(jù)源,其空間分辨率為5 m×20 m,時(shí)間分辨率為12 d。所用SAR數(shù)據(jù)共111景,圖1中紫色矩形框表示其覆蓋范圍(軌道號(hào)為128),監(jiān)測(cè)時(shí)段選取2014年10月—2019年12月。同時(shí)選取歐空局提供的精密軌道星歷數(shù)據(jù)(POD Precise Orbit Ephemerides)和精密軌道輔助文件(Precise Orbit Auxiliary File)。DEM采用由美國(guó)宇航局(NASA)提供的SRTM數(shù)據(jù),其分辨率為30 m×30 m。

隨后引入研究區(qū)GPS速度場(chǎng)來校正InSAR初始結(jié)果,使用的GPS觀測(cè)數(shù)據(jù)來源于中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)(CMONOC)二期工程以及國(guó)際GNSS服務(wù)提供的部分IGS觀測(cè)數(shù)據(jù)、星歷數(shù)據(jù)和表文件數(shù)據(jù)。利用GAMIT[12]進(jìn)行單日松弛解處理,再利用GLOBK[13]軟件將其與SOPAC (Scripps Orbital and Position Analysis Center)提供的ITRF2014全球框架解合并,并通過卡爾曼濾波得到站點(diǎn)在ITRF2014參考框架下的位置和速度,最后再將其轉(zhuǎn)化到鄂爾多斯塊體區(qū)域參考框架下。

2 研究方法

2.1 PS-InSAR基本原理

在SAR影像中,每個(gè)像素的回波是其相關(guān)地面分辨單元中所有散射體的相干和,而數(shù)據(jù)獲取過程中散射體間相對(duì)運(yùn)動(dòng)和雷達(dá)視角的改變會(huì)使相應(yīng)散射體返回能量以不同方式求和,從而引發(fā)影像去相關(guān)[14]。當(dāng)分辨單元內(nèi)某一散射體返回能量遠(yuǎn)大于其他散射體時(shí),對(duì)應(yīng)像素主要受這一個(gè)散射體影響,則去相關(guān)程度會(huì)大大減小,這種散射體即被定義為永久散射體(Persistent scatterer/Permanent scatterer,PS)[14]。人居環(huán)境下常見的PS點(diǎn)有建筑物屋頂、道路等人造地物,人跡稀少環(huán)境下則多為特定方向的裸露巖石。這些目標(biāo)的散射特性較為穩(wěn)定,雷達(dá)回波信號(hào)反射能力很強(qiáng),具有較高信噪比,幾乎不受時(shí)間和空間基線影響,能在較長(zhǎng)時(shí)間內(nèi)保持較高相干性[15]。

PS-InSAR技術(shù)本質(zhì)上仍然屬于D-InSAR,其基本原理是:假設(shè)研究區(qū)內(nèi)有N幅SAR影像,基于空間基線最小、時(shí)間基線最小、多普勒質(zhì)心頻率基線最小等原則選定其中一幅為主影像,其余為輔影像[16],配準(zhǔn)、干涉后可形成N-1幅差分干涉圖,再利用振幅離差閾值、相位離差閾值、小波相位分析等方法從差分干涉序列SAR影像中選取PS點(diǎn),對(duì)這些PS點(diǎn)進(jìn)行時(shí)間序列分析得到其形變量,進(jìn)而得到整個(gè)研究區(qū)的地表形變場(chǎng)。

2.2 趨勢(shì)改正原理

PS-InSAR得到的初始視距向形變速率存在一個(gè)斜坡趨勢(shì)[17],引入測(cè)區(qū)內(nèi)GPS站點(diǎn)速率對(duì)其進(jìn)行改正,同時(shí)也可以給InSAR結(jié)果一個(gè)參考框架。

首先,計(jì)算每個(gè)GPS站周圍圓形(半徑約400 m)范圍內(nèi)InSAR的反距離加權(quán)LOS向速率和入射角。其次,通過式(1)將GPS站速率(VN、VE、VU)轉(zhuǎn)化合成到雷達(dá)視線向(Vlos)[18],并求出其與對(duì)應(yīng)InSAR視距向加權(quán)速率的差值:

Vlos=VNsinθsinα-VEsinθcosα+VUcosθ

(1)

式中:θ和α分別代表雷達(dá)入射角和衛(wèi)星飛行方位角,對(duì)于Sentinel-1升軌,取α=-15°[18],θ即取上步算得加權(quán)入射角。

接著,使用最小二乘法令式(2)中的R2最小以求出趨勢(shì)參數(shù)m1、m2、m3,變形后的m1、m2、m3則由式(3)求出:

(2)

(3)

式中:loni、lati和ΔVi分別代表第i個(gè)GPS站點(diǎn)的經(jīng)度、緯度及其GPS與InSAR視距向速率差值。

最后,移除InSAR視距形變率初始值的線性趨勢(shì)。

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

3.1 數(shù)據(jù)處理

根據(jù)時(shí)間基線和垂直基線最優(yōu)原則選取2017年10月22日的SAR影像為主影像,其余為輔影像,組成110組干涉對(duì)。干涉對(duì)信息如圖2所示,干涉對(duì)空間基線最大為139 m,符合要求。

圖2 時(shí)空基線分布Fig.2 Distribution of temporal-spatial baselines

基于ISCE平臺(tái)[19]進(jìn)行預(yù)處理,再基于StaMPS平臺(tái)[5]采用StaMPS方法來提取PS點(diǎn),最終得到2 884 074個(gè)PS點(diǎn)(振幅離差閾值為0.4),其在城鎮(zhèn)地區(qū)分布較密集,野外山區(qū)分布較稀疏(圖3)?；谶x定的PS點(diǎn)進(jìn)行形變反演,依次去除殘余地形和大氣延遲等影響,解算出各點(diǎn)的形變速率[圖3(a)]。以區(qū)域內(nèi)7個(gè)GPS點(diǎn)三維速率為約束,去除InSAR初始結(jié)果中的線性趨勢(shì)后,得到修正的鄂爾多斯塊體區(qū)域參考框架下的PS點(diǎn)LOS向形變速率[圖3(b)]。

圖3 InSAR LOS向形變速率改正前后對(duì)比圖(圖中暗紅色橢圓為袁道陽(yáng)等[20]推測(cè)的1125年蘭州M7.0地震極震區(qū))Fig.3 InSAR LOS deformation rate before and after calibration

3.2 結(jié)果分析

將改正后的InSAR結(jié)果裁剪出研究區(qū)范圍,得到圖4(a)。由圖4(a)可知,相對(duì)于鄂爾多斯地塊,整個(gè)區(qū)域基本處于LOS沉降狀態(tài),沉降嚴(yán)重區(qū)主要分布在黃河流域及其支流和公路附近,沉降最顯著的區(qū)域位于三條峴鄉(xiāng)下莊村附近(圖4中橢圓區(qū)域),而非研究熱點(diǎn)黑方臺(tái)。將沉降最顯著區(qū)域——下莊沉降區(qū)放大,得到圖(4b)。由圖(4b)可知,下莊沉降區(qū)表現(xiàn)為橢圓形,其南側(cè)是洮河,西、北、東側(cè)均為居民聚居地,如紅峴子、大臺(tái)子、下莊、青和等村落,LOS沉降速率>4 mm/a,面積約25 km2,體量遠(yuǎn)大于黑方臺(tái);沉降最大值點(diǎn)距下莊、青和等村約2 km,距洮河約4 km,距永靖縣城約12 km,距劉家峽水庫(kù)約20 km,距黑方臺(tái)約20 km,距蘭州市區(qū)約40 km。

由圖4(b)可見,沉降速率>10 mm/a區(qū)域也表現(xiàn)為橢圓狀。為了進(jìn)一步分析下莊沉降區(qū)的變形特征,我們沿該橢圓的長(zhǎng)軸(BB′)和短軸(AA′)各作一條剖面(圖5)。如圖5所示,灰色圓點(diǎn)為L(zhǎng)OS向形變速率剖面,綠色實(shí)線為對(duì)應(yīng)的地表高程剖面。根據(jù)形變速率剖面可知,沉降速率>8 mm/a的長(zhǎng)軸約4 km,形成一個(gè)較大的沉降漏斗[圖5(b)]。中心沉降點(diǎn)P1的LOS向形變速率為-14.4 mm/a,約位于橢圓右焦點(diǎn)。結(jié)合高程剖面,山頂點(diǎn)P2和P3的沉降速率約為0 mm/a和-1 mm/a,分別向西和向南延伸至山下達(dá)到最大沉降值,說明若滑坡被觸發(fā),可能會(huì)向南方滑動(dòng),對(duì)東南的居民聚集地造成巨大危害,甚至可能堵塞南側(cè)的洮河。

圖4 研究區(qū)LOS向形變速率Fig.4 LOS deformation rate of the study area

圖5 沉降剖面Fig.5 Subsidence profile

圖3和圖4(a)中疊加了袁道陽(yáng)等[20]推斷的1125年蘭州M7.0地震極震區(qū),可看出下莊沉降區(qū)處于極震區(qū)內(nèi),因此推斷該滑坡隱患點(diǎn)可能與1125年蘭州地震有關(guān)。結(jié)合圖6(a)和圖6(b)中山底點(diǎn)P1與山頂點(diǎn)P2、山頂點(diǎn)P3間的時(shí)間序列差值可看出,2017年7月30日—8月23日沉降差值分別累計(jì)近20 mm和30 mm,之后未回到原位置,產(chǎn)生瞬態(tài)滑移,引起永久變形。2018年8月雖然也存在這樣的趨勢(shì),但差值并非持續(xù)增大,后續(xù)基本回到原位置,故推測(cè)每年8月沉降差值的增大或與降水有關(guān)。2017年8月的差值突增則可能是受到降水和2017年8月8日九寨溝MS7.0地震的雙重影響。

圖6 時(shí)間序列Fig.6 Time series of some special points

通過實(shí)地現(xiàn)場(chǎng)調(diào)查,沉降最大點(diǎn)P1處未見到明顯滑坡形態(tài),僅有一些邊坡失穩(wěn)現(xiàn)象。通過對(duì)沉降區(qū)東南方下莊村居民進(jìn)行走訪,得知不少居民家房屋有裂縫產(chǎn)生,最為明顯的是P4點(diǎn)某居民家墻上產(chǎn)生了約1 cm寬的裂縫(圖7)。房屋建設(shè)已近十年,裂縫約產(chǎn)生于2017年,與時(shí)間序列顯示的形變特征相符[圖6(c)],證明了InSAR結(jié)果的可信度。

圖7 P4點(diǎn)下莊村某居民家墻上的裂縫Fig.7 A crack on the wall of one resident’s house in Xiazhuang village (Point P4)

3.3 精度評(píng)定

為驗(yàn)證InSAR結(jié)果的精度,將GPS臺(tái)站的N、E、U分量合成LOS向與InSAR結(jié)果進(jìn)行對(duì)比,將結(jié)果列于表1。由表1可知其較差基本較小,證明InSAR結(jié)果精度較好。其中G095和G360兩點(diǎn)GPS和InSAR結(jié)果的差值較大,可能是因?yàn)檫@兩點(diǎn)接近圖幅邊緣,受到的約束較差。

表1 GPS與InSAR視距向速率對(duì)比Table 1 Comparison between LOS velocity of GPS and InSAR

最后,對(duì)比本研究得到的黑方臺(tái)地區(qū)的LOS向速率和前人結(jié)果以驗(yàn)證InSAR結(jié)果的可靠性。前人對(duì)黑方臺(tái)地區(qū)形變的研究較為豐富,其中王晨興等[21]得到黑方臺(tái)LOS向速率為-4～-13 mm/a,曾珠[22]得出LOS向區(qū)域均值為-2.7 mm/a,史緒國(guó)等[23]、朱文峰等[24]和張毅等[10]得到的LOS向結(jié)果最大值分別為-60 mm/a、-70 mm/a和-64 mm/a。經(jīng)ArcGIS統(tǒng)計(jì),本研究得到的黑方臺(tái)地區(qū)視線向形變速率為-0.05～-10.09 mm/a,均值為-4.19 mm/a[圖4(a)],與王晨興和曾珠的結(jié)果較為相近,與朱文峰和張毅的結(jié)果相差較大。這可能是因?yàn)榍皟烧吆捅狙芯恳粯硬捎肞S-InSAR方法,而后三者采用SBAS-InSAR方法;在研究期間發(fā)生滑坡的區(qū)域會(huì)失相干,長(zhǎng)時(shí)間序列PS-InSAR方法無法在這些區(qū)域提取到PS點(diǎn),這可能導(dǎo)致較大滑動(dòng)速率觀測(cè)值的丟失[25]。

4 結(jié)論

本文利用PS-InSAR技術(shù)處理111景覆蓋黃河流域劉家峽—蘭州段的Sentinel-1 C-SAR數(shù)據(jù),采用GPS數(shù)據(jù)對(duì)其進(jìn)行趨勢(shì)改正后獲得研究區(qū)2014年10月—2019年12月間的年平均地表LOS向形變場(chǎng);通過對(duì)比發(fā)現(xiàn)InSAR結(jié)果與GPS的較差小于2.3 mm/a,且與同期前人研究較為符合,驗(yàn)證了研究結(jié)果的有效性。

研究區(qū)沉降主要分布在黃河及其支流沿線和交通干線附近,最大沉降區(qū)位于永靖縣三條峴鄉(xiāng)下莊村附近,其規(guī)模和量級(jí)都大于研究熱點(diǎn)黑方臺(tái)。如果被外力(如強(qiáng)降雨、周邊強(qiáng)震等)觸發(fā)滑坡,對(duì)周圍村落和附近河流、交通都將造成嚴(yán)重?fù)p失,應(yīng)予以重視。通過歷史地震和時(shí)間序列分析,此沉降區(qū)可能與1125年蘭州M7.0地震有關(guān),2017年九寨溝M7.0地震可能對(duì)其產(chǎn)生了觸發(fā)作用。

本文通過長(zhǎng)時(shí)間序列PS-InSAR監(jiān)測(cè)獲取了大范圍的滑坡監(jiān)測(cè),但會(huì)丟失一些較大的快速形變,后續(xù)可以考慮和短時(shí)間序列PS-InSAR、SBAS-InSAR或機(jī)載/地基SAR相結(jié)合。本文的結(jié)果是LOS向形變場(chǎng),只能較為粗糙地分析形變,后續(xù)考慮利用升降軌和GPS聯(lián)合解算出該區(qū)域的三維速率場(chǎng)。

致謝:感謝ASF DAAC(https://asf.alaska.edu/)為本文提供的歐空局哥白尼Sentinel-1數(shù)據(jù)[201410—201912];感謝國(guó)家重大科學(xué)工程“中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)”(http://www.neiscn.org)為本文提供的GPS觀測(cè)資料;感謝IGS數(shù)據(jù)中心(http://www.igs.org/)為本文提供的全球IGS站資料;感謝中國(guó)地震臺(tái)網(wǎng)中心國(guó)家地震科學(xué)數(shù)據(jù)中心(http://data.earthquake.cn)提供的數(shù)據(jù)支撐。本文采用ISCE、StaMPS、GAMIT/GLOBK、ArcGIS等軟件進(jìn)行數(shù)據(jù)處理,部分圖件由GMT軟件繪制,在此一并表示感謝！