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

Abstract

The invention discloses a method for extracting road data based on rectangular neighborhood analysis, which includes a data preprocessing step, a step of screening road candidate points, a step of extracting a road center line and eliminating flat area point sets, and a road connection step. The invention mainly solves the problem of extracting road data according to the difference between road topological structure characteristics and surrounding environment when there is no reflection intensity or the reflection intensity cannot better reflect the ground object attributes, or there is no aerial image.

Description

A road data extraction method based on rectangular neighborhood analysis

technical field

The invention relates to a road data extraction method based on rectangular neighborhood analysis.

Background technique

The extraction of buildings and roads has always been the focus and difficulty of airborne LiDAR (Light Detection And Ranging) data feature extraction. In terms of road extraction, there are more researches on the extraction of road networks in urban areas, and less research on road extraction methods in wild areas. Among them, airborne LiDAR is a system that integrates three technologies of laser ranging, global positioning system (GPS) and inertial navigation system (INS), and is used to acquire data and generate accurate three-dimensional terrain (DEM); airborne LiDAR ( Or point cloud data) is the point data on the surface of an object collected by a 3D laser scanner. Each point has three-dimensional coordinates in space, and some also have information such as color and anti-color intensity; Non-ground points such as noise, trees, and houses; point cloud organization refers to managing point clouds according to a certain data format, which is convenient for data query and display; road extraction refers to automatically extracting road points from point cloud data. Existing road extraction methods mainly include road extraction methods based on strong reflection information, road extraction methods based on airborne LiDAR data and remote sensing images, etc. Commonly used road data extraction methods based on point clouds, reflection intensity, and images are as follows:

(1) Automatic extraction method of structured roads from vehicle-mounted laser scanning data

Vehicle-mounted laser scanning refers to setting up a laser scanner on a mobile platform such as a car for dynamic vehicles. The difference from airborne laser scanning is that the range of airborne laser scanning is larger and the speed is faster, while the data obtained by vehicle-mounted laser scanning are all Along the road, so extraction is easier.

Structured roads refer to urban roads or expressways, etc. Both sides of these roads have features higher than road points. Road extraction is completed through steps such as scanning line extraction, ground point cloud filtering, automatic road boundary extraction, and road boundary tracking. For the specific method, please refer to "Automatic Extraction Method of Structured Roads from Vehicular Laser Scanning Data" in "Journal of Surveying and Mapping", Volume 42, Issue 2, 2013.

(2) Road extraction based on LiDAR echo information

First, a method (such as layer-by-layer encryption TIN method) is used to extract or generate a digital terrain model (DTM), and road information is extracted from the DTM according to the echo information of the point cloud, and finally the noise is deleted through the filtering algorithm of searching isolated points point. The specific method is detailed in "Surveying and Mapping Science", Volume 36, Issue 2, "Road Extraction Based on LiDAR Echo Information" in 2011.

(3) Road extraction based on airborne LiDAR roughness index and echo intensity

First, the normalized digital surface model and roughness index are derived from the point cloud data, and then multi-resolution segmentation is performed on the registered multi-source data, and then the roughness index, echo intensity, road description factor and other features are used for classification , and finally, remove road noise and obtain an accurate road network. For details, see "Road Extraction Based on Airborne LiDAR Roughness Index and Echo Intensity" in Journal of Surveying and Mapping Science and Technology, Volume 30, Issue 1, 2013.

(4) Urban road network extraction based on airborne LiDAR and high-resolution remote sensing images

First, the two are accurately registered, and then the vegetation information and building information are removed by the method of pseudo road information removal to obtain the basic road outline, and then the centerline of the road is extracted by using the morphological refinement algorithm, and finally programmed in ArcGIS and Matlab Under the environment, the improved road pruning algorithm (IRT) is realized, and the smooth and coherent urban road network is obtained by using this algorithm for road pruning. For details, see "Urban Road Network Extraction Based on Airborne LiDAR and High Resolution Remote Sensing Images" in "Remote Sensing Technology and Application", Volume 28, Issue 4, 2013.

