基于非对称幅度插值的二维来波方位估计

陈筱月, 马晓川,*, 李 璇

(1. 中国科学院声学研究所水下航行器信息技术重点实验室, 北京 100190;2. 中国科学院大学电子电气与通信工程学院, 北京 100049)

0 引 言

来波方位(direction of arrival, DOA)估计作为阵列信号处理的研究热点之一,在目标跟踪,声源定位以及无线通信等多个领域当中具有广泛的应用[1-4]。常规波束形成(conventional beamforming, CBF)[5]是最为经典的DOA估计算法,其通过阵列加权得到指向性波束,调整加权向量使波束指向不同方位从而得到方位谱,方位谱峰值对应的方位即真实来波方位。基于无失真约束与波束输出噪声最小准则,Capon[6]提出了最小方差无失真方法(minimum variance distortionless response, MVDR),利用接收信号协方差矩阵自适应的调节波束加权向量,提高了估计结果的准确性。子空间类算法[7-8]利用噪声子空间与信号子空间正交的特性,实现了超分辨。近年来,一些基于稀疏优化、解卷积以及神经网络的DOA算法被提出。稀疏优化类算法构造过完备基矩阵将DOA估计视作优化问题,并通过对解向量进行稀疏先验约束从而获得稀疏解[9-11];解卷积后处理算法基于逆滤波的原理,消除点扩散函数造成的方位谱扩展,从而得到理想方位谱[12-13];神经网络类算法将各个角度分为不同的类别,然后利用深度神经网络等方法进行分类,从而获得估计结果[14-15]。

上述大多方法获得高精度DOA估计结果的前提是进行密集的方位搜索。然而,密集的方位搜索会导致计算量的升高,特别是对于二维DOA来说[16],计算量会按平方急剧膨胀。而对于实际的声纳雷达系统,由于实时性要求及硬件系统的限制,难以通过密集的方位扫描来获得高精度的DOA结果。因此,实际工程当中一般会采用幅度插值[17]的方法得到来波方位,该方法基于CBF,但与传统密集谱峰搜索不同的是,该方法仅在少数方位上计算波束输出结果,然后利用不同的函数对输出结果进行方位谱拟合,插值得到来波方位。其中,最常用的拟合函数有线性函数、二次函数、高斯函数[18]。研究者针对幅度比较测向法的误差成因以及估计精度的提高做了大量相关工作。顾敏剑[19]使用高斯函数对波束图进行近似,理论分析了系统噪声,通道幅度特性不一致以及波束轴角指向对估计结果的影响。在高斯噪声下,文献[20]利用泰勒级数展开得到了振幅比较单脉冲雷达的DOA估计均方误差的解析表达式,避免了使用蒙特卡罗计算均方误差带来的高计算量。文献[21]中自适应地使用相邻或交替天线进行DOA测量,有效的降低了由于噪声与来波方位偏离波束指向所造成的误差。文献[22]中利用神经网络模型对一阶泰勒展开引起的非线性误差进行补偿,得到最终的DOA结果,提高了估计精度。文献[17]中考虑波束偏离0°时,左右波束的不对称,利用波束图-3 dB点作为求解条件,通过非对称函数进行方位谱拟合,得到了更好的估计结果。文献[23]中,提出了一种基于雷达模式高斯斜率修正的优化算法,取得了比传统算法更好的精度。而文献[24]直接在真实方向图基础上进行公式推导,进而对目标角度进行求解。上述工作大多均在一维进行,针对二维情况,文献[25]通过蜂窝网结构分布的二维多波束,使用二元二次方程进行曲面拟合,得到了较准确的测向结果。

本文首先通过对二维平面阵方位谱的推导,将二维幅度插值DOA估计解耦为两个相互独立的一维幅度插值DOA估计,从而避免了文献[25]中使用多元高次方程进行曲面拟合。对于解耦后的一维幅度插值DOA估计,考虑波束指向偏离中心时左右波束不对称,若使用对称函数进行插值会带来误差,所以使用了不同的非对称函数进行幅度插值,并且在角度求解的过程中,区别于文献[17]中使用-3 dB波束宽度作为求解条件,本文提出以左右波束斜率比作为求解条件,提高了算法的性能表现。针对二维幅度插值中出现的角度失配问题,提出了以迭代插值的方式对估计结果进行修正,进一步地提高了角度估计精度。最后,通过仿真实验分析了不同插值方法的误差曲线以及角度失配,信噪比对插值精度的影响,并通过对湖试数据的处理对算法进行了验证,本文算法能够应用于对实时性有要求的系统中进行高精度单目标方位估计。

