Ship Wake Detection in SAR Images via Sparse Regularization

In order to analyze synthetic aperture radar (SAR) images of the sea surface, ship wake detection is essential for extracting information on the wake generating vessels. One possibility is to assume a linear model for wakes, in which case detection approaches are based on transforms such as Radon and Hough. These express the bright (dark) lines as peak (trough) points in the transform domain. In this article, ship wake detection is posed as an inverse problem, which the associated cost function including a sparsity enforcing penalty, i.e., the generalized minimax concave (GMC) function. Despite being a nonconvex regularizer, the GMC penalty enforces the overall cost function to be convex. The proposed solution is based on a Bayesian formulation, whereby the point estimates are recovered using a maximum a posteriori (MAP) estimation. To quantify the performance of the proposed method, various types of SAR images are used, corresponding to TerraSAR-X, COSMO-SkyMed, Sentinel-1, and Advanced Land Observing Satellite 2 (ALOS2). The performance of various priors in solving the proposed inverse problem is first studied by investigating the GMC along with the $L_{1}$ , $L_{p}$ , nuclear, and total variation (TV) norms. We show that the GMC achieves the best results and we subsequently study the merits of the corresponding method in comparison to two state-of-the-art approaches for ship wake detection. The results show that our proposed technique offers the best performance by achieving 80% success rate.


