A Comprehensive Horizon‐Picking Method on Subbottom Profiles by Combining Envelope, Phase Attributes, and Texture Analysis

We propose a new method that considers the envelope, phase attributes, and texture analysis of the subbottom profile to automatically obtain continuous and accurate horizon picking. This method overcomes the shortcomings of the traditional methods in envelope and phase automatic horizon picking in terms of accuracy and continuity, and in efficiency when it comes to manual picking. Under the constraint of the envelope threshold from the envelope image and the texture segmentation from the phase image, the accurate envelope horizons are picked from the envelope image. Through mean filtering and template enhancement, the fine phase horizons are picked from the phase image. Finally, by combining the two kinds of picked horizons, more continuous, accurate, and finer horizons are obtained. The proposed method was tested on a real data set.


Introduction
The subbottom profile (SBP) plays a significant role in ocean engineering and marine science research (Fakiris et al., 2019;Maroni et al., 2001). The SBP image can be analyzed for waterway dredging, pipeline locating, and more. Different subbottom layers displayed by the SBP have significant applications in marine science, such as mechanical sediment behavior and analysis of geo-acoustic characteristics. Since picked horizons allow the interpreter to identify and correlate lithologic interfaces and sequence stratigraphic boundaries, horizon picking is a basis to identify sediment layers and to build a sedimentary database. (Shilts & Clague, 1992;Wunderlich et al., 2005).
Manual or semiautomatic methods are often used to pick horizons. Generally, strong echoes are reflected from the interface between two sediments form obvious sediment layers in the SBP image. Manual picking can easily find these layers in SBP images based on this characteristic. Manual picking is the classic method and is frequently adopted in the horizon picking; however, it is time consuming and depends on the experience of the interpreters.
Many scholars have developed automatic horizon-picking algorithms to improve the efficiency of horizon picking. Bondár (1992) used an edge detector to find sediment layers. This detection method is efficient in horizon picking, but it obtains discrete horizons because it ignores the lateral continuity of layers. Maroni et al. (2001) took the distribution of reflector intensity into consideration to set the picking threshold and applied the multiresolution method to pick horizons. It is easy to provide smooth and continuous reflectors using this method by using a filtered low-resolution image. In addition, it is feasible to obtain accurate locations and details of the layers by the intersection operation between the reflectors and the high-resolution image. Unfortunately, the method may suffer challenges when the reflectors are not composed of high-intensity pixels in the SBP image. Multiscale rotated Haar-like features filtering was used for horizon picking by Fakiris et al. (2018), which can be helpful in picking reflectors with low intensities. However, the picked reflectors cannot be continuous. Patel et al. (2008) proposed automatic seismic horizon-picking algorithm by considering the data as a height field, but the discontinuities often appear in the picked horizons for ignoring the lateral continuity of layers. Goldner et al. (2015) proposed a new tracking algorithm for picking 2-D seismic horizons based on the shortest paths in directed acyclic graphs. The tracking algorithm can balance the global and local information of seismic data to correctly track horizons but is time consuming (Lomask, 2007). The methods described above only use envelope data. In actuality, the cosine of the instantaneous phase (hereinafter referred to as the cosine phase) is also an important piece of information from the raw recorded data, and it can be used for forming the finer subbottom structure image. Therefore, an alternative method has been proposed that takes the cosine phase into consideration in the SBP horizon picking. Dossi et al. (2015) applied an attribute-based autopicking algorithm, that is, an algorithm that uses the cosine phase to automatically pick horizons. Using this method, the picked horizons have good continuity but also include many unimportant or useless horizons due to only considering the lateral phase continuity. Forte et al. (2016) improved Dossi's method to take the instantaneous amplitude (hereinafter referred to as envelope) information into consideration. The improvement meets problems in the determination of the appropriate threshold used for removing these unimportant horizons, but it provides a good strategy for obtaining the ideal horizons by the combination of the cosine phase and envelope. By computing the spatial organization of reflections (or any other events) in terms of their continuity and homogeneity, Zhao et al. (2017) introduced texture attributes successively to the Ground Penetrating Radar (GPR) profile interpretation to improve the identification of complex stratigraphic features. This research indicates that texture attributes can be used to identify continuous and homogeneity reflections. These methods mentioned above are used for both electromagnetic GPR and seismic data processing but they can be improved further in some aspects.
In this paper, we propose a comprehensive horizon-picking method by combining the envelope, phase attributes, and texture analysis to obtain continuous, accurate, and fine horizons. The remainder of this paper is structured as follows: First of all, the preprocessing of the SBP data and the acquisition of the envelope data and phase data are described. Then, these data are converted into the corresponding envelope image and phase image, respectively, and the characteristics of the two images are analyzed. After that, the comprehensive horizon-picking method, including the envelope horizon picking, the phase horizon picking, and the combination of both the horizon-picking results, is proposed. Actual applications verified the proposed method, and the performances of the proposed method are also discussed further. In addition, some conclusions and suggestions are drawn out of the experiments.

