CN106056591A - 一种融合光谱图像和激光雷达数据进行城市密度估计方法 - Google Patents

Abstract

一种融合光谱图像和激光雷达数据进行城市密度估计方法，涉及数字图像处理领域。本发明要为解决现有城市密度估计方法中，存在利用的2维指标单一、3维指标少，无法合理、全面评价城市密度的问题。本发明方法按以下步骤进行：1、获取多/高光谱图像和激光雷达数据，分别对两种数据源进行预处理，利用激光雷达数据生成数字表面模型；2、在多/高光谱图像上提取光谱信息和空间信息，在激光雷达数据上提取高度信息；3、将提取得到的光谱信息、空间信息和高度信息输入到分类器中，得到分类主题图；4、用分类主题图和激光雷达提供的高度信息进行城市密度指标计算，最终生成城市密度主题图。本发明可应用于数字图像处理领域。

Description

一种融合光谱图像和激光雷达数据进行城市密度估计方法

技术领域

本发明属于数字图像处理领域，属于数字图像处理领域，尤其涉及一种融合光谱图像和激光雷达数据进行城市密度估计方法。

背景技术

多/高光谱图像(多光谱或高光谱图像)能提供目标丰富的光谱和空间信息，激光雷达能提供目标精确的高度信息。从所提供的信息角度上可以看出多/高光谱图像和激光雷达数据具有信息互补优势。联合两者数据源在农业、林业、城市发展规划等领域得到了广泛应用。

城市作为人类发展的产物，随着人口的增加，科技进步，在不断扩大。城市占地面积的增加，意味着耕地面积的减少。世界各国纷纷采用城市垂直方向发展策略来缓解这一矛盾，因此不管在城市的商业中心，还是居民区都出现了大量高楼大厦。使得仅仅利用传统的空间二维城市发展密度指标对城市密度进行评价不在合理，必须考虑城市垂直方向的高度信息。这促使了联合利用多/高光谱图像提供的光谱和空间二维信息和激光雷达数据提供的垂直方向高度信息进行城市密度估计。

为了计算城市密度指标，必须首先获取地物覆盖主题图。如果是三维(3-D)指标还必须知道对应的高度信息。用于计算城市密度指标的感兴趣区域可以是一个像素、一个网格、指定半径的圆或者别的自定义形状。如果采用单个像素作为感兴趣区域，那么只能计算二维(2-D)城市密度指标，通过光谱解混技术在亚像素水平计算特定地物所占的比例。在以别的感兴趣区域作为指标计算对象时，研究者们提出了许多2-D指标来评价城市密度：如不透水表面面积、建筑物覆盖率、植被覆盖度等；相对于2-D指标，3-D指标很少，如植被体积、建筑物体积。

目前在城市密度估计应用中存在利用的2-D指标单一、3-D指标少的问题，无法合理、全面评价城市密度。针对存在的这一问题，本发明联合利用多/高光谱图像和激光雷达数据计算能合理评价城市密度的3-D指标，并联合2-D和3-D指标对城市密度进行更合理评价。

发明内容

本发明为解决现有城市密度估计方法中，存在利用的2-D指标单一、3-D指标少，无法合理、全面评价城市密度的问题，而提出一种融合光谱图像和激光雷达数据进行城市密度估计方法。

一种融合光谱图像和激光雷达数据进行城市密度估计方法，按以下步骤进行：

一、获取多/高光谱图像和激光雷达数据，分别对两种数据源进行预处理，利用激光雷达数据生成数字表面模型(DSM)，并在两种数据源中选择控制点，对上述两种数据源进行配准；

其中，对光谱图像预处理为辐射矫正和几何矫正；对激光雷达数据预处理为奇异点剔除和图像栅格化；

二、在多/高光谱图像上提取光谱信息和空间信息，在激光雷达数据上提取高度信息：

其中光谱信息包括原始光谱波段、归一化植被指数和归一化建筑指数，空间信息包括通过使用均值、方差、形态学和Gabor空间滤波器生成的空间特征；在激光雷达数据上提取高度信息为归一化数字表面模型；

三、将提取得到的光谱信息、空间信息和高度信息输入到分类器中，得到分类主题图；

四、联合利用分类主题图和激光雷达提供的高度信息进行城市密度指标计算，最终生成城市密度主题图；

