[发明专利]一种基于多发多收阵列的稀疏波达方向角度估计方法

权利要求书

1.一种基于多发多收阵列的稀疏波达方向角度估计方法，其特征在于，该方法包括以下步骤：

S1建立信号模型

考虑一个由G个发射天线和M个接收天线组成的MIMO阵列，收到来自不同方向θ＝[θ₁ θ₂… θ_K]^T的K个目标的回波信号，(·)^T表示矩阵的转置；实际中用到的MIMO阵列，发射天线采用不同的正交编码，因此接收天线通过解码区分开来自不同发射天线的回波信号；于是，来自第i,i＝1,2…G个发射天线的阵列接收信号模型写为

y_i(t)＝Ax_i(t)+n_i(t) (1)

其中下标i用于指代某一序列或者矢量中的第i个元素，y_i(t)表示t时刻接收阵列所接收到的来自第i个发射天线的信号，表示字典矩阵，表示阵列的导向矢量，阵列的导向矢量为远场条件下的平面波入射到阵列所形成的接收信号矢量，它是一个复向量，维数等于阵列的阵元个数；字典矩阵A中的组成了一个关于空间入射角的一种网格划分，N表示该划分网格点的总数，网格划分越细，目标角度估计的分辨率越高；

接收天线阵列采用互质阵的排布方式，设阵元间最小间距为d，则该阵列在以下的位置处存在阵元：0:Pd:(Q-1)Pd，以及Qd:Qd:(2P-1)Qd，其中0:Pd:(Q-1)Pd表示从参考位置每隔Pd放置一个阵元，直到(Q-1)P的位置，Q:Q:(2P-1)Q的定义与之相同，P和Q分别为互质的整数，且PQ；对于这样的接收天线阵列，导向矢量可以表示成：

上式的通项中h代表阵元的坐标位置，即h∈{0:P:(Q-1)P}∪{Q:Q:(2P-1)Q}，λ表示阵列发射的电磁波波长，π为圆周率常数，j表示虚数符号；

公式(1)中，x_i(t)表示阵列接收到来自第i个发射天线的目标信号矢量，它的元素与字典矩阵中的每一列一一对应，当字典矩阵某一列所对应的入射角度存在目标时，与之相对应的x_i(t)的该元素等于目标信号的复数值；当其所对应的角度不存在目标时，x_i(t)的该元素值为0；n_i(t)表示系统的加性噪声；

S2信号协方差矩阵的估计与向量化

由于MIMO阵列具有G个不同的发射信号，因此一共可以形成G个不同的协方差矩阵和长矢量；按照下式计算得到第i个发射信号所对应的的协方差矩阵R_i

上面的公式中(·)^H表示矩阵的共轭转置，L表示用来计算R_i所采用的信号时刻数，该值越大，根据上式所计算的协方差矩阵R_i越接近真实值；

通过对协方差矩阵向量化可以得到

r_i＝vec(R_i) (4)

上式中vec(·)表示矩阵的向量化算子，它是将矩阵按照列由左至右的顺序依次堆叠成矢量；

S3利用互质阵处理方法进行空间自由度拓展

构造投影矩阵并按照下式计算得到投影后等效的线阵接收信号：

通过上面的投影运算，接收信号的维数实质上被降低了，即r_i的维数小于对于不同的发射阵元i，上式所采用的投影矩阵相同；

S4采用近邻算子法来得到目标波达方向估计

S4.1计算投影后等效的线性字典矩阵

不同于S1中的字典矩阵A，上式中是投影后的字典矩阵，它等效于一个由均匀线阵所对应的字典矩阵；对于S1所建立信号模型中所描述的互质阵列排布，所等效的均匀线阵维数等于2PQ+1，等效于一个从-PQd的位置到PQd的位置间隔为d的均匀线阵；因此，中每一列是一个均匀线阵的导向矢量，其维数为2PQ+1，即

S1中的信号模型仅仅是针对当目标的实际波达方向落在字典矩阵的网格点上的情况，当实际波达方向不落在网格点上时，采用S1中的模型所得到的角度估计存在固定的系统误差，这一误差来源于字典矩阵所采用的角度网格是离散的，因此S1中的模型所得到的角度估计只能是实际角度的一种近似；为此，本步骤通过泰勒展开公式来建立目标实际的导向矢量与由S1中网格化模型所得到的近似导向矢量之间的关系；对于字典矩阵计算其一阶导数矩阵

此处引入一阶导数矩阵是为了解决目标不在字典矩阵网格上的问题；

的每一列可以按照下式计算得到

