CN103218787B - 多源异构遥感影像控制点自动采集方法 - Google Patents

Abstract

一种多源异构遥感影像控制点自动采集方法，包括：从经过几何精纠正处理的多个数据源获取遥感影像；对所述获取的遥感影像进行分析归纳，统计并分析其属性信息；对所述获取的遥感影像进行优化格网设计，根据影像分辨率、幅宽将影像分为不同尺寸的格网；利用Wallis变换对每一个格网内的影像进行预处理；对经Wallis变换处理后的影像进行增强，综合运用Moravec算子、Harris算子、Forstner算子、SUSAN算子、尺度不变特征及最稳定极值区域检测算法提取特征点；以控制点影像片采集窗口来遍历统计特征点的个数，选择特征点数量最大的窗口为采集区域；针对所述选择的采集区域进行遥感控制点影像片的裁切保存，同时获取其对应区域的数字高程模型数据。

Description

多源异构遥感影像控制点自动采集方法

技术领域

本发明涉及一种遥感影像控制点采集方法，更具体而言涉及一种多源异构遥感影像控制点自动采集方法。

背景技术

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

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

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

所谓正射纠正，就是基于GCP及DEM信息，对在进行航天、航空摄影时，由于无法保证摄影瞬间航摄像机的绝对水平，得到的影像是一个倾斜投影的像片，像片各个部分的比例尺不一致问题进行影像倾斜纠正，同时对光学成像相机中心投影成像的方式，由地面的高低起伏在像片上形成投影差的问题进行投影差的改正，消除各种变形的纠正过程。

在遥感影像正射纠正系统的构建中，GCP的获取是最关键的步骤，GCP的精度直接影响校正后结果的几何精度。

传统的摄影测量中，GCP的获取是很麻烦的：在像片上选择明显的地物点，再到野外通过三角测量或其它大地测量方法来获取这些点的精确空间位置，由此得到一组控制点对。这种传统的GCP采集方法，虽然精度较高，但是采集周期长，作业成本高。

此外，目前，一些遥感图像处理软件，在满足一定精度的应用要求下，可以从相应比例尺的地图上量取控制点。但是，通过地形图资料获取的GCP存在地形图现势性、人工判读误差等问题。

发明内容

为了应对常规控制点获取技术中周期长、应用效率低、可重复利用率低的缺点，本发明提供了一种多源异构遥感影像控制点自动采集方法，该方法包括：从经过几何精纠正处理的多个数据源获取遥感影像；对所述获取的遥感影像进行分析归纳，统计并分析其属性信息，该属性信息包括影像格式、影像分辨率、坐标系统、成图比例尺、影像时相；对所述获取的遥感影像进行优化格网设计，根据影像分辨率、幅宽将影像分为不同尺寸的格网；利用Wallis变换对每一个格网内的影像进行预处理；对经Wallis变换处理后的影像进行增强，综合运用Moravec算子、Harris算子、Forstner算子、SUSAN算子、尺度不变特征及最稳定极值区域检测算法提取特征点；以控制点影像片采集窗口来遍历统计特征点的个数，选择特征点数量最大的窗口为采集区域；针对所述选择的采集区域进行遥感控制点影像片的裁切保存，同时获取其对应区域的数字高程模型数据。

通过该方法能够从多源异构影像自动提取遥感控制影像点，提高了控制数据获取的效率和精度，为遥感影像的自动化、智能化处理提供了基础资料。

附图说明

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

图1是根据本发明实施例的遥感影像控制点自动采集方法的流程图；

图2a和图2b分别示出了Wallis变换前的影像和Wallis变换后的影像；

图3是Moravec算子的示意图；

图4是尺度空间生成的示意图；

图5是空间极值点检测的示意图；

图6是特征描述符生成的示意图；

图7示意性地示出了在极值区域ε的MSER检测中拟合的椭圆区域；

图8示出了对于示例性图像的SIFT特征点描述。

具体实施方式

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

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

