Student Member, IEEE Andreas Arnold-Bos

Member, IEEE Ali Khenchaf
email: khenchal@ensieta.fr

Member, IEEE Arnaud Martin
email: martinar@ensieta.fr

Bistatic radar imaging of the marine environment. Part II: simulation and results analysis

Keywords: Marine surveillance systems, bistatic radar, bistatic scattering, radar simulation, SAR imagery

We present a bistatic, polarimetric and real aperture Marine Radar Simulator (MaRS) producing pseudo-raw radar signal. The simulation takes the main elements of the environment into account (sea temperature, salinity, wind speed). Realistic sea surfaces are generated using a two-scales model on a semi-deterministic basis, so as to be able to incorporate the presence of ship wakes. Then, the radar acquisition chain (antennas, modulation, polarization) is modeled, as well as the movements of the sensors on which uncertainties can be introduced, and ship wakes. The pseudo raw, temporal signals delivered by MaRS are further processed using, for instance, bistatic synthetic aperture beamforming. The scene itself represents the sea surface as well as ship wakes. The main points covered here are the scene discretization, the ship wake modeling and the computational cost aspects. We also present images simulated in various monostatic and bistatic configurations and discuss the results. This paper follows "Bistatic radar imaging of the marine environment. Part I: theoretical background", where much of the theory used here is recalled and developed in detail.

I. INTRODUCTION

Marine radar simulators are not a new trend, and are already well established in the literature. Those tools are valuable to validate image formation models, processing tools such as synthetic aperture focusing algorithms, and other post-processing tools to retrieve valuable information on the environment. With the renewed interest in bistatic imaging, there is a growing need of a simulation tool adapted to these configurations. Such a tool could help to understand the imaging process, predict "interesting" bistatic configurations which are not well known yet, benchmark post-processing algorithms such as synthetic aperture radar (SAR) focusing algorithms in the bistatic case, etc. To the best of our knowledge, while many detailed simulations exist in the monostatic case, no such simulator exists in the bistatic case. Developing such a tool requires to extend the theory behind each element of the acquisition chain for bistatic configurations. This aspect received detailed treatment in "Radar imaging of the marine environment. Part 1: theoretical background" [START_REF] Arnold-Bos | Bistatic radar imaging of the marine environment. Part I: theoretical background[END_REF]. These theoretical elements are now re-used here to describe the implementation of a polarimetric, bistatic Marine Radar Simulator (MaRS). This This work was supported by a grant awarded by the Regional Council of Brittany.

simulation also features ship wakes, which is another topic covered here.

Generally speaking, several categories of simulators can be established (see for instance [START_REF] Zurk | Comparison of actual and simulated synthetic aperture radar image spectra of ocean waves[END_REF] for a review). First, there are tools simulating the aspect of clutter in a radar image [START_REF] Antipov | Analysis of sea clutter data[END_REF], [START_REF] Radoi | Analysis and simulation of sea clutter at high range resolution and low grazing angles[END_REF]. Those are useful to develop, say, target detection and tracking algorithms in a noisy environment, but these are not radar simulators per se since they only try to model the aspect of the final image. Second, many publications center around the simulation of stripmap synthetic aperture radar (SAR) images. Because of their resolution and often large coverage, those images are indeed of particular operational interest for instance to detect oil spills or ship wakes [START_REF] Rey | Application of Radon transform techniques to wake detection in Seasat-A SAR images[END_REF]- [START_REF] Courmontagne | An improvement of ship wake detection based on the Radon transform[END_REF]. Since a large coverage is desirable without much computation cost, SAR image simulators do not necessarily emulate all steps of the real acquisition process and deliver synthetic images directly. The general outline of such simulators has been described by Franceschetti [START_REF] Franceschetti | SARAS: A Synthetic Aperture Radar (SAR) raw signal simulator[END_REF]. First, the scene reflectivity is modeled for a given incidence angle, taking shadows into account. The resulting reflectivity map is then convoluted with an appropriate modulation transfer function (MTF) which takes in particular the synthetic antenna irradiation pattern into account. Noise is also introduced at the physical level to introduce speckle. The convolution can be efficiently done in the Fourier domain using the fast Fourier transform (FFT) which is the point of such a process, since it is computationally efficient. The marine case has also been explored, either with a clean sea [START_REF] Franceschetti | On ocean SAR raw signal simulation[END_REF], with oil spills [START_REF] Franceschetti | SAR raw signal simulation of oil slicks in ocean environments[END_REF], or ship wakes [START_REF] Oumansour | Multifrequency SAR observation of a ship wake[END_REF], [START_REF] Cochin | Simulator of ocean scenes observed by polarimetric SAR[END_REF]. In the marine case, the modulation transfer function is more complex since the scene moves; as such, it is generally divided into three parts, the first and predominant one accounting for the radar cross-section (RCS) modulation due to tilt, the two other ones accounting for the velocity bunching process and non-linear hydrodynamic interactions. This method is, as we mentioned, fast, and is also valuable since it provides insight to perform the model inversion and retrieve e.g. the sea spectrum. The problem is that most often, a perfect platform motion is considered, yet atmospheric turbulence can defocus an airborne SAR and taking this into account with a FFTbased algorithm when simulating SAR images is difficult. A solution has however been recently proposed [START_REF] Franceschetti | Sar sensor trajectory deviations: Fourier domain formulation and extended scene simulation of raw signal[END_REF], [START_REF]Efficient simulation of airborne sar raw data of extended scenes[END_REF] for ground imaging for limited platform turbulence.

The last category is raw signal radar simulators; these emulate each element of the acquisition chain explicitly, and output raw signal as acquired by a real aperture radar. This signal can be further processed exactly as if it had been acquired on a real radar. The advantage of such a configuration is that the temporal evolution of the sensors' position is computed freely by the user and noise can be introduced on the position to simulate, for instance, the effect of atmospheric turbulence on the aircraft carrying the antenna. The periodic transmission of a pulse can then be simulated, with the shape of the pulse being also a parameter. The scene itself can be described by a list of facets, the position of which is updated when a pulse is transmitted. For each facet, the bistatic radar equation is solved and the contribution of each individual facet (a chirp appropriately attenuated, de-phased and frequencyshifted) is then added to a buffer representing the received raw signal. The disadvantage of the raw signal computation is the computational cost which explains why it has rarely been used in the past, and only has begun to be used recently [START_REF] Mori | A time-domain raw signal simulator for interferometric sar[END_REF].

