一种基于直线统计特征的大幅面遥感影像区域分类方法

摘要

本发明涉及一种基于直线统计特征的大幅面遥感影像区域分类方法。包括如下步骤：步骤1，均匀网格划分；步骤2，提取直线；步骤3，提取直线支持域；步骤4，椭圆拟合；步骤5，重新定义直线；步骤6，提取直线的长度属性；步骤7，计算直线的统计特征；步骤8，分类。充分利用不同类型的遥感影像区域具有不同的直线统计特征，能够准确地区分出大幅面遥感影像中包含的荒地、乡村、城郊和城市密集区等类型，解决了在遥感影像中盲目地搜索建筑物的问题，有利于提高建筑物检测的准确性和效率，达到完全自动化的效果。

说明书

技术领域

本发明涉及一种遥感影像处理领域，具体说是一种基于直线统计特征的大幅面遥感影像 区域分类方法。

背景技术

传统的基于像元的遥感影像区域分类方法只处在图像工程中的“图像处理”阶段，已经 不能满足现代遥感技术发展的要求。对于高分辨率遥感影像，单个像元一般不能反映真实的 地理对象，像元之间的拓扑关系也非常有限。由于基于像元的处理方法忽略了纹理、上下文 与形状等空间特征，即便将软分类器、亚像元分类方法与光谱解混技术引入其中，仅基于图 像光谱信息的分类结果也难有大幅改善。在对高分辨率遥感影像进行分类时，基于像元的方 法会导致更为严重的“椒盐效应”，从而影响分类结果的准确性。虽然传统目视判读的遥感信 息提取方法精度较高，但需要投入大量的人力和时间，己不能满足海量数据信息的处理要求。 当前，从高分辨率遥感影像中提取地理信息的手段还比较落后，较突出的问题仍然是“数据 海量、信息不足、知识难求”。另外，处在信息社会中的人们对信息的时效性越来越重视，这 有力地促进了遥感数据处理方法朝着半自动化、自动化的方向发展。因此，高效的、智能化 的遥感影像区域分类势必成为当前乃至今后遥感信息处理的研究热点。

发明内容

本发明提供了一种基于直线统计特征的大幅面遥感影像区域分类方法，可克服目前遥感 影像中建筑物提取准确率低下的问题，充分利用不同类型的遥感影像区域具有不同的直线统 计特征，能够准确地区分出大幅面遥感影像中包含的荒地、乡村、城郊和城市密集区等类型， 无需人工干预，自动化程度高。

为实现本发明的目标所采用的技术方案是：方法包括以下步骤：

步骤1：对输入遥感影像image1按均匀网格进行划分，以划分后的子图像I₁、I₂、…、I_n作为最小处理单元,n为划分后的子图像的总数；

步骤2：对子图像I₁、I₂、…、I_n分别运用直线提取算法提取直线；

步骤3：分别对子图像I₁、I₂、…、I_n中满足邻接关系的直线合并为一个直线支持域；

步骤4：用傅里叶描述子将直线支持域拟合成椭圆；

步骤5：提取椭圆的长轴，并重新定义为直线；

步骤6：提取直线的长度属性；

步骤7：计算直线的统计特征，包括直线长度的均值和熵；

步骤8：利用分类器对步骤7中的直线的统计特征进行分类。

所述的直线提取算法采用基于梯度的方法，x和y方向的梯度G_x(x,y)和G_y(x,y)分别利用以 下两个公式计算：

G_x(x,y)＝{e(x)h(y)*f(x,y)}

G_y(x,y)＝{e(y)h(x)*f(x,y)}

像素点(x，y)的梯度计算公式为：

$G(x,y)=\arctan \frac{G_x(x,y)}{G_y(x,y)}$

其中，e(k)为边缘滤波器，h(k)为映射滤波器，k的取值为x或y，对应的公式如下：

$e\left(k\right)=\frac{e^{-\mathrm{\alpha k}}}{\alpha }[-\frac{\cos (\mathrm{\alpha \beta k}+\pi /2)-\beta \sin (\mathrm{\alpha \beta k}+\pi /2)}{1+{\beta }^2}]$

h(k)＝e^-αkcos(αβk+π/2)

其中α为尺度参数，β为影像分辨率参数。

所述的邻接关系的确定方法为：选取一个阈值T(阈值T设定为5)，如果两条直线之间 的最短距离小于T，则判定这两条直接是邻接的，否则为非邻接的。

所述的用傅里叶描述子将直线支持域拟合成椭圆的方法为：利用傅里叶级数展开和欧拉公式， 对复周期函数进行变化得到如下傅里叶系数：

$a_n=\frac{1}{T}[\underset{k=0}{\overset{T-1}{\Sigma }}x\left(k\right)\cos \left(\frac{2\mathrm{\pi nk}}{T}\right)+\underset{k=0}{\overset{T-1}{\Sigma }}y\left(k\right)\sin \left(\frac{2\mathrm{\pi nk}}{T}\right)]$

${\beta }_n=\frac{1}{T}[\underset{k=0}{\overset{T-1}{\Sigma }}y\left(k\right)\cos \left(\frac{2\mathrm{\pi nk}}{T}\right)+\underset{k=0}{\overset{T-1}{\Sigma }}x\left(k\right)\sin \left(\frac{2\mathrm{\pi nk}}{T}\right)]$

其中x(k)和y(k)分别代表直线支持域边界上的点的实部和虚部。

