CN102103202A - Semi-supervised classification method for airborne laser radar data fusing images - Google Patents

Abstract

The invention relates to the technical field of airborne laser radar data processing, in particular to a semi-supervised classification method for airborne laser radar data with image fusion. Based on the semi-supervised concept, the present invention uses the rough classification results of point cloud data to extract high-precision training sample data for the classification of high-resolution images. The post-processing process removes false building points based on the target complexity and refines the classification results. , and integrate the multiple features of the LiDAR point cloud for cross-validation, and finally realize the fine classification of the airborne LiDAR data, which is a fusion classification method with high reliability and high classification accuracy. The present invention achieves a good effect of classifying point cloud areas with high vegetation and low vegetation without using near-infrared data.

Description

A semi-supervised classification method for airborne lidar data based on image fusion

technical field

The invention relates to the technical field of airborne laser radar data processing, in particular to a semi-supervised classification method for airborne laser radar data with image fusion.

Background technique

Airborne LiDAR is a new type of active aerial remote sensing earth observation technology, which can directly obtain the spatial three-dimensional point cloud information of the target. With the acceleration of urbanization, it is of great significance to use airborne LiDAR technology to achieve high-efficiency extraction of urban ground object information. The most basic and key technology is the classification of LiDAR point cloud data. Points classified as bare ground are used for the generation of digital ground models, providing basic data for topographical mapping, engineering surveying, environmental planning, etc.; points classified as buildings and vegetation can be used to improve the accuracy of DTM models and 3D digital cities. Building model reconstruction, urban green space research, etc. However, since the point cloud data provided by LiDAR cannot directly obtain the semantic information (material and structure, etc.) Handling becomes more difficult. The existing airborne LiDAR point cloud processing algorithm is not capable of automatic interpretation of complex urban area geomorphology. In practice, the time spent on quality control and manual operation occupies a considerable proportion of the entire data processing time. Therefore, it is necessary to design an automatic, efficient, and robust LiDAR point cloud classification and modeling algorithm.

Existing studies have shown that due to the lack of corresponding texture and semantic information, the classification and intelligent processing of ground objects using airborne laser scanning data alone has great limitations, and the classification accuracy for complex scenes is not high, which cannot meet the actual classification requirements. Handle application requirements.

Contents of the invention

Aiming at the limitations of the above-mentioned classification of single-source remote sensing data, especially the lack of classification accuracy of single laser radar data, the purpose of the present invention is to provide a semi-supervised classification method for airborne laser radar data that fuses images, using semi-supervised method to train Samples are classified by fusion of high-resolution images and airborne LiDAR data, and finally achieve the purpose of high-precision classification of point cloud data.

To achieve the above object, the present invention adopts the following technical solutions:

The original lidar data denoising step, which uses the K nearest neighbor sphere denoising algorithm to remove the noise existing in the point cloud, and interpolates to generate a digital surface model DSM;

A high-precision digital ground model DEM generation step, this step uses the iterative triangulation network progressive encryption filtering method to obtain the laser radar points on the ground for the denoised laser radar data, and interpolates to generate high-precision digital ground model DEM data;

nDSM data generation step, this step subtracts original DSM data and DEM data, obtains nDSM data;

Coarse classification of lidar data. In this step, for the non-ground point set obtained by iterative triangulation network filtering, the point cloud is first segmented by elevation information, and then the medium-high vegetation and buildings in the point cloud data are obtained by using local attribute estimation and elevation constraints. The initial category information of two categories such as objects;

The image classification training sample extraction step based on LiDAR data. In this step, for the point cloud after rough classification, the elevation value is first assigned according to the category information, and the corresponding grid data is generated by gridding; secondly, the seed point is manually selected and grown through the seed. The sample area is obtained by automatic growth using the above method; the sample information of high vegetation, buildings, and bare ground is semi-automatically obtained through the above method, and used as a high-precision training sample for image classification;

Combined with the nDSM mask classification step, this step uses nDSM data to generate a mask for the area with an elevation of 0 on the registered high-resolution image;