In our case, we desire to simulate bistatic configurations, where the transmitter and receiver are separate. Contrarily to monostatic radar, bistatic radar can be used in much more versatile configurations. Coastal radars could be used to illuminate a scene, the receiver being aboard an aircraft; the transmitter could also be aboard a satellite; bistatic SAR (BiSAR) imaging could be performed with the transmitter and receiver flying on non-parallel flight tracks [START_REF] Comblet | Bistatic radar imaging system for sea surface target detection[END_REF], and in general, the configurations are too numerous to be enumerated. This diversity is unfortunately against efficiency: while it is probably possible to generalize fast SAR image simulators for certain BiSAR configurations, the possibility of doing this for all cases (synthetic aperture or real aperture radar) is less clear. Besides, having pseudo-raw signal also is a way to test and improve BiSAR focusing algorithms or other post-processing algorithms, something that cannot be done if the synthetic image is directly simulated. Finally, we show that, with the advent of faster computers, raw radar signal simulation begins to be a viable alternative and the way to more precise simulations, with relatively high resolutions (of the metric order of magnitude) and sufficiently decent coverage for actual applications (500×500 m) with affordable simulation times.

This paper is structured as follows. First, the general simulation workflow is presented. Section III covers in detail aspects concerning the the scene digital elevation map. In particular, considerations on the discretization steps of the scene are exposed in subsection III-C. The generation of a ship wake is brought up in section IV. The last part is devoted to sample simulation results and the analysis of the computation complexity.

II. GENERAL SIMULATION DESCRIPTION

The main elements interacting in the image formation process as reviewed in Part I of our work are shown in figure 1. The additional feature is the presence of ship wakes, which are taken into account since they are very visible on SAR images [START_REF] Stapleton | Ship wakes in radar imagery[END_REF], [START_REF] Reed | Ship wakes and their radar images[END_REF]. Indeed, they can stretch over long distances (sometimes ten kilometers) and last for a long time; as a consequence, they are often investigated as the primary means to detect ships [START_REF] Rey | Application of Radon transform techniques to wake detection in Seasat-A SAR images[END_REF], [START_REF] Courmontagne | An improvement of ship wake detection based on the Radon transform[END_REF], [START_REF] Griffin | Ship hull characteristics from surface wake synthetic aperture radar (SAR) images[END_REF]. We do not cover the simulation of a ship's radar image and the ship/sea electromagnetic interaction in this paper. The main steps of the simulation algorithm are presented in Table I. Next to each step, the sections where each point is discussed in Part I and this paper are mentioned. 1): partially simulated, non-linear wave interactions not represented, first approximation only for the turbulent wake; ( 2) not yet simulated. Some relations are omitted for the sake of brevity.

As can be seen, the simulation closely follows the main steps of the radar acquisition chain presented in Part I, with the bistatic radar equation being computed for each facet of a moving digital elevation map (DEM) representing the sea surface. This facet-based description is an important aspect since the sea surface can locally be modified by the presence of a ship wake. It is therefore impossible to represent the sea image purely on a statistical basis since some determinism must be kept. The antennas are modeled either analytically for simple apertures (rectangular, circular, elliptic), but the possibility of using real gain patterns is also available, those patterns being readable from a file stored on disk.

III. HANDLING THE SURFACE DIGITAL ELEVATION MAP

The description of the scene as a DEM raises several questions, such as i) how to generate this surface, ii) how to make it evolve in time, and iii) more crucially, what discretization step and scene surface to use, which is a recurrent topic. This latter point is especially important since a too coarse step yields non-physical results; on the other hand, a finer mesh increases computation costs. We show in particular that the minimum width of the scene evolves as a function of the square of the wind speed.

A. Surface generation

At the beginning of the simulation (t = 0), a random sea is generated using a given sea spectrum, as described at section III of Part I. We know the 2-D power spectral density (PSD) of the sea height map, S, which is a function of the 2-D ocean wave vector K and physical parameters like the wind. A common practice is then to generate a random sea for time t = 0 by filtering random noise [START_REF] Toporkov | Numerical simulations of scattering from time-varying, randomly rough surfaces[END_REF]. The practical algorithm is the following. First we generate a matrix N of random complex numbers with both real and imaginary parts uniformly distributed between 0 and 1; the size of the matrix is equal to the dimension of the DEM. Then we multiply, on a term-by-term basis, the square root of the PSD by N:

Z t=0 = √ SU (1)
After an inverse Fourier Transform, the DEM is obtained:

z t=0 (x, y) = cF -1 [Z t=0 ] (x, y) (2)
Constant c is a normalization factor; its value depends on the implementation of the Fast Fourier Transform (FFT). With the FFTW package for instance, c = ∆K d where ∆K is the sampling step for the spatial wave number of the sea, and d is the dimension of the transform (e.g. d = 2 for a 2-D sea). Non-linearities are here left out of the model, but the fast computation of sea surface taking into account some non-linearities have begun to be investigated recently, e.g. by Toporkov [START_REF] Toporkov | Numerical simulations of scattering from time-varying, randomly rough surfaces[END_REF] and Saillard et al. [START_REF] Saillard | Sea surface probing with L-band Doppler radar : experiment and theory[END_REF].

B. Evolution of the elevation map in time

Once the sea map z t=0 is known, the sea map at a given time t can be deduced from it by multiplying its Fourier transform Z t=0 by a phase factor exp(-jωt). Assuming the sea depth is infinite, the temporal pulsation ω of an individual wave is linked to the modulus K of the spatial wave vector by the following dispersion relation [START_REF] Timman | Water waves and ship hydrodynamics, ser. Mechanics of Fluids & Transport processes[END_REF]:

ω 2 = g 0 K (3)
The array containing the ω values needs only be computed once. Also, it is convenient to save and re-use Z t=0 to perform the de-phasing in the Fourier plane so as to cut down on computation cost. The orbital speeds can similarly be computed efficiently within the linear framework by using the Fourier transform.

C. Choosing the sampling steps

