一种少阵元近场宽带信号源参数估计方法

摘要

本发明提提供一种利用较少阵元对较多的近场宽带信号源的DOA和距离参数进行估计的方法，充分利用宽带信号的特点，在较少阵元的情况下，利用频率族中阵元输出的互相关构造出满足现有窄带MUSIC算法的，含有方位角和距离参数信息的托普利茨矩阵，最后采用MUSIC算法来估计方位角和距离参数，实现对较多的近场宽带信号源进行定位。本发明可利用较少阵元对较多的近场宽带信号源进行定位，且不需参数配对，角度预估和宽带聚焦，运算量更小，便于实际应用。

说明书

技术领域

本发明应用于阵列信号处理领域的近场宽带目标定位技术。

背景技术

基于传感器阵列的辐射源定位是阵列信号处理的一个重要研究内容，它在雷达，声纳，无线通信，地震学以及射电天文学等众多领域中都有着广泛的应用。通常情况下，当信源位置与接收阵列距离较远时，在接收端可以将目标发射信号看成是一个平面波，目标的位置可以通过信源的方位角（DOA）来确定。传统的高分辨DOA估计方法均是基于远场模型下得到的。然而，当信源距离接收阵列较近即信源位于近场条件范围时，平面波的假设不再成立，信号以球面波的形式通过阵列，此时需同时估计出信源的DOA和距离参数才能实现定位，基于远场假设条件下得到的方法不再有效。因此近年来，对近场源参数估计的研究逐渐引起国内外学者的关注，成为阵列信号处理领域中一个新的热点。在过去十几年来，MUSIC（Multiple Signal Classification多信号分类）算法，极大似然(MLE)，二阶统计量和加权线性预测等方法都被用来对近场窄带目标定位。

随着现代信号处理技术的发展，宽带信号在阵列信号处理系统中的应用越来越普遍。目前为止，已经有许多学者对宽带信号的DOA估计进行了研究，并取得了显著的研究成果。近年来，各国相关学者针对近场宽带信源定位问题也相继提出了一系列的定位方法，如最大似然法，二维MUSIC算法，路径跟踪算法，高阶ESPRIT算法等。这些算法或多或少存在着一些不足，例如最大似然方法虽然具有最好的估计性能，但是计算量很大，限制了其在实际工程中的应用。二维MUSIC算法中的ISSM算法在低信噪比条件下估计性能较差，且计算量大，不能估计相干信号源，CSSM类算法通过应用聚焦矩阵，对不同频率点的数据进行聚焦，得到单一频率点（参考频点）的数据，从而计算出信号协方差矩阵，再应用传统窄带子空间方法计算出方位角，但该算法需要对信源位置进行预估，且最终定位结果对预估值敏感。高阶ESPRIT算法计算量较大，并且往往存在参数配对或孔径损失等问题。

上述的近场宽带信源定位方法都是建立在阵元个数个数大于或等于信号源个数条件下的。比如，MUSIC算法通过将接收信号变换到频域，在估计中心频点，并将宽带信号聚焦到中心频点上（以形成窄带），再通过接收向量得到一个自相关矩阵，通过该自相关矩阵得到信号子空间与噪声子空间从而计算得到信源的方位角DOA。MUSIC算法构造的自相关矩阵的秩需要等于信号源的个数，这样，就需要阵元个数个数大于或等于信号源个数。

显然要求阵元个数个数大于或等于信号源个数无法适合于用户日益增长的移动通信等环境中，一般来讲此时的用户是远远大于阵元个数的，而且从物理条件或者经济成本考虑，阵元数通常也是较有限的。然而当信源数大于阵元数时，阵列流型矩阵各列之间不再线性独立，信源协方差矩阵会发生秩亏缺，则此时阵列信号协方差矩阵的大特征值数目小于信源数，不能构成信号子空间，传统的子空间方法无法再使用。

发明内容

本发明所要解决的技术问题是，提供一种利用较少阵元对较多的近场宽带信号源的DOA和距离参数进行估计的方法。

