CN102609904A - 双变量非局部平均滤波x射线图像消噪方法 - Google Patents

Abstract

双变量非局部平均滤波X射线图像消噪方法，本发明特征是，方法为：1)模糊消噪窗的选择方法；2)双变量模糊自适应非局部平均滤波算法。本发明的有益效果为，为了更好的消除工业X射线扫描图像中存在的未知量子噪声的影响，提出了将不易处理的量子噪声模型转为常见的高斯加性噪声模型，运用模糊运算选择滤波器窗口的大小，寻找使误差函数最小的相关权值矩阵的双变量模糊自适应非线性平均滤波的X射线图像消噪方法。在本发明中，引入粒子群优化滤波参数，进而局部重建权值矩阵，降低了局部相关性对样本数据的影响，提高了算法收敛速度，提高了工业X射线扫描图像去噪处理的速度和精度，适用于对噪声模型不确定的X射线扫描图像的处理。

Description

双变量非局部平均滤波X射线图像消噪方法

技术领域

本发明属于X射线图像消噪领域，具体涉及一种模糊自适应参数调整的双变量非局部平均滤波(Fuzzy Adaptive Non Local means，简写为FANL means)算法的X射线图像的消噪方法。

背景技术

随着工业X射线探伤技术的不断发展，对X射线扫描图像的质量也提出了越来越多的要求，这就要求有效的消除实时检测过程中产生的噪声信息。由于X射线扫描图像灰度区间比较窄、缺陷边缘模糊、图像噪声多、缺陷特征有时被淹没的特点，影响了根据X射线图像对被检测工件进行分析和评价的效果。随着X射线采集设备的不断更新，新的X射线扫描机能更加全面、准确的描述工业图像的信息，这就导致图像中含有更多的像素，其像素间的冗余量大大增加。因此，通过模糊自适应优化算子去除信息中不相关的量，并保持其内部结构的不变性是非局部平均消噪亟待解决的关键技术。

非局部平均滤波算法利用冗余信息进行消噪处理，它是对传统局部消噪模型的一个革新，其主要思想不是用图像中单个像素的灰度值进行比较，而是对该像素周围的整个灰度的分布状况进行比较，根据整个灰度分布的相似性来估计权值，其本质是将原图像像素间的相关性映射到象空间进行处理。Buades指出应充分利用冗余信息为消噪服务，在Buades近年研究的理论分析和实验结果表明，传统非局部平均滤波消噪算法在主客观性能上都优于常见的图像消噪算法，如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等，它起源于邻域滤波算法，是对邻域滤波算法的一种推广，其权值根据像素周围整个区域灰度分布的相似性得到，在降低图像噪声的同时具有很强的保持图像空间分辨率的能力。利用传统的非局部平均滤波算法处理复杂图像时，其计算量较大，处理速度慢，尤其是在处理较大图像时，此问题更加突出；此外，该方法会在图像的平滑区域引入人工伪影，图像变得模糊，空间分辨率受到影响。

模糊自适应参数调整的双变量非局部平均滤波的X射线图像的消噪方法是利用模糊算法实现滤波参数选择窗口的自适应调整，以粒子群优化算子的概念和理论为基础，当有多个优化滤波参数需要选择时，利用模糊控制规则进行最佳滤波参数的选择，通过选择优化窗开启来完成进化搜索，获得更好处理效果，更快的收敛速度和全局寻优的能力。但到目前为止，还没有人将该算法应用于非局部平均优化方面。

针对上述问题，本发明对传统NL-means方法进行改进，算法从数据间的本征特性出发，寻找使重建误差最小的权重值，达到使误差函数最小的目的。同时，智能优化算法的本质决定了权重的确定不依赖于数据点间的距离，有效地降低图像邻域像素之间的相关性，运算复杂度降低，节省了计算时间，提高了算法运行速度，以适应工业生产任务中的检测需要。

发明内容

本发明对工业X射线扫描图像中产生的随机噪声模型的消噪情况进行了研究，采用基于模糊自适应参数调整的双变量非局部平均滤波快速消噪方法实现了工业X射线扫描图像的有效消噪处理。

模糊自适应参数调整有效地对X射线分块图像并行处理，降低了图像块之间的相关性，实现X射线扫描图像的快速消噪处理，同时，发明中引入x，y双变量保证了图像消噪过程的位置不变性，并且，由于其采用粒子群优化方法也可以获得较好的收敛性。因而，本发明利用粒子群优化算法从数据间的本征特性出发，寻找在X射线扫描图像中使重建误差最小的权重值，从而使得工业X射线消噪图像在保证精度的条件下获得较快的处理速度。

以下对本发明方法做进一步的说明，具体内容如下：

双变量非局部平均滤波X射线图像消噪方法，本发明特征是，方法为：

1).模糊消噪窗的选择方法

非局部平均滤波算法有一个前提假设：采样数据所在局部空间是线性的，即每个采样点与它的近邻点有相似关系，通过权重贡献值线性表示；

该算法的学习目标是：在低维空间中充分利用像素间的冗余关系，根据灰度分布的相似性来设置每个邻域中的权重，即假设镶入的映射窗在局部是线性的条件下，最小化不相关像素，重构原图像；

设c(x，y)为X射线扫描图像函数，r(x，y)为理想图像函数，n(x，y)为噪声图像函数，x，y为图像像素点的直角坐标系下的坐标，则有：

c(x，y)＝r(x，y)+n(x，y)                 (1)

现需要找一个模板来确定c(x，y)与r(x，y)各个元素之间的相关性，近而有效消除n(x，y)；

c(x，y)图像中像素点x，y对应的灰度值用c(i)来表示，则有

$\underset{i=1}{\overset{I}{\mathrm{\Σ}}}c\left(i\right)=\underset{y=1}{\overset{Y}{\mathrm{\Σ}}}\underset{x=1}{\overset{X}{\mathrm{\Σ}}}c(x,y)+\mathrm{\γ}---\left(2\right)$

其中，I＝X*Y为图像的大小尺寸，i为图像点x，y中的灰度像素，i∈(1，2....I)，x∈(1，2....X)，y∈(1，2....Y)c(i)为X射线扫描图像中第i点的灰度值，c(x，y)反映图像中点x，y处的灰度信息，γ为图像灰度修正因子，因此，通过调节每一点γ值，来调节X射线扫描图像的局部对比度，使去噪后的图像灰度分布适宜；

①随机噪声模型转换

X射线成像***图像降质的主要原因是***随机噪声，在X射线扫描图像上表现为水平或垂直的条纹，研究表明，X射线的产生以及与物质的相互作用产生的噪声，其分布认为是状态离散、时间连续的马尔可夫过程，在时间上和空间上都满足泊松随机过程，若独立增量过程Δc(t)，其增量的概率分布服从泊松分布，则有：

$p\left(k\right)\≈\frac{{\mathrm{\λ}}^k}{k!}{e}^{-\mathrm{\λ}}---\left(3\right)$

其中，k为随机变量增量Δc(t)出现的次数，k∈(1，2...K)，当k较大时，根据泊松分布公式计算误差十分复杂的实际应用中很不方便，因此，本文中采用对噪声图像进行斯蒂令(Stirling)近似公式，将泊松噪声转化为高斯白噪声；

斯蒂令(Stirling)近似公式认为，当k较大时，k！由公式(4)近似求得：

$k!\≈\sqrt{2\mathrm{\πk}}{k}^k{e}^{-k}---\left(4\right)$

此时，噪声从泊松分布转化为高斯分布，则X射线噪声分布表达式为：

