引言

水声通信技术是在水下进行无线信息传输的重要手段。然而, 水下声传播环境复杂, 边界的多次反射或物体的散射会造成典型的多径传播, 导致接收端形成多个具有不同时延和幅度衰减的发射波形副本的叠加^[1]。此外, 由于水中声速仅为1500 m/s, 收发双方之间以及反射界面的相对运动, 会产生显著的多普勒效应。上述影响的共同作用, 导致了水声信道的时延–多普勒双扩展现象, 给水声信道估计、水声定位、水声通信等带来了挑战^[2]。为了更为精准的获取信道特征, 实现高可靠的水声通信和高精度的水声定位, 对时延和多普勒参数同时进行估计是十分必要的。

匹配滤波(MF)是进行时延–多普勒参数估计的经典方法, 它通过估计接收信号与已知发送信号之间的互模糊函数(CAF)实现^[3,4]。Yu等通过在宽带模糊函数基础上引入分数阶傅里叶变换方法, 提升了标准匹配滤波算法的性能^[5]。然而, 受限于时间带宽积, 基于互模糊函数算法的分辨率有待提高。基于稀疏表示的信道脉冲响应估计近年来得到了广泛关注, 水声信道天然的稀疏特性, 也加速了相关技术的发展。以正交匹配追踪(OMP)为代表的压缩感知算法, 广泛应用于水声信道参数估计^[6,7]。Qu等考虑到大字典矩阵下的计算复杂度, 提出了一种两级稀疏信道估计技术, 大大减少了传统正交匹配追踪算法在网格上搜索的候选对象数量, 降低了计算复杂度^[8]。伍飞云等通过引入非均匀范数约束, 解决了传统算法容易进入局部最优和二维搜索导致运算复杂度高等问题^[9]。Sun在传统正交匹配追踪的基础上, 引入了Gram-Schmidt正交过程, 避免了逆运算, 提高了时延–多普勒的估计精度和效率^[10]。Cang等基于稀疏表示理论和解卷积思想, 提高了时延估计的分辨率^[11]。Wei等在传统算法中引入动态调整机制, 实现了稀疏度自适应, 显著降低了计算量^[12]。

然而, 上述方法都是基于高斯噪声模型背景展开的, 在极地环境中, 由于冰面的破裂、碰撞、挤压和摩擦等运动, 会造成大量的高频次、大幅度冲激噪声^[13,14]。这类噪声的概率密度函数具有较高的拖尾厚度, 无法使用高斯分布对其进行有效的表征。基于广义中心定理提出的对称$\alpha$稳定(${\text{S}}\alpha {\text{S}}$)分布是用来表示脉冲噪声的一种简单且常见的模型^[15,16]。该分布下, 噪声不具有二阶及以上高阶特征矩, 因此之前提到的基于高斯噪声模型假设的传统估计方法的性能将严重下降^[17]。

针对脉冲噪声下的水声信道参数估计, 研究人员提出了多种优化算法以提高脉冲噪声下的估计性能。黄健等采用基于分数低阶协方差谱的广义互相关函数, 实现了鲁棒的时延估计^[18]。陈梦等分析了脉冲噪声导致估计算法失效的原因, 并设计了一种新型分段非线性幅值变换函数^[19]。Dong等通过在传统匹配滤波器中引入${L_p}$范数约束, 有效地抑制了脉冲噪声干扰, 并利用实测脉冲噪声对算法的鲁棒性进行了验证^[20,21], 但该算法求解复杂度较高。Tian等采用基于M估计的鲁棒正交匹配追踪(RobOMP)算法, 实现了脉冲噪声下的信道估计, 并提高了后续信号检测能力^[22]。然而, 上述方法仅考虑了信道的时延特性, 冰下平台的运动和极地区域常见的冰面漂移等现象, 会造成时延–多普勒耦合效应, 导致相关算法的时延估计性能不佳^[10]。因此, 为了准确获取水声信道参数, 需要进一步对时延–多普勒参数开展联合估计。

针对上述问题, 为了实现极地脉冲噪声环境中水声信道参数的准确稳健估计, 本文结合最大相关熵准则(MCC)和${L_1}$范数约束, 提出了一种改进的正交匹配追踪算法。首先, 利用MCC对于脉冲噪声的抑制特性, 实现原子基的有效选择; 其次, 利用${L_1}$范数重构目标函数, 实现脉冲干扰下的参数求解。同时, 通过引入基于分布式迭代优化策略的交替方向乘子法(ADMM)加快了求解速度。数值仿真和基于中国第九次北极科考冰下实测脉冲噪声数据处理结果验证了算法的有效性和鲁棒性。

