Seawall terrain extraction method and system fusing vehicle-mounted laser point cloud and tilt model

Abstract

The application provides a kind of sea embankment terrain extraction method and system for fusing vehicle-mounted laser point cloud and tilt model, which comprises the following steps: forming a sea embankment range according to a sea embankment vector map and satellite imagery; dividing the sea embankment range into a vehicle-mounted laser radar collection range and an unmanned aerial vehicle oblique photogrammetry range; obtaining vehicle-mounted laser point cloud collected by vehicle-mounted laser radar, and pre-processing and precision optimizing the vehicle-mounted laser point cloud; obtaining oblique photogrammetry image data collected by a five-puzzle camera carried by an unmanned aerial vehicle, and generating an oblique photogrammetry model and an oblique three-dimensional model point cloud therefrom; extracting edge line elements and corner feature points in the common overlapping area of the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud, and performing point cloud fusion based on the edge line elements and corner feature points; filtering the fused point cloud to obtain sea embankment terrain point cloud data; and constructing a resolution-adaptive digital elevation model for different regions of the sea embankment dam based on the sea embankment terrain point cloud data, which can accurately represent the structural characteristics of the sea embankment.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of computer, and in particular to a sea embankment terrain extraction method and system fusing vehicle-mounted laser point cloud and inclined model. BACKGROUND

[0002] As an important coastal protection project, the accurate acquisition of the terrain data of the sea embankment is crucial for coastal safety monitoring, disaster assessment and engineering maintenance. The traditional terrain extraction method mainly relies on a single data source: ① Vehicle-mounted laser radar is suitable for flat areas such as embankment roads, and can obtain high-precision point cloud, but there are blind areas on the steep side of the sea embankment (such as revetment and wave dissipation structure), and it is limited by the vehicle traffic conditions; ② Inclined photogrammetry can cover the vertical surface of the sea embankment by reconstructing a three-dimensional model from multi-angle images, but the model accuracy is limited by the image resolution, and there are voids in weak texture areas (such as cement pavement), and the point cloud density is insufficient.

[0003] There are significant defects in fusing vehicle-mounted laser point cloud and inclined model data using existing methods:

[0004] ① Low data registration accuracy. The control point layout does not consider the linear features of the sea embankment, resulting in inconsistent spatial reference of vehicle-mounted point cloud and inclined model;

[0005] ② Poor fusion effect. Traditional point cloud fusion ignores the feature structures such as the edge and corner points of the sea embankment, resulting in elevation jumps at the connection between the embankment top and the vertical surface;

[0006] ③ Inaccurate terrain classification. The general filtering algorithm cannot distinguish between the wall surface of the sea embankment, vegetation and ground noise points, affecting the accuracy of the DEM. SUMMARY

[0007] In view of the above problems, the present application is proposed in order to provide a sea embankment terrain extraction method and system fusing vehicle-mounted laser point cloud and inclined model which overcomes the above problems or at least partially solves the above problems.

[0008] In one aspect of the present application, a sea embankment terrain extraction method fusing vehicle-mounted laser point cloud and inclined model is provided, the method comprising:

[0009] Extending a buffer zone to both sides of the sea embankment center line in the satellite image according to the sea embankment vector map to form a sea embankment range;

[0010] Dividing the sea embankment range into a vehicle-mounted laser radar collection range and an unmanned aerial vehicle inclined photogrammetry range;

[0011] Obtaining vehicle-mounted laser point cloud collected by the vehicle-mounted laser radar, performing unified spatial reference coordinate conversion and denoising processing on the vehicle-mounted laser point cloud, and performing precision optimization on the denoised vehicle-mounted laser point cloud;

[0012] Obtain the oblique photography image data and POS data collected by the unmanned aerial vehicle equipped with a five-pen camera, perform unified space reference coordinate conversion on the POS data, and generate an oblique photography model and an oblique three-dimensional model point cloud according to the converted POS data and the oblique photography image data;

[0013] Extract edge line elements and corner feature points in the common overlapping area of the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud, and fuse the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner feature points to obtain a fused point cloud;

[0014] Filter the fused point cloud to screen out dike terrain point cloud data;

[0015] According to the point cloud average density of the dike terrain point cloud data and the planar width of the dike dam, a digital elevation model with adaptive resolution is constructed for different regions of the dike dam as structural features.

[0016] Another aspect of the present application also provides a dike terrain extraction system for fusing vehicle-mounted laser point cloud and oblique model, the system comprises:

[0017] A dike range drawing module is configured to extend a buffer zone on both sides of the dike center line in the satellite image according to the dike vector map to form a dike range.

[0018] A measurement range division module is configured to divide the dike range into a vehicle-mounted laser radar collection range and an unmanned aerial vehicle oblique photogrammetry range.

[0019] A vehicle-mounted laser point cloud processing module is configured to obtain vehicle-mounted laser point cloud collected by a vehicle-mounted laser radar, perform unified space reference coordinate conversion and denoising processing on the vehicle-mounted laser point cloud, and perform precision optimization on the denoised vehicle-mounted laser point cloud.

[0020] An oblique photography data processing module is configured to obtain oblique photography image data and POS data collected by an unmanned aerial vehicle equipped with a five-pen camera, perform unified space reference coordinate conversion on the POS data, and generate an oblique photography model and an oblique three-dimensional model point cloud according to the converted POS data and the oblique photography image data.

[0021] A point cloud fusion module is configured to extract edge line elements and corner feature points in the common overlapping area of the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud, and fuse the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner feature points to obtain a fused point cloud.

[0022] A fused point cloud filtering module is configured to filter the fused point cloud to screen out dike terrain point cloud data.

[0023] The sea dike terrain extraction module is used for constructing a resolution-adaptive digital elevation model according to the point cloud average density of the sea dike terrain point cloud data and the planar width of the sea dike dam to take different regions of the sea dike dam as structural features.

[0024] Another aspect of the present application also provides a computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor; the computer program, when executed by the processor, implements the steps of the sea dike terrain extraction method of fusing the vehicle-mounted laser point cloud and the inclined model as described above.

[0025] Another aspect of the present application also provides a computer program product having a computer program stored thereon, the computer program, when executed by a processor, implements the steps of the sea dike terrain extraction method of fusing the vehicle-mounted laser point cloud and the inclined model as described above.

[0026] The sea dike terrain extraction method and system of fusing the vehicle-mounted laser point cloud and the inclined model provided by the embodiments of the present application, in the data acquisition link, perform unified spatial reference coordinate conversion and denoising processing on the vehicle-mounted laser point cloud and the inclined photographic image data, and perform precision optimization on the denoised vehicle-mounted laser point cloud, in the multi-source fusion stage, by extracting edge line elements and corner feature points of the public overlapping region in the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud, fusing the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner feature points to obtain fused point cloud, realizing automatic weighted fusion of the point cloud, filtering the fused point cloud to obtain sea dike terrain point cloud data related to the sea dike terrain, and constructing a resolution-adaptive digital elevation model according to the point cloud average density of the sea dike terrain point cloud data and the planar width of the sea dike dam to take different regions of the sea dike dam as structural features, which can accurately represent the structural features of the sea dike, provide sub-meter data support for coastal safety monitoring and disaster assessment, and fill the technical gap of the current specification for steep slope terrain measurement.

[0027] The above description is only a summary of the technical solutions of the present application, in order to more clearly understand the technical means of the present application, the specific embodiments of the present application can be implemented according to the content of the specification, and in order to make the above and other purposes, features and advantages of the present application more obvious and easy to understand, the following specific embodiments of the present application are described. BRIEF DESCRIPTION OF DRAWINGS

[0028] Various other advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments. The accompanying drawings are included to provide a description of the preferred embodiments and are not meant to limit the present application. Furthermore, the same reference numerals are used throughout the several drawings to represent similar components. In the drawings:

[0029] Figure 1A flowchart of the sea dike terrain extraction method fusing vehicle-mounted laser point cloud and an inclined model is provided for an embodiment of the present application.

[0030] Figure 2 A structural block diagram of the sea dike terrain extraction system fusing vehicle-mounted laser point cloud and an inclined model is provided for an embodiment of the present application. DETAILED DESCRIPTION

[0031] Exemplary embodiments of the present disclosure will be described in greater detail below with reference to the accompanying drawings. Although exemplary embodiments of the present disclosure are shown in the drawings, it is understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided so that the present disclosure can be more thoroughly understood, and the scope of the present disclosure can be accurately conveyed to those skilled in the art.

[0032] Those skilled in the art can understand that, unless specifically stated, the singular forms "a", "an", "said" and "the" used herein also include the plural forms. It should be further understood that the phrase "comprising" used in the specification of the present application means that the features, integers, steps, operations, elements and / or components exist, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or groups thereof.

[0033] Those skilled in the art can understand that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as that generally understood by those skilled in the art to which the present application belongs. It should also be understood that terms such as those defined in a general dictionary should be understood to have meanings consistent with those in the context of the prior art, and should not be interpreted in an idealized or overly formal sense unless specifically defined.

[0034] Figure 1 A flowchart of the sea dike terrain extraction method fusing vehicle-mounted laser point cloud and an inclined model of one embodiment of the present application is schematically shown. Referring to Figure 1 , the sea dike terrain extraction method fusing vehicle-mounted laser point cloud and an inclined model of an embodiment of the present application specifically includes the following steps:

[0035] S1, according to the sea dike vector map, a buffer zone is extended to both sides in the satellite image with the sea dike center line as the reference to form a sea dike range.

[0036] In this embodiment, according to the collected sea dike vector map, a high-resolution satellite image is superimposed in GIS (Geographic Information System) software, a buffer is extended to both sides with the sea dike as the center to draw the range on both sides, to form a sea dike range, i.e. a sea dike measurement area, and is converted into an ESRI shape file.

[0037] Then control point layout, plane control measurement, elevation control measurement, conversion relationship calculation and image control point brushing measurement are carried out. Among them, the control point layout is specifically: according to the linear range of the seawall, the control points are laid out on both sides of the seawall within 1 kilometer according to the layout distance; in the complex terrain range of the seawall, the number of control points is encrypted in the overlapping area of vehicle-mounted measurement and tilt measurement, and the number of control point positions is greater than or equal to 4, and the average side length L is obtained by the seawall mileage and the number of laid points. Determine L and the distance between adjacent points, when L is less than 2000 meters or the distance between adjacent points is less than 1000 meters or greater than 4000 meters, the control point position needs to be re-laid. The plane control measurement is specifically: selecting a standard initial point in the range near the seawall, observing the control point statically to obtain static base station data, and then processing the baseline to solve the seawall control point in the industry to obtain the plane coordinates of the control point. The elevation control measurement is: selecting an elevation initial point in the range near the seawall, performing leveling measurement to obtain observation data, and calculating the leveling route closure error in the industry to obtain the elevation coordinates of the control point. The elevation control measurement is specifically: selecting an elevation initial point in the range near the seawall, performing leveling measurement to obtain observation data, and calculating the leveling route closure error in the industry to obtain the elevation coordinates of the control point. The conversion relationship calculation is specifically: arranging the latitude and longitude of the control points and the two sets of data of the results, calculating the plane conversion relationship to obtain seven parameters, and calculating the elevation anomaly model to obtain the elevation conversion relationship. The image control point brushing measurement is specifically: setting the length of the square side, spraying two white squares with white paint, and the shape is two squares intersecting at a diagonal, as the model of the image control point. In the seawall measurement area, the image control points are uniformly laid according to the image control point layout distance and the precision check point layout distance. The GNSS receiver parameters are set in the field, the GNSS-RTK mode is adopted, the antenna of the GNSS receiver is placed on the image control point, and the geodetic coordinates and the basic control coordinate system results of the image control point are obtained. The length of the square side is 0.3 meters, the image control point layout distance is 500m, and the precision check point layout distance is 2000m.

