Ship Wake Detection in SAR Images via Sparse Regularization

TL;DR

This paper introduces a novel sparse regularization method for detecting ship wakes in SAR images, outperforming existing techniques with an 80% success rate by formulating the problem as a convex inverse problem with a GMC penalty.

Contribution

It proposes a new inverse problem approach using GMC regularization for ship wake detection in SAR images, demonstrating superior performance over state-of-the-art methods.

Findings

GMC regularization achieves the best detection results.

The proposed method outperforms existing approaches.

80% success rate in diverse SAR images.

Abstract

In order to analyse 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 paper, 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 non-convex 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 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 ALOS2. The performance of various priors in solving the proposed inverse problem is first studied by investigating the GMC along with the L1, Lp, 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.

Figures

Tables

Table 1. TABLE I: Details for Real SAR Images.

						Wake		Visible Wakes
Image	Satellite	Frequency Band	Polarization	Resolution		Number		Turbulent	Narrow1	Narrow2	Kelvin1	Kelvin2
Image 1	TerraSAR-X	X (8-12 GHz)	HH	$3\times 3$ m		1.1		✓	✓	✓	✓	✓
						1.2		✓	✓	✗	✓	✗
						1.3		✓	✓	✗	✓	✗
						1.4		✓	✓	✗	✗	✗
Image 2	TerraSAR-X	X (8-12 GHz)	HH	$3\times 3$ m		2.1		✓	✓	✗	✓	✓
						2.2		✓	✓	✗	✗	✗
						2.3		✓	✓	✗	✓	✗
						2.4		✓	✓	✗	✓	✗
						2.5		✓	✓	✓	✓	✓
Image 3	TerraSAR-X	X (8-12 GHz)	VV	$3\times 3$ m		3.1		✓	✓	✗	✓	✗
						3.2		✓	✓	✗	✓	✗
						3.3		✓	✓	✗	✓	✗
Image 4	ALOS2	L (1-2 GHz)	HH	$3\times 3$ m		4.1		✓	✓	✗	✗	✗
Image 5	ALOS2	L (1-2 GHz)	HH	$3\times 3$ m		5.1		✓	✓	✗	✗	✗
						5.2		✓	✓	✗	✗	✗
						5.3		✓	✓	✗	✓	✗
Image 6	COSMO-SkyMed	X (8-12 GHz)	VV	$3\times 3$ m		6.1		✓	✓	✓	✓	✗
						6.2		✓	✓	✗	✗	✗
						6.3		✓	✓	✗	✓	✗
Image 7	COSMO-SkyMed	X (8-12 GHz)	HH	$3\times 3$ m		7.1		✓	✓	✗	✓	✗
Image 8	Sentinel-1	C (4-8 GHz)	VV	$10\times 10$ m		8.1		✓	✓	✗	✗	✗
						8.2		✓	✓	✗	✓	✗
						8.3		✓	✓	✗	✗	✗
Image 9	Sentinel-1	C (4-8 GHz)	VV	$10\times 10$ m		9.1		✓	✓	✗	✗	✗
Image 10	Sentinel-1	C (4-8 GHz)	VV	$10\times 10$ m		10.1		✓	✓	✗	✗	✗
						10.2		✓	✓	✗	✗	✗
						10.3		✓	✓	✗	✗	✗
Image 11	Sentinel-1	C (4-8 GHz)	VV	$10\times 10$ m		11.1		✓	✓	✗	✗	✗

Table 2. TABLE II: Descriptions OF all parameters, priors and performance comparison metrics.

Expression		Description
$M$		Image dimension. ($M\times M$)
$k$		Parameter defines GMC, $L_1$, $L_p$, TV and nuclear priors for
		1, 2, 3, 4 and 5, respectively.
${\lambda }_k$		Regularisation constant for prior $k$.
$\gamma$		Parameter for GMC which controls convexity. Set to 0.9.
${ϵ}^{\left(i\right)}$		Error criterion term for iteration $i$ in Algorithm 2. Refer to
		(22)
