A multi-beam data outlier detection method based on adaptive progressive morphological filter

Abstract

This invention discloses a method for outlier detection in multi-beam data based on adaptive progressive morphological filtering. The method employs a closed-open-open combined average morphological operator for filtering, extracting a local terrain slope tolerance parameter. Based on this parameter, the elevation threshold is adaptively and dynamically adjusted during iterative morphological filtering. Through multi-level filtering windows and iterative calculations, the morphological reference plane gradually approximates the actual seabed until the iteration termination condition is met. Then, the terrain slope is extracted and a global judgment threshold is calculated. The depth deviation between the actual measured depth and the reference elevation is compared based on spatial location, thereby achieving accurate removal of outliers from the original discrete point cloud. By establishing a local dynamic adjustment and global decision mechanism for the filtering threshold, the method effectively avoids the filtering and under-filtering defects caused by fixed parameters, improving denoising accuracy while ensuring high fidelity of the actual seabed micro-topography.

Description

Technical Field

[0001] This invention relates to a method for quality control and post-processing of multibeam bathymetry data, and more particularly to a method for detecting outliers in multibeam bathymetry data based on adaptive progressive morphological filtering. Background Technology

[0002] Multi-beam echo sounders (MBES) possess full-coverage, high-precision, and high-resolution detection capabilities, making them a mainstream technology for modern marine surveying, seabed topography investigation, and underwater engineering exploration. With the development of acoustic technology, the amount of point cloud data acquired by multi-beam equipment has increased rapidly, enabling detailed characterization of seabed micro-topography features. However, in actual data acquisition, due to the extremely complex underwater environment (such as ambient noise, complex sea conditions, fish schools, wakes, and bubble layers), as well as the systematic and random errors of the measuring equipment itself, the raw multi-beam echo sounder data inevitably contains a large number of outliers. The presence of these outliers seriously affects the accuracy of subsequent digital elevation model (DEM) construction. Therefore, the detection and removal of outliers is crucial in the post-processing stage of multi-beam data.

[0003] Currently, methods for handling outliers in multibeam bathymetry data mainly include manual interactive editing, traditional statistical algorithms, spatial clustering-based algorithms, and classical mathematical morphological filtering. However, these methods all have certain limitations in practical applications. While manual interactive editing offers high accuracy and preserves realistic terrain well, it is time-consuming, labor-intensive, and inefficient when dealing with massive amounts of multibeam data. Furthermore, it heavily relies on the experience of operators, leading to significant subjectivity in the results. Traditional statistical algorithms, such as trend surface filtering and robust estimation, are typically based on the assumption of local flatness and perform well in areas with gentle terrain. However, in areas with complex underwater topography, these algorithms often fail to accurately distinguish between terrain undulations and noise, easily misclassifying real terrain features such as protruding reefs as outliers and filtering them out, resulting in a distorted final terrain model. In addition, spatial clustering-based outlier detection algorithms are also applied to point cloud denoising. However, multibeam bathymetry data is massive and its spatial distribution density is extremely uneven. Clustering algorithms not only face bottlenecks such as high computational complexity and low processing efficiency, making it difficult to meet the needs of rapid processing of massive data; at the same time, since complex underwater noise often forms local high-density point cloud clusters, clustering algorithms are very likely to misjudge them as real terrain features, and globally fixed distance or density thresholds are difficult to adapt to drastic spatial changes in point cloud density.

[0004] Therefore, to better handle complex terrain, Progressive Morphological Filter (PMF) has been gradually introduced into multibeam bathymetry data processing in recent years. This method separates terrain from noise by progressively increasing the filter window size and elevation threshold. However, the window size and elevation threshold of traditional PMF algorithms are mostly globally fixed parameters or simple linear growth patterns set empirically. When facing seabed terrain with dramatic undulations and diverse slopes, fixed parameter models lack the ability to adaptively perceive local terrain features. Specifically, setting the threshold too small will lead to over-smoothing of terrain features, producing a "filtered wave" phenomenon; setting the threshold too large will result in a large amount of residual noise in the data, leading to an "under-filtering" problem. In summary, to overcome the contradiction between "preserving terrain" and "removing noise" in the processing of complex underwater terrain, there is an urgent need in this field for a progressive morphological outlier detection method that can adaptively and dynamically adjust the filter parameters according to local terrain undulation features, so as to achieve efficient and high-precision automated cleaning of multibeam bathymetry data. Summary of the Invention

