Spatial Statistics and Digital Image Analysis (1991)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

8 Spatial Signal-Processing in Raciars and T. T. Kadota AT&T Bell Laboratories S.1 Introduction Sonars Radars and sonars are used for detecting and tracking targets. The surveil- lance radars and sonars typically employ arrays of sensors (or radiators) placed on the ground or in the ocean. A planar array is used for sonar to detect seismic explosions ant! for radar to track distant flying objects, and a linear array is used for sonar to detect distant underwater objects. The data thus obtained are in the form of a set of time-series that are relatecl by the spatial configuration of the array. The task of the radar and sonar systems is to process these data to detect the signal transmitted from a target and estimate the signal parameters related to the target location and velocity. Typically, the data contain noise and interfering signals besides the target signal, and the "signal-processing" (processing of these data) requires suppression of the noise and interference and enhancement of the signal. The literature on spatial signal-processing is enormous: the IEEE Trans- actions on Acoustic, Speech and Signal Processing, the IEEE Transactions on Aerospace Electronics, and the Journa1t of the American Statistical As- sociation, just to name a few. Covering the entire area and providing an adequate survey is beyond! the scope of this report. Insteaci, by using a sim- ple example, we give a list of how statistical decision and estimation theory is used on this form of multidimensional data to derive a signal-processing algorithm and indicate how to extend the basic approach to more complex problems in reality. 147

148 S.2 Detection ant} Estimation Problem The example we have chosen is the signal-processing of underwater acoustic data for detecting a narrow-band signal transmitted from a distant source and determining the direction of its arrival. The signal detection problem is traditionally cast as a problem of testing a nub hypothesis Bo (signal absent) against an alternative hypothesis B~ (signal present) as follows (Helstrom, 1986~: It) dt + dwelt), drj(t) j=1,...,J;O<t<T, (Ho) (H1) (8.1) where T is an observation time interval, drafty is the (incremental) acoustic data recorded at the josh sensor at time t, Aft) is the acoustic signal received at the Ah sensor at time t, dwj~t) is the (incremental) background noise at the Ah sensor at time t (assumed to be white Gaussian, independent from sensor to sensor), and Aft) is the interference (or adclitional noise) at the Ah sensor at time t. We assume that the background noise w and the interference z are mutually independent. Adopting the Neyman-Pearson formulation, we obtain the likelihood ra- tio between the two hypotheses Ho and B~ (i.e., the Radon-Nikodym deriva- tive between the two probability measures induced by the random fields of (~.1~. It can be expressed as sj(t) At + zj(t) At ~ dwj(t), dPl cry _ Ez expits + z, r) - 2 ITS + Z~2] dPo En exp[(z,r) - Arid where (~.2) {drj(t), O < t < T. j = 1,. . ., J), {sj(t), O < t < T. j = 1,. , J}, {Zj(t), o ~ ~ ~ T. i = l, .. . ~ J} ~ and Ez denotes the expectation with respect to the probability distribution of the z process and the inner product and the norm are defined by (a, r) = J T ~ | sj(t)drj(t)' llsll2 j=l o = J T As, | sj2(t) dt j=1 o

149 z 41/~4 ~ - ~ 4~/ FIGURE 8.1: Linear array of uniformly spaced sensors. 8.3 Using Linear Arrays of Uniformly Spaced Sensors The most widely studied case is the one with a linear array of equally spaced sensors (see Figure 8.1) and a narrow-band planewave signal with a known carrier frequency, which may be represented by skits = Resects exp[-iw~t—rj)] = Reaj(?,b)sett) exp(-i~t), (~.3) where Re denotes the real part of what follows, {sect)' O < t < T) is the complex signal envelope function, and ~ = the carrier (angular) frequency such that AT ~ 2;r, Tj = (;— 1~4 cos ¢/c = the time delay of the signal planewave arrival at the Josh sensor relative to the first, = the distance between two adjacent sensors,

