Analytic Expressions for Radar Sea Clutter WSSUS Scattering Functions †

Bello’s stochastic linear time-varying system theory has been widely used in the wireless communications literature to characterize multipath fading channel statistics. In the context of radar backscatter, this formulation allows for statistical characterization of distributed radar targets in range and Doppler using wide-sense stationary uncorrelated scattering (WSSUS) models. WSSUS models separate the channel from the effect of the waveform and receive filter, making it an ideal formulation for waveform design problems. Of particular interest in the radar waveform design community is the ability to suppress unwanted backscatter from the earth’s surface, known as clutter. Various methods for estimating WSSUS system functions have been studied in the literature, but to date no analytic expressions for radar surface clutter range-Doppler scattering functions exist. In this work we derive a frequency-selective generalization of the Jakes Doppler spectrum model, which is widely used in the wireless communications literature, adapt it for use in radar problems, and show how the maximum entropy method can be used to extend this model to account for internal clutter motion. Validation of the spectral and stationarity properties of the proposed model against a subset of the Australian Ingara sea clutter database is performed, and good agreement is shown.


Introduction
Random linear time-varying (LTV) system theory was first comprehensively described by Bello [2] and has widely been used in the wireless communications field ever since, particularly to model the multipath fading of mobile radio channels [3,4]. In particular, wide-sense stationary uncorrelated scattering (WSSUS) models, which are a subset of the category of random LTV systems, are the most common in the literature. As well, more recent work in this field has focused on nonstationary LTV communications channels [5][6][7], of which vehicle-to-vehicle (V2V) communication is a prime example [8]. A separate set of models that are distinct, but can be related to Bello's LTV theory are the Clarke/Jakes Doppler spectrum class of models [9,10], originally only applicable for flat fading (i.e., where the symbol time is much larger than the multipath delay spread). This model is also ubiquitous in the wireless communications literature, and many extensions have been proposed, such as for varying geometries [11] and more accurate fading statistics [12]. It is also commonly coupled with Bello's LTV theory (which is general enough to model frequency-selective fading where the multipath delay spread is much larger than the symbol duration) to model the Doppler component of multipath scattering [13]. Empirical [14] and analytic [15] expressions characterizing range-Doppler spreading of communication channels using LTV theory are commonplace.
Outside of communications applications, random LTV system models have been primarily used for sonar target detection and waveform design problems [16][17][18], although a few recent papers have applied this method to radar target detection problems as well [19,20]. Much of the literature on the subject treats the problem agnostically, treating radar and sonar as the same problem [21][22][23][24][25][26]. In both radar and sonar, the goal is to detect (usually) man-made targets of interest while suppressing thermal noise, which is internal to the receiver, and clutter, which are unwanted signal-dependent returns from the natural environment. Waveform design problems expressed in this form require characterization of the delay-Doppler "scattering functions" of the target and clutter, and the result is an optimized transmit waveform/receive filter pair.
In addition to waveform design problems, clutter Doppler spectrum characterization is useful for optimized moving target indication (MTI) and space-time adaptive processing (STAP) filtering, which optimize the receiver only. Specifically, the range-Doppler scattering function can be used to predict pulse-to-pulse correlations in the clutter return to suppress it using eigenfilter techniques [27,28]. An up-to-date summary of recent work on Doppler spectrum modeling can be found in [29]. Most papers on Doppler spectrum modeling do not use Bello's LTV formulation; the effect of the transmit waveform, measurement system, and clutter are usually combined. This is appropriate if the goal is to generate realistic clutter samples for simulation purposes or if the goal is to design an MTI/STAP filter because the processing is done on receive. However, waveform design requires the clutter to be partitioned separately from the transmit signal, which creates a motivation to "translate" these Doppler spectrum results to the "language" of random LTV system theory so that these Doppler spectrum models can be utilized in other problem domains.
Because scattering function estimation is a fundamental component of waveform design problems, it is a topic that has been studied for decades [30][31][32][33][34][35][36][37][38][39][40][41]. However, analytic expressions for scattering functions of sea clutter for airborne radars have not been reported in the literature except for some initial work by the author [1,42,43]. In this work we derive such analytic expressions for the case when localized internal clutter motion (ICM) is small relative to radar platform motion using a frequency-selective version of the Jakes model. Extensions to this model are then proposed to connect previous work on Doppler spectrum modeling to random LTV theory. As well, we demonstrate how an estimation of spatially-local ICM spectra can be cast as a probability density estimation problem, for which solutions can be found using the Jaynes maximum entropy (MaxEnt) method [44] and directional statistics [45]. The spectrum prediction method is then validated against a subset of the Australian Ingara medium grazing angle clutter dataset [46], and good agreement is shown.

