Reconstructing Seabed Topography from Side-Scan Sonar Images with Self-Constraint

Abstract

To obtain the high-resolution seabed topography and overcome the limitations of existing topography reconstruction methods in requiring external bathymetric data and ignoring the effects of sediment variations and Side-Scan Sonar (SSS) image quality, this study proposes a method of reconstructing seabed topography from SSS images with a self-constraint condition. A reconstruction model is deduced by Lambert’s law and the seabed scattering model. A bottom tracking method is put forward to get the along-track SSS towfish heights and the initial seabed topography in the SSS measuring area is established by combining the along-track towfish heights, towfish depths and tidal levels obtained from Global Navigation Satellite System (GNSS). The complete process of reconstructing seabed topography is given by taking the initial topography as self-constraint and the high-resolution seabed topography is finally obtained. Experiments verified the proposed method by the data measured in Zhujiang River, China. The standard deviation of less than 15 cm is achieved and the resolution of the reconstructed topography is about 60 times higher than that of the Digital Elevation Model (DEM) established by bathymetric data. The effects of noise, suspended bodies, refraction of wave in water column, sediment variation, the determination of iteration termination condition as well as the performance of the proposed method under these effects are discussed. Finally, the conclusions are drawn out according to the experiments and discussions. The proposed method provides a simple and efficient way to obtain high-resolution seabed topography from SSS images and is a supplement but not substitution for the existing bathymetric methods.

1. Introduction

High-resolution seabed topography has important applications in offshore oil exploration [1], underwater target detection [2], marine eco-environmental protection [3,4], underwater navigation and localization [5] and many other fields. Acoustic waves move quite efficiently through water and the ability to travel over such great distances allows remote sensing in a water environment, which is helpful for ocean observation and monitoring. Devices that use acoustic waves in such applications fall under the family of instruments known as sonars. Currently, the most widely used sonars in the bathymetric field are single-beam bathymetric system (SBS) and multibeam bathymetric system (MBS). Using SBS to obtain high-resolution seabed terrain needs dense surveying lines, which are low-efficiency and high-cost. Using MBS, one can carry out a large-scale bathymetric measurement, but the interval of surveying points in a ping will be enlarged and the bathymetric resolution will be decreased with the increases of water depth and beam incident angle [6,7,8]. Side-scan sonar (SSS) is widely used to reflect seabed surface features. The SSS transducer built in a towfish is usually towed behind a surveying vessel by a cable, it emits a wide-angle beam and receives the seabed echoes at fixed time intervals to form the seabed image. The resolution of a SSS image is about 20–100 times larger relative to that of MBS bathymetric topography. Consequently, SSS images can vividly reflect the seabed targets, topography and sediment reflectivity [9,10]. If the micro-topography features can be determined from SSS images, it will be helpful to obtain the high-resolution seabed topography in deep sea. In addition, in the shallow waters, although the MBS can achieve resolution of decimeter scale, the reconstruction method from SSS images can also provide an efficient and low-cost supplement to obtain high-resolution seabed topography.

Shape from shading (SFS) is an efficient technique to obtain the three-dimensional (3D) shape from the intensity variation of the two-dimensional (2D) image. In SFS, the reconstruction of the 3D shape is based on the illumination theory and the imaging mechanism [11,12,13]. As early as 1999, Bell et al. studied the SSS imaging mechanism, analyzed the effects of acoustic wave incident directions on the seabed texture presentations in SSS images, and qualitatively discussed the relationship between the incident directions and the seabed texture images [14]. The research implies that the seabed topography reconstruction is related to the incident direction. On this basis, many scholars have developed studies on reconstructing 3D seabed topographies using 2D SSS images. Johnson et al. studied the 3D shape reconstruction of seabed topography by different seabed scattering models, and concluded that Lambert’s law is appropriate [15,16].

To get absolute seabed topography from the SSS image, many scholars developed different seabed reconstruction methods [17,18,19,20,21,22]. These studies can be classified into two categories according to whether the reconstruction needs the external bathymetric data. Under the support of external bathymetric data, Lange et al. obtained the seabed topography of every ping by a propagation method [17]. The method obtained the position and seabed height of the first echo of each ping using the external bathymetric data firstly, and then built seabed elevation maps from SSS waterfall images according to the propagation principle of acoustic wave in water. Johnson et al. used external bathymetric data to build an initial terrain for the reconstruction process from SSS images [15,16]. Bikonis K et al. combined MBS and SSS data to build the statistical curves, namely the seabed backscattering coefficient angular dependence curves to reconstruct the seabed relief in a flat experimental area. As stated by Bikonis et al., the method needed to be further verified and improved for wide applications [18]. Zhao and Wang et al. studied the relationship between the reconstructed seabed topography and the bathymetric data, and discovered the strong correlation between them. A large-scale seabed topography inversion was carried out by Lambert’s law, and the absolute seabed topography was obtained through adjusting the inversion result by the correlation model [19,20]. The above studies need the bathymetric data as an external constraint, and thus the reconstruction results may be easily affected by the data quality [19]. Moreover, an extra SBS or MBS measurement has to be performed. Without the external bathymetric data, Dura et al. analyzed the propagation method and pointed out that the method simplifies the reconstruction of seabed topography, but it is more sensitive to noise and abnormal ping measurement, and thus they proposed a linear method to reconstruct the textured seabed from synthetic SSS images [21]. Using the linear method, Dura et al. reconstructed the relative shapes of sand ripples [21]. The quantitative evaluation of the reconstruction showed that the linear method had better robustness than the propagation method [21]. The conclusion was achieved by simulation experiments and it needs to be further verified by the actual data. Besides, Coral et al. used an assumed flat initial seabed to reconstruct the absolute seabed topography from the SSS image by combining the expectation-maximization and multiresolution algorithms, but they obtained the reconstructed small-size shapes with large quantization errors [22].

