A multi-beam point cloud filtering method based on adaptive terrain trend

Abstract

The present invention relates to the field of multi-beam bathymetry technology, and in particular to a multi-beam point cloud filtering method for adaptive terrain trends. The method comprises the following steps: obtaining original multi-beam Ping data, extracting the attributes of each beam point in each Ping, and obtaining a beam point attribute data set; performing multi-threshold condition screening on the beam point attribute data set, and eliminating abnormal beam points to obtain a preliminary filtered point cloud data set; performing linear interpolation on the inter-Ping gaps in the eliminated area of ​​the preliminary filtered point cloud data set to obtain a complete point cloud data set; extracting the central beam point data of three adjacent survey lines based on the complete point cloud data set, and performing coordinate transformation and fitting modeling on the central beam point data to obtain a long-wave term function model. The present invention accurately eliminates anomalies, repairs gaps, and corrects system errors through an automated process, thereby achieving high-precision restoration and real-time processing of multi-beam point cloud data, and improving the consistency and reliability of bathymetry results.

Description

Technical Field

[0001] The present invention relates to the technical field of multi-beam bathymetry, and in particular to a multi-beam point cloud filtering method for adaptive terrain trends. Background Art

[0002] Multibeam bathymetry systems have become the mainstream of current bathymetry equipment. Due to their comprehensive coverage, full water depth, high accuracy, and high resolution, they are widely used in marine scientific research and engineering applications, including underwater topography and seafloor geomorphology surveys, and underwater target detection. However, after calculating the spatial beam incidence angle, ray tracking, and beam footprint, the multibeam bathymetry system initially calculates the coordinates and water depth of each ping beam point in a specified geodetic coordinate system and vertical datum. These data still contain systematic errors caused by gross measurement errors and sound velocity errors. Furthermore, in complex ocean environments, both sound velocity measurement errors and installation errors can systematically affect the depth calculation of the measured points, often resulting in a "crying face" or "smiling face" effect on the measured topography. Traditional multibeam bathymetry data processing typically involves manual removal of gross error points and sound points with large systematic errors. However, manual data post-processing is impractical for real-time monitoring. Summary of the Invention

[0003] Based on this, it is necessary for the present invention to provide a multi-beam point cloud filtering method that is adaptive to terrain trends to solve at least one of the above technical problems.

[0004] To achieve the above object, a multi-beam point cloud filtering method for adaptive terrain trend includes the following steps:

[0005] Step S1: Obtain the original multi-beam Ping data, extract the attributes of each beam point in each Ping, and obtain the beam point attribute dataset;

[0006] Step S2: Perform multi-threshold screening on the beam point attribute dataset and remove abnormal beam points to obtain a preliminary filtered point cloud dataset;

[0007] Step S3: performing linear interpolation on the inter-Ping gaps in the removed area of ​​the preliminary filtered point cloud dataset to obtain a complete point cloud dataset;

[0008] Step S4: extracting the central beam point data of the three adjacent survey lines based on the complete point cloud data set, and performing coordinate transformation and fitting modeling on the central beam point data to obtain a long-wave term function model; extracting the edge beam point data based on the complete point cloud data set, establishing an edge beam point error fitting function model, and obtaining a short-wave term function model;

[0009] Step S5: performing depth correction on the edge beam points in the complete point cloud dataset based on the long-wave term function model and the short-wave term function model to obtain an adaptive trend filtered point cloud dataset;

[0010] Step S6: Build a grid model based on the adaptive trend filtering point cloud data set, perform depth sounding error statistics on the grid points, and obtain the filtering accuracy evaluation result.

[0011] This method can automatically address systematic errors and gross errors in multi-beam bathymetry systems caused by sound velocity errors, equipment installation errors, and extreme angles of incidence without manual intervention. It effectively addresses technical challenges such as low response efficiency, poor real-time processing capabilities, and the erroneous deletion of terrain details in traditional manual filtering methods. By constructing a unified format for beam point attribute datasets and implementing multi-threshold screening, this method improves the targetedness and accuracy of anomaly removal while maintaining structural integrity. Linear interpolation of inter-ping points restores the continuity of point cloud data, avoiding terrain fragmentation caused by missing points. Furthermore, by combining the long-wave and short-wave trends reflected by the central and edge beam points, an error trend function is constructed and corrected in different regions. This suppresses systematic errors in the edge beams while preserving detailed seafloor topography, addressing the "crying face" and "smiling face" terrain distortion phenomena. The grid model and error assessment mechanism constructed on this basis provide a basis for objective quality verification of point cloud data, ensuring the final results are highly stable, consistent, and accurate. This method is particularly suitable for marine surveying and seabed engineering scenarios where the accuracy and consistency of bathymetry results are extremely demanding. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Other features, objects and advantages of the present invention will become more apparent upon reading the detailed description of non-limiting embodiments thereof made with reference to the following drawings:

[0013] Figure 1 A schematic flow chart of the steps of a multi-beam point cloud filtering method for adaptive terrain trend according to the present invention;

[0014] Figure 2 for Figure 1 Detailed step flow diagram of step S1;

[0015] Figure 3 for Figure 1 Detailed step flow diagram of step S3;

[0016] Figure 4 This is a topographic trend map caused by multi-beam system errors;

[0017] Figure 5 This is a schematic diagram based on trend surface correction;

[0018] Figure 6It is the processed local terrain detail map;

[0019] Figure 7 This is a table of accuracy evaluation results in an embodiment of the present invention;

[0020] Figure 8 is an error histogram in an embodiment of the present invention;

[0021] in, Figure 5 In the figure, L1, L2, and L3 are three adjacent survey lines, represented by red, yellow, and blue respectively. The three dotted lines represent the lines connecting the seabed points measured by the three survey lines. DETAILED DESCRIPTION

[0022] The following is a clear and complete description of the technical method of the present invention in conjunction with the accompanying drawings. It is obvious that the embodiments described are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative efforts are within the scope of protection of the present invention.

[0023] In addition, the accompanying drawings are merely schematic illustrations of the present invention and are not necessarily drawn to scale. Identical reference numerals in the figures denote identical or similar parts, and thus repetitive descriptions thereof will be omitted. Some of the block diagrams shown in the accompanying drawings are functional entities that do not necessarily correspond to physically or logically separate entities. These functional entities may be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor and / or microcontroller approaches.

