基于改進粒子群優化算法的牙齒正畸路徑規劃方法

徐曉強，秦品樂，曾建朝

（中北大學大數據學院，太原 030051）

（*通信作者電子郵箱zjc@nuc.edu.cn）

0 引言

正畸（Orthodontics）又稱牙齒矯正，是指利用各種牙齒矯治器對患者牙齒的位置進行重新排列，使其恢復到與正常牙齒形態相對一致的過程［1］。牙齒矯治器包括傳統有托槽矯治器、無托槽隱形矯治器等。近年來，隨著無托槽隱形矯治器相較于前者具有更舒適、美觀、節約矯治時間、容易清潔等特點［2］，并且患者通過無托槽隱形矯治系統可以提前看到矯治完畢后的牙齒具體排列情況。因此，越來越多的患者選擇無托槽隱形矯治器進行牙齒矯正，這使得運用計算機技術對按數學模型排列的虛擬牙齒正畸技術成為研究的熱點。復雜的牙齒正畸路徑規劃作為上述研究熱點的一個重要步驟，既要保證矯治路徑具有可行性，同時也要保證矯治結果具有優化性。專業地說，牙齒正畸路徑規劃是保證患者的全部牙齒矯治過程在符合正畸學的前提下，使多顆牙齒在經過數個矯治階段后達到最終理想位置，并且在該過程中相鄰牙齒間不發生碰撞，因此需要尋找一條盡可能減少患者矯治時間、移動距離、旋轉角度等因素的最優路徑作為最終牙齒正畸移動路徑。文獻［3］提取了人工牙三維三角網格數據的特征點建立了牙齒選擇系統，通過計算機數控技術完成了全口義齒的虛擬排牙方法。文獻［4］提出單顆牙齒自動分割方法優化了傳統方法，通過搜索分割路徑、圖像形態，擬合牙弓線，匹配最優路徑確定牙縫邊界路徑，細化封閉分割輪廓完成整個牙齒的分割。此方法對于嚴重畸形的牙齒具有較好的正畸效果。文獻［5］提出了一種基于數字圖像與計算機斷層掃描模擬的虛擬近似真實正畸牙齒移動，并通過給牙齒特定位置施加力來減少平移及旋轉次數。文獻［6］根據壓力和張力與特定的信號傳導相關因子，建立了局部梯度來重塑牙周，以更快地移動牙齒，減少治療時間和依賴時間的不良后果。

隨著對遺傳算法、神經網絡算法等啟發式算法的深入研究，發現上述算法大多存在計算量大、操作復雜、算法運行時間過長等問題，越來越多的研究者提出運用粒子群優化（Particle Swarm Optimization，PSO）算法來解決路徑規劃的相關問題。PSO 算法是Eberhart等［7］通過模擬鳥類行為規律，提出的一種全局優化算法。PSO 算法因其結構簡單、效率較高等優點已經廣泛應用于各類優化問題［8-9］，近年來對改進粒子群算法在路徑規劃領域的研究陸續被提出。文獻［10］提出了斥力場下粒子群，利用柵格法進行初步規劃以及根據不同的障礙物所構建數學模型、適應度函數，從而得到了移動機器人最優路徑。文獻［11］在光柵化的基礎上構建了包括隨機地表模型、山地地形模型的數學模型，并分析了路徑威脅因素，引入遺傳算法的思想改進了粒子群算法最終得到了低威脅的干擾機路徑。文獻［12］提出了一種基于高斯參數更新規則和粒子重新初始化的改進粒子群算法，從縮短路徑長度及提高路徑光滑度兩方面優化了路徑規劃問題。文獻［13］通過個體粒子進化速度與群體離散度來調整粒子群的慣性權重，以此避免其過早收斂，并利用自然選擇方法改進量子行為粒子群的優化算法，將這兩種優化的算法引入移動機器人路徑規劃，并通過實驗證明方法的可行性。文獻［14］以ER3A-C60 機器人作為研究對象，建立了運動學和動力學模型，用參數合適的PSO算法求得了最小綜合誤差得到了一條最優軌跡。

