一种多策略改进遗传算法的医学图像配准方法

秦泽青，叶志伟，王泽松，徐 川，曹 羽

1(湖北工业大学 计算机学院，武汉 430068) 2(湖北大学 计算机与信息工程学院，武汉 430062) 3(湖南大学 信息科学与工程学院，长沙 410012) E-mail：weizhiye121@163.com

1 引 言

医学影像可以分为解剖成像和功能成像[1]，前者分辨率较高，能够提供解剖形态信息，后者分辨率较低，但可以很好地提供功能代谢等信息.医学图像配准是指通过调整几何变换的参数，从而使不同类型医学影像的共同部分达到空间位置上一致的过程[2，3]，可以更全面地提供病灶的信息，对现代医疗有着重大的意义[4].本质上，医学图像配准是一种参数最优化问题[5]，可以使用优化算法进行处理.

传统的优化方法已被应用于医学图像配准问题，然而，由于图像离散化，噪声等因素的影响，优化过程容易陷入局部最优解.近年来，基于进化算法和其他元启发式方法的医学图像配准方法展现出了良好的应用潜力[6]，配准精确度高，健壮性强，在很大程度上克服了传统方法的不足[7].

遗传算法(Genetic Algorithm，GA)是经典的进化算法，基于进化论模拟了自然界中的染色体交叉，碱基对变异等现象，有着良好的全局搜索能力，在许多领域得到了广泛的应用[8].但是，基于遗传算法的医学图像配准方法仍然存在着局部搜索能力较差且容易早熟收敛等不足[9].

为了改善上述问题，本研究从随机机制，种群多样性以及遗传算子3个方面出发，提出多种策略对遗传算法进行改进：1)设计一种基于混沌系统与频度记忆的随机策略，利用Logistic混沌映射产生随机数序列并结合频度记忆对随机个体的生成进行了干预，实现了随机数和随机个体的可再现性与可控性；2)提出一种动态的基因多样性控制策略，通过控制进化过程中个体间基因的差异，保证了种群的多样性；3)提出一种基于杂种优势的进化策略，模拟三系法杂交水稻协同进化育种模式设计新的遗传算子.最后结合多分辨率策略，利用改进的遗传算法搜索多模态医学图像配准的最优参数并进行实验仿真.结果表明本文提出的策略有效地提升了遗传算法的搜索能力并降低了早熟收敛的概率.

2 相关理论

2.1 图像配准

图像配准指的是通过特定的几何变换，从而使两个或多个图像的共同部分达到空间位置上一致的过程[2，3].在这个过程中，作为参照模板固定的图像叫做参考图像，被几何变换调整的图像叫做移动图像.图像配准可以分为基于灰度的配准方式以及基于特征的配准方式.从效率来看，前者只需要提取图像的特征点、边缘线等突出信息，时耗较小，而后者需要提取图像的所有信息，效率较低；从复杂程度来看，前者需要对图像进行特征提取等预处理，较为繁琐，而后者无需任何预处理，十分简洁；从应用范围来看，基于特征的图像配准只适用于有着明显特征的图片，而基于灰度的图像配准适用于所有类型的图像[2].本文采用基于灰度的图像配准方式，包含5个步骤：初始化参数过程、几何变换过程、评价过程、优化过程、输出过程.配准流程如图1所示.

图1 基于灰度的图像配准流程图Fig.1 Intensity-based image registration flow chart

图像配准作为一种有效的信息整合方法，已经在医学图像引导手术(image-guided surgery)[4]等方面得到了广泛的应用，可以被视作参数优化问题.基于传统优化算法的图像配准方法已取得了一定的成效，例如鲍威尔方法(Powell′s method)，下降单纯型法(downhill simplex method)，高斯牛顿法(Gauss-Newton)，梯度下降法(gradient descent)等[2]，但由于图像的离散化，噪声等负面因素的影响，这些方法非常容易陷入局部最优解[3].