1.   基本模型与理论

1.1   时延–多普勒双扩展信道模型

水声信道可视为线性时变信道, 信道响应$h\left( {\tau ,t} \right)$可以采用时变多径冲击模型描述为^[23]

其中, $\delta (\cdot)$表示冲击函数, $Q$是多径数目, ${a_q}\left( t \right)$和${\tau _q}\left( t \right)$分别为第$q$条路径上的幅度衰减和路径时延。幅度的衰减主要是由界面反射和介质的吸收引起的; 多途路径的不一致性导致了接收信号副本对应时延的差异。

同时, 由于水下声传播速度较慢, 信号收发双方之间的相对运动会导致显著的多普勒现象。针对这一问题, 假设发射信号$s\left( t \right)$的持续时间为$T$, 在一个同步周期$T$内, 幅度衰减${a_q}\left( t \right)$和时间延迟${\tau _q}\left( t \right)$可进行如下近似表示

式(2)中假设同步周期内路径上的幅值衰减是不变的。式(3)中${\tau _q}^\prime$为路径的初始时延因子, ${d_q}^\prime$为第$q$条路径上与发射源和接收器之间相对运动速度相关的系数。根据射线声学理论, ${d_q}^\prime$的大小取决于传播路径上收发双方之间的相对径向速度${v_q}$与声传播速度$c$的比值^[10], 即

其中, $\boldsymbol{v}$为运动矢量, $\| \cdot \|$为取模值, ${\theta _q}$为路径$q$所对应声线的掠射角。将式(2)和式(3)代入式(1)中, 可进一步将时变多径信道模型描述为

1.2   噪声声压统计模型

高斯噪声是最为常见的一种海洋背景噪声, 然而冰层的大范围形变、破裂、挤压、摩擦和融化等过程会产生大量的脉冲声, 导致单一的高斯噪声模型已无法有效的表征极地噪声背景。相关研究表明, 脉冲噪声可以使用${\text{S}}\alpha {\text{S}}$分布进行声压统计建模, 其特征函数可表示为

其中, $w$为表示声压的随机变量; $\alpha \in \left( {0,2} \right]$是表征概率密度函数拖尾厚度的特征指数, $\alpha$越大, 偏离中值的样本个数越少, 脉冲特性越弱; 当$\alpha = 2$时, ${\text{S}}\alpha {\text{S}}$分布退化为高斯分布; $\gamma \in \left( {0, + \infty } \right)$是表示分布范围的尺度参数, 其作用与方差类似, 用于表征噪声的起伏程度。

1.3   相关熵理论

相关熵可以用来度量两个随机变量$X$和$Y$之间的相似性, 其定义为^[24]

其中, ${\rm E}[\cdot]$表示期望, ${k}_{\sigma }(\cdot)$表示核函数, $\sigma$为核长。高斯核函数是使用最为广泛的一种核函数, 可表示为

实际上, 由于随机变量的联合概率密度函数难以获取, 一般可以通过有限的采样数据$\left\{ {\left( {{x_i},{y_i}} \right)} \right\}_{i = 1}^N$计算相关熵的估计值, 其中$N$为采样点数, 即

根据相关熵的有界性${M_\sigma }\left( {X,Y} \right) \leqslant {M_\sigma }\left( 0 \right)$和对称性${M_\sigma }\left( {X,Y} \right) = {M_\sigma }\left( {Y,X} \right)$, 可以得到相关熵诱导距离(Correntropy Induced Metric, CIM):

即CIM越小, 两个变量之间的相似程度越高。基于CIM距离测度, 最大相关熵准则被定义为

由式(8)可以看出, 相关熵本质上是一种受限的“局部测量”方法, 当误差大于由核长控制的某一阈值时, 相关熵趋近于0, 因此, 只有处于约束范围内的样本才能够有效地参与相似性度量的计算, 根据这一特性, 可以有效地抑制脉冲噪声。

2.   基于最大相关熵准则和${L_1}$范数约束的时延–多普勒估计方法

假设发射信号为$s\left( t \right)$, 则经过多径传播后的接收信号$y\left( t \right)$可表示为

其中, $\otimes$表示卷积, $h\left( {\tau ,t} \right)$为式(5)所示的时变多径冲击信道响应, $w\left( t \right)$为式(6)所表征的脉冲噪声。接收信号可以进一步写为

为了便于表述, 定义等效时延${\tau _q} = {{{\tau _q}^\prime } \mathord{\left/ {\vphantom {{{\tau _q}^\prime } {\left( {1 + {d_q}^\prime } \right)}}} \right. } {\left( {1 + {d_q}^\prime } \right)}}$, 等效多普勒因子${d_q} = 1 + {d_q}^\prime$, 接收信号可简化为

