用于皮膚病激光治療的脈沖光熱輻射模型

胡飛凡,徐晨辰,張浩,李東,陳斌,應朝霞

(1. 西安交通大學動力工程多相流國家重點實驗室,710049,西安; 2. 陜西省皮膚病防治所,725099,西安)

鮮紅斑痣(port wine stain, PWS)是一種常見的皮膚病,多發于面頸部等暴露位置[1],通常由皮膚表面/深層的血管擴張或異常增生引起,嚴重影響患者的外觀和身心健康?；谏锝M織的選擇性光熱作用[2],激光可以實現皮下病變血管的選擇性損傷并最大限度地保留正常組織,是目前臨床上血管性皮膚病治療的有效方法。

在皮膚病激光治療過程中,醫生首先對患者進行問診,觀察患處情況,通過對病情的評估選定治療激光類型,然后根據臨床經驗及醫學文獻中常用的激光參數,再結合相關理論知識,最終設定激光治療儀器的具體參數(脈寬、能量等)。目前,受限于皮膚組織結構的不確定性,僅憑醫生觀測進行參數選取,治療效果非常粗糙,常常會加重不良反應,如燙傷、術后疼痛、恢復時間延長及嚴重并發癥等,導致臨床完全治愈率已接近10年維持在低于20%的水平[3]。因此,迫切需要對皮膚組織進行內部溫度分布與結構的無損測量,為激光熱療參數的選取提供科學依據,有效指導醫生開展激光治療,從而提升治療效果。

近年來,脈沖光熱輻射法(pulsed photothermal radiometry, PPTR)[4-8]應用廣泛,作為一種非接觸式技術,它利用紅外熱像儀測量暴露于脈沖輻射(激光)下的樣品表面熱輻射信號(溫升信號)[4],結合有效的重建算法,反演被測皮膚樣品初始時刻的溫度分布和內部結構。PPTR法無需接觸患者皮膚,能避免對組織造成創傷,其測量信號由各個深度的溫度信號線性疊加至表面而成,其形式對于溫度反演來說為第一類Fredholm積分方程[9]。由于具備這種數學形式的問題往往是嚴重不適定的,這意味著重建的溫升分布通常不唯一,且對測量噪聲十分敏感。因此,設計反演算法時需要對問題采取有效的正則化策略。Mahdawi等采用V-循環多重網格法進行迭代,通過數值算例驗證了近似解的準確性和計算速度,得到更快的收斂性[10]。隨后,他們又改進了求解反問題的Landweber迭代法,并將其應用于求解熱傳導方程的反邊值問題[11]。與此同時,三維反演算法也由于計算機科學技術的飛速發展、計算性能的提升及實際工程問題的需求引起了廣泛關注,算法設計也是多種多樣。謝正超等采用基于廣義交叉驗證的Tikhonov正則化和截斷奇異值分解法,分別對煤粉鍋爐的爐內火焰進行了三維溫度場重建,發現 Tikhonov正則化的重建結果誤差較截斷奇異值分解法要小[12]。劉曉蔚等針對爆炸中溫度場的重建難以取得較高精度等問題,在聯合代數重建算法(SART)迭代系數中添加懲罰項,實現了更高質量的重建[13]。目前,脈沖光熱輻射法無損測量溫度分布的效果不斷提升,但鮮有學者將其應用于激光手術的溫度及結構探測輔助開展激光治療。

本文圍繞常見血管性皮膚病的激光治療這一主題開展研究,以鮮紅斑痣皮膚病為例,采用脈沖光熱輻射法構建皮膚傳熱理論模型,選取合適的正則化策略[14-15]開發溫升反演算法,實現了鮮紅斑痣患者皮下的溫升分布及結構重建,并開展體外實驗驗證其有效性,最后將其應用于臨床實驗。研究結果表明:在皮膚激光治療中,該應用模型可科學選取激光參數,減少醫生憑經驗選取參數導致的不良反應,能夠實現個性化治療,從而提升了治療的有效性。

1 數值模型構建

