An underwater terrain point cloud filtering and three-dimensional modeling method

Abstract

The present application relates to the technical field of image data processing and three-dimensional space modeling, and discloses a kind of underwater terrain point cloud filtering and three-dimensional modeling method, comprising: according to acoustic pulse timing and heading azimuth, projection is generated two-dimensional profile point array by sounding point cloud;Conversion into double-channel geometric feature tensor containing local point density and height gradient, utilize convolution target detection network to predict bounding box and eliminate noise point;According to the elevation mean and standard deviation of shallow point subset, the representative water depth of grid cell is calculated;Based on the strong constraint of non-convex closed topological boundary ring, triangular mesh is constructed, and the external triangle is removed by using the parity of the intersection of the center of mass ray, and a three-dimensional underwater digital elevation model is generated, the present application decouples suspended noise and micro-landform slope, avoids the influence of grid homogenization on characteristic shallow points, and realizes the automatic reconstruction of non-convex boundary continuous three-dimensional surface.

Description

Technical Field

[0001] This invention relates to the field of image data processing and three-dimensional spatial modeling technology, specifically to an underwater terrain point cloud filtering and three-dimensional modeling method. Background Technology

[0002] Three-dimensional digital elevation models of underwater topography are crucial foundational data for channel dredging and maintenance, as well as spatial analysis of underwater engineering projects. Current surveying operations typically involve a surveying platform equipped with an acoustic bathymetry system collecting underwater point cloud sequences along the survey line. After tidal and sound velocity corrections, a three-dimensional discrete point cloud is generated, serving as the basis for generating the digital elevation model and reconstructing three-dimensional surfaces. However, during actual water body exploration, factors such as suspended sediment, plankton, cavitation microbubbles, and multipath refraction of sound waves introduce numerous discrete noise points and elevation jumps into the original bathymetry sequence. Furthermore, in areas with significant undulations, such as channel slopes and deep channels, the actual terrain itself exhibits noticeable elevation changes over short distances. In general filtering algorithms, outliers are usually identified based on statistical elevation thresholds in three-dimensional space. However, this makes it difficult to accurately distinguish between suspended noise points in water bodies and the actual elevation changes on steep slope sections. When the filtering threshold is set too strictly, it is easy to smooth out the actual underwater slope breaks. When the filtering threshold is set too wide, it is easy to retain suspended noise points, which can create false terrain in subsequent modeling. In the process of discrete point cloud downsampling and regularized grid reconstruction, conventional grid algorithms usually use the arithmetic mean of the elevations of each point in the grid cell as the representative elevation. When there are many deep water points and few shallow points that obstruct navigation in the grid, the elevations of shallow points are easily masked by the average result, resulting in the omission of critical shallow point extreme values ​​in the regularized model.

[0003] In the processing chain from discrete point clouds to 3D reconstruction of complex underwater terrain, in addition to conventional spatial filtering and raster smoothing methods, the industry has attempted to introduce algorithms such as voxel isosurface extraction and local mesh refinement to construct 3D terrain. For example, Chinese invention patent CN120876766B discloses a method for digitally displaying marine terrain based on 3D modeling. This method constructs a 3D spatial density field by interpolating multi-source data, extracts isosurfaces using an improved MarchingCubes algorithm to generate an initial terrain model, and refines high-curvature areas using local 3D Delaunay triangulation. However, this type of solution has limitations when applied to high-noise underwater acoustic sounding points. In cloud conditions, spatial interpolation of the voxel density field can easily mix and homogenize suspended clutter in the water body with the real terrain boundary, making it difficult to decouple water noise and micro-topographic slopes during the isosurface extraction stage. In survey areas with non-convex polygonal boundaries, such as channel turning sections and irregularly shaped harbor basins, general spatial triangulation algorithms usually establish topological connections based on the convex hull of all discrete points. Since the convex hull cannot accurately fit the concave corner of the boundary, false triangles that cross the outer boundary of the actual water area are easily generated in the concave corner region. If these regions are then geometrically split and spliced ​​manually, it is easy to cause discontinuities in the topological connections of the triangulation network or geometric gaps in the boundary, thereby affecting the continuity of the three-dimensional surface reconstruction and the accuracy of spatial volume integration.

[0004] Therefore, the technical problems to be solved by this invention are: how to distinguish between floating noise and real terrain slopes and retain terrain boundaries in the feature space of a two-dimensional image, how to retain local control shallow point extrema during rasterization downsampling, and how to automatically construct a continuous three-dimensional triangular mesh surface that fits the actual boundary under the strong constraint of non-convex topological boundaries. Summary of the Invention

[0005] To address the problems in the background technology, the specific technical solution adopted is as follows: an underwater terrain point cloud filtering and 3D modeling method, comprising the following steps: Step S101: Obtain the set of three-dimensional discrete sounding points underwater in the test area. Extract the sounding points within the same transmitted pulse cross section based on the acoustic pulse timing and instantaneous heading azimuth. Calculate the projection distance of the sounding point perpendicular to the heading as the abscissa and the original elevation as the ordinate to generate a two-dimensional sounding profile spatial point array. Step S102: Convert the two-dimensional sounding profile spatial point array into a dual-channel geometric feature tensor containing local point density channels and height gradient channels. Input the convolutional target detection network to predict the abnormal noise bounding box, remove the sounding points that fall into the abnormal noise bounding box, and output a pure underwater terrain point agglomeration composed of pure underwater terrain points. Step S103: Divide the space into regular grid cells. Using the preset design reference elevation surface as a reference, extract shallow points with positive relative elevation differences to form a shallow point subset. Calculate the representative water depth value of the grid cell based on the elevation mean and elevation standard deviation of the shallow point subset, and generate a regularized water depth matrix containing regularized water depth feature points. Step S104: Extract the outer contour of the survey area to construct a non-convex closed topological boundary ring. Initialize and construct a planar triangular network using pure underwater topographic points and regularized water depth feature points. Insert the non-convex closed topological boundary ring as a strong constraint line segment and flip the intersecting diagonals to make them common edges. Eliminate the outer triangles based on the parity of the number of intersections of the centroid rays. Map the internal triangle topology to three-dimensional space to construct a three-dimensional underwater digital elevation model.

[0006] Preferably, generating the two-dimensional depth sounding profile spatial lattice in step S101 includes: step S1011, determining the transducer in the... The horizontal center coordinates and instantaneous heading azimuth angle at the trigger of the first transmission pulse; Step S1012, for the first... For any discrete sounding point of each transmitted pulse cross section, the lateral projection distance of the sounding point in the direction perpendicular to the flight direction is calculated based on the horizontal center coordinate and the instantaneous azimuth angle and used as the abscissa of the profile, and the original elevation is used as the ordinate of the profile; in step S1013, a two-dimensional sounding profile spatial point array is generated on the basis of a single transmitted pulse cross section, and the two-dimensional sounding profile spatial point array corresponding to each transmitted pulse cross section is generated pulse by pulse along the flight path according to the transmission pulse time sequence for feeding into the convolutional target detection network.

[0007] Preferably, step S102, converting to a dual-channel geometric feature tensor containing a local point density channel and a height gradient channel, includes: step S1021, dividing the local two-dimensional profile coordinate system into local pixel units with a spatial resolution of 0.05m; step S1022, counting the number of discrete sounding points falling into each pixel unit, and constructing a first feature channel representing the local point density; step S1023, calculating the first-order difference between the representative elevation of the point set inside each column of pixel units and the elevation of the adjacent column, and constructing a second feature channel representing the normalized height gradient; step S1024, combining the first feature channel and the second feature channel to generate a dual-channel two-dimensional geometric feature tensor, and inputting the dual-channel two-dimensional geometric feature tensor into a single-stage convolutional detection network containing a feature pyramid, outputting anomaly noise bounding boxes and their classification confidence scores; comparing the classification confidence scores with a preset judgment threshold, when the classification confidence score is greater than or equal to the preset judgment threshold, sounding points falling within the geometric range of the corresponding bounding box are removed, and the preset judgment threshold ranges from 0.80 to 0.90.

[0008] Preferably, step S104, inserting a non-convex closed topological boundary ring as a strongly constrained line segment and flipping the intersecting diagonals to make them common edges, includes: step S1041, extracting the actual outer boundary line segments of the survey area to construct an ordered closed point sequence; step S1042, using pure underwater topographic points and regularized water depth feature points as input discrete node sets, constructing an initial standard Delaunay triangulation on a two-dimensional plane; step S1043, traversing each strongly constrained line segment in the non-convex closed topological boundary ring one by one, searching for the triangle edges in the initial standard Delaunay triangulation that geometrically intersect with the current strongly constrained line segment, and performing local diagonal flipping on a pair of common triangles containing the intersecting edges, until each strongly constrained line segment is completely embedded in the triangulation to become a common edge of adjacent triangles.

