1.   引言

基于压缩感知(Compressive Sensing, CS)理论的波达方向(Direction Of Arrival, DOA)估计方法，在近十年里相当流行，该方法最早由Malioutov等人^[1]提出。相对于离散化的方位角数目来说，感兴趣的信源数目是少量的。基于这个事实，DOA估计被描述成了一个具有稀疏约束的欠定线性系统，称为压缩波束形成(CompreSsive BeamForming, CSBF)^[2]。与传统的常规波束形成(Conventional BeamForming, CBF)相比，CSBF即使在单快拍下也具有稳健的高分辨精度。因此，CSBF被广泛应用于地震波反演、声学、雷达以及无线通信等领域。

CSBF的经典策略是施加${l_1}$范数的稀疏约束^[3]，因为在稀疏信号处理中${l_1}$范数具有特殊的意义，它是${l_0}$范数的凸松弛近似，并且凸问题的可靠求解相对容易^[4]。然而，自然的稀疏度量仍是${l_0}$范数，使用${l_1}$范数近似，会导致由于稀疏性不足，使得CSBF的性能下降^[5]。使用非凸的惩罚函数约束，具有比${l_1}$范数更强的稀疏性，可以提高算法的性能^[6,7]。但是，非凸的惩罚函数在求解过程中需要解决一系列的平滑问题，不能直接最小化惩罚函数，可能出现局部最小化问题^[8]。为了增强稀疏性且同时避免局部最小化问题，Selesnick^[9]提出一种极小极大凹惩罚(Minimax Concave Penalty, MCP)函数，在趋近${l_0}$范数稀疏性的同时保持了整体惩罚函数的凸性质。Yang等人^[10]将MCP函数应用到CSBF中，并证明了其优于${l_1}$范数约束方法的强大性能。

然而在低信噪比下，文献[10]的算法DOA估计性能不佳，多快拍联合估计可以提高低信噪比下的DOA估计结果。传统多快拍方法，如最小方差无失真方法(Minimum Variance Distortion-less Response, MVDR)^[11]和多重信号分类方法(MUltiple SIgnal Classification, MUSIC)^[12]等，对快拍的数量有着最低限制，且在低信噪比下性能不佳。基于${l_1}$范数的多快拍压缩波束形成(Multiple Snapshot Compressive BeamForming based on ${l_1}$ norm, ${l_1}$-MCSBF)^[13]虽然对快拍数没有要求，但在低信噪比下，由于稀疏性不足，会导致伪源的产生。

基于MCP约束的多快拍压缩波束形成(Multiple Snapshot Compressive BeamForming based on MCP, MCP-MCSBF)可以提供比${l_1}$-MCSBF的更好的性能。但是，求解MCP-MCSBF的困难在于需要证明MCP-MCSBF具备整体凸性质的前提条件，因为MCP函数本身是非凸的，只有满足整体凸性质，才能收敛到全局最小值，这也是其与其他非凸稀疏约束相比的优势所在。另外，类似文献[11]的常规多快拍压缩波束形成求解方法，是对每一个快拍单独处理并将结果求均值，其计算成本随着快拍数目的增大而增大。

因此对于上述问题，本文的贡献在于证明了MCP-MCSBF的整体凸性质条件，提出了一种低计算成本的MCP-MCSBF求解方法，实现了增强稀疏性的同时，在低信噪比情况下获得更稳定的DOA估计结果。本文提出的MCP-MCSBF首先对阵列接收的观测矩阵数据进行SVD分解，并将其投影到信号子空间上；然后再根据Boyd等人^[14]的凸优化理论，将投影后的观测矩阵所有列加权叠加成一个新的向量；最后，再将这个新的向量代入计算，得出精确的DOA估计结果。

2.   单快拍 DOA估计

2.1   信号模型