$p\left(k\right)=\frac{1}{\sqrt{2\mathrm{\πh}}}{e}^{[-\frac{{(h-k)}^2}{2h}]}---\left(5\right)$

则原图像符合高斯白噪声模型，其数学期望E由式(6)所求：

$E{\left|\right|c(x_i,y_i)-c(x_j,y_j)\left|\right|}_{2a}^2=E{\left|\right|r(x_i,y_i)-r(x_j,y_j)\left|\right|}_{2a}^2+{2\mathrm{\σ}}^2---\left(6\right)$

其中，

为像素i，j所在邻域的的高斯加权欧式距离，a＞0为标准的高斯卷积核，噪声方差为σ，σ的大小决定模糊消噪窗口的尺寸；

②噪声水平估计

噪声水平与模块窗口大小和滤波器参数有着密切的关系，因此，要达到较好的去噪效果，就需要估计一幅图像中的噪声水平。对于高斯白噪声，均值为0，只需估计其标准差参数；

高斯噪声估计方法有很多种，常用的有基于图像块和基于滤波器的两种方法。基于图像块的方法将图像分为许多小块，计算其各自方差，以若干最小方差的均值作为估计结果；基于滤波器的方法先对图像进行一次平滑，再计算原图与平滑后图像的差别，由此差值图像估计噪声水平。本发明中先利用分块法将图像分成许多子块，利用并行滤波的方法对每个子块滤波，再利用粒子群优化算法找出子块中最大噪声水平，作为整体噪声水平估计结果；

Buades近年研究的理论分析和实验结果指出NL-Means算法在主客观性能上都优于常见的图像消噪算法，如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等，它的权值根据像素周围整个区域灰度分布的相似性得到，在降低图像噪声的同时具有很强的保持图像空间分辨率的能力；

③模糊噪声滤波模板窗

为了提高计算效率，选择两窗口并行运算，通过设置两个窗口尺寸，一个是像素邻域窗尺寸A×A，另一个是像素邻域窗搜索范围的窗口大小B×B，即在B×B大小的区域里面选择像素的邻域大小为A×A，执行自适应非局部平均(Adaptive Non Local means，简记作ANL-means)算法，B×B的窗口在A×A大小的区域里滑动，根据区域的相似性确定区域中心像素灰度的贡献权值。

2).双变量模糊自适应非局部平均滤波算法

FANL means算法中最核心的问题是利用高斯加性噪声水平求出使重建误差最小的相关局部重建权值矩阵，然而，该算法是针对局部图像块进行操作，研究者都采用与欧氏距离有关的变量来定义该权值矩阵，通过距离的远近来衡量影响的大小，这使得该算法对样本中的噪声很敏感，此外该算法收敛速度不够快。双变量模糊自适应非局部平均滤波算法建立在模糊消噪窗选择基础上，将每个小图像块进行非局部平均滤波去噪处理，使得每一个模块中都求出一对最优化滤波参数，并利用粒子群优化算法实现最优参数的更新操作，从而实现了目标的优化求解。所以，本发明利用粒子群优化算法从滤波器的最佳参数出发，寻找使X射线扫描图像特征向量降维重建误差最小的权重值，从而获得X射线扫描图像的有效去噪方法。

本发明在图像子块中进行ANL means滤波处理，根据样本点与其临近点的相关程序构造为粒子群优化算法的个体，然后粒子群的个体在寻优的过程中找到全局最优速度，确定全局最优位置，最终得到最相关近邻局部重建权值矩阵，从而有效地去掉不相关的冗余信息，实现X射线的有效消噪；

①双变量非局部平均滤波

设c(x，y)是一幅观测图像，即X射线检测图像，n是均值为0，方差为σ²的加性噪声，则把输入图像r(x，y)通过非局部平均(Non Local means，简记作NL means)算法映射到观测空间中，(x，y)定义一个二维有界区域，(x，y)∈R²，x_i，y_i与邻域点x_j，y_j之间的相关值NL(c(x_i，y_i))由(6)式求出：

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\∫e^{-\frac{(G_a*\left({|c(x_i,y_i)-c(x_j,y_j)|}^2\right)\left(0\right))}{h^2}}c(x_j,y_j)\mathrm{dxdy}---\left(6\right)$

其中变量为(x，y)，

为标准方差为a的高斯卷积核，h为滤波参数，主要由图像的噪声标准差决定，D(x_i，y_i)为NL means变换系数，它根据坐标(x，y)的相对位置的灰度值求出：

$D(x_i,y_i)={\∫e}^{\frac{(G_a*\left({\left|c(x_i,y_i)c(x_t,y_t)\right|}^2\right)\left(0\right))}{h^2}}{x}^{\′}\left(i\right)y^{\′}\left(i\right)\mathrm{di}---\left(7\right)$

对于数字化X射线扫描图像在离散域表示为c(i)，其像素点i∈(x，y)，C＝{c(i)，i∈I}，I为图像中像素集合，则像素点NL means算法根据最相关近邻局部重建权值矩阵转化为式(8)的形式求取像素点i与邻域之间的相关值NL(c(i))用相关近邻局部重建权值w(i，j)与其邻近点的灰度值c(j)来表示：

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\∈I}{\mathrm{\Σ}}w(i,j)c\left(j\right)---\left(8\right)$

其中，式(2)中得，点(x_i，y_i)处的灰度值为c(i)，由c(x，y)改变修正系数γ转化，其邻近点(x_j，y_j)的灰度值为c(j)，I为X射线检测图像的像素总点数，X为c(x，y)的宽度值，Y为c(x，y)的高度值；

其中，i＝{1，2...i...I}＝{(1，1)，(1，2)，...(xi，yi)...(X，Y)}，此时，相关近邻局部重建权值w(i，j)表示为：

