Multi-beam installation deviation overall calibration method based on prior submarine topography

Abstract

The invention discloses a multi-beam installation deviation overall calibration method based on prior submarine topography, and belongs to the field of multi-beam data processing. The method mainly comprises the following steps: preprocessing multi-beam original data, performing point cloud registration to obtain matching point pairs, and resolving the installation deviation by total least squares. In submarine topography measurement operation of a multi-beam sounding system, after a transducer is disassembled and reassembled, installation deviation calibration needs to be carried out again, but an existing calibration method is complex in operation and low in efficiency. According to the method, the multi-beam installation deviation can be effectively calibrated, the operation process is simplified, the operation efficiency is improved, and the operation cost is reduced. A new solution is provided for the installation deviation calibration problem of the multi-beam sounding system, and the method has practical engineering application value.

Description

Technical Field

[0001] The present invention belongs to the field of multi-beam data processing and relates to a multi-beam installation deviation overall calibration method based on a priori seabed topography. The invention provides a new solution to the installation deviation calibration problem of a multi-beam bathymetric system and has practical engineering application value. Background Art

[0002] Seabed topography measurement is an important basic work to ensure navigation safety and the development of marine resources. The emergence of the multi-beam echo sounder system (MBES) has transformed the measurement method from traditional "point" and "line" measurement to strip-type and full-coverage "surface" measurement. It has the advantages of high precision and high efficiency and has now become the main instrument for seabed topography measurement. As a complex system engineering, the success of a measurement task depends not only on the measurement process, but also on the various preparations before the measurement operation begins and the data processing after the measurement is completed.

[0003] Installation deviation calibration is an important task before measurement. Installation deviation is caused by geometric offset between navigation module and bathymetric module, including misalignment of center points and inconsistency of three-axis pointing. When there is installation deviation in multi-beam system, the measured seabed terrain will be tilted, translated and rotated, and the terrain distortion parallel to the track will be generated when adjacent strips are spliced. Therefore, the installation deviation must be calibrated before measurement, otherwise the bathymetric data obtained will be inaccurate or even meaningless.

[0004] At present, the installation deviation calibration methods of multi-beam bathymetry systems can be roughly divided into two categories, namely, separate calibration and overall calibration. Separate calibration is a more intuitive method. It uses the seabed topography with certain characteristics to obtain the transducer installation offset parameters in a certain order. This method uses the changes in the seabed topography obtained by adjacent survey lines in the measurement area to calculate the installation deviation parameters in turn. After multiple iterative calculations, when the deviation meets the set limit, the iteration is stopped. This process is called Patch Test. Although this type of method is simple to operate, it has the following problems: First, the parameters are not obtained in the optimal sense, so their estimation accuracy is subject to certain limitations; second, since the parameter estimation is not completed at the same time, it is difficult to determine the correlation between the estimated values ​​of each parameter; third, independent targets or specific terrain are required to detect deviations, which places high requirements on the seabed topography of the calibration water area and the design of the measurement carrier route. In view of the limitations of the separate calibration method, scholars have carried out research on the overall calibration method of multi-beam installation deviation. The overall calibration method takes the measured value as the observation value, combines the prior seabed topography information in the survey area, and solves all the offset parameters according to a certain optimal criterion. The overall calibration method effectively solves the optimality and correlation problems of parameter estimation in separation calibration and can significantly improve operational efficiency. However, in practical applications, this method needs to rely on reliable seabed topography information to ensure the calibration effect, and it is difficult to promote it as a general method. Therefore, it is necessary to find suitable application scenarios to give full play to the advantages of the overall calibration method.

[0005] Based on an in-depth study of the impact of multi-beam installation deviations and calibration methods, the present invention aims to solve the problem of recalibrating installation deviations after reinstallation when the survey ship must dismantle the multi-beam transducer under certain environments (such as bad weather or complex sea conditions). This process is not only cumbersome to operate, but also involves the investment of a large amount of human and material resources. To solve this problem, the present invention proposes a method for calibrating installation deviations based on a priori terrain. This method takes into account the actual situation that both bathymetric and navigation data have measurement errors, constructs the installation deviation calibration problem as a weighted total least squares problem, and solves the problem, thereby achieving an overall solution to the installation deviation of the multi-beam system. Summary of the invention

[0006] The purpose of the present invention is to address the problem that the installation deviation of the multi-beam bathymetric system needs to be recalibrated after the transducer is disassembled and reinstalled during the seabed topography measurement operation, but the existing calibration methods are complicated to operate and inefficient. A new method for overall calibration of the installation deviation using a point cloud registration algorithm and a total least squares adjustment method is proposed. This method takes into account the fact that the heading will change the rigid transformation direction caused by the installation deviation when performing point cloud registration. At the same time, it takes into account the influence of navigation positioning, attitude, and sound speed errors when calculating the installation deviation, and is more complete in theory. It can not only effectively calibrate the multi-beam installation deviation, but also simplify the operation process, improve operation efficiency, and reduce operation costs. It provides a new solution to the installation deviation calibration problem of the multi-beam bathymetric system, and has practical engineering application value.

[0007] The technical solution of the present invention is as follows:

[0008] The overall calibration method of multi-beam installation deviation based on prior seabed topography mainly includes: preprocessing of multi-beam raw data, obtaining matching point pairs through point cloud registration, and solving installation deviation by overall least squares. The specific implementation steps are as follows:

[0009] Step 1: Multi-beam raw data preprocessing

[0010] The idea of ​​the present invention is to eliminate or weaken the influence of installation deviation in the measured terrain based on accurate prior terrain. The specific method is to obtain matching point pairs through point cloud registration, and then solve the installation deviation parameters based on the information of the matching point pairs. Therefore, it is first necessary to obtain the measured terrain point cloud data without installation deviation calibration and the accurate prior terrain point cloud data through depth measurement point position reduction. This process can be carried out in Caris software, and the specific content will not be repeated.

[0011] Step 2: Point cloud registration to obtain matching point pairs

[0012] After obtaining the terrain point cloud data, the prior terrain point cloud is used as the target point cloud and the measured terrain point cloud is used as the measured point cloud. First, the measured point cloud is grouped according to the heading, and then each group of point cloud data is registered with the target point cloud, and the registration results of each group are integrated to form the final registration result. The following is a detailed introduction to this process:

[0013] (1) Point cloud data grouping

[0014] The measured point cloud data is grouped at fixed angle intervals, and each fixed angle interval is divided into an interval starting from the minimum heading. In each heading interval, the measured point cloud data is regarded as an independent data set, and its installation deviation direction is set to be consistent, so that accurate registration operations can be implemented.

[0015] Preferably, the fixed angle is set to 5°.

[0016] The installation deviation is specifically manifested as a three-dimensional rigid transformation of the transducer coordinate system relative to the survey ship coordinate system. When the survey ship's heading changes, the x-axis direction of the survey ship coordinate system will also deflect, resulting in a corresponding change in the direction of the rigid transformation. For example, roll deviation mainly causes the tilt of the terrain in the direction perpendicular to the track. When the heading changes, the direction of the terrain tilt will also change. At this time, if the measurement point clouds of different headings are directly mixed together and aligned with the target point cloud, the alignment results will be inaccurate due to the differences in the rigid transformation directions of the point clouds of each heading. Therefore, it is necessary to group the measurement point clouds according to the heading so that the measurement point cloud under each heading has the same relative structure as the target point cloud to ensure the effectiveness of the alignment process.

[0017] The present invention attempts to group the headings using multiple grouping intervals based on the heading angle. Studies have found that if the grouping interval is set too small, the amount of point cloud data in each group will be insufficient, the distribution of points in the point cloud will be uneven, it will be unfavorable for extracting point cloud features, and effective registration will not be possible, while increasing the probability of matching conflicts; conversely, if the grouping interval is too large, it will lead to significant differences in the rigid transformation direction, thereby reducing the accuracy of the registration results. After comparison, it was finally decided to group the measured point cloud data at intervals of 5°, and to divide an interval every 5° starting from the minimum heading. In each heading interval, the measured point cloud data is regarded as an independent data set, assuming that the installation deviation direction is consistent, and accurate registration operations can be implemented.

[0018] (2) Use the RANSAC algorithm and FPFH features for point cloud coarse registration

[0019] Affected by many factors such as complex marine environment, ship noise, and unreasonable manual operation, there are inevitably a large number of outliers in the multi-beam measurement point cloud data. These outliers may cause the ICP algorithm to fall into a local optimal solution, seriously affecting the point cloud registration results. The FPFH (fast point feature histograms) feature can effectively capture the local geometric information of the point cloud, thereby providing more representative features for point cloud matching. The RANSAC algorithm is an iterative method for estimating the parameters of a mathematical model from a set of observation data containing outliers. Its basic idea is to randomly select the minimum number of data points, use them to estimate the model parameters, and then determine how many data points are suitable for the model. Using the FPFH feature as a reference and the RANSAC algorithm to estimate the initial transformation matrix of the point cloud can not only find the correspondence of the point cloud on a larger scale, but also exclude most of the outliers, thereby improving the robustness and convergence of the subsequent ICP. Therefore, it is necessary to use the FPFH feature for RANSAC point cloud coarse registration before ICP fine registration.

[0020] The first step is to use the k-neighborhood method to calculate the normal vector of each point in the point cloud:

[0021] The target point cloud and the measured point cloud are processed as follows: for each point I, select its k nearest neighboring points, denoted as N(I) = {I 1 ,I 2 ,…,I k}. Calculate the covariance matrix C based on the selected neighborhood points I , for C I Perform eigenvalue decomposition to obtain the eigenvalue and its corresponding eigenvector ξ 1 , 2 , 3 The normal vector is the eigenvector corresponding to the minimum eigenvalue. The normal vector is the basis for FPFH feature calculation.

[0022] The second step is to calculate the FPFH feature histogram of each point and its neighborhood points. It only needs to calculate the feature elements between point I and its neighborhood points:

[0023] Let the normal vector of point I be n I , each neighboring point I of point I a The normal vector is n Ia , neighborhood point I a The distance vector from point I is v Ia , a=1,2,…,k. Then n I and n Ia 、n I and v Ia The angles are

[0024]

