CN103337052B - 面向宽幅遥感影像的自动几何纠正方法 - Google Patents

Abstract

基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，包括：确定待纠正影像的地理范围；针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索所有符合要求的控制点；进行自动匹配选点，确定用于几何纠正的控制点对；判断所述匹配的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回所述检索步骤；基于所述匹配的控制点对，构建不规则三角网并建立像素坐标与大地坐标的转换关系;基于小面元微分纠正方法进行几何纠正，得到经过纠正的数字正射影像。该方法实现了控制点影像自动检索和自动匹配，可减少控制点收集和选取的时间；采用了小面元几何纠正方法，可提高大幅宽影像几何纠正的精度。

Description

面向宽幅遥感影像的自动几何纠正方法

技术领域

本发明涉及一种遥感图像纠正方法，更具体而言涉及一种基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法。

背景技术

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

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

遥感图像与其他类图像的一个显著区别是，它是一种空间数据，具有空间地理位置信息。在应用遥感图像之前，必须将其投影到需要的地理坐标系中。因此，遥感图像的几何纠正处理是遥感信息处理过程中一个重要的环节，也是后续遥感影像应用的基础。

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

X=f_x(x,y,g)

Y=f_y(x,y,h)

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

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

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

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

环境减灾卫星全称中国环境与灾害监测预报小卫星星座，是我国首个专用于环境与灾害监测预报的小卫星星座，可实现灾害与环境大范围、全天候、全天时的动态监测，使我国综合减灾与环境保护工作更加科学化、现代化，为国民经济和社会持续稳定发展提供重要保障。整个星座采用分布实施的战略方案进行建设和完善，其中第一阶段构建两颗光学小卫星和一个合成孔径雷达小卫星组成的“2+1”星座；第二阶段构建四颗光学小卫星和四颗合成孔径雷达小卫星组成的“4+4”星座。现阶段已经成功发射HJ-1A、1B星。HJ-1A卫星搭载了宽覆盖的多光谱CCD相机和超光谱成像仪（HSI），HJ-1B卫星搭载了宽幅盖的多光谱CCD相机和红外相机（IRS）。

几何纠正是卫星影像进行实用化的基础，环境减灾卫星由于幅宽大（2个CCD可达720km的覆宽）、重返周期短（48小时），而且是作为灾害快速响应使用的业务化卫星，因此对几何纠正的自动化程度、处理效率等方面有着独特的需求。但环境减灾卫星影像覆盖范围大，难以用合适的多项式准确表达物点与像点之间的几何关系；姿态数据记录频率太低且姿态测量精度不高，影像畸变较大，难以通过全局的有理函数或基于姿轨参数的严格成像模型进行准确纠正。另一方面影像跨度大，地球曲率、像素大小不一致等引起的变形明显，而且在不同地区变形大小也不一致。

发明内容

根据本发明的实施例，提供了一种基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，包括：确定待纠正影像的地理范围；针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索所有符合要求的控制点；进行自动匹配选点，确定用于几何纠正的控制点对；判断所述匹配的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回所述检索步骤；基于所述匹配的控制点对，构建不规则三角网并建立像素坐标与大地坐标的转换关系；采用小面元微分纠正方法进行几何纠正，得到经过纠正的数字正射影像。

根据本发明的遥感图像纠正方法，实现了控制点影像自动检索和自动匹配、误匹配点自动剔除，从而可以减少控制点收集和选取的时间；采用了小面元几何纠正方法，从而可以提高大幅宽影像几何纠正的精度；采用了基于OpenMP的单机多核并行快速几何纠正计算，从而可以提高处理的效率，以满足减灾应用快速响应需求。

附图说明

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

图1是根据本发明的遥感图像几何纠正方法的流程图；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

图16是根据本发明实施例的分块纠正的示意图。

具体实施方式

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

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

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

图1示意性地示出了根据本发明实施例的遥感图像几何纠正方法的流程图。在步骤S1，确定待纠正影像的地理范围；在步骤S3，针对所确定的待纠正影像的地理范围，在控制点影像数据库中检索出所有符合要求的控制点；在步骤S5，进行自动匹配选点，确定用于几何纠正的控制点对；在步骤S7，判断所获得的控制点对的数量和分布是否满足几何纠正要求，如果是，则进入下一个步骤，如果不是，则返回步骤S3继续进行检索；在步骤S9，基于所匹配的控制点对，构建不规则三角网并建立像素坐标与大地坐标的转换关系；在步骤S11，基于小面元微分纠正方法进行几何纠正，得到经过纠正的数字正射影像。

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

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

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

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

控制点影像数据库构建

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

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

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

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

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

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

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

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

