CN103177417A - GPGPU (general purpose computing on graphics processing unit) based mathematical-morphology LiDAR (Light detection and ranging) point cloud quick-filtering method - Google Patents

Abstract

The invention relates to a GPGPU (general purpose computing on graphics processing unit) based mathematical-morphology LiDAR (Light detection and ranging) point cloud quick-filtering method. The method includes taking a CPU (central processing unit) as a main unit to complete task processing and series computing with high logicality, taking a GPU(graphics processing unit) as a coprocessor to complete parallel task processing with high thread parallelization, utilizing the GPU parallelization technology to realize establishment of a virtual regular grid, gross error reject and mathematical-morphology filtering of the virtual grid of original LiDAR point cloud data, and judge attributes of all points so as to complete filtering. By the method adopting the GPU-based parallelization technology, mathematical morphology operation during filtering is accelerated, operation time is shortened effectively, and filtering is realized quickly.

Description

Mathematical morphology LiDAR point cloud fast filtering method based on GPGPU

Technical field

The invention belongs to the Photogrammetry and Remote Sensing technical field, relate to airborne LiDAR(laser radar, LightDetection And Ranging) measuring technique and GPGPU (general-purpose computations graphic process unit, Gernel PurposeComputing on Graphics Processing Units), relate in particular to mathematical morphology LiDAR point cloud fast filtering method based on GPGPU.

Background technology

The filtering of LiDAR cloud data represents being distributed in of LiDAR system acquisition discrete point on the topographical surface different target, irregular distribution is processed and identified, most basic work is to distinguish with the point (as buildings, vehicle, trees and vegetation etc., being called non-ground point) that drops on other targets of ground being distributed in ground point (being called ground point).Based on traditional Mathematical Morphology theory, the LiDAR cloud data is carried out filtering and compare additive method, owing to not needing to resolve equation, thereby speed is fast, and efficient is higher.But, traditional method may need that original point cloud data is carried out the graticule mesh interpolation and process, existing Mathematical Morphology Method is difficult to according to Different Ground target and different terrain characteristics and automatically adjusts the size of filter window, and require the problem of less and dynamic window for less small hill or hillside landform, traditional mathematical Morphology Algorithm is difficult to take into account, and adopt the serial computing technology, make the consuming time longer of filtering method.

1.GPGPU technology

Traditionally, only GPU is applied to processing graphics plays up calculation task to people, does not fully take into account GPU and has high calculated performance, and this is undoubtedly the significant wastage to computational resource.2003, along with the programmability of GPU improves constantly, the Objective Concept american computer drawing professional association of GPGPU proposed, and made us GPU can be used for the calculating in field beyond graph rendering.This is GPU development milestone the most of historical significance.Also be subject to the restriction of various factors due to the core number of CPU, as technical barriers such as cost, heat radiations.And in recent years,

The performance of GPU has had significant raising, and its general-purpose computations ability is considerably beyond CPU.In the face of the demand of consumer to the image fast processing, Programmable GPU has developed into a kind of high-performance calculation machine platform with outstanding computing power.

2. Airborne LiDAR Technology

Airborne LiDAR Technology refers to the technology of a kind of height integrated laser range finding, dynamic GPS difference and inertial navigation attitude determination.Wherein, laser ranging be used for to be measured the laser radar signal transmitted-reference and is put distance between the laser pin point of ground; Dynamic difference GPS is used for determining the locus of laser radar signal transmitted-reference point; Inertial navigation is used for measuring the primary optical axis attitude parameter of scanister.Synchronous by this three, as to coordinate work has realized directly obtaining of three-dimensional coordinate on a surface target.Compare photogrammetric measurement, it be the active type measure mode, ageing strong, the laser pulse penetration capacity is strong, operating efficiency is high, production cost is low.

Filtering and supporting quality control thereof are that in airborne LiDAR Data Post generating digital elevation model (DEM) process, most critical is also a step the most consuming time, have almost accounted for the time of aftertreatment flow process 60%～80%.

Summary of the invention

The purpose of this invention is to provide a kind of mathematical morphology LiDAR point cloud fast filtering method based on GPGPU, solve conventional serial filtering method long problem consuming time.

For achieving the above object, mathematical morphology LiDAR point cloud fast filtering method technical scheme based on GPGPU of the present invention is as follows: the method is carried out the strong task of logicality with CPU as main frame and is processed and serial computing, carry out the parallel task of height thread parallelization take GPU as coprocessor and process, the step of filtering method is as follows:

(1) the original LiDAR cloud data of input is in the GPU video memory, and then GPU operates below the LiDAR cloud data is distributed to the different threads executed in parallel: calculate the graticule mesh ranks coordinate of LiDAR cloud data, and data allocations is arrived each self-corresponding graticule mesh;

(2) make each block correspondence in the GPU framework calculate the interior data of a graticule mesh, corresponding each the laser pin point that calculates in each graticule mesh of each thread in block is made comparisons according to certain judgment criterion and threshold value, will be judged as the rough error point greater than the laser pin point of threshold value and reject;

(3) according to the graticule mesh after step (2) excluding gross error point, all laser pin points according to its planimetric coordinates, are assigned in each self-corresponding regular grid, set up thus a dummy rules graticule mesh that can index each laser pin point.

(4) with the data copy in described dummy rules graticule mesh to the video memory space, and bind texture memory, different threads is distributed in the some position of all laser pin points in described dummy rules graticule mesh uses respectively and be of a size of l ₁Window each data that are tied to texture memory are carried out mathematical morphology open operator, and judge whether data corresponding to operation result are ground laser pin points, and filtering non-ground laser pin point;

(5) window size and the discrepancy in elevation threshold parameter that import into according to the CPU end increase described window size at the GPU end, and the operation of repeating step (4) is until window size is greater than maximum buildings;

(6) according to the difference of elevation through the laser pin point in the dummy rules graticule mesh after the laser pin point in the graticule mesh after step (2) and process step (4), whether the laser pin point in determining step (2) graticule mesh has and orders the same attribute through the dummy rules graticule mesh laser pin after step (4), so travel through all discrete laser pin points, remove the non-ground laser pin point in graticule mesh in step (2), complete filtering.

In step (1), the minimum height value of each grid unit as representative point, if there is no laser pin point data in grid unit, is got the height value of its neighbor point.

In step (1), the computing formula of graticule mesh ranks coordinate is $\left\{\begin{array}{c}I=(X-X_{\min })/\mathrm{mc}\\ J=(Y-Y_{\min })/\mathrm{mc}\\ c=1/\sqrt{n}\end{array}\right.,$ Wherein, (I, J) is graticule mesh ranks coordinates; (X, Y) is original point cloud planimetric coordinates; X _min, Y _minIt is respectively the minimum value of horizontal ordinate and ordinate in a cloud; C is sampling interval; M is constant, the multiple of expression sampling interval; N is that unit area inner laser pin is counted out.

In step (2), coarse error criterion is: if | V _p-U|〉T, laser pin point p is judged as the rough error point; Otherwise do not think that p is rough error, wherein the elevation difference V of laser pin point p _p=M _p-H _p, M _pTypical value for laser pin point p; H _pHeight value for laser pin point p; Standard deviation $\mathrm{SD}=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(V_i-U)}^2/(n-1)},$ Arithmetic mean

N is total number of laser pin point, and T is threshold value.

Described T is set to 3 times of SD.

In described step (4), the criterion of judgement ground laser pin point is: establish dh _p,1Poor for the elevation before and after laser pin point p opening operation, dh _{T, 1}Be the difference of elevation threshold value of opening operation for the first time, if dh _p,1≤ dh _T,1, the p point is judged as ground laser pin point; Otherwise be non-ground laser pin point; dh _{Max (t), 1}For topographical surface maximum elevation before and after opening operation poor, if dh _{T, 1}Dh _{Max (t), 1}, the part of ground laser pin point can be retained, but not topocentric part can be removed.

