Spectral Signature Calculations and Target ... - Semantic Scholar

1430

IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 55, NO. 4, AUGUST 2006

Spectral Signature Calculations and Target Tracking for Remote Sensing Mark B. Yeary, Senior Member, IEEE, Yan Zhai, Student Member, IEEE, Tian-You Yu, Shamim Nematifar, Student Member, IEEE, and Alan Shapiro

Abstract—Enhanced tornado detection and tracking can prevent loss of life and property damage. The research weather surveillance radar (WSR)-88D locally operated by the National Severe Storms Laboratory (NSSL) in Norman, OK, has the unique capability of collecting massive volumes of time-series data over many hours, which provides a rich environment for evaluating our new postprocessing algorithms. With the advent of more memory and computing power, new state-of-the-art algorithms can be explored. In this paper, an approach of identifying tornado vortices in Doppler spectra is proposed and investigated through the use of neural networks. Once the coordinate of the tornado has been established, the research question becomes the following: Can we apply target tracking algorithms to a volume of radar data to make estimations about where the tornado is going? In recent years, particle filters have attracted great attention in several research communities. These filters are used in problems where time-varying signals must be processed in real time, and the objective is to estimate various unknowns of the signals and to detect events described by the signals. The standard solutions of such problems in many applications are based on the Kalman or extended Kalman filters. In situations when the models that describe the behavior of the system are highly nonlinear and/or the noise that distorts the signals is non-Gaussian, the Kalman-filter-based algorithms provide solutions that may be far from optimal. Here, the path of the tornado follows a path that may be described by a set of nonlinear equations. To estimate the path, the particle filter will provide the better estimates. In addition to the single WSR-88D sensor designs, data fusion and tracing designs are also given for a four-node remote sensor network in central Oklahoma. By incorporating the data from each of the sensors, improvements in tracking are illustrated. The particle-filtering algorithms are especially effective in a networked system of sensors when they are in a data-fusion setting.

Manuscript received June 15, 2005; revised April 20, 2006. This work was supported in part by the NOAA under CSTAR Grant NA17RJ1227, by the National Science Foundation’s Course, Curriculum, and Laboratory Improvement Program under Grant NSF-0410564, and by the Engineering Research Centers Program of the National Science Foundation under NSF Award Number 0313747. This center is comprised of many collaborators and four universities: the University of Massachusetts at Amherst (lead institution); the University of Oklahoma, Norman; Colorado State University, Fort Collins; and the University of Puerto Rico at Mayaguez. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect those of the National Science Foundation. This was accomplished in association with the Adaptive Digital Signal and Array Processing and the Doppler Weather Radar Signal Processing courses. M. B. Yeary, Y. Zhai, T.-Y. Yu, and S. Nematifar are with the DSP and Embedded Systems Laboratory, Department of Electrical and Computer Engineering, University of Oklahoma, Norman, OK 73019 USA (e-mail: yeary@ ieee.org; [email protected]; [email protected]; [email protected]). A. Shapiro is with the School of Meteorology, University of Oklahoma, Norman, OK 73019 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TIM.2006.876574

Index Terms—Digital signal processing, National Weather Radar Testbed (NWRT), particle filters, radar measurements, real-time sensor instrumentation, remote sensing, sensor networks, spectral signature calculations, state estimation, weather surveillance radar (WSR)-88D (KOUN) radar.

I. I NTRODUCTION AND P ROPOSED A PPROACH A. Remote Sensing: Weather Surveillance Radar (WSR) About one third of the nation’s $10 trillion economy is sensitive to climate variability and weather [1]–[3]. Networked remote sensor systems, such as the ubiquitous WSR-88D (WSR—1988 Doppler) in Fig. 1, comprise the approximately 150 weather radars in the continental U.S. that are in a national network to provide the bulk of the nation’s weather information [4], [5]. It should also be mentioned that there are 47 more terminal Doppler weather radars (TDWRs) that are similar to the WSR-88D at airports throughout the nation. A prominent feature of the TDWR is that it can probe inside storms and measure dangerous wind shifts, such as those linked to windshear and downburst, which pose a threat to the aircraft during takeoff and landing [6]. Long-term warnings have improved greatly over the last five years and are now being used for critical decision making [7]. Further improvements are being aimed at providing longer warning lead times before severe weather events, better quantification of forecast uncertainties in hurricanes and floods, and tools for integrating probabilistic forecasts with other data sets [7]. WSRs, particularly the WSR-88D, have shown to be an important tool to observe severe and hazardous weather remotely, and to provide operational forecasters prompt information of rapidly evolving phenomena [8]. One of the most extreme weather phenomena is the tornado that can produce destructive wind speeds as high as 140 m/s [9]. Early and accurate detection of tornadic vortices can increase the lead time for tornado warning to prevent loss of life and property damage. A key component of the current tornado-detection algorithms on WSR-88D radars is a search for the presence of strong localized azimuthal shear of the radial winds. However, because the radar sample volume increases with distance from the radar, the shear signature deteriorates as the range of the tornado increases. In this paper, an independent approach of identifying the tornado vortices in Doppler spectra is proposed and investigated. The Doppler spectrum reveals the weighted velocity distribution within the radar volume. The Level II data (reflectivity, mean Doppler velocity, and spectrum width) are obtained by the first

0018-9456/$20.00 © 2006 IEEE

YEARY et al.: SPECTRAL SIGNATURE CALCULATIONS AND TARGET TRACKING FOR REMOTE SENSING

1431

Fig. 1. Hardware architecture of the operational system. The research WSR-88D (KOUN) operated by the NSSL in Norman, OK, has the unique capability of collecting massive volumes of Level I time-series data over many hours. A unique feature of the KOUN is that it has a dual polarization capability and hence, the two receivers. It can change the modes from dual to single polarization (picture courtesy Alan Zahari, NSSL, NOAA).