受上述改進算法的啟發，針對虛擬口腔正畸治療系統中牙齒移動路徑規劃問題，本文提出了一種基于正態分布的簡化均值粒子群的牙齒正畸路徑規劃方法。首先建立了單顆牙齒及整體牙齒的數學模型，并將牙齒正畸路徑規劃問題轉化為帶約束的優化問題；其次，引入正態分布及均值粒子群，結合簡化粒子群的位置項，提出了一種基于正態分布的簡化均值粒子群優化（Simplified Mean Particle Swarm Optimization based on Normal distribution，NSMPSO）算法；最后，構造了包含平移路徑長度、旋轉角度、碰撞檢測以及牙齒在單階段的移動量、旋轉量的高安全性的適應度函數，實現了牙齒正畸移動路徑的規劃。在三大基準測試函數上將NSMPSO 算法與其他三種粒子群算法進行對比，結果表明，本文改進的算法在收斂速度和收斂精度上均有明顯的提高，通過Matlab 的仿真實驗表明了利用該數學模型和改進算法求得的最優路徑安全可靠，可以為醫生輔助診斷提供有效的幫助。

1 建立數學模型

1.1 單顆牙齒及整體牙齒的數學模型

本文以分割后的牙齒為研究對象。首先，采用有向包圍盒（Oriented Bounding Box，OBB）對分割后的單顆牙齒建立包圍盒，以包圍盒的中心點作為牙齒的特征點建立局部坐標系，X、Y、Z軸的方向根據單顆牙齒包圍盒的長寬高來確定。圖1為尖牙、側切牙、后磨牙在三維數學模型下單顆牙齒包圍盒的局部坐標系。然后，以全部牙齒為目標，選取兩顆后磨牙、兩顆尖牙的特征點所確定的平面作為基準面來建立牙齒整體的世界坐標系。圖2 為整體牙齒在三維數學模型下的世界坐標系。最后，用每顆牙齒的包圍盒及中心點代替牙齒進行整個矯正過程的展示，可以更直觀、形象地了解每個階段內的牙齒位置的變化情況。在某一矯正階段內，由于全部牙齒的移動和旋轉是同時進行的，則可將該階段的牙齒移動路徑離散化表示為一系列路徑點的連線。

圖1 單顆牙齒的三維坐標系Fig.1 Three-dimensional coordinate system of single tooth

圖2 整體牙齒數學模型Fig.2 Mathematical model of whole teeth

1.2 牙齒移動路徑規劃的數學建模

將n顆牙齒的移動路徑規劃分為m個階段，m由每個階段可移動、旋轉的最大值及醫生經驗共同決定。牙齒在矯治過程中，在滿足單階段內患者可以承受的牙齒最大移動量、旋轉量并且無碰撞發生的前提下，確定每個階段的所有牙齒的位置信息，最終規劃出一條牙齒從初始位置到最終位置的最短路徑，即為可執行的牙齒移動最優路徑。因此可將牙齒正畸路徑規劃問題等價為帶約束的優化問題：其中：n為牙齒個數，m為矯治階段數。(xij，yij，zij)表示第i顆牙齒在第j個矯治階段結束后的位置坐標，目標函數之一f1表示牙齒平移路徑長度，是通過計算空間中兩點間的歐氏距離構建的。計算患者所有牙齒在全部矯治階段內的移動量之和，值越小表示牙齒在矯治過程中需要移動的量越小，從而減少了患者的矯治時間。目標函數f2表示牙齒旋轉角度，(αij，βij，γij)為第i顆牙齒第j個階段時包圍盒的三邊與世界坐標系的三軸之間的夾角。單顆牙齒包圍盒的局部坐標軸X1、Y1、Z1在全部牙齒建立的世界坐標軸X2、Y2、Z2的投影分別為X1X2、Y1Y2和Z1Z2，組成的三個夾角即為α、β、γ，計算患者全部牙齒在整個矯治階段內的旋轉量之和即為目標函數二。根據正畸學及醫生臨床經驗可知，每顆牙齒在某個矯治階段內移動范圍不能大于0.5 mm、旋轉角度不能超過3°。考慮到每顆牙齒在規定范圍內進行移動或者旋轉之后，可能出現某顆牙齒與左右相鄰的牙齒發生碰撞擠壓導致矯治失敗的情況，因此在每個矯治階段檢測牙齒間是否發生碰撞至關重要。綜上，用g1、g2、g3表示在對目標函數f1、f2進行優化時的約束條件。約束條件g1是限制每顆牙齒在前后兩個階段的移動量小于每個階段單顆牙齒能夠承受的最大移動量；約束條件g2是限制每顆牙齒在前后兩個階段的旋轉量小于人類能夠承受的每個階段的最大旋轉量；約束條件g3是檢測矯治過程中相鄰牙齒是否碰撞，式中ci表示判斷第i顆牙齒與其左右相鄰的牙齒是否發生碰撞。如果兩顆牙齒沒有發生碰撞，則ci等于0；否則ci不為0。因此使得約束條件g3中的fc等于0，則滿足牙齒在矯治過程中無碰撞發生。