Extraction of Instantaneous Attributes
The most common attributes are instantaneous amplitude, phase and frequency (Taner et al., 1979). The instantaneous amplitude (namely envelope) is applied to evaluate the reflection energy, and the instantaneous phase can reflect the continuity of the layer (Barnes, 1996). The instantaneous frequency is the first derivative of the instantaneous phase, which can depict the changes of lithological characters (Barnes, 2007). These instantaneous attributes can be obtained by Hilbert transformation from the raw SBP echo signals. The recorded trace is treated like a real part of the complex trace, and the imaginary part can be obtained after the Hilbert transformation (Dossi et al., 2015). If a recorded signal is h(t), then where b h t ð Þ is the transformation result; g(t) is the complex trace; and A(t) is the envelope, that is, the reflection strength. The instantaneous phase θ(t) can be obtained by The recorded signal and its transformation result can also be represented as Then, the cosine phase, that is, cosine of the instantaneous phase can be calculated as follows:

SBP Images and Their Characteristics 2.2.1. The Envelope Image and Its Characteristics
The envelope data can be converted to a gray level according to the varying range of the envelope (Haralick & Shanmugam, 1973;Zhao et al., 2017). To display the details clearly, it is necessary to use more gray levels to describe the reflection. However, the envelope distributions of the received echoes are unbalanced. The reflection from the seabed is too large while the reflection from the subbottom is low. Thus, there is a precision loss when attempting to directly convert the envelope value to grayscale. Therefore, the mean and standard deviation of the envelope data are calculated, that is, μ and δ, and, to minimize the loss, μ + 3δ is set as the upper bound and μ − 3δ is set as the lower bound in the transformation from the envelope data to the envelope image. [μ − 3δ, μ + 3δ] corresponds to [0,255]. The transformation can be carried out as follows: where G (i,j) is the gray level corresponding to the envelope data Enr; (i, j) represents the ith sample (row) of jth ping (column) in SBP data; and Bound up and Bound low are μ + 3δ and μ -3δ, respectively. By the transformation model, an envelope image ( Figure 1a) is formed by the data of the surveying line.
An envelope is the most common attribute used in horizon picking (Ding et al., 2012;Zhao et al., 2018), but this attribute is sensitive to interference phenomena that induce lateral variations. When the acoustic wave propagates in a medium, it attenuates with the distance. A thicker overlying formation can cause a weaker envelope of the below formation's reflection. Therefore, lateral variations of layer thickness produce gradual lateral changes in the envelope below this layer (Taner et al., 1979). In addition, intrinsic attenuation, which is related to the properties of the materials, also affects the envelope. The envelope variations are also associated with the measuring situation, sediment types, and distributions of the mediums, for example, a strong envelope exists in sand or gas accumulation areas (Garcia-Gil et al., 2002). All of these bring difficulties to horizon picking just based on the envelope image, and can cause discontinuous picked horizons. However, horizon picking from the envelope image has its advantages in ignoring unimportant horizons whose envelopes are relatively low. The picking is easy to distinguish layer reflections from nonlayer reflections since the reflections are relatively strong when the acoustic waves reach large impedance changes.

10.1029/2019EA000680
Earth and Space Science 2.2.2. Cosine Phase Image and Its Characteristics As shown in equation (2), the instantaneous phase removes the envelope information, and therefore it will not be influenced by lateral variations and other issues affecting the envelope horizon picking (Taner et al., 1979). Due to its characteristic, the instantaneous phase can be used to improve the lateral continuity of the picked layers. The instantaneous phase, as defined in equation (2), is wrapped between 180°and −180°, and therefore it cannot be easily processed (Forte et al., 2016). This problem can be solved by using the cosine phase. The cosine phase is smooth because it removes the 180°discontinuities, as defined in equation (4) (Forte et al., 2016). Similarly, a cosine phase image ( Figure 1b) can be formed by converting the cosine phase data into 0-255 gray level. Besides, the cosine phase image displays clear texture features, such as the continuity and homogeneity of the gray level in the horizons (Figure 1b), which can be used to distinguish horizons from the cosine phase image as well as to pick horizons.
The cosine phase can highlight the lateral continuity of any coherent event. A ping may contain long record at a limited measurement depth. In the record, a few signals come from layers' reflections while many signals come from other interference. All of these highlighted reflections are picked when trying to pick the local maximums as the horizon points. Moreover, a single reflector is made by several phases. Therefore, it is necessary to filter out those phases that can be considered either unrelated to reflectors or distorted by interference or noise. It is difficult to distinguish the useful and useless reflections in the cosine phase image directly but is easy to diagnose them by their intensities in the envelope image because of their differences in envelopes. Therefore, comprehensive horizon picking is proposed in this paper by combining the envelope image and the cosine phase image.