The step of removing false building points after classification is because the classification results based on the joint nDSM mask relying solely on spectral information will lead to misclassification of building categories due to excessive reliance on spectral information, so the use of shape index and Complexity calculation removes non-building category data, and removes laser points that are misjudged as building points;

The step of laser point classification based on the image classification result and point cloud multiple feature cross-validation, this step first uses the image classification result obtained by the joint nDSM mask classification step and the classified pseudo building point removal step to classify the point cloud Assignment; secondly, use the multiple features of the point cloud (intensity mean, dispersion) and DEM data to re-verify the classification assignment results, correct the misclassified points of the point cloud classification, and finally classify the point cloud data into bare ground, low vegetation, high 4 categories such as vegetation and buildings.

The step of generating the high-precision digital terrain model DEM further includes the following sub-steps:

① Carry out median filter processing on the original data to eliminate extremely low points in the data (noise points with very low elevation);

② Construct the outer rectangle of the data, the elevation values of the four vertices of the outer rectangle are set according to the nearest neighbor criterion, and then triangulate the outer rectangle, and use it as the initial terrain surface model;

③Grid organization of the data, the grid should be slightly larger than the size of the largest building, where the lowest point in each grid is the initial ground point, and the selected initial ground point is added to the triangulated irregular network (TIN);

④ Calculate the distance from each point to the triangle where it is located and the angle between it and the three vertices of the triangle. If the calculated value is less than the preset threshold condition, add it to the irregular triangulation;

⑤ Repeat ④ until no new points are added to the irregular triangulation;

⑥ Generate DEM by interpolation.

The classification step of the joint nDSM mask further comprises the following sub-steps:

① Based on nDSM information, divide the high-resolution image registered with the point cloud into a high-resolution area and a flat area, and use the building and high-vegetation samples obtained in the LiDAR data-assisted image classification training sample extraction step, in Areas whose elevation is 0 or whose elevation is within a specified threshold are classified using maximum likelihood estimation;

②Based on nDSM information, within the image range where the elevation area is 0 or lower than a certain threshold, artificially select surface vegetation samples, and at the same time use the high-precision surface samples obtained in the LiDAR data-assisted image classification training sample extraction step, using Maximum likelihood estimation was performed for classification to obtain high-accuracy classifications of bare earth surfaces with spectral information for green low vegetation categories.

The present invention has the following advantages and positive effects:

1) The accuracy of the samples used for image classification obtained by the present invention through LiDAR rough classification is very high.

2) On the premise of not using near-infrared data, the present invention achieves a good effect of classifying point cloud areas with high vegetation and low vegetation.

Description of drawings

Fig. 1 is a schematic diagram of the relationship between noise points and surrounding neighborhoods in the present invention.

Fig. 2 is a schematic diagram of iterative triangulation filtering in the present invention.

Fig. 3 is a flow chart of point cloud fusion classification based on semi-supervised classification and high-resolution images provided by the present invention.

Detailed ways

The present invention provides a semi-supervised classification method for airborne lidar data that fuses images. Based on the semi-supervised concept, the rough classification results of point cloud data are used to extract high-precision training sample data for the classification of high-resolution images, and The multiple features of the LiDAR point cloud are fused for cross-validation, and finally the fine classification of the airborne LiDAR data is realized.

Below in conjunction with accompanying drawing, the present invention will be further described with specific embodiment:

A semi-supervised classification method for airborne laser radar data provided by the invention comprises the following steps:

(1) Original lidar data denoising:

If the point cloud data has extremely low points and air points that are significantly lower or higher than the surrounding environment, it will greatly affect the accuracy of the post-processing algorithm, so these noise points should be removed before data processing. This method detects and removes the noise in the point cloud by establishing a K-nearest neighbor sphere. First, the data point set is divided into spatial grids. It is assumed that there is a space sphere, and the current measurement point is used as the center of the sphere, and the radius is respectively taken from the measurement point to the cube grid where it is located. 6 face distance. Take the space sphere with the smallest radius, and perform K-nearest neighbor search in the grid that interferes with it. If the established search termination principle is satisfied, the search is terminated; otherwise, take a space sphere with a larger radius to establish the K of the point to be determined. adjacent to the ball. In the process of noise processing on the point cloud, it mainly depends on the distance between the point to be determined and the point in the established K-nearest neighbor sphere to determine whether the point to be determined is noise. (As shown in Figure 1).

