宽带相干信号源DOA估计方法

摘要

本发明属于信号处理领域的范畴，为解决分布式宽带相干信号源的DOA估计问题，同时能够对分布式点源模型中的信号源DOA和由于散射、多径等因素生成的扩展角度进行估计。本发明采用的技术方案是，宽带相干信号源DOA估计方法，步骤如下：步骤一：首先要构建传感器阵列的信号模型；步骤二：对获取到的信号进行分数阶傅立叶变换FrFT；步骤三：对调频域中的信号模型进行简化步骤四：利用TLS-ESPRIT算法进行DOA估计，两个传感器子阵列的信号模型简化后，分别求出它们各自的协方差矩阵，并利用奇异值分解实现信号子空间与噪声子空间的分离，最后得到全部信号源的DOA估计值。本发明主要应用于信号处理场合。

说明书

技术领域

本发明属于信号处理领域的范畴，涉及到利用传感器阵列对多个信号源的波达方向进行精确估计。具体讲,涉及宽带相干信号源DOA估计算法。

背景技术

信号源DOA估计是阵列信号处理领域研究的主要内容之一，在声纳、雷达、地震探测和通信系统中都有着广泛的应用。

人们对DOA估计的研究已经有了很长的一段时间，形成了一系列经典的DOA估计算法。传统的波束形成算法是最早的DOA估计算法之一，它的实现比较简单，但分辨率较低且抗噪声能力较弱。之后，人们提出将子空间分解技术运用到DOA估计中，以多重信号分类算法(MUSIC)和基于旋转不变技术的信号参数估计算法(ESPRIT)这两种算法最为经典和常用，它们都是通过协方差矩阵分解得到信号子空间和噪声子空间，构造功率谱函数来对信号源的到达角进行估计，算法的精确度较高，并且环境适应能力较强。其中ESPRIT算法由于无需预先知道传感器信息，因此使用起来更加的方便，并衍生出一系列的改进算法，比如LS-ESPRIT、TLS-ESPRIT等。但上面介绍的几种算法一般只适合用于窄带信号的处理。

在宽带阵列信号处理中，通常是基于点源信号模型进行空间参数估计，点源信号的能量总是沿着信号的DOA方向传输。在信号传输过程中，由于空间散射以及多径的影响，生成的角度扩展在很多情况下是不能忽视的。因此我们需要采取某种方法能够同时对分布式信号模型中的到达角以及扩展角度进行精确估计。

通常我们可以采用贝叶斯方法和高阶统计量方法进行分布式宽带信号源空间参数估计，但是这些方法在信号源数目较多和点源散射严重等复杂场景中不太适用。实际上MUSIC算法通过改进可以将信号处理由窄带扩展到宽带领域，且能够解决信号散射带来的影响，但是这种算法在位置矢量时变的情况下不再适用，而实际环境中，大多数情况下位置矢量都是时变的。根据时变位置矢量的空间特征，可以采用传统分布式信号源参数估计算法(DSPE)进行处理，但是这种算法要求源信号数目不能超过传感器数目，且信号源必须独立，使得其应用具有很大的局限性。

发明内容

为克服现有技术的不足，本发明旨在解决分布式宽带相干信号源的DOA估计问题，同时能够对分布式点源模型中的信号源DOA和由于散射、多径等因素生成的扩展角度进行估计。本发明采用的技术方案是，宽带相干信号源DOA估计方法，步骤如下：

步骤一：首先要构建传感器阵列的信号模型，假设空间中有Q个信号源入射到一对均匀线性排列的传感器阵列，每组阵元的个数为P，其中P＜Q，阵元间距为d₁,d₂，则t时刻每一组子阵列中第i个传感器处的接收信号模型表示为：

式中，表示传感器处的位置矢量，是第i个信号源的角密度，其中ψ_i为包含中心角和扩展角的未知位置参数，是均值为0，方差为ο²的高斯白噪声，它表示信号x_i处的噪声，为入射角度；

