Point cloud filtering method based on Gaussian mixture model with multi-level curvature constraints

Abstract

The present invention proposes a Gaussian mixture model point cloud filtering method based on multi-level curvature constraints. The method mainly includes three steps: (1) According to the law of large numbers, first calculate the average distance from each measurement point to its K nearest neighbor points, and then eliminate various abnormal measurement values ​​in the point cloud based on the principle that the normal measurement value will not be greater than M times the standard deviation. (2) Secondly, the ground seed points obtained by the local minimum are interpolated using the thin plate spline function (TPS), and then the points whose height difference with the interpolation surface is greater than the curvature threshold are eliminated. This process is to iteratively approximate a reference surface by continuously reducing the interpolation grid units and gradually increasing the curvature threshold. (3) Finally, calculate the relative height difference of each measurement value, apply the binary Gaussian mixture model (GMM) to classify the ground object and the ground, and then use the expectation maximization algorithm (EM) to obtain the optimal solution of the model parameters, and judge whether the point belongs to the ground object based on the calculated posterior probability to complete the point cloud filtering.

Description

Technical Field

[0001] The present invention belongs to the fields of laser radar application, computer vision and geographic information system, and mainly includes aspects such as point cloud denoising, surface grid construction and point category determination according to posterior probability in the ground filtering process of airborne LiDAR point cloud. Background Art

[0002] The airborne LiDAR system measures the time of the laser echo to obtain the distance from the target to the center of the instrument. As long as the geographic reference system is given, the system can quickly obtain high-precision point cloud data. Compared with the traditional method of obtaining point clouds by photogrammetry, the advantages of the LiDAR system are shorter processing time and higher measurement accuracy. At present, LiDAR technology has been extended to many fields, such as digital terrain model construction, three-dimensional building model construction, road surface extraction, and estimation of forest station parameters. In many applications, it is an essential task to distinguish (filter) ground measurements from non-ground measurements in the original point cloud. Similarly, the construction of an accurate digital elevation model (DEM) also requires the removal of non-ground points in advance, and then the raster grid is calculated by the geoscientific interpolation method.

[0003] Airborne LiDAR technology has revolutionized the way of obtaining high-precision DEM, mainly because of its high efficiency, high accuracy and low price in obtaining point clouds. Compared with digital photogrammetry technology, airborne LiDAR provides a high-precision terrain mapping solution for large areas. Many researchers have focused on how to use LiDAR point cloud data to obtain terrain information. In order to obtain DEM, a necessary step is to filter out non-ground point clouds, including measurements from bridges, buildings, vegetation, etc. It is worth noting that although airborne LiDAR technology can obtain a large amount of data in a short time, its filtering and manual quality control account for 60%-80% of the entire post-processing time. Therefore, how to reduce the time complexity of the point cloud data processing algorithm and reduce its parameter adjustment time is the main problem of current research work.

[0004] Dozens of ground filtering algorithms have been proposed in the past two decades. However, filtering is still a major challenge in point cloud data processing in terms of processing time and processing accuracy. In 2004, Sithole and Vosselman conducted comprehensive comparative experiments on 8 classic filtering algorithms on rural and urban test area data. They pointed out that all algorithms can perform well on relatively flat terrain. However, when the terrain becomes steep and accompanied by complex objects (such as bridges, low shrubs, etc.), the filters will fail and perform unsatisfactory. Overall, surface-based filters perform better in most scenes. This is mainly because this type of filter uses more point cloud data context information than other types of filters. Since then, the 15 test samples they selected in rural and urban areas have been used as standard data sets for testing ground filtering algorithms.

[0005] The ground filtering algorithms of airborne LiDAR point clouds can be mainly divided into four categories: slope-based filters, mathematical morphology-based filters, surface interpolation-based filters, and segmentation strategy-based filters. Among them, filters based on mathematical morphology and surface interpolation are widely used methods. In particular, the commercial software TerraScan integrates the classic progressive triangular densification filter (PTD) algorithm, which can obtain terrain measurements with relatively high accuracy. However, a major drawback of this algorithm is that it requires complex parameter adjustments, and because the triangular facets will be severely deformed on special terrains, it cannot meet the ground filtering accuracy requirements of any terrain. The slope-based filter mostly assumes that the slope between terrain point measurements will be smaller than that of non-terrain points, and filters out non-terrain points by giving a slope threshold. Since the assumptions are mostly applicable to relatively flat terrains, and most terrain conditions are complex and changeable, the application of this type of method in practical work is limited. The basic idea of ​​the filter based on the segmentation strategy is to cluster first and then classify. The clustering algorithm such as the region growing method is used to merge the measurement points according to the homogeneity principle, and then each segment is classified into ground segment and non-ground segment given an angle or distance threshold to achieve the filtering purpose. However, this method usually requires the selection of a suitable clustering method to ensure the internal homogeneity of the segment, and is accompanied by complex parameter setting problems, so it is rarely used in filtering work.