Described definite difference of elevation threshold value dh _T,kFormula be $d{h}_{T,k}=\left\{\begin{array}{cc}d{h}_0& w_k\≤3\\ s(w_k-w_{k-1})c+d{h}_0& w_k>3\\ d{h}_{\max }& d{h}_{T,k}>d{h}_{\max }\end{array}\right.,$

dh ₀Be initial difference of elevation threshold value; C is a cloud average headway; dh _maxBe maximum elevation difference limen value, the poor dh of maximum elevation of landform before and after the k time iteration _{Max (t), k}With window size w _kWith the pass of terrain slope s be

In described step (5), the increase mode of window size is w _k=2kb+1 or w _k=2b ^k+ 1, wherein, k is iterations; w _kBe the k time window size; B is the home window size.

In described dummy rules graticule mesh, each graticule mesh has and only has a laser pin point.

The beneficial effect that the present invention reaches: the mathematical morphology LiDAR point cloud fast filtering method based on GPGPU of the present invention, utilization is based on the concurrent technique of GPU, to original serial filtering algorithm the changing that walk abreast, take full advantage of the powerful calculated performance of GPU, accelerate the mathematical morphological operation in filtering, effectively shorten and calculate working time, realize quick filter; Utilized Airborne LiDAR Technology, because the laser pulse penetration capacity is strong, avoided the occlusion issue of vegetation in aviation image institute corresponding ground, and production cost has been low.

Description of drawings

Fig. 1 is CUDA isomery programming model figure;

Fig. 2 is the mathematical morphology filter procedure chart;

Fig. 3 is the filtering process flow diagram;

Fig. 4-1, Fig. 5-1, Fig. 6-1, Fig. 7-1, Fig. 8-1, Fig. 9-1, Figure 10-1, Figure 11-1, Figure 12-1, Figure 13-1, Figure 14-1, Figure 15-1,16-1, Figure 17-1, Figure 18-1st, the gradual change hypsometric map that in the experiment of ISPRS sample, the original point cloud generates;

In Fig. 4-2, Fig. 5-2, Fig. 6-2, Fig. 7-2, Fig. 8-2, Fig. 9-2, Figure 10-2, Figure 11-2, Figure 12-2, Figure 13-2, Figure 14-2, Figure 15-2,16-2, Figure 17-2, Figure 18-2nd, ISPRS sample experiment based on GPGPU parallel algorithm filtering result;

Fig. 4-3, Fig. 5-3, Fig. 6-3, Fig. 7-3, Fig. 8-3, Fig. 9-3, Figure 10-3, Figure 11-3, Figure 12-3, Figure 13-3, Figure 14-3, Figure 15-3,16-3, Figure 17-3, Figure 18-3rd, reference filtering result in the experiment of ISPRS sample.

Embodiment

Calculate for the large-scale parallel under the Heterogeneous Computing resource, NVIDIA company released GPU general-purpose computations products C UDA (unified calculation equipment framework, Compute Unified Device Architecture) in 2007.This method is take CUDA isomery programming model as the basis, and this programming model of following paper: as Fig. 1, as main frame (HOST), GPU is as coprocessor (co-processer) or equipment (Device) with CPU for this model.Generally there are a main frame and a plurality of equipment in a system.In model, Each performs its own functions for main frame and equipment: main frame is responsible for the strong task of logicality and is processed and serial computing, and highly the parallel task of thread parallelization is given equipment.Carry out computing work with GPU, CPU compares with use, mainly contains several benefits:

Display chip has larger memory bandwidth usually.For example, the GTX580 of NVIDA has and surpasses

1) memory bandwidth of 180GB/s, and the memory bandwidth of the CPU of present high-order is only in the 35GB/s left and right.

2) display chip has more performance element.Have 128 as GeForce 8800GTX " StreamProcessors ", frequency is 1.35GHz.Cpu frequency is usually higher, and is many but the number of performance element will lack.

3) compare with high-end CPU, the price of video card is comparatively cheap.For example a present GTX 460 comprise the price of 1G internal memory and 2.99GHz four core CPU prices similar.

4) GPU-CPU associated treatment framework has become one of classic supercomputing platform, with traditional pure CPU cluster disposal system, the GPU supercomputer cost performance, take up an area the aspects such as space, power consumption advantage very obvious.

The committed step of this method comprises the mathematical morphology filter of the elimination of rough difference of original LiDAR cloud data, the foundation of discrete point cloud data dummy rules graticule mesh, virtual graticule mesh, the attribute that judgement is had a few.Concrete implementation step is as follows:

1. the foundation of discrete point cloud data dummy rules graticule mesh

Set up virtual graticule mesh and do not do the interpolation processing, its purpose is all laser pin points are assigned in each self-corresponding according to planimetric coordinates, and computing formula represents wherein that with reference to formula (1) parameter m of sampling interval multiple is set to 1.With the point that is assigned at first in each piece representative point as this piece.Because each piece records laser pin periods all in this piece, so it has kept original point cloud information fully.

