一种非相干分布宽带源的角度和频率参数估计方法
技术领域
本发明涉及信号处理
技术领域
，具体涉及一种非相干分布宽带源的角度和频率参数估计方法。背景技术
宽带信号具有分辨率高、信息量大、抗干扰能力强等特点，已经广泛应用于雷达、声呐、无线通信等领域。大多数宽带源的波达角(DOA)估计是基于点源假设。然而，源信号多径或散射传播引起的角分布是不可忽略的。大部分分布源的DOA估计方法是基于窄带假设。对于分布宽带源来说，随着频率带宽增大，基于窄带假设方法的DOA估计偏差也会增大。因此，分布宽带源的DOA估计问题需要更多的关注和研究。
近年来，一些分布宽带源的DOA估计方法已经被提出来。最大似然(ML)方法，通过重建角度分布的形状来估计角度参数。协方差匹配(CM)方法，一种基于频域快拍的匹配方法。然而这两种方法由于多维搜索维度随源数增加，使得计算负担较大。参数多项式方法，一种基于信号协方差时延的自回归模型方法。这三种方法都假定分布宽带源的宽带信号的频率分布已知，且频率范围相同。而实际应用中，宽带信号可能拥有不同的频率范围。
一些研究者也提出了一种将分数阶傅立叶变换(FRFT)算法与分布源参数估计器(DSPE)算法相结合的新方法，即FRFT-DSPE方法(参考文献Yu J,Zhang L,LiuK.Coherently Distributed Wideband LFM Source Localization[J].IEEE SignalProcessing Letters,2014,22(4):504-508.)。虽然该方法无需已知频率分布，但与上述方法一样，需要已知角度分布的参数化模型，而实际应用中确切的角度分布模型往往是未知的；同时该方法能估计源的角度分布，然而并不能有效估计频率分布，事实上频率也是目标源重要的特征参数。也有学者将稀疏贝叶斯学习(SBL)方法应用于分布宽带源，虽然该方法不需要已知角度和频率的分布模型，然而它要求分布宽带源在空间域中必须是稀疏的。而实际上分布源的稀疏性会随角度扩展而降低。
现有分布宽带源的角度分布估计方法，存在需要已知角度的参数化模型和频率分布以及计算复杂度高的问题，难以满足复杂场景下非相干分布宽带源的角度分布估计。因此提供一种无需已知角度和频率分布模型，又能降低计算复杂度，并能同时估计宽带分布源的角度和频率分布的方法是十分必要的。
发明内容
本发明的目的是针对现有分布宽带源的角度分布估计方法，存在需要已知角度的参数化模型和频率分布以及计算复杂度高的问题，基于低秩矩阵恢复公开发明了一种非相干分布宽带源的角度和频率参数估计方法。
本发明的目的至少通过如下技术方案之一实现。
一种非相干分布宽带源的角度和频率参数估计方法，包括如下步骤：
S1、通过接收阵列获取远场非相干分布宽带源信号，估计非相干分布宽带源信号的协方差并将其向量化；
S2、建立向量化的非相干分布宽带源信号协方差模型，根据协方差模型中的方向导数矩阵的行重复性和奇异性，同时对方向导数矩阵和步骤S1中估计的向量化的非相干分布宽带源信号协方差进行降维，以减少后续计算过程的复杂度；
S3、将步骤S2中向量化的非相干分布宽带源信号协方差模型公式化为秩最小化问题，结合交替方向乘子法(ADMM)和迭代加权核范数(IRNN)算法对秩最小化问题求解，得到角度-频率联合分布矩阵的估计量；
S4、根据步骤S3中得到的角度-频率联合分布矩阵的估计量，结合离散网格估计器估计出非相干分布宽带源的角度和频率分布的关键参数。
进一步地，步骤S1中，首先设置一个包括L个阵元的接收阵列，阵元间距为源信号最大频率对应的半波长；假设有K个远场非相干分布宽带源信号入射到接收阵列，则接收阵列在t时刻的输出信号为：
式中t＝1,2,...,Q，这里Q为信号采样的快拍数；f和Г分别是源信号的频率和频率范围；θ和Θ分别是源信号在空间传播的方位角和角度范围；dSk(f)是第k个源信号频率谱的一个测量；γk(θ,t)是第k个源信号在角度方向传播的复增益；是接收阵列关于(θ,f)的方向导数，j为虚数单位，(·)T为转置，τl(θ)为源信号传播到第l个阵元相对于传播到参考阵元(通常设置第一个阵元为参考阵元)的时延差，l＝1,2,...,L；n(t)为高斯白噪声；
根据输出信号x(t)，得到接收阵列输出信号的协方差为：
式中E(·)表示期望；(·)H表示共轭转置；Rs与分别是源信号与高斯白噪声的协方差；R的估计量为 为高斯白噪声方差，I为单位向量，对进行特征值分解，可将估计量设定为的最小特征值；由公式(2)可得到非相干分布宽带源信号的协方差估计量为：
向量化的为：
式中vec(·)表示将矩阵逐列堆叠为列向量，这里为复数集合。
进一步地，步骤S2的具体步骤如下：
S2.1、根据步骤S1中接收阵列输出信号x(t)，建立向量化的非相干分布宽带源信号协方差模型rs，并进一步构造方向导数矩阵
S2.2、利用方向导数矩阵的行重复性和奇异性，同时对方向导数矩阵和步骤S1中估计的向量化的非相干分布宽带源信号协方差进行降维。
4、根据权利要求3所述的一种非相干分布宽带源的角度和频率参数估计方法，其特征在于，步骤S2.1中，首先，对于非相干分布宽带源，假设不同源的信号不相关，同一源中来自不同角度的复增益不相干，同时源信号中不同频率的信号也不相干；由此，根据x(t)模型，可得出非相干分布宽带源信号的协方差为：
式中aθ,f为方向导数a(θ,f)的简化形式； 为第k个非相干分布宽带源信号传播复增益的角度分布，pk(f)为第k个非相干分布宽带源信号的频率分布；为K个非相干分布宽带源信号的角度-频率联合分布；
其次，将非相干分布宽带源信号的协方差公式(5)的积分式用求和式代替；用Pk表示当θ＝θ1,θ2,...,θM和f＝f1,f2,...,fN时pk(θ,f)的离散化，M和N分别是角度范围Θ和频率范围Γ离散化的数量，其中Pk(θm,fn)＝pk(θm,fn)表示矩阵Pk的第(m,n)个元素的值为pk(θm,fn)；同时用P表示pθ,f的离散化，这里为实数集合，那么P表示K个非相干分布宽带源的角度-频率联合分布矩阵；P是本发明待估计的未知量，P的求解是估计角度和频率参数的关键；
根据离散化，非相干分布宽带源信号的协方差公式(5)可改写为：
式中分别是Aθ，f和pθ，f在离散点(θm，fn)的值；
向量化的Rs可表示为：
式中是的向量化；pθ，f＝vec(P)，为矩阵P的向量化，具体如下：
最终，根据向量化的Rs模型公式(7)，可得出要构造的方向导数矩阵 这里为复数集合，具体如下：
进一步地，步骤S2.2中，首先，由于在离散点(θm，fn)处的方向导数a(θm,fn)向量的元素成等比数列，使得中存在相同的元素，那么方向导数矩阵存在相同的行；因此可删除方向导数矩阵中重复的行，同时可删除估计的向量化的非相干分布宽带源信号协方差中对应的行，以降低和的维度；
其次，方向导数矩阵是奇异的，可对和进行二次降维；其过程包括：
S2.2.1、对进行奇异值分解(SVD)，即∑A、UA和UA分别是进行SVD的奇异值对角矩阵、左奇异矩阵和右奇异矩阵；
S2.2.2、保留奇异值主成分；设置一个接近于1的奇异值主成分比例ε，当奇异值对角矩阵∑A对角线上前(q+1)个奇异值之和与所有奇异值之和的比例首次大于ε时，保留前q个较大的奇异值，删除其余较小的奇异值，即主成分的奇异值对角矩阵变为主成分对应的左奇异矩阵和右奇异矩阵更新为
S2.2.3、对方向导数矩阵和向量化的非相干分布宽带源信号协方差进行二次降维，令令
进一步地，步骤S3的具体步骤如下：
S3.1、根据步骤S2中向量化的非相干分布宽带源信号协方差模型Rs的角度-频率联合分布矩阵P具有低秩属性，将向量化的非相干分布宽带源信号协方差模型公式(7)化为秩最小化问题；
S3.2、结合交替方向乘子法(ADMM)和迭代加权核范数(IRNN)算法对秩最小化问题求解，估计角度-频率联合分布矩阵P。
进一步地，步骤S3.1中，首先，由于角度与频率是不相关的两个量，待估计量角度-频率联合分布矩阵P可表示为其中，pk(θ)＝[pk(θ1)，...，pk(θM)]T，pk(f)＝[pk(f1)，...,pk(fN)]T分别为角度和频率分布的离散向量，
由于向量的秩为1，可知矩阵Pk的秩也为1，从而P的秩小于等于K；而源数K通常小于离散化数量M和N，因此角度-频率联合分布矩阵P是具有低秩属性；
根据角度-频率联合分布矩阵P的低秩属性，将向量化的非相干分布宽带源信号协方差模型公式(7)化为秩最小化问题，具体如下：
式中rank(·)表示秩函数；
秩最小化问题是一个NP难问题，随着低秩矩阵恢复理论的发展，秩函数可以被核范数代替；用估计量替代rs，公式(10)更新为拉格朗日软约束形式：
式中λ为拉格朗日乘子，通常λ∈(0,1)；核范数||P||*为rank(P)非凸替代，表示P的核范数为P的奇异值之和，其中σi为P进行SVD的第i个奇异值；‖·‖2为2范数。
进一步地，步骤S3.2中，公式(11)只能获得角度-频率联合分布矩阵P的一个次优解，结合交替方向乘子法(ADMM)和迭代加权核范数(IRNN)算法对秩最小化问题进一步优化，获得角度-频率联合分布矩阵P的一个全局最优解，具体如下：
S3.2.1、根据交替方向乘子法(ADMM)将公式(11)写成分布式最小化问题：
式中为引入的未知量；公式(12)写成增广拉格朗日形式：
式中和β分别是交替方向乘子法(ADMM)中的双变量和惩罚参数；<·>表示內积算子；‖·‖F表示F范数；
交替方向乘子法(ADMM)分为四步迭代：
S3.2.1.1、更新P：
式中pq，Zq，Yq，βq分别表示P,Z,Y,β在第q次迭代中的估计值；
S3.2.1.2、更新Z：
S3.2.1.3、更新Y：
Yq+1＝Yq+βq(Pq+1-Zq+1)； (42)
S3.2.1.4、更新β：
βq+1＝min(βmax,ξβq)； (43)
βmax是迭代过程中设置的ξ最大值，ξ是更新算子，ξ＞1；
S3.2.2、为了求得P的最优解，利用迭代加权核范数(IRNN)算法解决公式(14)关于P的最小化子问题，求解过程如下：
考虑核范数的一个非凸替代：
式中，是Lipschitz连续可导函数；gδ(x)为非凸替代函数，当gδ(x)＝1-e-x/δ时能更好近似秩函数；在迭代过程中δ设置以较大值开始，δ取100～500，并在每次迭代中使δ＝δ/ρ，ρ＞1，避免初始δ过小使迭代陷入局部最小化；
将公式(18)更新为IRNN最小化形式：
其中μ＞L(f)确保收敛，L(f)是的Lipschitz常数；是gδ(x)的导数；
秩最小化问题获得最优解为：
其中Gq＝U∑VT是Gq的SVD；是广义软边界算子，
S3.2.3、公式(15)关于Z的最小化子问题具有封闭解，求解过程如下：
简化公式(15)中矩阵的F范数为向量的2范数，Z的最小化子问题简化为：
其封闭解为：
Zq+1通过vec(Zq+1)向量化逆操作重组得到；
因此，根据公式(10)初始化变量P0，Z0＝P0，Y0＝0，β0；设置λ，βmax，ξ，δ，ρ，μ；设置d为迭代停止阈值；ADMM中4个迭代变量P，Z，Y，β根据公式(20)、公式(22)、公式(16)和公式(17)迭代更新，当迭代收敛使得||Pq+1-Pq||F/||Pq||F＜d时，输出P作为的角度-频率联合分布矩阵的估计量
进一步地，步骤S4中，根据步骤S3中得到的角度-频率联合分布矩阵的估计量利用基于功率谱密度矩估计的离散网格估计器(参考文献Shahbazpanahi S,Valaee S,Gershman A B.A covariance fitting approach to parametric localization ofmultiple incoherently distributed sources[M].IEEE Press,2004.)，估计出非相干分布宽带源的关键参数，包括角度分布的中心DOA和角度扩展频率分布的中心频率和频率带宽具体如下：
式中,带宽定义为频率分布标准差的倍；Θk和Γk分别为第k个非相干分布宽带源信号感兴趣的角度和频率范围；对于角度和频率分布，当分布为高斯分布或拉普拉斯分布时，η＝1；当分布为均匀分布时，η＝3。
与现有非相干分布宽带源的角度分布估计方法相比，本发明具有以下优点：
(1)本发明无需已知非相干分布宽带源的角度分布参数化模型，能直接估计源的角度分布。
(2)本发明无需已知非相干分布宽带源信号的频率分布，能同时估计频率分布。
(3)本发明无需多维搜索和大规模信号重构，复杂度大为降低。
附图说明
图1是本发明分布宽带源的信号模型图；
图2是本发明参数估计方法的流程图；
图3是本发明仿真实际的角度-频率联合分布矩阵图；
图4是本发明仿真估计的角度-频率联合分布矩阵图；
图5是本发明中心DOA估计的RMSE随信噪比的变化图；
图6是本发明角度扩展估计的RMSE随信噪比的变化图；
图7是本发明中心频率估计的RMSE随信噪比的变化图；
图8是本发明频率带宽估计的RMSE随信噪比的变化图；
图9是本发明方法运行时间随阵元数的变化图。
具体实施方式
下面结合实施例对本发明作进一步详细的描述，本实施例说明了本发明的一种使用方式，但本发明的实施方式不限于此。
实施例：
一种非相干分布宽带源的角度和频率参数估计方法，如图2所示，包括如下步骤：
S1、通过接收阵列获取远场非相干分布宽带源信号，估计非相干分布宽带源信号的协方差并将其向量化，具体步骤如下：
图1是本发明分布宽带源(简称‘DW源’)的信号模型，首先设置一个包括L个阵元的接收阵列，阵元间距为信号最大频率对应的半波长；假设有K个远场非相干分布宽带源信号入射到接收阵列，θ0k，分别是第k个源的角度分布的中心DOA和角度扩展；
接收阵列在t时刻的输出信号为：
式中t＝1,2,...,Q，这里Q为信号采样的快拍数；f和Γ分别是源信号的频率和频率范围；θ和Θ分别是源信号在空间传播的方位角和角度范围；dSk(f)是第k个源信号频率谱的一个测量；γk(θ,t)是第k个源信号在角度方向传播的复增益；是接收阵列关于(θ,f)的方向导数，j为虚数单位，(·)T为转置，τl(θ)为源信号传播到第l个阵元相对于传播到参考阵元(通常设置第一个阵元为参考阵元)的时延差，l＝1,2，...,L；n(t)为高斯白噪声；
估计得到输出信号的协方差为：
式中(·)H表示共轭转置。
由公式(2)可得到非相干分布宽带源信号的协方差估计量为：
式中为高斯白噪声方差，I为单位向量，对进行特征值分解，可将估计量设定为的最小特征值；
向量化的为：
式中vec(·)表示将矩阵逐列堆叠为列向量，
S2、建立向量化的非相干分布宽带源信号协方差模型，根据协方差模型中的方向导数矩阵的行重复性和奇异性，同时对方向导数矩阵和步骤S1中估计的向量化的非相干分布宽带源信号协方差进行降维，以减少后续计算过程的复杂度，具体步骤如下：
S2.1、建立向量化的非相干分布宽带源信号协方差模型rs，构造方向导数矩阵其中rs表示为
将角度范围Θ和频率范围Γ离散化为[θ1,θ2，...，θM]和[f1，f2，...，fN]，这里M和N分别是角度范围和频率范围离散化的数量；式中是在离散点(θm，fn)的值；aθ，f为方向导数a(θ，f)的简化形式；是K个非相干分布宽带源信号的角度-频率联合分布pθ，f在离散点(θm,fn)的值；pθ,f＝vec(P)，为非相干分布宽带源信号的角度-频率联合分布矩阵P的向量化，具体如下：
这里是本发明待估计的未知量，P的求解是估计角度和频率参数的关键；
构造方向导数矩阵这里为复数集合，具体如下：
S2.2、利用方向导数矩阵的行重复性和奇异性，同时对方向导数矩阵和步骤S1中估计的向量化的非相干分布宽带源信号协方差进行降维；
首先，删除方向导数矩阵中重复的行，同时中对应的行，以降低和的维度；
其次，对和进行二次降维；其过程包括：
S2.2.1、对进行奇异值分解(SVD)，即∑A、UA和UA分别是进行SVD的奇异值对角矩阵、左奇异矩阵和右奇异矩阵；
S2.2.2、保留奇异值主成分；设置一个接近于1的奇异值主成分比例ε，当奇异值对角矩阵∑A对角线上前(q+1)个奇异值之和与所有奇异值之和的比例首次大于ε时，保留前q个较大的奇异值，删除其余较小的奇异值，即主成分的奇异值对角矩阵变为主成分对应的左奇异矩阵和右奇异矩阵更新为
S2.2.3、对方向导数矩阵和向量化的非相干分布宽带源信号协方差进行二次降维，令令
S3、将步骤S2中向量化非相干分布宽带源信号协方差rs模型公式化为秩最小化问题，结合ADMM和IRNN算法对秩最小化问题求解，从而得到的角度-频率联合分布矩阵的估计量，具体步骤如下：
S3.1、根据角度-频率联合分布矩阵P的低秩属性，可将rs模型公式化为秩最小化问题，具体如下：
式中rank(·)表示秩函数；
S3.2、结合交替方向乘子法(ADMM)和迭代加权核范数(IRNN)算法对秩最小化问题求解，估计角度-频率联合分布矩阵P；
根据ADMM方法可将公式(8)写成分布式最小化问题
式中λ为拉格朗日乘子，通常λ∈(0,1)；核范数||P||*为rank(P)非凸替代，这里表示P的核范数为P的奇异值之和，其中σi为P进行SVD的第i个奇异值；‖·‖2为2范数；为引入的未知量；公式(9)可写成增广拉格朗日形式：
式中和β分别是ADMM方法中的双变量和惩罚参数；<·>表示內积算子；‖·‖F表示F范数；
ADMM方法分为四步迭代：
S3.2.1、更新待估计矩阵P：
U∑VT＝Gq； (13)
式中q表示迭代次数；公式(11)中，gδ(x)＝1-e-x/δ是IRNN方法中核范数的替代函数；表示函数对x的求导。这里δ设置以较大值开始，让ρ＞1，每次迭代使δ＝δ/ρ，避免初始δ过小使迭代陷入局部最小化；公式(12)中μ＞L(f)确保迭代收敛，L(f)是的Lipschitz常数。公式(13)中U∑VT是Gq的SVD；
S3.2.2、更新未知量Z：
Zq+1可以通过vec(Zq+1)向量化逆操作重组得到；
S3.2.3、更新双变向量Y：
Yq+1＝Yq+βq(Pq+1-Zq+1)； (17)
S3.2.4、更新惩罚参数β：
βq+1＝min(βmax,ξβq)； (18)
βmax是迭代过程中β的最大值，ξ是更新算子，ξ＞1；
因此，根据公式(8)初始化变量P0，Z0＝P0，Y0＝0，β0；设置λ，βmax，ξ，δ，ρ，μ；设置d为迭代停止阈值；ADMM中4个迭代变量P，Z，Y，β根据公式(15)-公式(18)迭代更新，当迭代收敛使得||Pq+1-Pq||F/||Pq||F＜d时，输出P作为的角度-频率联合分布矩阵的估计量
S4、根据步骤S3中得到的角度-频率联合分布矩阵的估计量利用基于功率谱密度矩估计的离散网格估计器(参考文献Shahbazpanahi S,Valaee S,Gershman A B.Acovariance fitting approach to parametric localization of multipleincoherently distributed sources[M].IEEE Press,2004.)，估计出非相干分布宽带源的关键参数，包括角度分布的中心DOA和角度扩展频率分布的中心频率和频率带宽具体如下：
式中,带宽定义为频率分布标准差的倍；Θk和Γk分别为第k个非相干分布宽带源信号感兴趣的角度和频率范围；对于角度和频率分布，当分布为高斯分布或拉普拉斯分布时，η＝1；当分布为均匀分布时，η＝3。
本实施例中，本发明的效果可以通过以下的仿真结果进一步说明，仿真实验条件如下：
远场两个不同角度和频率分布的非相干分布宽带源入射到阵元数L＝30的阵列上。其中第一个源的角度服从高斯分布，频率服从均匀分布；第二个源的角度和频率都服从均匀分布。那么两个源的角度-频率联合分布如下：
这里pk(θ，f)为第k个源的角度-频率联合分布，它可以组成为K个非相干分布宽带源信号的角度-频率联合分布设置第一个源的角度分布的中心DOA和角度扩展为频率分布的中心频率和带宽为(f01，B1)＝(350Hz，100Hz)；第二个源的角度分布参数设置为 频率分布参数设置为(f02,B2)＝(150Hz,100Hz)。角度和频率的观测范围分别设置为Θ＝[20°,50°]，Γ＝[90Hz,420Hz]。联合分布矩阵的角度和频率的离散分辨率分别设置为1°和10Hz。快拍数Q为104,信噪比SNR为10dB。对于迭代程序，初始化δ＝100，并通过δ＝δ/2更新；初始化惩罚参数为β＝0.3，更新参数为ξ＝1.2，βmax＝105；设置λ＝0.3，μ＝103，ε＝1-106，同时迭代停止阈值d＝10-4。
如图3显示了仿真的实际角度-频率联合分布矩阵，图4是通过本发明方法估计的联合分布矩阵，对比两图可知通过本发明方法估计的联合分布矩阵与真实的有较好的一致性。体现了本发明无需已知角度和频率分布，就能估计非相干分布宽带源的角度-频率联合分布的可行性。
2.1关键参数的估计性能随信噪比变化：
为了体现本发明方法有更高的参数估计精度，根据估计的角度-频率联合分布，进一步估计了联合分布的关键参数，并与现有FRFT-DSPE方法和SBL方法比较了参数估计的性能。在仿真中，设置阵元数L＝20，快拍数Q＝100，信噪比SNR由-10dB变化到10dB，其他参数保持与上例相同。进行100次蒙特卡洛实验，分析本发明方法参数估计的RMSE随信噪比SNR的变化。如图5和图6分别显示了角度分布的中心DOA和角度扩展的估计性能。由图可知，本发明方法角度分布参数估计性能接近克拉美罗下界；同时比现有方法有更高估计精度，特别是在低信噪比时。由于FRFT-DSPE方法不对分布宽带源的频率参数进行估计，本发明方法与SBL方法比较了在频率参数估计的性能。如图7和图8所示，与SBL方法相比，本发明方法在频率参数估计方面也有更好的估计精度。
2.2计算复杂度随阵元数量的变化：
为了体现本发明方法有更低的计算复杂度，申请人在仿真实验中分析了发明方法运行时间随阵元数的变化。保持其他参数不变，设置信噪比SNR＝10dB，阵元数L从10增大到50。如图9所示，相比于需要多维搜索的FRFT-DSPE方法和需要大规模信号重构的SBL方法，本发明方法有更低的计算复杂度。
S5、通过步骤S3和S4估计的频率分布和参数信息，用于确定源目标的信号类型，可以进一步识别目标；估计的角度分布和参数信息，用于目标定位，可确定目标方位。
本发明的上述算例仅为详细地说明本发明的计算步骤和计算性能，而并非是对本发明的实施方式的限定。其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
- 上一篇：石墨接头机器人自动装卡簧、装栓机
- 下一篇：一种测向方法和装置
