楊正榮，喜文飛,2，史正濤，肖 波，周定義

（1.云南師范大學(xué)地理學(xué)部，云南 昆明 650500；2.云南省高校高烈度地震山區(qū)交通走廊工程地質(zhì)病害早期快速判識(shí)與防控重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650093；3.昆明理工大學(xué)國(guó)土資源工程學(xué)院，云南 昆明 650093；4.云南交通職業(yè)技術(shù)學(xué)院公路與建筑工程學(xué)院，云南 昆明 650500）

0 引言

庫(kù)岸滑坡是深切割高山峽谷型庫(kù)岸常見(jiàn)的破壞形式，多集中分布于我國(guó)西南山區(qū)[1-2]。庫(kù)岸滑坡由地表內(nèi)外營(yíng)力相互作用而形成，人類工程建設(shè)及庫(kù)區(qū)水位變化則使其演化特征更為突出，通常表現(xiàn)為庫(kù)區(qū)蓄水之后庫(kù)岸下緣坡體巖層軟化，從而引起上部庫(kù)岸的形變破壞[3-4]。庫(kù)岸滑坡受多方面因素影響，具有成因復(fù)雜、類型多樣和危害巨大等特點(diǎn)[5-7]。在水電站庫(kù)區(qū)，蓄水和泄洪等因素導(dǎo)致的水位變化直接影響庫(kù)岸滑坡的穩(wěn)定性，庫(kù)岸滑坡一旦失穩(wěn)會(huì)誘發(fā)一系列次生災(zāi)害，破壞區(qū)域生態(tài)系統(tǒng)，毀壞庫(kù)區(qū)壩體和發(fā)電設(shè)施[8]，嚴(yán)重威脅庫(kù)區(qū)上下游居民生命財(cái)產(chǎn)安全。因此，對(duì)水電站庫(kù)岸滑坡進(jìn)行變形監(jiān)測(cè)具有重要意義。

常規(guī)監(jiān)測(cè)手段已難以識(shí)別和監(jiān)測(cè)大面積的庫(kù)岸滑坡形變，相對(duì)傳統(tǒng)的精密水準(zhǔn)測(cè)量、全球衛(wèi)星導(dǎo)航系統(tǒng)（Global Navigation Satellite System,GNSS）和光學(xué)遙感技術(shù)，合成孔徑雷達(dá)干涉測(cè)量（Interferometric Synthetic Aperture Radar，InSAR）技術(shù)因其大范圍、全天時(shí)、全天候、高精度和高分辨率等特點(diǎn)，已成功應(yīng)用于庫(kù)岸滑坡災(zāi)害識(shí)別監(jiān)測(cè)分析中[9-11]。國(guó)內(nèi)外學(xué)者利用InSAR 技術(shù)在庫(kù)岸滑坡災(zāi)害的應(yīng)用方面做了大量研究，康亞等[8]利用三類InSAR 產(chǎn)品和DEM 數(shù)據(jù)對(duì)金沙江流域?yàn)鯑|德水電站段進(jìn)行滑坡早期識(shí)別，成功探測(cè)到多處未知和已知的滑坡體；徐帥等利用墨蘭指數(shù)對(duì)SBAS-InSAR 技術(shù)獲取的形變點(diǎn)進(jìn)行空間域異常值分析和聚類處理，成功識(shí)別出三峽庫(kù)區(qū)巫山—奉節(jié)段高概率潛在滑坡范圍[12]；王振林等[13]利用SBAS-InSAR 技術(shù)提取雅礱江流域錦屏一級(jí)水電站庫(kù)區(qū)左岸邊坡的形變特征信息，推斷出大幅水位上升是誘發(fā)滑坡復(fù)活的主要因素；朱同同等[14]結(jié)合時(shí)序InSAR 技術(shù)和GPS 觀測(cè)值分析了降雨和蓄水對(duì)三峽庫(kù)區(qū)樹坪滑坡變形的影響；史緒國(guó)等[15]聯(lián)合分布式目標(biāo)與點(diǎn)目標(biāo)的時(shí)序InSAR 技術(shù)對(duì)三峽庫(kù)區(qū)藕塘滑坡進(jìn)行穩(wěn)定性監(jiān)測(cè)；Tantianuparp 等[16]聯(lián)合多種SAR 數(shù)據(jù)和PS-InSAR 技術(shù)對(duì)三峽巴東進(jìn)行滑坡探測(cè)，并將PS 點(diǎn)時(shí)序形變與水位時(shí)間變化進(jìn)行初步相關(guān)性分析；Zhou 等[17]利用時(shí)序InSAR 技術(shù)發(fā)現(xiàn)三峽庫(kù)區(qū)木魚堡滑坡變形主要發(fā)生在水庫(kù)漲落期和高水位期；Liu等[18]利用SBAS-InSAR 技術(shù)對(duì)三峽巴東地區(qū)進(jìn)行滑坡探測(cè)并分析季節(jié)性滑坡運(yùn)動(dòng)與水位變化之間的相關(guān)性，上述研究證明了時(shí)序InSAR 技術(shù)在庫(kù)岸滑坡監(jiān)測(cè)中的可靠性，可以對(duì)水電站庫(kù)岸滑坡變形進(jìn)行有效分析。白鶴灘水電站地處四川和云南交界，自2021年4月開(kāi)始蓄水，庫(kù)區(qū)水位由660 m 升至825 m，上升幅度達(dá)165 m；水電站運(yùn)行期間最低水位765 m，最高水位825 m，升降水位差60 m，最大庫(kù)容達(dá)256×108m3[19-20]。庫(kù)區(qū)地形起伏較大、斷裂構(gòu)造發(fā)育，加之蓄水引起的水位變化直接影響庫(kù)岸潛在滑坡的變形趨勢(shì)，對(duì)水電站基礎(chǔ)設(shè)施和上下游居民生命財(cái)產(chǎn)安全造成潛在威脅[21-22]。因此，亟需對(duì)白鶴灘水電站庫(kù)岸潛在滑坡進(jìn)行變形分析。

