CN106249196A - 三分量声矢量传感器稀疏阵列四元数解模糊方法 - Google Patents

Abstract

三分量声矢量传感器稀疏阵列四元数解模糊方法，阵列接收K个不同频率远场、互不相关窄带声信号，两组接收数据都按照同阵元的同一次采样的z方向振速加乘以x方向振速加乘以y方向振速数据叠加构成四元数一次采样数据；对全阵列接收数据自相关矩阵特征分解并利用子空间理论得到阵列导向矢量估计值；由阵列导向矢量估计值计算子阵导向矢量并由子阵导向矢量间的旋转不变关系得到信号方向余弦的粗略估计值；利用方向余弦粗略估计值解稀疏空域导向矢量的相位模糊，从而得到高精度的信号到达角估计；本发明通过四元数子阵的构造，获得了俯仰角和方位角解耦的旋转关系矩阵，简化了到达角估计且可以进行全方位角度估计，解决了稀疏阵列的相位模糊。

Description

三分量声矢量传感器稀疏阵列四元数解模糊方法

技术领域

本发明属于声信号处理技术领域，尤其涉及三分量声矢量阵列的测向方法。

背景技术

声信号处理是信息处理研究领域中最为活跃的分支之一，传统的声传感器是由无指向性的声压传感器构成，它仅能测量声压强度；随着水声工程等应用的不断发展，对声系统处理提出了越来越高的要求，声矢量传感器技术应运而生并得到不断发展，声矢量传感器是由传统的无指向性的声压传感器和具有偶极子指向性的质点振速传感器按照几何、质量和相位三中心合一的原则构成，声矢量传感器能测量声场中某一点的声压强度和质点振速信息。

传统的长矢量模型处理方法是将声矢量传感器输出的各个分量按照一定顺序构成一个长的复数矢量，并没有考虑声矢量传感器各输出分量之间的正交关系，因而不能充分体现声矢量传感器技术的优越性。四元数具有多个虚部，因而一个四元数可以描述声矢量传感器的多个分量，基于四元数的算法能够利用多个虚部间的代数运算关系对声矢量传感器阵列信号在振速矢量域上的数据分量进行关联，从而可以更加有效地对阵列局部矢量信息加以利用。

基于四元数理论的声矢量传感器测向估计具有重要的意义，现有的方法只是简单的将各个分量叠加到一起，然后利用子阵空域导向矢量间的平移不变关系估计信号的波达方向，没有充分发挥声矢量传感器的矢量特性，本发明要充分利用各分量间的旋转不变关系，通过解模糊处理提高参数估计精度，解决了稀疏阵列的相位模糊问题。

发明内容

针对上述技术的不足和缺陷，本发明提出了一种三分量声矢量传感器稀疏阵列四元数解模糊方法，该方法简化了阵列结构、扩大了阵列孔径且增加了布阵的灵活性，通过对稀疏阵列空域导向矢量进行相位解模糊处理提高了到达角的估计精度。

实现本发明的技术思路是：将阵列接收数据构造成四元数数据，计算阵列接收数据自相关矩阵，通过数据自相关矩阵特征分解及相应的矩阵运算得到阵列导向矢量，由阵列导向矢量重构子阵导向矢量，利用重构后两个子阵导向矢量间的旋转不变关系，得到粗略而无模糊的信号到达角估计，然后利用该粗略估计值解阵列空域导向矢量的相位模糊，从而得到精确而无模糊的信号到达角估计；

本发明的目的是提供一种三分量声矢量传感器稀疏阵列四元数到达角估计方法，通过解模糊处理提高到达角的估计精度；

为了实现上述目的，本发明采取如下的技术解决方案：

三分量声矢量传感器稀疏阵列四元数解模糊方法，包括以下步骤：K个不同频率的远场、互不相关窄带声源信号同时入射到半径为R的稀疏圆形阵列上，

步骤一、将M个由x轴、y轴、z轴方向的振速传感器构成的三分量声矢量传感器组成一个圆形接收阵列，M个阵元特性完全相同，所有x轴振速传感器分量相互平行，所有y轴振速传感器分量相互平行，所有z轴振速传感器分量相互平行；将阵列直接输出的信号和延时ΔT后的输出信号进行N次同步采样得到两组采样数据X和Y，X和Y均为3M×N的矩阵，K＜M-1，采样数据X和Y的各个组成分量叠加构成四元数数据Z₁和Z₂；

其中，表示第m个阵元的x轴方向质点振速传感器输出信号的第n次采样数据，表示第m个阵元的y轴方向质点振速传感器输出信号的第n次采样数据，表示第m个阵元的z轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的x轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的y轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的z轴方向质点振速传感器输出信号的第n次采样数据，m＝1，...，M表示阵元数，n＝1，...，N表示采样数；

将第一组采样数据X和第二组采样数据Y都按照同阵元的同次快拍数据中的z轴方向质点振速传感器接收数据加上四元数虚数单位乘以x轴方向质点振速传感器接收数据，再加上四元数虚数单位乘以y轴方向质点振速传感器接收数据，叠加构成两组四元数N次同步采样接收数据矩阵Z₁和Z₂：

