CN106249225A - 稀疏圆形声矢量传感器阵列四元数esprit参数估计方法 - Google Patents

Abstract

稀疏圆形声矢量传感器阵列四元数ESPRIT参数估计方法，阵列接收K个不同频率的远场、互不相关窄带声信号，阵列接收的两组复数域数据都按照同阵元的同一次采样数据的声压及三个方向相互垂直的振速数据叠加构成四元数一次采样数据；对全阵列接收数据自相关矩阵特征分解及相关运算得到阵列导向矢量的估计值；由阵列导向矢量估计值重构子阵导向矢量并由子阵导向矢量间的旋转不变关系得到信号方向余弦的粗略估计值；利用方向余弦粗略估计值确定空域导向矢量的相位周期模糊数，利用最小二乘法得到信号到达角的精确估计值；本发明保持了声矢量传感器各分量间的正交特性，具有更好的模型误差鲁棒性，解决了稀疏声矢量传感器圆形阵列的相位模糊问题。

Description

稀疏圆形声矢量传感器阵列四元数ESPRIT参数估计方法

技术领域

本发明属于信号处理技术领域，尤其涉及声矢量传感器阵列四元数ESPRIT到达角估计方法。

背景技术

声矢量传感器是一种新型的声源信号测向设备，声矢量传感器技术近年来受到国内外许多学者的关注，声矢量传感器由传统无指向性的声压传感器与具有偶极子指向性的质点振速传感器构成，声矢量传感器能够同步、共点测量声场中某一点的声压强度和质点振速，由于声矢量传感器具有的偶极子指向性与频率无关，使得它更适合用于低频信号的测量与检测；

声矢量传感器信号处理是阵列信号处理的一个新领域，在声学工程特别是水声工程和空间声学工程中得到了广泛的应用。声矢量传感器技术在水声工程的应用领域有水声警戒声纳、拖曳线列阵声纳、舷侧共形阵声纳、水雷声引信、鱼雷探测声纳、多基地声纳、水下潜器的导航定位、分布式传感器网络等；声矢量传感器技术在空间声学工程的应用领域有战场警戒探测直升机和隐形飞机、声强测量、声功率测量及噪声源识别等；

声矢量传感器输出的各个分量按照一定顺序构成一个长的复数矢量，这就是声矢量传感器的长矢量模型。长矢量模型只是简单地将声矢量传感器的输出各分量排列成一个长矢量，没有考虑声矢量传感器各输出分量之间的正交关系，因此，长矢量模型不能充分体现声矢量传感器技术的优越性。四元数的正交结构更适合描述声矢量传感器各个输出分量间固有的正交性，研究基于四元数理论的声矢量传感器阵列信号处理不仅具有重要的理论价值，而且具有重要的工程应用价值；

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

发明内容

针对上述技术的不足和缺陷，本发明提出了一种基于四元数的稀疏圆形声矢量传感器阵列波达方向估计方法，该方法增加了布阵的灵活性，且通过相位解模糊处理提高了到达角的估计精度。

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

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

声矢量传感器阵列的四元数ESPRIT到达角估计方法，包括以下步骤：

K个不同频率的远场、互不相关窄带声源信号同时入射到声矢量传感器均匀圆形阵列上，所述阵列的阵元是由一个声压传感器和三个空间轴向垂直的振速传感器组成的声矢量传感器，

步骤一、对M个声矢量传感器组成的圆形阵列接收信号进行N次采样得到第一组采样数据X，阵列的接收信号延时ΔT后进行同步采样，得到第二组采样数据Y，X和Y均为4M×N的矩阵，K＜M-1，采样数据X和Y的各个组成分量叠加构成四元数数据Z₁和Z₂；

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

将第一组采样数据X和第二组采样数据Y都按照同阵元的声压强度、x轴方向质点振速传感器、y轴方向质点振速传感器和z轴方向质点振速传感器的同次快拍数据叠加构成第一组接收四元数数据矩阵Z₁和第二组接收四元数数据矩阵Z₂：

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

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

第二组采样数据Y构成的四元数数据矩阵Z₂＝B₂P+N₂，其中，B₂＝B₁Φ第二组采样数据对应的阵列导向矢量矩阵，为时延矩阵，N₂为四元数表示的噪声，信号和噪声互不相关；