To simplify the reconstruction and get accurate seabed topography without the constraint of external bathymetric data, this paper proposes a novel method to reconstruct seabed topography from SSS images with the SSS self-constraint which can be obtained by combining the SSS towfish height from the SSS waterfall image, the SSS towfish depth provided by a built-in pressure sensor and GNSS tidal level observed by the GNSS receiver mounted in surveying vessel. The paper is organized as follows. Section 2 deduces the reconstruction model using SSS images. Section 3 solves the model using Newton iteration algorithm, gives a method to extract the along-track depth from SSS waterfall images, and proposes the process of reconstructing 3D seabed topography using SSS images by taking the extracted along-track depth as self-constraint. Section 4 designs experiments to verify the proposed method and theoretically analyzes the results and assesses the accuracies of reconstructed seabed topographies. Section 5 discusses the influencing factors on the proposed method. Section 6 draws the conclusions according to the experiments and the discussions.

2. Reconstruction Model of 3D Seabed Topography

The acoustic waves emitted from the SSS transducer propagate through sea water to seabed, and the echoes from the seabed are received by the transducer at fixed time intervals. SSS images are formed by converting the backscatter strength (BS) to the gray level. The BS or gray level is related to the incident acoustic energies and directions, the seabed albedo and the topography gradients [9]. It is generally believed that the seabed is mainly comprised of diffuse scattering characteristics, and Lambert’s law is an appropriate model to describe the scattering process [23,24,25,26,27]. In optics, Lambert’s law says that the radiant intensity or luminous intensity observed from an ideal diffusely reflecting surface or ideal diffuse radiator is directly proportional to the cosine of the angle θ between the direction of the incident light and the normal vector [28]. This law is also suitable for acoustics measurement. A lot of coarse surfaces which obey Lambert’s law are called Lambertian and exhibit Lambertian reflectance. When given the parameter of incident direction, the Lambert’s law can be used to redistribute the acoustic energy [23]. Figure 1 describes the relationship between the BS and the seabed topography gradients [23].

E_s is proportional to the cosine of the incident angle θ, and can be expressed as

$$E_s=u{E}_i\cos \theta \cos \delta dS$$

(1)

where, u is the albedo, other parameters have the same meanings as shown in Figure 1. Let A = E_icosδdS, Equation (1) can be simplified as

$$E_s=Au\cos \theta$$

(2)

Equation (2) quantitatively describes the relationship between BS and the topography gradient. Figure 2 shows the situation of starboard transducer receiving echo under towfish reference frame (TFS). TFS is a right-hand coordinate system. The origin point, x, y and z axles of TFS are defined as the center of SSS transducer, the ping scanning line direction, the tracking direction and the direction perpendicular to the xoy surface.

Figure 2 simplifies the SSS transducer as a point source, and treats the propagation path as a straight line. After time varying gain (TVG) correction and radiation correction, the acoustic waves hitting on the seabed surface have the same incident energy [29,30,31]. For the homogeneous seabed, u is often set as a constant. Define I_norm = E_s/Au

$$I_{norm}(x,y)=\cos (\theta (x,y))$$

(3)

where, I_norm is called the normalized intensity, (x, y) are the coordinates of point g, θ reflects the topography gradient variations. The relationship of incident angle θ, incident vector $\stackrel{\rightharpoonup }{r}$ and the normal vector $\stackrel{\rightharpoonup }{N}$ can be expressed as

$$\cos (\theta (x,y))=-\frac{\overrightarrow{r_{xy}}·\overrightarrow{N_{xy}}}{|\overrightarrow{r_{xy}}|·|\overrightarrow{N_{xy}}|}$$

(4)

$$\overrightarrow{r_{xy}}=（\cos {\tau }_{xy}\sin {\phi }_{xy},\sin {\tau }_{xy}\sin {\phi }_{xy},\cos {\phi }_{xy}）$$

(5)

$$\begin{array}{c}{\overrightarrow{N}}_{xy}=(-p_{xy},-q_{xy},1)\\ p_{xy}=\frac{\partial z}{\partial x}=z_{x,y}-z_{x-1,y},q_{xy}=\frac{\partial z}{\partial y}=z_{x,y}-z_{x,y-1}\end{array}$$

