变形控制是高堆石坝工程设计、施工、科研和管理各方共同关注的问题^[1-2]. 变形监测作为变形控制的重要环节，是掌握大坝工作性态、评价大坝运行状况的重要措施，在工程建设和运行管理的全生命周期中不可或缺. 要准确地评估堆石坝的工作性态须掌握坝体应力变形的时空分布和演化，目前常用的方法是有限元数值模拟^[3-5]. 常用的堆石体本构模型存在缺陷和参数取值问题，导致有限元数值模拟结果与实际情况存在明显差异^[6-8]. 如何有效地利用变形监测数据重构堆石坝的变形场以进行工作性态评估，是堆石坝变形控制的关键环节之一.

为了掌握堆石坝的工作性态，须对其内部和外观进行定期的变形监测，其形变增量一般在毫米级，因此对监测精度要求较高. 常规监测技术，如水管式沉降仪、引张线式水平位移计、水准测量等，具有监测精度高、结果可靠、使用经验丰富等优点，但其只能在重点部位布设分散的监测点，按照以点代面、以局部代整体的方法采集离散的变形数据. 合成孔径雷达干涉（interferometric synthetic aperture radar，InSAR）技术^[9-11]具有全天候、全天时、空间分辨率高、不受恶劣天气影响等优点，极大地弥补了常规监测技术的不足. 然而，InSAR技术的观测结果为地表真实形变在雷达视线向的投影，难以反映地物的真实变形特征及演化规律^[12]；InSAR技术易受大气和时空失相关的影响，导致其单点监测精度略低于常规监测技术的. 因此，将新型与常规监测技术进行优势互补，可能是解决堆石坝变形监测难题的可行方案之一. 但由于新型监测技术与常规监测技术的观测机制迥异，多源监测数据的时空分辨率、精度差异较大，传统的数据处理方法在处理该类数据时有较大的局限性，目前尚缺乏合理的多源数据融合方法^[13].

本研究以水布垭面板堆石坝为例，利用InSAR技术获取水布垭面板堆石坝在2007年2月至2009年1月的视线向与方位向的大范围面状变形，充分考虑InSAR技术与水准监测的特点，基于集合卡尔曼滤波（ensemble Kalman filter，EnKF）将两者进行数据融合，实现优势互补；通过将多维度数据联合解算，重构水布垭面板堆石坝的外观变形场，以全面掌握堆石坝的整体变形性态.

与常规监测技术以点代面不同，InSAR技术具有空间连续的优势，时空连续的面状观测方式极大地丰富了堆石坝的外观变形数据. 不过，InSAR技术容易受到大气和时空失相关的影响，并且InSAR观测噪声在时间域上具有高斯性^[14]. 为了提高InSAR观测的精度，同时弥补常规监测技术空间离散的不足，采用集合卡尔曼滤波算法将分解到垂直向的InSAR观测结果与水准测量结果进行融合，实现2种监测技术的优势互补.

数据融合的概念于20世纪70年代被提出^[15]，后快速发展，被广泛应用到军事、工业、商业决策等领域，迄今仍是很多领域的研究热点，比如人工智能和大数据领域. 因其涵盖内容广泛、界限确定较困难，不同学者对数据融合有不同见解. 针对堆石坝变形监测的实际需求，本研究对数据融合的定义如下：将多层次观测体系获取的多源数据通过数据同化算法实现智能融合，获取比单一监测方式更高质量的结果（例如精度和分辨率），融合过程强调数据间的互补性，以实现对堆石坝变形性态的准确把握.

集合卡尔曼滤波算法以贝叶斯理论为框架，利用集合的思想求解误差方差矩阵，能够解决复杂的非线性高斯问题，对非线性较强、不连续的动态模型也能得到较好的同化结果，是数据同化领域的经典算法^[16]. 集合卡尔曼滤波根据背景场以及观测场集合的误差特征分布对其施加扰动，通过统计集合扰动场得到背景误差协方差矩阵，进而得到分析方程的增益矩阵，求取融合数据分析场的同化值以及相应的误差协方差，通过同化值实现多源数据的融合. 将通过不同监测手段获取的变形数据分别作为背景场数据与观测场数据，对背景场与观测场施加高斯扰动，得到背景场与观测场集合：

(1) $ {{\boldsymbol{X}}^{\text{b}}} = {{M}} {\boldsymbol{x}} + {\boldsymbol{q}};\;{\boldsymbol{q}}\sim N\left( {{\boldsymbol{0}},{\boldsymbol{Q}}} \right) \text{，} $

