CN103731189B - 共形天线阵动态子阵划分方法及波达方向估计方法 - Google Patents

Abstract

本发明公开一种共形天线阵动态子阵划分方法及波达方向估计方法，先确定k个初始聚类中心；根据各阵元与这些聚类中心的相似度（距离），分别将他们分配给与其最相似的子阵聚类（聚类中心所代表的）中。然后再计算每个新聚类的聚类中心（该聚类中心所有对象的均值）。不断重复这个过程直到标准测度函数开始收敛为止。本方法解决了天线载体在DOA估计中会发生转动的问题。

Description

共形天线阵动态子阵划分方法及波达方向估计方法

技术领域

本发明属于阵列信号处理领域，具体涉及一种共形天线阵动态子阵划分方法及波达方向估计方法。

背景技术

阵列信号处理将一组传感器按一定分布方式布置在空间不同位置上，形成传感器阵列，用传感器阵列来接收空间目标信号，并对接收到的信号进行处理。它作为信号处理的一个重要分支，在通信、雷达、声纳、地震勘探和射电天文等领域内获得了广泛应用和迅速发展。自适应空域滤波和空间谱估计是阵列信号处理最主要的两个研究方向，其中空间谱估计，又称为波达方向（Direction Of Arrival，DOA）估计，是指对空间信号波达方向分布进行的超分辨估计，其主要目的是估计信号的空域参数和信源位置。基于共形天线的DOA估计研究是信号处理和电磁场与微波技术学科的交叉结合，是一门新兴技术，也是未来提高通讯性能的主要方法之一。

共形天线是一种和物体外形保持一致的天线，它的应用可以在保证天线性能的条件下简化安装，节省空间，最大限度地减少对天线阵载体空气动力学性能的影响，因此关于共形天线阵的研究近年来日益得到关注。由于共形天线一般安装在曲面上，载体的曲率变化和遮挡效应会使共形天线阵在扫描时，阵列单元难以保持相同的方向图形状，甚至让部分阵元接收不到某一方向的入射信号。然而，目前大多数高性能的谱估计算法，如MUSIC类算法、ESPRIT类算法，都假设阵列导向矢量完整（即所有阵元都能接收到某个入射信号），所以，经典的DOA估计方法不能直接应用在共形阵列上。

针对上述问题，国内外学者提出了多种解决方案。首先是杨鹏等人在“MUSIC算法在柱面共形天线阵DOA估计中的应用研究”一文中提出将整个阵列进行子阵划分，克服圆柱体遮挡效应对MUSIC算法估计性能的影响，成功的在柱面共形阵上实现了信号的DOA估计。与子阵分割MUSIC算法进行对比，杨永建等人的“基于MUSIC算法的圆柱共形阵DOA估计”一文在考虑载体遮挡效应的同时，结合阵元的方向性，通过对导向矢量进行重构，解决了导向矢量与噪声子空间不完全正交这一问题。可这些工作均以天线载体与入射信源保持相对静止为前提条件，在实际工程应用中，如飞机、气球等空中平台的共形天线阵列，由于探测需要或者外力干扰等因素，天线载体很难一直固定在同一位置，此时上述方法很难达到预期效果。

发明内容

本发明所要解决的是天线载体在波达方向估计时发生转动的问题，而提供一种共形天线阵动态子阵划分方法及波达方向估计方法。

为解决上述问题，本发明是通过以下技术方案实现的：

一种共形天线阵动态子阵划分方法，包括如下步骤：

步骤1，根据共形天线载体的形状特点和工程应用要求，将该共形天线载体所覆盖的目标空域分为k块聚类区；

步骤2，选择每块聚类区的共形载体表面中心作为该块聚类区的初始聚类中心；

步骤3，检测共形天线载体是否发转动，若转动，则根据共形天线载体的转动情况对共形天线载体表面上的所有阵元进行实时定位；若不转动，继续用前面时刻的子阵划分状态；