150 the velocity of sound propagation, = the assumed direction of the signal planewave arrival, and aj(~) = expti~j—1)—tacos ¢], j = 1, e ~ ~ ~ J. ~ (~3~4) Then, in the absence of the interference z, the logarithm of the likelihood ratio becomes log,dpl~r) = (s,r)—Use = Re ~a1~t,b) sets) drj~t) - 2 So , (8 5) where drj is the envelope function of the narrow-band representation of the data, namely, drift) = Re alrj(~) exp(-i~t) . Suppose a planewave arrives in the direction of jot Then the signal part of the data-dependent term, the first term of (~.5), is proportional to Re a(?~la(?~'O) = Re 1—exptit7~a?(cos To - cos ¢~/c] (~.6) 1 - exp[~wcl~cos To - cos ?,b)/c] where Ably = fad (I), . . ., a, is the direction (or steering) column vector in the Redirection and ~ denotes the complex conjugate transpose. Equation (~.6) is in the form of a "main beam" centered at ~ and "side lobes" on each side of the main beam as To varies from -~r/2 to ~r/2. Hence, the process- ing (of the data r) described by (8.6) is called "beam forming" (Steinberg, 1976~. On the other hand, if the actual direction of the planewave arrival To is fixed and the assumed direction fib is varied, (~.6) attains the maximum at ~ = To, namely, when the beam is steered at the signal source. Thus, detection of the planewave signal and estimation of its direction of arrival are done by varying ~,b (steering the beam) from—~r/2 to ~r/2 to find the maximum of (~.6) and comparing it to a preassigned threshold determined by the false-alarm probability (according to the Neyman-Pearson criterion) (Helstrom, 1986~. Instead of steering a single main beam, one can place many beams to fill the angular sector ~—~r/2, ~r/2) by providing many direc- tion vectors a(?,km), m = 1,...,M. Then by comparing the magnitude of (~.6) for each lbm, instead of varying fib, we effectively accomplish the same task of detection and estimation. In the presence of the interference z, some modification to the beam- forming is necessary. For any interference to be effective against the signal, it must have energy at or near the carrier frequency (otherwise, it can be

151 simply filtered out). For the sake of simplicity, suppose we have one sinu- soidal interference arriving in the ,,6' direction with a Gaussian distributed amplitude, namely, Aft) = Reuexp(-u~- (j- l~cicosl,b'/c]) = Reuzjexp(-i~t), (~.7) where the second equality defines z; and u is a complex Gaussian variable with mean O and variance i32. This represents a planewave arriving through a multipath medium causing the Rayleigh fading. By carrying out the expec- tation with respect to z in (~.2) (i.e., with respect to u), the data-dependent term of the log likelihood ratio becomes Re (s, r - EO{z~r)) = Re (Se, art—I + 22 ably) a(¢ babe ~ ~ r) J / = where Redid+ R)-is,r`), sj(~) = aj(~)seLt), and Eo~z~r) is the conditional expectation of z given r under Ho, and R = ,(32a(¢')a(¢')~. The first member of (~.~) has an obvious interpre- tation: the optimum processing is to make the least-mean-square-error es- timate of the interference and subtract the estimate from the data before the beamforming. The second member, on the other hand, shows how the conventional beamformer is to be modified due to the presence of the inter- ference. By recalling that a(¢') is the steering vector in the direction of the interference, the modified beamformer has a considerably reduced output in the direction of the interference, thus acquiring the term, null-steering (Steinberg, 1976; Gabriel, 1976; Friediander and Porat, 1989~. In practice, the interfering source is not known a priori and its covari- ance matrix R must be estimated. The estimation may be done beforehand or simultaneously with the detection operation, assuming that the direction of the signal arrival is known (which is the case with the fixed multibeam scheme). This simultaneous method is referred to as the adaptive beam- forming and is implemented by attaching a variable gain (or weight) and a variable time-delay (which are adjusted as data are obtained) to the out- put of each sensor. Of course, such an adjustment must be done rapidly so that accurate signal detection and direction-of-arrival estimation can be accomplished. The iterative methods of adjusting and their convergence characteristics have been extensively studied. Monzingo and Miller (1980)

152 is one of the comprehensive textbooks on the subject. A benchmark paper series edited by Haykin (1980) has many important papers, including those of Gabriel, Applebaum, Widrow, Griffiths, and Owsley. 3.4 Using General Arrays of Sensors The results presented above can be generalized as follows (Kadota en cl Rot main, 1977~: first, instead of a linear array of equally spaced sensors, we can consider a general three-dimensional array (or configuration) with the coordinates if;, A, (j ), j = 1, . . ., J. Then the planewave arriving in the (8, ~b)-direction, where ~ and i,) are the elevation and the azimuthal angles, incurs at the JO sensor the phase shift expressed by aj(w, 8, ¢) = exp pi—(A cos ~ cos ~ + 77j cos ~ sin ~ + (j sin 8~] . (~.9) Next, instead of a single planewave, we consider a signal consisting of Ii planewaves, each having a different frequency ~k, k = 1,...,](, and each arriving in M different directions (69m,?/)m), m = 1,. . . ,M. For convenience, we assume that cokT is an integral multiple of Or for every k. Also, rather than a "slowly varying" envelope function sets), we consider a complex Gaus- sian variable (independent of time) as the amplitude of each planewave. Thus, the signal at the JO sensor is now given by K M sj(t) = ~ ~ Re Ukmaj(Wk,(9m'?/)m) exp(—inlet), (8.10) k=1 m=1 where (nkm), ukm = ukm+inkm, are complex, zero-mean, Gaussian variables with E UkmUk~m' = E Ukmuk~m' = Pkk~mm, ~ k, k = 1, . . ., h; m, m' = 1, . . ., M . We allow some As to be zero since not all frequency components have all M arrival directions. The interference is now generalized to vj(t) = ~ Re vje exp (-i T ) , (8.11) where {v;e}, vj' = vje + iVje, are complex, zero-mean, Gaussian variables. That is, for each j, vj(t) is a-discretized version of the spectral representation

153 of a general stationary noise. Then the data-dependent part of the log- likelihood ratio takes the following quadratic form in the data: (x, (! + V)-~S(! + V + S)-ix), (~.12) where X — (X1 In,. ,X(K_l)J+l, ~XKJ, X1. Id,. ,X(K_l)J+l, ~XKJ)' X(~l)J+i = (T) | cos~ktdrj(t), ~(~l)J+j = (T) | sinwktdrj(t)' [ ~ A ] with the ((k—l)J + j' (k' - 1)J + j,)th elements of S and S given respectively by the real and the imaginary parts of M ~ Pkk~mm~ai(Wk'6tm'(m)a:~(wki~f~mi,lbmI), j,; =1~, J; k,k =1,,X, m,m'=1 and V=[O Vat (V) (k—l ) J+ j,(k'—l ) J+ jl = E Vjl Vj'tikk' = E vJevj~e~kkl, i_ ANT 21r The J sensors constitute spatial samplers of the available (acoustic) data and their configuration specifies the pattern of spatial sampling. This sampling pattern is incorporated into the covariance matrices S and V to influence the detection statistic (8.123 which specifies the data-processing algorithm. Although the linear array (with or without the equal spacing) is the most common configuration, due primarily to the ease of implementation, the sampling pattern can be considered as a factor with respect to which the detection and estimation performance can be optimized. In fact, we show next an interesting example of this application. So far, we have assumed that the sensor positions are rigidly fixed and their coordinates are known a priori. Although this is the case with the phased-array radars and seismographic sensors, for underwater-acoustic sen- sors the exact positions in the ocean are difficult to determine and calibra- tion of the array becomes necessary. One way to deal with this problem

154 is to mode} the deviation (or fluctuation) of the sensor position (from the presumed value) as an additional noise and incorporate into the optimum processor (for detection and estimation) the sensitivity of the performance to this noise. For example, the array gain (the output signal-to-noise ratio of an array processor for a given direction of signal arrival) may be maximizer! under the constraint that the array-gain sensitivity to the sensor-position noise be kept below a given level (Cox et al., 1987~. An alternative is to devise the sensor configuration so as to make these test statistics immune to the sensor-position fluctuation. The ESPRIT (Estimation of Signal Param- eters by Rotational Invariance Techniques) method (Roy and Kailath, 1989) forms pairs of sensors to create an array of doublets such that it consists of two identical subarrays where one is a translate of the other. Suppose the signal consists of M planewaves with complex, zero-mean, Gaussian am- plitudes, having the same frequency ~ arriving from M directions (O. tbm), m = 1, . . ., M. We further assume for simplicity that the interference is ab- sent. Suppose we have already detected the signal and our goal is to estimate the M arrival directions ~m' m = 1,,M, which are specified relative to the axis of the doublet (the displacement vector). Denote the data from the two subarrays of sensors by two (~/2~-vectors x and y, assuming ~7 to be even, x = (~1, - ~ -, :~7/2 ), pi = ( T ) 1; exploit ~ drj (~ ~ ' be. ~3) ~ 2 ~ 2~( expticot) dr~+j~t), j = 1, . . ., J/2 . Y=(Y~,.~.,Y7/2), Ye= Then R== =Exx* =AUA*+! R=y = E by* = AURA* ~ where A' U. and 4} are ~7 x M, M x M, and M x M matrices respectively and specified by (A)jm = a;~77 O7 Abe) 7 (U)mm' = 2Pmm, 7 (q>)mm' = exp (i c sin (m) [mm' j = 1, ,J, m = 1, ,M, where ~ is the distance between the two paired sensors. Assuming U to be nonsingular, we observe that the determinant of R=X -·—yR=y = AUNT- ~ *)A* (8.14)

155 vanishes if and only if y=exp~—~ singes), m=l,...,M. This fact can be used to estimate hem, m = I, . . ., M, as follows: regard Rex and Ray as measured covariance matrices. For example, we might subdivide the observation interval T into N equal subintervals ((n—1)T/N,nT/N), n = 1,..., N. where N = ~T/~2;r), replace the integration limits in (~.13) by (n- 1~)T/N and nT/N, n = I,...,N, and denote the integrals by own) and yj~n), j = l, . . ., J; n = 1, . . ., N. Then, put N N RXX = N ~ x~n)x*(n), Ray = N ~ x(n)y*(n) n=1 r~=1 Now substitute these empirical matrices into the left-hand side of (~.14) and find M minima of the absolute value of the determinant as ~ moves on the unit circle centered at the origin of the complex plane. Substitute the M i- v~;u~ Cal 1 ~ ~QlL~;ll~ L~ L11~= 111;1~;111~ ~1~= ~1 v ~ 1V1 vim ~ ~ ~ ~—1 ~ . . . ~ M. Observe that the knowledge of the sensor positions incorporated into A and of the signal powers Pmm' is not required. Thus, this method of estimating the signal arrival directions is free of the costly array calibration. The price to be paid for this is that the two subarrays must be identical, with one being a translate of the other. ~~ ~~ I ~ ~~ +~ +~ ~~— —~^ ~1~^ the— ~/~ ~ — I S.5 Future Research Considerations The assumption that both the signal and the interference plus noise be Gaussian fields is primarily for mathematical convenience since the problem then is completely treatable by linear operators in Hilbert spaces, and Gaus- sian fields are the simplest class of the second-order random fields. How- ever, there are evidences, especially in the case of the ocean acoustics, that the probability distributions of the interference fields considerably deviate from the Gaussian distribution (Middleton, 1987~. Some simple analytical examples, such as the "contaminated Gaussian" distribution (Martin and Schwartz, 1971), have been proposed for the one-dimensional i.i.~. time se- ries. Although the non-Gaussian interference makes the analytical solution to the optimum processing problem infeasible, some suboptimum processing methods are explored in special cases (Monzingo and Miller, 1980~. Since it

156 is unrealistic to completely specify the probability distribution of the inter- ference, a robust method, such as the min-max solution (Huber, 1981), has been sought. The results so far are restricted to the one-dimensional time series having independent identical distributions (Kassam and Poor, 1985), and generalization to the higher dimensional case with dependent distribu- tions should be sought. Another area of investigation is the case where the interference is a nonstationary and inhomogeneous random field, such as a transient disturbance. In this case, one might use a semideterministic cri- terion rather than the totally probabilistic Neyman-Pearson criterion, and estimate (maximum likelihood) the interference z in (~.2) rather than av- erage with respect to its probability distribution. One practical problem in dealing with multidimensional data is computational complexity. Even if there is an explicit algorithm for the optimum signal-processing, the com- plex~ty may be too prohibitive to justify its use. Thus, a trade-off between the detection-estimation performance and the computational complexity, or the cost of processing the data, must be considered. Study of this trade-off is another area of useful research in the future. Bibliography [1] Cox, H., R. M. Zeskind, and M. M. Owen, Robust adaptive beamform- ing, IEEE Trans. Acoust., Speech, Signal Process. 35 (1987), 1365- 1375. [2] FriedIander, B., and B. Porat, Performance analysis of a nub steering al- gorithm based on direction-of-arrival estimation, IEEE Trans. Acoust., Speech, Signal Process. 37 (1989), 461-466. [3] Gabriel, W. F., Adaptive arrays An introduction, [EKE Proc. NatI. Aerosp. Electron. Conf. 64 (1976), 239-272.

157 [4] Haykin, S., ea., Array Processing Applications to Radar, Benchmark papers in electrical engineering and computer science 22, Dowden, Hutchinson and Ross, Stroudsburg, PA, 1980. [5] Helstrom, C. W., Statistical Theory of Signal Detection, 2nd ea., Perg- amon Press, Oxford, Englan(l, 1986, 87-95. [6] Huber, P. J., Robust Statistics, John Wiley and Sons, New York, 1981. [7] Kadota, T. T., and D. M. Romain, Optimum detection of Gaussian signal fields in the multipath-anisotropic noise environment and nu- merical evaluation of detection probabilities, IEEE Trans. Inf. Theory 23 (1977), 167-178. [~] Kassam, S. A., and H. V. Poor, Robust techniques for signal processing: A survey, IEEE Proc. NatI. Aerosp. Electron. Conf. 73 (1985), 433-481. [9] Martin, R. D., and S. C. Schwartz, Robust detection of a known signal in nearly Gaussian noise, IEEE Trans. Inf. Theory 17 (1971), 50-56. [10] Middleton, D., Channel modeling and threshold signal processing in underwater acoustics: An analytical overview, [EKE ]. Oceanic Eng. 12 (1987), 4-28. [11] Monzingo, R. A., and T. W. Miller, Introdluction to Adaptive Arrays, John Wiley and Sons, New York, 1980. [12] Roy, R., and T. Ka~lath, ESPRIT Estimation of Signal Parameters via Rotational Invariance Techniques, IEEE Trans. Acoust., Speech, Signal Process. 37 (1989), 984-995. [13] Steinberg, B. D., Principles of Aperture and Array System Design, John Wiley and Sons, New York, 1976. [14] Wegman, E. J., S. C. Schwartz, and J. B. Thomas, Topics in Non- Gaussian Signal Processing, Springer-VerIag, New York, 1989.