(2) $ {{\boldsymbol{Y}}^0} = {{H}} {\boldsymbol{y}} + {\boldsymbol{r}};\;{\boldsymbol{r}}\sim N\left( {{\boldsymbol{0}},{\boldsymbol{R}}} \right) . $

式中： ${{\boldsymbol{X}}^{\text{b}}}$、 ${{\boldsymbol{Y}}^0}$为背景场、观测场的状态量，x、y为通过不同监测手段获取的变形监测值，M为模型算子，H为观测算子； $ {\boldsymbol{q}} $、 $ {\boldsymbol{r}} $为数学期望为0，方差为Q、R的高斯白噪声；Q、R为背景场与观测场的误差协方差.

计算增益矩阵系数：

(3) $ {\boldsymbol{K}} = \frac{{{{\boldsymbol{P}}^{\text{b}}}{{\boldsymbol{H}}^{\text{T}}}}}{{{\boldsymbol{H}}{{\boldsymbol{P}}^{\text{b}}}{{\boldsymbol{H}}^{\text{T}}} + {\boldsymbol{R}}}} . $

式中： $ {\boldsymbol{P}}^{\mathrm{b}} $为样本数据的背景场误差协方差.

(4) $ {\boldsymbol{P}}_i^{\text{b}}{\text{ = }}\frac{1}{{n - 1}}\sum\limits_{i = 1}^n {\left( {{\boldsymbol{x}}_i^{\text{b}} - \overline {{\boldsymbol{x}}_i^{\text{b}}} } \right)} {\left( {{\boldsymbol{x}}_i^{\text{b}} - \overline {{\boldsymbol{x}}_i^{\text{b}}} } \right)^{\text{T}}}, $

(5) $ \overline {{\boldsymbol{x}}_i^{\text{b}}} = \frac{1}{{n - 1}}{\boldsymbol{x}}_i^{\text{b}}. $

式中： ${\boldsymbol{x}}_i^{\text{b}}$为背景场状态变量集合中第i个样本数据的预报值， $\overline {{\boldsymbol{x}}_i^{\text{b}}} $为 ${\boldsymbol{x}}_i^{\text{b}}$的平均值，n为样本数.

融合数据的分析场以及相应的误差协方差表达式如下：

(6) $ {\boldsymbol{x}}_i^{\text{a}} = {\boldsymbol{x}}_i^{\text{b}} + {\boldsymbol{K}}\left( {{\boldsymbol{Y}}_i^0 - {\boldsymbol{Hx}}_i^{\text{b}}} \right) \text{，} $

(7) $ {\boldsymbol{P}}_i^{\text{a}} = \frac{1}{{n - 1}}\sum\limits_{i = 1}^n {\left( {{\boldsymbol{x}}_i^{\text{a}} - \overline {{\boldsymbol{x}}_i^{\text{a}}} } \right)} {\left( {{\boldsymbol{x}}_i^{\text{a}} - \overline {{\boldsymbol{x}}_i^{\text{a}}} } \right)^{\text{T}}}. $

式中： ${\boldsymbol{x}}_i^{\text{a}}$为分析场状态变量集合中第i个样本数据的分析同化值即最优变形量， $\overline {{\boldsymbol{x}}_i^{\text{a}}} $为 ${\boldsymbol{x}}_i^{\text{a}}$的平均值.

(8) $ \overline {{\boldsymbol{x}}_i^{\text{a}}} = \frac{1}{{n - 1}}{\boldsymbol{x}}_i^{\text{a}}. $

通过误差分析，如果得到的同化值不满足预计精度，则对同化系统的参数组合同时进行调整与寻优，重新进行数据同化过程，直到符合精度要求并输出分析同化值以完成此轮数据融合工作，再进行下一时刻的融合计算，最终获得优势互补、高时间分辨率与高空间分辨率有机统一的多源变形监测数据融合结果，其技术路线图如图1所示.

要重构堆石坝的外观变形场，至少需要3个不同方向的变形数据或等效的先验信息^[17-18]. 视线向的时间序列InSAR观测与垂直向的水准测量提供了二维的变形信息. 多孔径InSAR（multiple aperture InSAR，MAI）技术显著提高了InSAR技术观测地表方位向变形的精度^[19]，为重构堆石坝外观变形场提供了方位向大范围面状变形信息. 多孔径InSAR技术通过方位向公共频谱滤波将干涉对中2幅合成孔径雷达（synthetic aperture radar，SAR）影像分别拆分为前视和后视影像，基于前、后视影像分别生成前视和后视干涉图，通过差分前、后视干涉图获得地表的方位向形变. 与Offset-tracking技术利用SAR影像的幅度信息不同，多孔径InSAR技术基于干涉相位信息，在高相干区域可以获取比Offset-tracking精度更高的方位向形变^[20]. 利用3个不同方向的变形监测结果重构堆石坝的外观变形场，其技术路线如图2所示.