2 用NSMPSO算法實現牙齒正畸路徑規劃

2.1 粒子群算法

在D維搜索空間中，隨機初始化z個粒子，這些粒子可視為優化問題的可行解，并根據預先設定的適應度函數來評價粒子所在位置的好壞。粒子i的當前位置為Xi=(xi1，xi2，…，xiD)，當前速度為Vi=(vi1，vi2，…，viD)。在每次迭代中，pi=(pi1，pi2，…，pid，…，piD)表示粒子迄今為止經歷的最優位置，pg=(pg1，pg2，…，pgd，…，pgD)表示整個種群迄今為止經歷的最優位置。粒子的速度和位置的更新公式如下：

其中：下標i=1，2，…，z，下標d=1，2，…，D，t是當前迭代次數，w是慣性權重。c1和c2是學習因子，r1、r2～U（0，1）的隨機數。為了使粒子在搜索空間內充分搜索，在迭代過程中限定粒子速度向量上下限為［vmin，vmax］，位置向量上下限為［xmin，xmax］。

2.2 基于正態分布的簡化均值粒子群算法

在D維搜索空間中，粒子均按照一定的速度和方向運動，因此忽略粒子在運動過程中的體積與質量，可將其看作一個點。在連續解空間過程中，粒子逐漸向全局最優解附近靠攏，可將多維可行解空間看作服從多元正態分布。在真實的迭代過程中，粒子在可行解空間的分布不完全意義上服從正態分布，相比之下仍存在一定的偏離，但鑒于這一部分的偏差的存在對整體實驗結果的影響不大，可忽略不計，并且絕大多數粒子的分布是服從正態分布的，所以引入正態分布作為真實粒子在搜索空間分布的近似值是合理可行的。因此，本文引入正態分布和均值粒子群［15］的思想，結合簡化粒子群［16］中的位置項，提出了一種基于正態分布的簡化均值粒子群（NSMPSO）算法。改進后粒子的位置公式如式（4）所示：

其中w為慣性權重。式（4）將簡化粒子群的位置項擴展為四項：第一項表示向粒子當前位置學習項。第二、三項中用粒子個體最優位置、種群全局最優位置的線性組合取代粒子的個體最優位置和全局最優位置，用慣性權重調整影響大小。第四項為正態分布項，其中，μ表示粒子當前位置、粒子個體最優位置及全局最優位置三者的平均值，根據式（5）計算；σ表示粒子當前位置、粒子個體最優位置及全局最優位置三者的標準差，根據式（6）計算；l1和l2是［0，1］范圍內均勻產生的隨機數，K是由指數函數和三角函數組成的常數項，用來控制標準差的影響程度，根據式（7）計算。為了避免出現標準差為0的情況，規定當標準差閾值t小于1E-05時，標準偏差的值取平均值的一半。綜上所述，本文提出的NSMPSO 算法增強了種群的多樣性，加入正態分布項進一步改善了PSO 算法在解決高維問題時求解精度不高、容易陷入局部最優等缺點。

