Multi-beam coordinate auxiliary side-scan image ground feature accurate positioning method and system

Abstract

The application relates to the technical field of marine surveying and mapping, and discloses a multi-beam coordinate auxiliary side-scan image ground feature accurate positioning method and system. The method comprises the following steps: converting an original side-scan image to a ground coordinate system by taking the side-scan sonar towing depth and slant range record as constraints to obtain a geometrically compensated side-scan image; calculating a resolution proportionality coefficient according to the ratio of the actual ground resolutions of two types of sonars, and respectively constructing a multi-beam scale space and a side-scan scale space; giving a credibility weight to a candidate feature point according to the depth effectiveness of the multi-beam sounding grid, and screening an effective feature point set; taking a high-weight key point as a priority sampling object, solving a homography matrix H through an RANSAC algorithm, mapping multi-beam geographic coordinates to a side-scan image, and obtaining an accurate positioning result. In the positioning of underwater obstacles such as wrecks, the positioning error is reduced from 5.9 meters to 2.5 meters.

Description

Technical Field

[0001] This application relates to the field of marine surveying technology, and in particular to a method and system for accurate positioning of ground features in multibeam coordinate-assisted side-scan images. Background Technology

[0002] Underwater obstacle detection is a crucial prerequisite for marine engineering, channel dredging, and shipwreck salvage operations. Multibeam sonar and side-scan sonar are currently the two main methods for underwater obstacle detection. Multibeam bathymetry systems use multi-angle beam arrays to perform full-coverage scanning of the seabed, acquiring three-dimensional bathymetry data with high-precision geographic coordinates. However, limited by the physical aperture of the transducer array, its image resolution is relatively low, and significant echo-free voids exist in areas obscured by obstacles, making it difficult to reveal the fine details of targets. Side-scan sonar operates in a towed manner, generating high-resolution two-dimensional images by recording the intensity of acoustic reflections from the seabed. It can clearly show the contour details and texture features of obstacles such as shipwrecks. However, towed sonar is affected by factors such as water currents, hull sway, and insufficient accuracy of its own attitude sensor, resulting in significant positioning errors in side-scan images, making them unsuitable for direct use in precise obstacle coordinate measurement.

[0003] Existing methods for fusing multibeam and side-scan sonar images primarily employ wavelet transform, affine transform, or the Demons non-rigid registration algorithm to register the two types of images, transferring the geographic coordinate information of the multibeam images to the side-scan images. However, these methods treat multibeam images as ordinary grayscale images, neglecting the depth validity information attached to each grid cell in the multibeam data. They fail to distinguish the different impacts of measured effective areas, interpolated filled areas, and echo-free areas on registration accuracy, leading to the inclusion of false feature points from echo-free hole areas in the registration calculation, introducing systematic errors. Furthermore, the actual ground resolution of multibeam images and side-scan images can differ by up to 20 times. Existing methods use the same initial scale parameters for both types of images when constructing the scale space, resulting in inconsistent physical ground scales at the same layer number in the two scale spaces. Feature points of the same ground feature in the two images cannot be detected and described at the same physical scale, fundamentally limiting the matching success rate.

[0004] Due to the dual problems of resolution differences and unutilized depth information between multibeam images and side-scan images, existing methods relying solely on image grayscale features for registration introduce a large number of unreliable candidate points during the feature point extraction stage. Consequently, they cannot effectively distinguish between high-confidence and low-confidence feature point pairs during the matching stage, ultimately resulting in insufficient accuracy of the solved coordinate transformation matrix due to interference from mismatched points. Further analysis reveals that the slant-range imaging mechanism of side-scan images causes nonlinear geometric distortion between image pixel coordinates and ground coordinates. If feature points are directly extracted in the slant-range coordinate system and matched with multibeam orthophoto images, the geometric correspondence of feature points will produce systematic deviations in the image edge regions. This causes the linear fitting residual of the homography matrix to accumulate as the image width increases. These three problems overlap and collectively restrict the positioning accuracy of existing methods. Summary of the Invention

[0005] This application provides a method and system for accurate localization of ground features in multibeam coordinate-assisted side-scan images. It solves the problems in existing multibeam and side-scan sonar image registration methods, such as scale space misalignment due to resolution differences, false feature points mixed into the registration calculation due to ignoring depth validity information, and systematic deviations in geometric correspondence due to uncompensated side-scan slant range distortion. It reduces the localization error of underwater obstacles such as shipwrecks from 5.9 meters to 2.5 meters.

[0006] In a first aspect, this application provides a method for precise localization of ground features in multi-beam coordinate-assisted side-scan images, the method comprising: Step S1: Using the drag depth and slant range records of the side-scan sonar as constraints, transform the coordinates of each pixel in the original side-scan image from the slant range coordinate system to the ground coordinate system to obtain the geometrically compensated side-scan image. Step S2: Calculate the resolution scaling factor λ based on the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter, construct the Gaussian difference pyramids for the multibeam image and the geometrically compensated side-scan image respectively to obtain the multibeam scale space and the side-scan scale space. Step S3: Extract candidate feature points in the multibeam scale space and the side scan scale space respectively. Assign confidence weights to the candidate feature points based on the depth validity of each point in the multibeam sounding grid. Filter the valid feature point set using the confidence weights. Step S4: Prioritize sampling based on high-weight points in the effective feature point set, solve the homography matrix H between the multibeam image and the geometrically compensated side scan image using the RANSAC algorithm, and map the multibeam geographic coordinates to the side scan image via the homography matrix H to obtain accurate positioning results.

[0007] Secondly, this application provides a system for precise location of ground features from multi-beam coordinate-assisted side-scan images, the system comprising: The conversion module is used to convert the coordinates of each pixel in the original side-scan image from the slant range coordinate system to the ground coordinate system, based on the drag depth and slant range records of the side-scan sonar, to obtain a geometrically compensated side-scan image. The construction module is used to calculate the resolution scaling factor λ based on the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter, Gaussian difference pyramids are constructed for the multibeam image and the geometrically compensated side-scan image respectively, so as to obtain the multibeam scale space and the side-scan scale space. The filtering module is used to extract candidate feature points in the multibeam scale space and the side scan scale space respectively, assign confidence weights to the candidate feature points according to the depth validity of each point of the multibeam sounding grid, and obtain a set of effective feature points by filtering with the confidence weights. The mapping module is used to prioritize sampling objects with high-weight points in the effective feature point set, solve the homography matrix H between the multibeam image and the geometrically compensated side scan image using the RANSAC algorithm, and map the multibeam geographic coordinates to the side scan image via the homography matrix H to obtain accurate positioning results.

[0008] Thirdly, a device for precise positioning of ground features in multi-beam coordinate-assisted side-scan images is provided, comprising: a memory and at least one processor, wherein the memory stores instructions; the at least one processor invokes the instructions in the memory to cause the device for precise positioning of ground features in multi-beam coordinate-assisted side-scan images to execute the above-described method for precise positioning of ground features in multi-beam coordinate-assisted side-scan images.

[0009] Fourthly, a computer-readable storage medium is provided, wherein instructions are stored therein, which, when executed on a computer, cause the computer to perform the above-described method for precise location of ground features from multi-beam coordinate-assisted side-scan images.

