一种基于粒子群优化的心电信号自适应非局部均值降噪方法

摘要

本发明公开了基于粒子群优化的心电信号自适应非局部均值降噪方法，包括以下步骤：(1)对原始心电信号进行非局部均值滤波，得到滤波结果；(2)采用粒子群算法对每一信号点i自适应确定非局部均值滤波方法中的两个关键参数，即衰减参数和搜索窗口；(3)利用自适应确定的上述两个参数，对原始心电信号进行非局部均值滤波，获得最终的降噪信号。

说明书

技术领域

本发明属于医学信号降噪领域，更具体地，涉及一种医学超声图 像降噪方法。

背景技术

心电信号采集自病人体表，是一种无创性的检测手段。因此，心 电信号采集过程中受到基线漂移、工频干扰、肌电信号等多种噪声的 干扰，影响后续临床诊断的准确性。因此这些噪声能够被适当的去除， 对心电信号诊断具有决定性的意义。

现有的心电信号滤波算法主要包括有限冲击响应(Finite ImpulseResponse,FIR)滤波、无限冲击响应(InfiniteImpulse Response,IIR)滤波、中值滤波、小波滤波和非局部均值降噪等。其 中，非局部均值方法因充分利用了信号的冗余性，较其它方法能更好 地抑制噪声的不利影响。但传统非局部均值降噪方法，对不同的心电 信号波段，都采用相同的衰减系数，不能较好地适应心电信号具有明 显周期性和区域性的特点，导致难以兼顾均匀区域的平滑和细节信息 的保护。

发明内容

针对现有技术中存在的问题，本申请提供的是一种基于粒子群优 化的心电信号自适应非局部均值降噪方法，简称ANLM。其中通过对 衰减参数和搜索窗口进行研究和设计，与现有降噪方法相比，ANLM 具有更好的恢复性能。

为实现上述目的，按照本发明的一个方面，提供了一种基于粒子 群优化算法的心电信号自适应非局部均值降噪方法,该方法包括以下 步骤：

步骤1对原始心电信号进行非局部均值滤波，得到预滤波结果。

对原始信号O进行非局部均值滤波，得到预滤波信号M，该信号 中每个信号点M_i按式(1)计算

$M_i=\underset{j\in \Omega }{\Sigma }{\omega }_{(i,j)}·O_j---\left(1\right)$

其中，O_j代表原始信号O中位于点j处信号的幅值，Ω表示搜 索窗，ω_(i,j)是衡量原始信号O中分别以i和j为中心的两信号块O(N_i) 和O(N_j)之间的相似度，其计算方法如(2)式所示：

${\omega }_{(i,j)}=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O\left(N_i\right)-O\left(N_j\right)|{|}_2^2}{h^2})---\left(2\right)$

其中，O(N_i)和O(N_j)分别代表信号O中分别以i和j为中心的 两信号块，C_i为归一化参数，表示L₂范数，G表示高斯内核，* 表示卷积操作符号，h代表固定的衰减参数。

步骤2利用粒子群算法来对每一信号点i优化衰减参数和搜索 窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适 应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数 和搜索窗口。

(2-1)初始化粒子群。

以搜索窗口Ω_i和衰减参数h_i为优化对象，在二维目标搜索空间中 初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位 置向量x_k＝(S_k,h_k),(k＝1,2,…,N)，S_k为搜索窗Ω_i的半径。而第k个粒 子的移动速度对应的二维的向量为v_k＝(v_k1,v_k2),(k＝1,2,…,N)。

(2-2)对于信号O，用逐信号点的均方误差的期望E[MSE(M_i)]作 为优化衰减参数h_i和搜索窗Ω_i的适应度函数去计算每个粒子的适应 度Fit(k)。考虑到均方误差的期望E[MSE(M_i)]越小，即表示对应粒子 的搜索窗大小Ω_i和衰减参数h_i降噪效果越好。故适应度Fit(k)采用下 式(3)计算：

Fit(k)＝-E[MSE(M_i)](3)

其中均方误差的期望E[MSE(M_i)]如下式(4)计算：

E[MSE(M_i)]＝bias²(M_i)+var(M_i)(4)

其中偏差平方项bias²(M_i)和方差项var(M_i)计算公式分别如式 (5)和式(6)所示：

