邢学敏 贺跃光 闻德保 周访滨

1 长沙理工大学特殊环境道路工程湖南省重点实验室，长沙市万家丽南路960号，410014

2 长沙理工大学交通运输工程学院，长沙市万家丽南路960号，410014

3 中南大学地球科学与信息物理学院，长沙市麓山南路932号，410083

永久散射体差分干涉测量技术（PSInSAR）广泛应用于区域地表形变的探测中［1－3］。然而，受区域相干性的限制，可能导致一些低相干区域内探测不到足够数目的PS点。人工角反射器技术（CR）是通过在研究区域内安装人工角反射器，通过对这些点上的相位进行建模分析，以求解其地表变形速率。人工角反射器安装灵活，当研究区域内探测不到PS点或当PS点密度不足时可发挥较强的优势［4］。此外，在CR 点上的GPS接收机还可获取CR 点的高程及坐标信息，一方面可大大避免在差分干涉处理过程中由于地理编码而引入的误差，另一方面可作为PS 网络解缠的起算数据，避免由于人为选取假设稳定点而引入的不确定性。本文将CRInSAR 技术与PSInSAR融合起来用于解算地表线性变形速率。设计了详细的算法流程，分别利用模拟实验及真实数据实验对算法进行验证。

1 PSInSAR 算法

假设有N＋1幅时间序列影像，利用二轨差分干涉处理方法可对应生成N幅差分干涉图。选取一定数量的PS点后，将任意相邻的两个PS点进行连接，可组建成PS 基线网络。针对网络中任意一条PS基线（假设包含第i和第j个PS点），其相位差可以表示为［5－6］：

其中m为干涉对序号；分别为两PS点的线性变形速率增量和高程改正增量；分别为相邻两PS点在第m幅干涉图中对应的相位增量和整周模糊度增量；为高程相位改正系数，其中分别表示时空基线，θ为影像的入射角，Ri表示雷达传感器与PS点之间的距离；表示残余相位，主要由大气延迟相位、非线性形变相位和噪声组成［7］。

针对上式未知参数Δvi，j、ΔδHi，j及的求解可利用LAMBDA算法进行［8］，ΔδHi，j及Δvi，j的估值结果可作为空间维相位解缠的观测数据，利用间接平差的方法求解出各PS 点上的绝对线性变形速率v及高程改正值δH。

2 CRInSAR 算法

利用CRInSAR 技术进行形变监测时，需要在影像上先提取出CR 点的行列号信息，然后建立如下函数关系模型［9］：

式中p表示CR 点序号，m表示干涉对序号，标有p和m的各分量均表示第p个CR 点与参考CR点的对应参数差；为干涉相位差；为整周模糊度差；为地形引起的相位差；vp为相对于参考点的形变速率差；Tm为时间基线；是残余相位，与（1）式中ε含义相同；地形相位，其中λ是雷达波长，Bm为垂直基线长度，Rp为CR 点与卫星位置间距离，ΔHp是CR 点与参考点的相对高程。

实际计算中，假设有N＋1 幅SAR 影像，有M＋1个CR 点，选取一个CR 点为参考点，一幅SAR 影像为主影像，可对应生成N幅干涉图。当CR 点间距在1 000 m 范围内时，大气相位影响可忽略，因此，可不考虑残余相位影响［10］。而均为已知量，则在式（2）中未知参数仅为整周模糊度增量Δk、形变速率v。利用LAMBDA 算法，可逐步分离出未知参数，进而得出最终的CR 点上的形变速率值vp。

3 PSInSAR 与CRInSAR 融合解算算法

利用LAMBDA 算法获取了相邻PS点的线性变形速率增量及高程改正增量后，即可作为输入数据，用于实现任意PS 点上绝对线性变形速率及高程改正值的求解，这一过程又被称为空间维相位解缠。传统空间维解缠会受主观选取稳定点的制约，带有很大的不确定性［11－12］。我们预先在研究区域内安装CR 点，获取其线性变形速率和高程改正值作为PS基线网络的空间维解缠的约束数据，利用间接平差方法求解各PS 点的绝对线性变形速率值和高程改正值，进而实现PS基线网络的空间维解缠。图1为融合算法中CR点与PS点相对位置示意图，以线性变形速率为例，针对任一条基线边，我们都可以将Δvi，j作为间接平差函数模型的观测值，线性变形速率参数vi、vj则为待求未知参数。间接平差函数模型可表示为：

图1 PS与CR 点相对位置示意图Fig.1 Diagrammatic sketch of reference position between CRs and PSs

