由于数据采集时间、大气干扰、光谱组合、传感器平台高度等差异和系统因素，多源遥感影像间在使用严格成像模型精确定位后仍然存在相对误差。因此，在检测、识别多源遥感影像所包含的地物信息等任务中，首要任务之一就是通过影像配准技术解决这个问题。目前，根据特征提取和描述算法的不同，国际上主流的影像配准方法主要分为基于灰度和基于特征的配准两类。例如，文献^[1]提出了具有尺度和旋转不变性的SIFT特征检测算子；文献^{[2, 3]}提出了具有尺度旋转不变性的点要素提取和特征描述算子SURF、ORB；文献^[4]将相关系数应用到多时相、多传感器遥感影像的配准中，使配准精度达到子像素级别；文献^{[5, 6]}通过取得互信息测度极大值的变形改正位置，解决了光谱特征差异显著的多源异构遥感影像的配准。但由于以上各种算法的局限性，现在仍然没有一个通用的算法能解决多尺度、多时相、多姿态、多光谱、多传感器等情况下的影像旋转、尺度畸变、辐射畸变、平移等问题，且存在效率与效果之间的矛盾^{[7, 8, 9, 10]}。本文通过对特征提取方法、特征匹配和约束策略、几何模型参数估计与粗差采样等环节算法的综合比较与研究，设计了一种利用相位一致性的尺度不变特征变换(phase congruency based on scale invitation feature transform，PC-SIFT)的分层逐级的多源遥感影像配准方案，以解决多源光学卫星影像配准过程中由辐射和几何畸变引起的同名特征提取不足、特征描述不稳定、匹配错误率高、配准精度低等问题。

本文研究主要包括4个部分：比较现有特征提取方法在应对几何畸变、辐射畸变及噪声干扰方面的能力，提出改进特征向量描述的方法；比较常见的向量型和区域型特征相似性检验方法，寻找适合PC-SIFT特征描述向量的相似性检验方法；寻找改进影像间存在的严重辐射畸变时，能有效提高PC-SIFT特征匹配结果正确率的约束策略；试验不同的参数求解与估计方法，选取合适的最优模型选择策略。

在多源遥感影像配准中，不仅需要多影像信息的参与和支持，而且由于构成SIFT算子特征描述符的核心要素——梯度的强度变化线性无关，导致从约束条件或特征空间着手的SIFT改进算法，不能从根本上解决特征描述对辐射非线性畸变抵抗力差的问题。因此，本文利用相位一致性对影像纹理描述稳定的特点，引入相位一致性代替SIFT中的梯度用于影像结构描述，提出了一种抗影像辐射畸变的特征描述新方法——PC-SIFT。PC-SIFT的具体改造过程如下。

1) PC-SIFT主方向确定。包括以下几个关键过程：① 通过搜索高斯差分影像(DOG)像素极值得到所有SIFT特征点,按尺度值进行统计，找出频率最高的尺度 σ及其对应的高斯滤波影像G_σ。② 根据所有特征点按尺度σ_i与σ的比值关系，换算得到在影像G_σ中的等效窗口S_i，其计算公式为S_i=[log₂(σ_i/σ)+1]×R\5C_{ori_sig}\5σ。同理，像素中心位置(x,y)也按上述比例关系换算为影像G_σ上的等效位置(x_i,y_i)。R为常数，为保证用于梯度直方统计的窗口足够大，通常取3.0，C_{ori_sig}代表对窗口内像素值进行距离加权运算的高斯核标准差比例因子。③ 对影像G_σ分不同方向按式(1)执行全局相位一致性计算，其中,Log-Gabor滤波器的方向角与统计采样角间距Δ_ori保持一致，最终得到由k=π/Δ_ori(取整)层相位一致性特征影像构成的集合C_PC。④ 以②和③的结果为基础，对所有任一特征点P(x,y)逐一做基于相位一致性特征值的高斯距离加权统计。

