CN103761739A - 一种基于半全局能量优化的影像配准方法 - Google Patents

Abstract

一种基于半全局能量优化的影像配准方法，包括根据地面采样距离对待配准影像进行预处理，使得参考影像和预处理后的待配准影像的分辨率相同；给定同名点搜索范围，针对参考影像相对于预处理后的待配准影像的二维视差场构建能量函数，用多个方向优化路径上的匹配代价来近似能量函数的匹配代价；进行匹配代价的聚合、二维视差场的计算及后处理和影像重采样，得到最终的配准影像。本发明技术方案实现过程简单，支持自动化运行，能够在没有任何几何变形的先验信息的情形下，有效地对待配准影像进行纠正，避免了配准模型选择的盲目性，并且能够处理较为严重的影像局部变形。

Description

一种基于半全局能量优化的影像配准方法

技术领域

本发明涉及遥感影像技术领域，尤其是涉及一种基于能量优化的影像配准方法。

背景技术

影像配准是将不同时间或不同传感器所获得的同一地区的影像归化到统一的坐标系中的操作，它是许多图像处理操作的先决条件，如影像融合、变化检测、统计模式识别、三维重构、地图修正等（孙家抦，2003）。传统的影像配准方法主要分为如下4个步骤：

1)影像特征提取，例如Harris角点特征（Harris，1988）、SUSAN角点特征（Smith，1997）、尺度不变特征（Scale-Invariant Features）（Lowe，2004）、最大稳定极值区域（MaximallyStable Extremal Regions）（Matas，2004）、加速稳健特征（Speeded Up Robust Features）（Bay，2006）等；

2)特征匹配及粗差剔除；

3)计算影像变换模型，如仿射变换、多项式、不规则三角网、薄板样条等；

4)根据变换模型对影像进行重采样。

然而，特征匹配本身会存在一定的误差。简单的影像变换模型，如仿射变换、低阶多项式等，可以通过变换参数的稳健估计使配准精度少受匹配误差的影响，但难以处理较复杂的影像变形。较为复杂的影像变换模型，如不规则三角网、薄板样条等，可以处理较复杂的影像变形，但必须有足够多的同名特征，且局部配准精度受匹配误差的直接影响。另外，影像变换模型的选择通常根据经验而定，带有一定的盲目性。

近十年来，密集匹配有了突破性的进展，涌现了一大批优秀的算法。Scharstein和Szeliski（2002）对密集匹配方法进行了分类和评估，他们将这些方法分为局部方法和全局方法两类。局部方法只利用对应点周围的局部影像信息，效率较高，但是缺少整体的约束，难以处理弱纹理、重复纹理、遮挡等问题，往往无法得到高精度的结果。全局方法通过对视差图D构建一个带平滑项的能量函数，将立体匹配问题转化为一个能量最小化问题。全局方法虽然能取得较好的匹配结果，但是往往需要进行多次迭代，计算耗时大。

最近，SGM算法以其出色的表现和较高的效率，在立体匹配、三维重建、移动机器人、驾驶员辅助系统中获得了巨大的成功（Hirschmüller，2011）。SGM算法采用多个方向的平滑约束来近似全局的平滑约束，计算量小，精度高，并且能够有效地处理遮挡和弱纹理的区域。另外，SGM算法还能够通过单指令多数据流（Single Instruction Multiple Data，SIMD）、图形处理单元（Graphic Processing Unit，GPU）和现场可编程门阵列（Field-Programmable GateArray，FPGA）等进行软硬件的加速，以满足一些实时性的需求。但尚未出现将SGM算法运用于影像配准的相应技术方案。另外，原始SGM算法是在核线影像上进行的，即影像不存在上下视差，只考虑地表高程起伏引起的水平视差，所以它的搜索范围是一维的；而对于影像配准来说，搜索范围往往是二维的，所以无法将SGM算法直接运用于影像配准。

发明内容

针对现有技术的不足，本发明设计了一种基于半全局能量优化的影像配准方法。

本发明技术方案提供一种基于半全局能量优化的影像配准方法，包括以下步骤：

步骤1，根据地面采样距离对待配准影像进行预处理，使得参考影像和预处理后的待配准影像的分辨率相同；

步骤2，将参考影像作为左影像，将预处理后的待配准影像作为右影像，给定同名点搜索范围，包括x方向视差的取值范围[dx_min,dx_max]，y方向视差的取值范围[dy_min,dy_max]；

步骤3，对参考影像相对于预处理后的待配准影像的二维视差场(Dx,Dy)构建能量函数如下，

$E(\mathrm{Dx},\mathrm{Dy})=\underset{p}{\mathrm{\Σ}}C(p,{\mathrm{Dx}}_p,{\mathrm{Dy}}_p)+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{1}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||=1]$