如图3所示为SAR成像几何示意图，SAR卫星飞行方向称为方位向，雷达视线向与方位向相互垂直，地距向为视线向在水平面的投影方向. 地面任意一点的视线向变形矢量可以分解为垂直向和地距向的变形矢量，同时在水平面视角中，地距向变形矢量由南北向和东西向的变形矢量组成，由此可以建立雷达视线向变形D_LOS和地表三维变形的关系^[12]：

(9) $\begin{split} {{\boldsymbol{D}}_{{\rm{LOS}}}} = \;&{{\boldsymbol{D}}_{\rm{U}}} \cos \;\theta - {{\boldsymbol{D}}_{\rm{E}}} \sin \;\theta \sin \;\left(\alpha - {{3{\text{π}} }}/{2}\right) -\\ \;&{{\boldsymbol{D}}_{\rm{N}}} \sin\; \theta \cos \;\left(\alpha - {{3{\text{π}}}}/{2}\right). \end{split} $

式中：D_U、D_E、D_N分别为垂直向、东西向和南北向的变形矢量； $ \theta $为雷达入射角； $ \alpha $为卫星飞行方位角，即正北方向和卫星飞行方向的顺时针夹角； $\alpha - {{3{\text{π}} }}/{2}$为正北方向与地距向的夹角. 对于本研究选用的升轨ALOS PALSAR数据, 入射角 $ \theta $≈38°，方位角 $\alpha $≈349°.

对该地面而言，多孔径InSAR技术获取的方位向变形D_AZI位于由东西向和南北向组成的二维平面中，根据几何关系可以建立方位向变形和地表三维变形的关系：

(10) $ {{\boldsymbol{D}}_{{\rm{AZI}}}} = - {{\boldsymbol{D}}_{\rm{E}}}\cos \;\left(\alpha - {{3 {\text{π}} }}/{2}\right) + {{\boldsymbol{D}}_{\rm{N}}} \sin \;\left(\alpha - {{3{\text{π}}}}/{2}\right) . $

在垂直向、东西向、南北向的三维变形中，垂直向的变形由分解到垂直方向的InSAR测量结果和水准测量数据进行数据融合后的结果表示：

(11) $ {{\boldsymbol{D}}_{\rm{U}}} = {{\boldsymbol{D}}_{\rm{{EnKF}}}} . $

将式（9）、（10）写成矩阵形式，利用联合解算模型求解东西向和南北向分量的变形：

(12) $ {\boldsymbol{L}}{\text{ = }}{\boldsymbol{BX}}, $

(13) $ {\boldsymbol{B}}\text=\left[\begin{array}{cc}-\mathrm{sin}\;\theta \;\mathrm{sin}\;\left(\alpha -{3{\text{π}} }/{2}\right)& -\mathrm{sin}\;\theta \;\mathrm{cos}\;\left(\alpha -{3{\text{π}} }/{2}\right)\\ -\mathrm{cos}\;\left(\alpha -{3{\text{π}} }/{2}\right)& \mathrm{sin}\;\left(\alpha -{3{\text{π}} }/{2}\right)\end{array}\right] . $

式中： ${\boldsymbol{L}} = {\left[ {{{\boldsymbol{D}}_{{\rm{LOS}}}}},\;{{{\boldsymbol{D}}_{{\rm{AZI}}}}} \right]^{\text{T}}}$， ${\boldsymbol{X}} = {\left[ {{{\boldsymbol{D}}_{\rm{E}}}},\;{{{\boldsymbol{D}}_{\rm{N}}}} \right]^{\text{T}}}$，B为系数矩阵.

水布垭面板堆石坝位于清江中游河段、湖北省恩施土家族苗族自治州巴东县境内，上距恩施市117 km，下距清江第2梯级隔河岩电站92 km，坝顶高程409 m，坝轴线长度608 m，最大坝高233.2 m，是世界最高的混凝土面板堆石坝. 水库正常蓄水位400 m，相应库容43.12亿m³，校核洪水位高程404 m，总库容45.8亿m³. 工程总工期7年9个月，2000年前为筹建期，2001年1月—2002年10月为施工准备期，2002年11月—2007年6月为主体工程施工期，2007年7月—2008年9月为工程完建期. 2006年10月大坝坝体填筑至高程405 m，并通过蓄水验收，2007年4月水库开始集中蓄水，2007年7月首台机组发电，2008年9月4台机组全部并网发电，2008年11月2日水库水位达到399.51 m. 本研究选取2007年2月—2009年1月的监测数据进行数据融合. 该时间段正值水布垭面板堆石坝的初期蓄水阶段，坝体变形速率较大，更有利于展示数据融合方法的优势.