其中城市密度指标包括2-D和3-D指标，2-D指标为植被覆盖率和人造表面覆盖率；3-D指标为建筑物占地面积与使用面积之比和建筑物集成指数。

本发明包括以下有益效果：

1、由于利用了激光雷达数据提供的高度信息，可以设计合理的3-D指标，克服的目前存在的3-D指标缺少的问题；

2、进行城市密度评价时，本发明联合利用了2-D和3-D指标对城市密度进行更合理评价，克服了传统使用单一指标进行评价存在的不合理问题。

附图说明

图1为本发明所述方法流程示意图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合图1和具体实施方式对本发明作进一步详细的说明。

具体实施方式一、本实施方式所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，按以下步骤进行：

一、获取多/高光谱图像和激光雷达数据，分别对两种数据源进行预处理，利用激光雷达数据生成数字表面模型(DSM)，并在两种数据源中选择控制点，对上述两种数据源进行配准；

其中，对光谱图像预处理为辐射矫正和几何矫正；对激光雷达数据预处理为奇异点剔除和图像栅格化；

二、在多/高光谱图像上提取光谱信息和空间信息，在激光雷达数据上提取高度信息：

其中光谱信息包括原始光谱波段、归一化植被指数和归一化建筑指数，空间信息包括通过使用均值、方差、形态学和Gabor空间滤波器生成的空间特征；在激光雷达数据上提取高度信息为归一化数字表面模型；

三、将提取得到的光谱信息、空间信息和高度信息输入到分类器中，得到分类主题图；

四、联合利用分类主题图和激光雷达提供的高度信息进行城市密度指标计算，最终生成城市密度主题图；

其中城市密度指标包括2-D和3-D指标，2-D指标为植被覆盖率和人造表面覆盖率；3-D指标为建筑物占地面积与使用面积之比和建筑物集成指数。

本实施方式包括以下有益效果：

1、由于利用了激光雷达数据提供的高度信息，可以设计合理的3-D指标，克服的目前存在的3-D指标缺少的问题；

2、进行城市密度评价时，本实施方式联合利用了2-D和3-D指标对城市密度进行更合理评价，克服了传统使用单一指标进行评价存在的不合理问题。

具体实施方式二、本实施方式是对具体实施方式一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明，步骤一中所述的控制点为道路交叉点，或建筑物的拐点。

具体实施方式三、本实施方式是对具体实施方式一或二所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明，步骤一中所述的激光雷达数据奇异点是指由于空中飞行的鸟、漂浮的垃圾和风筝等物体造成的高度明显高于地物的点云，采用直方图统计方法进行激光雷达奇异点云数据剔除。

具体实施方式四、本实施方式是对具体实施方式一至三之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明，步骤二中所述的归一化数字表面模型按如下步骤计算：

1、基于移动二次曲面拟合的激光雷达图像滤波，将数字表面模型中的点集分为了地面点与非地面点两部分，具体实施步骤为：

a、选取合适的滤波窗口尺寸m×n，其中m和n的大小为10²米数量级；

b、在该窗口中选取高度最低的10个点，作为初始种子点，输入到地面点集P中，此10个点一定是地面上物体构成的点；

c、利用点集P中的点进行在二次曲面拟合，拟合所涉及的函数方程为：

$Z_i=c_0+c_1{x}_i+c_2{y}_i+c_3{x}_i{y}_i+c_4{x}_i^2+c_5{y}_i^2$

其中，(x_i,y_i)为第i个点在图像中的坐标，Z_i为该点对应的高度值；

依次将P集合内的点输入，得到一系列的方程组，在最小二乘准则下解出各系数c₀,c₁,...,c₅，从而确定曲面方程；

d、利用得到的曲面方程对其他点的高度进行预测，若预测值与实际值的差大于阈值T，则判定该点为非地面点；反之，则该点为地面点，然后将该点加入到地面点集P中，重新计算各系数，得到新曲面，如此重复直至所有点都判定完毕；

e、移动窗口至图像其它位置，完成整幅图像的滤波；

2、归一化数字表面模型的生成：

采用反距离加权插值的算法进行高度内插，插值时对距待插值点近的地面点给予较大的权值，首先以内插点为中心，确定适当数目N₀的最近邻点作为源采样点，假设内插点为S₀(x₀,y₀)，采样点为Q_i(x_i,y_i,z_i),i＝(1,2,...,N₀)，反距离加权平均插值的数学表达式如下：