[0024] It should be understood that although the terms "first," "second," and the like may be used herein to describe various elements, these elements should not be limited by these terms. These terms are used solely to distinguish one element from another. For example, a first element may be referred to as a second element, and similarly, a second element may be referred to as a first element, without departing from the scope of the exemplary embodiments. The term "and / or" as used herein includes any and all combinations of one or more of the listed associated items.

[0025] To achieve this, please refer to Figures 1 to 8 The present invention provides a multi-beam point cloud filtering method for adaptive terrain trend, the method comprising the following steps:

[0026] Step S1: Obtain the original multi-beam Ping data, extract the attributes of each beam point in each Ping, and obtain the beam point attribute dataset;

[0027] Step S2: Perform multi-threshold screening on the beam point attribute dataset and remove abnormal beam points to obtain a preliminary filtered point cloud dataset;

[0028] Step S3: performing linear interpolation on the inter-Ping gaps in the removed area of ​​the preliminary filtered point cloud dataset to obtain a complete point cloud dataset;

[0029] Step S4: extracting the central beam point data of the three adjacent survey lines based on the complete point cloud data set, and performing coordinate transformation and fitting modeling on the central beam point data to obtain a long-wave term function model; extracting the edge beam point data based on the complete point cloud data set, establishing an edge beam point error fitting function model, and obtaining a short-wave term function model;

[0030] Step S5: performing depth correction on the edge beam points in the complete point cloud dataset based on the long-wave term function model and the short-wave term function model to obtain an adaptive trend filtered point cloud dataset;

[0031] Step S6: Build a grid model based on the adaptive trend filtering point cloud data set, perform depth sounding error statistics on the grid points, and obtain the filtering accuracy evaluation result.

[0032] Preferably, step S1 includes the following steps:

[0033] Step S11: Using a survey vessel equipped with a multi-beam bathymetry system to perform underwater bathymetry scanning of the survey area to obtain raw multi-beam Ping data;

[0034] Step S12: Decode and structurally analyze the original multi-beam Ping data, extract the Ping number, timestamp, heading, attitude, sound speed, each beam number and three-dimensional coordinate information, and obtain the Ping structure data set;

[0035] Step S13: performing spatial inversion calculation on each beam point in the Ping structure data set based on the attitude and sound speed to obtain a beam point spatial coordinate data set;

[0036] Step S14: performing attribute extraction processing on the beam point spatial coordinate data set to obtain a beam point attribute record set, where the attributes include the depth value, beam angle, beam number, quality factor, and Ping number of each beam point;

[0037] Step S15: constructing a beam point attribute data set in a unified format based on the beam point attribute record set.

[0038] In an embodiment of the present invention, a survey vessel equipped with a Norbit iWBMS multi-beam bathymetry system is first used to perform a full-coverage underwater scanning operation within a set survey area along equally spaced survey lines. The survey vessel's speed is maintained in the range of 1.2 to 1.5 knots, the survey line spacing is set to 20 meters, the equipment operating frequency is fixed at 400kHz, and the transmit beam angle is set to 140° to ensure that the beam covers the edge areas on both sides of the survey line. During the operation, the system continuously records Ping sequence data, collecting no less than 10 sets of Ping information per second, and sequentially forms raw multi-beam Ping data. Subsequently, a decoding operation is performed on the raw multi-beam Ping data, and a Norbit dedicated data unpacking program is used to perform structured parsing on the recorded binary data stream according to the Ping unit structure format. The Ping number, timestamp, measurement heading, attitude information (pitch angle, roll angle, heading angle), sound velocity profile data, all beam numbers in each Ping, and their corresponding initial ranging data are extracted one by one. Through this structured parsing process, a Ping structure data set is formed. Based on the P The ing structured dataset is constructed. The ranging data, attitude information, and sound speed information of each beam point are used as input. The spatial position inversion calculation of all beam points is performed one by one using the three-dimensional spatial rotation transformation relationship between the beam launch angle and the attitude angle. A fixed sound speed layered interpolation model is used in the calculation process. The sound speed profile is measured by the CTD instrument before the survey area starts. The inversion calculation uses the matrix transformation method to solve the three-dimensional coordinates of each beam point in the measurement coordinate system to form the beam point spatial coordinate dataset. After obtaining the beam point spatial coordinate dataset, an attribute extraction operation is performed. For each beam point, the calculated vertical depth value, the beam launch angle recorded in the original Ping, the beam number, the signal quality factor, and the Ping number to which it belongs are extracted. These records are organized into structured records according to a unified field format to form a beam point attribute record set. Finally, all records in the beam point attribute record set are summarized and sorted and merged according to the "Ping number-beam number" dual index structure to construct a beam point attribute dataset in a unified format.

[0039] Preferably, step S13 includes the following steps:

[0040] Step S131: extracting the original bathymetric data of the current beam point from the Ping structure data set to obtain the bathymetric data of the beam point;

[0041] Step S132: performing attitude compensation calculation on the beam point sounding data based on the attitude to obtain attitude compensated sounding data;

[0042] Step S133: performing sound ray tracking correction on the attitude-compensated sounding data according to the sound velocity to obtain sound velocity-corrected sounding data;

[0043] Step S134: Perform spatial inversion calculation based on the sound velocity corrected bathymetric data to obtain a beam point spatial coordinate data set.