[0005] This invention proposes an adaptive progressive morphological filtering method for outlier detection in multibeam data. This method effectively solves the problems of false detection of negative outliers in multibeam data by existing progressive morphological filtering methods and the inability of fixed algorithm parameters to adapt to complex and varied terrain changes, thereby improving the adaptability and accuracy of multibeam outlier detection.

[0006] To achieve the above objectives, the present invention adopts the following technical solution:

[0007] An adaptive progressive morphological filtering method for outlier detection in multi-beam data, comprising the following steps:

[0008] Step 1: Rasterize the multibeam point cloud data. Based on the spatial location information of the multibeam bathymetry data, divide it into a regular two-dimensional data grid, and use the minimum elevation of the bathymetry points contained in each grid as the elevation value of the current grid. Holes may appear during the grid division process. The nearest neighbor interpolation method is used to interpolate and fill the elevation values ​​of the holed grids, forming the minimum surface grid of the multibeam point cloud, denoted as . .

[0009] Step 2: Determine the morphological filtering operator. To overcome the problem that commonly used opening operators cannot effectively identify multi-beam negative outliers, a closed-open-open combined averaging operator is adopted. Morphological filtering is applied to the point cloud mesh.

[0010] Step 3: Calculate the terrain slope tolerance parameter. To improve the algorithm's adaptability to different terrain change trends, the terrain slope tolerance parameter is calculated in real time based on the multibeam point cloud data to be processed. In the filtering window for Below, the closed-open and open-closed combined average operator is adopted ( ) point cloud mesh Morphological filtering is performed to obtain a smoother terrain mesh surface that reflects terrain transformation trends; then, the grid is traversed to calculate the terrain slope tolerance parameter for local areas. .

[0011] Step 4: Perform progressive morphological filtering. A single filtering window cannot accommodate outlier information at different scales. Therefore, a progressively increasing filtering window is used to perform multi-level iterative operations on the surface grid, with the filtering window growing linearly. In the... In the next iteration, with Morphological structural elements are constructed for the dimensions, and morphological operations are performed on the elevation surfaces retained from the previous iteration to filter out convex elevation points at the current window scale, thereby obtaining the 1st... The reference smoothing datum plane for the next iteration.

[0012] Step 5: Marking Suspected Outliers. First, based on the terrain slope tolerance parameters pre-calculated in Step 3... Real-time calculation of corresponding grid nodes Local adaptive elevation difference threshold Its calculation model is as follows:

[0013]

[0014] In the formula, For the first The local adaptive elevation difference threshold is obtained from the next iteration. The initial height difference threshold is the height difference threshold from the previous iteration. This is typically set as the nominal depth sounding error limit for a multibeam echo sounding system. and These are the filter window sizes for the current round and the previous round, respectively. The horizontal physical resolution of the surface grid.

[0015] Compare the first Next and first The elevation difference at the corresponding horizontal position on the secondary reference smoothing datum surface, if the absolute value of the elevation difference at a certain node in the grid is greater than the local adaptive elevation difference threshold. If the absolute value of the elevation difference is less than or equal to the threshold, it is determined to be a real terrain point at that scale, and its original elevation value is retained to participate in the next iteration.

[0016] Step 6: Determine whether to terminate the iterative filtering. Considering that an excessively large filtering window might smooth out real terrain features, and to improve the algorithm's computational efficiency, an iteration termination condition is designed based on the adjacent elevation differences in the iteration results to avoid over-filtering due to unreasonable filtering window settings. Therefore, the iteration termination condition is set as follows:

[0017]

[0018] in, This is the depth error threshold. During the iteration process, if the absolute value of the elevation difference between adjacent filtering results is less than the depth error threshold... If the filtering process is considered to have stabilized, the iteration process ends. Simultaneously, if the structure window reaches its set maximum value, the iterative filtering also terminates, and the output consists of all marked outlier points from the iteration process. Otherwise, proceed to step 4.

