Source Localization and Sensing A Nonparametric Iterative Adaptive Approach Based on Weighted Least -hec.pdf

I. INTRODUCTION

Source Localization and

Sensing: A Nonparametric

Iterative Adaptive Approach

BasedonWeightedLeast

Squares

The goal of array processing is to estimate the

locations and waveforms of sources by combining

the received data from multiple sensors so that the

desired signal is enhanced, while the unwanted

signals, such as interference and noise, are suppressed.

In active sensing applications such as radar/sonar

range-Doppler imaging, the aim is to find targets

present in a region of interest. For channel estimation

in communications, the aim is to estimate the non-zero

channel taps and their Doppler shifts, which are then

fed to the subsequent equalizer for symbol detection.

Low numbers of snapshots and low signal-to-noise

ratios (SNR) are among the many challenges that

array processing systems frequently face. Another

challenge is the presence of nearby sources, in terms

of location or Doppler, since closely spaced sources

are harder to discriminate. Herein both passive array

processing and active sensing are considered.

The most basic approach to array processing

is the classical delay-and-sum (DAS) method,

in which the received signal from each sensor is

weighted and delayed so as to focus on different

points in space. However, this method suffers from

low resolution and high sidelobe levels. There is a

vast amount of literature on methods that provide

superior performance over the DAS approach when

certain assumptions are met [1]. The well-known

standard Capon beamformer (SCB) [2] and multiple

signal classification (MUSIC) [3, 4] methods provide

superresolution when the sources are uncorrelated and

the number of snapshots is high. Many extensions

to these methods have been proposed to deal with

modelling errors, such as steering vector mismatches.

(See, e.g., [5]—[9].) However, none of these methods

is able to cope with very low snapshot numbers,

coherent or highly correlated sources, or severe noise.

Only a few snapshots are available when the

environment being sensed by the array is stationary

for a short duration of time. Moreover, to avoid

smearing, i.e., losing resolution because of wide

main beamwidths in the array response, averaging

can only be done over a small bandwidth [10, 11].

Therefore the number of available snapshots, which is

directly related to the time-bandwidth product, can be

very small, sometimes as small as 3, for applications

such as underwater array processing. Furthermore,

as discussed later, the data models for single-input

single-output (SISO) radar/sonar range-Doppler

imaging and multi-input single-output (MISO) channel

estimation in communication problems are similar

to the model used in array processing with a single

snapshot and an arbitrary array geometry (dictated by

the probing waveforms).

The array processing problem has been carried into

the sparse signal representation area by noticing that

the number of actual sources is usually much smaller

TARIK YARDIBI, Student Member, IEEE

JIAN LI, Fellow, IEEE

University of Florida

PETRE STOICA, Fellow, IEEE

Uppsala University

Sweden

MING XUE, Student Member, IEEE

University of Florida

ARTHUR B. BAGGEROER, Fellow, IEEE

Massachusetts Institute of Technology

Array processing is widely used in sensing applications

for estimating the locations and waveforms of the sources in

a given field. In the absence of a large number of snapshots,

which is the case in numerous practical applications, such

as underwater array processing, it becomes challenging to

estimate the source parameters accurately. This paper presents

a nonparametric and hyperparameter, free-weighted, least

squares-based iterative adaptive approach for amplitude and

phase estimation (IAA-APES) in array processing. IAA-APES can

work well with few snapshots (even one), uncorrelated, partially

correlated, and coherent sources, and arbitrary array geometries.

IAA-APES is extended to give sparse results via a model-order

selection tool, the Bayesian information criterion (BIC). Moreover,

it is shown that further improvements in resolution and accuracy

can be achieved by applying the parametric relaxation-based

cyclic approach (RELAX) to refine the IAA-APES&BIC

estimates if desired. IAA-APES can also be applied to active

sensing applications, including single-input single-output

(SISO) radar/sonar range-Doppler imaging and multi-input

single-output (MISO) channel estimation for communications.

Simulation results are presented to evaluate the performance of

IAA-APES for all of these applications, and IAA-APES is shown

to outperform a number of existing approaches.

Manuscript received December 11, 2007; revised July 21 and

November 4, 2008; released for publication November 16, 2008.

IEEE Log No. T-AES/46/1/935952.

