CN103713276A - 基于最小互熵谱分析的波达方向估计方法 - Google Patents

Abstract

本发明公开了基于最小互熵谱分析的波达方向估计方法。步骤包括：发射换能器发射信号；用接收水听器阵接收发射的声信号的回波；对接收到的回波信号进行平面波模型的建模，然后运用最小互熵谱估计的方法进行处理，得到空间功率谱估计值；对所述的空间功率谱估计值进行分析，峰值处对应横坐标即为目标所在估计角度。本发明只需要一次较少阵元的阵列采样数据即可得到具有高分辨率的谱分析结果。进一步采用倒谱法，通过逆FFT变换提高了最小互熵谱分析算法的收敛速度。本方法较常规空间谱估计方法有更高的分辨力和更小的运算量，能够对阵列信号进行实时处理。算法不依赖于预先估计的信源数目，具有较好的宽容性，较高的分辨力以及极低的旁瓣电平。

Description

基于最小互熵谱分析的波达方向估计方法

技术领域

本发明涉及雷达、声纳及无线通信领域，特别涉及基于最小互熵谱分析的波达方向估计方法。

背景技术

波达方向（Directional of Arrival,DOA）估计是阵列信号处理领域的一个重要研究方向，在通信、雷达和声纳领域都有广泛的应用。在信号定位时，一般需要确定信号的二维到达角（方位角和俯仰角），即2D-DOA。学者们基于空间谱估计提出了各种有效且实用的算法，比如多重信号分类算法（MultipleSignal Classification,MUSIC和旋转不变算法（Estimation of SignalParameter via Rotational Invariance Techniques,ESPRIT）等经典算法来实现2D-DOA估计。然而这些算法的计算量较大，很难调和高分辨率和计算时间长这两者之间的矛盾。Malioutov等人针对均匀线阵的角度估计问题，提出了一种基于接收数据奇异值分解的简化算法。杨雪亚等结合稀疏解问题和二维DFT，提出一种DOA的迭代算法。然而这些算法虽然降低了计算时间，却存在数据盲区的问题。

发明内容

本发明的目的是针对阵列信号的实时处理，提供基于最小互熵谱分析的波达方向估计方法，在只有较少阵元数的阵列采样（快拍）数据情况下得到高分辨率的谱分析结果。

基于最小互熵谱分析的波达方向估计方法，包括：

（1.1）发射换能器发射信号；

（1.2）用接收水听器阵接收发射的声信号的回波；

（1.3）对接收到的回波信号进行平面波模型的建模，然后运用最小互熵谱估计的方法进行处理，得到空间功率谱估计值；

（1.4）对所述的空间功率谱估计值进行分析，峰值处对应横坐标即为目标所在估计角度。

在步骤（1.2）中，所述的接收水听器阵的布阵要求满足半波长布阵，当所述的声信号为宽带信号时，需满足最小波长的半波长布阵要求。

在步骤（1.3）中，所述的最小互熵谱估计包含以下步骤：

（3.1）最小互熵算法是估计一个真实的概率分布，使得它与给定的先验概率分布之间的相对熵最小；

（3.2）利用功率谱密度和自相关函数互为傅立叶变换的关系，用方便计算的自相关函数值来估计功率谱密度，而自相关函数值的数量有限，那么根据最小互熵算法来外推不能通过计算得到的自相关函数值，从而得到更加准确的功率谱密度；

（3.3）利用驾驶协方差矩阵，获得回波信号的高分辨率算法空间功率谱估计，取空间功率谱估计值的峰值对应的横坐标，即为目标的方位估计值。

步骤（1.2）中，所述的接收水听器阵接收到的回波信号只存在较小的多普勒扩展。

步骤（1.3）中，所述的空间功率谱估计要求探测水域混响小。

与现有技术相比，本发明的有益效果是：

传统的DOA估计方法如多重信号分类算法(MUSIC)需要对阵列接收数据的协方差矩阵进行特征分解，并在全空域进行谱峰搜索，运算量巨大，尤其是二维DOA估计方法还存在稳健性较差的问题。

