Obtaining High-Resolution Seabed Topography and Surface Details by Co-Registration of Side-Scan Sonar and Multibeam Echo Sounder Images

Abstract

Side-scan sonar (SSS) is used for obtaining high-resolution seabed images, but with low position accuracy without using the Ultra Short Base Line (USBL) or Short Base Line (SBL). Multibeam echo sounder (MBES), which can simultaneously obtain high-accuracy seabed topography as well as seabed images with low resolution in deep water. Based on the complementarity of SSS and MBES data, this paper proposes a new method for acquiring high-resolution seabed topography and surface details that are difficult to obtain using MBES or SSS alone. Firstly, according to the common seabed features presented in both images, the Speeded-Up Robust Features (SURF) algorithm, with the constraint of image geographic coordinates, is adopted for initial image matching. Secondly, to further improve the matching performance, a template matching strategy using the dense local self-similarity (DLSS) descriptor is adopted according to the self-similarities within these two images. Next, the random sample consensus (RANSAC) algorithm is used for removing the mismatches and the SSS backscatter image geographic coordinates are rectified by the transformation model established based on the correct matched points. Finally, the superimposition of this rectified SSS backscatter image on MBES seabed topography is performed and the high-resolution and high-accuracy seabed topography and surface details can be obtained.

1. Introduction

Seabed topography and morphology provide fundamental information for marine scientific research and ocean engineering [1,2,3,4]. Currently, the most widely used sensors in conducting ocean surveying and mapping are side-scan sonar (SSS) and multibeam echo sounder (MBES) [5,6].

As for the SSS, there are two types of operations: towing operation and hull-mounted [5]. The first type is more popular given that the second type of operation may degrade the SSS backscatter image quality because of the platform movements caused in the surveying process [5,7]. When towing, the SSS is usually installed in a towfish and towed by a cable behind a surveying vessel (Figure 1a). It emits a wide-angle beam and receives thousands of seabed echoes at fixed time intervals to form a high-resolution seabed image [5,8]. This kind of operation can minimize the effects of platform movements on the SSS backscatter images. More importantly, as the towfish is near the seabed, it can be used for deep seafloor observation [5]. Thus, the SSS backscatter image geographic coordinates are not accurate because of towing, potential currents dragging and flat bottom assumption [9,10,11,12]. By using the Ultra Short Base Line (USBL) or Short Base Line (SBL) in the surveying process, the accuracy of towfish positions can be improved. However, these two devices are not used in most cases because of high cost. Moreover, three-dimensional (3D) topography is not available by adopting SSS [5,13]. Although the emergence of interferometric SSS (ISSS) can acquire bathymetry, the measurement quality is very sensitive to underwater noise and seriously degrades for echoes close to nadir [14]. In addition to the ISSS, there also exist some other kinds of sonar systems that can acquire bathymetry data, such as Ping DSP 3DSS, Edgetech 6205 and Klein HydroChart 3500 [15], but they can only be used for hydrographic surveys in restricted areas, such as roadsteads or lakes [15]. The output from 3DSS is in the form of an intensity point cloud, which can be post-processed either as bathymetry, backscatter seabed image or a combination of both. However, this kind of data has been found less capable in the detection of seabed targets [16]. As for the Edgetech 6205 or Klein HydroChart 3500, they are mainly used in shallow waters, less than 35 m or 50 m respectively [17,18]. These new sonar systems have restrictions and their performances in the field of bathymetry measurement are usually inferior to those obtained by MBES [6].

MBES is designed for high-accuracy bathymetric work and modern MBES can receive hundreds of echoes from the seabed for one ping and provide a highly detailed backscatter image simultaneously. The MBES is usually installed on a surveying vessel (Figure 1b). Although it can be used for large-scale bathymetric measurement, the interval of footprints in a ping will be enlarged and the bathymetric resolution will be decreased as the water depth and beam incident angle grow. When producing a MBES image, backscatter data are extracted from the time series traces contained within a beam are recorded [19,20]. The central point of a beam is optimally positioned with minimal geometric distortions and the image pixels are distributed around it. With seabed topography, vessel attitudes and the refraction of acoustic waves considered, the positions of MBES image pixels can be accurately computed [19]. By superimposing the MBES image on the seabed topography, the comprehensive seabed topography and surface features can be presented together [21,22,23]. However, the resolution, signal to noise ratio (SNR) and intensity contrast of MBES images are lower than those of SSS backscatter images [6]. The complementarity of SSS and MBES data provides a way to obtain detailed seabed features by superposing a two-dimensional (2D) SSS backscatter image onto 3D MBES bathymetric terrain.

However, because SSS backscatter image geographic coordinates are inaccurate, it is a challenge to conduct the superimposition of these two categories of data directly. Much research has been carried out in this field. LeBas et al. [24] created synthetic image from the seabed topography and adopted the hierarchical chamfer matching method to match the SSS backscatter and synthetic images. According to the matched result, the SSS backscatter image geographic coordinates were rectified and the obtained SSS backscatter image was superposed on the seabed topography. Yang et al. [25] used the similarity between MBES topography isobaths and SSS backscatter image contours to carry out image matching and the following superimposition operation. The two methods need sufficient topological variations on the seafloor for creating distinct edges or feature points in sonar images and may fail when dealing with a flat seabed with various sediments. According to similar seabed topography characteristics, textures, targets and sediment distributions being reflected on both SSS and MBES images, Zhao et al. [26] adopted the improved Speeded-Up Robust Features (SURF) algorithm to extract common feature points and conducted image registration, then rectified the SSS backscatter image geographic coordinates, referring to the MBES image based on the matched results. To obtain more correct matches, they conducted seabed classification for SSS and MBES images and used the classification images for image matching followed by a superimposition operation [27]. This method can improve matching performance but its results depend heavily on the seabed classification accuracy. The above methods ignore the imaging mechanism differences between SSS and MBES images, which often lead to inaccurate descriptions of the feature points and result in inaccurate image matching. Considering the limitations of the existing methods, this paper proposes a new SSS and MBES image matching method for acquiring high-resolution and high-accuracy seabed topography and surface details, which can overcome the limitations of adopting a single MBES or SSS for seabed mapping. The remainder of this paper is organized as follows: Section 2 describes the superimposition method in detail; Section 3 designs the experiments to verify the proposed method; Section 4 analyzes the results; Section 5 discusses the proposed method; and Section 6 presents the conclusions according to the experiments and discussions.