$+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{2}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||>1]$

式中，p为图像中的任一像点，Dx_p为像点p在x方向的视差，Dy_p为像点p在y方向的视差；C(p,Dx_p,Dy_p)为像点p在x方向视差为Dx_p，y方向视差为Dy_p时的匹配代价；||·||表示向量的范数；N_p为像点p的邻域；q为像点p邻域内的任一像点，Dx_q为像点q在x方向的视差，Dy_q为像点q在y方向的视差；T[·]在方括号中的表达式为真时取1，为假时取0；P₁为对视差变化量等于1时的惩罚系数；P₂为对视差变化量大于1时的惩罚系数，并且满足P₂＞P₁＞0；

步骤4，用多个方向优化路径上的匹配代价来近似步骤3所得能量函数的匹配代价，每条优化路径上按下式计算匹配代价，

L_r(p,dx,dy)＝C(p,dx,dy)+min(L_r(p-r,dx,dy),

L_r(p-r,dx-1,dy)+P₁,L_r(p-r,dx+1,dy)+P₁,

L_r(p-r,dx,dy-1)+P₁,L_r(p-r,dx,dy+1)+P₁,

$\underset{i,j}{\min }{L}_r(p-r,i,j)+P^2)-\underset{i,j}{\min }{L}_r(p-r,i,j)$

式中，r表示传播方向，p表示方向r优化路径上的任一像点；p-r表示像点p在方向r上的前一个像点；dx为x方向视差，取值范围为[dx_min,dx_max]；dy为y方向视差，取值范围为[dy_min,dy_max]；L_r(p,dx,dy)表示像点p在方向r优化路径上当x方向视差为dx，y方向视差为dy时的路径匹配代价；C(p,dx,dy)为像点p在x方向视差为dx，y方向视差为dy时的匹配代价；

为像点p-r在x方向视差范围[dx_min,dx_max]、y方向视差范围[dy_min,dy_max]内最小的路径匹配代价，这里的i,j也表示x方向视差和y方向的视差；

步骤5，匹配代价的聚合，包括按下式将各个方向优化路径上的匹配代价相加，得到匹配代价聚合结果S(p,dx,dy)，

$S(p,\mathrm{dx},\mathrm{dy})=\underset{r}{\mathrm{\Σ}}{L}_r(p,\mathrm{dx},\mathrm{dy})$

步骤6，二维视差场的计算，包括对每个像素p，在使得匹配代价聚合结果S(p,dx,dy)最小的二维视差值(dx_p,dy_p)附近进行局部曲面拟合，求得拟合曲面的极小值作为子像素的二维视差，得到整个二维视差场；

步骤7，对步骤6所得二维视差场进行后处理；

步骤8，影像重采样，包括利用步骤7所得二维视差场，通过数字微分纠正的反解法对预处理后的待配准影像进行纠正，得到最终的配准影像。

本发明依据二维视差场中的局部平滑约束，对二维视差场构建能量函数，并采用半全局优化方法求解，利用解得的二维视差对影像进行微分纠正，实现待配准影像和参考影像的精确配准。通过对原始算法的能量函数进行相应的改进，将其推广至二维情形，将影像配准问题变成一个能量最小化问题，并借用SGM密集匹配方法来解决，可以有效地避免配准模型选择的盲目性。并且，本发明技术方案实现过程简单，支持自动化运行，能够在没有任何几何变形的先验信息的情形下，有效地对待配准影像进行纠正，避免了配准模型选择的盲目性，并且能够处理较为严重的影像局部变形。

附图说明

图1为本发明实施例的多个方向的动态规划示意图。

具体实施方式

本技术方案实现过程简单，支持采用计算机软件方式实现自动化运行。以下结合附图和实施例详细说明本发明技术方案。

实施例包括以下步骤：

步骤1，根据地面采样距离（Ground Sampling Distance，GSD）对待配准影像进行预处理，使得参考影像和预处理后的待配准影像的分辨率相同。一般参考影像的分辨率较高，待配准影像的分辨率较低，在以下步骤中，实施例以这种情形为例进行说明，预处理时对待配准影像上采样，将参考影像作为算法输入的左影像，将上采样后的待配准影像作为右影像。若参考影像的分辨率较低，待配准影像的分辨率较高，则预处理时对待配准影像下采样，将参考影像作为算法输入的左影像，将下采样后的待配准影像作为右影像，后续处理过程相同。

步骤2，将参考影像作为算法输入的左影像，将上采样后的待配准影像作为右影像，给定同名点搜索范围，记x方向搜索范围（即x方向视差的取值范围）为[dx_min,dx_max]，y方向搜索范围（即y方向视差的取值范围）为[dy_min,dy_max]。