[0010] The technical solution provided in this application uses the tow depth and slant range records of side-scan sonar as constraints to transform the pixel coordinates of the original side-scan image from the slant range coordinate system to the ground coordinate system, obtaining a geometrically compensated side-scan image. This fundamentally eliminates the nonlinear geometric distortion introduced by the slant range imaging mechanism of side-scan sonar, allowing subsequent feature point extraction to be performed in the ground coordinate system, avoiding systematic deviations in geometric correspondence caused by the accumulation of slant range distortion in the image edge region. Based on this, a resolution scaling factor λ is calculated according to the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter of the geometrically compensated side-scan image, Gaussian difference pyramids are constructed for each type of image. This ensures that the multibeam scale space and the side-scan scale space correspond to the same physical ground scale at the same group number and layer number, thereby guaranteeing that features with the same name are detected and described at physically comparable scale layers in both scale spaces. This solves the problem of feature point scale mismatch caused by the significant difference in resolution between the two types of sonar.

[0011] In the feature point selection stage, this application classifies candidate feature points into three categories—measured valid, interpolated filled, and echoless—based on the echo record status of each grid cell in the multibeam bathymetry grid, and assigns corresponding confidence weights. This introduces the depth validity information attached to the multibeam data into the feature point selection process, eliminating unreliable candidate points from echoless hole areas before they enter the matching calculation, fundamentally preventing false feature points from interfering with the coordinate transformation matrix solution. In this application, the RANSAC algorithm prioritizes feature points with high confidence weights for sampling, ensuring that the minimum sampling set for each iteration is primarily composed of measured valid feature points from the geographically most accurate region. This guarantees the initial fitting quality of the candidate homography matrix, accelerates inlier convergence, and results in a smaller residual in the final homography matrix H during the mapping process from multibeam geographic coordinates to the side-scan image. The reliability of the accurate positioning results obtained after mapping multibeam geographic coordinates to the side-scan image via the homography matrix H is significantly improved. Attached Figure Description

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

[0013] Figure 1 This is a schematic diagram of an embodiment of the method for precise localization of ground features using multibeam coordinate-assisted side-scan images in this application. Figure 2 This is a schematic diagram of the overall process of the method for precise localization of ground features in multi-beam coordinate-assisted side-scan images in the embodiments of this application; Figure 3The image shown is the original multibeam image obtained by multibeam sonar scanning the shipwreck area in the embodiments of this application. Figure 4 The images shown are the original side-scan images obtained by the side-scan sonar in the embodiments of this application, which scan the same shipwreck area. Figure 5 This is a schematic diagram illustrating the generation process of the Gaussian pyramid and the difference-of-Gaussian pyramid in the embodiments of this application; Figure 6 This is a schematic diagram of the step-by-step downsampling results during the construction of a Gaussian pyramid from a multibeam image in the embodiments of this application; Figure 7 This is a schematic diagram of DOG-scale spatial extremum point detection in an embodiment of this application; Figure 8 This is a schematic diagram illustrating the difference between discrete space extrema and the true extrema in continuous space in the embodiments of this application; Figure 9 This is a schematic diagram of the SIFT feature descriptor construction process in an embodiment of this application; Figure 10 This is a schematic diagram of the distribution of candidate feature points after extraction in the side-scan scale space in an embodiment of this application; Figure 11 This is a schematic diagram of the distribution of candidate feature points after extraction in the multi-beam scale space in an embodiment of this application. Figure 12 This is a schematic diagram of the image matching results after solving the homography matrix H using the confidence weighted RANSAC algorithm in this embodiment of the application; Figure 13 This is a schematic diagram of the matching results of the classic SIFT algorithm on the same multibeam image and side scan image dataset in the embodiments of this application; Detailed Implementation This application provides a method and system for precise location of ground features using multibeam coordinate-assisted side-scan images. The terms "first," "second," "third," "fourth," etc. (if present)," in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0014] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the method for precise localization of ground features using multi-beam coordinate-assisted side-scan images in this application includes: Step S1: Using the drag depth and slant range records of the side-scan sonar as constraints, transform the coordinates of each pixel in the original side-scan image from the slant range coordinate system to the ground coordinate system to obtain the geometrically compensated side-scan image. Specifically, the slant range coordinate system is a coordinate system established by the side-scan sonar with the towed fish as the origin and along the sound wave propagation path. Its horizontal axis represents the slant range of sound wave propagation rather than the horizontal distance, resulting in near-field targets being compressed and far-field targets being stretched in the image. The tow depth refers to the actual depth of the towed fish in the water, recorded in real time by the shipborne depth sensor. The slant range is recorded as the one-way sound wave propagation distance corresponding to each pixel. The horizontal distance value is calculated from the slant range and tow depth using the Pythagorean theorem. Invalid fill regions refer to blank grids caused by pixel position rearrangement after coordinate transformation. Masking these regions is used to exclude them in the subsequent feature point extraction stage, avoiding the detection of false extrema at locations without effective grayscale values.

[0015] Step S2: Calculate the resolution scaling factor λ based on the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter, construct the Gaussian difference pyramids for the multibeam image and the geometrically compensated side-scan image respectively to obtain the multibeam scale space and the side-scan scale space. Specifically, the beam footprint size refers to the actual area formed when the sonar beam illuminates the seabed. Due to the limitations of transducer array spacing, the beam footprint of multi-beam sonar is typically much larger than the pulse width of the linear array of side-scan sonar. The ratio of their actual ground resolution at the same water depth is λ. Using λ as the initial scale parameter for geometrically compensated side-scan images means that the side-scan images are physically aligned with the multi-beam images from the first layer of Gaussian blur. This ensures that the two difference pyramids correspond to the same physical ground scale at the same group number and layer number, thus making the extracted extrema points at that scale comparable across images.

[0016] Step S3: Extract candidate feature points in the multibeam scale space and the side scan scale space respectively. Assign confidence weights to the candidate feature points based on the depth effectiveness of each point in the multibeam sounding grid. Filter the valid feature point set by the confidence weights. Specifically, depth validity refers to whether each grid cell in the multibeam echo sounding grid contains a bathymetry value derived from actual seabed echoes. Measured valid cells are weighted 1.0 because real topographic measurements exist at that location, and the extracted feature point coordinates directly correspond to geographic coordinates. Interpolated filler cells are weighted 0.5 because the bathymetry value at that location is interpolated from surrounding measured points, resulting in lower coordinate accuracy than actual measurements. Cells without echoes are weighted 0 and directly eliminated because there are no topographic constraints at that location, making the feature point coordinates unreliable. Side-scan candidate feature points are weighted inversely proportional to distance because side-scan images themselves do not carry depth information; the reliability of their feature points is indirectly reflected by the density of their neighboring measured cells in the multibeam echo sounding grid. Closer proximity indicates more complete multibeam coverage in the area.

[0017] Step S4: Prioritize sampling based on high-weight points in the effective feature point set. Solve the homography matrix H between the multibeam image and the geometrically compensated side scan image using the RANSAC algorithm. Map the multibeam geographic coordinates to the side scan image via the homography matrix H to obtain accurate positioning results.

[0018] Specifically, the priority sampling targets are feature points with a weight value of no less than 1.0 because these points come from the effective area of ​​multibeam scanning, and their planar coordinates have the highest geographic accuracy. Using these points as the sampling source results in the smallest residual in the candidate homography matrix and the fastest convergence. RANSAC requires at least 4 sets of corresponding image points in each iteration because the homography matrix H contains 8 degrees of freedom. Each set of corresponding image points provides 2 equations, and 4 sets constitute the minimum solvable condition for the system of equations. A distance threshold is used to determine inliers; this threshold is set to 4 pixels, derived from the statistical results of the standard error after multibeam and side-scan image registration in the original experimental data.