We denote L the width of the sea surface we simulate, and n the number of facets the side is divided in. The choice of these parameters is important: a high number of facets entails a higher simulation cost, and careless setting of L and n can entail non-physical results in the simulation. The sampling of the height map is directly linked to the sampling of the sea spectrum, by the following relations:

∆K = 2π L (4)
K max = π n L ( 5
)
where ∆K is the sampling step of the sea spectrum and K max is the maximum frequency taken from the spectrum (above K max , the map is represented statistically and not in a deterministic way). In this part, we determine indicative values for the lower bound of L (written L min ) and n L . 1) Minimal sampling step ∆K of the spectrum: The low frequency peak corresponding to the dominant swell frequency spans on a thin wavenumber interval. The energy in this interval must be captured well enough when generating the sea; otherwise, the sea would appear as having wrong roughness characteristics. There are three possibilities to do this:

• use a fixed-step sampling step ∆K of the spectrum, with a step approximation of the PSD of the sea (that is, when generating the sea, the energy is captured through a simple rectangle integration of the PSD); this is easy, but the size of the rectangles must be adapted to the variations of the spectrum, otherwise a significant fraction of the energy is lost;

• use a variable-step sampling, while still retaining the rectangle integration rule. The problem is to do the inverse Fourier transform afterwards: doing a fast Fourier transform in this case is harder to do;

• use a fixed-step sampling of the spectrum, but allocate for each frequency sample, the average value of the spectrum in the box. This requires to integrate the spectrum in the box, which must be done numerically. Still, if the sampling step of the spectrum is too low, this eventually means that the sea surface is approached by a (nearly) monochromatic surface, which can or cannot be acceptable depending on what is intended. The first solution is usable in most situations if ∆K is chosen carefully enough, as we shall see. Indeed, it so happens that all common gravity spectra (Pierson, JONSWAP, Elfouhaily) have the form:

f (K, U ) = a 3 K 3 exp - a 4 .g 2 0 K 2 .U 4 (6)
where a 3 and a 4 are adequate scalars. Function f has the following property:

∀a, K, U > 0, f (K.a 2 , U ) = f (K, a.U )/a 6 (7)
which means that if the wind changes, the logarithmic bandwidth of the spectrum does not change; if the wind U is multiplied by a the logarithmic curve of the spectrum is merely translated along its +∞ asymptote by vector [-2 log a; 6 log a]. This means that in the log space, the bandwidth does not change. This can be noticed in figure 4 of Part 1. It can also be seen in figure 2 in this paper, which also presents the notations introduced in this part. Suppose we know the -3 dB interval

[K -3 dB,1 (U ref ), K -3 dB,2 (U ref )]
for a reference wind speed U ref , the corresponding interval at a wind speed U is therefore:

K -3 dB,1 (U ) = K -3 dB,1 (U ref ) K peak (U ) K peak (U ref ) (8)
K -3 dB,2 (U ) = K -3 dB,2 (U ref ) K peak (U ) K peak (U ref ) (9)
where K peak (U ) is the wavenumber corresponding to the maximum of the spectrum. So as to capture the energy of the low-frequency (gravity waves) spike, it is possible to take the sampling step ∆K smaller to, say, α percent of the -3 dB bandwidth of the spectrum (with respect to its maximum).

Then by using eq. ( 4), it follows that:

L min = α 2π K -3 dB,2 (U ref ) -K -3 dB,1 (U ref ) U U ref 2 (10)
This means that L min of the sea surface follows a quadratic dependence on the wind speed U . So as to have an idea of the order of magnitude of L min , using Fung and Lee's gravity spectrum (equation 25 of Part 1) is convenient since it is simple. The peak is analytically found at:

K peak (U ) = 0.702g 0 /U 2 [rad/m] (11)
If we arbitrarily set the reference wind speed U ref at 10 m/s, we find numerically:

K -3 dB,1 (10) = 4.55 × 10 -2 rad/m (12) K -3 dB,2 (10) = 1.22 × 10 -1 rad/m (13)
Finally, if α = 25% (which is a good compromise):

L min = 3.28U 2 (14)
This gives values of 82.1, 328 and 739 meters for L min with wind speeds equal to 5, 10 and 15 m/s (at 19.5 m AMSL) respectively. A numerically similar result should be obtained with other spectra. Those dimensions are perfectly acceptable in terms of computation cost for most of the configurations, except perhaps for higher wind speeds. Still, the linear model being increasingly approximate as wind increases, high wind speeds should be avoided anyway. It is worth noting that the rule on L min is indicative, not imperative.

2) Upper wavenumber bound K max : The upper wavenumber bound K max directly determines the facet density, n/L. Ideally, we would like to perform a coherent sum of the contributions from each facet of the sea surface, because this would be more rigorous. To do a coherent sum, all structures represented statistically on a given facet must have a typical height h such that the phase difference of the waves reflected by two points of the surface is much smaller than 2π. If the Gravity spectrum at U ref = aU (here a>1)

Spectrum at U Gravity spectrum translation [-2log a; 6log a] K -3 dB,1 (U ref ) K -3 dB,2 (U ref ) K peak (U ref ) K -3 dB,1 (U) K -3 dB,2 (U) K peak (U) K c (U) (-1 -p ) /2 (1 -3 p ) /2 log K log 5 @ (K)
capillary spectrum Fig. 2. Notations introduced in this part for the analysis of sea spectra. The Fung and Lee capillary spectrum is used and its asymptotes are drawn.

Rayleigh criterion is used, the maximum de-phasing is π/2 which means that:

h < λ 8 cos θ g i ..
.where θ g i is the incidence angle with respect to the average plane of the surface. This relation is written in the monostatic case, since we do not need more here. We call h Rayleigh = λ/8. Using the significant wave height definition, which is one among others, the typical height contributed by waves of wavenumber greater than K is:

h(K) = 4 +∞ K S 1d (K).dK (15)
It is possible either analytically or numerically to invert this relation, i.e. to obtain a relationship which allows to get K from h(K) and in particular K Rayleigh from h(K Rayleigh ).

Then, if waves of wavenumber lesser than K Rayleigh are to be represented deterministically by the DEM, equation ( 5) can be used to get the facet density in facets per meter.

Here, any sea spectrum could be used, but since only an order of magnitude is really needed, we opt here for an analytical approach with a simple spectrum for capillary waves, namely Philips' as modified by Fung and Lee (equation 26 of Part 1), which we recall here:

S 1d (K) = a 0 (1 + 3 K2 ) K(1 + K2 ) -(p+1)/2 (16)
where a 0 = 0.875(2π) p-1 g

(1-p)/2 0 and K = K/3.63 [rad/m]. Integrating this equation analytically according to (15) is not trivial. It is however easy to see that in logarithmic space, the spectrum is nearly linear by parts. Thus we have:

S 1d (K) ≈ a 0 K-p+1 2 if K < K c ( 17
)
S 1d (K) ≈ 3a 0 K 1-3p 2 if K ≥ K c ( 18
)
where K c is the wavenumber where the two lines intersect in the log-space:

K c = K m 3 1 p-1 (19)
Those equations are easier to integrate and, barring the difference between the real spectrum and the linear asymptotes (which is small enough to be neglected), we get after integration:

h(K) = 4 2a 0 K m 9(p -1) 9 K 1-p 2 -2 √ 3 if K < K c h(K) = 4 2a 0 K m p -1 K 3 2 (1-p) if K ≥ K c (20)
If K = K c , the height is equal to:

h(K c ) = 4 2a 0 K m 9(p -1) √ 3 (21)
If we equate those expressions to h Rayleigh , we get:

K Rayleigh = K m h 2 Rayleigh (p -1) 32a 0 K m 2 3(1-p) if h Rayleigh ≤ h c K Rayleigh = K m h 2 Rayleigh (p -1) 32a 0 K m + 2 9 √ 3 2 (1-p) else (22)
Numerically, we get the results provided at figure 3, where the Rayleigh height is the one obtained for a null incidence angle. In the case of sea spectra, it appears that the typical width of a facet, such that rugosity is below the Rayleigh height, would be of the order of the carrier wavelength. It is very doubtful that using the Small Perturbations Model on a small facet would be make any physical sense, since it was initially designed for infinite surfaces, or at the very least, surfaces large in front of the electromagnetic wavelength.

For this reason, smaller facet densities need to be used and an incoherent sum shall be performed. Using notations of equation ( 11), the complex amplitude scattering matrix S is approximated by taking the term-by-term square root of the power scattering matrix Σ 0 . Then, since the sum is incoherent, an additional random phase φ(x, y) taken uniformly between 0 and 2π is added to the signal reflected by that facet. Note that the summation is incoherent spatially, but there is some degree of coherence in time, at least for short durations; otherwise, SAR imaging of the ocean would be impossible. Thus, φ(x, y) was frozen in time for our simulations since the duration a facet was illuminated was under 0.5 s. However, φ(x, y) has no reason a priori not to depend on time for longer illuminations.

Still, the facet density cannot be too small either. So as to observe a local tilt effect, the facet density must be higher than the final (eventually synthetic) radar image resolution. This image resolution depends on the choice of the radar parameters and when synthetic aperture is to be used, it is the resolution after focusing that matters. However, if the effects of speckle are to be represented in a convincing manner, the number of scatterers per resolution cell must be large enough so that speckle becomes developed. It can be easily verified that when summing the complex signal coming from six scatterers with the same mean amplitude but a random phase, then amplitude already follows a Rayleigh distribution very closely. This is a good indicative value for the number of facets per unit of surface.

IV. SHIPS AND THEIR WAKES

Ship wakes are one of the most conspicuous features generated by ships on radar images. Since ship wakes last for often many hours, they leave a trail that can extend kilometers. The image of the wake is larger than the echo of the ship itself and is even visible from space with radars, thus making them useful for ship detection. We discuss here the possibility of incorporating their presence in the simulated scene. This will locally modify the elevation map and the surface properties.

A. Physical phenomena

The wake structure can be divided into Kelvin parts and non Kelvin parts [START_REF] Reed | Ship wakes and their radar images[END_REF].

The Kelvin wake is characterized by the far-field pattern resulting from the irrotational flow around the submerged body and can be accurately predicted within the inviscid fluid theory. The theory shows that the ship-generated Kelvin waves are stationary in the ship's referential, and can be divided into two types (cf. [START_REF] Newman | Marine Hydrodynamics[END_REF], section 6.10). Transverse waves travel mostly along the xb axis (which is aligned with the ship's main axis and pointing astern); they also go at the same speed V as the ship. Their wavelength is L t = 2πV 2 /g 0 (where g 0 is the acceleration of gravity). The second set is a system of divergent waves which have a speed V d = V. cos(θ) and a wavelength L d = 2πV 2 /g 0 cos 2 (θ), where θ is the angle from xb to the wave vector. Both transverse waves and divergent waves mostly live in a cone which, according to the theory, has an aperture of 38.54 • . Note that the relations quoted here are only valid if the sea depth is infinite.

The non-Kelvin features of the wake are created by viscous friction forces. They include some phenomena such as the breaking bow and stern waves, the vortices created by the ship's boundary layers and the turbulence created by the propellers [START_REF] Reed | Ship wakes and their radar images[END_REF]. Also, the turbulent flow around the ship and the flows from breaking waves have the effect of locally compressing the thin surface film into a thin layer which breaks short wavelength (capillary) waves behind the ship and "irons" the sea where the ship passed. Since this phenomenon has approximately a one-hour life, the calm zone can stretch for kilometers. When the volume of water presents a sharp gradient of temperature i.e. when the thermocline is shallow, those vortices can generate another set of perturbations on the surface, called internal waves [START_REF] Stapleton | Ship wakes in radar imagery[END_REF].

As pointed out in the literature [START_REF] Stapleton | Ship wakes in radar imagery[END_REF], [START_REF] Reed | Ship wakes and their radar images[END_REF], [START_REF] Aksnes | SAR detection of ship and ship wakes[END_REF]- [START_REF] Zilman | Kelvin and v-like ship wakes affected by surfactants[END_REF], both Kelvin and non-Kelvin phenomena are responsible for three typical wake signatures in radar images. First, the borders of the Kelvin triangle can sometimes be visible and appear as bright lines. Next, a second set of bright arms can appear and form an angle of typically 7 to 12 • . Much debate has been going on about the origin of these bright lines. The current consensus seems to be that the inner bright arms are created by waves entering into Bragg resonance with the radio wave, which are tilted by the wake waves in such a manner that they reflect the electromagnetic wave particularly well towards the receiver [START_REF] Reed | Ship wakes and their radar images[END_REF]. Finally, dark returns behind the ship (the socalled "dark V") can be visible, especially in L-band radar and are caused by the flattening of the sea as discussed above. The visibility of all these features is not guaranteed and varies with the frequency, the wind speed and the acquisition geometry [START_REF] Hennings | Radar imaging of kelvin arms of ship wakes[END_REF].