(2) High-precision DEM generation

The ground point set is obtained by iterative triangulation filtering, and the grid DEM is obtained by interpolating the ground point set. The key step is iterative triangulation filtering, and the steps are: ① Carry out median filtering on the original data to eliminate extremely low points (noise points with very low elevation) in the data (this step has been completed in the first step); ② Construct the data The surrounding rectangle of the surrounding rectangle, the elevation values of the four vertices of the surrounding rectangle are set according to the nearest neighbor criterion, and then the surrounding rectangle is triangulated and used as the initial terrain surface model (a in Figure 2); ③ for The data is organized in a grid, and the grid should be slightly larger than the size of the largest building. The lowest point in each grid is the initial ground point, and the selected initial ground point is added to the triangulated irregular network (TIN); ④ Calculation The distance from each point to the triangle where it is located and the angle between it and the three vertices of the triangle, if the calculated value is less than the preset threshold condition, it will be added to the irregular triangular network (as shown in Figure 2 b); ⑤Repeat ④ until no new points are added to the irregular triangulation network (c in Figure 2); ⑥ Interpolate to generate DEM.

The principle of iterative triangulation progressive encryption filtering algorithm is shown in Figure 2.

(3) nD SM (normalized digital surface model) data generation

The normalized digital surface model (nDSM), that is, the data obtained after the algebraic difference operation between the digital surface model (DSM) and DEM (ie DSM-DEM), nDSM can directly reflect the height information of ground objects, and alleviate the impact of terrain fluctuations. Elevation effects caused by ground features. This kind of treatment is usually suitable for the terrain coverage areas with drastic changes of topography.

According to the above principles, the denoised point cloud data is interpolated to generate a grid DSM, and the ground point data obtained by iterative triangulation filtering is used to interpolate to generate a grid DEM, and then the two data are subtracted to obtain nDSM data.

(4) Automatic rough classification of point cloud data

Segment the non-ground point data obtained by filtering, and segment the point cloud in the same plane into the same segment. Since the building points are obviously higher than the surrounding points, and most of the roof surfaces of the buildings change less, so for the filtered laser points, the elevation difference between the segmentation segment and the surrounding ground of the segment, and the segmentation segment The degree of local variation of the described surface (see Sections 4.1 and 4.2) is used to identify whether the segment belongs to the building area. On the basis of removing the building area, the elevation threshold is used to identify the point set of the vegetation area. After the automatic coarse classification process is completed, the original LiDAR point cloud will be roughly classified into point sets of three target categories: ground, buildings, and vegetation.

4.1. Normal Estimation

为p的邻域的质心，即Let the neighborhood of the sample point p be (p1, p2, ..., pk),

is the centroid of the neighborhood of p, ie

$\stackrel{\‾}{p}=\frac{1}{k}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{p}_i$ Formula 1

Since each point in the point cloud has three components x, y, and z, the covariance matrix of point p is a 3×3 matrix, which can be defined as

$C={\left[\begin{array}{c}{p}_1-\stackrel{\‾}{p}\\ \mathrm{\Λ}\\ p_k-\stackrel{\‾}{p}\end{array}\right]}^T\left[\begin{array}{c}{p}_1-\stackrel{\‾}{p}\\ \mathrm{\Λ}\\ p_k-\stackrel{\‾}{p}\end{array}\right]$ Formula 2

By accumulating sample points in the neighborhood of point p to the centroid In the square distance of the three component directions, the covariance matrix C can describe the statistical characteristics of the distribution of these sample points.

Consider the eigenvector problem

C · v _j = λ _j · v _j formula 3

Since C is a symmetric positive semidefinite matrix, all the eigenvalues should be real values, and the eigenvectors _vj constitute the vertical coordinate system, and correspond to the three main components of the sample point set in the neighborhood respectively. What the eigenvalue measures is the change of the sample point p _i (i=1, 2, . . . , k) in the neighborhood along the direction of the corresponding eigenvector.