[0038] The application effectively solves the uneven distribution of linear engineering control points by the dynamic encryption rule of control points (more than 4 points in complex terrain area and less than 2000m of L), avoids the precision fluctuation caused by traditional equal interval layout, and improves the plane precision by more than 30%. The seven-parameter conversion model (Bursa-Wolf) is combined with the quadratic surface fitting of elevation anomaly to realize the millimeter-level conversion from WGS84 coordinate system to engineering coordinate system, the plane conversion residual is less than 0.02m, and the fitting error in the elevation anomaly is less than 2cm.

[0039] S2, the range of the seawall is divided into a vehicle-mounted laser radar collection range and an unmanned aerial vehicle tilt photogrammetry range.

[0040] In this embodiment, the seawall range can be divided into the vehicle-mounted laser radar collection range and the unmanned aerial vehicle oblique photogrammetry collection range according to the seawall road information. Specifically, the boundary range of the seawall range is determined, high-resolution satellite images and DEM data are combined, the road width, terrain undulation and road traffic condition beside the seawall are analyzed, when the above conditions meet the vehicle-mounted laser radar vehicle driving, the red driving line is drawn, otherwise the green driving line is drawn, and the green driving line is extended to the red driving line by 500 meters. The red driving line is taken as the seawall range collected by the vehicle-mounted laser radar; the green driving line is buffered by 100 meters to form a planar vector data, which is taken as the collection range of the unmanned aerial vehicle oblique photogrammetry.

[0041] S3, acquiring the vehicle-mounted laser point cloud collected by the vehicle-mounted laser radar, performing unified space reference coordinate conversion and denoising processing on the vehicle-mounted laser point cloud, and performing precision optimization on the denoised vehicle-mounted laser point cloud.

[0042] S4, acquiring the oblique photogrammetry image data and POS data collected by the unmanned aerial vehicle carrying five-puzzle cameras, performing unified space reference coordinate conversion (plane coordinate conversion and elevation coordinate conversion) on the POS data, generating an oblique photogrammetry model and an oblique three-dimensional model point cloud according to the converted POS data and the image data;

[0043] S5, extracting the edge line elements and the corner feature points in the common overlapping area of the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud, and fusing the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner feature points to obtain a fused point cloud.

[0044] S6, filtering the fused point cloud to screen out seawall terrain point cloud data.

[0045] S7, constructing a resolution-adaptive digital elevation model according to the point cloud average density of the seawall terrain point cloud data and the plane width of the seawall dam, and taking the structural features as different regions of the seawall dam.

[0046] The sea embankment terrain extraction method and system provided by the embodiment of the present application fuse the vehicle-mounted laser point cloud and the inclined model, in the data acquisition link, the vehicle-mounted laser point cloud and the inclined photographic image data are subjected to unified space reference coordinate conversion and denoising processing, and the vehicle-mounted laser point cloud after denoising is subjected to precision optimization, in the multi-source fusion stage, the edge line elements and the corner feature points of the public overlapping area in the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud are extracted, the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud are fused based on the edge line elements and the corner feature points to obtain fused point cloud, the automatic weighted fusion of the point cloud is realized, the fused point cloud is filtered to screen out the sea embankment terrain point cloud data related to the sea embankment terrain, and the resolution adaptive digital elevation model is constructed for different regions of the sea embankment dam according to the point cloud average density of the sea embankment terrain point cloud data and the planar width of the sea embankment dam, the structure characteristics of the sea embankment are accurately characterized, the sub-meter data support is provided for the coastal safety monitoring and disaster assessment, and the technical blank of the current specification in the steep slope terrain measurement is filled. The problems of low registration accuracy, fusion edge jump, terrain classification distortion and the like in the traditional sea embankment terrain extraction method caused by a single data source are solved.

[0047] In the embodiment of the present application, the vehicle-mounted laser point cloud collected by the vehicle-mounted laser radar in step S3 includes: collecting satellite data by using a static observation method on the control point of the GNSS receiver half an hour before the vehicle-mounted laser radar is collected, and closing the receiver half an hour after the vehicle-mounted laser radar is collected to obtain static data. The pulse frequency, scanning frequency and vehicle driving speed of the vehicle-mounted laser radar device are set, the vehicle is driven according to the planned route, the laser radar transmitter is turned on, the laser beam is emitted and the reflected light is received, the laser distance data is obtained by measuring the flight time of the laser pulse; GNSS (Global Navigation Satellite System) and IMU (Inertial Measurement Unit) data are obtained at the same time to form POS data; and the collected distance data and POS data are stored in the local storage device of the vehicle-mounted computer in real time.

[0048] In the embodiment of the present application, the unified space reference coordinate conversion and denoising processing of the vehicle-mounted laser point cloud in step S3 specifically includes: solving the POS data in combination with the GNSS static base station observation file and the coordinate data of the control point in the WGS84; converting the distance data from the polar coordinates to the laser radar rectangular coordinate system, to the vehicle coordinate system, to the geodetic coordinate system and to the control survey coordinate system; and then denoising the vehicle-mounted laser point cloud data by using the distance threshold denoising method and the standard deviation denoising method. The specific implementation process includes.

[0049] (1) POS data post-difference;

[0050] Combined with GNSS static base station observation file, control point WGS84 coordinate data (B, L, H), the POS is post-differential calculated to obtain accurate trajectory data, including time, longitude, latitude, height, roll angle, pitch angle and heading angle.

[0051] (2) Polar coordinate to laser radar rectangular coordinate conversion of distance data;

[0052] In the laser radar coordinate system, the distance data is converted from the polar coordinate system point cloud (R, θ, φ) to the laser radar rectangular coordinate system point cloud (X3, Y3, Z3) combined with the emission angle of the laser beam and other information, wherein R, θ and φ represent the distance between the laser radar and the target object, the emission angle of the laser beam in the horizontal and vertical directions, respectively.

[0053] (3) Conversion of laser radar rectangular coordinate system to vehicle coordinate system;

[0054] The rectangular coordinate system (X3, Y3, Z3) in the laser radar coordinate system is converted to the vehicle coordinate system point cloud (X V , Y V , Z V ) through the rotation matrix (R V ) and the translation vector (T V ) obtained by external parameter calibration, wherein V represents the vehicle, and (X V , Y V , Z V ) represents the coordinate value in the vehicle coordinate system;

[0055] .

[0056] (4) Conversion of vehicle coordinate system to geodetic coordinate system;

[0057] According to the vehicle attitude angle and position information provided by the POS system, the point cloud coordinates (X V , Y V , Z V ) in the vehicle coordinate system are converted to the geodetic coordinate system, first a rotation matrix (R , ω, κ) is constructed according to the attitude angle, and then combined with the position (X0, Y0, Z0) in the geodetic coordinate system, the geodetic coordinate system point cloud coordinates (X g , Y g , Z g ) are calculated;

[0058] .

[0059] (5) Conversion of geodetic coordinate system to control survey coordinate system;

[0060] Based on the seven parameters calculated in the control survey, the point cloud data in the geodetic coordinate system is transformed to the control survey coordinate system; then, based on the elevation anomaly fitting model, the point cloud in the control survey coordinate system is transformed from the geodetic elevation system to the normal elevation system.

[0061] (6) Point cloud denoising; Point cloud denoising consists of two steps:

[0062] The first step is distance threshold denoising. For each point, a radius threshold R0 is set, and the number of points within this radius is counted. If the number of points is lower than the threshold, the point is considered noise and removed. For a point P in the point cloud... i (x) i ,y i ,z i Using it as the center, calculate the other points P. j (x) j ,y j ,z j ) to point P i distance d ij Determine point P. i To other points P j distance d ij Relationship with radius search threshold R0; when d ij If the value is greater than R0, it is not included in the count of point cloud numbers N0 that satisfy the relation; when d ij Points less than or equal to R0 are counted in the number of point clouds satisfying the relationship, N0. The number of point clouds satisfying the relationship, N0, is compared with a set threshold T0. When N0 is less than T0, point P is... i Points marked as noise are removed from the point cloud; points P are retained when N0 is greater than or equal to T0. i Repeat the above steps until all points in the point cloud have been traversed. The calculation formula is:

[0063] ;

[0064] The second step is standard deviation denoising.

[0065] Iterate through each point in the vehicle-mounted laser point cloud data and calculate the mean of the point cloud in the x, y, and z dimensions. , , ;

[0066] ;

[0067] ;

[0068] ;

[0069] Calculate the standard deviation of the point cloud in the x, y, and z dimensions. , , , set the appropriate multiple k0, the threshold range of x dimension is , the threshold range of y dimension is , the threshold range of z dimension is . Wherein k0 value range can be 1 meter to 5 meters.

[0070] ;

[0071] ;

[0072] ;

[0073] Iterate each point P (x, y, z) in the point cloud, judge whether its coordinate value is within the threshold range of the corresponding dimension. When the x coordinate of a certain point out of range, or the y coordinate out of range, or the z coordinate out of range, the point is regarded as a noise point and removed.

[0074] The number of point clouds in the vehicle-mounted laser point cloud data.

[0075] In the embodiment of the application, the laser point cloud precision optimization step further includes: taking the image control point coordinate as the center, setting a radius threshold to obtain an image control point set within a preset intensity threshold range, identifying the corner point by using an adaptive gray threshold and tension algorithm to obtain the image control point corner point coordinate; then obtaining the point cloud image control point center point according to the corner point coordinate, and forming a corresponding relationship with the image control point; combining the corresponding relationship and the algorithm, calculating the transformation parameter.

[0076] Specifically, in the embodiment of the application, the precision optimization of the denoised vehicle-mounted laser point cloud in step S3 includes the following steps not shown in the accompanying drawings:

[0077] S31, taking the denoised vehicle-mounted laser point cloud data as the original vehicle-mounted point cloud.

[0078] S32, taking the image control point coordinate pre-arranged in the seawall range as the center, and obtaining the image control point set within the preset intensity threshold range according to the preset radius threshold range in the original vehicle-mounted point cloud.

[0079] In this embodiment, the implementation process of step S32 further includes: taking the image control point coordinate pre-arranged in the seawall range as the center, extracting the laser point cloud data Plocal within the preset radius threshold R1, wherein the radius threshold R1 can be 2 meters. Calculate the intensity mean value With standard deviation , according to the intensity mean value With standard deviation Dynamically calculate the minimum intensity threshold I min And the maximum intensity threshold I max ;

[0080] ;

[0081] ;

[0082] ;

[0083] ;

[0084] I is the intensity value of the laser point cloud point; N1 is the number of laser point clouds within the radius threshold; n1 is a coefficient, n1 = 1.5-2.5.

[0085] Traverse the laser point cloud data within the preset radius threshold R2 range, if the intensity value of the current point is less than or equal to I max And greater than or equal to I min , the current point is added to the preset point set, and finally the image control point set Q is formed; wherein the radius threshold R2 can be 1 meter.

[0086] S33, using adaptive gray threshold and tension algorithm to identify corner points, obtaining image control point corner coordinates. The implementation process of step S33 specifically includes: using adaptive gray threshold and tension algorithm to identify corner points, obtaining image control point corner coordinates, including: two-dimensional projection of the image control point set Q, for each two-dimensional point after projection, according to the point cloud density and local terrain complexity, dynamically adjusting the adaptive window, calculating the gray gradient And And structure tensor M z , using structure tensor M z To construct corner point response function R z , and based on the corner point response function R z Determine whether the current two-dimensional point is an image control point corner point, and finally determine six two-dimensional image control point corners.