$w(i,j)=\frac{1}{Z\left(i\right)}{e}^{{\left|\right|c(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(9\right)$

$Z\left(i\right)=e^{-{\left|\right|c(\mathrm{NL}(x_{i,},y_i)-c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(10\right)$

其中，像素点i与邻近点j之间相关近邻局部重建权值w(i，s)建立式(9)的关系，Z(i)为权重加权转化系数，

为点i与j所在的两个邻域的基于灰度级的高斯加权欧氏距离，权重值w(i，s)主要由参数a和h决定，a为标准的高斯卷积核的标准差，h为双变量NL means的滤波参数，a越大，权值就越小，表明像素x_i，y_i距离中心像素越远，a的大小由选定像素邻域的窗口大小A×A决定；随着滤波参数h的增大，人工伪影逐渐减弱，但空间分辨率也会下降，滤波参数h由噪声的标准差决定；

②粒子群优化

设在第p个子模块中得到的局部最优参数为a_p和h_p，为局部最优解向量，利用粒子群最优算法，每个粒子能够通过一定规则估计自身位置的适应值；每个粒子能够记住自己当前所找到的局部最好位置z_p即 $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right],$ 粒子的局部位置为z_pq，z_p∈z_pq＝{z_p1，z_p2，....，z_pq..z_pQ}，z_pq为局部群种中粒子的位置，Q为种群中包含的粒子数，此外还记住群体中所有个体找到的最好位置，称为全局最优z_g，即 $\left[\begin{array}{c}a\\ h\end{array}\right],$ 选出种群中最好的 $\left[\begin{array}{c}a\\ h\end{array}\right]$ 解向量；

粒子群优化算法指出所有粒子都有一个速度决定他们的方向和距离，称为局部最优速度，用v_p表示；粒子们追随当前的最优粒子在解空间中搜索，在每次迭代中，粒子根据以下式子更新位置和速度：

z_p+1＝z_p+v_p+1                    (11)

v_p+1＝w₁v_p+g₁(s_pq-z_p)+g₂(s_gq-z_p)        (12)

式中，Q是目标搜索空间的维数，P表示种群中个体的数量，计算s_p当前位置适应值，v_p为粒子p的飞行速度，s_p∈(s_p1，s_p2，s_pq，...s_pQ)为粒子迄今为止搜索到的最优位置，s_pq为种群中粒子经过的位置，s_p为速度为v_p时搜索到的最优位置；s_g＝(s_g1+s_g1+...s_gq+...s_gQ)为整个粒子群迄今为止搜索到的最优位置和，s_gq为粒子群当前记录的最佳位置，g₁和g₂为自学习因子，w₁为惯性权重；

③适用度评价函数

适用度评价函数是衡量个体优劣的标志，其作用是在分块子图模块中找出整体最佳参数a和h。本发明根据工业X射线扫描图像噪声模型的不确定性作为个体模型的特殊性，所采用的适用度评价函数为：

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2---\left(13\right)$

式中，u_pq表示粒子群中第p代的第q个个体，D表示种群中个体的数量，该算法优化的具体步骤如下：

①初始化种群：选择阈值ε和最大迭代次数u_max，初始迭代次数u＝0，粒子位置范围为z_min～z_max，粒子速度范围为v_max～v_min；

②对种群中的个体进行测量：测量每个粒子的适应值u_pq，选取此次群体最优适应值与历史群体最优适应值进行比较，得到当前为止的群体最优适应值u_p，且u_p∈{u_pq|u_p1，u_p2，...u_pq，...u_pQ}，并存取相应位置的值zg∈{z_pq|z_p1，z_p2，...z_pq，...z_pQ}；每个粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值u_g∈{u_gq|u_g1，u_g2，...u_gq，...u_gQ}，并存取相应位置z_g∈{z_gq|z_g1，z_g2，...z_gq，...z_gQ}；

③评价种群U_p，并保留最优解；

④粒子群操作：如果U_p＞ε，顺序执行，否则结束循环；

⑤根据规则，更新v_p，z_p；

⑥改变进化代数：进化代数加1，如仍未满足最大进化代数P_max算法转至本节步骤②继续进行。

本发明的有益效果为，为了更好的消除工业X射线扫描图像中存在的未知量子噪声的影响，提出了将不易处理的量子噪声模型转为常见的高斯加性噪声模型，运用模糊运算选择滤波器窗口的大小，寻找使误差函数最小的相关权值矩阵的双变量模糊自适应非线性平均滤波的X射线图像消噪方法。在本发明中，引入粒子群优化滤波参数，进而局部重建权值矩阵，降低了局部相关性对样本数据的影响，提高了算法收敛速度，提高了工业X射线扫描图像去噪处理的速度和精度，适用于对噪声模型不确定的X射线扫描图像的处理。

附图说明

图1为粒子群优化算法流程图；

图2为工业X射线扫描图像去噪原理图；

图3为不同噪声水平下模糊窗口类型的隶属函数示意图。

具体实施方式

双变量非局部平均滤波X射线图像消噪方法，本发明方法为，

1).模糊消噪窗的选择方法

非局部平均滤波算法有一个前提假设：采样数据所在局部空间是线性的，即每个采样点与它的近邻点有相似关系，通过权重贡献值线性表示；

该算法的学习目标是：在低维空间中充分利用像素间的冗余关系，根据灰度分布的相似性来设置每个邻域中的权重，即假设镶入的映射窗在局部是线性的条件下，最小化不相关像素，重构原图像；

设c(x，y)为X射线扫描图像函数，r(x，y)为理想图像函数，n(x，y)为噪声图像函数，x，y为图像像素点的直角坐标系下的坐标，则有：

c(x，y)＝r(x，y)+n(x，y)                    (1)

现需要找一个模板来确定c(x，y)与r(x，y)各个元素之间的相关性，近而有效消除n(x，y)；

c(x，y)图像中像素点x，y对应的灰度值用c(i)来表示，则有

$\underset{i=1}{\overset{I}{\mathrm{\Σ}}}c\left(i\right)=\underset{y=1}{\overset{Y}{\mathrm{\Σ}}}\underset{x=1}{\overset{X}{\mathrm{\Σ}}}c(x,y)+\mathrm{\γ}---\left(2\right)$

其中，I＝X*Y为图像的大小尺寸，i为图像点x，y中的灰度像素，i∈(1，2....I)，x∈(1，2....X)，y∈(1，2.....Y)c(i)为X射线扫描图像中第i点的灰度值，c(x，y)反映图像中点x，y处的灰度信息，γ为图像灰度修正因子，因此，通过调节每一点γ值，来调节X射线扫描图像的局部对比度，使去噪后的图像灰度分布适宜；

①随机噪声模型转换

X射线成像***图像降质的主要原因是***随机噪声，在X射线扫描图像上表现为水平或垂直的条纹，研究表明，X射线的产生以及与物质的相互作用产生的噪声，其分布认为是状态离散、时间连续的马尔可夫过程，在时间上和空间上都满足泊松随机过程，若独立增量过程Δc(t)，其增量的概率分布服从泊松分布，则有：

$p\left(k\right)\≈\frac{{\mathrm{\λ}}^k}{k!}{e}^{-\mathrm{\λ}}---\left(3\right)$

其中，k为随机变量增量Δc(t)出现的次数，k∈(1，2...K)，当k较大时，根据泊松分布公式计算误差十分复杂的实际应用中很不方便，因此，本文中采用对噪声图像进行斯蒂令(Stirling)近似公式，将泊松噪声转化为高斯白噪声；

斯蒂令(Stirling)近似公式认为，当k较大时，k！由公式(4)近似求得：

$k!\≈\sqrt{2\mathrm{\πk}}{k}^k{e}^{-k}---\left(4\right)$

此时，噪声从泊松分布转化为高斯分布，则X射线噪声分布表达式为：

$p\left(k\right)=\frac{1}{\sqrt{2\mathrm{\πh}}}{e}^{[-\frac{{(h-k)}^2}{2h}]}---\left(5\right)$

则原图像符合高斯白噪声模型，其数学期望E由式(6)所求：

$E{\left|\right|c(x_i,y_i)-c(x_j,y_j)\left|\right|}_{2a}^2=E{\left|\right|r(x_i,y_i)-r(x_j,y_j)\left|\right|}_{2a}^2+{2\mathrm{\σ}}^2---\left(6\right)$

其中，

为像素i，j所在邻域的的高斯加权欧式距离，a＞0为标准的高斯卷积核，噪声方差为σ，σ的大小决定模糊消噪窗口的尺寸；

②噪声水平估计

噪声水平与模块窗口大小和滤波器参数有着密切的关系，因此，要达到较好的去噪效果，就需要估计一幅图像中的噪声水平。对于高斯白噪声，均值为0，只需估计其标准差参数；

高斯噪声估计方法有很多种，常用的有基于图像块和基于滤波器的两种方法。基于图像块的方法将图像分为许多小块，计算其各自方差，以若干最小方差的均值作为估计结果；基于滤波器的方法先对图像进行一次平滑，再计算原图与平滑后图像的差别，由此差值图像估计噪声水平。本发明中先利用分块法将图像分成许多子块，利用并行滤波的方法对每个子块滤波，再利用粒子群优化算法找出子块中最大噪声水平，作为整体噪声水平估计结果；

Buades近年研究的理论分析和实验结果指出NL-Means算法在主客观性能上都优于常见的图像消噪算法，如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等，它的权值根据像素周围整个区域灰度分布的相似性得到，在降低图像噪声的同时具有很强的保持图像空间分辨率的能力；

③模糊噪声滤波模板窗

为了提高计算效率，选择两窗口并行运算，通过设置两个窗口尺寸，一个是像素邻域窗尺寸A×A，另一个是像素邻域窗搜索范围的窗口大小B×B，即在B×B大小的区域里面选择像素的邻域大小为A×A，执行自适应非局部平均(Adaptive Non Local means，简记作ANL-means)算法，B×B的窗口在A×A大小的区域里滑动，根据区域的相似性确定区域中心像素灰度的贡献权值。

2).双变量模糊自适应非局部平均滤波算法

FANL means算法中最核心的问题是利用高斯加性噪声水平求出使重建误差最小的相关局部重建权值矩阵，然而，该算法是针对局部图像块进行操作，研究者都采用与欧氏距离有关的变量来定义该权值矩阵，通过距离的远近来衡量影响的大小，这使得该算法对样本中的噪声很敏感，此外该算法收敛速度不够快。双变量模糊自适应非局部平均滤波算法建立在模糊消噪窗选择基础上，将每个小图像块进行非局部平均滤波去噪处理，使得每一个模块中都求出一对最优化滤波参数，并利用粒子群优化算法实现最优参数的更新操作，从而实现了目标的优化求解。所以，本发明利用粒子群优化算法从滤波器的最佳参数出发，寻找使X射线扫描图像特征向量降维重建误差最小的权重值，从而获得X射线扫描图像的有效去噪方法。

本发明在图像子块中进行ANL means滤波处理，根据样本点与其临近点的相关程序构造为粒子群优化算法的个体，然后粒子群的个体在寻优的过程中找到全局最优速度，确定全局最优位置，最终得到最相关近邻局部重建权值矩阵，从而有效地去掉不相关的冗余信息，实现X射线的有效消噪；

①双变量非局部平均滤波

设c(x，y)是一幅观测图像，即X射线检测图像，n是均值为0，方差为σ²的加性噪声，则把输入图像r(x，y)通过非局部平均(Non Local means，简记作NL means)算法映射到观测空间中，(x，y)定义一个二维有界区域，(x，y)∈R²，x_i，y_i与邻域点x_j，y_j之间的相关值NL(c(x_i，y_i))由(6)式求出：

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\∫e^{-\frac{(G_a*\left({|c(x_i,y_i)-c(x_j,y_j)|}^2\right)\left(0\right))}{h^2}}c(x_j,y_j)\mathrm{dxdy}---\left(6\right)$

其中变量为(x，y)，

为标准方差为a的高斯卷积核，h为滤波参数，主要由图像的噪声标准差决定，D(x_i，y_i)为NL means变换系数，它根据坐标(x，y)的相对位置的灰度值求出：

$D(x_i,y_i)={\∫e}^{\frac{(G_a*\left({\left|c(x_i,y_i)c(x_t,y_t)\right|}^2\right)\left(0\right))}{h^2}}{x}^{\′}\left(i\right)y^{\′}\left(i\right)\mathrm{di}---\left(7\right)$

对于数字化X射线扫描图像在离散域表示为c(i)，其像素点i∈(x，y)，C＝{c(i)，i∈I}，I为图像中像素集合，则像素点NL means算法根据最相关近邻局部重建权值矩阵转化为式(8)的形式求取像素点i与邻域之间的相关值NL(c(i))用相关近邻局部重建权值w(i，j)与其邻近点的灰度值c(j)来表示：

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\∈I}{\mathrm{\Σ}}w(i,j)c\left(j\right)---\left(8\right)$

其中，式(2)中得，点(x_i，y_i)处的灰度值为c(i)，由c(x，y)改变修正系数γ转化，其邻近点(x_j，y_j)的灰度值为c(j)，I为X射线检测图像的像素总点数，X为c(x，y)的宽度值，Y为c(x，y)的高度值；

其中，i＝{1，2...i...I}＝{(1，1)，(1，2)，...(x_i，y_i)....(X，Y)}，此时，相关近邻局部重建权值w(i，j)表示为：

$w(i,j)=\frac{1}{Z\left(i\right)}{e}^{{\left|\right|c(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(9\right)$

$Z\left(i\right)=e^{-{\left|\right|c(\mathrm{NL}(x_{i,},y_i)-c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(10\right)$

其中，像素点i与邻近点j之间相关近邻局部重建权值w(i，s)建立式(9)的关系，Z(i)为权重加权转化系数，为点i与j所在的两个邻域的基于灰度级的高斯加权欧氏距离，权重值w(i，s)主要由参数a和h决定，a为标准的高斯卷积核的标准差，h为双变量NL means的滤波参数，a越大，权值就越小，表明像素x_i，y_i距离中心像素越远，a的大小由选定像素邻域的窗口大小A×A决定；随着滤波参数h的增大，人工伪影逐渐减弱，但空间分辨率也会下降，滤波参数h由噪声的标准差决定；

②粒子群优化

设在第p个子模块中得到的局部最优参数为a_p和h_p，为局部最优解向量，利用粒子群最优算法，每个粒子能够通过一定规则估计自身位置的适应值；每个粒子能够记住自己当前所找到的局部最好位置z_p即 $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right],$ 粒子的局部位置为z_pq，z_p∈z_pq＝{z_p1，z_p2....，z_pq..z_pQ}，z_pq为局部群种中粒子的位置，Q为种群中包含的粒子数，此外还记住群体中所有个体找到的最好位置，称为全局最优z_g，即 $\left[\begin{array}{c}a\\ h\end{array}\right],$ 选出种群中最好的 $\left[\begin{array}{c}a\\ h\end{array}\right]$ 解向量；

粒子群优化算法指出所有粒子都有一个速度决定他们的方向和距离，称为局部最优速度，用v_p表示；粒子们追随当前的最优粒子在解空间中搜索，在每次迭代中，粒子根据以下式子更新位置和速度：

z_p+1＝z_p+v_p+1                (11)

v_p+1＝w₁v_p+g₁(s_pq-z_p)+g₂(s_gq-z_p)        (12)

式中，Q是目标搜索空间的维数，P表示种群中个体的数量，计算s_p当前位置适应值，v_p为粒子p的飞行速度，s_p∈(s_p1，s_p2，s_pq，...s_pQ)为粒子迄今为止搜索到的最优位置，s_pq为种群中粒子经过的位置，s_p为速度为v_p时搜索到的最优位置；s_g＝(s_g1+s_g1+...s_gq+...s_gQ)为整个粒子群迄今为止搜索到的最优位置和，s_gq为粒子群当前记录的最佳位置，g₁和g₂为自学习因子，w₁为惯性权重；

③适用度评价函数

适用度评价函数是衡量个体优劣的标志，其作用是在分块子图模块中找出整体最佳参数a和h。本发明根据工业X射线扫描图像噪声模型的不确定性作为个体模型的特殊性，所采用的适用度评价函数为：

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2---\left(13\right)$

式中，u_pq表示粒子群中第p代的第q个个体，D表示种群中个体的数量，该算法优化的具体步骤如下：

①初始化种群：选择阈值ε和最大迭代次数u_max，初始迭代次数u＝0，粒子位置范围为z_min～z_max，粒子速度范围为v_max～v_min；

②对种群中的个体进行测量：测量每个粒子的适应值u_pq，选取此次群体最优适应值与历史群体最优适应值进行比较，得到当前为止的群体最优适应值u_p，且u_p∈{u_pq|u_p1，u_p2，...u_pq，...u_pQ}，并存取相应位置的值z_g∈{z_pq|z_p1，z_p2，...z_pq，...z_pQ}；每个粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值u_g∈{u_gq|u_g1，u_g2，...u_gq，...u_gQ}，并存取相应位置z_g∈{z_gq|z_g1，z_g2，...z_gq，...z_gQ}；

③评价种群U_p，并保留最优解；

④粒子群操作：如果U_p＞ε，顺序执行，否则结束循环；

⑤根据规则，更新v_p，z_p；

⑥改变进化代数：进化代数加1，如仍未满足最大进化代数P_max算法转至本节步骤②继续进行。

本发明针对电场工业设备无损检测过程中X射线扫描图像灰度区间比较窄、缺陷边缘模糊、图像噪声多、缺陷特征有时被淹没等特点进行去噪研究，分别以水平量子噪声，垂直量子噪声，随机噪声，过增强处理产生的噪声为例进行研究，输入BMP、JPG、PNG等多种格式图像任意大小的扫描图像。首先由用户根据去噪指标要求提出自选参数(目标图像)，如用户不提供则由***自动判别，并将输入的图像的量子噪声转换为高斯噪声，然后，对高期噪声进行等级划分，依据噪声模型水平的分布利用模糊控制规则设置不同的窗口，并在每个子窗口中设置不同的滤波参数，利用并行粒子群优化算法进行全局最优滤波参数的选择，近而对整幅X射线扫描图像进行FANL去噪算法，并用最小均方误差(Minimum Mean SquaredError，即MMSE)、峰值信噪比(peak signal noise ratio，即PSNR)及去噪所用时间t评价***去噪性能。***结构框图如图2所示，具体实施方法如下：

1.噪声水平建模与模板窗尺寸选择

X射线成像***图像降质的主要原因是***随机噪声。研究表明，X射线的产生以及与物质的相互作用，在时间上和空间上都满足泊松随机过程。对于快速X射线成像***，由于曝光时间短，X射线所产生的量子噪声更为突出，严重影响了图像的质量。因此，对输入的X射线扫描图像，由于其噪声模型未知，因此，在原图像基础上加入均值μ₀为0，初始方差σ₀为0.02的低噪声水平在高斯白噪声，其中，μ₀，σ₀为加入高斯噪声的初始值，为一常值，此时，利用斯蒂令(Stirling)公式(4)对X射线扫描加高斯均匀噪声后的图像函数模型进行近似，将泊松噪声转化为高斯白噪声，再求出转换后图像的总体高斯白噪声的方差σ，σ取值影响***性能。

根据模糊噪声滤波模板窗选用原则，需要对检测后的X射线图像加高斯噪声处理后的图像进行模糊分块并行计算各局部最优参数，因此，需要设置两个窗口，像素邻域窗尺寸A×A，邻域窗搜索范围的窗口B×B，常见的噪声较大的图像，一般取A＝7，B＝23，对于低噪声的图像取A＝3，B＝7基本就能满足降噪。因此，噪声水平来设置非局部平均消噪的模板窗口大小，对于输入的噪声先加入参数为(u₀，σ₀)的高斯白噪声后转换噪声模型，根据转换后图像函数中的高期噪声水平σ的不同分为五个等级，即低噪声σ＜0.2，较低噪声0.2＜σ＜0.5，中等噪声0.5＜σ＜0.8，较高噪声0.8＜σ＜1和高噪声σ＞1五个论域区间，对于不同的噪声设置不同的处理模板窗，在每个模板窗里自适应地调整滤波参数，不同水平下噪声模糊窗口的类型水平如图3中所示。

2.模板窗中局部最优参数选择

(1)局部参数a

非局部平均滤波中像素点i与邻近像素点j的权重值w(i，j)主要由参数a和h决定，a为标准的高斯卷积核，h为双变量NL means的滤波参数，高斯卷积核a越大，权值就越小，表明像素x_i，y_i距离中心像素越远，a的大小由选定像素邻域的窗口大小A×A决定；a的取值与图像本身的平滑程度也有一定关系，当a值较大时，包含了较多的噪声点，但同时也会将图像中本应与当前像素取值不同的点包含进来，因此a值并非越大越好，a常取邻域窗口大小的1～10倍的范围内；同时兼顾不同噪声水平的影响，取噪声标准差σ的10～15倍作为参考值，考虑到所有与当前像素加权绝对差值不超过aσ的像素点。因此，当选择的最优参数不满足以上条件时，则舍弃已选择的局部最优参数值。

(2)局部参数h

研究得出，滤波参数h与噪声方差σ²有近似线性正比关系，并受到噪声图像方差的影响，为了不损失噪声之外的信息，当两个像素的加权距离大于阈值h_β时(β为滤波器设定的截止条件)，权值w(i，j)应接近于零。

当加噪处理后图像中像素点的分布符合高斯分布时，其噪声权值函数w(i，j)也符合均值为0，标准差为h/2的典型高斯函数，因此，通过改变h值来调节权值w(i，j)的分布，则w(i，j)满足高斯分布的特点，与当前像素点及其邻域像素灰度相同的点在叠加了的高斯噪声(σ)后会与当前像素值有一些差异，这些差异决定了i，j两点间基于灰度级的高斯加权欧氏距离值，有一定比例的像素落在与当前考虑像素值的加权绝对差值不超过噪声标准差β倍的范围之内。设定忽略权值和与计算加权距离时忽略噪声像素百分比，当权值距离大于标准偏差某一倍数的所有权值之和估计出来，设定忽略的权值和，求得到这一倍数值。现假设像素值的加权绝对差值不超过噪声标准差的倍数为β，则不被忽略的距离应满足：

${\left|\right|c\left({\mathrm{NL}}_i\right)-c\left({\mathrm{NL}}_j\right)\left|\right|}_{2a}^2\≤\frac{\mathrm{\βh}}{\sqrt{2}}---\left(14\right)$

则得滤波参数与噪声水平的近似关系为：

$h=\sqrt{2}{\mathrm{\β\σ}}^2---\left(15\right)$

X射线图像灰度级数为256，因此，采用下式估计h值：

$h=[\sqrt{2}\mathrm{\β}+\frac{\max \left(O({\mathrm{\σ}}_c-{\mathrm{\σ}}_0)\right)}{10}]{\mathrm{\σ}}^2---\left(16\right)$

σ_c为检测图像的平均标准差，σ₀为参考标准差阈值，当采用式(16)中的参数对噪声图像进行滤波时，达到近似最优的效果。h参数值过小会使许多噪声没有被考虑；h参数值过大，会平滑掉一些超出噪声范围的图像差异，从而造成图像过于平滑、丢失边缘等细节信息。对于灰度为256级的X射线扫描图像，小于σ₀的标准差对h值影响较小，而当噪声水平较高于σ₀时对h影响较大。

3.基于粒子群优化的滤波参数解向量选择

对X射线图像进行消噪处理时，首先要对整幅图像进行分块，再利用ANL means算法求出每个子块的局部最优参数解向量。由于X射线机扫描得到中的噪声模型未知，因此对不同X射线图像分块的个数有所不同，每个图像子块中的局部最优参数也不尽相同，因此，采用粒子群优化算法在X射线扫描图像中寻找全局最优参数解向量。

本发明在高斯噪声模型转换后的X射线扫描图像子块中根据各像素点灰度值的相关性来提取特征信息，其中，邻域窗尺寸A×A，搜索范围的窗口大小B×B，此时，不仅相邻点之间是相关的，不相邻点在设定范围内也存在着一定的相关性，权重值w(i，j)根据像素点i，j之间的相关性得到，在每个子块中表现出求出每个子块的局部最优参数a_i，h_i，再利用粒子群优化算法求出全局最佳参数a，h。

设输入X射线图像分为P个子块，考虑到像素间的相关性，每个子块之间也存在着一定的相关性，因此，每一个子块看作待选最优参数空间的一个样本点c_p，计算每一个样本点(特征向量)c_p的k个最近邻点，计算该点与其他P-1个样本点之间的距离，将距离排序，选择前k个与c_p最近的点作为其邻近点。

由每个样本点的近邻点计算出该样本点的局部重建权值矩阵 $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right].$ 用每个特征向量的近邻点对该特征向量进行重建，求取使重建误差最小的近邻局部重建权值矩阵。

定义重建拟合误差函数如下：

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2$

其中，c_q(q＝1，2，...，k)为c_p样本中的q个近邻点，u_pq是c_p与c_q之间的权值。

在满足以下两个约束条件时，通过最小化误差函数得到局部重建权值矩阵，即由样本点的近邻点，构造出最优W矩阵使误差函数值达到最小。

1)每一个数据点c_q都只能由它的相关点来表示，若c_q不是相关点，则u_pq＝0；

2)权值矩阵的每一行的和为1，即满足归一化约束

在求取最优U矩阵的过程中，U为重建拟合误差函数解向量构成的误差函数，本发明中以样本点与其相关点的重建权值向量集合作为粒子群算法的初始粒子，然后粒子群中的个体在寻优的过程中找到全局最优位置和最佳速度，最终得到近邻局部重建权值矩阵U。其中，u_pq表示种群中第p代的第q个个体，P表示种群中个体的数量。该算法优化的具体步骤如下所示，其算法流程如图1所示。

式中，u_pq表示粒子群中第p代的第q个个体，D表示种群中个体的数量，该算法优化的具体步骤如下所示，其算法流程如图1所示。

1)初始化种群：选择阈值ε和最大迭代次数u_max，初始迭代次数u＝0，粒子位置范围为z_min～z_max，粒子速度范围为v_max～v_min；

2)对种群中的个体进行测量：测量每个粒子的适应值u_pq，选取此次群体最优适应值与历史群体最优适应值进行比较，得到当前为止的群体最优适应值u_p，且u_p∈{u_pq|u_p1，u_p2，...u_pq，...u_pQ}，并存取相应位置的值z_g∈{z_pq|z_p1，z_p2，...z_pq，...z_pQ}；每个粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值u_g∈{u_gq|u_g1，u_g2，...u_gq，...u_gQ}，并存取相应位置z_g∈{z_gq|z_g1，z_g2，...z_gq，...z_gQ}。

3)评价种群U_p，并保留最优解；

4)粒子群操作：如果U_p＞ε，顺序执行，否则结束循环；

