CN104484647B - 一种高分辨率遥感图像云高度检测方法 - Google Patents

Abstract

本发明提出了一种高分辨率遥感图像云高度检测方法，输入同一目标的两张不同角度遥感图像即可检测出云高度。首先，将输入图像分割为地面、云以及云和地面三个部分；然后，利用SIFT匹配算法计算出地面以及云的匹配结果，并用RANSAC算法剔除地面误匹配点，得到两幅图像地面点的坐标变换关系；利用地面点的变换关系，本发明给出了一种利用循环剔除云误匹配点的方法，可以准确得到云对地面的投影偏差；最后，根据投影偏差和已知的两幅图的成像角度，利用数学模型即可得到云的高度。本发明有效地解决因卫星姿态引起的不同角度图像中同一朵云的旋转、仿射等变换所带来的匹配精度下降的问题，提高了云高度检测的精度。

Description

一种高分辨率遥感图像云高度检测方法

技术领域

本发明涉及遥感图像处理、立体视觉、特征匹配等技术领域，可应用于多角度卫星成像云高度探测以及卫星对特定目标成像能力预测等系统。

背景技术

云作为卫星遥感图像的一个重要组成部分，它的处理一直以来都是遥感图像处理的重点。云的高度信息是云的诸多显著特性之一。一方面，它可以用于短期的天气预测以及长期的气候研究。另一方面，云常常会遮挡地面目标，影响卫星对地面目标的观测能力。通过云的高度信息的探测，可以快速判断卫星能否对特定目标进行观测，使卫星能够更加高效地工作并快速获得有效信息。

利用立体视觉的方法探测云高度的具体思路在20世纪80年代被学者提出，但受当时卫星获取和处理立体视觉数据能力的限制，该方法没能成功应用。而随着卫星平台以及计算机处理数据能力的提升，近年来相关的研究成果得以发表，并付诸实践。MISR系统(多角度成像分光光度计，Multi-Angle Imaging Spectro-Radiometer)就是成功案例之一。

利用不同角度的图像探测云高度的关键在于对不同图像中的云进行特征匹配。特征匹配的准确度直接影响着云高度探测的精度。SIFT算法(尺度不变特征变换，ScaleInvariant Feature Transform)被广泛应用于图像配准和图像拼接等领域，它能够在大量数据中进行快速、准确的匹配。利用SIFT算法，能够很好地解决特征匹配的效率和精度问题。

发明内容

本发明的目的是结合立体视觉的思路和SIFT算法，从不同角度高分辨率卫星遥感图像中快速探测云的高度。

本发明提出利用SIFT算法探测高分辨率遥感图像云高度的方法，该方法包括如下步骤：

步骤1，输入同一目标在两个不同角度下的卫星遥感图像，并进行预处理；将其中一幅作为参考图像，另一幅作为检测图像。

步骤2，将参考图像分为地面、云以及云与地面过渡三个区域。

步骤3，利用SIFT算法，将两幅图进行特征提取和匹配。

步骤4，剔除地面特征误匹配点，并以地面为基准计算出两幅图的空间坐标变换关系。

步骤5，剔除云特征误匹配点，计算出两幅图中云对地面的投影偏差。

步骤6，根据步骤5的结果以及已知的两幅图成像时的角度，利用数学模型，计算出云的高度。

本发明将SIFT算法应用于高分辨率遥感图像云高度探测，可以有效地解决因卫星姿态引起的不同角度图像中同一朵云的旋转、仿射等变换所带来的匹配精度下降的问题，从而提高云高度探测的精度。此外，本发明可以不依赖人为参数选择，对任意输入的两幅视场大面积重叠、云量适中的高分辨率卫星遥感图像，只需知道它们的成像角度，即可自适应地探测出云高度。基于以上优点，本发明可以应用于多角度卫星成像云高度探测以及卫星对特定目标成像能力预测等系统。

附图说明

图1：本发明方法流程图；

图2：本发明方法的数学模型示意图；

图3：(a)成像视场的俯视示意图；(b)参考图像；(c)检测图像；

图4：(a)剔除误匹配点前地面特征匹配结果示意图；(b)剔除误匹配点后地面特征匹配结果示意图；(c)剔除误匹配点前云特征匹配结果示意图；(d)剔除误匹配点后云特征匹配结果示意图；