$\left\{\begin{array}{c}I=(X-X_{\min })/\mathrm{mc}\\ J=(Y-Y_{\min })/\mathrm{mc}\\ c=1/\sqrt{n}\end{array}\right.---\left(1\right)$

Wherein, (I, J) is graticule mesh ranks coordinates; (X, Y) is original point cloud planimetric coordinates; X _min, Y _minIt is respectively the minimum value of horizontal ordinate and ordinate in a cloud; C is sampling interval; M is constant, the multiple of expression sampling interval; N is that unit area inner laser pin is counted out.

This step is as follows based on the parallel improvement of GPU:

Discrete points data is copied in the GPU video memory, distribute to different threads and carry out respectively, to reach the purpose of executed in parallel.Namely allow the parallel graticule mesh numbering of finding the solution the place of each discrete point position of all GPU computing units.

2. the elimination of rough difference of original LiDAR cloud data

Adopt the Method of Data Organization of virtual graticule mesh, discrete point cloud is carried out piecemeal process, make comparisons according to certain judgment criterion and the threshold value that calculates for each the laser pin point in piece, will be judged as the rough error point greater than the laser pin point of threshold value.

For laser pin point p, can obtain an elevation difference V _p, computing formula is shown below:

V _p＝M _p-H _p （2）

Wherein, M _pTypical value for laser pin point p; H _pHeight value for laser pin point p.

The elevation difference of having a few all will be used for the counting statistics value, as the basis that determines final threshold value.Arithmetic mean U and standard deviation S D calculate according to formula (3) and formula (4) respectively:

$U=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i/n---\left(3\right)$

$\mathrm{SD}=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(V_i-U)}^2/(n-1)}---\left(4\right)$

Wherein, n is total number of laser pin point.

Threshold value T is traditionally arranged to be 3 times of SD, and for arbitrary laser pin point p, its coarse error criterion is: if | V _p-U|〉T, laser pin point p is judged as the rough error point; Otherwise do not think that p is rough error.

This step is as follows based on the parallel improvement of GPU:

According to the thread enterprise schema of CUDA " grid-block-thread ", the corresponding data of calculating in a graticule mesh of each block, each thread is corresponding to be calculated.By utilizing shared storage to store the cloud data of cloud data corresponding to each block in a certain virtual graticule mesh, guarantee that in same graticule mesh, thread reads virtual grid points cloud data high-speed, realize fine-grained parallel.

3. the mathematical morphology filter of virtual graticule mesh

At first Mathematical Morphology Method is made original discrete point cloud data rule gridding and is processed.Namely choose minimum height value in each grid unit, if there is no laser pin point data in grid unit, get the height value of its neighbor point.

The graticule mesh that obtains after processing is carried out mathematical morphology filter, and method as shown in Figure 2.

(1) use is of a size of l ₁Window carry out opening operation for the first time, dh _p,1Poor for the elevation before and after laser pin point p opening operation, dh _T,1Be the difference of elevation threshold value of opening operation for the first time.If dh _p,1≤ dh _T,1, the p point is judged as ground laser pin point; Otherwise be non-ground laser pin point.dh _{Max (t), 1}For topographical surface maximum elevation before and after opening operation poor, if dh _T,1Dh _{Max (t), 1}, the part of landform shown in figure can be retained

(2) in next iteration, window size increases to l ₂, again carry out mathematical morphology open operator on previous filtering result, and do similar above some determined property.dh _{Max (t), 2}For opening operation front and back topographical surface maximum elevation is poor for the second time, if dh _{Max (t), 2}Less than the difference of elevation threshold value dh of opening operation for the second time _T,2, floor portion branch is retained dh _{Min (b), 2}For the minimum difference of elevation of building object point before and after opening operation, if dh _{Min (b), 2}Dh _T,2, the building object point can effectively be filtered out.

(3) continue to increase window size, so iterate, until window size is greater than maximum buildings.

1) filtering parameter determines

Use Mathematical Morphology class hour, choosing of filtering parameter is extremely important, and it directly affects the quality of filter effect.Filtering parameter is divided into window size and difference of elevation threshold value.

2) window size is determined

A kind of method is to utilize the linear window size that increases of formula (1).

w _k＝2kb+1 （6）

Wherein, k is iterations; w _kBe the k time window size; B is the home window size.The advantage of the method is to keep well the continually varying terrain feature.

Another kind method is that index increases window size, and formula is w _k=2b ^k+ 1 (7)

Its advantage is to reduce iterations.