5)根据规则，更新v_p，z_p；

6)改变进化代数：进化代数加1，如仍未满足最大进化代数P_max算法转至步骤2)继续进行。

按照以上步骤即求取最优近邻局部重建权值矩阵U，为保持上步求取的权值U(P)不变，算法流程第2)步执行的约束条件为：

1)位置更新操作：z_p+1＝z_p+v_p+1；

2)速度更新操作：v_p+1＝w₁v_p+g₁(s_pq-z_p)+g₂(s_gq-z_p)

式中，q∈{q₁，q₂...Q}，Q是目标搜索空间的维数，P是样本中数据点的个数，计算s_p当前位置适应值，v_p为粒子p的飞行速度，s_p∈(s_p1，s_p2，s_pq，...s_pQ)为粒子迄今为止搜索到的最优位置，s_pq为种群中粒子经过的位置，s_p为速度为v_p时搜索到的最优位置；s_g＝(s_g1+s_g1+...s_gq+...s_gQ)为整个粒子群迄今为止搜索到的最优位置和，s_gq为粒子群当前记录的最佳位置，g₁和g₂为自学习因子，w₁为惯性权重。因此，最小化拟合误差函数中取f(u)的最小的Q个特征值对应的向量为列向量组成矩阵[a_p，h_p]^T，则U的列向量即为Q维空间的向量表示。

