A Novel Horizon Picking Method on Sub-Bottom Profiler Sonar Images

Traditional manual horizon picking is time-consuming and laborious, while automatic picking methods often suffer from the limited scope of their applications and the discontinuity of picked results. In this paper, we propose a novel method for automatic horizon picking from sub-bottom profiles (SBP) by an improved filtering algorithm. First, a clear and fine SBP image is formed using an intensity transformation method. On this basis, a novel filtering method is proposed by improving the multi-scale enhancement filtering algorithm to obtain clear horizons from an SBP image. The improvement is performed by applying a vertical suppression weighting term based on the form of logistic function, which is constructed by using the eigenvectors from the Hessian matrix. Then, the filtered image is segmented using a threshold method, and the horizon points in the SBP image are picked. After that, a horizon linking method is applied, which uses the horizon directions to refine the picked horizon points. The proposed method has been verified experimentally, and accurate and continuous horizons were obtained. Finally, the proposed method is discussed and some conclusions are drawn.


Introduction
Sub-bottom profilers are designed to obtain the images of underwater SBP for waterway dredging, marine scientific research, and other tasks [1][2][3][4]. Horizon picking from SBP images is one of the most crucial steps in these applications and is the basis for analyzing attributes of sediment interfaces and building a sedimentary database [5]. Manual or semi-automatic methods are often used for horizon picking from SBP images but are time-consuming and laborious. Thus, many scholars have studied automatic horizon picking methods.
Bondár [6] used an edge detector to find sediment interfaces. This detection method is efficient for horizon picking but does not take the lateral continuity of interfaces into account and often results in discrete horizons. Several efforts have been made to remedy this shortcoming; for example, Maroni et al.
[5] applied a multi-resolution method to pick horizons and provide smooth and continuous reflectors, and Theuillon et al. [7] used a Markov Random Field (MRF) to preserve the continuity of picked horizons by taking the neighboring information into consideration. These methods improve the continuity of picked horizons but are easily influenced by the intensity imbalance in an SBP image. Different interfaces in an SBP image have different mean intensities, and this causes difficulties in the selection of appropriate intensity threshold, as well in the accurate and continuous horizon picking. horizon picking. Fakiris et al. [4] presented a novel sediment interface enhancement approach to solve the intensity imbalance problem by using Haar-like characteristics. Nevertheless, an often-occurring problem is the inconsistent interface boundaries among pings, which is induced by discontinuous reflectors due to a lack of appropriate image smoothing methods.
Horizon picking combining complex attributes has also drawn much attention [8,9]. Since this method introduces cosine phase attributes, it is quite helpful in highlighting all reflections and picking horizons [10]. This cosine phase attribute can preserve the continuity of the interface. However, in many cases, sub-bottom profilers cannot provide full-waveform data [11,12], and the cosine phase cannot be obtained. Thus, the application of this method has been limited. Furthermore, the cosine phase data is influenced easily by noise. Local phase has been used to highlight reflections when lacking full-waveform data [10,13], but it is not robust enough, and some pseudo structures may be highlighted too.
Recently, Wang et al. [14] presented a method that applied the Frangi's filtering method, known as multi-scale enhancement filtering, to highlight the horizontal characteristics of sediment interfaces in SBP images. However, the enhancement may lose some details of dipping interfaces and cause damage to the original interface structures. Therefore, Wang's method is useful for horizon picking from monotonically varying horizontal interfaces, but not from complex ones.
These existing methods can obtain discrete horizons due to intensity imbalance and may perform poorly since they ignore the inherent characteristics of interfaces. To automatically obtain continuous and accurate horizons from SBP images, in this paper we propose a novel horizon picking method based on the improved filtering method. The remainder of the paper is organized as follows. The proposed horizon picking method is detailed in Section 2 and tested in Section 3. A concise and precise description of the experimental results are given in Section 4. The paper closes with discussions and conclusions in Section 5 and Section 6. The flowchart of the proposed method is given in Figure 1.

Building Fine SBP Images
The proposed approach starts with forming the SBP images. Generally, the intensity distributions of the received echoes are imbalanced. For example, the reflections from the seabed may be too strong while those from the sub-bottom may be relatively weak. Thus, it is unavoidable that the presentation of sub-bottom structures in an SBP image in terms of intensity variations is weakened when attempting to directly convert the intensities to 0~255 grayscale.
Thus, a transformation is proposed as follows:

Building Fine SBP Images
The proposed approach starts with forming the SBP images. Generally, the intensity distributions of the received echoes are imbalanced. For example, the reflections from the seabed may be too strong while those from the sub-bottom may be relatively weak. Thus, it is unavoidable that the presentation of sub-bottom structures in an SBP image in terms of intensity variations is weakened when attempting to directly convert the intensities to 0~255 grayscale.
Thus, a transformation is proposed as follows: where G(i, j) is the gray level corresponding to the intensity Int.(i, j), and Int.(i, j) represents the intensity value of the ith sample (row) of jth ping (column). ub is the upper boundary, and lb is the low boundary, where gray level range [lb, ub] corresponds to range [0,255]. It is essential to determine the values of ub and lb, as an appropriate parameter setting can highlight the interface structures. To determine ub and lb more effectively, the potential interface structures, which are very likely the sampling points with reflections from interfaces, are obtained first. Then, by calculating the mean µ 0 and standard deviation δ 0 of intensities of these potential interface structures, and setting ub as µ 0 + 3δ 0 and lb as µ 0 − 3δ 0 , a better SBP image can be built using Equation (1).
In the above, the key is how to obtain potential interface structures. At first, the mean and standard deviation, namely, µ and δ, of the intensities of a survey line are calculated. Then, setting [lb, ub] as [µ − δ, µ +3 δ], and using Equation (1), an initial SBP image is formed. The potential interface structures are obtained based on the initial SBP image through a threshold horizon picking method described in Section 2.5, and µ 0 and δ 0 are finally calculated using the intensities of the potential interface structures. The final built SBP is shown in Figure 2, from which the details of sediments can be distinguished well.
Remote Sens. 2020, 12 It is essential to determine the values of ub and lb, as an appropriate parameter setting can highlight the interface structures. To determine ub and lb more effectively, the potential interface structures, which are very likely the sampling points with reflections from interfaces, are obtained first. Then, by calculating the mean μ0 and standard deviation δ0 of intensities of these potential interface structures, and setting ub as μ0 + 3δ0 and lb as μ0 − 3δ0, a better SBP image can be built using Equation (1).
In the above, the key is how to obtain potential interface structures. At first, the mean and standard deviation, namely, μ and δ, of the intensities of a survey line are calculated. Then, setting [lb, ub] as [μ − δ, μ +3δ], and using Equation (1), an initial SBP image is formed. The potential interface structures are obtained based on the initial SBP image through a threshold horizon picking method described in Section 2.5, and μ0 and δ0 are finally calculated using the intensities of the potential interface structures. The final built SBP is shown in Figure 2, from which the details of sediments can be distinguished well.