对于相干信号源，它的角密度用进行代换，其中s(t)是一个反映信号源特性的随机变量，为表示信号空间分布的角信号密度，它是一个关于角度θ的复数方程；

步骤二：对获取到的信号进行分数阶傅立叶变换FrFT得到

式中，为传感器阵列的位置矢量矩阵，S_i(α,u)为第i个信号源，N_x(α,u)为传感器阵列处的噪声；

根据分数阶傅立叶变换的旋转特性，调频域与能量集中域是相互垂直的，信号在这两个域中的旋转角度分别表示为

α_d＝arctanμ₀

α_e＝-arccotμ₀

式中，μ₀为信号的调频率，由此可知两个旋转角度满足α_d＝α_e+π/2，因此得到

X(α_d,u)＝F[X(α_e,u)]

式中，F(·)表示分数阶傅里叶变换，上式表示信号在调频域的分数阶傅里叶变换等于其在能量集中域分数阶傅里叶变换的分数阶傅里叶变换，在实际测量中，根据信号的能量值通过峰值搜索的方法进行信号提取，即为能量集中域信号X(α_e,u)，通过分数阶傅立叶变换将其转化到调频域中得到信号X(α_d,u)，从而求得信号源的旋转角度α_d；

步骤三：对调频域中的信号模型进行简化

定义参数b(α_d,ψ_i；u)为

因此，传感器阵列信号模型在调频域中表示为

进而得到

X(α_d,u)＝B(α_d,ψ；u)·S(α_d,u)+N_x(α_d,u)

式中，

B(α_d,ψ；u)＝[b(α_d,ψ₁；u),...,b(α_d,ψ_Q；u)]

S(α_d,u)＝[s₁(α_d,u),...,s_Q(α_d,u)]^T

而另外一组子阵列在调频域中的表达式为

Y(α_d,u)＝C(α_d,ψ；u)S(α_d,u)+N_y(α_d,u)

式中，

C(α_d,ψ；u)＝[c(α_d,ψ₁；u),...,c(α_d,ψ_Q；u)]

对于均匀线性阵列，传感器阵列的位置矢量矩阵为

而其中每个传感器处的位置矢量在调频域中表示为

其中，表示信号到达第k个传感器处的时间延时，c为光速；

由于第i个信号的估计角度与其实际角度θ_0i相近，即因此利用泰勒级数对上式在θ处进行展开

进一步表示为

其中

由于式中为关于角度θ_0i的归一化对称参数，且因此M_n,i中奇数项为0，得到：

$b(\alpha ,{\psi }_i)=\underset{k=0}{\overset{\infty }{\Sigma }}\frac{A^{\left(2k\right)}(\alpha ,{\theta }_{0i})}{\left(2k\right)!}{M}_{2k,i}$

同样的c(α_d,ψ_i)简化后变为

$c({\alpha }_d,{\psi }_i)\approx e^{jf\left(\theta \right)}\underset{k=0}{\overset{\infty }{\Sigma }}\frac{A^{\left(2k\right)}({\alpha }_d,{\theta }_{0i})e^{jf\left({\theta }_{0i}\right)}}{\left(2k\right)!}{M}_{2k,i}$

从而可得到c(α,ψ_i)与b(α,ψ_i)满足如下关系

步骤四：利用TLS-ESPRIT算法进行DOA估计，两个传感器子阵列的信号模型简化后，分别求出它们各自的的协方差矩阵，并利用奇异值分解实现信号子空间与噪声子空间的分离，最后得到全部信号源的DOA估计值。

对步骤一中的信号模型进行处理，可分为如下几个步骤：

1)、将得到的信号利用分数阶傅里叶变换转换到调频域中处理，表达式为

X(α_d,u)＝B(α_d,ψ；u)S(α_d,u)+N_x(α_d,u)

则信号的分数阶相关阵为

R_x(α_d,u)＝E[X(α_d,u)2X(α_d,u)^H]S(α_d,u)+N_x(α_d,u)

＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

