電子計算機斷層掃描（CT）引導下經皮肺穿刺是獲取肺部病變組織以確診肺癌的最有效的手段之一，規劃科學合理的穿刺路徑是避免穿刺并發癥，減少患者痛苦及穿刺死亡率的重要環節。本文提出了一種基于多級約束的肺部穿刺路徑規劃方法：首先利用患者CT建立胸部數字化模型；然后在以腫瘤病灶為中心的理想球面上構造斐波那契網格采樣，得到待選路徑集合；最后根據臨床穿刺準則提出了一種多級約束策略，并結合方向包圍盒層次樹（OBBTree）算法和帕累托（Pareto）優化算法，獲取最優穿刺路徑。模擬仿真實驗結果證明了該算法的有效性，能夠規避物理和生理障礙，可以作為醫生選擇穿刺路徑的輔助手段。
引用本文： 孫峰輝, 裴宏亮, 楊逸飛, 樊慶文, 李曉歐. 一種基于多級約束的肺部穿刺路徑規劃方法. 生物醫學工程學雜志, 2022, 39(3): 462-470. doi: 10.7507/1001-5515.202112029 復制
引言
肺癌是最常見的惡性腫瘤之一，病理活檢是肺癌診斷的金標準，經皮肺穿刺是提取病理組織的可靠手段。經皮肺穿刺是利用穿刺針穿透皮膚組織、肌肉組織和肺部組織達到腫瘤病灶，提取病變組織以進行病理檢查，確診是否為惡性腫瘤的一種常用方法[1-5]。在穿刺過程中，需要避開骨骼，肺部氣管支氣管樹分支和較大的肺血管，否則會引起氣胸、血胸、血氣胸、肺出血、大咯血、肺部感染、病灶播散，嚴重時甚至導致患者死亡[6-7]。數據調查表明，氣胸和肺出血是經皮肺穿刺手術的主要并發癥，其中氣胸的發生率為10%～40%，肺出血的發生率為26%～33%[8]。為了降低并發癥的發生率，根據腫瘤病灶的位置設計合理的穿刺路徑是肺穿刺領域的重要研究內容之一[9]。
目前，確定經皮肺穿刺路徑主要依靠醫生借助電子計算機斷層掃描（computed tomography，CT）圖像憑經驗手工完成。隨著計算機圖像處理和可視化技術的發展，通過CT圖像獲得最優的安全穿刺路徑成為了可能。穿刺手術路徑規劃算法的研究對象主要分為柔性針和剛性針兩類。柔性針路徑規劃算法借助柔性針本身可以改變進針方向的特點，試圖尋找復雜的曲線路徑以避開障礙，主要包括逆運動學反解法、人工勢場法和數值法三類[10]。柔性針路徑規劃算法的計算較為復雜，在實際穿刺過程中也需要構建符合實際情況的針—組織交互模型，同時考慮到柔性針本身的特點，目前仍處于實驗驗證階段。
相比于柔性針，剛性針路徑規劃算法的穿刺路徑為直線，是目前人工穿刺和機器人穿刺的主要手段。Kimura等[11]利用多平面重建（multiplanar reconstruction，MPR）技術重建患者CT數據集的特殊截面，引導穿刺針對處于復雜區域的腫瘤病灶進行斜向穿刺，但是穿刺路徑的規劃依舊需要人工完成。為了實現全自動剛性針路徑規劃， Villard等[12-13]、Baegert等[14-16]和Seitel等[17]以患者器官與組織的面繪制模型為基礎，針對臨床約束條件分別利用三角形網格曲面繪制方法對可行進針區域進行篩選與優化；Schumann等[18-21]以腫瘤病灶為原點，通過圓柱投影對不同優化目標生成二維約束地圖，采取加權乘積、手工交互或帕累托優化（Pareto optimization）的方法選擇穿刺路徑；張睿等[22]基于改進的立方體映射算法，在患者CT圖像的包圍盒表面進行二值化分類，得到滿足臨床約束條件的可行穿刺區域，再針對不同的臨床穿刺準則進行加權整合，提供若干穿刺路徑供醫生選擇。但是，上述剛性針路徑規劃算法往往會將數字化人體三維模型投影到某個特殊曲面，這個曲面可能是CT數據集的特殊截面、皮膚模型、圓柱面或立方體面，再基于曲面上的信息規劃穿刺路徑。其本質都是對原三維模型復雜空間關系的映射過程，不可避免地會發生信息失真，使投影到曲面的圖像發生畸變，難以確定最優穿刺路徑。另一方面，Gao等[23]以腫瘤病灶為中心，通過構建球面采樣點的方法直接在三維模型上進行路徑篩選，再根據入刺深度、目標區域可視大小等標準對穿刺路徑進行評分，避免了三維模型向曲面映射時所可能產生的信息失真問題。但是作者并未說明球面采樣點的具體生成方式，也未討論如何確定采樣點的數目，同時，不完整的皮膚模型使“入刺深度”這一評價標準僅適用于部分穿刺路徑，存在一定的缺陷。
針對以上問題，本文提出了一種基于多級約束的肺部穿刺路徑規劃方法。在患者數字化人體三維模型的基礎上，首先以腫瘤病灶為中心建立理想球面并在其上利用斐波那契網格（Fibonacci lattice）生成采樣點，其中每一個采樣點到病灶中心均構成一條可能的穿刺路徑；然后根據臨床穿刺準則對上述路徑進行約束，主要分為物理約束、生理約束和數理約束三部分：先根據體位、工具、設備等物理條件的限制排除大部分路徑；再利用方向包圍盒層次樹（oriented bounding box tree，OBBTree）模擬穿刺過程以規避風險組織和器官，篩選出滿足生理約束條件的穿刺路徑；最后在數理約束階段根據臨床穿刺準則對篩選出的穿刺路徑進行量化評價，并利用帕累托優化算法得到最優穿刺路徑。
1 方法
本文所提算法的流程如圖1所示。在胸部數字化模型上以腫瘤病灶為中心，穿刺針長度為半徑構建理想球面，通過斐波那契網格采樣得到離散點點集S0，其中每一點與病灶中心的連線均構成一條待選路徑；首先根據體位、工具、設備對穿刺路徑的限制完成物理約束篩選，得到路徑集合S1；再采用OBBTree算法模擬穿刺過程，篩選出滿足生理約束的穿刺路徑集合S2；最后采用帕累托優化算法對滿足物理和生理約束的穿刺路徑進行量化評價，得到備選路徑集合S3，選擇其中評分最高的路徑作為最優穿刺路徑s。
圖1
基于多級約束的肺部穿刺路徑規劃算法流程圖
Figure1.
Flow chart of the algorithm for lung puncture path planning based on the multi-level constraint
1.1 斐波那契網格采樣
斐波那契網格是由解析公式生成的一組點集，其中的點沿著密繞螺旋排列，在螺旋軸線方向上，每個點與上一點之間的角度差均構成黃金分割比，從宏觀上看，這些點所構成的點集在球面或圓平面上以基本各向同性的方式均勻分布[24]。利用其在球面或圓平面上點集分布的均勻性質，斐波那契網格已經在球面不規則圖形面積測量[24]、顏色量化與處理[25]以及相控陣天線布板設計[26]等方面得到應用。因此可以利用斐波那契網格對球面進行離散化采樣，以此確定待測穿刺路徑集合，記為點集M。如式(1)所示，設球面半徑為R，點集M中的元素總數為N，對于點集M中的任一點m(m = {(xn, yn, zn), n ∈ [1, N]})有：
 |