为了便于叙述，本文考虑远场理想均匀水平线阵模型。假设已知声速为1500 m/s，存在$K$个未知窄带源${s_k},k \in \{ 1,2,\cdots,K\}$从不同的方向角度入射。水平线阵阵元数为M。从源到每一个阵元的传播延迟由转向矢量描述，${\boldsymbol{a}}({\theta _k}) = [1,{{\rm{e}}^{{\rm{j}}(2\pi /\lambda )d1\sin {\theta _k}}}, {{\rm{e}}^{{\rm{j}}(2\pi /\lambda )d2\sin {\theta _k}}},\cdots,{{\rm{e}}^{{\rm{j}}(2\pi /\lambda )d(M - 1)\sin {\theta _k}}}]^{\rm{T}}$，$d$为两个相邻阵元间距，$\lambda$表示波长，${\theta _k}$代表第$k$个源的角度。

令阵列观测信号为${\boldsymbol{y}}$，并同样定义源信号${\boldsymbol{s}}$与加性噪声信号${\boldsymbol{n}}$，就可以将阵列观测信号建模为

其中，矩阵${\boldsymbol{A}}({\theta _k}) = [{\boldsymbol{a}}({\theta _1}),{\boldsymbol{a}}({\theta _2}),\cdots,{\boldsymbol{a}}({\theta _k})]$称为阵列流形矩阵。将感兴趣的角度空间${\theta _i} \in [ - {\text{9}}{{\text{0}}^{\text{o}}}{\text{,9}}{{\text{0}}^{\text{o}}}]$分解为$N$个可能的源方向，通常$N \gg K$。因此，源的数目是稀疏的，DOA估计问题就转化成稀疏重构问题

其中，${\boldsymbol{A}} \in {\mathbb{C}^{M \times N}}$，$\varepsilon \ge 0$是控制稀疏程度与拟合度的均衡。目前解决式(2)的方法有很多，包括快速迭代收缩阈值算法(Fast Iterative Shrinkage Thresholding Algorithm, FISTA)^[15]，利用MATLAB中的CVX工具箱的内点法求解^[16]等。但是，${l_1}$范数是${l_0}$范数一种非常松散的近似，所描述的稀疏结构并不令人满意，而稀疏性更强的非凸惩罚函数约束容易陷入局部最小化问题，所以有必要引入一个更加稀疏，且整体保持凸函数性质的惩罚函数。

2.2   非凸稀疏约束的压缩波束形成

MCP函数是由Huber函数^[17]的一个多元推广来定义的^[9]。其公式表达如式(3)所示

其中，${\boldsymbol{s}}$表示需要重构源信号向量，${\boldsymbol{u}}$表示引入的辅助变量，${\boldsymbol{\varSigma}}$为加权矩阵。将高维MCP函数${\phi _{\boldsymbol{\varSigma}} }({\boldsymbol{s}};{\boldsymbol{\varepsilon}} )$引入式(2)，就得到了基于MCP函数约束的压缩感知波束形成模型

MCP函数${\phi _{\boldsymbol{\varSigma}} }({\boldsymbol{s}};{\boldsymbol{u}})$是非凸的，因此可以比${l_1}$范数提供更强的稀疏性约束。根据参考文献[9]，当加权矩阵${\boldsymbol{\varSigma}}$满足${{\boldsymbol{\varSigma}} ^{\rm{T}}}{\boldsymbol{\varSigma}} \le {1/\varepsilon }{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}$时，式(4)的模型是凸的，通常加权矩阵${\boldsymbol{\varSigma}}$设置为${\boldsymbol{\varSigma}} = \sqrt {{\gamma /\varepsilon }} {\boldsymbol{A}}$。其中，$0 \le \gamma \le 1$是一个常数，实际中常使用$0.5 \le \gamma \le 0.8$。

将设置的固定加权矩阵${\boldsymbol{\varSigma}} = \sqrt {{\gamma/\varepsilon }} {\boldsymbol{A}}$代入式(4)，可得

