量子遗传算法在口腔种植定位中的研究与应用

黄俊华，唐 平，陈松龄，江小平，梁英蓬

(1.广东工业大学自动化学院，广州 510006；2.中山大学附属第一医院口腔科，广州 510080；3.广东工业大学医院，广州 510006)

1 概述

口腔种植体通过外科手术植入人体缺牙牙骨中完成缺牙修复。传统的口腔种植外科手术术前设计都是基于X线平片(根尖片和全景片)，临床操作中常因设计估计不足而发生牙槽骨意外穿孔等严重并发症，而且往往给上部结构修复带来困难，导致种植修复的最终失败[1]。随着医学图像三维可视化等计算机辅助技术的发展，根据CT影像数据制作的手术导板增加了手术定位的准确性和操作的可控性。国内外分别推出VISIT导航系统和 SIM/plant等口腔种植软件[2-3]。文献[1]通过重建CT影像以人工交互的方式测量牙槽骨相关数据，但该方法工作步骤繁琐，容易造成因主观因素而导致的定位参数不准确。

结合量子计算和遗传算法的量子遗传算法有较好的全局搜索能力，已广泛地应用于解决组合优化和多属性决策规划问题。文献[4]提出用于解决背包问题的遗传量子算法(Genetic Quantum Algorithm, GQA)，种群个体的多样性是基于量子比特和量子态叠加特性表示的。文献[5]提出一种利用量子染色体相位比较来更新旋转角的量子遗传算法。文献[6]提出一种基于量子位 Bloch球面的量子遗传算法(Bloch Quantum Genetic Algorithm, BQGA)，该算法中种群的染色体有并行的3条基因链，优化性能优于普通量子进化算法。

现有的口腔种植系统智能程度不高，而且鲜有运用人工智能解决种植手术的学术报道。基于上述分析，为了提高种植手术的精度以及种植系统的智能处理能力，针对口腔医学螺旋CT数据量较大的特点，本文提出一种改进的 Bloch量子遗传算法(Improved Bloch Quantum Genetic Algorithm, IBQGA)，用于完成口腔种植体的智能设计。

2 量子遗传算法

2.1 基本量子遗传算法

基本量子遗传算法可以从以下方面进行说明：

(1)量子比特的表示

量子比特是量子计算的基本单位。基于Bloch球面的量子比特可以用φ和θ一对相位角度表示，其对应的向量表达式如下：

其中，|·>为量子表示符号。

本文采用基于Bloch球面的量子比特作为BQGA种群个体的基因位。设pi为基于Bloch球面的染色体，则其染色体编码为：

其中，i是进化种群中第i个染色体；φij、θij为第i个染色体第j号基因的一组相位角度；i=1,2,…,m，j=1,2,…,n，m是种群规模，n是量子位数。

(2)量子旋转门

在量子遗传算法中，对量子比特状态实施遗传操作是由量子旋转门来实现的，由酉矩阵与种群个体的量子比特编码运算实现。对于基于Bloch球面的量子比特，相应的量子旋转门是3×3的酉矩阵U：

其中，Δφ和Δθ都是量子旋转门的旋转角度，可以控制变换后的量子比特所处的状态；旋转角Δφ、Δθ的正负和大小控制着进化算法的收敛速度，如果幅值过大，会导致种群早熟；若幅值过小，会使收敛速度减慢。

文献[6]中的 BQGA算法充分利用量子算法的多态性，但实际进化的是染色体中的某一基因链，造成最终优化群体中只有小部分优化个体。文献[7]提出一种引入进化稳定策略的量子竞争决策算法，用以解决旅行商问题。该算法利用竞争力函数评价竞争者的竞争力，从而提高算法的全局优化能力。

2.2 改进算法

针对上述 BQGA存在的问题，结合本文应用，现提出下列3种改进方法：

(1)把种群分为不同的进化群体，即特征群体，使各特征群体独自进化。

