基于電流逆灌注法的經(jīng)絡(luò)可視化研究

王學(xué)民 羅 鳴 馬曉雷 周 鵬 王忠海 馬思宇

(天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072)

引言

經(jīng)絡(luò)是中醫(yī)理論的基本概念之一，是指人體機(jī)能活動(dòng)聯(lián)絡(luò)、調(diào)節(jié)和反應(yīng)的體系。經(jīng)絡(luò)作為中醫(yī)的最重要的內(nèi)容和基礎(chǔ)之一，研究其實(shí)質(zhì)、用現(xiàn)代科學(xué)的方法去揭示其存在的形式和規(guī)律，是中醫(yī)與現(xiàn)代科學(xué)接軌的重要途徑［1］。雖然國(guó)內(nèi)外已有很多學(xué)者在研究經(jīng)絡(luò)的實(shí)質(zhì)，但是還沒(méi)有確鑿、完善的實(shí)驗(yàn)和理論能夠證明和說(shuō)明經(jīng)絡(luò)。經(jīng)絡(luò)的實(shí)質(zhì)是什么?這一問(wèn)題一直困擾醫(yī)學(xué)界，也是生命科學(xué)和醫(yī)學(xué)界迫切需要解決的重要學(xué)術(shù)問(wèn)題。因此，需要借助現(xiàn)代科學(xué)技術(shù)，從各個(gè)方面不斷地研究經(jīng)絡(luò)的現(xiàn)象和特性，將中西醫(yī)更好地結(jié)合起來(lái)，從而更深入地研究經(jīng)絡(luò)的本質(zhì)。

隨著經(jīng)絡(luò)研究的不斷深入，研究者從各個(gè)角度對(duì)經(jīng)絡(luò)進(jìn)行了研究，使人們對(duì)經(jīng)絡(luò)的一些特殊現(xiàn)象有了更進(jìn)深入的認(rèn)識(shí)，如經(jīng)絡(luò)的低流阻通道特性、經(jīng)絡(luò)的循經(jīng)感傳現(xiàn)象、經(jīng)絡(luò)的電阻抗特性等［2－4］。其中，經(jīng)絡(luò)的低流阻通道特性可以表述為:經(jīng)絡(luò)的組織間質(zhì)在沿經(jīng)絡(luò)的方向上，流導(dǎo)較大或流阻較小;在靠近經(jīng)絡(luò)的方向上，流導(dǎo)越來(lái)越大，流阻越來(lái)越小。此假設(shè)也可稱(chēng)為經(jīng)絡(luò)的低流阻假設(shè)。這些現(xiàn)象都是對(duì)經(jīng)絡(luò)這一客觀事實(shí)的進(jìn)一步確認(rèn)，也為經(jīng)絡(luò)物理模型的建立奠定了理論基礎(chǔ)。

本研究在經(jīng)絡(luò)的循經(jīng)感傳現(xiàn)象和經(jīng)絡(luò)的低流阻特性的理論基礎(chǔ)上，首次提出了一種經(jīng)絡(luò)研究的新方法。采用Comsol軟件，建立經(jīng)絡(luò)的簡(jiǎn)化三維物理模型，并用信息流注入的方法，研究信息流在經(jīng)絡(luò)和非經(jīng)絡(luò)上的傳導(dǎo)特性。在此基礎(chǔ)之上，采用動(dòng)態(tài)成像原理及截?cái)嗥娈愔捣纸庹?guī)化算法，在一個(gè)橫截面上進(jìn)行邊界電壓值的測(cè)量，根據(jù)邊界電壓值的差值，重建經(jīng)絡(luò)模型的電導(dǎo)率分布，為經(jīng)絡(luò)的可視化研究提供了有利的科學(xué)依據(jù)。

1 經(jīng)絡(luò)模型建立

研究表明，人體各組織具有不同的電特性［5］，人體各組織的電導(dǎo)率如表1所示。根據(jù)經(jīng)絡(luò)的生物物理特性，經(jīng)絡(luò)作為一種具有低流阻特性的通道，可以運(yùn)行組織液，也可以使組織中的化學(xué)物質(zhì)通過(guò)這一通道進(jìn)行運(yùn)輸和交換，一些物理量(如電流、電磁波、壓力、熱等)也可循著這條通道進(jìn)行傳播?？梢约僭O(shè)經(jīng)絡(luò)對(duì)這些物理量的傳輸效果更明顯，在電特性方面具有較高的電導(dǎo)率。