[0044] In an embodiment of the present invention, the original bathymetric data of each beam point is first extracted from the Ping structure data set, which includes field information such as Ping number, beam number, beam launch angle, ranging value and sound intensity quality factor. The ranging value is the initial distance data calculated by the Norbit iWBMS system based on the echo delay and the current sound speed, and the unit is meter. The system transmission frequency is 400kHz, and the Ping interval is controlled within 0.1 second. After the extraction is completed, the attitude compensation is performed on each beam point based on the attitude data (including pitch angle, roll angle, and heading angle) synchronously recorded by the system. The original ranging direction is converted from the transducer coordinate system to the geographic reference coordinate system using a three-dimensional rotation matrix transformation method, wherein the pitch angle range is controlled within ±10° and the roll angle range is controlled within ±5°. The output after compensation is attitude-compensated bathymetric data, and the original beam number information is still retained. The attitude-compensated bathymetric data is then processed according to the sound velocity profile data of the survey area. Acoustic ray tracking correction processing is performed. The sound velocity profile data is measured by the SBE19Plus CTD instrument. The sound velocity variation curve is constructed using layered linear interpolation. The beam propagation path is fitted with a multi-segment broken line based on Snell's law. The sound velocity gradient in each segment does not exceed the variation range of 2m / s. The corrected distance value is used as the sound velocity corrected bathymetry data. Finally, the sound velocity corrected bathymetry data is converted into three-dimensional spatial coordinates. Specifically, the corrected distance value and corresponding launch angle of each beam point are converted into polar coordinates. Then, combined with the current geographical location and attitude angle of the survey ship, the three-dimensional position coordinates (X, Y, Z) of each beam point in the measurement coordinate system are calculated using the formula from polar coordinates to rectangular coordinates. The X direction points to the navigation direction of the survey line, the Y direction points to the lateral direction of the hull, and the Z direction points vertically downward. The output is the beam point spatial coordinate dataset.

[0045] Preferably, in step S2, performing multi-threshold condition screening on the beam point attribute data set and eliminating abnormal beam points includes:

[0046] Set the minimum depth and maximum depth thresholds, maximum beam angle threshold, maximum width threshold, maximum beam point slope and quality factor optional range based on the operating water depth range and regional terrain of the survey area;

[0047] Determine the depth range of each beam point in the beam point attribute dataset, and remove beam points whose depth values ​​are not within the minimum depth and maximum depth thresholds to obtain a depth screening dataset;

[0048] The beam angle range of each beam point in the beam point attribute dataset is judged, and beam points with beam angles greater than the maximum beam angle threshold are eliminated to obtain an angle screening dataset;

[0049] For each beam point in the beam point attribute dataset, the horizontal distance from the beam point to the central beam point of the ping is calculated, and beam points whose horizontal distance exceeds the maximum width threshold are eliminated to obtain the width screening dataset.

[0050] The depth slope of each beam point in the beam point attribute dataset and its adjacent beam points is calculated, and beam points with slopes greater than the maximum beam point slope are eliminated to obtain a slope screening dataset.

[0051] Perform quality factor interval judgment on each beam point in the beam point attribute dataset, and eliminate beam points with quality factors below the lower limit or above the upper limit to obtain a quality screening dataset;

[0052] Perform an adjacency check on each beam point in the beam point attribute dataset, and remove isolated beam points with missing adjacent ping numbers or left and right adjacent beam numbers to obtain an adjacency screening dataset.

[0053] The depth filtering dataset, angle filtering dataset, width filtering dataset, slope filtering dataset, quality filtering dataset and adjacency filtering dataset are cross-fused, and the beam points that meet all the filtering conditions are retained to obtain the preliminary filtered point cloud dataset.

[0054] In an embodiment of the present invention, the process of performing multi-threshold conditional screening on a beam point attribute dataset first sets a specific screening parameter range, the operating water depth range is set to between 5 meters and 60 meters, the maximum threshold of the beam angle is set to 65 degrees, the maximum Ping width is limited to a horizontal distance of 25 meters around the central beam point of the survey line, the maximum threshold of the beam point depth slope is limited to 0.35, and the quality factor range is set to between 80 and 255; first, a depth screening operation is performed, and the depth value field of each beam point in the beam point attribute dataset is traversed in turn, and beam points with depths less than 5 meters or greater than 60 meters are eliminated, and a depth-screened dataset is output; then an angle screening operation is performed, and the beam angle field in each beam point record is read and compared with the set threshold of 65 degrees, and beam points with beam angles greater than 65 degrees are eliminated, and an angle-screened dataset is output; then a width screening operation is performed, and the central beam point with the beam number in the middle of each Ping is first determined as a reference point, and the horizontal distance of other beam points in the Ping to the central beam point in the horizontal coordinate direction is calculated, and the horizontal distance is calculated using the Euclidean distance formula. Projection length. If the distance is greater than 25 meters, the beam point is removed and the width screening data set is output. Then the slope screening process is performed. For the beam point sequence in each Ping, the center difference method is used to calculate the ratio of the depth difference between each beam point and its two adjacent beam points to the beam spacing as the slope value. If the calculated result is greater than 0.35, the beam point is removed and the slope screening data set is output. Then the quality factor screening operation is performed. The quality factor field of each beam point is read. If its value is less than 80 or greater than 255, then The point is removed and the quality-filtered dataset is output. Finally, the adjacency filtering operation is performed. For each beam point in the beam point attribute dataset, the neighboring points in the two preceding and following Pings and the two left and right beam numbers are queried in both the Ping number and the beam number directions to see if there are any. If the number of neighboring points is less than 3, the beam point is considered an isolated point and removed, and the adjacency-filtered dataset is output. The six types of filtered datasets are number-aligned and cross-compared, and only the beam points that have not been removed in the six filtering results are retained, finally forming the preliminary filtered point cloud dataset.

[0055] Preferably, step S3 includes the following steps:

[0056] Step S31: Based on the high-frequency sampling capability of the multi-beam bathymetry system and the Ping structure information, beam number matching is performed on the voids in the preliminary filtered point cloud dataset to obtain a set of beam point pairs to be interpolated;

[0057] Step S32: performing spatial linear fitting calculation on each pair of beam points in the set of interpolation beam point pairs, obtaining interpolation coefficients, and obtaining an interpolation parameter set;

[0058] Step S33: performing depth recovery on each interpolation target point in the interpolation parameter set, constructing a complete spatial coordinate record, and obtaining an interpolation point cloud dataset;

[0059] Step S34: The preliminary filtered point cloud dataset and the interpolated point cloud dataset are jointly merged to construct a complete point cloud dataset.