[0087] In this embodiment, for each point (x, y) of the two-dimensional projection of the image control point set Q, according to the point cloud density and local terrain complexity, the adaptive window is dynamically adjusted, and the gray gradient , And structure tensor M z .

[0088] ;

[0089] ;

[0090] ;

[0091] ;

[0092] ;

[0093] ;

[0094] , is the gradient of Sobel operator in x and y direction;

[0095] m, n are the row index number and column index number of the image pixel point;

[0096] I is the gray value or intensity value of the point;

[0097] is the gray value or intensity value of the point (x, y);

[0098] is the weight function, which is used to calculate the similarity between two pixel points G p (current pixel) and G q (another pixel in the neighborhood);

[0099] is the adaptive window, that is, the size and weight distribution of the window are controlled by the current pixel, ranging from 0.5 meters to 5 meters;

[0100] The corner response function R z is calculated again using the structure tensor M z . The average value of the corner response function R z of all points is calculated, and the response threshold T z is set to M z * . The relationship between the value of R z and the response threshold T z is determined. When R z > T z , the point is a gage point corner, otherwise, the point is not a gage point corner. Wherein k z is a constant greater than 1, M z = 1.5-3.

[0101] ;

[0102] ;

[0103] ;

[0104] is the structure vector M z .

[0105] is the trace of the structure vector M z .

[0106] k z is an empirical constant, generally taking a value between 0.04 and 0.06;

[0107] The six corner points in the two-dimensional projection are inversely calculated back to the control point coordinate system, and the six sets of corner point coordinates are divided into three sets of coordinate sets, three sets of center point coordinates are calculated, and the average value of the three pairs of center point coordinates is obtained as the image control point coordinate value. All image control points of the vehicle-mounted laser point cloud are compared with the laser point cloud data, and the image control point coordinate values of the problematic image control points are modified to obtain a new image control point set Q'.

[0108] S34, calculate the point cloud image control point center point according to the image control point corner point coordinates, select a target image control point from the image control point set measured in the field which matches the distance of the point cloud image control point center point, and construct a corresponding point pair between the point cloud image control point center point and the target image control point.

[0109] In this embodiment, the process of corresponding point searching is specifically as follows: selecting a point q i in the image control point set Q' of the laser point cloud, judging the distance between q i and p i in the image control point set P measured in the field; when the distance between q i and p i is less than or equal to the set distance threshold, q i and p i are one-to-one corresponding points; when the distance between q i and p i is greater than the set distance threshold, q i and p i are unrelated points, and the iteration is continued. The set distance threshold can be 0.5 meters to 2 meters.

[0110] S35, establish a corresponding relationship equation according to each corresponding point pair, and solve the transformation parameters between the point cloud image control point center point and the target image control point according to the corresponding relationship equation.

[0111] In this embodiment, according to each corresponding point (p i , q i ), an equation is established, and a target function is constructed, and the target function is respectively about R 1z and t 1zPartial derivative is taken and the partial derivative is set to zero to solve unknown coefficients of rotation matrix R 1z and translation vector t 1z The centroid of the laser point cloud image control point and the centroid of the field measured image control point are obtained, and the translation vector t 1z is obtained.

[0112] ;

[0113] ;

[0114] ;

[0115] ;

[0116] is the centroid of the field measured image control point set;

[0117] is the centroid of the laser point cloud image control point set;

[0118] is the objective function; is the number of corresponding points;

[0119] The translation vector equation is brought into the objective function to obtain a function about R 1z , and the optimal rotation matrix R 1z is obtained. 1z According to the equation of the translation vector t 1z , the translation vector t 1z is obtained.

[0120] ;

[0121] ;

[0122] ;

[0123] ;

[0124] ;

[0125] is the coordinate of the field measured image control point set after translation;

[0126] is the coordinate of the laser point cloud image control point set after translation;

[0127] H is the covariance matrix;

[0128] U and V are orthogonal matrices;

[0129] is a diagonal matrix.

[0130] S36, coordinate transform each point in the original vehicle-mounted point cloud according to the transformation parameters, and obtain a transformed point cloud.

[0131] In this embodiment, the implementation process of the point cloud transformation is as follows: the calculated rotation matrix R 1z and the translation vector t 1z are applied to each point q in the original vehicle-mounted laser radar point cloud Q to perform coordinate transformation, and an optimized point cloud q' is obtained, and the coordinate transformation formula is , wherein is a point in the optimized point cloud q'.

[0132] The transformed point cloud is used as a new original vehicle-mounted point cloud, and steps S32-S36 are repeatedly executed until the iteration termination condition is met, and the precision-optimized vehicle-mounted laser radar point cloud data is obtained, and the iteration termination condition is that the changes of each transformation parameter between two iterations are all less than the corresponding preset threshold value, or the coordinate difference between the point clouds before and after transformation is less than the preset error value (i.e. the value of the target function E(R 1z , t 1z ) is less than the corresponding error value, etc.).

[0133] The embodiment of the application also includes the step of laser point cloud precision checking. The specific implementation is as follows: the interval of the precision checking points is set, the landmark corner points near the seawall are collected in the field as the precision checking points, the precision checking points of the optimized transformed vehicle-mounted laser point cloud are manually selected, the field precision checking points and the precision checking points of the vehicle-mounted laser point cloud are compared and analyzed, the root mean square error and the average error are calculated, when the root mean square error and the average error are both less than or equal to the set threshold value, the laser point cloud precision meets the requirements; when the root mean square error and the average error are both greater than the set threshold value, the vehicle-mounted laser point cloud data is re-optimized.

[0134] The application adopts an adaptive denoising technology, a double denoising strategy of the distance threshold method (eliminating isolated points) combined with the standard deviation method (eliminating drift points), the noise elimination rate is > 95%, and the integrity of the seawall edge features is reserved. High-precision point cloud registration, automatic corner extraction algorithm based on intensity threshold and structure tensor, solves the registration problem in weak texture area. Dynamic precision optimization mechanism, iterative closest point (ICP) transformation parameter adaptive convergence, ensures the point cloud precision, meets the 1:500 topographic map specification.

[0135] In the embodiment of the present application, the obtaining of the oblique photography image data collected by the five-patch camera carried by the unmanned aerial vehicle in step S4 comprises: setting the flight height, the heading overlap and the lateral overlap according to the image resolution and the focal length of the camera, and generating the flight route according to the collection range of the seawall oblique photogrammetry, and obtaining the image data and the POS data of five different directions by using the five-patch camera carried by the unmanned aerial vehicle. The original image resolution range is 2-5 cm, and the flight height range is 150-500 m; the heading overlap is 70%-85%, and the lateral overlap is 70%-85%.

[0136] In the embodiment of the present application, the POS data is subjected to unified space reference coordinate conversion, and the relationship between the image position information and the image is matched, comprising:

[0137] In the data preprocessing step, the height anomaly model formula is used to convert the geodetic height system in the POS data into the normal height system; the corresponding relationship between the exterior orientation elements in the POS data and the original image is checked, and when there is a certain exterior orientation element in the POS data and there is no image in the original image, the exterior orientation element in the POS data is deleted; when there is no certain exterior orientation element in the POS data and there is an image in the original image, the image element is deleted.

[0138] In the coordinate system setting step, the corresponding ellipsoid is selected in combination with the coordinate system in the basic control measurement, and the central meridian of the seawall oblique aerial triangulation is set according to the central meridian of the control measurement.

[0139] In the embodiment of the present application, the generation of the oblique photography model and the oblique three-dimensional model point cloud according to the denoised oblique photography image data in step S4 comprises the following steps not shown in the accompanying drawings:

[0140] S41, corner point and edge feature extraction is performed on the oblique photography image data, SIFT algorithm is used to describe the key points of the image subjected to the corner point and edge feature extraction in different spatial scales, and the feature descriptor of the image is obtained.

[0141] The present application adopts the non-local mean denoising algorithm to perform denoising processing on the input image; the FAST algorithm and the Harris algorithm are used to identify the corner points, and then the gradient amplitude and direction of the image are calculated, and the double-threshold detection and edge connection are used to obtain the edge detection result; the image subjected to the corner point and edge feature extraction is input into the SIFT-Net to obtain the feature descriptor of the image. Specifically, it comprises:

[0142] (1) image preprocessing;

[0143] The input image is subjected to denoising processing, the non-local mean denoising algorithm is used to calculate the noise of each pixel, and the image grayscale value is mapped to the [0, 1] interval through grayscale normalization.

[0144] ;

[0145] G0 is the current pixel position, i.e. (x, y);

[0146] G1 is another pixel position within the neighborhood;

[0147] Ω is the neighborhood range;

[0148] I(G1) is the pixel value at y;

[0149] w(G0, G1) is the weight function;

[0150] C(G0) is the normalization constant;

[0151] h is a parameter to control the strength of denoising.

[0152] ;

[0153] I(x, y) is the original image pixel value,

[0154] I min and I max are the minimum and maximum gray scale values of the original image, respectively,

[0155] I norm (x, y) is the normalized pixel value.

[0156] (2) Feature extraction;

[0157] Corner extraction: the normalized image data uses the FAST algorithm to detect corners based on the local pixel gray scale change, and filters out candidate corners with higher probability; then, the Harris algorithm is combined to accurately locate the candidate corners, and the change of the correlation function in different directions is calculated to determine the corner.

[0158] Edge extraction: first, a Gaussian filter is used to smooth the image.

[0159] Edge extraction: then, the gradient amplitude and direction of the image are calculated, the gradient amplitude M(x, y) and direction θ(x, y) are calculated by finite difference of first-order partial derivative, non-maximum suppression is performed to retain the points with the largest local gradient and suppress other non-edge points. Then, double-threshold detection and edge connection are used, by setting a high threshold and a low threshold, points with gradient amplitude greater than the high threshold are regarded as strong edge points, points less than the low threshold are suppressed, and points between the two thresholds are judged as edge points according to their connectivity with strong edge points, to obtain the preliminary edge detection result.

[0160] ;

[0161] ;

[0162] G x (x, y) and G y (x, y) are the gradients of the image in x and y directions respectively.

[0163] For the preliminary edge detection result, the direction gradient histogram algorithm is used to further enhance the edge features and optimize the direction and intensity information of the edge.

[0164] (3) Feature matching

[0165] The image after corner and edge feature extraction is input into SIFT-Net, the key points are detected at different scales by constructing a scale space, and the descriptors of the key points are calculated, and a convolutional neural network is loaded. During the forward propagation process of the convolutional neural network, the image successively passes through each convolutional layer and pooling layer. At each layer, the features of the image are continuously extracted and abstracted. The extracted features are mapped to a fixed length vector space through a fully connected layer to obtain the feature descriptor of the image.

[0166] S42, according to the feature matching result and the ground control point information, the bundle adjustment algorithm is used to determine the exterior orientation elements of the aerial triangulation of the oblique photography image data. This step realizes not only the calculation of the exterior orientation elements of the aerial triangulation, but also the calculation of the three-dimensional coordinates of the encryption points. The feature matching result includes the corner points and the edge features; the aerial triangulation is calculated through the corner points and the edge features of each image. The encryption points are obtained by encrypting the points at any position in the image.

[0167] In this embodiment, according to the feature matching result and the ground control point information, the bundle adjustment algorithm is used to adjust the exterior orientation elements (including rotation and translation parameters) of each image and the three-dimensional point coordinates, so that the sum of the square errors between the projected image point coordinates and the actually observed image point coordinates is minimized, the position and attitude of the image in the three-dimensional space are determined, and the three-dimensional coordinates of the encryption points are obtained.