文章聯(lián)合2019年7月3日至2021年7月28日的150景升降軌Sentinel-1 SAR 數(shù)據(jù)集，采用SBAS-InSAR 技術(shù)獲取白鶴灘水電站庫(kù)區(qū)雷達(dá)視線方向（Line of sight，LOS）形變時(shí)間序列，在分析地表形變時(shí)間演化規(guī)律和空間分布特征的基礎(chǔ)上，結(jié)合無(wú)人機(jī)野外調(diào)查，分析白鶴灘水電站庫(kù)岸潛在滑坡的變形特征，重點(diǎn)研究蓄水因素對(duì)庫(kù)岸潛在滑坡變形趨勢(shì)的影響。

1 SBAS-InSAR 技術(shù)

小基線集InSAR（Small Baseline Subset InSAR，SBASInSAR）技術(shù)最早由Berardino 和Lanari 等[23]提出，該方法通過(guò)組合數(shù)據(jù)的方式獲得一系列短空間基線差分干涉圖，這些差分干涉圖能較好地克服空間失相關(guān)現(xiàn)象。SBAS-InSAR 技術(shù)利用奇異值分解（SVD）法求解形變速率，將被較大空間基線分開(kāi)的孤立SAR 數(shù)據(jù)進(jìn)行連接，進(jìn)一步提高觀測(cè)數(shù)據(jù)的時(shí)間采樣率[23]。該方法可以有效減弱大氣效應(yīng)，降低相位噪聲和誤差[24]，其基本原理及流程如下：假定已獲取覆蓋同一區(qū)域的按時(shí)間序列排序的N+1幅SAR 影像：

根據(jù)干涉組合規(guī)則，生成M幅干涉圖且M應(yīng)當(dāng)滿足：

假設(shè)以t0作 為影像獲取起始時(shí)刻且t0時(shí)刻影像覆蓋區(qū)域位移為0，則在去除軌道誤差、平地效應(yīng)及地形相位的影響后，第i(1≤i≤M)幅影像某像素的干涉相位可表示為：

式中：Δφidef——斜距向形變產(chǎn)生的相位；

Δφitopo——地形相位；

Δφiatm——大氣延遲引起的相位；

Δφinoise——相干噪聲造成的相位。

