運用增廣矩陣束方法稀布優(yōu)化平面陣

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

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

引 言

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

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

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

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

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

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

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

圖1 任意陣元的參考坐標

根據電磁波的疊加原理，陣因子可寫為[7]

(1)

式中:k=2π/λ，λ為工作波長;

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

(2)

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

(3)

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

1.2 最小陣元數目估計

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

(4)

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

任一采樣點處的值為

(5)

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

然后，由方向圖的采樣數據構造增廣矩陣Xe.

(6)

式中：

Xm=

(7)

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

(8)

對增廣矩陣Xe進行奇異值分解(SVD),表達式為

(9)

min=min{KL,(2N+2-K)(2N+2-L)}是增廣矩陣Xe較小的維數；Us、Σs、Vs包含N個主特征值和主特征向量，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的秩等于非零奇異值的數目，一般由N個陣元組成的陣列則有N個非零奇異值，即σi>0(i=1,…,N)，當i>N+1時，σi=0，因此，Σn為0.式(9)可以化簡為

(10)

文獻[8]的結果表明，重要奇異值的數目要小于總的陣元數，也就是說一些不重要的陣元的貢獻可以由其他重要陣元的組合代替，因此可以舍棄不重要的奇異值來獲得增廣矩陣Xe的低秩逼近矩陣XQ，這個低秩逼近矩陣對應著更少的陣元組成的新陣列.通常的處理方法是將這些小的奇異值設為0，即：

(11)

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

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

(12)

ε是一個很小的正數，ε的選擇取決于重構方向圖和期望方向圖的逼近程度.

1.3 求解特征值yi和zi

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

1.3.1 提取特征值yi

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

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

(13)

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

(U2-yU1)v=0，

(14)

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

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

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

(15)

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

1.3.2 提取特征值zi

為提取另一組特征值集合{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)是一個KL×1的向量，且除了第i行為1外，其他皆為0.

用P左乘Us，則得

UsP=PUs.

(17)

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

(U2P-zU1P)v=0,

(18)

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

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

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

(19)

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

1.3.3 對特征值yi和zi進行配對

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

(20)

通過標量β將矩陣F1和F2線性組合，并對其對角化分解的變換矩陣為T.

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

(21)

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

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

(22)

P1和P2每一行最大元素的位置構成向量p1和p2.p1中第k個位置所對應的特征值和p2中第k個位置對應的特征值構成正確的特征值對.

文獻[12]應用類ESPRIT算法對特征值配對，得到的只是特征值的近似值，此配對方法可得到更精確的特征值.

1.4 求解陣元位置和激勵

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

(23)

(24)

(25)

陣元激勵可通過下式求得

(26)

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

(27)

(28)

(29)

式中：

(30)

(31)

(32)

1.5 算法流程

1) 對期望平面陣的方向圖進行采樣，并由采樣點數據構造增廣矩陣Xe，如式(6)所示.

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

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

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

5) 由式(24)和(26)計算重構陣列陣元位置和激勵.

2 仿真實例

為說明增廣矩陣束方法的有效性，本文給出以下兩個實例來減小期望方向圖的陣元數目，并使重構陣列保持原陣列的特性.

例1： 綜合激勵為1的矩形平面陣

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

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

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

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

表1 均勻陣元的位置和非均勻陣列陣元的位置與激勵

圖3 方向圖的截面圖

圖4 陣元位置圖

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

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

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

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

圖6 方向圖的截面圖

圖7 陣元位置圖

切比雪夫平面陣的陣元位置(激勵)0．25，0．75(0．6854)0．75，0．75(0．2285)0．25，0．25(0．9008)0．75，0．25(0．6854)非均勻陣列的陣元位置(激勵)-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 結 論

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

[1] 陳客松,何子述.平面稀布天線陣列的優(yōu)化算法[J].電波科學學報,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].電波科學學報,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] 周云鐘,陳天麒.多信號極化與到達角估計算法[J].電波科學學報,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)