-
摘要: 在压缩感知理论中，设计好的稀疏重构算法是一个比较重要，同时也是一个具有挑战性的问题.稀疏重构的基本目标是用较少的数据样本，通过解一个优化问题完成信号或者图像重构.关于稀疏重构过程，一个重要的研究方向是在数据受噪声干扰的情况下，如何高效快速地重建原信号.本文提出了基于共轭梯度最小二乘法（Conjugate gradient least squares，CGLS）和最小二乘QR分解（Least squares QR，LSQR）的联合优化的匹配追踪算法.该算法采用Alpha散度来测量CGLS和LSQR之间的离散度（差异度），并通过离散度来选择最优的解序列.实验分析表明基于CGLS和LSQR的联合优化的匹配追踪算法在压缩采样的信号受噪声干扰情况下具有较好的恢复能力.Abstract: For compressed sensing theory, to design a good sparse reconstruction algorithm is a challenge. The basic purpose of sparse reconstruction is to implement the signal or image reconstruction by solving optimization problem on the condition of fewer data samples. For the sparse reconstruction, an important aspect is how to reconstruct original signal when data is contaminated by noise. In this article, we present a matching pursuit algorithm, in which the conjugate gradient least squares (CGLS) method combines the least squares QR (LSQR) method to solve the optimization. This algorithm uses alpha divergence to measure dispersion between CGLS and LSQR, then selects optimization solution sequence in light of the dispersion. Experiment shows that the matching pursuit algorithm has excellent reconstruction performance when the signal of compressed sampling is contaminated.
-
Key words:
- Compressed sensing /
- matching pursuit /
- sparse reconstruction /
- noise
-
单晶炉是利用直拉法生产集成电路材料的主要设备, 在其产业链中处于首要地位.在单晶炉直拉硅单晶的工作过程中, 为保证其品质, 要求硅单晶的直径必须保持恒定[1-3].直径控制系统依据直径测量值, 通过调整热场温度和拉速对直径进行控制[4].但由于单晶炉内高温、负压的环境, 无法对硅单晶直径进行直接测量.目前国内外大部分学者和单晶炉厂商主要利用CCD摄像机拍摄硅单晶周围形成的高亮光环, 然后对图像进行分割和边缘提取, 最后通过测量光环内圆的直径, 从而获得硅单晶直径的测量值[5].由于高亮光环和晶体分界线较为模糊, 图像直方图峰值较多, 使得准确分割光环图像较为困难.在传统的硅单晶直径测量中, 一般由操作人员根据经验和观察通过设置单阈值进行分割.随着低纳米线距集成电路对硅材料要求的不断提高和大尺寸晶体生长过程中图像对比度、亮度等发生变化较明显, 传统方法不能达到测量要求.因此, 如何准确分割高亮光环图像, 以满足硅单晶直径测量的需求, 成为硅单晶直径检测中的一个关键问题.
目前, 图像分割方法主要有阈值分割法[6]、聚类分割法[7]、基于区域的分割方法[8]和基于形态学分水岭的分割方法[9]等.阈值分割法以其简单、易实现的特性被广泛应用于工业和医学等领域[10-11]．为提高阈值分割法的鲁棒性, 国内外学者提出了很多改进方法, 文献[12]提出了一种基于二维直方图的Otsu分割方法, 同时考虑了图像灰度值分布和它们邻域内平均灰度值分布, 提高了图像分割的抗噪性; 文献[13]和文献[14]分别提出了一种基于Otsu和最大熵的二维直方图区域斜分方法, 降低了区域错分率; 文献[15]将粒子群优化算法应用到二维直方图区域分割方法中, 提高了分割的速度和准确性; 文献[16]将二维直方图区域斜分方法扩展为多阈值分割, 简化了距离测度函数, 并且应用到回转窑火焰图像分割中.以上方法虽对二维直方图区域分割方法提出了改进, 但分割的测度函数仅基于最大类间方差或最大熵的一种, 分割精度难以满足硅单晶直径检测的需求.近年来随着多目标群智能优化方法的提出[17-21], 使得同时使用最大类间方差和最大熵作为测度函数以进一步提高图像分割的准确性成为可能.
本文提出了一种基于多目标人工鱼群算法的二维直方图区域斜分多阈值分割方法, 以提高对硅单晶直径测量图像高亮光环的分割精度, 从而提高硅单晶直径检测的精度.首先设计了一种多目标人工鱼群算法, 改进了快速构造Pareto非劣解集的方法, 然后以最大类间方差和最大熵同时作为测度函数, 搜索最优的二维直方图区域斜分分割阈值.最后利用基准函数对多目标人工鱼群算法进行有效性验证, 并且利用硅单晶直径检测图像对本文所提出的二维直方图区域斜分多阈值分割方法进行验证.实验结果表明, 本文方法能够提高高亮光环的分割精度, 对提高硅单晶直径检测精度具有重要意义.
1. 问题描述
1.1 硅单晶直径检测
在单晶炉工作的过程中, 由于炉室内处于高温、负压的环境, 不能直接对硅单晶直径进行测量, 需要在炉室侧上方用CCD摄像机对硅单晶周围形成的高亮光环进行拍摄, 然后对图像进行分割, 测量光环内圆的直径, 从而得到硅单晶的直径, 测量原理如图 1所示.因此, 对硅单晶图像的高亮光环进行分割是测量硅单晶直径的前提, 分割算法的准确性直接影响直径测量的精度. CCD拍摄的图像如图 2所示, 图像中主要包含盖板、导流筒、硅单晶、硅熔液和高亮光环五个主要区域.由于硅单晶图像直方图存在多个峰值, 如图 3所示, 且高亮光环、硅单晶底部与硅熔液灰度值接近且较高, 分界线较模糊, 一般的单阈值分割方法(如Otsu法和最大熵分割方法)不能单独将高亮光环分割出来, 而是将灰度值较低的盖板和硅单晶顶部区域分割为背景区域, 将灰度值较高的导流筒、高亮光环、硅单晶底部与硅熔液分割为目标区域.随着二维直方图区域分割和多阈值分割方法的提出, 使得图像分割结果的抗噪性和精度均有了很大提高[12-16]．
1.2 二维直方图区域斜分方法
二维直方图区域直分法最早由文献[12]提出, 以图像像素灰度值与该像素邻域像素灰度平均值为二元组构建出图像二维直方图, 再利用阈值将直方图分为A、B、C、D四个区域, 如图 4所示.该算法明显的不足是忽略了C和D区域, 区域错分率较高[13-14].文献[13]和文献[14]分别提出了二维直方图区域斜分方法, 利用垂直于主对角线的直线g=-f + T, T ∈ [0, 2L -2](式中, T为最佳阈值, L为灰度级数)将图 4中的直方图分为E和F两个区域, 如图 5所示.该方法在选取阈值时考虑了所有区域, 降低了区域错分率, 分割后的二维图像如下式
$ b(x, y) = \left\{ {\begin{array}{*{20}{l}} {0, }&{f(x, y) + g(x, y) \le T}\\ {1, }&{f(x, y) + g(x, y) > T} \end{array}} \right. $
(1) 文献[13]为简化计算测度函数, 依据下式对横纵坐标进行了转换
$ \left\{ \begin{array}{l} x'=\frac{\sqrt{2}}{2}(x+y) \\ y'=-\frac{\sqrt{2}}{2}(x-y)+\frac{\sqrt{2}}{2}(L-1) \end{array} \right. $
(2) 设p(x', y')为(x', y')发生的概率, 选取下式为距离测度函数, 即二维直方图的最大类间方差[13]
$ \begin{align} & {\sigma}_{b}=\omega_{0}[(\boldsymbol{\mu}_{0}-\boldsymbol{\mu}_{z})(\boldsymbol{\mu}_{0}-\boldsymbol{\mu}_{z})^{\rm {T}}]+\notag\\ &\qquad\omega_{1}[(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{z})(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{z})^{\rm {T}}] \end{align} $
(3) 式中, 两个区域出现的概率为
$ \begin{array}{l} {\omega _0} = \sum\limits_{x' = 0}^t {\sum\limits_{y' = 0}^{\sqrt 2 (L - 1)} p } (x', y')\\ {\omega _1} = \sum\limits_{x' = t + 1}^{\sqrt 2 (L - 1)} {\sum\limits_{y' = 0}^{\sqrt 2 (L - 1)} p } (x', y') \end{array} $
两个区域对应的均值向量为
$ \begin{array}{l} {\boldsymbol{\mu} _0} = {({\mu _{0x'}}, {\mu _{0y'}})^{\rm{T}}} = \\ \;\;\;\;\;(\frac{1}{{{\omega _0}}}\sum\limits_{x' = 0}^t {\sum\limits_{y' = 0}^{\sqrt 2 (L - 1)} {x'} } p(x', y'), \\ \;\;\;\;\;\;\frac{1}{{{\omega _0}}}\sum\limits_{x' = 0}^t {\sum\limits_{y' = 0}^{\sqrt 2 (L - 1)} {y'} } p(x', y'){)^{\rm{T}}}\\ {\boldsymbol{\mu} _1} = {({\mu _{1x'}}, {\mu _{1y'}})^{\rm{T}}} = \\ \;\;\;\;\;\;\;(\frac{1}{{{\omega _1}}}\sum\limits_{x' = t + 1}^{\sqrt 2 (L - 1)} {\sum\limits_{y' = 0}^{\sqrt 2 (L - 1)} {x'} } p(x', y'), \\ \;\;\;\;\;\;\;\;\frac{1}{{{\omega _1}}}\sum\limits_{x' = t + 1}^{\sqrt 2 (L - 1)} {\sum\limits_{y' = 0}^{\sqrt 2 (L - 1)} {y'} } p(x', y'){)^{\rm{T}}} \end{array} $
二维直方图总的均值向量为
$ \begin{align*} \begin{aligned} \boldsymbol{\mu}_{z}=&(\mu_{zx'}, \mu_{zy'})^{\rm {T}}=\notag\\ &\Big(\sum^{\sqrt{2}(L-1)}_{x'=0}\sum^{\sqrt{2}(L-1)}_{y'=0}x'p(x', y'), \\ &\sum^{\sqrt{2}(L-1)}_{x'=0}\sum^{\sqrt{2}(L-1)}_{y'=0}y'p(x', y')\Big)^{\rm {T}} \end{aligned} \end{align*} $
由于σb为一个矩阵, 一般都选用该矩阵的迹为测度函数[12-13], 如下式
$ \begin{align} \begin{aligned} {\rm tr}{\sigma}_{b}(t)=&\omega_{0}[(\mu_{0x'}-\mu_{zx'})^{2}+(\mu_{0y'}-\mu_{zy'})^{2}]+\\ &\omega_{1}[(\mu_{1x'}-\mu_{zx'})^{2}+(\mu_{1y'}-\mu_{zy'})^{2}] \end{aligned} \end{align} $
(4) 文献[14]采用变换前的原坐标系, 设p(x, y)为(x, y)发生的概率, 选取最大Shannon熵作为阈值选取测度函数, 如下式所示
$ \begin{align} t^{*}=\arg\max_{0<t\leq2L-2}\{\min_{0<t\leq2L-2}\{H_{0}(t), H_{B}(t)\}\} \end{align} $
(5) 式中
$ \begin{array}{l} {H_0}(t) = \log [{P_0}(t)] + \frac{{{h_0}(t)}}{{{P_0}(t)}}\\ {H_B}(t) = \log [1-{P_0}(t)] + \frac{{{h_{LL}} - {h_0}(t)}}{{1 - {P_0}(t)}}\\ {h_{LL}} = - \sum\limits_{x = 0}^{L - 1} {\sum\limits_{y = 0}^{L - 1} p } (x, y)\log [p(x, y)] \end{array} $
当0 < t≤L -1时,
$ \begin{array}{l} {P_0}(t) = \sum\limits_{x = 0}^t {\sum\limits_{y = 0}^{t - x} p } (x, y)\\ {h_0}(t) = - \sum\limits_{x = 0}^t {\sum\limits_{y = 0}^{t - x} p } (x, y)\log [p(x, y)] \end{array} $
当L -1 < t≤2L -2时,
$ \begin{array}{l} {P_0}(t) = 1 - \sum\limits_{x = t - L + 1}^{L - 1} {\sum\limits_{y = t - x}^{L - 1} p } (x, y)\\ {h_0}(t) = {h_{LL}} + \sum\limits_{x = t - L + 1}^{L - 1} {\sum\limits_{y = t - x}^{L - 1} p } (x, y)\log [p(x, y)] \end{array} $
上述两种方法的测度函数仅适用于单阈值分割, 不能满足分割硅单晶直径检测图像的需求, 需将其拓展为多阈值分割, 然后选取分割结果中最亮的区域为高亮光环.随着分割方法拓展到多阈值, 确定多个分割阈值的计算复杂度呈指数增长, 穷举法搜索最佳阈值的运算量变得非常巨大, 所以, 利用智能优化算法搜索最佳阈值成为必然的选择[16].近几年随着多目标智能优化算法的提出[17-21], 使得同时使用两种测度函数作为目标以提高图像阈值分割的精度成为可能.为此, 本文提出一种基于多目标人工鱼群算法(Multi-objective artificial fish swarm algorithm, MAFSA)的二维直方图区域斜分多阈值图像分割方法, 同时使用两种测度函数作为目标函数, 用于提高硅单晶直径检测图像的分割精度, 以弥补使用一种测度函数分割精度的不足.
2. 基于MAFSA的改进二维直方图区域斜分方法
2.1 测度函数
由于式(4)和式(5)仅能用于单阈值分割方法中, 本文现将式(4)和式(5)扩展到多阈值分割中.因为坐标转换以后, 虽然计算阈值较为简便, 但搜索变换后的直方图边界复杂度却增加了, 所以, 本文为同时计算最大类间方差和Shannon熵方便, 统一使用原坐标系.设图像被分割为$K$个区域, 分别为C0, C1, …, CK-1, 分割阈值为t=[t0, t1, …, tK-2]T, ,将式(4)拓展到多阈值分割情况时为
$ \begin{align} %\begin{aligned} {\rm tr}{\sigma}_{b}({\textbf{t}})=\sum_{i=0}^{K-1}\omega_{i}[(\mu_{ix}-\mu_{zx})^{2}+(\mu_{iy}-\mu_{zy})^{2}] %\end{aligned} \end{align} $
(6) 式中
$ \begin{align*} \omega_{i}=\sum_{(x, y)\in C_{i}}p(x, y) \end{align*} $
$ \begin{align*} \mu_{ix}=\frac{1}{\omega_{i}}\sum_{(x, y)\in C_{i}}xp(x, y) \end{align*} $
$ \begin{align*} \mu_{iy}=\frac{1}{\omega_{i}}\sum_{(x, y)\in C_{i}}yp(x, y) \end{align*} $
$ \begin{align*} \mu_{zx}=\sum_{x=0}^{L-1}\sum_{y=0}^{L-1}xp(x, y) \end{align*} $
$ \begin{align*} \mu_{zy}=\sum_{x=0}^{L-1}\sum_{y=0}^{L-1}yp(x, y) \end{align*} $
同理, 将式(5)拓展到多阈值分割情况下为
$ \begin{array}{l} {\boldsymbol{t}^*} = \\ \quad \arg \mathop {\max }\limits_{0 < t \le 2L - 2} \left\{ {\mathop {\min }\limits_{0 < t \le 2L - 2} \{ {H_0}(\boldsymbol{t}), \cdots, {H_{K - 1}}(\boldsymbol{t})\} } \right\} \end{array} $
(7) 式中
$ \begin{array}{l} {H_i}(\boldsymbol{t}) = \log [\sum\limits_{(x, y) \in {C_i}} p (x, y)] - \\ \;\;\;\;\;\;\;\;\;\;\frac{{\sum\limits_{(x, y) \in {C_i}} p (x, y)\log [p(x, y)]}}{{\sum\limits_{(x, y) \in {C_i}} p (x, y)}}, \\ \;\;\;\;\;\;\;\;i = 0, 1, \cdots, K - 1 \end{array} $
2.2 多目标人工鱼群算法
人工鱼群算法是由国内学者提出的一种模仿鱼群觅食的单目标群智能优化算法, 通过觅食、追尾、聚群和随机四种行为更新状态[1, 22].现将其扩展成为多目标优化算法, 需要解决: 1)快速构造Pareto非劣解集; 2)为人工鱼选择向导; 3)维护外部档案集.根据本文的分割问题, 对多目标人工鱼群算法的优化问题做如下定义:目标函数为maxF(xi)=[f1(xi), f2(xi)], i=1, …, N, xi=[xi1, …, xiD]T为解向量, N为人工鱼数量, D为解空间维数, S为解空间, 两条人工鱼之间的距离为dij =‖xi -xj‖, V isual表示人工鱼的感知范围, Step表示人工鱼移动的步长, Nf表示人工鱼感知范围内的人工鱼数目, L表示人工鱼最大觅食次数, rand()表示产生0~1之间的随机数.现给出本文使用Pareto最优的相关定义[18]:
定义1 (Pareto支配).对于任意两个D维向量u和v, 当且仅当任意分量ui ≥vi, i=1, …, D且至少存在uj > vj, i ∈ 1, …, D, 称作u支配v, 记作u $ \prec $ v.
定义2 (Pareto非劣解).对于x* ∈ S, 当且仅当不存在x ∈ S使得x $ \prec $ x*, 称x*为Pareto非劣解.
定义3 (Pareto非劣解集).所有Pareto非劣解组成的集合称为Pareto非劣解集.
定义4 (Pareto前沿). Pareto非劣解集对应的目标函数值所组成的集合称为Pareto前沿.
2.2.1 改进快速构造Pareto非劣解集方法
在多目标优化问题中, Pareto非劣解的选取和非劣解集的构造是至关重要的环节[17].目前大部分算法非劣解的选取都采用将备选解与其他解逐个比较, 若备选解不被其他解支配, 则为非劣解[18].逐个比较的方法虽然能够确定一个解是否为Pareto非劣解, 但算法较复杂.为此, 文献[23]提出了一种快速构造非劣解集的方法, 在判断一个解是否为非劣解时不再与其他所有解进行比较, 而是首先搜索到第一个目标函数的最优值, 然后将该目标函数值对应的第二个目标函数设置为基准值, 其他解只跟基准值比较.本文基于该方法进行改进, 加入了对多个解的第一个目标函数相等情况时的考虑, 用于构造多目标人工鱼群算法的Pareto非劣解集.现定义目标函数个数为2, 需进行比较的解一共有M个, 最大化目标函数, 构造步骤如下:
步骤1.按照第一个目标函数的函数值进行快速降序排列, 同时对每一个人工鱼按排序后的顺序进行1到M标号, 并且对第一个目标函数值相同的人工鱼进行标记.暂定第1个解为非劣解, 将其第二个目标函数值作为基准值, 令i=2．
步骤2.首先判断第i条人工鱼的第一个目标函数值跟当前基准值对应的第一个目标函数值是否相同, 若相同, 转到步骤3;否则转到步骤4.
步骤3.比较第i条人工鱼的第二个函数值与基准值大小, 如果大于基准值, 则从Pareto非劣解集中剔除当前基准值对应的非劣解, 并且将第i条人工鱼的第二个目标函数值定义为基准值, 同时将该人工鱼暂时加入Pareto非劣解集中, 然后转到步骤5;否则直接转到步骤5．
步骤4. 比较第i条人工鱼的第二个函数值与基准值大小, 如果大于基准值, 将该人工鱼的第二个目标函数值定义为基准值, 将第i条人工鱼暂时加入Pareto非劣解集中, 然后转到步骤5;否则直接转到步骤5.
步骤5.令i=i + 1, 如果i≤M, 转到步骤2;否则终止比较, 输出最终Pareto非劣解集.
2.2.2 人工鱼群算法四种行为与状态更新
人工鱼群算法通过觅食、追尾、聚群和随机四种行为更新状态, 随着目标函数的增多和Pareto非劣解定义的引入, 觅食、追尾和聚群三种行为不能直接用于多目标优化.本文定义一个拥挤度距离Dis[24], 人工鱼同时根据Pareto支配规则和拥挤度距离判断人工鱼是否向新解位置移动.当人工鱼感知范围内人工鱼数量为Nf, 拥挤度距离Dis计算步骤如下:
步骤1.初始化每个人工鱼的拥挤度距离为0, 设置当前目标数j=1.
步骤2. 对第j个目标进行降序排列, 设第i条人工鱼排序后对应的位置为i', 排序后的第1条和最后1条人工鱼的拥挤度为无穷大, 令i=2.
步骤3.从第i条人工鱼按照下式计算拥挤度距离
$ \begin{align} Dis(i)=Dis(i)+|\frac{f_{j}(i'+1)-f_{j}(i'-1)}{\max(f_{j})-\min(f_{j})}| \end{align} $
(8) 步骤4.令i=i + 1, 判断i是否大于Nf, 若是, 转到步骤5;否则转到步骤3.
步骤5.令j= j +1, 判断j 是否大于2, 若是, 结束计算, 输出最后结果; 否则转到步骤2.
多目标人工鱼群优化算法的四种行为分别为:
1) 觅食行为.假设人工鱼的当前状态为x, 根据下面两式选取一个新的状态x', 如果F(x')支配F(x), 则 x移动到x', 否则人工鱼继续根据下面两式选择新的状态.如果人工鱼连续L次都没有满足更新状态的条件, 则执行随机行为.
$ {\mathit{\boldsymbol{x}}_V} = \mathit{\boldsymbol{x}} + Visual \cdot {\rm{rand}}\left( {} \right) $
(9) $ \mathit{\boldsymbol{x'}} = \mathit{\boldsymbol{x}} + \frac{{{\mathit{\boldsymbol{x}}_V} - \mathit{\boldsymbol{x}}}}{{\left\| {{\mathit{\boldsymbol{x}}_V} - \mathit{\boldsymbol{x}}} \right\|}} \cdot step \cdot {\rm{rand}}() $
(10) 2) 追尾行为.假设人工鱼的当前状态为x, 搜索其感知范围内的非支配解并计算拥挤度距离, 拥挤度距离最大的解为xgbest, 如果Dis(xgbest) > Dis(x), 则根据下式进行状态更新, 否则执行聚群行为.
$ \mathit{\boldsymbol{x'}} = \mathit{\boldsymbol{x}} + \frac{{{\mathit{\boldsymbol{x}}_{gbest}} - \mathit{\boldsymbol{x}}}}{{\left\| {{\mathit{\boldsymbol{x}}_{gbest}} - \mathit{\boldsymbol{x}}} \right\|}} \cdot step \cdot {\rm{rand}}() $
(11) 3) 聚群行为.假设人工鱼的当前状态为x, 在其感知范围内所有人工鱼的中心位置为xc, 如果F(xc)支配F(x), 并且Dis(xc) > Dis(x), 则根据下式进行状态更新, 否则执行觅食行为.
$ \mathit{\boldsymbol{x'}} = \mathit{\boldsymbol{x}} + \frac{{{\mathit{\boldsymbol{x}}_c} - \mathit{\boldsymbol{x}}}}{{\left\| {{\mathit{\boldsymbol{x}}_c} - \mathit{\boldsymbol{x}}} \right\|}} \cdot step \cdot {\rm{rand}}() $
(12) 4) 随机行为.人工鱼根据式(9)随机更新到一个新的状态.
2.2.3 外部档案集维护策略
外部档案集的更新和维护是保证鱼群在更新状态时能够得到具有更好分布性和多样性的Pareto前沿.本文采用基于格密度的维护策略[18, 23], 对每次外部档案集的每维目标空间分割成10块, 则整个目标空间分割成100块, 将每条人工鱼的格密度定义为该人工鱼所在区域所包含的人工鱼数.当每迭代一次以后, 将该次得到的Pareto非劣解集与外部档案集合并, 并且按照第2.2.1节中改进快速构造Pareto非劣解集方法将被支配的解剔除, 然后判断最后的解的数量是否大于外部档案集的最大数量Km, 如果大于Km, 计算每个解的格密度, 剔除格密度最大的解, 直到解的数量等于外部档案集的最大数量.
2.3 基于MAFSA的二维直方图区域斜分方法步骤
基于MAFSA的二维直方图区域斜分方法算法步骤如下:
步骤1.初始化.对人工鱼群算法参数初始化, 包括人工鱼的数量N, 人工鱼的初始位置xi=[ti, 0, ti, 1, …, ti, K-2]T, i=1, …, N, ,感知范围Visual初值, 感知范围最小值Visualmin, 最大觅食次数L, 人工鱼移动的步长Step初值, 步长最小值Stepmin, 最大迭代次数T, 外部档案集的最大数量Km, 外部档案集为空, 目标函数f1为式(6), f2为式(7).
步骤2.根据式(6)和式(7)分别计算目标函数f1和f2, 根据第2.2.1节中选取非劣解的方法构造Pareto非劣解集, 根据第2.2.3节中方法将Pareto非劣解集与外部档案集合并.
步骤3. 搜索每条人工鱼感知范围内的人工鱼数Nf, 如果Nf < 2, 执行觅食行为; 否则执行追尾行为.
步骤4. 根据式(6)和式(7)分别计算目标函数f1和f2, 根据第2.2.1节中选取非劣解的方法构造本次迭代Pareto非劣解集, 根据第2.2.3节中方法将Pareto非劣解集与外部档案集合并.
步骤5.根据下式更新人工鱼的感知范围Visual和人工鱼移动的步长Step.
$ \phi =\exp(-30\ast (t/T)^{5}) $
(13) $ Step=Step\cdot\phi+Step_{\rm{min}} $
(14) $ Visual=Visual\cdot\phi+Visual_{\rm{min}} $
(15) 步骤6. 如果达到最大迭代次数T, 输出外部档案集, 并且输出格密度最小的解作为最优解; 否则转到步骤3.
3. 实验与结果分析
为验证本文所提多目标人工鱼群算法和二维直方图区域斜分多阈值分割方法的有效性, 设计了2个蒙特卡洛对比实验:第一个是通过标准测试函数集ZDT对多目标人工鱼群算法进行验证, 第二个是通过具有代表性的两幅不同阶段的硅单晶直径检测图像对本文所提二维直方图区域斜分多阈值分割方法进行验证.整个实验在Intel Xeon(R) 2.00GHz(6*12核), 128GB内存, 操作系统为Linux 64的高性能计算集群上进行, 编程环境为Matlab R2013a, 为保证实验在同条件下执行, 所有程序采用单核串行执行, 独立重复实验200次.
3.1 多目标测试函数实验
为验证本文所提多目标人工鱼群算法的有效性, 采用常用的ZDT标准函数集进行测试, ZDT标准函数集如表 1所示[21], 并且与NSGA-Ⅱ[17]、p-optimality criteria evolutionary algorithm (p-OCEA)[25]两种多目标进化算法和多策略差分进化的元胞粒子群算法(Cellular multi-objective particle swarm based on multi-strategy differential evolution, MPSOCell)[20]、基于Pareto熵的多目标粒子群优化算法(Multiobjective particle swarm optimization based on pareto entropy, peMOPSO)[21]两种改进多目标粒子群优化算法进行对比.性能指标采用反转世代距离(Inverted generational distance, IGD)[21]和运行时间, IGD是度量真实Pareto前沿与算法得到的近似Pareto前沿的距离指标, 定义如下
$ IGD=\frac{1}{|P|}\sum_{i=1}^{|P|}Dist_{i} $
(16) $ Dis{t_i} = {\min _{j = 1, \cdots, |A|}}\sqrt {\sum\limits_{m = 1}^M {\left( {\frac{{{f_m}({\mathit{\boldsymbol{p}}_i}) - {f_m}({\mathit{\boldsymbol{a}}_j})}}{{f_m^{{\rm{max}}} - f_m^{{\rm{min}}}}}} \right)} } $
(17) 表 1 ZDT标准测试函数集Table 1 ZDT standard test function set函数标号 目标函数 变量数 变量范围 采样点数 ZDT1 30 xi ∈ [0, 1] 1 000 ZDT2 30 xi ∈ [0, 1] 1 000 ZDT3 30 xi ∈ [0, 1] 1 000 ZDT4 10 xi ∈ [0, 1]
xi ∈ [-5, 5]
i=2, …, n1 000 ZDT6 10 xi ∈ [0, 1] 1 000 式中, P为真实Pareto非劣解集, A是算法得到的近似Pareto非劣解集, fm max和fm min分别为P中第m个目标的最大值和最小值, m=1, 2为目标数, pi ∈ P, i=1, 2, ·…, |P|, aj ∈ A, j=1, 2, …, |A|, |P|和|A|分别为真实Pareto非劣解集采样样本数和算法得到的Pareto非劣解个数．
算法参数设置如下:所有算法的初始种群大小N=100, 最大迭代次数T=1000, 最大外部档案集数Km=100, 初值范围和种群维数对应表 1的变量范围和变量数.本文方法的其他参数设置如下:人工鱼的感知范围Visual=100, 感知范围最小值Visualmin=2, 最大觅食次数L=10, 人工鱼移动的步长为初值范围, 步长最小值Stepmin=0.01.另外四种对比方法的其他参数设置均与参考文献一致.
实验结果如表 2所示, 由表 2可得, 在对前3个测试函数的实验中, 本文方法、NSGA-Ⅱ、p-OCEA和MPSOCell的IGD较为稳定且接近, 而peMOPSO在ZDT2函数上的IGD略大.对于多模态ZDT4函数, MPSOCell的IGD较大, 本文方法与其他三种方法较为接近.而对于ZDT6函数, 由于ZDT6函数具有非均匀且不连续的特点, NSGA-Ⅱ、p-OCEA、MPSOCell和peMOPSO的IGD显著增大, 本文方法的IGD仅略微增大, 优于其他四种方法.综合前4个测试函数的结果, 本文方法与NSGA-II和p-OCEA的IGD和运行时间较为稳定, MPSOCell和peMOPSO的IGD和运行时间均波动较大.由于本文采用了快速构造Pareto非劣解的方法, 对5个测试函数的运行时间明显小于其他四种方法.因此, 根据上述实验结果分析可得, 本文方法在对Pareto前沿搜索的稳定性、准确性和实时性要优于其他四种方法.
表 2 多种算法对比实验结果Table 2 Experimental results of different algorithms对比方法 NSGA-Ⅱ p-OCEA peMOPSO MPSOCell 本文方法 测试函数 IGD 时间(s) IGD 时间(s) IGD 时间(s) IGD 时间(s) IGD 时间(s) ZDT1 0.0123 7 131.89 0.0826 6 288.28 0.0100 5 689.27 0.0172 2 599.12 0.0088 986.87 ZDT2 0.0207 6 367.91 0.0977 6 350.72 0.4995 3 190.81 0.0088 2 598.90 0.0211 829.74 ZDT3 0.0060 6 167.51 0.0532 6 531.65 0.0087 2 008.90 0.0116 1 285.60 0.0197 987.56 ZDT4 0.0354 6 443.92 0.0191 6 325.01 0.0138 1 107.43 7.5031 1 740.71 0.0597 551.46 ZDT6 96.8043 6 391.91 95.9782 6 656.29 71.1176 1 405.68 86.0420 917.85 0.4865 651.83 3.2 硅单晶直径检测图像分割实验
为验证本文所提出的二维直方图区域斜分多阈值分割方法分割硅单晶直径检测图像的有效性, 将本文方法应用到TDR-150型单晶炉直径检测工程实验的图像中. TDR-150型单晶炉为极大规模集成电路用硅单晶的全自动生长设备, 投料量为170kg, 控制直径为300mm, 等径长度为850mm, 测量原理如图 1所示, TDR-150型单晶炉和传感器装置如图 6所示, 直径测量采集图像分别采用晶体等径长度为10mm和等径长度为400mm的两幅图, 如图 7所示.
图 6 TDR-150单晶炉直径检测装置实物图Fig. 6 Photos of TDR-150 single crystal silicon furnace diameter detecting device3.2.1 分割实验结果
为验证本文方法的有效性, 采用直分Otsu法[12](下文称方法1)、斜分Otsu法[13](下文称方法2)和改进斜分Otsu法[16](下文称方法3)三种方法进行对比实验, 由于前两种方法仅适用于单阈值分割, 本文分别将其利用人工鱼群算法扩展为多阈值分割, 而方法3的测度函数选用文献[16]改进的测度函数, 搜索算法采用人工鱼群算法.
实验算法参数设置如下:所有方法的人工鱼数N=30, 感知范围Visual=1000, 感知范围最小值Visualmin=2, 最大觅食次数L=10, 人工鱼移动的步长初值为初值范围, 步长最小值Stepmin=0.01, 最大迭代次数T=100, 方法1初值范围为[0, 255], 其他三种方法的初值范围为[0, 510], 本文方法中最大外部档案集数Km=100.
为能够更明显地比较对高亮光环的分割结果, 本文将每个实验结果用两种形式表示, 图 8为在等径10mm时四种方法多阈值分割结果和高亮光环局部放大图像, 图 9为在等径10mm时四种方法将多阈值分割结果以高亮光环为目标、其余区域为背景的二值化图像表示.由图 8和图 9可知, 由于方法1仅考虑主对角线附近的区域, 错分率较高, 高亮光环区域中噪声点较多, 并且将硅熔液的部分区域错分到高亮光环的区域中; 方法2采用了区域斜分分割方法, 分割效果较方法1明显有所提高, 但因为晶线的亮度与高亮光环接近, 将晶线与高亮光环分割为一个区, 方法3引入了后处理方法, 效果与方法2较接近, 略优于方法2;本文方法由于采用了多目标优化的方法, 能够较为准确地将高亮光环分割出来.图 8和图 9表明本文方法在等径10mm时的分割精度优于其他三种方法.
图 8 等径10 mm时四种方法的多阈值分割结果Fig. 8 Results of different methods for multi-threshold image segmentation at 10 mm length in body growth
图 9 等径10 mm时四种方法的二值化图像Fig. 9 Binarization image of different methods at 10 mm length in body growth等径400mm的实验结果如图 10和图 11所示, 图 10为四种方法多阈值分割结果和高亮光环局部放大图像, 图 11为四种方法将以高亮光环为目标、其余区域为背景的二值化图像表示.从图 10可看出, 在等径400mm时, 随着硅单晶棒长度增加, 图像的区域也发生了变化, 晶线对分割结果的影响基本消除.方法1、方法2和方法3的分割结果都没有出现晶线的错分区域.但方法1仍将硅熔液的部分区域错分到高亮光环的区域中; 而方法2、方法3和本文方法结果基本一致, 都能分割出高亮光环.由图 11可得, 方法1在硅熔液区域包含一部分噪声点, 方法2和方法3在硅熔液区域也包含少量噪声, 方法3的噪声点数量少于方法2;而本文方法将高亮光环能够准确分割出来, 基本不含噪声点. 图 10和图 11表明了本文所提基于多目标人工鱼群优化的二维直方图区域斜分多阈值分割方法在等径400mm时的有效性.
图 10 等径400 mm时四种方法的多阈值分割结果Fig. 10 Results of different methods for multi-threshold image segmentation at 400 mm length in body growth
图 11 等径400 mm时四种方法的二值化图像Fig. 11 Binarization image of different methods at 400 mm length in body growth3.2.2 分割效果评价
为避免对以上图像分割结果视觉分析的主观性, 本文选用信噪比、分类误差率两个量化指标对分割精度进行评价[16], 同时利用运行时间对算法实时性进行评价.信噪比计算公式如下.
$ PSNR=20\lg\left(\frac{255}{RMSE}\right) $
(18) 其中,
$ RMSE=\sqrt{\frac{\sum_{i=1}^{m}\sum_{j=1}^{n}(I(i, j)-\hat{I}(i, j))^{2}}{mn}} $
式中, I和$\hat{I}$分别为模板图像(一般由人工根据先验知识分割所得)和算法分割结果图像, 图像大小为m × n, 式(18)说明分割结果与模板图像越相似, 信噪比结果越大, 分割结果越准确.
分类误差率定义如下
$ ME\!=\!1\!-\!\frac{|I_{0}\cap \hat{I}_{0}|\!+\!|I_{1}\cap \hat{I}_{1}|\!+\!\cdots\!+\!|I_{K-1}\cap \hat{I}_{K-1}|}{|I|} $
(19) 式中, $\hat{I}_{i}$为分割后第i个区域, 0≤i≤K -1, $|\bullet|$表示区域$\bullet$的势.式(19)说明分类误差率越小, 分割结果越接近模板图像, 分割效果越准确.
表 3和表 4为各方法分割结果的性能指标.由表 3和表 4可得, 在等径10mm和等径400mm时, 方法1的信噪比最小, 分类误差率最大, 而方法2采用区域斜分法, 分类误差率、信噪比和运行时间较方法1具有较大改善, 方法3改进了区域斜分法的测度函数, 降低了计算维度, 运行时间明显降低, 信噪比和分类误差率都略优于方法2.本文方法由于采用了多目标优化的方法, 信噪比和分类误差率均有较明显的改善, 信噪比最大, 分类误差率最小, 但算法运行时间较方法2和方法3略长.由于硅单晶生长是一个缓慢的过程, 本文方法的运行时间能够满足直径检测和控制工艺的实时性需求. 表 3和表 4与上节实验中的图像分割结果分析基本一致, 充分验证了本文方法的有效性和上节实验的结论, 对提高硅单晶直径检测精度具有重要意义.
表 3 等径10 mm时各种方法分割结果的性能指标Table 3 Results of efficiency among different methods for image segmentation at 10 mm length in body growth方法 PSNR (dB) 分类误差率（％) 运行时间（s) 方法1 12.8480 0.8008 11.0274 方法2 30.0228 0.1849 4.6387 方法3 31.0546 0.1772 1.8061 本文方法 35.1589 0.0626 7.6084 表 4 等径400 mm时各种方法分割结果的性能指标Table 4 Results of efficiency among different methods for image segmentation at 400 mm length in body growth方法 PSNR (dB) 分类误差率（％) 运行时间（s) 方法1 13.9686 0.6600 12.1390 方法2 32.1443 0.1042 4.6601 方法3 32.0341 0.0944 2.5306 本文方法 34.8116 0.0826 7.5028 4. 结论
为提高硅单晶的直径检测图像的分割精度, 本文提出了以最大类间方差和最大熵同时作为测度函数, 利用多目标优化的方法对图像二维直方图进行多阈值区域斜分.
设计了一种多目标人工鱼群算法, 改进了快速构造Pareto非劣解集的方法, 并将该方法和选择群体进化向导、外部档案集维护的方法引入到人工鱼群算法中.通过标准测试函数集进行测试, 并且与NSGA-Ⅱ、p-OCEA、MPSOCell、peMOPSO进行对比实验, 验证该算法具有较高的搜索精度.
提出一种基于多目标人工鱼群算法的二维直方图区域斜分多阈值分割方法, 将多目标人工鱼群算法引入二维直方图区域斜分的多阈值搜索中.分别利用等径10mm和等径400mm时的硅单晶直径检测图像进行测试.同时以PSNR和分类误差率作为衡量标准, 与二维直方图区域直分法、二维直方图区域斜分Otsu法和改进二维直方图区域斜分Otsu法对比.实验结果表明, 由于多目标优化方法的引入, 虽然略微增加了算法运行时间, 但使得图像分割的信噪比和分类误差率有较大的改善, 对提高硅单晶直径检测精度具有重要意义.
-
图 1 OMP、StOMP、CoSaOMP、ROMP和COCLMP在噪声干扰下对稀疏信号的重构效果
Fig. 1 Reconstruction effect of sparse signals include noise for OMP, StOMP, CoSaOMP, ROMP, COCLMP
图 2 噪声干扰的不同采样次数下OMP、StOMP、CoSaOMP、ROMP、COCLMP算法重构结果
Fig. 2 Reconstruction result include noise under different sampling number for OMP, StOMP, CoSaOMP, ROMP, COCLMP
图 3 噪声干扰下OMP、StOMP、CoSaOMP、ROMP、COCLMP算法对图像的重构效果
Fig. 3 Image reconstruction effect include noise for OMP, StOMP, CoSaOMP, ROMP, COCLMP
-
[1] Moshtaghpour A, Jacques L, Cambareri V, Degraux K, De Vleeschouwer C. Consistent basis pursuit for signal and matrix estimates in quantized compressed sensing. IEEE Signal Processing Letters, 2016, 23(1):25-29 doi: 10.1109/LSP.2015.2497543 [2] Zhao W Q, Beach T H, Rezgui Y. Efficient least angle regression for identification of linear-in-the-parameters models. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 2017, 473(2198):Article No.20160775 doi: 10.1098/rspa.2016.0775 [3] Lee H C, Song B, Kim J S, Jung J J, Li H H, Mutic S, et al. An efficient iterative CBCT reconstruction approach using gradient projection sparse reconstruction algorithm. Oncotarget, 2016, 7(52):87342-87350 http://europepmc.org/abstract/MED/27894103 [4] 许志强.压缩感知.中国科学:数学, 2012, 42(9):865-877 https://www.wenkuxiazai.com/doc/3851fb1eff00bed5b9f31d8b.htmlXu Zhi-Qiang. Compressed sensing:a survey. Science China:Mathematics, 2012, 42(9):865-877 https://www.wenkuxiazai.com/doc/3851fb1eff00bed5b9f31d8b.html [5] 方红, 杨海蓉.贪婪算法与压缩感知理论.自动化学报, 2011, 37(12):1413-1421 http://www.aas.net.cn/CN/abstract/abstract17639.shtmlFang Hong, Yang Hai-Rong. Greedy algorithms and compressed sensing. Acta Automatica Sinica, 2011, 37(12):1413-1421 http://www.aas.net.cn/CN/abstract/abstract17639.shtml [6] Cohen A, Dahmen W, DeVore R. Orthogonal matching pursuit under the restricted isometry property. Constructive Approximation, 2017, 45(1):113-127 doi: 10.1007/s00365-016-9338-2 [7] Do T T, Gan L, Nguyen N, Tran T D. Sparsity adaptive matching pursuit algorithm for practical compressed sensing. In: Proceedings of the 42th Asilomar Conference on Signals, Systems and Computers. Pacific Grove, CA, USA: IEEE, 2008. 581-587 [8] Needell D, Vershynin R. Signal recovery from incomplete and inaccurate measurements via regularized orthogonal matching pursuit. IEEE Journal of Selected Topics in Signal Processing, 2010, 4(2):310-316 doi: 10.1109/JSTSP.2010.2042412 [9] Lee D. MIMO OFDM channel estimation via block stagewise orthogonal matching pursuit. IEEE Communications Letters, 2016, 20(10):2115-2118 doi: 10.1109/LCOMM.2016.2594059 [10] Davenport M A, Needell D, Wakin M B. Signal space CoSaMP for sparse recovery with redundant dictionaries. IEEE Transactions on Information Theory, 2013, 59(10):6820-6829 doi: 10.1109/TIT.2013.2273491 [11] Giryes R, Elad M. RIP-based near-oracle performance guarantees for SP, CoSaMP, and IHT. IEEE Transactions on Signal Processing, 2012, 60(3):1465-1468 doi: 10.1109/TSP.2011.2174985 [12] Ramirez-Giraldo J C, Trzasko J, Leng S, Yu L, Manduca A, McCollough C H. Nonconvex prior image constrained compressed sensing (NCPICCS):theory and simulations on perfusion CT. Medical Physics, 2011, 38(4):2157-2167 doi: 10.1118/1.3560878 [13] Babaie-Kafaki S, Ghanbari R. A hybridization of the Hestenes-Stiefel and Dai-Yuan conjugate gradient methods based on a least-squares approach. Optimization Methods and Software, 2015, 30(4):673-681 doi: 10.1080/10556788.2014.966825 [14] Shaw C B, Prakash J, Pramanik M, Yalavarthy P K. Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography. Journal of Biomedical Optics, 2013, 18(8):Article No.80501 doi: 10.1117/1.JBO.18.8.080501 [15] Paige C C, S M A. LSQR:an algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software, 1982, 8(1):43-71 doi: 10.1145/355984.355989 [16] Cichocki A, Zdunek R, Amari S I. Csiszár's divergences for non-negative matrix factorization: family of new algorithms. In: Proceedings of the 2006 International Conference on Independent Component Analysis and Blind Signal Separation. Charleston, SC, USA: Springer, 2006. 32-39 [17] Cichocki A, Amari S I. Families of Alpha-Beta-and Gamma-divergences:flexible and robust measures of similarities. Entropy, 2010, 12(6):1532-1568 doi: 10.3390/e12061532 [18] 余南南, 邱天爽.压缩传感条件下红外和可见光图像融合技术的研究.信号处理, 2012, 28(5):692-698 http://cdmd.cnki.com.cn/Article/CDMD-10183-2008126444.htmYu Nan-Nan, Qiu Tian-Shuang. Fusion technology of infrared and visible images in compressive sensing. Signal Processing, 2012, 28(5):692-698 http://cdmd.cnki.com.cn/Article/CDMD-10183-2008126444.htm [19] 王琴, 沈远彤.基于压缩感知的多尺度最小二乘支持向量机.自动化学报, 2016, 42(4):631-640 http://www.aas.net.cn/CN/abstract/abstract18849.shtmlWang Qin, Shen Yuan-Tong. Multi-scale least squares support vector machine using compressive sensing. Acta Automatica Sinica, 2016, 42(4):631-640 http://www.aas.net.cn/CN/abstract/abstract18849.shtml 期刊类型引用(4)
1. 郭俊锋，胡婧怡，王智明. 基于压缩感知的缺失机械振动信号重构新方法. 振动与冲击. 2024(10): 197-204 . 
百度学术2. 田占峰，裴世建，于强. 复杂岩溶地区跨孔地震CT反演成像研究——以贵阳在建地铁3号线为例. CT理论与应用研究. 2024(05): 561-567 . 百度学术3. 焦帆，曹余庆，谢莉. 时滞和阶次未知的双率采样输出误差系统辨识. 控制与决策. 2024(09): 3006-3012 . 百度学术4. 王金甲，张玉珍，夏静，王凤嫔. 多层局部块坐标下降法及其驱动的分类重构网络. 自动化学报. 2020(12): 2647-2661 . 本站查看其他类型引用(5)
-
下载:
下载:
计量
- 文章访问数: 2699
- HTML全文浏览量: 453
- PDF下载量: 706
- 被引次数: 9
