Seafloor topography scanning method and system based on underwater robot

Abstract

The present application relates to the technical field of seabed topography mapping, and particularly relates to a seabed topography scanning method and system based on an underwater robot. The present application solves the technical problem that the prior art is difficult to accurately and efficiently locate the specific position of a hidden groove in a seabed pipe pile array under the premise of maintaining low energy consumption and large-scale sparse scanning. The method comprises: obtaining discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, and position information of rigid obstacles in the seabed topography; constructing a connected graph according to the discrete acoustic reflection intensity; determining an injection boundary region and an outflow boundary region based on the position information of the rigid obstacles; executing a maximum flow algorithm with the capacity value of each node in the connected graph as a constraint to determine a flow distribution matrix; and determining a navigation direction for guiding the underwater robot according to the flow distribution matrix. The present application is used in the scene of seabed topography mapping of an underwater robot.

Description

Technical Field

[0001] This invention relates to the field of seabed topographic mapping technology, specifically to a seabed topographic scanning method and system based on an underwater robot. Background Technology

[0002] With the large-scale construction of marine engineering projects such as photovoltaic power generation in shallow sea areas, numerous grid-like arrays of solid pipe piles are deployed on the seabed, supporting the floating photovoltaic modules. Natural ocean currents accelerate locally as they pass through the gaps between adjacent pipe piles, continuously eroding the seabed sediment. This causes the pipe pile bases to be hollowed out, forming depressions and trenches, seriously threatening the long-term stability of the pipe piles. Therefore, operation and maintenance units need to periodically scan the seabed topography of the pipe pile array area to determine the specific location of the trenches and monitor their elevation.

[0003] Currently, in marine engineering operations, underwater robots equipped with multibeam echo sounders typically scan the seabed along pre-set routes. To control battery consumption and complete a single inspection of a large-area pipe pile array, the underwater robot is often programmed to perform sparse scanning along parallel routes with large spacing. However, when the sonar beam crosses over a trench, the large grazing angle and multipath effects at the terrain edges severely interfere with the time-of-arrival (TOA) algorithm at the bottom of the sonar, resulting in the acquisition of a missing geometric elevation data of the trench bottom in the 3D point cloud. If the scanning frequency is forcibly increased or the route spacing is reduced to obtain the elevation of the deep trench, the underwater robot's limited battery power will be quickly depleted, making it impossible to complete the inspection of the entire area. In other words, existing technologies, while maintaining low-energy consumption and large-area sparse scanning, struggle to accurately and efficiently locate the specific positions of potentially hazardous trenches within a seabed pipe pile array. Summary of the Invention

[0004] To address the technical problem of existing methods struggling to accurately and efficiently locate hidden trenches in subsea pipe pile arrays while maintaining low-energy, large-area sparse scanning, the present invention aims to provide a seabed topographic scanning method and system based on an underwater robot. The specific technical solution adopted is as follows: In a first aspect, the present invention provides a seabed topography scanning method based on an underwater robot. The method includes: acquiring discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, and the location information of rigid obstacles in the seabed topography; the discrete acoustic reflection intensity is used to characterize the looseness of the seabed substrate; a connected graph is constructed based on the discrete acoustic reflection intensity; nodes in the connected graph have capacity values; the capacity values ​​are determined based on the acoustic reflection intensity, with nodes having higher capacity values ​​for lower acoustic reflection intensity; injection boundary regions and outflow boundary regions are determined based on the location information of the rigid obstacles; a maximum flow algorithm is executed with the capacity value of each node in the connected graph as a constraint to determine a flow distribution matrix; the flow distribution matrix is ​​used to characterize the degree of flow aggregation; and a navigation course is determined based on the flow distribution matrix to guide the underwater robot.

[0005] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: projecting discrete acoustic reflection intensity onto a two-dimensional grid to generate a grid average intensity matrix; each grid cell in the grid average intensity matrix corresponds to a node in a connected graph; normalizing the acoustic reflection intensity of each grid cell in the grid average intensity matrix according to the intensity extrema in the grid average intensity matrix to determine the normalization difference of each grid cell; converting the normalization difference of each grid cell into a capacity value of each node in the connected graph according to a preset nonlinear mapping relationship; the nonlinear mapping relationship allows nodes with smaller normalization differences to obtain higher capacity values; connecting each node with a capacity value according to its adjacency relationship, and assigning the capacity value of the connecting node to the connecting edge to generate a connected graph.

[0006] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: identifying adjacent first rigid obstacles and second rigid obstacles; the first rigid obstacle and the second rigid obstacle are a pair of rigid obstacles with a horizontal distance of less than a preset threshold within the seabed operating area; using the position of the first rigid obstacle as a reference, and according to the fitting radius of the first rigid obstacle and a preset search range, delineating a first peripheral area as an injection boundary area; using the position of the second rigid obstacle as a reference, and according to the fitting radius of the second rigid obstacle and a preset search range, delineating a second peripheral area as an outflow boundary area.

[0007] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: establishing the source and sink of the maximum flow algorithm; connecting the source to each node within the injection boundary region and connecting each node within the outflow boundary region to the sink; running a preset maximum flow algorithm with the source and sink as inputs and the capacity value of each node as a constraint; monitoring the number of algorithm iterations and the number of active nodes in the network, and determining that the algorithm has converged and stopping the push when the number of iterations reaches a preset upper limit or there are no active nodes in the network; generating a flow distribution matrix based on the cumulative flow value of each node in the converged state.

[0008] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: performing skeleton extraction processing on the flow distribution matrix to obtain a grid sequence with a width of a single grid and nodes sequentially adjacent; each grid in the grid sequence has a grid index; inversely mapping the grid index in the grid sequence to a two-dimensional coordinate set of the seabed operating area; performing linear fitting on the two-dimensional coordinate set to generate a navigation heading.

[0009] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: determining the geometric normal direction of the line connecting the first rigid obstacle and the second rigid obstacle; traversing each node in the flow distribution matrix, taking the currently traversed node as the current node, establishing a comparison neighborhood window of a preset size for the current node, and determining the nodes within the comparison neighborhood window that intersect with the spatial projection of the geometric normal direction as nodes to be compared; comparing the flow value of the current node with the flow values ​​of all nodes to be compared, and if the flow value of the current node is less than the flow value of any node to be compared, then clearing the flow value of the current node to zero; after traversing all nodes, connecting each node with a non-zero flow value in series according to the adjacency relationship to generate a grid sequence with a width of a single grid and nodes sequentially adjacent to each other.

[0010] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: obtaining the reference origin coordinates and grid edge length when constructing the connected graph; for each grid cell in the grid sequence, determining the two-dimensional coordinates of the geometric center of each grid cell based on its horizontal and vertical indices, the reference origin coordinates, and the grid edge length; using the least squares method to linearly fit the two-dimensional coordinates of the geometric center of each grid cell, and determining the navigation heading based on the angle between the fitted line and the reference direction.

[0011] In conjunction with the first aspect mentioned above, in one possible implementation, the method specifically includes: receiving the echo signal from a multibeam echo sounder and generating a basic three-dimensional point cloud set; each point in the basic three-dimensional point cloud set contains three-dimensional spatial coordinates and discrete acoustic reflection intensity; extracting the three-dimensional point cloud of a rigid obstacle from the basic three-dimensional point cloud set using a geometric fitting algorithm, and determining the horizontal center coordinates and fitting radius of the rigid obstacle; and determining the horizontal center coordinates and fitting radius of the rigid obstacle as the position information of the rigid obstacle.

[0012] In conjunction with the first aspect mentioned above, in one possible implementation, the method further includes: generating control commands based on the navigation heading and the position coordinates of the first rigid obstacle and the second rigid obstacle; the control commands are used to adjust the underwater robot's heading to the navigation heading and control the underwater robot to fly along the navigation heading at a preset altitude, perform vertical detection, and acquire a high-frequency encrypted point cloud; receiving the high-frequency encrypted point cloud acquired by the underwater robot based on the control commands, and superimposing and fusing the high-frequency encrypted point cloud with the basic three-dimensional point cloud set in the same coordinate system to generate a three-dimensional seabed topographic model.

