An optimal point cloud density prediction method based on airborne lidar

Abstract

The present invention discloses an optimal point cloud density prediction method based on airborne laser radar, comprising: S1, obtaining high-density original point clouds and orthophotos through airborne laser radar, and performing preprocessing to obtain a data source; S2, using the data source to make a data set for training; S3, using the data set to train an ICP-NN model; S4, inputting the high-density original point clouds and orthophotos to be processed into the trained ICP-NN model to obtain a prediction interval for the optimal point cloud density. The model trained by the present application can directly provide a collection interval for airborne laser radar point cloud data collection in complex mountainous areas.

Description

Technical Field

[0001] The present invention relates to the field of point cloud processing, and in particular to an optimal point cloud density prediction method based on airborne laser radar. Background Art

[0002] Airborne LiDAR technology can quickly and accurately obtain surface data over a large area. Through multiple echo data collection methods, it can penetrate vegetation cover to generate a three-dimensional terrain model of the real surface. Therefore, it is widely used in the field of geological disaster investigation. However, there are still many problems with the point cloud density of airborne LiDAR field data collection. The Digital Elevation Model (DEM) constructed by the point cloud data obtained by LiDAR is crucial in the manual visual interpretation of geological disasters and the automatic identification of cracks or disaster scope by machine learning. Its accuracy and visualization quality are directly related to the accuracy of identification and analysis. The size of the point cloud density directly determines the quality of the DEM. It is difficult to accurately judge whether the ground point cloud density can meet the accuracy requirements of DEM construction after the original point cloud collected is processed based on the experience of data collectors. On the one hand, if the density of laser point cloud required for collection is high, data collectors need to set a higher laser point frequency and heading overlap, or a lower flight speed, which will increase the cost of flight, data storage and processing; on the other hand, too low point cloud density will cause local scratching and missing micro-geomorphological elements (such as cracks, gullies, steep slopes, depressions, etc.) when TIN interpolation generates DEM. In addition, from the equipment level, due to the positioning accuracy limitations of the airborne LiDAR system, when the point cloud density is too high and exceeds the positioning accuracy, a large number of noise points will be generated.

[0003] In order to clarify the appropriate point cloud density value, various countries have formulated corresponding airborne LiDAR data collection specifications, but they are not sufficient to support a variety of collection scenarios. The United States Geological Survey (USGS) published the Lidar Base Specification for the surveying and mapping industry in August 2012. Its 2024 update stipulates that the quality of the constructed 3D elevation model must be no less than QL2 level, and the corresponding point cloud density is 2pts / m 2 The highest quality level QL0 requires a point cloud density of 8pts / m 2, the absolute vertical accuracy (RMSE) of the DEM should be no greater than 0.1m. However, this specification does not specify how to determine the quality level according to the survey needs, and implementers cannot infer the original point cloud density that should be collected on site only by the ground point cloud density required to construct the DEM. As one of the countries with frequent geological disasters, the Geospatial Information Authority of Japan (GSI) formulated a drone public survey manual in 2016. After being updated in 2020, the manual stipulates that the point cloud density collected by lidar in places with high vegetation influence needs to reach 20~200pts / m2, but the recommended value range is too wide to accurately guide on-site personnel to collect point cloud density according to different needs. China's existing airborne lidar data acquisition technical specifications were formulated by the surveying and mapping industry in 2011, and stipulated the point cloud density for constructing DEM at different scales. The maximum point cloud density required in the regulations is not less than 16 points / square meter, which is to meet the recognition needs of contour lines in the surveying and mapping field, and does not consider the recognition needs of micro-geomorphological elements under vegetation coverage. At present, there is still a lack of scientific and complete guidance on the density of airborne LiDAR point cloud acquisition for remote sensing surveys of geological hazards.