Extraction of the Local Phase
Some subbottom profilers can only provide envelope data (Park et al., 2019). The local phase of a monogenic signal in this situation can be obtained and used for forming the local phase image to highlight the reflections, as shown in Figure 2b. The local phase in monogenic signal analysis is usually used for detecting linear structures and is fit for horizon picking (Felsberg & Sommer, 2002;Hidalgo-Gato & Barbosa, 2017). Felsberg and Sommer (2000) introduced the following filters in frequency domain: The spatial representation of H (H = (H 1 ,H 2 )) is the convolution kernel of the Riesz transform and it reads as follows: For a 2-D signal f(x 1 ,x 2 ), the Riesz transform is given by A signal f(x 1 ,x 2 ) can be combined with its Riesz transform according to Then, the multidimensional generalization of the analytic signal is obtained, which is called the monogenic signal. Then, a bandpass filter is applied for the Riesz transform. A good choice for the bandpass filter is a lognormal filter. The radial bandpass filter with the transfer function is used as follows: where w 0 is the normalized radius from center of frequency plane and f 0 is the center frequency of the filter and wl is the center wavelength expressed in pixel units. The bandwidth parameter sigma is usually set as 0.74 or 0.55 or 0.41 (corresponding to one, two, and three octaves). In this paper, sigma was chosen to be 0.55, which is a modest choice and would be appropriate for most situations. In addition, wl is determined by the thickness of the layer (the mean vertical distance of the two adjacent minimum points in the envelope data). Then, the following filter is used: . Then the monogenic scale space is constructed by a set of bandpass filters as follows as Furthermore, a new Riesz transform is given by Then, the local phase ϕ(x 1 , x 2 ) is given by Similarly, the local phase can be converted into the local phase image ( Figure 2b). The local phase image removes the influence of intensity variations with penetration depth and enhances all of the reflections (Hidalgo-Gato & Barbosa, 2017), and thus it has the same characteristics as the cosine phase image. The characteristics of the local phase imply that the cosine phase horizons can be replaced with the local phase horizons when the cosine phase data is lacking. The code for converting the envelope to the local phase was proved by Chris Bridge (Bridge & Noble, 2015) and can be downloaded from https://github.com/CPBridge/ monogenic_signal_matlab.

Envelope Horizon Picking
The horizon that is picked based on the envelope image is called the envelope horizon. The picking method keeps the peaks (local maximums) of a ping as the horizon points. Before the work, it is necessary to smooth the envelope image to enhance the correlation between neighboring pings and reduce the reflector interruptions. After the smoothing, to avoid the false subbottom structures caused by noise, only the peaks whose envelope or gray level are larger than a given threshold will be selected as the preliminary envelope horizon points. The final envelope horizons are obtained from these selected peaks by using the constraint of texture segmentation from the phase image.

Envelope Threshold Horizon Picking
A ping-mean filtering method, that is, replacing the gray level of a ping with the mean gray level of their neighbor pings, is applied to smooth an envelope image. The method can be depicted as follows: where the mean gray level G new (i, j) of jth ping and ith row is obtained by averaging the gray level of its left neighbor n pings and its right neighbor n pings. 2n + 1 should be no greater than the horizontal resolution (Pinson et al., 2008) (7 pings, about 3.5 m in this paper). In this paper, n was set to 3.
A new image is formed after the filtering. A threshold for the envelope horizon picking also needs to be determined. The peaks (local maximums) whose gray levels are larger than a given threshold will be selected. The gray level distributions of different envelope images are not exactly identical, but they are related to a Rayleigh distribution. Therefore, the threshold can be determined by a Rayleigh distribution (Maroni et al., 2001). The probability density function (PDF) of the Rayleigh distribution f(x) is shown in equation (16) (Rashwan, 2013) along with the corresponding cumulative distribution function, F(x) is given by equation (16) (Rashwan, 2013), then the boundary location can be obtained where 95% of the data distributed as the last expression in equation (16) shows.
where δ represents the scale parameter and x in the last expression of equation (16) corresponds to the location where F(x) equals 95%. Figure 3a shows three simulated Rayleigh distributions, and Figure 3b displays the gray level histogram of the envelope image and its fitting curve obtained by equation (17). The fitting curve is approximated to the PDF curve. In an envelope image, the pixels of the noise and the medium interior reflections make up 95% of the whole envelope image whose gray level is low and is distributed in [g 0 , g 0 + 2.45δ] while that of the reflections from the horizons is high and more than g 0 + 2.45δ (g 0 is defined in equation (17)). The characteristics imply that we can determine the threshold by finding the boundary location of the fitting curve of the gray level histogram where 95% of the data are distributed.
If f(g) is a fitting function of the gray level histogram of an envelope image, then where g 0 is the start position (the first nonzero value) of f(g), δ represents the scale parameter, and k is used as a parameter to adjust f(g) to coincide with the gray level. g m corresponds to the maximum of the distribution. g 0 and g m can be found by a traversal search along the x axis. As shown in Figure 3b, δ = g m − g 0 . Then, the envelope threshold can be given as g 0 + 2.45δ. By the threshold, these local maxima, whose gray levels are greater than the threshold, are picked as preliminary envelope horizon points in a ping. A preliminary envelope horizon image is constructed using these local maxima of all pings by setting the gray levels of these local maxima as 255 and the others as 0.