(6)

where, τ_xy is the tilt angle between the plane component of incident direction and x-axis, p_xy and q_xy are topography gradients along x and y directions, z is topography height, z_x−1,y and z_x,y−1 are the neighbor heights of current pixel (x, y). In TFS, τ_xy is zero. ϕ_xy can be computed by towfish height and initial seabed terrain.

$$\begin{array}{c}{\phi }_{xy}=\arctan (\frac{x}{\mathrm{abs}(h_{xy})})\\ h_{xy}=h_x+(zi{n}_{xy}-zi{n}_x)\end{array}$$

(7)

where, x is the horizontal distance, h_xy is the vertical distance from the point g to the xoy surface, h_x is towfish height, zin_xy is the initial depth value of point g, zin_x is the initial depth value of the point beneath the SSS transducer. Other parameters have the same meanings as in Figure 2.

Combining Equations (3)–(7) to compute z, the reconstructed topography can be obtained.

3. Reconstructing Seabed Topography from SSS Image

3.1. Solution of Reconstruction Model

To take Equations (4)–(7) into Equation (3), Equation (8) is established

$$I_{norm}(x,y)=\frac{\cos {\phi }_{xy}-(z_{x,y}-z_{x-1,y})\sin {\phi }_{xy}}{\sqrt{1+{(z_{x,y}-z_{x-1,y})}^2+{(z_{x,y}-z_{x,y-1})}^2}}$$

(8)

Transform the Equation (8) into Equation (9) as

$$\begin{array}{ll}0& =f(I_{norm}(x,y),z_{x,y},z_{x-1,y},z_{x,y-1})\\ & =I_{norm}(x,y)-\frac{\cos {\phi }_{xy}-(z_{x,y}-z_{x-1,y})\sin {\phi }_{xy}}{\sqrt{1+{(z_{x,y}-z_{x-1,y})}^2+{(z_{x,y}-z_{x,y-1})}^2}}\end{array}$$

(9)

By Taylor series expansion, Equation (9) is depicted as

$$0=f^{n-1}(I_{norm}(x,y),z_{x,y}^{n-1},z_{x-1,y}^{n-1},z_{x,y-1}^{n-1})+\Gamma (1)+\Gamma (2)+\cdots +\Gamma (n)+\cdots$$

(10)

where, fⁿ⁻¹ is the (n − 1)th iteration result or the nth initial value, Г(1), Г(2), …, Г(n) are the 1st − nth order terms, zⁿ and zⁿ⁻¹ are the seabed heights achieved by nth and (n − 1)th iterations, z_x,y, z_x−1,y and z_x,y−1 are the seabed heights of pixel (x, y), (x − 1, y) and (x, y − 1). If assuming that seabed topography changes slowly, the linear part, Г(1) in Equation (10) plays the most important role in solving the reconstruction model [32]. By Taylor series expansion of up to first order terms, the linear approximation of the Equation (10) is obtained as

$$\begin{array}{ll}0& =f(I_{norm}(x,y),z_{x,y},z_{x-1,y},z_{x,y-1})\\ & \approx f^{n-1}(I_{norm}(x,y),z_{x,y}^{n-1},z_{x-1,y}^{n-1},z_{x,y-1}^{n-1})+\Gamma (1)\\ & =f^{n-1}(I_{norm}(x,y),z_{x,y}^{n-1},z_{x-1,y}^{n-1},z_{x,y-1}^{n-1})+a_{x,y}(z_{x,y}-z_{x,y}^{n-1})+a_{x-1,y}(z_{x-1,y}^{}-z_{x-1,y}^{n-1})+a_{x,y-1}(z_{x,y-1}-z_{x,y-1}^{n-1})\end{array}$$

(11)

where, $a_{x,y}=\frac{\partial f(I_{norm}(x,y),z_{x,y},z_{x-1,y},z_{x,y-1})}{\partial z_{x,y}^{n-1}}$, a_x−1,y and a_x,-y−1 have the similar formation as a_x,y.

The Equation (11) can also be transformed as

$$\begin{array}{ll}0& =f(I_{norm}(x,y),z_{x,y},z_{x-1,y},z_{x,y-1})\\ & =a_{x,y}{z}_{x,y}+a_{x-1,y}{z}_{x-1,y}^{}+a_{x,y-1}{z}_{x,y-1}-a_{x,y}{z}_{x,y}^{n-1}-a_{x-1,y}{z}_{x-1,y}^{n-1}-a_{x,y-1}{z}_{x,y-1}^{n-1}\\ & +f^{n-1}(I_{norm}(x,y),z_{x,y}^{n-1},z_{x-1,y}^{n-1},z_{x,y-1}^{n-1})\end{array}$$

(12)

To solve the Equation (12) by Newton iteration algorithm [33], Equation (13) is given as