[0006] Most of the current filtering algorithms require complex parameter adjustments based on the specific test area to achieve the required accuracy. However, this is not very friendly to many surveying and mapping industry workers, because they need to have a certain understanding of the algorithm mechanism and require a lot of experimental adjustment parameters. For example, when using a filter based on mathematical morphology, the algorithm often assumes that the slope s is a constant. But this assumption cannot be met in actual work, because the terrain slope is always changing under the action of internal and external forces on the surface. Using a constant s value may lead to errors in the calculation of the height threshold dT, and ultimately produce a large number of omission errors or commission errors, which will lead to catastrophic errors in the subsequent DEM construction. Summary of the invention

[0007] The purpose of the present invention is to propose a novel Gaussian mixture model point cloud filtering method based on multi-level curvature constraints to address the shortcomings of existing ground filtering algorithms, such as complex parameter adjustment and the need for later interactive verification and correction of filtering results.

[0008] The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints includes the following steps:

[0009] Step 1: Input the airborne LiDAR point cloud. According to the law of large numbers, first calculate the average distance from each point cloud to its K nearest neighbor points. Then, based on the principle that a point cloud greater than M times the standard deviation is a gross error, remove the point clouds that do not meet the conditions to obtain the point cloud to be processed.

[0010] Step 2: Build a grid index for the area where the point cloud to be processed is located, and interpolate the ground seed points obtained by the local minimum method using the thin plate spline function (TPS), remove the points whose height difference with the interpolation surface is greater than the curvature threshold, and obtain the ground points to be classified;

[0011] Step three, calculate the relative height difference ΔH of each ground point to be classified, apply the binary Gaussian mixture model (GMM) to classify the ground objects and the ground, use the expectation maximization algorithm (EM) to obtain the optimal solution of the model parameters, and judge whether the point belongs to the ground object based on the calculated posterior probability to complete the point cloud filtering.

[0012] The specific process of step one is:

[0013] (1) Input the airborne LiDAR point cloud, and according to the distribution of the point cloud in the study area, the user specifies the outlier denoising algorithm (SOR) parameters, which include the number of neighboring points K and the standard deviation multiplication factor M;

[0014] (2) According to the law of large numbers, we traverse the average distance of each point cloud to its K nearest neighbors, and then eliminate those points that do not meet the standard deviation based on the principle that those points that are greater than M times the standard deviation are gross errors. Point cloud to be processed is obtained; among them, is the average distance from the i-th point cloud to its K nearest neighbors, the KNN average distance Standard Deviation n is the number of input airborne LiDAR point clouds.

[0015] The specific process of step 2 is as follows:

[0016] (1) First, the initial ground seed point set G is obtained by the local minimum method for the point cloud to be processed in step 1. seed Given a grid size gw, construct a grid index, calculate the true plane coordinate value (x, y) of the center of the grid cell according to the grid row and column number (m, n), and then minimize the thin plate spline energy function E f Based on the optimization principle, N neighboring ground seed points of the grid point to be interpolated are selected as control points;

[0017] (2) The control point weights are set to satisfy the constraint that the interpolation function f(x, y) has a second-order partial derivative. The unknown coefficients in the interpolation function f(x, y) are solved according to the Lagrange multiplier method, and then the center elevation of the grid unit is interpolated.

[0018] (3) Given a curvature threshold parameter t, a surface grid pyramid R1 is constructed from gw and t, and a 3×3 mean moving window is used to pass through each grid where the point cloud to be processed is located to obtain a smoothed surface grid R new , according to R new Elevation value z new (m,n) and then calculate the average curvature c(m,n) = z new (m,n)+t; traverse the point cloud to be processed, if the jth point cloud to be processed P j (x j ,y j ,z j )’s elevation value z j If it is greater than or equal to c(m,n), it is added to the non-ground point, otherwise it is taken as the first ground point to be classified and enters (4);

[0019] (4) Let gw = gw / 2 and t = t+0.1, construct the surface grid pyramid R2, and traverse the first ground points to be classified to select the second ground points to be classified according to the method in (3), and enter (5);

[0020] (5) Let gw = gw / 2 and t = t+0.1, construct the surface grid pyramid R3, and traverse the second ground point to be classified according to the method in (3) to select the third ground point to be classified. The first to third ground points to be classified are collectively referred to as ground points to be classified.

[0021] In (1), the real plane coordinate values ​​(x, y) are converted according to the following formula:

[0022]

[0023] In the formula, (x min ,y min ) is the minimum coordinate value in the x and y directions of the minimum bounding box of the input airborne LiDAR point cloud.

[0024] Among them, the interpolation function in (2) is:

[0025]