其中，φ為定值，取
，n為1～N范圍內的整數，(xn, yn, zn)表示由n確定的m點坐標。
斐波那契網格的構成需要理想球面的半徑R與元素總數N作為輸入。根據穿刺針長度可以確定半徑R的大小，但在穿刺路徑的計算中往往并不關注元素總數N的多少，而更關注穿刺路徑分布的疏密程度，即點集中點與點之間的最近距離大小，記為分辨率δ。因此就需要得到分辨率δ與元素總數N的關系，以便在給出δ后可以直接估算出所需的N值。假設以各個點為圓心，δ/2為半徑在球面上構成一系列圓盤，在理想的絕對均勻條件下，這一系列圓盤在數學上形成了對球體的最密封裝，對于元素總數N的估計如式(2)所示：
 |
其中，SM表示所有圓盤的總面積；Sm表示單個圓盤的面積；S表示整個球面的面積。由于S略大于SM，故利用此公式計算得出的N也略大于理想的元素總數，因此該公式可以在保證分辨率要求的情況下，用于對元素總數N的估計。
以上斐波那契網格采樣所生成的是均勻分布于整個球面上的點集，從點集中的點m到穿刺目標點所形成的路徑中，有的并不能滿足臨床穿刺的要求。因此就需要對斐波那契網格采樣所產生的點集進行多級約束篩選，最終得到最優穿刺路徑。
1.2 基于臨床穿刺準則的多級約束設計
臨床穿刺準則是穿刺醫生在實踐中總結出的經驗標準，其中規定了穿刺操作的限制條件以及穿刺所要達到的理想目標。使用計算機技術進行肺部穿刺路徑規劃，需要對臨床穿刺準則進行量化分析。本文將臨床穿刺準則歸納為三個部分：物理約束、生理約束和數理約束。
1.2.1 物理約束
物理約束是指在肺部穿刺時由于體位、工具、設備的原因對穿刺路徑造成的限制，在規劃路徑時必須予以規避。雖然理論上醫生可以在任意方向實施任意距離的穿刺，但是由于患者在穿刺設備上所采取的體位、穿刺針長度和穿刺角度等因素的影響，客觀上限制了穿刺路徑的選擇。物理約束主要可分為三類：
（1）患者在穿刺設備上所采取的體位。CT引導經皮肺穿刺的常用體位包括仰臥位、俯臥位、左側臥位及右側臥位[27]。通過對腫瘤病灶位置的分析，醫生可以初步確定采用何種體位，由于設備的遮擋關系，不同體位所允許的穿刺角度區間也不同。以仰臥位為例，可以將胸部前方區域作為合適的穿刺區域。
（2）穿刺針本身的長度限制。在肺穿刺手術中選擇路徑的長度一定不會超過穿刺針的長度。假設穿刺針穿刺時可達的最大長度為R（考慮到醫生在臨床操作時的實際情況，保險起見R會比穿刺針的真實長度略短），即不必考慮超過R范圍的穿刺路徑，因為它們超過了穿刺針的最大作業范圍，無法成功完成穿刺。如式(3)所示，假設皮膚表面上任意點P(x, y, z)與選定的穿刺目標點O(x0, y0, z0)之間的總長度為d，則有：
 |
當d ≤ R時，可以認為從該點到穿刺目標點對應的路徑滿足了穿刺針長度約束，予以保留，否則剔除該點。對于滿足穿刺針長度約束的穿刺路徑，從皮膚到腫瘤病灶的長度越小，對患者造成的損傷越小，穿刺風險越低。在斐波那契網格采樣過程中，已經以腫瘤病灶為中心，穿刺針長度為半徑構建了理想球面，其上每一采樣點到病灶中心的距離均為R。結合數字化人體模型，只有當采樣點位于患者皮膚外時，對應路徑從入刺點到病灶中心的距離d在R范圍內，予以保留，否則剔除該點。
（3）穿刺針與皮膚的夾角限制。在穿刺針與患者皮膚作用的過程中，兩者之間的夾角需要大于臨床常規閾值（20°）以避免滑針[28]。如式(4)所示，設皮膚表面上任意點P有單位法向量N(i, j, k)，選定的穿刺目標點與備選點之間構成單位方向向量Z(a, b, c)，可以求得穿刺針進針角度θ為：
 |