[0004] In fact, most of the current research is aimed at the ground point cloud that constitutes the DEM, and no appropriate point cloud density has been proposed. At present, there are few studies on the density of point clouds collected to guide production. It is a relatively complex issue to set the ideal optimal ground point cloud density. First, the laser penetration depends greatly on the performance of the lidar sensor itself; second, the vegetation coverage greatly affects how many of the emitted point clouds can reach the ground; finally, the elevation change in the survey area will cause the real ground to be at a different distance from the sensor, that is, when the sensor is closer to the ground, the number of laser point clouds received per unit area will be greater; the slope change in the survey area will affect the size of the area receiving the laser, that is, when the slope is steeper, the number of laser point clouds received on the steep surface will be less. Therefore, there are many factors to consider when inferring the point cloud density that needs to be collected on site from the ideal ground point cloud density for constructing DEM. This application introduces a machine learning method. Taking into account the influence of forest canopy closure, slope, and height difference in the survey area, the ideal optimal ground point is used as input to predict the collection point cloud density interval that should be designed in the drone-borne LiDAR data collection stage. Because the vegetation cover, slope, etc. in the survey area fluctuate around the mean value, uncertainty measures should be provided for the predicted values. Summary of the invention

[0005] In view of the above-mentioned deficiencies in the prior art, the present invention provides a method for solving the problem that the existing point cloud density prediction method lacks uncertainty measurement.

[0006] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is: an optimal point cloud density prediction method based on airborne laser radar, comprising:

[0007] S1, obtain high-density original point cloud and orthophoto through airborne laser radar, and perform preprocessing to obtain data source;

[0008] S2, using the data source to create a data set for training;

[0009] S3. Use the dataset to train the ICP-NN model.

[0010] S4. Input the high-density original point cloud and orthophoto to be processed into the trained ICP-NN model to obtain the prediction interval of the optimal point cloud density.

[0011] Further: S1 includes:

[0012] S11, collect data through airborne laser radar to obtain high-density original point cloud and orthophoto;

[0013] S12, importing POS data into orthophotos, and using aerial triangulation to obtain target points;

[0014] S13, constructing a network according to the target points to obtain a processed orthophoto;

[0015] S14, performing band stitching and band adjustment on the high-density original point cloud to obtain a processed point cloud;

[0016] S15, automatically classifying the processed point cloud and filtering out noise to obtain a ground point cloud;

[0017] S16, constructing a digital elevation model based on the ground point cloud;

[0018] S17. Use the processed orthophoto and digital elevation model as data sources.

[0019] Further: S2 includes:

[0020] S21, extracting the number of point clouds and the number of ground point clouds in the data source;

[0021] S22, calculating the surface vegetation coverage of the data source;

[0022] S23. Create a data set for training based on the geographical environment of the data source, the surface vegetation coverage, the number of point clouds, and the number of ground point clouds.

[0023] Furthermore: in S22, vegetation density is used to characterize the surface vegetation coverage.

[0024] Further: S3 includes:

[0025] S31, training the neural network model through the data set to obtain a trained neural network model;

[0026] S32. Use the inductive commonality prediction algorithm to assist the neural network and obtain the ICP-NN model.

[0027] Further: In S4, the method for the ICP-NN model to obtain the prediction interval of the optimal point cloud density includes:

[0028] S41. Set the prediction interval execution degree of the ICP-NN model , calculate the unlabeled data x The prediction set ;

[0029] S42, according to the prediction set The constraint of Transformed into the prediction interval of the optimal point cloud density.

[0030] Further: In S41, the prediction set The expression is:

[0031]

[0032] in, H (.) is the Heaviside step function, For data points The NCM value of and For two machine learning models of the same structure and type, and is the weight set in the trained neural network model, For data Z j The corresponding NCM value is Data Z j Corresponding input variables and target labels, data Z j is an independent and identically distributed data set Any data in x is the unlabeled data, y For unlabeled data x The corresponding true target value is R is a real number, l and k Both represent range nodes for a data group.

[0033] Further: the independent and identically distributed data set middle, is the training set, is the calibration set, For the test set.

