面向震害信息提取的多源遥感图像自动配准

杨德贺，陈宜金，杨 溪，吕京国，张子昕，张 帅

(1.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083;2.北京建筑工程学院，北京 100044)

0 引言

多源遥感图像已应用于地震灾害监测、国土资源勘查和运动目标识别等各种领域，而图像配准是遥感应用的重要基础。单一传感器图像的几何配准就是通过寻找2 景图像的相同或相似区域的偏移量，对待配准图像进行重采样，受噪声和灰度变化的影响比较大。多传感器获取的图像具有互补性和冗余性，包含了更加丰富的信息，可反映地物特点的全面性;因此，多源遥感图像的几何配准是利用相似性测度计算图像间的变换参数，将同一区域的多景图像变换到同一坐标系统下，实现像元的最佳匹配。将配准后的图像应用于地震等级评估中，可解决地震救援和震害评估对基础资料需求的问题。地震灾害信息提取是指利用多源遥感图像的光谱信息或纹理特征判断地震灾害信息变化程度。

在特征定位、变化检测和立体视觉等不同应用中，需要进行亚像元级的图像配准，而像元级的配准算法(例如:空域互相关配准［1］、傅里叶变换配准［2］等)很难满足亚像元级图像配准的高精度定位要求。目前，学术界倾向于研究亚像元级配准方法，包括插值方法［3］、微分法［4］和最大互信息法［5］等。插值法取决于插值的质量;微分法对光照的变化敏感;而最大互信息法则需要人工干预选取特征点。上述方法均不适合多源遥感图像的自动化配准。

针对此问题，本文提出了多源遥感图像自动配准算法和分类流程，由以下几个步骤组成:①提取影像特征点以获取同名点，采用随机抽样一致性(random sample consensus，RANSAC)算法剔除误匹配点;②通过对数极坐标变换，求解旋转和缩放因子;③采用频域相位相关算法计算图像间的偏移量，利用插值算法对待配准图像重采样，实现待配准图像与基准图像的几何配准;④利用地震前图像和与之配准的地震后图像进行重点区域特征提取;⑤面向特征信息提取结果，开展基于特征的震害信息提取。

1 算法

1.1 Moravec 特征点提取

Moravec 是点特征提取的经典算子［6］，在以每一个像元为中心的影像窗口(w，w)中，计算对角线方向4个相邻像元灰度差的平方和作为感兴趣值;给定一个经验阈值t(一般取整景图像的灰度均值)，将感兴趣值大于阈值的点作为候选点;选取一定大小的窗口，从候选点中选择一个最大的感兴趣值作为特征点。该算子具有计算量小、不丢失灰度信息等优点，由(w，w)和t 决定特征点的数量。

RANSAC 算法是一个以迭代法来估计参数的数学模型［7］，应用于基于特征的图像配准中，具有较高的剔除误匹配点的能力。首先，它依据设定的容许误差，将所有的初步匹配点对分为内点和外点;然后，采用比较准确的内点来进行参数估计，剔除误匹配点。该算法具有较高的可靠性、配准精度和鲁棒性等优点。点匹配与误匹配点剔除的技术流程如图1 所示。

图1 特征点匹配流程Fig.1 Flow chart of matching process for feature points

以地震前全色波段遥感图像为基准图像、地震后高分辨率数码图像为待配准图像，采用Moravec特征提取算子进行特征点提取(图2)。

图2 Moravec 特征点提取效果Fig.2 Extraction of feature points with Moravec

1.2 对数极坐标变换

对数极坐标变换是从笛卡尔坐标系变换到对数极坐标系，将图像匹配中笛卡尔坐标系下的旋转和缩放变换转化为对数极坐标下的平移变换［8］。笛卡尔坐标系可表示为

式中:z 为笛卡尔坐标空间;(x，y)为点坐标;ρ 为点的极径;θ 为点的极角。

对数极坐标系可表示为

映射变换将图像在轴向缩放和旋转上的变化转变为对数极坐标上的水平和垂直方向位移，从而实现从一个圆形区域的图像映射成一个矩形区域的图像，即为对数极坐标的二维不变性。缩放和旋转可表示为

式中:Δρ 为放大倍数;Δθ 为旋转角度。

水平和垂直方向上的缩放和旋转表示为

图3 示出对数极坐标变换的结果。

图3 对数极坐标变换Fig.3 Logarithmic polar coordinate transformation

1.3 相位相关与插值

首先，利用相位相关和插值算法对2 景图像进行插值，以提高图像分辨率［9］;然后，进行傅里叶变换，利用频域相位相关的多相位分解求出图像之间的位移关系，则可得到图像之间的亚像元级配准结果。