$$\begin{array}{cc}{z}_{x,y}^n& =z_{x,y}^{n-1}-f(z_{x,y}^{n-1})/\frac{\partial f(z_{x,y}^{n-1})}{\partial z_{x,y}^{}}\\ & =z_{x,y}^{n-1}+\frac{1}{k^{n-1}}(-f(z_{x,y}^{n-1}))\end{array}$$

(13)

$$k^{n-1}=\frac{\sin {\phi }_{xy}}{\sqrt{1+{(z_{x,y}^{n-1}-z_{x-1,y}^{n-1})}^2+{(z_{x,y}^{n-1}-z_{x,y-1}^{n-1})}^2}}-\frac{(z_{x,y}^{n-1}-z_{x-1,y}^{n-1}+z_{x,y}^{n-1}-z_{x,y-1}^{n-1})(\cos {\phi }_{xy}-(z_{x,y}^{n-1}-z_{x-1,y}^{n-1})·\sin {\phi }_{xy})}{\sqrt{{(1+{(z_{x,y}^{n-1}-z_{x-1,y}^{n-1})}^2+{(z_{x,y}^{n-1}-z_{x,y-1}^{n-1})}^2)}^3}}$$

By iteration operation, the final topographic height zⁿ_xy is obtained and the seabed topography is reconstructed. In the reconstruction, the initial topography and the iteration terminal condition are required. Given the initial seabed topography z⁰_x,y, the final seabed topography can be reconstructed after iteration operation. The iteration can be stopped by a given threshold, ε. ε can be determined from the convergence curve of the iteration, which is shown in Section 5.5.

$$z_{x,y}^n-z_{x,y}^{n-1}\le \epsilon$$

(14)

3.2. Bottom Tracking and Initial Seabed Topography

Equation (13) shows that the seabed topography can be reconstructed from SSS images if the initial topography z⁰ is given. z⁰ is often provided by external bathymetric data, which means an extra measurement must be done. To simplify the work, the following proposes a method to obtain z⁰ by bottom tracking and combining towfish height, depth and tidal level.

In the SSS measurement, the first received echo is from the SSS transducer nadir. During the process of transmitting the beam and receiving the first echo, the transducer records nothing, and thus a blank area is formed in the SSS waterfall image. The width of the blank area represents the SSS towfish height [34,35]. If the first echo can be detected, the towfish height can also be obtained. Because of receiving nothing, BSs of the blank area are far lower than those of the seabed area. Therefore, the boundary point of two areas in a ping scanning line, namely the seabottom point, can be detected by calculating the BS differences of adjacent pixels and finding the maximum difference.

$$i_0=\max \{\Delta G_i=G_{i+1}-G_i,i=1,2,3,\dots \}$$

(15)

or

$$i_0=\max \{\Delta B{S}_i=B{S}_{i+1}-B{S}_i,i=1,2,3,\dots \}$$

where, i₀ is the pixel number or position of detected seabottom point in a ping sequences. G (BS) denotes gray level (backscatter strength), ΔG_i (ΔBS_i) is the gray level (BS) difference of i + 1 and i pixels. Then, the towfish height h_towfish can be calculated by

$$h_{\mathrm{towfish}}=\frac{S}{W}{i}_0$$

(16)

where, W and S are the width and the maximum slant range of the waterfall image.

To combine the towfish depth from the built-in pressure sensor and the tidal level measured by GNSS [36,37], the seabed height can be obtained by

$$z=T-(h_{\mathrm{towfish}}+D_{\mathrm{towfish}})$$

(17)

where, z is the seabed height, D_towfish is the towfish depth and T is the tidal level.

For the SSS image of a surveying line, the along-track seabed topographic sequence can be obtained by combining the towfish height derived from the waterfall image by the bottom tracking, towfish depth from the built-in pressure sensor and the tidal level measured by GNSS mounted on surveying ship. To reconstruct the seabed topography of the surveying line, the initial topography can be created by the method of triangulation with linear interpolation for the along-track seabed topographic sequences of the surveying line and those of its two adjacent surveying lines. Because the along-track topographic data is sparse in the coverage of the three SSS surveying lines, the created initial topography is coarse and used as the initial terrain in Equation (13) for reconstructing the seabed topography of the surveying line from the SSS image. Adopting a similar process to the above, the initial topographies can also be obtained for seabed reconstructions of other surveying lines. The process is depicted in Figure 3, which shows the initial topography is obtained from SSS waterfall images instead of the extra bathymetric data. Therefore, the initial topography or the constraint condition in the reconstruction is called self-constraint, which simplifies the reconstruction while saving costs.

3.3. Assessment

To assess if the system works, the evaluation of the reconstruction result needs to be implemented through computing the biases between the reconstructed data and the real bathymetric one. Using Equation (18), the maximum, minimum, mean and standard deviation of the biases can be computed and the probability distribution function (PDF) curves can be drawn. The obtained parameters and curves are used to assess the reconstruction method.

$$STD=\sqrt{\frac{{\displaystyle \sum _{i=1}^n\Delta D\Delta D}}{n}}，\Delta D=z-z_{\mathrm{real}}$$

(18)