是这样的一个平面，它通过质心点

且使得点p的邻接点到达该平面的平方距离和最小。也可以认为平面T(x)是曲面在点p处的切平面的逼近。因此，向量v₀可近似的看成是逼近点p处的曲面法线n_p，向量v1和v2则生成了曲面在点p处的切平面。Assuming that λ ₀ ≤λ ₁ ≤λ ₂ , the following conclusions can be drawn: the plane

is a plane that passes through the centroid point

And make the sum of the square distances of the adjacent points of point p to the plane minimum. The plane T(x) can also be thought of as an approximation of the tangent plane of the surface at point p. Therefore, the vector v ₀ can be approximated as the surface normal n _p at the point p, and the vectors v1 and v2 generate the tangent plane of the surface at the point p.

4.2 Curvature Estimation

Based on the 4.1 normal estimation method, the curvature of the point on the surface is estimated by using the normal of the sample point in the neighborhood. Assuming λ ₀ ≤λ ₁ ≤λ ₂ , the measurement is the change of the neighborhood of point p along the normal direction of the surface, that is, the degree of deviation of adjacent points from the tangent plane Tp. The overall degree of deviation of adjacent points, that is, the sum of the squares of distances between adjacent points p _i and the centroid can be given by the following formula:

$\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{|p_i-\stackrel{\‾}{p}|}^2={\mathrm{\λ}}_0+{\mathrm{\λ}}_1+{\mathrm{\λ}}_2$ Formula 4

Therefore, under the condition that the neighborhood size is k, the surface change at point p can be defined as

${\mathrm{\σ}}_k\left(p\right)=\frac{{\mathrm{\λ}}_0}{{\mathrm{\λ}}_0+{\mathrm{\λ}}_1+{\mathrm{\λ}}_2}$ Formula 5

If σ _k (p)=0, it means that all points are on the tangent plane Tp. When the changes of these points are the same in all directions, the surface change reaches 1/3 of its maximum value. The surface variation will vary with the size of the selected neighborhood. When the value of the neighborhood is larger, the estimated surface changes more sharply, and when the value of the neighborhood is smaller, the change of the surface is flatter.

(5) Image classification training sample extraction based on LiDAR data

After the rough classification of the point cloud, according to the category information, the elevation is assigned and then the grid is generated to generate the corresponding raster data. The raster data is displayed through the elevation, and a point in a certain category area in the raster image is manually interpreted and selected as the seed point , perform the seed point growth method, and reach any pixel in the area through eight directions of up, down, left, right, upper left, lower left, upper right, and lower right to obtain the sample area. Through this method, the sample information of ground category, building category, and high vegetation category will be accurately collected.

(6) Classification of joint nDSM masks.

a: Based on nDSM information, divide the high-resolution image registered with the point cloud into a high-resolution area and a flat area. Using the buildings and high-vegetation samples obtained in step (5), when the elevation is greater than 0 or the elevation is greater than the specified Image regions within the threshold are classified using maximum likelihood estimation. Usually when image data is only used, artificial ground and buildings, grassland and tall vegetation tend to interfere with each other, resulting in misclassification of categories. The accuracy of vegetation classification results has been improved.

: Based on nDSM information, within the image range where the elevation area is 0 or lower than a certain threshold, artificially select surface vegetation samples, and use the high-precision surface samples obtained in step (5) to classify using maximum likelihood estimation to obtain high Accuracy Classification of Bare Surface and Spectral Information for Green Low Vegetation Classes.

(7) Removal of false building points after classification

Since classification based solely on spectral information can easily lead to building category misclassifications, these misclassifications need to be corrected.

Buildings can often be interpreted as geometrically uniform, meaningful objects. Therefore, the non-building points misclassified as buildings can be removed according to the geometric characteristics of the target, such as the shape index (such as formula 6). The size of building objects is usually expressed by area or perimeter. The area of its area can be represented by the number of pixels that make up the area; its perimeter is obtained by calculating the length of the boundary line, that is, calculating the length of the boundary curve. The shape index is defined as Equation 6, where S is the target area and P is the perimeter. Non-building objects usually have a smaller shape index, and the more complex the shape of a building object, the larger the value of the shape index.