根據臨床常規閾值，當θ ≥ 20°時，可以認為穿刺針穿刺該點時不會滑針，予以保留，否則剔除該點。對于滿足穿刺針角度約束的穿刺路徑，它的方向向量與皮膚法向量之間的夾角越小，滑針的概率越小，穿刺風險越低。
穿刺針長度約束和角度約束示意圖如圖2所示，紅色區域為腫瘤病灶，其中心為穿刺目標點O，橙色線條為皮膚的橫截面。假設經過斐波那契網格采樣產生三條路徑，其中路徑1（綠色線條）和路徑2（紅色線條）的穿刺深度分別為d1和d2，穿刺角度為θ1和θ2，可以看出路徑1更好地滿足了物理約束條件；路徑3（藍色線條）由于長度不足，未能穿出皮膚，需要予以剔除。
圖2
穿刺針長度約束與角度約束
Figure2.
Length constraint and angle constraint of puncture needle
除了上述總結的約束類型，為了方便醫生更好地操作工具和設備、觀察穿刺路徑，還可以針對具體手術環境設置額外的物理約束。
1.2.2 生理約束
生理約束是指為了大幅度降低穿刺時氣胸、肺出血等并發癥的發生概率而對穿刺路徑進行的限制。這就需要穿刺針在到達腫瘤病灶的過程中規避骨骼、肺部大氣管、大血管、心臟等風險組織和器官。在重建的數字化人體模型基礎上，可以將骨骼、氣管、血管、心臟、肝臟等組織和器官視為穿刺路徑所不能通過的障礙，再進行生理約束的計算。如式(5)所示，假設點集中有任意點P(x, y, z)，選定的穿刺目標點為O(x0, y0, z0)，可以確定穿刺路徑的單位方向向量Z(a, b, c)為：
 |
其中，d為點P與點O之間的總長度。從穿刺目標點O(x0, y0, z0)沿著向量Z的方向進行路徑搜索，如果該路徑被重要器官或組織遮擋則剔除該路徑，否則繼續沿著原方向進行搜索，如果可以到達備選點P(x, y, z)，則證明該路徑滿足重要器官、組織避障約束，最終予以保留。對于滿足生理約束的穿刺路徑，它與最近障礙之間的直線距離越大，成功規避風險組織和器官的概率越大，穿刺風險越低。
1.2.3 數理約束
數理約束是在滿足前兩級約束的基礎上，為了進一步降低穿刺風險，提高穿刺可靠性，運用數學工具從可能的穿刺路徑中篩選最優路徑的過程。主要內容包括：①依照臨床穿刺準則對穿刺路徑進行量化評價。②根據帕累托優化結果得出最優穿刺路徑。
經過物理約束和生理約束篩選的穿刺路徑在理論上均滿足穿刺手術的基本要求，但是不同路徑對約束條件的滿足程度各不相同。有的路徑雖然滿足了物理、生理約束，但是距離風險組織和器官依然很近，考慮到實際穿刺時的操作誤差，此類路徑的手術風險過大，需要予以剔除。因此，為了能夠選出最優穿刺路徑，需要依照臨床穿刺準則對穿刺路徑進行量化評價。
根據物理約束和生理約束，定義從皮膚入刺點到病灶的路徑總長L、路徑與皮膚法向量夾角角度Θ和路徑與風險組織器官的最近距離D作為穿刺路徑的評價標準。其中L和Θ的大小可以根據穿刺路徑與皮膚的空間關系直接得出，為了消除量綱差異，使用離差標準化方法對上述兩個標準進行歸一化處理，如式(6)所示，對路徑i有：
 |