$Z_{S_0}=\underset{i=1}{\overset{N_0}{\Σ}}{\mathrm{\λ}}_i{Z}_i$

其中为内插点的高度估计值；λ_i为第i个采样点的权值；Z_i为第i个采样点的高度值，λ_i和Z_i由以下公式求得：

d_i为第i个采样点到内插点的距离：

$d_i=\sqrt{{(x_i-x_0)}^2+{(y_i-y_0)}^2}$

其中p为正整数，当p取值为2时，插值的效果最好；

完成插值过程后，即生成数字地形模型(DTM)，然后用DSM与DTM相减就是归一化数字表面模型。

具体实施方式五、本实施方式是对具体实施方式一至四之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明，步骤三中所述的分类器为决策树、支持向量机或神经网络。

具体实施方式六、本实施方式是对具体实施方式一至五之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明，步骤四中所述的植被覆盖率按如下公式进行计算：

$VF=\frac{A_t+A_g}{A_{AOI}}$

其中，VF为植被覆盖率，A_t为感兴趣区域内被树木覆盖的区域面积，A_g为感兴趣区域内被草和灌木覆盖的区域面积，A_AOI为感兴趣区域的面积；

所述的人造表面覆盖率按如下公式进行计算：

$ASC=\frac{A_b+A_i}{A_{AOI}}$

其中，ASC为人造表面覆盖率，A_b为感兴趣区域内被建筑物覆盖的区域面积，A_i为感兴趣区域内被道路、广场等人造地表所占面积；

所述的建筑物占地面积与使用面积之比按如下公式计算：

$iFAR=\frac{A_{ab}}{A_{fl}}$

其中，A_ab为建筑物占地面积，A_fl为建筑物使用面积，这里假设每层楼高3米，楼层数目N按公式计算：计算，其中H_B指楼高，表示向下取整；建筑物使用面积按公式计算：其中为每层楼的面积；

所述的建筑物集成指数按如下公式计算：

$BA=\frac{\frac{A_b}{A_{AOI}}}{Median\left(D_b\right)}-\frac{Median\left(iFAR\right)}{N_b}$

其中，D_b为感兴趣区域内每两栋建筑物中心之间的距离，Median代表取中值，N_b为感兴趣区域内的建筑物数目。

具体实施方式七、本实施方式是对具体实施方式一至六之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明，步骤四中所述生成城市密度主题图的具体内容为：每一个像素为计算中心，250米为半径的圆形区域为对应中心像素的感兴趣区域进行城市密度(UD)估计，按如下公式进行计算：

UD＝(BA+ASC)-(VF+iFAR)

其中，将植被覆盖率、人造表面覆盖率、建筑物占地面积与使用面积之比和建筑物集成指数代入城市密度估计公式计算前都归一化为0到1，最终得到的UD值区间为[-2,2]，值越大代表城市密度越大。

Claims (7)

1.一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于按以下步骤进行：

一、获取多/高光谱图像和激光雷达数据，分别对两种数据源进行预处理，利用激光雷达数据生成数字表面模型DSM，并在两种数据源中选择控制点，对上述两种数据源进行配准；

其中，对光谱图像预处理为辐射矫正和几何矫正；对激光雷达数据预处理为奇异点剔除和图像栅格化；

二、在多/高光谱图像上提取光谱信息和空间信息，在激光雷达数据上提取高度信息：

其中光谱信息包括原始光谱波段、归一化植被指数和归一化建筑指数，空间信息包括通过使用均值、方差、形态学和Gabor空间滤波器生成的空间特征；在激光雷达数据上提取高度信息为归一化数字表面模型；

三、将提取得到的光谱信息和空间信息输入到分类器中，得到分类主题图；

四、联合利用分类主题图和激光雷达提供的高度信息进行城市密度指标计算，最终生成城市密度主题图；

其中城市密度指标包括2-D和3-D指标，2-D指标为植被覆盖率和人造表面覆盖率；3-D指标为建筑物占地面积与使用面积之比和建筑物集成指数。

2.如权利要求1所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于步骤一中所述的控制点为道路交叉点，或建筑物的拐点。

3.如权利要求1或2所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于步骤一中所述的激光雷达数据奇异点是指由于空中飞行的鸟、漂浮的垃圾和风筝等物体造成的高度明显高于地物的点云，采用直方图统计方法进行激光雷达奇异点云数据剔除。