本发明为解决上述技术问题所提供的技术方案是，接收阵列中阵元的数量至少为5个，接收阵列中心的阵元设为参考阵元，包括以下步骤：

（一）对阵列中各阵元的接收信号x_p(t)进行离散傅里叶变换得到近场宽带信号源在频率域的阵列信号模型X_p(f)，之后进入步骤（二）和步骤（五）；

$>X_p\left(f\right)=\underset{k=1}{\overset{K}{\Sigma }}{S}_k\left(f\right)e^{{\mathrm{j\tau }}_{\mathrm{pk}}\left(f\right)}+N_p\left(f\right)$>

其中，p表示各阵元的编号，参考阵元的编号为p=0，与参考阵元为中心向左为负方向，向右为正方向；K为近场非相关宽带信号源的个数，S_k(f)表示第k个信号在频率f上的频谱，N_p(f)表示在频率f上的第p个附加噪声，附加噪声N_p(f)为与信号不相关的零均值空间白噪声，τ_pk(f)表示第k个信号入射到参考阵元相对于第p个阵元在频率f上的相位差；

（二）将近场宽带信号的频率带宽[f_min,f_max]分解成2N+1个频率族，N为正向频率段数的最大值，计算频率族中呈中心对称的阵元输出的含有方位角信息的互相关r_p(f₀+nΔf)：

$>r_p\left(n\right)\stackrel{\Delta }{=}{r}_p(f_0+\mathrm{n\Delta f})=\frac{1}{L}\underset{l=1}{\overset{L}{\Sigma }}{X}_p(f_0+\mathrm{n\Delta f},l)X_{-p}^{*}(f_0+\mathrm{n\Delta f},l);$>

其中，p=1,2；n=-N,…,-1,0,1,…,N；l=1,2…L，L表示频域快拍数；(·)^*表示复数共轭运算；f₀为中心频率，Δf为频率间隔，X_p(f₀+nΔf,l)表示将观测时间划分为L子段，然后对频率f₀+nΔf处的阵列接收信号进行离散傅里叶变换得到的频域信号，表示定义为；

（三）利用步骤（二）中计算得到的(2N+1)个互相关构造托普利茨矩阵R_p，1；

（四）利用步骤（三）得到的托普利茨矩阵R_p，1，使用MUSIC算法计算近场非相关宽带信号源的方位角参数；

（五）将近场宽带信号的频率带宽[f_min,f_max]分解成2N+1个频率族，N为正向频率段数的最大值，计算频率族中呈中心对称的阵元输出的含有方位角和距离参数的信息的互相关r_p(nΔf)：

$>r_{\mathrm{p\Delta f}}\left(n\right)\stackrel{\Delta }{=}{r}_p\left(\mathrm{n\Delta f}\right)=\frac{1}{L}\underset{l=1}{\overset{L}{\Sigma }}{X}_p(f_0+\mathrm{n\Delta f},l)X_{-p}^{*}(f_0-\mathrm{n\Delta f},l);$>

（六）利用步骤（五）中计算得到的(2N+1)个互相关构造托普利茨矩阵R_p,2；

（七）利用步骤（六）得到的托普利茨矩阵R_p，2以及步骤（四）得到的方位角参数，使用MUSIC算法计算近场非相关宽带信号源的距离参数。

本发明充分利用宽带信号的特点，提供了一种少阵元近场宽带信号源方位角和距离二维参数联合估计新方法，在较少阵元的情况下（最少5个阵元），利用频率族中阵元输出的互相关构造出满足现有窄带MUSIC算法的，含有方位角和距离参数信息的托普利茨矩阵（Toeplitz矩阵），最后采用MUSIC算法来估计方位角和距离参数，实现对较多的近场宽带信号源进行定位。

本发明的有益效果是，可利用较少阵元对较多的近场宽带信号源进行定位，且不需参数配对，角度预估和宽带聚焦，运算量更小，便于实际应用。

附图说明

图1：本发明方法流程图。