S4.2参数初始化：需要初始化的参数包括近邻算子系数μ，稀疏度约束系数以及Lipschitz常量

近邻算子系数μ的选取需要结合实际的需要，由用户根据自己的系统环境来定，近邻算子系数μ越小，得到的波达方向角估计精度就越高，但越小的近邻算子系数就意味着算法需要迭代更多的次数，即更长的时间才能获得估计结果；

稀疏度约束系数需要根据所估计目标的数量以及所设定网格点的疏密程度来定，目标数量越小，可以设的越大；

Lipschitz常量的取值直接与μ有关，假设MIMO阵列发射了G组不同的正交信号，那么该常量的取值应当满足

值越大，算法收敛需要的时间越长；

S4.3通过迭代的方式得到目标的波达方向估计，需要迭代的中间变量包括z_k和和z_k是与信号直接相关的迭代变量，它们均是维数为2N的矢量，z_k和中的下标k表示第k次迭代，则代表第k次迭代更新前迭代步长，和z_k表示第k次迭代更新前的值，具体步骤如下：

S4.3.1初始化迭代变量：令初始的这里1_2N表示一个元素全为1的向量，其维数为2N；令中间变量迭代步长的初始值迭代计数变量k＝1；

S4.3.2计算当前的误差梯度：令它表示当前的z所对应的模型与接收数据之间的误差，按照下式计算该函数的梯度：

这里的(·)^H表示矩阵的共轭转置；

S4.3.3计算稀疏度约束函数的梯度：定义稀疏度约束函数为h(z_k)＝||z_k||_2,1，这里的表示z_k的联合范数，z_k(i)和z_k(i+N)分别表示z_k的第i项和第i+N项，h(z_k)的Moreau包络记作h_μ(z_k)，按照下面的表达式计算函数h_μ(z_k)的梯度

在上面的表达式中，prox_μh(z_k)表示函数h_μ(z_k)的近邻算子，它的第i项和第i+N项如下所示

上面的表达式中，prox_μh(z_k)_i表示prox_μh(z_k)的第i项，prox_μh(z_k)_i+N表示prox_μh(z_k)的第i+N项；依次令i＝1,2,…N，按照(13)进行计算就可以得到完整的prox_μh(z_k)；

S4.3.4更新通过下面的表达式来更新对的估计

上式中的P_F函数表示投影函数，其输入的自变量为为了给出该函数的数学表达式，将其自变量拆分成前N个元素和后N个元素，并分别记为s和p，并且记s(i)和p(i)分别表示输入向量s和p的第i项，那么根据s和p，函数P_F输出结果的第i项和第i+N项可以写成：

上式中和分别表示函数输出结果的第i项和第i+N项，r表示目标波达方向角度与最近网格点角度之间相差的最大绝对值，参数|·|表示取绝对值；

S4.3.5更新按照如下的表达式更新

上式中等号左边的表示第k次迭代更新后的值，将用于第k+1次迭代，等号右边的表示更新前的值，在第一次迭代中即为初始化的值；

S4.3.6更新z_k+1：按照如下的表达式更新z_k+1

这里等号左边的z_k+1表示第k次迭代更新后的值，表示第k次迭代更新前的值，也即前一次迭代更新后的值；

S4.3.7收敛判断：在更新完z_k+1之后，需要判断迭代是否需要继续，判断的条件如下：如果||z_k+1-z_k||₁＞ε₀并且k+1＜ε₁，那么回到S4.3.2，令k＝k+1，重新开始迭代,；否则退出迭代，执行下一步；阈值参数ε₀和ε₁由用户决定，ε₀越小，ε₁越大，则得到的估计精度越高，但所花费的时间就越多；反之则估计精度降低，但花费的时间较少；

S5确定目标来波方向和数量

目标来波方向和数量的确定通过设定门限ε₂来进行判定得到；记算法输出结果z_k+1的前N项为z_f，后N项为z_b；在算法的输出结果中，z_f的物理含义是信号的幅度，它的各个元素与字典矩阵中的各列相对应，令z_f(i)表示z_f的第i项，那么z_f(i)的值就代表了方向上目标的幅度，z_b是一个与目标实际波达方向相关的物理量，通过它可以求解出目标实际波达方向与网格点之间的差值，通过如下方法来确定来波方向和数量：如果z_f(i)的值超过了门限ε₂，则可判定z_f(i)所对应的来波方向上存在目标，其波达方向角为z_b(i)表示z_b的第i个元素；如果z_f(i)的值低于门限ε₂，则不存在目标；按照这种门限检测的方法在得到了目标来波方向的同时也就确定了目标的数量。