[0009] Preferably, step S104, which involves eliminating external triangles based on the parity of the number of intersections of the centroid rays, includes: step S1044, calculating the geometric centroid coordinates of each triangle in the planar triangulation on the two-dimensional projection plane; step S1045, emitting a horizontal ray from each geometric centroid along the positive direction of the horizontal coordinate axis, and traversing and counting the total number of geometric intersections between the horizontal ray and all line segments in the non-convex closed topological boundary ring; step S1046, when the total number of geometric intersections is even, determining that the corresponding triangle is located outside the geometric boundary of the non-convex survey area and eliminating it; when the total number of geometric intersections is odd, determining that the corresponding triangle is located inside the geometric boundary of the non-convex survey area and retaining it; mapping the two-dimensional topological structure of all retained internal triangles to three-dimensional space, and assigning geodetic water depth and elevation values ​​to each grid vertex to generate a three-dimensional irregular triangulation surface model as a three-dimensional underwater digital elevation model.

[0010] Preferably, after constructing the three-dimensional underwater digital elevation model, the method further includes a dynamic extraction step of cross-sectional morphology: Step S105, in the constructed three-dimensional underwater digital elevation model, specify the starting point and ending point of the profile to generate the profile baseline segment, solve the spatial intersection sequence of the profile baseline segment and the boundary edge of each spatial triangle in the three-dimensional triangular network, arrange the spatial intersection sequence in linear ascending order according to the horizontal projection distance from the starting point of the profile, and interpolate between the spatial intersection points to generate a continuous curve of the transverse elevation change of the channel, and compare the continuous curve of the transverse elevation change of the channel with the design section slope line to output the verification data of the navigable bottom width and slope ratio.

[0011] Preferably, after constructing the three-dimensional underwater digital elevation model, the method also includes a multi-stage topographic evolution differential analysis step: Step S106, obtain the baseline three-dimensional topographic surface model of the same survey area, and perform collinear vertical projection difference between the three-dimensional underwater digital elevation model constructed in the current stage and the baseline three-dimensional topographic surface model to construct a spatial continuous differential field reflecting the changes in scour and deposition. When the differential field value is greater than zero, it is marked as a sediment deposition and backfilling area. When the differential field value is less than zero, it is marked as a water flow scour or excavation removal area. Based on the differential value gradient, a spatial scour and deposition color patch distribution map is generated.

[0012] Preferably, after constructing the three-dimensional underwater digital elevation model, the method further includes an engineering earthwork volume integration step: Step S107, extract all three-dimensional micro-triangles contained within the specified engineering calculation area, for any m-th three-dimensional micro-triangle, calculate the projected area of ​​the micro-triangle on the horizontal plane based on the planar coordinates of its three vertices, obtain the excavation depth of the three vertices relative to the design excavation reference plane, multiply the projected area by one-third of the sum of the excavation depths of the three vertices as the earthwork volume of the micro-triangular prism element, perform algebraic summation of the earthwork volumes of all micro-triangular prism elements within the specified engineering calculation area, and solve for the total engineering earthwork excavation volume of the area.

[0013] Preferably, after generating the regularized water depth matrix, a navigation water depth safety early warning step is also included: Step S108, setting a waterway maintenance water depth reference threshold, traversing and searching abnormal grid nodes in the regularized water depth matrix whose water depth values ​​are less than the waterway maintenance water depth reference threshold, extracting the plane coordinates and water depth elevation of the abnormal grid nodes, comparing the water depth elevation of the abnormal grid nodes with the preset design reference elevation surface to calculate the superelevation difference and plane influence range of local navigation obstruction shallow points, and outputting a spatial early warning mark indicating the location of the waterway under-dredging obstacle.

[0014] The beneficial effects of this invention are: 1. In underwater topographic point cloud filtering and 3D modeling, by converting the 3D discrete point cloud into a 2D profile orthogonal to the heading, and combining local point density and height gradient for abnormal noise identification, it is possible to more accurately distinguish between suspended clutter, bubble interference and elevation changes of the real terrain. This reduces the accidental deletion of channel slope breaks and deep channel boundaries during the filtering process. While filtering out noise, it can better preserve the original geometric continuity and orientation of the real underwater terrain, and improve the reliability of the point cloud data used in subsequent 3D terrain reconstruction.

[0015] 2. In regular spatial grid processing, the representative water depth of the grid cell is determined by the elevation characteristics of the shallow point subset. This allows the representative water depth to better reflect the actual situation of local shallow water points, avoiding the situation where a large number of deep water points are averaged to cover up a small number of navigation-obstructing shallow points. As a result, the extreme values ​​of local navigation-controlling shallow points can be preserved during point cloud downsampling and regularization reconstruction, reducing the possibility of key shallow points being smoothed out and providing more reliable topographic data for subsequent water depth safety verification.

[0016] 3. By constraining the boundary of the triangular mesh through a non-convex closed topological boundary ring and eliminating triangles outside the actual survey area boundary, the problem of cross-boundary triangles at the concave corners of the non-convex boundary can be reduced. Compared with relying on manual splitting and splicing of the survey area, this method can reduce the impact of topological tearing, geometric gaps and boundary distortion, making the generated three-dimensional irregular triangular mesh surface more continuous and better fit the actual boundary of the complex non-convex survey area.

[0017] 4. In engineering spatial analysis and comprehensive applications, the constructed three-dimensional underwater digital elevation model and regularized water depth matrix can continuously extract cross-sectional elevations along spatial intersections and accurately verify bottom width and slope ratio. Through multi-stage projection difference and color gradient rendering, the evolution and distribution of scour and sedimentation can be intuitively presented. Combined with the volume integral of the micro-element triangular prism unit, the calculation accuracy of earthwork excavation in non-convex irregular areas can be improved. Furthermore, through threshold comparison, proactive safety warnings can be provided for local navigation-obstructing shallow points. This comprehensively improves the operational efficiency and solution reliability of engineering measurement, morphological verification, and safety maintenance in complex water areas. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a flowchart of the underwater point cloud feature filtering and digital elevation modeling process of the present invention; Figure 2 This is a diagram illustrating the application architecture of underwater topographic data transfer and spatial analysis in this invention. Figure 3 This is a schematic diagram of the spatial distribution of noise bounding boxes in the two-dimensional cross-sectional coordinate system of this invention; Figure 4 This is an exploded diagram of the core logic control elements for constructing the digital elevation surface of this invention. Detailed Implementation

[0020] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0021] A method for filtering and 3D modeling underwater terrain point clouds includes the following steps: Step S101: Obtain the set of three-dimensional discrete sounding points underwater in the test area. Extract the sounding points within the same transmitted pulse cross section based on the acoustic pulse timing and instantaneous heading azimuth. Calculate the projection distance of the sounding point perpendicular to the heading as the abscissa and the original elevation as the ordinate to generate a two-dimensional sounding profile spatial point array. Step S102: Convert the two-dimensional sounding profile spatial point array into a dual-channel geometric feature tensor containing local point density channels and height gradient channels. Input the convolutional target detection network to predict the abnormal noise bounding box, remove the sounding points that fall into the abnormal noise bounding box, and output a pure underwater terrain point agglomeration composed of pure underwater terrain points. Step S103: Divide the space into regular grid cells. Using the preset design reference elevation surface as a reference, extract shallow points with positive relative elevation differences to form a shallow point subset. Calculate the representative water depth value of the grid cell based on the elevation mean and elevation standard deviation of the shallow point subset, and generate a regularized water depth matrix containing regularized water depth feature points. Step S104: Extract the outer contour of the survey area to construct a non-convex closed topological boundary ring. Initialize and construct a planar triangular network using pure underwater topographic points and regularized water depth feature points. Insert the non-convex closed topological boundary ring as a strong constraint line segment and flip the intersecting diagonals to make them common edges. Eliminate the outer triangles based on the parity of the number of intersections of the centroid rays. Map the internal triangle topology to three-dimensional space to construct a three-dimensional underwater digital elevation model.

[0022] Preferably, generating the two-dimensional depth sounding profile spatial lattice in step S101 includes: step S1011, determining the transducer in the... The horizontal center coordinates and instantaneous heading azimuth angle at the trigger of the first transmission pulse; Step S1012, for the first... For any discrete sounding point of each transmitted pulse cross section, the lateral projection distance of the sounding point in the direction perpendicular to the flight direction is calculated based on the horizontal center coordinate and the instantaneous azimuth angle and used as the abscissa of the profile, and the original elevation is used as the ordinate of the profile; in step S1013, a two-dimensional sounding profile spatial point array is generated on the basis of a single transmitted pulse cross section, and the two-dimensional sounding profile spatial point array corresponding to each transmitted pulse cross section is generated pulse by pulse along the flight path according to the transmission pulse time sequence for feeding into the convolutional target detection network.

[0023] Preferably, step S102, converting to a dual-channel geometric feature tensor containing a local point density channel and a height gradient channel, includes: step S1021, dividing the local two-dimensional profile coordinate system into local pixel units with a spatial resolution of 0.05m; step S1022, counting the number of discrete sounding points falling into each pixel unit, and constructing a first feature channel representing the local point density; step S1023, calculating the first-order difference between the representative elevation of the point set inside each column of pixel units and the elevation of the adjacent column, and constructing a second feature channel representing the normalized height gradient; step S1024, combining the first feature channel and the second feature channel to generate a dual-channel two-dimensional geometric feature tensor, and inputting the dual-channel two-dimensional geometric feature tensor into a single-stage convolutional detection network containing a feature pyramid, outputting anomaly noise bounding boxes and their classification confidence scores; comparing the classification confidence scores with a preset judgment threshold, when the classification confidence score is greater than or equal to the preset judgment threshold, sounding points falling within the geometric range of the corresponding bounding box are removed, and the preset judgment threshold ranges from 0.80 to 0.90.