图2：近场宽带信号接收阵列模型。

图3：两个近场宽带信号源方位角DOA的均方根误差（RMSE）随信噪比变化曲线图。

图4：两个近场宽带信号源距离参数的归一化均方根误差（RNMSE）随信噪比变化曲线图。

图5：本发明方法对多个信号源的估计能力仿真图。图5（a）和（b）分别给出五个和六个近场非相关宽带信号源其入射角分别为[-20°-10°0°10°20°]和[-30°-20°-10°0°10°20°]的伪谱图。

具体实施方式

下面结合附图和具体实施例进一步阐述本发明。这些实施例应理解为仅用于说明本发明而不用于限制本发明的保护范围。在阅读了本发明记载的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等效变化和修饰同样落入本发明权利要求所限定的范围。

实施例1

本发明的DOA和距离参数估计性能仿真：

实施例1的方法如附图1所示，接收阵列为如附图2所示的由5个阵元组成的非均匀线阵，且阵列中心所在的位置设为接收阵列的相位参考点，2个非相关的近场宽带信号的方位角参数以及距离参数分别为[θ₁,r₁]，[θ₂,r₂]，其中[θ₁,r₁]＝[-10°,4]，[θ₂,r₂]＝[5°,2]，宽带信号[f_min,f_max]的相对于中心频率的归一化频谱范围为[0.8,1.2]，即其带宽为中心频率f₀的40%，将宽带信号源分成21个频率族（2N+1），即正向（或反向）频率段数的最大值N＝10，相邻频率族之间的相差为中心阵元为参考阵元，参考阵元右边的第1个阵元与参考阵元的间距d₁以及参考阵元右边的第2个阵元与参考阵元的间距d₂分别为d₂=c/2Δf。其中，r_min=min{r₁,r₂}=2，c为信号传播媒质中的传播速度。

近场非相关宽带信号的信源个数K=2，观测时间划分为L个时间子段，L=16。信噪比SNR从-3dB到15dB变化，进行1000次蒙特卡罗实验。

实施例1中距离参数的估计性能用归一化均方根误差（RNMSE）衡量，方位角DOA的估计性能则用均方根误差（RMSE）衡量，采用的窄带子空间方法为谱MUSIC方法。

DOA和距离参数估计方法包括以下步骤：

（一）对阵列接收信号x_p(t)进行离散傅里叶变换（DFT）得到近场宽带信号源在频率域的阵列信号模型X_p(f)为：

$>X_p\left(f\right)=\underset{k=1}{\overset{K}{\Sigma }}{S}_k\left(f\right)e^{{\mathrm{j\tau }}_{\mathrm{pk}}\left(f\right)}+N_p\left(f\right)$>

其中，本实施例中频域输出对应p=0,±1,±2，p=0表示参考阵元，p=-1，p=1分别表示在参考阵元左右两边第1个阵元，p=-2，p=2分别表示在参考阵元左右两边第2个阵元，S_k(f)表示阵元p在频率f上的第k个信号频谱，N_p(f)表示在f频率族的第p个附加噪声，且为与信号不相关的零均值空间白噪声，τ_pk(f)表示第k个信号入射到参考阵元相对于第p个阵元在频率族f上的相位差，对相位差τ_pk(f)进行Fresnel（菲涅尔）近似后，忽略掉二次项即可得到τ_pk(f)≈ω_pk(f)+φ_pk(f)，其中ω_pk(f)为只含有方位角参数的电角度参数，θ_pk(f)为同时含有方位角和距离参数的电角度参数。

$>{\omega }_{\mathrm{pk}}\left(f\right)=\left(\begin{array}{cc}-p2\mathrm{\pi f}{d}_1\sin \left({\theta }_k\right)/c& p=\pm 1\\ -\mathrm{p\pi f}{d}_2\sin \left({\theta }_k\right)/c& p=\pm 2\\ 0& p=0\end{array}\right)$>

