1.   引言

分布式自适应估计是一种多节点协作的信息处理方式，即传感器网络中的各个节点均通过自适应迭代方式参与计算，并按照特定的协作策略与邻居节点进行信息交互，从而实现对感兴趣参数的有效估计。在各种节点协作策略中扩散策略^[1]更具灵活性和适应性，适于实现大规模网络参数的自适应估计。因而，在救灾管理、精确农业、电力系统建设等众多领域^[2-4]得到了广泛应用。

根据输入信号的特性不同，扩散式自适应估计算法可分为白输入信号算法和有色输入信号算法。最早提出的扩散式自适应估计算法—DLMS(Diffusion Least Mean Square)及其改进算法^[5-7]就是在白输入信号的假设之下得到的，对于有色输入信号该类算法的性能退化严重。为此，各种针对有色输入信号的算法相继出现，扩散式仿射投影算法(Diffusion Affine Projection Algorithm, DAPA)^[8,9]是其中之一，它是在DLMS的基础上通过对输入数据的重用来保证在输入信号存在相关性时仍快速收敛。为了便于算法分析，DAPA算法简单地假设系统加性噪声为符合大数定理的高斯分布。但是，实际应用中往往会遇到由大气、同信道干扰、节点失效等造成的具有脉冲特性的非高斯噪声^[10-16]，这类噪声在幅度上的强脉冲性会导致传统基于高斯假设的2阶统计量算法发生性能退化甚至完全失效，而DAPA算法恰恰是基于2阶统计量的，故其不适用于非高斯噪声。抑制非高斯噪声的方法大致可分为两类：信息论法和稳健统计量法。信息论法依据熵是局部相似性的度量而采用最大相关熵准则(Maximum Correntropy Criterion, MCC)^[12,13]和最小误差熵^[14]抑制具有脉冲性的异常值即非高斯噪声。稳健统计量法是目前主流的非高斯噪声抑制方法，用于扩散式自适应估计的稳健统计量主要是M估计量和L1范数，该类方法通过求取某个非2阶统计函数的最小化来获得最优估计^[15-19]。实际上，在稳健估计理论中还有另一种对抗异常值非常有效的估计量：R(Rank based estimator)估计量，其代价函数是通过对误差值进行排序的评分函数来定义的。R估计量的优点是计算量小且抗异常值能力强^[20]，其中的Wilcoxon范数近年来得到了信号处理领域的广泛关注，目前的研究主要集中在单节点的参数估计上：Majhi等人^[21]将Wilcoxon范数作为代价函数提出Wilcoxon最小均方误差算法(Wilcoxon Least Mean Square, WLMS)，实现了稳健的系统辨识算法；Ban等人^[22]提出了基于Wilcoxon范数的仿射投影算法(Wilcoxon Affine Projection Algorithm, WAPA)，通过最小化加权Wilcoxon范数来克服APA算法的不足；文献[23]提出符号回归Wilcoxon LMS算法(Sign-regressor Wilcoxon LMS)，获得了更快的收敛速度；文献[24]在块最小二乘算法(Block Recursive Least Square, BRLS)中加入Wilcoxon范数，通过QR分解减少其计算量。Wilcoxon范数在分布式自适应网络的应用起步较晚，最早出现的是由Kumar等人^[25]于2016年提出的扩散最小Wilcoxon范数(Diffusion Minimum Wilcoxon Norm, DMWN)算法，该文分别在不同非高斯条件下对算法进行了验证并给出了理论分析。

受上述文献启发，考虑到非高斯噪声对DAPA类算法的影响，本文基于Wilcoxon范数对原代价函数进行了改进，得到了新的迭代方程，同时为提升算法的适应能力，基于最优步长的上界提出了变步长策略，即通过建立步长与误差信号之间的非线性关系打破固定步长的局限性，在收敛初始阶段，采用大步长加快收敛速度；当接近收敛阶段，保持步长不变。仿真实验表明，本文新算法在非高斯噪声和高斯噪声环境下均表现出良好的性能。