上式中，Hcr为CR 点上的线性变形速率，可通过§2介绍的算法预先计算得出；δHcr为CR 点上的高程改正结果，可以通过式δHcr＝Hcr－Hdem进行计算。其中Hcr表示CR 点上的GPS接收机获取的高程值，Hdem表示CR 点对应的外部DEM 高程结果。由此，可利用最小二乘原理将研究区域内所有PS点上的线性变形速率及高程改正参数求解出来。根据上述讨论，将PSIn－SAR与CRInSAR 融合解算算法流程图表示在图2中。

图2 融合算法流程Fig.2 Flowchart of combined calculation algorithm

4 模拟实验

4.1 实验流程

为验证前述算法流程，设计一套模拟实验。首先利用peaks函数来模拟线性变形速度场，如图3 所示。模拟速度场的总面积为100 像素×100像素，区域内以无形变为主，变形区域有下沉、抬升，呈对称分布。从中随机选取200个像素作为参考PS 点，12 个像素作为CR 点。为增大PS点密度，再随机选取600个点位作为扩展PS点。对参考PS点和扩展PS点共同建立起分级PS基线网，类似于大地测量中分级控制网的方式（图4）。参考PS点类似于大地测量控制网中一级控制点，而扩展PS点为二级控制点，需要将参考PS点与其邻近的扩展PS点相连接构成基线边。高程改正值采用高斯模型进行模拟，其均值为0m，标准偏差为±3m。相位整周模糊度利用随机整数模拟器模拟。这样，PS点与CR 点的模拟真实速率值和模拟真实高程改正值均已给出。

图3 模拟线性变形速率场Fig.3 The simulated linear velocity field

图4 参考PS点与扩展PS点位置示意图Fig.4 Distribution of reference PSs and extended PSs

由式（1），差分干涉相位可以根据上面介绍的各模拟分量组建而成，其中涉及到的影像参数可直接在SAR 影像的头文件中读取（如时间基线、空间基线、高程相位转换系数）。由于现有数据主要为ALOS 卫星影像，因此本实验中选取的是PALSAR 影像各类参数。设置随机噪声均方差为0.5rad。

完成各项模拟工作后，我们利用融合算法求解各PS 点上绝对下沉速率及高程改正参数值。首先，对所有PS及CR 点进行基线网络布设。采用Delaunay算法将这200 个PS 点与12 个CR点共同组成的212个一级控制点建立成参考基线网络，共生成608条基线边。再将600个扩展PS点与212个一级控制点构建成扩展PS 网，在布网过程中限制二级网络基线边长不超过5像素，则共生成2 295条基线边。如图5所示，图中三角形表示CR 点位置。然后，依据式（1）构建基线边上任意两个PS点的差分相位，并利用LAMBDA算法对式（1）进行时间维相位解缠，求解出所有基线边上的线性变形速率增量及高程改正值增量。最后，将CR 点上预先给出的线性变形速率及高程改正值模拟值作为基线网络空间维解缠的起算数据，根据间接平差方法，求解出所有PS点上线性变形速率和高程改正结果。

图5 模拟的参考基线网和扩展基线网Fig.5 Simulated network of referenced PSs and extended PSs including CRs

4.2 实验结果

模拟实验预先给出了所有PS点的线性变形速率和高程改正模拟值，将其作为真值，与利用融合解算方法求解出的线性变形速率和高程改正计算值进行比较，验证算法的精度。对参考PS 点和扩展PS点上的线性变形速率计算值与模拟真值的差异作为其对应精度，将精度分布直方图表示在图6中。从图6明显看出，无论参考PS点还是扩展PS 点，线性变形速率计算值与模拟真值差异总体分布在－4～＋4mm 之间，说明其与模拟真实结果十分接近。根据统计，无论参考网还是扩展网，都有接近60%的PS点线性变形速率差异分布在－1～＋1mm，仅有个别PS点的差异超过10mm。计算得出，参考网点的线性速率差异均方根误差为±0.37mm／a。

图6 参考网PS点及扩展网PS点线性变形速率计算值与模拟真值差异分布直方图Fig.6 Statistics bars of Linear velocity rates errors at the reference PSs and extended PSs