three moments of the distribution [12]. Therefore, a characteristic spectrum for tornadoes may still facilitate their identification when the shear signature becomes difficult to identify. Doppler spectra of tornadoes have a distinct character that set them apart from other spectra. Flat and bimodal tornado spectral signatures (TSS) have been shown using simulated data [10], [11]. Although the history of tornado measurements is long, there have been only a few successes in obtaining spectra data. This is largely because neither the technology to process spectra, nor the technology to record voluminous amounts of time-series data were available. The present-day radar technology and computer resources are now advanced enough to study TSS in a systematic manner. As a pulsed Doppler radar intermittently radiates electromagnetic energy, it also receives energy that is reflected back to it from various targets. These targets may be rain drops, hail stones, birds, aircraft, and other meteorological and nonmeteorological reflections. When the energy is received at the radar, these analog signals are digitized, and the result is known as Level 1 data. The research WSR-88D (KOUN) operated by the National Severe Storms Laboratory (NSSL) is a pulsed Doppler radar that has the unique capability of collecting massive volumes of Level I time-series data over many hours. Doppler spectra can be obtained by processing these Level I data after the event. Doppler spectra from a tornado outbreak in central Oklahoma on May 10, 2003 are presented, as illustrated later in this paper.

B. Doppler Spectra and Spectral Analysis The mean Doppler velocity v¯(R0 ) is defined by the following equation: ∞ v¯(R0 ) =

vSn (R0 , v)dv

(1)

−∞

where Sn (R0 , v) is the Doppler spectrum normalized by the total power, where R0 is the center of the radar volume, and where v is the radial velocity component. Therefore, the mean Doppler velocity is the first moment of the spectrum, which is estimated from one realization (or probe volume). The current tornado-detection algorithms rely on the difference of the mean Doppler velocity between adjacent radar volumes in the azimuth. Tornado spectra observed by the NSSL research WSR-88D (KOUN) located at Norman, OK, on May 10, 2003, are analyzed. An F2–F3 tornado was reported approximately during the interval of time 9:30–10:00 PM (central time) starting south of Edmond, OK, (more details of the tornado can be found from NOAA National Climate Data Center (NCDC), http://www.ncdc.noaa.gov/oa/ncdc.html). Level I time-series data were collected during the entire period of tornado. As mentioned in this section, the Level I data are produced by the received reflected energy. That is, the energy is reflected back to the radar when it bounces back from an object.

1432


Fig. 2. PPI of the WSR-88D reflectivity at 9:43 PM, May 10, 2003. The small box outlines the location of the region of interest. The data for the spectral calculations in Fig. 7 originate from this box. The small white dots indicate the tornado-damage path from the ground survey.

These objects or targets are typically referred to as scatterers, since the energy is scattered back to the radar. Here, a target is considered to be a water droplet within a tornado. The Level I data are represented in a format called in-phase (I) and quadrature (Q). A one-dimensional (1-D) Fourier transform can be used to transform these time-domain I and Q samples into the frequency domain or Doppler-velocity spectrum. Three spectral moments of weather phenomena from the I and Q data are 1) weather signal power that is an indication of reflectivity, 2) mean Doppler velocity, and 3) spectrum width. Hence, reflectivities and mean Doppler velocities were obtained every half degree in azimuth using the autocovariance method [12]. This method yields results that are similar to the spectral analysis, but it is more computationally efficient since only the first three moments are being computed. As the WSR-88D rotates to capture an image of the atmosphere, several distinct elevation angles are used to do this. The lowest two or three elevations are typically used to search for tornados, as we have done here. To display this information, a plan position indicator (PPI), also known as a polar plot indicator, is used. In other words, it is a display of consecutive radials of data on a cone at constant elevation. The PPI of the reflectivity from the lowest elevation angle (0.439◦ ) at 9:43 PM, May 10, 2003, is given in Fig. 2. A well-defined hook signature and strong azimuthal shear were observed at 6 km west and 39 km north of the radar, indicating the existence of a tornado. As noted by Brown et al. [13], a hook-shaped reflectivity feature that is typically on the right rear flank of a severe storm indicates the presence of a mesocyclone that has the potential of having a tornado form within it. As discussed in [11] and [14], from simulation and real data, it is clear that a tornado can produce broad and flat spectra that are similar to the white-noise spectrum. However, significant signal power still remains in the tornado spectrum. Fig. 3 depicts the relative location of the radar and the path of the mesocyclone (prelude to a tornado).

Fig. 3. WSR-88D is located at the origin of the coordinate system for this single-sensor design. The path of the tornado begins at approximately 40 km north of the radar.

C. Tracking-Algorithm Development Once the location of the tornado has been established, tracking algorithms can be developed to estimate the movement of the vortex and to provide feedback into the tornadolocation-decision process to improve accuracy. These tracking algorithms are based on the emerging field of Monte Carlo sequential importance sampling (SIS), which is also nicknamed particle filtering (PF), as described later in Section III. The PF algorithms are especially effective in a networked system of sensors when they are in a data-fusion setting. By definition, the basic idea of a particle filter is to approximate the probability density function (pdf) of the estimated states by a large set of randomly generated samples (particles), which will ultimately render better estimates of the desired states, as depicted in Fig. 4. These filters are used in problems where time-varying signals must be processed in real time, and the objective is to estimate various unknowns of the signals and to detect events described by the signals. The standard solutions of such problems in many applications are based on the Kalman or extended Kalman filters (EKFs). In situations when the problems are highly nonlinear or the noise that distorts the signals is nonGaussian, the Kalman filters provide solutions that may be far from optimal, since the Kalman filter is designed to update the posterior distribution’s mean and variance. However, this distribution is estimated by a set of many discrete samples with the particle filter. Finally, the ultimate duty of the filtering process, be it the Kalman filter or the particle filter, is to estimate the state of a variable. Several comparisons are made at this time, based on our studies. In our nonlinear scalar example given in [15], 1000 Monte Carlo runs were implemented to produce the statistical performance, i.e., the root mean-square error (rmse), of a system designed to estimate the state of a nonlinear state variable. Fig. 5 depicts three curves. These are the EKF, a standard particle filter known as the bootstrap filter (bootstrap), and a particle filter based on state partitioning and parallel EKFs (PF-SP-PEKFs). The improvement of PF-SP-PEKF over EKF in terms of the mean rmse is approximately 78%, and more details are available in [15].


Fig. 4.

1433