Mathematical Background
For a monostatic radar of arbitrary polarization, the baseband backscattered signal y(t) will be modeled as the response of an LTV system h(τ, t) to a transmitted signal x(t): where τ represents downrange delay (i.e., "fast time"), t is absolute time ("slow time"), and w(t) is thermal noise. In the literature, h(τ, t) is commonly referred to as the time-varying impulse response. In this work, h will include the transceiver, antenna, atmospheric propagation, and backscattering environment, so x(t) is essentially the output of the transmit D/A, and y(t) is the input to the receive A/D, prior to matched filtering. Polarization considerations are contained in the antenna pattern and surface scattering radar cross-section (RCS). Variations in h(τ, t) with respect to t represent channel-induced Doppler shifts. In the case where the entire scenario is static, h(τ, t) = h(τ) (i.e., h is LTI), and the integral in Equation (1) reduces to a convolution. In applying the LTV model to a radar system, we use the following assumptions: • The LTV impulse response h(τ, t) is only valid over a single coherent processing interval (CPI) of N transmitted pulses • Over a single CPI, the range to target(s) is roughly constant • Relative motion produces a Doppler shift on the carrier only and does not introduce time dilation of the pulse. The conditions for this assumption to be valid are that BT < 0.1c/ṙ, where BT is the time-bandwidth product of the waveform andṙ is the maximum range rate [47] (p. 382).
We will model h(τ, t) as a wide-sense stationary uncorrelated scattering (WSSUS) process characterized by a set of related correlation functions and show how these are a function of sensor geometry, system parameters, and environmental conditions. WSSUS processes, which are a subset of all random LTV systems, are the most widely used model in channel characterization problems [3,13]. In the following sections, we will illustrate the properties of WSSUS processes and validate the core assumptions for the radar sea clutter problem.

WSSUS Processes
To apply the WSSUS assumption to a radar scattering environment, we will express the impulse response in the following form: where F −1 is the inverse Fourier operator, ρ is Doppler shift, and η(τ, ρ) is known as the "delay-Doppler spread function" [13,31]. This function is a stochastic description of the range-Doppler map of the targets. Substituting Equation (2) into Equation (1) yields: In Equation (3) it is clear to see that y(t) can be expressed as the superposition of delayed (by τ), frequency-shifted (by ρ), and scaled (by the complex gain η(τ, ρ)dτdρ) copies of x(t). It should be noted that this definition is general enough to include all scatterers in the radar's field of view-clutter as well as useful targets. In this case the overall delay-Doppler spreading function is: due to the linearity of the system. To characterize the behavior of h as a stochastic process, we will assume it is zero mean and compute its second-order statistics: Here will we apply the two fundamental assumptions of a WSSUS system: (1) The system is wide-sense stationary (WSS), implying that the autocorrelation depends only on ∆t along the t axis, and (2) the scattering at different lags τ are uncorrelated (US): The function A h (τ, ∆t) is known as the system correlation function [3,48]. By taking the Fourier transform of A h (τ, ∆t) with respect to ∆t, we can view the Doppler spectrum of the return as a function of delay, referred to in the literature as the "scattering function" S h (τ, ρ): The scattering function can be thought of as the true range-Doppler map of the scattering environment independent of the waveform used to probe it. It can also be shown that the autocorrelation function of the delay-Doppler spread function η for a WSSUS system is: It is sometimes easier to work with this form when deriving new expressions. This result shows that under a WSSUS model, the scattering not only at different ranges, but also different Doppler shifts are uncorrelated.
The full set of WSSUS system functions are shown in Figure 1, all of which are related by Fourier transforms of the different temporal variables. Knowledge of any one of these functions is a sufficient description of the second-order statistics of the system. These functions are standard tools in wireless communication channel modeling, and the diagram shown in Figure 1 can be found in virtually any textbook on the subject [3]. Figure 1. Wide-sense stationary uncorrelated scattering (WSSUS) system function relationships. Knowledge of any one of these functions grants complete knowledge of the second-order statistics of the system. Note that nowhere in this discussion have we specified a distribution for the samples of h(τ, t). If we assume that h is a circularly-symmetric complex Gaussian random process, then its distribution is completely specified by the WSSUS system functions described previously, and the magnitude envelope |y| of the return will be Rayleigh distributed [3]. However, in many cases a distribution with heavier tails, such as the Weibull or K distribution, is more appropriate, particularly for low grazing angles or when the radar can resolve individual sea spikes [49]. The results derived in the remainder of this paper, however, are valid regardless of the distribution of h. Note, however, that if h is non-Gaussian, higher-order statistics (third-order and above) are required to uniquely specify all of its properties.

