Ocean Waves in K-distributed Clutter

—Two examples of low grazing angle radar sea clutter, both well described by the compound K-distribution model, are studied. Pulse Doppler processing is applied to obtain two dimensional range-time textures for the intensity, centroid and width of the Doppler spectrum. The first example exhibits a monochromatic swell pattern, allowing phase averaging to be applied to the textures. The second example has a more typical ocean wave spectrum. The intensity textures are gamma distributed, consistent with the compound K-distribution model, but the Doppler spectrum centroid and width textures are also found to be gamma distributed. Based on this analysis, a new method for simulation of coherent radar sea clutter is proposed, where separate memoryless nonlinear transformations are applied to a simulated water surface to generate the spatially and temporally varying intensity, centroid and width of the Doppler spectrum. The method builds on the evolving Doppler spectrum model for radar sea clutter simulation and established methods for simulation of water surfaces.


I. INTRODUCTION
Good models of sea clutter are needed for radar design and performance evaluation [1], which can include operations research using human-in-the-loop simulations, run in real time.It is important that these models accurately represent the statistics of real sea clutter.For coherent radar, the Doppler spectrum of the clutter must be modeled, including its spatial and temporal variation.There are two main approaches to modeling sea clutter, physical modeling of electromagnetic scattering from the sea surface and empirical statistical modeling of recorded data.Good results have been obtained from physical modeling [2−4], but it is a technically demanding and computationally intensive undertaking, so far limited to one dimensional representations of the sea surface.The empirical approach is more practical for on-demand simulation with varying environmental conditions, radar configurations and geometries.A current example of this approach is the evolving Doppler spectrum model [1, 5−8].
The evolving Doppler spectrum model is based on a compound Gaussian representation of the clutter, where the fast-varying clutter speckle has a Gaussian amplitude distribution which is compounded with another distribution that models the slowly varying clutter texture.A variety of texture distributions have been employed [1].If the texture is gamma distributed, the overall clutter distribution is the K-distribution [8].The Doppler spectrum is obtained as the Fourier transform of the clutter speckle over a coherent processing interval (CPI).The power, centroid and bandwidth of the Doppler spectrum all vary over longer timescales and over range.In the evolving Doppler spectrum model, correlated sequences of random variates are used to represent the spatial or slow temporal variation in the clutter power and the Doppler centroid and bandwidth.The spatial or temporal autocorrelation function (ACF) is required to produce these correlated random sequences.Recently the evolving Doppler spectrum model has been applied to clutter modeling for a scanning radar [9].In this application a two-dimensional (2D) representation of the clutter texture is required, with both spatial and temporal variation.However, the 2D range-time ACF of the water surface is difficult to model, because each component in the wave spectrum propagates at a different phase velocity, so the correlations in the sea surface in range and time are not separable.Hence a different approach is needed to simulate the clutter texture.
This paper describes a careful study of two examples of radar sea clutter from the South African Council for Scientific and Industrial Research (CSIR) Fynmeet radar [10].The first example was chosen because it exhibits an almost monochromatic swell pattern.This allows phase averaging to be applied to the data, a technique more commonly used in wave tank experiments [11].The second example was chosen firstly because it has a more typical ocean wave spectrum and secondly because analysis of this particular data set led to the evolving Doppler spectrum model [5].Young et al. [12] performed a sophisticated three dimensional (3D) spectral analysis of data from a scanning marine radar to extract the wave spectrum.2D analysis is sufficient for a real aperture radar operated in staring mode, because the wave pattern is only resolved in range and time, not in cross-range, due to the low azimuth resolution [13], [14].In many studies of sea clutter from real aperture radars, only one-dimensional analysis is performed, either on a time sequence from a single range bin, or on a single range profile averaged over a short time duration.Ing et al. [15] developed a 2D texture estimation method with the aim of predicting the sea clutter texture for adaptive radar signal processing.Here the aim is to simulate representative clutter textures, not to replicate or predict a particular texture example.Hence the subspace estimation methods described in [15] are not employed, but the entire data block is processed coherently using 2D Fourier transforms, enabling the wave spectrum and modulation transfer function (MTF) to be estimated from the data.
The results of the data analysis are used to inform a new approach to sea clutter simulation for coherent radar.In the proposed method, Fourier synthesis is used to generate a timevarying range profile of the sea surface.The surface slope is generated in the same way.Fourier synthesis for simulation of water surfaces is a well established method in marine engineering and computer graphics imagery [16−20], that has also been used in the physical modeling of radar sea clutter [2−4].The novel aspect of the method presented here is the application of a memoryless nonlinear transformation S. Bocquet is with the Defence Science and Technology Group, 506 L o r i m e r S t ., F i s h e r m a n s B e n d , V I C 3 2 0 7 , A u s t r a l i a (email: stephen.bocquet@dst.defence.gov.au).
(MNLT) to the surface slope and displacement, to produce textures that represent the power, centroid and bandwidth of the clutter Doppler spectrum.Although Zhang et al. [21] applied a MNLT to produce a log normal texture from a simulated sea surface, they did not attempt to simulate the Doppler spectrum.
The paper is organised as follows.Section II describes the data analysis.In Section III the simulation method is described and illustrated with an example.Section IV contains a discussion of the analysis results and simulation method.Conclusions are drawn in Section V.