Texture Segmentation
Most of the horizon-picking methods are based on envelope data. However, some uninteresting points, such as these strong echoes from medium interior such as gas accumulation areas can be picked up when only considering characteristics of layers in the envelope. If the horizon picking focuses on the reflections from layers, then texture analysis provides a way of removing these uninteresting strong echoes (Haralick & Shanmugam, 1973). As mentioned in section 2.2.2, the texture features from the phase image can be used for discriminating between actual horizon points and nonhorizon points. The horizons points can be distinguished from the strong echoes through the characteristics of the continuity and homogeneity of gray level in layers.
There are many texture features, such as gray level co-occurrence matrix (GLCM) (Baraldi & Parmiggiani, 1995;Zhao et al., 2017), Gabor texture (Kamarainen et al., 2006), and Local Binary Pattern (LBP) (Ojala et al., 2002). Among these texture descriptors, the image segmentation using LBP features is easy to realize. In addition, LBP considers the relationship of the center pixel value with the surrounding ones, which is associated with the characteristics of the layer area and nonlayer area in a phase image. The LBP operator labels the pixels of an image by comparing the gray level of a pixel and those of its 3 × 3 pattern (Ojala et al., 2002). Then, the LBP of a pixel can be calculated as follows: where the center pixel's gray level is g c and the surrounding pixels' values are g p (p = 0,1,2, … 7). P = 8. s is the discrimination function, given as Based on the above, texture segmentation can be depicted as follows: 1. Divide the phase image into subareas of the same size. 2. Calculate the 3 × 3 LBP of every pixel and use the histogram (the number of bins is 8) of LBPs of a square subarea as a feature vector to describe the subarea, and then use the feature vectors of all the subareas to classify the phase image into two categories, that is, the layer parts with high gray levels and the nonlayer parts with low gray levels, through a k-means algorithm (Lomask, 2007). 3. Set the gray level of the layer part as 255 and that of the nonlayer part as 0.
To guarantee that the texture characteristics are efficiently found from the cosine phase image, the height of subarea is set as the mean vertical distance between two adjacent zero cross points in the cosine phase data to make the gray level of subarea have good consistency and the ability to distinguish the layer areas from the cosine phase image. The adjacent zero cross points are such that one is in the up neighborhood of the cosine phase peak and the other is in the down neighborhood, and both of their values are equal to 0. The adjacent zero cross points are closed to the positions of the preliminary envelope horizon points that were picked in section 3.1.1; thus, they can be found by referring to these positions. The height of the subarea is finally obtained by averaging the vertical distances between the two adjacent zero cross points in all of pings. It should be noted that the texture segmentation adopts the cosine phase image when both of cosine phase and the envelope attributes are provided, or when only the envelope attribute of the local phase image is provided. When texture segmentation adopts the local phase image, two adjacent minimum points are used to calculate the height of the subarea. Moreover, the aim of texture segmentation is to remove the uninteresting strong echoes. If the strong echoes are from interesting targets, then the texture segmentation can be ignored in the horizon picking, and the preliminary envelope horizons picked in section 3.1.1 are treated as the final envelope horizons.

Final Envelope Horizon
The final envelope horizon points are picked by the intersection operation applied between the preliminary envelope horizons and the texture segmentation. The intersection operation is as follows: where g (i, j) represents a pixel at the ith row and jth column and has gray level g. pAh and TS mean the preliminary envelope horizon and the texture segmentation, respectively. From equation (20), the revised preliminary envelope horizon points are obtained. Since the envelope maximums are close to the phase maximums in a subarea, the envelope horizon points should be near the upper boundary of the texture segmentation (Forte et al., 2016), and the vertical distances between the upper boundary and the phase maximums should not be larger than the height of the subarea used in the texture segmentation. Therefore, the envelope horizon points whose vertical distance to the nearest upper boundary is larger than the height of the subarea should be removed from the revised preliminary envelope horizon points. Then, the final envelope horizons are obtained.

Phase Horizon Picking
The phase horizon is obtained by determining positions of the peaks in the phase image (cosine phase image or local phase image). Peaks are picked from each ping as the horizon points and they are tracked to form the horizon segments.

Horizon Points Picking
To highlight layers and suppress noise, a template enhancement algorithm is put forward in this paper (Figure 4).
where R is the enhancing result of a cosine phase image Cimage; Size is the size of the enhancement template M. M is composed of three same-size parts. If the three parts are assigned −1, 2, and −1, respectively, then the fundamental form of M is [−1 … 2 … −1 …] T . If the center part of M corresponds to a layer, then Size decides the enhancement of characteristics of the layer. Therefore, the size of the median part of M is set as the mean vertical distance between two adjacent zero cross points in cosine phase data, which is defined in section 3.1.2. Specifically, Size is set to be three times the mean vertical distance.
A new cosine phase image is formed after the enhancement. The peak that is the local maximum in a small neighbor (usually three samples in a ping) will be accepted as a phase horizon point. Through a similar process, all of the cosine phase horizon points are obtained. In the algorithm, these peaks are enhanced and layer areas are highlighted by adopting the mean vertical distance to define the size of the enhancement template. Thus, the performance of horizon picking is improved.
The above stated process can be used to pick phase horizons from a cosine phase image. For a local phase image, the operation of image enhance can be ignored because the local phase image was enhanced when forming the image through the envelope data. The process of picking horizons from a local phase image is the same as that from a cosine phase image.