${\mathrm{bias}}^2\left(M_i\right)={\left(\frac{\underset{j\in {\Omega }_i}{\Sigma }(M_j-M_i)·{\omega }_{(i,j)}}{\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{(i,j)}}\right)}^2---\left(5\right)$

$\mathrm{var}\left(M_i\right)=\frac{{\sigma }^2\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{(i,j)}^2}{{\left(\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{(i,j)}\right)}^2}---\left(6\right)$

其中M_i是预滤波信号M中信号点i对应的幅度值,σ为信号O 噪声标准差。

(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置 p_best，并利用(3)式计算局部最优位置的适应度即个体最优fp_best(k)：

p_best＝(S_k,h_k),(k＝1,2,…,N)(7)

用对应每个粒子的适应度Fit(k)和个体极值fp_best(k)比较，如果 Fit(k)大于fp_best(k)，则用Fit(k)替换掉fp_best(k)；

(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置 g_best，并利用(3)式计算其适应度即全局极值fg_best：

g_best＝(S_g,h_g)(8)

对每个粒子，用它的适应度值Fit(k)和全局极值fg_best比较，如果 Fit(k)大于fg_best，则用Fit(k)替换掉fg_best；

(2-5)根据下式(9)和(10)更新粒子的速度和位置；

v_k＝T*v_k+c₁r₁(p_best-x_k)+c₂r₂(g_best-x_k)(9)

x_k＝x_k+v_k(10)

其中c₁和c₂是学习因子，r₁和r₂是介于[0,1]之间的随机数，T是 惯性权重。

如果满足结束条件(误差足够好或到达最大循环次数)退出，否则 返回(2-2)。

(2-6)粒子群优化完之后，会得到群体最优的位置g_best＝(S_g,h_g)， 于是就得到了逐信号点最优化的搜索窗S＝S_g，最优化的衰减参数 h＝h_g。

步骤3用步骤2粒子群算法所得逐信号最优的搜索窗半径S_g和 衰减参数h_g对原始信号O进行逐信号点的非局部均值降噪，得到降噪 之后的信号M′。该降噪信号中每一信号点i对应的幅度M_i′按式(11) 计算：

$M_i^{\prime }=\underset{j\in {\Omega }_i^{\prime }}{\Sigma }{\omega }_{(i,j)}^{\prime }·O_j---\left(11\right)$

其中Ω_i′表示信号中以信号点i为中心的搜索窗，其窗口大小为 步骤3得到的对应信号点i最优的S_g。ω′_(i,j)是衡量预处理信号M中 分别以i和j为中心的两信号块O′(N_i)和O′(N_j)之间的相似度，其计 算公式如下(12)：

${\omega }_{(i,j)}^{\prime }=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O^{\prime }\left(N_j\right)-O^{\prime }\left(N_i\right)|{|}_2^2}{{h_g}^2})---\left(12\right)$

O'(N_i)和O'(N_j)分别代表预处理信号M中以信号点i和j为中 心的信号块，h_g是对应信号点i最优的衰减参数。

优选地，所述心电信号降噪方法，步骤1非局部均值预滤波搜索 窗Ω窗口大小在800至1200之间，衰减参数h在15至20之间。为 了实现较好的噪声抑制效果，得到噪声少的预滤波信号，搜索窗Ω窗 口大小优选1000，衰减参数h优选为16。

优选地，所述心电信号降噪方法，步骤(2-1)所述的N取值在 15至30之间，为了兼顾粒子群算法优化效果和时间效率,优选N为 20。

优选地，所述心电信号降噪方法，其特征在于，考虑到心电信号 的细节尖峰区和信号平坦区的特征，步骤(2-1)所述搜索窗大小S_k的 搜索范围为500-2000，衰减参数h_k的搜索范围为5-30。

优选地，所述心电信号降噪方法，步骤(2-5)考虑到粒子群算法 优化迭代的效果和效率，所述惯性权重优选T＝0.7298，学习因子c₁和 c₂优选为c₁＝c₂＝1.4962，最大迭代次数优选为50。

优选地，为了提高步骤3自适应逐信号非局部均值降噪效果， 得到噪声较少的最终降噪效果，所述信号块O′(N_i)和O′(N_j)大小优选 10。

总体而言，按照本发明的上述技术构思与现有技术相比，主要具 备以下的技术优点：

1、在细节保留和噪声滤除的平衡性方面具有先天优势，同时考 虑到心电信号具有明显周期性和区域性的特点，对于心电信号的不同 波段区域，利用粒子群算法自适应的调整非局部均值降噪的核心参数， 即衰减系数和搜索窗大小；

2、能在有效抑制心电信号中噪声的同时更好地保护信号细节信 息。

附图说明

图1是本发明基于粒子群优化的心电信号自适应非局部均值降 噪方法的流程图；

图2是本发明粒子群算法步骤2的流程图；

图3(a)是模拟心电信号；

图3(b)是噪声方差为60时的原始信号；

图3(c)是对噪声方差为60的原始信号实施例1ANLM方法的降 噪结果；

图3(d)是对噪声方差为60的原始信号实施NLM方法的降噪结果；

图3(e)是对噪声方差为60的原始信号实施IIR滤波方法的降噪 结果；

图3(f)是对噪声方差为60的原始信号实施FIR滤波方法的降噪 结果；

图3(g)是对噪声方差为60的原始信号实施中值滤波方法的降噪 结果；

图3(h)是对噪声方差为60的原始信号实施小波变换滤波方法的 降噪结果；

图4(a)是真实心电信号；

图4(b)是实施例2ANLM方法的降噪结果；

图4(c)是实施NLM方法的降噪结果；

图4(d)是实施IIR滤波方法的降噪结果；

图4(e)是实施FIR滤波方法的降噪结果；

图4(f)是实施中值方法的降噪结果；

图4(g)是实施小波变换滤波方法的降噪结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合 附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描 述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外， 下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此 之间未构成冲突就可以相互组合。

本发明提供的心电信号降噪方法，简称ANLM，如图1所示，包 括以下步骤：

步骤1对原始心电信号进行非局部均值滤波，得到滤波结果。

对原始信号O进行非局部均值滤波，得到预滤波信号M，该信号 中每个信号点M_i按式(1)计算

$M_i=\underset{j\in \Omega }{\Sigma }{\omega }_{i,j}·O_j---\left(1\right)$

其中，O_j代表原始信号O中位于点j处信号的幅值，Ω表示搜 索窗，其大小选择为1000，ω_i,j是衡量信号O中分别以i和j为中心 的两信号块O(N_i)和O(N_j)之间的相似度，其计算方法如(2)式所示：

${\omega }_{i,j}=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O\left(N_i\right)-O\left(N_j\right)|{|}_2^2}{h^2})---\left(2\right)$

其中，C_i为归一化参数，表示L₂范数，G表示高斯内核，* 表示卷积操作符号，h代表固定的衰减参数，其优选为16。

步骤2利用粒子群算法来对每一信号点i优化衰减参数和搜索 窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适 应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数 和搜索窗口。

(2-1)初始化粒子群。

以搜索窗口Ω_i和衰减参数h_i为优化对象，在二维目标搜索空间中 初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位 置向量x_k＝(S_k,h_k),(k＝1,2,…,20)，S_k为搜索窗Ω_i的半径。而第k个 粒子的移动速度对应的二维的向量为v_k＝(v_k1,v_k2),(k＝1,2,…,20)。

优选的，S_k的搜索范围为[500,2000]，h_k的搜索范围为[5,30]。

(2-2)对于信号O，用逐信号点的均方误差的期望E[MSE(M_i)]作 为优化衰减参数h_i和搜索窗Ω_i的适应度函数去计算每个粒子的适应 度Fit(k)。考虑到均方误差的期望E[MSE(M_i)]越小，即表示对应粒子 的搜索窗大小Ω_i和衰减参数h_i降噪效果越好。故适应度Fit(k)采用下 式(3)计算：

Fit(k)＝-E[MSE(M_i)](3)

其中均方误差的期望E[MSE(M_i)]如下式(4)计算：

E[MSE(M_i)]＝bias²(M_i)+var(M_i)(4)

其中偏差平方项bias²(M_i)和方差项var(M_i)计算公式分别如式 (5)和式(6)所示：

${\mathrm{bias}}^2\left(M_i\right)={\left(\frac{\underset{j\in {\Omega }_i}{\Sigma }(M_j-M_i)·{\omega }_{i,j}}{\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}}\right)}^2---\left(5\right)$

$\mathrm{var}\left(M_i\right)=\frac{{\sigma }^2\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}^2}{{\left(\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}\right)}^2}---\left(6\right)$