2. Method

To get high-resolution seabed topography and surface backscatter data, the proposed method is shown in Figure 2.

Firstly, both SSS and MBES data are processed to form the geocoded images.

Secondly, the SSS and MBES image matching based on common features is conducted, which includes the initial matching with the SURF algorithm with image geographic coordinate constraint and finer matching by using dense local self-similarity (DLSS) descriptors.

Thirdly, the random sample consensus (RANSAC) algorithm is used for removing mismatches and then the SSS backscatter image geographic coordinates are rectified according to the relationship model established by the correct matched points.

Finally, according to the consistent geographic coordinates, the rectified SSS backscatter image is superposed on the MBES bathymetric terrain.

The proposed method is described in detail in the following paragraphs and the image matching program was written using Matlab in a computer with the i7, 3.40GHz Intel Core and 8.00 GB RAM.

2.1. SSS and MBES Data Processing

The original SSS data were first decoded to form waterfall images ping by ping. In a complete process of SSS data, bottom tracking, radiometric correction, geometric correction and geocoding,and so forth, are included [5,28,29,30,31,32]. This SSS data processing was conducted by self-developed software and the processing procedure is shown in Figure 3a. Because thousands of echoes exist in a ping scanning line, a high-resolution SSS backscatter image could be obtained after the above data processing. By referring to Figure 1a, the positions of SSS towfish were reckoned by combining vessel-mounted GPS positioning results, cable length and the towfish heading, and thus became inaccurate due to the towing operation, the flat bottom assumption, towing distance and bearing controlled by the towing speed and the currents. Therefore, the geographic coordinates of high-resolution SSS backscatter image needed to be rectified.

The original MBES data were first decoded to get the bathymetric data, backscatter data, GPS data, attitude data, sound velocity, and so forth. Quality control for these raw data was done first. Then, calculation of bathymetric point was conducted which includes attitude correction, sound ray tracing, absolute coordinates obtaining of bathymetric point, and so forth. After getting the locations of beam footprints, the authors conducted the radiometric correction, angular response correction and geocoding for the MBES backscatter data, and could then obtain the seabed image. Finally, after tidal correction for the sound ray tracing results, the absolute 3D coordinates of each bathymetric point were obtained and thus the seabed topography could be achieved. The procedure of MBES backscatter and bathymetric data processing was displayed in Figure 3b. To avoid gaps on the MBES image, the survey vessel velocity should be determined appropriately. According to the MBES data processing procedure and MBES measurement principle (Figure 1b), the geographic coordinates of bathymetric and backscatter data were accurate. However, the seabed topography and image obtained by using MBES were of low resolution, as only hundreds of echoes exist in one ping. As for the MBES backscatter, there were three types of acquired data: single beam intensity, individual beam time series and integrated time series. Although there exist a few to dozens of snippets in each beam footprint, the resolution of MBES image was still far lower than that of the SSS backscatter image. Thus, the MBES image is usually used for reflecting the large-scale seabed surface features instead of high-resolution visualization.

2.2. SSS and MBES Image Matching

Though with different characteristics (e. g. resolution, position accuracy and intensity contrast), both SSS and MBES images can reflect seabed features of the same area. Thus, the SSS and MBES image matching can be conducted by searching for common features in both images. To get accurate matched points, the matching method is proposed, which includes initial matching using the SURF algorithm with image geographic coordinate constraint, and the finer matching by using dense local self-similarity (DLSS) descriptors. The initial image matching step is to detect image feature points in the first place and use the SURF descriptors for feature-based image matching with the constraint of image coordinates. Since the SURF algorithm is not robust to nonlinear intensity differences between SSS and MBES images, some mismatches may arise in the initial matching step. To obtain more accurate matched points, the DLSS descriptor was used in the finer matching step and a template matching strategy was adopted. This two-step image matching process is described in detail in the following paragraphs.

2.2.1. Initial Image Matching Using SURF Algorithm with Image Geographic Coordinate Constraint

The SURF algorithm, which is invariant to rotation, scale, brightness and contrast [33], was adopted for detecting candidate feature points in SSS and MBES images to produce the SURF descriptors. By concatenating a 4D descriptor vector V of each sub-region, which underlies intensity structure for the square region of 4 × 4, the SURF descriptor of length 64 can be obtained [33].

V=(∑_dx,∑_dy,∑|_dx|,∑|_dy|)

where, d_x and d_y are respectively the Haar wavelet responses in horizontal and vertical directions. Since using the SURF descriptors for multi-model image matching directly may lead to many mismatches, some valid constraints can be added for specific applications [34,35].

For geocoded SSS and MBES images, the accuracies of their geographic coordinates are different and the geographic coordinates of the MBES image are more accurate than those of SSS because the SSS surveying is easily affected by the towing operation and many other factors. This inconsistency of image positional accuracy leads to the fact that the geographic coordinates of the same features presented in both images will be different. This kind of difference is within a range and can be determined by SSS location error, which is discussed in detail in Section 5.1. Thus, when the distances of the detected feature points in both images are too large, they will certainly not be the correct matches. Accordingly, the SURF algorithm for SSS and MBES image matching can be conducted when the distances of detected feature points in both images are within a threshold Dis, which can reduce the mismatches compared with using the SURF algorithm directly. With this constraint, the initial feature-based image matching was conducted as shown in Equation (2).