[0060] In the embodiment of the present invention, for the inter-ping structure missing areas in the preliminary filtered point cloud data set, first, based on the high-frequency sampling capability of the Norbit iWBMS multi-beam bathymetry system and the Ping transmission rate of more than 10 times per second, the Ping structure information is used to traverse the Ping number sequence one by one, and the beam number in each Ping is checked for integrity. By judging whether the same beam number exists continuously in adjacent Pings, the point position where a beam number is missing in the Ping number sequence is identified, and the coordinate information of the beam point with the same beam number in the two adjacent Pings before and after is recorded, extracted as the corresponding point pair of the beam number, and output as the set of beam point pairs to be interpolated; then, linear fitting processing is performed on the three-dimensional spatial coordinates of each pair of beam points in the set, and a vector linear interpolation method is used to construct a straight line connecting the two points before and after, and the distance proportional interpolation formula is used. ,in and is the depth value of the front and rear beam points, is the Ping distance of the current interpolation point relative to the previous point, The total Ping distance between the two points is calculated, and the depth value of the point where the current missing Ping number is located is calculated respectively to obtain the target depth parameter after interpolation, and the interpolation parameter set is output; after obtaining the interpolation parameter set, according to the Ping number, beam number, interpolation depth value and the emission direction vector corresponding to the beam angle of the interpolation target point, combined with the original attitude record, the spatial coordinate recovery operation is performed. First, the interpolation depth value and beam angle are converted into polar coordinate form, and then the polar coordinate value is converted into the rectangular coordinate value in the measurement coordinate system using the attitude angle and geographic reference coordinates of the survey ship at the corresponding Ping moment. The three-dimensional spatial coordinates of the interpolation target point are calculated to obtain the interpolated point cloud dataset; finally, the interpolated point cloud dataset and the original preliminary filtered point cloud dataset are combined in a dual-index structure according to the order of Ping number and beam number, and the missing beam position is inserted to maintain the continuity of the Ping structure and the integrity of the beam number, so as to construct a complete point cloud dataset.

[0061] Preferably, step S32 includes the following steps:

[0062] Step S321: extracting the three-dimensional spatial coordinate information and the corresponding Ping number for each pair of beam points in the interpolation beam point pair set, constructing a beam point pair spatial coordinate record, and obtaining a spatial point pair dataset;

[0063] Step S322: extracting the horizontal coordinate position of the target interpolation point for each group of beam points in the spatial point pair data set, calculating its relative position scale coefficient, and obtaining an interpolation coefficient data set;

[0064] Step S323: performing correlation matching on the interpolation coefficient dataset and the spatial point pair dataset, performing one-dimensional depth linear interpolation, and obtaining an interpolation depth value set;

[0065] Step S324: Synchronously record the target Ping number and beam number information for the interpolation depth value set and the interpolation coefficient data set to construct an interpolation parameter set.

[0066] In the embodiment of the present invention, a three-dimensional coordinate extraction operation is first performed on each pair of beam points in the interpolation beam point pair set. The X, Y, and Z three-dimensional spatial coordinate values ​​of the two points and the corresponding Ping number and beam number information are read respectively. The data are recorded in a structured format as a set of point pair data containing a "front point-back point" spatial relationship. At the same time, the increasing relationship between the front point and the back point in each point pair in the Ping number is clarified, and a complete beam point pair spatial coordinate record is constructed, and a spatial point pair data set is output. On this basis, the horizontal coordinate position of the target interpolation point is further extracted for each group of beam points in the spatial point pair data set, that is, the relative position between the Ping number of the interpolation target point and the Ping number of the two end points, and the relative position proportional coefficient is calculated. ,in Indicates the interpolation target Ping number, and are the Ping numbers of the front and rear beam points respectively. This coefficient is used to measure the linear distance ratio of the target point between the two points. All calculation results are organized and output as an interpolation coefficient dataset. Next, the interpolation coefficient dataset is matched with the spatial point pair dataset in a one-to-one correspondence. A one-dimensional linear interpolation operation is performed on the Z-direction depth value of each point pair using the formula ,in and are the depth values ​​of the front and rear beam points respectively, The calculated depth of the target interpolation point is output after group calculation. Finally, each interpolation result in the interpolation depth value set is synchronized with the corresponding interpolation coefficient data to record the target Ping number, beam number, interpolation coefficient and depth value calculated by linear interpolation, and a unified integration is constructed to construct the interpolation parameter set. It is also clear that the three-dimensional coordinates of each interpolation point still need to be combined with the original beam angle and attitude information in the subsequent steps for coordinate inversion conversion.

[0067] It is particularly important that step S33 includes the following steps:

[0068] Step S331: extracting the interpolation depth value, interpolation lateral coordinate position, Ping number and beam number of each target point in the interpolation parameter set, and constructing an interpolation point positioning information set;

[0069] Step S332: performing attitude constraint and sound ray direction calculation on each interpolation point based on the interpolation point positioning information set and the preset original survey vessel trajectory data to obtain the three-dimensional positioning parameters of the interpolation point in the measurement coordinate system;

[0070] Step S333: performing geodetic coordinate system transformation on the three-dimensional positioning parameter set in the interpolation point measurement coordinate system to generate three-dimensional space coordinates that conform to the measurement area coordinate reference, and constructing a three-dimensional coordinate data set of the interpolation point;

[0071] Step S334: Structurally encapsulate the interpolation point three-dimensional coordinate dataset, combine the corresponding Ping number and beam number information, construct a standard interpolation point cloud record, and obtain the interpolation point cloud dataset.

