An Improved Forward-Looking Sonar 3D Visualization Scheme of Underwater Objects

This paper presents an improved scheme for forward-looking sonar (FLS) visualization, including an echo acquisition method, a global threshold based on local outlier factor and a reconstruction method. Based on the model of line target in the shallow water, the paper presents an echo acquisition method with impulse response convolution and ray tracing, and considering the reverberation. Then, a transform matrix is presented for the coordinate transformation to obtain the point cloud. After that, the paper presents the global threshold based on local outlier factor (LOF) method to solving the problem of noise sensitivity in the global thresholding. Finally, the paper presents a novel reconstruction method, including the coarse step via random sample consensus (RANSC) and the refined step via Power crust. Both simulation and experimental results show that the scheme in this paper produced superior result compared with the state-of the art.


I. INTRODUCTION
In the complicated underwater environment, the detection of the underwater target and the acoustic imaging of targets are crucial to the marine exploration, and have received increasing attention in recent years [1]- [9]. Unlike the twodimensional imaging, the real-time three-dimensional (3D) imaging is expensive and difficult for implementation, which limits the practical utilization of 3D imaging [5], [10], [11].
To promote the practical application of 3D imaging sonar, the simulation of object echo is a way to verify the algorithm with low cost compared with sea trial. Palmese and Trucco [9] simulates the real scene of buried target, taking the target itself, seabed and sediment volume into account. However, it did not consider the interference of reverberation. In addition, Tang [14] simulates the echo waveform by highlight model, and the reverberation is not be considered either.
Highlight model indicate the echo of a complex target can be equivalent to assemble of several scatters in high frequency situation [12], [14]- [16]. Due to the roughness and the uneven The associate editor coordinating the review of this manuscript and approving it for publication was Zhan-Li Sun . of the seafloor in the real environment, reverberation can not be neglected. Generally, the seabed reverberation model includes ray model and normal model, Ellis and Dale [18], and ray model is more efficient and intuitive. The reverberation signal is obtained by calculation of envelope method in [17]. This method can simplify the simulation process, but can not describe actual situation.
To make better visual result, the segmentation and denoising are important in the visualization scheme. In recent years, many researchers have proposed various denoising and segmentation methods. Actually, in terms of noise reduction and quality improvement, FIR filter and thresholding have proved to allow a good compromise between computational complexity and performances compared with statistical modeling [4]. The fixed thresholding method has been widely applied in underwater 3D imaging [19]. Although this method can effectively eliminate the clusters, it also remove echoes of real targets. Compared to the fixed threshold method, the ostu method [20] which is a kind of global thresholding method, can obtain more appropriate threshold, but it is sensitive to the noise and there are still outliers after the threshold discrimination.  In this paper, the impulse response convolution method and the ray tracing method based on ray model are first employed to simulate the echo of the line target under the shallow water with the reverberation. And the influence and inhibition of reverberation are described. Then, the density method, which is called local density outlier in the data mining field, is proposed to solve the problem of the outlier after global threshold discrimination. In the end, combining the Power crust proposed by Amenta et al. [21] with the random sample consensus (RANSC), we can reconstruct the surface of the target reliably. Both simulation and sea trial verify the effectiveness of the scheme in the paper.
The remainder of this paper is organized as follows: The environment modeling is demonstrated In section II, which including target model, receiving array distribution and reverberation model, the formula derivation of signal form of object and reverberation are also included in this section. Section III describes the data processing procedures in detail, including global thresholding segmentation based on LOF method and power crust based on RANSC. Section IV focuses on the performance analysis and the experiment results. Finally, the conclusions are drawn in Section V. Figure 1 depicts the model of the line target with smooth surface, which is fixed to the wedge-shaped seabed. The diameter of the line target is less than the resolution of the sonar. And the receiving array is a plane array on the plane XOY. Suppose that the object surface can be scattered into N p points. The position of the i-th scatter is N i = (x Ni , y Ni , z Ni ), and the corresponding distance to the coordinate origin is n i = |N i |. The receiving array is depicted in Figure 2. Then the Fourier transform of the pressure of the N p point scatters at r = (x r , y r , 0) in receiving array can be expressed as:

II. ENVIRONMENT MODELING
where B i is the Rayleigh's scattering process for an unfixed rigid particle whose dimensions is very much less than the wavelength, Q(f ) is the Fourier transform of the emitted pulse q(t); a i is the radius of the point scatter; ρ is the densities of the propagating medium of seawater; ρ i is the densities of the propagating medium of scatter; k is the ompressibility of the propagating medium; k i is the ompressibility of the scatter. Assume that the cosine of the angle between r − N i and N i is close to 1. Figure 3 shows the model of ray acoustics reverberation. It shows that the sonic waves emitted by point source reach the scatters through different propagation paths, and then return to the receivers via different paths. The Reverberation at time t is the sum of the all the scattering intensities contributed at that time [22].
To obtain the signal of target and reverberation, the paper considers the multi-path of the sonic ray. The model of ray reverberation is derived based on the amplitude, propagation time and grazing angle calculated by BELLHOP. BELLHOP is a beam tracing model for predicting acoustic pressure fields in ocean environments [26], which uses the point scattering model to discretize the seafloor into I uniformly distributed scattering points corresponding to different distances. The impulse response function of each scatter can be expressed as: Scattering strength is important in reverberation simulation. They reflect boundary roughness effects and determine scattering energy and pattern. The paper combined (3) with the physical scattering model shown by Figure 4. Then the impulse response function with scattering characteristics is written as: where N is the number of scatter ring micro-units at one point on the object. The amplitude of the units obeys the normal distribution Amp ∼ N (0, 1), and the phase ϕ n obeys the average distribution. Then, the reverberation Rev i (t) can be derived as: Superimposing the reverberation of I points to acquire the reverberation Rev(t), and the echo can be obtained as combination of reverberation and target signal: Similarly, the target signal can be derived by using the point scattering model as reverberation generation. And the echo of the target is calculated according to (5) and regarding the target as nonuniform discrete grids, where τ mn = (L mn − L 11 ) /c is the time delay, L mn is the distance from the basic scatter N b to the (m, n) array element, (x mn , y mn , z mn ) is the position of the (m, n) array element. A m and A n are the amplitude of incident sonic ray and scattering sonic ray, respectively. θ m and θ n are the grazing angle of incident sonic ray and scattering sonic ray.

III. DATA PROCESSING A. 3D SPATIAL REPRESENTATION
The source emits the acoustic pulses and produces 128×128 beams are shown in Figure 5. And the diagram and geometry of plane array is shown in Figure 6. The 3D point cloud is acquired according to the sensing characteristics and geometry of plane array. Suppose that R k is the range of the k th slice in the scan scope, q is the direction vector of signal echo as expressed in (9), β is the elevation angle in the vertical plane and α is the azimuth angle in the horizontal plane. Then, the point cloud is expressed as below, where k p α,β is a 3D point at the azimuth angle α, the elevation β and the range k th . T (α, β) is the transformation matrix, which is defined as below:

B. DENOISING AND SEGMENTATION
Sonar point cloud is restricted by the background noise and the reverberation. They induce the target blur and the poor object identification. And it is necessary to extract the target from the background via denoising and segmentation before the 3D reconstruction. Generally, to eliminate the noise, a fixed threshold A is set and the elements which is greater than the threshold A are determined to be target. However, the fixed threshold is difficult to determine the optimal threshold. The global thresholding set an initial estimated threshold value and update the value by some rule, it can obtain a more accurate threshold by setting convergence error and the number of iterations.
To solve the problems of the noise sensitive and outliers, this section presents an adaptive global thresholding method based on local outlier factor (LOF). Outliers are defined by Hawkins as points that differ significantly from other data in the datasets. Currently, the outlier detection is mainly used in the data mining (DM) field and is mainly divided into four categories, i.e., outlier detection based on depth, outlier detection based on cluster outlier detection, outlier detection based on density and outlier detection based on range [23]. The LOF is a kind of outlier detection method based on density. The main idea of the LOF determine outlier based on the density of each point and the density of neighboring points [24]. It depends on how isolated the object is with respect to the surrounding neighborhood. The flow chart of the global thresholding based on LOF is depicted in Figure 6.
We can calculate the comparison factor LOF via (12), and can also detect the outlier through it. The parameter MinPts specifies a minimum number of objects. If the factor value VOLUME 7, 2019  is close to 1, the density of p is close to the density of neighborhood, and p belongs to the same cluster. If the factor value is less than 1, the density of p is higher than the density of neighborhood, and p belongs to dense point cluster. If the factor value is greater than 1, the density of p is less than the density of neighborhood, and p is an outlier. And the specific threshold for detection depends on the maximum of LOF.