[0024] Preferably, step S104, inserting a non-convex closed topological boundary ring as a strongly constrained line segment and flipping the intersecting diagonals to make them common edges, includes: step S1041, extracting the actual outer boundary line segments of the survey area to construct an ordered closed point sequence; step S1042, using pure underwater topographic points and regularized water depth feature points as input discrete node sets, constructing an initial standard Delaunay triangulation on a two-dimensional plane; step S1043, traversing each strongly constrained line segment in the non-convex closed topological boundary ring one by one, searching for the triangle edges in the initial standard Delaunay triangulation that geometrically intersect with the current strongly constrained line segment, and performing local diagonal flipping on a pair of common triangles containing the intersecting edges, until each strongly constrained line segment is completely embedded in the triangulation to become a common edge of adjacent triangles.

[0025] Preferably, step S104, which involves eliminating external triangles based on the parity of the number of intersections of the centroid rays, includes: step S1044, calculating the geometric centroid coordinates of each triangle in the planar triangulation on the two-dimensional projection plane; step S1045, emitting a horizontal ray from each geometric centroid along the positive direction of the horizontal coordinate axis, and traversing and counting the total number of geometric intersections between the horizontal ray and all line segments in the non-convex closed topological boundary ring; step S1046, when the total number of geometric intersections is even, determining that the corresponding triangle is located outside the geometric boundary of the non-convex survey area and eliminating it; when the total number of geometric intersections is odd, determining that the corresponding triangle is located inside the geometric boundary of the non-convex survey area and retaining it; mapping the two-dimensional topological structure of all retained internal triangles to three-dimensional space, and assigning geodetic water depth and elevation values ​​to each grid vertex to generate a three-dimensional irregular triangulation surface model as a three-dimensional underwater digital elevation model.

[0026] Preferably, after constructing the three-dimensional underwater digital elevation model, the method further includes a dynamic extraction step of cross-sectional morphology: Step S105, in the constructed three-dimensional underwater digital elevation model, specify the starting point and ending point of the profile to generate the profile baseline segment, solve the spatial intersection sequence of the profile baseline segment and the boundary edge of each spatial triangle in the three-dimensional triangular network, arrange the spatial intersection sequence in linear ascending order according to the horizontal projection distance from the starting point of the profile, and interpolate between the spatial intersection points to generate a continuous curve of the transverse elevation change of the channel, and compare the continuous curve of the transverse elevation change of the channel with the design section slope line to output the verification data of the navigable bottom width and slope ratio.

[0027] Preferably, after constructing the three-dimensional underwater digital elevation model, the method also includes a multi-stage topographic evolution differential analysis step: Step S106, obtain the baseline three-dimensional topographic surface model of the same survey area, and perform collinear vertical projection difference between the three-dimensional underwater digital elevation model constructed in the current stage and the baseline three-dimensional topographic surface model to construct a spatial continuous differential field reflecting the changes in scour and deposition. When the differential field value is greater than zero, it is marked as a sediment deposition and backfilling area. When the differential field value is less than zero, it is marked as a water flow scour or excavation removal area. Based on the differential value gradient, a spatial scour and deposition color patch distribution map is generated.

[0028] Preferably, after constructing the three-dimensional underwater digital elevation model, the method further includes an engineering earthwork volume integration step: Step S107, extract all three-dimensional micro-triangles contained within the specified engineering calculation area, for any m-th three-dimensional micro-triangle, calculate the projected area of ​​the micro-triangle on the horizontal plane based on the planar coordinates of its three vertices, obtain the excavation depth of the three vertices relative to the design excavation reference plane, multiply the projected area by one-third of the sum of the excavation depths of the three vertices as the earthwork volume of the micro-triangular prism element, perform algebraic summation of the earthwork volumes of all micro-triangular prism elements within the specified engineering calculation area, and solve for the total engineering earthwork excavation volume of the area.

[0029] Preferably, after generating the regularized water depth matrix, a navigation water depth safety early warning step is also included: Step S108, setting a waterway maintenance water depth reference threshold, traversing and searching abnormal grid nodes in the regularized water depth matrix whose water depth values ​​are less than the waterway maintenance water depth reference threshold, extracting the plane coordinates and water depth elevation of the abnormal grid nodes, comparing the water depth elevation of the abnormal grid nodes with the preset design reference elevation surface to calculate the superelevation difference and plane influence range of local navigation obstruction shallow points, and outputting a spatial early warning mark indicating the location of the waterway under-dredging obstacle.

[0030] Example 1: In the turning section of the main channel and the concave corner waters of irregularly shaped berths in the port area, the water contains suspended sediment particles and cavitation microbubbles from ships. The channel slope has a sharply undulating slope line, and there are local isolated shallow obstacles at the bottom. The outer contour of the survey area presents a non-convex boundary shape containing a concave polygonal outline. The measurement platform, equipped with a multibeam echo sounder, navigates along the preset survey line to collect underwater three-dimensional discrete echo sounder points after corrections for tide level, sound velocity, roll, pitch, and heave attitudes. ,in, For Gaussian projection plane coordinates, To represent water depth and geodetic elevation as negative values, the acoustic emission pulse trigger time and inertial navigation sampling time sequence are aligned based on a unified timestamp, and linear interpolation is performed on the heading data of adjacent sampling points to determine the measurement platform at the [missing information]. The horizontal geometric center coordinates when triggered by an acoustic emission pulse and instantaneous heading azimuth ,in Represent the acoustic emission pulse number; extract the set of sounding points within the same emission pulse cross section, and establish a local two-dimensional sounding profile coordinate system u,v with the survey line heading as the normal and the direction perpendicular to the heading as the transverse; for the part belonging to the first... Arbitrary discrete sounding points of each transmitted pulse cross section Calculate its lateral projection coordinates in the local two-dimensional sounding profile coordinate system. The original elevation As the longitudinal coordinate of the cross section The corresponding two-dimensional sounding profile spatial point array is generated on the basis of a single transmitted pulse cross section. Then, it is converted into a dual-channel geometric feature tensor pulse by pulse along the flight path according to the transmission pulse time sequence and then fed into the convolutional target detection network to predict the abnormal noise bounding box.

[0031] In the local two-dimensional sounding profile coordinate system u,v, local pixel units are divided with a spatial resolution of 0.05m. This 0.05m is the preferred spatial resolution for shallow water high-density sounding conditions. In actual operation, it can also be adaptively matched and adjusted according to the water depth of the test area and the beam footprint spacing. The number of discrete sounding points falling into each pixel unit is counted to construct the local point density channel as the first feature channel. When no depth point falls within a pixel unit, the local point density channel is set to zero. The first-order difference between the representative elevation of the point set within each column of pixel units and the elevation of the adjacent columns is calculated and normalized to construct the normalized height gradient channel as the second feature channel. A dual-channel two-dimensional geometric feature tensor is generated by combining the local point density channel and the normalized height gradient channel. Specifically, when constructing the normalized height gradient channel, for each column of pixel grid cells arranged horizontally in the local coordinate system, the mean ordinate of all sounding points in that column is extracted vertically as the representative elevation of the cross section at that horizontal position. When no sounding point falls into a column, the linear interpolation value of the representative elevation of the adjacent column is taken as the elevation of that column. Then, the absolute value of the central difference between the representative elevation of the current column and the representative elevation of the horizontally adjacent column is calculated, and divided by the spatial horizontal step size to obtain the local point density channel. The longitudinal slope of the terrain is calculated and normalized by division using a preset local maximum slope threshold. A dimensionless gradient scalar is then assigned to fill all pixel grid cells containing sounding points in that column, thus establishing a horizontally continuously varying terrain slope gradient feature in a two-dimensional raster space, decoupled from the vertical position of discrete points. A dual-channel two-dimensional geometric feature tensor is input into a single-stage convolutional object detection network containing a feature pyramid. Multi-scale receptive fields decouple shallow discrete patch features from deep global terrain trend features. The network detection head then predicts and outputs anomaly noise bounding boxes. The classification confidence score of the water noise category to which it belongs. ,in Indicates the bounding box number of the abnormal noise; compares the classification confidence score. With preset judgment threshold Among them, the preset judgment threshold Set to 0.85, its value ranges from 0.80 to 0.90; when When the noise falls within the corresponding abnormal noise bounding box, it is determined that the noise falls within the bounding box. Sounding points within the coverage area are marked and removed as noise points caused by suspended matter and air bubbles. Sounding points not falling within the bounding box are retained as topographic points. The output is a pure underwater topographic point cluster that maintains the geometric continuity of the slope polygon and the deep trench boundary. .