4.相似性度量

本发明中将输入的工业X射线扫描图像中的量子噪声转换为高斯噪声，近而对转换后的带有高斯噪声的图像利用FANL算法进行快速去噪。因此，对检测的带有不确定量子噪声的工业X射线检测图像和加入低噪声水平的噪声模型转换后的带有高斯噪声的X射线图像进行相似度测量，本发明采用如下相似测度方法：

(1)互信息度量

互信息作为衡量两幅图像配准的相似性测度函数，当互信息的值越大，表示两幅度图像越相似。

假设检测的带有不确定量子噪声的工业X射线检测图像为图像CR，转换后的带有高斯噪声的X射线图像为图像C，C按选择邻域窗分块后的子图集为C₁，C₂..C_a..C_A，C按搜索邻域窗分为C₁，C₂..C_b.C_B，CR，C为任意实数，则CR，C的边缘概率密度分别为P_CR(cr)和P_C(c)，它们的联合概率分布则表示为P_AB(a，b)，则用互相关信息表示两个信号之间的相关程度，则有

$I(\mathrm{CR},C)=\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log \frac{P_{\mathrm{CRC}}(\mathrm{cr},c)}{P_{\mathrm{CR}}\left(\mathrm{cr}\right)P_C\left(c\right)}---\left(17\right)$