[0013] Secondly, the present invention provides a seabed topography scanning system based on an underwater robot, used to implement the seabed topography scanning method based on an underwater robot described in the first aspect above; the system includes: an acquisition module, used to acquire discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, and the position information of rigid obstacles in the seabed topography; the discrete acoustic reflection intensity is used to characterize the looseness of the seabed sediment; a connected graph construction module, used to construct a connected graph based on the discrete acoustic reflection intensity; the nodes in the connected graph have capacity values; the capacity values ​​are determined based on the acoustic reflection intensity, and the lower the acoustic reflection intensity, the higher the capacity value of the node; a boundary determination module, used to determine the injection boundary region and the outflow boundary region based on the position information of the rigid obstacles; a deduction module, used to execute a maximum flow algorithm with the capacity value of each node in the connected graph as a constraint to determine the flow distribution matrix; the flow distribution matrix is ​​used to characterize the degree of flow aggregation; and a navigation determination module, used to determine the navigation course for guiding the underwater robot based on the flow distribution matrix.

[0014] The present invention has the following beneficial effects: This invention acquires discrete acoustic reflection intensity and rigid obstacle location information obtained from sparse scanning by an underwater robot. Reflection intensity is used to characterize the looseness of the seabed substrate, thereby constructing a connected graph where node capacity is negatively correlated with reflection intensity. Based on this, injection and outflow boundaries are determined according to obstacle locations. A maximum flow algorithm is executed with node capacity as a constraint to deduce the flow distribution matrix, and finally, the navigation heading is determined based on this matrix. This series of steps transforms the traditional geometric mapping problem relying on elevation data into a topology optimization problem based on acoustic intensity. While maintaining large-spacing sparse scanning (thus avoiding rapid battery depletion), it objectively locates the convergence direction of scour trenches through graph theory maximum flow deduction, bypassing the elevation data failure dilemma caused by large grazing angles and multipath effects. This achieves low-energy, high-efficiency spatial location identification of hidden seabed trenches, providing an accurate navigation benchmark for subsequent local resurveys. Thus, it solves the technical problem of existing methods struggling to accurately and efficiently locate the specific location of hidden trenches in subsea pipe pile arrays while maintaining low-energy, large-area sparse scanning. Attached Figure Description

[0015] To more clearly illustrate the technical solutions and advantages 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.

[0016] Figure 1 This is a flowchart illustrating a seabed topography scanning method based on an underwater robot, provided in one embodiment of the present invention. Figure 2 This is a schematic diagram of the system architecture of a seabed topography scanning system based on an underwater robot, provided in one embodiment of the present invention. Detailed Implementation

[0017] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of the seabed topography scanning method and system based on an underwater robot proposed according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0018] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0019] The specific solution of the seabed topography scanning method and system based on underwater robots provided by the present invention will be described in detail below with reference to the accompanying drawings.

[0020] Please see Figure 1 The diagram illustrates a flowchart of a seabed topography scanning method based on an underwater robot, according to an embodiment of the present invention. The method includes the following steps S101-S106, which will be described in detail below.

[0021] S101. Obtain the discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, as well as the location information of rigid obstacles in the seabed topography.

[0022] In one possible implementation, after entering the water, the underwater robot establishes a coordinate system reference origin and moves along a pre-set parallel route with large spacing. The robot's multibeam sonar emits fan-shaped sound waves towards the seabed, receiving echo signals from the seabed sediment and underwater structures, and analyzing them to generate a basic three-dimensional point cloud set. Each discrete data point in the basic three-dimensional point cloud set contains three-dimensional spatial coordinates and discrete acoustic reflection intensity, where the discrete acoustic reflection intensity is measured in decibels and is used to characterize the looseness of the seabed sediment at the corresponding location; the lower the value, the more severe the erosion of sediment and the softer the seabed sediment.

[0023] For example, discrete points conforming to the geometric characteristics of a vertical cylinder are extracted from the basic 3D point cloud set using the Random Sample Consensus (RANSAC) geometric fitting algorithm and aggregated into a solid pile foundation point cloud. For each independent cylinder, the intersection of its central axis and the horizontal reference plane of the nautical chart is calculated to generate the horizontal center coordinates of the pile foundation; simultaneously, the fitting radius parameter of the cylinder is extracted according to the cylinder fitting equation. The horizontal center coordinates and fitting radius of each pile foundation obtained in this way constitute the positional information of rigid obstacles in the seabed topography. At the same time, the discrete acoustic reflection intensity carried by each discrete data point in the basic 3D point cloud set, along with its corresponding horizontal coordinates, are saved as the raw data used in subsequent steps.

[0024] S102. Construct a connected graph based on the discrete acoustic reflection intensity.

[0025] In this connected graph, each node has a capacity value; the capacity value is determined based on the acoustic reflection intensity, and the lower the acoustic reflection intensity, the higher the capacity value of the node.

[0026] In one possible implementation, the discrete acoustic reflection intensity is first projected onto a two-dimensional planar grid, with each grid cell corresponding to a node in the connected graph. Then, based on the acoustic reflection intensity distribution of all grid cells, the global intensity extremum is determined, and the intensity of each grid cell is normalized to obtain a dimensionless normalized difference. Next, using a preset nonlinear mapping relationship (e.g., a negative exponential function), the normalized difference is converted into a capacity value for each node. This mapping relationship ensures that nodes with smaller normalized differences (i.e., lower original acoustic reflection intensity and softer substrate) receive higher capacity values. Finally, the nodes are connected according to the adjacency relationship between grid cells, and each connecting edge is assigned a capacity value to the node it connects to, thereby generating a complete connected graph.

[0027] S103. Based on the location information of rigid obstacles, determine the injection boundary region and the outflow boundary region.

[0028] In one possible implementation, the horizontal spacing between all pile foundations is first calculated. Pairs of adjacent pile foundations with a spacing less than a preset threshold (e.g., 1.5 times the design spacing of the pipe piles) are extracted and treated as a pair of adjacent rigid barriers, designated as the first rigid barrier and the second rigid barrier, respectively. Then, using the horizontal center coordinates of the first rigid barrier as the center, and its fitted radius plus a preset search range (e.g., the side length of a grid) as the radius, an annular perimeter region is delineated. All grids within this region are merged into an injection node set; this region is the injection boundary region. Similarly, using the horizontal center coordinates of the second rigid barrier as a reference, a second perimeter region is delineated according to the same rules, serving as the outflow boundary region.

[0029] S104. Using the capacity value of each node in the connected graph as a constraint, execute the maximum flow algorithm to determine the flow distribution matrix.

[0030] The flow distribution matrix is ​​used to characterize the degree of flow aggregation.

[0031] One possible implementation involves taking the constructed connected graph and its node capacities as input, directly running a pre-defined maximum flow algorithm, using the capacity of each node as edge capacity constraints, to deduce the distribution of traffic in the network. After the algorithm completes, the total traffic flowing through each node is counted, generating a traffic distribution matrix that corresponds one-to-one with the nodes in the connected graph. Each element in this matrix represents the degree of traffic concentration at the corresponding node location.

[0032] For example, the maximum flow algorithm receives flow from the nodes corresponding to the injection boundary region, propagates it node by node within the connected graph according to capacity constraints, and finally flows out from the nodes corresponding to the outflow boundary region. High-capacity nodes (corresponding to loose substrate) attract more flow, while low-capacity nodes (corresponding to dense substrate) restrict flow. After the algorithm completes, the high-value regions in the flow distribution matrix indicate the locations where flow accumulates, corresponding to the orientation of the seabed scour channels.

[0033] S105. Based on the flow distribution matrix, determine the navigation course used to guide the underwater robot.

