基于张量链分解的稀疏立方阵列波达方向估计方法

1.本发明属于阵列信号处理技术领域，尤其涉及基于稀疏阵列虚拟域张量的统计信号处理技术，具体是一种基于张量链分解的稀疏立方阵列波达方向估计方法，可用于目标定位。


背景技术：

2.相比于传统遵循奈奎斯特采样速率的均匀阵列，稀疏阵列具有大孔径、高分辨率的优势，能够突破传统均匀阵列波达方向估计在估计性能与成本开销上的性能瓶颈。具有系统化结构的稀疏阵列基于接收信号的二阶统计量推导增广虚拟阵列，从而基于虚拟域二阶等价信号实现奈奎斯特匹配的波达方向估计。嵌套互质阵列作为一种典型的稀疏阵列架构，能够推导出连续的虚拟阵列，因此得到广泛应用。然而，传统的虚拟域信号处理方法将接收信号表示成矢量，并通过矢量化接收信号的协方差矩阵推导虚拟域二阶等价信号；随着实际应用中稀疏阵列维度的不断扩张，多维稀疏阵列的接收信号涵盖多维度时空信息，而这种矢量化信号的处理方法将破坏多维接收信号的原始结构，造成严重的性能损失。
3.为了保留多维接收信号的结构化信息，张量作为一种多维的数据类型，开始被应用于阵列信号处理领域，用于表征涵盖多维时空信息的接收信号。现有面向多维稀疏阵列的张量信号处理方法通过推导二阶虚拟域张量，将对应的高维虚拟域协方差张量表示为canonical polyadic(cp)模型，对其进行直接canonical polyadic分解以实现高精度、高分辨的波达方向估计。然而，canonical polyadic分解在处理高维张量的时候，存在收敛速度慢、计算复杂度高的问题。而张量链分解作为另外一种张量分解形式，能够将高维张量高效分解为若干低维度核张量和矩阵的链式连接形式，进而针对低维核张量做后续处理，能够在保障计算效率的前提下有效提取信号特征，从而实现信道估计、谐波检索、高光谱图像去噪等具体功能。由此可见，张量链分解在提升稀疏阵列波达方向估计的处理效率方面具有极大潜力。针对高维虚拟域协方差张量，如何运用张量链分解设计一种低计算复杂度的波达方向估计方法仍是一个亟待解决的问题。


技术实现要素：