即接收信号是发射信号分别经历$Q$个不同时延、多普勒缩放以及幅度衰减之后的总和。

对式(14)进行离散化处理, 并定义如下的参数集合:

其中, ${d_q} = {d_0} + (q - 1){\varDelta _d}$, ${\tau _q} = {\tau _0} + (q - 1){\varDelta _\tau }$, ${d_0}$和${\tau _0}$表示多普勒因子和时延估计的最小可能值, ${\varDelta _d}$和${\varDelta _\tau }$分别表示对应的搜索精度, ${N_d}$和${N_\tau }$表示可能的多普勒因子和时延个数。因此, 接收信号${\boldsymbol{y}}$的离散形式可写为

其中, ${\boldsymbol{y}} = {[{y_1},{y_2},\cdots,{y_L}]^{\text{T}}}$, ${\boldsymbol{w}} = {[{w_1},{w_2},\cdots,{w_L}]^{\text{T}}}$分别表示接收信号和噪声的离散采样, $L$为接收信号长度。发射信号${\boldsymbol{s }}$经不同多普勒因子作用并循环移位构成的${N_d}$个Toeplitz矩阵共同组成字典矩阵${\boldsymbol{G}}$, 可表示为

其中, ${N_s}$为发送信号长度; ${s_m}\left( {{d_n}} \right)$表示对发射信号$s$的第$m$个采样点施加多普勒因子${d_n}$的作用。${\boldsymbol{h}}$是${N_d}{N_\tau }$维的待估计向量:

可以通过求解式(17)对时延–多普勒参数进行估计。时延多普勒参数共有${N_d}{N_\tau }$个组合作为候选, 对应于待估计信道${\boldsymbol{h}}$的长度。由于待求解未知量数目远大于观测数量$L$, 所以方程是严重病态的, 具有无穷多组解。根据式(5)可知, 水声信道具有稀疏特性, 即${\boldsymbol{h}}$中大量元素为0, 因此可以利用稀疏表示理论求解该欠定方程, 实现对稀疏向量${\boldsymbol{h }}$的重构, 从而进一步找到与稀疏解相对应的多普勒因子和时延参数。

压缩感知(CS)是实现稀疏估计的一种有效方法, 其中以OMP为代表的贪婪追踪算法在高斯噪声下的双扩展信道估计中已经得到了广泛研究, 其将式(17)表示为${L_2}$范数约束下的优化问题:

其中, ${\Vert \cdot \Vert }_{0}$和${\Vert \cdot \Vert }_{2}$分别为${L_0}$范数和$L_2$范数, $Q$为预设的稀疏度, 即信道多径数目。其具有表1所示的典型实现流程。

表  1  典型OMP算法实现流程

输入	接收信号${\boldsymbol{ y}}$, 字典矩阵${\boldsymbol{G}}$, 稀疏度$Q$
初始化	迭代次数索引$i = 1$, 残差${{\boldsymbol{e}}_i} = {\boldsymbol{y}}$, 原子索引集合${\boldsymbol{J }}= \emptyset$, ${\boldsymbol{h}}$的估计值$\widehat {\boldsymbol{h}} = {\mathbf{0}}$
迭代	1) 计算残差与各原子的内积, 并选择内积最大的原子, $j=\arg\max_j\left|\boldsymbol{c}_j^{\mathrm{T}}\boldsymbol{e}_i\right|$(其中$\boldsymbol{c}_j$表示字典矩阵${\boldsymbol{G}}$中第$j$个原子, 即${\boldsymbol{G}}$的第$j$列), 并将索引$j$放入集合${\boldsymbol{ J}}$中, 即${\boldsymbol{J}} = {\boldsymbol{J}} \cup \left\{ j \right\}$; 2) 计算信道的估计值$\widehat {\boldsymbol{h}} = {\left( {{{\boldsymbol{G}}_J}^{\text{T}}{{\boldsymbol{G}}_J}} \right)^{ - 1}}{{\boldsymbol{G}}_J}^{\text{T}}{\boldsymbol{y}}$, 这里的${{\boldsymbol{G}}_J}$表示由1)中挑选出来的原子$\boldsymbol{c}_j\left(j\in\boldsymbol{J}\right)$组成的基原子矩阵; 3) 更新$i = i + 1$, 如果$i < Q$, 返回步骤1), 否则输出$\widehat {\boldsymbol{h}}$。
输出	信道估计值$\widehat {\boldsymbol{h}}$

下载: 导出CSV 

| 显示表格