互相关与熵有着密切的关系，式(17)写成式(18)的形式：

I(CR，C)＝H(CR，C)＝H(CR)+H(C)-H(CR，C)＝H(C)H(C/CR)    (18)

H(A)和H(B)分别为两幅图像的熵。一幅图像CR的信息熵，它与各像素点之前的关系表示为如下式的形式：

$H\left(\mathrm{CR}\right)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{pi}\log \mathrm{pi}---\left(19\right)$

i∈(1，2，....n)表示信号A的所含有的像素数，设H(CR，C)表示图像CR，C的联合熵，H(CR/C)，H(C/CR)分别表示图像的条件熵，则有：

$H\left(\mathrm{CR}\right|C)=-\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log P_{\mathrm{CRC}}\left(\mathrm{cr}\right|c)---\left(20\right)$

$H(\mathrm{CR},C)=-\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log P_{\mathrm{CRC}}(\mathrm{cr},c)---\left(21\right)$

设m(e)为原检测图像的灰度值在点s处的灰度值，e为图像中的像素点，Γ(T(e))为原图像燥声模型转换后在点e处进行Γ变换后的灰度值，灰度值大小设置为0～255之间，同理，Γ(e)与m(T(e))的灰度值范围是0～255，当C在CR上移动时，重叠部分形成二维直方图，记作两幅图组成的灰度对(Γ(e)，m(T(e)))的个数为h_α(Γm)，则两幅相关联的图像CR和C的联合分布概率为

边缘分布概率

和

于是互相关函数为：

$I(\∂)=\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC},\∂}(\mathrm{cr},c){\log }_2\frac{P_{\mathrm{CRC},\∂}(\mathrm{cr},c)}{P_{\mathrm{CR},\∂}\left(\mathrm{cr}\right)P_C\left(c\right)}---\left(22\right)$

最优的互信息参数

有

${\∂}^{*}=\mathrm{arc}\max I\left(\mathrm{\α}\right)---\left(23\right)$

分别求出邻域窗C₁，C₂..C_a..C_A与搜索窗口C₁，C₂..C_b.C_B的最优互信息参数

和

则最优值 ${{\∂}_o}^{*}=\max \{{\∂}_a^{*},{\∂}_b^{*}\}.$

(2)余弦距离度量

为了考察两幅度图像各像素位置之间的相关度，利用直方图的交求余弦距离法的方法找出图像位置上的相关关系。