Refereeing of this contribution was handled by W. Koch.

This work was supported in part by the Office of Naval Research

(ONR) under Grants N00014-07-1-0193, N00014-07-1-0293,

and N00014-01-1-0257, the Army Research Office (ARO)

under Grant W911NF-07-1-0450, the National Aeronautics and

Space Administration (NASA) under Grant NNX07AO15A, the

National Science Foundation (NSF) under Grants CCF-0634786,

ECS-0621879, and ECS-0729727, the Swedish Research Council

(VR), and the European Research Council (ERC).

Opinions, interpretations, conclusions, and recommendations are

those of the authors and are not necessarily endorsed by the United

States Government.

Authors’ addresses: T. Yardibi, J. Li, M. Xue, Dept. of Electrical

and Computer Engineering, University of Florida, PO Box 116130,

Gainesville, FL 32611, E-mail: (li@dsp.ufl.edu); P. Stoica, Dept.

of Information Technology, Uppsala University, Uppsala, Sweden;

A. B. Baggeroer, Depts. of Mechanical Engineering & Electrical

Engineering and Computer Science, Massachusetts Institute of

Technology, C ambridge, MA 02139.

0018-9251/10/$26.00 ° 2010 IEEE

IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 46, NO. 1 JANUARY 2010

425

AUTHORIZED LICENSED USE LIMITED TO: IEEE XPLORE. DOWNLOADED ON MAY 13,2010 AT 11:53:44 UTC FROM IEEE XPLORE. RESTRICTIONS APPLY.

TABLE I

Notation used in the Text

k¢k 0 ` 0 -norm

k¢k 1 ` 1 -norm