[0032] Divide the horizontal two-dimensional projection plane into steps with size... A regular spatial grid system, defining grid cells. The coordinates of the center of the plane are The spatial projection of the search space falls into the grid cell. A subset of clean sounding points within the range Based on the preset design reference elevation surface For reference, a pre-set design datum elevation surface is established. Take the negative geodetic elevation corresponding to the design navigable bottom elevation of the channel or the engineering excavation reference surface, for example, -15.00m, and calculate the elevation difference of each sounding point within the grid cell relative to the design reference surface. Extract to satisfy Shallow points form a subset of shallow points Determine the number of elements contained in the shallow subset. When there are shallow points within a grid cell and the number of shallow points is... When calculating the arithmetic mean of the elevations of the shallow point subset, ,when When calculating the standard deviation of shallow point samples... ,when Take the standard deviation of shallow point samples According to the formula Calculate the representative water depth values ​​for each grid cell, where extreme values ​​retain the sensitivity coefficient. The value is set to 0.25, ranging from 0.10 to 0.50, to statistically shift the elevation of the grid nodes towards shallow water extremes. This extreme value retention calculation is based on the statistical derivation principle of the one-sided confidence upper bound. Under the premise that the shallow water sounding point samples follow an approximately normal distribution, the grid represents the water depth towards the shallow water extremes by introducing linear weighting of the sample standard deviation. Through engineering measurements and calibration, when the extreme value retention sensitivity coefficient is below 0.10, the standard deviation adjustment weight is too small, and the prominent extreme characteristics of isolated shallow points will still be diluted by the arithmetic mean of deep water points, leading to missed measurements of shallow points in the navigation channel. When the extreme value retention sensitivity coefficient is above 0.50, it excessively amplifies the contribution of random measurement disturbances and gross errors to the grid elevation, resulting in false topographical features in flat seabed areas. Limiting it to between 0.10 and 0.50 and calibrating it to 0.25 can both fully capture the navigation-obstructing shallow point extremes and effectively suppress local false elevation jumps caused by random noise. When there are no shallow points within the grid cell and... At that time, the arithmetic mean of the elevations of all clean sounding points within the grid cell is taken. As the representative water depth value of the grid cell; when the number of depth measurement points M=0 in the grid cell, the initial value is calculated by neighborhood inverse distance weighted interpolation, and the regularized water depth matrix that preserves the spatial extremum of shallow points of local obstacles in the channel is output.

[0033] In the navigation depth safety early warning step, after generating a regularized depth matrix, a navigation maintenance depth benchmark threshold is set according to the channel maintenance requirements. Abnormal grid nodes in the regularized depth matrix with depth values ​​less than the navigation maintenance depth benchmark threshold are traversed and retrieved. The planar coordinates and depth elevations of the abnormal grid nodes are extracted. The depth elevations of the abnormal grid nodes are compared with the preset design benchmark elevation surface. The superelevation difference of local navigation-obstructing shallow points is calculated, and the planar influence range is determined by the range of the corresponding regular spatial grid unit on the horizontal projection surface. Spatial early warning markers indicating the location of under-dredging obstacles in the channel are output. The preset design benchmark elevation surface is determined based on the channel design navigation bottom elevation. The design excavation benchmark surface in step S107 is determined based on the engineering dredging design excavation elevation. When there is no superelevation requirement, it coincides with the preset design benchmark elevation surface; when there is a superelevation requirement, it is shifted downwards according to the superelevation thickness. The navigation maintenance depth benchmark threshold in step S108 is determined based on the minimum allowable depth required for navigation safety. All three benchmarks are independent and have clear sources.

[0034] Based on the sequential connection of the boundary coordinate files of the engineering survey area, extract the polyline nodes of the actual outer contour of the engineering survey area to construct a non-convex closed topological boundary loop containing the concave polygon contour. Each boundary segment All are strongly constrained line segments that cannot be crossed. The boundary segment number is used as the input discrete node set, with filtered pure underwater topographic points and regularized water depth feature points as input. An initial standard Delaunay triangulation is constructed on a two-dimensional plane. In the actual construction of the input discrete node set, a spatial partitioning decoupling and fusion strategy is adopted: within the regular spatial grid cell, all original pure sounding points in the grid cell that has been selected for the shallow point subset are removed, and the calculated regularized water depth feature points are used as the unique control node input within that grid range; while in flat grids where shallow points have not been extracted and in areas with severe undulations such as channel slope breaks, the original pure underwater topographic points are completely retained as discrete input nodes. This eliminates the elevation numerical conflicts and dense redundancy of multiple nodes under the same horizontal spatial coordinates, ensuring the topological uniformity and three-dimensional surface smoothness of the initial triangulation in the planar projection domain. The non-convex closed topological boundary loops are traversed one by one. Each strongly constrained line segment in Search the planar triangulation network for line segments that are currently under strong constraints. For triangles whose sides intersect in spatial geometry, perform a local diagonal flip on the pair of triangles containing the intersecting sides, until all strongly constrained line segments are reached. The triangulation is fully embedded and becomes a shared edge of adjacent triangles. If, during the diagonal flipping process, the combined quadrilateral formed by a pair of shared-edge triangles containing intersecting edges is a non-convex quadrilateral, or if the new diagonal after flipping still overlaps with the strongly constrained line segment, the global flipping is paused. Instead, the geometric intersection of the strongly constrained line segment and the edge of the intersecting triangle is added to the triangulation as a temporary insertion node. The intersecting triangle is locally subdivided into sub-triangles with the intersection point as the vertex. Local convexity detection and diagonal flipping are then performed on the subdivided shared-edge diagonals. After the strongly constrained line segment becomes a shared edge of adjacent triangles in the triangulation, the temporary insertion node is removed through the collinear redundant edge merging algorithm to restore the original geometric integrity of the strongly constrained line segment, avoiding flipping deadlock and topological tearing. The geometric centroid coordinates of each triangle in the planar triangulation on the two-dimensional projection plane are calculated. A horizontal ray is emitted from each geometric centroid along the positive direction of the horizontal coordinate axis. The horizontal ray and the non-convex closed topological boundary loop are traversed and statistically analyzed. The total number of geometric intersections of all line segments is calculated. When the total number of geometric intersections is even, the corresponding triangle is determined to be outside the geometric boundary of the non-convex survey area and is topologically removed. When the total number of geometric intersections is odd, the corresponding triangle is determined to be inside the geometric boundary of the non-convex survey area and is retained. The two-dimensional topological structure of all retained internal triangles is mapped to the three-dimensional elevation space, and each grid vertex is assigned a geodetic water depth elevation value to generate a three-dimensional irregular triangular mesh surface model as a three-dimensional underwater digital elevation model.

[0035] In the multi-period topographic evolution differential analysis step, a baseline three-dimensional topographic surface model of the same survey area is obtained. The three-dimensional underwater digital elevation model constructed in the early stage, namely the three-dimensional irregular triangular mesh surface model, will be used. Compared with the baseline three-dimensional terrain surface model Collinear vertical projection differences are used to construct a spatially continuous difference field reflecting changes in scouring and silting. When the difference field value At that time, it was determined to be a siltation and backfilling area, when the difference field value When an area is identified as being eroded by water flow or excavated, a spatial erosion and sedimentation color distribution map is generated based on differential numerical gradient rendering. Specifically, a color mapping table is set with equal intervals of 0.1m erosion and sedimentation thickness. For sedimentation areas with differential values ​​greater than zero, a warm hue gradient is applied increasing with sedimentation thickness; for erosion or excavation areas with differential values ​​less than zero, a cool hue gradient is applied increasing with erosion depth. During this process, since water depth and geodetic elevation are uniformly represented using a negative value system, the current elevation is greater than the baseline elevation. Corresponding to the physical state of increased absolute elevation and shallower water depth, an algebraic difference of 0 meters between two periods of elevation directly indicates the thickness of sediment deposition. Conversely, when the current period elevation is less than the baseline period elevation, corresponding to the physical state of decreased absolute elevation and deeper water depth, an algebraic difference of 0 meters between the two periods of elevation directly indicates the depth of water scouring or excavation. The absolute value of the algebraic difference is the actual scouring and deposition thickness per unit projected area, thus ensuring the logical consistency between the topographic evolution judgment criteria and the physical scouring and deposition mechanism under the negative elevation system.

[0036] In the earthwork volume integration step, extract all three-dimensional infinitesimal triangles contained within the specified non-convex engineering calculation region. ,in Represents the index of a three-dimensional infinitesimal triangle; for any nth triangle... A three-dimensional infinitesimal triangle, based on the planar coordinates of its three vertices. Calculate the horizontal projected area Obtain the excavation depth of each of the three vertices relative to the designed excavation reference plane. Calculate the earthwork volume of the infinitesimal triangular prism element. The earthwork volume of all U infinitesimal triangular prism elements within the computational region is algebraically summed to output the total earthwork excavation volume for the region. .

[0037] In the dynamic extraction step of cross-sectional morphology, the starting point of the section is specified in the three-dimensional irregular triangular mesh surface model. and the finish line Generate a baseline segment for the profile, calculate the sequence of spatial intersection points between the baseline segment and the boundary edges of each spatial triangle in the 3D triangular network, arrange the sequence of spatial intersection points in linear ascending order according to the horizontal projection distance from the profile starting point A, interpolate between the spatial intersection points to generate a continuous curve of lateral elevation change of the waterway, and compare the continuous curve of lateral elevation change of the waterway with the design section slope line to output the verification data of navigable bottom width and slope ratio.