2) PC-SIFT特征描述向量构造。包括以下几个关键过程:① 与上述主方向确定中的①、②基本相同，仅需按照公式F_i=[log₂(σ_i/σ)+1]×F对原始SIFT统计窗口大小F按尺度比值关系做出换算，其中,F=2^1/2(D+1)\5W_Bin=2^1/2(D+1)\5(C_{des_fctr}\5σ)。C_{des_fctr}是常数，WBin是基本统计单元子窗口尺寸对应的像素单位长度，D是影像行(列)方向上规则排列的基本统计单元个数。② 利用上述主方向确定阶段得到的集合C_PC，对任意特征点P，在以其位置为中心、F_i为尺寸大小的窗口内，对C_PC上的每一个像素Q(m,n)，按照P的主方向角φ和统计采样角间距Δ_dsc，计算得到与像素梯度值方向角等效的θ角(θ∈O={φ,φ+Δ_dsc ,φ+2Δ_dsc…φ+(N-1)Δ_dsc}，N是统计空间的维数，对θ大于或等于π进行反向运算，保证θ∈[0,π)，O为当前统计方向)。③ 根据集合O中的θ值，找出C_PC中Log-Gabor滤波方向角与θ最近的三层影像，然后使用对应位置Q(m,n)上的三个计算值进行二次曲线拟合。

通过以上改进后，按式(2)计算某一像素Q(m,n) 对应于特征向量的位置及贡献值，并按上述过程循环遍历邻域窗口内的所有像素，实现PC-SIFT特征描述构造。

式中，W_o(x,y)为频率扩展权重因子；A_no(x,y)和Φ_no(x,y)为图像在给定滤波器尺度n和方向o在点(x,y)的振幅和相位;ΔΦ(x,y)为加权的平均相位; 符号为当正值时取本身，否则取0；T为噪声阈值;ε是避免除零的常数。式(2)中，由于特征向量 V 的计算公式中加入了主方向参照，同时考虑了窗口大小的尺度改变，很大程度消除了旋转形变和尺度畸变的影响，使得算法在抵抗仿射变形方面具有极高的稳定性。

本文针对PC-SIFT特征提取算法，选取最近距离比值法作为PC-SIFT特征描述向量的相似性检验方法，以其特征描述的尺度和主方向为主要关注对象，检验匹配结果集中所有特征点对的尺度比例和主方向旋转相对关系的正确性，并将影像频率域相位信息作为特征初始匹配结果正确性判断的二次约束策略。

特征匹配和约束策略的几个关键过程为： ① 使用最近欧氏距离测度值 ω₀和第二欧氏距离测度值ω₁的比例r=ω₀/ω₁作为新的测度，来检测出所有在信息表达上一致或相似的特征要素组合,生成初始匹配集M_i。② 在采用较小阈值作为初始匹配集测度判别条件的情况下，对M_i利用尺度和旋转约束，进一步获取正确匹配集。③在②的基础上对匹配集进行平移约束，且为进一步改善参数估计的精度水平，使用互功率谱影像上最强冲击响应位置附近的强度值建立与像素位置相关的函数模型进行极值内插。

本文为提高模型参数估计和粗差检校后模型选择的适用性和参数估计阶段的自动化水平，在采用RANSAC方法进行粗差剔除和模型优化运算后，利用相容率、覆盖率和整体离散度三个几何畸变模型评价指标决定最优目标函数的选择。

1) 相容率 (R_cst)。由集合C_in的样本数n_in与集合C的样本数n进行比值运算获得，其中,C为匹配同名点坐标集，C_in为对C采用RANSAC运算后得到的优化后模型参数具有一致性表达的正确同名点集合。R_cst值越大，表示对应的目标函数M对集合C的拟合程度越好。

2) 覆盖率(R_cvr)。由集合Cin所有样本的空间占有范围Ain与集合C所有样本的空间占有范围A进行比值运算获得，空间占有范围是指二维平面内包围所有参考样本点的最小凸多边形的面积。R_cvr值越大，表示相对应的目标函数M对影像的整体控制性越高。图 1 (a)展示了某点集经RANSAC计算后A_in和A的关系示意，图 1(a)中，内多边形为A_in，外多边形为A。

图 1 覆盖率比值范围和整体离散度计算TIN结构图 Fig. 1 Diagram of Coverage Ratio Scope and Discrete Degree Overall Computing with TIN Structure

图选项

3) 整体离散度(D_dsc)。由集合C_in中所有参考样本点彼此之间的欧氏距离平均值计算获得，而某一样本点与其邻近样本点之间的邻近关系由 集合C_in中所有内点参考样本点构建的不规则三 角网TIN来确定，统计组成TIN的所有边长取大小，D_dsc可由式(3)表示：

其中，N表示组成TIN的边的个数;L_diag代表包围TIN最小矩形的对角线长度(像素单位)；E_i为任意两点确定的欧氏距离； 为所有E_i的统计均值。图 1 (b)展示了用于整体离散度计算的TIN结构。