Phase Horizon Segment Tracking
After obtaining the phase horizon points, a phase horizon image is formed by setting 255 as the gray level of these horizon points and 0 for other points. Horizon tracking is performed for the phase horizon image to obtain the continuous phase horizons. The process is as follows: 1. Traverse each pixel of a row from left to right. 2. Come to a horizon point, if there is no horizon point in its 3 × 3 left neighborhood, then In Item 1 the current horizon point is set as the starting point, and an attempt is made to find next horizon point in the starting point's 3 × 3 right neighborhood. In Item 1 if there is no horizon point, then calculate the local orientation of the existing horizon points and find horizon points along the orientation ( Figure 5). If there is a horizon point along this orientation near the starting point and the distance between the starting point and the found point is smaller than a set value (this value can be determined by referring to the horizontal resolution) and the vertical distance is smaller than the mean vertical distances between two adjacent zero cross points, then the starting point can be connected with the found point, and they can be stored together in the same array. Each array is constructed as a horizon segment and is labeled with a number. Thus, every point in the same segment has the same label. The found point is set as a new starting point and (1) is repeated to find the next horizon point in the starting point's 3 × 3 right neighborhood. Otherwise, return to 1) and 2).
In Item 2 if there is a horizon point, then the horizon point is connected to the starting point and they are stored together in the same array. Additionally, the horizon point is set as a new starting point, and (1) is repeated to find the next horizon point in the starting point's 3 × 3 right neighborhood.
3. If there is a horizon point in the current point's 3 × 3 left neighborhood, continue to traverse this row until a horizon point is found and then return to 2). 4. Traverse all the rows in the phase horizon image to obtain the horizon segments.
After the above tracking, the horizon segments in the phase image are obtained. In Item 1, the local orientation needs to be calculated. As shown in Figure 5, the black points are the horizon points, the blue point is the current analysis point, and the green point is the found point. The local orientation can be obtained by designing a rectangle box centering on the current analysis point (the blue point). The length of the box m is defined as 1.2 times the lateral resolution (9 pixels in this study), and the width n is defined as 1.2 times the mean thickness of layers (5 pixels in this study). The initial orientation is set to 0°(the blue box shown in Figure 5), and then the pixels within the area with a width from i 0 − n/2 to i 0 + n/2 and the length from j 0 − m/2 to j 0 + m/2 are included if the location of the current point is (i 0 , j 0 ). The box is rotated with an angle interval of 5°to make the most of the points located in the box (the red box in Figure 5), and then the orientation of the long axis of the red box is also the local orientation.

Comprehensive Horizon Picking
The envelope horizon has an accurate location but poor continuity. The phase horizon compensates for the shortcomings of the envelope horizon but causes a large number of horizons marking strong and weak reflections alike (Barnes, 2007). The two kinds of horizons have strong complementarity; therefore, a comprehensive horizon-picking method is proposed as follows: 1. Obtain the envelope and phase (cosine phase or local phase) data and convert them to the corresponding images. 2. Pick the envelope horizons and the phase horizons as in sections 3.1 and 3.2. 3. Construct a set of the envelope horizons, EnvPoint, as {(i k , j k )} k = 1, 2…n . (i k , j k ) is the location of an envelope horizon point in the envelope image, and n is the number of envelope horizon points. 4. Obtain the corresponding phase horizon points using EnvPoint as a constraint through Gray phase I; J ð Þ ¼ 255; where Gray phase denotes the phase horizon image, (I, J) is the pixel at Ith row and Jth column in Gray phase , (i k , j k ) is the element in the EnvPoint (k = 1,2,3 … n.), and r is the same as height of the subarea in section 3.1.2. After the process given in section 3.2.2, the label of every phase horizon point can be visited automatically, and therefore the labels of the phase horizon points that satisfy equation (22) can be obtained. After removing duplicate labels, all the remaining labels are used for constructing a new set, Label as {label 1 , label 2 , … , label n }. Then, the phase horizon image is traversed, and the label of every phase horizon point is visited. Once a visited label belongs to Label, the corresponding phase horizon point is kept, and the horizon segments are constructed with the phase horizon points that have the same label. After traversing the whole phase horizon image, the phase horizon segments whose length is shorter than that of the segment in its vertical neighborhood (the size of the vertical neighborhood is Sr) are deleted. In addition, some tiny horizon segments are also ignored. The process of comprehensive picking is shown in Figure 6. Sr needs to be determined in the process. There are several different phase horizons close to the same envelope horizon point within the range of layer thickness because each reflector has several phases and the envelope horizon points have time shifting induced by filtering. Therefore, Sr is set as the layer thickness in this paper.

