1. Introduction

Since a couple of years, the number of illegal acts for taking control of maritime vessels has increased. For search-and-rescue reasons, it is suitable to find efficient sensor systems for detecting vessels. Vessel candidates for illegal acts are often commercial vessels with great dimensions. Such vessels (and all vessels with length greater than 45 m) are required to set up an Automatic Identification System (AIS) used by vessel traffic services for identifying and locating vessels. However, when a ship is shutting down its AIS due to illegal acts or material defects, there are no means to detect it in open sea.

Spaceborne sensors are a valuable tool for noncooperative ship detection when the AIS is not operational. Two classes of spaceborne sensors exist: radar and electro-optical (EO) sensors. As compared to radar sensors, EO sensors have lower cost, lower payload, and better computational processing load. To have a high revisiting time, a constellation of LEO micro-satellites is used. Micro-satellites limit the sensor payload to a few kilograms. Currently, EO sensors are then the best candidate for spaceborne applications. To have day-night capabilities, infrared (IR) sensors are used.

Optimum design of such sensors implies to be capable of simulating the evolution of sensor performance as a function of sensor or scene parameters before manufacturing the sensor. Sensor performance is often expressed using Receiver Operating Characteristic (ROC) curves representing the evolution of the probability of detection with respect to the probability of false alarms. So far, these curves are computed using results of detection algorithms applied to the image received by the sensor. This implies the simulation of these images and the choice of detection algorithms. For our application, since the payload is very limited, the Ground Sampling Distance (GSD) is large (about 100 m). Hence, ship detection cannot solely be pixel-based. This indeed leads to an important rate of false alarms. One possible solution is to detect wakes behind the ship. At large GSD, the turbulent wake is the most visible. It appears bright in optical images (Figure 1) and dark in long-wave IR (LWIR) images (Figure 2). Computing the evolution of ROC curves with sensor or scene parameters is then computationally intensive.

Figure 1: Example panchromatic (optical) SPOT image of moving ships. Source: [2]. Turbulent wake appears bright.

Figure 2: Example thermal infrared (LWIR) LANDSAT image of moving ships. Turbulent wake appears dark. Colors correspond to normalized radiance received by the sensor.

This paper proposes a methodology having real-time capabilities for helping sensor manufacturers in optimizing the design of new EO sensors for maritime (ship detection) applications. This implies to be able to test, in real-time, the effect of sensor or scene parameters on ROC curves. In IR, near real-time simulators exist [1]. However, they were designed for airborne applications, for which ship detection is done using the contrast between ship and sea background pixels. This cannot be used for spaceborne EO sensors with large GSD. Hence, to our knowledge, there are no real-time tools available for simulating performance of spaceborne EO sensors with large GSD in a maritime environment.

Our approach is based on the one described in [3, 4], where real-time capabilities are obtained by computing ROC curves from a model of the probability density function (pdf) of pixels contained in the image. This avoids simulating the image received by the sensor. In [3, 4], this idea was developed for land-cover scene modeling using hyperspectral sensors. In a maritime environment, a very first attempt to model sea pixels with a pdf was described in [5] for LWIR airborne sensors. To our knowledge, such methodology has not yet been considered for ship detection. This is the subject of the present paper. Our real-time statistical methodology is described in the case of a mid-wave IR (MWIR) sensor. The result is a simple simulator that produces ROC curves in real-time. The proposed statistical methodology can be applied to other EO sensors and even to radar sensors, if appropriate models for the pdfs are used. Such a tool can be very useful for ship detection using Synthetic Aperture Radar (SAR), for which the simulation of SAR images is time consuming [6, 7].

We emphasize that our aim is not to provide a very accurate, validated simulator. Hence, in this paper, performance of the proposed tool is not deeply examined and this tool is not validated using real data. This will be the subject of further research. Our aim is only to propose a real-time methodology for assessing EO sensor performance and to illustrate this methodology by the design of a simple simulator for ship detection using MWIR sensors. Remember that this methodology is inspired from [3, 4].

Section 2 describes the wakes generated behind a moving ship. Section 3 defines ROC curves. Section 4 presents models used for the sea surface and for the turbulent wake. Sections 5 and 6 explain the model of the signal received at the sensor. Section 7 presents the real-time statistical simulator. Section 8 studies its performance. Section 9 concludes.