典型OMP算法进行脉冲噪声下的信号处理存在两个问题。一是内积是在${L_2}$空间中进行定义的, 其并不具备抗离群值的能力, 因此不适合用于表征脉冲噪声下向量间的相关性(相似度)。二是由于除高斯分布($\alpha = 2$时)外的${\text{S}}\alpha {\text{S}}$分布不存在二阶矩, 因此采用${L_2}$范数约束构建的信道参数求解目标函数在脉冲噪声下失效。错误的信道估计会传递到下一次的原子基选择与参数估计中, 导致错误的连续发生。基于上述问题, 本文提出了一种基于MCC和${L_1}$范数约束的改进OMP方法, 利用相关熵作为相似性度量准则, 增强脉冲噪声下原子选择的鲁棒性; 采用${L_1}$范数重新构建目标函数, 增强信道估计的准确性。表2中描述了算法的具体流程。

表  2  基于MCC和${L_1}$范数约束的改进OMP算法实现流程

输入	接收信号${\boldsymbol{ y}}$, 字典矩阵${\boldsymbol{ G}}$, 稀疏度$Q$
初始化	迭代次数索引$i = 1$, 残差${{\boldsymbol{e}}_i} = {\boldsymbol{y}}$, 原子索引集合${\boldsymbol{J }}= \emptyset$, ${\boldsymbol{h}}$的估计值$\widehat {\boldsymbol{h}} = {\mathbf{0}}$;
迭代	1) 计算残差与各原子的相关熵, 并选择相关熵最大的原子, 即$j=arg\max_j\left|M_{\sigma}\left(\boldsymbol{c}_j^T,\boldsymbol{e}_i\right)\right|$, 并将索引$j$放入集合${\boldsymbol{J}}$中, 即${\boldsymbol{J}} = {\boldsymbol{J }}\cup \left\{ j \right\}$; 2) 计算信道的估计值$\widehat {\boldsymbol{h }}= \arg \mathop {\min }\limits_{\boldsymbol{h}} {\left\| {{\boldsymbol{y}} - {{\boldsymbol{G}}_{\boldsymbol{J}}}{\boldsymbol{h}}} \right\|_1}$(这里采用ADMM方法求解, 下文中有详细说明); 3) 更新$i = i + 1$, 如果$i < Q$, 返回步骤1), 否则输出$\widehat {\boldsymbol{h}}$。
输出	信道估计值$\widehat{\boldsymbol{ h}}$

下载: 导出CSV 

| 显示表格

相关熵对信号的“局部相似性”进行度量, 能够有效去除局部异常值的影响。同时, 相关熵的核长对噪声并不十分敏感, 无需${\text{S}}\alpha {\text{S}}$分布的先验知识, 相比于常规的基于低阶统计量参数选取的方法, 具有更好的使用前景^[25]。

在表2的步骤3)中, 为了抑制异常脉冲值的影响, 采用${L_1}$范数对式(22)中的优化目标进行重构, 并得到如下的优化目标:

该问题可以采用标准的优化工具如CVX工具箱进行求解。但该方法效率较低, 求解时间过长, 不利于信号的实时处理, 下面通过引入ADMM实现全局最优解的高效获取。

ADMM解决优化问题的思路是将难以求解的全局问题拆分为多个容易求解的局部子问题, 并进行迭代优化和交替求解。对式(23)引入辅助变量${\textit{z}}$, 可以到如下的目标函数

其中, $\varepsilon$为正则化参数, 用于控制待求解向量${\boldsymbol{h}}$的稀疏度在整个优化函数中的权值。文献[26]中的试验表明, 当$\varepsilon$的取值在0.4~0.8之间时, 上述问题可以取得较好的结果。式(24)对应的增广拉格朗日形式为

其中, ${\boldsymbol{u}}$为对偶变量, $\rho > 0$是与增广相关的惩罚系数。ADMM通过优化下面三个子代价函数, 对变量${\boldsymbol{ z }}$、${\boldsymbol{h }}$和${\boldsymbol{u}}$进行交替求解

式(26)和式(27)的求解问题, 本质上是${L_2} - {L_1}$问题, 这里通过对式(27)所示目标函数的二次项进行线性化进行近似求解, 在${{\boldsymbol{q}}^k}$处有

其中, $d\left( {{{\boldsymbol{h}}^k}} \right) = {{\boldsymbol{G}}_{\boldsymbol{J}}}^{\rm{T}}\left( {{{\boldsymbol{G}}_{\boldsymbol{J}}}{{\boldsymbol{h}}^k} - {{\boldsymbol{q}}^k}} \right)$, ${{\boldsymbol{q}}^k} = {\boldsymbol{y}} + {{\boldsymbol{z}}^{k + 1}} + {{{{\boldsymbol{u}}^k}} /\rho }$, $\kappa > 0$为近端参数。通过线性化近似, $h$的求解问题可用软阈值运算符表示为