where, ΔD is the bias between the reconstructed data and the real bathymetric one, STD is the standard deviation, n is the number of points used in the assessment process.

3.4. Process of Reconstructing 3D Seabed Topography

The reconstruction process can be described in detail below (Figure 4):

(1)

Detect the seabottom points from SSS waterfall images to get the towfish heights, and then obtain the initial seabed topography by that depicted in Figure 3.

(2)

Preprocessing the SSS image by sound ray tracing, slant distance correction, radiometric distortion correction, gray equalization, etc. [29,30,31].

(3)

Compute the topography gradients by Equation (6), the angle ϕ by Equation (7), and the seabed topography by Equation (13) using the self-constraint.

(4)

Execute iteration until the difference of two adjacent iterations is less than the given threshold ε. In the iteration, the next calculation uses the last result as the initial condition, and step (2) and step (3) are repeated until the difference meets Equation (14).

(5)

Assess the reconstruction by referring to the actual bathymetry data.

4. Experiments and Analysis

4.1. Experimental Area and Data Preparation

To verify the proposed method, an experiment was carried out in Zhujiang River with 10–14 m water depth, obvious seabed features and homogeneous sediment. In the experiment, Edgetech 4125 with 400 kHz sampling frequency, 2 cm lateral resolution, 0.6 m longitudinal resolution and 0.3° horizontal beam width, was adopted for the SSS measurement. The built-in pressure sensor with the accuracy of 0.5% depth was used to record the SSS towfish depth. Three surveying lines with the lengths of 1020 m, 900 m and 540 m were designed. The adjacent line interval and the SSS swath width are 200 m, and the overlapping ratio between adjacent lines is 50%. To assess the topography reconstructed by the proposed method, 20 surveying lines with 40 m interval were also accomplished by HY1600 with 208 kHz operating frequency, 8° beam angle, ±(0.01 m + 0.1% depth) bathymetric accuracy and 0.5 m along-track sampling interval in the area. In the measurement, the tidal level was also synchronously recorded. Raw SSS data were recorded in *.xtf files. To decode these files and deal with those by radiometric correction, slant range correction, geocoding and image mosaics, the image with the resolution of 0.6 m is formed as shown in Figure 5a. The features in the image are marked as ①–⑫. The bathymetric data are also processed by quality control, draft correction and tidal correction, and the Digital Elevation Model (DEM) is built as shown in Figure 5b. Features ①–⑫ are displayed clearly in Figure 5b. The results indicate the SSS image can reflect seabed features, and they also imply that the 3D shapes of these features can be recovered from a 2D SSS image. The corresponding 2D colored images are shown in Appendix A, which can also be used to compare.

4.2. Reconstructing Seabed Topography

To depict the reconstruction process in detail, Strip I is chosen. Before reconstructing the seabed topography, Strip I needs to be preprocessed. Firstly, the raw XTF file is decoded to form the waterfall image as shown in Figure 6a. Being affected by the propagation loss and absorption loss, the BSs will decrease with the increase of propagation range. TVG can compensate these losses to some extent. Being affected by beam pattern and TVG residual (or imperfect TVG correction), the radiometric distortion degrades the quality of the SSS images and affects the reconstruction process. A radiometric correction method for side-scan sonar images in consideration of seabed sediment variation is adopted before using the SSS images to reconstruct seabed topography [31]. After the correction, the effects of TVG and beam pattern on the topography reconstruction can be weakened. Relative to these in Figure 6a, the gray levels change more evenly, and the seabed features display more clearly in Figure 6b after the TVG and radiation correction. In the corrected image, the effect of acoustic incident direction is removed. After the bottom tracking from Figure 6a, the along-track towfish heights are achieved, and the slant-range correction is carried out and shown in Figure 6c. It can be found that the water column area disappears from the waterfall image, and the actual seabed image scanned by SSS transducer is obtained. Then, the image is geocoded and shown in Figure 6d. Due to homogeneous sediment in the measurement area, the effect of seabed sediment can be ignored, and the images only reflect the seabed topography features.

In the above preprocessing, the bottom tracking is very important to obtain towfish height and the initial seabed topography. In Figure 6a, it can be found that the blank area or water column area differs distinctly from the seabed areas in gray levels. According to the characteristic, the bottom tracking can be carried out. For a ping of Strip I as shown in Figure 7a, the gray level sequence is firstly extracted as shown in Figure 7b. Then, to calculate the gray level differences of any two adjacent echoes and find the maximum difference by Equation (15), the seabottom point can be detected as shown in Figure 7c, and the towfish height can be computed by Equation (16). To detect the seabottom points in other pings, the sequence of along-track towfish heights is obtained as shown in Figure 8. To avoid the effect of the water-column noise, the bottom tracking is smoothed by moving averaging for every 20 pings. Comparing sequences of the towfish heights and the corresponding seabed heights from bathymetric data, as shown in Figure 8, a strong similarity can be found between them except the various biases caused by the variations of tidal levels and towfish depths. Using the towfish heights, the slant-range correction can be performed and this is shown in Figure 6c. By geocoding, the image is obtained as shown in Figure 6d and it will be used to reconstruct the seabed topography. By combining the towfish heights, towfish depths and tidal levels, the along-track seabed height sequence is obtained and shown in Figure 8. It can be seen that the seabed height has high consistency with that from the SBS bathymetric data. Comparing the two sets of seabed heights, the statistic parameters of the biases are listed in