+B(α_d,ψ；u)E[S(α_d,u)·N(α_d,u)^H]+E[N_X(α_d,u)·S_X(α_d,u)^H]B(α_d,ψ；u)^H

由于信号与噪声不相关且各阵元噪声互不相关，则

R_x(α_d,u)＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

＝B(α_d,ψ；u)E_SB(α_d,ψ；u)^H+E_N

根据上述公式可分解得到信号子空间E_S和噪声子空间E_N，将上述得到的接收信号在调频域中进行旋转，每个信号得到不同的旋转角度，同时可以得到其在(α,u)域中的能量谱；

2)、在能量谱中进行峰值搜索，得到所有源的数目；

3)、由于调频率不同的源信号具有不同的旋转角度，则其在能量集中域中具有不同的空间表示，但在这个域中，只有一个信号能量值达到最大，而其他的信号和噪声能量相对比较分散，根据这个特性，利用带通滤波器将其分离出来，带通滤波器采用矩形窗函数，它的带宽为相邻信号峰值之间距离的最大值，能够用于消除其他信号源的影响；

4)、根据X(α_d,u)＝F[X(α_e,u)]，对滤波后的信号进行分数阶傅里叶变换，将其转化到调频域中，则分数阶相关阵变为利用特征值分解可分别得到其信号子空间E_S和噪声子空间E_N。

利用TLS-ESPRIT算法进行DOA估计，具体步骤为：

采用基于ESPRIT算法的改进算法TLS-ESPRIT进行DOA以及扩展角度的估计，首先定义一个矩阵Z(α,u)

$Z(\alpha ,u)=\left(\begin{array}{c}X(\alpha ,u)\\ Y(\alpha ,u)\end{array}\right)$

求得其协方差矩阵的信号子空间矩阵E_s(α,u)，根据阵列信号的接收模型，它可表示为E_s(α,u)＝[E₁(α,u)E₂(α,u)]^T，从而可知存在一个列满秩矩阵T，使得

E₁(α,u)＝A₁(α,u)T(α,u)

E₂(α,u)＝A₂(α,u)T(α,u)＝A₁(α,u)Φ(α,u)T(α,u)

式中，由于E₁,E₂都是满秩的，因此存在唯一的一个非奇异矩阵Ψ(α,u)＝T^-1(α,u)Φ(α,u)T(α,u)使得

E₂(α,u)＝E₁(α,u)Ψ(α,u)

在实际环境中，由于噪声以及其他环境因素的影响，用和来代替E₁(α,u)和E₂(α,u)，利用TLS估计对Ψ(α,u)进行处理

${\hat{\Psi }}_{TLS}(\alpha ,u)=-V_{1,2}(\alpha ,u)V_{2,2}(\alpha ,u)$

式中，V_1,2(α,u)，V_2,2(α,u)通过特征分解求出

$\left(\begin{array}{c}{\hat{E}}_1(\alpha ,u)\\ {\hat{E}}_2(\alpha ,u)\end{array}\right)\left[\begin{array}{cc}{\hat{E}}_1(\alpha ,u)& {\hat{E}}_2(\alpha ,u)\end{array}\right]=\left(\begin{array}{cc}{V}_{1,1}(\alpha ,u)& V_{1,2}(\alpha ,u)\\ V_{2,1}(\alpha ,u)& V_{2,2}(\alpha ,u)\end{array}\right)L(\alpha ,u)\left(\begin{array}{cc}{V}_{1,1}^H(\alpha ,u)& V_{1,2}^H(\alpha ,u)\\ V_{2,1}^H(\alpha ,u)& V_{2,2}^H(\alpha ,u)\end{array}\right)$

进而将其中包含的和估计出来，得到估计的位置矢量