本文以简单耦合的方式将三个指标联系在一起，且通过上述指标定义可以发现，当R_cst和R_cvr取值越大、D_dsc取值越小时，相应模型对数据集描述的适用性越好。依据这一关系，可构造公式R= (R_cst×R_cvr)/ D_dsc作为评价模型适用性的最终指标。通过RANSAC的迭代运算，不断记录最大R值与最优模型M_i，直到模型集中所有模型参与运算并完成评价值R的比较，输出最大R值对应的模型类型、参数和最大一致集，完成最优模型选择。

本文针对多源遥感光学卫星影像的配准问题，在分析和改进影像配准中的特征提取与描述、匹配和参数估计等关键步骤的基础上，提出了如图 2所示的基于PC-SIFT的多源遥感影像配准总体实施方案。

图 2 配准方案 Fig. 2 Registration Scheme

图选项

表 1为实验数据基本信息，为体现影像在辐射方面的差异，避免光谱组成上的线性相关性，目标影像和参考影像均选自光谱采样范围不同且非邻近波段。实验影像均为使用传感器成像模型初级几何约束或定标后产品，影像像素灰度的量化级别均为0~255，目标和参考影像间的覆盖率水平平均为80%左右，目标影像的投影和坐标系统未知。

本文采用的硬件环境为：Inter Core Duo CPU 2.53 GHz，内存4.0 G，硬盘60 G；软件环境为：开发和运行系统平台为Windows 7 Professional 32位操作系统，软件设计和开发环境Visual Studio 2010，开发语言VC++ 2010；外部辅助功能函数和数据接口库包括OpenCV 2.3、GDAL1.7和GSL。本文实验过程中，主要阶段各项参数配置为：① 影像准备阶段。影像金字塔的采样像素比例K=2，层级数N的大小根据影像图幅范围和K值的计算关系自动确定，且应保证顶层金字塔影像任一维(列或行方向)的长度值(像素单位)不小于500。② 粗匹配阶段。DOG尺度空间阶数为3，每一阶中高斯金字塔影像层数为6，初始高斯卷积核标准差为1.6，检测特征点对比度阈值为0.2，特征向量维数为128，Log-Gabor滤波方向数为18，倍频数为1，最小波长为3个像素单位长度，径向高斯卷积核标准差大小为0.55，角向角分辨率与高斯卷积核标准差比例系数为1.2；用于特征向量相似性检验的最近距离比阈值为0.65；相位相关运算窗口大小为256×256，相对旋转角评价阈值为2°，尺度比评价阈值为1.044 3，相对平移量评价阈值为12(单位为像素长度)。③ 模型参数估计与粗差检校阶段。初始内点正确率为0.25，随机抽取样本终止条件阈值为0.01，一致集搜索和判定的总体均方根误差阈值为2.0(单位为像素)；用于自适应选择的模型集合包括仿射形变、二次多项式、透视形变和二次投影形变模型。④ 逐级匹配阶段。梯度相关矩阵计算的窗口大小为7，Sobel扩展核大小为3，特征自由检测系数为0.04；同名特征点搜索半径为30×(2N-1-i)，其中,i为金字塔影像层级序数；最小二乘单点匹配运算窗口大小为11，最大迭代运算次数为8，迭代终止最大中误差阈值为1.0，同名点判定的相关系数阈值为0.80。⑤ 最终模型参数估计与粗差检校阶段。执行过程与第③阶段类似，RANSAC方法配置参数与③ 阶段保持相同，由于观测量变为影像像素灰度，因此,最终用于一致集判定的总体均方根误差阈值为10.0(单位为灰度级)；排除二次投影形变模型，用于自适应选择的模型降为3种。⑥ 重采样阶段。选择双线性内插法来重采样目标影像灰度值。

2.3 结果分析

图 3为4组实验数据在未进行几何校正重采样前，经过粗匹配和分级匹配处理运算后得到的同名点分布、局部卷帘与配准结果全局套合效果图。其中，第一组实验数据中，采用可见光蓝色波段目标影像与近红外波段参考影像，两幅影像上地物不仅在辐射变化上趋势各异，而且不具有线性相关，影像表现为全局刚性形变，同时，同名区域的初始对应关系未知，难以直接使用区域匹配方法。利用PC-SIFT结合相位约束的方法一定程度上解决了辐射畸变条件下特征检测重复率低、描述稳定性差的问题，为初始几何形变关系的预测提供了具有较高正确率的匹配同名点坐标观测集合。匹配得到433个控制点，分布较为均匀，最终校正模型自适应选择为二次多项式。从实验效果图可以看出，河道和陆地交界处较陆地不稳定，但匹配影像整体和局部几乎没有偏移和形变。影像配准总体均方根误差(RMSE)在1.27个像素左右，单点最大中误差为1.93个像素。