However, the laser echo intensity is directly proportional to the reflectivity of ground objects, that is, the greater the reflectivity, the higher the laser intensity value, but the echo intensity recorded by the laser radar system is an uncorrected value, which will be affected by atmospheric attenuation, laser incident Influenced by many factors such as angle, and the correction of laser echo intensity is still a relatively difficult subject, it is difficult for laser echo intensity to truly reflect the reflectivity information of ground objects. The extraction effect of the road extraction method based on remote sensing images mainly depends on the clarity of the images, but the remote sensing images collected under cloudy or insufficient visibility often cannot reflect the attributes of ground features well, so the classification is also more difficult.

Contents of the invention

The purpose of the present invention is to provide a method for extracting road data based on rectangular neighborhood analysis. Under the condition of not using reflection intensity and remote sensing images, by analyzing the structural characteristics of roads in hilly forest areas, the rectangular neighborhood judgment analysis method is used for road extraction. extraction.

The technical scheme of the present invention is: a kind of road data extraction method based on rectangular neighborhood analysis, it comprises:

The first step: data preprocessing, calculate point cloud spacing, establish regular grid organization, and then use fast mathematical morphology filtering method based on regular grid to filter, and eliminate non-ground points;

Step 2: Screen candidate road points, preliminarily determine the road point set, narrow the judgment range, search for points within its neighborhood for each point, and compare their elevation values. If they are all within the threshold range, they will be used as road candidate points ;

Step 3: Extract the road centerline, remove the flat area point set, calculate the elevation variance value within the rectangular range of different azimuths for each road candidate point, and judge whether the current point is the road center point according to the difference between the variance values ;

The fourth step: road connection, connect and smooth the discontinuous road.

On the basis of the above-mentioned technical solutions, the following subsidiary technical solutions are further included:

The first step includes:

(1) First read the point cloud data range (xmin, xmax, ymin, ymax), and count the number of points pt_Num;

(2) Calculate point cloud spacing:

$\mathrm{<span\; class="notranslate">\; pt</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; dis</span>}<span\; class="notranslate">\; =</span>\frac{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; max</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; min</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}\mathrm{<span\; class="notranslate">\; max</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}\mathrm{<span\; class="notranslate">\; min</span>}<span\; class="notranslate">\; )</span>}{\sqrt{\mathrm{<span\; class="notranslate">\; pt</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; num</span>}}}$

(3) According to the point cloud spacing pt_dis and calculation needs, establish a regular grid covering the point cloud range (xmin, xmax, ymin, ymax);

(4) distribute the airborne LiDAR data into individual grids according to the coordinates of the points,