I. INTRODUCTION
A CCURATE characterization of sea surface condition is not only important in isolation but also in the detection and characterization of ship wakes. These provide key information for tracking (illegal) vessels and are also useful in classifying the characteristics of the wake generating vessel. Until recently, one of the main factors hampering research into sea surface modeling was the lack of data of sufficiently high resolution (pixels need to be typically smaller than a few meters) and accuracy. SAR technologies have however shown remarkable progress in recent years, and the availability of remotely sensed data of the Earth and sea surface is continuously growing. Several European missions (e.g., the Italian COSMO-SkyMed, the German TerraSAR-X, or the British NovaSAR mission) have developed a new generation of satellites exploiting SAR to provide spatial resolutions previously unavailable from space-borne remote sensing. In all SAR images, bright areas correspond to high radar cross section per unit area and dark areas to low radar cross section. When imaging sea surface, high returns are caused by enhanced surface roughness or large wave steepness, either via specular reflection or Bragg scattering. On the other hand, specular reflections are usually caused by breaking waves. Sea waves become visible in SAR images due to the Bragg scattering of SAR electromagnetic waves by small-scale capillary and gravity-capillary waves that propagate on the surface of larger waves [1]. These include swells, coherent Kelvin waves, random sea waves, as well as their mixture [2], [3]. SAR images of moving ships exhibit some characteristic patterns that are directly determined by different ship wake formations. These are typically considered to fall in one of three categories: 1) turbulent wakes; 2) surface waves created by ships; and 3) ship-generated internal waves. Ship-generated surface waves can, in turn, be split into two subcategories [4], [5], the first being the short (centimeter scale) waves, while the second includes the (decameter scale) waves forming the classical Kelvin wake system [6]. The former is observed in SAR images through the Bragg scattering mechanism and appear as bright, narrow V-wakes due to the resonant interaction of the transmitted radar waves and ocean surface waves [7]. The two-scale composite model has been employed to simulate SAR images of rough sea surfaces with embedded Kelvin wake structures in [8], while Fujimura et al. [9] performed a validation study for simulated ship wakes via computational dynamics.
Since ship wakes can be modeled as linear structures, corresponding detection methods are mostly based on linear feature extraction approaches, such as the Hough or Radon transforms, both of which create bright peaks in the transform domain for bright lines and troughs for dark lines. Due to its high computational cost, the Hough transform has attracted less interest than the Radon transform [10]. Thanks to the lower computational complexity of its inverse, the Radon transform is widespread in ship wake detection applications and has been first utilized by Murphy [11]. The Radon transform has, however, a couple of drawbacks, e.g., bright pixels belonging to ships may cause false detections [12]. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ To address this, enhancing the Radon domain information is a common practice in the literature. Rey et al. [13] have proposed a method that combines the Radon transform and the Wiener filtering to increase the detectability of the peaks in the Radon domain. Tunaley [14] used a method that restricts the search area in the Radon domain. Eldhuset [15] proposed an automatic ship and wake detection method, whereby the detection performance is characterized by the number of lost and false wakes (wakelike features). Based on ship beam and speed estimation, Zilman et al. [16] have applied an enhancement operation to the Radon transform of the observed (noisy) image. Courmontagne [17] used a method based on a combination of Radon transform and stochastic matched filtering for wake detection. Kuo and Chen [18] have proposed a ship wake detection method using wavelet correlators. Zilman et al. [8] proposed a SAR image simulator, including ship wakes, and studied the performance of their ship wake detection method previously published in [16]. Recent studies based on L 2 regularized logistic regression [19] and low-rank plus sparse decomposition [20] also addressed the ship wake detection problem in SAR images. Graziano et al. [10], [12], [21], [22] have proposed wake detection methods that deal directly with the noisy image without performing any preliminary enhancement. The authors have first created ship-centeredmasked image tiles and performed a restricted area search in their Radon representation. The restricted area is the area that lies between two sine waves, the peak points of which have been selected using real ancillary data, local map, ship clusters, and the local traffic statistics.
The inverse problem formulation for line detection has been first proposed by Aggarwal and Karl [23]. The main advantage of formulating line detection as an inverse problem is the subsequent use of a regularization framework, which allows the incorporation of prior information about the object of interest. Anantrasirichai et al. [24] have further investigated the inverse problem formulation for B-line detection in lung ultrasound images. Although the inverse problem solution creates an enhanced image, this step is different from operations like despeckling, which require statistical assumptions [25], [26] or statistical model selection [27] and is well studied in [28]- [32].
In this article, we propose a novel approach to ship wake detection in SAR images, which is based on an inverse problem formulation. The main contribution of this article is to propose an innovative approach based on sparse regularization to obtain the Radon transform of the image, in which the linear features are enhanced. The solution to the inverse problem involves Bayesian methodology, which leads to a MAP estimator. Our proposed cost function uses the GMC penalty of Selesnick [33] as regularization term and investigates its merit in comparison to the TV, nuclear, L 1 , and L p norms. The use of the GMC penalty demonstrates the advantages of nonconvex sparse regularization while allowing the cost function to remain convex. We use SAR images of the sea surface from different sources, including TerraSAR-X, COSMO-SkyMed, SENTINEL-1, and ALOS2 to test the performance of the proposed method. For the ship wake detection step, we use a method that performs detection in the Radon domain as proposed in [10] and [22]. In this study, contrary to [10] and [22], the proposed method detects ship wakes, independently of real information about the ships and the environment, by selecting the required parameters directly from the observed SAR images. MAP estimates for the images are obtained using the FB method and the TwIST algorithm [34].
The rest of this article is organized as follows: the image formation model for ship wakes identification, the MAP estimation, GMC regularization, and priors employed are discussed in Section II. The detection algorithms and SAR data sets are presented in Sections III and IV, respectively. Experimental studies and results are provided in Section V. Section VI concludes this article with a brief summary.

A. Ship Wakes and Image Formation Model
As discussed briefly in Section I, a moving ship in deep sea typically creates three different types of wakes. The central dark streak is called turbulent wake. Generally, this central dark streak is surrounded by two bright arms called narrow Vwake, which lies either side of the turbulent wake within the half-angle from 1.5 • up to 4 • [10], [16]. Finally, two outer arms are known as Kelvin wake and limit the signatures of the moving ship on each side of the turbulent wake with a maximum half-angle of 19.5 • . The Kelvin wake half-angle can be somewhat smaller in real SAR images, e.g., between 10 • and 19.5 • [2], [8].
For the purpose of this study, we consider each arm of narrow-V and Kelvin wakes as a separate wake and refer to them as the first and second wakes in the order of their detection. Consequently, we attempt to detect five linear wake-corresponding structures, which are the turbulent wake, the first and second narrow V, and the first and second Kelvin wakes. In most cases, not all five wakes are however visible in SAR images as is generally the case with one of the narrow V and/or Kelvin wakes. In Fig. 1, two example SAR images are depicted. In particular, in the upper left image, all five wakes are visible, whereas the one on the upper right does not have a visible Kelvin wake.
Since we model ship wakes as linear features, the SAR image formation can be expressed in terms of its Radon transform as is the image in Radon domain, W refers to the additive noise, and the operator C = R −1 is the inverse Radon transform. X (r, θ) represents lines as a distance r from the center of Y and an orientation θ from the horizontal axis of Y .
In image processing applications, to compute an integral of the intensities of image Y (i, j ) over the hyperplane, which is perpendicular to θ leads to the Radon transform X (r, θ) of the given image Y . It can also be defined as a projection of the image along the angles, θ . Hence, for a given image Y , the general form of the Radon transform where δ(·) is the Dirac-delta function. The inverse Radon transform (Y = C X) of the projected image X can be obtained from the filtered backprojection [35] algorithm as where v is the radius in Fourier transform, and F [·] and F −1 [·] refer to the forward and inverse Fourier transforms, respectively. In this article, we use discrete operators R and C as described in [36].