Experiment and Analysis
An SBP survey was carried out in Shenzhen, China. The seabed sediments in the study area mainly consisted of muddy sand and sand, and the water depth varied from 5 to 15 m. The C-Boom Subbottom 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 adopted, which has a maximum penetration depth of about 60 m in a soft seabed. 99 SBP surveying lines were recorded for a total of 479 km.
A line was extracted to test the proposed methodology. Before the horizon picking, frequency filtering was performed to the raw data to improve the signal-to-noise ratio as well as the image performance to be formed subsequently. Then, the envelopes and cosine phases were calculated, and the corresponding envelope image (Figure 7a) and cosine phase image (Figure 7b) were constructed. In Figure 7a, the strong gray level marked with ellipsoids and the weak gray level marked with rectangles can be found easily. Many factors may cause the change of the gray level of a layer, such as changes of the thickness of the overlying formation, the ship attitude, and the intrinsic attention of sediment, which leads to lateral variations of the gray level in the layers.
Compared with Figure 7a, the cosine phase image (Figure 7b) displays continuous and clear horizons. Furthermore, homogeneous distributions of gray level in the layers were achieved (the ellipsoids in Figure 6. Flowchart of the proposed comprehensive picking method.

10.1029/2019EA000680
Earth and Space Science Figures 7c and 7d). However, the image also enlarged the low-intensity reflections from the medium interior. These phenomena are shown in the red rectangles of Figures 7c and 7d.
As shown above, it is necessary to combine the envelope image and cosine phase image for horizon picking. The proposed comprehensive picking method was adopted. First, the envelope horizon points were picked by the constraint of envelope threshold. The envelope threshold was determined according to equation (17). The gray level histogram of the envelope image and its fitting curve are shown in Figure 3b. g 0 , which is the g where the first nonzero value occurs in the histogram, equaled 100; the g corresponding to the maximum of the distribution, g m , equaled 115; and the corresponding δ was 15. Therefore, the envelope threshold was defined as 136. Then, all the local maximums in each ping could be found by sliding the 5-pixel window along the corresponding column in the envelope image, and the maximums more than the threshold were set as preliminary envelope horizons ( Figure 8a). As shown, a lot of the reflections from the medium interior appear in Figure 8a (such as those in the red rectangles of Figure 8a). Texture segmentation was applied on the cosine phase image to remove these useless horizon points, and the segmentation was used to revise the preliminary envelope horizons. Figure 8b shows the segmentation result. The result shown in Figure 8b is rough and unclear relative to that in Figure 7a, which is caused by the mechanism of the texture segmentation. Even so, it can be found that these useless reflections from the medium interior (red rectangles in Figures 8a and 8b) disappear in the texture segmentation image. Figure 8a was revised after the texture segmentation by the intersection operation between the preliminary envelope horizon and the segmentation result. The revised result, that is, the envelope horizons, are shown in Figure 8c. Figure 8c is accurate and clear relative to Figures 8a and 8b, which shows that the texture segmentation is necessary and the proposed envelope horizon picking is efficient.
The threshold used in envelope horizon picking ensures that major layers are picked, but there still are many interruptions in the result. In addition, some horizons may be missed. Moreover, a long envelope horizon may be composed of several discrete horizon segments with different lengths. When applying a threshold for the lateral extension of the picked horizons in Figure 8a, the segments with short lengths are removed and the long horizon is broken up into several discrete segments. These problems show that it is necessary to combine it with cosine phase horizon picking to improve the horizon picking.
A 3 × 3 mean filtering and a 1 × 12 template enhancement were applied to the cosine phase image before the phase horizon picking. The local maximums, that is, the cosine phase horizon points, were found by sliding a 5-pixel window along each column in the image and the result is shown in Figure 8d. As shown, the horizons in Figure 8d are clearer and more abundant than those in Figure 8c, and the missing horizons in the envelope horizon picking were obtained in the cosine phase horizon picking. The result can be attributed to the cosine phase enhancement for all the reflections. Furthermore, too many useless horizons were also picked, which increases the complexity of the image interpretation.
The horizons are accurate but interruptions and missing items exist in the envelope horizon image (Figure 8c). In addition, there were clear and continuous horizons that exist as part of the useless horizons in the cosine phase image (Figure 8d). The proposed comprehensive picking method was applied that considered the complementarity of the two kinds of horizons. By combining Figures 8c and 8d, the final horizons were determined and shown in Figure 8e. A grand total of 287 s were required to process this surveying line (the code was written in Matlab). In Figure 8e, the outliers in the red rectangle of Figure 8a have been removed. In addition, comparing the ellipses areas in Figure 8c with the same area in Figure 8e that displays continuous horizon shows that the shortcomings on the horizon interruptions caused by envelope picking were overcome by the proposed method. Moreover, in the circular areas in Figure 8c, only a few layers were mapped on Figure 8c, whereas these layers were described in detail on the final result obtained by the comprehensive method. Finally, the final horizons and the envelope image were superimposed, as shown in Figure 9a (the SBPs were processed by gray stretching to achieve a better display). The other surveying lines were also processed by the proposed comprehensive horizon-picking method. Figures 9b-9f randomly show the results of the five lines. The picked horizons in Figure 9 accurately and clearly display the continuous horizon and describe the fine subbottom structure, which verifies the efficiency of the proposed horizon picking. In Figures 9e and 9f, multiple reflections have been removed that refer to twice the TWT (two-way travel time) of the seabed. In Figure 9c, some reflections with low gray level were not picked, which may be caused by the inaccurate envelope threshold or the lateral extension filtering. This can be improved by adjusting the envelope picking threshold or the filtering threshold.

