一种有效的机载双频干涉地形高程重建方法

刘华有 郑明洁 张 衡* 王 宇 秦小芳

①(中国科学院电子学研究所 北京 100190)

②(中国科学院大学 北京 100049)

1 引言

合成孔径雷达干涉(Interferometry Synthetic Aperture Radar, InSAR)技术，是广泛应用在地理测量中的一种雷达遥感技术。利用同一场景的两幅合成孔径雷达(Synthetic Aperture Radar, SAR)复数据，提取相位信息来获取地面的高程信息[1,2]。可以应用于地形测绘、形变监测、海洋应用、以及地球动力学应用[3–5]。

传统的单波段InSAR系统只需两幅SAR图像形成一幅干涉图，即可实现高程反演，其处理流程[1]主要包括复图像配准、干涉相位去平地、干涉相位滤波、相位解缠、数字高程模型(Digital Elevation Model, DEM)重建等步骤。虽然单波段InSAR系统数据获取相对容易，处理流程相对简单[3]，但是单波段相位解缠必须满足Iton假设[6]，无法实现地形起伏比较大的高程反演。90年代以后多通道InSAR技术快速发展，多通道InSAR技术主要包括多频率干涉和多基线干涉两种。1994年，Xu Wei等人[7]提出了中国剩余定理(Chinese Reminder Theorem,CRT)、投影法和线性组合法3种多通道联合相位解缠算法。2000年，Pascazio 和Schirinzi[8]提出将大带宽的SAR数据进行分割来获取多频率干涉数据，并且使用最大似然估计(Maximum Likelihood,ML)来实现多频率联合相位解缠。2003年，Hoge等人[9]提出了一种基于频域奇异值分解(Singular Value Decomposition, SVD)的配准方法，可以很高效的实现图像亚像素级配准。2004年，Ferraiuolo等人[10]提出了一种高斯能量项的最大后验估计(Maximum A Posterior, MAP)的算法来多通道联合高程重建，该方法使用马尔科夫随机场来建模地形高程的先验分布，并通过仿真的多频率干涉验证了算法的有效性。2009年，Ferraioli等人[11]提出了一种基于图割法的变分能量项最大后验估计(Maximum Posteriori Estimation of Total Variational Energy Term, TV-MAP)的多通道联合高程重建方法。Ferraiuolo和 Meglio等人[12]深入分析最大似然估计和最大后验估计等两种多通道联合相位解缠算法，证明了最大后验估计法鲁棒性优于最大似然估计。2011年，于翰雯等人[13]提出了一种基于聚类分析的多通道联合相位解缠绕算法，该方法将具有相同模糊矢量的像素点聚成一类，随后逐类进行联合解缠。Deledalle等人[14]提出了一种非局部参数估计(Nonlocal Interferogram Estimation, NL-InSAR)的方法，可以实现对干涉相位、相干系数、幅度图的精确估计。2013年，袁志辉[15]提出了一种基于闭合形式鲁棒性中国剩余定理(CRT)的多通道联合高程重建算法。2016年，曾涛[16]等人提出了结合最大似然估计相位梯度和枝切法的多频率联合相位解缠算法。2017年，斯奇等人[17]提出了一种基于最大后验框架的聚类分析多基线干涉SAR高度重建算法。

双频InSAR技术利用不同频段的两幅干涉图进行联合高程反演，可以有效地提取地形复杂区域的高程信息。为实现高精度、强鲁棒性的高程反演，双频干涉数据处理需要解决不同频段之间的非相干配准、干涉相位降噪、联合相位解缠等一系列问题。针对这些问题，结合现有图像配准、相位滤波、多通道InSAR技术等研究进展[18,19]，本文提出了一种有效的双频干涉SAR地表高程重建方法。该方法对常规处理流程中的关键步骤进行了改进，首先采用NL-InSAR技术对幅度图、相干系数、干涉相位通过迭代进行精确估计，利用各个波段的幅度信息来实现不同波段的干涉相位的配准。然后采用聚类分析技术对联合解缠相位标记有效点和噪点，并利用有效点对联合解缠相位进行均值滤波。通过这些改进，提高了高程重建效果与高程反演精度。目前该方法已经应用在实际机载双频InSAR系统中，取得了较好的结果。

2 高程重建算法

本文针对常规的双频干涉处理流程中的关键步骤进行了改进，改进后的流程图如图1所示，主要包括频段内图像配准、去平地相位、干涉参数精确估计、频段间干涉相位配准、双频联合相位解缠、解缠相位噪点剔除等步骤。下面通过比较分析，将给出各个步骤的处理方法，并在机载实测双频干涉处理中验证其有效性。

2.1 图像配准

图像配准是通过调整主、辅SAR复图像，使得两幅图像上相同的点对应地面场景上的同一点。配准的步骤如下[20,21]：

