城市绿地具有重要的生态服务功能，相当于自然调节器，通过一系列的生态效应，综合调节着城市环境^[1-2]，能够吸收各种有害气体，起到净化环境的作用，还能够有效缓解城市热岛效应. 然而，随着经济的快速发展和城市化进程的加快，城市重复性建设、人为踩踏、非法占用、养护不当、尾气排放等因素造成城市绿地破坏，致使城市绿地生态结构单一、生态服务功能降低，从而导致城市环境问题日益严重^[3]，迫切需要采取科学有效的城市绿地管理措施. 提取城市绿地分布及其变化信息对城市绿地变化监测具有重要意义，可为城市绿地管理提供可靠的基础数据，有利于维护城市环境，促进城市可持续发展.

现有的航天、航空和地面遥感技术已广泛运用于提取城市绿地分布及其变化信息^[4-7]. 在遥感图像数据的获取过程中，城市地面的起伏变化和航拍视角变化造成城市绿地遥感图像间产生非刚性畸变、重叠度低等问题，从而引起同一场景遥感图像（如：多时相、多视角）不在同一坐标系，不能直接进行变化检测. 遥感图像配准可以通过2幅图像的特征点对应关系，将2幅图像校正到同一坐标系，是后续变化检测的前提条件.

遥感图像配准被广泛运用于提取城市绿地分布及其变化信息的关键技术，当前的遥感图像配准方法主要包括：1）基于区域的配准方法；2）基于特征的配准方法^[8]. 其中，基于区域的配准方法主要是利用图像的灰度信息，建立2幅图像之间的相似性度量，再统计相似性度量值最大或最小变换模型的参数值，以达到配准图像的目的；而基于特征的图像配准方法是通过提取图像中的某些特征，如点、线、面等，将这些特征与被匹配的图像进行比较分析，从而得到匹配结果. 相比前者，基于特征的方法不受光照、旋转影响，计算信息量少且效率较高. 本文主要研究城市绿地遥感图像间具有非刚性畸变及重叠度低等问题的图像配准. 同时，特征点具有普遍性以及易提取性，因此，本研究从基于特征点的角度进行配准研究. 苏珊（small univalue segment assimilating nucleus，SUSAN）算子^[9]、哈里斯（Harris）算子^[10]、莫拉维克（Moravec）算子^[11]与尺度不变特征变换（scale-invariant feature transform，SIFT）算子^[12-13]等是目前常用的几种图像特征点提取算法. 其中SIFT算子不仅具有尺度、旋转、视角、光照不变性，对目标的运动、遮挡、噪声等因素也保持较好的匹配性，受到领域内学者们的青睐.

在特征点配准问题上，Myronenko等^[14]提出了著名的一致性点漂移（coherent points drift，CPD）算法，该算法从概率密度评估的角度解决特征点匹配问题，能够很好地计算中心特征点的误差和缺失. 该算法还利用快速高斯变换（fast Gaussian transform）^[15]和矩阵低秩逼近（Low-rank matrix approximation）技术^[16]，减少了算法的计算复杂度，提升了计算速度. 为进一步提高遥感图像配准精度，梁栋等^[17]利用了非下采样轮廓波变换（non-subsampled contourlet transform，NSCT）在图像分解上的灵活性和SIFT算法在特征描述上的有效性进行遥感图像配准. Jian等^[18]提出了一种基于高斯混合模型（Gaussian mixture model，GMM）的非刚性点集配准算法. Liu等^[19]提出了有限空间顺序约束（restricted spatial order constraints，RSOC）算法，该方法采用鲁棒的图匹配技术来去除错误的匹配. Zhao等^[20]提出了基于尺度不变特征变换（SIFT）和区域互信息优化的配准方法. Yang等^[21]提出了全局和局部混合距离-薄板样条（global and local mixture distance-thin plate spline，GLMDTPS）配准算法，该算法主要将全局和局部结构特征差异看作线性分配问题进行处理. 张闻宇等^[22]针对传统局部不变特征的景象匹配算法冗余点多、实时性差、抗几何变换不突出的情况，提出了基于CenSurE-star的无人机景象匹配算法.

绿地区域地形设计复杂（如：凹凸面地形、土丘、坡地）、局部布局相似（局部区域主要呈现块状、带状、楔形和混合式），以及图像获取时航拍视角变化（例如：多颗卫星不同时间航拍同一地点或无人机因飞行高度、风向、人为操控等问题引起的视角变化）造成城市绿地遥感图像间产生非刚性畸变、重叠度低等问题. 使用上述方法进行配准时会产生以下问题：1）由于绿地区域局部布局相似，仅使用基于灰度信息的算法（如SIFT算法）或特征点集的几何结构特征（如CPD、GLMDTPS等），会出现特征点辨识度不高，图像中多个区域特征点集特征相似等问题，导致图像配准失败；2）由于绿地区域地形设计复杂和航拍视角变化，加剧了遥感图像间的非刚性畸变，只能使用非刚性变换模型进行配准，然而非刚性变换模型具有较高的自由度，在配准城市绿地遥感图像过程中容易产生不适定问题（ill-posed problem），造成配准精度降低.