$\left(\begin{array}{c}{\hat{A}}_1(\alpha ,u)={\hat{E}}_1(\alpha ,u){\hat{T}}^{-1}(\alpha ,u)={\hat{E}}_2(\alpha ,u){\hat{T}}^{-1}(\alpha ,u){\hat{\Phi }}^{-1}(\alpha ,u)\\ =\frac{1}{2}\{{\hat{E}}_1(\alpha ,u){\hat{T}}^{-1}(\alpha ,u)+{\hat{E}}_2(\alpha ,u){\hat{T}}^{-1}(\alpha ,u){\hat{\Phi }}^{-1}(\alpha ,u)\}\end{array}\right)$

根据因此对其中包含的角度信息进行估计；然后结合信号源在能量集中域的峰值信息可以求解出对应DOA估计值所包含的扩展角度信息。本发明的特点及有益效果是：

1.本发明在信号源数目未知的情况下能够对各个信号源DOA值进行精确估计；

2.本发明对信号模型进行简化处理，减少运算量；

3.本发明克服了传统DOA估计算法中当传感器数目少于源数目时不能进行精确DOA估计的问题；

4.本发明利用改进的TLS-ESPRIT算法进行信号源DOA估计，提高了估计精度。

附图说明：

图1是本发明的方向流程图；

图2是传感器阵列分布示意图；

图3是不同信号变换到调频域中的示意图；

图4是利用本发明中算法将信号变换到调频域的频谱图；

图5是使用带通滤波器进行信号分离的示意图。

具体实施方式

本发明的实现过程如下：

步骤一：首先要构建传感器阵列的信号模型，假设空间中有Q个信号源入射到一对均匀线性排列的传感器阵列，每组阵元的个数为P，则t时刻每一组子阵列中第i个传感器处的接收信号模型可表示为

式中，表示传感器处的位置矢量，是第i个信号源的角密度，其中ψ为包含中心角和扩展角的未知位置参数，是均值为0，方差为ο²的高斯白噪声，它表示信号x_i处的噪声，为入射角度。

对于相干信号源，它的角密度可以用进行代换，其中s(t)是一个反映信号源特性的随机变量，为表示信号空间分布的角信号密度，它是一个关于角度θ的复数方程。

步骤二：对获取到的信号进行分数阶傅立叶变换(FrFT)得到

根据分数阶傅立叶变换的旋转特性，调频域(Dechirping domain)与能量集中域(Energy-concentrated domain)是相互垂直的。信号在这两个域中的旋转角度满足如下关系式

α_d＝arctanμ₀

α_e＝-arccotμ₀

式中，μ₀为信号的调频率。由此可知两个旋转角度满足α_d＝α_e+π/2。因此可以得到

X(α_d,u)＝F[X(α_e,u)]式中，F(·)表示分数阶傅里叶变换，上式表示信号在调频域的分数阶傅里叶变换等于其在能量集中域分数阶傅里叶变换的分数阶傅里叶变换。在实际测量中，我们可以根据信号的能量值通过峰值搜索的方法进行信号提取，即为能量集中域信号X(α_e,u)，通过分数阶傅立叶变换将其转化到调频域中得到信号X(α_d,u)，从而求得信号源的旋转角度α_d。

步骤三：对调频域中的信号模型进行简化。

定义参数b(α,ψ_i；u)为

因此，传感器阵列信号模型在调频域中可表示为

$X({\alpha }_d,u)=\underset{i=1}{\overset{Q}{\Sigma }}b({\alpha }_d,{\psi }_i;u)S_i({\alpha }_d,u)+N_x({\alpha }_d,u)$

进而可以得到

X(α_d,u)＝B(α_d,ψ；u)·S(α_d,u)+N_x(α_d,u)

式中，

B(α_d,ψ；u)＝[b(α_d,ψ₁；u),...,b(α_d,ψ_Q；u)]

S(α_d,u)＝[s₁(α_d,u),...,s_Q(α_d,u)]^T

而另外一组子阵列在调频域中的表达式为

Y(α_d,u)＝C(α_d,ψ；u)·S(α_d,u)+N_y(α_d,u)

式中，

C(α_d,ψ；u)＝[c(α_d,ψ₁；u),...,c(α_d,ψ_Q；u)]

对于均匀线性阵列，传感器阵列的位置矢量矩阵为

