地形高程自适应的星载InSAR图像配准方法

张金强， 索志勇， 李真芳， 保　铮

(西安电子科技大学 雷达信号处理国家重点实验室，陕西 西安 710071)

星载干涉合成孔径雷达(Interferometric Synthetic Aperture Radar， InSAR)技术利用不同观测几何下，获取的两幅合成孔径雷达(Synthetic Aperture Radar，SAR)复图像间的干涉相位反演来获取地表三维地形[1-3]．为了从SAR复图像对中获取干涉相位信息，需要对其进行配准处理，以确保用于干涉处理的图像像素对应于地面同一散射单元[4-5]．图像配准精度将影响InSAR系统的最终产品精度[6]．因此，研究稳健高精度InSAR图像配准方法具有重要意义．已有的InSAR图像配准方法主要可分为两大类：第1类方法首先利用基于图像数据[7-8]或图像特征[9-10]的方法估计控制点处配准偏移量，然后根据二维二元低阶多项式模型计算其他像素处的配准偏移量．该类方法没有考虑地形高程对配准偏移量的影响，当基线较长或观测场景内地形复杂度较高时，利用二维二元低阶多项式模型拟合配准偏移量的精度较低[11]；当利用控制点处配准偏移量估计多项式系数时，控制点数量和位置分布对系数估计精度的影响较大；第2类方法首先利用几何配准法[12]计算每一像素的初始配准偏移量，然后利用控制点处高精度配准偏移量根据二维二元一阶多项式模型校正初始配准偏移量[11]，从而获取全场景高精度配准偏移量．该类方法需要进行逐像素几何配准处理，运算效率较低．

针对上述问题，笔者提出了地形高程自适应的降维图像配准方法．首先，给出了主辅图像间的配准偏移量计算模型；然后，为了追踪配准偏移量随地形高程的变化，根据配准偏移量计算模型提出了包含地形高程项的二维三元一阶多项式模型．利用该模型不仅可以精确拟合主辅图像间的配准偏移量，而且能够避免逐像素几何配准处理；最后，通过降维处理将需要利用控制点处配准偏移量估计的多项式系数减少为两个常数项系数，提高了系数估计的稳健性．实测数据处理结果验证了文中方法的精确性和稳健性．

1　星载InSAR降维图像配准方法

1.1　配准偏移量计算模型

主辅图像间的配准偏移量计算流程包括两步:首先，通过正向定位获取主图像像素对应的目标三维位置;然后，通过反向定位获取上述目标在辅图像上的坐标．针对主图像上某一像素(am，rm)，将成像参数和轨道数据代入下列3式：

所示的距离-多普勒模型[11-13]求解目标三维位置，即正向定位，其中，R表示雷达斜距；t表示方位时间；P(t)= (Px，Py，Pz)T，表示卫星位置矢量；T= (Tx，Ty，Tz)T，表示目标三维位置；fdi表示成像多普勒频率；V(t)= (Vx，Vy，Vz)T，表示卫星速度矢量；λ表示雷达波长；h表示目标高程；P(t)和V(t)随时间变化；Re和Rp分别表示赤道和极地半径，且Rp= (1-f)(Re+h)，f为参考椭球扁率；P、V和T均定义在世界测地系统(World Geodetic System-1984，WGS-84)坐标系下．

利用目标三维位置T、辅图像成像参数和轨道数据，首先根据式(2)计算该目标对应的辅图像方位时间t，然后根据式(1)计算时刻t目标到雷达的斜距R，最后根据

计算该目标在辅图像上的二维坐标(as，rs)，其中，t0表示图像方位向第1个像素对应的方位时间，fp表示脉冲重复频率，R0表示图像距离向第1个像素对应的雷达斜距，fs表示脉冲采样频率，c表示光速．上述过程称为反向定位．

1.2　配准偏移量拟合模型

由配准偏移量计算模型可知，目标在辅图像上的方位和距离坐标随其在主图像上的方位、距离坐标和高程变化．传统配准方法假设目标高程对配准偏移量的影响可以忽略，用含有方位和距离坐标项的二维二元低阶多项式模型拟合配准偏移量随主图像方位和距离坐标的变化．当主辅图像间的基线较长或观测场景内地形复杂度较高时，地形高程对配准偏移量的影响不能忽略．为了追踪配准偏移量随地形高程的变化，可用

所示的包含地形高程项的二维三元一阶多项式模型拟合主辅图像间的坐标关系，其中，d0～d3和g0～g3表示拟合系数．

1.3　系数估计