$>{\phi }_{\mathrm{pk}}\left(f\right)=\left(\begin{array}{cc}\mathrm{\pi f}{d}_1^2{\cos }^2\left({\theta }_k\right)/c{r}_k\left(f\right)& p=\pm 1\\ \mathrm{\pi f}{d}_2^2{\cos }^2\left({\theta }_k\right)/{\mathrm{cr}}_k\left(f\right)& p=\pm 2\\ 0& p=0\end{array}\right)$>

式中，c是信号传播媒质中的传播速度，θ_k是第k个信号的方位角DOA，如附图2所示，d₁和d₂分别为参考阵元右边的第1个阵元和第2个阵元所在位置与参考阵元之间的相对距离，r_k(f)=r_kc/f，其中r_k为频率族f上第k个信源的距离参数。

（二）假设宽带信号在给定的归一化频率带宽[0.8,1.2]内有相同的功率谱,且关于中心频率f₀对称，其中归一化中心频率即有S_k(f₀+nΔf)=S_k(f₀-nΔf)k=1,…,K,n=1,…,N，K为信源总数，N为正向（或反向）频率段数的最大值。不失一般性，假定每个宽带信号可以分解成2N+1个频率族，其相邻之间的频率相差为Δf=0.02。

基于以上假设，计算与中心频率f₀对称的阵元的接收信号进行DFT后的频域输出X₁(f₀+nΔf)及X_-1(f₀+nΔf)以及X₂(f₀+nΔf)和X_-2(f₀+nΔf)之间的互相关为：

$>r_p(f_0+\mathrm{n\Delta f})\stackrel{\Delta }{=}E\left\{X_p(f_0+\mathrm{n\Delta f})X_{-p}^{*}(f_0+\mathrm{n\Delta f})\right\}$>

$>=\underset{k=1}{\overset{K}{\Sigma }}{P}_k(f_0+\mathrm{n\Delta f})e^{j2{\omega }_{\mathrm{pk}}(f_0+\mathrm{n\Delta f})}$>p=1,2；n=-N,…,-1,0,1,…,N

其中E{·}和(·)^*分别表示数学期望以及复数共轭运算。P_k(f₀+nΔf)是第k个信号在频率族f₀+nΔf处的功率谱。由此可以看出，在感兴趣的频率段内信号功率谱相同即P_k(f₀+nΔf)与频率族f无关，因此互相关r_p(f₀+nΔf)可以重写为：

$>r_p\left(n\right)\stackrel{\Delta }{=}{r}_p(f_0+\mathrm{n\Delta f})=\underset{k=1}{\overset{K}{\Sigma }}{P}_{k{f}_0}{e}^{j2{\omega }_{\mathrm{pk}}\left(\mathrm{n\Delta f}\right)}\stackrel{\Delta }{=}\underset{k=1}{\overset{K}{\Sigma }}{P}_{k{f}_0}{e}^{\mathrm{jn}{\alpha }_{\mathrm{pk}}};$>p=1,2；n=-N,…,-1,0,1,…,N式中，可看作互相关r_p(n)的幅值，含方位角信息的参数$>{\alpha }_{\mathrm{pk}}\stackrel{\Delta }{=}2{\omega }_{\mathrm{pk}}\left(\mathrm{\Delta f}\right),$>表示定义为。

（三）利用正向2（p=1,2）个阵元的(2N+1)个频率族的互相关构造如下的Toeplitz托普利茨矩阵R_p，1：

其中，矩阵A_p,1＝[a_p,1(θ₁),a_p,1(θ₂),…,a_p,1(θ_K)]是一个(N+1)×K维的范德蒙矩阵，所以当θ_i≠θ_j(i≠j)时，A_p,1是个列向量之间线性无关的列满秩矩阵且其第k列可以表示为a_p，1(θ_k)=[1,exp(-jα_pk),…,exp(-jNα_pk)]^T；$>D\left(P\right)\stackrel{\Delta }{=}\mathrm{diag}\left(P\right)=\mathrm{diag}\{P_{1f_0},P_{2f_0},···,P_{K{f}_0}\}$>是个对角矩阵。由k=1，2,…,K，可以看出所以D(P)的秩为K。通过以上描述可知，矩阵R_p，1的秩等于近场宽带信号源的个数。