Simulation Geometry
The simulation geometry is shown in Figure 2, where the airborne radar is at an altitude H above the surface traveling with velocity We will assume a coordinate system fixed to the phase center of the radar antenna, such that the antenna is at the origin, the X-axis is pointed parallel to the surface and in the plane of symmetry of the aircraft, the Y-axis is pointed out the right side of the aircraft and parallel to the surface, and the Z-axis is pointing down. In this work, we will assume that v always lies in the X-Z plane (i.e., no crabbing). Because we are seeking to model the scattering characteristics as a function of t and τ, we will need to express all spatial quantities in terms of these temporal variables.
surface airborne radar at origin X Z Y H Figure 2. Flat earth geometry in the X-Z plane. Because the goal is to ultimately produce a range-Doppler map, we need to express all quantities in terms of downrange delay τ and absolute time t.
We will model the radar signal using an approach similar to Barrick [50,51], i.e., as a spherical wave emanating from the source located at the origin as is shown in Figure 3. Each small segment of the wave reflects off of successive isorange rings on the sea surface; the total return is therefore the superposition of the returns from each isorange ring. The incremental power gain as a function of delay can be obtained using the radar range equation: where: • P T = transmit power • dP R = incremental received power from isorange ring at delay τ • φ, θ = azimuth and depression angles relative to platform • α = grazing angle, which equals θ in a flat earth model

