InTech uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Earth and Planetary Sciences » "Land Applications of Radar Remote Sensing", book edited by Francesco Holecz, Paolo Pasquali, Nada Milisavljevic and Damien Closson, ISBN 978-953-51-1589-2, Published: June 11, 2014 under CC BY 3.0 license. © The Author(s).

Chapter 1

Adaptive Speckle Filtering in Radar Imagery

By Edmond Nezry
DOI: 10.5772/58593

Article top

Overview

TerraSAR-X 4-looks SAR images (HH and VV polarisations) acquired on August 10, 2007, over Oberpfaffenhofen, Germany (Credits: Astrium and InfoTerra; filtered images: Privateers NV, using ®SARScape).
Figure 1. TerraSAR-X 4-looks SAR images (HH and VV polarisations) acquired on August 10, 2007, over Oberpfaffenhofen, Germany (Credits: Astrium and InfoTerra; filtered images: Privateers NV, using ®SARScape).
RADARSAT-2 1-look SAR image (HH polarisation) acquired on April 30, 2009, over Capetown, South Africa (Credits: Canadian Space Agency and MDA; filtered images: Privateers NV, using ®SARScape).
Figure 2. RADARSAT-2 1-look SAR image (HH polarisation) acquired on April 30, 2009, over Capetown, South Africa (Credits: Canadian Space Agency and MDA; filtered images: Privateers NV, using ®SARScape).
Series of 18 (3x6) COSMO-SkyMed 1-look spotlight SAR images (HH and VV polarisations) acquired between February 2 and September 30, 2009, over Perito Moreno in Patagonia, Argentina (Credits: original images: Italian Space Agency; filtered images: Privateers NV, using ®SARScape).
Figure 3. Series of 18 (3x6) COSMO-SkyMed 1-look spotlight SAR images (HH and VV polarisations) acquired between February 2 and September 30, 2009, over Perito Moreno in Patagonia, Argentina (Credits: original images: Italian Space Agency; filtered images: Privateers NV, using ®SARScape).
ALOS-PALSAR 6-looks polarimetric SAR imagery acquired on June 30, 2006 over Bavaria, Germany (Credits: JAXA and MITI; filtered images: Privateers NV, using ®SARScape).
Figure 4. ALOS-PALSAR 6-looks polarimetric SAR imagery acquired on June 30, 2006 over Bavaria, Germany (Credits: JAXA and MITI; filtered images: Privateers NV, using ®SARScape).

Adaptive Speckle Filtering in Radar Imagery

Edmond Nezry1

1. Introduction

Historically, foundations in the domain of radar speckle properties have been laid down from the 1940’s to the 1980’s. Decisive theoretical advances were made by teams led by Professor Fawwaz Ulaby at the University of the Michigan (USA), by Professor Christopher Oliver at the Defence Research Agency in Great Malvern (UK), and by Professor Keith Raney at the Canadian Centre of Remote Sensing (Canada) since the 1970’s. Then, the domain of speckle filtering in SAR images matured in the period 1976-2000, mostly under the impulsion of Dr. Jong Sen Lee of the Naval Research Center, Washington D.C. (USA). Since 1986, the team led by Dr. Armand Lopès at the Centre d'Etude Spatiale des Rayonnements in Toulouse (France), has carried out and then inspired the development of the most efficient speckle filters existing today. Since 2000, with speckle filters having reached a satisfactory level of performance, no significant advances have been made. Nevertheless, in this period, the use of speckle filters in a wide range of applications using SAR imagery has become generalized.

A radar wave can be considered, with a good approximation, as plane, coherent and monochromatic. It is emitted by an antenna towards a target. The target backscatters partially the radar wave in the direction of a receiving antenna. In the vast majority of spaceborne Synthetic Aperture Radars (SAR), a single antenna assumes the two functions of emission and reception (monostatic radar).

The complete radar measurement is the combination of the horizontally (H) and vertically (V) linearly polarised radar waves, at emission and at reception after backscattering by the observed target. After signal calibration, this measurement, affected by noise, enables to restitute for each resolution cell a polarimetric backscattering matrix S:

S = [Spq] = (SHHSHVSVHSVV)
(1)

whose coefficients Spq=|Spq|.exp(jpq) are the complex backscattered field amplitudes for p-transmit and q-received wave polarisations. They relate the backscattered wave vector Ed to the incident wave vector Ei: Ed = exp(j.k.r).S. Ei. This complete radar measurement is said "fully polarimetric", often abbreviated in "polarimetric".

Let consider the interaction of the radar wave with an extended random surface, i.e. a surface containing a sufficiently great number of scatterers per resolution cell, with no preponderant scatterer with respect to the others. For whatever configuration pq of polarisation, the signal Spq received from such a surface by the antenna becomes, after quadratic detection in intensity Ipq:

Ipq = |Spq|2 = Spq.Spq*(where*denotes complex conjugation)
(2)

This detected signal intensity I is proportional in average to the radar "backscattering coefficient" σ°. The backscattering coefficient σ°=4.π.|Spq|2 is the average radar cross-section per surface unit [1]. σ°, expressed in m2/m2, is a dimensionless physical quantity. It is a physical property of the sensed surface, which depends principally on its roughness, its dielectric properties, its geometry, and the arrangement of its individual scatterers.

Carrying the radiometric information with regard to the sensed target, σ° is a function of the frequency of the radar wave, of its angle of incidence upon the target, and of the configuration of polarisation. In terms of physical meaning, the radar backscattering coefficient is analogous to the bidirectional reflectance in the domain of optical wavelengths: σ° # 4 cos2θ, where is the incidence angle of illumination on the sensed target.

Nevertheless, detected radar images look visually very noisy, exhibiting a very characteristic salt-and-pepper appearance with strong tonal variations from a pixel to the next. Indeed, since radar imaging systems are time-coherent, radar measurements over random rough surfaces are corrupted by "speckle" noise due to the random modulation of waves reflected by the elementary scatterers randomly located in the resolution cell. Then, coherent summation of the phases of elementary scatterers within the resolution cell results in a random phase of the complex pixel value.

This speckle "noise" makes both photo-interpretation and the estimation of σ° extremely difficult. Actually, speckle is a physical phenomenon, which is inherent to all coherent imaging systems (radar, lidar, sonar, echography). In most remote sensing applications using radar/SAR imagery, speckle is generally considered a very strong noise that must be energically filtered to obtain an image on which classic and proven information extraction techniques could be further applied, in particular the techniques used for optical imagery acquired in the visible and near-infrared part of the electromagnetic spectrum.

Therefore, speckle filtering and radar reflectivity restoration are among the main fields of interest in radar images processing for remote sensing. Speckle filtering is a pre-processing aiming at the restoration of σ° value in the first place. This pre-processing must account for both the particular properties of the speckle, and those of extended imaged targets (often called "clutter"). It must also account for the radar imaging system that has sensed the target and for the processor that has generated the images. For stationary targets of infinite size, speckle filtering is equivalent to a simple smoothing using a moving processing (averaging) window. An ideal filter must nevertheless avoid image degradation through excessive smoothing of the signal. To this end, it must respect structural image information (road and water networks, etc.), and the contours of radiometric entities. In addition, it must also respect the particular texture of the clutter, in forested or urban environments for example. Last, it must also identify the strong local radiometric variations due to the presence of strong scatterers (often artificial in nature) from those due to spatially rapid speckle fluctuations.

Therefore, an ideal speckle filter must satisfy to the following specifications:

  1. Preserve accurately the local mean value of the radar reflectivity (i.e. the quantity actually measured by the radar, which is proportional to σ°) to enable, for example, the comparison of radar reflectivities in the framework of a multitemporal analysis of radar acquisition series.

  2. Smooth as much as possible homogeneous image areas and therefore reduce the speckle to increase the Equivalent Number of Looks (ENL) of the radar image (cf. § 2.1.4). The minimum ENL depends on the desired radiometric accuracy. For example, a 1 dB accuracy with 90% confidence level (i.e. less than 25% variation of the radar intensity) requires an ENL value around 230, and a 2 dB accuracy with 90% confidence level (i.e. less than 60% variation of the radar intensity) needs an ENL value around 40.

  3. Preserve texture as much as possible where it exists in the image (forests, non-homogeneous fields, etc.) to avoid confusions among radiometrically similar areas exhibiting different texture. Therefore, a speckle filter must be able to discriminate heterogeneity effects due to texture from those due to speckle.

  4. Both preserve and denoise image structures (contours, lines) as well as the quasi-deterministic responses due to corner reflector effects within strongly textured areas such as urban environments. Indeed, the energy of artificial radar reflectors responses must be preserved to enable radiometric calibration, in particular when calibration targets are dispersed in the radar image.

  5. Minimise, and whenever possible prevent loss in useful spatial resolution during the speckle filtering process.

2. Statistical properties of speckle and texture in radar images

In this section, the statistical properties of speckle in images produced by coherent imaging systems such as imaging radars, lidars or sonars, are exposed. Since a good speckle filter must restore the texture of the scene imaged by the radar, the statistical properties of texture in radar images are examined as well. This analysis intentionally restrains to the first order statistical properties, since only these are generally used by the estimation techniques involved in speckle reduction methods. Explicit use of second order statistical properties of both the speckle and the imaged scene in the filtering process is adressed in Section 4.

2.1. Speckle in SAR images

First, let consider a natural target with radar backscattering coefficient σ° which remains constant from a resolution cell to the next. Such a target, said stationary, "homogeneous", or "textureless", may be found in a medium-resolution radar image within cultivated fields or forest cut areas for example. Spatial variations of the radar signal backscattered by such a target are therefore only due to speckle and the statistical properties of the radar signal are those of the speckle, which depends only on the radar system and the image production system.

2.1.1. First order statistics for "1-look" complex radar images

The radar imaging system is linear, spatially invariant, and can be characterised at each image pixel (x,r), where x is the azimuth and r is the radial distance with respect to the flight path of the sensor’s carrier, by a complex impulse response h(x,r), which defines the resolution cells. The 3dB widths of the impulse response, often used to quantify the spatial resolutions δ(x) in arimuth and δ(r) in radial distance.

If the number of individual scatterers within a resolution cell is large enough and none of them has absolute predominance in scattering the radar wave [2] [3], speckle can be modelled in output of the radar processor by a circular complex gaussian random process uc(x,r) at the corresponding image pixel of coordinates (x,r). The complex amplitude (complex radar signal, in a complex radar image) in output of the radar processor can be expressed as Ac(x,r)=a(x,r)+j.b(x,r), where a(x,r) is proportional to the field backscattered by the surface element located in (x,r). In the literature, the terms a(x,r) and b(x,r) are generally called i(x,r) ("in-phase" term) and q(x,r) ("in-quadrature" or "out-of-phase" term), respectively.

These complex data in output of the radar processor are called a "1-look complex" image or equivalently, a "Single-Look-Complex" (SLC) image.

For a homogeneous area where A(x,r)=A0 (constant), i.e. in the presence of speckle only, the complex amplitude has the form [4]:

Ac(x,r) = [A(x,r).uc(x,r)]*h(x,r) = A0. uc(x,r)*h(x,r) = A0. V(x,r)
(3)

where A=A0 is the amplitude of the backscattered wave and * is the convolution operator. The random process uc(x,r) represents the circular white gaussian (by application of the central limit theorem over a large number of individual scatterers within the resolution cell) complex process responsible for the speckle. Then, V(x,r) is a correlated complex gaussian process characterising the "gaussian fully developed speckle", which behaves as multiplicative "noise" with a very high spatial frequency dependent on h(x,r).

Owing to the particular character of the speckle phenomenon in coherent imagery, which is the case of radar imagery, information extraction from a radar image results in a statistical estimation problem. It is therefore mandatory to have an as complete as possible statistical speckle model. The statistical model of the fully developed speckle has proven perfectly adapted to SAR images, at the spatial resolutions / wavelength frequencies combinations actually used in radar remote sensing. Note that this speckle model emphasises the importance of wave phase information under the condition that the phase has not been lost during signal detection.

Goodman hypotheses [2] [5] enable to calculate speckle statistics before (i.e. in a complex radar image) and after radar signal quadratic detection in amplitude A or in intensity I. Taking this hypotheses into account, one obtains:

E[i]=E[q]=0           and           Var[i]=Var[q]=σi2 =σq2 =σ2
(4)

where E[.] denotes the mean expectation and Var[.] denotes the variance. Therefore, the mean backscattered intensity R is expressed as:

R=E[I]=E[A2]=E[i2+q2]=2σ2
(5)

R=E[I] is proportional to the radar backscattering coefficient σ°. R can be estimated by R≈<I>, where <.> denotes the averaging operator applied in a neighbourhood of the pixel under consideration. Then, from R, it is possible to retrieve the σ° value through the calibration parameters of the radar image. Problems related to radar image calibration are discussed in detail in [6] [7] [8]. R=E[I] can be estimated locally in a radar image through the spatial averaging <I> of the intensities of a number N of image pixels in a spatial neighbourhood of the image pixel of interest [9]. This estimate is the unbiased Maximum Likelihood estimate with minimal variance [6].

Besides, the in-phase and the out-of-phase components i and q of the complex radar signal are decorrelated, E[i.q]=0, and therefore independent of each other. It results of the above considerations that the probability density functions (pdf) of i and q are Gaussian distributions:

P(i) =1/2π.σ2. exp(i2/2σ2)     and    P(q) =1/2π.σ2. exp(q2/2σ2)
(6)

with a phase of the complex radar signal ϕ=Arctg(q/i), which is uniformly distributed in the interval [0, 2π].

2.1.2. First order statistics for "1-look" detected radar images

Since i and q are independent variables, their joint pdf is P(i,q)=P(i).P(q). The pdf Pu(I) of the speckle in intensity is obtained by a simple change of variable I=a2+b2 in Equation (6). Then, Pu(I) results being an exponential pdf:

Pu(I) = (1/E[I]) . exp(I/E[I])        for      I>0
(7)

It is important to note that, since E[I] is proportional to the radar reflectivity R, Pu(I) is the pdf of I conditional to the value of R, therefore Pu(I)=Pu(I/R).

The two first first-order moments of this pdf can be expressed as a function of its standard-deviation as follows:

mean   E[I] = 2σ2    and     standard-deviation   σI = 2σ2
(8)

To characterise the strength of speckle in radar images, it is convenient to consider the normalised coefficient of variation of the intensity, CI, i.e. the standard-deviation of the radar signal intensity I over its mean value:

E[CI]  =  σI/E[I] =  1
(9)

The value of the coefficient of variation of speckle only, sometimes called "contrast index", is a constant for every type of radar image product. Equation (9) means that, in a radar image, the dispersion (variance) of radiometry increases as the mean signal backscattered by the target increases. This justifies in part the qualification of "multiplicative noise" given to the speckle.

Clearly, with a signal-to-noise ratio of 1, the radiometric degradation due to speckle makes very difficult the discrimination of two homogeneous targets with very different radar backscattering coefficients. As an example, a theoretical computation [10] demonstrates that two textureless target classes of homogeneous radar reflectivities (i.e. exhibiting only speckle) radiometrically separated by 2.5 dB present a probability of confusion of 40% in a 1-look SAR image.

Therefore, a first radiometric enhancement is needed to achieve a reduction of the coefficient of variation of the speckle over homogeneous areas. It corresponds to an enhancement of the signal-to-noise ratio and to a preliminary speckle "noise" reduction.

2.1.3. First order statistics in "multilook" images

A first method of speckle reduction consists in averaging incoherently M independent samples of the intensity I (or "looks") obtained from the same target:

I = (1/M) . k=1MIk with each of the Ik distributed according to Equation (7)
(10)

The goal of this method is to reduce speckle enough to make radar image photo-interpretation possible. Indeed, experience has shown that M values of the order of 3 or 4 enable a photo-interpreter to use a SAR image [1]. Such values have been adopted for most spaceborne SARs (ERS, Almaz, JERS-1, Radarsat, Envisat, ALOS: 3-looks, Seasat, SIR-B: 4-looks), including recent ones (TerraSAR-X, Cosmo-Skymed, Sentinel-1, etc.).

This operation is realised by splitting the Doppler bandwidth of the backscattered signal into M sections. This is equivalent to divide the length of the synthetic antenna into M sub-antennas. Independent processing of each section results in an independent image called a "look". Nevertheless, this operation results in a degradation by a factor M of the spatial resolution in azimuth of each look [1].

Over the same target, the mean intensity resulting from this operation remains the same mean intensity of each of the individual looks. If the M individual looks were independent, the standard-deviation is divided by M. Thus, the coefficient of variation of the speckle measured over a homogeneous area of the intensity multilook image becomes:

CI (homogeneous area)  =  Cu =  1/M
(11)

It is important to note that multilook radar image formation is at the expense of the spatial resolution in azimuth. In practice, the value of M does never exceed a few units (less than 16 for airborne SARs, and in general only 3 or 4 for spaceborne SARs). This remains insufficient to improve satisfactorily the signal-to-noise ratio of a radar image. Indeed, in our example of two target classes with homogeneous radar reflectivities radiometrically separated by 2.5 dB, the probability of confusion of 40% in a 1-look SAR image reduces only to 33% in a 4-look image.

If the individual looks are uncorrelated with each other, the pdf of the speckle, which is the sum of M independent exponential distributions, becomes a χ2 distribution with 2M degrees of freedom:

Pu(I) = Pu(I/R) = (M/E[I])M / (M1)! . exp(M.I/E[I]) . I(M1)
(12)

2.1.4. "Equivalent Number of Looks" of a SAR image

If the M sections of Doppler bandwidth used to produce the individual looks overlap, the averaged samples, and therefore the individual looks, are correlated. The coefficient of variation of the speckle Cu, measured over a perfectly homogeneous image area will therefore be always superior to what could have been expected from Equation (11). In this situation, speckle strength results as if the number of independent looks were equal to:

L = 1/Cu2   <  M
(13)

The L value, which is generally a non-integer value, is called the "Equivalent Number of Looks", or "ENL", of the multilook radar image.

Hence, the pdf of the speckle in intensity can be approximated, by extension of Equation (12), for whatever value of L, by a Gamma distribution with parameters E[I] and L:

Pu(I) = Pu(I/R) =(L/E[I]) LΓ(L) . exp(L.I/E[I]) . I(L1)
(14)

which is rigorously equivalent to Equation (12) when L=M (L integer) for M independent looks.

All speckle models for multilook images (M correlated or uncorrelated averaged looks, i.e. L equivalent looks, with M L) use this approximated distribution. Note that the pdf of the speckle as formulated in Equation (14) is slightly inexact if the M looks are correlated, as it is the case in general. Therefore, the pdf of the speckle for the average of M correlated looks is close to, but is not exactly, a Gamma pdf [11].

Let us consider that the M looks (I1,I2,...,IM) are correlated to each other with correlation coefficients ρij, and that the Ii are distributed according to the same exponential marginal pdf (1-look, cf. Equation (7)) with same parameter R=E[I]=E[Ii], for whatever i. This last condition can be actually fulfilled by re-weighting the individual looks to make them contain an identical mean energy E[I].

For the intensity image resulting of the averaging of the M looks, I=(∑Ii)/M, the set of looks correlations ρij is taken into account to compute the ENL value L (L<M), which is related to the coefficient of variation of the speckle in intensity, Cu, by:

Cu2 =  1/L = [ 1 + (1/M) .  i=1Mj=1i1σij2] / M
(15)

The exact pdf of the average intensity I of M correlated looks, is extremely difficult to establish mathematically in a closed-form (for 2-looks, cf. [11]), and, when established, to manipulate. Nevertheless, Equation (14) is a very close and therefore satisfactory approximation when using the appropriate ENL value L calculated using Equation (15).

2.1.5. The logarithmic (homomorphic) transformation

Since the radar backscattering coefficient σ° is generally expressed in dB in physics, and with the goal to make speckle an additive noise that would be independent on the level of the radar signal, some authors [12] [13] [14] chose to use a neperian logarithmic (homomorphic) transformation:

D =  Log (I)
(16)

Arsenault and April [15] [17] [19] demonstrated that, after this transformation, the pdf of the speckle for a multilook radar image, Equation (14), becomes a Fisher-Tippett distribution:

Pu(D) =  LL/ Γ(L) .  exp[L.(DD0)] . exp{L . exp[(DD0)]}
(17)

with D0=Log(E[I])

The first order statistical moments of this distribution are as follows [15]:

1) Mean <D>:

<D> = D0  Γ'(L)/Γ(L) + Log(L)
(18)

Equation (18) shows that the logarithmic transformation causes a signal distortion, increasing with decreasing number L of independent speckle samples taken into account, and with decreasing R (therefore, σ°) values. There is therefore a tendency to systematically underestimate the value of σ° [15] [16] [18].

2) Variance σD2 :