Characteristics of SBP Sediment Interface
A novel filtering method is proposed by improving the multi-scale enhancement filtering algorithm to obtain fine horizons in an SBP image. The improvement is performed by constructing a vertical suppression weighting term based on the form of a logistic function using the eigenvectors from the Hessian matrix. Then, the filtered image is segmented using the threshold method, and the horizon points in the SBP image are picked.
Sub-bottom profilers generate vertical section images of ocean sub-bottom sediments based on the propagation of low-frequency pulses through the water and the sub-bottom. The sub-bottom measurement is carried out ping by ping. When the impedance changes are large, the echo reflections are relatively strong [5,10]. Thus, the acoustic impedance contrast interfaces in the sub-bottom can be considered to be SBP image pixels with high intensity [14,15].
From the above, it can be derived that: 1. A sediment interface has lateral continuity; therefore, the pixels with high gray level representing interfaces in adjacent pings are continuous along lateral directions. Then the sediment interfaces can be approximated as line-like structures [14]. 2. Adjacent pings have different attitudes, incident energy, and signal-noise ratio (SNR); therefore the intensity imbalance often occurs along pings in SBP images.

Filtering Line-Like Sructures
The multi-scale enhancement filtering method is suitable for the enhancement of line-like structures [14]. The multi-scale enhancement filter is designed based on the Hessian matrix and uses the eigenvalues of the Hessian matrix to describe the characteristics of targets [16]. As described by Frangi et al. [16] and Wang et al. [14], the Hessian matrix can capture the intensity structure of the local neighborhood, including line-like structures. In multi-scale enhancement filtering algorithm,

Characteristics of SBP Sediment Interface
A novel filtering method is proposed by improving the multi-scale enhancement filtering algorithm to obtain fine horizons in an SBP image. The improvement is performed by constructing a vertical suppression weighting term based on the form of a logistic function using the eigenvectors from the Hessian matrix. Then, the filtered image is segmented using the threshold method, and the horizon points in the SBP image are picked.
Sub-bottom profilers generate vertical section images of ocean sub-bottom sediments based on the propagation of low-frequency pulses through the water and the sub-bottom. The sub-bottom measurement is carried out ping by ping. When the impedance changes are large, the echo reflections are relatively strong [5,10]. Thus, the acoustic impedance contrast interfaces in the sub-bottom can be considered to be SBP image pixels with high intensity [14,15].
From the above, it can be derived that:

1.
A sediment interface has lateral continuity; therefore, the pixels with high gray level representing interfaces in adjacent pings are continuous along lateral directions. Then the sediment interfaces can be approximated as line-like structures [14].

2.
Adjacent pings have different attitudes, incident energy, and signal-noise ratio (SNR); therefore the intensity imbalance often occurs along pings in SBP images.

Filtering Line-Like Sructures
The multi-scale enhancement filtering method is suitable for the enhancement of line-like structures [14]. The multi-scale enhancement filter is designed based on the Hessian matrix and uses the eigenvalues of the Hessian matrix to describe the characteristics of targets [16]. As described by Frangi et al. [16] and Wang et al. [14], the Hessian matrix can capture the intensity structure of the local neighborhood, including line-like structures. In multi-scale enhancement filtering algorithm, while obtaining the Hessian matrix, the differentiation is defined as a convolution with derivatives of Gaussians: The original image is I(x). γ is usually set to unity, σ is the scale parameter, while x = (i, j) represents the location of a pixel (i and j represent the horizontal and vertical coordinates of a pixel). G(x, σ) is given by: The Hessian matrix corresponding to the pixel at x, is given by: where G ii (x), etc., are the second-order partial derivatives of G(x, σ). Let λ 1 , λ 2 be the eigenvalues of matrix H(x), λ 2 ≥ λ 1 ; a descriptor can be constructed using these eigenvalues. Three cases are considered: 1.
If both λ 1 and λ 2 are small, the analyzed region has approximately constant intensities.