根据本发明的实施例，提供了一种多源异构遥感控制点影像片自动采集方法。通过多源异构遥感影像的数据源分析，综合影像格式、影像分辨率、坐标系统、成图比例尺等属性信息，构建遥感影像特征点自动提取和描述算法，利用关键点的邻域像素的梯度模值和方向分布特性，自动采集遥感影像控制点，以解决遥感影像控制点人工采集费时费力、特征不适合机器匹配且实时性差的问题。

根据本发明的技术方案，多源指数据来源于不同的飞行平台和传感器；异构指不同的成像方式，如框幅式单中心投影、推扫式多中心投影、雷达成像等；控制点影像片，是指为了满足多源遥感影像的几何纠正处理，通过机器自动匹配或者辅之以人工选择，按照相应的规范要求，从已经经过精确几何纠正处理的正射影像上裁剪的具有地理坐标信息、纹理信息的影像块，该影像块的每个像素都具有地理坐标，都可以作为控制点使用。

图1示出了根据本发明实施例的遥感图像控制点采集方法的流程。

如图1所示，在步骤S01，获取遥感图像采集数据，该数据来自于经过几何精纠正处理的多个数据源，这些数据源例如，航空影像、资源系列卫星影像、SPOT系列卫星影像、雷达卫星影像等。

在步骤S02，对在步骤S01获取的遥感影像进行分析归纳，统计其影像格式、影像分辨率、坐标系统、成图比例尺、影像时相等属性信息，从而对原始采集数据（原始影像）进行综合信息分析

具体而言，尽管如前所述用于遥感影像控制点采集的经纠正影像是具有地理坐标的影像，但是，由于生产单位不同、数据源不同（传感器不同）、时相不同、坐标系信息不同，需要对这些属性信息进行归纳，并在采集的过程中将这些信息与所采集的控制点影像片相关联，以便于输入和更新控制点影像数据库。

在步骤S03，进行优化格网设计，根据影像的分辨率、幅宽等信息将影像分为不同尺寸的格网，例如，3×3、5×5或7×7等，从而可以提高所提取的遥感影像控制点的空间分布的合理性。

网格的设计尺寸主要是由具体应用的需求所确定的，一般按照理论要求一般选择3×3的既可以满足需求，但是根据测绘数据处理的要求一般会选择冗余条件，以便于检核，如有相关要求即需要选择大的格网；同时航空影像和航天影像也有不同的要求，一般航空影像处理需要的控制点较多。

在步骤S04，利用Wallis变换对每一个格网内的影像进行预处理。Wallis滤波器在计算影像的局部灰度方差和均值时使用平滑算子，所以其可以在增强影像有用信息的同时抑制噪声，提高影像的信噪比，使影像中存在的模糊的纹理模式得到增强，有利于特征提取。

Wallis变换可以表示为：

$f(x,y)=[g(x,y)-m_g]\frac{{\mathrm{cs}}_f}{{\mathrm{cs}}_g+(1-c)s_f}+{\mathrm{bm}}_f+(1-b)m_g---\left(1\right)$

式（1）中，g(x,y)为原始影像的灰度值；f(x,y)为Wallis变换结果影像的灰度值；m_g为原始影像的局部灰度均值；s_g为原始影像的局部灰度方差值；m_f为结果影像局部灰度均值的目标值；s_f为结果影像的局部灰度方差值的目标值；c∈[0,1]，为影像方差的扩展常数；b∈[0,1]，为影像的亮度系数，当b→1时，影像均值被强制到m_f，当b→0时，影像均值被强制到m_g。

式（1）也可表示为：

f(x,y)=g(x,y)r₁+r₀                       （2）

式（2）中，r₀=bm_f+(1-b-r₁)m_g中的参数r₁,r₀分别为乘性系数和加性系数。

Wallis变换的实现过程如下：

1）对数字图像分为互不重叠矩形区域，每区域的尺度对应于要增强的纹理模式的尺度；

2）计算每一块的灰度均值和方差；

