基于卫星平台的时间序列InSAR分析方法具有高分辨率、大范围、厘米级地表形变场获取的能力^[1]，被认为是当前最具前景的形变监测技术^{[2, 3]}。该技术是针对传统InSAR技术易受时间和几何去相关以及大气扰动影响而提出的^[4]，推动了InSAR技术在地表形变监测中的广泛应用^[3]。

SAR传感器固有的侧视成像方式使得干涉相位只对视线方向(line of sight，LOS)上的形变敏感^[5]。因此，对于重复轨道获取的单独一组SAR数据集，当研究目标出现在几何畸变(如阴影)的发生区域时，无法提取有效信息来完成形变监测；或者目标的形变方向与卫星飞行方向平行，则在该数据集的视线向上无法得到很好的反映。这些都造成了研究区域形变信息的缺失，限制了时间序列InSAR分析技术对真实地表形变信息的精确反演。

随着越来越多SAR卫星的运行，联合具有不同成像几何参数的多卫星平台获取的SAR数据是当前该领域研究的趋势^{[6, 7]}，能为区域目标信息的获取提供相互补充。结合对目标的多次观测信息可对真实的地表形变做出准确估计，并为以上问题提供一种有效的解决方案。

1 多平台InSAR数据集联合估计原理

本文的多平台InSAR数据集联合估计方法结合了时间序列InSAR分析技术和带权最小二乘估计来提取无偏估计的地表形变量。该方法主要包含两部分：第一，多个SAR数据集分别进行时间序列分析，通过提取在时间序列上具有稳定散射特性的高相干点进行相位分析。依据不同相位源的时空特性来获得在观测时间段内，经过地形以及轨道、大气误差改正^[4]的形变速率场；第二，根据多平台数据之间不同的观测几何，建立有效的形变约束模型，然后根据模型利用最小二乘算法对多平台观测结果优化，最终得到具有最优统计特性的形变估计值。

本文将首先阐述联合估计的算法原理，分析具体的处理中由于引入了具有不同的观测几何、时间以及空间分辨率的SAR数据需要考虑的关键问题，并出处理方案。

1.1   联合估计算法

在局部参考系统内，LOS方向监测形变量R可分解为三个方向(南北N、东西E、垂直U)的分量，在同一个参考系内设定单位矢量 $\boldsymbol{\hat u} = \left[ {{u_{\rm{N}}}\; {u_{\rm{E}}}\; {u_{\rm{U}}}} \right]$ ，它的方向从地面指向卫星。

本文的目标是从观测量R中提取地表三个方向上的真实形变量v=[v_N v_E v_U]^T。根据式(1) 至少需要获取3个不在同一个平面的观测量。以单个点目标为例，假设从3个不同的视角以升轨或降轨来观测地面上的同一个点x，定义观测量R=[R_x¹ R_x² R_x³]，每一个R_xⁱ(i=1, 2, 3) 是从不同方向上获取的LOS观测量，则：

式中，R是LOS方向的观测量，对所有像素是已知的；A是一个与SAR系统观测几何相关的3×3矩阵; v是目标点与参考点的速率差。

在一个小区域(5×5 km²)内，针对城市区域形变，以沉降为主引入一个假设：在研究区域内不同目标点的视线向观测量的差异主要来源于垂直方向的不同。这样的假设是基于对实验区历史沉降情况的了解^{[1, 8]}以及为了简化计算。因此在小区域内，不同数据集测量获得的水平方向形变量应是一致的，那么当获取了三个数据集对M个点的观测量时，有如下简化模型：

式中，P是观测量的权阵；σ是m个观测点在3个数据集上观测量的标准差；f是利用标准差建立观测的随机模型。对观测量建立约束方程，获得地表真实形变量的最优无偏估值：

1.2   关键问题分析

(1) 观测几何

根据式(6) 可知，真实形变量估值的获取取决于系数阵是否可逆。系数阵是SAR传感器观测的几何参数的函数，即视线向每个观测量R_i和相应的真实地表形变量[v_N, v_E, v_U]之间存在关系：

式中，θ、φ是入射角和卫星飞行方向。系数阵满秩的条件是不同SAR数据获取时入射角和飞行方向的差异性。