[0072] In the embodiment of the present invention, first, the interpolation depth value, lateral coordinate position, Ping number and beam number information of each target point in the interpolation parameter set are extracted. The lateral coordinate position is the lateral offset value along the survey line direction calculated according to the Ping interval and the target Ping number position, and is stored in a structured record format to construct an interpolation point positioning information set. Then, based on the interpolation point positioning information set, the attitude data (including pitch angle, roll angle, and heading angle) recorded with the corresponding Ping number in the original survey ship trajectory data and the dead reckoning result are combined. The beam angle corresponding to the beam number to which the interpolation point belongs and the interpolation depth value are used to perform spatial direction constraints. The three-dimensional rotation matrix and trigonometric function transformation method are used to calculate the X, Y, and Z three-dimensional positioning parameters of each interpolation point in the measurement coordinate system (shipborne coordinate system), where the X axis is the navigation direction, the Y axis is the lateral direction, and the Z axis is the vertical depth. The input is The three-dimensional positioning parameter set of the interpolation point in the measurement coordinate system is obtained; then, the coordinate transformation from the measurement coordinate system to the geodetic coordinate system is performed on the three-dimensional positioning parameter set in sequence, and the reference position and attitude angle recorded by the ship-borne GNSS system are used as the conversion reference. The local three-dimensional coordinates in the measurement coordinate system are converted into global three-dimensional coordinates in the WGS84 reference through the Euler angle transformation method. The accuracy of the conversion result is controlled, and the error is limited to within ±0.05 meters. The three-dimensional spatial coordinates that conform to the projection reference of the survey area are generated and output as the interpolation point three-dimensional coordinate dataset; finally, the interpolation point three-dimensional coordinate dataset is structured and packaged, and the three-dimensional coordinate value of each point is integrated and stored with its corresponding Ping number, beam number, interpolation depth value, beam angle, attitude angle and lateral projection position, and they are uniformly arranged in the order of Ping number and beam number to construct a standard interpolation point cloud record, which is output as the interpolation point cloud dataset.

[0073] Preferably, in step S4, extracting the central beam point data of three adjacent survey lines based on the complete point cloud data set, and performing coordinate transformation and fitting modeling on the central beam point data includes:

[0074] Perform index matching on the Ping data of three adjacent survey lines selected from the complete point cloud dataset, extract the central beam point in each Ping, and obtain the central beam point dataset;

[0075] Perform coordinate system transformation on the central beam point data set to establish a local coordinate system with a certain point on the middle survey line as the origin, the starboard side as the y-axis, and the vertical downward direction as the z-axis, and obtain the transformed central beam point coordinate set;

[0076] The least squares polynomial fitting is performed based on the central beam point coordinate set to construct a long-wave term function model of the terrain trend.

[0077] In this embodiment of the present invention, data from three adjacent survey lines are first selected from the complete point cloud dataset based on the order of the survey line number and the Ping number. Each survey line is set to contain no less than 500 Ping records, ensuring that the lateral spacing between survey lines is less than 1.5 meters. A Ping index sequence is established for each of the three selected survey lines. After synchronization and alignment based on the Ping number, the central beam point data is extracted from each Ping record. The central beam point is defined as the beam point located in the middle of the beam number sequence. The three-dimensional spatial coordinates and corresponding Ping number of this point are extracted to construct a central beam point dataset. The structure contains the X, Y, and Z coordinate values ​​of each point, as well as its survey line number and Ping number information. Subsequently, A coordinate system transformation operation is performed on the central beam point data set. A central beam point in the middle area of ​​the middle survey line is selected as the origin of the local coordinate system. The survey line direction of the origin is defined as the X-axis, the starboard direction of the hull is defined as the Y-axis, and the vertical downward direction is the Z-axis. The original coordinates are rotated and translated by constructing a direction cosine matrix. All central beam point coordinates are converted to a local coordinate system based on the origin to obtain the transformed central beam point coordinate set. All coordinate data are kept in meters and the accuracy is controlled within 0.01 meters. On this basis, the least squares fitting method is used to perform polynomial surface fitting on the relationship between the Z value and the X and Y coordinates. The fitting function is set to , solve the normal equations for the X, Y, and Z data in the central beam point coordinate set to obtain the fitting coefficients to , which is used to represent the overall trend deformation of the terrain in the survey area. The surface fitting function is the long-wave function model of the terrain trend.

[0078] Preferably, in step S4, extracting edge beam point data based on the complete point cloud data set and establishing an edge beam point error fitting function model includes:

[0079] Extract edge beam points from three adjacent measurement lines in the complete point cloud dataset and construct an edge beam point dataset;

[0080] Perform local coordinate transformation on the edge beam point data set to construct the edge beam point local coordinate set;

[0081] Based on the local coordinate set of the edge beam point, the m-order polynomial fitting is performed using the least squares method to construct the shortwave term function model of the terrain error.

[0082] In the embodiment of the present invention, three adjacent survey lines are first selected from the complete point cloud dataset, with the lateral spacing between the survey lines not exceeding 1.5 meters. All Ping numbers in each survey line are traversed, and the leftmost and rightmost n beam points (n ranges from 3 to 5) in each Ping are extracted as edge beam points according to the beam number index. The extracted content includes the three-dimensional spatial coordinates, Ping number, beam number and survey line number of each point, and is organized into an edge beam point dataset. The edge beam point dataset is then transformed into a coordinate system using the same system as the central beam point in step S4. The local coordinate system is constructed in a consistent way, that is, a central beam point in the middle of the middle survey line is selected as the origin, the direction along the survey line is defined as the X axis, the direction to the starboard of the hull is the Y axis, and the vertical downward is the Z axis. Based on the local coordinate system, the three-dimensional coordinates of the edge beam points are rotated and translated, and the edge beam points are uniformly mapped to the local coordinate system to form a local coordinate set of the edge beam points; then, the data in the coordinate set are subjected to terrain error fitting according to the spatial distribution relationship, with the X and Y coordinates as independent variables and the Z coordinate as the dependent variable, and the m-order two-dimensional polynomial surface function is set as follows: ,in is the polynomial order, and the recommended value is 3 to 4. The normal equation is constructed using the least squares principle to minimize the sum of squares of the fitting residuals of all points, and the coefficients of each polynomial are obtained by matrix solution. , and the final two-dimensional polynomial function is defined as the shortwave term function model of terrain error.

[0083] Preferably, step S5 includes the following steps:

[0084] Step S51: Perform local coordinate transformation on the edge beam points in the complete point cloud data set to obtain the local coordinate values ​​of the edge beam points. ;

[0085] Step S52: Based on the shortwave term function model, the The error trend calculation is performed on the value to obtain the depth detail change of the edge beam point;

[0086] Step S53: Based on the long-wave function model, the The trend surface estimation is performed on the value to obtain the trend prediction depth value;

[0087] Step S54: Calculating a target corrected depth value based on the depth detail variation and the trend predicted depth value;

[0088] Step S55: performing depth correction on each edge beam point using the target corrected depth value to construct a corrected edge beam point data set;

