唐 斌 鄭美燕 陳客松 吳宏剛 劉先攀

(1.電子科技大學(xué)電子工程學(xué)院，四川 成都 611731；2.成都航空職業(yè)技術(shù)學(xué)院航空電子工程系，四川 成都 610100；3.中國(guó)民用航空局第二研究所，四川 成都 610041；4.電子科技大學(xué)航空航天學(xué)院，四川 成都 611731)

引 言

近幾十年來(lái)，平面天線在雷達(dá)、通信、衛(wèi)星電視接收等方面得到廣泛的應(yīng)用.通常情況下陣元數(shù)目決定一個(gè)系統(tǒng)的復(fù)雜度和成本，因此使用盡可能少的陣元達(dá)到系統(tǒng)要求是陣列設(shè)計(jì)的重要問(wèn)題.

綜合非均勻平面陣列的陣元激勵(lì)是一個(gè)線性問(wèn)題，然而綜合陣元相位和位置是一個(gè)包含多個(gè)未知量的高度非線性優(yōu)化問(wèn)題.一種方案是對(duì)均勻間隔的陣列進(jìn)行稀疏化設(shè)計(jì)，得到稀疏陣列；另一種方案是對(duì)陣元隨機(jī)稀布，使陣元在稀布過(guò)程中有更大的自由度，可以獲得比稀疏陣列更好的性能，稱為稀布陣列.目前已經(jīng)有許多綜合稀布陣列的方法，包括優(yōu)化算法(如遺傳算法(Genetic Algorithms,GA)[1]、微分進(jìn)化算法(Differential Evolution Algorithm,DEA)[2]、粒子群優(yōu)化方法(Particle Swarm Optimization,PSO)[3]、模擬退火方法[4]、矩陣束方法(Matrix Pencil Method,MPM))和其他綜合技術(shù)等.其中，GA、DEA和PSO適合求解全局最優(yōu)解，但是非常耗時(shí)；MPM已成功應(yīng)用于可分離型分布的平面陣列的稀布綜合中，然而可分離型平面陣列只能保證兩個(gè)主平面內(nèi)方向圖的最佳特性，并不能保證在任一平面內(nèi)方向圖都是最佳的.

使平面陣產(chǎn)生的方向圖在每一剖面上都是最佳方向圖，其關(guān)鍵是采用不可分離型分布的平面陣.如何將MPM擴(kuò)展應(yīng)用到不可分離型分布的陣列綜合中，本文提出一種新的方法-增廣矩陣束方法(Matrix Enhancement and Matrix Pencil, MEMP)[5].

目前還未見(jiàn)將MEMP應(yīng)用到平面陣列綜合中的報(bào)道，本文對(duì)此展開(kāi)研究.首先對(duì)期望的方向圖進(jìn)行采樣得到離散的數(shù)據(jù)集合，并由采樣點(diǎn)數(shù)據(jù)根據(jù)隔斷和堆放的過(guò)程構(gòu)造增廣矩陣，對(duì)此增廣矩陣進(jìn)行奇異值分解，確定在誤差允許范圍內(nèi)所需的陣元數(shù)目；然后基于廣義特征值分解分別求解兩組特征值，并根據(jù)文獻(xiàn)[6]的配對(duì)方法實(shí)現(xiàn)兩組特征值的正確配對(duì)；最后在最小二乘準(zhǔn)則的條件下求得稀布面陣的陣元位置和激勵(lì).仿真試驗(yàn)分別優(yōu)化激勵(lì)為1的平面陣(等激勵(lì)陣列)和切比雪夫平面陣(非等激勵(lì)陣列)，使用盡可能少的陣元逼近均勻分布陣列的方向圖，并保持原陣列的特性，仿真結(jié)果證明了該方法的有效性.

1 增廣矩陣束用于平面陣列綜合

1.1 陣列的最優(yōu)化模型

