DOA Estimation Technology Based on Khatri-Rao Product for 3D Forward-looking Sonar

1.   引言

近年来，3维前视声呐作为一种重要的水下探测系统，在海底地形测绘、石油管道检测、沉船打捞、水下航行器避障、水下机器视觉等领域具有广泛的应用^[1-3]。目前，国外多个公司已研制出实用化的3维前视声呐设备。从功能原理上，国际上把3维前视声呐主要分为两类：(1)相控扫描式3维前视声呐，如TriTech公司的Eclipse^[4]，采用“T”型阵列结构，以水平向波束形成、垂直向扫描的方式完成3维探测，这种探测方式效率非常低；(2)多波束3维前视声呐，如CodaOctopus公司的Echoscope^4G[5]，采用2维方形接收面阵设计，采用2维波束形成技术完成3维探测，然而系统实现成本非常高；相比之下，Reson公司的SeaBat 7130^[6]采用2维矩形接收面阵设计，水平向采用波束形成技术，垂直向采用相干技术即可完成3维探测，相比Echoscope^4G，该方案具有较低的系统成本。

3维前视声呐作为一种2维方位主动探测声呐，回波信号具有强相干的特点，也即来自不同方位回波信号的皮尔森相关系数约为1。因而在实际工程应用中，在考虑改善目标方位分辨性能的同时，还需考虑2维方位估计失效和信号相干的问题。围绕提高方位分辨能力这一目标，陈建等人^[7]采用4阶累积量进行虚拟阵列孔径扩展，有效地改善了多目标分辨性能，然而4阶累积方法运算量非常大；杨永晶等人^[8]采用Kronecker积构造虚拟协方差矩阵，相比4阶累积方法计算复杂度有一定程度降低，然而在时间序列上的Kronecker积处理引入大量冗余信息，系统实时性不够友好；韩佳辉和谭伟杰等人^[9,10]直接对协方差矩阵进行向量化操作，基于Khatri-Rao积构造观测信号模型，相比Kronecker积处理，避免了时间维度引入的冗余，而仅仅对阵列孔径进行扩展，进一步降低了计算复杂度。此外，本文所述3维前视声呐还存在信号相干的问题，4阶累积量、Kronecker积等进行虚拟阵列扩展将出现协方差矩阵秩亏损的问题。为了恢复矩阵的秩，现阶段主要有降维处理和非降维处理两类方法进行解相干操作。其中降维处理方法以空间平滑类算法为代表，然而这类算法是以阵列孔径损失作为代价^[11]；非降维处理方法有效避免了阵列孔径损失^[12]。基于Khatri-Rao积方法在进行去冗余操作时隐含协方差矩阵平滑过程，有效避免了阵列孔径损失。值得注意的是，自协方差矩阵将抵消回波信号模型中垂直角时延分量(见本文3.1节)，应用于前视声呐的2维方位估计方法将失效，因而，上述4阶累积、Kronecker积、Khatri-Rao虚拟孔径扩展技术不能直接应用于文献[13,14]中所述的2维方位估计方法。

围绕提升目标分辨能力这一目标，本文重点解决现有阵列扩展技术存在2维目标方位估计失效、协方差矩阵秩亏损的问题，综合考虑2维稀疏面阵的结构特点，首先基于不同子阵构造互协方差矩阵，然后进行Khatri-Rao积扩展，提出单测量向量Khatri-Rao方法(Single Measured Vector Khatri-Rao, SVM-KR)和多测量向量Khatri-Rao方法(Multiple Measured Vector Khatri-Rao, MVM-KR)，并进行了理论推导、计算机仿真和工程实现验证。通过仿真实验，验证了所提出方法在目标分辨能力方面优于文献[13,14]中所述的单测量向量波束形成(Single Measured Vector Beamforming, SVM-B)方法，计算复杂度低于4阶累积量和Kronecker积虚拟阵列扩展方法；通过水池目标探测试验，验证了声呐系统和相关算法的工程实用性。

2.   问题阐述

2.1   稀疏阵列配置

自研前视声呐接收阵列配置如图1所示。接收阵列由3个均匀线性子阵构成，子阵间阵元间距分别为$5d$和${\rm{3}}d$，其中每个子阵由$M$个阵元构成，为了避免相位模糊阵元间距为半波长$d{\rm{ = }}\lambda /{\rm{2}}$^[15,16]。因此，该阵列由$M \cdot I$个阵元组成，$I{\rm{ = }}3,M = 48$。根据文献[1,13,17]，该2维稀疏面阵具有下面的优势：(1)适当改善目标方位估计分辨能力；(2)增强方位估计方法对于噪声的鲁棒性；(3)避免相位模糊问题。

