中国现有各型水库9.8万余座（不含港、澳、台地区），总库容超过9 300亿m^3[1]。大多数水库兴建于20世纪50~70年代，受限于当时的经济、技术条件以及数十年间工程老化等问题的影响，水库病险情况曾普遍存在，最多时有3万余座水库处于病险状态。病险水库除险加固任务具有艰巨性、复杂性和长期性的特点，根据专项规划和治理，全国范围内将逐步消除病险水库。

形变监测是水库大坝安全评价的基础和重要组成。现阶段水库形变监测多是在重点部位布设分散、不连续的单个监测点，利用全站仪、全球导航卫星系统等测量手段按以点代面、以局部代替整体的方法采集离散点形变信息进行数据处理，这种方式难免会遗漏一些重大的安全隐患。遥感技术的飞速发展为水利工作者提供了新型工具和技术手段。星载合成孔径雷达干涉测量（interferometric synthetic aperture radar, InSAR）作为近年来迅速发展的空间大地测量新技术，能够全天时、全天候、大面积同步获取地形信息及地表形变信息，与传统测量手段相比，有效缩短了观测周期，扩大了空间范围，提高了监测精度^[2-4]。

国内外已有学者利用InSAR技术进行水利工程形变监测的探索。廖明生等^[5]基于永久散射体（persistent scatterer, PS）干涉测量的时序InSAR技术对三峡大坝和三峡库区滑坡体进行了稳定性监测，详尽分析了三峡大坝坝体在建设过程中及工程完工后蓄水期的变形，其监测结果与实际情况相吻合。Voege等^[6]利用升降轨数据对挪威斯瓦特湖坝开展了基于小基线集技术的形变监测，结果显示，时序InSAR可以探测到坝体的局部变形情况。裴媛媛等^[7]利用时序合成孔径雷达（synthetic aperture radar, SAR）影像提取了上海长江口南岸和杭州湾北侧堤坝的沉降速率。Milillo等^[8]利用多时相InSAR技术对伊拉克摩苏尔大坝的持续失稳情况进行了调查，指出2004—2010年间大坝变形速率较大，2012—2014年间失稳情况得到控制，但自2014年8月大坝停止注浆作业后，形变再次加速。现有研究多集中于对超大型水利设施、枢纽工程的监测，将星载时序InSAR作为地面常规监测手段的补充，而对其他大中型水库大坝缺乏关注。从方法上而言，大多数研究是将时序InSAR方法直接应用于水库库区形变监测，并未考虑环境的复杂性和独特性，一定程度上限制了其发展和应用。本文综合考虑水库大坝特殊的环境特征和自身特点，研究适合实际需求的时序InSAR分析方法，并以山东省广南水库为研究区开展实验研究，证实了方法的有效性。

1 面向水库研究区时序InSAR方法

相比于城市地面沉降监测等应用领域，水库特殊的环境特征给时序InSAR数据处理造成了一定困难，主要表现在：（1）水利工程施工、运营周期长，且库区地表多生长茂密植被，面临因时空失相干严重而导致SAR影像无法精确配准的问题；（2）水库研究区常伴随大面积水体覆盖，SAR系统侧视成像水域范围内后向散射强度低，经典PS选取方法得到的空间分布较为稀疏，而简单放宽选取条件可能引起大量误报；（3）经典PS算法中非缠绕相位的恢复本质上是基于先验形变模型的时间维一维解缠方法，由于库区地质条件复杂，难以预先建立形变模型，形变参数解算存在困难。

针对上述现实问题，本文对常规时序InSAR数据处理流程进行改进，通过应用适当的策略，包括SAR影像传递配准方法、基于自适应阈值分割的PS选取方法以及时空维相位解缠策略等，满足水库复杂环境下时空连续、高精度形变监测要求。改进后的数据处理流程如图 1所示。

图  1  面向水库研究区的时序InSAR数据处理流程

Figure  1.  Flowchart of Time Series InSAR Data Processing Used for Reservoirs Study Area

下载: 全尺寸图片 幻灯片

1.1   SAR影像传递配准方法

SAR影像精配准依赖于影像间的相干关系^[9]，若需要配准的两景影像因时间间隔长、基线条件不良等因素导致失相干严重，则常规配准方法可能导致配准失败。相对于常规配准方法，选取1景或多景与两幅待配准影像都具有较好相干性的辅助影像（集），通过几何关系和数学模型进行传递配准，可以从一定程度上解决这类问题。选择高质量的辅助影像成为传递配准法的关键。