4.本发明的目的在于针对现有技术的不足，提供一种基于张量链分解的稀疏立方阵列波达方向估计方法。
5.本发明的目的是通过以下技术方案实现的：一种基于张量链分解的稀疏立方阵列波达方向估计方法，包含以下步骤：
6.(1)接收端使用m
xmymz
+n
x
nyn
z-1个物理天线阵元，按照三维嵌套互质立方阵列的结构进行架构；该嵌套互质立方阵列分解为一个均匀立方子阵列和一个稀疏立方子阵列其中包含m
x
×my
×mz
个天线阵元，阵元间距d为入射窄带信号波长λ的一半，即d＝λ/2，包含n
x
×
ny×
nz个天线阵元，在x轴方向、y轴方向和z轴方向上的间距分别为m
x
d、myd和mzd；
7.(2)假设有k个来自方向的远场窄带非相关信号源，θk和分别为第k个入射信号源的方位角和俯仰角，k＝1，2，
…
，k，将嵌套互质立方阵列中均匀立方子阵列的t个采样快拍信号在第四个维度进行叠加后，得到一个四维张量信号建模为：
[0008][0009]
其中，sk＝[s
k，1
，s
k，2
，
…
，s
k，t
]
t
为对应第k个入射信号源的多快拍采样信号波形向量，[
·
]
t
表示转置操作，表示矢量外积，为与各信号源相互独立的噪声张量，和分别为在x轴、y轴和z轴方向上的导引矢量，表示为：
[0010][0011][0012][0013]
其中，
[0014]
稀疏立方子阵列的接收信号用四维张量信号表示为：
[0015][0016]
其中，为与各信号源相互独立的噪声张量，和分别为在x轴、y轴和z轴方向上的导引矢量，表示为：
[0017][0018][0019][0020]
通过求四维张量信号和的互相关统计量，得到二阶互相关张量
[0021][0022]
其中，表示第k个入射信号源的功率，表示六维互相关噪声张量，《
·
，
·
》r表示两个张量沿着第r维度的张量缩并操作，e[
·
]表示取数学期望操作，(
·
)
*
表示共轭操作；这里，六维互相关噪声张量的第(1，1，1，1，1，1)个元素取值为噪声功率而其他元素取值均为0；
[0023]
(3)定义维度集合通过对互相关张量进行维度合并的张量变换，得到一个三维张量
[0024][0025]
其中，其中，通过在指数项上形成差集数组，构造出一个大小为m
xnx
×my
ny×mz
nz的虚拟均匀立方阵列的虚拟均匀立方阵列表示克罗内克积；噪声张量的第(1，1，1)个元素取值为而其他元素取值均为0；将三维张量中的元素进行排列，以对应中虚拟阵元的位置，得到三维虚拟域张量
[0026][0027]
其中，
[0028]
和分别为虚拟均匀立方阵列沿x轴、y轴和z轴上的导引矢量，噪声张量的第(m
xnx-m
x
+1，myn
y-my+1，mzn
z-mz+1)个元素取值为而其他元素取值均为0；
[0029]
(4)计算对应于虚拟域张量的虚拟域协方差张量的虚拟域协方差张量表示为如下canonical polyadic模型：
[0030]
[0031]
其中，对角张量的对角线元素为b
x
＝＝为canonical polyadic因子矩阵，为噪声张量，
×i表示沿着第i维度的张量-矩阵内积；
[0032]
(5)通过张量链分解，将虚拟域协方差张量表示为矩阵与三维核张量的链式连接形式，即：
[0033][0034]
其中，为首矩阵，为首矩阵，为四个三维核张量，为尾矩阵，表示沿着张量第r1维度与张量第r2维度的张量缩并操作；
[0035]
(6)对核张量进行canonical polyadic分解，表示为：
[0036][0037]
其中，对角张量的对角线元素为幅度因子[λ1，λ2，
…
，λk]，为对应虚拟域协方差张量第一和第二维度的基变换矩阵，为canonical polyadic因子矩阵by的估计值；由此，的第一个canonical polyadic因子矩阵b
x
通过张量链首矩阵g1的基变换估计得到，即：
[0038][0039]
其中(
·
)-1
表示矩阵求逆操作；从估计的canonical polyadic因子矩阵和中提取角度信息，得到二维波达方向估计结果
[0040]
进一步地，步骤(1)所述的三维嵌套互质立方阵列结构具体描述为：均匀立方子阵列中天线阵元在三维坐标系上的位置坐标为{(m
x
d，myd，mzd)|m
x
＝0，1，
…
，m
x-1，my＝0，1，
…
，m
y-1，mz＝0，1，
…
，m
z-1}；稀疏立方子阵列中天线阵元在三维坐标系上的位置坐标为{(n
xmx
d，n
ymy
d，n
zmz
d)|n
x
＝0，1，
…
，n
x-1，ny＝0，1，
…
，n
y-1，nz＝0，1，
…
，n
z-1}；{m
x
，n
x
}、{my，ny}、{mz，nz}分别为一对互质整数；将和按照坐标系原点位置阵元重叠的方式进行子阵列组合，由于和的阵元排布满足互质数条件，因此除坐标系原点位置以外的其余位置阵元均不重叠，获得实际包含m
xmymz
+n
x
nyn
z-1个天线阵元的嵌套互质立方阵列。
[0041]
进一步地，步骤(2)所述的二阶互相关张量推导，在实际中，通过计算四维张量信号和的采样互相关统计量近似得到，即采样互相关张量
[0042]
[0043]
进一步地，步骤(5)中，虚拟域协方差张量的张量链分解按如下步骤进行计算：沿着第一维度对进行展开，得到对进行截断奇异值分解：
[0044][0045]
其中，为奇异值矩阵，为右奇异值矩阵；将重构为一个辅助矩阵其中reshape{
·
}表示将一个矩阵按指定维度大小重构为另一矩阵或张量的操作；通过对c1使用截断奇异值分解得到：
[0046]
c1＝u2λ
2v2
，
[0047]
其中，为左奇异值矩阵，为奇异值矩阵，为奇异值矩阵，为右奇异值矩阵，通过重构u2得到核张量即重复上述步骤，核张量和尾矩阵g6利用左奇异值矩阵{ur，r＝3，4，5，6}依次对应生成。
[0048]
进一步地，步骤(6)中，从估计的canonical polyadic因子矩阵和中提取参数μk和vk，表示为：
[0049][0050][0051]
其中，与分别表示的第m行元素与的第n行元素，∠(
·
)表示对一个复数取幅角操作；基于(μk，vk)与二维波达方向的对应关系，得到二维波达方向估计的闭式解为：
[0052][0053][0054]
本发明的有益效果是：本发明主要解决现有方法中通过canonical polyadic(cp)分解直接处理高维虚拟域协方差张量时计算复杂度高的问题，其实现步骤是：构建三维嵌套互质立方阵列；嵌套互质立方阵列接收信号的张量建模；通过互相关张量变换构造对应于虚拟均匀立方阵列的虚拟域张量；推导虚拟域协方差张量；对虚拟域协方差张量进行张
量链分解；从张量链首矩阵和核张量中提取cp因子矩阵并获取波达方向估计结果。本发明通过张量链分解对高维虚拟域协方差张量进行降维表征，进而从低维度矩阵和核张量中高效提取角度信息，实现了低复杂度、高精度的波达方向估计，可用于目标定位。
附图说明
[0055]
图1是本发明的总体流程框图。
[0056]
图2是本发明所构建的嵌套互质立方阵列结构示意图。
[0057]
图3是本发明所提方法在不同信噪比条件下的波达方向估计精度性能比较图。
[0058]
图4是本发明所提方法在不同采样快拍数条件下的波达方向估计精度性能比较图。
[0059]
图5是本发明所提方法在不同信噪比条件下的运行时间比较图。
[0060]
图6是本发明所提方法在不同采样快拍数条件下的运行时间比较图。
具体实施方式
[0061]
以下参照附图，对本发明的技术方案作进一步的详细说明。
[0062]
为了解决现有方法通过canonical polyadic分解直接处理高维虚拟域协方差张量时计算复杂度高的问题，本发明提出了一种基于张量链分解的稀疏立方阵列波达方向估计方法，通过对虚拟域协方差张量进行张量链分解，从降维张量链首矩阵和核张量中提取canonical polyadic因子矩阵，实现高效的波达方向估计。参照图1，本发明的实现步骤如下：
[0063]
步骤1：构建三维嵌套互质立方阵列。接收端使用m
xmymz
+n
x
nyn
z-1个物理天线阵元，按照三维嵌套互质立方阵列的结构进行架构，如图2所示，该嵌套互质立方阵列分解为一个均匀立方子阵列和一个稀疏立方子阵列其中包含m
x
×my
×mz
个天线阵元，阵元间距d为入射窄带信号波长λ的一半，即d＝λ/2，中天线阵元在三维坐标系上的位置坐标为{(m
x
d，myd，mzd)|m
x
＝0，1，
…
，m
x-1，my＝0，1，
…
，m
y-1，mz＝0，1，
…
，m
z-1}；包含n
x
×
ny×
nz个天线阵元，这些天线阵元在x轴方向、y轴方向和z轴方向上的间距分别为m
x
d、myd和mzd，在三维坐标系上的位置坐标为{(n
xmx
d，n
ymy
d，n
zmz
d)|n
x
＝0，1，
…
，n
x-1，ny＝0，1，
…
，n
y-1，nz＝0，1，
…
，n
z-1}；{m
x
，n
x
}、{my，ny}、{mz，nz}分别为一对互质整数；将和按照坐标系原点位置阵元重叠的方式进行子阵列组合，由于和的阵元排布满足互质数条件，因此除坐标系原点位置以外的其余位置阵元均不重叠，获得实际包含m
xmymz
+n
x
nyn
z-1个天线阵元的嵌套互质立方阵列；
[0064]
步骤2：嵌套互质立方阵列接收信号的张量建模。假设有k个来自方向的远场窄带非相关信号源，θk和分别为第k个入射信号源的方位角和俯仰角，k＝1，2，
…
，k，将嵌套互质立方阵列中均匀立方子阵列的t个采样快拍信号在第四个维度进行叠加后，可以得到一个四维张量信号建模为：
[0065][0066]
其中，sk＝[s
k，1
，s
k，2
，
…
，s
k，t
]
t
为对应第k个入射信号源的多快拍采样信号波形向量，[
·
]
t
表示转置操作，表示矢量外积，为与各信号源相互独立的噪声张量，和分别为在x轴、y轴和z轴方向上的导引矢量，表示为：
[0067][0068][0069][0070]
其中，
[0071]
类似地，稀疏立方子阵列的接收信号可以用四维张量信号的接收信号可以用四维张量信号表示为：
[0072][0073]
其中，为与各信号源相互独立的噪声张量，和分别为在x轴、y轴和z轴方向上的导引矢量，表示为：
[0074][0075][0076][0077]
通过求四维张量信号和的互相关统计量，得到二阶互相关张量
[0078][0079]
其中，表示第k个入射信号源的功率，表示六维互相关噪声张量，《
·
，
·
》r表示两个张量沿着第r维度的张量缩并操作，e[
·
]表
示取数学期望操作，(
·
)
*
表示共轭操作。这里，六维互相关噪声张量的第(1，1，1，1，1，1)个元素取值为噪声功率而其他元素取值均为0。在实际中，通过计算张量信号和的采样互相关统计量近似得到，即采样互相关张量
[0080][0081]
步骤3：通过互相关张量变换构造对应于虚拟均匀立方阵列的虚拟域张量。由于互相关张量中包含了均匀立方子阵列和稀疏立方子阵列的空间信息，通过合并中表征同一方向空间信息的维度，可以使两个子阵列对应相同方向的导引矢量在指数项上形成差集数组，从而构造三维的虚拟均匀立方阵列。具体地，互相关张量的第1、4维度表征x轴方向的空间信息，第2、5维度表征y轴方向的空间信息，第3、6维度表征z轴方向的空间信息；为此，定义维度集合通过对互相关张量进行维度合并的张量变换，得到一个三维张量
[0082][0083]
其中，其中，通过在指数项上形成差集数组，构造出一个大小为m
xnx
×my
ny×mz
nz的虚拟均匀立方阵列的虚拟均匀立方阵列表示kronecker积。噪声张量的第(1，1，1)个元素取值为而其他元素取值均为0。将三维张量中的元素进行排列，以对应中虚拟阵元的位置，即可得到三维虚拟域张量得到三维虚拟域张量
[0084][0085]
其中，其中，和分别为虚拟均匀立方阵列沿x轴、y轴和z轴上的导引矢量，噪声张量的第(m
xnx-m
x
+1，myn
y-my+1，mzn
z-mz+1)个元素取值为而其他元素取值均为0；
[0086]
步骤4：推导虚拟域协方差张量。计算对应于虚拟域张量的虚拟域协方差张量表示为如下canonicalpolyadic模型：
[0087][0088]
其中，对角张量的对角线元素为的对角线元素为的对角线元素为的对角线元素为为canonical polyadic因子矩阵，为噪声张量，
×i表示沿着第i维度的张量-矩阵内积；
[0089]
步骤5：对虚拟域协方差张量进行张量链分解。为了确保对虚拟域协方差张量中的canonical polyadic因子进行高效提取，不同于传统方法对其直接进行canonical polyadic分解的操作，通过张量链分解对其进行降维处理。具体而言，将虚拟域协方差张量表示为矩阵与三维核张量的链式连接形式，即：
[0090][0091]
其中，为首矩阵，为首矩阵，为四个三维核张量，为尾矩阵，表示沿着张量第r1维度与张量第r2维度的张量缩并操作。上述虚拟域协方差张量的张量链分解可按如下步骤进行计算：沿着第一维度对进行展开，得到进而对进行截断奇异值分解(truncated singular value decomposition)：
[0092][0093]
其中，为奇异值矩阵，为右奇异值矩阵。将重构为一个辅助矩阵其中reshape{
·
}表示将一个矩阵按指定维度大小重构为另一矩阵或张量的操作。通过对c1使用截断奇异值分解可得到：
[0094]
c1＝u2λ
2v2
，
[0095]
其中，为左奇异值矩阵，为奇异值矩阵，为奇异值矩阵，为右奇异值矩阵，通过重构u2可得到核张量即重复上述步骤，核张量和尾矩阵g6可利用左奇异值矩阵{ur，r＝3，4，5，6}依次对应生成；
[0096]
步骤6：从张量链首矩阵和核张量中提取cp因子矩阵并获取波达方向估计结果。由于canonical polyadic因子矩阵b
x
与by可用于估计信号源的二维波达方向，因此将虚拟域协方差张量通过张量链分解得到的首矩阵g1和核张量转换为b
x
与by。具体过程为：对核张量进行canonical polyadic分解，表示为：
[0097][0098]
其中，对角张量的对角线元素为幅度因子[λ1，λ2，
…
，λk]，为对应虚拟域协方差张量第一和第二维度的基变换矩阵，为canonical polyadic因子矩阵by的估计值。由此，的第一个canonical polyadic因子矩阵b
x
可以通过张量链首矩阵g1的基变换得到，即：
[0099][0100]
其中(
·
)-1
表示矩阵求逆操作。进而，从和中估计μk和νk，表示为：
[0101][0102][0103]
其中，与分别表示的第m行元素与的第n行元素，∠(
·
)表示对一个复数取幅角操作。基于(μk，vk)与二维波达方向的对应关系，得到二维波达方向估计的闭式解为：
[0104][0105][0106]
下面结合仿真实例对本发明的效果做进一步的描述。
[0107]
仿真实例：采用嵌套互质立方阵列接收入射信号，其参数选取为m
x
＝my＝mz＝3，n
x
＝ny＝nz＝4，即架构的嵌套互质立方阵列共包含m
xmymz
+n
x
nyn
z-1＝90个物理阵元。所构造的虚拟均匀立方阵列共有12
×
12
×
12个虚拟阵元。将本发明所提基于张量链分解的稀疏立方阵列波达方向估计方法与传统基于canonical polyadic分解的方法和基于旋转不变子空间(estimating signal parameter via rotational invariance techniques，esprit)的方法进行对比。假定有2个入射信号，每次试验的信号方位角和俯仰角均在[15
°
，60
°
]范围内随机生成。在采样快拍数t＝150的条件下，绘制均方根误差(root-mean-square error，rmse)随信噪比(signal-to-noise ratio，snr)变化的性能对比曲线，如图3所示；在信噪比snr＝-5db的条件下，绘制rmse随采样快拍数t变化的性能对比曲线，如图4所示。进一步地，对比上述方法的计算效率，绘制计算机运行时间随snr变化的对比曲线，如图5所示；并绘制计算机运行时间随采样快拍数t的对比曲线，如图6所示。从对比结果看出，相比于传统基于esprit的方法，本发明所提方法通过构建虚拟域张量，充分利用所对应高维虚拟域协方差
张量的结构化信息，具备更优的波达方向估计性能，且基于张量链分解对虚拟域协方差张量进行高效特征提取，具备更短的运行时间；相比基于canonical polyadic分解的方法，本发明所提方法在虚拟域协方差张量的张量链分解过程中最大程度降低信息损失，不仅表现出更高的估计精度，还使得计算复杂度大幅度降低，显著提升了计算效率。
[0108]
综上所述，本发明通过对高维虚拟域协方差张量进行张量链分解，并设计了相应canonical polyadic因子的高效提取策略，实现了低复杂度、高精度的稀疏立方阵列波达方向估计。
[0109]
以上所述仅是本发明的优选实施方式，虽然本发明已以较佳实施例披露如上，然而并非用以限定本发明。任何熟悉本领域的技术人员，在不脱离本发明技术方案范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案做出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何的简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。