2 景图像存在的水平和竖直方向的平移矢量，通过插值放大一定的倍数，以解决频域相位相关非整数不能求解的问题。对2 景插值前图像的归一化互功率谱进行傅里叶逆变换，得到单位脉冲响应函数。它采用二维sinc 函数来近似，通过解求sinc 函数，求得相位相关峰值处坐标，即可解出平移矢量。常用的插值算法有最邻近插值法、双线性插值法和三线性插值法。图4 示出一维sinc 函数。

图4 sinc 二维空间(左)和三维空间显示(右)Fig.4 Display of sinc in 2D (left)and 3D(right)

1.4 基于特征的图像分类

应用于地物分类的图像含有不同类型的地物，需要选择一定的特征代表不同类型的地物;以一定的方法将特征所对应的空间划分为不同的子空间，使图像中代表各类地物的像元归为相对应的子空间。因此，在进行遥感图像分类时，既要考虑影像的灰度信息，也要顾及影像的结构信息，故需要选择不同的影像特征参数来反映这些要求。目前，常用的影像特征参数包括均值、方差、能量、对比度、熵、相关性和均质度等［10］。采用光谱信息与结构信息等多种特征相结合的方法，可以充分利用影像灰度及其分布的信息，改善图像分类的效果。

2 实验

2.1 数据源

实验数据为四川省玉树地区2009年5月Quick-Bird 全色波段遥感图像和2010年6月ADS40 机载高分辨率数字图像。选取其中1个波段的图像进行几何配准实验，配准精度为0.1个像元;并进行灰度变化和抗噪性实验，测试两者对配准结果的影响。实验证明，基于傅里叶变换的配准方法具有较高的鲁棒性和较强的抗噪性。在旋转和缩放的情况下，本文提出的自动配准方法具有较高的稳定性。分类后的变化检测实验表明，稳定的几何配准对信息变化检测有较高的鲁棒性。图5 示出本次实验的流程。

图5 实验流程Fig.5 Flow chart of experiment

2.2 配准实验

以震前QuickBird 全色图像为基准图像(图6(a))，以震后ADS40 航摄图像为待配准图像(图6(b))，使用本文提出的自动配准方法进行多源遥感图像的几何配准实验，获得与QuickBird 全色图像配准的ADS40 航摄图像(图7)，用于配准的特征点均方根误差RMSE=0.96 像元(表1)。

图6 地震前QuickBird(左)和地震后ADS40(右)图像Fig.6 QuickBird image before earthquake (left)and ADS40 image after earthquake (right)

图7 配准后的ADS40 图像Fig.7 ADS40 image after registration

表1 部分特征点误差Tab.1 Error of some feature points (像元)

不同时相的多源遥感图像，因受大气、太阳入射角、视角及地形的影响，导致地震前后影像灰度变化大，影响到特征匹配的效果，所以在几何配准前要对地震后图像进行相对辐射校正(即辐射归一化)。经灰度调整后，采用Moravec 对2 景图像提取匹配的特征点;利用RANSAC 剔除误匹配点对，对匹配的特征点对进行对数极坐标变换;解求变换点对的旋转角和缩放系数，进行相位相关求解偏移量;采用双线性插值重采样完成配准过程。

2.3 灰度变化和噪声对配准的影响

通过灰度线性变化函数D=aT+b(其中a，b 为灰度变化因子)，使原来的灰度值T 变为D。由于噪声对图像的配准存在影响，故利用衡量噪声强弱的信噪比的对数函数来加入噪声。图8 示出灰度变化后和加入椒盐噪声后的图像。

图8 灰度变化后(左)和加入椒盐噪声后(右)图像Fig.8 Images after gray change(left)and adding noise(right)

从表2 可以看出，灰度变化和噪声对配准误差的影响都比较小。

表2 灰度变化和噪声对配准误差的影响Tab.2 Influence of gray change and noise to registration error

2.4 配准误差对分类精度的影响

对自动配准后的2 期遥感图像进行了分类比较。结果表明，当配准误差较大时，分类面积的平均变化率会较大，但也会稳定在一定范围;当配准误差较小时，分类面积的平均变化率较小，分类精度较高。表3 示出不同配准误差下的地物分类变化率。

表3 不同配准误差下的分类变化率Tab.3 Change rate of classification under different registration errors

2.5 震害信息提取

本文采用分类后处理的方式提取地震灾害信息。针对所提取的重点目标区域，以多特征统计向量为分类标准。首先对地震前后的重点区域分别进行基于特征的分类，然后根据2 期图像的分类结果进行变化检测(图9)。

图9 震害信息提取Fig.9 Extraction of earthquick damage informatiom

震害信息的提取结果详见表4。

表4 震害信息提取结果Tab.4 Extraction results of earthquick damage information (像元)