Evolution of the moving target’s pdf and its adaptive particle approximation.

(an input layer, a hidden layer consisting of 40 neurons, and an output layer) provided us with the most efficient result in terms of training time and detection. A. Adaptive Algorithms for Feedforward–Backpropagation (FFB) Neural Networks

Fig. 5.

Kalman and particle-filter performance comparisons.

D. Project Genesis To quote a recent report by the National Science Board titled Science and Engineering Infrastructure for the 21st Century [16] “fueled by exponential growth in computer power, communication bandwidth and data storage, the Nation’s research infrastructure is increasingly characterized by interconnected, distributed systems of hardware, solfware, information bases, and expert systems.” The authors of this paper believe in the same philosophy, which serves as the unifying theme of Sections I-A–C. It is noted that this paper builds upon and improves upon our IEEE-International Multimedia Teleconferencing Consortium (IMTC) 2005 conference paper [17]. In brief, this new paper includes a more detailed description of the tracking-algorithm development, more laboratory results, expanded theoretical details, and comprehensive models. II. A PPLYING N EURAL N ETWORKS TO THE T ORNADO -D ETECTION P ROBLEM The task at hand was to train a neural network to recognize the broad and flat Doppler spectra associated with tornados. Our team examined a diverse class of artificial neural networks for the most efficient candidate. In particular, we looked at their architecture, training time, and the generalization ability of the network (how well a given network responds to new inputs). After exploring a variety of network architectures and training algorithms, our team decided that a feedforward–backpropagation neural network with three layers

An FFB network is a generalization of the network architecture of Fig. 6. Hornik et al. [18] have shown that a threelayer FFB neural network with sigmoid transfer functions in the hidden layer and linear transfer functions in the output layer can approximate any continuous function to an arbitrary precision. The scaled-conjugate-gradient (SCG) training algorithm was employed to adjust the weights in the network. SCG is a powerful extension of a class of optimization techniques known as the CG methods. The idea behind the training is to adjust the weights of a neural network such that a desired mapping from the input to the output is achieved. The squared mapping error is defined as F (Wk ) = e2 (k) = [t(k) − a(k)]2 , where Wk is the vector of network weights and biases, e(k) is the error signal, t(k) is a vector of desired target outputs, and a(k) is a vector current outputs of the network at the kth iteration. In order to minimize the squared mapping error, CG algorithms traverse the error surface along a set of conjugate directions. A set of vectors p1 , p2 , . . . , pn constructs a conjugate system if there is a nonsingular symmetric N × N matrix A, such that pT i Apj = 0, for i = j, i = 1, 2, . . . , n. Scales [19] and Gill et al. [20] have shown that the exact minimum of any quadratic function can be found by moving along a set of conjugate directions p1 , p2 , . . . , pn . Such a conjugate system can be constructed by first defining p1 = −g1 , then, subsequently, pk = −gk + βk pk−1 . The variable gk = ∇F (W)|W=Wk is the gradient evaluated at the old weights, i.e., Wk , k = 2, . . . , n, and βk is given by βk =

gkT gk . T g gk−1 k−1

(2)

CG methods adaptively update the weights in a sequence of steps by moving along the conjugate directions Wk+1 = Wk + αk pk . The variable αk is called the learning rate parameter and determines how far the movements occur along each search direction. For a quadratic function, αk can be analytically calculated to minimize the performance function along the search path as defined by the conjugate directions [21]. Thus αk = −

gkT gk . T gk−1 Ak gk−1

(3)

1434


Fig. 6. Feedforward/backpropagation network that was employed. Thirty-two input nodes were followed by 40 hidden neurons. The CG algorithm was employed to train the 1280 weights of the first hidden layer, which are denoted by W1 . Similarly, the CG algorithm was employed to train the 40 second-layer weights, which are denoted by W2 .

The variable Ak = ∇2 F (W) is the Hessian matrix for a quadratic function. It is noticed that, in many practical applications, the performance function is a higher order function of hundreds of weight parameters. However, a local quadratic approximation can be achieved by applying the Taylor expansion formula [21]. Finally, the SCG algorithm employs an approximation of the Hessian matrix to reduce the computational time and complexity (see [22] for more details). B. Neural Network Training and Results After atmospheric testing, the WSR-88D Level I time-series data consisted of thousands of range gates of data, each gate providing 32 complex I and Q samples taken in a 1◦ beam width (overlapping every 1/2◦ ). The Doppler spectrum of each vector of 32 samples is calculated by taking the Fourier transform of the I and Q time-series samples, taking its magnitude squared, and expressing it in units of decibels. Further preprocessing of the data was performed by passing each Doppler spectrum through a four-point lowpass filter. The smoothed Doppler spectra form the input feature vectors that are fed into our network. Training consisted of the spectra of 60 range gates. A “range gate” denotes a specific atmospheric volume that was sampled at a particular range. Eight of the range gates were defined as tornadic, while the others were not. After much experimentation in the training process, this preprocessing procedure enhanced the network’s ability to discriminate between targets and nontargets (tornadic regions or not). A subset of the resulting data is shown in Fig. 7. Each spectrum has 32 points, and the radar’s unambiguous velocity is 32 m/s. The network was able to adapt itself to the training dataset in less than 1 min. Subsequent presentation of the input datasets to the trained network resulted in eight out of eight recognition of the tornadic patterns. The system was designed to have a low false-alarm rate during the training phase of the neural network. Following the training, during the pattern-recognition phase of the network, the network achieved a more reasonable falsealarm rate while providing a high accuracy on the targets of interest. It is noted that determining the probability of detection

is not straightforward since scientists typically do not have access to the damage path on the ground. For more details about the probability of detection, see [23] and [24]. After training the neural network to positively detect tornadic patterns on the training data in Fig. 7, the next step is to test the ability of our trained neural network to generalize and detect tornadolike patterns on new datasets consisting of spectral data from hundreds of range gates. Our team acquired a set of four data files, containing I and Q samples, recorded within a halfhour time span during the May 10, 1999 tornado outbreak in central Oklahoma. Each file contained 32 complex I and Q samples taken in a 1◦ beam width (overlapping every 1/2◦ ) over an azimuthal angle of 20◦ and a 10-km range, resulting in 1600 vectors of length 32 each. In each case, our trained neural network was able to successfully scan through these massive volumes of data and pinpoint a very good estimate of the location of each tornado in less than a second. Such expedient performance is due to the fact that the data flow through the network in a feedforward scheme, and all the weight coefficients are fixed (no feedback or parameter adjustment is required). Finally, as a practical consideration, we note that a hardware implementation of this network, with computations performed within neurons operating in parallel, can provide very high realtime performance [25].