σD2 = ( π2/6) . [ γ+ Γ'(L)/Γ(L)]2  +  2  k=1L2 [ (1/(Lk). j=1Lk1(1/j) ]  
(19)

Lopès [20] has shown that for N independent L-look radar data samples (NL independent speckle samples), the standard-deviation of the samples averaged after logarithmic transformation, and then retransformed into intensity, is always significantly larger than the standard-deviation σI of the same samples averaged in intensity. Thus, the logarithmic transformation degrades both the measurement accuracy of the backscattering coefficient σ° and its local range of fluctuations due to local scene texture.

2.2. Texture in SAR images

Texture concerns the spatial distribution of grey levels. It is important in the analysis of SAR images for a wide range of applications (e.g. [21] [22] [23] [24] etc.). In most of these applications, not only the radiometric, but also the textural information must be retrievable after adaptive speckle filtering of the SAR images.

2.2.1. Texture of the imaged scene

As seen above, within a homogeneous area image of a detected radar image, one can consider the observed speckle I as originating from a real non-gaussian random process u(x,r), with unit-mean <u>=1, and proportional to the radar signal I (multiplicative noise). A simplification of Equation (3) enables to write, for every pixel (x,r) located within a homogeneous image area (i.e.R(x,y)=E[R]=R, (x,y)):

I(x,r) = u(x,r) . E[I(x,r)]= u(x,r) . R
(20)

In most remote sensing applications, it is reasonable to consider the imaged scene as an arrangement of discrete objects superimposed to a background with mean reflectivity E[R]. The imaged scene is composed of classes of non-homogeneous (R(x,y)E[R]) objects characterised by the statistical properties and parameters of the variable R [25]:

If Tj is the random variable which represents the spatial fluctuations of the reflectance (let consider in the following its analog in radar images, the denoised/speckle-filtered radar reflectivity) R(x,r) around its mean expectation E[R(x,r)], in the terrain class j to which the pixel at coordinates (x,r) belongs, the following results are obtained:

  1. R(x,r) is proportional to the σ° value of the resolution cell containing pixel (x,r);

  2. E[R(x,r)] is proportional to σ° in a neighbourhood of pixel (x,r) belonging to class j;

  3. E[Tj]=1 for every class j;

R(x,r) = E[R(x,r)] . Tj
(21)

The spatial structure of actual imaged scenes, called "texture", induces a mesurable spatial structure (after the properties of  Tj) in the images of these scenes. This is the case, for both optical and radar images. Indeed, a simple visual interpretation of a radar image reveals radiometric spatial variations at a longer scale than those due to speckle only. These variations originate from the spatial fluctuations of the radar backscattering coefficient σ° within a given terrain class, and are affected by the presence of speckle.

To characterise a given class j, one must be able to restore  Tj, that means to separate the respective effects of speckle and texture in the spatial variations of the radar intensity signal. To this end, the multiplicative fully-developed speckle model exposed above will be used. It has been shown [26] [3] that, even when the scene is not stationary, speckle still remains multiplicative, eventually correlated, and independent on the imaged scene as long as Goodman’s conditions [2] [5] are realised.

2.2.2. First order statistics of texture in SAR images

2.2.2.1. Spatially uncorrelated speckle

Considering the scene texture model of [25] in Equation (21), it is possible [27] [28] to model the radar image intensity at pixel of coordinates (x,r) by generalising the multiplicative speckle model for a wide-sense stationary scene (Equation (20)) as follows:

I(x,r) = R(x,r) . u(x,r) = (E[R(x,r)].Tj) . u(x,r)
(22)

This relationship is valid as long as the spatial variations of R(x,r) happen at length scales longer than the size of the resolution cell of the radar imaging system.

The variance of the intensity I within a given target class is computed using:

σI2 = E[I2]E2[I] = E2[R].(E[T2.u2]E2[T.u])
(23)

Since the fluctuations of σ°, and therefore those of R, are independent on the speckle, one gets [27] [28]:

(σI/E[I])2 = (σR/E[R])2 .(σu/E[u])2 + (σR/E[R])2 + (σu/E[u])
(24)

and:

E[I] = E[R.u] = E[R].E[u] = E[R]
(25)

Therefore, E[R] is locally estimated in the radar image by averaging pixel intensities in a neighbourhood of pixel (x,r):

Introducing the coefficients of variation (CI=σI/<I> for the radar intensity, CR=σR/<R> for the imaged scene, Cu=σu/<u> for the speckle), Equations (24) and (25) lead to the important result [27] [28]:

CR2  =  (CI2  Cu2) / (1 + Cu2)
(27)

which characterises scene texture in terms of heterogeneity, with Cu2 = CuI2 = 1/L for a homogeneous/textureless area of a L-look intensity radar image. Indeed, Equation (27) shows that the more heterogeneous is the scene, the easier its texture can be restored.

For a better description of the scene, one must use the pdfs of its diverse classes (distribution of the random variable R or σ°). Their knowledge may result, either from a priori knowledge, or from direct estimation in the image. If the general form of the pdf is known, its parameters can be estimated from the image data. In all cases, if one can establish the pdf of R, PR(R), for a given terrain class in the scene, the unconditional pdf of the corresponding radar image intensity I is:

PI(I) =   0Pu(I) . PR(R) . dR   
(28)

where Pu(I) is the pdf of the speckle for multilook images (cf. Equation (14)). Equation (28) shows, as does also Equation (27), that the respective contributions of scene texture and speckle micro-texture can be separated in the radar image.

2.2.2.2. Spatially correlated speckle

To produce a SAR image, speckle samples at the output of the SAR processor are correlated through the SAR impulse response to obtain a sampling rate (pixel size) of about half the spatial resolution of the radar sensor, thus avoiding aliasing effects. Using the multiplicative noise model and Equations (24) and (25) [27] [28], the correlation coefficient ρI(Δx,Δy) between two pixels of intensities I1 and I2 separated by (Δxy) in the SAR image:

ρI(Δx,Δy) = (<I1(x,y)<I1>> . <I2(x,y)<I2>>) / (σI1 .σI2)
(29)

can be related to the scene correlation coefficient ρR(Δx,Δy) and to the speckle correlation coefficient ρu(Δx,Δy) by [10]:

<CI>2. ρI(Δx,Δy) = <CR>2.<Cu>2.ρR(Δx,Δy).ρu(Δx,Δy) + <Cu>2.ρu(Δx,Δy) + <CR>2.ρR(Δx,Δy)
(30)

Therefore, if ρu(Δx,Δy) >0 for Δx>0 and/or Δy>0, which is the case in correctly sampled SAR images, the local statistics estimated in the neighbourhood of a given pixel will have to be corrected of the effect of speckle spatial correlation properties (cf. Section 4).

2.2.3. Local statistics computation and consequences of the speckle model

Using the simplified multiplicative noise model for the speckle u (cf. Equation (22)), the first-order non-stationary statistics of the scene, <R> and σR2, can be deduced locally (cf. Equations (26) and (27)) from those of the radar image intensity, <I> et σI2.

In practice, CI is estimated in the radar image, and its locally estimated value can be inferior to Cu. Since σR2 and CR2 are positive quantities, finding cases where CI2<Cu2 must be attributed to the limited size of the neighbourhood of the pixel under processing over which <I> and σI2 are estimated [29]. Independently of the scene model PR(R), the multiplicative speckle model fixes therefore an inferior threshold CImin=Cu to the possible values of the coefficient of variation CI.

Below the CImin value, speckle filters based on the use of local statistics, it means all adaptive speckle filters, are no longer valid. In this case, the image area under processing must be considered homogeneous and textureless.

2.2.4. The coefficient of variation as a heterogeneity indicator

The coefficient of variation CI is a heterogeneity measure particularly well appropriate to the case of radar images and for isotropic textures. Nevertheless, if the source of heterogeneity presents a particuliar orientation (contour, line...), additional detectors able to identify or retrieve this orientation are needed.

CI increases with scene heterogeneity. Depending on the heterogeneity of the area under processing, it enables to discriminate among the following situations:

  1. homogeneous areas with neither texture, nor structures.

  2. areas with low or moderate heterogeneity where the presence of structures is improbable.

  3. strongly heterogeneous areas where one must search for the possible presence of structures (contour, line, strong scatterers).

3. Adaptive speckle filters for single-channel detected radar images

This section is dedicated to a theoretical analysis of the most common and more efficient speckle filtering techniques developed for Synthetic Aperture Radar (SAR) images, and of their corresponding estimation techniques. These filters use variants of the statistical speckle model exposed in the preceding Section 1. They also use diverse statistical estimators to restore the radar reflectivity (in the sense of the CEOS’s (Committee for Earth Observation Satellites) "radar brightness"). These statistical estimators are, either Minimum Mean Square Error estimators (e.g. Lee et al. filter, Kuan et al. filter, etc...), or autoregressive estimators (e.g. Frost et al. filter), or Bayesian estimators (e.g. Gamma-Gamma and DE-Gamma MAP filters). These estimators are discussed, and their behaviours are analysed.

3.1. Wiener Method: the Frost et al. filter (1982) [30]

3.1.1. Theoretical development

The SAR image model considered by [30] is as follows:

I(t) = [R(t) . u(t)] * h(t)
(31)

where R(x,y) is a wide-sense stationary random process decribing the radar reflectivity of the observed scene at pixel of coordinates t=(x,y) located in a homogeneous area of the radar image. u(t) is the multiplicative noise due to speckle, modelled by a real white stationary random process with a Gamma pdf (cf. Equation (14)). h(t) is the impulse response function of the system. This model is a simplification of the fully developed speckle model exposed above (cf. Equation (3); [31] [32].

To estimate the radar reflectivity R(t) of the noise-free image, Frost et al. apply a linear filter whose impulse response m(t) minimises the mean quadratic error (MMSE estimator):

 ε= E[{ R(t) I(t)*h(t) }2]
(32)

The MMSE least-squares solution is valid for homogeneous areas for which R(t) can be represented by a stationary random process. To deduce m(t), [30] assume that the transfer function H(f) is constant over a certain bandwidth. This leads to an uncorrelated multiplicative speckle model:

I(t) = R(t) . u(t)
(33)

where u(t) is a pseudo-white noise, independent of the signal, also used by other authors [33] [34]. The multiplicative speckle model is deduced from Equation (3) assuming that the bandwidth of the signal R(t) is much narrower that that of the linear filter h(t).

The impulse response of the filter is calculated adopting an autoregressive process model for R(t) with an exponential autocorrelation function, a classic hypothesis for this family of models [35] [36]. Later studies showed that the choice of an exponential autocorrelation function is appropriate to SAR scenes [37].

These hypotheses enable to define an optimal MMSE filter (Wiener filter) with impulse response m(t) :

R^=  t  neighborhood of N pixelst=(x,y)m(t) . I(t)
(34)

with:

m(t)  =  K1 . exp[ K . CI2 . d(t) ]
(35)

where K is the filter parameter, and d(t) is the distance between the pixel located in t and the pixel under processing. K1 is a normalisation constant granting the preservation of the mean value of the radar reflectivity:

K1 =  (1/N).  t neighborhood of N pixelst=(x,y)m(t)
(36)

The final Frost et al. filter (Equation (36)) is not exactly a Wiener filter. It is designed using Yaglom’s method [38]. In this method, contrarily to that of Wiener, one searches for the frequential characteristic, not for the impulse transfer function susceptible not to exist. This function is chosen empirically taking into account the hypotheses it must fit.

3.1.2. Behaviour of the Frost et al. filter

The behaviour of the filter depends on the values of the locally observed coefficient of variation, i.e. on the local heterogeneity:

  • Extremely homogeneous areas, for which CI = 0,     then:     R^ = <I>

  • Very strong scatterers (extreme heterogeneity), where CI , then:      R^ = I

  • Between these two situations, in textured natural areas, when CI value increases, neighbour pixels located far away from the pixel under processing are given less weight and the neighbourhood on which the weighted mean value is actually estimated narrows.

Considering the requirements for an ideal speckle filter, the filter should:

  • Restore the mean intensity value over homogeneous areas, where CICu. A consequence is that K might be close to 0, unless the value of Cu tends towards 0, which is obviously not possible for usual radar images.

  • Preserve the observed intensity when the value of CI is very high but obviously always finite. For high CI values, the only pixel to weigh might be the pixel under processing, and K might therefore tend towards infinity since CI has always a finite value.

Since it is impossible to fulfil simultaneously these two conditions (cf. [39] [40], one must make a difficult choice between, either an efficient speckle reduction, or the preservation of image structures and texture.

3.2. Locally adaptive linear MMSE estimators: Lee and Kuan filters

3.2.1. Kuan et al. filter (1985) [41] [32]

The radar image I(t) is modelled as a function of the radar reflectivity image R(t) to be restored and of a white uncorrelated noise n(t) with 0 mean and dependent on R(t). This is an additive noise model different of the multiplicative speckle model used in Equation (33):

I(t) = R(t) + n(t)
(37)

Kuan et al. [41] introduce a scene model R(t) with non-stationary mean and non-stationary variance (NMNV model), where the non-stationary mean describes the general structure of the image, whereas the non-stationary variance characterises texture and the presence of local structures. With this model, the linear filter minimising the quadratic mean error (LMMSE) has the general form [32]:

R^LMMSE =  E[R]  +  CRI. CI1. ( I <I>)
(38)

where CRI and CI are the non-stationary spatial covariance matrices of the image. These matrices are diagonal if, as it is assumed in the model, n(t) is a white noise and the image model is NMNV. Then, Equation (38) takes a scalar form [31]:

R^ =    E[R] + ( I <I>) . σR2 / (σR2 +σn2 )
(39)

The local statistics of the scene R are deduced from those of the intensity I. The Kuan et al. filter is obtained by replacing locally:

E[R]   by  <I>      and      σR2  by  [σR2 <I2>.σn2 ] / [1+σn2]  
(40)

Since the variance of the additive noise σn2 is numerically equal to the coefficient of variation of the speckle Cu in the multiplicative model, Equation (40) is, in practice, exactly equivalent to Equations (26) and (27) established for the multiplicative speckle model.

It is noteworthy that, since the image model does not assume independence of the noise n and the signal R, the filter is theoretically able to take into account the effects of an eventual dependence between speckle and the radar reflectivity, when speckle is no longer fully developed as in the case of strong scatterers.

3.2.2. Lee filter (1980) [33] [42] [43] [44] [45]

Lee uses the unit-mean uncorrelated multiplicative speckle modele (cf. Equation (33)). A linear approximation is done by developing Equation (33) (I as a function of u for pixel located in t) in a first-order Taylor series with respect to u(t):

I(t) = R(t). E[u(t)] +  E[R(t)] .[u(t)E[u(t)] ]        with    E[u(t)] = 1
(41)

This approximation enables to transform Equation (33) into a weighted sum of the signal and of a noise independant on the signal. The linear MMSE estimator, and Equations (26) and (27) that are consequences of the multiplicative speckle model, enable to establish the Lee filter [33], historically the first speckle filter designed to be adaptive to local radar image statistics:

R^ = <I> + ( I <I>). (σI2 σu2)/σI2
(42)

It is noteworthy that, assuming the independence of noise and signal in the model used by the Kuan et al. filter, this filter become identical to the Lee filter [41].

3.2.3. Behaviour of the Kuan et al. and Lee filters

The linear MMSE speckle filters of Lee and of Kuan et al., based on the use of the local statistics of the observed intensity can be written under the same general form:

R^(t) =  I(t). W(t) + <I(t)>. [1W(t)]
(43)

with:

W(t) = (1[Cu2/CI2(t)])/(1+Cu2)         for the Kuan filter
(44)

and

W(t) = 1  [Cu2 / CI2(t)]    for the Lee filter
(45)

Both methods perform a linearly weighted average of the local mean of the intensity and of the observed pixel intensity. In both cases, weights depend on the ratio of the coefficients of variation of the observed intensity and of the noise [7] [40]. As for the Frost et al. filter, the locally observed heterogeneity governs the behaviour of these filters:

  • Homogeneous areas where CI=Cu, therefore R^=<I>

  • Very stong scatterers (extreme heterogeneity) where CI →∞, therefore R^=I

  • Between these two extreme situations, the weight of the observed value of the pixel intensity increases with the heterogeneity of its neighbourhood.

Nevertheless, note that the weight on the estimated mean intensity remains significant for high values of the coefficient of variation. Thus, responses of small impulse targets are cut-off by these filters and image structures suffer some amount of smoothing. To correct this drawback, Lopès et al. [40] have proposed a very efficient enhancement to the Lee and Kuan et al. filters.

Inversely, when CI <Cu, some weighting on the observed intensity reappears, resulting in an amplification of the noise [13]. In this case, the local mean of its neighbour pixels must be assigned to the pixel under processing as filtered value.

The linear MMSE filters differ from the Wiener filters by the fact that the A Priori mean and variance, <R> and σR2, are estimated using the local statistics of the original image I, and not implicitly estimated from an assumed autocorrelation model (cf. [30], Equation (35)). The goal of this approach is to obtain a better adaptivity to the local properties of the scene R, which verifies in pratice.

It is remarkable that the often used MMSE estimation ([33] [43] [44] [45] and [41] [32], among others) is nothing else than the mean of the conditional A Posteriori pdf P(R/I) [46]. Note that the complete evaluation of the terms enabling to obtain the A Posteriori pdf P(R/I) requires a description of PI(I). This last pdf, which depends on the non-linear transformation leading from R to I is, in general a K-distribution. The solution adopted by Lee and Kuan et al. consists in a forced linearization of the problem assuming a Gaussian PI(I) [33] [42] [44], which is unjustified for low (< 3) ENL values.

As an effect, this linearization restrains the validity of the linear MMSE estimator to situations where the noise level is not too high (multilook images with high values of L), or where the local heterogeneity CI is not too high. These considerations justify one of the modifications proposed to the Lee and Kuan et al. filters par Lopès et al. [40]: an upper heterogeneity threshold CImax, above which the observed intensity of the pixel under processing is preserved.

3.3. Bayesian Maximum A Posteriori (MAP) speckle filters

3.3.1. Originality of the MAP approach in the case of SAR images

The Bayesian MAP approach [47] [48] consists in characterising the imaged scene and the speckle "noise" by their statistical description, using their associated pdfs.

In the Bayesian perspective, the theory of probabilities is extended to the logics of probabilistic inference. Probabilities are seen as a relationship between a formal hypothesis and a possible conclusion. This relationship corresponds to a certain degree of rational credibility and is limited only by the extreme relationships of certitude and impossibility [49]. The classic deductive logics considering only these extreme relationships (à la Sherlock Holmes: "When all impossible hypotheses have been excluded, what remains, even if unlikely, corresponds to reality") is nothing but a particular case of this general development [50].

As a general consequence, the theory of probabilities cannot base on the sole concepts of classic logics (frequentist as in the MMSE or in Wiener’s approaches). In particular, the relationship of probability cannot be defined in terms of certitude, since certitude is viewed as a particular case of probability. The frequentist definition of probabilities based on relations of certitude related to the knowledge of a number of parameters (<R> for example, cf. § 3.3.5 below) is therefore no longer sufficient [50].

In this context, probabilities are used to describe stochastically an incomplete information on a global phenomenon (here, the radar reflectivity and the superimposed speckle), rather than to describe only the noise randomness that corrupts its comprehension. Probability relationships are viewed as being conditional to the context. This way, the pdf of the speckle Pu(I) (cf. § 2.1.4, Equation (14)) is therefore formulated under the form P(I/R). Taking into account a reasonable A Priori statistical model for the radar reflectivity P(R), one must also take into account – which is new, with respect to the methods precedently exposed – the probability of R, given the information obtained through observation (the intensites I measured by the radar), formulated under the form of an A Posteriori pdf P(R/I). This approach has the great advantage to enable the characterisation of speckled radar image formation while easing the description of non-linear effects.

The least error-cost inference mechanism leading from the observed intensity I to the A Posteriori deduced radar reflectivity R through Bayes’s theorem [47] [48], allows a rigorous combination of the A Priori knowledge P(R) and of the new knowledge provided by the observation I:

P(R/I)  =   P(I/R) . P(R) / P(I)
(46)

P(R/I) depends on the pdf of R introduced as A Priori information about the scene to restore. Thus, the estimate is influenced by the A Priori statistical knowledge about R or, by default, the hypotheses made about R. The noisiest are the radar data, the less the A Priori information would contribute to estimate R [46].

In theory, the MAP method enables to avoid direct estimation of the mean of the conditional A Posteriori pdf, which is necessary to the MMSE estimation. This feature is of great interest in the resolution of non-linear problems where the evaluation of the conditional mean is, either difficult, or impossible [51].

The Maximum A Posteriori (MAP) estimation of R corresponds to the maximum value of P(R/I), i.e. to the mode of the A Posteriori pdf P(R/I). The mode and the mean of P(R/I) coincide if P(R/I) is symmetrical, in particular in the Gaussian case. Since the available knowledge is:

  • P(I/R), the conditional pdf of speckle, which is a Gamma distribution (cf. Equation (14)),

  • P(R), the A Priori unconditional pdf of the imaged scene, which is an assumed statistical distribution (not necessarily Gaussian),

there is no particular justification to prefer the MMSE estimation by minimisation of the mean quadratic error implicitly assuming a symmetrical P(R/I), rather than to use this more complete available knowledge through the Bayesian estimation technique, which is less subject to arbitrary assumptions [46].

It is of interest to notice that the "Maximum Likelihood" estimation is nothing more than a particular case of Maximum A Posteriori, in the situation where the conditional A Posteriori and A Priori pdf’s are equal: P(R/I)=P(I/R), with P(R)=P(I). Note that the Maximum Likelihood estimation of the scene R is the radar image I itself, which is of interest in the case of strong scatterers whose radar response is deterministic.

Once the first-order statistical speckle model P(I/R) is either known, or reasonably approximated, (cf. Equation (14)), the form of the MAP estimate depends principally on the form of P(R).

The NMNV scene model [32] enables to solve locally, in an analytic manner, the estimation of R. The parameters of the pdf P(R) are locally estimated in a neighbourhood of the pixel under processing. This way, the complicated recursive forms of mathematical resolution as exposed by Hunt [52], Kuan et al. [53] and Geman & Geman [54] are avoided. This approach is particularly adapted to the high spatial resolution of airborne and modern spaceborne SARs. Indeed, it allows a good preservation of scene texture in an adaptive way by taking into account the local fluctuations of σ° in the A Priori model, while being at the same time computationally efficient.

3.3.2. Computation of the MAP estimate

Since the logarithm function is a monotonically increasing function, Bayes’s formula ([47] [48]; Equation (46)) can be rewritten as:

Log[P(R/I)] = Log[P(I/R)]  Log[P(I)] + Log[P(R)]
(47)

which gives the local MAP estimation of R when Log[P(R/I)]  is maximum, i.e. when its first-order derivative with respect to R is locally equal to 0:

/R Log[P(I/R)]+/R Log[P(R)]=/R Log[P(R/I)]+/R Log[P(I)]  
(48)

with:

/ R Log[P(I)]=0 because P(I) does not depend on the local fluctuations of R.

P(R/I) reaches its maximum when:

/R Log[P(R/I)]= 0   for   R=R^MAP,   where R^MAPis the MAP estimation of R
(49)

Then, the general equation of the MAP speckle filter becomes, locally:

/R Log[P(I/R)]+/R Log[P(R)]= 0    for   R =  R^MAP
(50)

The first term of Equation (50), the Maximum Likelihood term, accounts for the effects of the whole imaging system on the radar image and describes the detected radar intensity once the speckle statistical model is known. The second term, the Maximum A Priori term, represents the prior statistical knowledge with regard to the imaged scene.

In the Bayesian approach, probabilities are used to describe incomplete information rather than randomness. As Equation (50) shows, in the Bayesian inference process, induction is influenced by the prior expectations allowed by the prior knowledge of P(R) [46]. In addition, the non-linear system and scene effects are implicitly taken into account by the restoration process. Therefore, MAP speckle filtering can be considered as a controlled restoration of R, where A Priori knowledge controls the inference restoration process and allows an accurate estimation of the radar backscattering coefficient σ°.

The pdf of the speckle in intensity for a L-looks radar image, P(I/R), is a Gamma distribution with parameters R and L (cf. Equation (14)). The Maximum Likelihood term for a L-looks radar image is then equal to:

/R Log[P(I/R)]  =  L . ( I/R2  1/R )
(51)

The Maximum A Priori term must be calculated according to the scene model chosen as A Priori knowledge, hypothesis, or belief.

3.3.3. Gaussian distributed imaged scene: The Gaussian-Gamma MAP filter (1987)

The hypothesis of Gaussian-distributed scene has been adopted as a natural hypothesis by a large number of authors who had worked, either on images from optical sensors, or on images from passive/active microwave sensors. These authors are comforted in this hypothesis by both the force of habit and by the mathematical ease in manipulating a Gaussian distribution.

Kuan et al. [32] have developed a MAP filter (Gaussian-Gamma MAP filter, with a Gaussian-distributed scene and Gamma-distributed speckle) for radar images under this hypothesis. Though, the hypothesis of a Gaussian-distributed scene, although widely spread in the literature, is inappropriate. Indeed, this hypothesis assumes implicitly the theoretical possibility of negative σ° values-which has no physical sense-in the extreme case of a large variance with a low mean value, therefore needing a regularisation of the filter behaviour in such a situation.

Therefore, one must preferably take into account a positive pdf as a realistic scene model. For reasons that are both experimental and theoretical, in the case of natural extended targets as it is most often the case dealt with in remote sensing, a Gamma distributed scene model is better appropriate.

3.3.4. Gamma distributed imaged scene: The Gamma-Gamma MAP filter (1990)

Natural textures observed, either by coherent radar imagery, or by incoherent optical imagery are due to a common contribution corresponding to the variability in spatial distribution of the objects within the scene. Even if the interaction mechanisms between the electromagnetic wave and the observed medium are very different in either case, the natural arrangement of the scene makes the second-order statistics very similar in either kind of imagery [55] [56]. At the scale of a large number of resolution cells, the pdfs of the cross-section variables corresponding to either mechanism belong to the same family of distributions, at least for the high radar frequencies in bands Ku, X, and C for which wave penetration into natural media is limited. This point is more arguable for radar bands L and P.

In a wide range of radar backscattering situations, the Gamma distribution is experimentally the one that best fits, not only the distribution of the radar backscattering coefficient [21] [57], but also the distribution of radiometries observed in incoherent optical images [55]. This scene model has been successfully used also for radar images of the sea [58] [59] [60].

The local pdf of a scene statistically described by a Gamma distribution, has the form:

P(R) =  (α/E[R])α/Γ(α) . exp(α.R/E[R]) . Rα1
(52)

with α= 1/CR2. The parameter is called the "heterogeneity coefficient".

Note that assuming a Gamma-distributed R, and by performing the integration in Equation (28), the pdf P(I)=P(I/R) of the intensity is a K-distribution [57]. Introduced in 1976 [58] by British researchers of the RSRE (later DRA) to describe the non-Gaussian properties of waves backscattered by objects within a radar resolution cell, the K-distribution has been theoretically recognised as the pdf of the intensity I backscattered by a rough non-stationary surface [61] such as most natural scenes observed by a radar.

Nevertheless, in an illustration of the so-called "Cromwell’s rule" [62], even hard A Priori conviction that the scene presents a Gamma-distributed texture must not be insensitive to counter-evidence. Therefore, the complete MAP filter is a set of three filters adapted to diverse situations locally encountered in a radar image: the application of this or that filter is decided depending on the degree of heterogeneity of the image part under processing, that is on the locally estimated value of the coefficient of variation CI. This may eventually imply the determination of thresholds on CI calculated as a function of a user-defined probability of false alarm with respect to the local presence of texture or of strong scatterers (cf. [63] [64] [10]). Note that these considerations apply to all MAP speckle filters, and can be extended to all other adaptive speckle filters.

3.3.4.1. Textured areas

Assuming that the pdf of the scene is a Gamma distribution, the Maximum A Priori term is locally equal to: / R Log[P(R)]=(α-1)/R-α/<R>. Once E[R] is estimated locally by <I> in an ergodic stationary neighbourhood of the pixel under processing in the radar image, the equation of the Gamma-Gamma MAP filter ([65] [10] [63] [64]) is:

α.R2 + (1+Lα).<I>.R  L.I.<I>  =  0
(53)

This second-degree equation admits only one real positive solution R in the interval ranging between the mean intensity <I>=<R> and its observed value I. Therefore, the Gamma-Gamma MAP estimate of the radar reflectivity of the pixel under processing, is:

R^=<I>.(αL1) +  <I>2.(αL1)2 + 4α.L.I.<I> 2α
(54)

The integration of a heterogeneity/texture detector based on the coefficient de variation and of specific detectors (ratio-of-amplitudes – RoA; [66] [7]) for contours, linear structures and strong scatterers in the filtering process is described by the general algorithm presented in [63] [64] [10]. The integration of texture and structure detectors using second-order statistics (autocorrelation functions) of both the speckle and the radar reflectivity of the scene [67] is presented in Section 4. In all cases, mage areas identified as textured are filtered using Equation (54).

3.3.4.2. Homogeneous areas

In the particular case of a perfectly homogeneous (textureless) scene, with CICu, the radar reflectivity R is a constant (R=E[R]), and can be statistically represented by a Dirac distribution:

P(R)  =  δ(R)
(55)

This distribution is the limit of Gamma distributions when tends towards+∞. In such a case, the MAP estimate is equal to the local mean intensity:

This case is taken into account when the local statistics calculated in the neighbourhood of the processed pixel show a nearly perfect homogeneity of the scene.

3.3.4.3. Strong scatterers and impulsional targets

The other extreme case regards strong scatterers, when speckle is no longer fully developed. The considerations that led to consider the Gamma pdf as an A Priori model for a textured natural scene are no longer valid: we no longer have any A Priori information about the scene. In this situation, the information content of every grey level of the image being A Priori the same, the pdf of the imaged scene can be represented by an uniform distribution:

P(R) =  1/(Rmax  Rmin)
(57)

with undetermined extreme values Rmin and Rmax, which we may consider equal to respectively 0 and+∞. Thus, the Maximum A Priori term of the MAP equation becomes / R[Log(P(R))]=0, and the MAP estimate is:

If the resolution cell contains only one isolated strong scatterer, the response (value I of the pixel) of this scatterer is deterministic and must therefore be preserved. This situation leads to the same conclusion and is therefore treated similarly.

3.3.5. Scene pdf estimated from the data: The Gamma/Distribution-Entropy MAP filter (1998)

SAR images of dense tropical forest, urban areas, or very strong and rapidly varying topography often show very strong or mixed textures. This is also the case of high-and very-high spatial resolution SAR images. In these situations, it may be hazardous to make an assumption about the probability density function of the radar reflectivity.

Indeed, the MAP technique does not account for any uncertainty in the value of the parameters of the A Priori pdf chosen as a Gamma distribution once it has been locally estimated on a given image area. Hence, in the presence of mixed (forests with underlying structures, for example) or rapidly varying (strongly textured area located on strong slopes, or very-high spatial resolution radar images, for example) texture, the MAP estimator will underestimate the variance of the predictive distribution. Indeed, this predictive distribution can hardly take into account the fact that it results of a compound of a mix of different distributions.

In this context, the A Priori knowledge with regard to the observed scene can hardly be an analytical first order statistical model, chosen on the base of prior scene knowledge. However, to retrieve local statistical scene knowledge directly from SAR image data, Datcu & Walessa [68] [69] proposed to introduce the local entropy of the radar reflectivity, S(R), as a measure of local textural disorder. This concept originates from the theory of information. S(R) is estimated on a neighbourhood (Npix pixels) of the pixel under processing:

S(R) = k=1Npix[Rk. log(Rk)]   
(59)

Because the radar reflectivities Rk are non-negative and exp[S(R)]/Z is normalised, and since S(R) is a measure of the spread/dispersion of the radar reflectivities of the scene, characterising its heterogeneity, Equation (59) can be treated as a pdf whose entropy is S(R) [70]:

P(R) = (1/Z. exp[S(R)] =  (1/Z). exp(k=1Npix[Rk. log(Rk)] )
(60)

For a single detected SAR image, the conditional pdf of the speckle can be, as long as speckle is fully developed, modelled as a Gamma distribution:

P(I/R) = (L/R)L/Γ(L) . exp(L.I/R) . IL1
(61)

Incorporating these scene and speckle models, the Gamma/Distribution-Entropy MAP (Gamma-DE MAP) speckle filter for single-channel detected SAR images is the solution of the following equation [71]:

L.I  L.R  R2.k=1Npix[log(Rk)1/Ln(10)]  = 0 
(62)

The radar reflectivites Rk in the neighbourhood of the pixel under processing are pre-estimated by a first speckle filtering pass.

Note that the local DE MAP estimation of R is contrained by the value of the entropy S(R) retrieved from image data. Since this fixes upper/lower bounds to the local entropy S(R), the DE MAP speckle filter combines both the Bayesian MAP and the Maximum/Minimum Entropy estimation techniques.

The DE-MAP filters adapt to a much larger range of textures than the other MAP filters ([10] [11] [72] [73]) developed under the assumption of K-distributed SAR intensity (i.e. Gamma-distributed scene). In particular, these filters are of particular interest in the case of very high-resolution SAR images and strongly textured scenes.

Compared to those of the other MAP filters, performances in terms of speckle reduction are identical. However, texture restoration and structures or point targets preservation, identical for moderate textures, are superior in strongly textured areas. These filters have proven a remarkable efficiency in operational remote sensing (cf. [74]).

3.3.6. Behaviour of the MAP filters

The local MAP estimate is the mode of the local A Posteriori pdf of the radar reflectivity R, under the hypothesis made for P(R). The MAP estimation is therefore probabilistic: it corresponds to the most probable value of radar reflectivity. This justifies to chose the local mean observed intensity as filtered value within areas identified as perfectly homogeneous (cf. § 2.4.4; [64] [10). Excellent speckle reduction and preservation of the mean value of the radar reflectivity are therefore granted over such areas.

In the presence of texture, the estimator takes into account the non-linear effects in getting from R to I via the imaging system. As for the Lee and Kuan et al. filters, the local A Priori mean and variance, <R> et σR2, are estimated using the first-order non-stationary local statistics of the original image I (cf. Equations (26) and (27)), thus allowing adaptivity to the local isotropic properties of the scene R. At the scale of a terrain parcel, the MAP filter must theoretically restore the whole pdf of the radar reflectivity. This includes its mean E[R] and its variance (equivalent to parameter α), as well as the spatial relationships between pixels if one introduces speckle and scene second-order statistics in the filtering process (cf. [10] [75] [67]).

Nevertheless, in the presence of structures, the NMNV model may be used abusively on an inappropriate neighbourhood of the pixel under processing and therefore the local statistics <R> and σR are no longer those of an ergodic process. Thus, in order to use the filter in the conditions required by the NMNV model, i.e. in the statistical situation it has been designed for, one must introduce, as one might also do for the Lee and Kuan et al. filters, structure detectors. These detectors contribute to select a neighbourhood of the pixel under processing composed of pixels pertaining to the same thematic class. They may be, either geometrical improvements based on the use of ratio (RoA) detectors ([7] [66] [63] [64] [10]), or detectors based on the use of the second-order statistics of both the speckle and the scene ([67]; cf. Section 4). The desired result is a better preservation of image texture, structures, and responses of strong scatterers in the speckle filtered radar image.

4. Using speckle and scene spatial correlation to preserve structural information in SAR images

The assumptions usually made by adaptive speckle filters ([76] [30] [32] [10]) with regard to the first order statistics of the speckle are:

  • the multiplicative speckle model;

  • the pdf of the speckle which reflects only incompletely the spatial correlation of the speckle through the ENL.

  • The assumptions regarding the imaged scene are:

  • Local wide-sense stationarity and ergodicity;

  • A formal Non-stationary Mean Non-stationary Variance (NMNV) model accounting for the spatial variation of the first order statistical properties [32];

  • For the Maximum A Posteriori (MAP) filter, the general form of the pdf of the scene radar cross-section [10].

The main consequences of these assumptions are that:

  • The formal NMNV model justifies local treatment (processing window) using the local statistics;

  • Assuming spatially uncorrelated speckle leads to the design of scalar (single-point) filters;

  • Local adaptivity to tonal properties (mean intensity <I> and scene radar reflectivity R) is achieved by controlling the filtering process through the local coefficient of variation CR of the radar reflectivity R of the scene, through the local coefficient of variation CI of the original SAR image intensity I [27] [28].

The major drawbacks of these assumptions are that:

  • The coefficients of variation are statistically [10] sensitive to texture and speckle strength [40], but do not provide direct information on spatial correlation properties and texture directionality. Locally estimated, they can also be biased if speckle samples are not independent (cf. § 2.2.2), as it is often the case in SAR images. Indeed, speckle can have strong spatial correlation properties, a situation the filters are not designed to deal with;

  • Compliance with the ergodicity and stationarity hypothesis requires a preliminary identification of the structural elements of the scene to fully exploit the NMNV model and to correctly estimate the local mean reflectivity around the pixel under processing;

In this Section, local second order properties, describing spatial relationships between pixels are introduced into single-point speckle adaptive filtering processes, in order to account for the effects of speckle and scene spatial correlations. To this end, texture measures originating from the local autocorrelation functions (ACF) are used to refine the evaluation of the non-stationary first order local statistics as well as to detect the structural elements of the scene [75].

4.1. Effects of the spatial correlation on the speckle

In practice, the usual single-point filters do preserve texture, only due to the spatial variation of the local first order statistics, using the NMNV formal model [32]. Fairly good preservation of textural properties and structural elements (edges, lines, strong scatterers) can be achieved by associating constant false alarm rate structure detectors such as the directional Ratio-Of-Amplitudes (RoA) detectors; ([7] [66] [10] [63] [64]) to the speckle filtering process. The combination of detection and speckle filtering allows to preserve scene structures, and to enhance scene texture estimation on a shape adaptive neighborhood of the pixel under consideration. When the conditions for which they have been developed (especially spatially uncorrelated or low-correlated speckle) are fulfilled, the best single-point adaptive speckle filters and their structure retaining associated processes based on RoA detectors retain enough scene texture to allow its use as an additional discriminator ([23] [24] [63] [64]).

However, the performances of the usual single-point filters, even when refined with associated RoA structure detectors, degrade when the actual spatial correlation of speckle samples becomes significant [77], i.e. generally when SAR images are far too much oversampled [31]. In fact, the uncorrelated speckle model approximation is seldom justified, according to sampling requirements of the radar signal and to the SAR system and processor which have produced the image. In multilook SAR images, spatial correlations are introduced by a series of weighting functions: complex coherent weighting related to data acquisition (antenna pattern, pre-sum filter and Doppler modulation in azimuth), coding of transmitted pulses in range, and looks formation (useful Doppler bandwidth selection, extraction window in azimuth). In practice, they are related to sampling and scale requirements in the design of SAR image products, up to the extreme case of signal oversampling [31].

As an example of the effects of speckle spatial correlation, the case of structure detection using the RoA detectors [10] implemented in the Gamma-Gamma MAP filter is illustrated in [67]. Detection is successfully performed when the correlation of speckle is low between pixels. However, the performance of RoA detectors substantially degrades when the spatial correlation of speckle samples becomes significant. It has been shown theoretically in [67] that spatial speckle correlation results in an increasing detection of non-existing structures, i.e. in the increase of the probability of false alarm (Pfa). The expected 1% Pfa for uncorrelated speckle raises to about 25% if the correlation of adjacent pixels is 0.7 in both range and azimuth (cf. ERS: about 0.55 in range and 0.65 in azimuth). In addition, the spatial speckle correlation also decreases the probability of detection. The combination of these two effects causes artefacts in the filtered image (visual “crumbled paper” effect), making photo-interpretation uneasy and reducing substantially the performance of further automatic image processing like segmentation [78] and classification [23] [24].

Therefore, second order statistical properties must also be considered, including both scene texture and resolution/sampling related properties, to achieve a more complete restoration of the radar reflectivity.

4.2. A possible solution: Spatial vector filtering

A possible solution is the implementation of vector (multiple-points) filters ([32] [79]) where the spatial covariance matrices of the speckle and of the scene are taken into account. The development of a filtering method using second order statistics, i.e. an estimation of the complete ACF of the scene through the speckled image ACF, results in a vectorial equation giving rise to a multiple-points filter, such as the LMMSE vector filter developed theoretically by Kuan et al. [32]. A practical implementation of multiple-point filters has been first suggested by Lopès et al. [10], and then detailed by Lopès & Séry [79]. Such filters have good spatial memory, but they result in a complicated implementation and very heavy computations.

4.3. Single-point speckle filtering using spatial second order statistics

To avoid the mathematical complexity and the heavy computational burden of multiple-point filters, an alternative solution is to introduce an appropriate description of both speckle correlation properties and spatial relations between resolution cells into a single-point filter.

Second order statistics have explicitly been used in the past, following the scheme proposed by Woods & Biemond [80] implemented later by Quelle & Boucher [81] in the adaptive Frost speckle (single-point) filter [30]. The filter, which belongs to the family of Wiener filters, is established using Yaglom's method [38], where a frequency characteristic is determined, instead of an impulse transfer function, which actually exists only within wide-sense stationary areas. This frequency characteristic is determined in the speckled image at a distance of one pixel only in range and azimuth, thus mainly limiting the description of correlation properties to those of the speckle.

Considering the NMNV model [32], a much better solution is to introduce locally estimated second order statistics to refine the computation of the local NMNV first order statistics. In this way, speckle correlation properties can be taken into account for filtering. In addition, the local mean radar reflectivity E[R] would be estimated, taking better account of the textural properties of the imaged scene (natural variability of the radar reflectivity, including spatial relations between resolution cells). On a given textured class of the scene, when scene correlation length is smaller than the processing window size, the mean radar reflectivity E[R] will then be modulated by the scene ACF, i.e. by a function of the correlation coefficients of the radar reflectivity in all possible directions. When the correlation length of the imaged scene is longer than the processing window size, the estimate for the non-stationary mean radar reflectivity tends towards to the classical Maximum Likelihood estimate E[R]=<I>.

4.4. Local ACF’s and texture fields

Spatial relations between pixels are well described by the intensity ACF, defined by a set of correlation coefficients z), where Δz=(Δa, Δr), the normalised ACF of an intensity SAR image, {ρI(Δz)}, is a composition of the underlying scene ACF, {ρR(Δz)}, convoluted with an overlap function depending on the point spread function (PSF).

The local estimates for the normalised intensity ACF, {ρI(Δz)}, are, either deduced from the spatial autocovariance, which is computed on the domain of interest D (N pairs of pixels separated from each other by Δz) as follows:

Co^vI(Δz) = (1/N.D[I(z+Δz)<I>] .[I(z)<I>]
(63)

or, equivalently, directly computed within the domain of interest D as follows:

ρ^I(Δz) =D([I(z+Δz)<I>].[I(z)<I>])D[I(z)<I> ]2
(64)

These estimates have the attractive property that their mean square error is generally smaller than that of other estimators; considerations on the estimation accuracy can be found in Rignot & Kwok [82].

The contribution of scene texture must be separated from that of the speckle for all displacements Δz: the scene ACF is deduced from the intensity image ACF, using the following transformation [28]:

ρ^R(Δz) = [ 1 + ρ^I(Δz).CI2] [1 + |Gc(Δz)|2/L] 
(65)

where L is the equivalent number of independent looks, CI is the local coefficient of variation of the image intensity, and Gc(Δz) is the normalised ACF of the individual 1-look complex amplitudes (i.e. the ACF of the SAR imaging system).

Gc(Δz) depends only on the SAR complex PSF [83]. It is realistic to adopt an exponential form for the ACF of the speckle for the separated looks detected in intensity:

|Gc(Δz)| = exp [(Δa.dx/rx+Δr.dy/ry)]
(66)

where dx, dy are the pixel dimensions, and rx, ry are the spatial resolutions in azimuth and range directions. However, it is preferable to use directly the true ACF of the speckle in the actual multilook SAR image, if it is available, or if its estimation is possible.

The computation of a local non-stationary estimate of E(R) is performed by a simple convolution in the domain of interest around the pixel under consideration of the normalised ACF of the textured scene {ρR(Δz)} with the intensity:

E^(R) = D[ρ^R(Δz).I(z+Δz)]D[ρ^R(Δz)]  
(67)

for all ρR(Δz) within D, and as long as the following condition applies:

[ρ^R(z) /ρ^R(z=(0,0))]   >  1/e
(68)

This means that, along a direction z, the pixels whose correlation ρ^R(z+Δz) to the pixel under process is less than 1/e are ignored for the computation of the first order NMNV local statistics through Equation (67). This way, the algorithm achieves powerful speckle reduction in homogeneous areas, preservation of the textural properties (heterogeneity and directionality) wherever they exist, and correct structure detection.

The implementation of this operator, which estimates the local non-stationary E[R] is inspired from the notion of "texture fields" introduced by Faugeras & Pratt [84]. It is important to notice the following properties of the operator:

  1. in homogeneous areas, Equation (67) acts as a smoothing operator: E[R]<I>, which is the Maximum Likelihood estimate of E[R] for a perfectly homogeneous neighbourhood. (The local energy scaling factor in the denominator of Equation (67) ensures the preservation of the absolute level of the radar reflectivity).

  2. for a strong scatterer, it acts as a Laplacian: E[R]I(z).

  3. if a perfect step edge between two homogeneous areas is present in the neighbourhood (correlations=1,0,-1), it acts as a gradient and automatically stops weighting further pixels. This indicates structure detection capability;

  4. in all other cases, i.e. for a textured scene, the E[R] estimate is a weighted average of the neighbourhood, modulated by the scene ACF in all possible directions.

This concept of "texture fields", based on the analysis of second order statistical properties, generalises image processing concepts such as filtering and detection, which are usually considered as distinct. It provides the possibility to perform simultaneously structure detection and image pre-filtering. Structural elements such as edges are detected (value of the estimate of the correlation coefficient becoming less than 1/e) when the distance from the pixel under processing (central pixel of the processing window) to the edge is reached, going in the direction of the edge. Therefore, the position of the edge is known as soon as the edge enters the “field of view” of the processing window.

The local mean radar reflectivity estimation takes now into account the textural properties of the imaged scene (spatial variability of R between resolution cells, directionality) as well as the spatial correlation properties of the speckle. Introduction of the non-stationary estimates of E[R] and CR into the scalar equation of a single-point speckle filter improves the restoration of R fluctuations in the filtered radar image [75]. As speckle and scene correlations are explicitly accounted for, this method can be considered as being complete.

The whole process acts as an adaptive focusing process in addition to the filtering process, since emphasis is put on the restoration and the enhancement of the small-spatial-scale fluctuations of the radar reflectivity. Full profit is taken of the useful resolution offered by the compound sensor and processing systems. Thus, speckle related features and thin scene elements (scene short-scale texture and thin structures) are automatically differentiated. The latter are restored with enhanced spatial accuracy, according to the local statistical context.

There is no geometrisation of the scene structural elements, as when using templates-based detectors used in structure retaining speckle filters ([76] [63] [64] [10]). The neighbourhood on which the NMNV local statistics are estimated is delimited by the estimates of ρR(Δz) for all possible Δz orientations. The only limitation to the potential extension of the domain D is the spatial extent of the available ACF of the speckle.

5. MAP speckle filters for series of detected SAR images

Since the launch of ERS-1 satellite in 1991 temporal series of calibrated SAR images have been made available. This has stimulated the development of multichannel ("vector") filters especially dedicated to filter the speckle in a series of images, taking into account the correlation of the SAR signal between image acquisitions, thus opening the way to furher developments of change detection techniques specific to SAR series of images.

Interest in multichannel adaptive speckle filtering arises from their ability to combine both multi-image diversity through the exploitation of correlation between speckle and scene between images, and spatial diversity through the locally adaptive averaging of pixel values in the spatial domain. The objectives of this combination are the simultaneously achievement of both, a better restoration of existing texture, and a stronger speckle smoothing in textureless areas.

Although numerous work with regard to multi-channel speckle filtering in SAR images had already being done in a remote past, the introduction of A Priori knowledge or A Priori guess which implies the use of Bayesian methods in the processing of multi-SAR’s images, multi-date SAR images or SAR and optical images began only in 1992. Bayesian MAP vector speckle filters developed for multi-channel SAR images incorporate statistical descriptions of the scene and of the speckle in multi-channel SAR images. These models are able to account for the scene and system effects, which result in a certain amount of correlation between the different channels.

To account for the effects due to the spatial correlation of both the speckle and the scene in SAR images, estimators originating from the local ACF’s (cf. Section 4) are incorporated to these filters to refine the evaluation of the non-stationary first order local statistics on an ergodic wide-sense neighbourhood of the pixel under processing. The goal is to improve the restoration of scene textural and structural properties and the preservation of the useful spatial resolution in the filtered SAR image.

5.1. Multichannel vector MAP speckle filtering

In the case of multi-channel detected SAR images, let define the vector quantities of interest: I is the speckled intensity vector available in the actual SAR data; R is the radar reflectivity vector, which is the quantity to restore. The MAP filtering method bases on Bayes’s theorem [47] [48] in its matricial form:

P(R/I)= P(I/R) . P(R) / P(I)
(69)

where P(I/R) is the joint conditional pdf of the speckle.

For each individual detected SAR image i, P(Ii/Ri) is known to be well approximated by a Gamma distribution for multilook intensity images (cf. § 2.1.4; Equation (14)).

P(R) is the joint pdf of the radar reflectivity, introduced as statistical A Priori information in the restoration process. To describe the first order statistical properties of natural scenes, it has already been shown that, in most situations, for each individual image i, a Gamma distribution (cf. § 3.3.4) is an appropriate choice for P(Ri). Nevertheless, when in doubt about the reliability of such a choice, P(Ri) can be estimated directly from the data in the corresponding radar image i.

For multi-channel detected SAR images, MAP filtering is a vector filtering method. For every channel i, the posterior probability is maximum if the following condition is verified:

[Ln(P(I/R))]/Ri+[Ln(P(R))]/Ri  = 0      for      Ri= R^i MAP
(70)

The MAP speckle filtering process acts as a data fusion process, since the information content of the whole image series is exploited to restore the radar reflectivity in each individual image. Among other advantages, this allows a better detection and restoration of thin scene details.

5.2. Multichannel scene models

To describe the first order statistical properties of a natural scene viewed at the spatial resolution and at the wavelength frequency of a radar remote sensing sensor, it has been shown that a Gamma pdf would be a suitable representation (cf. § 3.3.4).

However, to describe the first order statistical properties of a natural scene as viewed by diverse radar sensors (different physical scene features), or at different dates (scene evolution over time), there is no analytic multivariate Gamma pdf available under closed-form. Therefore, for the sake of mathematical tractability, a multivariate Gaussian pdf is used as analytic "ersatz" of a multi-channel scene statistical model:

P(R) =[(2π)N|CovR|]1/2. exp[(R<R>)t.CovR1.(R<R>)]  
(71)

CovR is the local covariance matrix of the imaged scene in all image channels. For each pixel location, CovR is estimated from the local covariance matrix of the intensities observed in all image channels CovI, using the multiplicative speckle model (cf. § 2.2.2; [27] [28]).

This statistical scene model is convenient to preserve the mathematical tractability of the problem. In addition, the Gaussian model, commonly used to describe the statistical properties of natural scenes in processing of optical imagery, is appropriate as a first-order statistical description of moderate, not too heterogeneous, textures.

However, when scene texture becomes very heterogeneous, this multivariate Gaussian scene model is not any longer reliable, and in the absence of specific prior knowledge about the distribution of scene texture, its model has to be directly estimated from the SAR image.

5.3. Speckle uncorrelated between SAR images: The Gamma-Gaussian MAP filter

Let consider the case of a series of N images originating from different SAR systems (different frequencies, incidence angles, or spatial resolution for example, but same spatial sampling). In such a case, it is justified to consider that speckle is independent between the N image channels.

Under this assumption, since each one of the SARs senses slightly different physical properties of the same scene, the joint conditional pdf of the fully developed speckle, P(I/R), can reasonably be modelled as a set of N independent Gamma distributions P(Ii /Ri):

P(Ii/Ri) =(Li /Ri)Li /Γ(Li). exp(LiIi /Ri) .Ii(Li1)
(72)

where the Li parameters are the ENLs of the individual SAR images.

The first MAP filtering algorithm for multi-channel detected SAR images results in a set of N scalar equations, with N independent speckle models and a coupled scene model (Equation (71)). The coupled scene model is justified by the fact that the physical properties of the scene do contribute, but not in an identical way, to the formation of the images provided by the different sensors. The set of equations describing the Gamma-Gaussian MAP filter for multi-channel detected SAR images is as follows ([73] [85]):

Li(Ii/Ri21/Ri) (1i)t.CovR1.(R<R>)(R<R>)t.CovR1.(1i)  1/2 Tr[CovR1.CovR/Ri]+(R<R>)t.CovR1.CovR/Ri.CovR1.(R<R>) = 0
(73)

where (1i) is a vector where all components, but the ith, are equal to zero. The set of Equations (73) is solved numerically.

The introduction of coupling between the scene statistical representations leads to a data fusion process taking profit of the correlation between texture as it is observed in all the images in the series. Indeed, replacing speckle noise model by optical noise (or film grain noise, or any other appropriate noise model) in some of the N image channels, this filter adapts also easily to the case of multi-channel optical and SAR images.

5.4. Speckle correlated between SAR images: The Gaussian-Gaussian MAP filter

To filter series of images originating, either from the same SAR system operating in repeat-pass mode, or from different SAR systems with relatively close properties, a second speckle filtering algorithm has been developed. The different SARs may be close in terms of frequency, angle of incidence, spatial resolution and sampling, geometry of images, with difference in polarisation configuration only, or small differences in incidence angle, for example. In such cases, speckle correlation between individual SAR images must be taken into account to deal optimally with system effects on the SAR images series.

Taking into consideration speckle correlation between image channels, the joint conditional pdf of the fully developed speckle P(I/R) should theoretically be a multivariate Gamma pdf. Nevertheless, since there is no analytic multivariate Gamma pdf available in closed-form, another reasonable choice for P(I/R) should be done for the sake of mathematical tractability.

To solve this problem, Lee’s assumption [33] is adopted: for multilook SAR images (more than 3 looks), the joint conditional pdf of the speckle can be reasonably approximated by a Gaussian distribution. Therefore, for convenience, P(I/R) is approximated by a Gaussian multivariate distribution in the case of multilook SAR images:

P(I/R) = [(2π)N|CovS|]1/2. exp[(IR)t.CovS1.(IR)]
(74)

where CovS is the covariance matrix of the speckle between images of the series. In practice, CovS is estimated within the most homogeneous area one can identify in the series of detected SAR images.

It is noteworthy that CovS takes into account possible different Li values of the respective ENLs of the SAR images in the series as well as energy unbalance between images. This enables to compose the series combining images from different SARs operating at the same wavelength, or acquired by the same SAR at different times or in different configurations of polarization.

At this point, the second MAP filtering algorithm for multi-channel detected multilook SAR images results in a set of N scalar equations. The set of equations describing the Gaussian-Gaussian MAP filter for multi-channel detected multilook SAR images is ([73] [85]):

t(1i).CovS1.(IR) +(IR)t.CovS1.(1i)  1/2 Tr[CovR1.CovR/Ri]+(R<R>)t.CovR1.CovR/Ri.CovR1.(R<R>)(1i)t.CovR1.(R<R>) t(R<R>).CovR1.(1i) = 0
(75)

where (1i) is a vector where all components, but the ith, are equal to zero. Again, this set of equations is solved numerically.

Including a coupled speckle model (Equation (71)) and a coupled scene model (Equation (74)), this speckle filter takes profit of both speckle and scene texture correlations, thus restoring the radar reflectivity through a complete data fusion process.

5.5. Prior knowledge gathered from SAR image series: DE-MAP filters

In the presence of very strong or mixed textures, eventually combined with the presence of strong topographic relief, it may be hazardous to make an A Priori assumption about the pdf of the radar reflectivity. This situation is often the case in SAR images of dense tropical forest or of urban environments located in slanted terrain with rapidly varying slopes and counter-slopes.

This is also often the case in high-and very-high spatial resolution SAR images, where strong texture is omnipresent. In addition, at such very fine scale, the textural properties of the scene vary strongly within the timeframe elapsed between successive image acquisitions by the same SAR, and exhibit strong differences when imaged by different SARs or using different SAR configurations in terms of wavelength, polarisation, or angle of incidence, for example.

5.5.1. Speckle correlated between SAR images: The Gaussian-DE MAP filter

As mentioned precedently (cf. § 5.4), for a series of N multilook (Li>3 for each image i) detected SAR images, the joint conditional pdf of the speckle, P(I/R), can be modelled as a multivariate Gaussian distribution [33] when speckle is correlated between image channels with covariance matrix CovS:

P(I/R)=[(2π)N.|CovS|]1/2. exp[t(IR).CovS1.(IR)]
(76)

The entropy of the scene texture [68] [69] [70] is introduced (on a sample of Npix pixels of an ergodic neighbourhood of the pixel under processing) as follows:

S(Ri) = kNpix[Rik. log(Rik)]     for the ithchannel
(77)

Because the radar reflectivities Rk are non-negative and exp(S(Ri))/Z is normalized, and since S(Ri) is a measure of the spread/dispersion of radar reflectivities in image channel i, Equation (77) can be treated as a pdf whose entropy is S(Ri) [70]:

P(Ri) = (1/Z. exp(S(Ri)) = (1/Z) . exp(kNpix[Rik.log(Rik)] )
(78)

Under this assumption, the Gaussian/Distribution-Entropy MAP (Gaussian-DE MAP) filter for multi-channel multilook detected SAR images (N image channels) results in a set of N coupled scalar equations of the form [71]:

(1i)t.CovS1.(IR)+(IR)t.CovS1.(1i) Ri2.kNpix[log(Rik)1/Ln(10)]= 0
(79)

where (li) is a vector where all components, but the ith, are equal to zero.

To estimate P(Ri), the radar reflectivities Rik in a neighbourhood of the pixel under processing in the ith SAR image are pre-estimated in a first speckle filtering pass; indeed, the iterative estimation process does already converge after a first application of Equation (79).

Note that this filter does not need to introduce implicitly scene texture correlation. Indeed, the restoration of the radar reflectivity in each image channel bases only on the A Priori knowledge retrieved from the speckled image channel. Nevertheless, the filter takes profit of speckle correlation between image channels. An application of this filter is illustrated in Figure 1.

media/image185_w.jpg

(a) color composition of the unfiltered HH (red) and VV (green and blue) amplitudes. (b) color composition of the speckle filtered HH and VV amplitudes obtained using the Gaussian-DE MAP filter for series of multilook detected SAR images (cf. § 5.5.1).

Figure 1.

TerraSAR-X 4-looks SAR images (HH and VV polarisations) acquired on August 10, 2007, over Oberpfaffenhofen, Germany (Credits: Astrium and InfoTerra; filtered images: Privateers NV, using ®SARScape).

5.5.2. Speckle correlated between SAR images: The Gamma-DE MAP filter

Note that if speckle is not correlated between the N image channels, the Gamma/Distribution-Entropy MAP (Gamma-DE MAP) filter for multi-channel detected SAR images (N channels) results in the resolution of a set of N independent (uncoupled) scalar equations similar to Equation (62) (cf. § 3.3.5).

5.6. Vector MAP filtering and control systems

A number of advantages and improvements these filters offer is listed in this section. Most of these advantages arise from the use of the covariance matrices of the speckle and of the imaged scene (i.e. from the use of coupled models) in a Bayesian reconstruction of the radar reflectivity.

5.6.1. Imaging system’s effects and preservation of high-resolution

In a series of SAR (but not only...) images, the resolution cells never overlap perfectly between the different individual images. The covariance matrix of the speckle CovS contains the information about these overlaps, and the filtering process accounts for these overlaps in order to preserve, and even improve the spatial resolution in the filtered SAR images, thus allowing us to restore the thinnest scene structures.

In addition, the simultaneous attempt to detect images structures and targets in all radar image channels improves the probability to detect such structures. As shown in Figure 1, such an improvement is already very substantial using only two SAR images.

5.6.2. Potential for change detection

In series of radar images acquired over time, the covariance matrix of the scene CovR contains the information about the temporal evolution of the imaged area and enables the detection of changes in scene physical properties (scene temporal evolution).

In a series of SAR images acquired by (not too...) different SAR’s, the covariance matrix of the scene enables the detection of more aspects of the scene, as it is viewed by more SAR sensors.

5.6.3. MAP speckle filters and control systems

Let consider the mathematical form of the set of Equations (73), (75) and (79). Indeed, their formulation is that of a control system, since both equations can be rewritten under the form of Riccati’s [86] continuous time algebraic matricial equations:

 A . X  X .tA  Q + X .tC . P1.C . X = 0
(80)

Equation (10) represents the optimal state controlled reconstruction at constant gain of linear invariant processes (R and the textures of the channels) perturbed by white noises (speckle, pixel spatial mismatch between channels), from observed evolving state variables (I and CovI) [87]. The scene A Priori model acts as a command, and the covariance matrices act as controls. It is also noticeable that the MAP elaboration of information generates automatically the feedback process, which enables to control the filtering process.

In addition, it is highly remarkable that Riccati’s theorem [86] stipulates the existence of a unique positive definite solution for Equation (10). Therefore, this property holds also for the MAP Equations (73), (75) and (79).

The extension of this noise filtering technique to the case of SAR and optical images set is straightforward, as mentioned in § 5.3. It is also noticeable that, extending this method to the case of complex SAR images, this technique presents a very interesting potential for super-resolution. Such methodologies can be found in the literature (cf. [46], for example).

A major interest of control systems is that they offer wide possibilities for the choice and design of additional commands (statistical and physical models) for further data exploitation. In this view, speckle filtering should be regarded as a first step of integrated application oriented control systems.

6. MAP speckle filters for complex SAR data

SLC images are complex valued images. Hence, every image pixel is assigned a complex value, its complex amplitude Ac=i+j.q=I.exp(j.ϕ), which corresponds to an intensity I=i2+q2 and a phase ϕ=Arctg(q/i).

Whereas the intensity I carries information related to the radar reflectivity R and the backscattering coefficient σ° of the imaged scene, the phase of a complex pixel of a single SLC image is totally random and does not carry any information about the scene when speckle is fully developed. (It is noteworthy that only, phase differences between SLCs acquired in interferometric conditions, or phase differences between configurations of polarisation in SAR polarimetric data, carry information about the imaged scene).

Therefore, the quantity of interest to restore in SLC images through speckle filtering is the radar reflectivity R. Thus, whereas the input of the speckle filters consists in complex radar data, their output consists in radar reflectivity images, i.e. speckle filtered intensity images.

6.1. Single-look complex (SLC) SAR image: The CGs-DE MAP filter

In SLC radar data the spatial correlation of the speckle can be dealt with by taking into account the spatial covariance matrix of the speckle in the filtering process. Indeed, the joint pdf of a local sample X (Npix-dimention vector) of 1-look spatially correlated speckle is given by Equation (7) (cf. § 2.1.2), can be written as [88] [89] [90] [91]:

Pu(X) = P(X/R) = [1/(πNpix.|CS|)] .exp [Xt*.CS1.X]
(81)

where CS is the spatial covariance matrix of the complex speckle, estimated within the most homogeneous part of the SLC radar image.

In high and very-high spatial resolution SLC SAR images, which often exhibit strong textural properties, it becomes difficult to consider any theoretical statistical model as a reliable A Priori knowledge about scene texture. In such a situation, it seems reasonable to retrieve statistical scene knowledge directly from the SAR image data [68] [69]Introducing the local entropy of the radar reflectivity as A Priori scene knowledge [70] (cf. § 3.3.5), the Complex-Gaussian/Distribution-Entropy MAP (CGs-DE MAP) filter for an SLC radar image is the solution of the following equation [71]:

(1/Npix) . tX*.CS1.X R  R2.k=1Npix[log(Rk)1/Ln(10)] = 0 
(82)

where Npix is the effective number of neighbour pixels taken into account in the computation of local statistics. The radar reflectivites Rk are pre-estimated in this neighbourhood of the pixel under processing by a first pass of speckle filtering in intensity.

An example of application of this filter is shown in Figure 2.

media/image190_w.jpg

(a) unfiltered 1-look SAR image amplitude. (b) filtered SAR image amplitude obtained using the Complex-Gaussian-DE MAP filter for SLC SAR image (cf. § 6.1).

Figure 2.

RADARSAT-2 1-look SAR image (HH polarisation) acquired on April 30, 2009, over Capetown, South Africa (Credits: Canadian Space Agency and MDA; filtered images: Privateers NV, using ®SARScape).

6.2. Separate complex looks, or series of SLC SAR images

Interest in adaptive speckle filtering of SLC image series arises from their ability to combine both temporal or spectral diversity, and spatial diversity. The objective is to obtain a multilook radar reflectivity image, where additional speckle reduction is obtained in the spatial domain by averaging pixel values through locally adaptive filtering.

The final objective is to obtain very strong speckle smoothing in textureless areas and to reach an ENL high enough (i.e. ENL>>30, and preferably ENL>150) to enable information extraction using classic techniques, or combined use with optical imagery, while better restoring existing or characteristic scene texture.

6.2.1. Compount statistical speckle model for series of SLC SAR images

A SAR sensor can produce series of SLC SAR images in three different situations:

  1. L separate complex looks are usually extracted from the useful Doppler band of the radar signal, usually to produce a L-looks multilook image, but the individual looks may be conserved at this stage of SAR image production. Let remind that this operation is done at the cost of a loss in azimuth spatial resolution (cf. § 2.1.3). These looks generally overlap, in the sense that they have in common a part of the Doppler spectrum of the original full-resolution 1-look signal, in order to preserve as completely as possible the total amount of backscattered radar signal. This overlap introduces a certain amount of speckle correlation between looks extracted from adjacent portions of the Doppler spectrum.

  2. In the case of a non-interferometric 1-look complex SAR data set, the individual SAR SLC images acquired independently of each other show low or inexistent speckle correlation between images.

  3. However, in the case of interferometric 1-look complex SAR data, speckle correlation is close to 1 (corresponding to perfect correlation) or very high, since speckle correlation between SAR data acquisitions is exactly the fundamental requirement for interferometry.

As a consequence of these three possible configurations of a series of SLC images, a speckle filter dealing with the whole series must be able to take into account the complex correlation of speckle between the different images. The local complex covariance matrix CovS describes completely both, the radar intensities in each SLC image, and the correlation between images.

Considering the whole set of SLC images, the measurement vector for each pixel is X={yn}, where yn=in+j.qn. When speckle is fully developed, the (in, qn) are statistically independent random processes. However, the yn are correlated complex Gaussian random processes with a joint conditional pdf given by [88]:

P(X/CovS) = exp(tX*.CovS1.X) / (πL.|CovS|)
(83)

where CovS=E[tX. X*] is the complex covariance matrix of the speckle between SLC images, |CovS| is the determinant of CovS, and (CovS-1/|CovS|) is the inverse complex correlation matrix of the speckle. Theoretically, it is possible to compute the correlation matrix of the speckle from the sensor’s and SAR processor’s parameters. But, since these parameters are in general not easily available to SAR images users, CovS is estimated in practice over the most homogeneous area identified in the images.

6.2.2. Gamma distributed scene model: The Complex-Gaussian–Gamma MAP filter

In the case of separate complex looks corresponding to the same SAR data acquisition, or of series of SLC images acquired over time by the same SAR over a time-invariant scene, it is quite reasonable to assume that the textural properties of the scene are similar in all the SLC’s.

Therefore, scene texture may be locally modelled for all SLCs in the series, as well as in the detected multilook image formed by incoherent pixel-per-pixel averaging, by the same Gamma pdf with parameters E[R] and α (cf. § 3.3.4, Equation (52)). Locally, E[R] is estimated, either by <I> or using Equation (67), and α is estimated by 1/CR2 (cf. § 2.2.2, Equations (26) and (27)), in an ergodic wide-sense stationary neighbourhood of the pixel under processing in the (less noisy) multilook radar image.

In this situation, the Complex-Gaussian/Gamma MAP (CGs-Gamma MAP) filter for separate complex looks (L 1-look SAR images) has the following expression [11]

R^=<I>.(αL1) + <I>2.(αL1)2+ 4α.<I>.(tX*.(CovS1/|CovS|).X) 2α
(84)

Note that when the SLC images or the separate looks are independent, i.e. when the non-diagonal terms of CovS are equal to zero, then tX*.(CovS-1/|CovS|).X=L.<I>, and the filter is equivalent to the Gamma-Gamma MAP filter for single detected SAR image (cf. § 3.3.4). Therefore, independent SLC images can be incoherently summed to produce an L-looks intensity image that is filtered using the Gamma-Gamma MAP filter (Equation (54)) with the same α and R local parameters, for the same final result.

Finally, it is noteworthy that the ENL L’ of the multilook radar image can be computed from the complex correlation coefficients ρmn between images m and n, computed from the covariance matrix of the speckle, i.e. CovS estimated over the most homogeneous area identified in the images [11]:

L = L . (1 +(2/L) .n=1L1m=n+1L|ρmn|2)1
(85)

This relationship may prove useful if SLC images are incoherently summed to produce a L’-looks intensity image.

6.2.3. Prior knowledge gathered from the data: The Complex-Gaussian-DE MAP filter

High and very-high spatial resolution SLC SAR images often exhibit strong textural properties. As exposed above, in such a situation, statistical scene knowledge is estimated from the radar data [68] [69]Introducing the local entropy of the radar reflectivity as A Priori scene knowledge [70] (cf. § 3.3.5), the Complex-Gaussian/Distribution-Entropy MAP (CGs-DE MAP) filter for an SLC image series (L images) is the solution of the following equation [71]:

tX*.CovS1.X L.Ri Ri2.k=1Npix[log(Rik)1/Ln(10)] = 0
(86)

where CovS is the covariance matrix of the speckle between the different SLC images [11] and Npix is the effective number of neighbour pixels taken into account in the computation of local statistics. The radar reflectivites Rk are pre-estimated in this neighbourhood of the pixel under processing by a first pass of multichannel speckle filtering of the L images in intensity (using Equation (84) for example).

An application of this filter to a series of 18 (3x6) SLC SAR images is shown in Figure 3.

media/image195_w.jpg

(a) color composition of three unfiltered 6-looks SAR images, each one obtained by pixel-per-pixel averaging of 6 individual 1-look images. (b) color composition of three speckle filtered images, each one obtained using the Complex-Gaussian-DE MAP filter for series (6 images) of SLC SAR images (cf. § 6.2.3).

Figure 3.

Series of 18 (3x6) COSMO-SkyMed 1-look spotlight SAR images (HH and VV polarisations) acquired between February 2 and September 30, 2009, over Perito Moreno in Patagonia, Argentina (Credits: original images: Italian Space Agency; filtered images: Privateers NV, using ®SARScape).

7. MAP speckle filters for polarimetric radar data

7.1. Polarimetric radar data representation, and polarimetric speckle model

A polarimetric radar system produces, for each pixel location, a scattering measurement matrix, which is the scattering matrix S of the corresponding surface (cf. Section 1, Equation (1)) corrupted by speckle noise.

The first polarimetric speckle filter ever developed [92] resulted in an optimal summation of the intensities IHH, IHV, IVV of the polarisation channels. The resulting image is called the "improved span" image:

ISpan= IHH+(1+|ρHHVV|2) .IHV. E(IHH)/E(IHV)+ IVV.E(IHH)/E(IVV)
(87)

where ρHHVV is the complex correlation coefficient between SHH and SVV.

However, Equation (87) clearly shows that the ENL achieved in the span image will barely reach 3, which remains largely insufficient to meet the noise reduction requirements of remote sensing applications. Besides, the physical interpretation of ISpan, which is a mixture of HH, VV, HV radar reflectivities is questionable. In addition, scene information contained in the polarimetric diversity of the radar measurement as well as in the complex correlations between polarimetric channels is lost in the process.

Lee et al. [93] developed a polarimetric optimally combining the polarisation channels to restore the radiometric information in the HH, VV, and HV channels. This polarimetric speckle filter has been shown to preserve polarimetric signatures [94].

Nevertheless, the ultimate aim of polarimetric speckle filtering is to restore the full polarimetric scattering matrix S, or at least the radar reflectivities in the HH, HV, VV configurations of polarisation, while achieving a very high ENL value (>100) in the filtered data. To achieve these objectives, fully polarimetric adaptive speckle filters take profit of both the entire polarimetric diversity (i.e. the full scattering matrix) and the spatial diversity (i.e. local statistics around the location under processing) [95].

7.1.1. Polarimetric covariance matrix

Under the assumption of reciprocity, SHV=SVH, the speckle-free (denoised) scattering matrix S can be transformed into its covariance matrix CS, as follows:

CS=(|SHH|2SHH.SHV*SHH.SVV*SHH*.SHV|SHV|2 SHV.SVV* SHH*.SVVSHV*.SVV|SVV|2)
(88)

This representation puts emphasis on the radar reflectivities RHH=|SHH|2, RHV=|SHV|2 and RVV=|SVV|2, as well as the covariances between HH, HV and VV configurations of polarisation.

If the polarimetric measurement is made by a monostatic radar system, the covariance matrix ΣS of the actual polarimetric radar measurement is the speckle-corrupted version of CS.

It is noteworthy that it is possible to produce multilook polarimetric data (L-looks) from the original 1-look data. For L-looks polarimetric data obtained by either spatial averaging (generally in azimuth) or spectral multilooking (cf. § 1.1.4), ΣS multilook is defined by:

ΣS multilook = [ (1/L).1L(Spq.Spqt*) ]      for all combinations p,q of H,V
(89)

7.1.2. Polarimetric vector, degrees of coherence, phase differences

For convenience, the polarimetric covariance matrix is often expressed, without changing its overall information content, under the form of a real valued vector, called the "polarimetric feature vector", or "polarimetric vector" X:

tX = ( |SHH|2, |SHV|2, |SVV|2, Re(SHH.SVV*), Im(SHH.SVV*),Re(SHH.SHV*), Im(SHH.SHV*), Re(SVV.SVV*), Im(SVV.SVV*) )   
(90)

Two quantities of great interest in polarimetric SAR applications can be obtained from the polarimetric vector:

1) The Phase Differences ΔφKL-MN between radar returns in configurations of polarisation KL and MN (where K,L,M,N ∈ {H,V} are preferably considered in most applications), which are computed as follows:

ΔφKL-MN= Arg(SKL.SMN*)= Arctan[Im(SKL)/Re(SKL)]Arctan[Im(SMN)/Re(SMN)] 
(91)

2) The (complex) Degrees of Coherences γKL-MN, which measure wave coherence between radar returns in two configurations of polarisation KL and MN, and express textural properties of an extended wide-sense stationary scene, are defined as follows:

γKL-MN= E[Re(SKL.SMN*)]+j. E[Im(SKL.SMN*)](E[|SKL|2].E[|SMN|2] 
(92)

However, in the actual polarimetric radar measurement, the observations ΔφKL-MN of the phases differences, and the observations ρKL-MN of the complex degrees of coherence are the speckle-corrupted versions of Arg(SKL.SMN*) and γKL-MN, respectivly. Therefore, these quantities, which represent important contributions of radar polarimetry must also be speckle filtered, justifying the development of fully-polarimetric speckle filters that do not limit themselves to the restoration of radar reflectivities.

7.1.3. Polarimetric speckle model

In the case of polarimetric radar data, ΣS is the actually observed polarimetric covariance matrix, and CS is the unspeckled polarimetric covariance matrix, i.e. the quantity to restore through speckle filtering. In the reciprocal case, and for low look correlation, the conditional pdf of ΣS is a complex Wishart distribution of the form [72]:

P(ΣS/CS)=(det ΣS)L3.LL3π3.Γ(L).Γ(L1).Γ(L2).(det CS)L . exp[Tr(L.CS1.ΣS)]  
(93)

where L is the ENL of the "improved span" image defined by Equation (87) [92].

7.2. Restoration of the whole polarimetric vector / covariance matrix

7.2.1. Gamma distributed scene model: The Wishart-Gamma MAP filter

Using physical backscattering models, assuming (as a rough approximation) that texture is identical in all polarizations, we get the following approximation [72]:

CS=μ . E[CS]
(94)

where μ is the scalar texture parameter equal to the normalized number of scatterers within the resolution cell, and E(CS) is the mean covariance matrix [92].

P(μ) = αα/Γ(α) . exp(α.μ) . μα1     and     E[μ]= 1
(95)

The characterisation of scene heterogeneity by only one textural random variable μ bases on the assumption that textural properties are identical in all configurations of polarisation. However, this has been shown to be inexact, both experimentally [96] and theoretically [97]. Therefore, μ should be regarded as an average textural characterisation of the scene observed in all possible configurations of polarisation. Besides the advantage of estimating α=1/ CR2 (using Equation (27)) in a less noisy multilook image is an additional justification to do it in the "improved span" image defined by Equation (87) [92].

Introducing the first-order statistical models for fully-polarimetric speckle (Equation (93)) and the polarimetric texture parameter μ in the MAP equation (Equation (50), cf. § 3.3.2), the fully-polarimetric Wishart-Gamma MAP filter restores the value of μ as:

μ^=(αL1) + (αL1)2+ 4.α.L.Tr(E[CS]1.ΣS) 2α
(96)

where Tr(.) denotes the trace of a matrix.

Finally, the restored (speckle filtered) version of the covariance matrix CS is obtained using the maximum likelihood estimator described in [72]:

C^S= μ^ .E[CS]
(97)

7.2.2. Prior knowledge gathered from polarimetric SAR data: The Wishart-DE MAP filter

In high and very-high spatial resolution polarimetric SAR data, strong and/or mixed textural properties justify to estimate statistical scene knowledge from the data themselves, rather than assuming a A Priori theoretical model. Assuming that textural properties are identical in all configurations of polarisation, the entropy constraint on scene texture ([68] [69] [70]) becomes:

P(μ) = (1/μ). exp(k[μk.log(μk)]           and     E[μ] = 1
(98)

In this case, the complex Wishart/Distribution-Entropy MAP (CW-DE MAP) filter for polarimetric multilook SAR data is expressed as [71]:

L.Tr(E[CS]1.ΣS)  L.μμ2.k[log(μk)1/Ln(10)] = 0
(99)

E[CS] is obtained using the maximum likelihood estimator (Equation (97)) described in Lopès et al. [72].

media/image211_w.jpg

(a) color composition of the unfiltered HH (red), HV (green), and VV (blue) amplitudes. (b) color composition of the speckle filtered HH, HV, and VV amplitudes obatained using the fully polarimetric Wishart-DE MAP filter (cf. § 7.2.2). (c) color composition of the unfiltered HH-VV (red), HH-HV (green) and HV-VV (blue) degrees of coherence. (d) color composition of the HH-VV, HH-HV, and HV-VV degrees of coherence obtained from the Wishart-DE MAP speckle filtered polarimetric vector (cf.cf. § 7.2.2). (e) Unfiltered HH-VV phase difference. (f) HH-VV phase difference obtained from the Wishart-DE MAP speckle filtered polarimetric vector (cf.cf. § 7.2.2).

Figure 4.

ALOS-PALSAR 6-looks polarimetric SAR imagery acquired on June 30, 2006 over Bavaria, Germany (Credits: JAXA and MITI; filtered images: Privateers NV, using ®SARScape).

The texture parameters μ k in the neighbourhood of the pixel under processing are pre-estimated over an ergodic neighbourhood of the pixel under processing by a first speckle filtering pass using the Wishart-Gamma MAP filter (cf. § 7.2.1).

Although these fully polarimetric MAP filters assume identical texture properties in the HH, HV, and VV channels, which has been shown both experimentally [96] and theoretically [97] inexact, they have nevertheless been shown to preserve polarimetric signatures [98]. Fully polarimetric speckle filtering is illustrated in Figure 4.

Acronyms

ACF: spatial AutoCorrelation Function (of a signal)

DE: Distribution-Entropy (statistical distribution)

ENL: Equivalent Number of Looks (of a radar image)

LMMSE: Linear Minimum Mean Square Error (estimation)

MAP: Maximum A Posteriori (statistical estimation)

MMSE: Minimum Mean Square Error (estimation)

NMNV: Non-stationary Mean Non-stationary Variance (model)

pdf: Probability Density Function (random variable distribution)

Pfa: Probability of False Alarm (detection)

PSF: Point Spread Function (of an imaging system)

RoA: Ratio Of Amplitudes (structure detector)

SAR: Synthetic Aperture Radar (imaging sensor)

SLC: Single-Look-Complex (radar image)

References

1 - F.T. Ulaby, A.K. Fung and R.K. Moore, 1982: "Microwave Remote Sensing", Vol.1 & 2, Artech House Publishers, 1982.
2 - J.C. Dainty, 1975: "Statistical properties of laser speckle patterns", in Laser Speckle and Related Phenomena, Springer Verlag, Berlin (Germany), 1975.
3 - M. Tur, K.C. Chin and J.W. Goodman, 1982: "When is speckle noise multiplicative?", Applied Optics, Vol.21, n°7, pp.1157-1159, April 1982.
4 - R.K. Raney, 1983: "Transfert function for partially coherent SAR system", IEEE Trans. on AES, Vol.AES-19, n°5, pp.740-750, September 1983.
5 - J.W. Goodman, 1976: "Some fundamental properties of speckle", J. Opt. Soc. Am., Vol.66, n°11, pp.1145-1150, November 1976.
6 - A. Lopès and R. Touzi, 1986: "Extraction of the backscatter coefficient of agricultural fields from an airborne SAR image", Proc. of IGARSS'86, Zürich (Switzerland), ESA SP-254, Vol.3, pp.621-627, 8-11 September 1986.
7 - R. Touzi, 1988: "Analyse d'images radar en télédétection: améliorations radiométriques, filtrage du speckle et détection des contours", Ph.D dissertation, Université Paul Sabatier (Toulouse III, France), n°258, 4 March 1988.
8 - A.L. Gray, P.W. Vachon, C.E. Livingstone and T.I. Lukowski, 1990: "Synthetic aperture radar calibration using reference reflectors", IEEE Trans. on GRS, Vol.GE-28, n°3, pp.374-383, May 1990.
9 - V.S. Frost and K.S. Shanmugan, 1983: "The information content of synthetic aperture radar images of terrain", IEEE Trans. on AES, Vol.AES-19, n°5, pp.768-774, September 1983.
10 - A. Lopès, E. Nezry, R. Touzi and H. Laur, 1993: "Structure detection and statistical adaptive speckle filtering in SAR images", Int. J. Rem. Sens., Vol.14, n°9, pp.1735-1758, June 1993.
11 - A. Lopès, E. Nezry, S. Goze, R. Touzi and G. Aarbakke-Solaas, 1992: "Adaptive prossessing of multilook complex SAR images", Proc. of IGARSS'92, Clear Lake (USA), Vol.2, pp.890-892, 26-29 May 1992.
12 - H.H. Arsenault and M. Levesque, 1984: "Combined homomorphic and local-statistics processing for restoration of images degraded by signal dependant noise", Applied Optics, Vol.23, n°6, pp.845-850, March 1984.
13 - H.H. Arsenault, 1987: "Information extraction from images degraded by speckle", Proc. of IGARSS'87, Ann Arbor (USA), Vol.2, pp.1317-1320, 18-21 May 1987.
14 - Ping Fan Yan and C.H. Chen, 1986: "An algorithm for multiplicative noise in wide range", Traitement du signal, Vol.3, n°2, pp.91-96, 1986.
15 - H.H. Arsenault and G. April, 1976: "Properties of speckle integrated with a finite aperture and logarithmically transformed", J. Opt. Soc. Am., Vol.66, n°11, pp.1160-1163, November 1976.
16 - E. Nezry, 1988: "Evaluation comparative des méthodes de filtrage du speckle sur les images radar (SAR)", Rapport de DEA, CESR - Université Paul Sabatier (Toulouse, France), June 1988.
17 - H.H. Arsenault and G. April, 1985: "Information content of images degraded by speckle noise", Proc. of SPIE, Vol.556, International conference on speckle, pp.190-195, 1985.
18 - D.H. Hoekman, 1991: "Speckle ensemble statistics of logarithmically scaled data", IEEE Trans. on GRS, Vol.GE-29, n°1, pp.180-182, January 1991.
19 - H.H. Arsenault and G. April, 1986: "Information content of images degraded by speckle noise", Optical Engineering, Vol.25, n°5, pp.662-666, May 1986.
20 - A. Lopès, 1983: "Etude expérimentale et théorique de l'atténuation et de la rétrodiffusion des micro-ondes par un couvert de blé. Application à la télédétection", Ph.D. dissertation, Université Paul Sabatier (Toulouse III, France), n°852, 25 October 1983.
21 - H. Laur, 1989: "Analyse d'images radar en télédétection: discriminateurs radiométriques et texturaux", Ph.D. dissertation, Université Paul Sabatier (Toulouse III, France), n°403, 23 March 1989.
22 - G. Kattenborn, E. Nezry and G. DeGrandi, 1993: "High resolution detection and monitoring of changes using ERS-1 time series", Proc. of the 2nd ERS-1 Symposium, 11-14 Oct. 1993, Hamburg (Germany), ESA SP-361, Vol.1, pp.635-642.
23 - E. Nezry, E. Mougin, A. Lopès, J.-P. Gastellu-Etchegorry and Y. Laumonier, 1993: "Tropical vegetation mapping with combined visible and SAR spaceborne data", Int. J. Rem. Sens., Vol.14, n°11, pp.2165-2184, 20 July 1993.
24 - E. Nezry, A. Lopès, D. Ducros-Gambart, C. Nezry and J.S. Lee, 1996: "Supervised classification of K-distributed SAR images of natural targets and probability of error estimation", IEEE Trans. on GRS, Vol.34, n°5, pp.1233-1242, September 1996.
25 - D.L.B. Jupp, A.H. Strahler and C.E. Woodcock, 1988: "Autocorrelation and regularization in digital images : I. Basic theory", IEEE Trans. on GRS, Vol.GE-26, n°4, pp.463-473, July 1988.
26 - G. April and H.H. Arsenault, 1984: "Non-stationary image-plane speckle statistics", J. Opt. Soc. Am. A, Vol.1, n°7, pp.738-741, July 1984.
27 - F.T. Ulaby, A.K. Fung and R.K. Moore, 1986: "Microwave Remote Sensing", Vol.3, Artech House Publishers, 1986.
28 - F.T. Ulaby, F. Kouyate, B. Brisco and T.H. Lee Williams, 1986: "Textural information in SAR images", IEEE Trans. on GRS, Vol. GE-24, n°2, pp-235-245, March 1986.
29 - H. Laur, T. Le Toan and A. Lopès, 1987: "Textural segmentation of SAR images using first order statistical parameters", Proc. of IGARSS'87, Ann Arbor (USA), Vol.2, pp.1463-1468, 18-21 May 1987.
30 - V.S. Frost, J.A. Stiles, K.S. Shanmugan and J.C. Holtzman, 1982: "A model for radar images and its application to adaptive digital filtering of multiplicative noise", IEEE Trans. on PAMI, Vol.PAMI-4, n°2, March 1982.
31 - R.K. Raney and G.J. Wessels, 1988: "Spatial considerations in SAR speckle simulation", IEEE Trans. on GRS, Vol.GE-26, n°5, pp.666-672.
32 - D.T. Kuan, A.A. Sawchuk and T.C. Strand, 1987: "Adaptive restoration of images with speckle", IEEE Trans. on ASSP, Vol.ASSP-35, n°3, pp.373-383, March 1987.
33 - J.S. Lee, 1980: "Digital image enhancement and noise filtering by use of local statistics", IEEE Trans. on PAMI, Vol.PAMI-2, pp.165-168, March 1980.
34 - J.S. Lim and H. Nawab, 1981: "Techniques for speckle noise removal", Optical Engineering, Vol.20, n°3, pp.472-480, May 1981.
35 - I. Ghikhman and A. Skorokhod, 1977-1980: "Introduction à la théorie des processus aléatoires", Chapt. 5, pp.290-305, Editions MIR, Moscow (USSR), in Russian 1977, in French 1980.
36 - D.S. Zrnic, 1977: "Mean power estimation with a recursive filter", IEEE Trans. on AES, Vol.AES-13, n°3, pp.281-289, May 1977.
37 - G. Horgan, 1994: "Choosing weight functions for filtering SAR", Int. J. Rem. Sens., Vol.15, n°5, pp.1053-1164.
38 - A. Yaglom, 1955: "Introduction to the theory of random stationary functions", in Russian, Ouspekhi Matematicheskoie Naouk, Vol. 7, n°5, pp.3-168, 1955.
39 - G.C. Deane, S. Quegan, et P.N. Churchill, A. Hendry, A. Evans and J.W. Trewett, 1986: "Land-use features detection in SAR images", ESA SP-257, 1986.
40 - A. Lopès, R. Touzi and E. Nezry, 1990: "Adaptive speckle filters and scene heterogeneity", IEEE Trans. on GRS, Vol.GE-28, n°6, pp.992-1000, November 1990.
41 - D.T. Kuan, A.A. Sawchuk, T.C. Strand and P. Chavel, 1985: "Adaptive noise smoothing filter for images with signal-dependant noise", IEEE Trans. on PAMI, Vol.PAMI-7, n°2, pp.165-177, March 1985.
42 - J.S. Lee, 1981: "Speckle analysis and smoothing of synthetic aperture radar images", Computer Graphics and Image Processing, n°17, pp.24-32, 1981.
43 - J.S. Lee, 1983: "A simple speckle smoothing algorithm for synthetic aperture radar images", IEEE Trans. on SMC, Vol.SMC-13, n°1, pp.85-89, January 1983.
44 - J.S. Lee, 1986: "Speckle suppression and analysis for synthetic aperture radar images", Optical Engineering, Vol.25, n°5, pp.636-643, May 1986.
45 - J.S. Lee, 1987: "Statistical modelling and suppression of speckle in synthetic aperture radar images", Proc. of IGARSS'87, Ann Arbor (USA), Vol.2, pp.1331-1339, 18-21 May 1987.
46 - S.P. Luttrell, 1991: "The theory of bayesian super-resolution of coherent images: a review", Int. J. Rem. Sens., Vol.12, n°2, pp.303-314, 1991.
47 - T. Bayes and R. Price, 1763: "An essay towards solving a problem in the doctrine of chance. By the late Rev. Mr. Bayes, communicated by Mr. Price, in a letter to John Canton, A.M.F.R.S.", Philosophical Trans. of the Royal Society of London, Vol.53, pp.370–418.
48 - P.-S. Laplace, 1774: "Mémoire sur la probabilité des causes par les événements", Mémoires de l'Académie royale des Sciences de Paris (Savants étrangers), Vol.4, pp.621-656.
49 - J.M. Keynes, 1929: "A treatise on probability", Editions Macmillan, London (UK), 1929.
50 - R.P. Cox, 1946: "Probability, frequency and reasonable expectation", American Journal of Physics, Vol.14, n°1, pp.1-13, January-February 1946.
51 - H.L. Van Trees, 1968: "Detection, Estimation, and Modulation Theory", Editions J. Wiley & Sons, New York, 1968.
52 - B.R. Hunt, 1977: "Bayesian methods in non-linear digital image restoration", IEEE Trans. on Computers, Vol.C-26, n°3, pp.219-229, March 1977.
53 - D.T. Kuan, A.A. Sawchuk, T.C. Strand and P. Chavel, 1982: "MAP Speckle reduction filter for complex amplitude speckle images", The Proceedings of the Pattern Recognition and Image Processing Conference, pp.58-63, IEEE, 1982.
54 - S. Geman and D. Geman, 1984: "Stochastic relaxation, Gibbs distribution, and the Bayesian restoration of images", IEEE Trans. on PAMI, Vol.PAMI-6, n°6, pp.721-741, November 1984.
55 - J.R. Garside and C.J. Oliver, 1988: "A comparison of clutter texture properties in optical and SAR images", Proc. of IGARSS'88, Edinburgh (Scotland), ESA SP-284, Vol.2, pp.1249-1255, 13-16 September 1988.
56 - G. Edwards, R. Landry and K.P.B. Thompson, 1988: "Texture analysis of forest regeneration sites in high-resolution SAR imagery", Proc. of IGARSS'88, Edinburgh (Scotland), ESA SP-284, Vol.3, pp.1355-1360, 13-16 September 1988.
57 - A. Lopès, H. Laur and E. Nezry, 1990: "Statistical distribution and texture in multilook and complex SAR images", Proc. of IGARSS'90, Washington D.C. (USA), Vol.3, pp.2427-2430, 20-24 May 1990.
58 - E. Jakeman and P.N. Pusey, 1976: "A model for non-Rayleigh sea echo", IEEE Trans. on AP, Vol.AP-24, n°6, pp.806-814, November 1976.
59 - C.J. Oliver, 1984: "A model for non-Rayleigh scattering statistics", Optica Acta, Vol.31, n°6, pp.701-722, 1984.
60 - C.J. Oliver, 1988: "Representation of radar sea clutter", IEE Proceedings, Vol.135, Pt. F, n°6, pp.497-500, December 1988.
61 - E. Jakeman and R.J.A. Tough, 1987: "Generalized K distribution: a statistical model for weak scattering", J. Opt. Soc. Am. A, Vol.4, n°9, pp.1764-1772, September 1987.
62 - D. Lindley, 1991: "Making Decisions", 2nd edition, J. Wiley & Sons Inc. Ed., p.104, 1991.
63 - E. Nezry, A. Lopès and R. Touzi, 1991: "Detection of structural and textural features for SAR images filtering", Proc. of IGARSS'91, Espoo (Finland), Vol.4, pp.2169-2172, 3-6 June 1991.
64 - E. Nezry, 1992: "Restauration de la réflectivité radar pour l'utilisation conjointe des images radar et optiques en télédétection", Ph.D. dissertation, Université Paul Sabatier (Toulouse III, France), n°1236, 237 p., 3 July 1992.
65 - A. Lopès, E. Nezry, R. Touzi and H. Laur, 1990: "Maximum a posteriori speckle filtering and first order texture models in SAR images", Proc. of IGARSS'90, Washington D.C. (USA), Vol.3, pp.2409-2412, 20-24 May 1990.
66 - R. Touzi, A. Lopès and P. Bousquet, 1988: "A statistical and geometrical edge detector for SAR images", IEEE Trans. on GRS, Vol.GRS-26, n°6, pp.764-773, November 1988.
67 - E. Nezry, G. De Grandi and M. Leysen, 1995: "Speckle and scene spatial statistical estimators for SAR image filtering and texture analysis: Some application to agriculture, forestry and point targets detection", Proc. of SPIE, Vol.2584, pp.110-120, 25-29 September 1995.
68 - M. Datcu and M. Walessa, 1997: "Maximum entropy methods for despeckleling and resampling of synthetic aperture radar images of rough terrain", Proc. of SPIE, Vol.3217, pp.76-83, September 1997.
69 - M. Walessa and M. Datcu, 1997: "Bayesian approach to SAR image reconstruction", Proc. of IGARSS’97, Vol.2, pp.767-770, 03-08 August 1997.
70 - A.K. Jain, 1989: "Fundamentals of digital image processing", Chapt. 8, pp.316-317, Prentice Hall Ed., Englewood Cliffs (NJ), 1989.
71 - E. Nezry and F. Yakam-Simen, 1999: "Five new distribution-entropy MAP speckle filters for polarimetric SAR data, and for single or multi-channel detected and complex SAR Images", Proc. of SPIE, Vol.3869, pp.100-105, 20-24 September 1999.
72 - A. Lopès, S. Goze and E. Nezry, 1992: "Polarimetric speckle filters for SAR data", Proc. of IGARSS'92, Houston (TX), Vol.1, pp.80-82, 26-29 May 1992.
73 - E. Nezry, F. Zagolski, A. Lopès and F. Yakam-Simen, 1996: "Bayesian filtering of multi-channel SAR images for detection of thin structures and data fusion", Proc. of SPIE, Vol.2958, pp.130-139, September 1996.
74 - F. Yakam-Simen, E. Nezry and J. Ewing, 1998: "The legendary lost city 'Ciudad Blanca' found under tropical forest in Honduras, using ERS-2 and JERS-1 SAR imagery", Proc. of SPIE, Vol.3496, pp.21-28, September 1998.
75 - E. Nezry, H.G. Kohl and H. De Groof, 1994: "Restoration and enhancement of textural properties in SAR images using second order statistics", Proc. of SPIE, Vol.2316, pp.115-124, Rome (Italy), 26-30 September 1994.
76 - J.S. Lee, 1981: "Refined filtering of image noise using local statistics", Computer Graphics and Image Processing, n°15, pp.380-389, 1981.
77 - R.G. Caves, P.J. Harley and S. Quegan, 1992: "Edge structures in ERS-1 and airborne SAR data", Proc. of IGARSS'92, Clear Lake (TX), Vol.2, pp.1117-1119, 26-29 May 1992.
78 - R.P.H.M. Schoenmakers, 1995: "Segmentation of Remotely Sensed Imagery", Ph.D. dissertation, Katholieke Universiteit Nijmegen (NL), EUR 16087 EN, Office for Official Publications of the E.C., Luxembourg, 168 p., 13 Sept. 1995.
79 - A. Lopès and F. Séry, 1994: "The LMMSE polarimetric Wishart and the LMMSE textural vector speckle filters for SAR images". Proc. of IGARSS'94, Pasadena (USA), Vol.4, pp.2143-2145, 8-12 August 1994.
80 - J.W. Woods and J. Biemond, 1982: "Comments on 'A model for radar images and its application to adaptive digital filtering of multiplicative noise'", IEEE Trans. on PAMI, Vol.PAMI-4, n°2, pp.168-169.
81 - H.-Ch. Quelle and J.-M. Boucher, 1990: "Combined use of parametric spectrum estimation and Frost-algorithm for radar speckle filtering", Proc. of IGARSS'90, Washington D.C. (USA), Vol.1, pp.295-298, 20-24 May 1990.
82 - E. Rignot and R. Kwok, 1993: "Characterization of spatial statistics of distributed targets in SAR data", Int. J. Rem. Sens., Vol.14, n°2, pp.345-363.
83 - C.J. Oliver, 1985: "Correlated K-distributed clutter models", Optica Acta, Vol.32, n°12, pp.1515-1547.
84 - O.D. Faugeras and W.K. Pratt, 1980: "Decorrelation methods of texture feature extraction", IEEE Trans. on PAMI, Vol.PAMI-2, n°4, pp.323-332.
85 - E. Nezry, F. Yakam-Simen, F. Zagolski and I. Supit, 1997: "Control systems principles applied to speckle filtering and geophysical information extraction in multi-channel SAR images", Proc. of SPIE, Vol.3217, pp.48-57, September 1997.
86 - J. Riccati, 1724: "Animadversiones in aequationes differentiales secundi gradus", Actorum Eruditorum, quae Lipsiae publicantur, Supplementa, Vol.8, pp.66-73.
87 - P. Lancaster and L. Rodman, 1995: "Algebraic Riccati equations", Oxford University Press, 1995.
88 - J.W. Goodman, 1985: "Statistical optics", J. Wiley & Sons Inc. Ed., New York (USA), 1985.
89 - R. Fjørtoft, A. Lopès, J. Bruniquel and P. Marthon, 1998: "Optimal edge detection in SLC images with correlated speckle", Proc. of IGARSS'98, Seattle (USA), Vol., pp., 6-10 July 1998.
90 - R. Fjørtoft, J.-C. Cabada, A. Lopès and P. Marthon, 1998: "Optimal edge detection and segmentation of SLC SAR images with spatially correlated speckle", Proc.of SPIE, Vol.3497, Barcelona (Spain), 21-25 September 1998.
91 - R. Fjørtoft, A. Lopès, J. Bruniquel and P. Marthon, 1999: "Optimal edge detection and edge localization in complex SAR images with correlated speckle", IEEE Trans. on GRS, Vol.37, pp.2272-2281, September 1999.
92 - L.M. Novak and M.C. Burl, 1990: "Optimal speckle reduction in polarimetric SAR imagery", IEEE Trans. on AES, Vol.26, pp.293-305, March 1990.
93 - J.S. Lee, M.R. Grunes and S.A. Mango, 1991: "Speckle reduction in multipolarization multifrequency SAR imagery", IEEE Trans. on GRS, Vol.GE-29, n°4, pp.535-544, July 1991.
94 - J.S. Lee, G. De Grandi, M.R. Grunes and E. Nezry, 1996: "Polarimetric signature preservation in SAR speckle filtering", Proc. of IGARSS'96, Lincoln (USA), Vol.3, pp.1574-1576, 27-31 May 1996.
95 - R. Touzi and A. Lopès, 1994: "The principle of speckle filtering in polarimetric SAR imagery", IEEE Trans. on GRS, Vol.32, nr.5, pp.1110-1114, 1994.
96 - G. De Grandi, E. Nezry and G. Kattenborn, 1993: "Experimental characterization of spatial statistics in polarimetric multifrequency airborne SAR images", Proc. of IGARSS'93, Tokyo, (Japan), Vol.2, pp.806-812, 18-21 August 1993.
97 - G. De Grandi, J.S. Lee, D. Schuler and E. Nezry, 2003: "Texture and speckle statistics in polarimetric SAR synthesized images", IEEE Trans. on GRS, Vol.41, n°9, pp.2070-2088, September 2003.
98 - E. Nezry and F. Yakam-Simen, 1999: "On the preservation of polarimetric signatures and polarimetric texture signatures by fully polarimetric MAP speckle filters", Proc. of IGARSS'99, Vol.3, 3 p., Hamburg (Germany), 28 June-2 July 1999.