[0168] S43, according to the exterior orientation elements of aerial triangulation, dense matching is performed on the oblique photography image data, a dense oblique three-dimensional model point cloud is generated according to a preset laser point cloud density threshold, a mesh three-dimensional surface model is constructed by using the oblique three-dimensional model point cloud, texture information of the oblique photography image data is mapped to the mesh three-dimensional surface model by using a texture mapping technology, and an oblique photography model with real texture is generated. In this embodiment, according to the exterior orientation elements of aerial triangulation, dense matching is performed on the image, a laser point cloud density threshold is set, and a dense oblique three-dimensional model point cloud is generated. A triangular mesh is constructed by using the dense point cloud, a mesh three-dimensional surface model is obtained, and texture information of the image is mapped to the three-dimensional model by using a texture mapping technology, so that a three-dimensional real scene model with real texture is generated.

[0169] In S43, the mesh three-dimensional surface model is constructed by using the oblique three-dimensional model point cloud, the texture information of the oblique photography image data is mapped to the mesh three-dimensional surface model by using the texture mapping technology, and the oblique photography model with real texture is generated. Specifically, the following steps not shown in the accompanying drawings are included:

[0170] S431, an octree structure is used to construct a spatial index for the oblique three-dimensional model point cloud. The boundary range of the dense point cloud data in the x, y and z directions is obtained, a root node of a cube containing all the point clouds is created, the root node is divided into eight sub-nodes, and the number of point clouds in each sub-node is checked. A sub-node point cloud number threshold is set, if the number of point clouds in the sub-node is greater than the set threshold, recursive division is continued until the number of point clouds in each sub-node is less than or equal to the threshold or the maximum division layer is reached. (Different scene models have different numbers of sub-nodes, mainly including small scene, medium scene and large scene.)

[0171] S432, adjacent point searching is performed according to the spatial index, and a triangular mesh is constructed based on the searched candidate adjacent points. The triangular mesh is optimized based on the removal of narrow triangles and the edge exchange optimization principle.

[0172] In this embodiment, the steps of triangular mesh construction include:

[0173] Seed point selection: a point is randomly selected from the point cloud as a seed point, and the seed point is added to the initial triangular mesh.

[0174] Adjacent point searching: the spatial index is used to search for adjacent points within a certain distance range of the seed point. The distance between the seed point P1(X1, Y1, Z1) and the adjacent point P2(X2, Y2, Z2) is calculated, and the points with a relatively close distance are selected as candidate adjacent points according to the laser point cloud density threshold. When the distance between two points is less than or equal to the density threshold, the candidate adjacent point can be selected; when the distance between two points is greater than the density threshold, the candidate adjacent point is directly discarded.

[0175] Triangle mesh growth: select two points from the candidate neighboring points, form a triangle with the seed point, obtain the circumcircle equation of the triangle, then traverse other points in the point cloud, check whether other points are located inside the circumcircle.When there are other points falling within the circle, the triangle mesh does not meet the requirements, is directly discarded, and iteration calculation from the candidate neighboring points is continued;When there is no other point falling within the circle, the triangle mesh is directly added to the triangle mesh set of the seed point, the smallest angle in each triangle mesh is obtained from the triangle mesh set of the seed point, is sorted, and the triangle mesh with the maximum value of the smallest angle is selected as the final triangle mesh of the seed point;The final triangle mesh of the seed point is added to the point cloud triangle mesh set, and other points in the final triangle mesh of the seed point are selected as seed points, and iteration calculation is continued.

[0176] In the embodiment, the steps of triangle mesh optimization include:

[0177] Remove long and narrow triangle: calculate the edge length and internal angle of the triangle, and delete or adjust the long and narrow triangle with too short edge length or too small internal angle.Calculate the triangle surface, combine the triangle area and edge length to judge the shape of the triangle, and when the area is too small and the edge length ratio difference is too large, it is determined as a long and narrow triangle

[0178] Edge exchange optimization: traverse each edge in the triangle mesh, check the quadrilateral composed of adjacent triangles, and judge whether the edge exchange condition is met.If the minimum internal angle of two triangles can be increased after the edge exchange, the edge exchange operation is performed to optimize the quality of the triangle mesh.

[0179] ;

[0180] ;

[0181] ;

[0182] ;

[0183] ;

[0184] S A The triangle area is A, B and C are three angles of the triangle, and a, b and c are three sides of the triangle.

[0185] The application adopts spatial index octree partition (dynamic adjustment of child node threshold) combined with Delaunay triangle mesh local optimization, and the long and narrow triangle rejection rate is 100%, and the model geometric distortion rate is reduced by 70%.

[0186] S433, the texture information of the oblique photography image data is mapped to the optimized triangle mesh through the texture mapping technology, and the oblique photography model with real texture is obtained.

[0187] In this embodiment, for the optimized triangular net, one of the triangles is selected, the vertices of which are V1(x1, y1, z1), V2(x2, y2, z2) and V3(x3, y3, z3), and the corresponding texture coordinates are T1(u1, v1), T2(u2, v2) and T3(u3, v3). The barycentric coordinates (α1, β1, θ1) of a point P(x, y, z) in the triangle are calculated, and the texture coordinates (u, v) of the point P(x, y, z) are calculated, which are used to specify the corresponding position coordinates of the vertex on the three-dimensional model in the texture image. According to the texture coordinate value, the pixel value of the point in the texture is collected and assigned to the corresponding point in the triangular net.

[0188] ;

[0189] ;

[0190] ;

[0191] .

[0192] The embodiment of the application also includes the step of model precision checking. Specifically, the interval of the model precision checking points is set, the feature points near the seawall with small ups and downs are collected as the precision checking points in the field, and the point records of the precision checking points are made; the precision checking points are unfolded in the indoor work, the precision checking point positions of the inclined model are selected, the precision checking points of the field work and the precision checking points of the inclined model are compared and analyzed, the root mean square error and the average error are calculated, when the root mean square error and the average error are less than or equal to the model precision threshold, the precision of the inclined model meets the requirements; when the root mean square error or the average error is greater than the model precision threshold, the feature points and the inclined aerial triangulation are recalculated.

[0193] The embodiment of the application also includes the step of orthographic image production. The DEM of the regular grid is generated, and the elevation value of each grid is recorded; the corresponding position (x, y) of the image point on the ground is calculated by using the DEM and the image exterior orientation elements obtained by aerial triangulation, and the gray value of the position is assigned to the corresponding image element of the orthographic image; a plurality of corrected orthographic images are spliced to obtain the orthographic image. Specifically, the point cloud of the inclined three-dimensional model is interpolated to generate the digital elevation model of the regular grid, and the elevation value Z(x, y) of each grid is recorded;

[0194] The corresponding position (x, y) of the image point on the ground is calculated by using the DEM and the image exterior orientation elements obtained by aerial triangulation, and the gray value of the position is assigned to the corresponding image element of the orthographic image;

[0195] ;

[0196] ;

[0197] (x0, y0) is the image principal point coordinate;

[0198] f is the camera focal length;

[0199] (X S ,Y S ,Z S ) is the photograph center coordinate;

[0200] (x, y) is the image point coordinate;

[0201] (X W ,Y W ,Z W ) is the object point coordinate;

[0202] a i ,b i ,c i (i = 1, 2, 3, 4) is the element of the rotation matrix;

[0203] For the overlapping area, the weight is assigned according to the photograph center distance, the multiple corrected orthographic images are spliced, and the orthographic image is obtained.

[0204] ;

[0205] ;

[0206] ;

[0207] m i ,n i are integer coordinates of adjacent pixels in the original image;

[0208] I(m i ,n i ) is the gray value of the original image at (m i ,n i );

[0209] N is the number of images participating in fusion;

[0210] d k is the distance from the current point to the center of the kth image;

[0211] is a small constant to prevent division by zero.

[0212] This invention uses the covariance matrix to calculate the curvature of each point in the point cloud; it fits the quadratic surface equation, calculates the coefficient matrices of its first and second fundamental forms, and obtains the principal curvatures k1 and k2; it uses the curvature to extract edge and corner features from the vehicle-mounted seawall point cloud and the UAV tilted point cloud. In a specific embodiment, step S5, extracting the edge line features and corner feature points of the common overlapping area in the tilted 3D model point cloud and the vehicle-mounted laser point cloud, includes the following steps not shown in the attached figures:

[0213] S501. For each point in the tilted 3D model point cloud and the vehicle-mounted laser point cloud, obtain its neighborhood point set, construct the covariance matrix C' of the neighborhood points, and calculate the curvature k of the neighborhood points based on the covariance matrix C'.

[0214] In this embodiment, the neighborhood points are first calculated: for each point p in the point cloud... i Determine its k P Nearest neighbor or at P i The set of N neighborhood points within a radius r of the center Q .

[0215] Then, construct the covariance matrix: [construct the neighborhood points N]. Q The coordinate data is centered, that is, the coordinates of each point are subtracted from the coordinates of the center point P. i We obtain the centered coordinate data from the coordinates; then we calculate the covariance matrix C' based on this centered coordinate data.

[0216] ;

[0217] Where: n p is the number of neighboring points; i and j are the values ​​of x and y, respectively; It is the kth P The coordinates of the neighboring points in the i-th direction; It is the kth P The coordinates of the neighboring points in the j-direction; It is the mean of the coordinates along the i-th direction; It is the average of the coordinates in the j-direction.

[0218] Finally, the covariance moment C' matrix is ​​decomposed into eigenvalues ​​λ. 11 , λ 21 and λ 31 ;Calculate the curvature k based on the eigenvalues.

[0219] .

[0220] S502, for each point and neighborhood point in the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud, surface fitting is performed to obtain a local surface equation, and a first principal curvature k1, a second principal curvature k2 and respective principal curvature direction vectors are calculated according to the local surface equation .

[0221] In this embodiment, the local surface fitting specifically includes that for the point P and the neighborhood points thereof, a fitting equation is used, and coefficients are determined through a least square equation, so as to obtain a local surface equation.

[0222] ;

[0223] wherein a5, b5, c5, d5, e5 and f5 are coefficients in the equation.

[0224] In this embodiment, the calculation of the principal curvature direction specifically includes:

[0225] According to the obtained fitting quadratic surface equation, a coefficient matrix of a first fundamental form and a second fundamental form is calculated, and then a characteristic value equation is solved to obtain a first principal curvature k1 and a second k2 and a principal curvature direction vector .

[0226] ;

[0227] ;

[0228] ;

[0229] ;

[0230] wherein I is the first fundamental form; is a first-order partial derivative of the surface with respect to x; is a first-order partial derivative of the surface with respect to y;

[0231] ;

[0232] ;

[0233] ;

[0234] ;

[0235] ;

[0236] k are the principal curvatures k1 and k2 respectively; , , is a second-order partial derivative of the surface Z with respect to x and y.

[0237] S503, when the curvature k of the vehicle-mounted laser point cloud and the inclined model point cloud data is greater than or equal to a preset first curvature threshold, it is determined that the point is located in an edge region, and is classified into a vehicle-mounted laser point cloud potential region and an inclined model point cloud potential region; when the curvature k of the vehicle-mounted laser point cloud and the inclined model point cloud data is less than the preset curvature threshold, it is determined that the point is located in a non-edge region.

[0238] S504, according to a preset neighborhood search radius range, the curvature k of the neighborhood points of each point in the potential region is checked, the number of neighborhood points with curvature greater than the curvature threshold in the neighborhood is judged, and when the ratio of the number of neighborhood points with curvature greater than the first curvature threshold to the total number of neighborhood points is greater than or equal to a preset edge point critical proportion, the current point is determined to be an edge feature point, otherwise the current point is a non-edge feature point. The curvature threshold can be 0.3-0.8, and the edge point critical proportion can be 0.4-0.6.