（1）本发明只需要一次快拍数据，即能在信噪比和阵元数都很小的情况下得到较高的空间谱分辨率，体现了该算法在DOA估计中的优势。

（2）采用倒谱法实现的最小互熵算法很大程度上降低了计算量，能够对阵列信号进行实时处理。

（3）不依赖于预先估计的信源数目，而且能够高分辨力辨识它们，具有较好的宽容性。

附图说明

图1是本发明的原理图；

图2是发射换能器和8元接收水听器阵及目标的示意图；

图3是DOA为-45°和60°时不同样本数下MCE和MUSIC算法的空间谱曲线；

图4是DOA为-5°和5°时不同样本数下MCE和MUSIC算法的空间谱曲线；

图5是莫干山湖上实验声速剖面图；

图6是不同样本数下MCE和MUSIC算法的空间谱曲线。

具体实施方式

下面结合附图和具体实例对本发明做进一步的描述。

本发明所采取的具体技术方案，步骤包括：

（1）发射换能器发射信号；

（2）用接收水听器阵接收发射的声信号的回波；

（3）对接收到的回波信号进行平面波模型的建模，然后运用最小互熵谱估计的方法进行处理，得到高分辨率的DOA估计；

（4）对空间功率谱估计值进行分析，其峰值处对应横坐标即为目标所在估计角度。

利用本发明的基于最小互熵谱分析的波达方向估计方法包括如下步骤：

（1）对于N元均匀线阵，阵元间距为d，d≤λ/2，λ为波长，在远场、平面波假设下有M个窄带点源信号以θ_k(k＝1,2,...,M)方向入射，则阵列接收的快拍数信号可表示为X(t)＝A(θ)S(t)+N(t)，其中X(t)_N×1为快拍数据矢量，N(t)_N×1为阵列噪声矢量，阵列噪声假定为空时均独立的高斯白噪声，其均值为0，方差为σ²，S(t)为入射信号复幅度矢量，θ＝[θ₁,...θ_M]为信源方位矢量，A(θ)为驾驶向量矩阵，且A(θ)＝[a(θ₁),a(θ₂),...,a(θ_M)]其中，a(θ_k)＝[1,e^jβk,...,e^j(N-1)βk]^T,k＝1,2,...,M为第k个信源的驾驶矢量，β_k＝-2πdsin(θ_k)λ。阵列的协方差矩阵R定义为R＝E[X(t)X^H(t)]＝AR_sA^H+σ²I，其中，R_s＝E[S(t)S^H(t)]为入射信号的协方差矩阵，I为单位阵。

（2）根据谱分析的空时等效性原理，对于一个给定阵列，不同的入射角对应不同的频率分量。通过对空间谱峰的搜索，就可以获得空间信号的方向信息，从而实现测向。为满足均匀线阵无侧向模糊，通常d≤λ/2，则M个信号的频率为|-dsin(θ_k)/λ|≤1/2，其满足采样定理中最小频率值为1，因此x(n)所组成的离散采样序列满足采样定理的要求。由于噪声为高斯白噪声，且各阵元间的噪声相互独立，其互相关函数为0，频谱为平坦谱，故可以得到M个谱峰。根据谱峰对应的频率即可求出信号的到达角

θ_k＝arcsin(-f_kλ/d)  （1）

（3）最小互熵算法从严格的互熵（也称相对熵或交叉熵）的定义出发，用Lagrange乘数法求其最小优化时的解。设

为真实概率密度，p是我们关于的估计的先验密度。约束条件

其中x∈D，连续情况下互熵定义为

式（2）（3）并不能完全确定

如何从上述约束条件的解的集合Ω里找到的估计，最小互熵为这个推断问题提供了一个一般性的解决方法：对于满足限制条件的所有密度，我们选择的后验q要和先验p有着最小互熵，即

（4）根据拉格朗日乘数法来求最小互熵分布，设

那么令

联合式（2）和（3）得到互熵最小的后验密度

的估计q具有形式

$q\left(x\right)=p\left(x\right)\exp (-\mathrm{\λ}-\underset{k=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\β}}_k{g}_k\left(x\right))---\left(5\right)$

其中λ和β_k均为拉格朗日乘数，由步骤（3）中的约束条件式（2）和（3）决定。

