0.   引言

库岸滑坡是深切割高山峡谷型库岸常见的破坏形式，多集中分布于我国西南山区^[1-2]。库岸滑坡由地表内外营力相互作用而形成，人类工程建设及库区水位变化则使其演化特征更为突出，通常表现为库区蓄水之后库岸下缘坡体岩层软化，从而引起上部库岸的形变破坏^[3-4]。库岸滑坡受多方面因素影响，具有成因复杂、类型多样和危害巨大等特点^[5-7]。在水电站库区，蓄水和泄洪等因素导致的水位变化直接影响库岸滑坡的稳定性，库岸滑坡一旦失稳会诱发一系列次生灾害，破坏区域生态系统，毁坏库区坝体和发电设施^[8]，严重威胁库区上下游居民生命财产安全。因此，对水电站库岸滑坡进行变形监测具有重要意义。

常规监测手段已难以识别和监测大面积的库岸滑坡形变，相对传统的精密水准测量、全球卫星导航系统（Global Navigation Satellite System, GNSS）和光学遥感技术，合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar，InSAR）技术因其大范围、全天时、全天候、高精度和高分辨率等特点，已成功应用于库岸滑坡灾害识别监测分析中^[9-11]。国内外学者利用InSAR技术在库岸滑坡灾害的应用方面做了大量研究，康亚等^[8]利用三类InSAR产品和DEM数据对金沙江流域乌东德水电站段进行滑坡早期识别，成功探测到多处未知和已知的滑坡体；徐帅等利用墨兰指数对SBAS-InSAR技术获取的形变点进行空间域异常值分析和聚类处理，成功识别出三峡库区巫山—奉节段高概率潜在滑坡范围^[12]；王振林等^[13]利用SBAS-InSAR技术提取雅砻江流域锦屏一级水电站库区左岸边坡的形变特征信息，推断出大幅水位上升是诱发滑坡复活的主要因素；朱同同等^[14]结合时序InSAR技术和GPS观测值分析了降雨和蓄水对三峡库区树坪滑坡变形的影响；史绪国等^[15]联合分布式目标与点目标的时序InSAR技术对三峡库区藕塘滑坡进行稳定性监测；Tantianuparp等^[16]联合多种SAR数据和PS-InSAR技术对三峡巴东进行滑坡探测，并将PS点时序形变与水位时间变化进行初步相关性分析；Zhou等^[17]利用时序InSAR技术发现三峡库区木鱼堡滑坡变形主要发生在水库涨落期和高水位期；Liu等^[18]利用SBAS-InSAR技术对三峡巴东地区进行滑坡探测并分析季节性滑坡运动与水位变化之间的相关性，上述研究证明了时序InSAR技术在库岸滑坡监测中的可靠性，可以对水电站库岸滑坡变形进行有效分析。白鹤滩水电站地处四川和云南交界，自2021年4月开始蓄水，库区水位由660 m升至825 m，上升幅度达165 m；水电站运行期间最低水位765 m，最高水位825 m，升降水位差60 m，最大库容达256×10⁸ m^3[19-20]。库区地形起伏较大、断裂构造发育，加之蓄水引起的水位变化直接影响库岸潜在滑坡的变形趋势，对水电站基础设施和上下游居民生命财产安全造成潜在威胁^[21-22]。因此，亟需对白鹤滩水电站库岸潜在滑坡进行变形分析。

文章联合2019年7月3日至2021年7月28日的150景升降轨Sentinel-1 SAR数据集，采用SBAS-InSAR技术获取白鹤滩水电站库区雷达视线方向（Line of sight，LOS）形变时间序列，在分析地表形变时间演化规律和空间分布特征的基础上，结合无人机野外调查，分析白鹤滩水电站库岸潜在滑坡的变形特征，重点研究蓄水因素对库岸潜在滑坡变形趋势的影响。

1.   SBAS-InSAR技术

小基线集InSAR（Small Baseline Subset InSAR，SBAS-InSAR）技术最早由Berardino和Lanari等^[23]提出，该方法通过组合数据的方式获得一系列短空间基线差分干涉图，这些差分干涉图能较好地克服空间失相关现象。SBAS-InSAR技术利用奇异值分解（SVD）法求解形变速率，将被较大空间基线分开的孤立SAR数据进行连接，进一步提高观测数据的时间采样率^[23]。该方法可以有效减弱大气效应，降低相位噪声和误差^[24]，其基本原理及流程如下：

假定已获取覆盖同一区域的按时间序列排序的$ N + 1 $幅SAR影像：

根据干涉组合规则，生成M幅干涉图且M应当满足：