近年来，进化计算和元启发式算法在图像配准问题上受到了广泛地关注.相较于传统优化算法，这些方法求解速度快、寻优能力强、计算成本低.Valsecchi等人[10]提出了一种改进后的散射搜索方法(Scatter Search*)，采用双层设计(2-tier design)同时考虑了最优子集的质量和多样性，并按照图像配准的特征对散射搜索的基本模板进行改进，最后将改进后的散射搜索法与多分辨率进行结合，表现出较好的性能.Bermejo等人[11]首次将珊瑚礁优化算法(Coral Reef Optimization)用于医学图像配准问题，基于5种算子结合的进化模式为珊瑚设计了附着层，在基于灰度和特征的医学图像配准实验中都表现出了较高的配准精度和较强的健壮性.Chen等人[12]将生物地理学优化算法(Biogeography-based Optimization)首次用于基于灰度的多模医学图像配准问题，并结合精英学习策略对其做出了改进，提升了所得解集的质量并提高了收敛的效率.

上述算法都在医学图像配准问题上表现出了良好的寻优能力.然而，这些方法一定程度上还存在容易陷入局部最优解的问题，收敛性也没有得到充分证明.作为经典的最优化算法，遗传算法具有着简单高效的特性，其收敛性也得到了证明，本文拟通过多种策略对遗传算法做出改进，使其更加适用于医学图像配准问题.

2.2 几何变换

几何变换分为刚体变换和非刚体变换两类，后者又包括仿射变换以及透视变换等.利用几何变换可以对图像进行旋转、平移、缩放、剪切等操作.为了减少计算代价，本文采用自由度最低的三维刚体变换.此变换包含3个旋转参数rx、ry、rz和3个平移参数tx、ty、tz.

2.3 遗传算法

遗传算法是一种基于生物界自然选择和遗传机制构建的随机搜索算法，由美国的Holland教授提出.这种算法来源于达尔文的进化论、魏茨曼的物种选择学说和孟德尔的群体遗传学说，模拟了自然界生物繁育过程中的染色体交叉和基因变异等现象，有着全局优化性、梯度信息不依赖性、简单易施性等特点，在很多领域有着广泛的应用[8].作为一种自适应全局优化算法，遗传算法使用二进制编码，其基本执行过程为：

步骤1.确定种群规模N，交叉概率Pc以及变异概率Pm，设定停止准则；

步骤2.按照某一随机机制初始化N个个体的基因；

步骤3.计算各个体的适应度，按照选择算子选择亲代；

步骤4.按照交叉算子遵循概率Pc执行交叉；

步骤5.按照变异算子执行概率Pm执行变异；

步骤6.计算各个体的适应度，按照选择算子进行淘汰；

步骤7.终止检验.若已经满足停止准则，则输出目前种群中适应度最大的个体基因作为最优解；否则转步骤3.

针对遗传算法的研究主要包含3个部分：随机机制(Random Mechanisms)、遗传算子(Genetic Operators)、种群多样性(Diversity of Population).其中遗传算子包括了选择、交叉和变异3种基本形式.

在医学图像配准问题中，遗传算法存在着编码方式受限、局部搜索能力较差、容易早熟收敛等不足[9].图2表明了这些不足与本文主要研究内容之间的关系.

图2 本文主要的研究内容Fig.2 Topics of this work

3 提出的方法

为了让遗传算法更加适用于基于灰度的医学图像配准问题，本研究从随机机制、种群多样性、遗传算子3个研究方向对标准的遗传算法提出多种策略进行改进并用于医学图像配准.

3.1 随机机制

遗传算法是一种随机搜索算法，随机数和随机个体的产生机制在很大程度上决定了搜索的精度和速度.在标准遗传算法的实现中，大多使用类如randomize()或rand()等随机函数，存在着效率低，不可再现等问题.本文将混沌系统和频度记忆结合，实现了随机数和随机个体的可再现性与可控性.

3.1.1 混沌系统

混沌系统是一种复杂的非线性动力学系统，产生的混沌序列具有着遍历性、类随机性、确定可再生性等特点[13].混沌系统被广泛应用于保密通信、数据安全等领域中并取得了良好的效果.本文采用混沌系统中简单的一维Logistic混沌映射来生成随机数，如式(1)所示：