水布垭面板堆石坝的外观变形监测主要由大坝下游面布设的11个观测房组成，分别位于235 m（WS01）、265 m（WS02）、300 m（WS03、WS04、WS05）、335 m（WS06、WS07、WS08）和370 m（WS09、WS10、WS11）高程处，如图4所示. 图中，d为沉降，H为水位。11个监测点处（WS01~WS11）地面水准测量的沉降历时曲线表明，在蓄水阶段，随着库水位上升，测点沉降明显增大. 在约6个月后，各个测点的变形速率逐渐减缓并趋于稳定.

SAR影像选用覆盖水布垭地区的11景L波段ALOS PALSAR升轨影像（2007年2月28日—2009年1月18日），基本参数如表1所示. 此外，选用覆盖研究区域的美国宇航局航天飞机雷达地形测量（shuttle radar topography mission，SRTM）数据作为外部数字高程模型（digital elevation model，DEM）以去除干涉图的地形相位.

时间序列InSAR技术根据主影像的数量，可以大致分为以永久散射体（persistent scatterer，PS）为代表的单主影像时间序列InSAR技术^[21]和以短基线集（small baseline subset，SBAS）为代表的多主影像时间序列InSAR技术^[22]. 为了充分利用SAR影像数据，提高观测数据的时间采样率，避免出现干涉对时空基线较长而产生去相干的影响，采用SBAS-InSAR技术观测水布垭面板堆石坝视线向的变形. 该方法通过设定时空基线阈值将获取的数据进行组合配对得到一系列短基线差分干涉图，能够较好地克服空间去相关现象. 在求解形变速率时，采用奇异值分解法（singular value decomposition, SVD），将研究区域不同时间的形变连接成时间序列的形变场. 将ALOS PALSAR FBD模式下获取的影像数据进行重采样，获得与FBS模式下相同分辨率的数据，并将数据裁剪至合适的范围，如图5所示，可以清楚地看到清江河道、水布垭面板堆石坝坝体轮廓、溢洪道等.

通过设置时间、空间基线阈值，生成干涉对. 时间基线和空间基线分布如图6所示. 图中，B为空间基线. 利用轨道信息与外部 DEM 数据，对干涉对进行差分处理，去除平地效应与地形效应，得到差分干涉图. 为了进一步抑制噪声的影响，采用Goldstein滤波方法消除差分干涉图中的噪声，并利用三维最小费用流方法进行相位解缠^[23]，通过Delaunay三角网的其他高相干离散点来帮助解缠低相干性区域. 最后，执行标准的SBAS反演步骤获得水布垭面板堆石坝视线向的变形速率.

为了获得时间序列的方位向形变，结合多孔径InSAR和SBAS-InSAR技术的方法，根据与SBAS方法相同的规则改进多孔径InSAR技术的数据处理流程，如图7所示. 将SAR影像数据集划分为以较小的时间和空间基线为特征的干涉组合，利用多孔径InSAR方法生成MAI干涉图，并应用SBAS反演步骤得到方位向时间序列的变形信息.

将InSAR观测结果进行地理编码，使其位于WGS84坐标系. 为了更直观地展示坝体变形结果，进行矢量裁剪处理. 如图8（a）所示为由时间序列InSAR技术观测的水布垭面板堆石坝2007年2月—2009年1月年平均沉降速率v_u（由视线向分解得到）；如图8（b）所示为由多孔径InSAR技术观测的同时间段水布垭面板堆石坝卫星方位向的年平均变形速率v_AZI. 虽然多孔径InSAR技术干涉图的质量较全孔径干涉图的低且更易受时空失相关的影响，可能由于噪声而出现一些形变异常的区域，但是其依然直观地呈现出了堆石坝下游面的整体变形趋势，符合大坝变形的一般规律. 堆石坝下游面的外观变形基本成对称分布，左右岸发生指向坝体最大断面的变形，且左岸变形大于右岸变形，这是由于水布垭面板堆石坝的右岸岸坡比左岸岸坡平缓^[24]. 该时间段正值水布垭面板堆石坝的初期集中蓄水期，坝体发生明显的沉降变形，其下游坝面中上部高程的沉降较大，最大年平均沉降速率约为163 mm/a，向下游和沿坝轴向向两岸逐渐递减.

