一种均匀面阵下基于快速收敛平行因子的相干信源测向方法与流程

本发明涉及阵列信号处理技术领域，具体涉及一种均匀面阵下基于快速收敛平行因子的相干信源测向方法。

背景技术：

由于信号在传播中会受到多径效应影响，不可避免会有相干信号产生。研究相干信号和非相干信号一样，都是阵列信号处理中非常重要的内容。相干信号的主要研究内容就是考虑如何去相干。去相干方法最重要的就是在进行空间谱估计前通过一定的方法对已亏损的矩阵秩进行恢复，矩阵秩的亏损是由接收多个相干信号引起的。具体的去相干方法分为两大类：一是空间平滑技术，二是通过阵列的移动或频率的平滑方法等。空间平滑的核心思想是在谱估计前对协方差矩阵进行预处理，使得协方差矩阵的秩恢复到信源个数，然后再利用非相干信号中的相应算法进行信号的角度参数估计。

在传统的相干信源doa估计算法中，使用较为广泛的算法是fbss-esprit和fbss-pm算法，但是这两种算法的信号角度估计性能有限，而传统的paralind分解通过使目标矩阵收敛的方式能够获得较高的估计精度，却有着相当大的计算复杂度。

技术实现要素：

发明目的：为了克服现有技术的不足，本发明通过结合ss-pm算法和paralind算法的方式，可以通过ss-pm算法进行初估计的方式有效提高paralind分解的收敛速度，而且最后的角度估计性能和传统paralind算法十分接近。

技术方案：本发明所述的均匀面阵下基于快速收敛平行因子的相干信源测向方法，该方法包括：

(1)根据均匀面阵中信源的方向矢量建立阵列信号的数学模型；

(2)利用空间平滑结构计算信号的协方差矩阵；

(3)根据所述协方差矩阵，通过ss-pm算法进行信号角度参数初始估计；

(4)构造初始化方向阵，利用paralind分解使得方向矩阵收敛；

(5)估计二维doa。

进一步地，包括：

步骤(1)中，所述阵列信号的数学模型建立过程为：

(11)根据远场目标中非相干目标和相干目标表示均匀面阵中x轴和y轴上信源的方向矢量ax(θk,φk)和ay(θk,φk)；θk为第k个信源对应的仰角，φk表示第k个信源对应的方位角，1≤k≤k，k为信源总数。

(12)计算x轴上m个阵元对应的方向矩阵为ax和y轴上n个阵元对应的方向矩阵为ay；

(13)根据各个子阵的接收信号，确定整个均匀面阵中信源为相干信源时的接收信号向量，所述信源为相干信源的接收信号向量由信源矢量、高斯白噪声、ax、ay以及相关系数矩阵表示。

进一步地，包括：

步骤(2)中，所述计算信号的协方差矩阵具体包括：

(21)根据相干信源的接收信号向量得到第(m，n)个子阵的接收信号向量ymn；

(22)表示第(m，n)个子阵的协方差矩阵rmn，进而计算全部子阵的均值

进一步地，包括：

步骤(3)中，所述通过ss-pm算法进行信号角度参数初始估计，具体包括：

(31)对所述全部子阵的均值进行分块，得到其中，其中，m为面阵中x轴上阵元个数，n为y轴上阵元个数，为复数矩阵的符号，上标即为几行几列的矩阵；

(32)传播算子p的估计值的最小二乘解为：

分别构造矩阵a1，a2，分别表示为：

a1，a2之间相差了一个旋转因子φy，即a2＝a1φy，其中：

其中，d为两个相邻阵元的间距，λ是波长，d1(·),d2(·),...,dm-1(·)表示矩阵a1从第一行到第m-1行构造的对角矩阵，d2(·),d3(·),...,dm(·)表示矩阵a2从第二行到第m行构造的对角矩阵；定义矩阵：

构造矩阵ex＝e(1：m(n-1)，：)，ey＝e(m+1：nm，：)；