3）灰度均值和方差的目标值分别设定为127和40-70之间的数值，其中后者可以根据区域尺度的减小而减小以防止大量像素的灰度值被饱和（即落于[0,255]之外）；

4）计算出每一块的Wallis变换的乘性常数r1和加性常数r0。

5）由于各矩形区域不重叠，所以数字图像的任一像素的系数由双线性内插得到，并根据公式（2）计算变换后的灰度值f(x,y)。

图2a和图2b分别示出了Wallis变换前的影像和Wallis变换后的影像。

在步骤S05，对经Wallis变换处理后的影像进行增强。

在每一景影像上，由于获取该影像的传感器、分辨率、时相等均不同，所以影像反映出的特征点也不同（所谓特征点例如房屋角点、高亮点等），这些特征点是后期人工刺点和进行机器自动匹配的基础，特征点的多少直接反应匹配点数量和成功率。为保证采集的控制点影像片内包含足够的影像特征点，即，既包含足够传统灰度匹配的特征点，又包含足够的满足光照、尺度、旋转等不变特征的特征点，综合运用Moravec算子、Harris算子、Forstner算子、SUSAN算子、尺度不变特征及最稳定区域检测等多种算法提取特征点，即是先利用Moravec算子、Harris算子、Forstner算子、SUSAN算子提取特征点，然后利用尺度不变特征及最稳定区域算法来检测所提取到的特征点，记录特征点个数，并且利用描述算法对提取出的特征点进行具体描述，以确定其特征的精确性和稳定性。

a.Moravec算子检测步骤：

首先，计算各格网区域影像的兴趣值IV。在以像素(c,r)为中心的w×w的影像窗口中（w限定窗口的大小，一般选用奇数，如5×5的窗口），计算图3所示四个方向相邻像素灰度差的平方和：

$V_1=\underset{i=-k}{\overset{k-1}{\mathrm{\Σ}}}{(g_{c+i,r}-g_{c+i+1,r})}^2$

$V_2=\underset{i=-k}{\overset{k-1}{\mathrm{\Σ}}}{(g_{c+i,r+i}-g_{c+i+1,r+i+1})}^2$

$V_3=\underset{i=-k}{\overset{k-1}{\mathrm{\Σ}}}{(g_{c,r+i}-g_{c,r+i+1})}^2$

$V_4=\underset{i=-k}{\overset{k-1}{\mathrm{\Σ}}}{(g_{c+i,r-i}-g_{c+i+1,r-i-1})}^2---\left(3\right)$

其中k=INT(w/2)，g表示灰度值。取其中最小者为该像素(c,r)的兴趣值，即

IV_c,r=min(V₁,V₂,V₃,V₄)        （4）

接下来，给定一个经验阈值，将兴趣值大于该阈值的点（即归一化后的兴趣值计算窗口的中心点）作为候选点。阈值的选择应该以候选点中包括所需要的特征点而又不含过多的非特征点为原则。

之后，选取候选点中的极值点作为特征点。在一定大小的窗口内（可不同于兴趣值计算窗口，例如5×5像元，7×7像元等），将候选点中兴趣值不是最大者均去掉，仅保留兴趣值最大者，该像素即为一个特征点。通过该极值点筛选操作，可以抑制局部非最大。

b.Harris算子的计算步骤如下：

首先，求出影像上所有像素点在x，y方向的梯度g_x,g_y，即对每个像素点的灰度值进行一阶差分运算。

g_x=g_i,j-g_i,j+1;g_y=g_i,j-g_i+1,j     （5）

接下来，确定一个n×n大小的窗口，利用上式（5）生成n×n大小的高斯卷积模板。其中，高斯模板的σ取0.3-0.9之间的值。

$G(x,y)=e^{\frac{x^2+y^2}{{2\mathrm{\σ}}^2}}---\left(6\right)$

然后，利用生成的高斯模板对梯度值g_x,g_y进行高斯滤波，并计算强度值I。

之后，选取局部极值点为特征点，如果I大于阈值0.6，则将该点取为特征点，并且进行排序保存。