在三維空間任意排列的陣列如圖1所示.此陣列由N個(gè)陣元組成，位于(ri,θi,φi)第i個(gè)單元的激勵(lì)記為Ri，每個(gè)陣元均為全向輻射元.

圖1 任意陣元的參考坐標(biāo)

根據(jù)電磁波的疊加原理，陣因子可寫(xiě)為[7]

(1)

式中:k=2π/λ，λ為工作波長(zhǎng);

cosαi=sinθsinθicos(φ-φi)+cosθcosθi,

(2)

0≤θ≤π,0≤φ≤2π分別表示方位角和俯仰角.針對(duì)本文的平面陣，式(1)可簡(jiǎn)化為

(3)

式中:dx和dy分別為第i個(gè)陣元的橫坐標(biāo)和縱坐標(biāo);u=sinθcosφ；v=sinθsinφ.

1.2 最小陣元數(shù)目估計(jì)

MEMP是在誤差允許范圍內(nèi)，使用盡可能少的陣元形成新的平面陣列來(lái)逼近期望方向圖.因此，最優(yōu)化問(wèn)題的數(shù)學(xué)描述為

(4)

對(duì)期望的方向圖從u=-1、v=-1到u=1、v=1進(jìn)行均勻采樣，則um=mΔ=m/N，其中m=-N,…,0,…,N；vn=nΔ=n/N，其中n=-N,…,0,…,N.共有(2N+1)(2N+1)個(gè)采樣點(diǎn).

任一采樣點(diǎn)處的值為

(5)

式中yi=ejkdixΔ=ejkdix/N和zi=ejkdiyΔ=ejkdiy/N.

然后，由方向圖的采樣數(shù)據(jù)構(gòu)造增廣矩陣Xe.

(6)

式中：

Xm=

(7)

式中：x(m,n)=f(m-N,n-N)；Xe是Hankel塊矩陣；Xm是Hankel矩陣.參數(shù)K和L的選擇滿足條件：

(8)

對(duì)增廣矩陣Xe進(jìn)行奇異值分解(SVD),表達(dá)式為

(9)

min=min{KL,(2N+2-K)(2N+2-L)}是增廣矩陣Xe較小的維數(shù)；Us、Σs、Vs包含N個(gè)主特征值和主特征向量，Un、Σn、Vn包含剩余特征值和特征向量.具體為：

Us=[u1,u2,…,uN]，

Σs=diag{σ1,σ1,…,σN}，

Vs=[v1,v2,…,vN]，

Un=[uN+1,uN+2,…,umin]，

Σn=diag{σN+1,σN+2,…,σmin}，

Vn=[vN+1,vN+2,…,vmin].

式中，σ1≥σ2≥…≥σmin.增廣矩陣Xe的秩等于非零奇異值的數(shù)目，一般由N個(gè)陣元組成的陣列則有N個(gè)非零奇異值，即σi>0(i=1,…,N)，當(dāng)i>N+1時(shí)，σi=0，因此，Σn為0.式(9)可以化簡(jiǎn)為

(10)

文獻(xiàn)[8]的結(jié)果表明，重要奇異值的數(shù)目要小于總的陣元數(shù)，也就是說(shuō)一些不重要的陣元的貢獻(xiàn)可以由其他重要陣元的組合代替，因此可以舍棄不重要的奇異值來(lái)獲得增廣矩陣Xe的低秩逼近矩陣XQ，這個(gè)低秩逼近矩陣對(duì)應(yīng)著更少的陣元組成的新陣列.通常的處理方法是將這些小的奇異值設(shè)為0，即：

(11)

式中，ΣQ=diag{σ1,σ1,…,σQ,0,…,0},Q≤N.

在實(shí)際的陣列綜合中，Q的最小值可以通過(guò)下式確定[9]

(12)

ε是一個(gè)很小的正數(shù)，ε的選擇取決于重構(gòu)方向圖和期望方向圖的逼近程度.

1.3 求解特征值yi和zi