1 信号模型

本文所使用的三维坐标系统[26]如图1所示,O为三维坐标系的原点,在此坐标系统中定义空间角为Ω=(θ,φ),其中θ为Ω与平面yOz的夹角,范围为[-90°,90°],而φ为Ω与平面xOz的夹角,范围为[-90°,90°]。在上述定义下,θ与φ实际上无法同时取到各自的边界值,两者之间存在着如下的限制关系:

图1 三维坐标系示意图Fig.1 Diagram of the three-dimensional coordinate system

|θ|+|φ|≤90°

(1)

即超出式限制的空间角Ω是不存在的。

二维平面阵位于xOy平面上,阵元个数为M,阵元坐标为pm=[pmx,pmy,pmz]T,则阵元坐标矩阵为P=[p1,p2,…,pM],在下文中,为了仿真的计算方便,使用的为四元矩形阵,且阵元间距d等于半波长,但是实际上本文方法对于任意二维平面阵均适用。现一窄带平面波s(t)由远场入射至二维平面阵上,其频率为f,波长为λ,声速为c,入射角度为Ω0=(θ0,φ0),其入射方向的单位方向矢量v0为

(2)

则阵列接收数据x(t)可以写为

x(t)=a(θ0,φ0)s(t)+n(t)

(3)

式中:n(t)为加性高斯白噪声;a(θ0,φ0)为(θ0,φ0)方向上的导向矢量,其表达式为

(4)

式(3)为连续形式,实际中会经过采样得到其离散形式如下:

x(n)=a(θ0,φ0)s(n)+n(n),n=1,2,…,N

(5)

式中:N为快拍数;x(n)∈CM×N;a(θ0,φ0)∈CM×1;s(n)∈C1×N;n(n)∈CM×N。

2 二维幅度插值分解

二维平面阵通过对采样得到的阵元接收信号进行加权求和,即可实现波束形成,对期望方向的信号进行输出。常规波束形成加权向量表达式如下:

(6)

波束输出信号为

y(θ,φ,n)=ω(θ,φ)Hx(n)

(7)

进而可以计算出其功率为

(8)

为了应对上述情况,实际系统中多采用幅度插值的方法进行角度估计。其原理如图2所示,图中的曲面是按照0.1°的间隔对整个空间进行扫描得到的方位幅度谱,可近似认为是一个连续的曲面(这里需要注意的是,尽管根据式(1),图中某些角度并不存在,但是为了之后插值与理解的方便,仍然计算了这些角度的波束输出值),红色星号所对应的位置(20°,10°)为通过波束密集扫描估计的来波方位。幅度插值法仅在少数的空间角度计算波束输出,如图2中的品红色圆点对应的角度,该操作相当于对真实的方位幅度谱进行空间采样。这里的空间采样间隔则与阵列在水平竖直方向上的波束宽度有关。对于矩形阵来讲,其在水平竖直两个方向上的分辨率相同,所以在这两个方向上的采样间隔也相同,如图3(a)所示。但对于任意阵,如2×4的五元L形阵,如图4(a)所示,其在水平竖直两个方向上的分辨率不同,因此在这两个方向上的采样间隔也不一样,在分辨率更高水平的方向上则需要更密的采样间隔,以免波束主瓣落在两个采样点之间,一般要求采样间隔小于主瓣宽度的一半。之后利用曲面函数对最大采样点及其周围多个采样点进行曲面拟合,插值得到的曲面最大值所对应的角度则为估计的来波角度。该方法避免了对二维空间进行密集扫描,有效的降低了计算量。

图2 二维幅度插值原理示意图Fig.2 Two-dimensional amplitude interpolation principle diagram

图3 二维方位幅度谱及其一维切片(矩形阵)Fig.3 Two-dimensional azimuth amplitude spectrum and its one-dimensional section (rectangular array)