[0034] In one possible implementation, the high-value regions in the flow distribution matrix are distributed in a continuous band, and the spatial extension direction of this band corresponds to the direction of the trenches formed by the scouring of seabed sediment. Spatial analysis is performed on the band to determine the extension direction of its centerline, and then the angle between this direction and a reference direction (such as geographic north or the robot's initial heading) is calculated. This angle is used as the navigation heading to guide the underwater robot.

[0035] For example, each node in the flow distribution matrix corresponds to a cumulative flow value. High-value nodes are arranged continuously in space, forming a flow band with a certain width. First, the geometric center line of the flow band is identified to eliminate lateral width interference, obtaining a continuous trajectory with a single pixel width. Then, the discrete grid positions on this trajectory are restored to the real geographic coordinates within the seabed operating area. Finally, linear fitting is performed on the restored coordinate set to obtain an optimal straight line, and the direction angle of this line is the navigation heading. After receiving this heading, the underwater robot can adjust its attitude and fly stably along the trench.

[0036] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment obtains discrete acoustic reflection intensity and rigid obstacle location information obtained by sparse scanning of an underwater robot, and uses reflection intensity to characterize the looseness of the seabed substrate. Then, a connected graph with a negative correlation between node capacity and reflection intensity is constructed. Based on this, the injection boundary and outflow boundary are determined according to the obstacle location. The maximum flow algorithm is executed with node capacity as a constraint to deduce the flow distribution matrix, and finally, the navigation heading is determined based on this matrix. This series of steps transforms the traditional geometric mapping problem relying on elevation data into a topology optimization problem based on acoustic intensity. While maintaining large-spacing sparse scanning (thus avoiding rapid battery depletion), it can objectively locate the aggregation direction of scour trenches through graph theory maximum flow deduction, bypassing the elevation data failure dilemma caused by large grazing angles and multipath effects. This achieves low-energy, high-efficiency spatial location identification of hidden trenches on the seabed, providing an accurate navigation benchmark for subsequent local supplementary surveys. Thus, it solves the technical problem that existing methods, while maintaining low-energy, large-area sparse scanning, are unable to accurately and efficiently locate the specific location of hidden trenches in seabed pipe pile arrays.

[0037] In one possible implementation, the process of constructing a connected graph based on discrete acoustic reflection intensity can be specifically implemented through the following steps S201-S204, which will be explained in detail below.

[0038] S201. Project the discrete acoustic reflection intensity onto a two-dimensional grid to generate a grid-averaged intensity matrix.

[0039] In one possible implementation, the side length of the square grid is first determined based on the single-beam footprint diameter of the multibeam bathysound at the current operating water depth. Then, a reference origin for the two-dimensional grid coordinate system is established, forcibly coinciding with the reference origin of the three-dimensional absolute coordinate system established when the underwater robot enters the water. For each discrete seabed data point (point cloud excluding rigid obstacles), its horizontal coordinates are extracted, and its horizontal and vertical indices are calculated according to the grid division formula. The calculation method is: subtract the reference origin coordinates from the horizontal coordinates, divide by the grid side length, and then round down. Rounding down means taking the largest integer not greater than the value in parentheses.

[0040] For example, if a discrete data point has a horizontal index of 3 and a vertical index of 5, its reflection intensity is stored in the grid register at index (3, 5). After all discrete points on the natural seabed have been projected, each grid cell is traversed, the number of reflection intensity samples stored in its register is counted, and all samples are summed. The summation result is divided by the number of samples to obtain the arithmetic mean of the acoustic reflection intensity of that grid cell, which is denoted as an element of the grid average intensity matrix. For empty grid cells where no data points fall due to sparse scanning, no value is assigned yet, to be reserved for subsequent interpolation processing. The resulting two-dimensional grid average intensity matrix has each element representing the average acoustic reflection intensity of the seabed area covered by the corresponding grid cell, reflecting the looseness of the seabed sediment in that area.

[0041] For example, the first i Horizontal coordinates of data points Horizontal index With vertical index Satisfy the following formula:

[0042]

[0043] in, This is the horizontal index of the grid, representing the column number of the data point in the two-dimensional grid after projection; This is the vertical index of the grid, representing the row number of the data point in the two-dimensional grid after projection; This is the floor function, which rounds down to the nearest integer less than or equal to the value within the parentheses. and The origin of the two-dimensional grid coordinate system is the reference point, which coincides with the horizontal coordinate of the underwater robot's entry point into the water. The side length of the square grid is equal to the diameter of the single-beam footprint of the multibeam echo sounder at the current water depth.

[0044] To calculate the lateral offset distance of a data point relative to the reference origin, a positive value indicates that the point is east of the origin, and a negative value indicates that it is west of the origin; To calculate the longitudinal offset distance of a data point relative to the reference origin, a positive value indicates that the point is north of the origin, and a negative value indicates that it is south of the origin. Divide the offset distance by the grid side length to obtain the relative position in units of the grid side length, and round it to obtain the integer index of the grid to which the point belongs. and It performs independent horizontal and vertical operations, mapping any point on a two-dimensional plane to a unique grid cell.

[0045] S202. Based on the intensity extreme values ​​in the grid average intensity matrix, normalize the acoustic reflection intensity of each grid cell in the grid average intensity matrix to determine the normalization difference of each grid cell.

[0046] In one possible implementation, each grid cell in the average intensity matrix is ​​traversed. During the traversal, the intensity value of each grid cell is compared one by one, and the maximum and minimum values ​​encountered are recorded. After the traversal is completed, the recorded values ​​are the global maximum and global minimum values. For each grid cell, its acoustic reflection intensity is subtracted from the global minimum value to obtain the difference. This difference is then divided by the difference between the global maximum and global minimum values ​​(to prevent the denominator from being zero, a very small positive parameter adjustment coefficient with the same dimensions as the difference, such as 0.001, can be added). The result is the normalized difference degree of that grid cell. This normalized difference degree quantitatively characterizes how close the current grid cell's seabed looseness is to the loosest seabed in the seabed operating area. The closer the value is to 0, the closer the seabed is to the loosest state; the closer the value is to 1, the closer it is to the densest state.

[0047] S203. Based on the preset nonlinear mapping relationship, the normalized difference of each grid cell is converted into the capacity value of each node in the connected graph.

[0048] In one possible implementation, the capacity value characterizes the upper limit of flow allowed through a node. The transformation rule assigns a higher capacity value to nodes corresponding to grid cells with smaller normalized difference. Specifically, since smaller normalized difference represents lower original acoustic reflection intensity and looser substrate, which is more easily eroded and penetrated by water flow, nodes corresponding to such grid cells need to be assigned higher capacity values ​​to attract more flow aggregation in the subsequent maximum flow algorithm. The nonlinear mapping relationship uses a monotonically decreasing exponential function, which maps the normalized difference to a capacity value between 0 and 1. When the normalized difference approaches 0, the capacity value approaches 1, and when the normalized difference approaches 1, the capacity value approaches 0. A preset intensity sensitivity coefficient is introduced during the mapping process to control the sensitivity of the mapping: the larger the sensitivity coefficient, the more drastic the change in capacity value will occur due to a small change in the normalized difference; the smaller the sensitivity coefficient, the smoother the mapping curve.

[0049] For example, the specific conversion formula for converting normalized variance to capacity values ​​is as follows:

[0050] in, For grid The node capacity value is dimensionless. Its range is (0, 1], representing the maximum capacity of the grid to allow virtual fluid to pass through. The larger the value, the softer the seabed sediment at that location, and the easier it is to be eroded and hollowed out by the water flow. It is an exponential function with the natural constant e as its base (natural exponential function), used to convert negative linear combinations into positive capacity values; The preset intensity sensitivity coefficient (e.g., 3) is based on the seabed sediment characteristics of the operating area. Specifically, it is pre-calibrated according to the median particle size of the sediment. The larger the intensity sensitivity coefficient, the steeper the slope of the mapping curve. Small changes in the normalized difference will lead to more drastic changes in the capacity value, thus amplifying the differences in the looseness of the seabed sediment between different grids. The smaller the intensity sensitivity coefficient, the flatter the mapping curve, the less sensitive the capacity value is to changes in the normalized difference, and the smaller the difference in capacity value between each grid. For grid The average acoustic reflection intensity within the grid represents the density of the sediment in the corresponding seabed area. The lower the value, the looser the sediment (and the more severe the erosion). This represents the minimum acoustic reflection intensity of all grids within the seabed operating area. This represents the maximum acoustic reflection intensity of all grids within the seabed operating area; These are parameter tuning coefficients, and their values ​​should be extremely small positive numbers (e.g., 0.001) to avoid denominators of 0.

[0051] molecular Calculate the difference between the current grid strength and the global minimum strength to characterize the degree of deviation of the grid relative to the softest substrate; Calculate the global intensity range (maximum minus minimum), and add the parameter tuning coefficients. Ensure the denominator is not zero; Taking the negative sign and multiplying it by the sensitivity coefficient yields the exponent of the exponential part. The negative sign reverses the positive-negative relationship: the larger the original score (density), the more negative the exponential part becomes, resulting in a smaller exponential function output. Therefore, the final capacity value is negatively correlated with the normalized difference, that is, negatively correlated with the density of the substrate and positively correlated with the softness. This value represents the ability of a certain grid location on the seabed to conduct virtual fluid. The higher the value, the softer the seabed at that location, and the easier it is for the water flow to erode and penetrate it, thus attracting more flow accumulation in the subsequent maximum flow algorithm. This value reflects the vulnerability of the silt being eroded.

[0052] For example, in sea areas with fine sandy bottoms, the small particle size of the sediment makes the acoustic reflection intensity highly sensitive to changes in the looseness of the bottom, with a significant difference in intensity between soft and dense bottoms. In this case, setting the intensity sensitivity coefficient to a larger value (e.g., 3.0 to 5.0) can fully utilize this difference, making the capacity value close to 1 in the loose region and close to 0 in the dense region, thus forming a clear flow concentration zone in the maximum flow projection. In sea areas with coarse sand or gravel bottoms, the larger bottom particles result in a smaller variation in acoustic reflection intensity with the degree of looseness. Using an excessively large sensitivity coefficient can lead to noise amplification and increased misjudgment. In this case, setting the intensity sensitivity coefficient to a smaller value (e.g., 2.0 to 3.0) makes the mapping curve smoother and avoids drastic fluctuations in capacity values ​​due to small measurement errors. For mixed bottoms with uneven particle size distribution, the coefficient value that best matches the projected scour location can be selected through on-site sampling or historical data calibration.

[0053] S204. Connect each node with a capacity value according to its adjacency relationship, and assign the capacity value of the connected node to the connecting edge to generate a connected graph.

[0054] In one possible implementation, each node corresponds to a grid cell in a two-dimensional mesh. These grid cells have spatial adjacency relationships, including horizontal, vertical, and diagonal adjacency, i.e., eight-neighborhood adjacency. For each node, connections are established between it and nodes in its eight neighboring directions. When establishing these connections, each edge is assigned a capacity value equal to the capacity value of the node it points to. In other words, the capacity of an edge is determined by the capacity value of its endpoint (the target node), not the capacity value of its originating node. In this way, traffic flowing along edges in the network is limited by the capacity value of the target node.

[0055] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment further projects the discrete acoustic reflection intensity onto a two-dimensional grid to generate a grid average intensity matrix, uses the global intensity extremum to normalize the intensity of each grid to obtain a normalized difference degree, and then converts the normalized difference degree into a node capacity value according to a preset nonlinear mapping relationship. Finally, the nodes are connected according to the adjacency relationship and the edges are assigned capacity values ​​to generate a connected graph. Through gridded projection and normalization processing, the dimensional differences caused by different seabed operating areas and different sonar ranges are eliminated, making the extrapolation process universal; and the nonlinear mapping (such as negative exponential transformation) successfully corresponds low reflection intensity (soft bottom) to high capacity and high reflection intensity (dense bottom) to low capacity, realizing the physical property reversal from static acoustic parameters to the capacity boundary of the graph theory network, providing the maximum flow algorithm with constraints that conform to the hydrodynamic scour law, and ensuring that the extrapolation results can truly reflect the continuous fragile zone of sediment erosion.

[0056] In one possible implementation, the process of determining the injection boundary region and the outflow boundary region based on the location information of the rigid obstacle can be specifically implemented through the following S301-S303, which will be described in detail below.

[0057] S301. Determine the adjacent first rigid obstacle and second rigid obstacle.

[0058] Among them, the first rigid obstacle and the second rigid obstacle are a pair of rigid obstacles in the seabed operation area with a horizontal distance of less than a preset threshold.

[0059] In one possible implementation, the horizontal center coordinates of each rigid obstacle are first read, and the horizontal straight-line distance between any two obstacles is calculated. Then, the calculated distance is compared with a preset threshold, which is determined based on the design spacing of the pile array, typically 1.2 to 1.8 times the design spacing, to ensure that only actually adjacent pile pairs with a potential for scour are included in subsequent processing. When the horizontal distance between two obstacles is less than the preset threshold, they are determined to be an adjacent obstacle pair, with one designated as the first rigid obstacle and the other as the second rigid obstacle. The roles of the first and second rigid obstacles are physically symmetrical and can be interchanged according to the water flow direction or actual simulation needs, but each pair of adjacent piles will independently participate in the subsequent boundary area delineation.

[0060] S302. Based on the position of the first rigid obstacle, and according to the fitting radius of the first rigid obstacle and the preset search range, a first surrounding area is defined as the injection boundary area.

[0061] In one possible implementation, the horizontal center coordinates and fitting radius of the first rigid obstacle are first obtained. The fitting radius is directly provided by a preceding geometric fitting algorithm (such as RANSAC cylindrical fitting), reflecting the actual physical dimensions of the obstacle. Then, a preset search range is determined, typically a grid side length (i.e., the diameter of the sonar footprint), to add a buffer zone outside the fitting radius, covering the substrate area around the pile foundation that may be affected by water flow. A circular region (approximately a polygonal region in the actual meshed representation) is delineated with the horizontal center coordinates of the first rigid obstacle as the center and the sum of the fitting radius and the preset search range as the radius. All grid cells covered by this circular region are identified as the first peripheral region. This first peripheral region is the injection boundary region, representing the entry point for flow into the connected graph network from this region.

[0062] Optionally, in practical engineering applications, the preset search range can be appropriately adjusted according to factors such as the diameter of the pile foundation, the water flow velocity, and the characteristics of the seabed. For example, for large-diameter pile foundations or strong-current sea areas, the influence range of the water flow around the pile is larger, and the search range can be appropriately expanded to 1.5 or 2 times the grid side length; for small-diameter pile foundations or weak-current sea areas, it can be reduced to 0.5 times the grid side length. This disclosure does not limit this.

[0063] S303. Based on the position of the second rigid obstacle, and according to the fitting radius of the second rigid obstacle and the preset search range, delineate the second surrounding area as the outflow boundary area.

[0064] In one possible implementation, the horizontal center coordinates and fitting radius of the second rigid obstacle are first obtained. The fitting radius is also given by the preceding geometric fitting algorithm. The preset search range is consistent with step S302, typically taking one grid side length. A circular region is delineated with the horizontal center coordinates of the second rigid obstacle as the center and the sum of the fitting radius and the preset search range as the radius. All grid cells covered by this circular region are identified as the second peripheral region. This second peripheral region is the outflow boundary region, representing the exit region where flow exits from the connected graph network and converges to the sink.

[0065] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment determines adjacent rigid obstacle pairs (i.e., adjacent pipe piles) with a horizontal spacing of less than a preset threshold within the seabed operating area. Based on the positions of these two obstacles and according to their fitting radii and preset search ranges, the surrounding areas are delineated as injection boundary areas and outflow boundary areas. By using the actual rigid obstacles on the seabed as physical anchor points for the start and end points of scour, the gaps between pipe piles in the hydrodynamic environment are transformed into source and sink boundaries in graph theory algorithms. This allows the maximum flow extrapolation to naturally focus on the narrow areas most prone to local scour. Simultaneously, the introduction of fitting radii and search ranges effectively eliminates the rigid body regions of the obstacles themselves, avoiding invalid nodes in the flow injection, improving the specificity and computational efficiency of the extrapolation, and ensuring that the high-value areas of the flow distribution matrix accurately correspond to the trench locations with the softest seabed sediment between the pipe piles.

[0066] In one possible implementation, the process of using the capacity value of each node in the connected graph as a constraint to execute the maximum flow algorithm and determine the flow distribution matrix can be specifically implemented through the following S401-S405, which will be explained in detail below.

[0067] S401. Establish the source and sink nodes for the maximum flow algorithm.

[0068] In one possible implementation, using the defined injection and outflow boundary regions as input, two special nodes are established outside the connected graph, serving as the source (flow injection point) and sink (flow outflow point) of the maximum flow algorithm, respectively. Specifically, a new node is created outside the connected graph as the source. This source does not belong to any existing node in the connected graph; its role is to serve as the sole origin of the flow. Simultaneously, another new node is created outside the connected graph as the sink, serving as the sole final destination of the flow. The source and sink nodes themselves do not possess capacity properties; they are only used to define the starting and ending boundaries of the flow in the maximum flow algorithm.

[0069] It is understandable that the source and sink nodes themselves do not possess capacity attributes; they merely serve as boundary nodes for flow injection and outflow. In this way, multiple nodes that were originally scattered within the injection boundary region are unified into a common source node, and multiple nodes that were originally scattered within the outflow boundary region are unified into a common sink node, enabling the maximum flow algorithm to perform deductions between a single source and sink node.

[0070] S402. Connect the source point to each node within the injection boundary region, and connect each node within the outflow boundary region to the sink point.

[0071] One possible implementation involves iterating through all nodes within the injection boundary region. For each node, a directed edge is created from the source node to that node, with an infinite capacity (or no capacity limit) to ensure unimpeded flow from the source node to any node within the injection boundary region. Similarly, all nodes within the outflow boundary region are iterated through. For each node, a directed edge is created from that node to the sink node, also with an infinite capacity, to ensure unimpeded flow from any node within the outflow boundary region to the sink node.

[0072] S403. Using the source and sink as inputs and the capacity of each node as a constraint, run the preset maximum flow algorithm.

[0073] In one possible implementation, the established source and sink nodes, the constructed connected graph, and the capacity of each node are taken as input. A pre-defined maximum flow algorithm is run to calculate the maximum flow value from the source node, through the connected graph network, and finally to the sink node, and the flow distribution on each node after algorithm convergence is obtained. Specifically, the capacity of each edge in the connected graph is set to the capacity value of the node pointed to by the edge (i.e., the capacity of the target node). The source node is used as the flow source, and the sink node is used as the flow destination. The pre-defined maximum flow algorithm is called for iterative calculation. During the algorithm's operation, flow is continuously pushed until the total flow from the source node to the sink node can no longer be increased.

[0074] For example, a pre-defined maximum flow algorithm can employ a push and relabel algorithm. Starting from the source, the algorithm pushes traffic to nodes within the injection boundary region connected to the source, and then, based on the capacity constraints of each edge, pushes traffic node by node towards the sink. When all paths from the source to the sink reach their capacity limits, the algorithm converges, outputting a set of flow values ​​corresponding to the edges of the connected graph. By accumulating the inflow traffic to each node, the total flow value at each node can be obtained; these flow values ​​reflect the degree of traffic aggregation in the network.

[0075] S404. Monitor the number of algorithm iterations and the number of active nodes in the network. When the number of iterations reaches the preset limit or there are no active nodes in the network, determine that the algorithm has converged and stop pushing.

[0076] In one possible implementation, the number of iterations and the number of active nodes are monitored in real time during the execution of the maximum flow algorithm. When the number of iterations reaches a preset upper limit, or when there are no more active nodes in the network, the algorithm is considered to have converged, and traffic pushing is immediately stopped. The preset upper limit can be set based on the total number of network nodes, for example, taking the square of the total number of nodes (the logic behind this setting is that the worst-case time complexity of the push and remark maximum flow algorithm is O(V²E) (where V is the number of nodes and E is the number of edges). For a typical grid-structured connected graph, the number of edges E and the number of nodes V are linearly related, so the number of iterations in the worst case is approximately on the order of V³. Taking V² as the preset upper limit can effectively identify and terminate abnormal situations while ensuring the normal operation of the algorithm); active nodes refer to nodes that still have remaining traffic to push or are highly marked as needing updates.

[0077] S405. Generate a flow distribution matrix based on the cumulative flow value of each node under convergence.

[0078] In one possible implementation, for each node in the connected graph, the total traffic received, forwarded, and retained by that node after the algorithm converges is recorded; this total is the traffic accumulation value for that node. The traffic accumulation values ​​of all nodes are then filled into the corresponding positions of a two-dimensional matrix according to their corresponding grid cell positions (horizontal and vertical indices) to form a traffic distribution matrix.

[0079] For example, assuming the connected graph corresponds to a 100x100 grid matrix, the flow distribution matrix is ​​also a 100x100 two-dimensional array. The matrix element located in the i-th row and j-th column has a value equal to the accumulated flow value of the corresponding node in the connected graph. For nodes that have never received any flow, their accumulated value is zero. The resulting flow distribution matrix shows high-value regions corresponding to areas of concentrated flow, indicating continuous areas with loose substrate and easy erosion.

[0080] The technical solution provided by the above embodiments can bring at least the following beneficial effects: After establishing the source and sink of the maximum flow algorithm, this embodiment connects the source to each node in the injection boundary region and each node in the outflow boundary region to the sink. The preset maximum flow algorithm is run with the capacity of each node as a constraint, and the number of algorithm iterations and the number of active nodes in the network are monitored. When the number of iterations reaches the preset upper limit or there are no active nodes in the network, convergence is determined and the push stops. Finally, a flow distribution matrix is ​​generated based on the cumulative flow value of each node under the convergence state. By introducing the iteration upper limit and active node monitoring mechanism, infinite loops or high-remarking deadlocks under extreme terrain (such as two piles blocked by hard rock veins) are effectively prevented, ensuring stable convergence of the algorithm. The statistical analysis of the cumulative flow value of each node after convergence allows the high-value areas in the distribution matrix to objectively present the compression and aggregation pattern of flow under limited capacity, providing a high-quality data foundation for subsequent skeleton extraction and demonstrating the quantitative indication capability of the most severely eroded continuous fragile zone.

[0081] In one possible implementation, the process of determining the navigation course for guiding the underwater robot based on the flow distribution matrix can be specifically implemented through the following steps S501-S503, which will be described in detail below.

[0082] S501. Perform skeleton extraction processing on the flow distribution matrix to obtain a grid sequence with a width of a single grid and nodes that are sequentially adjacent.

[0083] Each grid in the grid sequence has a grid index.

[0084] In one possible implementation, the goal of skeleton extraction is to retain the local peak nodes representing the center line of the trench in the flow distribution matrix, while removing redundant nodes that diverge at the edges, so that the remaining effective nodes form a continuous path with no branches, no width, and single pixel level in the two-dimensional grid.

[0085] For example, a non-maximum suppression method based on neighborhood comparison along the normal direction can be used. First, the geometric normal direction of the line connecting adjacent rigid obstacles corresponding to the current flow distribution matrix is ​​determined. Then, each grid node is traversed, and a comparison neighborhood window of a preset size (e.g., 3×3 or 5×5) is established with the current node as the center. Nodes within the window that intersect with the spatial projection of the normal direction are extracted as comparison nodes. If the flow value of the current node is less than the flow value of any comparison node, or equal to the flow value of a comparison node on one side of a preset priority direction, the flow value of the current node is reset to zero. After traversal, all nodes with non-zero flow values ​​are local peak nodes. These nodes are then connected end-to-end according to an eight-neighbor relationship to obtain a grid sequence with a width of a single grid and nodes sequentially adjacent.

[0086] S502. Reverse map the grid indices in the grid sequence to a set of two-dimensional coordinates of the seabed operating area.

[0087] In one possible implementation, the coordinates of the reference origin and the grid edge length are first obtained when constructing the two-dimensional grid. The reference origin coincides horizontally with the reference origin of the three-dimensional absolute coordinate system established when the underwater robot enters the water. The grid edge length is equal to the single-beam footprint diameter of the multibeam bathysound sonar at the current operating water depth. For each grid cell in the grid sequence, its horizontal and vertical indices represent the cell's number in the matrix, corresponding to the position of the lower left vertex of the cell. To obtain the geographic coordinates of the geometric center of the grid cell, 0.5 (half a grid offset) is added to the index value, multiplied by the grid edge length, and then the reference origin coordinates are added. This process is repeated for each grid cell in the grid sequence, recording the geometric center coordinates of each cell sequentially to form a two-dimensional coordinate point set. The points in this set are arranged in spatial order, reflecting the actual geographic location of the trench centerline.

[0088] S503. Perform linear fitting on the two-dimensional coordinate set to generate the navigation heading.

[0089] In one possible implementation, the points in the two-dimensional coordinate set are arranged in spatial order according to the trench centerline, but these points exhibit a jagged distribution due to grid discretization. To obtain smooth and stable heading commands, a least-squares method is used to fit a straight line to these points, calculating an optimal straight line that minimizes the sum of the squares of the vertical distances from all points to this line. Then, the angle between this fitted straight line and a reference direction (e.g., geographic north or the initial heading of the underwater robot) is calculated; this angle is the navigation heading used to guide the underwater robot.

[0090] The technical solution provided in the above embodiments can bring at least the following beneficial effects: This embodiment performs skeleton extraction processing on the flow distribution matrix to obtain a grid sequence with a width of a single grid and nodes sequentially adjacent (each grid has a grid index). Then, the grid index is inversely mapped to a two-dimensional coordinate set of the seabed operation area. Finally, linear fitting is performed on this coordinate set to generate a navigation heading. The skeleton extraction compresses the planar flow distribution into a single-line feature, eliminating lateral width redundancy and branch interference, ensuring that the extracted sequence is strictly located on the center line of the trench. The inverse mapping restores the abstract graph theory index to real geographic coordinates, providing the robot with a directly understandable metric position. The linear fitting further eliminates the inflection points of the discrete grid sequence, generating a smooth and stable straight heading, enabling the underwater robot to fly stably along the trench direction, providing an accurate navigation reference for subsequent vertical supplementary measurements, and realizing a seamless conversion from virtual simulation results to physical motion control commands.

[0091] In one possible implementation, the process of performing skeleton extraction on the flow distribution matrix to obtain a grid sequence with a width of a single grid and nodes sequentially adjacent can be specifically implemented through the following S601-S604, which will be explained in detail below.

[0092] S601. Determine the geometric normal direction of the line connecting the first rigid obstacle and the second rigid obstacle.

[0093] In one possible implementation, the horizontal center coordinates of the first rigid obstacle and the second rigid obstacle are first obtained. A vector pointing from the first obstacle to the second obstacle is then calculated; the direction of this vector is the direction of the connecting line. In a two-dimensional plane, a straight line has two normal directions, pointing to opposite sides of the line. This step selects one of these normal directions as the reference direction for subsequent neighborhood comparisons. This direction is used to search for nodes to be compared along a direction perpendicular to the trench direction during skeleton extraction.

[0094] S602. Traverse each node in the flow distribution matrix, take the currently traversed node as the current node, establish a comparison neighborhood window of a preset size for the current node, and determine the nodes in the comparison neighborhood window that intersect with the spatial projection of the geometric normal direction as the nodes to be compared.

[0095] In one possible implementation, a square neighborhood window of a preset size is defined centered on the currently traversed node, with the window side length being an odd number (such as 3 or 5). Then, a straight line is drawn along the geometric normal direction and passes through the current node. All grid nodes that this line passes through within the window are extracted and used as the comparison nodes for the current node.

[0096] Understandably, since the current node is located at the center of the window, it is necessary to ensure that the window is symmetrical about the current node in both horizontal and vertical directions so that it can extend evenly to both sides along the normal direction. Odd-numbered side lengths (such as 3, 5, 7) just satisfy this symmetry. If the window is too small (such as 1), it cannot cover adjacent nodes and loses the meaning of comparison. If the window is too large (such as 7 or more), it will introduce too many nodes far away from the current node. These nodes may exceed the effective range of the trench cross section, leading to misjudgment of retaining non-local peaks and increasing the computational burden. A 3×3 window can cover one node each above, below, left, and right in the normal direction, which is suitable for scenarios with narrow trench width and concentrated flow distribution. A 5×5 window can cover a wider range and is suitable for scenarios with larger trench width or more dispersed flow distribution. In engineering practice, it can be pre-selected according to the historical scour width of the operating sea area or the pixel resolution of the flow distribution matrix.

[0097] S603. Compare the current node's traffic value with the traffic values ​​of all nodes to be compared. If the current node's traffic value is less than the traffic value of any node to be compared, then clear the current node's traffic value to zero.

[0098] In one possible implementation, the flow value of the current node is retrieved and compared one by one with the flow values ​​of all nodes to be compared. If the flow value of the current node is less than the flow value of any node to be compared, it means that the current node is not the flow peak on the normal cross section (i.e., it is located at the edge of the trench or on a downhill slope), and the flow value of the current node is set to zero; if the flow value of the current node is greater than or equal to the flow values ​​of all nodes to be compared, its original flow value is retained.

[0099] S604. After traversing all nodes, connect each node with a non-zero flow value in series according to the adjacency relationship to generate a grid sequence with a width of a single grid and nodes sequentially adjacent to each other.

[0100] In one possible implementation, after all nodes in the flow distribution matrix have completed the aforementioned comparison and zeroing operations, only nodes with non-zero flow values ​​remain in the matrix. These nodes are the local peak nodes along the normal direction. These nodes are then connected end-to-end according to spatial adjacency relationships (e.g., eight-neighborhood adjacency) to form a continuous sequence. In this sequence, adjacent nodes are adjacent to each other in the grid, and the width of the sequence is strictly one grid cell, with no branches or discontinuities.

[0101] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment specifically defines the steps of skeleton extraction processing: determining the geometric normal direction of the line connecting adjacent rigid obstacles, traversing each node in the flow distribution matrix, establishing a comparison neighborhood window of a preset size for the current node, extracting nodes within the window that intersect with the spatial projection of the normal direction as nodes to be compared, and if the flow value of the current node is less than any node to be compared, then its flow value is cleared to zero. Finally, all non-zero flow nodes are connected in series according to their adjacency relationship to generate a single-line continuous grid sequence. By comparing neighborhoods based on the normal direction, local peak nodes at the centerline of the trench can be accurately preserved, while secondary peaks and noise flow values ​​diverging at the edges are cleared to zero. This processing method not only ensures the continuity of the centerline but also eliminates the interference of width and branches, making the width of the generated grid sequence strictly a single grid, providing high-quality single-line skeleton data for subsequent coordinate restoration and heading fitting, and significantly improving the accuracy and robustness of trench centerline extraction.

[0102] In one possible implementation, the process of reverse mapping the grid indexes in the grid sequence to a two-dimensional coordinate set of the seabed operating area can be specifically implemented through the following S701-S703, which will be explained in detail below.

[0103] S701. Obtain the coordinates of the origin and the grid edge length when constructing the connected graph.

[0104] In one possible implementation, two pre-set and saved parameters are read from the underwater robot's internal storage unit: the reference origin coordinates and the grid edge length. The reference origin coordinates are the origin of the absolute coordinate system established on the horizontal plane when the underwater robot enters the water. This origin coincides with the horizontal projection of the entry point in the three-dimensional spatial coordinate system and has been calibrated and stored before the operation begins. The grid edge length is equal to the single-beam footprint diameter of the multibeam echo sounder at the current operating water depth. This value is pre-determined and stored in the system before the operation based on the sonar equipment's technical parameters and water depth conditions.

[0105] S702. For each grid cell in the grid sequence, determine the two-dimensional coordinates of the geometric center of each grid cell based on its horizontal and vertical indices, the coordinates of the reference origin, and the grid side length.

[0106] In one possible implementation, since the grid index is obtained by rounding down and projecting, the index value corresponds to the position of the lower left corner vertex of the grid cell (or the corner point closest to the origin), rather than the center of the grid cell. To obtain the geographic coordinates of the geometric center of the grid cell, 0.5 (i.e., an offset of half the grid side length) needs to be added to the index value to shift the index from the corner point to the center; then, this result is multiplied by the grid side length to obtain the metric offset distance from the reference origin to the grid center point; finally, the coordinates of the reference origin are added to obtain the absolute coordinates of the geometric center of the grid cell. For the horizontal coordinate, the calculation formula is: the horizontal coordinate of the reference origin plus (horizontal index plus 0.5) multiplied by the grid side length; the same processing is performed for the vertical coordinate, and the calculation is performed sequentially for each grid cell in the grid sequence to obtain the two-dimensional coordinates of the geometric center of each cell.

[0107] S703. The least squares method is used to linearly fit the two-dimensional coordinates of the geometric center of each grid cell, and the navigation heading is determined based on the angle between the fitted line and the reference direction.

[0108] For example, suppose the two-dimensional coordinate point set obtained through inverse mapping contains 50 points, roughly distributed in a northeast-southwest direction. After fitting using the least squares method, a straight line is obtained, which has an angle of 135 degrees relative to true north. This angle is used as the navigation heading output to the flight control module of the underwater robot, allowing the robot to adjust its yaw angle to 135 degrees and fly stably along the centerline of the trench. Even if the point set is more discrete, the least squares method can still provide a statistically optimal straight line that reflects the overall trend of the trench and avoids frequent heading fluctuations caused by local curvature or noise.

[0109] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment specifically defines the steps of generating a two-dimensional coordinate set through reverse mapping: obtaining the reference origin coordinates and grid edge length when constructing the connected graph; for each grid cell in the grid sequence, determining the two-dimensional coordinates of the geometric center of the cell based on its horizontal index, vertical index, reference origin coordinates, and grid edge length; then using the least squares method to linearly fit all geometric center coordinates; and determining the navigation heading based on the angle between the fitted line and the reference direction. By calculating the geometric center (implying the offset compensation of half a grid), the systematic positional deviation caused by the downward rounding operation during forward projection is accurately eliminated, and the graph theory index is accurately restored to the real metric coordinates; the least squares method fitting of the line further smooths the discrete point set, eliminates errors, and obtains the optimal straight heading angle, ensuring that the underwater robot yaw controller can stably and oscillate to execute navigation commands, significantly improving the trajectory tracking accuracy and engineering practicality of local rescanning.

[0110] In one possible implementation, the process of obtaining the discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, as well as the position information of rigid obstacles in the seabed topography, can be specifically implemented through the following S801-S803, which will be described in detail below.

[0111] S801 receives the echo signal from the multibeam sonar and generates a basic three-dimensional point cloud set.

[0112] In one possible implementation, a multibeam echo sounder mounted on an underwater robot transmits sound waves to the seabed and receives echo signals from the seabed sediment and underwater structures. The sonar receiver converts the echo signals into digital sampled data, calculates the three-dimensional spatial coordinates of each echo based on the time of arrival and angle of arrival, and records the energy integral of the echo as the discrete acoustic reflection intensity, thereby generating a basic three-dimensional point cloud set. Each point in this point cloud set contains three-dimensional spatial coordinates (x-coordinate, y-coordinate, and elevation) and the corresponding discrete acoustic reflection intensity.

[0113] S802. Extract the 3D point cloud of the rigid obstacle from the basic 3D point cloud set using a geometric fitting algorithm, and determine the horizontal center coordinates and fitting radius of the rigid obstacle.

[0114] In one possible implementation, since the subsea pipe pile has the geometric characteristics of a vertical cylinder, the Random Sampling Consensus (RANSAC) algorithm can be used to set the cylindrical surface model. Several points are randomly sampled from the point cloud set, and the cylinder parameters (axial direction, central axis position, radius) are fitted. Then, the points in the point cloud that conform to the model are statistically analyzed, and these points are grouped into a solid pile foundation point cloud. For each identified independent cylinder, the intersection of its central axis and the nautical chart horizontal reference surface (i.e., the water depth datum) is calculated to obtain the horizontal center coordinates. Simultaneously, the radius parameter is extracted from the fitted cylindrical model as the fitted radius for that pile foundation.

[0115] S803 receives the high-frequency encrypted point cloud obtained by the underwater robot based on control commands, and determines the horizontal center coordinates and fitting radius of the rigid obstacle as the position information of the rigid obstacle.

[0116] In one possible implementation, the horizontal center coordinates and fitting radius of each rigid obstacle obtained from the geometric fitting algorithm in the previous step are associated and stored as the location information of that rigid obstacle. Specifically, for each identified pile foundation, its horizontal center coordinates and fitting radius are bound as a data record and stored in the underwater robot's memory or storage unit.

[0117] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment defines the specific method for obtaining discrete acoustic reflection intensity and rigid obstacle location information: receiving the echo signal from a multibeam bathymetry sonar to generate a basic three-dimensional point cloud set (each point contains spatial coordinates and discrete acoustic reflection intensity), extracting the three-dimensional point cloud of the rigid obstacle from the point cloud through a geometric fitting algorithm, and determining its horizontal center coordinates and fitting radius as the location information of the rigid obstacle. By directly obtaining the required data from the original sonar echo, seamless connection between data acquisition and subsequent processing is achieved, without the need for additional sensors or manual annotation, reducing system complexity and deployment costs; the geometric fitting algorithm can accurately separate columnar obstacles from the mixed point cloud and calculate the precise center and radius, providing reliable parameter basis for subsequent delineation of injection and outflow boundary areas, ensuring the adaptability and robustness of the entire method to complex seabed environments.

[0118] In one possible implementation, it is also necessary to generate a three-dimensional seabed topographic model. This process can be specifically implemented through the following S901-S902, which will be explained in detail below.

[0119] S901. Generate control commands based on the navigation heading and the position coordinates of the first rigid obstacle and the second rigid obstacle.

[0120] Among them, the control commands are used to adjust the underwater robot's course to the navigation course, and to control the underwater robot to fly at a preset altitude along the navigation course, perform vertical detection, and acquire high-frequency encrypted point clouds.

[0121] In one possible implementation, a complete control command is generated using the determined navigation heading and the position coordinates of the adjacent rigid obstacle pair (the first rigid obstacle and the second rigid obstacle) as input. This command guides the underwater robot to adjust its heading to the navigation heading and fly along that heading to the area above the trench between the first and second rigid obstacles. Specifically, the control command contains at least two fields: target heading angle and target position range. The target heading angle is the navigation heading (angle value relative to the reference direction) determined in the previous steps, and the target position range is calculated based on the horizontal center coordinates of the first and second rigid obstacles. For example, the midpoint of the line connecting the two coordinates can be taken as the starting position for the robot to perform vertical detection.

[0122] S902. The high-frequency encrypted point cloud and the basic three-dimensional point cloud are superimposed and fused in the same coordinate system to generate a three-dimensional seabed landform model.

[0123] In one possible implementation, a high-frequency encrypted point cloud acquired through local vertical detection and a basic 3D point cloud set generated by an initial sparse scan are used as inputs. These two sets are then superimposed and fused within the same absolute coordinate system to generate a complete 3D seabed topographic model. Specifically, since both scans are performed by the same underwater robot in a single underwater operation and are based on the same absolute coordinate system origin established at the time of entry, the two point clouds naturally share the same spatial reference system, eliminating the need for complex point cloud registration calculations. During fusion, each point in the high-frequency encrypted point cloud is directly appended to the basic 3D point cloud set. For points with overlapping spatial locations (i.e., both sparse and encrypted points exist within the same grid cell), the elevation data of the encrypted points is preferentially retained, or the depth values ​​of the encrypted points are used to replace the missing or invalid elevation values ​​in the sparse points. After fusion, the merged complete point cloud set is rendered using surface meshing to generate a visualized 3D seabed topographic model.

[0124] The technical solution provided by the above embodiments can bring at least the following beneficial effects: After determining the navigation heading, this embodiment generates control commands based on the navigation heading and the position coordinates of the first and second rigid obstacles, adjusts the underwater robot's heading to the navigation heading, and controls the robot to fly along the heading at a preset height to perform vertical detection to obtain high-frequency encrypted point clouds. Finally, the high-frequency encrypted point clouds and the basic three-dimensional point clouds are superimposed and fused in the same coordinate system to generate a complete three-dimensional seabed topographic model. By using the precise heading obtained through deduction to guide the robot to perform local vertical detection only directly above the trench, the vertical incident attitude fundamentally avoids side wall occlusion and multipath effects, allowing the wave arrival time calculation algorithm to return to normal and successfully obtaining high-precision elevation data of the bottom of the deep trench. At the same time, only local high-frequency detection is performed on the hidden danger area instead of global high-frequency scanning, enabling precise positioning under the condition of limited battery power. Finally, through the fusion of the same source coordinate reference, both the large-scale outline of the global sparse scan is preserved, and the fine depth information of key parts is supplemented, forming a complete three-dimensional seabed topographic model, providing reliable data support for the stability assessment of pipe piles.

[0125] Please see Figure 2 This document illustrates a schematic diagram of the system architecture of a seabed topography scanning system 200 based on an underwater robot, provided in an embodiment of the present invention. This system implements the aforementioned seabed topography scanning method based on an underwater robot. The system includes: an acquisition module 201, used to acquire discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, and the location information of rigid obstacles in the seabed topography; the discrete acoustic reflection intensity is used to characterize the looseness of the seabed sediment; a connected graph construction module 202, used to construct a connected graph based on the discrete acoustic reflection intensity; nodes in the connected graph have capacity values; the capacity values ​​are determined based on the acoustic reflection intensity, with nodes having higher capacity values ​​for lower acoustic reflection intensity; a boundary determination module 203, used to determine the injection boundary region and the outflow boundary region based on the location information of the rigid obstacles; a deduction module 204, used to execute a maximum flow algorithm constrained by the capacity value of each node in the connected graph to determine the flow distribution matrix; the flow distribution matrix is ​​used to characterize the degree of flow aggregation; and a navigation determination module 205, used to determine the navigation course for guiding the underwater robot based on the flow distribution matrix.

[0126] The technical solution provided by the above embodiments can bring at least the following beneficial effects: This embodiment realizes the method steps in the above embodiments through the collaborative work of the acquisition module, connected graph construction module, boundary determination module, inference module, and navigation determination module. By decoupling the functions of data acquisition, graph theory graph construction, boundary delineation, maximum flow inference, and navigation generation through modular design, each module has a clear division of labor and a clear data flow. Under the premise of maintaining the low energy consumption and large spacing sparse scanning of the underwater robot, it can accurately infer the trench direction using acoustic reflection intensity data and guide local vertical supplementary measurement. Overall, it solves the technical contradictions of traditional methods that cannot locate hidden dangers due to elevation failure or exhaust power due to global high-frequency detection, significantly reducing operation and maintenance costs and improving the intelligence level and engineering practicality of marine engineering structure safety monitoring.

[0127] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0128] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A seabed topographic scanning method based on an underwater robot, characterized in that, The method includes: The discrete acoustic reflection intensity is obtained from the seabed topography determined by the underwater robot scanning, as well as the location information of rigid obstacles in the seabed topography; the discrete acoustic reflection intensity is used to characterize the looseness of the seabed substrate. A connected graph is constructed based on the discrete acoustic reflection intensities; the nodes in the connected graph have capacity values; the capacity values ​​are determined based on the acoustic reflection intensities, with nodes having higher capacity values ​​for nodes with lower acoustic reflection intensities. Based on the location information of the rigid obstacle, the injection boundary region and the outflow boundary region are determined; Using the capacity value of each node in the connected graph as a constraint, the maximum flow algorithm is executed to determine the flow distribution matrix; the flow distribution matrix is ​​used to characterize the degree of flow aggregation. Based on the flow distribution matrix, determine the navigation course used to guide the underwater robot; The step of constructing a connected graph based on the discrete acoustic reflection intensities includes: The discrete acoustic reflection intensity is projected onto a two-dimensional grid to generate a grid-averaged intensity matrix; each grid cell in the grid-averaged intensity matrix corresponds to a node in the connected graph. Based on the intensity extreme values ​​in the average intensity matrix, the acoustic reflection intensity of each grid cell in the average intensity matrix is ​​normalized to determine the normalization difference of each grid cell. According to a preset nonlinear mapping relationship, the normalized difference of each grid cell is converted into the capacity value of each node in the connected graph; the nonlinear mapping relationship makes the node with the smaller normalized difference obtain a higher capacity value. Each node with a capacity value is connected according to its adjacency relationship, and the capacity value of the connected node is assigned to the connecting edge to generate the connected graph.