而其中每个传感器处的位置矢量在调频域中表示为

其中，表示信号到达第k个传感器处的时间延时，为入射角度，c为光速。

由于第i个信号的估计角度与其实际角度θ_0i相近，即因此可利用泰勒级数对上式在θ处进行展开

进一步表示为

$b(\alpha ,{\psi }_i)=\underset{n=0}{\overset{\infty }{\Sigma }}\frac{A^{\left(n\right)}(\alpha ,{\theta }_{0i})}{n!}{M}_{n,i}$

其中

由于式中为关于角度θ_0i的归一化对称参数，且因此M_n,i中奇数项为0，可得

$b(\alpha ,{\psi }_i)=\underset{k=0}{\overset{\infty }{\Sigma }}\frac{A^{\left(2k\right)}(\alpha ,{\theta }_{0i})}{\left(2k\right)!}{M}_{2k,i}$

同样的c(α_d,ψ_i)简化后变为

$c({\alpha }_d,{\psi }_i)\approx e^{jf\left(\theta \right)}\underset{k=0}{\overset{\infty }{\Sigma }}\frac{A^{\left(2k\right)}({\alpha }_d,{\theta }_{0i})e^{jf\left({\theta }_{0i}\right)}}{\left(2k\right)!}{M}_{2k,i}$

从而可得到c(α,ψ_i)与b(α,ψ_i)满足如下关系

对信号模型进行简化处理后求出其相关函数，并利用奇异值分解得到信号子空间和噪声子空间，并进行峰值搜索和滤波处理，将滤波后信号转换到调频域，再进行子空间分解。

步骤四：利用TLS-ESPRIT算法进行DOA估计。根据步骤三中滤波后得到的信号子空间，利用TLS-ESPRIT算法进行信号源DOA和扩展角度的估计。

本发明提出了一种基于TLS-ESPRIT算法的宽带相干信号源DOA估计算法，具体的实施过程可以分为三个部分。下面结合附图对本发明进行更完整、更详细的描述。

第一部分为信号模型的构建。

假设空间中有Q个相干信号源入射到一对均匀线性排列的传感器阵列，每组阵元的个数为P，传感器的阵列分布如附图2所示，则t时刻每一组子阵列中第i个传感器处的接收信号模型可表示为

传统的DOA估计算法在传感器阵元少于信号源数目的情况下不再适用，本发明为了解决这个问题，因此在参数设计时设定P＜Q。

第二部分主要介绍的是信号模型处理，可分为如下几个步骤。

1、将得到的信号利用分数阶傅里叶变换转换到调频域中，并利用上述步骤三中的方法对信号模型进行简化处理，得到表达式

X(α_d,u)＝B(α_d,ψ；u)S(α_d,u)+N_x(α_d,u)

则信号的分数阶相关阵为

R_x(α_d,u)＝E[X(α_d,u)·X(α_d,u)^H]S(α_d,u)+N_x(α_d,u)

＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

+B(α_d,ψ；u)E[S(α_d,u)·N(α_d,u)H]+E[N_X(α_d,u)·S_X(α_d,u)^H]B(α_d,ψ；u)^H

由于信号与噪声不相关且各阵元噪声互不相关，则

R_x(α_d,u)＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

＝B(α_d,ψ；u)E_SB(α_d,ψ；u)^H+E_N

将上述简化处理后得到的信号在调频域中进行旋转，每个信号得到不同的旋转角度，如图3所示，同时可以得到其在(α,u)域中的能量谱。

2、如附图4所示，在能量谱中进行峰值搜索，可以得到所有源的数目。

3、由于调频率不同的源信号具有不同的旋转角度，则其在能量集中域中具有不同的空间表示。但在这个域中，只有一个信号能量值达到最大，而其他的信号和噪声能量相对比较分散。根据这个特性，利用带通滤波器可将其分离出来。如附图5所示，带通滤波器采用矩形窗函数，它的带宽为相邻信号峰值之间距离的最大值，能够用于消除其他信号源的影响。

