From Maxwell's Equations to Polarimetric SAR Images: A Simulation Approach

A new electromagnetic approach for the simulation of polarimetric SAR images is proposed. It starts from Maxwell's equations, employs the spectral domain full-wave technique, the moment method, and the stationary phase method to compute the far electromagnetic fields scattered by multilayer structures. A multilayer structure is located at each selected position of a regular rectangular grid of coordinates, which defines the scene area under imaging. The grid is determined taking into account the elementary scatter size and SAR operational parameters, such as spatial resolution, pixel spacing, look angle and platform altitude. A two-dimensional separable “sinc” function to represent the SAR spread point function is also considered. Multifrequency sets of single-look polarimetric SAR images are generated, in L-, C- and X-bands and the images are evaluated using several measurements commonly employed in SAR data analysis. The evaluation shows that the proposed simulation process is working properly, since the obtained results are in accordance with those presented in the literature. Therefore, this new approach becomes suitable for carrying out theoretical and practical studies using polarimetric SAR images.


Introduction
Retrieval of targets' biophysical and geophysical parameters is one of the main goals of microwave remote sensing, having been the subject of numerous studies. Understanding the electromagnetic reflective properties of targets is a key to the correct interpretation of microwave imaging data. In this sense an image simulator might become a powerful tool for remote sensing researchers, since the use of simulated images may improve considerably the knowledge on several synthetic aperture radar (SAR) applications.
The applicability of synthesized images ranges from theoretical considerations to practical problems. For instance, from SAR simulated images it is possible to develop dedicated algorithms for filtering, segmenting or classifying images. Simulated images can also be used in remote sensing inversion techniques, in the identification of the scattering mechanisms intrinsic to a set of pixels, and in sensor calibration procedures, among others.
There are several ways to synthesize an image; for instance, in [1][2][3] a statistical technique is used to generate amplitude SAR data. The employed distributions are derived from multiplicative models and are associated with the homogeneity degree of each target [4]. However, to perform the simulation process as realistic as possible, the polarimetric SAR image simulation problem is addressed in this work from the point of view of electromagnetic modeling. The employed simulation approach is based on the computation of the far electric field scattered by a multilayer structure excited by a plane wave, using the moment method technique [5].
According to [6], the next challenging era of satellite programmes are the Cartwheel satellites, consisting of a transmitter and several small receivers for a specific purpose. This technique is completely based on the bistatic behavior of radar waves for different targets. Therefore, there is currently a need to perform such bistatic radar measurements for various kinds of targets and to develop a less complex model for the retrieval of the target characteristics. There are currently in the literature few works that have used the polarimetric analysis of bistatic radar data to retrieve target parameters, such as [7]. Although the SAR image simulation process proposed here can be used in either monostatic or bistatic SAR framework without any additional complexity, this paper assumes the monostatic framework to simulate polarimetric SAR images.
The paper is organized as follows. A general theory behind the electromagnetic model is outlined in Section 2. In this section the moment method is also developed and applied to a particular structure. The generation of polarimetric SAR images is detailed in Section 3, where SAR images for a simple multilayer structure are simulated. In Section 4, the simulated polarimetric images are evaluated based on intrinsic properties of amplitude and polarimetric SAR data; results of different classification approaches are also analyzed. Finally, the conclusions are drawn in Section 5.

Electromagnetic Model
The electromagnetic model is based on the determination of the electromagnetic fields scattered by a multilayer planar structure that is excited by plane waves. The structure under analysis is composed of N+2 isotropic, linear and homogenous layers stacked up in z direction. The layers are assumed to be unbounded along the x and y directions. The lower layer, having complex permittivity ε g and complex permeability µ g , is denoted as ground layer and occupies the negative-z region. The next N layers are characterized by thickness ℓ n , complex permittivity ε n and complex permeability µ n , where 1 ≤ n ≤ N.
The planar interface z = d N separates the N-th layer from free space (the upper layer). Metallic patches, which behave as scattering elements, are printed at arbitrary positions on each one of the N+1 interfaces of the structure. The development is based on a global right-handed rectangular coordinate system located on the top of the ground layer (interface z = 0) and lying on the xy-plane. The geometry of the planar multilayer structure is depicted in Figure 1.

