CN109886895A - Local fit filtering method - Google Patents
Local fit filtering method Download PDFInfo
- Publication number
- CN109886895A CN109886895A CN201910151831.1A CN201910151831A CN109886895A CN 109886895 A CN109886895 A CN 109886895A CN 201910151831 A CN201910151831 A CN 201910151831A CN 109886895 A CN109886895 A CN 109886895A
- Authority
- CN
- China
- Prior art keywords
- point
- echo
- data
- grid
- ground
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
The present invention provides local fit filtering methods, specific steps include: that initial ground seed point obtains, the acquisition of encryption ground seed point, merges initial ground seed point and encryption ground seed point, and ground seed point cloud data are exported, seed point cloud data filtering program in three-dimensional laser ground is completed in conjunction with Weighted Average Algorithm finally by minimum range filtering.Filtering algorithm of the invention not only can be reduced data volume but also can guarantee filter quality.
Description
Technical field
The present invention relates to technical field of data processing, especially local fit filtering method.
Background technique
The laser point cloud filtering of broad sense refers to from a certain data acquisition system, and the number useful to research is extracted according to correlation principle
According to rejecting the process of extraneous data.For the point cloud data acquired in the laser radar system, filtering is generally referred to laser
The earth's surfaces non-ground points such as artificial construction of structures, forest cover in point cloud data are rejected, to obtain the mistake of practically millet cake cloud
Journey.
The problem of laser radar scanning system is in terms of hardware and the system integration has been basically solved with the development of science and technology, but
Its Data Post technology is that current LiDAR technology can't widely applied key problem also in conceptual phase.?
There are two major issue in LiDAR data last handling process, one is that influence of the systematic error to three-dimensional laser pin point coordinate should
How to eliminate, another is exactly the filtered classification of three-dimensional point cloud.For LiDAR 3 d-dem laser point converges P, system energy
Directly acquiring the D coordinates value Pi (xi, yi, zi) of each laser footpoint, (i=1,2,3 ... n), and it is strong to additionally include echo
The additional informations such as degree, echo times.It is converged in laser point and contains ground point, culture point and rough error point in P, and LiDAR data
Post-processing main purpose be exactly to obtain digital terrain model, this just need to isolate by certain method atural object point set and
Rough error point set only retains ground point collection.Atural object point set and rough error point set are removed, the process for obtaining digital terrain model is exactly radar
Point cloud data filtering.
Current filtering method is based on height value, by choosing ground seed point or seed region compared with neighborhood mostly
Perhaps gradient etc. achievees the effect that filtering also with good grounds echo strength, more echo informations or several method knot for elevation, height difference
It closes to complete filtering operation.In recent years, the research of filtering method achieves certain achievement, but that there are still filtering accuracies is not high,
The problems such as data volume is big does not have blanket filtering algorithm yet, becomes restriction laser radar surface sweeping technology and further develops
Bottleneck.
Summary of the invention
In view of the foregoing drawbacks and problem, the object of the present invention is to provide local fit filtering methods, solve existing laser
Filtering algorithm precision in radar system Data Post is not high, error is big, and data volume is big, so as to cause data-handling efficiency
It is low, and cannot achieve blanket technical problem.
In order to achieve the above object, the invention provides the following technical scheme:
A kind of local fit filtering method, specifically includes the following steps:
Step 1: initial ground seed point obtains: grid is established in region according to single and last echo point cloud density,
Grid size is determined by the minimum amount threshold that must contain in cloud density and single grid；It is selected respectively in each grid
Take minimum point as initial ground seed point；
Step 2: the acquisition of encryption ground seed point: grid being re-established according to current point cloud density, for each
The initial ground point O determined in step 1i(xi, yi, zi), search for single and last echo point cloud undetermined in its neighborhood
Rj(xj, yj, zj), wherein contiguous range can be arranged according to ground seed point grid is obtained；
According to the initial ground point Oi(xi, yi, zi) and its neighborhood in point to be located Rj(xj, yj, zj) coordinate calculates its slope
Degree；If the gradient is in threshold range, then it is assumed that the last echo point is that encryption ground seed point is recognized if the gradient exceeds threshold value
It is point to be judged for the last echo point；
Step 3: after obtaining encryption ground seed point by step 2, by itself and the initial ground seed point that is obtained in step 1
Merge, by traversal check, removal repeat point, export ground seed point cloud data and have not determined single echo point cloud,
Last echo point cloud data；
Step 4: by the output ground seed point cloud data and the single echo point cloud having not determined, last echo point
Cloud data are filtered by minimum range and complete seed point cloud data filtering program in three-dimensional laser ground in conjunction with Weighted Average Algorithm；
Wherein, the implementation procedure of the minimum range filtering algorithm specifically:
Step 101: acquisition data simultaneously pre-process data: then progress abnormality value removing first carries out zero-mean
Secondly processing is removal trend item parts, be finally that extracting cycle item forms time series data；
Step 102: time series data pretreated in step 1 being modeled and adjusts model parameter: specifically including: building
Erection system model: the parameter Estimation of time series linear model refers on the basis of identification obtains model classification and order, asks
The numerical value of autoregressive coefficient in model and sliding average coefficient out；The applicability of testing model: time series linear model
Applicability, which is examined, refers to whether be applicable in by the model that above method determines with the data detection of sample, i.e. inspection residual sequence is
No is white noise sequence；
Step 103: Kalman filter being designed to each system model and application is based on the limited of minimum vector distance
Model algorithm carries out online switching in real time to model；It specifically includes: Kalman filtering is designed to the system model that each is established
Device, given system noise and observation noise battle array on the basis of the system model picked out；To each system model according to base
Carry out the accuracy of assessment system model in minimization vector distance criterion fitting evaluation function；It will be to be selected according to the evaluation function
The filter result of the highest system model of accuracy is exported as final result in the model selected；
Wherein, the Weighted Average Algorithm specifically:
Step 201: assuming that there is N number of unit element, carrying out operation by data weighted average first；
Step 202: if input value is K, lucky unit element N-K, N-K+1 ... .N-1, N are selected, then counting
Number device C1 work, count results add 1, otherwise keep original count result；
Step 203: when counter C1 count value reaches preset value, another counter C2 work, count results add
1, while counter C1 is set to 0, otherwise counter C2 keeps original count result；
Step 204: if the count value of counter C2 reaches identical element number of packages N, counter C2 count results are set to
0；
Step 205: barrel shifter is using the count value of counter C2 as shift amount, by the operation knot in the first step
Fruit ring shift right obtains new result and rotation is gone to select N number of unit element.
Further, in step 1, the initial ground seed point is obtained specifically:
(1) all last echo point clouds and single echo point cloud are traversed first, obtain the X in point cloud datamax、Ymax、
Xmin、Ymin；
(2) then pass through last echo point cloud and single echo point cloud quantity N1, Experimental Area area S, obtain point Yun Ping
Equal density p1, ρ1=N1/ S determines grid size by a cloud averag density；
(3) Experimental Area is divided according to the grid size of setting；
(4) point cloud data in each grid is detected, point cloud data is ranked up by height value, obtains elevation minimum
Value；
(5) all grid are traversed, the smallest cloud of height value in each grid is exported as initial ground seed point.
Further, in step 2, the acquisition of the encryption ground seed point specifically:
(1) according to initial ground seed point and single echo undetermined, last echo point cloud quantity N2, Experimental Area area S,
Obtain a cloud averag density ρ2, ρ2=N2/ S redefines grid size by a cloud averag density；
(2) initial ground seed point cloud is traversed；
(3) for each initial ground seed point cloud, all single echoes undetermined and last echo in grid are searched for
Point cloud, calculates separately the gradient of seed point and each point to be located cloud in grid；
(4) gained value of slope and gradient threshold value comparison will be calculated, if being less than gradient threshold value, output is encryption ground seed
Point wouldn't then be handled if more than gradient threshold value；
(5) all grid are traversed, obtain all encrypting ground seed point.
Further, in the step 2, the calculation formula of the gradient are as follows:
In formula, xi, yi, ziFor initial ground point three-dimensional coordinate, xj, yj, zjFor point to be located three-dimensional coordinate.
Further, in step 4, the filter specifically:
(1) input data: non-ground seed point cloud number in input ground seed point cloud data and single echo, last echo
According to；
(2) all data, the X found out respectively are traversedmax、Ymax、Xmin、YminAnd point cloud sum N3, find out test block
The minimum outsourcing rectangular area S in domain, thus finds out the packing density ρ of discrete point3=N3/S；
(3) each point to be located institute is determined by seed point cloud number needed for ground seed dot density, Weighted Average Algorithm
The grid side length L for needing to search for；
(4) by each point to be located Qi (xi, yi, zi), horizontal plane is projected to, then centered on each point to be located, is built
Grid is found, then searches for the ground seed point in its neighborhood by grid side length；
(5) the maximum height difference dh in the grid centered on point to be located between the seed point of ground is calculated；
If meeting: 1) height difference is even variation and is located on two opposite side of grid, 2) set interpolated error as h meters of △,
Calculate ds=(△ h*L)/dh, there are the actual range of ground seed point and point to be located is small in the grid centered on point to be located
In or equal to ds；
Minimum distance method is then used, directly replaces weighting flat with the height value of the ground seed point nearest apart from point to be located
Equal height value；
If not: then using weighted mean method interpolation: the n ground kind searched for each point to be located Qi (xi, yi, zi)
(j=1,2 ... n), its height value is obtained the meter of interpolation elevation ZQ, ZQ using Weighted Average Algorithm by son point Mj (xj, yj, zj)
Calculate formula are as follows:
Wherein, the elevation that Zj is j-th point in n ground seed point is j-th point in n ground seed point of weight,
It is defined as the inverse of point-to-point transmission horizontal distance square, it may be assumed that
(6) compare the height value Z of point to be located QiElevation Z is weighted and averaged with its neighborhood seed pointQIf Zi-ZQGreater than height difference threshold
Value, then it is assumed that Q point is vegetation point, if not, it is believed that the point to be located is ground point；
(7) single and last echo point cloud for judging each point to be located, are classified as two class of vegetation point and ground point.
Further, in step 3 further include total data to acquisition, handled by following steps:
S301, the total data that the step 3 is obtained and geographical location carry out standardization matching treatment；
S302, the data of the point of proximity acquisition in matched treated data and peripheral location preset range are carried out pair
Than being labeled as falling point when difference is more than threshold value between correlation data；
S303, all falling points are stored in database profession convenient for later retrieval analysis, inquiry, are shown.
Local fit filtering method of the invention carries out echo free to all data first, will return for the first time with intermediate time
Wave is classified as vegetation point cloud, then divides two-stage to choose ground seed point last and single echo point cloud, is finally with point to be located
Ground seed point is searched at center, completes to filter by minimum range combination Weighted Average Algorithm.Only to last and single echo point
Cloud filtering, not only can be reduced data volume but also can guarantee filter quality, and substantially increased the operation efficiency of filtering.
Detailed description of the invention
In order to more clearly explain the embodiment of the invention or the technical proposal in the existing technology, to embodiment or will show below
There is attached drawing needed in technical description to be briefly described, it should be apparent that, the accompanying drawings in the following description is only this
Some embodiments of invention for those of ordinary skill in the art without creative efforts, can be with
It obtains other drawings based on these drawings.
Fig. 1 is local fit filtering method step schematic diagram of the invention.
Specific embodiment
Below in conjunction with attached drawing of the invention, technical solution of the present invention is clearly and completely described, it is clear that institute
The embodiment of description is only a part of the embodiment of the present invention, instead of all the embodiments.Based on the embodiments of the present invention,
Every other embodiment obtained by those of ordinary skill in the art without making creative efforts, belongs to this hair
The range of bright protection.
Laser point cloud data is the superposition of more echo datas, and each echo information is all a kind of embodiment of earth's surface information, repeatedly
Echo information is a big advantage of LiDAR system.In more echo informations of LiDAR, the data volume of usually intermediate time echo than
For the first time with the data volume much less of last echo, three times and the above echo data amount is than twice return much smaller number.
The echo information of LiDAR can be divided into single echo (single echo) according to whether laser beam occurs multiple reflections
With multiecho (multiple-echo), multiecho can be divided into echo (first echo), last echo (last for the first time again
Echo the centre time echo) and between first and last time.Single echo refers to that same laser beam receives only an echo information.
For earth's surface, building top surface etc. only exist the object to be measured in individual reflection face in vertical distribution, usually only
Single echo.Multiecho refers to that same laser beam touches the various pieces of object to be measured until depleted of energy, every part point
Other reflection echo information and successively recorded by system and the echo information that is formed.Laser radar system is anti-by receive first
Penetrate signal and be known as echo for the first time, the last one the reflection signal received is known as last echo, between for the first time with last echo
Between reflection signal be known as intermediate echo.For vegetative coverage region, building facade or edge etc. in spatial vertical
There are for the object to be measured of multiple reflectings surface, usually occurring multiecho in distribution, usual first and last time echo can exist than
Biggish difference, and intermediate echo relative different is smaller.Therefore, that the echo information of LiDAR is divided into multiple subsets is more advantageous
In being filtered.
1, by each echo overlay analysis:
First last echo information overlay analysis: echo almost all is superimposed on the top of last echo for the first time, illustrates for the first time
A possibility that echo is substantially the non-ground points such as Vegetation canopy, trunk, and last echo is ground point is maximum.
First, last, single respectively with intermediate echo information overlay analysis:
Echo is superimposed with intermediate echo information for the first time: intermediate time echo data almost all is located under echo information for the first time
Just and it is close to echo information for the first time, can release echo for the first time is the echo information that Vegetation canopy reflects mostly, intermediate time
Echo is mostly the echo information that the limb of short vegetation or tall and big trees is reflected.
Last echo is superimposed with intermediate echo information: it is close that last echo point cloud density is significantly greater than intermediate echo point cloud
Degree, intermediate time echo almost all is contained among last echo, and part last echo point cloud is located above intermediate echo,
Most last echo is located at below intermediate echo, and last echo can be speculated for most of ground return and small part vegetation branch
It is dry to reflect obtained echo information.
Single echo is superimposed with intermediate echo information: single echo point cloud density with last echo as, hence it is evident that be greater than
Intermediate echo point cloud density, and intermediate echo almost all is contained among single echo, unlike, it is mostly single
Secondary echo point cloud is located above intermediate echo, and small part single echo is located at below intermediate echo.It will thus be seen that single
Echo is the echo information reflected by small part ground, most vegetation Malabar Pied Hornbill and part vegetation limb.
2, the selection based on first, last time echo information filtering
LiDAR data is to need to determine to return based on echo for the first time or last based on the major issue that more echoes filter
Wave filtering.Because laser beam can penetrate vegetation, it is more likely in research using last echo as ground point, but last echo
It is also likely to be that other short vegetation generate, is easy therefore to cause error.
Because LiDAR system integration error and multipath effect, there may be more scattered in last echo-signal
Low spot, although carrying out the rejecting of rough error and low spot before filtering, excessive low spot also will affect filtering accuracy, and echo for the first time
Middle low spot data volume is fewer than last echo, so most of can select to be filtered based on echo-signal for the first time.Based on for the first time
The filtering method of echo-signal does not make full use of the penetration capacity of laser, can from the corresponding type of ground objects of the more echoes of following table
It arrives, echo may be that vegetation or building may be blocked due to large area vegetation in the forest zone that forest is intensive and lack ground for the first time
Millet cake, to cause the loss of DTM precision.
Also there is the research of comprehensive utilization first and last twice return filtering, this method can be in morphologic filtering and linear prediction
It is accomplished in the parts filtering algorithms such as filtering, and it is not suitable for other algorithm, such as it is applied to irregular triangle network
Filtering, will not only enhance filter effect, cause biggish error due to also filtering statistical data because changing, also increase data
Amount, reduces data-handling efficiency instead.
In summary analyze: the Project Areas intensive for vegetative coverage, is more likely ground based on last echo information
In contrast the fundamemtal phenomena of information, is based on last in the case where can guarantee has enough ground laser footpoints in Project Areas
It is a good selection that echo, which is filtered LiDAR data,.Intermediate, Er Qieji relatively fewer with echo information amount for the first time
All vegetation area point cloud datas improve filtration efficiency to reduce data volume, direct with echo information for the first time by intermediate time
It is classified as non-ground points collection, filtering is no longer participate in and calculates.
Therefore, the present invention will be filtered based on last echo.To optimize the filter effect based on last echo, the present invention
It is classified using grid and chooses ground seed point.
As shown in Figure 1, specific filtering algorithm the following steps are included:
Step S1: initial ground seed point obtains: according to single and last echo point cloud density the built-in legislate in region then
Grid chooses minimum point as initial ground seed point respectively in each grid.
Specifically:
(1) all last echo point clouds and single echo point cloud are traversed first, obtain the X in point cloud datamax、Ymax、
Xmin、Ymin；
(2) then pass through last echo point cloud and single echo point cloud quantity N1, Experimental Area area S, obtain point Yun Mi
Spend ρ, ρ=N1/ S determines regular grid size by cloud density；
(3) Experimental Area is divided according to the regular grid size of setting；
(4) point cloud data in each grid is detected, point cloud data is ranked up by height value, obtains elevation minimum
Value；
(5) all grid are traversed, the smallest cloud of height value in each grid is exported as initial ground seed point.
Step S2: the acquisition of encryption ground seed point: the initial ground point O determined in step 1 for eachi
(xi, yi, zi), search for single and last echo point cloud R undetermined in its neighborhoodj(xj, yj, zj), wherein contiguous range can basis
Obtain the setting of ground seed point grid；
According to the initial ground point Oi(xi, yi, zi) and its neighborhood in point to be located Rj(xj, yj, zj) coordinate calculates its slope
Degree；If the gradient is in threshold range, then it is assumed that the last echo point is that encryption ground seed point is recognized if the gradient exceeds threshold value
It is point to be judged for the last echo point.
Specifically:
(1) initial ground seed point cloud is traversed；
(2) for each initial ground seed point cloud, all single echoes undetermined and last in regular grid are searched for
Echo point cloud calculates separately the gradient of seed point and each point to be located cloud in grid；
(3) gained value of slope and gradient threshold value comparison will be calculated, if being less than gradient threshold value, output is encryption ground seed
Point wouldn't then be handled if more than gradient threshold value；
(4) all grid are traversed, obtain all encrypting ground seed point.
Wherein, the calculation formula of the gradient are as follows:
In formula, xi, yi, ziFor initial ground point three-dimensional coordinate, xj, yj, zjFor point to be located three-dimensional coordinate.
Step S3: after obtaining encryption ground seed point by step 2, by itself and the initial ground seed that is obtained in step 1
Point merges, and is checked by traversal, removal repeats point, the single echo point for exporting ground seed point cloud data and having not determined
Cloud, last echo point cloud data.
Wherein, step 3 further includes the total data to acquisition, is handled by following steps:
S301, the total data that step 3 is obtained and geographical location carry out standardization matching treatment；
S302, the data of the point of proximity acquisition in matched treated data and peripheral location preset range are carried out pair
Than being labeled as falling point when difference is more than threshold value between correlation data；
S303, all falling points are stored in database profession convenient for later retrieval analysis, inquiry, are shown.
Step S4: the filter of three-dimensional laser ground seed point cloud data is completed into minimum range filtering in conjunction with Weighted Average Algorithm
Wave-path sequence.
Minimum range filtering algorithm and Weighted Average Algorithm are explained first.
(1) minimum range filtering algorithm:
The implementation procedure of algorithm are as follows:
Step 1: acquiring data and being pre-processed to data:
Abnormality value removing is carried out first, then carries out zero-mean processing, is secondly removal trend item parts, is finally to extract
Periodic term forms time series data.
Step 2: being modeled to time series data pretreated in step 1 and adjusting model parameter:
2.1 establish system model: the parameter Estimation of time series linear model, which refers to, obtains model classification and rank in identification
On the basis of number, the numerical value of the autoregressive coefficient and sliding average coefficient in model is found out；
The applicability of 2.2 testing models: the applicability of time series linear model examines the data detection referred to sample
Whether it is applicable in by the model that above method determines, i.e. whether inspection residual sequence is white noise sequence.
Step 3: designing Kalman filter and application based on the finite module for minimizing vector distance to each system model
Type algorithm carries out online switching in real time to model:
3.1 pairs of each system models established design Kalman filter, on the basis of the system model picked out
Given system noise and observation noise battle array；
3.2 pairs of each system models are fitted evaluation function according to based on minimization vector distance criterion come assessment system mould
The accuracy of type；
3.3 according to the above-mentioned evaluation function provided by the filtering knot of the highest system model of accuracy in model to be selected
Fruit exports as final result.
(2) Weighted Average Algorithm:
Step 1: assuming there is N number of unit element, operation is carried out by traditional data Weighted Average Algorithm first；
Step 2: if when input value is K, lucky unit element N-K, N-K+1 ... .N-1, N are selected, then counting
Device C1 work, count results add 1, otherwise keep original count result；
Step 3: another counter C2 work, count results add 1 when counter C1 count value reaches preset value,
Counter C1 is set to 0 simultaneously, otherwise counter C2 keeps original count result；
Step 4: counter C2 count results are set to 0 if the count value of counter C2 reaches identical element number of packages N；
Step 5: barrel shifter is using the count value of counter C2 as shift amount, by the operation knot in the first step
Fruit ring shift right obtains new result and rotation is gone to select this N number of unit element.
The specific work process of step S4 are as follows:
(1) input data: non-ground seed point cloud number in input ground seed point cloud data and single echo, last echo
According to；
(2) all data, the X found out respectively are traversedmax、Ymax、Xmin、YminAnd point cloud sum N3, find out test block
The minimum outsourcing rectangular area S in domain, thus finds out the packing density ρ of discrete point3=N3/S；
(3) each point to be located institute is determined by seed point cloud number needed for ground seed dot density, Weighted Average Algorithm
The grid side length L for needing to search for；
(4) by each point to be located Qi (xi, yi, zi), horizontal plane is projected to, then centered on each point to be located, is built
Grid is found, then searches for the ground seed point in its neighborhood by grid side length；
(5) the maximum height difference dh in the grid centered on point to be located between the seed point of ground is calculated；；
(6) compare the height value Z of point to be located QiElevation Z is weighted and averaged with its neighborhood seed pointQIf Zi-ZQGreater than height difference threshold
Value, then it is assumed that Q point is vegetation point, if not, it is believed that the point to be located is ground point；
(7) single and last echo point cloud for judging each point to be located, are classified as two class of vegetation point and ground point.
Wherein, specific filter specifically includes that
(1) according to Experimental Area area S, total ground seed point number N2, obtain ground seed point density p2=N2/S；
(2) each point to be located institute is determined by seed point cloud number needed for ground seed dot density, Weighted Average Algorithm
The grid side length L for needing to search for, such as to reach has 8 ground seed points in averagely each grid, then
(3) firstly for each point to be located Qi(xi, yi, zi), it is projected into horizontal plane, then with each point to be located
Centered on, regular grid is established, then the ground seed point in its neighborhood is searched for by grid side length；
(4) the maximum height difference d in the regular grid centered on point to be located between the seed point of ground is calculatedh,
If meeting: 1) height difference is even variation and is located on two opposite side of grid, 2) set interpolated error as h meters of △,
Calculate ds=(△ h*L)/dh, there are the actual range of ground seed point and point to be located is small in the grid centered on point to be located
In or equal to ds；
Minimum distance method is then used, directly replaces weighting flat with the height value of the ground seed point nearest apart from point to be located
Equal height value；
If not: then using weighted mean method interpolation: the n ground kind searched for each point to be located Qi (xi, yi, zi)
(j=1,2 ... n), its height value is obtained the meter of interpolation elevation ZQ, ZQ using Weighted Average Algorithm by son point Mj (xj, yj, zj)
Calculate formula are as follows:
Wherein, the elevation that Zj is j-th point in n ground seed point is j-th point in n ground seed point of weight,
It is defined as the inverse of point-to-point transmission horizontal distance square, it may be assumed that
(6) compare the height value Z of point to be located QiElevation Z is weighted and averaged with its neighborhood seed pointQIf Zi-ZQGreater than height difference threshold
Value, then it is assumed that Q point is vegetation point, if not, it is believed that the point to be located is ground point；
(7) judge each single undetermined and last echo point cloud, be classified as two class of vegetation point and ground point.
Precision analysis evaluation finally is carried out to using filtering algorithm result of the invention, " refusing true " error is 3.62%, " is received
Puppet " error is 1.45%, overall error 2.16%, and filtering overall effect is good.In addition, the present invention during filtering, has
29144 are put cloud with intermediate time for the first time, i.e. the 29% of total amount of data is not engaged in filtering algorithm, illustrates that this research filtering method exists
Guarantee that filtration efficiency has been effectively ensured in the case where filter quality.
Wherein, " refuse true " error, which refers to, is considered error caused by non-ground points for the ground point in reference data, " receives
Puppet " error, which refers to, is considered that error caused by ground point, overall error are " refusing true " errors for the non-ground points in reference data
The synthesis of " receive puppet " error, is able to reflect the accuracy of algorithm.
The above description is merely a specific embodiment, but scope of protection of the present invention is not limited thereto, any
Those familiar with the art in the technical scope disclosed by the present invention, can easily think of the change or the replacement, and should all contain
Lid is within protection scope of the present invention.Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
Claims (6)
1. a kind of local fit filtering method, it is characterised in that: specifically includes the following steps:
Step 1: initial ground seed point obtains: grid, grid are established in region according to single and last echo point cloud density
Size is determined by the minimum amount threshold that must contain in cloud density and single grid；It is chosen respectively most in each grid
Low spot is as initial ground seed point；
Step 2: the acquisition of encryption ground seed point: grid being re-established according to current point cloud density, for each
The initial ground point O determined in step 1i(xi, yi, zi), search for single and last echo point cloud R undetermined in its neighborhoodj
(xj, yj, zj), wherein contiguous range can be arranged according to ground seed point grid is obtained；
According to the initial ground point Oi(xi, yi, zi) and its neighborhood in point to be located Rj(xj, yj, zj) coordinate calculates its gradient；If
The gradient is in threshold range, then it is assumed that the last echo point is encryption ground seed point, if the gradient exceeds threshold value, then it is assumed that the end
Secondary echo point is point to be judged；
Step 3: after obtaining encryption ground seed point by step 2, it being closed with the initial ground seed point obtained in step 1
And checked by traversal, removal repetition point, export ground seed point cloud data and the single echo point cloud having not determined, end
Secondary echo point cloud data；
Step 4: by the ground seed point cloud data exported in step 3 and the single echo point cloud, the last echo that have not determined
Point cloud data is filtered by minimum range and completes seed point cloud data filtering journey in three-dimensional laser ground in conjunction with Weighted Average Algorithm
Sequence；
Wherein, the implementation procedure of the minimum range filtering algorithm specifically:
Step 101: acquisition data simultaneously pre-process data: then progress abnormality value removing first carries out zero-mean processing,
Secondly it is removal trend item parts, is finally that extracting cycle item forms time series data；
Step 102: time series data pretreated in step 1 being modeled and adjusts model parameter: being specifically included: establishing system
System model: the parameter Estimation of time series linear model refers on the basis of identification obtains model classification and order, finds out mould
The numerical value of autoregressive coefficient and sliding average coefficient in type；The applicability of testing model: time series linear model is applicable in
Property examine refer to the data detection of sample by above method determine model whether be applicable in, i.e., inspection residual sequence whether be
White noise sequence；
Step 103: Kalman filter and application are designed based on the finite model for minimizing vector distance to each system model
Algorithm carries out online switching in real time to model；It specifically includes: Kalman filter is designed to the system model that each is established,
Given system noise and observation noise battle array on the basis of the system model picked out；To each system model according to based on minimum
Change the accuracy that vector distance criterion fitting evaluation function carrys out assessment system model；According to the evaluation function by mould to be selected
The filter result of the highest system model of accuracy is exported as final result in type；
Wherein, the Weighted Average Algorithm specifically:
Step 201: assuming that there is N number of unit element, carrying out operation by data weighted average first；
Step 202: if input value is K, lucky unit element N-K, N-K+1 ... .N-1, N are selected, then counter
C1 work, count results add 1, otherwise keep original count result；
Step 203: when counter C1 count value reaches preset value, another counter C2 work, count results add 1, together
When counter C1 is set to 0, otherwise counter C2 keeps original count result；
Step 204: if the count value of counter C2 reaches identical element number of packages N, counter C2 count results are set to 0；
Step 205: barrel shifter follows the operation result in the first step using the count value of counter C2 as shift amount
Ring moves to right to obtain new result and rotation is gone to select N number of unit element.
2. local fit filtering method according to claim 1, it is characterised in that: in step 1, the initial ground seed
Point obtains specifically:
(1) all last echo point clouds and single echo point cloud are traversed first, obtain the X in point cloud datamax、Ymax、Xmin、
Ymin；
(2) then pass through last echo point cloud and single echo point cloud quantity N1, Experimental Area area S, obtain a cloud averag density
ρ1, ρ1=N1/ S determines grid size by a cloud averag density；
(3) Experimental Area is divided according to the grid size of setting；
(4) point cloud data in each grid is detected, point cloud data is ranked up by height value, obtains elevation minimum value；
(5) all grid are traversed, the smallest cloud of height value in each grid is exported as initial ground seed point.
3. local fit filtering method according to claim 1, it is characterised in that: in step 2, encryption ground seed
The acquisition of point specifically:
(1) according to initial ground seed point and single echo undetermined, last echo point cloud quantity N2, Experimental Area area S, obtain
Point cloud averag density ρ2, ρ2=N2/ S redefines grid size by a cloud averag density；
(2) initial ground seed point cloud is traversed；
(3) for each initial ground seed point cloud, all single echoes undetermined and last echo point cloud in grid are searched for,
Calculate separately the gradient of seed point and each point to be located cloud in grid；
(4) gained value of slope and gradient threshold value comparison will be calculated, if being less than gradient threshold value, output is encryption ground seed point,
If more than gradient threshold value, then wouldn't handle；
(5) all grid are traversed, obtain all encrypting ground seed point.
4. local fit filtering method according to claim 1, it is characterised in that: in the step 2, the meter of the gradient
Calculate formula are as follows:
In formula, xi, yi, ziFor initial ground point three-dimensional coordinate, xj, yj, zjFor point to be located three-dimensional coordinate.
5. local fit filtering method according to claim 1, it is characterised in that: in step 4, the filter is specific
Are as follows:
(1) input data: non-ground seed point cloud data in input ground seed point cloud data and single echo, last echo；
(2) all data, the X found out respectively are traversedmax、Ymax、Xmin、YminAnd point cloud sum N3, find out Experimental Area most
Thus small outsourcing rectangular area S finds out the packing density ρ of discrete point3=N3/S；
(3) it is determined required for each point to be located by seed point cloud number needed for ground seed dot density, Weighted Average Algorithm
The grid side length L of search；
(4) by each point to be located Qi (xi, yi, zi), horizontal plane is projected to, then centered on each point to be located, establishes lattice
Net, then search for by grid side length the ground seed point in its neighborhood；
(5) the maximum height difference dh in the grid centered on point to be located between the seed point of ground is calculated；
If meeting: 1) height difference is even variation and is located on two opposite side of grid, 2) interpolated error is set as h meters of △, it calculates
Ds=(△ h*L)/dh out, be less than in the grid centered on point to be located there are the actual range of ground seed point and point to be located or
Person is equal to ds；
Minimum distance method is then used, directly replaces weighted average high with the height value of the ground seed point nearest apart from point to be located
Journey value；
If not: then using weighted mean method interpolation: the n ground seed point searched for each point to be located Qi (xi, yi, zi)
(j=1,2 ... n), and the calculating that its height value is obtained interpolation elevation ZQ, ZQ using Weighted Average Algorithm is public by Mj (xj, yj, zj)
Formula are as follows:
Wherein, the elevation that Zj is j-th point in n ground seed point is j-th point in n ground seed point of weight, definition
For the inverse of point-to-point transmission horizontal distance square, it may be assumed that
(6) compare the height value Z of point to be located QiElevation Z is weighted and averaged with its neighborhood seed pointQIf Zi-ZQGreater than height difference threshold value, then
Think that Q point is vegetation point, if not, it is believed that the point to be located is ground point；
(7) single and last echo point cloud for judging each point to be located, are classified as two class of vegetation point and ground point.
6. local fit filtering method according to claim 1, which is characterized in that step 3 further includes the whole to acquisition
Data are handled by following steps:
S301, the total data that step 3 is obtained and geographical location carry out standardization matching treatment；
S302, the data of matched treated data and the point of proximity acquisition in peripheral location preset range are compared,
Falling point is labeled as when difference is more than threshold value between correlation data；
S303, all falling points are stored in database profession convenient for later retrieval analysis, inquiry, are shown.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910151831.1A CN109886895A (en) | 2019-02-28 | 2019-02-28 | Local fit filtering method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910151831.1A CN109886895A (en) | 2019-02-28 | 2019-02-28 | Local fit filtering method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109886895A true CN109886895A (en) | 2019-06-14 |
Family
ID=66930003
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910151831.1A Pending CN109886895A (en) | 2019-02-28 | 2019-02-28 | Local fit filtering method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109886895A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110660027A (en) * | 2019-08-28 | 2020-01-07 | 青岛秀山移动测量有限公司 | Laser point cloud continuous profile ground filtering method for complex terrain |
| CN114548277A (en) * | 2022-02-22 | 2022-05-27 | 电子科技大学 | Method and system for fitting ground points and extracting crop height based on point cloud data |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101969307A (en) * | 2010-08-20 | 2011-02-09 | 浙江大学 | Improved data weighed averaging algorithm and device |
| CN102679984A (en) * | 2012-05-29 | 2012-09-19 | 北京理工大学 | Finite model filtering method based on vector distance minimizing criterion |
-
2019
- 2019-02-28 CN CN201910151831.1A patent/CN109886895A/en active Pending
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101969307A (en) * | 2010-08-20 | 2011-02-09 | 浙江大学 | Improved data weighed averaging algorithm and device |
| CN102679984A (en) * | 2012-05-29 | 2012-09-19 | 北京理工大学 | Finite model filtering method based on vector distance minimizing criterion |
Non-Patent Citations (1)
| Title |
|---|
| 彭丽: "植被密集的陡坡区机载 LiDAR数据滤波方法研究", 《中国优秀硕士学位论文全文数据库农业科技辑》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110660027A (en) * | 2019-08-28 | 2020-01-07 | 青岛秀山移动测量有限公司 | Laser point cloud continuous profile ground filtering method for complex terrain |
| CN110660027B (en) * | 2019-08-28 | 2023-03-31 | 青岛秀山移动测量有限公司 | Laser point cloud continuous profile ground filtering method for complex terrain |
| CN114548277A (en) * | 2022-02-22 | 2022-05-27 | 电子科技大学 | Method and system for fitting ground points and extracting crop height based on point cloud data |
| CN114548277B (en) * | 2022-02-22 | 2023-09-08 | 电子科技大学 | Method and system for ground point fitting and crop height extraction based on point cloud data |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wu et al. | Individual tree crown delineation using localized contour tree method and airborne LiDAR data in coniferous forests | |
| CN105488770B (en) | A kind of airborne laser radar point cloud filtering method of object-oriented | |
| Gibbs et al. | Approaches to three-dimensional reconstruction of plant shoot topology and geometry | |
| Mian et al. | Three-dimensional model-based object recognition and segmentation in cluttered scenes | |
| CN110473178A (en) | A kind of open defect detection method and system based on multiple light courcess fusion | |
| CN102103202B (en) | Semi-supervised classification method for airborne laser radar data fusing images | |
| CN109118265A (en) | Commercial circle determines method, apparatus and server | |
| CN106780551B (en) | A kind of Three-Dimensional Moving Targets detection method and system | |
| Bienert et al. | Voxel space analysis of terrestrial laser scans in forests for wind field modelling | |
| CN105389849B (en) | Vehicle collision angle resolved systems based on three-dimensional reconstruction | |
| CN109767464A (en) | A kind of point cloud registration method of low Duplication | |
| CN109886895A (en) | Local fit filtering method | |
| CN109933539A (en) | A kind of Software Defects Predict Methods based on principal component analysis and combination sampling | |
| Schilling et al. | Tree topology representation from TLS point clouds using depth-first search in voxel space | |
| CN113869629A (en) | Laser point cloud-based power transmission line safety risk analysis, judgment and evaluation method | |
| CN106295498A (en) | Remote sensing image target area detection apparatus and method | |
| CN108038081A (en) | The landslide disaster logistic regression analysis of feature based function space filter value | |
| CN110363299A (en) | Space reasoning by cases method towards delamination-terrane of appearing | |
| CN113450269A (en) | Point cloud key point extraction method based on 3D vision | |
| CN113947770B (en) | Method for identifying object placed in different areas of intelligent cabinet | |
| Wijayanto et al. | Machine learning approaches using satellite data for oil palm area detection in Pekanbaru City, Riau | |
| CN105956544A (en) | Remote sensing image road intersection extraction method based on structural index characteristic | |
| Perrin et al. | Multisensor fusion in the frame of evidence theory for landmines detection | |
| Gong et al. | Point cloud segmentation of 3D scattered parts sampled by RealSense | |
| Ma et al. | Land covers classification based on Random Forest method using features from full-waveform LiDAR data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190614 |