（5）假设时域信号为

其中a_k和b_k为随机变量，f_k为频率，且频率f_k处的功率为

我们用联合概率密度来描述这个随机过程，其中x＝x₁,x₂,...,x_N。那么频率f_k处概率密度为

的功率谱为

把P_k作为

的先验估计，那么把（7）作为概率密度

的先验估计的形式，那么我们可以从M+1个自相关函数R(t_r)中获得有关的新信息，

其中t₀＝0。由于后验功率谱的估计S_k＝∫x_kq(x)dx  （9）

联合式（5），（7），（8）和（9）得

$S_k=\frac{1}{\frac{1}{P_k}+\underset{r=0}{\overset{M}{\mathrm{\Σ}}}2{\mathrm{\β}}_r\cos \left(2\mathrm{\π}{t}_r{f}_k\right)}---\left(10\right)$

β_r是由（6）作为约束条件的拉格朗日常数（用S_k代替

）。具体求解时可以采用大范围收敛的迭代法——连续延拓法求解。

（6）Shore提出的最小互熵算法较为复杂，后来有学者提出了一种新的最小互熵谱估计方法——倒谱法，它首先应用Tzannes的互熵概念推导最小互熵谱公式，证明了Lagrange因子λ_k是被测信号先验信息倒谱和后验信息倒谱之差，即 $\mathrm{IFFT}\left[\underset{k=0}{\overset{\∞}{\mathrm{\Σ}}}{\mathrm{\λ}}_k\cos \mathrm{\ω}{\mathrm{\τ}}_k\right]=\frac{1}{2}\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{\mathrm{\λ}}_k{\mathrm{\δ}}_{n-k}=-{\mathrm{\δ}}_n-\mathrm{IFFT}\left[\log S\left(\mathrm{\ω}\right)\right]+\mathrm{IFFT}\left[\log P\left(\mathrm{\ω}\right)\right]---\left(11\right)$

其中，δ_n-k为单位冲激函数，当n=k是为1，否则为0。IFFT表示逆FFT变换，则IFFT[logS(ω)]为信号序列的倒谱，IFFT[logP(ω)]为信号先验序列的倒谱。

而对于一个因果、稳定和最小相位序列x(n)，和其倒谱

之间存在如下递推关系

$\hat{x}\left(n\right)=\left\{\begin{array}{cc}& \\ 0,& n<0\\ \log \left[x\left(0\right)\right],& n=0\\ \frac{x\left(n\right)}{x\left(0\right)}-\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}\left(\frac{k}{n}\right)\hat{x}\left(n\right)\frac{x(n-k)}{x\left(0\right)},& n>0\end{array}\right.---\left(12\right)$

其中，先验序列倒谱

由给定的先验功率谱P(ω)直接计算。则

$\left.\begin{array}{c}{\mathrm{\&lambda;}}_0=-2-2\hat{x}\left(0\right)+2\hat{p}\left(0\right)\\ {\mathrm{\&lambda;}}_k=\left\{\begin{array}{c}-2\hat{x}\left(k\right)+2\hat{p}\left(k\right),k=1,...,M\\ 0,k>M+1\end{array}\right.\end{array}\right\}---\left(13\right)$

根据公式（11），采用倒谱法，即使在阵元数量很少的情况下也能得到分辨率较高的频谱，采用步骤（2）中的式(1)即可计算信号的波达方向。

实施例1

如图1所示，本发明所使用的探测装置包括干端和湿端，干端包括信号发射机、功率放大器、信号接收机和信号处理器，湿端包括发射换能器和一个8元接收水听器阵，此时发射阵阵元个数为M=1,接收阵阵元个数为N=8。

需要说明的是，虽然图1为方便说明而以单个发射换能器和8元接收水听器为示例，但本发明的发射换能器阵中的发射换能器可以是M元，接收水听器阵中的接收水听器可以是N元，其中，M的取值与信号发射机所发送的相互正交的正交信号的数量相同，且M为≥1的正整数；N是接收阵元的个数，接收阵元的个数是大于1的正整数，即N为≥1的正整数。