B. Bayesian Inference in Inverse Problems
Assume that the aim is to extract information about the unknown signal x given an observation signal y. Then, suppose a statistical model which defines the relation of y to x can be expressed. This statistical model is referred to as the likelihood p(y|x). Obtaining x from y might contain numerous uncertainties and thus extracting information belonging to x would be ill-posed and problematic. The Bayesian inference framework helps to reduce these uncertainties by employing knowledge on x, namely, the prior of x, which can promote structural properties such as sparsity.
Hence, incorporating this prior in conjunction with the observed statistical model creates the knowledge on x given y, namely, the posterior distribution p(x|y) via Bayes' theorem where the denominator p(y|x) p(x)dx is the marginal likelihood p(y), which is not related to x and assumed to be constant according to x. Since the aim is to extract x, we can write the unnormalized posterior distribution as Generally, posterior distributions are assumed to be log-concave as where F(x) is a convex function. Estimating x directly using the posterior density can be intractable, especially for high-dimensional cases. Hence, the MAP estimator maximizes the posterior p(x|y) to obtain a point estimatex. The MAP estimator for the unknown x iŝ In cases when F(x) is not convex, i.e., the posterior will no longer be log-concave, the minimizer will not be convex but nevertheless can be solved via proximal operators [37].
According to the SAR image formation defined in (1), the posterior distribution of the desired line image X with respect to the observed noisy SAR image Y can be written using (6) as Assuming an independent identically distributed (i.i.d.) standard normal noise case, the likelihood distribution p(Y |X) can be expressed as In this article, we assume priors to be of exponential form as p(X) ∝ exp{−λψ(X)}. Thus, the cost function in (8) becomes where λ is the scale parameter, namely, the regularization constant.

C. Generalized Minimax Concave Penalty
The GMC regularization has been proposed by Selesnick [33] as a sparsity enforcing penalty in inverse problems. It exploits the advantages of using nonconvex penalties as well as preserving the convexity of the cost function. It is based on L 1 -norm and the generalized Huber function.
In particular, the MC penalty in the univariate case is The relation between the MC penalty and the Huber function s(t) can be written as For a scalar, b = 0, the scaled MC penalty, ψ b (t) can be written as scaled Huber function, s b (t) as All these definitions are based on the univariate case, but they can be adapted to the multivariate case. If we assume a scaling matrix, B, then the generalized Huber function, S B (t) can be defined as [33] Combining (14) and (15) gives the GMC penalty function as The scaling matrix B should be selected in relation to the inverse Radon operator C to provide the convexity of the cost function as where λ 1 is the scale parameter of the GMC prior and γ is a parameter, which controls the nonconvexity. Note that for 0 ≤ γ ≤ 1, B ensures the convexity of the cost function. The nominal range of 0.5 ≤ γ ≤ 0.9 should be used for better performance [33]. Thus, the GMC sparse prior can be written as As GMC regularization does not have an explicit formulation, the minimization problem with the cost function in (11) can be solved using proximal operators. Thus, we rewrite the cost function as which leads to a minimax optimization problem The solution to this problem can be obtained using the FB algorithm [33]. The FB algorithm for GMC regularized cost functions requires only a couple of computational steps and soft-thresholding, which is the proximal operator for L 1 -based regularization. The FB algorithm to solve (21) is given in Algorithm 1. For the purpose of this study, we chose the maximum number of iterations MaxIter = 1000, which was experimentally set. The algorithm stops when the error term reaches 10 −3 . This error term (i) is calculated for iteration i as