B. Simulating the wake 1) Modeling the Kelvin wake:

To model the Kelvin wake, we assume the body's hull to be described by function y b = Y (x b , z b ). The body is also assumed to be thin so as to verify Michell's thin-ship theory [START_REF] Michell | The wave-resistance of a ship[END_REF]. Within the linear framework with inviscid fluids, the surface elevation Z w induced by a moving body at a point P(x b , y b ) in the body's frame can be seen as the linear superposition of sinusoidal waves traveling at various angles θ, as long as P is in the far field (that is, more than about two ship lengths away):

Z w (P) = R + π 2 -π 2 A(θ)e -iΩ(θ) dθ (23)
where A(θ) is the complex free wave function or Kochin function, Ω(θ) = k(θ)[x b cos θ + y b sin θ] is a phase function and k(θ) is a wave number. In infinite water depth, which we assume here, this is equal to:

k(θ) = k 0 sec 2 θ ( 24
)
with k 0 = g 0 /V 2 . As pointed out by Tuck et al. [START_REF] Tuck | Sea wave pattern evaluation, part i report, primary code and test results (surface vessels)[END_REF], when the ship has no transom stern, the amplitude function can be simplified into

A(θ) = - 2i π k 2 0 sec 4 (θ)[P (θ) + Q(θ)] (25)
with:

P (θ) = x b F (x b , θ) cos(k 0 x b sec θ)dx b ( 26
)
Integral Q is similar, using a sin() instead of a cos(). Function F is itself an integral:

F (x b , θ) = z b Y (x b , z b ) exp(k 0 z b sec 2 θ)dz b (27)
It is important to keep in mind that the theory as exposed above is only valid at large distances from the ship. The formulas can be adapted for near fields and ships with transom sterns but the result is more computation-intensive since some convenient approximations cannot be made. Radar images are often at a resolution low enough so that only the far field matters much. Furthermore, the above model does not degenerate in the near field, in the sense that the heights obtained with it are "realistic", even if they are not truly exact, and bring satisfactory radar images. As a consequence, we feel that any refinement here is probably very superfluous for most applications.

The computation of the height map induced by the ship wakes amounts to solving equation [START_REF] Fung | A semi-empirical sea-spectrum model for scattering coefficient estimation[END_REF]. When the ship hull is simple enough, for instance in the case where it is a parabolic Wigley hull, equation ( 23) can be solved analytically by using the method of the stationary phase. This gives an explicit formula for the elevation for a point P living in the Kelvin cone [START_REF] Oumansour | Modélisation de la rétrodiffusion des sillages de navire en imagerie radar polarimétrique[END_REF]. This was the approach we used in a preliminary version of our code [START_REF] Arnold-Bos | A versatile bistatic & polarimetric marine radar simulator[END_REF]. However, the formula is not valid outside of the cone (where the divergent waves, even if they are very small, still exist). With the analytic method, the elevation map presents a sharp discontinuity at the boundary of the cone, which create artefacts (that is, artificially strong returns) in the radar image.

In order to capture the continuous decay of the Kelvin waves outside of the cone, and to be able to generalize the computation to real ship hulls, we followed Tuck, Lazauskas and Scullen [START_REF] Tuck | Sea wave pattern evaluation, part i report, primary code and test results (surface vessels)[END_REF] and implemented a more precise ship wake computation module. The process begins by the numerical evaluation of F for all sections x b , then P and Q for all propagation angles θ. P and Q need only be evaluated once, and can be tabulated for later use. Then, integral ( 23) is evaluated for all points P. Special care must be taken for the evaluation of these integrals, due to the presence of the sec θ terms, which diverge near π/2. As suggested in [START_REF] Tuck | Sea wave pattern evaluation, part i report, primary code and test results (surface vessels)[END_REF], we used a Filon-like quadrature scheme to compute F , P and Q and take into account the contribution of the rapidly oscillating terms near π/2. The evaluation of P and Q is even less trivial, since the integrand is rapidly oscillating as well. To avoid using a too fine discretization step, it is possible to emulate the stationary phase method: the integration is done only on a window centered around the two points corresponding to a stationary phase; the window is tapered by an exponential fade factor as exposed in [START_REF] Tuck | On ship wave patterns and their spectra[END_REF].

In the simulations presented in this paper, we used the standard DTMB hull model 5415. Numerical models of this hull are freely available and well-studied in the literature; see for instance [START_REF] Reed | Ship wakes and their radar images[END_REF]. Figure 4 shows the surface elevation induced by this hull at two typical speeds. At the lower Froude number Fr = 0.143, transverse waves tend to be more visible than at the higher Froude number (Fr = 0.287). Those results are comparable to those found in [START_REF] Reed | Ship wakes and their radar images[END_REF].

An important note needs to be given here. The simulations presented here are in X-band, with radio wavelength of the order of 3 cm; Kelvin-induced waves of this wavelength (which are below the facet size) are very small and are implicitly approached or masked by the ocean capillary waves. More attention was devoted to X-band radar in our work, but a word must be said about L or C-band radar. For these radio bands, the Bragg waves have a much longer wavelength and those Kelvin waves that enter into Bragg resonance have a much higher amplitude which cannot be approached by the ocean spectrum only. Yet the wavelength of these waves is smaller than the tile size; furthermore, the sum of the contribution of adjacent tiles is computed incoherently. This means that a local spectrum modification must be devised to take into account Bragg interference within a single facet. This would require to take the Fourier transform of the wake amplitude. Such an approach was proposed in [START_REF] Zilman | Kelvin and v-like ship wakes affected by surfactants[END_REF], [START_REF] Zilman | Radar backscatter of a V-like ship wake from a sea surface covered by surfactants[END_REF], but for a simplified ship hull only.