图7为真实PS点线性变形速率场与计算获取的PS 点线性变形速率场的对比效果图，从中可以看出二者颜色分布非常接近，但仍有部分区域存在细微的颜色差异。如图中红色方框标出的A 区内，计算结果中有约一半的点出现隆起，颜色表现为桔黄色，对应的变形速率为0.02～0.04 cm／a，而模拟真值（图7（b））中这一区域内只有少量点出现隆起，大部分点表现为无形变。在B 区内，计算结果（图7（a））中出现了零星分布的下沉点，而模拟真值结果（图7（b））中整个区域内所有点都表现为隆起。出现这种情况主要是由于在模拟过程中加入了噪声，另外，在算法求解过程中没有考虑非线性形变及大气延迟等因素的影响。在实验过程中获取的高程改正值结果均方误差为±0.5m，说明利用这一算法进行求解，可以通过增加高程改正的方式修正低精度的外部DEM，且修正结果精度为亚m级，可以减少在外部DEM 订购方面的成本，甚至可以不使用外部DEM 去除地形信息，直接获取待求PS 点上的高程信息及线性变形速率。

图7 计算获取的PS点线性变形速率与模拟真值对比Fig.7 Comparison of linear velocities at all PSs with the simulated real value

5 结 语

本文将CRInSAR 与PSInSAR 两种时间序列InSAR 地表变形监测技术融合，用于解算地表线性变形速率，设计了详细的算法流程。算法总体思想是利用CR 点上获取的线性变形速率结果和高程改正值作为PS基线网络的空间维解缠间接平差函数模型的起算数据，利用最小二乘原理求解出所有PS点上的线性变形速率与高程改正参数估值。这一算法可在一定程度上避免人为选取假设稳定点的不确定性，融合了CR 与PS两种技术的优势，可应用于PS 点较为稀少且缺少先验变形信息的研究区域。为验证算法的可行性，本文设计并实现了模拟实验流程。从模拟实验的结果可以看出，这一融合算法解算线性变形速率可达到mm级精度，具备一定的可行性及可靠性。下一步将重点研究不同CR 约束网型对融合算法精度的影响，指导CR 点的安装与布设，以起到对PS基线网络的最佳约束效果。

［1］Ferretti A，Savio G，Barzaghi R，et al.Submillimeter Accuracy of InSAR Time Series：Experimental Validation［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（5）：1 142－1 153

［2］Crosetto M，Monserrat O，Iglesias R，et al.Persistent Scatterer Interferometry：Potential，Limits and Initial Cand X－Band Comparison［J］.Photogrammetric Engineering and Remote Sensing，2010，76（9）：1 061－1 069

［3］Osmanoĝlu B，Dixon T H，Wdowinski S，et al.Mexico City Subsidence Observed with Persistent Scatterer InSAR［J］.International Journal of Applied Earth Observation and Geoinformation，2011，13（1）：1－12

［4］Xia Y，Kaufmann H，Guo X.Landslide Monitoring in the Three Gorges Area Using D－InSAR and Corner Reflectors［J］.Photogrammetric Engineering and Remote Sensing，2004，70（4）：1 167－1 172

［5］Ferretti A，Prati C，Rocca F.Permanent Scatterers in SAR Interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing，2001，39（1）：8－20

［6］Ferretti A，Prati C，Rocca F.Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing，2000，38（5）：2 202－2 212

［7］Kampes B M，Hanssen R F.Ambiguity Resolution for Permanent Scatterer Interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing， 2004， 42（11）：2 446－2 453

［8］Xia Y，Kaufmann H，Guo X.Differential SAR Interferometry Using Corner Reflectors［C］.Geoscience and Remote Sensing Symposium，2002

［9］Liu G，Luo X，Chen Q，et al.Detecting Land Subsidence in Shanghai by PS－networking SAR Interferometry［J］.Sensors，2008，8（8）：4 725－4 741

［10］陈强，丁晓利，刘国祥，等.雷达干涉网络的基线识别与解算方法［J］.地球物理学报，2009，52（9）：2 230－2 236（Chen Qiang，Ding Xiaoli，Liu Guoxiang，et al.Baseline Recognition and Parameter Estimation of Persistent－scatterer Network in Radar Interferometry［J］.Chinese Journal of Geophysics，2009，52（9）：2 230－2 236）

［11］陈强，丁小利，刘国祥.永久散射体雷达差分干涉应用于区域地表沉降探测［J］.地球物理学报，2007，50（3）：737－743（Chen Qiang，Ding Xiaoli，Liu Guoxiang.Radar Differential Interferometry Based on Permanent Scatterers and Its Application to Detecting Regional Ground Subsidence［J］.Chinese Journal of Geophysics，2007，50（3）：737－743）

［12］Xia Y.CR－Based SAR－Interferometry for Landslide Monitoring［C］.Geoscience and Remote Sensing Symposium，2008