对已有的1：10000尺度的DOM（数字正射影像图，DigitalOrthophotoMap）、DEM（数字高程模型，DigitalElevationModel）成果图进行整理，选取现势性强、纹理清晰、特征明显的地区，如道路交叉口、桥梁、田埂角点等地区按一定大小（如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{{(a+\mathrm{\&beta;})}^2}{\mathrm{a\&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}$

θ(x,y)=αtan2((L(x,y+1)-L(x,y-1))²+(L(x+1,y)-L(x-1,y)))²

上式为(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+\&CenterDot;\&CenterDot;\&CenterDot;+{(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+\&CenterDot;\&CenterDot;\&CenterDot;{+\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)若则对于,k≠j,

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;}}}{x}_j\&Element;\stackrel{\&OverBar;}{R}{w}_i{u}_{\mathrm{ij}}^m{d}_{\mathrm{ij}}^2-\underset{j=0}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i(\underset{i=0}{\overset{c}{\mathrm{\&Sigma;}}}{x}_j\&Element;\stackrel{\&OverBar;}{R}w{u}_{\mathrm{ij}}-1)$

上式的极值条件为：

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

对上式进行计算可得：

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

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

$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_{l,j}=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_{\mathrm{lj}}}}=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。

目标函数 $J_m=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{d}_{\mathrm{ij}}^2-a\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m,$ 给出了u_ij的更新公式

$u_{\mathrm{ij}}=1/\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{((d_{\mathrm{ij}}^2-a)/(d_{\mathrm{kj}}^2-a))}^{1/(m-1)}.$ 根据此公式，以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.若则停止，否则，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）。

小面元微分纠正

小面元在遥感中主要用于同一区域两幅不同影像间的配准。首先，如前所述，在参考影像上自动提取特征点作为影像配准的控制点，通过影像匹配获得同名点对，然后由这些同名点对构成不规则三角网，再以小三角形面元为单位进行微分纠正得到精确纠正的影像。

根据本发明的实施例，不规则三角网为狄洛尼（Delaunay）三角网。

狄洛尼三角网：

区域D上有n个离散点Pi(Xi,Yi)(i=1,2,…,n)，若将D用一组直线段分成n个互相邻接的多边形，满足：

1）每个多边形内含且仅含一个离散点；

2）D中任意一点P'(X',Y')若位于Pi所在的多边形内，则满足：

$\sqrt{{(X^{\&prime;}-X_i)}^2+{(Y^{\&prime;}-Y_i)}^2}<\sqrt{{(X^{\&prime;}-X_j)}^2+{(Y^{\&prime;}-Y_j)}^2}(j\&NotEqual;i)$

若P’在与所在的两多边形的公共边上，则：

$\sqrt{{(X^{\&prime;}-X_i)}^2+{(Y^{\&prime;}-Y_i)}^2}<\sqrt{{(X^{\&prime;}-X_j)}^2+{(Y^{\&prime;}-Y_j)}^2}(j\&NotEqual;i)$

这样的多边形称为泰森多边形。用直线段连接每两个相邻多边形内的离散点而生成的三角网称为狄洛尼三角网。

对三角网中的对每一对三角形，设为ΔP₁P₂P₃和ΔP₁'P₂'P₃'，利用其三个顶点的对应坐标(X_i,Y_i)，(x_i,y_i)，i=1,2,3，求解仿射变换：