WSSUS System Function Derivations
To create a model for the scattering function S h (τ, ρ), we need to characterize the signal return as a function of delay and Doppler. To do so, we will consider the Doppler spectrum generated from the backscatter from a single isorange ring and then apply the principle of superposition to obtain the total response from all ranges. The approach taken in this section can be considered a generalization of the Clarke model for flat fading [9,12], applied to modeling radar surface clutter.
We will assume that the impulse response from this ring is the superposition of returns from N equiangular patches in azimuth with random amplitudes and phases. We will also assume that the scatterers on this ring are located at delayτ and write the incremental impulse response of this thin ring as follows: where φ n = 2πn/N is the azimuth angle to patch n, da(φ n ) is the infinitesimal amplitude gain at angle φ n , ρ(φ n ) is the Doppler shift at angle φ n , and γ n is a random phase shift. The incremental autocorrelation function of the impulse response atτ is then given by: We will also assume that the scatterer amplitudes da(φ n ) are mutually independent, and we will assume that the amplitude da(φ n ) is independent of the Doppler shift ρ(φ n ), therefore: Note that E |da(φ n )| 2 is simply the backscattered power gain from scatterer n, therefore by using Equation (9) and scaling it to account for the fact that the surface area of each patch is smaller by a factor of 1/N, we obtain: Note that N = 2π/∆φ, where ∆φ is the angular spacing between patches, which upon substituting in Equation (13) yields: Combining this with Equation (12) and taking the limit as N → ∞ yields Note that the expression E e j2πρ φ,θ ∆t is the characteristic function of the random angular frequency 2πρ φ,θ [52]. In this work, we will define the characteristic function k(∆t|φ, θ) as: where p(ρ|φ, θ) is the probability density function (pdf) of the random frequency ρ φ,θ . Substituting this expression in Equation (15) yields: One interpretation of the function p(ρ|φ, θ) that is used in the wireless communications literature is that the pdf of the random Doppler shift can be thought of as a normalized power spectral density (PSD) [3,8]. Thus if one has some model of the local Doppler spectrum, this information is accounted for in Equation (16).
To obtain the autocorrelation function for all delays τ we will use the US property that the return from each isorange ringτ is uncorrelated and thus we can integrate overτ to obtain the total response: It is clear from Equation (18) that A h (τ, ∆t) is therefore: To find the clutter scattering function S h (τ, ρ), we take the Fourier transform of A h (τ, ∆t) with respect to ∆t: The expression in Equation (20) is significant because it gives an analytical expression for the full clutter spectrum, not just the mainlobe clutter, for any antenna pattern and provides a mechanism for supplying a priori information about the local Doppler spectra as a function of look angle. In this integral it can be seen that the antenna pattern has the effect of performing a weighted average of the local Doppler spectra over azimuth.

Important Special Cases
In this section, we will use Equation (20) to derive analytic expressions for the scattering function for several useful cases.

No ICM
In the degenerate case where there is no ICM, i.e., the surface motion is small relative to the platform motion, then the localized Doppler shifts ρ φ,θ are deterministic and only depend on the look angle relative to platform motion. This means that the local Doppler spectra are each just an impulse, and the scattering function Equation (20) reduces to [42]: where the nominal Doppler shift β(φ, θ) due to platform motion is given as: where ρ X = 2v X /λ and ρ Z = 2v Z /λ are the Doppler contributions due to motion in the X − Z plane, ρ X (τ) = ρ X cos(θ(τ)) and ρ Z (τ) = ρ Z sin(θ(τ)) are the effective Doppler contributions in the observation direction, and φ (τ, ρ) = arccos((ρ − ρ Z (τ))/ρ X (τ)) is the nominal azimuth angle. The result in Equation (21) is a frequency-selective generalization of the Clarke/Jakes spectrum for airborne radar geometries [9,10]. A constant-τ cut of Equation (21) without the antenna pattern is plotted in Figure 4, with G = 1 and ρ Z = 0. The singularity of the scattering function at ρ − ρ Z (τ) ≈ ρ X (τ) is extremely critical for modeling nose-aspect clutter, as it serves to narrow the Doppler spread of the antenna pattern. For side-looking antennas, however, the singularities are irrelevant as the Jakes spectrum is flat near  (21) with the antenna pattern removed (i.e., G = 1) and no vertical motion (ρ Z = 0). The true scattering function will be windowed by the antenna pattern to focus on a specific region of ρ-a side-looking antenna will be focused near ρ = 0, whereas a nose-aspect antenna will be focused near ρ = ρ X .

