V
OLUME
20 N
OVEMBER
2003JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY
q 2003 American Meteorological Society 1433
Whitening in Range to Improve Weather Radar Spectral Moment Estimates.
Part I: Formulation and Simulation
S
EBASTIA
´
N
M. T
ORRES
Cooperative Institute for Mesoscale Meteorological Studies, University of Oklahoma, and National Severe Storms Laboratory,
Norman, Oklahoma
D
US
ˇ
AN
S. Z
RNIC
´
National Severe Storms Laboratory, Norman, Oklahoma
(Manuscript received 29 May 2002, in final form 14 March 2003)
ABSTRACT
A method for estimation of spectral moments on pulsed weather radars is presented. This scheme operates
on oversampled echoes in range; that is, samples of in-phase and quadrature-phase components are collected at
a rate several times larger than the reciprocal of the transmitted pulse length. The spectral moments are estimated
by suitably combining weighted averages of these oversampled signals in range with usual processing of samples
(spaced at the pulse repetition time) at a fixed range location. The weights in range are derived from a whitening
transformation; hence, the oversampled signals become uncorrelated and, consequently, the variance of the
estimates decreases significantly. Because the estimate errors are inversely proportional to the volume scanning
times, it follows that storms can be surveyed much faster than is possible with current processing methods, or
equivalently, for the current volume scanning time, accuracy of the estimates can be greatly improved. This
significant improvement is achievable at large signal-to-noise ratios.
1. Introduction
Doppler weather surveillance radars probe the at-
mosphere and retrieve spectral moments for each res-
olution volume in the surrounding space. In computing
these moments, it is customary to average signals from
many pulses to reduce the statistical uncertainty of the
estimates. With such processing, the variance reduction
of averaged estimates is inversely proportional to the
equivalent number of independent samples, which de-
pends on the correlation between samples
r
and the total
number of averaged samples (or pulses) M (Walker et
al. 1980). The number of samples available for aver-
aging is determined by the pulse repetition time T
s
and
the dwell time, which is usually controlled by the re-
quired azimuthal resolution. If averaging along sample
time is not enough to keep estimation errors below ac-
ceptable limits, it is possible to trade range resolution
for estimate accuracy by averaging a few samples along
range time. From an operational point of view, we are
faced with conflicting requirements. On the one hand,
large estimation errors restrict the applicability of
weather surveillance radars for precise quantification
Corresponding author address: Dr. Sebastia´n M. Torres, NSSL,
1313 Halley Circle, Norman, OK 73069.
E-mail: sebastian.torres@noaa.gov
and identification of weather phenomena. On the other
hand, the need for faster updates between volume scans
calls for faster antenna rotation rates, which limit the
number of samples available for each resolution volume.
As mentioned before, the number of samples is inversely
related to the variance of estimates.
Several solutions have been proposed to reduce spec-
tral moment errors in weather surveillance radars. In the
quest for finding better estimators of spectral moments,
Zrnic´ (1979) showed that maximum likelihood (ML)
estimators yield errors one order of magnitude less than
those obtained with conventional autocovariance meth-
ods. Later, Frehlich (1993) improved Zrnic´’s results and
derived simplified expressions to test new estimators
based on the ML approach. Due to the complexity of
ML estimators, researchers focused on ways to simplify
spectral moment estimators by assuming knowledge of
some of the underlying parameters of the weather signal.
Bamler (1991) computed the Cramer–Rao lower bound
(CRLB) for Doppler frequency estimates assuming both
the correlation (or spectrum) of samples and signal-to-
noise ratio (SNR) are known. Later, Chornoboy (1993)
obtained an optimal estimator for Doppler velocity that
is simpler than ML formulations, but again, the SNR
and the spectrum width were assumed to be known.
Summarizing, ML estimators provide better accuracy
compared with classic estimators (Zrnic´ 1979) and are
1434 V
OLUME
20JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY
only moderately complex if the spectrum width is
known a priori. However, this last assumption restricts
their applicability because the correlation coefficient of
the weather echo along sample time is not known and
must be estimated. That is, to estimate spectral mo-
ments, the joint distribution of signal power, mean
Doppler velocity, and spectrum width would need to be
calculated, which turns out to be computationally very
intensive. Dias and Leita˜o (2000) proposed ML esti-
mators that do not exhibit the characteristic computa-
tional burden. They derived a nonparametric method to
obtain ML estimates of spectral moments from weighted
sums of autocovariance estimates; the only assumption
is that the power spectral density of the underlying pro-
cess is bandlimited.
Schulz and Kostinski (1997) suggested that knowl-
edge of the correlation coefficient of weather-echo time
series could improve the variance of spectral moment
estimates. By means of the whitening transformation
(i.e., a linear operation that decorrelates echo samples),
they devised estimators that theoretically achieve the
CRLB. Further, Koivunen and Kostinski (1999) ex-
plored practical aspects of the whitening transformation
on estimation of signal power. A principal obstacle to
the application of the whitening transformation on the
time series is the need to know (or estimate) the cor-
relation coefficient of the signal. These authors suggest
ways to overcome this difficulty and call for an exper-
imental study to verify the technique. Using the whit-
ening approach, Frehlich (1999) investigated the per-
formance of ML estimators of spectral moments under
the assumption of a known spectrum width. Although
Frehlich concluded that the estimator derived by Schultz
and Kostinski cannot be applied to the case of finite
SNR, Kostinski and Koivunen (2000) showed that the
problem is not in the estimator but in the Gaussian as-
sumption for the Doppler spectrum. In their work, they
suggested simple numerical recipes intended to avoid
the so-called Gaussian anomaly.
Range oversampling to improve the estimates of per-
iodograms was explored by Urkowitz and Katz (1996).
They acknowledged that the periodograms for each
range location are correlated; they computed the equiv-
alent number of uncorrelated range samples and showed
that the variance reduction is not optimal. Strauch and
Frehlich (1998) considered oversampled signals in
range to estimate Doppler velocity within one pulse.
This approach fails because the phase shift within a
pulse is smaller than the uncertainty of the estimates;
the authors point out that simple averaging is not enough
to achieve the required equivalent number of indepen-
dent samples. Acquisition and processing of samples
over finer range scales was investigated also in the con-
text of pulse compression. Pulse compression can be
applied to increase the equivalent number of indepen-
dent samples by averaging high-resolution estimates in
range (Mudukutore et al. 1998). However, most ground-
based weather surveillance radars do not use pulse com-
pression due to the required larger transmission band-
widths.
Range oversampling and whitening can be used to
increase the equivalent number of independent samples
without increasing the transmission bandwidth. The
main advantage of this technique is that the whitening
transformation is derived from a known correlation
function. In a rather short but important study, Dias and
Leita˜o (1993) derived an iterative technique to obtain
ML estimates of spectral moments. They consider both
the time and space variables and assume a known cor-
relation in range. Implicit in their solution is whitening
of the signals in range; thus, it is likely the earliest
application of this technique to radar remote sensing.
More recently, Fjørtoft and Lope`s (2001) proposed a
method for estimating the reflectivity in synthetic ap-
erture radar (SAR) images with correlated samples (pix-
els). The method is based on a modified whitening trans-
formation that exhibits low computational complexity
and is suitable for oversampled data.
This paper describes an application of the whitening
transformation in range that increases the equivalent
number of independent samples while keeping the dwell
time constant with no significant degradation of the
range resolution. Obtaining more independent samples
reduces the estimate errors at the same antenna rotation
rate or speeds up volume scans while keeping the errors
at previous levels. Our work has been inspired by Schulz
and Kostinski (1997) but shares some common elements
with the work of Dias and Leita˜o (1993). In simulation
studies we consider a receiver with a large bandwidth
and a perfect rectangular pulse. We use autocovariance
processing in sample time and examine in detail the
effects of white noise. At the end we briefly describe
some other practical aspects and trade-offs.
2. Why whitening?
Current implementations of spectral moment esti-
mators use a simple method of averaging samples in
range at the expense of degradation in range resolution.
Simple averaging, however, does not yield the best per-
formance when the observations are correlated. Schulz
and Kostinski (1997) computed the variance bounds for
reflectivity estimates and demonstrated that they do not
depend on the correlation structure of the observations.
Hence, it can be inferred that it is not the correlation
between observations that limits the accuracy of a given
estimator but the way those observations are used to
compute the estimates. Therefore, it is reasonable to
think that knowledge of the correlation coefficient
r
(mT
s
) could be used to formulate estimators that attain
the CRLB. In their work, Schulz and Kostinski proposed
whitening the time series data along sample time to
produce uncorrelated samples. The main drawback for
this technique is that the whitening transformation de-
pends on other meteorological parameters such as the
spectrum width. Dias and Leita˜o (1993) derived man-
N
OVEMBER
2003 1435TORRES AND ZRNIC
´
F
IG
. 1. Depiction of sampling/oversampling in range and process-
ing of the signals. (a) Samples in range with spacing equal to the
pulse length
t
; standard processing to obtain correlation estimates is
indicated; (b) oversampling in range; (c) zoomed presentation of ov-
ersampled range locations where range samples to be whitened with
matrix
W
are indicated; (d) processing of whitened samples to obtain
estimates of correlations in range and average of these estimates in
range to reduce the statistical errors.
ageable approximate solutions to the ML estimators of
spectral moments for signals that are oversampled in
range. Their solution, although not immediately trans-
parent to readers, amounts to whitening the samples in
range and applying Fourier transforms for estimating
spectral moments.
We suggest combining the whitening transformation
of samples in range with autocovariance processing in
sample time and thus improving the spectral moment
estimates. The proposed processing increases the equiv-
alent number of independent samples in a simple manner
while the sacrifice in range resolution is minimal and
the transmission bandwidth is not broadened. While the
correlation of samples separated by T
s
needs to be es-
timated for each particular case (it depends on the me-
teorological conditions being observed), samples spaced
in range exhibit a correlation coefficient that can be
exactly computed a priori; the underlying assumption
here is that the mean echo power changes very little
over the averaging interval in range. We elaborate more
on this fundamental constraint later. By exactly knowing
the correlation coefficient, it is possible to apply the
whitening transformation without worrying about the
pitfalls originating from an estimated quantity. As a re-
sult, the equivalent number of independent samples be-
comes equal to the number of available samples, and
the variance reduction through averaging is maximized.
Maximization of the equivalent number of independent
samples leads to the following.
• For the same uncertainty as that obtained with cor-
related samples, faster scan rates are possible, as the
total number of samples M for a resolution volume is
determined by the pulse repetition time and the dwell
time. Rapid acquisition of volumetric radar data has
significant scientific and practical ramifications. For
example, observations at minute intervals are required
to understand the details of vortex formation and de-
mise near the ground. Even faster rates of volumetric
data are required to determine the presence of trans-
verse winds. Fast update rates would also yield more
timely warnings of impending severe weather phe-
nomena such as tornadoes and strong winds.
• For the same scanning rates, lower uncertainties can
be obtained, making the use of polarimetric variables
feasible for accurate rainfall estimation and hydro-
meteor identification.
With the advent of digital receivers (Brunkow 1999),
oversampling is indeed feasible (Ivic´ 2001). Therefore,
it is possible to maintain the same current radar capa-
bilities (as with a digital matched filter in classical pro-
cessing) while adding, in parallel, a set of more reliable
estimates obtained from whitened oversampled range
data.
3. The whitening transformation
The procedure (as depicted in Fig. 1) begins with
oversampling in range so that there are L samples during
the pulse duration
t
(i.e., oversampling by a factor of
L). Assume that a range of depth c
t
(where c is the
speed of light) is uniformly filled with scatterers, which
is a common occurrence for relatively short pulses. For
convenience, the contribution from the resolution vol-
ume to the received sampled complex voltage V(nT
s
)
5 I(nT
s
) 1 jQ(nT
s
) at a fixed time delay nT
s
can be
decomposed into subcontributions s(l
t
o
, nT
s
) from L
contiguous elemental shells or ‘‘slabs,’’ each c
t
/2L
thick. Each of these is an equivalent scattering center.
For simplicity
t
o
and T
s
are dropped hereafter so the
indices l and n indicate range-time increments
t
o
(sam-
pling time) and sample-time increments T
s
(pulse rep-
etition time), respectively. The voltages s(l, n) are iden-
tically distributed, complex, Gaussian random variables,
where the real and imaginary parts have variances
s
2
,
and the average power of s(l, n)is 5 2
s
2
. A pulse
2
s
s
with an arbitrary envelope shape p(l) induces weighting
to the contributions from contiguous slabs such that the
composite voltage after synchronous detection is
V(l, n) 5 I(l, n) 1 jQ(l, n)
L21
5 s(l 1 i, n)p(L 2 1 2 i) , h(l), (1)
O
[]
i50
where the , denotes convolution and h(l) is the impulse
response of the receiver filter. The summation in (1) is
1436 V
OLUME
20JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY
a correct representation for the composite voltage V(l,
n) given that we consider ‘‘scattering centers,’’ and each
of these implicitly admits integration. Then, as a gen-
eralization of (4.40) in Doviak and Zrnic´ (1993), the
correlation of range samples is
(R)
2
R (l) 5
s
[p (l) , p*(2l)],
Vsmm
(2)
where the modified pulse envelope p
m
is given by p
m
(l)
5 p(l) , h(l). Hence, the correlation coefficient of range
samples is
(R)
r
V
(R)
R (l) p (l) , p*(2l)
Vmm
(R)
r
(l) 55 . (3)
V
L21
(R)
R (0)
V
2
p (l9)
O
m
l950
In practice, can be evaluated by attenuating the
(R)
r
V
transmitted pulse, injecting it directly into the receiver,
and oversampling the result to obtain the modified pulse
envelope p
m
. Introducing p
m
into (3) produces ; this
(R)
r
V
needs to be done only once for a given pulse shape and
receiver bandwidth.
The procedure for implementing the whitening trans-
formation follows Koivunen and Kostinski (1999) and
is listed here for completeness (and reader convenience).
Define the Toeplitz–Hermitian normalized correlation
matrix as
(R)
C
V
(R)(R)
1
r
(1) · · ·
r
(L 2 1)
VV
(R)(R)
r
(1)* 1 · · ·
r
(L 2 2)
VV
(R)
C
5 .
V
__5_
(R)(R)
r
(L21)*
r
(L 2 2)* · · · 1
VV
(4)
Because this matrix is positive semidefinite (Therrien
1992), it can be decomposed as
(R)
T
C
5
H
*
H
,
V
(5)
where the superscript T indicates matrix transpose and
* is the usual complex conjugation operation. Any
H
that satisfies (5) is called a square root of (Faddeev
(R)
C
V
and Faddeeva 1963) and is the inverse of a whitening
transformation matrix,
21
W
5
H
, (6)
which, if applied to the range samples, produces L un-
correlated random variables with identical power (Kay
1993).
In general, the decomposition of the correlation matrix
is not unique and many well-known methods can be ap-
plied to generate different whitening transformations. Two
such methods are the eigenvalue decomposition (Therrien
1992) and Cholesky decomposition, which is equivalent
to Gram–Schmidt orthogonalization (Therrien 1992; Pa-
poulis 1984). In the eigenvalue decomposition method,
is represented as 5
U
L
U
*
T
, where L is a diagonal
(R)(R)
CC
VV
matrix of eigenvalues of , and
U
is the unitary matrix
(R)
C
V
whose columns are the eigenvectors of . With this
(R)
C
V
decomposition, an expression of the same form as (5) can
be obtained as 5
H
*
H
T
5 (
U
*L
1/2
)*(
U
*L
1/2
)
T
, where
(R)
C
V
L
1/2
is a diagonal matrix with the square roots of the
eigenvalues on the diagonal. With this decomposition,
H
5
U
*L
1/2
, and
W
is obtained as
W
5
H
21
5 L
21/2
U
T
,
which is the Mahalanobis transformation (Tong 1995). In
the case of Cholesky (or triangular) decomposition, the
correlation matrix is factored as 5
TT
*
T
, where the
(R)
C
V
matrix
T
is a lower triangular matrix. Here, 5
H
*
H
T
(R)
C
V
5 (
T
*)*(
T
*)
T
, and
H
5
T
*; hence, the whitening matrix
W
5
H
21
5 (
T
*)
21
is also lower triangular. A possible
advantage of lower triangular
W
matrices is that whitening
can proceed in a pipeline manner; that is, computations
can start as soon as the first sample is taken and progress
through subsequent samples. Non-lower-triangular
W
ma-
trices require the presence of all data before a whitened
sample can be computed.
In the following sections, we denote with X(l, n) the
sequence of whitened samples obtained from V(l, n) for
a fixed sample time nT
s
as
L21
X(l, n) 5 wV(j, n); l 5 0,1,...,L 2 1, (7)
O
l,j
j50
where w
l,j
are the entries of the whitening matrix. Al-
ternatively, the previous equation can be written using
matrix notation as
X 5
W
V ,
nn
(8)
where V
n
5 [V(0, n), V(1, n),...,V(L 2 1, n)]
T
and
X
n
5 [X(0, n), X(1, n),...,X(L 2 1, n)]
T
. It is important
to note that regardless of the method used to decompose
, the whitening procedure is given by (8), and even
(R)
C
V
though the whitening matrices may be different, the
results in terms of data decorrelation are statistically
equivalent.
4. The noise enhancement effect
The presence of noise is inherent in every radar sys-
tem; therefore, it is necessary to analyze the perfor-
mance of the whitening transformation, under noisy
conditions. Let V 5 V
S
1 V
N
, where the subscripts S
and N stand for signal and noise components, respec-
tively. When applying the whitening transformation,
both signal and noise are similarly affected:
X 5
W
V 5
W
V 1
W
V 5 X 1 X .
S NSN
(9)
For simplicity, we dropped the subscript ‘‘n’’ that is
used to indicate sample time. From (9), we can see that
the signal is whitened and the noise, which was white
prior to the whitening transformation, becomes colored.
Let us apply the transformation matrix to the data (8)
and compute the range-time correlation for the ran-
(R)
R
X
dom vector X using the expectation operation E[.] as
(R)
TT
T
R
5 E[X*X ] 5
W
*E[V*V ]
W
. (10)
X
The correlation matrix of V (signal plus noise) is given
by S
V
1 N
V
I
, where is given in (4);
I
is the L-
(R)(R)
CC
VV
SS
by-L identity matrix; S
V
is the signal power; and N
V
is
N
OVEMBER
2003 1437TORRES AND ZRNIC
´
the noise power. Then, substituting the expression above
for the correlation matrix of V in (10), using (6), dis-
tributing the matrix product, and using (5),
(R)(R)
T
21 21
R
5 (
H
)*[S
C
1 N
I
](
H
)
XV
VV
S
21
T
5 S
I
1 N (
HH
*) . (11)
VV
It is thus evident that the range-time correlation of X
is the sum of a signal ( ) and a noise ( ) component,
(R)(R)
RR
XX
SN
where 5 S
V
I
and 5 N
V
(
H
T
H
*)
21
. By definition,
(R)(R)
RR
XX
SN
X
S
is white because its correlation is a diagonal matrix.
Additionally, all components have identical power S
V
because is a scalar multiple of the identity matrix.
(R)
R
X
S
On the other hand, the noise becomes colored, and its
mean power after whitening can be computed by av-
eraging the powers of the individual components of X
N
,
which correspond to the diagonal elements of . Then,
(R)
R
X
N
1 N
V
(R)
21
T
N 5 tr[
R
] 5 tr[(
HH
*) ], (12)
X
X
N
LL
where tr(.) is the matrix trace operation. Therefore, the
noise enhancement factor (NEF) defined as NEF 5 N
X
/
N
V
can be obtained from (12) as
(R)
21 21
NEF 5 L tr{[
C
]},
V
S
(13)
where we used the cyclic property of the trace and the
matrix decomposition in (5). For an ideal system (i.e.,
a system with a rectangular transmitted pulse and radar
bandwidth much larger than the reciprocal of pulse
width) (l) 5 1 2 | l |/L for |l | , L, and the trace
(R)
r
V
in (13) can be computed using (A32) to obtain NEF 5
L
2
(L 1 1)
21
. Note that an extra L factor should be added
if comparing with the noise power in the classical pro-
cessing. To effectively oversample by a factor of L, an
L-times larger bandwidth than the reciprocal of the pulse
width is needed; this increases the noise power by the
same factor. Under these considerations, the effective
noise increase over the matched filter case (for an ideal
system) becomes L
3
(L 1 1)
21
.
The trade-off between noise enhancement (radar sen-
sitivity) and variance reduction makes the whitening
transformation useful in cases of relatively large SNR.
For weather surveillance radars, the SNR of signals
from storms is large and the effects of noise when using
the whitening transformation are negligible. For ex-
ample, a 3 mm h
21
rain in the Weather Surveillance
Radar-1988 Doppler (WSR-88D) produces an SNR of
37 dB at 50 km. Clearly, for measuring light rain, the
noise enhancement by whitening would not be a prob-
lem to the full 230 km of required coverage. However,
the threshold for display is set at an SNR of 6 dB mainly
to observe snow, which typically has smaller reflectiv-
ity. This SNR falls in the range where noise enhance-
ment dominates and may preclude the use of whitening.
A solution to the noise enhancement problem is to
relax the whitening requirements and select a transfor-
mation such that the output noise power is also mini-
mized. A transformation that is optimized based on the
minimum mean-square error (MMSE) criterion accom-
plishes the desired goal but requires a priori knowledge
of the SNR at every range location (Ebbini et al. 1993).
Alternatively, we can look at the same problem in terms
of the eigenvalues of the range-time correlation matrix.
The ability to limit the gain of the whitening transfor-
mation to reduce the noise enhancement effect arises
from the relation between the eigenvalues of a corre-
lation matrix and the corresponding power spectral den-
sity. That is, the range spanned by the power spectral
density matches closely the range of eigenvalues (John-
son and Dudgeon 1993). Accordingly, by limiting the
span of eigenvalues, it is possible to place a bound on
the gain of the transformation (Torres 2001). The anal-
ysis of these and other suboptimal techniques is a subject
for further study.
5. Spectral moment estimators
The estimation of spectral moments using a whitening
transformation on oversampled data is performed in
three steps. First, oversampled data in range are whit-
ened as discussed in section 3. Then, the sample-time
autocorrelation at lags zero and one are estimated for
each range location, and these estimates are averaged
(in range) to reduce the standard errors. Finally, these
improved correlation estimates are used to compute
power, Doppler velocity, and Doppler spectrum width
with the usual algorithms (Doviak and Zrnic´ 1993).
Whitening-transformation-based (WTB) estimators for
the spectral moments and the theoretically derived equa-
tions for their variances (appendix A) are presented next.
These theoretical equations are used to compare and
validate simulation results. Further, they clearly show
the interplay of the various variables in the reduction
of variances.
a. Signal power estimator
The WTB power estimator for oversampled signals
in noise is given by
L21 M21
1
2
ˆˆ
S 5 S 5 |X(l, m)| 2 N(NEF), (14)
OO
X
(WTB)
LM
l50 m50
where L is the oversampling factor, M is the number of
pulses, N is the noise power, NEF is the noise enhance-
ment factor in (13), and X(l, n) is the whitened over-
sampled weather signal as in (7). Using the results in
the appendix for an ideal system, the normalized stan-
dard deviation of WTB signal power estimates is ob-
tained from (A34) as
ˆ
SD{S}
1112LN
(WTB)
51
12
[
SLL1 1 S
ÏM 2
s
Ï
p
vn
1/2
2
2
L(3L 1 2L 2 3) N
1 , (15)
2
12
]
2(L 1 1) S