图4 二维方位幅度谱及其一维切片(2×4 L形阵)Fig.4 Two-dimensional azimuth amplitude wpectrum and its one-dimensional section (2×4 L-shaped array)

然而,进行二维曲面拟合会面临复杂的多元高次方程组[25],求解较为困难,且面对任意阵型时,并不一定能够找到合适的拟合函数,从而会导致估计结果误差增加。因此,本文提出将二维幅度插值解耦为两个不同方向(θ方向,φ方向)的一维幅度插值,下面对该方式的可行性进行验证。对于式(8),在忽略噪声项后得到的幅度表达式如下:

(9)

进一步将式(6)代入式(9)可以得到

(10)

然而,θ0与φ0是所需要估计的量,实际上其具体值事先未知,仅知道在方位幅度谱上的等间隔采样值(即图3(a)与图4(a)中的黑色下三角点与品红色圆点处),而其中黑色下三角点的交叉处则是这些采样点中的最大幅度所对应的角度,记作(θm,φm)。因此,将二维曲面插值分解为(θ,φm)与(θm,φ)两个切面(黑色下三角点指示)上的一维曲线插值。图3(b)、图3(c)与图4(b)、图4(c)分别给出了四元矩形阵与2×4的五元L形阵的方位幅度谱及其在(θ0,φ0)的两个方向上(红色上三角)与(θm,φm)的两个方向上(黑色下三角)的一维切片。可以看出,对于矩形阵而言,其方位幅度谱在(θ0,φ0)的两个方向上与(θm,φm)的两个方向上的一维切片虽然存在幅度上的差异,但是最大值均所对应的角度相同。并且对于矩形阵来说,在任意角度的两个方向上的一维切片均满足此条性质。然而,对于其他任意二维平面阵,该条性质并不一定满足,如图4所示,L形阵方位幅度谱在(θ0,φ0)的两个方向上与(θm,φm)的两个方向上的一维切片极值点并不相同,若利用(θm,φm)的两个方向上采样点进行插值,则会产生误差,在本文中这种误差被称作角度失配误差。为了使算法能适用于任意二维平面阵,后续章节提出了使用多次迭代插值的方法来降低该项误差。

3 一维幅度插值

图5 φ=φm的一维切面Fig.5 One-dimensional section of φ=φm

3.1 对称幅度插值

3.1.1 对称线性插值

(11)

(12)

(13)

3.1.2 对称二次插值

二次插值则是用二次函数对方位幅度谱进行拟合,二次函数的表达式如下:

(14)

(15)

3.1.3 对称高斯插值

高斯插值则是用高斯函数对方位幅度谱进行拟合,其表达式如下:

(16)

(17)

3.2 非对称幅度插值

上述方法均认为对真实二维方位幅度谱进行切片得到的一维方位幅度谱是关于其极大值左右对称的。然而,事实上,这一点仅当一维方位幅度谱的极大值位于0°的时候满足,当其极大值偏离0°的时候,左右两侧将不再对称。如图6所示,为真实来波位于不同方位时的一维方位幅度谱,可以看出,随着来波角度偏离0°越严重,方位幅度谱的左右不对称性就越来越严重。因此,如果用对称函数进行插值拟合,会给估计结果带来较大的误差,为了解决这个问题,下面提出使用非对称函数进行插值拟合,以提高估计精度。

图6 不同来波方位下的方位幅度谱Fig.6 Azimuth amplitude spectrum in different wave directions

3.2.1 非对称线性插值

非对称线性插值认为一维方位幅度谱在其极大值左右两侧的斜率不同,因此线性拟合函数如下:

(18)

将α、β、γ与相应的波束幅度响应Uα、Uβ、Uγ(Uα≥Uγ)代入式(18)可以得到方程组:

(19)

(20)

(21)

(22)

(23)

3.2.2 非对称二次插值

非对称二次插值使用两个不同的二次函数对一维方位幅度谱极值点左右两侧进行拟合,两个二次函数仅在二次项系数上存在差别,其插值函数如下:

(24)