矩陣束方法[10-11]求解特征值是通過(guò)構(gòu)造兩個(gè)矩陣求其廣義特征值而得到，利用文獻(xiàn)[6]中的配對(duì)算法估計(jì)出(yi,zi)對(duì).

1.3.1 提取特征值yi

低秩矩陣XQ獲得后，求解特征值yi即是求解下式的廣義特征值：

(XQ,f-yXQ,l)v=0，

(13)

式中：XQ，f和XQ，l分別由XQ去掉前L列和后L列得到.等效于求解下列廣義特征值問(wèn)題

(U2-yU1)v=0，

(14)

式中：U2和U1由UQ分別去掉前L行和后L行得到，UQ是式(10)Us的Q個(gè)左特征向量.

因此，矩陣束U2-yU1可以化為

(U2-yU1)=E1(Yd-yI)Ta，

(15)

式中，Yd是由特征值{yi,i=1,…,Q}組成的對(duì)角矩陣.

1.3.2 提取特征值z(mì)i

為提取另一組特征值集合{zi,i=1,…,Q}，引入置換矩陣P

P= [pT(1),pT(1+L), …,pT(1+(K-1)L),

pT(2),pT(2+L), …,pT(2+(K-1)L),

…

…

pT(L),pT(L+L),…,pT(1+(K-1)L)]T.

(16)

矩陣P的元素p(i)是一個(gè)KL×1的向量，且除了第i行為1外，其他皆為0.

用P左乘Us，則得

UsP=PUs.

(17)

由式(14)可知，求解特征值z(mì)i等效于求解下式的廣義特征值

(U2P-zU1P)v=0,

(18)

式中：U2P和U1P由UQP分別去掉前K行和后K行得到，UQP是式(17)UsP的Q個(gè)左特征向量.

因此，矩陣束U2P-zU1P可以化為

(U2P-zU1P)=E1P(Zd-zI)Tb，

(19)

式中,Zd是由特征值{zi,i=1,…,Q}組成的對(duì)角矩陣.

1.3.3 對(duì)特征值yi和zi進(jìn)行配對(duì)

由式(15)和(19)可得：

(20)

通過(guò)標(biāo)量β將矩陣F1和F2線性組合，并對(duì)其對(duì)角化分解的變換矩陣為T(mén).

βF1+(1-β)F2=T-1DT.

(21)

由T、Ta和Tb求得兩組置換矩陣P1和P2：

P1=T-1Ta，P2=T-1Tb.

(22)

P1和P2每一行最大元素的位置構(gòu)成向量p1和p2.p1中第k個(gè)位置所對(duì)應(yīng)的特征值和p2中第k個(gè)位置對(duì)應(yīng)的特征值構(gòu)成正確的特征值對(duì).

文獻(xiàn)[12]應(yīng)用類ESPRIT算法對(duì)特征值配對(duì)，得到的只是特征值的近似值，此配對(duì)方法可得到更精確的特征值.

1.4 求解陣元位置和激勵(lì)

一旦求出特征值yi和zi，陣元位置可以通過(guò)文獻(xiàn)[8]的式(13)求出：

(23)

(24)

(25)

陣元激勵(lì)可通過(guò)下式求得

(26)

矩陣R的對(duì)角線上的元素即是Ri(i=1,…,Q).式中:

(27)

(28)

(29)

式中：

(30)

(31)

(32)

1.5 算法流程

1) 對(duì)期望平面陣的方向圖進(jìn)行采樣，并由采樣點(diǎn)數(shù)據(jù)構(gòu)造增廣矩陣Xe，如式(6)所示.

2) 對(duì)增廣矩陣Xe進(jìn)行SVD,計(jì)算出其奇異值和左特征向量Us.

3) 根據(jù)式(12)確定逼近期望方向圖所需的最小陣元數(shù)目Q.

4) 由式(15)和(19)分別提取特征值yi和zi，并按照式(21)對(duì)特征值yi和zi進(jìn)行配對(duì).