式(5)具有凹凸结构，因此它在鞍点$({{\boldsymbol{s}}^ * },{{\boldsymbol{u}}^ * })$处具有最优的估计。本文提出基于块近端梯度下降算法框架^[18]的方法解决这个问题，实际上是两个FISTA变种，核心的思想是保持一个变量不变的同时，优化另一个变量。因此需要保持变量${\boldsymbol{s}}$不变，首先优化${\boldsymbol{u}}$

其中，$q$表示迭代次数。式(6)可以改写为

其中，${{\boldsymbol{L}}_{\boldsymbol{u}}} = \left\| {{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}} \right\|$是关于变量${\boldsymbol{u}}$利普希茨常数^[10]。为了加快收敛，参照FISTA算法^[15]将变量${\boldsymbol{u}}$的迭代重新设计为${\tilde {\boldsymbol{u}}^q} = {{\boldsymbol{u}}^q} + {{w}}_u^q({{\boldsymbol{u}}^q} - {{\boldsymbol{u}}^{q - 1}})$，并代入式(7)中可得

因此，式(8)的解表示为^[19]

接下来保持变量${\boldsymbol{u}} = {{\boldsymbol{u}}^{q + 1}}$不变，类似地优化变量${\boldsymbol{s}}$

对优化变量${\boldsymbol{s}}$执行类似变量${\boldsymbol{u}}$的操作可得

其中，${{\boldsymbol{L}}_{\boldsymbol{s}}} = (1 - \gamma )\left\| {{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}} \right\|$表示变量${\boldsymbol{s}}$的利普希茨常数，${\tilde {\boldsymbol{s}}^q} = {{\boldsymbol{s}}^q} + w_s^q({{\boldsymbol{s}}^q} - {{\boldsymbol{s}}^{q - 1}})$是变量${\boldsymbol{s}}$迭代的外推点。同上可得，式(11)的解表示为

根据参考文献[15]的加速策略，上文所提到的权值$w_u^q$和$w_s^q$可以被迭代的设置为

并且，权值${w^q}$的更新迭代遵循初始值${t_0} = 1$时

因为MCP函数的整体凸性质，所以向量${\boldsymbol{u}}$和${\boldsymbol{s}}$的初始值没有特别的限制，这里都设置为${ {\textit{0}}}$向量。算法的迭代次数可以任意设置，迭代的停止一个简单规则是利用两个相邻迭代向量的变化，${{{\left\| {{{\boldsymbol{s}}^{q + 1}} - {{\boldsymbol{s}}^q}} \right\|}_2}}/ {{{\left\| {{{\boldsymbol{s}}^q}} \right\|}_2}} \le {\rm{tol}}$，其中，${\rm{tol}}$为一个较小的常数。具体的算法步骤如表1所示。

表 1  基于MCP函数约束的压缩波束形成算法步骤

输入：观测数据${\boldsymbol{y}}$，参数$\varepsilon$，超参数$\gamma$，最大迭代次数$Q$，终止条 　　　　件${\rm{tol}} = {10^{ - 3} }$。参数向量${t_0} = 1$，${{\boldsymbol{u}}^0}$和${{\boldsymbol{s}}^0}$的初始值均设为 　　　　${ {\textit{0}}}$向量，且${{\boldsymbol{u}}^0} = {{\boldsymbol{u}}^1}$, ${{\boldsymbol{s}}^0} = {{\boldsymbol{s}}^1}$。
(1) 当$q \le Q$时：
(2) 　　根据式(14)和式(13)更新$w_u^q$，
(3) 　　获得变量${\tilde {\boldsymbol{u}}^q}$，并按照式(9)更新变量为${{\boldsymbol{u}}^{q + 1} }$，
(4) 　　根据式(14)和式(13)更新$w_s^q$，
(5) 　　同理获得变量${\tilde {\boldsymbol{s}}^q}$，并按照式(12)更新变量为${{\boldsymbol{s}}^{q + 1} }$。
(6) 　　如果${{\boldsymbol{s}}^{q + 1} }$和${{\boldsymbol{s}}^q}$的变化$\le {\rm{tol}}$：
(7) 　　　那么停止迭代，且最优估计${{\boldsymbol{s}}^ * } = {{\boldsymbol{s}}^{q + 1} }$，
(8) 　　　找到${{\boldsymbol{s}}^ * }$的峰值，并将其对应到DOA估计的角度${{\boldsymbol{\theta}} ^ * }$。
(9) 　　结束
(10) 结束
输出：DOA估计的角度${{\boldsymbol{\theta}} ^ * }$。