Electromagnetic Fields in the Structure
The electromagnetic fields in a multilayer structure are determined through the methodological approach described in [8]. According to this methodology, which employs the spectral domain full-wave technique, the structure is treated as a boundary value problem, where the induced electric surface current densities on the metallic patches are the virtual sources of the scattered fields. Since the layers of the structure are free of sources, the Maxwell's equations for the n-th layer, assuming time dependence of the form where, for free space and the ground layer the index n is equal to 0 and g, respectively, ω is the angular frequency and the vectors E n (x, y, z), H n (x, y, z), D n (x, y, z), and B n (x, y, z) denote the complex electric field, magnetic field, electric flux density and magnetic flux density, respectively (bold face letters represent vectors). Using the following constitutive relations for each of these media the wave equations for the n-th layer is written as where n n n k ε µ ω 2 2 = gives the wave number in the n-th layer. The wave equations can be solved in the spectral domain using the double Fourier transform. In this paper, the Fourier transform pair is defined as where the ) , , ( z y x F function represents the fields E n (x, y, z) or H n (x, y, z). Application of the double Fourier transform to (7) and (8) yields a differential equation system whose general solution, in terms of the fields components, is given by where ) , ( are the amplitudes of the transformed field components, k x and k y are the spectral variables, γ n is the propagation constant in the n-th layer, ϑ = x, y or z, and Im( ) means the imaginary-part function. The τ variable, which defines the wave propagation direction, can assume values 1 or 2. Only the former value, representing propagation in the positive-z direction, occurs in the upper layer (free space). For the ground layer, on the other hand, τ equals 2, i.e., a wave propagating in the negative-z direction. For the confined layers, however, both values of τ will occur.
Interesting relations among the amplitudes of the transformed fields are derived by introducing the inverse Fourier transform of (11) and (12) in the Maxwell's curl equations, such that the amplitudes of the transversal components (x and y directions) are written as functions of the amplitude of the longitudinal ones (z direction). By enforcing the boundary conditions for the electromagnetic fields at each interface a set of 4N+4 equations with an equal number of unknowns is obtained. The analytical solution of this system leads to the spectral Green's functions. These functions, jointly with the transformed superficial density currents, allow the determination of the transformed fields at any point of the multilayer structure. The transformed electromagnetic field components are expressed by Green's functions in the n-th layer, which relate the ς (ς = x or y) components of the electric and magnetic fields to the transformed superficial density current j υς (k x , k y ) located at the interface d υ , with υ ∈ {g, 1, 2, …, N}. Note that d g = 0.

Moment Method -MoM
Once the Green's functions are derived, the next step is to set up integral equations constrained to the required boundary conditions. The integral equation is a statement of the boundary condition requiring that the total electric field tangential to the each of the perfectly conducting surfaces is zero [9]. That is, where S υ are the conducting surfaces, E s (x, y, d υ ) denotes the scattered field excited by the current on S υ , E i (x, y, d υ ) stands for the incident electric field and E r (x, y, d υ ) identifies the field that is reflected by the multilayer structure in the absence of patches. The E i (x, y, d υ ) and E r (x, y, d υ ) fields define the excitation mechanism of the structure, which in this analysis is due to an elliptically polarized plane wave at an arbitrary incidence angle. The currents induced on the conducting surfaces by these fields are unknown. To solve the electric field integral equation (16), with the unknown surface currents, the moment method is applied. This method is one the most popular numerical techniques used to analyze the radiation and scattering from complex structures. In the MoM, first the surface current is linearly expanded in a set of basis functions with unknown coefficients   where the left sides of (18) and (19) define the [V] matrix and the double integrals are related to the [Z] matrix.
The double integrations in equation (18) and (19) must be performed numerically, usually in a very inefficient and time-consuming way. In order to improve the computation efficiency some mathematical simplifications are employed. These simplifications include the evaluation of the even and the odd properties of the Green's functions, the change of the coordinate system (rectangular to polar) and the asymptotic extraction technique.
The far electromagnetic fields scattered by the multilayer structure are computed based on asymptotic expressions, which are derived from the stationary phase method [10]. The electric far field, using the stationary phase method, is given by in a spherical coordinate system, where the intrinsic impedance of free space is represented by η 0 , k xe = k 0 sinθ cosφ and k ye = k 0 sinθ sinφ are the stationary phase points, k 0 is the wave number of the excitation wave and r characterizes the distance between the receiving antenna and the target. Notice that from the knowledge of the electric far field it is possible to calculate the scattering matrix elements.