5) 由式(24)和(26)計(jì)算重構(gòu)陣列陣元位置和激勵(lì).

2 仿真實(shí)例

為說(shuō)明增廣矩陣束方法的有效性，本文給出以下兩個(gè)實(shí)例來(lái)減小期望方向圖的陣元數(shù)目，并使重構(gòu)陣列保持原陣列的特性.

例1： 綜合激勵(lì)為1的矩形平面陣

設(shè)有一7×7的矩形平面陣，其陣元均勻分布在矩形柵格上，陣元間距dx=dy=λ/2，方向圖如圖2(a)所示.首先對(duì)此方向圖進(jìn)行采樣，共有(2×49+1)(2×49+1)個(gè)采樣點(diǎn)，并由這些采樣點(diǎn)數(shù)據(jù)按照隔斷和堆放的過(guò)程構(gòu)造增廣矩陣.矩陣束參數(shù)K=L=50.由式(12)可知，當(dāng)ε=10-6所需的陣元數(shù)為Q=36，然后基于廣義特征值分解提取兩組特征值，并應(yīng)用配對(duì)算法對(duì)特征值進(jìn)行配對(duì)，最后在最小二乘準(zhǔn)則的約束下求得稀布面陣的陣元位置和激勵(lì).根據(jù)上述流程求得重構(gòu)陣列的方向圖如圖2(b)所示.圖3對(duì)比了期望方向圖和重構(gòu)方向圖的切面圖.

(a) 均勻平面陣的方向圖

(b) 重構(gòu)平面陣的方向圖圖2 綜合激勵(lì)為1的矩形平面陣的方向圖

由圖可知，非均勻分布的陣列只需要36個(gè)陣元便可精確重構(gòu)均勻分布時(shí)需要49個(gè)陣元才能產(chǎn)生的方向圖，此例可節(jié)省27%的陣元.圖4給出了均勻分布和稀布后的陣元位置圖.表1列出了均勻分布的陣元位置以及稀布后的陣元位置和激勵(lì).因陣元是對(duì)稱分布的，只給出了第一象限內(nèi)陣元的位置和激勵(lì).

表1 均勻陣元的位置和非均勻陣列陣元的位置與激勵(lì)

圖3 方向圖的截面圖

圖4 陣元位置圖

例2： 綜合切比雪夫平面陣

(a) 切比雪夫平面方向圖

(b) 重構(gòu)平面陣方向圖圖5 綜合切比雪夫平面陣的方向圖

設(shè)有一4×4的切比雪夫平面陣，其陣元均勻分布在矩形柵格上，陣元間距dx=dy=λ/2，要求其環(huán)狀副瓣的電平為-20 dB.在此例中，共有(2×16+1)(2×16+1)個(gè)采樣點(diǎn)，增廣矩陣束參數(shù)K=L=17.按照例1的步驟求得稀布平面陣的陣元位置和激勵(lì).圖5是切比雪夫方向圖和稀布平面陣方向圖.旁瓣電平為-16.504 2 dB.圖6對(duì)比了重構(gòu)方向圖和期望方向圖的切面圖，進(jìn)一步增加陣元數(shù)目也不會(huì)改善方向圖的特性.雖然旁瓣電平有小幅抬高(3.495 8 dB)，但可以節(jié)省43.75%的陣元.圖7給出了均勻切比雪夫平面陣的陣元位置和稀布后的陣元位置圖.表2列出了均勻分布條件下陣元的位置和激勵(lì)以及稀布后的陣元位置和激勵(lì).

圖6 方向圖的截面圖

圖7 陣元位置圖

切比雪夫平面陣的陣元位置(激勵(lì))0．25，0．75(0．6854)0．75，0．75(0．2285)0．25，0．25(0．9008)0．75，0．25(0．6854)非均勻陣列的陣元位置(激勵(lì))-0．5819，0．5815(0．8882)0．0049，0．6404(1．3491)0．5845，0．5871(0．8724)-0．6362，-0．0036(1．3614)0，0(1．7599)0．6362，0．0036(1．3614)-0．5845，-0．5871(0．8724)-0．0049，-0．6404(1．3491)0．5819，-0．5815(0．8882)