图  1  2维稀疏接收阵列

下载: 全尺寸图片 幻灯片

2.2   信号模型

如图2所示，在快拍号$t$，考虑$K$个远场窄带相干回波信号同时入射2维接收阵列，第$i$行第$m$列阵元接收回波信号模型可表述为^[16]

图  2  3维笛卡尔坐标系下的目标方位表述

下载: 全尺寸图片 幻灯片

其中，${\left\langle {x\left( t \right)} \right\rangle _{i,m}}$为第$i$行第$m$列阵元接收回波信号；$s(t){\rm{ = }}{{\rm{e}} ^{{\rm{j}}2\pi ft/{f_{\rm{s}}}}}$为主动探测信号；$f$和${f_{\rm{s}}}$分别为工作频率和采样频率；$n(t)$为加性复高斯白噪声且与探测信号$s(t)$统计独立；${\theta _k}$和${\eta _k}$分别为水平角和垂直角；${{c}}$为声速；$\lambda = {{c}}/f$为波长；阵元间距为半波长$d{\rm{ = }}\lambda /{\rm{2}}$；$L$为快拍数。本文中，第1行第1列阵元设置为相位参考点。令${\omega _k} = 2\pi d\cos {\eta _k}\cos {\theta _k}/\lambda$, ${\phi _k} =$$2\pi d\sin {\eta _k}/\lambda$, ${\left\langle {{\tau _k}} \right\rangle _{i,m}} \approx m{\omega _k}{\rm{ + }}{\gamma _i}{\phi _k}$，由最小角定理，${\omega _k} = 2\pi d\cos {\vartheta _k}/\lambda$, ${\tau _k}$中第1项和第2项分别为空间角${\vartheta _k}$和垂直角${\eta _k}$的1阶多项式。式(2)中

第$i$个子阵接收回波信号矩阵形式为

其中，${{{X}}_i}{\rm{ = }}{\left[ {{{{x}}_{i,1}}\;{{{x}}_{i,2}}\; ··· \;{{{x}}_{i,M}}} \right]^{\rm{T}} } \in {\mathbb{C}^{M \times L}}$, ${{{x}}_{i,m}} =$${\left[{{{\left\langle {x\left( 1 \right)} \right\rangle }_{i,m}}\;{{\left\langle {x\left( 2 \right)} \right\rangle }_{i,m}}\; ··· \;{{\left\langle {x\left( L \right)} \right\rangle }_{i,m}}} \right]^{\rm{T}} }\in {\mathbb{C}^{L \times 1}}$, ${{S}} \in {\mathbb{C}^{K \times L}}$, ${{{N}}_i} \in {\mathbb{C}^{M \times L}}$与${{{X}}_i}$结构类似，不再展开赘述；${{{A}}_i}{\rm{ = }}$$\left[ {{{{a}}_{i,1}}\;{{{a}}_{i,2}}\; ··· \;{{{a}}_{i,K}}} \right] \in {\mathbb{C}^{M \times K}}$，令${u_k} = {{\rm{e}} ^{{\rm{j}}{\omega _k}}}$, ${v_k} = {{\rm{e}} ^{{\rm{j}}{\phi _k}}}$, 导向矢量可写作${{{a}}_{i,k}}{\rm{ = }}{\left[ {{u_k}v_k^{{\gamma _i}}\;u_k^2v_k^{{\gamma _i}}\; \cdots \;u_k^Mv_k^{{\gamma _i}}} \right]^{\rm{T}} } \in {\mathbb{C}^{M \times 1}}$, ${{\tilde{ A}}_i}{\rm{ = }}\left[ {{{{\tilde{ a}}}_{i,1}}\;{{{\tilde{ a}}}_{i,2}}\; ··· \;{{{\tilde{ a}}}_{i,K}}} \right] \in {\mathbb{C}^{M \times K}}$, ${{\tilde{ a}}_{i,k}}{\rm{ = }}{\left[ {{u_k}\;u_k^2\; ··· \;u_k^M} \right]^{\rm{T}} }$$\in {\mathbb{C}^{M \times 1}}$，也即${{{A}}_i}{\rm{ = }}v_k^{{\gamma _i}}{{\tilde{ A}}_i}$, ${{{a}}_{i,k}}{\rm{ = }}v_k^{{\gamma _i}}{{\tilde{ a}}_{i,k}}\;$。