[0034] Further: In said S42, the prediction set The constraints are:

[0035]

[0036] Among them, Pr(.) represents the joint distribution.

[0037] Further: The prediction interval of the optimal point cloud density is ;

[0038] in, , , q For all data in the calibration set Quantile.

[0039] The beneficial effects of the present invention are:

[0040] 1. It provides uncertainty measurement, which enables decision makers to better control risks;

[0041] 2. Improve the interval prediction algorithm and use it for the design of point cloud density collected by airborne laser radar, and be able to specify the prediction interval of the confidence level;

[0042] 3. A dataset based on geographical environment, vegetation density, ground point cloud density and collected point cloud density was established, which can be used to train other point cloud density prediction models. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 This is the process of the optimal point cloud density prediction method based on airborne lidar.

[0044] Figure 2 Schematic diagram of data collection and preprocessing.

[0045] Figure 3 Schematic diagram of the ICP-NN model.

[0046] Figure 4 This is a canopy density distribution diagram for a certain study area in the embodiment.

[0047] Figure 5 This is the loss curve of the ICP-NN model.

[0048] Figure 6 This is the MPIW distribution diagram of the ICP-NN model under different confidence levels CL. DETAILED DESCRIPTION

[0049] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0050] like Figure 1 As shown, in one embodiment of the present invention, an optimal point cloud density prediction method based on airborne laser radar is provided, comprising:

[0051] S1, obtain high-density original point cloud and orthophoto through airborne laser radar, and perform preprocessing to obtain data source;

[0052] S2, using the data source to create a data set for training;

[0053] S3. Use the dataset to train the ICP-NN model.

[0054] S4. Input the high-density original point cloud and orthophoto to be processed into the trained ICP-NN model to obtain the prediction interval of the optimal point cloud density.

[0055] like Figure 2 As shown, the S1 includes:

[0056] S11, collect data through airborne laser radar to obtain high-density original point cloud and orthophoto;

[0057] S12, importing POS data into orthophotos, and using aerial triangulation to obtain target points;

[0058] S13, constructing a network according to the target points to obtain a processed orthophoto;

[0059] S14, performing band stitching and band adjustment on the high-density original point cloud to obtain a processed point cloud;

[0060] S15, automatically classifying the processed point cloud and filtering out noise to obtain a ground point cloud;

[0061] S16, constructing a digital elevation model based on the ground point cloud;

[0062] S17. Use the processed orthophoto and digital elevation model as data sources.

[0063] The S2 includes:

[0064] S21, extracting the number of point clouds and the number of ground point clouds in the data source;

[0065] S22, calculating the surface vegetation coverage of the data source;

[0066] S23. Create a data set for training based on the geographical environment of the data source, the surface vegetation coverage, the number of point clouds, and the number of ground point clouds.

[0067] In S22, the vegetation density is used to characterize the surface vegetation coverage.

[0068] Specifically, the number of collected point clouds and the number of ground point clouds per unit area were extracted using the PointCloud Density Calculator toolbox in the FME software;

[0069] The calculation method of vegetation canopy density is:

[0070] In ArcGIS, by subtracting DSM (processed point cloud) from DEM (digital elevation model), the grid with a height greater than the threshold is assigned a value of 1, indicating that there is vegetation cover within the grid range, otherwise it is assigned a value of 0. Then, the number of grids with a value of 0 and 1 in the statistical unit is calculated and normalized.

[0071] The geographical environment includes slope and elevation. The extraction method is: use ArcGIS software to extract the slope and elevation according to statistical units, and subtract the minimum elevation value from each elevation value to obtain the elevation difference, which is used to represent the degree of elevation fluctuation in the survey area.

[0072] The S3 includes:

[0073] S31, training the neural network model through the data set to obtain a trained neural network model;

[0074] S32. Use the inductive commonality prediction algorithm to assist the neural network and obtain the ICP-NN model.