3 結(jié) 論

提出了一種基于增廣矩陣束(MEMP)方法的減小最小陣元數(shù)目、求解陣元位置和設(shè)計(jì)激勵(lì)幅度的平面陣列綜合方法，與基于隨機(jī)優(yōu)化的算法相比，基于增廣矩陣束的方法是一種非迭代算法，適合于要求窄波束、低副瓣陣列的設(shè)計(jì).另外，與基于矩陣束方法的可分離型分布的平面陣列綜合相比，增廣矩陣束方法可用于不可分離型分布的平面陣的綜合，從而保證方向圖在每一剖面的最佳特性.本文的探討豐富了增廣矩陣束方法在稀布平面陣列綜合中的應(yīng)用，為其他種類的稀布陣列綜合提供了有益的提示，也為其在工程應(yīng)用中提供了有價(jià)值的參考.

[1] 陳客松,何子述.平面稀布天線陣列的優(yōu)化算法[J].電波科學(xué)學(xué)報(bào),2009,24(2):193-198.

CHEN Kesong,HE Zishu. Synthesis approach for sparse plane arrays[J]. Chinese Journal of Radio Science,2009,24(2):193-198.(in Chinese)

[2] KURUP D G, HIMDI M. Synthesis of uniform amplitude unequally spaced antenna arrays using the differential evolution algorithm[J]. IEEE Trans on Antenna and Propagation,2003,51(9): 2210-2217.

[3] DELIGKARIS K.V,ZAHARIS Z.D.Thinned planar array design using boolean PSO with velocity mutation[J].IEEE Trans on Magnetics,2009, 45(3):1490-1493.

[4] 郭陳江,張 鋒,丁 君,等.基于循環(huán)差集與模擬退火法的陣列綜合[J].電波科學(xué)學(xué)報(bào),2007,22(6):962-964.

GUO Chenjiang,ZHANG Feng,DING Jun,et al. Array synthesis using cyclic difference sets and simulated annealing[J].Chinese Journal of Radio Science,2007,22(6):962-964.(in Chinese)

[5] HUA Yingbo.Estimating Two-dimensional frequencies by matrix enhancement and matrix pencil[J].IEEE Trans on Sgnal Processing,1992,40(9):2267-2280.

[6] ROUQUETTE S, NAJIM M.Estimation of frequencies and damping factors by two dimensional ESPRIT type methods[J].IEEE Trans on Signal Processing,2001,49(1):237-245.

[7] CHEN D K.Optimization techniques for antenna arrays[J].Processing of the IEEE,1971,59(12):1664-1674.

[8] MILLER E K,GOODMAN D M.A pole-zero modeling approach to linear array synthesis I:the unconstrained solution[J].Radio Science,1983,18(1):57-69.

[9] LIU Yanhui, NIE Zaiping, LIU Qinghuo. Reducing the number of elements in a linear antenna array by the matrix pencil method[J].IEEE Trans on Antennas and Propagation,2008,56(9):2955-2962.

[10] SAPKAR T K,PEREIRA O.Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[J].IEEE Antennas and Propagation Magazine,1995,37(1):48-54.

[11] HUA Yingbo, SAPKAR T K. Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[J].IEEE Trans on Acoustics Speech and Signal Processing,1990,38(5):814-824.

[12] 周云鐘,陳天麒.多信號(hào)極化與到達(dá)角估計(jì)算法[J].電波科學(xué)學(xué)報(bào),1997,12(2):220-224.

ZHOU Yunzhong,CHEN Tianqi.Angles of arrival polarizations estimation[J].Chinese Journal of Radio Science,1997,12(2):220-224.(in Chinese)