Gradient Projection Algorithm for Total Variation Image Restoration by Adaptive Steplength Selection Rules
-
摘要: 针对图像去噪问题，本文基于全变差对偶公式提出一个新的梯度投影算法.算法采用改进的非单调线搜索和自适应BB（Barzilai-Borwein）步长，有效地改善了Chambolle梯度投影算法收敛慢的缺点.数值结果表明新算法优于一些已有的梯度投影算法.Abstract: We propose a new gradient projection algorithm for image denoising based on the dual of total variation. The new method exploits nonmonotone line-search and adaptive steplength selection based on strategies for alternation of the well-known Barzilai-Borwein rules. The proposed method is much faster than the Chambolle's gradient projection algorithm. Numerical results illustrate the efficiency of this method.
-
随着科技的巨大发展, 更多的产业开始实现半机械化甚至全机械化, 创造更多效益的同时却给人们带来了更多噪音困扰, 这些噪音已经严重影响了人们的身心健康, 特别是对老人和婴儿, 在某种程度上, 它已经成为了一种隐形的杀手.
大部分机械噪声都是由其内部的旋转机构引起的, 这种噪声的频率主要集中在低频段, 传统的被动噪声控制(Passive noise control, PNC)系统不能有效地对其进行抑制, 这时主动噪声控制(Active noise control, ANC)系统应运而生, 它对低频噪声作用效果显著[1]. ANC系统是一种基于声波相交干涉原理的电声学装置, 通过引入一个与待处理噪声(初级噪声源)幅度相同、频率相同、相位相反的噪声(次级噪声源), 使其达到降噪或者消噪的目的[2].
随着ANC系统的诞生, 大量的结构和算法也相应地被提出, 其中基于滤波-X最小均方(Filtered-X least mean square, FXLMS)算法的结构在ANC研究领域中最受欢迎而且也最有效[3-4]. ANC系统按照噪声的频谱特性可以分为窄带ANC (Narrowband ANC, NANC)和宽带ANC (Broadband ANC, BANC).这两种系统对其相应带宽的噪声消噪效果明显, 但是应对宽窄带混合的噪声时, 会出现``火花''现象, 造成系统不稳定(这里是利用BANC系统处理宽窄带混合噪声, 因为NANC系统不能处理宽带噪声).生活中常见的噪声基本都具有宽窄带混合特性, 例如风扇噪声、电机噪声以及切割机噪声等, 为此, 文献[5]中提出了一种混合ANC (Hybrid ANC, HANC)系统, 理想情况下可以完整分离混合噪声中的宽带分量和窄带分量, 然后分别利用BANC和NANC对其相应分量进行单独处理, 这种系统在处理宽窄带混合噪声上取得了巨大的进步, 同时也是处理宽窄带混合噪声最典型的结构.但是这种HANC系统主要存在两个问题：1)该系统中NANC子模块的输入频率是由非声学传感器(转速计、加速度计等)获取得到, 因此很容易与真实的频率出现偏差, 即频率失调(Frequency mismatch, FM), 这种现象会造成混合系统的不稳定, 使得噪声处理效果变差[6-7].文献[8-10]给出了FM的几种解决方法, 这些方法主要应对10%以内的FM，而且随着FM的增大, 系统收敛速度会变慢, 稳态误差也会随之增大, 当FM大于10%时, 这些方法不再适用; 2)该系统应对非平稳的混合噪声时效果不佳.
根据HANC系统处理宽窄带混合噪声的思想：预先分离出混合噪声中的宽带分量和窄带分量，然后对各分量进行独立的ANC处理, 本文提出了一种基于总体平均经验模态分解(Ensemble empirical mode decomposition, EEMD)的ANC系统.本文是第一次将EEMD算法应用到ANC领域, 并且提出了实时EEMD算法, 通过实时EEMD算法和ANC系统的有效结合, 自适应地分离出混合噪声中的宽带分量和窄带分量, 独立进行ANC处理后不仅可以消除平稳的混合噪声, 而且对非平稳的混合噪声也有明显的抑制效果.与HANC系统相比, 该系统没有FM的风险, 而且混合噪声的处理效果更好.
EEMD是经验模态分解(Empirical mode decomposition, EMD)的改进形式. EMD是1998年由Huang等提出的一种时频数据分析方法[11]，可以认为是一种自适应的二进制带通滤波器组(由多个带通滤波器组成, 这些带通滤波器的形状固定, 但是相邻两个滤波器的频率范围相互重叠, 而且一个是另外一个的一半或者二倍)[12]. EMD算法具有自适应多尺度分解特性, 利用信号本身的时间尺度特征, 通过一定的筛选条件将一个复杂信号中不同时间尺度的波动或趋势逐级分解为若干个固有模态函数(Intrinsic mode function, IMF)和一个余量[13-14], 每个分量的主要频带范围都不同, 因此该方法可以实现混合噪声中宽窄带分量的分离, 而且不需要估计窄带频率, 因此不受FM的影响. EMD的提出主要是为了解决非线性、非平稳问题, 可以实现对原始信号和数据的平稳化处理. EMD~经常和~Hilbert~变化结合使用, 通常称为希尔伯特黄变换(Hilbert-Huang transform, HHT).传统的EMD有时会有模态混叠(即没有合适的时间特征尺度对不同的模态成分进行有效地分离, 使得在一个固有模态函数中含有不只一个模态分量, 从而不能清楚地反映信号的内在固有性质).模态混叠现象一旦产生, 会对信号的后续分解产生影响, 使经过EMD分解后得到的结果不能完整地反映真实的物理过程[15-17].受Flandrin和Gledhill在噪声辅助分析方面研究的启发, 2009年, Huang等提出了总体平均经验模态分解(EEMD)[15], 该方法利用白噪声频谱的均匀分布, 使不同尺度的信号自动地映射到合适的参考尺度上, 可以得到合适的时间尺度对信号进行分解, 成功地避免模态混叠的影响. EEMD算法的应用条件和EMD算法相同, 但是需要的计算量更大[11].
通过仿真表明, 提出的系统在不引入FM的情况下也能很好分离出混合噪声中的宽带分量和窄带分量, 从而解决混合噪声问题, 同时由于EEMD对信号的平稳化处理(消除骑行波, 平滑不均匀幅度, 使其更加对称), 当噪声中出现非平稳变化时, 该系统也能保持较高的系统稳定性.
1. 传统宽窄带混合噪声控制系统
由宽窄带分量组成的原始混合噪声信号可以表示为:
$ \begin{equation} \begin{array}{l} {x_r}(n)= {x_n}(n)+ {x_w}(n)=\\ \;\;\;\;\;\;\;\;\;\, \sum\limits_{i = 1}^q {\left[ {{a_{r, i}}{x_{a_i}}(n)+ {b_{r, i}}{x_{b_i}}(n)} \right]} + {x_w}(n) \end{array} \end{equation} $
(1) $ \begin{equation} {x_{a_i}}(n)= \cos({\omega _i}n), {\rm{ }}{x_{b_i}}(n)= \sin({\omega _i}n) \end{equation} $
(2) 其中${x_n}(n)$是一个包含q种频率成分的窄带噪声信号, $\left\{ {{a_{r, i}}, {b_{r, i}}} \right\}_{i = 1}^q$是其傅里叶系数. ${x_w}(n)$表示宽带噪声信号, 一般是均值为零、方差为$\sigma _w^2$的白噪声或有色噪声.
为了解决传统BANC在处理宽窄带混合噪声时的不稳定现象, Xiao等提出了一种混合主动噪声控制(HANC)系统[5], 其框图如图 1所示, 主要包含NANC子系统、BANC子系统以及正弦噪声消除器(Sinusoidal noise canceller, SNC)子系统三个部分, 其中NANC子系统和BANC子系统分别针对噪声中的窄带分量和宽带分量, SNC子系统用来去除混合噪声中的窄带部分, 使进入BANC中的噪声都为宽带噪声. P(z)为初级路径的传递函数
$ \begin{equation} P(z)= \sum\limits_{t = 0}^{T - 1} {{p_t}{z^{ - t}}} \end{equation} $
(3) 真实的次级路径S(z)及其估计${\hat S}(z)$一般都用滑动平均(Moving average, MA)模型表示
$ \begin{equation} S(z)= \sum\limits_{m = 0}^{M - 1} {{s_m}{z^{ - m}}} \end{equation} $
(4) $ \begin{equation} {\hat S}(z)= \sum\limits_{m = 0}^{\hat M - 1} {{{\hat s}_m}{z^{ - m}}} \end{equation} $
(5) 其中$T-1, M-1, {\hat M}-1$和${p_t}, {s_m}, {{\hat s}_m}$分别为其阶数和系数.
混合噪声经过SNC子系统后, 理想情况下会变成一个纯净的宽带信号x(n), 然后利用BANC子系统对其进行处理, 产生次级噪声中的宽带分量yB(n). "Sync signal"通过实时的信号传感器(如加速度计)获得, 用来计算窄带信号的估计频率.次级噪声中的窄带分量yN(n)由NANC子系统生成.初级噪声p0(n)则由混合噪声经过初级通道P(z)直接产生, 此时误差信号表示为
$ \begin{equation} e(n)= {p_0}(n)- \sum\limits_{m = 0}^{M - 1} {{s_m}y(n-m)} + {v_p}(n) \end{equation} $
(6) $ \begin{equation} {{p_0}(n)= \sum\limits_{t = 0}^{T-1} {{p_t}} {x_r}(n - t)} \end{equation} $
(7) $ \begin{equation} y(n)= {y_B}(n)+ {y_N}(n) \end{equation} $
(8) 其中vp(n)为消噪点的背景噪声.
HANC系统虽然解决了混合噪声造成的系统不稳定现象, 但是同时也引入了FM问题, 给该系统带来了新的不稳定因素，FM的定义如下:
$ \begin{equation} \Delta {\omega _i} = \frac{{{\omega _i} - {{\hat \omega }_i}}}{{{\omega _i}}} \times 100\, \% \end{equation} $
(9) 其中${\omega _i}$为真实频率, ${\hat \omega _i}$为其估计值.
由于FM的影响, HANC系统中的SNC子系统不仅无法准确地去除混合噪声中的窄带分量, 而且还会引入新的窄带分量, 此时所有的窄带成分和噪声中的宽带成分最后都被送入到BANC子系统中同时进行处理, 由于混合噪声中窄带分量的增多, 使得整个系统变得更加不稳定, "火花"现象也更加严重.
2. 本文提出的系统
解决宽窄带混合噪声的关键是分离出其中的宽带分量和窄带分量, 并对它们单独进行处理, 才能避免相互之间的影响, 传统的HANC系统是一种精确的分离方式, 即将窄带分量和宽带分量进行完整分离, 然后对其单独处理.该方法虽然可以解决混合噪声问题, 但是却会引入FM现象, 严重影响ANC系统的处理效果.为此本文提出一种宽窄带分量模糊分离的方法, 利用EEMD算法对宽窄带混合噪声进行分解, 得到其所有IMF分量, 即将混合噪声按频率从高到低分解到各个IMF分量中, 此时宽带分量和窄带分量被分解到不同的IMF中.由于每个分量包含的是某一段频率的信号, 而不是某个确定频率信号, 所以称该分离方式为模糊分离.该方法不需要估计窄带分量的频率, 因此也不受这种估计误差带来的FM影响.
本文中提出的基于总体平均经验模态分解的主动噪声控制(EEMDANC)系统如图 2所示, 主要包含两个部分: EEMD算法和若干个BANC子系统(具体数量由被处理噪声的复杂程度决定), 利用EEMD算法将混合噪声分解成若干个分量, 然后通过BANC子系统对每个分量进行独立消噪.
2.1 实时EEMD算法处理
EEMD算法的基本原理为:给信号中加入高斯白噪声将信号和噪声组成一个``总体'', 加入信号的白噪声将遍布整个时频空间, 与被滤波器组分离的不同尺度分量相一致.当信号加在这些一致分布的白噪声背景上时, 不同尺度的信号, 自动地映射到合适的参考尺度上.而且由于零均值白噪声的特性, 加入的次数越多, 将这些多次分解的结果取平均后, 白噪声最终将被互相抵消而消除.总体平均的结果即被当作真实信号.加噪声的次数越多, 平均后的结果所含的剩余白噪声越小, 结果就越接近真实值. EEMD算法是EMD算法的改进形式, 同样假设任何信号都是由若干个IMF组成的, 每一个IMF必须满足下面两个条件[13-14]： 1)在整个信号上, 极值点的个数和过零点的个数相差不大于1; 2)在任意点处, 上下包络的均值为0.
EEMD算法的实现步骤如下:
1)对待分解的宽窄带混合噪声xr(n)重复加入白噪声NE次, 每次加入的白噪声都不同, 幅度用vi表示, 一般为0.1～0.4.通过在混合噪声中添加白噪声来平滑脉冲干扰, 使脉冲干扰在EMD分解过程中融入到白噪声模式中, 抑制脉冲干扰与其他固有模态之间的混叠.
$ \begin{equation} {x_i}(n)= {x_r}(n)+ {\nu _i}{\upsilon _i}(n) \end{equation} $
(10) 2)对每次加入白噪声后的混合噪声进行EMD分解.
a)首先确定加入白噪声后的混合噪声xi(n)的极大值和极小值, 然后通过三次样条插值函数拟合出极大值包络emax (n)和极小值包络emin(n).
b)计算整个包络线的均值, 用mi, 1(n)表示, 并得到xi(n)与mi, 1(n)的差值hi, 1(n), hi, 1(n)为去掉低频的新信号.
$ \begin{equation} {m_{i, 1}}(n)= \frac{{{e_{\max }}(n)+ {e_{\min }}(n)}}{2} \end{equation} $
(11) $ \begin{equation} {h_{i, 1}}(n)= {x_i}(n)- {m_{i, 1}}(n) \end{equation} $
(12) c) ${h_{i, 1}}(n)$一般不满足IMF定义的两个条件, 需要重复上述过程多次, 直到其满足IMF的定义, 此时得到第一个IMF分量${c_{i, 1}}(n)$.但是实际中, 上下包络的均值无法为零, Huang把标准偏差系数(Standard deviation, SD)作为停止准则, 即当SD满足下式时, 就认为包络的均值满足IMF的均值为零的条件[14].
$ \begin{equation} SD = \sum\limits_{n = 1}^{NT} {\left[{\frac{{{{\left| {{h_{i, k- 1}}(n)-{h_{i, k}}(n)} \right|}^2}}}{{{h_{i, k-1}}(n)}}} \right]} \le \varepsilon \end{equation} $
(13) 式中NT为待分解信号的总长度, ε为筛分门限, 一般取值在0.2到0.3之间.
d)用加入白噪声后的混合噪声信号减去${c_{i, 1}}(n)$, 得到一个去掉高频成分的余量${r_{i, 1}}(n)$.判断${r_{i, 1}}(n)$是否满足EMD的停止准则，如果满足则EMD结束, 否则将${r_{i, 1}}(n)$看成一个新的信号重复步骤a)～d). EMD停止准则的判断标准通常有两种:一种是当余量为常数或单调函数时, 从中不能再筛分出其他IMF时停止; 另一种是当最后一个IMF分量或余量的幅度小于预设值时停止.最后得到每组加入白噪声后信号的EMD结果为
$ \begin{equation} {x_i}(n)= \sum\limits_{j = 1}^N {{c_{i, j}}(n)} + {r_{i, N}}(n) \end{equation} $
(14) 其中N为IMF的总个数, ${r_{i, N}}(n)$是一个趋势项, 表示信号的平均趋势或均值.
3)将上述每次分解得到的IMF分量及其余量对应求平均, 最终得到的IMF分量和余量分别为
$ \begin{equation} {c_j}(n)= \frac{1}{{NE}}\sum\limits_{i = 1}^{NE} {{c_{i, j}}(n)} \end{equation} $
(15) $ \begin{equation} {r_N}(n)= \frac{1}{{NE}}\sum\limits_{i = 1}^{NE} {{r_{i, N}}(n)} \end{equation} $
(16) 此时, 就得到原始混合噪声${x_r}(n)$的EEMD结果
$ \begin{equation} {{x_r}(n)= \sum\limits_{j = 1}^N {{c_j}(n)} + {r_N}(n)} \end{equation} $
(17) 为保证信号分解过程中的实时性, 这里采用实时的EEMD算法, 即通过一个窗口W对输入信号从左到右进行实时分割, 分割得到的部分利用EEMD对其分解.同时为了确保分解得到的每组IMF具有连续性, 每次分割信号的窗口是连续移动的, 整个过程的示意图如图 3所示.
2.2 BANC子系统的FXLMS算法处理
利用实时EEMD对输入的混合噪声进行分段分解, 每段的长度都为W, 每段都可以被分解成N个IMF分量和一个余量, 再把得到的IMF分量和余量分别看成新的噪声信号, 利用FXLMS算法对其进行计算处理, 得到其各自的次级噪声.总的次级噪声为
$ \begin{equation} y(n)= \sum\limits_{j = 1}^{N + 1} {{y_j}(n)} \end{equation} $
(18) 令${{\pmb c}_{N+1}}(n) = {{\pmb r}_N}(n)$, 则
$ {y_j}(n) = {\pmb W}_j^{\rm{T}}(n){{\pmb c}_j}(n),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1 \le j \le N + 1 $
(19) 其中
$ \begin{equation} {{\pmb W}_j}(n) = [{w_{j, 0}}, \, {w_{j, 1}}, \, \cdots, \, {w_{j, {L_j} -1}}{{\rm{]}}^{\rm{T}}} \end{equation} $
(20) $ \begin{equation} {{\pmb c}_j}(n) = [{c_j}(n), \, {c_j}(n-1), \, \cdots , \, {c_j}{\rm{(}}n-W + 1{\rm{)}}{{\rm{]}}^{\rm{T}}} \end{equation} $
(21) 其中$\left\{ {{L_j}} \right\}_{j = 1}^{N + 1}$和$\left\{ {{w_{j, k}}} \right\}_{k = 0}^{{L_j} - 1}$分别为各控制器的长度和系数, 一般控制器的长度要小于分割窗口$W$.
利用FXLMS算法对$\left\{ {{{\pmb W}_j}(n)} \right\}_{j = 1}^{N + 1}$进行更新
$ \begin{equation} {{\pmb W}_j}(n + 1) = {{\pmb W}_j}(n){\rm{ + }}{\mu _j}(n){\hat S^{\rm{T}}}{{\pmb c}_j}(n)e(n) \end{equation} $
(22) $ \begin{equation} e(n)= p(n)- \sum\limits_{j = 1}^{N + 1} {{\pmb c}_j^{\rm{T}}(n)S{{\pmb W}_j}(n)} + {v_p}(n) \end{equation} $
(23) $${\hat S}$和${S}$都为非对称托普利兹矩阵. ${\hat S}$的第一列为$ {[{{\hat {\pmb s}}^{\rm{T}}}, {\bf 0}_{W-\hat M}^{\rm{T}}]^{\rm{T}}}$, 第一行为$[{\hat {\pmb s}_0}, {\bf 0}_{{L_j}-1}^{\rm{T}}]$. S的第一列为$ {[{{\pmb s}^{\rm{T}}}, {\bf 0}_{W-M}^{\rm{T}}]^{\rm{T}}}$, 第一行为$[{{\pmb s}_0}, {\bf 0}_{{L_j}-1}^{\rm{T}}]$.
$ \begin{equation} {\hat {\pmb s}} = {[{\hat s_0}, {\hat s_1}, \cdots, {\hat s_{\hat M -1}}]^{\rm{T}}} \end{equation} $
(24) $ \begin{equation} {\pmb s}= {[{s_0}, {s_1}, \cdots, {s_{M-1}}]^{\rm{T}}} \end{equation} $
(25) ${\mu _j}(n)$为一可变的步长, 用下面的表达式表示
$ \begin{equation} {\mu _j}(n)= {\mu _c}(j){\gamma _j}(n) \end{equation} $
(26) $ \begin{equation} {\gamma _j}(n)= \beta \left[{\frac{1}{{1 + {{\rm e}^{-\alpha \left| {{\delta _j}(n)e(n)} \right|}}}}-\xi } \right] \end{equation} $
(27) $ \begin{equation} {\delta _j}(n)= \frac{{\sum\limits_{i = 1}^W {{\left|{{c_{j, i}}{(n)}}\right|}^2} }}{{\sum\limits_{j = 1}^{N + 1} {\sum\limits_{i = 1}^W {{\left|{{c_{j, i}}{(n)}}\right|}^2} } }} \end{equation} $
(28) 式(27)中的$\alpha $是控制函数形状的常数, $\beta $是控制函数范围的常数, ${\delta _j}(n)$表示每一个分量在每段信号中占的能量比, ${\delta _j}(n)e(n)$则近似表示每个分量产生的误差.
2.3 系统稳定性分析
EEMDANC系统是由多个BANC子系统组成, 只有当每个子系统都稳定时, EEMDANC系统才能稳定.将式(23)代入式(22)，并对其两端取期望
$ \begin{equation} \begin{array}{l} {\rm{E}}\left\{ {{{\pmb W}_j}(n + 1)} \right\} = {\rm{E}}\left\{ {{{\pmb W}_j}(n)} \right\} + {\mu _j}(n){{\hat S}^{\rm{T}}}{\rm{E}}\left\{ {{{\pmb c}_j}(n)p(n)} \right\}-\\ {\mu _j}(n){\rm{E}}\left\{ {{{\hat S}^{\rm{T}}}\sum\limits_{i = 1}^{N + 1} {{{\pmb c}_j}(n){{\pmb c}_i^{\rm{T}}(n)}S{{\pmb W}_i}(n)} } \right\}+\\ {{\hat S}^{\rm{T}}}{\rm{E}}\left\{ {{{\pmb c}_j}(n){v_p}(n)} \right\} \end{array} \end{equation} $
(29) 输入信号向量与权重向量一般相互独立, 则
$ \begin{align} {\rm{E}}\left\{ {{{\hat S}^{\rm{T}}}\sum\limits_{i = 1}^{N + 1} {{{\pmb c}_j}(n){\pmb c}_i^{\rm{T}}(n)S{{\pmb W}_i}(n)} } \right\}=~~~~~~~~~~~~~~~~~~~~~~~~~~~\nonumber\\ {{\hat S}^{\rm{T}}}\sum\limits_{i = 1}^{N + 1} {{\rm{E}}\left\{ {{{\pmb c}_j}(n){\pmb c}_i^{\rm{T}}(n)} \right\}} S{\rm{E}}\left\{ {{{\pmb W}_i}(n)} \right\} \end{align} $
(30) 同时由于EMD分解得到的各分量之间相互正交且与${v_p}(n)$相互独立[11], ${\rm{E}}\left\{ {{v_p}(n)} \right\} = 0$, 式(29)可进一步化简为
$ \begin{equation} \begin{array}{l} {\rm{E}}\left\{ {{{\pmb W}_j}(n + 1)} \right\} = \left( {I - {\mu _j}(n){{\hat S}^{\rm{T}}}{R_j}S} \right){\rm{E}}\left\{ {{{\pmb W}_j}(n)} \right\}+\\ {\mu _j}(n){{\hat S}^{\rm{T}}}{{\pmb D}_j} \end{array} \end{equation} $
(31) 其中
$ \begin{equation} {{\pmb D}_j} = {\rm{E}}\left\{ {{{\pmb c}_j}(n)p(n)} \right\} \end{equation} $
(32) $ \begin{equation} {R_j} = {\rm{E}}\left\{ {{{\pmb c}_j}(n){{\pmb c}_j}^{\rm{T}}(n)} \right\} \end{equation} $
(33) 要想${\rm{E}}\left\{ {{{\pmb W}_j}(n)} \right\}$最终趋于稳定值, 则
$ \begin{equation} {W_{{j_\infty }}} = {\lim _{n \to \infty }}{\rm{E}}\left\{ {{{\pmb W}_j}(n)} \right\} = {\left( {{{\hat S}^{\rm{T}}}{{R}_j}S} \right)^{ - 1}}{{\hat S}^{\rm{T}}}{{\pmb D}_j} \end{equation} $
(34) 权重向量的误差定义为
$ \begin{equation} {{\pmb v}_j}(n)= {{\pmb W}_j}(n)- {{\pmb W}_{{j_\infty }}} \end{equation} $
(35) 将上式代入式(31), 进一步化简得
$ \begin{equation} {\rm{E}}\left\{ {{{\pmb v}_j}(n)} \right\} = {\left( {I - {\mu _j}(n){{\hat S}^{\rm{T}}}{R_j}S} \right)^n}{\rm{E}}\left\{ {{{\pmb v}_j}(0)} \right\} \end{equation} $
(36) 假设${\hat S} = S$, 则${{\hat S}^{\rm{T}}}{R_j}S = {S^{\rm{T}}}{R_j}S$为对称矩阵, 可以表示为
$ \begin{equation} {S^{\rm{T}}}{R_j}S = {Q_S}{\Lambda}{Q_S^{\rm{T}}} \end{equation} $
(37) 其中${\Lambda }$为${S^{\rm{T}}}{R_j}S$矩阵特征值组成的对角矩阵.将上式代入到式(36)中, 由${Q_S^{\rm{T}}}{Q_S} = {I}$化简得
$ \begin{equation} {\rm{E}}\left\{ {{{\pmb v}_j}(n)} \right\} = {Q_S}{\left( {I - {\mu _j}(n){\Lambda }} \right)^n}{Q_S^{\rm{T}}}{\rm{E}}\left\{ {{{\pmb v}_j}(0)} \right\} \end{equation} $
(38) 上式收敛的必要条件为, 即
$ \begin{equation} 0 < \frac{{{\mu _j}(n)}}{2} < \frac{1}{{{\lambda _{\max }}}} \end{equation} $
(39) 式中${\lambda _{\max }}$为矩阵${S^{\rm{T}}}{R_j}S$的最大特征值.
此时${{\pmb v}_j}(n)= \mathop {\lim }_{n \to \infty }{\rm{E}}\left\{ {{{\pmb v}_j}(n)} \right\} = {\bf 0}$, 从而
$ \begin{equation} {{\pmb W}_j}(n)={{\pmb W}_{{j_{opt} }}} = {\left( {{S^{\rm{T}}}{R_j}S} \right)^{ - 1}}{S^{\rm{T}}}{{\pmb D}_j} \end{equation} $
(40) 系统的代价函数为
$ \begin{equation} \begin{array}{l} J({W})= {\rm{E}}\left\{ {{e^2}(n)} \right\}=\\ \sigma _p^2 + \sigma _v^2 - 2\sum\limits_{j = 1}^{N + 1} {{{\pmb W}_j^{\rm{T}}}(n){S^{\rm{T}}}{{\pmb D}_j}}\, +\\ \sum\limits_{j = 1}^{N + 1} {{{\pmb W}_j^{\rm{T}}}(n){S^{\rm{T}}}{R_j}S{{\pmb W}_j}(n)} \end{array} \end{equation} $
(41) 其中$W = [{{\pmb W}_1}, {{\pmb W}_2}, \cdots, {{\pmb W}_{N + 1}}]$, $\sigma _p^2$是初级噪声$p(n)$的方差, $\sigma _v^2$是背景噪声${v_p}(n)$的方差.当${{\pmb W}_j}(n)$满足式(40)时, $J({W})$达到其最优值
$ \begin{equation} \begin{array}{l} J({W})={J_{\min }}=\\ \sigma _p^2 + \sigma _v^2 - {\pmb D}_j^{\rm{T}}\sum\limits_{j = 1}^{N + 1} {S{{\pmb W}_{{j_{opt} }}}} \end{array} \end{equation} $
(42) 此时EEMDANC系统达到稳定状态.
3. 仿真
在仿真1和仿真2中, 信号的采样频率都为2 000 Hz.初级路径$P(z)$和次级路径$S(z)$都是截止频率为$0.4\pi$ (400 Hz)的低通滤波器, 次级路径训练的信号为白噪声, 最终取其稳态时的1 000个点的均值作为次级路径估计的系数.次级路径估计时的参数设置如表 1所示.
表 1 次级路径估计参数设置Table 1 Parameters of secondary path estimation参数名称 参数值 初级路径$P(z)$的长度 ${\rm{40 (}}L{\rm{ = 41)}}$ 次级路径$S(z)$的长度 ${\rm{20 (}}M{\rm{ = 21)}}$ 次级路径估计${\hat S}(z)$的长度 ${\rm{30 ( }}\hat M{\rm{ = 31)}}$ 计算步长 0.001 训练输入信号的方差 $\sigma _{xs}^2(n)= 1$ 训练观测噪声的方差 $\sigma _{vs}^2(n)= 0.1$ 训练次数 10 000 3.1 仿真1:平稳变化的混合噪声
在该仿真过程中, 混合噪声中窄带分量的频率为$\omega = [{0.03\;}{0.06\;}0.09]\pi$ (因为采样率为2 000 Hz, 该频率也可表示为(30 Hz 60 Hz 90 Hz)), 其中各频率的傅里叶系数分别为${a_{r, i}} = [2\;1\;0.5], {b_{r, i}} = [{\rm{-}}1\;{\rm{-}}0.5\;0.1]\, (1 \le i \le 3)$, 宽带分量为一均值为0、方差为1的白噪声.分别利用HANC系统以及提出的EEMDANC系统对该混合噪声进行处理.
根据实际情况, HANC在窄带频率估计时会出现FM现象, 本次仿真中, FM分别为1%、5%、10%, 对应的BANC子系统控制器步长都为0.0002, 其他参数与文献[5]相同, 在本文提出的EEMDANC系统中, 实时EEMD算法每次分解的分量个数都为5 (4个IMF分量和1个余量), 分割信号的窗口长度为$W=101$, 各分量控制器的固定步长分别为${\mu _c}{\rm{ = }}\left[ {0.0004\;0.0009\;0.00009\;0.00006\;0.0003} \right]$, 步长中其他参数分别为$\alpha = 300, \beta = 0.4, \xi = 0.1$, 各分量主控制器的阶数都为51, 即$\left\{ {{L_j}} \right\}_{j = 1}^{N + 1} = 51$, 在该仿真条件下, 不同FM时HANC系统的消噪结果分别如图 4中的(a)、(b)、(c)所示.由于本文提出的EEMDANC系统没有频率估计, 不受FM的影响, 最终的消噪结果如图 4(d)所示.该次仿真的点数都为30 000, 独立运行次数为40.
图 4 平稳宽窄带混合噪声系统性能对比Fig. 4 Comparisons between systems$'$ performance for stationary mixture noise of wideband and narrowband3.2 仿真2:非平稳变化的混合噪声
分别利用HANC系统和提出的EEMDANC系统对非平稳的混合噪声进行消噪处理.仿真的参数设置为:总仿真点数为50 000, 独立运行次数为40次.在25 000点时混合噪声中窄带分量的幅值和频率同时发生突变:幅度由${a_{r, i}}, {b_{r, i}}$变为${\rm{ - }}{a_{r, i}}, {\rm{ - }}{b_{r, i}}$, 频率由$[0.03\;0.06\;0.09]\pi$ (30 Hz 60 Hz 90 Hz)变为$[0.05\;0.08\;0.12]\pi$ (50 Hz 80 Hz 120 Hz).宽带分量保持不变.在本次仿真中, HANC的BANC子系统主控制器的步长为0.00025, 提出的EEMDANC系统参数和HANC系统的其他参数设置与仿真1中一样.仿真结果如图 5所示.
图 5 非平稳宽窄带混合噪声的系统性能对比Fig. 5 Comparisons between systems$'$ performance for non-stationary mixture noise of wideband and narrowband3.3 仿真3:实际录制的混合噪声
在该仿真中, 被处理噪声为实时录制的某地下室某房间排风扇产生的噪音, 录制设备为NI公司的PXI-4472数据采集卡以及SKC公司的NP21麦克风.实际测量噪声的时域图和频域图分别如图 6(a)和图 6(b)所示, 信号的采样频率为22 kHz, 测得窄带分量的频率成分主要包含11 Hz ($0.001\pi$)和44Hz ($0.004\pi$). EEMDANC系统控制器的固定步长为${\mu _c} = \left[ {0.03\;0.01\;0.01\;0.00009\;0.001} \right]$, HANC系统中BANC子系统的控制器步长为0.0007, 其他参数设置与仿真1中一样.消噪后的结果分别如图 7和图 8所示, 其中图 7为时域消噪结果, 图 8为频域消噪结果.
3.4 仿真分析
1)当HANC系统引入FM后, 该系统变得不再稳定, "火花"现象又会重新出现, 而且变得更加严重, 随着FM增大, "火花"的个数也逐渐增多.出现这种现象的原因是真实频率和估计频率之间存在误差, 使得真正的窄带频率和估计的窄带频率都被送到了BANC子系统中进行处理, 由于窄带频率个数的增多, 导致"火花"现象更加严重.当FM减小时, 估计频率和真实频率接近, 产生的"火花"基本发生在同一处, 所以"火花"个数减少，随着FM增大, 真实频率和估计频率差距变大, 各自产生的"火花"现象也逐渐分离, "火花"个数也随之增多.为了便于分析EEMDANC系统, 这里给出平稳的混合噪声24950～25050点经过EEMD后时域图和频域图, 分别如图 9和图 10所示.从图 10中可以看出, 宽窄带混合噪声被分解成5个不同频率段的分量: IMF1的带宽为300 Hz以上; IMF2的带宽为100～600 Hz; IMF3的带宽为50～300 Hz; IMF4的带宽为20～150 Hz; 余量$r4(n)$的带宽为0～60 Hz.这也充分说明EEMD类似于一个自适应二进制带通滤波器组.仿真中窄带信号的频率分别为30 Hz 60 Hz和90 Hz, 因此经过EEMD后窄带分量集中在IMF4和余量$r4(n)$中, 宽带分量集中在IMF1～IMF3中, 通过这种分解, 同样有效分离了宽带分量和窄带分量, 避免了宽窄带混合噪声同时处理带来的"火花"现象, 而且该系统不需要估计窄带频率, 因此也不存在FM带来的影响.仿真结果图 4(d)也充分证明提出的EEMDANC系统在不需要估计窄带频率的条件下依然可以很好地处理宽窄带混合噪声.由于实际中FM现象很难避免, 而且随着非声学传感器的老化, 这种现象会更加严重, 所以提出的EEMDANC系统比HANC系统具有更好的实用性.
图 9 平稳宽窄带混合噪声的EEMD时域图Fig. 9 Time domain EEMD view of stationary mixture noise of wideband and narrowband
图 10 平稳宽窄带混合噪声的EEMD频域图Fig. 10 Frequency domain EEMD view of stationary mixture noise of wideband and narrowband2)通过观察图 5发现:当混合噪声发生非平稳变化时, HANC系统的收敛曲线发生迟滞现象, 不能及时跟踪控制, 需要经过重新调整才能再一次达到稳定状态, 而提出的EEMDANC系统却可以正常工作, 不受非平稳变化的影响. 图 11给出了非平稳变化混合噪声24950～25050点的EEMD图, 原始混合噪声在25000点时出现突变, EEMDANC系统将这种突变分解到4个IMF分量和一个余量中. EEMD算法主要有两个优点: 1)消除骑行波; 2)平滑不均匀幅度, 使其更加对称.这两个优点被称为EEMD的平稳化处理, 本文就是利用这种平稳化处理过程来解决ANC系统中的非平稳现象. 图 12为混合噪声的Hilbert谱(上图)和经过EEMD分解后各分量的Hilbert谱(下图), 从图中可以看出分解得到的各分量频率随着时间基本都是连续变化的, 而且每个分量的能量变化都比较均匀, 这说明EEMD分解后的各分量基本都趋于平稳, 然后对这些平稳分量进行独立的ANC处理, 就会减少甚至消除原始混合噪声非平稳变化带来的影响.
图 11 非平稳宽窄带混合噪声的EEMD图Fig. 11 EEMD view of non-stationary mixture noise of wideband and narrowband
图 12 混合噪声和EEMD后各分量的Hilbert谱Fig. 12 Hilbert spectrum of the mixture noise and the components after EEMD3)通过观察图 7和图 8, 提出的EEMDANC系统对排风扇噪声的消噪效果要明显优于HANC系统.这是由于实际排风扇噪声中窄带分量的频率复杂, 而传统的HANC系统不能同时对所有的频率成分进行跟踪控制, 只能对其主要频率进行作用, 而且部分窄带分量的频率会随着风扇的转速不均匀实时变化, 使得HANC系统不能对其频率做到精确估计.本文提出的EEMDANC系统是对信号频率的一种模糊分解, 不需要对频率进行精确估计, 因此也不受频率估计不完整和频率估计误差的影响, 同时这种信号分解方式也简化了宽带分量, 优化了系统对宽带噪声的消噪效果.
4. 总结
本文主要研究参考噪声信号具有宽窄带混合特性的ANC系统, 首先介绍了传统BANC系统以及Xiao等提出的混合ANC系统, 探讨了这两种方法处理混合噪声的可能性及其存在的风险.同时提出一种基于总体平均经验模态分解(EEMD)的ANC系统, 分析了其降低宽窄带混合噪声的机理.通过大量的仿真分析, 提出的EEMDANC系统不仅能有效消除宽窄带混合噪声, 而且也不会引入FM的影响, 同时当噪声出现非平稳变化时, 也能保证较高的系统稳态性.
-
图 2 Boat(1 024×1 024)梯度投影的相对范数vs. CPU时间(秒)
Fig. 2 Relative norm of projected gradient vs. CPU time (s) of Boat (1 024×1 024)
图 3 Cameraman( $Tol\, =\, 10^{-4}$ )和Barbara( $Tol\, =\, 10^{-6})$ MGPSSABB的去噪结果
Fig. 3 Denoising results of MGPSSABB of Cameraman ( $Tol\, =\, 10^{-4}$ ) and Barbara ( $Tol\, =\, 10^{-6}$ )
图 4 四个算法的对偶间隙vs. CPU时间(秒) (Tol=10−6)
Fig. 4 Duality gap vs. CPU time(s) for four algorithms (Tol=10−6)
表 1 数值结果
Table 1 Numerical results
Image Algorithm Niter Time(s) RMSE Shape(128×128) Chambolle 737 2.06 0.0739 Chambolle1 592 1.86 0.0738 MChambolle 196 0.81 0.0412 GPSSABB 182 0.70 0.0414 MGPSSABB 130 0.61 0.0414 Lena(256×256) Chambolle 352 5.79 0.0829 Chambolle1 272 4.73 0.0827 MChambolle 129 3.62 0.0221 GPSSABB 152 4.13 0.0214 MGPSSABB 85 2.34 0.0214 Peppers(512×512) Chambolle 389 35.80 0.0525 Chambolle1 298 26.27 0.0520 MChambolle 124 19.77 0.0520 GPSSABB 124 19.67 0.0079 MGPSSABB 92 15.79 0.0079 Boat(1024×1024) Chambolle 388 141.80 0.0078 Chambolle1 300 117.05 0.0078 MChambolle 140 101.66 0.0078 GPSSABB 156 107.90 0.0078 MGPSSABB 97 67.16 0.0078 
下载: 导出CSV
表 2 Cameraman (256 × 256)迭代步数和CPU时间(秒)
Table 2 Number of iterations (Iter) and CPU time (s) of Cameraman (256 × 256)
Algorithms Tol=10−2 Tol=10−3 Tol=10−4 Tol=10−6 Iter Time(s) Iter Time(s) Iter Time(s) Iter Time(s) Chambolle 27 0.45 163 3.18 822 16.04 14625 271.89 GPCL 32 0.59 114 1.98 549 10.19 10116 187.28 GPLS 18 0.53 132 5.38 607 26.22 12100 550.32 GPBBM 20 0.56 128 3.24 596 16.99 11032 295.31 GPBBM(3) 20 0.58 72 1.78 292 7.53 3186 79.15 SQPBBM 14 0.50 43 1.95 167 7.66 2532 105.18 GPBBsafe 16 0.33 47 1.05 165 4.23 2398 73.68 GPBBNM 16 0.34 48 1.06 162 3.42 2974 62.19 GPABB 16 0.45 47 1.23 179 4.68 2276 49.45 GPSSABB 13 0.41 48 1.36 146 4.31 1372 38.50 MGPSSABB 13 0.37 47 1.33 129 3.67 678 18.90
下载: 导出CSV
表 3 Barbara (512 × 512)迭代步数和CPU时间(秒)
Table 3 Number of iterations (Iter) and CPU time (s) of Barbara (512 × 512)
Algorithms Tol=10−2 Tol=10−3 Tol=10−4 Tol=10−6 Iter Time(s) Iter Time(s) Iter Time(s) Iter Time(s) Chambolle 27 3.04 131 14.27 536 60.47 8583 939.73 GPCL 24 2.78 79 9.31 331 37.10 5854 641.63 GPLS 36 6.61 90 20.90 319 77.84 5532 1316.20 GPBBM 20 3.03 84 13.32 340 52.10 6096 872.62 GPBBM(3) 20 2.89 74 11.31 229 34.94 2588 365.88 SQPBBM 14 2.75 33 7.71 106 23.20 1530 325.25 GPBBsafe 15 2.09 40 5.76 118 18.75 1494 290.38 GPBBNM 15 2.42 41 5.55 116 15.69 1467 194.15 GPABB 14 2.43 35 5.91 125 21.23 1123 180.51 GPSSABB 12 2.03 34 5.97 104 19.25 878 159.95 MGPSSABB 14 2.39 37 6.26 90 15.55 431 86.31
下载: 导出CSV
-
[1] Rudin L I, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena, 1992, 60(1-4): 259-268 doi: 10.1016/0167-2789(92)90242-F [2] Aubert G, Kornprobst P. Mathematical Problems in Image Processing. Berlin: Springer, 2002. [3] Chan T, Shen J H. Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods. Philadelphia: SIAM, 2005. [4] 王诗言, 于慧敏.基于全变分的运动分割模型及分裂Bregman算法.自动化学报, 2015, 41(2): 396-404 http://www.aas.net.cn/CN/abstract/abstract18618.shtmlWang Shi-Yan, Yu Hui-Min. Motion segmentation model based on total variation and split Bregman algorithm. Acta Automatica Sinica, 2015, 41(2): 396-404 http://www.aas.net.cn/CN/abstract/abstract18618.shtml [5] 张瑞, 冯象初, 王斯琪, 常莉红.基于稀疏梯度场的非局部图像去噪算法.自动化学报, 2015, 41(9): 1542-1552 http://www.aas.net.cn/CN/abstract/abstract18729.shtmlZhang Rui, Feng Xiang-Chu, Wang Si-Qi, Chang Li-Hong. A sparse gradients field based image denoising algorithm via non-local means. Acta Automatica Sinica, 2015, 41(9): 1542-1552 http://www.aas.net.cn/CN/abstract/abstract18729.shtml [6] Chan T F, Golub G H, Mulet P. A nonlinear primal-dual method for total variation-based image restoration. SIAM Journal on Scientific Computing, 1999, 20(6): 1964-1977 doi: 10.1137/S1064827596299767 [7] Chambolle A. An algorithm for total variation minimization and applications. Journal of Mathematical Imaging and Vision, 2004, 20(1): 89-97 http://cn.bing.com/academic/profile?id=2144275666&encoded=0&v=paper_preview&mkt=zh-cn [8] Chambolle A. Total variation minimization and a class of binary MRF models. In: Proceedings of the 5th International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition. St. Augustine, FL, USA: Springer, 2005. 136-152 [9] 任福全, 邱天爽.基于二阶广义全变差正则项的模糊图像恢复算法.自动化学报, 2015, 41(6): 1166-1172 http://www.aas.net.cn/CN/abstract/abstract18691.shtmlRen Fu-Quan, Qiu Tian-Shuang. Blurred image restoration method based on second-order total generalized variation regularization. Acta Automatica Sinica, 2015, 41(6): 1166-1172 http://www.aas.net.cn/CN/abstract/abstract18691.shtml [10] Zhang B X, Zhu Z B, Li S A. A modified spectral conjugate gradient projection algorithm for total variation image restoration. Applied Mathematics Letters, 2014, 27: 26-35 doi: 10.1016/j.aml.2013.08.006 [11] Xiao Y H, Wu S Y, Qi L Q. Nonmonotone Barzilai-Borwein gradient algorithm for l1 regularized nonsmooth minimization in compressive sensing. Journal of Scientific Computing, 2014, 61(1): 17-41 doi: 10.1007/s10915-013-9815-8 [12] Loris I, Bertero M, De Mol C, Zanella R, Zanni L. Accelerating gradient projection methods for l1-constrained signal recovery by steplength selection rules. Applied and Computational Harmonic Analysis, 2009, 27(2): 247-254 doi: 10.1016/j.acha.2009.02.003 [13] Huang Y K, Liu H W, Zhou S. An efficient monotone projected Barzilai-Borwein method for nonnegative matrix factorization. Applied Mathematics Letters, 2015, 45: 12-17 doi: 10.1016/j.aml.2015.01.003 [14] Xiao Y H, Jin Z F. An alternating direction method for linear-constrained matrix nuclear norm minimization. Numerical Linear Algebra with Applications, 2012, 19(3): 541-554 doi: 10.1002/nla.v19.3 [15] Yu G H, Qi L Q, Dai Y H. On nonmonotone Chambolle gradient projection algorithms for total variation image restoration. Journal of Mathematical Imaging and Vision, 2009, 35(2): 143-154 doi: 10.1007/s10851-009-0160-3 [16] Dai Y H, Fletcher R. Projected Barzilai-Borwein methods for large-scale box-constrained quadratic programming. Numerische Mathematik, 2005, 100(1): 21-47 doi: 10.1007/s00211-004-0569-y [17] Bertsekas D P. Nonlinear Programming. Nashua: Athena Scientific, 1999. [18] Frassoldati G, Zanghirati G, Zanni L. New adaptive stepsize selections in gradient methods. Journal of Industrial and Management Optimization, 2008, 4(2): 299-312 doi: 10.3934/jimo [19] Barzilai J, Borwein J M. Two-point step size gradient methods. IMA Journal of Numerical Analysis, 1988, 8(1): 141-148 doi: 10.1093/imanum/8.1.141 [20] Birgin E G, Martínez J M, Raydan M. Nonmonotone spectral projected gradient methods on convex sets. SIAM Journal on Optimization, 2000, 10(4): 1196-1211 doi: 10.1137/S1052623497330963 [21] Grippo L, Lampariello F, Licidi S. A nonmonotone line search technique for Newton's method. SIAM Journal on Numerical Analysis, 1986, 23(4): 707-716 doi: 10.1137/0723046 [22] Zhu M Q, Wright S J, Chan T F. Duality-based algorithms for total-variation image restoration. Computational Optimization and Applications, 2010, 47(3): 377-400 doi: 10.1007/s10589-008-9225-2 [23] Ekeland I, Témam R. Convex Analysis and Variational Problems. Philadelphia: SIAM, 1999. 期刊类型引用(2)
1. 张绍荣，朱志斌，冯宝，余天佑，李智. 基于组稀疏贝叶斯逻辑回归运动想象脑电信号分类模型的通道选择与分类新算法. 仪器仪表学报. 2019(10): 179-191 . 
百度学术2. 石曼曼，李雷. 一种改进的OMP算法及其在图像重构上的应用. 计算机技术与发展. 2018(02): 94-97 . 百度学术其他类型引用(3)
-
下载:
计量
- 文章访问数: 2287
- HTML全文浏览量: 517
- PDF下载量: 1724
- 被引次数: 5