3) the difference of elevation threshold value is determined

The difference of elevation threshold value can be determined by terrain slope, and the poor dh of maximum elevation of the k time iteration front and back landform _{Max (t), k}With window size w _kThere is following relation with terrain slope s

$s=\frac{d{h}_{\max \left(t\right),k}}{(w_k-w_{k-1})/2}---\left(8\right)$

Therefore, determine difference of elevation threshold value dh _T,kFormula be

$d{h}_{T,k}=\left\{\begin{array}{cc}d{h}_0& w_k\≤3\\ s(w_k-w_{k-1})c+d{h}_0& w_k>3\\ d{h}_{\max }& d{h}_{T,k}>d{h}_{\max }\end{array}\right.---\left(9\right)$

Wherein, dh ₀Be initial difference of elevation threshold value; C is a cloud average headway; dh _maxBe maximum elevation difference limen value.

This step is as follows based on the parallel improvement of GPU:

A. whole virtual Grid square copied to the video memory space and binds texture memory, accelerating changing the access speed of data;

B. for the morphology opening operation, first being tied to texture memory data are carried out erosion operation, operation result is stored in global storage, then result data is tied to texture memory, carries out dilation operation, and acquired results is the result of an opening operation.

C. when carrying out the burn into dilation operation, different threads is distributed in the computing of whole regional difference position carried out respectively, to reach the purpose of executed in parallel.Should suitable thread enterprise schema be set according to the actual hardware environment, guarantee the high usage of computational resource.

D. utilize the fence function " syncthreads " guarantee thread synchronization, namely in judging area, whether all opening operations all are finished.

E. increase window at the GPU end, iterate, until meet the condition that stops iteration.After carrying out above improvement, can carry out simultaneously mathematical morphology filter to reach the purpose of Accelerating running to a large amount of cloud datas.

4. judge the attribute of having a few

Judge according to the difference of elevation of virtual graticule mesh inner laser pin point and its representative point whether this laser pin point has with representative and order the same attribute.So all discrete laser pin points of traversal, complete filtering.

This step is as follows based on the parallel improvement of GPU:

With the original discrete point cloud data that is stored in the GPU video memory, distribute to different threads and carry out respectively, to reach the purpose of executed in parallel.Namely allow the parallel graticule mesh numbering of finding the solution the place of each discrete point position of all GPU computing units.

For whether checking this method can produce a desired effect, this method is tested.

1) experimental data

Fig. 4-1, Fig. 5-1, Fig. 6-1, Fig. 7-1, Fig. 8-1, Fig. 9-1, Figure 10-1, Figure 11-1, Figure 12-1, Figure 13-1, Figure 14-1, Figure 15-1, 16-1, Figure 17-1, Figure 18-1 is the gradual change hypsometric map that the original point cloud of ISPRS15 sample data experiment generates, Fig. 4-2, Fig. 5-2, Fig. 6-2, Fig. 7-2, Fig. 8-2, Fig. 9-2, Figure 10-2, Figure 11-2, Figure 12-2, Figure 13-2, Figure 14-2, Figure 15-2, 16-2, Figure 17-2, in Figure 18-2nd, ISPRS sample experiment based on GPGPU parallel algorithm filtering result, Fig. 4-3, Fig. 5-3, Fig. 6-3, Fig. 7-3, Fig. 8-3, Fig. 9-3, Figure 10-3, Figure 11-3, Figure 12-3, Figure 13-3, Figure 14-3, Figure 15-3, 16-3, Figure 17-3, Figure 18-3rd, reference filtering result in the experiment of ISPRS sample.

2) experimental result and analysis

A. filter effect evaluation

Above-mentioned 1) described in, each figure belongs to qualitative evaluation for the evaluation of this method filter effect, because only compared the gradual change hypsometric map before and after filtering, is the validity of further verification method, and the below will do evaluation to this filtering method with quantitative method.Concrete grammar is the total false rate of each sample after statistical filtering, and compares with filtering method in the ISPRS test report.

Claims (9)

1. mathematical morphology LiDAR point cloud fast filtering method based on GPGPU, it is characterized in that, the method is carried out the strong task of logicality with CPU as main frame and is processed and serial computing, carry out the parallel task of height thread parallelization take GPU as coprocessor and process, the step of filtering method is as follows:

(1) the original LiDAR cloud data of input is in the GPU video memory, and then GPU operates below the LiDAR cloud data is distributed to the different threads executed in parallel: calculate the graticule mesh ranks coordinate of LiDAR cloud data, and data allocations is arrived each self-corresponding graticule mesh;

(2) make each block correspondence in the GPU framework calculate the interior data of a graticule mesh, corresponding each the laser pin point that calculates in each graticule mesh of each thread in block is made comparisons according to certain judgment criterion and threshold value, will be judged as the rough error point greater than the laser pin point of threshold value and reject;

(3) according to the graticule mesh after step (2) excluding gross error point, all laser pin points according to its planimetric coordinates, are assigned in each self-corresponding regular grid, set up thus a dummy rules graticule mesh that can index each laser pin point.

(4) with the data copy in described dummy rules graticule mesh to the video memory space, and bind texture memory, different threads is distributed in the some position of all laser pin points in described dummy rules graticule mesh uses respectively and be of a size of l ₁Window each data that are tied to texture memory are carried out mathematical morphology open operator, and judge whether data corresponding to operation result are ground laser pin points, and filtering non-ground laser pin point;

(5) window size and the discrepancy in elevation threshold parameter that import into according to the CPU end increase described window size at the GPU end, and the operation of repeating step (4) is until window size is greater than maximum buildings;

(6) according to the difference of elevation through the laser pin point in the dummy rules graticule mesh after the laser pin point in the graticule mesh after step (2) and process step (4), whether the laser pin point in determining step (2) graticule mesh has and orders the same attribute through the dummy rules graticule mesh laser pin after step (4), so travel through all discrete laser pin points, remove the non-ground laser pin point in graticule mesh in step (2), complete filtering.

2. method according to claim 1 is characterized in that: in described step (1), the minimum height value of each grid unit as representative point, if there is no laser pin point data in grid unit, is got the height value of its neighbor point.

3. method according to claim 1, is characterized in that, in described step (1), the computing formula of graticule mesh ranks coordinate is $\left\{\begin{array}{c}I=(X-X_{\min })/\mathrm{mc}\\ J=(Y-Y_{\min })/\mathrm{mc}\\ c=1/\sqrt{n}\end{array}\right.,$ Wherein, (I, J) is graticule mesh ranks coordinates; (X, Y) is original point cloud planimetric coordinates; X _min, Y _minIt is respectively the minimum value of horizontal ordinate and ordinate in a cloud; C is sampling interval; M is constant, the multiple of expression sampling interval; N is that unit area inner laser pin is counted out.

4. method according to claim 1, is characterized in that, in described step (2), coarse error criterion is: if | V _p-U|〉T, laser pin point p is judged as the rough error point; Otherwise do not think that p is rough error, wherein the elevation difference V of laser pin point p _p=M _p-H _p, M _pTypical value for laser pin point p; H _pHeight value for laser pin point p; Standard deviation

Arithmetic mean

N is total number of laser pin point, and T is threshold value.

5. method according to claim 4, it is characterized in that: described T is set to 3 times of SD.

6. method according to claim 1 is characterized in that: in step (4), the criterion of judgement ground laser pin point is: establish dh _p,1Poor for the elevation before and after laser pin point p opening operation, dh _T,1Be the difference of elevation threshold value of opening operation for the first time, if dh _p,1≤ dh _T,1, the p point is judged as ground laser pin point; Otherwise be non-ground laser pin point; dh _{Max (t), 1}For topographical surface maximum elevation before and after opening operation poor, if dh _T,1Dh _{Max (t), 1}, the part of ground laser pin point can be retained, but not topocentric part can be removed.

7. method according to claim 6, is characterized in that: described definite difference of elevation threshold value dh _T,kPublic dh ₀Be initial difference of elevation threshold value; C is a cloud average headway; dh _maxBe maximum elevation difference limen value, the poor dh of maximum elevation of landform before and after the k time iteration _{Max (t), k}With window size w _kWith the pass of terrain slope s be $s=\frac{d{h}_{\max \left(t\right),k}}{(w_k-w_{k-1})/2}.$

8. the described method of any one according to claim 1-7, is characterized in that, in described step (5), the increase mode of window size is w _k=2kb+1 or w _k=2b ^k+ 1, wherein, k is iterations; w _kBe the k time window size; B is the home window size.

9. method according to claim 4, is characterized in that, in described dummy rules graticule mesh, each graticule mesh has and only have a laser pin point.

Images