[0239] S505, the edge feature points obtained through the spatial position relationship and the geometric continuity are connected to form different edge line elements.

[0240] S506, the first principal curvature k1, the second principal curvature k2 and the included angle of the direction vector of the principal curvature of each point in the edge point are calculated, when the included angle is greater than or equal to a preset angle threshold and the first principal curvature and the second principal curvature of the point are both greater than a second curvature threshold, the current point is determined to be a corner point feature point, otherwise the current point is a non-corner point feature point. The angle threshold can be 70-90 degrees.

[0241] In the embodiment of the application, the fusion of the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner point feature points in step S5 to obtain the fused point cloud, comprising:

[0242] S51, the fusion weights of each point cloud in the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud are calculated, the point cloud number with fusion weight greater than a preset weight threshold is named as a first character, the point cloud number with fusion weight less than or equal to the preset weight threshold is named as a second character, the point cloud data with the number named as the first character is taken as target point cloud data, and the point cloud data with the number named as the second character is taken as source point cloud data.

[0243] In the embodiment of the application, the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud in different formats are converted into a unified format, the point cloud fusion data is set to be in las format; the fusion weight is calculated according to the data density and the accuracy level, the laser point cloud plane accuracy coefficient A1, the height accuracy coefficient B1 and the point density coefficient C1 are set, A1+B1+C1=1; the vehicle-mounted laser point cloud point density ρ 11 , the plane accuracy s 11 and the height accuracy h 11 are obtained, and the laser point cloud point density ρ21 , plane accuracy s 21 , and height accuracy h 21。

[0244] ;

[0245] .

[0246] Meanwhile, according to the fusion weight, the vehicle-mounted laser point cloud and the inverse calculation laser point cloud of the tilt model in the multi-source data are numbered. For a large value of the fusion weight, the point cloud number is named as the first character, which can be selected as 1. Conversely, for a small value of the fusion weight, the point cloud number is named as the second character, which can be selected as 2.

[0247] S52, edge line elements and corner feature points of a common overlapping area of the vehicle-mounted laser point cloud and the three-dimensional model point cloud are acquired, and the start point coordinates, end point coordinates and length of each edge line element are acquired; the edge matching relationship of the two kinds of point clouds is judged according to the start point distance, end point distance and length of the edge line elements of the vehicle-mounted laser point cloud and the three-dimensional model point cloud, and two edge line elements of the two kinds of point clouds matched in the edge are formed into a corresponding point pair set of the vehicle-mounted laser point cloud and the three-dimensional model point cloud.

[0248] In the embodiment, the implementation of the edge matching relationship calculation includes: the matching weight is calculated through the start point distance, end point distance and length of the edge line elements of the two kinds of point clouds, and the matching relationship of different point clouds is determined. When the matching weight is greater than or equal to the weight threshold value, two edge lines of different point clouds can be matched and form a corresponding point pair set; when the matching weight is less than the weight threshold value, two edge lines of different point clouds are irrelevant.

[0249] ;

[0250] ;

[0251] ;

[0252] (x 1O , y 1O , z 1O ) is the start point coordinates of the vehicle-mounted laser point cloud;

[0253] (x 2O , y 2O , z 2O ) is the end point coordinates of the vehicle-mounted laser point cloud;

[0254] L1 is the length of the vehicle-mounted laser point cloud line element;

[0255] (x 3O , y 3O , z 3O) is the start coordinate of the inclined three-dimensional model point cloud;

[0256] (x 4O , y 4O , z 4O ) is the end coordinate of the inclined three-dimensional model point cloud;

[0257] L2 is the length of the inclined model point cloud line element;

[0258] S q is the matching weight of the vehicle-mounted laser point cloud and the three-dimensional model point cloud.

[0259] S53, according to the corresponding point pair set of the vehicle-mounted laser point cloud and the three-dimensional model point cloud, the conversion parameters for converting the source point cloud into the target point cloud are calculated, and the coordinate conversion is performed on the source point cloud data according to the conversion parameters, so that the optimized source point cloud data is obtained.

[0260] In this embodiment, through the corresponding point pair set (P v , Q s ) of the vehicle-mounted laser point cloud and the inclined three-dimensional point cloud, the equation is established, wherein is a point in the target point cloud, is a point in the source point cloud. At the same time, the objective function is constructed, and the partial derivatives of the objective function with respect to R S and t S are solved, and the partial derivatives are set to zero to solve the unknown coefficients of the rotation matrix R S and the translation vector t S . At the same time, the centroid of the source point cloud and the centroid of the target point cloud are obtained, and the equation of the translation vector t S is obtained.

[0261] ;

[0262] ;

[0263] ;

[0264] ;

[0265] wherein t s is a translation vector; is the centroid of the target point cloud; is the centroid of the source point cloud; is the objective function. is the number of point cloud pairs in the corresponding point pair set of the vehicle-mounted laser point cloud and the inclined three-dimensional point cloud.

[0266] The translation vector equation is brought into the objective function, a function about R s is obtained, and the optimal rotation matrix Rs , according to the equation of translation vector t s , the translation vector t s is obtained.

[0267] ;

[0268] ;

[0269] ;

[0270] ;

[0271] ;

[0272] wherein, is the coordinate of the target point cloud after translation; is the coordinate of the source point cloud after translation; is the covariance matrix; and is an orthogonal matrix; is a diagonal matrix.

[0273] The obtained rotation matrix R s and translation vector t s are applied to each point in the source point cloud Q s , a coordinate transformation is performed to obtain the optimized point cloud , and the coordinate transformation formula is , wherein is the point in the optimized source point cloud .

[0274] The optimized source point cloud data is taken as new source point cloud data, and the above steps S52-S53 are repeatedly executed until the iteration termination condition is met, the precision-optimized source point cloud data is obtained, and the iteration termination condition is that the change of each conversion parameter between two iterations is less than the corresponding preset threshold value or the coordinate difference of the point cloud before and after transformation is less than the preset error value;

[0275] S54, the precision-optimized source point cloud data is fused with the target point cloud data to obtain a fused point cloud.

[0276] In this embodiment, the transformed source point cloud data is fused with the target point cloud to obtain a fused point cloud. In the overlapping area of the seawall, the seawalls are overlapped with each other; in the non-overlapping seawall area, according to the adjustment parameter at the boundary of the overlapping area, the source point cloud is linearly changed within the buffer radius on both sides of the seawall to obtain a new fused point cloud.

[0277] ;

[0278] ;

[0279] ;

[0280] 、 、 , is the adjustment value in the xyz direction at different distances within the buffer zone; 、 、 is the adjustment parameter in the xyz direction at the boundary of the overlap area; is the buffer zone radius; d0 is the distance to the boundary of the overlap area.

[0281] The present application accurately captures the seawall ridge line feature and edge point feature through curvature driving feature extraction and principal curvature direction analysis, which is 35% higher than the traditional gradient method. The elastic fusion mechanism combines the rigid transformation of the overlap area and the linear gradual change of the non-overlap area to reduce the elevation mutation at the fusion joint. The point cloud is automatically weighted, dynamically balancing the advantages of high-precision laser point cloud and high-density inclined point cloud, and the fusion point cloud plane precision reaches 0.05m, and the density is >200 points / m 2 .

[0282] In the embodiment of the present application, the step S6 of filtering the fusion point cloud to screen out the seawall terrain point cloud data includes the following steps not shown in the drawings:

[0283] S61, point cloud initial classification: block the fusion point cloud and select seed points in each point cloud block, construct an initial triangular net according to the seed points of all point cloud blocks, and classify the seed points into the ground point category; traverse the fusion point cloud, calculate the distance d1 of each point to the nearest initial triangle and the included angle α formed between the nearest initial triangle vertex, judge the relationship between α and the corresponding iteration angle threshold T α and the relationship between d1 and the corresponding iteration distance threshold T d , if d1 is less than or equal to the iteration distance threshold T d and α is less than or equal to the iteration angle threshold T α , the current point is inserted into the initial triangular net and classified into the ground point category, and the next point is continued to be traversed; if d1 is greater than the iteration distance threshold T d or α is greater than the iteration angle threshold T α , the current point is directly skipped, and the next point is continued to be traversed, and the finally obtained point cloud classified into the ground point category is taken as the coarse ground point data.

[0284] In the embodiment, the point cloud initial classification is realized by the following steps, specifically including:

[0285] Point cloud classification: the unified fusion point cloud data is classified as default class;

[0286] Chunking and selecting seed points: first, the fused point cloud survey area is chunked according to the preset building size (usually the diagonal size of the building). Then the lowest point of each block in the survey area is selected as the seed point, and all seed points construct an initial triangular net, and the seed points are classified into ground. The maximum building size is set to 80-200 meters.

[0287] Construction and optimization of triangular net: traverse the default class points in the point cloud, calculate the distance d1 of each point to the nearest triangle and the angle a formed between the nearest triangle vertex, judge the relationship between a and d1 and the set threshold value; when d1 is less than or equal to the set iteration distance threshold T d and a is less than or equal to the set iteration angle threshold T α , the point is inserted into the triangular net and classified into ground, and the traversal continues; when d1 is greater than the set iteration distance threshold T d or a is greater than the set iteration angle threshold T α , it is directly skipped and the traversal continues. Through the traversal of the point cloud, the coarse ground point data is obtained. The distance threshold T d is 0.5-3 meters, and the iteration angle threshold T α is 3-10 degrees.

[0288] ;

[0289] ;

[0290] is the nearest triangular net placed in the plane equation;

[0291] is the coordinate of the default point p in the fused point cloud;

[0292] d1 is the distance d1 of the default point p to the nearest triangle;

[0293] ;

[0294] ;

[0295] ;

[0296] is the projection P' coordinate of P on the triangular net plane;

[0297] is the nearest triangular net vertex;

[0298] Alpha is the included angle between the p point in the default and the nearest triangular vertex.

[0299] S62, the coarse ground point data is filtered for noise points to obtain optimized ground point data with noise points removed. Specifically, according to the fused point cloud data, a seawall surface range line is finely drawn. A distance threshold method is used to remove noise (first step of filtering operation), including: traversing the coarse ground point data, calculating the first elevation difference of each point in the coarse ground point data and the first elevation of each point within a first preset search radius range of the current point; when the first elevation difference is greater than a first preset height difference threshold H1, the current point is classified into a noise category, and when the first elevation difference is less than or equal to the first preset height difference threshold H1, the ground point category is maintained, and the next point is continued to be traversed, and finally the ground point data with coarse noise points removed is obtained. The first search radius is 0.5 meters-2 meters, and the first height difference threshold is 0.5 meters-3 meters. Then a fitting plane distance method is used to remove noise (second step of filtering operation), including: traversing the ground point data with coarse noise points removed, obtaining all neighborhood points within a second preset search radius range of the current point for each point therein, constructing a fitting plane equation according to the neighborhood point coordinates of the current point, calculating the first distance h1 of the current point to the fitting plane, the second distance h2 of the neighborhood points to the fitting plane, and the distance average h3 of the neighborhood points to the fitting plane, if the first distance h1 is less than or equal to the distance average h3*limit and the first distance h1 is less than or equal to the second distance h2, the ground point category of the current point is maintained, and the next point is continued to be traversed; if the first distance h1 is greater than the distance average h3*limit or the first distance h1 is greater than the second distance h2, the current point is classified into a noise category, and the next point is continued to be traversed, and finally the optimized ground point data with secondary noise points removed is obtained; wherein limit is a preset limit coefficient. The second search radius is 0.5 meters-2 meters, the second height threshold is 0.10 meters-1 meter, and the limit coefficient is 2-4. h1: the vertical distance from the current point to the fitting plane. h2: the set of vertical distances from all neighborhood points to the same fitting plane. h3: the average of the vertical distances from all neighborhood points to the fitting plane, representing the overall fluctuation degree of the local area.