[0075] like Figure 3 As shown, in S4, the method for obtaining the prediction interval of the optimal point cloud density by the ICP-NN model includes:

[0076] S41. Set the prediction interval execution degree of the ICP-NN model , calculate the unlabeled data x The prediction set ;

[0077] S42, according to the prediction set The constraint of Transformed into the prediction interval of the optimal point cloud density.

[0078] In S41, the prediction set The expression is:

[0079]

[0080] in, H (.) is the Heaviside step function, For data points The NCM value of and For two machine learning models of the same structure and type, and is the weight set in the trained neural network model, For data Z j The corresponding NCM value is Data Z j Corresponding input variables and target labels, data Z j is an independent and identically distributed data set Any data in x is the unlabeled data, y For unlabeled data x The corresponding true target value is R is a real number, l and k Both represent range nodes for a data group.

[0081] The independent and identically distributed data set middle, is the training set, is the calibration set, For the test set.

[0082] In S42, the prediction set The constraints are:

[0083]

[0084] Among them, Pr(.) represents the joint distribution.

[0085] The prediction interval of the optimal point cloud density is ;

[0086] in, , , q For all data in the calibration set Quantile.

[0087] In one embodiment of the present invention, a study area is selected. The study area is located in a typical tectonic erosion low mountain landform area, and the surface vegetation types are varied, including shrubs, trees, grasslands, etc. Partial sliding occurred in the middle of the study area due to continuous rainfall. The landslide has clear boundaries, multiple cracks on the rear edge, multiple gullies on the front edge, and steep slopes and secondary sliding in the middle. Because there are many micro-geomorphological elements and rich geological information in the area, it is suitable as a study area.

[0088] During the data collection phase, a multi-rotor drone equipped with a DV-LiDAR10 laser radar module was used for data collection. During the collection, the laser radar point frequency was adjusted to 200Hz, the lateral overlap was set to 35%, the flight altitude and flight speed were 250m and 10m / s respectively, so that the average collection point cloud density reached 144.84 points / square meter. This embodiment adopts a collection point cloud density far higher than that used in conventional geological disaster remote sensing survey operations to ensure that a wider range of collection point cloud density values ​​can be obtained.

[0089] In order to obtain the required ground point density, the collection density needs to be adjusted. Given the complexity of the survey area during the collection process, the relationship between ground point density and collection point density is nonlinear, not one-to-one, and changes with the changes in the terrain and vegetation in the survey area. According to previous work, in the process of collecting point clouds, the obtained ground point density will also be affected by canopy density, height difference, and slope. Therefore, we propose to use neural networks to estimate the collection density, and use canopy density, height difference, slope, and ground point density as input neural network features to calculate the mapping between ground point density and collection density in different geographical environments. This method involves using a database containing terrain data, real ground point density, and corresponding collection density to train a neural network for prediction. When using it, the user needs to first estimate the geographical environment data (canopy density, height difference, slope), and calculate the optimal ground point density, and then input the obtained data into the trained neural network to obtain the collection density that meets the expected ground point density.

[0090] When using the Inductive Conformal Predictor (ICP) framework to assist a neural network in generating high-quality prediction intervals, you first need to train a general machine learning regression model that can predict label values, and then train a machine learning model of the same structure and type to predict the error of the previous model.

[0091] There are many advantages to using ICP. First, it can theoretically ensure that the PICP of the generated prediction interval meets the confidence level set by the user. Second, compared with the bootstrap method, which usually requires training a large number of models greater than or equal to 50, ICP has very mild requirements on the distribution of data and only needs to train two models, which has higher operating efficiency.

[0092] In this embodiment, the two indicators for measuring the prediction interval are the mean prediction interval width (MPIW) and the prediction interval coverage probability (PICP), which are calculated as follows:

[0093]

[0094] A high-quality prediction interval should meet the following requirements: PICP should be as close as possible to the confidence level preset by the user, and MPIW should be as small as possible.