c.Forstner算子的计算步骤如下：

首先，计算各像素的Robert梯度：

$g_u=\frac{\mathrm{\δg}}{\mathrm{\δu}}=g_{i+1,j+1}-g_{i,j}$

$g_v=\frac{\mathrm{\δg}}{\mathrm{\δv}}=g_{i,j+1}-g_{i+1,j}---\left(7\right)$

接下来，计算l×l窗口中灰度的协方差矩阵：

$Q=N^{-1}={\left[\begin{array}{cc}\mathrm{\Σ}{g}_u^2& \mathrm{\Σ}{g}_u{g}_v\\ \mathrm{\Σ}{g}_v{g}_u& \mathrm{\Σ}{g}_v^2\end{array}\right]}^{-1}---\left(8\right)$

其中，

$\mathrm{\Σ}{g}_u^2=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\Σ}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\Σ}}}{(g_{i+1,j+1}-g_{i,j})}^2$

$\mathrm{\Σ}{g}_v^2=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\Σ}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\Σ}}}{(g_{i,j+1}-g_{i+1,j})}^2$

$\mathrm{\Σ}{g}_u{g}_v=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\Σ}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\Σ}}}(g_{i+1,j+1}-g_{i,j})(g_{i,j+1}-g_{i+1,j})$

然后，计算兴趣值q与w：

$w=\frac{1}{\mathrm{trQ}}=\frac{\det N}{\mathrm{trN}}$

$q=\frac{4\det N}{{\left(\mathrm{trN}\right)}^2}---\left(9\right)$

其中，detN代表矩阵N的行列式；trN代表矩阵N的直迹。

之后，确定特征点：

如果兴趣值大于给定的阈值，则该像元为特征点。阈值为经验值，可使用下式（10）作为参考公式：

其中，为权的平均值；w_c为权的中值。

当q>T_q且w>T_w时，选该像元为特征点。

d.SUSAN(Small univalue segment assimilating nucleus)算子是一种基于灰度的特征点获取方法，适用于图像中边缘和角点的检测。所谓角点，是指二维图像亮度变化剧烈的点或图像边缘曲线上曲率极大值的点，这些点在保留图像图形重要特征的同时，可以有效地减少信息的数据量。SUSAN角点检测的步骤如下：

首先，用SUSAN模板在图像上滑动，将模板内部每个图像像素点的灰度值分别与模板中心像素的灰度值作比较，用下式来计算模板中各点的相似度。

上式（11）中，为模板内属于USAN区域（即，吸收核同值区，Univalue Segment Assimilating Nucleus,USAN)。）的相似比较函数，是模板中心在图像中的位置；是模板内除中心之外的任意一个点的位置；的灰度；的灰度；t为灰度差的阈值。通常采用更稳定可靠的指数函数形式，如下式（12）所示：

$c(\stackrel{\→}{r},{\stackrel{\→}{r}}_0)=e^{-{\left(\frac{I\left(\stackrel{\→}{r}\right)-I\left({\stackrel{\→}{r}}_0\right)}{t}\right)}^6}---\left(12\right)$

接下来，用下式（13）计算图像中每一点的USAN区域大小。

$n\left({\stackrel{\→}{r}}_0\right)=\underset{\stackrel{\→}{r}\∈C\left({\stackrel{\→}{r}}_0\right)}{\mathrm{\Σ}}c(\stackrel{\→}{r},{\stackrel{\→}{r}}_0)---\left(13\right)$

然后，用下式（14）的角点响应函数来计算角点初始响应：

其中，几何阈值g控制角点的生成质量，一般取3×n_max/4，n_max是所能达到的最大值。g的大小主要影响特征点的尖锐性，g越小，所检测到的角度越尖锐。灰度差阈值t用来控制角点密度，t越小，检测的图像的细微差别越多，同样角点数越多。

之后，通过USAN区域的质心与模板中心的距离去除伪角点。对角点响应图像中的正值所对应的点，如果不满足以下（i）和（ii）两点要求，则剔除该点。

