CN103218783B - 基于控制点影像数据库的卫星遥感图像快速几何纠正方法 - Google Patents

Abstract

基于控制点影像数据库的卫星遥感图像快速几何纠正方法，包括：确定待纠正影像的地理范围；针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索所有符合要求的控制点；进行自动匹配选点，确定用于几何纠正的控制点对；判断匹配的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回所述检索步骤；根据传感器类型和成像特点，选择几何纠正模型；基于匹配的控制点对和选择的几何纠正模型建立像素坐标与大地坐标的转换关系；基于数字高程模型和所述转换关系进行几何纠正，得到经过纠正的数字正射影像。本图像纠正方法，将影像自动匹配技术引入到控制点的自动选取中，可实现遥感影像的自动几何纠正处理。

Description

基于控制点影像数据库的卫星遥感图像快速几何纠正方法

技术领域

本发明涉及一种遥感图像纠正方法，更具体而言涉及一种基于控制点影像数据库的卫星遥感图像快速几何纠正方法。

背景技术

随着遥感技术的发展，特别是遥感传感器技术的不断发展，通过遥感技术所获得的遥感影像或数据的用途越来越广。目前，遥感数据的应用范围已经扩展到社会信息服务领域，例如，广泛应用于测绘、农业、林业、地质矿产、水文与水资源、环境监测、自然灾害、区域分析与规划、军事、土地利用等方面。具有精确地理编码的遥感影像可以为土地、规划、环保、农业、林业、海洋等不同的领域提供各自需要的地物特征和信息。

在通过卫星或航空平台等飞行平台获取遥感图像数据或其他数据时，会受到天气、日光、遮挡等外在因素的影响，同时，在数据采集时飞行平台的高度、姿态会发生变化，因此，在进行遥感图像拍摄时往往会造成图像平移、旋转、缩放等问题。此外，根据光学成像原理，相机成像时是按照中心投影方式成像的，因此地面上的高低起伏在成像时就会导致投影差的存在。上述因素综合，会造成遥感影像的误差，例如倾斜误差、投影误差等。因此，在使用这些遥感影像/数据之前需要对所获得的原始遥感影像进行正射纠正。

传统的遥感影像正射纠正一般包括：首先通过外业测量或者从已有的地形图资料采集地面控制点（GroundControlPoint，GCP）及数字高程模型（DigitalElevationModel，DEM）信息；然后将这些信息导入专业的遥感或者数字摄影测量***；接下来由加载了GCP和DEM信息的***对遥感影像进行正射纠正。

地面控制点GCP(GroundControlPoints)是卫星遥感影像几何纠正和地理定位时重要的基准数据源。在遥感影像几何纠正处理过程中，为达到一定的纠正精度，一定数量的地面控制点是必不可少的，通过控制点的物方坐标和相应的像点坐标构建影像成像模型并求解模型参数或对已有成像模型进行优化补偿求解补偿参数，最终建立成像过程中物方与像方的正确转换关系。

传统作业过程中的控制点一般采用全野外测量的模式，需要经过“收集已有控制资料、测区踏勘、选点埋石、外业测量、内业整理”等一系列复杂的工序，虽然随着GPS-RTK等先进测量技术的发展，外业测量的工作量和复杂度大大减少，但必要的控制点野外测量工作仍不可避免。

此外，由于卫星遥感影像覆盖范围大（以资源三号卫星为例，覆盖范围为50公里×50公里），要想获得均匀分布的地面控制点，通常需要在几百公里乃至上千平方公里的范围内施测，其外业工作量之大、耗费人力和物力之巨是不言而喻的。此外，在地震、洪水、泥石流等自然灾害多发区域或原始森林、沼泽、沙漠等人迹罕至地区，测量人员常常无法进入实地进行测量。

即使控制点是通过人工选取遥感影像和地形图同名点来获得，也存在效率低下、精度难以保障等问题。而对同一地区的不同时相或不同传感器影像进行纠正，又会出现重复选点的情况。

发明内容

根据本发明的实施例，提供了一种基于控制点影像数据库的卫星遥感图像快速几何纠正方法。该方法包括：确定待纠正影像的地理范围；针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索所有符合要求的控制点；进行自动匹配选点，确定用于几何纠正的控制点对；判断所述匹配的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回所述检索步骤；根据传感器类型和成像特点，选择相应的几何纠正模型；基于所述匹配的控制点对和所述选择的几何纠正模型建立像素坐标与大地坐标的转换关系；基于数字高程模型和像素坐标与大地坐标之间的转换关系进行几何纠正，得到经过正射纠正的数字正射影像。

根据本发明实施例的遥感图像纠正方法，将影像自动匹配技术引入到控制点的自动选取中，可以实现遥感影像的自动几何纠正处理。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍，显而易见地，下面描述中的附图仅仅涉及本发明的一些实施例，而非对本发明的限制。

图1是根据本发明的基于控制点影像数据库的卫星遥感图像快速几何纠正方法的流程图；

图2a-2c是示例性的控制点影像片与资源三号卫星影像的对比图；

图3是控制点影像片检索的一种示意性流程图；

图4是控制点自动匹配方法的一种示意性流程图；

图5示出了不同尺度空间下的二维影像的一组例子；

图6示出了高斯差分尺度空间（DoG）影像的一组例子；

图7是DoG尺度空间局部极值检测的示意图；

图8是由关键点邻域梯度信息生成特征向量的示意图；

图9示出了不同影像由关键点邻域梯度信息生成特征向量图的一组例子；

图10是基于SIFT描述子的匹配算法的流程的示意图；

图11是包括了误匹配点剔除步骤的匹配算法流程的示意图；

图12示出了线阵推扫式成像传感器的成像方式；

图13示出了几何纠正中的像点与地面点的数学关系；

图14示出了严格成像模型可能用到的坐标系；

图15是成像方向与地球椭球的交会的示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例的附图，对本发明实施例的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于所描述的本发明的实施例，本领域普通技术人员在无需创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

除非另作定义，此处使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明专利申请说明书以及权利要求书中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。同样，“一个”或者“一”等类似词语也不表示数量限制，而是表示存在至少一个。

根据本发明实施例的基于控制点影像数据库的卫星遥感图像快速几何纠正方法，是以控制点影像数据库的建立为基础，通过控制点影像数据库的检索、控制点自动匹配以及目标图像纠正来实现的。所述控制点影像数据库是控制点影像片的集合。而控制点影像片是传统控制点的拓展，也就是说，用具有地理坐标信息的正射影像DOM块和相应的数字高程模型DEM块取代单一的控制点作为数据库的数据单元，从而生成控制点影像数据库。

图1示意性地示出了根据本发明实施例的遥感图像几何纠正方法的流程图。在在步骤S1，确定待纠正影像的地理范围；在步骤S3，针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索出所有符合要求的控制点；在步骤S5，进行自动匹配选点，确定用于几何纠正的控制点对；在步骤S7，判断所获得的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回步骤S3继续进行检索；在步骤S9，根据传感器类型和成像特点，选择相应的几何纠正模型类别；步骤S11，基于控制点对、几何纠正模型建立像素坐标与大地坐标的转换关系；在步骤S13，基于数字高程模型(DEM)和像素坐标与大地坐标之间的转换关系,对原始影像进行几何纠正，得到经过正射纠正的数字正射影像。

可选地，在进行步骤S1之前，对原始待纠正影像进行影像粗纠正。该粗纠正可以采用非基于控制点影像数据库的常规方法，例如，采用原始影像的元数据信息中记录的四个角点大地坐标和影像坐标，构建一次多项式模型，完成概略的几何纠正，获得没有经控制点改正的纠正影像。通过粗纠正可有效预测控制点影像片对应的原始影像的范围，从而缩小影像匹配的搜索范围，减少匹配时间。因此，采用该影像粗纠正可以提高本发明的几何纠正的速度。