[0019] Step 7: Generate a terrain reference surface. Considering that the progressive morphological filtering process essentially smooths the terrain mesh, and may lose some detailed terrain surface features while marking suspected outliers, the terrain surface mesh output by the filter is not used directly. Instead, the marked suspected outlier points are placed on the original point cloud mesh. The mesh is then culled, and the resulting mesh holes are interpolated to form a new reference point cloud mesh, denoted as . .

[0020] Step 8: Outlier Detection and Removal. First, calculate the outlier judgment threshold based on the reference point cloud mesh. Calculate terrain slope information The size of the filtering window at the end of the iterative filtering process. This is to support neighborhood calculation. An outlier detection threshold can be calculated. :

[0021]

[0022] in, This is a scaling factor, typically set to 1.2-1.5 based on the acoustic reflection characteristics of the substrate and the complexity of the terrain in the survey area. The horizontal physical resolution of the point cloud mesh; This is the size of the filter window when the iterative filtering terminates.

[0023] Horizontal coordinates based on actual depth sounding points In point cloud mesh Interpolation values ​​provide in-depth reference information. By comparing with actual measured depth information To identify outliers, the identification criteria are as follows:

[0024]

[0025] in, For actual depth measurement points Actual measured depth information at the location; For depth sounding points In the reference terrain grid Interpolation values ​​provide reference depth information; For depth sounding points The outlier threshold is calculated at the specified location.

[0026] This enables the detection and removal of outliers in multi-beam data.

[0027] As a further improvement of this invention, in order to overcome the limitations of morphological opening operations that easily lead to deformation of the reference terrain and unidirectional elevation bias, this method uses a closed-open-open-closed combined average morphological operator to pre-smooth the original point cloud elevation grid before extracting terrain features and in the subsequent progressive morphological iterative filtering process. This effectively offsets the systematic deviation caused by unidirectional operations and achieves the preservation of the reference terrain contour.

[0028] As a further improvement of the present invention, in order to solve the problem that traditional morphological filtering is prone to generating filtered waves and under-filtering in complex underwater terrain, step 4 breaks the limitation of the global fixed threshold, and substitutes the terrain slope tolerance parameter as an adaptive weight coefficient into the progressive iterative formula. The elevation difference threshold is dynamically adjusted according to the local terrain undulation characteristics, so as to remove high-frequency clutter while protecting the real terrain protrusion.

[0029] As a further improvement of the present invention, in order to ensure the spatial distribution characteristics and original measurement accuracy of the original multibeam discrete sounding data, in step 8, the macroscopic terrain slope is extracted based on the reference surface after all noise is removed after the iteration is terminated to establish a global decision standard. The reference depth and the decision threshold are mapped back to the original discrete sounding points for point-by-point verification through spatial interpolation technology, thereby achieving accurate removal of discrete point cloud outliers.

[0030] Advantages compared to existing technologies:

[0031] This invention overcomes the limitations of traditional morphological filtering, which relies on fixed parameters and whose traditional opening operators easily lead to asymmetric shifts in terrain contours and systematic elevation reductions. By introducing a closed-opening-opening-closeding combined averaging morphological operator, it effectively offsets the systematic elevation deviation caused by unidirectional calculations, achieving faithful contour reproduction of the baseline terrain. Furthermore, by combining terrain slope tolerance parameters and terrain slope, a local adaptive dynamic adjustment mechanism for the filtering threshold and an outlier detection standard are constructed. This method fundamentally overcomes the false positives and false negatives that easily occur in complex underwater terrain. It not only improves the computational efficiency and automation level of outlier cleaning for massive bathymetry data but also achieves an effective balance between accurate removal of multi-scale spatial noise and high fidelity of real, minute terrain features, effectively improving the construction accuracy of seabed digital elevation models.

[0032] To address the issues of false positives and false negatives in traditional progressive morphological filtering methods due to fixed parameters when dealing with drastic terrain undulations, this paper proposes an adaptive progressive morphological outlier detection method suitable for complex underwater terrain by improving the setting of morphological filtering operators and filtering thresholds based on multibeam local terrain features. This method enhances the denoising accuracy and terrain fidelity of the method.