步骤4，计算每个阵元到各个虚拟聚类中心的距离，并把阵元分配到距离最近的虚拟聚类中心所代表的子阵聚类中；

步骤5，所有阵元分配完成后，将每个子阵聚类中的所有阵元的坐标之和除以阵元数以获得每个子阵聚类的几何聚类中心；

步骤6，计算每个阵元到各个几何聚类中心的距离，并把阵元分配到距离最近的几何聚类中心所代表的子阵聚类中；

步骤7，重新划分之后，若步骤6所得子阵聚类划分结果的标准测度函数未收敛，则需要返回步骤2重新选定另一点作为该块聚类区的新的虚拟聚类中心，并重新进行子阵聚类划分，直到标准测度函数收敛为止；若步骤6所得子阵聚类划分结果的标准测度函数开始收敛，则输出步骤6所得子阵聚类划分结果；上式中，p代表天线载体上任意一阵元对象的坐标；m_i表示第i个聚类中心坐标。

在步骤4和6中，若某一阵元与2个或2个以上的聚类中心的距离相等，则把这个阵元同时分配到这些2个或2个以上的聚类中心所代表的聚类中。

一种共形天线阵波达方向估计方法，包括如下步骤：

步骤1，根据共形天线载体的形状特点和工程应用要求，将该共形天线载体所覆盖的目标空域分为k块聚类区；

步骤2，选择每块聚类区的共形载体表面中心作为该块聚类区的初始聚类中心；

步骤3，检测共形天线载体是否发转动，若转动，则根据共形天线载体的转动情况对共形天线载体表面上的所有阵元进行实时定位；若不转动，继续用前面时刻的子阵划分状态；

步骤4，计算每个阵元到各个虚拟聚类中心的距离，并把阵元分配到距离最近的虚拟聚类中心所代表的子阵聚类中；

步骤5，所有阵元分配完成后，将每个子阵聚类中的所有阵元的坐标之和除以阵元数以获得每个子阵聚类的几何聚类中心；

步骤6，计算每个阵元到各个几何聚类中心的距离，并把阵元分配到距离最近的几何聚类中心所代表的子阵聚类中；

步骤7，重新划分之后，若步骤6所得子阵聚类划分结果的标准测度函数未收敛，则需要返回步骤2重新选定另一点作为该块聚类区的新的虚拟聚类中心，并重新进行子阵聚类划分，直到标准测度函数收敛为止；若步骤6所得子阵聚类划分结果的标准测度函数开始收敛，则输出步骤6所得子阵聚类划分结果；上式中，p代表天线载体上任意一阵元对象的坐标；m_i表示第i个聚类中心坐标；

步骤8，分别在每个子阵聚类上进行波达方向估计。

上述方案中，所述步骤8具体为：

步骤8.1，在同一个坐标系中依据建立各子阵聚类的接收数据模型；上式中，X代表共形阵列天线的接收数据矢量，S代表信号矢量，N代表噪声矢量，代表阵列导向矢量矩阵，由构成，表示第i个信号从方向入射至共形阵列的导向矢量，分别为第i个入射信号的方位角和俯仰角，K为入射信源个数，

步骤8.2，从子阵接收数据模型的协方差矩阵中提取出信号子空间U_S和噪声子空间U_N；上式中，R代表子阵接收数据矢量的协方差矩阵，X_n代表第n次快拍的接收数据，L代表快拍数，上标H代表共轭转置；

步骤8.3，使变化，根据高分辨率波达方向估计算法即MUSIC算法的谱估计公式进行谱峰搜索，通过寻求谱估计公式的峰值完成第k个子阵的波达方向估计，其峰值对应的角度即为信源的入射方向；上式中，表示共形天线覆盖领域内各方向的空间谱，表示任意一方向的共形阵导向矢量，U_N是由从子阵接收数据矢量协方差矩阵R中提取出来的噪声子空间，上标H代表共轭转置。