针对上述问题，提出一种基于图像混合特征的城市绿地遥感图像配准算法：1）为了解决仅使用基于灰度信息或几何结构特征对城市绿地遥感图像配准造成的问题，构建基于混合特征的高斯混合模型（mixture-feature Gaussian mixture model, MGMM），并将全局几何结构特征、局部几何结构特征与灰度信息特征（SIFT特征）融入到MGMM中，形成特征互补，从而得到可靠的对应关系；2）为了解决城市绿地遥感图像间非刚性畸变、重叠度低问题带来的不适定问题，在 ${{\rm L}_2}{\rm E}$ 的能量方程中引入几何结构约束项，以增强非刚性变换的稳定性.

提出的遥感图像配准算法的流程图如图1所示. 使用SIFT算法进行遥感图像特征点的提取，主要按以下步骤来介绍遥感图像的配准过程.

1）基于SIFT的特征点提取；

2）基于混合特征的SIFT特征点配准；

3）图像配准：假设，从待配准图像 ${{{I}}_ {\rm{s}}}$ 和参考图像 ${{{I}}_ {\rm{r}}}$ 中提取图像特征点集分别为 $\{ {{{x}}_{{i}}}, i = 1, 2, \cdots , n\} $ 和 $\{ {{{y}}_{{j}}}, j = 1, 2,\cdots, m\} $, ${{{x}}_{{i}}}$ 和 ${{{y}}_{{j}}}$ 分别为提取的遥感图像 ${{{I}}_ {\rm{s}}}$ 和 ${{{I}}_ {\rm{r}}}$ 上的任意特征点.

本研究使用的特征点为SIFT算法中提取的SIFT特征点. 给定待配准图像 ${{{I}}_ {\rm{s}}}$ 和参考图像 ${{{I}}_ {\rm{r}}}$，分别将提取自 ${{{I}}_ {\rm{s}}}$ 和 ${{{I}}_ {\rm{r}}}$ 的SIFT特征点集设为源点集Y和目标点集X. 此外，基于SIFT的特性，将其作为特征使用，后面将详细介绍.

特征点集配准的过程主要包括以下2个步骤：1）在待配准图像特征点集Y和参考图像特征点集X之间进行对应关系评估；2）利用1）中对应关系建立空间映射函数不断更新 ${\hat Y}$ （形变后的源点集）的位置，直到 ${\hat Y}$ 与参考图像特征点集X达到最大程度匹配.

构建能够稳健地描述遥感图像几何结构特征的混合特征描述子，并以高斯混合模型为基础，构建基于几何结构特征和灰度信息的高斯混合模型来进行对应关系评估.

1）特征差异描述子.

SIFT距离特征差异描述子：基于SIFT特征描述子的特性，使用128维的特征描述子以实现更为鲁棒的特征点配准. 将SIFT特征距离定义为

(1) ${s_{ij}} = {\left\| {{{{u}}_{{i}}} - {{{v}}_{{j}}}} \right\|^2}.$

式中：s的维度为 $m \times n$， ${{{u}}_{{i}}}$ 和 ${{{v}}_{{j}}}$ 分别为 ${{{x}}_{{i}}}$ 和 ${{{y}}_{{j}}}$ 的特征描述子.

全局几何结构特征差异描述子：使用 ${{{x}}_{{i}}}$ 和 ${{{y}}_{{j}}}$ 之间的欧氏距离，将全局几何结构特征差异描述子定义为

(2) ${g_{ij}} = {\left\| {{{{x}}_{{i}}} - {{{y}}_{{j}}}} \right\|^2}.$

局部几何结构特征差异描述子：假设 ${{{x}}_{{i}}}$ 和 ${{{y}}_{{j}}}$ 中每个点和其所在特征点集中的k个最近相邻点以距离升序构造一个局部片段，如图2所示. $N{({{{x}}_{{i}}})_k}$ 和 $N{({{{y}}_{{j}}})_k}$ 分别为 ${{{x}}_{{i}}}$ 和 ${{{y}}_{{j}}}$ 的第k个最近相邻点， $T(N{({{{x}}_{{i}}})_k}, {{{y}}_{{j}}})$ 为点 ${x_i}$ 的第k个最近相邻点 $N{({{{x}}_{{i}}})_k}$ 沿向量 $ {{{{x}}_{{i}}}{{{y}}_{{j}}}} $ 进行平移操作，并将其表示为

(3) $T(N{({{{x}}_{{i}}})_k}, {{{y}}_{{j}}}){{ = }}N{({{{x}}_{{i}}})_k}{{ + (}}{{{y}}_{{j}}} - {{{x}}_{{i}}}{{) }}.$

由此，通过计算2个局部片段中对应序列点之间的欧氏距离之和，即局部不变距离特征，得到特征点集x与y之间的局部特征差异描述子，将其定义为

(4) ${l_{ij}} = \sum\limits_{k = 1}^K\; {{{\left\| {T(N{{({{{x}}_{{i}}})}_k}, {{{y}}_{{j}}}) - N{{({{{y}}_{{j}}})}_k}} \right\|}^2}}.$

其中，K为相邻点的数量.

2）基于混合特征的高斯混合模型.