k¢k 2 ` 2 -norm

k¢k F Frobenius norm of a matrix

¯ the Hadamard (elementwise) matrix product

tr( ¢ ) trace of a matrix

( ¢ ) T

transpose of a vector or matrix

( ¢ ) H

conjugate transpose of a vector or matrix

Fig. 1. Far-field linear array.

than the number of potential source points that can

be considered [12—29]. Sparsity-based techniques

have also been used in spectral estimation [30], image

processing, and array design [31—33], among many

other application areas. Sparse signal representation

algorithms can deal with a few snapshots (even one).

However, they may require large computation times

and the fine tuning of one or more hyperparameters.

Application areas, ranging from passive source

localization to active radar/sonar range-Doppler

imaging and channel estimation for communications

are addressed here. Section II considers passive

sensing applications and describes a weighted

least squares-based iterative adaptive approach

for amplitude and phase estimation (IAA-APES).

The algorithm is named IAA-APES herein since

its derivation resembles that of the amplitude and

phase estimation (APES) algorithm [34—36]. (See

Section IIC.) This name also distinguishes it from the

maximum likelihood (ML)-based iterative adaptive

approach (IAA-ML) discussed in the Appendix to

further motivate IAA-APES. Using the Bayesian

information criterion (BIC) [37, 38], IAA-APES can

be extended to yield point source estimates. (This

approach is referred to as IAA-APES&BIC.) Next, a

parametric relaxation-based cyclic approach, namely

RELAX [39, 40], is discussed as a way to further

refine the results of IAA-APES&BIC. (This approach

is referred to as IAA-APES&RELAX.) Section III

presents the data models for the aforementioned active

sensing applications and emphasizes the similarities

to the passive sensing case. IAA-APES is evaluated

via comprehensive simulations in Section IV, and the

IAA-APES’ performance is compared with that of a

number of existing approaches.

Notation : We denote vectors and matrices by

boldface lowercase and boldface uppercase letters,

respectively. The k th component of a vector x is

written as x k .The k th diagonal element of a matrix

P is written as P k . See Table I for other symbols and

their meanings.

used to estimate the desired source characteristics.

This section first introduces the data model for such

applications and then the IAA-APES algorithm in this

model.

A. Data Model

Consider the wavefield generated by K sources

located at Μ ,where Μ =[ Μ 1 , Μ 2 , ::: , Μ K ]and Μ k 1 are the

location parameters of the k th signal, k =1, ::: , K .In

the narrowband, multi-snapshot case, the M£ 1 array

output vector of an M element array in the presence

of additive noise can be represented as [5, 43]

y ( n )= A ( Μ ) s ( n )+ e ( n ), n =1, ::: , N (1)

where N is the number of snapshots, A ( Μ )isthe

M£K steering matrix defined as A ( Μ )=[ a ( Μ 1 ),

a ( Μ 2 ), ::: , a ( Μ K )], and s ( n )=[ s 1 ( n ), s 2 ( n ), ::: , s K ( n )] T ,

n =1, ::: , N , is the source waveform vector at time n .

The array steering matrix has different expressions,

depending on the array geometry and on whether the

source is in the near-field or far-field of the array. For

instance, the steering vector corresponding to the k th

source for a far-field linear array (where Μ k represents

the impinging angle of source k in this case) is given

a ( Μ k )=[ e ¡j (2 ¼f=c 0 ) x 1 cos( Μ k ) , ::: , e ¡j (2 ¼f=c 0 ) x M cos( Μ k ) ] T

(2)

where f is the center frequency, c 0 is the wave

propagation velocity, and x m is the position of the

m th sensor, m =1, ::: , M . (See Fig. 1.) Note that for a

planar array or near-field sources, only the expressions

for the steering vectors have to be modified; the

algorithms that are presented can be applied without

any modifications, since a ( Μ ) is assumed to be a

known function of Μ . The number of sources, K ,is

usually unknown; hence, here, K is considered to be

the number of scanning points in the region. In other

words every point of a predefined grid that covers

the region of interest is considered as a potential

source whose power is estimated. Consequently,

II. PASSIVE SENSING

In passive sensing applications such as

aeroacoustic noise measurements [41] and underwater

acoustic measurements [42], an array of sensors is

1 With a slight abuse of notation, we do not use bold font for Μ k ,

k =1, ::: , K , which might be multi-dimensional, for simplicity.

426

IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 46, NO. 1 JANUARY 2010

AUTHORIZED LICENSED USE LIMITED TO: IEEE XPLORE. DOWNLOADED ON MAY 13,2010 AT 11:53:44 UTC FROM IEEE XPLORE. RESTRICTIONS APPLY.

K will be much larger than the actual number of

sources present, and only a few components of

f s ( n ) g will be non-zero. This is the main reason why

sparse algorithms can be used in array processing

applications.

[26]. Fuchs [27, 28] uses a sparsity-constrained

deconvolution approach that assumes the sources

are uncorrelated and that the number of snapshots

is large. The sparsity-constrained solution is

obtained with a LASSO or BP type of algorithm.

Reference [29] introduces two hyperparameter

free deconvolution algorithms exploiting sparsity:

a sparsity-based extension to the deconvolution

approach for the mapping of acoustic sources

(DAMAS) [45] (which is similar to [27] and [28]

andwidelyusedinpractice)andasparsitybased

covariance-matrix fitting approach. Extensions

to the correlated source case are also provided.

However, the methods in [29] are based on the

sample covariance matrix, and hence, these

methods do not work well with a limited number of

snapshots.

B. Related Work

We focus our attention on array processing

algorithms exploiting sparsity, which have

gained noticeable interest recently. Sparse signal

representation aims at finding the sparsest s ,such

that y = As is satisfied, i.e., to minimize k s k 0 ,such

that y = As ,where A is known and y is measured.

The problem in its original form is a combinatorial

problem and is nondeterministic polynomial-time

(NP) hard, making it impractical [12]. Fortunately,

when s is sufficiently sparse [12—15], k s k 0 can

be replaced by k s k 1 , which leads to a convex

optimization problem that can be solved much more

easily by using, for instance, the least absolute

shrinkage and selection operator (LASSO) [16] or

basis pursuit (BP) [17] algorithms. Alternatively, the

focal underdetermined system solution (FOCUSS)

algorithm [18], which is derived using Lagrange

multipliers, can be used to iteratively solve the

sparse problem. A Bayesian approach, such as

sparse Bayesian learning (SBL) [19, 20] or the

approach in [21], can also be used to estimate s .

The algorithms mentioned up to this point are for

the single-snapshot case only. The extensions of

FOCUSS and SBL to the multiple-snapshot case

are M-FOCUSS [22] and M-SBL [23], respectively.

Another algorithm, the ` 1 -SVD (singular value

decomposition) algorithm [24, 25] is similar