假设以$ {t_0} $作为影像获取起始时刻且$ {t_0} $时刻影像覆盖区域位移为0，则在去除轨道误差、平地效应及地形相位的影响后，第$ i\left( {1 \leqslant i \leqslant M} \right) $幅影像某像素的干涉相位可表示为：

式中：$ \Delta {\varphi _{{i_{\rm{def}}}}} $——斜距向形变产生的相位；

$ \Delta {\varphi _{{i_{\rm{topo}}}}} $——地形相位；

$ \Delta {\varphi _{{i_{\rm{atm}}}}} $——大气延迟引起的相位；

$ \Delta {\varphi _{{i_{\rm{noise}}}}} $——相干噪声造成的相位。

利用最小二乘或者奇异值分解（SVD）对m个解缠相位进行三维时空相位解缠即可获得不同SAR时刻对应的时序形变速率。

2.   研究区概况和研究数据

本文以四川省与云南省交界白鹤滩水电站库区作为研究区域，如图1所示。研究区长约30.38 km，宽约11.65 km，总面积353.93 km²，地处横断山脉东北部、青藏高原东南边缘，区域内断裂构造发育，构造运动强烈，河谷深切，山体陡峻，地震频发^[25-27]。最高海拔3556 m，最低海拔520 m，高差达3036 m，地势陡峭，致使该区存在大量滑坡、崩塌和泥石流等地质灾害隐患。

图  1  研究区位置

Figure  1.  Location of study area

下载: 全尺寸图片 幻灯片

形变监测数据选用从欧州航天局(European Space Agency, ESA)免费下载的150景C波段Sentinel-1雷达影像（其中升轨数据50景，降轨数据100景并在每个时间点上下两景拼接），升降轨数据覆盖区域如图2所示。时间跨度为2019年7月3日至2021年7月28日，极化方式为VV，成像方式为IW，数据参数如表1所示。为提高影像轨道精度，引入POD精密定轨星历数据。使用日本宇宙航空研究开发机构（Japan Aerospace Exploration Agency，JAXA）发布的ALOS WORLD 3D 30 m空间分辨率的数字高程模型（Digital Elevation Model，DEM），用于去除地形相位影响，如图3所示。

图  2  SAR卫星影像覆盖范围

Figure  2.  SAR satellite image coverage

下载: 全尺寸图片 幻灯片

表  1  Sentinel-1A数据参数

Table  1.  Sentinel-1A data parameters

轨道方向	成像模式	极化方式	波长	波段	入射角/(°)
升轨	IW	VV	5.63	C	39.44
降轨	IW	VV	5.63	C	39.28

下载: 导出CSV 

| 显示表格

图  3  研究区DEM

Figure  3.  Digital elevation model of study area

下载: 全尺寸图片 幻灯片

3.   SBAS-InSAR技术数据处理

采用SBAS-InSAR技术，选取经镶嵌、配准和裁剪后的100景Sentinel-1A斜距单视复数（Single Look Complex，SLC）影像（升降轨数据各50景），根据时间基线和垂直基线最优原则，升轨和降轨数据分别以日期为20191216和20200204的影像作为超级主影像。设置时间基线阈值180d，空间基线为临界基线阈值的50%，共生成654和888对干涉像对。为抑制斑点噪声，设置多视数为1∶4，采用Minimum Cost Flow 解缠方法和Goldstein滤波方法进行干涉处理，将组合干涉对经过配准，调整删除不理想的数据后生成干涉图，研究区部分较理想的干涉图如图4所示。

图  4  研究区部分较理想的干涉图

注：（a）、（b）、（c）为升轨数据干涉图，（d）、（e）、（f）为降轨数据干涉图

Figure  4.  Ideal interference patterns in the study area

下载: 全尺寸图片 幻灯片

经过轨道精炼和重去平，利用最小二乘法和奇异值矩阵分解进行形变反演，然后估算和去除大气相位，得到研究区时间序列形变信息，对时序信息地理编码后获取研究区2019年7月3日至2021年7月28日LOS方向的形变结果。如图5所示，形变速率为正值表示靠近卫星，负值表示远离卫星。对比图5（a）、（b）研究区形变结果可知，降轨数据集探测的形变信息较为丰富，主要集中在库区西岸，最大LOS向形变速率−61.425 mm/a；升轨数据集仅在库区东岸部分区域形变较为明显，最大LOS向形变速率为91.426 mm/a。升降轨数据集形变信息不一致的原因是白鹤滩水电站库区两岸地形起伏较大，山势陡峭险峻，而升轨数据飞行方向大致沿东南向西北，雷达视线方向位于右侧，降轨数据则与之相反，故利用InSAR探测形变过程中阴影、叠掩和透视收缩等几何畸变现象严重。

