CN110796741B - Airborne laser sounding point cloud filtering method based on bidirectional cloth simulation - Google Patents
Airborne laser sounding point cloud filtering method based on bidirectional cloth simulation Download PDFInfo
- Publication number
- CN110796741B CN110796741B CN201910902029.1A CN201910902029A CN110796741B CN 110796741 B CN110796741 B CN 110796741B CN 201910902029 A CN201910902029 A CN 201910902029A CN 110796741 B CN110796741 B CN 110796741B
- Authority
- CN
- China
- Prior art keywords
- point
- filtering
- laser
- node
- submarine
- 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.)
- Expired - Fee Related
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 101
- 239000004744 fabric Substances 0.000 title claims abstract description 59
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000004088 simulation Methods 0.000 title claims abstract description 35
- 230000002457 bidirectional effect Effects 0.000 title claims abstract description 30
- 238000012876 topography Methods 0.000 claims abstract description 27
- 230000002159 abnormal effect Effects 0.000 claims abstract description 18
- 238000012937 correction Methods 0.000 claims abstract description 16
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 20
- 230000033001 locomotion Effects 0.000 claims description 10
- 230000001133 acceleration Effects 0.000 claims description 9
- 238000010276 construction Methods 0.000 claims description 9
- 238000001514 detection method Methods 0.000 claims description 8
- 230000000694 effects Effects 0.000 claims description 4
- 238000012546 transfer Methods 0.000 claims description 4
- 230000005489 elastic deformation Effects 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 239000000463 material Substances 0.000 claims 1
- 238000013507 mapping Methods 0.000 abstract description 3
- 230000005540 biological transmission Effects 0.000 abstract description 2
- 238000005259 measurement Methods 0.000 description 3
- 241000251468 Actinopterygii Species 0.000 description 1
- 241000195493 Cryptophyta Species 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000007667 floating Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009401 outcrossing Methods 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000013535 sea water Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T19/00—Manipulating 3D models or images for computer graphics
- G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2200/00—Indexing scheme for image data processing or generation, in general
- G06T2200/04—Indexing scheme for image data processing or generation, in general involving 3D image data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Remote Sensing (AREA)
- Architecture (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Image Processing (AREA)
Abstract
The invention discloses an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation, which belongs to the technical field of ocean mapping, and comprises the steps of firstly acquiring airborne laser sounding data and calculating the coordinates of a submarine laser point; then eliminating the negative abnormal data of the seabed by constructing a transmission type iteration trend surface and expressing the continuity of the seabed topography; and finally, constructing a filtering surface according to the bidirectional cloth simulation correction model to finish airborne laser sounding point cloud filtering. The method realizes high-precision and strong-robustness filtering of the airborne laser sounding point cloud, and effectively solves the problems that the existing method can not identify negative abnormal data of the airborne laser sounding and is easy to generate excessive filtering.
Description
Technical Field
The invention belongs to the technical field of ocean mapping, and particularly relates to an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation.
Background
The airborne laser sounding system has the characteristics of high measurement accuracy, high measurement point density, high working efficiency, strong maneuverability, measurement continuity and the like, is particularly suitable for rapid detection of complex terrains such as shallow water areas, island reef nearby areas and the like, and can realize seamless splicing of coastline water and underwater terrains. Has important significance for meeting urgent demands of ocean, mapping, water conservancy, traffic, outcrossing, navy and other aspects. The airborne laser sounding system adopts 532nm blue-green laser with stronger water body penetrating ability as a submarine detection means, when the blue-green laser reaches the sea surface, one part of laser beams are reflected along the original path, and the other part of laser beams penetrate the sea surface to enter the water body. However, when a large amount of floating algae, fish shoals or submarine vegetation exist in the shallow sea water body, laser beams entering the water body cannot reach the sea bottom and are reflected in advance, so that a large amount of irregular noise points are generated, and the authenticity expression of the submarine topography is affected. The submarine topography is a direct basis for understanding the external shape of the earth, the movement of submarine structures and the submarine evolution, and is also important basic data in the aspects of ocean economic development, ocean scientific research, ocean military application and the like. Therefore, in order to extract the real submarine topography, it is necessary to study an effective submarine laser point cloud filtering method.
The current mainstream filtering methods are mainly divided into multi-beam water depth data filtering and land airborne laser scanning data filtering, but the applicability of the methods to airborne laser sounding data is not verified, and the methods related to airborne laser sounding point cloud data filtering are relatively few, so that an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation is necessary to be provided, and high-precision and strong-robustness filtering of airborne laser sounding point clouds is realized.
Disclosure of Invention
Aiming at the technical problems in the prior art, the invention provides the airborne laser sounding point cloud filtering method based on bidirectional distribution simulation, which is reasonable in design, solves the problems that the airborne laser sounding negative abnormal data cannot be identified and excessive filtering is easy to generate in the prior art, fully reserves the characteristics of seafloor micro-topography, ensures the integrity of seafloor topography data, and has better applicability and robustness for the airborne laser sounding data.
In order to achieve the above purpose, the present invention adopts the following technical scheme:
an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation comprises the following steps:
step 1: acquiring airborne laser sounding point cloud data, and calculating the coordinates of the submarine laser points;
step 2: fitting by adopting a quadric LM (Levenberg-Marquardt) algorithm, constructing a transfer iteration trend surface, removing submarine negative abnormal data, and expressing submarine continuity;
step 3: and constructing a bidirectional distribution simulation correction model, and performing bidirectional filtering of the submarine laser point cloud according to the correction model to obtain a final filtering surface.
Preferably, in step 2, the construction method of the transitive iterative trend surface is as follows:
step 2.1: quadric LM algorithm fitting: assuming that a point P is taken as a center, r is taken as a radius, searching all submarine laser points in a given radius area, and a quadric surface model formed by fitting all laser points in the area is shown as a formula (1):
wherein, (x) j ,y j ,z j ) For the point set N j Geographic coordinates of the laser point within (j=1, 2, …, k) in the local coordinate system; a. b, c, d, e and f are quadric fitting parameters;
step 2.2: at point P 1 Taking r as a radius as a center, searching all laser points in the area and performing quadric surface LM fitting; calculating the water depth average value of laser points in curved surfaces
Standard deviation sigma of water depth H Find and +.>
The difference is greater than 2 sigma H The laser points are marked as negative abnormal points and stored into a suspicious point set, and the water depth values of all suspicious points are calculated into a fitted quadric surface to ensure a discrete point set N 1 The method does not contain coarse differences and ensures the integrity of sounding data;
step 2.3: find and point P 1 The nearest point, designated as point P 2 Find (N) 2 -N 1 ∩N 2 ) Point concentration
The difference is greater than 2 sigma H The laser points are marked as negative abnormal points and stored into a suspicious point set, and the water depth values of all suspicious points are calculated into a fitted quadric surface to ensure a discrete point set N 2 The method does not contain coarse differences;
step 2.4: find and point P 2 The nearest point, excluding the completed point, is denoted as point P 3 And (2.3) repeating the step until all submarine laser points are traversed, storing no laser points into a suspicious point set, and obtaining the final submarine topography trend surface when the fitted topography trend surface reaches a stable state.
Preferably, in step 3, the specific construction method of the bidirectional cloth simulation correction model is as follows:
step 3.1: setting filtering parameters: the grid resolution GR, namely the distance between two adjacent nodes; the stiffness parameter IMRI controls the softness of the simulated cloth; the maximum iteration number MI is used for judging whether the simulation process is terminated; judging classification thresholds of the sea bottom points and the non-sea bottom points according to the self-adaptive distance threshold H;
step 3.2: respectively performing forward filtering and backward filtering on the submarine laser points, namely respectively inverting the submarine laser points and not inverting the submarine laser points, and placing the simulated cloth right above the laser points and keeping a certain distance from the laser points; projecting the distribution node and the laser point to the same horizontal plane, performing grid-networking, determining the corresponding relation between the distribution node and the laser point, and recording the elevation values of the distribution node and the laser point;
step 3.3: the simulated cloth is enabled to do free falling motion, and the position and the stress effect of the cloth node meet the formula (2):
wherein m represents node Z l Is the mass of (3); x represents Z l Position at time t; g represents gravitational acceleration; c e Representing the elastic deformation coefficient of the node spring; e represents Z l A set of neighbors;
representing node Z l And node Z n Distance at time t>
Representing node Z l And node Z n Original distance between>
Representing Z l Pointing to Z n Is a unit vector of (a);
step 3.4: in the process of simulating the descending of cloth, performing collision detection on cloth nodes;
assuming that the collision state is defined as C, the collision constraint is determined by the position and velocity constraints of the cloth node, and there are two time steps:
where dt is the time step; XC (t), XC (t+dt) and XC (t+2dt) are the positions of the collision at times t, t+dt and t+2dt, respectively; XC ' (t), XC ' (t+dt) and XC ' (t+2dt) are the movement velocities of the collisions at times t, t+dt and t+2dt, respectively; XC "(t) and XC" (t+dt) are the motion accelerations of the collision at times t and t+dt, respectively; Δxc "(t) and Δxc" (t+dt) are corrected accelerations of the collision at times t and t+dt, respectively, and thus equation (4) can be obtained:
wherein Δxc "(t) can be converted into additional impact forces acting on the distribution node;
step 3.5: judging the position relation between the distribution node and the submarine laser points on the basis of collision detection;
if the elevation value of the cloth node is H cp Less than the elevation value H of the submarine laser spot sp Moving the node to the corresponding laser point location and marking the node as an immovable point; otherwise calculate the node inside the nodeDisplacement to be adjusted under the action of traction force;
step 3.6: repeating the steps 3.3 to 3.5 until the iteration times of the nodes are larger than the set maximum iteration times MI, and ending the cloth simulation process;
step 3.7: calculating the height difference delta H between the distribution node and the corresponding submarine laser point, and judging the size relation between delta H and the self-adaptive distance threshold H; if delta H is less than H, the node is reserved as a submarine point, otherwise, the node is judged to be a non-submarine point and deleted;
step 3.8: establishing a forward filtering surface F F (P T IMRI, H) and inverse filtering plane F I (P T ,IMRI,H)；
Defining a forward filtering surface as the simulated cloth formed after forward filtering, and defining a reverse filtering surface as the simulated cloth formed after reverse filtering; let the forward filtering surface be F F (P T IMRI, H), on which the distribution node is denoted P F (x F ,y F ,z F ) The method comprises the steps of carrying out a first treatment on the surface of the The reverse filtering surface is F I (P T IMRI, H), on which the distribution node is denoted P I (x I ,y I ,z I ) The method comprises the steps of carrying out a first treatment on the surface of the The two filter surfaces are represented as shown in equation (5):
wherein F is F (P T IMRI, H) and F I (P T IMRI, H) is the forward and reverse filtering face; p (P) T (x T ,y T ,z T )(z T ＝f(x T ,y T ) A submarine laser spot); the IMRI is a rigidity parameter for controlling the softness of the cloth; h is a distance threshold value for judging that the distribution node is a submarine point and a non-submarine point; (x) F ,y F ,z F ) And (x) I ,y I ,z I ) For distributing the node P F And P I Is satisfied with z F ＝f F (x F ,y F ) And z I ＝f I (x I ,y I ) The method comprises the steps of carrying out a first treatment on the surface of the T represents the sea floor topography;
step 3.9:establishing a submarine filtering surface F (P) T ,IMRI,H)；
Since some nodes of the forward and reverse filtering surfaces will be moved to the same subsea laser spot position, the subsea filtering surface F (P T IMRI, H) can be expressed as forward filtering plane F F (P T IMRI, H) and inverse filtering plane F I (P T IMRI, H) as shown in equation (6):
F(P T ,IMRI,H)＝F F (P T ,IMRI,H)∪F I (P T ,IMRI,H) (6)。
the invention has the beneficial technical effects that:
compared with the prior art, the invention utilizes a quadric surface LM algorithm to construct a transfer type iteration trend surface, eliminates negative abnormal data of the seabed, and constructs a bidirectional cloth simulation filtering model by combining a cloth simulation technology and bidirectional filtering so as to obtain a final seabed filtering surface, thereby solving the problem that excessive filtering is easy to generate in the prior art, effectively improving the filtering precision of the airborne laser sounding point cloud, and enhancing the robustness of sounding point cloud filtering.
Drawings
Fig. 1 is a flow chart of an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation.
FIG. 2 is a schematic diagram of a transitive iterative trend surface construction in the present invention.
FIG. 3 is a flow chart of the construction of the simulation correction model of the bidirectional cloth in the invention.
Detailed Description
The invention is described in further detail below with reference to the attached drawings and detailed description:
the invention provides an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation, and the flow of the filtering method is shown in figure 1.
The method comprises the following steps:
step 1: and acquiring airborne laser sounding data, and calculating the coordinates of the submarine laser points.
Acquiring sounding laser data by using an airborne LiDAR sounding system; the effective sounding laser point is the submarine laser point after attitude correction, multi-source data fusion and refraction correction.
Step 2: and (5) quadric surface fitting, and constructing a transitive iteration trend surface.
The natural topography of the sea floor is continuously variable, not abrupt, despite the gradient, as the sea dynamics erode the sea floor irregularities. Due to the fact that secondary reflection and the like of the submarine echo occur, negative abnormal points lower than the real submarine topography inevitably occur in actual submarine topography, however, in a cloth simulation filtering algorithm, filtering processing cannot be performed when negative abnormality exists. Therefore, it is necessary to establish a submarine trend surface, so that negative abnormal data can be removed, and continuity of submarine topography can be expressed. The construction of the transitive iterative trend surface is shown in fig. 2.
In a further embodiment, the step 2 specifically includes the following steps:
step 2.1: quadric LM algorithm fitting: assuming that a point P is taken as a center, r is taken as a radius, searching all submarine laser points in a given radius area, and fitting all laser points in the area to form a quadric surface model as follows:
wherein, (x) j ,y j ,z j ) For the point set N j Geographic coordinates of the laser point within (j=1, 2, …, k) in the local coordinate system; a. b, c, d, e and f are quadric fitting parameters, and fitting parameter values can be obtained based on an LM algorithm;
step 2.2: at point P 1 Taking r as a radius for centering, searching all laser points in the area and performing quadric LM fitting. Calculating the water depth average value of laser points in curved surfaces
Standard deviation sigma of water depth H Find and +.>
The difference is greater than 2 sigma H The laser points are marked as negative abnormal points and stored into a suspicious point set, and the water depth values of all suspicious points are calculated into a fitted quadric surface to ensure a discrete point set N 1 The method does not contain coarse differences and ensures the integrity of sounding data;
step 2.3: find and point P 1 The nearest point, designated as point P 2 Due to (N) 1 ∩N 2 ) The dot concentration of (2) no longer contains coarse differences, and is defined by N 2 The quadric surface fitting accuracy is higher, and the (N 2 -N 1 ∩N 2 ) Point concentration
The difference is greater than 2 sigma H The laser points are marked as negative abnormal points and stored into a suspicious point set, and the water depth values of all suspicious points are calculated into a fitted quadric surface to ensure a discrete point set N 2 The method does not contain coarse differences;
step 2.4: find and point P 2 The nearest point (excluding the completed point) is denoted as point P 3 And (2.3) repeating the step until all submarine laser points are traversed, storing no laser points into a suspicious point set, and obtaining the final submarine topography trend surface when the fitted topography trend surface reaches a stable state.
In the concrete implementation, the construction of the transmission type iteration trend surface not only effectively eliminates the negative abnormal points on the sea bottom, but also reduces the water depth values of the negative abnormal points to the fitted trend surface, thereby ensuring the integrity of sounding data and providing guarantee for smooth implementation of bidirectional cloth simulation filtering.
Step 3: and constructing a bidirectional cloth simulation correction model to obtain the submarine filtering surface.
The submarine topography is complex and changeable, when the submarine topography has more obvious positive topography (such as cliffs, submarine volcanoes, mud volcanic and reef) or negative topography (such as notching, notching and pit), the submarine laser point cloud filtering can generate more serious excessive filtering phenomenon by adopting a forward distribution simulation method, and the integrity of the submarine topography is destroyed. Therefore, the bidirectional cloth simulation correction model is adopted to carry out the filtering of the seabed laser point cloud, and the problem of excessive filtering of the point cloud is effectively avoided. The specific flow of the construction of the bidirectional cloth simulation correction model is shown in fig. 3.
In a further embodiment, the step 3 specifically includes the following steps:
step 3.1: setting filtering parameters: the grid resolution GR, namely the distance between two adjacent nodes; the stiffness parameter IMRI controls the softness of the simulated cloth; the maximum iteration number MI is used for judging whether the simulation process is terminated; judging classification thresholds of the sea bottom points and the non-sea bottom points according to the self-adaptive distance threshold H;
step 3.2: respectively performing forward filtering and backward filtering on the submarine laser points, namely respectively inverting the submarine laser points and not inverting the submarine laser points, and placing the simulated cloth right above the laser points and keeping a certain distance from the laser points; projecting the distribution node and the laser point to the same horizontal plane, performing grid-networking, determining the corresponding relation between the distribution node and the laser point, and recording the elevation values of the distribution node and the laser point;
step 3.3: the simulated cloth is enabled to do free falling movement, and the position and the stress effect of the cloth node meet the formula (2);
wherein m represents node Z l Is the mass of (3); x represents Z l Position at time t; g represents gravitational acceleration; c e Representing the elastic deformation coefficient of the node spring; e represents Z l A set of neighbors;
representing node Z l And node Z n Distance at time t>
Representing node Z l And node Z n Original distance between>
Representation ofZ l Pointing to Z n Is a unit vector of (a).
Step 3.4: in the process of simulating the descending of cloth, performing collision detection on cloth nodes;
assuming that the collision state is defined as C, the collision constraint is determined by the position and velocity constraints of the cloth node, and there are two time steps:
where dt is the time step; XC (t), XC (t+dt) and XC (t+2dt) are the positions of the collision at times t, t+dt and t+2dt, respectively; XC ' (t), XC ' (t+dt) and XC ' (t+2dt) are the movement velocities of the collisions at times t, t+dt and t+2dt, respectively; XC "(t) and XC" (t+dt) are the motion accelerations of the collision at times t and t+dt, respectively; Δxc "(t) and Δxc" (t+dt) are corrected accelerations of the collision at times t and t+dt, respectively, and thus equation (4) can be obtained:
wherein Δxc "(t) can be converted into additional impact forces acting on the distribution node;
step 3.5: based on collision detection, judging the position relationship between the distribution node and the submarine laser points:
if the height value H of the cloth node is satisfied cp Less than the elevation value H of the submarine laser spot sp Moving the node to the corresponding laser point location and marking the node as an immovable point; otherwise, calculating the displacement of the node which needs to be adjusted under the action of internal force (pulling force inside the node).
Step 3.6: repeating the steps 3.3-3.5 until the iteration times of the nodes are larger than the set maximum iteration times MI, and ending the cloth simulation process.
Step 3.7: calculating the height difference delta H between the distribution node and the corresponding submarine laser point, and judging the size relation between delta H and the self-adaptive distance threshold H; if ΔH < H, the node is reserved as a submarine point, otherwise, the node is judged as a non-submarine point and deleted.
Step 3.8: establishing a forward filtering surface F F (P T IMRI, H) and inverse filtering plane F I (P T ,IMRI,H)：
Defining a forward filtering surface as the simulated cloth formed after forward filtering, and defining a reverse filtering surface as the simulated cloth formed after reverse filtering; let the forward filtering surface be F F (P T IMRI, H), on which the distribution node is denoted P F (x F ,y F ,z F ) The method comprises the steps of carrying out a first treatment on the surface of the The reverse filtering surface is F I (P T IMRI, H), on which the distribution node is denoted P I (x I ,y I ,z I ) The method comprises the steps of carrying out a first treatment on the surface of the The two filtering surfaces can be expressed as
Wherein F is F (P T IMRI, H) and F I (P T IMRI, H) is the forward and reverse filtering face; p (P) T (x T ,y T ,z T )(z T ＝f(x T ,y T ) A submarine laser spot); the IMRI is a rigidity parameter for controlling the softness of the cloth; h is a distance threshold value for judging that the distribution node is a submarine point and a non-submarine point; (x) F ,y F ,z F ) And (x) I ,y I ,z I ) For distributing the node P F And P I Is satisfied with z F ＝f F (x F ,y F ) And z I ＝f I (x I ,y I ) The method comprises the steps of carrying out a first treatment on the surface of the T represents the sea floor topography;
step 3.9: establishing a submarine filtering surface F (P) T ,IMRI,H)：
Since some nodes of the forward and reverse filtering surfaces will be moved to the same subsea laser spot position, the subsea filtering surface F (P T IMRI, H) can be expressed as forward filtering plane F F (P T IMRI, H) and inverse filtering plane F I (P T Union of IMRI, H), i.e
F(P T ,IMRI,H)＝F F (P T ,IMRI,H)∪F I (P T ,IMRI,H) (6)。
In the concrete implementation, the submarine filtering surface is constructed according to the bidirectional cloth simulation correction model, so that non-submarine points are effectively filtered, the submarine micro-topography characteristics are fully reserved, and the integrity of submarine topography data is ensured. Therefore, the submarine filtering surface can accurately express the actual submarine topography, and the bidirectional cloth simulation filtering method has better applicability and robustness for airborne laser sounding data.
In summary, the invention provides an airborne laser sounding point cloud filtering method based on bidirectional cloth simulation, which comprises the following steps: on the premise of obtaining the submarine laser point coordinates, adopting a quadric LM algorithm to fit and construct a transfer type iteration trend surface, removing submarine negative abnormal data, and expressing the submarine continuity; and constructing a bidirectional distribution simulation correction model, and performing bidirectional filtering on the submarine laser point cloud according to the correction model to obtain a final submarine filtering surface. The method effectively solves the problems that the existing method can not identify negative abnormal data of airborne laser sounding and is easy to generate excessive filtering, and realizes high-precision and strong-robustness filtering of the airborne laser sounding point cloud.
It should be understood that the above description is not intended to limit the invention to the particular embodiments disclosed, but to limit the invention to the particular embodiments disclosed, and that the invention is not limited to the particular embodiments disclosed, but is intended to cover modifications, adaptations, additions and alternatives falling within the spirit and scope of the invention.
Claims (2)
1. An airborne laser sounding point cloud filtering method based on bidirectional cloth simulation is characterized in that: the method comprises the following steps:
step 1: acquiring airborne laser sounding point cloud data, and calculating the coordinates of the submarine laser points;
step 2: adopting quadric LM algorithm fitting, constructing a transfer type iteration trend surface, removing negative abnormal data of the seabed, and expressing the continuity of the seabed;
step 3: constructing a bidirectional distribution simulation correction model, and performing bidirectional filtering of the submarine laser point cloud according to the correction model to obtain a final filtering surface;
the specific construction method of the bidirectional cloth simulation correction model is as follows:
step 3.1: setting filtering parameters: the grid resolution GR, namely the distance between two adjacent nodes; the stiffness parameter IMRI controls the softness of the simulated cloth; the maximum iteration number MI is used for judging whether the simulation process is terminated; judging classification thresholds of the sea bottom points and the non-sea bottom points according to the self-adaptive distance threshold H;
step 3.2: respectively performing forward filtering and backward filtering on the submarine laser points, namely respectively inverting the submarine laser points and not inverting the submarine laser points, and placing the simulated cloth right above the laser points and keeping a certain distance from the laser points; projecting the distribution node and the laser point to the same horizontal plane, performing grid-networking, determining the corresponding relation between the distribution node and the laser point, and recording the elevation values of the distribution node and the laser point;
step 3.3: the simulated cloth is enabled to do free falling motion, and the position and the stress effect of the cloth node meet the formula (2):
wherein m represents node Z l Is the mass of (3); x represents Z l Position at time t; g represents gravitational acceleration; c e Representing the elastic deformation coefficient of the node spring; e represents Z l A set of neighbors;
representing node Z l And node Z n Distance at time t>
Representing node Z l And node Z n Original distance between>
Representing Z l Pointing to Z n Is a unit vector of (a);
step 3.4: in the process of simulating the descending of cloth, performing collision detection on cloth nodes;
assuming that the collision state is defined as C, the collision constraint is determined by the position and velocity constraints of the cloth node, and there are two time steps:
where dt is the time step; XC (t), XC (t+dt) and XC (t+2dt) are the positions of the collision at times t, t+dt and t+2dt, respectively; XC ' (t), XC ' (t+dt) and XC ' (t+2dt) are the movement velocities of the collisions at times t, t+dt and t+2dt, respectively; XC "(t) and XC" (t+dt) are the motion accelerations of the collision at times t and t+dt, respectively; Δxc "(t) and Δxc" (t+dt) are corrected accelerations of the collision at times t and t+dt, respectively, and thus equation (4) can be obtained:
wherein Δxc "(t) can be converted into additional impact forces acting on the distribution node;
step 3.5: judging the position relation between the distribution node and the submarine laser points on the basis of collision detection;
if the elevation value H of the material distribution node cp Elevation value H smaller than the undersea laser spot sp Moving the node to the corresponding laser point location and marking the node as an immovable point; otherwise, calculating the displacement of the node to be adjusted under the action of the traction force in the node;
step 3.6: repeating the steps 3.3 to 3.5 until the iteration times of the nodes are larger than the set maximum iteration times MI, and ending the cloth simulation process;
step 3.7: calculating the height difference delta H between the distribution node and the corresponding submarine laser point, and judging the size relation between delta H and the self-adaptive distance threshold H; if delta H is less than H, the node is reserved as a submarine point, otherwise, the node is judged to be a non-submarine point and deleted;
step 3.8: establishing a forward filtering surface F F (P T IMRI, H) and inverse filtering plane F I (P T ,IMRI,H)；
Defining a forward filtering surface as the simulated cloth formed after forward filtering, and defining a reverse filtering surface as the simulated cloth formed after reverse filtering; let the forward filtering surface be F F (P T IMRI, H), on which the distribution node is denoted P F (x F ,y F ,z F ) The method comprises the steps of carrying out a first treatment on the surface of the The reverse filtering surface is F I (P T IMRI, H), on which the distribution node is denoted P I (x I ,y I ,z I ) The method comprises the steps of carrying out a first treatment on the surface of the The two filter surfaces are represented as shown in equation (5):
wherein F is F (P T IMRI, H) and F I (P T IMRI, H) is the forward and reverse filtering face; p (P) T (x T ,y T ,z T ) Z is a submarine laser spot T ＝f(x T ,y T ) The method comprises the steps of carrying out a first treatment on the surface of the The IMRI is a rigidity parameter for controlling the softness of the cloth; h is a distance threshold value for judging that the distribution node is a submarine point and a non-submarine point; (x) F ,y F ,z F ) And (x) I ,y I ,z I ) For distributing the node P F And P I Is satisfied with z F ＝f F (x F ,y F ) And z I ＝f I (x I ,y I ) The method comprises the steps of carrying out a first treatment on the surface of the T represents the sea floor topography;
step 3.9: establishing a submarine filtering surface F (P) T ,IMRI,H)；
Since some nodes of the forward and reverse filtering surfaces will be moved to the same subsea laser spot position, the subsea filtering surface F (P T IMRI, H) can be expressed as forward filtering plane F F (P T IMRI, H) and inverse filtering plane F I (P T IMRI, H) as shown in equation (6):
F(P T ,IMRI,H)＝F F (P T ,IMRI,H)∪F I (P T ,IMRI,H) (6)。
2. the airborne laser sounding point cloud filtering method based on bidirectional cloth simulation according to claim 1, wherein the method is characterized in that: in step 2, the construction method of the transitive iterative trend surface is as follows:
step 2.1: quadric LM algorithm fitting: assuming that a point P is taken as a center, r is taken as a radius, searching all submarine laser points in a given radius area, and a quadric surface model formed by fitting all laser points in the area is shown as a formula (1):
wherein, (x) j ,y j ,z j ) For the point set N j Geographic coordinates of the laser points in the local coordinate system, j=1, 2, …, k; a. b, c, d, e and f are quadric fitting parameters;
step 2.2: at point P 1 Taking r as a radius as a center, searching all laser points in the area and performing quadric surface LM fitting; calculating the water depth average value of laser points in curved surfaces
Standard deviation sigma of water depth H Find and +.>
The difference is greater than 2 sigma H The laser points are marked as negative abnormal points and stored into a suspicious point set, and the water depth values of all suspicious points are calculated into a fitted quadric surface to ensure a discrete point set N 1 The method does not contain coarse differences and ensures the integrity of sounding data;
step 2.3: find and point P 1 The nearest point, designated as point P 2 Find (N) 2 -N 1 ∩N 2 ) Point concentration
The difference is greater than 2 sigma H The laser points are marked as negative abnormal points and stored into a suspicious point set, and the water depth values of all suspicious points are calculated into a fitted quadric surface to ensure a discrete point set N 2 The method does not contain coarse differences;
step 2.4: find and point P 2 The nearest point, excluding the completed point, is denoted as point P 3 And (2.3) repeating the step until all submarine laser points are traversed, storing no laser points into a suspicious point set, and obtaining the final submarine topography trend surface when the fitted topography trend surface reaches a stable state.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910902029.1A CN110796741B (en) | 2019-09-24 | 2019-09-24 | Airborne laser sounding point cloud filtering method based on bidirectional cloth simulation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910902029.1A CN110796741B (en) | 2019-09-24 | 2019-09-24 | Airborne laser sounding point cloud filtering method based on bidirectional cloth simulation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110796741A CN110796741A (en) | 2020-02-14 |
| CN110796741B true CN110796741B (en) | 2023-04-25 |
Family
ID=69438592
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910902029.1A Expired - Fee Related CN110796741B (en) | 2019-09-24 | 2019-09-24 | Airborne laser sounding point cloud filtering method based on bidirectional cloth simulation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110796741B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111369436A (en) * | 2020-02-27 | 2020-07-03 | 山东科技大学 | Airborne LiDAR point cloud rarefying method considering multi-terrain features |
| CN113156393B (en) * | 2021-03-29 | 2022-05-27 | 山东科技大学 | Airborne laser sounding broken wind wave sea surface model construction method |
| CN114820400B (en) * | 2022-07-01 | 2022-09-23 | 湖南盛鼎科技发展有限责任公司 | Airborne LiDAR point cloud ground point filtering method |
| CN115880189B (en) * | 2023-02-22 | 2023-05-30 | 山东科技大学 | Multi-beam point cloud filtering method for submarine topography |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015039375A1 (en) * | 2013-09-17 | 2015-03-26 | 中国科学院深圳先进技术研究院 | Method and system for automatically optimizing quality of point cloud data |
| CN107798657A (en) * | 2017-10-30 | 2018-03-13 | 武汉珞珈新空科技有限公司 | A kind of vehicle-mounted laser point cloud filtering method based on circular cylindrical coordinate |
| CN109299739A (en) * | 2018-09-26 | 2019-02-01 | 速度时空信息科技股份有限公司 | The method that vehicle-mounted laser point cloud is filtered based on the surface fitting of normal vector |
| CN110047036A (en) * | 2019-04-22 | 2019-07-23 | 重庆交通大学 | Territorial laser scanning data building facade extracting method based on polar coordinates grid |
| CN110119438A (en) * | 2019-04-23 | 2019-08-13 | 东华理工大学 | Airborne LiDAR point cloud filtering method based on Active Learning |
-
2019
- 2019-09-24 CN CN201910902029.1A patent/CN110796741B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015039375A1 (en) * | 2013-09-17 | 2015-03-26 | 中国科学院深圳先进技术研究院 | Method and system for automatically optimizing quality of point cloud data |
| CN107798657A (en) * | 2017-10-30 | 2018-03-13 | 武汉珞珈新空科技有限公司 | A kind of vehicle-mounted laser point cloud filtering method based on circular cylindrical coordinate |
| CN109299739A (en) * | 2018-09-26 | 2019-02-01 | 速度时空信息科技股份有限公司 | The method that vehicle-mounted laser point cloud is filtered based on the surface fitting of normal vector |
| CN110047036A (en) * | 2019-04-22 | 2019-07-23 | 重庆交通大学 | Territorial laser scanning data building facade extracting method based on polar coordinates grid |
| CN110119438A (en) * | 2019-04-23 | 2019-08-13 | 东华理工大学 | Airborne LiDAR point cloud filtering method based on Active Learning |
Non-Patent Citations (1)
| Title |
|---|
| 基于激光点云约束的多波束低掠射波束几何改正;卜宪海;《CNKI硕士电子期刊》;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110796741A (en) | 2020-02-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110796741B (en) | Airborne laser sounding point cloud filtering method based on bidirectional cloth simulation | |
| CN102446367B (en) | Method for constructing three-dimensional terrain vector model based on multi-beam sonar submarine measurement data | |
| CN107314768B (en) | Underwater terrain matching auxiliary inertial navigation positioning method and positioning system thereof | |
| CN110906935B (en) | Unmanned ship path planning method | |
| CN103292792B (en) | Actual measurement SVP reconstruction method suitable for submarine detection and pseudo-landform processing | |
| CN108489491A (en) | A kind of Three-dimensional Track Intelligent planning method of autonomous underwater vehicle | |
| CN109752727B (en) | Airborne LiDAR depth sounding sea air interface refraction correction method | |
| CN102419436A (en) | Multibeam data processing method based on total propagation error filter | |
| CN115167465B (en) | Unmanned submersible vehicle three-dimensional path planning method based on artificial potential field grid method | |
| CN111665846A (en) | Water surface unmanned ship path planning method based on rapid scanning method | |
| CN115145315A (en) | Unmanned aerial vehicle path planning method suitable for chaotic environment and with improved A-star algorithm | |
| CN111523077A (en) | Correction calculation method for terrain in middle-far area in coastal zone land gravity measurement | |
| CN113156393B (en) | Airborne laser sounding broken wind wave sea surface model construction method | |
| KR101271402B1 (en) | Interpolation method of erosion-based fractal river channel and computer readable media using the same | |
| CN104680583B (en) | A kind of method that sea-floor relief automatically generates | |
| CN107554719B (en) | A kind of ship load measurement method based on Sonar system | |
| CN105388460A (en) | Indoor underwater target positioning method based on genetic algorithm | |
| CN115031733A (en) | Multi-beam sounding active synchronous positioning and mapping method for underwater robot | |
| CN106908036B (en) | A kind of AUV multi-beam Bathymetric Data patterning process based on local offset | |
| CN113239520B (en) | Near-water three-dimensional underwater topography environment modeling method | |
| CN113012286B (en) | Method for constructing road DEM based on dense point cloud data | |
| CN113124873B (en) | UUV multi-index constraint three-dimensional route planning method based on marine environment information | |
| KR101522203B1 (en) | Method for generating synthesis environment using warm eddy/internal wave/submarine topography | |
| CN114137566A (en) | Method for inverting rock movement parameters of ponding subsidence basin | |
| CN105931268A (en) | Mean shift tracking method based on scale adaption in UUV underwater recovery process |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20230425 |