根据本发明实施例的上述匹配算法，在一定误差范围内，可实现快速自动化选点操作。可选地，在步骤S5中，采用人工辅助筛选与自动匹配相结合，可以提高控制点的选择质量，进而提高几何纠正的精度。

另外，可选地，在步骤S13之后，进行纠正结果精度检查，判断纠正结果是否满足要求，如果不满足要求，则对控制点的选择结果进行调整。该检查结果可作为纠正结果的几何精度评价参考指标。

以下分部分描述根据本发明的纠正方法的具体技术内容。

控制点影像数据库构建

现代遥感技术的快速发展，遥感图像获取的周期与精度也逐渐得到了提高，这为控制点影像数据库的建立提供了可靠的资料基础。而计算机自动识别技术的发展，则为控制点数据库的高效使用提供了技术基础。根据本发明的实施例，建立了控制点影像数据库，将控制点的属性信息和影像信息进行统一建库管理，实现了“一次建库，部分更新，多次使用”的目标。同时，将影像自动匹配技术引入到控制点的自动选取中，实现遥感影像的自动或半自动几何纠正处理。

在一个实施例中，每一个控制点影像片都包含两种数据：图像数据和属性数据。其中属性数据用来描述地理位置信息，包括四方面内容：

①描述GCP地理位置信息，如三维坐标X、Y、Z；

②描述地理坐标的一些必要的辅助信息,如所采用的坐标系、投影方式、椭球参数等；

③描述控制点影像的辅助信息,如传感器的类型、波段、图像宽高、图像分辨率等；

④GCP选取的特征描述,如道路的交叉口或桥梁中心点，这些信息可以作为查询的附属条件。

控制点影像数据库对大量控制点影像片采用数据库的方式进行存储、管理和服务。传统的方法在数据存储管理上是把控制点的属性和图像数据割裂开来，在数据库中存储的是点对应图像的文件指针，而图像数据则在数据库外部以文件方式单独存储，这种方式破坏了控制点信息的完整性和数据库的安全性，极易由于文件的误删除而丢失控制点的图像信息。根据本发明的控制点影像数据库采用将图像与属性结合在一起的数据存储形式，即对控制点影像片对应的图像信息采用大二进制对象BLOB(BinaryLargeObject)类型并作为表结构的一个字段与控制坐标、椭球类型、投影方式等属性进行一体化存储管理。

对已有的1：10000尺度的DOM（数字正射影像图，DigitalOrthophotoMap）、DEM（数字高程模型，DigitalElevationMode）成果图进行整理，选取现势性强、纹理清晰、特征明显的地区，如道路交叉口、桥梁、田埂角点等地区按一定大小（如512×512像素）采集影像片，从DOM中获得平面坐标信息、椭球、投影等信息，同时从对应的DEM数据中获取相应区域的高程值，然后将获得图像信息和属性信息进行统一入库存储。

图2a是从某地的DOM影像上裁剪出的一个控制点影像片（512×512像素，分辨率为1米）；图2b为同一地区的资源三号的正视全色影像（分辨率为2.1米），图2c为同一地区的资源三号正视多光谱影像（分辨率为5.8米）。对照图2b和图2c可以很直观地从图2a中选取相应的同名点作为控制点。

此外，控制点影像数据库的覆盖范围和存储量到一定规模后，如果按常规的物理顺序进行逐一检索，则会耗费大量时间，不利于实际应用，为了快速检索出控制片，在实际应用中，需要对控制点数据库按地理坐标进行分区存储。在分区过程中，根据地形起伏、地物的复杂程度以及具体的应用，对感兴趣的地区及控制点分布较密的地区，减小划分范围，而对于次要地区或大面积水域等扩大划分范围，最终保证每一区域内的控制点数量基本保持一致，加快检索的速度。

控制点影像数据库的检索

检索功能是衡量一个数据库***的重要技术指标。对于控制点影像数据库来说，在进行几何纠正时,能否根据待纠正影像快速检索出可利用的控制点影像是应用***建设的成败关键之一。

图3是步骤S3中控制点影像片检索的一种示意性流程图。在步骤S301，基于步骤S1中所估算的概略地理位置范围，进行基于目标区域中心点经纬度的检索；在步骤S303，根据控制点影像片的属性信息进行筛选；在步骤S305，进行基于内容的高级检索；在步骤S307，获得所需的控制点影像片。下面分别介绍控制点影像片检索的步骤。

首先，根据待纠正影像的轨道参数预测其概略的地理范围，进而估算出目标区域内控制点的大致地理位置范围。由于预测的概略位置有误差、加之待纠正遥感影像存在几何变形，一般在估算目标区域的地理位置范围时都会参考一个估计的误差半径R，其值一般为控制点影像大小的2至3倍。其中，假设该范围的左上角点经纬度坐标为（L1，B1），右下角点经纬度坐标为（L2，B2），待查询的目标区域影像中心的坐标为（L0，B0）；则L1=L0-R，L2=L0+R，B1=B0-R，B2=B0+R。

基于目标区域影像的中心经纬度的检索即是基于所估算的概略地理位置的检索。控制点影像数据库记录表中的数据项L和B，为基于概略位置的控制点影像检索提供了基本条件。由上述坐标关系可知，可建立“L1≤Li&&Li≤L2&&B2≤Bi&&Bi≤B1”的检索表达式进行关系检索。

基于目标影像属性检索，即根据目标影像的分辨率、传感器类型、成像时间等，筛选出可利用的控制点影像片，有效缩小检索结果范围。

根据本发明的实施例，首先估算待纠正遥感影像目标区域的概略地理位置范围，通过基于目标区域中心点经纬度的检索来缩小检索的范围；然后判断可用的控制点影像片的分辨率范围、时相及传感器类型等，进一步缩小检索结果的范围。具体的检索方式有：基于给定的坐标范围检索，基于传感器、时相及分辨率等属性信息综合检索，基于内容（需求分布特征、颜色特征、形状特征、纹理特征）的检索技术。按照这种检索流程对控制点影像片检索之后，满足条件的控制点影像数量会大大减少，基本上可以满足要求。

控制点对的匹配

遥感影像的纠正过程中最关键的步骤是同名控制点的选取，这是决定影像纠正自动化的关键。而根据已有的控制点影像数据进行自动匹配是实现影像的自动纠正基础。

控制点影像片的一个优势在于既有地理信息又有影像纹理信息，因此可以采用影像自动匹配算法找到控制点影像片对应的待纠正影像的像点位置，然后取出控制点影像片的地理坐标信息就可组成一个包括了物方坐标和像方坐标的控制点对。

根据本发明的遥感图像几何纠正方法，针对控制点影像片和待纠正影像在时相、分辨率、成像角度等方面都不一致的特点，采用了基于SIFT的几何不变特征提取与匹配算法，并采用粗糙模糊C-均值方法等方法对误匹配点进行自动剔除，最后采用经典的最小二乘匹配算法进行子像素级的精细化匹配。

具体流程如图4所示，在步骤S501，根据控制点的坐标信息、待纠正影像的元数据信息以及成像模型等计算对应的像点的初始坐标，然后按控制点影像片的大小从待纠正影像上切割出待搜索影像块；在步骤S503，利用Sift算法对控制点影像片以及切割后的待搜索影像块进行匹配，获得初步的匹配结果信息；在步骤S505，采用粗糙模糊C-均值方法和几何约束等方法进行误匹配点的剔除，保留可靠的准确匹配点对；在步骤S507，利用最小二乘算法对匹配结果进行精匹配，得到亚像素级的匹配精度；在步骤S509，将匹配成功的控制点对按规定的格式输出到包括控制点号、物方坐标、像方坐标的控制点信息文件中。

如图1所示，对于待纠正影像覆盖范围内的多个控制点影像片进行上述匹配操作，产生多个控制点对。如果控制点对数量不够，则通过继续检索获得在覆盖待纠正影像范围内的更多的控制点影像片，从而对控制点影像片的密度或数量进行增加，然后自动匹配的方式获得更多的控制点对，当控制点对的数量和分布满足几何纠正要求时，就可以进行自动化的几何纠正处理。下面给出根据本发明实施例的控制点影像片匹配相关的具体技术内容。