III. T ARGET T RACKING U SING P ARTICLE F ILTERS IN N ETWORKED S ENSOR S YSTEMS Once the coordinates of the tornado have been established via our neural network algorithms, tracking algorithms based on particle filters can be developed to estimate the movement of the vortex. A. Introduction to the Particle Filter The PF technique has drawn an increasing attention during the recent years because of its superior ability to deal with processes that have nonlinear models and/or contain nonGaussian noise sources [15], [26]–[30]. PF is a sequential Monte Carlo method based on the concept of importance sampling and Bayesian theory. For the sake of completeness, in this section, we briefly recall the PF technique. Consider the following nonlinear system: xt = f (xt−1 ) + vt−1 yt = h(xt ) + wt

(4)

where xt and yt denote the states and the observations of the system at time t. The process noise and the observation noise are represented by vt and wt . Both f (·) and h(·) are nonlinear functions. The filtering objective is to simulate the posterior distribution p(x0:t |y1:t ) via a particle filter. According to Bayes’ theorem, this distribution is given below: p(x0:t |y1:t ) =

p(y1:t |x0:t )p(x0:t ) . p(y1:t |x0:t )p(x0:t )dx0:t

(5)


1435

Fig. 7. Laboratory data: Samples of the Doppler spectra from the box illustrated in Fig. 2. Every column has seven range gates corresponding to 1/2◦ angular sampling. Here, a search area is defined by a distance of 39.125 to 40.625 km and spanning 5.10◦ to 12.02◦ in angle. The independent axis of each graph above spans from −32 to +32 m/s. The dependent variable 20 log10 |V | is plotted on the vertical axis. This axis has a span of −70 to 10 dB. Thirty-two digital samples were employed to calculate each of the spectral plots above. As described in the paper, tornadic activity is denoted by spectra that are relatively flat and have a high signal-to-noise ratio (SNR). This is exemplified by the regions defined by 8.5◦ and 9.0◦ , spanning from 39.375 to 40.125 km.

By applying the SIS, the posterior distribution can be approximated by a set of weighted samples shown previously in Fig. 4 as follows: p(x0:t |y1:t ) ≈

Ns

(i) ω ˜ (i) (t)δ x0:t − x0:t

(6)

i=1

where Ns denotes the number of particles (samples). The (i) variable ω ˜ t is called the normalized importance weight and is given by (i)

ω (i) ω ˜ t = N t s

j=1

where (i) ωt

(j)

(7)

ωt

(i) (i) (i) p yt |xt p xt |xt−1 (i) = ωt−1 (i) (i) q xt |x0:t , y1:t

(8)

(i)

and ωt denotes the importance weight associated with the ith particle at time t. The particles (samples), which are denoted by (i) xt , are drawn from the proposal distribution (i) (i) (9) xt ∼ q xt |x0:t−1 , y0:t . A good proposal distribution should provide an enough support from the true posterior distribution of system states, and it should be easy to draw samples from this distribution. As our previous results have indicated [15], and as depicted in

Fig. 8. As depicted here, the pdf of a state variable evolves over time. Carefully modeling this pdf is essential to state estimation.

Fig. 8, the pdf of a state variable as it evolves, especially when the underlying process is nonlinear, deviates from the Gaussian assumptions of the Kalman filter. The problem of tornado tracking is extremely difficult to solve because of the nonlinearities and the chaotic behavior of these phenomena. Thus, the particle filter is an excellent tool. B. Target Tracking With the Particle Filter: Nonmaneuvering Case Target tracking is a classic state-estimation problem, and is especially challenging due to the nonlinear observation model.

1436


In this problem, the target is usually represented by a state vector defined as follows:

where (xt , yt ) is the target location in a Cartesian coordinates at time t. The variables vxt and vyt denote the target velocities along the x-axis and y-axis, respectively. The discrete time model for the kinematics of a nonmaneuvering target is given as [31]

presents additional tracking challenges, as compared to nonmaneuvering approaches. When tracking a target that can perform maneuvers, the estimations generated by a bootstrap filter based on the aforementioned model [(11) and (14)] are not always accurate enough [31], [33]. The problem of maneuvering target tracking is often referred to as a jump Markov process [31], [34], in which the system is assumed to operate according to one model from a finite set of hypothetical models, known as regimes or modes. The state process model of this tracking problem is defined as

xt = F · xt−1 + Γ · vt−1

(11)

xt = Frt (xt−1 ) + Γ · vt−1

(12)

where rt is the regime variable, which is governed by a regime transition matrix defined as   π11 π12 · · · πNr 1  π21 π22 · · · ·    (16) Π= . . ..  . .. ..  .. .  π1Nr · · · · · · πNr Nr

xt = [xt

yt

vxt

vyt ]T

(10)

where 

1 0 F= 0 0

0 T 1 0 0 1 0 0



0 T  0 1



2

T /2  0 Γ=  T 0



0 T /2   0  T 2

and where T is the sampling rate. The process noise vt is a 2 × 1 noise vector with a distribution vt ∼ N (0, Qv ), where