[0025] According to the angle calculation result of formula (1), fill the FPFH feature histogram and perform normalization. The contribution FPFH of each neighborhood point can be expressed as

[0026] FPFH(I a )=(θ nn ,θ nv ) (2)

[0027] The third step is to use the RANSAC algorithm to obtain the optimal transformation matrix:

[0028] Set a threshold, and calculate the difference in the FPFH values ​​of the points in the two groups of point clouds based on the FPFH feature histogram obtained in the second step. If it is less than the threshold, the two points are taken as matching point pairs. Perform bidirectional matching in this way to obtain a set of matching point pairs: randomly select a set of matching point pairs from the matching point pair set, calculate the rigid transformation matrix from the measured point cloud to the target point cloud based on the selected matching point pairs, and apply the rigid transformation matrix to the measured point cloud to obtain the transformed measured point cloud. For all matching point pairs, use the transformed measured point cloud to compare with the target point cloud, and calculate the number of inliers within a given threshold. Inliers refer to points in the transformed measured point cloud whose distance to the corresponding point in the target point cloud is less than the set threshold. Repeat the process of random selection, measurement point cloud transformation, and calculation of the number of inliers until the iteration ends. Select the rigid transformation matrix with the largest number of inliers in all iterations as the coarse registration result and pass it to the ICP algorithm as the initial rigid transformation matrix.

[0029] (3) ICP algorithm realizes precise point cloud registration

[0030] The RANSAC algorithm is used for coarse registration to obtain the initial transformation matrix, which provides a reasonable initial value for the subsequent ICP algorithm. The ICP algorithm is a classic algorithm for point cloud registration. It continuously minimizes the distance between the measured point cloud and the target point cloud to find the best rigid transformation (translation and rotation) so that the measured point cloud is aligned with the target point cloud as much as possible. On the basis of coarse registration, the ICP algorithm can be used to eliminate subtle differences between point clouds and ultimately achieve the purpose of precise alignment. The ICP algorithm mainly includes three steps: estimating corresponding points, minimizing errors, and iterative optimization. All the measured point clouds mentioned in this step are the latest measured point clouds after transformation.

[0031] The first step is to estimate the corresponding points. For each point in the measured point cloud, the nearest neighbor method is used to solve the point closest to the target point cloud as its corresponding point. However, due to the large number of points, the time complexity of this operation is very high. In order to improve the efficiency of the algorithm, when solving the corresponding point of a point in the measured point cloud, it is not necessary to calculate the distance to each point in the target point cloud. A threshold can be set. When the distance is less than the threshold, it is used as a corresponding point. The threshold has a decisive influence on the result of point cloud registration. If the threshold is selected too high, the root mean square error of the inliers in the matching result will increase, thereby causing mismatching; conversely, if the threshold is set too low, the matching degree will be too small, and it will not be possible to obtain enough matching point pairs for calculation. Since the matching degree will be affected by the scale of the target point cloud, it is not suitable as a unified standard for selecting the threshold. Therefore, the threshold is determined based on the root mean square error (RMSE) of the inliers. Under the premise of ensuring that the RMSE of the inliers does not exceed a certain value (set according to experience, preferably 0.2), a larger threshold is preferred, but its upper limit shall not exceed 1. Construct two sets of corresponding point sets based on the corresponding points obtained:

[0032]

[0033] where p w (w=1,2,…,c) is the coordinate of the point in the measured point cloud, q w (w=1,2,…,c) is p w The coordinates of the point in the corresponding target point cloud.

[0034] The second step is to minimize the error. After obtaining the corresponding point set, the ICP algorithm is used to solve the rigid transformation R ICP (rotation matrix) and t ICP (translation vector), R ICP (rotation matrix) and t ICP (translation vector) is applied to the currently updated measured point cloud so that the difference between the transformed measured point cloud and the target point cloud is as small as possible. Let the optimized objective function be the sum of square errors of the distances between corresponding points after transformation, that is,

[0035]

[0036] After processing formula (4), the final objective function form is as follows:

[0037]

[0038] in, are the centroids of the two groups of point clouds, q′ w , p′ w is the coordinate of the point after de-centroiding. Analyzing formula (5), we find that for the latter Regardless of R ICP Whatever the value, there is always a corresponding t ICP Make this term equal to 0. Therefore, we only need to solve R that minimizes the value of the previous term. ICP , and then let the latter term equal to 0 to solve for t ICP , the objective function J can be minimized. The optimal rotation matrix R * for

[0039]

[0040] make Performing SVD decomposition on K yields K = U K Σ K V K T , where Σ K is a diagonal matrix of singular values.

[0041]

[0042] From the properties of SVD decomposition, we know that V K , U K 、V K T R ICP U K are all orthogonal matrices. Since the element value of the orthogonal matrix is ​​less than or equal to 1, when V K T R ICP U K When is the unit matrix, formula (7) takes the maximum value.

[0043]

[0044] Among them, t * is the optimal translation vector.