采用PPTR方法重建皮膚結構輔助激光治療的理論模型如圖1所示。首先,以紅外熱像儀作為表面輻射信號的測量工具,實現皮膚結構的非接觸測量,其原理為采用基于脈沖激光照射后皮膚表面溫度的衰減變化曲線,重建輻射照射結束瞬間的皮膚溫升分布,該溫升分布是激光輻射能在皮膚組織內傳輸、吸收以及皮膚組織導熱共同作用的結果[15]。根據生物組織選擇性光熱效應[2]理論,鮮紅斑痣血液層對不同波長的激光具有選擇性吸收效應,若選取的激光波長為鮮紅斑痣血液層吸收激光峰值所對應的波長,且激光作用時間小于血液層的熱弛豫時間,則該重建皮膚內的初始溫升分布可準確反映皮膚內血管層的結構信息(如血管層厚度、深度等)。同時,由于人體表皮層對激光同樣具有較強的吸收作用,因此該重建溫度也可反映表皮層的厚度等結構信息。隨后,將探測到的皮膚結構信息根據生物傳熱模型進行計算機建模、計算和優化,得到最優的激光治療參數來輔助醫生進行激光治療,從而提升了治療的精準性、安全性和高效性。

圖1 皮膚激光治療中的脈沖光熱輻射模型Fig.1 Model of pulsed photothermal radiation in dermal laser therapy

1.1 皮膚組織傳熱模型的建立

根據PWS的病癥特點[1],可知病變的毛細血管分布較為密集,如圖2(a)所示。假定真皮中雜亂分布的PWS血管在真皮內分布均勻,可把PWS病變區簡化為具有一定厚度的層狀區域,記為PWS層。這種做法對血管分布的處理具有一定的合理性,也可極大地降低生物組織內光和熱傳播計算的難度,提高了計算效率。此外,由于脈沖激光照射時間遠小于熱傳導所需時間,且皮膚激光手術的徑向尺度遠大于軸向深度尺度,因此可把該問題簡化成一個半無限大平板在半空間中、無內熱源的一維非穩態導熱問題。根據以上分析,可建立包含黑色素表皮層、真皮層、真皮層中PWS層在內的多層皮膚模型,如圖2(b)所示。

(a)鮮紅斑痣患者

將圖2所示模型假設為一個具有對流邊界的、常物性、非穩態、無內熱源的一維半無限大平板導熱問題,其溫度分布符合如下表達式

(1)

式中:ρ為密度;c為比熱容;z為模型的深度;t為時間;λ為熱導率;T(z,t)為溫度分布函數。

皮膚表面的Robin邊界條件可表示如下

(2)

式中:h為表面傳熱系數;T∞為皮膚周圍的環境溫度。

由于紅外熱像儀捕捉的是溫升信號,令熱擴散率a=λ/ρc,溫升分布函數ΔT(z,t)=T(z,t)-T∞,則式(1)、式(2)描述的定解問題可簡化為

(3)

(4)

式中:ν為熱損失系數,定義為h/λ,即表面傳熱系數與熱導率之比。

對于式(3)、式(4)組成的定解問題,可采用格林函數法[16]求解,以獲得溫升分布函數ΔT(z,t),數學解析表示如下

(5)

式中:ΔT(z′, 0)為初始溫升分布函數,用于描述處于位置z′的平面熱源;erfc(u)為互補誤差函數,可寫為

至此,通過建立一維模型和采用格林函數法求解正問題,獲得了圖2(b)所示多層皮膚組織內溫升分布函數的解析式,其積分項中包含的初始溫升分布函數即是待重建項。

1.2 導熱反問題的建立

由于PPTR方法以紅外熱像儀作為輻射信號的測量工具,因此需要基于紅外熱像儀的信號測量原理來建立PPTR信號幅度與初始溫升分布的關系。紅外熱像儀測量的表面信號幅度實際上是深度方向的皮下各層組織能量在皮膚表面的沉積,即在脈沖激光照射后,各深度組織輻射力的變化在皮膚表面的線性疊加。以深度為z的組織層為例,若其在激光照射下的溫升分布函數為ΔT(z,t),假設該溫升遠小于環境溫度,則深度為z處組織輻射力的變化量ΔEz滿足

(6)

式中:ε為實際物體的發射率(黑度);σ為黑體輻射常數,σ=5.67×10-8W·m-2·K-4。

從式(6)可以看出,輻射力的變化量與溫升成正比。根據比爾-郎伯定律,光譜輻射能在具有吸收性的介質中傳播時將按指數規律衰減。因此,當深度為z處的輻射力變化量傳播到皮膚表面時,衰減為