其中，ex＝e(1：m(n-1)，：)表示取e的1到m(n-1)行，ey＝e(m+1：nm，：)表示取e的m+1到nm行；

ex、ey可以表示成

其中，为k×k的满秩矩阵；

可得

ey＝ext^-1φyt＝exψ

其中，ψ＝t^-1φyt，ex和ey张成相似的子空间，且矩阵φy的对角元素为ψ的特征值；

(33)根据最小二乘准则，得出，构造的矩阵ψ的估计

对进行特征值分解得到φy的估计值利用的特征向量，得到矩阵t的估计值在无噪声模型下：

其中，п为置换矩阵；

由于与φy的特征值相同，对特征值分解得到uk＝sinθksinφk的估计值通过如下得到：

其中，是矩阵的第k个特征值的估计值，angle(.)表示取复数的相角；

对矩阵e重构，得到

由e′构造矩阵ex′＝e′(1：n(m-1)，：)，ey′＝e′(n+1：mn，：)，其中ex′＝e′(1：n(m-1)，：)表示取e′的1到n(m-1)行，ey′＝e′(n+1：mn，：)表示取e′的n+1到mn行，d1(·),d2(·),...,dn(·)表示矩阵ax从第一行到第n行构造的对角矩阵；

定义

则

进一步得到

(ex′)+ey′＝πφxπ^-1

在无噪声影响时

其中，∏为置换矩阵；

因此，由下式可以得到vk＝cosθksinφk的估计值

其中，εk是矩阵(ex′)+ey′的第k个对角线元素，angle(·)为取复数的相角；表示对仰角估计，为对方位角的估计；

进一步地，包括：

步骤(4)中，所述构造初始化方向阵，利用paralind分解使得方向矩阵收敛，具体包括：

(41)根据和作为已知角带入步骤(11)中得到x轴上m个阵元对应的方向矩阵a′x和y轴上n个阵元对应的方向矩阵a′y，处理后的输出为：

x′＝a′y⊙a′x]γs^t

其中，s表示信源矩阵，γ表示相关矩阵。

(42)将a′y，a′x，γ，s作为paralind分解的初始矩阵，则

定义矩阵表示为

其中，y＝(γs)^t(a′y⊙a′x)^t是无噪接收信号；

(43)对于最小二乘目标函数为：

无噪声模型下，接收信号表示为：

(44)根据最小二乘目标函数更新步骤，将a'x取共轭的最小二乘更新为a′x满秩，是非奇异的，因此得到为：

进一步：

则最小二乘估计可由取复数共轭得到；

(45)根据最小二乘目标函数更新步骤，a′y的估计最小二乘更新为：

进一步地，包括：

步骤(5)中所述估计二维doa的具体步骤包括：

(51)由获得方向矩阵ay和ax精确估计矩阵和后，用和分别代表和的第k列，则：

其中：αxk和ayk分别是ax和ay的第k列，u和v是系数矩阵，nxk和nyk是噪声；

(52)构造两个范德蒙矩阵asx∈c^m×p，asy∈c^n×p，p>>m，p>>n，p表示可能信源数量，m和n分别表示面阵的行和列数；

其中，g是一个采样矢量，矩阵asx和asy可以看作是完备字典；

(53)因此，和可表示为：

xs和ys可以通过范数约束条件得到：

(54)提取出xs，ys中最大模元素的位置作为索引，分别记成px和py，在asx和asy中找到相对应的列，即可得到g(px)和g(py)，它们就是sinθksinφk和sinθkcosφk的估计；

(55)定义γk＝g(px)+jg(py)，仰角和方位角精确值估计可以通过下式得到：