Four-Layer Structure
The approach described above is now applied to a structure, consisting of four layers (N = 2), as illustrated in Figure 2. A flat electric dipole of infinitesimal thickness is selected to represent the metallic patch. The dipole, printed on the interface z = d 2 (d 2 = ℓ 1 + ℓ 2 ) and oriented along the x direction, has dimensions 2a and 2b (a >> b) in the x and the y directions respectively. For this structure, the amplitudes of the transformed far field components (in the z direction) in free space are given by where α 0 = γ 0 d 2 , α 1 = γ 1 d 1 , and α = γ 2 (d 2 −d 1 ). Notice that only the amplitudes  0z (k x , k y ) and  0z (k x , k y ) are necessary to compute the electric far field, as shown in (20). In this particular case, the surface current along the y direction is neglected since the dipole width is considered to be very thin.
with the factors A 1 , A 2 , A 3 and  pm given, respectively, by )] Equation (31) is obtained from the modeling of the surface current density as the summation of piecewise-linear subdomain basis functions (rooftop functions) taking into account the edge condition. In this equation the sinc(.) is defined as sinc(α) = sin(α)/α, J 0 (.) stands for the zero-order Bessel function of the first kind and ∆x = 2a/(M+1).
From (27) to (31) it is noted that the first double integral of [Z pm ] is dependent on the operating frequency, whereas the second one is not. These are the well-known Sommerfeld integrals, which exhibit singularities in the form of branch points and poles; as such, their computation requires careful attention. The poles (generally complex) correspond to surface and leaky waves that can be excited in the layers. According to [11] the number of poles and their locations depends on the thickness of the layers, their relative electric permittivity and the wave number. For a multilayer structure and depending on the thickness of the layers, the Green's functions might present hundreds of poles, making their integration a formidable task. As an example, Figures 3 and 4 illustrate the real and the imaginary parts of the Green's function of a four-layer structure for two different operating frequencies: 1.25 GHz and 9.6 GHz, respectively. The electric parameters that characterize this structure are: ℓ 1 = ℓ 2 = 263.82 mm, ε r1 = ε r2 = 2.33, tanδ 1 = tanδ 2 = 1.2×10 -4 , ε rg = 5.0 and tanδ g = 2.0×10 -1 . Notice that the layers 1 and 2 have the same electric characteristics, and the variations of the real and the imaginary parts of the Green's function become more numerous and more abrupt as the frequency increases. These graphics illustrate that neither the analytical nor the numerical treatment of this kind of function is an easy task.
The first double integral of (27) could be computed by a singularity extraction method, however this technique requires the calculation of residues and Cauchy principal values at the singularity points. A major problem in computing the poles' contribution is finding the accurate location of all poles in a region. Since the number of poles and their locations are not known beforehand [9], the use of a deformed path to compute the integrations seems to be an efficient way to avoid this problem. Therefore, in this work a parabolic path was chosen to compute the integrations over the interval [0, B] and then proceeding along the Re(β)-axis from B to ∞, assuming there are no singularities in the latter interval. The deformed contour P, which is defined by P 1 and P 2 paths, is depicted in Figure 5, where P 1 represents the parabolic path and P 2 the path along the real axis. The advantage of contour P is avoiding numerical integration near the poles. Further, no knowledge of the number of poles and their locations is required. Generally, the values of parameters A and B are, respectively, about 0.1k 0 and  According to [12] the scattering of a flat electric dipole printed on the interface z = d 2 of a fourlayer structure can be accurately characterized by rooftop subdomain basis functions with twenty expansion modes and taking into account the edge condition. After that, the integrations are computed using the 96-point Gauss-Legendre quadrature rule, by truncating the upper limit of the β variable at 100,000 for the frequency-independent integral and at 50,000 for the frequency-dependent one, and by applying an additional subdivision technique, consisting of 10 subintervals, to compute the integral of the β variable. In this work, such MoM setup was used to compute the scattering of the printed dipole for simulation process purposes.

