第二章 几何配准

第二章几何配准图像配准是遥感数字图像处理的重要内容，用于将不同时间、不同空间、不同视觉、不同场景下，不同传感器或者不同成像条件下的两幅或多幅图像进行叠加、拼接、对准、匹配等操作，以校正图像之间的平移、缩放、旋转、偏扭等几何差异和灰度差异。通常可将图像配准看作是特征空间、配准模型及相似性测度等几个要素的基本组合。特征空间指从图像中提取用于匹配的特征，配准模型是用于描述配准图像和参考图像之间空间关系的变换或映射函数，相似性测度用来估计待匹配特征之间的相似程度。由此可见，图像配准就是利用控制点数据，建立两幅之间的变换关系，将配准图像中的所有像元逐个变换到纠正空间中去。这样，如何实现地面控制点的自动提取，成了自动几何配准欲解决的关键问题，也是本文研究的重点所在。2.1特征提取图像特征是由于景物的物理与几何特性使图像中局部区域的灰度产生明显变化而形成的。因此，可以利用信息量的大小来判断图像特征的存在与否，然后提取特征点。提取出的特征点要求数量足够，尽量均匀分布在整幅图像，只有这样才能保证后续工作的正确率和准确度。我们将提取点特征的算子称为有利算子或兴趣算子，即运用某种算法从图像提取有利于某种应用目的或对之感兴趣的点。好的特征点的提取算法应该满足以下条件〔18〕：第一，特征点的重复度高，即在不同条件下拍摄的两幅图像中抽取的特征点具有较高的一致性与稳定性；第二，特征点所包含的信息量大，也就是说被提取的点很突出，能够有效地反映图像的主要结构信息；此外，特征点的提取算法还应具有旋转、平移的不变性，并且在发生小的尺度变化和透视形变时，能够检测出相同位置特征点，只有这样才能保证两幅图像的有效配准。比较知名的有利算子有Moraverc算子、Forstner算子、Hannah算子等。Moraverc算子以方差为判断依据，它对窗口中各个方向的方差进行统计，取最小方差都超过阈值的点为兴趣点，计算相对简单，但提取出的特征点精度相对不高;Hannah算子实质上与Forstner算子相似，而且计算较之复杂。考虑到AVHRR影像的分辨率和环境变化监测的需要，本文选用Forstner算子来提取影像的点特征。Forstner算子通过计算各像素的Robert'梯度和以像素（c，r）为中心的一个窗口的灰度协方差矩阵，在影像中寻找具有尽可能小而接近圆的误差椭圆的点作为特征点。由于Forstner算子比较复杂，可首先用一简单的差分算子提取初选点，然后采用Forstner算子在3X3窗口计算兴趣值，并选择备选点最后提取的极值点为特征点【16】。其具体实现步骤如下：（1）利用差分算子提取初选点差分算子为：计算像素（c,r）在上下左右四个方向的灰度差分绝对值dgl,dg2，dg3，dg4。dg1—gc,r-gc+1,rdg2—g-gc,rc,r+1>dg3—g-gc,rc-1,rdg4—gc,r-gc,r-1M=min{dgl,dg2,dg3,dg4}若对给定的阈值T,M>T,则(c,r)为一初选点，否则(c,r)不是特征点。也就是四个方向的差分绝对值有任意两个大于阈值，则该像素有可能是一特征点。在以初选点为中心的3X3窗口中逐像元计算Robert's梯度以组成协方差Q=N-1gg矩阵NQ=N-1gguvg2式中：c^k-c^k-1r^k-1g2=ui—c-kj—r-k(gi+1,j+1-g)2