SIFT算法

SIFT算法基于图像特征尺度选择的思想，建立图像的多尺度空间，在不同尺度下检测到同一个特征点，确定特征点位置的同时确定其所在尺度，以达到尺度抗缩放的目的，此外，该算法剔除一些对比度较低的点以及边缘响应点，并提取旋转不变特征描述符以达到抗仿射变换的目的。该算法主要包含：（1）建立尺度空间，寻找候选点；（2）精确确定关键点位置，剔除不稳定点；（3）确定关键点梯度的模及方向；（4）提取特征描述符。

①尺度空间的生成

尺度空间理论最早出现于计算机视觉领域时其目的是模拟图像数据的多尺度特征。Koendetink在文献中证明高斯卷积核是实现尺度变换的唯一变换核，而Lindeberg等人则进一步证明高斯核是唯一的线性核。

二维高斯函数定义如下：

$G(x,y,\mathrm{\σ})=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)}/{2\mathrm{\σ}}^2$

σ代表了高斯正态分布的方差。

一幅二维图像，在不同尺度下的尺度空间表示可由图像与高斯核卷积得到：

L(x,y,σ)=G(x,y,σ)*I(x,y)

图5示出了不同尺度空间下的二维影像的一组例子。

为了有效的在尺度空间检测到稳定的关键点，提出了高斯差分尺度空间（DoGscale-space）。利用不同尺度的高斯差分核与图像卷积生成：

D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ)

DoG算子计算简单，是尺度归一化的LoG算子的近似。

然后构建图像金字塔，图像金字塔共O组，每组有S层，下一组的图像由上一组图像降采样得到。图6示出了高斯差分尺度空间（DoG）影像的一组例子。在图6中，对尺度空间octave，原始影像经过多次高斯卷积运算，产生一系列设定的尺度空间的影像。在右边的DoG影像是通过临近的高斯滤波后的影像进行差分运算产生的。在每一阶之后，高斯影像做因子为2的降采样，并重复进行该过程。

②空间极值点检测

为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如图7所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

③关键点位置确定与不稳定点剔除

a)关键点精确位置确定

利用尺度空间函数D(x,y,σ)的泰勒二次展开式进行最小二乘拟合，通过计算拟合曲面的极值来进一步确定关键点的精确位置和尺度。关键点最终的坐标和尺度可以精确到子像素级。

用泰勒公式展开D(x,y,σ)，则采样点原点为：

$D\left(X\right)=D+\frac{{\∂D}^T}{\∂X}X+\frac{1}{2}{X}^T\frac{{\∂}^{2}D}{{\∂X}^2}X$ （其中Χ=(x,y,σ)^T）

对X求导，并令其为零，即：，便可求得采样原点的位置

为： $\hat{X}=-\frac{{\∂}^2{D}^{-1}}{{\∂X}^2}\frac{\∂D}{\∂X}$

即为 $\left[\begin{array}{ccc}\frac{{\∂}^{2}D}{{\∂\mathrm{\σ}}^2}& \frac{{\∂}^{2}D}{\∂\mathrm{\σy}}& \frac{{\∂}^{2}D}{\∂\mathrm{\σx}}\\ \frac{{\∂}^{2}D}{\∂\mathrm{\σy}}& \frac{{\∂}^{2}D}{{\∂y}^2}& \frac{{\∂}^{2}D}{\∂\mathrm{yx}}\\ \frac{{\∂}^{2}D}{\∂\mathrm{\σx}}& \frac{{\∂}^{2}D}{\∂\mathrm{yx}}& \frac{{\∂}^{2}D}{{\∂x}^2}\end{array}\right]\left[\begin{array}{c}\mathrm{\σ}\\ y\\ x\end{array}\right]=-\left[\begin{array}{c}\frac{\∂D}{\∂\mathrm{\σ}}\\ \frac{\∂D}{\∂y}\\ \frac{\∂D}{\∂x}\end{array}\right]$

b)低对比度剔除

由，如果|D(X)|<0.03，则该点对比度较低，剔除。

c)边缘响应的去除

一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率通过一个2x2的Hessian矩阵H求出：

$H=\left[\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{xy}}& D_{\mathrm{yy}}\end{array}\right]$

导数由采样点相邻差估计得到。

D的主曲率和H的特征值成正比，令α为最大特征值，β为最小的特征值，则：

Tr(H)=D_xx+D_yy=α+β

Det(H)=D_xxD_yy-(D_xy)²=αβ

令α=γβ，则：

$\frac{{\mathrm{Tr}\left(H\right)}^2}{\mathrm{Det}\left(H\right)}=\frac{{(\mathrm{\&alpha;}+\mathrm{\&beta;})}^2}{\mathrm{\&alpha;\&beta;}}=\frac{{(\mathrm{r\&beta;}+\mathrm{\&beta;})}^2}{{\mathrm{r\&beta;}}^2}=\frac{{(r+1)}^2}{r}$

的值在两个特征值相等的时候最小，随着r的增大而增大。因此，为了检测主曲率是否在某域值r下，只需检测：

$\frac{{\mathrm{Tr}\left(H\right)}^2}{\mathrm{Det}\left(H\right)}<\frac{{(r+1)}^2}{r}$

取r＝10。

④关键点梯度模及方向计算

利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性。

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x,y+1)-L(x,y-1))}^2}$

$\mathrm{\&theta;}(x,y)=\mathrm{\&alpha;}\tan 2{\left((L(x,y+1)-L(x,y-1)\right)}^2+{(L(x+1,y)-L(x-1,y)))}^2$

上式为(x，y)处梯度的模值和方向公式。其中L所用的尺度为每个关键点各自所在的尺度。

⑤特征描述符生成

图8示出了由关键点邻域梯度信息生成特征向量的过程。图8左部分的中心点为当前关键点的位置。首先，将坐标轴旋转为关键点的方向，以确保旋转不变性。接下来以关键点为中心取8×8的窗口。在图8中，每个小格代表关键点邻域所在尺度空间的一个像素，箭头方向代表该像素的梯度方向，箭头长度代表梯度模值，图中圆圈代表高斯加权的范围（越靠近关键点的像素梯度方向信息贡献越大）。然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图8右部分所示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。

在计算过程中，为了增强匹配的稳健性，可选地，对每个关键点使用4×4共16个种子点来描述，这样对于一个关键点就可以产生128个数据，即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。

图9示出了不同影像由关键点邻域梯度信息生成特征向量图的一组例子。

基于SIFT描述子的遥感影像匹配

在对单独的一幅图像进行上述的特征提取和特征描述以后，就得到了该图中所有的特征及其描述子，设为图像1（即实时图像），其特征点数量为m。要做到将两幅图像匹配起来，也就是得到两幅图像中匹配的像素点（这里指特征点）。首先，要对另外一幅图像（图像2，即参考图像）进行完全相同的特征提取和特征描述过程，得到数量为n的特征；其次就是要在图像1中m个特征点和图像2中n个特征点中搜索出r对匹配的特征点，其中r≤m，r≤n，而且为了能够鲁棒地计算出图像1和图像2之间的几何关系，应该保证r≥8，如果不能寻找到满足条件的r，就需要调整匹配的精度，比如降低匹配特征点之间的相似性要求，使得能获得更多的匹配特征点对。

特征点的匹配实质上是其描述子的匹配，特征点的描述子其实是对该特征进行了一个定量的描述，以能够应用于匹配算法。特征点描述子的匹配其实是在描述空间中进行的，比如SIFT特征描述子为128维的向量，因此SIFT描述子是在128维空间中进行匹配的。在描述空间中，特征点描述子的匹配程度则是以距离来测定的，距离最近的两个描述子一般就代表了匹配的一对特征点。而在描述空间中，一般有以下两种距离的定义：