构造全阵列接收数据矩阵其中，是全阵列接收数据对应的阵列导向矢量矩阵，是全阵列接收数据四元数噪声矩阵；

步骤二、计算全阵列接收数据矩阵Z的数据自相关矩阵R_z，对数据自相关矩阵进行四元数特征分解，得到第一组四元数数据对应的阵列导向矢量估计值第二组四元数数据对应的阵列导向矢量估计值和全阵列接收数据对应的阵列导向矢量矩阵估计值

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

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

对数据自相关矩阵R_z进行四元数矩阵特征分解，得到由K个大特征值对应的特征向量构成的信号子空间E_s和2M-K个小特征值对应的特征向量构成的噪声子空间E_n，根据子空间原理，存在K×K的非奇异矩阵T，满足E_s＝BT，分别取E_s的前M行及后M行组成信号子空间矩阵E₁和E₂，由信号子空间的定义，B₁、B₂与E₁、E₂之间满足E₁＝B₁T，E₂＝B₂T＝B₁ΦT，则有其中

矩阵不是厄米特矩阵不能直接进行四元数的特征分解，需要转换为复数自伴随矩阵其中Ψ＝Ψ_r+iΨ_i，i是复数域虚数单位矢量，Ψ_r和Ψ_i分别是Ψ的实部和虚部，和分别是Ψ_i和Ψ_r的共轭矩阵，对C_Ψ进行特征分解其中，c_k是特征矢量，v_k是特征值，Ψ的右特征值也是v_k，对应的特征矢量其中，是维数为J×2J的还原矩阵，J＝M，I_J代表维数为J×J的单位阵；从而由第二次特征分解得到阵列导向矢量的估计值和延时ΔT后的阵列导向矢量估计值及全阵列导向矢量估计值以及时延矩阵的估计值其中，为Ψ的特征向量矩阵；

步骤三、由时延矩阵估计值估计信号频率，由导向矢量重构子阵导向矢量和根据两子阵导向矢量之间的旋转不变关系矩阵Ω得到信号方位角和俯仰角的粗略估计值，进而得到方向余弦的粗略估计值；

由时延矩阵估计值可得到频率的估计值

${\hat{f}}_k=\frac{\arg \left({\widehat{\mathrm{\Φ}}}_{kk}\right)}{2\mathrm{\π}\mathrm{\Δ}T}$

由导向矢量重构子阵导向矢量和即：且与满足的关系，Ω＝diag[Ω₁₁，…，Ω_KK]，由关系式可获得声波到达角粗略估计值和其中，arg(·)表示取相位；

由到达角的粗略估计值可得到方向余弦的粗略而无模糊的估计值

步骤四、由导向矢量矩阵的估计值得到归一化空域导向矢量的估计值根据方向余弦的粗略估计值解空域导向矢量的相位周期性模糊，从而得到信号方位角和俯仰角的精确估计值；

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

通过这个优化问题得到模糊数矢量r_opt，其中是由步骤三中方向余弦的粗略估计值得到的粗略相位矩阵估计值，则精确相位矩阵估计值为根据计算第k个声波信号的x轴方向和y轴方向余弦的精确估计值和其中，[W]^#＝[(W)^H W]^-1(W)^H是位置矩阵W的伪逆矩阵，位置矩阵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{\θ}}}_k=arcsin\left(\sqrt{{\widetilde{\mathrm{\α}}}_k^2+{\widetilde{\mathrm{\β}}}_k^2}\right)\\ {\widetilde{\mathrm{\φ}}}_k=\arctan \frac{{\widetilde{\mathrm{\α}}}_k}{{\widetilde{\mathrm{\β}}}_k}\end{array}\right.;$

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

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

四元数方法能够更好地保持声矢量传感器的矢量特性，具有比长矢量方法更好的性能，耦合误差更小。本发明方法采用基于四元数的ESPRIT算法联合利用子阵间的旋转不变关系和空域导向矢量估计入射信号的多参数，与只利用空域导向矢量的非稀疏阵四元数方法相比，可以更好的体现声矢量传感器各组成分量的正交特性，提高参数估计的精度，增加了布阵的灵活性。

附图说明

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

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

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

图3为旋转不变算法到达角估计散布图。

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

图5为仿真实验的俯仰角估计标准偏差与信噪比的关系图。

图6为仿真实验的方位角估计标准偏差与信噪比的关系图。