有益效果：本发明采用的技术方案与现有技术相比，具有以下技术效果：(1)该算法通过先对信号角度参数利用ss-pm算法进行初估计，然后利用估计的参数初始化paralind分解的承载矩阵的方式；(2)运算复杂度相对传统paralind算法更低；(3)能有效地用于相干信号的二维doa估计，同时可以得到信号的相干关系矩阵；(4)能够得到配对的方位角和仰角估计，无需额外的参数配对过程；(5)本发明算法的角度估计性能接近于传统的paralind算法，同时优于ss-pm算法和ss-esprit算法。

附图说明

图1是本发明所述的均匀面阵示意图；

图2是本发明所述的均匀面阵的空间平滑结构示意图；

图3a为信噪比为20db时本发明的估计方法得到的角度估计散点图的散点图，图3b为信噪比为5db时本发明的估计方法得到的角度估计散点图的散点图；

图4是本发明所述的算法的角度估计性能在不同快拍数条件下的对比图；

图5是本发明所述的算法的角度估计性能在不同阵元数条件下的对比图；

图6是本发明所述的算法和传统paralind算法的复杂度在相同阵列结构和相同快拍数条件下的对比图；

图7是本发明所述的算法和fbss-pm、fbss-esprit、paralind四种算法的角度估计性能在相同阵列结构和相同快拍数条件下的对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，并不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

符号表示：本发明中用(·)^t表示矩阵转置，大写字母x表示矩阵，小写字母x(·)表示矢量，表示kronecker积，⊙表示khatri-rao积，angle(·)表示取复数的相角，abs(·)表示取复数的模，符号上带尖号表示估计，符号上带星号表示取共轭，()⁺表示矩阵的广义逆矩阵。

本发明中涉及的均匀面阵结构如图1所示，该面阵共有m×n个阵元，均匀分布，两个相邻元素的间距是d，d≤λ/2(λ是波长)。假设空间有k个信源入射到此均匀面阵上，其二维波达方向为(θk，φk)，k＝1,2,…,k，其中，θk，φk分别代表第k个信源的仰角和方位角。定义uk＝sinθksinφk，vk＝cosθksinφk。根据ss-pm的旋转不变性得到角度参数的初始估计，后利用paralind分解不断更新目标矩阵至收敛，最后得出相干信源信号角度参数估计值。本例中基于快速收敛平行因子的相干信源测向方法的具体实现如下：

步骤1：根据均匀面阵中信源的方向矢量建立阵列信号的数学模型：

假定k个远场目标中包括k1个非相干目标以及(k-k1)个相干目标，则均匀面阵中x轴和y轴上信源的方向矢量分别为：

x轴上m个阵元对应的方向矩阵为ax＝[ax(θ1，φ1)，ax(θ2，φ2)，…，ax(θk，φk)]，具体表示为：

y轴上n个阵元对应的方向矩阵为ay＝[ay(θ1，φ1)，ay(θ2，φ2)，…，ay(θk，φk)]，具体表示为

图1中面阵中子阵1的接收信号为

x1(t)＝axs(t)+n1(t)

式中：ax＝[ax(θ1，φ1)，ax(θ2，φ2)，…，ax(θk，φk)]是子阵1的方向矩阵，n1(t)是子阵1的加性高斯白噪声。是信源矢量，为复数矩阵的符号，上标即为几行几列的矩阵。

第n个子阵的接收信号为：

xn(t)＝axφ^n-1s(t)+nn(t)

式中：nn(t)是第n个子阵的加性高斯白噪声。可得整个面阵的接收信号为

也可以表示为：

x(t)＝[ay⊙ax]s(t)+n(t)

当信源为相干信源时：

x(t)＝(ay⊙ax)γs(t)+n(t)

其中，γ为相关系数矩阵。

步骤2：利用空间平滑结构求信号的协方差矩阵：

图2所示的m×n的矩形面阵分为若干个重叠大小为p×q的矩形子面阵。由步骤1得第(m，n)个子阵的接收信号向量为：

其中，1≤m≤m-p+1，1≤n≤n-q+1，a为第(1，1)个子阵的方向矩阵和表示k×k阶对角矩阵φx和φy的n次幂，