下载: 导出CSV 

| 显示表格

3.   多快拍DOA估计

使用多次观测数据的联合DOA估计，比单快拍具有更稳定的估计结果^[20]。本文近似地认为源信号的DOA在连续的多个快拍数据里都是相同的。将阵列接收到的多快拍测量数据建模成一个矩阵${\boldsymbol{Y}} \in {\mathbb{C}^{M \times L}}$

其中，$L$表示快拍数目，${\boldsymbol{N}} \in {\mathbb{C}^{M \times L}}$为噪声矩阵，矩阵${\boldsymbol{S}} \in {\mathbb{C}^{N \times L}},{\boldsymbol{S}} = [s(1),s(2),\cdots,s(L)]$表示出相同的行稀疏性。

对多快拍观测数据矩阵${\boldsymbol{Y}} \in {\mathbb{C}^{M \times L}}$进行SVD分解^[21]可得

其中，${\boldsymbol{U}}$表示$M \times M$的酉矩阵，${\boldsymbol{B}}$表示对角线上具有非负实数的对角线$M \times L$矩阵，${\boldsymbol{V}}$表示$L \times L$的酉矩阵。通过估计${\boldsymbol{B}}$中较大的奇异值个数$\hat K$，并且得到${\boldsymbol{V}}$中与奇异值对应的信号子空间${{\boldsymbol{V}}_{\hat K}} = [{\boldsymbol{V}}(1), {\boldsymbol{V}}(2),\cdots,{\boldsymbol{V}}(\hat K)]$。因为噪声与信号子空间正交，通过正交投影，从而将多快拍观测数据矩阵截断，得到降噪后的观测矩阵${{\boldsymbol{Y}}_{\hat K}} \in {\mathbb{C}^{M \times \hat K}}$

其中，${{\boldsymbol{D}}_{\hat K}} = {[{{\boldsymbol{I}}_{\hat K}}{\text{ }}{{ {\textit{0}}}}]^{\rm{T}}}$，${{\boldsymbol{I}}_{\hat K}}$是$\hat K \times \hat K$的单位阵。

此时基于MCP约束的多快拍压缩波束形成的模型可以表示为

其中，${\phi _{\boldsymbol{\varSigma}} }({{\boldsymbol{S}}_{\hat K}};{{\boldsymbol{U}}_{\hat K}}) = {\left\| {{{\boldsymbol{S}}_{\hat K}}} \right\|_{2,1}} - \mathop {\min }\limits_{{{\boldsymbol{U}}_{\hat K}}} \left\{ {{\left\| {{{\boldsymbol{U}}_{\hat K}}} \right\|}_{2,1}} + {1/ 2} \left\| {\boldsymbol{\varSigma}} ({{\boldsymbol{S}}_{\hat K}} - {{\boldsymbol{U}}_{\hat K}}) \right\|_{\rm{F}}^2 \right\}$。那么式(18)可以进一步写成

显然${\phi _{\boldsymbol{\varSigma}} }({{\boldsymbol{S}}_{\hat K}};{{\boldsymbol{U}}_{\hat K}})$是非凸的，所以式(19)也是非凸的。为了消除局部最小值，保证式(19)的整体凸性质，定理1给出了约束条件。

