运用增广矩阵束方法稀布优化平面阵

唐 斌 郑美燕 陈客松 吴宏刚 刘先攀

(1.电子科技大学电子工程学院，四川 成都 611731；2.成都航空职业技术学院航空电子工程系，四川 成都 610100；3.中国民用航空局第二研究所，四川 成都 610041；4.电子科技大学航空航天学院，四川 成都 611731)

引 言

近几十年来，平面天线在雷达、通信、卫星电视接收等方面得到广泛的应用.通常情况下阵元数目决定一个系统的复杂度和成本，因此使用尽可能少的阵元达到系统要求是阵列设计的重要问题.

综合非均匀平面阵列的阵元激励是一个线性问题，然而综合阵元相位和位置是一个包含多个未知量的高度非线性优化问题.一种方案是对均匀间隔的阵列进行稀疏化设计，得到稀疏阵列；另一种方案是对阵元随机稀布，使阵元在稀布过程中有更大的自由度，可以获得比稀疏阵列更好的性能，称为稀布阵列.目前已经有许多综合稀布阵列的方法，包括优化算法(如遗传算法(Genetic Algorithms,GA)[1]、微分进化算法(Differential Evolution Algorithm,DEA)[2]、粒子群优化方法(Particle Swarm Optimization,PSO)[3]、模拟退火方法[4]、矩阵束方法(Matrix Pencil Method,MPM))和其他综合技术等.其中，GA、DEA和PSO适合求解全局最优解，但是非常耗时；MPM已成功应用于可分离型分布的平面阵列的稀布综合中，然而可分离型平面阵列只能保证两个主平面内方向图的最佳特性，并不能保证在任一平面内方向图都是最佳的.

使平面阵产生的方向图在每一剖面上都是最佳方向图，其关键是采用不可分离型分布的平面阵.如何将MPM扩展应用到不可分离型分布的阵列综合中，本文提出一种新的方法-增广矩阵束方法(Matrix Enhancement and Matrix Pencil, MEMP)[5].

目前还未见将MEMP应用到平面阵列综合中的报道，本文对此展开研究.首先对期望的方向图进行采样得到离散的数据集合，并由采样点数据根据隔断和堆放的过程构造增广矩阵，对此增广矩阵进行奇异值分解，确定在误差允许范围内所需的阵元数目；然后基于广义特征值分解分别求解两组特征值，并根据文献[6]的配对方法实现两组特征值的正确配对；最后在最小二乘准则的条件下求得稀布面阵的阵元位置和激励.仿真试验分别优化激励为1的平面阵(等激励阵列)和切比雪夫平面阵(非等激励阵列)，使用尽可能少的阵元逼近均匀分布阵列的方向图，并保持原阵列的特性，仿真结果证明了该方法的有效性.

1 增广矩阵束用于平面阵列综合

1.1 阵列的最优化模型

在三维空间任意排列的阵列如图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是在误差允许范围内，使用尽可能少的阵元形成新的平面阵列来逼近期望方向图.因此，最优化问题的数学描述为

(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个阵元才能产生的方向图，此例可节省27%的阵元.图4给出了均匀分布和稀布后的阵元位置图.表1列出了均匀分布的阵元位置以及稀布后的阵元位置和激励.因阵元是对称分布的，只给出了第一象限内阵元的位置和激励.

表1 均匀阵元的位置和非均匀阵列阵元的位置与激励

图3 方向图的截面图

图4 阵元位置图

例2： 综合切比雪夫平面阵

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

(b) 重构平面阵方向图图5 综合切比雪夫平面阵的方向图

设有一4×4的切比雪夫平面阵，其阵元均匀分布在矩形栅格上，阵元间距dx=dy=λ/2，要求其环状副瓣的电平为-20 dB.在此例中，共有(2×16+1)(2×16+1)个采样点，增广矩阵束参数K=L=17.按照例1的步骤求得稀布平面阵的阵元位置和激励.图5是切比雪夫方向图和稀布平面阵方向图.旁瓣电平为-16.504 2 dB.图6对比了重构方向图和期望方向图的切面图，进一步增加阵元数目也不会改善方向图的特性.虽然旁瓣电平有小幅抬高(3.495 8 dB)，但可以节省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)方法的减小最小阵元数目、求解阵元位置和设计激励幅度的平面阵列综合方法，与基于随机优化的算法相比，基于增广矩阵束的方法是一种非迭代算法，适合于要求窄波束、低副瓣阵列的设计.另外，与基于矩阵束方法的可分离型分布的平面阵列综合相比，增广矩阵束方法可用于不可分离型分布的平面阵的综合，从而保证方向图在每一剖面的最佳特性.本文的探讨丰富了增广矩阵束方法在稀布平面阵列综合中的应用，为其他种类的稀布阵列综合提供了有益的提示，也为其在工程应用中提供了有价值的参考.

[1] 陈客松,何子述.平面稀布天线阵列的优化算法[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] 郭陈江,张 锋,丁 君,等.基于循环差集与模拟退火法的阵列综合[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)