2. Wakes behind Moving Ships

If a ship is moving, wakes are generated behind it. These wakes are observed for any vessel speed and dimensions and can persist for hours and grow several tens of kilometers long, making it a feature which can easily be detected using spaceborne sensors. It can also provide information on the vessel's heading, speed, and potentially its hull dimensions, which makes it a very desirable feature for detection and tracking purposes. Therefore, wake detection is often used either in combination with or even instead of other ship detection methods.

A ship produces two types of wakes [8]. The turbulent wake, a zone of reduced sea surface roughness which appears as a long bright (optical sensors) or dark (LWIR sensors) streak behind the ship, bounded by a v-wake, and the Kelvin wake, a system of ripples occurring inside a cone of 39 degrees originating at the ship's bow. The Kelvin wake's wave spectrum can be analyzed for determining the ship's speed and heading, and its dimensions. Figure 3 shows a typical wakes pattern.

Figure 3: Different types of wakes appearing behind a moving vessel.

2.1. Kelvin Wake

The Kelvin wake consists of two systems of ripples, the transverse and divergent waves. These systems [9] are bounded by two cusp-lines separated by an angle of deg. On the cusp-line, a wave propagates with a wavelength depending upon the ship speed : with being the gravity constant.

2.2. Turbulent Wake

The turbulent wake is a zone of high-frequency low-amplitude waves behind the ship's stern. It behaves like a flat but rough surface, therefore contrasting with its surroundings. Hence, the physical quantities of interest are the width and the length of the wake. The turbulent wake's width depends upon ship dimensions, more specifically its beam (width) , and its length . We have where is the distance from the ship stern. Here and are derived from an empiric approximative formula for the turbulent wake width at four ship lengths. Experimental data show that is a good approximation, though can vary between and . In general, is a good approximation which varies very little for common ship designs [9], resulting in further simplification: .

The wake length is a more difficult problem and depends upon sea state. The turbulent wake is caused by water displacement due to the ship's hull and propulsion system. This water displacement has a kinetic energy decreasing according to [10]. As long as this kinetic energy is significantly larger than the energy of the top water layers, the turbulent wake remains detectable. Typically, turbulent wakes exist during a long period of time. Their length is typically a few kilometers. Figure 4 shows example simulated turbulent wake widths.

Figure 4: Width of the turbulent wake as a function of the distance behind the ship for various ship dimensions.

3. Definition of ROC Curves

ROC curves are an important signal processing tool for assessing the performance of a sensor or an algorithm. They rely on the definition of a probability density function (pdf) for the signal and the noise [11].

3.1. Signal and Noise Pdfs

In our application, the target signal is the turbulent wake radiance, and the noise signal is the open-sea radiance. Each signal is characterized by a pdf. We thus have two pdfs representing the statistical distribution of the wake signal and of the open-sea signal. They are, respectively, denoted and , where is the level of the signal displayed by the sensor.

3.2. ROC Curves and Detection Algorithm

The detection algorithm works as follows. The value of each pixel in the image received at the sensor is a realization of either or of (or a mix of both pdfs). The mean of each pdf is denoted and , respectively. To perform the detection, we apply a threshold to the pixels in the signal image. If , all pixels greater than are classified as target pixels and other pixels as noise pixels. However, among target pixels, some of them are noise pixels and thus correspond to false alarms. Below, we describe how to evaluate the rate of false alarms.

For a given , we can define a probability of detection and a probability of false alarms . If , and are given by

Hence, represents the probability that an open-sea pixel is classified as a wake pixel and represents the probability that a wake pixel is effectively classified as a wake pixel. and are represented graphically in Figure 5. Hence, for each , we have one and one . ROC curves are obtained by plotting versus for all possible values of . We can repeat the reasoning if . These curves serve as basis for discussing sensor performance: for a given , should be as high as possible. Below, we describe a model for and .

Figure 5: Probability of detection and probability of false alarms.

4. Sea and Turbulent Wake Surface Models

Finding and implies to compute the signal received at each pixel in the detector plane of the spaceborne sensor. There are mainly two classes of pixels: open-sea and wake pixels, respectively, containing open-sea and turbulent wake radiances. We first describe how the geometrical models of the sea surface and of the turbulent wake surface are obtained.

4.1. Open-Sea Surface Modeling

Our model is based on the model presented in [12, 13]. In realistic sea surface models, we consider three classes of waves: () capillarity waves with small wavelength ( cm) influenced by viscosity and surface tension, () gravity waves that are wind-driven waves with wavelength  cm and smaller than a few meters, () swells being waves with great wavelength, that is, is greater than a few meters (these waves originate due to the presence of wind. However, they remain active for a long time after the wind has blown), () choppy waves appearing for high wind speed and introducing nonlinearities in the sea surface model (they are the starting point of breaking waves and of the apparition of foam). We only consider gravity waves and swells.

To obtain the sea surface model, we divide the sea surface in small facets. Then, vertical displacements are applied to these facets. These displacements are obtained by modeling the sea surface as a superposition of linear plane waves [14]. A plane wave is given as where is the wave amplitude, is time, is the position vector, is the phase, and is the wave vector given by , where is the wave number where is the wavelength. is the direction of propagation of the plane wave. Gravity waves are modeled as a superposition of a great number of plane waves. Each wave is characterized by a value for , , and . The wave height at location and time is found by integrating the plane waves over the entire space spanned by . We thus have where Hence, is the inverse Fourier Transform (FT) of . Modeling gravity waves is done by specifying a model for and for . Modeling swells is done in the same way. The only difference is the model for and . For gravity waves, in the case where capillarity waves can be neglected, we have the dispersion relationship [12] where is the gravity constant. is modeled as a random process (RP) that determines the random character of wind-generated waves. Here, is modeled as a Gaussian RP with zero mean and unit variance. The model of depends upon wind speed and wind direction . We can write as where and is the power spectrum often given by the Pierson-Moskowitz spectrum [14], that is, where , , , and , where is the wind speed at 19.5 m above the sea level. There exist other spectra that are tailored to a particular sea [14].

In practice, in (4) is computed using the 2D inverse FFT (IFFT). Indeed, by discretizing as , where and and as , where and , (4) becomes

The length of the patch where the IFFT is computed is given by . The periodicity of the IFFT can be used to replicate the 's in both spatial directions. Hence, we can compute sea heights for extended surfaces at an acceptable computation cost. Figure 6 shows examples of sea surface heights generated with the previous model.

Figure 6: Examples simulated color-coded sea heights for a wind speed of 11 m/s: (a) gravity waves and (b) gravity waves and swells. -axis and -axis are labeled in meters. Color indicates sea heights (in meters). Sea height zero is the mean sea level.

Only considering gravity waves and swells for modeling sea surface is valid for low sea states. For high sea states (typically 5), breaking waves appear due to gravity. These waves are not handled in this model. The presence of breaking waves only modifies the model for ; the principles of the method remain unchanged.

4.2. Turbulent Wake Surface Model

A turbulent wake is modeled as a very rough flat surface [9, 10]. Hence, we model this wake as a flat sea. This flat sea is divided into microfacets (to simulate turbulences), the orientation of these micro-facets being uniformly distributed between and to simulate surface roughness.

In Section 5, we see that sea water emissivity (resp., reflectivity) goes down (resp., up) as the angle of arrival of the optical beam on a sea facet increases. Hence, wakes can be distinguished from open-sea thanks to a change in the emissivity (or reflectivity) between wake and open-sea pixels. For optical sensors, the wake appears bright (Figure 1) due to a higher value (higher reflectivity) of the sun glint for wake than for open-sea pixels. For LWIR sensors, the wake appears dark due to a reduction in the emissivity of the sea surface in the wake compared to its value for open-sea pixels. For MWIR sensors, there is a competition between reflection (sun glint and sky irradiance) on sea facets and self-emission of sea facets. This is discussed further below.

5. Radiance Received at the Sensor

Below, we present a model for computing the radiance received at the entry of the sensor. This model can be applied to open-sea and wake pixels.

5.1. Radiance at the Sea Surface (One Sea Facet)

We first describe the method for computing the radiance leaving one sea facet . The radiance leaving for wavelength is computed using the following equation [12, 13, 15]:

We describe below a real-time model for each term in (9). Figure 7 defines useful variables relative to .

Figure 7: Useful variables for a sea facet .

5.1.1. Emitted Radiance

In (9), represents the radiance emitted by due to its nonzero temperature. It is computed using Planck's law [16], that is, where is the open sea water emissivity at , is the absolute open sea surface temperature, and is the blackbody radiance [16]. if , zero otherwise. The variation of with is mainly due to the variation of with the elevation angle of the optical beam that goes to the sensor. Neglecting the dependence upon wavelength, we have [17]

For wake pixels, is the mean of for . Hence, . Hence, for wake facets, for most values of .

5.1.2. Sky Radiance

In (9), is the irradiance produced by the sky. It is present at any time. There are two models. The first model as a blackbody at sky temperature (depending upon weather parameters) [18, 19], that is, where is the sea reflectivity and is a visibility factor representing the portion of the sky hemisphere seen by . This model is realistic under clear-sky conditions. The second model uses MODTRAN [20]. However, it is computationally intensive. To have a real-time model, we use the blackbody model.

5.1.3. Solar Irradiance

The solar extraterrestrial radiation not back-scattered to space reaches the ground in two ways. The radiation reaching the ground directly is the beam irradiance. The scattered radiation reaching the ground is the diffuse irradiance. Below, we assume clear-sky conditions. (For images containing clouds, we assume that cloud masking algorithms [21, 22] have been applied prior to ship detection.) The beam irradiance incident on a surface of 1  on the earth's ground is where where is a proportionality constant and is obtained from MODTRAN and accounts for propagation through the atmosphere. is the spectral radiance of the sun computed either using a blackbody at  K or using MODTRAN (more accurate). In (9), corresponds to the reflected beam solar irradiance, that is, the solar glint. Assuming that is a Lambertian (diffuse) reflector (diffuse solar glint), we have

is the diffuse irradiance reflected by . Since we consider clear-sky conditions, we neglect since it is a small fraction of . To reduce computation time, at all zenith angles are precomputed and the value corresponding to a given zenith angle is obtained by interpolation of the pre-computed 's.

5.2. Radiance at Sea Level (One Pixel)

The radiance leaving an open-sea pixel corresponding to the IFOV of the sensor is where is the area of and is a shadowing coefficient smaller than one if the satellite is not at zenith. Indeed, in this case, some sea facets are shadowed by other sea facets. Modeling of implies to resort to ray-tracing algorithms (highly time-consuming). Here, we use a simplified, but realistic expression [12], that is, where is the error function, , where is the satellite look angle and is the RMS slope of the facets [23], that is, , where is the average wind speed at  m above sea level.

5.3. Radiance at the Spaceborne Sensor

To obtain the radiance arriving at the entrance of the spaceborne sensor, we multiply by the solid angle of the sensor (using the radius of the entrance pupil and the satellite height ). We obtain where represents the radiance received on the path between the sea surface and the sensor and is the atmospheric transmittance. For MWIR sensors, represents the radiance emitted by the atmosphere on the path between and the sensor. It can then be modeled as the integral of a blackbody with height-dependent temperature. We then approximate using a blackbody at a temperature being the mean of the air temperature along the path to the sensor.

6. Signal Displayed by the Sensor

We describe the model for converting to the signal displayed at each pixel of the sensor.

6.1. Model for the Displayed Signal (No Noise)

is transferred by the sensor optics to the detector focal plane where the image is formed. The spectral irradiance at the entry of a detector located on the optical axis is related to by the camera equation [15] where is the optical system transmittance (often % and nearly flat), is the -number. Then, the detector converts collected photons in an electrical current (photo-electric effect). The efficiency of this transformation is . Next, is spectrally filtered by the spectral response . The resulting signal is the integration of over the spectral interval corresponding to the bandpass of the detector. To increase the SNR, the signal is temporally integrated over a time interval (integration time) specified by . Hence, Here, we assume that and for all . is expressed in Coulomb (). Dividing by the electrical charge of an electron, we get the number of electrons collected by the detector, that is, .