软阈值算子定义为^[27]

其通过保留元素大于门限的部分与元素的原始符号, 并将其余元素用0取代, 以保证所求解的稀疏性。同样${\boldsymbol{z}}$的解可表示为

对于惩罚系数$\rho$, 其最优值的设定目前没有较为简单的方法^[28], 通常采用经验值, 文献^[29]给出了一种适用于ADMM算法的修正方案, 它通过一个渐进系数, 在每次迭代中对$\rho$值进行更新, 从而保证算法的求解效果。近端参数$\kappa$用于控制每次迭代计算的步长, 当选择较大值时, 会加速算法收敛, 但会损失精度, 反之精度提升, 但收敛速度变慢^[30]。此外, 还需要设置相邻两次估计的差值门限$\xi$对ADMM算法的迭代过程进行控制, 同时, 设置最大迭代次数$K$防止算法的过度迭代。ADMM求解算法的具体实施步骤如表3所示。

表  3  ADMM算法求解流程

输入	原子索引集合${{\boldsymbol{G}}_{\boldsymbol{J}}}$, 残差矩阵${\boldsymbol{y }}$, 设置超参数$\varepsilon$, $\rho$, $\kappa$, $K$和$\xi$
初始化	迭代次数索引$k = 1$, 待估计中间变量${{\boldsymbol{z}}^k} = 0$, ${{\boldsymbol{h}}^k} = 0$, ${{\boldsymbol{u}}^k} = 0$
迭代	1) 根据式(26)计算${{\boldsymbol{z}}^{k + 1}}$; 2) 根据式(27)计算${{\boldsymbol{h}}^{k + 1}}$; 3) 根据式(28)计算${{\boldsymbol{u}}^{k + 1}}$; 4) 更新$k = k + 1$, 如果${\left\| {{{\boldsymbol{h}}^k} - {{\boldsymbol{h}}^{k - 1}}} \right\|^2} < \xi$或$k > K$, 则输出结果$\widehat {\boldsymbol{h}} = {{\boldsymbol{h}}^{k + 1}}$, 否则返回步骤1)。
输出	信道估计值$\widehat {\boldsymbol{h}}$

下载: 导出CSV 

| 显示表格

3.   数值仿真及实验数据分析

3.1   极地实测脉冲噪声统计特性分析

为了验证对于脉冲噪声采用${\text{S}}\alpha {\text{S}}$分布进行建模的合理性, 首先对采集的真实北极冰下脉冲噪声进行拟合分析。冰下环境噪声采集实验于中国第九次北极科学考察期间(2018年8月18—19日)进行, 在位于(84°09′47.10′′N, 167°15′44.384′W)的长期冰站开展。采集设备由8个间隔10 m的自容式水听器组成, 覆盖了冰下10 m到80 m的范围, 接收设备采样率为50 kHz, 水听器工作频带范围为10 Hz~20 kHz, 灵敏度为$- 170\;{\text{dB}} \pm 2\;{\text{dB}} \; {\text{re}} \; {\text{1}}\,\text{V} / \mu {\text{Pa}}$。图1(a)是第5个采集设备采集到的北极噪声时域波形, 可以明显地看到许多大能量值异常脉冲, 其对应的频谱如图1(b)所示。可知北极冰下脉冲噪声具有丰富的频率成分, 涵盖了几乎所有常用频段, 因此无法通过简单的滤波手段消除脉冲噪声的干扰。

图  1  实测北极冰下脉冲噪声及分析结果 (a)北极冰下脉冲噪声的时域波形; (b)脉冲噪声的频谱图

下载: 全尺寸图片 幻灯片

图2是图1(a)中北极冰下噪声的概率密度函数以及采用高斯分布和${\text{S}}\alpha {\text{S}}$分布的拟合结果。这里采用Fama和McCuloch等提出的基于顺序统计量的分位数方法^[31,32]对参数$\alpha$和$\gamma$进行估计, 高斯分布的估计是基于均值和方差的。图2(a)中蓝色实线表示真实的脉冲噪声分布; 红色实线表示${\text{S}}\alpha {\text{S}}$分布的拟合结果, 其中$\alpha {\text{ = 1}}{\text{.44}}$和$\gamma {\text{ = 1087}}{\text{.97}}$; 粉色实线表示高斯分布拟合结果。从图中可以看出, 冰下脉冲噪声具有严重的拖尾现象, ${\text{S}}\alpha {\text{S}}$分布具有更好的拟合效果, 证明了采用${\text{S}}\alpha {\text{S}}$分布对冰下脉冲噪声进行建模的合理性, 为了更直观地展示拟合效果, 图2(b)中将脉冲声压幅值大小限制在$[-{\text{10}}^{\text{4}}, {\text{10}}^{\text{4}}]$。