其中M_i是预滤波信号M中信号点i对应的幅度值,σ为信号O 噪声标准差。

(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置 p_best，并利用(3)式计算局部最优位置的适应度即个体最优fp_best(k)：

p_best＝(S_k,h_k),(k＝1,2,…,20)(7)

用对应每个粒子的适应度Fit(k)和个体极值fp_best(k)比较，如果 Fit(k)大于fp_best(k)，则用Fit(k)替换掉fp_best(k)；

(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置 g_best，并利用(3)式计算其适应度即全局极值fg_best：

g_best＝(S_g,h_g)(8)

对每个粒子，用它的适应度值Fit(k)和全局极值fg_best比较，如 果Fit(k)大于fg_best，则用Fit(k)替换掉fg_best；

(2-5)根据下式(9)和(10)更新粒子的速度和位置；

v_k＝T*v_k+c₁r₁(p_best-x_k)+c₂r₂(g_best-x_k)(9)

x_k＝x_k+v_k(10)

其中c₁和c₂是学习因子，优选为1.4962，r₁和r₂是介于[0,1]之 间的随机数，T是惯性权重，优选为0.7298。

如果满足最大循环次数退出，否则返回(2-2)。其最大循环次数 设置为50。

(2-6)粒子群优化完之后，会得到群体最优的位置g_best＝(S_g,h_g)， 于是就得到了逐信号点最优化的搜索窗S＝S_g，最优化的衰减参数 h＝h_g。

步骤3用步骤2粒子群算法所得逐信号最优的搜索窗半径S_g和 衰减参数h_g对原始信号O进行逐信号点的非局部均值降噪，得到降噪 之后的信号M′。该降噪信号中每一信号点i对应的幅度M_i′按式(11) 计算：

$M_i^{\prime }=\underset{j\in {\Omega }_i^{\prime }}{\Sigma }{\omega }_{i,j}^{\prime }·O_j---\left(11\right)$

其中Ω_i′表示信号中以信号点i为中心的搜索窗，其窗口大小为 步骤3得到的对应信号点i最优的S_g。ω′_i,j是衡量预处理信号M中分 别以i和j为中心的两信号块O′(N_i)和O′(N_j)之间的相似度，其计算 公式如下(12)：

${\omega }_{i,j}^{\prime }=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O^{\prime }\left(N_j\right)-O^{\prime }\left(N_i\right)|{|}_2^2}{{h_g}^2})---\left(12\right)$

O′(N_i)和O′(N_j)分别代表预处理信号M中以信号点i和j为中 心的信号块，其信号块大小为10，h_g是对应信号点i最优的衰减参 数。

以下为实施例：

实施例1

对图3(a)所示的模拟心电信号加入方差为20、30、40、50、60 的高斯噪声，得到图3(b)原始信号O,接着采用本专利提出的心电信 号降噪方法进行降噪处理：

步骤1对原始心电信号图3(b)进行非局部均值滤波，得到滤波 结果。

对原始信号O进行非局部均值滤波，得到预滤波信号M，该信号 中每个信号点M_i按式(1)计算

$M_i=\underset{j\in \Omega }{\Sigma }{\omega }_{i,j}·O_j---\left(1\right)$

其中，O_j代表原始信号O中位于点j处信号的幅值，Ω表示搜 索窗，其大小选择为1000，ω_i,j是衡量信号O中分别以i和j为中心 的两信号块O(N_i)和O(N_j)之间的相似度，其计算方法如(2)式所示：

${\omega }_{i,j}=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O\left(N_i\right)-O\left(N_j\right)|{|}_2^2}{h^2})---\left(2\right)$

其中，C_i为归一化参数，表示L₂范数，G表示高斯内核，* 表示卷积操作符号，h代表固定的衰减参数，对于模拟信号优选为16。 C_i为归一化参数，以确保成立，采用下式计算

$C_i=\underset{j\in \Omega }{\Sigma }\exp (-\frac{G*\left|\right|O\left(N_i\right)-O\left(N_j\right)|{|}_2^2}{h^2})---\left(3\right)$

步骤2利用粒子群算法来对每一信号点i优化衰减参数和搜索 窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适 应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数 和搜索窗口。

(2-1)初始化粒子群。

以搜索窗口Ω_i和衰减参数h_i为优化对象，在二维目标搜索空间中 初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位 置向量x_k＝(S_k,h_k),(k＝1,2,…,20)，S_k为搜索窗Ω_i的半径。而第k个 粒子的移动速度对应的二维的向量为v_k＝(v_k1,v_k2),(k＝1,2,…,20)。

优选的，S_k的搜索范围为[500,2000]，h_k的搜索范围为[5,30]。

(2-2)对于信号O，用逐信号点的均方误差的期望E[MSE(M_i)]作 为优化衰减参数h_i和搜索窗Ω_i的适应度函数去计算每个粒子的适应 度Fit(k)。考虑到均方误差的期望E[MSE(M_i)]越小，即表示对应粒子 的搜索窗大小Ω_i和衰减参数h_i降噪效果越好。故适应度Fit(k)采用下 式(3)计算：

Fit(k)＝-E[MSE(M_i)](3)

其中均方误差的期望E[MSE(M_i)]如下式(4)计算：

E[MSE(M_i)]＝bias²(M_i)+var(M_i)(4)

其中偏差平方项bias²(M_i)和方差项var(M_i)计算公式分别如式 (5)和式(6)所示：

${\mathrm{bias}}^2\left(M_i\right)={\left(\frac{\underset{j\in {\Omega }_i}{\Sigma }(M_j-M_i)·{\omega }_{i,j}}{\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}}\right)}^2---\left(5\right)$

$\mathrm{var}\left(M_i\right)=\frac{{\sigma }^2\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}^2}{{\left(\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}\right)}^2}---\left(6\right)$

其中M_i是预滤波信号M中信号点i对应的幅度值,σ为信号O 噪声标准差。

(2-3)用下式表示每个粒子迄今为止搜索到的局部最优位置 p_best，并利用(3)式计算局部最优位置的适应度即个体最优fp_best(k)：

p_best＝(S_k,h_k),(k＝1,2,…,20)(7)

用对应每个粒子的适应度Fit(k)和个体极值fp_best(k)比较，如果 Fit(k)大于fp_best(k)，则用Fit(k)替换掉fp_best(k)；

(2-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置 g_best，并利用(3)式计算其适应度即全局极值fg_best：

g_best＝(S_g,h_g)(8)

对每个粒子，用它的适应度值Fit(k)和全局极值fg_best比较，如 果Fit(k)大于fg_best，则用Fit(k)替换掉fg_best；

(2-5)根据下式(9)和(10)更新粒子的速度和位置；

v_k＝T*v_k+c₁r₁(p_best-x_k)+c₂r₂(g_best-x_k)(9)

x_k＝x_k+v_k(10)

其中c₁和c₂是学习因子，优选为1.4962，r₁和r₂是介于[0,1]之 间的随机数，T是惯性权重，优选为0.7298。

如果满足最大循环次数退出，否则返回(2-2)。其最大循环次数 设置为50。

(2-6)粒子群优化完之后，会得到群体最优的位置g_best＝(S_g,h_g)， 于是就得到了逐信号点最优化的搜索窗S＝S_g，最优化的衰减参数 h＝h_g。

步骤3用步骤2粒子群算法所得逐信号最优的搜索窗半径S_g和 衰减参数h_g对原始信号O进行逐信号点的非局部均值降噪，得到降噪 之后的信号M′。该降噪信号中每一信号点i对应的幅度M_i′按式(11) 计算：

$M_i^{\prime }=\underset{j\in {\Omega }_i^{\prime }}{\Sigma }{\omega }_{ij}^{\prime }·O_j---\left(11\right)$

其中Ω_i′表示信号中以信号点i为中心的搜索窗，其窗口大小为 步骤3得到的对应信号点i最优的S_g。ω′_i,j是衡量预处理信号M中分 别以i和j为中心的两信号块O′(N_i)和O′(N_j)之间的相似度，其计算 公式如下(12)：

${\omega }_{i,j}^{\prime }=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O^{\prime }\left(N_j\right)-O^{\prime }\left(N_i\right)|{|}_2^2}{{h_g}^2})---\left(12\right)$

O′(N_i)和O′(N_j)分别代表预处理信号M中以信号点i和j为中 心的信号块，其信号块大小为10，h_g是对应信号点i最优的衰减参 数。

为了说明本专利提出的ANLM方法的优越性，我们将传统NLM 方法、IIR滤波、FIR滤波、中值滤波、小波变换滤波这五种算法作 为对比算法。