Earth and Space Science
To quantitative evaluate the horizon picking, classic detection accuracy metrics, such as Precision, Recall, F measure and Accuracy, were used to assess the picked results. Precision denotes the ratio of sampling points correctly recognized as horizon points to the total number of sampling points recognized as horizons points; Recall means the ratio of sampling points correctly recognized as horizon points to the total number of horizon points that were picked manually; F measure is the weighted average of Precision and Recall; and Accuracy is the ratio of sampling points correctly recognized as horizon points, or correctly recognized as nonhorizon points, to total number of sampling points. Taking the horizons that were picked manually from Figure 7a as the reference (Figure 10), the horizons in Figure 8e were evaluated by the principle that the horizon point was diagnosed to be correct as the vertical distance between the diagnosed horizon point and its reference was less than five pixels. The four statistic accuracies are listed in Table 1. Precision was relatively low because many tiny layers were missing in manual picking results, while the automatic picking results were finer. Recall was relatively high since only a few horizons picked manually were missing in the automatically picked results. Accuracy was the highest because there were so many nonhorizon points, and neither the manual method nor the automatic method could pick these sampling points. These statistic parameters show that the proposed method had good performance in horizon picking. To justify the proposed method, five bore hole samplings were performed. Due to the specific engineering requirement, these samplings are sparse and shallow, and locate along the SBP line shown in Figure 9f. It can be seen from Figure 9f To further test the proposed methodology, the SBP data recorded by Edge Tech 3100P was also used for the experiments. The subbottom profiler, Edge Tech 3100P, has a maximum penetration depth of about 80 m in soft seabed and 6 m in hard sediment. It transmits an FM pulse that is linearly swept from 4 to 24 kHz, whose sampling interval is 0.032 ms. The Edge Tech 3100P Profiler carried out 96 SBP surveying lines at a total of 540 km and only the envelope data were provided by this profiler.
A surveying line was extracted (Figure 11a), and the main reflectors were picked by the envelope picking method (Figure 11b). Many interruptions marked with a red rectangle appear in the picked result. This phenomenon is like what is shown in Figure 8a and it cannot be avoided when only using the intensity data due to gradual lateral changes in gray level. According to section 2.2.3, the local phase image (Figure 11c) was formed, and the phase horizons ( Figure 11d) were picked out from the image. Figure 11d includes many unimportant horizon points. By combining the picked envelope horizons and phase horizons, the final horizons were obtained and are shown in Figure 11e. Compared with Figures 11b and 11d, the horizons have good continuity, accurate location, and fine distributions. The final horizons and the envelope image were

10.1029/2019EA000680
Earth and Space Science superimposed and shown in Figure 12a. In addition, the other surveying lines were also processed by the proposed comprehensive horizon-picking method. Figures 12b-12d randomly show the results of three lines.
Similar evaluation criteria to what was used in the last experiment were also adopted. The horizons picked manually from Figure 11a are shown in Figure 13 and were used as the reference, and the horizons shown in Figure 11e were assessed. The statistic parameters are listed in Table 2. Compared with Table 1, a similar result was achieved, which further shows that the proposed comprehensive horizon-picking method achieved good performance even when only the envelope data was proved by the profiler. Dossi et al. and Forte et al. (Dossi et al., 2015) systematically analyzed the cosine phase data and the envelope data in a synthetic signal obtained by the summation of two cosine functions with frequencies of 0.875 Hz and 0.625 Hz, and drew out the following characteristics:

Envelope Horizon and Phase Horizon
1. It is clear that the zero crossings of the raw recorded amplitude data and the cosine phase data coincide in time. Therefore, there is a one-to-one correspondence between the raw recorded amplitude peak and the cosine phase peak ( Figure 14). There are several phase peaks that belong to the same envelope maximum. Therefore, a maximum in the envelope does not necessarily correspond to a peak in the raw recorded amplitude data. Thus, the envelope cannot be used directly to obtain useful information about the raw recorded amplitude while cosine phase can. This phenomenon shows the advantage of the cosine phase, which also implies that the horizons picked from the cosine phase are more accurate. 2. The envelope can be used to quickly locate the strongest reflections in a data set. These characteristics show the complementarity between envelope data and cosine phase data. Due to the complementarity, the envelope data can be used to quickly locate the strongest reflections and the strongest reflections can be taken as constraints to construct phase horizons as well as to obtain accurate, continuous, and fine horizons by combining the envelope horizons and the phase horizons. Compared with horizon picking from single envelope data (Zhao et al., 2018)), the comprehensive horizon picking proposed in this paper introduces the phase horizon and thus improves the continuity of the picked horizons. Dossi et al. and Forte et al. also used the cosine phase data in horizon picking (Dossi et al., 2015;Forte et al., 2016). Compared with their work, the proposed method takes the distribution of the SBP image gray level into consideration and provides the criteria in the determination of envelope threshold. In addition, the texture segmentation is given for refining the envelope horizon obtained by threshold segmentation. Moreover, the enhancement and the orientation pattern ensure continuous tracking and accurate phase horizons. Furthermore, the local phase is introduced into the horizon picking as lack of cosine phase data, which expands the applications of the proposed method. Therefore, the comprehensive horizon picking proposed in this paper provides a better strategy for horizon picking. The above conclusions were well validated by experiments. Recently, a pure imaging processing method has been proposed by Wenbo Wang to pick horizons on seismic profile (Wang et al., 2018).Wang's method can be further improved to make it more appropriate for the SBP horizon picking when only envelope data are provided, and can also be used for forming a new comprehensive horizon-picking method by combining with the proposed phase horizon-picking method.  Figure 14. An example of the attribute analysis of a synthetic signal.