图  5  研究区视线向形变速率

Figure  5.  Line-of-sight deformation rate of the study area

下载: 全尺寸图片 幻灯片

4.   试验结果与分析

4.1   研究区库岸典型潜在滑坡选取

对升轨和降轨数据获取的研究区形变结果进行综合解译，升轨数据库岸形变区域解译结果如图6所示，共选取库岸形变较大区域4处。结合无人机野外调查结果，发现典型潜在滑坡2处，分别用H1和H2表示；非滑坡形变区2处，分别用X1和X2表示，升轨数据详细解译结果如表2所示。

图  6  升轨潜在滑坡解译及实地考察结果

Figure  6.  Interpretation and field investigation results of potential landslide in ascending orbit

下载: 全尺寸图片 幻灯片

表  2  升轨数据库岸形变区域解译结果列表

Table  2.  List of interpretation results of shore deformation region in orbit lifting database

编号	形变区域名称	最大形变速率/（mm·a−1）	形变区域类别
H1	观音岩	19.846	潜在滑坡
H2	鱼坝	18.537	潜在滑坡
X1	六城村	76.259	非滑坡形变
X2	半坡	55.947	非滑坡形变

下载: 导出CSV 

| 显示表格

降轨数据库岸形变区域解译结果如图7所示，共选取库岸形变较大区域6处。结合无人机野外调查结果，发现典型潜在滑坡4处，分别用H3至H6编号；非滑坡形变区2处，分别用X3和X4表示，降轨数据详细解译结果如表3所示。

表  3  降轨数据库岸形变区域解译结果列表

Table  3.  List of interpretation results of shore deformation region in orbit descent database

编号	形变区域名称	最大形变速率/（mm·a−1）	形变区域类别
H3	观音岩	10.726	潜在滑坡
H4	清水沟	17.605	潜在滑坡
H5	鱼坝	19.326	潜在滑坡
H6	大湾子	15.888	潜在滑坡
X3	六城村	48.871	非滑坡形变
X4	半坡	61.425	非滑坡形变

下载: 导出CSV 

| 显示表格

图  7  降轨潜在滑坡解译及实地考察结果

Figure  7.  Interpretation and field investigation results of potential landslide in descending orbit

下载: 全尺寸图片 幻灯片

对比升轨和降轨数据解译结果可以看出，非滑坡形变区X1、X2与X3、X4分别相同，潜在滑坡H1、H2与H3、H5相互对应。另外，降轨数据还解译出除上述区域以外的潜在滑坡H4和H6，同一时间段不同轨道SAR数据集探测的形变结果能够相互对应，从侧面验证了本文InSAR结果的准确性，但受时间、空间失相干因素和几何畸变影响，升降轨形变信息有所差异，说明升降轨结合的方式能够有效弥补仅利用单一轨道识别结果不全面、不准确的缺陷，提升库岸潜在滑坡灾害识别和监测的准确性和有效性。

4.2   库岸典型潜在滑坡变形分析

结合4.1节升降轨数据集库岸潜在滑坡解译结果，本文选取H1、H2、H4和H6四处典型潜在滑坡进行变形分析，分别在各滑坡形变结果中选取特征点，引入研究区降雨数据，绘制特征点在蓄水前后的时序形变曲线，并结合无人机野外调查结果分析库岸典型潜在滑坡的变形特征。

H1滑坡地处观音岩，位于沿江公路东岸，滑坡形变速率如图8（a）所示。滑坡整体形变速率范围为−10.726～15.433 mm/a，分别选取滑坡体上缘和下缘特征点A、B与降雨数据构建时序形变曲线如图8（b）所示，特征点A和B时序形变速率波动趋势大致相同，每年雨季形变速率较旱季明显增加。2019年10月—2020年5月形变速率减小，2021年4月后形变速率增大，同比增加约16 mm/a。

图  8  H1潜在滑坡形变特征

Figure  8.  H1 potential landslide deformation characteristics

下载: 全尺寸图片 幻灯片

经实地勘察，该滑坡坡体上缘为自然坡体，坡体下缘已进行边坡加固，故在2019年10月至2020年5月间B点较A点形变速率变化相对稳定。受降雨因素影响，坡体在雨季滑动速率增大。2021年4月至5月，降雨量几乎为零，水电站蓄水导致库区水位上升，库岸下缘受到江水侵蚀改变坡体上下缘间的平衡关系，使该滑坡体形变量增大。