[0038] Example 2: This example selects a surface measurement platform equipped with a multibeam echo sounder. Based on measured underwater acoustic point cloud data, a semi-physical simulation test case is constructed to perform three-dimensional spatial modeling verification on the turning section of the main channel of the port area and the concave corner waters of irregular berths. The water in the measurement area is mixed with suspended sediment particles and cavitation microbubbles from ships. The channel slope has a sharply undulating slope line, and there are local isolated shallow point obstacles at the bottom. The outer contour of the measurement area presents a non-convex boundary shape containing a concave polygonal outline. The measurement platform collects the original underwater three-dimensional discrete echo sounder point cloud. ,in, For Gaussian projection plane coordinates, To represent water depth and geodetic elevation as negative values, the operating frequency of the multibeam echo sounder system was set to 400 kHz, and the transmission opening angle was limited to 1 degree x 2 degrees. During the initial signal acquisition process, Gaussian white noise with a signal-to-noise ratio of 15 dB and power frequency interference harmonics at a frequency of 50 Hz were superimposed to simulate electromagnetic and acoustic disturbances in a real underwater acoustic detection environment, determining the measurement platform's position on the [missing information - likely a specific location or stage]. The instantaneous horizontal center coordinates at the moment of triggering an acoustic emission pulse and instantaneous heading azimuth Extract the set of sounding points within the same transmitted pulse cross section, and establish a local two-dimensional sounding profile coordinate system u,v with the survey line heading as the normal and the direction perpendicular to the heading as the lateral direction; calculate the lateral projection coordinates of the discrete sounding points in this coordinate system. The original elevation As the longitudinal coordinate of the cross section A two-dimensional sounding profile spatial lattice is generated using a single transmitted pulse cross section as a unit. Then, it is converted into a two-channel two-dimensional geometric feature tensor pulse by pulse along the flight path according to the transmission pulse time sequence and input into a single-stage convolutional target detection network.

[0039] In the local two-dimensional sounding profile coordinate system u,v, local pixel units are divided with a spatial resolution of 0.05m. The number of discrete sounding points falling into each pixel unit is counted, and a local point density channel is constructed as the first feature channel. Calculate the first-order difference between the representative elevation of the point set within each column of pixels and the elevation of the adjacent columns, and normalize it to construct a normalized height gradient channel as the second feature channel. The two feature channels are combined to generate a dual-channel two-dimensional geometric feature tensor. This tensor is then input into a single-stage convolutional object detection network. Multi-scale receptive fields are used to extract shallow discrete patch features and deep global terrain trend features. The network's detection head then predicts bounding boxes for abnormal noise. and its classification confidence score Compare classification confidence scores With preset judgment threshold Remove clutter and output a pure underwater terrain point cluster. To verify the influence between the average of isolated shallow point extrema and the prominence of false shallow points caused by noise residue during the discrete point regularization downsampling process, four parallel experimental groups were set up: Comparison group one used a spatial grid division algorithm that combined global elevation threshold denoising with arithmetic mean; comparison group two used a single-stage convolutional target detection network for filtering, and retained the sensitivity coefficient of the extrema in the grid reconstruction stage. The value is set to 0.05, which is lower than the lower limit. In comparison, sample group three uses a single-stage convolutional target detection network for filtering, and retains the sensitivity coefficient for extreme values. The value was set to 0.65, which is higher than the upper limit. The experimental group adopted a complete grid reconstruction logic based on the preset design reference elevation surface. Extracting shallow point subsets Calculate the arithmetic mean of shallow point elevations. Standard deviation of shallow sample According to the formula The representative water depth values ​​of the computational grid cells are merged, and the extreme values ​​are preserved by the sensitivity coefficient. Set the calibration value to 0.25, which is within the range of 0.10 to 0.50.

[0040] Divide the horizontal two-dimensional projection plane into steps with size... The regular spatial grid system selects six continuous grid cells within the area of ​​frequent shallow obstacles characteristic of the waterway. The plane centers are sequentially labeled G1 to G6, and the pre-set design reference elevation surface is... -15.00m was used as the reference depth. The absolute elevations of the true shallow points of these six grid cells were verified by high-density sonar. The values ​​of G1 to G6 were -14.62m, -14.85m, -14.71m, -14.58m, -14.92m and -14.68m, respectively. The representative water depth values ​​for each sample group were calculated for the same set of grid cells, and the elevation jump deviation scalar was obtained according to the same test procedure. Within the same grid cell, three repeated independent static scans were performed using a multibeam system with a fixed heading. The absolute value of the difference between the maximum and minimum values ​​of the representative water depth values ​​of the grid cell in the three independent bathymetry results was taken as the elevation jump deviation scalar, which was used to characterize the stability of the reconstructed elevation under power frequency interference and Gaussian random fluctuation conditions.

[0041] After global elevation threshold denoising and arithmetic mean grid reconstruction, the water depth values ​​represented by grid cells G1 to G6 in the comparison sample group are -15.16m, -15.28m, -15.28m, -15.04m, -15.43m, and -15.12m, respectively, with corresponding elevation jump deviation scalar values ​​of 0.27m, 0.33m, 0.25m, 0.35m, 0.26m, and 0.30m, respectively. Compared with the aforementioned absolute elevation of the actual shallow points, arithmetic mean smoothing makes the water depth values ​​represented by the grid cells more than 0.40m lower than the actual shallow point elevation, and the deep water background point set averages the extreme values ​​of the shallow points. At the same time, this sample group did not use multi-scale receptive field feature extraction, and the average jump error corresponding to power frequency harmonics and Gaussian perturbations reached 0.29m.

[0042] Compared to sample group two, after filtering in a single-stage convolutional object detection network, the sensitivity coefficient for retaining extreme values ​​is improved. With a value set to 0.05, the representative water depth values ​​for the same six grid cells are -14.94m, -15.16m, -14.98m, -14.93m, -15.15m, and -15.03m, respectively. The elevation jump deviation scalars are 0.03m, 0.04m, 0.02m, 0.04m, 0.03m, and 0.03m, respectively. After front-end filtering, the jump error of each grid cell is reduced to within 0.04m; however, Values ​​below 0.10 are considered lower limits; shallow sample standard deviation. The weights in the fusion calculation are too low, and the water depth values ​​represented by the grid cells still produce a navigational obstruction deviation of more than 0.20m compared with the corresponding actual shallow point elevations.

[0043] The third comparison sample also uses a single-stage convolutional target detection network filter to retain the sensitivity coefficients of the extreme values. Set to 0.65, the grid cells G1 to G6 represent water depths of -14.49m, -14.63m, -14.56m, -14.33m, -14.80m, and -14.49m respectively, with elevation jump deviation scalars of 0.21m, 0.25m, 0.18m, 0.24m, 0.19m, and 0.22m respectively. The reason for the higher-than-0.50 value is that although the front-end single-stage convolutional object detection network has removed suspended clutter points in the water, The linear amplification represents the random measurement fluctuation component of the standard deviation of shallow samples in three repeated scans, when When the standard deviation reaches 0.65, the elevation disturbance introduced by the weighted standard deviation term has exceeded the stability gain brought by the convolutional filter, resulting in an increased impact of local high-frequency topographic undulations on the water depth calculation represented by the grid cells, and the residual elevation jump error is in the range of 0.18m to 0.25m.

[0044] The experimental group used a single-stage convolutional object detection network to remove clutter points within the abnormal noise bounding boxes and retained the sensitivity coefficients by retrieving extrema. Set it to 0.25 within the range of 0.10 to 0.50, and then use the arithmetic mean of shallow point elevations. Standard deviation of shallow sample The representative water depth values ​​of the fused computational grid cells, G1 to G6, are -14.63m, -14.89m, -14.73m, -14.61m, -14.93m, and -14.72m, respectively. The corresponding elevation jump deviation scalars are 0.03m, 0.04m, 0.02m, 0.03m, 0.04m, and 0.03m, respectively. Front-end clutter removal reduces the impact of false elevation jumps on shallow point statistics. Subsequent values ​​are calculated using the arithmetic mean of the shallow point elevations. Standard deviation of shallow sample The water depth was calculated by jointly participating in the calculation of the water depth represented by the grid cells; it was compared one by one with the absolute elevation of the actual shallow points in the six grid cells. The residual between the output water depth value represented by the grid cells and the actual shallow point elevation was controlled within 0.04m. The data of the four sample groups showed a corresponding convergence trend as the extreme value retention sensitivity coefficient changed: when When the value is below the lower limit of 0.05, there is still a significant under-detection bias in the shallow point extreme values; when When the elevation jump error is increased to 0.65, exceeding the upper limit, the error increases again; When the value is set to 0.25 within the range of 0.10 to 0.50, the convolutional target detection network first removes clutter points within the abnormal noise bounding box, and then calculates the water depth represented by the grid cells using the arithmetic mean of shallow point elevations and the standard deviation of shallow point samples, forming a continuous numerical processing process from front-end noise filtering to regular spatial grid reconstruction.