if_Di≤Disotherwise then SURF matchingRefused

where, D_i is the distance between the candidate matched points. The determination of threshold Dis is also described in detail in Section 5.1.

2.2.2. Finer Image Matching Using DLSS Descriptors

Because the SSS backscatter reflects relative backscatter levels and the MBES images are usually normalized to dB scale, there is no linear relationship between their backscatter intensity levels. Meanwhile, as the resolution, SNR and intensity contrast of SSS backscatter image are higher than those of the MBES image, the same feature points of the seabed may be presented more clearly in SSS backscatter image than in the MBES image. As a result, some feature points can be detected in the SSS backscatter image but not in the MBES image. The SURF algorithm may detect incorrect matches when dealing with the non-linear intensity differences because it mainly focuses on the intensity variations of feature points and spatial distributions of their gradient information [33,34,35,36]. Moreover, as the number of detected feature points in the two images were different, the finally obtained number of matched points was limited. To get better matched points, the shape properties of distinct targets, topography changes or sediment distributions in these images can be used.

The DLSS descriptor can represent these shape properties and is robust against significant nonlinear intensity differences by integrating the local self-similarity (LSS) descriptors of multiple local image regions in a dense sampling grid [37]. LSS was calculated by correlating the image patch centered at a point q with a larger surrounding image region, resulting in a correction surface S_q(x,y) as below [38].

