CN102663680A - 基于面特征的图像几何校正方法 - Google Patents

Abstract

基于面特征的图像几何校正方法，提出了点到多边形、多边形到多边形的距离、距离矢量的定义，并给出了求解点到多边形距离矢量、多边形到多边形距离矢量的算法，在此基础上将面状地物作为控制特征，建立了一种新的图像几何校正方法。该方法实用性很强，与影像的具体成像模型无关，适用于各种遥感影像、医学影像以及普通照片。

Description

基于面特征的图像几何校正方法

技术领域

本发明涉及一种图像几何校正方法，能够应用在遥感、测绘、图像处理、机器人视觉等领域。

技术背景

遥感图像几何校正是遥感图像进一步应用的基础。遥感影像在成像过程中受到各种复杂因素的影响，使其产生几何变形。图像几何校正的实质就是用数学模型来描述物方空间坐标系中的地面点坐标与它在图像平面上像点坐标之间的几何关系，从而消除几何畸变，生成具有地理参考的影像。因此，遥感图像高精度几何定位和几何纠正的关键是建立精确的传感器成像的数学模型。

传统的几何校正模型是基于控制点建立的。在实际应用中，当能够准确获取地面控制点时，采用地面控制点进行几何校正能够较好地改善图像的定位精度；然而在很多困难地区的卫星影像上，如沙漠、山区等，很难辨认出准确的物理点特征，因此难以获得足够数量和足够准确的地面控制点，而线特征、面特征则更易获得(如已知的道路、水体等)。在以不同分辨率影像作为参考数据源时，点的位置坐标也难以准确确定。事实上，就特征提取而言，提取有意义的面特征比提取有意义的点特征更容易，而且，无论自然环境还是人工环境中都存在丰富的面状地物，如建筑物、运动场、公园、湖泊等。参考影像、数字线划地图(DLG)、GIS矢量数据中存在大量的面特征，使用面特征进行几何校正能够充分地利用这些数据。另一方面，由于面特征是由许多点组成的，个别点的误差对于控制面本身来说影响并不大，因此，基于面特征的几何校正方法比基于点特征的几何校正方法具有更强的容错能力。

本发明提出了一种全新的基于面特征的图像几何校正方法，使用该方法能够获得高精度几何校正图像，该方法实用性很强，与影像的具体成像模型无关，适用于各种遥感影像、医学影像以及普通照片的几何校正与配准。

发明内容

本发明提出了一种点到多边形、多边形到多边形的距离度量方法，在此基础上将面状地物作为控制特征，建立一种基于面特征的几何校正模型，用于图像几何校正。

基于面特征的几何校正模型中关键的技术是利用面特征建立和求解误差方程，其中误差方程的建立需要计算控制面多边形与像平面多边形之间的距离，而多边形之间的距离是以多边形上的点到多边形的距离为基础的。因此，下面依次介绍点到多边形的距离计算方法、多边形到多边形的距离计算方法和基于面特征的几何校正模型建立方法。

1、点到多边形的距离计算方法

定义1：点p到多边形A的距离定义为点p到多边形A边界的最小距离，特别地，当点p在多边形A内部时，距离为0。

记多边形的边界点集合为

任意点p到多边形A的距离可表示为：

其中，x为多边形A的边界点，|x-p|为点p到点x的欧氏距离，

表示|x-p|的最小值。

欲求点p到多边形A边界的最小距离，首先分别计算点p到多边形A各条边的距离，然后得到其中最小值。其中，点p到多边形A某条边l的距离分为两种情况：过点p作线段l的垂线，若垂足在线段l上，则点p到垂足之间的距离为点p到l的距离；若垂足不在线段l上，则点p到线段l的两个端点的距离中较小值为点p到l的距离。

定义2：若(1.1)式中，距离值在x＝a处取得，则定义点p到点a的矢量为点p到多边形A的距离矢量，记为ρ(p，A)。特别地，如果ρ(p，A)＝0，则ρ(p，A)＝0。

算法1：点p到多边形A的距离矢量算法(算法流程见附图1)：

步骤1：令ρ为一个很大的值MAX，令l为多边形A的第一条边；

步骤2：计算点p到线段l的距离矢量d及距离d；