$MaxIter$		Maximum number of iterations for Algorithms 1 and 2.
		Set to 1000.
${\psi }_k\left(X\right)$		Negative logarithm of prior $k$. Please refer to (19), (24), (26),
		(29) and (31) for $k=1,2,3,4,5$, respectively.
$p$		$L_p$ norm order. .
$A$		Maximum sine wave amplitude for wake detection process.
		Set to $M/10$. Refer to Section III.
Kelvin wake searching range		An interval of $\pm [{10}^o,{20}^o]$ on either side of turbulent wake.
Narrow V-wake searching range		An interval of $\pm [0^o,4^o]$ on either side of turbulent wake.
True Positive (TP)		Correct confirmation of visible wakes.
True Negative (TN)		Correct discard of invisible wakes
False Positive (FP)		False detection of invisible wakes
		FP includes mislocated wake confirmations.
False Negative (FN)		False discard of visible ship wakes.
N		Total number of all possible wakes.
		Equal to 140 (28 Images $\times$ 5 possible wakes)
Sensitivity		TP/(TP+FN)
Specificity		TN/(TN + FP)
% Accuracy		100(TP+TN)/N
$F_1$		2TP/(2TP + FP + FN)
LR+		Sensitivity/(1-Specificity)
Youden’s $J$		Sensitivity + Specificity - 1

Table 3. TABLE III: Detection performance of different priors over all data sets.

				TP	TN	FP	FN	Sensitivity	Specificity	% Accuracy	$F_1$	LR+	Youden’s $J$
GMC				49.29%	30.71%	17.86%	2.14%	95.83%	63.24%	80.00%	0.83	2.61	0.59
$L_1$				35.00%	28.57%	32.14%	4.29%	89.09%	47.06%	63.57%	0.66	1.68	0.36
$L_{02}$				40.00%	26.43%	32.86%	0.71%	98.25%	44.58%	66.43%	0.70	1.77	0.43
$L_{05}$				37.86%	30.00%	30.71%	1.43%	96.36%	49.41%	67.86%	0.70	1.90	0.46
$L_{08}$				33.57%	32.14%	29.29%	5.00%	87.04%	52.33%	65.71%	0.66	1.83	0.39
TV				35.71%	34.29%	26.43%	3.57%	90.91%	56.47%	70.00%	0.70	2.09	0.47
Nuclear				33.57%	29.29%	32.14%	5.00%	87.04%	47.67%	62.86%	0.64	1.66	0.35

Table 4. TABLE IV: Detection performance comparison of methods over all data sets.

				TP	TN	FP	FN	Sensitivity	Specificity	% Accuracy	$F_1$	LR+	Youden’s $J$
Proposed (GMC)				49.29%	30.71%	17.86%	2.14%	95.83%	63.24%	80.00%	0.83	2.61	0.59
Graziano et. al. [10]				33.57%	35.71%	26.43%	4.29%	88.68%	57.47%	69.29%	0.69	2.09	0.46
Log-Hough [23]				38.57%	25.71%	33.57%	2.14%	94.74%	43.37%	64.29%	0.68	1.67	0.38

Table 5. TABLE V: Detection performance comparison of methods for each data sets separately.

				TerraSAR-X
				TP	TN	FP	FN	Sensitivity	Specificity	% Accuracy	$F_1$	LR+	Youden’s $J$
Proposed (GMC)				61.67%	21.67%	13.33%	3.33%	94.87%	61.90%	83.33%	0.88	2.49	0.57
Graziano et. al. [10]				35.00%	23.33%	33.33%	8.33%	80.77%	41.18%	58.33%	0.63	1.37	0.22
Log-Hough [23]				48.33%	15.00%	33.33%	3.33%	93.55%	31.03%	63.33%	0.73	1.36	0.25
				ALOS2
				TP	TN	FP	FN	Sensitivity	Specificity	% Accuracy	$F_1$	LR+	Youden’s $J$