[0095] like Figure 4 As shown, in this embodiment, the stand canopy density is obtained by the canopy density estimation method of the two-dimensional canopy height model, where 0-0.2 is low canopy density, 0.2-0.69 is medium canopy density, and 0.7 and above is high canopy density (such as Figure 4 b). Actual vegetation coverage compared to DOM Figure 4 a) In contrast, evergreen broad-leaved forests are basically identified as high canopy density, coniferous forests are identified as medium canopy density, and sparse low shrubs, fruit trees and roads are identified as low canopy density. The overall inversion effect is good. At the same time, the height difference, slope value, collection point cloud density and ground point cloud density values ​​of the study area are extracted.

[0096] In this embodiment, after removing NoData, negative values ​​and outliers from the acquired data sources, 94970 sets of data remain as data sets, of which 40% are used as training sets, 20% as calibration sets, 20% as test sets, and the remaining 20% ​​are used to adjust hyperparameters; the slope, elevation difference, canopy density and ground point cloud density of the study area are used as feature values, and the point cloud density is collected as the output value. The neural model composed of 30 hidden layer neurons is used as the regression prediction model, and the learning rate is 10 -4 , the number of iterations reaches 1000, the model fits well, and the loss curve of the training model is as follows Figure 5 .

[0097] Finally, the prediction interval width calculated by the ICP-NN model on the test set is as follows Figure 6 , the figure x The coordinates are the actual density values ​​of the collected point cloud. yThe axis represents the predicted value. When the predicted value is completely equal to the true value, the red prediction curve will coincide with the true value represented by the line with a slope of 45°. Figure 6 The point cloud density collected is 100pts / m 2 ~175pts / m 2 The prediction error is extremely small within the range, and the point cloud density is higher than 300pts / m 2 There will be a large prediction error. And as the confidence level CL increases, the interval width also increases. The MPIW and PICP values ​​under different confidence levels CL are shown in Table 1:

[0098] Table 1 MPIW and PICP of ICP-NN on the test set

[0099] Confidence Level CL MPIW PICP 0.5 36.00 0.497 0.6 43.73 0.601 0.7 52.34 0.705 0.8 61.62 0.798 0.9 70.33 0.899

[0100] Quantifying the uncertainty of point predictions can help decision makers better control risks. Similar algorithms have been used to assist medical decision-making and autonomous driving. In this embodiment, the interval prediction algorithm is improved and used for the design of airborne LiDAR acquisition point cloud density. The conformal algorithm has been tested many times and is considered to be an effective method for providing uncertainty reporting. This algorithm can specify the prediction interval of the confidence level (when there are enough samples, PICP is very close to the set confidence level CL).

[0101] The ICP-NN model is used to take the height difference of the measurement area, forest canopy density, slope, and ideal optimal ground point cloud density as feature values. The predicted value is the collection point cloud density. From the results, the collection point cloud density is 100~175pts / m 2 The predicted value and the actual value fit very well. When the point cloud density is 50~100pts / m 2 or 175~250pts / m 2 The deviation between the predicted result and the actual value gradually increases and is less than 50pts / m 2 Or higher than 250pts / m 2 The prediction effect is poor. Analysis of MPIW distribution shows that it is less than 100 pts / m 2 The predicted point cloud density is always larger than the actual value, and is higher than 250pts / m 2 The predicted value is always lower than the actual value. The reason for analyzing the data set and model may be that the data collected by the researchers is less than 100 pts / m 2 Below and 250pts / m 2 The density of the collected point cloud is relatively low, so the amount of data in the training and test data sets is not large. Therefore, the model pays more attention to 100~250pts / m during training. 2With the increase of confidence, MPIW tends to widen. When the confidence CL is 80%, MPIW reaches 61.62, which means that when predicting the density of the collected point cloud, the probability of covering the true value with 80% must be guaranteed, and the predicted value needs to be expanded by about 30 pts / m. 2 The current wide range is mainly due to 100pts / m 2 Below and 250pts / m 2 The above-mentioned low point cloud density prediction accuracy can be reduced by collecting individual points less than 100 pts / m 2 and above 250 pts / m 2 By supplementing the database with point cloud data, the length of the prediction interval can be further optimized. In short, the current ICP-NN algorithm has a good prediction effect within the commonly used collection point cloud range.