xn+1=r·xn·(1-xn)

(1)

其中：系统控制参数r∈(0，4]，x∈(0，1).

Logistic映射的动态行为与系统控制参数r密切相关，不同的r可以使Logistic映射处于周期性或混沌状态.研究[14]表明，当r∈(3.5699456，4]时，Logistic映射进入混沌状态，产生的混沌序列{xn|n=0，1，2，3，…}呈现出非周期性，具有良好的随机分布特征且不收敛.在r=4时，Logistic映射的输入输出都在(0，1)区间内，混沌轨迹最为明显，具有极好的遍历特性.本文采用r=4时的Logistic混沌映射来产生随机序列，对应的映射表达式如式(2)所示：

xn+1=4·xn·(1-xn)

(2)

3.1.2 频度记忆

在禁忌搜索的实现中，频度记忆(frequency-based memory)被广泛应用于长期记忆策略[15]，这种方法可以自适应地调整随机个体生成的概率.本文利用频度记忆对遗传算法中随机个体的生成过程进行干预，使其分布更加均匀，从而提升随机个体的可控性.基于频度记忆的随机个体生成策略分为3步：

步骤1.将随机个体参数的搜索空间S从小到大划分为4个大小相等的区间S={s1，s2，s3，s4}，并赋予每个区间一个频率计数器freqi，初始值记为0；

步骤2.随机生成一个属于[0，4]的整数来选择具体的区间.随机整数的生成遵循Logistic混沌映射产生的混沌序列，但是其生成概率与freqi的值成反比，生成后对应的freqi值加1；

步骤3.在被选择的区间中利用Logistic混沌映射生成随机个体的参数值.

3.2 种群多样性

早熟收敛是指超级个体的出现使得种群一致向其靠拢的现象.由于拥有相对良好的性状，超级个体的基因会很快在种群中占有绝对的比例，从而导致多样性迅速降低，整个群体丧失进化能力.为了降低这一现象发生的概率，本研究提出了一种动态的基因多样性控制器.设个体基因的维度为n，我们要求进化过程中子代基因x与亲代基因y的距离不得小于多样性阈值Diversity_Threshold，如式(3)所示：

(3)

由于在进化的过程中群体的需求是变化的，所以多样性控制器的阈值应该是动态的.自然界中物种的进化现象表明，群体进化初期个体之间性状差异较大，而中后期为了适应于特定的环境，群体的性状趋于一致[16].所以进化的早期不需要多样性控制或者说只需要很弱的控制，过早的控制可能会丢失宝贵的临近解；然而在进化的中后期则需要加大控制的力度来避免早熟收敛.基于上述思想，多样性控制器的阈值Diversity_Threshold是随着种群的进化逐渐增大的.假设种群的进化次数为n，最大差异限度为Diversity_Max.Diversity_Threshold是自变量处于[1，n]，函数值处于[0，Diversity_Max]上的线性函数.当种群进化次数为1时Diversity_Threshold=0，当进化次数为n时Diversity_Threshold=Diversity_Max.

3.3 遗传算子

为了进一步学习自然界中优秀的繁育机制并扩大基因差异化对种群进化的正向影响，受到杂种优势的启发，本研究对标准的遗传算子进行了改进.杂种优势是指杂合体在一种或多种性状上优于亲本的一种生物学现象，在农业中得到了广泛的应用[17].中国的三系杂交水稻[18]是杂种优势的典型应用案例，在很大程度上解决了粮食短缺的问题.三系是指保持系、恢复系以及不育系，在繁育的过程中水稻种子被人为地划分进这三系.位于不同系的水稻种子通过杂交、自交等方式相互学习优良的性状并繁育后代.相较于传统遗传算法的双亲学习模式，中国的三系杂交水稻这种三系协同进化的模式更好地发挥群体智能的优势并在很大程度上保证了群体基因的差异.本研究模拟了中国三系杂交水稻的育种机制式，改进了标准遗传算法的交叉、变异以及选择算子.具体的步骤如下：