2.3 構造適應度函數

構造一個合理的、安全性高的適應度函數不僅可以作為衡量優化得出的牙齒正畸移動路徑優劣程度的標準，同時也是檢測本文改進的粒子群算法在收斂速度及收斂精度方面是否有所提高的判斷準則。根據牙齒正畸移動路徑規劃問題的特點，本文對它的評價取決于5 個方面，分別是平移路徑長度、旋轉角度、是否碰撞以及牙齒在單階段的移動、旋轉量。

平移路徑長度評價函數對應目標函數f1，旋轉角度評價函數對應目標函數f2，用式（1）中的g1、g2表示對牙齒在矯治過程中某一階段移動距離和旋轉角度的限制約束。本文對約束條件g3通過添加常數項進一步改進，函數f3作為評價是否發生碰撞的標準，定義如下：

其中，S表示給定的一個較大常數項，其作用是當相鄰的牙齒發生碰撞時，較大的常數S可以使得整個f3的值明顯增大。如果兩顆牙齒沒有發生碰撞，則f3=0。

綜上，結合f1、f2、f3、g1、g2定義本文算法的適應度函數為：

其中，i=1，2，…，z；δ、η、λ、μ、θ分別為路徑長度、旋轉角度、是否碰撞以及牙齒在單階段內的移動、旋轉量在適應度函數中的權重，且為大于0 的實數。式（9）中f1和f2是基于目標函數1 和2 構建的，設置其權重略大于f3、g1、g2的權重。分析適應度函數Fi可知，函數值越小說明得到的牙齒正畸路徑越優。

2.4 NSMPSO算法在路徑規劃中的實現

2.4.1 解的編碼問題

PSO 算法采用實數編碼，一個粒子代表著所有牙齒在全部階段內的矯治結果，也就是說每個粒子對應著一個可行解。由于每個粒子的位置由每顆牙齒在每個矯治階段中平移所確定的坐標(xij，yij，zij)和旋轉所確定的(αij，βij，γij)組成。牙齒個數為n，矯治階段數為m，因此粒子的位置是6×n×m維變量。這樣，粒子的編碼結構如下所示：

2.4.2 算法實現的具體步驟

步驟1 獲取患者牙齒數據，提取牙齒特征點作為初始位置坐標，本文采用β函數來擬合理想牙弓曲線，從而確定牙齒矯正結束后的最終位置坐標［17］，公式如下：

其中：W是患者兩個磨牙特征點之間的距離，H是患者中切牙特征點到兩顆磨牙特征點之間連線的最短距離，e是將患者尖牙的特征點坐標代入式中所確定。

步驟2 按照2.4.1節對粒子進行編碼，初始化種群并設定相關參數，其中包括患者牙齒個數n，需要矯治階段數m，種群大小z，每個粒子維數D，最大迭代次數Tmax。初始化粒子的位置，公式如下：

其中，r為［0，1］區間的隨機數。

步驟3 根據式（9）計算每個粒子的適應度值，并更新粒子的個體最優值和全局最優值。

步驟4 根據式（5）計算粒子當前位置、粒子個體最優位置及全局最優位置的平均值，根據式（6）計算三者的標準差，根據式（7）計算常數項K。添加正態分布項，對種群中的每個粒子根據式（4）進行位置更新。

步驟5 判斷算法迭代次數是否達到所設定的最大迭代次數，若是，停止迭代，輸出牙齒正畸移動最優路徑；否則，返回步驟3）繼續執行。

算法流程如圖3所示。

圖3 NSMPSO算法流程Fig.3 Flowchart of NSMPSO algorithm

3 實驗環境及分析

3.1 測試函數