D. Alternative Sparse Priors
In this article, we also investigate other types of priors, which are known to be sparse and common in optimization problems. The first of which is Laplace (or L 1 norm) prior, which is log-concave and given by  9: 10: v (i) = soft(u (i) , μλ) 11: i + + 12: while (i) > 10 −3 or i < MaxIter We further investigate the (nonconvex) L p -norm-based prior, which is where 0 < p < 1.
Another important type of prior is the TV norm, T V (·), which can be defined as [38] T V (X) = ∇ X 1 (27) where ∇ is the 2-D discrete gradient operator. Hence, the TV norm-based prior is given by The nuclear norm is an important type of prior for data with low-rank and low singular values. In particular, X * is the nuclear norm of X and is defined as the sum of its singular values [38]. Here, we define the nuclear norm prior as Replacing ψ(X) in (11) with the previously defined priors determines different MAP estimators aŝ where k = 2, 3, 4, 5.
Minimizations of (32) are subsequently carried out with the TwIST algorithm [34]. The TwIST algorithm is given in Algorithm 2 where MaxIter and the stopping criterion are the same as in the GMC case.
The operator λ k (·) in Algorithm 2 is a shrinkage/thresholding/denoising function. In this article, we used shrinkage/thresholding operators in minimizers with L 1 and L p for λ k (·) as also used in [24] Algorithm 2 TwIST Algorithm 1: Input: Ship-centered-masked SAR image Y 2: Output: Radon image X 3: Set: λ k > 0 and α = 1.96 (as in [34]) 4: Set: X (1) = λ k (X (0) ) and i = 1 5: do 6: Obtain λ k (w (i) ) for ψ k 8: i + + 10: while (i) > 10 −3 or i < MaxIter and the denoising operators for the TV and nuclear norm priors where prox λ k ψ k (·) is the Moreau proximal operator for ψ k (·). The soft thresholding operation is the proximal operator for L 1 -norm, whereas for L p -norm, the proximal operator is computed with an iterative algorithm called generalized soft thresholding (GST) [24], [39]. The proximal operator for nuclear norm is obtained via singular value soft thresholding as in [38], and for TV norm, it is efficiently computed by using Chambolle's method in [40].