图5：剔除云特征误匹配点过程示意图。

具体实施方式

现结合实施例和附图，对本发明具体步骤进行进一步描述：

首先，对数学模型进行说明。对于高分辨率卫星，其视场角通常较小，如图3(a)，所以可以将卫星的物方成像光线视为平行光，而且成像地面可以视为平面；同时，两次成像的间隔很短，云的移动可以忽略不计。图2为云上一点沿两次成像光轴对地面的投影示意图。设云对地面的垂直投影点为坐标原点，向右为正方向，光轴与地面夹角定义为与负方向夹角；其中h为云高，α₁、α₂为光轴与地面夹角，δ₁、δ₂为云对地面投影点与原点的偏差(向右为正、向左为负)，|Δ|为两次成像的投影偏差大小。如图2所示，可以将两次成像分为三种情况：光轴与地面夹角均为锐角、光轴与地面夹角分别为锐角和钝角以及光轴与地面夹角均为钝角。三种情况的投影偏差大小|Δ|均可表示为则云高特别地，当α₁＝90°时，h＝|tanα₂·Δ|，当α₂＝90°时，h＝|tanα₁·Δ|。由于α₁和α₂可以通过成像时的卫星姿态得出，因此只需求出投影偏差Δ即可得到云的实际高度。在实际应用中，投影偏差可能沿图像的任意方向。为便于说明，在本发明的实施例中，假设投影偏差只发生在图像的水平方向。

图1为本发明的流程图，其具体步骤如下：

步骤1，输入同一目标在两个不同角度下的卫星遥感图像，并进行预处理；将其中一幅作为参考图像，另一幅作为检测图像。

预处理方法：均匀抽取图像中适量的行和列，以增加算法运行速度，记下原图尺寸与抽取图像尺寸的比值ratio。这里抽取行和列，而不是取平均，所以图像地面分辨率不变，检测精度不会受到影响。

处理后其中一幅作为参考图像，另一幅作为检测图像。两幅图像可由同一卫星在飞临目标上空时拍摄，也可由不同卫星拍摄，但两次拍摄时间应较短，以排除风使云移动或变形所造成的干扰。图3(b)、图3(c)分别为本实施例的参考图像和检测图像，原图大小为6144×2304，地面分辨率为3米每像素，预处理后所得图像大小1536×576，比值ratio＝4。

步骤2，将参考图像分为地面、云以及云与地面过渡三个区域；

现有方法很难既快速又精确地将地面与云区分开，所以后续剔除错误点的步骤必不可少。因此，将参考图像分为三个区域能够减小地面以及云的误判断概率，为后续处理提供方便条件。本发明利用K均值聚类的方法，其具体步骤如下：

(1)选取3个灰度值作为初始聚类中心；

为保证算法的自适应性，利用最大方差法自动选取聚类中心，步骤如下：

a.计算整幅图像灰度均值其中p_i为灰度级i的像素数在图像中所占的比例，G_max为最大灰度级；设一灰度阈值t，分别计算低灰度与高灰度像素平均灰度和

b.从0至G_max循环t，计算方差其中，

c.找出方差最小的t，作为云与地面过渡类的初始聚类中心t/2作为地面的初始聚类中心t+(G_max-t)/2作为云的初始聚类中心

(2)计算图像中每个像素灰度与这3个聚类中心的距离：并根据每个像素的d_ij,m的大小，将像素(i,j)划分至d_ij,m最小的那一类m。

(3)根据新的聚类结果，分别计算3类中像素的灰度均值作为新的聚类中心。

(4)循环步骤(2)到(3)，直至3个聚类中心的变化均小于某一阈值T。对于灰度图像，其灰度间隔一般为1，所以本发明实施例中，设置阈值T＝1。结果中灰度均值从小到大的三类依此为地面、云与地面过渡和云。

步骤3，利用SIFT算法，将两幅图进行特征提取和匹配；

该步骤所用SIFT算法，具体步骤如下：

(1)构建高斯差分(Difference of Gaussian,DOG)尺度空间

设图像为I(x,y)，则DOG尺度空间的每一层可以表示为相邻两层高斯模糊图像的差值，即

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

其中，为尺度可变高斯函数，*代表卷积运算。

(2)在构建的每一层DOG尺度空间上提取极值点

针对DOG空间上的每一点，在邻域内，检测该点是否为极值点。这里的邻域包括同一层内的8邻域点和相邻两层的9×2共26个点。