[0089] Step S56: Merge the corrected edge beam point dataset with the unmodified central and non-edge beam points in the complete point cloud dataset, reconstruct the beam point structure of all Pings, and obtain an adaptive trend filtered point cloud dataset.

[0090] In the embodiment of the present invention, based on the complete point cloud data set obtained by the Norbit iWBMS multi-beam bathymetry system, the edge beam points at both ends of each survey line are first identified, and a local coordinate transformation operation is performed on their three-dimensional coordinate data. A central beam point in the middle of the survey line is set as the origin, and a local coordinate system is established, in which the starboard direction is the y-axis and the vertical downward direction is the z-axis. The edge beam points are rotated by two-dimensional coordinates. Perform transformation to local coordinates Then, based on the edge beam error trend function obtained by the least squares fitting method in the early stage, the error of each edge beam point is input. Value, directly calculate its corresponding error trend depth value through polynomial function and with the original depth Subtract and get the corresponding depth detail change This process relies on the original sounding distribution and position parameters of the edge beam points, without the need to introduce high-order models; then, based on the terrain trend function model established by the central beam point, the same The value is used as input to calculate the trend surface prediction depth. , which represents the reference depth of the ideal seabed topography at that location; and Perform algebraic summation to obtain the corrected depth value of each edge beam point ;Will Value replaces the original Value, combined with the horizontal coordinate of the original point And longitudinal distance information, reconstruct the new coordinates of each edge beam point in the local coordinate system, and transform it back to the three-dimensional coordinates in the original measurement area coordinate system with the local coordinate transformation parameters to form a corrected edge beam point data set; finally, the corrected edge beam point data set is paired with the uncorrected central beam point and non-edge beam point data sets in the complete point cloud data set for numbering and structural reconstruction to ensure that the Ping number and beam number information are consistent, and all beam points are merged into structured point cloud data in the original measurement order to complete the construction of the adaptive trend filtering point cloud data set. The input range should be limited to [-50m, 50m]. The calculation needs to ensure that the interpolation points are within the range of the existing edge data set, and all polynomial functions should be limited to third order or less to ensure the numerical stability and continuity of the fitted terrain trend and error trend in the local area.

[0091] Preferably, step S6 includes the following steps:

[0092] Step S61: performing spatial range analysis on the adaptive trend filtering point cloud data set to determine the measurement area boundary range and resolution parameters, and obtain a grid division parameter set;

[0093] Step S62: performing spatial projection and discrete mapping on the adaptive trend filtering point cloud dataset based on the grid division parameter set, constructing an equidistant regular grid model, and obtaining a grid model dataset;

[0094] Step S63: Comparing and analyzing the results of different survey lines in the overlapping area of ​​the grid model dataset, calculating the sounding difference of the same grid point by locating the repeated coverage area of ​​the same grid point in multiple survey lines, and obtaining a sounding difference dataset;

[0095] Step S64: performing statistical analysis on the sounding difference data set, calculating the difference mean, absolute average, maximum, minimum and sounding root mean square error, and obtaining an error statistical index table;

[0096] Step S65: performing frequency distribution statistics on the depth difference data set, constructing error histogram data, and obtaining error distribution layer data;

[0097] Step S66: Generate a filtering accuracy evaluation result based on the error statistical indicator table and the error histogram data set.

[0098] In the embodiment of the present invention, first, statistical operations are performed on all the beam point coordinates in the adaptive trend filtering point cloud data set obtained by step S56, and the coordinates of the beam points are extracted. and The minimum and maximum values ​​in the direction are used to determine the coordinates of the lower left corner and upper right corner of the survey area boundary, and the two-dimensional spatial boundary range of the survey area is constructed based on this. The grid spacing parameters are set in combination with the density of survey line layout in the operation area and the desired spatial resolution. The side length of the grid unit is set to 0.1 meter to meet Figure 3 The terrain detail accuracy requirements shown in the figure are used to obtain the grid division parameter set; based on the above boundary and spacing parameters, the two-dimensional coordinate divisible mapping method is used to map the grid of each beam point. The coordinates are discretized, specifically the Subtract border The minimum value is divided by the grid side length and rounded down to get the column number. Similarly, the coordinates are used to obtain row numbers to form unique grid numbers, so that all beam points are assigned to equally spaced regular grid cells to construct a grid model dataset; the same grid cells covered by multiple survey lines in the grid model dataset are compared and processed, and the beam points belonging to each grid number are grouped according to the Ping number, and the point pairs belonging to different survey lines are paired one by one. The depth values ​​of each pair of beam points are extracted to perform difference calculation, and the results are used as the bathymetric difference of the grid points, and the bathymetric difference dataset is obtained by summarizing them; the bathymetric difference dataset is processed using basic statistical methods, and the arithmetic mean, absolute mean, maximum difference, minimum difference and root mean square error (RMS) of all differences are calculated in turn, where the mean reflects the offset direction of the system error, and the absolute mean and RMS are used to evaluate the overall filtering accuracy. The calculation methods are the square root of the average of all differences, the average of the absolute values ​​of the differences, and the average of the square values, and the results form an error statistical index table; all bathymetric differences are numerically grouped, and equal-width intervals are constructed with an interval width of 0.01 meters. Each difference is classified into the corresponding group according to the interval, and the number of samples in each interval is counted to form a difference frequency distribution table. Error histogram data is generated with frequency as the vertical axis and difference as the horizontal axis to construct error distribution layer data; finally, the error statistical index table and error distribution layer data are combined to analyze whether the error is symmetrical, whether the distribution is concentrated, and whether it is approximately normal, and to judge the consistency of the filtered point cloud data set between different survey lines, and then form an accuracy evaluation report to confirm that the systematic error and gross error have been effectively weakened during the implementation process, while retaining the true trend and detailed structure of the seabed topography.

[0099] It is particularly important that step S62 includes the following steps:

[0100] Step S621: parsing the grid division parameter set, extracting the measurement area boundary coordinate range and grid resolution parameters, constructing a two-dimensional grid index structure, and obtaining a grid index mapping template;

[0101] Step S622: performing coordinate discretization processing on each beam point in the adaptive trend filtering point cloud dataset, determining the corresponding grid row and column number based on its position, and obtaining the grid belonging dataset;