（i）对于那些与正确角点相关的USAN区域，其质心位置远离模板的中心位置，质心的计算公式为：

$\stackrel{\→}{r}\left({\stackrel{\→}{r}}_0\right)=\frac{\underset{\stackrel{\→}{r}}{\mathrm{\Σ}}\stackrel{\→}{r}c(\stackrel{\→}{r},{\stackrel{\→}{r}}_0)}{\underset{\stackrel{\→}{r}}{\mathrm{\Σ}}c(\stackrel{\→}{r},{\stackrel{\→}{r}}_0)}---\left(15\right)$

（ii）对模板内位于从模板中心指向USAN区域质心的直线上的所有像素都必须是USAN区域的一部分。

再后，通过非极大值抑制来寻找角点。

因为，在角点附近区域，响应函数可能都具有较高的值，而只有在局部有最大值的位置才是角点的正确位置，所以为了得出角点的正确位置必须进行非极大值抑制。

e.尺度不变特征变换（Scale Invariant Feature Transform，SIFT），用于描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。由于SIFT特征与影像的大小和旋转无关，对于光线、噪声、微小视角改变的容忍度也相当高，因而利用少量的SIFT物体特征就足以计算出位置与方位，适合在海量数据库中快速准确匹配。

SIFT特征提取步骤主要包括：建立尺度空间，选取候选点，确定关键点精确位置，删除不稳定点，关键点方向分配以及关键点描述。

（1）尺度空间的生成

首先，如下式（16）定义二维高斯函数：

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

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

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

L(x,y,σ)=G(x,y,σ)*I(x,y)    （17）

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

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

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

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

（2）空间极值点检测

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

（3）关键点位置确定

利用尺度空间函数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---\left(19\right)$ （其中Χ=(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]---\left(20\right)$

然后，进行低对比度剔除

由如果|D（X）|<0.03，则该点对比度较低，将其剔除。

（4）边缘响应的去除

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

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

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

具体而言，在数字影像中，例如，导数的计算可以用差分予以近似估计， $D_x=\frac{\&PartialD;g}{\&PartialD;x}=g_{i+1,j+1}-g_{i,j};$ $D_y=\frac{\&PartialD;g}{\&PartialD;y}=g_{i,j+1}-g_{i+1,j},$ 其中g表示灰度值，i,j表示像素的行列号）。

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

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

Det(H)=D_xxD_yy-(D_xy)²=αβ     （21）

令α=γβ，则：

$\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}{r{\mathrm{\&beta;}}^2}=\frac{{(r+1)}^2}{r}---\left(22\right)$

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

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

取r＝10。

对候选点进行筛选后，还要确定所确定的控制点的梯度模、方向及特征描述。

（5）关键点梯度模及方向计算

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

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

θ(x,y)=α tan 2((L(x,y+1)-L(x,y-1))²+(L(x+1,y)-L(x-1,y)))²     （24）

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

（6）特征描述符生成

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

f.最稳定极值区域（Maximally Stable Extremal Region，MSER）检测

MSER是图像中的一些极值区域，而这样的极值区域ε有着一些非常好的特性。首先，图像本身灰度的单调变化并不会改变这些极值区域ε，因为极值区域是通过比较区域内像素点和区域边界上像素点的灰度值而生成的，只是一种相对关系的结果，而图像灰度的单调变化并不会改变像素点之间的相对灰度大小关系，因此也就不会改变生成的ε了。这种特性保证了在光照线性变化以后仍能检测出同样的极值区域ε，从而消除了光照变化对图像匹配的影响。其次，对图像进行连续的几何变换以后，图像中像素点的拓扑关系并没有改变，即像素点会从一个连通区域变换到同一个连通区域。因此即使对图像进行了仿射变换甚至连续的非线性扭曲以后，得到的变换以后的区域ε'仍然将会是与ε相匹配的极值区域。最后，由于所检测到的极值区域的数量总是小于图像中所有像素点的总数量（当且仅当每个像素点都是一个极值区域时极值区域的数量等于所有像素点的总数量），这保证了整个检测过程不会有太高的运算复杂度。