其中，表示由第m个阵元的z轴、x轴及y轴方向质点振速传感器输出信号的第n次采样数据叠加构成第m个阵元的第n次采样的四元数数据，表示由第m个阵元的x轴、y轴z轴方向质点振速传感器接收信号延时ΔT后的第n次采样数据叠加构成第m个阵元延时ΔT后的第n次采样的四元数数据；

第一组采样数据X构成的四元数数据矩阵Z₁＝B₁P+N₁，其中，P＝[p₁，p₂，…，p_K]为K个互不相关入射声波的声强，N₁为四元数噪声，为M×K的导向矢量矩阵，A₁＝[A₁₁，…，A_1k，…，A_1K]表示z轴方向的振速矢量构成的子阵导向矢量，A₂＝[A₂₁，...，A_2k，...，A_2K]表示y轴和x轴方向的振速矢量构成的子阵导向矢量，其中 a_1k＝cosθ， 是M个传感器与位于坐标原点处的传感器之间相位差构成的空域导向矢量，θ_k∈[0，π/2]和φ_k∈[-π，π]是第k个声源信号的俯仰角和方位角，表示第m个传感器的角坐标；

第二组采样数据Y构成的四元数数据矩阵Z₂＝B₂P+N₂，其中，B₂是第二组采样数据对应的阵列导向矢量矩阵，N₂为四元数噪声，两组数据的导向矢量间满足：B₂＝B₁Φ，为时延矩阵；

步骤二、构造全阵列接收数据矩阵Z＝BP+N，其中，Z＝[Z₁；Z₂]，B＝[B₁；B₂]是全阵列接收数据对应的阵列导向矢量矩阵，N＝[N₁；N₂]是全阵列接收数据四元数噪声矩阵；计算全阵列接收四元数数据矩阵Z的自相关矩阵R_z，对数据自相关矩阵进行四元数特征分解，得到信号子空间E_s，按照阵列导向矢量矩阵相同的分块模式对信号子空间进行分块，通过矩阵块运算以及伪伴随矩阵特征分解得到阵列导向矢量估计值及全阵列导向矢量矩阵估计值

$R_z=\frac{1}{N}\[{\mathrm{ZZ}}^H\]={\mathrm{BR}}_P{B}^H+{\mathrm{\σ}}^{2}I,$

其中，为入射信号的自相关函数，σ²为噪声的方差，I为2M×2M单位矩阵，(·)^H表示转置复共轭操作；

对数据自相关矩阵R_z进行四元数矩阵特征分解，得到信号子空间E_s，信号子空间与阵列导向矢量的关系为：E_s＝BT，其中T为K×K的非奇异矩阵，分别取E_s的前M行及后M行组成信号子空间矩阵E₁＝B₁T和E₂＝B₂T＝B₁ΦT，则有其中

将矩阵转换为复数自伴随矩阵其中Ψ_r和Ψ_i的分别是Ψ的实部和虚部，和分别是Ψ_i和Ψ_r的共轭矩阵，是复数的虚数单位，对C_Ψ进行特征分解其中，V_k是特征矢量，U_k是特征值，Ψ的右特征值也是U_k，对应的特征矢量其中，Φ_M＝[I_M，-iI_M]是维数为M×2M的还原矩阵，I_M代表维数为M×M的单位阵；从而由第二次特征分解得到阵列导向矢量的估计值和延时ΔT后的阵列导向矢量估计值及全阵列导向矢量估计值其中，为Ψ的特征向量矩阵；

步骤三、重构子阵导向矢量矩阵的估计和根据步骤一中四元数的构造模式、四元数阵列导向矢量B₁的结构特点以及阵列导向矢量B₁和子阵导向矢量矩阵A₁和A₂关系得到两个子阵导向矢量矩阵的估计和通过矩阵运算得到和的旋转不变关系矩阵Ω，Ω是关于俯仰角和方位角解耦的对角矩阵，由Ω得到信号方位角和俯仰角的粗略估计值，进而得到方向余弦的粗略估计值，该粗略估计值将用于后续的稀疏阵列空域导向矢量矩阵解相位模糊；

根据导向矢量计算子阵导向矢量和根据步骤一中阵列导向矢量B₁和子阵导向矢量矩阵A₁和A₂的关系知：且与满足的关系，Ω＝diag[Ω₁₁，…，Ω_KK]是子阵导向矢量和间的旋转不变关系矩阵，是方位角和俯仰角的解耦函数，由Ω_kk可以很容易地获取声波到达角粗略的估计值和由到达角的粗略估计值可得到方向余弦的粗略而无模糊的估计值

步骤四、获取归一化空域导向矢量的估计值导向矢量矩阵的估计值的每一列对该列的第一个元素进行归一化处理得到归一化空域导向矢量的估计值根据步骤三中得到的方向余弦粗略估计值解空域导向矢量的相位周期性模糊，从而得到信号方位角和俯仰角的精确估计值；

