體外沖擊波是治療骨關節炎的一種有效方法, 研究骨結構對體外沖擊波治療的影響以及進一步提高體外沖擊波治療療效都具有重要的應用價值。本文建立了沖擊波治療骨關節炎的有限元模型, 并制作人工骨結構及高速攝像系統對模型加以驗證, 證明所建立模型的有效性。最后利用所建的模型, 預測了沖擊波治療骨關節時的能量分布。結果表明骨結構的存在會使得沖擊波治療時能量產生折射, 且折射角與骨結構的位置相關。有限元方法結合磁共振/計算機斷層掃描成像可以更好地保障體外沖擊波治療骨關節炎的治療劑量及精度。
引用本文： 張仕年, 王蕭峰, 章東. 骨關節炎體外沖擊波治療的能量分布研究. 生物醫學工程學雜志, 2015, 32(2): 300-304. doi: 10.7507/1001-5515.20150055 復制
引言
骨關節疾病包括退行性關節炎、滑囊炎、滑膜炎、頸椎病、腰椎病、肩周炎、骨質增生、風濕性關節炎、類風濕性關節炎、股骨頭壞死等，多數屬于慢性病，藥物治療時間長，見效慢，有些甚至無效，容易產生耐藥性。傳統治療技術(如磁療、熱療、推拿等)只能緩解癥狀，不能從根本上解決問題。手術方法對患者有損傷，且費時費力[1-2]。而體外沖擊波治療則可以有效解決這些難題，通過直接刺激成骨細胞和間質細胞，使其快速修復，快速消除炎癥、緩解疼痛，具有無手術損傷且短期內達到治療目的的優勢[3-5]。
在體外沖擊波治療骨關節炎過程中，如果沖擊波的聲壓足夠大，體液中會產生空化泡，并在沖擊波的負壓相中快速生長，然后劇烈崩潰，這種現象稱為穩態空化[6]。很多研究者指出空化效應是體外沖擊波治療的主要機制[7-8]。盡管體外沖擊波治療技術從80年代起就在臨床上得以應用，但研究者僅建立了簡單的數值模型以研究沖擊波在水及組織中的傳播[9]。這些簡單模型沒有考慮到骨關節炎治療中骨結構的存在，而骨結構的存在會改變沖擊波的治療焦域及能量分布，甚至會激發復雜的模式轉換，影響治療精度。因此更好地理解沖擊波治療骨關節炎時的能量分布對提高其治療療效有積極意義。
本文基于有限元商用軟件，建立了體外沖擊波治療骨關節炎的有限元模型。為驗證有限元模型的有效性，采用高速相機實驗記錄體外沖擊波治療儀治療人工骨結構時產生的空化泡，并與仿真結果對比。最后，基于有限元模型計算了沖擊波治療骨關節炎時的能量分布。
1 有限元模型及實驗方法
1.1 有限元模型
本文沖擊波治療骨關節炎的有限元模型是基于線性聲傳播理論[10-11]。沖擊波脈沖的形成是由于聲壓由靜壓突然上升到最高壓，波速由于聲壓的上升而增加，波形畸變形成沖擊波形。如圖 1所示，沖擊波治療器中的火花塞位于半橢園焦點F1以產生沖擊波脈沖。骨結構處于橢圓的第二個焦點F2，圖中T表示骨結構的表面。焦點F1處的沖擊波可以用如下公式表示，即
| ${P_s}\left(t \right)=2{P_a}{{\rm{e}}^{-\alpha t}}\cos \left({2\pi ft + \pi /3} \right), $ |
其中α=3.5×105，f=50 kHz，Pa=50 MPa。
圖1
體外沖擊波治療骨關節炎圖示
Figure1.
Schematic illustration of the electrohydraulic lithotripsy system
沖擊波在水、組織及骨結構中的傳播采用如下聲波方程：
| $\rho \frac{{{\partial ^2}\boldsymbol{v}}}{{\partial {\boldsymbol{t}^2}}}=\left({\lambda + 2\mu } \right)\nabla \left({\nabla \cdot \boldsymbol{v}} \right)-\mu \nabla \times \left({\nabla \times \boldsymbol{v}} \right), $ |
其中ρ為傳播介質的密度，v為聲速，λ、μ為拉密系數。骨結構中的彈性聲波方程為
| $\begin{array}{l} {\sigma _{xx}}=\lambda \left({{\varepsilon _{xx}} + {\varepsilon _{yy}} + {\varepsilon _{zz}}} \right) + 2\mu {\varepsilon _{xx}}\\ {\sigma _{yy}}=\lambda \left({{\varepsilon _{xx}} + {\varepsilon _{yy}} + {\varepsilon _{zz}}} \right) + 2\mu {\varepsilon _{yy}}\\ {\sigma _{zz}}=\lambda \left({{\varepsilon _{xx}} + {\varepsilon _{yy}} + {\varepsilon _{zz}}} \right) + 2\mu {\varepsilon _{zz}}, \\ {\sigma _{yz}}=\mu {\varepsilon _{yz}}\\ {\sigma _{zx}}=\mu {\varepsilon _{zx}}\\ {\sigma _{xy}}=\mu {\varepsilon _{xy}} \end{array}$ |
其中σij(i, j=x, y, z)為彈性矢量，且εij(i, j=x, y, z)=為彈性張量，位移uij(i, j=x, y, z)與聲速相關。骨表面分布可由CT掃描圖像獲得。數值模擬中水及骨結構參數如下：ρwater=1 000 kg/m3, λwater=2 190.4 MPa, μwater=0 MPa, ρABS=1 210 kg/m3, λABS=3 264 MPa, μABS=1 032 MPa。有限元網格的設定與材料的ρ、λ和μ相關。本文采用動態網格以優化沖擊波的傳播波形。有限元模擬中對橢圓反射器的邊界采用硬邊界條件，而其它水邊界采用吸收邊界。在水與骨結構連接處需滿足垂直速度與應力張量連續，即vtissue, n=vbone, n以及σtissue=σbone。
1.2 實驗方法
實驗中基于CT掃描的圖像，采用人工骨制作樣品，置于一個可調節位置的臺上并浸沒于水中。類似于臨床治療情況，人工骨與沖擊波傳播方向相對。沖擊波治療儀(Dornier HM-3 Lithotripter, 美國)，工作電壓為22 kV，采用2 Hz脈沖重復頻率。數值仿真及實驗中為考察骨結構對沖擊波治療的影響，以3 mm為步長，改變骨結構表面與F2焦點的距離。為建立穩態的空化場，每次實驗前，火花塞需連續20次點火。空化場采用高速相機(Imacon 200, 美國)拍攝，該相機幀率達到2×107幀/s，每次曝光可得到14幅高分辨率圖像，幀間隔為50 ns~20 ms。圖像采集后，采用MatLab軟件進行分析。
2 結果與討論
2.1 數值仿真結果
為方便與實驗結果驗證，在有限元模擬中采用二維(xz平面)仿真。圖 2給出了沿z方向移動骨結構時彈性張量的折射變化情況。圖中暗區表示沖擊波傳播1.5μs時150幅圖像累積的彈性張量的最大值，表示了空化泡出現的區域。當骨表面遠離F2焦點(即z=-12 mm)時，沖擊波的傳播是自由場的情況；而當骨結構接近沖擊波治療器的焦點后，沖擊波聲場的折射角變大。
圖2
靶點T與焦點F2之間的距離變化時，有限元方法得到的最大彈性張量
圖中星號表示焦點
the asterisk and solid point represent the shock focus
2.2 實驗結果
為驗證有限元仿真的正確性，本文采用高速相機拍攝沖擊波治療時產生的空化泡，圖 3為高速相機得到的系列圖像。由圖可見骨結構位置對沖擊波聲場的影響。每幅圖像中從左邊開始出現的空化泡表示沖擊波的傳播路徑。當骨結構與聲焦點的相對位置由z=-12 mm改變至9 mm時，空化泡群區域的形狀出現彎曲向上，偏離骨結構表面。這種折曲現象與數值仿真結果相符，表明當沖擊波接近骨表面時聲能量向上產生折射。
圖3
靶點T與焦點F2之間的距離變化時，高速相機系統得到的空化泡分布
(a)-12 mm; (b)-9 mm; (c)-6 mm; (d)-3 mm; (e) 0 mm; (f) 3 mm; (g) 6 mm; (h) 9 mm
Figure3. Cavitation bubbles using high-speed imaging technology for various vertical distances between F2 and T point(a)-12 mm; (b)-9 mm; (c)-6 mm; (d)-3 mm; (e) 0 mm; (f) 3 mm; (g) 6 mm; (h) 9 mm
假設骨結構由線性彈性材料構成，其彈性勢能, 其中E為楊氏模量，且有及。根據上面的有限元模型，可以計算沖擊波治療時的能量分布。圖 4(a)給出了當靶點T位于F2的-6 mm處時骨中的能量分布。由圖可見，面對沖擊波源的骨表面具有較高的能量，并且能量隨著在骨中透射深度的增加而發散。在靶點T與F2的垂直距離z變化時，圖 4(b)給出了以靶點T為圓心、半徑5 mm的圓上彈性勢能的積分隨時間變化曲線。由圖可見，自沖擊波治療器中的火花塞發射脈沖，經過約0.18 ms可使骨區域的聲能量達到峰值。當骨結構遠離F2時(z=-9 mm)，由于靶點T點與焦點偏離很多，傳播至骨結構的聲能量很小；當骨結構沿著z方向移動，靶點T接近焦點，越來越多的聲能量傳至骨結構；當靶點T與F2近似重合時(z=3 mm)，聲能量達到峰值。需要指出的是，當靶點T與F2完全重合時(即z=0 mm)，聲能量并非最大，這是由于骨結構表現的不規則性造成的，也因此說明根據有限元數值模擬可以更準確地了解聲能量的傳播途徑及其在骨結構中的分布情況，從而更好地指導臨床治療。
圖4
骨結構中能量的有限元模擬
(a)靶點
(a) calculated energy density distribution, as the vertical distance between
2.3 討論
本文從數值仿真及實驗兩方面研究沖擊波治療骨關節炎時的聲傳播。圖 2及圖 3得到的結果均表明沖擊波脈沖會受到骨結構的影響。當骨結構在z方向對于焦點F2向上移動時，沖擊波會向上折射且折射角會改變。
盡管數值模擬與實驗結果的比較驗證了有限元模型能較好地預測沖擊波的傳播及骨結構存在時產生的折射，但從圖 2與圖 3的比較可看出，測量結果比計算結果更明顯，這是因為：①有限元計算中網格一般要小于1/6波長，但網格變小需極大地增加計算時間；②本文的有限元模型是基于線性聲學假設，但沖擊波治療中高壓引起的強烈的非線性效應必須考慮；③此外，周圍組織的吸收、非均勻性等的考慮將會進一步完善有限元模型。
3 結束語
體外沖擊波治療骨關節炎的主要機制是空化效應，利用此效應能疏通微血管、松解關節和軟組織粘連; 并且沖擊波的高能量也可緩解疼痛。該作用機制表明了沖擊波能夠在骨科廣泛應用，對骨肌系統的很多疾病都有較好的療效。但如何定位體外沖擊波治療的沖擊點以及治療劑量的選擇是取得良好療效的關鍵。本文的主要目的在于研究骨結構的存在對沖擊波聲場的影響。在線性聲波傳播理論基礎上建立了有限元模型，數值仿真結果表明骨結構的存在會使得空化微泡場折射并且折射角會隨著目標點與焦點的位置而改變。高速相機系統的實驗證實了該有限元模型的有效性。并且該有限元模型可以有效地預測骨表面的沖擊波能量分布，可以更為有效地治療骨肌損傷。
引言
骨關節疾病包括退行性關節炎、滑囊炎、滑膜炎、頸椎病、腰椎病、肩周炎、骨質增生、風濕性關節炎、類風濕性關節炎、股骨頭壞死等，多數屬于慢性病，藥物治療時間長，見效慢，有些甚至無效，容易產生耐藥性。傳統治療技術(如磁療、熱療、推拿等)只能緩解癥狀，不能從根本上解決問題。手術方法對患者有損傷，且費時費力[1-2]。而體外沖擊波治療則可以有效解決這些難題，通過直接刺激成骨細胞和間質細胞，使其快速修復，快速消除炎癥、緩解疼痛，具有無手術損傷且短期內達到治療目的的優勢[3-5]。
在體外沖擊波治療骨關節炎過程中，如果沖擊波的聲壓足夠大，體液中會產生空化泡，并在沖擊波的負壓相中快速生長，然后劇烈崩潰，這種現象稱為穩態空化[6]。很多研究者指出空化效應是體外沖擊波治療的主要機制[7-8]。盡管體外沖擊波治療技術從80年代起就在臨床上得以應用，但研究者僅建立了簡單的數值模型以研究沖擊波在水及組織中的傳播[9]。這些簡單模型沒有考慮到骨關節炎治療中骨結構的存在，而骨結構的存在會改變沖擊波的治療焦域及能量分布，甚至會激發復雜的模式轉換，影響治療精度。因此更好地理解沖擊波治療骨關節炎時的能量分布對提高其治療療效有積極意義。
本文基于有限元商用軟件，建立了體外沖擊波治療骨關節炎的有限元模型。為驗證有限元模型的有效性，采用高速相機實驗記錄體外沖擊波治療儀治療人工骨結構時產生的空化泡，并與仿真結果對比。最后，基于有限元模型計算了沖擊波治療骨關節炎時的能量分布。
1 有限元模型及實驗方法
1.1 有限元模型
本文沖擊波治療骨關節炎的有限元模型是基于線性聲傳播理論[10-11]。沖擊波脈沖的形成是由于聲壓由靜壓突然上升到最高壓，波速由于聲壓的上升而增加，波形畸變形成沖擊波形。如圖 1所示，沖擊波治療器中的火花塞位于半橢園焦點F1以產生沖擊波脈沖。骨結構處于橢圓的第二個焦點F2，圖中T表示骨結構的表面。焦點F1處的沖擊波可以用如下公式表示，即
| ${P_s}\left(t \right)=2{P_a}{{\rm{e}}^{-\alpha t}}\cos \left({2\pi ft + \pi /3} \right), $ |
其中α=3.5×105，f=50 kHz，Pa=50 MPa。
圖1
體外沖擊波治療骨關節炎圖示
Figure1.
Schematic illustration of the electrohydraulic lithotripsy system
沖擊波在水、組織及骨結構中的傳播采用如下聲波方程：
| $\rho \frac{{{\partial ^2}\boldsymbol{v}}}{{\partial {\boldsymbol{t}^2}}}=\left({\lambda + 2\mu } \right)\nabla \left({\nabla \cdot \boldsymbol{v}} \right)-\mu \nabla \times \left({\nabla \times \boldsymbol{v}} \right), $ |
其中ρ為傳播介質的密度，v為聲速，λ、μ為拉密系數。骨結構中的彈性聲波方程為
| $\begin{array}{l} {\sigma _{xx}}=\lambda \left({{\varepsilon _{xx}} + {\varepsilon _{yy}} + {\varepsilon _{zz}}} \right) + 2\mu {\varepsilon _{xx}}\\ {\sigma _{yy}}=\lambda \left({{\varepsilon _{xx}} + {\varepsilon _{yy}} + {\varepsilon _{zz}}} \right) + 2\mu {\varepsilon _{yy}}\\ {\sigma _{zz}}=\lambda \left({{\varepsilon _{xx}} + {\varepsilon _{yy}} + {\varepsilon _{zz}}} \right) + 2\mu {\varepsilon _{zz}}, \\ {\sigma _{yz}}=\mu {\varepsilon _{yz}}\\ {\sigma _{zx}}=\mu {\varepsilon _{zx}}\\ {\sigma _{xy}}=\mu {\varepsilon _{xy}} \end{array}$ |
其中σij(i, j=x, y, z)為彈性矢量，且εij(i, j=x, y, z)=為彈性張量，位移uij(i, j=x, y, z)與聲速相關。骨表面分布可由CT掃描圖像獲得。數值模擬中水及骨結構參數如下：ρwater=1 000 kg/m3, λwater=2 190.4 MPa, μwater=0 MPa, ρABS=1 210 kg/m3, λABS=3 264 MPa, μABS=1 032 MPa。有限元網格的設定與材料的ρ、λ和μ相關。本文采用動態網格以優化沖擊波的傳播波形。有限元模擬中對橢圓反射器的邊界采用硬邊界條件，而其它水邊界采用吸收邊界。在水與骨結構連接處需滿足垂直速度與應力張量連續，即vtissue, n=vbone, n以及σtissue=σbone。
1.2 實驗方法
實驗中基于CT掃描的圖像，采用人工骨制作樣品，置于一個可調節位置的臺上并浸沒于水中。類似于臨床治療情況，人工骨與沖擊波傳播方向相對。沖擊波治療儀(Dornier HM-3 Lithotripter, 美國)，工作電壓為22 kV，采用2 Hz脈沖重復頻率。數值仿真及實驗中為考察骨結構對沖擊波治療的影響，以3 mm為步長，改變骨結構表面與F2焦點的距離。為建立穩態的空化場，每次實驗前，火花塞需連續20次點火。空化場采用高速相機(Imacon 200, 美國)拍攝，該相機幀率達到2×107幀/s，每次曝光可得到14幅高分辨率圖像，幀間隔為50 ns~20 ms。圖像采集后，采用MatLab軟件進行分析。
2 結果與討論
2.1 數值仿真結果
為方便與實驗結果驗證，在有限元模擬中采用二維(xz平面)仿真。圖 2給出了沿z方向移動骨結構時彈性張量的折射變化情況。圖中暗區表示沖擊波傳播1.5μs時150幅圖像累積的彈性張量的最大值，表示了空化泡出現的區域。當骨表面遠離F2焦點(即z=-12 mm)時，沖擊波的傳播是自由場的情況；而當骨結構接近沖擊波治療器的焦點后，沖擊波聲場的折射角變大。
圖2
靶點T與焦點F2之間的距離變化時，有限元方法得到的最大彈性張量
圖中星號表示焦點
the asterisk and solid point represent the shock focus
2.2 實驗結果
為驗證有限元仿真的正確性，本文采用高速相機拍攝沖擊波治療時產生的空化泡，圖 3為高速相機得到的系列圖像。由圖可見骨結構位置對沖擊波聲場的影響。每幅圖像中從左邊開始出現的空化泡表示沖擊波的傳播路徑。當骨結構與聲焦點的相對位置由z=-12 mm改變至9 mm時，空化泡群區域的形狀出現彎曲向上，偏離骨結構表面。這種折曲現象與數值仿真結果相符，表明當沖擊波接近骨表面時聲能量向上產生折射。
圖3
靶點T與焦點F2之間的距離變化時，高速相機系統得到的空化泡分布
(a)-12 mm; (b)-9 mm; (c)-6 mm; (d)-3 mm; (e) 0 mm; (f) 3 mm; (g) 6 mm; (h) 9 mm
Figure3. Cavitation bubbles using high-speed imaging technology for various vertical distances between F2 and T point(a)-12 mm; (b)-9 mm; (c)-6 mm; (d)-3 mm; (e) 0 mm; (f) 3 mm; (g) 6 mm; (h) 9 mm
假設骨結構由線性彈性材料構成，其彈性勢能, 其中E為楊氏模量，且有及。根據上面的有限元模型，可以計算沖擊波治療時的能量分布。圖 4(a)給出了當靶點T位于F2的-6 mm處時骨中的能量分布。由圖可見，面對沖擊波源的骨表面具有較高的能量，并且能量隨著在骨中透射深度的增加而發散。在靶點T與F2的垂直距離z變化時，圖 4(b)給出了以靶點T為圓心、半徑5 mm的圓上彈性勢能的積分隨時間變化曲線。由圖可見，自沖擊波治療器中的火花塞發射脈沖，經過約0.18 ms可使骨區域的聲能量達到峰值。當骨結構遠離F2時(z=-9 mm)，由于靶點T點與焦點偏離很多，傳播至骨結構的聲能量很小；當骨結構沿著z方向移動，靶點T接近焦點，越來越多的聲能量傳至骨結構；當靶點T與F2近似重合時(z=3 mm)，聲能量達到峰值。需要指出的是，當靶點T與F2完全重合時(即z=0 mm)，聲能量并非最大，這是由于骨結構表現的不規則性造成的，也因此說明根據有限元數值模擬可以更準確地了解聲能量的傳播途徑及其在骨結構中的分布情況，從而更好地指導臨床治療。
圖4
骨結構中能量的有限元模擬
(a)靶點
(a) calculated energy density distribution, as the vertical distance between
2.3 討論
本文從數值仿真及實驗兩方面研究沖擊波治療骨關節炎時的聲傳播。圖 2及圖 3得到的結果均表明沖擊波脈沖會受到骨結構的影響。當骨結構在z方向對于焦點F2向上移動時，沖擊波會向上折射且折射角會改變。
盡管數值模擬與實驗結果的比較驗證了有限元模型能較好地預測沖擊波的傳播及骨結構存在時產生的折射，但從圖 2與圖 3的比較可看出，測量結果比計算結果更明顯，這是因為：①有限元計算中網格一般要小于1/6波長，但網格變小需極大地增加計算時間；②本文的有限元模型是基于線性聲學假設，但沖擊波治療中高壓引起的強烈的非線性效應必須考慮；③此外，周圍組織的吸收、非均勻性等的考慮將會進一步完善有限元模型。
3 結束語
體外沖擊波治療骨關節炎的主要機制是空化效應，利用此效應能疏通微血管、松解關節和軟組織粘連; 并且沖擊波的高能量也可緩解疼痛。該作用機制表明了沖擊波能夠在骨科廣泛應用，對骨肌系統的很多疾病都有較好的療效。但如何定位體外沖擊波治療的沖擊點以及治療劑量的選擇是取得良好療效的關鍵。本文的主要目的在于研究骨結構的存在對沖擊波聲場的影響。在線性聲波傳播理論基礎上建立了有限元模型，數值仿真結果表明骨結構的存在會使得空化微泡場折射并且折射角會隨著目標點與焦點的位置而改變。高速相機系統的實驗證實了該有限元模型的有效性。并且該有限元模型可以有效地預測骨表面的沖擊波能量分布，可以更為有效地治療骨肌損傷。