[0045] The third step is iterative optimization. * and t * Apply to the current updated measured point cloud, update the measured point cloud, and calculate the registration error. If the current error is less than the predetermined tolerance range relative to the previous round, the registration is considered to have converged. If not, repeat steps 1 to 3 for a new round of iterations until the convergence condition is met or the maximum number of iterations is reached. Each round of iteration will gradually bring the measured point cloud closer to the target point cloud, ultimately achieving high-precision point cloud registration.

[0046] After completing the point cloud registration, matching point pairs can be obtained according to the corresponding point indexes. The matching point pairs are substituted into the model of the present invention to solve the installation deviation parameters and complete the installation deviation calibration.

[0047] Step 3: Overall least squares solution for installation deviation

[0048] In traditional least squares, an implicit premise is that the observed value contains errors, while the coefficients are considered to be error-free. For the overall calibration problem of installation deviation, it is no problem to consider the water depth as a measured value to contain errors, but the coefficients are composed of the position and attitude of the survey ship, which inevitably have errors (especially attitude errors). Obviously, this problem has obvious error-in-variable (EIV) characteristics, that is, there are errors in the independent variables. Taking this feature into account, the corresponding solution should be the total least squares solution. The following is a detailed description of the installation deviation solution method based on the total least squares model.

[0049] (1) Coordinate system establishment and parameter definition

[0050] In order to facilitate the description of various variables and their mutual relations, the present invention introduces the local horizontal coordinate system, the ship surveying coordinate system and the depth measurement center coordinate system. The origin b of the ship surveying coordinate system is the center of the surveying ship, the x-axis is parallel to the keel line direction of the surveying ship, pointing to the bow is positive, the y-axis points to the starboard is positive, and the z-axis is perpendicular to the bxy plane and constitutes a right-handed coordinate system. The origin l of the local horizontal coordinate system is located on the depth level surface, the X-axis and the Y-axis coincide with the X-axis and the Y-axis of the local UTM coordinate system, and the Z-axis is perpendicular to the lXY plane and constitutes a right-handed coordinate system. The origin m of the depth measurement center coordinate system is the geometric center of the transducer, the x'axis is the center line of the transmitting transducer, pointing to the bow is positive, the y'axis is the center line of the receiving transducer, pointing to the starboard is positive, and the z'axis is perpendicular to the mx'y'plane and constitutes a right-handed coordinate system.

[0051] The offset between the geometric center of the transducer and the center of the survey ship is defined as the translation offset parameter (TOC), and the position vector Indicated by, where the subscript represents the vector from point b to point m, and the superscript represents the coordinate of the vector in the b system; the inconsistency between the three-axis directions of the survey ship coordinate system and the sounding center coordinate system is defined as the rotation offset parameter (ROC), and the rotation matrix is ​​used The upper and lower subscripts indicate that the rotation matrix is ​​used to transform the coordinates of a vector in the m system into the coordinates of the vector in the b system. These two examples illustrate the meaning of the upper and lower subscripts of the position vector r and the rotation matrix R. The same applies when the upper and lower subscripts change, and will not be repeated below.

[0052] (2) Establishment of installation deviation calculation model

[0053] The idea of ​​this model is to use the seabed target point with known position and depth as the benchmark, and use the idea of ​​adjustment to calculate the installation deviation parameters that meet the highest accuracy conditions. Given that the current GNSS second pulse (1pps) synchronization technology has been widely used in hydrographic surveying practice, the delay factor will not be considered in this model. Assume that the position vector of the target point t is At a certain moment, the beam footprint covers the target point, and the position vector of the measurement carrier is The carrier's posture matrix is TOC is ROC is The position of the target point relative to the center of the transducer is Then there is the following functional relationship:

[0054]

[0055] Reorganize equation (9) and according to the properties of the rotation matrix, we get

[0056]

[0057] in As the location of the target point is known; and It needs to be calculated based on the measurement data of the MBES navigation module and the water level observation data, and its value contains measurement errors; It is obtained by using the sound velocity profile to track the sound line. The calculation result contains the sound velocity error. The least square method sets equation (10) to contain an error only on the left side of the equal sign. Obviously, this is not in line with reality. Therefore, we reorganize equation (10) to obtain

[0058]

[0059] Defining Observables and

[0060]

[0061] in, for The observation quantity, i = 1, 2, ..., n represents the i-th target point, e i and ε i is the measurement error. Substituting equation (12) into equation (11), we have

[0062]

[0063] Let the cost function be the weighted sum of squares of measurement errors, that is

[0064]

[0065] Where P i and Q i e i and ε i The problem now becomes a cost function minimization problem that satisfies the constraints of equation (13), and the Lagrange multiplier λ is introduced. i , construct the following Lagrange function

[0066]

[0067] According to the first-order necessary condition for minimizing φ, we have

[0068]

[0069] Combining equation (16) and equation (13), we get

[0070]

[0071] Where W i =P i +Q i, Substituting formula (17) into formula (14), we have

[0072]

[0073] according to The first-order necessary condition for minimization is

[0074]

[0075] Solved

[0076]

[0077] Substituting equation (20) into equation (18), we get

[0078]

[0079] in

[0080]

[0081] Expand (21) and separate Items, there are

[0082]

[0083] Since the first two terms on the right side of equation (23) are It doesn't matter, so Minimizing is equivalent to maximizing the following function

[0084]