[0102] Step S623: Aggregate the grid-attributed data sets, organize all beam points under the same grid number into grid unit data, and construct an initial grid data unit set;

[0103] Step S624: Perform depth averaging processing on each grid cell in the initial grid data cell set, calculate the representative depth value, construct a regular grid model, and obtain a grid model data set.

[0104] In the embodiment of the present invention, first, the plane coordinates of all beam points are extracted based on the adaptive trend filtering point cloud data set that has completed the filtering correction. , among which The minimum and maximum values ​​of the coordinates are used as the boundary range of the survey area in the east-west direction. The minimum and maximum values ​​of the coordinates are used as the boundary range in the north-south direction to form a two-dimensional boundary coordinate frame of the survey area. Combined with the operating accuracy requirements of the multi-beam system and the goal of retaining terrain details, the grid resolution is set to 0.1 meters. The total number of rows and columns is calculated based on the boundary box and resolution, and a two-dimensional grid index structure is constructed in a column-first manner starting from the lower left corner, where the columns are numbered 0, 1, 2... from west to east, and the rows are numbered 0, 1, 2... from south to north to generate a grid index mapping template; for each beam point in the adaptive trend filtering point cloud dataset, the grid unit number to which it belongs is determined as follows: the grid unit number of the point is assigned to the grid cell. Coordinate subtraction Divide the minimum boundary by the grid side length and round it up to get the column number , Coordinate subtraction Divide the minimum boundary by the grid side length and round it up to get the row number ,use The grid number of the beam point is used as the grid number, and a grid affiliation dataset containing the grid number information of all beam points is finally obtained. Based on the grid affiliation dataset, beam points with the same grid number are clustered and organized into the same grid unit, forming a key-value mapping structure between the grid number and its subordinate beam point list, and output as the initial grid data unit set in the order of the grid number. The depth values ​​of all beam points in each grid unit in the initial grid data unit set are summed and divided by the number of points to calculate the representative depth value of the grid, ensuring that the depth result of each grid unit is actually representative and not dominated by local outliers. The calculation process uses a standard floating-point average calculation method and sets an outlier rejection rule (for example, if the number of points with depth values ​​exceeding the median ±0.5 meters in the grid accounts for more than 50%, the grid will not participate in the averaging process). Finally, a regular grid model is constructed with the average depth value corresponding to each grid number, and the grid model dataset is output.

[0105] Therefore, no matter from which point of view, the embodiments should be regarded as illustrative and non-restrictive, and the scope of the present invention is not limited by the above description. Therefore, it is intended that all changes that fall within the meaning and scope of the equivalent elements of the application documents are included in the present invention.

[0106] The foregoing description is intended only to provide specific embodiments of the present invention, which will enable those skilled in the art to understand and implement the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not intended to be limited to the embodiments shown herein, but is to be construed in the widest possible manner consistent with the principles and novel features disclosed herein.

Claims

1. A multi-beam point cloud filtering method for adaptive terrain trend, characterized by: Applied to a survey vessel equipped with a multi-beam bathymetric system, the method comprises the following steps: Step S1: Obtain the original multi-beam Ping data, extract the attributes of each beam point in each Ping, and obtain the beam point attribute dataset; Step S2: Perform multi-threshold condition screening on the beam point attribute dataset and remove abnormal beam points to obtain a preliminary filtered point cloud dataset. The multi-threshold condition screening includes depth screening, angle screening, width screening, slope screening, quality factor screening, and adjacency screening. Step S3: performing linear interpolation on the inter-Ping gaps in the removed area of ​​the preliminary filtered point cloud dataset to obtain a complete point cloud dataset; Step S4: extracting the central beam point data of the three adjacent survey lines based on the complete point cloud data set, and performing coordinate transformation and fitting modeling on the central beam point data to obtain a long-wave term function model; extracting the edge beam point data based on the complete point cloud data set, establishing an edge beam point error fitting function model, and obtaining a short-wave term function model; Step S5: performing depth correction on the edge beam points in the complete point cloud dataset based on the long-wave term function model and the short-wave term function model to obtain an adaptive trend filtered point cloud dataset; Step S6: Build a grid model based on the adaptive trend filtering point cloud data set, perform depth sounding error statistics on the grid points, and obtain the filtering accuracy evaluation result.

2. The multi-beam point cloud filtering method for adaptive terrain trend according to claim 1, characterized in that: Step S1 includes the following steps: Step S11: Using a survey vessel equipped with a multi-beam bathymetry system to perform underwater bathymetry scanning of the survey area to obtain raw multi-beam Ping data; Step S12: Decode and structurally analyze the original multi-beam Ping data, extract the Ping number, timestamp, heading, attitude, sound speed, each beam number and three-dimensional coordinate information, and obtain the Ping structure data set; Step S13: performing spatial inversion calculation on each beam point in the Ping structure data set based on the attitude and sound speed to obtain a beam point spatial coordinate data set; Step S14: performing attribute extraction processing on the beam point spatial coordinate data set to obtain a beam point attribute record set, where the attributes include the depth value, beam angle, beam number, quality factor, and Ping number of each beam point; Step S15: constructing a beam point attribute data set in a unified format based on the beam point attribute record set.

3. The multi-beam point cloud filtering method according to claim 2, characterized in that: Step S13 includes the following steps: Step S131: extracting the original bathymetric data of the current beam point from the Ping structure data set to obtain the bathymetric data of the beam point; Step S132: performing attitude compensation calculation on the beam point sounding data based on the attitude to obtain attitude compensated sounding data; Step S133: performing sound ray tracking correction on the attitude-compensated sounding data according to the sound velocity to obtain sound velocity-corrected sounding data; Step S134: Perform spatial inversion calculation based on the sound velocity corrected bathymetric data to obtain a beam point spatial coordinate data set.