_Sq(x,y)=exp(−SS_Dq(x,y)max(_varnoise,_varauto(q))

where, SSD_q(x,y) is the sum of square differences between image patch intensities, var_noise is a constant corresponding to intensity variations and var_auto is the maximal SSD of all patches within a very small neighborhood of q. Then, the S_q(x,y) was transformed into a log-polar representation and partitioned into bins. The maximal correction value of each bin was selected to generate the LSS descriptor associated with the pixel q. In order to be invariant to the intensity variations of the image, the LSS was linearly stretched into the range [0, 1] [38]. An example of the LSS descriptor is shown in Figure 4, which indicates that the LSS descriptors of feature points are similar when the points share a common layout.

By collecting the LSS descriptors from spatial regions within a dense overlapping grid covering a defined square area, the DLSS was produced as a combined feature descriptor [37]. Using the DLSS descriptors, a template matching strategy was conducted based on the initial image matching results.

First, the DLSS descriptors of initial matched points in the SSS backscatter image were produced and a searching area centered at the initial matched points in the MBES image with the radius r was defined. Meanwhile, the DLSS descriptors of every pixel in this defined area were calculated.

Secondly, calculating the similarity between the DLSS descriptors of the SSS backscatter image feature points and those of the points still unprocessed in the above defined searching area in the MBES image. The normalized cross correlation (NCC) of DLSS descriptors as shown in Equation (4) was adopted as the similarity metric.

NCC=∑i=1n(DLSS_1i−mean(DLSS1))(DLSS_2i−mean(DLSS2))∑i=1n^{(DLSS_1i−mean(DLSS1))2}∑i=1n^{(DLSS_2i−mean(DLSS2))2}

where, DLSS1 and DLSS2 are the DLSS descriptors of points in SSS and MBES images, N is the dimension of DLSS descriptor.

Finally, by choosing maximum NCCs and selecting the corresponding points in the defined area of the MBES image, the finer matches could be obtained. After conducting the above process for all the initial matched points in the SSS backscatter and MBES images, more correct matched points were obtained.

2.3. SSS Backscatter Image Geographic Coordinates Rectification Based on the Matched Points

Even after the finer matching step, mismatches still existed. The reason is that, affected by the measurement environment (e.g., the currents, the surveying vessel velocity and attitude changes), the seabed features presented in the SSS and MEBS images may not be identical. To eliminate these mismatches, the random sample consensus (RANSAC) algorithm was adopted. Since the RANSAC algorithm was conducted by selecting some existing matched points to build a geometric model [39] and the local geometric distortions in different parts of the SSS backscatter image were different, the SSS backscatter image could be segmented into different blocks according to feature distributions and towfish attitudes and the RANSAC algorithm could be conducted separately in each block [40,41,42]. When segmenting the SSS backscatter image, enough matched points should be contained in each block and the split lines are supposed to lie on the edge of the targets [41,42]. Meanwhile, to ensure the continuity between adjacent blocks, an overlapping ratio (e.g., one seventh of the block) can be set between them. To further ensure the remaining matches were correct, a geometric transformation model was established by using the matched points after application of the RANSAC algorithm and the coordinates of feature points in the SSS backscatter image were rectified. By referring to the matched feature points in the MBES image, the corresponding feature points in the rectified SSS backscatter image were considered as correct matches when their coordinate errors were less than double standard deviations. After this process, the correct matched points were obtained.

Based on the correct matched points between SSS and MBES images in each block, a geometric transformation model could be built to rectify the geographic coordinates of the SSS backscatter image by referring to those of the MBES image. Since the geographic distortions in local regions of the SSS backscatter image can be modeled by the affine function [43], the affine transformation model as shown in Equation (5) can be built within each block.

_xm=_a0+_{a1 xs}+_{a2 ysym}=_b0+_{b1 xs}+_{b2 ys}

where, (x_s,y_s) and (x_m,y_m) are geographic coordinates of matched feature points in SSS and MBES images, a₀, a₁, a₂ and b₀,b₁,b₂ are the model parameters to be calculated. After obtaining the geometric transformation models within each segmented block, the SSS backscatter image geographic coordinates could be rectified.

2.4. Superimposition of Rectified SSS Backscatter Image on MBES Terrain

After the rectification, the geographic coordinates of the SSS backscatter image were inconsistent with those of MBES image. Thus, the rectified SSS backscatter image could be superposed on the MBES terrain. This process was conducted by following steps:

(1) Based on the matched points, the geographic coordinates of the SSS backscatter image were rectified by referring to those of the MBES image.

(2) The 3D seabed topography was obtained by using the regular square grid method, which has the same resolution as the SSS backscatter image.

(3) According to the geographic coordinates, each grid point of the 3D topography was attributed with the corresponding grey value of SSS backscatter image.

After this process, the high-resolution 3D seabed topography and surface details can be presented together.

3. Experiments and Results 3.1. Experimental Data

To validate the proposed method, an experiment was carried out in Meizhou Bay, China with the area of 1.4 × 0.3 km. In this water area, the depth ranges from 10 to 14 m; seabed sampling showed sediments mainly contain silt and sand. In this experiment, EdgeTech 4100P with the slant range of 112 m and 7502 sampling points for each ping was adopted for the SSS measurement. For the MBES measurement, EM 3002 with a maximum angular sector of 130° and 131 beams was adopted. The MBES data were collected using the Seabed Information System (SIS) and restored in “*.all” files. The SSS data were collected using the EdgeTech Discover Software and restored in “*.xtf” files. After collecting these data, both SSS and MBES data were processed by self-developed software. The processing of SSS and MBES data is shown in Figure 3. After data processing, the geocoded SSS backscatter image, MBES image and seabed topography were created, as shown in Figure 5a–c. The coordinates are in WGS84. Considering the sensitivity of geographical data, the coordinates have been subtracted by a constant. Meanwhile, the MBES and SSS backscatter images and the MBES topography were resampled in the same ground sample distance (GSD) of 1m. Even with the same GSD, the SSS backscatter image could reflect clearer contours of seabed surface features because there exist many more sampling points for one ping of the SSS measurement than that of the MBES measurement. For example, the pipeline is presented more clearly in the SSS backscatter image than in the MBES image. Meanwhile, the topography shown in Figure 5c mainly reflects the 3D seabed undulations but with less detail.

Thus, the superimposition of the high-resolution SSS backscatter image on high-accuracy seabed terrain can take advantage of both datasets, which can capture both the 3D seabed topography and the detailed surface features together.

3.2. Experimental Procedure and Results

According to the proposed method depicted in Section 2, the SURF algorithm was first used for detecting feature points in SSS and MBES images. In a unified geographical framework, the image geographic coordinates can serve as the constraint in the initial image matching step as described in Section 2.2.1. This matching result is shown in Figure 6b. Compared with the matching results shown in Figure 6a, which were obtained by the classical SURF algorithm without constraint, the initial matching result with constraint seems better. The reason is that the SURF algorithm is not robust when dealing with the nonlinear intensity difference between SSS and MBES images. Meanwhile, the features existed on the seabed are relatively simple and the detected edge and point features can be represented by similar SURF descriptors, which may decrease the distinctiveness of SURF descriptors. Because of the two factors, even when two feature points in the SSS and MBES images are distant, they may be regarded as matches when they are represented by similar SURF descriptors. When using the image geographic coordinates as constraint in the initial matching step, the initial image matching is only conducted within a valid distance, which avoids mismatches when the geographic coordinates of two feature points are too far away.

Moreover, since there are many more sampling points for SSS measurement than for MBES measurement, the seabed features in the MBES image may not be as clear as those in the SSS backscatter image. As a result, some feature points may be detected in the SSS backscatter image but cannot be presented in the MBES image. Consequently, it was found that several detected feature points in the SSS backscatter image are matched to one feature point in the MBES image, particularly in zoomed area (b1) of Figure 6b. This several-to-one matching problem is obviously incorrect, which will be settled in the finer image matching step.

As image geographic coordinates are used as a constraint in the initial image matching step, the correct matched points in the MBES image are supposed to be located within an area centered at the initial matched ones. Thus, the optimization search operation described in Section 2.2.2 was conducted and the better matched points are shown in Figure 6c. In this process, a template matching strategy using DLSS descriptors was conducted. Compared with the SURF descriptor, the DLSS descriptor is robust to the nonlinear intensity difference between SSS and MBES images because it reflects the shape properties of feature points and their surrounding areas. Meanwhile, this finer image matching using a template matching strategy can also help find the same number of feature points in the MBES image as these in the SSS backscatter image. As a result, more correct matched points can be obtained by adopting the finer image matching step. Compared with Figure 6b, obtained from the initial image matching step, the several-to-one matching phenomena was removed after the finer matching step, which is clearly shown in the zoomed area (c1) of Figure 6c. This improvement of matching performance is also shown in Table 1. Table 1 lists the numbers of all the matches obtained by using SURF, initial matching and finer matching and that of the correct matches after the application of the RANSAC algorithm.

Even if the finer image matching step is followed, there still exist mismatches in Figure 6c. To eliminate theses mismatches, the SSS backscatter image was segmented into blocks and the RANSAC algorithm was adopted in each block. In the segment process, it can be seen that there exist three main independent areas with abundant seabed features. Thus, the SSS backscatter image was segmented into three blocks according to the feature distributions and towfish attitudes. Meanwhile, the overlapping ratio of one seventh between adjacent blocks was set and the split lines lay on the edge of isolated targets, which can be seen in the red rectangles in Figure 6 and Figure 7. The image matching results after this process are shown in Figure 7 and Table 1. It can be seen that even after using the RANSAC algorithm, some mismatches still exist in Figure 7a, because too many mismatches occurred due to the use of the SURF algorithm directly and the RANSAC algorithm is sensitive to the high outlier rate. As a result, these mismatches cannot be effectively eliminated by using the RANSAC algorithm. This phenomenon also proves that adding the constraint of image geographic coordinates in the initial matching step was necessary, which improved the matching ratio from 8% to 55% as shown in Table 1. Even so, mismatches still exist after the initial step, which can be seen in the zoomed area (b1) of Figure 7b. After the finer image matching step, the matching ratio rose to 86%. More importantly, these matches are correct, as shown in the zoomed area (c1) of Figure 7c. This is because the combination of both the initial and finer image matching not only makes full use of the intensity variations and gradient information of feature points, but also considers their geographic coordinates and the shape properties of the area around them.

Using the correct matched points that were finally obtained in each block, the transformation model shown in Equation (5) could be established. According to the obtained transformation model, the SSS backscatter image geographic coordinates could be rectified following the method in Section 2.3. The accuracy of the rectified SSS backscatter image geographic coordinates is discussed in Section 4.

Superimposing the rectified SSS backscatter image on seabed topography, the comprehensive presentation of 3D topography and surface details is shown in Figure 8. Compared with Figure 5c, Figure 8 can reflect the topography undulations in combination with highly detailed SSS backscatter imagery.

4. Evaluation and Analysis

To assess the accuracy of the rectified SSS backscatter image geographic coordinates, the interior checking and exterior checking were conducted. The interior checking involves calculating the geographic coordinate biases of matched feature points used to build the transformation model before and after the rectification, while the exterior checking involves calculating the matches not used to build the transformation model. Meanwhile, the statistical results of these biases were also calculated. In this experiment, two-thirds of matched points were used for establishing the transformation model and the interior checking; the remaining one-third of matched points were consequently used for the exterior checking. The statistical results and probability distribution function (PDF) curves of each checking are shown in Table 2, Table 3, Figure 9 and Figure 10. In Table 2 and Table 3, “Raw” means the coordinate biases of matched feature points in the raw SSS and MBES images; “Rectified” means those in the rectified SSS and MBES images; “dy” means those of the east direction, and “dx” means those of the north direction.

From Table 2, Table 3, Figure 9 and Figure 10, it can be seen that obvious systemic errors and large standard deviations exist between the raw geographic coordinates of matched points. Table 2 and Table 3 show the systemic errors in the two directions are 4.32 m and 5.98 m for interior checking and 5.06 m and 6.17 m for exterior checking. After using the transformation model, these systemic errors were removed. Meanwhile, the standard deviations became smaller, which can be seen in Table 2, Table 3, Figure 9 and Figure 10. These results prove that the accuracy of the SSS backscatter image geographic coordinates improved after the rectification. Moreover, the accuracy of both checkings were consistent with each other, which indicates that the proposed method is robust in building the transformation model.

5. Discussions 5.1. Determination of the Threshold Dis and Searching Radius

In the initial image matching process, the image geographic coordinates served as the constraint and the threshold Dis were determined by SSS location error. Usually, SSS is installed in a towfish and towed by a cable behind a surveying vessel for seabed surveying, which is shown in Figure 1a. During this process, the towfish positions can be estimated by Equation (6).

_{xtowfishytowfish}=_{xvesselyvessel}+L·cosα·cos(A+π)L·cosα·sin(A+π)

where, (x_towfish, y_towfish) and (x_vessel, y_vessel) are respectively the coordinates of the towfish and vessel tow point, L is the cable length, α is the angle between the cable and horizontal direction, A is the vessel heading. Affected by the vessel velocity variation, the accuracies of vessel positions and heading data, the towfish position error d_towfish can be determined by Equation (7).

_dtowfish=d_xtowfish ²+d_ytowfish ²=d_xvessel ²+d_yvessel ²+^L2·^sin2α·d^α2+^L2·^cos2α·d^A2

where, (dx_towfish, dy_towfish) and (dx_vessel, dy_vessel) are respectively the towfish and tow point position errors in x and y directions, dα is the variation range of the angle α, dA is the heading difference between the vessel and towfish. For most SSS measurements in offshore waters, the vessel position error is about 1~3m; the length of the cable is about 10~20m; the angle α is about 30° and its variation range is about 10°; the heading difference between the vessel and towfish is about 1°–5° [44]. After substituting these values into Equation (7), the d_towfish ranges roughly from 10 to 100 m. The above parameters will become bigger with the increase of surveying vessel size in open sea. In this experiment, the sea was calm and the vessel’s course was steady. Thus, the d_towfish was about 18 m, which served as a referenced threshold in the initial image matching. To determine the optimal threshold, different thresholds in Equation (2) were used by referring to d_towfish and the matched results are shown in Table 4, which shows that the detected matched points increased as the threshold became larger but the number of correct ones peaked when the threshold was 20. Thus, the initial matched result was obtained with this threshold and followed by the finer matching.

In the finer matching process, the template matching strategy was used and the searching radius for better matches was determined. A small searching radius is not enough for finding better matches while a large one will cost more time and the matched results may also not be the optimal. During the finer matching process, different searching radiuses were used and the matched results are shown in Table 5. This indicated that 20 m could be adopted as a suitable searching radius in the finer matching step as most correct matches were obtained when using this value.

5.2. Method Repeatability: Applicability of the Proposed Method in another Situation

To further verify the performance of the proposed method in other waters when adopting different kinds of sonars, another experiment was conducted in BeiBu Gulf. In this area, the depth ranges from 40 to 50 m; seabed sampling shows that sediments mainly consist of sand, which has two forms of existence: sand waves and ripples. Klein 3000 and EM710 were separately adopted for SSS and MBES measurement. After processing the obtained measurement data—the seabed topography is shown in Figure 11—the SSS and MBES images were presented in Figure 12. It can be seen that the seabed topography and the MBES image can only reflect the sand waves but the SSS backscatter image can also capture the small scale sand ripples clearly. This difference can be seen in areas (a1) and (b1) of Figure 11. This difference can also be seen in zoomed areas (a1) and (b1) of Figure 12.

Based on the common sand wave features reflected by both images, the proposed image matching method described in Section 2.2 was performed and the matched result shown in Figure 12 proved that the method can also be effective when using other kinds of measurement instruments in a different water area. An area with sharper edges is zoomed in for better visualization, as shown in Figure 12c. Using the correct matches that were finally obtained to rectify the SSS backscatter image geographic coordinates, the obtained SSS backscatter image was superposed on the seabed topography. The comprehensive presentation of 3D topography and surface details is shown in Figure 13. Compared with Figure 11, it can reflect both the obvious sand waves and detailed sand ripples.

5.3. Sonar Frequencies

The acoustic waves emitted by a sonar propagate through sea water to seabed and the backscattered echoes are received by the transducer to form seabed backscatter images. When using different frequencies of acoustic wave, the impacts of surface and volume scatter are different [45,46,47]. Taking the sandy sediments as an example, penetration depth is limited to 12mm for 500 kHz, while 100 kHz may penetrate 90mm into the subsurface [48]. In addition, the acoustic frequency affects sonar backscatter image texture presentations due to changing surveying parameters, such as footprint sizes and different sensitivity to seafloor features [23,45]. As a result, the presentation of seabed features may be different when using different SSS and MBES frequencies [23,45,47]. As a result, finding matches might become difficult.

At specific frequencies, the sonar backscatter images can always reflect seabed features and the relative positions of these features are stable. Therefore, based on the geometrical structure formed by the association of three detected features, image matching can still be conducted [43]. Meanwhile, deep learning has also been applied in the sonar image processing field [49,50] and some related image matching methods [51] can be used for sonar image matching.

Recently, multi-frequency SSS and MBES have been applied, such as EM2040 and Edgetech 4200SP. As for multi-frequency sonar backscatter data, they can be rendered as acoustically colorful images [23,46]. The formed colorful image can reflect more seabed features compared with those obtained from monochromatic sonar with a single center frequency. Consequently, more features can be detected from the colorful image and can be used for image matching.

6. Conclusions

Despite some new sonar systems appearing that can provide co-registered seabed images and topography, they can only be used in restricted waters. The SSS and MBES are still the most effective sonars for obtaining seabed images and topography, especially for deep seafloor observation. Based on the complementarity of SSS and MBES data, this paper proposes a new method for acquiring high-resolution seabed topography and surface details by combining SSS and MBES data. Through taking the image geographic coordinates as the constraint when using the SURF algorithm for initial image matching, the authors have obtained more correct initial matched points compared to that without constraint. Then, the finer matching step is conducted by adopting a template matching strategy, which uses the DLSS descriptor to reflect the shape properties of the area centered feature points. The combination of the initial and finer matching steps can help improve the matching performance and the percentage of correct matched points rose to 86%. Based on the obtained matched points and considering the difference of geographic distortions in different parts of the SSS backscatter image, the SSS backscatter image was segmented into several blocks and the geometric transformation model within each block can be established to rectify its geographic coordinates. Subsequently, the rectified SSS backscatter image can be superposed on the topography, which can reflect both the 3D seabed undulations in combination with a high detailed SSS backscatter imagery.

Experiments have verified this method. By this method, the high-resolution and high-accuracy seabed topography and surface details can be represented together, which is meaningful for understanding and interpreting seabed topography. Meanwhile, this paper discusses the accuracy of the reckoned SSS positions and uses it as a reference threshold in the image matching process. In addition, this paper discusses the impact of sonar frequency on the sonar backscatter image and provides some useful suggestions when dealing with multi-frequency sonar image matching.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

	Matches	Correct Matches	Percentage
SURF	341	28	8%
Initial matching	231	128	55%
Fine matching	231	198	86%

	Bias	Max	Min	Mean	Std.
Raw	dy	7.79	−19.88	−4.32	3.67
	dx	8.57	−17.87	−5.98	4.45
Rectified	dy	9.61	−7.69	0.01	2.95
	dx	9.43	−7.35	0.00	3.17

	Bias	Max	Min	Mean	Std.
Raw	dy	1.68	−19.93	−5.06	3.41
	dx	8.47	−15.79	−6.17	4.83
Rectified	dy	7.97	−7.12	−0.48	2.58
	dx	9.88	−7.28	−0.27	3.53

Threshold (m)	Detected Matches	Correct matches	Percentage
10	109	92	84%
20	231	128	55%
30	285	120	42%
40	317	102	32%
50	327	100	30%

Searching Radius (m)	Detected Matches	Correct Matches	Percentage
10	231	167	72%
20	231	198	86%
30	231	174	75%
40	231	162	70%
50	231	146	63%

Author Contributions

Conceptualization, X.S. and J.Z.; Data curation, X.S. and J.Z.; Formal analysis, X.S., J.Z. and H.Z.; Funding acquisition, J.Z. and H.Z.; Methodology, X.S. and J.Z.; Supervision, J.Z. and H.Z.; Visualization, X.S., J.Z. and H.Z.; Writing—original draft, X.S.; Writing—review & editing, X.S., J.Z. and H.Z.

Funding

This research was funded by the National Natural Science Foundation of China, grant numbers 41576107, 41376109 and 41606114.

Acknowledgments

The authors would like to thank the editors and anonymous reviewers for their comments and suggestions. The data used in this study were provided by the Guangzhou Marine Geological Survey Bureau. The authors are grateful for their support.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Mayer, L.; Jakobsson, M.; Allen, G.; Dorschel, B.; Falconer, R.; Ferrini, V.; Lamarche, G.; Snaith, H.; Weatherall, P. The Nippon Foundation—GEBCO Seabed 2030 Project: The Quest to See the World’s Oceans Completely Mapped by 2030. Geosciences 2018, 8, 63.

2. DeSanto, J.B.; Sandwell, D.T.; Chadwell, C.D. Seafloor geodesy from repeated sidescan sonar surveys. J. Geophys. Res. Solid Earth 2016, 121, 4800–4813.

3. Powers, J.; Brewer, S.K.; Long, J.M.; Campbell, T. Evaluating the use of side-scan sonar for detecting freshwater mussel beds in turbid river environments. Hydrobiologia 2015, 743, 127–137.

4. Lurton, X. Forty years of progress in multibeam echosounder technology for ocean investigation. J. Acoust. Soc. Am. 2017, 141, 3948.

5. Blondel, P. The Handbook of Sidescan Sonar; Springer: Berlin/Heidelberg, Germany; New York, NY, USA, 2009; ISBN 978-3-540-42641-7.

6. Lurton, X.; Jackson, D. An Introduction to Underwater Acoustics, 2nd ed.; Springer-Praxis: New York, NY, USA, 2008; ISBN 3540429670.

7. Tamsett, D.; Hogarth, P. Sidescan sonar beam function and seabed backscatter functions from trace amplitude and vehicle roll data. IEEE J. Ocean. Eng. 2015, 411, 155–163.

8. Bell, J.M.; Linnett, L.M. Simulation and analysis of synthetic sidescan sonar images. IEE Proc. Radar Sonar Navig. 1997, 144, 219–226.

9. Cobra, D.T.; Oppenheim, A.V.; Jaffe, J.S. Geometric distortions in side-scan sonar images: A procedure for their estimation and correction. IEEE J. Ocean. Eng. 1992, 17, 252–268.

10. Clarke, J. Dynamic Motion Residuals in Swath Sonar Data: Ironing out the Creases. Int. Hydrogr. Rev. 2003, 4, 6–23.

11. Cervenka, P.; Moustier, C.D.; Lonsdale, P.F. Geometric corrections on sidescan sonar images based on bathymetry: Application with SeaMARC II and Sea Beam data. Mar. Geophys. Res. 1995, 17, 217–219.

12. Cervenka, P.; de Moustier, C. Postprocessing and corrections of bathymetry derived from sidescan sonar systems: Application with SeaMARC II. IEEE J. Ocean. Eng. 1994, 19, 619–629.

13. Coiras, E.; Petillot, Y.; Lane, D.M. Multiresolution 3-D reconstruction from side-scan sonar images. IEEE Trans. Image Process. 2007, 16, 382–390.

14. Fezzani, R.; Zerr, B.; Mansour, A.; Legris, M.; Vrignaud, C. Fusion of Swath Bathymetric Data: Application to AUV Rapid Environment Assessment. IEEE J. Ocean. Eng. 2019, 44, 111–120.

15. Stateczny, A.; Gronska, D.; Motyl, W. HydroDron—New step for professional hydrography for restricted waters. BGC Geomat. 2018, 226–230.

16. Crawford, A.; Connors, W. Performance Evaluation of a 3-D Sidescan Sonar for Mine Countermeasures. In Proceedings of the OCEANS 2018 MTS, Charleston, SC, USA, 22–25 October 2018.

17. Ai, Y.; Armstrong, S.; Fleury, D. Evaluation of the Klein Hydrochart 3500 Interferometric Bathymetry Sonar for Noaa Sea Floor Mapping. In Proceedings of the OCEANS 2015 MTS, Washington, DC, USA, 18–21 May 2015.

18. Brisson, L.; Wolfe, D.; Staley, M. Interferometric swath bathymetry for large scale shallow water hydrographic surveys. In Proceedings of the Canadian Hydrographic Conference, St. John’s, NL, Canada, 14–17 April 2014.

19. Schimel, A.C.; Beaudoin, J.; Parnum, I.; LeBas, T.P.; Schmidt, V.; Keith, G.; Ierodiaconou, D. Multibeam sonar backscatter data processing. Mar. Geophys. Res. 2018, 39, 121–137.

20. Lucieer, V.; Roche, M.; Degrendele, K.; Malik, M.; Dolan, M.; Lamarche, G. User expectations for multibeam echo sounders backscatter strength data-looking back into the future. Mar. Geophys. Res. 2018, 39, 23–40.

21. LeBas, T.P.; Huvenne, V. Acquisition and processing of backscatter data for habitat mapping–comparison of multibeam and sidescan systems. Appl. Acoust. 2009, 70, 1248–1257.

22. Mitchell, G.A.; Orange, D.L.; Gharib, J.J.; Kennedy, P. Improved detection and mapping of deepwater hydrocarbon seeps: Optimizing multibeam echosounder seafloor backscatter acquisition and processing techniques. Mar. Geophys. Res. 2018, 39, 323–347.

23. Fakiris, E.; Blondel, P.; Papatheodorou, G.; Christodoulou, D.; Dimas, X.; Georgiou, N.; Kordella, S.; Dimitriadis, C.; Rzhanov, Y.; Geraga, M.; et al. Multi-Frequency, Multi-Sonar Mapping of Shallow Habitats—Efficacy and Management Implications in the National Marine Park of Zakynthos, Greece. Remote Sens. 2019, 11, 461.

24. LeBas, T.P.; Mason, D.C. Automatic registration of TOBI side-scan sonar and multi-beam bathymetry images for improved data fusion. Mar. Geophys. Res. 1997, 19, 163–176.

25. Yang, F.; Wu, Z.; Du, Z.; Jin, X. Co-registering and fusion of digital information of multibeam sonar and side-scan sonar. Geomat. Inf. Sci. Wuhan Univ. 2006, 31, 740–743.

26. Zhao, J.; Wang, A.; Guo, J. Study on fusion method of the block image of MBS and SSS. Geomat. Inf. Sci. Wuhan Univ. 2013, 38, 287–290.

27. Yan, J. Acquisition and superposition of the high-quality measurement information of multibeam echo sonar. Acta Geod. Et Cartogr. Sin. 2019, 48, 400.

28. Fakiris, E.; Papatheodorou, G. Quantification of regions of interest in swath sonar backscatter images using grey-level and shape geometry descriptors: The TargAn software. Mar. Geophys. Res. 2012, 33, 169–183.

29. Wang, A.; Zhang, H.; Wang, X.; Shang, X. Processing Principles of Side-scan Sonar Data for Seamless Mosaic Image. J. Geomat. 2017, 42, 26–29.

30. Cervenka, P.; de Moustier, C. Sidescan sonar image processing techniques. IEEE J. Ocean. Eng. 1993, 18, 108–122.

31. Ye, X.; Yang, H.; Li, C.; Jia, Y.; Li, P. A Gray Scale Correction Method for Side-Scan Sonar Images Based on Retinex. Remote Sens. 2019, 11, 1281.

32. Capus, C.G.; Banks, A.C.; Coiras, E.; Ruiz, I.T.; Smith, C.J.; Petillot, Y.R. Data correction for visualisation and classification of sidescan SONAR imagery. IET Radar Sonar Navig. 2008, 2, 155–169.

33. Bay, H.; Ess, A.; Tuytelaars, T.; Gool, L. Speeded-Up Robust Features (SURF). Comput. Vis. Image Underst. 2008, 110, 346–359.

34. Sedaghat, A.; Mohammadi, N. High-resolution image registration based on improved SURF detector and localized GTM. Int. J. Remote Sens. 2019, 40, 2576–2601.

35. Tola, E.; Lepetit, V.; Fua, P. Daisy: An efficient dense descriptor applied to wide-baseline stereo. IEEE Trans. Pattern Anal. Mach. Intell. 2009, 32, 815–830.

36. Li, J.; Hu, Q.; Ai, M.; Zhong, R. Robust feature matching via support-line voting and affine-invariant ratios. ISPRS J. Photogramm. Remote Sens. 2017, 132, 61–76.

37. Ye, Y.; Shen, L.; Hao, M.; Wang, J.; Xu, Z. Robust optical-to-SAR image matching based on shape properties. IEEE Geosci. Remote Sens. Lett. 2017, 14, 564–568.

38. Shechtman, E.; Irani, M. Matching local self-similarities across images and videos. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Minneapolis, MN, USA, 17–22 June 2007; pp. 1–8.

39. Fischler, M.A.; Bolles, R.C. Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. Commun. ACM 1981, 24, 381–395.

40. Chailloux, C.; Caillec, J.; Gueriot, D.; Zerr, B. Intensity-Based Block Matching Algorithm for Mosaicing Sonar Images. IEEE J. Ocean. Eng. 2011, 36, 627–645.

41. Zhao, J.; Wang, A.; Zhang, H.; Wang, X. Mosaic method of side-scan sonar strip images using corresponding features. IET Image Process. 2013, 7, 616–623.

42. Wang, A.; Zhao, J.; Guo, J.; Wang, X. Elastic Mosaic Method in Block for Side Scan Sonar Image Based on Speeded-Up Robust Features. Geomat. Inf. Sci. Wuhan Univ. 2018, 43, 697–703.

43. Daniel, S.; Leannec, F.L.; Roux, C.; Solaiman, B.; Maillard, E.P. Side-Scan Sonar Image Matching. IEEE J. Ocean. Eng. 1998, 23, 245–259.

44. Yang, L.; Jiao, Y.; Xu, J. Underwater target positioning technology of side scan sonar based on ultra short baseline. China Harb. Eng. 2017, 37, 6–9.

45. Feldens, P.; Schulze, I.; Papenmeier, S.; Schönke, M.; Schneider von Deimling, J. Improved Interpretation of Marine Sedimentary Environments Using Multi-Frequency Multibeam Backscatter Data. Geosciences 2018, 8, 214.

46. Tamsett, D.; McIlvenny, J.; Watts, A. Colour Sonar: Multi-Frequency Sidescan Sonar Images of the Seabed in the Inner Sound of the Pentland Firth, Scotland. J. Mar. Sci. Eng. 2016, 4, 26.

47. Hughes Clarke, J.E. Multispectral Acoustic Backscatter from Multibeam, Improved Classification Potential. In Proceedings of the US Hydrographic Conference, National Harbor, MD, USA, 16–19 March 2015; pp. 1–18.

48. Huff, L.C. Acoustic Remote Sensing as a Tool for Habitat Mapping in Alaska Waters. In Marine Habitat Mapping Technology for Alaska; Reynolds, J.R., Greene, H.G., Eds.; Alaska Sea Grant, University of Alaska Fairbanks: Fairbanks, AK, USA, 2008; pp. 29–46.

49. Chen, J.L.; Summers, J.E. Deep neural networks for learning classification features and generative models from synthetic aperture sonar big data. Acoust. Soc. Am. J. 2016, 140, 3423.

50. Song, Y.; He, B.; Liu, P.; Yan, T. Side scan sonar image segmentation and synthesis based on extreme learning machine. Appl. Acoust. 2019, 146, 56–65.

51. Chen, Z.; Liu, L.; Sa, I.; Ge, Z.; Chli, M. Learning context flexible attention model for long-term visual place recognition. IEEE Robot. Autom. Lett. 2018, 3, 4015–4022.