式(6)和式(7)所示的二维三元一阶多项式模型存在8个未知系数．传统配准方法利用控制点处配准偏移量估计上述未知系数，控制点数量和位置分布对系数估计精度影响较大．为了降低控制点对系数估计的影响，提高系数估计的稳健性，文中给出了利用主辅图像的成像参数和轨道数据计算一次项系数的公式．一次项系数计算公式的推导包括三步:首先，根据正向定位模型，推导目标三维位置相对于主图像方位时间、雷达斜距和目标高程的一阶偏导数;然后，根据反向定位模型，推导目标对应辅图像方位时间和雷达斜距相对于目标三维位置的一阶偏导数;最后，结合上述两步结果得到一次项系数的计算公式．

根据正向定位模型，将式(1)～式(3)中目标三维位置分别对方位时间、雷达斜距和目标高程求一阶偏导数，得到

其中，Am表示主卫星加速度矢量．其他参数与1.1节的参数含义相同，下标m表示主图像对应的参数．

根据反向定位模型，将式(1)和式(2)中方位时间和雷达斜距分别对目标三维位置求一阶偏导数，得到

其中，As表示辅卫星加速度矢量．其他参数与1.1节的参数含义相同，下标s表示辅图像对应的参数．

根据式(4)～式(5)和式(8)～式(14)，得到目标在辅图像上的二维坐标相对于其在主图像上的二维坐标和高程的一阶偏导数为

其中，d3和g3的单位为 像素/m．上述偏导数随主图像二维坐标和目标高程变化．由于星载SAR的雷达斜距较大[14]，在一定观测场景范围内，上述偏导数的变化对配准偏移量的影响可以忽略．因此，在实测数据处理过程中，利用观测场景中心处目标对应的主辅图像的成像参数和轨道数据，根据式(15)～式(20)计算d1～d3和g1～g3．

图1 文中方法的处理流程

1.4　处理流程

如图1所示，地形高程自适应的降维图像配准方法处理流程包括如下步骤：

(1) 利用观测场景中心处目标对应的主辅图像的成像参数和轨道数据，根据式(15)～式(20)计算d1～d3和g1～g3．

(2) 在主图像上选取一定数量的控制点，采用基于图像数据或图像特征的方法估计控制点处配准偏移量．

(3) 根据式(6)和式(7)，补偿控制点处配准偏移量的一次分量，利用剩余配准偏移量估计d0和g0．

(4) 根据式(6)和式(7)，计算其他像素处配准偏移量．

文中方法需要输入主图像像素对应的目标高程，其可以利用先验数字高程模型(Digital Elevation Model，DEM)通过正向定位获取．

2　实测数据处理与分析

利用TerraSAR-X分别于2010年3月8日和2010年3月19日对某地区进行重复航过观测任务时，获取的两幅SAR复图像来验证文中方法，取2010年3月8日获取的SAR复图像为主图像．主要系统参数如表1所示，SAR幅度图像如图2(a)所示．在图2(a)方框所示的区域内选择12个特显点以定量评估配准精度，如图2(b)所示．图2(c)为观测场景对应的航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission，SRTM)DEM，已将其投影到主图像斜距平面，观测场景幅宽为 15 km× 15 km，观测场景内目标高程变化为 695 m．

表1　TerraSAR-X主要系统参数

图2　TerraSAR-X重复航过数据

利用图1所示的处理流程对主辅图像进行配准处理，在主图像内均匀选择100个控制点．控制点处配准偏移量采用实相关函数法[7]估计，图像子块的大小为 64像素× 64像素．配准后主辅图像间的干涉相位图如图3所示．利用文中方法配准后的主辅图像可以获得清晰的干涉相位图，验证了文中方法的有效性．

为了验证文中方法的精确性和稳健性，利用蒙特卡罗实验对比分析文中方法与二阶多项式方法(利用二维二元二阶多项式模型拟合配准偏移量，并且多项式系数均利用控制点处配准偏移量估计)在不同控制点数量下的配准性能．蒙特卡罗实验次数为100次．在主图像内均匀选择900个控制点，控制点处配准偏移量采用实相关函数法估计，图像子块的大小为 64像素× 64像素．每次实验从上述900个控制点中随机选择一定数量的控制点估计多项式拟合系数．配准误差利用图2(b)所示的12个特显点进行估计．对配准后主辅图像上的同名特显点进行200倍插值，将峰值间的距离向和方位向偏移量作为该特显点处的配准误差．图4给出了二阶多项式方法和文中方法配准误差均方根值随控制点数量的变化．由图4可知，二阶多项式方法配准误差随控制点数量变化较大，而文中方法配准误差随控制点数量变化较小．例如，当控制点数量为10个时，二阶多项式方法的方位向和距离向配准误差分别为0.79像素和0.86像素，文中方法的方位向和距离向配准误差分别为0.05像素和0.07像素; 当控制点数量为100个时，二阶多项式方法的方位向和距离向配准误差分别为0.05像素和0.12像素，文中方法的方位向和距离向配准误差分别为0.04像素和0.05像素．由表1可知，主辅图像间的垂直基线大于水平基线，导致目标高程对距离向配准偏移量的影响大于方位向配准偏移量，所以二阶多项式方法的距离向配准误差大于方位向配准误差．文中方法可以追踪配准偏移量随目标高程的变化，所以文中方法的方位向和距离向配准误差近似相同，并且可以满足0.1像素的配准精度要求．上述结果验证了文中方法的精确性和稳健性．