表1 人體各組織的電導(dǎo)率Tab.1 Conductivity of human tissues

因此，為了更方便地研究經(jīng)絡(luò)，本研究采用Comsol軟件，建立了模擬人體臂部大腸經(jīng)的簡(jiǎn)化物理模型，模型直徑為10 cm、高為15 cm。如圖1(a)所示，圓柱為人體臂部的簡(jiǎn)化模型，小圓柱體為經(jīng)脈的低流阻通道，其余部分為肌肉組織，電導(dǎo)率值分別選擇0.7和1。在圖1(b)中，小圓柱頂部的電極為注入電流激勵(lì)電極，底部的圓圈為圍繞圓柱底部均勻分布的地電極，中間為8個(gè)測(cè)量電極。本研究采用動(dòng)態(tài)成像的方法，分別從經(jīng)絡(luò)口和非經(jīng)絡(luò)口注入電流，根據(jù)兩種方式的差值來(lái)研究經(jīng)絡(luò)的特性參數(shù)。

圖1 人體臂部的簡(jiǎn)化模型。(a)人體臂部簡(jiǎn)化模型;(b)模型電極分布Fig.1 Simplified model of human arm.(a)simplified modelof human arm diagram;(b)electrode distribution model

正問(wèn)題就是已知物場(chǎng)的分布和敏感場(chǎng)的初始及邊界條件，求取電磁場(chǎng)分布。通常采有限元法進(jìn)行求解，本研究采用Comsol軟件進(jìn)行電磁場(chǎng)有限元仿真。圖2為物理模型的剖分圖，通過(guò)Comsol軟件對(duì)所建模型進(jìn)行了三維剖分，以便求出其靈敏度矩陣。圖3為在經(jīng)絡(luò)口注入電流時(shí)的電壓分布，可以很明顯地看到在靠近所設(shè)經(jīng)絡(luò)區(qū)域，電導(dǎo)率明顯比其他區(qū)域要高。

圖2 物理模型的三維剖分Fig.2 Three-dimensional physical model of subdivision

圖3 物理模型測(cè)量截面的電壓分布Fig.3 Physicalmodel of the voltage distribution measured cross section

由于反問(wèn)題的計(jì)算采用基于靈敏度定理的方法，因此需要構(gòu)建靈敏度矩陣J，這也是正問(wèn)題所需解決的一個(gè)重要問(wèn)題。靈敏度矩陣是基于場(chǎng)域內(nèi)電導(dǎo)率的微小變化，而引起邊界測(cè)量電壓變化的原理，實(shí)現(xiàn)非線(xiàn)性問(wèn)題的線(xiàn)性化處理。因此，線(xiàn)性化物理模型可表示為

式中，J為線(xiàn)性化系統(tǒng)矩陣(靈敏度矩陣)，Δσ為在經(jīng)絡(luò)口和非經(jīng)絡(luò)口注入電流時(shí)對(duì)應(yīng)的電導(dǎo)率變化，ΔU為電導(dǎo)率變化引起的邊界電壓變化。

對(duì)于已經(jīng)建立的模型，為了減少問(wèn)題的病態(tài)性，采用方形網(wǎng)格進(jìn)行靈敏度矩陣求解以及逆問(wèn)題計(jì)算。將成像區(qū)域分成一個(gè)個(gè)小區(qū)域，每一個(gè)區(qū)域電導(dǎo)率為常量，經(jīng)過(guò)一系列的推導(dǎo)［6］，可以得到靈敏度矩陣的計(jì)算公式為

式中，Vdm為第d種和第m種類(lèi)型的電流激勵(lì)條件下的電勢(shì)差，σk為第k個(gè)區(qū)域的特征表達(dá)式(電導(dǎo)率)，Ed和Em分別為第 d種和第m種類(lèi)型的電流激勵(lì)條件下的電場(chǎng)強(qiáng)度分布。

由式(2)獲得了該模型的靈敏度矩陣，為經(jīng)絡(luò)模型反問(wèn)題的求解奠定了基礎(chǔ)。

2 TSVD正則化方法