图7为仿真实验的俯仰角估计成功概率随信噪比的关系图。

图8为仿真实验的方位角估计成功概率随信噪比的关系图。

具体实施方式

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

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

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

步骤一、对M个声矢量传感器组成的圆形阵列接收信号进行N次采样得到第一组采样数据X，阵列的接收信号延时ΔT后进行同步采样，得到第二组采样数据Y，X和Y均为4M×N的矩阵，K＜M-1，采样数据X和Y的各个组成分量叠加构成四元数数据Z₁和Z₂；

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

将第一组采样数据X和第二组采样数据Y都按照同阵元的声压强度、x轴方向质点振速传感器、y轴方向质点振速传感器和z轴方向质点振速传感器的同次快拍数据叠加构成第一组接收四元数数据矩阵Z₁和第二组接收四元数数据矩阵Z₂：

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

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

第二组采样数据Y构成的四元数数据矩阵Z₂＝B₂P+N₂，其中，B₂＝B₁Φ第二组采样数据对应的阵列导向矢量矩阵，为时延矩阵，N₂为四元数表示的噪声，信号和噪声互不相关；

构造全阵列接收数据矩阵其中，是全阵列接收数据对应的阵列导向矢量矩阵，是全阵列接收数据四元数噪声矩阵；

步骤二、计算全阵列接收数据矩阵Z的数据自相关矩阵R_z，对数据自相关矩阵进行四元数特征分解，得到第一组四元数数据对应的阵列导向矢量估计值第二组四元数数据对应的阵列导向矢量估计值和全阵列接收数据对应的阵列导向矢量矩阵估计值

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

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

对数据自相关矩阵R_z进行四元数矩阵特征分解，得到由K个大特征值对应的特征向量构成的信号子空间E_s和2M-K个小特征值对应的特征向量构成的噪声子空间E_n，根据子空间原理，存在K×K的非奇异矩阵T，满足E_s＝BT，分别取E_s的前M行及后M行组成信号子空间矩阵E₁和E₂，由信号子空间的定义，B₁、B₂与E₁、E₂之间满足E₁＝B₁T，E₂＝B₂T＝B₁ΦT，则有其中

矩阵不是厄米特矩阵不能直接进行四元数的特征分解，需要转换为复数自伴随矩阵其中Ψ＝Ψ_r+iΨ_i，i是复数域虚数单位矢量，Ψ_r和Ψ_i分别是Ψ的实部和虚部，和分别是Ψ_i和Ψ_r的共轭矩阵，对C_Ψ进行特征分解其中，c_k是特征矢量，v_k是特征值，Ψ的右特征值也是v_k，对应的特征矢量其中，是维数为J×2J的还原矩阵，J＝M，I_J代表维数为J×J的单位阵；从而由第二次特征分解得到阵列导向矢量的估计值和延时ΔT后的阵列导向矢量估计值及全阵列导向矢量估计值以及时延矩阵的估计值其中，为Ψ的特征向量矩阵；

步骤三、由时延矩阵估计值估计信号频率，由导向矢量重构子阵导向矢量和根据两子阵导向矢量之间的旋转不变关系矩阵Ω得到信号方位角和俯仰角的粗略估计值，进而得到方向余弦的粗略估计值；

由时延矩阵估计值可得到频率的估计值

${\hat{f}}_k=\frac{\arg \left({\widehat{\mathrm{\Φ}}}_{kk}\right)}{2\mathrm{\π}\mathrm{\Δ}T}$

由导向矢量重构子阵导向矢量和即：且与满足的关系，Ω＝diag[Ω₁₁，…，Ω_KK]，由关系式可获得声波到达角粗略估计值和其中，arg(·)表示取相位；

由到达角的粗略估计值可得到方向余弦的粗略而无模糊的估计值

步骤四、由导向矢量矩阵的估计值得到归一化空域导向矢量的估计值根据方向余弦的粗略估计值解空域导向矢量的相位周期性模糊，从而得到信号方位角和俯仰角的精确估计值；

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