其中，Li、Θi是當前路徑的長度和夾角，Lmin、Θmin是所有有效路徑長度和夾角的最小值，Lmax、Θmax是所有有效路徑長度和夾角的最大值。由定義可知，L(i)和Θ(i)的值越小，對應路徑的評價越高，其風險越小。
而求解各路徑的D值涉及到穿刺路徑與風險組織器官之間的復雜空間關系，很難進行直接的計算，因此需要采用其他方法間接地對D值進行量化。設通過生理約束篩選的點集為A，未能通過生理約束篩選的點集為 
。相比于直接計算穿刺路徑的D值，找出點集A中與點集 
的距離在δ以內的點更加容易。因此，可以通過這種方式得到點集A所代表的穿刺路徑中D值小于δ的部分，記為A1。以此類推，將A1合并到 
中，重復上述計算步驟，可以得到A2、A3、A4、
，在這些集合中穿刺路徑的D值越來越大，風險越來越低。對點集A1、A2、A3、A4、…所代表的穿刺路徑分別賦予D值1、2、3、4、
，即完成了對D值的量化。由定義可知，D(i)的值越大，對應路徑的風險越小。
在對每一條穿刺路徑進行量化評價后，需要依據評分高低選擇出最優穿刺路徑。但是不同的評價標準往往存在相互制約的關系，在某一標準上表現較好的穿刺路徑，在其他標準上可能具有較差的表現。為了解決上述問題，可以為不同的評價標準分別設置權重系數，將經由權重系數調整的評分加和作為穿刺路徑的總體評價。這種基于權重的評價方法可以直接選出最優穿刺路徑，但是需要人為對權重系數進行設置，而權重系數的設置并沒有一個明確的標準，這為穿刺路徑的規劃帶來了主觀因素的影響。
為了避免上述方法的缺點，可以采用帕累托優化的方法得出最優穿刺路徑[29]，通過此方法可以在不額外設置系數的前提下求解得到位于帕累托前沿面上的解，這些位于帕累托前沿面上的解之間無法進一步區分優劣，而在多個評價標準中均優于非帕累托前沿面上的解。
2 實驗與結果
2.1 實驗數據
肺部穿刺路徑規劃算法必須以患者胸部的數字化模型為基礎，該模型應包括肺穿刺需要規避的重要組織和器官，數字化模型的精度直接影響到穿刺路徑規劃的效果。目前重建胸部數字化模型的素材主要來自于患者的CT數據集，隨著計算機技術，特別是神經網絡的發展，針對CT數據集的自動分割算法已經能夠對包括肺部精細化結構在內的胸部組織器官進行高效的三維重建，可以滿足本文所提算法的需求，在獲取CT影像后及時構建模體并進行手術[30-32]。本文研究重點在于肺部穿刺路徑規劃算法，因此依然采用人工分割的方式建立數字化人體模型。本實驗所用數據來源于美國國家癌癥研究所發布的公開數據集肺部圖像數據庫聯盟與圖像數據庫資源計劃（lung image database consortium and image database resource initiative，LIDC-IDRI）（網址為：https://wiki.cancerimagingarchive.net/display/Public/LIDC-IDRI）[33]。數字化人體模型的建立過程如圖3所示，采用CT圖像的像素尺寸為0.75 mm × 0.75 mm，層厚1 mm，切片數341層，層分辨率為512 × 512像素，利用醫學圖像處理軟件3Dslicer4.11（The Slicer Community，美國）重建了包括患者皮膚、骨骼、肺（肺實質、肺氣管、肺血管）、心臟和部分肝在內的三維模型，最后裝配成為患者胸部數字化模型，數字化模型中的黑色區域為模擬腫瘤的位置[34]。
圖3
數字化人體模型的建立過程
Figure3.
The process of building the digital human body model
另外，在原有CT圖像坐標系O’-ijk中，CT圖片左上角為坐標原點，向右為i軸正方向，向下為j軸正方向，向后為k軸正方向。為了更方便地表示穿刺路徑并減少計算量，需要建立以病灶中心O(x0, y0, z0)為原點的肺穿刺路徑規劃坐標系O-xyz，兩個坐標系可以互相轉換。
2.2 斐波那契網格采樣的創建
以患者胸部數字化模型為基礎，在肺穿刺路徑規劃坐標系O-xyz中進行斐波那契網格采樣。為此需要確定三個參數：理想球面的球心，理想球面的半徑R以及分辨率δ。理想球面的球心就是腫瘤病灶的中心O(x0, y0, z0)，其半徑R可由穿刺針長度確定，穿刺精度決定了分辨率δ的大小。以人工穿刺為例，在CT引導下經皮肺穿刺活檢中使用的穿刺針長度約為21 cm[28]，針尖偏移距離誤差約為5 mm[35]。根據以上數據，有R = 210 mm，δ = 5 mm。
2.3 基于多級約束的肺部穿刺路徑規劃
2.3.1 物理約束篩選
物理約束可以進一步縮小穿刺路徑的可行區域，其結果如圖4所示，相比于整個理想球面，由綠色區域代表的穿刺針可行區域大大減少，在該區域內的斐波那契網格采樣點到穿刺目標點之間的連線均構成一條待選穿刺路徑。
圖4
滿足物理約束的路徑集合
Figure4.
The set of paths satisfying the physical constraints
2.3.2 生理約束篩選
物理約束篩選的結果構成一離散化點集，以其中的點為一端，穿刺目標點為另一端模擬穿刺過程，利用OBBTree算法篩選出滿足生理約束的穿刺路徑。OBBTree是一種多邊形有向包圍盒的二叉樹數據結構和算法，用于物體之間的快速相交測試[36]。利用基于OBBTree的快速相交測試算法，通過判斷路徑和重要器官或組織是否存在交點，可以選擇出規避風險區域的穿刺路徑。如圖5所示，由綠色區域代表的穿刺針可行區域進一步減少，經過放大可以看到，點集輪廓（黑色邊緣）能夠完全反映風險區域對腫瘤病灶的遮擋范圍（紅色邊緣），完成了對風險區域的規避。
圖5
滿足生理約束的路徑集合
Figure5.
The set of paths satisfying the physiological constraints
2.3.3 數理約束篩選
對于從生理約束中篩選出的大量路徑，需要通過數理約束對其進行量化評價，得出最優穿刺路徑。在各項評價標準中，從皮膚入刺點到病灶的路徑總長L和路徑與皮膚法向量夾角角度Θ比較容易計算得出；各路徑與風險組織器官的最近距離D如圖6所示，可以看到，隨著D值的增大，穿刺路徑與風險組織器官的距離同時增加，穿刺風險降低。
圖6
各路徑與風險組織器官的最近距離
Figure6.
Nearest distance of each pathway to the risk tissue and organ
考慮到D值較小的路徑過于靠近風險組織器官，因此選擇在D值較大的路徑中采用帕累托優化方法篩選最優穿刺路徑。如圖7所示，x軸表示從皮膚入刺點到病灶的路徑總長評分L(i)，y軸表示路徑與皮膚法向量夾角角度評分Θ(i)，兩者均由式(6)計算得出。針對不同的D值得到的帕累托前沿面顯示為紅線，在帕累托前沿面上的點可以視為該D值下的最優點。可以看到，在帕累托前沿面p1、p2和p3中，p1上的穿刺路徑具有最短的穿刺距離和最好的穿刺角度，但是距離風險組織器官較近；相反，p3上的穿刺路徑距離風險組織器官最遠，但穿刺距離和穿刺角度都稍遜于p1和p2上的穿刺路徑。
圖7
不同D值下的帕累托前沿面
Figure7.
Pareto front surface for different D values
將處于各帕累托前沿面上的點所代表的穿刺路徑表示在患者數字化人體模型上，如圖8所示，位于p1上的2條穿刺路徑標記為淺藍色，位于p2上的2條穿刺路徑標記為深藍色，位于p3上的2條穿刺路徑標記為紅色，這六條路徑之間沒有明確的優劣關系，而在三個評價標準上均優于其他穿刺路徑，醫生可以在上述六條路徑中進行選擇，大大減少了主觀因素的影響。
圖8
路徑規劃的最終結果
Figure8.
The final result of path planning
2.4 實驗結果
為了驗證算法的穩定性，基于不同體型、年齡、性別患者的CT數據建立了另外9套數字化人體模型，分別為模型1～模型9，隨機設定病灶原點，并運用本文所提出的方法進行肺部穿刺路徑規劃。如圖9所示，最終篩選出的穿刺路徑標記為黑色，其數目在4～12條之間，平均在每個模型中可以計算得到7條穿刺路徑。在子圖模型2、模型3、模型4、模型6、模型8中，篩選出的穿刺路徑具有一定的聚集性，表示最優穿刺路徑即在此區域內；在子圖模型5、模型7、模型9中，篩選出的穿刺路徑比較分散，表示在已有的穿刺標準下有多個可能的穿刺方案；在子圖模型1中，篩選出的穿刺路徑總體具有分散的特點，而在個別區域又有一定的聚集性。實驗結果證明，該算法具有一定的穩定性，數量適中的穿刺路徑可以在幫助醫生確定穿刺范圍的同時，提供一定的選擇方案。
圖9
基于不同患者模型的路徑規劃結果
Figure9.
Path planning results based on different patient models
3 總結
本文提出了一種基于多級約束的肺部穿刺路徑規劃算法，在患者數字化人體模型的基礎上，首先在以腫瘤病灶為中心的理想球面上構造斐波那契網格采樣生成穿刺路徑，根據體位、工具、設備對穿刺造成的限制完成物理約束的篩選，利用OBBTree算法篩選滿足生理約束條件的穿刺路徑，最后利用帕累托優化算法得到了有代表性的肺部穿刺路徑。經過實驗驗證，該算法可以規避骨骼、氣管、血管等風險組織，能夠選擇出相對安全的穿刺路徑供醫生選擇，同時也為未來實現機器人穿刺提供了技術支撐。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：孫峰輝主要負責算法程序設計、實驗數據分析以及論文初稿撰寫；裴宏亮主要負責實驗方法設計和實驗結果可視化；楊逸飛主要負責數據整理與管理；樊慶文主要負責提供實驗指導，數據分析指導和論文審閱修訂；李曉歐主要負責提供醫學指導以及對文章的醫學內容作批評性審閱。
引言
肺癌是最常見的惡性腫瘤之一，病理活檢是肺癌診斷的金標準，經皮肺穿刺是提取病理組織的可靠手段。經皮肺穿刺是利用穿刺針穿透皮膚組織、肌肉組織和肺部組織達到腫瘤病灶，提取病變組織以進行病理檢查，確診是否為惡性腫瘤的一種常用方法[1-5]。在穿刺過程中，需要避開骨骼，肺部氣管支氣管樹分支和較大的肺血管，否則會引起氣胸、血胸、血氣胸、肺出血、大咯血、肺部感染、病灶播散，嚴重時甚至導致患者死亡[6-7]。數據調查表明，氣胸和肺出血是經皮肺穿刺手術的主要并發癥，其中氣胸的發生率為10%～40%，肺出血的發生率為26%～33%[8]。為了降低并發癥的發生率，根據腫瘤病灶的位置設計合理的穿刺路徑是肺穿刺領域的重要研究內容之一[9]。
目前，確定經皮肺穿刺路徑主要依靠醫生借助電子計算機斷層掃描（computed tomography，CT）圖像憑經驗手工完成。隨著計算機圖像處理和可視化技術的發展，通過CT圖像獲得最優的安全穿刺路徑成為了可能。穿刺手術路徑規劃算法的研究對象主要分為柔性針和剛性針兩類。柔性針路徑規劃算法借助柔性針本身可以改變進針方向的特點，試圖尋找復雜的曲線路徑以避開障礙，主要包括逆運動學反解法、人工勢場法和數值法三類[10]。柔性針路徑規劃算法的計算較為復雜，在實際穿刺過程中也需要構建符合實際情況的針—組織交互模型，同時考慮到柔性針本身的特點，目前仍處于實驗驗證階段。
相比于柔性針，剛性針路徑規劃算法的穿刺路徑為直線，是目前人工穿刺和機器人穿刺的主要手段。Kimura等[11]利用多平面重建（multiplanar reconstruction，MPR）技術重建患者CT數據集的特殊截面，引導穿刺針對處于復雜區域的腫瘤病灶進行斜向穿刺，但是穿刺路徑的規劃依舊需要人工完成。為了實現全自動剛性針路徑規劃， Villard等[12-13]、Baegert等[14-16]和Seitel等[17]以患者器官與組織的面繪制模型為基礎，針對臨床約束條件分別利用三角形網格曲面繪制方法對可行進針區域進行篩選與優化；Schumann等[18-21]以腫瘤病灶為原點，通過圓柱投影對不同優化目標生成二維約束地圖，采取加權乘積、手工交互或帕累托優化（Pareto optimization）的方法選擇穿刺路徑；張睿等[22]基于改進的立方體映射算法，在患者CT圖像的包圍盒表面進行二值化分類，得到滿足臨床約束條件的可行穿刺區域，再針對不同的臨床穿刺準則進行加權整合，提供若干穿刺路徑供醫生選擇。但是，上述剛性針路徑規劃算法往往會將數字化人體三維模型投影到某個特殊曲面，這個曲面可能是CT數據集的特殊截面、皮膚模型、圓柱面或立方體面，再基于曲面上的信息規劃穿刺路徑。其本質都是對原三維模型復雜空間關系的映射過程，不可避免地會發生信息失真，使投影到曲面的圖像發生畸變，難以確定最優穿刺路徑。另一方面，Gao等[23]以腫瘤病灶為中心，通過構建球面采樣點的方法直接在三維模型上進行路徑篩選，再根據入刺深度、目標區域可視大小等標準對穿刺路徑進行評分，避免了三維模型向曲面映射時所可能產生的信息失真問題。但是作者并未說明球面采樣點的具體生成方式，也未討論如何確定采樣點的數目，同時，不完整的皮膚模型使“入刺深度”這一評價標準僅適用于部分穿刺路徑，存在一定的缺陷。
針對以上問題，本文提出了一種基于多級約束的肺部穿刺路徑規劃方法。在患者數字化人體三維模型的基礎上，首先以腫瘤病灶為中心建立理想球面并在其上利用斐波那契網格（Fibonacci lattice）生成采樣點，其中每一個采樣點到病灶中心均構成一條可能的穿刺路徑；然后根據臨床穿刺準則對上述路徑進行約束，主要分為物理約束、生理約束和數理約束三部分：先根據體位、工具、設備等物理條件的限制排除大部分路徑；再利用方向包圍盒層次樹（oriented bounding box tree，OBBTree）模擬穿刺過程以規避風險組織和器官，篩選出滿足生理約束條件的穿刺路徑；最后在數理約束階段根據臨床穿刺準則對篩選出的穿刺路徑進行量化評價，并利用帕累托優化算法得到最優穿刺路徑。
1 方法
本文所提算法的流程如圖1所示。在胸部數字化模型上以腫瘤病灶為中心，穿刺針長度為半徑構建理想球面，通過斐波那契網格采樣得到離散點點集S0，其中每一點與病灶中心的連線均構成一條待選路徑；首先根據體位、工具、設備對穿刺路徑的限制完成物理約束篩選，得到路徑集合S1；再采用OBBTree算法模擬穿刺過程，篩選出滿足生理約束的穿刺路徑集合S2；最后采用帕累托優化算法對滿足物理和生理約束的穿刺路徑進行量化評價，得到備選路徑集合S3，選擇其中評分最高的路徑作為最優穿刺路徑s。
圖1
基于多級約束的肺部穿刺路徑規劃算法流程圖
Figure1.
Flow chart of the algorithm for lung puncture path planning based on the multi-level constraint
1.1 斐波那契網格采樣
斐波那契網格是由解析公式生成的一組點集，其中的點沿著密繞螺旋排列，在螺旋軸線方向上，每個點與上一點之間的角度差均構成黃金分割比，從宏觀上看，這些點所構成的點集在球面或圓平面上以基本各向同性的方式均勻分布[24]。利用其在球面或圓平面上點集分布的均勻性質，斐波那契網格已經在球面不規則圖形面積測量[24]、顏色量化與處理[25]以及相控陣天線布板設計[26]等方面得到應用。因此可以利用斐波那契網格對球面進行離散化采樣，以此確定待測穿刺路徑集合，記為點集M。如式(1)所示，設球面半徑為R，點集M中的元素總數為N，對于點集M中的任一點m(m = {(xn, yn, zn), n ∈ [1, N]})有：
|
其中，φ為定值，取，n為1～N范圍內的整數，(xn, yn, zn)表示由n確定的m點坐標。
斐波那契網格的構成需要理想球面的半徑R與元素總數N作為輸入。根據穿刺針長度可以確定半徑R的大小，但在穿刺路徑的計算中往往并不關注元素總數N的多少，而更關注穿刺路徑分布的疏密程度，即點集中點與點之間的最近距離大小，記為分辨率δ。因此就需要得到分辨率δ與元素總數N的關系，以便在給出δ后可以直接估算出所需的N值。假設以各個點為圓心，δ/2為半徑在球面上構成一系列圓盤，在理想的絕對均勻條件下，這一系列圓盤在數學上形成了對球體的最密封裝，對于元素總數N的估計如式(2)所示：
|
其中，SM表示所有圓盤的總面積；Sm表示單個圓盤的面積；S表示整個球面的面積。由于S略大于SM，故利用此公式計算得出的N也略大于理想的元素總數，因此該公式可以在保證分辨率要求的情況下，用于對元素總數N的估計。
以上斐波那契網格采樣所生成的是均勻分布于整個球面上的點集，從點集中的點m到穿刺目標點所形成的路徑中，有的并不能滿足臨床穿刺的要求。因此就需要對斐波那契網格采樣所產生的點集進行多級約束篩選，最終得到最優穿刺路徑。
1.2 基于臨床穿刺準則的多級約束設計
臨床穿刺準則是穿刺醫生在實踐中總結出的經驗標準，其中規定了穿刺操作的限制條件以及穿刺所要達到的理想目標。使用計算機技術進行肺部穿刺路徑規劃，需要對臨床穿刺準則進行量化分析。本文將臨床穿刺準則歸納為三個部分：物理約束、生理約束和數理約束。
1.2.1 物理約束
物理約束是指在肺部穿刺時由于體位、工具、設備的原因對穿刺路徑造成的限制，在規劃路徑時必須予以規避。雖然理論上醫生可以在任意方向實施任意距離的穿刺，但是由于患者在穿刺設備上所采取的體位、穿刺針長度和穿刺角度等因素的影響，客觀上限制了穿刺路徑的選擇。物理約束主要可分為三類：
（1）患者在穿刺設備上所采取的體位。CT引導經皮肺穿刺的常用體位包括仰臥位、俯臥位、左側臥位及右側臥位[27]。通過對腫瘤病灶位置的分析，醫生可以初步確定采用何種體位，由于設備的遮擋關系，不同體位所允許的穿刺角度區間也不同。以仰臥位為例，可以將胸部前方區域作為合適的穿刺區域。
（2）穿刺針本身的長度限制。在肺穿刺手術中選擇路徑的長度一定不會超過穿刺針的長度。假設穿刺針穿刺時可達的最大長度為R（考慮到醫生在臨床操作時的實際情況，保險起見R會比穿刺針的真實長度略短），即不必考慮超過R范圍的穿刺路徑，因為它們超過了穿刺針的最大作業范圍，無法成功完成穿刺。如式(3)所示，假設皮膚表面上任意點P(x, y, z)與選定的穿刺目標點O(x0, y0, z0)之間的總長度為d，則有：
|
當d ≤ R時，可以認為從該點到穿刺目標點對應的路徑滿足了穿刺針長度約束，予以保留，否則剔除該點。對于滿足穿刺針長度約束的穿刺路徑，從皮膚到腫瘤病灶的長度越小，對患者造成的損傷越小，穿刺風險越低。在斐波那契網格采樣過程中，已經以腫瘤病灶為中心，穿刺針長度為半徑構建了理想球面，其上每一采樣點到病灶中心的距離均為R。結合數字化人體模型，只有當采樣點位于患者皮膚外時，對應路徑從入刺點到病灶中心的距離d在R范圍內，予以保留，否則剔除該點。
（3）穿刺針與皮膚的夾角限制。在穿刺針與患者皮膚作用的過程中，兩者之間的夾角需要大于臨床常規閾值（20°）以避免滑針[28]。如式(4)所示，設皮膚表面上任意點P有單位法向量N(i, j, k)，選定的穿刺目標點與備選點之間構成單位方向向量Z(a, b, c)，可以求得穿刺針進針角度θ為：
|
根據臨床常規閾值，當θ ≥ 20°時，可以認為穿刺針穿刺該點時不會滑針，予以保留，否則剔除該點。對于滿足穿刺針角度約束的穿刺路徑，它的方向向量與皮膚法向量之間的夾角越小，滑針的概率越小，穿刺風險越低。
穿刺針長度約束和角度約束示意圖如圖2所示，紅色區域為腫瘤病灶，其中心為穿刺目標點O，橙色線條為皮膚的橫截面。假設經過斐波那契網格采樣產生三條路徑，其中路徑1（綠色線條）和路徑2（紅色線條）的穿刺深度分別為d1和d2，穿刺角度為θ1和θ2，可以看出路徑1更好地滿足了物理約束條件；路徑3（藍色線條）由于長度不足，未能穿出皮膚，需要予以剔除。
圖2
穿刺針長度約束與角度約束
Figure2.
Length constraint and angle constraint of puncture needle
除了上述總結的約束類型，為了方便醫生更好地操作工具和設備、觀察穿刺路徑，還可以針對具體手術環境設置額外的物理約束。
1.2.2 生理約束
生理約束是指為了大幅度降低穿刺時氣胸、肺出血等并發癥的發生概率而對穿刺路徑進行的限制。這就需要穿刺針在到達腫瘤病灶的過程中規避骨骼、肺部大氣管、大血管、心臟等風險組織和器官。在重建的數字化人體模型基礎上，可以將骨骼、氣管、血管、心臟、肝臟等組織和器官視為穿刺路徑所不能通過的障礙，再進行生理約束的計算。如式(5)所示，假設點集中有任意點P(x, y, z)，選定的穿刺目標點為O(x0, y0, z0)，可以確定穿刺路徑的單位方向向量Z(a, b, c)為：
|
其中，d為點P與點O之間的總長度。從穿刺目標點O(x0, y0, z0)沿著向量Z的方向進行路徑搜索，如果該路徑被重要器官或組織遮擋則剔除該路徑，否則繼續沿著原方向進行搜索，如果可以到達備選點P(x, y, z)，則證明該路徑滿足重要器官、組織避障約束，最終予以保留。對于滿足生理約束的穿刺路徑，它與最近障礙之間的直線距離越大，成功規避風險組織和器官的概率越大，穿刺風險越低。
1.2.3 數理約束
數理約束是在滿足前兩級約束的基礎上，為了進一步降低穿刺風險，提高穿刺可靠性，運用數學工具從可能的穿刺路徑中篩選最優路徑的過程。主要內容包括：①依照臨床穿刺準則對穿刺路徑進行量化評價。②根據帕累托優化結果得出最優穿刺路徑。
經過物理約束和生理約束篩選的穿刺路徑在理論上均滿足穿刺手術的基本要求，但是不同路徑對約束條件的滿足程度各不相同。有的路徑雖然滿足了物理、生理約束，但是距離風險組織和器官依然很近，考慮到實際穿刺時的操作誤差，此類路徑的手術風險過大，需要予以剔除。因此，為了能夠選出最優穿刺路徑，需要依照臨床穿刺準則對穿刺路徑進行量化評價。
根據物理約束和生理約束，定義從皮膚入刺點到病灶的路徑總長L、路徑與皮膚法向量夾角角度Θ和路徑與風險組織器官的最近距離D作為穿刺路徑的評價標準。其中L和Θ的大小可以根據穿刺路徑與皮膚的空間關系直接得出，為了消除量綱差異，使用離差標準化方法對上述兩個標準進行歸一化處理，如式(6)所示，對路徑i有：
|
其中，Li、Θi是當前路徑的長度和夾角，Lmin、Θmin是所有有效路徑長度和夾角的最小值，Lmax、Θmax是所有有效路徑長度和夾角的最大值。由定義可知，L(i)和Θ(i)的值越小，對應路徑的評價越高，其風險越小。
而求解各路徑的D值涉及到穿刺路徑與風險組織器官之間的復雜空間關系，很難進行直接的計算，因此需要采用其他方法間接地對D值進行量化。設通過生理約束篩選的點集為A，未能通過生理約束篩選的點集為 。相比于直接計算穿刺路徑的D值，找出點集A中與點集 的距離在δ以內的點更加容易。因此，可以通過這種方式得到點集A所代表的穿刺路徑中D值小于δ的部分，記為A1。以此類推，將A1合并到 中，重復上述計算步驟，可以得到A2、A3、A4、，在這些集合中穿刺路徑的D值越來越大，風險越來越低。對點集A1、A2、A3、A4、…所代表的穿刺路徑分別賦予D值1、2、3、4、，即完成了對D值的量化。由定義可知，D(i)的值越大，對應路徑的風險越小。
在對每一條穿刺路徑進行量化評價后，需要依據評分高低選擇出最優穿刺路徑。但是不同的評價標準往往存在相互制約的關系，在某一標準上表現較好的穿刺路徑，在其他標準上可能具有較差的表現。為了解決上述問題，可以為不同的評價標準分別設置權重系數，將經由權重系數調整的評分加和作為穿刺路徑的總體評價。這種基于權重的評價方法可以直接選出最優穿刺路徑，但是需要人為對權重系數進行設置，而權重系數的設置并沒有一個明確的標準，這為穿刺路徑的規劃帶來了主觀因素的影響。
為了避免上述方法的缺點，可以采用帕累托優化的方法得出最優穿刺路徑[29]，通過此方法可以在不額外設置系數的前提下求解得到位于帕累托前沿面上的解，這些位于帕累托前沿面上的解之間無法進一步區分優劣，而在多個評價標準中均優于非帕累托前沿面上的解。
2 實驗與結果
2.1 實驗數據
肺部穿刺路徑規劃算法必須以患者胸部的數字化模型為基礎，該模型應包括肺穿刺需要規避的重要組織和器官，數字化模型的精度直接影響到穿刺路徑規劃的效果。目前重建胸部數字化模型的素材主要來自于患者的CT數據集，隨著計算機技術，特別是神經網絡的發展，針對CT數據集的自動分割算法已經能夠對包括肺部精細化結構在內的胸部組織器官進行高效的三維重建，可以滿足本文所提算法的需求，在獲取CT影像后及時構建模體并進行手術[30-32]。本文研究重點在于肺部穿刺路徑規劃算法，因此依然采用人工分割的方式建立數字化人體模型。本實驗所用數據來源于美國國家癌癥研究所發布的公開數據集肺部圖像數據庫聯盟與圖像數據庫資源計劃（lung image database consortium and image database resource initiative，LIDC-IDRI）（網址為：https://wiki.cancerimagingarchive.net/display/Public/LIDC-IDRI）[33]。數字化人體模型的建立過程如圖3所示，采用CT圖像的像素尺寸為0.75 mm × 0.75 mm，層厚1 mm，切片數341層，層分辨率為512 × 512像素，利用醫學圖像處理軟件3Dslicer4.11（The Slicer Community，美國）重建了包括患者皮膚、骨骼、肺（肺實質、肺氣管、肺血管）、心臟和部分肝在內的三維模型，最后裝配成為患者胸部數字化模型，數字化模型中的黑色區域為模擬腫瘤的位置[34]。
圖3
數字化人體模型的建立過程
Figure3.
The process of building the digital human body model
另外，在原有CT圖像坐標系O’-ijk中，CT圖片左上角為坐標原點，向右為i軸正方向，向下為j軸正方向，向后為k軸正方向。為了更方便地表示穿刺路徑并減少計算量，需要建立以病灶中心O(x0, y0, z0)為原點的肺穿刺路徑規劃坐標系O-xyz，兩個坐標系可以互相轉換。
2.2 斐波那契網格采樣的創建
以患者胸部數字化模型為基礎，在肺穿刺路徑規劃坐標系O-xyz中進行斐波那契網格采樣。為此需要確定三個參數：理想球面的球心，理想球面的半徑R以及分辨率δ。理想球面的球心就是腫瘤病灶的中心O(x0, y0, z0)，其半徑R可由穿刺針長度確定，穿刺精度決定了分辨率δ的大小。以人工穿刺為例，在CT引導下經皮肺穿刺活檢中使用的穿刺針長度約為21 cm[28]，針尖偏移距離誤差約為5 mm[35]。根據以上數據，有R = 210 mm，δ = 5 mm。
2.3 基于多級約束的肺部穿刺路徑規劃
2.3.1 物理約束篩選
物理約束可以進一步縮小穿刺路徑的可行區域，其結果如圖4所示，相比于整個理想球面，由綠色區域代表的穿刺針可行區域大大減少，在該區域內的斐波那契網格采樣點到穿刺目標點之間的連線均構成一條待選穿刺路徑。
圖4
滿足物理約束的路徑集合
Figure4.
The set of paths satisfying the physical constraints
2.3.2 生理約束篩選
物理約束篩選的結果構成一離散化點集，以其中的點為一端，穿刺目標點為另一端模擬穿刺過程，利用OBBTree算法篩選出滿足生理約束的穿刺路徑。OBBTree是一種多邊形有向包圍盒的二叉樹數據結構和算法，用于物體之間的快速相交測試[36]。利用基于OBBTree的快速相交測試算法，通過判斷路徑和重要器官或組織是否存在交點，可以選擇出規避風險區域的穿刺路徑。如圖5所示，由綠色區域代表的穿刺針可行區域進一步減少，經過放大可以看到，點集輪廓（黑色邊緣）能夠完全反映風險區域對腫瘤病灶的遮擋范圍（紅色邊緣），完成了對風險區域的規避。
圖5
滿足生理約束的路徑集合
Figure5.
The set of paths satisfying the physiological constraints
2.3.3 數理約束篩選
對于從生理約束中篩選出的大量路徑，需要通過數理約束對其進行量化評價，得出最優穿刺路徑。在各項評價標準中，從皮膚入刺點到病灶的路徑總長L和路徑與皮膚法向量夾角角度Θ比較容易計算得出；各路徑與風險組織器官的最近距離D如圖6所示，可以看到，隨著D值的增大，穿刺路徑與風險組織器官的距離同時增加，穿刺風險降低。
圖6
各路徑與風險組織器官的最近距離
Figure6.
Nearest distance of each pathway to the risk tissue and organ
考慮到D值較小的路徑過于靠近風險組織器官，因此選擇在D值較大的路徑中采用帕累托優化方法篩選最優穿刺路徑。如圖7所示，x軸表示從皮膚入刺點到病灶的路徑總長評分L(i)，y軸表示路徑與皮膚法向量夾角角度評分Θ(i)，兩者均由式(6)計算得出。針對不同的D值得到的帕累托前沿面顯示為紅線，在帕累托前沿面上的點可以視為該D值下的最優點。可以看到，在帕累托前沿面p1、p2和p3中，p1上的穿刺路徑具有最短的穿刺距離和最好的穿刺角度，但是距離風險組織器官較近；相反，p3上的穿刺路徑距離風險組織器官最遠，但穿刺距離和穿刺角度都稍遜于p1和p2上的穿刺路徑。
圖7
不同D值下的帕累托前沿面
Figure7.
Pareto front surface for different D values
將處于各帕累托前沿面上的點所代表的穿刺路徑表示在患者數字化人體模型上，如圖8所示，位于p1上的2條穿刺路徑標記為淺藍色，位于p2上的2條穿刺路徑標記為深藍色，位于p3上的2條穿刺路徑標記為紅色，這六條路徑之間沒有明確的優劣關系，而在三個評價標準上均優于其他穿刺路徑，醫生可以在上述六條路徑中進行選擇，大大減少了主觀因素的影響。
圖8
路徑規劃的最終結果
Figure8.
The final result of path planning
2.4 實驗結果
為了驗證算法的穩定性，基于不同體型、年齡、性別患者的CT數據建立了另外9套數字化人體模型，分別為模型1～模型9，隨機設定病灶原點，并運用本文所提出的方法進行肺部穿刺路徑規劃。如圖9所示，最終篩選出的穿刺路徑標記為黑色，其數目在4～12條之間，平均在每個模型中可以計算得到7條穿刺路徑。在子圖模型2、模型3、模型4、模型6、模型8中，篩選出的穿刺路徑具有一定的聚集性，表示最優穿刺路徑即在此區域內；在子圖模型5、模型7、模型9中，篩選出的穿刺路徑比較分散，表示在已有的穿刺標準下有多個可能的穿刺方案；在子圖模型1中，篩選出的穿刺路徑總體具有分散的特點，而在個別區域又有一定的聚集性。實驗結果證明，該算法具有一定的穩定性，數量適中的穿刺路徑可以在幫助醫生確定穿刺范圍的同時，提供一定的選擇方案。
圖9
基于不同患者模型的路徑規劃結果
Figure9.
Path planning results based on different patient models
3 總結
本文提出了一種基于多級約束的肺部穿刺路徑規劃算法，在患者數字化人體模型的基礎上，首先在以腫瘤病灶為中心的理想球面上構造斐波那契網格采樣生成穿刺路徑，根據體位、工具、設備對穿刺造成的限制完成物理約束的篩選，利用OBBTree算法篩選滿足生理約束條件的穿刺路徑，最后利用帕累托優化算法得到了有代表性的肺部穿刺路徑。經過實驗驗證，該算法可以規避骨骼、氣管、血管等風險組織，能夠選擇出相對安全的穿刺路徑供醫生選擇，同時也為未來實現機器人穿刺提供了技術支撐。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：孫峰輝主要負責算法程序設計、實驗數據分析以及論文初稿撰寫；裴宏亮主要負責實驗方法設計和實驗結果可視化；楊逸飛主要負責數據整理與管理；樊慶文主要負責提供實驗指導，數據分析指導和論文審閱修訂；李曉歐主要負責提供醫學指導以及對文章的醫學內容作批評性審閱。