If the imaged scene is a point source, the image produced at the detector is a blurred point due to diffraction. The resulting image is called the Point Spread Function (PSF) . For any other imaged scene, the signal at each pixel on the detector plane is given by a convolution of with , that is, For real systems, the PSF also includes nonideal effects. With each effect, a PSF is associated. The global PSF is the convolution of all PSFs. Typical nonidealities are the following. First, the optics induce blurring by the optical PSF as explained above. The image formed by the optics may move during the integration time; this introduces image motion PSF (also called smearing PSF). High-frequency (resp., low frequency) vibrations of the satellite also imply a degradation of the signal. We then associate to these vibrations a jitter (resp., pointing) PSF. The detector also adds additional blurring due to the detector PSF. Finally, the detected signal is further degraded by the electronics PSF.

Computation of (21) is computationally intensive. One alternative is to compute the FT of the PSFs. The convolution becomes a product. The FT of a PSF is called a Modulation Transfer Function (MTF). Hence, where (resp., ) is the FT (resp., inverse FT) of . In practice, is a discrete function since the number of detectors in the detector plane is finite. represents the current received at detector located at position . Similarly, is the number of electrons collected at .

6.2. Inclusion of Noise

So far, the proposed model for does not include noise present in the detector. In EO sensors, the most important noise sources are the following: () the photon (shot) noise associated with the nonequilibrium conditions in a potential energy barrier of a photovoltaic detector through which a dc current flows; () the thermal (Johnson) noise associated to fluctuations in the voltage current caused by the thermal motion of charge carriers in resistive materials, () the multiplexer (read out) noise.

Each noise is modeled as a random process (RP) with zero mean and variance (photon noise), (thermal noise), or (multiplexer noise). Models for these variances can be found in [24]. Each noise is expressed in number of electrons. The total detector noise variance is then the sum of , , and , so that

There exist two other noise sources (the quantization noise and bit errors) [3]. However, they are not considered here. Notice that these noises only modify the value of ; the reasoning remains unchanged. The signal displayed by the sensor at detector is then where is a realization of the zero-mean Gaussian RP with variance .

7. Real-Time Simulator

Evaluating sensor performance implies first to simulate for all detectors in the detector plane. This is computationally intensive due to the inclusion of the PSF (or MTF). Hence, evaluating sensor performance using this approach is not possible in real-time. Below, we propose an efficient, real-time strategy.

First, observe that, for an image in the open-sea, we have three classes of pixels: () pixels only composed of open-sea radiance, () pixels only composed of wake radiance, and () mixed pixels composed partially of open-sea radiance and of wake radiance. For each class of pixels, we propose below an RP for the received signal. Hence, instead of simulating the entire image, we only have to find a model for the pdf of the three classes of pixels. Indeed, the entire image is found by considering realizations of these three RPs. To summarize, we propose to reduce the computation of the entire image to the computation of three pdfs, one for each class of pixels. ROC curves are then obtained as discussed in Section 3.

7.1. Probability Density Function for

We first consider the pdf of an open-sea pixel. Then, we consider a wake pixel and finally, a mixed pixel.

7.1.1. Open-Sea Pdf

The signal corresponding to an open-sea pixel is denoted and is given by (20) using the geometric model of Section 4.1. mainly depends on satellite position , sun location , and wind speed . Consider that , , and are fixed. Consider a great open-sea area divided in small planar facets for which we compute sea heights (see Section 4). We then compute the received for each facet, and we plot the corresponding histogram. This gives an idea of the pdf of for open-sea pixels. Results are shown for various sea states (various wind speeds) in Figure 8 for MWIR sensors. Figure 8 shows normalized 's, denoted as , that is, , obtained as where and are, respectively, the minimum and the maximum values of for all sea facets.

Figure 8: Histogram of 's and corresponding beta distributions.

Histograms of all have the shape of a beta statistical distribution. The pdf of this distribution has two free parameters and and is given by where where is the gamma function. To find the that best fits the 's, we estimate and using the mean and the variance of the 's. We have [25] Hence, the pdf of the 's is , where is replaced by using (25) and and are, respectively, replaced by and , that is, The reason why the open-sea 's can be modeled as a beta distribution is currently not well understood.

7.1.2. Turbulent Wake Pdf