LOF MinPts
And the definition of LOF depends on the locally reachable density (lrd). Intuitively, the lrd of p is the inverse of the average reachability distance of p. Then the lrd of point p is defined as: The diagram of reachability distance of p is depicted in Figure 8. The reachability distance takes the maximum of d k (o) and d(p, o). And d k (o) is the k-th distance from o to the k-th closet points, d(p, o) is the distance from p to o. N MinPts (p) is k-distance neighborhood of p, which contains all points whose distance from p is not greater than d k (p).

C. RECONSTRUCTION
To further optimize the number of points for the reconstruction, the scheme in the paper includes the coarse reconstruction based on RANSC and the refined reconstruction based on power crust.
Power crust is presented by Amenta et al. [21], whose main idea is to calculate the Medial Axis (MA) and to obtain the object surface by inverse transformation. Power crust can be seemed as a way to extract triangular surface based on Delaunay. Medial axis is the central set of the largest polar sphere.
Assuming that S is sample point cloud on the object surface, the procedures of the power crust method can be described as follows [21]: 1) Constructing the Delaunay triangulation of S to generate triangulation mesh, and finding the voronoi vertices; 2) Choosing positive and negative poles ±p of each sample point; 3) Constructing the set of pole balls B p depicted in Figure 9. The set includes interior balls and exterior balls, and the center of the interior ball circle is the point on the MA. Then calculating the power diagram Pow(B p ), which is the weighted voronoi diagram; 4) Labeling its infinite outer pole with O and the opposite inner pole with I , and inserting both poles in the queue; 5) If the queue is nonempty, a labeled pole p is removed from the queue, and the Pow(B p ) of each unlabeled adjacent point q relative to p is examined. If the voronoi ball surrounding q intersects the voronoi ball of p at an angle of more than π/4, then q is set to the same label as p and is inserted in the queue. For sample point s, q is a pole of s. If the opposite point q of the pole o is unlabeled, o is set to the opposite label of q and is inserted into the queue; 6) Separating the cells of pole. Cells set labeled I is medial axis transform (MAT) and cells set labeled O is power crust.
The random sample consensus (RANSC) [25] method get effective sample data by calculating the mathematical model of the data parameter. The RANSC is usually used in the computer vision to calculate the concatenated transformation matrices. It selects the optimal fitting results iteratively via voting. This paper make surface fitting by RANSC and reserve interior point, and then obtain the coarse shape of object before the Power crust. RANSC can reduce the point   number to some extent, and the point number reduction depends on the distribution of points and the shape of the target.

A. ECHO SIMULATON
The simulation parameters of model shown by Figure 1 are set as follows: the water depth is 8 meters; the seabed tilt angle is 4 degrees; the object is about L = 2m length, 1 meter under water, at an angle of 45 degrees to sea level. The open angle of sonar scanning range up and down is ±22.5 degrees; the number of array elements is M × N = 48 × 48, and the distance between the array elements is 0.05 meters; the sample frequency is f s = 500kHz; the center frequency of the narrow-band pulse signal is f = 300kHz. The acoustic velocity distribution is constant. The unit of seabed volume absorption is taken as dB/m; seabed interface elect acoustic elastic half space; bottom scattering coefficient is µ = 0.027; open angle of source emission is [−45 • , 45 • ]. Figure 10 shows the seabed reverberation signal. Figure 11 shows the actual echo under the model. Figure 12 shows three-dimensional imaging results.