4.如权利要求3所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于步骤二中所述的归一化数字表面模型按如下步骤计算：

4.1、基于移动二次曲面拟合的激光雷达图像滤波，将数字表面模型中的点集分为了地面点与非地面点两部分，具体实施步骤为：

a、选取合适的滤波窗口尺寸m×n，其中m和n的大小为10²米数量级；

b、在该窗口中选取高度最低的10个点，作为初始种子点，输入到地面点集P中，此10个点一定是地面上物体构成的点；

c、利用点集P中的点进行在二次曲面拟合，拟合所涉及的函数方程为：

Z_i＝c₀+c₁x_i+c₂y_i+c₃x_iy_i+c₄x_i ²+c₅y_i ²

其中，(x_i,y_i)为第i个点在图像中的坐标，Z_i为该点对应的高度值；

依次将P集合内的点输入，得到一系列的方程组，在最小二乘准则下解出各系数c₀,c₁,...,c₅，从而确定曲面方程；

d、利用得到的曲面方程对其他点的高度进行预测，若预测值与实际值的差大于阈值T，则判定该点为非地面点；反之，则该点为地面点，然后将该点加入到地面点集P中，重新计算各系数，得到新曲面，如此重复直至所有点都判定完毕；

e、移动窗口至图像其它位置，完成整幅图像的滤波；

4.2、归一化数字表面模型的生成：

采用反距离加权插值的算法进行高度内插，插值时对距待插值点近的地面点给予较大的权值，首先以内插点为中心，确定适当数目N₀的最近邻点作为源采样点，假设内插点为S₀(x₀,y₀)，采样点为Q_i(x_i,y_i,z_i),i＝(1,2,...,N₀)，反距离加权平均插值的数学表达式如下：

$Z_{S_0}=\underset{i=1}{\overset{N_0}{\Σ}}{\mathrm{\λ}}_i{Z}_i$

其中为内插点的高度估计值；λ_i为第i个采样点的权值；Z_i为第i个采样点的高度值，λ_i和Z_i由以下公式求得：

d_i为第i个采样点到内插点的距离：

$d_i=\sqrt{{(x_i-x_0)}^2+{(y_i-y_0)}^2}$

其中p为正整数；

完成插值过程后，即生成数字地形模型DTM，然后用DSM与DTM相减就是归一化数字表面模型。

5.如权利要求4所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于步骤三中所述的分类器为决策树、支持向量机或神经网络。

6.如权利要求5所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于步骤四中所述的植被覆盖率按如下公式进行计算：

$VF=\frac{A_t+A_g}{A_{AOI}}$

其中，VF为植被覆盖率，A_t为感兴趣区域内被树木覆盖的区域面积，A_g为感兴趣区域内被草和灌木覆盖的区域面积，A_AOI为感兴趣区域的面积；

所述的人造表面覆盖率按如下公式进行计算：

$ASC=\frac{A_b+A_i}{A_{AOI}}$

其中，ASC为人造表面覆盖率，A_b为感兴趣区域内被建筑物覆盖的区域面积，A_i为感兴趣区域内被道路、广场人造地表所占面积；

所述的建筑物占地面积与使用面积之比按如下公式计算：

$iFAR=\frac{A_{ab}}{A_{fl}}$

其中，A_ab为建筑物占地面积，A_fl为建筑物使用面积，这里假设每层楼高3米，楼层数目N按公式计算：计算，其中H_B指楼高，表示向下取整；建筑物使用面积按公式计算：其中为每层楼的面积；

所述的建筑物集成指数按如下公式计算：

$BA=\frac{\frac{A_b}{A_{AOI}}}{Median\left(D_b\right)}-\frac{Median\left(iFAR\right)}{N_b}$

其中，D_b为感兴趣区域内每两栋建筑物中心之间的距离，Median代表取中值，N_b为感兴趣区域内的建筑物数目。

7.如权利要求6所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法，其特征在于步骤四中所述生成城市密度主题图的具体内容为：每一个像素为计算中心，250米为半径的圆形区域为对应中心像素的感兴趣区域进行城市密度UD估计，按如下公式进行计算：

UD＝(BA+ASC)-(VF+iFAR)

其中，将植被覆盖率、人造表面覆盖率、建筑物占地面积与使用面积之比和建筑物集成指数代入城市密度估计公式计算前都归一化为0到1，最终得到的UD值区间为[-2,2]，值越大代表城市密度越大。