a)欧氏（Euclidean）距离。即对于p维空间中的两点x和y（均为p维向量），有它们的欧氏距离定义为：

$d_E(x,y)=\sqrt{{(x_1-y_1)}^2+...+{(x_p-y_p)}^2}=\sqrt{{(x-y)}^T(x-y)}$

式中x=(x₁,…,x_p)^T,y=(y₁,…,y_p)^T。欧氏距离可应用于基于直方图的描述子，即描述子的每一维都有相同的权重，比如SIFT描述子，GLOH描述子和PCA-SIFT描述子等。

b)马氏（Mahalanobis）距离。若描述子的每一维都有不同的权重，则需要使用马氏距离测量它们之间的距离，此时设描述子的权重向量为s=(s₁,…,s_p)^T，则P维空间中的两点x和y之间的马氏距离为：

$d_M(x,y)=\sqrt{{\left(\frac{x_1-y_1}{s_1}\right)}^2+...+{\left(\frac{x_p-y_p}{s_p}\right)}^2}=\sqrt{{(x-y)}^T{D}^{-1}(x-y)}$

式中马氏距离可以应用于方向可控滤波描述子，差分不变描述子，矩不变描述子和复数滤波描述子等。

在根据本发明的实施例中，使用SIFT特征描述子进行匹配，其将欧氏距离作为评判描述子匹配程度的标准。图10示出了根据本发明一个实施例的基于Sift描述子的匹配算法流程。图10中方框内的部分表示SIFT特征点的提取和特征描述过程。

根据图10所示的算法，在输出匹配结果之前还要通过“最近距离与第二近距离比较选择法”来选择匹配点。具体而言，首先设定一个合理的阈值t，然后比较A与n个点中所有点的距离，然后找出找出与A距离最小的点B以及与A距离第二小的点C，并且若d_ABd_AC<t时，才认为A与B是合理正确的匹配点。这样做的好处是当A有很多相似的匹配时，即n个点中与A最小的距离和与A第二小的距离相差不是很多时，会认为这不是一个合理的匹配。只有当最小的距离与第二小的距离相差很多时，即距离最小的点“遥遥领先”其它点时，才认同它，因为这基本保证了这是一个非常稳定的正确匹配而不是一个模棱两可的匹配。

误匹配点剔除

在得到了实时图像和参考图像的匹配特征点对后，基本上已经实现了图像匹配的任务，但是匹配的特征点只能代表图像对之间的局部关系，有限数量的匹配特征点并不能完全反映图像对之间的全局关系。在对应的匹配点对中，可能存在误匹配或是匹配误差较大的点对，由于这些点对的存在将影响匹配精度，为基于匹配点文件的应用造成影响，所以在初始匹配结束后，探索误匹配点对的剔除方法也非常重要。进一步，误匹配点自动剔除能有效自动化地提高匹配的准确率。

图11示出了根据本发明一个实施例的包括了误匹配点剔除步骤的匹配算法流程。在图11中，首先根据前述SIFT算法基于几何不变特征进行特征点匹配，然后对于匹配结果进行筛选删除，其中涉及到粗糙模糊C-均值方法，随机抽样一致性检测方法（RANSAC方法），以及多项式拟合的最小二乘方法。下面分别描述这些方法步骤。

如前所述，利用SIFT算法分别提取实时影像与参考影像的特征向量信息。根据两幅图像的SIFT特征向量，采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量，具体原则是最近邻(NN)和第二近邻(SCN)的距离之比(NN/SCN)最小作为两幅图像中关键点的相似性判定度量。之后，利用粗糙模糊C-均值方法进行误匹配点的初步剔除。

误匹配点剔除的粗糙模糊C-均值方法

聚类分析的基本思想非常朴素、直观和简单，它是根据各个待分类的模式特征相似程度进行分类的，相似的归为一类，不相似的作为另外一类。聚类分析包括两个基本内容：模式相似性的度量和聚类算法。模糊C-均值（FuzzyC-Means）算法是一种以最小类内平方误差和为聚类准则，计算每个样本属于各模糊子集（聚类）的隶属度，通过目标函数极小化的必要条件之间的Pickard迭代来实现算法。

定义1：设X={x₁,x₂,x₃,…，x_n}为待分类对象集，对于第i类w_i，其质心为v_i，定义 为可以肯定属于w_i类的对象集合，则：

a)若，则对于 $\&ForAll;k\&Element;\{\mathrm{1,2},...,n\},k\&NotEqual;j,x_j\&NotElement;{\stackrel{\&OverBar;}{R}w}_i;$

b)同时为了保证A_i不至于取得过小，需要满足对于，则至少存在k∈{1,2,…,n}，使得

其中A_i成为上近似限，上近似限刻画了所有可能属于第i类的对象的边界，若某个对象不属于上近似限所界定的范围，则它属于这个类的负域，即完全不属于这个类。

定义2：粗糙模糊C-均值(RoughFuzzyC-Means,RFCM)算法的目标函数为：

$J_m(U,V)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}}_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}{u}_{\mathrm{ij}}^m{d}_{\mathrm{ij}}^2$

约束条件为：

●u_ij∈[0,1]

● $0<\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}<N$

● ${\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}}_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}{u}_{\mathrm{ij}}^m=1$

利用拉格朗日乘数法，可得到无约束的准则函数：

$F=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}{d}_{\mathrm{ij}}^2-\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i\left(\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}-1_{x_j\&Element;\stackrel{\&OverBar;}{R}w}\right)$

上式的极值条件为：

$\frac{\&PartialD;F}{{\&PartialD;u}_{\mathrm{ij}}}=0,\frac{\&PartialD;F}{{\&PartialD;\mathrm{\&lambda;}}_j}=0$

对上式进行计算可得：

$u=1/\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{(d_{\mathrm{ij}}^2/d_{\mathrm{kj}}^2)}^{\frac{1}{m-1}}_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}$

质心计算公式不变公式为：

$v_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{x}_j/\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m$

从RFCM算法可以得到以下两个性质：

● $x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i\&DoubleLeftRightArrow;u_{\mathrm{ij}}=1;$

● $u_{\mathrm{ij}}^{\mathrm{FCM}}\&le;u_{\mathrm{ij}}^{\mathrm{RFCM}}$

RFCM算法的主要思想是把属于某个类的对象分成了肯定的、可能的和否定的三个集合，以所有可能的对象的最小类内平方误差和为聚类准则进行聚类。RFCM算法和FCM算法最大的不同在于，它认为x_j属于w_i的隶属度u_ij的计算只与上近似中包含x_j的类有关，若某个类w_i的上近似中不包含x_j，则这个类对x_j的隶属度是没有任何贡献的。

RFCM算法的目标函数是寻找各个聚类最小类内距离平方和，它排除了超球体以外的不可能属于类的对象。对于u_ij的计算，若对于类w₁，对象x_j仅仅位于类的上近似，内，则，即u_1j=1，若x_j属于和的交集，则u_1j，u_2j的计算仅仅与类w₁，和w₂有关，而与w₃无关，即：

$u_{1j}=1/{(1+\frac{{\left|\right|x_j-v_1\left|\right|}^2}{{\left|\right|x_j-v_2\left|\right|}^2})}^{\frac{1}{m-1}},$ $u_{1j}=1/{(1+\frac{{\left|\right|x_j-v_2\left|\right|}^2}{{\left|\right|x_j-v_1\left|\right|}^2})}^{\frac{1}{m-1}}$

若x_j不在内，则u_1j=0。

目标函数，给出了u_ij的更新公式。根据此公式，以v_j为中心半径为α的超球体内的u_ij为1。

基于上述数理分析，根据本发明实施例的几何约束下的粗糙模糊C-均值方法误匹配点剔除步骤主要包括：

根据生成的匹配点对及影像的分辨率信息，计算影像间的几何约束条件：

●同分辨率情况下，求对应点对的斜率k及距离S；

●不同分辨率情况下，求对应点对连线延长线的交点坐标（X，Y）。

