義齒應力分布數據是評價義齒工藝參數合理性的重要判據，而源于對頜牙的咬合力是計算該數據的關鍵載荷條件。為了盡可能地提高義齒應力分布數據的準確性，本文以正中咬合下的下頜第一磨牙全瓷冠為研究對象，采用一種能夠反映咬合實際狀況的咬合力加載方式。首先對標準牙頜模型中下頜第一磨牙咬合面形貌進行光柵掃描與三維重建，同時根據預備體制備工藝進行結合面的曲面造型，并運算得到義齒與預備體幾何模型；其次，借助咬合動力學分析軟件（OFA）函數庫開發參數化的咬合分析程序，并使用遺傳算法對正中頜位進行優化；最后，基于咬合優化結果設計網格模型優化工況，同時以尖窩接觸特征設計參考工況，將以上數據均導入有限元分析軟件 Abaqus 中進行應力分布的對比分析。本文研究結果表明，遺傳算法適用于咬合優化問題的解決；相較于參考工況，咬合優化后的義齒具有更小的應力最大值和更均勻的應力分布特征。基于以上研究結果證明，本文所提出的正中咬合優化方法能夠進一步提升義齒在有限元模型中的受力準確性，也為人體其他關節的受力分析提供了新思路。
引用本文： 秦文龍, 叢明, 任翔, 劉冬. 基于正中咬合優化的磨牙全瓷冠應力分析. 生物醫學工程學雜志, 2020, 37(5): 802-808, 817. doi: 10.7507/1001-5515.201902013 復制
引言
正中咬合（centric occlusion，CO）是牙列間最均勻、最廣泛、最緊密與最穩定的接觸，且伴隨著咀嚼循環中最大咬合力的發生，其作為關鍵咬合狀態影響著義齒的力學表現[1-2]。目前，有限元分析法已經成為評價義齒正中咬合狀態力學表現的重要途徑，可用于獲取感興趣部位的應力、應變和受力等數據[3]。考慮到載荷是影響有限元分析結果的基礎性因素[4]，因此盡可能建立與正中咬合狀態相吻合的載荷施加方式是提高義齒有限元建模準確性的關鍵，這就要求研究設計的技術方案能夠同時解決好正中咬合重建與咬合力設計這兩項內容。
由于不同的正中頜位（centric position，CP）將產生不同的咬合分析結果[5]，因此對于不確定的、初始的正中頜位有進行咬合重建的必要。并且，正中咬合的重建質量將在很大程度上決定咬合力設計的質量。常見的方法是，對人體或體外正確頜位關系采用光學[6-7]、探針[8]掃描記錄的方式在虛擬模型中建立正中咬合關系。目前，對咬合模型進行掃描已經成為牙頜模型掃描儀的標準操作步驟[9]，但該方法不適用于咬合關系不明的場景（例如：牙頜缺失）[10]。當無法直接獲取正中咬合狀態時，可以采用計算機輔助的頜位搜尋技術輔助研究。在咬合動力學分析軟件咬合指紋分析系統（occlusal fingerprint analyser，OFA）（ZiLoX IT Inc.，德國）中，具備基于接觸面積最大的頜位遍歷功能，但卻存在計算步數過多、解集質量不佳的問題[11]。同時，未見將正中咬合重建視為優化問題進行求解的報道。
咬合力設計，其目的是解決如何將咬合分析結果（例如：虛擬模型的接觸點云、實體模型的試紙壓痕等）設計映射到有限元網格模型。目前，咬合力加載方式包括集中力與均勻分布力兩種。其中，集中力的載荷施加方法，通常是將集中載荷以某個方向角施加在功能牙尖或牙脊上，因為凸起的形貌特征更容易成為碰撞的先發位置[12]。另外，根據尖窩接觸的特征，頰窩、近中窩和中央窩也常作為載荷施加點位[13-14]。然而，這種集中力的施加方式忽略了正中咬合復雜的接觸點分布特征，局部容易產生失真的應力集中現象。越來越多的研究使用均勻分布力作為載荷施加方式[15]。咬合力被施加在理想的中央窩區域或頰尖舌斜面，同時力的方向與牙齒長軸呈 45° 夾角。然而，此時的接觸點位置及接觸面邊界并非源于咬合分析結果，更多是經驗性的學知識，因此往往是規則的、理想的，與實際的接觸細節存在差異。為了克服這一不足，研究人員開始根據咬合動力學分析軟件或咬合試紙分析所得非規則的與實際吻合的接觸面，設計均勻分布式咬合力[16-19]，但是仍存在接觸位置映射至網格模型失準的情況。
基于以上原因，目前的正中咬合力加載方式同時存在正中咬合重建質量有限與咬合力加載位置失準的不足。為了解決這一問題，本文針對下頜第一磨牙全瓷冠有限元模型，創新性地提出一種基于接觸面積最大的正中咬合優化模型，并根據咬合優化結果設計基于網格模型的咬合力加載方式。通過本文研究，期望可以為義齒有限元模型提供新的咬合力加載方式，實現應力計算準確度的進一步提升，另一方面也可以為人體其他關節的受力分析提供新思路。
1 下頜第一磨牙全瓷冠光柵掃描與三維重建
1.1 光柵掃描與三維重建
牙頜模型是根據不同牙科教學或實驗目的仿照不同類型人體牙列的工業模型，其中標準牙頜模型是具有正常咬合面形貌的仿真牙列。本文選取大連醫科大學口腔醫學院提供的上下頜右側第一磨牙（D16FE-500HGSF-QF，Nissin. Inc，日本）作為人體磨牙的代替物進行掃描。具體掃描步驟如下：① 對掃描樣品進行熒光粉噴涂處理；② 使用光柵掃描設備（AutoScan DS100，先臨 3D，中國）分別對單顆牙進行掃描，記錄牙齒形貌，其掃描精度 0.015 mm。借助逆向工程軟件 Geomagic Studio（Geomagic Inc.，美國）對掃描獲取的牙齒表面點云模型進行多邊形破損修復和精確造面，實現咬合面、鄰面、頰面和舌面的三維重建，如圖 1 所示。
圖1
下頜第一磨牙全瓷冠三維模型
Figure1.
Three dimensional model of the mandibular first molar with all-ceramic crown prosthesis
1.2 預備體建模
根據全瓷冠基牙制備流程[20]，右側下頜第一磨牙預備工藝可分為肩部預備、鄰面和頰舌面預備、咬合面預備、斜面預備和光順打磨等 5 個步驟，涉及的工藝參數包括肩部寬度 a1、鄰面預備深度 a2、頰面預備深度 a3、舌面預備深度 a4、咬合面預備深度 a5 和過渡斜面預備深度 a6。這里 a1、a2 與 a4 取 1 mm；a3 與 a5 取 1.5 mm；a6 取 0.5 mm。采用基于面的等距方法[21]，參考下頜磨牙的表面形貌，利用三維建模軟件 CREO 3.0（PTC Inc.，美國）進行預備體結合面參數化曲面造型。
（1）肩部預備：以齦緣曲線為掃描軌跡使用可變剖截面掃描命令，并在剖截面內控制肩部寬度參數 a1。
（2）鄰面和頰舌面預備：依據不同預備深度參數 a2、a3 與 a4 構造空間曲線，作為鄰面和頰舌面預備面的內部曲線，聯合肩部內環，使用曲面構造命令形成該預備面。
（3）咬合面預備：提取咬合面特征點（包括頰尖、舌尖和中央窩特征等），通過造型命令建立與咬合面特征相一致的近似曲面，并通過使用曲面偏移命令，控制咬合面預備深度參數 a5。
（4）斜面預備：以咬合面與軸面的交線為掃描軌跡使用可變剖截面掃描命令，并在剖截面內控制斜面預備深度參數 a6。
（5）光順打磨：使用圓角命令，實現交線的光順。
對磨牙表面和預備體結合面進行布爾運算與實體化操作，獲取義齒與預備體的三維模型，如圖 1 所示。
2 正中咬合優化
2.1 優化模型設計
如圖 2 所示，考慮到正中咬合為上下頜最廣泛的接觸型式，因此符合實際咬合狀況的下頜位姿應滿足接觸面積最大。以下頜位姿為設計變量，接觸面積最大為優化目標，建立正中咬合優化問題。以下頜點云模型幾何中心 O 為原點建立下頜連體坐標系與基坐標系。下頜模型的位姿使用連體坐標系相對于基坐標系的 6 個參數描述，即：V（x，y，z，α，β，γ），其中，（x，y，z）是點 O 的平移量，（α，β，γ）是歐拉角，依次為繞基坐標系進行 x 軸翻滾、y 軸俯仰和 z 軸偏航。已知，位移（x，y，z）取值范圍為［? 2，2］ mm，角度（α，β，γ）取值范圍為［? 2°，2°］。
圖2
基于接觸面積的目標函數定義
Figure2.
Definition of the objective function based on contact area
采用基于圖形空間的八叉樹-包圍球（octree-spheres）算法，對上下頜點云模型進行實時的碰撞檢測。檢測結果使用碰撞指標 h 表示，當發生碰撞時，h 取 1；反之，取 0。具體地，將當前點云模型均勻分割為 8 份規則的空間單元格；設計包含各空間單元格的最小球體作為包圍球；對 8 個包圍球進行碰撞檢測，排除未發生碰撞的包圍球，保留發生碰撞的包圍球（即該區域為可能碰撞區域）；對碰撞包圍球內點云模型重復上述流程。當排除所有包圍球時代表兩者之間未發生碰撞；當保留的包圍球達到半徑下限閾值時代表發生碰撞，并將該狀態定義為碰撞狀態。
對于非碰撞狀態，計算上頜點云與下頜各節點之間的最近距離di，如式（1）所示：
 |