[0033] In summary, the improved adaptive PMF method resolves the sharp contradiction between preserving topographic features and eliminating high-frequency clutter in complex seabed environments. Attached Figure Description

[0034] Figure 1 The flowchart shows a multi-beam data outlier detection method based on adaptive progressive morphological filtering.

[0035] Figure 2 1. Actual measured multibeam bathymetry data;

[0036] Figure 3 1 is the multibeam bathymetry data filtered by the method of the present invention;

[0037] Figure 4 2. Actual measured multibeam bathymetry data;

[0038] Figure 5 2. This is the multibeam bathymetry data filtered by the method of the present invention. Detailed Implementation

[0039] Example 1, refer to Appendix Figure 1-5 An adaptive progressive morphological filtering method for outlier detection in multi-beam data is described in the attached algorithm flowchart. Figure 1 As shown, the specific steps to implement the method are as follows:

[0040] 1) Acquisition and preprocessing of multibeam bathymetry data

[0041] Acquire raw multibeam bathymetry point cloud data and extract the position and depth information of each measuring point, i.e., the three-dimensional spatial coordinates. Based on the point cloud density and terrain features of the survey area, an appropriate grid resolution is set, and the discrete, scattered point cloud data is regularized into a two-dimensional elevation grid. To address the data gaps caused by uneven point cloud distribution during the gridding process, the nearest neighbor interpolation method is used to interpolate and fill the elevation values ​​of the gapped grids, forming the minimum surface grid for the multibeam point cloud, denoted as […]. .

[0042] 2) Determine the morphological filtering operator

[0043] Since single opening or closing operations can easily cause asymmetric morphological shifts in terrain contours, this method employs a combined averaging operator to effectively offset the systematic biases introduced by unidirectional morphological operations. This maintains the volume conservation and boundary features of the baseline terrain while filtering out positive and negative high-frequency random noise. The specific mathematical model of this combined averaging operator is expressed as follows:

[0044]

[0045] in, For data to be processed, For structural elements, select a rectangular structural element here; and They are respectively Opening and closing and closed-open form operators, and Let represent the opening and closing operations, respectively, and define them as follows:

[0046]

[0047] in, and These represent the erosion and dilation operations, respectively, and are defined as follows:

[0048]

[0049] 3) Calculate terrain slope tolerance parameters

[0050] To improve the algorithm's adaptive perception of different underwater terrain change trends and avoid interference with the extraction of real terrain features caused by local anomalies in the original multibeam data, this invention calculates the terrain slope tolerance parameter based on macroscopic terrain trends in real time based on the multibeam point cloud grid data to be processed.

[0051] In morphological filtering window size Below, the closed-open and open-closed combined average operator is adopted ( ) on the original point cloud elevation grid surface Morphological pre-filtering is performed by calculating the results of the first-close-then-open operation and the first-open-then-close operation, and then averaging the two to obtain a relatively smooth terrain grid reference surface that effectively suppresses local high-frequency clutter and presents the true macroscopic undulations of the terrain. Based on this smoothed terrain grid reference surface, the terrain slope tolerance parameter is further calculated. The specific process is as follows:

[0052] First, to quantify the severity of local terrain undulations, arbitrary target nodes on the terrain grid surface are smoothed. Construct a local slope calculation window centered on [the target area]. Traverse all neighboring nodes within the window. Calculate the spatial elevation change rate between the target node and each neighboring node, and extract the maximum value as the local terrain slope of the target node. Its calculation model is as follows:

[0053]

[0054] In the formula, and These represent the elevation values ​​of the target node and the neighboring nodes on the smooth grid surface, respectively; the denominator is the Euclidean distance between the two nodes on the horizontal projection plane.

[0055] Secondly, because a few underwater man-made structures or extremely strong residual random noise may cause the local calculated slope to tend to infinity, in order to ensure the stability of the adaptive filtering system, it is necessary to adjust the calculated local terrain slope. By applying truncation constraints, the final terrain slope tolerance parameters are obtained. An empirical threshold for the maximum topographic slope of the survey area is introduced. Its truncation model is as follows:

[0056]

[0057] The terrain slope tolerance parameter obtained through the above steps This method not only filters out high-frequency noise from the original data from interfering with slope calculations but also fully preserves the macroscopic slope characteristics of the actual terrain. In subsequent progressive morphological filtering iterations, this parameter will be used as an adaptive weighting coefficient in the elevation difference threshold calculation formula: in areas with large, volatile parameter values, the algorithm will automatically relax the elevation difference threshold to protect the actual terrain; in flat seabed areas with smaller parameter values, the algorithm will automatically tighten the elevation difference threshold to ensure that minor clutter is eliminated. This makes the calculated terrain slope tolerance parameter more robust and reliable, providing a reference for the adaptive dynamic adjustment of the elevation difference threshold in subsequent progressive filtering processes.

[0058] 4) Progressive morphological filtering

[0059] A single filtering window cannot accommodate outlier information at different scales in the original multibeam data. Therefore, this invention employs a gradually increasing filtering window to perform multi-level iterative operations on the surface grid. The filtering window grows linearly, and its iterative growth model is as follows:

[0060]

[0061] in, Indicates the first The morphological filtering window size during each iteration (in units of grid nodes, and the value is an odd number greater than zero). Indicates the size of the initial morphological filter window; The current iteration number ( ); This is the linear increment step size of the filtering window in each iteration.

[0062] In the In the next iteration, with Morphological structural elements are constructed for the dimensions, and morphological operations are performed on the elevation surfaces retained from the previous iteration to filter out convex elevation points at the current window scale, thereby generating the 1st generation. The reference smoothing datum plane for the next iteration.

[0063] 5) Marking suspected outliers

[0064] To achieve faithful processing of complex underwater terrain, after obtaining the reference smoothing datum surface, the algorithm is based on the terrain slope tolerance parameters pre-calculated in step 3. Real-time calculation of corresponding grid nodes Local adaptive elevation difference threshold Its calculation model is as follows:

[0065]

[0066] In the formula, For the first The local adaptive elevation difference threshold is obtained from the next iteration. The initial height difference threshold is the height difference threshold from the previous iteration. This is typically set as the nominal depth sounding error limit for a multibeam echo sounding system. and These are the filter window sizes for the current round and the previous round, respectively. The horizontal physical resolution of the surface grid is given. This computational model achieves local spatial adaptive dynamic adjustment of the filtering threshold: in In areas with significant topographic relief, the cumulative elevation difference threshold is rapidly amplified as the window size increases, effectively preventing real seabed protrusions from being misjudged as outliers; while... In smaller, flat seabed areas, the growth of the threshold is strictly suppressed, ensuring that residual large-scale clutter is accurately identified.

[0067] Finally, the previous elevation surface was compared with the first... The secondary reference smoothing datum surface at the target node The absolute value of the elevation difference at a point and the adaptive elevation difference threshold at that point Perform a one-by-one comparison. If the absolute value of the actual elevation difference is greater than... If so, the original data point is determined to be a "suspected outlier".

[0068] 6) Determine whether to terminate iterative filtering.

[0069] Considering that an excessively large filtering window might smooth out real terrain features, and to improve the algorithm's computational efficiency, an iteration termination condition is designed based on the adjacent elevation differences in the iteration results to avoid over-filtering due to an unreasonable filtering window setting. Therefore, the iteration termination condition is set as follows:

[0070]

[0071] in, This is the depth error threshold. During the iteration process, if the absolute value of the height difference between adjacent filtering results is less than the depth error threshold... If the filtering process is considered to have stabilized, the iteration process ends. Simultaneously, if the structure window reaches its set maximum value, the iterative filtering also terminates, and all marked outlier points from the iteration process are output. Otherwise, proceed to step 4.

[0072] 7) Generate terrain reference surface