$\left\{\begin{array}{c}\mathrm{XAxisNum}=\mathrm{int}((x_{\max }-x_{\min })/d)+1\\ \mathrm{YAxisNum}=\mathrm{int}(({\mathrm{the\; y}}_{\max }-{\mathrm{the\; y}}_{\min })/d)+1\end{array}\right.$ (Formula 2)

GridNum＝XAxisNum*YAxisNum (Formula 3)

XAxisNum and YAxisNum are the number of grids in the x and y directions, and GridNum is the total number of grids;

The grid of any point p _i ( _xi , y _i ) in the airborne LiDAR data is:

$\left\{\begin{array}{c}{x}_i=\mathrm{int}((x_i-x_{\min })/d)\\ Y_i=\mathrm{int}(({\mathrm{the\; y}}_i-{\mathrm{the\; y}}_{\min })/d)\\ G_i=\mathrm{int}((x_i-x_{\min })/d)+x_i\&times;\mathrm{int}(({\mathrm{the\; y}}_i-{\mathrm{the\; y}}_{\min })/d)\end{array}\right.$ (Formula 4)

Among them, p _i is the number of grids in the X direction where p _i is located, Yi is the number of grids in the Y direction where Pi is located, and Gi is the number of grids recorded in the one-dimensional array where Pi is located;

Carry out non-ground filtering on the airborne LiDAR data, retain ground point data including road points, and use regular grid as the unit of mathematical morphology filtering. The result of grid erosion is the minimum elevation value in the w×w grid, and the expansion The result is the maximum elevation value in the w×w grid, where w represents the size of the window, and the operation of formula 5 is performed:

$\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)

where A represents the target set, B represents the structural element, AΘB represents corrosion, In the formula, zi _′ and zi _′ ′ are the elevation values after grid cell erosion and expansion respectively, and Grid _i (z _min ) is the minimum elevation value in the grid.

The first step further includes:

(1) Calculate the number of grids from the grid spacing, and allocate memory space for the grid, store the points one by one in the corresponding grid unit according to the plane coordinates, and compare with the minimum elevation value saved in the grid unit , record the minimum elevation value after comparison;

(2) Traverse each grid, compare the minimum elevation value in the surrounding w×w grids, use the first formula of formula 5, take the minimum value as the elevation value of the current grid after corrosion, and set the maximum value of the grid The elevation value is set to the elevation value after corrosion;

(3) Traversing through each grid, comparing the corrosion elevation values in the surrounding w×w grids, using the second formula of formula 5, taking the maximum value as the elevation value after the expansion of the current grid, and making the grid minimum The elevation value is set to the inflated elevation value;

(4) Traversing each grid, calculating the difference between the elevation value of each point in the grid and the inflated elevation value, when its absolute value is less than the set threshold, the point is a ground point, otherwise it is a non-ground point.

The second step includes:

(1) For each point p _i (x _i , y _i ), calculate its grid by formula 4;

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

(3) Let the point set in the rectangular neighborhood of point p _i (xi _, y _i ) be

P _{i_near} ＝{p _j (x _j ,y _j )|x _j ∈(xi _-r , _xi +r),y _j ∈(y _i -r,y _i +r)} (where j is the point The point number or subscript is different from the current point pi), so the point set P _{i_near} in the surrounding grid is searched with the current point pi as the center;

(4) Calculate the sum ΣH _i of the absolute value of the elevation difference between each point in P _{i_near} and p _i (xi _, y _i ), if ΣH _i is less than the threshold, the current point is a road candidate point;

$\mathrm{\&Sigma;}{h}_i=\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}|z_j-z_i|$ (Formula 6)

zj is the elevation value of the jth point in P _{i_near} , zi is the elevation value of p _i ( _xi , y _i ), and n is the total number of points in P _{i_near} ;

(5) Repeat the above steps until the judgment of all points is completed, and the reserved road candidate points are set as P _c .

The third step includes:

(1) For each road candidate point p _i , calculate the rectangular range points p _i,r1 , p _i,r2 , p _i,r3 , p _i,r4 according to formula 7, that is, the four vertices of the rectangle, p _i,r1 It coincides with the current point, and the other points are arranged counterclockwise, where dθ is the rotation angle of the rectangle, and divides 360 degrees into n intervals;

