CN113484844A

CN113484844A - Point cloud classification method driven by LiDAR full-waveform control decomposition

Info

Publication number: CN113484844A
Application number: CN202110976992.1A
Authority: CN
Inventors: 扆亮海
Original assignee: Individual
Current assignee: Individual
Priority date: 2021-08-24
Filing date: 2021-08-24
Publication date: 2021-10-08

Abstract

According to the point cloud classification method driven by LiDAR full-waveform control decomposition, wavelet factors are estimated from an original waveform one by adopting a layer-by-layer extraction method, estimated waveform characteristic factors are optimized by adopting a dynamic global algorithm, characteristic factors obtained by waveform decomposition are extracted and high-density point cloud is generated, and finally point cloud data are classified by adopting an SVM (support vector machine) based on the waveform characteristic factors; first, LiDAR full waveform data preprocessing: an algorithm suitable for removing LiDAR full-waveform random noise is provided, and an addition and subtraction alternative iterative correction denoising method is improved, so that the processing result effect is better; second, LiDAR full waveform data decomposition: introducing a dynamic global algorithm in the characteristic factor optimization, and adding control in the optimization process of the characteristic factor optimization to ensure that the optimization result does not deviate from the reality; and thirdly, point cloud classification based on waveform characteristic factors. The airborne small-light-spot full-waveform LiDAR rapid and accurate waveform decomposition and extraction are realized, and the precise and efficient classification of LiDAR point cloud is realized.

Description

Point cloud classification method driven by LiDAR full-waveform control decomposition

Technical Field

The application relates to a LiDAR full-waveform point cloud classification method, in particular to a point cloud classification method driven by LiDAR full-waveform control decomposition, and belongs to the technical field of LiDAR point cloud classification.

Background

The airborne LiDAR is a remote sensing data acquisition technology which is rapidly developed in recent years, can quickly and directly acquire a space three-dimensional coordinate, and is widely applied to the fields of topographic survey, urban modeling, power line extraction and the like due to the characteristics of high precision, high efficiency, low cost and the like.

The traditional discrete airborne LiDAR system only records the first echo and the last echo reflected by a ground object, the multi-echo system records 4 to 5 echoes at most, the two systems only can provide three-dimensional point coordinates and related intensity information for a user, the user can perform applications such as ground object classification, DEM generation and the like by using the three-dimensional coordinates and the intensity information through a filtering algorithm, the aspects such as ground object classification, urban three-dimensional reconstruction, forest evaluation and the like are performed according to geometric information, only the geometric information of point clouds is used, or the intensity information is added at most, and the precision of the ground object classification is greatly limited.

Commercial onboard small spot full waveform LiDAR systems are capable of recording the entire backscatter echo waveform of a scatterer at very small sampling intervals, and the user can extract more information from the full waveform data. The full waveform data can reflect the vertical distribution of the ground features, the geometric attributes and the physical attributes of the ground features are revealed, and the information extracted from the full waveform data is an important feature for classifying the ground features. The processing of the full waveform data can not only extract parameters (echo times, pulse width, amplitude and the like) reflecting the vertical structure and physical characteristics of ground objects, but also generate three-dimensional point cloud with higher density and higher precision.

Since the advent of airborne small-spot full-waveform LiDAR systems, the enormous amount of information it contains has attracted a great deal of research, development, and application. The current processing method of waveform data mainly comprises a threshold value method, a waveform decomposition method and a deconvolution method. The threshold method considers the signal with amplitude larger than a given value in the echo as a valid echo, and then calculates the position and amplitude of the target by adopting a local maximum value method. The threshold method is simple to implement, can quickly process waveforms, but only can effectively process simple waveforms, has a great problem on the accuracy of the processing result of complex echoes, and has no essential difference from the traditional method for acquiring the discrete point cloud by the LiDAR system in principle and does not give full waveform data. The waveform decomposition method considers the LiDAR full waveform data as superposition of a plurality of Gaussian waves, then decomposes each sub-waveform, extracts parameters (wave crest size, wave crest position and half wave width) of the sub-waveform, and the wave crest position of the sub-waveform is the reflection point position of multiple reflection.

Gaussian decomposition is firstly used for processing satellite-borne large-spot LiDAR full waveform data, the earliest relatively complete waveform data waveform decomposition algorithm is used for processing waveform data of a satellite-borne large-spot laser altimeter LVIS, and then an algorithm for extracting echo basic information of the laser altimeter is developed on the basis of the Gaussian decomposition by the Yang Heng and the like, so that the number of sub-waveforms of the waveform data and parameters such as the peak position, half wave width, amplitude and the like of each sub-waveform are obtained. Current methods of waveform decomposition for airborne LiDAR full waveform data include: a hierarchical gaussian function waveform fitting method based on nonlinear least squares, a gauss-newton method, an Expectation Maximization (EM) algorithm, a waveform decomposition method of an improved EM algorithm, and the like.

The prior art of ground feature classification using airborne LiDAR waveform data includes: the method comprises the steps of carrying out deconvolution processing on full waveform data, dividing the data into ground points and non-ground points by a method of reserving the last sub-waveform, classifying the data by utilizing the waveform data, and further carrying out waveform decomposition based on the data, obtaining the peak position of the sub-waveform by utilizing the waveform decomposition, calculating high-density three-dimensional point cloud data, and obtaining information reflecting the reflection characteristics of ground objects such as the peak size, the half wave width and the like. The initial research was limited to the classification of terrain using high density point clouds computed using waveform decomposition parameters.

In the prior art, unsupervised classification is carried out on vegetation through decomposition of full waveform data of LiDAR, but parameters obtained through Gaussian decomposition are not used in the classification process, and the structure of the vegetation is analyzed by only using dense point cloud calculated by waveform parameters, so that the purpose of classification is achieved. Then, the feature factors extracted by waveform decomposition are gradually added in the research of ground feature classification by using waveform data, and the method is mostly used for classification research of vegetation areas. On the basis of waveform decomposition, the trees in the forest are classified by utilizing the wave width and the intensity value through a K-means clustering method and an EM algorithm respectively, and the result proves that the effect is better than that of classifying by only using the first echo/the last echo. The prior art uses waveform decomposition parameters to classify different types of regions, including unsupervised classification, supervised classification, and classification based on building decision trees. The supervised classification and the unsupervised classification are similar to the remote sensing image classification, the adopted algorithms comprise a maximum likelihood method, an EM algorithm, a support vector machine and the like, most of research applications perform feature selection while classifying the ground features, and the influence of different waveform parameters on different ground feature classifications is analyzed. Decision tree classification is then based on statistics of waveform decomposition parameters. However, due to the limitation of various conditions, the LiDAR point cloud classification method in the prior art has various problems, and a method which can achieve higher level of point cloud classification accuracy and efficiency is lacked.

In summary, the LiDAR point cloud classification method in the prior art has disadvantages, and the difficulties and problems to be solved in the present application mainly focus on the following aspects:

firstly, the traditional discrete airborne LiDAR system only records the first echo and the last echo reflected by a ground object, the multi-echo system records 4 to 5 echoes at most, the two systems only can provide three-dimensional point coordinates and related intensity information for a user, and only utilizes the geometric information of point cloud or adds the intensity information at most aiming at the aspects of ground object classification, urban three-dimensional reconstruction, forest evaluation and the like of the geometric information, so that the precision of ground object classification is greatly limited; the processing of the full waveform data can not only extract parameters reflecting the vertical structure and physical characteristics of ground features, but also generate three-dimensional point cloud with higher density and higher precision, and compared with the discrete point cloud data acquired by the multi-echo airborne LiDAR, the waveform data of the full waveform airborne LiDAR can provide more target information, but also put higher requirements on the analysis and information extraction of the full waveform data, but the prior art cannot extract waveform decomposition parameters quickly and accurately, and cannot obtain accurate point cloud results in the subsequent processing;

secondly, for the removal of the background noise of LiDAR waveform data, in the prior art, ten sampling points before and after the waveform are mainly regarded as data without ground echo information, and the average value of the data is regarded as the background noise removal. When the number of the echoes is more in complex ground conditions, the part exceeding the sampling number of the system is not recorded, and the sampling point of the waveform digitization end contains ground echo information; when high ground objects exist, sampling points at which waveform digitization starts also contain echo data of the ground objects, if the average value of the points is simply taken as background noise to be removed from the original waveform, the removed background noise is larger than an actual value, and weak ground echoes are mistaken for noise, so that the waveform after denoising is inconsistent with the actual value, the final decomposition result is influenced, and further the calculation and classification result of point cloud coordinates after denoising is influenced, and the prior art lacks a method for removing full waveform background noise;