III. SHIP WAKE DETECTION
Our detection algorithm includes three important steps as shown with different colored rectangles in Fig. 2:  1) preprocessing and inverse problem solution (blue rectangles); 2) detection of wakes in the Radon domain (gray rectangles); and 3) validation in the spatial domain (green rectangle). The blue rectangle within the red dashed-line shape in Fig. 2 constitutes the main contribution of this article, while the detection/validation steps are inspired by the method in [10] with all changes explained in detail in the sequel.
The preprocessing step consists of two steps including the creation of ship-centered and masked images and the inverse problem to obtainX MAP as discussed in the previous sections. For the masking operation by Graziano et al. [10], instead of replacing the area of the ship location with the mean intensity as is done in this article, only the unmasked pixels have been taken into account and the area containing the ship and land returns is ignored (for details see [10, p. 4]).
Upon obtaining the estimateX MAP , the first detection step is to restrict the peak/trough searching area in the Radon domain (r, θ). Since the image is centered on the ship, we ensure that the peaks/trough of all possible wakes are located between two sine waves [10], [14] with: |r | ≤ A sin θ where the sine wave peak point A refers to the maximum azimuth shift. Although the antenna velocity and the ship velocity along the slant range could be calculated using the slant range in the presence of field data, we chose to select A depending on the number of pixels in the image azimuth dimension M.
The value A has a crucial importance in the detection process. Selecting a large value for A increases the searching area, which is increasing the possibility of misdetecting noisy peaks as wakes. Conversely, selecting it as a small value may cause the ship wakes to fall outside the search area. Here, we use A = M/10, which we have shown to lead to the best results for a range of [M/15, M/5] in [41].
In the literature, it has been widely stated that the most visible wake is the turbulent wake [4], [15] and as stated in Section II-A, in most cases, it can be surrounded by bright narrow V-wake. Hence, the next step is detecting a peak/trough pair in the restricted area ofX MAP that corresponds to the turbulent and one arm of the narrow V-wake, respectively. A window of size 4 • for detecting turbulent/narrow V-wake pair is then scanned. The pair which maximizes the difference in amplitudes is selected as turbulent/narrow V-wake pair. Following the detection of the turbulent/narrow V-wake pair, the other arm of the narrow V-wake is searched on the other side of the turbulent wake within the same window size. The point with the maximum amplitude is selected as the second narrow V-wake.
As discussed in the previous sections, the Kelvin wake is outer signature of a moving vessel. In order to detect both Kelvin arms, both sides of the detected turbulent wake are searched with a half-angle window of size 10 • starting from ±10 • to ±20 • [2], [8] and points with maximum amplitudes on each side are then selected as Kelvin arms.
Up to this point, candidate wakes are detected by using the inverse problem solutionX MAP , and the next step consists in validating the detected wakes, which includes removal of half-lines and confirmation of detected wakes via measure indexes. The validation step is performed in the image domain since it has been stated in [21] that validation in the Radon domain might lead to erroneous results. Hence, points detected in the Radon domain are instead transformed into lines in the image domain via the inverse Radon transform.
The half of the lines is first removed. This has been referred to as the 180 • ambiguity problem in [10]. To solve this ambiguity, only the detected turbulent wake is used. The average intensity over the line representing the turbulent wake is calculated, and the half-line having lower average (since the turbulent wake is a dark line) is selected as the un-confirmed half-line corresponding to the turbulent wake. Half-lines belonging to the other detected wakes are selected if they are located in ±45 • on either side of the unconfirmed turbulent wake.
Confirmation of the candidate half-lines is then performed using a measure index F I , which is calculated as: F I =Ī w /Ī − 1 and interpreted as a measure of the difference between the average intensity over the unconfirmed wake,Ī w , and the average intensity of the image,Ī . The F I index is positive for bright wakes and negative only for the turbulent wake. Moreover, deciding a margin will help to reduce the possibility of false confirmations. Thus, we assumed a margin of 10% after a trial-error procedure (see Section V-A). Detected half-lines, which do not follow: F I < 0 for turbulent wake F I > 0.1 for narrow-V and Kelvin wakes (35) are discarded, whereas the remaining half-lines are confirmed.

IV. SAR DATA SETS
In this section, we describe the data sets employed in the experimental part of this article. These consist of both simulated and real SAR images.

A. Simulated SAR Images of The Sea Surface
In order to generate simulated SAR images of the sea surface, a numerical SAR image simulation software was developed based on previous studies in [8], [42], and [43]. The first part of the simulations consists in sea surface modeling, where the linear theory of surface waves was used. The ship wake modeling (Kelvin wake) considers ships as rigid bodies moving in inviscid incompressible fluid. Here, basic parameters belonging to ship such as length, beam, draft, and the Froude number are used to model different types of wakes.
The SAR imaging mechanism is usually described via RAR and specific SAR imaging. The RAR imaging is represented as an NRCS backscattering with VV and HH signal polarization modes and two linear modulations of tilt and hydrodynamic. Specific SAR imaging is based on a velocity bunching modulation (for details see [8], [44], and references therein).

B. Real SAR Data
SAR images from four different satellite platforms, namely TerraSAR-X, COSMO-SkyMed, Sentinel-1, and ALOS2, are employed. All TerraSAR-X images [45] are X-band Stripmap products with 3-m resolution for both azimuth and range directions. We have used three images, two of which are in HH polarization with the third in VV polarization. Of the COSMO-SkyMed data sets, we have used 2 X-band images including ship wake patterns. These are Stripmap products and their resolution is 3 m for both azimuth and range directions, in different polarizations. The Sentinel-1 images utilized in this article are C-band SAR images and have 10-m resolution for both azimuth and range directions. We have used four images, all of which are in VV polarization. Finally, we have used two ALOS2 images that are L-band and in HH polarization. ALOS2 images have 3-m resolution. From all images, we have selected 28 different ship wake patterns (28 different ships with corresponding wakes) by visual inspection and used them for experimental analysis.
We created ship centered image tiles and used them as input for the ship wake detection procedure discussed in the previous section. We assume that ship locations for the selected ships are known. In practice, this step can be replaced with ship detection techniques such as CFAR approaches (see [46]).
All ships in the created ship-centered image tiles are first masked to remove the bright spots in the Radon domain resulting from the ship itself. This allows us to discriminate bright points as possible ship wakes in the Radon domain. Following the preprocessing operations, the proposed inverse problem-based ship wake detection algorithm has been tested using the 28 aforementioned image tiles. In Table I, each image tile is shown with their visible wakes and corresponding details. As it can be seen, only two image tiles (1.1 and 2.6) out of 28 have all five possible types of ship wakes. The remaining images have less than five wakes, and the performance of the proposed method and the reference methods are evaluated as true detections of ship wakes. True detection in experimental analysis implies detecting visible wakes and discarding invisible wakes.