2) Modeling the dark V: The origin of the turbulent dark V is still not well understood, yet it is often a very prominent feature of radar images since the turbulence decays very slowly with distance x from the stern (typically in a power x -4/5 fashion) [START_REF] Reed | Ship wakes and their radar images[END_REF]. The width W of the dark V is given by a semiempirical relation developed by Zilman, Zapolski and Marom [START_REF] Zilman | The speed and beam of a ship from its wake's SAR images[END_REF], which has been shown to have a good agreement with experiments:

W (x) = w0 x0L B 1/α B (α-1)/α x 1/α ( 28
)
where B is the beam of the ship, L its length, x0 ≈ 4, and w0 ≈ 4. α has a predicted value of 5 but may vary between in 4 and 5 in practice. This relation only describes the locus of the dark V but not its surface characteristics. Since internal turbulence have the effect of damping capillary waves behind the ship, the energy in the capillary part of the spectrum must be diminished. In the present state of the simulation as it is presented in the paper, we opted for a very simple and qualitative method, by pretending that the wind has blown more softly above the zone where the dark V lies. Without changing the height map, the reflectivity of those points in the dark V is computed as explained above, with the following changes. A slope probability corresponding to a slick sea, such as the one given by Cox & Munk, is used to compute the specular part. For the diffuse part, the Small Perturbations Method is used with a spectrum where a smaller wind speed is used and manually tuned to yield realistic looking values. The method gives reasonably good visual results, as demonstrated in the images given in the next section, at a very cheap computational cost. However, this is a first approach that does by no means reflect the true mechanisms at work. Other, more sophisticated spectrum modifications must be locally incorporated to address these phenomena (see for instance [START_REF] Reed | Ship wakes and their radar images[END_REF] for a spectrum modification caused by turbulent sources within the wake or [START_REF] Skop | Modification of a directional wave number spectra by surface currents[END_REF] for a spectral modification induced by a surface current generated by a ship). Similarly, the effects of stratified fluids with different layers of water with different characteristics (temperature, salinity) are known to change the surface elevation and models and a few models have been developed [START_REF] Radko | Ship waves in a stratified fluid[END_REF].

V. EXPERIMENTS AND RESULTS

A. Monostatic and bistatic SAR images 1) Configurations:

We present four configurations in this paper (see figure 5): two airborne monostatic SAR configurations, with an incidence angle of 45 • and 63.4 • ; as well as two bistatic SAR configurations. The first bistatic configuration is a compromise between the two aforementioned monostatic SAR configurations, with two aircraft flying on a parallel path. This kind of configuration has been experimented by the French Office National d' Études et Recherches Aérospatiales (ONERA) and the German Zentrum für Luft-und Raumfahrt (DLR) in a joint experiment in 2004 [START_REF] Dubois-Fernandez | Analysis of bistatic scattering behavior of natural surfaces[END_REF], albeit with slightly different radars characteristics and above the ground only. A second bistatic configuration is an hypothetical cooperation scheme between a coastal radar and an airborne radar. In all the configurations, the two radars are aiming at the same point, at every instant.

The characteristics we used to simulate the sensors as well as the environment are summed up in tables II and III. In the present case, the radar characteristics have been more or less inspired by the Furuno FAR 28x7 series of coastal radar with its close-range setting (save for the pulse repetition frequency). The wake is created by a DTMB 5415 hull scaled so that it has a length of 100 m; the ship goes at a speed of 4.5 m/s so as to enhance the visibility of the transverse waves. The resulting elevation map is shown on figure 6.

Figure 7 shows the elementary contribution to the bistatic radar equation [START_REF] Oumansour | Multifrequency SAR observation of a ship wake[END_REF] of Part 1 for each facet. Since the transmitter has a very low grazing angle, many facets are hidden and do not actually reflect anything. The superposition of the antenna lobes is clearly visible here, as well as the secondary lobes, which do not overlap perfectly since the transmitter and the receiver are separate. Once the raw signal has been obtained, it is fed to a post-processing chain which performs the range compression and the azimuthal compression so as to obtain the synthetic aperture image. Given the present configuration, a simple Range-Doppler algorithm is sufficient to perform the beam sharpening; a description of this algorithm in the bistatic case may be found in e.g. [START_REF] Comblet | Bistatic synthetic aperture radar imaging: Theory, simulations and validations[END_REF]. The images have been focused for stationary targets. The intensity histogram of the images is then clipped down so that the intensity corresponding to the 98% percentile becomes the maximum.

The images obtained in the monostatic configuration are shown in figures 8(a) and 8(b). They compare favorably with those published in the literature (see for instance [START_REF] Stapleton | Ship wakes in radar imagery[END_REF], [START_REF] Reed | Ship wakes and their radar images[END_REF]); however most sensors are either HH-polarized or VVpolarized; images of ship wakes acquired in a cross-polarized configuration are still very uncommon. The simulation suggests that an excellent contrast can be obtained with crosspolarized images, but that remains to be compared to real images as they become available. Also, since the strongest diffuse returns from the water are obtained with an VVpolarization, this configuration is probably the worst for wake imaging, since the specular returns of the distinctive wake waves are lost in the sea clutter. As a consequence, only the dark, turbulent V is clearly visible. Judging by other simulations not presented here, we found all these conclusions to generally hold for any configuration (bistatic or monostatic), frequency and sea state.

The image obtained for configuration 3 is shown in figure 8(c) and is not much different qualitatively from the images obtained in configuration 1 and configuration 2. In a sense, configuration 3 is only weakly bistatic since the transmitter and the receiver are not widely apart; the received image is thus a compromise between the two monostatic images and does not yield much additional information on the scene. On the other hand, configuration 4 (figure 8(d)) is more of a real bistatic nature; the use of a coastal radar implies the reflection be totally diffuse since the incidence angle is large for the transmitter. Also, since the back of the waves is not in the line of sight of the transmitter, a better contrast is obtained; this is valid for all four polarization channels. This fourth configuration shows one of the advantages of bistatic imaging: using an aircraft to acquire the return signal allows for using a synthetic aperture algorithm which yields an increased resolution (here, circa 2.25 m at the center of the image) compared to what would be obtained with the bare coastal radar (here, about 27.5 m at the center of the image); at the same time, the excellent contrast of low grazing angles is retained. On the downside, this configuration requires the coastal radar antenna be servoed on the aircraft's antenna, which is incompatible with passive, non-cooperative radar applications.