i,j-gi+1,j)2c+k-gi+1,j)2g2— (gv i,j+1i—c-kj—r-kuv送刃(g-g)(g-g)uvi+1,j+1 i,j i,j+1 i+1,ji—c-k j—r-kk—INT(l/2),l为窗口大小

兴趣值：q—航，其中DetN代表矩阵N之行列式，trN为矩阵N之迹。可以证明，q即像素椭圆的圆度:(a2-b2)2q—1-(a2+b2)2a与b为椭圆之长、短半轴。如果a,b中任一个为零，则q=0,表明该点可能位于边缘上；如果a=b，则q=1,表明为一圆。确定阈值Tq,若q>Tq,则该像素为以备选点，按以下原则确定其权。0 (q<Tq)w—\DetN/ 、(q>Tq)ItrN(4)以权值为依据，大者为特征点。选取一适当窗口中的极值点为特征点，即选取窗口中权最

2.2特征匹配在提取出数量足够且分布均匀的特征点之后，就要在待匹配图像上搜寻与之对应的同名点(匹配点)，即确定两幅图像中对应同一实物的像点对。匹配点确定的方法很多，常用的有相关函数、协方差函数、相关系数、差平方和、差绝对值之和等。相关系数是标准化的协方差函数，协方差函数除以两信号的方差即得相关系数。它能消除两个灰度阵列之间的线性畸变，是一个最常用的相似性测度。如上图，在基准图像中以特征点为中心选取一大小为nXn的矩阵TW作为目标窗，在待配准图像中同名点的预测位置取一大一些窗口SW(mXm)作为搜索窗，TW在SW中每移动一个像元位置可求得一个相关系数，当移动完搜索窗SW中所有的位置时得到一个相关系数矩阵[P],选取具有最大相关系数对应的位置上的点作为匹配点。假设与目标窗口灰度矩阵相应的灰度函数为g(x,y)，(x,y)gTW，与搜索窗灰度矩阵相应的灰度函数为g'(x',y')，(x',y')gTW'，则g(x,y)与g'(x',y)的相关系数为() c(p,q)p(p,q)二CC(p,q)ggg'g'其中：C(P，q)= A {g(x,y)-E[g(x,y)]}{g'(x+p,y+p)-E[g'(x+p,y+q)]}dxdy(x,y)gTWC=J!{g(x,y)—E[g(x,y)]}2dxdygg(x,y)gTWC=g'g'(x+p,y+q)gTW'JJ{g'(x+C=g'g'(x+p,y+q)gTW'其中：1g(x,y)gTWE[g(x,y)]g(x,y)gTWITWI为目标窗TW的面积

1E[1E[g'(x+p，y+q)]二阿ag(x,y)e7W'g'(X+p,y+q)dxdy若p(p,q)>p(p,q)(p丰p,q丰q)，则p,q为搜索区影像相对于目标区000000影像的位移参数。由离散灰度数据对相关系数的估计为p(c,p(c,r)二工艺(g-g)(g' -g')i,j i+r,j+cr,cg二丄工Kgn2 i,ji=ij=1a=aa=a，Ni

i=1b= 工bNii=1g'=——g'r,c n2 i+r,j+ci=1j=1若p(c,r)>p(c,r)(c丰c,丫丰r)，则c,r为搜索区影像相对于目标区影像000000位移的行、列数参数。下面从几何角度分析相关系数。为了表示的方便，假定A和B为一维阵列,且A和B包含N个像元灰度值，则A和B分别表示N维空间中的一个矢量。B=(b,b2…bt)1NA=a—a，B=b—bii ii如果A和B为相同图像，则A=B.协方差：C=工AB=AB=ABcos0AB ii 11i=1其中||A||和||B|分别表示矢量A和矢量B的模，即欧式范数：0为矢量A和B的夹角。下图以最简单的二维矢量形式表示了矢量A和矢量B之间的相关系数的几何意义。

几何意义。相关系数：p=iA侖=cmax'相关系数越大，矢量A和矢量B的夹角9越小，相似性程度越高。相关系数虽然计算比较复杂，但却是一种较好的匹配测度。通过上面的分析不难发现，相关系数取值范围满足：P<1。经过特征匹配得到一定数量的同名点对，但不是所有的同名点对都是所要求的控制点。此时可按照相关系数的大小，将同名点对从大到小排列，选取前面相关系数较大一部分作为控制点。2.3配准模型的建立通过特征匹配确定了控制点对，然后利用这些点对之间的关系确定几何变换模型。常用的方法有三角形线性法和多项式配准模型等。三角形线性法的前提是假定小区域内几何畸变是线性的，同时，为了提高纠正精度，划分区域应越小越好，但是实际工作中并不能保证GCP分布均匀，特别是针对空间分辨率低的AVHRR数据，三角形线性法并不常用。多项式配准模型是实践中经常使用的一种方法。多项式模型的基本思想是回避成像的空间几何过程，直接对图像变形本身进行数学模拟。多项式模型将配准过程中的总体形变看作是平移、缩放、旋转、仿射、偏扭、弯曲以及更高次的基本变形的综合作用的结果，因而配准图像和参考图像对应点之间的坐标关系可以用一个适当的多项式来表达。多项式纠正的实质是用二元多项式来描述几何畸变，其数学模型为：x二f(u,v)二另乞auivjx ij< i=0j=0y=f(u,v)=另艺buivj

y ij、 i=0j=0式中i=0,1,2,...,n、j=0,1,2,...n—i,为二元多项式的次数；x,y为原始图像空间中的像元坐标；u,v为纠正图像空间网格坐标；n为二元多项式次数；a,ijb为二元多项式系数。ij利用提取出来的控制点数据运用最小二乘法进行曲面拟合求系数，假设得到的控制点对数为M(M>(n+1)(n+2)/2),令：艺(Xmm=1=Min艺(Xmm=1=Minauiv j)ijmmi=0j=0对个系数求一阶偏导并令其为0可以得到关于未知系数的联立线性方程组，解方程组可得到多项式系数，从而建立了纠正空间与标准空间位置转换的数学模型。根据这个数学模型当已知纠正空间某一位置时，可反求出其在畸变空间的位置。本文的重点在于控制点的选取，而控制点数据是用来求取最佳系数，以建立有效的多项式配准模型，从而将畸变空间中所有的像元更好的变换到基准空间中去。由于纠正后图像和原始图像的形状、大小、方向都不一样。所以在纠正过程实施之前，必须首先确定新图像的大小范围。

根据已建立的多项式配准模型求出原始图像四个角点(a,b,c,d)在纠正图像中的对应点(a',b',c',d')的坐标(Xa',Ya')(Xb',Yb')(Xc',Yc')(Xd',Yd')，求出最大值和最小值，即：X1=min(Xa',Xb',Xc',Xd')X2=max(Xa',Xb',Xc',Xd')

Y1=min(Ya',Yb',Yc',Yd')

Y2=max(Ya',Yb',Yc',Yd')然后，根据精度要求，在新图像范围内划分网格，每个网格就是一个像元，任意一个像元由其行列号唯一确定，新图像的行数M=(Y2-Y1)/^Y+1;新图像的列数N=(X2-X1)/AX+1。2.4灰度的重采样纠正后得到的新图像的每一个像元，根据变换函数，可以得到它在原始图像上的位置。如果求得的位置为整数，则该位置处的像元灰度就是新图像的灰度值。但是计算后的点多半不在原图的像元中心处，因此必须重新计算新位置的灰度值，通常有三种方法：最近邻法、双向线性内插法和三次卷积法。最近邻法图像中两相邻点的距离为1,即行间距厶x=1,列间距△『=】，取与所计算点(x,y)周围相邻的4个点，比较它们与被计算点的距离，哪个点距离最近，就取哪个的亮度值作为(x,y)点的亮度值f(x,y)。设该最近邻点的坐标为(k,l),则k=Integer(x+0.5)l=Integer(y+0.5)有：f(x,y)=f(k,l)最邻近内插法以距内插点最近的观测点的像元值为所求的像元值。该方法最大可产生0.5个像元的位置误差,优点是不破坏原来的像元值,处理速度快。双线性内插法以实际位置临近的4个像元值,确定输出像元的灰度值。公式为：为f.P.iif(x，y)二专Zp.ii=1式中f(x,y)为输出像元灰度值，f为邻i近点i的灰度值。p为邻近点对投影i点的权重(pi=1/d.，d.表示邻近点到投影点的距离，最近者权重最大)。双线性内插法使用内插点周围的4个观测点的像元值,对所求的像元值进行

线性内插。缺点是破坏了原来的数据，但具有平均化的滤波效果三次卷积法增加邻点来获得最佳插值函数，以实际位置临近的16个像元值，确定输出像元的灰度值。公式为：f-1z-lJ+2J-17+2a(x,y)M 》( 3计1J-l汁1j+2i+2y-i计2J+2艺f.p.iiJ=1 艺piii=1式中f(x,y)为输出像元灰度值，f为邻近点i近点i的灰度值。p为邻近点对投影

i点的权重(pi=1/J,d.表示邻近点到投影点的距离，最近者权重最大)。因这种三次多项式内插过程实际上是一