CN102353930B - 一种高精度测向阵列结构设计方法 - Google Patents

Abstract

本发明涉及一种用于波达估计的阵列结构设计方法。通过对阵元位置的优化以实现非等距阵列的布置，包括：确定最小阵元间距d_min，最大阵元位置d_max，所用阵元数N；计算阵列可调布阵距离D＝d_max-(N-1)d_min；对由性能参数G和旁瓣或栅瓣高度来表征的阵列进行优化设计；输出现有最优结果等几大步骤。本发明方法通过优化设计，在保证阵元间距不小于某一特定值条件下，使测向阵列在测向精度与空间谱伪峰高度之间达到最优化，选用该方法设计的合适阵列可以在较少阵元条件下同时实现较低的空间谱伪峰与较高测向精度，并能使阵元间距突破半波长限制，以减小阵元间互耦效应，提高测向精度，也为增大阵元孔径探测更微弱信号预留一定空间，提供了一种行之有效的优化设计方法。

Description

一种高精度测向阵列结构设计方法

技术领域

本发明属于以发射信号可被无方向性接收机检测到并确定发射机相关的方向位置技术领域，涉及一种用于波达方向(Direction of Arrival，DOA)估计的阵列结构设计方法。 

背景技术

波达方向(Direction of Arrival，DOA)估计在电子对抗、水声对抗及新一代无线通信中扮演着非常重要的角色。常用的波达方向估计方法有针对于单辐射源的干涉仪测向与针对于多辐射源的空间谱估计法。这些方法都是依赖于阵列的，阵列结构见附图1所示。干涉仪阵列一般是3～5元非等距阵列，通过对辐射源号到达各阵元时的相位差解模糊估计辐射源的方向。常用的空间谱估计阵列是等距线阵。为了在180°范围内实现无模糊波达方向估计，阵元间距必须小于半波长。这样的阵列主要有如下不足：(1)阵元间互耦效应较大，对波达方向估计不利；(2)阵元尺寸如孔径较大时阵列难以安装，无法使用大孔径阵元接收更微弱的信号；(3)欲实现高精度波达方向估计需增大阵列有效孔径，需要更多的阵元，造价昂贵。为了在大阵列孔径下减小阵元数目，人们已设计出了最小冗余阵列、最大连续延迟阵列、最小间隙阵等非等距线阵，但这些非等距阵列的最小阵元间距仍为半波长，且其它阵元间距为半波长的整数倍，这样的阵列仍无法克服等距线阵的前两个不足。 

发明内容

针对上述现有技术状况，本发明的目的在于，提供一种高精度测向阵列结构设计方法，该方法通过对干涉仪阵列作一定的扩展，可以设计出最小阵元间距大于半波长的非等距阵列，增大阵列孔径并减小阵元间互耦效应，以提高阵列的测向精度，也为增大阵元孔径以探测更微弱信号预留一定的空间。 

现将本发明构思及方法步骤叙述如下： 

本发明的基本构思是通过对阵元位置的优化以实现非等距阵列的布置，阵元位置的优化基于以下技术基础： 

假设阵列是N元线阵，各阵元的位置为d_i，i＝0，1，…，N-1，其中d₀为首阵元位置，其值为0。由空间谱估计理论可知，阵列波达方向估计中常用的MUSIC算法的测向协方差矩阵为： 

$C_{\mathrm{MU}}=\frac{{\mathrm{\σ}}^2}{L}{(H\·I)}^{-1}\mathrm{Re}[H\·{\left(A^H\mathrm{UA}\right)}^T]{(H\·I)}^{-1}---\left(1\right)$

式中，σ²为接收信号中的噪声功率，L为快拍数，I为单位矩阵，A为阵列流形矩阵， <·>表示Hadamard乘积运算，上标H表示矩阵的共轭转置， 

A^HUA＝P^-1+σ²P^-1(A^HA)^-1P^-1                (2) 

$H=D^H{U}_n{U}_n^{H}D=D^H[1-A{\left(A^{H}A\right)}^{-1}{A}^H]D$

D＝[d(β₁)d(β₂)…d(β_m)] 

$d\left({\mathrm{\&beta;}}_i\right)=\frac{\mathrm{d\&alpha;}\left(\mathrm{\&beta;}\right)}{\mathrm{d\&beta;}}{|}_{\mathrm{\&beta;}={\mathrm{\&beta;}}_i}$

其中，β_i＝sinθ_i，θ_i为第i个辐射源的入射角，α(β)为β方向的导向矢量，m为辐射元数目。 

当仅有一个辐射源时，方向的估计方差能说明阵列的测向性能，此时阵列流形矩阵A就为α(β₁)。MUSIC算法的对β的估计方差为： 