Polarimetric SAR Image Simulation
To exemplify the process of polarimetric SAR image simulation, the four-layer structure presented in Section 2.3 is considered. The parameters that can be varied in this structure are: the thickness of each confined layer, the dielectric characteristics of each layer (except for free space), and the size, orientation and location of the dipoles. Meanwhile, to define an image region having similar electromagnetic characteristic in its pixels, only dipole orientation is varied whereas the other parameters are constant.
Image generation begins with the definition of a regular rectangular grid of coordinates over the scene to be imaged. The grid size is determined by the dipole size and by SAR operational parameters, such as spatial resolution, pixel spacing and the extension of the area to be imaged. The coordinates of the rectangular grid determine the possible positions that any structure can occupy. Mutual interaction among the dipoles can be avoided by the definition of a guard band around each one. In Figure 6 is illustrated a rectangular grid showing a few selected dipole positions (blue points) and a zooming part depicting the guard band (gray circle) and the dipole orientation. In the simulation process it is assumed that there is no dominant scatter in the image. In addition, according to [13] six elementary scatters are sufficient, in practice, to the components in-phase (I) and quadrature (Q) of received backscattered signal have Gaussian distribution with zero mean and equal standard deviation σ. Considering that the multilayer structure is excited by (vertically and horizontally) polarized plane waves, the far electric field scattered by the structure is computed based on the electromagnetic model and the stationary phase method. In the spherical coordinate system, this field is given by (20), where it refers to the scattering element whose phase center is at the origin of the coordinate system. Each individual field scattered by a structure will contribute to the total return observed in each resolution cell. It is important to mention that due to the location of each scatter element phase center will appear an additional phase factor in its scattered field. This phase factor is given by exp{i k 0 d m cos ζ m }, where d m is the Euclidean distance to phase center of the m-th element and cos ζ m = cos ξ m sinθ cosφ + sin ξ m sinθ cosφ, with ξ m the angle measured counterclockwise from the x-axis to the distance line.
In order to form an image pixel the radar return is calculated by the vector summation of the individual fields weighted by a separable two-dimensional sinc(.) function (32). This function is introduced in the image generation process to represent the SAR spread point function, which takes into account the spatial correlation of the contiguous pixels in SAR image. The weight values are estimated based on a neighborhood of each resolution cell. The two-dimensional separable sinc(.) function is given by where x and r represent, respectively, the terrain azimuth and range coordinates, h 0 is a proportionality constant and the required spatial resolutions (in the azimuth and range directions) for the SAR image are expressed by δ a and δ r , respectively. Note that the xy-plane coincides with the azimuth-ground range plane.