定理1　定义函数${\boldsymbol{F}}({{\boldsymbol{S}}_{\hat K}}) = {1/2} \left\| {{\boldsymbol{Y}}_{\hat K}} - {\boldsymbol{A}}{\boldsymbol{{S}}_{\hat K}} \right\|_{\rm{F}}^2 + \varepsilon {\phi _{\boldsymbol{\varSigma}} }({{\boldsymbol{S}}_{\hat K}};{{\boldsymbol{U}}_{\hat K}})$，当加权矩阵${\boldsymbol{\varSigma}}$满足${{\boldsymbol{\varSigma}} ^{\rm{T}}}{\boldsymbol{\varSigma}} \le {1/ \varepsilon }{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}$时，${\boldsymbol{F}}({{\boldsymbol{S}}_{\hat K}})$是凸函数。

证明　由式(19)可知，${\boldsymbol{F}}({{\boldsymbol{S}}_{\hat K}})$可以写成

根据$|| \cdot ||_{\rm{F}}^2$的计算公式，可以将式(20)中的矩阵F范数转化为

对于向量$|| \cdot ||_2^2$的计算公式，式(21)可以转化为

将式(22)展开并合并同类项可得

其中，$j$表示矩阵${{\boldsymbol{S}}_{\hat K}} \in {\mathbb{C}^{N \times \hat K}}$的第$j$个快拍，$\hat K$是快拍数目。最后一项$G({{\boldsymbol{S}}_{\hat K}};{{\boldsymbol{U}}_{\hat K}})$是$|| \cdot ||_{\rm{F}}^2$展开之后${{\boldsymbol{S}}_{\hat K}}$和${{\boldsymbol{U}}_{\hat K}}$的相关项，它是一组凸函数的最大值^[22]，因此也是凸的。所以，当条件${{\boldsymbol{\varSigma}} ^{\rm{T}}}{\boldsymbol{\varSigma}} \le {1/\varepsilon }{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}$满足时，函数${\boldsymbol{F}}({{\boldsymbol{S}}_{\hat K}})$是凸的。证毕

同2.2节，将加权矩阵${\boldsymbol{\varSigma}}$设置为${\boldsymbol{\varSigma}} = \sqrt {{\gamma/\varepsilon }} {\boldsymbol{A}}$，式(23)可以改写为

显然，可以参考2.2节的单快拍算法，将式(24)看成两个交替求解的多快拍稀疏约束问题求解。但问题在于传统求解多快拍稀疏约束的算法，如多快拍交替方向法(Multiple snapshots Alternate Direction Method, M-ADM)^[23]等，本质上只是对每一个快拍单独运算并求和平均。这样的做法存在两个问题：(1)本质上仍未脱离单快拍运算的范畴，只能一定程度上减少估计结果的误差，对于高噪声环境下导致的稀疏约束能力不足，从而产生的估计偏差，并没有得到真正的解决；(2)计算成本与$\hat K$成正比。

本文提出MCP-MCSBF算法将降维后的观测矩阵${{\boldsymbol{Y}}_{\hat K}}$通过加权求和，进一步表示成$M \times 1$维向量${{\boldsymbol{y}}_{\hat K}}$，然后将其代入单快拍求解算法中求解。且可以通过定理2证明，当权数${\hat w_1} + {\hat w_2} + \cdots + {\hat w_{\hat K}} = 1$，${{\boldsymbol{\varSigma}} ^{\rm{T}}}{\boldsymbol{\varSigma}} \le {1/ \varepsilon }{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}$时，目标函数仍具备整体凸性质。

定理2　令加权向量$\hat {\boldsymbol{w}} = {[{\hat w_1},{\hat w_2},\cdots,{\hat w_{\hat K}}]^{\rm{T}}}, {\hat w_i} > 0, i = 1,2,\cdots,\hat K$，定义

