A Multiscale Time-Varying Approach to Moving Source Tracking

A Multiscale Time-Varying Approach to Moving Source Tracking

Copyright © IFAC Adaptive Systems in Control and Signal Processing. Grenoble. France. 1992 A MULTISCALE TIME-VARYING APPROACH TO MOVING SOURCE TRACKI...

1MB Sizes 3 Downloads 118 Views

Copyright © IFAC Adaptive Systems in Control and Signal Processing. Grenoble. France. 1992

A MULTISCALE TIME-VARYING APPROACH TO MOVING SOURCE TRACKING o. Zugmeyer and J.P. Le Cadre IRISAlCNRS. Campus de BeaU/itu. 35042 ReMes Cedex. France

Abstract: After a simplified model of the spatial frequencies of a moving source has been defined. spatio-temporal methods including the (unknown) motion model in the source one are derived.This study is situated in the general context of passive processing and deals with the estimation of the sourco motion parameters in order to enhance detection and tracking of weak moving sources. The proposed methods are based. for a large part. on multiscale analysis of a state-space modelling of the data. In order to achieve a coherent fusion of the state estimators. the interpolation procedures constitute a basic tool.The interpolation may be considered either on the observations or on the state vector

1

means of markovian models . The smoothness constraints are related with the hyperparameters defining the state of these markovian models. Indeed, any real source trajectory, even a close one, can be described by such a model. This model is then related to a time varying temporal AR model of the spatio-temporal file {y(i, m)} . The time varying parameters of this AR model is then estimated by means of Kalman filtering and smoothing. Afterwards, the method is extended to the multidimensional (multifile) analysis. These extensions rely heavily on interpolation procedures. Furthermore, this formalism allows us to combine the advantages of the spatial processings with those of the multiscale time-varying analysis. The paper is organized as follows : in a first time, a time-varying AR model of each file is considered as well as the estimation of the paramters defining this model. For that purpose, the Kalman filtering plays a major role. The hyperparameters are then estimated by maximizing a likelihood functional. This approach is then extended to multi file estimation and more generally to spatio-temporal analysis.

Introduction

For numerous applications (e.g. passive sonar) large integration times are required in order to improve the performances of array processors in terms of detection. angular resolution. tracking .. . It is then necessary to incorporate the (unknown) source motion models into the source ones and to define array processing using basically this extended source model. For that purpose. the following steps must then be considered :

2

• derivation of a simplified model of source motion, • incorporation of the motion model into the source one,

For the rest of the paper, we shall assume that the sensor array is linear and constituted of p equispaced sensors. The angle 0 denotes the bearing. defined w.r.t . the array axis. the parameter>' is the wavelength. Then the spatial frequency is defined by : (1) k = cos 0/>.

• spatio-temporal analysis (source motion analysis) . Classical array processing arc based upon a short time analysis itself followed by steps of source tracking, target motion analysis and data assocation . All these methods are based upon the maximization of a spatial contrast functional. However, for numerous practical applications (e.g. passive sonar), a large amount of spatio-temporal data is available. It seems, then, possible to use all these data in order to separate, detect and track the sources by using their respective trajectories. In this spatio-temporal approach, the concept of source trajectory, and more precisely its modelling, replaces the instantaneous source bearings in the (classical) spatial analysis. For that purpose, the basic analysis is monodimensional and closely, related to time series analysis . However, using a unique formalism, the basic method can be extended to the following steps: • multidimensional estimation of the parameters

A spatio-temporal model of the data

For a moving source, the following temporal model of kt will be considered : (2)

Note that eq. 2 amounts to replace locally the exact value of kt+1 by its first order expansion . As the temporal changes of k t are usually very small, this approximation is (locally) quite acceptable. Furthermore, our main objective consists in estimating the parameters kt associated with the various sources. The spatial frequencies corresponding to s sources are described by the following model :

k, k •. t +1

• multiscale coherent analysis,

=

k•. t

+ k•.t (3)

• multifrequency analysis, • spatio-temporal analysis (after beamforming),

k •.t +.

• optimization of the spatio-temporal analysis.

=

k•. t

+ k.,t

The vectors X t represent the consecutive array snapshots at a given frequency (omitted). Then, the estimated CSM

The time varying spatial frequencies are, then, modelled by

481

(Cross Spectral Matrix)

Note that rank

(hI)

RI

is defined as follows:

k

~ XIX;

lE (w(t, m)) = 0

is equal to I. Generally, this is not

k

