針對肌肉動態收縮情況下的高密度表面肌電(sEMG)信號，本文提出了一種基于空間位置的sEMG信號分解方法。首先，根據肌肉運動單元(MU)在各個通道上的波形相關性，提取發放時刻，然后利用肌肉MU的空間位置分類發放時刻，最后得到MU發放序列。仿真結果表明，分類后單個MU發放序列準確率大于91.67%。針對實際sEMG信號，通過“二源法”找到同一個MU發放序列的準確率達到（88.3 ± 2.1）%以上。本文為動態sEMG信號分解提供了一種新思路。
引用本文： 何金保, 管冰蕾, 黃凱, 駱再飛. 一種肌肉動態收縮的高密度表面肌電信號分解新方法. 生物醫學工程學雜志, 2021, 38(6): 1081-1086. doi: 10.7507/1001-5515.202104020 復制
引言
表面肌電信號（surface electromyogram，sEMG）是由肌肉運動單元（motor unit，MU）產生的動作電位在皮膚表面疊加形成的，MU的發放頻率、波形幅值等參數反映了神經肌肉系統狀態，通過MU信息有助于理解神經肌肉工作機制和控制策略。與針電極等插入式電極相比，sEMG信號有著無創、簡單的優點。因此，高密度sEMG信號在臨床醫學、人機功效學、康復醫學以及體育科學等方面均有重要的實用價值[1]。
為了從高密度sEMG信號中獲取MU信息，國內外學者提出了不同的分解方法[2-5]，比較典型有基于模板匹配的模式識別方法和“盲源分離”方法。模式識別方法對采樣頻率要求較高，而且使用的電極數量較少，不利于獲取MU空間信息。“盲源分離”方法應用于高密度sEMG信號，更有利于提取MU，并且能得到MU空間信息。雖然模式識別方法和“盲源分離”方法已經比較成熟，但是它們主要應用于肌肉靜態收縮，即肌肉恒力收縮的情況。
由于肌肉動態收縮時，肌肉力時變、肌纖維長度、傳導速度等參數會變化，導致MU的波形畸變，應用于肌肉靜態收縮的分解方法不能直接應用于肌肉動態收縮情況。目前有少數學者對肌肉動態收縮sEMG信號進行了研究，Glaser等[6]提出卷積核分解方法，該方法將sEMG信號分段，采用線性函數得到發放時刻，但是發放時刻分類仍然存在問題。Glaser等[7]還提出了基于MU波形非傳輸分量的分解算法，但由于傳輸分量與非傳輸分量的分離存在難度，實現比較困難。DeLuca等[8]使用機器學習算法分解sEMG信號，但是MU波形變化影響因素眾多，難以預測，所以識別仍然較困難。
綜上所述，針對動態肌肉收縮sEMG信號分解問題，本文提出基于MU空間位置的分解方法，以仿真和實際的sEMG信號驗證本文方法，以期該方法可有效實現肌肉動態收縮sEMG 信號的分解。本文的方法為神經肌肉系統的生理研究提供了新的思路，有助于深入了解神經肌肉的控制機理，同時也為人機交互系統、臨床診斷等提供了有效的工具。
1 研究方法
1.1 構建sEMG信號模型
肌肉動態收縮情況下，高密度sEMG信號的數學模型方程如式（1）所示[3-4]：
 |
其中，X(n)為高密度sEMG信號向量，e(n)為零均值白噪聲向量，
(n)是N個源序列
的擴展形式，H是肌肉動態收縮下混合矩陣，它由
長度為P個采樣點的通道響應所組成，?H是肌肉動態收縮下的干擾矩陣。本文將?H
(n)看作噪聲干擾，采用“盲源分離”方法提取發放時刻。
1.2 分解方法
將式(1)兩側同時乘以向量f，則可以得到如式(2)所示：
 |
采用變換矩陣Xc(n)替換X(n)，在不考慮噪聲影響的情況下，第i個MU的發放序列τi可近似表示為如式(3)所示:
 |
其中，Xc(n)是與矩陣X(n)有較強列相關性的矩陣，保留sEMG信號中的時域特征，f是第i個MU的發放時刻。Xc(n)可以采用奇異值分解(singular value decomposition, SVD)獲得，即用SVD方法得到的右酉矩陣構造，式(3)可以寫為如式（4）所示：
 |
其中，矩陣U由SVD方法得到的右酉矩陣，即X(n) = VDUT，V、U分別是左、右酉矩陣，D矩陣除了主對角線上的元素以外全為0，主對角線上的每個元素都稱為奇異值。矩陣U相對于原始信號矩陣X(n)，具有較強的穩定性和抗噪性，既增強了信號的時域特征，又降低了計算負荷[9-16]。
式（4）得到的發放序列τi中會存在較多干擾，很可能會出現不同MU的發放時刻在同一個序列中，因此，正確的發放時刻分類十分重要。本論文采用基于MU的空間位置來分類MU發放序列，首先根據發放時刻的雙極性波形確定MU在表面電極上的平面位置，然后根據發放時刻波形的幅值擬合確定MU的深度，最后通過MU平面位置和深度實現該發放時刻的分類。
MU平面位置的確定是利用MU在終板區兩側傳輸過程中雙極性波形互為相反的特點，從而得到MU平面位置。
MU的深度可由表面電極上檢測的發放波形幅值估計[5, 11-13]，將肢體當作圓柱體得到如式（5）、式（6）所示關系：
 |
 |
其中，h、Amp分別對應MU深度、MU波形幅值，Lth是肌纖維垂直方向上的Amp所對應電極到(1/2)Amp所對應電極的弧長，α、β、d、γ均為常數，Lth的值通過擬合方法得到。綜上，動態分解算法流程如圖1所示。
圖1
動態分解算法流程圖
Figure1.
Flow diagram of the dynamic decomposition algorithm
2 仿真驗證
2.1 仿真模型
纖維動作電位是sEMG信號的基礎，sEMG信號就是大量纖維動電位在皮膚表面疊加形成的。本文從纖維動作電位角度建立模型，生成sEMG信號。
纖維動作電位采用三極子代替[11,13]，如式（7）～式（8）所示：
 |
 |
其中，P1、P2、P3表示三極子，a、b是三極子間距離。
兩組三極子從終板區向肌纖維兩側傳播[11,13]，將肢體看作圓柱體，以皮膚表面作為x-z平面，第i根纖維在電極上電壓 
，如式(9)所示：
 |