2. The seabed topography scanning method based on an underwater robot according to claim 1, characterized in that, The step of determining the injection boundary region and the outflow boundary region based on the position information of the rigid obstacle includes: Identify adjacent first and second rigid obstacles; the first and second rigid obstacles are pairs of rigid obstacles with a horizontal distance of less than a preset threshold within the seabed operating area; Based on the position of the first rigid obstacle, a first surrounding area is defined as the injection boundary area according to the fitting radius of the first rigid obstacle and the preset search range. Based on the position of the second rigid obstacle, and according to the fitting radius of the second rigid obstacle and the preset search range, a second surrounding area is defined as the outflow boundary area.

3. The seabed topography scanning method based on an underwater robot according to claim 2, characterized in that, The step of constraining the capacity value of each node in the connected graph and executing the maximum flow algorithm to determine the flow distribution matrix includes: Establish the source and sink nodes for the maximum flow algorithm; Connect the source point to each node within the injection boundary region, and connect each node within the outflow boundary region to the sink point; Using the source and sink nodes as inputs and the capacity of each node as a constraint, a preset maximum flow algorithm is run. The algorithm is monitored for the number of iterations and the number of active nodes in the network. When the number of iterations reaches the preset limit or there are no active nodes in the network, the algorithm is determined to have converged and the push is stopped. The flow distribution matrix is ​​generated based on the cumulative flow value of each node under convergence.

