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

Abstract

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.

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.

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

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:

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:

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;

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. 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. 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

Fig. 1 is a schematic flow chart of the method of the present invention.

detailed description

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:

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:

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;

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. 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. 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. 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. Select an appropriate filter window size m×n, wherein the size of m and n is on the order of 10 ² meters;

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. Use the points in the point set P to fit the quadratic surface. The function equation involved in the fitting is:

${\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 4</span>}}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 5</span>}}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

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;

Input the points in the P set in turn to obtain a series of equations, and solve the coefficients c ₀ , c ₁ ,...,c ₅ under the least squares criterion to determine the surface equation;

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. Move the window to other positions in the image to complete the filtering of the entire image;

2. Generation of normalized digital surface model:

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 ₀ nearest neighbor points is determined as the source sampling point. , assuming that the interpolation point is S ₀ (x ₀ ,y ₀ ), the sampling point is Q _i ( _xi ,y _i , _zi ), i=(1,2,...,N ₀ ), inverse distance weighted The mathematical expression for mean interpolation is as follows:

${\mathrm{<span\; class="notranslate">\; Z</span>}}_{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}$

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:

d _i is the distance from the ith sampling point to the interpolation point:

${\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\sqrt{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}$

Where p is a positive integer, when the value of p is 2, the interpolation effect is the best;

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:

$\mathrm{<span\; class="notranslate">\; V</span>}\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; g</span>}}}{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; I</span>}}}$

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:

$\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; S</span>}\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; b</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; I</span>}}}$

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:

$\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; f</span>}\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; R</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; b</span>}}}{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}\mathrm{<span\; class="notranslate">\; l</span>}}}$

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:

$\mathrm{<span\; class="notranslate">\; B</span>}\mathrm{<span\; class="notranslate">\; A</span>}<span\; class="notranslate">\; =</span>\frac{\frac{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; b</span>}}}{{\mathrm{<span\; class="notranslate">\; A</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; I</span>}}}}{\mathrm{<span\; class="notranslate">\; m</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; b</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; m</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; f</span>}\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; R</span>}<span\; class="notranslate">\; )</span>}{{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; b</span>}}}$

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.

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)

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 10²The 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:

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

Wherein, (x_i,y_i) it is i-th point coordinate in the picture, Z_iFor 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 squares₀, c₁,...,c₅, 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 N₀Nearest neighbor point as source sampling point, it is assumed that interpolated point is S₀ (x₀,y₀), sampled point is Q_i(x_i,y_i,z_i), i=(1,2 ..., N₀), the mathematic(al) representation of inverse distance-weighting average interpolation is such as Under:

$Z_{S_0}=\underset{i=1}{\overset{N_0}{\&Sigma;}}{\mathrm{\&lambda;}}_i{Z}_i$

WhereinHeight Estimation value for interpolated point；λ_iWeights for ith sample point；Z_iFor the height value of ith sample point, λ_iAnd Z_iTried to achieve by below equation:

d_iDistance for ith sample point to interpolated point:

$d_i=\sqrt{{(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:

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

Wherein, VF is vegetation coverage, A_tFor the region area covered by trees in area-of-interest, A_gFor in area-of-interest The region area covered by grass and shrub, A_AOIArea for area-of-interest；

Described artificial surfaces coverage rate calculates as follows:

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

Wherein, ASC is artificial surfaces coverage rate, A_bFor the region area covered by building in area-of-interest, A_iFor 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:

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

Wherein, A_abFor building floor space, A_flFor building usable floor area, it is assumed here that every floor is high 3 meters, floor number N Calculate by formula:Calculate, wherein H_BRefer 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:

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

Wherein, D_bFor the distance between two solitary building centers every in area-of-interest, Median represents and takes intermediate value, N_bEmerging 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.