II. DATA ANALYSIS
The data analyzed in this paper were recorded during a 2006 trial conducted by the CSIR at the Overberg test range in South Africa [22].The Fynmeet radar used in the trial was vertically polarized on transmit and receive (VV), with a 9 GHz transmission frequency, variable pulse repetition frequency (PRF) and 15 m range resolution.The antenna 3 dB beam width was 1.8° in both azimuth and elevation.Radar and environmental parameters for the two data sets used are given in Table I.

A. Monochromatic swell example
The data set CFC18-012 exhibits almost monochromatic ocean swell.Fig. 1 is a range-time-intensity (RTI) image of the data.The data is well described by the compound Kdistribution model, as shown in Fig. 2, where the probability density function (PDF) of the intensity is compared to a K+noise distribution with a shape parameter ν = 0.40, estimated from the first two intensity moments [5].The z log z estimator for the K+noise distribution [23] was also tried, but this gave a slightly inferior fit to the right tail, with a shape estimate of 0.32.In both cases the noise power was estimated from the Doppler spectrum, as described below.Pulse Doppler processing was applied with a coherent processing interval (CPI) of N = 128 pulses (51 ms).The data block was divided into 1180 CPIs in each of 48 range bins, giving a total of 56640 individual Doppler spectra.A −55 dB Dolph-Chebyshev window was applied in the time domain, then a fast Fourier transform was used to obtain the Doppler spectra.Averaged over the entire data set, the Doppler spectrum has a Gaussian shape (Fig. 3), as does its Fourier transform, the short time correlation (Fig. 4).The Gaussian fits in Fig. 3 and  4 have standard deviations of 29.9 Hz and 6.2 ms respectively.The noise power was estimated from the background power level in the Doppler spectrum, corrected for the Dolph-Chebyshev window weighting, giving a clutter to noise ratio (CNR) of 12.88 dB.The pair of peaks in the power spectrum at Doppler frequencies of ±900 Hz are radar artifacts.
The intensity I m , centroid C m and width σ m of each individual Doppler spectrum were estimated as where f r is the PRF and S (n, m) is the power spectral density (PSD) in the n th Doppler bin of the m th spectrum.The sums were restricted to L = N/8 =16 Doppler bins in the centre of the spectrum, starting at bin n 0 , as indicated by the vertical lines in Fig. 3.This removal of the exoclutter part of the spectrum reduces the noise power included in I m , and increases the CNR by 9 dB.The intensity I m from the individual spectra is the local clutter and noise power.The PDF of I m is shown in Fig. 5, with fits to the gamma, inverse gamma, inverse Gaussian and log normal distributions.Each fit includes a constant level representing the noise power, indicated by the dotted vertical line at −22.1 dB in Fig. 5.The gamma distribution clearly provides the best fit to the data.The estimated shape parameter from this fit is ν = 0.24.Texture images can be formed from the intensity, centroid and width of the individual Doppler spectra (Fig. 6).It is noticeable that the width is strongly anticorrelated with the intensity.This is due to the very low CNR in the wave troughs, where the spectra consist mainly of noise.Fig. 7 shows the autocorrelation function (ACF) of the intensity texture, which consists of equally spaced parallel ridges, due to the nearly monochromatic wave spectrum.This is confirmed by the PSD of the texture, which is dominated by a single pair of peaks (Fig. 8, top).The PSD has two-fold rotational symmetry, so the negative frequency and wavenumber quadrant at lower left is a rotated copy of the top right quadrant.A second pair of peaks corresponding to the second harmonic is also visible.The bottom plot in Fig. 8 is explained below.Fig. 9 shows the ACF along the time and range axes.The unmodified texture ACF has a 'bell bottomed' appearance with rounded troughs and sharp peaks.A correlated gamma distributed random process can be produced with the application of a MNLT to a correlated Gaussian random process.For transformation of zero mean unit variance Gaussian variates ξ to unit mean gamma variates x with shape parameter ν, the MNLT takes the form where erfc(•) is the complementary error function and Tough and Ward [24] derived the relationship between the correlations in the two processes, refer to [24] or [8] for details.This relationship can be applied to transform the ACF of the gamma distributed texture to the corresponding Gaussian process ACF.The transformed ACF is much more symmetrical about zero (Fig. 9(a) and (c)).Instead of transforming the ACF, the inverse of the MNLT (2) may be applied directly to the texture data, before calculating the ACF.The inverse of ( 2) is the MNLT Fig. 9(b) and (d) shows the ACF calculated in this way, with results that are even more symmetrical.Fig. 10 shows the cumulative distribution function (CDF) of the texture.As shown in the left hand graph, this follows a gamma distribution with shape parameter ν = 0.24, indicated by the dashed curve.If a modified texture 1 + ξ is constructed using the MNLT (3), this follows a normal distribution for intensities greater than 1, as shown in the right hand graph.The dashed curve is the normal CDF for unit mean and unit variance.The lower plot in Fig. 8 shows the PSD calculated from this modified texture.The second harmonic peaks visible in the upper plot in Fig. 8 have almost disappeared.The range and time ACFs from the modified texture can be fitted to the ACF for a windowed cosine wave, where t is the time lag, ω is the angular frequency and τ is the time duration of the window.The sine terms in (4) arise from interference between the positive and negative frequency components of the PSD; these terms are negligible if τ is large compared to the wave period.In the case of the range ACF, range lag is substituted for time lag, wave number for angular frequency and range extent for time duration in (4).These fits are shown in Fig. 9(b) and (d).The fitted wavelength is 121.8 m and the period is 11.58 s.The fits are in fact the average of two components (4) with wavelengths and periods slightly above and below these values, indicating that the wave is not quite monochromatic.There is also a small delta function contribution at zero lag, corresponding to a white noise component in the texture PSD.This is due to random fluctuations in the texture, not thermal noise from the radar.
The almost monochromatic wave allows phase averaging to be applied to the data.A graduated circular shift is applied to the texture in each range bin so that the wave crests are aligned in time (Fig. 11).First the wave period is rounded to the nearest whole number of CPIs, giving an adjusted value of 11.57 s, then the data is trimmed in time so that it consists of a whole number of wave periods.The white lines in Fig. 11 mark the edges of the data block prior to circular shifting.The aligned data can then be partitioned into single wave period blocks and averaged (Fig. 12).Doppler frequencies f D can be converted to velocities V D via the relationship V D = λ t f D /2, where λ t is the radar wavelength.The phase averaged Doppler centroid (Fig. 12, top) provides a measure of the long wave orbital velocity V Orb , which is proportional to the long wave displacement η.For a long wave with angular frequency ω L , V Orb = ω L η [11].The mean Doppler shift of −0.71 Hz or −0.012 ms −1 has been subtracted; the negative sign here indicates motion away from the radar, probably due to the offshore wind.The data has been fitted to two Fourier components, shown with dark blue curves.The principal component (dashed curve) has the long wave period, but there is also a second component with half this period (dotted curve).The two components are not in phase.The amplitude of the principal component is 0.62 m and the amplitude of the second component is 0.14 m.
The RMS surface displacement from the fit is 0.45 m.If the sea surface elevation is described by a Gaussian distribution, the envelope is Rayleigh distributed and the significant wave height is 4 times the standard deviation of the surface elevation [25].Hence the RMS surface displacement gives an estimate of 1.8 m for the significant wave height, which may be compared with the measurement of 2.09 m from the wave buoy (Table I).Alternatively, estimates of 0.70 m for the RMS surface displacement and 2.8 m for the significant wave height are obtained from the Doppler centroid texture without phase averaging.Regardless of which significant wave height estimate is used, the wave steepness of ~0.02 is very low.
The phase averaged intensity is concentrated in a fraction of the wave period (Fig. 12, center).The model curve is a Gaussian with a standard deviation of 1.4 s, which gives a good fit to the data.The long wave displacement has been overlaid, showing that nearly all the backscattered power comes from the front slope of the wave.There is almost no clutter in the region plotted with a thick black line, where the local grazing angle is negative.The geometric shadow extends a bit further, as shown by the black dashes superimposed on the displacement curve.If the MNLT is applied to the intensity data before phase averaging, the intensity variation is more sinusoidal (Fig. 12, bottom).The intensity variation can now be fitted by two Fourier components, shown as dark blue curves.The two components are almost in phase.After application of the MNLT, the intensity variation is roughly proportional to the long wave gradient, but the two are not quite in phase.The antenna elevation of −1.19° (Table I) corresponds to a gradient of −0.0208, so the local grazing angle is negative where the gradient is less than this value, as indicated by the black segment of the gradient curve.
The modulation transfer function (MTF) relates the fluctuations in backscattered power to the long wave slope [26].It is estimated by cross-correlating the power with either the wave height or the horizontal component of the long wave orbital velocity.The cyclic variations in the Doppler shift over range and time are due to the long wave orbital velocity, so the MTF can be obtained from the cross-correlation of the intensity and Doppler textures.For the CFC18-012 data set the MTF has magnitude 44.5 and phase 71°.The MTF has a single value because the long wave is monochromatic in this example.As discussed by Plant [26], in general the MTF is a function of ocean wave frequency and other parameters.It is usually estimated by incoherently averaging a large number of Doppler spectra, with the coherence being an important measure of the data quality.Here the entire block of data has been processed coherently, so the coherence is exactly 1, but it needs to be recognized that this data set corresponds to just one particular set of environmental conditions at one time and location.It is not possible to make a meaningful comparison between the MTF estimated here for a grazing angle of 1°and theory.The MTF theory is based on the composite model, which breaks down at low grazing angles, because the tilt modulation in the theory diverges as the grazing angle approaches zero.The phase difference between the intensity and the principal component of the Doppler centroid fit in Fig. 12 is 70°, which almost exactly matches the MTF phase.
Fig. 13 is a scatter plot of the spectrum width σ m vs CNR, which shows that the broadest spectra have low CNR, while almost all the spectra with a CNR above 15 dB have widths less than 50 Hz.This is confirmed by the spectrum width PDF (Fig. 14).The width distribution is much narrower if the low CNR spectra are excluded.The low CNR spectra consist mainly of white noise, hence the estimated width is large.The width estimated from (1) for a spectrum consisting entirely of noise is σ noise ≈ f r L/(2√3̅ N) = 90 Hz.The model fit to the width distribution for CNR > 10 dB, shown in Fig. 14, is a gamma distribution offset by half the Doppler bin width, i.e. f r / (2N) = 9.8 Hz.This offset is the minimum spectrum width that can be estimated with the chosen CPI.The mean width of the high CNR spectra is 24.8 Hz, with a standard deviation of 6.8 Hz.