thirdly, the noise distribution in LiDAR signals is very complex, signal denoising is unrealistic only by selecting a cut-off frequency through a traditional digital filter, the waveform of the denoised signals is similar to that of a wavelet function through discrete wavelet transform, distortion is serious, the wavelet multi-scale decomposition reconstruction calculation amount is large, and if the LiDAR is used for denoising echoes for thousands of times, the consumed time space is too large; for the denoising of the space domain, if a better denoising effect is needed, prior knowledge of LiDAR signals needs to be added when a support vector machine is trained, the LiDAR waveforms are not superimposed by Gaussian distribution, the denoising result is not ideal, the traditional image space domain processing filtering algorithm does not need the prior knowledge of the waveforms, but the waveforms are distorted while denoising, and the prior art lacks a method for removing LiDAR full-waveform random noise in a targeted manner;

fourthly, after the number of sub-waveforms and the initial values of the waveform characteristic factors are estimated, the estimated initial values of the characteristic factors need to be optimized, Gaussian waves fitted by the sub-waveforms are enabled to be closest to original waveform data, the waveform decomposition problem of airborne LiDAR full waveform data belongs to a multi-dimensional nonlinear optimization problem, algorithms commonly used for solving the problems comprise a gradient method, a quasi-Newton method, an LM algorithm and an EM algorithm, the LM algorithm overcomes the problem that search fails when a Gaussian-Newton algorithm Jacobian matrix is a non-column matrix, and the LM algorithm is easy to obtain a local optimal solution; the evolution method adopts multi-point parallel search, generates new individuals through crossing and variation in each iteration process, and continuously enlarges the search range, so that the evolution algorithm is easy to search out the global optimal solution rather than the local optimal solution; the prior art lacks a method to optimize the characteristic factors for LiDAR full waveform decomposition.

Disclosure of Invention

In order to solve the problems, the method comprises the steps of estimating wavelet factors from original waveforms one by adopting a layer-by-layer extraction method based on Gaussian attributes of LiDAR full waveform data, optimizing the estimated waveform characteristic factors by adopting a dynamic global algorithm, extracting characteristic factors obtained by waveform decomposition and generating high-density point cloud, and finally classifying point cloud data by adopting an SVM based on the waveform characteristic factors, wherein the method comprises the following steps: first, LiDAR full waveform data preprocessing: based on the characteristics of LiDAR full waveform data and noise thereof, an algorithm suitable for removing LiDAR full waveform random noise is provided, and an addition and subtraction alternative iterative correction denoising method is improved, so that the processing result effect is better; second, LiDAR full waveform data decomposition: estimating initial characteristic factors of LiDAR waveform data by adopting a one-by-one extraction method, filtering the initial characteristic factors based on an error theory, introducing a dynamic global algorithm in characteristic factor optimization, and adding control in the optimization process so that the optimization result does not deviate from the reality; thirdly, point cloud classification based on waveform characteristic factors: on the basis of analyzing the ground feature reflection characteristics reflected by the waveform characteristic factors, an SVM classifier is adopted, the point cloud is classified by utilizing waveform parameters obtained by waveform decomposition and point cloud coordinates calculated by utilizing the waveform parameters, in the training process of the classifier, the punishment coefficient and the kernel function parameters are optimized by utilizing a dynamic global algorithm by taking classification precision as a fitness function, the rapid and accurate waveform decomposition and extraction of the airborne small-spot full-waveform LiDAR are realized, and the accurate and efficient classification of the LiDAR point cloud is realized.

In order to realize the technical characteristics, the technical scheme adopted by the application is as follows:

a point cloud classification method driven by LiDAR full-waveform control decomposition is characterized in that full-waveform data are subjected to waveform decomposition based on Gaussian attributes of airborne LiDAR full-waveform data, waveform characteristic factors are estimated and extracted one by one, a dynamic global algorithm is adopted to optimize the characteristic factors, LiDAR waveform characteristic factors are extracted and high-density point cloud is generated, and finally the point cloud is classified by utilizing waveform decomposition parameters and point cloud height through a support vector machine; the method mainly comprises the following steps:

first, LiDAR full waveform data preprocessing: based on the characteristics of noise in the LiDAR full waveform data, a method for removing background noise and random noise in the LiDAR full waveform data is provided, on the basis of verifying the denoising superiority of the LiDAR waveform data by the addition and subtraction alternative correction denoising method, the addition and subtraction alternative correction denoising method is improved, the improved addition and subtraction alternative iterative correction denoising method is provided, the random noise is removed, the shape of the waveform is kept, and the signal-to-noise ratio of the waveform data is improved;

second, LiDAR full waveform decomposition: the waveform decomposition process comprises pre-estimating sub-waveform characteristic factors and optimizing characteristic factors, and in the pre-estimating characteristic factor part, the characteristic factors of the sub-waveforms are continuously pre-estimated from original waveform data by adopting a one-by-one extraction method until the residual waveforms only contain noise points; the optimized waveform characteristic factor introduces a dynamic global algorithm and controls the change of a peak value in the optimization process, so that the maximum peak value can not emit violent change in the optimization process; finally, generating high-density point cloud data by using parameters obtained by waveform decomposition;

thirdly, point cloud classification based on waveform characteristic factors: on the basis of analyzing the ground feature reflection characteristics reflected by the waveform characteristic factors, an SVM classifier is adopted, the point cloud is classified by utilizing the waveform characteristic factors obtained by waveform decomposition and point cloud coordinates calculated by utilizing waveform parameters, in the classifier training process, classification accuracy is taken as a fitness function, a penalty coefficient and kernel function parameters are optimized by utilizing a dynamic global algorithm, wherein the kernel function selects a radial basis function, the dynamic global algorithm is adopted for optimizing the classification parameters in the SVM classifier training process, and finally the classifier trained by sample data is used for classifying the point cloud data of the whole measuring area.

A point cloud classification method driven by LiDAR full-waveform control decomposition, further an improved addition and subtraction alternative iterative correction denoising method: for a one-dimensional signal z, it is represented as a column vector (z)₁，…，z_m)^tThe simple form of filtering is expressed as equation 1:

z_i′＝z_i+aΔz_iformula 1

Wherein

a is a scale factor greater than 0 and less than 1, and the equation is written in matrix form, as in equation 2:

z ═ I-aW z formula 2

Wherein W matrix is as follows:

replacing the matrix I-aW with a non-shrinking function on the matrix W, the following result:

z ═ g (W) z formula 4

If iterate M times, the result is as follows:

z^N＝g(W)^N formula 5

The W matrix is symmetrical and has real eigenvalue and eigenvector, and the real eigenvalue of the W matrix is set to be more than or equal to 0 and less than or equal to W₁≤w₂≤…≤w_m

N is less than or equal to n, and the corresponding feature vectors are v1, …, v_mFormula 4 is represented by formula 6:

where g is the filter kernel, h_iIs a coefficient, satisfies that after M iterations, when w belongs to [0, n ∈]Time, low frequency component g (w)_i)^M1, high frequency component g (w)_i)^N0, the kernel function takes the form of equation 7:

g (w) ═ 1-aw) (1-cw) formula 7

Where c is a new negative scale factor and c < -a, which is equivalent to a further approximation step after gaussian smoothing with a positive scale factor by equation 1:

z_i′＝z_i+cΔz_iformula 8

Since g (0) ═ 1, a + c < 0, a threshold value w is sought_QLet g (w)_Q) 1, a, c satisfy formula 9:

when the Gaussian smoothing is carried out by using the formula 1, the subtraction operation of the formula 9 is added, and the two calculations are alternately carried out, so that the defect of the Gaussian filtering is overcome.

A point cloud classification method driven by LiDAR full-waveform control decomposition, further an improved alternative iterative correction denoising method:

the addition and subtraction alternate correction denoising method modifies a kernel function on the basis of Gaussian filtering, adopts an addition and subtraction alternate calculation method to inhibit contraction generated by the Gaussian filtering, but because c is less than-a, the addition and subtraction degrees are different, the denoising result can be influenced, the addition and subtraction alternate correction denoising method is further improved, multiple alternate operations are performed on the addition and subtraction alternate correction denoising, and the following operations are performed when the operand is an odd number:

z_i′＝z_i+aΔz_iwhen i is an odd number;

z_i′＝z_i+cΔz_iwhen i is an even number;

when the operand is an even number, the following operations are performed:

z_i′＝z_i+aΔz_iwhen i is an even number;

z_i′＝z_i+cΔz_iwhen i is an odd number;

each processing of the plus-minus alternative correction denoising method is weighted correction with a single direction, and correction directions of adjacent points are different, so that although the difference between the point and the neighborhood average value is small after multiple times of processing, the difference between the adjacent two points is likely to be large, and the noise elimination strength is insufficient; the improved alternative iterative correction denoising method adopts alternative addition and subtraction correction denoising along with the increase of the processing times, namely, in the processing of adjacent times, the modification of the same point is in different directions, so that the deformation introduced by the method can be counteracted by the iterative processing, and the denoising effect is further improved.

The LiDAR full-waveform control decomposition driving point cloud classification method comprises the following steps of: LiDAR full waveform data is formed by the superposition of several single Gaussian waves, represented by equation 10:

where k is the number of single Gaussian waves, g_j(z) is the probability density function of the jth Gaussian distribution, d_jIs the size of the wave crest, c_jIs the mean value, i.e. the peak position, σ_jThe standard deviation is half wave width, the number of the sub-waveforms, the size of the wave peak, the position of the wave peak and the optimal value of the half wave width characteristic factor are calculated by waveform decomposition, so that the superposition result of the sub-waveforms is closer to the original waveform;

LiDAR full waveform decomposition procedure: firstly, preprocessing waveform data, removing background noise and random noise in the waveform data, namely the content of the previous part, then estimating the number of sub-waveforms and characteristic factors of each sub-waveform, wherein the characteristic factors comprise the size of a wave crest, the position of the wave crest and half wave width, and finally optimizing the estimated initial value of the characteristic factor by using a dynamic global adjustment method to obtain an accurate sub-waveform characteristic factor.

The point cloud classification method driven by LiDAR full-waveform control decomposition is further characterized in that an initial characteristic factor is estimated: estimating the number of sub-waveforms and the characteristic factor of each sub-waveform;

the method for estimating the initial parameters of the LiDAR waveforms which are extracted one by one comprises the following specific processes: firstly, the original waveform data is preprocessed, background noise and random noise are removed to obtain new waveform data L_mDetecting L_mThe maximum value of the wave-shape is taken as the wave peak value of the sub-waveformD, taking the position of the maximum value as the peak position e of the sub-waveform, and in order to avoid the influence of random noise as much as possible, D needs to be more than 3 times of the error in the random noise, namely more than 3 sigma_mOtherwise, the waveform data L is considered_mIf the position r of the sampling point with half of the wave peak value D is known, the noise remains_gThe half-wave width σ of the sub-waveform can be found due to the discreteness of the waveform data sampling points, and r_gIs not necessarily exactly at the sampling point, so r cannot be directly obtained_gR is obtained by statistically analyzing the sampled data on both sides of the peak_gThen, calculating an approximate value of the sigma;

r_gthe analysis calculation obtaining method comprises the following steps: firstly judging the waveform data on the right side of the wave crest, if the value of the w-1 waveform data sampling point is greater than D/2 and the value of the w +1 waveform data sampling point is less than D/2, then r_gThe approximate value of the wave data is the position of the w-th wave data sampling point, then the approximate value of sigma is calculated, then the wave data on the left side of the wave crest is judged, if the value of the w-1 th wave data sampling point is less than D/2 and the value of the w +1 th wave data sampling point is more than D/2, then r_gThe approximate value of the sigma is calculated by taking the position of the w-th waveform data sampling point; when the waveform data on the two sides are judged to be finished, if a waveform data sampling point meeting the condition is not searched on one side, the r searched on the other side is selected_gTo calculate an approximation of σ as a final result; if waveform data sampling points on two sides of a wave crest can be searched to meet the condition r_gThen r searched from the left and right sides is compared_gAnd taking a smaller sigma value as a final sigma approximate value according to the calculated sigma approximate value:

after calculating the approximate value of sigma, recording (D, e, sigma) as the initial values of three characteristic factors of the estimated sub-waveform, subtracting the sub-waveform from the original waveform data, repeating the estimation steps for the rest waveform data, and continuously estimating the estimated sub-waveform from the original waveform dataSubtracting from the waveform data until the residual waveform does not satisfy L_max＞3σ_mAnd (4) indicating that only noise remains in the current waveform data, and finishing the estimation of the sub-waveform characteristic factor.

The point cloud classification method driven by LiDAR full-waveform control decomposition further comprises the following steps of dynamically and globally adjusting and optimizing characteristic factors: according to the method, the optimization process of the waveform characteristic factors is improved, based on the characteristic factor estimation method which is extracted one by one, the estimated first sub-waveform is a pulse with the largest wave peak value, namely the strongest echo signal, the strongest signal is less influenced by noise, the detection result is closer to the real situation, and the estimated sub-waveform is more easily influenced by the noise due to the weaker signal with the smaller wave peak value, so that the change of the size and the position of the wave peak of the first sub-waveform is controlled in the process of optimizing the characteristic factors by adopting a dynamic global algorithm, and the optimization result of the first sub-waveform cannot deviate from the actual situation due to too large change.

The point cloud classification method driven by LiDAR full-waveform control decomposition is characterized in that the specific process of dynamic global optimization of characteristic factors comprises the following steps:

firstly, initializing a characteristic factor: determining the size N of the population, and randomly selecting an initial population:

Z(r)＝(Z1(r)，Z2(r)，…Z_m(r)) formula 12

Wherein Z_i(r) represents the ith individual in the r-th generation, r is 0, Z_i(0) Is an m-dimensional vector;

and a second process, characteristic factor group evolution: for each individual Z in Z (r)_i(r) performing the following operations:

firstly, mutation operation: for randomly selecting two individuals z from characteristic factor group_p1，z_p2The following operations are carried out:

u_ij(r+1)＝z_best，j(r)+G(z_p1，j(r)-z_p2，j(r)) formula 13

Wherein u is_ij(r +1) is the jth component of the ith individual of the (r +1) th generation, z_best，j(r) is the best individual vector in the r generation, G isVariation factor, z_p1，j(r)-z_p2，j(r) is a difference vector, i, p1, p2 are different from each other and p1, p2 are [1, N]Any two random integers,

i

1, 2, and N,

j

1, 2, m;

and step two, cross operation: increasing the diversity of the characteristic factor group, and the calculation process is as follows:

wherein S is cross probability, and S is more than 0 and less than or equal to 1, randk_ijIs [0, 1 ]]Random (i) is [1, n ]]Random integer between, cross-operation guarantees v_ijAt least one component of (r +1) is composed of u_ij(r + 1);

selecting operation: to judge Z_iWhether it can become a member of the next generation or not, vector V is calculated_iAnd a target vector Z_iThe fitness of (c) is compared when V_iFitness ratio Z of_iHigh is selected as the offspring, otherwise Z is directed_iAs a child, the selection operation is calculated in the following manner:

wherein g (Z)_i(r)) is the fitness of the ith generation of the ith individual;

and a third process, namely termination judgment: let the new population resulting from process two be:

Z(r+1)＝(Z₁(r+1)，Z₂(r+1)，…，Z_M(r +1)) formula 16

The most optimal individual in Z (r +1) is designated as Z_best(r +1), when the result meets the precision requirement or the whole evolution process has reached the maximum algebra, the operation is terminated and Z is added_bestAnd (r +1) is output as an optimal solution, otherwise, r is made to be r +1, and the process is switched to a second process.

And (3) a point cloud classification method driven by LiDAR full-waveform control decomposition, and further generating high-density point cloud data: the characteristic factors of the sub-waveforms are extracted through waveform decomposition of LiDAR waveform data, the characteristic factors comprise peak positions, peak sizes and half wave widths, a point data recording format comprises three parameters X (r), Y (r) and Z (r), the three parameters calculate space three-dimensional coordinates of points on corresponding waveforms, and the space coordinates of the points on the waveforms are calculated by a formula 17:

x, Y, Z is the spatial location of the point calculated from the sub-waveform feature factors, X₀，Y₀，Z₀Is the position of the starting point, i.e. the three-dimensional coordinates of the point data, X (r), Y (r), Z (r) are the velocity vector of the laser signal, r is the time relative to the starting point, i.e. the peak position of the sub-waveform, in picograms seconds, the product of X (r), Y (r), Z (r), and r is the displacement of the laser pulse from the starting time to the time at which the peak is detected, the position of the starting point plus this displacement yields the coordinates of the point corresponding to the sub-waveform, and the units of X, Y, Z are the units of the LAS data coordinate system.

The point cloud classification method driven by LiDAR full-waveform control decomposition is further based on the point cloud classification of waveform characteristic factors: the method comprises the steps of classifying point cloud data based on a support vector machine, classifying sample data by using four kernel functions in a LibSVM toolbox respectively, selecting a radial basis function with high speed and high classification precision as a kernel function of an SVM, optimizing classification parameters by using a dynamic global algorithm of the application and using the classification precision as a fitness function, giving training samples and corresponding classes, optimizing classifier parameters by using the dynamic global algorithm of the application and using the classification precision as the fitness function, and classifying the data of the whole test area by using the classifier.

Compared with the prior art, the innovation points and advantages of the application are as follows:

first, conventional discrete airborne LiDAR typically records only the first and last echoes reflected by terrain, and multi-echo systems record up to 4 to 5 echoes, both of which provide the user with three-dimensional point coordinates and associated intensity information. The airborne small-spot full-waveform LiDAR system can record the whole backscatter echo waveform of a scatterer at a very small sampling interval, and a user can extract more information by performing autonomous processing analysis on full-waveform data. According to the method, based on Gaussian attributes of airborne LiDAR full waveform data, waveform decomposition is carried out on full waveform data, a one-by-one extraction method is adopted for estimating waveform characteristic factors, a dynamic global algorithm is adopted for optimizing the characteristic factors, LiDAR waveform characteristic factors are extracted and high-density point clouds are generated, and finally the point clouds are classified by utilizing waveform decomposition parameters and point cloud elevation through a support vector machine; through a series of innovative improvements, rapid and accurate waveform decomposition and extraction of airborne small-spot full-waveform LiDAR are realized, and accurate and efficient classification of LiDAR point cloud is realized;

secondly, the first step LiDAR full waveform data preprocessing of the application provides a method for removing background noise and random noise in LiDAR full waveform data based on the characteristics of noise in the LiDAR full waveform data, improves an addition and subtraction alternative correction denoising method on the basis of verifying the advantage of the addition and subtraction alternative correction denoising method in denoising the LiDAR waveform data, provides an improved addition and subtraction alternative iterative correction denoising method, performs experiments on real waveform data, compares the effects of a traditional image processing filtering algorithm, an addition and subtraction alternative correction denoising filtering method and the improved alternative iterative correction denoising method in removing random noise of the LiDAR waveform data, points out that the reduction of a peak and the increase of a half-wave width are caused when the LiDAR waveform is denoised by the traditional image processing algorithm, proves the advantage of the improved addition and subtraction alternative iterative correction denoising method, can effectively remove random noise and simultaneously keep the shape of the waveform, the signal-to-noise ratio of the waveform data is improved, and the subsequent processing of the waveform data is facilitated;

thirdly, the second step LiDAR full waveform decomposition of the application comprises pre-estimating sub-waveform characteristic factors and optimizing characteristic factors, and in the pre-estimating characteristic factor part, the application adopts a one-by-one extraction method to continuously estimate the characteristic factors of the sub-waveforms from the original waveform data until the residual waveforms only contain noise points; the optimized waveform characteristic factor introduces a dynamic global algorithm and controls the change of a peak value in the optimization process, so that the maximum peak value can not emit violent change in the optimization process; finally, high-density point cloud data are generated by using parameters obtained by waveform decomposition, and compared with point cloud provided by a system, the density and the quality are greatly improved;

fourthly, in the third step of the method, based on the point cloud classification of the waveform characteristic factors, on the basis of analyzing the ground object reflection characteristics reflected by the waveform characteristic factors, an SVM classifier is adopted, the point cloud is classified by utilizing the waveform characteristic factors obtained by waveform decomposition and point cloud coordinates calculated by utilizing waveform parameters, the classification precision is taken as a fitness function, a penalty coefficient and a kernel function parameter are optimized by utilizing a dynamic global algorithm, the classification precision of test sample data by the SVM classifier in the training process of the SVM classifier is up to 97.9%, finally the point cloud data of the whole measuring area is classified by utilizing the classifier trained by sample data, the classification result is analyzed, and the precise classification high efficiency of the LiDAR point cloud is obtained, and particularly, the classification precision of the ground, trees and buildings is greatly improved.

Drawings

FIG. 1 is a graph comparing the background noise removal effect of different methods of the present application.

Fig. 2 is a schematic diagram of denoising effects of the methods on real waveform data when a window takes 5.

FIG. 3 is a detail diagram of the addition and subtraction alternate correction denoising filter and the improved alternate iterative correction denoising effect.

FIG. 4 is a diagram showing PSNR comparison after denoising of real waveform data of each method.

FIG. 5 is a comparison graph of waveform characteristic parameters before and after de-noising of real waveform data of each method.

FIG. 6 is a flow chart of the present LiDAR full waveform decomposition.

FIG. 7 is a detailed flow diagram of a method for estimating initial parameters of a one-by-one extracted LiDAR waveform of the present application.

FIG. 8 is a comparison graph of estimation of waveform characteristic factors and the optimization results of LM algorithm.

FIG. 9 is a comparison graph of the point cloud provided by the system and calculated by the waveform decomposition factor.

FIG. 10 is a comparison graph of the density of the corresponding point clouds of FIG. 9 (unit: dot/square meter).

FIG. 11 is a flowchart of the present application for point cloud classification based on waveform feature factors.

Fig. 12 is a schematic diagram of a point cloud data classification result according to the whole embodiment of the present application.

FIG. 13 is a schematic diagram of the point cloud data of the embodiment of FIG. 12 separately displaying three types of ground features.

Detailed description of the invention

The technical solution of the point cloud classification method driven by LiDAR full-waveform control decomposition provided by the present application will be further described below with reference to the accompanying drawings, so that those skilled in the art can better understand the present application and can implement the same.

Conventional discrete airborne LiDAR typically records only the first and last echoes reflected by terrain, and multi-echo systems record up to 4 to 5 echoes, both of which provide the user with three-dimensional point coordinates and associated intensity information. The airborne small-spot full-waveform LiDAR system can record the whole backscatter echo waveform of a scatterer at a very small sampling interval, and a user can extract more information by performing autonomous processing analysis on full-waveform data.

According to the method, based on Gaussian attributes of airborne LiDAR full waveform data, waveform decomposition is carried out on full waveform data, a one-by-one extraction method is adopted for estimating waveform characteristic factors, a dynamic global algorithm is adopted for optimizing the characteristic factors, LiDAR waveform characteristic factors are extracted and high-density point clouds are generated, and finally the point clouds are classified by utilizing waveform decomposition parameters and point cloud elevation through a support vector machine; the method mainly comprises the following steps:

first, LiDAR full waveform data preprocessing: the method for removing background noise and random noise in the full waveform data of the LiDAR is improved on the basis of verifying the advantage of the denoising method for the LiDAR waveform data, the improved denoising method for the addition-subtraction alternating correction is provided, the experiment is carried out on the real waveform data, the effect of the traditional image processing filtering algorithm, the addition-subtraction alternating correction denoising filtering and the improved denoising method for the LiDAR waveform data random noise removal is compared, the reduction of the peak and the increase of the half-wave width can be caused when the LiDAR waveform is denoised by the traditional image processing algorithm, and the advantage of the improved denoising method for the addition-subtraction alternating iterative correction is proved, the shape of the waveform can be kept while random noise is effectively removed, the signal-to-noise ratio of waveform data is improved, and subsequent processing of the waveform data is facilitated;