B. Analysis of speckle

It is interesting to analyze the characteristics of speckle as obtained in those images. Traditional models of speckle in the marine environment include the Rayleigh distribution, the Weibull distribution, and the K-distribution [START_REF] Ward | Compound representation of high resolution sea clutter[END_REF]; a good review of models can be found in [START_REF] Antipov | Analysis of sea clutter data[END_REF] and [START_REF] Dong | Distribution of X-band high resolution and high grazing angle sea clutter[END_REF]. Of those distributions, only the Rayleigh and the K-distribution have a physical ground. The former results from a coherent sum of multiple scatterers within the same resolution cell, with the same average amplitude but a random phase. The latter results from the composition of a Gamma probability distribution, which denotes spatially and temporally slow varying features directly linked to the texture of the sea in the radar image, and a rapidly varying component modeled by a Rayleigh distribution. This is well suited to the marine environment; as such, the K-distribution has received a great deal attention. The probability density distribution of the amplitude z of a pixel under the K distribution is given by:

p(z) = 2b Γ(ν) bz 2 ν K ν-1 (bz) ( 29
)
where b is a scale factor, ν controls the shape of the distribution, and K ν-1 is the modified Bessel function of the second sort of order ν -1. When ν is large (>10), the law rapidly converges to the Rayleigh distribution; when ν is small (¡2), the distribution is more spiky, that is, the image has a nearly uniform amplitude but has some very bright pixels. We computed the histogram of the power distribution of the images in the absence of wakes; then for each reference distribution (Rayleigh, K or Weibull) we estimated the parameters of that distribution by a 2-D optimization algorithm where the mean square distance between the measured histogram and the reference probability density function, was minimized. Finally, the acceptability of the law was evaluated using the Kolmogorov-Smirnov test. It appeared that the K law was the most appropriate model for the speckle in the simulated image, since it always passes the test with less than a 10 -5 failure probability. For each simulation, the shape parameter ν is consistently lower for HH polarizations and higher for VV polarizations: radar images in HH polarization are more spiky than VV images. This is consistent with literature [START_REF] Antipov | Analysis of sea clutter data[END_REF], [START_REF] Dong | Distribution of X-band high resolution and high grazing angle sea clutter[END_REF].

We present in figure 9 the values of ν computed with surfaces of 500 × 1500 pixels (ground resolution 1 m in azimuth and 0.75 m in range) at 10 GHz, an incidence angle θ g i of 89 • which emulates the working conditions of a coastal radar (SAR processing is only used to have a constant viewing direction on a large surface). Finally, there are an average of 3.1 facets per final resolution cell; and we checked that results do not change significantly when the facet density is increased. The environment conditions are those of table III.

This configuration was chosen so as to compare our results to the empirical model developed by Ward, Baker and Watts [START_REF] Ward | Maritime surveillance radar. part 1: Radar scattering from the ocean surface[END_REF] in the monostatic case:

log ν = 2 3 log(90 -θ g i ) + 5 8 log ∆ a ∆ r 4.2 -k m -k p (30)
where θ g i is the incidence angle in degrees (80 to 89.1 • ), ∆ a and ∆ r are the azimuth and distance resolution (in meters), k p = 1 for VV polarization and 1.7 for HH polarization; k m is a wind-dependency term set Watts and Wicks [START_REF] Watts | Empirical models for prediction in Kdistribution radar sea clutter[END_REF] experimentally as:

k m = 1 3 cos(2ψ) (31)
where ψ is the wind direction with respect to the radar direction (0 • is downwind, 180 • is upwind.) The mean level for ν in our simulations is slightly higher than what could be predicted with the model (ν = 0.0657 ± 0.0423 in HH, ν = 0.1324 ± 0.0851 in VV); this means that the clutter is slightly more Rayleigh-like than what is expected, especially for the VV channel. The discrepancy between simulation and experiment can perhaps be explained by the fact that non-linearities and breaking waves have not been taken into account in the model. However, the evolution with the wind direction is consistent with the model proposed by Watts and Wicks. It is interesting to relate this evolution to that of the RCS of the sea with the wind direction, since the evolution is the same; see for instance figure 8 in Part I [START_REF] Arnold-Bos | Bistatic radar imaging of the marine environment. Part I: theoretical background[END_REF]. Generally speaking, the wind-dependence remains the same for smaller incidence angles, the clutter in the VV channel becoming more Rayleigh like as the incidence decreases; this is found both with our simulation and with experiments [START_REF] Stacy | Polarimetric analysis of fine resolution X-band SAR sea clutter data[END_REF]. Comparing the clutter to experiments in the bistatic is harder since few actual measures exist in the literature; yet some tendencies can be already be seen in the two bistatic configuration exposed before. Configuration 3 (parallel BiSAR) has speckle characteristics that are a compromise to the two SAR configurations 1 and 2; configuration 4 (hybrid coastal-airborne BiSAR) has characteristics that are numerically comparable with what is observed with a monostatic coastal radar in our simulations, which agrees with intuition.

C. A word about time-wise computation performances

Computation time is often a major concern when designing a simulator; generally, a trade-off must be found between the computation time and the accuracy of the models used to simulate the images. In order to keep performances up, the main routines of MaRS have been coded in C++, but are then called via an integrated Lua scripting-language interpreter. The code is portable on any POSIX-compliant machine. Surprisingly enough, the computation time is much below what we imagined it would be when we began to code our simulator.

For a typical session, the following steps are done. First, a sea digital elevation map (DEM) is produced for a surface of n facets; the computation time is theoretically dominated by the complexity of the FFT, which is a O(n. log n), but the hidden constant behind n calls to the functions computing the spectra is much more important. The sea and its Fourier Transform can be stored on disk for another session. Then, a DEM is produced for the ship wake. As described above, the inner integrals P and Q are first precomputed over the ship hull, which is itself discretized on n x × n z points. The integrals are computed for n θ values of θ between 0 and π/2. Then the height is evaluated for the n facets of the DEM. This process is essentially linear in computation time, and surprisingly short as compared to other hydrodynamic codes. The wake DEM needs only be computed once, since it is stationary in time. Finally, as many pulses as necessary are "shot" on the surface. For each of the m transmitted pulses, the contribution of the n facets is computed. This is (obviously enough) the stage which takes most of the computation time, and where optimization is much needed. Trigonometric function accounted for as much as 30 % total computation time, and their tabulation led for instance to an optimization by a factor 8 of the innermost loop. The computation times for a single processor and complexities are summarized in table IV.

