CN104050473A - Road data extraction method based on rectangular neighborhood analysis - Google Patents

Abstract

The invention discloses a road data extraction method based on rectangular neighborhood analysis. The method comprises the steps of data pre-processing, road candidate point screening, road center line extracting, flat terrain point set removing and road connecting. The road data extraction method based on the rectangular neighborhood analysis mainly aims to carry out road data extraction according to the difference between the road topological structure characteristics and the ambient environment under the condition that reflection strength is unavailable, or reflection strength cannot reflect ground feature attributes well, or aerial images are unavailable.

Description

A kind of road data extracting method of analyzing based on rectangular neighborhood

Technical field

The present invention relates to a kind of road data extracting method of analyzing based on rectangular neighborhood.

Background technology

Buildings and road extraction are the Focal point and difficult points that airborne LiDAR (Light Detection And Ranging) data characteristics is extracted always, aspect road extraction, mainly more around downtown roads net extraction research, less for method for extracting roads research in region, field.Wherein airborne LiDAR is that the range finding of a kind of set laser, GPS (GPS) and three kinds of technology of inertial navigation system (INS) are in the system of one, and for obtaining data and generating accurate three-dimensional landform (DEM); Airborne LiDAR (or claim cloud data) is the body surface point data that gathered by three-dimensional laser scanner, each some three-dimensional coordinate that has living space, have also with information such as color and inverse intensity; Point cloud filtering refers to the non-ground points such as the noise removed in some cloud, trees, house; Point cloud tissue refers to a cloud is managed according to certain data layout, is convenient to inquiry and the demonstration of data; Road extraction refers to automatically extracts road waypoint from cloud data.Existing method for extracting roads mainly contains method for extracting roads, the method for extracting roads based on airborne LiDAR data and remote sensing image etc. based on the strong information of reflection, and the conventional road data extracting method based on a cloud, reflection strength and image is as follows:

(1) the structured road extraction method of Vehicle-borne Laser Scanning data

Vehicle-borne Laser Scanning refers to and laser scanner is erected to the enterprising Mobile state vehicle of the mobile platforms such as automobile, be that with airborne laser scanning difference the scope of airborne laser scanning is larger, speed is faster, and the data that Vehicle-borne Laser Scanning obtains are all along road, it is comparatively easy therefore to extract.

Structured road refers to downtown roads or highway etc., and there is the feature that exceeds waypoint the both sides of these roads, by the extraction that sweep trace extracts, the step such as the filtering of ground point cloud, road boundary are extracted automatically, road boundary tracking completes road.Concrete grammar refers to " mapping journal " 42 2 phases of volume in 2013 " the structured road extraction methods of Vehicle-borne Laser Scanning data ".

(2) road extraction based on LiDAR echo information

First adopt a kind of method (as: successively encrypting TIN method) to extract or generating digital relief block (DTM), from DTM, extract road information according to a cloud echo information, finally delete noise spot wherein by the filtering algorithm of search isolated point.Concrete grammar refers to " mapping science " 36 2 phases of volume " based on the road extraction of LiDAR echo information " in 2011.

(3) road extraction based on airborne LiDAR roughness index and echo strength

First, by the derivative normalization digital surface model of cloud data and roughness index, then the multi-source data after registration is carried out to multi-resolution segmentation, and then use the features such as roughness index and echo strength, the road description factor to classify, finally, remove road noise, and obtain road frame net accurately.Refer to " surveying and drawing scientific and technical journal " 30 1 phases of volume in 2013 " based on the road extraction of airborne LiDAR roughness index and echo strength ".

(4) urban road network based on airborne LiDAR and high-resolution remote sensing image extracts

First both being carried out to accuracy registration then utilizes the method for pseudo-road information removal respectively the removal such as vegetation information and building information to be obtained to basic road profile, recycling form thinning algorithm extracts the center line of road, finally under ArcGIS and Matlab programmed environment, has realized improved road Pruning Algorithm (IRT) and has utilized this algorithm to carry out road to prune and obtained level and smooth and coherent urban road network.Refer to " remote sensing technology and application " 28 4 phases of volume in 2013 " urban road network based on airborne LiDAR and high-resolution remote sensing image extracts ".

But the reflectivity of return laser beam intensity and ground object target is proportional, be that the larger laser intensity value of reflectivity is higher, but the echo strength of laser radar system record is the value through overcorrect not yet, can be subject to the impact of the many factors such as atmospheric attenuation, laser incident angle, and the correction of return laser beam intensity so far or a more difficult problem, therefore return laser beam intensity is difficult to reflect truly the reflectivity information of atural object.The extraction effect of the road extraction method based on remote sensing image depends primarily on the readability of image, but the remote sensing image gathering in the situation that cloudy or visibility is inadequate often can not react atural object attribute well, therefore also difficulty comparatively of classification.