[0073] Considering that the progressive morphological filtering process essentially smooths the terrain mesh, it may lose some detailed terrain surface features while marking suspected outliers. Therefore, instead of directly using the terrain surface mesh output by the filter, the marked suspected outlier points are placed on the original point cloud mesh. The mesh is then culled, and the resulting mesh holes are interpolated to form a new reference point cloud mesh, denoted as . .

[0074] 8) Outlier Detection and Removal

[0075] After the aforementioned multi-level progressive morphological iterative filtering, a smooth reference topographic mesh surface that closely approximates the real seabed is obtained. Macroscopic topographic slope information is then calculated based on this reference mesh surface. It should be noted that, unlike the topographic slope tolerance parameter calculated in step 3 based on the initial pre-smoothed mesh and local windows, which aims to dynamically control the increase of the elevation difference threshold in a single iteration, the topographic features extracted in this step are intended to serve as the basis for identifying outliers in the original point cloud. Therefore, this macroscopic topographic slope... Based on the reference surface after all noise has been removed, and the maximum filter window size at the termination of iterative filtering. This is obtained as a supporting neighborhood calculation.

[0076] The specific calculation method is as follows: taking any node on the surface of the reference terrain grid... Centered on, in a size of Within a large-scale neighborhood, calculate the rate of change of spatial elevation between the central node and each neighboring node, and take the maximum value as the terrain slope at that node. .

[0077] Based on the aforementioned macroscopic terrain slope information, the final outlier judgment threshold at each node in the grid is further calculated. Its calculation model is as follows:

[0078]

[0079] in, This is a scaling factor, typically set to 1.2-1.5 based on the acoustic reflection characteristics of the substrate and the complexity of the terrain in the survey area. The horizontal physical resolution of the point cloud mesh; This is the final filter window size.

[0080] After obtaining the gridded decision threshold, point-by-point detection is performed on the original discretely distributed multibeam data. This is based on the horizontal coordinates of the actual sounding points. By employing methods such as bilinear interpolation within the reference terrain grid, the reference depth information for the corresponding location is calculated. Similarly, the outlier judgment threshold at that point is obtained. .

[0081] Finally, by comparing the reference depth information with the actual measured depth information... To perform the final identification of outliers, the identification criteria are as follows:

[0082]

[0083] If the absolute difference between the measured depth and the reference depth of an actual sounding point is greater than the threshold, the sounding point is identified as an outlier and removed from the original point cloud; if it is less than or equal to the threshold, it is considered true seabed depth data and retained. By iterating through all original sounding points and applying the above identification criteria, efficient detection and accurate removal of outliers in multibeam bathymetry data are achieved, ultimately outputting cleaned, high-quality multibeam bathymetry point cloud data.

[0084] The embodiments of the present invention are described in detail below. These embodiments are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0085] To verify the effectiveness of this invention, the algorithm was simulated using the MATLAB platform. The multibeam bathymetry data used was measured depth point cloud data from a lake, with four measurement areas selected, each containing 20,000 bathymetry points. The parameter settings for the method of this invention are shown in Table 1.

[0086] Table 1. Parameter settings for the improved PMF algorithm

[0087] parameter Parameter value Preprocessing structure window Maximum structure window Depth error threshold Height scaling factor

[0088] The filtering result is as follows Figure 2-5 As shown, where, Figure 2 and Figure 4 This is the raw multibeam bathymetry data. Figure 3 and Figure 5 This is the filtered multibeam bathymetry data.

[0089] Figure 2 The underwater topography exhibits sloping characteristics, and its outliers display a characteristic arc-shaped distribution, commonly seen in the data of the same ping or several consecutive pings. This is usually caused by the sidelobe effect, meaning that the secondary beams of the multibeam system also receive reflected signals. These signals may originate from the topographic edges or other seabed features, but do not represent the true underwater topography. After filtering using the method of this invention, a total of 262 outlier points were detected. The multibeam data after removing outliers is as follows: Figure 3 As shown, the method of the present invention can accurately detect outliers in sloping terrain, reducing the rate of missed detections and false detections.