2 σv 0 Qv = . (13) 0 σv2 Here, the measurement of the tracking system is the angle between the observation platform (the coordinate origin) and the location of the target. It is also known as the direction of arrival (DOA), which is given by xt (14) zt = arctan + wt . yt It is obvious that the observation model, given in the above equation, of the tracking problem involves strong nonlinearities. In addition, there is no direct measurement of the current velocities. Historically, the EKF has been used to solve the tracking problem. However, it was found in the current literature that the estimates from EKF are not only easily to diverge from the true target trajectory but also always bare large estimation errors [31], [32]. Recently, PF has emerged as a promising target tracking algorithm due to its accuracy and robustness when dealing with nonlinear systems and non-Gaussian noises. The bootstrap filter is a special case of particle filter, in which the filter takes the (i) (i) state transition prior p(xt |xt−1 ) as the proposal distribution (i) (i) (i) (i) p(xt |x0:t , y1:t ) so that (8) is reduced to ωt ∝ p(yt |xt ). In other words, the importance weights are proportional to the system likelihood. The bootstrap filter is easy to implement and computationally efficient. As indicated in [31] and [32], the bootstrap filter gives better estimates than EKF in the problem of tracking nonmaneuvering targets. C. Target Tracking With the Multiple-Model Particle Filter: Maneuvering Case The target to be tracked in this example is a maneuvering target. By definition, this implies that the target’s x and y velocities are nonconstant as it makes its path. This aspect

(15)

where πij = Pr(rt = j|rt−1 = i), (i, j ∈ 1, . . . , Nr ), and Nr j=1 πij = 1. The variable Nr is the total number of possible regimes. In addition, the observation model is the same as (14). Usually, dynamic models for a maneuvering target can be divided into two categories: 1) two-dimensional (2-D)/threedimensional (3-D) motion models and 2) coordinate-uncoupled maneuver models [34]. In this paper, we only consider the 2-D horizontal motion model [31] that is composed of three different regimes defined as 1) constant-velocity (CV) model, in which the target has a rectilinear motion. This model is the same as the model discussed in the nonmaneuvering case, i.e., F1 (xt ) = F, where F is given in (11); 2) clockwise coordinated turn (CT) model; and 3) counterclockwise CT model. For r = 2 and 3, these two CT models are given as follows:   (r) (r) 1−cos Ωt T sin Ωt T 1 0 − (r) (r)   Ωt Ωt     (r) (r) 1−cos Ωt T sin Ωt T   r 0 1   (17) (r) (r) F (xt ) =  Ωt Ωt    (r) (r) 0 0 cos Ωt T − sin Ωt T    (r)

0 0

sin Ωt T

(r)

cos Ωt T

(r)

where Ωt denotes the model-conditioned turning rate given as (2)

Ωt

(3)

Ωt

=

=

am

(18)

vx2t + vy2t −am vx2t + vy2t

.

(19)

The variable am in the above equation represents the maneuver acceleration. Fig. 9 illustrates the possible turn result from various maneuver accelerations with the target initial condition given as: xt = [0 0 1 0]T . As indicated by (17)–(19), the two CT models have strong nonlinearities. The interacting multiple model-EKF (IMM-EKF) is a standard way to track


1437

(i)

where xt and zt are defined in (10) and (14). This equation is similar to a single-model PF given in (8), except that the regime (i) variable rt is incorporated. After this, the particle weights are normalized followed by a resampling step. For the sake of completeness, the MMPF algorithm is summarized in Table I. The multimodel bootstrap filter is a special case of MMPF, in (i) (i) (i) which the proposal distribution of the PF, q(xt |xt−1 , rt , zt ), (i) (i) (i) is taken as the state transition prior p(xt |xt−1 , rt ), i.e, in this case, the particle weights only depend on the system’s (i) (i) likelihood p(zt |xt , rt ), as indicated in the previous section. D. Remote-Sensor Geographical Layout

Fig. 9. Possible turns with various maneuver accelerations. Movements to the left are represented by positive turning rates, while movements to the right are represented by negative values. The values above correspond to the variable am in (18) and (19).

maneuvering target. The basic idea of IMM-EKF is to use a separate EKF for each hypothetical model; then, the final estimate is the weighted sum of all the components. The filter weight is a function of the regime transition probabilities and the filter weight at previous time index. The approximation error of IMM-EKF comes from two sources. First of all, each EKF approximates a nonlinear system by linearizing it about a certain operating point. This procedure introduces the linearization error. Second, the IMM approximates the exponentially growing Gaussian mixture with a finite Gaussian mixture that result in an approximation error. These two kinds of error may cause the filter to diverge from the true target trajectory in some scenarios. More recently, various multiple model PF (MMPF) algorithms have been proposed as superior alternatives to IMMEKF to perform the nonlinear filtering with switching dynamic models [33], [35]–[37]. As a novel extension of a single-model particle filter, the MMPF generates state estimates based on all regimes. More specifically, at each iteration, a set regime (i) s samples {rt }N i=1 are generated based on the regime transition matrix Π defined in (16) and the regime samples at previous (i) s time {rt−1 }N i=1 such that (i) (i) Pr rt = j|rt−1 = i = πij

(20)

With the advances in networking technology, a new paradigm of networked radars is on the verge of becoming a reality. The many advantages offered by these networks of radars are manifested by their ability to offer the opportunity to potentially revolutionize the science of radar meteorology. The NSF sponsored center for a collaborative adaptive sensing of the atmosphere (CASA) is a prime example of emerging distributed collaborative and adaptive systems [38], [39]. The University of Oklahoma is a partner of CASA initiative, with three other universities. The vision of the CASA is to revolutionize our ability to observe the lower troposphere through distributed collaborative adaptive sensing (DCAS), improving the ability to detect, understand, and predict severe storms, floods, and other atmospheric and airborne hazards. As depicted in Fig. 10, the system will cover a 7000-km2 region of Southwestern Oklahoma, as described in [40]. The CASA team has engineered the main components of the low-power radar sensing node, including the transmitter, digital receiver, processor, and pedestal. As of the time of this writing, the four-node system has been installed and is being calibrated. The Oklahoma statewide Internet Service Provider (ISP), known as OneNet, provides the high-speed communication infrastructure (towers and associated ISP) that supports the new test-bed project. Since the four-node network is in the final construction phase, it is noted that the results below are based on the simulated data. In this paper, we adopt the multimodel bootstrap filter using multiple sensors to track a tornado as it moves across central Oklahoma. In this tracking scenario, we applied separate multimodel bootstrap filters using DOA measurements from each sensor, then, the overall estimation is given by a weighted sum of these individual estimates. The sensor weight is chosen as the normalized reflectivity received by each radar, which is proportional to the inverse of the distance squared. The algorithm of MMPF using multiple sensors is illustrated in Fig. 11. IV. E XPERIMENTAL R ESULTS