文中用到的各种运算及其含义为：${\rm{col}}\{ {{{x}}_1}\; {{{x}}_2} ···\; {{{x}}_N}\}$表示将列向量${{{x}}_1} \;{{{x}}_2}\;···\; {{{x}}_N}$依次纵向组成一个列向量，${[ \cdot ]^{\rm{T}}}$表示求矩阵的转置，${\rm{sign}} ( \cdot )$表示符号运算，${\rm{E}}[ \cdot ]$表示求期望，${\rm{diag}} \{ \cdot \}$表示产生对角矩阵，${\lambda _{\max }}( \cdot )$为矩阵的最大特征值，${{{I}}_{NM}}$为$N \times M$的单位矩阵，${{X}} \otimes {{C}}$表示求两矩阵X和C的克罗内克积，$|| \cdot ||_2^2$为向量的欧氏范数的平方，$| \cdot |$表示求绝对值，${\rm{med}}( \cdot )$表示求中值，$\triangleq$表示定义。

2.   问题描述

2.1   系统模型

考虑一个由$N$个节点组成的分布式网络，所有与节点$k$直接相连的邻居节点集合记为${N_k}$(包括节点$k$自身在内)，其中的任意两个邻居节点$l,k$通过融合系数${c_{lk}}$实现信息交互。网络中待估计的未知参数向量${{{w}}_0}$是$M \times 1$维向量，${{{u}}_k}(i)$是$M \times 1$的输入向量，则节点$k$在$i$时刻的局部测量值为$\left\{ {{d_k}(i),{{{u}}_k}(i)} \right\}$, $k = 1,2, ··· ,N$，它们之间的线性关系为

其中，$M$为滤波器长度，${{{u}}_k}(i) = [{u_k}(i)\ {u_k}(i - 1)\ ··· {u_k}(i - M + 1)]^{\rm{T}}$是$M \times 1$维的输入向量。假设${{{u}}_k}(i)$和${\eta _k}(i)$相互独立，${\eta _k}(i)$为背景噪声，在实际场景中它可能是符合高斯分布的高斯噪声，也可能是不符合高斯分布的非高斯脉冲噪声。本文采用混合高斯模型^[13,16,19]来描述非高斯噪声：

其中，${\chi _k}(i)$与${v_k}(i)$是独立同分布的零均值高斯白噪声，方差分别为$\sigma _{\chi ,k}^2$和$\sigma _{v,k}^2$($\sigma _{v,k}^2 \gg \sigma _{\chi ,k}^2$), ${m_k}(i) = {b_k}(i) {v_k}(i)$表示非高斯噪声，${b_k}(i)$为伯努利过程，其概率密度函数为$P[{b_k}(i) = 1] = {p_k}$, $P[{b_k}(i) = 0] = 1 - {p_k}$, ${p_k}$表示在节点$k$脉冲性噪声发生的概率。相应地，${\eta _k}(i)$的概率密度函数为

其中，$\sigma _{\eta ,k}^2 = {p_k}\sigma _{s,k}^2 + (1 - {p_k})\sigma _{\chi ,k}^2$, $\sigma _{s,k}^2 = \sigma _{\chi ,k}^2 + \sigma _{v,k}^2$, $\sigma _{m,k}^2 = {p_k}\sigma _{v,k}^2$, ${\zeta _{im}} = \dfrac{{\sigma _{m,k}^2}}{{\sigma _{\chi ,k}^2}}$表示非高斯噪声与高斯噪声在${\eta _k}(i)$中所占比例。本文算法的目标就是在这种非高斯噪声模型下，利用$\left\{ {{d_k}(i),{{{u}}_k}(i)} \right\}$估计${{{w}}_0}$。

2.2   传统的DAPA算法

扩散式策略有两种实现形式，一是先融合再进行自适应的CTA(Combine-Then-Adapt)模式，二是先进行自适应再融合的ATC(Adapt -Then- Combine)模式。二者的基本结构相同，但融合和自适应的顺序不同导致ATC模式比CTA模式能够更快速地遍历所有的节点，因此本文采用ATC模式^[1]。