将2组点集的配准问题考虑为概率密度估计问题. 给定 ${d_{\rm e}}$ 维上的2个点集X和Y（即目标点集和源点集），将点集 ${{Y}} = \{ {y_j}|j = 1, 2,\cdots, M\} $ 中的点看作GMM中高斯成员的质心，将点集 $X = \left\{ {{x_i}|i = 1,2,} \right. $ $\left. { \cdots ,N} \right\}$ 中的点看成由GMM产生的独立同分布的数据点. 在配准过程中，使GMM的质心保持一定的邻域拓扑结构，逐渐向数据点移动，质心与数据点间的距离越近，质心和数据点成为对应点的概率越高. 将特征点集X在高斯混合模型中的概率密度函数（probability densitjy function，PDF）定义为

(5) $p({{{x}}_{{i}}}) = (1 - \xi )\sum\nolimits_{j = 1}^m {\left( {{C_{ij}}\;p({{{x}}_{{i}}}|{{{y}}_{{j}}})} \right)} + \xi/n.$

其中，

(6) $p({{{x}}_{{i}}}|{{{y}}_{{j}}}) = {(2\pi {\sigma ^2})^{{{-{d_{\rm e}}} / 2}}}\exp \left[ {{{{{\left\| {{{{x}}_{{i}}} - {{{y}}_{{j}}}} \right\|}^2\Big/{\left(2\sigma ^2\right)}}}}{{}}} \right].$

用于表示第i个数据点由第j个高斯成员函数生成的概率. $\sigma^2 $为高斯成员的协方差， ${C_{ij}}$ 为各个高斯成员的先验概率，即 ${C_{ij}} = {m^{{{ - 1}}}}$ 且满足 $\displaystyle\sum\nolimits_{j = 1}^m {{C_{ij}}} = 1$. 在GMM模型中，先验概率人为设定. TPSMPM^[23]中指出，由于实际的先验参数未知，保守地将每个先验概率赋予相同的值是合理的选择.

当存在冗余点（Outlier）时，曲线拟合和特征提取都将变得特别困难，且冗余点会对转换模型造成严重的影响，从而产生不可靠的对应关系. 因此添加一个概率为1/n的均匀分布，用来处理噪声和冗余点， $\xi $ 为该均匀分布的权重系数，且满足 $0< \xi <1$.

由于 $p({{{x}}_{{i}}}|{{{y}}_{{j}}})$ 在进行对应关系评估时，仅考虑点集X到点集Y之间的欧氏距离特征，当距离特征相似时，会得到错误的对应关系. 例如：假设 ${{{x}}_{{ a}}}$ 与 ${{{x}}_{{ b}}}$ 具有完全不同的局部结构，若两者到 ${{{y}}_{{\rm cer}}}$ 的欧氏距离相等，则有 $p({{{x}}_{{ a}}}|{{{y}}_{{j}}}) = p({{{x}}_{{ b}}}|{{{y}}_{{j}}})$，这意味着 ${{{x}}_{{ a}}}$ 和 ${{{x}}_{{ b}}}$ 错误地得到了与 ${{{y}}_{{\rm cer}}}$ 对应的相同概率.

为了解决上述问题，提出能够融合上述多种图像特征的基于混合特征的高斯混合模型（MGMM）：

(7) $p({{{x}}_{{i}}}) = (1 - \xi )\sum\limits_{j = 1}^m \;\left[{{M_{ij}}} f({{{x}}_{{i}}}|{{{y}}_{{j}}})\right] + \xi/{n}.$

式（7）不仅易于建模且计算方便，能够合理地表示点集间的概率分布，具有较强的鲁棒性. 式（7）中各个参数的详细介绍如下. 首先，MGMM的密度函数

(8) $f({{{x}}_{{i}}}|{{{y}}_{{j}}}) = {(2\pi {\sigma ^2})^{{{-{{{d}}_{{\rm e}}}} / 2}}}\exp\; [ - ({G_{ij}} + \alpha {L_{ij}})].$

其中， ${G_{m \times n}}$ 和 ${L_{m \times n}}$ 分别用于表示全局和局部结构差异，其次， ${M_{ij}}$ 为先验概率（同上述的 ${C_{ij}}$），根据式（1）得到的SIFT距离，N（或M）个点就会有N（或M）维度的矩阵，对其每一行进行归一化（令和为1，然后再用1去减）后，就可以充当先验，此处用于表示图像的灰度信息差异. M_ij的先验概率被定义为

${M_{ij}} = 1 - {\tilde s_{ij}},$

在MGMM中用于表示各个高斯成员的比重. 再次， ${G_{ij}} = {{{g_{ij}}} / ({2{\sigma ^2}})}$ 为全局几何结构特征差异， ${\sigma ^2}$ 为MGMM的协方差. 最后， ${L_{ij}} = \alpha {l_{ij}}$ 为局部几何结构特征差异， $\alpha $ 为常量. 本算法采用贝叶斯法则评估点集之间的对应关系，得到二维高斯混合模型的后验概率如下：

(9) ${p_{ij}} = \frac{{{M_{ij}}f({{{x}}_{{i}}}|{{{y}}_{{j}}})}}{{\displaystyle\sum\nolimits_{k = 1}^m {\left( {{M_{ik}}f({{{x}}_{{i}}}|{{{y}}_{{j}}})} \right) + \xi /{n}} }}.$

通过式（10），评估出源点集与目标点集的对应概率矩阵 ${P_{m \times n}}$，用于实现后述配准过程中点集之间对应关系的评估.