考虑到不同的植入点对种植方案的差异，故以植入点为分类标准，把种群细分为不同的特征群体。每个个体在自己所属的特征群体中进化：即与同一特征群体中的其他个体比较，执行遗传操作和交叉操作。

(2)针对不同的特征群体及特征群体的进化程度，实施自适应调整进化步长的遗传操作。

对于量子遗传算法的遗传操作，文献[4-5]通过遍历查询表生成量子旋转门角度，此方法适合二进制编码的量子进化算法。虽能推广到一般的量子算法中，但涉及的多路判断影响算法的效率。在不同特征群体的基础上，相应进化特征群体中各个体的植入角度。本文采用如下自适应调整旋转角度的生成方法：

其中，φb和 θb为最优解中量子位概率幅所对应的相位；φc和θc是当前解中量子位概率幅所对应的相位；若A=0时，sgn(A)取正负号均可；若B=0时，sgn(B)取正负号均可。Δφ0和Δθ0同为种群进化的步长。fij是第j代种群某特征群体中个体i的适应值。fjmax和fjmin分别为当代此群体的最大优适应值和最小适应值。

(3)引入随机错位的交叉操作

在许多遗传算法的应用中容易发现，当种群进化到一定程度时，种群中各染色体相似或相同[8]，种群的多样性减弱，难以突破局部最优所带来的影响。而适当的交叉操作能保证遗传算法良好的搜索性能，有利于保持种群的多样性。考虑到经典遗传算法中，二进制编码的染色体执行交叉操作的特点是，分别随机地选择两父代的某一个基因位作为交叉基因位，父代之间的交叉基因位相互交换，生成2个新的子代染色体，从而子代染色体隐含父代染色体的特征信息。针对医学影像数据数据量大的特点，本文算法引入随机的错位交叉操作。随机的错位交叉操作的算法如下：

(1)随机生成一个0～1之间的数rand。

(2)由父代 V1、V2 分别按下式生成子代 V1’、V2’：

3 种植体定位参数优化算法设计

3.1 口腔种植定位问题的描述

设口腔种植定位问题有一个解集 SP={S1, S2,…,Sn}，n为解的总数，而每一个解Si都有一个评价种植方案 Si优劣的综合评价 criticism(Si)。每一个解Si={Positioni, Anglei, Li, Ri}，i=1, 2,…, n；每个种植方案有如下特征属性：植入点位置 Position，植入方向Angle，种植体长度L以及种植体半径R。解中的植入点位置是主特征属性，可根据各个解的植入点属性的差异，把解集细分不同的小集合SPj，SPj即为特征群体。口腔种植体定位设计就是要在解集SP中寻找一个这样的Si，Si的各项特征属性使其对应的综合评价criticism(Si)最好。

口腔种植手术中种植体定位，由种植体的半径、长度、植入点与方位角度决定。上述4个种植体信息构成具体的种植方案，是口腔种植体定位模型的参数。假设已知植入点的坐标(x0, y0, z0)，植入方向角度ρ和 ω，种植体的半径和长度分别为 r和 L；则可建立以植入点(x0, y0, z0)为原点的球面坐标系，分别以植入方向角ρ和ω为此球面坐标系的2个角度坐标，以种植体长度L为此球面坐标系的长度坐标，口腔种植体定位模型如图1所示。建立的球面坐标系是以种植植入点(x0, y0, z0)为原点的种植方案搜索空间。

图1 口腔种植体定位模型示意图

下列数学语言可以描述口腔种植体定位模型：

口腔的影像数据量庞大且包含噪声信息，因而不能直接作为搜索空间，还需要考虑定位模型和原数据空间的坐标转换。故需要提取缺牙区域的影像数据作为原数据空间。通过数字图像处理的方法去噪，简化区域数据；对提取的数据进行坐标变换，映射到搜索空间。因为口腔内牙列和牙槽骨的分布及形态具有一定的特征，把口腔牙列分为左上/下和右上/下 4个部分(分别记为I、IV、II和III区)。只需按照下式进行适当的坐标变换，口腔中缺牙区域的数据就能映射到相应的搜索空间，解决口腔种植体的定位问题：