利用最小二乘或者奇異值分解（SVD）對(duì)m 個(gè)解纏相位進(jìn)行三維時(shí)空相位解纏即可獲得不同SAR 時(shí)刻對(duì)應(yīng)的時(shí)序形變速率。

2 研究區(qū)概況和研究數(shù)據(jù)

本文以四川省與云南省交界白鶴灘水電站庫(kù)區(qū)作為研究區(qū)域，如圖1所示。研究區(qū)長(zhǎng)約30.38 km，寬約11.65 km，總面積353.93 km2，地處橫斷山脈東北部、青藏高原東南邊緣，區(qū)域內(nèi)斷裂構(gòu)造發(fā)育，構(gòu)造運(yùn)動(dòng)強(qiáng)烈，河谷深切，山體陡峻，地震頻發(fā)[25-27]。最高海拔3 556 m，最低海拔520 m，高差達(dá)3 036 m，地勢(shì)陡峭，致使該區(qū)存在大量滑坡、崩塌和泥石流等地質(zhì)災(zāi)害隱患。

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

形變監(jiān)測(cè)數(shù)據(jù)選用從歐州航天局(European Space Agency,ESA)免費(fèi)下載的150 景C 波段Sentinel-1 雷達(dá)影像（其中升軌數(shù)據(jù)50 景，降軌數(shù)據(jù)100 景并在每個(gè)時(shí)間點(diǎn)上下兩景拼接），升降軌數(shù)據(jù)覆蓋區(qū)域如圖2所示。時(shí)間跨度為2019年7月3日至2021年7月28日，極化方式為VV，成像方式為IW，數(shù)據(jù)參數(shù)如表1所示。為提高影像軌道精度，引入POD 精密定軌星歷數(shù)據(jù)。使用日本宇宙航空研究開(kāi)發(fā)機(jī)構(gòu)（Japan Aerospace Exploration Agency，JAXA）發(fā)布的ALOS WORLD 3D 30 m空間分辨率的數(shù)字高程模型（Digital Elevation Model，DEM），用于去除地形相位影響，如圖3所示。

圖2 SAR 衛(wèi)星影像覆蓋范圍Fig.2 SAR satellite image coverage

圖3 研究區(qū)DEMFig.3 Digital elevation model of study area

表1 Sentinel-1A 數(shù)據(jù)參數(shù)Table 1 Sentinel-1A data parameters

3 SBAS-InSAR 技術(shù)數(shù)據(jù)處理

采用SBAS-InSAR 技術(shù)，選取經(jīng)鑲嵌、配準(zhǔn)和裁剪后的100 景Sentinel-1A 斜距單視復(fù)數(shù)（Single Look Complex，SLC）影像（升降軌數(shù)據(jù)各50 景），根據(jù)時(shí)間基線和垂直基線最優(yōu)原則，升軌和降軌數(shù)據(jù)分別以日期為20 191 216 和20 200 204 的影像作為超級(jí)主影像。設(shè)置時(shí)間基線閾值180d，空間基線為臨界基線閾值的50%，共生成654 和888 對(duì)干涉像對(duì)。為抑制斑點(diǎn)噪聲，設(shè)置多視數(shù)為1∶4，采用Minimum Cost Flow 解纏方法和Goldstein 濾波方法進(jìn)行干涉處理，將組合干涉對(duì)經(jīng)過(guò)配準(zhǔn)，調(diào)整刪除不理想的數(shù)據(jù)后生成干涉圖，研究區(qū)部分較理想的干涉圖如圖4所示。

圖4 研究區(qū)部分較理想的干涉圖Fig.4 Ideal interference patterns in the study area