图  2  脉冲噪声概率密度函数拟合结果 (a)原始拟合结果图; (b)拟合结果局部放大图(将幅值限制在$[-{\text{10}}^{\text{4}}, {\text{10}}^{\text{4}}]$)

下载: 全尺寸图片 幻灯片

3.2   数值仿真分析

在数值仿真中, 发射信号$s\left( t \right)$采用广泛使用的线性调频(LFM)信号, 起始频率为${f_0}$, 脉冲持续时间为$T$的LFM信号可表示为

其中, $k$为线性调频速率。$s\left( t \right)$在时刻$t$的瞬时频率为$f\left( t \right) = {f_0} + kt$, 发射信号的带宽为$B = f\left( T \right) - f\left( 0 \right) = kT$。仿真所用LFM信号起始频率为1 kHz, 带宽为2 kHz, 采样率为8 kHz, 信号持续时间为100 ms。

图3(a)中是仿真使用的声速剖面, 其来自于2018年采集的真实北极数据, 采集地点的经纬度为76.6233N,163.8606W, 其具有北极海域典型的双声道特性。仿真中, 声源在180 m深度以6 m/s的恒定速度做水平运动, 接收器固定在200 m深度, 接收器与声源之间的水平距离为3 km。接收信号通过将发射信号与模型化的信道冲击响应(CIR)进行卷积生成, 信道冲激响应由Bellhop声线传播模型生成。图3(b)是仿真得到的声线轨迹, 图3(c)中是对应的信道冲激响应。

图  3  基于真实极地声速剖面的Bellhop仿真结果 (a)某北极冰站CTD声速剖面图; (b)模拟的声线轨迹; (c)仿真的信道冲击响应

下载: 全尺寸图片 幻灯片

仿真过程中采用式(6)表述的${\text{S}}\alpha {\text{S}}$分布对脉冲噪声进行仿真生成。高斯噪声背景下, 通常采用信噪比(SNR)表征信号和噪声的相对强度, 其定义为

其中, $P\mathrm{_S}$表示信号功率, $P_{\mathrm{N}}$表示噪声功率。但由于除$\alpha = 2$之外的${\text{S}}\alpha {\text{S}}$分布不存在二阶矩, 即$P_{\mathrm{N}} \to \infty$。此时, 传统信噪比在脉冲噪声下失去意义, 因此这里采用广义信噪比(Generalized SNR, GSNR)表征信号和噪声的相对强度, 其定义为

其中, $\gamma$为前文提到的${\text{S}}\alpha {\text{S}}$分布的尺度参数。

为了衡量不同方法的性能, 采用均方误差作为评价标准, 时延和相对速度的均方误差分别为

其中, $N_{\mathrm{mc}}$是独立重复实验的次数, ${\widehat \tau _n}$和${\widehat v_n}$为第$n$次实验中的估计值, ${\tau _n}$和${v_n}$分别为对应的真值。

3.2.1   典型方法性能对比

从信道的仿真结果中可以提取不同声线对应的掠射角${\theta _q}$, 并根据式(4)得到不同路径上的多普勒因子。这里只考虑图1中的前五条路径, 信道的具体参数如表4所示, 之后的实验均采用此信道。

表  4  仿真信道参数

路径编号	1	2	3	4	5
幅值	1	0.54	0.99	0.71	0.44
时延 (s)	2.079	2.111	2.137	2.226	2.247
相对速度 (m/s)	6	5.4	5.8	4.6	5.4
多普勒因子	0.0040	0.0036	0.0039	0.0031	0.0036

下载: 导出CSV 

| 显示表格

根据式(4), 以相对速度的估计值表征多普勒估计结果。仿真过程中, 相对速度的搜索空间设置为0~8 m/s, 搜索精度设置为0.2 m/s。时延的搜索范围设置为0~200 ms, 搜索精度${\varDelta _\tau } = {1 \mathord{\left/ {\vphantom {1 {{f_s} = }}} \right. } {{f_s} = }}0.125{\text{ ms}}$。因此, 相对速度(多普勒)维度${N_d} = 41$, 时延维度${N_\tau } = 1600$, 接收信号长度$L = 2400$, 字典矩阵$G$的规模为$\left[ {L,{N_d}{N_\tau }} \right] = \left[ {2400,65600} \right]$。仿真过程中添加特征指数$\alpha = 1.5$, GSNR=0 dB的脉冲噪声, 并假设稀疏度$Q = 5$已知。图4是添加脉冲噪声前后的时域波形, 从图4(b)可以看到, 接收信号的尾段受到了较为强烈的脉冲离群值干扰。

