CN110109053A - 一种未知声速环境下快速doa估计方法 - Google Patents
一种未知声速环境下快速doa估计方法 Download PDFInfo
- Publication number
- CN110109053A CN110109053A CN201910259714.7A CN201910259714A CN110109053A CN 110109053 A CN110109053 A CN 110109053A CN 201910259714 A CN201910259714 A CN 201910259714A CN 110109053 A CN110109053 A CN 110109053A
- Authority
- CN
- China
- Prior art keywords
- matrix
- axis
- signal
- formula
- array
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
- G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明公开了一种未知声速环境下快速DOA估计方法,步骤如下:采用任意交叉线阵结构进行一维DOA估计;对互协方差矩阵R进行奇异值分解,其奇异值从大到小排列,前K个奇异值对应的左奇异向量张成的空间为信号子空间,这K个奇异向量组成矩阵Us,将Us划分为矩阵U1和矩阵U2;求出目标信号在x轴上的波达方向角θxi以及在y轴上的波达方向角θyi的表达式;将目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式;根据两条线阵的几何关系推导出θxi和θyi的关系式,再最终求解出目标信号的波达方向角。该发明在得到声速无关的DOA估计表达式的同时解决了DOA估计处理过程中的参数配对问题,减小了DOA估计处理的复杂度。
Description
技术领域
本发明涉及目标定位技术领域,具体涉及一种基于任意夹角的二维均匀线阵在未知声速环境进行的快速水下一维波达方向估计的方法。
背景技术
波达方向估计(DOA估计)在众多领域已得到广泛应用,而一维水下DOA估计就是指在水面放置传感器阵列利用阵列信号处理技术来对水下声源目标进行一维波达方向估计的方法。与空间目标DOA估计有所不同的是,水下声源目标的DOA估计采用声波作为传播载体,由于声波信号在水下环境传播时,水声信道中的各种障碍物及崎岖不平的海底造成的声波散射作用,导致了信号的衰减。除了水声环境造成信号的快速衰减,水下DOA估计面临的另一个问题就是声速影响。由于河流和海洋等水下环境复杂且不稳定,声波的速度随位置和时间而变化,水下DOA算法的估计精度受到很大影响,当实际声速偏离预先设定速度,估计精度将因此降低。
现有研究已提出了利用不同构型的二维阵列构造二维角度关系,对声速因子进行消除的一维水下目标DOA估计方法。这类方法提高了未知声速环境下的水下目标DOA估计精度,如专利申请201811241541.8和201811338421.X等。但是在多目标声源的情况下,由于这类方法采用了全组合遍历对多组参数进行配对,使得算法复杂度非常高,装置的估计实时性得不到保证。如何在消除声速影响的同时降低算法的算法复杂度成了一个急需解决的问题。
发明内容
本发明的目的是为了解决现有技术中的上述缺陷,提供一种未知声速环境下快速DOA估计方法,通过对任意夹角的二维均匀线阵的接收信号进行处理,在DOA波达方向估计中消除声速这个因子,从而消除水下声速不确定性对目标定位精度的影响。同时由于采用了快速配对方法,算法复杂度得到很大的下降,有利于在实际测量中快速估计目标。同时在估计过程中采用双线阵接收信号的互协方差矩阵进行处理,本发明无需对多条线阵的接收信号进行独立估计再进行配对,有效减小了算法复杂度。
本发明的目的可以通过采取如下技术方案达到:
一种基于未知声速环境进行的快速水下一维波达方向估计的方法,所述的估计方法步骤如下:
S1、用于一维DOA估计的任意交叉线阵结构如图3所示,与常规一维DOA估计方法采用单条均匀线阵不同,本文所提的方法采用两条均匀线阵,每条线阵上各自分布有M个阵元,两条线阵接口处有一个公共阵元,阵元间距为d。两条线阵之间的夹角设为δ且建立坐标系,设一条线阵所在直线为x轴,另一条线阵所在直线设为y轴。由于水下DOA估计大部分情况是在水面位置对水面下方的目标进行定位,因此本文只考虑半个平面空间内目标源信号的情况,即x轴的上半平面空间。目标信号满足窄带条件,即当信号延迟远小于带宽倒数时,延迟作用相当于使基带信号产生一个相移。假设目标信号的个数为K,目标信号的中心频率为fi,i=1,2,...,K,各个信号入射路径的声速定义为ci,i=1,2,...,K,目标信号在x轴上的波达方向角设为θxi,i=1,2,...,K,在y轴上的波达方向角设为θyi,i=1,2,...,K。目标信号的待估计的波达方向角设为θi且θi=θxi,i=1,2,...,K。第i个信号与x轴线阵的夹角为αi且αi∈[0,π]。当有K个远场窄带相互独立的信号入射到图3所示阵列,x轴和y轴阵列接收到的信号可以写成如下的矢量形式:
X(t)=AxS+Nx (1)
Y(t)=AyS+Ny (2)
其中S是一个K×1维的源信号矩阵,另外Nx和Ny则是M×1维的噪声矩阵,Ax和Ay分别为x轴和y轴阵列的M×K阶方向矩阵,写成矢量形式有:
其中a(θxi)和a(θyi)分别是x轴和y轴阵列的第i个声源的导向向量,i=1,2,...,K。
对于两条线阵各自接收的信号X(t)和Y(t),求它们的互协方差矩阵Rxy=E{X(t)YH(t)},然后将互协方差矩阵Rxy划分为子矩阵Rxy,1和Rxy,2,最后将子矩阵Rxy,1和Rxy,2组合成新的互协方差矩阵R。
S2、对矩阵R进行奇异值分解,其奇异值从大到小排列,前K个奇异值对应的左奇异向量张成的空间为信号子空间,这K个奇异向量组成矩阵Us。将Us划分为矩阵U1和矩阵U2。
S3、对矩阵进行特征值分解,得到K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK,并由这K个特征向量组成的矩阵估计出非奇异阵T,根据K个特征值求出目标信号在x轴上的波达方向角θxi的表达式,以及求出目标信号在y轴上的波达方向角θyi的表达式。
S4、将步骤S3得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式(33)。
S5、根据两阵列的几何关系推导出θxi和θyi的关系式,再最终求解出目标信号的波达方向角。
进一步地,所述的步骤S1中对接收信号矩阵进行处理得到新的互协方差矩阵R的过程如下:
根据线阵夹角δ和αi将目标源信号入射区域分成4部分:当入射信号与x轴的夹角αi∈[0,δ)时,信号入射区域设为区域①,时为区域②,时为区域③,时为区域④。
设入射信号与x轴法线的夹角为xni且与y轴法线的夹角为yni且
以原点处的阵元为参考阵元,对于x轴阵列,当入射信号最先到达参考阵元时,则波达方向角θxi为正值,此时波达方向角等于信号与阵列法线的夹角,即θxi=xni。否则当信号最后到达参考阵元时,波达方向角θxi为负值,且θxi=-xni。对于y轴阵列也有相同的结论。因此有:
sinθyi=sin(θxi-δ) (5)
得到X(t)和Y(t)的M×M阶互协方差矩阵为:
由于噪声为零均值白噪声,各噪声之间相互统计独立且都与目标信号统计独立,所以式(6)的后三项都为零,公式(6)改写为:
其中Rs=E{S(t)SH(t)}为信源部分的互协方差矩阵,由各信源是统计独立的,且Rs是对角矩阵。
取Ay的前(M-1)行和后(M-1)行分别设为Ay1和Ay2,即:
由Ay的表达式推出:
Ay2=Ay1ΩH (9)
其中
接下来将矩阵Rxy的前M-1列划分为子阵Rxy,1,后M-1列划分为子阵Rxy,2,即:
结合式(9)进一步可以得出:
将Rxy,1和Rxy,2按下式组合成新的2M×(M-1)阶矩阵,将此矩阵定义为新的互协方差矩阵R,
进一步地,所述的步骤S2中对矩阵R进行奇异值分解和排列,得到矩阵U1和矩阵U2的过程如下:
根据式(12)定义新的方向矩阵为:
则式(12)可以重写为:
对R进行奇异值分解:
其中U和V分别表示R的左右奇异矢量,且它们都是酋矩阵。Σ为对角矩阵,其对角线元素表示奇异值且它们从大到小排列,在无噪声的情况下,前K个奇异值大于零,剩余的奇异值为零。在有噪声影响的情况下,前K个奇异值也远远大于剩余奇异值。因此可以用前K个奇异值对应的奇异矢量构成信号子空间Us,剩余(M-1-K)个奇异值对应的奇异矢量构成噪声子空间Un。Σs和Σn对角线上的元素分别表示信号子空间和噪声子空间的奇异值,Vs和Vn分别为信号子空间和噪声子空间的右奇异矢量。
由奇异值分解性质可知:
RVn=ΣnUn (16)
两边都取共轭转置并右乘Un得到:
由式(15)和式(17)可知,由于噪声之间的独立性,矩阵R不受噪声的影响,因此噪声子空间的奇异值都为0,即Σn=0。因此式(17)改写为:
Vn HRHUn=0 (18)
又Vn为矩阵R奇异值分解得到的奇异值向量矩阵,所以Vn是满秩矩阵,由满秩矩阵性质可解出式(18)得:RHUn=0,并结合式(14)代入式(18)可得:
由于是满秩的,因此:
BHUn=0 (20)
由于U=[Us,Un]为酋矩阵,所以有:
结合式(20)和(21)可知,方向矩阵B与矩阵R的信号特征向量组成的子矩阵Us所张成的子空间相同,因此存在一个非奇异矩阵T使得:
Us=BT (22)
由式(13)可以得知方向矩阵B为2M×K阶矩阵,由两个M×K的矩阵Ax和AxΩ组成,易知Ax具有范德蒙德结构。定义两个子阵B1和B2,B1由Ax的前M-1行和AxΩ的前M-1行组成,B2由Ax的后M-1行和AxΩ的后M-1行组成,即:
观察式(13)和(23)很显然可以得到:
B2=B1Ψ (24)
其中Ψ=diag{ψ1,ψ2,...,ψK}且Ψ叫做旋转矩阵其对角线元素为相位旋转算子,易知通过求解旋转算子ψi就可以求出θxi。将Us按照相同的方法划分为U1和U2,即:
由Us=BT可得到:
进一步地,所述的步骤S3中对矩阵进行特征值分解,进行处理后求出目标信号进行处理后求出目标信号在在x轴上的波达方向角θxi的表达式和y轴上的波达方向角θyi的表达式的过程如下:
由式(24)和(26)可以推出:
U2=B1ΨT=B1TT-1ΨT=U1T-1ΨT (27)
其中T-1为T的逆矩阵,定义为矩阵U1的Moore-Penrose广义逆,公式(27)两边左乘有:
则矩阵为矩阵Ψ的相似变换,所以它们有相同的特征值,而Ψ由于是对角矩阵因此其特征值就是其对角线元素所以可以求解出的K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK。最后可以求出θxi的表达式为:
观察公式(28)可知,矩阵T-1由的K个特征向量α1,α2,...,αK组成,则可得到估计值
由(22)可以得出方向矩阵B的估计值为:
观察B的展开式(11),可以得到如下表达式:
其中表示矩阵的第k行第i列的元素。
进一步地,所述的步骤S4中对得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式的过程如下:
此时求出的结果还跟声速和频率有关,还需进一步消除声速和频率的影响。
用公式(32)和(29)相除可以得到:
进一步地,所述的步骤S5中根据两阵列的几何关系推导出θxi和θyi的关系式,最终求解出目标信号的波达方向角的过程如下:
由于只对矩阵进行了一次特征值分解,由特征值得到上式的分母,特征值对应的特征向量矩阵得到上式的分子,因此无需进行参数配对。再结合由两线阵之间的几何关系得出的公式(5)可以求出θxi的表达式如下:
根据第i个信号的波达方向角θi就是其在x轴上的波达方向角θxi,得到最终的一维波达方向估计结果为:
本发明相对于现有技术具有如下的优点及效果:
1、本发明公开的DOA估计方法不仅仅适用于L型阵列,还适用于其他任意夹角的二维均匀线阵,在阵列构型的选择上具有更大的灵活空间。
2、本发明与利用传统的一维DOA算法进行水下目标波达方向估计的方法相比更具有实用性,估计精确度也更高。传统的一维DOA算法通常假定声速为一个常量,而在实际的复杂水下环境中,声速往往是不断变化的,如果把其当成一个常量来进行计算的话,会导致较大的误差。本发明通过阵列夹角与波达方向角之间的关系消去了声速这个变量,使得最后的运算结果与声速无关,从而提高了估计精度。
3、本发明方法的时间复杂度主要集中在3个部分:一次互协方差矩阵构造、一次奇异值分解和一次特征值分解,其它部分的计算复杂度可以忽略。设采样快拍数为L,则两个线阵的接收信号X(t)和Y(t)的阶数都是M×L,因此构造互协方差矩阵的复杂度为O(M 2L)。矩阵R的阶数为2M×(M-1),因此对它进行奇异值分解的复杂度为O((M-1)3)。矩阵的阶数为K×K,对其进行特征值分解的复杂度为O(K 3),因此总的复杂度为
O(M 2L+(M-1)3+K3);
相较于其他未知声速环境的DOA估计方法,本发明提出的DOA估计方法时间复杂度更低,而且当目标信源数越多时,本发明所提方法的优势越明显,有利于目标的实时DOA估计。
4、本发明方法的实施装置在传统的测量装置上进行了改进,使用夹角可调节的均匀线阵,可行性强,安装简单。除此之外,现代处理器计算处理能力的不断提高,这使得本发明所使用的处理器等芯片的集成度高,并且计算能力强,从而保证了本发明的可行性。
附图说明
图1是本发明装置的硬件结构模块图;
图2是本发明中阵列的接收阵元与处理器连接示意图;
图3是本发明所用的任意夹角二维均匀线阵模型示意图;
图4是子线阵3的旋转连接示意图;
图5是x轴均匀线阵的接收信号模型示意图;
图6是本发明公开的未知声速环境下快速DOA估计方法的流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例一
本发明采用任意夹角的二维线阵,窄带目标声源为S,中心频率为f。声波入射方向,即方位角和仰角可表示为θk和
如附图6所示,本实施例中基于水下DOA估计方法包括以下步骤:
S1、用于一维DOA估计的任意交叉线阵结构如图3所示,与常规一维DOA估计方法采用单条均匀线阵不同,本文所提的方法采用两条均匀线阵,每条线阵上各自分布有M个阵元,两线阵接口处有一个公共阵元,阵元间距为d。两条线阵之间的夹角设为δ且建立坐标系,设一条线阵所在直线为x轴,另一条线阵所在直线设为y轴。由于水下DOA估计大部分情况是在水面位置对水面下方的目标进行定位,因此本文只考虑半个平面空间内目标源信号的情况,即x轴的上半平面空间。目标信号满足窄带条件,即当信号延迟远小于带宽倒数时,延迟作用相当于使基带信号产生一个相移。当有K个远场窄带相互独立的信号入射到图3所示阵列,x轴和y轴阵列接收到的信号可以写成如下的矢量形式:
X(t)=AxS+Nx (1)
Y(t)=AyS+Ny (2)
其中S是一个K×1维的源信号矩阵,另外Nx和Ny则是M×1维的噪声矩阵,Ax和Ay分别为x轴和y轴阵列的M×K阶方向矩阵,写成矢量形式有:
其中a(θxi)和a(θyi)分别是x轴和y轴阵列的第i个声源的导向向量,i=1,2,...,K。
假设目标信号的个数为K,目标信号的中心频率为fi,i=1,2,...,K,各个信号入射路径的声速定义为ci,i=1,2,...,K,目标信号在x轴上的波达方向角设为θxi,i=1,2,...,K,在y轴上的波达方向角设为θyi,i=1,2,...,K。目标信号的待估计的波达方向角设为θi且θi=θxi,i=1,2,...,K。第i个信号与x轴线阵的夹角为αi且αi∈[0,π]。根据线阵夹角δ和αi将目标源信号入射区域分成4部分:当入射信号与x轴的夹角αi∈[0,δ)时,信号入射区域设为区域①,时为区域②,时为区域③,时为区域④。
设入射信号与x轴法线的夹角为xni且与y轴法线的夹角为yni且
当信号从区域①入射时,有则有:
以原点处的阵元为参考阵元,对于x轴阵列,当入射信号最先到达参考阵元时,则波达方向角θxi为正值,此时波达方向角等于信号与阵列法线的夹角,即θxi=xni。否则当信号最后到达参考阵元时,波达方向角θxi为负值,且θxi=-xni。对于y轴阵列也有相同的结论。因此有:
sinθyi=sin(θxi-δ) (5)
X(t)和Y(t)的M×M阶互协方差矩阵为:
由于噪声为零均值高斯白噪声,各噪声之间相互统计独立且都与目标信号统计独立所以式(6)的后三项都为零,公式(6)改写为:
其中Rs=E{S(t)SH(t)}为信源部分的互协方差矩阵,由各信源是统计独立的,且Rs是对角矩阵。
取Ay的前(M-1)行和后(M-1)行分别设为Ay1和Ay2,即:
由Ay的表达式推出:
Ay2=Ay1ΩH (9)
其中
接下来将矩阵Rxy的前M-1列划分为子阵Rxy,1,后M-1列划分为子阵Rxy,2,即:
结合式(9)进一步可以得出:
将Rxy,1和Rxy,2按下式组合成新的2M×(M-1)阶矩阵,将此矩阵定义为新的互协方差矩阵R,
S2、对互协方差矩阵R进行奇异值分解,其奇异值从大到小排列,前K个奇异值对应的左奇异向量张成的空间为信号子空间,这K个奇异向量组成矩阵Us。将Us按照式(36)的方式划分为矩阵U1和矩阵U2。根据式(12)定义新的方向矩阵为:
则式(12)可以重写为:
对互协方差矩阵R进行奇异值分解:
其中U和V分别表示R的左右奇异矢量,且它们都是酋矩阵。Σ为对角矩阵,其对角线元素表示奇异值且它们从大到小排列,在无噪声的情况下,前K个奇异值大于零,剩余的奇异值为零。在有噪声影响的情况下,前K个奇异值也远远大于剩余奇异值。因此可以用前K个奇异值对应的奇异矢量构成信号子空间U s,剩余(M-1-K)个奇异值对应的奇异矢量构成噪声子空间Un。Σs和Σn对角线上的元素分别表示信号子空间和噪声子空间的奇异值,Vs和Vn分别为信号子空间和噪声子空间的右奇异矢量。
由奇异值分解性质可知:
RVn=ΣnUn (16)
两边都取共轭转置并右乘Un得到:
由式(6)和式(7)可知,由于噪声之间的独立性,矩阵R不受噪声的影响,因此噪声子空间的奇异值都为0,即Σn=0。因此式(17)改写为:
Vn HRHUn=0 (18)
又Vn为互协方差矩阵R奇异值分解得到的奇异值向量矩阵,所以Vn是满秩矩阵,由满秩矩阵性质可解出式(18)得:
RHUn=0 (19)
结合式(13)代入式(19)可得:
由于是满秩的,因此BHUn=0。由于U=[Us,Un]为酋矩阵,所以有:
结合式(20)和(21)可知,方向矩阵B与矩阵R的信号特征向量组成的子矩阵Us所张成的子空间相同,因此存在一个非奇异矩阵T使得:
Us=BT (22)
由式(13)可以得知方向矩阵B为2M×K阶矩阵,由两个M×K的矩阵Ax和AxΩ组成,易知Ax具有范德蒙德结构。定义两个子阵B1和B2,B1由Ax的前M-1行和AxΩ的前M-1行组成,B2由Ax的后M-1行和AxΩ的后M-1行组成,即:
观察式(22)和(23)很显然可以得到:
B2=B1Ψ (24)
其中Ψ=diag{ψ1,ψ2,...,ψK}且Ψ叫做旋转矩阵其对角线元素为相位旋转算子,易知通过求解旋转算子ψi就可以求出θxi。将Us按照相同的方法划分为U1和U2,即:
由Us=BT可得到:
S3、对矩阵进行特征值分解,得到K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK,并由这K个特征向量组成的矩阵估计出非奇异阵T,根据K个特征值求出目标信号在x轴上的波达方向角θxi的表达式,求出目标信号在y轴上的波达方向角θyi的表达式。
由式(24)和(26)可以推出:
U2=B1ΨT=B1TT-1ΨT=U1T-1ΨT (27)
定义为矩阵U1的Moore-Penrose广义逆,公式(27)两边左乘有:
则矩阵为矩阵Ψ的相似变换,所以它们有相同的特征值,而Ψ由于是对角矩阵因此其特征值就是其对角线元素所以可以求解出的K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK。最后可以求出θxi的表达式为:
此时求出的结果还跟声速和频率有关,还需进一步消除声速和频率的影响。观察公式(28)易知矩阵T-1由的K个特征向量α1,α2,...,αK组成,则得到估计值
由(22)可以得出方向矩阵B的估计值为:
观察B的展开式(13),可以得到如下表达式:
其中表示矩阵的第k行第i列的元素。
S4、将步骤S3得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相比得到一个与声速无关的表达式。用公式(32)和(29)相除可以得到:
S5、根据两阵列的几何关系推导出θxi和θyi的关系式,再结合式(33)最终求解出目标信号的波达方向角。由于只对矩阵进行了一次特征值分解,由特征值得到上式的分母,特征值对应的特征向量矩阵得到上式的分子,因此无需进行参数配对。再结合由两线阵之间的几何关系得出的公式(5)可以求出θxi的表达式如下:
因为第i个信号的波达方向角θi就是其在x轴上的波达方向角θxi,所以最终的一维波达方向估计结果为:
实施例二
本实例公开了一种基于未知声速环境的快速水下波达方向估计装置,所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块,如图1和图2所示。
数据处理与控制模块由一对A/D、D/A转换器和一个处理器组成,是整个装置的核心部分,其它所有模块都与它直接相连。它可以控制发射模块,使发射模块发射我们指定的信号;可以控制接收模块的夹角可调线阵,使均匀线阵1保持固定,均匀线阵2以连接点为中心进行自由旋转,可转至设定值;同时能够对接收模块传过来的信号进行处理,通过本发明的算法计算出波达方向角,然后将结果传输至输出模块。
接收模块包括2个以均匀间距摆放的超声波探头阵列,步进电机和步进电机驱动电路。步进电机是将电脉冲信号转变为角位移或线位移的开环控制电机,当步进电机驱动电路收到一个脉冲信号,它就驱动步进电机按设定的方向转动固定的角度,称为步距角。所以可以通过使数据处理与控制模块发射一定数量的脉冲信号来达到期望的角度值。如图3所示,使均匀线阵1排布于坐标系x轴上保持固定,因为接收模块会放置在水中,所以固定支架采用塑料材质以增大浮力。均匀线阵2安装到步进电机上,可由步进电机带动旋转,从而达到夹角调节的目的,图4为均匀线阵2与步进电机的连接旋转示意图,如图所示,步进电机通过旋转转子与均匀线阵2连接在一起,以控制均匀线阵2的旋转。
发射模块由一个阻抗匹配电路和一个超声波发射探头组成,通过D/A转换器与处理器相连,能够根据处理器发出的指令发射指定的信号。
输出模块由一个USB接口和一个显示器组成,并且与数据处理与控制模块和电源模块相连。它能够提供人机交互,将数据处理与控制模块中处理好的数据通过USB接口输出到外部装置或者在显示器上显示出来。
电源模块由一个电源组成,并且与数据处理与控制模块、发射模块、接收模块和输出模块相连。它能够为这些模块供电。
本发明装置的主要工作流程如下:在实测过程中根据我们想要发射的信号参数,通过数据处理与控制模块输入对应的参数,使处理器产生相应的数字信号,然后通过D/A转换后传给发射模块,超声波发射探头就能产生我们需要的信号并进行发射。均匀线阵1和均匀线阵2之间的夹角值δ可以通过数据处理与控制模块进行设定,处理器发送特定的脉冲信号到步进电机驱动电路,然后驱动步进电机转动至我们需要的角度。接收模块中的接收阵列收到从目标声源反射回来的信号后将其通过A/D转换成数字信号后发送给处理器,然后处理器根据本发明提供的算法计算出结果。最后数据处理与控制模块将计算结果传给输出模块,输出模块将结果通过USB接口传给外部设备或者通过显示器显示出来。电源模块为所有其它模块供电。
实施例三
本实例公开了一种基于未知声速环境的快速水下波达方向估计装置,所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块,如图1和图2所示。
数据处理与控制模块可以用DSP芯片实现(如:TI公司TMS320VC5509A型号的DSP芯片),此DSP芯片可实现A/D转换和D/A转换的功能,并能够实现三维均匀线阵的旋转算子和最终波达方向的计算;
接收模块中的步进电机采用东芝公司的23HY6606-CP型号的电机,此步进电机的步距角为1.8度,步进电机驱动电路采用东芝公司的TC78S600FTG型芯片。此外接收模块还使用1个固定均匀直线阵列和1自由旋转的均匀直线阵列,其中每个阵列包括多个超声接收探头,并且数量相同,2个均匀阵列按图3所示组装;发射模块使用一个超声波发射探头;输出模块使用一个USB接口和一个LCD显示屏。图1即为本发明所述装置的硬件结构模块图。
本发明的主要工作步骤具体如下:
步骤T1、按图2连接好具体装置,其中接收模块中的每个均匀线阵中的阵元个数M统一定为8。利用数据处理与控制模块发送指令,控制超声发射探头发射超声信号s(t);海水中声速范围大致为1430m/s-1550m/s,则取最小声速为1430m/s,可以求出最小半波长为7.15cm。所以设置两个均匀线阵的平均间距为5cm,即第一个阵元和最后一个阵元相隔35cm。任意两相邻线阵之间的距离必须小于7.15cm,在满足此限制条件下可以任意选取阵元间距,这里2个均匀线阵的间距都取为4cm。线阵1与线阵2之间的夹角设置线阵夹角值为72°。在数据处理与控制模块设定线阵夹角值,首先将均匀线阵夹角δ转为72°。在水下放置K=3个目标声源,目标信号的中心频率分别为f=10,15,20kHz,各个信号入射路径的声速分别为c=1450,1500,1550m/s,入射的波达方向角分别为(40°,60°,70°)。
步骤T2、对超声接收探头线阵接收到的目标声源信号进行采样;均匀线阵1接收到的信号为x1(t),x2(t),…,x8(t),均匀线阵2接收的信号为y1(t),y2(t),…,y8(t),共采样接收200次,并将接收到的信号传递给数据处理与控制模块进行分析处理。
步骤T3、信号在处理模块中的分析处理步骤具体如下:
1)对于两条线阵各自接收的信号X(t)和Y(t),求它们的互协方差矩阵Rxy=E{X(t)YH(t)},然后将互协方差矩阵Rxy划分为子矩阵Rxy,1和Rxy,2,最后将Rxy,1和Rxy,2按组合成新的矩阵R。
2)对矩阵R进行奇异值分解,其奇异值从大到小排列,前K个奇异值对应的左奇异向量张成的空间为信号子空间,这K个奇异向量组成矩阵Us。将Us划分为矩阵U1和矩阵U2。
3)对矩阵进行特征值分解,得到K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK,并由这K个特征向量组成的矩阵估计出非奇异阵T,根据K个特征值按式(29)求出目标信号在x轴上的波达方向角θxi的表达式,按式(32)求出目标信号在y轴上的波达方向角θyi的表达式。
4)将以上步骤3)得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相比得到一个与声速无关的表达式(33)。
5)根据两阵列的几何关系推导出θxi和θyi的关系式,再结合式(33)最终求解出目标信号的波达方向角。
根据本发明中快速DOA估计方法,估计出的二维波达方向角(40.15°,59.89°,70.23°),对目标估计达到了预期精度,说明估计结果正确,本发明方法及装置可行。
综上所述,上述实施例提出一种基于任意交叉线阵的无需配对的声速无关快速一维DOA估计方法与装置,利用对交叉线阵的接收信号进行联合处理,在得到声速无关的DOA估计表达式的同时解决了DOA估计处理过程中的参数配对问题,减小了DOA估计处理的复杂度。
上述实施例为本发明较佳的实施方式,但本发明的实施方式并不受上述实施例的限制,其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化,均应为等效的置换方式,都包含在本发明的保护范围之内。
Claims (6)
1.一种未知声速环境下快速DOA估计方法,其特征在于,所述的估计方法包括以下步骤:
S1、采用任意交叉线阵结构进行一维DOA估计,该任意交叉线阵结构采用两条均匀线阵,每条线阵上各自分布有M个阵元,两条线阵接口处有一个公共阵元,阵元间距为d,两条线阵之间的夹角设为δ且建立坐标系,设一条线阵所在直线为x轴,另一条线阵所在直线设为y轴,假设只考虑半个平面空间内目标源信号的情况,即x轴的上半平面空间,假设目标信号满足窄带条件,即当信号延迟远小于带宽倒数时,延迟作用相当于使基带信号产生一个相移。假设目标信号的个数为K,目标信号的中心频率为fi,i=1,2,...,K,各个信号入射路径的声速定义为ci,i=1,2,...,K,目标信号在x轴上的波达方向角设为θxi,i=1,2,...,K,在y轴上的波达方向角设为θyi,i=1,2,...,K,目标信号的待估计的波达方向角设为θi且θi=θxi,i=1,2,...,K,第i个信号与x轴线阵的夹角为αi且αi∈[0,π],当有K个远场窄带相互独立的信号入射到上述任意交叉线阵结构,x轴和y轴阵列接收到的信号写成如下的矢量形式:
X(t)=AxS+Nx (1)
Y(t)=AyS+Ny (2)
其中S是一个K×1维的源信号矩阵,Nx和Ny是M×1维的噪声矩阵,Ax和Ay分别为x轴和y轴阵列的M×K阶方向矩阵,写成矢量形式有:
其中a(θxi)和a(θyi)分别是x轴和y轴阵列的第i个声源的导向向量,i=1,2,...,K,
对于两条线阵各自接收的信号X(t)和Y(t),求它们的互协方差矩阵Rxy=E{X(t)YH(t)},然后将互协方差矩阵Rxy划分为子矩阵Rxy,1和Rxy,2,最后将子矩阵Rxy,1和Rxy,2组合成新的互协方差矩阵R;
S2、对互协方差矩阵R进行奇异值分解,其奇异值从大到小排列,前K个奇异值对应的左奇异向量张成的空间为信号子空间,这K个奇异向量组成矩阵Us,将Us划分为矩阵U1和矩阵U2;
S3、对矩阵U1 +U2进行特征值分解,得到K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK,并由这K个特征向量组成的矩阵估计出非奇异阵T,根据K个特征值求出目标信号在x轴上的波达方向角θxi的表达式,以及求出目标信号在y轴上的波达方向角θyi的表达式;
S4、将步骤S3得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式;
S5、根据两条线阵的几何关系推导出θxi和θyi的关系式,再最终求解出目标信号的波达方向角。
2.根据权利要求1所述的一种未知声速环境下快速DOA估计方法,其特征在于,所述的步骤S1中对接收信号矩阵进行处理得到新的互协方差矩阵R的过程如下:
根据线阵夹角δ和αi将目标源信号入射区域分成4部分:当入射信号与x轴的夹角αi∈[0,δ)时,信号入射区域设为区域①,时为区域②,时为区域③,时为区域④;
设入射信号与x轴法线的夹角为xni且与y轴法线的夹角为yni且
以原点处的阵元为参考阵元,对于x轴阵列,当入射信号最先到达参考阵元时,则波达方向角θxi为正值,此时波达方向角等于信号与阵列法线的夹角,即θxi=xni,否则当信号最后到达参考阵元时,波达方向角θxi为负值,且θxi=-xni,对于y轴阵列也有相同的结论,因此有:
sinθyi=sin(θxi-δ) (5)
得到X(t)和Y(t)的M×M阶互协方差矩阵为:
由于噪声为零均值白噪声,各噪声之间相互统计独立且都与目标信号统计独立,将式(6)的后三项都为零,式(6)改写为:
其中Rs=E{S(t)SH(t)}为信源部分的互协方差矩阵,由各信源是统计独立的,且Rs是对角矩阵;
取Ay的前(M-1)行和后(M-1)行分别设为Ay1和Ay2,即:
由Ay的表达式推出:
Ay2=Ay1ΩH (9)
其中
接下来将矩阵Rxy的前M-1列划分为子阵Rxy,1,后M-1列划分为子阵Rxy,2,即:
结合式(9)得出:
将Rxy,1和Rxy,2按下式组合成新的2M×(M-1)阶矩阵,将此矩阵定义为新的互协方差矩阵R,
3.根据权利要求2所述的一种未知声速环境下快速DOA估计方法,其特征在于,所述的步骤S2中对互协方差矩阵R进行奇异值分解和排列,得到矩阵U1和矩阵U2的过程如下:
根据式(12)重构出新的方向矩阵B为:
则式(12)重写为:
对R进行奇异值分解:
其中U和V分别表示R的左右奇异矢量,且它们都是酋矩阵,Σ为对角矩阵,其对角线元素表示奇异值且它们从大到小排列,在无噪声的情况下,前K个奇异值大于零,剩余的奇异值为零,在有噪声影响的情况下,前K个奇异值也远远大于剩余奇异值,因此用前K个奇异值对应的奇异矢量构成信号子空间Us,剩余(M-1-K)个奇异值对应的奇异矢量构成噪声子空间Un,Σs和Σn对角线上的元素分别表示信号子空间和噪声子空间的奇异值,Vs和Vn分别为信号子空间和噪声子空间的右奇异矢量;
由奇异值分解性质可知:
RVn=ΣnUn (16)
两边都取共轭转置并右乘Un得到:
由式(15)和式(17)可知,由于噪声之间的独立性,互协方差矩阵R不受噪声的影响,因此噪声子空间的奇异值都为0,即Σn=0,式(17)改写为:
Vn HRHUn=0 (18)
又Vn为矩阵R奇异值分解得到的奇异值向量矩阵,所以Vn是满秩矩阵,由满秩矩阵性质可解出式(18)得:RHUn=0,并结合式(14)代入式(18)可得:
由于是满秩的,因此:
BHUn=0 (20)
由于U=[Us,Un]为酋矩阵,所以有:
结合式(20)和(21)可知,方向矩阵B与互协方差矩阵R的信号特征向量组成的子矩阵Us所张成的子空间相同,因此存在一个非奇异矩阵T使得:
Us=BT (22)
由式(13)可以得知方向矩阵B为2M×K阶矩阵,由两个M×K的矩阵Ax和AxΩ组成,得出Ax具有范德蒙德结构,定义两个子阵B1和B2,B1由Ax的前M-1行和AxΩ的前M-1行组成,B2由Ax的后M-1行和AxΩ的后M-1行组成,即:
由式(13)和(23)得到:
B2=B1Ψ (24)
其中Ψ=diag{ψ1,ψ2,...,ψK}且Ψ叫做旋转矩阵其对角线元素为相位旋转算子,通过求解旋转算子ψi求出θxi,将Us按照相同的方法划分为U1和U2,即:
由Us=BT可得到:
4.根据权利要求3所述的一种未知声速环境下快速DOA估计方法,其特征在于,所述的步骤S3中对矩阵U1 +U2进行特征值分解,进行处理后求出目标信号在在x轴上的波达方向角θxi的表达式和y轴上的波达方向角θyi的表达式的过程如下:
由式(24)和(26)推出:
U2=B1ΨT=B1TT-1ΨT=U1T-1ΨT (27)
其中T-1为T的逆矩阵,定义U1 +为矩阵U1的Moore-Penrose广义逆,公式(27)两边左乘U1 +有:
U1 +U2=T-1ΨT (28)
则矩阵U1 +U2为矩阵Ψ的相似变换,所以它们有相同的特征值,而Ψ由于是对角矩阵因此其特征值就是其对角线元素求解出U1 +U2的K个特征值λ1,λ2,...,λK和对应的特征向量α1,α2,...,αK,最后求出θxi的表达式为:
根据式(28)可知,矩阵T-1由U1 +U2的K个特征向量α1,α2,...,αK组成,则可得到估计值
由(22)可以得出方向矩阵B的估计值为:
由方向矩阵B的展开式(11)得到如下表达式:
其中表示矩阵的第k行第i列的元素。
5.根据权利要求4所述的一种未知声速环境下快速DOA估计方法,其特征在于,所述的步骤S4中对得到的目标信号在y轴上的波达方向角表达式和在x轴上的波达方向角表达式相除得到一个与声速无关的表达式的过程如下:
用公式(32)和公式(29)相除得到:
6.根据权利要求5所述的一种未知声速环境下快速DOA估计方法,其特征在于,所述的步骤S5中根据两条线阵的几何关系推导出θxi和θyi的关系式,最终求解出目标信号的波达方向角的过程如下:
通过结合由两条线阵之间的几何关系得出的公式(5)求出θxi的表达式如下:
根据第i个信号的波达方向角θi就是其在x轴上的波达方向角θxi,得到最终的一维波达方向估计结果为:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910259714.7A CN110109053B (zh) | 2019-04-02 | 2019-04-02 | 一种未知声速环境下快速doa估计方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910259714.7A CN110109053B (zh) | 2019-04-02 | 2019-04-02 | 一种未知声速环境下快速doa估计方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110109053A true CN110109053A (zh) | 2019-08-09 |
CN110109053B CN110109053B (zh) | 2021-01-19 |
Family
ID=67484906
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910259714.7A Expired - Fee Related CN110109053B (zh) | 2019-04-02 | 2019-04-02 | 一种未知声速环境下快速doa估计方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110109053B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112327244A (zh) * | 2020-10-22 | 2021-02-05 | 中国电子科技集团公司第五十四研究所 | 一种基于l型阵列的二维非相干分布式目标参数估计方法 |
CN113504505A (zh) * | 2021-06-02 | 2021-10-15 | 华南理工大学 | 一种适用于低信噪比环境下的一维doa估计方法 |
CN113504504A (zh) * | 2021-06-04 | 2021-10-15 | 华南理工大学 | 一种水下高精度一维doa估计方法 |
CN113640737A (zh) * | 2021-07-27 | 2021-11-12 | 哈尔滨工程大学 | 一种基于二维功率分布的少阵元阵列高分辨方位估计方法 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102608565A (zh) * | 2012-03-23 | 2012-07-25 | 哈尔滨工程大学 | 一种基于均匀圆阵列的波达方向估计方法 |
CN101431354B (zh) * | 2007-11-09 | 2013-03-27 | 中兴通讯股份有限公司 | 一种波达角估计方法 |
CN106154219A (zh) * | 2015-04-22 | 2016-11-23 | 常熟海量声学设备科技有限公司 | 一种新的水声目标方位估计方法 |
US20170090016A1 (en) * | 2015-09-25 | 2017-03-30 | Texas Instruments Incorporated | Method for Joint Antenna-Array Calibration and Direction of Arrival Estimation for Automotive Applications |
CN108594166A (zh) * | 2018-04-19 | 2018-09-28 | 广东工业大学 | 一种二维波达方向估计方法及装置 |
US20190049579A1 (en) * | 2017-08-09 | 2019-02-14 | Sony Corporation | PERFORMANCE OF A TIME OF FLIGHT (ToF) LASER RANGE FINDING SYSTEM USING ACOUSTIC-BASED DIRECTION OF ARRIVAL (DoA) |
CN109521392A (zh) * | 2018-10-24 | 2019-03-26 | 华南理工大学 | 基于非圆信号和l型线阵的水下一维doa估计方法和装置 |
-
2019
- 2019-04-02 CN CN201910259714.7A patent/CN110109053B/zh not_active Expired - Fee Related
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101431354B (zh) * | 2007-11-09 | 2013-03-27 | 中兴通讯股份有限公司 | 一种波达角估计方法 |
CN102608565A (zh) * | 2012-03-23 | 2012-07-25 | 哈尔滨工程大学 | 一种基于均匀圆阵列的波达方向估计方法 |
CN106154219A (zh) * | 2015-04-22 | 2016-11-23 | 常熟海量声学设备科技有限公司 | 一种新的水声目标方位估计方法 |
US20170090016A1 (en) * | 2015-09-25 | 2017-03-30 | Texas Instruments Incorporated | Method for Joint Antenna-Array Calibration and Direction of Arrival Estimation for Automotive Applications |
US20190049579A1 (en) * | 2017-08-09 | 2019-02-14 | Sony Corporation | PERFORMANCE OF A TIME OF FLIGHT (ToF) LASER RANGE FINDING SYSTEM USING ACOUSTIC-BASED DIRECTION OF ARRIVAL (DoA) |
CN108594166A (zh) * | 2018-04-19 | 2018-09-28 | 广东工业大学 | 一种二维波达方向估计方法及装置 |
CN109521392A (zh) * | 2018-10-24 | 2019-03-26 | 华南理工大学 | 基于非圆信号和l型线阵的水下一维doa估计方法和装置 |
Non-Patent Citations (1)
Title |
---|
LIANG ZHANG: "Direction-of-arrival estimation for far-field", 《ELECTRONICS LETTERS》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112327244A (zh) * | 2020-10-22 | 2021-02-05 | 中国电子科技集团公司第五十四研究所 | 一种基于l型阵列的二维非相干分布式目标参数估计方法 |
CN112327244B (zh) * | 2020-10-22 | 2022-06-24 | 中国电子科技集团公司第五十四研究所 | 一种基于l型阵列的二维非相干分布式目标参数估计方法 |
CN113504505A (zh) * | 2021-06-02 | 2021-10-15 | 华南理工大学 | 一种适用于低信噪比环境下的一维doa估计方法 |
CN113504505B (zh) * | 2021-06-02 | 2023-11-03 | 华南理工大学 | 一种适用于低信噪比环境下的一维doa估计方法 |
CN113504504A (zh) * | 2021-06-04 | 2021-10-15 | 华南理工大学 | 一种水下高精度一维doa估计方法 |
CN113504504B (zh) * | 2021-06-04 | 2023-06-20 | 华南理工大学 | 一种水下高精度一维doa估计方法 |
CN113640737A (zh) * | 2021-07-27 | 2021-11-12 | 哈尔滨工程大学 | 一种基于二维功率分布的少阵元阵列高分辨方位估计方法 |
CN113640737B (zh) * | 2021-07-27 | 2022-06-21 | 哈尔滨工程大学 | 一种基于二维功率分布的少阵元阵列高分辨方位估计方法 |
Also Published As
Publication number | Publication date |
---|---|
CN110109053B (zh) | 2021-01-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110109053A (zh) | 一种未知声速环境下快速doa估计方法 | |
Hawkes et al. | Effects of sensor placement on acoustic vector-sensor array performance | |
CN105607033B (zh) | 基于正交均匀线阵的水下波达方向估计方法及系统 | |
CN109581275B (zh) | 基于非圆信号和三维正交阵的二维水下doa估计方法和装置 | |
CN108008348B (zh) | 基于可调夹角均匀线阵的水下波达方向估计方法及装置 | |
USRE45823E1 (en) | System and method of acoustic doppler beamforming | |
CN108535682B (zh) | 一种基于旋转非均匀双l阵的水下二维doa估计方法及装置 | |
CN107942284B (zh) | 基于二维正交非均匀线阵的水下波达方向估计方法与装置 | |
CN108414967A (zh) | 基于夹角可调双l阵的水下二维波达方向估计方法与装置 | |
CN104931929B (zh) | 基于线阵综合声速补偿的近场波达方向估计方法及装置 | |
CN106249196A (zh) | 三分量声矢量传感器稀疏阵列四元数解模糊方法 | |
CN109407048B (zh) | 基于非圆信号和夹角可调阵的水下doa估计方法与装置 | |
CN109521392B (zh) | 基于非圆信号和l型线阵的水下一维doa估计方法和装置 | |
CN109884580A (zh) | 水下一维doa估计方法和装置 | |
CN109581274B (zh) | 基于夹角可调三维阵的非圆信号水下doa估计方法和装置 | |
CN112098938B (zh) | 一种基于六元锥矢量阵的水声目标降维匹配声场定位方法 | |
CN208000373U (zh) | 基于夹角可调双l阵的水下二维波达方向估计装置 | |
CN111722232A (zh) | 一种具备三维定位能力的多波束成像声纳实时信号处理装置 | |
CN209514042U (zh) | 基于非圆信号和夹角可调阵的水下doa估计装置 | |
CN209640474U (zh) | 基于夹角可调三维阵的非圆信号水下doa估计装置 | |
CN209640473U (zh) | 基于非圆信号和三维正交阵的二维水下doa估计装置 | |
CN110824484B (zh) | 一种基于恒模算法的阵元位置估计方法 | |
CN113504504B (zh) | 一种水下高精度一维doa估计方法 | |
Jin et al. | A Single-snapshot Compressed Sensing Technique for Source Localization in Underwater Transportation | |
CN114142900B (zh) | 一种基于lcmv准则的多通道自适应波束形成方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210119 |
|
CF01 | Termination of patent right due to non-payment of annual fee |