second, LiDAR full waveform decomposition: the LiDAR full waveform number is formed by superposing a plurality of single Gaussian waves, the waveform decomposition obtains the characteristic factors of the single Gaussian waves, including the size of the peak, the position of the peak and the half wave width, the waveform decomposition process comprises pre-estimating sub-waveform characteristic factors and optimizing characteristic factors, in the pre-estimating characteristic factor part, the application adopts a one-by-one extraction method to continuously pre-estimate the characteristic factors of the sub-waveforms from the original waveform data until the rest waveforms only contain noise points; the optimized waveform characteristic factor introduces a dynamic global algorithm and controls the change of a peak value in the optimization process, so that the maximum peak value can not emit violent change in the optimization process; finally, generating high-density point cloud data by using parameters obtained by waveform decomposition, and comparing the high-density point cloud data with the density of the point cloud provided by the system;

thirdly, point cloud classification based on waveform characteristic factors: on the basis of analyzing the ground feature reflection characteristics reflected by the waveform characteristic factors, an SVM classifier is adopted, the point cloud is classified by utilizing the waveform characteristic factors (wave crest size, wave crest position and half wave width) obtained by waveform decomposition and point cloud coordinates calculated by utilizing waveform parameters, in the classifier training process, classification precision is taken as a fitness function, a penalty coefficient and kernel function parameters are optimized by utilizing a dynamic global algorithm, wherein the kernel function is a radial basis function, the dynamic global algorithm is adopted for optimizing the classification parameters in the SVM classifier training process, the classification precision of test sample data by the SVM classifier reaches 97.9%, finally, the classifier trained by sample data is used for classifying the point cloud data of the whole test area, and classification results are analyzed.

First, LiDAR full waveform data preprocessing

Removing background noise

Background noise in LiDAR waveform data is removed, data processing amount can be reduced, and operation speed is improved. For the removal of the background noise of LiDAR waveform data, in the prior art, front and back ten sampling points of a waveform are mainly regarded as data without ground echo information, and the average value of the data is regarded as background noise removal.

When the number of the echoes is more in complex ground conditions, the part exceeding the sampling number of the system is not recorded, and the sampling point of the waveform digitization end contains ground echo information; when high ground objects exist, sampling points at which waveform digitization starts also contain echo data of the ground objects, if the mean value of the points is simply taken as background noise to be removed from the original waveform, the removed background noise is larger than an actual value, and weak ground echoes are mistaken for noise, so that the denoised waveform is not consistent with the actual value, the final decomposition result is influenced, and the calculation and classification result of the point cloud coordinates is influenced. As shown in the waveform data in fig. 1(a), the first ten sampling points contain echo data, and the background noise estimated by using the first ten sampling points and the last ten sampling points is large, so that the situation that the echo intensity is a negative value can occur when the background noise is removed by using the method, the estimated characteristic factor after the influence is influenced, in order to avoid the phenomenon, the minimum value in the waveform data is taken as the background noise to be removed from the original waveform, and the result of removing the background noise by using the method is shown in fig. 1 (b).

(II) removing random noise

The LiDAR waveform data denoising method based on the conventional image processing filtering algorithm (median filtering, Gaussian filtering and the like) can possibly cause the peak value and the wave width of the waveform to be changed remarkably, and the filtering result directly influences the accuracy and the reliability of waveform parameter calculation, so that the research and development of the filtering method capable of effectively keeping the LiDAR waveform data characteristics have important significance and practical value.

The distribution of noise in LiDAR signals is quite complex, and it is impractical to denoise a signal by selecting a cutoff frequency using only conventional digital filters. The signal waveform after de-noising by discrete wavelet transform is similar to the wavelet function waveform, the distortion is serious, the wavelet multi-scale decomposition reconstruction calculation amount is large, and if the method is used for de-noising echo of thousands of times of LiDAR, the consumed time space is too large, so that the performance of the method can be used for solving the practical problem only after further research and improvement.

As for the denoising of the space domain, if a better denoising effect is needed, the prior knowledge of the LiDAR signals is added when the support vector machine is trained, if the LiDAR waveforms are not the superposition of Gaussian distribution, the denoising result is not ideal, and in fact the LiDAR waveform signals are not Gaussian distribution. The traditional image space domain processing filtering algorithms do not need the prior knowledge of the waveform, but the algorithms can distort the waveform while denoising.

1. Improved addition-subtraction alternative iterative correction denoising method

(1) Addition and subtraction alternative correction denoising method

For a one-dimensional signal z, it is represented as a column vector (z)₁，…，z_m)^tThe simple form of filtering is expressed as equation 1:

z_i′＝z_i+aΔz_iformula 1

Wherein

z ═ I-aW z formula 2

Wherein W matrix is as follows:

z ═ g (W) z formula 4

If iterate M times, the result is as follows:

z^N＝g(W)^N formula 5

The W matrix is symmetrical and has real eigenvalue and eigenvector, and the real eigenvalue of the W matrix is set to be more than or equal to 0 and less than or equal to W₁≤w₂≤…≤w_mN is less than or equal to n, and the corresponding feature vectors are v1, …, v_mFormula 4 is represented by formula 6:

g (w) ═ 1-aw) (1-cw) formula 7

z_i′＝z_i+cΔz_iformula 8

when the gaussian smoothing is performed by using the formula 1, the addition operation is performed, which causes the image deformation, and if the subtraction operation of the formula 9 is performed, the two calculations are performed alternately, which overcomes the defect of the gaussian filtering.

(2) Improved alternative iteration correction denoising method

z_i′＝z_i+aΔz_iwhen i is an odd number;

z_i′＝z_i+cΔz_iwhen i is an even number;

when the operand is an even number, the following operations are performed:

z_i′＝z_i+aΔz_iwhen i is an even number;

z_i′＝z_i+cΔz_iwhen i is an odd number;

each processing of the addition and subtraction alternative correction denoising method is weighted correction with a single direction, and correction directions of adjacent points are different, so that although the difference between the point and the neighborhood average value is small after multiple processing, the difference between the adjacent two points is possibly large, and the noise elimination strength is insufficient. The improved alternative iterative correction denoising method adopts alternative addition and subtraction correction denoising along with the increase of the processing times, namely, in the processing of adjacent times, the modification of the same point is in different directions, so that the deformation introduced by the method can be counteracted by the iterative processing, and the denoising effect is further improved.

2. De-filtering experiment

In order to evaluate the denoising effect of each filtering algorithm, experiments and comparison are carried out by using real waveform data. The adopted real waveform data is collected by a Leica ALS80 system and stored in an LAS1.4 format, the sampling number of each waveform is 128, the sampling interval is 1 nanosecond, and FIG. 2 shows the result of filtering and denoising the real LiDAR waveform data by respectively using wavelet denoising, Gaussian filtering, plus-minus alternative correction denoising and improved alternative iterative correction denoising when a filtering window is 5

As can be seen from fig. 2, the mean filtering and the gaussian filtering cause significant waveform distortion and the peak value is severely reduced. FIG. 3 shows the denoising effect of the addition and subtraction alternative correction denoising filtering method and the improved alternative iteration correction denoising method by amplifying the peak part in FIG. 2, the denoised waveform of the improved alternative iteration correction denoising method is closer to the original waveform, in order to quantitatively evaluate the denoising effect, a peak signal to noise ratio (PSMR) is introduced to compare the two methods, the peak signal to noise ratio of the real waveform after denoising is shown in FIG. 4, and the peak signal to noise ratio of the waveform after denoising by the improved alternative iteration correction denoising method is larger. The characteristic factors before and after the de-noising of the real waveform are shown in FIG. 5, the wavelet de-noising, the median filtering and the Gaussian filtering cause more serious wave crest shrinkage, the addition and subtraction alternate correction de-noising filtering method only enables the wave crest size to have smaller shrinkage, and the improved alternate iterative correction de-noising method better maintains the wave crest size.

Two, LiDAR full waveform decomposition

LiDAR full waveform data is formed by the superposition of several single Gaussian waves, represented by equation 10:

where k is the number of single Gaussian waves, g_j(z) is the probability density function of the jth Gaussian distribution, d_jIs the size of the wave crest, c_jIs the mean value, i.e. the peak position, σ_jThe standard deviation is half wave width, the number of the sub-waveforms, the size of the wave peak, the position of the wave peak and the optimal value of the half wave width characteristic factor are calculated by waveform decomposition, so that the superposition result of the sub-waveforms is closer to the original waveform.