将α、β、γ与相应的波束幅度响应Uα、Uβ、Uγ(Uα≥Uγ)代入式(24)可以得到方程组:

(25)

(26)

A1=ηq(β)-1-M1ηq(β)+M1

A2=1-ηq(β)+M2ηq(β)-M2

B1=2β-2ηq(β)α+2M1ηq(β)α-2M1γ

B2=2ηq(β)α-2γ-2M2ηq(β)β+2M2γ

C1=ηq(β)α2-β2-M1ηq(β)α2+M1γ2

C2=γ2-ηq(β)α2+M2ηq(β)β2-M2γ2

(27)

(28)

3.2.3 非对称高斯插值

非对称高斯插值使用两个不同的高斯函数对一维方位幅度谱极值点左右两侧进行拟合,使用的非对称高斯插值函数如下:

(29)

将α、β、γ与相应的波束幅度响应Uα、Uβ、Uγ(Uα≥Uγ)代入式(29)可以得到方程组:

(30)

(31)

M3=[ln(Uα/Uβ)]/[ln(Uα/Uγ)]

M4=ln(Uγ/Uα)/ln(Uγ/Uβ)

A3=-ηg(β)+1+M3ηg(β)-M3

A4=-1+ηg(β)-M4ηg(β)+M4

B3=-2β+2ηg(β)α-2M3ηg(β)α+2M3γ

B4=-2ηg(β)α+2γ+2M4ηg(β)β-2M4γ

C3=-ηg(β)α2+β2+M3ηg(β)α2-M3γ2

C4=-γ2+ηg(β)α2-M4ηg(β)β2+M4γ2

(32)

(33)

4 仿真实验

(34)

下面将对二维幅度插值方法的估计误差进行讨论,并分析不同因素对算法估计结果精度的影响。

4.1 不同插值方式的影响

仿真使用的阵列为四元矩形阵并按照半波长布阵,来波信噪比为20 dB。在与z轴夹角小于50°的区域当中,按照10°的间隔进行幅度采样,使用上面提出的6种插值方法分别进行角度估计,并通过100次蒙特卡罗实验得到该区域内的二维误差曲面。图7给出了非对称高斯插值的估计误差曲面,为了方便对不同方法的估计精度进行观察与对比,之后不依次展示各个方法的估计误差曲面,而是展示各方法二维误差曲面在φ=0°处的切面,如图8所示。

图7 非对称高斯插值二维误差曲面Fig.7 Two-dimensional error surface of asymmetric Gaussian interpolation

图8 不同插值方法方位估计误差曲线(φ=0°切面)Fig.8 Azimuth estimation error curves of different interpolation methods (φ=0° section)

可以看出,随着入射角度逐渐偏离0°,估计误差呈现振荡上升。对称插值方法中,由于线性函数与方位谱的相似程度较差,其估计结果误差最大,而二次函数,高斯函数与方位谱具有较好的相似度,其拟合结果更加准确。而非对称插值相对对称插值,线性函数的整体估计误差并没有明显的改变,只是其误差曲线的极值点出现了移动,二次函数与高斯函数的误差出现了较为明显的下降,得到了更加精确的估计结果。进一步可以发现,非对称二次插值与非对称高斯插值的误差曲线的极小值与幅度采样角度一致,而误差曲线的极大值出现在两个相邻采样角度中间,即当入射信号方位在幅度采样角度附近的时候可以得到更小的估计误差。

4.2 角度失配的影响

由图8(b)可以看出,非对称二次插值与非对称高斯插值的误差曲线极大值点出现在两个相邻的幅度采样角度中点处,极小值点则刚好出现在幅度采样角度上。这是因为当真实角度刚好位于两个相邻幅度采样点中点的时候,真实角度离两侧的幅度采样点距离最远,所产生的角度失配最严重,误差较高;而入射角度刚好位于幅度采样点处时,几乎没有角度失配产生,所以误差很小。于是本文提出多次迭代插值,在每次迭代插值中,以上一次的插值结果为基础建立新的采样网格,这样,真实来波方位距离新的采样网格点更近,因此产生的失配误差就越小,下面分别对矩形阵与任意二维阵的迭代插值方式进行介绍。