分别将由式（1）、（2）和（4）所求得的SIFT距离特征 ${s_{ij}}$、欧氏距离特征 ${g_{ij}}$、局部距离特征 ${l_{ij}}$ 代入式（9）中，得到如下后验概率公式：

(10) ${p_{ij}}^* = \frac{{(1 - {{\tilde s}_{ij}}){ E}}}{{\displaystyle\sum\nolimits_{k = 1}^m {\left( {(1 - {{\tilde s}_{ij}}){ E}} \right) + {{(2\pi {\sigma ^2})}^{{{{d_{\rm e}}}}/{2}}}{\xi }/{n}} }}.$

其中， ${\tilde s_{ij}} \in [0, 1]$ 为归一化后的SIFT特征距离差异， ${{E}} = \exp\; \left[ { - \left( {\displaystyle\frac{{{{\left\| {{x_i} - {y_j}} \right\|}^2}}}{{2{\sigma ^2}}} + \displaystyle\sum\nolimits_{k = 1}^K {\left( {{{\left\| {T(N{{({x_i})}_k},{y_j}) - N{{({y_j})}_k}} \right\|}^2}} \right)} } \right)} \right],$ K为局部距离特征中最近相邻点的数量，使用大小为 $m \times n$ 的对应概率矩阵 ${{ P}^*}$ 左乘目标点集X，可以得到对应目标点集 ${{{ X}}^*}$，即

(11) ${{{ X}}^{{*}}}{{ = }}{{{{ P}}}^{{*}}}{{ X }}.$

从而实现点集之间对应关系的评估. 由式（10）可知，当3个特征中有某1个存在相似，就能根据其余特征的差异而得到可靠的对应关系评估.

（1）参数评估：采用 ${\rm {L_2}E}$ ^[24]估计子（一种统计学中用于最小化密度之间距离的鲁棒性估计子）进行参数优化.

（2）基于运动一致性的几何结构约束：图像发生非刚性形变，导致2个点集之间的映射关系过于自由. 因此，本研究引入基于运动一致性^[25]的几何结构约束以解决这种不适定问题（ill-posed problem）. 在运动感知中，有一些重要的现象涉及到一致性，其中，速度一致性是对运动进行平滑的一种特殊方式. 图3给出了在2种不同的速度场约束下分别得到的对应关系.

上述运动一致性理论写作如下形式：

(12) $R(T) = {\left\| T \right\|^2}.$

这意味着在SIFT特征点之间融入了运动一致性约束，并抑制了整个变换过程中速度场不一致的变换.

（3）变换建模：在再生核希尔伯特空间（reproducing kernel hilbert space，RKHS）^[26]中对非刚性转换T进行建模. 首先选用高斯核来定义RKHS，形式如下：

(13) $ \varPhi ({y_i}, {y_j}) = \exp\; \left( - \frac{1}{{2{\beta ^2}}}{\left\| {{y_i} - {y_j}} \right\|^2}\right).$

式中： $\beta $ 为常数，用于控制空间平滑度，矩阵 $ \varPhi $ 的大小为 $m \times m$. 根据里斯表示定理（Resiz representation theorem），得到最优变换函数T的形式如下：

(14) $ T(Y) = Y + { \varPhi} { W}. $

其中， ${{ W}_{m \times {d_{\rm e}}}} = {\{ {w_1}, {w_2},\cdots, {w_m}\} ^{\rm T}}$ 为需要求得的系数矩阵. 可以通过迭代求解基于几何结构约束的 ${\rm {L_2}E}$ 能量方程：

$Q(T, {\rho ^2}) = {\rm {L_2}E}(T, {\rho ^2}) + \frac{\lambda }{2}R(T).$

其中， ${\rho ^2}$ 为用来表示源点集和目标点集之间偏差的密度模型参数， ${\rm {L_2}E}( \cdot )$ 表示 ${\rm {L_2}E}$ 估计， $R( \cdot )$ 表示图像非刚性变换过程中基于Tikhonov正则化原理的几何结构约束，用于克服图像变换过程中的不适定问题， $\lambda $ 用于控制约束强度，从而得到该系数矩阵W. 上述的非刚性转换 $T(Y)$ 等同于在其原始位置加入位移函数，即

(15) $T(Y) = Y + V(Y).$

由于 $R(T)$ 在仿射变换下是不变的，则 $R(T)$ 等价于 $R(V)$，即R（T）=R（Y+V（Y））=R（V）. 本研究用V代替T并将基于几何结构约束的 ${L_2}E$ 能量方程重写为如下形式：

(16) $Q({ W}, {\rho ^2}) = \frac{1}{{{2^{{d_{\rm e}}}}{{({\text π}\rho )}^{{{{d_{\rm e}}}}/{2}}}}}{ B} + \lambda \ {\rm tr}\ ({ W \varPhi W}).$

其中，

${B}{{ = }}\exp \left[ {{{{{\left\| {{x_i}^{{*}} - \displaystyle\sum\nolimits_{j = 1}^m { \varPhi ({y_i}, {y_j}){w_j}} } \right\|}^2}}}\Bigg/{{\left(2{\rho ^2}\right)}}} \right] ,$

常数 $\lambda $ 用于控制约束项的强度， ${\rm tr}( \cdot )$ 为矩阵的迹. 系数矩阵W可由对式（16）求相应的偏导得到：