Proposed (GMC)				40.00%	35.00%	25.00%	0.00%	100.00%	58.33%	75.00%	0.76	2.40	0.58
Graziano et. al. [10]				30.00%	40.00%	30.00%	0.00%	100.00%	57.14%	70.00%	0.67	2.33	0.57
Log-Hough [23]				45.00%	35.00%	20.00%	0.00%	100.00%	63.64%	80.00%	0.82	2.75	0.64
				COSMO-SkyMed
				TP	TN	FP	FN	Sensitivity	Specificity	% Accuracy	$F_1$	LR+	Youden’s $J$
Proposed (GMC)				45.00%	40.00%	10.00%	5.00%	90.00%	80.00%	85.00%	0.86	4.50	0.70
Graziano et. al. [10]				40.00%	40.00%	15.00%	5.00%	88.89%	72.73%	80.00%	0.80	3.26	0.62
Log-Hough [23]				35.00%	30.00%	30.00%	5.00%	87.50%	50.00%	65.00%	0.67	1.75	0.38
				Sentinel-1
				TP	TN	FP	FN	Sensitivity	Specificity	% Accuracy	$F_1$	LR+	Youden’s $J$
Proposed (GMC)				37.50%	37.50%	25.00%	0.00%	100.00%	60.00%	75.00%	0.75	2.50	0.60
Graziano et. al. [10]				30.00%	50.00%	20.00%	0.00%	100.00%	71.43%	80.00%	0.75	3.50	0.71
Log-Hough [23]				22.50%	35.00%	42.50%	0.00%	100.00%	45.16%	57.50%	0.51	1.82	0.45

Equations

$$Y=\mathcal{C}X+W,$$

$$X(r,\theta )={\int }_{\text{mathds}{R}^2}Y(i,j)\delta (r-i\cos \theta -j\sin \theta )didj$$

$$g(r,\theta )$$

$$Y(i,j)$$

$$p(x|y)=\frac{p(y|x)p(x)}{\int p(y|x)p(x)dx}$$

$$p(x|y)\propto p(y|x)p(x).$$

$$p(x|y)\propto \exp \{-F(x)\}$$

$${\hat{x}}_{\text{MAP}}=\arg \underset{x}{\max }p(x|y)=\arg \underset{x}{\min }F(x).$$

$$p(X|Y)\propto p(Y|X)p(X).$$

$$p(Y|X)\propto \exp \{-\Vert Y-\mathcal{C}X{\Vert }_2^2\}.$$

$$F(X)\propto \Vert Y-\mathcal{C}X{\Vert }_2^2+\lambda \psi (X),$$