3.   提出的算法

3.1   1维空间角估计方法

根据线阵几何性质，空间域的导向波束具有锥形形状^[18]，也即沿同一锥面入射的信号具有相同的空间角，如图3所示为空间角${\vartheta _1}$和${\vartheta _2}$的导向波束示意图。需要说明的是，根据第2节式(2)，如果构造常规协方差矩阵，也即${{{R}}_i} = {{{X}}_i}{{X}}_i^{\rm{H}}/L$，其中${{{A}}_i}{{S}}{{{S}}^{\rm{H}} }{{A}}_i^{\rm{H}}{\rm{ = }}v_k^{{\gamma _i}}{{\tilde{ A}}_i}{{S}}{{{S}}^{\rm{H}} }{\tilde{ A}}_i^{\rm{H}}v_k^{ - {\gamma _i}} = {{\tilde{ A}}_i}{{S}}{{{S}}^{\rm{H}} }{\tilde{ A}}_i^{\rm{H}}$，时延${\left\langle {{\tau _k}} \right\rangle _{i,m}}$中的${v_k}$分量将被抵消不能开展后续垂直角估计任务。本文综合考虑了上述2维面阵的特点，定义下面特殊的协方差矩阵

图  3  导向波束的3维波束图案

下载: 全尺寸图片 幻灯片

对上述协方差矩阵进行向量化操作

其中，$\otimes$表示Kronecker积。定义投影矩阵

其中，${{{e}}_{M + m}}$为单位向量，也即向量中除第$M + m$个元素为1其余皆为0。进而有

这里，${{B}} = \left[ {{{{b}}_1}\;{{{b}}_2}\; ··· \;{{{b}}_K}} \right] \in {\mathbb{C}^{2M - 1 \times K}}$为导向矩阵，${{{b}}_k} = {\left[ {u_k^{M - 1}\;u_k^{M - 2}\; ··· \;1\; \cdots \;u_k^{2 - M}\;u_k^{1 - M}} \right]^{\rm{T}} }$$\in {\mathbb{C}^{2M - 1 \times 1}}$, ${\mu _{21}}{\rm{ = }}v_k^{{\gamma _2} - {\gamma _1}}$, ${\mu _{32}}{\rm{ = }}v_k^{{\gamma _3} - {\gamma _2}}$, $\odot$表示Khatri-Rao积。显然地，相比SMV-B方法中的导向矩阵${{{A}}_i}$，阵列孔径大大增加，有望增强算法目标分辨能力。由于${{P}}$为满秩可逆矩阵，有

由于矩阵${{{P}}^{{\rm{ - }}1}}$左乘${\rm{vec}} \left( {{{{R}}_{21}}} \right)$近似等效为协方差矩阵平滑，${{{P}}^{{\rm{ - }}1}}\left( {{{A}}_1^* \otimes {{{A}}_2}} \right){\rm{vec}} \left( {{{{R}}_{\rm{s}}}} \right) \triangleq {{{P}}^{{\rm{ - }}1}}\left( {{{A}}_1^* \otimes {{{A}}_2}} \right)$${\rm{vec}} \left( {{{{\tilde{ R}}}_{\rm{s}}}} \right)$，这里${{\tilde{ R}}_{\rm{s}}}$为对角矩阵，也即${\rm{diag}} \left( {{{{\tilde{ R}}}_{\rm{s}}}} \right) =$${{q}} \in {\mathbb{C}^{K \times 1}}$, ${\rm{diag}} \left( \cdot \right)$为提取矩阵的主对角线元素构成向量的函数。根据Kronecker积和Khatri-Rao积的关系，当矩阵${{\tilde{ R}}_{\rm{s}}}$为对角矩阵时，有

进而，联系式(9)有

显然，根据式(12)可用1维方位估计算法完成空间方位估计${\hat \vartheta _k}$。下面依据Vernier法完成垂直角估计。

3.2   基于Vernier法^[17,19]的垂直角估计方法

由第2节可知，子阵间阵元间距大于半波长，进而目标回波相位差可以在$\left[ { - \pi ,\pi } \right]$范围内展开

其中，