Side-Looking Antenna, Level Flight Path
In many airborne radars, the antenna broadside is pointed at φ = −π/2. In this case, we will create a new angular variable ψ = φ + π/2 to represent angular deviation from broadside. The antenna pattern and local Doppler spectrum with respect to this coordinate will be denotedG andp, respectively. If we assume a level flight path (i.e., ρ Z = 0), then: which, for an antenna with beamwidth of less than 20 • , is a very good approximation. Substituting into Equation (20) to compute the scattering function yields: where the second equality is due to the fact that the antenna gain is nearly zero near the edges of the angular limit, so extending the limits of the integral to infinity does not change the result, and the third equality converts the integral from the angular domain to the Doppler domain. If we simplify the integrals by defining G(β|τ) =G(β/ρ X (τ), θ(τ)) and p(ρ|β, τ) =p(ρ|β/ρ X (τ), θ(τ)), which are just the antenna pattern and local Doppler spectrum projected into delay-Doppler space, we can see that the scattering function is just a linear transformation of the antenna pattern with the local Doppler spectrum acting as the kernel function: Note that p(ρ|β, τ) is the Doppler spectra observed by a moving platform, meaning that its center frequency is being modulated by the platform Doppler shift β(φ, θ). If we assume that p does not change much relative to β, which is a reasonable assumption due to the fact that radar antennas usually have a very small azimuthal beamwidth, we can express p in terms of a "baseband" Doppler spectrum, i.e., the Doppler spectrum caused only by ICM that would be seen by a stationary radar with an infinitesimal beamwidth. The ICM spectrum is usually what is discussed in papers on Doppler spectrum modeling (e.g., [28,29,53]). We will denote this quantity as b(ρ|τ) and note that under the narrow beamwidth assumption, p(ρ|β, τ) ≈ b(ρ − β|τ), therefore: If we note that the coefficients outside the brackets only depend on τ, we can simplify Equation (27): The relationship in Equation (27) is a common model of the relationship between ICM spectrum, antenna pattern, and observed Doppler spectrum used in side-looking airborne radars (e.g., [54]); this derivation explicitly shows the conditions that are required for this model to be accurate.
We can also define a baseband analogue to the characteristic function of Equation (16) as: and note that for a side-looking radar: where A (platform) h is the autocorrelation due to platform motion only. Because k serves the function of windowing the correlation function (which in this case is also a covariance because the random process is zero mean) in ∆t, it is referred to in the STAP literature as a covariance matrix taper (CMT) [27,28].

Arbitrary Orientation
In general, for a narrow-azimuthal beamwidth antenna, the scattering function will be: where S (no ICM) h is the scattering function defined in Equation (21), which includes the Jakes spectrum peaks, and b is the ICM spectrum.

Output Time-Frequency Power Distribution
We will assume the matched-filter output is the correlation of N pulses against an infinite pulse train input, represented by the periodic ambiguity function (PAF) |χ NT |, given by [55]: where x(t) is the normalized unit-energy input pulse train and T r is the pulse repetition interval (PRI). The pulse repetition frequency (PRF) F r = 1/T r . The PAF has the following properties that are relevant to our application: and: One important result from Green's work [56] was that the time-frequency power distribution is proportional to the convolution of the ambiguity function with the WSSUS system scattering function: where E x is the energy of the pulse train. Note that since |χ NT | 2 is periodic along the τ axis, the result of the convolution integral of Equation (34) is periodic along this axis as well. If we assume the support of S h is limited to [0, QT r ] for some positive integer Q along the τ axis, then the output power distribution can be written as a circular convolution: where represents circular convolution and * represents ordinary convolution. The problem was cast into this form so that discrete implementations can use the Fast Fourier Transform (FFT) to efficiently compute the periodic output power distribution along the τ axis without zero padding.

Maximum Entropy Prior
As it has been noted in Section 2.4, b(ρ|τ) can be interpreted simultaneously as either (a) a PSD corresponding to a small patch of sea surface, or (b) a pdf of a random Doppler shift introduced by surface motion. Thus the spectrum characterization problem is reduced to a problem in prior probability density estimation.
Choosing a prior distribution is an important component of Bayesian estimation, and thus, there is a wide selection of literature available on the subject. If one has prior knowledge of the Doppler distribution from measurements or physical calculations, they can immediately compute the appropriate clutter taper using Equation (16). However, this is not usually the case, as the Doppler spectra for any given operating wavelength λ may depend on a multitude of factors, such as sea state, wind speed, temperature, salinity, etc.
In the absence of such detailed information, one tool for selecting an appropriate prior is the principle of maximum entropy (MaxEnt) [44]. It is based on the objectivist Bayesian philosophy that probabilities represent a state of knowledge rather than a degree of belief, and thus the selection of a prior distribution should be based on objective criteria, such as known expectations or known support of the random variable, in such a way that the amount of "assumed" information is minimized.
Formally, this is expressed as choosing some distribution function b(x) with support (a, b) of the random variable X that maximizes the differential entropy H: subject to the constraints imposed by the known "testable information" F i : The solution to this optimization problem is found using the calculus of variations: where c 0 is a normalization constant, m(x) is a partition function that is constant in the support region x ∈ (a, b) and zero elsewhere, and λ i are the Lagrange multipliers. The shape of the resulting distribution depends on the amount of known "testable information".