[0026] Where a1, a2, and a3 are coefficients that reflect the overall trend of the interpolation surface; w s Represents the sth control point c s The relevant weight coefficient; is the radial basis kernel function associated with the thin plate spline function (TPS), (x s ,y s ) is c s The x and y coordinate values ​​of , N is the number of control points;

[0027] The second-order partial derivative constraints of f(x,y) are:

[0028]

[0029] In the formula, z s c s Coordinate value in the z direction;

[0030] According to the Lagrange multiplier method, the following linear equations are used to solve the unknown coefficients in f(x,y):

[0031]

[0032] In the formula, the element in the Ith row and Jth column of the matrix U is I,J∈[1,2,…,N],(x I ,y I ) and (x J ,y J ) are the first control point c I and the Jth control point c J The x and y coordinate values ​​of

[0033] Finally, the calculation formula for the elevation value z of the unit to be interpolated is: z=f(x,y).

[0034] The specific process of step three is:

[0035] (1) Traverse all the points to be classified, calculate the relative elevation of each point to be classified, substitute ΔH into the binary Gaussian mixture model (GMM), calculate the response of the Gaussian mixture model to the relative elevation of each point to be classified, and apply the expectation maximization algorithm (EM) to iteratively update and calculate the parameter value α of each sub-model l ,μ l ,σ l 2 ,l=1,2; among them, the kth point to be classified q k The relative height ΔH k =H–Hgrid, H is the elevation value of the point to be classified, Hgrid is the elevation of the surface grid unit of R3 where the point to be classified is located, α l is the mixing coefficient of the model and is greater than or equal to 0, μ l Yes k The average value of the relative elevation; σ l Yes k The standard deviation of the relative elevation;

[0036] (2) Since it is assumed that the relative elevations of ground points and non-ground points conform to Gaussian distribution, the conditional probability estimates containing latent variables are obtained based on the calculated Gaussian parameter values ​​of each sub-model, and then the posterior probability P(G|q k ), if the posterior probability P(G|q k ) is greater than or equal to 0.5, then the kth point to be classified q k belongs to the ground point and is used to construct the digital elevation model DEM; if the posterior probability P(G|q k ) is less than 0.5, then q k It is a non-ground point and is discarded;

[0037] Among them, the parameter value α of each sub-model is updated and calculated according to the following formula: l ,μ l ,σ l 2 ,l=1,2:

[0038]

[0039]

[0040]

[0041] In the formula, ΔH k is the ground point q to be classified k Relative height; P(G|q k ) is the point to be classified q k The posterior probability of belonging to a ground point; P(NG|q k) is the point to be classified q k The posterior probability of belonging to a non-ground point; N c is the total number of point clouds to be processed.

[0042] In (2), the point to be classified q is calculated according to the following formula: k The posterior probability P(G|q k ) and the posterior probability P(NG|q k ):

[0043]

[0044]

[0045] In the formula, P(G) is the point to be classified q k The prior probability of belonging to the ground point; P(NG) is the point to be classified q k The prior probability of belonging to a non-ground point; P(q k ) is calculated by the total probability formula, P(q k )=P(q k |G)P(G)+P(q k |NG)P(NG).

[0046] Each point to be classified q k The conditional probability of:

[0047]

[0048] In the formula, α l is the mixing coefficient of the model and is greater than or equal to 0; ΔH k is the pending point q k The relative elevation of; φ() is the Gaussian density function, θ l =(μ l , σ l ) is its parameter, μ l is the point to be classified q k The average value of relative elevation; σ l is the point to be classified q k The standard deviation of the relative elevation.

[0049] Apply the above formula to calculate the conditional probability P(q k |G), only one set of parameters (α1, μ1, σ1) is used:

[0050] P(q k |G)=α1φ(ΔH k |μ1,σ1)

[0051] Similarly, calculate the conditional probability P(q k|NG) uses another set of parameters (α2, μ2, σ2):

[0052] P(q k |NG)=α2φ(ΔH k |μ2,σ2)

[0053] The Gaussian density function is:

[0054]

[0055] The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints of the present invention is an adaptive filtering method, which combines the advantages of statistical methods and surface interpolation methods. The algorithm provides a relatively general set of parameters, which can be convenient for users who need filtering work. The filtering method of the present invention adopts a surface interpolation strategy, constructs a surface pyramid through a thin plate spline function (TPS), avoids the misclassification of terrain points caused by changes in local terrain slope, and ensures that the ground object points can be effectively distinguished and eliminated through curvature constraints; in addition, the relative height difference is calculated for each measurement point, which can effectively overcome the misclassification problem caused by the absolute elevation of the measurement point on the steep terrain, such as the cliff, being higher than the absolute elevation of the tree measurement point on the flat terrain, and then uses a binary Gaussian mixture model (GMM) to complete the filtering by probability discrimination. The solution of the model is automatically optimized by the expectation maximization algorithm (EM), avoiding excessive subjective judgment, and reducing a lot of manual parameter adjustment time. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 is a method framework diagram of an embodiment of the present invention;