根據(jù)邊界電壓的測(cè)量值，通過(guò)適當(dāng)?shù)乃惴ǎ蟮媚Ｐ蛢?nèi)的阻抗分布或阻抗變化，這在電磁場(chǎng)分析中被稱(chēng)為逆問(wèn)題，而對(duì)經(jīng)絡(luò)特性參數(shù)的圖像重建就是對(duì)反問(wèn)題的求解。對(duì)于經(jīng)絡(luò)特性參數(shù)的重建，需要借助一定的先驗(yàn)知識(shí)，主要包括經(jīng)絡(luò)組織的結(jié)構(gòu)特征以及一般組織的電導(dǎo)率分布信息。

通過(guò)計(jì)算得到靈敏度矩陣J，分析這個(gè)矩陣的相關(guān)參數(shù)，可知它的條件數(shù)為 cond(A)=1.48×1028很大，因此矩陣J是一個(gè)病態(tài)矩陣，求解式(1)是一個(gè)不適定問(wèn)題。奇異值分解是分析不適定問(wèn)題的一種非常有效的方法，而截?cái)嗥娈愔捣纸?TSVD)正則化方法已經(jīng)成功解決了許多不適定問(wèn)題，因此本研究采用截?cái)嗥娈愔捣纸獾恼齽t化方法。

根據(jù)所求的靈敏度矩陣，采用截?cái)嗥娈愔档恼齽t化方法［7－10］，用一個(gè)良態(tài)的虧秩矩陣 B來(lái)近似A。所謂良態(tài)的虧秩矩陣，是指這個(gè)矩陣的條件數(shù)

不是很大。

定義:設(shè) A∈ Rm×n的奇異值分解由 A =給出，給定 ε ，取 k(ε)滿(mǎn)足

即 σ1≥ … ≥ σk(ε)≥ σk(ε)+1≥ … ≥ σn，則

稱(chēng)為A的截?cái)郤VD分解(truncated singular value decomposition，TSVD)，k(ε)為截?cái)鄥?shù)。在不引起混淆的情況下，記Ak(ε)為Ak，實(shí)質(zhì)上是通過(guò)控制ε濾掉對(duì)矩陣性質(zhì)影響較大的小奇異值。當(dāng)A是從有噪聲的觀測(cè)數(shù)據(jù)得到時(shí)，這種方法能夠起到很好的作用。因此，可以將問(wèn)題轉(zhuǎn)化為求解

其最小二乘解為

這就是TSVD正則化方法。

為了確定k(ε)，采用范數(shù)比的方法。

3 結(jié)果

為了更好地研究靈敏度矩陣的性質(zhì)和算法的比較，先假設(shè)經(jīng)絡(luò)在圓形場(chǎng)域的中間，其電導(dǎo)率的分布如圖4所示，中心部分的電導(dǎo)率高于周?chē)糠值碾妼?dǎo)率。采用動(dòng)態(tài)成像原理，分別從經(jīng)絡(luò)口和非經(jīng)絡(luò)口(靠近經(jīng)絡(luò)處)注入頻率為100 Hz、有效值為10 mA的激勵(lì)電流，通過(guò)截?cái)嗥娈愔捣纸庹?guī)化算法，重建模型電導(dǎo)率分布。

圖4 模型電導(dǎo)率分布Fig.4 Conductivity distribution model

根據(jù)范數(shù)比的定義計(jì)算范數(shù)比，有

由于截?cái)鄥?shù)與噪聲是相關(guān)的，在不同的參數(shù)下相關(guān)系數(shù)是不一樣的，因此選擇υ(k)≥α=0.999 9的最小整數(shù) k=19作為截?cái)鄥?shù)，成像效果如圖5所示。

圖5 k=19時(shí)的圖像重建結(jié)果Fig.5 The image reconstruction when k=19

然后，在計(jì)算得到的測(cè)量數(shù)據(jù)中分別加入1%和5%的隨機(jī)噪聲，成像效果如圖6所示。圖6(a)是在所測(cè)得的數(shù)據(jù)中加入1%隨機(jī)噪聲后重建所得，圖6(b)是在所測(cè)得的數(shù)據(jù)中加入5%隨機(jī)噪聲后重建所得，兩幅圖的電導(dǎo)率分布基本沒(méi)有變化，且圓心部位的電導(dǎo)率都明顯比其他區(qū)域的電導(dǎo)率要高?？梢?jiàn)，該截?cái)鄥?shù)能夠起到較好的圖像重建效果，使靈敏度矩陣J的條件數(shù)減小，起到抑制噪聲的作用。

圖6 加入不同水平噪聲時(shí)的圖像重建結(jié)果。(a)1%;(b)5%Fig.6 Image reconstruction by adding different random noises.(a)1%;(b)5%