(17) $\frac{{\partial Q}}{{\partial { W}}} = \frac{{2{ \varPhi} [{{U}} \odot ({{E}} \otimes {{1}})]}}{{m{\rho ^2}{{(2\pi {\rho ^2})}^{{{{d_{\rm e}}}}/{2}}}}} + 2\lambda { W}.$

其中， ${{ U}_{m \times {d_{\rm e}}}} = { \varPhi} { W} - { X^*}$， ${{{E}}_{m \times 1}} = \exp \left\{ {{{{\rm diag}({{U}}{{{U}}^{\rm T}})} / (\rm {2{\rho ^2}}}} )\right\}$， ${{\bf{1}}_{1 \times {d_{\rm e}}}}$ 为所有元素均是1的行向量， ${\rm diag}\,{\rm tr}\;( \cdot )$ 为对角矩阵， $ \odot $ 和 $ \otimes $ 分别表示为矩阵的Hadamard乘积和Kronecker乘积. 由于式（17）并非线性系统，采用基于梯度的数值优化技术进行求解. 此外，本文使用关于协方差 ${\rho ^2}$ 的另一个确定性退火技术来改善数值求解的收敛性，即通过 ${\rho ^2} = k{\rho ^2}$ （k为常数且k $ \epsilon$（0，1））将一个数值较大的 ${\rho ^2}$ 逐渐减小，该方案将确保式（16）的最终收敛.

求得非刚性形变系数矩阵W之后，可以根据式（15）得到新的Y，然后通过 ${\sigma ^2} \leftarrow \varepsilon {\sigma ^2}$ 对MGMM的协方差进行退火，接着进入新一轮的对应关系评估，继续进行配准，直到达到最大迭代次数，配准结束.

基于源点集和经过变换的源点集构造映射函数 $F= \left\{ {{y_i}, {{\hat y}_i}} \right\}_{i = 1}^m$，从而对图像进行配准. 为避免输出图像发生孔和（或）重叠的现象（由离散化和四舍五入所引起的配准错误），本文采用反推法（Backward approach）^[27]来获取变换后的图像 ${I_{\rm{t}}}$，即通过使用参考图像 ${I_{\rm{r}}}$ 的网格化图像以及特征点集之间的映射关系来指导针对待配准图像 ${I_{\rm{s}}}$ 的像素重采样处理.

设源点集为 ${\hat Y}$，目标点集为Y_q，使用TPS函数建立图像非刚性变换模型，并将TPS核定义为

(18) $k\left( {{{{\hat y}}_{{i}}}, {{{\hat y}}_{{j}}}} \right) = {\left\| {{{{\hat y}}_{{i}}} - {{{\hat y}}_{{j}}}} \right\|^2}\log\; \left\| {{{{\hat y}}_{{i}}} - {{{\hat y}}_{{j}}}} \right\|.$

由此，TPS转换模型可以由以下公式得到：

(19) ${{\varDelta}} = {\left( {\begin{array}{*{20}{c}} {{k}}&{{Q}} \\ {{{{Q}}^ {\rm{T}}}}&{{O}} \end{array}} \right)^{ - 1}}{\left( {{{Y_{\rm q}\; 0\; 0 \;0}}} \right)^ {\rm{T}}}.$

其中，TPS模型 ${{\varDelta}}$ 为（m+3）×3的矩阵，O为3×3的零矩阵，Q为图像像素点的齐次坐标矩阵，其第i行写为 $\left( {1, {{{\hat y}}_{ia}}, {{{\hat y}}_{ib}}} \right)$，其中， ${{\hat y}_{ia}}$ 和 ${{\hat y}_{ib}}$ 为像素 ${{\hat y}_i}$ 在图像上的坐标. 对图像像素点进行二维索引化处理，可以得到网格 ${\varOmega ^{\rm{r}}}_{z \times 2} = {\left\{ {{\tau ^{\rm{r}}}_1, {\tau ^{\rm{r}}}_2,\cdots, {\tau ^{\rm{r}}}_z} \right\}^{\rm{T}}}$，其中， $z = {N'}_{\rm w} \times {N'}_{\rm h}$表示图像大小. 若令网格点集 ${\varOmega ^{\rm{r}}}$ 为源点集， ${{\varDelta}}$ 为TPS转换模型，则转换后的网格图像可由以下公式得到：

(20) ${\hat \varOmega ^{\rm{r}}}_{z \times 3} = \left( {{k'}{Q'}} \right)\varDelta.$

通过 ${\hat \varOmega ^{\rm{r}}} \leftarrow \left( {{{\hat \varOmega }^{\rm{r}}}_{\left( {., 1} \right)} {{\hat \varOmega }^{\rm{r}}}_{\left( {., 2} \right)}} \right)$ 将网格齐次坐标恢复为二维笛卡尔坐标，式（20）中， ${k'} = $ ${\left\| {{\tau ^{\rm{r}}}_i - {{\hat y}_j}} \right\|^2}\log\; $ $\left\| {{\tau ^{\rm{r}}}_i - {{\hat y}_j}} \right\| $ 是大小为z×m的TPS核，大小为z×3的矩阵 ${{{{ Q}}}^{{'}}}$ 是网格图像 ${\varOmega ^{\rm{r}}}$ 的齐次坐标，第i行表示为 $\left( {1, {{\left( {{\tau ^{\rm{r}}}_i} \right)}_a}, {{\left( {{\tau ^{\rm{r}}}_i} \right)}_b}, } \right)$.