4、根据X(α_d,u)＝F[X(α_e,u)]，对滤波后的信号进行分数阶傅里叶变换，将其转化到调频域中，则分数阶相关阵变为利用特征值分解可分别得到其信号子空间E_S和噪声子空间E_N。

第三部分为利用TLS-ESPRIT算法进行DOA估计。

传统的分布式信号源参数估计方法(DSPE)是基于噪声特征矢量最小化来进行参数估计的，即对于要估计的空间分布参数ψ，我们可以通过下式求出其估计值

$\hat{\psi }=\arg \underset{\psi }{max}\frac{1}{\left|\right|b^H\left(\psi \right)E_n\left|\right|}$

式中，

这种算法虽然能够对信号空间参数进行估计，但是它的分辨率较低，且信源数目会对估计精度产生较大的影响，因此在本发明中我们采用基于ESPRIT算法的改进算法TLS-ESPRIT进行DOA以及扩展角度的估计。

在这里首先定义一个矩阵Z(α,u)

$Z(\alpha ,u)=\left(\begin{array}{c}X(\alpha ,u)\\ Y(\alpha ,u)\end{array}\right)$

根据第二部分可求得其协方差矩阵的信号子空间矩阵E_s(α,u)，可表示为E_s(α,u)＝[E₁(α,u)E₂(α,u)]^T。

从而可知存在一个列满秩矩阵T，使得

E₁(α,u)＝A₁(α,u)T(α,u)

E₂(α,u)＝A₂(α,u)T(α,u)＝A₁(α,u)Φ(α,u)T(α,u)

式中，由于E₁,E₂都是满秩的，因此存在唯一的一个非奇异矩阵Ψ(α,u)＝T^-1(α,u)Φ(α,u)T(α,u)使得

E₂(α,u)＝E₁(α,u)Ψ(α,u)

在实际环境中，用和来代替E₁(α,u)和E₂(α,u)。

利用TLS估计对Ψ(α,u)进行处理

${\hat{\Psi }}_{TLS}(\alpha ,u)=-V_{1,2}(\alpha ,u)V_{2,2}(\alpha ,u)$

式中，V_1,2(α,u)，V_2,2(α,u)可通过特征分解求出

$\left(\begin{array}{c}{\hat{E}}_1(\alpha ,u)\\ {\hat{E}}_2(\alpha ,u)\end{array}\right)\left[\begin{array}{cc}{\hat{E}}_1(\alpha ,u)& {\hat{E}}_2(\alpha ,u)\end{array}\right]=\left(\begin{array}{cc}{V}_{1,1}(\alpha ,u)& V_{1,2}(\alpha ,u)\\ V_{2,1}(\alpha ,u)& V_{2,2}(\alpha ,u)\end{array}\right)L(\alpha ,u)\left(\begin{array}{cc}{V}_{1,1}^H(\alpha ,u)& V_{1,2}^H(\alpha ,u)\\ V_{2,1}^H(\alpha ,u)& V_{2,2}^H(\alpha ,u)\end{array}\right)$

进而可以将其中包含的和估计出来。

可得到估计的位置矢量

$\left(\begin{array}{c}{\hat{A}}_1(\alpha ,u)={\hat{E}}_1(\alpha ,u){\hat{T}}^{-1}(\alpha ,u)={\hat{E}}_2(\alpha ,u){\hat{T}}^{-1}(\alpha ,u){\hat{\Phi }}^{-1}(\alpha ,u)\\ =\frac{1}{2}\{{\hat{E}}_1(\alpha ,u){\hat{T}}^{-1}(\alpha ,u)+{\hat{E}}_2(\alpha ,u){\hat{T}}^{-1}(\alpha ,u){\hat{\Phi }}^{-1}(\alpha ,u)\}\end{array}\right)$

根据上文中因此可对其中包含的角度信息进行估计。然后结合信号源在能量集中域的峰值信息可以求解出对应DOA估计值所包含的扩展角度信息。