4. The seabed topography scanning method based on an underwater robot according to claim 2, characterized in that, Determining the navigation course for guiding the underwater robot based on the flow distribution matrix includes: The flow distribution matrix is ​​subjected to skeleton extraction processing to obtain a grid sequence with a width of a single grid and nodes that are sequentially adjacent; each grid in the grid sequence has a grid index; The grid indices in the grid sequence are inversely mapped to a set of two-dimensional coordinates of the seabed operating area; A linear fit is performed on the set of two-dimensional coordinates to generate the navigation heading.

5. The seabed topography scanning method based on an underwater robot according to claim 4, characterized in that, The process of performing skeleton extraction on the flow distribution matrix to obtain a grid sequence with a width of a single grid and nodes sequentially adjacent to each other includes: Determine the geometric normal direction of the line connecting the first rigid obstacle and the second rigid obstacle; Traverse each node in the flow distribution matrix, take the currently traversed node as the current node, establish a comparison neighborhood window of a preset size for the current node, and determine the nodes in the comparison neighborhood window that intersect with the spatial projection of the geometric normal direction as nodes to be compared. The current node's traffic value is compared with the traffic values ​​of all nodes to be compared. If the current node's traffic value is less than the traffic value of any node to be compared, the current node's traffic value is cleared to zero. After traversing all nodes, each node with a non-zero flow value is connected in series according to its adjacency relationship to generate a grid sequence with a width of a single grid and nodes sequentially adjacent to each other.