[0085] make Perform singular value decomposition (SVD) on H, and we have

[0086]

[0087] where σ 1 , σ 2 , σ 3 is the singular value of matrix H, U and V are orthogonal matrices. Substituting equation (25) into equation (24), we have

[0088]

[0089] To maximize g, while ensuring the rotation matrix The determinant is equal to 1, so

[0090]

[0091] Solve ROC according to formula (27), and substitute formula (27) into formula (20) to obtain TOC. At this point, the calculation of transducer installation deviation is completed.

[0092] Beneficial effects of the method of the present invention:

[0093] In the seabed topography measurement operation of the multi-beam bathymetry system, the transducer needs to be recalibrated for installation deviation after disassembly and reinstallation, but the existing calibration method is complicated and inefficient. The method of the present invention can not only effectively calibrate the multi-beam installation deviation, but also simplify the operation process, improve the operation efficiency and reduce the operation cost. It provides a new solution to the installation deviation calibration problem of the multi-beam bathymetry system and has practical engineering application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 The figure is a flow chart of the method of the present invention.

[0095] Figure 2 Schematic diagram of the coordinate system and target point position.

[0096] Figure 3 The color map of the seabed topography in the test area, where (a) is the prior seabed topography measured in the first voyage, and (b) is the seabed topography with installation deviation measured in the second voyage.

[0097] Figure 4 A three-dimensional seafloor topography map of the selected area.

[0098] Figure 5 A map of the seafloor topography measured along the calibration line.

[0099] Figure 6 Schematic diagram of the position relationship between the measured point cloud and the target point cloud before and after point cloud registration, where (a) is a schematic diagram of the position relationship between the measured point cloud and the target point cloud before point cloud registration, and (b) is a schematic diagram of the position relationship between the measured point cloud and the target point cloud after point cloud registration.

[0100] Figure 7 The figures are statistical diagrams of the water depth discrepancies of the measurement data of two voyages on the same survey line before and after the installation deviation calibration (i.e., the external conformity) and the water depth discrepancies of the overlapping areas of adjacent survey lines of the second voyage before and after the calibration (i.e., the internal conformity), wherein (a) is the external conformity of the water depth data before the installation deviation calibration, (b) is the external conformity of the water depth data after the installation deviation calibration by the Patch Test method, (c) is the external conformity of the water depth data after the installation deviation calibration by the method of the present invention, (d) is the internal conformity of the water depth data before the installation deviation calibration, (e) is the internal conformity of the water depth data after the installation deviation calibration by the Patch Test method, and (f) is the internal conformity of the water depth data after the installation deviation calibration by the method of the present invention. DETAILED DESCRIPTION

[0101] In order to make the problems solved, the methods adopted and the effects achieved by the model of the present invention clearer, the present invention is further described in detail below in conjunction with the accompanying drawings and experiments. It is understood that the specific experiments described herein are only used to explain the present invention, rather than to limit the present invention. In addition, it should be noted that, for the convenience of description, only the parts related to the present invention are shown in the accompanying drawings, rather than all the contents.

[0102] 1. Experimental introduction

[0103] The test area was selected in a certain sea area, where the average water depth is about 40m. The Seabat T-20P multi-beam bathymetric system was used to measure the seabed topography. The transducer operating frequency was 400kHz, the number of single Ping beams was 256, the beam footprint angle was 1°×1°, the maximum equidistant beam opening angle was 140°, the maximum equiangular opening angle was 160°, the maximum range was 300m, and the highest Ping rate reached 50ping / s. Two voyages were designed to measure the same survey line. Installation deviation calibration was performed in one voyage; in the second voyage, the transducer was removed and reinstalled, and installation deviation calibration was not performed. The bathymetric data of the two voyages were corrected for navigation attitude, sound ray tracking, and water level (the bathymetric data of the first voyage was also calibrated for installation deviation), and outliers and edge beam data with low accuracy were eliminated to obtain the coordinates of the bathymetric points in the local horizontal coordinate system. Figure 3 Color maps of the seafloor topography obtained from the two measurements.

[0104] contrast Figure 3 In (a) and (b), it can be observed that there are obvious differences in the seabed topography measured in the two voyages, indicating that the installation deviation has a significant impact on the bathymetric data. In this experiment, the seabed topography measured in the first voyage is used as the prior topography to calibrate the installation deviation of the second voyage. Figure 3 It can be found that the topographic features of area 1 are relatively complex and the terrain is undulating. This terrain condition is conducive to accurately identifying the installation deviation. To verify the effectiveness of the method of the present invention, the seabed topography of area 1 measured in the first voyage is used as the priori terrain to calibrate the installation deviation of the second voyage. To evaluate the robustness of the algorithm of the present invention, a certain number of outliers are artificially introduced into the bathymetric data. Finally, the following is obtained: Figure 4 Three-dimensional seafloor topography of the area shown.