发射阵和接收阵的布阵示意图如图2所示，接收阵的阵间距为7.5cm,以发射阵元为参考原点，那么发射驾驶向量a_t(θ)和接收驾驶向量a_r(θ)分别可以用式（14）和式（15）表示，

a_t(θ,f_n)＝exp(-j2πd_tf_nsin(θ)/c),d_t＝[050]^T/100 （14）

a_r(θ,f_n)＝exp(j2πd_rf_nsin(θ)/c),  （15）

d_r＝[-30 -22.5 -15 -7.5 0 7.5 15 22.5]^T/100;

式（14）和式（15）中，f_n为发射信号对应第n个采样点的频率,c为水的声速，实验的实际声速如图3声速剖面所示，θ为目标可能存在的所有方位，T表示矩阵转置运算。

在此系统上进行仿真。图3显示了不同样本数下MCE和MUSIC算法的空间谱曲线，其中来波方向为-45°和60°。容易看出在样本数为1，即只有一次快拍数据时，MUSIC算法失去了有效性；而样本数为1时，MCE算法却能很好地识别出波达方向。仿真结果还表明，达到和MCE相当的辨识度，MUSIC算法需要的样本数至少为16。由于MUSIC算法需要对阵列接收数据的协方差矩阵进行特征分解，并在全空间进行谱峰搜索，所以需要较长的计算时间。仿真结果表明，在CPU为Inter(R)Core(TM)i7，主频为3.07GHz的计算机上，得到大致相同的谱分析结果，MCE算法所需要的时间（0.0069秒）仅为MUSIC算法所需要时间（0.74秒）的1/100左右。

图4显示了来波方向为-5°和5°时MCE和MUSIC算法的空间谱曲线。容易看出，即使只采用了一次快拍数据，MCE算法能够识别两个位置很近的来波角度，而尽管增加样本数，MUSIC算法却不能很好地辨识来波方向。

实施例2

采用浙江大学水声实验室2013年4月在浙江省湖州市莫干山湖的湖试实验数据，阵元数N=8，样本数不定，信源数为1，信源信号为PCW信号，信号持续时间10ms，频率6kHz，采样频率48kHz，来波方向为0°。实验当天声速梯度剖面如图5所示，数据处理结果如图6。容易看出，对于实际采样信号而言，采用一次快拍数据的MCE算法能够很好的辨识出来波方向，而MUSIC算法却不能，这和前面的仿真结果吻合较好。

Claims (5)

1.一种基于最小互熵谱分析的波达方向估计方法，其特征在于，包括：

（1.1）发射换能器发射信号；

（1.2）用接收水听器阵接收发射的声信号的回波；

（1.3）对接收到的回波信号进行平面波模型的建模，然后运用最小互熵谱估计的方法进行处理，得到空间功率谱估计值；

（1.4）对所述的空间功率谱估计值进行分析，峰值处对应横坐标即为目标所在估计角度。

2.根据权利要求1的方法，其特征在于，在步骤（1.2）中，接收水听器阵的布阵要求满足半波长布阵，当所述的声信号为宽带信号时，需满足最小波长的半波长布阵要求。

3.根据权利要求1的方法，其特征是，在步骤（1.3）中，所述的最小互熵谱估计包含以下步骤： 

（3.1）最小互熵算法是估计一个真实的概率分布，使得它与给定的先验概率分布之间的相对熵最小；

（3.2）利用功率谱密度和自相关函数互为傅立叶变换的关系，用方便计算的自相关函数值来估计功率谱密度，而自相关函数值的数量有限，那么根据最小互熵算法来外推不能通过计算得到的自相关函数值，从而得到更加准确的功率谱密度；

（3.3）利用驾驶协方差矩阵，获得回波信号的高分辨率算法空间功率谱估计，取空间功率谱估计值的峰值对应的横坐标，即为目标的方位估计值。

4.根据权利要求1的方法，其特征在于，步骤（1.2）中，接收水听器阵接收到的回波信号只存在较小的多普勒扩展。

5. 根据权利要求1的方法，其特征在于，步骤（1.3）中，所述的空间功率谱估计要求探测水域混响小。

Images