通过提取三组系数：(α_-1,β_-1),(α₀,β₀)和(α₁,β₁)，可以得到椭圆的三个参数：

中心：c_xy＝(α₀,β₀)

长度：$L=2[\sqrt{{\alpha }_1^2+{\beta }_1^2}+\sqrt{{\alpha }_{-1}^2+{\beta }_{-1}^2}]$

方向：$\theta =\frac{\arctan \left(\frac{{\beta }_1}{{\alpha }_1}\right)+\arctan \left(\frac{{\beta }_{-1}}{{\alpha }_{-1}}\right)}{2}$

所述的直线的统计特征中直线长度的熵的计算方法为：构建一个由50个bin的直方图，每 个bin的宽度为4个像素，利用以下公式进行计算：

$E=-\underset{i=1}{\overset{50}{\Sigma }}h\left(i\right){\log }_2\left(h\left(i\right)\right)$

所述的分类器采用帕尔森窗分类器。

本发明的有益效果是：解决了在遥感影像中盲目地搜索建筑物的问题，有利于提高建筑 物检测的准确性和效率，达到完全自动化的效果。

附图说明

图1是本发明的总体处理流程图。

具体实施方式

下面结合附图详细描述本发明的具体实施方式。

在步骤101，输入的待处理遥感影像为Quick bird高空间分辨率遥感影像image1，尺寸 为8000×8000，且已经进行辐射校正和几何校正等预处理。

在步骤102，将遥感影像image1按均匀网格划分为400个400×400的子图像。

在步骤103，分别对每个400×400的子图像，采用基于梯度的方法计算梯度，x和y方向的 梯度G_x(x,y)和G_y(x,y)分别利用以下两个公式计算：

G_x(x,y)＝{e(x)h(y)*f(x,y)}

G_y(x,y)＝{e(y)h(x)*f(x,y)}

像素点(x，y)的梯度计算公式为：

$G(x,y)=\arctan \frac{G_x(x,y)}{G_y(x,y)}$

其中，e(k)为边缘滤波器，h(k)为映射滤波器，k的取值为x或y，对应的公式如下：

$e\left(k\right)=\frac{e^{-\mathrm{\alpha k}}}{\alpha }[-\frac{\cos (\mathrm{\alpha \beta k}+\pi /2)-\beta \sin (\mathrm{\alpha \beta k}+\pi /2)}{1+{\beta }^2}]$

h(k)＝e^-αkcos(αβk+π/2)

其中α为尺度参数，取1.8，β为影像分辨率参数，取70。

最后，排除梯度幅值小于10的像元，用于抑制噪声。

在步骤104，对于每个400×400的子图像在步骤103中提取的直线，选取阈值T＝5，计 算两条直线之间的最短距离L_min，如果L_min小于T，则判定这两条直接是邻接的，否则为非邻 接的，将所有邻接的直线合并为一个直线支持域。

在步骤105，用傅里叶描述子将步骤104得到的直线支持域拟合成椭圆：利用傅里叶级数 展开和欧拉公式，对复周期函数进行变化得到如下傅里叶系数：

$a_n=\frac{1}{T}[\underset{k=0}{\overset{T-1}{\Sigma }}x\left(k\right)\cos \left(\frac{2\mathrm{\pi nk}}{T}\right)+\underset{k=0}{\overset{T-1}{\Sigma }}y\left(k\right)\sin \left(\frac{2\mathrm{\pi nk}}{T}\right)]$

${\beta }_n=\frac{1}{T}[\underset{k=0}{\overset{T-1}{\Sigma }}y\left(k\right)\cos \left(\frac{2\mathrm{\pi nk}}{T}\right)+\underset{k=0}{\overset{T-1}{\Sigma }}x\left(k\right)\sin \left(\frac{2\mathrm{\pi nk}}{T}\right)]$

其中x(k)和y(k)分别代表直线支持域边界上的点的实部和虚部。

通过提取三组系数：(α_-1,β_-1),(α₀,β₀)和(α₁,β₁)，可以得到椭圆的三个参数：

中心：c_xy＝(α₀,β₀)

长度：$L=2[\sqrt{{\alpha }_1^2+{\beta }_1^2}+\sqrt{{\alpha }_{-1}^2+{\beta }_{-1}^2}]$

方向：$\theta =\frac{\arctan \left(\frac{{\beta }_1}{{\alpha }_1}\right)+\arctan \left(\frac{{\beta }_{-1}}{{\alpha }_{-1}}\right)}{2}$

在步骤106，将步骤105得到的椭圆用其长轴表示，即将椭圆拟合成直线。

在步骤107，计算直线的统计特征，包括直线长度的均值和熵。

直线长度的均值的计算公式为：

$M=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{L}_i$

其中，N为子图像中最后拟合直线的数量，L_i为第i条直线的长度

直线长度的熵的计算方法为：构建一个由50个bin的直方图，每个bin的宽度为4个像素， 利用以下公式进行计算：

$E=-\underset{i=1}{\overset{50}{\Sigma }}h\left(i\right){\log }_2\left(h\left(i\right)\right)$

在步骤108，采用帕尔森窗分类器进行分类。

在步骤109，输出分类结果：在输入遥感影像image1中，将分类结果用伪彩色进行叠加， 其中，用黄色表示荒地，绿色表示乡村，蓝色表示城郊，红色表示城市密集区。