$\begin{array}{c}\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{{\mathrm{<span\; class="notranslate">\; l</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{{\mathrm{<span\; class="notranslate">\; l</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>\end{array}\right.\end{array}$

k＝0,1,2...18 (Formula 7)

Among them, p _{i, r1} (x) is the x coordinate value of p _{i, r1 , p i, r1} (y) is the y coordinate value of p _{i, r1} , and the rest are the same, l is the length of the rectangle, w is the width of the rectangle, dθ is the rotation angle interval of the rectangle, and k is the number of rotations or the number of iterations;

(2) According to the neighborhood judgment method in step 4, search for the neighboring point P _{i_near} of point p _i , and then judge whether P _{i_near} is formed by p _i,r1 , p _i,r2 , p _i,r3 , p _i,r4 Within the rectangular range, save the point set P _{i_r,k} within the range, and calculate its elevation variance

Elevation mean: Where z _j is the elevation value of the jth point in the point set P _{i_r,k} , and m is the number of points in the point set P _{i_r,k} ;

Elevation variance value: ${\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 the above two steps, calculate the elevation variance value of points in the neighborhood of all rectangles (k=0,1,2...18) of point p _i , and obtain the elevation variance value set Compare the elevation variance values in the set Ψ respectively, and record the three values with the smallest, second smallest, and third smallest elevation variance values and and the corresponding rectangle numbers k ₁ , k ₂ , k ₃ , if the following two equations are satisfied at the same time, the current point p _i is recorded as the road center point;

$\left\{\begin{array}{}\end{array}\begin{array}{c}\frac{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; \&psi;</span>}\\ \mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}<span\; class="notranslate">\; \&le;</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}\end{array}\right.$

where ψ is the road roughness threshold set by the user, and μ is the road turning threshold set by the user;

(4) Repeat the above steps until the judgment of all points is completed, and record the road center point set as P _k .

The fourth step includes:

According to the calculated road candidate point set P _c and road center point set P _k , the road candidate point set may contain a large number of points in flat areas such as parking lots, sports fields, and roofs. These points are eliminated by the following methods:

(1) Divide the point sets P _c and P _k into standard grids, and use codes 1 and 2 to indicate whether the point type is a road candidate point or a road center point;

(2) For each point in the point set _Pc , search for a point with an attribute of 2 within the r/2 range around it, if there is, keep the point, otherwise it will be judged as a non-road point, and the parameter r/2 represents the width of the road The ratio of one-half of to the grid spacing.

Advantage of the present invention is:

The present invention mainly solves the problem of extracting roads according to the differences between road topological structure characteristics and the surrounding environment when there is no reflection intensity or the reflection intensity cannot better reflect the ground object attributes, or there is no aerial image. Only point cloud data is used for road extraction in hilly areas, which reduces the dependence on original data.

Description of drawings

The present invention will be further described below in conjunction with accompanying drawing and embodiment:

Fig. 1 is the division schematic diagram that airborne LiDAR data is distributed to each grid by the present invention;

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

Fig. 3 is a schematic diagram of regular grid division of the present invention;

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

Fig. 5 is a schematic diagram of rectangular range points of the present invention;

Fig. 6 is a schematic diagram of a rotating rectangle in the present invention.

Detailed ways

Embodiment: Referring to Figures 1-6, the present invention provides a method for extracting road data based on rectangular neighborhood analysis, especially an airborne LiDAR road extraction method based on rectangular neighborhood analysis, the overall framework as follows:

The first step: data preprocessing, first calculate the point cloud spacing, establish a regular grid organization, and then use a fast mathematical morphology filtering method based on a regular grid to filter out non-ground points;

Step 2: Screen candidate road points, preliminarily determine the road point set, narrow the judgment range, search for points within its neighborhood for each point, and compare their elevation values. If they are all within the threshold range, they will be used as road candidate points .

Step 3: extract the road centerline, calculate the elevation variance value within the rectangular range of different azimuths for each road candidate point, and judge whether the current point is the road center point according to the difference between the variance values.

Step 4: Eliminate non-road points.

The detailed steps are as follows:

1. Estimate the point cloud distance pt_dis according to the range and quantity of the airborne LiDAR data, where the airborne LiDAR data range refers to the range covered by the point cloud, which is generally represented by the coordinates of the four corners;

(1) First read the airborne LiDAR data range (xmin, xmax, ymin, ymax) through the airborne LiDAR, and the number of statistical points pt_num;

(2) Calculate the point cloud distance pt_dis:

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

2. Carry out regular grid division on the airborne LiDAR data to facilitate the search for neighboring points. The specific steps and instructions are as follows:

(1) According to the point cloud spacing pt_dis and calculation requirements, a regular grid covering the point cloud range (xmin, xmax, ymin, ymax) is established. The grid spacing d is preferably 2 to 3 times the point cloud spacing pt_dis in the present invention.

(2) Distribute the airborne LiDAR data into individual grids according to the coordinates of the points.

$\left\{\begin{array}{c}\mathrm{XAxisNum}=\mathrm{int}((x_{\max }-x_{\min })/d)+1\\ \mathrm{YAxisNum}=\mathrm{int}(({\mathrm{the\; y}}_{\max }-{\mathrm{the\; y}}_{\min })/d)+1\end{array}\right.$ (Formula 2)

GridNum=XAxisNum*YAxisNum (Formula 3)

XAxisNum and YAxisNum are the number of grids in the x and y directions, and GridNum is the total number of grids.

The grid of any point p _i ( _xi , y _i ) in the airborne LiDAR data is:

$\left\{\begin{array}{c}{x}_i=\mathrm{int}((x_i-x_{\min })/d)\\ Y_i=\mathrm{int}(({\mathrm{the\; y}}_i-{\mathrm{the\; y}}_{\min })/d)\\ G_i=\mathrm{int}((x_i-x_{\min })/d)+x_i\&times;\mathrm{int}(({\mathrm{the\; y}}_i-{\mathrm{the\; y}}_{\min })/d)\end{array}\right.$ (Formula 4)

Among them, Xi is the number of grids in the X direction where Pi is located, Yi is the number of grids in the Y direction where Pi is located, and Gi is the number of grids recorded in the one-dimensional array where Pi is located. As shown in Figure 1, the data in the figure is calculated as XAxisNum =11, YAxisNum=7, GridNum=77, Xi=7, Yi=4, Gi=4*11+7=51, indicating that the Pi point is in the 51st grid.

3. Perform non-ground filtering on the airborne LiDAR data through the basic operation of mathematical morphology and its extended operation, and retain the ground point data including road points.

Mathematical morphological filtering with a regular grid as the unit, the grid erosion result is the minimum elevation value in the w×w grid, and the dilation result is the maximum elevation value in the w×w grid, where w represents the size of the window, that is, Grid size, as shown in Figure 2, which is a schematic diagram of 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)

Among them, Formula 5 is a typical regular mathematical morphology formula, A represents the target set, B represents the structural element, AΘB represents corrosion, In the formula, zi _′ and zi _′ ′ are the elevation values of the grid unit after corrosion and expansion respectively, Grid _i (z _min ) is the minimum elevation value in the grid, and w is the size of the window.

Specific steps are as follows:

(1) Calculate the number of grids from the grid spacing, and allocate memory space for the grid, store the points one by one in the corresponding grid unit according to the plane coordinates, and compare with the minimum elevation value saved in the grid unit (The initial default minimum elevation value is a very large value), record the minimum elevation value after comparison;

(2) Traverse each grid, compare the minimum elevation value in the surrounding w×w grids, use the first formula of formula 5, take the minimum value as the elevation value of the current grid after corrosion, and set the maximum value of the grid The elevation value is set to the elevation value after corrosion, specifically comparing the minimum elevation value in the surrounding grid within a certain range, retaining the minimum elevation value, and assigning it to the variable of the maximum elevation value of the current grid;

(3) Traversing through each grid, comparing the corrosion elevation values in the surrounding w×w grids, using the second formula of formula 5, taking the maximum value as the elevation value after the expansion of the current grid, and making the grid minimum The elevation value is set to the elevation value after expansion. Specifically, compare the corrosion elevation values (results calculated in the previous step) in the grid within a certain range around, keep the minimum elevation value, and assign it to the minimum value of the current grid. The variable of elevation value;

(4) Traversing each grid, calculating the difference between the elevation value of each point in the grid and the height value after expansion (that is, the minimum elevation value in the grid), when its absolute value is less than the set threshold, the point is a ground point, otherwise it is a non-ground point.

After the above four steps of processing, the filtering processing of the airborne LiDAR data can be quickly completed.

4. Estimate the width road_width of the road, and select the road candidate points in the neighborhood.

(1) For each point p _i (x _i , y _i ), calculate its grid by formula 4;

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

(3) Let the point set in the rectangular neighborhood of point p _i (xi _, y _i ) be

P _{i_near} ＝{p _j (x _j ,y _j )|x _j ∈(xi _-r , _xi +r),y _j ∈(y _i -r,y _i +r)} (where j is the point The point number or subscript, in order to distinguish it from the subscript i) of the current point pi, so that the current point pi is used as the center to search for the point set P _{i_near} in the surrounding grid, as shown in Figure 4, with arrows pointing to The point of is p _i ( _xi , y _i ), when R=3, the point set of P _{i_near} is the marked point.

(4) Calculate the sum ΣH _i of the absolute value of the elevation difference between each point in P _{i_near} and p _i (xi _, y _i ), if ΣH _i is smaller than the threshold, the current point is a road candidate point.

$\mathrm{\&Sigma;}{h}_i=\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}|z_j-z_i|$ (Formula 6)

zj is the elevation value of the jth point in P _{i_near} , zi is the elevation value of p _i ( _xi , y _i ), and n is the total number of points in P _{i_near} .

(5) Repeat the above steps until the judgment of all points is completed, and record the reserved road candidate point set as P _c .

5. Search for road centerline

For each point cloud, a rectangular search neighborhood is established, and the length l of the rectangle is preferably 2 to 3 times the width w.

(1) For each road candidate point p _i , calculate the rectangular range points p _i,r1 , p _i,r2 , p _i,r3 , p _i,r4 according to formula 7, that is, the four vertices of the rectangle, p _i,r1 Coincident with the current point, all the other points are arranged counterclockwise, wherein dθ is the rotation angle of the rectangle, and 360 degrees are divided into n intervals. The present invention takes dθ=20 degrees, i.e. n=18, as shown in Figure 4.

$\begin{array}{c}\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{{\mathrm{<span\; class="notranslate">\; l</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{{\mathrm{<span\; class="notranslate">\; l</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>\end{array}\right.\end{array}$

k＝0,1,2...18 (Formula 7)

Among them, p _{i, r1} (x) is the x coordinate value of p _{i, r1 , p i, r1} (y) is the y coordinate value of p _{i, r1} , and the rest are the same, l is the length of the rectangle, w is the width of the rectangle, dθ is the rotation angle interval of the rectangle, and k is the number of rotations or the number of iterations.

(2) According to the neighborhood judgment method in step 4, search for the neighboring point P _{i_near} of point p _i , and then judge whether P _{i_near} is formed by p _i,r1 , p _i,r2 , p _i,r3 , p _i,r4 Within the rectangular range, save the point set P _{i_r,k} within the range, and calculate its elevation variance

Elevation mean: Where z _j is the elevation value of the jth point in the point set P _{i_r,k} , and m is the number of points in the point set P _{i_r,k} ;

Elevation variance value: ${\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 the above two steps, calculate the elevation variance value of points in the neighborhood of all rectangles (k=0,1,2...18) of point p _i , and obtain the elevation variance value set Compare the elevation variance values in the set Ψ respectively, and record the three values with the smallest, second smallest, and third smallest elevation variance values and and the corresponding rectangle numbers k ₁ , k ₂ , k ₃ , if the following two equations are satisfied at the same time, the current point p _i is recorded as the road center point.

$\left\{\begin{array}{}\end{array}\begin{array}{c}\frac{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; \&psi;</span>}\\ \mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}<span\; class="notranslate">\; \&le;</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}\end{array}\right.$

In the formula, ψ is the road roughness threshold set by the user, which is taken as 0.5 in the present invention, and μ is the road turning threshold set by the user, which is taken as 1 in the present invention, that is, the maximum turning angle is 20 degrees.

(4) Repeat the above steps until the judgment of all points is completed, and record the road center point set as P _k .

6. Elimination of non-road points

According to the road candidate point set P _c and the road center point set P _k calculated in steps 4 and 5 respectively, the road candidate point set may contain a large number of points in flat areas such as parking lots, sports fields, and roofs. These points are divided by the following method remove.

(1) Using a method similar to step 2, divide the point sets P _c and P _k into standard grids, and use codes 1 and 2 to indicate whether the point type is a road candidate point or a road center point.

(2) For each point in the point set _Pc , search for points with an attribute of 2 within the r/2 range around it, if there is, keep the point, otherwise it is judged as a non-road point. The parameter r/2 represents the ratio of 1/2 of the road width to the grid spacing, and r/2 is taken as 2 in the present invention.

Of course, the above-mentioned embodiments are only to illustrate the technical conception and characteristics of the present invention, and its purpose is to enable those skilled in the art to understand the content of the present invention and implement it accordingly, and not to limit the protection scope of the present invention. All equivalent changes or modifications made according to the spirit of the main technical solutions of the present invention shall fall within the protection scope of the present invention.

Claims (6)

1. A road data extraction method based on rectangular neighborhood analysis, characterized in that it comprises: The first step: data preprocessing, calculate point cloud spacing, establish regular grid organization, and then use fast mathematical morphology filtering method based on regular grid to filter, and eliminate non-ground points; Step 2: Screen candidate road points, preliminarily determine the road point set, narrow the judgment range, search for points within its neighborhood for each point, and compare their elevation values. If they are all within the threshold range, they will be used as road candidate points ; Step 3: Extract the road centerline, remove the flat area point set, calculate the elevation variance value within the rectangular range of different azimuths for each road candidate point, and judge whether the current point is the road center point according to the difference between the variance values ; The fourth step: road connection, connect and smooth the discontinuous road. 2. a kind of road data extraction method based on rectangular neighborhood analysis according to claim 1, is characterized in that: described first step comprises: (1) First read the point cloud data range (xmin, xmax, ymin, ymax), and count the number of points pt_Num; (2) Calculate point cloud spacing: $\mathrm{<span\; class="notranslate">\; pt</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; dis</span>}<span\; class="notranslate">\; =</span>\frac{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; max</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; min</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}\mathrm{<span\; class="notranslate">\; max</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}\mathrm{<span\; class="notranslate">\; min</span>}<span\; class="notranslate">\; )</span>}{\sqrt{\mathrm{<span\; class="notranslate">\; pt</span>}<span\; class="notranslate">\; \_</span>\mathrm{<span\; class="notranslate">\; num</span>}}}$ (3) According to the point cloud spacing pt_dis and calculation needs, establish a regular grid covering the point cloud range (xmin, xmax, ymin, ymax); (4) distribute the airborne LiDAR data into individual grids according to the coordinates of the points, $\left\{\begin{array}{c}\mathrm{XAxisNum}=\mathrm{int}((x_{\max }-x_{\min })/d)+1\\ \mathrm{YAxisNum}=\mathrm{int}(({\mathrm{the\; y}}_{\max }-{\mathrm{the\; y}}_{\min })/d)+1\end{array}\right.$ (Formula 2) GridNum＝XAxisNum*YAxisNum (Formula 3) XAxisNum and YAxisNum are the number of grids in the x and y directions, and GridNum is the total number of grids; The grid of any point p _i ( _xi , y _i ) in the airborne LiDAR data is: $\left\{\begin{array}{c}{x}_i=\mathrm{int}((x_i-x_{\min })/d)\\ Y_i=\mathrm{int}(({\mathrm{the\; y}}_i-{\mathrm{the\; y}}_{\min })/d)\\ G_i=\mathrm{int}((x_i-x_{\min })/d)+x_i\&times;\mathrm{int}(({\mathrm{the\; y}}_i-{\mathrm{the\; y}}_{\min })/d)\end{array}\right.$ (Formula 4) Among them, X _i is the number of grids in the X direction where P _i is located, Y _i is the number of grids in the Y direction where P _i is located, and G _i is the number of grids recorded in the one-dimensional array where P _i is located; Carry out non-ground filtering on the airborne LiDAR data, retain ground point data including road points, and use regular grid as the unit of mathematical morphology filtering. The result of grid erosion is the minimum elevation value in the w×w grid, and the expansion The result is the maximum elevation value in the w×w grid, where w represents the size of the window, and the operation of formula 5 is performed: $\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) where A represents the target set, B represents the structural element, AΘB represents corrosion, In the formula, zi _′ and zi _′ ′ are the elevation values after grid cell erosion and expansion respectively, and Grid _i (z _min ) is the minimum elevation value in the grid. 3. a kind of road data extraction method based on rectangular neighborhood analysis according to claim 2, is characterized in that: described first step further comprises: (1) Calculate the number of grids from the grid spacing, and allocate memory space for the grid, store the points one by one in the corresponding grid unit according to the plane coordinates, and compare with the minimum elevation value saved in the grid unit , record the minimum elevation value after comparison; (2) Traverse each grid, compare the minimum elevation value in the surrounding w×w grids, use the first formula of formula 5, take the minimum value as the elevation value of the current grid after corrosion, and set the maximum value of the grid The elevation value is set to the elevation value after corrosion; (3) Traversing through each grid, comparing the corrosion elevation values in the surrounding w×w grids, using the second formula of formula 5, taking the maximum value as the elevation value after the expansion of the current grid, and making the grid minimum The elevation value is set to the inflated elevation value; (4) Traversing each grid, calculating the difference between the elevation value of each point in the grid and the inflated elevation value, when its absolute value is less than the set threshold, the point is a ground point, otherwise it is a non-ground point. 4. a kind of road data extraction method based on rectangular neighborhood analysis according to claim 3, is characterized in that: described second step comprises: (1) For each point p _i (x _i , y _i ), calculate its grid by formula 4; (2) Calculate the search radius, $r=\frac{\mathrm{the\; road}\_\mathrm{width}}{d};$ (3) Let the point set in the rectangular neighborhood of point p _i (xi _, y _i ) be P _{i_near} ＝{p _j (x _j ,y _j )|x _j ∈(xi _-r , _xi +r),y _j ∈(y _i -r,y _i +r)} (where j is the point The point number or subscript is different from the current point pi), so the point set P _{i_near} in the surrounding grid is searched with the current point pi as the center; (4) Calculate the sum ΣH _i of the absolute value of the elevation difference between each point in P _{i_near} and p _i (xi _, y _i ), if ΣH _i is less than the threshold, the current point is a road candidate point; $\mathrm{\&Sigma;}{h}_i=\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\&Sigma;}}}|z_j-z_i|$ (Formula 6) z _j is the elevation value of the jth point in P _{i_near} , z _i is the elevation value of p _i ( _xi , y _i ), n is the total number of points in P _{i_near} ; (5) Repeat the above steps until the judgment of all points is completed, and the reserved road candidate points are set as P _c . 5. a kind of road data extraction method based on rectangular neighborhood analysis according to claim 4, is characterized in that: described 3rd step comprises: (1) For each road candidate point p _i , calculate the rectangular range points p _i,r1 , p _i,r2 , p _i,r3 , p _i,r4 according to formula 7, that is, the four vertices of the rectangle, p _i,r1 It coincides with the current point, and the other points are arranged counterclockwise, where dθ is the rotation angle of the rectangle, and divides 360 degrees into n intervals; $\begin{array}{c}\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{{\mathrm{<span\; class="notranslate">\; l</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{{\mathrm{<span\; class="notranslate">\; l</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; w</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; l</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>\\ {\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; cos</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; d\&theta;</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&pi;</span>}<span\; class="notranslate">\; )</span>\end{array}\right.\end{array}$ k＝0,1,2...18 (Formula 7) Among them, p _{i, r1} (x) is the x coordinate value of p _{i, r1 , p i, r1} (y) is the y coordinate value of p _{i, r1} , and the rest are the same, l is the length of the rectangle, w is the width of the rectangle, dθ is the rotation angle interval of the rectangle, and k is the number of rotations or the number of iterations; (2) According to the neighborhood judgment method in step 4, search for the neighboring point P _{i_near} of point p _i , and then judge whether P _{i_near} is formed by p _i,r1 , p _i,r2 , p _i,r3 , p _i,r4 Within the rectangular range, save the point set P _{i_r,k} within the range, and calculate its elevation variance Elevation mean: Where z _j is the elevation value of the jth point in the point set P _{i_r,k} , and m is the number of points in the point set P _{i_r,k} ; Elevation variance value: ${\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 the above two steps, calculate the elevation variance value of points in the neighborhood of all rectangles (k=0,1,2...18) of point p _i , and obtain the elevation variance value set Compare the elevation variance values in the set Ψ respectively, and record the three values with the smallest, second smallest, and third smallest elevation variance values and and the corresponding rectangle numbers k ₁ , k ₂ , k ₃ , if the following two equations are satisfied at the same time, the current point p _i is recorded as the road center point; $\left\{\begin{array}{}\end{array}\begin{array}{c}\frac{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; min</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; \&psi;</span>}\\ \mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}<span\; class="notranslate">\; \&le;</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}\end{array}\right.$ where ψ is the road roughness threshold set by the user, and μ is the road turning threshold set by the user; (4) Repeat the above steps until the judgment of all points is completed, and record the road center point set as P _k . 6. a kind of road data extraction method based on rectangular neighborhood analysis according to claim 5, is characterized in that: described 4th step comprises: According to the calculated road candidate point set P _c and road center point set P _k , the road candidate point set may contain a large number of points in flat areas such as parking lots, sports fields, and roofs. These points are eliminated by the following methods: (1) Divide the point sets P _c and P _k into standard grids, and use codes 1 and 2 to indicate whether the point type is a road candidate point or a road center point; (2) For each point in the point set _Pc , search for a point with an attribute of 2 within the r/2 range around it, if there is, keep the point, otherwise it will be judged as a non-road point, and the parameter r/2 represents the width of the road The ratio of one-half of to the grid spacing.