在步骤4和6中，若某一阵元与2个或2个以上的聚类中心的距离相等，则把这个阵元同时分配到这些2个或2个以上的聚类中心所代表的聚类中。

本发明将K-means聚类算法应用到子阵划分技术中，其核心思想是把共形天线载体表面的n个阵元数据对象划分为k个聚类（即k个子阵），使每个子阵聚类中的阵元到该聚类中心的平方和最小。与现有技术相比，本发明解决了空中平台上天线载体发生转动导致各阵元与信源位置产生相对移动这一实际问题。同时，该方法适用于各种几何形态的天线载体，不影响经典DOA估计算法的分辨力和估计精度。

附图说明

图1为基于K-means算法的动态子阵划分方法的流程图。

图2为球面载体阵列天线模型图。

图3为球面阵信源入射示意图。

图4为单元子阵仿真对比图。（a）为固定子阵划分，（b）为动态子阵划分。

图5为球面共形阵列天线DOA估计仿真图。

具体实施方式

一种共形天线阵波达方向估计方法，包括共形天线阵动态子阵划分的过程和在各子阵聚类上进行波达方向估计的过程。其中共形天线阵动态子阵划分的过程为本发明的核心点，而在各子阵聚类上进行波达方向估计的过程与现有技术相同或相近似。

①上述共形天线阵动态子阵划分的过程如下（参见图1）：

步骤1，根据共形天线载体的形状特点和工程应用要求，将该共形天线载体所覆盖的目标空域分为k块聚类区。

k块子阵聚类的阵元会随着共形阵列天线载体的转动发生相应调整，聚类块的空间位置是固定不变的，划分方案由共形阵列天线载体的形态特点和其覆盖的待测空域范围制定，方案确定后需保证每个子阵的所有阵元可以接收到该子阵对应空域任何一个方向入射来的信号。一般而言，若天线载体为平面立体图形，即该载体为n面时，则n≤k≤2n；若天线载体为曲面立体图形，则要求k值使每块分界线在同一圆周平面的切线夹角≤90°。以图2所示球面共形阵为例，在已建立的坐标系中，将该球面阵覆盖的空域分为四块，分别为：x≥0,y≥0,z≥0，x≥0,y≤0,z≥0，x≤0,y≤0,z≥0，x≤0,y≥0,z≥0，划分的子阵聚类则为相应空域块内的球面阵元。

步骤2，选择每块聚类区的共形载体表面中心作为该块聚类区的初始聚类中心。若某一阵元与2个或2个以上的聚类中心的距离相等，则把这个阵元同时分配到这些2个或2个以上的聚类中心所代表的聚类中。

初始聚类中心为一虚拟点，该处可以有阵元，也可以没有，初始聚类中心是由步骤1确定的聚类块决定，选取各子阵聚类块边界范围的中点坐标值；

步骤3，检测共形天线载体是否发转动，若转动，则根据共形天线载体的转动情况对共形天线载体表面上的所有阵元进行实时定位。若不转动，继续用前面时刻的子阵划分状态。

步骤4，计算每个阵元到各个虚拟聚类中心的距离，并把阵元分配到距离最近的虚拟聚类中心所代表的子阵聚类中。

步骤5，所有阵元分配完成后，将每个子阵聚类中的所有阵元的坐标之和除以阵元数以获得每个子阵聚类的几何聚类中心。

步骤6，计算每个阵元到各个几何聚类中心的距离，并把阵元分配到距离最近的几何聚类中心所代表的子阵聚类中。若某一阵元与2个或2个以上的聚类中心的距离相等，则把这个阵元同时分配到这些2个或2个以上的聚类中心所代表的聚类中。

步骤7，重新划分之后，若步骤6所得子阵聚类划分结果的标准测度函数未收敛，则需要返回步骤2重新选定另一点作为该块聚类区的新的虚拟聚类中心，并重新进行子阵聚类划分，直到标准测度函数收敛为止；若步骤6所得子阵聚类划分结果的标准测度函数开始收敛，则输出步骤6所得子阵聚类划分结果；上式中，p代表天线载体上任意一阵元对象的坐标；m_i表示第i个聚类中心坐标。