6. The seabed topography scanning method based on an underwater robot according to claim 4, characterized in that, The step of inversely mapping the grid indices in the grid sequence to a set of two-dimensional coordinates of the seabed operating area includes: Obtain the coordinates of the origin and the grid edge length when constructing a connected graph; For each grid cell in the grid sequence, the two-dimensional coordinates of the geometric center of each grid cell are determined based on its horizontal and vertical indices, the coordinates of the reference origin, and the grid side length. The least squares method is used to linearly fit the two-dimensional coordinates of the geometric center of each grid cell, and the navigation heading is determined based on the angle between the fitted line and the reference direction.

7. The seabed topography scanning method based on an underwater robot according to claim 6, characterized in that, The acquisition of discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, and the location information of rigid obstacles in the seabed topography, includes: The echo signal from the multibeam sonar is received to generate a basic three-dimensional point cloud set; each point in the basic three-dimensional point cloud set contains three-dimensional spatial coordinates and discrete acoustic reflection intensity. From the basic 3D point cloud set, the 3D point cloud of the rigid obstacle is extracted by a geometric fitting algorithm, and the horizontal center coordinates and fitting radius of the rigid obstacle are determined. The horizontal center coordinates and fitting radius of the rigid obstacle are used to determine the position information of the rigid obstacle.