则第(m，n)个子阵的协方差矩阵表示成：

故二维空间平滑的协方差矩阵可由全部子阵的均值代替

其中，ms＝m-p+1，ns＝n-q+1。

步骤3：根据协方差矩阵，通过ss-pm算法进行信号角度参数初始估计：

对步骤2中求得的分块可得

其中，

传播算子p的估计值的最小二乘解为：

分别构造矩阵a1，a2

a1，a2之间相差了一个旋转因子φy，即a2＝a1φy，其中

定义矩阵

构造矩阵ex＝e(1：m(n-1)，：)，ey＝e(m+1：nm，：)。其中ex＝e(1：m(n-1)，：)表示取e的1到m(n-1)行，ey＝e(m+1：nm，：)表示取e的m+1到nm行。

ex、ey可以表示成

其中，为k×k的满秩矩阵。

可得

ey＝ext^-1φyt＝exψ

其中，ψ＝t^-1φyt，至此可知，ex和ey张成相似的子空间，且矩阵φy的对角元素为ψ的特征值。

根据最小二乘准则，ψ的估计可由下式得出

对进行特征值分解得到φy的估计值利用的特征向量，得到了矩阵t的估计值

在无噪声模型下

其中，π为置换矩阵。

由于与φy的特征值相同，对特征值分解得到uk＝sinθksinφk的估计值通过如下得到：

其中，是矩阵的第k个特征值，angle(·)表示取复数的相角。

然后对矩阵e重构，得到

由e′构造矩阵ex′＝e′(1：n(m-1)，：)，ey′＝e′(n+1：mn，：)。其中ex′＝e′(1：n(m-1)，：)表示取e′的1到n(m-1)行，ey′＝e′(n+1：mn，：)表示取e′的n+1到mn行。

定义

则

进一步得到

(ex′)+ey′＝πφxπ^-1

在无噪声影响时

其中，п为置换矩阵。

因此由下式可以得到vk＝cosθksinφk的估计值

其中εk是矩阵(ex′)+ey′的第k个对角线元素，angle(.)为取复数的相角。

步骤4：构造初始化方向阵，利用paralind分解使得方向矩阵收敛：

通过ss-pm算法处理得到仰角和方位角的估计值，然后将和作为已知角，代入步骤1中得到x轴上m个阵元对应的方向矩阵a'x和y轴上n个阵元对应的方向矩阵a'y，处理后的输出为：

x′＝[a′y⊙a′x]γs^t

将a′y，a′x，γ，s作为paralind分解的初始矩阵，则

定义可以表示为

其中，y＝(γs)^t(a'y⊙a'x)^t是无噪接收信号。

上面的y表达式数据模型可以看做是一个paralind模型,下面的推导部分为关于paralind分解的推导与证明。

对于最小二乘目标函数为

无噪声模型下，接收信号表示为：

根据步骤4最小二乘目标函数更新步骤，将a'x取共轭的最小二乘更新为

进一步：

则最小二乘估计可由取复数共轭得到；

根据最小二乘目标函数更新步骤，a′y的估计最小二乘更新为：

步骤5：完成二维doa估计：

由获得方向矩阵ay和ax精确估计矩阵和后，用和分别代表和的第k列，则：

其中：αxk和αyk分别是ax和ay的第k列，u和v是系数矩阵，nxk和nyk是噪声；

构造两个范德蒙矩阵asx∈c^m×p，asy∈c^n×p，p>>m，p>>n，p表示可能信源数量，m和n分别表示面阵的行和列数；

其中，g是一个采样矢量，矩阵asx和asy可以看作是完备字典；

因此，和可表示为：

xs和ys可以通过范数约束条件得到：

提取出xs，ys中最大模元素的位置作为索引，分别记成px和py，在asx和asy中找到相对应的列，即可得到g(px)和g(py)，它们就是sinθksinφk和sinθkcosφk的估计；