通过这个优化问题得到模糊数矢量r_opt，其中是由步骤三中方向余弦的粗略估计值得到的粗略相位矩阵估计值，则精确相位矩阵估计值为根据计算第k个声波信号的x轴方向和y轴方向余弦的精确估计值和其中，[W]^#＝[(W)^HW]^-1(W)^H是位置矩阵W的伪逆矩阵，位置矩阵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{\θ}}}_k=arcsin\left(\sqrt{{\widetilde{\mathrm{\α}}}_k^2+{\widetilde{\mathrm{\β}}}_k^2}\right)\\ {\widetilde{\mathrm{\φ}}}_k=\arctan \frac{{\widetilde{\mathrm{\α}}}_k}{{\widetilde{\mathrm{\β}}}_k}\end{array}\right.;$

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

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

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

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

现有技术的旋转不变算法就是利用四元数子阵间的旋转不变关系进行俯仰角和方位角的估计；仿真实验采用旋转不变算法与本发明方法的解模糊算法进行对比，如图3至图8所示；

图3和图4为信噪比是10dB时，旋转不变算法和解模糊算法到达角估计的散布图，从图3至图4可以看出旋转不变算法的俯仰角和方位角在真值附近的很大范围内扰动，而本发明方法的俯仰角和方位角在真值附近较小范围内扰动，旋转不变算法的到达角估计精度低于本发明的解模糊算法的到达角估计精度的原因在于本发明方法通过稀疏布阵，增加了阵列孔径，并通过解模糊处理提高了到达角的估计精度；

从图5和图6可以看出，在信噪比区间内本发明方法的标准偏差都低于旋转不变算法的标准偏差，也就是估计值在均值附近的更小范围内容扰动；在信噪比为0dB时，本发明方法在方位角和俯仰角估计标准偏差跟旋转不变算法相比分别小3.5°和0.5°；

当估计值与真值的偏差小于1度认为估计成功，本发明方法的俯仰角和方位角的估计成功概率高于旋转不变算法，特别是0dB时，从图7可以看出，本发明方法的俯仰角成功概率超过了95％，而旋转不变算法的成功概率低于60％；从图8可以看出，本发明方法的方位角成功概率达到了59％，而旋转不变算法的成功概率仅是17％；

因此从参数估计的散布图、标准偏差和成功概率上可以看出本发明方法比旋转不变方法具有更高的参数估计性能；旋转不变算法没有利用阵列孔径信息，虽然可以估计到达角，但精度不高，稀疏阵列排布扩大了阵列孔径，但引起了相位模糊，本发明算法首先利用旋转不变关系求得粗略估计值，然后利用粗略估计值解相位模糊，从而大大提高了参数估计精度，精度的提高是因为充分利用了阵列孔径信息和子阵间的旋转不变关系；

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

Claims (1)

1.稀疏圆形声矢量传感器阵列四元数ESPRIT参数估计方法，其特征在于：

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

稀疏圆形声矢量传感器阵列四元数ESPRIT参数估计方法的步骤如下：阵列接收K个不同频率的远场、互不相关的窄带声波入射信号，

步骤一、对M个声矢量传感器组成的圆形阵列接收信号进行N次采样得到第一组采样数据X，阵列的接收信号延时ΔT后进行同步采样，得到第二组采样数据Y，X和Y均为4M×N的矩阵，K＜M-1，采样数据X和Y的各个组成分量叠加构成四元数数据Z₁和Z₂；

$X=\left[\begin{array}{ccccc}{x}_1\left(1\right),& ...,& x_1\left(n\right),& ...,& x_1\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .\\ x_m\left(1\right),& ...,& x_m\left(n\right),& & x_m\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .& ...,& .,& ...,& .\\ x_M\left(1\right),& ...,& x_M\left(n\right),& ...,& x_M\left(N\right)\end{array}\right],Y=\left[\begin{array}{ccccc}{y}_1\left(1\right),& ...,& y_1\left(n\right),& ...,& y_1\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .,& ...,& .,& ...,& .\\ y_m\left(1\right),& ...,& y_m\left(n\right),& & y_m\left(N\right)\\ .& & .& & .\\ .& & .& & .\\ .& ...,& .,& ...,& .\\ y_M\left(1\right),& ...,& y_M\left(n\right),& ...,& y_M\left(N\right)\end{array}\right]$

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