ΔEz0=ΔEzexp(-μaz)∝ΔT(z,t)exp(-μaz)

(7)

式中:μa為組織吸收系數。

由此,表面輻射能量沉積ΔE0,即各深度輻射力在表面的線性疊加,可表示如下

(8)

理想PPTR信號幅度ΔS(t)是由測量的表面輻射能量沉積信號ΔE0經過光電轉換處理而來,Milner等研究中推導的理想PPTR信號幅度ΔS(t)公式表達如下[8]

(9)

式中:Cd為紅外測量系統決定的比例常數。

將正問題推導所得ΔT(z,t)的關系式,代入式(9)中并對z積分,可得到理想PPTR信號幅度關于初始溫升分布的數學表達,形式如下

[erfcx(u+)-erfcx(u1)]}dz′

(10)

式中:erfcx(u)為指數互補誤差函數,表達式可寫為

erfcx(u)=exp(u2)erfc(u)

其中u1、u-、u+分別為

至此,通過推導構建了理想PPTR信號幅度ΔS(t)與初始溫升分布函數ΔT(z′, 0)之間的數學關系。將被積函數中除初始溫升之外的其他項整理成核函數K(z′,t),寫為

(11)

結合式(10)、式(11),整理可得

(12)

式(12)即為典型的第一類Fredholm積分方程。對于第一類Fredholm方程的反問題,由于該類問題解的嚴重不適定性[13],使得ΔT(z′, 0)幾乎不存在解析解,因此只能尋求數值解法。接下來,本文基于3種常用的正則化策略,即截斷奇異值分解(TSVD)法[17-18]、Tikhonov正則化[19]以及共軛梯度(CG)法[20],對上述方程進行求解。

2 PPTR數值模擬結果分析

針對1.2小節建立的PPTR皮膚結構模型反問題開展仿真模擬,基于上述常用3種算法,定義了基于自定義截斷參數選取方法的截斷奇異值分解(CTP-TSVD)法、基于L曲線法選取正則化參數[21]的Tikhonov正則化(L-Tikhonov)法以及給定迭代次數的共軛梯度(GNI-CG)正則化算法。3種算法的主函數形式相似,仿真流程可概括為以下5個步驟。

步驟1通過網格劃分,將核函數離散成核矩陣。

步驟2除了GNI-CG 算法,其余都需要對核矩陣進行奇異值分解,分別獲得左、右奇異矩陣及奇異值向量。

步驟3假設初始溫升分布函數為ΔTS(z, 0),代入至線性方程組ΔS(t)=KΔT,獲得PPTR理想溫升信號ΔS(t),然后基于給定信噪比向ΔS(t)添加白噪聲,從而獲得含擾動的測量溫升信號ΔS′(t)。由于建立的熱傳導模型中,表皮層、PWS層對脈沖激光的輻射吸收均較強,因而經激光照射后,兩層的中心位置具有最大初始溫升,然后分別向兩邊衰減。因此在激光照射結束瞬間,人體皮膚內的溫升分布函數ΔTS(z, 0),可認為符合如下超高斯分布

(13)

式中:ΔT1m、ΔT2m分別為表皮層、PWS層的最大溫升;z1、z2分別為表皮層、PWS層沿z軸方向的中心位置深度(簡稱中心深度);d1、d2分別為表皮層、PWS層的厚度;w1、w2分別為表皮層、PWS層的輪廓系數。通過調整上述參數,可模擬激光照射結束瞬間人體皮膚內的溫升分布。

步驟4調用正則化參數策略,選取最優正則化參數[18-19],執行正則化方法重建初始溫升分布ΔTC(z, 0)。

步驟5繪制重建初始溫升ΔTC(z, 0)、原初始溫升ΔTS(z, 0)曲線,通過比較兩者差異,并基于ΔTC(z, 0)、ΔTS(z, 0)兩條曲線之間的差異在ΔTS(z, 0)中的占比,建立如下誤差評估函數

(14)

參照文獻[7]選取仿真實驗參數,通過改變初始溫升分布的表皮層厚度d1、PWS層的厚度d2及中心深度z2,來分析3種算法的綜合表現。