其中，i 為下頜節點編號；Qi 為下頜節點 Pi 在上頜點云中的最近點，使用鄰近點搜索函數（nearest neighbor search，NNS）求解[22]。當節點最近距離的最小值不大于檢測距離，這部分節點定義為接觸節點，構成接觸點云，該狀態定義為接觸狀態；反之，定義為間隙狀態。
最后，根據不同咬合狀態，定義咬合優化目標 G，如式（2）所示：
 |
其中，min（·）為最小值函數；d0 為檢測距離，取 0.2 mm；S 為接觸點云的面積，如式（3）所示：
 |
其中，k 為接觸點云構成的三角面片編號；Ak、Bk 和 Ck 為組成三角面片 k 的三個節點坐標。
2.2 目標函數自動計算
咬合動力學分析軟件 OFA 內置點云的碰撞檢測、距離計算、面積計算等函數。該軟件可以計算出下頜模型沿咀嚼路徑起始點 D、若干中間迭代點 E 和終止點 F 運動時，點云模型間的實時接觸面積，如圖 3 所示。計算結果以規則化的表格呈現迭代步數 N 與接觸面積 SN 之間的關系。
圖3
咬合動力學分析軟件 OFA 計算流程
Figure3.
Calculating workflow in occlusal analysis software OFA
由于軟件無法脫離咀嚼路徑而直接輸出單個位姿的接觸面積，根據給定位姿參數 V 設計如下咀嚼路徑：定義起始點 D 位于下頜點云幾何中心 O，不設置中間迭代點 E，構造輔助參考點（位于起始點 D 的上方 0.01 mm）作為路徑終止點 F。利用軟件函數庫與項目庫構造腳本，自動實現起始點與終止點的位姿參數修改、咬合迭代分析和規則化表格輸出。
上述表格的規則化特點包括：當處于間隙狀態時，迭代步數 N 為空；當處于接觸狀態時，迭代步數 N 有非零值，輸出起始點 D 和接觸面積 S；當處于碰撞狀態時，軟件內部將不進行迭代運算，迭代步數 N 為零。據此使用數值分析軟件 Matlab R2016（MathWorks Inc.，美國）實現面積數據的提取與目標函數的輸出。
2.3 遺傳算法優化
鑒于遺傳算法作為智能優化算法具有可全局搜索且應用范圍廣的特點，因此本文嘗試使用其解決咬合分析這類復雜的優化問題。首先對設計變量進行編碼獲得初始種群，其中種群數量設置為 50。終止迭代條件為最大適應度函數（即目標函數 G 的相反數）變化范圍小于 0.001。前代種群如未滿足終止迭代條件，則進行基因選擇、交叉與變異操作，產生新種群，再次計算新種群是否滿足迭代終止條件。如此反復，最終輸出最佳種群并解碼出最佳設計變量。
在未進行優化前，初始頜位下接觸面積為 5.627 mm2。遺傳算法通過 3 050 次咬合迭代，輸出收斂的最優接觸面積為 15.477 mm2，如圖 4 所示。遺傳優化后得到的最佳下頜位姿參數為：［? 0.279，0.020，0.097，0.030，0.021，? 0.008］。同時采用咬合動力學分析軟件 OFA 自帶的頜位遍歷功能作為評價遺傳算法優化效果的對照項。試算過程顯示，該功能模塊通過 19 683 次咬合迭代，輸出最佳接觸面積為 12.918 mm2。對比發現，遺傳算法的優化效果明顯優于頜位遍歷功能，所得接觸面積提高 19.8%，且迭代步數減少 84.5%。
圖4
遺傳算法優化迭代過程
Figure4.
Iteration process of genetic algorithm
2.4 優化結果評價
使用咬合優化得到的最佳下頜位姿，在數值分析軟件 Matlab R2016 中修正下頜點云模型，并利用鄰近點搜索函數輸出與咬合動力學分析軟件 OFA 可視化結果相一致的節點位置。此時點云模型接觸點的數目為 717 個，集中在頰尖頰斜面與舌尖頰斜面。接觸區域可劃分為 4 個區域，包括接觸區域 1（橙色）、接觸區域 2（藍色）、接觸區域 3（綠色）和接觸區域 4（黃色），分別覆蓋 107、787、230 和 629 個三角單元，如圖 5 所示。
圖5
最佳下頜位姿在不同軟件中的接觸表現
Figure5.
Occlusal performance of the optimum mandibular position in different software
3 有限元建模與分析
3.1 工況設置
使用網格劃分軟件 HyperMesh 14.0（Altair Engineering Inc.，美國）對下頜第一磨牙全瓷冠以 0.2 mm 為全局單元尺寸進行網格劃分，形成有限元網格模型。根據口腔生物力學文獻給出的數據，設置義齒預備體與氧化鋯全瓷冠的物理參數：預備體密度為 2 140 kg/m3、彈性模量為 24.5 GPa、泊松比為 0.31；氧化鋯全瓷冠密度為 6 100 kg/m3、彈性模量為 220 MPa、泊松比為 0.30[23-24]。將下頜磨牙牙頸平面內的節點位移設置為 0，作為位移邊界條件。
為了進一步解決點云模型與網格模型節點不一致的問題，這里進行模型間接觸節點的映射。以點云模型接觸節點為出發點，使用最鄰近搜索算法尋找網格模型中的鄰近點。進一步刪除重復節點后，得到網格模型接觸節點的數目為 384 個。設計正中咬合法向力合值為 150 N，換算得到咬合節點力 0.39 N。利用節點鄰近點求解其所在曲面的近似法向量作為咬合力具體加載方向，如圖 6 所示。
圖6
網格模型的優化工況與參考工況設計
Figure6.
Design of optimized case and referenced case in finite element model
為了對咬合優化產生的影響進行討論，根據尖窩接觸的理想接觸點（中央窩、遠中窩、近中頰尖、遠中頰尖與遠中尖）設計參考工況。使用鄰近點搜索函數以理想接觸點為中心，擴展至周邊 77 個點（含自身），此時接觸點總數為 385，保持與優化工況基本一致。
3.2 應力結果討論
將網格模型結合載荷條件導入到有限元分析軟件 Abaqus 11.9（3DS Inc.，美國），運算得到如圖 7 所示的應力云圖。整體來看，新的咬合力加載方式顯著影響了應力分析結果。兩種工況義齒外表面應力較大區域與咬合位置基本保持一致，且最大應力位置存在明顯不同。對于優化工況，應力最大值為 12.12 MPa，發生在義齒舌側肩部；應力較大區域集中在頰尖舌斜面與頰斜面、舌尖頰斜面以及牙尖之間的溝形貌。對于參考工況，應力最大值為 16.02 MPa，發生在中央窩位置；應力較大區域集中在頰尖、中央窩等位置。
圖7
磨牙全瓷冠應力分布
Figure7.
Stress contour of the all-ceramic crown in the molar
觀察咬合節點位置圖可以發現，優化工況多為平面受力，而參考工況多為曲面受力，這使得前者相較于后者，具有更廣的咬合力傳遞空間。最終的應力云圖也顯示出，平面受力時應力分布更加均勻、連通，而曲面受力時應力分布更加集中、分散。此外，曲面受力呈現特殊的力傳遞特性，表現為：當曲面（如牙尖形貌）為凸形時，此時的最高點應力較小；而當曲面為凹形（如窩形貌）時，此時的最低點應力較大。
4 結論
本文考慮正中咬合優化對有限元模型載荷的影響，進行了下頜第一磨牙全瓷冠的三維建模、正中咬合優化與有限元分析，得到以下結論：遺傳算法可用于求解咬合優化問題，效率與性能明顯高于咬合動力學分析軟件 OFA 自帶的頜位遍歷功能；使用鄰近點搜索函數能夠實現咬合位置在點云模型與網格模型間接觸節點位置的映射，復用該函數可進一步求解接觸節點受力方向；最后，優化工況與參考工況具有明顯不同的應力極值與分布特征，主要原因在于前者多為平面受力，后者多為曲面受力。
在正中咬合重建的過程中，本文選取的優化目標是檢測距離為 0.2 mm 的咬合接觸面積，盡管優化效果良好，但仍值得進一步討論。首先，由于不同的檢測距離將呈現不同的咬合接觸面積，因此優化效果對檢測距離的敏感度有待后續確認。其次，正中咬合一系列定性特征[1-2]不僅體現在接觸面積上，同時也體現在接觸間隙、接觸力等定量咬合指標上，這使得以其他咬合指標為優化目標成為可能。
在咬合力設計的過程中，本文是在非規則接觸區域上使用均勻分布特征進行設計。根據咬合力分析設備的測試結果[25]，咬合力不僅加載區域是不規則的，而且大小也呈現非均勻分布。為了進一步提高咬合力載荷的仿生性，非均勻分布特征將在后續的研究中通過點云間相對距離換算來實現。
正中咬合優化不僅可以應用在單牙接觸的場景，對于全口牙列而言同樣具備應用的空間。此外，類似的接觸區域優化也可進一步應用于其他人體關節的受力分析，如顳下頜關節、膝關節和髖關節等。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
正中咬合（centric occlusion，CO）是牙列間最均勻、最廣泛、最緊密與最穩定的接觸，且伴隨著咀嚼循環中最大咬合力的發生，其作為關鍵咬合狀態影響著義齒的力學表現[1-2]。目前，有限元分析法已經成為評價義齒正中咬合狀態力學表現的重要途徑，可用于獲取感興趣部位的應力、應變和受力等數據[3]。考慮到載荷是影響有限元分析結果的基礎性因素[4]，因此盡可能建立與正中咬合狀態相吻合的載荷施加方式是提高義齒有限元建模準確性的關鍵，這就要求研究設計的技術方案能夠同時解決好正中咬合重建與咬合力設計這兩項內容。
由于不同的正中頜位（centric position，CP）將產生不同的咬合分析結果[5]，因此對于不確定的、初始的正中頜位有進行咬合重建的必要。并且，正中咬合的重建質量將在很大程度上決定咬合力設計的質量。常見的方法是，對人體或體外正確頜位關系采用光學[6-7]、探針[8]掃描記錄的方式在虛擬模型中建立正中咬合關系。目前，對咬合模型進行掃描已經成為牙頜模型掃描儀的標準操作步驟[9]，但該方法不適用于咬合關系不明的場景（例如：牙頜缺失）[10]。當無法直接獲取正中咬合狀態時，可以采用計算機輔助的頜位搜尋技術輔助研究。在咬合動力學分析軟件咬合指紋分析系統（occlusal fingerprint analyser，OFA）（ZiLoX IT Inc.，德國）中，具備基于接觸面積最大的頜位遍歷功能，但卻存在計算步數過多、解集質量不佳的問題[11]。同時，未見將正中咬合重建視為優化問題進行求解的報道。
咬合力設計，其目的是解決如何將咬合分析結果（例如：虛擬模型的接觸點云、實體模型的試紙壓痕等）設計映射到有限元網格模型。目前，咬合力加載方式包括集中力與均勻分布力兩種。其中，集中力的載荷施加方法，通常是將集中載荷以某個方向角施加在功能牙尖或牙脊上，因為凸起的形貌特征更容易成為碰撞的先發位置[12]。另外，根據尖窩接觸的特征，頰窩、近中窩和中央窩也常作為載荷施加點位[13-14]。然而，這種集中力的施加方式忽略了正中咬合復雜的接觸點分布特征，局部容易產生失真的應力集中現象。越來越多的研究使用均勻分布力作為載荷施加方式[15]。咬合力被施加在理想的中央窩區域或頰尖舌斜面，同時力的方向與牙齒長軸呈 45° 夾角。然而，此時的接觸點位置及接觸面邊界并非源于咬合分析結果，更多是經驗性的學知識，因此往往是規則的、理想的，與實際的接觸細節存在差異。為了克服這一不足，研究人員開始根據咬合動力學分析軟件或咬合試紙分析所得非規則的與實際吻合的接觸面，設計均勻分布式咬合力[16-19]，但是仍存在接觸位置映射至網格模型失準的情況。
基于以上原因，目前的正中咬合力加載方式同時存在正中咬合重建質量有限與咬合力加載位置失準的不足。為了解決這一問題，本文針對下頜第一磨牙全瓷冠有限元模型，創新性地提出一種基于接觸面積最大的正中咬合優化模型，并根據咬合優化結果設計基于網格模型的咬合力加載方式。通過本文研究，期望可以為義齒有限元模型提供新的咬合力加載方式，實現應力計算準確度的進一步提升，另一方面也可以為人體其他關節的受力分析提供新思路。
1 下頜第一磨牙全瓷冠光柵掃描與三維重建
1.1 光柵掃描與三維重建
牙頜模型是根據不同牙科教學或實驗目的仿照不同類型人體牙列的工業模型，其中標準牙頜模型是具有正常咬合面形貌的仿真牙列。本文選取大連醫科大學口腔醫學院提供的上下頜右側第一磨牙（D16FE-500HGSF-QF，Nissin. Inc，日本）作為人體磨牙的代替物進行掃描。具體掃描步驟如下：① 對掃描樣品進行熒光粉噴涂處理；② 使用光柵掃描設備（AutoScan DS100，先臨 3D，中國）分別對單顆牙進行掃描，記錄牙齒形貌，其掃描精度 0.015 mm。借助逆向工程軟件 Geomagic Studio（Geomagic Inc.，美國）對掃描獲取的牙齒表面點云模型進行多邊形破損修復和精確造面，實現咬合面、鄰面、頰面和舌面的三維重建，如圖 1 所示。
圖1
下頜第一磨牙全瓷冠三維模型
Figure1.
Three dimensional model of the mandibular first molar with all-ceramic crown prosthesis
1.2 預備體建模
根據全瓷冠基牙制備流程[20]，右側下頜第一磨牙預備工藝可分為肩部預備、鄰面和頰舌面預備、咬合面預備、斜面預備和光順打磨等 5 個步驟，涉及的工藝參數包括肩部寬度 a1、鄰面預備深度 a2、頰面預備深度 a3、舌面預備深度 a4、咬合面預備深度 a5 和過渡斜面預備深度 a6。這里 a1、a2 與 a4 取 1 mm；a3 與 a5 取 1.5 mm；a6 取 0.5 mm。采用基于面的等距方法[21]，參考下頜磨牙的表面形貌，利用三維建模軟件 CREO 3.0（PTC Inc.，美國）進行預備體結合面參數化曲面造型。
（1）肩部預備：以齦緣曲線為掃描軌跡使用可變剖截面掃描命令，并在剖截面內控制肩部寬度參數 a1。
（2）鄰面和頰舌面預備：依據不同預備深度參數 a2、a3 與 a4 構造空間曲線，作為鄰面和頰舌面預備面的內部曲線，聯合肩部內環，使用曲面構造命令形成該預備面。
（3）咬合面預備：提取咬合面特征點（包括頰尖、舌尖和中央窩特征等），通過造型命令建立與咬合面特征相一致的近似曲面，并通過使用曲面偏移命令，控制咬合面預備深度參數 a5。
（4）斜面預備：以咬合面與軸面的交線為掃描軌跡使用可變剖截面掃描命令，并在剖截面內控制斜面預備深度參數 a6。
（5）光順打磨：使用圓角命令，實現交線的光順。
對磨牙表面和預備體結合面進行布爾運算與實體化操作，獲取義齒與預備體的三維模型，如圖 1 所示。
2 正中咬合優化
2.1 優化模型設計
如圖 2 所示，考慮到正中咬合為上下頜最廣泛的接觸型式，因此符合實際咬合狀況的下頜位姿應滿足接觸面積最大。以下頜位姿為設計變量，接觸面積最大為優化目標，建立正中咬合優化問題。以下頜點云模型幾何中心 O 為原點建立下頜連體坐標系與基坐標系。下頜模型的位姿使用連體坐標系相對于基坐標系的 6 個參數描述，即：V（x，y，z，α，β，γ），其中，（x，y，z）是點 O 的平移量，（α，β，γ）是歐拉角，依次為繞基坐標系進行 x 軸翻滾、y 軸俯仰和 z 軸偏航。已知，位移（x，y，z）取值范圍為［? 2，2］ mm，角度（α，β，γ）取值范圍為［? 2°，2°］。
圖2
基于接觸面積的目標函數定義
Figure2.
Definition of the objective function based on contact area
采用基于圖形空間的八叉樹-包圍球（octree-spheres）算法，對上下頜點云模型進行實時的碰撞檢測。檢測結果使用碰撞指標 h 表示，當發生碰撞時，h 取 1；反之，取 0。具體地，將當前點云模型均勻分割為 8 份規則的空間單元格；設計包含各空間單元格的最小球體作為包圍球；對 8 個包圍球進行碰撞檢測，排除未發生碰撞的包圍球，保留發生碰撞的包圍球（即該區域為可能碰撞區域）；對碰撞包圍球內點云模型重復上述流程。當排除所有包圍球時代表兩者之間未發生碰撞；當保留的包圍球達到半徑下限閾值時代表發生碰撞，并將該狀態定義為碰撞狀態。
對于非碰撞狀態，計算上頜點云與下頜各節點之間的最近距離di，如式（1）所示：
|
其中，i 為下頜節點編號；Qi 為下頜節點 Pi 在上頜點云中的最近點，使用鄰近點搜索函數（nearest neighbor search，NNS）求解[22]。當節點最近距離的最小值不大于檢測距離，這部分節點定義為接觸節點，構成接觸點云，該狀態定義為接觸狀態；反之，定義為間隙狀態。
最后，根據不同咬合狀態，定義咬合優化目標 G，如式（2）所示：
|
其中，min（·）為最小值函數；d0 為檢測距離，取 0.2 mm；S 為接觸點云的面積，如式（3）所示：
|
其中，k 為接觸點云構成的三角面片編號；Ak、Bk 和 Ck 為組成三角面片 k 的三個節點坐標。
2.2 目標函數自動計算
咬合動力學分析軟件 OFA 內置點云的碰撞檢測、距離計算、面積計算等函數。該軟件可以計算出下頜模型沿咀嚼路徑起始點 D、若干中間迭代點 E 和終止點 F 運動時，點云模型間的實時接觸面積，如圖 3 所示。計算結果以規則化的表格呈現迭代步數 N 與接觸面積 SN 之間的關系。
圖3
咬合動力學分析軟件 OFA 計算流程
Figure3.
Calculating workflow in occlusal analysis software OFA
由于軟件無法脫離咀嚼路徑而直接輸出單個位姿的接觸面積，根據給定位姿參數 V 設計如下咀嚼路徑：定義起始點 D 位于下頜點云幾何中心 O，不設置中間迭代點 E，構造輔助參考點（位于起始點 D 的上方 0.01 mm）作為路徑終止點 F。利用軟件函數庫與項目庫構造腳本，自動實現起始點與終止點的位姿參數修改、咬合迭代分析和規則化表格輸出。
上述表格的規則化特點包括：當處于間隙狀態時，迭代步數 N 為空；當處于接觸狀態時，迭代步數 N 有非零值，輸出起始點 D 和接觸面積 S；當處于碰撞狀態時，軟件內部將不進行迭代運算，迭代步數 N 為零。據此使用數值分析軟件 Matlab R2016（MathWorks Inc.，美國）實現面積數據的提取與目標函數的輸出。
2.3 遺傳算法優化
鑒于遺傳算法作為智能優化算法具有可全局搜索且應用范圍廣的特點，因此本文嘗試使用其解決咬合分析這類復雜的優化問題。首先對設計變量進行編碼獲得初始種群，其中種群數量設置為 50。終止迭代條件為最大適應度函數（即目標函數 G 的相反數）變化范圍小于 0.001。前代種群如未滿足終止迭代條件，則進行基因選擇、交叉與變異操作，產生新種群，再次計算新種群是否滿足迭代終止條件。如此反復，最終輸出最佳種群并解碼出最佳設計變量。
在未進行優化前，初始頜位下接觸面積為 5.627 mm2。遺傳算法通過 3 050 次咬合迭代，輸出收斂的最優接觸面積為 15.477 mm2，如圖 4 所示。遺傳優化后得到的最佳下頜位姿參數為：［? 0.279，0.020，0.097，0.030，0.021，? 0.008］。同時采用咬合動力學分析軟件 OFA 自帶的頜位遍歷功能作為評價遺傳算法優化效果的對照項。試算過程顯示，該功能模塊通過 19 683 次咬合迭代，輸出最佳接觸面積為 12.918 mm2。對比發現，遺傳算法的優化效果明顯優于頜位遍歷功能，所得接觸面積提高 19.8%，且迭代步數減少 84.5%。
圖4
遺傳算法優化迭代過程
Figure4.
Iteration process of genetic algorithm
2.4 優化結果評價
使用咬合優化得到的最佳下頜位姿，在數值分析軟件 Matlab R2016 中修正下頜點云模型，并利用鄰近點搜索函數輸出與咬合動力學分析軟件 OFA 可視化結果相一致的節點位置。此時點云模型接觸點的數目為 717 個，集中在頰尖頰斜面與舌尖頰斜面。接觸區域可劃分為 4 個區域，包括接觸區域 1（橙色）、接觸區域 2（藍色）、接觸區域 3（綠色）和接觸區域 4（黃色），分別覆蓋 107、787、230 和 629 個三角單元，如圖 5 所示。
圖5
最佳下頜位姿在不同軟件中的接觸表現
Figure5.
Occlusal performance of the optimum mandibular position in different software
3 有限元建模與分析
3.1 工況設置
使用網格劃分軟件 HyperMesh 14.0（Altair Engineering Inc.，美國）對下頜第一磨牙全瓷冠以 0.2 mm 為全局單元尺寸進行網格劃分，形成有限元網格模型。根據口腔生物力學文獻給出的數據，設置義齒預備體與氧化鋯全瓷冠的物理參數：預備體密度為 2 140 kg/m3、彈性模量為 24.5 GPa、泊松比為 0.31；氧化鋯全瓷冠密度為 6 100 kg/m3、彈性模量為 220 MPa、泊松比為 0.30[23-24]。將下頜磨牙牙頸平面內的節點位移設置為 0，作為位移邊界條件。
為了進一步解決點云模型與網格模型節點不一致的問題，這里進行模型間接觸節點的映射。以點云模型接觸節點為出發點，使用最鄰近搜索算法尋找網格模型中的鄰近點。進一步刪除重復節點后，得到網格模型接觸節點的數目為 384 個。設計正中咬合法向力合值為 150 N，換算得到咬合節點力 0.39 N。利用節點鄰近點求解其所在曲面的近似法向量作為咬合力具體加載方向，如圖 6 所示。
圖6
網格模型的優化工況與參考工況設計
Figure6.
Design of optimized case and referenced case in finite element model
為了對咬合優化產生的影響進行討論，根據尖窩接觸的理想接觸點（中央窩、遠中窩、近中頰尖、遠中頰尖與遠中尖）設計參考工況。使用鄰近點搜索函數以理想接觸點為中心，擴展至周邊 77 個點（含自身），此時接觸點總數為 385，保持與優化工況基本一致。
3.2 應力結果討論
將網格模型結合載荷條件導入到有限元分析軟件 Abaqus 11.9（3DS Inc.，美國），運算得到如圖 7 所示的應力云圖。整體來看，新的咬合力加載方式顯著影響了應力分析結果。兩種工況義齒外表面應力較大區域與咬合位置基本保持一致，且最大應力位置存在明顯不同。對于優化工況，應力最大值為 12.12 MPa，發生在義齒舌側肩部；應力較大區域集中在頰尖舌斜面與頰斜面、舌尖頰斜面以及牙尖之間的溝形貌。對于參考工況，應力最大值為 16.02 MPa，發生在中央窩位置；應力較大區域集中在頰尖、中央窩等位置。
圖7
磨牙全瓷冠應力分布
Figure7.
Stress contour of the all-ceramic crown in the molar
觀察咬合節點位置圖可以發現，優化工況多為平面受力，而參考工況多為曲面受力，這使得前者相較于后者，具有更廣的咬合力傳遞空間。最終的應力云圖也顯示出，平面受力時應力分布更加均勻、連通，而曲面受力時應力分布更加集中、分散。此外，曲面受力呈現特殊的力傳遞特性，表現為：當曲面（如牙尖形貌）為凸形時，此時的最高點應力較小；而當曲面為凹形（如窩形貌）時，此時的最低點應力較大。
4 結論
本文考慮正中咬合優化對有限元模型載荷的影響，進行了下頜第一磨牙全瓷冠的三維建模、正中咬合優化與有限元分析，得到以下結論：遺傳算法可用于求解咬合優化問題，效率與性能明顯高于咬合動力學分析軟件 OFA 自帶的頜位遍歷功能；使用鄰近點搜索函數能夠實現咬合位置在點云模型與網格模型間接觸節點位置的映射，復用該函數可進一步求解接觸節點受力方向；最后，優化工況與參考工況具有明顯不同的應力極值與分布特征，主要原因在于前者多為平面受力，后者多為曲面受力。
在正中咬合重建的過程中，本文選取的優化目標是檢測距離為 0.2 mm 的咬合接觸面積，盡管優化效果良好，但仍值得進一步討論。首先，由于不同的檢測距離將呈現不同的咬合接觸面積，因此優化效果對檢測距離的敏感度有待后續確認。其次，正中咬合一系列定性特征[1-2]不僅體現在接觸面積上，同時也體現在接觸間隙、接觸力等定量咬合指標上，這使得以其他咬合指標為優化目標成為可能。
在咬合力設計的過程中，本文是在非規則接觸區域上使用均勻分布特征進行設計。根據咬合力分析設備的測試結果[25]，咬合力不僅加載區域是不規則的，而且大小也呈現非均勻分布。為了進一步提高咬合力載荷的仿生性，非均勻分布特征將在后續的研究中通過點云間相對距離換算來實現。
正中咬合優化不僅可以應用在單牙接觸的場景，對于全口牙列而言同樣具備應用的空間。此外，類似的接觸區域優化也可進一步應用于其他人體關節的受力分析，如顳下頜關節、膝關節和髖關節等。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