经过多次迭代插值,即可以一步步地降低角度失配误差,使估计结果逼近真实值。下面通过仿真对提出的方法进行验证,仿真条件与第4.1节所使用的仿真条件相同,图9给出了非对称插值方法经过一次迭代插值后得到的误差曲线,*代表迭代插值结果。

图9 单次非对称插值误差估计曲线与一次迭代非对称插值误差估计曲线对比Fig.9 Comparison of the error estimation curve of single asymmetric interpolation with that of one-iteration asymmetric interpolation

由仿真可以看出,通过以第一次插值结果为基准,重新进行幅度采样并插值得到的估计误差有了明显的降低。实际上,对于矩形阵而言,只需要进行一次迭代插值就能够得到理想的估计结果,而对于任意二维阵,其角度失配更加严重,需要经过多次迭代插值才能够得到理想的估计结果。

4.3 不同求解条件的稳定性对比

从图10可以发现,常规波束形成的方位谱仅在高信噪比的情况下与波束图有较高的一致性,而随着信噪比的逐渐降低,两者之间的一致性也出现了明显的下降。在第4.2节进行非对称幅度插值的过程当中,为了对方程组进行求解,文献[17]引入了波束图的-3 dB点作为求解条件,本文则是定义了波束图的左右波束斜率比作为求解条件。而幅度插值实际上是针对方位谱进行插值,在低信噪比下,由于真实方位谱会相对波束图出现展宽,此时使用基于波束图的额外求解条件进行求解会导致最终的插值结果出现误差。尽管如此,由于无法事先知道方位谱,所以只能使用基于波束图的求解条件。这里想要说明的是,本文所使用的求解条件波束图斜率比对信噪比的变化并不敏感,其相对于-3 dB波束宽度更加稳定。

图10 不同信噪比下二维方位谱切片(φ=0°)Fig.10 Section of two-dimensional azimuth amplitude spectrum at different signal to noise ratios (φ=0°)

图11给出了不同来波方位下,方位谱的-3 dB宽度随信噪比的变化,图12给出了不同来波方位下,方位谱的斜率比随信混比的变化曲线。

图11 不同DOA下方位谱的-3 dB宽度随信噪比的变化曲线Fig.11 Variation curve of -3 dB width of azimuth spectrum with signal to noise ratio under different wave directions

图12 方位幅度谱的斜率比随信噪比的变化曲线Fig.12 Curve of slope ratio of the azimuth spectrum varies with the signal to noise ratio

从图11可以明显看出,方位谱的-3 dB宽度随着信噪比的降低逐渐变大,甚至当信噪比低到某个程度之后,方位谱的-3 dB宽度已经不存在了,并且对于小孔径的阵列来说,当来波方位偏离0°之后,其方位谱的一侧也不存在-3 dB宽度。因此,使用波束图的-3 dB点作为求解条件进行非对称插值,其估计结果受信噪比的影响较为严重,虽然文献[17]中通过粗估信噪比,并根据信噪比对插值过程中左右波束宽度进行修正,但这样仍然会存在较大的误差。而本文定义的斜率比随信噪比的降低变化较小,相对稳定,这代表虽然低信噪比下方位谱相对于波束图有所展宽,但是在这个过程中斜率比几乎不改变,因此使用波束图的斜率比进行非对称插值,能够在低信噪比下获得较高的插值精度。

4.4 信噪比的影响

接下来以10°作为采样间隔,在不同的信噪比下计算了不同插值方法在与图1中z轴之间空间夹角小于50°范围内的平均估计误差,经过100次蒙特卡罗计算得到的误差曲线如图13所示。

图13 不同插值方法估计误差随信噪比的变化曲线Fig.13 Estimation error curves of different interpolation methods varies with the signal to noise ratio

在图13(a)中,线性插值由于其与真实方位谱拟合较差,所以具有较高的估计误差;而二次函数插值与高斯插值方法的性能相近,在各个信噪比下均低于线性插值的误差,且在考虑真实方位谱的非对称性后其误差有明显的下降。同时,对比图13(b)可以发现,在经过一次迭代插值之后,对称线性插值误差存在略微下降,而非对称线性插值误差下降更加明显;对称高斯插值与对称二次插值误差反而有所升高,这是因为图8(a)中这两种方式的误差曲线的极小值点并不是在采样点附近,而是在两个采样点之间;而非对称二次插值与非对称高斯插值在低信噪比下误差相对单次插值改变不明显,但是在较高信噪比下(10 dB以上),出现了下降,取得了上述所有方法中最优秀的估计效果。