设 ${\varOmega ^{\rm{s}}}$ 为获取自参考图像的网格图像，遵循下式：

(21) ${\hat \varOmega ^{*}} = {\hat \varOmega ^{\rm{r}}} \cap {\varOmega ^{\rm{s}}}.$

得到针对待配准图像 ${I_{\rm{s}}}$ 进行像素重采样的坐标范围 ${\hat \varOmega ^*}$，据此得到变换后的图像 ${I_{\rm{t}}}$，并将非采样区域的像素填充为黑色. 在重采样过程中，本文使用双三次插值（Bicubic interpolation）算法来提高变换后图像 ${I_{\rm{t}}}$ 的平滑度，更确切地说，通过双三次插值所创建的图像 ${I_{\rm{t}}}$ 上的每一个像素均是由该像素及其附近（4×4）窗口内的邻近像素值确定.

（1）单一特征与混合特征配准精度比较.

由于航拍视角的变化造成城市绿地遥感图像间产生非刚性畸变、重叠性低等问题，增加了图像配准的难度. 有效地选择图像特征点集，不但会减少图像配准的计算复杂度，还会提升图像配准精度. 使用SIFT算法提取城市绿地遥感图像的特征点. 为了证明使用混合特征的MGMM的优点，比较使用单一特征和混合特征的MGMM点集配准过程，对配准过程进行可视化展示，如图4所示。其中，t为迭代次数，点集X和Y提取自本文测试数据所用的遥感图像对，图中的横纵坐标表示特征点的空间坐标. 在第一次迭代中，可以看到单一特征配准最初评估所得的对应点集位于区域的质心，然而，使用MGMM进行配准，其分布更接近真正的目标点集. 随着迭代次数增加，使用单一特征配准得到的假想目标点集仍分布于区域质心周围，相比之下，MGMM迭代更新得到的源点集位置更优. 使用混合特征配准相对于单一特征配准不仅配准效果更佳且迭代次数减少. 由图4中矩形框可以看出，基于单一特征的点集配准存在许多圆圈（假想目标点集）缺失，而在所提出的基于图像混合特征的配准中，几乎所有的圆圈都正确地分布到目标位置，冗余点位置基本没有任何圆圈，说明冗余点得到了有效排除.

（2）不同混合特征配准精度比较.

由上述（1）中的介绍可知，城市绿地遥感图像间具有非刚性畸变、重叠性低等问题，增加了配准的难度. 所提方法是将SIFT距离特征、欧氏距离特征、局部距离特征进行融合，形成特征互补，而在实验部分，对比方法中GLMDTPS也使用了混合特征，为了进一步验证本文方法的精度，比较使用不同混合特征（如本文算法和GLMDTPS）的点集配准过程，并对配准过程进行了可视化展示，如图5所示。虽然GLMDTPS也使用了混合特征，但是结果存在着明显的偏差（见图5中方框），而本文的方法却表现出良好的性能. GLMDTPS利用基于欧氏距离的全局特征差异描述子和局部距离特征差异描述子建立一一对应关系，可以用于具有不同程度形变、旋转图像的配准. 从图5中可以看出，GLMDTPS算法在迭代初期分布比较均匀，但由于一一对应关系，GLMDTPS无法处理冗余点，而城市绿地遥感图像间具有非刚性畸变、重叠度低等问题，虽然其在迭代初期点集分布均匀，但最终配准效果却不理想. 本文的方法有效地解决了这一问题，表现出良好的性能.

（3）混合特征与基于图像灰度信息方法配准精度比较.

实验部分的对比方法SIFT算法是基于图像灰度信息进行匹配的（SIFT算法不具有迭代过程，无法展示其点集迭代配准过程），由于城市绿地遥感图像很多局部区域灰度信息近似，从而造成误配，而本文的方法使用混合特征，形成特征互补，配准效果更好. 两者的配准性能可以参考实验部分的所提算法和SIFT算法的配准结果.

（4） $\rm {L_2}E$ 的鲁棒性.

$\rm {L_2}E$ 估计子具有较强的鲁棒性，其与最大似然估计（maximum likelihood estimation, MLE）相比，更适用于分析包含冗余点的大规模数据集，而 $\rm MLE$ 在提取的特征点包含大量冗余点时会产生严重偏差. 假设参数模型 $p\left( {z|u} \right)$，对 $\rm {L_2}$ 距离进行最小化，即

　　 ${u^*} = \arg \mathop {\min }\limits_u \int [ p(z|u) -$ $ p(z|{u_0})]^2 {\rm d}z$，

其中，参数u的真实值 ${u_0}$未知，省略与u无关的常数项 $\int {p{{(z|{u_0})}^2}} {\rm d}z$. 基于 $\rm {L_2}$ 距离，参数u可以由 $\rm {L_2}E$ 估计子对模型 $p\left( {z|u} \right)$ 最小化来估计得到：

${u}^{*{\rm {L_2}E}} = \arg \mathop {\min }\limits_u$ $ \left[ {\int {p{{(z|{u_0})}^2}} {\rm d}z - \displaystyle\frac{2}{m}\sum\nolimits_{i = 1}^m {p(z|{u_0})} } \right]$