图 3 4组影像数据实验效果图 Fig. 3 Effect of Four-Group Experiments

图选项

第二组实验数据中，采用近红外通道目标影像和可见光绿光波段参考影像，在灰度值上与第一组数据有相同的变化形式，且实验影像投影方式未知，并且，由于参考影像产品自带6个仿射变形参数，导致影像间除在全局存在仿射变形外，还存在因投影差异引起的未知几何形变。通过PC-SIFT方法配准，匹配得到298对特征点，空间分布相对均匀，最终影像校正模型自适应选择为二次多项式。从实验效果图可以看出，除河道泥沙堆积小岛轮廓因江水冲刷引起形变和桥梁显示出1~2个像素偏移外，整体和局部配准后影像对齐较为准确。影像配准RMSE在1.69个像素左右，单点最大中误差1.97个像素。第三组实验数据中，采用红外波段目标影像和全色参考影像，空间分辨率相差最大，且由于参考影像包含了整个可见光范围，波段相关性较其他组实验数据高，影像在全局范围内反映为平面刚性几何形变。通过PC-SIFT方法配准，匹配得到控制点225对，分布较其他实验组稀疏，最终影像校正模型自适应选择为透视变形，这也反映出影像上部地势起伏明显而下部地势平坦引起了估计参数的小范围调整。从实验效果图可以看出，雪线沿山脊变化趋势和位置高度一致，且人工地物也没有因模型类型发生扭曲或变形，影像叠合度较高。影像配准RMSE在1.13个像素左右，单点最大中误差为1.81个像素。第四组实验数据中，采用近红外通道目标影像和可见光蓝色波段参考影像，其辐射差异与第一、二组数据类似，但因影像分辨率很高，高程引起的几何畸变与遮挡影响较大，不同地物表现的光谱曲线变化趋势与形式各异，导致实验数据几何与辐射畸变形式复杂。通过PC-SIFT方法配准，匹配得到控制点641对，与其他三组实验相比分布均匀性最差，最终影像校正模型自适应选择为仿射模型。从实验效果图可以看出，除桥梁等地物有微小视差外，影像整体间对齐较好。实验结果表明，本文方案对二级产品的中分辨率遥感影像配准的平面精度在1.8个像素以内，不考虑地形起伏和人工地物高程影响，对二级产品的高分辨率遥感影像配准的平面精度在1.5个像素内。但在高分辨率影像中，因为没有三维高程观测量辅助配准校正，在配准后的结果影像和参考影像中选取人工检查点计算残差会出现较大困难，且二维几何形变关系模型对高分辨率影像的局部控制有限，所以本文提供的配准方案中的特征点匹配方法适用于高分辨率影像，而估计模型参数的过程需要采用三维几何形变模型和高程观测量。

为进一步验证本文所提PC-SIFT配准方案的精度和准确性，将上述实验结果与传统遥感影像处理软件ENVI配准模块的输出结果进行对比，ENVI配准模块采用基于局部网格划分与区域搜索的最小二乘单点匹配方法，粗差检校和校正模型选择两个过程需要人工参与判断。为保证该模块配准结果与本文实验结果在容许残差精度上尽可能一致，选择实验中自适应选择结果相同的模型进行参数估计，并采用手工排除最大中误差的方式将单控制点中误差精度控制在2个像素内。从表 2可以看出，两种方案在最终控制点总体匹配精度上水平相似，说明本文配准方案与常规基于灰度匹配的配准方法具有相一致的精度；但在控制点获取总数上，除第四组实验外，ENVI配准模块高于本文方案，这也说明影像纹理信息对基于特征匹配的方法十分重要，而灰度匹配的方法仅与格网大小和同名点搜索区域大小相关。

本文针对多源遥感光学卫星影像配准中的问题，首先提出了一种基于影像局部相位分析、兼具几何和辐射畸变抵抗性的特征提取方法；接着，针对该方法提出了以相位相关理论估计局部几何形变参数为基础的SIFT型特征匹配策略；然后，基于随机采样一致性粗差检校方法，提出了一种在参数估计过程中校正模型的自适应选择方法； 最后综合以上研究成果，在特征匹配和灰度匹配相结合、以影像金字塔结构分级匹配策略为指导的基础上，设计实现了一套完整、流程化的多源遥感光学卫星影像自动配准处理与实施方案，并分析评价了配准后的影像精度，验证了本文所提配准策略的稳定性和可靠性。