B. Wave spectrum example
We now turn to a different data set, CFC17-001, which exhibits a more typical ocean wave spectrum.Fig. 15 is a RTI plot of this data set.Watts [5] studied this data set and used it to construct what has become known as the evolving Doppler spectrum model.He analyzed data from individual range bins.Here we apply pulse Doppler processing to the entire data block, as described above for CFC18-012.The CFC17-001 data was recorded with a 5 kHz PRF, so a 256 pulse CPI was used to retain the same 51 ms CPI duration.The texture is formed from the intensity of the individual Doppler spectra.The 2D Fourier transform of the texture is shown in Fig. 16.This is often referred to as an ω -k plot, but note that the axes in Fig. 16 are cyclic frequency f = ω/(2π) and reciprocal wavelength 1/λ = k/(2π).Almost all of the power lies along the gravity wave dispersion curve.In water of depth d, the gravity wave dispersion relation between angular frequency ω and wavenumber k is where g is the acceleration due to gravity.The solid white curve is a fit to this relation, with an estimated water depth of 15 m.The dashed curve indicates the deep water dispersion relation, ω 2 = k g.The wave spectrum can be extracted from the texture PSD by summing the intensity that lies along the dispersion curve.
The frequency resolution is Δf = 1/τ, where τ = 32.6 s is the duration of the data recording, and the wavenumber resolution is Δk = 2π/R, where R = 1425 m is the range extent of the data.The PSD image consists of a series of rectangular patches lying along the dispersion curve.These overlap in frequency, because the frequency resolution is relatively coarse due to the short data recording.Hence it is best to sum over frequency, including only the power that lies within ±Δf of the dispersion curve.This gives the wave spectrum as a function of wavenumber, which can be mapped to frequency by interpolation.The resulting wave spectrum is shown in Fig. 17.The smooth curve, shown for comparison, is a Bretschneider spectrum [27], with the same modal frequency f p = 0.075 Hz as the data.The Bretschneider spectrum has the same form as the Pierson-Moskowitz spectrum, but it does not assume a fully developed sea, so the parameter f p is not necessarily related to wind speed [28].The model spectrum does not match the data very well, but one would not expect to obtain a good estimate of the wave spectrum from 33 s of data.The modulation transfer function (MTF) can also be calculated from the data.The magnitude and phase of the MTF are shown in Fig. 17, as a function of long wave frequency.The MTF peak magnitude is at the modal frequency f p and the phase is 90° at f p , indicating that the strongest radar return is from the front face of the waves.
The mean Doppler spectrum and the Doppler centroid and width distributions can all be fitted to offset gamma distributions, as shown in Fig. 18.Parameters from these fits are listed in Table II.Note that an offset gamma distribution with offset γ, scale β and shape ν has mean γ + βν and standard deviation (SD) β√ν.The mean Doppler spectrum parameters were estimated by nonlinear least squares fitting on a log scale; parameters for the centroid and width distributions were estimated using the z log z estimator for the K+noise distribution [23], with the number of looks set to infinity.Although the fit is good in each case, fitting the PSD on a log scale overestimates the centroid by 10 Hz.Spectra with a CNR below 10 dB were excluded from the width distribution, to remove the bias from thermal noise in the radar.Fig. 19 shows the time evolution of the Doppler spectrum in a single range bin.The power in the spectrum is intermittent and quasi-periodic.The spectrum centroid and width also vary, with the periodic variation being more regular on the negative Doppler frequency side and more erratic on the positive Doppler frequency side.This behavior is reflected in the asymmetry of the mean Doppler spectrum from the whole data block (Fig. 18, 20).If ζ is the clutter speckle plus noise power, the variance in the local clutter power z is where 〈〉 denotes the mean.The variance spectrum shows the variance in the local clutter power as a function of Doppler frequency.This is also asymmetric, due to the greater variability in the spectrum on the positive Doppler side.The smooth curves in Fig. 20 are a fit to the model proposed by Miller [29].In this model a Gaussian-shaped spectrum is modulated in both power and centre frequency.The power in the spectrum is gamma distributed, while the centre frequency has an offset or non-central gamma distribution.The coefficient of variation (CV) spectrum (Fig. 21) is another way to show the variation as a function of Doppler frequency.
The CV is the standard deviation in the clutter plus noise power, divided by the mean power.Thermal noise from the radar has a CV of 1, as seen in the extremities of the spectrum.The smooth curve in Fig. 21 is calculated from the model fit to the power and variance spectra in Fig. 20.