步骤3，对参考影像相对于上采样后的待配准影像的二维视差场(Dx,Dy)构建能量函数E(Dx,Dy)：

$E(\mathrm{Dx},\mathrm{Dy})=\underset{p}{\mathrm{\Σ}}C(p,{\mathrm{Dx}}_p,{\mathrm{Dy}}_p)+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{1}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||=1]$

$+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{2}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||>1]$

式中，p为图像中的任一像点，Dx_p为像点p在x方向的视差，Dy_p为像点p在y方向的视差；C(p,Dx_p,Dy_p)为像点p在x方向视差为Dx_p，y方向视差为Dy_p时的匹配代价；||·||表示向量的范数；N_p为像点p的邻域；q为像点p邻域内的任一像点，Dx_q为像点q在x方向的视差，Dy_q为像点q在y方向的视差；T[·]在方括号中的表达式为真时取1，为假时取0；P₁为对视差变化量等于1时的惩罚系数；P₂为对视差变化量大于1时的惩罚系数，并且满足P₂＞P₁＞0。

本步骤考虑基于全局平滑约束的视差场能量最小化模型，根据视差的局部平滑性假设，对视差场D构建一个带平滑项的能量函数E(D)，将密集同名点的匹配问题转化为一个能量最小化问题：

E(D)＝E_data(D)+λE_smooth(D)

式中，E_data(D)为数据项，表示左右影像在视差场为D时所有像点匹配代价的总和；E_smooth(D)为平滑项，是基于局部视差变化的平滑性假设而引入的项，它对局部视差不连续进行了惩罚；λ为惩罚因子，表示对视差不连续的容忍程度，惩罚因子越大，视差的局部平滑性就越明显。

在本发明方法中，即对图像中的每个像素p，将其在x方向视差为Dx_p，y方向视差为Dy_p时的匹配代价进行累加得到的数据项。

$E_{\mathrm{smooth}}\left(D\right)=\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{1}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||=1]+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{2}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||>1],$

即对当前像素p与其领域像素q存在视差变化时的惩罚项。

但是全局模型被证明是一个NP-hard问题，无法在多项式时间内精确求解出来，所以需要对模型进行一定的近似。

步骤4，匹配代价的计算。采用半全局的思路，即用多个方向优化路径上的带平滑约束的匹配代价来近似步骤3能量函数中带全局平滑约束的匹配代价。与原始SGM算法的优化路径相同，采用8邻域或16领域的优化路径，如图1所示16领域的优化路径，从16个方向传播到像点P。每个像点的匹配代价只考虑从这几个方向上的所有像素迭代传播过来的平滑约束。

实施例中每条优化路径上按下式计算匹配代价：

L_r(p,dx,dy)＝C(p,dx,dy)+min(L_r(p-r,dx,dy),

L_r(p-r,dx-1,dy)+P₁,L_r(p-r,dx+1,dy)+P₁,

L_r(p-r,dx,dy-1)+P₁,L_r(p-r,dx,dy+1)+P₁,

$\underset{i,j}{\min }{L}_r(p-r,i,j)+P_2)-\underset{i,j}{\min }{L}_r(p-r,i,j)$

式中，r表示传播方向；p表示方向r优化路径上的任一像点；p-r表示像点p在方向r上的前一个像点；dx为x方向视差，取值范围为[dx_min,dx_max]；dy为y方向视差，取值范围为[dy_min,dy_max]；L_r(p,dx,dy)表示像点p在方向r优化路径上当x方向视差为dx，y方向视差为dy时的路径匹配代价；C(p,dx,dy)为像点p在x方向视差为dx，y方向视差为dy时的匹配代价；

为像点p-r在x方向视差范围[dx_min,dx_max]、y方向视差范围[dy_min,dy_max]内最小的路径匹配代价，这里的i,j也表示x方向视差和y方向的视差，为了避免与式中特定的视差值dx,dy混淆，用不同的字母表示以示区分。

步骤5，匹配代价的聚合，即将各个方向优化路径上的匹配代价相加，得到匹配代价聚合结果S(p,dx,dy)：

$S(p,\mathrm{dx},\mathrm{dy})=\underset{r}{\mathrm{\Σ}}{L}_r(p,\mathrm{dx},\mathrm{dy})$

最终，视差场D的半全局的能量函数就变成了：

$E\left(D\right)=\underset{p}{\mathrm{\Σ}}S(p,\mathrm{dx},\mathrm{dy})=\underset{p}{\mathrm{\Σ}}\underset{r}{\mathrm{\Σ}}{L}_r(p,\mathrm{dx},\mathrm{dy})$

这就是基于半全局平滑约束的视差场能量最小化模型，并且该模型可以通过半全局优化算法有效地进行求解。