(3)针对步骤(2)中所找到的极值点，去除其中不稳定的点和DOG局部曲率非常不对称的点，保留的点即为SIFT特征点

对于空间尺度函数D(x)，其泰勒展开式为对上式求导，并令其为零，可得精确位置将其带入原式，只取前两项，可得DOG空间的极值点取值为可以认为的极值点为低对比度的和不稳定的点，将其剔除。

为剔除局部曲率不对称的点，计算保留下的极值点处的Hessian矩阵矩阵H的迹和行列式值分别记为Tr(H)＝D_xx+D_yy和若则表示该点处的局部曲率非常不对称，将其剔除，那么保留下来的点即为良好的SIFT特征点。

(4)确定特征点的方向和大小，并生成特征点描述子

利用图像的局部特征得到每一个特征点的方向和大小。对每一个特征点的邻域，计算得其梯度方向为以及该方向上的梯度大小为在0～360°范围内，以每45°一个柱，总共8个柱，做出每一个特征点的梯度方向直方图，其中的峰值代表该特征点的主方向。

对每一个特征点，将坐标轴旋转至该点的主方向上，以确保旋转不变性。再取特征点的16×16邻域，在邻域内每4×4小区域中分别计算8个梯度方向的累加值，将这个累加值向量作为一个种子点。这样，邻域内将有4×4个种子点，将这些种子点连接起来，可以得到一个4×4×8的特征向量，作为该特征点的描述子。

(5)将参考图像和检测图像的SIFT特征点进行匹配

对参考图像的每一个特征点，利用描述子，计算与检测图像中的所有特征点的欧氏距离，找出欧氏距离最近的两个点作为候选点。求出最近距离与次近距离的比值，如果该比值小于0.5，则认为匹配成功，检测图像中的最近距离点为该特征点的匹配点。记下这两个相匹配的点在各自图像中的坐标(x_1k,y_1k)和(x_2k,y_2k)，其中k代表第k个特征点。

步骤4，剔除地面特征误匹配点，并以地面为基准计算出两幅图的空间坐标变换关系；

根据步骤2中区域分割的结果，提取出地面的特征匹配点。利用随机抽样一致性(RANdom SAmple Consensus,RANSAC)算法，剔除误匹配点，提取出精确匹配点，并计算出检测图像与参考图像的坐标变换关系：

其中(x_1g,y_1g)和(x_2g,y_2g)分别为参考图像和检测图像的地面特征点在各自图像坐标系下的坐标，H_3×3为单应性矩阵，代表两坐标系的变换关系，包括平移、旋转、仿射等。图4(a)、图4(b)分别为剔除错误点前后的地面特征匹配结果。

步骤5，剔除云特征误匹配点，计算出两幅图中云对地面的投影偏差；

根据步骤2中区域分割的结果，提取出云的特征匹配点，这些特征点中包括实际云点和将地面高灰度点误认为云的点(下文称之为伪云)。因此，云特征点存在以下两类误匹配情况：1.伪云之间SIFT特征匹配正确的点：虽然它们匹配正确，但由于它们的匹配关系与实际的云不同，因此会影响对云特征点匹配关系的判断；2.云与云、伪云与伪云以及云与伪云之间SIFT特征匹配错误的特征点。

正如前文在数学模型中所说，由于云与地面存在高度差，导致两次成像时云对地面的投影产生偏差，所以地面之间以及云之间的特征匹配关系是不同的。因此，可以利用它们匹配关系的差异剔除误匹配云点，得到正确的云特征匹配关系。对于正确匹配云点，其投影偏差应集中分布在一个特定的点周围，且这个点较远离0；对于误匹配情况1，它们的投影偏差为0；对于误匹配情况2，它们的投影偏差呈随机分布。提取云匹配关系的具体步骤如下：

(1)将检测图像中的云特征点利用地面特征匹配关系进行坐标变换，将其变换至参考图像坐标系下：其中(x_2c,y_2c)和(x'_1c,y'_1c)分别为检测图像的云特征点在检测图像和参考图像坐标系下的坐标，为地面匹配关系矩阵的逆矩阵。求出两个坐标方向上云对地面的投影偏差像素数δ_pixel＝(x'_1c,y'_1c)-(x_1c,y_1c)，其中(x_1c,y_1c)为与检测图像云特征点相匹配的参考图像特征点的坐标。