Summary of the invention

The present invention seeks to: a kind of road data extracting method of analyzing based on rectangular neighborhood is provided, do not adopting under the condition of reflection strength and remote sensing image, by analyzing the architectural feature of territory, forest land, mound road, adopt rectangular neighborhood judgement to carry out the extraction of road.

Technical scheme of the present invention is: a kind of road data extracting method of analyzing based on rectangular neighborhood, and it comprises:

The first step: data pre-service, calculation level cloud spacing, sets up regular grid tissue, then adopts the quick mathematical morphology filter method of rule-based grid to carry out filtering, rejects non-ground point;

Second step: screening road candidate point, tentatively determine road point set, dwindle determination range, to the point within the scope of its neighborhood of each point search, and compare its height value, if all in threshold range, as road candidate point;

The 3rd step: whether extract road axis, reject smooth region point set, to each road candidate point, calculate the interior elevation variance yields of rectangular extent of different orientations, be road center point according to the current point of the diversity judgement between variance yields;

The 4th step: road connects, connects the road being interrupted, level and smooth.

On the basis of technique scheme, further comprise following attached technical scheme:

The described first step comprises:

(1) first read cloud data scope (xmin, xmax, ymin, ymax), the quantity pt_Num of statistics point;

(2) calculation level cloud spacing:

$\mathrm{pt}\_\mathrm{dis}=\frac{(x\max -x\min )*(y\max -y\min )}{\sqrt{\mathrm{pt}\_\mathrm{num}}}$

(3) according to a cloud spacing pt_dis and calculating needs, set up a regular grid that covers some cloud scope (xmin, xmax, ymin, ymax);

(4) according to point coordinate by airborne LiDAR data allocations in each grid,