（四）利用步骤（三）中所求的Toeplitz矩阵R_p，1，应用窄带MUSIC算法求得DOA，方法如下：

首先对矩阵R_p，1进行奇异值分解（SVD），则零奇异值对应的左奇异矢量构成的(N+1)×(N+1-K)即(11×9)维矩阵U_N1为矩阵R_p,1的噪声子空间。

其次，进行如下谱峰搜索得到信号的方位角DOA：

$>P_1\left({\theta }_k\right)=\frac{1}{a_{p,1}^H\left({\theta }_k\right)U_{N1}{U}_{N1}^H{a}_{p,1}\left({\theta }_k\right)}$>

其中a_p，1(θ_k)=[1,exp(-jα_pk),…,exp(-jNα_pk)]^T，含方位角信息的参数ω_pk如步骤（一）中所述，Δf为频率间隔，(·)^T和(·)^H分别表示转置和共轭转置，P₁(θ)为功率谱。

（五）计算与中心频率f₀对称的阵元的接收信号进行DFT后的频域输出X₁(f₀+nΔf)和X_-1(f₀-nΔf)或者X₂(f₀+nΔf)和X_-2(f₀-nΔf)之间的互相关为：

$>r_p\left(\mathrm{n\Delta f}\right)\stackrel{\Delta }{=}E\left\{X_p(f_0+\mathrm{n\Delta f})X_{-p}^{*}(f_0-\mathrm{n\Delta f})\right\}$>

$>=\underset{k=1}{\overset{K}{\Sigma }}{P}_k\left(\mathrm{n\Delta f}\right)e^{j2{\omega }_{\mathrm{pk}}\left(f_0\right)}{e}^{\mathrm{j\pi }{d}_p^2{\cos }^2\left({\theta }_k\right)4f_0\mathrm{n\Delta f}/c^2{r}_k}$>p=1,2n=-N,…,-1,0,1,…,N式中E{·}和(·)^*分别表示数学期望以及复数共轭运算，f₀为中心频率，Δf为频率间隔，由步骤（二）可知功率谱与频率无关，即有$>P_k\left(\mathrm{n\Delta f}\right)\stackrel{\Delta }{=}E\left\{S_k(f_0+\mathrm{n\Delta f})S_k^{*}(f_0-\mathrm{n\Delta f})\right\}=P_k(f_0+\mathrm{n\Delta f})=P_k(f_0-\mathrm{n\Delta f})=P_k,$>则以上互相关可以重写为

$>r_{\mathrm{p\Delta f}}\left(n\right)\stackrel{\Delta }{=}{r}_p\left(\mathrm{n\Delta f}\right)=\underset{k=1}{\overset{K}{\Sigma }}{P}_{k{f}_0}{e}^{\mathrm{j\pi }{d}_p^2{\cos }^2\left({\theta }_k\right)4f_0\mathrm{n\Delta f}/c^2{r}_k}\stackrel{\Delta }{=}\underset{k=1}{\overset{K}{\Sigma }}{P}_{k{f}_0}{e}^{\mathrm{jn}{\beta }_{\mathrm{pk}}};$>p=1,2；n=-N,…,-1,0,1,…,N其中

$>P_{k{f}_0}\stackrel{\Delta }{=}{P}_k{e}^{j2{\omega }_{\mathrm{pk}}\left(f_0\right)},$>$>{\beta }_{\mathrm{pk}}\stackrel{\Delta }{=}\pi d_p^2{\cos }^2\left({\theta }_k\right)4f_0\mathrm{\Delta f}/c^2{r}_k,$>式中d_p表示第p个阵元所在位置与参考阵元之间的相对距离，c是信号传播媒质中的传播速度，θ_k是第k个信号的方位角DOA，r_k为第k个信源的距离参数。