根据导向矢量矩阵计算归一化空域导向矢量其中，表示的第k列，表示的第k列的第1个元素；对空域导向矢量取相位得到相位矩阵arg(·)表示取相位，是有模糊的相位估计值；

由方向余弦的粗略估计值得到稀疏空域导向矢量相位矩阵的粗略估计值通过解这个优化问题得到模糊数是r_opt，为阵列空域导向矢量的精确相位估计值，根据计算第k个声波信号的x轴方向的方向余弦的估计值和y轴方向的方向余弦估计值 其中，[W]^#＝[(W)^HW]^-1(W)^H是位置矩阵W的伪逆矩阵，r_opt是模糊数矢量，是方向余弦的精确估计值，位置矩阵W为

$W=\frac{2\mathrm{\π}R}{\mathrm{\λ}}\left[\begin{array}{cc}0& 1\\ sin\mathrm{\Δ}& \cos \mathrm{\Δ}\\ .& .\\ .& .\\ .& .\\ sin\[(M-1)\mathrm{\Δ}\]& cos\[(M-1)\mathrm{\Δ}\]\end{array}\right]$

其中，Δ＝2π/M。

根据方向余弦的精确估计值和得到信号到达角的精确估计值：

$\left\{\begin{array}{c}{\widetilde{\mathrm{\&theta;}}}_k=\arcsin \left(\sqrt{{\widetilde{\mathrm{\&alpha;}}}_k^2+{\widetilde{\mathrm{\&beta;}}}_k^2}\right)\\ \left\{\begin{array}{cc}{\widetilde{\mathrm{\&phi;}}}_k=\arctan \frac{{\widetilde{\mathrm{\&alpha;}}}_k}{{\widetilde{\mathrm{\&beta;}}}_k}& {\widetilde{\mathrm{\&beta;}}}_k\&GreaterEqual;0\\ {\widetilde{\mathrm{\&phi;}}}_k=\mathrm{\&pi;}+\arctan \frac{{\widetilde{\mathrm{\&alpha;}}}_k}{{\widetilde{\mathrm{\&beta;}}}_k}& {\widetilde{\mathrm{\&beta;}}}_k<0\end{array}\right.\end{array}\right.;$

前述步骤中的m＝1，...，M，M为阵列的阵元数，n＝1，...，N，N为采样次数， 为四元数的3个虚数单位。

本发明的阵列为稀疏圆形阵列，声矢量传感器的阵元由三个空间轴向垂直的振速传感器构成，M个完全相同的阵元均匀分布在圆环上，坐标原点位于圆环的圆心。

本发明将三分量声矢量传感器的接收数据构造成四元数数据，利用两个四元数子阵间的旋转不变关系估计方位角和俯仰角，本发明的突出优点是旋转不变关系矩阵的对角线元素是俯仰角和方位角的解耦函数，大大方便了俯仰角和方位角的估计且可以进行全方位角度估计，通过四元数的解模糊处理提高了参数的估计精度，与现有技术的长矢量方法相比，本发明方法不需要参数配对运算，降低了算法复杂度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中需要使用的附图做简单介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例声矢量传感器阵列的示意图。

图2为本发明方法的流程图。

图3为长矢量不解模糊算法到达角估计散布图。

图4为四元数不解模糊算法到达角估计散布图。

图5为四元数解模糊算法到达角估计散布图。

图6为仿真实验的俯仰角估计均方根误差与信噪比的关系图。

图7为仿真实验的方位角估计均方根误差与信噪比的关系图。

图8为仿真实验的到达角估计均方根误差与信噪比的关系图。

图9为仿真实验的到达角估计成功概率随信噪比的关系图。

具体实施方式

为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明实施例，并配合所附图示，做详细说明如下。

图1所示为本发明的声矢量传感器阵列的示意图。如图1所示，声矢量传感器阵列为圆环形阵列，其阵元是由三个空间轴向垂直的振速传感器构成，M个阵元在半径为R的圆环上均匀分布，第1个阵元位于x轴上，沿圆周逆时针方向分别是第1，...，M个阵元，坐标原点位于圆环的圆心，第m个阵元与x轴正向的夹角其中，R＞＞λ_min/2，λ_min为入射信号的最小波长；

结合图2，本发明的三分量声矢量传感器稀疏阵列四元数解模糊方法，包括以下步骤：K个不同频率远场、互不相关窄带声源信号同时入射到声矢量传感器阵列上，

步骤一、将M个由x轴、y轴、z轴方向的振速传感器构成的三分量声矢量传感器组成一个圆形接收阵列，M个阵元特性完全相同，所有x轴振速传感器分量相互平行，所有y轴振速传感器分量相互平行，所有z轴振速传感器分量相互平行；将阵列直接输出的信号和延时ΔT后的输出信号进行N次同步采样得到两组采样数据X和Y，X和Y均为3M×N的矩阵，K＜M-1，采样数据X和Y的各个组成分量叠加构成四元数数据Z₁和Z₂；

其中，表示第m个阵元的x轴方向质点振速传感器输出信号的第n次采样数据，表示第m个阵元的y轴方向质点振速传感器输出信号的第n次采样数据，表示第m个阵元的z轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的x轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的y轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的z轴方向质点振速传感器输出信号的第n次采样数据，m＝1，...，M表示阵元数，n＝1，...，N表示采样数；

将第一组采样数据X和第二组采样数据Y都按照同阵元的同次快拍数据中的z轴方向质点振速传感器接收数据加上四元数虚数单位乘以x轴方向质点振速传感器接收数据，再加上四元数虚数单位乘以y轴方向质点振速传感器接收数据，叠加构成两组四元数N次同步采样接收数据矩阵Z₁和Z₂：

其中，表示由第m个阵元的z轴、x轴及y轴方向质点振速传感器输出信号的第n次采样数据叠加构成第m个阵元的第n次采样的四元数数据，表示由第m个阵元的x轴、y轴z轴方向质点振速传感器接收信号延时ΔT后的第n次采样数据叠加构成第m个阵元延时ΔT后的第n次采样的四元数数据；

第一组采样数据X构成的四元数数据矩阵Z₁＝B₁P+N₁，其中，P＝[p₁，p₂，…，p_K]为K个互不相关入射声波的声强，N₁为四元数噪声，为M×K的导向矢量矩阵，A₁＝[A₁₁，…，A_1k，…，A_1K]表示z轴方向的振速矢量构成的子阵导向矢量，A₂＝[A₂₁，...，A_2k，...，A_2K]表示y轴和x轴方向的振速矢量构成的子阵导向矢量，其中 a_1k＝cosθ， 是M个传感器与位于坐标原点处的传感器之间相位差构成的空域导向矢量，θ_k∈[0，π/2]和φ_k∈[-π，π]是第k个声源信号的俯仰角和方位角，表示第m个传感器的角坐标；

第二组采样数据Y构成的四元数数据矩阵Z₂＝B₂P+N₂，其中，B₂是第二组采样数据对应的阵列导向矢量矩阵，N₂为四元数噪声，两组数据的导向矢量间满足：B₂＝B₁Φ，为时延矩阵；

步骤二、构造全阵列接收数据矩阵Z＝BP+N，其中，Z＝[Z₁；Z₂]，B＝[B₁；B₂]是全阵列接收数据对应的阵列导向矢量矩阵，N＝[N₁；N₂]是全阵列接收数据四元数噪声矩阵；计算全阵列接收四元数数据矩阵Z的自相关矩阵R_z，对数据自相关矩阵进行四元数特征分解，得到信号子空间E_s，按照阵列导向矢量矩阵相同的分块模式对信号子空间进行分块，通过矩阵块运算以及伪伴随矩阵特征分解得到阵列导向矢量估计值及全阵列导向矢量矩阵估计值

$R_z=\frac{1}{N}\&lsqb;{\mathrm{ZZ}}^H\&rsqb;={\mathrm{BR}}_P{B}^H+{\mathrm{\&sigma;}}^{2}I,$

其中，为入射信号的自相关函数，σ²为噪声的方差，I为2M×2M单位矩阵，(·)^H表示转置复共轭操作；

对数据自相关矩阵R_z进行四元数矩阵特征分解，得到信号子空间E_s，信号子空间与阵列导向矢量的关系为：E_s＝BT，其中T为K×K的非奇异矩阵，分别取E_s的前M行及后M行组成信号子空间矩阵E₁＝B₁T和E₂＝B₂T＝B₁ΦT，则有其中

将矩阵转换为复数自伴随矩阵其中Ψ_r和Ψ_i的分别是Ψ的实部和虚部，和分别是Ψ_i和Ψ_r的共轭矩阵，是复数的虚数单位，对C_Ψ进行特征分解其中，V_k是特征矢量，U_k是特征值，Ψ的右特征值也是U_k，对应的特征矢量其中，Φ_M＝[I_M，-iI_M]是维数为M×2M的还原矩阵，I_M代表维数为M×M的单位阵；从而由第二次特征分解得到阵列导向矢量的估计值和延时ΔT后的阵列导向矢量估计值及全阵列导向矢量估计值其中，为Ψ的特征向量矩阵；

步骤三、重构子阵导向矢量矩阵的估计和根据步骤一中四元数的构造模式、四元数阵列导向矢量B₁的结构特点以及阵列导向矢量B₁和子阵导向矢量矩阵A₁和A₂关系得到两个子阵导向矢量矩阵的估计和通过矩阵运算得到和的旋转不变关系矩阵Ω，Ω是关于俯仰角和方位角解耦的对角矩阵，由Ω得到信号方位角和俯仰角的粗略估计值，进而得到方向余弦的粗略估计值，该粗略估计值将用于后续的稀疏阵列空域导向矢量矩阵解相位模糊；

根据导向矢量计算子阵导向矢量和根据步骤一中阵列导向矢量B₁和子阵导向矢量矩阵A₁和A₂的关系知：且与满足的关系，Ω＝diag[Ω₁₁，…，Ω_KK]是子阵导向矢量和间的旋转不变关系矩阵，是方位角和俯仰角的解耦函数，由Ω_kk可以很容易地获取声波到达角粗略的估计值和由到达角的粗略估计值可得到方向余弦的粗略而无模糊的估计值

步骤四、获取归一化空域导向矢量的估计值导向矢量矩阵的估计值的每一列对该列的第一个元素进行归一化处理得到归一化空域导向矢量的估计值根据步骤三中得到的方向余弦粗略估计值解空域导向矢量的相位周期性模糊，从而得到信号方位角和俯仰角的精确估计值；

根据导向矢量矩阵计算归一化空域导向矢量其中，表示的第k列，表示的第k列的第1个元素；对空域导向矢量取相位得到相位矩阵arg(·)表示取相位，是有模糊的相位估计值；

由方向余弦的粗略估计值得到稀疏空域导向矢量相位矩阵的粗略估计值通过解这个优化问题得到模糊数是r_opt，为阵列空域导向矢量的精确相位估计值，根据计算第k个声波信号的x轴方向的方向余弦的估计值和y轴方向的方向余弦估计值 其中，[W]^#＝[(W)^HW]^-1(W)^H是位置矩阵W的伪逆矩阵，r_opt是模糊数矢量，是方向余弦的精确估计值，位置矩阵W为

$W=\frac{2\mathrm{\&pi;}R}{\mathrm{\&lambda;}}\left[\begin{array}{cc}0& 1\\ sin\mathrm{\&Delta;}& \cos \mathrm{\&Delta;}\\ .& .\\ .& .\\ .& .\\ sin\&lsqb;(M-1)\mathrm{\&Delta;}\&rsqb;& cos\&lsqb;(M-1)\mathrm{\&Delta;}\&rsqb;\end{array}\right]$

其中，Δ＝2π/M。

根据方向余弦的精确估计值和得到信号到达角的精确估计值：

$\left\{\begin{array}{c}{\widetilde{\mathrm{\&theta;}}}_k=\arcsin \left(\sqrt{{\widetilde{\mathrm{\&alpha;}}}_k^2+{\widetilde{\mathrm{\&beta;}}}_k^2}\right)\\ \left\{\begin{array}{cc}{\widetilde{\mathrm{\&phi;}}}_k=\arctan \frac{{\widetilde{\mathrm{\&alpha;}}}_k}{{\widetilde{\mathrm{\&beta;}}}_k}& {\widetilde{\mathrm{\&beta;}}}_k\&GreaterEqual;0\\ {\widetilde{\mathrm{\&phi;}}}_k=\mathrm{\&pi;}+\arctan \frac{{\widetilde{\mathrm{\&alpha;}}}_k}{{\widetilde{\mathrm{\&beta;}}}_k}& {\widetilde{\mathrm{\&beta;}}}_k<0\end{array}\right.\end{array}\right.;$

前述步骤中的m＝1，...，M，M为阵列的阵元数，n＝1，...，N，N为采样次数， 为四元数的3个虚数单位。

本发明利用两组同步采样数据构造全阵列接收数据并计算自相关矩阵，对自相关矩阵进行四元数特征分解，并根据子空间理论得到阵列导向矢量的估计，根据阵列导向矢量重构两个子阵导向矢量，利用两个子阵导向矢量间的关系得到信号二维到达角的粗略估计值，进而得到x轴方向和y轴方向的方向余弦粗略估计值，通过对空域导向矢量估计值进行相位解模糊，运算得到x轴方向和y轴方向的精确方向余弦估计值，从而得到信号二维到达角的精确估计值。

本发明的效果可以通过以下的仿真结果进一步说明：

仿真实验条件如下：采用半径为R＝2.5λ的均匀圆形阵列作为接收阵列，M＝7个声矢量传感器均匀分布于圆周上。两个不同频率的远场、互不相关窄带声源信号入射到该接收阵列上，入射信号的参数分别为(θ₁，φ₁)＝(10°，70°)，(θ₂，φ₂)＝(25°，45°)，进行512次快拍，200次独立蒙特卡罗试验。

长矢量不解模糊方法是指圆形阵列半径R＝0.5λ时基于最小二乘估计到达角的方法，四元数不解模糊方法是指利用子阵的旋转不变关系求到达角估计的方法，本发明专利的四元数解模糊方法是指利用子阵旋转不变关系的粗略估计值解稀疏圆形阵列的相位模糊，从而得到高精度的到达角估计的方法；

到达角估计均方根误差是指的大小，到达角估计成功概率是指在200次独立试验中俯仰角和方位角估计值满足关系式的实验次数占总实验次数的百分比；其中，和是估计值，θ₀和φ₀是真值，mean是指求平均，N是实验次数，和是指第i次实验的估计值；

上述所述的三种算法的仿真对比结果如图3至图9所示；

图3，图4和图5为信噪比是10dB时，长矢量不解模糊算法，四元数不解模糊算法和四元数解模糊算法的到达角估计散布图，从图4至图5可以看出本发明的四元数解模糊算法的到达角估计精度高于四元数不解模糊算法的到达角估计精度，这是因为本发明方法利用了稀疏阵列的孔径信息，提高了参数估计精度；从图3至图4可以看出利用子阵旋转不变关系估计到达角的四元数不解模糊算法性能优于利用非稀疏阵空域导向矢量估计到达角的长矢量方法；

从图6，图7和图8可以看出，在信噪比区间内本发明方法的均方根误差估计性能在三种算法中是最好的，因为它即利用子阵间的旋转不变关系，又利用了稀疏阵列孔径信息，四元数不解模糊方法只利用子阵间的旋转不变关系，长矢量不解模糊方法只利用非稀疏阵列的孔径信息；从图8可以看出，在信噪比为0dB时，长矢量不解模糊算法，四元数不解模糊算法和本发明方法的到达角均方根误差分别为2°，1.3°和0.25°；

从图9可以看出，本发明方法的到达角估计成功概率高于四元数不解模糊算法和长矢量不解模糊算法，特别是0dB时，本发明方法的成功概率超过了99％，而四元数不解模糊算法的成功概率为50％，长矢量不解模糊算法的成功概率低于30％；

因此从参数估计的散布图、均方根误差和成功概率上可以看出本发明方法比四元数不解模糊算法和长矢量不解模糊算法具有更高的参数估计性能；

以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (1)

1.一种三分量声矢量传感器稀疏阵列四元数解模糊方法，其特征在于：

所述声矢量传感器阵列的阵元由三个空间轴向垂直的振速传感器组成，M个阵元均匀分布在半径为R的圆环上，第1个阵元位于x轴上，沿圆环逆时针方向分别是第1，...，M个阵元，坐标原点位于圆环的圆心，第m个阵元与x轴正向的夹角为其中，R＞＞λ_min/2，λ_min为入射信号的最小波长；

三分量声矢量传感器稀疏阵列四元数解模糊方法的步骤如下：阵列接收K个不同频率的远场、互不相关的窄带声波入射信号，

步骤一、将M个由x轴、y轴、z轴方向的振速传感器构成的三分量声矢量传感器组成一个圆形接收阵列，M个阵元特性完全相同，所有x轴振速传感器分量相互平行，所有y轴振速传感器分量相互平行，所有z轴振速传感器分量相互平行；将阵列直接输出的信号和延时ΔT后的输出信号进行N次同步采样得到两组采样数据X和Y，X和Y均为3M×N的矩阵，K＜M-1，采样数据X和Y的各个组成分量叠加构成四元数数据Z₁和Z₂；

$X=\left[\begin{array}{ccccc}{x}_{1v_x}\left(1\right),& ...,& x_{1v_x}\left(n\right),& ...,& x_{1v_x}\left(N\right)\\ x_{1v_y}\left(1\right),& ...,& x_{1v_y}\left(n\right),& ...,& x_{1v_y}\left(N\right)\\ x_{1v_z}\left(1\right),& ...,& x_{1v_z}\left(n\right),& ...,& x_{1v_z}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .,\\ x_{{\mathrm{mv}}_x}\left(1\right),& ...,& x_{{\mathrm{mv}}_x}\left(n\right),& ...,& x_{{\mathrm{mv}}_x}\left(N\right)\\ x_{{\mathrm{mv}}_y}\left(1\right),& ...,& x_{{\mathrm{mv}}_y}\left(n\right),& ...,& x_{{\mathrm{mv}}_y}\left(N\right)\\ x_{{\mathrm{mv}}_z}\left(1\right),& ...,& x_{{\mathrm{mv}}_z}\left(n\right),& ...,& x_{{\mathrm{mv}}_z}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .,\\ x_{{\mathrm{Mv}}_x}\left(1\right),& ...,& x_{{\mathrm{Mv}}_x}\left(n\right),& ...,& x_{{\mathrm{Mv}}_x}\left(N\right)\\ x_{{\mathrm{Mv}}_y}\left(1\right),& ...,& x_{{\mathrm{Mv}}_y}\left(n\right),& ...,& x_{{\mathrm{Mv}}_y}\left(N\right)\\ x_{{\mathrm{Mv}}_z}\left(1\right),& ...,& x_{{\mathrm{Mv}}_z}\left(n\right),& ...,& x_{{\mathrm{Mv}}_z}\left(N\right)\end{array}\right],Y=\left[\begin{array}{ccccc}{y}_{1v_x}\left(1\right),& ...,& y_{1v_x}\left(n\right),& ...,& y_{1v_x}\left(N\right)\\ y_{1v_y}\left(1\right),& ...,& y_{1v_y}\left(n\right),& ...,& y_{1v_y}\left(N\right)\\ y_{1v_z}\left(1\right),& ...,& y_{1v_z}\left(n\right),& ...,& y_{1v_z}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .,\\ y_{{\mathrm{mv}}_x}\left(1\right),& ...,& y_{{\mathrm{mv}}_x}\left(n\right),& ...,& y_{{\mathrm{mv}}_x}\left(N\right)\\ y_{{\mathrm{mv}}_y}\left(1\right),& ...,& y_{{\mathrm{mv}}_y}\left(n\right),& ...,& y_{{\mathrm{mv}}_y}\left(N\right)\\ y_{{\mathrm{mv}}_z}\left(1\right),& ...,& y_{{\mathrm{mv}}_z}\left(n\right),& ...,& y_{{\mathrm{mv}}_z}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .,\\ y_{{\mathrm{Mv}}_x}\left(1\right),& ...,& y_{{\mathrm{Mv}}_x}\left(n\right),& ...,& y_{{\mathrm{Mv}}_x}\left(N\right)\\ y_{{\mathrm{Mv}}_y}\left(1\right),& ...,& y_{{\mathrm{Mv}}_y}\left(n\right),& ...,& y_{{\mathrm{Mv}}_y}\left(N\right)\\ y_{{\mathrm{Mv}}_z}\left(1\right),& ...,& y_{{\mathrm{Mv}}_z}\left(n\right),& ...,& y_{{\mathrm{Mv}}_z}\left(N\right)\end{array}\right]$

其中，表示第m个阵元的x轴方向质点振速传感器输出信号的第n次采样数据，表示第m个阵元的y轴方向质点振速传感器输出信号的第n次采样数据，表示第m个阵元的z轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的x轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的y轴方向质点振速传感器输出信号的第n次采样数据，表示延时ΔT后第m个阵元的z轴方向质点振速传感器输出信号的第n次采样数据，m＝1，...，M表示阵元数，n＝1，...，N表示采样数；

将第一组采样数据X和第二组采样数据Y都按照同阵元的同次快拍数据中的z轴方向质点振速传感器接收数据加上四元数虚数单位乘以x轴方向质点振速传感器接收数据，再加上四元数虚数单位乘以y轴方向质点振速传感器接收数据，叠加构成两组四元数N次同步采样接收数据矩阵Z₁和Z₂：

$Z_1=\left[\begin{array}{ccccc}{z}_{11}\left(1\right),& ...,& z_{11}\left(n\right),& ...,& z_{11}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .\\ z_{1m}\left(1\right),& ...,& z_{1m}\left(n\right),& & z_{1m}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .& ...,& .,& ...,& .\\ z_{1M}\left(1\right),& ...,& z_{1M}\left(n\right),& ...,& z_{1M}\left(N\right)\end{array}\right],Z_2=\left[\begin{array}{ccccc}{z}_{21}\left(1\right),& ...,& z_{21}\left(n\right),& ...,& z_{21}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .\\ z_{2m}\left(1\right),& ...,& z_{2m}\left(n\right),& & z_{2m}\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .& ...,& .,& ...,& .\\ z_{2M}\left(1\right),& ...,& z_{2M}\left(n\right),& ...,& z_{2M}\left(N\right)\end{array}\right]$

其中，表示由第m个阵元的z轴、x轴及y轴方向质点振速传感器输出信号的第n次采样数据叠加构成第m个阵元的第n次采样的四元数数据，表示由第m个阵元的x轴、y轴z轴方向质点振速传感器接收信号延时ΔT后的第n次采样数据叠加构成第m个阵元延时ΔT后的第n次采样的四元数数据；

第一组采样数据X构成的四元数数据矩阵Z₁＝B₁P+N₁，其中，P＝[p₁，p₂，...，p_K]为K个互不相关入射声波的声强，N₁为四元数噪声，为M×K的导向矢量矩阵，A₁＝[A₁₁，...，A_1k，...，A_1K]表示z轴方向的振速矢量构成的子阵导向矢量，A₂＝[A₂₁，...，A_2k，...，A_2K]表示y轴和x轴方向的振速矢量构成的子阵导向矢量，其中 a_1k＝cosθ， 是M个传感器与位于坐标原点处的传感器之间相位差构成的空域导向矢量，θ_k∈[0，π/2]和φ_k∈[-π，π]是第k个声源信号的俯仰角和方位角，表示第m个传感器的角坐标；

第二组采样数据Y构成的四元数数据矩阵Z₂＝B₂P+N₂，其中，B₂是第二组采样数据对应的阵列导向矢量矩阵，N₂为四元数噪声，两组数据的导向矢量间满足：B₂＝B₁Ф，为时延矩阵；

步骤二、构造全阵列接收数据矩阵Z＝BP+N，其中，Z＝[Z₁；Z₂]，B＝[B₁；B₂]是全阵列接收数据对应的阵列导向矢量矩阵，N＝[N₁；N₂]是全阵列接收数据四元数噪声矩阵；计算全阵列接收四元数数据矩阵Z的自相关矩阵R_z，对数据自相关矩阵进行四元数特征分解，得到信号子空间E_s，按照阵列导向矢量矩阵相同的分块模式对信号子空间进行分块，通过矩阵块运算以及伪伴随矩阵特征分解得到阵列导向矢量估计值及全阵列导向矢量矩阵估计值

$R_z=\frac{1}{N}\&lsqb;{\mathrm{ZZ}}^H\&rsqb;={\mathrm{BR}}_P{B}^H+{\mathrm{\&sigma;}}^{2}I,$

其中，为入射信号的自相关函数，σ²为噪声的方差，I为2M×2M单位矩阵，(·)^H表示转置复共轭操作；

对数据自相关矩阵R_z进行四元数矩阵特征分解，得到信号子空间E_s，信号子空间与阵列导向矢量的关系为：E_s＝BT，其中T为K×K的非奇异矩阵，分别取E_s的前M行及后M行组成信号子空间矩阵E₁＝B₁T和E₂＝B₂T＝B₁ФT，则有其中

将矩阵转换为复数自伴随矩阵其中Ψ_r和Ψ_i的分别是Ψ的实部和虚部，和分别是Ψ_i和Ψ_r的共轭矩阵，是复数的虚数单位，对C_Ψ进行特征分解其中，V_k是特征矢量，U_k是特征值，Ψ的右特征值也是U_k，对应的特征矢量其中，Ф_M＝[I_M，-iI_M]是维数为M×2M的还原矩阵，I_M代表维数为M×M的单位阵；从而由第二次特征分解得到阵列导向矢量的估计值和延时ΔT后的阵列导向矢量估计值及全阵列导向矢量估计值其中，为Ψ的特征向量矩阵；

步骤三、重构子阵导向矢量矩阵的估计和根据步骤一中四元数的构造模式、四元数阵列导向矢量B₁的结构特点以及阵列导向矢量B₁和子阵导向矢量矩阵A₁和A₂关系得到两个子阵导向矢量矩阵的估计和通过矩阵运算得到和的旋转不变关系矩阵Ω，Ω是关于俯仰角和方位角解耦的对角矩阵，由Ω得到信号方位角和俯仰角的粗略估计值，进而得到方向余弦的粗略估计值，该粗略估计值将用于后续的稀疏阵列空域导向矢量矩阵解相位模糊；

根据导向矢量计算子阵导向矢量和根据步骤一中阵列导向矢量B₁和子阵导向矢量矩阵A₁和A₂的关系知：且与满足的关系，Ω＝diag[Ω₁₁，...，Ω_KK]是子阵导向矢量和间的旋转不变关系矩阵，是方位角和俯仰角的解耦函数，由Ω_kk可以很容易地获取声波到达角粗略的估计值和由到达角的粗略估计值可得到方向余弦的粗略而无模糊的估计值

步骤四、获取归一化空域导向矢量的估计值导向矢量矩阵的估计值的每一列对该列的第一个元素进行归一化处理得到归一化空域导向矢量的估计值根据步骤三中得到的方向余弦粗略估计值解空域导向矢量的相位周期性模糊，从而得到信号方位角和俯仰角的精确估计值；

根据导向矢量矩阵计算归一化空域导向矢量其中，表示的第k列，表示的第k列的第1个元素；对空域导向矢量取相位得到相位矩阵arg(·)表示取相位，是有模糊的相位估计值；

由方向余弦的粗略估计值得到稀疏空域导向矢量相位矩阵的粗略估计值通过解这个优化问题得到模糊数是r_opt，为阵列空域导向矢量的精确相位估计值，根据计算第k个声波信号的x轴方向的方向余弦的估计值和y轴方向的方向余弦估计值其中，[W]^#＝[(W)^HW]^-1(W)^H是位置矩阵W的伪逆矩阵，r_opt是模糊数矢量，是方向余弦的精确估计值，位置矩阵W为

$W=\frac{2\mathrm{\&pi;}R}{\mathrm{\&lambda;}}\left[\begin{array}{cc}0& 1\\ \sin \mathrm{\&Delta;}& \cos \mathrm{\&Delta;}\\ .& .\\ .& .\\ .& .\\ sin\&lsqb;(M-1)\mathrm{\&Delta;}\&rsqb;& cos\&lsqb;(M-1)\mathrm{\&Delta;}\&rsqb;\end{array}\right]$

其中，Δ＝2π/M。

根据方向余弦的精确估计值和得到信号到达角的精确估计值：

$\left\{\begin{array}{c}{\widetilde{\mathrm{\&theta;}}}_k=\arcsin \left(\sqrt{{\widetilde{\mathrm{\&alpha;}}}_k^2+{\widetilde{\mathrm{\&beta;}}}_k^2}\right)\\ \left\{\begin{array}{cc}{\widetilde{\mathrm{\&phi;}}}_k=\arctan \left(\frac{{\widetilde{\mathrm{\&alpha;}}}_k}{{\widetilde{\mathrm{\&beta;}}}_k}\right)& {\widetilde{\mathrm{\&beta;}}}_k\&GreaterEqual;0\\ {\widetilde{\mathrm{\&phi;}}}_k=\mathrm{\&pi;}+\arctan \left(\frac{{\widetilde{\mathrm{\&alpha;}}}_k}{{\widetilde{\mathrm{\&beta;}}}_k}\right)& {\widetilde{\mathrm{\&beta;}}}_k<0\end{array}\right.\end{array}\right.;$

前述步骤中的m＝1，...，M，M为阵列的阵元数，n＝1，...，N，N为采样次数， 为四元数的3个虚数单位。