其中，${\rm{angle}} \left( \cdot \right)$为计算复数相位角的函数。进而由式(13)得到目标回波方位的候选子集为

其中，$\left\lfloor \cdot \right\rfloor$为向下取整函数。根据互质特性，由候选子集${{{g}}_p}$和${{{g}}_q}$可得到目标回波方位估计

然而，式(16)成立建立在相位差$\Delta {\varphi _{21}}$和$\Delta {\varphi _{32}}$无偏的前提下，定义垂直方位估计的正则化表达式

3.3   基于最小角定理的水平角估计方法

对于水平角，由空间角估计${\hat \vartheta _k}$和垂直角估计${\hat \eta _k}$，进而依据最小角定理可得到水平角估计

至此，完成目标3维方位估计(包括水平角${\hat \theta _k}$和垂直角${\hat \eta _k}$)的相关算法，总结相关算法步骤如下：(1)由式(12)进行空间角估计${\hat \vartheta _k}$；(2)由式(14)计算子阵间相位差$\Delta {\varphi _{21}}$和$\Delta {\varphi _{32}}$；(3)由式(15)展开得到垂直角候选子集${{{g}}_p}$和${{{g}}_q}$；(4)由式(17)解算得到垂直角估计${\hat \eta _k}$；(5)结合步骤1和步骤4结果，由式(18)得到水平角估计${\hat \theta _k}$。

3.4   虚拟扩展阵列特性分析

(1)虚拟阵元接收信号相关性分析。如本文所述3维前视声呐，接收回波信号具有相干的特点，协方差矩阵的秩为1，应用本文所提方法进行虚拟阵列扩展，扩展后的虚拟阵元接收信号仍具有相干的特点。经过Khatri-Rao积向量化投影变换，等效为矩阵平滑，得到一个满秩协方差矩阵，可恢复为非相干信号模型^[9,10]。

(2)阵列多目标分辨能力分析。阵列波束形成的行为类似于空域滤波器，其阵列波束图案中的主瓣对应于滤波器的通带，表征了阵列对一定入射角度间隔的多目标分辨能力，阵列束宽越窄意味着其分辨能力越强。由式(12)可知，本文涉及的KR子空间方法导向矩阵${{B}} \in {\mathbb{C}^{2M - 1 \times K}}$理论孔径为$(2M - 1) \cdot d$，而SVM-B方法导向矩阵${{{A}}_n} \in {\mathbb{C}^{M \times K}}$等效阵元数为$M \cdot d$，KR子空间方法阵列孔径约为SVM-B方法的2倍。进而，SVM-B方法和KR子空间方法波束宽度可分别由式(19)、式(20)表述

其中，$\beta$为–3 dB加权窗系数，例如当窗函数为矩形窗时，$\beta \approx 0.443$。${\rm{arc}} \sin ( \cdot )$函数在$\left[ {0,\pi /2} \right]$区间内单调递增，我们有${\rm{BW}}_{\rm{KR},-3\;{\rm{dB}}}<{\rm{BW}}_{\rm{B},-3\;{\rm{dB}}}$，也即KR子空间方法阵列孔径得到扩展，其波束宽度减小，多目标分辨能力得到增强。下面两节开展仿真和水池试验进一步验证：(a) 多目标分辨能力；(b) 2维目标方位估计能力。

4.   仿真实验及分析

4.1   仿真实验

考虑两个窄带相干回波信号，水平角分别为90°和92°，垂直角分别为–7°和–5°，根据式(18)空间角分别为90°和91.99°。信噪比–5 dB，采用单测量向量4阶累计方法(Single Measured Vector Fourth-Order Cumulant, SMV-FOC)^[7]、单测量向量Kronecker积方法(Single Measured Vector Kronecker Product, SMV-KP)^[8]、SMV-B^[13,14]、本文第3节所述SMV-KR和MMV-KR进行空间角估计，SMV-FOC, SMV-KP, SMV-B和SMV-KR方法中$L = 1$, MMV-KR方法中$L = 3$，1000次独立的蒙特卡洛仿真实验水平角、垂直角估计误差棒图如图4、图5所示，其中15次输出波束图案如图6所示。

图  4  目标1(水平角90°垂直角–7°)仿真估计结果误差棒图

下载: 全尺寸图片 幻灯片

图  5  目标2(水平角92°垂直角–5°)仿真估计结果误差棒图

下载: 全尺寸图片 幻灯片