[0045] Example 3: The current surface measurement platform is equipped with a multibeam echo sounder and an inertial navigation positioning system. It advances along a preset route in the widened area of ​​the main channel bend and irregularly shaped waterways of the port area, collecting underwater three-dimensional discrete sounding point aggregates after real-time tide level correction and layering correction of the measured sound velocity profile in the measurement area. ,in For Gaussian projection plane coordinates, The water depth and geodetic elevation are expressed as negative values. The transducer is calculated in the [missing information - likely a specific step or condition]. Instantaneous horizontal geometric center coordinates at the moment of acoustic emission pulse triggering and instantaneous heading azimuth Extract the set of sounding points within the same transmitted pulse cross section, and establish a local two-dimensional sounding profile coordinate system u,v with the survey line heading as the normal and the direction perpendicular to the heading as the transverse direction; for the points belonging to the first... Arbitrary discrete sounding points of each transmitted pulse cross section Calculate the transverse projection coordinates based on the coordinate projection transformation relationship. The original elevation As the longitudinal coordinate of the cross section A two-dimensional sounding profile spatial lattice is generated using a single transmitted pulse cross-section as the unit. Each pulse along the flight path is then converted into a dual-channel two-dimensional geometric feature tensor according to the transmission pulse sequence and input into a single-stage convolutional target detection network. In a Gaussian projection plane rectangular coordinate system, the vertical axis points due north, and the horizontal axis points due east. The instantaneous heading azimuth angle is the angle between the north direction and the heading propulsion direction, rotated clockwise. The unit direction vector along the heading propulsion direction consists of cosine and sine components. The transverse orthogonal projection unit vector perpendicular to the right of the heading consists of sine and negative cosine components. The planar position vector of the sounding point relative to the transducer center is projected onto this transverse orthogonal axis to obtain the transverse projection distance perpendicular to the survey line heading, and the two-dimensional sounding profile coordinate transformation is completed accordingly.

[0046] During the transformation of the spatial point matrix of the two-dimensional sounding profile, the local two-dimensional sounding profile coordinate system u,v is divided into local pixel grid units with a spatial resolution of 0.05m. The total number of discrete sounding points falling into each pixel grid unit is counted, and a local point density channel is constructed as the first feature channel. The first-order central difference quotient representing the elevation between adjacent column pixel grid cells is calculated along the lateral projection coordinate u direction. After normalization by local maximum slope threshold division, a normalized height gradient channel is generated as the second feature channel. The two feature channels are concatenated along the depth dimension to form a dual-channel two-dimensional geometric feature tensor. This dual-channel two-dimensional geometric feature tensor is then input into a single-stage convolutional object detection network. The network backbone uses multiple residual layers with depth-separable convolutional kernels to extract shallow high-resolution features and deep large receptive field features. The feature pyramid aggregates deep continuous slope geometric features and shallow discrete patch location information from top to bottom. A multi-scale detection head then regresses and predicts the bounding boxes for abnormal noise on the feature layer. The classification confidence score of the water noise category to which it belongs. When training this single-stage convolutional object detection network, at least 3,000 historical two-dimensional sounding profile geometric feature tensors, containing typical water body suspended clutter, cavitation microbubbles, and underwater slope breakage features, covering at least three typical waterways and port basin survey areas, were pre-collected as the training dataset. At least 8,000 real two-dimensional bounding boxes and noise category labels were manually labeled for isolated suspended noise clusters and dense clutter regions (including at least 3,500 samples of water body suspended sediment clutter, at least 2,500 samples of cavitation microbubbles, and at least 2,000 samples of underwater slope breakage features). An intersection-over-union ratio (IoU) greater than 0.50 was used as the positive sample candidate matching criterion, and an IoU less than 0.30 was used as the negative sample candidate matching criterion. The negative sample criterion is used to balance the positive and negative sample candidate matching at a ratio of 1:3. Network training employs a multi-task loss function for joint optimization. This loss function consists of a smoothed first-order absolute error regression loss (based on the differences in horizontal coordinates, vertical elevation, width, and height between the predicted and ground truth boxes) and a cross-entropy classification loss (based on the predicted probability of noise categories), weighted 1:1. The input training set tensor is augmented through random horizontal flipping and minor coordinate translation. The initial learning rate is set to 0.001, and iterative training using mini-batch gradient descent is employed until the total loss value stabilizes below 0.01 and shows no decrease for 10 consecutive training epochs. Subsequently, the network weight parameters are determined, and online forward inference is performed, with a preset decision threshold. The threshold is set to 0.85, with an applicable range of 0.80 to 0.90. This threshold is calibrated based on experimental results regarding the filtration rate of suspended solids in water and the fidelity of terrain. When the preset threshold is below 0.80, the network classifies some continuous micro-topographical bounding boxes with large slopes as noise, and the false rejection rate of slope break feature points rises to over 3.5%. When the preset threshold is above 0.90, the recall ability for weakly scattered micro-cavitation bubbles and edge discrete sediment clutter decreases, and the noise residue rate exceeds 8.0%. Within the value range of 0.80 to 0.90, the filtration rate of suspended noise in water remains stable above 98.2%, and the false rejection rate of terrain slope break feature points is controlled within 0.3%. 0.85 is used as a compromise between noise rejection and terrain fidelity. When this happens, the object will fall into the corresponding abnormal noise bounding box. Sounding points within the coverage area are marked as suspended sediment or microbubble disturbance noise and removed. Sounding points not falling within the bounding box are retained as topographic points. The output is a pure underwater topographic point cluster that maintains the geometric continuity of the channel slope polygon and the deep channel boundary. .

[0047] In the process of reconstructing a regular spatial grid model from a cluster of underwater topographic points in pure water, the step size is divided on the horizontal two-dimensional projection plane. For 2m, For a regular spatial grid system of 2m, determine each grid cell. The coordinates of the center of the plane are The retrieved plane projection falls into the grid cell. A subset of clean sounding points within the range Based on the preset design reference elevation surface For reference, calculate the elevation difference between each sounding point within the grid cell and the design datum. Extract the depth measurement points with a difference greater than zero to form a shallow point subset. Count the number of shallow points When the number of shallow spots When calculating the arithmetic mean of the elevations of the shallow point subset, ,when Calculate the standard deviation of shallow point samples ,when Take the standard deviation of shallow point samples According to the formula Calculate the representative water depth value of the grid cell, where the extreme values ​​retain the sensitivity coefficient. Set to 0.25, its calibration value range is 0.10 to 0.50; when there are no shallow points in the grid cell and K=0, take the arithmetic mean of the elevations of all clean sounding points in the grid cell. As a representative water depth value; when there is no landing point within the grid cell and the total number of sounding points M=0, with the current grid center as the center, in Within the search radius, adjacent valid grid cells are retrieved. Inverse distance weighted interpolation is used to assign initial water depth values ​​and output a regularized water depth matrix. When performing inverse distance weighted interpolation, the inverse square root of the horizontal distance from the center of the current blank grid to the center of the adjacent valid grid is used as the interpolation weight, and each weight is normalized and then weighted averaged. When the number of adjacent valid grid cells retrieved within the 6m search radius is 0, the three valid grid cells with the closest Euclidean distance to the center of the current grid are searched throughout the entire survey area. The same inverse square distance weight is used to perform a weighted average, and the resulting value is assigned to the current grid as the initial water depth value, so that the missing nodes in the regularized water depth matrix obtain the corresponding water depth values.

[0048] During the continuous construction of the 3D mesh with non-convex boundary constraints, the polyline nodes of the harbor basin corners and the actual outer contour of the concave berthing area are extracted to construct a non-convex closed topological boundary loop containing the concave polygon contour. Each boundary segment All are treated as strongly constrained line segments; using pure underwater topographic points and regularized water depth feature points as input discrete node sets, an initial standard Delaunay triangulation is constructed on a two-dimensional plane, and non-convex closed topological boundary loops are traversed one by one. Strongly constrained line segments The process involves retrieving triangle edges from the initial standard Delaunay triangulation that intersect the current strongly constrained line segment spatially. Local diagonal flipping is applied to the pair of triangles sharing the intersecting edge until the strongly constrained line segment is fully embedded in the triangulation and becomes a shared edge of adjacent triangles. Subsequently, the geometric centroid coordinates of each triangle in the planar triangulation on the two-dimensional projection plane are calculated. A horizontal ray is emitted from each geometric centroid along the positive direction of the horizontal coordinate axis. When the ray coincides with a vertex of the boundary polygon, an application is applied to the ray direction vector. The minute angular perturbation of 0.000001 rad is used to traverse and statistically count the total number of geometric intersections between the horizontal ray and all line segments in the non-convex closed topological boundary ring. In a double-precision floating-point computing environment, the perturbation in the horizontal and vertical coordinate calculations caused by this minute angular deflection of 0.000001 rad is greater than the limit of floating-point machine computing precision. In an engineering survey area with a maximum span of 5000m, the geometric lateral displacement caused by this deflection angle at the farthest end of the survey area does not exceed 0.005m, which is less than the minimum actual distance of 0.5m between adjacent vertices of the boundary polygon. This avoids geometric misjudgment caused by the ray crossing adjacent non-intersecting line segments due to the minute angular deflection. When the total number of geometric intersections is even, the corresponding triangle is determined to be outside the geometric boundary of the non-convex survey area and is topologically removed. When the total number of geometric intersections is odd, the corresponding triangle is determined to be inside the geometric boundary of the non-convex survey area and is retained. The two-dimensional topological structure of all retained internal triangles is mapped to the three-dimensional elevation space, and the corresponding geodetic water depth elevation value is configured for each grid vertex to generate a three-dimensional irregular triangular mesh surface model as a three-dimensional underwater digital elevation model.