2.
If one eigenvalue is high and the other low, there are small gray level changes in one direction and significant changes in the orthogonal direction, which indicates that there is a line-like structure. 3.
If both λ 1 and λ 2 are high, their gray level changes in both directions, which indicates the presence of a corner.
Frangi et al. [16] used the geometric characteristics derived from the Hessian matrix to construct descriptors suitable for line-like structures. The constructed descriptors are shown below: where R B is the line-like structure descriptor, and S is the intensity descriptor. When the local area is a line-like structure, R B has a very large value. Assuming that the sediment interfaces are the foreground and the non-interfaces are the background, S is low in the background because the eigenvalues are small in the background. This phenomenon helps us distinguish interface structures from the background. Then, a filter based on above two descriptors, R B and S, is designed as: where IM is the filtering result for the SBP image, σ represents the Gaussian scale parameter and c is a threshold that controls the sensitivity of the filter to the measures R B and S. The filtered result IM(σ) is obtained by calculating Equation (7). In Equation (7), R B and S are two matrixes and every element in these matrixes can be obtained according to eigenvalues of the Hessian matrix. The eigenvalues of every element are obtained by eigen-decomposition of the Hessian matrix in Equation (4).
Remote Sens. 2020, 12, 3322 5 of 20 By combining the two descriptors, a line-like structure filter can be built. β is set as 0.5, as suggested by Frangi et al. [16]. c depends on the grayscale range of the image, and value of half of the maximum Hessian norm has been proven to work in most cases [16,17].
It is worth noting that the second-order derivative of a Gaussian kernel has a similar structure as the interface area. The second derivative of a Gaussian kernel at scale σ generates a probe kernel that measures the contrast between the regions inside and outside the range (−σ, σ) in the direction of the derivative (as shown in Figure 3, in which parameter σ is set as 10). As Figure 3 shows, the second-order derivative of a Gaussian kernel, the range (−σ, σ), corresponds to the interface area, while the other intervals correspond to the non-interface area. Since convolution is conducted during multi-scale enhancement filtering, after the convolution operation in multi-scale enhancement filtering, a high value will be obtained at the locations of interface structures, and the interface structures are enhanced.
Remote Sens. 2020, 12,3322 5 of 21 of the derivative (as shown in Figure 3, in which parameter σ is set as 10). As Figure 3 shows, the second-order derivative of a Gaussian kernel, the range (−σ, σ), corresponds to the interface area, while the other intervals correspond to the non-interface area. Since convolution is conducted during multi-scale enhancement filtering, after the convolution operation in multi-scale enhancement filtering, a high value will be obtained at the locations of interface structures, and the interface structures are enhanced. Because (−σ, σ) corresponds to the interface area, the interval length, 2σ, also equals the thickness of the interface. A multi-scale method using a set of scale parameters is applied for the determination of the final scale and the final result. The final result corresponding to the largest response can be obtained using the equation below (the value set of σ0 is {σmin, σ1, σ2, ..., σmax}):  Figure 4b,d show the results of SBP images filtered using multi-scale enhancement filtering method. The reflections from sediment interfaces are enhanced and the other discrete ones are filtered. However, it is noticeable that not only the interface structures are enhanced, but so are some vertical structures, marked with rectangles in Figure 4.
As depicted before, the sub-bottom profiler obtains the interface structure along pings, thus the noise caused by the measuring instrument is characterized by a vertical distribution. Additionally, the imbalanced emission energies or the different measuring attitudes in pings also produce a vertical boundary (area pointed by the green arrows in Figure 4a). These vertical boundaries can also be enhanced after filtering. These enhanced vertical structures will induce the erroneous picking of the horizon. Moreover, biases can also be induced in determining the positions of horizons. All of these factors show the necessity of suppressing vertical structures.  Because (−σ, σ) corresponds to the interface area, the interval length, 2σ, also equals the thickness of the interface. A multi-scale method using a set of scale parameters is applied for the determination of the final scale and the final result. The final result corresponding to the largest response can be obtained using the equation below (the value set of σ 0 is {σ min , σ 1 , σ 2 , . . . , σ max }): Figure 4b,d show the results of SBP images filtered using multi-scale enhancement filtering method. The reflections from sediment interfaces are enhanced and the other discrete ones are filtered. However, it is noticeable that not only the interface structures are enhanced, but so are some vertical structures, marked with rectangles in Figure 4.
Remote Sens. 2020, 12,3322 5 of 21 of the derivative (as shown in Figure 3, in which parameter σ is set as 10). As Figure 3 shows, the second-order derivative of a Gaussian kernel, the range (−σ, σ), corresponds to the interface area, while the other intervals correspond to the non-interface area. Since convolution is conducted during multi-scale enhancement filtering, after the convolution operation in multi-scale enhancement filtering, a high value will be obtained at the locations of interface structures, and the interface structures are enhanced. Because (−σ, σ) corresponds to the interface area, the interval length, 2σ, also equals the thickness of the interface. A multi-scale method using a set of scale parameters is applied for the determination of the final scale and the final result. The final result corresponding to the largest response can be obtained using the equation below (the value set of σ0 is {σmin, σ1, σ2, ..., σmax}): Figure 4b,d show the results of SBP images filtered using multi-scale enhancement filtering method. The reflections from sediment interfaces are enhanced and the other discrete ones are filtered. However, it is noticeable that not only the interface structures are enhanced, but so are some vertical structures, marked with rectangles in Figure 4.
As depicted before, the sub-bottom profiler obtains the interface structure along pings, thus the noise caused by the measuring instrument is characterized by a vertical distribution. Additionally, the imbalanced emission energies or the different measuring attitudes in pings also produce a vertical boundary (area pointed by the green arrows in Figure 4a). These vertical boundaries can also be enhanced after filtering. These enhanced vertical structures will induce the erroneous picking of the horizon. Moreover, biases can also be induced in determining the positions of horizons. All of these factors show the necessity of suppressing vertical structures.  As depicted before, the sub-bottom profiler obtains the interface structure along pings, thus the noise caused by the measuring instrument is characterized by a vertical distribution. Additionally, the imbalanced emission energies or the different measuring attitudes in pings also produce a vertical boundary (area pointed by the green arrows in Figure 4a). These vertical boundaries can also be enhanced after filtering. These enhanced vertical structures will induce the erroneous picking of the horizon. Moreover, biases can also be induced in determining the positions of horizons. All of these factors show the necessity of suppressing vertical structures.

The Verticality Descriptor
Besides the line-like feature, interfaces have another feature type, namely, the reflectors may be non-vertical. Many scholars, for example, Maroni [5] and Theuillon [7], have taken the non-vertical characteristic into consideration to avoid wrong interface structure picking in the multi-resolution filtering method and MRF method.
As described above, using multi-scale enhancement filtering method directly may lead to the enhancement of wrong structures, namely the vertical structures, and produce wrong horizons during horizon picking. To overcome the problem, it is necessary to improve multi-scale enhancement filtering method by suppressing these vertical structures. In this paper, we propose a weighted vertical suppression by taking the non-vertical characteristics of interfaces into account. If we define V as the weighted vertical suppression term, then Equation (7) is rewritten as: V has a suppressing function on the vertical structures and preserves non-vertical structures. As described in Section 2.3, the multi-scale enhancement filtering method only uses the eigenvalue descriptors. As for vertical suppression, it is better to use the eigenvectors for determining whether a structure is vertical or not. In fact, the eigenvectors have been obtained already when calculating the line-like descriptors. Therefore, a verticality descriptor can also be constructed simultaneously by leveraging the information of the eigenvectors.
Assuming that v 1 and v 2 are the two eigenvectors of the interface area, v 2 , which corresponds to the smaller eigenvalue, points to the principal direction of a local area. We set v 2 = (L x , L y ). Then, a verticality descriptor is constructed as: where L x and L y are the coordinates of v 2 in the x and y directions, respectively, and ϕ is the verticality descriptor. The smaller ϕ is, the more vertical the structure is.