②在各子阵聚类上进行波达方向估计的过程如下：

步骤8.1，在同一个坐标系中依据建立各子阵聚类的接收数据模型；上式中，X代表共形阵列天线的接收数据矢量，S代表信号矢量，N代表噪声矢量，代表阵列导向矢量矩阵，由构成，表示第i个信号从方向入射至共形阵列的导向矢量，分别为第i个入射信号的方位角和俯仰角，K为入射信源个数，

步骤8.2，从子阵接收数据模型的协方差矩阵中提取出信号子空间U_S和噪声子空间U_N；上式中，R代表子阵接收数据矢量的协方差矩阵，X_n代表第n次快拍的接收数据，L代表快拍数，上标H代表共轭转置；

步骤8.3，使变化，根据高分辨率波达方向估计算法即MUSIC算法的谱估计公式进行谱峰搜索，通过寻求谱估计公式的峰值完成第k个子阵的波达方向估计，其峰值对应的角度即为信源的入射方向；上式中，表示共形天线覆盖领域内各方向的空间谱，表示任意一方向的共形阵导向矢量，U_N是由从子阵接收数据矢量协方差矩阵R中提取出来的噪声子空间，上标H代表共轭转置。

本实施例以一置于空中平台的球面阵为研究对象，如图2所示，该球面阵一共包含25个阵元，分为3层，第一层阵元数为1,第二层阵元数为8,第三层阵元数为16，每层阵元均匀分布在相应的圆周上。

首先，依照图2所示方向建立坐标系，将该球面阵覆盖的空域分为4个部分，每部分选择一恰当的参考点（见表1聚类中心坐标）作为相应子阵的虚拟聚类中心，然后检测球面载体是否发生转动，若发生，根据其转动特点对载体表面的阵元进行定位。最后将球面各阵元作为数据对象依照标准测度函数将阵元进行聚类划分。表1为初始位置球面阵应用K-means聚类算法的子阵划分情况，表2为球面阵转动后的子阵划分情况，表3为球面阵转动后的子阵划分情况。从三个表中可以看到，本发明可以根据空中平台上的球面共形载体转动情况精确地对载体上的所有阵元进行即时子阵划分。

表1球面阵转动前子阵划分状态

表2球面阵转动23°子阵划分状态

表3球面阵转动45°后子阵划分状态

基于本发明的转动球面共形天线阵DOA估计中的效果可以通过以下实验仿真进行说明。

第一步，建立球面共形阵接收数据模型，图3所示为球面阵信源入射示意图，其坐标轴正方向如图所示，设球面载体的半径为r，图中“×”代表球面阵任意一阵元，θ为阵元所在点与z轴正方向所成的角，为该阵元所在点与x轴正方向所成的角，第m个阵元的位置为以坐标原点为参考点，波长为λ₀的N个不相关同频窄带信号分别从方向入射到球面共形天线阵，其中φ_i和分别表示第i个入射信号的方位角和俯仰角，R_0i和R_mi分别为该信号到参考点和第m个阵元的距离，令τ_mi为第i个信号到第m个阵元相对于该信号到参考点的时延，则有：

${\mathrm{\τ}}_{\mathrm{mi}}=\frac{R_{0i}-R_{\mathrm{mi}}}{c}=-\frac{\mathrm{\Δ}{R}_i}{c}$

导向矢量为：

式中：c为光速。

当K个信号分别从这些方向入射到第k个子阵时，该子阵的接收数据可以表示为：X、S、N分别为子阵接收数据矢量、信号矢量和噪声矢量。

第二步，根据上文K个接收信号矢量得到协方差矩阵的估计值，即X_n表示相应子阵的第n次快拍接收数据。对协方差矩阵进行特征值分解，得到R=UΣU^H，其中U是由R的特征向量构成，Σ由R的特征值构成。