Table 1. Statistical parameters of seabed height biases.

Max. (m)	Min. (m)	Mean (m)	STD (±m)
0.28	−0.36	0.00	0.13

, which shows that the along-track depth data achieved by the proposed method has the standard deviation of 0.13 m and can be used to build the initial seabed topography.

The similar process is used for Strip II and III, the sequences of along-track seabed heights are obtained and the seabed DEM can be established by a linear interpolation. The DEM as shown in Figure 9 will be used as the initial seabed topography in the reconstruction of Strip I. Using the initial seabed topography shown in Figure 9 as self-constraint, the seabed topography is reconstructed from the image as shown in Figure 6d by the process depicted in Figure 4. The reconstruction result is shown in Figure 10a. The reconstructed terrain grid size is 0.6 m × 0.6 m, which is the same as the image pixel resolution representing the actual. The resolution of terrain grid data is improved significantly relative to the bathymetric terrain with terrain grid size of 40 m × 40 m. It can be seen that the topography vividly presents the six features as shown in Figure 6d.

For convenient comparison, an interpolation is performed for the SBS bathymetric data to get the same data resolution as that of the reconstructed topography. Two DEMs in Figure 10 are established by using the reconstruction data and the interpolated bathymetric data. Taking the bathymetric DEM as reference, the reconstructed DEM in Figure 10a has high similarity with the bathymetric one in Figure 10b on the whole trend, but reflects more refined seabed features. To assess the capability of the reconstruction in reflecting microtopography, frequency analysis is adopted for the two DEMs. The amplitude-frequency curves are shown in Figure 11a. It can be seen that the reconstruction has higher amplitude than the bathymetric data at every frequency, especially at mid-high frequency band. The amplitude-frequency curves show that the reconstruction has better performance in reflecting the microtopography. According to variations of the two curves, the cut-off frequency as 0.02 Hz (Figure 11a) is chosen for dividing the DEMs into high-frequency and low-frequency parts. The high-frequency part and low-frequency part differences of two DEMs are shown in Figure 11b,c. It can be seen that the low-frequency part differences are less than 0.1 m and the high-frequency differences vary mainly between 0.0–0.2 m and then between 0.2–0.4 m at some specific positions, which shows that the reconstructed DEM has consistent terrain tendency with the bathymetric one and can reflect the seabed features with the height more than 0.2 m at least.

To assess the accuracy of reconstruction, the external bathymetric data is used as a reference, and the biases of the reconstructed seabed topography at these bathymetric points are calculated. The statistical results of the biases are listed in

Table 2. Statistical parameters of the reconstruction biases of Strip I relative to the bathymetric data.

Max.	Min.	Mean	Standard Deviation
0.30	−0.43	0.00	0.12

. It can be seen that the maximum, minimum, mean and standard deviation of the biases are 0.30 m, −0.43 m, 0.00 m and 0.12 m and the proportion of absolute biases less than 0.30 m is 99.22%. The reconstruction achieves nearly the same accuracy as the bathymetric data. According to the reconstruction accuracy and the ability to reflect small-size terrain features, the reconstruction can reflect the seabed feature with the height more than 0.2 m at least and achieve the standard deviation better than 0.15 m.

To analyze how the reconstruction accuracies vary with the across-track distances, the bathymetric data of the eight SBS surveying lines symmetrically distributed on both sides of the SSS towfish tracking line in Figure 12a are used as the reference to assess the reconstruction accuracies at the 8 lines. The reconstruction errors at these 8 lines are calculated and the statistical parameters and PDF curves of the biases are shown in

Table 3. Statistical parameters of the reconstruction biases at different bathymetric lines.

	Line	1	2	3	4	5	6	7	8
Bias
Max./m		0.39	0.39	0.12	0.39	0.26	0.38	0.22	0.34
Min./m		−0.10	−0.25	−0.09	−0.20	−0.30	−0.29	−0.35	−0.10
Mean/m		0.12	0.04	0.01	0.07	0.00	0.01	0.03	0.10
Std/m		0.10	0.10	0.04	0.10	0.09	0.11	0.10	0.08

and Figure 12b.

It can be seen from Figure 12 and Table 3 that the variation ranges and standard deviations of reconstruction biases at different across-track distance lines are closed to those shown in Table 2, but the systematic biases appear in the reconstruction results of line 1 and 8 because the corresponding mean biases are respectively 0.12 m and 0.10 m. The problem appearing at the outmost across-track locations may be induced by imperfect radiometric distortion correction and is discussed in Section 5.2.

4.3. Seabed Topography Reconstruction of Surveying Area