以距离向为例，传递配准法所涉及的主影像、副影像p及辅助影像q之间的空间关系为：

式中，α和β为多项式系数；(R，A)为配准联系点在影像距离向和方位向的坐标；对于配准多项式阶次d，所需求解的未知数有(d + 1)(d + 2)个。当d = 2时，式（1）展开为：

引入辅助影像的传递配准方法在配准点及联系点选取、互相关系数计算、偏移量平差及影像重采样等方面均与常规配准方法相同^[10]。由于加入了辅助影像，需要确定的未知参数个数增加，但相比于精配准过程中数百个配准窗口上的亚像素级偏移量，仍有大量多余观测可供平差和多项式拟合。传递配准方法相对于常规配准法有更好的适用性，尤其适合低相干SAR影像对的配准。

1.2   基于自适应阈值分割的PS选取方法

基于PS的时序InSAR方法首先从SAR影像中提取具有稳定反射特性的高质量散射体，然后仅对具有高相干特性的像元进行差分干涉、相位建模和形变参数解算，最终获得高可靠性的地表微小形变信息^[11]。如何在影像上识别稳健的PS目标是时序InSAR技术的基础和关键。水库库区地物类型复杂，常伴随大面积水体覆盖，既要顾及PS的探测概率，即让尽量多的真实PS被成功选取；又要保障PS探测结果的可靠性，即选取结果中只能包含少量失相干目标。针对这一实际问题，基于SAR振幅图像统计特征，提出利用自适应阈值分割结果进行约束，综合分析振幅及干涉相位信息，以得到足够密度的高质量PS目标。

二值分割将图像分割为背景和目标两类，此时图像灰度直方图可视为组合背景和目标像素的混合分布。由SAR影像的成像特征^[12]及中心极限定律可知，该混合分布的两个分量均服从高斯分布。给定判断标准：

式中，μ和σ分别为均值与方差，下标0和1分别代表目标和背景；参数f需根据灰度分布特征确定。

对于输入的时域平均SAR振幅图像，首先根据其灰度直方图概率密度分布确定初始阈值，随后按照一定的规则对图像分块，计算块内灰度值的方差，并将其与方差阈值进行比较。若块内像素灰度值方差较小，可认为该子块内所有像素为同一地物类型，对比块内灰度均值与图像总体灰度均值，以判定该子块为陆地或水体；对于非同类子块，利用Otsu方法^[13]选取局部图像分割阈值，迭代处理直至满足式(3)的条件。通过上述自适应阈值分割，可得到水体、陆地两类分割结果。

基于自适应阈值分割的水陆分离结果，仅对SAR图像上被分类为陆地的部分利用经典振幅离差指数阈值法^[11]和时域相干性测度法^[14]进行稳定性分析，可有效避免水域区域PS的误判。

1.3   时空维相位解缠

经典PS算法是利用干涉图时间序列上两点间的缠绕相位差估计形变参数^[15]，这种相位恢复的思路是一种基于模型的时间维一维解缠策略。水库库区涉及范围广，形变情况复杂，难以预先建立形变模型。对干涉对进行逐一解缠，再将结果统一归算到时间维上的方法是一种简化处理，但时序InSAR相位是时域上一维、空间域上二维的三维数据，充分利用其三维特性进行相位解缠将取得更高精度的参数估计结果。Hooper等^[16]提出了三维解缠理论框架，给出了无多环路相位中断平面的理想情况下似L^∞范数解。

现有的解缠算法均基于如下假设，即邻近像元的相位差绝对值小于π。干涉相位组成部分中，地形相位、形变相位及大气相位在空间上具有相似性，对于这部分空间相关的信号而言，当空间采样率足够高时，上述假设前提是成立的；但信号中存在的空间不相关的部分，即使空间采样间隔足够小，仍有可能存在邻近像元绝对相位差大于π的情况。空间不相关的信号中，影响最大的成分是视角误差，主影像中空间非相关的误差部分也将被估计。从缠绕的干涉相位φ_int中去除空间不相关的部分，仅对剩余部分进行解缠，即：