經(jīng)過(guò)軌道精煉和重去平，利用最小二乘法和奇異值矩陣分解進(jìn)行形變反演，然后估算和去除大氣相位，得到研究區(qū)時(shí)間序列形變信息，對(duì)時(shí)序信息地理編碼后獲取研究區(qū)2019年7月3日至2021年7月28日LOS 方向的形變結(jié)果。如圖5所示，形變速率為正值表示靠近衛(wèi)星，負(fù)值表示遠(yuǎn)離衛(wèi)星。對(duì)比圖5（a）、（b）研究區(qū)形變結(jié)果可知，降軌數(shù)據(jù)集探測(cè)的形變信息較為豐富，主要集中在庫(kù)區(qū)西岸，最大LOS 向形變速率-61.425 mm/a；升軌數(shù)據(jù)集僅在庫(kù)區(qū)東岸部分區(qū)域形變較為明顯，最大LOS 向形變速率為91.426 mm/a。升降軌數(shù)據(jù)集形變信息不一致的原因是白鶴灘水電站庫(kù)區(qū)兩岸地形起伏較大，山勢(shì)陡峭險(xiǎn)峻，而升軌數(shù)據(jù)飛行方向大致沿東南向西北，雷達(dá)視線方向位于右側(cè)，降軌數(shù)據(jù)則與之相反，故利用InSAR 探測(cè)形變過(guò)程中陰影、疊掩和透視收縮等幾何畸變現(xiàn)象嚴(yán)重。

圖5 研究區(qū)視線向形變速率Fig.5 Line-of-sight deformation rate of the study area

4 試驗(yàn)結(jié)果與分析

4.1 研究區(qū)庫(kù)岸典型潛在滑坡選取

對(duì)升軌和降軌數(shù)據(jù)獲取的研究區(qū)形變結(jié)果進(jìn)行綜合解譯，升軌數(shù)據(jù)庫(kù)岸形變區(qū)域解譯結(jié)果如圖6所示，共選取庫(kù)岸形變較大區(qū)域4 處。結(jié)合無(wú)人機(jī)野外調(diào)查結(jié)果，發(fā)現(xiàn)典型潛在滑坡2 處，分別用H1 和H2 表示；非滑坡形變區(qū)2 處，分別用X1 和X2 表示，升軌數(shù)據(jù)詳細(xì)解譯結(jié)果如表2所示。

圖6 升軌潛在滑坡解譯及實(shí)地考察結(jié)果Fig.6 Interpretation and field investigation results of potential landslide in ascending orbit

表2 升軌數(shù)據(jù)庫(kù)岸形變區(qū)域解譯結(jié)果列表Table 2 List of interpretation results of shore deformation region in orbit lifting database

降軌數(shù)據(jù)庫(kù)岸形變區(qū)域解譯結(jié)果如圖7所示，共選取庫(kù)岸形變較大區(qū)域6 處。結(jié)合無(wú)人機(jī)野外調(diào)查結(jié)果，發(fā)現(xiàn)典型潛在滑坡4 處，分別用H3 至H6 編號(hào)；非滑坡形變區(qū)2 處，分別用X3 和X4 表示，降軌數(shù)據(jù)詳細(xì)解譯結(jié)果如表3所示。

圖7 降軌潛在滑坡解譯及實(shí)地考察結(jié)果Fig.7 Interpretation and field investigation results of potential landslide in descending orbit

表3 降軌數(shù)據(jù)庫(kù)岸形變區(qū)域解譯結(jié)果列表Table 3 List of interpretation results of shore deformation region in orbit descent database

對(duì)比升軌和降軌數(shù)據(jù)解譯結(jié)果可以看出，非滑坡形變區(qū)X1、X2 與X3、X4 分別相同，潛在滑坡H1、H2 與H3、H5 相互對(duì)應(yīng)。另外，降軌數(shù)據(jù)還解譯出除上述區(qū)域以外的潛在滑坡H4 和H6，同一時(shí)間段不同軌道SAR 數(shù)據(jù)集探測(cè)的形變結(jié)果能夠相互對(duì)應(yīng)，從側(cè)面驗(yàn)證了本文InSAR 結(jié)果的準(zhǔn)確性，但受時(shí)間、空間失相干因素和幾何畸變影響，升降軌形變信息有所差異，說(shuō)明升降軌結(jié)合的方式能夠有效彌補(bǔ)僅利用單一軌道識(shí)別結(jié)果不全面、不準(zhǔn)確的缺陷，提升庫(kù)岸潛在滑坡災(zāi)害識(shí)別和監(jiān)測(cè)的準(zhǔn)確性和有效性。

4.2 庫(kù)岸典型潛在滑坡變形分析