目前的SAR卫星均以近极地轨道收集数据，即其飞行方向与北方向夹角约为10°(降轨角度约为170°，升轨约为350°)，获取数据的入射角一般在20°~45°之间。通过对式(7) 求偏导，可以近似估计InSAR技术对不同方向形变监测的敏感度。从式(8) 结果来说，观测量对垂直方向的形变最敏感，东西方向次之，南北方向最差。同一个方向上的形变，不同平台数据监测的敏感度存在差异，敏感度的差异性正是真实形变得以恢复的契机。在联合估计中，应充分考虑各个平台数据监测精度，本文通过修改式(5) 中的权阵来实现。

权阵描述平差问题中的观测量可靠性以及相互间统计相关性质可由观测量的方差-协方差确定。假设不同数据估计的LOS形变量均为独立观测量，以估计的形变量的标准偏差来确定观测量的权重，这是一种简单有效的选取方法。此外，还有许多采用验后估计来确定平差中采用的观测量权阵^[9]，如胡俊^[7]采用最小二乘准则和Helmert方差分量估计融合升降轨道获取的方位向和距离向形变，来反演Bam地震地表三维形变场。

(2) 时间周期

当前主要的星载SAR数据，ENVISAT ASAR为35 d，ALOS PALSAR为46 d，TerraSAR-X为11 d。考虑到不同卫星重访的时间周期不同，难以实现多个卫星平台在同一时刻获取特定研究区域数据。时间采样率的差异会造成某个观测时间点上数据出现缺失，法方程的系数阵出现秩亏，从而无法求解真实形变量。因此，本文采用在时间段内的线性形变(形变速率)进行分析，在一个较长共同覆盖时间段内，探测目标的行为应该具有一致的总体趋势。同时，目标在时间段内的演变行为可作为目标机制解译与灾害防治的重要信息。

(3) 空间分辨率

卫星飞行轨道的不同和观测时的传感器参数差异带来的是物理特性和相干点密度的区别，导致不同数据集间的配准困难。因此，针对不同数据的配准，需要保证配准的数据集有相同的空间尺度。

针对升降轨数据的融合，已有相应的实现方法^{[6, 7]}。值得注意的是，由于SAR侧视成像的特性，在升降轨道数据集的融合时，应剔除非公共的相干点(例如位于高耸建筑物不同墙面的相干点)，可通过在配准同时综合考虑相干点幅度值和高程值来完成。

本文的处理策略是，首先对不同数据集的高相干点进行局部最大窗口滤波，保留各个数据集窗口内最强散射点^[10]；其次，高相干点数据集的融合不是简单的合并，不同平台的观测数据需要经过投影计算到相同的地图坐标系下，同时为保证数据集间的空间一致性，用克里金内插方法从不同数据集获取规则格网数据，滤除高程值差异过大的点对；最后，通过地理编码数据进行最终的配准^[6]。

2 试验及分析

2.1   试验区及数据集

以上海为试验区，上海积累了丰富的SAR数据(例如18景TerraSAR-X、16景ENVISAT ASAR、20景ALOS PALSAR)，我们前期对上海已展开了研究^{[1, 8]}，对其沉降情况较了解，加之该地区完备的水准监测能验证融合算法的精确性。数据覆盖范围如图 1, 其参数如表 1。

图  1  TerraSAR-X、ALOS PALSAR和ENVISAT ASAR的数据范围(三角点：水准点数据监测点位)

Figure  1.  Tracks of TerraSAR-X, ALOS PALSAR and ENVISAT ASAR Data (Triangular Point: Location of Levelling Point)

下载: 全尺寸图片 幻灯片

表  1  实验用SAR数据的基本参数

Table  1.  Basic Parameters of Experimental SAR Data

卫星	时间	方向	θ/(°)	φ/(°)
TerraSAR-X	2009-03-28~2010-01-30	降轨	25.7	190.72
ALOS PALSAR	2007-01-07~2010-07-18	升轨	36.8	347.21
ENVISAT ASAR	2009-01-26~2009-09-13	升轨	22.1	346.80

下载: 导出CSV 

| 显示表格

2.2   多平台时间序列反演结果分析

3个试验数据集在共同观测时间段(2009年~2010年)的InSAR反演结果如图 2所示。整体结果在沉降区域分布上有很好的一致性，即在共同覆盖区域内，3种数据探测到的沉降区域是一致的。

图  2  典型实验区的沉降速率

Figure  2.  Subsidence Velocity of Typical Experiment Area

下载: 全尺寸图片 幻灯片