8. The seabed topography scanning method based on an underwater robot according to claim 7, characterized in that, The method further includes: Based on the navigation heading and the position coordinates of the first rigid obstacle and the second rigid obstacle, control commands are generated; the control commands are used to adjust the underwater robot's heading to the navigation heading, and to control the underwater robot to fly along the navigation heading at a preset altitude, perform vertical detection, and acquire high-frequency encrypted point clouds; The system receives a high-frequency encrypted point cloud obtained by an underwater robot based on control commands, and then overlays and fuses the high-frequency encrypted point cloud with the basic three-dimensional point cloud set in the same coordinate system to generate a three-dimensional seabed topographic model.

9. A seabed topography scanning system based on an underwater robot, characterized in that, The system is used to implement the seabed topographic scanning method based on an underwater robot as described in any one of claims 1-8, and the system includes: The acquisition module is used to acquire discrete acoustic reflection intensity determined by the underwater robot scanning the seabed topography, as well as the position information of rigid obstacles in the seabed topography; the discrete acoustic reflection intensity is used to characterize the looseness of the seabed substrate. A connected graph construction module is used to construct a connected graph based on the discrete acoustic reflection intensity; the nodes in the connected graph have capacity values; the capacity values ​​are determined based on the acoustic reflection intensity, and the lower the acoustic reflection intensity, the higher the capacity value of the node; The boundary determination module is used to determine the injection boundary region and the outflow boundary region based on the position information of the rigid obstacle; The deduction module is used to execute the maximum flow algorithm with the capacity value of each node in the connected graph as a constraint to determine the flow distribution matrix; the flow distribution matrix is ​​used to characterize the degree of flow aggregation. The navigation determination module is used to determine the navigation course for guiding the underwater robot based on the flow distribution matrix.

Images