結(jié)合4.1 節(jié)升降軌數(shù)據(jù)集庫(kù)岸潛在滑坡解譯結(jié)果，本文選取H1、H2、H4 和H6 四處典型潛在滑坡進(jìn)行變形分析，分別在各滑坡形變結(jié)果中選取特征點(diǎn)，引入研究區(qū)降雨數(shù)據(jù)，繪制特征點(diǎn)在蓄水前后的時(shí)序形變曲線，并結(jié)合無(wú)人機(jī)野外調(diào)查結(jié)果分析庫(kù)岸典型潛在滑坡的變形特征。

H1 滑坡地處觀音巖，位于沿江公路東岸，滑坡形變速率如圖8（a）所示?；抡w形變速率范圍為-10.726～15.433 mm/a，分別選取滑坡體上緣和下緣特征點(diǎn)A、B 與降雨數(shù)據(jù)構(gòu)建時(shí)序形變曲線如圖8（b）所示，特征點(diǎn)A 和B 時(shí)序形變速率波動(dòng)趨勢(shì)大致相同，每年雨季形變速率較旱季明顯增加。2019年10月—2020年5月形變速率減小，2021年4月后形變速率增大，同比增加約16 mm/a。

圖8 H1 潛在滑坡形變特征Fig.8 H1 potential landslide deformation characteristics

經(jīng)實(shí)地勘察，該滑坡坡體上緣為自然坡體，坡體下緣已進(jìn)行邊坡加固，故在2019年10月至2020年5月間B 點(diǎn)較A 點(diǎn)形變速率變化相對(duì)穩(wěn)定。受降雨因素影響，坡體在雨季滑動(dòng)速率增大。2021年4月至5月，降雨量幾乎為零，水電站蓄水導(dǎo)致庫(kù)區(qū)水位上升，庫(kù)岸下緣受到江水侵蝕改變坡體上下緣間的平衡關(guān)系，使該滑坡體形變量增大。

H2 滑坡地處魚壩村，金沙江支流末端?；滦巫兯俾嗜鐖D9（a）所示，形變較大值處于坡體中上部，形變范圍在-19.326～8.254 mm/a。在坡體兩側(cè)分別選擇特征點(diǎn)C、D 結(jié)合降雨數(shù)據(jù)構(gòu)建時(shí)序形變曲線見(jiàn)圖9（b），特征點(diǎn)C 和D 形變速率變化趨勢(shì)基本一致，與降雨數(shù)據(jù)呈現(xiàn)一定相關(guān)性，雨旱兩季形變速率差異較小。2020年1月至10月間，形變速率逐漸減小，2021年1月后形變速率振蕩變化，2021年4月之后，形變速率增加值超過(guò)10 mm/a。

圖9 H2 潛在滑坡形變特征Fig.9 H2 potential landslide deformation characteristics

經(jīng)野外實(shí)地調(diào)查，H2 滑坡滑面自上而下呈倒“V”字形。由于隧道工程尚未完工，附近仍伴有部分工程活動(dòng)，故在2020年雨季坡體滑動(dòng)速率對(duì)降雨因素響應(yīng)較弱，表現(xiàn)為形變速率逐漸減小。從圖9（b）可以看出，2021年4月以后，特征點(diǎn)C 和D 形變速率相對(duì)之前有所增加，此時(shí)降雨量較小，說(shuō)明蓄水導(dǎo)致的庫(kù)區(qū)水位抬升也對(duì)遠(yuǎn)離河道的坡體產(chǎn)生影響。

H4 滑坡體地處清水溝，位于庫(kù)區(qū)西岸?；滦巫兯俾嗜鐖D10（a）所示，滑坡整體形變范圍為-17.605～9.012 mm/a，形變速率較大區(qū)域位于坡體中部。由圖10（b）特征點(diǎn)與降雨數(shù)據(jù)構(gòu)建的時(shí)序形變曲線可知，特征點(diǎn)E 呈振蕩變化趨勢(shì)，雨旱兩季形變速率差異明顯。2020年8月后形變速率急劇增大，2021年4月之后形變速率相比同期增加約17 mm/a。