步骤1.利用随机机制初始化n个遗传个体X的基因，并按照适应度函数值从高到低将遗传个体等量地划分入保持系{X1，X2，…Xm}、恢复系{Xm+1，Xm+2，…X2m}和不育系{X2m+1，X2m+2，…X3m}，每系中的个体数量m=⎣n/3」；

步骤2.从保持系和不育系中分别选择亲代进行杂交，对于保持系的个体Xi选择不育系的个体Xn-i与其配对.相较于传统遗传算法的单点交叉模式，本研究采用BLX-α交叉算子来对亲代基因进行学习.设被选中的亲代为x、y，对于他们基因的的每一个维度i有xi，yi∈[ai，bi]，其中ai是维度i取值范围的下界，bi是取值范围的上界.子代基因zi是在区间[LEFTi，RIGHTi]上服从均匀分布的随机数.LEFTi是区间下界，见式(4).RIGHTi是区间上界，见式(5)：

LEFTi=max{ai，min(xi，yi)-αdi}

(4)

RIGHTi=min{max(xi，yi)+αdi，bi}

(5)

其中：亲代基因第i维的距离di=|xi-yi|，α为杂交率.

如果子代的适应度优于不育系亲代，则后者的基因将会被直接替换；

步骤3.恢复系的个体充分学习保持系中表现最好个体的基因.在这一步，基于随机机制选择两个恢复系的个体，并选择保持系中适应度最好的个体.这些个体将进行交叉操作完成自交.设来自保持系的个体为Zbest，恢复系中的个体为x、y，且对于基因的的每一个维度i有xi，yi∈[ai，bi]，其中ai是维度i取值范围的下界，bi是取值范围的上界.子代基因的第i维从Fgene和Sgene两个区间中的一个遵循随机机制生成，Fgene和Sgene的表达式如式(6)、式(7)所示：

Fgene=[max{ai，Zbesti-d1i·β}，min{Zbesti+d1i·β，bi}](6)

Sgene=[max{ai，Zbesti-d2i·β}，min{Zbesti+d2i·β，bi}](7)

其中：d1i=|Zbesti-xi|，d2i=|Zbesti-yi|，β是自交率.

恢复系中的个体变异的概率为Pm，假设会发生基因变异的个体为x，根据高斯分布概率密度函数构造高斯变异函数，具体表达式如式(8)所示：

(8)

其中：Best是当前的最佳解，σ是高斯分布的标准偏差.由于三维刚体变换有6个参数，i是在区间[1，6]中随机生成的整数.

受到高斯突变的干扰，更新后的个体可以逃脱最优解的引力，从而跳出局部最优值.当前解的第i个位置的值如式(9)所示：

xi=xi·GMF(x，i)

(9)

3.4 具体应用

3.4.1 遗传个体编码

本研究采用实数编码的策略来代替二进制编码.对三维刚体变换的参数进行编码，组成遗传个体的基因.设种群的大小为N，编码后的遗传个体如式(10)所示：

Xi=(xi，1，xi，2，xi，3，xi，4，xi，5，xi，6)，i=1，2，3，…，N

(10)

其中xi，j代表三维刚体变换的第j个参数.三维刚体变换共有6个参数，分别是3个旋转参数和3个平移参数.参照文献[10，11]，旋转参数rx、ry、rz的取值范围是[-90，90]，平移参数tx、ty、tz的搜索空间为[-30，30].

3.4.2 适应度函数

本研究利用归一化互信息度量遗传个体的适应度.互信息可以看作是一个随机变量中包含的关于另一个随机变量的信息量，被广泛用于测量医学图像配准效果[3].设图像X，Y的联合分布概率为p(x，y)，p(x)和p(y)是他们的边缘分布概率.互信息I(x，y)是联合分布p(x，y)与乘积分布概率p(x)p(y)的相对熵，如式(11)所示：

=H(X)+H(Y)-H(X，Y)

(11)

其中：H(X)，H(Y)分别是X和Y的信息熵，H(X，Y)是他们的联合熵.

适应度函数的计算如式(12)所示，配准效果越好时归一化互信息值NMI越趋近于1，反之则越趋近于0.

(12)