Simulated Images
Multifrequency sets of single-look polarimetric SAR images have been generated in the L-, C-and Xbands, corresponding to 1.25, 5.3 and 9.6 GHz respectively. The acquisition geometry is particularized for a monostatic sensor flying at an altitude of 6,000 m (airborne platform altitude) and 35º grazing angle imaging a 290 m × 290 m area terrain. For this imaging geometry, the look angle between near-and farrange changes less than 1 o , that is a small variation around of the look angle at the center swath width. az rg The 3.0 m spatial resolution and 2.8 m pixel spacing were set in the range and the azimuth directions. In the simulation process the elementary scatterer was represented by a 50 mm × 1 mm electric dipole printed on the interface z = d 2 (see the structure in Figure 2). A 9 × 9 pixels square window was used to determine the neighborhood for the estimates of sinc(.) function weights.
The simulated images are based on a phantom image (an idealized cartoon model), which contains five different regions. The phantom image is depicted in Figure 7a. The analysis that follows is based on a set of twelve 10 × 10 pixels square samples for each image region, as shown in Figure 7b   Differentiation among the image regions is based on the local orientation of the dipoles and the electric characteristics of the layers. The dipole local orientation is relative to the azimuth-axis (az), and the coordinate system origin is located at the center of the image. The main characteristics of each region of the structure are summarized in Table 1, where TR stands for 'totally random'. It is assumed that the magnetic permeability of all layers is µ 0 and that the confined layers (layers 1 and 2 of the structure) have the same electric characteristics and thickness (ℓ i ), equivalent to 1.1 L-band wavelengths. As observed from Table 1, the main difference among regions A, B and C is the local orientation, whereas the relative electric permittivity and the loss tangent of the confined layers distinguish regions C, D and E.  [14] in the simulation process.
The first part of the SAR data analysis carried out is visual inspection (qualitative analysis). It shows that all images present a granular appearance typical of speckle noise. Mainly in the HV and VV channels, the boundary between some regions is unclear, due to the speckle noise effect. The mean backscatter plots, shown in Figure 11, where the error bars represent one standard deviation, reinforce the visual perception relative to the regions discrimination. From these graphics it can be seen that region A of the HH channel presents the lowest backscatter level in all bands, becoming the only region that can be clearly distinguishable from the others. In general, each region's backscatter depends on the frequency band. In the L-band, regions C and D of the HH channel have similar backscatter mean values, off by about 0.28 dB, which makes extremely difficult the separation between these regions. The largest difference between the mean backscatter levels among regions B, C and D of the HV channel is about 1.73 dB, for which the visual distinction between regions is already problematic. For the VV channel, regions A, D and E, as well as regions B and C, present the same separability issue. This fact can be confirmed by visual inspection of the HV and the VV channels, where distinction between regions becomes a hard task, except for region A of the HV channel.
In the C-band, the ability to discriminate the image regions is an issue for regions C and D in the HH channel, for the pairs of regions A-C and D-E in the HV channel, and for regions B and C in the VV channel. The largest difference between the mean backscatter levels for these regions reaches 1.94 dB. For the HH channel in X-band, all regions are visually distinguishable since the boundaries between any two contiguous regions can be clearly identified. This statement is not true however for regions B and C of the HV channel, as well as for regions A and E and for regions B and C of the VV channel, where the largest difference between the mean backscatter levels reaches 2.5 dB. Consequently, it can be stated that, for any variation in the backscatter mean levels that is less than 2.5 dB, the two corresponding regions will not be visually distinguishable.

Image Analysis
The purpose of this section is to carry out a quantitative analysis in the simulated images aiming at the validation of the simulation methodology. The analysis will be performed through statistical tests and feature extraction from SAR amplitude and polarimetric data. An application employing the simulated data is also shown using two classification procedures.