Definitions of the Key Points
After defining ϕ, a soft threshold function, such as the Sigmoid Function, is needed to better preserve non-vertical structures. ϕ and V are the independent and dependent variables of the soft threshold function. The complex form of the Sigmoid function is called the logistic function; it has the same features as the sigmoid function and is suitable for complex situations [18]. Thus, the logistic function is used in this paper. The general form of the logistic function is expressed as: where A, c and k are the parameters of the logistic function [18]. ϕ and V are related to x and y, respectively. As shown in Figure 5, the logistic curve is S-shaped and can be divided into five areas at four key locations depending on the speed of curve change. These are the flat area (−∞, a), the initially ascending area (a, b), the ascending area (b, c), the initially flattening area (c, d); and the flat area (d, +∞). d and a are defined as key points, where the weighted vertical suppressing V starts to play a powerful role in the filtering operation. Thus, the locations of d and a determine the effect of vertical suppression and must be defined accurately.
The locations of a and d are related to the derivatives of the logistic function, initially, the first, the second, and the third derivatives of the logistic function are calculated as: Then, as shown in Figure 5, the black points at both sides of the y-axis represent the points where the third derivatives equal to 0 and the second derivatives reach their maximums. Therefore, the coordinates of the two black points can be calculated using: Third, since the curvature is approximated by the second derivative, the curvature reaches its maximum at the two black points. In Figure 5, the two red circles represent the curvature circle of the two black points and the red points are the centers of red circles. The abscissa values of the red points are defined as the locations where the logistic curve starts increase and starts to become significantly flat, namely the values of a and d, respectively. Thus, it is necessary to determine the centers of the red circles. There are several steps to achieve this.
First, we assume that the red lines are tangent to the black curve, while the blue lines are perpendicular to the red lines. Then, the equations of the blue lines are: d and a are defined as key points, where the weighted vertical suppressing V starts to play a powerful role in the filtering operation. Thus, the locations of d and a determine the effect of vertical suppression and must be defined accurately.
The locations of a and d are related to the derivatives of the logistic function, initially, the first, the second, and the third derivatives of the logistic function are calculated as: Then, as shown in Figure 5, the black points at both sides of the y-axis represent the points where the third derivatives equal to 0 and the second derivatives reach their maximums. Therefore, the coordinates of the two black points can be calculated using: Third, since the curvature is approximated by the second derivative, the curvature reaches its maximum at the two black points. In Figure 5, the two red circles represent the curvature circle of the two black points and the red points are the centers of red circles. The abscissa values of the red points are defined as the locations where the logistic curve starts increase and starts to become significantly flat, namely the values of a and d, respectively. Thus, it is necessary to determine the centers of the red circles. There are several steps to achieve this. First, we assume that the red lines are tangent to the black curve, while the blue lines are perpendicular to the red lines. Then, the equations of the blue lines are: Second, taking the points of maximum curvature (black points) as the centers of the circles (marked with black dotted lines), and the maximum curvatures as the radiuses, the black dashed circles are drawn. The equations of these circles are: Finally, these circles intersect the blue lines at the red points, so the abscissa values of the red points, namely a and d, can be defined using: Thus, the key points a, d are defined as:

Determination of the Final Vertical Suppression Model
To build the vertical suppression model, the relation between a, d and parameter k need to be determined first. Initially, the vertical suppression model is as shown below, setting c and A in Equation (11) to 1 for simplicity: In Equation (19), a lower V value corresponds to stronger vertical suppression. The parameter k and x determine what kind of vertical structure is suppressed.
Second, parameter k is related to the key points a and d, and thus the relations between k and a, d should be formulated. Many experiments have verified that it is appropriate to consider that the effect of vertical suppression is significant when the verticality descriptor ϕ lies in the range [0 • , 5 • ], and is not significant at all in the range [45 • , 90 • ]. This implies that the key locations of the vertical suppression model are at 5 • (corresponding to a) and 45 • (corresponding to d). Then, through Equations (17) and (18) and setting c and A as 1, a and d can be revised as: where m is the offset item. The additional offset term is equivalent to moving the logistic function curve in the positive direction of the x-axis to ensure that the key points of the curve are positive. By adding the offset term, x in Equation (19) equals ϕ − m. Through Equations (20) and (21), the relations between k and a, d are formulated. After calculating k and m, given a vertical structure, we can determine to what extent it should be suppressed using Equation (19). Lastly, it is necessary to determine the value of k and m according to Equations (20) and (21). These two equations are nonlinear systems; thus, it is difficult to solve them directly. An intelligent optimization algorithm called the particle swarm optimization (PSO) method is used to deal with this problem [19]. Let D be the dimension of a particle, while X i = (x i1 , x i2 , x i3 , . . . , x iD ) T represents the current position of the ith particle, p i = (p i1 , p i2 , p i3 , . . . , p iD ) T represents the optimal position that the ith particle has reached and v i1 = (v i1 , v i2 , v i3 , . . . , v iD ) T is the ith particle's velocity. The evolutionary pattern of PSO algorithm is shown in Equation (22): where t represents the number of iterations. c 1 and c 2 are known as the cognitive and social parameters, and c 1 = c 2 = 2 is commonly used. r 1 and r 2 are random numbers in [0, 1]; p gd means the particle with the best position; ω is the inertial weight, which is set as 1. In this paper, the dimension of the particle D is 2, because we only have two parameters.
The cost function of the optimization method is shown in (23). f cost is defined as the fitness value of each particle. If p besti represents the position corresponding to the best fitness value of the ith particle, and g besti represents the position corresponding to the best fitness value of all particles in the particle swarm, then the PSO algorithm is the process of searching for g besti . Through iterative computation, k was finally calculated as 0.69, and m was calculated as 0.44.
Then, the final vertical suppression model is formulated as below: where V is the suppression value, and ϕ is the verticality descriptor calculated using Equation (10).

Horizon Picking and Linking
Horizon picking is conducted based on the filtered SBP image. The threshold method is often adopted in horizon picking. In this method, we assume that the gray level of pixels belongs to a interface area is higher than a given threshold [20]. For the picking, the following rule is adopted: where µ p and δ p are the mean and standard deviation of the pixel gray levels in one ping; G p (m) is the gray level of the mth sampling point in this ping; and h is the threshold, which is usually set to 2. If the gray level of a pixel satisfies Equation (25), then the pixel will be accepted as a horizon point. After detecting all pings, a new image is formed by setting the accepted horizon points to 255, and all others to 0. Therefore, a binary image is obtained, which is called the horizon image. Horizon points that are close to each other in a 3 × 3 neighborhood form a horizon segment, such as the one shown by the black and blue points in Figure 6. Horizon linking is necessary for obtaining long horizon segments and analyzing the physical properties of the horizon. The simplest method of horizon linking is to analyze the endpoint of each horizon segment in a neighborhood. In Figure 6, the blue point is the endpoint of a horizon segment, and the green point is the point that requires linking. The red arrow represents the direction information of the blue point. If a horizon point, such as the green point in Figure 6, lies on the orientation (θ), is closest to the blue point and the horizontal distance (D h ) between the blue point and the green point is smaller than a set value (this value can be determined by referring to the horizontal resolution) and the vertical distance (D v ) is smaller than 2σ in this area (σ is given by Equation (8), while 2σ corresponds to the thickness of the interface), the blue point can be connected with the green point.
Remote Sens. 2020, 12, 3322 10 of 21 been calculated during horizon linking. As described in Section 2.4.1, the direction information of every pixel has already been obtained through the calculation of the eigenvectors.

Experiment and Analysis
In this section, the performance of the proposed method was verified and evaluated using data acquired by different sub-bottom profilers in different areas. EdgeTech 3100P profiler (made by EdgeTech, West Wareham, the United States) and C-Boom profiler (made by C-Products, Hong Kong, China) were used in ShenZhen, and the surveyed areas are marked with red rectangle in Figure  7a. C-Boom was also used in Bohai Bay, and the surveyed areas are marked with red rectangles in Figure 7b. Parasound 70 (made by Teledyne Marine company, North Falmouth, the United States) was used in the South China sea, and the surveyed area was not shown in this paper to maintain secrecy.