V. EXPERIMENTAL RESULTS
The proposed method was tested from three different perspectives using both simulated and real data. 1) We first used simulated SAR images of the sea surface including ship wakes. 2) Subsequently, we conducted experiments to determine the best choice of prior for solving the inverse problem in (1). 3) Finally, we compared the best method tested at 2) to state-of-the-art approaches for ship wake detection. The performance comparison was carried out in terms of receiver operation characteristics (ROC) as well as other measures, which are specifically described in Table II. Sensitivity quantifies the number of true positive (TP) in proportion to the total number of detections and specificity does the same for true negatives (TNs). The percentage accuracy shows the correct detection percentage over TP and TN values and was used to assess the ship wake detection performance of the method directly. The F 1 score is also a metric, which shows the accuracy of the methods. Positive likelihood ratio (LR+) shows how likely (or suitable) the utilized model is to detect correct classes. The last metric used in this article is Youden's J index, which shows the success of the test with values between 0 and 1. The value 0 implies that the  I   DETAILS FOR REAL SAR IMAGES   TABLE II DESCRIPTIONS OF ALL PARAMETERS, PRIORS, AND PERFORMANCE COMPARISON METRICS "test is unsuccessful," whereas 1 implies that "the test is successful."

A. Wake Detection in Simulated SAR Images
In the first set of experiments, we tested the proposed method with simulated SAR images of the sea surface. We generated two SAR images with VV polarization and a 4-m/s wind speed, 50 m size of ship, which is moving with a velocity of 9 m/s. In addition, both images are with different ship orientation, which effects the visibility of the wakes [47]. Inverse problem solution with all employed priors is obtained, and the procedure described above for detecting ship wakes is modified just to detect Kelvin arms instead of all five wakes as simulated images only contain one or two visible Kelvin arms.
A comparison study has been carried out for both images, and all visible wakes are successfully detected for all types of priors. Moreover, we apply different F index values to determine the most suitable value for F in the confirmation step, and we conclude that the F index value of 0.1 is enough to discard wrong detections, which follows our assumption for real SAR images in this article. In Fig. 3, simulated images and the wake detection results for GMC, L 1 , and TV are depicted. On examining the figure, it is clear that all three visible Kelvin wakes are detected, and the invisible one in Fig. 3(a) is discarded.

B. Inverse Problem Solution for Various Priors
In the second set of experiments, we tested the proposed ship wake detection framework in the context of using different types of priors. The results are presented in Table III. The performance metrics in Table III correspond to all 28 real SAR images in our data set. Examining the values in Table III, it is obvious to state that the inverse problem approach based on the GMC achieves the highest performance metrics among all methods tested. The percentage accuracy is around 10% higher than the TV, which is the second best.
ROC curves for all priors are depicted in Fig. 4, which demonstrates ship wake detection performance via TPR (or Sensitivity) versus FPR (or 1-Specificity). Examining the plots in Fig. 4, it can apparently be seen that the proposed   GMC-based inverse problem method outperforms all other priors in terms of ROC analysis.
In Fig. 5, we show visual results to assess the ship wake detection performance. In Fig. 5(a), the image for wake 2.2 is shown. There are three visually detected wakes corresponding to turbulent wake, one arm of narrow V-wake, and one Kelvin arm. Specifically, in Fig. 5, for wake image 2.2, GMC and TV detect two out of three visible wakes and discard all invisible wakes, which leads to an 80% detection accuracy. For the rest of the priors, detection accuracy is lower, e.g., for L 0.5 and L 1 , it is 60%.