[0105] To verify the correctness of the method of the present invention, two calibration lines J1 and J2 were designed in an area containing flat terrain and characteristic terrain. In the second voyage, two measurements were made along one of the calibration lines at the same speed and opposite heading, and one measurement was made along each of the two calibration lines at the same speed and heading. The seabed terrain swept by the survey ship when measuring along the calibration line is shown in Fig. Figure 5After the measurement is completed, the installation deviation is calculated using the Patch Test method and the calculation results are used as a reference. Since the multi-beam measurement system uses a positioning receiver that can transmit 1pps signals and 1pps data synchronization is performed before the measurement, there is no need for delay calibration. The bathymetric data of the two survey lines are also used for statistical analysis of water depth discrepancies in the overlapping area.

[0106] 2. Test content

[0107] The bathymetric data with horizontal coordinates in area 1 were extracted from the two voyages respectively. The navigation, attitude and surface sound speed information of the second voyage were extracted using the Command_line tool of the Caris software. These auxiliary information were then matched with the beam through timestamp synchronization and linear interpolation methods to provide the necessary input parameters for model calculation. Finally, the bathymetric data set of the second voyage and the bathymetric data set of the first voyage were formed.

[0108] The bathymetric data set of the first voyage was directly used as the target point cloud in the point cloud registration process. The bathymetric data set of the second voyage was grouped according to the heading. The minimum heading recorded in the data was 155.6°, and the data set was divided into two groups: 155.6°-160.6° and 160.6°-165.6°. The first group had a total of 13,925 points, and the second group had a total of 29,655 points.

[0109] After completing the grouping of the measurement point cloud, the point cloud registration algorithm of the present invention is used to match the corresponding points between each group of measurement point cloud and the target point cloud. The first group of measurement point cloud and the target point cloud are registered, and the matching threshold is set to 0.55m, and 13865 pairs of matching point pairs are obtained, with a matching degree of 99.57% and an inlier root mean square error of 0.154; the second group of measurement point cloud and the target point cloud are registered, and the matching threshold is set to 0.55m, and 29392 pairs of matching point pairs are obtained, with a matching degree of 99.11% and an inlier root mean square error of 0.191. The spatial relationship between the measurement point cloud and the target point cloud before and after registration is shown in the figure. Figure 6 As shown. Figure 6 In (a) and (b), it can be observed that the measured point cloud and the target point cloud after registration have a high degree of overlap, and the method proposed in the present invention achieves a good registration effect.

[0110] The two sets of matching point pairs are combined to obtain a matching point data set, which is substituted into the overall least squares installation deviation calculation model to calculate the installation deviation parameter estimation. The calculation results are then compared with the calculation results of the Patch Test, and the results are shown in Table 1. As can be seen from Table 1, the calculation results of the method of the present invention are basically consistent with those of the Patch Test, thereby verifying the correctness of the method of the present invention.

[0111] Table 1 is a diagram showing the analysis of the calculation results of the installation deviation of the method of the present invention and the Patch Test

[0112]

[0113] In the practice of water depth measurement, the comparison of the sounding values ​​of the overlapping area of ​​adjacent survey lines in the same measurement is called internal conformity, and the comparison of the sounding values ​​of the same area in different measurements is called external conformity. The internal conformity accuracy and external conformity accuracy are widely used to evaluate the quality of water depth measurement data. In order to further verify the effectiveness of the method of the present invention under complex terrain conditions, the calculation results of the method of the present invention and the Patch Test were respectively applied to the installation deviation calibration of the second voyage, and the water depth discrepancies of the measurement data of the two voyages of the same survey line before and after the installation deviation calibration (i.e., external conformity) and the water depth discrepancies of the overlapping area of ​​adjacent survey lines in the second voyage before and after the calibration (i.e., internal conformity) were statistically analyzed. The statistical results are shown in Figure 7 .

[0114] observe Figure 7 From (a), (b) and (c) in Figure 7 Compared with (a), the depth discrepancy bar graph in (c) has become higher and narrower, the standard deviation of the discrepancy value has dropped from 2.143 before calibration to 0.031, and the mean absolute value of the discrepancy value has dropped from 1.792 before calibration to 0.005, and the discrepancy value is closer to 0 as a whole. This result shows that the difference between the bathymetric data of the second voyage and the bathymetric data of the first voyage after installation deviation calibration has been significantly reduced, thus verifying the effectiveness of this research method under complex terrain conditions. Figure 7 From the analysis of (b) in the figure, we can find Figure 7 The distribution law of the non-conforming values ​​in (c) and (b) is basically the same, indicating that the two methods have equivalent calibration effects on installation deviation. However, compared with the Patch Test method, the method of the present invention has higher operating efficiency and lower operating cost, and has significant advantages in practical applications.

[0115] observe Figure 7 (d), (e) and (f) in the figure can reach the same conclusion. Figure 7 In (d) and (f), the standard deviation of the discrepancy values ​​in the overlapping area of ​​adjacent survey lines decreased from 4.919 to 0.143, and the mean absolute value of the discrepancy values ​​decreased from 0.91 to 0.039, and the discrepancy values ​​as a whole were closer to 0. This shows that after using the method of the present invention to calibrate the installation deviation, the problem of abnormal splicing of the bathymetric data of adjacent survey lines has been greatly improved. Figure 7 From the analysis of (e) in the figure, we can find Figure 7 The distribution rules of the non-conforming values ​​in (f) and (e) are basically the same, and the method of the present invention can achieve the same effect as the Patch Test method.