通過(guò)野外調(diào)查可知，H4 滑坡屬于臨江大型沖溝，溝面呈褶皺形態(tài)，目前尚未發(fā)育為真正意義的滑坡。圖10（b）時(shí)序形變曲線在2020年雨季后呈梯度下降趨勢(shì)，主要原因是降水沖刷溝壑表面使沖溝坡面向下滑動(dòng)。2021年4月之后相比同期形變速率明顯增加，此時(shí)受降雨影響微弱，說(shuō)明該滑坡體對(duì)水位變化有較強(qiáng)響應(yīng)，原本裸露的坡體下緣遭受江水侵蝕，下緣坡體在動(dòng)水壓力作用下土壤結(jié)構(gòu)趨向松散狀態(tài)，上緣沖溝體失穩(wěn)，自然產(chǎn)生向下形變。

圖10 H4 潛在滑坡形變特征Fig.10 H4 potential landslide deformation characteristics

H6 滑坡位于庫(kù)區(qū)西岸大灣子隧道，滑動(dòng)面處于隧道臨江一側(cè)，其形變速率如圖11所示，整體形變速率為-15.888～16.326 mm/a,選取滑坡體中部特征點(diǎn)F 與降雨數(shù)據(jù)建立時(shí)序形變曲線（圖11），特征點(diǎn)F 形變速率整體呈波動(dòng)趨勢(shì)，在2020年雨季形變速率較旱季增速明顯。2021年4月之后，形變速率緩慢增大，較同期增加約16 mm/a。

圖11 H6 潛在滑坡形變特征Fig.11 H6 potential landslide deformation characteristics

經(jīng)野外實(shí)地考察，發(fā)現(xiàn)H6 滑坡已發(fā)育且有部分滑動(dòng)痕跡，在坡體頂端還發(fā)育有一定程度的裂縫（圖11），所以在雨季降水沖刷坡面且沿裂縫滲入坡體改變其土體應(yīng)力結(jié)構(gòu)，使坡體產(chǎn)生較大形變。由圖11可以看出，2021年4—5月間，降雨量幾乎為零，但形變速率變化明顯，說(shuō)明該坡體對(duì)水位抬升具有較強(qiáng)響應(yīng)，降水沿裂縫進(jìn)入坡體內(nèi)部，促進(jìn)了破裂面的貫通，而水位抬升致使坡體下緣和滑動(dòng)面軟化，降低其抗剪強(qiáng)度，降雨和水位抬升的共同作用可能使H6 滑坡進(jìn)一步發(fā)育，后續(xù)應(yīng)當(dāng)對(duì)該滑坡進(jìn)行重點(diǎn)監(jiān)測(cè)。

5 結(jié)論

本文聯(lián)合升降軌Sentinel-1 SAR 數(shù)據(jù)，采用SBASInSAR 技術(shù)并結(jié)合無(wú)人機(jī)野外調(diào)查數(shù)據(jù)，分析白鶴灘水電站庫(kù)岸潛在滑坡的變形特征，得到以下結(jié)論：

（1）白鶴灘水電站庫(kù)區(qū)LOS 方向形變速率為-90.959～91.426 mm/a，受蓄水因素影響，各庫(kù)岸典型潛在滑坡形變速率明顯加快，蓄水前后形變平均增速達(dá)10 mm/a 以上；

（2）白鶴灘水電站庫(kù)岸潛在滑坡對(duì)水位變化具有較強(qiáng)響應(yīng)，蓄水量增加是當(dāng)前庫(kù)岸潛在滑坡發(fā)育的關(guān)鍵性誘因，水位抬升之后潛在滑坡形變速率變化明顯，在降雨和蓄水等因素共同作用下，白鶴灘水電站庫(kù)岸潛在滑坡存在失穩(wěn)風(fēng)險(xiǎn)；

（3）降軌數(shù)據(jù)集探測(cè)的形變信息較為豐富，主要集中在庫(kù)區(qū)西岸，而升軌數(shù)據(jù)集僅在庫(kù)區(qū)東岸部分區(qū)域形變較為明顯，故聯(lián)合升降軌SAR 數(shù)據(jù)能有效克服僅利用單一軌道導(dǎo)致的幾何畸變等問(wèn)題，使水電站庫(kù)岸潛在滑坡變形監(jiān)測(cè)更加準(zhǔn)確、全面。