#### M-Estimators of Roughness and Scale for -Modelled SAR Imagery

EURASIP Journal on Applied Signal Processing
M-Estimators of Roughness and Scale for 0 -Modelled SAR Imagery GA
Oscar H. Bustos 0 1 2
María Magdalena Lucini 0 1 2
Alejandro C. Frery
0 Centro de Informática, Universidade Federal de Pernambuco , CP 7851, 50732-970 Recife - PE , Brasil
1 Facultad de Matemática, Astronomía y Física, Universidad Nacional de Córdoba , Ciudad Universitaria, 5000 Córdoba , Argentina
2 Facultad de Matemática, Astronomía y Física, Universidad Nacional de Córdoba , Ciudad Universitaria, 5000 Córdoba , Argentina
The GA0 distribution is assumed as the universal model for multilook amplitude SAR imagery data under the multiplicative model. This distribution has two unknown parameters related to the roughness and the scale of the signal, that can be used in image analysis and processing. It can be seen that maximum likelihood and moment estimators for its parameters can be influenced by small percentages of “outliers”; hence, it is of outmost importance to find robust estimators for these parameters. One of the best-known classes of robust techniques is that of M-estimators, which are an extension of the maximum likelihood estimation method. In this work we derive the M-estimators for the parameters of the GA0 distribution, and compare them with maximum likelihood estimators with a Monte-Carlo experience. It is checked that this robust technique is superior to the classical approach under the presence of corner reflectors, a common source of contamination in SAR images. Numerical issues are addressed, and a practical example is provided.
and phrases; inference; likelihood; M-estimators; Monte-Carlo method; multiplicative model; speckle; synthetic aperture radar; robustness
1. INTRODUCTION
The statistical modeling of synthetic aperture radar (SAR)
imagery has provided some of the best tools for the processing
and understanding of this kind of data. Among the statistical
approaches the most successful is the multiplicative model.
This model offers a set of distributions that, with a few
parameters, are able to characterize most of SAR data. This
model is presented, for instance, in [1], and extended in [2].
This extension is a general and tractable set of
distributions within the multiplicative model, used to describe every
kind of SAR return. It was then called a universal model, and
its properties are studied in [3, 4, 5].
In this paper, the problem of estimating the parameters of
this extension, namely of the GA0 distribution, is studied for
the single look (the noisiest) situation. Two typical
estimation situations arise in image processing and analysis, namely
large and small samples, being the latter considered in this
work. The small samples problem arises in, for instance,
image filtering where with a few observations within a window
a new value is computed.
Statistical inference with small samples is subjected to
many problems, mainly bias, large variance, and sensitivity
to deviations from the hypothesized model. This last issue is
also a problem when dealing with large samples.
Robustness is a quite desirable property for estimators,
since it allows their use even in situations where the quality of
the input data is not perfect [6, 7, 8]. Most image processing
and analysis procedures (filtering, classification,
segmentation, etc.) use spatial data. Even experienced users are unable
to guarantee that all the input data are free of spurious values
and/or structures.
A situation where this occurs is in the presence of corner
reflectors. These devices, which are essential for data
calibration, produce a return which is quite higher than the rest
of the image. If data from a corner reflector enter a
nonrobust estimation procedure, the results may be completely
unreliable. Figure 1 shows man-made corner reflectors, the
white spot regularly placed horizontally along the middle of
the image, possibly buildings, on a SAR image where crops
predominate.
This paper presents classical estimators, namely those
based on sample moments and the maximum likelihood ones,
and derives robust M-estimators for the single look GA0 model.
Once these estimators are derived, several situations are
compared by means of a Monte-Carlo experience. These
estimators are then applied to SAR imagery, and it is shown that
the robust approach exhibits a good performance even in the
presence of corner reflectors.
M-estimators have been mainly used for symmetric data,
being [7, 8] two very relevant exceptions. In this application
they are computed, implemented, and assessed for speckled
imagery that, as will be seen in the next section, follow laws
whose densities can be highly asymmetric. Since speckle noise
also appears in ultrasound B-scan, sonar and laser images,
the procedures here presented have potential application in
all these techniques.
2. NOTATION AND MODEL DEFINITION
In the following, 1A will denote the indicator function of the
set A, that is,
1A(x) =
1 if x ∈ A,
0 else.
(
1
)
The single-look GA0(α, γ) distribution is characterized by
the following probability density function:
−2αx
f x, (α, γ) = γ(1 + (x2/γ))1−α 1(0,+∞)(x), −α, γ, x > 0.
(
2
)
Multilook, intensity, and complex versions of this distribution
can be seen in [2]. The polarimetric (multivariate) extension
of the G0 distribution is presented in [9].
The parameter α in (
2
) describes the roughness of the
target, being small values (say α < −15) usually associated
to homogeneous targets, like pasture, values ranging in the
[−15, −5] interval usually observed in heterogeneous clutter,
like forests, and big values (−5 < α < 0 for instance)
commonly seen when extremely heterogeneous areas, as urban
spots, are imaged. The γ parameter is related to the scale, in
the sense that if ZA is a GA0(α, 1)-distributed random
vari(
3
)
(
4
)
(
6
)
(
7
)
able then Z = √γZ obeys a GA0(α, γ) law and, therefore, it is
related to the brightness of the scene. Both parameters are
essential when characterizing targets, and in image processing
techniques involving statistical modeling.
In the following it will be assumed that α < 0 and γ > 0.
Calling
∂ ln f x, (α, γ) ,
s1 x, (α, γ) = ∂α
∂ ln f x, (α, γ) ,
s2 x, (α, γ) = ∂γ
the score function of the GA0(α, γ) distribution is given by
s x, (α, γ) =
s1 x, (α, γ)
s2 x, (α, γ)
This will be used to compute an estimator based on order
statistics. In [5] a relation between a more general version of
the GA0(α, γ) law and Snedekor’s F distribution is shown to
allow writing (
5
) as a function of the cumulative distribution
function of the latter.
The moments of a GA0(α, γ) distribution are given by
E Xk =
∞
γk/2 kΓ (k/2)Γ (−α − k/2)
2Γ (−α)
k
if − α > 2 ,
else.
This distribution can be derived as the square root of
the ratio of two independent random variables: one
obeying a unitary-mean exponential distribution (which conveys
the speckle noise in one look) and one following a Γ (α, γ)
distribution, related to the unobserved ground truth, the
backscatter. Densities of the GA0(α, γ) distribution are shown
in Figure 2, for α ∈ {−1, −3, −9} and the scale parameter
γ = γα adjusted to have unitary mean, that is, using (
6
) and
k = 1, one derives
4 Γ (−α)
γα = π Γ (−α − 1/2)
Following Barndorff-Nielsen and Blaesild [10], it is
interesting to see these densities as log probability functions,
particularly because the GA0 law is closely related to the class
of hyperbolic distributions [11]. Figure 3 shows the densities
0
of the GA(−2.5, 7.0686/π ) and the Gaussian distribution
N (1, 4(1.1781 − π /4)/π ) distributions in semilogarithmic
scale, along with their mean value. The parameters were
(
9
)
(
10
)
1.0
0.8
chosen, using (
6
), so that these distributions have equal mean
and variance. The different decays of their tails in the
logarithmic plot are evident: the former behaves logarithmically,
while the latter decays quadratically. This behavior ensures
the ability of the GA0 distribution to model data with extreme
variability.
3. INFERENCE FOR THE GA0 MODEL
It can be seen that the maximum likelihood estimators of α
and γ based in X1, . . . , XN , represented here by αˆML,N and
γˆML,N , respectively, are given by
1 N
αˆML,N = − N k=1 ln 1 +
Xk
Assume α < −1/2 in order to have random variables
with finite mean. For each j = 1/2, 1, define the jth order
moment as mj,N = N−1 kN=1 Xkj . Using (
6
) it can be seen
that the half and first order moments estimators of α and γ
based on X1, . . . , XN , denoted here by αˆMOM,N and by γˆMOM,N ,
respectively, are given by the solution of
m1/2,N = γˆM1/O4M,N
m1,N = γˆM1/O2M,N
Γ − αˆMOM,N − 1/4 Γ (1/4)
4Γ − αˆMOM,N
√π Γ − αˆMOM,N − 1/2
2Γ − αˆMOM,N
Moment estimation is extensively used in remote
sensing applications [2], mainly because it is inexpensive from
the computational point of view and analytically tractable in
most situations. In the presence of corner reflectors, though,
severe numerical instabilities were observed.
A comparison among different estimators for
roughness parameter of the GA0 distribution was performed in [4]
through a Monte-Carlo experience. No contamination was
considered, and the αˆML estimator was the best one in
almost all cases with respect to the mean square error (mse)
and the closeness to the true value (bias) criteria, assuming
γ is known. No analytical comparison is possible, due to the
complexity of the estimators.
In previous works (see [12], for instance) it is shown that
moment and maximum likelihood estimators have many
optimal properties when the observations, x1, . . . , xN are the
outcome from independent, identically distributed random
variables, with common density f (·, (α, γ)). Among these
properties, maximum likelihood estimators are
asymptotically unbiased, that is, if X1, . . . , XN is a sequence of
independent, identically distributed random variables with
common density f (·, (α, γ)) then, when N → ∞, it holds that
αˆMOM,N → α, γˆMOM,N → γ, αˆML,N → α, and γˆML,N → γ.
Nevertheless, good performance is neither warranted with
finite samples nor if the sample does not obey precisely the
hypothesis. A common departure from classical hypothesis is
the contamination by a percentage of “outliers,” that is, when
some of the observed data come from a different distribution.
“Corner reflectors” can be considered outliers in SAR
imagery. These are physical artifacts in the sensed area that
return most of the power they receive. The image in these areas
is dominated by the biggest possible values admitted by the
storage characteristics, and their effect is typically limited to a
few pixels. Corner reflectors are either placed on purpose, for
image calibration, due to man-made objects, such as highly
reflective urban areas, or the result of double-bounce
reflection [1].
Since the routines that compute both ML- and
Mestimators require initial guesses, the unstable (and, therefore,
unacceptable) behavior of moment estimators demanded the
use of another technique. A procedure based on analogy using
order statistics and moments [13], defined in the following
for α < −1, was used.
If X is GA0(α, γ) distributed the median of the scaled
random variable Y = X/E(X) can be computed using (
5
)
Q2 = √π
This scaled median does not depend on γ, so αˆ can be
estimated using the sample median q2(y), where y = (xi/x¯)i,
using standard numerical tools. An estimate of γ using the
first-order moment can then be computed as
4 Γ (−αˆ)
γˆ = π m12,N Γ (−αˆ − 1/2)
This estimator of α derived from (
11
) and the one computed
for γ through (
12
) will be called mixed estimator for (α, γ),
and denoted (αˆMIX, γˆMIX).
In the following, ML- and M-estimators will be assessed
in two situations: the pure model, when no contamination is
present, and cases where outliers simulating a corner reflector
enter the data.
The contamination model here considered is defined
as a sequence of independent, identically distributed
random variables X1, . . . , XN with common distribution
function F (·, (α, γ), , z) given by
F x, (α, γ), , z = (1 − )F x, (α, γ) + δz(x),
(
13
)
where δz(x) = 1[z,+∞)(x), with z a “very big” value with
respect to most of the values typically assumed by a random
variable with distribution GA0(α, γ). Equation (
13
) describes
a contamination that occurs at random with probability ,
that is, in average N out of N samples will be “very big”
values (corner reflectors), while N − N samples will obey
0
the GA(α, γ) distribution (the background). The flexibility
of the GA0(α, γ) distribution will allow the modeling of any
kind of background, namely crops, forests, and urban areas.
The contamination value z will be chosen as a factor of the
mean value of the underlying distribution GA0(α, γ).
Other contamination hypothesis may include different
distributions for the departure of the model, and/or spatial
dependence among observations.
In order to quantitatively assess the sensitivity of
MLestimators in our case, using the strong law of large numbers
on (
8
) and assuming γ = 1, it is immediate that
lim αˆML,N (z) = − EF(·,(α,γ),ε,z) ln 1 + X2
N→∞
1
= − ε ln(1 + z2) − (1 − ε)(1/α)
−1
Figure 4 shows, for α = −5, this limit as a function of ε,
the proportion of outliers in the sample, for several values of
the outlier z. Three values of z are shown for ε ∈ [0, 0.1].
It can be seen that a sample that suffers from a small
contamination of ε = 0.02 leads to the wrong conclusion that
extremely heterogeneous targets are under observation, since
the estimated value will be around −3 whilst the true value
is −5 (corresponding to an heterogeneous target). The
bigger the contamination the worse this behavior. This influence
on the ML-estimator is noticed regardless the value of z in
the chosen range. This result justifies a careful analysis of the
performance of estimators for a single representative value of
the contamination z, and for several values of ε, as presented
in Section 5.
4.
M-ESTIMATORS
These estimators offer a useful alternative when small
proportions of values may be far from the bulk of data. A good survey
on these estimators and their use in practical situations can
be seen in [7]. The difficulty in our case arises due to the fact
that the underlying distribution is asymmetric. Inspired in
the work by Marazzi and Ruffieux [8], that robustified
maximum likelihood estimators for the parameters of the gamma
distribution with success, we will compute M-estimators for
the parameters of the GA0(α, γ) distribution.
One has to devise ways to prevent a large influence of
outliers in the likelihood equations. M-estimators propose
the use of ψ functions that truncate the score of influential
observations in the likelihood equations.
Let b1 and b2 be two real positive numbers. We will call
M-estimators of parameters α and γ based on the sample
X1, . . . , XN , the statistics αˆM,N , and γˆM,N such that
= 0,
= 0,
(15)
N
k=1
N
k=1
(
14
)
ψb1 s1 Xk, αˆM,N , γˆM,N
− c1 αˆM,N , γˆM,N , b1
ψb2 (s2(Xk, αˆM,N , γˆM,N
− c2 αˆM,N , γˆM,N , b2
where s1 and s2 are given in (
4
), and ψb1 and ψb2 are
ψbi (x) =
bi
−bi if x ≤ bi,
x
if − bi < x < bi,
if x ≥ bi,
(16)
with x ∈ R and i = 1, 2. Note that making these functions
the identity and c1 = c2 ≡ 0, equations (15) reduce to the
familiar system of likelihood equations. Figure 5 illustrates
the ψ6(x) function, where its effect on the data becomes
clear: it truncates the score of those values above and below
its threshold. In this manner, this function prevents abnormal
(too small and too big) data having excessive influence on the
estimates.
The functions ci : Θ × (0, +∞) → R in (15) are defined in
such way that (αˆM,N , γˆM,N )N is a sequence of asymptotically
unbiased estimators of (α, γ), that is,
ψb1 s1 x, (α, γ) − c1 α, γ, b1 f x, (α, γ) dx = 0,
ψb2 s2(x, (α, γ) − c2 α, γ, b2 f x, (α, γ) dx = 0,
(17)
for any (α, γ) ∈ Θ, b1 > 0 and b2 > 0. The computation of
these functions is presented in Appendices A and B.
10
5
0
−5
−10
−10
−5
5
10
The values b1 and b2, called “tuning parameters” are
chosen in such a way that the efficiency of the M-estimators is not
unacceptably poor with respect to the maximum likelihood
ones. Thus, we will choose b1 and b2 such that
V αˆML
V αˆM
e1,
V γˆML
V γˆM
e2,
(18)
where V denotes variance, with 0.9 ≤ ei ≤ 0.975 for
example. This restriction imposes that the variance of the robust
estimators will not surpass those of the maximum likelihood
ones in more than a certain factor. The computation of these
parameters is done in Appendix C.
In order to assess the finite sample behavior of the
proposed estimators consider the situation where N observations
from independent identically distributed GA0(α, γ) random
5
10
Contamination Value
15
variables are contaminated by L observations taking the value
z. The contaminated sample z = (z1, . . . , zN , z, . . . , z), where
the brackets denote the N + L outliers, will be used to
estimate α and γ. Depending on the true parameters, on the
contamination (the value z and the number of outliers L)
and on the sample size N, it is expected that any non-robust
estimator will be mislead by the departure from the
hypothesized model. Figure 6 shows the maximum likelihood and
M-estimators for the situation where N = 77 pure
observations and L = 4 equal contaminated values are used. It can
be seen that the former suffers from both the departure from
the true model and from numerical instabilities (the peak
around the value 10). It is noticeable that the robust
estimator remains closer to the true value (the dash-dot line) than
the other in the presence of contamination. The bigger the
value of the contamination the further maximum likelihood
will go from the true value, while the M-estimator will render
the same estimate. The same behavior was observed for the
estimators of γ. The values employed for N, L, α and z are
representative of the kind of problem at hand: filtering
images over a heterogeneous area (α = −3) with 9 × 9 windows
(N + L = 81) where some observations (L = 4) are due to a
corner reflector (z varying and taking large values).
5. RESULTS
A Monte-Carlo experience was performed to compare
maximum likelihood (ML) and M-estimators (M), since analytical
comparisons are unfeasible. GA0 deviates can be obtained in
two ways, namely (a) by the inversion of the cumulative
distribution function, given in (
5
), or (b) by the constructive
definition of this distribution, that is, returning z = y/x
where y is a sample of the unitary mean exponential
distribution and x is a sample of the Γ (−α, γ) distribution, X
independent of Y . Alternative (b) was chosen in this
experience. Since the parameter γ is a scale factor, all the studied
situations will use the parameter that makes the expected
value 1, that is, γ = γα = 4(Γ (−α)/Γ (−α − 1/2))2/π .
We studied four models that will be identified as Models
1, 2, 3, and 4: Model 1 (the “ pure” model) consisting of
samples of the F (·, (α, γα)) distribution; Model 2 consisting of
samples taken from the F (·, (α, γα), , z) distribution with
= 0.01 and z = 15 , Model 3 is a collection of samples
obtained from the F (·, (α, γα), , z) distribution with = 0.05
and z = 15, while Model 4 is a collection of samples obtained
from the F (·, (α, γα), , z) distribution with = 0.10 and
z = 15. In this manner, the pure situation and three levels
of contamination are examined. The contamination value z
is held fixed, as discussed in Section 3, as fifteen times the
expected value of the underlying model which is big enough
to describe corner reflectors.
In each of these models we deal with three situations,
namely with: Situation 1 (α = −1.0), Situation 2 (α = −3),
and Situation 3 (α = −6). These situations cover
representative targets in practice: extreme heterogeneity, heterogeneity
and homogeneity.
In order to determine the number of replications an
empirical precision criterion was used. One hundred replications
(1 ≤ r ≤ R = 100) with samples of size N = 9 × 9 = 81
were performed. One estimator for α is computed for each
replication r , say αˆ(r ). The final number of replications
R = R + , ≥ 1 integer, is the smallest number of
samples that allows forming an asymptotic confidence interval
at the 95% level on the mean of the R estimators, that is,
2 × 1.96sR/√R < 0.5 where sR is the standard deviation of
the sample {αˆ(
1
), . . . , αˆ(R)}. This is computed for every
estimator for α and every estimator for γ in every model and
every situation. The biggest number of replications R (the
worst case) is then used. This procedure led to R = 101 in
every model and every situation.
The following tables show the results obtained in this
study for each situation: for every model the sample mean of
the estimator over the R replications Eˆ(ζ) = R−1 rR=1 ζ(r ),
the mean square error mse = Eˆ(ζ −θ)2 = Vˆ(ζ)+(Eˆ(ζ)−θ)2,
and the absolute relative bias B(ζ) = θ−1|Eˆ(ζ) − θ|, where
θ is the true value of the parameter and ζ is the estimator
under study.
In all the tables it can be seen that M-estimators are almost
as good as ML-estimators when there is no contamination,
while in the presence of outliers M-estimators exhibit smaller
absolute relative bias and smaller mean squared error. As α is
smaller (the observed region is more homogeneous),
estimation becomes less reliable. This is probably due to the shape
of the likelihood function, that becomes flat and, therefore, is
hard to find its maximum location. The bigger the proportion
of contamination the worse the behavior of both estimators,
but M-estimators remain consistently closer to the true value
than ML-estimators, as expected.
As an application, the image shown in Figure 7 is
analyzed. Windows of size 9 × 9 are used to estimate the
background parameters over a trajectory that spans from the
upper left corner to the lower right corner. This trajectory passes
through the three corner reflectors of the image, which are
clearly seen as white spots.
Figure 8 shows the result of computing both ML (solid
line) and M-estimators (dashes) over 9 × 9 windows in every
position of the trajectory in Figure 7. The dotted line shows
the background value of α, and the corner reflectors are
located in the positions marked with dash-dots. It is clear that
the M-estimator is consistently closer than the ML-estimator
to the true value at positions where there is contamination.
Besides that, it is always possible to compute M-estimators,
whereas ML-estimator fails to converge in 5 out of 45
positions (the gaps in the solid line). Similar results were observed
estimating the scale parameter γ.
6. CONCLUSIONS AND FUTURE WORK
There are numerical problems with the computation of these
estimators for certain parameters. For small samples there is
often no solution, being this situation more critical for small
values of α, that is, for homogeneous areas. This issue will be
further investigated in future works.
As it can be seen from the tables, in the pure model
(Model 1), the behavior of the robust estimators is as good
as the maximum likelihood ones, as it is expected. In the
contaminated situations (Models 2 and 3), we can see that
estimation by moments or maximum likelihood methods is
poor.
It is relevant to notice that under a very small
contamination or a very small deviation from the model the behavior
of the classical estimators is not reliable.
This work will continue computing the M-estimators for
the multilook case, that is, for the n > 1 situation and for
polarimetric (multivariate) data. An alternative approach
using alpha-stable distributions [14] is also possible, but at the
expense of loosing the interpretability of the parameters that
stem from the multiplicative model.
APPENDICES
A. COMPUTATION OF c1
Let α < 0. The problem consists of finding c1(α, 1) that
satisfies I1(α, c1, b) = 0, where I1(α, c, b) = 0∞ ψb((1/α) +
ln(1 + u) − c)((−α)/(1 + u)1−α)du is the equation that has
0
to be solved in order to define the M-estimators for the GA
distribution (equations (15)).
Lemma A.1. Let α < 0, then
(
1
) If exp(2αb) > bα + 1 then there exists c0, with −b +
α−1 < c0 < b + α−1 that satisfies I1(α, c0, b) = 0. In this case
c0 is αc0 = exp(αb + αc0 − 1).
(
2
) If exp(2αb) < bα + 1 then there exists c0, with c0 ≥
b + α−1 that satisfies I1(α, c0, b) = 0, and it is given by c0 =
α−1(ln(−α) + ln(b) − ln(e−αb − eαb) + 1).
Proof. The function G(c) = I1(α, c, b) is a continuous
monotone decreasing function, then
lim G(c) = b,
c→−∞
lim G(c) = −b.
c→∞
(A.1)
Therefore, there exits c0 that satisfies G(c0) = I1(α, c0, b) =
0. Denoting A = A(α, c) = exp(−b − (1/α) + c) − 1 and
B = B(α, c) = exp(b − α−1 + c) − 1, then A < B and one can
write
ψb
1
α + ln(1 + u) − c
−b if u < A,
1
= α + ln(1 + u) − c if A < u < B,
b if u > B.
If c < −b + α−1, then B < 0, and then I1(α, c, b) = b.
So, there is no c0 < −b + α−1 satisfying I1(α, c0, b) = 0. If
−b + α−1 ≤ c < b + α−1, then A < 0 < B and, therefore,
I1(α, c, b) =
exp(αb + αc − 1)
α
− c.
As I1(α, −b + α−1, b) = b > 0, then c0 > −b + α−1, and
If I1(α, −b + (1/α), b) < 0, then there exits −b + α−1 <
c0 < b + α−1 such that I1(α, c0, b) = 0. This occurs when
exp(2αb) > bα + 1 and c0 must satisfy
Let α < 0 and γ > 0. As for the computation of c1, the
problem of finding the function c2 consists of finding c2 =
c2(α, γ, b2) that satisfies
where I2(α, γ, c, b2) = 0∞ ψγb2 (−α − (1 − α)/(1 + u) −
c)(−α)/(1 + u)1−αdu. In this situation there are more cases
than in the computation of c1 due to the presence of the
parameter γ.
Lemma B.1. Let α < 0 and γ > 0. In order to find the function
c2 that satisfies (B.1) one has to consider the following cases:
(
1
) if 1 − α − 2γb2 > 0 then
(a) if −(1 − α − 2γb2)1−α ≤ (γb2 − 1)(1 − α)1−α
then there exists c2 such that c2 ≤ (γb2 −
1), and it is given by the solution of c2 =
−((−γb2 − α − c2)/(1 − α))1−α,
(b) if −(1 − α − 2γb2)1−α > (γb2 − 1)(1 − α)1−α
then c2 is given by
(i) the solution of γb2(1 − α)1−α = (γb2 −
α − c2)1−α − (−γb2 − α − c2)1−α if
(2γb2/(1 − α))1−α − γb2 < 0, or
(ii) c2 = γb2 − α − (γb2)1/(1−α)(1 − α) if
(2γb2/(1 − α))1−α − γb2 ≥ 0,
(
2
) if 1 − α − 2γb2 ≤ 0 then c2 satisfying (B.1) is given by
(a) the solution of c2 = −((−γb2 − α − c2)/
(1 − α))1−α, if γb2 < −α, or
(b) c2 = 0, if γb2 > −α and γb2 > 1, or
(c) c2 = γb2 − α − (γb2)(1−α)−1 (1 − α) if γb2 > −α
and γb2 ≤ 1.
Proof. As the function G(c) = I2(α, γ, c, b) is a continuous
monotone decreasing function, and
lim G(c) = b2,
c→−∞
lim G(c) = −b2,
c→∞
(B.2)
then there exists c2 that satisfies G(c2) = I2(α, γ, c2, b2) = 0.
− b1.
(A.6)
Case 1. Suppose that 1−α−2γb2 > 0, so γb2−1 < −α−γb2.
If c ≤ γb2 − 1 and A = (1 − γb2 + c)/(γb2 − α − c), then
c < −α − γb2, and
(B.3)
(B.4)
(B.6)
1 − α
ψγb2 − α − 1 + u − c
−γb2
γb2
= −α − 11 +− uα − c if A ≤ u ≤ B,
if u < A,
if u > B,
where B = (1 + γb2 + c)/(−γb2 − α − c); then
I2 α, γ, c, b2 = −
−γb2 − α − c
1 − α
1−α
− c.
(B.5)
Taking c = γb2 − 1 one has that I2(α, γ, γb2 − 1, b2) =
−(1 − 2γb2/(1 − α))1−α − (γb2 − 1), therefore, if −(1 − α −
2γb2)1−α ≤ (γb2 − 1)(1 − α)1−α one can say that there
exists c2(α, γ, b2) satisfying (B.1), given by the solution of
c2 = −((−γb2 − α − c2)/(1 − α))1−α.
For the other situation, if −(1 − α − 2γb2)1−α > (γb2 −
1)(1 − α)1−α one has to look for c2 in the interval (γb2 −
1, −α − γb2]. In this case,
1 − α
ψγb2 − α − 1 + u − c
1 − α
−γb2 if − α − 1 + u − c < −γb2,
1 − α 1 − α
= −α − 1 + u − c if − γb2 < −α − 1 + u − c < γb2,
γb2 if − α − 11 +− uα − c > γb2.
I2 α, γ, c, b2 = −γb2 +
γb2 − α − c2
1 − α
1−α
−
−γb2 − α − c2
1 − α
1−α
.
Replacing c in the above equation by −α − γb2 one
obtains I2(α, γ, −α − γb2, b2) = (2γb2/(1 − α))1−α − γb2.
Then, if (2γb2/(1 − α))1−α − γb2 ≤ 0, so there exists
γb2 − 1 < c2 ≤ −α − γb2 satisfying (B.1) and it is given
by the solution c2 that makes equation (B.6) zero.
But, when (2γb2/(1 − α))1−α − γb2 > 0 one has to seek
for the solution of (B.1) in the interval (−α − γb2, +∞). In
this case
γb2 − α − c
1 − α
1−α
Hence, c2 = γb2 − α − (γb2)(1−α)−1 (1 − α).
Case 2. Consider the situation 1 − α − 2γb2 ≤ 0, so that
γb2 − 1 ≥ −α − γb2. If c ≤ −α − γb2 then one has that
−((−γb2 − α − c2)/(1 − α))1−α.
I2(α, γ, c, b2) = −((−γb2 − α − c)/(1 − α))1−α − c. Taking
c = −α − γb2, then I2(α, γ, −α − γb2, b2) = α + γb2. As
I2(α, γ, −α − γb2, b2) < 0 if γb2 < −α, one can say that
c2 satisfying equation (B.1) is given by the solution of c2 =
But, if γb2 ≥ −α, one can take −γb2 − α < c ≤ γb2 −
1. In this case one has that I2(α, γ, c, b2) = −c, therefore
I2(α, γ, γb2 − 1, b2) = 1 − γb2. Thus, if 1 < γb2 one can say
that the solution of equation (B.1) is c2 = 0.
If γb2 ≤ 1, then one takes c > γb2 − 1, leading to
Hence, c2 = γb2 − α − (γb2)1/(1−α)(1 − α) is the solution
of (B.1).
where
Therefore, replacing matrix (C.4) and the inverse of
matrix (C.5) in equation (C.1), we have
V (α, γ)ML
=
The definition of efficient M-estimators requires finding
tuning parameters b1 and b2 ensuring that the relative
asymptotic variance (with respect to the variance of maximum
likelihood estimators) satisfies V (αˆML)/V (αˆM) e1 and
V (γˆML)/V (γˆM) e2 with, for instance, 0.9 ≤ ei ≤ 0.975.
V (αˆML), V (γˆML) and V (αˆM), V (γˆM) are the diagonal
elements of the asymptotic covariance matrix of the
MLand M-estimators of parameters α and γ, here denoted by
V ((αˆ, γˆ)ML) and V ((αˆ, γˆ)M), respectively.
C.1. Computation of V ((α, γ)ML)
We denote
V (α, γ)ML = M s, F (α, γ) −1J(α, γ)M s, F (α, γ) −T ,
where F (α, γ) is the cumulative distribution function of
0
a random variable having GA(α, γ, 1) distribution, s =
s(x, (α, γ)) is the vector of score functions,
J(α, γ) =
s x, (α, γ) s x, (α, γ) T dF x, (α, γ)
(C.2)
is the Fisher information matrix, and
∂
∂(α, γ)
s x, (α, γ) dF x, (α, γ) .
M s, F (α, γ) = −
It can be seen that
J(α, γ) =
α−2
γ(1 − α) −1
γ(1 − α) −1
−α (2 − α)γ2 −1 ,
and that
M s, F (α, γ) = J(α, γ)
=
α−2
γ(1 − α) −1
γ(1 − α) −1
−α (2 − α)γ2 −1 .
.
(C.6)
(C.7)
,
C.2. Computation of V ((α, γ)M)
We define
V (α, γ)M = M Ψ , F (α, γ) −1Q(α, γ)M Ψ , F (α, γ) −T ,
(C.1)
ψb1 s1(x; (α, γ) − c1 α, γ, b1)
ψb2 s2(x; (α, γ) − c2 α, γ, b2
Ψ x, (α, γ) dF x, (α, γ) ,
Ψ = Ψ x; (α, γ) =
∂
M Ψ , F (α, γ) = − ∂α
Q(α, γ) =
Ψ x, (α, γ) Ψ x, (α, γ) T dF x, (α, γ) .
(C.8)
The matrices M and Q are computed taking into account
all the cases that define the functions ci, i = 1, 2, and solving
explicitly the integrals.
ACKNOWLEDGEMENT
The authors are grateful to Conicor, SeCyT, CNPq, and Vitae
for the partial support of this research.
(C.3)
(C.4)
(C.5)
Alejandro C. Frery received the degree in
Electronic Engineering from the
Universidad de Mendoza, Argentina, in 1987. His
MSc degree was in Applied Mathematics
(Statistics) from the Instituto de Matemática
Pura e Aplicada (IMPA, Rio de Janeiro,
Brazil, 1990) and his Ph.D. was in Applied
Computing from the Instituto Nacional de
Pesquisas Espaciais (INPE, São José dos
Campos, Brazil, 1994). He is currently a lecturer the Computer
Science Institute of the Universidade Federal de Pernambuco, Brazil.
His research areas include stochastic models for image formation
and analysis, as well as stochastic simulation.
[1] C. Oliver and S. Quegan , Understanding Synthetic Aperture Radar Images, Artech House, Boston, 1998 .
[2] A. C. Frery , H. - J. Müller , C. C. F. Yanasse , and S. J. S. Sant'Anna , “ A model for extremely heterogeneous clutter , ” IEEE Transactions on Geoscience and Remote Sensing , vol. 35 , pp. 648 - 659 , 1997 .
[3] J. Jacobo-Berlles , M. Mejail , and A. C. Frery , “ The GA0 distribution as the true model for SAR images ,” in Simpósio Brasileiro de Computação Gráfica e Processamento de Imagens, SBC, IEEE, Campinas, SP , Brazil, 1999 .
[4] M. Mejail , J. Jacobo-Berlles , A. C. Frery , and O. H. Bustos , “ Parametric roughness estimation in amplitude SAR images under the multiplicative model ,” Revista de Teledetección, vol. 13 , pp. 37 - 49 , 2000 .
[5] M. E. Mejail , A. C. Frery , J. Jacobo-Berlles , and O. H. Bustos , “ Approximation of distributions for SAR images: Proposal, evaluation and practical consequences , ” Latin American Applied Research , vol. 31 , pp. 83 - 92 , 2001 .
[6] O. H. Bustos , “ Robust statistics in SAR image processing , ” in Proceedings of the First Latino-American Seminar on Radar Remote Sensing: Image Processing Techniques , pp. 81 - 89 , ESA, Paris, 1996 .
[7] F. R. Hampel , E. M. Ronchetti , P. J. Rousseeuw , and W. A. Stahel , Robust Statistics: The Approach Based on Influence Functions , Wiley, New York, 1986 .
[8] A. Marazzi and C. Ruffieux , “ Implementing M-estimators of the gamma distribution ,” in Robust Statistics, Data Analysis , and Computer Intensive Methods, R. Helmut, Ed., vol. 109 of Lecture Notes in Statistics, Springer-Verlag, Berlin, 1996 .
[9] A. C. Frery , A. H. Correia , C. D. Rennó , C. C. Freitas , J. JacoboBerlles, M. E. Mejail , and K. L. P. Vasconcellos , “ Models for synthetic aperture radar image analysis,” Resenhas (IME-USP) , vol. 4 , no. 1 , pp. 45 - 77 , 1999 .
[10] O. E. Barndorff-Nielsen and P. Blaesild , “ Hyperbolic distributions and ramifications: Contributions to theory and applications,” in Statistical distributions in scientific work, C. Taillie and B. A. Baldessari , Eds., pp. 19 - 44 , Reidel, Dordrecht, 1981 .
[11] A. C. Frery , C. C. F. Yanasse , and S. J. S. Sant'Anna , “ Alternative distributions for the multiplicative model in SAR images,” in Quantitative remote sensing for science and applications , IGARSS'95 Proc. , pp. 169 - 171 , IEEE, Florence, July 1995 .
[12] P. K. Sen and J. M. Singer , Large Sample Methods in Statistics: An Introduction with Applications , Chapman and Hall, London, 1993 .
[13] C. F. Manski , Analog Estimation Methods in Econometrics, vol. 39 of Monographs on Statistics and Applied Probability, Chapman & Hall, New York, 1988 , available at http://elsa.berkeley.edu/books/analog.html.
[14] C. L. Nikias and M. Shao , Signal Processing with Alpha-Stable Distributions and Applications , Wiley, New York, 1995 .