where Ns is the sample size. Then, the state particles are drawn (i) (i) (i) from the proposal q(xt |xt−1 , rt , zt ), and the particle weights are calculated as follows: (i)

ωt

(i) (i) (i) (i) (i) p zt |xt , rt p xt |xt−1 , rt (i) = ωt−1 (i) (i) (i) q xt |xt−1 , rt , zt

(21)

In this section, we apply the multiple-sensor MMPF to track a tornado trajectory using four sensors in a geometrical arrangement similar to the positions of the CASA radars in Fig. 10. In addition, a single-sensor single-model bootstrap filter and a single-sensor MMPF are also implemented for a purpose of comparison. The estimation results of these three tracking algorithms are shown in Fig. 12. In this figure, the true

1438


TABLE I MMPF ALGORITHM

trajectory is denoted by the plus signs (+). The estimates from the single-sensor single-model bootstrap filter are represented by the line with a cross (− × −), while estimates from MMPF with single sensor and with multiple sensors are depicted by the upward triangle line (−−) and the line with dots (−•−), respectively. In addition, each PF algorithm uses 5000 particles. As indicated in the figure, although the single-model singlesensor bootstrap filter can follow the general direction of the target, it failed to detect the small turns made by the target. On the other hand, the multiple-sensor MMPF algorithm is able to detect the maneuvers made by the target and outperforms the single-sensor MMPF. V. F UTURE V ISIONS

Fig. 10. Four-node remote sensor network in central Oklahoma. By incorporating the data from each of the sensors, dramatic improvements in the tracking algorithm can be made.

Enhanced tornado tracking can prevent loss of life and property damage. Moreover, tracking other targets of national importance will continue to evolve and gain attention. Therefore, multifunction radar capabilities are needed, that is, new radars will also need the provisions to track noncooperating airplanes, toxic plumes from industrial accidents, and the like. One example of homeland interest is as follows. In the United States, airport capacity has increased by only 1% in the past ten years, while air traffic increased 37% during that time,


1439

Fig. 11. MMPF tracking algorithm using the DOA measurements from multiple sensors.

Fig. 12. Zoom-in view of Fig. 3 illustrating the tracking results of three algorithms: Single-sensor single-model bootstrap filter (− × −), single-sensor MMPF (−−), and multisensor MMPF (−•−). In addition, the true target location is represented by plus signs (+). It is clearly demonstrated from this figure that the MMPF algorithm with multiple sensors yields very accurate estimates.

1440


been detected, its location is stored, and these coordinates are employed in a stochastic state-estimation algorithm to predict the future movement of the tornado. Rather than relying on the traditional Kalman filter to provide tracking estimates, a Monte Carlo sequential state-estimation technique known as PF was employed. Based on ever expanding computing power, the particle filter allows the posterior probability distribution of the state that defines the tornado’s location to be approximated by a set of many samples, rather than approximating its mean and variance as the Kalman filter does. This provides more reliable estimates of the future path of the vortex. Particle filters work well not only in the single-sensor environment but in a multisensor setting as well. As we have illustrated in this paper, by employing a carefully designed data-fusion center that relies on weighted inputs from several sensors or radars, a comprehensive PF algorithm may be developed to estimate the position of a moving target, such as a tornado, with a high degree of fidelity. Fig. 13. Plus marks indicate the true locations or measurements. The line with the stars indicates the estimates with the extended Kalman filter. The line with the triangles denotes the output of a particle filter. Finally, the line with the circles indicates our estimates with a new method, based on PF, state partitioning, and a bank of extended Kalman filters (known as PF-SP-PEKF, see Section I-C and [42] for more details).

as reported by the American Society of Civil Engineers [41]. It is clear that America’s infrastructure is aging. Providing discipline specific solutions will be extremely costly. However, by working together, a diverse group of scientists and engineers can develop a system of networked radars that can be used in a multifunction capacity to provide both weather and target-tracking data. This strategic alliance greatly reduces costs, while providing enormous benefits to the public. The electronically steerable beams of the new multifunction phased array radar at the National Weather Radar Testbed (NWRT) in Norman, OK, have the ability to detect a variety of airborne targets. Fig. 13 depicts aircraft detections and tracking results at the Oklahoma City International Airport, based on the data taken on September 15, 2005 [42]. VI. C ONCLUDING R EMARKS The main focus of our research has been to develop an optimum method of identifying tornadic patterns in the Doppler spectrum of a large ensemble of weather data and subsequently employ this identification data into a tracking algorithm. To this end, we examined a novel tornado-prediction technique based on artificial neural networks (as pattern recognizers). Traditional tornado-detection techniques rely on the difference of mean Doppler velocities between adjacent radar volumes. According to these methods, strong localized gate-togate shears (or azimuthal shears) are an indication of a tornado. However, this approach may break down when the tornado is far from the radar since the azimuthal resolution degrades with increasing distance. Tornadic vortices demonstrate a distinct pattern in the Doppler spectrum that differs from ordinary weather signatures. An FFB artificial neural network as a pattern recognizer has the potential adaptability required to detect a wide range of tornadic patterns. After the tornado has