步骤3：如果d＜ρ，则令ρ＝d，ρ＝d；

步骤4：如果l是多边形A的最后一条边，则转到步骤5；否则，令l为多边形A的下一条边，并转到步骤2；

步骤5：输出点p到多边形A的距离矢量ρ。

2、多边形到多边形的距离计算方法

定义3：多边形A到多边形B的距离为多边形A所有边界点到多边形B的距离中最大者。

多边形A到多边形B的距离可以表示为：

$\mathrm{\ρ}(A,B)=\underset{p\∈\∂A}{\max }\mathrm{\ρ}(p,B)---\left(1.2\right)$

其中，p为多边形A的边界点，

表示ρ(p，B)的最大值。

利用该定义计算多边形之间的距离矢量时间复杂度较大，下面给出一个命题，可以简化多边形距离的计算。

命题1：多边形A到多边形B的距离为多边形A所有顶点到多边形B的距离中的最大者。

附图2说明了多边形和边界点的区别。可以从数学上证明命题1的正确性。

借助命题1，如果记多边形的顶点集合为

则多边形A到多边形B的距离可以表示为：

$\mathrm{\ρ}(A,B)=\underset{p\∈\tilde{A}}{\max }\mathrm{\ρ}(p,B)=\underset{p\∈\tilde{A}}{\max }\left(\underset{x\∈\∂B}{\min }\right|x-p\left|\right)---\left(1.3\right)$

其中max和min分别表示最大值和最小值，|x-p|为点p到点x的欧氏距离。

定义4：若(1.3)式中，多边形A到多边形B的距离值在p＝a，x＝b处取得，则定义点a到点b的矢量为多边形A到多边形B的距离矢量，记为ρ(A，B)。特别地，如果ρ(A，B)＝0，则ρ(A，B)＝0。

算法2：多边形A到多边形B的距离矢量算法(算法流程见附图3)：

步骤1：令ρ＝0，令p为多边形A的第一个顶点；

步骤2：计算点p到多边形B的距离矢量d及距离d(利用算法1进行计算)；

步骤3：如果d＞ρ，则令ρ＝d，ρ＝d

步骤4：如果p是多边形A的最后一个顶点，则转到步骤5；否则，令p为多边形A的下一个顶点，并转到步骤2；

步骤5：输出多边形A到多边形B的距离矢量ρ。

3、基于面特征的几何校正模型建立方法

经典的成像几何模型是用来建立地面点三维空间坐标与相应像点二维平面坐标之间的关系，一般地，通用模型可表示为：