[0116] Finally, it should be noted that the above experiments are only used to illustrate the method scheme of the present invention, rather than to limit it.

Claims

1. A multi-beam installation deviation overall calibration method based on a priori seabed topography, characterized in that: The steps include: Step 1: Multi-beam raw data preprocessing Obtain the measured terrain point cloud data and accurate prior terrain point cloud data without installation deviation calibration by calculating the depth measurement point position; Step 2: Point cloud registration to obtain matching point pairs After obtaining the terrain point cloud data, the prior terrain point cloud is used as the target point cloud and the measured terrain point cloud is used as the measured point cloud. First, the measured point cloud is grouped according to the heading, and then each group of point cloud data is registered with the target point cloud, and the registration results of each group are integrated to form the final registration result. Step 3: Solve the installation deviation using the overall least squares method.

2. According to the multi-beam installation deviation overall calibration method based on a priori seabed topography according to claim 1, the step 2 is specifically as follows: (1) Point cloud data grouping The measured point cloud data are grouped at fixed angle intervals. Starting from the minimum heading, each fixed angle interval is divided into an interval. In each heading interval, the measured point cloud data are regarded as an independent data set, and the installation deviation direction is set to be consistent, so that accurate registration operations can be implemented. (2) Use the RANSAC algorithm and FPFH features for point cloud coarse registration The first step is to use the k-neighborhood method to calculate the normal vector of each point in the point cloud: The target point cloud and the measured point cloud are processed as follows: for each point I, select its k nearest neighboring points, denoted as N(I)={I1,I2,,I k }; Calculate the covariance matrix C based on the selected neighborhood points I , for C I Perform eigenvalue decomposition to obtain the eigenvalue and its corresponding eigenvectors ξ1, ξ2, ξ3; the normal vector is the eigenvector corresponding to the minimum eigenvalue; The second step is to calculate the FPFH feature histogram of each point and its neighborhood points: Set point The normal vector of I is n I , each neighboring point I of point I a The normal vector is n Ia , neighborhood point I a The distance vector from point I is v Ia , a=1,2,,k; then n I and n Ia 、n I and v Ia The angles are According to the angle calculation result of formula (1), the FPFH feature histogram is filled and normalized; the contribution FPFH of each neighborhood point can be expressed as FPFH(I a )=(θ nn ,i nv )(2) The third step is to use the RANSAC algorithm to obtain the optimal transformation matrix: Set a threshold, and calculate the difference in FPFH values ​​between the points in the two groups of point clouds based on the FPFH feature histogram obtained in the second step. If it is less than the threshold, the two points are taken as matching point pairs, and bidirectional matching is performed in this way to obtain a set of matching point pairs: randomly select a set of matching point pairs from the matching point pair set, calculate the rigid transformation matrix from the measured point cloud to the target point cloud based on the selected matching point pairs, and apply the rigid transformation matrix to the measured point cloud to obtain the transformed measured point cloud; for all matching point pairs, use the transformed measured point cloud to compare with the target point cloud, and calculate the number of inliers within a given threshold; inliers refer to points in the transformed measured point cloud whose distance to the corresponding point in the target point cloud is less than the set threshold; repeat the process of random selection, measurement point cloud transformation and inlier number calculation until the iteration ends; select the rigid transformation matrix with the largest number of inliers in all iterations as the coarse registration result, and pass it to the ICP algorithm as the initial rigid transformation matrix; (3) ICP algorithm realizes precise point cloud registration The ICP algorithm finds the best rigid transformation by continuously minimizing the distance between the measured point cloud and the target point cloud, so that the measured point cloud is aligned with the target point cloud as much as possible; all the measured point clouds mentioned in this step are the latest measured point clouds after transformation; The first step is to estimate the corresponding points. For each point in the measured point cloud, the nearest point in the target point cloud is solved by the nearest neighbor method as its corresponding point. When solving the corresponding point of a point in the measured point cloud, it is not necessary to calculate the distance to each point in the target point cloud. A threshold can be set. When the distance is less than the threshold, it is regarded as the corresponding point. The threshold is determined based on the root mean square error of the inliers. Under the premise of ensuring that the RMSE of the inliers does not exceed a certain value, a larger threshold is selected, but its upper limit shall not exceed 1; two sets of corresponding point sets are constructed based on the corresponding points obtained: where p w (w=1,2,,c) is the coordinate of the point in the measured point cloud, q w (w=1,2,,c) is p w The coordinates of the points in the corresponding target point cloud; The second step is to minimize the error; after obtaining the corresponding point set, the ICP algorithm is used to solve the rigid transformation R ICP (rotation matrix) and t ICP (translation vector), R ICP (rotation matrix) and t ICP (translation vector) is applied to the currently updated measured point cloud so that the difference between the transformed measured point cloud and the target point cloud is as small as possible; let the optimized objective function be the sum of square errors of the distances between corresponding points after transformation, that is, After processing formula (4), the final objective function form is as follows: in, are the centroids of the two groups of point clouds, q′ w , p′ w is the coordinate of the point after de-centroiding; analyzing formula (5), it is found that for the latter Regardless of R ICP Whatever the value, there is always a corresponding t ICP Make this term equal to 0; therefore, we only need to solve R that minimizes the value of the previous term. ICP , and then let the latter term equal to 0 to solve for t ICP , the objective function J can be minimized; the optimal rotation matrix R * for make Performing SVD decomposition on K yields K = U K Σ K V K T , where Σ K is a diagonal matrix composed of singular values; then From the properties of SVD decomposition, we know that V K , U K 、V K T R ICP U K are all orthogonal matrices; since the element value of the orthogonal matrix is ​​less than or equal to 1, when V K T R ICP U K When is the unit matrix, formula (7) takes the maximum value; thus we can get Among them, t * is the optimal translation vector; The third step is iterative optimization; the calculated R * and t * Apply to the current updated measured point cloud, update the measured point cloud, and calculate the registration error; if the current error is less than the predetermined tolerance range relative to the previous round, the registration is considered to have converged; if not, repeat steps 1 to 3 and perform a new round of iterations until the convergence condition is met or the maximum number of iterations is reached; each round of iteration will gradually bring the measured point cloud closer to the target point cloud, and ultimately achieve high-precision point cloud registration; After completing the point cloud registration, matching point pairs can be obtained according to the corresponding point indexes. The matching point pairs are substituted into the model of the present invention to solve the installation deviation parameters and complete the installation deviation calibration.

