Disclosure of Invention
The invention aims to overcome the defects of the prior art, provides a river channel erosion and deposition measuring method based on three-dimensional laser scanning, realizes river channel erosion and deposition measurement by combining a three-dimensional laser scanning technology, and has the advantages of high efficiency, small workload, high precision, low risk and the like.
The purpose of the invention is realized by the following technical scheme: a river channel erosion and deposition measuring method based on three-dimensional laser scanning comprises the following steps:
s1, river terrain measurement: determining the plane arrangement of a scanner erection station and a target according to the geometric shape of a measurement area, and acquiring point cloud data which can meet the erosion-deposition deformation measurement requirement and completely cover the scanning area;
s2, point cloud data registration: firstly, multi-site point cloud registration of a first period is carried out, preliminary registration is carried out through a target, registration errors are further reduced through a closest point iterative algorithm, and after completion, point cloud is exported and named by date; then, performing multi-period point cloud registration, using a first-period point cloud coordinate system as an engineering coordinate system, registering a second-period point cloud through a closest point iterative algorithm, and finishing deriving the second-period point cloud, wherein the second-period point cloud is named by date;
s3, calculating and visualizing erosion-deposition deformation of the river channel: through the altitude variation value of the same point of the plane coordinate in different periods of contrast, the erosion and deposition variable quantity of the river channel is obtained: because plane coordinates of point clouds in different periods are not coincident, difference calculation cannot be directly carried out, firstly, deleting the point clouds in regions outside an original point cloud river channel, and extracting a river channel point cloud boundary; then, generating a plane coordinate matrix with equal spacing, and deleting points outside the boundary to obtain a new matrix; then, obtaining a height coordinate corresponding to each point through interpolation, and storing the height coordinates in a matrix form; and finally, calculating the height difference of the corresponding points, and displaying the junction in the form of a deformed cloud picture or a contour map.
Wherein the step S1 includes the following substeps:
s101, determining an effective measurement distance rho of a scanner according to the following formula:
N=1/(ρ×h×((tan(θ+dθ)-tanθ)×dα);
wherein θ is arccos (h/ρ), N is the minimum point cloud density, h is the height from the center of the instrument to the riverbed, ρ is the effective measurement distance of the instrument, and da and do are the horizontal resolution and the vertical resolution of the scanner;
s102, acquiring the geometric shape of the measured river reach by using the satellite image, measuring the data of the width and the length of the river reach, and calculating the maximum distance L between scanner standsmax;
Wherein d is 1/2 of the river reach width;
arranging the erection stations on the image according to the maximum distance between the erection stations, and alternately arranging the stations on two banks to avoid instrument scanning blind areas;
s103, erecting an instrument at each preset station position, and arranging at least 3 targets between two different testing stations; and then scanning at different stations until the measurement of the whole area is completed, and acquiring point cloud data which meets the erosion and deposition deformation measurement requirement and completely covers the scanning area.
Wherein the step S2 includes the following substeps:
s201, importing first-stage point cloud data, and carrying out point cloud denoising and vegetation removal; taking the scanner coordinate system of the first station as an engineering coordinate system, taking the target coordinate in the first station as an engineering coordinate, and firstly, calculating a coordinate conversion parameter by solving the following formula:
A=RA'+ΔT
in the formula:
r comprises three independent parameters alpha, beta,
Rotation angles around x, y, z axes, respectively;
Δ X, Δ Y, Δ Z are three translation parameters;
engineering coordinates for the target;
scanner coordinates for the target;
the simultaneous solution is carried out to obtain six parameters of alpha, beta,
Delta X, delta Y and delta Z, completing the initial registration of the second station to the engineering coordinate system, and then reducing the registration error by adopting a closest point iterative algorithm to realize the final registration; and repeating the registration of the 3 rd, 4 th,.., n stations;
s202, importing second-stage point clouds, denoising and removing vegetation, importing the first-stage registered point clouds in the same project, deleting the point clouds in the river channel, and reserving points of which the external topography of the river channel is not changed; secondly, registering the clouds of the first station in the second period to the clouds of the first station by using the external area of the riverway overlapped by the point clouds in the two periods by adopting a closest point iterative algorithm; finishing the registration of the remaining point clouds in the second period according to the same method;
s203, deleting the point clouds in the region outside the river channel from the two-stage point clouds after registration, reserving the rest point clouds according to the convex polygon region, respectively exporting point cloud coordinate information in a txt format, and naming by using dates.
The invention has the beneficial effects that: the river channel erosion and deposition measurement method is combined with a three-dimensional laser scanning technology to achieve river channel erosion and deposition measurement, and has the advantages of being fast in data acquisition, high in measurement accuracy, simple and convenient to operate and small in workload.
Detailed Description
The technical solutions of the present invention are further described in detail below with reference to the accompanying drawings, but the scope of the present invention is not limited to the following.
The invention adopts a three-dimensional laser scanner to obtain the river channel point cloud information. The measurement objects are exposed riverbed and bank slope. Because the emitted laser is refracted on the water surface, the intensity of the return signal received by the instrument is very weak, and the terrain below the water surface is difficult to measure, so the measurement is recommended to be carried out in a dry season or when the exposed area of a riverbed is large in the specific implementation process.
As shown in fig. 1, a river channel erosion and deposition measuring method based on three-dimensional laser scanning includes the following steps:
s1, river terrain measurement: according to the geometric shape of the measurement area, determining the plane arrangement of the scanner erection station and the target, and acquiring point cloud data which can meet the erosion-deposition deformation measurement requirement and completely cover the scanning area: in the embodiment of the application, the measurement area is selected as a downstream river channel of a certain sluice engineering, and the downstream river channel is straight. The measuring instrument adopts a Ruigger VZ-2000i laser scanner, the distance measurement reaches 2000m furthest, the speed measurement reaches 3000 ten thousand points/minute, the maximum horizontal angle resolution is 0.0007 degrees, and the maximum vertical angle resolution is 0.0015 degrees;
the step S1 includes the following sub-steps:
s101, considering that the farther the relative distance between the scanner and the scanner is, the smaller the point cloud density is, the effective measurement distance of the scanner is determined according to the minimum point cloud density required by deformation observation before measurement. For a river channel with large surface relief, a larger point cloud density can be selected, for a river channel with a flat surface, a smaller point cloud density can be selected, and the effective measurement distance rho of the scanner is determined according to the following formula:
N=1/(ρ×h×((tan(θ+dθ)-tanθ)×dα);
wherein θ is arccos (h/ρ), N is the minimum point cloud density, h is the height from the center of the instrument to the riverbed, ρ is the effective measurement distance of the instrument, and da and do are the horizontal resolution and the vertical resolution of the scanner;
in the above embodiment, the minimum point cloud density required for observation is 1000 points/m2. The height h from the center of the instrument to the river bed surface is the height of the embankment plus the height of a tripod, and is about 20 m. The horizontal resolution of the scanner was set to 0.004 °, and the vertical resolution was set to 0.002 °. Converting the angle into radian when calculating the effective measuring distance; for the convenience of solution, firstly, a rho value is assumed to obtain a point cloud density N, and if the point cloud density is close to the set minimum point cloud density, the assumed rho value is used as an approximate value of the effective observation distance. Assuming that ρ is 200m, N is calculated to be about 1100 points/m2. Therefore, the effective measuring distance is 200 m;
s102, acquiring the measured river by using satellite imagesMeasuring the geometric shape of the section, measuring the width and length data of the river section, and calculating the maximum spacing L of the scanner erecting stationsmax;
Wherein d is 1/2 of the river reach width;
arranging the erection stations on the image according to the maximum distance between the erection stations, and alternately arranging the stations on two banks to avoid instrument scanning blind areas;
in the above embodiment, the river reach width is about 300m, the length is about 650m, etc. measured from the satellite image, Lmax is calculated to be 132m, 12 stations are arranged at equal intervals of 130m on both banks, and 6 stations are respectively arranged on the right bank and the left bank;
s103, erecting an instrument at each preset station position, and arranging at least 3 targets between two different testing stations; and then scanning at different stations until the measurement of the whole area is completed, and acquiring point cloud data which meets the erosion and deposition deformation measurement requirement and completely covers the scanning area.
In the above embodiment, 3 targets were placed between two different stations. For example, when scanning is performed at the second station, 3 targets are arranged between the first and second stations and between the second and third stations. And after the scanning is finished, if the missing area is found, immediately performing the additional measurement. Likewise, the second stage scanning is done at the same site location.
S2, point cloud data registration: firstly, multi-site point cloud registration of a first period is carried out, preliminary registration is carried out through a target, registration errors are further reduced through a closest point Iteration (ICP) algorithm, and point cloud is exported and named by date after completion; then, multi-period point cloud registration is carried out, a first-period point cloud coordinate system is used as an engineering coordinate system, a second-period point cloud is registered through a closest point iterative algorithm, the second-period point cloud is derived, and the method is named by date:
specifically, in the embodiment of the present application, the step S2 includes the following sub-steps:
s201, a project is newly built in point cloud processing software, point cloud data in the same period are imported, point cloud denoising is carried out, and vegetation is removed. Taking the scanner coordinate system of the first station as an engineering coordinate system, taking the target coordinate in the first station as an engineering coordinate, and firstly, calculating a coordinate conversion parameter by solving the following formula:
A=RA'+ΔT
in the formula:
r comprises three independent parameters alpha, beta,
Rotation angles around x, y, z axes, respectively;
Δ X, Δ Y, Δ Z are three translation parameters;
engineering coordinates for the target;
scanner coordinates for the target;
the simultaneous solution is carried out to obtain six parameters of alpha, beta,
Delta X, delta Y and delta Z, completing the initial registration of the second station to the engineering coordinate system, and then reducing the registration error by adopting a closest point iterative algorithm to realize the final registration; and repeating the registration of the 3 rd, 4 th,.., n stations; in the above embodiment, n is 12.
S202, building a new project, importing second-stage point clouds, denoising and removing vegetation, importing the first-stage registered point clouds in the same project, deleting the point clouds in the river channel, and reserving points of which the external topography of the river channel is not changed; then, registering the clouds of the first station in the second stage to the clouds of the first stage by using an external region of the riverway overlapped by the clouds of the two stages by adopting a closest point Iteration (ICP) algorithm; finishing the registration of the remaining point clouds in the second period according to the same method;
s203, deleting the point clouds in the region outside the river channel from the two-stage point clouds after registration, reserving the rest point clouds according to the convex polygon region, respectively exporting point cloud coordinate information in a txt format, and naming by using dates.
S3, calculating and visualizing erosion-deposition deformation of the river channel: through the altitude variation value of the same point of the plane coordinate in different periods of contrast, the erosion and deposition variable quantity of the river channel is obtained: because plane coordinates of point clouds in different periods are not coincident, difference calculation cannot be directly carried out, firstly, deleting the point clouds in regions outside an original point cloud river channel, and extracting a river channel point cloud boundary; then, generating a plane coordinate matrix with equal spacing, and deleting points outside the boundary to obtain a new matrix; then, obtaining a height coordinate corresponding to each point through interpolation, and storing the height coordinates in a matrix form; and finally, calculating the height difference of the corresponding points, and displaying the junction in the form of a deformed cloud picture or a contour map.
Specifically, in the embodiment of the present application, the step S3 includes:
s301: and reading txt format point cloud in matlab software, storing the point cloud in a 4-row matrix, and naming the point cloud as pointclosed _1 and pointclosed _2, wherein the first row is the serial number of the point, and the second, third and fourth rows are the x, y and z coordinates of the point. And then, a delaunay algorithm is used for generating a triangular net on the basis of the first-stage point cloud, and the side and the triangular surface are numbered respectively. All the edges are stored in a 3-column matrix line _1, the first column is the serial number of the edge, and the second and third columns are the serial numbers of two points included in the edge. All triangle faces are stored in a 4-column matrix tri _1, the first column is the number of the triangle face, and the second, third and fourth are the serial numbers of three sides of each triangle face. And counting the serial numbers of the edges in the matrix for storing the triangular surface information, wherein the serial number is the inner edge of the triangular network when the serial number appears for 2 times, and the serial number is the edge on the boundary of the triangular network when the serial number appears for 1 time. And extracting point cloud coordinates of points forming the outer boundary edge, storing the point cloud coordinates as a matrix bound _1, and reordering the point sequence in the matrix counterclockwise.
S302: firstly, a meshgrid function is utilized to generate grid coordinate matrixes x _ coord and y _ coord, x-direction coordinates and y-direction coordinates of grid points are stored respectively, and coordinate increment is 0.03. And then, interpolating the vertical coordinates of the grid points on the basis of the point cloud pointclose _1 and pointclose _2 by using a gritdata function, wherein the interpolation method can be a bilinear interpolation method, a cubic spline interpolation method, a bicubic interpolation method and the like. And obtaining interpolated z-coordinate matrixes z1_ coord and z2_ coord, and subtracting the former from the latter to obtain a z-coordinate change value matrix delz _ coord.
S303: the reshape function is then applied to straighten the x _ coord, y _ coord, delz _ coord into column vectors, which are combined and stored as a matrix gridcoord. And then, judging the positions of the points in the gridcoord and the boundary bound _1 by using an inpolygon function, and storing the points inside the boundary in a matrix in _ gridcoord.
S304: and exporting the coordinates in the in _ gridcoord as a txt file, and drawing a erosion deformation cloud picture or a contour map through importing drawing software, as shown in fig. 2 and 3.
The sorting method for S032 comprises the following steps: firstly, calculating barycentric coordinates according to plane coordinates of boundary points, then generating vectors from the barycenter to all the boundary points, calculating included angles between all the vectors and the x direction, and finally sequencing the boundary points according to the sequence of the angles from small to large.
The invention is based on three-dimensional laser scanning, thus improving the data acquisition speed; the measurement precision is high, the error generated by contact measurement is reduced, and the precision of laser scanning can reach mm level; compared with the traditional single-point measurement, the deformation value of any point of the river channel plane can be obtained; the method is simple and convenient to operate, has small workload, and can be realized by using a point cloud registration function of software and a simple matlab statement; the visualization of results can be realized through a software drawing function, and the distribution condition and the deformation amount of the erosion and deposition deformation area of the river channel are vividly displayed.
The foregoing is a preferred embodiment of the present invention, it is to be understood that the invention is not limited to the form disclosed herein, but is not to be construed as excluding other embodiments, and is capable of other combinations, modifications, and environments and is capable of changes within the scope of the inventive concept as expressed herein, commensurate with the above teachings, or the skill or knowledge of the relevant art. And that modifications and variations may be effected by those skilled in the art without departing from the spirit and scope of the invention as defined by the appended claims.