3.4.3 算法流程

本文算法的执行流程如图3所示.

图3 改进后的遗传算法的流程图Fig.3 Flow chart of improved genetic algorithm

4 实 验

本节通过一系列实验来探究改进后的遗传算法(以下标记为GA*)在医学图像配准问题上的搜索能力以及早熟收敛的状况.

4.1 数据介绍和多分辨率策略

本研究选用来自Retrospective Image Registration Evaluation Project (RIRE)(1)https://www.insight-journal.org/rire/index.php提供的公开数据集中6号患者脑部的CT成像与MR-T1成像、MR-T2成像以及MR-PD成像进行配准.其中CT成像的大小为512×512×30voxels，MR成像的大小为256×256×26voxels.图4(a)、(d)、(g)分别是MR-PD、MR-T1、以及MR-T2成像，作为配准中的参考图像(Reference image)；图4 (b)、(e)、(h)是CT成像，作为配准中的浮动图像(Moving image).图4(c)、(f)、(i)显示了这3组图像初步对齐的效果.

图4 数据集可视化Fig.4 Visualization of datasets

本研究采用多分辨策略降低优化过程陷入局部最优解的概率.在进行配准之前首先对原图像进行下采样和高斯平滑处理，从而构建双层图像金字塔.金字塔的第1层是低分辨率的图像，第2层是未经处理的原始图像.本研究首先对第1层的低信息图像进行配准，记录最优解，并将最优解用于第2层配准参数的初始化.为了防止算法陷入局部最优解，优化算法需要在低分辨率图像上运行多次.

4.2 实验设计与参数设置

研究[6，7]表明，相较于传统的优化算法，进化计算和元启发式优化算法在图像配准问题上有着更佳的效果.本文选取5种元启发式优化算法与GA*进行对比，分别为改进后的散射搜索算法(SS*)[10]、珊瑚礁优化算法 (CRO)和带有附着层的珊瑚礁优化算法(CRO-SL)[11]、生物地理学优化算法(BBO)和基于精英学习策略改进的生物地理学优化算法(BBO-EL)[12].

本研究的配准流程如图5所示.

图5 配准执行流程Fig.5 Image registration process

所有的优化算法的参数设定严格与相关研究[10-12]保持一致.为了保证公平性，参照文献[11]，在图像金子塔第1层上的配准时间限制为180秒，在第2层的时间限制为280秒.第1层最优解集的大小由算法具体执行的效率决定.本研究将多样性控制阈值Diversity_Max设置为0.015.在基于灰度的图像配准实验中，文献[11]将所有优化算法的交叉概率和混合因子分别设置为0.5和0.3，所以GA*的杂交率α=0.5，自交率β=0.3.参照文献[9]，GA*的变异率Pm=0.08.

4.3 实验结果与分析

4.3.1 固定时间下的配准效果研究

本研究在给定的时间限制下对所给图像进行了配准.将MR-T1成像作为参考图像，CT成像作为浮动图像，采用GA*、SS*、CRO、CRO-SL、BBO以及BBO-EL这6种优化算法各自独立地进行了20次配准，并将归一化互信息值(NMI)作为衡量指标，实验结果如表1所示.

表1 多模态医学图像配准20次结果Table 1 Results of multi-model image registration

实验结果的描述性统计如表2所示.

表2 配准结果的描述性统计Table 2 Statistical description of registration results

从表1可以看出，相较于其他算法，GA*在所给多模医学图像配准的场景上总体表现最佳，在20次实验中有11次精度最高，占比55%.对表2的数据进行分析，从20次配准的NMI平均值来看，GA*达到了0.1839，而其余算法都在0.183以下.这说明GA*在所给数据集上的搜索能力较强，在给定的时间内寻优精度最高.除此之外，本研究参照所有算法20次配准的最差值对是否错误配准进行评判.去除CRO的极端值，其余算法最差值的平均值为0.1767，本研究将配准结果低于0.1767的视为误配.GA*的误配率为0%，和BBO-EL同时达到了最低，远小于CRO的50%.这说明GA*早熟收敛的概率相对较低.