B. MULTIPATH INHIBITION
There is mirror image existed in the imaging result which is caused by reverberation and multi-path. This would generate more disturbance points and change the target contour features and position information, making it more difficult for the object reconstruction and the recognition operation. The false target outside the water is depicted in Figure 13. Since the sound ray have different paths back to the receiving array, some paths would lead to the target above the sea surface or below the seabed. Based on the inverse analysis of sound ray or prior information in the experiment, the false target outside the seawater can be removed by restricting the points within a certain geometric region. Figure 14 shows the result after multipath inhibition. It shows that the proposed method is effective for the mirror false target suppression and the outline of the line object is clearer. However, the point number are declined drastically.
To verify the efficiency of the model and the multi-path inhibition, the sea trial was carried under the same environment. The data was collected from MaShan acoustic laboratory in Wuxi, Jiangsu Province. The depth of water is 2.2m, and the direction range is 50m. Receiving array has 64 × 128 beams and then the open angle of direction is 80 degree, while the open angle of pitch is 30 degree. There are 2048 distance slices in the range 0-50 meters. The point cloud and reconstruction results are depicted in Figure 15. The experimental results are consistent with those obtained by using simulated echoes. Furthermore, the results of multi-path inhibition for the experimental data are depicted in Figure 16, in which the contour of the line object is clearer.

C. VISUALIZATION RESULTS
To show the validity of data processing scheme, tank experiment has been carried in the tank and the diagram is shown in Figure 17. The depth of sonar system is 2m and it is 5.3m from the noise elimination wall, the distance between object VOLUME 7, 2019    and sonar is 3.4m. The object is a metal circular tube with a certain camber, whose length is 2.05m and diameter is 0.03m, whole body with white cloth wrapped, end with screw thread, as depicted in Figure 18. The receiver array has 48 × 48 elements and has 128 × 128 beams, and has 2048 slices in the range of 0-50m. According to the beams diagram and transformation equation, the point cloud of the metal tube  in the tank experiment is depicted in Figure 19. Because of noise and multi-path caused by the tank wall, it is difficult to recognize the target from the point cloud.
This section compares the fixed threshold segmentation, the global threshold segmentation and the global threshold segmentation based on LOF. The value of fixed threshold is set as 0.9 multiples of maximum intensity. The cycle number is 100 and the predefined error is 0.5 in the global thresholding. Besides, the K in the LOF is selected automatically, which is 0.9 multiples of points number. And points higher than 0.03 multiple max LOF is outlier.
As shown in Figure 21, although the global threshold has less outlier than fixed threshold, it still has some isolated point which does not belong to the object. Consequently, the global threshold based on LOF presented in this paper has good performance at the cost of sharp reduction in the point number.
After the segmentation, two-step reconstruction is carried out by the RANSC and the power crust. The RANSC fitting is shown in by 22. From the first step, the basic contour of the object is emerged, and interior point around the plane calculated processed as the input of the power crust. The results of the Medial axis and power crust are depicted in Figure 23. After refined reconstruction, the contour and size of the object is clearer compared to the photo depicted in Figure 18.  In order to validate the stability and feasibility of scheme, the paper repeats the data processing 10 times to get the size and position information of the object in the tank experiment, and the statistical analysis between the output of the scheme and the real object are shown in table 1.
From the statistical analysis in table 1, comparing with the real sizes 2.05m and 3.4m, the length and the range of the object are 1.16m and 3.505m, which can be accepted and means that the results of the scheme is reliable. And the compression ratio in this scheme based on point number after  processing is the lowest, which means the faster speed of processing.

V. CONCLUSION
In this paper, the 3D beam forming is first used to obtain the raw point cloud of the environment. Then, a novel denosing and segmentation method is employed for the raw point cloud to realize multi-path inhibition and the segmentation. Finally, the object reconstruction is carried out via RANSC and power crust. The simulation and experimental results indicate that the model and the scheme mentioned in the paper could be trusted in shape and size compared with the real object, and the realization has low cost and high flexibility in the meantime. Moreover, compared to other methods, the proposed method can obtain a more satisfactory contour since the disturb points are delete and the two-step reconstruction is more reliable. VOLUME 7, 2019 In the future, to avoid discontinuities and holes in the target surface, we can try to use interpolation method to obtain dense point cloud.