定义γk＝g(px)+jg(py)，仰角和方位角精确值估计可以通过下式得到：

由于估计矩阵的列都是自动匹配的，所以仰角和方位角也是自动匹配的，该方法无需谱峰搜索和对协方差矩阵进行特征值分解，运算复杂度低；可实现信号角度参数中仰角和方位角配对，避免了传统方法的额外配对。

本发明的方法运算复杂度分析如下：

设均匀面阵的阵元数为m×n，非相干信源数为k1，相干信源数为k-k1，信源总数为k，快拍数为l，本算法的主要复杂度包括：计算ss-pm算法所需要的复杂度和计算paralind分解所需要的复杂度。在ss-pm算法中，构造协方差矩阵平均值所需复杂度为o(mnlp²q²)，特征值分解所需复杂度为o(p³q³)，对和求解所需要的复杂度为o(mn(mn-k)k+k²(mn-k)+2k²(m-1)n+4k³+k²mn+(n-1)mk)而在paralind分解中，所需要的单次分解复杂度为：o(n1(mnl(2k²k1²+2kk1+2k+k1))其中，n1为paralind分解在本发明的方法中的迭代次数。则该算法总复杂度约为：

o(mnlp²q²+p³q³+mn(mn-k)k+k²(mn-k)+2k²(m-1)n+4k³+k²mn+(n-1)mk+n1(mnl(2k²k1²+2kk1+2k+k1))

而传统paralind算法的复杂度约为：

n2(mnl(2k²k1²+2kk1+2k+k1)+mn(2k1²+k²+kk1+2k)+k²(ln+l+n²+n+m+2)+k³k1³+2k³+k1³+kk1l(n+1))+2k²(m+n)+6k²)。由于传统parlaind分解需要的迭代次数n2＞＞n1，所以本发明算法的复杂度远低于传统的paralind算法，这一点从图6也可以得到验证。

图3为当3个相干信源入射均匀面阵时，利用本发明的估计方法得到的角度估计散点图。由图可以看出算法可以有效地估计出信号角度参数。其中，入射信号的角度参数(φ，θ)为：(φ1，θ1)＝(15°，10°)，(φ2，θ2)＝(25°，30°)，(φ3，θ3)＝(35°，50°)；paralind模型的大小为：m×n×l(m＝n＝10，l＝200)，空间平滑子阵大小为p＝q＝3，图3中两个图对应的信噪比为snr＝20db和snr＝5db。

图4是本发明算法角度估计性能在不同快拍下的曲线图。快拍数增加，即采样数据增多。由图可以得出，算法的角度估计性能随着快拍数增加变得更好。其中，入射信号的角度参数：(φ1，θ1)＝(15°，10°)，(φ2，θ2)＝(25°，30°)，(φ3，θ3)＝(35°，50°)；paralind模型的大小为：m×n×l(m＝n＝10)，空间平滑子阵大小为p＝q＝3。

图5是本发明算法角度估计性能在不同阵元下的曲线图。阵元数增加，即分集增益增加。由图可以得出，算法的角度估计性能随着阵元数增加变得更好。其中，入射信号的角度参数：(φ1，θ1)＝(15°，10°)，(φ2，θ2)＝(25°，30°)，(φ3，θ3)＝(35°，50°)；paralind模型的大小为：m×n×l(m＝10，l＝200)，空间平滑子阵大小为p＝q＝3。

图7是本发明算法和fbss-pm、fbss-esprit、paralind四种算法的角度估计性能在相同阵列结构和相同快拍数条件下的曲线图。由图可以得出，本发明算法的角度估计性能优于fbss-pm、fbss-esprit算法，且和传统paralind算法估计性能接近。入射信号的角度参数：(φ1，θ1)＝(15°，10°)，(φ2，θ2)＝(25°，30°)，(φ3，θ3)＝(35°，50°)；paralind模型的大小为：m×n×l(m＝10，l＝200)空间平滑子阵大小为p＝q＝3。