[0057] Figure 2 1 is a schematic diagram of removing statistical outliers according to an embodiment of the present invention, wherein (a) is a schematic diagram before removal, (b) is a schematic diagram of determining an outlier set, and (c) is a schematic diagram after removal;

[0058] Figure 3 Schematic diagram of a multi-level curvature constrained surface grid according to an embodiment of the present invention, wherein (a) is a digital surface model (DSM) constructed by using test sample point cloud data provided by the third committee of ISPRS, (b) is a first-level surface grid R1 in the surface grid pyramid constructed by the embodiment of the present invention, (c) is a second-level surface grid R2 in the surface grid pyramid constructed by the embodiment of the present invention, and (d) is a third-level surface grid R3 in the surface grid pyramid constructed by the embodiment of the present invention;

[0059] Figure 4: is a measured point cloud filtering result diagram of an embodiment of the present invention, wherein (a) is a digital surface model (DSM) constructed using the original LiDAR point cloud, (b) is a digital elevation model (DEM) obtained using the Gaussian mixture model point cloud filtering method based on multi-level curvature constraints of the present invention; (c) is a digital elevation model (DEM) obtained using the classical mathematical morphology filtering method. DETAILED DESCRIPTION

[0060] The present invention will be further described in detail below in conjunction with specific embodiments and with reference to the accompanying drawings.

[0061] According to the law of large numbers, the present invention first calculates the average distance from each measurement point to its K nearest neighbor points, and then removes various abnormal measurement values ​​in the point cloud based on the principle that the normal measurement value will not be greater than M times the standard deviation; secondly, the ground seed points obtained by the local minimum are interpolated using the thin plate spline function (TPS), and then the points whose height difference with the interpolation surface is greater than the curvature threshold are removed. This process is to iteratively approximate a reference surface by continuously reducing the interpolation grid units and gradually increasing the curvature threshold; finally, the relative height difference of each measurement value is calculated, and the binary Gaussian mixture model (GMM) is applied to classify the ground objects and the ground, and then the expectation maximization algorithm (EM) is used to obtain the optimal solution of the model parameters, and the point cloud filtering is completed based on the calculated posterior probability to determine whether the point belongs to the ground object, so as to effectively reduce the time of manual parameter adjustment and improve the accuracy and efficiency of airborne LiDAR point cloud filtering.

[0062] As attached Figure 1 As shown in the figure, the Gaussian mixture model point cloud filtering method based on multi-level curvature constraints of the present invention includes three parts: ① determining the points to be processed in the point cloud according to the law of large numbers; ② multi-level curvature constraints based on thin plate spline interpolation; ③ relative elevation calculation and Gaussian mixture model classification. The specific implementation steps are:

[0063] Step 1: Determine the points to be processed in the point cloud according to the law of large numbers.

[0064] There are two types of noise points in the airborne LiDAR point cloud data: high outliers and low outliers. The KNN neighborhood search in the present invention can detect and eliminate these two types of noise.

[0065] The schematic diagram of determining the points to be processed in the point cloud according to the law of large numbers is shown in the attached figure. Figure 2 As shown in (a) to (c) in the figure, the specific process is as follows:

[0066] (1) Input the airborne LiDAR point cloud, and according to the distribution of the point cloud in the study area, the user specifies the outlier denoising algorithm (SOR) parameters, which include the number of neighboring points K and the standard deviation multiplication factor M;

[0067] (2) According to the law of large numbers, we traverse the average distance of each point cloud to its K nearest neighbors, and then eliminate those points that do not meet the standard deviation based on the principle that those points that are greater than M times the standard deviation are gross errors. Point cloud to be processed is obtained; among them, is the average distance from the i-th point cloud to its K nearest neighbors, the KNN average distance Standard Deviation n is the number of input airborne LiDAR point clouds.

[0068] Assume that the point cloud q to be calculated i (x i ,y i , z i ) to its nearest K neighbor point p m (x m ,y m , z m ), (m=1,2,...,K) is d im , then the calculation formula is as follows:

[0069]

[0070] Then its KNN average distance The calculation of can be given by the following formula:

[0071]

[0072] In the formula, the value of K is generally specified by the user.

[0073] According to The KNN average distance d can be calculated mean And the standard deviation stddev can be given by the statistical calculation formula:

[0074]

[0075] Where n is the total number of input airborne LiDAR point clouds.

[0076] Therefore, let the standard deviation multiply the coefficient M = 3, then the point cloud to be processed q i The screening criteria are as follows:

[0077]

[0078] If the point satisfies the above condition, it will be added to the pending point set Q; if the point does not satisfy the above condition, it will be added to the outlier point set O and discarded.