4. The multi-beam point cloud filtering method according to claim 1, characterized in that: In step S2, the beam point attribute dataset is screened using multiple threshold conditions, and abnormal beam points are removed, including: Set the minimum depth and maximum depth thresholds, maximum beam angle threshold, maximum width threshold, maximum beam point slope and quality factor optional range based on the operating water depth range and regional terrain of the survey area; Determine the depth range of each beam point in the beam point attribute dataset, and remove beam points whose depth values ​​are not within the minimum depth and maximum depth thresholds to obtain a depth screening dataset; The beam angle range of each beam point in the beam point attribute dataset is judged, and beam points with beam angles greater than the maximum beam angle threshold are eliminated to obtain an angle screening dataset; For each beam point in the beam point attribute dataset, the horizontal distance from the beam point to the central beam point of the ping is calculated, and beam points whose horizontal distance exceeds the maximum width threshold are eliminated to obtain the width screening dataset. The depth slope of each beam point in the beam point attribute dataset and its adjacent beam points is calculated, and beam points with slopes greater than the maximum beam point slope are eliminated to obtain a slope screening dataset. Perform quality factor interval judgment on each beam point in the beam point attribute dataset, and eliminate beam points with quality factors below the lower limit or above the upper limit to obtain a quality screening dataset; Perform an adjacency check on each beam point in the beam point attribute dataset, and remove isolated beam points with missing adjacent ping numbers or left and right adjacent beam numbers to obtain an adjacency screening dataset. The depth filtering dataset, angle filtering dataset, width filtering dataset, slope filtering dataset, quality filtering dataset and adjacency filtering dataset are cross-fused, and the beam points that meet all the filtering conditions are retained to obtain the preliminary filtered point cloud dataset.

5. The multi-beam point cloud filtering method according to claim 1, characterized in that: Step S3 includes the following steps: Step S31: Based on the high-frequency sampling capability of the multi-beam bathymetry system and the Ping structure information, beam number matching is performed on the voids in the preliminary filtered point cloud dataset to obtain a set of beam point pairs to be interpolated; Step S32: performing spatial linear fitting calculation on each pair of beam points in the set of interpolation beam point pairs, obtaining interpolation coefficients, and obtaining an interpolation parameter set; Step S33: performing depth recovery on each interpolation target point in the interpolation parameter set, constructing a complete spatial coordinate record, and obtaining an interpolation point cloud dataset; Step S34: The preliminary filtered point cloud dataset and the interpolated point cloud dataset are jointly merged to construct a complete point cloud dataset.

6. The multi-beam point cloud filtering method according to claim 5, characterized in that: Step S32 includes the following steps: Step S321: extracting the three-dimensional spatial coordinate information and the corresponding Ping number for each pair of beam points in the interpolation beam point pair set, constructing a beam point pair spatial coordinate record, and obtaining a spatial point pair dataset; Step S322: extracting the horizontal coordinate position of the target interpolation point for each group of beam points in the spatial point pair data set, calculating its relative position scale coefficient, and obtaining an interpolation coefficient data set; Step S323: performing correlation matching on the interpolation coefficient dataset and the spatial point pair dataset, performing one-dimensional depth linear interpolation, and obtaining an interpolation depth value set; Step S324: Synchronously record the target Ping number and beam number information for the interpolation depth value set and the interpolation coefficient data set to construct an interpolation parameter set.

7. The multi-beam point cloud filtering method according to claim 1, characterized in that: In step S4, the central beam point data of three adjacent survey lines are extracted based on the complete point cloud data set, and coordinate transformation and fitting modeling are performed on the central beam point data, including: Perform index matching on the Ping data of three adjacent survey lines selected from the complete point cloud dataset, extract the central beam point in each Ping, and obtain the central beam point dataset; Perform coordinate system transformation on the central beam point data set to establish a local coordinate system with a certain point on the middle survey line as the origin, the starboard side as the y-axis, and the vertical downward direction as the z-axis, and obtain the transformed central beam point coordinate set; The least squares polynomial fitting is performed based on the central beam point coordinate set to construct a long-wave term function model of the terrain trend.

8. The multi-beam point cloud filtering method according to claim 1, characterized in that: In step S4, edge beam point data is extracted based on the complete point cloud data set, and an edge beam point error fitting function model is established, including: Extract edge beam points from three adjacent measurement lines in the complete point cloud dataset and construct an edge beam point dataset; Perform local coordinate transformation on the edge beam point data set to construct the edge beam point local coordinate set; Based on the local coordinate set of the edge beam point, the m-order polynomial fitting is performed using the least squares method to construct the shortwave term function model of the terrain error.

9. The multi-beam point cloud filtering method according to claim 1, characterized in that: Step S5 includes the following steps: Step S51: Perform local coordinate transformation on the edge beam points in the complete point cloud data set to obtain the local coordinate values ​​of the edge beam points. ; Step S52: Based on the shortwave term function model, the The error trend calculation is performed on the value to obtain the depth detail change of the edge beam point; Step S53: Based on the long-wave function model, the The trend surface estimation is performed on the value to obtain the trend prediction depth value; Step S54: Calculating a target corrected depth value based on the depth detail variation and the trend predicted depth value; Step S55: performing depth correction on each edge beam point using the target corrected depth value to construct a corrected edge beam point data set; Step S56: Merge the corrected edge beam point dataset with the unmodified central and non-edge beam points in the complete point cloud dataset, reconstruct the beam point structure of all Pings, and obtain an adaptive trend filtered point cloud dataset.

10. The multi-beam point cloud filtering method according to claim 1, wherein: Step S6 includes the following steps: Step S61: performing spatial range analysis on the adaptive trend filtering point cloud data set to determine the measurement area boundary range and resolution parameters, and obtain a grid division parameter set; Step S62: performing spatial projection and discrete mapping on the adaptive trend filtering point cloud dataset based on the grid division parameter set, constructing an equidistant regular grid model, and obtaining a grid model dataset; Step S63: Comparing and analyzing the results of different survey lines in the overlapping area of ​​the grid model dataset, calculating the sounding difference of the same grid point by locating the repeated coverage area of ​​the same grid point in multiple survey lines, and obtaining a sounding difference dataset; Step S64: performing statistical analysis on the sounding difference data set, calculating the difference mean, absolute average, maximum, minimum and sounding root mean square error, and obtaining an error statistical index table; Step S65: performing frequency distribution statistics on the depth difference data set, constructing error histogram data, and obtaining error distribution layer data; Step S66: Generate a filtering accuracy evaluation result based on the error statistical indicator table and the error histogram data set.

Images