式中，W {∙}为缠绕算子；$\text{ }\!\!\Delta\!\!\text{ }\hat{\phi }_{\theta }^{\text{nc}}$和$\phi _{m}^{\text{nc}}$分别表示空间不相关的视角误差和主影像误差；；ϕ_def为卫星视线方向上的形变；ϕ_atm表示两次SAR成像时刻的大气相位差异；Δϕ_orb是由轨道不确定性造成的相位差；$\text{ }\!\!\Delta\!\!\text{ }\phi _{\theta }^{\text{corr}}$为视角误差中空间相关部分造成的地形相位残余；；Δϕ_noise为对应的噪声项。

通过二维空间滤波从干涉相位中去除视角误差和主影像误差中空间不相关的部分后，在InSAR时间序列的时间维上进行相位恢复，随后利用该结果作为空间上二维解缠的优化初值。这不是严格意义上的三维解缠，但相比于三维解缠理论框架中的似L^∞范数解，该方法更具有普适性，尤其是数据中存在多环路相位中断平面的情况。二维相位解缠有最小费用流、最小统计费用流^[17]等诸多方法。

2 研究区域概况及SAR数据情况

山东省东营市地处黄河最下游，境内客水资源较丰沛，但黄河来水在年际、年内分配不均，无法满足生产生活对水资源的需求。修建平原水库，利用黄河水调剂当地水资源，做到“丰蓄枯用”、“冬蓄春用”。广南水库位于东营市东南，濒临渤海，在支脉沟以北、广利河以南、防潮大堤以西退海滩涂的开阔地带。为优化胶东半岛水资源配置，缓解潍坊、青岛地区的水资源供需矛盾，于2017年开工的黄水东调应急工程将从广南水库向潍坊北部地区调水。

广南水库围坝为碾压式均质土坝，坝顶高程7.0 m，平均库底高程2.0 m，设计水位5.0 m，占地39 km²，库容1.14亿m³，属大（2）型水库。本研究选取30景ENVISAT ASAR影像为数据源，时空基线详细信息见表 1。图 2显示研究区地理位置，其中黑色矩形框是该组影像整景覆盖范围，蓝色矩形对应广南水库研究区，截取的区域在影像坐标系下为1 000×1 000像素，底图显示1″空间分辨率的SRTM（shuttle radar topography mission）高程信息，该区域平均高程不超过5 m。

表  1  用于时序InSAR分析的ENVISAT ASAR影像时空基线表

Table  1.  Detailed Information of Temporal and Spatial Baselines of ENVISAT ASAR Images Used for Time Series InSAR Analysis

影像获取时间	时间基线/d	垂直基线/m
2007-02-01	-805	-161
2007-03-08	-770	687
2007-04-12	-735	36
2007-05-17	-700	73
2007-06-21	-665	276
2007-07-26	-630	207
2007-10-04	-560	115
2007-11-08	-525	408
2008-01-17	-455	166
2008-02-21	-420	-114
2008-03-27	-385	551
2008-10-23	-175	451
2008-11-27	-140	-48
2009-02-05	-70	-81
2009-03-12	-35	725
2009-04-16	0	0
2009-05-21	35	186
2009-06-25	70	549
2009-07-30	105	201
2009-09-03	140	517
2009-10-08	175	111
2009-11-12	210	496
2009-12-17	245	-97
2010-01-21	280	474
2010-02-25	315	-97
2010-04-01	350	638
2010-05-06	385	453
2010-06-10	420	458
2010-07-15	455	143
2010-09-23	525	457

下载: 导出CSV 

| 显示表格

图  2  研究区地理位置示意图

Figure  2.  Sketch Map of Geographical Location of the Study Area

下载: 全尺寸图片 幻灯片

3 水库大坝时序InSAR监测结果分析

广南水库围坝内坡为现浇混凝土板护坡，外坡大部分为草皮护坡。由于东坝外坡是临海面，为防止外坡受潮汐的冲击侵蚀，在外坡坡脚以上至4.0 m高程范围内采用浆砌石护坡和混凝土护坡。坝体大部分为混凝土浇筑，具有较好的电磁波反射特性。莱州湾是风暴潮频发区，为防止潮浪侵袭，莱州湾南岸粉砂淤泥质海岸筑有人工土石海堤。防潮堤、广南水库库岸护坡及湖心岛等人工地物上均提取到较多高质量的PS目标，居民区和企业分布的地区也探测到一定数目的稳定点，水产养殖池堤等位置由于变化较大，仅在部分区域获取到零星点。广南水库西侧为内陆侧背海面，护坡为混凝土浇筑，受潮汐影响较小，选取分布于该护坡上77个PS的平均相位作为基准，时间维以第1景影像的获取时间（2007-02-01）为计算基准，经时空维相位解缠、数字高程模型（digital elevation model，DEM)及其他相关误差项的估计，得到各PS上对应的相对变化情况，卫星视线方向形变时间序列如图 3所示。