Known Mean, Unknown Variance
We often know the wind speed v w and direction φ w but no other information about the shape of the Doppler prior. If we start with the following assumptions:

•
The mean wave speed (and hence mean Doppler shift) is proportional to the wind speed. For example, it is commonly assumed that the wave speed is 1/8th the wind speed [57], and • The wind velocity vector has no Z component, and • The waves move in the same direction as the wind, then the following mathematical statements can then be constructed to express this testable information: There are cases when these statements are not true, such as when there is a sudden change of wind direction, but for a fully developed wind-wave these are not particularly controversial statements [49]. The MaxEnt prior that satisfies these criteria can be shown to be: where u(·) is the Heaviside step function. This is an exponential distribution with respect to the wind direction. The intuition behind this is that, while the mean Doppler shift will be proportional to the wind speed, there will be some wave components that move much faster, but none that move in the opposite direction of the wind.

Known Mean and Variance
If the mean E [ρ] = µ ρ and variance E (ρ − E [ρ]) 2 = σ 2 ρ of the Doppler frequency are known, then the MaxEnt prior is a truncated Gaussian: In general µ ρ = µ ρ and σ 2 ρ = σ 2 ρ due to the truncation of the tails of the distribution by m(ρ). The Doppler spectra model given in §3.8.2 of [49] can be considered a version of this with a → −∞ and b → ∞, in which case µ ρ = µ ρ and σ ρ = σ ρ .

Distribution Comparison
The claimed exponential shape of the MaxEnt Doppler spectrum may seem puzzling to experienced readers who have experience with Doppler spectra being commonly modeled as having a Gaussian shape. The reader is reminded that Equation (39) is not the spectrum the radar "sees"; the shape of the spectrum a (side-looking) radar will see is determined by the convolution (along ρ) of the Doppler spectra b(ρ|τ) with the antenna pattern G 2 (ρ|τ), which is then convolved (along τ and ρ) with the waveform ambiguity squared |χ NT (τ, ρ)| 2 . Each of these processing steps alters the resolution of the radar.
As an illustrative example, we will consider the Doppler spectra modeled in §3.8.2 of [49]. Since the radar is stationary, G 2 (ρ|τ) ∝ δ(ρ). The predicted spectrum in this scenario is described using the following parameters: Overall clutter to noise ratio (CNR) = 20 dB The computed MaxEnt spectra for both the unknown and known variance cases are shown in Figure 5. The predicted power distribution P(ρ) = |χ NT (0, ρ)| 2 * b(ρ) is shown in Figure 6. It can be seen in Figure 5 that the specification of the variance narrows the Doppler pdf b(ρ) significantly, but if too short of a CPI is used, as is the case in Figure 6, the observed spectrum after matched filtering will not change significantly.   Note that in this case the predicted spectrum when the variance is known becomes waveform-limited because σ ρ F r /N.