（六）正向2（p=1,2）个阵元的(2N+1)个频率族的互相关构造如下的Toeplitz矩阵R_p，2：

式中A_p，2=[a_p，2(θ₁,r₁),a_p，2(θ₂,r₂),…,a_p，2(θ_K,r_K)]；

a_p，2(θ_k,r_k)=[1,exp(-jβ_pk),…,exp(-jNβ_pk)]^T，$>{\beta }_{\mathrm{pk}}\stackrel{\Delta }{=}\pi d_p^2{\cos }^2\left({\theta }_k\right)4f_0\mathrm{\Delta f}/c^2{r}_k,$>$>D\left(P\right)\stackrel{\Delta }{=}\mathrm{diag}\left(P\right)=\mathrm{diag}\{P_{1f_0},P_{2f_0},···,P_{K{f}_0}\};$>

其中θ_k是第k个信号的方位角DOA，r_k为第k个信源的距离参数，β_pk为包含方位角和距离信息的参数，d_p表示第p个阵元所在位置与参考阵元之间的相对距离，c是信号传播媒质中的传播速度，f₀为中心频率，Δf为频率间隔，表示互相关r_pΔf(n)的幅值。

（七）结合步骤（六）中所求的Toeplitz矩阵R_p,2以及步骤（四）中求得的方位角DOA，应用窄带MUSIC算法即可求得近场宽带信号源的距离参数，方法如下：

（1）对矩阵R_p,2进行奇异值分解（SVD），则零奇异值对应的左奇异矢量构成的(N+1)×(N+1-K)维矩阵U_N2为矩阵R_p,2的噪声子空间。

（2）进行如下谱峰搜索得到

$>P_2\left({\beta }_{\mathrm{pk}}\right)=\frac{1}{a_{p,2}^H\left({\beta }_{\mathrm{pk}}\right)U_{N2}{U}_{N2}^H{a}_{p,2}\left({\beta }_{\mathrm{pk}}\right)}$>

式中p=1,2，P₂(β_pk)为功率谱，第k个列向量a_p，2(β_pk)=[1,exp(-jβ_pk),…,exp(-jNβ_pk)]^T，其中，β_pk为包含方位角和距离信息的参数，(·)^T和(·)^H分别表示转置和共轭转置，K为近场宽带信号源个数。

（3）利用下式计算近场宽带信号源的距离参数：

$>r_k=4\pi d_p^2{\cos }^2\left({\theta }_k\right)f_0\mathrm{\Delta f}/c^2{\beta }_{\mathrm{pk}}$>

式中p=1,2，d_p表示第p个阵元所在位置与参考阵元之间的相对距离，c是信号传播媒质中的传播速度，θ_k是第k个信源的方位角DOA，f₀为中心频率，Δf为频率间隔，r_k为第k个信源的距离参数。

图3表示方位角DOA的RMSE随SNR=-3dB到SNR=15dB变化曲线图。图4表示距离参数的RNMSE随SNR=-3dB到SNR=15dB变化曲线图。从图3可以看出，DOA的RMSE即使在信噪比（SNR）很低的情况下也非常小，例如，当SNR=-2dB时，RMSE仅仅0.2°。从图4可以看出，距离参数越小，RNMSE越低，即估计性能越好。

实施例2：

本发明对于多个信号源的估计能力：

实施例2的方法如附图1所示，5个近场非相关宽带信号源其入射角分别为[-20°-10°0°10°20°]，SNB=5dB，6个近场非相关宽带信号源其入射角分别为[-30°-20°-10°0°10°20°]，SNB=10dB，其余仿真条件与实施例1的相同，改变仿真条件后再次执行实施例1的步骤即可得到图5，其中图5（a）和（b）分别给出了五个和六个近场非相关宽带信号源其入射角分别为[-20°-10°0°10°20°]和[-30°-20°-10°0°10°20°]的伪谱图。从图5中可以看出所提算法可以估计超过阵元个数的近场宽带信号，而这是传统的近场宽带空间谱算法是无法达到的。