从图9 和表4 可以看出，地震前、后房屋、裸地和道路3 种类别的区分明显，分类效果较好;地震后的总体分类精度为68%，说明地震后的分类结果与地震前的实际类型相差较大;地震后的裸地类别变化相对较小，地震后道路和房屋变化相对较大。尽管研究区域中的房屋损毁情况比较大，但是对房屋的分类较为准确。震害信息变化图(图9(c))表明，对多源遥感图像进行几何配准、分类后，能够检测出房屋、裸地、道路等不同类别的损毁程度，可以为震害损毁分析与评价提供重要信息。

3 结论

本文提出了多源遥感图像自动配准的算法，并通过震害信息检测对配准的效果进行了验证。结果表明:

1)解决基准图像和待配准图像的旋转、缩放以及配准自动化问题是当前多源遥感图像几何配准的难点。几何配准精度的高低也决定了图像分类效果的好坏，并会对变化检测的结果产生影响。

2)本文提出的面向地震灾害信息提取的多源遥感图像自动配准算法，经遥感图像数据的几何配准实验证明，该算法综合了Moravec 算子、RANSAC算子、对数极坐标变换和相位相关等算法的优点，实现了多源遥感像亚像元级的自动配准，较好地解决了灰度变化和噪声对几何配准的影响以及基准图像和待配准图像的旋转和缩放等问题。

3)将用本文方法配准后的震前和震后图像应用于地震灾害信息提取，震害类型变化信息统计结果表明，较高的几何配准精度提高了震害的分类精度，说明本文提出的高精度自动配准方法和基于特征的分类可以较好地应用于地震震害评估、海洋水色反演等领域的实时变化检测中。

本文的不足之处为:①因获取遥感图像的时相及传感器等方面的不同，会出现实际未变化地物却存在灰度变化的情况，所以本文的几何配准方法研究应采用相对辐射校正处理后的遥感图像;②本文涉及的算法比较复杂，下一步工作将研究如何运用并行计算方法提高计算效率。

［1］吴金鹿，陈绍炜，孙 磊.一种基于互相关匹配的空域目标提取方法［J］.科学技术与工程，2011，11(32):7948－7951.Wu J L，Chen S W，Sun L.An airspace object extraction method based on cross－correlation matching［J］.Science Technology and Engineering，2011，11(32):7948－7951.

［2］林 茂.基于改进的曲线傅里叶变换图像配准研究［J］.计算机仿真，2011，28(10):265－268.Lin M.Improved digital image contour curve Fourier descriptor algorithm［J］.Computer Simulation，2011，28(10):265－268.

［3］詹 毅.基于泰勒展开式的图像插值方法［J］.计算机工程，2012，38(13):202－204.Zhan Y.Image interpolation method based on Tyler expansion［J］.Computer Engineering，2012，38(13):202－204.

［4］张晓芬.基于偏微分方程的图像配准技术研究［D］.河北:华北电力大学，2009:32－39.Zhang X F.Research on image registration technique based on PDE’s［D］.Hebei:North China Electric Power University，2009:32－39.

［5］程有娥.基于角点特征和最大互信息的图像配准［J］.计算机系统应用，2012，21(6):186－190.Cheng Y E.Image registration based on mutual information and corner points［J］.Computer Systems and Applications，2012，21(6):186－190.

［6］刘晓龙，李英成.基于多源遥感影像融合的影像匹配技术［J］.测绘科学，2007，32(3):59－61.Liu X L，Li Y C.Image matching based on multi－source remotesensing image fusion［J］.Science of Surveying and Mapping，2007，32(3):59－61.

［7］周剑军，欧阳宁，张 彤，等.基于RANSAC 的图像拼接方法［J］.计算机工程与设计，2009，30(24):5692－5694.Zhou J J，Ouyang N，Zhang T，et al.Image mosaic method based on RANSAC［J］.Computer Engineering and Design，2009，30(24):5692－5694.

［8］雷 凯，刘艳滢，王延杰，等.基于特征点的对数极坐标变换图像配准算法［J］.光学技术，2007，32(3):435－437.Lei K，Liu Y Y，Wang Y J，et al.Log polar transformation based on feature points for image registration［J］.Optical Technique，2007，32(3):435－437.

［9］刘卫光，崔江涛，周利华.插值和相位相关的图像亚像素配准方法［J］.计算机辅助设计与图形学学报，2005，17(6):1273－1277.Liu W G，Cui J T，Zhou L H.Subpixel registration based on interpolation and extension phase correlation［J］.Journal of Computer－aided Design and Computer Graphics，2005，17(6):1273－1277.

［10］贾永红.数字图像处理［M］.武汉:武汉大学出版社，2003:182－186.Jia Y H.Digital image processing［M］.Wuhan:Wuhan University Press，2003:182－186.