[0019] Figure 2 This is a schematic diagram of the overall process of the method for accurate localization of ground features in multi-beam coordinate-assisted side-scan images in this application embodiment. It shows the complete processing flow from input image through Gaussian pyramid construction, Gaussian difference pyramid construction, key point localization, principal direction calculation, SIFT feature descriptor generation to feature vector output, and based on the feature sample set of multi-beam image and side-scan image through random sampling, homography matrix parameter calculation, interior point statistics to output homography matrix parameters and interior point set.

[0020] In one specific embodiment, step S1 includes: Based on the drag depth record of the side-scan sonar and the slant distance record of each pixel, the corresponding horizontal distance value of each pixel in the original side-scan image is calculated to obtain the slant distance-horizontal distance mapping table. Based on the slant range-horizontal distance mapping table, the slant range coordinates of each pixel in the original side scan image are replaced with horizontal distance coordinates. The gray values ​​of the pixels after coordinate replacement are resampled by bilinear interpolation to obtain the ground range domain side scan image. Based on the track coordinates and heading angle records of the side-scan sonar, the ship coordinate system coordinates of each pixel in the ground range domain side-scan image are transformed to the geodetic coordinate system to obtain the georegistered side-scan image. Invalid fill regions caused by coordinate transformation in georegistered side scan images are masked and marked, and the georegistered side scan images with mask marking are used as geometrically compensated side scan images.

[0021] Specifically, the tow depth refers to the vertical distance between the side-scan sonar towed fish and the water surface, recorded in real time by the ship's pressure sensor and stored frame by frame during navigation. The slant range refers to the one-way propagation distance of the sound waves emitted by the towed fish from the seabed target to the return, obtained by multiplying the round-trip time of the sound wave corresponding to each pixel by the speed of sound in water and then dividing by two. The horizontal distance value is calculated by the slant range and tow depth using the Pythagorean theorem, that is, the horizontal distance is equal to the square of the slant range minus the square of the tow depth, and the positive square root of the result. The slant range-horizontal distance mapping table is a lookup table that stores the above calculation results pixel by pixel. After replacing the slant range coordinates of each pixel with the horizontal distance coordinates based on this mapping table, the spacing of each pixel in the horizontal direction is no longer uniform. The gray values ​​of the pixels originally uniformly sampled need to be redistributed to a uniformly spaced ground coordinate grid. Bilinear interpolation resampling is to take the weighted average gray value of the four original pixels surrounding each position in the new grid. The weights are determined by the bidirectional distance between the new position and the center of the four original pixels, thus obtaining the ground range domain side-scan image.

[0022] Track coordinates refer to the geodetic latitude and longitude coordinates of the side-scan sonar towed fish at each sampling time, recorded by the shipborne differential GPS and aligned with the sonar timestamp; heading angle refers to the angle between the towed fish's direction of travel and true north, recorded by the shipborne attitude sensor. The ship's coordinate system has the towed fish as the origin, the direction of travel as the vertical axis, and the perpendicular direction of travel as the horizontal axis. To transform the ship's coordinate system coordinates of each pixel in the ground distance domain side-scan image to the geodetic coordinate system, the horizontal distance and longitudinal track distance of each pixel need to be multiplied by the sine and cosine components of the heading angle, respectively, and then superimposed onto the towed fish's geodetic coordinate system to obtain the geodetic latitude and longitude coordinates of that pixel, i.e., the georegistered side-scan image. After coordinate transformation, some grid positions become blank because no original pixels fall into them. These positions are marked as invalid filling areas. The mask markings are stored in the form of a binary mask matrix. Valid pixel positions are assigned a value of 1, and invalid positions are assigned a value of 0. During subsequent feature point extraction, extreme point detection is only performed in areas where the mask value is 1. The georegistration side scan image after mask marking is completed is used as the geometric compensation side scan image.

[0023] Figure 3The image shows the original multibeam image obtained by multibeam sonar scanning the shipwreck area in this embodiment of the application. The green area in the image is the normal sounding coverage area of ​​the multibeam sonar, and the black area is the echo-free void area caused by the shipwreck blocking the sound waves. This void area is a visual representation of the echo-free unit in the multibeam sounding grid and is the data source basis for assigning a weight of 0 to the echo-free unit in the depth effectiveness classification grid.

[0024] Figure 4 The image shown is the original side-scan image obtained by the side-scan sonar scanning the same shipwreck area in the embodiment of this application. The yellow area in the image represents the distribution of acoustic reflection intensity on the seabed. The outline of the shipwreck and its acoustic shadow are clearly visible, and the image has rich texture details. This image is the original side-scan image input for the slant range-horizontal distance mapping table calculation and bilinear interpolation resampling processing in step S1.

[0025] In one specific embodiment, step S2 includes: Based on the beam footprint size of the multibeam sonar and the side-scan sonar, the actual ground resolution of the two under the same water depth is calculated, and the ratio of the actual ground resolution is used as the resolution scaling factor λ. Using the resolution scaling factor λ as the initial scale parameter of the geometrically compensated side-scan image and 1 as the initial scale parameter of the multibeam image, the multibeam image and the geometrically compensated side-scan image are subjected to progressive Gaussian blurring and downsampling respectively to obtain the multibeam Gaussian pyramid and the side-scan Gaussian pyramid. Perform differential operations on adjacent layers in the same group of the multibeam Gaussian pyramid to obtain the multibeam differential pyramid; perform differential operations on adjacent layers in the same group of the side-scan Gaussian pyramid to obtain the side-scan differential pyramid. By using the multibeam differential pyramid and the side-scan differential pyramid as their respective extreme point detection spaces, the multibeam scale space and the side-scan scale space are obtained.

[0026] Specifically, the beam footprint size refers to the lateral width of the acoustic illumination area formed on the seabed surface when the sonar beam illuminates the seabed. The beam footprint of multi-beam sonar is determined by the physical aperture of the transducer array and the water depth, with a lateral resolution of approximately 40 centimeters at a water depth of 40 meters. Side-scan sonar uses linear pulse imaging, and its lateral resolution is determined by the pulse width and sound velocity, with a lateral resolution of approximately 2 centimeters at the same water depth. Dividing the actual ground resolution of the side-scan sonar by the actual ground resolution of the multi-beam sonar yields the resolution scaling factor λ, which is approximately 20 under the aforementioned parameters. λ is used as the initial scale parameter for the geometrically compensated side-scan image. This means that when applying the first layer of Gaussian blur to the geometrically compensated side-scan image, the standard deviation of the Gaussian kernel is set to λ times the reference scale, while when applying the first layer of Gaussian blur to the multi-beam image, the standard deviation is set to 1 times the reference scale, with the reference scale uniformly set to 0.5 pixels. This setting ensures that both images are blurred at the same physical ground scale from the first layer, and the physical ground sampling intervals corresponding to the same group and layer numbers of the two pyramids are consistent.