$\left\{\begin{array}{c}X=a_0+a_{1}x+a_{2}y\\ Y=b_0+b_{1}x+b_{2}y\end{array}\right.$

可得系数a₀,a₁,a₂,b₀,b₁,b₂。

然后按上式将待纠正图像上的三角形ΔP₁'P₂'P₃'纠正成与目标图像对应的三角形ΔP₁P₂P₃。

虚拟控制点:

由于构建的狄洛尼三角网不能完全覆盖原始影像，而是一个不规则的多边形，因此为了加强影像纠正后的完整性，在本发明的纠正方案中，可选地，引入原始影像的四个角点作为虚拟控制点。

按二次多项式的反解模型，建立原始影像坐标与地面坐标之间的对应关系，如下：

$\left\{\begin{array}{c}X=a_{00}+a_{10}u+a_{01}v+a_{20}{u}^2+a_{11}\mathrm{uv}+a_{02}{v}^2\\ Y=b_{00}+b_{10}u+b_{01}v+b_{20}{u}^2+b_{11}\mathrm{uv}+b_{02}{v}^2\end{array}\right.$

根据已匹配的控制点对的信息，计算多项式的系数。之后，计算四个角点（0，0），（0，img_width）,(img_height,0),(img_height,img_width)对应的地面坐标。将求出的四个控制点及与之对应的原始像点坐标作为虚拟的控制点对加入到已有的经匹配产生的控制点对的集合中，一起构建三角网，然后按小面元微分纠正的方法对原始影像进行纠正，最终得到一幅完整的经纠正后的影像。

在小面元微分纠正过程中，避开了宽幅遥感影像成像模型复杂、几何变形因素多等问题，将整幅影像的纠正划归为一个个三角面元的几何纠正，从而保证了几何纠正的精度，这个的前提就是有足够的控制点可供选择，而控制点影像库中密集的控制点影像则提供了有力的支撑。

基于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。在此基础上，为了考虑对大图像处理的需要，可以对原始图像进行分块操作，即将原图像分为若干块，每一块看作一幅图像单独进行纠正然后写入磁盘，最后合成。图16是根据本发明实施例的分块纠正的示意图。如图16所示，将待纠正的原始图像划分成如图16中左侧图所示的小块，其四个角点为a，b，c和d；然后进行成像模型变换，得到右侧图，在该右侧图中，原来的四个角点a，b，c和d对应变换为a’，b’，c’和d’，原始待纠正影像中的p点变换到纠正后的p’点。

以下给出了根据本发明实施例的基于控制点影像数据库的面向宽幅影像的自动几何纠正方法应用的例子，以及得到的实验结果。

在该示意性实施例中，采用了环境减灾卫星在中国中心部地区四景多光谱CCD影像数据，控制点影像库来源于TM的DOM影像上采集的数据。

根据待纠正影像成像地理范围，从控制点影像库中检索出该区域的控制点片数据集合，然后采用自动匹配的方式在待纠正影像上找出像点坐标，最后通过小面元微分法完成几何纠正。

表一大幅宽的环境减灾卫星遥感图像的几何纠正时间统计表

每景从自动匹配选点到完成几何纠正，总耗时不超过30分钟。纠正结果精度如下表所示，可以满足环境减灾卫星进行业务化应用的需求。

表二几何纠正结果精度统计表

通过控制点影像库，检索出待纠正影像覆盖内的控制点，采用自动影像匹配和误匹配点剔除的方法可以快速地获得控制点对，然后利用控制点通过严格成像模型，可以实现环境减灾星多光谱CCD影像的快速几何纠正，通过检查，纠正的精度结果可以达到2个像素左右，其精度可以满足环境星用于环境监测、灾害预报评估等方面的实际应用需求。另外从表一可以看出，在不改变硬件性能的前提下，如果采用并行化处理，可以有效提高计算的效率。每景环境卫星影像从影像匹配到几何纠正，总时间也不超过30分钟，非常适合于该卫星的大规模自动化生产和业务化快速应急救灾服务。

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

Claims (7)

1.一种基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，其特征在于，包括：

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

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

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

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

基于所述匹配的控制点对，构建不规则三角网并建立像素坐标与大地坐标的转换关系；

采用小面元微分纠正方法进行几何纠正，得到经过纠正的数字正射影像，

在所述构建不规则三角网中，将待纠正影像的四个角点作为虚拟控制点，按照如下二次多项式的反解模型，建立待纠正影像坐标与地面坐标之间的对应关系：

$\left\{\begin{array}{c}X=a_{00}+a_{10}u+a_{01}v+a_{20}{u}^2+a_{11}uv+a_{02}{v}^2\\ Y=b_{00}+b_{10}u+b_{01}v+b_{20}{u}^2+b_{11}uv+b_{02}{v}^2\end{array}\right.$

其中a₀₀，a₁₀，a₀₁，a₂₀，a₁₁，a₀₂，b₀₀，b₁₀，b₀₁，b₂₀，b₁₁，b₀₂为多项式的系数；根据已匹配的控制点对的信息，计算上述多项式的系数；之后，计算所述四个角点对应的地面坐标，并分别与对应的待纠正影像像点坐标构成控制点对；将由所述虚拟控制点和与之对应的待纠正影像像点构成的控制点对加入到已有的经匹配产生的控制点对的集合中；然后基于该控制点对的集合来构建所述不规则三角网。

2.根据权利要求1所述的基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，其特征在于，所述检索控制点包括：

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

基于所述估算的概略地理位置范围，对所述控制点影像数据库中存储的控制影像片进行基于目标区域中心点经纬度的检索；

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

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

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

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

5.根据权利要求1所述的基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，其特征在于，所述自动匹配包括：

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

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

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

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

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

6.根据权利要求1所述的基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，其特征在于，在所述纠正之后，进行纠正结果精度检查，判断纠正结果是否满足要求，如果不满足要求，则对控制点的检索结果进行调整。

7.根据权利要求1所述的基于控制点影像数据库的面向宽幅遥感影像的自动几何纠正方法，其特征在于，在所述纠正中，进行分块的几何纠正，并采用OpenMP进行单机多核并行处理。