步骤6，二维视差场的计算。对每个像素p，在使得S(p,dx,dy)最小的二维视差值(dx_p,dy_p)附近进行局部曲面拟合（张祖勋，1997），求得拟合曲面的极小值，即子像素的二维视差。最终得到整个二维视差场。

步骤7，对步骤6所得二维视差场的后处理。后处理的目的是去除一些明显不一致的离群值，现有技术中常见的后处理有3*3窗口的中值滤波、小斑块滤除、左右一致性检查、局部平面拟合、空洞填补等等，具体实施时可根据具体情况选用。

步骤8，影像重采样。利用步骤7所得二维视差场，通过数字微分纠正的反解法（张祖勋，1997）对上采样后的待配准影像进行纠正，得到最终的配准影像。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (1)

1.一种基于半全局能量优化的影像配准方法，其特征在于，包括以下步骤：

步骤1，根据地面采样距离对待配准影像进行预处理，使得参考影像和预处理后的待配准影像的分辨率相同；

步骤2，将参考影像作为左影像，将预处理后的待配准影像作为右影像，给定同名点搜索范围，包括x方向视差的取值范围[dx_min,dx_max]，y方向视差的取值范围[dy_min,dy_max]；

步骤3，对参考影像相对于预处理后的待配准影像的二维视差场(Dx,Dy)构建能量函数如下，

$E(\mathrm{Dx},\mathrm{Dy})=\underset{p}{\mathrm{\Σ}}C(p,{\mathrm{Dx}}_p,{\mathrm{Dy}}_p)+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{1}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||=1]$

$+\underset{q\∈N_p}{\mathrm{\Σ}}{P}_{2}T\left[\right||{({\mathrm{Dx}}_p,{\mathrm{Dy}}_p)}^T-{({\mathrm{Dx}}_q,{\mathrm{Dy}}_q)}^T||>1]$

式中，p为图像中的任一像点，Dx_p为像点p在x方向的视差，Dy_p为像点p在y方向的视差；C(p,Dx_p,Dy_p)为像点p在x方向视差为Dx_p，y方向视差为Dy_p时的匹配代价；||·||表示向量的范数；N_p为像点p的邻域；q为像点p邻域内的任一像点，Dx_q为像点q在x方向的视差，Dy_q为像点q在y方向的视差；T[·]在方括号中的表达式为真时取1，为假时取0；P₁为对视差变化量等于1时的惩罚系数；P₂为对视差变化量大于1时的惩罚系数，并且满足P₂＞P₁＞0；

步骤4，用多个方向优化路径上的匹配代价来近似步骤3所得能量函数的匹配代价，每条优化路径上按下式计算匹配代价，

L_r(p,dx,dy)＝C(p,dx,dy)+min(L_r(p-r,dx,dy),

L_r(p-r,dx-1,dy)+P₁,L_r(p-r,dx+1,dy)+P₁,

L_r(p-r,dx,dy-1)+P₁,L_r(p-r,dx,dy+1)+P₁,

$\underset{i,j}{\min }{L}_r(p-r,i,j)+P^2)-\underset{i,j}{\min }{L}_r(p-r,i,j)$

式中，r表示传播方向，p表示方向r优化路径上的任一像点；p-r表示像点p在方向r上的前一个像点；dx为x方向视差，取值范围为[dx_min,dx_max]；dy为y方向视差，取值范围为[dy_min,dy_max]；L_r(p,dx,dy)表示像点p在方向r优化路径上当x方向视差为dx，y方向视差为dy时的路径匹配代价；C(p,dx,dy)为像点p在x方向视差为dx，y方向视差为dy时的匹配代价；为像点p-r在x方向视差范围[dx_min,dx_max]、y方向视差范围[dy_min,dy_max]内最小的路径匹配代价，i,j也表示x方向视差和y方向的视差；

步骤5，匹配代价的聚合，包括按下式将各个方向优化路径上的匹配代价相加，得到匹配代价聚合结果S(p,dx,dy)，

$S(p,\mathrm{dx},\mathrm{dy})=\underset{r}{\mathrm{\Σ}}{L}_r(p,\mathrm{dx},\mathrm{dy})$

步骤6，二维视差场的计算，包括对每个像素p，在使得匹配代价聚合结果S(p,dx,dy)最小的二维视差值(dx_p,dy_p)附近进行局部曲面拟合，求得拟合曲面的极小值作为子像素的二维视差，得到整个二维视差场；

步骤7，对步骤6所得二维视差场进行后处理；

步骤8，影像重采样，包括利用步骤7所得二维视差场，通过数字微分纠正的反解法对预处理后的待配准影像进行纠正，得到最终的配准影像。

Images