一种互质线阵中基于级联的子空间拟合角度估计方法
技术领域
本发明涉及通信
技术领域
，尤其涉及一种互质线阵中基于级联的子空间拟合角度估计方法。背景技术
在传统的来向角(Direction of Arrival,DOA)估计算法中，阵列中相邻阵元的间距必须不大于入射信号半波长，以免产生角度模糊。但此时存在相邻阵元互耦严重、阵列孔径小、分辨率差和测向精度低等不足，且无法满足实际环境需求。在这种背景下阵元间距大于半波长的稀疏阵列被提出并开始推广应用。互质阵列是稀疏阵列的一种主要形式，通过子阵的互质关系能够有效消除测向模糊，与阵元数相等的均匀线阵相比具有更大的阵列孔径，以及更高的自由度，从而在谱估计精度和分辨率上更有优势。
由于互质阵列的设计要求阵元间距是半波长的整数倍，而当阵元间距大于半波长时会有模糊现象产生。子空间拟合算法因其具有较高的估计精度被广泛的应用于阵列信号处理和通信工程领域，Viberg M等提出了均匀阵列中多信源下的子空间拟合算法，包括信号子空间拟合算法和噪声子空间拟合算法，该算法能实现近似于最大似然的DOA估计性能。
发明内容
为了解决现有技术存在的上述问题，本发明的目的在于提供一种互质线阵中基于级联的子空间拟合角度估计方法，解决了如何将均匀阵列空间谱算法应用于互质阵的角度估计中，突破了传统均匀阵阵元间距不大于半波长的限制，充分发挥互质阵列的优势，在提升阵列参数估计性能的同时贴近实际场景，能够从阵列接收信号中获得更多关于所接收信号的信息，使其实际应用性更强。
为了解决上述技术问题，本发明采用的技术方案如下：一种互质线阵中基于级联的子空间拟合角度估计方法，包括以下步骤：
步骤1：建立互质阵列信号的数学模型；
步骤2：根据建立的互质阵列信号的数学模型求互质阵列信号的协方差矩阵；
步骤3：对求得的协方差矩阵进行特征值分解，得到信号子空间；
步骤4：根据建立的阵列信号的数学模型中的传播算子得到互质阵列各子阵的模糊角度估计值；
步骤5：根据互质阵列各子阵的模糊角度估计值恢复角度的所有估计值，包括真实值和模糊值，再通过互质线阵两个子阵列的互质特性消除所有角度估计值中的模糊角度，得到真实的初始角度估计值；
步骤6：根据得到的阵列的初始角度估计值和信号子空间，重构第1个方向矩阵且在设定的搜索区间根据拟合函数进行局部搜索，得到第1个入射信号角度θ1的估计值
步骤7：根据得到的第1个入射信号角度的估计值重构第2个方向矩阵且在设定的搜索区间根据拟合函数进行局部搜索，得到第2个入射信号角度θ2的估计值
步骤8：依次类推，直到根据得到的K-1个入射信号角度的估计值，重构第K个方向矩阵再通过设定的搜索区间进行局部搜索，得到第K个入射信号角度θK的估计值
作为本发明的进一步改进，所述步骤1具体包括：
所述互质线阵结构由两个均匀线阵组成，包括子阵1和子阵2，所述子阵1含有M个阵元，子阵2含有N个阵元，其中M和N互为质数，且M＜N，由于两个子阵在(0,0)点重合，所以整个阵列的阵元数目表示为T＝M+N-1，子阵1的阵元间距是d1＝Nλ/2，子阵2的阵元间距为d2＝Mλ/2,其中λ表示波长；
设信源数为K，假设有K个非相关的窄带远场信号入射到该阵列中，其中 K＜min{M,N}，入射方向Θ表示为Θ＝[θ1,θ2,…,θK]，
其中，θ1表示第1个入射信号角度，θ2表示第2个入射信号角度，θK表示第K个入射信号角度，θ1,θ2,…,θK统一记作θk，表示第k个入射信号角度，且θk∈[0,π/2]，k＝1,2,…,K，子阵i的信号模型可表示为，
xi(t)＝Ais(t)+ni(t) 公式1
式中，i＝1,2，s(t)为窄带远场非相干信号，具体表示为s(t)＝[s1(t),s2(t),…,sK(t)]T，其中，t＝1,2,…,L，L表示快拍数，ni(t)为子阵i均值为0方差为的加性、白高斯噪声矢量，且独立于信源，Ai为第i个子阵的方向矩阵，具体表示为Ai＝[ai(θ1),ai(θ2),…,ai(θK)]，ai(θ1)为子阵i第1个入射信号的方向矢量，ai(θ2)为子阵i第2个入射信号的方向矢量，ai(θK)为子阵i第K个入射信号的方向矢量，ai(θ1),ai(θ2),…,ai(θK)统一记为ai(θk)，表示子阵i第k个信号的方向矢量，其中k＝1,2,…,K，分别得到子阵1的方向矢量a1(θk)和子阵2的方向矢量 a2(θk)，具体表示为
式中，j为复数虚部，(·)T为矩阵的转置操作；
通过将两个子阵的接收信号堆叠，得到整个阵列的总接收信号，
其中，A1为第1个子阵的方向矩阵，A2为第2个子阵的方向矩阵， n(t)表示均值为0，方差为的加性、白高斯噪声矢量，具体表示为s(t) 为入射信源。
作为本发明的进一步改进，所述步骤2具体包括：
所述阵列信号的协方差矩阵具体表示为，
式中，x为公式4中的信号矩阵，L为快拍数，xH为信号矩阵x的伪逆操作。
作为本发明的进一步改进，所述步骤3具体包括：
对协方差矩阵进行特征值分解，可得
式中，Us为由对应K个最大特征值的特征矢量构成的信号子空间，Ds是由K个最大特征值构成的对角矩阵，Un为其余特征矢量构成的噪声子空间，Dn是由其它所有特征值构成的对角矩阵，(·)H为矩阵的伪逆操作。
作为本发明的进一步改进，所述步骤4具体包括如下步骤：
步骤4-1，根据建立的互质信号的数学模型对子阵1的方向矩阵进行分块，得到方向矩阵A1的前K行和(M-K)行矩阵，分别表示为A11和A12，假设A11是K×K的非奇异矩阵，则 A12可以通过A11的线性变化得到，即A12＝P1cA11，其中P1c为子阵1的传播算子，为维度，定义中间矩阵P1为
式中，I1K是一个K×K的单位矩阵，得到
步骤4-2，再分别用P1a和P1b表示P1的前M-1行和后M-1行，用A1a和A1b表示A1的前 M-1行和后M-1行，可得
则存在关系式定义由于A11为满秩矩阵，Ψ1r为Φ1r的相似变换，Φ1r为包括特征值的对角矩阵，因此Ψ1r和Φ1r有相同的特征值，对Ψ1r进行特征值分解即可得到Φ1r的对角线元素δ1,k，从而得到基于子阵1的接收信号的模糊角度估计值用正弦函数具体表示为其中，angle(·)表示相位函数，k＝1,2,…,K；
步骤4-3，相同的方式，得到基于子阵2的接收信号的模糊角度估计值用正弦函数具体表示为，其中，k＝1,2,…,K。
作为本发明的进一步改进，所述步骤5具体包括：
根据正弦函数的特性可知真实值和模糊值之间的相位差是即
其中，ki表示整数，i＝1,2，
k1＝-(N-1),…,-1,1,…,(N-1)，k2＝-(M-1),…,-1,1,…,(M-1)，
d1＝Nλ/2，d2＝Mλ/2,其中λ表示波长，
对于每一个子阵，分别存在
步骤5-1，根据互质阵列各子阵的模糊角度估计值恢复角度的所有估计值；
式中，为子阵1的K个信号的角度真实值和模糊值的集合，具体表示为
为子阵2的K个信号的角度真实值和模糊值的集合，具体表示为
步骤5-2，根据互质阵列的互质特性消除所有角度估计值中的模糊角度，得到阵列真实的初始角度估计值，
在两个备选集合和中选择距离最近的角度作为估计值，
式中，为子阵2所有的估计值的集合，Ω＝1,2,…,N，为子阵1所有的估计值的集合，
最后，得到阵列的真实的初始角度估计值具体表示为，
式中，k＝1,2,…,K。
作为本发明的进一步改进，所述步骤6具体包括如下步骤：
步骤6-1，由阵列的真实的初始角度估计值得到初始角度估计值的集合具体表示为得到第1个入射信号的方向矩阵具体表示为，
式中，a(θ)为待测角度的方向矢量，为利用第2个初始估计角度的重构的方向矢量，为第K个初始估计角度的重构的方向矢量，统一记作表示第k个初始估计角度的重构的方向矢量，分别表示为，
其中，k＝2,3,…,K；
步骤6-2，根据第1个方向矩阵和信号子空间，且设定第一搜索区间进行局部搜索，得到第1个入射信号角度的估计值
其中，第一搜索区间为Δθ为设定系数，为Frobenius范数，为信号子空间。
作为本发明的进一步改进，所述步骤7具体包括如下步骤：
步骤7-1，根据第1个入射信号角度的估计值且保持不变，得到第2个入射信号的方向矩阵
式中，为第1个初始估计角度的重构的方向矢量，为第3个初始估计角度的重构的方向矢量，为第K个初始估计角度的重构的方向矢量，统一记作表示第k个初始估计角度的重构的方向矢量，其中，k＝3,4,…,K；
步骤7-2，根据第2个方向矩阵且设定第二搜索区间进行局部搜索，得到第2个入射信号角度θ2的估计值为，
其中，第二搜索区间为
作为本发明的进一步改进，所述步骤8具体包括如下步骤：
步骤8-1，由此类推，根据第K-1个入射信号角度的估计值得到第K个方向矩阵
步骤8-2，再根据第K个方向矩阵且设定第三搜索区间进行局部搜索，得到第K个入射信号角度θK的估计值
其中，第三搜索区间为
本发明与现有技术相比，具有如下有益效果：
1、本发明一种互质线阵中基于级联的子空间拟合角度估计方法，通过子阵的互质关系能够有效消除测向模糊，相较于传统均匀阵列，能够获得更好的角度估计性能。
2、本发明一种互质线阵中基于级联的子空间拟合角度估计方法，能突破传统均匀阵列阵元间距不大于半波长的限制，拓展阵列孔径。
附图说明
图1是互质线阵结构示意图；
图2是本发明算法(LC-SF)估计得到角度参数的散点图；
图3是分别采用本发明算法(LC-SF)与传统的信号子空间算法(S-SF)、噪声子空间拟合算法(N-SF)和全谱峰搜索算法(TSS)在不同阵元数下的复杂度对比图；
图4是分别采用本发明算法(LC-SF)与传统的信号子空间算法(S-SF)、噪声子空间拟合算法(N-SF)和全谱峰搜索算法(TSS)在不同快拍数下的复杂度对比图；
图5是分别采用本发明算法(LC-SF)与传统的信号子空间算法(S-SF)、噪声子空间拟合算法(N-SF)和全谱峰搜索算法(TSS)在不同搜索步长下的复杂度对比图；
图6是本发明互质阵的(LC-SF)算法与均匀阵的(LC-SF)算法角度估计性能在不同快拍数下的对比图；
图7是本发明算法(LC-SF)的角度估计性能在不同信噪比下的对比图；
图8是本发明算法(LC-SF)的角度估计性能在不同快拍数下的对比图；
图9是本发明算法(LC-SF)实施方法流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明涉及一种互质线阵中基于级联的子空间拟合角度估计方法，本发明的方法运算复杂度分析如下：
本小节对所设计的LC-SF算法、信号子空间算法S-SF，噪声子空间拟合算法N-SF以及全谱峰搜索算法TSS算法的复杂度进行了比较。
对于LC-SF算法，其计算复杂度表示为
O((M2+N2)L+M3+N3+2[(M-1)2+(N-1)2]K+(M-1)3+(N-1)3+[(M-1)+(N-1)]K2]+ G1[4K(M2+N2)+(M3+N3)])，其中G1为级联机制的搜索次数，G1＝K·2Δ/ds，ds＝0.001表示搜索步长，O表示计算复杂度；
对于S-SF算法，复杂度表示为O((M2+N2)L+M3+N3+G2[4K(M2+N2)+(M3+N3)])，其中G2为全局搜索次数，G2＝(90/ds)K；
对于N-SF算法，复杂度表示为O((M2+N2)L+M3+N3+G2[M(M-K)K+N(N-K)K])；
对于TSS算法，复杂度表示为O((M2+N2)L+M3+N3+G3[M(M-K)+N(N-K)])，其中， G3＝90/ds。
表1-1给出了不同算法在CLA阵列下的计算复杂度比较，表1-2给出了不同算法在一次迭代下的运行时间，该测试是在Intel(R)Xeon(R)CPU [email protected]的MATLAB R2015b环境下进行。
仿真参数设置如下：M＝3,N＝4，K＝2，L＝100，ds＝0.01，SNR＝10dB。从表1-2 中可以看出所提算法LC-SF比其他算法拥有更低的计算复杂度。
表1-1不同算法的计算复杂度比较
表1-2不同算法的运算时间比较
图2为当2个信号入射到互质线阵时，利用本发明的估计方法得到的角度估计的散点图，由图可以看出本发明的算法可以有效地估计出角度。其中，入射信号的角度为θ1＝30°,θ2＝50°,互质线阵的两个子阵阵元数分别为M＝5，N＝6，快拍数L＝200，信噪比 SNR＝-5dB。
图3是本发明的LC-SF算法在不同阵元数下的复杂度比较图。由图可以得出，随着阵元数N的增多，复杂度增加。相较于传统的S-SF算法、N-SF算法和TSS算法，本发明的算法能够大大降低计算复杂度。其中，快拍数L＝200，步长ds＝0.01。
图4为本发明的LC-SF算法在不同快拍数下的复杂度比较图。相较于传统的S-SF算法、 N-SF算法和TSS算法，本发明的LC-SF算法能够实现更低的计算复杂度。其中，互质线阵的两个子阵阵元数分别为M＝5，N＝6，步长ds＝0.01。从表1-1可以看出快拍数L仿真设置时最大为500，与快拍数L相乘的部分最多为十的三次方数量级，相对于搜索步长ds部分的12次方数量级，微乎其微，从图4中看到随着快拍数L的增加，各种算法的复杂度的变化很微小，所以图中几乎是平的。
图5为本发明的LC-SF算法在不同步长下的复杂度比较图。由图可以得出，随着搜索步长ds的减小，复杂度增加。相较于传统的S-SF算法、N-SF算法和TSS算法，本发明的 LC-SF算法能够大大降低计算复杂度。其中，互质线阵的两个子阵阵元数分别为M＝5，N＝6，快拍数L＝200。
图6是本发明互质阵的LC-SF算法(对应图中的LC-SF for CLA)与均匀阵的LC-SF算法 (对应图中的LC-SF for ULA)角度估计性能在相同算法不同阵列且不同快拍数L下的比较图。随着快拍数L增加，即采样数据增多。由图可以得出，互质阵中的LC-SF算法的角度估计性能明显好于均匀阵中的LC-SF算法，其中，入射信号的角度分别为θ1＝25°,θ2＝35°,互质线阵的两个子阵阵元数分别为M＝5，N＝6，信噪比SNR＝5dB。
图7是本发明的LC-SF算法角度估计性能在不同信噪比下的比较图。由图可以得出，算法的角度估计性能随着信噪比SNR的增加都变得更好，同时可以看到，本发明的LC-SF算法在估计性能上近似于传统的S-SF算法和N-SF算法，且优于ESPRIT算法和PM算法。其中，入射信号的角度分别为θ1＝25°,θ2＝35°,互质线阵的两个子阵阵元数分别为 M＝5，N＝6，快拍数L＝200。
图8是本发明的LC-SF算法角度估计性能在不同快拍数下的比较图。由图可以得出，算法的角度估计性能随着块拍数的增加都变得更好，同时可以看到，本发明的LC-SF算法在估计性能上近似于传统的S-SF算法和N-SF算法，且优于ESPRIT算法和PM算法。其中，入射信号的角度分别为θ1＝25°,θ2＝35°,互质线阵的两个子阵阵元数分别为M＝5，N＝6，快拍数L＝200。
本发明涉及一种互质线阵中基于级联的子空间拟合角度估计方法，包括以下步骤：
步骤1：建立互质阵列信号的数学模型；
所述互质线阵结构由两个均匀线阵组成，包括子阵1和子阵2，所述子阵1含有M个阵元，子阵2含有N个阵元，其中M和N互为质数，且M＜N，由于两个子阵在(0,0)点重合，所以整个阵列的阵元数目表示为T＝M+N-1，子阵1的阵元间距是d1＝Nλ/2，子阵2的阵元间距为d2＝Mλ/2,其中λ表示波长；
设信源数为K，假设有K个非相关的窄带远场信号入射到该阵列中，其中 K＜min{M,N}，入射方向Θ表示为Θ＝[θ1,θ2,…,θK]，
其中，θ1表示第1个入射信号角度，θ2表示第2个入射信号角度，θK表示第K个入射信号角度，θ1,θ2,…,θK统一记作θk，表示第k个入射信号角度，且θk∈[0,π/2]，k＝1,2,…,K，子阵i的信号模型可表示为，
xi(t)＝Ais(t)+ni(t) 公式1
式中，i＝1,2，s(t)为窄带远场非相干信号，具体表示为s(t)＝[s1(t),s2(t),…,sK(t)]T，其中， t＝1,2,…,L，L表示快拍数，ni(t)为子阵i均值为0方差为的加性、白高斯噪声矢量，且独立于信源，Ai为第i个子阵的方向矩阵，具体表示为Ai＝[ai(θ1),ai(θ2),…,ai(θK)]，ai(θ1)为子阵i第1个入射信号的方向矢量，ai(θ2)为子阵i第2个入射信号的方向矢量，ai(θK)为子阵i第K个入射信号的方向矢量，ai(θ1),ai(θ2),…,ai(θK)统一记为ai(θk)，表示子阵i第k个信号的方向矢量，其中k＝1,2,…,K，分别得到子阵1的方向矢量a1(θk)和子阵2的方向矢量 a2(θk)，具体表示为
式中，j为复数虚部，(·)T为矩阵的转置操作；
通过将两个子阵的接收信号堆叠，得到整个阵列的总接收信号，
其中，A1为第1个子阵的方向矩阵，A2为第2个子阵的方向矩阵， n(t)表示均值为0，方差为的加性、白高斯噪声矢量，具体表示为s(t) 为入射信源。
步骤2：根据建立的互质阵列信号的数学模型求互质阵列信号的协方差矩阵；
阵列信号的协方差矩阵具体表示为，
式中，x为公式4中的信号矩阵，L为快拍数，xH为信号矩阵x的伪逆操作。
步骤3：对求得的协方差矩阵进行特征值分解，得到信号子空间；
对协方差矩阵进行特征值分解，可得
式中，Us为由对应K个最大特征值的特征矢量构成的信号子空间，Ds是由K个最大特征值构成的对角矩阵，Un为其余特征矢量构成的噪声子空间，Dn是由其它所有特征值构成的对角矩阵，(·)H为矩阵的伪逆操作。
步骤4：根据建立的阵列信号的数学模型中的传播算子得到互质阵列各子阵的模糊角度估计值；
步骤4-1，根据建立的互质信号的数学模型对子阵1的方向矩阵进行分块，得到方向矩阵A1的前K行和(M-K)行矩阵，分别表示为A11和A12，假设A11是K×K的非奇异矩阵，则 A12可以通过A11的线性变化得到，即A12＝P1cA11，其中P1c为子阵1的传播算子，为维度，定义中间矩阵P1为
式中，I1K是一个K×K的单位矩阵，得到
步骤4-2，再分别用P1a和P1b表示P1的前M-1行和后M-1行，用A1a和A1b表示A1的前 M-1行和后M-1行，可得
则存在关系式定义由于A11为满秩矩阵，Ψ1r为Φ1r的相似变换，Φ1r为包括特征值的对角矩阵，因此Ψ1r和Φ1r有相同的特征值，对Ψ1r进行特征值分解即可得到Φ1r的对角线元素δ1,k，从而得到基于子阵1的接收信号的模糊角度估计值用正弦函数具体表示为其中，angle(·)表示相位函数，k＝1,2,…,K；
步骤4-3，相同的方式，得到基于子阵2的接收信号的模糊角度估计值用正弦函数具体表示为，其中，k＝1,2,…,K。
步骤5：根据互质阵列各子阵的模糊角度估计值恢复角度的所有估计值，包括真实值和模糊值，再通过互质线阵两个子阵列的互质特性消除所有角度估计值中的模糊角度，得到真实的初始角度估计值；
根据正弦函数的特性可知真实值和模糊值之间的相位差是即
其中，ki表示整数，i＝1,2，
k1＝-(N-1),…,-1,1,…,(N-1)，k2＝-(M-1),…,-1,1,…,(M-1)，
d1＝Nλ/2，d2＝Mλ/2,其中λ表示波长，
对于每一个子阵，分别存在
对于传统的ULA结构而言，由于阵元间距小于半波长，k1和k2只能为0，无模糊现象产生，然而，由于阵元间距的增大，存在多对k1和k2使得公式11和公式12成立；
步骤5-1，根据互质阵列各子阵的模糊角度估计值恢复角度的所有估计值；
式中，为子阵1的K个信号的角度真实值和模糊值的集合，具体表示为
为子阵2的K个信号的角度真实值和模糊值的集合，具体表示为
步骤5-2，根据互质阵列的互质特性消除所有角度估计值中的模糊角度，得到阵列真实的初始角度估计值，
由于M和N为互质整数对，因此唯一存在k1＝k2＝0使得公式15成立，意味着可以通过两个子阵的估计值得到阵列的真实的初始角度估计值，实际中，由于噪声的存在，不可能直接找到无误差的完全重合的估计，因此，在两个备选集合和中选择距离最近的角度作为估计值，
式中，为子阵2所有的估计值的集合，Ω＝1,2,…,N，为子阵1所有的估计值的集合，
最后，得到阵列的真实的初始角度估计值具体表示为，
式中，k＝1,2,…,K。
步骤6：根据得到的阵列的初始角度估计值和信号子空间，重构第1个方向矩阵且在设定的搜索区间根据拟合函数进行局部搜索，得到第1个入射信号角度θ1的估计值
步骤6-1，由阵列的真实的初始角度估计值得到初始角度估计值的集合具体表示为得到第1个入射信号的方向矩阵具体表示为，
式中，a(θ)为待测角度的方向矢量，为利用第2个初始估计角度的重构的方向矢量，为第K个初始估计角度的重构的方向矢量，统一记作表示第k个初始估计角度的重构的方向矢量，分别表示为，
其中，k＝2,3,…,K；
步骤6-2，根据第1个方向矩阵和信号子空间，且设定第一搜索区间进行局部搜索，得到第1个入射信号角度的估计值
其中，第一搜索区间为Δθ为设定系数，为Frobenius范数，为信号子空间，Δθ是一个特别小的值。
步骤7：根据得到的第1个入射信号角度的估计值重构第2个方向矩阵且在设定的搜索区间根据拟合函数进行局部搜索，得到第2个入射信号角度θ2的估计值
步骤7-1，根据第1个入射信号角度的估计值且保持不变，得到第2个入射信号的方向矩阵
式中，为第1个初始估计角度的重构的方向矢量，为第3个初始估计角度的重构的方向矢量，为第K个初始估计角度的重构的方向矢量，统一记作表示第k个初始估计角度的重构的方向矢量，其中，k＝3,4,…,K；
步骤7-2，根据第2个方向矩阵且设定第二搜索区间进行局部搜索，得到第2个入射信号角度θ2的估计值为，
其中，第二搜索区间为
步骤8：依次类推，直到根据得到的K-1个入射信号角度的估计值，重构第K个方向矩阵再通过设定的搜索区间进行局部搜索，得到第K个入射信号角度θK的估计值
步骤8-1，由此类推，根据第K-1个入射信号角度的估计值得到第K个方向矩阵
步骤8-2，再根据第K个方向矩阵且设定第三搜索区间进行局部搜索，得到第K 个入射信号角度θK的估计值
其中，第三搜索区间为
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。
- 上一篇：石墨接头机器人自动装卡簧、装栓机
- 下一篇：一种基于嵌套阵的循环平稳信号直接定位方法