[0300] S63, the coarse ground point data with noise points removed in the optimized ground point data is reduced in iterative distance threshold T d and iterative angle threshold T α , and the point cloud classification step is performed again to obtain fused point cloud data containing new ground. The distance threshold T d is 0.5 meters-1 meter, and the iterative angle threshold T α is 2 degrees-5 degrees.

[0301] The seawall specific range line is drawn again, and a second iterative distance threshold T d2 and a second iterative angle threshold Tα2 , the block size, i.e., the maximum building size L, and performing a point cloud initial classification step on the point cloud data within the seawall range to classify seawall top surface point data and two side wall surface point data and classify them as ground point data;

[0302] S64, obtaining a normal vector of each point in the ground point data: setting a neighborhood point set for each point in the ground point data, calculating a geometric center of the neighborhood point set, constructing a covariance matrix of the neighborhood point set according to the geometric center of the neighborhood point set, and performing eigenvalue decomposition on the covariance matrix to obtain eigenvalues and eigenvectors of the covariance matrix, and selecting an eigenvector corresponding to the smallest eigenvalue as the normal vector of the current point;

[0303] In this embodiment, the step of obtaining the normal vector of each point in the seawall ground point data comprises:

[0304] Determine the neighborhood point cloud. For the point p i containing ground in the fusion point cloud, set the neighborhood number k i , calculate the distance between other points q i and p i , and take the k i closest points as the neighborhood point set N i .

[0305] ;

[0306] Where (x p , y p , z p ) is the coordinate of p i ; (x q , y q , z q ) is the coordinate of q i . is the distance between two points p i and q i .

[0307] Calculate the centroid. Calculate the geometric center of the neighborhood point set N i , i.e., the centroid C O .

[0308] ;

[0309] ;

[0310] ;

[0311] ;

[0312] Where, is the neighborhood point set N inumber of midpoints; p j is the jth point in the set of neighborhood midpoints; is p j 's coordinate.

[0313] Compute the covariance matrix. Construct the covariance matrix Cov i of the set of neighborhood points N

[0314] ;

[0315] ;

[0316] ;

[0317] ;

[0318] ;

[0319] ;

[0320] ;

[0321] where, is the variance in the x-dimension; is the variance in the y-dimension; is the variance in the z-dimension; is the covariance between the x and y dimensions; is the covariance between the x and z dimensions; is the covariance between the y and z dimensions.

[0322] Eigen decomposition. Perform an eigen decomposition of the covariance matrix Cov to obtain eigenvalues (including , , ) and corresponding eigenvectors , , .

[0323] ;

[0324] ;

[0325] ;

[0326] where, is the eigenvalue; is the identity matrix.

[0327] Normal vector determination. The normal vector n i to point P corresponding eigenvector . Traverse the fusion point cloud containing ground, and obtain the normal vector of each point.

[0328] S65, triangulate the points in the ground point data, obtain the three corner points of the triangle in which each point falls, construct the plane equation of the triangle, and calculate the plane distance from the current point to the triangle based on the plane equation.

[0329] In this embodiment, the ground points in the laser point cloud of the previous step are triangulated to obtain triangulation data; the plane equation of the triangle in which the points in the laser point cloud fall is constructed; and the plane distance from the point to the triangle is calculated. Specifically, it includes:

[0330] (1) Ground TIN network construction. Triangulate the ground points in the laser point cloud of the previous step to obtain triangulation data.

[0331] (2) Plane equation establishment. Obtain the three corner points of the triangle in which the point P i (x7, y7, z7) in the laser point cloud falls, and calculate the plane equation coefficients of the triangle.

[0332] ;

[0333] (3) Distance calculation. Calculate the plane distance d i from the point P T to the triangle.

[0334] .

[0335] S66, filter the points in the ground point data whose corresponding plane distance is greater than a preset seawall height distance threshold as noise points, and obtain the point cloud data near the seawall and the ground.

[0336] In this embodiment, a seawall height distance threshold is set, the points in the fusion point cloud are traversed, and the relationship between the distance value of the point and the seawall height distance threshold is judged; when the distance value of the point is greater than the seawall height distance threshold, the point is classified as a low vegetation class; when the distance value of the point is less than or equal to the seawall height distance threshold, the point is still the original class point, and the point cloud data near the seawall and the ground is obtained.

[0337] S67, traverse the point cloud data near the seawall and the ground, select the points located on the side vertical surface as seed points according to the normal vector of the points, and constantly add the adjacent points that meet the normal vector similarity and distance constraint conditions to the growth region according to a preset growth criterion until there is no point meeting the conditions, and the growth criterion of the growth region is that the normal vector of the seed point is similar to the normal vector of the point in the growth region and the distance between the seed point and the point in the growth region is less than a preset distance threshold.

[0338] In this embodiment, the growth criterion is determined according to the direction amount and distance of the point; the seed point is selected, and the adjacent points that meet the normal vector similarity and distance condition at the same time are continuously added to the growth region according to the growth criterion. Specifically, it includes:

[0339] Growth criterion determination. The angle between the normal vector of the seed point and the normal vector of the point in the grown region is calculated, and when the cosine value is greater than or equal to the normal vector threshold, the normal vectors of the two points are similar. At the same time, the distance between the seed point and the point in the grown region is calculated to be less than the distance threshold, which meets the distance condition.

[0340] Seed point P ii (x ii , y ii , z ii ) The normal vector of is , the normal vector of the point P jj (x jj , y jj , z jj ) in the growth region is , and the cosine value of the angle between the two normal vectors is calculated.

[0341] ;

[0342] ;

[0343] Some points located on the side vertical surface are selected as seed points. The normal vector of the side vertical surface point (x n , y n , z n ).

[0344] ;

[0345] Starting from the seed point, the adjacent points that meet the normal vector similarity and distance condition at the same time are continuously added to the growth region according to the growth criterion, and the point cloud data of the dike fusion processing is traversed until there is no more point that meets the condition to be added.

[0346] S68, the area of the growth region obtained by growth is counted, and when the area of the region is greater than or equal to the preset area threshold, the growth region is the dike side vertical surface point cloud, the point cloud in the growth region is classified as the side vertical surface category, and the point cloud classified into the side vertical surface category is taken as the side vertical surface point data;

[0347] In this embodiment, the area threshold is set, the area of the region point cloud obtained by growth is counted, when the area of the region is less than the area threshold, the region point cloud type is not processed; when the classification processing is performed, when the area of the region is greater than or equal to the area threshold, the region is the dike side vertical surface point cloud, and the region point cloud is classified as the side vertical surface wall category.

[0348] S69. Use the side elevation point data and the ground point data as the initial seawall terrain point cloud data related to the seawall terrain.

[0349] Furthermore, the embodiment of the present invention further includes: performing accuracy detection on the initial seawall terrain point cloud data. The accuracy detection of the initial seawall terrain point cloud data specifically includes: dividing the polygon area corresponding to the seawall range according to the seawall length to form grid blocks, loading lidar point cloud data, and obtaining the point cloud data corresponding to the grid blocks; loading the oblique photography model data and the point cloud data, and judging the category of the point cloud data corresponding to each grid block one by one according to the texture information and spatial structure information provided by the oblique photography model data, and screening out the point cloud of the category related to the seawall terrain from the point clouds of each sub-block to obtain the seawall terrain point cloud data. Specifically, screen out the point cloud of the category related to the "ji" - shaped terrain of the seawall from the point clouds of each sub-block.

[0350] In the embodiment of the present invention, based on the seawall range line, buffer analysis is performed to generate buffer polygons with a width of 50 meters on both sides of the seawall range line, and then the polygon area is further divided according to the seawall length of 200 meters to form regular small grid blocks. In professional point cloud processing software, combined with the generated regular small grid blocks, batch - load lidar point cloud data to obtain the corresponding small grid block point cloud data; load the oblique photography model data and the point cloud data, classify the point cloud data of each sub - block one by one, and judge the category of the point cloud by means of the texture information and spatial structure information provided by the oblique model. Check the point cloud on the side of the main body of the seawall and classify it as the "wall" category. Check the point cloud of the surrounding vegetation and classify it as the "low vegetation" category. Check the point cloud of the exposed ground and classify it as the "ground" category. Screen out the point cloud of the category related to the seawall terrain from the point clouds of all sub - blocks to obtain the seawall terrain point cloud data.

[0351] The terrain adaptive classification proposed by the present invention dynamically adjusts through iterative distance thresholds and angle thresholds. The correct classification rate of ground points is greater than 92%, and the misjudgment rate in the steep slope area drops by 60%. The protection of the facade structure is realized. The seawall side wall is accurately extracted through the normal vector - constrained region growth algorithm (cosine similarity > 0.95), and the integrity retention rate of the facade structure is 98%. Based on the collaborative quality inspection with the texture information provided by the oblique model, grid - based sub - division (200×50m) combined with the texture verification of the oblique model, the correction efficiency is increased by 80%, and the purity of the terrain point cloud is > 99%.

[0352] In the embodiment of the present invention, the step of constructing a resolution - adaptive digital elevation model for different regions of the seawall dam with structural features according to the point cloud average density of the seawall terrain point cloud data and the planar width of the seawall dam in step S7 includes steps not shown in the following drawings:

[0353] S71, determining the standard grid size of the digital elevation model (DEM) according to the point cloud average density of the seawall terrain point cloud data and the planar width of the seawall dam, taking the crest and slope of the seawall dam as a high-resolution area, taking the base of the seawall dam as a low-resolution area, and reducing the grid size in the high-resolution area. The crest is a high platform, the slope is a steep slope connecting the high platform, and the base is a gentle slope or platform transition zone at the bottom.

[0354] In this embodiment, the grid size, i.e., the resolution, of the DEM is determined according to the point cloud average density of the seawall terrain point cloud data and the planar width of the seawall dam. The local density of the point cloud is calculated, and for each point P i , a search sphere radius R with a variable radius is set as the center i , the number of points n in the sphere is counted i , and the local density is calculated.

[0355] ;

[0356] ;

[0357] According to the planar width of the seawall dam and the structural characteristics, the crest and slope of the seawall dam are taken as a high-resolution area, the base of the seawall dam is taken as a low-resolution area, the grid size is reduced in the high-resolution area, and the weight coefficient of the high-resolution area is w (w>1). The grid size S local =S global / w, where S global is the grid size preliminarily estimated according to the global average density.

[0358] S72, for the high-resolution area and the low-resolution area of the seawall dam, a triangular mesh is constructed by using different neighborhood radii, and the elevation values of the grid nodes in the corresponding area are calculated by linear interpolation according to the elevation of the vertices of the triangular mesh. For the overlapping area of the high-resolution area and the low-resolution area of the seawall, the linear interpolation results of the overlapping area are weighted and fused to obtain the elevation value of the overlapping area.

[0359] In this embodiment, a smaller search radius is set in the high-resolution area of the seawall, and a larger search radius is set in the low-resolution area of the seawall. A Delaunay triangular mesh is constructed in its neighborhood, and then the elevation of the grid node is calculated by linear interpolation according to the elevation of the vertices of the triangular mesh.

[0360] The overlapping range of the high-resolution area and the low-resolution area of the seawall is determined by a spatial analysis algorithm, and the interpolation results of the overlapping area are weighted and fused. For each grid node in the overlapping area, the grid size S high of the high-resolution area and the grid size S low of the low-resolution area are obtained, and the weight is calculated.