$I=\frac{\sqrt{S}}{P}$ Formula 6

(8) Laser point classification based on cross-validation of image classification results and point cloud multiple features

The ground object category map obtained after processing through (6) and (7) steps is used to reassign the category information of the point cloud data after rough classification, but relying solely on the category information provided by the category map will lead to the loss of buildings and high vegetation. Misclassification occurs in some areas, so when using the category map to classify the point cloud, it is necessary to use various features of the point cloud data to re-verify the category information of the point cloud, that is, by comprehensively considering features such as elevation, intensity, and dispersion to determine the class attribute of the final laser point.

a. Using elevation features to identify misclassified points in the feature category map

The elevation of buildings and high vegetation points should theoretically be much higher than their surrounding ground points, so it can be used as a criterion to distinguish ground points that are misclassified as buildings or high vegetation categories. Therefore, when a point is judged to be a building or a high vegetation category according to the category information of the category map, the coordinate value of the point should be used to interpolate in the DEM data, and according to the difference between the elevation obtained at the point and the elevation of the point itself Whether the difference is greater than a preset threshold is used to determine whether the point is a ground point. If it is greater than the specified threshold, then the point is far from the ground, which conforms to the judgment theory, and the original category information can be kept unchanged. On the contrary, the category information of this point should be a ground point, and the category of this point can be directly assigned as a bare ground object point.

b. Use the intensity feature to identify misclassified points in the object category map

The building intensity information is completely different from the high vegetation intensity. The intensity is used as a distinguishing standard, and the right and wrong are corrected for the high vegetation point category information of the building class. In the category map, the point cloud intensity values in the building area and high vegetation area are counted to obtain the respective average values BI and TI, and the edge points of the building area are traversed. If the intensity of a building category point is not within the set If it is within the specified building threshold range (usually the threshold is BI±BI/4), and at the same time it is within the high vegetation intensity threshold range (usually the threshold is TI), then the point is likely to be a high vegetation point, and the laser should be used Points are divided into high vegetation points, otherwise keep the original category information unchanged.

c. Use the dispersion feature to identify misclassified points in the object category map

The spatial dispersion of laser point cloud is an important clue to distinguish buildings and vegetation. Since buildings and bare ground are usually composed of planes, its degree of discreteness in space can be considered to be distributed along a two-dimensional surface (not necessarily horizontal), while points on trees are scattered in all directions in space. This discreteness can be obtained through The eigenvalues of the discrete matrix to analyze.

The specific method is as follows: search all adjacent points in the neighborhood of the laser point, and establish a 3×3 discrete matrix of the laser point in space, see formula 7:

$S_j=\underset{i=0}{\overset{\mathrm{no}}{\mathrm{\Σ}}}\left(v_i^T{v}_i\right)(j=\mathrm{0,1,2}.....m)$ Formula 7

为v_i的转置矩阵，M为激光点数。Where S _j is the 3×3 discrete matrix of the jth point, n is the number of adjacent points in the neighborhood of the jth point, and v _i is the spatial coordinate of the ith adjacent point of the jth point v _i =(x _i , y _i , z _i ),

is the transposition matrix of v _i , and M is the number of laser points.

Singular value decomposition of the discrete matrix S _j can obtain the three eigenvalues of the point matrix, and arrange the eigenvalues from small to large. Set three categories:

a. One eigenvalue is much larger than the other eigenvalue, then the point is marked as a plane

b. One eigenvalue is much larger than the other two eigenvalues, then the point is marked as an edge class

c. If the three eigenvalues are large enough, they are marked as spatially discrete.

Based on the above three criteria, traverse the vegetation points and use the dispersion to correct the misclassified building points as vegetation.

The above embodiments are only for the purpose of illustrating the present invention, rather than limiting the present invention. Those skilled in the relevant technical fields can also make various changes or modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent All technical solutions fall within the protection scope of the present invention.