然后，利用粗糙模糊C-均值方法(RFCM)，求（k,S）或（X,Y）的隶属度，对所有的匹配点对进行聚类分析，删除不同类的点对，只保留包含在肯定聚集集合内的点。

i.确定类数c(2≤c≤N)，参数m，初始矩阵，类的上近似边界A_i和一个适当小数，s=0；

ii.计算质心

iii.若则u_ij=0，否则更新

iv.若||U^s-U^s+1||<ε，则停止，否则，s=s+1，转到b。

如图11所示，如果经过粗糙模糊C-均值方法进行匹配点筛查后仍有超过设定数目（例如，7）的匹配点，则采用随机抽样一致性检测方法进行进一步的匹配点筛除。

随机抽样一致性检测方法（RANSAC方法）

a)RANSAC算法的计算量

RANSAC算法中，要求保证在一定的置信概率下，M组抽样中至少有一组抽样的数据全是内点（inliers）。利用下式可以求得满足要求的最小抽样数M。

1-(1-(1-ε)^m)^M=p

其中:ε为数据错误率（外点（outliers）在原始数据所占的比例），m为计算模型参数需要的最小数据量，P为置信概率。从上式可以看出，M和ε，m，P呈指数关系。下表示出了P=0.95时，M随m,ε变化的情况。

从上表可以看出，当模型比较复杂、ε较高时，M很大，直接造成RANSAC算法效率下降。设从原始数据中随机抽取一组抽样需要时间；由一组抽样计算模型参数需要时间t；用一个数据检验模型参数需时间平均为t，则用N个数据检验（全数据检验）需要时间N_t。因此RANSAC算法所需计算时间为：

T=M(t_s+t_c)+MN_t

其中，M(t_s+t_c)为M组抽样抽取及模型参数计算需要的时间，MN_t为M个模型参数检验需要的时间。

从该式看出，RANSAC算法需要的时间由两部分组成：

●M组抽样选择和模型参数估计需要的时间；

●模型参数检验需要的时间。在模型确定、数据错误率确定的情况下，

为了保证结果的置信概率，M是不能减少的。

因此，为了提高算法的效率，只能从减少参与检验的模型参数数量、减少模型参数检验需要的时间出发。

b)RANSAC算法步骤

根据置信概率P和数据错误率ε计算最小抽样数M；计算抽样对应的模型参数，用所有原始数据检验模型参数质量，获得每个模型参数的inliers数量；根据inliers数量和误差的方差来选择最优模型参数；用最优模型参数对应的inliers估计最终模型参数。

最小二乘法精匹配

经过RANSAC算法进行更进一步的匹配点选择后，还可以利用多项式拟合的最小二乘方法剔除拟合残差大于倍中误差的匹配点对。

如图11所示，如果匹配点数较小，例如不大于7，可选地，可以采用参数自适应SIFT算法来增加供筛选的匹配点。

卫星遥感影像由于受到天气、阳光、遮挡等外界因素的严重影响，并且存在因不同的成像时间、角度、距离等因素而导致的图像平移、旋转、缩放等问题以及各种传感器之间的成像模式相机参数轨道模型的不同，即使采用SIFT算法，在极端情况下也可能出现提取的特征点少或根本提取不到特征点的现象，造成误匹配或匹配失败。针对这种情况，可以采用改进的参数自适应SIFT影像匹配方法，根据不同的影像特点及质量，选取相应的策略并根据相关的权重定义规则确定特征点探测的相应阈值。

在SIFT特征提取及匹配方法中涉及的参数较多，为了验证各参数在特征提取及匹配过程中起到的作用，进行了匹配实验，主要通过修改某个待验证的参数的大小，即对该参数按照算法中理论变化规则进行调节，与此同时保持其他参数不变，通过匹配计算统计各参数大小时的匹配点个数、误匹配点个数、匹配时间等，并进行统计分析。

通过对实验的统计结果分析可知，其一般规律如下：

●当每阶（Octave）的采样间隔逐渐增大时，匹配点个数增加迅速，误匹配点数量没有改变，匹配时间增加；

●逐渐增加高斯卷积核σ，匹配点数整体成下降趋势，但当σ=1.6时结果最好；

●对比度阈值逐渐增大，匹配点个数逐渐减少，误匹配点数据逐渐减少；

●随着曲比率阈值的逐渐增大，匹配点个数逐渐增多，误匹配点数几乎不变；

●匹配门限阈值最近距离与次最近距离的比值增大，匹配点个数逐渐增加，误匹配点数亦逐步增加。

基于上述实验结果，提出了参数自适应SIFT算法。参数自适应SIFT算法的步骤是通过影像自身的信息计算自动确定各参数的阈值，从而通过改变参数起到增加提取的特征点数量、增加匹配点数量的目的。具体可以包括如下步骤：

a)计算输入影像的灰度平均值或对影像做自动色阶，然后将平均灰度除以10如果结果小于1.0，则减小对比度的阈值，增加提取的特征点数量；

b)在满足匹配精度的情况下，如果匹配点数小于7时，则增大曲比率阈值；

c)若经过a）b)两步提取的特征数量仍较少或是匹配数量较少时，在不影响匹配精度的情况下，适当放大最近距离与次最近距离的比例阈值。

如图1所示，在进行了控制点匹配后，根据传感器类型和成像特点，选择相应的几何纠正模型（步骤9），然后基于所匹配的控制点对和所选择的几何纠正模型建立像素坐标与大地坐标之间的转换关系（步骤11）,最后基于数字高程模型（DEM）和像素坐标与大地坐标之间的转换关系对原始影像进行几何纠正（步骤13）。下面主要描述几何纠正模型的选择和几何纠正的步骤。

几何纠正模型及其选择

线阵推扫式成像传感器是逐行以时序方式获取二维图像的。一般是先在像面上形成一条线图像，然后卫星沿着预定的轨道向前推进，逐条扫描后形成一幅二维影像，成像方式如图12所示。影像上每一行像元在同一时刻成像且为中心投影，整个影像为多中心投影。

在图12中：p_k为影像上任一像点，x_k为扫描线k上影像点的x坐标，c为传感器主距，O_k为扫描线k的投影中心，o_k为扫描线k的主点，l_k为扫描线k从投影中心O_k发出的光线。

几何纠正中一个最基本的问题就是要建立合理的遥感影像成像模型，所谓遥感影像的成像模型是指建立影像上的坐标(x,y)与其对应的地面点的大地坐标(X,Y,Z)之间的数学关系。即：

X=f_x(x,y,g)

Y=f_y(x,y,h)

其中g，h为其他因素的影响。

遥感影像的成像模型可以分为两大类：物理模型和通用模型。

物理模型是指考虑成像时造成影像变形的物理意义如地表起伏、地球曲率、大气折射、相机畸变等因素，然后利用这些物理条件建构成像几何模型，最为代表性的是共线条件方程为基础的严格成像模型。

通用模型不考虑成像机理，而是直接用数学函数来描述像点与物点的几何关系，它具有一般性、保密性、高效性的特点。常用的通用成像模型有多项式、直接线性变换、仿射变换、有理函数模型等；

在根据本发明的实施例中，以严格成像模型、有理函数模型和多项式模型为几何纠正的可选模型。根据传感器类型和成像特点，选择相应的几何纠正模型，例如，SPOT5卫星影像，一般选择严格成像模型，而对于IKONOS、P5、RapidEye等卫星影像则选择有理函数模型，无资轨参数或有理函数模型系数文件的，则选择多项式模型。

严格成像模型描述的是成像时像点与物方点一一对应的严格几何过程，以常用的共线条件方程为代表。根据本发明的一个实施例，依据SPOT5提供的高精度的卫星星历数据、姿态数据以及时间数据，通过成像过程，建立像点(l,p)与地面点(X,Y,Z)或之间的数学关系。图13示出了像点与地面点的数学关系。

严格成像模型所用到的坐标系可能包括以下坐标系，如图14所示：