3. According to the multi-beam installation deviation overall calibration method based on a priori seabed topography according to claim 1, the step three is specifically as follows: (1) Coordinate system establishment and parameter definition The local horizontal coordinate system, the ship surveying coordinate system and the bathymetric center coordinate system are introduced; the origin b of the ship surveying coordinate system is the center of the ship surveying, the x-axis is parallel to the keel line of the surveying ship, pointing to the bow is positive, the y-axis points to the starboard is positive, and the z-axis is perpendicular to the bxy plane and forms a right-handed coordinate system; the origin l of the local horizontal coordinate system is located on the depth level surface, the X-axis and the Y-axis coincide with the X-axis and the Y-axis of the local UTM coordinate system, and the Z-axis is perpendicular to the lXY plane and forms a right-handed coordinate system; the origin m of the bathymetric center coordinate system is the geometric center of the transducer, the x'axis is the center line of the transmitting transducer, pointing to the bow is positive, the y'axis is the center line of the receiving transducer, pointing to the starboard is positive, and the z'axis is perpendicular to the mx'y'plane and forms a right-handed coordinate system; The offset between the geometric center of the transducer and the center of the survey ship is defined as the translation offset parameter (TOC), and the position vector Indicated by, where the subscript represents the vector from point b to point m, and the superscript represents the coordinate of the vector in the b system; the inconsistency between the three-axis directions of the survey ship coordinate system and the sounding center coordinate system is defined as the rotation offset parameter (ROC), and the rotation matrix is ​​used The upper and lower subscripts indicate that the rotation matrix is ​​used to transform the coordinates of a vector in the m system into the coordinates of the vector in the b system. These two examples illustrate the meaning of the upper and lower subscripts of the position vector r and the rotation matrix R. The same applies when the upper and lower subscripts change, and will not be repeated below. (2) Establishment of installation deviation calculation model Assume the position vector of the target point t is At a certain moment, the beam footprint covers the target point, and the position vector of the measurement carrier is The carrier's posture matrix is TOC is ROC is The position of the target point relative to the center of the transducer is Then there is the following functional relationship: Reorganize equation (9) and according to the properties of the rotation matrix, we get in As the location of the target point is known; and It needs to be calculated based on the measurement data of the MBES navigation module and the water level observation data, and its value contains measurement errors; It is obtained by using the sound velocity profile to track the sound line. The calculation result contains the sound velocity error. The least square method sets equation (10) with only the error on the left side of the equal sign. Obviously, this is not in line with reality. Therefore, we reorganize equation (10) to obtain Defining Observables and in, for The observation quantity, i = 1, 2, ..., n represents the i-th target point, e i and ε i is the measurement error; Substituting formula (12) into formula (11), we have Let the cost function be the weighted sum of squares of measurement errors, that is Where P i and Q i e i and ε i The weight matrix; At this time, the problem becomes a cost function minimization problem that satisfies the constraints of formula (13), and the Lagrange multiplier λ is introduced i , construct the following Lagrange function According to the first-order necessary condition for minimizing φ, we have Combining equation (16) and equation (13), we get Where W i =P i +Q i , Substituting formula (17) into formula (14), we have according to The first-order necessary conditions for minimization are Solved Substituting equation (20) into equation (18), we get in Expand (21) and separate Items, there are Since the first two terms on the right side of equation (23) are It doesn't matter, so Minimizing is equivalent to maximizing the following function make Perform singular value decomposition (SVD) on H, and we have Where σ1, σ2, σ3 are the singular values ​​of the matrix H, and U and V are orthogonal matrices. Substituting equation (25) into equation (24), we have To maximize g, while ensuring that the rotation matrix The determinant is equal to 1, so According to formula (27), ROC is solved and formula (27) is substituted into formula (20) to obtain TOC. At this point, the calculation of transducer installation deviation is completed.

Images