Using a similar method and process to that depicted in Section 4.2, the seabed topographies of Strip II and III are reconstructed. Combining the reconstruction results of three strips, the seabed DEM of the surveying area is established and shown in Figure 13. Compared with the DEM as shown in Figure 5b, the reconstructed topography shows good consistence with the real one. However, because of the resolution difference of the SSS image (0.6 m) and the SBS bathymetric data (40 m), the reconstruction improves the topography resolution about 66 times and reflects clearer seabed features relative to the real DEM. Besides the seabed features ①–⑫ in Figure 5, some micro features are also displayed in Figure 13.

To assess the accuracy in the surveying area, the biases between the whole reconstruction topography and the bathymetric data are calculated, and the statistical results and PDF curves of the reconstruction biases are shown in

Table 4. Statistical parameters of the biases of the reconstructed topographies in different areas.

	Biases (m)	Max.	Min.	Mean	Standard Deviation
Area
Entire measurement area		0. 47	−0.54	0.00	0.12
Overlapping area of Strip I & II		0.44	−0.54	−0.08	0.14
Overlapping area of Strip II & III		0.49	−0.38	0.03	0.10

and Figure 14a. Compared with Table 2, the reconstruction accuracy in the area is close to that of Strip I. The results show the proposed method has good applicability and consistent accuracy in the reconstruction. Using Strip II and III, the proposed method can also be evaluated in terms of the reconstruction consistency of the overlapping area of adjacent strips. Comparing topography differences of the overlapping area, the statistical results of these biases are shown in Table 4 and the corresponding PDF curves are shown in Figure 14b,c. It can be seen that the biases obey the normal standard distribution and the proportion of the absolute biases less than 0.3 m is 98.66%. Statistical results prove that the reconstruction of the overlapping area between adjacent strips is consistent.

The same phenomenon that some big biases exist in the reconstruction of Strip I also happens in other strips. The reason for this could be attributed to the abnormal bottom tracking and the image noise caused by the water environment. Nonetheless, the statistical results in Table 4 show that the proposed method achieves the standard deviation less than 0.15 m, and has consistent accuracy in the reconstructions of all SSS strips.

5. Discussion

5.1. Noise in SSS Image

The proposed method reconstructs the seabed topography from SSS images with a self-constraint. Once the SSS transducer receives abnormal echoes from fish, suspended solids or particles in the water column, abnormal features or noise will be formed in the SSS image [9,35] and this will result in abnormal topography reconstruction. The problem can be avoided by filtering the SSS image and the reconstructed topography. By obeying discrete distribution, the noise in the SSS image can be easily removed by an average filtering or the threshold filtering. In the threshold filtering, the pixel whose gray level difference is more than 10 away from its surrounding pixels will be filtered. Although most noise can be removed from the SSS image, some residual noise still remains and is reflected in the reconstructed terrain. The noise from the suspended solids has strong backscatter strengths. Thus, the reconstruction of the noise will differ obviously from those of its surrounding terrain and can be removed by the topographic trend filtering [7].

5.2. Eeffects of Refraction of Waves, Towfish Depth and Cross-Track Distance

5.2.1. Refraction of Waves in the Water Column

In general, in SSS image processing, the refraction of the wave is ignored and the sound ray is regarded as a straight line for the calculation of echo position in a ping measurement, which results in inaccurate target position in a SSS image. To weaken the effect, the sound ray tracing is added in the reconstruction process (2) of Section 3.4.

In the sound ray tracing, two parameters, the propagation time and incident angle of an echo are required. The propagation time can be obtained by the equal-interval time sampling in a ping measurement. The incident angle can be approximately estimated by the towfish height and the across-track distance of the echo. After the sound ray tracing, the positions of the echoes or seabed features in SSS image can be determined accurately.

5.2.2. Towfish Depth and Across-Track Distance

The footprint of a beam varies with towfish depth and across-track distance. The larger the towfish depth and the across-track distance, the smaller the backscatter intensities are. This phenomenon leads to the disequilibrium of gray level in the SSS image and the inaccurate presentation of sediment reflectivity as well as the incorrect reconstruction of seabed topography. The problem can be solved by the radiometric distortion correction and the slant distance correction [29,30,31]. The radiometric distortion correction can remove the effects of towfish depth change. The integration of the slant distance correction and radiometric distortion correction can weaken the effect of various across-track distances. After these corrections, the quality of the SSS image is improved (Figure 6) relative to the raw waterfall image. Using the improved SSS image to reconstruct seabed topography, nearly consistent accuracies were obtained in different across-track distances (Figure 12, Table 3). We can also find slight accuracy decreases with the increase of the across-track distances, which shows that the effects of towfish depth and across-track distance were not removed completely by the radiometric distortion correction used in this paper and needs to be weakened by further improving the model.

5.3. Seabed Sediments