Amplitude Data
Within the SAR image processing community, the multiplicative model is widely used to describe statistically the data [4], [15]. From this modeling, it is well known that the amplitude SAR data (linear detection) from homogeneous areas obey a square root of gamma distribution, once for constant terrain backscatter (no texture) the observed variation is due to speckle noise. The square root of gamma is a two-parameter family of continuous distributions, having a scale parameter β and a shape parameter n.
This distribution can be written for every x > 0 by: where parameter n stands for the equivalent number of looks (ENL). The n and β parameters are estimated using an interactive procedure [21] based on the first and the second order moments of each selected sample. In order to test the hypothesis that the generated SAR data has a square root of gamma distribution, a χ 2 goodness-of-fit test was performed for all channels using only one sample of each image region. This sample has size of 1200 pixels since it was formed by grouping the training and test sample sets. The test was applied to amplitude data and their resulting p-values are listed in Table 2. The lowest p-value (11.51%) occurs for HV channel in C-band. Therefore, there is no evidence to reject the hypothesis that all regions for any channel in the three bands are homogeneous areas and exhibit a square root of gamma distribution. The p-values of χ 2 goodness-of-fit test for each channel in each band are shown in Table 2, where the largest p-values are typed in bold face. The data histogram (square symbols) and their corresponding fitted distribution (solid line) for these largest p-values are shown in Figures 12 to 14.  Another important quantity, commonly used within the SAR literature, is the equivalent number of looks, which can be estimated from the moments of the square root of gamma distribution. The estimated ENL was used as a quantitative measure for evaluating the one-look, using the individual samples of 100 pixels. The final estimates of ENL in each image region were computed as the average of those individual estimates, and are presented in Table 3. The estimated values are very close to one, as expected for a single look data. Under the linear detection and for single look data the ratio of the standard deviation and the expected value (the C v ) over homogeneous area is constant and equal to [(4-π)/π] 1/2 = 0.5227. This value can be obtained from the moments of the Rayleigh distribution, which is a particular case of the square root of gamma when ENL is equal to one. As a consequence, another simple way of evaluating the simulated data is to check whether the mean value (µ) and the standard deviation (σ) of the data holds the linear relationship σ = 0.5227×µ within homogeneous area. This analysis was performed by applying a simple linear regression model [16] to the estimated means and standard deviations for all samples in each channel that were properly fitted by a square root of gamma distribution. The fitted linear regression model states that Y = b 0 +b 1 ×X, where b 0 and b 1 are the estimated intercept and slope, respectively. In the case under analysis it is expected that the intercept should be zero and the slope should be equal to 0.5227. The adjusted linear model (Y = b 0 +b 1 ×X), where Y represents the standard deviations and X the mean values, and the expected linear model (Y = 0.5227×X) are depicted on Figures 15 to 17, respectively for L-, C-and X-bands in each channel. In these figures the straight line is represented by the expected linear model and the symbols follow the same legend adopted in Figure 11, i.e., the square, the triangle, the circle, the star and the diamond represent the regions A, B, C, D, and E, respectively. The obtained straight line for each channel is not shown in these figures because it is very close to the expected one. By observing the Figures 15 to 17 it can be noted that the estimated values for the intercept (b 0 ) are always around zero and, in general, the values of slope (b 1 ) are over estimated, but around the 0.5227 value.  Table 4. From the number of samples it is possible to note that few samples (approximately 3%) were not fitted well. Analyzing the p-values is observed that only one value lower than 5% is encountered for the slope test in HH channel of X-band, that is, for this case there is no evidence to accept the hypothesis of the slope to be equal to 0.5227.

Polarimetric Data
A major problem in analyzing polarimetric SAR data arises from the complexity of the scattering mechanisms that give rise to features in the different polarization parameters. A lot of work has been done for modeling polarimetric radar backscatter for various types of targets. In [17], for example, the relationship between the HH-VV polarization phase difference (PPD) at P-band and some forestry parameters, such as diameter at breast height, stand age, basal area and truck biomass for a pine plantation was analyzed. Similarly Durden et al. [18] reported the PPD returns show that the ground/tree interaction is important at P-band. Moreover, in [19], is presented a technique for unsupervised classification of scattering behavior by selecting the dominant scattering mechanism, based on PPD, where each pixel is classified as either an odd number of reflection (small PPD values), even number of reflection (large PPD values) or diffuse scattering. A PPD of 0 o means that the scattering mechanism is a single scattering, whereas a PPD of 180 o suggests that the scattering mechanism be a double-bounce scattering.
The PPD, given by (34), was used as quantitative measure to evaluate the generated images from the point of view of polarimetric data. The mean PPD values and its corresponding standard deviation (shown in Table 5) were computed for each image region in all bands, based on the twelve samples, as well as their histograms shown in Figures 18 to 20 for L-, C-and X-bands, respectively. In these figures the image regions are presented in alphabetic order. The histograms of PPD attribute appear to be symmetric for all bands and all type of image region. It was also performed a Kolmogorov-Smirnov goodness-of-fit test verifying that PPD values can be derived from a Gaussian distribution. The p-values obtained for these fits are shown in Table 5, all being greater than 30.56%. In Figures 18 to 20 a comparison with theoretical Gaussian distribution plotted over the histograms is also illustrated. Table 5. PPD mean and standard deviation and p-value for the goodness-of-fit test. The mean values of PPD attribute ranges from -17.91 o to 31.24 o , low values, which suggest that the major scattering mechanism is characterized by a single scattering (single bounce) for all image region. On the other hand, the standard deviations of PPD vary considerably with radar band and region type over the range 1.34 o to 82.98 o . The high values (consistently greater than 47.82 o ) for regions C, D, and E are indicators that these regions present structural heterogeneities. Due to the low values of standard deviations of the PPD the regions A and B can be seen as targets having uniform scattering properties (homogeneous structure).