图  4  信号时域波形 (a)无噪声接收信号时域波形; (b)加入脉冲噪声后接收信号时域波形

下载: 全尺寸图片 幻灯片

分别选取了GSOMP^[10], 基于${L_p}$范数的匹配滤波器(${L_p}$-MF) ^[18]和RobOMP^[20]三种方法进行对比。对于${L_p}$-MF算法, 选取$p = 1$与本文的${L_1}$范数约束相对应; 对于RobOMP算法, 采用具有代表性的“Bisquare”损失函数。对于所提方法的ADMM解算过程, 根据文献[29]给出的指导方案, 正则化参数$\varepsilon$设置为0.5, 惩罚系数$\rho$设置为5, 近端参数$\kappa$设置为100, 停止门限$\xi$和最大迭代次数$K$分别设置为${10^{ - 4}}$和30; 此外, 高斯核函数核长的选取也会影响算法的性能, 为了对参数进行合理选择, 首先对不同核长下的算法性能进行了仿真对比, 仿真过程中设置GSNR=0 dB。仿真结果如图5所示。

图  5  不同核长下本文算法性能仿真结果(GSNR = 0 dB) (a)时延均方误差变化曲线; (b)等效速度均方误差变化曲线

下载: 全尺寸图片 幻灯片

从仿真结果可以看出, 核长越小, 对于强脉冲噪声的抑制性能越好, 这是由于小核长对应于更小的“局部”范围, 但同时, 当脉冲噪声减弱时, 由于其有效数据利用范围过小, 导致性能下降; 反之, 当核长增大时, 系统在弱脉冲噪声下的效果提升, 但是对于强脉冲噪声的抑制能力变弱。通常, 冰下脉冲噪声的特征指数在1 ~ 2之间, 为了保证算法对于不同脉冲噪声的鲁棒性, 根据仿真结果将核长$\sigma$设置为2。

下面, 对不同算法的性能进行对比, 图6中是采用不同方法得到的时延–相对速度估计结果, 右侧二维图中红色圆圈表示添加的真值。${L_{\text{p}}}$-MF算法能够分辨前四条路径, 但无法准确的估计相对速度, 并且由于接收信号尾部脉冲噪声的存在, 基于相关计算的匹配滤波的输出结果被大能量异常值污染, 导致第五条路径被淹没, 无法被准确识别, 此外该方法只能确定全局最大值点, 对其他路径的局部最大值不易识别, 无法进行后续处理。由于假定稀疏度已知, 三种基于压缩感知的方法都展示了五条路径, 但由于脉冲噪声的影响, GSOMP算法完全失效, 不能提供准确的信息, RobOMP算法虽然有一定的脉冲噪声抑制能力, 但还是在第四条路径上出现了错误, 而本文提出的MCC-${L_1}$-OMP算法则展示了完全正确的估计结果。

图  6  不同方法得到的估计结果 (a)(b) ${L_{\text{p}}}$-MF方法($p = 1$); (c)(d) GSOMP方法; (e)(f) RobOMP方法; (g)(h) MCC-${L_1}$-OMP方法

下载: 全尺寸图片 幻灯片

3.2.2   统计性能分析

为了对比各种算法的统计性能, 分别进行了不同脉冲噪声强度以及不同广义信噪比下的仿真实验。由于${L_{\text{p}}}$-MF方法分辨率较低, 难以有效提取估计参数, 因此不参与之后的性能比较。

(1) 不同脉冲噪声强度下的仿真结果

对比不同算法在不同脉冲噪声强度下的表现, 仿真过程中接收信号的GSNR设为0 dB, 脉冲噪声特征指数$\alpha$从1增长到1.8, 步长为0.1, 每个脉冲噪声强度下进行1000次独立重复实验, 即${N_{{\text{mc}}}} = 1000$。图7给出了不同方法的仿真实验结果。可以看出, 在广义信噪比保持不变的情况下, 随着脉冲噪声强度的减弱, 脉冲噪声极大异常值出现的频次越来越低, 对于采用内积进行原子选择和采用LS准则进行求解带来的影响也逐渐降低, 因此GSOMP的估计性能不断提高。但由于该方法是针对高斯噪声环境设计的, 在脉冲噪声下均方误差仍停留在较高水平。RobOMP和MCC-${L_1}$-OMP都能够有效抑制脉冲噪声, 并且MCC-${L_1}$-OMP方法有着更为优异的结果, 在不同脉冲噪声强度下均能够有更好的性能表现。