to BP or LASSO, but this algorithm can work

with multiple snapshots. (For the single-snapshot

case, ` 1 -SVD becomes a LASSO and BP type

of method.) M-FOCUSS requires the tuning of

two hyperparameters, which might affect the

performance of the algorithm significantly. ` 1 -SVD

requires the tuning of a hyperparameter and an

estimate for the number of sources. Moreover,

implementing ` 1 -SVD requires convex optimization

software, such as SeDuMi [44]. M-SBL does

not require any hyperparameters. However,

M-SBL converges quite slowly in its original form

[19, 20, 23].

Besides the above algorithms, which are the

focus of our attention in the numerical examples,

there are other sparsity based approaches worth

mentioning. Reference [26] adds an additional spatial

sparsity regularizing term (an ` 2 -norm constraint)

to the ` 1 -norm constraint, and it minimizes a cost

function similar to that of ` 1 -SVD. However, this

method has two hyperparameters, assumes that

the source waveforms can be represented by a

sparse basis, and has high computational complexity

C. IAA-APES

IAA-APES is a data-dependent, nonparametric

algorithm based on a weighted least squares (WLS)

approach. Let P be a K£K diagonal matrix, whose

diagonal contains the power at each angle on the

scanning grid. Then P can be expressed as

P k = 1

js k ( n ) j 2 , k =1, ::: , K: (3)

n =1

Furthermore, define the interference (signals at angles

other than the angle of current interest Μ k ) and noise

covariance matrix Q ( Μ k )tobe

Q ( Μ k )= R ¡P k a ( Μ k ) a H ( Μ k )

( 4 )

where R = A ( Μ ) PA H ( Μ ). Then the WLS cost function

is given by (see, e.g., [34]—[36] and [43])

k y ( n ) ¡s k ( n ) a ( Μ k ) k 2 Q ¡ 1 ( Μ k )

(5)

n =1

a H ( Μ k ) Q ¡ 1 ( Μ k ) a ( Μ k ) , n =1, ::: , N: (6)

This looks like the result that would be obtained by

employing APES [34—36], but it is actually different

than APES since APES obtains Q ( Μ k ) from the data

by forming subapertures, while IAA-APES computes

Q ( Μ k ), as in (4). Moreover, IAA-APES is iterative, but

APES is not, and APES cannot be used with arbitrary

array geometries.

Using (4) and the matrix inversion lemma, (6) can

be written as

ˆ s k ( n )= a H ( Μ k ) R ¡ 1 y ( n )

a H ( Μ k ) R ¡ 1 a ( Μ k ) , n =1, ::: , N: (7)

YARDIBI ET AL.: SOURCE LOCALIZATION AND SENSING: A NONPARAMETRIC ITERATIVE ADAPTIVE APPROACH 427

AUTHORIZED LICENSED USE LIMITED TO: IEEE XPLORE. DOWNLOADED ON MAY 13,2010 AT 11:53:44 UTC FROM IEEE XPLORE. RESTRICTIONS APPLY.

where k x k 2 Q ¡ 1 ( Μ k ) = x H Q ¡ 1 ( Μ k ) x and s k ( n ) represents the

signal waveform at angle Μ k and at time n . Minimizing

(5) with respect to s k ( n ), n =1, ::: , N , yields

ˆ s k ( n )= a H ( Μ k ) Q ¡ 1 ( Μ k ) y ( n )

TABLE II

The IAA-APES Algorithm

than necessary reduction of the number of degrees

of freedom when some of the interfering sources are

coherent since the cancellation of multiple coherent

interfering sources would require only one DOF if

the correct structure of P were known. However, we

do not assume that the true structure of P is known.

Moreover, it is the diagonal structure of P assumed by

IAA-APES that makes the algorithm work properly

even for low number of snapshot cases and coherent

sources.

( a H ( Μ k ) a ( Μ k )) 2 N

P k =

j a H ( Μ k ) y ( n ) j 2 , k =1, ::: , K

n =1

repeat

R = A ( Μ ) PA H ( Μ )

for k =1, ::: , K

ˆ s k ( n )= a H ( Μ k ) R ¡ 1 y ( n )