第三步，按照特征值的大小顺序，把与该子阵区域信号个数K相等的大特征值和对应的特征矢量看作信号子空间U_S，把剩下的M_k-K（M_k为第k个子阵阵元数）个特征值和对应的特征矢量看作噪声子空间U_N，则有： $R=U_s{\mathrm{\Σ}}_s{U_s}^H+U_N{\mathrm{\Σ}}_N{U_N}^H.$

第四步，使变化，按照来计算谱函数，通过寻求峰值完成第k个子阵的DOA估计，其峰值对应的角度即为信源的入射方向。上式中，表示共形天线覆盖领域内各方向的空间谱，表示任意一方向的共形阵导向矢量，U_N是由从子阵接收数据矢量协方差矩阵R中提取出来的噪声子空间；

图4和图5是中心频率为1KHz的同频窄带信源，噪声选取为全向高斯白噪声，快拍数为1024，信噪比取20dB，信源入射角随机选取条件下按照上述步骤进行DOA估计的matlab仿真图。在球面共形阵相对于初始位置进行23度转动后，用固定子阵划分法（图4（a））和基于K-means算法的共形天线阵动态子阵划分方法（图4（b））的子阵单元DOA估计仿真图如图4所示，从该图中可以看到，在球面阵幅度稍小的转动情况下两种方案均能探测出信源的DOA，显然（b），即在应用基于K-means算法的共形天线阵动态子阵划分方法的子阵单元仿真图中，信源DOA估计精确度更高，谱峰较为尖锐；当球面载体发生稍大幅度的转动，45度甚至更大时，就不能保证对于某一入射信号，对应子阵单元的所有阵元可以一直接收到该信号。图5为基于K-means算法的共形天线阵动态子阵划分方法在整个球面阵上的应用仿真图，该仿真中一共有4个信源，入射方向分别为：(30°,10°)，(120°,60°)，(200°,35°)，(300°,75°)。从该仿真图中可以看到本发明应用到转动的球面共形阵列天线上可以准确地估计出多个信源DOA。

综上所述，本发明提出的基于K-means算法的共形天线阵动态子阵划分方法解决了实际工程应用中共形天线载体因各类因素发生转动这一问题，所选实施例对本发明技术方案的详细说明是示意性的而非限制性的，基于这样的理解，本领域技术人员对本发明所记载的技术方案进行修改或者对部分技术特征的等同替换并不能使相应技术方案的本质脱离本发明的精神和范围。

Claims (5)

1.共形天线阵动态子阵划分方法，其特征是，包括如下步骤：

步骤1，根据共形天线载体的形状特点和工程应用要求，将该共形天线载体所覆盖的目标空域分为k块聚类区；

步骤2，选择每块聚类区的共形载体表面中心作为该块聚类区的初始聚类中心；

步骤3，检测共形天线载体是否发生转动，若转动，则根据共形天线载体的转动情况对共形天线载体表面上的所有阵元进行实时定位；若不转动，继续用前面时刻的子阵划分状态；

步骤4，计算每个阵元到各个虚拟聚类中心的距离，并把阵元分配到距离最近的虚拟聚类中心所代表的子阵聚类中；

步骤5，所有阵元分配完成后，将每个子阵聚类中的所有阵元的坐标之和除以阵元数以获得每个子阵聚类的几何聚类中心；

步骤6，计算每个阵元到各个几何聚类中心的距离，并把阵元分配到距离最近的几何聚类中心所代表的子阵聚类中；

步骤7，重新划分之后，若步骤6所得子阵聚类划分结果的标准测度函数未收敛，则需要返回步骤2重新选定另一点作为该块聚类区的新的虚拟聚类中心，并重新进行子阵聚类划分，直到标准测度函数收敛为止；若步骤6所得子阵聚类划分结果的标准测度函数开始收敛，则输出步骤6所得子阵聚类划分结果；上式中，p代表天线载体上任意一阵元对象的坐标；m_i表示第i个聚类中心坐标；C_i表示第i个聚类，k表示聚类的块数。