Doppler Spectrum Modeling
To validate our framework, we will compare the modeled clutter spectra to empirically measured spectra. For our comparison, we will use the Ingara dataset, which is a medium grazing angle dataset containing measured returns from an airborne radar flying in a circular path at a speed of approximately 200 knots, operating in spotlight mode (i.e., illuminating the patch of sea at the center of the circle.) Data were collected at grazing angles from 15 to 45 degrees. The radar operates at a center frequency of 10.1 GHz using a linear frequency modulated (LFM) waveform with a pulse width of 20 µs and a bandwidth 200 MHz, leading to a range resolution of about 0.75 m.
The data in the flight that we will use for comparison was collected from an altitude of 0.5 nmi with the plane flying in a circular path of radius 1.9 nmi, leading to a grazing angle of approximately 15 degrees. This run contains the measured radar IQ returns over a 305 second interval, or about 1.4 revolutions around the circular path. A plot of the downrange video at each slow time is shown in Figure 7. Close inspection of this picture shows that the individual wave crests and troughs are clearly visible.
To test the ability of the WSSUS framework to predict Doppler spectra, we will compare the predicted spectra to the Doppler spectra observed while the radar is looking downwind. The observed range-Doppler map seen when the radar is looking downwind is shown in Figure 8. It can be seen that the spectrum is biased in the negative Doppler direction, as expected. Significant variation in the mean Doppler frequency in each range bin can be seen; this is likely due to the high range resolution of the radar. Because these minute fluctuations are essentially random and difficult to predict, we will take the ensemble average of the normalized spectra over τ. In this dataset, the elevation beampattern was removed as a preprocessing step, so it is expected that the Doppler spectra in each range bin will have similar statistics with the exception of path loss-induced amplitude decay. The ensemble average spectrum is shown in Figure 9, where it is labeled "empirical".
Because the Ingara radar is side-looking and is flying at a constant altitude, we will use the simplifications of the WSSUS model in Section 2.5.2 to compute the scattering function S h , which we will convolve with the ambiguity to obtain the delay-Doppler map P(τ, ρ). The azimuth antenna pattern and waveform parameterization were supplied with the dataset, so the only unknown in our prediction efforts is the ICM Doppler prior b(ρ|τ).  . Empirical ensemble average spectrum plotted versus WSSUS predicted spectra for the unknown variance (u.v.) and known variance (k.v.) cases. As well, the "standard" Gaussian clutter spectrum when the variance is known is plotted for reference. It can be seen that the Gaussian spectrum model vastly underestimates the clutter floor caused by antenna and waveform sidelobes.
We will use the MaxEnt procedure of Section 2.7.1 to estimate the ICM prior. For these models, it is necessary to predict mean and variance. Because the spectrum is periodic, we estimated the mean µ ρ Doppler shift using the circular mean [45]: where E P denotes expectation is taken with respect to the empirically observed spectrum P(ρ), normalized to unit area. From the data in Figure 9, µ ρ = −47.9 Hz, which corresponds to a radial wave speed of −0.71 m/s. We will assume that µ ρ accounts for both the mean ICM Doppler shift as well as antenna pointing errors. For estimating the variance σ 2 ρ of the prior, we need to account for the fact that, according to Equations (27) and (35), the spectral width is due to the convolution of the prior with the antenna with the waveform ambiguity, therefore each step "broadens" the spectrum. In the elementary case when all three functions are Gaussian, the observed spectral variance σ 2 obs is the sum of the component variances: where σ a and σ w are the RMS Doppler widths of the antenna pattern and waveform, respectively. This can be used to estimate the variance of the clutter prior: However, because the observed distribution P(ρ) and the waveform ambiguity cut |χ NT (0, ρ)| 2 are both periodic, the appropriate way to measure σ obs and σ w are using the circular standard deviation [45]: where E |χ| 2 is expectation with respect to the normalized ambiguity squared.
The predicted spectrum using the WSSUS framework with the MaxEnt priors, labeled "WSSUS u.v." for the unknown variance case, and "WSSUS k.v." for the known variance case, as well a plot of the "standard" Gaussian Doppler spectrum model for P(ρ), are shown in Figure 9. It can be seen that both predicted spectra capture the shape of the empirical spectrum; the specification of the variance produces a spectrum that is within approximately 1 dB of the empirical spectrum at all frequencies. The Gaussian model for P(ρ) produces excelled agreement for in the main lobe but does not predict the presence of the sidelobe clutter floor of approximately -29 dB. This could lead to excessively optimistic predictions of radar performance in the sidelobe clutter region.
It is expected that for a side-looking radar in level flight the clutter Doppler spectrum will be centered at 0 Hz. The negative frequency bias of the observed Doppler spectrum in Figure 9 could be due to: • Wave motion away from the radar, as described in Section 2. Some combination of the above effects.