其中，x、y、z为原数据空间中某点的3个直角坐标，sx、sy、sz为搜索空间下的坐标。Tn为属于n区的原数据空间映射到搜索空间的变换矩阵。

运用坐标变换和利用上面的定位模型，把种植方案的植入方向优化限制为2个植入角度的组合优化，且 2个角度的优化范围限制为 0°～90°，而不是 360°乃至空间中的搜索优化，能有效地减少进化算法搜索最优解的难度，有利于提高问题的优化效率。数据的提取与映射流程如图2所示。

图2 数据的提取与映射流程

3.2 量子染色体编码

一个染色体代表一个种植方案，分别由植入点、种植方向和种植体半径长度等参数决定。为了简化染色体结构，现考虑给定种植体半径和长度的情况，即每条染色体由植入点编码 1、植入点编码 2、角度编码1以及角度编码2这4个部分组成。每条染色体的量子比特均采用如式(1)表示的基于 Bloch球面的实数相位角度编码。

3.3 适应值函数

合理的种植体定位需要种植体外围有足够密度的牙骨支撑，以及植入方向符合生物力学、吻合度，此外，还需要考虑牙槽骨的骨密度、形态、宽度、高度和骨缺损等情况[9-10]。因此，本文算法的适应值函数 Fitness分为支撑部分 Fsup和方向部分 Fdir，Fsup评价此种植方案的种植体外围牙骨密度的大小；Fdir作为种植体对应的植入方向是否满足生物力学的评价。最优的种植方案有如下特点：适应值函数的支撑部分和方向部分的适应值总和最大，即max{Fsup+Fdir}。

Fsup作为评价种植体周围的牙骨灰度，假设植入坐标(x0, y0, z0)，角度ρ、ω，半径r和长度L已知，其数学表述如下：

其中，value(x, y, z)表示的是点(x, y, z)对应的骨灰度值；X、Y、Z为原数据空间的各维度的最大值。

种植方向的适应度Fdir，由2个定位方向角度来评价。牙长轴与种植体植入方向的夹角应在一定角度范围以内，为合理的种植方向，符合牙合力传导的规律[10]。定位方向的适应度可通过引入与植入方向和牙长轴的夹角相关的惩罚函数决定。如果植入方向和牙长轴的夹角较大，则使适应度值为0；若夹角在口腔医学允许范围之内，可增大对应的适应度值。

本文算法步骤如下：

Step1 随机初始化种群，并建立一个用以记录种植点、植入方向和最优适应值的空表 list。初始化迭代截止代数。

Step2 计算种群中各个体的适应值。然后根据植入点的不同划分特征群体，记入 list。其中，各特征群体需记录最优及最小适应值。

Step3 根据list表中当前特征群体的最优适应值以及最小适应值，分别按式(2)～式(4)更新群体中个体，直至当前特征群体更新完毕。

Step4 重复Step3，直至全部群体被更新为止。

Step5 更新当前代数，并返回 Step2，直到满足收敛条件为止。

4 实验结果与分析

为了证明本文算法在解决种植体定位问题具有优于经典优化算法的特性，下面将本文算法与 GA、BQGA进行比较。生成一个模拟CT医学影像数据的三维数组，作为测试种植体定位算法的实验数据。实验数据是一个随机生成的服从正态分布的大小为10×10×10的三维数组，仅设置 3个对角相邻，且数组标号递增的数组元素值为10，已知最优的坐标及对应的最优角度(π/4)和理论最优值(37)。下面利用GA、BQGA与本文算法测试上述测试数据，得出各优化算法的综合比较结果如表1所示，其中，qq表示平均角度；pp表示角度平均相对误差；kk表示平均优化值。

表1 各优化算法的综合比较结果

3种算法种群规模均为 200，截止进化代数均为200。GA中的交叉概率为 0.8，变异概率为 0.3，而BQGA与本文算法的交叉概率均为 0.6，变异概率均为 0.1。进化算法的优化值相对误差比较如图 3所示。