${\mathrm{\&sigma;}}_1^2=\frac{{\mathrm{\&sigma;}}^2}{L}{H}^{-1}\mathrm{Re}(P^{-1}+{\mathrm{\&sigma;}}^2{N}^{-1}{P}^{-2})H^{-1}\&Proportional;G^{-1}---\left(3\right)$

式中， 

$G=\underset{i=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{d}_i^2-\frac{1}{N}{\left(\underset{i=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{d}_i\right)}^2---\left(4\right)$

当G取得最大值时，阵列具有最高的测向精度。由理论分析可知，G取最大值时d_i取0或d_max，d_max是阵元设置的最大位置。这种取值方式必定会使阵列测向模糊，且阵元无法安装。 

因此，求取阵元位置d_i值时，应使G取值尽可能大，且还需考虑d_i-d_i-1(i＝1，2，…，N-1)之值大于等于设定门限d_min(如半波长)。更为重要的是，d_i取值时还要考虑伪峰高度。由MUSIC算法可知，空间谱估计式为： 

$P_{\mathrm{MUSIC}}\left(\mathrm{\&beta;}\right)=-10\log \left[\mathrm{\&alpha;}{\left(\mathrm{\&beta;}\right)}^H{U}_n{U}_n^H\mathrm{\&alpha;}\left(\mathrm{\&beta;}\right)\right]=-10\log [N^2-\mathrm{\&alpha;}{\left(\mathrm{\&beta;}\right)}^H{U}_s{U}_s^H\mathrm{\&alpha;}\left(\mathrm{\&beta;}\right)]---\left(5\right)$

式中，U_n为噪声子空间，U_s为信号子空间。 

当β在非波达方向时 与N²接近便出现伪峰。由于信号子空间是阵列流形矩阵A张成的子空间，所以避免伪峰的方法就是使β在非波达方向时导向矢量与信号子空间尽可能正交，即非波达方向导向矢量与波达方向导向矢量的内积绝对值r(Δβ)应尽可能为零。r(Δβ)的值为： 

r(Δβ)＝|α(β)^Hα(β+Δβ)|＝|α(Δβ)|                  (6) 

由式(6)可知，r(Δβ)实际上就是阵列的静态方向图。所以阵列静态方向图的旁瓣(栅瓣)高度决定了伪峰的高度，旁瓣(栅瓣)的高度接近于零时，伪峰高度较低，易于判断向，当旁瓣(栅瓣)接近于主瓣高度时，伪峰高度较高，不易判断波达方向。 

终上所述，一个高精度测向阵列应满足以下三个性能参数条件： 

(1)式(4)取得最大值 

(2)阵列静态方向图的旁瓣最低 

(3)阵元间距不小于设定值 

根据上述结论，在实际阵列设计中，以上三个条件之间相互约束，任何阵列都不可能同时满足以上三个条件。本发明是阵列优化设计方法，通过本方法可以保证所设计的阵列阵元间距不小于指定值，这些阵列结构是在工程应用中，可在所设计的阵列中根据设计要求选择所需的阵列结构。 

因此，本发明一种高精度测向阵列结构设计方法，其特征在于：在满足条件“阵元间距不小于设定值”且性能参数相同的所有阵列中，其静态方向图的旁瓣是最低的参数条件下通过对阵元位置的优化以实现非等距阵列的布置，求得一系列阵列结构，具体设计步骤是： 

步骤1：确定最小阵元间距d_min，最大阵元位置d_max，所用阵元数N； 

步骤2：计算阵列可调布阵距离D＝d_max-(N-1)d_min； 

步骤3：对由性能参数G和旁瓣高度来表征的阵列进行优化设计； 

步骤3.1：设置最大迭代次数M，设置N-2维待优化变量{v₁，v₂，...，v_N-2}，这些变量分别与阵元位置d_i，i＝1，2，...，N-2相对应，设置总迭代次数k＝0，设置现有最优结果数目L＝1，其对应的待优化变量为任意值，对应的阵列结构为任意值，其对应的阵列性能参数G为0，其对应的旁瓣高度为无穷大； 

步骤3.2：用粒子群算法或蚁群算法或遗传算法生成新的变量； 

步骤3.3：根据生成的新变量计算优化的阵元位置； 

步骤3.3.1：对所有待优化变量取绝对值，并计算总和s＝abs(v₁)+abs(v₂)+…+abs(v_N-2)，其中abs()为取绝对值函数； 

步骤3.3.2：计算各阵元间距Δd_i＝D×abs(v_i)/s+d_min，i＝1，2，，N-2； 

步骤3.3.3：计算各阵元位置d_i＝Δd₁+Δd₂+…+Δd_i，i＝1，2，N-2，令d₀＝0，d_N-1＝d_max； 

步骤3.4：由式(4)计算阵列性能参数G，由式(6)的主极大值与次极大值计算方向图旁瓣高度； 

步骤3.5：将步骤3.4取得的所有结果加入现有最优结果中，然后逐个结果对比，保留最优结果。 

步骤3.5.1：计算现有最优结果数目L，设置临时最优结果为空集，设置序号l＝1 

步骤3.5.2：如果现有最优结果中存在一个或多个结果的性能参数G大于第l个结果的性能参数，且它(们)的旁瓣高度低于第l个结果的旁瓣高度，转步骤3.5.4 

步骤3.5.3：将现有最优结果中第l个结果加入临时最优结果 

步骤3.5.4：令l＝l+1，如果l＜＝L转步骤3.5.2 

步骤3.5.5：用临时最优结果替换现有最优结果 

步骤3.6：如果觉得已取得最优结果过多，可采用小生境技术在保证多样性的同时缩减已取得的最优结果； 

步骤3.7：令k＝k+1，如果k＜M，转步骤3.2； 

步骤4：输出现有最优结果 

本方法中，步骤3.3保证了所设计阵列的阵元间距大于设定的最小阵元间距，步骤3.5确保各个已取得的最优结果不差于其它结果。 

本发明同现有技术相比，通过优化设计，在保证阵元间距不小于某一特定值条件下，使测向阵列在测向精度与空间谱伪峰高度之间达到最优化，选用该方法设计的合适阵列可以在较少阵元条件下同时实现较低的空间谱伪峰与较高测向精度，并能使阵元间距突破半波长限制，以减小阵元间互耦效应，提高测向精度，也为增大阵元孔径探测更微弱信号预留一定空间，为波达方向估计的阵列结构设计提供了一种行之有效的优化设计方法。 

附图说明

图1本发明高精度测向阵列结构设计方法阵列设计程序流程图 

图2非等距线阵结构示意图 

图3实施例1利用本发明设计的阵列性能关系图 

图4实施例2设计阵列1估计的空间谱 

图5实施例2设计阵列2估计的空间谱 

图6实施例2设计阵列3估计的空间谱 

图7实施例3中各阵列估计成功率关系图 

图8实施例3中对-45°方向的估计偏差关系图 

图9实施例3中对-30°方向的估计均方差关系图 

图10实施例3中对-45°方向的估计偏差关系图 

图11实施例3中对-30°方向的估计均方差关系图 

具体实施方式

实施例1 

参见图：3：根据本发明的阵列结构设计方法，按照本实施例1的要求，设计出一系列阵列结构，然后选择测向性能不同的三个阵列进行空间谱估计，比较其性能。假设需设计的阵列的阵元数N＝8，最大阵元位置d_max＝50λ(λ为信号波长)，最小阵元间距d_min＝1.5λ，设置最大迭代次数M，采用粒子群算法，按本发明具体设计步骤进行优化设计。由步骤4输出最优结果的性能关系如图3所示。图中，由步骤3得到的性能参数G为纵坐标，旁瓣高度为横坐标得到一曲线图。根据曲线上任意一点所对应的阵元位置d_i所确定的阵列测向性能为该点的G值和旁瓣高度值。该曲线可供设计人员对阵列的优化设计进行选择。 

实施例2 

参见图3～6：对所设计阵列测向性能进行仿真实验验证：假设有功率相等的两个信号从相近的两个方向θ₁＝-45°与θ₂＝-30°入射，信噪比为10dB(参见式(1)式(2))。采用MUSIC算法，利用实施例1设计出的三个阵列即：测向精度最高阵列(阵列1)、测向精度居中阵列(阵列2)，测向精度最差阵列(阵列3)由式(5)进行空间谱估计，设计出的阵列结构见下表所示，三个阵列估计的空间谱如附图3～5所示。 

表一阵列结构表 

阵元位置	d0	d1	d2	d3	d4	d5	d6	d7
									阵列1(×λ)	0	1.5016	3.0031	4.5807	45.4830	46.9950	48.4956	50
阵列2(×λ)	0	1.5060	3.0061	5.0088	37.6822	46.7049	48.4803	50
									阵列3(×λ)	0	1.5556	3.9477	11.9312	28.6029	44.9379	48.0375	50

从三个阵列估计出的空间谱可以看出，阵列1的伪峰较高，如果预先不知道只有两个辐射源，根据空间谱可能会认为存在多个辐射源；阵列2的局部空间谱较为尖锐，伪峰高度较低，但也会存在认为有多个辐射源的可能；阵列3的伪峰高度最低，通常能正确地判断出辐射源的方向与数目。 

实施例3 

参见图7～11：对表一中的三个阵列的目标分辨成功概率、分辨率、测向偏差与测向均方差性能进行测试。假设有功率相等的两个信号从相近的两个方向-45°与-30°入射，采用MUSIC算法进行空间谱估计，针对各阵列进行1000次独立的仿真，如果根据最高的两个谱峰能估计出两个辐射源，且估计偏差不大于2°，认为算法成功，对成功估计结果进行统计。各阵列估计成功率见附图6所示，各阵列估计偏差见附图7，8所示，各阵列估计均方差见图9，10所示。 

从该实施例可看出，在限定最小阵元间距情况下，阵列测向精度与伪峰高度不能同时达到最优，当测向精度高时，伪峰也较高，当测向精度低时，伪峰也较低。因此，在实际应用中应根据阵列布阵位置从该方法设计的阵列结构中根据需要选择合适的阵列结构。 

Claims (3)

1.一种高精度测向阵列结构设计方法，其特征在于：在满足条件“阵元间距不小于设定值”且性能参数相同的所有阵列中，通过对阵元位置的优化以实现非等距阵列的布置，求得一系列阵列结构，使其静态方向图的旁瓣最低，其具体步骤为： 

步骤1：确定最小阵元间距d_min，最大阵元位置d_max，所用阵元数N； 

步骤2：计算阵列可调布阵距离D=d_max-(N-1)d_min； 

步骤3：对由性能参数G和旁瓣或栅瓣高度来表征的阵列进行优化设计； 

步骤3.1：设置最大迭代次数M，设置N-2维待优化变量{v₁，v₂，...，v_N-2}，这些变量分别与阵元位置d_i，i=1，2，...，N-2相对应，设置总迭代次数k=0，设置现有最优结果数目L=1，其对应的待优化变量为任意值，对应的阵列结构为任意值，其对应的阵列性能参数G为0，其对应的旁瓣或栅瓣高度为无穷大； 

步骤3.2：用粒子群算法或蚁群算法或遗传算法生成新的变量； 

步骤3.3：根据生成的新变量计算优化的阵元位置； 

步骤3.4：由式(4)计算阵列性能参数G，由式(6)的主极大值与次极大值计算方向图旁瓣高度； 

r(Δβ)=|α(β)^Hα(β+Δβ)|=|α(Δβ)|   (6) 

式中：G为性能参数；d_i为阵元i位置；N为阵元数；r(Δβ)为阵列静态方向图；β＝sinθ，θ为辐射源的入射角；α(β)为β方向的导向矢量，上标H表示矩阵的共轭转置； 

步骤3.5：将步骤3.4取得的所有结果加入现有最优结果中，然后逐个结果对比，保留最优结果； 

步骤3.6：令k=k+1，如果k<M，转步骤3.2； 

步骤3.7：如果觉得已取得最优结果过多，可采用小生境技术在保证多样性的同时缩减已取得的最优结果；

步骤4：输出现有最优结果。 

2.根据权利要求1所述的一种高精度测向阵列结构设计方法，其特征在于：步骤3.3中所述的“根据生成的新变量计算优化的阵元位置”的具体步骤为： 

步骤3.3.1：对所有待优化变量取绝对值，并计算总和s=abs(v₁)+abs(v₂)+…+abs(v_N-2)，其中abs()为取绝对值函数； 

步骤3.3.2：计算各阵元间距Δd_i＝D×abs(v_i)/s+d_min，i=1，2，……，N-2； 

步骤3.3.3：计算各阵元位置d_i＝Δd₁+Δd₂+…+Δd_i，i=1，2，……，N-2，令d₀=0，d_N-1=d_max。 

3.根据权利要求2所述的一种高精度测向阵列结构设计方法，其特征在于：步骤3.5中所述的“将步骤3.4取得的所有结果加入现有最优结果中，然后逐个结果对比，保留最优结果”的具体步骤为： 

步骤3.5.1：计算现有最优结果数目L，设置临时最优结果为空集，设置序号l=1； 

步骤3.5.2：如果现有最优结果中存在一个或多个结果的性能参数G大于第l个结果的性能参数，且它或它们的旁瓣或栅瓣高度低于第l个结果的旁瓣或栅瓣高度，转步骤3.5.4； 

步骤3.5.3：将现有最优结果中第l个结果加入临时最优结果； 

步骤3.5.4：令l=l+1，如果l<=L转步骤3.5.2； 

步骤3.5.5：用临时最优结果替换现有最优结果。 

Images