III. SIMULATION
In the evolving Doppler spectrum model [1, 5−8], the texture ACF is used to describe the variation in spectrum intensity and centre frequency with range or slow time.This is a practical approach provided that the range and time variations are considered separately.It runs into difficulty if we require both range and time variation in the model, because in general the 2D range-time ACF of the water surface is complicated.Each component in the wave spectrum propagates at a different phase velocity, so the correlations in the sea surface in range and time are not separable.The solution to this problem is to generate a set of amplitudes from the wave spectrum, propagate them forward in time at their respective phase velocities, then use an inverse Fourier transform to obtain a realization of the sea surface profile as a function of time.Here we only consider the 2D case of time variation in the range profile, but the same approach can be used for the 3D case of a time-varying surface.This Fourier synthesis method for simulating the water surface has a long history, with applications in marine engineering and computer graphics imagery [16−20].
The surface displacement, or wave height h, at a range vector r and time t is expressed as the sum of Fourier components with wavevectors k and complex amplitudes h, The wavevectors k lie in the plane of the sea surface and the range vector r is along the radar boresight, so the dot product gives the wavevector components parallel to the boresight.For a grazing angle ψ, azimuth angle θ and wave direction θ w , k.r = k r cos ψ cos(θ − θ w ).In the example presented here, ψ ~ 1°, so cos ψ ≈ 1.If the radar is operated in staring mode, the range profile is a function of only the wavevector components along the fixed antenna azimuth.Waves propagating in other directions are filtered out due to the low resolution in azimuth [13].A scanning radar records information on all wave directions [12], so two spatial dimensions and time may need to be simulated, but here we only simulate one spatial dimension as a function of time.The initial set of complex amplitudes is produced as where ξ r and ξ i are independent normal variates with zero mean and unit variance and S k (k) is the wave spectral density as a function of wavenumber.In terms of wavenumber, the Bretschneider spectrum for deep water waves is where k p is the modal or peak wavenumber.Then at time t the Fourier amplitudes are where ω(k) is the dispersion relation.The surface displacement as a function of range and time is then obtained as the real part of the inverse Fourier transform of these amplitudes.The imaginary part provides a second, independent realization of the displacement, if multiple realizations are required.Tessendorf [18] and Mitchell [19] propose using conjugate symmetric amplitudes, so that the inverse Fourier transform is real, but it is not necessary to do this.Fig. 22 shows an example of the simulated surface displacement as a function of range and time.In order to make the simulation comparable with the CFC17-001 data, the dispersion relation for a water depth of 15 m has been used in (11).The deep water wave spectrum (10) was used, on the grounds that the waves originate in deep water, but the peak wavenumber k p /(2π) = 0.0065 m -1 was obtained by numerically inverting the dispersion relation (5) for f p = 0.075 Hz and a water depth of 15 m.The dispersion relation is recovered in t h e ω -k plot from the simulation, Fig. 16 (bottom).The dimensions of the rectangular patches along the dispersion curve reflect the range extent and time duration of the simulation.If a larger block of data is simulated, the spectral density will be more sharply concentrated on the dispersion curve.The dispersion curve obtained from large blocks of real data has a finite width, due to nonlinear hydrodynamics [30], so these simulations based on linear hydrodynamics are only realistic for short durations, up to about a minute.
The backscattered power is a function of the surface gradient, which can be computed at time t as the real part of an inverse Fourier transform with the amplitudes i k h(k, t).Since the displacement is normally distributed, the gradient will be also, but for K-clutter the texture is gamma distributed.Gamma variates are obtained from the application of an MNLT (2) with the desired shape parameter ν to the surface gradient.An example of the normalized texture intensity obtained in this way is shown in Fig. 22, with ν = 0.44.Doppler spectra can now be simulated as a function of range and slow time.The aim of these simulations is to produce spectra with the same general appearance and statistics as the real ones, but not to reproduce any particular data set exactly.The individual spectra are based on a Gaussian function which gives the power at a Doppler frequency f n in the m th spectrum, for an N pulse CPI, ) , (13)   where At each range and time, the spectrum intensity I m is obtained from the MNLT applied to the simulated surface gradient, while the centroid C m and width σ m are obtained from separate MNLTs applied to the simulated displacement.Different shape parameters are used in the MNLTs for the power, centroid and width.The Doppler centroid and width are much less variable than the power, so much larger shape parameters (Table II) are needed to produce the correct distributions.A phase shift of −30° was applied to the width, so that the spectrum width peaks after the wave crests.The complex amplitude spectra are simulated as (14)   where c is the CNR and χ(n) and ρ(n) are complex normal variates representing the clutter speckle and noise respectively.
A simulated Doppler-time-intensity plot for a single range bin is shown in Fig. 19 (bottom).Both the real and simulated spectra have a distinctive tilted appearance, whereby the intermittent spectral intensity shifts to higher frequencies as it increases, before subsiding.The mean Doppler and variance spectra (Fig. 20, right) and the CV spectrum (Fig. 21) are all reproduced fairly well in the simulation.However, it needs to be recognized that there is a lot of variability in the wave spectrum as observed or simulated over a relatively short range extent and time.This is true despite the large numbers of individual Doppler spectra (61248 from the data and 98304 in the simulation), which result from the high PRF and short CPI.In order to make the simulation example reproducible, the Mersenne twister random number generator was specified in Matlab, with the random seed set to 5.This particular example was chosen because the randomly generated components in the wave spectrum were similar to those in the data.The mean Doppler spectrum is insensitive to the choice of random seed, but the variance and CV spectra change appreciably with different random seeds, especially the height of the positive Doppler peak in the CV spectrum.Nevertheless, in all cases the simulation is able to reproduce the tilted appearance of the real spectra, something that is difficult to produce with other methods.