a H ( Μ k ) R ¡ 1 a ( Μ k ) , n =1, ::: , N

P k = 1

j ˆ s k ( n ) j 2

D. IAA-APES&BIC

n =1

end for

until (convergence)

In many applications it is desirable to obtain point

estimates rather than a continuous spatial estimate.

To achieve this sparsity we incorporate a model-order

selection tool, i.e., the BIC [37, 38], into IAA-APES.

Let P denote a set containing the indices of the peaks

selected from the IAA-APES spatial power spectrum

estimate. Also let I denote the set of the indices of

the peaks selected by the BIC algorithm so far. The

IAA-APES&BIC algorithm works as follows: first

the peak, from the set P , giving the minimum BIC is

selected. Then the second peak, from the set P¡I ,

which together with the first peak gives the minimum

BIC, is selected, and so on, until the BIC value

does not decrease anymore. 2 The IAA-APES&BIC

algorithm is summarized in Table III. BIC i ( ´ )is

calculated as follows (see [38])

TABLE III

The IAA-APES&BIC Algorithm

P : Set of peaks obtained from IAA-APES

I =Ø, ´ =1, quit=0, BIC old = 1

repeat

i 0 =argmin

i2P¡I

° 2

This avoids the computation of Q ¡ 1 ( Μ k ) for each

scanning point, i.e., K times. Moreover, f ˆ s k ( n ) g can

be computed in parallel for each scanning point,

which makes IAA-APES amenable to implementation

on parallel hardware. IAA-APES is summarized in

Table II. Since IAA-APES requires R ,whichitself

depends on the unknown signal powers, it has to

be implemented as an iteration. The initialization is

done by a standard DAS beamformer. Our empirical

experience is that IAA-APES does not provide

significant improvements in performance after

about 15 iterations. In IAA-APES, P and, hence

R , are obtained from the signal estimates of the

previous iteration and not from the snapshots as in

conventional adaptive beamforming algorithms, such

as SCB, which fails to work properly with coherent

or highly correlated sources, or few snapshots [43].

The Appendix provides an alternative derivation of

IAA-APES, based on the ML principle. IAA-APES

is shown to be an approximation of the IAA-ML

algorithm, which is locally convergent due to

cyclically maximizing the likelihood function. Hence,

the analysis in the Appendix provides an approximate

calculation that shows the local convergence of

IAA-APES.

Because P is assumed to be a diagonal matrix,

one degree of freedom (DOF) is needed to suppress

one interfering source. This may lead to a larger

BIC i ( ´ )=2 MN ln

y ( n ) ¡

a ( Μ j ) ˆ s j ( n )

n =1