其中，
，
是x、y方向上的電導率，
是z方向上的電導率，Pi表示第i個三極子，
表示三極子坐標，
表示表面電極坐標。上述模型較好地反映了各電極間電位相關性，是目前常用的生成sEMG信號的模型。
2.2 仿真結果
設計皮膚表面電極個數為11 × 11，每個電極間間隔9 mm。MU個數20，每個MU肌纖維數量30～50，MU深度在2～10 mm范圍內隨機分布，且分布在電極覆蓋范圍內，MU終板區分布在 ± 1 mm范圍內，肌纖維長度100～170 mm變化，發放頻率小于25 Hz且隨機，傳導速度4.0 m/s，20個MU參與發放，三極子電量參數為
在[18～28]范圍內變化，a、b的值分別取為2.1和7，電導率 
為0.063，K為6，手臂半徑為70 mm。仿真信號中增加均值為零的高斯白噪聲，三種信噪比(signal noise ratio, SNR)分別為5、15、25 dB，采樣頻率2 048 Hz。生成的其中某一通道sEMG信號及其放大如圖2所示。
圖2
仿真信號及放大
Figure2.
Simulated signal and its amplification
準確率（true pulse ratio，TPR）是分解得到的MU發放時刻與相應原始MU發放時刻的比值。計算仿真信號在不同SNR情況下，能找到所有MU的TPR：在5 dB時找到6個MU，TPR = 165/180 × 100% = 91.67%；15 dB時找到12個MU，TPR = 335/360 × 100% = 93.06%；25 dB時找到18個MU，TPR = 510/540 × 100% = 94.44%。具體每個MU的TPR如圖3所示。
圖3
每個MU的發放時刻準確率
Figure3.
Accuracy rate of the firing time of each MU
為了驗證有效性，仿真時，MU發放完全隨機，因為肌肉動態收縮時，MU的發放沒有規律性，所以隨機設置發放頻率是合理的。從分解效果來看，本文的算法可以應對MU隨機發放的情況。
如圖4所示是不同深度MU在皮膚表面得到的擬合弧長，從圖4可以看出，MU在不同深度時，最大幅值處到一半幅值處的弧長在1.2～6.0 mm之間，因此可以利用弧長來分類發放時刻。
圖4
不同深度MU的擬合弧長
Figure4.
The fitting arc length of MU with different depth
3 試驗驗證
3.1 信號采集及驗證方法
肌肉動態收縮試驗征集了8個健康人作為受試者，其中女性3名（年齡35～42歲,體重52～56 kg），男性5名（年齡分別為27～39歲,體重60～72 Kg）。在寧波市第一醫院采集數據，通過該院倫理審查，試驗前，告知受試者試驗過程及注意事項，簽訂有關受試者知情同意書。
試驗采用64通道的高密度sEMG信號測量儀（Refa64，TMSi BV，荷蘭）采集sEMG信號，參考電極置于手腕，采集頻率2 048 Hz，設置帶通濾波頻率10～500 Hz。64通道（8 × 8）高密度電極陣列極間距8.5 mm，電極直徑5 mm。試驗前需要清除右手肱二頭肌皮膚上角質層，把電極陣列貼在肱二頭肌肌腹部位，保證電極與皮膚的接觸良好。測試時，受試者自然站立，右手握1 kg啞鈴，從自然垂下位置向上作運動并反復數次。
為了驗證本文分解真實sEMG信號方法的可靠性，采用“二源法”驗證。將64 通道的信號按照電極序號的奇偶序號分成兩組，按照電極序號的奇偶序號來分，即第一組信號包括電極1、3、5、···、61、63上采集到的信號，第二組信號包括電極2、4、6、···、62、64上采集到的信號，然后利用本文的方法分別對兩組信號進行分解對比。
3.2 試驗結果及分析
試驗采集到的某一通道信號及其放大如圖5所示。
圖5
試驗采集的信號及放大
Figure5.
Real signal and its amplification
利用本文的方法對真實sEMG信號進行分解，取采集信號長度為3 s，8個受試者的分解結果如表1所示。式（3）中直接使用原始信號，即“未采用SVD分解方法”，8個受試者分別得到的11～13個MU；采用本文方法，8個受試者分別得到11～15個MU，并且受試者1、4、7、8找到2個平面位置相同、深度不同的MU。從表1中“未采用SVD分解方法”和“本文方法”找到的“MU個數”可以看出，本文方法采用SVD預處理以后，分解得到的MU個數稍多，也就是說，SVD預處理增強了sEMG信號的時域特征，這與其他學者的研究結論相符[14, 16]，并且本文方法利用MU空間位置分類發放時刻，可以找到平面位置重疊的MU，見表1中受試者1、4、7、8，這個是其它分類方法無法做到的。
表1
真實sEMG信號的分解結果
Table1.
Decomposition results of real sEMG signals
另外采用本文的方法，針對8個受試者，“二源法”找到的MU發放時刻與全部通道信號找到的MU發放時刻比例為（88.3 ± 2.1）%，該結果與文獻[8]中手臂肘關節屈曲/伸展運動時，對高密度sEMG信號分解得到結果基本相符。由于肌肉動態收縮的MU波形變化，文獻[8]采用的是基于波形的機器學習方法，所以該方法受到運動時MU波形變化影響相對較大。
如圖6所示是受試者1采用二源法找到的MU發放時刻對比圖，共13個MU。其中第1個MU和第13個MU的平面位置是重疊的，從弧長的角度分類得到，第1個MU擬合弧長是1.6 mm，第13個MU擬合弧長2.3 mm。從圖6可以看出，在肌肉動態收縮過程中，MU的發放頻率不規律，變化較大。由于篇幅限制，受試者2～受試者8的MU發放時刻對比圖不再給出。
圖6
“二源法”發放時刻對比(受試者1)
Figure6.
Comparison of firings with “two-source method”(subject 1)
受試者1采用“二源法”得到13個MU，兩組信號找到的相同MU發放時刻與全部通道信號找到的MU發放時刻比例為86%。通常從奇通道和偶通道找到MU的發放時刻比全部通道的要少，這是由于通道數的減少，使得分解得到的發放時刻有所減少，這與通道數越多得到的發放時刻越多結論相吻合。從已經找到的發放時刻來看，驗證了本文分解方法的可靠性和準確性。
4 結論
本文提出的針對肌肉動態收縮的sEMG信號分解方法，首先對原始sEMG信號進行SVD分解預處理，然后根據MU在各個通道上的波形相關性提取發放時刻，最后通過MU空間位置，對發放時刻進行分類。仿真及試驗結果表明本文的方法獲得了較好的分解性能，該方法可以應用在對動態sEMG 信號分解的相關研究中。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
致謝：衷心感謝美國休斯頓大學生物醫學工程系張迎春教授團隊的大力協助！
引言
表面肌電信號（surface electromyogram，sEMG）是由肌肉運動單元（motor unit，MU）產生的動作電位在皮膚表面疊加形成的，MU的發放頻率、波形幅值等參數反映了神經肌肉系統狀態，通過MU信息有助于理解神經肌肉工作機制和控制策略。與針電極等插入式電極相比，sEMG信號有著無創、簡單的優點。因此，高密度sEMG信號在臨床醫學、人機功效學、康復醫學以及體育科學等方面均有重要的實用價值[1]。
為了從高密度sEMG信號中獲取MU信息，國內外學者提出了不同的分解方法[2-5]，比較典型有基于模板匹配的模式識別方法和“盲源分離”方法。模式識別方法對采樣頻率要求較高，而且使用的電極數量較少，不利于獲取MU空間信息。“盲源分離”方法應用于高密度sEMG信號，更有利于提取MU，并且能得到MU空間信息。雖然模式識別方法和“盲源分離”方法已經比較成熟，但是它們主要應用于肌肉靜態收縮，即肌肉恒力收縮的情況。
由于肌肉動態收縮時，肌肉力時變、肌纖維長度、傳導速度等參數會變化，導致MU的波形畸變，應用于肌肉靜態收縮的分解方法不能直接應用于肌肉動態收縮情況。目前有少數學者對肌肉動態收縮sEMG信號進行了研究，Glaser等[6]提出卷積核分解方法，該方法將sEMG信號分段，采用線性函數得到發放時刻，但是發放時刻分類仍然存在問題。Glaser等[7]還提出了基于MU波形非傳輸分量的分解算法，但由于傳輸分量與非傳輸分量的分離存在難度，實現比較困難。DeLuca等[8]使用機器學習算法分解sEMG信號，但是MU波形變化影響因素眾多，難以預測，所以識別仍然較困難。
綜上所述，針對動態肌肉收縮sEMG信號分解問題，本文提出基于MU空間位置的分解方法，以仿真和實際的sEMG信號驗證本文方法，以期該方法可有效實現肌肉動態收縮sEMG 信號的分解。本文的方法為神經肌肉系統的生理研究提供了新的思路，有助于深入了解神經肌肉的控制機理，同時也為人機交互系統、臨床診斷等提供了有效的工具。
1 研究方法
1.1 構建sEMG信號模型
肌肉動態收縮情況下，高密度sEMG信號的數學模型方程如式（1）所示[3-4]：
|
其中，X(n)為高密度sEMG信號向量，e(n)為零均值白噪聲向量，(n)是N個源序列的擴展形式，H是肌肉動態收縮下混合矩陣，它由長度為P個采樣點的通道響應所組成，?H是肌肉動態收縮下的干擾矩陣。本文將?H(n)看作噪聲干擾，采用“盲源分離”方法提取發放時刻。
1.2 分解方法
將式(1)兩側同時乘以向量f，則可以得到如式(2)所示：
|
采用變換矩陣Xc(n)替換X(n)，在不考慮噪聲影響的情況下，第i個MU的發放序列τi可近似表示為如式(3)所示:
|
其中，Xc(n)是與矩陣X(n)有較強列相關性的矩陣，保留sEMG信號中的時域特征，f是第i個MU的發放時刻。Xc(n)可以采用奇異值分解(singular value decomposition, SVD)獲得，即用SVD方法得到的右酉矩陣構造，式(3)可以寫為如式（4）所示：
|
其中，矩陣U由SVD方法得到的右酉矩陣，即X(n) = VDUT，V、U分別是左、右酉矩陣，D矩陣除了主對角線上的元素以外全為0，主對角線上的每個元素都稱為奇異值。矩陣U相對于原始信號矩陣X(n)，具有較強的穩定性和抗噪性，既增強了信號的時域特征，又降低了計算負荷[9-16]。
式（4）得到的發放序列τi中會存在較多干擾，很可能會出現不同MU的發放時刻在同一個序列中，因此，正確的發放時刻分類十分重要。本論文采用基于MU的空間位置來分類MU發放序列，首先根據發放時刻的雙極性波形確定MU在表面電極上的平面位置，然后根據發放時刻波形的幅值擬合確定MU的深度，最后通過MU平面位置和深度實現該發放時刻的分類。
MU平面位置的確定是利用MU在終板區兩側傳輸過程中雙極性波形互為相反的特點，從而得到MU平面位置。
MU的深度可由表面電極上檢測的發放波形幅值估計[5, 11-13]，將肢體當作圓柱體得到如式（5）、式（6）所示關系：
|
|
其中，h、Amp分別對應MU深度、MU波形幅值，Lth是肌纖維垂直方向上的Amp所對應電極到(1/2)Amp所對應電極的弧長，α、β、d、γ均為常數，Lth的值通過擬合方法得到。綜上，動態分解算法流程如圖1所示。
圖1
動態分解算法流程圖
Figure1.
Flow diagram of the dynamic decomposition algorithm
2 仿真驗證
2.1 仿真模型
纖維動作電位是sEMG信號的基礎，sEMG信號就是大量纖維動電位在皮膚表面疊加形成的。本文從纖維動作電位角度建立模型，生成sEMG信號。
纖維動作電位采用三極子代替[11,13]，如式（7）～式（8）所示：
|
|
其中，P1、P2、P3表示三極子，a、b是三極子間距離。
兩組三極子從終板區向肌纖維兩側傳播[11,13]，將肢體看作圓柱體，以皮膚表面作為x-z平面，第i根纖維在電極上電壓 ，如式(9)所示：
|
其中，，是x、y方向上的電導率，是z方向上的電導率，Pi表示第i個三極子，表示三極子坐標，表示表面電極坐標。上述模型較好地反映了各電極間電位相關性，是目前常用的生成sEMG信號的模型。
2.2 仿真結果
設計皮膚表面電極個數為11 × 11，每個電極間間隔9 mm。MU個數20，每個MU肌纖維數量30～50，MU深度在2～10 mm范圍內隨機分布，且分布在電極覆蓋范圍內，MU終板區分布在 ± 1 mm范圍內，肌纖維長度100～170 mm變化，發放頻率小于25 Hz且隨機，傳導速度4.0 m/s，20個MU參與發放，三極子電量參數為在[18～28]范圍內變化，a、b的值分別取為2.1和7，電導率 為0.063，K為6，手臂半徑為70 mm。仿真信號中增加均值為零的高斯白噪聲，三種信噪比(signal noise ratio, SNR)分別為5、15、25 dB，采樣頻率2 048 Hz。生成的其中某一通道sEMG信號及其放大如圖2所示。
圖2
仿真信號及放大
Figure2.
Simulated signal and its amplification
準確率（true pulse ratio，TPR）是分解得到的MU發放時刻與相應原始MU發放時刻的比值。計算仿真信號在不同SNR情況下，能找到所有MU的TPR：在5 dB時找到6個MU，TPR = 165/180 × 100% = 91.67%；15 dB時找到12個MU，TPR = 335/360 × 100% = 93.06%；25 dB時找到18個MU，TPR = 510/540 × 100% = 94.44%。具體每個MU的TPR如圖3所示。
圖3
每個MU的發放時刻準確率
Figure3.
Accuracy rate of the firing time of each MU
為了驗證有效性，仿真時，MU發放完全隨機，因為肌肉動態收縮時，MU的發放沒有規律性，所以隨機設置發放頻率是合理的。從分解效果來看，本文的算法可以應對MU隨機發放的情況。
如圖4所示是不同深度MU在皮膚表面得到的擬合弧長，從圖4可以看出，MU在不同深度時，最大幅值處到一半幅值處的弧長在1.2～6.0 mm之間，因此可以利用弧長來分類發放時刻。
圖4
不同深度MU的擬合弧長
Figure4.
The fitting arc length of MU with different depth
3 試驗驗證
3.1 信號采集及驗證方法
肌肉動態收縮試驗征集了8個健康人作為受試者，其中女性3名（年齡35～42歲,體重52～56 kg），男性5名（年齡分別為27～39歲,體重60～72 Kg）。在寧波市第一醫院采集數據，通過該院倫理審查，試驗前，告知受試者試驗過程及注意事項，簽訂有關受試者知情同意書。
試驗采用64通道的高密度sEMG信號測量儀（Refa64，TMSi BV，荷蘭）采集sEMG信號，參考電極置于手腕，采集頻率2 048 Hz，設置帶通濾波頻率10～500 Hz。64通道（8 × 8）高密度電極陣列極間距8.5 mm，電極直徑5 mm。試驗前需要清除右手肱二頭肌皮膚上角質層，把電極陣列貼在肱二頭肌肌腹部位，保證電極與皮膚的接觸良好。測試時，受試者自然站立，右手握1 kg啞鈴，從自然垂下位置向上作運動并反復數次。
為了驗證本文分解真實sEMG信號方法的可靠性，采用“二源法”驗證。將64 通道的信號按照電極序號的奇偶序號分成兩組，按照電極序號的奇偶序號來分，即第一組信號包括電極1、3、5、···、61、63上采集到的信號，第二組信號包括電極2、4、6、···、62、64上采集到的信號，然后利用本文的方法分別對兩組信號進行分解對比。
3.2 試驗結果及分析
試驗采集到的某一通道信號及其放大如圖5所示。
圖5
試驗采集的信號及放大
Figure5.
Real signal and its amplification
利用本文的方法對真實sEMG信號進行分解，取采集信號長度為3 s，8個受試者的分解結果如表1所示。式（3）中直接使用原始信號，即“未采用SVD分解方法”，8個受試者分別得到的11～13個MU；采用本文方法，8個受試者分別得到11～15個MU，并且受試者1、4、7、8找到2個平面位置相同、深度不同的MU。從表1中“未采用SVD分解方法”和“本文方法”找到的“MU個數”可以看出，本文方法采用SVD預處理以后，分解得到的MU個數稍多，也就是說，SVD預處理增強了sEMG信號的時域特征，這與其他學者的研究結論相符[14, 16]，并且本文方法利用MU空間位置分類發放時刻，可以找到平面位置重疊的MU，見表1中受試者1、4、7、8，這個是其它分類方法無法做到的。
表1
真實sEMG信號的分解結果
Table1.
Decomposition results of real sEMG signals
另外采用本文的方法，針對8個受試者，“二源法”找到的MU發放時刻與全部通道信號找到的MU發放時刻比例為（88.3 ± 2.1）%，該結果與文獻[8]中手臂肘關節屈曲/伸展運動時，對高密度sEMG信號分解得到結果基本相符。由于肌肉動態收縮的MU波形變化，文獻[8]采用的是基于波形的機器學習方法，所以該方法受到運動時MU波形變化影響相對較大。
如圖6所示是受試者1采用二源法找到的MU發放時刻對比圖，共13個MU。其中第1個MU和第13個MU的平面位置是重疊的，從弧長的角度分類得到，第1個MU擬合弧長是1.6 mm，第13個MU擬合弧長2.3 mm。從圖6可以看出，在肌肉動態收縮過程中，MU的發放頻率不規律，變化較大。由于篇幅限制，受試者2～受試者8的MU發放時刻對比圖不再給出。
圖6
“二源法”發放時刻對比(受試者1)
Figure6.
Comparison of firings with “two-source method”(subject 1)
受試者1采用“二源法”得到13個MU，兩組信號找到的相同MU發放時刻與全部通道信號找到的MU發放時刻比例為86%。通常從奇通道和偶通道找到MU的發放時刻比全部通道的要少，這是由于通道數的減少，使得分解得到的發放時刻有所減少，這與通道數越多得到的發放時刻越多結論相吻合。從已經找到的發放時刻來看，驗證了本文分解方法的可靠性和準確性。
4 結論
本文提出的針對肌肉動態收縮的sEMG信號分解方法，首先對原始sEMG信號進行SVD分解預處理，然后根據MU在各個通道上的波形相關性提取發放時刻，最后通過MU空間位置，對發放時刻進行分類。仿真及試驗結果表明本文的方法獲得了較好的分解性能，該方法可以應用在對動態sEMG 信號分解的相關研究中。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
致謝：衷心感謝美國休斯頓大學生物醫學工程系張迎春教授團隊的大力協助！