导航坐标系（也称本体坐标系）O₁-X₁Y₁Z₁：是由卫星的姿态控制***（ACS）决定的，其中原点在卫星质心，X₁轴、Y₁轴、Z₁轴分别取卫星的三个主惯量轴。X₁轴沿着卫星横轴，Y₁轴沿着纵轴指向卫星飞行方向（不完全重合），Z₁轴按照右手系规则确定。卫星姿态测量在本体坐标系中进行，描述其空间姿态的三个参数是：俯仰（pitch）为绕本体坐标系X₁轴的旋转，滚动（roll）为绕本体坐标系Y₁轴的旋转，航偏（yaw）为绕本体坐标系Z₁轴的旋转。

(ψ_x)_p和(ψ_y)_p分别代表电荷耦合器件（CCD）第P个元器件成像时沿轨道方向和垂直轨道方向的倾角。

轨道坐标系：是导航坐标系与地面坐标系进行转化的一个过渡坐标系，轨道坐标系的原点在卫星质心，Z₂轴指向地心反向，Y₂轴在卫星轨道面上指向卫星运动的方向，X₂轴按照右手系规则确定。

O₂-X₂Y₂Z₂各个轴的方向由卫星在t时刻的位置和速度两个矢量方向确定，如下：

${\stackrel{\&RightArrow;}{Z}}_2=\frac{\stackrel{\&RightArrow;}{P}\left(t\right)}{\left|\right|\stackrel{\&RightArrow;}{P}\left(t\right)\left|\right|}$

${\stackrel{\&RightArrow;}{X}}_2=\frac{\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2}{\left|\right|\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2\left|\right|}$

${\stackrel{\&RightArrow;}{Y}}_2={\stackrel{\&RightArrow;}{Z}}_2\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{X}}_2$

地面坐标系：SPOT卫星的位置和速度都是定义在地球协议参考坐标系（ITRS），ITRF是ITRS的实现，是由IERS中心局的地球参考框架部负责建立和维护的，该参考框架的站坐标和速度是利用诸如甚长基线干涉测量（VLBI）、激光测卫等空间大地测量技术所采集的数据计算出来的，是一个不断变化的参考框架，目前采用的是ITRF90，为了方便计算，一般采用WGS84坐标系来代替ITRF框架确定的坐标***。这是因为，WGS84为地心地固坐标系，ITRF则为惯性坐标系，最终测图需要的是地固坐标系，且ITRS与WGS84相差仅在1米左右，从ITRF90到WGS84的转换参数如下表所示：

T1(m)	T2(m)	T3(m)	D(ppm)	R1(arcsec.)	R2(arcsec.)	R3(arcsec.)
							0.060	-0.517	-0.223	-0.011	0.0183	-0.0003	0.0070

ITRF90到WGS84的转换参数表

转化公式如下：

$\left(\begin{array}{c}{X}_{\mathrm{WGS}84}\\ Y_{\mathrm{WGS}84}\\ Z_{\mathrm{WGS}84}\end{array}\right)=\left(\begin{array}{c}{T}_1\\ T_2\\ T_3\end{array}\right)+(1+D)\left(\begin{array}{ccc}0& {-R}_3& R_2\\ R_3& 0& -R_1\\ -R_2& R_1& 0\end{array}\right)\left(\begin{array}{c}{X}_{\mathrm{ITRF}}\\ Y_{\mathrm{ITRF}}\\ Z_{\mathrm{ITRF}}\end{array}\right)$

WGS84坐标系是一个由全球地心参考坐标框架和一组相应的模型以及WGS84大地水准面所组成的测量参照系。该坐标系的原点位于地球质心，Z轴与IERS参考极指向相同，X轴指向IERS参考子午线与过原点且垂直于Z轴平面的交点，Y轴由右手系规则确定；

成像过程推导(l,p)→(X,Y,Z)如下：

ⅰ、计算l行p列像素所对应的成像时间

t=t_c+lsp×(l-l_c)

其中，t_c为中央扫描行的摄影时间，lsp为每行的扫描时间，l_c为中央扫描行的行号

ⅱ、内插出t时刻卫星的位置

$\stackrel{\&RightArrow;}{P}\left(t\right)=\underset{j=1}{\overset{8}{\mathrm{\&Sigma;}}}\frac{\stackrel{\&RightArrow;}{P}\left(t_j\right)\&times;\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}{\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}$

$\stackrel{\&RightArrow;}{V}\left(t\right)=\underset{j=1}{\overset{8}{\mathrm{\&Sigma;}}}\frac{\stackrel{\&RightArrow;}{V}\left(t_j\right)\&times;\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}{\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}$

ⅲ、通过p列的成像倾角(ψ_x)_p和(ψ_y)_p计算(l,p)在导航坐标系下的成像方向矢量 ${\stackrel{\&RightArrow;}{u}}_1^{\&prime;}=\left(\begin{array}{c}-\tan \left({\mathrm{\&psi;}}_Y\right)\\ +\tan \left({\mathrm{\&psi;}}_X\right)\\ -1\end{array}\right)$

${\stackrel{\&RightArrow;}{u}}_1=\frac{{\stackrel{\&RightArrow;}{u}}_1^{\&prime;}}{\left|\right|{\stackrel{\&RightArrow;}{u}}_1^{\&prime;}\left|\right|}$

ⅳ、计算轨道坐标系下的成像方向矢量

${\stackrel{\&RightArrow;}{Z}}_2=\frac{\stackrel{\&RightArrow;}{P}\left(t\right)}{\left|\right|\stackrel{\&RightArrow;}{P}\left(t\right)\left|\right|}$

${\stackrel{\&RightArrow;}{X}}_2=\frac{\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2}{\left|\right|\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2\left|\right|}$

${\stackrel{\&RightArrow;}{Y}}_2={\stackrel{\&RightArrow;}{Z}}_2\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{X}}_2$

其中，和分别是像素(l,p)成像时刻t卫星的位置和速度矢量内插出t时刻卫星的姿态角：