H2滑坡地处鱼坝村，金沙江支流末端。滑坡形变速率如图9（a）所示，形变较大值处于坡体中上部，形变范围在−19.326～8.254 mm/a。在坡体两侧分别选择特征点C、D结合降雨数据构建时序形变曲线见图9（b），特征点C和D形变速率变化趋势基本一致，与降雨数据呈现一定相关性，雨旱两季形变速率差异较小。2020年1月至10月间，形变速率逐渐减小，2021年1月后形变速率振荡变化，2021年4月之后，形变速率增加值超过10 mm/a。

图  9  H2潜在滑坡形变特征

Figure  9.  H2 potential landslide deformation characteristics

下载: 全尺寸图片 幻灯片

经野外实地调查，H2滑坡滑面自上而下呈倒“V”字形。由于隧道工程尚未完工，附近仍伴有部分工程活动，故在2020年雨季坡体滑动速率对降雨因素响应较弱，表现为形变速率逐渐减小。从图9（b）可以看出，2021年4月以后，特征点C和D形变速率相对之前有所增加，此时降雨量较小，说明蓄水导致的库区水位抬升也对远离河道的坡体产生影响。

H4滑坡体地处清水沟，位于库区西岸。滑坡形变速率如图10（a）所示，滑坡整体形变范围为−17.605～9.012 mm/a，形变速率较大区域位于坡体中部。由图10（b）特征点与降雨数据构建的时序形变曲线可知，特征点E呈振荡变化趋势，雨旱两季形变速率差异明显。2020年8月后形变速率急剧增大，2021年4月之后形变速率相比同期增加约17 mm/a。

图  10  H4潜在滑坡形变特征

Figure  10.  H4 potential landslide deformation characteristics

下载: 全尺寸图片 幻灯片

通过野外调查可知，H4滑坡属于临江大型冲沟，沟面呈褶皱形态，目前尚未发育为真正意义的滑坡。图10（b）时序形变曲线在2020年雨季后呈梯度下降趋势，主要原因是降水冲刷沟壑表面使冲沟坡面向下滑动。2021年4月之后相比同期形变速率明显增加，此时受降雨影响微弱，说明该滑坡体对水位变化有较强响应，原本裸露的坡体下缘遭受江水侵蚀，下缘坡体在动水压力作用下土壤结构趋向松散状态，上缘冲沟体失稳，自然产生向下形变。

H6滑坡位于库区西岸大湾子隧道，滑动面处于隧道临江一侧，其形变速率如图11所示，整体形变速率为−15.888～16.326 mm/a, 选取滑坡体中部特征点F与降雨数据建立时序形变曲线（图11），特征点F形变速率整体呈波动趋势，在2020年雨季形变速率较旱季增速明显。2021年4月之后，形变速率缓慢增大，较同期增加约16 mm/a。

图  11  H6潜在滑坡形变特征

Figure  11.  H6 potential landslide deformation characteristics

下载: 全尺寸图片 幻灯片

经野外实地考察，发现H6滑坡已发育且有部分滑动痕迹，在坡体顶端还发育有一定程度的裂缝（图11），所以在雨季降水冲刷坡面且沿裂缝渗入坡体改变其土体应力结构，使坡体产生较大形变。由图11可以看出，2021年4—5月间，降雨量几乎为零，但形变速率变化明显，说明该坡体对水位抬升具有较强响应，降水沿裂缝进入坡体内部，促进了破裂面的贯通，而水位抬升致使坡体下缘和滑动面软化，降低其抗剪强度，降雨和水位抬升的共同作用可能使H6滑坡进一步发育，后续应当对该滑坡进行重点监测。

5.   结论

本文联合升降轨Sentinel-1 SAR数据，采用SBAS-InSAR技术并结合无人机野外调查数据，分析白鹤滩水电站库岸潜在滑坡的变形特征，得到以下结论：

（1）白鹤滩水电站库区LOS方向形变速率为−90.959～91.426 mm/a，受蓄水因素影响，各库岸典型潜在滑坡形变速率明显加快，蓄水前后形变平均增速达10 mm/a以上；

（2）白鹤滩水电站库岸潜在滑坡对水位变化具有较强响应，蓄水量增加是当前库岸潜在滑坡发育的关键性诱因，水位抬升之后潜在滑坡形变速率变化明显，在降雨和蓄水等因素共同作用下，白鹤滩水电站库岸潜在滑坡存在失稳风险；

（3）降轨数据集探测的形变信息较为丰富，主要集中在库区西岸，而升轨数据集仅在库区东岸部分区域形变较为明显，故联合升降轨SAR数据能有效克服仅利用单一轨道导致的几何畸变等问题，使水电站库岸潜在滑坡变形监测更加准确、全面。