步骤1 主辅图像选取若干控制点，采用相关法计算控制点处粗偏移量，精度为1个像素。

计算控制点处的偏移量可以采用空域相关函数法的计算基于窗口的配准方法，分别在主图像s1和辅图像s2以控制点为中心建立匹配窗和搜索窗。滑动搜索窗求得窗口的相关函数值，相关函数最大的滑动位置即为图像的粗偏移量(u，v)；

步骤2 利用相关法求得图像粗偏移量，并修正辅图像控制点的位置，为得到亚像素级的精配准，可以通过空域相干系数法[21]、傅里叶变换法[21]、频域SVD法[9]求得控制点处的精确偏移量；

步骤3 在求得控制点处精确偏移量后，可以利用二项式进行拟合：

其中， (x，y)为主图像的控制点坐标，通过最小二乘法可以得拟合系数，再利用拟合系数得到辅图像整幅图像配准后精确的位置，通过线性插值把辅图像插值到求得的精确位置来实现精配准[21]。

下面采用某机载双频干涉系统飞行的数据为例来分析图像配准技术。

空域相关法、傅里叶变换法、频域SVD分解法这3种方法都具有较高的配准精度，本文分别采用了这3种方法对上述图像进行了配准，配准结果如图2所示。由于频域SVD分解法运算效率最高，因此机载双频干涉高重建过程中采用频域SVD法配准。

2.2 去平地

配准后的干涉相位由于平地相位的存在，使干涉相位缠绕量较大，干涉条纹过密，会影响后续的相位解缠。因此在相位解缠之前，需要去除平地相位，以便降低联合相位解缠的难度。去平地方法主要有频谱法和轨道法两种，频谱法是将干涉复图像变换到频域，找出频谱中最大值对应的频率并搬移到零频。轨道法是根据轨道参数计算平地相位。两种方法的处理结果如图3所示。由于在实测数据处理中频谱法并不能有效地去除平地相位，因此，本文采用轨道法去平地。

2.3 干涉参数精确估计

为了提高不同频段配准和联合解缠的精度，采用NL-InSAR技术分别对各个频段的干涉相位、相干系数、幅度进行精确估计[14]，再进行不同波段配准，可以有效提高不同频段之间干涉相位的配准精度，并且实现了对干涉相位去噪的同时较好地保持了干涉相位的纹理信息。图4给出了NLInSAR滤波前后的干涉相位对比情况。从图中可以看出，NL-InSAR在去除干涉相位存在的斑点噪声。NL-InSAR滤波前的幅度图和滤波后的幅度图如图5所示，对比表明NL-InSAR在保持SAR图像边缘细节信息的前提下能较好地去除SAR图像的斑点噪声。NL-InSAR估计前的相干系数如图6所示，可以看出，NL-InSAR估计后的相干系数的可靠性更强，可以有效辅助后续双频联合解缠。

2.4 频段间干涉图配准

交叉频段干涉相位配准精度往往较低，本文采用的方法为首先进行同频段图像配准、滤波、去平地，通过利用各自主图像的幅度信息拟合出多项式系数，进而来实现频段间干涉相位的配准。

由于不同频段的主图像幅度存在较多的斑点噪声，造成不同频段的主图像幅度图的相似性较差。由图5可知，NL-InSAR可以有效地对SAR图像的幅度图进行去噪。图7(a)为不进行NL-InSAR幅度滤波的C波段主图像幅度图，图7(b)为不进行NLInSAR滤波的X波段主图像幅度。可以看出不进行NL-InSAR滤波的不同波段主图像幅度图较多的斑点噪声，造成不同波段的主图像幅度图之间相似性较差。图8(a)为采用NL-InSAR对幅度滤波后的C波段主图像幅度，图8(b)为采用NL-InSAR对幅度滤波后的X波段主图像。NL-InSAR可以很好地去除SAR图像幅度图的斑点噪声，同时较好地保持了图像的边缘细节信息，提高了不同波段主图像幅度图之间的相似性。因此本文采用经过NL-InSAR滤波后的两个频段的主图像幅度拟合出式(1)的多项式系数，从而实现对两个频段的干涉相位图之间的配准。针对实测数据不同频段配准前后的不同波段之间相关系数图如图9所示。图9(a)所示的配准前X波段和C波段的相关系数较低，均值为0.46，图9(b)所示的不采用NL-InSAR对幅度进行滤波的配准后的相关系数相对较高，均值为0.90。图9(c)所示为采用NL-InSAR对幅度进行滤波的配准好的相关系数最高，均值为0.94，验证了本文提出的不同波段干涉相位配准方法的有效性。

2.5 多通道联合相位解缠