[0361] ;

[0362] ;

[0363] The elevation value z in the overlapping area is obtained by weighting the elevation value z interpolated in the high-resolution area and the elevation value z interpolated in the low-resolution area. high and the elevation value z interpolated in the low-resolution area. low The elevation value z in the overlapping area is obtained by weighting the elevation value z interpolated in the high-resolution area and the elevation value z interpolated in the low-resolution area.

[0364] ;

[0365] Wherein, S high and S low are the grid sizes of the high-resolution area and the low-resolution area respectively.

[0366] S73, constructing a digital elevation model in the range of the sea dike according to the elevation values of the high-resolution area, the low-resolution area and the overlapping area of the sea dike dam.

[0367] The method for extracting the sea dike terrain by fusing the vehicle-mounted laser point cloud and the tilt model provided by the application is a full-process solution for the problems of low registration accuracy, fusion edge jump and terrain classification distortion caused by traditional single data source (vehicle-mounted laser radar has blind area of steep slope shielding, hollow area of weak texture in tilt photography). The method first realizes the millimeter-level reference unification of linear engineering by dynamic encryption control network (more than 4 points in complex terrain area, Lmean<2000m), and makes the plane residual error less than 0.02m and the height fitting mean error less than 2cm through seven-parameter conversion (Boussinesq-Wolf model) combined with quadratic surface fitting of height anomaly. In the data acquisition link, the vehicle-mounted point cloud is denoised by distance threshold and standard deviation (noise elimination rate>95%) and the corner points are automatically extracted based on intensity threshold and structure tensor (weak texture registration success rate 92%), and the tilt model is improved in geometric fidelity through octree dynamic segmentation and Delaunay triangulation optimization (narrow triangle elimination rate 100%). The method realizes the automatic weighted fusion of point cloud through the innovation of curvature-driven feature extraction (main curvature analysis makes the ridge line recognition rate increased by 35%) and elastic fusion mechanism (rigid registration in overlapping area and linear gradual change in non-overlapping area, edge height difference<0.05m). The method ensures the purity of sea dike terrain point cloud>99% through terrain adaptive iteration and normal vector constraint region growing (facade structure retention rate 98%) combined with grid human-computer collaborative quality inspection (200x50m block+tilt model texture verification). The resolution adaptive DEM finally constructed can accurately represent the structural characteristics of the sea dike, provide sub-meter data support for coastal safety monitoring and disaster assessment, and fill the technical gap of steep terrain measurement in the current specification.

[0368] For the method embodiments, the methods will be described in terms of a collection of actions to be performed by a subject, but it will be understood that the methods described herein are not limited to the order of acts described. The methods can be performed in other orders or concurrently, and elements of one method can be incorporated into other methods. Further, the methods described herein are not limited to the specific examples described and should be understood to include any method that accomplishes the same, similar or equivalent results.

[0369] Another embodiment of the present application also provides a system for extracting sea embankment terrain by fusing vehicle-mounted laser point cloud and tilt model, which comprises functional modules for implementing the method for extracting sea embankment terrain by fusing vehicle-mounted laser point cloud and tilt model according to any one of the above. Figure 2 The structural block diagram of the system for extracting sea embankment terrain by fusing vehicle-mounted laser point cloud and tilt model according to another embodiment of the present application is schematically shown. Referring to Figure 2 The system for extracting sea embankment terrain by fusing vehicle-mounted laser point cloud and tilt model according to the embodiment specifically comprises a sea embankment range drawing module 201, a measurement range dividing module 202, a vehicle-mounted laser point cloud processing module 203, a tilt photography data processing module 204, a point cloud fusion module 205, a fused point cloud filtering module 206 and a sea embankment terrain extraction module 207, wherein:

[0370] The sea embankment range drawing module 201 is configured to extend a buffer zone to both sides of the sea embankment center line in the satellite image according to the sea embankment vector diagram to form a sea embankment range;

[0371] The measurement range dividing module 202 is configured to divide the sea embankment range into a vehicle-mounted laser radar collection range and a UAV tilt photography measurement range;

[0372] The vehicle-mounted laser point cloud processing module 203 is configured to acquire the vehicle-mounted laser point cloud collected by the vehicle-mounted laser radar, perform unified spatial reference coordinate conversion and denoising processing on the vehicle-mounted laser point cloud, and perform precision optimization on the denoised vehicle-mounted laser point cloud;

[0373] The tilt photography data processing module 204 is configured to acquire the tilt photography image data and POS data collected by the UAV-mounted five-patch camera, perform unified spatial reference coordinate conversion on the POS data, and generate a tilt photography model and a tilt three-dimensional model point cloud according to the converted POS data and the tilt photography image data;

[0374] The point cloud fusion module 205 is configured to extract edge line elements and corner feature points in the common overlapping area of the tilt three-dimensional model point cloud and the vehicle-mounted laser point cloud, and fuse the tilt three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner feature points to obtain a fused point cloud;

[0375] The fusion point cloud filtering module 206 is configured to filter the fusion point cloud to obtain seawall terrain point cloud data.

[0376] The seawall terrain extraction module 207 is configured to construct a resolution-adaptive digital elevation model according to the point cloud average density of the seawall terrain point cloud data and the planar width of the seawall dam, and take the different regions of the seawall dam as structural features.

[0377] For the system embodiments, the description is relatively simple because the system embodiments are basically similar to the method embodiments, and the relevant parts can be referred to the part of the description of the method embodiments.

[0378] The system embodiments described above are only illustrative, wherein the units described as separate components can or can not be physically separated, and the components shown as units can or can not be physical units, that is, can be located in one place or distributed on multiple network units. Part or all of the modules can be selected to achieve the purpose of the embodiment according to actual needs. Those skilled in the art can understand and implement without creative labor.

[0379] In addition, another embodiment of the present application further provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor; the computer program is executed by the processor to implement the steps of the seawall terrain extraction method of fusing the vehicle-mounted laser point cloud and the tilt model.

[0380] In addition, another embodiment of the present application further provides a computer program product, wherein the computer program product stores a computer program, and the computer program is executed by a processor to implement the steps of the seawall terrain extraction method of fusing the vehicle-mounted laser point cloud and the tilt model.

[0381] Those skilled in the art can understand that although some embodiments herein include certain features but not others included in other embodiments, the combination of features of different embodiments means to be within the scope of the present application and form different embodiments. For example, any one of the claimed embodiments can be used in any combination.

[0382] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, but not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement to part of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for extracting a seawall terrain by fusing a vehicle-mounted laser point cloud and an inclined model, characterized by, The method comprises: According to the seawall vector diagram, a buffer zone is extended to both sides of the seawall center line in the satellite image to form a seawall range; The seawall range is divided into a vehicle-mounted laser radar collection range and an unmanned aerial vehicle oblique photogrammetry range; Vehicle-mounted laser point clouds collected by the vehicle-mounted laser radar are obtained, unified spatial reference coordinate conversion and denoising processing are performed on the vehicle-mounted laser point clouds, and precision optimization is performed on the denoised vehicle-mounted laser point clouds; Oblique photography image data and POS data collected by the unmanned aerial vehicle carrying a five-patch camera are obtained, unified spatial reference coordinate conversion is performed on the POS data, and a tilt photography model and a tilt three-dimensional model point cloud are generated according to the converted POS data and the oblique photography image data; Edge line elements and corner feature points in the common overlapping area of the tilt three-dimensional model point cloud and the vehicle-mounted laser point cloud are extracted, and the tilt three-dimensional model point cloud and the vehicle-mounted laser point cloud are fused based on the edge line elements and the corner feature points to obtain fused point clouds; The fused point clouds are filtered to obtain seawall terrain point cloud data; According to the point cloud average density of the seawall terrain point cloud data and the planar width of the seawall dam, a resolution-adaptive digital elevation model is constructed for different regions of the seawall dam as structural features.

2. The method of claim 1, wherein, The precision optimization of the denoised vehicle-mounted laser point cloud comprises: S31, taking the denoised vehicle-mounted laser point cloud data as original vehicle-mounted point cloud; S32, taking the coordinates of the image control points arranged in the seawall range in advance as the center, and obtaining an image control point set within a preset intensity threshold range according to a preset radius threshold range in the original vehicle-mounted point cloud; S33, identifying the corner points by using an adaptive gray threshold and a tension algorithm to obtain the coordinates of the image control point corners; S34, calculating the point cloud image control point center point according to the coordinates of the image control point corners, selecting a target image control point from the image control point set measured in the field which matches the distance of the point cloud image control point center point, and constructing a corresponding point pair between the point cloud image control point center point and the target image control point; S35, establishing a corresponding relationship equation according to each corresponding point pair, and solving the transformation parameters between the point cloud image control point center point and the target image control point according to the corresponding relationship equation; S36, performing coordinate transformation on each point in the original vehicle-mounted point cloud according to the transformation parameters to obtain transformed point clouds; Taking the transformed point clouds as new original vehicle-mounted point clouds, steps S32 to S36 are repeatedly executed until an iteration termination condition is met, to obtain precision-optimized vehicle-mounted laser radar point cloud data, and the iteration termination condition is that the changes of each transformation parameter between two iterations are all less than a corresponding preset threshold value or the coordinate difference of the point clouds before and after transformation is less than a preset error value.

3. The method of claim 2, wherein, In the original vehicle-mounted point cloud, the image control point set within the preset intensity threshold range is obtained according to the preset radius threshold range with the image control point coordinates arranged in the seawall range as the center, including: extracting the laser point cloud data Plocal within the preset radius threshold R1 range with the image control point coordinates arranged in the seawall range as the center, calculating the intensity mean value of Plocal and the standard deviation , calculating the minimum intensity threshold I min and the maximum intensity threshold I max according to the intensity mean value and the standard deviation ; traversing the laser point cloud data within the preset radius threshold R2 range with the image control point coordinates as the center, if the intensity value of the current point is less than or equal to I max and greater than or equal to I min , the current point is added to the preset point set, and finally the image control point set Q is formed, R2 Adaptive gray threshold and tension algorithm is used to identify the corner, and the image control point corner coordinates are obtained, including: two-dimensional projection is performed on the image control point set Q, for each two-dimensional point after projection, the adaptive window is dynamically adjusted according to the point cloud density and local terrain complexity, the gray gradient of each two-dimensional point in x direction and y direction is calculated and and structural tensor M z , the structural tensor M z is used to construct the corner response function R z , and whether the current two-dimensional point is an image control point corner is judged based on the corner response function R z , and finally six two-dimensional image control point corners are determined.

4. The method of claim 1, wherein, According to the converted POS data and the oblique photography image data, the tilt photography model and the tilt three-dimensional model point cloud are generated, which comprises: Performing corner and edge feature extraction on the oblique photography image data, using a SIFT algorithm to describe the key points of the images subjected to the corner and edge feature extraction in different spatial scales to obtain feature descriptors of the images; According to the feature matching results of the corner and edge features and the ground control point information, using a bundle adjustment algorithm to determine the exterior orientation elements of the aerial triangulation of the oblique photography image data. According to the exterior orientation elements of aerial triangulation, dense matching is performed on the oblique photography image data, a dense oblique three-dimensional model point cloud is generated according to a preset laser point cloud density threshold, and a mesh three-dimensional surface model is constructed by using the oblique three-dimensional model point cloud; texture information of the oblique photography image data is mapped to the mesh three-dimensional surface model through a texture mapping technology, and an oblique photography model with real texture is generated.