[0090] Figure 4 In the dataset, outliers are located beneath a relatively flat terrain surface and show significant separation from the main dataset, exhibiting considerable dispersion. This distribution may be due to multipath reflection, where sound waves undergo multiple reflections before reaching the transducer, generating multiple echoes and causing some outliers to appear below the actual terrain surface. After filtering using the method of this invention, a total of 342 outlier points were detected. The multibeam dataset after removing outliers is shown below. Figure 5As shown, the method of this invention can effectively identify outliers with large amplitudes, and is also quite sensitive to the detection of outliers with small amplitudes. This result is due to the progressive window iterative calculation, which enables the method to extract terrain change trends more accurately and form a more accurate reference surface. Meanwhile, improvements to the morphological operators allow the algorithm to effectively identify both positive and negative outliers.

[0091] Example 2, refer to Appendix Figure 1-3 An adaptive progressive morphological filtering method for outlier detection in multi-beam data includes the following steps:

[0092] Step 1: Obtain the raw multibeam bathymetry point cloud data and regularize it into an elevation grid surface.

[0093] Step 2: Perform morphological processing on the elevation grid surface, and extract local terrain features based on the processed elevation grid surface. The local terrain features include terrain slope tolerance parameters.

[0094] When performing morphological processing on the elevation grid surface, the closed-open and open-closed combined average morphological operator is used, and its calculation model is as follows:

[0095]

[0096] in, For data to be processed, For structural elements, select a rectangular structural element here; and f Opening and closing and closed-open form operators, and Let represent the opening and closing operations, respectively, and define them as follows:

[0097]

[0098] in, and These represent the erosion and dilation operations, respectively, and are defined as follows:

[0099]

[0100] The extraction of local terrain features based on the processed elevation grid surface specifically includes: using any target node on the smoothed elevation grid surface... A local slope calculation window is constructed centered on the target node; the spatial elevation change rate between the target node and each neighboring node within the local slope calculation window is calculated, and the maximum value is extracted as the local terrain slope. And an empirical threshold for the maximum topographic slope of the survey area is introduced. The slope of the local terrain After truncation, the final terrain slope tolerance parameter is obtained. .

[0101] Step 3: Construct a progressive morphological filtering process with an increasing filter window size. In a single iteration: Based on the current filter window size, use morphological operators to process the elevation grid surface retained in the previous round to obtain a reference smoothing datum surface; dynamically calculate the adaptive elevation difference threshold corresponding to the grid node using the terrain slope tolerance parameter; compare the grid surface before iteration with the reference smoothing datum surface, and mark the elevation of the grid node in the current iteration as a suspected outlier based on the adaptive elevation difference threshold.

[0102] The calculation model for dynamically calculating the adaptive elevation difference threshold corresponding to the grid node using the terrain slope tolerance parameter is as follows:

[0103]

[0104] In the formula, For the first The local adaptive elevation difference threshold is obtained from the next iteration. The initial height difference threshold is the height difference threshold from the previous iteration. This is typically set as the nominal depth sounding error limit for a multibeam echo sounding system. and These are the filter window sizes for the current round and the previous round, respectively. The horizontal physical resolution of the surface grid.

[0105] Step 4: Determine whether the filtering termination condition has been met. If so, terminate the iteration and output the final retained terrain reference surface.

[0106] Step 5: Extract macroscopic terrain slope information based on the terrain reference surface and calculate the global outlier judgment threshold. Specifically, this includes:

[0107] With the nodes on the terrain reference surface Centered on the node, within the neighborhood of the maximum filter window size at the termination of iterative filtering, the maximum rate of change of spatial elevation between the center node and each neighboring node is calculated to obtain the macroscopic terrain slope at that node. The global outlier threshold is calculated using the following formula. :

[0108]

[0109] in, This is a scaling factor, typically set to 1.2-1.5 based on the acoustic reflection characteristics of the substrate and the complexity of the terrain in the survey area. The horizontal physical resolution of the point cloud mesh; This is the final filter window size.

[0110] Step 6: Based on the horizontal position of the original discrete multibeam echo sounder point, obtain its corresponding reference elevation and global outlier judgment threshold; compare the depth deviation between the actual measured depth and the reference elevation, and detect and remove outliers in the point cloud data according to the global outlier judgment threshold. Specifically, this includes: based on the horizontal coordinates of the actual multibeam echo sounder point, using a spatial interpolation algorithm, interpolating in the terrain reference surface and the global outlier judgment threshold grid respectively to obtain the reference depth information corresponding to the point and the mapped global outlier judgment threshold.