We consider a model for the pdf of the signal corresponding to a wake pixel given by (20) with the geometrical model of Section 4.2. We consider that , , and are fixed. In Section 4, we saw that corresponds to the radiance of a flat sea with important roughness. This roughness is modeled by dividing the wake pixel in microfacets with arbitrary orientation. (For computing , we consider that the wake surface temperature is equal to . However, in practice, [26]. Modeling this temperature difference is outside the scope of this paper.) The simplest model for thus considers a uniform pdf for the emissivity leading to a uniform distribution of . However, this is not realistic since having microfacets with arbitrary orientation is more probable than having microfacets with horizontal orientation. One solution is to use a beta distribution with high probability density near the signal corresponding to microfacets with orientation uniformly distributed between and (denoted as ) and a very small probability density near the signal corresponding to a flat sea (denoted as ). If , we have where and are such that and , with . Simulations show that and lead to a meaningful pdf. Figure 9 shows .

Figure 9: Pdf of : beta distribution with and .

7.1.3. Mixed Pdf

Some pixels, called mixed pixels and located at the edge of the wake, are composed of a portion of wake and a portion of open-sea (Figure 10). The signal corresponding to a mixed pixel is then where is the portion of wake signal in the pixel. We then have to find a model for the pdf of . Computing the analytical expression of the pdf of a linear combination of different beta distributions is challenging. In Section 7.2.3, we propose a method for computing this pdf. We use this method here with weights and The resulting pdf is then a beta distribution (see Section 7.2.3) given by where and are, respectively, the minimum and the maximum values of , obtained using the minimum and the maximum values of and . Coefficients and are obtained as explained in Section 7.2.3. An example of mixed pixel pdf is given in Figure 11.

Figure 10: Mixed pixels: (a) Edge and (b) overlapping mixed pixels.

Figure 11: Example pdf of a mixed pixel for , sea state , GSD of 100 m, and wake width of 40 m. Simulation corresponds to sunlight dominating (bright) wake.

7.2. Statistical Model for

Below, we describe the model for the pdf of for a mixed pixel. The approach is similar for open-sea and wake pixels. Finding a model for implies to include the effect of the PSF. We can either compute the convolution of the PSF with the image pixels or perform the FT of the image and multiply the resulting image by the MTF. Both methods are computationally intensive: they require the computation of the entire image. Below, we propose an efficient method to include the effect of the PSF without computing the entire image.

Moreover, to simulate the effect of changing the PSF (or MTF) of one particular nonideality on sensor performance, we should first be able to easily change the shape of the PSF (or MTF) and second to update ROC curves in real-time. Hence, we propose to represent each PSF (or MTF) with one scalar value: the MTF at Nyquist. This allows to rapidly update sensor performance.

7.2.1. MTF at Nyquist

If a sensor is looking at a scene, each detector of the sensor senses a pixel of size equal to the GSD. For a line of detectors, the values received at these detectors correspond to the sampling of a continuous signal corresponding to the radiance produced by all patches on the ground line corresponding to the detector line. Hence, the maximum frequency of the signal that can be sensed is . Signals with higher frequencies produce aliasing. Hence, the MTF at Nyquist frequency plays an important role in evaluating sensor performance. is often expressed using the detector size to be independent upon the GSD. Hence, . The MTF can then be characterized by one scalar value, that is, .

7.2.2. Model for the MTF

Each MTF is then described by its . For each nonideality , we have a value of , denoted . The global is the product of the 's. Sensor designers provide two MTF functions: the along-track and the across-track MTF. Both MTF are combined to give the 2D MTF. We thus have two , that is, the along-track () and the across-track (). We make the reasonable assumption that the 2D MTF is Gaussian [15]. This allows to compute the inverse FT analytically, saving computation time. Hence, where and are normalized frequency variables (with respect to ) and and are determined using and . Estimates and of and are and . is the inverse FT of , which is also a Gaussian, that is, where and are normalized detector locations (with respect to ). is used to compute in (21).

7.2.3. Model for

Consider a detector . Hence, discretizing integrals in (21), we obtain where and are the sets of pixel indexes centered on and for which has a nonnegligible value. With typical values of , the sizes of and are about to . Hence, is evaluated by summing signals of about to pixels, which is very efficient. Now, we describe a model for the pdf of . in (35) is a weighted sum of RPs. Indeed, each in (35) is the signal corresponding either to an open-sea pixel or to a wake pixel or to a mixed pixel that all are RPs. Hence, we can compute mean and variance of . We have where is the mean of . For , we have where is the variance of . Using and , we can model as a beta distribution with parameters and , respectively, given by (27) and (28), where and are, respectively, replaced by and . Hence, is