图3 进化算法的优化值相对误差比较

从表1的实验数据中，可以得出：GA的优化效率不如 BQGA和本文算法高。在 10次实验中，GA只有2次优化接近最优方案，即植入坐标优化为最优的植入坐标，总体优化值的平均相对误差是0.416 6；而 BQGA和本文算法均能收敛且优化值的平均相对误差均优于GA。在本文实验的GA中，种群中的坐标以及角度的优化信息都依赖于交叉以及变异算子。即使增大种群数量或进化代数，也难以保证GA能搜索到理论的最优坐标。而角度优化仅在坐标得以优化下讨论才有意义，因此，植入点坐标的优化程度决定了最终优化解的优劣。在种群进化过程中，如果种群中坐标及植入角度得到充分的优化，最终得到的优化解就越接近问题的最优解。与 BQGA相比，由于本文算法考虑到种群个体间的差异性，本文算法在角度平均相对误差及优化值的平均相对误差上均优于BQGA算法。

针对中山大学第一附属医院提供的 CT影像资料，结合本文研究的种植定体位模型，将本文算法应用于种植体方案设计优化上，完成种植体的智能生成，最终优化结果的矢状切面图如图4所示，其中，牙骨中的直线部分为种植体。

图4 本文算法优化种植体方案的矢状切面图

由图4可见，本文算法能正确生成缺牙区域的种植体种植方案。但就其植入方向等指标的合理性，仍有待提高与改进。造成这样的原因可能是图像处理的方法及其参数的选取设置不当造成的，也可能是量子遗传算法适应值函数的设置问题。

5 结束语

本文提出一种改进的量子遗传算法。利用基于Bloch球面的量子基因编码种群多样性，引入交叉操作，使该算法的全局搜索能力和种群的多样性得以保证。实验结果表明，该算法比GA和BQGA有较好的搜索能力和收敛性能，能解决种植体的定位问题。今后的研究内容主要有：如何合理设置进化算法的参数与约束，如种植体长度半径以及适应函数的定型等。

[1]康 璐, 黄远亮, 顾立栩, 等.计算机辅助设计软件在口腔种植外科的应用研究[J].口腔颌面外科杂志, 2008,18(4): 265-268.

[2]吴 婷, 廖文和, 戴 宁, 等.口腔种植导板计算机辅助技术研究[J].生物医学工程学杂志, 2011, 28(1): 1-6.

[3]王培志, 夏 露, 陈 宁, 等.种植导航模板的计算机辅助设计和制造[J].中国口腔种植学杂志, 2010, 15(3):128-130, 133.

[4]Han Kuk-Hyun, Kim Jong-Hwan.Genetic Quantum Algorithm and Its Application to Combinatorial Optimization Problems[C]//Proc.of IEEE Conference on Evolutionary Computation.[S.l.]: IEEE Press, 2000.

[5]李士勇.一种基于相位比较的量子遗传算法[J].系统工程与电子技术, 2010, 32(10): 2219-2222.

[6]李士勇, 李盼池.量子计算与量子优化算法[M].哈尔滨:哈尔滨工业大学出版社, 2009.

[7]刘 勇, 马 良, 宁爱兵.量子竞争决策算法及其在旅行商问题中的应用[J].计算机应用研究, 2010, 27(2):586-589.

[8]陈国龙, 陈火旺, 郭文忠, 等.基于随机错位算术交叉的遗传算法及其应用[J].模式识别与人工智能, 2004,17(2): 250-255.

[9]周 苗, 陈松龄, 陈建灵, 等.螺旋CT在牙种植术前评估和设计中的应用[J].实用医技杂志, 2005, 12(6):1562-1563.

[10]孙婷婷, 潘 瑾, 崔明军, 等.基于CT影像的牙种植模板相关的颌骨解剖学研究[J].口腔颌面外科杂志, 2007,17(3): 252-255.