IV. DISCUSSION
The clutter simulation method described by Ward, Tough and Watts [8,Chapter 4] is based on correlated sequences of gamma variates, obtained by the application of a MNLT to correlated normal variates.The approach is to first model the correlation in the gamma-distributed clutter texture, then deduce the ACF that must be applied to the input Gaussian sequence in order to produce this correlation in the output gamma sequence.With this approach in mind, Tough and Ward [24] devised a power series representation of the mapping between the correlation functions of the input Gaussian process and the output process.Application of this mapping to the texture correlation in the CFC18-012 data set clearly shows that it is the correlation in the underlying Gaussian process that is more regular and easier to model, not the correlation in the gamma-distributed texture.Moreover, rather than modifying the ACF, it is better to apply a MNLT to the texture first, so that the transformed data is approximately normally distributed.
The monochromatic wave pattern in the CFC18-012 data allowed phase averaging to be applied, so that the surface profile could be estimated from the Doppler centroid over one period or wavelength of the long wave.Phase averaging was also applied to the backscattered intensity, showing this to be concentrated on the forward slope of the long wave.If the intensity data is first transformed with a MNLT and then phase averaged, the transformed intensity is approximately proportional to the surface gradient.These observations motivate a different approach to the simulation of sea clutter.
It has long been known that the sea surface is well described by a Gaussian random surface [31], [25], provided that the wave steepness is not too large.Analysis of the CFC18-012 data set showed that the clutter texture could be approximately mapped to a Gaussian random surface with the application of a MNLT.This suggested that the textures required for a clutter simulation could be obtained with the application of a MNLT to a simulated sea surface.The relationship between the sea surface and the clutter texture also explains the difficulties encountered in modeling the 2D range-time ACF for the clutter texture.Water waves are dispersive, so the spatial and temporal correlations in the sea surface are not separable, but the Fourier synthesis method provides a solution to this problem.
The simulations described here are based on linear hydrodynamics.Whether or not this is satisfactory depends on the required fidelity and the temporal and spatial extent of the simulation.Stuhlmeier and Stiassnie [32] compared three methods of wave forecasting: linear dispersion, Stokes' corrected dispersion and nonlinear dispersion from numerical solution of the Zakharov equation.They found that departures from linear theory became apparent after a time ~ ϵ -1 T p , where ϵ is the wave steepness and T p is the peak wave period of the spectrum.So, for example, if T p ~ 6 s and ϵ ~ 0.1, linear dispersion should be satisfactory for a simulation lasting up to about a minute.
The proposed simulation method is easily able to reproduce the distinctive tilted appearance of the slow time-Doppler or range-Doppler images.McDonald and Cerutti-Maori [33] proposed a simulation method based on a mapping between the swell pattern and clutter 'events' to account for this, but their method does not work if the radar is looking directly up or down swell [34].The evolving Doppler spectrum model is able reproduce the tilted appearance in the simulated spectra by means of a correlation between Doppler and intensity, but only if this correlation is significant.In the case of VV polarization, the correlation is very small, regardless of look direction [6], [33], so the simulated spectra are not noticeably tilted in range or slow time, unlike the real spectra.In the method proposed here, instead of an explicit correlation, the Doppler-intensity correlation is local and implicit, due to the Doppler centroid and intensity both being generated from a common realization of the simulated ocean surface.Some modification to the proposed method will be needed for a cross-wave geometry, where the radar boresight is orthogonal to the wave direction, because as it stands, the method will not produce any slow time variation in the Doppler spectrum in this situation.Significant further development is needed for the simulation method proposed here to be practically useful, but the evolving Doppler spectrum model provides a blueprint for this.The evolving Doppler spectrum model began with simulation of clutter spectra for a ground-based radar with VV polarization, operated in staring mode [5].It was first extended to different wind and wave look directions [35] then to airborne radar with medium grazing angles and other polarizations [6], [7].Bimodal spectra were introduced at this stage, to include fast scattering.The model was then extended to airborne scanning radar and multiple phase centers [1], [9].Lastly, variations in radar resolution and integration time have been allowed for [36].Empirical model parameters need to be related to environmental conditions and geometry for the model to be useful.The model proposed here is based on a simulation of the ocean surface, so we can draw on existing simulation work for this [37], [20].The K-distribution shape parameter has already been related to environmental conditions and geometry for the evolving Doppler spectrum model [7], but the gamma distribution shape parameters used in the MNLTs for the Doppler spectrum centroid and width are new empirical parameters that will need to be modeled.
A single Doppler spectrum component is used in the model as presented here, which is satisfactory for VV polarization, but it is likely that a second component will be needed to represent breaking waves and sea spikes when the model is extended to horizontal polarization.The evolving Doppler spectrum, spatially and temporally limited clutter [33] and compound G2 [38] models all use two components to model the Doppler spectrum.The second component provides additional degrees of freedom for an accurate model, but it also complicates the model and introduces more empirical parameters that must be related to the environmental conditions if the model is to be useful [34].
The monochromatic swell seen in the CFC18-012 data set is unusual, but not unique.For example, the Starea4 data set from the IPIX radar, studied by Greco, Bordoni and Gini [39], appears very similar.The phase averaging used here could be applied to other data sets exhibiting monochromatic swell patterns.In particular, it would be useful to analyze horizontal and cross polarized data in order to model the sea spike and breaking wave components.The data analysis presented here was limited to vertical polarization, because the CSIR Fynmeet radar is vertically polarized on both transmit and receive.
It is interesting to compare the analysis results from the CFC18-012 data set with the numerical study of Johnson et al. [3].In their study, three methods for retrieval of the sea surface profile from low grazing angle radar measurements were compared using simulation.The three methods used (i) the normalized radar cross section (NRCS), (ii) interferometry and (iii) the Doppler centroid.The Doppler centroid gave the best results, with fairly accurate measurements even in shadowed regions of the surface profile.The NRCS method also provided measurements from shadowed regions, albeit less accurate than the Doppler centroid.Interferometry requires two vertically separated antennas, so it is not comparable with the Fynmeet radar data.The significant wave height measurement from the wave buoy is the only ground truth available for the CSIR data; nevertheless the Doppler centroid texture calculated from the CFC18-012 data set appears to give a good estimate of the sea surface profile.Phase averaging showed that the backscattered intensity is mainly from the forward slopes of the waves, with very weak returns from the wave troughs.Thus it is not surprising that the NRCS method does not perform so well in retrieval of the sea surface profile in a low grazing angle situation where there is significant shadowing.However, the fact that good estimates of the Doppler centroid can be obtained from shadowed regions shows that geometric shadowing is not a satisfactory model for the backscatter.This is in accord with both theory [40], [41] and much more extensive field measurements [42].
V. CONCLUSIONS Study of an example of low grazing angle sea clutter data with a monochromatic swell pattern showed that the clutter intensity texture could be transformed into a correlated Gaussian random field with the application of a MNLT.The correlation in this Gaussian random field was found to be much more regular and simpler to model than the correlation in the gamma distributed clutter texture.Similar textures were obtained from the Doppler centroid and width, with the Doppler centroid providing a measure of the surface displacement.A more typical ocean wave spectrum was obtained from 2D Fourier analysis of a second example of low grazing angle sea clutter data.In this case the mean Doppler spectrum and the Doppler centroid and width distributions could each be modeled with offset gamma distributions.
These observations led to the proposal of a new method for simulation of coherent radar sea clutter, where separate MNLTs are applied to a simulated water surface to generate the spatially and temporally varying intensity, centroid and width of the Doppler spectrum.The method builds on the evolving Doppler spectrum model and established methods for simulation of water surfaces.Further development will be needed for the method to be useful, in particular the empirical parameters in the model must be related to environmental conditions and geometry, but the evolving Doppler spectrum model provides a roadmap for this work.