$\left\{\begin{array}{c}\mathrm{XAxisNum}=\mathrm{int}((x_{\max }-x_{\min })/d)+1\\ \mathrm{YAxisNum}=\mathrm{int}((y_{\max }-y_{\min })/d)+1\end{array}\right.$ (formula 2)

GridNum=XAxisNum*YAxisNum (formula 3)

XAxisNum, YAxisNum are the grid number of x, y direction, and GridNum is grid sum;

Any point p in airborne LiDAR data _i(x _i, y _i) grid at place is:

$\left\{\begin{array}{c}{X}_i=\mathrm{int}((x_i-x_{\min })/d)\\ Y_i=\mathrm{int}((y_i-y_{\min })/d)\\ G_i=\mathrm{int}((x_i-x_{\min })/d)+X_i\×\mathrm{int}((y_i-y_{\min })/d)\end{array}\right.$ (formula 4)

Wherein p _ifor p _ithe grid number of the directions X grid number at place, the Y-direction grid number that Yi is Pi place, one-dimension array record that Gi is Pi place;

Airborne LiDAR data are carried out to the filtering of non-ground, retain and comprise waypoint in interior ground point data, mathematical morphology filter taking regular grid as unit, grid Corrosion results is the minimum height value in w × w grid, expansion results is maximum elevation value in w × w grid, wherein w represents the size of window, and carries out the computing of formula 5:

$\left\{\begin{array}{c}\mathrm{Grid}{\left(\mathrm{A\ΘB}\right)}_i=\mathrm{Grid}\left(z_i^{\′}\right)=\underset{{\mathrm{Grid}}_i\⋐w}{\min }\left({\mathrm{Grid}}_i\left(z_{\min }\right)\right)\\ \mathrm{Grid}{(A\⊕B)}_i=\mathrm{Grid}\left(z_i^{\″}\right)=\underset{{\mathrm{Grid}}_i\⋐w}{\max }\left({\mathrm{Grid}}_i\left(z_i^{\′}\right)\right)\end{array}\right.$ (formula 5)

Wherein A represents goal set, B representative structure element, and A Θ B represents corrosion, represent to expand, z in formula _i' and z _iheight value after the corrosion of ' ' be respectively grid cell and expansion, Grid _i(z _min) be elevation minimum value in grid.

The described first step further comprises:

(1) by mesh spacing computing grid quantity, and be grid storage allocation space, point is stored in corresponding grid cell by planimetric coordinates one by one, simultaneously with this grid cell in the minimum height value of having preserved compare, the minimum height value of record after relatively;

(2) travel through each grid, minimum height value in w × w grid relatively around, first formula of employing formula 5, getting minimum value is the height value after current grid corrosion, and this grid maximum elevation value is made as to the height value after corrosion;

(3) travel through each grid, w × w grid internal corrosion height value relatively around, second formula of employing formula 5, getting maximal value is the height value after current grid expands, and minimum this grid height value is made as to the height value after expansion;

(4) travel through each grid, the height value after asking for each point height value in grid and expanding poor, in the time that its absolute value is less than setting threshold, this point is ground point, otherwise is non-ground point.

Described second step comprises:

(1) for each some p _i(x _i, y _i) calculate its place grid by formula 4;

(2) calculate search radius, $r=\frac{\mathrm{road}\_\mathrm{width}}{d};$

(3) p that sets up an office _i(x _i, y _i) rectangular neighborhood in point set be

P _{i_near}={ p _j(x _j, y _j) | x _j∈ (x _i-r, x _i+ r), y _j∈ (y _i-r, y _i+ r) } (wherein j is period or perhaps the subscript of point, and distinguishes with current some pi), thus centered by current some pi, search is the point set P in grid around _{i_near};

(4) calculate P _{i_near}in each point and p _i(x _i, y _i) difference of elevation absolute value sum Σ H _iif, Σ H _ibe less than threshold value, current point is road candidate point;

$\mathrm{\Σ}{H}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}|z_j-z_i|$ (formula 6)

Zj is P _{i_near}in j point height value, zi is p _i(x _i, y _i) height value, n is P _{i_near}always count;

(5) repeat above step, until complete to judgement a little, and be P by the road candidate point set of reservation _c.

Described the 3rd step comprises:

(1) to each road candidate point p _i, calculate rectangular extent point p according to formula 7 _{i, r1}, p _{i, r2}, p _{i, r3}, p _{i, r4}, i.e. four of rectangle summits, p _{i, r1}overlap with current point, all the other press counterclockwise arrangement, and 360 degree are divided into n interval by the anglec of rotation that wherein d θ is rectangle;

$\begin{array}{c}\left\{\begin{array}{c}{p}_{i,r1}\left(x\right)=p_i\left(x\right)p_{i,r1}\left(y\right)=p_i\left(y\right)\\ p_{i,r2}\left(x\right)=p_i\left(x\right)+l*\cos (k*\mathrm{d\θ})p_{i,r2}\left(y\right)=p_i\left(y\right)+l*\sin (k*\mathrm{d\θ})\\ p_{i,r3}\left(x\right)=p_i\left(x\right)+\sqrt{l^2+w^2}*\cos (k*\mathrm{d\θ}+\arctan (w/l))p_{i,r3}\left(y\right)=p_i\left(y\right)+\sqrt{l^2+w^2}*\sin (k*\mathrm{d\θ}+\arctan (w/l))\\ p_{i,r4}\left(x\right)=p_i\left(x\right)+w*\cos (k*\mathrm{d\θ}+\mathrm{\π})p_{i,r4}\left(y\right)=p_i\left(y\right)+w*\sin (k*\mathrm{d\θ}+\mathrm{\π})\end{array}\right.\end{array}$

K=0,1,2 ... 18 (formula 7)

Wherein p _{i, r1}(x) be p _{i, r1}x coordinate figure, p _{i, r1}(y) be p _{i, r1}y coordinate figure, remaining in like manner, l is rectangle length, w is rectangle width, the anglec of rotation interval that d θ is rectangle, k be rotation number of times or be called iterations;

(2) according to the neighborhood determination methods of the 4th step, search p _ithe neighbor point P of point _{i_near}, then judge P _{i_near}whether at p _{i, r1}, p _{i, r2}, p _{i, r3}, p _{i, r4}in the rectangular extent forming, will in scope, save as point set P _{i_r, k}, and calculate its elevation variance

Elevation average: wherein z _jfor point set P _{i_r, k}in j point height value, m is point set P _{i_r, k}point number;

Elevation variance yields: ${\mathrm{\δ}}_{i,k}^2=\frac{\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{(z_j-{\stackrel{\‾}{z}}_{i,k})}^2}{m}$

(3) according to above two steps, calculation level p _iall rectangles (k=0,1,2 ... 18) the elevation variance yields of point in neighborhood, obtains the set of elevation variance yields relatively gather respectively the elevation variance yields in Ψ, record elevation variance yields minimum, second little and the 3rd three little values with with corresponding rectangle numbering k ₁, k ₂, k ₃if meet following two current some p of equation simultaneously _ibe designated as road center point;

$\left\{\begin{array}{}\end{array}\begin{array}{c}\frac{({\mathrm{\&delta;}}_{i,\min 1}^2-{\mathrm{\&delta;}}_{i,\min 1}^2)}{({\mathrm{\&delta;}}_{i,\min 2}^2-{\mathrm{\&delta;}}_{i,\min 1}^2)}<\mathrm{\&psi;}\\ 9-\mathrm{\&mu;}\&le;|k_1-k_2|\&le;9+\mathrm{\&mu;}\end{array}\right.$

In formula, ψ sets the coarse threshold value of road for user, and μ sets road turning threshold value for user;

(4) repeat above step, until complete to judgement a little, and remember that road-center point set is P _k.

Described the 4th step comprises:

According to the road candidate point set P calculating _cwith road-center point set P _k, road candidate point is concentrated the point that may contain the smooth regions such as a large amount of parking lots, sports ground, roof, by following method, these points is rejected:

(1) by point set P _cand P _kcarry out the division of specification grid, and 1 and 2 represent that kind is a little road candidate point or road center point with coding;

(2) to point set P _cin each point, search for its point that around attribute is 2 in r/2, retain this point if having, otherwise be just judged as non-road waypoint, parameter r/2 represents the ratio of 1/2nd and mesh spacing of road width.

Advantage of the present invention is:

The present invention mainly solves is not having reflection strength or reflection strength can not better react atural object attribute, or there is no the situations such as aviation image, the different road extraction of carrying out according to road topology design feature from surrounding environment.Only adopt cloud data to carry out the road extraction in territory, forest land, mound, the less dependence to raw data.

Brief description of the drawings

Below in conjunction with drawings and Examples, the invention will be further described:

Fig. 1 is the division schematic diagram that the present invention arrives airborne LiDAR data allocations each grid;

Fig. 2 is w × w of the present invention grid schematic diagram;

Fig. 3 is that regular grid of the present invention is divided schematic diagram;

Fig. 4 is neighborhood point set schematic diagram of the present invention;

Fig. 5 is rectangular extent point schematic diagram of the present invention;

Fig. 6 is rotation rectangle schematic diagram of the present invention.

Embodiment

Embodiment: shown in figure 1-6, the invention provides a kind of road data extracting method of analyzing based on rectangular neighborhood, especially a kind of territory, forest land, airborne LiDAR mound method for extracting roads of analyzing based on rectangular neighborhood, overall framework is as follows:

The first step: data pre-service, first calculation level cloud spacing, sets up regular grid tissue, then adopts the quick mathematical morphology filter method of rule-based grid to carry out filtering, rejects non-ground point;

Second step: screening road candidate point, tentatively determine road point set, dwindle determination range, to the point within the scope of its neighborhood of each point search, and compare its height value, if all in threshold range, as road candidate point.

The 3rd step: whether extract road axis, to each road candidate point, calculate the interior elevation variance yields of rectangular extent of different orientations, be road center point according to the current point of the diversity judgement between variance yields.

The 4th step: non-road waypoint is rejected.

Concrete detailed step is as follows:

1,, according to airborne LiDAR data area and quantity survey (surveying) point cloud spacing pt_dis, wherein airborne LiDAR data area is the scope of giving directions cloud to cover, generally represents with the coordinate at four angles;

(1) first read airborne LiDAR data area (xmin, xmax, ymin, ymax) by airborne LiDAR, the quantity pt_num of statistics point;

(2) calculation level cloud spacing pt_dis:

$\mathrm{pt}\_\mathrm{dis}=\frac{(x\max -x\min )*(y\max -y\min )}{\sqrt{\mathrm{pt}\_\mathrm{num}}}$ (formula 1)

2, airborne LiDAR data are carried out to regular grid division, are convenient to search for neighborhood point, concrete steps and being described as follows:

(1) according to a cloud spacing pt_dis with calculate needs, set up a regular grid that covers some cloud scope (xmin, xmax, ymin, ymax), the present invention gets mesh spacing d and is preferably 2～3 times of a cloud spacing pt_dis.

(2) according to point coordinate by airborne LiDAR data allocations in each grid.

$\left\{\begin{array}{c}\mathrm{XAxisNum}=\mathrm{int}((x_{\max }-x_{\min })/d)+1\\ \mathrm{YAxisNum}=\mathrm{int}((y_{\max }-y_{\min })/d)+1\end{array}\right.$ (formula 2)

GridNum=XAxisNum*YAxisNum (formula 3)

XAxisNum, YAxisNum are the grid number of x, y direction, and GridNum is grid sum.

Any point p in airborne LiDAR data _i(x _i, y _i) grid at place is:

$\left\{\begin{array}{c}{X}_i=\mathrm{int}((x_i-x_{\min })/d)\\ Y_i=\mathrm{int}((y_i-y_{\min })/d)\\ G_i=\mathrm{int}((x_i-x_{\min })/d)+X_i\&times;\mathrm{int}((y_i-y_{\min })/d)\end{array}\right.$ (formula 4)

The grid number of the one-dimension array record that wherein Xi is the directions X grid number at Pi place, Y-direction grid number that Yi is Pi place, Gi is Pi place, as shown in Figure 1, data in figure are calculated as XAxisNum=11, YAxisNum=7, GridNum=77, Xi=7, Yi=4, Gi=4*11+7=51, represents that in the 51st grid at Pi o'clock.

3, by fundamental operation and the extended arithmetic thereof of mathematical morphology, airborne LiDAR data are carried out to the filtering of non-ground, retain and comprise waypoint in interior ground point data.

Mathematical morphology filter taking regular grid as unit, grid Corrosion results is the minimum height value in w × w grid, expansion results is maximum elevation value in w × w grid, wherein w represents the size of window, namely sizing grid, as shown in Figure 2, its schematic diagram that is w=2.

$\left\{\begin{array}{c}\mathrm{Grid}{\left(\mathrm{A\&Theta;B}\right)}_i=\mathrm{Grid}\left(z_i^{\&prime;}\right)=\underset{{\mathrm{Grid}}_i\&Subset;w}{\min }\left({\mathrm{Grid}}_i\left(z_{\min }\right)\right)\\ \mathrm{Grid}{(A\&CirclePlus;B)}_i=\mathrm{Grid}\left(z_i^{\&Prime;}\right)=\underset{{\mathrm{Grid}}_i\&Subset;w}{\max }\left({\mathrm{Grid}}_i\left(z_i^{\&prime;}\right)\right)\end{array}\right.$ (formula 5)

Wherein formula 5 is typical regular mathematical morphology formula, and A represents goal set, B representative structure element, and A Θ B represents corrosion, represent to expand, z in formula _i' and z _iheight value after the corrosion of ' ' be respectively grid cell and expansion, Grid _i(z _min) be elevation minimum value in grid, the size that w is window.

Concrete steps are as follows:

(1) by mesh spacing computing grid quantity, and be grid storage allocation space, point is stored in corresponding grid cell by planimetric coordinates one by one, simultaneously with this grid cell in the minimum height value of having preserved compare (the minimum height value of initial default is a very large value), the minimum height value of record after relatively;

(2) travel through each grid, minimum height value in w × w grid relatively around, adopt first formula of formula 5, getting minimum value is the height value after current grid corrosion, and this grid maximum elevation value is made as to the height value after corrosion, specifically the minimum height value in grid in certain limit is around compared, retain minimum height value, and its assignment is given to this variable of maximum elevation value of current grid;

(3) travel through each grid, w × w grid internal corrosion height value relatively around, adopt second formula of formula 5, getting maximal value is the height value after current grid expands, and minimum this grid height value is made as to the height value after expansion, specifically the corrosion height value in grid in certain limit around (result that previous step is calculated) is compared, retain minimum height value, and its assignment is given to this variable of minimum height value of current grid;

(4) travel through each grid, the height value (i.e. minimum height value in this grid) after asking for each point height value in grid and expanding poor, in the time that its absolute value is less than setting threshold, this point is ground point, otherwise is non-ground point.

Through above four step processing, can complete fast the filtering processing to airborne LiDAR data.

4, estimate the width road_width of road, choose road candidate point in neighborhood.

(1) for each some p _i(x _i, y _i) calculate its place grid by formula 4;

(2) calculate search radius, $r=\frac{\mathrm{road}\_\mathrm{width}}{d};$

(3) p that sets up an office _i(x _i, y _i) rectangular neighborhood in point set be

P _{i_near}={ p _j(x _j, y _j) | x _j∈ (x _i-r, x _i+ r), y _j∈ (y _i-r, y _i+ r) (wherein j is period or perhaps the subscript of point, for the subscript that is different from current some pi i), thus centered by current some pi, search is the point set P in grid around _{i_near}, as shown in Figure 4, be p with the point of arrow points _i(x _i, y _i), when R=3, P _{i_near}point set is the point marking.

(4) calculate P _{i_near}in each point and p _i(x _i, y _i) difference of elevation absolute value sum Σ H _iif, Σ H _ibe less than threshold value, current point is road candidate point.

$\mathrm{\&Sigma;}{H}_i=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|z_j-z_i|$ (formula 6)

Zj is P _{i_near}in j point height value, zi is p _i(x _i, y _i) height value, n is P _{i_near}always count.

(5) repeat above step, until complete to judgement a little, and the road candidate point set of reservation is designated as to P _c.

5, search road axis

To each some cloud, set up rectangular search neighborhood, the length l of rectangle is preferably 2～3 times of width w.

(1) to each road candidate point p _i, calculate rectangular extent point p according to formula 7 _{i, r1}, p _{i, r2}, p _{i, r3}, p _{i, r4}, i.e. four of rectangle summits, p _{i, r1}overlap with current point, all the other press counterclockwise arrangement, and 360 degree are divided into n interval by the anglec of rotation that wherein d θ is rectangle, and the present invention gets d θ=20 degree, i.e. n=18, as shown in Figure 4.

$\begin{array}{c}\left\{\begin{array}{c}{p}_{i,r1}\left(x\right)=p_i\left(x\right)p_{i,r1}\left(y\right)=p_i\left(y\right)\\ p_{i,r2}\left(x\right)=p_i\left(x\right)+l*\cos (k*\mathrm{d\&theta;})p_{i,r2}\left(y\right)=p_i\left(y\right)+l*\sin (k*\mathrm{d\&theta;})\\ p_{i,r3}\left(x\right)=p_i\left(x\right)+\sqrt{l^2+w^2}*\cos (k*\mathrm{d\&theta;}+\arctan (w/l))p_{i,r3}\left(y\right)=p_i\left(y\right)+\sqrt{l^2+w^2}*\sin (k*\mathrm{d\&theta;}+\arctan (w/l))\\ p_{i,r4}\left(x\right)=p_i\left(x\right)+w*\cos (k*\mathrm{d\&theta;}+\mathrm{\&pi;})p_{i,r4}\left(y\right)=p_i\left(y\right)+w*\sin (k*\mathrm{d\&theta;}+\mathrm{\&pi;})\end{array}\right.\end{array}$

K=0,1,2 ... 18 (formula 7)

Wherein p _{i, r1}(x) be p _{i, r1}x coordinate figure, p _{i, r1}(y) be p _{i, r1}y coordinate figure, remaining in like manner, l is rectangle length, w is rectangle width, the anglec of rotation interval that d θ is rectangle, k be rotation number of times or be called iterations.

(2) according to the neighborhood determination methods of the 4th step, search p _ithe neighbor point P of point _{i_near}, then judge P _{i_near}whether at p _{i, r1}, p _{i, r2}, p _{i, r3}, p _{i, r4}in the rectangular extent forming, will in scope, save as point set P _{i_r, k}, and calculate its elevation variance

Elevation average: wherein z _jfor point set P _{i_r, k}in j point height value, m is point set P _{i_r, k}point number;

Elevation variance yields: ${\mathrm{\&delta;}}_{i,k}^2=\frac{\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}{(z_j-{\stackrel{\&OverBar;}{z}}_{i,k})}^2}{m}$

(3) according to above two steps, calculation level p _iall rectangles (k=0,1,2 ... 18) the elevation variance yields of point in neighborhood, obtains the set of elevation variance yields relatively gather respectively the elevation variance yields in Ψ, record elevation variance yields minimum, second little and the 3rd three little values with with corresponding rectangle numbering k ₁, k ₂, k ₃if meet following two current some p of equation simultaneously _ibe designated as road center point.

$\left\{\begin{array}{}\end{array}\begin{array}{c}\frac{({\mathrm{\&delta;}}_{i,\min 1}^2-{\mathrm{\&delta;}}_{i,\min 1}^2)}{({\mathrm{\&delta;}}_{i,\min 2}^2-{\mathrm{\&delta;}}_{i,\min 1}^2)}<\mathrm{\&psi;}\\ 9-\mathrm{\&mu;}\&le;|k_1-k_2|\&le;9+\mathrm{\&mu;}\end{array}\right.$

In formula, ψ sets the coarse threshold value of road for user, and the present invention is taken as 0.5, μ and sets road turning threshold value for user, and the present invention is taken as 1, and angle of turn is 20 degree to the maximum.

(4) repeat above step, until complete to judgement a little, and remember that road-center point set is P _k.

6, non-road waypoint is rejected

The road candidate point set P calculating respectively according to the 4th, 5 steps _cwith road-center point set P _k, road candidate point is concentrated the point that may contain the smooth regions such as a large amount of parking lots, sports ground, roof, by following method, these points is rejected.

(1) adopt and the similar method of the 2nd step, by point set P _cand P _kcarry out the division of specification grid, and 1 and 2 represent that kind is a little road candidate point or road center point with coding.

(2) to point set point set P _cin each point, search for its point that around attribute is 2 within the scope of r/2, retain this point if having, otherwise be just judged as non-road waypoint.Parameter r/2 represents the ratio of 1/2nd and mesh spacing of road width, and in the present invention, r/2 is taken as 2.

Certainly above-described embodiment is only explanation technical conceive of the present invention and feature, and its object is to allow person skilled in the art can understand content of the present invention and implement according to this, can not limit the scope of the invention with this.All equivalent transformation or modification that according to the present invention, the Spirit Essence of main technical schemes does, within all should being encompassed in protection scope of the present invention.

Claims (6)

1. a road data extracting method of analyzing based on rectangular neighborhood, is characterized in that it comprises:

The first step: data pre-service, calculation level cloud spacing, sets up regular grid tissue, then adopts the quick mathematical morphology filter method of rule-based grid to carry out filtering, rejects non-ground point;

Second step: screening road candidate point, tentatively determine road point set, dwindle determination range, to the point within the scope of its neighborhood of each point search, and compare its height value, if all in threshold range, as road candidate point;

The 3rd step: whether extract road axis, reject smooth region point set, to each road candidate point, calculate the interior elevation variance yields of rectangular extent of different orientations, be road center point according to the current point of the diversity judgement between variance yields;

The 4th step: road connects, connects the road being interrupted, level and smooth.

2. a kind of road data extracting method of analyzing based on rectangular neighborhood according to claim 1, is characterized in that: the described first step comprises:

(1) first read cloud data scope (xmin, xmax, ymin, ymax), the quantity pt_Num of statistics point;

(2) calculation level cloud spacing:

$\mathrm{pt}\_\mathrm{dis}=\frac{(x\max -x\min )*(y\max -y\min )}{\sqrt{\mathrm{pt}\_\mathrm{num}}}$

(3) according to a cloud spacing pt_dis and calculating needs, set up a regular grid that covers some cloud scope (xmin, xmax, ymin, ymax);

(4) according to point coordinate by airborne LiDAR data allocations in each grid,

$\left\{\begin{array}{c}\mathrm{XAxisNum}=\mathrm{int}((x_{\max }-x_{\min })/d)+1\\ \mathrm{YAxisNum}=\mathrm{int}((y_{\max }-y_{\min })/d)+1\end{array}\right.$ (formula 2)

GridNum=XAxisNum*YAxisNum (formula 3)

XAxisNum, YAxisNum are the grid number of x, y direction, and GridNum is grid sum;

Any point p in airborne LiDAR data _i(x _i, y _i) grid at place is:

$\left\{\begin{array}{c}{X}_i=\mathrm{int}((x_i-x_{\min })/d)\\ Y_i=\mathrm{int}((y_i-y_{\min })/d)\\ G_i=\mathrm{int}((x_i-x_{\min })/d)+X_i\&times;\mathrm{int}((y_i-y_{\min })/d)\end{array}\right.$ (formula 4)

Wherein X _ifor P _ithe directions X grid number at place, Y _ifor P _ithe Y-direction grid number at place, G _ifor P _ithe grid number of the one-dimension array record at place;

Airborne LiDAR data are carried out to the filtering of non-ground, retain and comprise waypoint in interior ground point data, mathematical morphology filter taking regular grid as unit, grid Corrosion results is the minimum height value in w × w grid, expansion results is maximum elevation value in w × w grid, wherein w represents the size of window, and carries out the computing of formula 5:

$\left\{\begin{array}{c}\mathrm{Grid}{\left(\mathrm{A\&Theta;B}\right)}_i=\mathrm{Grid}\left(z_i^{\&prime;}\right)=\underset{{\mathrm{Grid}}_i\&Subset;w}{\min }\left({\mathrm{Grid}}_i\left(z_{\min }\right)\right)\\ \mathrm{Grid}{(A\&CirclePlus;B)}_i=\mathrm{Grid}\left(z_i^{\&Prime;}\right)=\underset{{\mathrm{Grid}}_i\&Subset;w}{\max }\left({\mathrm{Grid}}_i\left(z_i^{\&prime;}\right)\right)\end{array}\right.$ (formula 5)

Wherein A represents goal set, B representative structure element, and A Θ B represents corrosion, represent to expand, z in formula _i' and z _iheight value after the corrosion of ' ' be respectively grid cell and expansion, Grid _i(z _min) be elevation minimum value in grid.

3. a kind of road data extracting method of analyzing based on rectangular neighborhood according to claim 2, is characterized in that: the described first step further comprises:

(1) by mesh spacing computing grid quantity, and be grid storage allocation space, point is stored in corresponding grid cell by planimetric coordinates one by one, simultaneously with this grid cell in the minimum height value of having preserved compare, the minimum height value of record after relatively;

(2) travel through each grid, minimum height value in w × w grid relatively around, first formula of employing formula 5, getting minimum value is the height value after current grid corrosion, and this grid maximum elevation value is made as to the height value after corrosion;

(3) travel through each grid, w × w grid internal corrosion height value relatively around, second formula of employing formula 5, getting maximal value is the height value after current grid expands, and minimum this grid height value is made as to the height value after expansion;

(4) travel through each grid, the height value after asking for each point height value in grid and expanding poor, in the time that its absolute value is less than setting threshold, this point is ground point, otherwise is non-ground point.

4. a kind of road data extracting method of analyzing based on rectangular neighborhood according to claim 3, is characterized in that: described second step comprises:

(1) for each some p _i(x _i, y _i) calculate its place grid by formula 4;

(2) calculate search radius, $r=\frac{\mathrm{road}\_\mathrm{width}}{d};$

(3) p that sets up an office _i(x _i, y _i) rectangular neighborhood in point set be

P _{i_near}={ p _j(x _j, y _j) | x _j∈ (x _i-r, x _i+ r), y _j∈ (y _i-r, y _i+ r) } (wherein j is period or perhaps the subscript of point, and distinguishes with current some pi), thus centered by current some pi, search is the point set P in grid around _{i_near};

(4) calculate P _{i_near}in each point and p _i(x _i, y _i) difference of elevation absolute value sum Σ H _iif, Σ H _ibe less than threshold value, current point is road candidate point;

$\mathrm{\&Sigma;}{H}_i=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|z_j-z_i|$ (formula 6)

Z _jp _{i_near}in j point height value, z _ip _i(x _i, y _i) height value, n is P _{i_near}always count;

(5) repeat above step, until complete to judgement a little, and be P by the road candidate point set of reservation _c.

5. a kind of road data extracting method of analyzing based on rectangular neighborhood according to claim 4, is characterized in that: described the 3rd step comprises:

(1) to each road candidate point p _i, calculate rectangular extent point p according to formula 7 _{i, r1}, p _{i, r2}, p _{i, r3}, p _{i, r4}, i.e. four of rectangle summits, p _{i, r1}overlap with current point, all the other press counterclockwise arrangement, and 360 degree are divided into n interval by the anglec of rotation that wherein d θ is rectangle;

$\begin{array}{c}\left\{\begin{array}{c}{p}_{i,r1}\left(x\right)=p_i\left(x\right)p_{i,r1}\left(y\right)=p_i\left(y\right)\\ p_{i,r2}\left(x\right)=p_i\left(x\right)+l*\cos (k*\mathrm{d\&theta;})p_{i,r2}\left(y\right)=p_i\left(y\right)+l*\sin (k*\mathrm{d\&theta;})\\ p_{i,r3}\left(x\right)=p_i\left(x\right)+\sqrt{l^2+w^2}*\cos (k*\mathrm{d\&theta;}+\arctan (w/l))p_{i,r3}\left(y\right)=p_i\left(y\right)+\sqrt{l^2+w^2}*\sin (k*\mathrm{d\&theta;}+\arctan (w/l))\\ p_{i,r4}\left(x\right)=p_i\left(x\right)+w*\cos (k*\mathrm{d\&theta;}+\mathrm{\&pi;})p_{i,r4}\left(y\right)=p_i\left(y\right)+w*\sin (k*\mathrm{d\&theta;}+\mathrm{\&pi;})\end{array}\right.\end{array}$

K=0,1,2 ... 18 (formula 7)

Wherein p _{i, r1}(x) be p _{i, r1}x coordinate figure, p _{i, r1}(y) be p _{i, r1}y coordinate figure, remaining in like manner, l is rectangle length, w is rectangle width, the anglec of rotation interval that d θ is rectangle, k be rotation number of times or be called iterations;

(2) according to the neighborhood determination methods of the 4th step, search p _ithe neighbor point P of point _{i_near}, then judge P _{i_near}whether at p _{i, r1}, p _{i, r2}, p _{i, r3}, p _{i, r4}in the rectangular extent forming, will in scope, save as point set P _{i_r, k}, and calculate its elevation variance

Elevation average: wherein z _jfor point set P _{i_r, k}in j point height value, m is point set P _{i_r, k}point number;

Elevation variance yields: ${\mathrm{\&delta;}}_{i,k}^2=\frac{\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}{(z_j-{\stackrel{\&OverBar;}{z}}_{i,k})}^2}{m}$

(3) according to above two steps, calculation level p _iall rectangles (k=0,1,2 ... 18) the elevation variance yields of point in neighborhood, obtains the set of elevation variance yields relatively gather respectively the elevation variance yields in Ψ, record elevation variance yields minimum, second little and the 3rd three little values with with corresponding rectangle numbering k ₁, k ₂, k ₃if meet following two current some p of equation simultaneously _ibe designated as road center point;

$\left\{\begin{array}{}\end{array}\begin{array}{c}\frac{({\mathrm{\&delta;}}_{i,\min 1}^2-{\mathrm{\&delta;}}_{i,\min 1}^2)}{({\mathrm{\&delta;}}_{i,\min 2}^2-{\mathrm{\&delta;}}_{i,\min 1}^2)}<\mathrm{\&psi;}\\ 9-\mathrm{\&mu;}\&le;|k_1-k_2|\&le;9+\mathrm{\&mu;}\end{array}\right.$

In formula, ψ sets the coarse threshold value of road for user, and μ sets road turning threshold value for user;

(4) repeat above step, until complete to judgement a little, and remember that road-center point set is P _k.

6. a kind of road data extracting method of analyzing based on rectangular neighborhood according to claim 5, is characterized in that: described the 4th step comprises:

According to the road candidate point set P calculating _cwith road-center point set P _k, road candidate point is concentrated the point that may contain the smooth regions such as a large amount of parking lots, sports ground, roof, by following method, these points is rejected:

(1) by point set P _cand P _kcarry out the division of specification grid, and 1 and 2 represent that kind is a little road candidate point or road center point with coding;

(2) to point set P _cin each point, search for its point that around attribute is 2 in r/2, retain this point if having, otherwise be just judged as non-road waypoint, parameter r/2 represents the ratio of 1/2nd and mesh spacing of road width.