JOURNAL OF

AERONAUTICS

Chinese Journal of Aeronautics, (2014),27(6): 1538-1543

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

cja@buaa.edu.cn www.sciencedirect.com

DOA estimation and mutual coupling calibration with the SAGE algorithm

Xiong Kunlai, Liu Zhangmeng *, Liu Zheng, Jiang Wenli

College of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China

Received 31 December 2013; revised 29 April 2014; accepted 19 June 2014 Available online 19 October 2014

Abstract In this paper, a novel algorithm is presented for direction of arrival (DOA) estimation and array self-calibration in the presence of unknown mutual coupling. In order to highlight the relationship between the array output and mutual coupling coefficients, we present a novel model of the array output with the unknown mutual coupling coefficients. Based on this model, we use the space alternating generalized expectation-maximization (SAGE) algorithm to jointly estimate the DOA parameters and the mutual coupling coefficients. Unlike many existing counterparts, our method requires neither calibration sources nor initial calibration information. At the same time, our proposed method inherits the characteristics of good convergence and high estimation precision of the SAGE algorithm. By numerical experiments we demonstrate that our proposed method outperforms the existing method for DOA estimation and mutual coupling calibration.

© 2014 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA.

Open access under CC BY-NC-ND license.

KEYWORDS

Array self-calibration; Convergence;

Direction of arrival estimation;

Mutual coupling;

Space alternating generalized

expectation-maximization

algorithm

1. Introduction

The problem of estimating the direction of arrival (DOA) of multiple narrowband signals plays an important role in many fields, including radar, wireless communications, seismology, sonar and so on. Most of the existing DOA estimation methods rely heavily on perfect knowledge of the array manifold. However, the array manifold is often affected by array imperfections in practice, such as mutual coupling, gain/phase uncertainty and sensor location error. The performance of

* Corresponding author. Tel.: +86 731 84573489. E-mail address: liuzhangmeng@nudt.edu.cn (Z. Liu). Peer review under responsibility of Editorial Committee of CJA.

the DOA estimation methods may deteriorate significantly without array manifold calibration.1-6

In this paper, we focus on the mutual coupling calibration and DOA estimation problem for deterministic signals. Various array calibration methods have been proposed with respect to the mutual coupling effect.7-14 In Ref. 7, an iterative least-mean-square algorithm is proposed to estimate the calibration matrix, but it requires preliminary calibration. In Ref. 8, Fabrizio and Alberto present an online mutual coupling compensation algorithm for uniform and linear arrays, which could estimate the coupling coefficients through an alternating minimization technique, but the convergence is not well guaranteed. The iterative algorithm in Ref. 9 compensates the mutual coupling effect and the gain/phase perturbation simultaneously. However, the convergence rate is low, and the computational cost is very expensive. In Ref. 10, a maximum-likelihood algorithm is presented for array calibration, which takes into consideration the array imperfections of mutual coupling,

http://dx.doi.org/10.1016/j.cja.2014.10.020

1000-9361 © 2014 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. Open access under CC BY-NC-ND license.

gain/phase perturbation, as well as sensor position errors. However, the method requires a set of calibration sources with known locations, which is hardly available in practice. Compared with the methods that require calibration sources or initial calibration information, self-calibration methods are more attractive under some circumstances. Liu et al.11 present a novel DOA estimation algorithm with uniform linear arrays in the presence of mutual coupling via blind calibration. Weiss and Friedlander12 present a self-calibration algorithm that compensates for mutual coupling and gain/phase perturbation. However, they require multi-dimensional parameter searching, and the convergence is not well guaranteed. In Ref. 13, an array self-calibration method to compensate for sensor position perturbation is presented by the SAGE algorithm. It remains robust in critical scenarios including large sensor position errors and closely located signals. But it only focuses on sensor position perturbations. In Ref. 14, a unified framework is proposed for DOA estimation and array self-calibration to all kinds of array perturbations, which can achieve good performance in scenarios with low SNR, limited snapshots and spatially adjacent signals, but it only aims at stochastic signals.