[0102] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. An optimal point cloud density prediction method based on airborne laser radar, characterized in that: include: S1, obtain high-density original point cloud and orthophoto through airborne laser radar, and perform preprocessing to obtain data source; S2, using the data source to create a data set for training; S3. Use the dataset to train the ICP-NN model. Among them, ICP stands for Inductive Conformal Predictor framework; S3 includes: S31, training the neural network model through the data set to obtain a trained neural network model; S32, using the inductive commonality prediction algorithm to assist the neural network, and obtaining the ICP-NN model; S4, inputting the high-density original point cloud and orthophoto to be processed into the trained ICP-NN model to obtain the prediction interval of the optimal point cloud density; In S4, the method for the ICP-NN model to obtain the prediction interval of the optimal point cloud density includes: S41. Set the prediction interval execution degree ε of the ICP-NN model and calculate the prediction set Γ of the unlabeled data x. ε (x); Prediction set Γ ε The expression of (x) is: Where H(.) is the Heaviside step function, A(x,y) = |ym(η;X)| / σ(θ;x) is the NCM value of the data point (x,y), m and σ are two machine learning models with the same structure and type, η and θ are the weight sets in the trained neural network model, and A(x j ,y j ) is the data Z j The corresponding NCM value, x j ,y j The data Z j Corresponding input variables and target labels, data Z j is an independent and identically distributed data set {Z1,Z2,...,Z k ,...,Z l ,...,Z n }, x is unlabeled data, y is the true target value corresponding to the unlabeled data x, R is a real number, l and k both represent the range nodes of the data group; S42, according to the prediction set Γ ε (x) constraint, the prediction set Γ ε (x) is converted into the prediction interval of the optimal point cloud density; Prediction set Γ ε The constraints on (x) are: Pr(y i ∈Γ ∈ (x i ))≥1-ε Among them, Pr(.) represents the joint distribution; The prediction interval of the optimal point cloud density is in, q is the ε quantile of all data in the calibration set.

2. The optimal point cloud density prediction method based on airborne laser radar according to claim 1 is characterized in that: The S1 includes: S11, collect data through airborne laser radar to obtain high-density original point cloud and orthophoto; S12, importing POS data into orthophotos, and using aerial triangulation to obtain target points; S13, constructing a network according to the target points to obtain a processed orthophoto; S14, performing band stitching and band adjustment on the high-density original point cloud to obtain a processed point cloud; S15, automatically classifying the processed point cloud and filtering out noise to obtain a ground point cloud; S16, constructing a digital elevation model based on the ground point cloud; S17. Use the processed orthophoto and digital elevation model as data sources.

3. The optimal point cloud density prediction method based on airborne laser radar according to claim 2 is characterized in that: The S2 includes: S21, extracting the number of point clouds and the number of ground point clouds in the data source; S22, calculating the surface vegetation coverage of the data source; S23. Create a data set for training based on the geographical environment of the data source, the surface vegetation coverage, the number of point clouds, and the number of ground point clouds.

4. The optimal point cloud density prediction method based on airborne laser radar according to claim 3 is characterized in that: In S22, the vegetation density is used to characterize the surface vegetation coverage.

5. The optimal point cloud density prediction method based on airborne laser radar according to claim 1, characterized in that: The independent and identically distributed data set {Z1, Z2, ..., Z k ,...,Z l ,...,Z n }, {Z1,Z2,...,Z k } is the training set, {Z k+1 ,...,Z l } is the calibration set, {Z l+1 ,...,Z n } is the test set.

Images