极值区域ε的MSER检测的具体过程如下：

首先，检测出图像中所有的极值区域ε，这个过程分为以下几步：

（1）将所有的像素点按照灰度值进行排列。由于每个像素点的灰度值都是在0-255的范围内，因此可以使用O(n)复杂度的计数排序算法来实现（其中n为像素点总数，下同）。

（2）在将所有像素点排序好以后，按照升序或者降序将这些像素点在图像中标记出来，在标记的过程中，使用union-find算法将一系列不断增大或者出现的连通区域和它们的面积保存下来，并可以按照树的数据结构进行保存，这样的运算复杂度为O(nα(n))，其中α为反Ackerman函数且对于有限的n均有α(n)≤4，因此它十分接近线性复杂度。而这样的一系列连通区域就是极值区域ε了。

极值区域找到后，进行拟合将其拟合为一个近似的椭圆区域，然后寻找拟合的椭圆区域的重心，如图7所示。然后，以该重心点坐标作为虚拟特征点，利用SIFT匹配算法中的描述方法进行特征点描述。图8示出了对于示例性图像的SIFT特征点描述。

在根据本发明的实施例中，例如，首先根据设计要求设定遥感影像控制点采集的数量，即控制点片的尺寸，例如，待采集影像的影像分辨率为0.05-1.0m，采集控制点影像片为512×512个像素；影像分辨率为1.0-5.0m，采集控制点影像片为1024×1024个像素；影像分辨率大于5.0m，采集控制点影像片为2048×2048个像素。然后，如前所述，以此为窗口到构建的相应格网范围内进行遍历，并统计该范围内的所有特征点的数目,包括综合运用Moravec算子、Harris算子、Forstner算子、SUSAN算子、尺度不变特征及最稳定区域检测算法提取的所有特征点，记为数组Feature_NUM(i)，之后，对Feature_NUM(i)数据进行排序。

在步骤S06，以控制点影像片采集窗口来遍历统计特征点个数，选择特征点数量最大的窗口为采集区域。

在步骤S07，针对所选择的采集区域进行遥感控制点影像片的裁切保存，同时获取其对应区域的数字高程模型数据。于是，采集得到所需的遥感影像控制点数据及对应的数字高程模型数据。

通过对已有的经几何精确纠正的影像进行控制点影像片的采集，能够获得可供多次反复使用的控制点信息，从而解决野外测量的点位坐标无法重复使用的问题。

根据本发明的遥感影像控制点自动采集技术，改善了现有基于人工和半自动的提取技术，保证了采集影像特征点的丰富性，提高了工作效率，为遥感影像自动化处理提供了基础。以上所述仅是本发明的示范性实施方式，而非用于限制本发明的保护范围，本发明的保护范围由所附的权利要求确定。

Claims (10)

1.一种多源异构遥感影像控制点自动采集方法，其特征在于，包括：

从经过几何精纠正处理的多个数据源获取遥感影像；

对所述获取的遥感影像进行分析归纳，统计并分析其属性信息，该属性信息包括影像格式、影像分辨率、坐标系统、成图比例尺、影像时相；

对所述获取的遥感影像进行优化格网设计，根据影像分辨率、幅宽将影像分为不同尺寸的格网；

利用Wallis变换对每一个格网内的影像进行预处理；

对经Wallis变换处理后的影像进行增强，综合运用Moravec算子、Harris算子、Forstner算子、SUSAN算子、尺度不变特征及最稳定极值区域检测算法提取特征点；

以控制点影像片采集窗口来遍历统计所述特征点的个数，选择特征点数量最大的窗口为采集区域；

针对所述选择的采集区域进行遥感控制点影像片的裁切保存，同时获取其对应区域的数字高程模型数据。