[0049] In the process of multi-stage topographic evolution differential analysis and engineering quantity calculation, a three-dimensional topographic surface model of the same survey area in the reference period is obtained. The three-dimensional underwater digital elevation model constructed in the early stage, namely the three-dimensional irregular triangular mesh surface model, will be used. Compared with the baseline three-dimensional terrain surface model Perform collinear vertical projection difference operations to construct a spatially continuous difference field. When the difference field value is greater than zero, it is marked as a sediment deposition and backfilling area; when the difference field value is less than zero, it is marked as a water scour or excavation removal area. A spatial scour and deposition color distribution map is generated based on the difference field value gradient, and all three-dimensional micro-triangles contained within the specified non-convex engineering calculation area are extracted. For any number of A three-dimensional infinitesimal triangle, based on the planar coordinates of its three vertices. Calculate the horizontal projected area Obtain the excavation depth of each of the three vertices relative to the designed excavation reference plane. Calculate the earthwork volume of the infinitesimal triangular prism element. The earthwork volumes of all U infinitesimal triangular prism elements within the engineering calculation area are algebraically summed to output the total earthwork excavation volume for the area. .

[0050] During the dynamic extraction of cross-sectional morphology, the starting point of the section is specified in the three-dimensional irregular triangular mesh surface model. and the finish line A baseline profile segment is generated, and the spatial intersection sequence between the baseline profile segment and the boundary edges of each spatial triangle in the 3D triangular mesh is calculated. The spatial intersection sequence is arranged in linear ascending order according to the horizontal projection distance from the profile starting point A. A continuous curve of lateral elevation change of the waterway is generated by interpolation between adjacent spatial intersection points, and this continuous curve is overlapped and compared with the design section slope line to output the verification data of navigation bottom width and slope ratio. In the specific calculation of the spatial intersection sequence, the baseline profile segment is first projected onto a horizontal two-dimensional plane, and the boundary segments of the triangular mesh that intersect with the baseline profile segment in two dimensions are retrieved one by one, and the coordinates of the two-dimensional intersection points in the plane are calculated. Then, the three-dimensional geodetic elevations of the two endpoints of the intersecting boundary segments are obtained, and one-dimensional linear elevation interpolation is performed according to the length ratio of the two-dimensional intersection points along the boundary segments to calculate the three-dimensional spatial elevations corresponding to the two-dimensional intersection points. Based on this, a continuous sequence of intersection points with three-dimensional spatial coordinates is generated, so that the profile elevation along the triangular mesh is consistent with the three-dimensional spatial coordinates. Boundary continuous extraction; when overlapping and comparing and outputting verification data, extract the left and right spatial intersections of the continuous curve of the transverse elevation change of the waterway and the design navigation bottom elevation line in the design section, and calculate the horizontal projection distance between the two intersections as the measured navigation bottom width value; extract the point set of the continuous curve of the transverse elevation change of the waterway in the left and right design slope projection sections, use the least squares method to perform linear fitting to obtain the slope, and calculate the measured slope ratio value based on the reciprocal of the fitted slope; calculate the width deviation between the measured navigation bottom width and the design bottom width, compare the measured slope ratio with the design section slope ratio, and when the measured navigation bottom width is greater than or equal to the design bottom width and the width deviation is within the engineering allowable limit, and the measured slope ratios on both the left and right sides are not steeper than the design slope ratio, the cross-section dimensions are deemed to be qualified; otherwise, it is deemed to be unqualified due to under-excavation or excessive slope ratio, and the verification data containing the measured navigation bottom width value, the measured slope ratio value, and the qualification / disqualification conclusion is output.

[0051] Example 4: This example combines Figures 1 to 4 This paper describes a method for filtering and 3D modeling underwater terrain point clouds, such as... Figure 1As shown, in step S101, a three-dimensional discrete sounding point set of the underwater area is obtained. Sounding points within the same transmitted pulse cross section are extracted based on the acoustic pulse timing and instantaneous heading azimuth. The projection distance of the sounding point perpendicular to the heading is calculated as the abscissa, and the original elevation is used as the ordinate to generate a two-dimensional sounding profile spatial point array. In step S102, the two-dimensional sounding profile spatial point array is converted into a dual-channel geometric feature tensor containing local point density channels and height gradient channels. This tensor is input into a convolutional target detection network to predict the bounding box of abnormal noise. Sounding points falling within the bounding box of abnormal noise are removed, and a pure underwater terrain point set composed of pure underwater terrain points is output. In step S103, a regular spatial grid is divided. The unit extracts shallow points with positive relative elevation differences from the preset design datum elevation surface to form a shallow point subset. Based on the elevation mean and elevation standard deviation of the shallow point subset, the representative water depth value of the grid unit is calculated, generating a regularized water depth matrix containing regularized water depth feature points. In step S104, the outer contour of the survey area is extracted to construct a non-convex closed topological boundary loop. A planar triangular network is initialized with pure underwater topographic points and regularized water depth feature points. The non-convex closed topological boundary loop is inserted as a strong constraint line segment, and the intersecting diagonals are flipped to make them common edges. The outer triangles are eliminated according to the parity of the number of intersections of the centroid rays. The internal triangle topology is mapped to three-dimensional space to construct a three-dimensional underwater digital elevation model.

[0052] like Figure 2 As shown, the underwater 3D discrete sounding point set is transformed by spatial projection to generate a 2D sounding profile spatial point matrix. The 2D sounding profile spatial point matrix is ​​filtered by target detection to output a pure underwater topographic point set. The pure underwater topographic point set is statistically calculated using shallow point elevation features to generate a regularized water depth matrix. The regularized water depth matrix is ​​retrieved by threshold comparison to output spatial warning markers. At the same time, the pure underwater topographic point set and the regularized water depth matrix are used together as the input discrete node set. Combined with the strong constraint boundary conditions provided by the non-convex closed topological boundary ring, a 3D underwater digital elevation model is constructed. The 3D underwater digital elevation model is further used to solve the spatial intersection sequence and interpolate to generate a continuous curve of the lateral elevation change of the channel. The total earthwork excavation volume of the region is output by performing the accumulation and calculation of the micro-element triangular prism unit. The spatial scour and siltation color distribution map is output by vertical projection differential rendering.

[0053] like Figure 3 As shown, within the plane bathymetry data distribution space starting from value 0, which is composed of the lateral projection distance u / m (perpendicular to the heading) and the original water depth elevation v / m, the continuous and densely distributed real underwater topographic bathymetry points in space jointly delineate and connect to form the channel slope break region. Isolated clutter points that are detached from this continuous elevation connectivity are defined by the corresponding geometric boundaries and fall within the spatial coverage of the anomalous noise bounding box BmB shown in the figure.

[0054] like Figure 4As shown, the global technical goal of constructing a three-dimensional underwater digital elevation model relies on the interactive support of four underlying logical branches. Among them, the two-dimensional sounding profile spatial point array, which serves as the spatiotemporal reference, is formed by the interweaving and convergence of acoustic pulse timing, instantaneous heading azimuth, projection distance, and original elevation. The pure underwater topographic point cluster, which is responsible for data purification, is reconstructed by the hierarchical screening logic of local pixel units, dual-channel geometric feature tensors, and classification confidence scores. The regularized water depth matrix, which represents the extreme characteristics of water depth, is limited by the joint calculation of the preset design reference elevation surface, regular spatial grid units, shallow point subsets, and representative water depth values. The three-dimensional irregular triangular mesh surface model, which serves as the final surface carrier, is strictly controlled by spatial geometric rules such as non-convex closed topological boundary loops, strongly constrained line segments, and the parity of the number of intersections of centroid rays.

[0055] Example 5: Before the multibeam echo sounder is mounted on a surface measurement platform to collect underwater topographic data, the spatial relative installation deviation parameters between the positioning antenna and the geometric center of the acoustic transducer are obtained. Based on the measured sound velocity profile and tidal observation sequence of the measurement area, the beam refraction bending and instantaneous water level fluctuations are corrected. After corrections for roll, pitch, and heave attitudes to eliminate hull sway errors, an underwater three-dimensional discrete echo sounder point set in a geodetic coordinate system is calculated and generated. ,in For Gaussian projection plane coordinates, To determine the water depth and geodetic elevation, extract the first... The horizontal geometric center coordinates of the transducer when triggered by an acoustic emission pulse and instantaneous heading azimuth Establish a local two-dimensional sounding profile coordinate system u,v with the survey line heading as the normal and the direction perpendicular to the heading as the lateral direction, and calculate the lateral projection coordinates. and the original elevation As the longitudinal coordinate of the cross section In the local two-dimensional sounding profile coordinate system, a channel containing local point density is constructed at a resolution of 0.05m. With normalized height gradient channels The dual-channel two-dimensional geometric feature tensor is input into a single-stage convolutional object detection network to predict anomalous noise bounding boxes. and its classification confidence score Compare with preset judgment threshold Remove suspended debris and bubble disturbance noise from the water body that fall within the abnormal noise frame, and output a pure underwater topographic point cluster that retains the geometric continuity of micro-topographic slopes. .