為了驗證本文提出的NSMPSO 算法在收斂精度及收斂速度方面上的提升，將NSMPSO 算法與基本粒子群優化（PSO）、均值粒子群優化（Mean Particle Swarm Optimization，MPSO）和動態調整慣性權重的簡化均值粒子群優化（Simplified Mean Particle Swarm Optimization with Dynamic adjustment of inertia weight，DSMPSO）算法進行對比。本實驗選取三大基準函數作為測試函數，其中：Sphere 函數是單峰函數，主要用于檢驗算法的收斂速度和求解精度；Ackley 函數和Griewank 函數是多峰函數，主要用于檢驗算法的全局搜索能力。測試函數描述如下。

1）Sphere測試函數：

其中，搜索區間為［-100，100］，當f1(x)=0時，Sphere 測試函數取得全局最優值。

2）Ackley測試函數：

其中，搜索區間為［-32，32］，當f2(x)=0時，Ackley 測試函數取得全局最優值。

3）Griewank測試函數：

其中，搜索區間為［-600，600］，當f3(x)=0時，Griewank 測試函數取得全局最優值。

3.2 參數設置

在實驗中，測試函數的維數D=30，最大迭代次數genmax=150，種群粒子數m=30。對于PSO 算法和MPSO 算法，w=0.8、c1=c2=2。對于DSMPSO算法，c1=c2=2，wmax=0.9、wmin=0.4、f1(x)==0.1、a=1、b=2。對于NSMPSO 算法，w取［0.4，0.9］區間的隨機數，c1=c2=1.494。

3.3 實驗結果分析

3.3.1 收斂性分析

本實驗為了驗證NSMPSO 算法的收斂性，分別與PSO、MPSO、DSMPSO 算法在三大測試函數上進行對比分析。四種算法迭代150次的收斂效果如圖4所示。

從圖4中可以看出，對于單峰Sphere 函數，基本PSO 算法由于全局搜索能力較弱，陷入局部最優導致算法過早停滯，在迭代150 次后未收斂，而MPSO 算法收斂速度較慢，DSMPSO、NSMPSO 算法均能較快找到全局最優解，相比之下，NSMPSO算法適應值小于其他三種算法，求解精度相對較高。對于多峰Ackley 函數，基本PSO 算法由于難以處理復雜的多峰函數導致算法停滯，MPSO 算法由于跳出局部最優能力較弱，經過140 次迭代后適應值開始趨于穩定，相比其他兩種改進算法，收斂速度過慢，而DSMPSO 和NSMPSO 算法均能在50 次內到達全局最優值或者達到最優值附近可接受的范圍內，分別在40 次、25 次左右開始趨于收斂。對于多峰Griewank 函數，基本PSO 算法在迭代150 次過后，整體適應值變化不大，沒有達到收斂精度。MPSO 算法在解決復雜多峰函數時穩定性不強，DSMPSO 算法以及NSMPSO 算法分別在迭代20 次、15 次左右后適應值穩定收斂，從圖4（c）可知，NSMPSO 算法找到全局最優解所需迭代次數最少，收斂速度最快。

綜上所述，本文提出的基于正態分布的簡化均值粒子群算法（NSMPSO）無論對于單峰函數還是對于多峰函數來說，在尋找最優解過程中均能有效地避免種群陷入局部最優解，并且能夠在較少的迭代次數內達到穩定收斂，比其他PSO 算法具有更快的收斂速度，在求解精度方面也較為出色，全局搜索能力更佳。

圖4 三個函數收斂曲線對比Fig.4 Convergence curve comparison of three functions

3.3.2 有效性分析

為了進一步驗證本文提出的NSMPSO 算法的有效性，選取某患者下牙頜14 顆牙齒作為矯治對象，選取牙齒特征點作為初始位置，根據理想牙弓曲線確定牙齒矯治結束時的位置。利用NSMPSO 算法在約束條件g1、g2、g3的約束下對目標函數f1和f2進行優化，從而得到每顆牙齒在每個階段時的坐標位置，然后在Matlab上進行仿真實驗。

