中国指挥与控制学会会刊 
信号波达方向估计(Direction of Arrival, DOA)是阵列信号处理领域重要的研究内容之一,其目的是对传感器阵列接收的数据进行处理,从而估计出目标方位。信号波达方向估计技术在现代军事和国民经济等领域均发挥了重要作用,并具有良好的应用前景。
DOA估计是声呐信号处理的基础。目前水下目标DOA估计的研究主要集中于均匀浅海声场环境,阵列信号建模以平面波声场模型为基础[1⇓⇓-4]。与浅海不同,典型的深海环境(深度大于2 000 m)[5]主要有5种声传播模式:表面波导模式、会聚区模式、海底反射模式、可靠声路径模式和表面波导泄漏模式。由于多途效应和深海声速分布不均匀等因素的影响,基于平面波声场模型的深海DOA估计存在误差,这一问题一直未得到很好解决。宫先仪[6]等人研究了多途信号、波导中的简正波等典型传播模型下的阵列信号处理方法。吴俊楠[7-8]根据简正波理论和射线理论分析表明在不同距离情况下,有贡献的简正波不同,与这类简正波对应的相速度也不同,参考声速选取不当会导致 DOA 估计出现误差。李辉、杨坤德[9-10]对拖曳线列阵在会聚区内对远场声源进行DOA估计时存在偏差进行了研究分析,发现方位估计误差随距离呈周期性变化,并与传播损失随距离的变化规律成正相关,从简正波和射线理论分析了深海环境下测向误差的成因,并提出一种稳健的测向误差修正方法。韩志斌[11]基于简正波理论和射线理论,指出了通过声线到达角或简正波相速度估计,合理设置波束形成权矢量,是消除水平接收阵测向误差的有效途径,并进一步提出采用线性拟合法快速预报目标所在距离的声线到达角,然后估计简正波相速度代替常规波束形成的参考声速,实现了测向误差的有效校正。
本文从射线声学理论出发,建立了多途信道下深海DOA估计的阵列信号模型,以及深海多径声线传播时间和相邻阵元声线传播时延差的表征方法,结合解相干DOA估计算法,提高了深海DOA估计性能,从根本上解决了深海DOA估计存在误差的问题,并通过仿真实验进行了验证。
1 阵列信号建模
1.1 基于平面波声场模型的阵列信号建模
1.1.1 窄带阵列信号模型[12]
设阵元数为M,K个窄带声源,对于远场平面波入射,第m个阵元的接收信号可以表示为
xm(t)= si(t-τmi)+nm(t),m=1,2…M
其中,si(t)为参考阵元接收到的第i个信源的信号,τmi是第i个信号到达第m个阵元时的相对于参考阵元的时延,nm(t)是第m个阵元的噪声。
对于窄带信号:
si(t-τmi)≈si(t) ,i=1,2,…K
其中,ω0=2πf=2π 为信号频率,c是声速,λ是波长。
将M个阵元在t时刻的接收信号排列成一个列矢量:
= +
将式(3)改写成矢量形式,得
X(t)=AS(t)+N(t)
其中,X(t)表示阵列M×1维快拍数据矢量,A是M×K维阵列流型矩阵,S(t)是空间信号的K×1维矢量,N(t)表达阵列M×1维噪声数据矢量,阵列流型矩阵A可以表示为
Α=[a1(ω0),a2(ω0),…,aK(ω0)]
其中,ai(ω0)= (i=1,2,…,K)为导向矢量,T表示转置,τmi是第i个信号到达第m个阵元时相对于参考阵元的时延。对于平面波声场中的线列阵,τm= ,其中,d为阵元间距,θ为入射信号方位。
1.1.2 相干信源的阵列信号模型[12]
对于随机的两个平稳入射信号si(t)和sk(t),它们间的相关系数定义为
γ=
根据Schwartz不等式,可知|γ|≤1,信号之间的相关性定义如下:
由定义式(7)可知,相干信号源之间只差一个复常数,即
si(t)=ais0(t),i=1,2,…,K
将公式(8)代入(4),可得到相干信号入射时的阵列信号模型:
X(t)=A s0(t)+N(t)=Aρs0(t)+N(t)
其中,ρ代表K×1维的复常数矢量。
以上为平面波声场模型下的阵列信号模型,平面波声场模型认为声速均匀分布,且只考虑直达波,目前的水下目标DOA估计大部分是基于平面波声场模型的,而实际的水下信道为多途信道,声速分布可能是非均匀的,如深海信道。
1.2 基于多途信道声场模型的阵列信号建模
深海的非均匀声速分布和多途传播信道使深海DOA估计不同于浅海,本文从射线声场模型出发,推导声线传播时间和相邻阵元声线传播时延差的表示式,并对多途信道的阵列信号进行建模,为深海DOA估计研究奠定基础。
1.2.1 射线声场模型与本征声线传播时延
根据射线声学理论[13],程函φ可以写成水平分量φ1(r)和垂直分量φ2(z)的线性叠加:
φ(r,z)=φ1(r)+φ2(z)=cosα0r+ dz+C
其中,a0是声线的初始掠射角,r,z分别是信源与接收阵列的水平距离与接收阵深度,n为折射率,C为常数。
平面问题的第l条射线声场表示式为[13]
pl(r,z)=Al(r,z)
其中,k0= 为波数,c(zs)为声源深度处的声速,Al为第l条声线的声压振幅,φl(r,z)表示第l条声线的程函。
根据叠加性原理,接收点的所有L条声线的总声压可表示为[13]
p(r,z)= pl(r,z)= Al(r,z)
其中,k0φl(r,z)即为第l条声线传播时延所产生的相位,到达第m个阵元的第l条声线的传播时延为
Δ = rm+ dz
其中, =k(zs)cos 为到达第m个阵元的第l条声线的水平波数, = 为第l条声线的垂直波数, 为到达第m个阵元的第l条声线的出射角,zs和zr分别是声源的深度与接收阵列深度,rm为信源到第m个阵元的水平距离。
根据式(13),相邻阵元间的传播时延差为
Δ = rm- rm+1+ dz
以均匀浅海声场环境为例,此时声速c=c0为常数,因此折射率n=n(z)=1,波数k= ,到达第m个阵元的第l条声线的水平波数kr,m=k(zs)cos = cos ,垂直波数kz,m= ,利用式(14)得到相邻阵元的传播时延差为
Δ =
上式对直达声线和任意多途声线均成立,不同的声线出射角不同。进一步分析表明,远场条件下,到达相邻阵元的第l条声线的出射角近似相等,式(15)可简化为
Δt=
式(16)即为常用的平面波声场条件下相邻阵元的时延差。对比式(14)和式(16),表明本文推导出的相邻阵元传播时延差表达式更具有普遍性,为解决深海DOA估计问题奠定了基础。
1.2.2 多途信道下的阵列信号建模
设阵元数为M,K个窄带声源,每个声源的发射信号经多径传播有Lk(k=1,2,…,K)条声线到达接收阵元,则第m个阵元的接收信号可以表示为
xm(t)= bk,lsk(t- )+nm(t)
以第一个阵元为参考阵元,令:
X=[x1(t),…,xM(t)]T,N=[n1(t),…,nM(t)]T
S=[s1(t),…,sK(t)]T,A=[A1(θ),…,AK(θ)]T
Ak(θ)=[Ak(θ1),…,Ak()],Ak(θl)=[1, ,…, ]T
B=diag(B1,…,BK),Bk=diag(bk,l,…, )
则M个阵元输出矢量为
X=ABS+N
其中, 表示到达第m个阵元的第k个信号源的第l条声线相对于参考阵元的时延差,可由式(14)计算得到。在进行DOA估计时,以此时延差结合式(20)构造导向矢量,可以提高DOA估计的性能,将在第四章中仿真验证。由此可见,相干多途信道下,如深海Munk声速分布信道,相邻阵元的时延差,除与阵元间距和目标方位有关外,还与声源深度和声源到各阵元的水平距离有关,DOA估计时导向向量的构造要考虑这些因素的影响,此时DOA估计实为多维参数优化问题。
2 DOA估计算法
2.1 MUSIC[14](Multiple Signal Classification,MUSIC)算法
MUSIC算法即多重信号分类算法,其基本思想是利用信号子空间和噪声子空间的正交性估计信号的波达方向,基本原理是对阵列接收数据的协方差矩阵进行特征值分解,得到信号子空间和噪声子空间,由于这两个子空间之间是相互正交的,可构造出“针状”空间谱从而突破瑞利限,提高角度分辨率。
对于平面波声场中的阵列信号模型,即式(4),其阵列数据的协方差为
R=E[XXH]=ARSAH+σ2I
其中,H表示矩阵的共轭转置,RS为信号的协方差矩阵,σ2为噪声功率。对R进行特征分解有
$\boldsymbol{R}=\boldsymbol{U}_{s} \Sigma_{s} \boldsymbol{U}_{s}^{\mathrm{H}}+\boldsymbol{U}_{N} \Sigma_{N} \boldsymbol{U}_{N}^{\mathrm{H}}$
式中,US对应信号子空间的特征矢量矩阵,$\sum$S是R的K个大特征值构成的对角阵,UN对应噪声子空间的特征矢量矩阵,$\sum$N是R的M-K个小特征值构成的对角阵。由于噪声子空间和信号子空间是相互正交的,有
a(θ)UN=0
MUSIC算法的基础即来自于式(25)。现实中接收数据长度有限,其协方差矩阵的最大似然估计为
= XXH
对其进行特征分解,得到噪声子空间的特征矢量矩阵 ,由于存在噪声,a(θ)与 并不完全正交,实际的方位估计是通过最小优化搜索实现的,即
θMUSIC= aH(θ) a(θ)
因此,MUSIC算法的空间谱估计为
PMUSIC=
传统的MUSIC算法只适用于窄带非相干源的DOA估计,在多途相干环境下,传统的MUSIC方法失效,因此本文采用一种重构子空间的改进MUSIC方法:RS-MUSIC,进行相干信道下的DOA估计。
设均匀线列阵有M个阵元,K个信源,由式(23)得到阵列数据协方差矩阵RX为
RX=E[X(t)XH(t)]
定义新的向量Y(t)=J(X(t))*,其中H表示矩阵的共轭转置,J是M×M阶反单位矩阵,因此Y(t)的协方差矩阵为
RY=E[Y(t)YH(t)]=J J
计算X(t)和Y(t)的互协方差矩阵,得
RXY=E[X(t)YH(t)]
计算Y(t)和X(t)的互协方差矩阵,得
RYX=E[Y(t)XH(t)]
根据式(29)、(30)、(31)、(32),利用X(t)和Y(t)各自的自协方差矩阵以及相互间的互协方差矩阵构造一个新的协方差矩阵:
RA=[RX RY RXY RYX]
RA是一个M×4M维的矩阵,对其进行奇异值分解得
RA=USVH
其中,U是一个M×M维的矩阵,S是一个M×4M维的矩阵,V是一个4M×4M维的矩阵,S=diag(λ1,λ2,…,λM)是特征值矩阵,且λ1≥λ2≥…≥λM。
根据式(29)(30)构建一个新的数据协方差矩阵R1:
R1= =
对其进行奇异值分解得
R1=U1S1(V1)H
其中,S1=diag(η1,η2,…,ηM)是特征值矩阵,且η1≥η2≥…≥ηM。
对RA和R1的特征值矩阵S和S1按照顺序比较对应的特征值,例如,比较λ1和η1,如果λ1>η1,则选择λ1作为新的特征值以构成新的特征值矩阵和新的特征向量矩阵,从新的特征向量矩阵中选择第K+1列到第M列重构噪声子空间Urecn,再利用式(28)得到RS-MUSIC算法的空间谱估计为
PRS-MUSIC=
3 仿真与结果分析
3.1 Bellhop声场模型的原理及使用方法
Bellhop声场模型是基于射线声学理论的信道模型,该模型通过高斯波束跟踪方法来计算水平非均匀环境中的声场。高斯波束跟踪方法的核心思想是将高斯强度的分布与水下声场中各条声线相联系,因此水下声场各条声线可用高斯声束的中心声线表示,这些声线能较平滑地穿过水下声场的焦散线区域,也能较平滑地过渡到水下声影区,有利于高频宽带信号在水平方向快速变化问题的解决,因此,基于Bellhop模型构建的水声信道与实际声场更一致。
Bellhop的基本原理是在射线中心坐标系求解波动方程,得到声场的声压解。
如图1 所示,单位正交向量t(s)、e1(s)、e2(s)构成一个新的坐标系(s,n1,n2),该坐标系称为射线中心坐标系。n1、n2分别表示沿e1(s)、e2(s)方向上的距离,e1(s)、e2(s)为垂直于声线r0(s)的平面上的两个单位正交向量,定义为:
=κ1(s)t(s), =κ2(s)t(s)
由波动方程及射线中心坐标系的拉普拉斯算子表达式,可以得到射线中心坐标系下的波动方程为
=
其中,h=1+n1κ1+n2κ2,(s,n1,n2)为射线中心坐标系下的坐标,c为声速,求解此波动方程即可得到声场的声压u。
3.1.2 Bellhop声场模型的使用方法[17]
在使用Bellhop模型时,首先需要根据实验环境及需求配置Bellhop输入文件,以环境文件.env为例,在配置.env文件时,需要根据场景需求在文件的对应行中,按一定格式设置输出文件;接着将配置好的输入文件如.env文件放在对应的信号处理MATLAB脚本文件(.m)的同一目录下,随后在MATLAB脚本文件中通过“bellhop xx.env”指令调用Bellhop模型,由Bellhop模型计算声场,得到描述声场特性的相应输出文件。具体的输入输出文件需要根据实际需要选取以下形式。
Bellhop 模型的输入文件:用于描述发送端与接收端之间的水声环境和几何特性。在最简单同时也是最典型的情况下,只有一个文件(*.env)。该文件作为环境文件,包含水声信道的声速分布以及有关海底海面情况的信息。如果海底特性随距离变化,还必须添加一个文件(*.bty)来确定海底形状。同理,当海面特性随距离变化时,需要添加一个文件(*.ati)来确定海面形状。此外,当声速分布随距离变化时,需要添加一个声速剖面分布文件(*.ssp),简称SSP(Sound Speed Profile, SSP)文件,来描述海洋声速的变化。最后,当需要指定一个由任意底部反射系数来表征海底时,可提供一个底部反射系数文件(.brc)来定义海底反射角度。同理,可通过提供顶部反射系数文件(*.trc)来定义海面的反射角度。
Bellhop 模型的输出文件:主要是通过Bellhop 模型程序对水声环境文件和相关选项的选择,产生不同的输出文件。当选择射线追踪选项时,就会产生一个包含从声源发出的一系列扇形声线的文件(*.ray);当选择本征射线选项时,扇形射线就只包含到达指定位置的接收阵元的声线;当选择幅度和传播时间选项时,则生成(*.arr)文件,该文件包含到达指定位置接收阵元的各条声线的时间延迟和振幅等信息,利用这些信息,可以计算出水声信道的冲击响应,再与发射信号进行卷积,即可得到阵元的接收信号。
3.2 相干多途信道下的空间谱
1)选择海洋环境为深海声场,Munk声速分布,海水深度H=5 000 m,海水密度ρ=1 024 kg/m3。海底声速1 600 m/s,海底密度是 1 024 kg/m3,吸收系数 0.35 dB/λ。采样频率为48 kHz,快拍数为4 096,快拍次数为4,均匀线列阵的阵元数为16,阵元间距取信号中心频率对应波长的一半,接收线列阵和信号源的深度均为300 m,信号源与均匀线列阵的首个阵元距离r=1 km,信噪比均为5 dB,设置直达波的入射方位角为60°,将这些参数输入海洋环境文件(*.env),利用 Bellhop 模型得到多途信道下的阵列接收信号。对于多途信道,阵列接收到的来自不同路径的信号是相干的;多途信道下的DOA估计,需要考虑相干源分辨的问题。考虑直达波和海面反射波,基于平面波声场阵列信号模型的RS-MUSIC 算法和基于本文多途信道阵列信号模型的RS-MUSIC 算法的空间谱如图2 所示。
其他条件不变,直达波角度设置为120°时,基于两种阵列信号模型的RS-MUSIC 算法的空间谱如图3 所示。
此外,图2 、图3 表明,与基于平面波声场模型的RS-MUSIC 算法相比,基于本文多途信道阵列信号模型的RS-MUSIC 算法的空间谱旁瓣更低,性能更好。
2)信号源与均匀线列阵的首个阵元距离改为10 km,其他条件同实验3.2(1)。直达波的入射方位角60°时基于平面波声场阵列信号模型的RS-MUSIC 算法和基于本文多途信道阵列信号模型的RS-MUSIC 算法的空间谱如图4 所示。
其他条件不变,直达波角度设置为120°时,基于两种阵列信号模型的RS-MUSIC 算法的空间谱如图5 所示。
3.3 相干多途信道下的DOA估计性能
实验条件同3.2(2),信噪比的变化范围是-40 dB到10 dB,分别对每组信噪比进行50次蒙特卡洛实验,深海Munk声速分布下的基于平面波声场阵列信号模型的RS-MUSIC 算法和基于本文多途信道阵列信号模型RS-MUSIC 算法的DOA估计偏差和均方根误差随信噪比变化的关系如图6 所示。
仿真表明,由于充分考虑了深海声场的声传播特性,基于本文多途信道阵列信号模型的DOA估计算法性能优于传统的基于平面波声场阵列信号模型的DOA估计算法;多途声线传播时间和相邻阵元多途声线传播时延差表征方法,以及多途信道下深海DOA估计阵列信号模型的建立,从根本上解决了深海DOA估计存在误差的问题。
4 结束语
为了提高深海DOA估计的性能,建立了深海多径声线传播时间和相邻阵元声线传播时延差的表征方法,以及基于相干多途信道的阵列信号模型,结合解相干DOA估计算法,提高了深海DOA估计精度。研究表明,深海DOA估计实质为多维参数优化问题;充分考虑海洋声场的声传播特性,可以从根本上解决深海DOA估计的误差问题。多个相干源的分辨和DOA估计,以及深海DOA估计多维参数优化的快速实现是下一步要研究的问题。