對(duì)于經(jīng)絡(luò)特性的研究，引用一些先驗(yàn)知識(shí)(包括內(nèi)部組織的結(jié)構(gòu)特性和組織的電導(dǎo)率信息)，建立了三維的經(jīng)絡(luò)物理模型，圖7為經(jīng)絡(luò)物理模型的原始電導(dǎo)率分布。采用電流激勵(lì)-電壓測(cè)量的方式，分別從經(jīng)絡(luò)口和非經(jīng)絡(luò)口注入電流，測(cè)量?jī)煞N方式的電壓值，通過(guò)差值研究經(jīng)絡(luò)的特性參數(shù)。根據(jù)動(dòng)態(tài)成像的方法，采用截?cái)嗥娈愔捣纸馑惴?，圖像重建效果如圖8所示，可以明顯看到經(jīng)絡(luò)部分高于周?chē)渌糠值内厔?shì)。

4 結(jié)論

圖7 經(jīng)絡(luò)物理模型原始電導(dǎo)率分布Fig.7 Physicalmodelofthe originalconductivity distribution Meridian

圖8 重建經(jīng)絡(luò)模型電導(dǎo)率分布Fig.8 Meridian reconstructed conductivity distribution model

本研究在經(jīng)絡(luò)的低流阻特性和生物物理特性的基礎(chǔ)上，合理地建立了人體臂部的簡(jiǎn)化數(shù)學(xué)模型。然后，提出了一種經(jīng)絡(luò)研究的新方法——縱向注入電流、在一個(gè)橫截面上進(jìn)行邊界電壓的測(cè)量。采用在經(jīng)絡(luò)口和非經(jīng)絡(luò)口注入電流的方法，研究電壓分布情況的差異。采用有限元法和Comsol軟件進(jìn)行了電磁場(chǎng)的仿真，得出了靈敏度矩陣公式。最后，采用動(dòng)態(tài)成像方法和截?cái)嗥娈愔捣纸獾恼齽t化算法，基于先驗(yàn)知識(shí)重建了經(jīng)絡(luò)模型的電特性分布(電導(dǎo)率)，為進(jìn)一步研究經(jīng)絡(luò)的本質(zhì)提供了新的方法。

［1］費(fèi)倫，魏瑚.經(jīng)絡(luò)學(xué)說(shuō)基礎(chǔ)理論的構(gòu)建及其學(xué)術(shù)地位的確立［J］.世界科學(xué)，2005，9:43－45.

［2］劉芳，黃光英.經(jīng)絡(luò)的生物物理學(xué)、化學(xué)特性的研究進(jìn)展［J］.針刺研究，2007，32(4):281－285.

［3］楊威生.低阻經(jīng)絡(luò)研究Ⅲ 對(duì)經(jīng)絡(luò)組織學(xué)本質(zhì)的推斷［J］.北京大學(xué)學(xué)報(bào)(自然科學(xué)版)，2008，44(2):227－280.

［4］王永紅.經(jīng)絡(luò)的實(shí)質(zhì)探析［J］.世界中醫(yī)藥，2010，5(4):236－237.

［5］王研.EIT固定參考點(diǎn)圖像重建和電極結(jié)構(gòu)優(yōu)化方法研究［D］.北京:中國(guó)協(xié)和醫(yī)科大學(xué)，2002:1－5.

［6］Polydorides N，Lionheart WRB.A Matlab toolkit for threedimensional electrical impedance tomography:a contribution to the electricalimpedance and diffuse opticalreconstruction software project［J］.Measurement Science and Technology，2002，13:1871－1883.

［7］聶守平，魏曉燕.數(shù)字圖像的奇異值分解［J］.南京師大學(xué)報(bào)，2001，24(1):59－61.

［8］黃小為，吳傳生，朱華平，基于奇異值分解建立的一種新的正則化方法［J］.數(shù)學(xué)物理學(xué)報(bào)，2005，25A(3):331－336.

［9］黃小為，吳傳生，李卓球.TSVD正則化方法的參數(shù)選取及數(shù)值計(jì)算［J］.華中師范大學(xué)學(xué)報(bào)，2006，40(2):154－157.

［10］董超峰，孫方裕.基本解方法求解各向異性材料中熱傳導(dǎo)方程的時(shí)間反向問(wèn)題［J］.浙江大學(xué)學(xué)報(bào)，2007，34(1):33－39.