Investigation of joint probability density function of high-resolution VV- and HH-polarized X-band sea backscatter obtained through direct numerical simulations

Joint statistical properties of range-resolved sea backscatter at different polarizations are investigated using large Monte Carlo ensembles of numerically generated data. The simulations are based on the first-principles boundary integral equation technique in the two-dimensional (2-D) space and produce noise-free sets of backscatter corresponding to well-defined wind conditions. This study focuses on the joint probability density function (PDF) of vertically (VV) and horizontally (HH) polarized clutter intensities, investigating the behavior of the function's orientation and shape with variations in the incidence angle and wind speed. Comparisons with the results generated using coherent formulation of the popular Two-Scale Model reveal strengths and weaknesses of the latter in reproducing detailed statistical descriptors of sea clutter such as the joint VV-HH PDF.


INTRODUCTION
Sea backscatter is present in many radar applications. Depending on system purpose, it can be either a source of useful remote sensing information, or a hindrance (clutter) that masks echoes from targets. In either case, good knowledge and understanding of the backscatter properties can lead to performance improvements or can even be an enabling factor. Due to randomness of ocean surface and its backscatter, statistical quantities such as magnitude probability density functions (PDFs) and their dependence on variety of radar and environmental parameters are often of interest. Historically, those detailed statistics were derived separately for individual transmit/receive polarization combinations [1]- [3]. As radar sensors increasingly adopt dual-polarization and full-polarization capabilities, more comprehensive properties of ocean backscatter become relevant. Such characteristics include joint probability distributions of signal magnitudes in various polarization channels.
Direct numerical simulations of surface backscatter that engage large Monte Carlo ensembles have proven to be a useful, high-fidelity tool [3]. The method, based on a firstprinciples boundary integral equation technique, is computationally intensive and is limited to the 2-D space. As a result, only co-polarized components of scattering matrix corresponding to signals both transmitted and received at vertical (VV) and horizontal (HH) polarizations can be obtained. Apart from this limitation, the numerical experiment is capable of producing well-calibrated noisefree dual-channel complex sea backscatter under controlled conditions, yielding both magnitude and polarimetric phase information [4], [5]. In this work, we focus on the joint magnitude PDF, trying to parametrize it and establish the behavior of the key parameters with incidence angles (including low grazing regime) and near-surface wind speed. The results are compared with the predictions of the coherent Two-Scale scattering model by Bass and Fuks [6].

NUMERICAL SIMULATION DETAILS
The numerical scattering simulations follow the scheme outlined in [3]. Sea-like surface profiles (in 2-D space) are generated as realizations of the Gaussian random process with the power spectral density modeled after the winddriven Elfouhaily wave spectrum [7]. Several values of wind speed up to 7 m/s are considered. To account for wave-wave interactions, a non-linear transformation by Creamer et al [8] is applied yielding, in particular, steeper crests and shallower troughs.
For each such surface profile the scattered electromagnetic field is calculated exactly, using a frequency-domain boundary integral equation technique for the induced surface current. The incident field at the surface acts as a forcing term and is assumed to be a tapered plane wave. The incidence angle (measured from the vertical) varies from 20 o to 85 o . Once the induced electric current on the surface is determined, the backscattered field can be calculated as its radiation effect, i.e. by applying the appropriate propagator. The frequency-domain formulation assumes that all fields and currents are time-harmonic. Range-resolving capability is achieved by repeating the simulations for a number of frequencies (1024 in all) covering the band of more than 1 GHz around =9.3 GHz ( =3.22 cm). Application of a spectral window followed by the inverse Fourier transform yields the synthesized surface response to a short pulse with a prescribed envelope [3]. The choice of the Harris-Nutall spectral weighing provides very low sidelobes and results in the pulse width of 0.4 m in magnitude (0.3 m in power). Both VV and HH-polarized returns are produced by the solution for the same surface.
These calculations are repeated for 10000 surface realizations providing robust Monte Carlo sets for statistical analysis. The generated radar echoes are easily corrected for the incident field taper and other factors that are well known from the simulations. The calibrated output, therefore, represents complex scattering amplitudes and that relate to the normalized radar cross section [3] as 〈| | 〉 . In this work, only absolute values | | and/or intensities | | will be of interest.
Scattering calculations with the coherent Two-Scale Model (TSM) are conducted within the same framework. For each surface profile, instead of solving the integral equation for induced surface current, one essentially uses the approximate expression for such a current as given in [6]. Apart from that, the calculations proceed as described above. This approach provides a rigorous test for an approximate model as no other simplifications (e.g. during analytical averaging) are introduced, and scattering from the same set of surfaces is considered. The TSM expressions in [6] allow incorporating the "shadowing factor" (being 0 or 1) to account for geometric blockage of rays by wave crests; both shadowed and non-shadowed cases were implemented.