(2)设定一个较小的阈值t₁，如果|δ_pixel|≤t₁，则认为其该特征点为误匹配情况1的点，将其剔除；保留|δ_pixel|＞t₁的点。实施例中t₁＝2。

(3)对所有保留下的点，求出它们投影偏差像素数的均值和两个方向上的标准差σ_x、σ_y。可以认为凡是在这个矩形范围外的点即为误匹配情况2的点，将其剔除；保留矩形范围内的点。

(4)重复步骤(3)，直至两次循环的均值的差值的两个分量均小于一个阈值t₂。图5为循环示意图，每个矩形的中心和半边长分别代表每次循环的投影偏差均值和两个方向上的标准差。对于阈值t₂的选择，计算匹配误差其中，mean代表均值运算，(x'_1g,y'_1g)为检测图像中地面点在参考图像坐标系下的坐标。如果匹配误差较小，如小于1个像素，则设置其等于1；如果匹配误差较大，则设置其等于匹配误差。

(5)求得即可表征两次成像云对地面投影偏差。

图4(c)、图4(d)分别为剔除误匹配点前后的云特征匹配结果。

步骤6，根据步骤5的结果以及已知的两幅图成像时的角度，利用数学模型，计算出云的高度。

根据数学模型，求出云高其中，GSD为图像的地面分辨率大小。特别地，当α₁＝90°时，h＝|tanα₂·Δ|，当α₂＝90°时，h＝|tanα₁·Δ|。实施例中，人为设置云高为4km，α₁＝25°，α₂＝30°，利用本发明方法测得云高为4.0002km。

Claims (2)

1.一种高分辨率遥感图像云高度检测方法，其特征在于：

步骤1，输入同一目标在两个不同角度下的卫星遥感图像，并进行预处理；将其中一幅作为参考图像，另一幅作为检测图像；所述的预处理的方法为均匀抽取输入图像中适量的行和列，并记下原图行列数与抽取图像行列数的比值ratio；

步骤2，将参考图像分为地面、云以及云与地面过渡三个区域；

步骤3，利用SIFT算法，将两幅图进行特征提取和匹配，具体流程如下：

(1)构建高斯差分尺度空间；

(2)在构建的每一层构建高斯差分尺度空间上提取极值点；

(3)针对所找到的极值点，去除其中不稳定的点和构建高斯差分局部曲率非常不对称的点，保留的点即为SIFT特征点；

(4)确定特征点的方向和大小，并生成特征点描述子；

(5)将参考图像和检测图像的SIFT特征点进行匹配；

步骤4，剔除地面特征误匹配点，并以地面为基准计算出两幅图的空间坐标变换关系，记下矩阵H_3×3代表地面匹配关系；

步骤5，剔除云特征误匹配点，计算出两幅图中云对地面的投影偏差，其流程如下：

(1)将检测图像中的云特征点变换至参考图像坐标系下：其中(x_2c,y_2c)和(x'_1c,y'_1c)分别为检测图像的云特征点在检测图像和参考图像坐标系下的坐标，为地面匹配关系矩阵的逆矩阵；求出两个坐标方向上云对地面的投影偏差像素数δ_pixel＝(x'_1c,y'_1c)-(x_1c,y_1c)，其中(x_1c,y_1c)为与检测图像云特征点相匹配的参考图像特征点的坐标；

(2)设定一个阈值t₁，如果|δ_pixel|≤t₁，将其剔除；

(3)对所有保留下的点，求出它们投影偏差像素数的均值和两个方向上的标准差σ_x、σ_y；对于在这个矩形范围外的点将其剔除，保留矩形范围内的点；

(4)重复步骤(3)，直至两次循环的投影偏差均值的差值的两个分量均小于一个阈值t₂；

(5)求得即可表征两次成像云对地面投影偏差；

步骤6，根据步骤5的结果以及已知的两幅图成像时的角度，利用数学模型，计算出云的高度，计算云的高度所用公式和数据如下：

云高其中，GSD为图像的地面分辨率大小，α₁、α₂为输入图像的成像角度。

2.根据权利要求1所述方法，其特征在于，阈值t₂的选择方法为：

计算匹配误差，若匹配误差小于1，则t₂＝1；否则，t₂等于匹配误差。