ACKNOWLEDGMENT The authors would like to thank the scientists who maintain, operate, and collect data from the KOUN radar at the National Severe Storms Laboratory, Norman, OK, for their dedicated work. R EFERENCES [1] National Research Council, The Atmospheric Sciences Entering the Twenty First Century, 1998, Washington, DC: National Academy. [2] Bureau of Economic Analysis. (2002). GDP by Industry, [Online]. Available: www.bea.gov/bea/dn2/gpo.htm [3] J. Dutton, “Opportunities and priorities in a new era for weather and climate services,” Bull. Amer. Meteorol. Soc., vol. 83, no. 9, pp. 1303–1311, Sep. 2002. [4] G. E. Klazura and D. A. Imy, “A description of the initial set of analysis products available from the NEXRAD WSR-88D system,” Bull. Amer. Meteorol. Soc., vol. 74, no. 7, pp. 1293–1311, 1993. [5] T. Crum and R. Alberty, “The WSR-88D and the WSR-88D operational support facility,” Bull. Amer. Meteorol. Soc., vol. 74, no. 9, pp. 1669–1687, 1993. [6] Terminal Doppler Weather Radar. [Online]. Available: www.raytheon. com/newsroom/photogal/tdwr_h.htm [7] National Research Council, Making Climate Forecasts Matter, 1999, Washington, DC: National Academy. [8] R. Serafin and J. Wilson, “Operational weather radar in the United States: Progress and opportunity,” Bull. Amer. Meteorol. Soc., vol. 81, no. 3, pp. 501–518, 2000. [9] R. Davies-Jones, R. Trapp, and H. Bluestein, “Tornadoes and tornadic storms,” in Severe Convective Storms. Boston, MA: Amer. Meteorological Soc., 2001, ch. 5, pp. 167–221. [10] D. S. Zrnić and R. J. Doviak, “Velocity spectra of vortices scanned with a pulsed-doppler radar,” J. Appl. Meteorol., vol. 14, no. 8, pp. 1531–1539, 1975. [11] T.-Y. Yu, D. S. Zrnić, A. Shapiro, and M. B. Yeary, “Feasibility of earlier tornado detection using doppler spectra,” in Proc. 31st Conf. Radar Meteorol., 2003, pp. 333–336. [12] R. J. Doviak and D. S. Zrnić, Doppler Radar and Weather Observations. San Diego, CA: Academic, 1993. [13] R. Brown et al., “Improved detection of severe storms using experimental fine-resolution WSR-88D measurements,” Weather Forecast., vol. 20, no. 1, pp. 3–14, Feb. 2005. [14] T.-Y. Yu, A. Shapiro, D. S. Zrnić, M. Foster, D. Andra, R. Doviak, and M. Yeary, “Tornado spectral signature observed by WSR-88D,” in Proc. 22nd Severe Local Storms Conf., Hyannis, MA, Oct. 4–8, 2004, pp. 1–5. Session 8b.1. [15] Y. Zhai and M. Yeary, “A nonlinear state estimation technique based on sequential importance sampling and parallel filter banks,” in Proc. IEEE Conf. Control Appl., Toronto, ON, Canada, 2005, pp. 2071–2075.


[16] National Science Board, “Science and engineering infrastructure for the 21st century,” National Science Foundation, Arlington, VA, Rep. NSB 02-190, 2003. [17] M. Yeary, T.-Y. Yu, S. Nematifar, and A. Shapiro, “Spectral signature calculations for remote sensing,” in Proc. IEEE-Instrum. Meas. Technol. Conf., May 2005, pp. 2071–2075. [18] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural Netw., vol. 2, no. 5, pp. 359–366, 1989. [19] L. Scales, Introduction to Non-Linear Optimization. New York: Springer-Verlag, 1985. [20] P. Gill, W. Murray, and M. Wright, Practical Optimization. New York: Academic, 1981. [21] M. Hagan, H. Demuth, and M. Beale, Neural Network Design. Boulder, CO: University Colorado Press, 1996. [22] M. Møller, “A scaled conjugate gradient algorithm for fast supervised learning,” Comput. Sci. Dept., Univ. Aarhus, Aarhus, Denmark, Tech. Rep. PB-339, Nov. 1990. [23] A. Witt, M. Eilts, G. Stumpf, E. Mitchell, J. Johnson, and K. Thomas, “Evaluating the performance of WSR-88D severe storm detection algorithms,” Weather Forecast., vol. 13, no. 2, pp. 513–518, Jun. 1998. [24] E. Mitchell, S. Vasiloff, G. Stumpf, A. Witt, M. Eilts, J. Johnson, and K. Thomas, “The national severe storms laboratory tornado detection algorithm,” Weather Forecast., vol. 13, no. 2, pp. 352–366, 1998. [25] B. E. Boser, E. Säckinger, J. Bromley, Y. LeCun, and L. Jackel, “Hardware requirements for neural network pattern classifiers: A case study and implementation,” IEEE Micro, vol. 12, no. 1, pp. 32–40, Feb. 1992. [26] A. Doucet, N. de Freitas, and N. Gordon, Sequential Monte Carlo Methods in Practice. New York: Springer-Verlag, 2001. [27] P. Djurić et al., “Particle filtering,” IEEE Signal Process. Mag., vol. 20, no. 5, pp. 19–38, Sep. 2003. [28] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” IEEE Trans. Signal Process., vol. 50, no. 2, pp. 174–188, Feb. 2002. [29] Y. Zhai and M. Yeary, “Implementing particle filters with MetropolisHastings algorithms,” in Proc. IEEE-Region 5 Conf., Apr. 2004, pp. 149–152. [30] A. Doucet, S. Godsill, and C. Andrieu, “On sequential Monte Carlo sampling methods for Bayesian filtering,” Stat. Comput., vol. 10, no. 3, pp. 197–208, 2000. [31] B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications. Boston, MA: Artech House, 2004. [32] N. Gordon and D. Salmond, “Novel approach to nonlinear/non-Gaussian Bayesian state estimation,” Proc. Inst. Elect. Eng.—Radar Signal Process., vol. 140, no. 2, pp. 107–113, Apr. 1993. [33] C. Musso, N. Oudjane, and F. LeGland, “Improving regularized particle filter,” in Sequential Monte Carlo Methods in Practice, A. Doucet, N. de Freitas, and N. Gordon, Eds. New York: Springer-Verlag, 2001. [34] X. Li and V. Jilkov, “Survey of maneuvering target tracking. Part 1: Dynamic models,” IEEE Trans. Aerosp. Electron. Syst., vol. 39, no. 4, pp. 1333–1364, Oct. 2003. [35] S. McGinnity and G. Irwin, “Multiple model bootstrap filter for maneuvering target tracking,” IEEE Trans. Aerosp. Electron. Syst., vol. 36, no. 3, pp. 1006–1012, Jul. 2000. [36] D. Angelova, T. Semerdjiev, V. Jilkov, and E. Semerdjiev, “Application of Monte Carlo method for tracking maneuvering target in clutter,” Math. Comput. Simul., vol. 1851, pp. 1–9, 2000. [37] A. Farina, B. Ristic, and D. Benvenuti, “Tracking a ballistic target: Comparison of several nonlinear filters,” IEEE Trans. Aerosp. Electron. Syst., vol. 38, no. 3, pp. 854–867, Jul. 2002. [38] Collaborative Adaptive Sensing of the Atmosphere, National Science Foundation Engineering Research Center. [Online]. Available: www. casa.umass.edu [39] D. McLaughlin, V. Chandrasekar, K. Droegemeier, S. Frasier, J. Kurose, F. Junyent, B. Philips, S. Cruz-Pol, and J. Colom, “Distributed collaborative adaptive sensing (DCAS) for improved detection, understanding, and prediction of atmospheric hazards,” in Proc. 9th Symp. IOAS-AOLS, Amer. Meteorol. Soc. Conf., San Diego, CA, 2005. [40] CASA, System Requirements Document, 2005. [41] American Society of Civil Engineers, Report Card for America’s Infrastructure, 2003. [Online]. Available: www.asce.org [42] Y. Zhai, M. Yeary, and J.-C. Noyer, “Target tracking in a sensor network based on particle filtering and energy-aware design,” in Proc. IEEEIMTC, Apr. 2006.