目前，多频段联合相位解缠算法大致有基于统计信号处理和非统计信号处理两类，本文仿真对比了非统计信号处理类的鲁棒性中国剩余定理(CRT)[15]，统计信号处理类的最大似然估计(ML)[8]和变分最大后验估计(TV-MAP)[11]。在实测数据处理中，双频联合相位解缠算法对算法鲁棒性要求较高，本文选用鲁棒性较强的变分最大后验估计作为双频联合解缠的处理方案。最大后验估计相位解缠就利用观测相位Ψ并且加入真实相位的先验分布P(Φ)来估计真实相位Φ，即：

假设所有像素点之间是相互独立的，那么：

对于N0个波段相位联合解缠情况下，相应似然函数为：

其中单个波段的观测相位的似然函数为：

我们可以使用马尔科夫随机场来描述图像的局部特性，根据MarKov-Gibbs等价性可以用Gibbsian分布来描述马尔科夫先验分布[8]：

变分能量项最大后验估计使用的是一种在图像领域常用的整体变分(TV)模型，其能量函数定义为：

式(7)中，β是一个固定的常量，用来调节真实相位的光滑性与观测相位和真实相位之间相似性的平衡的常数，Np为像素点p的邻域的集合，wp，q表征邻域连通性，若Np为4邻域模型，其值为1。

把式(7)代入式(6)结合式(2)–式(5)并取对数且添加负号，可以得到变分能量项最大后验估计的在TV模型下的相位估计公式为：

上述优化模型可以通过图割理论求解[11]，通过求解上述优化模型可以求解得到各个波段干涉相位的缠绕量和联合解缠绕相位。

对添加相同噪声的仿真数据，采用上述3种方法分别进行联合解缠结果如图10所示，用来仿真干涉相位的原始相位相减的残差图如图11所示。通过上述仿真分析，进一步验证了变分能量项最大后验估计用于双频联合解缠的鲁棒性。

2.6 聚类均值滤波

模糊矢量[13,22]定义为：

其中，Kc为C波段干涉相位的缠绕量，Kx为X波段干涉相位的缠绕量。将图像各个点的模糊矢量求和得到：

如图12所示，将两种标记结果相与，得到用于加权均值滤波的标记结果，对联合解缠结果进行如下加权均值滤波：

3 实测数据分析

下面采用实际飞行数据对本文方法进行验证，同时对比了单波段和双波段的高程重建效果。实测数据是由中国科学院电子学研究所航天微波遥感系统部提供的某场景双波段(C波段中心频率5.4 GHz和X波段中心频率9.6 GHz)的机载双频干涉数据。

实测数据包括同一场景的C波段和X波段主、辅SAR图像复数据，处理数据的场景如图16所示，为陕西澄城地区的某山区，其中地形起伏比较大区域为图中右上角的高架，不同波段配准之后的C波段干涉图和X波段干涉图分别如图17所示。

单波段使用最小代价流(MCF)解缠，处理结果如图18所示。使用最大似然估计方法双频联合解缠结果如图19所示，由于波段数较少，最大似然估计鲁棒性较差，完全无法实现正确的相位解缠。按照本文方法，双频联合解缠结果如图20所示。对图中B区域红框内地形起伏比较大的高架区域解缠结果放大对比如图21所示，从图中可以看出在地形起伏较大的区域，单波段也无法实现正确解缠，而使用本文提出的双频联合解缠算法可以正确解缠。

根据轨道参数反演的DEM如图22所示，图中下半部分为上半部分中沿虚线的地形走势，通过对比可知，无论是X波段还是C波段，都无法正确反演高架桥高程，只有联合X波段和C波段，才能正确反演。这也更加明确地说明了，在地形起伏比较大的区域，只有双频联合才可以有效地反演高程信息。

由于本文所用的实测数据场景中没有准确的地本文采用的方法是将谷歌光学地图面控制点，本文采用的方法是将谷歌光学地图作为主图像和SAR图像的幅度图进行配准，得到配准的多项式系数，利用拟合的多项式系数实现对DEM进行配准。在光学图像上选取5个控制点对DEM进行精度分析。单波段反演高程的均方根误差(Root Mean Square Error, RMSE)为4.5826 m，且存在部分地区高程反演失败，而双频联合反演的RMSE为1.3476 m，可以看出本文提出双频联合高程重建方法精度优于单波段高程重建。

4 结束语

本文提出了一种有效的双频干涉处理方法，分析并给出了算法中各个步骤的处理方案。针对双频联合处理中的双频联合配准、双频解缠绕等关键技术难题，给出了有效的解决方案，提高了双频联合反演DEM的精度和鲁棒性。采用常规单波段方法和本文方法对机载实际飞行数据处理，处理结果表明，单波段使用MCF进行相位解缠，无法实现大起伏地形区域的高程重建；而采用本文提出的处理流程能够准确地反演大起伏地形区域的高程，验证了本文所提出的双频干涉处理方法的有效性。