Horizon Picking on EdgeTech 3100P Data in ShenZhen
An SBP survey was carried out in the water area of Shenzhen, China. The sediments in the study area mainly consist of silt and sand, and the water depth varies from 5 m to 15 m. SBP data recorded using an EdgeTech 3100P was used for the experiments. The sub-bottom profiler EdgeTech 3100P transmits an FM pulse that is linearly swept from 4 to 24 kHz, with a sampling interval of 0.032 ms.
At first, a line surveyed using the EdgeTech 3100P was extracted to test the proposed methodology. No full-waveform data is provided by the EdgeTech 3100P, and the provided envelope data was converted into an SBP image. The transformation method has been described in Section 2.1. When converting the received data into an SBP image (Figure 8), a low ub value in Equation (1) can result in the loss of interface details. As is shown in Figure 8, the red arrow points to the seafloor, while the yellow arrow points to the sub seafloor sediment interface. These two interfaces marked with arrows cannot be distinguished well from sediment interiors. This phenomenon makes horizon Traversing the horizon image from top to down and from left to right, all horizon points are processed, and a final horizon segment can be obtained. It is worth noting that the proposed method is different from the previously presented method [20], as the direction information needs not have been calculated during horizon linking. As described in Section 2.4.1, the direction information of every pixel has already been obtained through the calculation of the eigenvectors.

Experiment and Analysis
In this section, the performance of the proposed method was verified and evaluated using data acquired by different sub-bottom profilers in different areas. EdgeTech 3100P profiler (made by EdgeTech, West Wareham, MA, USA) and C-Boom profiler (made by C-Products, Hong Kong, China) were used in ShenZhen, and the surveyed areas are marked with red rectangle in Figure 7a. C-Boom was also used in Bohai Bay, and the surveyed areas are marked with red rectangles in Figure 7b. Parasound 70 (made by Teledyne Marine company, North Falmouth, MA, USA) was used in the South China sea, and the surveyed area was not shown in this paper to maintain secrecy.
Remote Sens. 2020, 12, 3322 10 of 21 been calculated during horizon linking. As described in Section 2.4.1, the direction information of every pixel has already been obtained through the calculation of the eigenvectors.

Experiment and Analysis
In this section, the performance of the proposed method was verified and evaluated using data acquired by different sub-bottom profilers in different areas. EdgeTech 3100P profiler (made by EdgeTech, West Wareham, the United States) and C-Boom profiler (made by C-Products, Hong Kong, China) were used in ShenZhen, and the surveyed areas are marked with red rectangle in Figure  7a. C-Boom was also used in Bohai Bay, and the surveyed areas are marked with red rectangles in Figure 7b. Parasound 70 (made by Teledyne Marine company, North Falmouth, the United States) was used in the South China sea, and the surveyed area was not shown in this paper to maintain secrecy.

Horizon Picking on EdgeTech 3100P Data in ShenZhen
An SBP survey was carried out in the water area of Shenzhen, China. The sediments in the study area mainly consist of silt and sand, and the water depth varies from 5 m to 15 m. SBP data recorded using an EdgeTech 3100P was used for the experiments. The sub-bottom profiler EdgeTech 3100P transmits an FM pulse that is linearly swept from 4 to 24 kHz, with a sampling interval of 0.032 ms.
At first, a line surveyed using the EdgeTech 3100P was extracted to test the proposed methodology. No full-waveform data is provided by the EdgeTech 3100P, and the provided envelope

Horizon Picking on EdgeTech 3100P Data in Shenzhen
An SBP survey was carried out in the water area of Shenzhen, China. The sediments in the study area mainly consist of silt and sand, and the water depth varies from 5 m to 15 m. SBP data recorded using an EdgeTech 3100P was used for the experiments. The sub-bottom profiler EdgeTech 3100P transmits an FM pulse that is linearly swept from 4 to 24 kHz, with a sampling interval of 0.032 ms.
At first, a line surveyed using the EdgeTech 3100P was extracted to test the proposed methodology. No full-waveform data is provided by the EdgeTech 3100P, and the provided envelope data was converted into an SBP image. The transformation method has been described in Section 2.1. When converting the received data into an SBP image (Figure 8), a low ub value in Equation (1) can result in the loss of interface details. As is shown in Figure 8, the red arrow points to the seafloor, while the yellow arrow points to the sub seafloor sediment interface. These two interfaces marked with arrows cannot be distinguished well from sediment interiors. This phenomenon makes horizon picking difficult. The proposed method provides a way to select a better transform model. The key point is that the picked potential interface structures are used to estimate the intensity range of the interfaces. Then, a better transform model is built by referring to the intensity range. Thus, the characteristics of all interfaces are more pronounced, as shown in Figure 9a. Both the seabed and the sub-bottom in Figure 9a have relatively high gray levels, which is helpful for horizon picking.
Then, the proposed filtering method with vertical structures suppression was applied to the SBP image filtering. In the filtering, the set of scale parameters was set as {5, 7, 9, 11, 13} according to the thicknesses of sediment interfaces. Through multi-scale parameters setting, thicker interfaces such as sandy and thin interfaces such as silty are enhanced, and some sharp points are smoothed.
It is obvious that there are vertical structures marked with blue rectangles in Figure 9a. These vertical structures may be caused by high emission energies and attitude changes. If we apply the multi-scale enhancement filtering method directly, the vertical structures in the blue rectangle of Figure 9a will be enhanced. The enhancement result of the structures in the blue rectangle in Figure 9a is shown in Figure 4b. These enhanced vertical structures will induce mistaken horizon picking. The proposed method deals well with this problem by adding the vertical suppression term. Through the proposed method, the vertical structures have been suppressed as shown in Figure 9b, while the real interface structures are well preserved. As shown in Figure 9a, the area marked by the red rectangle is an underwater tunnel. Using the proposed method, this artificial structure is distinguished from these vertical structures caused by the measurement environment. The result verifies the good performance of the proposed method in preserving interface structures, as shown in Figure 9b.
Then, the horizon picking process described in Section 2.5 is applied on the filtered image. During picking, h in Equation (25) is set as 2. After the picking, the horizon linking is adopted for the picked horizons. Some tiny gaps were repaired, and the final picked horizons are shown in Figure 9c, while in Figure 9d they are shown superimposed on Figure 9a. It can be seen that the picked horizons are visually consistent with the original profile. It is worth noting that the sandy interfaces composed of discrete pixels with high intensities are displayed continuously in Figure 9d. Therefore, the proposed method overcomes the shortcoming of the traditional intensity threshold method that can only pick out a series of discrete horizons. Thus, the result verifies the good performance of the proposed method in obtaining continuous horizons.
Remote Sens. 2020, 12, 3322 11 of 21 characteristics of all interfaces are more pronounced, as shown in Figure 9a. Both the seabed and the sub-bottom in Figure 9a have relatively high gray levels, which is helpful for horizon picking. Then, the proposed filtering method with vertical structures suppression was applied to the SBP image filtering. In the filtering, the set of scale parameters was set as {5, 7, 9, 11, 13} according to the thicknesses of sediment interfaces. Through multi-scale parameters setting, thicker interfaces such as sandy and thin interfaces such as silty are enhanced, and some sharp points are smoothed. It is obvious that there are vertical structures marked with blue rectangles in Figure 9a. These vertical structures may be caused by high emission energies and attitude changes. If we apply the multi-scale enhancement filtering method directly, the vertical structures in the blue rectangle of Figure 9a will be enhanced. The enhancement result of the structures in the blue rectangle in Figure  9a is shown in Figure 4b. These enhanced vertical structures will induce mistaken horizon picking. Then, the horizon picking process described in Section 2.5 is applied on the filtered image. During picking, h in Equation (25) is set as 2. After the picking, the horizon linking is adopted for the picked horizons. Some tiny gaps were repaired, and the final picked horizons are shown in Figure  9c, while in Figure 9d they are shown superimposed on Figure 9a. It can be seen that the picked horizons are visually consistent with the original profile. It is worth noting that the sandy interfaces composed of discrete pixels with high intensities are displayed continuously in Figure 9d. Therefore, the proposed method overcomes the shortcoming of the traditional intensity threshold method that can only pick out a series of discrete horizons. Thus, the result verifies the good performance of the proposed method in obtaining continuous horizons.
To quantitatively evaluate the horizon picking result, classic detection accuracy metrics, including Precision, Recall, F-measure, and Accuracy [10], were used to assess the picked results. Taking manual horizon picking result as the reference (shown in Figure 10), the automatically picking horizons in Figure 9c were evaluated on the basis that a horizon point is deemed to be correct if the vertical distance between it and its reference was smaller than seven pixels. The statistical results are To quantitatively evaluate the horizon picking result, classic detection accuracy metrics, including Precision, Recall, F-measure, and Accuracy [10], were used to assess the picked results. Taking manual horizon picking result as the reference (shown in Figure 10), the automatically picking horizons in Figure 9c were evaluated on the basis that a horizon point is deemed to be correct if the vertical distance between it and its reference was smaller than seven pixels. The statistical results are shown in Table 1. It can be seen that an F-measure value of 85% and an accuracy of 99% were achieved using the proposed method.  Table 1. It can be seen that an F-measure value of 85% and an accuracy of 99% were achieved using the proposed method.