通过估算在含冗余点的正态分布N(5,1)干扰下，正态分布N(0,1)所产生的均值，对比L₂E和最大似然估计MLE的鲁棒性，如图6所示（图中垂直的虚线表示极值，纵坐标表示方程的值）。例如，一个大小为600的数据样本表示内点，其服从正态分布N(0,1)，而5组数量递增的冗余点n=30,60,120,180,300由正态分布N(5,1)生成. ${\rm {L_2}E}$ 在“0”附近有一个全局最小值，在“5”附近存在一个局部最小值，均符合内点H和冗余点B分布，并且随着冗余点的增加极值会变得更深. 与此相反，MLE随着冗余点的增多偏差增大，稳定性变差.

使用默认阈值1.5提取假定的对应点集和SIFT描述子. 常数 $\alpha $ 用于控制全局和局部几何结构特征差异之间的平衡. 常数 $\xi $ 作为冗余点均匀分布 ${1}/{n}$ 的权重. 本研究使用的2个确定性退货方案分别为： ${\sigma ^2} \leftarrow \varepsilon {\sigma ^2}$ 和 ${\rho ^2} = k{\rho ^2}$，在对非凸性（non-convexity）问题处理上作用明显. 在退火中，若不考虑效率问题，可以优先选择每次退火都乘以一个小于1的数，使得退火速度变慢. 常数 ${t\,^\sigma }$ 和 ${t\,^\rho }$ 用于确定本文方法和Matlab函数fminunc的最大迭代次数. $\beta $ 用于确定源点集邻域的大小. $\lambda $ 用于控制变换模型T中约束项的强度. 分别设置 $\alpha = 10$， $\xi = 0.3$， $\varepsilon = 0.9$， $k = 0.75$， $\beta = 2$ 和 $\lambda = 3$，并将 ${\sigma ^2}$ 和 ${\rho ^2}$ 分别初始化为1和0.05. 由于迭代需要终止条件，使用Matlab中函数fminunc的如下选项：{MaxIter=50}，即设置 ${t\,^\rho } = 50$， ${t\,^\sigma } = 100$. 在确定参数的过程中，使用实验数据集I和II（参考实验部分2.1节）进行实验. 在保证其他参数不变的情况下，改变某一参数值进行实验，从而求得参数不同情况下图像配准的均方根误差（root of mean square error，RMSE）和平均绝对误差（mean absolute error，MAE），此处将其作为判断的主要标准，根据RMSE和MAE的值确定最佳参数值. 其中2个参数（ $\beta $ 和 $\varepsilon $）的实验，如图7所示，从图中可以看出 $\beta = 2$ 和 $\varepsilon = 0.9$ 时，误差δ最小，性能最佳.

本实验所覆盖区域主要为昆明市. 研究区域内土地覆盖类型包括建筑物、道路、空地、水体和城市绿地. 为了保证实验的广泛性和鲁棒性，使用不同角度、不同时间的城市绿地遥感图像，根据遥感图像的类型，实验数据集主要划分成两大类：1）选取60组不同时相且存在视角变化的卫星遥感图；2）选取60组不同视角的小型无人机低空遥感图像. 所有卫星遥感图像来自谷歌地球，小型无人机低空遥感图像来自720云平台和大疆精灵4无人机航拍. 图片分辨率大小范围为620×480~1 300×890，实验中所用的数据集如表1所示. 实验中仿真平台硬件环境为：Intel（R）Core（TM）i5-3230M CPU，2.6 GHz，8 G内存的PC机；软件开发工具为 Visual Studio 2013与Matlab 2014b.

主要采用均方根误差和平均绝对误差评估图像配准精度. 首先，在待配准图像和参考图像上手动确定至少15对地面真值（Landmarks）. 为了减小在评估过程中所产生的误差，所选择的Landmarks都是均匀分布于容易识别的区域. 评估中的相关定义如下：