图  3  广南水库库区卫星视线方向形变时间序列

Figure  3.  Time Series of Deformation of Guangnan Reservoir Area in Line of Sight Direction

下载: 全尺寸图片 幻灯片

由图 3可见，广南水库西岸和湖心长廊较为稳定，无显著变形现象；莱州湾防潮堤受潮汐冲击侵蚀，水库东岸受海水内侵影响，在研究期内有形变发生，但总体情况可控；其余地区如水产养殖场和盐池等局部区域，由于受工农业生产活动影响，其变化较大。假设研究期间为均匀线性形变，各PS平均形变速率如图 4(a)所示，各点位线性形变速率标准差如图 4(b)所示，其中底图为SAR平均强度。将PS点及其平均形变速率信息叠加到GoogleEarth^TM上，以识别所对应的地物类型，如图 5所示。

图  4  广南水库库区平均形变速率与标准差分布图

Figure  4.  Distribution Maps of Averaged Deformation Rates and Standard Deviations of the Guangnan Reservoir Area

下载: 全尺寸图片 幻灯片

图  5  广南水库库区卫星视线方向平均形变速率图

Figure  5.  Averaged Deformation Rates of the Guangnan Reservoir Area in Line of Sight Direction with Optical Image Background

下载: 全尺寸图片 幻灯片

在探测到的所有PS点中，线性形变速率拟合标准差最大为3.7 mm/a，最小为0.3 mm/a，其中有超过91%的点位标准差小于2 mm/a。堤坝、库岸等混凝土浇筑的人工地物由于相干条件较好，反映在图 4(b)中表现为标准差较小，优于1 mm/a。由低洼闲置田改建的淡水养殖池塘等地物本身相干条件较差，且受工农业生产活动影响，形变情况并不规律，因此线性形变速率拟合的标准差较大。选取4类典型地物，即莱州湾防潮堤D、广南水库湖心人工长廊L、迎海面受潮汐影响的东侧库岸E和内陆侧不受潮汐影响的西侧库岸W（标识于图 4(a)中），分析其形变时间序列，如图 6所示。图 6中直线是利用线性函数拟合得到的形变趋势，形变参考时间为第1景影像的获取时间（2007-02-01）。

图  6  典型特征区域PS目标点形变时间序列

Figure  6.  Time Series of Deformation of PS Target Points in Typical Feature Areas

下载: 全尺寸图片 幻灯片

图 6显示，莱州湾防潮堤上典型点位D的形变速率约为5 mm/a，形变发生的主要原因是受海潮冲刷、侵蚀等外界作用影响。同样，受到潮汐影响的还有广南水库的东侧迎海面库岸，典型点位E上的平均形变速率约为3 mm/a，相比于防潮堤，东侧库岸受潮汐侵蚀作用稍小，因此形变速率也较之有所下降。广南水库内陆侧西侧库岸和湖心人工长廊由于不受海潮的影响，在研究期间均保持稳定，考虑到干涉测量误差影响，可认为西侧库岸和湖心人工建筑物无显著形变发生。

综上所述，莱州湾防潮堤、广南水库库岸及景区湖心人工建筑物由混凝土浇筑，研究期间干涉条件较好，可以从相关人工地物上提取到足量的高质量PS进行时间序列分析，获得有效的形变信息。

4 结语

表面形变监测结果是水库大坝风险分析与预警的基础，能够为防灾减灾提供决策支持。星载时序InSAR提供了一种遥感新方法，可以较低的成本进行大面积同步观测，获取高空间分辨率、高精度的形变信息。结合大坝、堤坝及库区实际应用环境的复杂性和特殊性，本文对时序InSAR方法的数据处理流程进行了相应改进，使其更加适应水电工程形变监测实际需求。广南水库应用实践表明，由混凝土浇筑的大坝、防潮堤等水工建筑物干涉条件较好，时序InSAR方法可探测到高质量的散射体目标，并提取形变时间序列，验证了时序InSAR分析方法对水库库区、大坝、防潮堤等水工建筑物进行形变监测的有效性。时序InSAR形变监测方法在水库安全状况普查及形变历史回溯等工作中具有巨大的应用潜力。