首先，牙齒在矯治過程中垂直方向上的移動量較小，因此本實驗忽略牙齒在垂直方向上的移動，將原本為一條三維空間的牙弓曲線投影到由患者兩顆磨牙、兩顆尖牙以及中切牙所確定的牙頜平面上，從而簡化成一條二維平面上的曲線。其次，患者在矯治過程中所有牙齒是同時進行移動、旋轉的，并且需要確保過程中牙齒間不能發生碰撞、擠壓，否則容易導致矯治失敗，這也是將碰撞檢測作為牙齒正畸路徑規劃的重要原因之一。本實驗利用OBB 進行碰撞檢測，在二維平面上用包圍框代替牙齒可以更直觀、形象地表示牙齒移動的可行性，確保牙齒正畸的每一階段都是合理移動的。最后，以理想牙弓曲線作為牙齒正畸的標準，牙齒正畸路徑規劃的結果如圖5所示。

圖5 牙齒初始位置與牙齒矯正結束位置對比Fig.5 Comparison of initial position of teeth and position of teeth after orthodontics

圖5 表示牙齒初始位置與牙齒矯正結束位置對比，其中實線包圍框表示牙齒初始位置，虛線包圍框表示牙齒矯正結束位置，包圍框的中心點即為每顆牙齒的特征點。為了方便敘述，從左到右依次對牙齒進行編號為1～14，可以看出該患者牙齒排列不齊，里出外進，在一塊狹小的區域內多個牙齒相互重疊遮擋。牙齒正畸路徑規劃過程如圖6 所示，其中，牙齒從初始位置到矯正結束一共需要7 步，每一階段矯正結束時的位置在圖6中表示，由于1 號和14 號第一磨牙、2 號和13 號第二磨牙相對移動、旋轉較小，因此在矯治第三階段時第一磨牙、第二磨牙已經達到理想位置，其他牙齒在經過七個階段的矯治后，牙齒之間排列整齊，沒有縫隙，10 號右側尖牙和4 號前磨牙、3 號前磨牙沒有之前被擠出的跡象，并且所優化出來的結果在牙齒矯治過程中，各個相鄰的包圍框沒有發生重疊部分，意味著牙齒之間沒有發生碰撞，說明本文提出的改進粒子群算法NSMPSO應用到牙齒正畸路徑規劃是合理可行的。

圖6 牙齒正畸路徑規劃結果Fig.6 Orthodontic path planning results

圖7 表示四種算法迭代150 次的牙齒正畸路徑規劃收斂效果對比。可以看出，采用PSO 和MPSO 算法求解牙齒正畸移動最優路徑時效果不佳，很難找到全局最優解。與DSMPSO 相比，本文提出的NSMPSO 算法有更快的牙齒正畸路徑規劃收斂速度，在迭代45 次左右算法開始收斂，表示NSMPSO 算法可以規劃出一條有效的牙齒正畸移動最短無碰撞路徑。

圖7 牙齒正畸路徑規劃收斂圖Fig.7 Orthodontic path planning convergence diagram

4 結語

針對虛擬口腔正畸治療系統中牙齒移動路徑規劃問題，本文提出了一種基于NSMPSO 算法的牙齒正畸路徑規劃方法。首先根據牙齒運動特性，將問題轉化成帶約束的多目標優化問題；其次，本文在簡化粒子群算法的基礎上，引入正態分布及均值粒子群的思想，提出了基于正態分布的簡化均值粒子群算法（NSMPSO）；最后，構造了包含平移路徑長度、旋轉角度、碰撞檢測以及牙齒在單階段的移動量、旋轉量的適應度函數，得到牙齒從開始位置到最終位置的最優路徑解。在三大測試函數上進行對比實驗，驗證了提出的NSMPSO 算法與其他粒子群算法相比具有最高的收斂速度及最高的求解精度。同時在Matlab中進行仿真實驗表明，利用該數學模型和改進算法求得的最優路徑安全可靠，可以為醫生輔助診斷提供有效的幫助，具有一定的實用價值。