水准测量不受大气噪声的影响，其单点监测精度高于时间序列InSAR技术，但受限于水准点的个数，不能反映坝体全域的变形情况. 堆石坝的变形具有较高的空间自相关性，在一定范围内，具有趋势一致性和量级渐变性的特点. 因此，首先在时间序列InSAR观测结果提供的变形范围、时空分辨率、发展趋势等先验信息的基础上，利用ArcGIS软件统计模块中的克里金插值方法将高精度地面水准测点变形监测结果进行空间插值，同化地面水准测量和InSAR观测的空间分辨率.在此基础上，将水准测量数据与分解到垂直向的时间序列InSAR观测结果进行数据融合. 如图9（a）所示为测点WS07~WS10的数据融合结果. 图中，数据融合Ⅰ、Ⅱ分别为以水准测量数据与InSAR观测数据作为观测场的数据融合结果. 可以看出，在基于集合卡尔曼滤波进行数据融合前、后，位移的变化趋势总体基本一致，但数值呈现一定的差异，同化数值处于InSAR观测与地面水准测量值之间，且数据融合后的结果更接近观测场. 如图9（b）所示为以水准测量数据为基准，2种数据融合方式的误差对比图. 可以看出，当以水准测量数据为观测场时，均方根误差较小. 因此，以下采用以水准测量数据为观测场的同化方式进行融合.

通过融合水准测量数据弥补了InSAR技术由于各种噪声导致的观测误差，提高了变形监测数据的精度和可靠性. 由于各种噪声的影响，InSAR技术得到的部分像元间变形不连续（见图8（a）），在基于集合卡尔曼滤波进行数据融合后，变形分布更加光滑连续（见图10（a））. 进一步评估数据融合的精度，分析数据融合前、后InSAR观测结果与地面水准监测值的相关性和误差（见图10（b））. 图中，v_EnkF为数据融合后的年平均沉降速率，v_In为InSAR观测结果，v_ws为水准测量结果. 可以看出，与InSAR观测结果相比，数据融合后的可决系数R²由0.84提高到0.96，均方根误差RMSE由18.63 mm/a 降为9.66 mm/a，表明数据融合后的InSAR观测精度有较大提高.

通过时间序列InSAR技术与多孔径InSAR技术获取水布垭面板堆石坝高空间分辨率的面状变形. 考虑变形范围、时空分辨率、发展趋势等先验信息，基于集合卡尔曼滤波将InSAR观测数据与水准测量数据进行融合，通过式（9）~（13）将不同方向的变形监测数据进行联合解算，建立各方向变形监测数据与地表三维变形的关系，重构坝体的外观变形场，如图11所示. 重构的外观变形场直观反映了坝体在2007年2月—2009年1月期间的整体变形状态。图中，颜色和箭头分别表示垂直年平均变形速率v_EnkF和水平年平均变化速率v_H. 此期间正值水布垭面板堆石坝初期集中蓄水阶段，坝体出现以沉降为主的明显变形，年平均沉降速率普遍大于水平变形速率，其中沉降最大值出现在下游坝面的顶部和中部，并呈现沿河向向下游和沿坝轴向向两岸逐渐减小的规律；水平位移沿顺河向与沉降有相似的分布规律，最大值同样出现在下游坝面的顶部和中部，沿顺河向和向两岸逐渐减小，左右岸发生指向对岸的相互变形，且左岸大于右岸（由于右岸岸坡较缓）.

（1）InSAR技术提供了有效的堆石坝外观变形观测方法，高空间分辨率的观测结果可以有效弥补堆石坝常规监测技术以点代面、以局部代整体方式的不足. 利用时间序列InSAR与多孔径InSAR技术获取的视线向和方位向变形可以展现堆石坝的整体变形性态.

（2）基于集合卡尔曼滤波的多源数据融合技术可以实现“大范围、低精度、高效率”的新型监测技术与“离散点、高精度、低效率”的常规监测技术之间的优势互补. 以地面水准测量数据作为参照，数据融合显著降低了InSAR技术由于各种噪声造成的观测误差，提高了InSAR技术的变形观测精度.

（3）通过多源监测数据重构堆石坝的外观变形场能够全面地掌握堆石坝的变形性态，解决单一监测方式可能导致对实际变形漏判或误判的问题. 重构的外观变形场展示了水布垭面板堆石坝在2007年2月—2009年1月的变形特征，由于正值初期集中蓄水阶段，坝体变形明显并符合堆石坝的一般性规律.