2.1 表皮層厚度對重建結果的影響

本小節仿真參數選取:a=0.11 mm2·s-1,μa=50 mm-1,ν=0.03 mm-1,h=20 W·m-2·K-1,ΔT1m=1 K,z1=0 mm,ΔT2m=1.5 K,z2=0.23 mm,d2=0.2 mm,w1=w2=1。參考Oltulu等的研究工作[22],保持以上參數不變,選取d1分別為 0.05、0.10、0.25、0.40 mm進行仿真模擬,得到的重建曲線如圖3所示。

(a)d1=0.05 mm

由圖3可見,在不同表皮層厚度下,3種逆算法都呈現出較好的溫升重建結果。隨著表皮層厚度增大,PWS層重建深度、溫升峰值、厚度(即半峰值展寬,簡稱展寬)的變化均較小,其中心深度變化不超過0.05 mm(5%),溫升峰值變化不超過0.2 K(9.5%),展寬變化不超過0.13 mm(13%)。淺層區域的溫升分布逼近精度不高的原因,很可能是由于未能考慮表面對流、輻射等原因對反演過程的影響?？偟膩砜?3種逆算法均可滿足不同厚度表皮層結構重建的需求。

2.2 PWS層厚度對重建結果的影響

本小節仿真參數選取:a=0.11 mm2·s-1,μa=50 mm-1,ν=0.03 mm-1,h=20 W·m-2·K-1,ΔT1m=1 K,z1=0 mm,d1=0.05 mm,ΔT2m=1.5 K,z2=0.23 mm,w1=w2=1。保持以上參數不變,選取d2分別為0.2、0.4、0.6、0.8 mm,所得的重建結果如圖4所示。

(a)d2=0.2 mm

由圖4可見,隨著PWS層厚度增加,3種逆算法的重建效果整體均較好。PWS層的重建深度、溫升峰值、厚度變化仍較小,其中心深度變化不超過0.1 mm(10%),溫升峰值變化不超過0.1 K(6.7%),展寬變化不超過0.12 mm(12%)。這主要是因為PWS層厚度的增加使得溫升信號占比增大,從而能更多地反映PWS層的信息,因此重建效果較好。

2.3 PWS層中心深度對重建結果的影響

由2.1和2.2小節可知,表皮層、PWS層厚度的改變對3種逆算法的重建效果影響不大。下面分析PWS層中心深度對重建結果的影響,仿真參數分別選取:a=0.11 mm2·s-1,μa=50 mm-1,ν=0.03 mm-1,h=20 W·m-2·K-1,ΔT1m=1 K,z1=0 mm,d1=0.05 mm,ΔT2m=1.5 K,d2=0.2 mm,w1=w2=1。保持以上參數不變,選取z2分別為0.23、0.43、0.63、0.83 mm,得到的重建曲線如圖5所示。

(a)z2=0.23 mm

由圖5可見,在PWS層中心深度為0.23 mm處,重建曲線與模擬初始溫升吻合度較高,但隨著PWS層中心深度增大,重建效果變差,具體表現為:①PWS層重建深度發生偏移,當中心深度增加至0.83 mm時,偏差最大達到了0.12 mm(14.5%);②PWS層溫升峰值降低,當中心深度增加至0.83 mm時,溫升峰值最大降低了0.4 K(26.7%);③重建PWS層展寬明顯增大,當中心深度增加至0.83 mm時,重建展寬增大了0.23 mm(23%)。由此可見,PPTR溫升反演技術的使用深度有限。然而對于大部分PWS層來說,其中心深度均保持在0.4 mm以內[23]。由圖5(b)可見,當PWS層中心深度為0.43 mm 時,深度方向反演的相對誤差較小,評價指數的最大值小于0.18,PWS層的重建深度、溫升峰值、厚度變化較小,其中心深度變化不超過0.05 mm(5%),溫升峰值變化不超過0.1 K(6.7%),展寬變化不超過0.08 mm(8%),均在可接受范圍內。

2.4 綜合性能評價

表1 各逆算法在不同表皮層厚度下的評價指數Table 1 Evaluation index of each inverse algorithm under different epidermal thicknesses

表2 各逆算法在不同PWS層厚度下的評價指數Table 2 Evaluation index of each inverse algorithm under different PWS layer thicknesses