$a_p\left(t\right)=a_p\left(t_i\right)+(a_p\left(t_{i+1}\right)-a_p\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

$a_r\left(t\right)=a_r\left(t_i\right)+(a_r\left(t_{i+1}\right)-a_r\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

$a_y\left(t\right)=a_y\left(t_i\right)+(a_y\left(t_{i+1}\right)-a_y\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

a_p(t_i)，a_r(t_i)，a_y(t_i)分别为t_i时刻卫星的姿态角，可以从辅助数据文件中获取。

将导航坐标系下的成像方向矢量转化为轨道坐标系下的成像方向矢量

${\stackrel{\&RightArrow;}{u}}_2=\frac{{\stackrel{\&RightArrow;}{u}}_2^{\&prime;}}{\left|\right|{\stackrel{\&RightArrow;}{u}}_2^{\&prime;}\left|\right|}$

${\stackrel{\&RightArrow;}{u}}_2^{\&prime;}=M_p\&CenterDot;M_r\&CenterDot;M_y\&CenterDot;{\stackrel{\&RightArrow;}{u}}_1$

$M_p=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(a_p\left(t\right)\right)& \sin \left(a_p\left(t\right)\right)\\ 0& -\sin \left(a_p\left(t\right)\right)& \cos \left(a_p\left(t\right)\right)\end{array}\right]$

$M_r=\left[\begin{array}{ccc}\cos \left(a_r\left(t\right)\right)& 0& -\sin \left(a_r\left(t\right)\right)\\ 0& 1& 0\\ \sin \left(a_r\left(t\right)\right)& 0& \cos \left(a_r\left(t\right)\right)\end{array}\right]$

$M_y=\left[\begin{array}{ccc}\cos \left(a_y\left(t\right)\right)& -\sin \left(a_y\left(t\right)\right)& 0\\ \sin \left(a_y\left(t\right)\right)& \cos \left(a_y\left(t\right)\right)& 0\\ 0& 0& 1\end{array}\right]$

ⅴ、计算地面坐标系下的成像方向矢量

${\stackrel{\&RightArrow;}{u}}_3=\left[\begin{array}{ccc}{\left(X_2\right)}_X& {\left(Y_2\right)}_X& {\left(Z_2\right)}_X\\ {\left(X_2\right)}_Y& {\left(Y_2\right)}_Y& {\left(Z_2\right)}_Y\\ {\left(X_2\right)}_Z& {\left(Y_2\right)}_Z& {\left(Z_2\right)}_Z\end{array}\right]\&CenterDot;{\stackrel{\&RightArrow;}{u}}_2$

其中，X₂,Y₂,Z₂为轨道坐标系O-X₂Y₂Z₂的三个坐标轴方向矢量ⅵ、地面点的定位：

如图15所示，观测矢量u₃与地面的交点为M，卫星的位置矢量为M与地心连接对应的矢量为。又因为M点在长轴为a+h,短轴为b+h的椭球面上，所以联立方程后有：

$\stackrel{\&RightArrow;}{O_{3}M}=\stackrel{\&RightArrow;}{P}\left(t\right)+\mathrm{\&mu;}\&times;u_3\&DoubleRightArrow;\left\{\begin{array}{c}X=X_P+\mathrm{\&mu;}\&times;{\left(u_3\right)}_X\\ Y=Y_P+\mathrm{\&mu;}\&times;{\left(u_3\right)}_Y\\ Z=Z_P+\mathrm{\&mu;}\&times;{\left(u_3\right)}_Z\end{array}\right.$

$\frac{X^2+Y^2}{A^2}+\frac{Z^2}{B^2}=1$

其中，

$\left\{\begin{array}{c}A=a+h\\ B=b+h\end{array}\right.$

结合前面的几个变换关系，写成类似共线条件方程的形式：

$u_1=\left[\begin{array}{c}-\tan {\left({\mathrm{\&psi;}}_y\right)}_p\\ +\tan {\left({\mathrm{\&psi;}}_x\right)}_p\\ -1\end{array}\right]=\frac{1}{u}{R}_{21}^{-1}{.R}_{32}^{-1}.\left[\begin{array}{c}X-X_P\\ Y-Y_P\\ Z-Z_P\end{array}\right]$

获得导航坐标系下的成像方向矢量与地面坐标系的观测方向的关系；投影中心为(X_P,Y_P,Z_P)，坐标矢量为(X,Y,Z)，投影系数为u，旋转矩阵为:

$R=R_{21}^{-1}[a_r\left(t\right),a_p\left(t\right),a_y\left(t\right)]R_{32}^{-1}[P\left(t\right),V\left(t\right)]=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]$

可得共线方程：

$\frac{r_{11}(X-X_P)+r_{12}(Y-Y_p)+r_{13}(Z-Z_p)}{r_{31}(X-X_P)+r_{32}(Y-Y_p)+r_{33}(Z-Z_p)}-\tan {\left({\mathrm{\&psi;}}_y\right)}_p=0$

$\frac{r_{21}(X-X_P)+r_{22}(Y-Y_p)+r_{23}(Z-Z_p)}{r_{31}(X-X_P)+r_{32}(Y-Y_p)+r_{33}(Z-Z_p)}+\tan {\left({\mathrm{\&psi;}}_x\right)}_p=0$

由此可以得到遥感影像成像过程中影像上的坐标(x,y)与其对应的地面点的大地坐标(X,Y,Z)之间的数学映射关系。

RFM模型（RationalFunctionModel，有理函数模型）是通用传感器模型中的一种，它的实质是通过更复杂的数学模型来描述严格成像过程，建立物方空间坐标与像方像点坐标的直接对应关系。RFM模型可以有效对卫星轨道参数、传感器属性等涉密参数进行封装，用户可以在不知道相关信息的前提下，直接根据RFM模型及参数就可解算出像点的空间坐标，从而进行各种应用。RFM模型一般采用RPC参数（RationalPolynomialCoefficients，有理多项式系数）来表示，RPC参数一般有80个。

a、利用RPC定位的RFM模型原理

RFM模型是将像点坐标(s,l)表示为以相应地面点空间坐标(φ,λ,h)(其中φ为经度，λ为纬度，h为高程)为自变量的多项式的比值：

$s=\frac{\mathrm{Nums}(U,V,W)}{\mathrm{Dens}(U,V,W)}$

$l=\frac{\mathrm{Numl}(U,V,W)}{\mathrm{Denl}(U,V,W)}$

其中，

$U=(\mathrm{\&phi;}-{\mathrm{\&phi;}}_0)/{\mathrm{\&phi;}}_s$ $V=(\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_0)/{\mathrm{\&lambda;}}_s$

$W(h-h_0)/h_s$

$s=(S-S_0)/S_s$

$l=(L-L_0)/L_s$

(U,V,W)是标准化后的地面点空间坐标(φ,λ,h)，(s,l)是标准化后的像点坐标(S,L)，a_i,b_i,c_i,d_i(i=1～20),φ₀,φ_s,λ₀,λ_s,h₀,h_s,S₀,S_s,L₀,L_s为相片的RPC参数,φ₀,λ₀,h₀，S₀，L₀为标准化平移参数，φ_s,λ_s,h_s,S_sL_s为标准化比例参数。它们与影像数据一起被提供。

Nums(U,V,W)=a₁+a₂V+a₃U+a₄W+a₅VU+a₆VW+a₇UW+a₈V²

+a₉U²+a₁₀W²+a₁₁VUW+a₁₂V³+a₁₃VU²+a₁₄VW²

+a₁₅V²U+a₁₆U³+a₁₇UW²+a₁₈V²W+a₁₉U²W+a₂₀W³

Dens(U,V,W)、Numl(U,V,W)、Denl(U,V,W)也同样可以表示为关于(U,V,W)的三次多项式的形式。

b、RPC参数的求解

将RFM模型变形为：

F_X=Num_s(P,L,H)-X*Den_s(P,L,H)=0

F_Y=Num_L(P,L,H)-Y*Den_L(P,L,H)=0

则误差方程为：

V=Bx-l,W

上式中，

$B=\left[\begin{array}{cccc}\frac{{\&PartialD;F}_X}{{\&PartialD;a}_i}& \frac{{\&PartialD;F}_X}{{\&PartialD;b}_j}& \frac{{\&PartialD;F}_X}{{\&PartialD;c}_i}& \frac{{\&PartialD;F}_X}{{\&PartialD;d}_j}\\ \frac{{\&PartialD;F}_Y}{{\&PartialD;a}_i}& \frac{{\&PartialD;F}_Y}{{\&PartialD;b}_j}& \frac{{\&PartialD;F}_Y}{{\&PartialD;c}_i}& \frac{{\&PartialD;F}_Y}{{\&PartialD;d}_j}\end{array}\right],$ (i=1…20,j=1…20)

$1=\left[\begin{array}{c}-F_X^0\\ -F_Y^0\end{array}\right]$

$x={\left[\begin{array}{cccc}{a}_i& b_j& c_i& d_j\end{array}\right]}^T$

W为权矩阵。

根据最小二乘平差原理，可将误差方程变换为法方程：

(B^TB)x=B^Tl

当用于解算RPC模型参数的采用二阶或者二阶以上时，RPC模型存在过度参数化，RPC模型中分母的变化非常剧烈，这样就导致设计矩阵(B^TB)的状态变差，设计矩阵变为奇异矩阵，导致最小二乘平差不能收敛。

采用谱修正迭代法求解法方程的方法，不论法方程呈良态、病态或秩亏，其解算程序均不需加任何当法方程呈良态时，经几次迭代就可收敛到精确解。当法方程呈病态时，收敛速度稍慢，且估计结果无偏。具体公式表述为：

设有(B^TB)x=B^Tl，将上面式子两边同时加上x，得到

(B^TB+E)x=B^Tl+x

其中，E为t阶单位矩阵，由于式子两边都含有未知参数x，所以只能采用迭代的方法求解，其迭代公式为：

x^(k)=(B^TB+E)^-1(B^Tl+x^(k-1))

多项式模型是回避成像的空间几何过程，直接采用一个次数适当的数学多项式来描述纠正前后图像相应点之间的坐标关系。

一般多项式纠正变换公式为：

u=a₀₀+a₁₀X+a₀₁Y+a₂₀X²+a₁₁XY+a₀₂Y²+a₃₀X³+a₂₁X²Y+a₁₂XY²+a₀₃Y³+...

v=b₀₀+b₁₀X+b₀₁Y+b₂₀X²+b₁₁XY+b₀₂Y²+b₃₀X³+b₂₁X²Y+b₁₂XY²+b₀₃Y³+...

式中，（u，v）为像点的像平面坐标；（X，Y）为其对应地面点的大地坐标；a_ij,b_ij为多项式的系数。这里多项式的阶数一般不大于三次，因为更高阶的多项式往往不能提高精度反而会引起参数的相关，造成模型定向精度的降低。

在基于匹配的控制点对和选择的几何纠正模型建立像素坐标与大地坐标的转换关系后，在步骤S13，基于数字高程模型和像素坐标与大地坐标之间的转换关系进行几何纠正，得到经过正射纠正的数字正射影像。这里的高程模型是从已有的基础测绘成果图中得到的。

另外，可选地，在步骤S13之后，进行纠正结果精度检查，判断纠正结果是否满足要求，如果不满足要求，则对控制点的分布和数量选择结果进行人工干预调整。

根据本发明的实施例，在实现上述基于控制点影像数据库的卫星遥感图像快速几何纠正方法时，可以采用基于OpenMP的几何纠正并行加速处理方法。

基于OpenMP的几何纠正并行加速处理

OpenMP的规范由SGI发起，它是一种面向共享内存以及分布式共享内存的多处理器多线程并行编程语言。OpenMP具有良好的可移植性，支持Fortran和C/C++编程语言，操作***平台方面则支持UNIX***以及Windows***。OpenMP的重要性在于，它能够为编写多线程程序提供一种简单的方法，无需程序员进行复杂的线程创建、同步、负载平衡和销毁工作。OpenMP对于For循环语句特别适用，它能过通过在原有程序中添加几句简单的语句就实现多核的并行处理，充分利用电脑的CPU资源，实现计算的加速。因此，针对分块的几何纠正采用OpenMP进行单机多核并行处理将非常有效，而且对硬件也不会有过多要求。

遥感影像在纠正过程实际上是从一个几何空间转换到另一几何空间，一旦建立了纠正模型，则这种转换是唯一的，只与位置有关，另外由于成像是连续的，纠正后影像也是连续的，因此可以通过分块的思想，通过计算四个角点转换前后的坐标，然后建立块内相对简单的仿射变换模型，从而简化纠正模型的计算，减少计算量。考虑到地形的影响，分块的大小是需要考虑的，理论上应该是地形起伏大的地方分块小，起伏小的地方分块大，而诸如SPOT5、P5等分辨率都在2.5米左右，所以最终分块大小定为15×15，而对于HJ_1A\1B分块大小则定为5×5。在此基础上，为了考虑对大图像处理的需要，可以对原始图像进行分块操作，即将原图像分为若干块，每一块看作一幅图像单独进行纠正然后写入磁盘，最后合成。

下表示出了针对不同影像类型进行基于OpenMP的几何纠正处理的例子。

从上表可以看出，对于严格成像模型，采用分块加并行的处理能显著提高计算效率，加速比在50左右，而对于有理函数模型和多项式其加速比也在4左右。

根据本发明实施例的遥感图像纠正方法，将影像自动匹配技术引入到控制点的自动选取中，可以实现遥感影像的自动或半自动几何纠正处理。具体而言，该方法基于控制点影像数据库，通过对于控制点影像库的查询检索、控制点影像自动匹配获取控制点匹配数据；再利用严格成像模型、有理函数模型、多项式模型等多种模型的几何纠正；实现了完整的遥感影像自动几何纠正解决方案，有效地利用了已有的控制点资源，提高了生产效率。同时，该方法可以通过单机多核并行快速计算，从而进一步提高遥感影像的纠正效率。

以上所述仅是本发明的示范性实施方式，而非用于限制本发明的保护范围，本发明的保护范围由所附的权利要求确定。

Claims (8)

1.一种基于控制点影像数据库的卫星遥感图像快速几何纠正方法，其特征在于，包括：

确定待纠正影像的地理范围；

针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索所有符合要求的控制点；

进行自动匹配选点，确定用于几何纠正的控制点对；

判断所述匹配的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回所述检索步骤；

根据传感器类型和成像特点，选择相应的几何纠正模型；

基于所述匹配的控制点对和所述选择的几何纠正模型建立像素坐标与大地坐标的转换关系；

基于数字高程模型和所述像素坐标与大地坐标之间的转换关系进行几何纠正，得到经过正射纠正的数字正射影像，

其中，对于SPOT5卫星影像选择严格成像模型，依据SPOT5提供的高精度的卫星星历数据、姿态数据以及时间数据，通过成像过程，建立ITRF90坐标系下的像点坐标与WGS84坐标系下的地面点坐标之间的数学转换关系，基于所述数字高程模型和所述像点坐标与所述地面点坐标之间的所述数学转换关系，通过基于OpenMP的单机多核并行分块处理来对所述SPOT5卫星影像进行几何纠正，

其中所述通过基于OpenMP的单机多核并行分块处理来对所述SPOT5卫星影像进行几何纠正包括：以15像素×15像素为单位对所述SPOT5卫星影像进行分块，通过计算每个图像块四个角点转换前后的坐标，建立用于该图像块的仿射变换模型，并行地对各个所述图像块单独进行纠正，之后再合成。

2.根据权利要求1所述的遥感图像快速几何纠正方法，其特征在于，所述检索包括：

估算待纠正影像的概略地理位置范围；

基于所述估算的概略地理位置范围，进行基于目标区域中心点经纬度的检索；

根据控制点影像片的属性信息进行筛选；

进行基于内容的高级检索。

3.根据权利要求2所述的遥感图像快速几何纠正方法，其特征在于，所述根据控制点影像片的属性信息进行筛选包括根据目标影像的分辨率、传感器类型、成像时间来筛选可利用的控制点影像片。

4.根据权利要求2所述的遥感图像快速几何纠正方法，其特征在于，所述基于内容的高级检索包括基于需求分布特征、颜色特征、形状特征、纹理特征来对控制点影像片进行检索。

5.根据权利要求1所述的遥感图像快速几何纠正方法，其特征在于，所述自动匹配包括：

根据控制点的坐标信息、待纠正影像的元数据信息以及成像模型计算对应的像点的初始坐标，然后按控制点影像片的大小从待纠正影像上切割出待搜索影像块；

利用Sift算法对控制点影像片以及切割后的待搜索影像块进行匹配，获得初步的匹配结果信息；

采用粗糙模糊C-均值方法和几何约束方法进行误匹配点的剔除，保留可靠的准确匹配点对；

利用最小二乘算法对匹配结果进行精匹配，得到亚像素级的匹配精度；

将匹配成功的控制点对按规定的格式输出到包括控制点号、物方坐标、像方坐标的控制点信息文件中。

6.根据权利要求5所述的遥感图像快速几何纠正方法，其特征在于，所述利用Sift算法进行匹配包括：

建立尺度空间，寻找候选点；

精确确定关键点位置，剔除不稳定点；

确定关键点梯度的模及方向；

提取特征描述符。

7.根据权利要求1所述的遥感图像快速几何纠正方法，其特征在于，在所述匹配中，采用人工辅助筛选与自动匹配相结合。

8.根据权利要求1所述的遥感图像快速几何纠正方法，其特征在于，在所述纠正之后，进行纠正结果精度检查，判断纠正结果是否满足要求，如果不满足要求，则对控制点的选择结果进行调整。