Comparison Between the Cosine Phase and the Local Phase
The raw recorded amplitude data in Figure 7a, which contains envelope data and cosine phase data, was extracted out for the following experiment to analyze the cosine phase and the local phase. A cosine phase image (Figure 15a) was formed by the cosine phase data, and a local phase image (Figure 15b) was constructed by the envelope data. Both of the two phase images removed the influence of intensity variations with penetration depth and enhanced all of the reflections. The cosine phase image showed a relatively clearer and sharper phase variations than the local phase image even though the local phase image was finer and clearer than Figure 7a. According to section 3.2, the cosine phase horizons and the local phase horizons were picked out from the two images and are shown in Figures 16a and 16b, respectively. The two sets of picked horizons were almost the same and display good continuity, which implies that the local phase data makes almost the same effects as the cosine phase data. The above analysis also provides a way to obtain accurate and continuous horizons by combining the local phase image under the situation where single envelope data provided.

Determination of the Envelope Threshold
To analyze the impacts of noise and multiple reflections on the gray level histogram of an envelope image, the gray level histograms of Figures 9c and 9f were calculated, as shown in Figures 17a and 17b, respectively. The histogram in Figure 17a has the fitting curve parameters g 0 and δ as 116 and 3, respectively. The corresponding fitting curve parameters in Figure 17b are 108 and 7. Multiple reflections exist in Figure 9f, while the noise and interference bring impacts, as shown in Figure 9c. The differences of these fitting curve parameters show that multiple reflections, noise, and interference influence the parameters of PDFs of different images but will not change their PDF curve shapes. This phenomenon also shows that the determination of the threshold in the envelope horizon picking by the PDF curve of the actual envelope image is appropriate.
The pixels of noise, the medium interior reflections, and weak reflections from the layer make up 95% of the whole envelope image. Through probability statistics, the 95% probability distribution corresponds to the gray level g 0 of the PDF curve boundary adding 2.45 times the scale parameter δ of the PDF curve. Therefore, the threshold was set to g 0 + 2.45δ. For this study, a proportion of 95% was adopted. If the proportion was decreased to 86.5%, then the threshold would become g 0 + 2.00δ and finer horizons would be picked out ( Figure 18a); however, interference would also appear in the picked horizons. On the contrary, if the proportion was increased to 99%, then the corresponding threshold would be g 0 + 3.00δ and only the main horizon in the image would be picked and some small layers might be ignored (Figure 18c). The experiment shows that the fineness of picked horizons depends on the envelope threshold and can be adjusted according to the actual requirements.

The Size of the Subarea in Texture Segmentation
In texture segmentation, the size of the subarea is decided by the mean vertical distances between two adjacent zero cross points. Sometimes, the size is not very appropriate due to acquiring the size by averaging all the vertical distances. A size that is too small or too large leads to the local feature being described inaccurately and some envelope horizon points to be deleted mistakenly as well as a discontinuity of the envelope horizon.
By referring to some multiscale features filtering method (Fakiris et al., 2019;Talukder & Casasent, 1998), a possible way to overcome the problem is to use different sizes of subareas to implement texture segmentation. By analyzing these different-scale segmentation results, the segmentation with the most consistency as the envelope horizons is desirable.

Limitations 5.5.1. Small Horizons
For this study, we adopted a 95% proportion to determine the envelope threshold for horizon picking. As shown in Figures 9 and 12, some small or weak reflections were not picked out. As discussed in section 5.3, these small or weak reflections can be obtained by setting a small envelope threshold to guarantee that the gray level of interesting reflections is larger than the threshold. Noise or interference may be introduced into the picked horizons through this operation, which can be weakened by a basic filtering algorithm, such as frequency filtering and amplitude recovery, to improve the signal-to-noise ratio as well as  the image performance and to find an appropriate threshold before the envelope horizon picking (Dossi et al., 2015).

Application of the Picked Horizons
This study focuses on the horizon picking from the SBP data but it does not consider further applications of the picked horizons. In practice, picking is related to its applications. For example, if we were interested in some small targets, an appropriate threshold is needed for finding the targets, and the polarity of reflections, as well as the determination reflection coefficients, are used for recognizing the target (Plets et al., 2009). In addition, by referring to the horizon distributions and the reflection coefficients, the picked horizons can also be used for the analysis of the subbottom layer. Therefore, the variations of reflection coefficients can be calculated in horizon picking and marked on the horizon image for stratigraphic interpretation.

Conclusions
The proposed comprehensive horizon-picking method takes the advantages of the envelope data and the phase data. This method provides the criteria for accurately determining the envelope threshold. Texture segmentation is adopted to refine the envelope horizon. Enhancement filtering and the orientation pattern are used to track the phase horizons, introduce the local phase as the lack of cosine phase data, and obtain continuous, accurate, and clear horizons automatically. Compared with the current horizon-picking method, the proposed method provides a more complete strategy in horizon picking.
The proposed method was applied to subbottom surveys. Two typical data sets, containing the envelope data and the cosine phase data surveyed by C-Boom and the envelope data obtained by the EdgeTech 3100P, were used for validating the proposed method. Almost all the horizons were picked continuously and accurately at a given threshold. Furthermore, the proposed method was also efficient in picking horizons from other SBP data.
There are two aspects of this research that should be further developed in future work. One is to improve the proposed method by referring to the application requirements of the picked horizons. The other is in applications of the picked results, such as finding the interesting targets in the subbottom and estimating the reflection coefficient as well as the distributions of the subbottom sediment types and so on.