表3 各逆算法在不同PWS層中心深度下的評價指數Table 3 Evaluation index of each inverse algorithm at different PWS layer central depths

表4 各逆算法的綜合評價指數Table 4 Comprehensive evaluation index of each inverse algorithm

綜上,L-Tikhonov正則化和GNI-CG正則化算法相較于CTP-TSVD算法的綜合評價指數更小,表明前兩者的綜合表現比CTP-TSVD算法更優秀。下一步,我們將選取最優的GNI-CG算法進行實驗驗證。

3 實驗研究

上節通過數值模擬驗證了算法的有效性,結合理論分析證明了采用PPTR法進行皮膚結構探測的可行性。本節將針對臨床應用開展體外仿體實驗,進一步驗證皮膚組織結構探測的有效性。采用PPTR方法對二附院患者進行皮膚組織的溫升分布求解及結構信息探測,將探測到的信息用于激光參數選取,從而達到輔助醫生激光治療的目的。

3.1 皮膚仿體制備

皮膚仿體的制備流程如下:首先將3 g瓊脂粉溶解在100 ml蒸餾水中,然后加入2 g脂肪乳劑增加散射性,接著將混合物加熱至沸點引發聚合,最后將適量聚合瓊脂溶液注入小型培養皿,制備表皮層,將適量聚合瓊脂溶液注入大型培養皿,制備真皮層,鮮紅斑痣PWS層為厚度50 μm的粉紅色聚乙烯箔,3層的結合在水浴中開展以防止相鄰層之間形成氣泡。在吸收層下放置薄膜熱電偶,測量激光照射后吸收層的溫度。圖6(a)為制備的皮膚仿體,圖6(b)為紅外相機拍攝的圖像。

(a)皮膚仿體

3.2 皮膚仿體實驗流程

采用脈沖激光器、紅外相機、計算機、皮膚仿體搭建如圖7所示的PPTR實驗臺。使用的激光波長為1 064 nm,能量為10 J/cm2;脈沖形式為單脈沖,寬度為1 ms,光束直徑為5 mm,入射角θ為45°。同時采用幀率為120 Hz、記錄波段長度為7.5～13.5 μm的紅外相機采集仿體表面溫度信號,記錄時間為1.5 s。將采集到的表面溫度信號傳輸至計算機,采用逆算法重建獲取皮下溫度分布。

(a)示意圖

3.3 皮膚仿體實驗結果

圖8給出了PWS層中心深度z2為0.2、0.4、0.6 mm時,溫升重建結果與實驗值的對比。圖9給出了PWS層深度、展寬、最大溫升的重建結果與實驗值對比。從圖8～圖9可以看出,PWS層中心深度越淺,重建效果越好,無論是PWS深度、展寬、最大溫升,都與實驗值更為接近。當PWS層中心深度為0.2 mm 時,PWS層重建深度與實驗值高度重合,誤差低于0.1%,展寬的重建誤差為4%,最大溫升的重建誤差為3%。隨著PWS層中心深度增加,重建深度仍與實驗值保持較高重合度,誤差最大為2.5%。隨著PWS層中心深度繼續增加,在0.8 mm處,重建展寬逐漸增大,與實驗值相差0.061 mm,誤差為122%;重建的最大溫升與實驗值相差4.15℃,誤差達到了25.8%。由此可見,隨著PWS層中心深度增加,PWS層深度、最大溫升的重建效果均較好,但重建展寬即PWS層厚度的誤差卻顯著增大,這主要是因為經過組織的吸收,傳輸到表面的深層輻射信號隨著PWS層中心深度的增加而減弱,在諸如環境噪音等因素的干擾下,重建精度減弱。因而在實際應用中,需要根據圖9(b)的誤差曲線對展寬進行修正。

(a)z2=0.2 mm

3.4 臨床實驗結果