Claims

1. A method for outlier detection in multi-beam data based on adaptive progressive morphological filtering, characterized in that, Includes the following steps: Step 1: Acquire the raw multibeam bathymetry point cloud data and regularize it into an elevation grid surface; Step 2: Perform morphological processing on the elevation grid surface, and extract local terrain features based on the processed elevation grid surface. The local terrain features include terrain slope tolerance parameters. Step 3: Construct a progressive morphological filtering process with an increasing filter window size. In a single iteration: Based on the current filter window size, use morphological operators to process the elevation grid surface retained in the previous round to obtain a reference smoothing datum surface; dynamically calculate the adaptive elevation difference threshold corresponding to the grid node using the terrain slope tolerance parameter; compare the grid surface before iteration with the reference smoothing datum surface, and mark the suspected outlier of the current iteration's grid node elevation based on the adaptive elevation difference threshold. Step 4: Determine whether the filtering termination condition has been met. If so, terminate the iteration and output the final retained terrain reference surface. Step 5: Extract macroscopic terrain slope information based on the terrain reference surface, and calculate the global outlier judgment threshold; Step 6: Based on the horizontal position of the original discrete multibeam echo sounder point, obtain its corresponding reference elevation and global outlier judgment threshold; compare the depth deviation between the actual measured depth and the reference elevation, and detect and remove outliers in the point cloud data according to the global outlier judgment threshold.

2. The method according to claim 1, characterized in that: When performing morphological processing on the elevation grid surface, the closed-open and open-closed combined average morphological operator is used, and its calculation model is as follows: in, For data to be processed, For structural elements, select a rectangular structural element here; and They are respectively Opening and closing and closed-open form operators, and Let represent the opening and closing operations, respectively, and define them as follows: in, and These represent the erosion and dilation operations, respectively, and are defined as follows:

3. The method according to claim 1, characterized in that: Step 2, which involves extracting local terrain features based on the processed elevation grid surface, specifically includes: using any target node on the smoothed elevation grid surface... A local slope calculation window is constructed centered on the target node; the spatial elevation change rate between the target node and each neighboring node within the local slope calculation window is calculated, and the maximum value is extracted as the local terrain slope. And an empirical threshold for the maximum topographic slope of the survey area is introduced. The slope of the local terrain After truncation, the final terrain slope tolerance parameter is obtained. .

4. The method according to claim 1, characterized in that: The calculation model for dynamically calculating the adaptive elevation difference threshold corresponding to the grid node using the terrain slope tolerance parameter is as follows: In the formula, For the first The local adaptive elevation difference threshold is obtained from the next iteration. The initial height difference threshold is the height difference threshold from the previous iteration. This is typically set as the nominal depth sounding error limit for a multibeam echo sounding system. and These are the filter window sizes for the current round and the previous round, respectively. The horizontal physical resolution of the surface grid.

5. The method according to claim 1, characterized in that: The step of extracting macroscopic terrain slope information based on the terrain reference surface and calculating the global outlier judgment threshold specifically includes: using nodes on the terrain reference surface... Centered on the node, within the neighborhood of the maximum filter window size at the termination of iterative filtering, the maximum rate of change of spatial elevation between the center node and each neighboring node is calculated to obtain the macroscopic terrain slope at that node. The global outlier threshold is calculated using the following formula. : in, This is a scaling factor, typically set to 1.2-1.5 based on the acoustic reflection characteristics of the substrate and the complexity of the terrain in the survey area. The horizontal physical resolution of the point cloud mesh; This is the final filter window size.

6. The method according to claim 1, characterized in that: Step 6 specifically includes: based on the horizontal coordinates of the actual multibeam bathymetry point, using a spatial interpolation algorithm, interpolating in the terrain reference surface and the global outlier judgment threshold grid respectively, to obtain the reference depth information corresponding to the point and the mapped global outlier judgment threshold.

Images