[0056] Divide the horizontal two-dimensional projection plane into steps with size... Regular spatial grid cells Based on the preset design reference elevation surface Extract the shallow point subset with positive relative elevation difference. Arithmetic mean of shallow point elevations Standard deviation of shallow sample Combined with extreme value retention sensitivity coefficient ,in accordance with Calculate the representative water depth value of the grid cell; extract the ordered polyline nodes of the outer contour of the irregular berth and the channel turning area, and construct a non-convex closed topological boundary loop. Each boundary segment As a strongly constrained line segment, it is inserted into the initial standard Delaunay triangulation formed by pure underwater topographic points and regularized water depth feature points. The strongly constrained line segment is then flipped locally through diagonal rotation. Become a common edge of adjacent triangles; calculate the geometric centroid coordinates of each triangle on the horizontal projection plane, and emit horizontal rays along the positive direction of the horizontal coordinate axis. When the horizontal ray coincides with a vertex of the boundary polygon, apply a vector to the ray direction. Tiny angular perturbations of rad to eliminate degenerate singular intersections, statistical horizontal rays and non-convex closed topological boundary loops. The total number of intersections is calculated, and external triangles with an even number of intersections are removed, while internal triangles with an odd number of intersections are retained. The two-dimensional topological structure of the retained internal triangles is mapped to the three-dimensional elevation space and assigned geodetic water depth and elevation values ​​to generate a three-dimensional irregular triangular network surface model that fits the actual non-convex physical boundary, which serves as the three-dimensional underwater digital elevation model.

[0057] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for filtering and 3D modeling underwater terrain point clouds, characterized in that, Includes the following steps: Step S101: Obtain the set of three-dimensional discrete sounding points underwater in the test area. Extract the sounding points within the same transmitted pulse cross section based on the acoustic pulse timing and instantaneous heading azimuth. Calculate the projection distance of the sounding point perpendicular to the heading as the abscissa and the original elevation as the ordinate to generate a two-dimensional sounding profile spatial point array. Step S102: Convert the two-dimensional sounding profile spatial point array into a dual-channel geometric feature tensor containing local point density channels and height gradient channels. Input the convolutional target detection network to predict the abnormal noise bounding box, remove the sounding points that fall into the abnormal noise bounding box, and output a pure underwater terrain point agglomeration composed of pure underwater terrain points. Step S103: Divide the space into regular grid cells. Using the preset design reference elevation surface as a reference, extract shallow points with positive relative elevation differences to form a shallow point subset. Calculate the representative water depth value of the grid cell based on the elevation mean and elevation standard deviation of the shallow point subset, and generate a regularized water depth matrix containing regularized water depth feature points. Step S104: Extract the outer contour of the survey area to construct a non-convex closed topological boundary ring. Initialize and construct a planar triangular network using pure underwater topographic points and regularized water depth feature points. Insert the non-convex closed topological boundary ring as a strong constraint line segment and flip the intersecting diagonals to make them common edges. Eliminate the outer triangles based on the parity of the number of intersections of the centroid rays. Map the internal triangle topology to three-dimensional space to construct a three-dimensional underwater digital elevation model.

2. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, Step S101, generating the two-dimensional sounding profile spatial lattice, includes: Step S1011, determining the transducer in the... The horizontal center coordinates and instantaneous heading azimuth angle at the trigger of the first transmission pulse; Step S1012, for the first... For any discrete sounding point of each transmitted pulse cross section, the lateral projection distance of the sounding point in the direction perpendicular to the flight direction is calculated based on the horizontal center coordinate and the instantaneous azimuth angle and used as the abscissa of the profile, and the original elevation is used as the ordinate of the profile; in step S1013, a two-dimensional sounding profile spatial point array is generated on the basis of a single transmitted pulse cross section, and the two-dimensional sounding profile spatial point array corresponding to each transmitted pulse cross section is generated pulse by pulse along the flight path according to the transmission pulse time sequence for feeding into the convolutional target detection network.

3. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, Step S102, converting to a dual-channel geometric feature tensor containing a local point density channel and a height gradient channel, includes: Step S1021, dividing the local two-dimensional profile coordinate system into local pixel units with a spatial resolution of 0.05m; Step S1022, counting the number of discrete sounding points falling into each pixel unit to construct a first feature channel representing the local point density; Step S1023, calculating the first-order difference between the representative elevation of the point set inside each column of pixel units and the elevation of the adjacent column to construct a second feature channel representing the normalized height gradient; Step S1024, combining the first and second feature channels to generate a dual-channel two-dimensional geometric feature tensor, and inputting the dual-channel two-dimensional geometric feature tensor into a single-stage convolutional detection network containing a feature pyramid, outputting anomaly noise bounding boxes and their classification confidence scores; comparing the classification confidence scores with a preset judgment threshold, when the classification confidence score is greater than or equal to the preset judgment threshold, sounding points falling within the geometric range of the corresponding bounding box are removed, and the preset judgment threshold ranges from 0.80 to 0.

90.

4. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, Step S104 involves inserting a non-convex closed topological boundary ring as a strongly constrained line segment and flipping the intersecting diagonals to make them common edges. This includes: Step S1041, extracting the actual outer boundary line segments of the survey area to construct an ordered closed point sequence; Step S1042, using pure underwater topographic points and regularized water depth feature points as input discrete node sets to construct an initial standard Delaunay triangulation on a two-dimensional plane; Step S1043, traversing each strongly constrained line segment in the non-convex closed topological boundary ring one by one, searching for the triangle edges in the initial standard Delaunay triangulation that geometrically intersect with the current strongly constrained line segment, and performing local diagonal flipping on a pair of common triangles containing the intersecting edges until each strongly constrained line segment is completely embedded in the triangulation to become a common edge of adjacent triangles.

5. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, Step S104, which involves eliminating external triangles based on the parity of the number of intersections of the centroid rays, includes: Step S1044, calculating the geometric centroid coordinates of each triangle in the planar triangulation on the two-dimensional projection plane; Step S1045, emitting a horizontal ray from each geometric centroid along the positive direction of the horizontal coordinate axis, and traversing and counting the total number of geometric intersections between the horizontal ray and all line segments in the non-convex closed topological boundary ring; Step S1046, when the total number of geometric intersections is even, determining that the corresponding triangle is located outside the geometric boundary of the non-convex measurement area and eliminating it; when the total number of geometric intersections is odd, determining that the corresponding triangle is located inside the geometric boundary of the non-convex measurement area and retaining it. The two-dimensional topology of all retained internal triangles is mapped to three-dimensional space, and each grid vertex is assigned a geodetic water depth elevation value, generating a three-dimensional irregular triangular mesh surface model as a three-dimensional underwater digital elevation model.

6. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, After constructing the three-dimensional underwater digital elevation model, the process also includes a dynamic extraction step of cross-sectional morphology: Step S105, in the constructed three-dimensional underwater digital elevation model, specify the starting point and ending point of the profile to generate the profile baseline segment, solve the spatial intersection sequence of the profile baseline segment and the boundary edge of each spatial triangle in the three-dimensional triangular network, arrange the spatial intersection sequence in linear ascending order according to the horizontal projection distance from the starting point of the profile, and interpolate between the spatial intersections to generate a continuous curve of the transverse elevation change of the channel, and compare the continuous curve of the transverse elevation change of the channel with the design section slope line to output the verification data of the navigable bottom width and slope ratio.

7. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, After constructing the three-dimensional underwater digital elevation model, the process also includes a multi-stage topographic evolution differential analysis step: Step S106, obtain the baseline three-dimensional topographic surface model of the same survey area, and perform collinear vertical projection difference between the three-dimensional underwater digital elevation model constructed in the current stage and the baseline three-dimensional topographic surface model to construct a spatial continuous differential field reflecting the changes in scour and deposition. When the differential field value is greater than zero, it is marked as a sediment deposition and backfilling area. When the differential field value is less than zero, it is marked as a water flow scour or excavation removal area. Based on the differential value gradient, a spatial scour and deposition color patch distribution map is generated.

8. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, After constructing the three-dimensional underwater digital elevation model, the process also includes an engineering earthwork volume integration step: Step S107, extract all three-dimensional micro-triangles contained within the specified engineering calculation area. For any m-th three-dimensional micro-triangle, calculate the projected area of ​​the micro-triangle on the horizontal plane based on the planar coordinates of its three vertices, obtain the excavation depth of the three vertices relative to the design excavation reference plane, multiply the projected area by one-third of the sum of the excavation depths of the three vertices as the earthwork volume of the micro-triangular prism element, and algebraically accumulate the earthwork volumes of all micro-triangular prism elements within the specified engineering calculation area to obtain the total engineering earthwork excavation volume of the area.

9. The underwater terrain point cloud filtering and 3D modeling method according to claim 1, characterized in that, After generating the regularized water depth matrix, the navigation water depth safety early warning step is also included: Step S108, set the waterway maintenance water depth reference threshold, traverse and search for abnormal grid nodes in the regularized water depth matrix whose water depth values ​​are less than the waterway maintenance water depth reference threshold, extract the plane coordinates and water depth elevation of the abnormal grid nodes, compare the water depth elevation of the abnormal grid nodes with the preset design reference elevation surface to calculate the superelevation difference and plane influence range of the local navigation obstruction shallow points, and output spatial early warning marks indicating the location of the waterway under-dredging obstacle.

Images