[0079] Step 2: Multi-level curvature constraints based on thin plate spline interpolation.

[0080] The schematic diagram of multi-level curvature constraint based on thin plate spline interpolation is shown in the attached figure. Figure 3 As shown in (a) to (d) in the figure, the process is as follows:

[0081] (1) First, the initial ground seed point set G is obtained by the local minimum method for the point cloud to be processed in step 1. seed Given a grid size gw, a grid index is constructed, and the true plane coordinate value (x, y) of the center of the grid unit is calculated according to the grid row and column number (m, n), and then the energy function E of the thin plate spline is minimized. f Based on the optimization principle, N neighboring ground seed points of the grid point to be interpolated are selected as control points;

[0082] Assuming that the row and column numbers of the grid to be interpolated are (m, n), they are converted into corresponding real plane coordinate values ​​(x, y) according to the following formula:

[0083]

[0084] Where gw is the grid size; (x min ,y min ) is the minimum coordinate value in the x and y directions of the minimum bounding box of the original point cloud.

[0085] Given a point, there are N = 12 control points c in 3D space s (x s ,y s , z s ), the thin plate spline function (TPS) passes through all control points c through the constructor f(x, y). s And minimize the bending energy function E f To complete the interpolation. In two-dimensional space, E f is defined as the integral of the square of the second-order partial derivatives and is given by:

[0086]

[0087] Where f is the interpolation function f(x, y); (dx, dy) is the differential of f(x, y) in the x, y directions.

[0088] (2) The control point weights are set to satisfy the constraint that the interpolation function f(x, y) has a second-order partial derivative. The unknown coefficients in the interpolation function f(x, y) are solved according to the Lagrange multiplier method, and then the center elevation of the grid cell is interpolated.

[0089] To be interpolated unit g idx For example, the true plane coordinates (x, y) are derived from formula (6), and the elevation z is calculated as follows:

[0090] First, we need to calculate the unknown parameters in f(x, y), and the calculation formula is given by the following formula:

[0091]

[0092] Where a1, a2, a3 are coefficients that reflect the overall trend of the interpolation surface; w s Represents and control point c s The relevant weight coefficient; is the radial basis kernel function associated with the thin plate spline function (TPS), (x s ,y s ) is the control point c s The x and y coordinate values ​​of f(x, y); at the same time, in order to satisfy the second-order partial derivative constraint of f(x, y), the following conditions need to be met.

[0093]

[0094] Then, according to the Lagrange multiplier method, the unknown coefficients in f(x, y) are solved using the linear equations. The calculation formula is as follows:

[0095]

[0096] In the formula, U is I, J∈[1, 2, ..., N] defines the matrix, N is the number of control points, (x I ,y I ), (x J ,y J ) are the x and y coordinates of the control point; and w is the coordinates of w s A is a column vector composed of a1, a2, and a3; accordingly, P and v can be expressed as follows:

[0097]

[0098] The elements in the matrix come from a set of control points c s (x s ,y s , z s ), (s=1, 2, ..., N), where N is the number of control points.

[0099] Finally, the interpolation unit g idx The calculation formula for the elevation value z is:

[0100] z=f(x,y) (12)

[0101] (3) Given the curvature threshold parameter t, construct the surface grid pyramid R1 from gw and t, and use a 3×3 mean shift window to pass through each grid where the point cloud to be processed is located, obtaining a smoothed surface grid R new , according to R new elevation value z new (m, n), and then calculate the mean curvature c(m, n) = z new (m, n) + t; traverse the point cloud to be processed. If the point P j (x j , y j , z j )'s elevation value z j is greater than or equal to c(m, n), then add it to the non-ground points; otherwise, take it as the first ground point to be classified and enter (4).

[0102] Obtain the initial ground seed point set G by the local minimum method seed , then given the grid size gw and the curvature threshold t, calculate the interpolated surface grid R1 by the formula z = f(x, y). Then traverse the point q j to be processed, and calculate its grid cell row and column numbers (m, n) by the formula (6). Use image convolution operation and a 3×3 mean shift window to process this grid cell, and calculate to obtain a smoothed surface R new , according to R new elevation value z new (m, n), and calculate the mean curvature c(m, n) according to the following formula:

[0103] c(m, n) = z new (m, n) + t (13)

[0104] In the formula, Z new (m, n) represents the smoothed surface R new grid cell where the row and column numbers (m, n) are located; t is the curvature threshold parameter.

[0105] Judge whether the elevation value z j (x j , y j , z j ) of the point q j to be processed satisfies the condition z j < c(m, n). If so, add it to the ground point set G to be classified; otherwise, add it to the non-ground point set NG.

[0106] (4) Let gw = gw / 2 and t = t + 0.1, construct the surface grid pyramid R2, and screen out the second ground point to be classified according to the method in (3) for the first ground point to be classified, and enter (5);