假设F和Q是图像C和CR图像的包含w个点的直方图，则它们之间的相交距离l用下式表示：

$l=\underset{w=1}{\overset{N}{\mathrm{\Σ}}}\min (F_w,Q_w)---\left(24\right)$

直方图的相交是指两个直方图在每个点中共有的像素数量。有时，该值还通过除以其中一个直方图中所有的像素数量来实现标准化，从而使其值处于[0，1]的值域范围，则有：

$l(F,Q)=\underset{w=1}{\overset{N}{\mathrm{\Σ}}}\min (F_w,Q_w)/\underset{w=1}{\overset{N}{\mathrm{\Σ}}}{Q}_w---\left(25\right)$

求得余弦距离为：

l(F，Q)＝F^T*Q/(||F||*||Q||)          (26)

其中，F和Q分别表示查询图像和数据库中图像的特征向量，||*|表示向量范数。计算得到的相似性度量值在[0，1]之间，该值越大，表示图像越相似。

对于组合特征图像，相似性度量定义为各个特征相似性度量的加权和。其公式为：

$l(F,Q)=\underset{w=1}{\overset{m}{\mathrm{\Σ}}}{u}_w*l_w(F,Q)---\left(27\right)$

表示由m个特征组合而成，其中u_w表示第w个特征的权重系数，它表示第w个特征的重要性，一般取各u_w相等，S_w(F，Q)表示第w个特征的相似性度量函数值。

5.算法的性能评价

本发明采用去噪性能评价中常用的最小均方误差MMSE(Minimum Means SquaredError)，峰值信噪比PSNR和响应时间t来客观评价去噪性能的好坏。

MMSE评价去噪后图像数据的变化程度，MMSE的值越小，说明去噪效果越精确。原X射线检测图像的噪声模型是量子噪声，将其转换后得出的高斯噪声模型中的噪声均方误差为σ₀，图像经过本发明中的方法消噪处理后，可知消噪后新图像各点像素的信息熵，将消噪后图像的像素熵与原X射线扫描图像各点的像素信息熵比较，求出整体的均方差σ₁，则MMSE＝min{σ₀，σ₁}，MMSE反映中去噪后图像各点平均消噪水平的情况。为了说明图像中差值较大点对去噪效果的影响，本发明中引入PSNR作为评价图像去噪效果好坏的另一个客观标准，用于衡量图像中最大值点与噪声之间的变化情况，已知检测得到X射线扫描图像为CR，CR经噪声模型转换得到图像为C，C消噪后得到图像R^*，理想情况下期望得到的不含噪声在图像为R，则PSNR定义如下所示：

$\mathrm{PSNR}=10\×\log \left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(28\right)$

$\mathrm{MMSE}=\min \left\{\right|\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}({\left(R^x\right)}^n-R^n)}{N}\left|\right\}---\left(29\right)$

其中，N为图像尺寸，n为图像中的像素点，n∈N表示图像中所含的像素数。PSNR和MMSE表明该去噪处理***的效果的性能指标，实际中根据用户需求，设置接受的PSNR和MMSE阈值。

在衡量去噪***性能指标中响应时间也是一个忽略的指标，本发明中由于采用分块后各模块并行计算，因而在处理时间在优于传统的NL means去噪算法，但考虑到精准度要求方面，本发明中利用x(n)，y(n)来描述像素点n在图像中的位置信息，因而能更精确的处理指定区域的去噪情况。

应用传统非局部平均滤波消噪算法处理复杂图像时，其计算量较大，处理速度慢，尤其是在处理较大图像时，此问题更加突出；此外，传统非局部平均滤波方法会在图像的平滑区域引入人工伪影，图像变得模糊，空间分辨率受到影响。本发明中提出的模糊双变量非局部平均滤波消噪算法利用分块的方式把复杂图像进行分块，再对子块进行自适应参数并行运算，在降低图像复杂度的同时，节省了计算时间，达到了较好的去噪性能。

Claims (1)

1.双变量非局部平均滤波X射线图像消噪方法，其特征是，方法为：

1).模糊消噪窗的选择方法

非局部平均滤波算法有一个前提假设：采样数据所在局部空间是线性的，即每个采样点与它的近邻点有相似关系，通过权重贡献值线性表示；

此算法在低维空间中充分利用像素间的冗余关系，根据灰度分布的相似性来设置每个邻域中的权重，即假设镶入的映射窗在局部是线性的条件下，最小化不相关像素，重构原图像；

设c(x，y)为X射线扫描图像函数，r(x，y)为理想图像函数，n(x，y)为噪声图像函数，x，y为图像像素点的直角坐标系下的坐标，则有：

c(x，y)＝r(x，y)+n(x，y)            (1)

现需要找一个模板来确定c(x，y)与r(x，y)各个元素之间的相关性，近而有效消除n(x，y)；

c(x，y)图像中像素点x，y对应的灰度值用c(i)来表示，则有

$\underset{i=1}{\overset{I}{\mathrm{\Σ}}}c\left(i\right)=\underset{y=1}{\overset{Y}{\mathrm{\Σ}}}\underset{x=1}{\overset{X}{\mathrm{\Σ}}}c(x,y)+\mathrm{\γ}---\left(2\right)$

其中，I＝X*Y为图像的大小尺寸，i为图像点x，y中的灰度像素，i∈(1，2...I)，x∈(1，2...X)，y∈(1，2.....Y)c(i)为X射线扫描图像中第i点的灰度值，c(x，y)反映图像中点x，y处的灰度信息，γ为图像灰度修正因子，因此，通过调节每一点γ值，来调节X射线扫描图像的局部对比度，使去噪后的图像灰度分布适宜；

①随机噪声模型转换

X射线成像***图像降质的主要原因是***随机噪声，在X射线扫描图像上表现为水平或垂直的条纹，X射线的产生以及与物质的相互作用产生的噪声，其分布是状态离散、时间连续的马尔可夫过程，在时间上和空间上都满足泊松随机过程，若独立增量过程Δc(t)，其增量的概率分布服从泊松分布，则有：

$p\left(k\right)\≈\frac{{\mathrm{\λ}}^k}{k!}{e}^{-\mathrm{\λ}}---\left(3\right)$

其中，k为随机变量增量Δc(t)出现的次数，k∈(1，2...K)，当k较大时，根据泊松分布公式计算误差十分复杂的实际应用中很不方便，因此，采用对噪声图像进行斯蒂令(Stirling)近似公式，将泊松噪声转化为高斯白噪声；

斯蒂令(Stirling)近似公式认为，当k较大时，k！由公式(4)近似求得：

$k!\≈\sqrt{2\mathrm{\πk}}{k}^k{e}^{-k}---\left(4\right)$

此时，噪声从泊松分布转化为高斯分布，则X射线噪声分布表达式为：

$p\left(k\right)=\frac{1}{\sqrt{2\mathrm{\πh}}}{e}^{[-\frac{{(h-k)}^2}{2h}]}---\left(5\right)$

则原图像符合高斯白噪声模型，其数学期望E由式(6)所求：

$E{\left|\right|c(x_i,y_i)-c(x_j,y_j)\left|\right|}_{2a}^2=E{\left|\right|r(x_i,y_i)-r(x_j,y_j)\left|\right|}_{2a}^2+{2\mathrm{\σ}}^2---\left(6\right)$

其中，

为像素i，j所在邻域的的高斯加权欧式距离，a＞0为标准的高斯卷积核，噪声方差为σ，σ的大小决定模糊消噪窗口的尺寸；

②噪声水平估计