鮮紅斑痣的激光臨床治療過程中,如波長、脈寬、入射能量等激光參數的選取對治療效果起著決定性的作用。以入射能量為例,若能量過低,病變血管無法得到有效清除;若能量過高,則會引發表皮燙傷等副作用。目前,由于缺乏有效的皮膚結構檢測技術,激光參數的選取主要憑借醫師臨床經驗,而無法根據每位患者具體的皮膚情況實現個性化治療。以美國Candela激光公司波長為595 nm的Vbeam激光治療儀為例,其可選取的參數范圍包括:脈寬為0.5～40 ms,入射能量為0～20 J/cm2,頻率為1～10 Hz。如何在如此多的參數中選出最優參數組合,對于醫生來說十分困難。因而臨床醫生迫切需要患者皮膚結構探測、激光參數智能化選取的技術和手段。本文提出的脈沖輻射法在臨床應用過程中,首先可基于該方法反演出個性化患者的皮膚組織結構,其次可基于該皮膚結構信息,利用已有的參數優化軟件幫助醫生選取治療參數,輔助開展激光治療。

實驗的具體實施過程如圖10～圖11所示。首先,用安全激光劑量2 J/cm2照射患者面部皮膚表面,如圖10(a),開展PPTR測量。然后,采用紅外熱像儀探測皮膚表面溫度,如圖10(b),提取并處理激光作用區域的溫度信號,如圖10(c)。最后,利用逆算法根據探測信號,反算皮下組織的初始溫升分布和皮膚組織結構信息,如圖11。

(a)患者面部

圖11 皮下初始溫升分布重建結果Fig.11 Reconstruction of the initial subcutaneous temperature rise distribution

圖11中,紫色區域A代表具有黑色素的表皮層,紅色區域B、C、D、E均代表病變血管。由于激光的作用,表皮層及PWS層在吸收激光能量后瞬時升溫,呈現高斯溫峰形狀。由于重建的血管深度較深,為避免展寬重建誤差過大,利用圖9(b)的誤差曲線對展寬結果進行了修正。修正完成后,根據圖11,可估算出表皮層厚度dA=0.18 mm;血管B位于zB=0.55 mm處,厚度dB=0.088 mm; 血管C位于zC=0.79 mm處,厚度dC=0.032 mm; 血管D 位于zD=0.91 mm處,厚度dD=0.021 mm; 血管E位于zE=0.98 mm處,厚度dE=0.01 mm。上述重建結果均符合患者面部特征下的皮膚結構參數范圍。此外,本模型對于皮下血管溫度的重建效果較好,有望實現激光熱療過程中的溫度在線監測。

基于PPTR方法進行皮膚結構探測,反演出患者的皮膚結構信息后,可采用激光參數優化軟件[24](圖12)進行治療參數優化,軟件的計算過程如下。

圖12 參數優化軟件界面Fig.12 Parameter optimization software

步驟1以經典的多層均勻模型為基礎,構建出患者的皮膚結構模型。

步驟2基于光傳輸模型[25],計算激光能量在患者皮膚內的傳播及分布情況。

步驟3基于生物傳熱模型[25],計算激光作用后皮膚組織的溫度分布。

步驟4基于熱損傷模型[26],以PWS層血管達到損傷閾值為條件,進一步優化如脈寬、能量等治療參數。

針對圖10(a)中的患者,利用脈沖光熱輻射模型反演的個性化皮膚結構如圖11所示,在激光參數優化軟件內輸入結構參數信息后,即可得到推薦的激光治療參數,如圖12所示。最終醫生可以此為依據,選取合適的激光參數對患者開展治療。

4 結 論

本文針對血管性皮膚病激光熱療中,激光參數選取無依據、臨床療效不佳的瓶頸問題,提出了一種基于脈沖光熱輻射法重建皮膚結構,智能化選取激光參數用以輔助激光治療的模型。采用脈沖光熱輻射法建立數值模型,利用MATLAB軟件開展仿真實驗,通過調整模型輪廓等參數,全面分析了CTP-TSVD、L-Tikhonov以及GNI-CG 3種正則化逆算法對皮下血管深度、厚度及溫度的重建效果,發現GNI-CG正則化算法的綜合性能最好,L-Tikhonov正則化算法綜合表現次之,CTP-TSVD算法的綜合表現最差。然而無論何種方法,對血管層深度及溫度的反演精度均較好,而對血管層厚度的反演精度稍差,需要采用實驗值進行修正。接著,基于模擬結果進一步開展體外實驗,驗證了模型的有效性。最后,將最優算法GNI-CG應用于臨床研究,為臨床醫生的激光光熱治療血管性皮膚病提供了定量化的數據支持。