Undoubtedly, when the width of the map increases by a factor a, the computation time increases by a factor a 2 , which makes it unsuitable for very large computations. There are however several strategies that can help to keep the computation time manageable. The problem is highly separable, both in time as in space. Indeed, the result corresponding to the surface at a time t + ∆t does not strictly depend on the result at time t, since the use of a linear model allows to compute the position and characteristics of each element of the scene directly from the scene at t = 0. The computation can therefore be parallelized on several computers. Also, the signal returned by one facet does not depend on another facet whatsoever, and the scene can be partitioned into sub-scenes which can each be processed by a different processor, only the access to the final received signal buffer needing to be shared.

VI. FINAL DISCUSSION AND FUTURE OUTLOOKS

The great advances in computing speed make it possible today to simulate radar raw signals in an affordable time while keeping a hand on each individual element of the chain and making fewer and fewer approximations. The simulated signal can then be fed to pre-existing post-processing algorithms, such as synthetic aperture image formation. The simulator can readily be specialized so as to represent coastal radars, airborne radars, or spaceborne radars, or a mixture of all these. However, even if raw radar simulation is possible, computations are still not instantaneous. Therefore, as far as BiSAR imaging is concerned, there is a clear niche for MTFbased simulations as explained in the Introduction, but this calls for the models to be extended for the envisaged bistatic configurations. Extending these will also provide invaluable insight for model inversion. Depending on what is sought to be inverted, some parts of our models must be improved. To that respect, we acknowledge the limits of our work in that the modeling of the turbulent ship wake and non-linearities in the scene must be improved; this is the focus of our current work.

The code we presented in part II is modular. As it has been designed, it accepts various spectra, surface characteristics (such as surface permittivity, roughness, etc), and scattering matrix computation subroutines. In fact, each facet has its own surface characteristics and points to a given scattering computation subroutine which can vary from facet to facet. It is therefore immediate to compute scenes with mixed characteristics, such as a sea surface with localized oil spills, etc, provided that the scattering properties are known. Additionally, the code is easily usable since all the basic routines are called using an embedded script language. It is therefore easy for someone to adapt the simulator for their own usage.

The simulator allows for having an extensive control on the environment and to generate a pseudo ground-truth. This is important when designing detection and tracking algorithms based on radar images. Interesting perspectives are that the the robustness of these algorithms will be able to be evaluated. Also, additional knowledge will be derived about which radar configuration brings the best results for a given algorithm.

All in all, this simulator is a first step towards answering two major questions concerning bistatic radar.

The first question concerns the operational benefits of bistatic configurations. The fact that the receiver is passive is already interesting per se, but other arguments in favor of bistatic radar have been put forth: increased resolution (in some cases), better images (again, in some cases), though the meaning of expressions "better images" and "some cases" remain vague, subjective and highly application-dependent.

The answer requires at worst, some example images to build the intuition; at best, they require metrics which can only be developed when a ground-truth is available. The immediate next step in our research will thus be to investigate in which cases bistatic configurations may be more interesting than monostatic configurations. On the contrary -and said very bluntly-cases where bistatic radar is not worth the trouble will also be apparent.

The second question concerns the operational requirements and problems to solve when bistatic imaging is performed, especially when the transceiver and the receiver do not cooperate. How will the final image look like when a transmitter of opportunity does not, for instance, use chirped pulses but continuous emissions, and more importantly, how can the image be improved? What happens when the two antennas are not exactly focused on the target? What happens when the location of the transmitter is not exactly known? These uncertainties and shortcomings can now be simulated and help designers develop new algorithms to handle those cases.

Fig. 1 .

 1 Fig. 1. Simulated configuration. The arrows indicate the dependencies (some of them are omitted for the sake of brevity). (1): partially simulated, non-linear wave interactions not represented, first approximation only for the turbulent wake; (2) not yet simulated. Some relations are omitted for the sake of brevity.

Fig. 3 .

 3 Fig.3. Facet density [facets/m] as a function of the radar carrier f 0 [Hz] so that waves of wavelength smaller than a facet develop a height h Rayleigh , for different wind speeds at 19.5 m AMSL. The asymptotic form of Fung and Lee's capillary spectrum is used[START_REF] Fung | A semi-empirical sea-spectrum model for scattering coefficient estimation[END_REF].

Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 4567 Fig. 4. Surface elevation for the DTMB 5415 hull, scaled for a length of 100 m; the speed is 4.5 m/s (8.7 kts, Fr = 0.143) on the left and 9 m/s (17.5 kts, Fr = 0.287) on the right.. All figures are given in meters.

Fig. 8 .

 8 Fig.8. Simulation results for the configurations provided in figure5

  determined

Fig. 9 .

 9 Fig.9. The shape parameter ν as a distribution of azimuth angle at a 1 degree grazing angle (monostatic configuration). The fit to a cosine curve has been done by least squares.

  Wake-induced DEM: P & Q for all θ O(n θ .nx.nz) 0.2 n θ = 300, nx = 50, nz ≤ 28 Wake-induced DEM: map computation O(n.n θ ) 45 n = 512 × 512, n θ = 300

	Operation	Complexity	Typical time (s) Values used to compute time
	Generation of sea DEM	O(n. log n)	5	n = 512 × 512
	Time-shifting of sea DEM	O(n. log n)	0.13	n = 512 × 512
	Single pulse	O(n)	6.6	n = 512 × 512
	TOTAL for m pulses		3497	m = 500

TABLE IV TYPICAL

 IV COMPUTATION COMPLEXITIES (NUMBER OF ARRAY ACCESSES) AND COMPUTATION TIME ON A FIRST-GENERATION PENTIUM IV @ 3 GHZ .

ACKNOWLEDGMENT

We wish to thank Arnaud Coatanhay and Michel Legris for many stimulating discussions, as well as the anonymous reviewers for their meticulous work and excellent advice.