ATC模式下扩散式自适应估计的节点$k = 1,2,···,N$的迭代过程为

式中，$f[{{{\varPhi}} }_{l}(i+1); l\in {N}_{k}]$表示节点$k$与其邻居节点$l$间的融合函数，${N_k}$表示节点$k$的所有邻居节点的集合；$\nabla {J_k}({{{w}}_k}(i)$是代价函数的梯度，根据代价函数的不同可以得到不同算法的迭代过程。

DAPA算法在DLMS算法的基础上对输入数据进行重组，将P(称为仿射投影阶数)个输入向量组合形成输入信号的$M \times P$维矩阵形式：${{{U}}_k}(i) = \left[ {{{{u}}_k}(i)\ {{{u}}_k}(i - 1)\ ··· \ {{{u}}_k}(i - P + 1)} \right]$；期望信号表示为${{{d}}_k}(i) \;=\; {{w}}_0^{\rm{T}}{{{U}}_k}(i) + {{{\eta}} _k}(i) = [{d_k}(i)\ {d_k}(i - 1)\ ··· \ {d_k} (i - P + 1)]$；噪声向量表示为${{{\eta}} _k}(i) = [{\eta _k}(i)\ {\eta _k}(i - 1) ··· \ {\eta _k}(i - P + 1)]$。传统的DAPA 算法在噪声符合高斯分布的假设下，通过求解式(5)的最优化问题来估计权向量${{{w}}_0}$。

由此可以看到，DAPA算法是一种2阶统计量方法，在非高斯噪声下必将失效。

3.   本文提出的算法

3.1   迭代方程的推导

Wilcoxon范数^[24]的定义通过评分函数$\varphi (u)$描述，$\varphi (u)$是一个递增函数，并满足$\displaystyle\int_0^1 {{\varphi ^2}(u) {\rm{d}}u < \infty }$，其自变量$u \in [0,1]$是实数。对于维数为$L$的向量${{v}} = {[{v_1}\ {v_2}\ ···\ {v_L}]^{\rm{T}}} \in {R^L}$，其Wilcoxon范数定义为

其中，${v_i}$表示向量${{v}}$的第$i$个分量，$R({v_i})$表示${v_i}$在新构成的升序排列的向量中对应的位置序号。为描述简便，定义$a({v_i}) = \varphi \left[\dfrac{{R({v_i})}}{{L + 1}}\right]$。

考虑到Wilcoxon的诸多优点，本文基于Wilcoxon范数推导新的迭代方程。后验误差向量：

由于冲击噪声样本对算法性能的影响发生在迭代过程中，因此为提高算法对非高斯噪声的鲁棒性，增加一个约束条件，从而最小化条件约束公式为

采用拉格朗日算子法求解上述优化问题，有

对式(9)求关于${{{w}}_k}(i)$的偏导数并令其为零，同时考虑式(6)的Wilcoxon范数定义，得到

对式(8)的约束条件取等号，并将式(10)代入，得到

继而可得

$\delta$为很小的数。

同时，考虑到Wilcoxon范数中评分函数的范围，本文将原来的评分函数：$\varphi (u) = \sqrt {12} (u - 0.5)$替换为$\varphi (u) = {\rm{sign}} (u - 0.5)$，从而得到本文的迭代方程为

其中，${c_{lk}}$称为融合参数，其取值满足条件$\displaystyle\sum\nolimits_{l \in {N_k}} {{c_{lk}}} = 1$。

为进一步提高算法的适应性，本文采用步长可变策略，具体的推导过程见3.2节。

3.2   变步长的推导

由式(13)定义中间估计误差${\tilde {{\varPhi}} _k}(i) = {{{w}}_0} - {{{\varPhi}} _k}(i)$，这里假设待估计向量${{{w}}_0}$是时不变的，则有

所以

这里，$\Delta ({\mu }_{k}(i))\triangleq -2{\mu }_{k}(i)$${\rm{E}}\left\{\dfrac{a[{{e}}_{k}^{{\rm{T}}}(i)]{{U}}_{k}^{{\rm{T}}}(i){\tilde{{{\varPhi}} }}_{k}(i)}{\sqrt{{\{{{U}}_{k}(i)a[{{e}}_{k}(i)]\}}^{{\rm{T}}}\{{{U}}_{k}(i)a[{{e}}_{k}(i)]\}+\delta }}\right\}+{\mu }_{k}{}^{2}(i)$是关于${\mu _k}(i)$的函数。为使式(15)所示的从$i$时刻到$i + 1$时刻的均方差达到最小，$\Delta ({\mu _k}(i))$必须进行最小化计算。将$\Delta ({\mu _k}(i))$重写为

式(16)中第1项分子部分的$a[{{e}}_k^{\rm{T}}(i)]{{{\eta}} _k}(i)$直接包含噪声成分，因此对其进一步讨论：

将$||{{{\eta}} _k}(i)|{|_1}$近似为其期望：${\rm{E}}[||{{{\eta}} _k}(i)|{|_1}]$，并考虑$\left| {{\eta _k}(i)} \right|$可看作为半正态分布^[26]，得到

因此，

将式(19)的右边定义为关于${\mu _k}(i)$的函数，并对其求关于${\mu _k}(i)$的导数，再令其为零，同时考虑到期望难以得到，首先设

式(20)中的${\sigma _{{\eta _k}}}$可由文献[26]的方法得到，然后采用滑动平均方法获得可变步长：

$\gamma$($0 \!<\! \gamma \!<\! 1$)是滑动系数，变步长的初值${\mu _k}(0){\rm{ = }}\sqrt {\dfrac{{\sigma _{{d_k}}^2}}{{\sigma _{{u_k}}^2N}}}$, $\sigma _{{d_k}}^2$和$\sigma _{{u_k}}^2$分别是输入信号和输出信号的功率。非高斯的冲激噪声样本是以小概率随机出现的，且具有很大的幅度，也就是说，强脉冲性的非高斯噪声样本不是经常出现的，同时考虑到在初始阶段${\mu _k}(i)$比较大，所以步长选择为${\beta _k}(i)$，从而步长会根据$\left| {a[{{e}}_k^{\rm{T}}(i - 1)]{{{e}}_k}(i - 1)} \right|$的变化而更新，而在接近稳定阶段就保持步长不变。当非高斯噪声出现时，步长保持不变，以避免脉冲噪声样本对步长参数的负面影响。

4.   实验仿真

为验证算法的有效性，这里将本文算法、DLMS^[5], DAPA^[8], DSELMS^[7], DMCC^[13], DLMP^[6]算法在系统辨识问题中所获的结果进行了对比分析。待估计的系统为${{{w}}_0}{\rm{ = rand(}}M,1)/{\rm{norm}} [{\rm{rand}}(M{\rm{,1)]}}$, ${\rm{rand(}} \cdot {\rm{)}}$表示标准的均匀分布函数，$M{\rm{ = }}20$，系统所在的分布式网络由20个节点组成，融合参数由Metropolis准则^[1]获取。本文分别对白信号和有色信号两种输入进行了仿真：白信号是均值为零，方差为$\sigma _g^2$的高斯白噪声；有色信号由上述白信号经过1阶AR系统$G(z) = \dfrac{1}{{(1 + 0.9{z^{ - 1}})}}$得到。仿真中加性噪声分别为：高斯噪声和非高斯噪声，高斯噪声为${\chi _k}(i)$，其方差为$\sigma _{\chi ,k}^2$; ${v_k}(i)$是与${\chi _k}(i)$独立同分布的零均值高斯白噪声，其方差为$\sigma _{v,k}^2$，二者相加组成非高斯噪声${\eta _k}(i)$，如式(2)所示。仿真中取${\zeta _{im}} = 10000$, ${p_k}{\rm{ = }}0.01$。同时，定义信噪比为

其中，$\sigma _{{d_0}{,_k}}^2$为${{u}}_k^{\rm{T}}(i){{{w}}_0}$的方差。设所有节点的信噪比均为20dB，其它参数设置为：$\gamma = 0.99$; $P = 4$; DMCC的核宽取2, DLMP中的p=1.8。文中采用NMSD(Network Mean Square Deviation) 作为算法收敛性能的衡量标准：

NMSD的值越小意味着所估计的向量越逼近未知系统，所有结果均为50次独立平均的结果。

4.1   输入为白信号

(1) 背景噪声为高斯噪声

图1所示为输入为白高斯信号，背景噪声为高斯噪声时各算法的性能，DLMS, DSELMS, DMCC, DLMP, DAPA算法的步长分别为0.0038, 0.0066, 0.2, 0.0072, 0.021。由图1可以看到，在这种条件下本文算法仍可以收敛，但是性能不如其他算法，原因是本文算法是针对有色输入和非高斯噪声而提出的。

图 1  输入为白高斯信号，各算法在高斯噪声下的NMSD曲线

下载: 全尺寸图片 幻灯片

(2)背景噪声为非高斯噪声

图2所示为在非高斯噪声下这7种算法的NMSD曲线。实验中设置DLMS, DSELMS, DMCC, DLMP, DAPA算法的步长分别为0.0038, 0.0066, 0.2, 0.0068, 0.008。可以看到，DLMS和DAPA算法已经失效，其它几种算法可以在非高斯噪声下较好工作，本文算法由于是针对相关输入信号而提出的，因此性能不如DSELMS和DMCC， 但与DSELMS相差不多，且远远优于DLMP。

图 2  输入为白高斯信号，各算法在非高斯噪声下的NMSD曲线

下载: 全尺寸图片 幻灯片

4.2   输入为有色信号

(1) 背景噪声为高斯噪声

图3所示为输入信号为有色信号，背景噪声为高斯噪声下各算法的性能曲线。DLMS, DSELMS, DMCC, DLMP, DAPA算法的步长分别为0.0023, 0.0076, 0.56, 0.0072, 0.008。可以看到，尽管背景噪声是高斯噪声，但DLMS算法由于并未考虑输入信号的相关性而导致性能较差，而DSELMS, DLMP, DMCC算法均是针对非高斯噪声而提出的，且没有输入信号的相关性，因此性能也不理想。DAPA算法采用了信号重用策略减少了信号相关性的影响，因此在高斯噪声下获得的性能最好。但是，同时也看到，本文算法优于DSELMS, DLMP及DMCC算法，与DAPA算法几乎相同，这说明本文算法在高斯噪声条件下也可很好地工作。

(2) 背景噪声为非高斯噪声

输入信号为有色信号，背景噪声为非高斯噪声，并为进一步说明本文算法优越性，考虑了系统突变的情况，即在$i = 10000$时，待估计系统由${{{w}}_0}$突变为${\rm{ - }}{{{w}}_0}$，各算法的跟踪性能如图4所示。实验中设置DLMS, DSELMS, DMCC, DLMP, DAPA算法的步长参数分别为0.0023, 0.0056, 0.75, 0.0071, 0.008。可以看到，DLMS算法和DAPA算法在非高斯噪声下失调，不能跟踪系统的变化，其他算法能够正常收敛，而本文算法在信道改变之后仍保持收敛，且稳态误差最小即抑制非高斯噪声的能力最强。

图 3  输入为有色信号，各算法在高斯噪声下的NMSD曲线

下载: 全尺寸图片 幻灯片

图 4  输入为有色信号，各算法在非高斯噪声下的跟踪性能

下载: 全尺寸图片 幻灯片

5.   结论

本文基于Wilcoxon范数改进了DAPA算法的代价函数，并考虑固定步长的局限，提出了具有较强鲁棒性的扩散式自适应参数估计新算法。仿真结果表明，本文的新算法不仅能有效抑制非高斯噪声的干扰，而且在高斯噪声下也能获得良好性能。在非高斯噪声下，与现有一些针对非高斯噪声而提出的方法相比，本文算法具有更好的性能。需要指出的是，本文算法的复杂度与文中所对比的算法要略高一些。