[0107] (5) Let gw = gw / 2 and t = t+0.1, construct the surface grid pyramid R3, and traverse the second ground point to be classified according to the method in (3) to select the third ground point to be classified. The first to third ground points to be classified are collectively referred to as ground points to be classified.

[0108] The specific steps for constructing surface grid pyramids R1, R2, and R3 are as follows:

[0109] First, the ground point set G to be classified is calculated using the initial gw and t by the process in (3); the updated ground seed point set is calculated by the local minimum method (the purpose of this step is to prevent the linear equation system from calculating ill-conditioned problems), and a new surface grid is established based on the updated point set. According to the curvature constraint condition of formula (13), the new ground point set G to be classified is obtained. new The process is repeated three times to complete the construction of the surface grid R1 at this level; then, let gw = gw / 2 and t = t+0.1 to construct the surface grid R2 at the next level, and use the local minimum method to calculate the updated ground seed point set G for the ground point set to be classified obtained in the previous step. new , a new surface grid is established based on the updated point set, and a new ground point set to be classified is obtained according to the curvature constraint of formula (13). This process is repeated three times to complete the construction of the surface grid R2 at this level; finally, let gw = gw / 2 and t = t+0.1 to construct the next level of surface grid R3, and the local minimum method is used again to calculate the updated ground seed point set G for the ground point set to be classified obtained in the previous step. new , a new surface grid is established based on the updated point set, and a new ground point set to be classified is obtained according to the curvature constraint condition of formula (13). This process is repeated three times to complete the construction of the surface grid R3 at this level.

[0110] Step 3: Relative elevation calculation and Gaussian mixture model classification.

[0111] (1) Traverse all the points to be classified, calculate the relative elevation of each point to be classified, substitute ΔH into the binary Gaussian mixture model (GMM), calculate the response of the Gaussian mixture model to the relative elevation of each point to be classified, and apply the expectation maximization algorithm (EM) to iteratively update and calculate the parameter value α of each sub-model l , μ l , σ l 2 , l = 1, 2; among them, the kth point to be classified q k The relative height ΔH k =H-Hgrid, H is the elevation value of the point to be classified, Hgrid is the elevation of the surface grid unit R3 where the point to be classified is located;

[0112] The calculation of relative elevation is to distinguish ground points from non-ground points more significantly. Because the general algorithm assumes that if the height difference of a point is greater than the height difference of its neighboring points, it will be classified as a non-ground point. However, this assumption will fail when encountering complex terrain such as steep cliffs. By calculating each point to be processed q k The relative height ΔH k , can effectively solve this problem, and its formula is as follows:

[0113] ΔH k =H-Hgrid (14)

[0114] Where H is the point to be processed q k The coordinate value in the z direction, Hgrid is the point to be processed q k The corresponding surface elevation value of R3(m,n).

[0115] ΔH k Substituting into the binary Gaussian mixture model (GMM), each point to be processed q can be expressed as k The conditional probability of:

[0116]

[0117] In the formula, α l is the mixing coefficient of the model and is greater than or equal to 0; ΔH k is the pending point q k The relative elevation of; φ() is the Gaussian density function, θ l =(μ l , σ l ) is its parameter. l is the pending point q k The average value of relative elevation; σ l is the pending point q k The standard deviation of the relative elevation; where the Gaussian density function is:

[0118]

[0119] Apply formula (15) to calculate the conditional probability P(q k |G), only one set of parameters (α1, μ1, σ1) is used:

[0120] P(q k |G)=α1φ(ΔH k |μ1,σ1) (17)

[0121] Apply formula (15) to calculate the conditional probability P(q k |NG) uses another set of parameters (α2, μ2, σ2):

[0122] P(qk |NG)=α2φ(ΔH k |μ2,σ2) (18)

[0123] Since there are both observed variables ΔH in the Gaussian mixture model k , there are also hidden variables G and NG, so the conventional maximum likelihood method cannot be used directly to estimate the parameters. At this time, the expectation maximization algorithm (EM) algorithm provides a solution to the approximate solution of the hidden variables. Usually, the expectation maximization algorithm contains four main steps:

[0124] i. Initialize input parameters, including α l , μ l , σ l , l = 1, 2.

[0125] First, we obtain the number of ground seed points and the total number of points to be processed N according to the multi-level curvature constraints. c , initialize α1 and α2. It should be noted that α1+α2=1. And μ1, μ2, σ1, σ2 are calculated by the relative elevation of the ground seed point and the non-ground seed point respectively.

[0126]

[0127] In the formula, ΔH k The meaning of is consistent with formula (15); N c is the total number of point clouds to be processed.

[0128] ii. Step E: Calculate the membership of the Gaussian mixture model for each point can be expressed as P(G|q k ) and P(NG|q k ).