Horizon Picking on C-boom Data in ShenZhen
The C-Boom sub-bottom profiler system, with a central frequency of 1.76 kHz, a sampling

Horizon Picking on C-boom Data in Shenzhen
The C-Boom sub-bottom profiler system, with a central frequency of 1.76 kHz, a sampling interval of 0.1 ms, and a data recording length of 90 ms, was also adopted in ShenZhen. The seabed sediments in the study area mainly consist of silty and sand, and the water depth varies from 5 m to 15 m.
The data provided by the C-Boom sub-bottom profiler is full-waveform data. First, envelope data was extracted from the full waveform data through the Hilbert transformation. The following processing steps were kept the same as before. The C-Boom profiler achieves deeper penetration compared with the EdgeTech 3100P, so deeper interfaces can be detected using this profiler, while the EdgeTech 3100P can detect more sub-bottom details [21]. Compared with the EdgeTech 3100P, the data collected using C-Boom has a low SNR. During this measurement, the data quality was reduced by poor measurement environment factors such as swift currents. Figure 11a shows this profile. In Figure 11a, there is much salt-and-pepper noise, which was caused by the poor surveying environment. This noise has high intensity and often leads to erroneous horizon points when adopting the traditional intensity threshold methods in horizon picking, which only consider the high-intensity characteristics of interfaces. The proposed method takes full advantage of the line-like structures of sediment interfaces and is thus robust to this noise. Figure 11a was filtered by the proposed method, and the result is shown in Figure 11b. It can be seen that the noise is almost completely removed. The result verifies the performance of the proposed method in picking the horizon from C-boom data.

Horizon Picking on C-boom Data in Bohai Bay
To test the effectiveness of the proposed method, another line was surveyed using C-Boom (Figure 12a) in Bohai Bay. The seabed sediments in the study area mainly consist of silty clay, silt, and muddy-silty clay, while the water depth varies from about 8 m to 30 m. Figure 12a was processed using the proposed method, and the continuous horizon (Figure 12c) was picked. The multiples cannot be distinguished from the first echoes, and thus the false horizons caused by the multiples are

Horizon Picking on C-boom Data in Bohai Bay
To test the effectiveness of the proposed method, another line was surveyed using C-Boom (Figure 12a) in Bohai Bay. The seabed sediments in the study area mainly consist of silty clay, silt, and muddy-silty clay, while the water depth varies from about 8 m to 30 m. Figure 12a was processed using the proposed method, and the continuous horizon (Figure 12c) was picked. The multiples cannot be distinguished from the first echoes, and thus the false horizons caused by the multiples are still picked. It is obvious that there are vertical structures in the SBP image ( Figure 12a, marked with blue Boxes). These vertical structures may be caused by high emission energies and overlay of multiple and primary echoes. When using multi-scale enhancement filtering method directly, the vertical structures will be enhanced. The mistakenly enhanced result of Box2 in Figure 12a is shown in Figure  4d. These mistakenly enhanced vertical structures will induce mistaken horizon picking.
The proposed method deals with vertical structures mistakenly enhancement well by adding the vertical suppression term. Through the proposed method, the vertical structures have been suppressed and the interfaces have been picked effectively, as shown in Figure 12c.

Horizon Picking on Parasound P70 Data in the South China Sea
To further verify the performance of the proposed method in deep-sea applications, an SBP line obtained using the Parasound 70 sub-bottom profiler in the South China Sea with the water depth ranging from 400 m to 2400 m was used to pick the sub-bottom horizons. Figure 13c shows the picked results. In the experiment, the SBP data was acquired using the secondary low frequency (4 kHz) produced by this profiler. It is obvious that there are vertical structures in the SBP image ( Figure 12a, marked with blue Boxes). These vertical structures may be caused by high emission energies and overlay of multiple and primary echoes. When using multi-scale enhancement filtering method directly, the vertical structures will be enhanced. The mistakenly enhanced result of Box2 in Figure 12a is shown in Figure 4d. These mistakenly enhanced vertical structures will induce mistaken horizon picking.
The proposed method deals with vertical structures mistakenly enhancement well by adding the vertical suppression term. Through the proposed method, the vertical structures have been suppressed and the interfaces have been picked effectively, as shown in Figure 12c.