j2fI[ig

+3 ´ ln(2 MN )

( 8 )

where ´ = jIj +1, jIj denotes the size of the set I , i is

the index of the current peak under consideration, and

f ˆ s j ( n ) g n =1 is the IAA-APES signal waveform estimate

corresponding to angle Μ j , j 2fI[ig . Note that the

second term on the right side of (8) does not matter

to peak selection; it matters only when (8) is used to

select the number of peaks to retain.

E. IAA-APES&RELAX

RELAX [39, 40] is a parametric cyclic algorithm

that requires an estimate of the number of sources in

the field. The results of IAA-APES&BIC can be used

to provide a good initial estimate for the last step of

RELAX. This approach is outlined in Table IV. The

RELAX iterations can be terminated when the norm

of the difference between two consecutive estimates

is smaller than a certain threshold (5 £ 10 ¡ 4 for the

examples considered herein). The maximization

2 Alternatively, the largest ´ peaks can be selected so that when

the ( ´ +1)st largest peak is added to I , the BIC value does not

decrease. This simplified version gives similar results to the one

described above in our numerical examples.

428

IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 46, NO. 1 JANUARY 2010

AUTHORIZED LICENSED USE LIMITED TO: IEEE XPLORE. DOWNLOADED ON MAY 13,2010 AT 11:53:44 UTC FROM IEEE XPLORE. RESTRICTIONS APPLY.

BIC i ( ´ )

if BIC i 0 ( ´ ) < BIC old

I = fI , i 0 g

BIC old =BIC i 0 ( ´ )

´ = ´ +1

else quit=1

until (quit=1)

TABLE IV

The IAA-APES&RELAX Algorithm

Μ 0 : Locations of the peaks obtained from IAA-APES&BIC

K 0 : Number of peaks obtained from IAA-APES&BIC

f ˆ s k ( n ) g : Corresponding waveforms obtained from

IAA-APES&BIC

repeat

for k =1, ::: , K 0

K 0

y k ( n )= y ( n ) ¡

a ( Μ 0 i ) ˆ s i ( n ), n =1, ::: , N

i =1

i6 = k

Μ 0 k =argmax Μ 0 2 R

j a H ( Μ 0 ) y k ( n ) j 2

n =1

M a H ( Μ 0 k ) y k ( n ), n =1, ::: , N

Fig. 2. Pulse compression for radar/sonar range-Doppler imaging.

end for

until (convergence)

step of the algorithm can be implemented without

much computational effort by means of a fine search

only around the peak values of the IAA-APES&BIC

result. 3 RELAX can be useful in estimating off-grid

sources accurately and for further improving the

IAA-APES waveform and angle estimates.

the Doppler interval of interest is divided into L bins,

and the Doppler frequency for the l th Doppler bin is

denoted as ! l . Then the M samples of the received

signal that is temporally aligned with the return from

the range bin of current interest r (see Fig. 2) can be

represented by

y r =

® r , l s ( ! l )+

® r + m , l J m s ( ! l )+ e r

III. ACTIVE SENSING

l =1

m = ¡M +1

m6 =0

l =1

(12)

for r =1, ::: , R ,where ® r , l (which is proportional to

the complex “voltage” radar-cross section (RCS)

of the corresponding target) denotes the complex

amplitude of the returned signal from the range bin

of current interest r and the l th Doppler bin f® m + r , l g

denote the complex amplitudes of the returned signals

from the adjacent range bins, and e r denotes the noise.

The reflections from nearby range bins are considered

to be clutter. The M£M shift matrix J m takes into

account the fact that the clutter returns from adjacent

range bins need different propagation times to reach

the radar/sonar receiver:

In active sensing applications, besides the

receiver, there are also one or more transmitters.

The radar/sonar range-Doppler imaging problem for

a SISO system is first investigated in this section.

Then the channel estimation problem for MISO

communications is discussed.

A. Range-Doppler Imaging

Pulse compression refers to the process of

transmitting a modulated pulse and then matched

filtering the returned signal, which arrives at the

antenna altered by complex coefficients that bear

target information [48].

1) Data Model : Consider a range-Doppler

imaging radar/sonar with a transmitted pulse

˜ s =[ ˜ s (1), ˜ s (2), ::: , ˜ s ( M )] T

Ã¡ m ¡!

0 ¢¢¢ 1 ¢¢¢ 0

. . . . . . . . . . . .

0 . . . 1

. . . . . . . . . . . .

0 ¢¢¢ 0 ¢¢¢ 0

(9)

5 = J ¡m (13)

where M is the pulse length. Let

s ( ! l )= ˜ s ¯ d ( ! l )

J m =

(10)

be the Doppler shifted signal, where

d ( ! l )=[1, e j! l , ::: , e j ( M¡ 1) ! l ] T , l =1, ::: , L

for m =0, ::: , M¡ 1. Equation (12) can be written as

(11)

y r = S ® r + e r

(14)

3 The fine search can be implemented efficiently using the fast

Fourier transform (FFT) [39, 40] or a derivative-free uphill search

method, such as the Nelder-Mead algorithm [46, 47]. Note that the

latter method is available in the MATLAB optimization toolbox

with the name of “fminsearch.”

where

S =[ s ( ! 1 ), s ( ! 2 ), ::: , s ( ! L )] (15)

® r =[ ® r ,1 , ® r ,2 , ::: , ® r , L ] T

(16)

YARDIBI ET AL.: SOURCE LOCALIZATION AND SENSING: A NONPARAMETRIC ITERATIVE ADAPTIVE APPROACH 429

AUTHORIZED LICENSED USE LIMITED TO: IEEE XPLORE. DOWNLOADED ON MAY 13,2010 AT 11:53:44 UTC FROM IEEE XPLORE. RESTRICTIONS APPLY.

ˆ s k ( n )= 1

M¡ 1

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika: