CN106056591A - Method for estimating urban density through fusion of optical spectrum image and laser radar data - Google Patents

Method for estimating urban density through fusion of optical spectrum image and laser radar data Download PDF

Info

Publication number
CN106056591A
CN106056591A CN201610356309.3A CN201610356309A CN106056591A CN 106056591 A CN106056591 A CN 106056591A CN 201610356309 A CN201610356309 A CN 201610356309A CN 106056591 A CN106056591 A CN 106056591A
Authority
CN
China
Prior art keywords
point
laser radar
building
area
radar data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201610356309.3A
Other languages
Chinese (zh)
Other versions
CN106056591B (en
Inventor
谷延锋
王青旺
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Heilongjiang Industrial Technology Research Institute Asset Management Co ltd
Original Assignee
Harbin Institute of Technology Shenzhen
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology Shenzhen filed Critical Harbin Institute of Technology Shenzhen
Priority to CN201610356309.3A priority Critical patent/CN106056591B/en
Publication of CN106056591A publication Critical patent/CN106056591A/en
Application granted granted Critical
Publication of CN106056591B publication Critical patent/CN106056591B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01JMEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J3/00Spectrometry; Spectrophotometry; Monochromators; Measuring colours
    • G01J3/28Investigating the spectrum
    • G01J3/2823Imaging spectrometer
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88Lidar systems specially adapted for specific applications
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01JMEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J3/00Spectrometry; Spectrophotometry; Monochromators; Measuring colours
    • G01J3/28Investigating the spectrum
    • G01J3/2823Imaging spectrometer
    • G01J2003/2826Multispectral imaging, e.g. filter imaging
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • G06T2207/10036Multispectral image; Hyperspectral image
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20081Training; Learning

Landscapes

  • Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Image Processing (AREA)

Abstract

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

A method for estimating urban density by fusing spectral images and lidar data relates to the field of digital image processing. The purpose of the present invention is to solve the problem that the existing urban density estimation method has a single 2-dimensional index and few 3-dimensional indexes, and cannot reasonably and comprehensively evaluate the urban density. The method of the present invention is carried out according to the following steps: 1, obtain multi/hyperspectral image and laser radar data, carry out pretreatment to two kinds of data sources respectively, utilize laser radar data to generate digital surface model; 2, extract on multi/hyperspectral image Spectral information and spatial information, and extract height information from the lidar data; 3. Input the extracted spectral information, spatial information and height information into the classifier to obtain a classified theme map; 4. Use the classified theme map and lidar to provide The height information of the city is used to calculate the urban density index, and finally generate the urban density theme map. The invention can be applied to the field of digital image processing.

Description

一种融合光谱图像和激光雷达数据进行城市密度估计方法A Method for Urban Density Estimation by Fusion of Spectral Imagery and LiDAR Data

技术领域technical field

本发明属于数字图像处理领域,属于数字图像处理领域,尤其涉及一种融合光谱图像和激光雷达数据进行城市密度估计方法。The invention belongs to the field of digital image processing, in particular to a method for estimating urban density by fusing spectral images and laser radar data.

背景技术Background technique

多/高光谱图像(多光谱或高光谱图像)能提供目标丰富的光谱和空间信息,激光雷达能提供目标精确的高度信息。从所提供的信息角度上可以看出多/高光谱图像和激光雷达数据具有信息互补优势。联合两者数据源在农业、林业、城市发展规划等领域得到了广泛应用。Multi/hyperspectral images (multispectral or hyperspectral images) can provide rich spectral and spatial information of targets, and lidar can provide accurate height information of targets. From the perspective of the information provided, it can be seen that multi/hyperspectral images and LiDAR data have complementary information advantages. Combining the two data sources has been widely used in agriculture, forestry, urban development planning and other fields.

城市作为人类发展的产物,随着人口的增加,科技进步,在不断扩大。城市占地面积的增加,意味着耕地面积的减少。世界各国纷纷采用城市垂直方向发展策略来缓解这一矛盾,因此不管在城市的商业中心,还是居民区都出现了大量高楼大厦。使得仅仅利用传统的空间二维城市发展密度指标对城市密度进行评价不在合理,必须考虑城市垂直方向的高度信息。这促使了联合利用多/高光谱图像提供的光谱和空间二维信息和激光雷达数据提供的垂直方向高度信息进行城市密度估计。As a product of human development, cities are constantly expanding with the increase of population and technological progress. The increase of urban area means the reduction of cultivated land. Countries around the world have adopted urban vertical development strategies to alleviate this contradiction. Therefore, a large number of high-rise buildings have appeared in both commercial centers and residential areas of the city. It is unreasonable to evaluate urban density only by using the traditional spatial two-dimensional urban development density index, and the height information in the vertical direction of the city must be considered. This motivates the joint utilization of spectral and spatial two-dimensional information provided by multi/hyperspectral imagery and vertical height information provided by lidar data for urban density estimation.

为了计算城市密度指标,必须首先获取地物覆盖主题图。如果是三维(3-D)指标还必须知道对应的高度信息。用于计算城市密度指标的感兴趣区域可以是一个像素、一个网格、指定半径的圆或者别的自定义形状。如果采用单个像素作为感兴趣区域,那么只能计算二维(2-D)城市密度指标,通过光谱解混技术在亚像素水平计算特定地物所占的比例。在以别的感兴趣区域作为指标计算对象时,研究者们提出了许多2-D指标来评价城市密度:如不透水表面面积、建筑物覆盖率、植被覆盖度等;相对于2-D指标,3-D指标很少,如植被体积、建筑物体积。In order to calculate the urban density index, the thematic map of ground object coverage must be obtained first. If it is a three-dimensional (3-D) indicator, the corresponding height information must also be known. The ROI used to calculate the urban density indicator can be a pixel, a grid, a circle with a specified radius, or another custom shape. If a single pixel is used as the region of interest, then only two-dimensional (2-D) urban density indicators can be calculated, and the proportion of specific ground features can be calculated at the sub-pixel level through spectral unmixing technology. When using other areas of interest as indicators to calculate objects, researchers have proposed many 2-D indicators to evaluate urban density: such as impermeable surface area, building coverage, vegetation coverage, etc.; compared to 2-D indicators , there are few 3-D indicators, such as vegetation volume, building volume.

目前在城市密度估计应用中存在利用的2-D指标单一、3-D指标少的问题,无法合理、全面评价城市密度。针对存在的这一问题,本发明联合利用多/高光谱图像和激光雷达数据计算能合理评价城市密度的3-D指标,并联合2-D和3-D指标对城市密度进行更合理评价。At present, in the application of urban density estimation, there are problems of single 2-D indicators and few 3-D indicators, which make it impossible to evaluate urban density reasonably and comprehensively. In view of this existing problem, the present invention jointly utilizes multi/hyperspectral images and lidar data to calculate 3-D indicators that can reasonably evaluate urban density, and combines 2-D and 3-D indicators to evaluate urban density more reasonably.

发明内容Contents of the invention

本发明为解决现有城市密度估计方法中,存在利用的2-D指标单一、3-D指标少,无法合理、全面评价城市密度的问题,而提出一种融合光谱图像和激光雷达数据进行城市密度估计方法。In order to solve the problem that the existing urban density estimation method has a single 2-D index and few 3-D indexes, and cannot reasonably and comprehensively evaluate the urban density, the present invention proposes a fusion spectral image and laser radar data for urban density estimation. Density Estimation Method.

一种融合光谱图像和激光雷达数据进行城市密度估计方法,按以下步骤进行:A method for estimating urban density by fusing spectral images and lidar data, the steps are as follows:

一、获取多/高光谱图像和激光雷达数据,分别对两种数据源进行预处理,利用激光雷达数据生成数字表面模型(DSM),并在两种数据源中选择控制点,对上述两种数据源进行配准;1. Obtain multi/hyperspectral images and lidar data, preprocess the two data sources respectively, use the lidar data to generate a digital surface model (DSM), and select control points in the two data sources, and the above two Data source for registration;

其中,对光谱图像预处理为辐射矫正和几何矫正;对激光雷达数据预处理为奇异点剔除和图像栅格化;Among them, the preprocessing of spectral images is radiation correction and geometric correction; the preprocessing of lidar data is singular point removal and image rasterization;

二、在多/高光谱图像上提取光谱信息和空间信息,在激光雷达数据上提取高度信息:2. Extract spectral information and spatial information on multi/hyperspectral images, and extract height information on lidar data:

其中光谱信息包括原始光谱波段、归一化植被指数和归一化建筑指数,空间信息包括通过使用均值、方差、形态学和Gabor空间滤波器生成的空间特征;在激光雷达数据上提取高度信息为归一化数字表面模型;Among them, the spectral information includes the original spectral band, normalized vegetation index and normalized building index, and the spatial information includes the spatial features generated by using the mean, variance, morphology and Gabor spatial filter; the height information extracted on the lidar data is Normalized digital surface model;

三、将提取得到的光谱信息、空间信息和高度信息输入到分类器中,得到分类主题图;3. Input the extracted spectral information, spatial information and height information into the classifier to obtain the classification theme map;

四、联合利用分类主题图和激光雷达提供的高度信息进行城市密度指标计算,最终生成城市密度主题图;4. Combine the classified theme map and the height information provided by the lidar to calculate the urban density index, and finally generate the city density theme map;

其中城市密度指标包括2-D和3-D指标,2-D指标为植被覆盖率和人造表面覆盖率;3-D指标为建筑物占地面积与使用面积之比和建筑物集成指数。Among them, urban density indicators include 2-D and 3-D indicators, 2-D indicators are vegetation coverage and artificial surface coverage; 3-D indicators are the ratio of building area to usable area and building integration index.

本发明包括以下有益效果:The present invention comprises following beneficial effect:

1、由于利用了激光雷达数据提供的高度信息,可以设计合理的3-D指标,克服的目前存在的3-D指标缺少的问题;1. Due to the use of the height information provided by the lidar data, a reasonable 3-D index can be designed to overcome the current lack of 3-D indexes;

2、进行城市密度评价时,本发明联合利用了2-D和3-D指标对城市密度进行更合理评价,克服了传统使用单一指标进行评价存在的不合理问题。2. When evaluating urban density, the present invention jointly utilizes 2-D and 3-D indicators to evaluate urban density more reasonably, and overcomes the unreasonable problem of traditional evaluation using a single index.

附图说明Description of drawings

图1为本发明所述方法流程示意图。Fig. 1 is a schematic flow chart of the method of the present invention.

具体实施方式detailed description

为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合图1和具体实施方式对本发明作进一步详细的说明。In order to make the above-mentioned purpose, features and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below in conjunction with FIG. 1 and specific embodiments.

具体实施方式一、本实施方式所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法,按以下步骤进行:Specific Embodiments 1. A method for estimating urban density by fusing spectral images and lidar data described in this embodiment is carried out in the following steps:

一、获取多/高光谱图像和激光雷达数据,分别对两种数据源进行预处理,利用激光雷达数据生成数字表面模型(DSM),并在两种数据源中选择控制点,对上述两种数据源进行配准;1. Obtain multi/hyperspectral images and lidar data, preprocess the two data sources respectively, use the lidar data to generate a digital surface model (DSM), and select control points in the two data sources, and the above two Data source for registration;

其中,对光谱图像预处理为辐射矫正和几何矫正;对激光雷达数据预处理为奇异点剔除和图像栅格化;Among them, the preprocessing of spectral images is radiation correction and geometric correction; the preprocessing of lidar data is singular point removal and image rasterization;

二、在多/高光谱图像上提取光谱信息和空间信息,在激光雷达数据上提取高度信息:2. Extract spectral information and spatial information on multi/hyperspectral images, and extract height information on lidar data:

其中光谱信息包括原始光谱波段、归一化植被指数和归一化建筑指数,空间信息包括通过使用均值、方差、形态学和Gabor空间滤波器生成的空间特征;在激光雷达数据上提取高度信息为归一化数字表面模型;Among them, the spectral information includes the original spectral band, normalized vegetation index and normalized building index, and the spatial information includes the spatial features generated by using the mean, variance, morphology and Gabor spatial filter; the height information extracted on the lidar data is Normalized digital surface model;

三、将提取得到的光谱信息、空间信息和高度信息输入到分类器中,得到分类主题图;3. Input the extracted spectral information, spatial information and height information into the classifier to obtain the classification theme map;

四、联合利用分类主题图和激光雷达提供的高度信息进行城市密度指标计算,最终生成城市密度主题图;4. Combine the classified theme map and the height information provided by the lidar to calculate the urban density index, and finally generate the city density theme map;

其中城市密度指标包括2-D和3-D指标,2-D指标为植被覆盖率和人造表面覆盖率;3-D指标为建筑物占地面积与使用面积之比和建筑物集成指数。Among them, urban density indicators include 2-D and 3-D indicators, 2-D indicators are vegetation coverage and artificial surface coverage; 3-D indicators are the ratio of building area to usable area and building integration index.

本实施方式包括以下有益效果:This embodiment includes the following beneficial effects:

1、由于利用了激光雷达数据提供的高度信息,可以设计合理的3-D指标,克服的目前存在的3-D指标缺少的问题;1. Due to the use of the height information provided by the lidar data, a reasonable 3-D index can be designed to overcome the current lack of 3-D indexes;

2、进行城市密度评价时,本实施方式联合利用了2-D和3-D指标对城市密度进行更合理评价,克服了传统使用单一指标进行评价存在的不合理问题。2. When evaluating the urban density, this embodiment uses 2-D and 3-D indicators to evaluate the urban density more reasonably, which overcomes the unreasonable problem of traditional evaluation using a single index.

具体实施方式二、本实施方式是对具体实施方式一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明,步骤一中所述的控制点为道路交叉点,或建筑物的拐点。Specific embodiment 2. This embodiment is a further description of the method for estimating urban density by fusing spectral images and lidar data described in specific embodiment 1. The control points described in step 1 are road intersections, or building turning point of things.

具体实施方式三、本实施方式是对具体实施方式一或二所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明,步骤一中所述的激光雷达数据奇异点是指由于空中飞行的鸟、漂浮的垃圾和风筝等物体造成的高度明显高于地物的点云,采用直方图统计方法进行激光雷达奇异点云数据剔除。Specific Embodiment Three. This embodiment is a further description of a method for estimating urban density by fusing spectral images and laser radar data described in specific embodiment one or two. The singular point of the laser radar data described in step one refers to Due to objects such as birds flying in the air, floating garbage, and kites, the height of which is significantly higher than the point cloud of ground objects, the histogram statistical method is used to remove the singular point cloud data of lidar.

具体实施方式四、本实施方式是对具体实施方式一至三之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明,步骤二中所述的归一化数字表面模型按如下步骤计算:Specific Embodiment 4. This embodiment is a further description of a method for estimating urban density by fusing spectral images and lidar data described in one of specific embodiments 1 to 3. The normalized digital surface model described in step 2 Calculate as follows:

1、基于移动二次曲面拟合的激光雷达图像滤波,将数字表面模型中的点集分为了地面点与非地面点两部分,具体实施步骤为:1. Based on the lidar image filtering based on moving quadratic surface fitting, the point set in the digital surface model is divided into two parts: ground points and non-ground points. The specific implementation steps are as follows:

a、选取合适的滤波窗口尺寸m×n,其中m和n的大小为102米数量级;a. Select an appropriate filter window size m×n, wherein the size of m and n is on the order of 10 2 meters;

b、在该窗口中选取高度最低的10个点,作为初始种子点,输入到地面点集P中,此10个点一定是地面上物体构成的点;b. In this window, select the 10 points with the lowest height as the initial seed points, and input them into the ground point set P. These 10 points must be points formed by objects on the ground;

c、利用点集P中的点进行在二次曲面拟合,拟合所涉及的函数方程为:c. Use the points in the point set P to fit the quadratic surface. The function equation involved in the fitting is:

ZZ ii == cc 00 ++ cc 11 xx ii ++ cc 22 ythe y ii ++ cc 33 xx ii ythe y ii ++ cc 44 xx ii 22 ++ cc 55 ythe y ii 22

其中,(xi,yi)为第i个点在图像中的坐标,Zi为该点对应的高度值;Among them, (x i , y i ) is the coordinate of the i-th point in the image, and Z i is the height value corresponding to the point;

依次将P集合内的点输入,得到一系列的方程组,在最小二乘准则下解出各系数c0,c1,...,c5,从而确定曲面方程;Input the points in the P set in turn to obtain a series of equations, and solve the coefficients c 0 , c 1 ,...,c 5 under the least squares criterion to determine the surface equation;

d、利用得到的曲面方程对其他点的高度进行预测,若预测值与实际值的差大于阈值T,则判定该点为非地面点;反之,则该点为地面点,然后将该点加入到地面点集P中,重新计算各系数,得到新曲面,如此重复直至所有点都判定完毕;d. Use the obtained surface equation to predict the height of other points. If the difference between the predicted value and the actual value is greater than the threshold T, the point is determined to be a non-ground point; otherwise, the point is a ground point, and then the point is added to In the ground point set P, recalculate each coefficient to obtain a new surface, and repeat until all points are determined;

e、移动窗口至图像其它位置,完成整幅图像的滤波;e. Move the window to other positions in the image to complete the filtering of the entire image;

2、归一化数字表面模型的生成:2. Generation of normalized digital surface model:

采用反距离加权插值的算法进行高度内插,插值时对距待插值点近的地面点给予较大的权值,首先以内插点为中心,确定适当数目N0的最近邻点作为源采样点,假设内插点为S0(x0,y0),采样点为Qi(xi,yi,zi),i=(1,2,...,N0),反距离加权平均插值的数学表达式如下:The inverse distance weighted interpolation algorithm is used for height interpolation. When interpolating, a larger weight is given to the ground point close to the point to be interpolated. First, the interpolation point is centered, and an appropriate number of N 0 nearest neighbor points is determined as the source sampling point. , assuming that the interpolation point is S 0 (x 0 ,y 0 ), the sampling point is Q i ( xi ,y i , zi ), i=(1,2,...,N 0 ), inverse distance weighted The mathematical expression for mean interpolation is as follows:

ZZ SS 00 == ΣΣ ii == 11 NN 00 λλ ii ZZ ii

其中为内插点的高度估计值;λi为第i个采样点的权值;Zi为第i个采样点的高度值,λi和Zi由以下公式求得: in is the estimated height of the interpolation point; λ i is the weight of the i-th sampling point; Z i is the height value of the i-th sampling point, and λ i and Z i are obtained by the following formula:

di为第i个采样点到内插点的距离:d i is the distance from the ith sampling point to the interpolation point:

dd ii == (( xx ii -- xx 00 )) 22 ++ (( ythe y ii -- ythe y 00 )) 22

其中p为正整数,当p取值为2时,插值的效果最好;Where p is a positive integer, when the value of p is 2, the interpolation effect is the best;

完成插值过程后,即生成数字地形模型(DTM),然后用DSM与DTM相减就是归一化数字表面模型。After the interpolation process is completed, a digital terrain model (DTM) is generated, and then the DSM is subtracted from the DTM to obtain a normalized digital surface model.

具体实施方式五、本实施方式是对具体实施方式一至四之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明,步骤三中所述的分类器为决策树、支持向量机或神经网络。Specific embodiment five. This embodiment is a further description of a method for estimating urban density by fusing spectral images and lidar data described in one of specific embodiments one to four. The classifier described in step three is a decision tree, Support vector machines or neural networks.

具体实施方式六、本实施方式是对具体实施方式一至五之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明,步骤四中所述的植被覆盖率按如下公式进行计算:Specific embodiment six, this embodiment is a further description of a method for estimating urban density by fusing spectral images and lidar data described in one of specific embodiments one to five, and the vegetation coverage described in step four is as follows Calculation:

VV Ff == AA tt ++ AA gg AA AA Oo II

其中,VF为植被覆盖率,At为感兴趣区域内被树木覆盖的区域面积,Ag为感兴趣区域内被草和灌木覆盖的区域面积,AAOI为感兴趣区域的面积;Wherein, VF is the vegetation coverage rate, A t is the area covered by trees in the region of interest, A g is the area covered by grass and shrubs in the region of interest, and A AOI is the area of the region of interest;

所述的人造表面覆盖率按如下公式进行计算:The artificial surface coverage is calculated according to the following formula:

AA SS CC == AA bb ++ AA ii AA AA Oo II

其中,ASC为人造表面覆盖率,Ab为感兴趣区域内被建筑物覆盖的区域面积,Ai为感兴趣区域内被道路、广场等人造地表所占面积;Among them, ASC is the artificial surface coverage rate, A b is the area covered by buildings in the area of interest, and A i is the area occupied by artificial surfaces such as roads and squares in the area of interest;

所述的建筑物占地面积与使用面积之比按如下公式计算:The ratio of the mentioned building area to the usable area is calculated according to the following formula:

ii Ff AA RR == AA aa bb AA ff ll

其中,Aab为建筑物占地面积,Afl为建筑物使用面积,这里假设每层楼高3米,楼层数目N按公式计算:计算,其中HB指楼高,表示向下取整;建筑物使用面积按公式计算:其中为每层楼的面积;Among them, A ab is the floor area of the building, and A fl is the usable area of the building. Here, it is assumed that each floor is 3 meters high, and the number of floors N is calculated according to the formula: Calculate, where H B refers to the height of the building, Indicates rounding down; the usable area of the building is calculated according to the formula: in is the area of each floor;

所述的建筑物集成指数按如下公式计算:The building integration index is calculated according to the following formula:

BB AA == AA bb AA AA Oo II Mm ee dd ii aa nno (( DD. bb )) -- Mm ee dd ii aa nno (( ii Ff AA RR )) NN bb

其中,Db为感兴趣区域内每两栋建筑物中心之间的距离,Median代表取中值,Nb为感兴趣区域内的建筑物数目。Among them, D b is the distance between the centers of every two buildings in the area of interest, Median represents the median value, and N b is the number of buildings in the area of interest.

具体实施方式七、本实施方式是对具体实施方式一至六之一所述的一种融合光谱图像和激光雷达数据进行城市密度估计方法的进一步说明,步骤四中所述生成城市密度主题图的具体内容为:每一个像素为计算中心,250米为半径的圆形区域为对应中心像素的感兴趣区域进行城市密度(UD)估计,按如下公式进行计算:Specific Embodiment 7. This embodiment is a further description of a method for estimating urban density by fusing spectral images and lidar data described in one of specific embodiments 1 to 6. The specific method for generating the urban density theme map described in step 4 The content is: each pixel is the calculation center, and a circular area with a radius of 250 meters is used to estimate the urban density (UD) of the region of interest corresponding to the central pixel, and the calculation is performed according to the following formula:

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

其中,将植被覆盖率、人造表面覆盖率、建筑物占地面积与使用面积之比和建筑物集成指数代入城市密度估计公式计算前都归一化为0到1,最终得到的UD值区间为[-2,2],值越大代表城市密度越大。Among them, the vegetation coverage rate, artificial surface coverage rate, ratio of building area to usable area, and building integration index are all normalized to 0 to 1 before being substituted into the urban density estimation formula, and the final UD value range is [-2,2], the larger the value, the greater the urban density.

Claims (7)

1. a fusion spectrum picture and laser radar data carry out city density method of estimation, it is characterised in that according to the following steps Carry out:
One, obtain many/high spectrum image and laser radar data, respectively two kinds of data sources are carried out pretreatment, utilize laser thunder Reach data genaration digital surface model DSM, and in two kinds of data sources, select control point, above two data source is joined Accurate;
Wherein, spectrum picture pretreatment is corrected and Geometry rectification for radiation;It is that singular point picks to laser radar data pretreatment Remove and image tiles;
Two, on many/high spectrum image, spectral information and spatial information are extracted, extraction elevation information on laser radar data:
Wherein spectral information includes original spectrum wave band, normalized differential vegetation index and normalization building index, and spatial information includes By the space characteristics using average, variance, morphology and Gabor spatial filter to generate;Laser radar data extracts Elevation information is normalization digital surface model;
Three, spectral information extraction obtained and spatial information are input in grader, obtain classification scheme figure;
Four, combine the elevation information utilizing classification scheme figure and laser radar to provide and carry out city density index calculating, the most throughout one's life Become city density thematic map;
Wherein city density index includes that 2-D and 3-D index, 2-D index are vegetation coverage and artificial surfaces coverage rate;3-D Index is building floor space and the ratio of usable floor area and building collection exponentially.
2. a kind of fusion spectrum picture as claimed in claim 1 and laser radar data carry out city density method of estimation, its It is characterised by that the control point described in step one is road junction, or the flex point of building.
3. a kind of fusion spectrum picture as claimed in claim 1 or 2 and laser radar data carry out city density method of estimation, It is characterized in that the laser radar data singular point described in step one refers to the bird due to airflight, floating rubbish and wind The height that the objects such as zither cause, apparently higher than the some cloud of atural object, uses statistics with histogram method to carry out laser radar singular point cloud number According to rejecting.
4. a kind of fusion spectrum picture as claimed in claim 3 and laser radar data carry out city density method of estimation, its It is characterised by that the normalization digital surface model described in step 2 calculates as follows:
4.1, lidar image filtering based on mobile Quadratic Surface Fitting, has been divided into ground by the point set in digital surface model Cake and non-ground points two parts, being embodied as step is:
A, choosing suitable filter window size m × n, wherein the size of m and n is 102The rice order of magnitude;
B, in this window, choose minimum 10 point, as initial seed point, be input in the point set P of ground, these 10 Point must be the point that on ground, object is constituted;
C, utilizing the point in point set P to carry out at Quadratic Surface Fitting, the functional equation involved by matching is:
Zi=c0+c1xi+c2yi+c3xiyi+c4xi 2+c5yi 2
Wherein, (xi,yi) it is i-th point coordinate in the picture, ZiFor the height value that this point is corresponding;
Point input in being gathered by P successively, obtains a series of equation group, solves each coefficient c under criterion of least squares0, c1,...,c5, so that it is determined that surface equation;
Other height put are predicted by the surface equation that d, utilization obtain, if the difference of predictive value and actual value is more than threshold value T, Then judge that this point is as non-ground points;Otherwise, then this point is ground point, is then joined in the point set P of ground by this point, recalculates Each coefficient, obtains new curved surface, be so repeated up to the most all judge complete;
E, moving window, to other position of image, complete the filtering of entire image;
4.2, the generation of normalization digital surface model:
The algorithm using inverse distance weighted interpolation carries out height interpolation, gives bigger to the ground point near away from interpolation point during interpolation Weights, first centered by interpolated point, determine proper number N0Nearest neighbor point as source sampling point, it is assumed that interpolated point is S0 (x0,y0), sampled point is Qi(xi,yi,zi), i=(1,2 ..., N0), the mathematic(al) representation of inverse distance-weighting average interpolation is such as Under:
Z S 0 = Σ i = 1 N 0 λ i Z i
WhereinHeight Estimation value for interpolated point;λiWeights for ith sample point;ZiFor the height value of ith sample point, λiAnd ZiTried to achieve by below equation:
diDistance for ith sample point to interpolated point:
d i = ( x i - x 0 ) 2 + ( y i - y 0 ) 2
Wherein p is positive integer;
After completing Interpolation Process, i.e. generating digital terrain model DTM, then subtracting each other with DSM Yu DTM is exactly normalization digital surface Model.
5. a kind of fusion spectrum picture as claimed in claim 4 and laser radar data carry out city density method of estimation, its It is characterised by that the grader described in step 3 is decision tree, support vector machine or neutral net.
6. a kind of fusion spectrum picture as claimed in claim 5 and laser radar data carry out city density method of estimation, its It is characterised by that the vegetation coverage described in step 4 calculates as follows:
V F = A t + A g A A O I
Wherein, VF is vegetation coverage, AtFor the region area covered by trees in area-of-interest, AgFor in area-of-interest The region area covered by grass and shrub, AAOIArea for area-of-interest;
Described artificial surfaces coverage rate calculates as follows:
A S C = A b + A i A A O I
Wherein, ASC is artificial surfaces coverage rate, AbFor the region area covered by building in area-of-interest, AiFor interested By road, artificial earth's surface, square occupied area in region;
Described building floor space is calculated as follows with the ratio of usable floor area:
i F A R = A a b A f l
Wherein, AabFor building floor space, AflFor building usable floor area, it is assumed here that every floor is high 3 meters, floor number N Calculate by formula:Calculate, wherein HBRefer to that building is high,Represent and round downwards;Building usable floor area is pressed formula and is calculated:WhereinArea for every floor;
Described building collection exponentially is calculated as follows:
B A = A b A A O I M e d i a n ( D b ) - M e d i a n ( i F A R ) N b
Wherein, DbFor the distance between two solitary building centers every in area-of-interest, Median represents and takes intermediate value, NbEmerging for sense Building number in interest region.
7. a kind of fusion spectrum picture as claimed in claim 6 and laser radar data carry out city density method of estimation, its It is characterised by that the particular content generating city density thematic map described in step 4 is: each pixel is calculating center, 250 meters The area-of-interest that border circular areas is corresponding center pixel for radius carries out city density UD estimation, counts as follows Calculate:
UD=(BA+ASC)-(VF+iFAR)
Wherein, by integrated with the ratio of usable floor area and building to vegetation coverage, artificial surfaces coverage rate, building floor space Index substitutes into before city density estimation formulas calculates and is all normalized to 0 to 1, and the UD value finally given is interval is [-2,2], and value is more Represent greatly city density the biggest.
CN201610356309.3A 2016-05-25 2016-05-25 A kind of fusion spectrum picture and laser radar data carry out city density estimation method Active CN106056591B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610356309.3A CN106056591B (en) 2016-05-25 2016-05-25 A kind of fusion spectrum picture and laser radar data carry out city density estimation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610356309.3A CN106056591B (en) 2016-05-25 2016-05-25 A kind of fusion spectrum picture and laser radar data carry out city density estimation method

Publications (2)

Publication Number Publication Date
CN106056591A true CN106056591A (en) 2016-10-26
CN106056591B CN106056591B (en) 2019-01-18

Family

ID=57175271

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610356309.3A Active CN106056591B (en) 2016-05-25 2016-05-25 A kind of fusion spectrum picture and laser radar data carry out city density estimation method

Country Status (1)

Country Link
CN (1) CN106056591B (en)

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106529484A (en) * 2016-11-16 2017-03-22 哈尔滨工业大学 Combined spectrum and laser radar data classification method based on class-fixed multinucleated learning
CN107092020A (en) * 2017-04-19 2017-08-25 北京大学 Merge the surface evenness monitoring method of unmanned plane LiDAR and high score image
CN107392130A (en) * 2017-07-13 2017-11-24 西安电子科技大学 Classification of Multispectral Images method based on threshold adaptive and convolutional neural networks
CN107657206A (en) * 2016-12-23 2018-02-02 航天星图科技(北京)有限公司 A kind of method based on remote sensing technology estimation forest coverage rate
CN108020324A (en) * 2017-11-07 2018-05-11 湖北理工学院 The filter detecting method of single beam laser superimposed pulse signal
CN108334819A (en) * 2017-01-17 2018-07-27 德尔福技术有限公司 Ground classifier system for automated vehicle
CN108846352A (en) * 2018-06-08 2018-11-20 广东电网有限责任公司 A kind of vegetation classification and recognition methods
CN110427857A (en) * 2019-07-26 2019-11-08 国网湖北省电力有限公司检修公司 A kind of transmission line of electricity geological disasters analysis method based on Remote Sensing Data Fusion Algorithm
CN111257882A (en) * 2020-03-19 2020-06-09 北京三快在线科技有限公司 Data fusion method and device, unmanned equipment and readable storage medium
CN111638185A (en) * 2020-05-09 2020-09-08 哈尔滨工业大学 Remote sensing detection method based on unmanned aerial vehicle platform
CN111679287A (en) * 2020-06-05 2020-09-18 中国科学院空天信息创新研究院 An Active Video Stereo Hyperspectral Imaging Method
CN111861170A (en) * 2020-07-07 2020-10-30 河南大学 A carbon emission spatial mapping method, density spatial distribution determination method and device
CN112130169A (en) * 2020-09-23 2020-12-25 广东工业大学 A point cloud-level fusion method of lidar data and hyperspectral images
CN112419402A (en) * 2020-11-27 2021-02-26 广东电网有限责任公司肇庆供电局 Positioning method and system based on multispectral image and laser point cloud
CN112859108A (en) * 2021-01-28 2021-05-28 中国科学院南京土壤研究所 Method for extracting under-forest vegetation coverage under complex terrain condition by using ground laser radar data
CN115906476A (en) * 2022-11-18 2023-04-04 国网湖北省电力有限公司经济技术研究院 A Calculation Method of Mountain Photovoltaic Power Generation

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5659391A (en) * 1996-01-26 1997-08-19 The United States Of America As Represented By The Secretary Of The Army Earth monitoring satellite system with combined infrared interferometry and photopolarimetry for chemical and biological sensing
CA2547359A1 (en) * 2003-11-26 2005-06-16 Florida Environmental Research Institute, Inc. Spectral imaging system
CN101866385A (en) * 2010-06-24 2010-10-20 浙江大学 A simulation and optimization method for surface temperature of target plots
CN101894210A (en) * 2010-06-24 2010-11-24 浙江大学 A Method of Acquiring Land-block Energy Saving Energy Based on 3D Scene Image

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5659391A (en) * 1996-01-26 1997-08-19 The United States Of America As Represented By The Secretary Of The Army Earth monitoring satellite system with combined infrared interferometry and photopolarimetry for chemical and biological sensing
CA2547359A1 (en) * 2003-11-26 2005-06-16 Florida Environmental Research Institute, Inc. Spectral imaging system
CN101866385A (en) * 2010-06-24 2010-10-20 浙江大学 A simulation and optimization method for surface temperature of target plots
CN101894210A (en) * 2010-06-24 2010-11-24 浙江大学 A Method of Acquiring Land-block Energy Saving Energy Based on 3D Scene Image

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
EMMANUEL BRATSOLIS等: "Automated Building Block Extraction and Building Density Classification Using Aerial Imagery and LiDAR Data", 《JOURNAL OF EARTH SCIENCE AND ENGINEERING》 *
冯凯: "基于多核学习的多/高光谱图像与激光雷达数据联合分类研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》 *
满其霞: "激光雷达和高光谱数据融合的城市土地利用分类方法研究", 《中国博士学位论文全文数据库基础科学辑》 *

Cited By (27)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106529484A (en) * 2016-11-16 2017-03-22 哈尔滨工业大学 Combined spectrum and laser radar data classification method based on class-fixed multinucleated learning
CN107657206A (en) * 2016-12-23 2018-02-02 航天星图科技(北京)有限公司 A kind of method based on remote sensing technology estimation forest coverage rate
CN108334819A (en) * 2017-01-17 2018-07-27 德尔福技术有限公司 Ground classifier system for automated vehicle
CN108334819B (en) * 2017-01-17 2021-12-03 安波福技术有限公司 Ground classifier system for automated vehicles
CN107092020A (en) * 2017-04-19 2017-08-25 北京大学 Merge the surface evenness monitoring method of unmanned plane LiDAR and high score image
CN107092020B (en) * 2017-04-19 2019-09-13 北京大学 Road roughness monitoring method based on UAV LiDAR and high-resolution images
CN107392130B (en) * 2017-07-13 2020-12-08 西安电子科技大学 Multispectral Image Classification Method Based on Threshold Adaptive and Convolutional Neural Networks
CN107392130A (en) * 2017-07-13 2017-11-24 西安电子科技大学 Classification of Multispectral Images method based on threshold adaptive and convolutional neural networks
CN108020324A (en) * 2017-11-07 2018-05-11 湖北理工学院 The filter detecting method of single beam laser superimposed pulse signal
CN108846352A (en) * 2018-06-08 2018-11-20 广东电网有限责任公司 A kind of vegetation classification and recognition methods
CN108846352B (en) * 2018-06-08 2020-07-14 广东电网有限责任公司 Vegetation classification and identification method
CN110427857A (en) * 2019-07-26 2019-11-08 国网湖北省电力有限公司检修公司 A kind of transmission line of electricity geological disasters analysis method based on Remote Sensing Data Fusion Algorithm
CN111257882B (en) * 2020-03-19 2021-11-19 北京三快在线科技有限公司 Data fusion method and device, unmanned equipment and readable storage medium
CN111257882A (en) * 2020-03-19 2020-06-09 北京三快在线科技有限公司 Data fusion method and device, unmanned equipment and readable storage medium
CN111638185A (en) * 2020-05-09 2020-09-08 哈尔滨工业大学 Remote sensing detection method based on unmanned aerial vehicle platform
CN111638185B (en) * 2020-05-09 2022-05-17 哈尔滨工业大学 Remote sensing detection method based on unmanned aerial vehicle platform
CN111679287B (en) * 2020-06-05 2023-01-17 中国科学院空天信息创新研究院 An Active Video Stereo Hyperspectral Imaging Method
CN111679287A (en) * 2020-06-05 2020-09-18 中国科学院空天信息创新研究院 An Active Video Stereo Hyperspectral Imaging Method
CN111861170A (en) * 2020-07-07 2020-10-30 河南大学 A carbon emission spatial mapping method, density spatial distribution determination method and device
CN111861170B (en) * 2020-07-07 2023-08-04 河南大学 A carbon emission spatial mapping method, density spatial distribution determination method and device
CN112130169A (en) * 2020-09-23 2020-12-25 广东工业大学 A point cloud-level fusion method of lidar data and hyperspectral images
CN112130169B (en) * 2020-09-23 2022-09-16 广东工业大学 A point cloud-level fusion method of lidar data and hyperspectral images
CN112419402A (en) * 2020-11-27 2021-02-26 广东电网有限责任公司肇庆供电局 Positioning method and system based on multispectral image and laser point cloud
CN112859108A (en) * 2021-01-28 2021-05-28 中国科学院南京土壤研究所 Method for extracting under-forest vegetation coverage under complex terrain condition by using ground laser radar data
CN112859108B (en) * 2021-01-28 2024-03-22 中国科学院南京土壤研究所 Method for extracting vegetation coverage under forests under complex terrain condition by using ground laser radar data
CN115906476A (en) * 2022-11-18 2023-04-04 国网湖北省电力有限公司经济技术研究院 A Calculation Method of Mountain Photovoltaic Power Generation
CN115906476B (en) * 2022-11-18 2023-09-01 国网湖北省电力有限公司经济技术研究院 A Calculation Method of Mountain Photovoltaic Power Generation

Also Published As

Publication number Publication date
CN106056591B (en) 2019-01-18

Similar Documents

Publication Publication Date Title
CN106056591A (en) Method for estimating urban density through fusion of optical spectrum image and laser radar data
CN108363983B (en) An urban vegetation classification method based on UAV images and reconstructed point clouds
CN102103202B (en) Semi-supervised classification method for airborne laser radar data fusing images
Asadi et al. Simulation of green roofs and their potential mitigating effects on the urban heat island using an artificial neural network: A case study in Austin, Texas
CN102054274B (en) Method for full automatic extraction of water remote sensing information in coastal zone
CN107609526A (en) Rule-based fine dimension city impervious surface rapid extracting method
CN103824077B (en) Urban impervious layer rate information extraction method based on multi-source remote sensing data
CN106780089B (en) Permanent basic farmland planning method based on neural network cellular automaton model
CN110135354B (en) Change detection method based on live-action three-dimensional model
CN104318270A (en) Land cover classification method based on MODIS time series data
CN109508585A (en) A method of urban function region is extracted based on POI and high-resolution remote sensing image
CN105118090A (en) Adaptive point-cloud filtering method for complex terrain structure
CN110222586A (en) A kind of calculating of depth of building and the method for building up of urban morphology parameter database
CN104463164A (en) Tree canopy structure information extraction method based on rib method and crown height ratio
CN110990511B (en) A method for classifying local climate zones considering two-dimensional and three-dimensional forms of cities
CN108052876A (en) Regional development appraisal procedure and device based on image identification
CN104517024A (en) Urban green degree space evaluation modeling method based on construction dimension
CN116561509A (en) Urban vegetation overground biomass accurate inversion method and system considering vegetation types
CN113468982B (en) Method, device and storage medium for classifying urban functional areas
CN106096129B (en) A kind of foot of the hill water surface scale analysis method calculated based on mountainous region charge for remittance
CN106844979B (en) Interactive combined design system and method for urban virtual park
CN104463971A (en) Method for creating green-degree spatial arrangement curve used for evaluating urban landscaping three-dimensional layout
Nelson et al. Spatial statistical techniques for aggregating point objects extracted from high spatial resolution remotely sensed imagery
CN116152668B (en) A Method of Obtaining Block-scale LCZ Based on Artificial Intelligence
Alphan et al. Mapping availability of sea view for potential building development areas

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20210118

Address after: Building 9, accelerator, 14955 Zhongyuan Avenue, Songbei District, Harbin City, Heilongjiang Province

Patentee after: INDUSTRIAL TECHNOLOGY Research Institute OF HEILONGJIANG PROVINCE

Address before: 150001 No. 92 West straight street, Nangang District, Heilongjiang, Harbin

Patentee before: HARBIN INSTITUTE OF TECHNOLOGY

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230221

Address after: 150027 Room 412, Unit 1, No. 14955, Zhongyuan Avenue, Building 9, Innovation and Entrepreneurship Plaza, Science and Technology Innovation City, Harbin Hi tech Industrial Development Zone, Heilongjiang Province

Patentee after: Heilongjiang Industrial Technology Research Institute Asset Management Co.,Ltd.

Address before: Building 9, accelerator, 14955 Zhongyuan Avenue, Songbei District, Harbin City, Heilongjiang Province

Patentee before: INDUSTRIAL TECHNOLOGY Research Institute OF HEILONGJIANG PROVINCE

TR01 Transfer of patent right