The flow of LiDAR full waveform decomposition is shown in FIG. 6. Firstly, preprocessing waveform data, removing background noise and random noise in the waveform data, namely the content of the previous part, then estimating the number of sub-waveforms and characteristic factors of each sub-waveform, wherein the characteristic factors comprise the size of a wave crest, the position of the wave crest and half wave width, and finally optimizing the estimated initial value of the characteristic factor by using a dynamic global adjustment method to obtain an accurate sub-waveform characteristic factor.

(one) estimating initial characteristic factor

The estimated characteristic factors comprise the number of estimated sub-waveforms and the characteristic factors of all sub-waveforms, and the specific flow of the method for estimating the initial parameters of the LiDAR waveforms which are extracted one by one is shown in figure 7. Firstly, the original waveform data is preprocessed, background noise and random noise are removed to obtain new waveform data L_mDetecting L_mThe maximum value in the random noise is used as the wave peak value D of the sub-waveform, the position of the maximum value is used as the wave peak position e of the sub-waveform, and in order to avoid the influence of random noise as much as possible, D needs to be more than 3 times of the error in the random noise, namely more than 3 sigma_mOtherwise, the waveform data L is considered_mIf the position r of the sampling point with half of the wave peak value D is known, the noise remains_gThe half-wave width σ of the sub-waveform can be found due to the discreteness of the waveform data sampling points, and r_gIs not necessarily exactly at the sampling point, so r cannot be directly obtained_gR is obtained by statistically analyzing the sampled data on both sides of the peak_gThen, an approximation of σ is calculated.

r_gThe analysis calculation obtaining method comprises the following steps: firstly judging the waveform data on the right side of the wave crest, if the value of the w-1 waveform data sampling point is greater than D/2 and the value of the w +1 waveform data sampling point is less than D/2, then r_gThe approximate value of the wave data is the position of the w-th wave data sampling point, then the approximate value of sigma is calculated, then the wave data on the left side of the wave crest is judged, if the value of the w-1 th wave data sampling point is less than D/2 and the value of the w +1 th wave data sampling point is more than D/2, then r_gThe approximate value of the sigma is calculated by taking the position of the w-th waveform data sampling point; when the waveform data on the two sides are judged to be finished, if a waveform data sampling point meeting the condition is not searched on one side, the r searched on the other side is selected_gTo calculate an approximation of σ as a final result; if the waveform data sampling points on both sides of the wave crest can beCan search to satisfy the condition r_gThen r searched from the left and right sides is compared_gAnd taking a smaller sigma value as a final sigma approximate value according to the calculated sigma approximate value:

after the approximation of σ is calculated, the three initial values of the characteristic factors of the predicted sub-waveform are recorded (D, e, σ), and this sub-waveform is subtracted from the original waveform data, as shown in step (c) of fig. 7. Then repeating the estimation steps for the residual waveform data, and continuously subtracting the estimated sub-waveform from the waveform data until the residual waveform does not satisfy L in step (a) of FIG. 7_max＞3σ_mAnd (4) indicating that only noise remains in the current waveform data, and finishing the estimation of the sub-waveform characteristic factor.

(II) dynamically and globally adjusting optimized characteristic factors

After estimating the number of the sub-waveforms and the initial values of the waveform characteristic factors, optimizing the estimated initial values of the characteristic factors to enable the Gaussian waves fitted by the sub-waveforms to be closest to the original waveform data. The optimization is to give an initial value of a characteristic factor and an objective function, the optimization of parameters is realized by minimizing the objective function, the problem of waveform decomposition of airborne LiDAR full waveform data belongs to a multi-dimensional nonlinear optimization problem, and algorithms commonly used for solving the problem comprise a gradient method, a quasi-Newton method, an LM algorithm and an EM algorithm. The LM algorithm solves the problem of search failure when the Jacobian matrix of the Gauss-Newton algorithm is a non-column matrix, is used for nonlinear programming and least square curve fitting, and is easy to obtain a local optimal solution. The evolution method adopts multi-point parallel search, generates new individuals through crossing and variation in each iteration process, and continuously enlarges the search range, so that the evolution algorithm is easy to search out the global optimal solution rather than the local optimal solution.

1. Optimization strategy of the application

Because the characteristic factor optimization is a pure mathematical calculation process, the phenomena that the fitting of an optimization result and an original waveform is good, but the difference between the optimization result and an actual situation is large, such as the serious reduction of a wave peak value and the large deviation of the wave peak position, for example, the wave peak position and the wave peak size of the waveform in fig. 8(a) are close to the estimated first sub-waveform, but the wave peak position and the wave peak size are greatly changed after the optimization of the LM algorithm in fig. 8(b), the reduction of the wave peak size and the half wave width are taken as a characteristic of the subsequent classification, and the classification precision can be influenced, and the deviation of the wave peak position can seriously influence the coordinate precision of the generated point cloud.

Aiming at the problems, the waveform characteristic factor optimization process is improved, based on the characteristic factor estimation method which is extracted one by one, the estimated first sub-waveform is a beam of pulse with the largest wave peak value, namely the strongest echo signal, the strongest signal is less influenced by noise, the detection result is closer to the real situation, and the later estimated sub-waveform is more easily influenced by noise due to the weaker signal with the smaller wave peak value, so that the change of the wave peak size and the wave peak position of the first sub-waveform is controlled in the process of optimizing the characteristic factor by adopting a dynamic global algorithm, and the optimization result of the first sub-waveform cannot deviate from the actual situation due to too large change.

2. Dynamic global optimization feature factor

The specific process of dynamically and globally optimizing the characteristic factors comprises the following steps:

Z(r)＝(Z1(r)，Z2(r)，…Z_m(r)) formula 12

u_ij(r+1)＝z_best，j(r)+G(z_p1，j(r)-z_p2，j(r)) formula 13

Wherein u is_ij(r +1) is the jth component of the ith individual of the (r +1) th generation, z_best，j(r) is the best individual vector in the r-th generation, G is a variation factor, z_p1，j(r)-z_p2，j(r) is a difference vector, i, p1, p2 are different from each other and p1, p2 are [1, N]Any two random integers,

i

1, 2, and N,

j

1, 2, m;

wherein g (Z)_i(r)) is the fitness of the ith generation of the ith individual;

Z(r+1)＝(Z₁(r+1)，Z₂(r+1)，…，Z_M(r +1)) formula 16

(III) generating high-density point cloud data

The characteristic factors of the sub-waveforms are extracted through waveform decomposition of LiDAR waveform data, the characteristic factors comprise peak positions, peak sizes and half wave widths, a point data recording format comprises three parameters X (r), Y (r) and Z (r), the three parameters calculate space three-dimensional coordinates of points on corresponding waveforms, and the space coordinates of the points on the waveforms are calculated by a formula 17:

x, Y, Z is the spatial location of the point calculated from the sub-waveform feature factors, X₀，Y₀，Z₀Is the position of the starting point, i.e. the three-dimensional coordinates of the point data (parameters X, Y, Z in the point recording format), X (r), Y (r), Z (r) are the velocity vector of the laser signal, r is the time relative to the starting point (the starting point time is 0), i.e. the peak position of the sub-waveform, in picoseconds, the product of X (r), Y (r), Z (r) and r is the displacement of the laser pulse from the starting time to the time when the peak is detected, the position of the starting point plus this displacement yields the coordinates of the point corresponding to the sub-waveform, and the units of X, Y, Z are the units of the LAS data coordinate system.

The selected LAS data is urban data, the point cloud data is provided with 264156 points in total, and the newly generated point cloud data after waveform decomposition comprises 426239 points. The point cloud provided by the system and the point cloud obtained by calculating the waveform decomposition parameter are shown in fig. 9(a) and 9 (b). In order to quantify the contrast density, sample areas with different types of ground features are selected from fig. 9(a) and 9(b) respectively to calculate the point cloud density. Fig. 10 compares the point cloud densities of the four sample areas in fig. 9, and it can be obtained from fig. 10 that since the ground does not include a three-dimensional structure, there is generally only one echo, and only when the laser foot is hit on a target such as a street lamp, an automobile, or the like, there are two or more echoes, so that the point cloud density increase amplitude on the ground is small; when the laser points are hit on the edge of a building, the laser points can be reflected by the roof and the ground of the building, laser signals can penetrate through vegetation canopies to reach trunks, if vegetation is not dense, part of the laser signals can reach the ground to form at least three echoes, and therefore the range of the point cloud density after waveform decomposition is increased greatly in a building area and a vegetation area.