C. Comparison to the State of the Art
In the third set of experiments, we used the best method from the previous section, i.e., the one based on the GMC and compared it to two state-of-the-art methods: 1) the ship wake detection method proposed by Graziano et al. [10] and 2) the log-regularized Hough transform-based method (Log-Hough) [23]. The results are presented in Tables IV and V. The performance of all three methods over all data sets is presented in Table IV in terms of the same metrics discussed in the previous section. The proposed method outperforms the reference methods by 10% in terms of accuracy and to various degree about the other performance metrics presented. Specifically, the TP value is higher than for the other methods, whereas the TN value is slightly lower than Graziano et al. [10]. (h) Nuclear. Yellow, green and red lines represent turbulent, narrow-V and Kelvin wakes, respectively. Ship wake image in (a) has only three visible wakes, which are turbulent, one narrow-V and one Kelvin wakes. Detecting these three wakes and discarding two unimaged wakes at the same time has 100% accuracy.  Table V for each data source, the proposed method achieves at least 75% detection performance for each satellite platform. However, for Sentinel-1, the method of Graziano et al. [10] achieves the best detection results. The lower performance of the proposed method for Sentinel-1 images can be explained in several ways. First, the resolution of Sentinel-1 images is less than that of the other platforms (i.e., 10 m). We believe that the low resolution affects the performance of the proposed method since the visibility of the wakes is reduced in low resolution. Besides, as high frequency (smaller wavelength), SAR sensors determine more surface scattering from rough surfaces, images in L-band (ALOS2) and C-band (Sentinel-2) include less sea surface details than X-band (TerraSAR-X and COSMO-SkyMed). This directly influenced the performance of the proposed method and determined the most accurate results to be obtained for X-band images.
When examining the TN values in Table V for Sentinel-1 images, the proposed method achieves 37.5%, whereas the method of Graziano et al. 50%, even though the TP value of the proposed method is the highest. Enhancing the image in Radon space using the proposed methodology increased the detectability of the visible wakes (TP values) in the images. However, in discarding the invisible wakes (TN values), Graziano's method is generally better than the pro-  [15] performance metrics that quantify false and lost wakes. In particular, FP and FN can be defined as the percentage of false and lost wakes, respectively. Examining FN and FP values, we can clearly see that the proposed method is substantially better than the reference methods in terms of FP results (false wakes), which is obvious especially in TerraSAR-X   Fig. 6, ROC curves for the methods tested in this section are presented. The superiority of the proposed method compared to existing ones is obvious. Graziano et al. [10] is the closest to the proposed method in terms of the ROC curve. The Log-Hough falls short compared to the proposed method and Graziano et al., which is not surprising since the log-regularization in [23] converges to an L p when p tends to 0. Example detection results for various images are presented in Fig. 7.

VI. CONCLUSION
In this article, we proposed a novel automatic approach for ship wake detection in SAR images of the sea surface acquired by various satellite platforms including TerraSAR-X, ALOS2, COSMO-SkyMed, and Sentinel-1. The proposed approach handles wake detection as an inverse problem to enhance information in the Radon domain to promote linear features. The solution to the corresponding optimization problem is obtained via a Bayesian formulation using a MAP estimator. The proposed method based on the GMC prior was first compared to various benchmark priors that are based on L 1 , L p , T V , and the nuclear norms. The GMC-based method was then compared to two state-of-the-art methods including Graziano et al. [10], [22] and Log-Hough [23]. The superiority of the GMC-based method has been clearly demonstrated in both sets of simulations with at least 10% accuracy gain over all data sets.
We conclude that enhancing sea SAR images in Radon space provides a more suitable platform for detecting the peaks/through compared to the direct approach by Graziano et al. Nevertheless, merely solving an inverse problem is not sufficient to obtain the most accurate results, since the choice of the prior has a crucial effect on the results. Indeed, only two priors (GMC and TV) lead to higher accuracy than Graziano et al.
Furthermore, the proposed approach differs from other approaches in the literature, which require the use of some ad hoc thresholding for enhancing the image. The main tunable parameter in our method is the regularization constant, λ. The need to adjust the scale parameter will be removed in future studies, which will consider a hierarchical Bayesian inference step. Investigating more complex sparsity enforcing priors in conjunction with nonconvex optimization algorithms is also one of our current endeavors.