5 湖试数据处理

为了进一步对本文提出的非对称幅度插值方位估计算法的有效性与估计精度进行验证,下面将对实际的湖试采集数据进行处理。数据于杭州千岛湖采集,千岛湖是一个人造湖,水底存在一定起伏,水面较为平静。实验区域湖深约80 m,无水底山体干扰,实验时接收阵列与目标均布放在水深20 m附近。其中,接收阵列为均匀半波长布阵四元矩形阵,阵列朝向斜向下,由于一些其他原因,四元矩形阵的右下角阵元没有采集到有效数据,因此实际上使用的阵列为三元的L形阵列。目标为一人工声源,发射信号为10 kHz的连续波(continuous wave, CW)信号。

对接收信号按照0.01°的间隔在整个空间进行密集波束形成扫描,得到二维方位谱如图14所示,可以发现峰值出现在 (-17.85°,43.51°)的位置,该处即为密集扫描估计得到的来波方位。而在二维方位谱的左下角也存在峰值,需要说明的是,根据图2坐标系的定义,该峰值所对应角度|θ|+|φ|>90°,实际上这个空间角并不存在,因此将其排除。以二维方位谱最大值所对应的空间角作为基准,采样间隔为10°,在表1中给出了不同的幅度插值方法的估计结果进行对比。表1中,密集波束扫描的扫描间隔为0.01°,*代表进行了4次迭代插值。

表1 不同插值方法的估计结果Table 1 Estimated results of different interpolation methods (°)

图14 密集扫描方位幅度谱Fig.14 Azimuth amplitude spectrum with intensive scanning

从表1的插值结果可以看出,在只进行单次插值的结果中,非对称插值的估计结果虽然略微优于对称插值,但是整体的估计结果与密集波束扫描得到的结果相差较大。这是因为采集数据使用的阵列为三元L形阵列而非矩形阵,而对于任意阵来说,在幅度采样最大值处将其二维幅度插值分解为两个方向上独立的一维幅度插值会因为角度失配产生误差,这一点已经在第3节中进行了详细论述。为了降低由角度失配造成的误差,使用了在第4.2节中提出的多次插值的方法,经过4次迭代插值后,可以发现各种插值方法的估计精度有了较为明显的提升,其中非对称二次插值与非对称高斯插值的估计精度提升最为明显,其估计结果已经相当的接近密集波束扫描的结果。而对称插值法在经过多次插值后,其估计结果精度依然较差,事实上继续增加迭代插值次数,对称插值算法的估计精度也不会有较好的提升,这是由于这类插值方法没有充分的考虑方位谱的左右不对称性所导致。另外需要说明的是,虽然多次插值会给计算量带来一定的提升,但是每次迭代插值只需要多计算5次波束输出,而这所需的计算量是相当小的,远小于以0.01°间隔进行空间扫描的计算量。

6 结束语

本文提出了一种任意阵列的二维非对称幅度插值方位估计算法。算法将二维幅度插值解耦为两个独立的一维幅度插值,然后利用左右波束斜率比作为求解条件,通过非对称函数对少量离散波束输出进行幅度插值得到来波角度。进一步的,针对插值过程中的角度失配误差提出迭代插值的方法进行修正。通过仿真,本文提出的算法相对其他插值方法估计误差更小,湖试数据的处理结果说明,本文方法能够以较小的计算量取得与密集扫描常规波束形成类似的估计精度。本文算法具有如下优点:一是采用插值的方法进行方位估计,计算量很小;二是对阵型没有特殊要求,且能够在阵列孔径较小的情况下取得精确的估计结果。因此,本文算法可以应用于一些对实时性有要求的小型航行器上,实现对单目标的高精度实时方位估计。但由于本文算法是基于常规波束形成进行插值,目前只局限于对单目标的估计,如何将算法扩展至多目标,是下一步的研究目标。