is orthogonall.!J projected on the a Toeplitz matrix. Then Toeplitz subspace by means of the classical algorithm (this matrix is denoted iZI ,T)'

with:

tr: trace of a matrix, Rim, ~ zm, . RI t: transposi tion,

and: and:

• ' Z·.) ri(t) = -1-. 17' ( Ill

P-l

tr : denoting the trace of a matrix and Z the shift matrix. Now, consider the time serie {i'm(t)L (m fixed), then the sin ratio can be improved by using a low-pass filtering since, the value of k is usually very small. The filtered time-series will be called {y(t,m)L ,m and are represented by figure I. Y 11. I)

y (t + I I)

Y (t. 2)

Y (I + 1.2)

Y (I. P .1)

Y (I + I. p. I)

Y(I+N.I) _

'file I'

Y(I+N.2) _

'filer

RI: exact covariance matrix of the data y(t, m) (6) at the instant t Note that all the parameters correspond to a narrow band analysis. Our main problem consists now in estimating the time-varying transition matrix Ft". At this point, it is worth noting that, thanks to the Cayley-Halmilton theorem, there exist s scalar coefficients such that the following equality holds: ,

Ft":,'

y (t + N. p . I r--- "file p. I'

=

L aj(t, m) Ft":.(,-j)

(7)

j=1

I

Consequently, denoting Z(t, m) the noise free output of the linear system (eq. 5), i.e. :

Figure 1: General organization of the data files {y(t, In)} At this point, it is worth not.ing the exact expression of )'m (t), i.e. :

Z(t, m) ~ h" X(t, m) eq. 5 yields directly:

, Z(t

o ::; In o}

::; P - 1 , 8m = 0 for m

'#

+ s, m) = L

aj(t, m)Z(t + s - j, m)

(8)

j=t

0 (1::; m ::; P - 1)

The eigenvalues of the transition matrix Ft":. are also the roots of the associated polynomial p.,m(Z) defined as follows: ,

and b2 source (I'csp. noise) power spectral densities

(4) According to eqs (3) and (1), the data can be locally described by a state-space model defined as follows:

P"m(z) = zs - Laj(t,m)z,-j

(9)

j=1

Using eqs. (5) and (8), the following equality holds:

X(t+l,m)

,

Ft":.X(t,m)

y(t, m) = L aj(t, m)y(t - j, m) X(t,m+l)

Fo,IX(t,m)

y(t,m)

h" X(t, m)

with :

+ w(t, m)

, n(t, m) = w(t, m) - L aj(t, m)w(t - j, m)

with :

For the rest of this paper, the above model is assumed to be valid and we shall now consider stochastic modelling of the time varying variations of the spatial frequencies. for that purpose, the approach of Kitagawa and Gersh [3J is instrumental. More precisely, define /((t) as the vector constituted of all the kj(t) and Am(t) the vector formed of the coefficients aj(t,m), i.e. :

Fo,. ~ diag (exp(2irrdk t (t)), .. ·,exp(2irrdk,(t)))

d

(10)

j=t

Ft. I ~ diag (exp (Zirrdkl(t)), ... ,exp (Zirrdk,(t)))

h"

+ n(t, m)

j=1

~

(1,1, .. ·,1) m,#O intersensor distance

(5) The above equation is fundamental for the sequel. The transition matrix Ft,l is associated with the instantaneous variations of the spatial frequencies {k i ( t)} :=t. Therefore the parameters {ki(t)} are assumed to take constant values )n all the elementary intervals 6.T ~ (t + 1) - t. The transi;ion matrix Fo,. is associated with the instantaneous spatial ~requencies {ki,d . finally, the estimation noise (of the data I(t, m)) is w(t, m) . The study of its statistical properties is 'undamental [1], [2], t.he corresponding results are summa, 'izcd below :

[((t) ~

(k

~(t)

[~I(t), .. ·, ~,(t)J'

~

l

(t),i.: 2 (t), ... ,k,(t))1

(11)

The components of the vector ~(t) are independent, cen, tered, gaussian random processes. Therefore, the vector ~ is gaussian, i.e. : ~ is N(O, r) with: (12)

482

Obviously, higher values of n may be considered, but for practical applications an order 2 (or 3) seems quite sufficient even for very close (to the array) sources. Furthermore, for numerous applications, the matrices J" J'_I, J'-2 are rather similar so that the eqs. (15) and (16) simply result in ;

The variances {Tn are unknown. They are called hyper parameters [3] since they model the smoothness of the temporal changes of the parameters kj(I). Finally, the "gradient" operator is defined recursi vel)' as follows ;

yrO[((I)

• n • n

/\"(t)

li

yrl[((I) ~

yro[((I)-yrOI((t-l)

yr2[((I) ~

yrl[((t)-yrl[((t-l)

Assuming that the component of the vector Am(t) are differentiable functions of the vector k (t) is defined as usually:

(. ") ',m I,)

li

=

A(t) = 2A(t - 1) - A(t - 2)

=3

A(I) = :JA(t - 1) - 3A(t - 2)

3)

+ J,(t)

The time varying-AR (15) model of the data y(l ,m) is defined conditionnaly on the following parameters, which according to Gersh [3]' will be called hyperparameters ;

(13)

J

+ J,(I) + A(I -

=2

• the order n of the difference equation, • the input variances

m)) 1 SI,).. Ss

T?

(eq. 15)

(Da;(I, akj(t)

• the variance a~(I) of the (colored) noise n(l,m).

yielding;

Generally, these parameters are unknown. As it will be seen later, they can be obtained by maximizing a likelihood functional ". However , for most of the practical situations, their values seem to not too critical.

(11) (the symbol U denotes the differential) Consider the fir st order difference approxi mat ion of D11", l I), I.e . :

3

then the eq. 13 yields direct.ly :

We are now concerned with the es timation of the vector O(t), defin ed by eq. 15. For that purpose, the Kalman filter is the basic tool and the basi c equations take the following form in this case ['1] :

V'I , I ... (I)

= J"m yrl k(l)

More generally, if a n - th order difference equation is considered for k(I) , then this modelling y ields directly for various values of n (the index m is omitted for that part and the symbol Id stands for the identity matrix) ;

Estimation of the time-varying AR model

[((i) = P(III - 1)·1/"(1)· [H(t). P(III - J)/i"(t) O(llt)

n = 1:

A(l) n

(Id

+ J,J,__\)

I)

+ [((I) · [y(t,m)

- H(l)O(tll - I)]

P(tll) = (Id - I«I)U(I)) P(llt - 1)

A(I -1) +J,cf>(I)

= 2: A(I)

= O(llt -

+ a~(I)]

A(I - I) - J,J,__ \ A(t - 2)

O(llt - 1) = F (I)O( III)

+ J , rt>(I)

n = 3;

+ lit) = F(I)P(llt)T(I) + r'(I)

P(I

A(I)

with:

n = 3;

-I) 11(t - 1) (Id + 21,J,_,

A(I)

1

+J,J'__ 2 A (t - 3)

(2J,J'__1+ J,J'__') A(I 1

2

+ J.cI>(l)

(15)

Consequently, according to eqs. 5,14, the data {y(l, m) can be modelled by the following state-space model:

l.

r'~GrG·

2)

and the following notations: I

f( (I);

Kalm a n gain matrix at time I

P(III): Extrapolated error covariance matrix

J Ott)

F(t)O(l - 1)

+ G(I)(I) 0(111) ; Estimated state at time I

1

llm (I) O(l)

y(t,m)

+ a(t , m)

O( I + Ill):

Extrapolated estimate of the state.

with:

Ott)

~ (A'(I) IA'(I -

unitialization of the Kalman filter

1)1· · ·IA'(t - n + 1)1)'

[((010 ) = [0, ... , 0],

(n - th order differe nce equation)

lI(t) ~ (y(t - 1, TJI)," y(1 - s, m), 101', ...

,101')

(16)

0(010)

= 1 F(I) = Id, • n = 2 •

P(OIO) = >.Id(>. "great") i(O)

71

Id

G(l)

= J,

+ J,J,__\ Id

= Id

Furthermore, note that the above algorithm requires the calculation of the Jacobian matrix J,. This filter can be associated with a smoother. The previous approach is limited to the monodimensional case i.e. the analysis of the data {y(t,m)}, (m fixed) . The basic idea will be extended to the multidimensional analysis

-J,J,__\

F(t) = (

= (A'(O) IA'(O)I" 'IA'(O))' (18)

The matrices F and G are straightforwardly deduced from cq. 11 , i.c. ;

o 483

by using a multiscale approach. The proposed methods relv

heavily on an interpolation procedure, presented below. By using eq. 5, t he following interpolation formulae are valid for the state of the "fil m (i.e. {y(t,m)L):

X

A symmetrical way consists ill interpolating the state. For that purpose, define an interpolated state ott), then one has :

0(')

'l'l~erefore, with the notations of eq. 8, the interpolated data Z(t, m) satisfy:

Z

~

(I + I:o,m) = Ft·X(t,m)

(I + ':>, m) = h*F(m·.'{(t,m)

and , consequently, there exist coefficients

aj

[

::::2) ) 0(t,2) ;, "l.t,d to the "fil," 2.

and :

(t, mol s.t. :

Y,

= lI(t) . Ott) + N(t) (22)

Ott)

(19)

(t j

(t, mol

- ~ j, m)

Y (t

1)

+ G,~(t)

Then, one has with classical notations [41 :

8,(+) = 0,(-) + k,· [Y(t) - H,·

So that finally the interpolated data Y(t, m) can be modelled by time-varying A R-models with a common set of coellicients {aj (t, 1I1u)) ; =1 whatever the spatial index TIl, i.e. :

y(t, m) = 2::j=1

= F(t)O(t -

0,(-)]

(}, = (ldlO) 0, = id· 0, and considering the covariance matrix P of the prediction error c( t) ; one obtai ns directly:

+ 'i(t, m) (20)

P,( +)

= (id - [(,H,) j>,( -) (id - [(,H,) + k,R,k;

with:

R,

The model is rdcrcnced w.r.t . the spatial index moo The eITect of the data interpolation is illustrated by the figure 2, below. y(t- l,mol

y(t, mol

Y(I ---.1!!!L , 0+ 2) mo+2

~ cov

(N,)

and

k,:

Kalman gain

(23)

"file" mo

_

_

_

-

-

- "file" mo + 1

-

-

-

"fUe"mo+2

y(l,mo+ J)

The markovian model (15) is then extended to the multidimensional case hy using the eq. 19, i.e. :

I

0(1) =

1

~(t)O(t -

Y(t) = lI(t)O(t)

1)

+ G,~(t)

Figure 2: Temporal interpolation of the data "files" y(t, m), for an AR (2), 3 files . The problem consists now in estimating the Kalman gain which minimizes the mean square error TT (P,( +I). The calculation is direct and gives:

+ ri(t)

k,

with: (0, F, G) defined as in eq. 15

=

P,(+) =

idP,(-)ll, · (II,P,(-)II;

+ R,)

(id-[(,il,)P,(-)fd"

(24)

Whatever the considered approach, the last point to be considered is the estimation of the hyperparameters n (order of the difference model) and r2/(j2. More precisely, one obtains the following expression of the likelihood functional L.

y' Y(t,m)

L (Jl, n/ {Y}) = P (Yd n~2 P (Y, IY1 , Y2 ,'"

(y (t, mo),'" ,y (t, md)

mo))' t - -;; , "' , Y-(t - s~ (}-. (mo)

,

Y,-I)

(21 )

furthermore, one has:

The covariance lIlatrix of the noise is deduced fmlll eq. 6. Practically, the data interpolation procedure is replaced by the procedure described below. The scalar Kalman filter equations (17) are replaced by vectorial equations. Since they are quite similar, they are omitted. It is worth noting that the interpolation procedure of the data is not critical thanks to the large number of temporal data on each file.

-I

P (Y, IYI' Y2,"', Y.-I) = c · det

(re (,)) exp (-c*(t)r~:)2{t))

c(t) bcing the innovation vector at the time t, i.e. :

r«,)=li,p(tlt-l)Ji;+(j2Jd

(25)

The above expression allow us to calculate thc likelihood functional conditionnally to Jl.

484

4

Simulation results

5

The methods exposed above will be now illustrated by simulation results. For ~ ll this section , the source trajectori es are fix ed a nd correspond to it linear motion for the two sources. The temporal evolut ions of the associated spatial frequ encies are represented in figure 4 in terms of spatial freque ncies. The snapshot numb0r is plotted on the x axis from 0 to 3600. The elementary inte rval be tween two snapshots is 6.T = I sec .. The iustantaneous values of k,(t = 1 to 3600) for the two sources are represented in figure 5, they a re far to be constant. The applicatioll of the I
are represented on the figure 6 by solid lines while the estimated values are plotted on the das hed and dotted lines. The algorithm is, in this case, the multidimensional one. The array is linear, constituted of 100 equispaced sensors. The considered valu es of m are comprised from m = 40 to 80. In order to judge th e importance o f the hy pe rpara meter J.I. (J.I. = (72 / T2) , th e res ults of the mult idimensional Kalman filte r (section 3) have been represented for various of the parameter J.I. (J.I. = 10-8,10-12 , 10-16). The better results are obtained for J.I. = 10- 12 which see ms the optimal value of the hype rparameter. Thi s is confirmed by the analysis of the likelihood functional , however its maximum is rath er rough . In order to illustrate the benefits of the multidimensional anal ysis, these result.s are compared with the results obtained on various data files {y(t, m)} (m = 40,50,80), presented on figure 7. One can see that the filters converge more slowly, th erefore the result s a re clearly enhanced by using th e multi pie analysis.

t+ !

y mo

ymo

t+

2

t+ 0

Conclusion

The aim of this p~pcr is to present an original and general framework for source tracking. On the opposite to class ical methods which rel y upon an updating of the spatial model, the source motion model have been included into the source one. The utilization of the Kalman filter is then quite natural. The modelling of any source trajectory by a markovian model is quite satisfying. The proposed methods can be applied to close and fast sources, it seems thus promissing for enhancing the source tracking.

References [IJ J .P. LE CADRE e t O. ZUGMEYER, "Sur les propric tes statistiques des covariances estimces apres rectification, appli cations". Traitement du Signal, vol. 9, nQ I, mars 1992.

[2J O . ZUGMEYER and J.P. LE CADRE, " On the optimization of spatio-temporal analysis for source motion estimation ". ICASSP 92, San Francisco, mars 1992.

[3J G. KITAGAWA a nd W. GERSII , "A smoothness pri ors time-varying AR coeffi cie nt modeling of non-stationary covariance time seri es". IEEE Trans . on Automatic Control, vol. AC-30, nQ 1, p. 48-56, janvier 198.5.

[4J A. GELD , "Applied Optimal Estimation ". A. Gelb Edt., MIT Press, 1989.

t+ 0 +!

ymO

filemo+!

filemo+c

Kt.mo.c Figure 3: Jumping interpolation procedure

485

kt estimate 15~---------=-----'-=-----~---' m=4O .... 80 mu= le-S

. ~ ,.

11.5 11

·0.5 JQUO

01)

·2

: f--; ,

(J

500

_ _ _

..3. ::,._ .. :,~"":..~.~._::,:,,

1000

1500

2000

3500

JOOO

2500

4000

kt estimate m =40.... 80

:J '0

-------

----

_--~ \ourte

211L - - · _ _ _ _ _ _ _ _~_ _~_ , /lOO

1500

~UO(l

mu= le - 12(oplimaJ)

10

source I

-,

50n

_0., - . ... .. .. - r _ _

~---------------------

IS .II.IO.()

--'.

11

-...

,I

SOURCETRAJECTO~

JI)

CO

~ .

--..

'-'.

.... -.----~ .- ..----

2

_,_~--~---.J

~500

3000

J(HI

35011

·5

0

Figure 4: Temporal evo lutions of the source spatial frequ encies (two sources)

500

1000

1500

2500

3000

J~OO

4000

kt estimate

l"'i :\ 10.#0 ((I

2000

m =40 .... 80 mu = le- 16

",:

" -..

.In" IJ

SPATIAL FREQUENCY RATES

.. -,

"

...

~ -:.

:,-:- -

U

·5 0

500

1000

1500

2000

2500

3000

3500

4000

HI

Figure 6: Ml!ltidilllen~ional estimation of'the spatial frequency rates kl(t) and ~'2(t) for two sources at - 15 dB, m = 16 40 to 80, It = 10- 8 , 10- 12 (opt imal), 10-

o ·2

~oun;e

I

-----_.-

-40'---.~ "-"---1-"-"- - ~-"--~I-"-'---~-5~ I)U--J~ UUO---J5-UO---J-'OO()

Figure 5: Temporal evolutions of the spatial frequency rates kl (t) and k1(t) for the trajectories of figure 4, xlO's

3'~~~--~--~--~m~=~~~----~--~~

: V\"~------".- _____=.~.::::~~:':

__._. __

..............•...•..

-I

-20;----S~00~--~1~000~--7IS~00~~2~000~--~~~00~--~3~000~---3S~00----4~000 3 xlO-s m=SO

mu

=le-I I (optimal)

3 xlO·s

Figure 7: M~nodilllensi,onal estimation of the spatial frequency rates k\(t) and k1(t) for two sources at • 15 dO for various values of the spatial index m = 40, 50, 80

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

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

-I O~---;;;:;;---~:----;-f,;~-:::::=-~=:----:-::'-:-:-==::.:::....----.l SOO

486

1000

1500

2000

~OO

3000

3500

4000