[0027] The Gaussian pyramid is set to 4 groups, with 5 layers in each group. The standard deviation scaling factor k of the Gaussian kernel between adjacent layers is taken as the positive root of the number of layers in each group, i.e., k equals 2 to the power of 1 / 5. Successive Gaussian blurring refers to convolving the image of the previous layer within the same group with a Gaussian kernel that increases by a factor of k in standard deviation. Successive downsampling refers to downsampling the final layer image of the previous group by a factor of 2 between groups as the input of the first layer of the next group, and so on, to obtain the multi-beam Gaussian pyramid and the side-scan Gaussian pyramid. Differential operation refers to subtracting the pixel grayscale values ​​of adjacent layers within the same group bit by bit to obtain the Gaussian difference image of each layer. The set of all adjacent difference images within a group constitutes the difference pyramid of that group. The four difference pyramids are merged to form the multi-beam difference pyramid and the side-scan difference pyramid. The aforementioned difference pyramid is directly used as the operation space for extreme point detection, namely the multi-beam scale space and the side-scan scale space. Subsequently, each pixel is compared with its 26 three-dimensional neighboring points in this space to detect extreme points.

[0028] Figure 5 This is a schematic diagram of the generation process of Gaussian pyramids and difference-of-Gaussian pyramids in the embodiments of this application. It shows the construction logic and inter-layer correspondence of generating corresponding difference-of-Gaussian pyramids by subtracting adjacent layers in the same group after the image is successively blurred and downsampled to form multiple sets of multi-layer Gaussian pyramids.

[0029] Figure 6This is a schematic diagram of the step-by-step downsampling results of the multibeam image during the construction of the Gaussian pyramid in this embodiment of the application. The diagram shows the original image and the three-level downsampled images corresponding to scaling factors of 0.5, 0.25 and 0.125, respectively, from the top left to the bottom right. It shows the process of the change of the texture of the ground features in each layer of the image after applying step-by-step Gaussian blur and downsampling to the multibeam image with the scale reduction, with the resolution scaling factor λ as the initial scale parameter. The diagram also shows the image content of the first layer of each group of the multibeam Gaussian pyramid.

[0030] In one specific embodiment, step S3 involves extracting candidate feature points in the multi-beam scale space and the side-scan scale space, respectively, including: Each pixel in the multi-beam scale space and the side-scan scale space is compared one by one with 26 neighboring pixels in the same group, the same layer and adjacent layers. Pixels whose gray values ​​are strictly greater than or strictly less than all 26 neighboring pixels are marked as initial extreme points, thus obtaining the initial extreme point set of multi-beam and the initial extreme point set of side-scan. For each point in the initial extreme point set of the multi-beam and the initial extreme point set of the side scan, the local surface of its scale space is fitted based on a three-dimensional quadratic function. The coordinates of the extreme points are corrected at the sub-pixel level according to the fitting offset. Points with contrast below a set threshold and edge response points with principal curvature ratio exceeding a set threshold are removed to obtain the precise feature point set of the multi-beam and the precise feature point set of the side scan. Taking each point in the multi-beam precise feature point set and the side-scan precise feature point set as the center, the pixel gradient direction histogram in its 3σ neighborhood is calculated. The peak direction of the histogram is taken as the main direction of the corresponding feature point, and the direction with more than 80% of the peak is taken as the auxiliary direction, thus obtaining the multi-beam feature point set with direction and the side-scan feature point set with direction. The directional multibeam feature point set and the directional side scan feature point set are concatenated into an 8-directional gradient histogram within a 4×4 sub-region of each point to form a 128-dimensional description vector. The 128-dimensional description vector is then L2 normalized to obtain multibeam candidate feature points and side scan candidate feature points, which together constitute candidate feature points.

[0031] Specifically, the 26 neighboring pixels are composed of the 8 adjacent pixels in the same layer as the current pixel, the 9 pixels in the corresponding position in the previous layer, and the 9 pixels in the corresponding position in the next layer, forming a three-dimensional cubic neighborhood. Three-dimensional quadratic function fitting refers to constructing a three-dimensional quadratic polynomial in the local gray-level distribution of the current initial extreme point as the center. By taking the partial derivative of this polynomial and setting it to zero, the continuous offset of the extreme point in spatial and scale coordinates is obtained. This offset is superimposed on the discrete coordinates of the initial extreme point to obtain the precise coordinates of the sub-pixel. The contrast threshold is set to 0.03, meaning that if the gray-level response value at the fitted extreme point is lower than 0.03, it is discarded. This value comes from the original parameters of the SIFT algorithm and is used to eliminate false extreme points with low response caused by noise. The principal curvature ratio threshold is set to 10, meaning that if the ratio of the two principal curvatures of the Hessian matrix exceeds 10, the corresponding point is identified as an edge response point and discarded. This value also comes from the original parameters of the SIFT algorithm and is used to eliminate unstable extreme points located at the image edges.

[0032] The 3σ neighborhood refers to a circular region centered on the precise feature point with a radius three times the Gaussian standard deviation σ corresponding to the scale layer of that point. The gradient directions of all pixels within this region are divided into 36 directional intervals, each covering 10 degrees. The cumulative gradient magnitude values ​​falling into each interval are used to construct a direction histogram. The direction corresponding to the maximum peak in the histogram is designated as the primary direction, and the directions corresponding to the remaining peaks exceeding 80% of the maximum peak are designated as secondary directions. This 80% ratio originates from the original setting of the SIFT algorithm and is used to assign multiple stable directions to a single feature point in rotational scenarios. The 4×4 sub-region refers to a 16×16 pixel neighborhood centered on the precise feature point and rotated and aligned along the primary direction. This neighborhood is evenly divided into 16 sub-regions (4 rows, 4 columns). Gradient magnitude histograms for 8 directional intervals are calculated within each sub-region, with each interval covering 45 degrees. Each of the 16 sub-regions generates an 8-dimensional vector, which is then concatenated to obtain a 128-dimensional descriptive vector. L2 normalization involves dividing the 128-dimensional description vector by its own Euclidean norm to make the normalized vector have a magnitude of 1. After normalization, components with values ​​exceeding 0.2 in each dimension are truncated to 0.2 before L2 normalization is performed again. The 0.2 truncation threshold comes from the original setting of the SIFT algorithm and is used to suppress the influence of nonlinear illumination changes on the description vector.

[0033] Figure 7 This is a schematic diagram of DOG-scale spatial extremum point detection in an embodiment of this application, showing the three-dimensional comparison range and spatial distribution relationship between the currently detected pixel and 26 neighboring pixels in the same group, layer, and adjacent layers above and below.

[0034] Figure 8 This is a schematic diagram illustrating the difference between discrete space extrema and the true extrema in continuous space in the embodiments of this application. Figure 8 This diagram illustrates the positional deviation between discrete spatial extrema and true continuous spatial extrema in the keypoint localization step of the SIFT algorithm. The curves in the figure represent the variation of pixel grayscale values ​​in continuous space. All points are represented by solid black dots. Points labeled as detected extrema are those directly detected on the discrete sampling grid. Due to the discrete sampling method used to store the image in the computer, the positions of these points have quantization errors and are not the true extrema of the curve. Points labeled as true extrema are obtained by fitting a surface to the neighborhood of the discrete extrema using a three-dimensional quadratic function, taking the partial derivative of the fitted polynomial, and setting it to zero. These points have sub-pixel-level accuracy. The positions of the two types of points do not overlap in the figure, and their spatial deviation intuitively reflects the quantization error introduced by discrete sampling and the necessity of sub-pixel-level correction.

[0035] Figure 9 This is a schematic diagram of the SIFT feature descriptor construction process in this application embodiment. (a) shows a circular neighborhood region centered on a keypoint and defined in the original image coordinate system (xy axis). The red arrow indicates the main direction of the keypoint, the green circle is the boundary of the neighborhood range, and the internal 4×4 grid is the sub-region to be statistically analyzed. This corresponds to the process of calculating the pixel gradient direction histogram within the 3σ neighborhood of the feature point and using the peak direction of the histogram as the main direction of the corresponding feature point. (b) shows the result after rotating and aligning the neighborhood region in (a) according to the main direction of the keypoint. That is, the entire 4×4 grid region rotates with the main direction, making the coordinate axes of the descriptor consistent with the main direction, thereby achieving rotation invariance. This corresponds to the process of defining the neighborhood region after rotating and aligning the feature point according to the main direction. (c) shows the process of defining the neighborhood region in the rotated and aligned coordinate system. The specific method of dividing the neighborhood into 4×4 sub-regions is shown in the red and green boxes, which respectively indicate the size parameters of the sub-regions. × (d) shows the result after calculating the 8-directional gradient histogram for each sub-region. The arrow combination in each cell is the visualization of the 8-directional gradient histogram for that sub-region. The red solid dots are the key point locations. Each of the 16 sub-regions generates an 8-dimensional vector, which is then concatenated to obtain a 128-dimensional description vector. This corresponds to the process of concatenating the 8-directional gradient histograms within the 4×4 sub-regions into a 128-dimensional description vector.

[0036] Figure 10This is a schematic diagram of the distribution of candidate feature points extracted in the side-scan scale space in the embodiments of this application. Each circle in the figure represents a candidate feature point in the side scan. The size of the circle corresponds to the scale parameter of the scale layer where the feature point is located. The direction of the arrow inside the circle corresponds to the main direction of the feature point. The number of candidate feature points in the figure is dense and widely distributed in various regions of the image, which reflects the characteristics of the side-scan image having rich texture and a large number of extractable feature points.

[0037] Figure 11 This is a schematic diagram showing the distribution of candidate feature points extracted in the multi-beam scale space in an embodiment of this application. The meanings of the circles in the diagram are as follows: Figure 10 Similarly, candidate feature points are mainly concentrated on the outline edge of the shipwreck target and in areas with significant changes in the surrounding terrain. There are few candidate feature points in the hollow areas. This distribution pattern directly reflects the spatial distribution of each unit category in the depth effectiveness classification grid in the multibeam image, and is the data basis for assigning confidence weights to multibeam candidate feature points according to depth effectiveness.

[0038] In one specific embodiment, step S3 assigns confidence weights to candidate feature points based on the depth validity of each point in the multibeam echo sounding grid, including: Based on the echo record status of each grid cell in the multibeam bathymetry grid, the grid cells are divided into three categories: measured effective cells, interpolated filling cells, and no-echo cells, thus obtaining the depth effectiveness classification grid. The planar coordinates of the multi-beam candidate feature points are projected onto the depth effectiveness classification grid. According to the category of the grid cell in which each feature point falls, the multi-beam candidate feature points are assigned values ​​one by one according to the rule that the measured effective cells are assigned a weight of 1.0, the interpolated filling cells are assigned a weight of 0.5, and the no-echo cells are assigned a weight of 0, so as to obtain the weight set of the multi-beam candidate feature points. The planar coordinates of the side-scan candidate feature points are projected onto the depth validity classification grid. Based on the grid distance between the grid cell where each feature point falls and the nearest measured valid cell, the weight value of each point is calculated inversely proportional to the distance. The side-scan candidate feature points are assigned values ​​one by one to obtain the side-scan candidate feature point weight set. The confidence weights of the candidate feature points are obtained by merging the weight set of the multi-beam candidate feature points with the weight set of the side-scan candidate feature points.

[0039] Specifically, a multibeam bathymetry grid refers to a grid data structure formed by dividing the seabed plane of the survey area into regular rectangular grids at fixed intervals, assigning each multibeam bathymetry point to its corresponding grid cell according to its planar coordinates, and recording its depth value and data source marker. Echo recording status refers to the data source type of each grid cell during the multibeam acquisition process: if the cell contains at least one actual echo return value, it is marked as a valid measured cell; if the cell has no actual measured value but the depth value is interpolated from surrounding measured cells, it is marked as an interpolation-filled cell; if the cell contains neither actual measured value nor any interpolated surrounding data, it is marked as a cell without echo. These three classifications together constitute the depth validity classification grid. When assigning weights to candidate feature points in a multibeam array, the measured effective unit is weighted 1.0 because there is a real acoustic echo constraint at this location, and the feature point's planar coordinates directly correspond to the geographic coordinates, resulting in the highest geographic accuracy. The interpolated fill unit is weighted 0.5 because the depth value at this location is calculated from nearby measured points, and although the coordinate accuracy is lower than the measured values, it still has some reference value. The no-echo unit is weighted 0 because there are no terrain constraints at this location, and the feature point coordinates do not have geographic reliability, so they are directly excluded.

[0040] The process of projecting side-scan candidate feature points onto the depth-effectiveness classification grid is as follows: The difference between the geodetic coordinates of each side-scan candidate feature point in the geometrically compensated side-scan image and the coordinates of the origin of the depth-effectiveness classification grid is calculated, then divided by the grid spacing. The integer part of the result is the grid row and column index where the feature point falls. Starting from this index, neighboring cells in the classification grid are traversed in breadth-first order to find the nearest measured effective cell. The grid distance *d* between the two cells is calculated; this grid distance refers to the Euclidean distance between the center points of two grid cells, expressed in units of grid spacing. The weight values ​​of the side-scan candidate feature points are calculated according to the inverse distance rule: the weight value equals 1 divided by *d* plus 1. When *d* is 0, the weight value is 1.0. As *d* increases, the weight value monotonically decreases, approaching 0. Adding 1 to the denominator avoids division by zero when *d* is 0 and maintains consistency with the 1.0 weight value for multi-beam measured effective cells. After calculating the side-scan candidate feature points one by one according to the above rules, the weight set of side-scan candidate feature points is obtained. This weight set is then merged with the weight set of multi-beam candidate feature points according to the feature point index to form a confidence weight that covers all candidate feature points.

[0041] In one specific embodiment, step S3, which involves selecting the effective feature point set by confidence weight, includes: Based on the confidence weight, multibeam candidate feature points with a weight value of 0 are directly removed from the candidate feature points, and side scan candidate feature points with a weight value of less than 0.5 are directly removed from the candidate feature points, resulting in a weighted filtered feature point set. For the 128-dimensional description vectors of multibeam candidate feature points and side-scan candidate feature points in the feature point set after weight filtering, calculate the Euclidean distance between each pair. Feature point pairs with a ratio of nearest neighbor Euclidean distance to second nearest neighbor Euclidean distance of less than 0.6 are marked as initial matching pairs to obtain the initial matching feature point set. Based on the confidence weight, the Euclidean distance ratio of each matching pair in the initial matching feature point set is adjusted by weight. Matching pairs whose adjusted Euclidean distance ratio is still lower than 0.6 are retained to obtain the weight-adjusted matching feature point set, which is then used as the effective feature point set.

[0042] Specifically, points with a weight value of 0 are removed from multibeam candidate feature points because a weight value of 0 corresponds to no echo cell, and there are no topographic survey constraints at this location. The feature point coordinates do not have geographic reliability, and retaining such points will introduce invalid image point pairs in subsequent matching. Points with a weight value lower than 0.5 are removed from side scan candidate feature points because a weight value lower than 0.5 means that the grid distance between the location of the side scan feature point and the nearest measured effective cell is greater than 1, that is, there is no adjacent multibeam measured coverage. The geographic coordinate reliability of this point is insufficient to support the subsequent coordinate mapping accuracy requirements. After weighted filtering, the calculation range for the pairwise Euclidean distance between the 128-dimensional descriptive vectors of multi-beam candidate feature points and side-scan candidate feature points in the feature point set is as follows: For each multi-beam candidate feature point, traverse all side-scan candidate feature points and calculate the Euclidean distance between their 128-dimensional descriptive vectors respectively. Sort them in ascending order of distance, take the one with the smallest distance as the nearest neighbor, and the one with the second smallest distance as the second nearest neighbor. The ratio of the nearest neighbor's Euclidean distance to the second nearest neighbor's Euclidean distance is the distance ratio of the multi-beam candidate feature point. The distance ratio threshold is set to 0.6, which comes from the original parameter setting of the SIFT algorithm. A ratio lower than 0.6 indicates that the distance difference between the nearest neighbor and the second nearest neighbor is sufficiently significant, and the match is unique. Feature point pairs that meet this condition are marked as initial matching pairs, forming the initial matching feature point set.

[0043] The specific calculation method for weight correction is as follows: For each matching pair in the initial matching feature point set, the product of the confidence weight values ​​of the multi-beam candidate feature points and the side-scan candidate feature points in the matching pair is taken as the joint weight of the matching pair. The original Euclidean distance ratio is divided by the joint weight to obtain the corrected Euclidean distance ratio. A larger joint weight indicates that the geographical confidence of the feature points at both ends of the matching pair is higher. After dividing by the joint weight, the ratio decreases, and the matching pair is more likely to pass the 0.6 threshold screening. A smaller joint weight indicates that the geographical confidence of at least one feature point is lower. After correction, the ratio increases, and the matching pair is more likely to be eliminated. The corrected threshold is still 0.6, consistent with the initial screening threshold, to ensure that the final retained matching pairs meet the requirements in both the description vector similarity and geographical confidence dimensions. All matching pairs that pass the corrected threshold screening constitute the weight-corrected matching feature point set, i.e., the effective feature point set.

[0044] In one specific embodiment, step S4 includes: Based on the confidence weight of each point in the effective feature point set, feature points with a weight value of not less than 1.0 are assigned to the priority sampling point set, and feature points with a weight value of less than 1.0 are assigned to the candidate sampling point set. The priority sampling point set is used as the first sampling source and the candidate sampling point set is used as the supplementary sampling source to obtain the weighted sampling sequence. In each iteration of the RANSAC algorithm, four pairs of image points with the same name are extracted from the effective feature point set based on the weighted sampling sequence and the candidate feature points of the multi-beam candidate feature points and the side scan candidate feature points. The candidate homography matrix is ​​solved based on the pair of image points with the same name. The projection error of each point to the candidate homography matrix is ​​calculated by traversing all point pairs in the effective feature point set. Point pairs with projection errors lower than a set distance threshold are marked as inliers. The number of inliers is counted. The candidate homography matrix with the largest number of inliers is updated as the current best homography matrix. The iteration is repeated for a set number of times K to obtain the homography matrix H. Substitute the geodetic coordinates of the multibeam candidate feature points in the effective feature point set into the homography matrix H, calculate the corresponding pixel coordinates of each point in the geometrically compensated side-scan image, and obtain the coordinate mapping result. The planar residuals between each mapping point and its corresponding side-scan candidate feature points in the coordinate mapping results are statistically analyzed. Mapping points with residuals exceeding a set threshold are removed, and the geodetic coordinates of the remaining mapping points are used as the accurate positioning results.

[0045] Specifically, the priority sampling point set partitioning threshold is set to 1.0, corresponding to points from actual measured effective units among the multibeam candidate feature points. The geodetic coordinates of these points are directly determined by multibeam measured echoes, resulting in the highest geographic accuracy. Therefore, the candidate homography matrix constructed using these points as the primary sampling source has the smallest residual. The weighted sampling sequence is constructed as follows: all points in the priority sampling point set are arranged at the beginning of the sequence, and all points in the candidate sampling point set are arranged in descending order of confidence weight at the end of the sequence. Each iteration prioritizes sampling from the beginning of the sequence; if there are fewer than four sets of points in the beginning, points from the end are used to supplement them. At least four pairs of corresponding image points are sampled in each iteration because the homography matrix H is a 3×3 homogeneous matrix with 8 effective degrees of freedom. Each pair of corresponding image points provides two linear equations, and four pairs constitute the minimum solvable linear equation system of eight equations. Based on this equation system, the candidate homography matrix is ​​obtained by solving the direct linear transformation method. Projection error refers to the Euclidean distance between the predicted side-scan coordinates obtained by substituting the coordinates of multi-beam candidate feature points into the candidate homography matrix and the actual side-scan candidate feature point coordinates. The distance threshold is set to 4 pixels, derived from the statistical median of the registration residuals between multi-beam and side-scan images in the original experimental data. Point pairs with a projection error of less than 4 pixels are marked as inliers. The number of iterations K is determined according to the RANSAC theoretical formula, as follows: Where p is the probability that at least one sample will be an interior point, taken as 0.99, and e is the initial estimate of the interior point rate, taken as 0.5. The calculated K is 72 times, but in practice it is set to K equal to 30 times, which achieves a balance between computational efficiency and solution stability.

[0046] The geodetic coordinates of multibeam candidate feature points in the effective feature point set are substituted into the homography matrix H in homogeneous coordinate form. This involves representing the geodetic coordinates as a three-dimensional homogeneous column vector, multiplying it by H on the left, and then dividing the resulting vector by its third component and taking the first two components. These first two components represent the corresponding pixel coordinates of the multibeam candidate feature point in the geometrically compensated side-scan image. The set of pixel coordinates after mapping all multibeam candidate feature points constitutes the coordinate mapping result. The planar residual refers to the Euclidean distance between the pixel coordinates of each mapped point and the actual pixel coordinates of its corresponding side-scan candidate feature point. This Euclidean distance is multiplied by the ground resolution of the geometrically compensated side-scan image to convert it into a ground distance residual in meters. A rejection threshold of 2.5 meters is set; mapped points with residuals exceeding 2.5 meters are considered unreliable points for coordinate correction and are rejected. The geodetic coordinates of the remaining mapped points, corrected by the homography matrix H, are then used to inversely calculate the geodetic coordinates from the side-scan image coordinates, which are then output as the precise positioning result.

[0047] Tables 1 and 2 are shown below: Table 1 is a summary of the setting values ​​of key parameters of the RANSAC algorithm in the embodiments of this application. The table lists the specific values ​​of four parameters: the number of iterations K is set to 30 times, the minimum number of model points is set to 4 pairs of corresponding image points, the inlier distance threshold is set to 4 pixels, and the inlier proportion threshold is set to 0.02. The above parameters are determined based on the statistical results of the registration residuals of the measured multibeam image and side scan image datasets. The minimum number of model points of 4 corresponds to the minimum solvable condition of the 8 degrees of freedom parameters of the homography matrix. The inlier distance threshold of 4 pixels corresponds to the statistical median of the registration residuals in the measured data. The number of iterations of 30 times meets the computational efficiency requirements of actual engineering operations while ensuring the stability of the solution.

[0048] Table 2 presents a statistical comparison of the quantitative accuracy of the RANSAC algorithm based on confidence weight and the classic SIFT algorithm on the same shipwreck area dataset in the embodiments of this application. The table uses multibeam measured coordinates as the benchmark and lists the north and east coordinate positioning results of each pair of corresponding image points for the two algorithms. The mean square error of the positioning for each algorithm is given in the last row. The mean square error of the method in this application is 2.5 meters, and the mean square error of the classic SIFT algorithm is 5.9 meters. The difference in the coordinate residuals of each group intuitively reflects the degree of improvement in the accuracy of homography matrix solution by introducing the confidence weight of depth validity.

[0049] Table 1 Distribution of key RANSAC-SIFT parameters Table 2 Quantitative Experiment Statistics of Two Types of Algorithms Figure 12 This is a schematic diagram of the image matching results after solving the homography matrix H based on the confidence weighted RANSAC algorithm in the embodiments of this application. The left side of the figure is the multi-beam image, and the right side is the side scan image. The red line represents the 10 pairs of image points with the same name that are finally retained after effective feature point set screening, weight correction matching and RANSAC in-point screening. The two ends of the line are the coordinate positions of the multi-beam candidate feature points and the side scan candidate feature points in their respective images.

[0050] Figure 13 This is a schematic diagram illustrating the matching results of the classic SIFT algorithm on the same multibeam image and side-scan image datasets in an embodiment of this application. The red lines in the diagram represent only three pairs of identical image points obtained by the classic SIFT algorithm. Figure 12 The 10 pairs of identical image points in the middle form an intuitive comparison. Combined with the quantitative statistical results in Table 2, which show that the mean square error has decreased from 5.9 meters to 2.5 meters, the improvement effect of the depth validity and credibility weight constraint introduced in this application on the matching results is verified.

[0051] The above describes the method for precise localization of ground features using multi-beam coordinate-assisted side-scan images in the embodiments of this application. The following describes the system for precise localization of ground features using multi-beam coordinate-assisted side-scan images in the embodiments of this application. One embodiment of the system for precise localization of ground features using multi-beam coordinate-assisted side-scan images in the embodiments of this application includes: The conversion module is used to convert the coordinates of each pixel in the original side-scan image from the slant range coordinate system to the ground coordinate system, based on the drag depth and slant range records of the side-scan sonar, to obtain a geometrically compensated side-scan image. The construction module is used to calculate the resolution scaling factor λ based on the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter, Gaussian difference pyramids are constructed for the multibeam image and the geometrically compensated side-scan image respectively, so as to obtain the multibeam scale space and the side-scan scale space. The filtering module is used to extract candidate feature points in the multibeam scale space and the side scan scale space respectively, assign confidence weights to the candidate feature points according to the depth validity of each point of the multibeam sounding grid, and obtain a set of effective feature points by filtering with the confidence weights. The mapping module is used to prioritize sampling objects with high-weight points in the effective feature point set, solve the homography matrix H between the multibeam image and the geometrically compensated side scan image using the RANSAC algorithm, and map the multibeam geographic coordinates to the side scan image via the homography matrix H to obtain accurate positioning results.

[0052] This invention also provides a device for precise location of ground features from multi-beam coordinate-assisted side-scan images. This device can be a server. The device includes a processor, memory, display screen, input device, network interface, and database connected via a system bus. The processor, designed as a computer, provides computational and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database stores the data corresponding to this embodiment. The network interface communicates with external terminals via a network connection. The computer program, when executed by the processor, implements the above-described method.

[0053] The present invention also provides a computer-readable storage medium, which can be a non-volatile computer-readable storage medium or a volatile computer-readable storage medium, wherein the computer-readable storage medium stores instructions that, when the instructions are executed on a computer, cause the computer to perform the steps of the method for precise localization of ground features in multi-beam coordinate-assisted side-scan images.

[0054] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0055] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a multi-beam coordinate-assisted side-scan image-based ground feature precise positioning device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

Claims

1. A method for precise ground feature localization from multi-beam coordinate-assisted side-scan images, characterized in that, The method includes: Step S1: Using the drag depth and slant range records of the side-scan sonar as constraints, transform the coordinates of each pixel in the original side-scan image from the slant range coordinate system to the ground coordinate system to obtain the geometrically compensated side-scan image. Step S2: Calculate the resolution scaling factor λ based on the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter, construct the Gaussian difference pyramids for the multibeam image and the geometrically compensated side-scan image respectively to obtain the multibeam scale space and the side-scan scale space. Step S3: Extract candidate feature points in the multibeam scale space and the side scan scale space respectively. Assign confidence weights to the candidate feature points based on the depth validity of each point in the multibeam sounding grid. Filter the valid feature point set using the confidence weights. Step S4: Prioritize sampling based on high-weight points in the effective feature point set, solve the homography matrix H between the multibeam image and the geometrically compensated side scan image using the RANSAC algorithm, and map the multibeam geographic coordinates to the side scan image via the homography matrix H to obtain accurate positioning results.

2. The method for precise localization of ground features in multi-beam coordinate-assisted side-scan images according to claim 1, characterized in that, Step S1 includes: Based on the drag depth record of the side-scan sonar and the slant distance record of each pixel, the corresponding horizontal distance value of each pixel in the original side-scan image is calculated to obtain the slant distance-horizontal distance mapping table. Based on the slant range-horizontal distance mapping table, the slant range coordinates of each pixel in the original side scan image are replaced with horizontal distance coordinates, and the gray values ​​of the pixels after coordinate replacement are resampled by bilinear interpolation to obtain the ground distance domain side scan image. Based on the track coordinates and heading angle records of the side-scan sonar, the ship coordinate system coordinates of each pixel in the ground distance domain side-scan image are transformed to the geodetic coordinate system to obtain the georegistered side-scan image. The invalid filled areas caused by coordinate transformation in the geo-registered side scan image are masked and marked, and the geo-registered side scan image after masking is used as the geometrically compensated side scan image.

3. The method for precise localization of ground features from multi-beam coordinate-assisted side-scan images according to claim 1, characterized in that, Step S2 includes: Based on the beam footprint size of the multibeam sonar and the side-scan sonar, the actual ground resolution of the two under the same water depth is calculated, and the ratio of the actual ground resolution is used as the resolution scaling factor λ. Using the resolution scaling factor λ as the initial scale parameter of the geometrically compensated side-scan image and 1 as the initial scale parameter of the multibeam image, the multibeam image and the geometrically compensated side-scan image are respectively subjected to progressive Gaussian blur and downsampling to obtain the multibeam Gaussian pyramid and the side-scan Gaussian pyramid. Perform differential operations on adjacent layers in the same group of the multibeam Gaussian pyramid to obtain a multibeam differential pyramid; perform differential operations on adjacent layers in the same group of the side-scan Gaussian pyramid to obtain a side-scan differential pyramid. By using the multibeam differential pyramid and the side-scan differential pyramid as their respective extreme point detection spaces, the multibeam scale space and the side-scan scale space are obtained.

4. The method for precise localization of ground features in multi-beam coordinate-assisted side-scan images according to claim 1, characterized in that, Step S3 involves extracting candidate feature points in both the multi-beam scale space and the side-scan scale space, including: Each pixel in the multi-beam scale space and the side-scan scale space is compared one by one with 26 neighboring pixels in the same group, the same layer and adjacent layers. Pixels whose gray values ​​are strictly greater than or strictly less than all 26 neighboring pixels are marked as initial extreme points, thus obtaining the multi-beam initial extreme point set and the side-scan initial extreme point set. For each point in the initial extreme point set of the multi-beam and the initial extreme point set of the side scan, the local surface of its scale space is fitted based on a three-dimensional quadratic function. The coordinates of the extreme points are corrected at the sub-pixel level according to the fitting offset. Points with contrast below a set threshold and edge response points with principal curvature ratio exceeding a set threshold are removed to obtain the precise feature point set of the multi-beam and the precise feature point set of the side scan. Taking each point in the multi-beam precise feature point set and the side-scan precise feature point set as the center, the pixel gradient direction histogram in its 3σ neighborhood is calculated. The peak direction of the histogram is taken as the main direction of the corresponding feature point, and the direction with a peak of more than 80% is taken as the auxiliary direction, thus obtaining the multi-beam feature point set with direction and the side-scan feature point set with direction. The directional multibeam feature point set and the directional side-scan feature point set are concatenated into an 8-directional gradient histogram within a 4×4 sub-region of each point to form a 128-dimensional description vector. The 128-dimensional description vector is then L2 normalized to obtain multibeam candidate feature points and side-scan candidate feature points, which together constitute candidate feature points.

5. The method for precise localization of ground features in multi-beam coordinate-assisted side-scan images according to claim 4, characterized in that, In step S3, confidence weights are assigned to the candidate feature points based on the depth validity of each point in the multibeam echo sounding grid, including: Based on the echo record status of each grid cell in the multibeam bathymetry grid, the grid cells are divided into three categories: measured effective cells, interpolated filling cells, and no-echo cells, thus obtaining the depth effectiveness classification grid. The planar coordinates of the multi-beam candidate feature points are projected onto the depth effectiveness classification grid. According to the category of the grid cell in which each feature point falls, the multi-beam candidate feature points are assigned values ​​one by one according to the rule that the measured effective cell is assigned a weight of 1.0, the interpolation filling cell is assigned a weight of 0.5, and the no-echo cell is assigned a weight of 0, so as to obtain the weight set of the multi-beam candidate feature points. The planar coordinates of the side-scan candidate feature points are projected onto the depth validity classification grid. Based on the grid distance between the grid cell where each feature point falls and the nearest measured valid cell, the weight value of each point is calculated inversely proportional to the distance. The side-scan candidate feature points are assigned values ​​one by one to obtain the side-scan candidate feature point weight set. The confidence weights of the candidate feature points are obtained by merging the weight set of the multi-beam candidate feature points with the weight set of the side-scan candidate feature points.

6. The method for precise localization of ground features in multi-beam coordinate-assisted side-scan images according to claim 5, characterized in that, The effective feature point set obtained by filtering using the confidence weight in step S3 includes: Based on the confidence weight, multi-beam candidate feature points with a weight value of 0 are directly removed from the candidate feature points, and side-scan candidate feature points with a weight value of less than 0.5 are directly removed from the candidate feature points, resulting in a weighted filtered feature point set. For the 128-dimensional description vectors of multibeam candidate feature points and side-scan candidate feature points in the weighted feature point set, calculate the Euclidean distance between each pair. Feature point pairs with a ratio of nearest neighbor Euclidean distance to second nearest neighbor Euclidean distance of less than 0.6 are marked as initial matching pairs to obtain the initial matching feature point set. Based on the confidence weight, the Euclidean distance ratio of each matching pair in the initial matching feature point set is adjusted by weight. Matching pairs whose adjusted Euclidean distance ratio is still lower than 0.6 are retained to obtain the weight-adjusted matching feature point set. The weight-adjusted matching feature point set is used as the effective feature point set.

7. The method for precise localization of ground features in multi-beam coordinate-assisted side-scan images according to claim 1, characterized in that, Step S4 includes: Based on the confidence weight of each point in the effective feature point set, feature points with a weight value of not less than 1.0 are assigned to the priority sampling point set, and feature points with a weight value less than 1.0 are assigned to the candidate sampling point set. The priority sampling point set is used as the first sampling source, and the candidate sampling point set is used as the supplementary sampling source to obtain the weighted sampling sequence. In each iteration of the RANSAC algorithm, four pairs of image points consisting of multi-beam candidate feature points and side-scan candidate feature points are extracted from the effective feature point set according to the weighted sampling sequence. The candidate homography matrix is ​​solved based on the pair of image points. The projection error of each point to the candidate homography matrix is ​​calculated by traversing all point pairs in the effective feature point set. Point pairs with projection errors lower than a set distance threshold are marked as inliers. The number of inliers is counted. The candidate homography matrix with the largest number of inliers is updated as the current best homography matrix. The iteration is repeated for a set number of times K to obtain the homography matrix H. Substitute the geodetic coordinates of the multibeam candidate feature points in the effective feature point set into the homography matrix H, calculate the corresponding pixel coordinates of each point in the geometrically compensated side-scan image, and obtain the coordinate mapping result; The planar residuals between each mapping point and its corresponding side-scan candidate feature point in the coordinate mapping result are statistically analyzed. Mapping points with residuals exceeding a set threshold are removed, and the geodetic coordinates of the remaining mapping points are used as the accurate positioning result.

8. A system for precise ground feature localization based on multi-beam coordinate-assisted side-scan images, characterized in that, A method for accurately locating ground features in a multi-beam coordinate-assisted side-scan image as described in any one of claims 1-7, wherein the accurate ground feature location system for the multi-beam coordinate-assisted side-scan image comprises: The conversion module is used to convert the coordinates of each pixel in the original side-scan image from the slant range coordinate system to the ground coordinate system, based on the drag depth and slant range records of the side-scan sonar, to obtain a geometrically compensated side-scan image. The construction module is used to calculate the resolution scaling factor λ based on the ratio of the actual ground resolution of the multibeam sonar to that of the side-scan sonar. Using λ as the initial scale parameter, Gaussian difference pyramids are constructed for the multibeam image and the geometrically compensated side-scan image respectively, so as to obtain the multibeam scale space and the side-scan scale space. The filtering module is used to extract candidate feature points in the multibeam scale space and the side scan scale space respectively, assign confidence weights to the candidate feature points according to the depth validity of each point of the multibeam sounding grid, and obtain a set of effective feature points by filtering with the confidence weights. The mapping module is used to prioritize sampling objects with high-weight points in the effective feature point set, solve the homography matrix H between the multibeam image and the geometrically compensated side scan image using the RANSAC algorithm, and map the multibeam geographic coordinates to the side scan image via the homography matrix H to obtain accurate positioning results.

9. A device for precise ground feature localization using multi-beam coordinate-assisted side-scan images, characterized in that, The method includes a memory and a processor, the memory storing a computer program that can run on the processor, and the processor executing the computer program to implement the method for precise localization of ground features in multibeam coordinate-assisted side-scan images according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is run by the processor, it causes the processor to perform a method for precise localization of ground features in a multi-beam coordinate-assisted side-scan image as described in any one of claims 1 to 7.

Images