2.根据权利要求1所述的影像控制点自动采集方法，其特征在于，所述数据源包括航空影像、资源系列卫星影像、SPOT系列卫星影像、雷达卫星影像中的一个或多个。

3.根据权利要求1所述的影像控制点自动采集方法，其特征在于，所述Wallis变换包括：

将数字图像分为互不重叠矩形区域，每区域的尺度对应于要增强的纹理模式的尺度；

计算每一块区域的灰度均值和方差；

将灰度均值和方差的目标值分别设定为127和40-70之间的数值；

计算出每一块所述区域的Wallis变换的乘性常数和加性常数；

由双线性内插计算所述数字图像的任一像素的系数，并计算变换后的灰度值。

4.根据权利要求1所述的影像控制点自动采集方法，其特征在于，所述综合运用Moravec算子、Harris算子、Forstner算子、SUSAN算子、尺度不变特征及最稳定极值区域检测算法提取特征点包括：先利用Moravec算子、Harris算子、Forstner算子、SUSAN算子提取特征点，然后利用尺度不变特征变换及最稳定极值区域算法来检测所提取到的特征点，记录特征点个数，并且利用描述算法对提取出的特征点进行具体描述。

5.根据权利要求4所述的影像控制点自动采集方法，其特征在于，所述Moravec算子的计算包括：

计算各格网区域影像的兴趣值；

给定一个阈值，将所述计算的兴趣值大于该阈值的点作为候选点；

选取所述候选点中的极值点作为特征点。

6.根据权利要求4所述的影像控制点自动采集方法，其特征在于，所述Harris算子的计算包括：

求出影像上所有像素点的梯度，即对每个像素点的灰度值进行一阶差分运算；

确定一个n×n大小的窗口，生成n×n大小的高斯卷积模板，其中，高斯模板的σ取0.3-0.9之间的值；

利用生成的高斯模板对所述梯度值进行高斯滤波，并计算强度值；

选取所述强度值的局部极值点为特征点，如果所述强度值大于阈值0.6，则将该点取为特征点，并且进行排序保存。

7.根据权利要求4所述的影像控制点自动采集方法，其特征在于，所述Forstner算子的计算包括：

计算各像素的Robert梯度；

根据所计算的Robert梯度值来计算l×l窗口中灰度的协方差矩阵；

利用所计算的协方差矩阵来计算兴趣值；

确定特征点，其中，如果兴趣值大于给定的阈值，则该像素为特征点。

8.根据权利要求4所述的影像控制点自动采集方法，其特征在于，所述SUSAN算子的计算包括：

用SUSAN模板在图像上滑动，将模板内部每个图像像素点的灰度值分别与模板中心像素的灰度值作比较，计算模板中各点的相似度；

基于所计算的相似度来计算图像中每一点的USAN区域大小；

基于所计算的USAN区域大小，用角点响应函数来计算角点初始响应；

通过所述USAN区域的质心与所述模板中心的距离去除伪角点；

在去除伪角点后，再通过非极大值抑制来寻找角点。

9.根据权利要求4所述的影像控制点自动采集方法，其特征在于，所述尺度不变特征变换包括：

建立尺度空间；

检测空间极值点作为候选点；

确定关键点位置；

去除边缘响应，删除不稳定点；

计算所述关键点的梯度模及方向；以及

生成特征描述符来描述所述关键点。

10.根据权利要求4所述的影像控制点自动采集方法，其特征在于，所述最稳定极值区域检测包括：

检测图像中所有的极值区域ε，其中，将所有的像素点按照灰度值进行排列，在将所有像素点排序以后，按照升序或者降序将这些像素点在图像中标记出来，在标记的过程中，使用union-find算法将一系列不断增大或者出现的连通区域和它们的面积保存下来，这样的一系列连通区域就是所述极值区域ε；

检测出所述极值区域后，进行拟合将其拟合为一个近似的椭圆区域，然后确定拟合的椭圆区域的重心；

以所述重心的坐标作为虚拟特征点，利用尺度不变特征变换算法中的描述方法进行特征点描述。