$\left\{\begin{array}{c}x=f_x(X,Y,Z,t)\\ y=f_y(X,Y,Z,t)\end{array}\right.---\left(1.4\right)$

其中，(X，Y，Z)表示控制点的地面坐标，(x，y)表示控制点在影像上的量测坐标，t＝(t₁，t₂，…t_n)^T表示传感器几何校正模型的参数。目前，常用的几何校正模型有共线方程[Wong，K.W.，1980，Basic Mathematics of Photogrammetry，in Manual of Photogrammetry 4thEdition，Chapter II，Editor in chief：C.C.Slama，ASP Publishers，Falls Church，USA，pp.37-101.]、仿射变换模型[Okamoto，A.，Orientation Theory of CCD Line-scanner Images.InternationalArchives of ISPRS，1988.27(B3)：p.9.]、多项式模型[DE LEEUW，A.J.，VEUGEN，L.M.M.，andVAN STOKKOM，H.T.C.，1988，Geometric correction of remotely-sensed imagery using groundcontrol points and orthogonal polynomials.International Journal of Remote Sensing，9，1751-1759.]、有理函数模型[Tao，C.Vincent and Yong Hu.A Comprehensive Study of the RationalFunction Model for Photogrammetric Processing.Photogrammetric Engineering & Remote Sensing，Vol.67，No.12，December 2001，pp.1347-1357.]等。

表1给出了具体成像模型与通用模型的关系，虽然在不同的遥感影像成像模型中，f_x、f_y、t具有不同含义，但是这些模型均可用(1.4)式来表示。

表1 具体的成像模型与通用模型

本发明提出的方法可以用于各种几何校正模型。

面特征用多边形来描述。控制面由地面多边形A′((X_i，Y_i，Z_i)，i＝1，2，…，m)和影像多边形

组成，其中m为地面多边形A′的顶点个数，(X_i，Y_i，Z_i)为地面多边形A′各顶点的坐标，m₀为影像多边形A₀的顶点个数，

为影像多边形A₀各顶点的坐标。

利用基于点的几何成像模型(1.4)可计算出地面多边形A′在影像上所成的像多边形A，A的m个顶点坐标分别为(x_i，y_i)，i＝1，2，…，m，多边形A和A₀应满足以下约束关系：

ρ(A，A₀)＝0                 (1.5)

(1.5)式表示多边形A到多边形A₀的距离矢量为0。由命题1和定义4可知，多边形A到多边形A₀的距离必在多边形A的某一顶点处取得，不妨记为

其中k为整数且1≤k≤m，则：

ρ(a_k，A₀)＝0              (1.6)

矢量ρ(a_k，A₀)具有x和y两个分量，可记ρ(a_k，A₀)＝(ρ_x，ρ_y)，其中A₀为常量，且根据(1.4)式有

$\left\{\begin{array}{c}{x}_k=f_x(X_k,Y_k,Z_k,t)\\ y_k=f_y(X_k,Y_k,Z_k,t)\end{array}\right.---\left(1.7\right)$

其中(X_k，Y_k，Z_k)为地面多边形A′的第k个顶点坐标，则(1.6)式等价于(1.8)式：

$\left\{\begin{array}{c}{\mathrm{\ρ}}_x(f_x(X_k,Y_k,Z_k,t),f_y(X_k,Y_k,Z_k,t))=0\\ {\mathrm{\ρ}}_y(f_x(X_k,Y_k,Z_k,t),f_y(X_k,Y_k,Z_k,t))=0\end{array}\right.---\left(1.8\right)$

对(1.8)式进行线性化可以得到控制面特征的误差方程：

$\left\{\begin{array}{c}-v_x=(\frac{{\∂\mathrm{\ρ}}_x}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{\∂{\mathrm{\ρ}}_x}{\∂f_y}\frac{\∂f_y}{\∂t})\mathrm{\Δt}-l_x\\ -v_y=(\frac{{\∂\mathrm{\ρ}}_y}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{{\∂\mathrm{\ρ}}_y}{\∂f_y}\frac{{\∂f}_y}{\∂t})\mathrm{\Δt}-l_y\end{array}\right.---\left(1.9\right)$

其中，v_x和v_y为随机误差，l_x＝-ρ_x(f_x，f_y)，l_y＝-ρ_y(f_x，f_y)，Δt＝(Δt₁，Δt₂，…Δt_n)^T表示t的改正向量。

方程式(1.9)可表示成矩阵形式：

V＝AΔt-L         (1.10)

其中V＝(v_x，v_y)^T，L＝(l_x，l_y)^T， $A=\left[\begin{array}{c}\frac{{\∂\mathrm{\ρ}}_x}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{\∂{\mathrm{\ρ}}_x}{\∂f_y}\frac{\∂f_y}{\∂t}\\ \frac{{\∂\mathrm{\ρ}}_y}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{{\∂\mathrm{\ρ}}_y}{\∂f_y}\frac{{\∂f}_y}{\∂t}\end{array}\right].$

这里ρ_x和ρ_y并没有解析形式，所以采用数值计算的方法来逼近各函数偏导，如：

$\frac{{\∂\mathrm{\ρ}}_x}{\∂f_x}=\underset{\mathrm{\Δ}{f}_x\→0}{\lim }\frac{{\mathrm{\ρ}}_x(f_x+\mathrm{\Δ}{f}_x,f_y)}{\mathrm{\Δ}{f}_x}---\left(1.11\right)$

一个控制面特征可以导出如(1.9)式的2个误差方程，因此利用不少于n/2个控制面特征即可求解出模型中的n个未知参数t，参数的求解使用LM(Levenberg-Marquardt)算法[J.J.Mor6，″The Levenberg-Marquardt algorithm：implementation and theory，″in：G.A.Watson，ed.，Numerical Analysis Dundee 1977，Lecture Notes in Mathematics 630(Springer-Verlag，Berlin，1977)pp.105-116.]。

具体实施过程

下文结合说明书附图4，以卫星遥感影像为例，对本发明所提出的基于面特征的图像几何校正方法的实施方式加以说明。本发明所述方法包含但不限于所举实例。该方法实施的主要步骤如下：

步骤1：读取卫星遥感影像数据；

步骤2：输入控制面特征；

步骤3：利用卫星影像参数获取模型参数的初始值，当前迭代次数置为0；

步骤4：若当前迭代次数小于最大迭代次数，转至步骤5；若当前迭代次数达到最大迭代次数，转至步骤8；

步骤5：依次对各个控制面特征建立误差方程；

步骤6：将各控制面特征建立的误差方程联立，使用LM算法求出未知参数的改正值向量，并更新参数的值；

步骤7：用求得的改正值向量的模与规定的限差比较，若小于限差，转至步骤8；否则，当前迭代次数增加1，转至步骤4；

步骤8：解算完成，得到影像模型参数；

步骤9：利用模型参数对图像进行几何校正。

算法流程见附图4。

附图说明

图1是点到多边形距离矢量的求解过程示意图。

图2是多边形边界点和顶点的示意图。

图3是多边形到多边形距离矢量的求解过程示意图。

图4是基于面特征的图像几何校正过程示意图。

本发明方法的具体应用

应用实例一：使用安徽地区一景Landsat5TM 30米分辨率的多光谱数据和1∶25万DEM，进行两组实验。第一组试验中，选取9个控制面特征，用本发明方法进行正射校正；第二组试验中，在第一组试验中的9个控制面特征上各选取一个控制点，用这9个控制点进行正射校正。分别选取20个分布均匀的检查点验证两次校正的精度，结果如表2所示。

表2 使用控制面和控制点的校正精度对比

从表2中可知，在控制点质量不高的情况下，利用传统的基于控制点的方法进行校正，校正结果误差很大；将控制点替换为控制面，采用本发明的方法进行校正能够大大提高校正的精度。试验表明控制面比控制点具有更好的稳定性，受到粗差的影响相对较小。

应用实例二：使用武汉地区一景Landsat5TM 30米分辨率的多光谱数据和1∶25万DEM，进行两组试验。第一组实验中，选取10个控制点(分布不均匀，缺一角)，进行正射校正；第二组实验中，在缺少控制点的区域加入1个控制面特征，结合第一组实验中的10个控制点，用本发明方法进行正射校正。分别选取20个检查点验证两次校正的精度。表3列出了两次试验中后10个检查点的残差，17、18、19、20为缺少控制点区域内的四个检查点。

表3 面特征补充控制点进行正射校正的精度

从表3中可知，第一组实验中，在缺少控制点的区域，检查点有一个30米左右的固定偏差；而在第二组实验中，该区域的检查点精度偏差在10米以内。试验表明本发明用于补充控制点进行图像校正的有效性。

Claims (4)

1.一种基于面特征的图像几何校正方法，其主要步骤是：

步骤1：读取待校正卫星遥感影像数据；

步骤2：输入控制面特征；

步骤3：利用卫星影像参数获取模型参数的初始值，当前迭代次数置为0；

步骤4：若当前迭代次数小于最大迭代次数，转至步骤5；若当前迭代次数达到最大迭代次数，转至步骤8；

步骤5：依次对各个控制面特征建立误差方程；

步骤6：将各控制面特征建立的误差方程联立，使用LM算法求出未知参数的改正值向量，并更新参数的值；

步骤7：用求得的改正值向量的模与规定的限差比较，若小于限差，转至步骤8；否则，当前迭代次数增加1，转至步骤4；

步骤8：解算完成，得到影像模型参数。

步骤9：利用模型参数对图像进行几何校正。

2.根据权利要求1所述的方法，其特征是，在步骤5中，利用控制面特征建立的误差方程的方法如下：

经典的成像几何模型是用来建立地面点三维空间坐标与相应像点二维平面坐标之间的关系，一般地，通用模型可表示为：

$\left\{\begin{array}{c}x=f_x(X,Y,Z,t)\\ y=f_y(X,Y,Z,t)\end{array}\right.---\left(2.1\right)$

其中，(X，Y，Z)表示控制点的地面坐标，(x，y)表示控制点在影像上的量测坐标，t＝(t₁，t₂，…t_n)^T表示传感器几何校正模型的参数；目前，常用的具体成像模型有共线方程、仿射变换模型、多项式模型、有理函数模型等，表1给出了具体成像模型与通用模型的关系，虽然在不同的成像模型中，f_x、f_y、t具有不同含义，但是这些模型均可用通用模型来表示；

表1 具体的成像模型与通用模型

面特征用多边形来描述，控制面由地面多边形A′((X_i，Y_i，Z_i)，i＝1，2，…，m)和影像多边形

组成，其中m为地面多边形A′的顶点个数，(X_i，Y_i，Z_i)为地面多边形A′各顶点的坐标，m₀为影像多边形A₀的顶点个数，为影像多边形A₀各顶点的坐标；

利用基于点的几何成像模型(2.1)可计算出地面多边形A′在影像上所成的像多边形A，A的m个顶点坐标分别为(x_i，y_i)，i＝1，2，…，m，多边形A和A₀应满足以下约束关系：

ρ(A，A₀)＝0            (2.2)

其中，ρ(A，A₀)表示多边形A到多边形A₀的距离矢量；由多边形到多边形距离的定义可知，多边形A到多边形A₀的距离必在多边形A的某一顶点处取得，记为

其中，

为多边形A的顶点集合，k为整数且1≤k≤m，则：

ρ(a_k，A₀)＝0           (2.3)

其中，ρ(a_k，A₀)表示点a_k到多边形A₀的距离矢量；矢量ρ(a_k，A₀)具有x和y两个分量，记为ρ(a_k，A₀)＝(ρ_x，ρ_y)，其中，A₀为常量，且根据(2.1)式有

$\left\{\begin{array}{c}{x}_k=f_x(X_k,Y_k,Z_k,t)\\ y_k=f_y(X_k,Y_k,Z_k,t)\end{array}\right.---(2.4)$

其中，(X_k，Y_k，Z_k)为地面多边形A′的第k个顶点坐标，则(2.3)式等价于(2.5)式：

$\left\{\begin{array}{c}{\mathrm{\ρ}}_x(f_x(X_k,Y_k,Z_k,t),f_y(X_k,Y_k,Z_k,t))=0\\ {\mathrm{\ρ}}_y(f_x(X_k,Y_k,Z_k,t),f_y(X_k,Y_k,Z_k,t))=0\end{array}\right.---(2.5)$

对上式进行线性化可以得到控制面特征的误差方程：

$\left\{\begin{array}{c}-v_x=(\frac{{\∂\mathrm{\ρ}}_x}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{\∂{\mathrm{\ρ}}_x}{\∂f_y}\frac{\∂f_y}{\∂t})\mathrm{\Δt}-l_x\\ -v_y=(\frac{{\∂\mathrm{\ρ}}_y}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{{\∂\mathrm{\ρ}}_y}{\∂f_y}\frac{{\∂f}_y}{\∂t})\mathrm{\Δt}-l_y\end{array}\right.---\left(2.6\right)$

其中，v_x和v_y为随机误差，l_x＝-ρ_x(f_x，f_y)，l_y＝-ρ_y(f_x，f_y)，Δt＝(Δt₁，Δt₂，…Δt_n)^T表示t的改正向量；

方程式(2.6)可表示成矩阵形式：

V＝AΔt-L      (2.7)

其中，V＝(v_x，v_y)^T，L＝(l_x，l_y)^T， $A=\left[\begin{array}{c}\frac{{\∂\mathrm{\ρ}}_x}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{\∂{\mathrm{\ρ}}_x}{\∂f_y}\frac{\∂f_y}{\∂t}\\ \frac{{\∂\mathrm{\ρ}}_y}{\∂f_x}\frac{\∂f_x}{\∂t}+\frac{{\∂\mathrm{\ρ}}_y}{\∂f_y}\frac{{\∂f}_y}{\∂t}\end{array}\right]$

这里ρ_x和ρ_y并没有解析形式，所以采用下式数值计算的方法来逼近各函数偏导：

$\frac{{\∂\mathrm{\ρ}}_x}{\∂f_x}=\underset{\mathrm{\Δ}{f}_x\→0}{\lim }\frac{{\mathrm{\ρ}}_x(f_x+\mathrm{\Δ}{f}_x,f_y)}{\mathrm{\Δ}{f}_x}---\left(2.8\right)$

3.根据权利要求2所述的方法，其特征是，多边形到多边形的距离和距离矢量的定义及计算方法如下：

定义1：多边形A到多边形B的距离为多边形A所有边界点到多边形B的距离中最大者；

多边形A到多边形B的距离可以表示为：

$\mathrm{\ρ}(A,B)=\underset{p\∈\∂A}{\max }\mathrm{\ρ}(p,B)---\left(3.1\right)$

其中，p为多边形A的边界点，ρ(p，B)表示点p到多边形B的距离，

表示ρ(p，B)的最大值；

利用该定义计算多边形之间的距离矢量时间复杂度较大，下面给出一个命题，可以简化多边形距离的计算；

命题1：多边形A到多边形B的距离为多边形A所有顶点到多边形B的距离中的最大者；

借助命题1，如果记多边形的顶点集合为

则多边形A到多边形B的距离可以表示为：

$\mathrm{\ρ}(A,B)=\underset{p\∈\tilde{A}}{\max }\mathrm{\ρ}(p,B)=\underset{p\∈\tilde{A}}{\max }\left(\underset{x\∈\∂B}{\min }\right|x-p\left|\right)---\left(3.2\right)$

其中，max和min分别表示最大值和最小值，|x-p|为点p到点x的欧氏距离；

定义2：若(3.2)式中，距离值在p＝a，x＝b处取得多边形A到多边形B的距离值，则定义点a到点b的矢量为多边形A到多边形B的距离矢量，记为ρ(A，B)，特别地，如果ρ(A，B)＝0，则ρ(A，B)＝0；

多边形A到多边形B的距离矢量算法如下：

步骤3-1：令ρ＝0，令p为多边形A的第一个顶点；

步骤3-2：计算点p到多边形B的距离矢量d及距离d；

步骤3-3：如果d＞ρ，则令ρ＝d，ρ＝d；

步骤3-4：如果p是多边形A的最后一个顶点，则转到步骤5；否则，令p为多边形A的下一个顶点，并转到步骤3-2；

步骤3-5：输出多边形A到多边形B的距离矢量ρ。

4.根据权利要求2、3所述的方法，其特征是，点到多边形的距离和距离矢量的定义及计算方法如下：

定义3：点p到多边形A的距离定义为多边形A所有边界点到点p的最小距离，特别地，当点p在多边形A内部时，距离为0；

记多边形的边界点集合为

任意点p到多边形A的距离可表示为：

其中，x为多边形A的边界点，|x-p|为点p到点x的欧氏距离，

表示|x-p|的最小值；

欲求点p到多边形A边界的最小距离，首先分别计算点p到多边形A各条边的距离，然后得到其中最小的一个。其中，点p到多边形A某条边l的距离分为两种情况：过点p作线段l的垂线，若垂足在线段l上，则点p到垂足之间的距离为点p到l的距离；若垂足不在线段l上，则点p到线段l的两个端点的距离中较小的一个为点p到l的距离；

定义4：若(4.1)式中，距离值在x＝a处取得，则定义点p到点a的矢量为点p到多边形A的距离矢量，记为ρ(p，A)；特别地，如果ρ(p，A)＝0，则ρ(p，A)＝0；

点p到多边形A的距离矢量算法如下：

步骤4-1：令ρ为一个很大的值MAX，令l为多边形A的第一条边；

步骤4-2：计算点p到线段l的距离矢量d及距离d；

步骤4-3：如果d＜ρ，则令ρ＝d，ρ＝d；

步骤4-4：如果l是多边形A的最后一条边，则转到步骤5；否则，令l为多边形A的下一条边，并转到步骤4-2；

步骤4-5：输出点p到多边形A的距离矢量ρ。

Images