当且仅当${\hat w_1} + {\hat w_2} + \cdots + {\hat w_{\hat K}} = 1$，${{\boldsymbol{\varSigma}} ^{\rm{T}}}{\boldsymbol{\varSigma}} \le {1/\varepsilon } {{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}$时，${\hat {\boldsymbol{s}}_{\hat K}}$是整体凸的。

证明　首先由2.2节已知，当${{\boldsymbol{\varSigma}} ^{\rm{T}}}{\boldsymbol{\varSigma}} \le {1 \mathord{\left/ {\vphantom {1 \varepsilon }} \right. } \varepsilon }{{\boldsymbol{A}}^{\rm{H}}}{\boldsymbol{A}}$时，式(4)中的目标函数$\hat s$是整体凸的；其次，定义集合${\boldsymbol{\varOmega}} = \left\{ {{{\boldsymbol{Y}}_i},i = 1,2,\cdots,\hat K} \right\}$，同时集合${\boldsymbol{\varOmega}}$中的元素组合加权叠加可以表示为集合${\boldsymbol{\varTheta}}$

根据文献[14]中1.1节第1部分可知，当${\hat w_1} + {\hat w_2} + \cdots + {\hat w_{\hat K}} = 1$时，集合${\boldsymbol{\varTheta}}$是包含集合${\boldsymbol{\varOmega}}$的凸集，因此${\hat {\boldsymbol{s}}_{\hat K}}$是整体凸的。证毕

本文所提MCP-MCSBF算法如表2所示。本文算法的收敛性证明可见文献[10,24]，因此这里不再重复说明。

表 2  MCP-MCSBF算法步骤

输入：多快拍观测数据矩阵${\boldsymbol{Y}}$。
(1) 根据式(16)和式(17)计算降噪后的观测矩阵${{\boldsymbol{Y}}_{\hat K} }$，
(2) 构造加权向量$\hat {\boldsymbol{w}}$，满足${\hat w_1} + {\hat w_2} + \cdots + {\hat w_{\hat K}} = 1$，
(3) 根据式(25)计算${{\boldsymbol{y}}_{\hat K} }$，
(4) 　将${{\boldsymbol{y}}_{\hat K} }$代入式(26)中得到全新的目标函数，
(5) 　根据表1的算法得到DOA估计角度${{\boldsymbol{\theta}} ^ * }$。
输出：DOA估计的角度${{\boldsymbol{\theta}} ^ * }$。

下载: 导出CSV 

| 显示表格

4.   仿真结果

本节主要探讨基于MCP约束的多快拍压缩波束形成的DOA估计性能。考虑了一个传感器数量$M = 24$且半波长分离的均匀线性阵列，源信号为载波频率3000 Hz的窄带信号，噪声为服从高斯分布点的白噪声信号。传感矩阵$A$定义在一个间隔为${{\text{1}}^{\text{o}}}$，从–90°～90°的角度网格中，并使用均方根误差(Root Mean Squared Error, RMSE)量化DOA估计的性能。

4.1   本文算法的性能

首先，考虑单快拍与提出多快拍算法之间的性能对比。假设有3个能级相同的源信号，它们的入射角分别为$[ - {\text{4}}{{\text{0}}^{\text{o}}}{\text{,}} - {\text{1}}{{\text{5}}^{\text{o}}}{\text{,2}}{{\text{0}}^{\text{o}}}]$。信噪比SNR=–10 dB，快拍数分别设置为1和50，算法的最大迭代次数$Q$设置为3000。图1展示了MCP约束的单快拍与多快拍压缩波束形成算法的性能。从图1(a)可以看出，由于较强噪声的存在，MCP约束的单快拍压缩波束形成算法产生了多个能量很强的伪源，且DOA估计结果存在偏差，估计性能较差。而MCP约束的多快拍压缩波束形成算法，在准确定位了3个源信号入射角度的同时，也不存在错误的峰值，如图1(b)所示。

图 1  信噪比–10 dB下的DOA估计结果

下载: 全尺寸图片 幻灯片

其次，比较MCP约束的传统多快拍与提出多快拍压缩感知算法的性能。这里的传统多快拍压缩感知算法指的是，通过对每一个快拍单独处理，然后将结果求和平均的方法，见文献[11,23]。仿真设置与上文相同，仿真结果如图2所示。图2(a)展示了传统多快拍压缩感知算法的估计性能。由于该算法只是对每个快拍分别计算然后求和平均，没有根本上解决强噪声背景下，稀疏性不足的问题，因此可以看到图2(a)中尽管估计出了准确的方位角度，但是存在两个较强的伪源，干扰了目标判断。图2(b)为提出多快拍算法的估计结果，展示了稳健而精确的方位估计性能。

图 2  信噪比–10 dB下多快拍的DOA估计结果

下载: 全尺寸图片 幻灯片

4.2   与传统多快拍DOA估计算法的性能对比

为了进一步验证MCP-MCSBF算法的性能。本文将其与CBF, MUSIC以及${l_1}$-MCSBF对比。首先，考察SNR为–10 dB下不同的算法的DOA估计结果。令3个能级相同的独立源信号，入射角分别为[–4°, 0°, 20°]，快拍数为50，仿真结果如图3所示。在图3(a)和图3(b)中由于–4°和0°的两个源间隔较小，CBF和MUSIC均无法有效地区分这两个目标源；而图3(c)中${l_1}$-MCSBF虽然能够区分这两个目标源，但由于${l_1}$-MCSBF的稀疏性不足，所以出现了许多的伪峰。MCP-MCSBF具有比${l_1}$-MCSBF更强的稀疏性，因此在图3(d)中，不仅能区分间隔较近的两个目标源，且除了估计的目标源外，不存在别的伪峰。

图 3  –10 dB下不同算法的DOA估计结果对比

下载: 全尺寸图片 幻灯片

其次，分别考虑独立源和相干源情况。SNR从–10 dB变化到10 dB，其余设置与上文相同，每次均进行400次蒙特卡罗仿真模拟，得到的RMSE结果如图4所示。图4(a)表示独立源情况下RMSE随SNR变化的情况。可以看出，随着SNR的增大，这4种算法的DOA估计精度都逐渐提高。因为CBF波束宽度近似为${{{{102}^\circ}}/{M = }}{4.25^\circ}$，大于$- {{\text{4}}^\circ}$和${{\text{0}}^\circ}$的两个源的角度间隔，所以CBF的估计精度要明显低于另外3种算法。MUSIC估计精度在高信噪比下较好，但在低信噪比下估计精度较差。${l_1}$-MCSBF和MCP-MCSBF估计精度相似，但在低信噪比下MCP-MCSBF具备的估计精度更高。图4(b)表示相干源情况下RMSE随SNR变化的情况。由于MUSIC算法在估计DOA过程中，涉及接收信号协方差求逆，因此相干源会导致DOA估计精度大大降低。CBF，${l_1}$-MCSBF以及MCP-MCSBF则不受相干源影响，仍具备很高的DOA估计精度。

图 4  不同算法的性能比较

下载: 全尺寸图片 幻灯片

然后，比较信噪比–10 dB，快拍数从1～100变化下的DOA估计结果，如图4(c)所示。这里考虑到阵元个数为24，因此MUSIC算法的快拍起点为30快拍。图中所有算法的RMSE都随着快拍数增多而下降，但是可以清楚地看到本文提出的MCP-MCSBF在任何快拍下都要优于CBF，MUSIC以及l₁-MCSBF。

最后，从SNR和角度间隔两个方面考察不同算法的角度分辨率。对于不同算法的DOA估计结果，如果满足$\sum\nolimits_{k = 1}^2 {\left| {{{\hat \theta }_{ik}} - \theta _{ik}^{{\text{true}}}} \right|} < \left| {\theta _2^{{\text{true}}} - \theta _1^{{\text{true}}}} \right|$，则认为成功区分了两个独立源，将成功分辨的次数比上蒙特卡罗仿真次数，得到成功分辨概率，可以用来表示算法的性能。

考虑SNR从–10 dB变化到10 dB，两个能级相同的独立源信号，入射角分别为$[{0^\circ}{\text{,}}{{\text{5}}^\circ}]$，快拍数为50，每次均进行400次蒙特卡罗仿真模拟。不同算法的成功分辨概率结果随SNR的变化如图5(a)所示。接下来考虑角度间隔${2^\circ}$从变化到${12^\circ}$的两个独立源，SNR为–10 dB。不同算法的成功分辨概率结果随角度间隔的变化如图5(b)所示。显然，当角度间隔较大时，所有算法的成功分辨概率都很高；但当角度间隔较小时，本文提出的MCP-MCSBF算法成功分辨概率要远优于其他4种算法。

图 5  成功分辨率随SNR和角度间隔的变化

下载: 全尺寸图片 幻灯片

5.   实测数据结果

上文已经在仿真中验证了算法的性能，但是数值的模拟并不能代表真实的水下环境。本节展示算法在水下实测数据中的优良性能。采用湖试数据作为本文提出的MCP-MCSBF算法的实测数据。湖面声速为1500 m/s，由UUV平台搭载的舷侧均匀水平线阵接收数据，UUV航行深度在水下6 m，阵元数目为24个，阵元之间的距离为0.1875 m。此次试验中存在两个声源信号，一个是位于水下6 m处的静止声源信号，发射频率为1～4 kHz宽带噪声信号，另一个为水面运动的小艇发出的辐射声源，且小艇与静止声源之间存在交叉。将每秒的数据分为50快拍，使用全频带的信息进行DOA估计。

图6表示CBF, MUSIC以及${l_1}$-MCSBF与MCP-MCSBF算法的时间方位历程图。图6(a)表示CBF的方位历程图，可以看出当两个声源逐渐接近时，由于CBF较宽的主瓣，将两个声源合并成了一个声源，且存在自噪声和干扰噪声的影响，导致目标声源不够突出；图6(b)表示已知精确的信号源个数这一先验信息下的MUSIC的方位历程图，MUSIC算法相对于CBF来说主瓣更窄，对自噪声和干扰噪声有了很好的抑制，因此目标声源更加突出，但是对背景噪声的抑制力不强，最低能量为–10 dB；图6(c)表示${l_1}$-MCSBF的方位历程图，显然${l_1}$-MCSBF相对于CBF以及MUSIC来说，在抑制自噪声和干扰噪声的同时，大大减弱了背景噪声，此时最低能量为–100 dB，增强了目标声源。但是由于稀疏性不足，所以方位历程图中存在较多的噪点；图6(d)表示MCP-MCSBF的方位历程图，与${l_1}$-MCSBF相比，由于本文提出算法增强了稀疏性，因此图中存在的噪点很少。且与传统CBF和MUSIC算法相比，目标声源更加清晰，对背景噪声的抑制也更强，此时最低能量为–150 dB。

图 6  方位历程图

下载: 全尺寸图片 幻灯片

通过湖试数据可以看出，本文方法在噪声背景下保持着稳健的DOA估计性能，抑制了旁瓣的产生；同时由于较窄的主瓣宽度，也提高了目标声源的分辨率，在目标声源角度间隔较近时，展现了良好的性能。

6.   结论

本文提出一种基于MCP约束的多快拍压缩波束形成算法，本方法通过非凸函数MCP提高了传统压缩波束形成算法的稀疏性，且由于MCP的整体凸性质，不存在局部最小值问题，保证了DOA估计的可靠性；同时利用多个快拍的数据联合估计，提高了DOA估计的稳定性。仿真分析证明了本文算法与CBF, MUSIC以及${l_1}$-MCSBF算法相比，具有更优的精确性和更高的角度分辨率，湖试结果也再一次证明了本文算法的有效性。