在指端選取血管密集區域并在其上方施加適當壓力可改善指端容積脈搏信號幅值。針對這一現象，本文運用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）建立指尖單支小動脈血管組織模型進行仿真。基于有限元法求解三維納維-斯托克斯（Navier-Stokes）方程，計算血液的速度場和壓力分布，對血管及周邊組織進行形變分析。基于朗伯比爾定律（Lambert Beer’s law）探究橫向光出射面區域及組織模型縱向壓縮位移對出射光強度信號的影響。研究結果表明，選取血管較密集區域可以增大單支血管組織模型的光強度信號幅值，減小其峰值；向組織施加一定深度的下壓會增大信號的幅值及峰值。最終，期望通過本文仿真研究結果結合以往實驗經驗，可為提高指端容積脈搏信號質量提供切實可行的方案。
引用本文： 趙鶴童, 周亮, 劉朝暉, 喬文龍, 孫笑敩, 姜樂, 呂媛媛. 有限壓力下皮膚組織模型的血管動力學及光學耦合分析. 生物醫學工程學雜志, 2022, 39(3): 527-536. doi: 10.7507/1001-5515.202106039 復制
引言
在生物醫學光子學領域，光電容積描記術（photoplethysmography，PPG）具備無創、便捷、低成本等優勢，是目前研究熱點之一[1]。PPG利用光電傳感器，檢測經過人體血液和組織吸收后的不同反射光強度，描記出血管容積在心動周期內的變化，能夠從得到的脈搏波形中計算出心率等生理信息。然而，以光學方法提取的脈搏波容易受到多種因素影響，其檢測精度仍有待提高。包括血管壁和組織彈性、血液動力黏度、測量區域血管密集程度差異以及外界環境施加給組織表層的壓力等因素，都會引起光電容積脈搏波波形的變化[2-3]。
PPG由Hertzman等[4]率先提出，是基于皮下組織血管容積變化導出光強度信號變化而得到的光電容積脈搏波的技術。Womersley[5]通過納維-斯托克斯方程（Navier-Stokes equations，NSE）計算出血流速度及血液對血管壁周邊壓力的分布。劉永樂[6]結合三維流固耦合方程及橈動脈模型內血液的流動，分析了泵血橈動脈初始邊界條件的表達方程。李頂立[7]提出以朗伯-比爾定律（Lambert-Beer law，LBL）從數學的角度來理解單支血管血液容積變化對光電容積脈搏波透射光強信號的影響。劉靜紈[8]在橈動脈及指脈之間建立了指端微循環容積脈搏波血流模型，首次把心血管雙彈性腔模型理論擴大并推廣到微循環領域，從理論上表示了PPG在微循環系統中檢測的容積脈搏波。王靜宜[9]基于蒙特卡羅方法（Monte-Carlo method，MCM）進行了光在人體組織內的傳輸分析。
對脈搏波的建模一般基于力學或電學進行研究分析[10]。此外，Sherwin等[11]首次提出了基于數值模擬時空變量的血管網絡模型。而后陸續提出多種優化方法以完善該模型，如低夏皮羅水壓重建法（the low-Shapiro hydro-static reconstruction，HR-LS）[12]、龍格-庫塔間斷伽遼金方法（Runge–Kutta discontinuous Galerkin methods，RKDG）[13]、長波納維-斯托克斯方程耦合法（long-wave Navier–Stokes equations coupled methods，LWNS）[14]。本文基于有限元數值模擬法、利用納維-斯托克斯方程及朗伯-比爾定律，分析在不同血管密度的區域和適當組織下壓深度的血管組織模型上，透射光能量強度信號的變化情況。
依據甲板生理結構[17]，本研究運用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）完成對指端甲板深叢層單支小動脈血管組織模型的建模工作，在軟件中配置力學、光學仿真相關的生理參數、計算方程及邊界條件。基于不連續試函數和基函數的有限元數值計算分析法，求解模型的三維納維-斯托克斯方程，得出數值解來解決計算流體動力學[15]（computational fluid dynamics，CFD）的流固耦合問題[16]，實現模擬“泵血”的動態過程。然后，本文進一步計算出血液的速度場和壓力分布，再對血管及周邊組織進行形變分析，得出模型幾何變化量。通過動力學仿真數據進行模型重建，而后基于朗伯-比爾定律模擬均勻單色光源與組織之間的光能量強度損耗的光學作用，實現對光學特性的仿真研究分析。本文從有限元數值仿真角度，驗證了采取適當壓力及選取血管較密集區，會改善脈搏波信號質量，或可為提高指端容積脈搏信號質量提供切實可行的方案。
1 基于單支血管的簡化組織模型的建立
在人體中，作為人體血液輸送載體的血管主要分三大類：動脈、靜脈和毛細血管網。血液由作為泵血源頭的心臟泵出后，依次經過動脈、毛細血管網、靜脈。動、靜脈均有大、中、小、微之分。在微循環中，血流通常是從小、微動脈流入毛細血管網，再流至微、小靜脈。血管的管徑會依次變化，其壁厚也會因血管膜的中膜內的平滑肌層數變化而改變。不同血管類型的管徑范圍及主要成分如表1所示。
表1
各類血管尺寸及管膜主要成分
Table1.
Sizes and main components of tubules of different types of blood vessels
指端甲板深叢層單支小動脈血管組織模型是參考郭家松等[17]對甲床動脈分布的研究，取指掌側固有橈動脈的血管尺寸參數建立的。該模型的建立旨在分析單支血管與周邊組織受周期泵血及外界下壓的影響而產生的形變量、出射光強度信號區域（血管密集程度）Lx及組織下壓深度L對出射光強度信號的影響。為便于后續仿真分析，組織設定為規則長方體，組織尺寸為3 mm × 2 mm × 2 mm，取管徑為460 μm，血液直徑為400 μm，壁厚30 μm。建立的指端單支小動脈血管組織模型如圖1所示。依據力學仿真計算每時刻模型隨“泵血”及組織下壓引起的血管擴張變化數據，并在重建二維光學模型時，分為血液、血管及組織兩部分區域進行分析。
圖1
指端單支小動脈血管組織模型
Figure1.
The model of single branch vascular tissue
2 基于有限元分析法對組織模型的力學特性進行分析
2.1 力學仿真分析流程
依據甲板生理結構、運用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）完成對指端甲板深叢層小動脈的模型的建模工作后，定義模型的組織、血管壁、血液的力學參數以及血液動力學壓力條件和相關邊界條件，即定義血液與血管假設無相對滑移、血管壁與組織位移連續、血流進出口壓力參數、模型外表面為輥支撐邊界，使得組織邊界不會隨血流產生垂直方向上的位移[6]。為考慮組織及血管向兩側擠壓時的力學結構變化，設定給予模型上表面適當下壓深度條件，分別為：0、0.1、0.2、0.3、0.4、0.5 mm；而后將模型劃分有限網格，利用不連續試函數和基函數的有限元法求解三維納維-斯托克斯方程，得出血液的速度場和壓力分布，再對血管及周邊組織進行形變分析[6]。
2.2 模型力學參數的簡化及定義
在對生物組織進行建模中，力學特征通常表征為非線性。其模型有多種，如粘彈性模型、超彈性模型、線彈性模型等。在指端處血管隨泵血產生的形變量微小，將力學應變視為小應變極限，模型中的組織、血管壁定義為線彈性模型。模型血管選取為指端血管深叢層中的小動脈區域。參考趙瑞剛等[18]給定的人體血管力學參數，小動脈中膜由1～4層平滑肌細胞組成，大動脈中膜中含有40～70層的彈性纖維形成的彈性膜，適當調整得出本仿真中所涉及到的組織、血管壁、血液的力學參數（密度ρ，拉梅參數λ、μ，動力黏度η），如表2所示。
表2
單支血管組織模型力學參數
Table2.
Mechanical parameters of single-branch vascular tissue model
依據彈性腔的相關理論及成年人心臟泵血生理參數來定義血液動力學輸入、輸出壓力條件。輸入壓力條件是由彈性腔理論推出的主動脈血液壓力條件進行放縮得出小動脈中的血流壓力函數。仿真主要探究選取出射區域（血液密集度）、下壓深度對出射光強度信號幅值及峰值（主峰）的影響。在設置壓力條件時，為降低計算量，不考慮血壓及脈搏波中引起的潮波及重搏波現象。
定義血液動力輸入、輸出壓力條件時，應先由彈性腔理論導出的主動脈血液入口壓力條件[19-20]。在彈性腔理論假定的腔體容積與血壓之間呈線性關系，靜脈血壓近似為零，心室射血瞬時完成（即 
）的情況下，得出主動脈血液入口壓力條件p主動脈(t)函數，如式(1)所示：
 |