[0129] iii. Step M: Calculate the parameters (α l , μ l , σ l 2 , l=1,2).

[0130] The M step in the iteration is to find the parameter value that maximizes the expected membership in the E step, α l , μ l , σ l 2 The calculation formula is as follows:

[0131]

[0132] In the formula, ΔH k The meaning of is consistent with formula (15).

[0133] iv. Determine the convergence conditions of the expectation-maximization algorithm.

[0134] The expectation-maximization algorithm (EM) repeats steps (i)-(iv) until the model parameters (α l , μ l , σ l 2 , l=1, 2) until it no longer changes.

[0135] (2) Since it is assumed that the relative elevations of ground points and non-ground points conform to Gaussian distribution, the conditional probability estimates containing latent variables are obtained based on the calculated Gaussian parameter values ​​of each sub-model, and then the posterior probability P(G|q k ), if the posterior probability P(G|q k ) is greater than or equal to 0.5, then the kth point to be classified q k belongs to the ground point and is used to construct the digital elevation model DEM; if the posterior probability P(G|q k ) is less than 0.5, then q k It is a non-ground point and is discarded;

[0136] After the model parameters are calculated by the expectation maximization algorithm (EM), the Gaussian mixture model formula (15-18) is used to directly calculate the conditional probability estimate P(q k |G), and the conditional probability estimate P(q k |NG); then, according to the Bayesian formula, the point to be processed q k The posterior probability P(G|q k ) and the posterior probability P(NG|q k ) is calculated as follows:

[0137]

[0138] Where P(G) represents the processing point q k The prior probability of belonging to the ground point; P(NG) represents the processing point q k The prior probability of belonging to a non-ground point; P(q k ) is calculated by the total probability formula, P(q k )=P(q k |G)P(G)+P(q k |NG)P(NG).

[0139] According to the condition P(G|q k )≥0.5, if it satisfies, the point to be classified q k is marked as a ground point; otherwise, it is classified as a non-ground point and discarded.

[0140] Therefore, the Gaussian mixture model point cloud filtering method based on multi-level curvature constraints of the present invention is a classification method based on probability theory and mathematical statistics, which completes the filtering of airborne LiDAR point clouds under the guidance of posterior probability. Figure 4 (a) to (c) in the figure show the filtering results of a measured point cloud sample of an ISPRS dataset, where (a) is a digital surface model (DSM) constructed using the original LiDAR point cloud, (b) is a digital elevation model (DEM) obtained using the Gaussian mixture model point cloud filtering method based on multi-level curvature constraints of the present invention; (c) is a digital elevation model (DEM) obtained using the classical mathematical morphology filtering method.

[0141] From the results, it can be seen that the digital elevation model (DEM) obtained by the classical mathematical morphology filtering method has poor filtering effect in areas with drastic terrain slope changes and mixed with low vegetation, resulting in many omission errors, and some ground point measurements are misclassified as non-ground points, resulting in erroneous DEM results. The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints in the present invention not only provides a set of more general algorithm parameter values, but also has a significant improvement in filtering accuracy.