Point cloud classification based on waveform characteristic factors

(I) Classification method

The method is characterized in that point cloud data are classified based on a support vector machine, sample data are classified by four kernel functions in a LibSVM toolbox respectively, and finally, a radial basis function with high speed and high classification precision is selected as a kernel function of the SVM, SVM classification parameters are selected by adopting the dynamic global algorithm, the classification precision is used as a fitness function, and classification parameters are optimized. The point cloud classification process based on the waveform characteristic factors is shown in fig. 11. Given training samples and corresponding categories, using the dynamic global algorithm of the application to optimize classifier parameters by taking classification precision as a fitness function, and then using the classifier to classify the data of the whole measuring area.

(II) experiments and analyses

1. Training of classifiers

The embodiment has 12213 sample areas selected, wherein the number of ground points 4229, the number of building points 6378, and the number of tree points 1606, because the dimensions of several characteristic factors are different, before training and classification, normalization operation is performed on five characteristic factor values of amplitude, peak position, half-wave width, intensity, and elevation, training sets are randomly selected from samples according to different proportions to train a classifier, initial values of a penalty coefficient S and a parameter σ of a radial basis function are randomly selected in an initialization process, a variable factor G is respectively selected as 0.55 and a cross probability S is selected as 0.9 for parameters of a dynamic global algorithm, then the overall classification accuracy is used as a fitness function to optimize S and σ, and finally the overall classification accuracy and the optimal classification parameter are obtained. The classification precision of the classifiers trained by training sets with different proportions is over 97.9 percent.

2. Point cloud classification

The point cloud data of the whole embodiment is classified by using the trained classifier, and the classification result is shown in fig. 12. It can be seen that the classification effect of the ground points and the non-ground points is very obvious, and the classification effect of the buildings and the trees is also improved, and in order to evaluate the classification effect of each ground feature, the classification results of the three ground features are displayed separately, as shown in fig. 13. The LiDAR point cloud classification obtained by analyzing the classification result is accurate and efficient, and particularly the classification accuracy of the ground, trees and buildings is greatly improved.

Claims

The point cloud classification method is characterized in that full-wave type data are subjected to waveform decomposition based on Gaussian attributes of airborne LiDAR full-wave type data, waveform characteristic factors are estimated and a one-by-one extraction method is adopted, the optimized characteristic factors adopt a dynamic global algorithm, LiDAR waveform characteristic factors are extracted and high-density point cloud is generated, and finally the point cloud is classified by utilizing waveform decomposition parameters and point cloud height through a support vector machine; the method mainly comprises the following steps:

first, LiDAR full waveform data preprocessing: based on the characteristics of noise in the LiDAR full waveform data, a method for removing background noise and random noise in the LiDAR full waveform data is provided, on the basis of verifying the denoising superiority of the LiDAR waveform data by the addition and subtraction alternative correction denoising method, the addition and subtraction alternative correction denoising method is improved, the improved addition and subtraction alternative iterative correction denoising method is provided, the random noise is removed, the shape of the waveform is kept, and the signal-to-noise ratio of the waveform data is improved;

second, LiDAR full waveform decomposition: the waveform decomposition process comprises pre-estimating sub-waveform characteristic factors and optimizing characteristic factors, and in the pre-estimating characteristic factor part, the characteristic factors of the sub-waveforms are continuously pre-estimated from original waveform data by adopting a one-by-one extraction method until the residual waveforms only contain noise points; the optimized waveform characteristic factor introduces a dynamic global algorithm and controls the change of a peak value in the optimization process, so that the maximum peak value can not emit violent change in the optimization process; finally, generating high-density point cloud data by using parameters obtained by waveform decomposition;

thirdly, point cloud classification based on waveform characteristic factors: on the basis of analyzing the ground feature reflection characteristics reflected by the waveform characteristic factors, an SVM classifier is adopted, the point cloud is classified by utilizing the waveform characteristic factors obtained by waveform decomposition and point cloud coordinates calculated by utilizing waveform parameters, in the classifier training process, classification accuracy is taken as a fitness function, a penalty coefficient and kernel function parameters are optimized by utilizing a dynamic global algorithm, wherein the kernel function selects a radial basis function, the dynamic global algorithm is adopted for optimizing the classification parameters in the SVM classifier training process, and finally the classifier trained by sample data is used for classifying the point cloud data of the whole measuring area.
2. The LiDAR full-waveform controlled decomposition driven point cloud classification method of claim 1, wherein the improved addition-subtraction alternating iterative correction denoising method: for a one-dimensional signal z, it is represented as a column vector (z)₁，…，z_m)^tThe simple form of filtering is expressed as equation 1:

z_i′＝z_i+aΔz_iformula 1

Wherein
a is a scale factor greater than 0 and less than 1, and the equation is written in matrix form, as in equation 2:

z ═ I-aW z formula 2

Wherein W matrix is as follows:

replacing the matrix I-aW with a non-shrinking function on the matrix W, the following result:

z ═ g (W) z formula 4

If iterate M times, the result is as follows:

z^N＝g(W)^Nformula 5

The W matrix is symmetrical and has real eigenvalue and eigenvector, and the real eigenvalue of the W matrix is set as0≤w₁≤w₂≤…≤w_mN is less than or equal to n, and the corresponding feature vectors are v1, …, v_mFormula 4 is represented by formula 6:

where g is the filter kernel, h_iIs a coefficient, satisfies that after M iterations, when w belongs to [0, n ∈]Time, low frequency component g (w)_i)^M1, high frequency component g (w)_i)^N0, the kernel function takes the form of equation 7:

g (w) ═ 1-aw) (1-cw) formula 7

Where c is a new negative scale factor and c < -a, which is equivalent to a further approximation step after gaussian smoothing with a positive scale factor by equation 1:

z_i′＝z_i+cΔz_iformula 8

Since g (0) ═ 1, a + c < 0, a threshold value w is sought_QLet g (w)_Q) 1, a, c satisfy formula 9:

when the Gaussian smoothing is carried out by using the formula 1, the subtraction operation of the formula 9 is added, and the two calculations are alternately carried out, so that the defect of the Gaussian filtering is overcome.
3. The LiDAR full-waveform controlled decomposition driven point cloud classification method of claim 2, characterized by an improved alternating iterative correction denoising method:

the addition and subtraction alternate correction denoising method modifies a kernel function on the basis of Gaussian filtering, adopts an addition and subtraction alternate calculation method to inhibit contraction generated by the Gaussian filtering, but because c is less than-a, the addition and subtraction degrees are different, the denoising result can be influenced, the addition and subtraction alternate correction denoising method is further improved, multiple alternate operations are performed on the addition and subtraction alternate correction denoising, and the following operations are performed when the operand is an odd number:

z_i′＝z_i+aΔz_iwhen i is an odd number;

z_i′＝z_i+cΔz_iwhen i is an even number;

when the operand is an even number, the following operations are performed:

z_i′＝z_i+aΔz_iwhen i is an even number;

z_i′＝z_i+cΔz_iwhen i is an odd number;

each processing of the plus-minus alternative correction denoising method is weighted correction with a single direction, and correction directions of adjacent points are different, so that although the difference between the point and the neighborhood average value is small after multiple times of processing, the difference between the adjacent two points is likely to be large, and the noise elimination strength is insufficient; the improved alternative iterative correction denoising method adopts alternative addition and subtraction correction denoising along with the increase of the processing times, namely, in the processing of adjacent times, the modification of the same point is in different directions, so that the deformation introduced by the method can be counteracted by the iterative processing, and the denoising effect is further improved.
4. The LiDAR full waveform control decomposition driven point cloud classification method of claim 1, wherein LiDAR full waveform decomposition: LiDAR full waveform data is formed by the superposition of several single Gaussian waves, represented by equation 10:

where k is the number of single Gaussian waves, g_j(z) is the probability density function of the jth Gaussian distribution, d_jIs the size of the wave crest, c_jIs the mean value, i.e. the peak position, σ_jThe standard deviation, namely the half wave width, the optimal values of the number, the wave peak size, the wave peak position and the half wave width characteristic factor of the sub-waveforms are solved by the waveform decomposition, so that the result of the sub-waveform superposition is closer to the result of the half wave widthOriginal waveform;

LiDAR full waveform decomposition procedure: firstly, preprocessing waveform data, removing background noise and random noise in the waveform data, namely the content of the previous part, then estimating the number of sub-waveforms and characteristic factors of each sub-waveform, wherein the characteristic factors comprise the size of a wave crest, the position of the wave crest and half wave width, and finally optimizing the estimated initial value of the characteristic factor by using a dynamic global adjustment method to obtain an accurate sub-waveform characteristic factor.
5. The LiDAR full-waveform controlled decomposition driven point cloud classification method of claim 1, wherein the estimated initial feature factors are: estimating the number of sub-waveforms and the characteristic factor of each sub-waveform;

the method for estimating the initial parameters of the LiDAR waveforms which are extracted one by one comprises the following specific processes: firstly, the original waveform data is preprocessed, background noise and random noise are removed to obtain new waveform data L_mDetecting L_mThe maximum value in the random noise is used as the wave peak value D of the sub-waveform, the position of the maximum value is used as the wave peak position e of the sub-waveform, and in order to avoid the influence of random noise as much as possible, D needs to be more than 3 times of the error in the random noise, namely more than 3 sigma_mOtherwise, the waveform data L is considered_mIf the position r of the sampling point with half of the wave peak value D is known, the noise remains_gThe half-wave width σ of the sub-waveform can be found due to the discreteness of the waveform data sampling points, and r_gIs not necessarily exactly at the sampling point, so r cannot be directly obtained_gR is obtained by statistically analyzing the sampled data on both sides of the peak_gThen, calculating an approximate value of the sigma;

r_gthe analysis calculation obtaining method comprises the following steps: firstly judging the waveform data on the right side of the wave crest, if the value of the w-1 waveform data sampling point is greater than D/2 and the value of the w +1 waveform data sampling point is less than D/2, then r_gThe approximate value of the wave data is the position of the w-th wave data sampling point, then the approximate value of sigma is calculated, then the wave data on the left side of the wave crest is judged, if the value of the w-1-th wave data sampling point is less than D/2And the value of the w +1 waveform data sampling point is greater than D/2, then r_gThe approximate value of the sigma is calculated by taking the position of the w-th waveform data sampling point; when the waveform data on the two sides are judged to be finished, if a waveform data sampling point meeting the condition is not searched on one side, the r searched on the other side is selected_gTo calculate an approximation of σ as a final result; if waveform data sampling points on two sides of a wave crest can be searched to meet the condition r_gThen r searched from the left and right sides is compared_gAnd taking a smaller sigma value as a final sigma approximate value according to the calculated sigma approximate value:

after the approximate value of sigma is calculated, recording (D, e, sigma) as three characteristic factor initial values of the estimated sub-waveform, subtracting the sub-waveform from the original waveform data, repeating the estimation steps for the residual waveform data, and continuously subtracting the estimated sub-waveform from the waveform data until the residual waveform does not meet L_max＞3σ_mAnd (4) indicating that only noise remains in the current waveform data, and finishing the estimation of the sub-waveform characteristic factor.
6. The LiDAR full-waveform control decomposition-driven point cloud classification method of claim 1, wherein the feature factors are dynamically and globally adjusted to optimize: according to the method, the optimization process of the waveform characteristic factors is improved, based on the characteristic factor estimation method which is extracted one by one, the estimated first sub-waveform is a pulse with the largest wave peak value, namely the strongest echo signal, the strongest signal is less influenced by noise, the detection result is closer to the real situation, and the estimated sub-waveform is more easily influenced by the noise due to the weaker signal with the smaller wave peak value, so that the change of the size and the position of the wave peak of the first sub-waveform is controlled in the process of optimizing the characteristic factors by adopting a dynamic global algorithm, and the optimization result of the first sub-waveform cannot deviate from the actual situation due to too large change.
7. The LiDAR full-waveform controlled decomposition driven point cloud classification method of claim 6, wherein the specific process of dynamically and globally optimizing the feature factors is:

firstly, initializing a characteristic factor: determining the size N of the population, and randomly selecting an initial population:

Z(r)＝(Z1(r)，Z2(r)，…Z_m(r)) formula 12

Wherein Z_i(r) represents the ith individual in the r-th generation, r is 0, Z_i(0) Is an m-dimensional vector;

and a second process, characteristic factor group evolution: for each individual Z in Z (r)_i(r) performing the following operations:

firstly, mutation operation: for randomly selecting two individuals z from characteristic factor group_p1，z_p2The following operations are carried out:

u_ij(r+1)＝z_best，j(r)+G(z_p1，j(r)-z_p2，j(r)) formula 13

Wherein u is_ij(r +1) is the jth component of the ith individual of the (r +1) th generation, z_best，j(r) is the best individual vector in the r-th generation, G is a variation factor, z_p1，j(r)-z_p2，j(r) is a difference vector, i, p1, p2 are different from each other and p1, p2 are [1, N]Any two random integers, i 1, 2, and N, j 1, 2, m;

and step two, cross operation: increasing the diversity of the characteristic factor group, and the calculation process is as follows:

wherein S is cross probability, and S is more than 0 and less than or equal to 1, randk_ijIs [0, 1 ]]Random (i) is [1, n ]]Random integer between, cross-operation guarantees v_ijAt least one component of (r +1) is composed of u_ij(r + 1);

selecting operation: to judge Z_iWhether or not to become next generationMember, will vector V_iAnd a target vector Z_iThe fitness of (c) is compared when V_iFitness ratio Z of_iHigh is selected as the offspring, otherwise Z is directed_iAs a child, the selection operation is calculated in the following manner:

wherein g (Z)_i(r)) is the fitness of the ith generation of the ith individual;

and a third process, namely termination judgment: let the new population resulting from process two be:

Z(r+1)＝(Z₁(r+1)，Z₂(r+1)，…，Z_M(r +1)) formula 16

The most optimal individual in Z (r +1) is designated as Z_best(r +1), when the result meets the precision requirement or the whole evolution process has reached the maximum algebra, the operation is terminated and Z is added_bestAnd (r +1) is output as an optimal solution, otherwise, r is made to be r +1, and the process is switched to a second process.
8. The LiDAR full-waveform controlled decomposition driven point cloud classification method of claim 1, generating high density point cloud data: the characteristic factors of the sub-waveforms are extracted through waveform decomposition of LiDAR waveform data, the characteristic factors comprise peak positions, peak sizes and half wave widths, a point data recording format comprises three parameters X (r), Y (r) and Z (r), the three parameters calculate space three-dimensional coordinates of points on corresponding waveforms, and the space coordinates of the points on the waveforms are calculated by a formula 17:

x, Y, Z is the spatial location of the point calculated from the sub-waveform feature factors, X₀，Y₀，Z₀Is the position of the starting point, i.e., the three-dimensional coordinates of the point data, X (r), Y (r), Z (r) are the velocity vectors of the laser signals, and r is the time relative to the starting point, i.e., the peak position of the sub-waveformThe product of the units of picosecond, X (r), Y (r), Z (r) and r is the displacement of the laser pulse from the starting time to the time when the peak is detected, the position of the starting point plus the displacement obtains the coordinates of the point corresponding to the wavelet, and the units of X, Y and Z are the units of the LAS data coordinate system.
9. The LiDAR full waveform control decomposition driven point cloud classification method of claim 1, wherein the point cloud classification based on waveform characterization factors: the method comprises the steps of classifying point cloud data based on a support vector machine, classifying sample data by using four kernel functions in a LibSVM toolbox respectively, selecting a radial basis function with high speed and high classification precision as a kernel function of an SVM, optimizing classification parameters by using a dynamic global algorithm of the application and using the classification precision as a fitness function, giving training samples and corresponding classes, optimizing classifier parameters by using the dynamic global algorithm of the application and using the classification precision as the fitness function, and classifying the data of the whole test area by using the classifier.