5. The method of claim 4, wherein, The mesh three-dimensional surface model is constructed by using the oblique three-dimensional model point cloud, the texture information of the oblique photography image data is mapped to the mesh three-dimensional surface model through the texture mapping technology, and the oblique photography model with real texture is generated, including: An octree structure is used to construct a spatial index for the oblique three-dimensional model point cloud; Adjacent point searching is performed according to the spatial index, and a triangular mesh is constructed based on the searched candidate adjacent points; the triangular mesh is optimized based on a long and narrow triangle removal optimization principle and an edge exchange optimization principle; The texture information of the oblique photography image data is mapped to the optimized triangular mesh through the texture mapping technology, and the oblique photography model with real texture is obtained.

6. The method of claim 1, wherein, Edge line elements and corner feature points in a common overlapping area of the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud are extracted, including: A set of neighborhood points of each point in the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud is obtained, a covariance matrix of the neighborhood points is constructed, and the curvature k of the neighborhood points is calculated according to the covariance matrix; For each point and neighborhood point in the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud, surface fitting is performed to obtain a local surface equation, and a first principal curvature k1, a second principal curvature k2, and respective principal curvature direction vectors are calculated according to the local surface equation ; When the curvature k of the vehicle-mounted laser point cloud and the oblique model point cloud data is greater than or equal to a preset curvature threshold, it is determined that the point is located in an edge region and belongs to a vehicle-mounted laser point cloud potential region and an oblique model point cloud potential region; when the curvature k of the vehicle-mounted laser point cloud and the oblique model point cloud data is less than the preset curvature threshold, it is determined that the point is located in a non-edge region; According to a preset neighborhood search radius range, the curvature k of the neighborhood points of each point in the potential region is checked, the number of neighborhood points whose curvature k is greater than a first curvature threshold is judged, and when the ratio of the number of neighborhood points whose curvature k is greater than the first curvature threshold to the total number of neighborhood points is greater than or equal to a preset edge point critical proportion, it is determined that the current point is an edge feature point; otherwise, the current point is a non-edge feature point; The edge feature points obtained through the spatial position relationship and the geometric continuity are connected to form different edge line elements; The first principal curvature k1, the second principal curvature k2 and the included angle of the direction vectors of the principal curvatures of each point in the edge points are calculated, and when the included angle is greater than or equal to a preset angle threshold and the first principal curvature and the second principal curvature of the point are both greater than a second curvature threshold, it is determined that the current point is a corner feature point; otherwise, the current point is a non-corner feature point.

7. The method of claim 1, wherein, The oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud are fused based on the edge line elements and the corner feature points to obtain a fused point cloud, including: S51, calculate the fusion weight of each point cloud in the inclined three-dimensional model point cloud and the vehicle-mounted laser point cloud, name the point cloud whose fusion weight is greater than a preset weight threshold as a first character, name the point cloud whose fusion weight is less than or equal to the preset weight threshold as a second character, take the point cloud data numbered as the first character as target point cloud data, and take the point cloud data numbered as the second character as source point cloud data; S52, obtain the edge line elements and corner feature points of the common overlapping area of the vehicle-mounted laser point cloud and the three-dimensional model point cloud, and obtain the start point coordinates, end point coordinates and length of each edge line element; determine the edge matching relationship of the two kinds of point clouds according to the start point distance, end point distance and length of the edge line elements of the vehicle-mounted laser point cloud and the three-dimensional model point cloud, and form a corresponding point pair set of the vehicle-mounted laser point cloud and the three-dimensional model point cloud from the two edge line elements of the two kinds of point clouds matched by edges; S53, calculate the conversion parameters for converting the source point cloud into the target point cloud according to the corresponding point pair set of the vehicle-mounted laser point cloud and the three-dimensional model point cloud, and perform coordinate conversion on the source point cloud data according to the conversion parameters to obtain optimized source point cloud data; take the optimized source point cloud data as new source point cloud data, and repeatedly execute the above steps S52~S53 until the iteration termination condition is met, to obtain precision-optimized source point cloud data, and the iteration termination condition is that the change of each conversion parameter between two iterations is less than the corresponding preset threshold value or the coordinate difference of the point cloud before and after transformation is less than a preset error value; S54, fuse the precision-optimized source point cloud data and the target point cloud data to obtain fused point cloud data.

8. The method of claim 1, wherein, Filtering the fused point cloud to obtain dike terrain point cloud data, comprising: Point cloud initial classification: the fusion point cloud is blocked and a seed point is selected in each point cloud block, an initial triangular net is constructed according to the seed points of all point cloud blocks, and the seed points are classified into the ground point category; the fusion point cloud is traversed, the distance d1 of each point to the nearest initial triangle and the included angle a between the nearest initial triangle vertex are calculated, the relationship between a and the corresponding iteration angle threshold T α and the relationship between d1 and the corresponding iteration distance threshold T d are judged, if d1 is less than or equal to the iteration distance threshold T d and a is less than or equal to the iteration angle threshold T α , the current point is inserted into the initial triangular net and classified into the ground point category, and the next point is continued to be traversed; if d1 is greater than the iteration distance threshold T d or a is greater than the iteration angle threshold T α , the current point is directly skipped, the next point is continued to be traversed, and the finally obtained point cloud classified into the ground point category is taken as the coarse ground point data; performing noise point filtering on the coarse ground point data to obtain optimized ground point data free of noise points; Lowering the iterative distance threshold T d and the iterative angle threshold T α and re-performing the point cloud initial classification step to obtain new ground point data Drawing a specific range line of the seawall, setting a second iteration distance threshold T matching the size of the seawall d2 , a second iteration angle threshold T α2 , a block size, and performing a point cloud preliminary classification step on the point cloud data in the seawall range to classify the seawall top surface point data and the two side wall point data into ground point data; obtaining the normal vector of each point in the ground point data: setting a neighborhood point set for each point in the ground point data, calculating the geometric center of the neighborhood point set, constructing a covariance matrix of the neighborhood point set according to the geometric center of the neighborhood point set, and performing eigenvalue decomposition on the covariance matrix to obtain the eigenvalues and eigenvectors of the covariance matrix, and selecting the eigenvector corresponding to the minimum eigenvalue as the normal vector of the current point; triangulation is constructed for the points in the ground point data, three corner points of the triangle in which each point falls are obtained, and a plane equation of the triangle is constructed, and the plane distance of the current point to the triangle is calculated based on the plane equation; filtering the points in the ground point data corresponding to the plane distance greater than a preset dike height distance threshold as noise points to obtain dike and ground near point cloud data; traverse the dike and ground near point cloud data, select the points located on the side vertical surface as seed points according to the normal vector of the points, and continuously add the adjacent points meeting the normal vector similarity and distance constraint conditions to the growing region according to a preset growing rule until there is no point meeting the conditions, and the growing rule of the growing region is that the normal vector of the seed point is similar to the normal vector of the point in the growing region and the distance between the seed point and the point in the growing region is less than a preset distance threshold; The area of the growth region point cloud obtained by growth is counted, and when the area of the region is greater than or equal to a preset area threshold, the growth region is a seawall side elevation point cloud, the point cloud in the growth region is classified into a side elevation category, and the point cloud classified into the side elevation category is taken as side elevation point data; The side elevation point data and the ground point data are taken as initial seawall terrain point cloud data related to the seawall terrain.

9. The method of claim 8, wherein, The method further comprises: The accuracy of the initial seawall terrain point cloud data is detected, including: dividing the polygon region corresponding to the seawall range according to the length of the seawall to form a grid block, loading the laser radar point cloud data to obtain the point cloud data corresponding to the grid block, loading the oblique photography model data and the point cloud data, and judging the category of the point cloud data corresponding to the grid block according to the texture information and spatial structure information provided by the oblique photography model data, and filtering out the point cloud of the category related to the seawall terrain from the point cloud of each sub-block to obtain the seawall terrain point cloud data.

10. The method of claim 8, wherein, The coarse ground point data is filtered to obtain optimized ground point data without noise, including: The coarse ground point data is traversed, the first elevation difference of each point in the coarse ground point data and each point in the first preset search radius range of the current point is calculated, when the first elevation difference is greater than the first preset elevation difference threshold H1, the current point is classified into a noise category, when the first elevation difference is less than or equal to the first preset elevation difference threshold H1, the ground point category is maintained, and the next point is continuously traversed, and finally the ground point data without coarse noise is obtained; The coarse ground point data is traversed, the first elevation difference of each point in the coarse ground point data and each point in the first preset search radius range of the current point is calculated, when the first elevation difference is greater than the first preset elevation difference threshold H1, the current point is classified into a noise category, when the first elevation difference is less than or equal to the first preset elevation difference threshold H1, the ground point category is maintained, and the next point is continuously traversed, and finally the ground point data without coarse noise is obtained; 11. The method of claim 1, wherein, According to the average density of the seawall terrain point cloud data and the planar width of the seawall dam, a resolution-adaptive digital elevation model is constructed for different regions of the seawall dam with structural features, including: According to the average density of the seawall terrain point cloud data and the planar width of the seawall dam, the standard grid size of the digital elevation model DEM is determined, the top of the seawall dam and the slope are taken as high-resolution regions, and the base of the seawall dam is taken as a low-resolution region, and the grid size is reduced in the high-resolution region; For the high-resolution area and the low-resolution area of the seawall dam, the triangular net is constructed by adopting different neighborhood radii, and the elevation values of the grid nodes in the corresponding area are calculated by linear interpolation according to the elevation of the vertex of the triangular net; for the overlapping area of the high-resolution area and the low-resolution area of the seawall, the linear interpolation results of the overlapping area are fused by weighting to obtain the elevation value of the overlapping area; The digital elevation model in the range of the seawall is constructed according to the elevation values of the high-resolution area, the low-resolution area and the overlapping area of the seawall dam.

12. A system for fusing vehicle-mounted laser point cloud and tilt model for seawall terrain extraction, characterized in that, The system comprises: A seawall range drawing module is configured to extend a buffer area on both sides of the seawall center line in the satellite image according to the seawall vector diagram to form a seawall range; A measurement range division module is configured to divide the seawall range into a vehicle-mounted laser radar collection range and an unmanned aerial vehicle oblique photogrammetry range; A vehicle-mounted laser point cloud processing module is configured to acquire the vehicle-mounted laser point cloud collected by the vehicle-mounted laser radar, perform unified spatial reference coordinate conversion and denoising processing on the vehicle-mounted laser point cloud, and perform precision optimization on the denoised vehicle-mounted laser point cloud; An oblique photography data processing module is configured to acquire oblique photography image data and POS data collected by the five-patch camera carried by the unmanned aerial vehicle, perform unified spatial reference coordinate conversion on the POS data, and generate an oblique photography model and an oblique three-dimensional model point cloud according to the converted POS data and the oblique photography image data; A point cloud fusion module is configured to extract edge line elements and corner feature points of a common overlapping area in the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud, fuse the oblique three-dimensional model point cloud and the vehicle-mounted laser point cloud based on the edge line elements and the corner feature points to obtain fused point cloud; A fused point cloud filtering module is configured to filter the fused point cloud to obtain seawall terrain point cloud data; A seawall terrain extraction module is configured to construct a resolution-adaptive digital elevation model according to the point cloud average density of the seawall terrain point cloud data and the planar width of the seawall dam, and take the structural features as different areas of the seawall dam.

13. A computer device, comprising: The computer program is stored on the memory and can be run on the processor; the computer program is executed by the processor to realize the steps of the method according to any one of claims 1-11.

14. A computer program product, characterised in that, The computer program product stores a computer program, and the computer program is executed by the processor to realize the steps of the method according to any one of claims 1-11.

Images