其中，Q為血管內血流充盈程度，R為體表外層組織容抗性，C為脈管的順暢性，Ts、Tpulse分別為心室收縮期和單次脈搏時間， Pd、Ps分別為舒張壓（diastolic blood pressure，DBP）和收縮壓（systolic blood pressure，SBP）。
血管內血流充盈程度與體表外層組織容抗性，兩種參數的系數關系，如式(2)、式(3)所示。通過式(2)、式(3)，可以將主動脈血液入口壓力條件用心室收縮期Ts、單次脈搏時間Tpulse、舒張壓Pd、收縮壓Ps參數進行表征。
 |
 |
以心率Hr = 80 次/min，舒張壓Pd = 75 mm Hg，收縮壓Ps = 115 mm Hg的成年人心臟泵血生理參數作為壓力條件的相關參數。泵血周期Tpulse = 60/Hr = 0.75 s，心室收縮期Ts = 0.15 s。血管開始擴張t = 0 s，血容積處于最大狀態t = Ts。將以上參數代入p主動脈(t)后，將分段函數連接處以連續一階導數的方式進行平滑處理。
在定義血液動力學輸入條件時，考慮到模型血管取自指端處深叢層中的小動脈區域血管。其管徑遠小于主動脈，處于遠離心臟的位置。將主動脈血液入口壓力條件縮小，將其時間取至3 s，循環4個周期，如圖2所示。模型血液輸入壓力條件定義如式(4)所示：
圖2
單支血管組織模型血液流動入口壓力條件
Figure2.
Inlet pressure conditions of blood flow in single-branch vascular tissue model
 |
在定義血液動力學輸出條件時，參考趙瑞剛等[18]給定的人體血管力學參數，可得出血管血壓損耗0.008 mm Hg/mm，該簡化模型出口為單出口，模型長度為3 mm，即可得模型出入口損耗壓強lossmodle = 0.024 mm Hg。模型血液出口壓力條件可以定義為集中參數函數，如式(5)所示：
 |