According to Section 2, the gray level of a pixel in the SSS image or the backscatter strength of a beam footprint is related to the seabed sediment and topographic variation. Because the experimental area is of homogeneous sediment, the SSS images only reflect the variation of seabed topography, and thus the impacts of seabed sediments on the reconstruction are ignored. The above experiments have verified that the reconstruction has good results under the situation of homogeneous seabed sediment. For further applications of the proposed method in inhomogeneous areas, the sediment variations need to be considered. The inhomogeneous sediments have different acoustic impedances and reflectivity albedos [38], which affects the backscatter intensities (Figure 15a) and the reconstruction by Equation (2). If the effect is ignored, the reconstruction will result in a systematic bias in the reconstructed topographies of two different sediment areas due to using the same reflectivity albedo (Figure 15(c1,d1)). To remove the impacts of sediment variations on the reconstruction, an unsupervised classification on the SSS images [31,39] or the sediment sampling [38,40] in the reconstructed area needs to be implemented before the reconstruction. To some extent, the unsupervised classification may provide a convenient way to obtain approximate distributions of different sediments in the area [31]. According to the classifications, the SSS image is divided into several parts (Figure 15b). In the separated homogeneous sediment areas, the reconstructions of seabed topographies are carried out by the proposed method. As shown in Equation (3), the seabed shape is obtained by the normalized intensity I_norm in this paper, and a given constant reflectivity albedo u will not bring effect on the reconstruction of the homogeneous sediment seabed areas. Therefore, the entire seabed topography (Figure 15(c2)) is finally obtained by combining the reconstructed seabed topographies of different sediment areas. Compared with the bathymetric topography (Figure 15c), the reconstruction by considering the sediment variations (Figure 15(c2)) achieved higher accuracy (mean bias of 0.01 m, standard bias of 0.06 m, Figure 15(d2)) than that by ignoring the variations (mean bias of 0.17 m, standard bias of 0.21 m, Figure 15(d1)).

5.4. Bottom Tracking and Initial Seabed Topography

The initial topography is used as the constraint of reconstructing absolute topography in the proposed method. Without the constraint, the proposed method can only get the relative topographic variations [19]. Relative to the existing methods [15,16,17,18,19,20,21,22], the proposed method obtains the initial topography by bottom tracking from SSS waterfall images, but not the external bathymetric data. Therefore, the constraint is named self-constraint in this paper. Being affected by the suspended solids in the water column, the bottom tracking may be affected and this results in inaccurate towfish height tracked by the mutation of gray levels depicted in Section 3.2. Under the situations as shown in Figure 16, the symmetric principle, namely the detected port towfish height equals the starboard one and the continuous seabed variation principle are used for diagnosing and repairing the abnormal [35]. If one side of the water column is covered by suspended solids, the detected towfish height of the covered side will be replaced by that of the uncovered one (Figure 16b). If the two sides are both covered, the abnormal towfish height can be repaired by interpolating correct ones of its neighbor pings (Figure 16c).

5.5. Determination of Iteration Termination Threshold

In the reconstruction, Equation (13) is solved by an iteration operation. As shown in Equation (14), the iteration operation is performed using an initial z⁰ and stopped by a given threshold ε. It is well known that the Netwon method converges to a stable result when an accurate starting value is known [32,33]. To analyze the importance of initial terrain in the iteration, an inappropriate initial terrain and an appropriate one are respectively used for the iteration and the corresponding convergence curves are shown in Figure 17a,b. It can be seen that the curve converges difficultly to a stable value within finite iterations when an inappropriate initial terrain is provided. Reversely, the curve easily reaches a stable result after the 42nd iteration when an appropriate initial terrain is given. This experiment demonstrates that it is better to provide an appropriate initial terrain in the iteration. The method of how to get the initial terrain is depicted in Section 3.2. Figure 17b displays the iteration from unstable to stable changes, which implies that the difference of adjacent iterations corresponding to the inflection point of the convergence curve can be set as ε. In the experiments depicted in Section 4, ε is set as 0.05 m. Similarly, the method for determining ε can also be used for the other seabed topography reconstructions from SSS images.

6. Conclusions

This study proposes a new method to reconstruct detailed seabed topography from side-scan sonar image with self-constraint. Without the need of external bathymetric data, the self-constraint or the initial topography, which is approximate to the real one, is formed by combining the towfish height from bottom tracking, towfish depth from built-in pressure sensor and tide level obtained from the SSS measurement. Therefore, the robustness of the reconstruction is improved and the reconstruction work is significantly simplified. Besides, the effects of sediment variations can be weakened by the unsupervised classification. The quality of the SSS image that is influenced by noise, suspended bodies, refraction of wave in water and the cross-track distance etc. can be improved by SSS image preprocessing. As a result, the accuracy of the reconstruction is improved. These conclusions were proved by the experiments. In the experiments, the reconstructed seabed topography reflected refined micro-terrain features and achieved the standard deviation less than 0.15 m. Moreover, the resolution of the reconstructed topography is about 66 times higher than that of the real one.

The propose method, which provides an approach to get detailed seabed topography from the high-resolution SSS image, is a supplement but not substitution for the existing bathymetric methods, and the proposed method has wide applications in obtaining fine seabed topography, especially in deep waters.