Region L-band p-value (%) C-band p-value (%) X-band p-value (%)
The polarimetric images analysis follows by deriving some polarimetric features from the standard Cloude-Pottier eigenvalue/eigenvector target decomposition [20]. These features are deduced from the decomposition applied on the average coherency matrix. Under the reciprocity assumption framework and a monostatic measurement the coherency matrix can be expressed as a linear combination of the outer products of three eigenvectors. This means that the average coherency matrix can be decomposed into a sum of three independent scattering mechanisms, since each eigenvector corresponds to one scattering matrix. The eigenvalue/eigenvector decomposition of the coherency matrix into elementary mechanisms (i.e. single, double and volume scattering) is employed in order to identify the global mean scattering mechanism. From the eigenvalue/eigenvector can be defined the α -angle, which ranges from 0° to 90° and is used to represent physical scattering mechanism. Furthermore, eigenvalues can be combined to form the anisotropy (A) and the entropy (H) parameters, scalar quantities ranging from 0 to 1; the former is a measure of the degree of randomness the scattering process and the latter is a complementary measure to H, related to secondary scattering mechanism. Low entropy (H ≈ 0) indicates a single scattering mechanism (isotropic scattering) while high entropy (H ≈ 1) indicates a totally random mixture of scattering mechanisms with equal probability and hence a depolarizing target.
From the eigenvalue analysis it was observed that for all image regions their coherency matrix has only one nonzero eigenvalue (coherency matrix with rank 1). It leads to a zero entropy, which complies with a deterministic scattering process (or pure target), characterizing a single scattering matrix equivalent descriptor. It means that the region does not depolarize the incident wave, and in this case the anisotropy is zero also.
A pointwise estimation was employed to form an α image in all bands and their histograms are illustrated in Figure 21. This behavior is analogous in all regions, indicating that they have similar scattering mechanisms. Note that the greatest peak (maximum occurrence) is around 45°, suggesting that there exists a high percentage of volumetric scattering mechanism in the image regions and low contribution of the surface scattering (α = 0 o ) and the double bounce scattering ( α = 90 o ). In addition, the similarity among the region leads to a low discriminatory capability based on α . The histograms reinforce the remark that each region has a deterministic scattering process, since as stated in [20], an α value equal to 45° characterizes a dipole scatter mechanism. This result is expected because all regions were created with the same dipole having only different local orientation angles.