图3 干涉相位图图4 二阶多项式方法和文中方法配准误差统计结果

为了直观证明文中方法相对于二阶多项式方法的精确性和稳健性，图5给出了控制点数量为10个时，100次试验的12个特显点处的配准误差分布．图5(a)～(d)分别为二阶多项式方法的方位向和距离向配准误差、文中方法的方位向和距离向配准误差．由图5可知，二阶多项式方法的方位向和距离向配准误差最大值分别为3.47像素和3.45像素，而文中方法的方位向和距离向配准误差最大值分别为0.17像素和0.30像素．

图5　二阶多项式方法和文中方法配准误差分布

3　结 束 语

针对长基线或复杂场景星载InSAR图像配准问题，文中提出了地形高程自适应的降维图像配准方法．首先，给出了主辅图像间的配准偏移量计算模型；然后，提出了包含地形高程项的二维三元一阶多项式模型．该模型可以追踪配准偏移量随地形高程的变化，提高了配准偏移量拟合精度；最后，利用主辅图像的成像参数和轨道数据计算多项式一次项系数，将需要利用控制点处配准偏移量估计的系数减少为两个常数项系数，提高了系数估计的稳健性．实测数据处理结果验证了文中方法的精确性和稳健性．

参考文献：

[1] YAGUE-MARTINEZ N, PRATS-IRAOLA P, GONZALEZ F R, et al. Interferometric Processing of Sentinel-1 TOPS Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(4): 2220-2234.

[2]LIANG C, FIELDING E J. Interferometry with ALOS-2 Full-aperture ScanSAR Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(5): 2739-2750.

[3]于瀚雯, 保铮. 利用L1范数的多基线InSAR相位解缠绕技术[J]. 西安电子科技大学学报, 2013, 40(4): 37-41.

YU Hanwen, BAO Zheng.L1-norm Method for Multi-baseline InSAR Phase Unwarpping[J]. Journal of Xidian University, 2013, 40(4): 37-41.

[4]靳峰, 冯大正. 一种快速准确的图像配准算法[J]. 西安电子科技大学学报, 2015, 42(6): 88-93.

JIN Feng, FENG Dazheng. Fast and Precise Image Registration Algorithm[J]. Journal of Xidian University, 2015, 42(6): 88-93.

[5]HAN Y, BOVOLO F, BRUZZONE L. Segmentation-based Fine Registration of Very High Resolution Multitemporal Images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(5): 2884-2897.

[6]KRIEGER G, MOREIRA A, FIEDLER H, et al. TanDEM-X: a Satellite Formation for High-resolution SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(11): 3317-3341.

[7]STONE H S, ORCHARD M T, CHANG E C, et al. A Fast Direct Fourier-based Algorithm for Subpixel Registration of Images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(10): 2235-2243.

[8]薛海伟, 冯大正. 一种快速的InSAR图像精配准方法[J]. 西安电子科技大学学报, 2016, 43(3): 172-178.

XUE Haiwei, FENG Dazheng. Fast Subpixel Registration Method for InSAR Images[J]. Journal of Xidian University, 2016, 43(3): 172-178.

[9]SUN K, LIU L M, TAO W B. Image Matching via Feature Fusion and Coherent Constraint[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(3): 289-293.

[10]FAN J W, WU Y, WANG F, et al. New Point Matching Algorithm Using Sparse Representation of Image Patch Feature for SAR Image Registration[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(3): 1498-1510.

[11]NITTI D O, HANSSEN R F, REFICE A, et al. Impact of DEM-assisted Coregistration on High-resolution SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(3): 1127-1143.

[12]SANSOSTI E, BERARDINO P, MANUNTA M, et al. Geometrical SAR Image Registration[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(10): 2861-2870.

[13]刘艳阳, 李真芳, 杨娟娟, 等. 分布式卫星InSAR目标定位近似闭式解[J]. 西安电子科技大学学报, 2012, 39(4): 87-93.

LIU Yanyang, LI Zhenfang, YANG Juanjuan, et al. Quasi-closed-form Solution for Distributed Satellites InSAR Geolocation[J]. Journal of Xidian University, 2012, 39(4): 87-93.

[14]李锦伟, 李真芳, 侯英龙, 等. 基于递推公式的星载SAR高效正向定位算法[J]. 电子与信息学报, 2014, 36(2): 409-414.

LI Jinwei, LI Zhenfang, HOU Yinglong, et al. A Novel Efficient Spaceborne SAR Geolocation Method Based on Recursion Formulae[J]. Journal of Electronics & Information Technology, 2014, 36(2): 409-414.