2.根据权利要求1所述的共形天线阵动态子阵划分方法，其特征是，在步骤4和6中，若某一阵元与2个或2个以上的聚类中心的距离相等，则把这个阵元同时分配到这些2个或2个以上的聚类中心所代表的聚类中。

3.共形天线阵波达方向估计方法，其特征是，包括如下步骤：

步骤1，根据共形天线载体的形状特点和工程应用要求，将该共形天线载体所覆盖的目标空域分为k块聚类区；

步骤2，选择每块聚类区的共形载体表面中心作为该块聚类区的初始聚类中心；

步骤3，检测共形天线载体是否发生转动，若转动，则根据共形天线载体的转动情况对共形天线载体表面上的所有阵元进行实时定位；若不转动，继续用前面时刻的子阵划分状态；

步骤4，计算每个阵元到各个虚拟聚类中心的距离，并把阵元分配到距离最近的虚拟聚类中心所代表的子阵聚类中；

步骤5，所有阵元分配完成后，将每个子阵聚类中的所有阵元的坐标之和除以阵元数以获得每个子阵聚类的几何聚类中心；

步骤6，计算每个阵元到各个几何聚类中心的距离，并把阵元分配到距离最近的几何聚类中心所代表的子阵聚类中；

步骤7，重新划分之后，若步骤6所得子阵聚类划分结果的标准测度函数未收敛，则需要返回步骤2重新选定另一点作为该块聚类区的新的虚拟聚类中心，并重新进行子阵聚类划分，直到标准测度函数收敛为止；若步骤6所得子阵聚类划分结果的标准测度函数开始收敛，则输出步骤6所得子阵聚类划分结果；上式中，p代表天线载体上任意一阵元对象的坐标；m_i表示第i个聚类中心坐标；C_i表示第i个聚类，k表示聚类的块数；

步骤8，分别在每个子阵聚类上进行波达方向估计。

4.根据权利要求3所述的共形天线阵波达方向估计方法，其特征是，上述步骤8具体为：

步骤8.1，在同一个坐标系中依据X＝A(φ,θ)S+N建立各子阵聚类的接收数据模型；上式中，X代表共形阵列天线的接收数据矢量，S代表信号矢量，N代表噪声矢量，A(φ,θ)代表阵列导向矢量矩阵，由a(φ_i,θ_i)(i＝1,2,…,n)构成，a(φ_i,θ_i)表示第i个信号从方向(φ_i,θ_i)入射至共形阵列的导向矢量，φ_i,θ_i分别为第i个入射信号的方位角和俯仰角，K为入射信源个数，A(φ,θ)＝[a(φ₁,θ₁),a(φ₂,θ₂),...,a(φ_K,θ_K)]；

步骤8.2，从子阵接收数据模型的协方差矩阵中提取出信号子空间U_S和噪声子空间U_N；上式中，R代表子阵接收数据矢量的协方差矩阵，X_n代表第n次快拍的接收数据，L代表快拍数，上标H代表共轭转置；

步骤8.3，使(φ,θ)变化，根据高分辨率波达方向估计算法即MUSIC算法的谱估计公式进行谱峰搜索，通过寻求谱估计公式的峰值完成第k个子阵的波达方向估计，其峰值对应的角度即为信源的入射方向；上式中，P_MUSIC(φ,θ)表示共形天线覆盖领域内各方向的空间谱，a(φ,θ)表示任意一方向(φ,θ)的共形阵导向矢量，U_N是由从子阵接收数据矢量协方差矩阵R中提取出来的噪声子空间，上标H代表共轭转置。

5.根据权利要求3所述的共形天线阵波达方向估计方法，其特征是，在步骤4和6中，若某一阵元与2个或2个以上的聚类中心的距离相等，则把这个阵元同时分配到这些2个或2个以上的聚类中心所代表的聚类中。