Variations with incidence angle and wind speed
Joint PDFs for signal intensities in VV and HH channels derived from the simulated data are shown in Fig. 1 for four incidence angles. The color scale is linear and stretches from zero to the maximum of the respective PDF in each panel. It is notable that when the arguments are expressed in dB units as in Fig. 1, the PDFs are obviously elongated along a straight line. Indeed the analysis using methods of differential geometry shows that the points on the "ridge" of a distribution follow a straight line remarkably well. These observations suggest describing the "2-dimensional" bivariate PDFs with a set of parameters that include the "orientation angle" of the ridge line (cf. Fig. 1) as well as the widths along and across that direction. Such parametrization allows succinct quantitative description of the PDF behavior with various factors. For example, Fig. 2 shows how the orientation of the PDF in the VV-HH plane changes with the incidence angle and the wind speed. One certainly notes the change in behavior of the solid curves around 70 o -80 o as the grazing regime is approached, with the downward trend for higher wind speeds (i.e. rougher surfaces) starting at smaller incidence angles. Fig. 3 examines the slices of the joint bivariate PDF along the direction of the "ridge" discussed above ("major axis" of the PDF footprint, shown in Fig. 1 as a dotted green arrow) and along the orthogonal direction ("minor axis" ). In fact, the shown profiles result from summation of a number of parallel slices to achieve more stable result: for minor transect average is taken over slices within 5-dB band along centered at the PDF peak, while the width of the similar band for the major transect result is taken to be one standard deviation along . For display convenience, both functions were normalized to a unit maximum value. Fig. 4 proceeds to investigate the behavior of the width of these PDF transects vs. incidence angle. This width is calculated in the standard-deviation sense (as the square root of the second central moment of a function in Fig. 3). The results for minor transect exhibit a steady growth as the incident angle increases. A higher wind speed generally yields a larger width parameter. However, near grazing (beyond 75 o ) this distinction fades. The PDF width along the "major axis" direction appears to exhibit a shallow minimum at intermediate incidence (50 o ), growing rapidly towards grazing. Again, a higher wind speed results in a greater width, but this time the differences become prominent for large  incidence angles. The curves in Fig. 3 are notably asymmetric, so such a parameter as skewness can also be of interest. Being proportional to a higher-order distribution moment (it is defined as the third central moment divided by the cube of standard deviation), the skewness is harder to estimate, so the plots in Fig. 5 may contain bigger statistical uncertainties. It is apparent that the minor-transect skewness generally grows with incident angle (although reversing the trend around 75 o -80 o , somewhat similar to the orientation angle in Fig 2). However, it is hard to reach a conclusion on any definitive trend with the wind speed here. For major transect, the skewness is negative and reaches minimum between 70 o and 80 o . The impact of the wind speed also becomes apparent towards grazing geometry, beyond 70 o . It is worth noting that the results for small incidence angles (20 o and even 30 o ) should be viewed with caution, since the joint PDF is very narrow (cf. Fig. 1) and its estimate is probably affected by the 0.5-dB size of the histogram bin.