将第一组采样数据X和第二组采样数据Y都按照同阵元的声压强度、x轴方向质点振速传感器、y轴方向质点振速传感器和z轴方向质点振速传感器的同次快拍数据叠加构成第一组接收四元数数据矩阵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个阵元的声压强度、x轴方向质点振速传感器、y轴方向质点振速传感器、z轴方向质点振速传感器输出信号的第n次采样数据叠加构成第m个阵元的第n次采样的四元数数据，表示由第m个阵元的声压强度、x轴方向质点振速传感器、y轴方向质点振速传感器、z轴方向质点振速传感器接收信号延时ΔT后的第n次采样数据叠加构成第m个阵元延时ΔT后的第n次采样的四元数数据；

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

第二组采样数据Y构成的四元数数据矩阵Z₂＝B₂P+N₂，其中，B₂＝B₁Ф第二组采样数据对应的阵列导向矢量矩阵，为时延矩阵，N₂为四元数表示的噪声，信号和噪声互不相关；

构造全阵列接收数据矩阵其中，是全阵列接收数据对应的阵列导向矢量矩阵，是全阵列接收数据四元数噪声矩阵；

步骤二、计算全阵列接收数据矩阵Z的数据自相关矩阵R_z，对数据自相关矩阵进行四元数特征分解，得到第一组四元数数据对应的阵列导向矢量估计值第二组四元数数据对应的阵列导向矢量估计值和全阵列接收数据对应的阵列导向矢量矩阵估计值

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

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

对数据自相关矩阵R_z进行四元数矩阵特征分解，得到由K个大特征值对应的特征向量构成的信号子空间E_s和2M-K个小特征值对应的特征向量构成的噪声子空间E_n，根据子空间原理，存在K×K的非奇异矩阵T，满足E_s＝BT，分别取E_s的前M行及后M行组成信号子空间矩阵E₁和E₂，由信号子空间的定义，B₁、B₂与E₁、E₂之间满足E₁＝B₁T，E₂＝B₂T＝B₁ФT，则有其中

矩阵不是厄米特矩阵不能直接进行四元数的特征分解，需要转换为复数自伴随矩阵其中Ψ＝Ψ_r+iΨ_i，i是复数域虚数单位矢量，Ψ_r和Ψ_i分别是Ψ的实部和虚部，和分别是Ψ_i和Ψ_r的共轭矩阵，对C_Ψ进行特征分解其中，c_k是特征矢量，v_k是特征值，Ψ的右特征值也是v_k，对应的特征矢量其中，是维数为J×2J的还原矩阵，J＝M，I_J代表维数为J×J的单位阵；从而由第二次特征分解得到阵列导向矢量的估计值和延时ΔT后的阵列导向矢量估计值及全阵列导向矢量估计值以及时延矩阵的估计值其中，为Ψ的特征向量矩阵；

步骤三、由时延矩阵估计值估计信号频率，由导向矢量重构子阵导向矢量和根据两子阵导向矢量之间的旋转不变关系矩阵Ω得到信号方位角和俯仰角的粗略估计值，进而得到方向余弦的粗略估计值；

由时延矩阵估计值可得到频率的估计值

${\hat{f}}_k=\frac{\arg \left({\widehat{\mathrm{\Φ}}}_{kk}\right)}{2\mathrm{\π}\mathrm{\Δ}T}$

由导向矢量重构子阵导向矢量和即：且与满足的关系，Ω＝diag[Ω₁₁，...，Ω_KK]，由关系式可获得声波到达角粗略估计值和其中，arg(·)表示取相位；

由到达角的粗略估计值可得到方向余弦的粗略而无模糊的估计值

步骤四、由导向矢量矩阵的估计值得到归一化空域导向矢量的估计值根据方向余弦的粗略估计值解空域导向矢量的相位周期性模糊，从而得到信号方位角和俯仰角的精确估计值；

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

通过这个优化问题得到模糊数矢量r_opt，其中是由步骤三中方向余弦的粗略估计值得到的粗略相位矩阵估计值，则精确相位矩阵估计值为根据计算第k个声波信号的x轴方向和y轴方向余弦的精确估计值和其中，[W]^#＝[(W)^HW]^-1(W)^H是位置矩阵W的伪逆矩阵，位置矩阵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{\θ}}}_k=\arcsin \left(\sqrt{{\widetilde{\mathrm{\α}}}_k^2+{\widetilde{\mathrm{\β}}}_k^2}\right)\\ {\widetilde{\mathrm{\φ}}}_k=\arctan \frac{{\widetilde{\mathrm{\α}}}_k}{{\widetilde{\mathrm{\β}}}_k}\end{array}\right.;$

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