1441

Mark B. Yeary (S’94–M’99–SM’03) received the B.S. (honors), M.S., and Ph.D. degrees from the Department of Electrical Engineering, Texas A&M University (TAMU), College Station, in 1992, 1994, and 1999, respectively. As a student at TAMU, he was a Charter Member and an Officer of the Engineering Scholars Program and received the Dean’s Outstanding Student Award. In 1995, he was with IBM, Austin, TX, as a Member of a microprocessor development team. As a graduate student, he was a Teaching and Research Assistant, for which he received the Outstanding Teaching Assistant award from the local IEEE Student Chapter two years in a row, and was also nominated to be an NSF/FIE 1998 New Faculty Fellow. After his graduation in 1999, he continued to be a member of the DSP group and a Lecturer in the Department of Electrical Engineering, TAMU. Since the fall of 2002, he has been an Assistant Professor with the Department of Electrical and Computer Engineering, University of Oklahoma, Norman. His research, teaching, and consulting interests are in the areas of customized embedded DSP systems and digital signal processing as applied to radar signal processing, digital communications, image processing, adaptive filter design, and real-time systems. His applied signal processing contributions are many, which include the design of all-digital system-ona-chip scheme for a Ka band radar and various target tracking-algorithm developments. Dr. Yeary is a member of the Tau Beta Pi and Eta Kappa Nu honor societies. He has served on numerous national and international review panels, and several are noted here. He was a Session Chairman of the 2001 IEEE International Symposium on Intelligent Signal Processing and Communication Systems. In 2001, he served on the steering committee of the 2001 Future Energy Challenge. He has also served as a Session Chairman at 2003 IEEEInternational Multimedia Teleconferencing Consortium (IMTC) and as a Technical Committee Member of the 2003 and 2004 IEEE-ICIP. He is serving as a Technical Committee Member for IMTC in 2006 and will do so again in 2007.

Yan Zhai (S’03) received the B.S. degree from the Tsinghua University, Beijing, China, in 1998 and the M.S. degree from the Department of Mechanical Engineering (systems and control), Oklahoma State University, Stillwater, in 2002. He is currently working toward the Ph.D. degree in electrical engineering at the School of Electrical and Computer Engineering, the University of Oklahoma, Norman. His primary research emphasis is in the area of stochastic signal processing with a focus on developing sequential importance sampling (SIS) [particle filtering (PF)] algorithms and applying them to nonlinear state-estimation applications, such as image and video processing. He is currently working as a Graduate Research Assistant and Teaching Assistant. Mr. Zhai is a member of the Tau Beta Pi and Phi Kappa Phi honor societies. He has written a variety of papers, and he has served as a reviewer for several national and international conferences. He is a recent recipient of the Sooner Heritage Scholarship at University of Oklahoma.

Tian-You Yu is an Assistant Professor with the School of Electrical and Computer Engineering, University of Oklahoma, Norman. His education at the University of Nebraska and his postdoc experience at the National Center for Atmospheric Research, Boulder, CO, provide a unique crossdisciplinary background of atmospheric research. He has many reviewed technical journal and conference papers in the areas of applications of signal processing techniques to radar problems and studies using atmospheric radars. In parallel with his technical strength, he has a passion for delivering high-quality education. He has developed and taught several undergraduate and graduate courses at the University of Oklahoma.

1442


Shamim Nematifar (S’04) received the Bachelor’s degree in mathematics from the University of Oklahoma, Norman, in May 2005. He is currently working toward the master’s degree in signal processing, computational, and applied mathematics. His research interests include machineintelligence-based systems, signal and image processing, and system identification and modeling. Mr. Nematifar is a winner of the IEEE Region 5 Student Paper Contest in 2005.

Alan Shapiro received the B.S. degree (with distinction) in atmospheric science from Cornell University, Ithaca, NY, in 1983 and the M.A. and Ph.D. degrees in geophysical fluid dynamics from The Johns Hopkins University, Baltimore, MD, in 1985 and 1987, respectively. From 1987 to 1989, he was a Postdoctoral Scientist at the National Meteorological Center (now one of the National Centers for Environmental Prediction), where he worked on the initialization of numerical hurricane-prediction models. From 1990 to 1995, he was a member of the research team at the Center for Analysis and Prediction of Storms (CAPS), where he worked on problems of radar data analysis, single-Doppler-velocity retrieval, and numerical weather prediction. He joined the faculty of the School of Meteorology, University of Oklahoma, Norman, in 1996, where he is currently an Associate Professor. His research interests range from the dynamics of natural convective flows, waves, vortices, and katabatic flows to the development of single- and multiple-Doppler radaranalysis techniques and tornado-detection algorithms. Dr. Shapiro is a member of the American Association for the Advancement of Science, American Geophysical Union, American Mathematical Society, American Meteorological Society, and the Society for Industrial and Applied Mathematics.