Horizon Picking on Parasound P70 Data in the South China Sea
To further verify the performance of the proposed method in deep-sea applications, an SBP line obtained using the Parasound 70 sub-bottom profiler in the South China Sea with the water depth ranging from 400 m to 2400 m was used to pick the sub-bottom horizons. Figure 13c shows the picked results. In the experiment, the SBP data was acquired using the secondary low frequency (4 kHz) produced by this profiler. Through the proposed method, the sediment interfaces with clear line-like structures were picked well. The profiles marked with blue rectangles in Figure 13 show acoustic blanks and acoustic turbidities. In some areas in the SBP image, there are large acoustic impedance differences because of gas seeps or carbonate rock, which cause the acoustic blank phenomenon [22]. Gas accumulation can also induce acoustic turbidities, which may be caused by the reflection and scattering of the acoustic energy by myriads of in situ gas bubbles, causing dark smearing on the SBP image.
These phenomena disturb the interface structure and weaken the line-like characteristics of interfaces [22]. In these disturbed areas, only the main interfaces are picked. The phenomenon may be beneficial to distinguish the abnormal area from other normal interface structures for finding gas seeps or some other feature areas.

Results
The superposition diagrams of the horizon and the original SBP images show that the picked horizons are visually consistent with the original profiles. To quantitatively evaluate the horizon picking, classic detection accuracy metrics, such as Precision, Recall, F-measure and Accuracy, were used to assess the picked results. Taking manual horizon picking result as the reference, the automatically picking horizons were evaluated. A precision of 86%, a recall of 85%, an F-measure Through the proposed method, the sediment interfaces with clear line-like structures were picked well. The profiles marked with blue rectangles in Figure 13 show acoustic blanks and acoustic turbidities. In some areas in the SBP image, there are large acoustic impedance differences because of gas seeps or carbonate rock, which cause the acoustic blank phenomenon [22]. Gas accumulation can also induce acoustic turbidities, which may be caused by the reflection and scattering of the acoustic energy by myriads of in situ gas bubbles, causing dark smearing on the SBP image.
These phenomena disturb the interface structure and weaken the line-like characteristics of interfaces [22]. In these disturbed areas, only the main interfaces are picked. The phenomenon may be beneficial to distinguish the abnormal area from other normal interface structures for finding gas seeps or some other feature areas.

Results
The superposition diagrams of the horizon and the original SBP images show that the picked horizons are visually consistent with the original profiles. To quantitatively evaluate the horizon picking, classic detection accuracy metrics, such as Precision, Recall, F-measure and Accuracy, were used to assess the picked results. Taking manual horizon picking result as the reference, the automatically picking horizons were evaluated. A precision of 86%, a recall of 85%, an F-measure value of 85%, and an accuracy of 99% were achieved using the proposed method as shown in Table 1. The accuracy value was the highest because there were so many non-horizon points, and neither the manual method nor the automatic method could pick these sampling points. The high values of precision, recall and F-measure showed that the proposed method had a good performance in horizon picking.

Scale Parameter
The selection of scale parameter σ is important in filtering on line-like structures. σ is strongly related to the sediment interface thicknesses. Different sediment interfaces have different thickness; thus, a multi-scale method is needed for the selection of appropriate parameters. It is most important to determine the minimum scale parameter (σ min in Equation (8)). A larger σ min will help avoid some unnecessary tiny interface structures. For example, during sub-bottom surveying for channel dredging, these tiny sediment interface structures can be ignored. However, a smaller σ min is required when these tiny sediment interfaces are used for the analysis of the sub-bottom geological structures.
As described in Section 2.3, it is appropriate to set σ as the interface thickness. Therefore, we define σ min as the size of the smallest structure that we want. In theory, we do not know the thickness of small horizons in advance. To set σ min , the intensity threshold method is used for picking horizons and σ min is defined as half of the mean thickness of all horizons. Through the above, σ min was set to 5, and the horizons using this value are shown in Figure 9c. The proposed method was also used to filter the SBP image of Figure 9a with σ min equal to 9, and the result is shown in Figure 14a. The horizon picked from this filtered result is shown in Figure 14b. It can be seen that some tiny interfaces are filtered, and the main interfaces have been retained. When compared with the horizons shown in Figure 9c, there are some interfaces in Figure 14b that were not picked. The proposed method was also used to filter the SBP image of Figure 9a with σ min set to 3, and the result is shown in Figure 15a,b. Compared with Figure 14b, it can be seen that many tiny interfaces were maintained in Figure 15b. value of 85%, and an accuracy of 99% were achieved using the proposed method as shown in Table  1. The accuracy value was the highest because there were so many non-horizon points, and neither the manual method nor the automatic method could pick these sampling points. The high values of precision, recall and F-measure showed that the proposed method had a good performance in horizon picking.

Scale Parameter
The selection of scale parameter σ is important in filtering on line-like structures. σ is strongly related to the sediment interface thicknesses. Different sediment interfaces have different thickness; thus, a multi-scale method is needed for the selection of appropriate parameters. It is most important to determine the minimum scale parameter (σmin in Equation (8)). A larger σmin will help avoid some unnecessary tiny interface structures. For example, during sub-bottom surveying for channel dredging, these tiny sediment interface structures can be ignored. However, a smaller σmin is required when these tiny sediment interfaces are used for the analysis of the sub-bottom geological structures.
As described in Section 2.3, it is appropriate to set σ as the interface thickness. Therefore, we define σmin as the size of the smallest structure that we want. In theory, we do not know the thickness of small horizons in advance. To set σmin, the intensity threshold method is used for picking horizons and σmin is defined as half of the mean thickness of all horizons. Through the above, σmin was set to 5, and the horizons using this value are shown in Figure 9c. The proposed method was also used to filter the SBP image of Figure 9a with σmin equal to 9, and the result is shown in Figure 14a. The horizon picked from this filtered result is shown in Figure 14b. It can be seen that some tiny interfaces are filtered, and the main interfaces have been retained. When compared with the horizons shown in Figure 9c, there are some interfaces in Figure 14b that were not picked. The proposed method was also used to filter the SBP image of Figure 9a with σmin set to 3, and the result is shown in Figure 15a,b. Compared with Figure 14b, it can be seen that many tiny interfaces were maintained in Figure 15b.
Thus, an appropriate σmin is helpful to avoid obtaining a large number of tiny horizons just a few pings long, as well as missing some important horizons.

Vertical Exaggeration
The direction φ calculated in Section 2.4.1 is the direction of the SBP image structure, not that of the real sediment interface. Different mapping scales used in the formation of the SBP image can result in different picking results. In this paper, the SBP image is formed by arranging the sequential samples of each ping along the depth direction and the successive pings along the sailing direction with navigation information. Thus, the SBP image is stretched longitudinally compared to reality, a phenomenon called vertical exaggeration. Because of this, the interface structures may have a similar appearance to vertical structures, and thus be subject to vertical structure suppression.
This situation can be handled by adjusting the value of key point a (in Equation (20)) in vertical suppression, which in this paper was set to 5°-m (m is a parameter in Equation (20)). Decreasing a means that the vertical exaggeration is mitigated, which is helpful to preserve some nearly vertical interfaces. In most cases, setting a to 5°-m is appropriate, and this has been verified through experiments in different situations.