为了充分验证联合算法的特性，选取数据共同覆盖范围内具有较明显特征(较大沉降速率)的区域作为联合估计的典型实验区，即图 2中红色矩形框区域，其具体细节见图 3。

图  3  典型实验区的沉降速率场

Figure  3.  Velocity Fields of Typical Experiment Area

下载: 全尺寸图片 幻灯片

图 3显示了3种数据结果有所差异。3种数据的分辨率、波长、轨道方向及入射角度的差异导致沉降量值具有不同的敏感度，同时也是导致图 3视线向监测结果差异的原因。

3种数据结果各有特点。ALOS PALSAR反演结果噪声点较多，掩盖了区域形变分布特征，影响目视解译，然而ALOS的穿透性使得其即使在植被覆盖区也能够获得可信相干点，该特点值得应用。ENVISAT ASAR和TerraSAR-X数据反演结果形变分布更为接近，但TerraSAR-X高分辨率、短重访周期等相干性保持更好，能探测到更多的点。

2.3   融合试验结果及评价

采用本文方法在典型实验区提取的多平台联合估计结果如图 4所示，整体的形变分布情况清晰，有效地抑制了噪声。

图  4  典型实验区的多平台融合垂直方向估计结果

Figure  4.  Subsidence Velocity of Typical Experiment Area by Multi-platform Fusion

下载: 全尺寸图片 幻灯片

除了对形变情况的目视解译效果的改善，本文更加关注该方法对监测可靠性的提高。为此，本文对比了一个目标点在观测的共同时间段(2009~2010年间)，3种数据集InSAR反演累积沉降量与融合结果，如图 5(a)所示。对于ALOS PALSAR数据而言，时间监测的密度有很大改善，并且能较好抑制L波段数据相位噪声带来的误差对融合结果的影响；对于ASAR数据，除了改善时间监测密度，也削弱了相关误差源对沉降量估计影响，使得融合结果更加平滑；对于TerraSAR-X数据，改善了X波段数据监测形变数量级偏小的问题。

图  5  多平台融合结果

Figure  5.  Velocity Comparison Between Fusion Results and Three Respective InSAR/levelling

下载: 全尺寸图片 幻灯片

同时，将本文方法获得的无偏估计速率值与获取的18个同期观测的水准控制点的监测数据进行对比, 具体数据如表 2所示。本文并没有使用外部数据对结果校正，仅仅对三种数据结果进行联合估计，结果与地面测量能较好吻合。可见，融合算法能显著提高监测的精度，其示意图如图 5(b)。

表  2  水准观测数据与InSAR反演结果以及多平台融合结果的对比(估计速率/(mm·a^-1)

Table  2.  Comparison Among Levelling Data and InSAR Derived Results and Multi-platform Fusion Results/(mm·a^-1)

	多平台	TerraSAR	ASAR	PALSAR
平均误差	2.70	4.87	5.64	13.10
标准差	3.19	5.49	6.45	18.44

下载: 导出CSV 

| 显示表格

图 5(b)中存在个别水准点与InSAR融合结果差异较大，这是因为相干点位与水准点位并非准确对应。InSAR监测的空间密度较高，一个水准点位附近会有多个高相干点的监测结果。对比验证时，在水准点位附近依据高相干点的地理编码结果进行缓冲区分析，提取水准点位一定距离的高相干点反演结果与水准监测结果对比。也就是说，在缓冲内的高相干点存在与水准点形变行为不同的高相干点，则会引起对比结果偏差。在大多情况下，一定空间距离分布内的高相干点是具有相似形变行为，但也存在局部的空间不均匀形变，如建筑物高楼形成的高相干点的监测结果会与周围地面高相干点有差异。

3 结语

本文结合时间序列分析方法和最小二乘平差原理，给出了一种适用于多平台SAR数据融合的算法，讨论了由于不同平台数据的引入带来的成像几何参数、空间和时间分辨率的差异融合问题，给出相应的解决方案。

本文引入多平台数据，实现了精确的沉降速率反演，提供了更可靠的监测数据。研究区域的联合结果较好的平衡了各种数据之间的差异，使得形变分布更加容易解译；与同期地面观测水准数据对比结果证明了该方法的有效性，其平均误差不超过3 mm，提高了反演的精度。

下一步的研究将在具有更复杂形变区域开展真实形变的三维反演，以拓展该方法的可用性。