Claims (3)

1. a kind of airborne lidar data semi-supervised classification method of fusion image, is characterized in that, comprises the following steps: The original lidar data denoising step, which uses the K nearest neighbor sphere denoising algorithm to remove the noise existing in the point cloud, and interpolates to generate a digital surface model DSM; A high-precision digital ground model DEM generation step, this step uses the iterative triangulation network progressive encryption filtering method to obtain the laser radar points on the ground for the denoised laser radar data, and interpolates to generate high-precision digital ground model DEM data; nDSM data generation step, this step subtracts original DSM data and DEM data, obtains nDSM data; Coarse classification of lidar data. In this step, for the non-ground point set obtained by iterative triangulation network filtering, the point cloud is first segmented by elevation information, and then the medium-high vegetation and buildings in the point cloud data are obtained by using local attribute estimation and elevation constraints. The initial category information of two categories such as objects; The image classification training sample extraction step based on LiDAR data. In this step, for the point cloud after rough classification, the elevation value is first assigned according to the category information, and the corresponding grid data is generated by gridding; secondly, the seed point is manually selected and grown through the seed. The sample area is obtained by automatic growth using the above method; the sample information of high vegetation, buildings, and bare ground is semi-automatically obtained through the above method, and used as a high-precision training sample for image classification; Combined with the classification step of nDSM mask, this step uses nDSM data to generate a mask on the registered high-resolution image to apply mask processing to the planar area whose elevation is lower than the specified elevation threshold and the high area whose elevation is higher than the specified elevation threshold, The step of removing false building points after classification is because the result of joint nDSM mask classification based on spectral information alone will lead to misclassification of building categories due to purely relying on spectral information, so the use of shape index and complexity Calculate and remove non-building category data, and remove laser points that are misjudged as building points; The step of laser point classification based on the image classification result and point cloud multiple feature cross-validation, this step first uses the image classification result obtained by the joint nDSM mask classification step and the classified pseudo building point removal step to classify the point cloud Assignment; secondly, use the multiple features of the point cloud (intensity mean, dispersion) and DEM data to re-verify the classification assignment results, correct the misclassified points of the point cloud classification, and finally classify the point cloud data into bare ground, low vegetation, high 4 categories such as vegetation and buildings. 2. the airborne lidar data semi-supervised classification method of fusion image according to claim 1, is characterized in that: The step of generating the high-precision digital terrain model DEM further includes the following sub-steps: ① Carry out median filter processing on the original data to eliminate extremely low points in the data (noise points with very low elevation); ② Construct the outer rectangle of the data, the elevation values of the four vertices of the outer rectangle are set according to the nearest neighbor criterion, and then triangulate the outer rectangle, and use it as the initial terrain surface model; ③Organize the data in a grid. The grid should be slightly larger than the size of the largest building. The lowest point in each grid is the initial ground point, and the selected initial ground point is added to the triangulated irregular network (TIN); ④ Calculate the distance from each point to the triangle where it is located and the angle between it and the three vertices of the triangle. If the calculated value is less than the preset threshold condition, add it to the irregular triangulation; ⑤ Repeat ④ until no new points are added to the irregular triangulation; ⑥ Generate DEM by interpolation. 3. the airborne lidar data semi-supervised classification method of fusion image according to claim 1 or 2, is characterized in that: The classification step of the joint nDSM mask further comprises the following sub-steps: ① Based on nDSM information, divide the high-resolution image registered with the point cloud into a high-resolution area and a flat area, and use the building and high-vegetation samples obtained in the LiDAR data-assisted image classification training sample extraction step, in Areas whose elevation is greater than 0 or whose elevation is greater than a specified threshold are classified using maximum likelihood estimation; ②Based on nDSM information, within the image range where the elevation area is 0 or lower than a certain threshold, artificially select surface vegetation samples, and at the same time use the high-precision surface samples obtained in the LiDAR data-assisted image classification training sample extraction step, using Maximum likelihood estimation was performed for classification to obtain high-accuracy classifications of bare earth surfaces with spectral information for green low vegetation categories.

Images