(22) ${\rm RMSE} = \sqrt {\frac{1}{{{m'}}}\sum\nolimits_{i = 1}^{{m'}} {({{{y}}_i}^t - {{{x}}_i}^t)} }.$

RMSE反映出选定地标与实际位置间的均方根误差，其中 ${m'}$ 为所选地标的总数， ${{{y}}_i}^t$ 为对应于 ${{{x}}_i}^t$ 的地标.

(23) ${\rm MAE} =\frac{1}{{{{m'}}}} {{\displaystyle\sum\nolimits_{i = 1}^{{}} {\left| {{{{y}}_i}^t - {{{x}}_i}^t} \right|} }}.$

MAE通过计算选定地标与实际位置间的平均绝对差异，来评估算法的配准精度.

本系列实验针对数据集I，主要用于实现对城市绿地变化的长时间监测. 选用RMSE和MAE对不同方法进行定量评估比较，评估结果如表2所示(“−”表示由于所得转换图像过于模糊而无法得到对应的地面真值,导致无法计算所求评估值RMSE和MAE；粗体字表示最好的结果)图8展示了来自数据集I的3组具有代表性的图像配准示例，并与其他4种流行的算法（SIFT、CPD、RSOC、GLMDTPS）进行了对比，结果表明本文算法在配准中均给出了更高的配准精度。图8中第一行展示了每对图像中的待配准图像和参考图像。每组配准方法用一行进行展示，左侧为转换图像，右侧为5×5的棋盘格。使用5×5的棋盘格将转换图像I_t与参考图像I_r交替展示，可视化最终配准结果。配准错误区域用矩形框标出。

在本组实验中，主要使用数据II中的小型无人机低空遥感图像，能够帮助实现对城市绿地信息的实时监测. 使用RMSE和MAE对本系列实验进行定量评估，其评估结果如表3所示(“−”表示由于所得转换图像过于模糊而无法得到对应的地面真值，导致无法计算所求评估值RMSE和MAE；粗体字表示最好的结果). “如图9所示为来自数据集II的3组具有代表性的小型无人机低空遥感图像配准结果，并与其他4种流行的算法（SIFT、CPD、RSOC、GLMDTPS）进行了对比，结果表明本文算法在配准中均给出了更高的配准精度。图9中第一行展示了每对图像中的待配准图像和参考图像。每组配准方法用一行进行展示，左侧为转换图像，右侧为5×5的棋盘格。使用5×5的棋盘格将转换图像I_t与参考图像I_r交替展示，可视化最终配准结果。配准错误区域用矩形框标出。

为了更加直观地了解本文方法的精度及鲁棒性，分别展示数据集I多时相、多视角的卫星遥感图像（见图10）和数据集II小型无人机低空遥感图像（见图11）的配准结果；使用RMSE和MAE对本文算法的精度进行定量评估，如表4和5所示. 实验结果表明：本文方法在所有实验中均保持较高的配准精度.

在城市绿地遥感图像的配准中，由于SIFT的优良特性，其被用于提取特征点，但SIFT算法仅考虑灰度信息进行匹配，而城市绿地通常局部布局相似，很多局部区域的灰度信息也相似，从而造成误配，因此使用SIFT进行城市绿地遥感配准时效果不理想，如图8中城市绿地卫星遥感图像和图9中城市绿地小型无人机低空遥感图像中第3行SIFT算法配准结果所示. CPD使用单一特征差异描述子、欧氏距离并利用高斯混合模型（GMM）建立对应关系，采用基于运动一致性的几何结构约束来调整点集之间的位移场，可以较好地处理噪声和缺失点等情形，但是，CPD也会受到基于单一特征的对应性评估限制. GLMDTPS利用基于欧氏距离的全局特征差异描述子和局部距离特征差异描述子的混合特征描述子建立对应关系，可以用于具有不同程度形变、旋转图像的配准，但是GLMDTPS对冗余点比较敏感. 与SIFT相比，CPD和GLMDTPS均在一定程度上提高了图像配准的精度，如图8中城市绿地卫星遥感图像和图9中城市绿地小型无人机低空遥感图像中第4行CPD、第6行GLMDTPS算法配准结果所示. RSOC采用图匹配技术，具有较高的鲁棒性，对错误匹配具有较强校正能力，其精度相对较高，如图8中城市绿地卫星遥感图像第二组（昆明洛龙公园）第5行RSOC算法配准结果所示，本组的2幅图像是多时相、多视角的卫星遥感图像，图像间存在非刚性畸变，容易出现误匹配，而RSOC可以去除误匹配，因此使用RSOC进行配准得到的结果相对较好，但是RSOC算法运行速度较慢.

本文算法以高斯混合模型（MGMM）为基础，并将全局几何结构特征、局部几何结构特征与灰度信息特征（SIFT特征）融入到MGMM中，形成特征互补，从而得到可靠的对应关系，因此可以较好地处理灰度信息相似、特征点集的几何结构特征相似、发生非刚性畸变、重叠度低的遥感图像，如图8中的第二组图像具有相似的灰度信息和特征点集几何结构，图9中第一组图像发生了非刚性畸变，重叠度低. 其他算法均未取得好的结果，而本文算法得到了非常精确的结果.

以上5种方法使用的主要特征包括：1）灰度信息特征；2）全局几何结构特征差异；3）局部几何结构特征差异；4）几何约束. SIFT只采用特征1）提取对应的特征点完成配准，而CPD和GLMDTPS还分别采用了特征2）、4）和2）、3），能够在一定程度上提高图像配准的精度. RSOC采用了鲁棒的图匹配（graph matching）技术去除误匹配，其性能相对较高. 本文算法配准结果优于其他配准结果，主要原因是采用了特征1）、2）、3）、4）特征，形成特征互补.

本文提出了一种基于图像混合特征的城市绿地遥感图像配准方法. 该研究的主要贡献如下：1）构建了基于混合特征的高斯混合模型（MGMM），并将全局几何结构特征、局部几何结构特征与灰度信息特征（SIFT特征）融入到MGMM中，形成特征互补，从而得到可靠的对应关系；2）为了防止遥感图像配准过程中出现由不适定问题（ill-posed problem）造成的错误匹配，在 ${\rm {L_2}E}$ 的能量方程中引入了几何结构约束项，从而实现准确的空间变化更新.

本研究分别使用卫星遥感图像、小型无人机低空遥感图像与4种同类算法（SIFT、CPD、RSOC、GLMDTPS）进行了对比实验. 结果证明：提出的算法可以有效解决城市绿地遥感图像配准时存在的灰度信息特征相似、特征点集的几何结构特征相似、图像间的非刚性畸变、重叠度低等问题，在城市绿地卫星遥感图像和无人机遥感图像配准中均表现出了较高的精度和鲁棒性.