S and RB in the Multi-Scale Filter
Using Equation (9), S and , as well as V, can be obtained. Since the effectiveness of V was illustrated in the experiments, the performances of S and RB are described in this section. Figure 16a shows the image of S, and the structures of sediment interfaces have been filtered and displayed well because of the specific structure of the second derivative of a Gaussian kernel in the multi-scale enhancement filter (as described in Section 2.3). However, some sediment interfaces having low intensity cannot be enhanced in S, as shown in the red rectangle area in Figure 16a.
can better show the line-like structures even low intensity occurs compared with S. The image of RB is shown in Figure 16b, and the red rectangle area in this figure has been well displayed. However, in the blue rectangle marked area in Figure 16b, some pseudo structures having line-like structures can also be enhanced in RB. These outliers have low intensities because they are reflections of weak acoustic impedance interfaces. Thus, only combining RB and S, as well as V, can obtain a better filtered result (as shown in Figure 9b). Thus, an appropriate σ min is helpful to avoid obtaining a large number of tiny horizons just a few pings long, as well as missing some important horizons.

Vertical Exaggeration
The direction ϕ calculated in Section 2.4.1 is the direction of the SBP image structure, not that of the real sediment interface. Different mapping scales used in the formation of the SBP image can result in different picking results. In this paper, the SBP image is formed by arranging the sequential samples of each ping along the depth direction and the successive pings along the sailing direction with navigation information. Thus, the SBP image is stretched longitudinally compared to reality, a phenomenon called vertical exaggeration. Because of this, the interface structures may have a similar appearance to vertical structures, and thus be subject to vertical structure suppression.
This situation can be handled by adjusting the value of key point a (in Equation (20)) in vertical suppression, which in this paper was set to 5 • -m (m is a parameter in Equation (20)). Decreasing a means that the vertical exaggeration is mitigated, which is helpful to preserve some nearly vertical interfaces. In most cases, setting a to 5 • -m is appropriate, and this has been verified through experiments in different situations.

S and R B in the Multi-Scale Filter
Using Equation (9), S and R B , as well as V, can be obtained. Since the effectiveness of V was illustrated in the experiments, the performances of S and R B are described in this section. Figure 16a shows the image of S, and the structures of sediment interfaces have been filtered and displayed well because of the specific structure of the second derivative of a Gaussian kernel in the multi-scale enhancement filter (as described in Section 2.3). However, some sediment interfaces having low intensity cannot be enhanced in S, as shown in the red rectangle area in Figure 16a. R B can better show the line-like structures even low intensity occurs compared with S. The image of R B is shown in Figure 16b, and the red rectangle area in this figure has been well displayed. However, in the blue rectangle marked area in Figure 16b, some pseudo structures having line-like structures can also be enhanced in R B . These outliers have low intensities because they are reflections of weak acoustic impedance interfaces. Thus, only combining R B and S, as well as V, can obtain a better filtered result (as shown in Figure 9b).  Figure 17a shows the result of local phase calculated based on Figure 9a (the parameters used were suggested by [10]). Local phase image can highlight structures with low intensities; however, some pseudo vertical structures may also be enhanced. Furthermore, when using Equation (25) to detect horizon points, outliers can be induced into the final horizon results. For example, a vertical structure is shown in Figure 17a (red rectangle marks). This vertical structure induces outlier in the final picked horizon points (marked by red rectangle in Figure 17b). Moreover, in the blue rectangle area in Figure 17b, the sediment interfaces which should be line-like structures become arc-like structures because of the interferences of the highlighted pseudo structures. Compared with the result in Figure 9c, the result based on local phase image is not robust enough.   Figure 17a shows the result of local phase calculated based on Figure 9a (the parameters used were suggested by [10]). Local phase image can highlight structures with low intensities; however, some pseudo vertical structures may also be enhanced. Furthermore, when using Equation (25) to detect horizon points, outliers can be induced into the final horizon results. For example, a vertical structure is shown in Figure 17a (red rectangle marks). This vertical structure induces outlier in the final picked horizon points (marked by red rectangle in Figure 17b). Moreover, in the blue rectangle area in Figure 17b, the sediment interfaces which should be line-like structures become arc-like structures because of the interferences of the highlighted pseudo structures. Compared with the result in Figure 9c, the result based on local phase image is not robust enough.  Figure 17a shows the result of local phase calculated based on Figure 9a (the parameters used were suggested by [10]). Local phase image can highlight structures with low intensities; however, some pseudo vertical structures may also be enhanced. Furthermore, when using Equation (25) to detect horizon points, outliers can be induced into the final horizon results. For example, a vertical structure is shown in Figure 17a (red rectangle marks). This vertical structure induces outlier in the final picked horizon points (marked by red rectangle in Figure 17b). Moreover, in the blue rectangle area in Figure 17b, the sediment interfaces which should be line-like structures become arc-like structures because of the interferences of the highlighted pseudo structures. Compared with the result in Figure 9c, the result based on local phase image is not robust enough.

Applicability of the Proposed Method
The intrinsic characteristic of an SBP image is that the sediment interface extends in a particular direction and is displayed as a line-like structure on the SBP image. Then, the directions of these line-like structures can be used for helping the horizon picking. Thus, based on this characteristic, the proposed method achieved excellent performance. When dealing with the full waveform data, the envelope data is first derived and then the proposed method is applied to the converted envelope data. Thus, the proposed method can be applied to both the full-waveform data and the envelope data. However, the instantaneous phase information of the full-waveform data was not taken into consideration, and only envelope data was used. As described by Dossi et al. [8], the instantaneous phase information provides some advantages in accuracy and continuity compared with the envelope data in horizon picking. The proposed method is therefore more suitable when only the envelope data is provided, or the instantaneous phase data is very noisy.

Conclusions
In this paper, we propose a novel horizon picking method based on the improved multi-scale enhancement filtering algorithm. The developed new horizon picking method effectively overcomes the limitations of the method based on the original multi-scale enhancement filtering algorithm and successfully obtained high-quality SBP horizons. In the proposed method, good SBP images were first formed. The improved filtering method enhances line-like structures well and weakens the effect of noise. The suppression of vertical structures maintains the actual horizon, which allows accurate horizon picking and linking. The measures guarantee the fine and accurate horizon picked using the proposed method. The performance of the proposed method has been verified and evaluated using data acquired by different sub-bottom profilers in different areas. In future work, we will further extend the method and combine it with the deep learning framework to get better results and to realize real-time SBP horizon picking.