Fitting model probability distributions
Having a few first distribution moments, one can attempt to fit PDF transects with some model probability densities [9]. An eventual goal may be to approximately represent the joint bivariate PDF as a product of two 1-D distributions. Looking at the shapes of the slices in Fig. 3, one obvious choice is the Skew normal distribution, cf. [9]. It is characterized by three parameters (location, scale and the skewness-controlling factor ) that can be calculated from mean, variance and skewness according to straightforward expressions [9]. However, the distribution cannot have the absolute value of skewness exceeding 0.9953, and as this threshold is approached the shape of the model PDF becomes increasingly lopsided. From Fig. 5 it is apparent that the skewness of the PDF slices can come close to or exceed this value, especially for the minor transect. This was a motivation to consider another model, the Gumbel distribution [9]. It may be more familiar than it sounds, as upon variable transformation from dBs to linear units it gives rise to the Weibull PDF. The distribution is defined by two parameters, location and scale (those again can be derived from the estimated mean and variance). Its skewness is constant and approximately equals 1.14. It is expected therefore that this model should provide better approximation when the observed skewness is high, e.g. from perhaps 65 o down to grazing for minor transects, cf. Fig. 5. In principle, one can switch between the two models once the estimated skewness reaches a certain threshold, say 0.85. Results of such fit are demonstrated in Fig. 3 (dashed lines). The major transect ( =-0.793) is approximated with the Skew normal distribution, while its minor counterpart ( =1.132) is fitted with the Gumbel PDF. The comparison is encouraging, but the "transition region" when | | is around 0.85-1 (too high for Skew normal and too low for Gumbel) may still be problematic. Other distributions that provide smoother control of the skewness should be considered.

Comparisons with the coherent Two-Scale Model
Joint VV-HH PDFs can be similarly evaluated from the Monte Carlo sets generated using the approximate Two-Scale Model expressions in lieu of rigorous integralequation solution. In general, the appearance and behavior of the joint PDFs are analogous to what is observed for their "exact" counterparts in Fig. 1. They are elongated along a tilted direction in the VV-HH plane, having a blade-thin shape for near-nadir incidences and becoming more spreadout towards grazing. (The "shadowed" results at low grazing Figure 3. Normalized cuts of the bivariate PDF along the "major" and "minor" axes of its footprint (cf. Fig. 1). Dashed lines show fits with Skew Normal (major transect) and Gumbel (minor transect) distributions. angles exhibit, expectedly, a local peak at very low intensity values that are consistent with the pulse sidelobe levels; it is interesting that the direct numerical simulations do not contain such deep shadowing signature). Examining the outlined distribution parameters such as orientation angle, major-and minor-axis width, etc. helps make a quantitative judgment of how well the approximate model captures all the relevant scattering mechanisms. Figs. 2, 4 and 5 show quantities derived from the TSM-based scattering simulations. First of all, the differences between non-shadowed (dashed line) and shadowed (dash-dotted line) results become apparent after 70 o , which is not unexpected; for some quantities (e.g. minor transect width in Fig. 4) any difference between the two curves can be noticed only much closer to grazing. In comparison, the deviations of the approximate results from their exact 5 m/s counterparts are generally more pronounced. Such discrepancies at low grazing angles are not surprising. Yet the differences that some parameters (e.g. orientation angle in Fig. 2 and minor transect width in Fig.4) exhibit within the supposed validity range of the TSM (say, down to 60 o ) are noteworthy.

CONCLUSION
The joint VV-HH PDF provides much more detailed description of the surface clutter than the customarily used single-polarization PDFs. Its knowledge can be consequential for polarimetric applications. Its shape and the ways it evolves with radar and environmental parameters may hold additional insight in the surface scattering mechanisms.
Our study indicates that the joint VV-HH distribution for sea backscatter acquires a simpler shape when the arguments (signal magnitudes) are expressed in dB. This shape can be characterized by a few parameters, and direct numerical simulations provide a precise tool to establish both qualitative and quantitative dependence of these descriptors on angle of incidence, near-surface wind speed and other factors. Such a study could eventually lead to development of an approximate analytical expression for the twodimensional VV-HH PDF for sea clutter. Comparisons with the predictions of the Two-Scale Model indicate that the latter may not be fully adequate in describing probabilistic properties of the sea backscatter even within its accepted validity range.

ACKNOWLEDGMENT
This work was supported by the US Naval Research Laboratory under the 6.1 and 6.