图  6  仿真实验输出波束图案

下载: 全尺寸图片 幻灯片

本文中取波束图案上两目标间幅值下降3 dB(0.707)作为分辨准则。由图6可见，SMV-B方法不能分辨出两个目标，而SMV-FOC, SMV-KP, SVM-KR和MMV-KR方法均能分辨两个目标，验证了本文所提方法在多目标分辨方面的优势，此外，上述4种方法估计均值均较为接近真值，验证了2维目标方位估计能力。相比SMV-B方法，SMV-FOC, SMV-KP, SVM-KR和MMV-KR方法具有更低的估计误差，这是由于估计误差随着阵列孔径的扩展而降低。进一步地，相比SVM-KR方法，MMV-KR方法具有更低的估计误差，这是由于估计误差随快拍数的增加而降低。

4.2   计算复杂度分析

上述5种算法的计算复杂度差异主要集中在1维空间角估计这一操作，本文重点对这一过程进行分析。计算复杂度由flops(floating-point operations per second)的次数近似表示，一个flops定义为一次浮点加法或乘法操作。SMV-B, SMV-FOC, SMV-KP, SMV-KR和MMV-KR方法的计算复杂度约为${\cal O}(3M \cdot B)$, ${\cal O}(4{M^5} + 8{M^4} - 2{M^3} + 13{M^2} - 8M +$$2 + 4MB \cdot (2M - 1))$, ${\cal O}(6{M^3} + 8{M^2} - 8M + 2 +$$4MB \cdot (2M-1))$, ${\cal O}(2(2M-1)$ $\cdot B+4{M^3})$和${\cal O}(4{M^3}+$$4{M^2} + 2B \cdot (2M - 1))$，因而计算复杂度方面所提出的算法优于SMV-FOC和SMV-KP。

5.   水池试验及分析

为进一步检验上述算法的有效性，本文依托哈尔滨工程大学研制的3维前视声呐于哈尔滨工程大学水声技术重点实验室信道水池开展实验研究。3维前视声呐关键参数如下。工作频率：150 kHz；发射波束宽度：60°×20°；接收波束宽度：2°×20°；阵列构成参数：$M{\rm{ = }}48$, $I{\rm{ = }}3$；探测信号：CW脉冲，脉宽0.2 ms。信道水池宽6 m，水深5 m，距离声呐基阵5 m处放置两个塑料球目标，声呐基阵置于水面以下2.5 m，目标置于水面以下2.2 m。水池实验场景示意图如图7所示。

图  7  水池实验场景示意图

下载: 全尺寸图片 幻灯片

如图8为上述SMV-B, SMV-FOC, SMV-KP, SVM-KR和MMV-KR方法应用波束形成技术输出的声呐图像，SMV-B方法不能很好地分辨两个目标，而SMV-FOC, SMV-KP, SMV-KR和MMV-KR能较好地分辨两个目标，为了更清楚地说明这一问题，从上述5种处理方法中抽取某一快拍下的波束图案如图9所示，SMV-B方法由于波束宽度较宽，不能有效分辨出两个目标，而SMV-FOC, SMV-KP, SMV-KR和MMV-KR方法由于扩展了阵列孔径，波束宽度较窄，能够很好地分辨两个目标。最后，应用上述3种方法得到目标估计结果如图10所示，SMV-FOC, SMV-KP, SVM-KR和MMV-KR方法能够对水池中目标水平角、垂直角进行估计，联系回波到达时间，可完成目标的3维定位。

图  8  波束形成输出声图像

下载: 全尺寸图片 幻灯片

图  9  某一快拍下的波束图案输出

下载: 全尺寸图片 幻灯片

图  10  2维方位估计结果

下载: 全尺寸图片 幻灯片

6.   结束语

本文依托2维SVM-B方位估计方法，基于Khatri-Rao积进行虚拟阵列扩展，提出1维空间角估计方法、基于Vernier法的垂直角估计方法和基于最小角定理的水平角方位估计方法。相比SVM-B方法原有的阵列流型，在保证2维方位估计方法有效的前提下，扩展了阵列孔径并改善了目标分辨能力。仿真实验验证了所提出的SVM-KR和MMV-KR方法在目标分辨性能优于SVM-B方法，且方位估计精度更高，计算复杂度低于SVM-FOC和SMV-KP，水池实验进一步验证了本文所提方法能够对多目标进行有效分辨，具备较强的工程实用性。