图  7  MSE随脉冲噪声强度变化曲线(GSNR=0)) (a)时延估计; (b)相对速度估计

下载: 全尺寸图片 幻灯片

(2) 不同广义信噪比下的仿真结果

本小节对比不同算法在不同广义信噪比下的性能表现, 仿真分别针对强脉冲噪声($\alpha = 1.2$), 中等脉冲噪声($\alpha = 1.5$), 弱脉冲噪声($\alpha = 1.8$)展开, 设置接收信号的广义信噪比从−5 dB增长到5 dB, 步长为1 dB, 每组实验同样进行1000次。图8给出了不同方法的实验结果。仿真结果纵向对比可以看出, 脉冲噪声强度越弱, 各种方法的估计性能越好; 同时, 随着GSNR的增加, 同等脉冲噪声出现频次下, 脉冲噪声异常值的强度逐渐降低, 因此三种方法的估计性能均有提升。但由于GSOMP没有对脉冲噪声进行针对性设计, 即使GSNR提升到较高的水平, 由于脉冲噪声的影响, 其估计结果仍处于失效状态。相比于RobOMP算法, 本文算法在GSNR=0 dB时, 基本能够实现对于时延参数的无误差估计, 同时对于相对速度的估计精度提升约一个数量级。

图  8  不同脉冲噪声强度下估计结果的MSE随广义信噪比(GSNR)变化曲线 (a) 时延估计, $\alpha = 1.2$; (b) 相对速度估计, $\alpha=1.2$; (c) 时延估计, $\alpha = 1.5$; (d) 相对速度估计, $\alpha=1.5$; (e) 时延估计, $\alpha = 1.8$; (f) 相对速度估计, $\alpha=1.8$

下载: 全尺寸图片 幻灯片

3.3   基于真实极地脉冲噪声的实验结果

基于图3(c)中的信道以及图1(a)中的实测噪声进行实测脉冲噪声下的仿真, 以验证所提算法的有效性, 仿真过程中选用的估计算法中的参数与前文相同。图9是无噪声接收信号和添加了实测北极冰下脉冲噪声信号的时域波形。由图9(b)可知, 脉冲噪声的能量远远超出了目标信号, 会严重影响时延–多普勒参数的估计精度。

图  9  信号时域波形 (a)无噪声接收信号时域波形; (b)添加实测北极冰裂噪声接收信号时域波形

下载: 全尺寸图片 幻灯片

图10(e)(f)是使用本文所提算法的时延–多普勒估计结果, 所提MCC-${L_1}$-OMP算法能够完全估计出五条多径的准确时延与多普勒信息, RobOMP算法在第五个多途路径上出现了严重错误, 并且在第一个多途路径上也有一定的偏差, GSOMP算法的性能则接近失效。

图  10  真实北极脉冲噪声下不同算法的估计结果 (a)(b) GSOMP方法; (c)(d) RobOMP方法; (e)(f) MCC-${L_1}$-OMP方法

下载: 全尺寸图片 幻灯片

为了验证所提算法的稳定性, 对图1(a)中所示的噪声以步长200 ms进行了切割, 共得到了200条脉冲噪声信号, 通过蒙特卡罗仿真实验, 得到了如表5所示的结果。从实验结果可以看出, 传统的GSOMP算法在真实脉冲噪声下失效, 本文所提算法能够有效抵抗脉冲噪声, 准确估计时延和多普勒信息, 参数估计准确性相较于RobOMP方法分别提升了58.8%和75.1%。

表  5  不同估计算法在实际脉冲噪声下的统计性能

方法	MSE-T	MSE-V
GSOMP	979.63	1.06
RobOMP	$2.3\times {10}^{-3}$	$8.15\times {10}^{-2}$
MCC-${L_1}$-OMP	$9.48\times {10}^{-4}$	$2.03\times {10}^{-2}$

下载: 导出CSV 

| 显示表格

4.   结论

为了提高极地脉冲噪声环境下水声信道时延–多普勒参数估计的准确性与稳健性, 本文提出了一种鲁棒的正交匹配追踪算法。该方法在传统正交匹配追踪框架内引入最大相关熵准则实现基底的准确挑选, 通过将信道估计问题重构为${L_1}$优化问题有效抑制脉冲噪声异常值对于参数求解的干扰, 并采用交替方向乘子法实现快速求解。数值模拟和基于真实北极冰下脉冲噪声的仿真结果表明, 与现有方法相比, 本文提出的MCC-${L_1}$-OMP算法具有更高的估计精度和更强的鲁棒性, 实现了最佳的均方误差性能。