On the other hand, among the existing DOA estimation methods, the maximum-likelihood (ML) ones achieve the best asymptotic performance. The main drawback of the ML approach is the high computational complexity caused by maximization of the likelihood function. To overcome this drawback, the expectation maximization (EM) algorithm has been derived for both deterministic and stochastic signal mod-els.15-18 The SAGE algorithm is a variation of the widely used EM algorithm, which updates subsets of parameters sequentially in one iteration. Through the augmentation scheme specified by the EM or SAGE algorithm, the complicated multidimensional searching involved in maximizing the likelihood functions can be simplified to one-dimensional searching. It has been proved in Ref. 19 that SAGE converges faster than EM while retaining the advantages in numerical simplicity and stability.

This paper develops a novel calibration method to compensate for unknown mutual coupling in uniform linear arrays. The method jointly estimates the DOA parameters and the coupling coefficients by using the SAGE algorithm. First, an array output model with unknown mutual coupling is presented to highlight the relationship between the array output and the mutual coupling coefficients. Then we derive a SAGE-based calibration algorithm that updates the DOA parameters and the coupling coefficients sequentially based on this model. Compared to the existing methods, our algorithm requires neither calibration sources nor initial calibration information. At the same time, simulation results demonstrate that it can achieve higher parameter estimation precision.

2. Perturbed array output formulation

Suppose that K unknown and deterministic narrowband signals impinge onto an M element uniform linear array (ULA) from directions of # = [#1, #2,..., #K]T, N snapshots are collected, and the output of the accurately calibrated array at time t is