本发明方法中先利用分块法将图像分成许多子块，利用并行滤波的方法对每个子块滤波，再利用粒子群优化算法找出子块中最大噪声水平，作为整体噪声水平估计结果；

③模糊噪声滤波模板窗

设置两个窗口尺寸，一个是像素邻域窗尺寸A×A，另一个是像素邻域窗搜索范围的窗口大小B×B，即在B×B大小的区域里面选择像素的邻域大小为A×A，执行自适应非局部平均(Adaptive Non Localmeans，简记作ANL-means)算法，B×B的窗口在A×A大小的区域里滑动，根据区域的相似性确定区域中心像素灰度的贡献权值；

2).双变量模糊自适应非局部平均滤波算法

本发明方法利用粒子群优化算法从滤波器的最佳参数出发，寻找使X射线扫描图像特征向量降维重建误差最小的权重值，从而获得X射线扫描图像的有效去噪方法；

本发明在图像子块中进行ANL means滤波处理，根据样本点与其临近点的相关程序构造为粒子群优化算法的个体，然后粒子群的个体在寻优的过程中找到全局最优速度，确定全局最优位置，最终得到最相关近邻局部重建权值矩阵，从而有效地去掉不相关的冗余信息，实现X射线的有效消噪；

①双变量非局部平均滤波

设c(x，y)是一幅观测图像，即X射线检测图像，n是均值为0，方差为σ²的加性噪声，则把输入图像r(x，y)通过非局部平均(Non Local means，简记作NL means)算法映射到观测空间中，(x，y)定义一个二维有界区域，(x，y)∈R²，x_i，y_i与邻域点x_j，y_j之间的相关值NL(c(x_i，y_i))由(6)式求出：

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\∫e^{-\frac{(G_a*\left({|c(x_i,y_i)-c(x_j,y_j)|}^2\right)\left(0\right))}{h^2}}c(x_j,y_j)\mathrm{dxdy}---\left(6\right)$

其中变量为(x，y)，

为标准方差为a的高斯卷积核，h为滤波参数，主要由图像的噪声标准差决定，D(x_i，y_i)为NL means变换系数，它根据坐标(x，y)的相对位置的灰度值求出：

$D(x_i,y_i)={\∫e}^{\frac{(G_a*\left({\left|c(x_i,y_i)c(x_t,y_t)\right|}^2\right)\left(0\right))}{h^2}}{x}^{\′}\left(i\right)y^{\′}\left(i\right)\mathrm{di}---\left(7\right)$

对于数字化X射线扫描图像在离散域表示为c(i)，其像素点i∈(x，y)，C＝{c(i)，i∈I}，I为图像中像素集合，则像素点NL means算法根据最相关近邻局部重建权值矩阵转化为式(8)的形式求取像素点i与邻域之间的相关值NL(c(i))用相关近邻局部重建权值w(i，j)与其邻近点的灰度值c(j)来表示：

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\∈I}{\mathrm{\Σ}}w(i,j)c\left(j\right)---\left(8\right)$

其中，式(2)中得，点(x_i，y_i)处的灰度值为c(i)，由c(x，y)改变修正系数γ转化，其邻近点(x_j，y_j)的灰度值为c(j)，I为X射线检测图像的像素总点数，X为c(x，y)的宽度值，Y为c(x，y)的高度值；

其中，i＝{1，2...i...I}＝{(1，1)，(1，2)，...(x_i，y_i)....(X，Y)}，此时，相关近邻局部重建权值w(i，j)表示为：

$w(i,j)=\frac{1}{Z\left(i\right)}{e}^{{\left|\right|c(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(9\right)$

$Z\left(i\right)=e^{-{\left|\right|c(\mathrm{NL}(x_{i,},y_i)-c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(10\right)$

其中，像素点i与邻近点j之间相关近邻局部重建权值w(i，s)建立式(9)的关系，Z(i)为权重加权转化系数，

为点i与j所在的两个邻域的基于灰度级的高斯加权欧氏距离，权重值w(i，s)主要由参数a和h决定，a为标准的高斯卷积核的标准差，h为双变量NL means的滤波参数，a越大，权值就越小，表明像素x_i，y_i距离中心像素越远，a的大小由选定像素邻域的窗口大小A×A决定；随着滤波参数h的增大，人工伪影逐渐减弱，但空间分辨率也会下降，滤波参数h由噪声的标准差决定；②粒子群优化

设在第p个子模块中得到的局部最优参数为a_p和h_p，为局部最优解向量，利用粒子群最优算法，每个粒子能够通过一定规则估计自身位置的适应值；每个粒子能够记住自己当前所找到的局部最好位置z_p即 $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right],$ 粒子的局部位置为z_pq，z_p∈z_pq＝{z_p1，z_p2....，z_pq..z_pQ}，z_pq为局部群种中粒子的位置，Q为种群中包含的粒子数，此外还记住群体中所有个体找到的最好位置，称为全局最优z_g，即 $\left[\begin{array}{c}a\\ h\end{array}\right],$ 选出种群中最好的 $\left[\begin{array}{c}a\\ h\end{array}\right]$ 解向量；

粒子群优化算法指出所有粒子都有一个速度决定他们的方向和距离，称为局部最优速度，用v_p表示；粒子们追随当前的最优粒子在解空间中搜索，在每次迭代中，粒子根据以下式子更新位置和速度：

z_p+1＝z_p+v_p+1            (11)

v_p+1＝w₁v_p+g₁(s_pq-z_p)+g₂(s_gq-z_p)        (12)

式中，Q是目标搜索空间的维数，P表示种群中个体的数量，计算s_p当前位置适应值，v_p为粒子p的飞行速度，s_p∈(s_p1，s_p2，s_pq，...s_pQ)为粒子迄今为止搜索到的最优位置，s_pq为种群中粒子经过的位置，s_p为速度为v_p时搜索到的最优位置；s_g＝(s_g1+s_g1+...s_gq+...s_gQ)为整个粒子群迄今为止搜索到的最优位置和，s_gq为粒子群当前记录的最佳位置，g₁和g₂为自学习因子，w₁为惯性权重；

③适用度评价函数

适用度评价函数是衡量个体优劣的标志，其作用是在分块子图模块中找出整体最佳参数a和h；根据工业X射线扫描图像噪声模型的不确定性作为个体模型的特殊性，所采用的适用度评价函数为：

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2---\left(13\right)$

式中，u_pq表示粒子群中第p代的第q个个体，D表示种群中个体的数量，该算法优化的具体步骤如下：

①初始化种群：选择阈值ε和最大迭代次数u_max，初始迭代次数u＝0，粒子位置范围为z_min～z_max，粒子速度范围为v_max～v_min；

②对种群中的个体进行测量：测量每个粒子的适应值u_pq，选取此次群体最优适应值与历史群体最优适应值进行比较，得到当前为止的群体最优适应值u_p，且u_p∈{u_pq|u_p1，u_p2，...u_pq，...u_pQ}，并存取相应位置的值z_g∈{z_pq|z_p1，z_p2，...z_pq，...z_pQ}；每个粒子通过比较历史上自身适应值大小得到自己当前为止的最优适应值u_g∈{u_gq|u_g1，u_g2，...u_gq，...u_gQ}，并存取相应位置z_g∈{z_gq|z_g1，z_g2，...z_gq，...z_gQ}；

③评价种群U_p，并保留最优解；

④粒子群操作：如果U_p＞ε，顺序执行，否则结束循环；

⑤根据规则，更新v_p，z_p；

⑥改变进化代数：进化代数加1，如仍未满足最大进化代数P_max算法转至本节步骤②继续进行。

Images