Validation of WSSUS Assumption
Another validation test that was performed using the Ingara dataset was to determine the validity of the WSSUS assumption. If the clutter impulse response h(τ, t) is non-WSS in t, then the autocorrelation A h (τ, ∆t) would not just depend on the time-difference ∆t, but also the absolute time t.
To characterize the change in the autocorrelation over time, we computed the autocorrelation of the data over a sliding window of width 1 s over the entire record length of 305 s in range bin 1. The start time of each window was advanced by 10 ms, leading to a 99% overlap between windows. We will denote this function A t (∆t). To characterize the autocorrelation at each time instant t, we measure the decorrelation time T c , defined implicitly by: i.e., T c is the time it takes for the autocorrelation centered at time t to decay to 1/e from its peak value.
The measured values of T c versus time are plotted in Figure 10. tV T c PV Figure 10. Clutter decorrelation times measured over time in range bin 1. Note that the decorrelation time remains roughly constant over several intervals much longer than a CPI. This provides support for the assumption that the clutter statistics are locally WSS for sufficiently long instants in time for the assumptions in Section 2.4 to apply.
We will use the stability of T c over time as a proxy to for the stationarity properties of the clutter. It is clear that the clutter is nonstationary over long periods of time, but for short instants of time it may be approximated as being locally WSS. Ironically, even researchers that study the nonstationary nature of clutter spectra must implicitly assume the clutter signal is locally WSS for them to be able to apply the Wiener-Khinchin theorem to estimate the clutter power spectrum in the first place! So the important question is not whether or not clutter is stationary, but how long does it remain locally stationary? Our modeling so far has assumed that the stationarity duration is larger than a CPI.
To test this assumption, we created a statistic called the "stability duration" of T c , which is defined as the region of time over which T c changes less than X%. (We used 10 percent, but this is admittedly arbitrary.) Because T c = T c (t) is itself a function of time, the stability duration was computed relative to the value of T c (t) at each timestep. A histogram of the stability times over the entire record length is shown in Figure 11. It can be seen that the median stability time is about 1.5 s, with the mean being even higher. Thus for any waveform with a PRI on the order of one millisecond or less (i.e., PRF > 1 kHz, which is very common) and integrating less than 1000 pulses, this is a perfectly reasonable assumption. This assumption may be violated for low-PRF waveforms with long integration times, such as those used in surveillance radars.
T c VWDELOLW\GXUDWLRQV &RXQW PHDQ PHGLDQ KLVWRJUDP Figure 11. Histogram of stability durations. It can be seen that the decorrelation time T c is roughly constant over intervals much greater than one CPI for medium-and high-PRF waveforms.

Discussion and Conclusions
Effective suppression of unwanted sea clutter returns requires an accurate characterization of their spatio-temporal statistics. We have created a WSSUS model for sea clutter that allows for the prediction of the clutter range-Doppler spectrum and provided a mechanism by which internal clutter motion may be accounted for and quickly estimated. Validation against the Ingara medium grazing angle dataset shows agreement between the WSSUS model and the ensemble average spectrum to within 1 dB at all frequencies.
Future work could apply this model to developing signal processing techniques such as adaptive transmit waveform design and improved MTI filtering to optimize signal detection in the presence of heterogeneous clutter.

Materials and Methods
The Ingara dataset, which was used for all experiments in this paper, is not publicly available; access is controlled by the Australian Defence Science and Technology (DST) Group [46]. Data analysis methods are described fully in the paper.