为客观评价各种方法的滤波性能，采用信噪比(SNR)和均方误 差(RMSE)作为评价指标，这两个指标其定义如下：

$SNR=10log\left(\frac{I^2}{{(I-In)}^2}\right)$

其中I为未被噪声污染的ECG信号，In为加高斯噪声的ECG信号

$\mathrm{RMSE}=\sqrt{\frac{\Sigma {(I-In)}^2}{n}}$

其中n为未加噪声ECG信号的长度。

表1给出了各种滤波算法对应仿真心电信号降噪结果的均方误 差和信噪比。从表1可看出，本专利提出的ANLM方法较其它比较 方法具有更低的RMSE和更高的SNR，这表明ANLM方法在所有的比 较方法中具有最优的恢复性能。

表1不同噪声水平下各种滤波算法对应仿真心电信号降噪结果 的信噪比和均方误差

从图3可以看出：IIR、FIR和中值滤波对应的心电降噪信号中 残留了较多的噪声，如图3(e)、图3(f)、图3(g)中包含平滑区的 长方形窗口所示；NLM和小波算法对心电图中的细节造成一定的破 坏，如图3(d)、图3(h)中包含尖峰的正方形窗口所示。相比之下， 我们提出的ANLM算法既能很好地保护边缘等信号细节，又能有效地 滤除平滑区的噪声，具有较其它比较算法更优的恢复性能。

实施例2

为了说明本算法的实用性，我们将其用于对儿科头部受伤病人在 500Hz采样率下采集到的实际心电图4(a)的降噪处理，采用本专利提 出的ANLM方法进行降噪处理的步骤如下：

步骤1对图4(a)原始心电信号进行噪声评估，得到噪声方差。

对原始信号O进行噪声评估，该信号的噪声方差σ′按下式计算

σ′＝1.4826*median(|R-median(R)|)

其中R＝{R₁,…,R_x}是通过手动选择的噪声信号O中同质区域的 局部残差集，并且$R_x=\frac{(2O_x-O_{x-1}-O_{x+1})}{\sqrt{6}}.$

步骤2对实际心电信号图4(a)进行非局部均值滤波，得到滤波 结果。

对图4(a)进行非局部均值滤波，得到预滤波信号M，该信号中 每个信号点M_i按式(1)计算

$M_i=\underset{j\in \Omega }{\Sigma }{\omega }_{i,j}·O_j---\left(1\right)$

其中，O_j代表原始信号O中位于点j处信号的幅值，Ω表示搜 索窗，其大小选择为1000，ω_i,j是衡量信号O中分别以i和j为中心 的两信号块O(N_i)和O(N_j)之间的相似度，其计算方法如(2)式所示：

${\omega }_{i,j}=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O\left(N_i\right)-O\left(N_j\right)|{|}_2^2}{h^2})---\left(2\right)$

其中，C_i为归一化参数，表示L₂范数，G表示高斯内核，* 表示卷积操作符号，h代表固定的衰减参数，对于实际信号其优选为 16。C_i为归一化参数，以确保成立，采用下式计算

$C_i=\underset{j\in \Omega }{\Sigma }\exp (-\frac{G*\left|\right|O\left(N_i\right)-O\left(N_j\right)|{|}_2^2}{h^2})$

步骤3利用粒子群算法来对每一信号点i优化衰减参数和搜索 窗口,将信号的偏差的平方和其方差的和作为适应度函数，使其自适 应于不同的信号点，进而得到每一个信号点对应的最优化的衰减参数 和搜索窗口。

(3-1)初始化粒子群。

以搜索窗口Ω_i和衰减参数h_i为优化对象，在二维目标搜索空间中 初始化N个粒子组成一个种群，其中第k个粒子表示为一个二维的位 置向量x_k＝(S_k,h_k),(k＝1,2,…,20)，S_k为搜索窗Ω_i的半径。而第k个 粒子的移动速度对应的二维的向量为v_k＝(v_k1,v_k2),(k＝1,2,…,20)。

优选的，S_k的搜索范围为[500,2000]，h_k的搜索范围为[5,30]。

(3-2)对于信号O，用逐信号点的均方误差的期望E[MSE(M_i)]作 为优化衰减参数h_i和搜索Ω_i的适应度函数去计算每个粒子的适应度 Fit(k)。考虑到均方误差的期望E[MSE(M_i)]越小，即表示对应粒子的 搜索窗大小Ω_i和衰减参数h_i降噪效果越好。故适应度Fit(k)采用下式 (3)计算：

