C++编程实现遥感影像几何纠正及效果分析
摘要 对于卫星遥感获取的影像作为空间数据，具有空间地理位置的概念。而影像在获取过程中总存在误差，因此，遥感影像的几何纠正是遥感信息处理中的一个重要环节，是减小遥感影像与地面真实形态差异的重要手段。本文就遥感影像几何纠正的C++程序实现算法进行研究和分析。
关键词 多项式系数 几何纠正 参照系统 坐标位置 灰度重采样
一、原理介绍
在遥感理论上,将遥感平台位置和运动状态、地形起伏、地球表面曲率、大气折射等遥感系统内外因素影响造成的遥感图像几何位置上的变化统称为几何畸变，也就是遥感图像在几何位置上与实际地面位置有差异。在图像上表现为像元行列分布不均匀，像元大小与地面大小对应不准确等。针对不同因素，可以采取相应的纠正方法。
本文所探讨的基于多项式的几何校正，适用于地面平坦,不考虑高程信息，或地面起伏较大而无高程信息的情况。方法机理是利用地面控制点的图像坐标和其同名点的地面坐标通过平差原理计算多项式中的系数，然后用该多项式对图像进行纠正。一般参照高分辨率的相应带坐标信息的标准影像，通过待纠正影像上和标准影像上的同名点建立多项式的对应关系，反算解求多项式系数，再利用解出的多项式对待纠正影像进行灰度重采样，实现遥感图像与实际地理图件间的配准，达到消减以及消除遥感图像的几何畸变的目的。
二、算法设计
1、算法概述
本文研究的几何纠正算法的数学基础是基于多项式的，首先通过寻找待纠正影像与对应参照坐标空间的同名点，将同名点的坐标代入多项式方程，通过平差计算得到待定多项式的系数，从而确定待纠正影像到参照系统的一种映射关系。
然后，利用解算的多项式方程求出纠正后影像的图幅范围以及与参照系统尺度的比例系数，新建该图幅大小的图像。依据同理解算的多项式方程的逆变换公式，以遍历新影像的方式，对应地从待纠正影像获取新影像的灰度值，从而实现影像的几何纠正。
2、流程图
从待纠正影像的输入，到生成纠正后影像的操作过程及其算法流程图如下：
图-1
三、实现方法与过程
影像的多项式纠正的原理是利用原图与参照系统的同名点的位置信息，解算二者之间的多项式对应关系，因此关键是二者之间控制点信息的获取。
这里采取了两种方案，一是手动采点，打开待纠正影像和参照的标准影像，通过手动获取同名点坐标信息解算多项式系数并灰度重采样；一是读入控制点坐标文件，直接解算多项式系数并进行影像的几何纠正。
以下是研究影像几何纠正的程序执行过程以及相应的C++实现代码：
1、读入待纠正影像
打开待纠正影像，获取影像的几何信息以及灰度信息。若是手动选点，还需打开参考影像，为了方便对照选取同名点，采用创建一幅新图像，同时显示两幅影像。
新建图像同时显示两幅影像时，其灰度赋值的主体实现代码如下：
int cols=pDoc->m_bmpfile.m_Cols; //获得原图的信息
int raws=pDoc->m_bmpfile.m_Rows;
new_bmpfile=pDoc->m_bmpfile; //复制原图
std_bmpfile.Load4File(dlg.m_path); //打开标准图像
int scols=std_bmpfile.m_Cols; //获得标准图像数据
int sraws=std_bmpfile.m_Rows;
int r,c;
int arow=raws>sraws?raws:sraws;
int acol=cols+scols+50;
all_bmp.CreateBmp(acol,arow,1); //创建图像实现两幅影像的预览
for(r=0;r<arow;r++)
for (c=0;c<acol;c++)
{
if(r<raws&&c<cols)
{
all_bmp.m_pImgDat[r*acol+c]=new_bmpfile.m_pImgDat[r*cols+c];
}
else if(r<sraws&&c>=cols+50)
{
all_bmp.m_pImgDat[r*acol+c]=std_bmpfile.m_pImgDat[r*scols+(c-cols-50)];
}
else
all_bmp.m_pImgDat[r*acol+c]=255;
}
2、获取同名点坐标
方式一：读入同名点数据文件
直接输入文件路径，读取事先确定好的同名点坐标文件，得到待纠正影像和参照系统上对应控制点的坐标值。
方式二：手动采集同名点
在一张新建影像上同时显示两张影像后，通过定义鼠标左键函数，实现采集并存储同名点的坐标信息。
鼠标左键函数的C++主体代码如下：
CPoint Org=GetScrollPosition(); //获得滚动视窗的滚动量
Group[pressnum].y=(point.x+Org.x); //将屏幕坐标转化为图像内的像素坐标并存储
Group[pressnum].x=pDoc->m_bmpfile.m_Rows-(point.y+Org.y);
pressnum++;
if((pressnum)%2!=0&&Group[pressnum-1].y>new_bmpfile.m_Cols) //纠错提示功能
{
AfxMessageBox("该选点无效！请在待纠正影像上选点！");
pressnum--;
}
else if((pressnum)%2==0&&Group[pressnum-1].y<(new_bmpfile.m_Cols+50)) //纠错提示
{
AfxMessageBox("该选点无效！请在参照影像上选点！");
pressnum--;
}
lable(pDoc->m_bmpfile,Group[pressnum-1]); //选点标记函数
3、解算多项式的正逆变换系数
多项式模型的次数可以有多个，这里以一次多项式系数的解求方法为例予以说明。
（1）输入控制点坐标构成系数矩阵
正变换C++代码：
for(i=0,j=0;i<Num;i+=2,j++) //给正变换系数矩阵A和常数矩阵L赋值
{
A[i][0]=1; A[i][1]=a[j].x; A[i][2]=a[j].y; A[i][3]=0; A[i][4]=0; A[i][5]=0;
A[i+1][0]=0; A[i+1][1]=0; A[i+1][2]=0; A[i+1][3]=1; A[i+1][4]=a[j].x; A[i+1][5]=a[j].y;
L[i][0]=sa[j].x; L[i+1][0]=sa[j].y;
}
逆变换C++代码：
for(i=0,j=0;i<Num;i+=2,j++) //给逆变换系数矩阵A和常数矩阵L赋值
{
A[i][0]=1; A[i][1]=sa[j].x; A[i][2]=sa[j].y; A[i][3]=0; A[i][4]=0; A[i][5]=0;
A[i+1][0]=0; A[i+1][1]=0; A[i+1][2]=0; A[i+1][3]=1; A[i+1][4]=sa[j].x; A[i+1][5]=sa[j].y;
L[i][0]=a[j].x; L[i+1][0]=a[j].y;
}
（2）最小二乘平差解算正逆变换系数
正变换C++代码：
p=dlg.zhuan(A,Num,6); //矩阵转置
q=dlg.MatrixMul(p,6,Num,A,Num,6); //矩阵相乘
rr=dlg.InverseMatrix(q,6); //矩阵求逆
s=dlg.MatrixMul(rr,6,6,p,6,Num);
t=dlg.MatrixMul(s,6,Num,L,Num,1); //得出平差后的解算结果
a0=t[0][0]; //赋值给正变换系数
a1=t[1][0];
a2=t[2][0];
b0=t[3][0];
b1=t[4][0];
b2=t[5][0];
逆变换C++代码：
p=dlg.zhuan(A,Num,6); //矩阵转置
q=dlg.MatrixMul(p,6,Num,A,Num,6); //矩阵相乘
rr=dlg.InverseMatrix(q,6); //矩阵求逆
s=dlg.MatrixMul(rr,6,6,p,6,Num);
t=dlg.MatrixMul(s,6,Num,L,Num,1); //得出平差后的解算结果
pa0=t[0][0]; //赋值给逆变换系数
pa1=t[1][0];
pa2=t[2][0];
pb0=t[3][0];
pb1=t[4][0];
pb2=t[5][0];
4、创建新图并依据多项式方程灰度重采样实现影像纠正
实现灰度重采样的主要C++代码如下：
int stand_drow=0,stand_dcol=0;
int stand_xrow=1000,stand_xcol=1000;
for(r=0;r<raws;r++) //确定纠正后影像的图形区大小
for (c=0;c<cols;c++)
{
r1=int(a0+a1*r+a2*c);
c1=int(b0+b1*r+b2*c);
if(r1>stand_drow)
stand_drow=r1; //r1从0开始的数组脚码，所以下面计算真正的行数时要加1
if(r1<stand_xrow)
stand_xrow=r1;
if (c1>stand_dcol)
stand_dcol=c1;
if(c1<stand_xcol)
stand_xcol=c1;
}
int stand_row=int((stand_drow-stand_xrow)/k+0.5)+1; //得到纠正后影像的行列数
int stand_col=int((stand_dcol-stand_xcol)/k+0.5)+1;
standbmp.CreateBmp(stand_col,stand_row,1); //创建新图
for(r=0;r<stand_row;r++)
for (c=0;c<stand_col;c++)
{
standbmp.m_pImgDat[r*stand_col+c]=0;
}
double X,Y;
for(r=0;r<stand_row;r++) //以遍历新影像的方式灰度重采样
for (c=0;c<stand_col;c++)
{
X=k*r+stand_xrow; //还原由待纠正影像解算的标准图的尺度
Y=k*c+stand_xrow;
r1=int(pa0+pa1*(r+stand_xrow)+pa2*(c+stand_xcol)); //反算回待纠正影像并获取器灰度值
c1=int(pb0+pb1*(r+stand_xrow)+pb2*(c+stand_xcol));
if(r1>=0&&r1<raws && c1>=0&&c1<cols)
standbmp.m_pImgDat[r*stand_col+c]=new_bmpfile.m_pImgDat[r1*cols+c1];
else
standbmp.m_pImgDat[r*stand_col+c]=0;
}
四、实验与分析
1、操作界面：
图-2
由操作界面可以看出：影像的几何纠正可以采用手动采集同名点和输入同名点文件两种方式进行。
2、手动采集同名点
图-3
定义鼠标左键函数只能获取以客户区右上角为原点的平面坐标，而影像的首行首列是以影像的左下角为起始点的，这就需要进行坐标系间的转换，才能得到当前影像的实际行列数，这时依此行列坐标给选择的同名点用十字圆形图标予以标记显示。最后根据左右影像在整张图中的显示位置，还原以各自影像左下角为首行首列时的同名点行列位置，并存储下来。
3、引入数据确定同名点
图-4
4、纠正效果
（1）手动选点方式的纠正效果：
图-5（标准图） 图-6（待纠正图）
图-7（一次纠正模型同名点分布）
图-8（使用一次纠正模型纠正后）
图-9（二次纠正模型同名点分布）
图-10（使用二次纠正模型纠正后）
由以上效果图可知，手动获取同名点的几何纠正算法影像的几何变形起到了纠正作用，当待纠正影像变形严重且分辨率不太差时，二次纠正模型较一次纠正模型的纠正效果好，但是从几何纠正的原理上来说，当待纠正影像分辨率不太高以致影响同名点的选取精度时，选取一次纠正模型较好。
几何纠正实际上是根据同名点在标准影像以及待纠正影像上的位置分布特征，使待纠正的整张影像按照标准影像上图形元素的分布，通过灰度重采样实现影像的几何纠正。
（2）数据输入方式的纠正效果：
图-11（纠正前）
图-12（一次纠正后） 图-13（二次纠正后）
数据引入方式与手动采集方式只是同名点的获取方式不同，几何纠正的原理和过程是一样的，这里的同名点数据固定为7组，一次和二次纠正模型的纠正效果相差不大，但是相对纠正前整张影像明显拉长。由于没有标准参照影像，参照对象是待纠正影像上控制点的地面绝对坐标，所以无法评价影像的纠正效果。
五、结束语
遥感影像的几何纠正的理论性虽然很强，但是其算法思想不难理解，我们必须深入理解其机理，这样才能对控制点的选择、控制点的分布和数量要有很好的把握，才能使纠正后的影像有较好的效果。而编程实现遥感影像的几何纠正可以帮助我们更深地理解几何纠正的本质，不仅是对理论知识的一次升华，更是为更熟练地操作遥感影像处理软件打下的的重要基础。
几何纠正的算法思想在此就不再赘述，主要列举一下在C++编程实现过程中应该注意的问题：
1、定义鼠标消息函数获取图像上所选同名点的行列位置时，系统记录鼠标位置的是CPoint类的对象point，point的值所在坐标系是以客户区左上角为原点，向下向左为坐标轴正轴方向的坐标系，而解算多项式系数需要得到的像素坐标是以影像的首行首列位置为原点，向上向右为坐标轴正方向的坐标系中的坐标。即获得的鼠标位置需要进行转换才能得到影像上同名点的对应像素坐标。
2、多项式系数解求过程中，正逆变换的多项式系数都需要求解，正变换用于根据待纠正影像的图幅确定新图像的图幅大小，反变换用于新图像到待纠正影像的灰度重采样。另外，在确定新图图幅大小的过程中，由待纠正影像根据正变换公式得到的尺寸是与其参照影像或者说参照坐标系统的尺寸一致，显然要得到的纠正影像应该与纠正前大小一致，所以事先要确定待纠正影像与参照影像或是参照坐标系统之间的尺寸比列关系。
3、实现视窗滚动功能以后，需要注意，鼠标位置获取的坐标一直是以客户区左上角为原点的，所以在手动选取影像同名点的过程中，为了避免因滑动滚动条而得到错误的坐标信息，应在鼠标消息函数里添加滚动量的代码，将其纠正为真实的坐标值。
4、灰度重采样时，根据对应关系遍历新建影像获取灰度信息的效果比遍历待纠正影像根据对应关系给新图赋灰度值更好，因为后者容易使纠正后的新影像出现因在待纠正影像中无对应的像素采集灰度而出现斜纹噪声。
总之，遥感影像几何纠正的编程实现首先需要深入理解的理论知识，然后是
编程中遇到实际问题的解决能力。编程练习对更深入地理解理论知识有着非常大的帮助，在经历的遇到问题、分析问题，解决问题的过程中，能学到之前未曾遇到过的问题，理论实践过程能让学习者受益匪浅。
参考文献
【1】王学平.遥感影像几何校正原理及效果分析.计算机应用与软件，2008.9，第25卷第9期
【2】孙家抦.遥感原理与运用.第二版.湖北武汉：武汉大学出版社,2009.6
成绩评定