$$\psi (t)=\{\begin{array}{ll}|t|-\frac{1}{2}{t}^2& |t|\le 1\\ \frac{1}{2}& |t|\ge 1.\end{array}$$

$$\psi (t)=|t|-s(t).$$

$${\psi }_b(t)=|t|-s_b(t).$$

$$S_B(t)=\underset{v}{\inf }\{\Vert v{\Vert }_1+\frac{1}{2}\Vert B(t-v){\Vert }_2^2\}.$$

$${\psi }_B(t)=\Vert t{\Vert }_1-S_B(t).$$

$$B=\sqrt{\frac{\gamma }{{\lambda }_1}}\mathcal{C}$$

$$p_1(X)$$

$$\propto \exp \{-{\lambda }_1\left(\Vert X{\Vert }_1-S_B(X)\right)\}.$$

$$F(X,v)=\Vert Y-\mathcal{C}X{\Vert }_2^2$$

$${\hat{X}}_{\text{MAP-GMC}}=\arg \underset{X}{\min }\underset{v}{\max }F(X,v).$$

$${ϵ}^{(i)}=\frac{\Vert X^{(i)}-X^{(i-1)}\Vert }{\Vert X^{(i-1)}\Vert }.$$

$$p_2(X)$$

$$\propto \exp \{-{\lambda }_2\Vert X{\Vert }_1\}.$$

$$p_3(X)$$

$$\propto \exp \{-{\lambda }_3\Vert X{\Vert }_p^p\}$$

$$TV(X)=\Vert \nabla X{\Vert }_1.$$

$$p_4(X)$$

$$\propto \exp \{-{\lambda }_{4}TV(X)\}.$$

$$p_5(X)$$

$$\propto \exp \{-{\lambda }_5\Vert X{\Vert }_{\ast }\}.$$

$${\hat{X}}_{\text{MAP}-{\psi }_k}=\arg \underset{X}{\min }\{\Vert Y-\mathcal{C}X{\Vert }_2^2$$

$${\Gamma }_{{\lambda }_k}(u)=\frac{|pro{x}_{{\psi }_k}^{{\lambda }_k}(u)|}{|pro{x}_{{\psi }_k}^{{\lambda }_k}(u)|+{\lambda }_k}u\text{for }k=2\text{ or }3,$$

Full text

I Introduction

Accurate characterisation of sea surface condition is not only important in isolation, but also in the detection and characterisation 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 modelling was the lack of data of sufficiently high resolution (pixels need to be typically smaller than few meters) and accuracy. Synthetic aperture radar (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 synthetic aperture radar (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. Whereas specular reflections are usually caused by breaking waves, the Bragg scatterers are short gravity waves or capillary waves depending on the wavelength of the SAR sensor itself. 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 which are directly determined by different ship wake formations. These are typically considered to fall in one of three categories: (i) turbulent wakes, (ii) surface waves created by ships and (iii) ship-generated internal waves. Ship-generated surface waves can in turn be split into two sub-categories [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 are 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. performed a validation study for simulated ship wakes via computational dynamics in [9].

Since ship wakes can be modelled 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 utilised by Murphy [11]. The Radon transform has however a couple of drawbacks, e.g. bright pixels belonging to ships may cause false detections [12]. To address this, enhancing the Radon domain information is common practice in the literature. Rey et al. [13] have proposed a method which combines the Radon transform and Wiener filtering to increase the detectability of the peaks in Radon domain. Tunaley [14] used a method which restricts the search area in Radon domain. Eldhuset [15] proposed an automatic ship and wake detection method whereby the detection performance is characterised by the number of lost and false wakes (wake-like 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}$ regularised logistic regression [19] and low-rank plus sparse decomposition [20] also addressed the ship wake detection problem in SAR images. Graziano et. al. [21, 12, 10, 22] have proposed wake detection methods which deal directly with the noisy image without performing any preliminary enhancement. The authors have first created ship-centred-masked 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 regularisation 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 the literature [28, 29, 30, 31, 32].

In this paper, 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 paper is to propose an innovative approach based on sparse regularisation 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 maximum a-posteriori (MAP) estimator. Our proposed cost function uses the generalised minimax concave (GMC) penalty of Selesnick [33] as regularisation term and investigates its merit in comparison to the total variation (TV), nuclear, $L_{1}$ and $L_{p}$ norms. The use of the GMC penalty demonstrates the advantages of non-convex sparse regularisation 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 which performs detection in the Radon domain as proposed in [10, 22]. In this study, contrary to [10, 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 *forward-backward *(FB) method and the two-step iterative shrinkage/thresholding (TwIST) algorithm [34].

The rest of the paper is organised as follows: the image formation model for ship wakes identification, the MAP estimation, GMC regularisation 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 the paper with a brief summary.

II Theoretical Preliminaries

II-A Ship Wakes and Image Formation Model

As discussed briefly in the previous section, 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 V-wake, which lies either side of the turbulent wake within the half angle from $1.5^{\circ}$ up to $4^{\circ}$ [16, 10]. Lastly, 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^{\circ}$. The Kelvin wake half angle can be somewhat smaller in real SAR images, e.g. between $10^{\circ}$ and $19.5^{\circ}$ [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 (1), the first (2) and second (3) narrow V - and the first (4) and second (5) 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 Figure 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 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

[TABLE]

where $Y$ is the $M\times M$ SAR image, $X(r,\theta)$ is the image in Radon domain, $W$ refers to the additive noise and the operator $\mathcal{C}=\mathcal{R}^{-1}$ is the inverse Radon transform. $X(r,\theta)$ represents lines as a distance $r$ from the centre of $Y$ and an orientation $\theta$ 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 $\theta$ leads to the Radon transform $X(r,\theta)$ of the given image $Y$. It can also be defined as a projection of the image along the angles, $\theta$. Hence, for a given image $Y$, the general form of the Radon transform ($X=\mathcal{R}Y=\mathcal{C}^{T}Y$) is

[TABLE]

where $\delta(\cdot)$ is the Dirac-delta function.

The inverse Radon transform ($Y=\mathcal{C}X$) of the projected image $X$ can be obtained from the filtered back-projection [35] algorithm as

[TABLE]

where $v$ is the radius in Fourier transform, and $\mathcal{F}\left[\cdot\right]$ and $\mathcal{F}^{-1}\left[\cdot\right]$ refer to forward and inverse Fourier transforms, respectively. In this paper, we use discrete operators $\mathcal{R}$ and $\mathcal{C}$ as described in [36].

II-B Bayesian Inference in Inverse Problems

Assume 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:

[TABLE]

where, the denominator $\int 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

[TABLE]

Generally, posterior distributions are assumed to be log-concave as

[TABLE]

where $F(x)$ is a convex function. Estimating $x$ directly using the posterior density can be intractable especially for high dimensional cases. Hence, the maximum a-posteriori (MAP) estimator maximises the posterior $p(x|y)$ to obtain a point estimate $\hat{x}$. The MAP estimator for the unknown $x$ is

[TABLE]

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

[TABLE]

Assuming an i.i.d. standard normal noise case, the likelihood distribution $p(Y|X)$ can be expressed as

[TABLE]

In this paper, we assume priors to be of exponential form as, $p(X)\propto\exp\{-\lambda\psi(X)\}$. Thus, the cost function in (8) becomes

[TABLE]

where $\lambda$ is the scale parameter, namely the regularisation constant.

II-C Generalised Minimax Concave (GMC) Penalty

The GMC regularisation has been proposed by Selesnick [33] as a sparsity enforcing penalty in inverse problems. It exploits the advantages of using non-convex penalties as well as preserving the convexity of the cost function. It is based on $L_{1}$ norm and the generalised Huber function.

In particular, the minimax concave (MC) penalty in the univariate case is

[TABLE]

The relation between the MC penalty and the Huber function $s(t)$ can be written as

[TABLE]

For a scalar, $b\neq 0$, the scaled MC penalty, $\psi_{b}(t)$ can be written as scaled Huber function, $s_{b}(t)$ as

[TABLE]

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 generalised Huber function, $S_{B}(t)$ can be defined as [33]

[TABLE]

Combining (14) and (15) gives the generalised MC (GMC) penalty function as

[TABLE]

The scaling matrix $B$ should be selected in relation to the inverse Radon operator $\mathcal{C}$ to provide the convexity of the cost function as

[TABLE]

where $\lambda_{1}$ is the scale parameter of the GMC prior and $\gamma$ is a parameter which controls the non-convexity. Note that for $0\leq\gamma\leq 1$, $B$ ensures the convexity of the cost function. The nominal range of $0.5\leq\gamma\leq 0.9$ should be used for better performance [33].

Thus, the GMC sparse prior can be written as

[TABLE]

As GMC regularisation does not have an explicit formulation, the minimisation problem with the cost function in (11) can be solved using proximal operators. Thus, we rewrite the cost function as

[TABLE]

which leads to a minimax optimisation problem

[TABLE]

The solution to this problem can be obtained using the forward-backward (FB) algorithm [33]. The FB algorithm for GMC regularised cost functions requires only a couple of computational steps and soft-thresholding, which is the proximal operator for $L_{1}$ based regularisation. 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 $\epsilon^{(i)}$ is calculated for iteration $i$ as

[TABLE]

II-D Alternative Sparse Priors

In this paper, we also investigate other types of priors, which are known to be sparse and common in optimisation problems. The first of which is Laplace (or $L_{1}$ norm) prior, which is log-concave and given by

[TABLE]

We further investigate the (non-convex) $L_{p}$ norm based prior, which is

[TABLE]

where $0<p<1$.

Another important type of prior is the TV norm, $TV(\cdot)$, which can be defined as [38]

[TABLE]

where $\nabla$ is the 2-D discrete gradient operator. Hence, the TV norm based prior is given by

[TABLE]

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:

[TABLE]

Replacing $\psi(X)$ in (11) with the previously defined priors determines different MAP estimators as

[TABLE]

where $k=2,3,4,5.$

Minimisations of (32) are subsequently carried out with the two-step iterative shrinkage/thresholding (TwIST) algorithm [34]. The TwIST algorithm is given in Algorithm 2 where $MaxIter$ and the stopping criterion $\epsilon$ are the same as in the GMC case.

The operator $\Gamma_{\lambda_{k}}(\cdot)$ in Algorithm 2 is a shrinkage/thresholding/denoising function. In this paper, we used shrinkage/thresholding operators in minimizers with $L_{1}$ and $L_{p}$ for $\Gamma_{\lambda_{k}}(\cdot)$ as also used in [24]

[TABLE]

and the denoising operators for the TV and nuclear norm priors,

[TABLE]

where $prox_{\psi_{k}}^{\lambda_{k}}(\cdot)$ is the Moreau proximal operator for $\psi_{k}(\cdot)$.

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 generalised soft thresholding (GST) [39, 24]. 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 coloured rectangles in Figure 2: 1) Pre-processing and inverse problem solution (blue rectangles), 2) detection of wakes in the Radon domain (grey rectangles) and 3) validation in the spatial domain (green rectangle).The blue rectangle within the red dashed-lined shape in Figure 2 constitutes the main contribution of this paper, while the detection/validation steps are inspired by the method in [10] with all changes explained in detail in the sequel.

The pre-processing step consists of two steps including the creation of ship-centred and masked images and the inverse problem to obtain $\hat{X}_{\text{MAP}}$ as discussed in the previous sections. For the masking operation in Graziano et. al. [10], instead of replacing the area of the ship location with the mean intensity as is done in this paper, only the unmasked pixels have been taken into account and the area containing the ship and land returns is ignored (for details see p.4 in [10]).

Upon obtaining the estimate $\hat{X}_{\text{MAP}}$, the first detection step is to restrict the peak/trough searching area in the Radon domain ($r,\theta$). Since the image is centered on the ship, we ensure that the peaks/trough of all possible wakes are located between two sine waves [14, 10] with: $|r|\leq A\sin\theta$ where the sine wave peak point $A$ refers to the maximum azimuth shift. Although antenna velocity and 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 mis-detecting 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 [15, 4] 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 of $\hat{X}_{\text{MAP}}$ that corresponds to the turbulent and one arm of the narrow V-wake, respectively. A window of size $4^{\circ}$ for detecting turbulent/narrow V-wake pair is then scanned. The pair which maximises 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, 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^{\circ}$ starting from $\pm 10^{\circ}$ to $\pm 20^{\circ}$ [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 solution $\hat{X}_{\text{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. Validation step is performed in 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 are first removed. This has been referred to as the $180^{\circ}$ 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 $\pm 45^{\circ}$ on either side of the un-confirmed turbulent wake.

Confirmation of the candidate half lines is then performed using a measure index $F_{I}$, which is calculated as: $F_{I}=\bar{I}_{w}/\bar{I}-1$ and interpreted as a measure of the difference between the average intensity over the un-confirmed wake, $\bar{I}_{w}$, and the average intensity of the image, $\bar{I}$. 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 (Please see Section V-A). Detected half lines which do not follow:

[TABLE]

are discarded, whereas the remaining half lines are confirmed.

IV SAR Data Sets

In this section, we describe the datasets employed in the experimental part of the paper. These consist of both simulated and real SAR images.

IV-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, 43]. The first part of the simulations consists in sea surface modelling, where the linear theory of surface waves was used. The ship wake modelling (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 real aperture radar (RAR) and specific SAR imaging. The RAR imaging is represented as a normalised radar cross section (NRCS) backscattering with vertical-vertical (VV) and horizontal-horizontal (HH) signal polarisation 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).

IV-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 three meters resolution for both azimuth and range directions. We have used 3 images, two of which are in HH polarisation with the third in VV polarisation. 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 three meters for both azimuth and range directions, in different polarisations. The Sentinel-1 images utilised in this paper are C-band SAR images and have ten meters resolution for both azimuth and range directions. We have used 4 images, all of which are in VV polarisation. Lastly, we have used 2 ALOS2 images which are L-band and in HH polarisation. ALOS2 images have three meters 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 constant false alarm rate (CFAR) approaches (see e.g. [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 pre-processing 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 2 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. i) We first used simulated SAR images of the sea surface including ship wakes. ii) Subsequently, we conducted experiments to determine the best choice of prior for solving the inverse problem in eq. (1). iii) Lastly, we compared the best method tested at (ii) 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 in proportion to the total number of detections and specificity does the same for true negatives. 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 utilised model is to detect correct classes. The last metric used in this paper is Youden’s $J$ index, which shows the success of the test with values between 0 and 1. The value 0 implies the “test is unsuccessful” whereas 1 implies “the test is successful”.

V-A Wake detection in simulated SAR images of sea surface

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 polarisation and a 4 m/s wind speed, 50 m size of ship, which is moving with 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 are 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 1 or 2 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 paper. In Figure 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 Figure 3-(a) is discarded.

V-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 corresponds 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 Figure 4 which demonstrates ship wake detection performance via true positive rate (TPR, or Sensitivity) vs. false positive rate (FPR, or 1 - Specificity). Examining the plots in Figure 4, it can apparently be seen that the proposed GMC based inverse problem method outperforms all other priors in terms of ROC analysis.

In Figure 5, we show visual results to assess the ship wake detection performance. In Figure 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 Figure 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%.

V-C Comparison to the state-of-the-art

In the third set of experiments, we used the best method from 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], ii) the log-regularised Hough transform based method (Log-Hough) [23]. The results are presented in Table IV and V.

The performance of all three methods over all data sets are presented in Table IV in terms of the same metrics discussed in 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].

When examining the values in Table V for each data source, the proposed method achieves at least 75% detection performance for each satellite platforms. However, for Sentinel-1, the method of Graziano et. al 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 are less than that of the other platforms (i.e. 10 meters). We believe that the low resolution affects the performance of the proposed method since 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- (ALOS2) and C- (Sentinel-2) bands 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 proposed method whereas it falls short in detecting the visible wakes. The values of FP and FN are also related to the Eldhuset’s [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 images where it results in 20% higher wake detection accuracy. For all the methods, FN values (the number of lost wakes) are relatively small and do not directly correlate with performance. Example detection results for various images are presented in Figure 7.

In Figure 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 ROC curve. The Log-Hough falls short compared to the proposed method and Graziano et. al, which is not surprising since the log-regularisation in [23] converges to an $L_{p}$ when $p$ tends to 0.

VI Conclusions

In this paper, 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 Radon domain to promote linear features. The solution to the corresponding optimisation 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 which are based on $L_{1}$, $L_{p}$, $TV$, 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 in 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 regularisation constant, $\lambda$. 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 non-convex optimisation algorithms is also one of our current endeavours.

List of Abbreviations

[TABLE]

Acknowledgment

We are grateful to the UK Satellite Applications Catapult for providing us the COSMO-SkyMed data sets employed in this study.