Fit(k)＝-E[MSE(M_i)](3)

其中均方误差的期望E[MSE(M_i)]如下式(4)计算：

E[MSE(M_i)]＝bias²(M_i)+var(M_i)(4)

其中偏差平方项bias²(M_i)和方差项var(M_i)计算公式分别如式 (5)和式(6)所示：

${\mathrm{bias}}^2\left(M_i\right)={\left(\frac{\underset{j\in {\Omega }_i}{\Sigma }(M_j-M_i)·{\omega }_{i,j}}{\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}}\right)}^2---\left(5\right)$

$\mathrm{var}\left(M_i\right)=\frac{{\left({\sigma }^{\prime }\right)}^2\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}^2}{{\left(\underset{j\in {\Omega }_i}{\Sigma }{\omega }_{i,j}\right)}^2}---\left(6\right)$

其中M_i是预滤波信号M中信号点i对应的幅度值,σ′是步骤1 所得预估的信号噪声标准差。

(3-3)用下式表示每个粒子迄今为止搜索到的局部最优位置 p_best，并利用(3)式计算局部最优位置的适应度即个体最优fp_best(k)：

p_best＝(S_k,h_k),(k＝1,2,…,20)(7)

用对应每个粒子的适应度Fit(k)和个体极值fp_best(k)比较，如果 Fit(k)大于fp_best(k)，则用Fit(k)替换掉fp_best(k)；

(3-4)用下式表示整个粒子群迄今为止搜索到的全局最优位置 g_best，并利用(3)式计算其适应度即全局极值fg_best：

g_best＝(S_g,h_g)(8)

对每个粒子，用它的适应度值Fit(k)和全局极值fg_best比较，如 果Fit(k)大于fg_best，则用Fit(k)替换掉fg_best；

(3-5)根据下式(9)和(10)更新粒子的速度和位置；

v_k＝T*v_k+c₁r₁(p_best-x_k)+c₂r₂(g_best-x_k)(9)

x_k＝x_k+v_k(10)

其中c₁和c₂是学习因子，优选为1.4962，r₁和r₂是介于[0,1]之 间的随机数，T是惯性权重，优选为0.7298。

如果满足最大循环次数退出，否则返回(3-2)。其最大循环次数 设置为50。

(3-6)粒子群优化完之后，会得到群体最优的位置g_best＝(S_g,h_g)， 于是就得到了逐信号点最优化的搜索窗S＝S_g，最优化的衰减参数 h＝h_g。

步骤4用步骤3粒子群算法所得逐信号最优的搜索窗半径S_g和 衰减参数h_g对原始信号O进行逐信号点的非局部均值降噪，得到降噪 之后的信号M′。该降噪信号中每一信号点i对应的幅度M_i′按式(11) 计算：

$M_i^{\prime }=\underset{j\in {\Omega }_i^{\prime }}{\Sigma }{\omega }_{i,j}^{\prime }·O_j---\left(11\right)$

其中Ω_i′表示信号中以信号点i为中心的搜索窗，其窗口大小为 步骤3得到的对应信号点i最优的S_g。ω′_(i,j)是衡量预处理信号M中 分别以i和j为中心的两信号块O′(N_i)和O′(N_j)之间的相似度，其计 算公式如下(12)：

${\omega }_{i,j}^{\prime }=\frac{1}{C_i}\exp (-\frac{G*\left|\right|O^{\prime }\left(N_j\right)-O^{\prime }\left(N_i\right)|{|}_2^2}{{h_g}^2})---\left(12\right)$

O′(N_i)和O′(N_j)分别代表预处理信号M中以信号点i和j为中 心的信号块，其信号块大小为10，h_g是对应信号点i最优的衰减参 数。

为了说明本专利提出的ANLM方法的优越性，我们将传统NLM 方法、FIR滤波、IIR滤波、中值滤波、小波变换滤波这五种算法作 为对比算法。从图4对比可以明显看出：ANLM滤波后的实际心电图 更平滑，在边缘区保留了最多的细节信息。上述对比证实了本专利提 出的算法在实际心电图降噪中的有效性和优越性。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例 而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任 何修改、等同替换和改进等，均应包含在本发明的保护范围之内。