[0142] The above description is only a specific implementation of the present invention, but the protection scope of the present invention is not limited thereto. Any person familiar with the technology can understand and think of any changes or substitutions within the technical scope disclosed by the present invention, which should be included in the scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A Gaussian mixture model point cloud filtering method based on multi-level curvature constraints, characterized in that: The following steps are involved: Step 1: Input the airborne LiDAR point cloud. According to the law of large numbers, first calculate the average distance from each point cloud to its K nearest neighbor points. Then, based on the principle that a point cloud greater than M times the standard deviation is a gross error, remove the point clouds that do not meet the conditions to obtain the point cloud to be processed. Step 2: A grid index is established for the area where the point cloud to be processed is located, and the ground seed points obtained by the local minimum method are interpolated using the thin plate spline function TPS, and the points whose height difference with the interpolation surface is greater than the curvature threshold are eliminated to obtain the ground points to be classified; Step 3: Calculate the relative height difference of each ground point to be classified, apply the binary Gaussian mixture model GMM to classify the ground object and the ground, use the expectation maximization algorithm EM to obtain the optimal solution of the model parameters, and judge whether the point belongs to the ground object based on the calculated posterior probability to complete the point cloud filtering; The specific process of step three is: (1) Traverse all the points to be classified, calculate the relative elevation of each point to be classified and substitute it into the binary Gaussian mixture model GMM, calculate the response of the Gaussian mixture model to the relative elevation of each point to be classified, and apply the expectation maximization algorithm EM to iteratively update and calculate the parameter value α of each sub-model l ,μ l ,σ l 2 ,l=1,2; among them, the kth point to be classified q k The relative height ΔH k =H–Hgrid, H is the elevation value of the point to be classified, Hgrid is the elevation of the surface grid unit of R3 where the point to be classified is located, α l is the mixing coefficient of the model and is greater than or equal to 0, μ l Yes k The average value of the relative elevation; σ l Yes k The standard deviation of the relative elevation; (2) Based on the calculated Gaussian parameter values ​​of each sub-model, the conditional probability estimate containing the hidden variables is obtained, and then the Bayesian formula and the total probability formula are used to solve q k The posterior probability P(G|q k ), if P(G|q k ) is greater than or equal to 0.5, then q k belongs to the ground point and is used to construct the digital elevation model DEM; if P(G|q k ) is less than 0.5, then q k It is a non-ground point and is discarded; in, P(q k )=P(q k |G)P(G)+P(q k |NG)P(NG),P(G) is q k The prior probability of belonging to a ground point, P(NG) is q k The prior probability of belonging to a non-ground point, P(q k |G)=α1φ(ΔH k |μ1,σ1),P(q k |NG)=α2φ(ΔH k |μ2,σ2), (1) The parameter value α of each sub-model is updated and calculated according to the following formula: l ,μ l ,σ l 2 ,l=1,2: In the formula, P(NG|q k ) is q k The posterior probability of belonging to a non-ground point; N c is the total number of point clouds to be processed; 2. The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints according to claim 1 is characterized in that: The specific process of step one is: (1) Input the airborne LiDAR point cloud, and according to the distribution of the point cloud in the study area, the user specifies the outlier denoising algorithm SOR parameters, which include the number of neighboring points K and the standard deviation multiplication factor M; (2) According to the law of large numbers, we traverse the average distance of each point cloud to its K nearest neighbors, and then eliminate those points that do not meet the standard deviation based on the principle that those points that are greater than M times the standard deviation are gross errors. Point cloud to be processed is obtained; among them, is the average distance from the i-th point cloud to its K nearest neighbors, the KNN average distance Standard Deviation n is the number of input airborne LiDAR point clouds.

3. The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints according to claim 1 is characterized in that: The specific process of step 2 is as follows: (1) First, the initial ground seed point set G is obtained by the local minimum method for the point cloud to be processed in step 1. seed Given a grid size gw, construct a grid index, calculate the true plane coordinate value (x, y) of the center of the grid cell according to the grid row and column number (m, n), and then minimize the thin plate spline energy function E f Based on the optimization principle, N neighboring ground seed points of the grid point to be interpolated are selected as control points; (2) The control point weights are set to satisfy the constraint that the interpolation function f(x, y) has a second-order partial derivative. The unknown coefficients in the interpolation function f(x, y) are solved according to the Lagrange multiplier method, and then the center elevation of the grid unit is interpolated. (3) Given a curvature threshold parameter t, a surface grid pyramid R1 is constructed from gw and t, and a 3×3 mean moving window is used to pass through each grid where the point cloud to be processed is located to obtain a smoothed surface grid R new , according to R new Elevation value z new (m,n) and then calculate the average curvature c(m,n) = z new (m,n)+t; traverse the point cloud to be processed, if the jth point cloud to be processed P j (x j ,y j ,z j )’s elevation value z j If it is greater than or equal to c(m,n), it is added to the non-ground point, otherwise it is taken as the first ground point to be classified and enters (4); (4) Let gw = gw / 2 and t = t+0.1, construct the surface grid pyramid R2, and traverse the first ground points to be classified to select the second ground points to be classified according to the method in (3), and enter (5); (5) Let gw = gw / 2 and t = t+0.1, construct the surface grid pyramid R3, and traverse the second ground point to be classified according to the method in (3) to select the third ground point to be classified. The first to third ground points to be classified are collectively referred to as ground points to be classified.

4. The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints according to claim 3 is characterized in that: (1) The real plane coordinates (x, y) are converted according to the following formula: In the formula, (x min ,y min ) is the minimum coordinate value in the x and y directions of the minimum bounding box of the input airborne LiDAR point cloud.

5. The Gaussian mixture model point cloud filtering method based on multi-level curvature constraints according to claim 3 is characterized in that: (2) Interpolation function: Where a1, a2, and a3 are coefficients that reflect the overall trend of the interpolation surface; w s Represents the sth control point c s The relevant weight coefficient; is the radial basis kernel function associated with the thin plate spline function TPS, (x s ,y s ) is c s The x and y coordinate values ​​of , N is the number of control points; The second-order partial derivative constraints of f(x,y) are: In the formula, z s c s Coordinate value in the z direction; According to the Lagrange multiplier method, the following linear equations are used to solve the unknown coefficients in f(x,y): In the formula, the element in the Ith row and Jth column of the matrix U is (x I ,y I ) and (x J ,y J ) are the first control point c I and the Jth control point c J The x and y coordinate values ​​of Finally, the calculation formula for the elevation value z of the unit to be interpolated is: z=f(x,y).

Images