表2的数据显示，GA*标准差为0.0017，明显低于大多数算法，仅高于BBO_EL的标准差0.0001，表现出了很强的鲁棒性.在20次迭代中，BBO_EL虽最高值低于GA*，但其最差值略高于GA*，解质量的稳定性略优.

4.3.2 固定代数下的配准效果研究

为了探究不同算法进化过程之间的差异，本研究进一步在固定的迭代次数下利用不同优化算法对所给数据集进行了配准实验.

4.3.1的实验已经表明，CRO-SL的配准效果要优于CRO且BBO-EL的配准效果要优于BBO，本实验主要针对SS*、CRO-SL、BBO-EL以及GA*展开进一步分析.实验严格按照图5所示的配准流程执行，在图像金子塔第1层上的配准时间限制为180秒，在第2层上取消了时间限制，迭代次数固定为100代，其余参数保持不变.本实验在所给数据集的3个多模配准场景上展开，4种算法随着迭代次数增加在原图像上的配准效果如图6所示.

图6 固定代数下的配准Fig.6 Image registration in fixed iteration

由于多分辨率策略的存在，不同的算法在图像金子塔第1层上表现的差异导致了在原图像上配准初始值的不同.从图6所示的3个场景可以看出，GA*的初始值均为最高，甚至高于CRO-SL迭代100次后的最终配准值.这说明多分辨率策略在很大程度上提高了本文算法的搜索能力，GA*在低分辨率图像上的表现最好，局部搜索能力较强，为高分辨率图像上的配准打下了很好的基础.

除此之外，GA*的全局收敛能力佳，可以在较少的迭代次数内接近或达到全局最优的配准效果，图6(a)和图6(c)中GA*在第10代的配准效果就已经超过了其他算法在迭代100次后所能达到的效果.最后，GA*的搜索不易陷入局部最优解，图6(a)中GA*在第10代已经达到了很好的配准效果，在第30代却仍然能够出现上升趋势，突破局部最优解.

表3给出了4种算法在3个场景下最终的配准参数、配准结果的NMI以及总耗时.从最终的配准结果来看，GA*在3个场景上的总体表现最好，时间性能略优.从总耗时来看，CRO-SL在最短的时间内完成了3个场景上的一百次迭代，其次是BBO-EL，GA*每一代的搜索的时间相对较长，有进一步提升的空间.

表3 4种算法的配准结果Table 3 Image registration results of four algorithms

以上实验表明，本研究提出的基于混沌映射和频度记忆的随机策略以及基于杂种优势的进化策略增强了遗传算法的搜索能力；基因多样性控制策略保证了种群的多样性，从而降低了遗传算法陷入局部最优解和早熟收敛的概率.相较于SS*、CRO、CRO-SL、BBO以及BBO-EL这5种优化算法，GA*在所给数据集上总体表现最佳.在给定时间限制的前提下GA*的性能有着较好的优势，而在固定进化代数的情况下GA*在早期低分辨率图像上的搜索效果较好，但在后期进化幅度不大.这说明相较于所对比的优化算法，GA*的主要优势在于可以在较短的时间内达到或接近全局最优解，适用于时间优先的应用场景.

总的来说，本文提出的基于多策略改进的遗传算法误配率低，搜索能力强，不易早熟收敛且有着很好的鲁棒性.

5 结 论

本文从随机机制、种群多样性以及遗传算子3个方面出发，提出多种策略来解决遗传算法在医学图像配准问题上局部搜索能力较差，容易早熟收敛等不足.并将多分辨率策略与改进后的遗传算法相结合，以NMI作为配准效果的测度函数，实现了医学图像的精确配准.实验结果表明，与SS*、CRO、CRO-SL、BBO以及BBO-EL相比，本文的方法可以较好地搜索到质量最高的解且不易早熟收敛，在很大程度上克服了标准遗传算法的不足，应用到医学图像配准问题上配准精度高、误配率低，有着良好的应用前景.未来我们将尝试把改进后的遗传算法与其他优化算法相结合，在不同的分辨率上发挥各自的优势，进一步提高医学图像配准的精度和效率.