Data Classification
Digital classification is one of the most extensively used tools in remote sensing applications. Using this tool the discriminatory capability of polarimetric generated images is quantitatively evaluated. The classification procedure used is based on the Iterated Conditional Modes (ICM) algorithm [21][22][23], that is a supervised procedure.
The ICM method is a contextual procedure that, in order to classify every pixel, uses both the observed value in the corresponding co-ordinate and the classification of the surrounding sites. In order to incorporate the context within a statistical framework, a Markovian model is used for the classes. This model is known in the literature as Potts model [24]. The ICM algorithm consists of the iterative improvement of the classification of the co-ordinate s, using the information of its return and the classes of its neighboring sites. This improvement is obtained by maximizing the a posteriori distribution of classes given the observation and the surrounding classes, which is given by: where f ξ (z s ) is the density associated to class ξ, which has radiometric value z s on co-ordinate s, β is a real parameter that quantifies the influence of the neighboring class and it is estimated iteratively and ∂ s is the set of neighboring co-ordinates around site s. The iterative technique stops according to the number of co-ordinates whose classification changes from one iteration to the next [21]. The expression of L(ξ) can be reduced to the Maximum Likelihood (ML) classifier and to Mode Filter, when β = 0 and β → ∞, respectively. For more details of the algorithm, the reader is referred to [21,22,25].
In order to evaluate the discriminatory capability of the five image regions two classification approaches based on the ICM classifier were applied. The first one takes into account the bivariate distribution of the HH and VV intensities channels developed in [26,27]. The distribution of the pair of intensities arises from the multivariate complex Wishart distribution [28], which models the covariance matrix of the multilook polarimetric data, since it was assumed that the speckle obeys a multivariate complex Gaussian law [29] and the terrain backscatter has a constant distribution (no texture). The second classification approach is based on two polarimetric attributes derived from the HH and VV backscattering coefficient. The former is called by polarization discrimination ratio (PDR) was proposed in [6] and the later is denoted here as polarimetric description square root (PDS) being proposed in [30]. Both attributes aim for the retrieval of soil moisture from SAR data, and are, respectively, given by  Tables 6 to 11.
The classification results can be considered excellent in both approaches and for all bands, since all regions could be well distinguished from the others with little confusion among them. The results show the sensibility of the HH and VV channels and the attributes PDR and PDS to variations in the electric characteristic of the regions as well as variations in the elementary scatter orientation. This fact leads to further theoretical studies, using the simulated polarimetric images, aiming the extraction of geophysical parameters. In general, L-band presents the worst results of classification followed by C-and X-bands. In L-band the greatest confusion is found between C and D regions, followed by the confusion between C and E regions. On the other hand the confusion between B and C regions is the largest one when using Cand X-bands. These results suggest that the L-band might not be so adequate as C-and X-bands to conduct studies of the sensitive of electric characteristic of a target to microwave frequencies.  .κ σ κ Table 7. Confusion matrix (# of pixels) using the bivariate HH-VV distribution (C-Band). .κ σ κ Table 8. Confusion matrix (# of pixels) using the bivariate HH-VV distribution (X-Band). .κ σ κ A two-sided statistical z-test was performed to evaluate the equality between all pair of κ values.

Reference Data Classification
The tests of equality of two pairs of κ , at a significance level of 95%, demonstrated that all classifications obtained using C-and X-bands can be considered statistically equal. In this way, for the used parameter characteristic of each region and classification approaches, C-and X-bands presented the same performance. In contrast, in L-band the bivariate HH-VV classification is statistically equal only to the classification obtained by PDR-PDS attributes for L-band, showing that the information gathered by HH and VV channels is not significantly changed when these channels are combined in the PDR and PDS attributes for classification purposes.  .κ σ κ Table 10. Confusion matrix (# of pixels) using the Gaussian bivariate distribution for PDR-PDS (C-band). .κ σ κ Table 11. Confusion matrix (# of pixels) using the Gaussian bivariate distribution for PDR-PDS (X-band). .κ σ κ

Conclusions
This paper presented an electromagnetic way to simulated polarimetric SAR images, starting from Maxwell's equations. Images were simulated with five different regions, and their electromagnetic characteristics were used to distinguish them. The generated images were evaluated according to several measurements commonly employed in SAR data analysis. Firstly, the evaluation analysis consisted of statistical tests to amplitude data, showing that the data are adequately fitted by a square root of gamma distribution, which is the characteristic distribution of the multilook amplitude SAR data. Secondly, the equivalent number of looks was estimated, proving that the simulated data have only one look as simulated. It was checked by using a simple regression linear model whether the mean value (µ) and the standard deviation (σ) of the data exhibit the linear relationship σ = 0.5227×µ within homogeneous areas. From regression linear analysis it can be concluded that this relationship holds for all simulated images. It was also analyzed the polarimetric content of the data based on the HH-VV polarization phase difference (PPD) and the standard Cloude-Pottier eigenvalue/eigenvector target decomposition. The polarimetric evaluation revealed that the simulated data can be used as polarimetric SAR data, since the results are in accordance with those found in the literature. Finally, it was assessed the discriminatory capability of the image regions by applying two classification approaches based on the ICM classifier. The classification results showed that C-and X-bands have greater discriminatory power than L-band, for the parameters used to describe each region. Consequently, from the evaluation results can be affirmed that the simulation process is adequate and the simulated polarimetric data are in good agreement with those produced by a SAR sensor. Furthermore this simulation process can used to improve the understanding of SAR data properties in different situation, images can be generated to do theoretical as well as practical studies in several SAR subjects. For example, studies of the sensitivity of the relative electric permittivity and loss tangent of a target to certain microwave frequency can be carried out, which is an important topic in the retrieval of soil moisture content from SAR data.