其中p0max是模型血流入口峰值。注：在定義血液輸出壓力條件時，也可以將其恒定置為0 Pa。其相應的血管位移會存在相應的失真，但仿真求解出的幾何特征變化趨勢大致不變。即，在滿足仿真幾何形變正常的情況下，依然可以得出與下文所述相同的結論。
由于血液分子的尺度是納米量級（10?9 m），血管管徑尺寸近似毫米量級（10?4 m），相差約5個量級，而血液內部血紅細胞尺寸為5 μm左右，也遠小于血管管徑。在空間角度，血液可以理解為微觀上充分大、宏觀上充分小的分子團集合。在時間角度上分析，泵血周期定義為0.75 s，泵血的時間尺度遠大于其對應的弛豫時間（即流體質點在失去平衡后通過分子間的碰撞最終達到新的平衡狀態所經歷的不平衡時間），可以理解為泵血過程微觀充分長、宏觀充分短，故在建立好的三維單支血管組織模型上滿足流體連續性假設條件。此外，將血液視為粘性流體，流動時存在內摩擦力，用動力黏度表征其大小。仿真時血流按照層流運動來分析，即各層流體僅做相對滑動，互不混合。在血液層流域中求解納維-斯托克斯方程，進而求解出模型中血液的速度場和壓力分布及血管壁內側受到的壓力載荷分布。
2.3 力學分析與結果
本研究選取6種組織下壓深度，即L = 0、0.1、0.2、0.3、0.4、0.5 mm進行力學仿真。為了避免首個周期數據出現初始特殊波動，在獲取穩定泵血循環過程中的數據時，故提取第三泵血周期數據（t = 1.5～2.25 s）進行分析。首先，將組織下壓深度L = 0 mm的數據進行處理，得出內外壁位移大小隨時間波動關系及血液中心點血液流速隨時間波動關系，如圖3所示。從圖3兩幅子圖中可看出，第二個周期后，位移及流速波動相對穩定。血液最高流速vmax = 0.193 5 cm/s，血最大位移量為0.135 mm，則血液最大半徑rmax = 0.335 mm與實際流速相似。血液只做層流運動時，要滿足雷諾數Re = ρblood·vmax·rmax/η<2 000[21]。將代入雷諾數公式求得雷諾數Re = 1.37。依照Lew[22]對血液流雷諾數的分類可知，此時為小雷諾數均勻發展流動。綜上，符合仿真起初血流為層流的假設。
圖3
無組織下壓時，泵血模型的位移量及血液流速
Figure3.
The displacement and blood velocity of the pump model without tissue pressure
在6種不同組織下壓深度L的作用下，橫、縱向管壁厚度在泵血周期內的變化情況及最小值如圖4所示。由橫、縱向血管厚度變化的兩幅圖可以看出，隨著壓力深度的增大，其橫向血管壁厚度逐漸變大，縱向血管壁逐漸變小。由厚度最小值變化的兩幅圖可以看出，血管壁厚最小值與組織下壓深度L呈線性關系，且當下壓深度L = 0.5 mm時，縱向管壁厚度出現負值，預示著管壁破裂。此外，血管壁厚度會隨著心臟泵血、血流壓力而波動，軸向血管半徑的變化趨勢與其壁厚變化趨勢相同。
圖4
橫向、縱向血管壁厚度隨泵血產生的變化
Figure4.
Changes of lateral and longitudinal vessel wall thickness with pumping
對比血液體積、血管壁體積和血液與血管壁總體積在不同組織下壓深度L下隨時間變化情況以及在t = 1.66 s時，三者隨組織下壓深度L變化的關系，如圖5所示。可以得出這三種體積在同一時刻的變化趨勢相同。由t = 1.66 s時隨組織下壓深度變化的關系可得，三者相對比，血管壁體積受組織下壓的影響小。橫、縱向血管尺寸（內、外半徑及壁厚）隨著組織下壓深度L變化程度大致相同。
圖5
血液、血管壁隨泵血的體積變化趨勢
Figure5.
The variation trend of blood and vascular wall with the volume of blood pumped
3 基于有限元分析法對組織模型的光學分析
3.1 光學仿真分析流程
由于各介質的吸收系數不同，血管的泵血搏動會引起透射光強度的波動。且血液壓力波動 
與透射光強度波動 
存在線性關系，即 
[7]。本研究在不考慮光散射損失的情況下，依據朗伯-比爾定律模擬均勻單色光源與組織之間的光學吸收作用。光學模型按區域分為血液、血管與組織兩部分來定義參數。其幾何參數是由動力學仿真計算出的每時刻模型血管擴張變化數據定義的，其光學參數根據理論皮膚組織光譜吸收系數得出。在軟件中設定朗伯-比爾定律計算方程以及模型初始邊界條件，包括定義入射強度、將模型外邊定義為透明表面等。該邊界條件使得邊界處不積蓄能量，且進入吸收介質的輻射光束的強度為零。再將模型劃分有限網格進行有限元計算。得出利用朗伯-比爾定律模擬均勻單色光源與組織之間光學作用的結果。
3.2 模型光學參數的簡化及定義
皮膚自表層往里可分為3層——表皮層（epidermis，EP）、真皮層（dermis，DE）、皮下組織層（subcutishypodermis，SH）。真皮層依照含血量分為4類——微循環層（papilmry dermis，PA）、上叢層（superficial arteriovenous plexus，SP）、網狀真皮層（reticular dermis，RD）、深叢層（deep arteriovenous plexus，DP）[23]。模型的血管按尺寸應處于深叢層，研究出射面光強度分布受組織下壓深度L的影響，參考周亮等[23]提出的理論皮膚組織光譜吸收系數，定義血液吸收系數μblood=5 cm?1、組織吸收系數μtunicae = 3 cm?1。設定入射沉積光束功率密度I0 = 5 W/m2，入射方向為垂直向下。入射沉積光束被組織吸收后轉化成熱能很小，故忽略熱能轉換量。模型中的熱能參數，如導熱系數、恒壓熱容、比熱率等均為仿真參考值，對仿真分析出射面光強信號無影響，可不考慮。
3.3 光學分析與結果
先將組織下壓深度L = 0 mm的數據進行處理，得出無組織下壓、峰值時刻（t=1.66 s）時，模型二維截面的光強度分布信息。再對比6種組織下壓深度L = 0、0.1、0.2、0.3、0.4、0.5 mm時，橫坐標直線上的強度分布，如圖6所示。可以看出隨著組織下壓深度L的增大，出射點的光強度逐漸變大，并保留原有的線性趨勢。
圖6
不同組織下壓深度時，模型光能量強度分布圖
Figure6.
Distribution of model light energy intensity under different pressing depths of tissues
由于在實際接收、檢測光強度信號時，只能選取一塊感興趣區域（region of interest，ROI），無法聚焦于某一點。在仿真中，要研究出射面一定區域的出射光強度平均值。以組織下壓深度L = 0 mm模型為例，當出射面以中心區域Lx為 0.2 mm開始，逐漸增大至2.0 mm時，光能量強度平均值與時間t的關系，如圖7所示。
圖7
光能量強度平均值與時間的關系
Figure7.
The relationship between average light intensity and time
提取光強度信號的幅值及峰值，對比選取的橫向出射面區域Lx對信號幅值、峰值的影響，如圖8所示。在單支血管模型中，要選取合適的橫向光出射面區域，獲取的信號幅值才能為最大，選取的區域過大或過小，都會降低信號幅值；選取的橫向光出射面區域越大，其峰值會越大。
圖8
無組織下壓時的信號幅值、峰值
Figure8.
Signal amplitude and peak value without down tissue pressure
模型選取不同組織下壓深度L時，當處于心臟膨脹期（t = 1.48 s）時，血管直徑最小，處于心臟收縮期（t = 1.66 s）時，血管直徑最大。通過處理上述模型的動力學仿真數據，可得出血液直徑波動數據如表3所示。
表3
不同下壓深度時，模型血液直徑數據
Table3.
Data of model blood diameter at different pressing depths
由表3可知，組織下壓深度L = 0 mm模型的血液直徑在0.557～0.703 mm之間波動。橫向光出射面區域Lx > 0.703 mm的幅值點處在圖8信號幅值圖曲線的下降沿處，即當橫向光出射面選取區域大于血液直徑波動最大值時，信號的幅值會隨著橫向光出射面Lx選取區域增大而減小，信號峰值會隨之增大而增大。實際應用中，選取的ROI區域有多支血管分布，選取的單支區域一定大于血管直徑波動最大值。可以推測，實際上在血管組織同等厚度情況下，血管相對密集的區域意味著選取的單支血管橫向光出射面區域Lx更接近血管直徑波動最大值。在此處進行提取光強度信號，幅值相對血管稀疏區（單支血管橫向光出射面大于血管直徑波動最大值）要大一些。
探究當橫向光出射面區域大于血管尺寸（Lx > 0.917 mm）時，信號幅值及峰值隨組織下壓深度L及橫向光出射面區域Lx的變化情況，如圖9所示。結合4幅子圖綜合分析可得，信號幅值及峰值均隨組織下壓深度L的增大而增大。且橫向光出射面區域越小，其幅值增大的效果更明顯。在模型仿真中，選取橫向光出射面區域Lx接近血液直徑波動最大值，且給予適當的下壓深度L，其出射面區域的光能量強度信號幅值和峰值均會增大。在實際透射組織厚度相同情況下，選取血液密集處，給予適當的壓力一定程度上可以提高信號質量。
圖9
不同組織下壓深度、橫向光出射面區域時的信號幅值、峰值
Figure9.
Signal amplitude and peak value of different tissue pressing depths and transverse light emitting surface area
當給予單支血管組織模型下壓深度時，血液區域及組織區域都會產生形變，從而影響出射面區域能量強度信號峰值、幅值。將血液區域的擠壓變形、組織區域的擠壓變薄兩種物理效果進行分離建模，探究血液形變、組織厚度形變二者對信號影響的程度。為貼近實際，僅考慮橫向光出射面區域大于血液直徑波動最大值的情況，即選取橫向光出射面區域Lx > 0.917 mm的數據進行處理分析。將原有模型數據與組織和血液分別受組織下壓深度L影響的模型數據結合在一起，對比不同模型下，組織下壓深度L與信號幅值、峰值之間的關系，如圖10所示，其中a模型為組織和血液均受組織下壓的影響，b模型為僅組織受影響，c模型為僅血液受影響。可以看出，血液區域擠壓形變對信號幅值及峰值的影響均不大，增大信號幅值的關鍵因素是組織區域的形變。實際中，血液和組織均受到下壓效果引起形變，而此時對信號的影響最大。
圖10
不同模型及深度時的信號幅值、峰值
Figure10.
Signal amplitude and peak value at different models and depths
4 結論與展望
本文利用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）根據組織血管的基本結構及泵血物理過程，參照郭家松等[17]對指尖甲床動脈的尺寸的研究，定義指端單支小動脈血管組織模型的血管尺寸；依據劉永樂[6]對血管組織動力學上的研究等，完成了模型在泵血周期內產生的幾何形變特征的動力學計算；依據李頂立[7]對單只血容積光電容積脈搏波的數學分析等，完成對模型的光學分析。本研究探究了出射面光強度信號的幅值、峰值受選取橫向光出射面區域和組織下壓深度的影響。從仿真角度得出：在模型組織上施加適當的下壓，會造成橫向血管壁變厚，縱向血管壁變薄；選取血管密集區域相對于血管稀疏區的信號幅值更高，峰值較低；對組織給予適當的下壓深度，會增大出射面光強信號的幅值及峰值。驗證了王焱等[2]提出的適當壓力會改善在動脈血液檢測中的脈搏波信號質量及余江軍等[3]提出的適當壓力對指尖自適應脈搏信號具有改善作用的理論。
本文的力學仿真屬于流體與固體耦合計算，為減少計算量，設定為單向作用耦合，血管壁及組織的形變后的幾何模型不會反過來影響血液流動，即無法分析按壓組織對血液流速的影響。在光學方面，依據朗伯-比爾定律考慮光損耗效果，未考慮光色散情況。以上簡化，一定程度上提高了仿真中的組織下壓操作對信號幅值及峰值的影響效果。在后續的研究中，可以設定雙向耦合物理場，建立多支血管組成血管網模型，來探究血管分布等因素對出射光信號及血液產生的影響。如調整相關設定：設定多支血管不同排列方式、多支血管泵血存在延時、多支血管分布過密集（受按壓后，血管相互抵觸阻礙形變）等。在雙向耦合場中，可以實現對心臟動脈瓣關閉，血液反流等因素形成的潮波、脈搏波現象的仿真；建立指端血管網模型，探究其出射光強度分布，對其ROI區域的提取進行自適應分析。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：由趙鶴童主要負責仿真的全過程及后續數據分析，由周亮、劉朝暉、喬文龍等負責實驗指導，數據分析指導，由喬文龍等負責論文審閱修訂。
引言
在生物醫學光子學領域，光電容積描記術（photoplethysmography，PPG）具備無創、便捷、低成本等優勢，是目前研究熱點之一[1]。PPG利用光電傳感器，檢測經過人體血液和組織吸收后的不同反射光強度，描記出血管容積在心動周期內的變化，能夠從得到的脈搏波形中計算出心率等生理信息。然而，以光學方法提取的脈搏波容易受到多種因素影響，其檢測精度仍有待提高。包括血管壁和組織彈性、血液動力黏度、測量區域血管密集程度差異以及外界環境施加給組織表層的壓力等因素，都會引起光電容積脈搏波波形的變化[2-3]。
PPG由Hertzman等[4]率先提出，是基于皮下組織血管容積變化導出光強度信號變化而得到的光電容積脈搏波的技術。Womersley[5]通過納維-斯托克斯方程（Navier-Stokes equations，NSE）計算出血流速度及血液對血管壁周邊壓力的分布。劉永樂[6]結合三維流固耦合方程及橈動脈模型內血液的流動，分析了泵血橈動脈初始邊界條件的表達方程。李頂立[7]提出以朗伯-比爾定律（Lambert-Beer law，LBL）從數學的角度來理解單支血管血液容積變化對光電容積脈搏波透射光強信號的影響。劉靜紈[8]在橈動脈及指脈之間建立了指端微循環容積脈搏波血流模型，首次把心血管雙彈性腔模型理論擴大并推廣到微循環領域，從理論上表示了PPG在微循環系統中檢測的容積脈搏波。王靜宜[9]基于蒙特卡羅方法（Monte-Carlo method，MCM）進行了光在人體組織內的傳輸分析。
對脈搏波的建模一般基于力學或電學進行研究分析[10]。此外，Sherwin等[11]首次提出了基于數值模擬時空變量的血管網絡模型。而后陸續提出多種優化方法以完善該模型，如低夏皮羅水壓重建法（the low-Shapiro hydro-static reconstruction，HR-LS）[12]、龍格-庫塔間斷伽遼金方法（Runge–Kutta discontinuous Galerkin methods，RKDG）[13]、長波納維-斯托克斯方程耦合法（long-wave Navier–Stokes equations coupled methods，LWNS）[14]。本文基于有限元數值模擬法、利用納維-斯托克斯方程及朗伯-比爾定律，分析在不同血管密度的區域和適當組織下壓深度的血管組織模型上，透射光能量強度信號的變化情況。
依據甲板生理結構[17]，本研究運用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）完成對指端甲板深叢層單支小動脈血管組織模型的建模工作，在軟件中配置力學、光學仿真相關的生理參數、計算方程及邊界條件。基于不連續試函數和基函數的有限元數值計算分析法，求解模型的三維納維-斯托克斯方程，得出數值解來解決計算流體動力學[15]（computational fluid dynamics，CFD）的流固耦合問題[16]，實現模擬“泵血”的動態過程。然后，本文進一步計算出血液的速度場和壓力分布，再對血管及周邊組織進行形變分析，得出模型幾何變化量。通過動力學仿真數據進行模型重建，而后基于朗伯-比爾定律模擬均勻單色光源與組織之間的光能量強度損耗的光學作用，實現對光學特性的仿真研究分析。本文從有限元數值仿真角度，驗證了采取適當壓力及選取血管較密集區，會改善脈搏波信號質量，或可為提高指端容積脈搏信號質量提供切實可行的方案。
1 基于單支血管的簡化組織模型的建立
在人體中，作為人體血液輸送載體的血管主要分三大類：動脈、靜脈和毛細血管網。血液由作為泵血源頭的心臟泵出后，依次經過動脈、毛細血管網、靜脈。動、靜脈均有大、中、小、微之分。在微循環中，血流通常是從小、微動脈流入毛細血管網，再流至微、小靜脈。血管的管徑會依次變化，其壁厚也會因血管膜的中膜內的平滑肌層數變化而改變。不同血管類型的管徑范圍及主要成分如表1所示。
表1
各類血管尺寸及管膜主要成分
Table1.
Sizes and main components of tubules of different types of blood vessels
指端甲板深叢層單支小動脈血管組織模型是參考郭家松等[17]對甲床動脈分布的研究，取指掌側固有橈動脈的血管尺寸參數建立的。該模型的建立旨在分析單支血管與周邊組織受周期泵血及外界下壓的影響而產生的形變量、出射光強度信號區域（血管密集程度）Lx及組織下壓深度L對出射光強度信號的影響。為便于后續仿真分析，組織設定為規則長方體，組織尺寸為3 mm × 2 mm × 2 mm，取管徑為460 μm，血液直徑為400 μm，壁厚30 μm。建立的指端單支小動脈血管組織模型如圖1所示。依據力學仿真計算每時刻模型隨“泵血”及組織下壓引起的血管擴張變化數據，并在重建二維光學模型時，分為血液、血管及組織兩部分區域進行分析。
圖1
指端單支小動脈血管組織模型
Figure1.
The model of single branch vascular tissue
2 基于有限元分析法對組織模型的力學特性進行分析
2.1 力學仿真分析流程
依據甲板生理結構、運用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）完成對指端甲板深叢層小動脈的模型的建模工作后，定義模型的組織、血管壁、血液的力學參數以及血液動力學壓力條件和相關邊界條件，即定義血液與血管假設無相對滑移、血管壁與組織位移連續、血流進出口壓力參數、模型外表面為輥支撐邊界，使得組織邊界不會隨血流產生垂直方向上的位移[6]。為考慮組織及血管向兩側擠壓時的力學結構變化，設定給予模型上表面適當下壓深度條件，分別為：0、0.1、0.2、0.3、0.4、0.5 mm；而后將模型劃分有限網格，利用不連續試函數和基函數的有限元法求解三維納維-斯托克斯方程，得出血液的速度場和壓力分布，再對血管及周邊組織進行形變分析[6]。
2.2 模型力學參數的簡化及定義
在對生物組織進行建模中，力學特征通常表征為非線性。其模型有多種，如粘彈性模型、超彈性模型、線彈性模型等。在指端處血管隨泵血產生的形變量微小，將力學應變視為小應變極限，模型中的組織、血管壁定義為線彈性模型。模型血管選取為指端血管深叢層中的小動脈區域。參考趙瑞剛等[18]給定的人體血管力學參數，小動脈中膜由1～4層平滑肌細胞組成，大動脈中膜中含有40～70層的彈性纖維形成的彈性膜，適當調整得出本仿真中所涉及到的組織、血管壁、血液的力學參數（密度ρ，拉梅參數λ、μ，動力黏度η），如表2所示。
表2
單支血管組織模型力學參數
Table2.
Mechanical parameters of single-branch vascular tissue model
依據彈性腔的相關理論及成年人心臟泵血生理參數來定義血液動力學輸入、輸出壓力條件。輸入壓力條件是由彈性腔理論推出的主動脈血液壓力條件進行放縮得出小動脈中的血流壓力函數。仿真主要探究選取出射區域（血液密集度）、下壓深度對出射光強度信號幅值及峰值（主峰）的影響。在設置壓力條件時，為降低計算量，不考慮血壓及脈搏波中引起的潮波及重搏波現象。
定義血液動力輸入、輸出壓力條件時，應先由彈性腔理論導出的主動脈血液入口壓力條件[19-20]。在彈性腔理論假定的腔體容積與血壓之間呈線性關系，靜脈血壓近似為零，心室射血瞬時完成（即 ）的情況下，得出主動脈血液入口壓力條件p主動脈(t)函數，如式(1)所示：
|
其中，Q為血管內血流充盈程度，R為體表外層組織容抗性，C為脈管的順暢性，Ts、Tpulse分別為心室收縮期和單次脈搏時間， Pd、Ps分別為舒張壓（diastolic blood pressure，DBP）和收縮壓（systolic blood pressure，SBP）。
血管內血流充盈程度與體表外層組織容抗性，兩種參數的系數關系，如式(2)、式(3)所示。通過式(2)、式(3)，可以將主動脈血液入口壓力條件用心室收縮期Ts、單次脈搏時間Tpulse、舒張壓Pd、收縮壓Ps參數進行表征。
|
|
以心率Hr = 80 次/min，舒張壓Pd = 75 mm Hg，收縮壓Ps = 115 mm Hg的成年人心臟泵血生理參數作為壓力條件的相關參數。泵血周期Tpulse = 60/Hr = 0.75 s，心室收縮期Ts = 0.15 s。血管開始擴張t = 0 s，血容積處于最大狀態t = Ts。將以上參數代入p主動脈(t)后，將分段函數連接處以連續一階導數的方式進行平滑處理。
在定義血液動力學輸入條件時，考慮到模型血管取自指端處深叢層中的小動脈區域血管。其管徑遠小于主動脈，處于遠離心臟的位置。將主動脈血液入口壓力條件縮小，將其時間取至3 s，循環4個周期，如圖2所示。模型血液輸入壓力條件定義如式(4)所示：
圖2
單支血管組織模型血液流動入口壓力條件
Figure2.
Inlet pressure conditions of blood flow in single-branch vascular tissue model
|
在定義血液動力學輸出條件時，參考趙瑞剛等[18]給定的人體血管力學參數，可得出血管血壓損耗0.008 mm Hg/mm，該簡化模型出口為單出口，模型長度為3 mm，即可得模型出入口損耗壓強lossmodle = 0.024 mm Hg。模型血液出口壓力條件可以定義為集中參數函數，如式(5)所示：
|
其中p0max是模型血流入口峰值。注：在定義血液輸出壓力條件時，也可以將其恒定置為0 Pa。其相應的血管位移會存在相應的失真，但仿真求解出的幾何特征變化趨勢大致不變。即，在滿足仿真幾何形變正常的情況下，依然可以得出與下文所述相同的結論。
由于血液分子的尺度是納米量級（10?9 m），血管管徑尺寸近似毫米量級（10?4 m），相差約5個量級，而血液內部血紅細胞尺寸為5 μm左右，也遠小于血管管徑。在空間角度，血液可以理解為微觀上充分大、宏觀上充分小的分子團集合。在時間角度上分析，泵血周期定義為0.75 s，泵血的時間尺度遠大于其對應的弛豫時間（即流體質點在失去平衡后通過分子間的碰撞最終達到新的平衡狀態所經歷的不平衡時間），可以理解為泵血過程微觀充分長、宏觀充分短，故在建立好的三維單支血管組織模型上滿足流體連續性假設條件。此外，將血液視為粘性流體，流動時存在內摩擦力，用動力黏度表征其大小。仿真時血流按照層流運動來分析，即各層流體僅做相對滑動，互不混合。在血液層流域中求解納維-斯托克斯方程，進而求解出模型中血液的速度場和壓力分布及血管壁內側受到的壓力載荷分布。
2.3 力學分析與結果
本研究選取6種組織下壓深度，即L = 0、0.1、0.2、0.3、0.4、0.5 mm進行力學仿真。為了避免首個周期數據出現初始特殊波動，在獲取穩定泵血循環過程中的數據時，故提取第三泵血周期數據（t = 1.5～2.25 s）進行分析。首先，將組織下壓深度L = 0 mm的數據進行處理，得出內外壁位移大小隨時間波動關系及血液中心點血液流速隨時間波動關系，如圖3所示。從圖3兩幅子圖中可看出，第二個周期后，位移及流速波動相對穩定。血液最高流速vmax = 0.193 5 cm/s，血最大位移量為0.135 mm，則血液最大半徑rmax = 0.335 mm與實際流速相似。血液只做層流運動時，要滿足雷諾數Re = ρblood·vmax·rmax/η<2 000[21]。將代入雷諾數公式求得雷諾數Re = 1.37。依照Lew[22]對血液流雷諾數的分類可知，此時為小雷諾數均勻發展流動。綜上，符合仿真起初血流為層流的假設。
圖3
無組織下壓時，泵血模型的位移量及血液流速
Figure3.
The displacement and blood velocity of the pump model without tissue pressure
在6種不同組織下壓深度L的作用下，橫、縱向管壁厚度在泵血周期內的變化情況及最小值如圖4所示。由橫、縱向血管厚度變化的兩幅圖可以看出，隨著壓力深度的增大，其橫向血管壁厚度逐漸變大，縱向血管壁逐漸變小。由厚度最小值變化的兩幅圖可以看出，血管壁厚最小值與組織下壓深度L呈線性關系，且當下壓深度L = 0.5 mm時，縱向管壁厚度出現負值，預示著管壁破裂。此外，血管壁厚度會隨著心臟泵血、血流壓力而波動，軸向血管半徑的變化趨勢與其壁厚變化趨勢相同。
圖4
橫向、縱向血管壁厚度隨泵血產生的變化
Figure4.
Changes of lateral and longitudinal vessel wall thickness with pumping
對比血液體積、血管壁體積和血液與血管壁總體積在不同組織下壓深度L下隨時間變化情況以及在t = 1.66 s時，三者隨組織下壓深度L變化的關系，如圖5所示。可以得出這三種體積在同一時刻的變化趨勢相同。由t = 1.66 s時隨組織下壓深度變化的關系可得，三者相對比，血管壁體積受組織下壓的影響小。橫、縱向血管尺寸（內、外半徑及壁厚）隨著組織下壓深度L變化程度大致相同。
圖5
血液、血管壁隨泵血的體積變化趨勢
Figure5.
The variation trend of blood and vascular wall with the volume of blood pumped
3 基于有限元分析法對組織模型的光學分析
3.1 光學仿真分析流程
由于各介質的吸收系數不同，血管的泵血搏動會引起透射光強度的波動。且血液壓力波動 與透射光強度波動 存在線性關系，即 [7]。本研究在不考慮光散射損失的情況下，依據朗伯-比爾定律模擬均勻單色光源與組織之間的光學吸收作用。光學模型按區域分為血液、血管與組織兩部分來定義參數。其幾何參數是由動力學仿真計算出的每時刻模型血管擴張變化數據定義的，其光學參數根據理論皮膚組織光譜吸收系數得出。在軟件中設定朗伯-比爾定律計算方程以及模型初始邊界條件，包括定義入射強度、將模型外邊定義為透明表面等。該邊界條件使得邊界處不積蓄能量，且進入吸收介質的輻射光束的強度為零。再將模型劃分有限網格進行有限元計算。得出利用朗伯-比爾定律模擬均勻單色光源與組織之間光學作用的結果。
3.2 模型光學參數的簡化及定義
皮膚自表層往里可分為3層——表皮層（epidermis，EP）、真皮層（dermis，DE）、皮下組織層（subcutishypodermis，SH）。真皮層依照含血量分為4類——微循環層（papilmry dermis，PA）、上叢層（superficial arteriovenous plexus，SP）、網狀真皮層（reticular dermis，RD）、深叢層（deep arteriovenous plexus，DP）[23]。模型的血管按尺寸應處于深叢層，研究出射面光強度分布受組織下壓深度L的影響，參考周亮等[23]提出的理論皮膚組織光譜吸收系數，定義血液吸收系數μblood=5 cm?1、組織吸收系數μtunicae = 3 cm?1。設定入射沉積光束功率密度I0 = 5 W/m2，入射方向為垂直向下。入射沉積光束被組織吸收后轉化成熱能很小，故忽略熱能轉換量。模型中的熱能參數，如導熱系數、恒壓熱容、比熱率等均為仿真參考值，對仿真分析出射面光強信號無影響，可不考慮。
3.3 光學分析與結果
先將組織下壓深度L = 0 mm的數據進行處理，得出無組織下壓、峰值時刻（t=1.66 s）時，模型二維截面的光強度分布信息。再對比6種組織下壓深度L = 0、0.1、0.2、0.3、0.4、0.5 mm時，橫坐標直線上的強度分布，如圖6所示。可以看出隨著組織下壓深度L的增大，出射點的光強度逐漸變大，并保留原有的線性趨勢。
圖6
不同組織下壓深度時，模型光能量強度分布圖
Figure6.
Distribution of model light energy intensity under different pressing depths of tissues
由于在實際接收、檢測光強度信號時，只能選取一塊感興趣區域（region of interest，ROI），無法聚焦于某一點。在仿真中，要研究出射面一定區域的出射光強度平均值。以組織下壓深度L = 0 mm模型為例，當出射面以中心區域Lx為 0.2 mm開始，逐漸增大至2.0 mm時，光能量強度平均值與時間t的關系，如圖7所示。
圖7
光能量強度平均值與時間的關系
Figure7.
The relationship between average light intensity and time
提取光強度信號的幅值及峰值，對比選取的橫向出射面區域Lx對信號幅值、峰值的影響，如圖8所示。在單支血管模型中，要選取合適的橫向光出射面區域，獲取的信號幅值才能為最大，選取的區域過大或過小，都會降低信號幅值；選取的橫向光出射面區域越大，其峰值會越大。
圖8
無組織下壓時的信號幅值、峰值
Figure8.
Signal amplitude and peak value without down tissue pressure
模型選取不同組織下壓深度L時，當處于心臟膨脹期（t = 1.48 s）時，血管直徑最小，處于心臟收縮期（t = 1.66 s）時，血管直徑最大。通過處理上述模型的動力學仿真數據，可得出血液直徑波動數據如表3所示。
表3
不同下壓深度時，模型血液直徑數據
Table3.
Data of model blood diameter at different pressing depths
由表3可知，組織下壓深度L = 0 mm模型的血液直徑在0.557～0.703 mm之間波動。橫向光出射面區域Lx > 0.703 mm的幅值點處在圖8信號幅值圖曲線的下降沿處，即當橫向光出射面選取區域大于血液直徑波動最大值時，信號的幅值會隨著橫向光出射面Lx選取區域增大而減小，信號峰值會隨之增大而增大。實際應用中，選取的ROI區域有多支血管分布，選取的單支區域一定大于血管直徑波動最大值。可以推測，實際上在血管組織同等厚度情況下，血管相對密集的區域意味著選取的單支血管橫向光出射面區域Lx更接近血管直徑波動最大值。在此處進行提取光強度信號，幅值相對血管稀疏區（單支血管橫向光出射面大于血管直徑波動最大值）要大一些。
探究當橫向光出射面區域大于血管尺寸（Lx > 0.917 mm）時，信號幅值及峰值隨組織下壓深度L及橫向光出射面區域Lx的變化情況，如圖9所示。結合4幅子圖綜合分析可得，信號幅值及峰值均隨組織下壓深度L的增大而增大。且橫向光出射面區域越小，其幅值增大的效果更明顯。在模型仿真中，選取橫向光出射面區域Lx接近血液直徑波動最大值，且給予適當的下壓深度L，其出射面區域的光能量強度信號幅值和峰值均會增大。在實際透射組織厚度相同情況下，選取血液密集處，給予適當的壓力一定程度上可以提高信號質量。
圖9
不同組織下壓深度、橫向光出射面區域時的信號幅值、峰值
Figure9.
Signal amplitude and peak value of different tissue pressing depths and transverse light emitting surface area
當給予單支血管組織模型下壓深度時，血液區域及組織區域都會產生形變，從而影響出射面區域能量強度信號峰值、幅值。將血液區域的擠壓變形、組織區域的擠壓變薄兩種物理效果進行分離建模，探究血液形變、組織厚度形變二者對信號影響的程度。為貼近實際，僅考慮橫向光出射面區域大于血液直徑波動最大值的情況，即選取橫向光出射面區域Lx > 0.917 mm的數據進行處理分析。將原有模型數據與組織和血液分別受組織下壓深度L影響的模型數據結合在一起，對比不同模型下，組織下壓深度L與信號幅值、峰值之間的關系，如圖10所示，其中a模型為組織和血液均受組織下壓的影響，b模型為僅組織受影響，c模型為僅血液受影響。可以看出，血液區域擠壓形變對信號幅值及峰值的影響均不大，增大信號幅值的關鍵因素是組織區域的形變。實際中，血液和組織均受到下壓效果引起形變，而此時對信號的影響最大。
圖10
不同模型及深度時的信號幅值、峰值
Figure10.
Signal amplitude and peak value at different models and depths
4 結論與展望
本文利用多物理場耦合仿真有限元分析軟件Comsol Multiphysics 5.6（Comsol，瑞典）根據組織血管的基本結構及泵血物理過程，參照郭家松等[17]對指尖甲床動脈的尺寸的研究，定義指端單支小動脈血管組織模型的血管尺寸；依據劉永樂[6]對血管組織動力學上的研究等，完成了模型在泵血周期內產生的幾何形變特征的動力學計算；依據李頂立[7]對單只血容積光電容積脈搏波的數學分析等，完成對模型的光學分析。本研究探究了出射面光強度信號的幅值、峰值受選取橫向光出射面區域和組織下壓深度的影響。從仿真角度得出：在模型組織上施加適當的下壓，會造成橫向血管壁變厚，縱向血管壁變薄；選取血管密集區域相對于血管稀疏區的信號幅值更高，峰值較低；對組織給予適當的下壓深度，會增大出射面光強信號的幅值及峰值。驗證了王焱等[2]提出的適當壓力會改善在動脈血液檢測中的脈搏波信號質量及余江軍等[3]提出的適當壓力對指尖自適應脈搏信號具有改善作用的理論。
本文的力學仿真屬于流體與固體耦合計算，為減少計算量，設定為單向作用耦合，血管壁及組織的形變后的幾何模型不會反過來影響血液流動，即無法分析按壓組織對血液流速的影響。在光學方面，依據朗伯-比爾定律考慮光損耗效果，未考慮光色散情況。以上簡化，一定程度上提高了仿真中的組織下壓操作對信號幅值及峰值的影響效果。在后續的研究中，可以設定雙向耦合物理場，建立多支血管組成血管網模型，來探究血管分布等因素對出射光信號及血液產生的影響。如調整相關設定：設定多支血管不同排列方式、多支血管泵血存在延時、多支血管分布過密集（受按壓后，血管相互抵觸阻礙形變）等。在雙向耦合場中，可以實現對心臟動脈瓣關閉，血液反流等因素形成的潮波、脈搏波現象的仿真；建立指端血管網模型，探究其出射光強度分布，對其ROI區域的提取進行自適應分析。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：由趙鶴童主要負責仿真的全過程及后續數據分析，由周亮、劉朝暉、喬文龍等負責實驗指導，數據分析指導，由喬文龍等負責論文審閱修訂。