where A(#) = [«(#1), «(#2),..., «(#K)] and «(tf^ — jM, ejuk2, ■■■, ejutM]T are the array responding matrix and vector, respectively, uk m — (2pdm sin #k)/k, where k represents the signal wavelength and d1, d2,..., dM are the distances of the M array sensors from the reference located on the array axes. sk(t) is the waveform of the kth signal at time t, s(t) = [s1(t), S2(t),..., SK(t)]T, v(t) = [v1(t), V2(t),..., VM(t)]T is the noise vector with variance r2/, xpf(t) — [xpf(t) , xpf(t) , ■■■ , xM (t)] and the superscript "pf" is the short of ''perturbation-free''.

When mutual coupling effects are taken into consideration, the array output is given as follows:

x(t) = HM + v(t) = A'(#)s(t) + v(t)

where A'(#) = [«'(#1), «'(#2),...,«'(#K)], and a'(#) is the perturbed array responding vector, «' (#) = C«(#) with C denoting the mutual coupling matrix. The mutual coupling matrix can be written more explicitly as

toeplitz^j1, c1 , c2 , ■■■ ,cP , 0TM_P_1)x1J ^ in ULA,10 where

cp is the coupling coefficient of the two sensors displaced by (p _ 1) times the inter-element spacing of the ULA, and cp is very small when p > P and cp is thus neglected.

Although the directional information of the incident signals is still reserved in the perturbed array outputs, it can hardly be extracted directly based on the array output formulation in Eq. (2). In order to highlight the perturbation-free signal components, and also make clearer the relationship between the array output and the perturbation parameters, the following model can be established for the perturbed array output.

x(t) = A (#)s(t) + v(t) = A(#)s(t) + Q(t)c + v(t)

where c = [c1, c2, ■■■, cP] stands for the coupling coefficient vector. And the following equation holds for Q(t) and c.

Q(t)c = [A'(#)- A(#)]s(t) = (C - lM)A(#)s(t)

xpf(t) = J2a(#k)sk(t) + v(t) = A(#)s(t) + v(t)

where /M stands for M dimension unitary matrix. The matrix of Q(t) can be expressed explicitly according to [Q(t)]:p = GpA(#)s(t) and Gp = 0C/0cp contains nonzero elements of 1 only on the ±p diagonals. Such an expression of Q(t) is concluded by taking the differentiations of both sides of Eq. (4) with respect to cp.

The perturbed array output formulation given in Eq. (3) reveals the relationship between the array output and the coupling coefficients, and such a relationship can be exploited to estimate the coupling coefficients by maximizing the log-likelihood of the array output x(t). Meanwhile, for a particular coupling coefficients estimate, the log-likelihood can be maximized to determine the signal directions. In the next section, we follow this guideline to calibrate the array imperfection and estimate the signal directions iteratively.

3. SAGE-based algorithm for DOA estimation and array calibration

Both the EM and SAGE algorithms have been applied to the problem of DOA estimation without array imperfections, and they show good convergence properties. The EM algorithm is a general numerical method for finding the maximum-likelihood estimates which is characterized with simple implementation and stable convergence. The SAGE algorithm is a

variation of the EM algorithm, which follows a more flexible optimization scheme and converges faster than the EM algorithm. Instead of estimating all the parameters at once in EM, SAGE breaks up the problem into several smaller ones and uses EM to update the parameter subsets associated with each decomposed problem. In order to solve the DOA estimation problem in the presence of mutual coupling, this section extends the existing research to calibrate the mutual coupling effect together with estimating the source directions by the SAGE algorithm.

Similar to the DOA estimation problem in calibrated arrays, we construct an augmented data vector as follows to separate different signal components contained in the array output based on the model of the perturbed array output in Eq. (3).

yk{t) = d{#k)sk{t) + v{t)

= a(#k)sk(t) + Qk (t)c + v(t) k = 1,2,..., K (5)

where Qk(t)c = [«'(#*) - «(#k)]sk(t), [Qk(t)]:,p = [«'(#k)-

«(#k)]sk(t) = Gk,pa(#k)sk(t), Gkp is a variation of Gp corresponding to the kth signal component, and the noise vector v(t) ~ N(0, r2lu).

In the framework of the SAGE algorithm, the likelihood of yi(t), J2(t),...,yKyt) is given by

p(jk(t)|Xk(t)) = |pr2lM|-1ex^ - - ||jk(t)-«'(#k)sk(t)||2

where Xk(t) is the unknown parameter set associated with yk(t), and it is a subset of X = {#, {s(tn)}N=1, c, r2}. Obviously, y1(t), y2(t),...,yK(t) are normally distributed with the same variances and different means. The log-likelihood of {yi(tn), yi(tn),..., jK(tn)}1i can be obtained by calculating

the mean of the log-likelihoods of {yk(tn)}№= 1 for different k.

L({Уl(tn), У2(tn), ... , jK(tn)}N=1) = £({yk(tn)}'=1)

= MN ln r2 + tt-^EE|LVk(tn)-a(#k)sk(tn)||2 (7)

According to Eq. (7), the log-likelihood of {Ji (tn), J2(tn),...,yK(tn)}'=1 can be written as

L({Уl(tn), y2(tn) ..., jK(tn)}'=l)

= KMNln r2 + J2Hyk(tn)-a'(#k)sk(tn)||

n=1 k=1 NK

= KMN ln r2 + ^||jk(tn)- a(#k)sk(tn)- Qk(tn)c||2

r2 n=1 k=1

The second expression of L({ji(tn), J2(tn),..., y^n)}"^) in Eq. (8) will be used for optimizing c.

The SAGE algorithm updates the unknown parameters with an iterative guideline. Each iteration of the SAGE algorithm consists of an E-step and an M-step. Given the estimate of the (q - 1)th iteration Xq-1, the qth iteration consists of the following steps.

Step 1. E-step Calculate

F(Xq; Xq-1) = E[L({yi (tn), J2(tn), ■■■, Jz(t„)}N= i) I {*&)£ i, Xq-

where the superscript q stands for the iteration index, and Eq. (9) is equivalent to computing the following conditional expectations,

Jk(t) = E[jk(t)|x(t), Xq-1]

= a(#k-1)sk-1(t) + x(t) - A'(cq-1, #q-1)sq-1(t) (10)

RJk = E

nY>k(tn)JH (tn)|{x(tn)}N=i, Xq-n=1

= k(tnW (tn ))H

In order to emphasize the dependence on the array perturbation parameters, we use A'(cq-1, #q-1) to denote the perturbed responding matrix in Eq. (10), what relies on the current estimates of the mutual coupling coefficients.

Step 2. M-step

The unknown parameters contained in X can be updated by maximizing F(Xq; Xq-1), and the updating schemes, which are derived by setting the differentiations of F(Xq; Xq-1) to those parameters to zero, are given as follows:

#k = argmax

(«' (cq-1,#))HRj a (cq-1,#)

||a ' (cq-

Sk(tn) =

(r2)q =

(a (cq-1, #k))Hyk(tn)

|a ' (cq-1,#kk)||

cq = {E

Y^^^n)-a (#k)sk(tn

(tn )Qk(tn

n=1 k=1 NK

(tn)(yk(tn)- a(#k)sk(tn))

As can be seen from Eq. (12) to Eq. (15), the update processes of #k and {sk(tn)}"=1 mainly depend on the corresponding augmented data {yk(tn)}"= and they do not matter

with {yi(tn),y2 (tn), ...,Уk-1(tn),Уk+1(tn), . .. , yK(tn)}N=l. As the parameters of c and r2 exist in the models of all the signal components, their update processes rely on {yT (tn),y2(tn),..., yK(tn)}"=1. The iteration process can be divided into two parts, the first one contains the sequential updating of #k and {sk(tn)}1 t for different k, and the second one updates c and r2 by using all the augmented data. Straightforward formulation of Qk(t) in the iterations is very complicated, thus we carry out some deeper analysis in the Appendix A on the update processes associated with this matrix to simplify the implementation of the iteration.

4. Simulation results

The simulations in this section will be carried out to demonstrate the performance of the proposed algorithms. The existing methods of the S-S method8 and the Y-L method,20 as well as the Cramer-Rao lower bound (CRLB),21 are also implemented for performance comparison.

We consider a uniform linear array of 5 sensors with interelement spacing of half a wavelength of the incident signals. Suppose that two equal-power independent signals impinge onto the array from directions of _90 and 25°. The coupling coefficient vector is c = [0.4 + 0.3j, 0.3 + 0.1j]T. The required initial DOA estimate #0 can be initialized by subspace methods such as MUSIC or the standard ML approach. The initial coupling coefficient vector estimate is chosen to be c0 = [0, 0]T. And s0(t) is obtained by computing the least squares estimation for # = #0 and c = c0. The algorithm is terminated if the relative increase in the data likelihood function is smaller than 10_4 or the number of iterations achieves 50.

The average root-mean-square error (RMSE) of the DOA estimates of the two signals is considered for statistical direction estimation precision evaluation, which is defined as

RMSE# =

£||#w - #w||/(KW)

where W denotes the number of Monte Carlo experiments. #w and are the estimated and true directions in the wth simulation.

Fig. 1(a) shows the RMSE of the DOA estimates of different methods against SNR computed via 500 Monte Carlo runs for each SNR (the number of snapshots is fix to 50). Fig. 1(b) shows the RMSE of the DOA estimates of different methods versus the number of snapshots (with SNR = 5 dB) computed via 500 Monte Carlo runs for each number of snapshots. We can observe that the proposed SAGE-based method can achieve higher estimation accuracy than S-S method and Y-L Method, and it improves with a similar speed as that of the CRLB.

The RMSE of the mutual coupling coefficients is used for array calibration precision evaluation, which is defined as

RMSEc =

Eiicw - c|i2/w

where c keeps constant in each scenario, and cw is the estimated coupling coefficient vector in the wth simulation.

Fig. 2(a) shows the RMSE of the mutual coupling coefficients versus SNR in the same simulations (the number of snapshots is fix to 50). Fig. 2(b) shows the RMSE of the mutual coupling coefficients estimates of different methods versus the number of snapshots (with SNR = 5 dB) computed

Fig. 1 RMSE of DOA estimates of different methods.

Fig. 2 RMSE of coupling coefficient estimates of different methods.

via 500 Monte Carlo runs for each number of snapshots. The results illustrate that SAGE-based Method achieves the best performance among the three methods, and its RMSE is very close to the CRLB in the considered scenarios.

5. Conclusion

In this letter, we propose a SAGE-based algorithm for DOA estimation and mutual coupling calibration in the case of uniform linear arrays and unknown determined signals. The implementation of the proposed method follows the standard SAGE iterations and thus has guaranteed convergence. The simulation results show that compared with the existing methods, the proposed one performs much better in both DOA estimation and array calibration precisions, especially when the SNR is low. The guideline that the new method follows can also be extended to the calibration of other kinds of array imperfections, such as gain/phase perturbation and sensor location error, but the straightforward extensions are excluded in this letter due to space limitation.

Acknowledgment

The work was supported by the National Natural Science Foundation of China (No. 61302141).

(hk)p = «H(#k)GHp-1 E

l>k(tn)s?(t)-«(#k)sk (tn )sH(t))

= «H(#k)G«£(R9k ,sk - a(#k)Rk) = aH(#k)GHp

«(#)(«' (c9-1,

a (cq-1, #k)ll

R?k a (c9-1, #k)

a (cq-1, #k)ll2

k ' it p

£sk(t)«H(#k)GHpi Gkp2 a(#k)sk(t)

£>(0^(0

flH(#k)GHpt Gk p «(#k)

(a (cq-1, #k))X a (cq-1, #k)

#k)) Ryt

lla (c(q-1), #k)ll2

aH (#k)GHp, Gk p2 a(#k)

where (hk)p represents the pth element of hk and (Rk)pip stands for the (p1, p2) th element of Rk.

Finally, Eqs. (10), (11) and Eqs. (A3)-(A5) can be substituted into Eqs. (12)-(15) to yield more convenient steps for updating the unknown parameters and implementing the iterations.

Appendix A

According to Eqs. (10), (11) and (13), one can obtain that

R9 = E

sk(tnKH(tn)|{*(tn)}Nl1, Xq1

(a'(c9-1, #^-1))HRyt "'(cq-1, #k)

ll«'(cq-1, tf9)!

Rq — e

Ryk ,sk " E

k(tn)skH(tn)|{x(tn)}l1, X9-1

R\k a'(c9-1, #k)

ll«'(cq-1, #k)ll2

Moreover, according to Eqs. (12), (13), (A1) and (A2), one can get that

Elljk(tn)-«' (#k )Sk(tn ,n=1

tr[R9k - a(cq-1;#k)(Ryk ,sk)H - ^ ,sk(a'(cq-1;#9))H +RJk a (c9-1, #k)(«' (c9-1, #k))H]

a (c9-1, #k)(«' (c9-1, tf9))1"

la' (c9-1, #k)ll2

Denote hk — N ^^N=1Qf (t»)(Vk(t»)_ «(#kH(t„))] and

Rk — N ! QH(tn)Qk(tn^ , then the elements of the vector

and the matrix can be computed according to [Qk(t)]:,p = Gk,p«(#k)sk(t) as follows:

References

1. Roller C, Wasylkiwskyj W. Effects of mutual coupling on superresolution DF in linear arrays. IEEE international conference on acoustics, speech, and signal processing; 1992 Mar 23-26; San Francisco, CA, USA: IEEE; 1992. p. 257-260.

2. Liu ZM, Huang ZT, Zhou YY. Bias analysis of MUSIC in the presence of mutual coupling. IET Signal Proc 2009;3(1):74-84.

3. Weiss AJ, Friedlander B. Mutual coupling effects on phase-only direction finding. IEEE Trans Antennas Propag 1992;40(5):535-41.

4. Yang B, He F, Jin J, Xiong HG, Xu GH. DOA estimation for attitude determination on communication satellites. Chin J Aeronaut 2014;27(3):670-7.

5. Liu A, Liao G, Zeng C, Yang Z, Xu Q. An eigenstructure method for estimating DOA and sensor gain-phase errors. IEEE Trans Signal Process 2011;59(12):5944-56.

6. Chen YM, Lee JH, Yeh CC, Mar J. Bearing estimation without calibration for randomly perturbed arrays. IEEE Trans Signal Process 1991;39(1):194-207.

7. Hung EKL. Matrix-construction calibration method for antenna arrays. IEEE Trans Aerosp Electron Syst 2012;36(3):819-28.

8. Fabrizio S, Alberto S. A novel online mutual coupling compensation algorithm for uniform and linear arrays. IEEE Trans Signal Process 2007;55(2):560-73.

9. Friedlander B, Weiss AJ. Direction finding in the presence of mutual coupling. IEEE Trans Antennas Propag 1991;39(3):273-84.

10. Ng BC, See CMS. Sensor-array calibration using a maximum-likelihood approach. IEEE Trans Antennas Propag 1996;44(6):827-35.

11. Liu Z, Huang Z, Wang FH, Zhou YY. DOA estimation with uniform linear arrays in the presence of mutual coupling via blind calibration. Signal Process 2009;89(7):1446-56.

12. Weiss AJ, Friedlander B. DOA and steering vector estimation using a partially calibrated array. IEEE Trans Aerosp Electron Syst 1996;32(3):1047-57.

13. Chung PJ, Wan S. Array self-calibration using SAGE algorithm. The 5th IEEE international conference on sensor array and

multichannel signal processing workshop; 2008 July 21-23; Darmstadt Germany: IEEE; 2008. p.165-9.

14. Liu ZM, Zhou YY. A unified framework and sparse Bayesian perspective for direction-of-arrival estimation in the presence of array imperfections. IEEE Trans Signal Process 2013;61(15): 3786-98.

15. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J R Stat B 1977;39(1):1-38.

16. Feder M, Weinstein E. Parameter estimation of superimposed signals using the EM algorithm. IEEE Trans Acoust Speech Signal Process 1988;36(4):477-89.

17. Miller MI, Fuhrmann DR. Maximum-likelihood narrow-band direction finding and the EM algorithm. IEEE Trans Acoust Speech Signal Process 1990;38(9):1560-77.

18. Mada KK, Wu HC, Iyengar SS. Efficient and robust EM algorithm for multiple wideband source localization. IEEE Trans Veh Technol 2009;58(6):3071-5.

19. Pei JC, Johann FB. Comparative convergence analysis of EM and SAGE algorithms in DOA estimation. IEEE Trans Signal Process 2001;49(12):2940-8.

20. Ye ZF, Liu C. On the resiliency of MUSIC direction finding against antenna sensor coupling. IEEE Trans Antennas Propag 2008;56(2):371-80.

21. Liu ZM. Conditional cramer-rao lower bounds for DOA estimation and array calibration. IEEE Signal Process Lett 2014;21(3): 361-4.

Xiong Kunlai is a Ph.D. student in the College of Electronic Science and Engineering at National University of Defense Technology . He received the B.S. and M.S. degrees in Electronic Engineering Institute in 2008 and 2011 respectively. His main research interests are statistical signal processing, blind separation, array processing.

Liu Zhangmeng is a lecturer in the College of Electronic Science and Engineering at National University of Defense Technology. His main research interests are statistical signal processing, array processing.

Liu Zheng is an associate professor in the College of Electronic Science and Engineering at National University of Defense Technology. He received his P.h.D. degree in 2005. His main research interests are integrated electronic warfare system and technology.