微波熱療溫度場(chǎng)仿真中SAR分布影響因素的研究

張曼，吳水才，胡芮，林嵐，高宏建

北京工業(yè)大學(xué) 生命科學(xué)與生物工程學(xué)院，北京 100124

微波熱療溫度場(chǎng)仿真中SAR分布影響因素的研究

張曼，吳水才，胡芮，林嵐，高宏建

北京工業(yè)大學(xué) 生命科學(xué)與生物工程學(xué)院，北京 100124

為了提高微波熱療溫度場(chǎng)仿真的精度，本文從以下3個(gè)方面研究了影響比吸收率（SAR）分布的因素：對(duì)實(shí)驗(yàn)中的溫度誤差做定量計(jì)算；確定擬合中的最佳時(shí)間段；提出一種SAR分布的優(yōu)化擬合方法。經(jīng)過(guò)大量計(jì)算與對(duì)比仿真結(jié)果后發(fā)現(xiàn)，當(dāng)實(shí)驗(yàn)中溫度誤差控制在3.5 ℃以內(nèi)，選取前60 s、70 s…110 s時(shí)間段進(jìn)行擬合對(duì)仿真結(jié)果影響不大，所提出的SAR加權(quán)最小二乘擬合法更加穩(wěn)定有效。

微波熱療；溫度場(chǎng)仿真；SAR擬合；加權(quán)最小二乘法

惡性腫瘤嚴(yán)重威脅著人類的健康與生命。根據(jù)世界衛(wèi)生組織對(duì)全世界死亡者病因的統(tǒng)計(jì)， 癌癥已居首位[1]。目前，手術(shù)切除是治療癌癥的首選方法，但術(shù)后創(chuàng)傷大，患者恢復(fù)慢，且病灶極易轉(zhuǎn)移。微波熱療作為一種微創(chuàng)治療腫瘤的方法，可在病灶原位進(jìn)行加熱以滅活腫瘤組織[2]，能夠克服手術(shù)治療的缺陷。

由于腫瘤的個(gè)體差異性很大，需要建立精確的微波熱療三維溫度場(chǎng)仿真來(lái)確定微波熱療手術(shù)中的最佳加熱參數(shù)[3]。該溫度場(chǎng)是生物組織吸收微波能量后所轉(zhuǎn)化的，所以建立精確溫度場(chǎng)分布的關(guān)鍵是準(zhǔn)確計(jì)算微波能量在人體組織內(nèi)的比吸收率（Specific Absorption Rate，SAR）分布。目前，對(duì)于不同類型的微波天線SAR分布的研究有很多[4-9]，主要方法有3種：第一種是根據(jù)Maxwell方程[10]求出輻射電磁場(chǎng)，通過(guò)理論推導(dǎo)和數(shù)值計(jì)算得到SAR分布；第二種是翟偉明[11]等報(bào)道的通過(guò)半經(jīng)驗(yàn)公式推導(dǎo)與測(cè)量數(shù)據(jù)修正相結(jié)合的方法來(lái)確定SAR分布；第三種是根據(jù)實(shí)驗(yàn)測(cè)量[12]的方法得到SAR分布。

實(shí)驗(yàn)測(cè)量的方法較其他兩種方法不僅有良好的理論支持，而且更能體現(xiàn)不同類型的天線實(shí)際的熱消融能力，因此具有更高的可靠性。但該方法在實(shí)驗(yàn)和擬合精準(zhǔn)性方面的研究較少。本文根據(jù)2450 MHz微波消融治療儀的大量實(shí)驗(yàn)，定量計(jì)算實(shí)驗(yàn)數(shù)據(jù)處理中溫度誤差和最佳擬合時(shí)間段的選取對(duì)仿真結(jié)果的影響，并提出一種擬合SAR分布的新方法。

1 原理與方法

1.1 實(shí)驗(yàn)原理

實(shí)驗(yàn)所用儀器為南京康有公司研制的KY-2000型2450MHz微波消融治療儀，測(cè)溫系統(tǒng)是基于熱電偶測(cè)溫原理的3497A型數(shù)據(jù)采集儀。實(shí)驗(yàn)體模采用四川大學(xué)江漢保等研制的微波體模配方[13]。為了精確分析體模的溫升情況，采用22根測(cè)溫針，按模具上所規(guī)劃的位置（布針示意圖見(jiàn)圖1）固定插入體模中。

體模中SAR根據(jù)Pennes生物傳熱方程計(jì)算[14]。它是當(dāng)前生物熱領(lǐng)域應(yīng)用最廣泛的模型，也是研究微波輻射與生物組織之間熱相互作用規(guī)律非常經(jīng)典的方程，即：

其中ρ、ρb為組織和血液的密度，單位為kg/m3；c、cb為組織和血液的比熱，單位為J/(kg ·℃)；T、Tb為組織和血液的溫度，單位為℃；λ為組織的導(dǎo)熱率，單位為J/（s ·m·℃）；ωb為血液灌注率，單位為kg/(m3· s)；qm為代謝率，單位為J/(m3· s)；qr為單位組織吸收的微波輻射能，單位為J/(m3·s)；

體模中，ωb、qm均為零，且在微波開(kāi)啟瞬間導(dǎo)熱項(xiàng)可忽略不計(jì)，則此方程可簡(jiǎn)化為位組織吸收的微波輻射SAR等于局部穩(wěn)定變化函數(shù)。因此，只需得到實(shí)驗(yàn)溫升初始階段斜率即可得到該點(diǎn)的SAR值。

1.2 研究方法

用實(shí)驗(yàn)測(cè)量方法計(jì)算SAR分布，具體求解流程見(jiàn)圖2。

從該流程圖可以看出，在實(shí)驗(yàn)及數(shù)據(jù)處理過(guò)程中存在一系列的問(wèn)題會(huì)影響最終溫度場(chǎng)仿真結(jié)果?？偨Y(jié)如下：

（1）實(shí)驗(yàn)中測(cè)溫針直徑為1.8 mm，且兩根測(cè)溫針之間的最短距離僅為5 mm（圖1），由于實(shí)驗(yàn)條件有限，測(cè)溫針的數(shù)量較多，分布密集，在插針過(guò)程中容易產(chǎn)生偏差[15]。并且實(shí)驗(yàn)環(huán)境可變性因素較多，易對(duì)初始條件產(chǎn)生影響，所以定量分析實(shí)驗(yàn)中數(shù)據(jù)測(cè)量的溫度誤差對(duì)結(jié)果的影響是很有必要的。為了研究溫度誤差的影響，對(duì)不同微波加熱功率進(jìn)行≥12組的實(shí)驗(yàn)，取其平均值。以該平均值所得擬合結(jié)果作為標(biāo)準(zhǔn)，計(jì)算在平均值中加入不同溫度誤差后對(duì)溫度場(chǎng)仿真的影響。

為了模擬實(shí)驗(yàn)中的誤差情況，對(duì)12原始組數(shù)據(jù)進(jìn)行分析，發(fā)現(xiàn)實(shí)驗(yàn)中數(shù)據(jù)點(diǎn)溫度誤差的方差不同，范圍為0.07～2.00 ℃，平均值為0.445 ℃。因此在平均值數(shù)據(jù)中加入方差為0.01～5.00 ℃的零均值高斯噪音，對(duì)所得結(jié)果進(jìn)行匯總，觀察不同的溫度誤差對(duì)標(biāo)準(zhǔn)溫度場(chǎng)仿真的影響。

（2）在對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行預(yù)處理時(shí)，主要目的是選擇實(shí)驗(yàn)中有效的數(shù)據(jù)點(diǎn)。由于實(shí)驗(yàn)參數(shù)設(shè)置為10 min，所得溫升變化曲線有較長(zhǎng)一段的線性較好，因此在求取斜率時(shí)我們普遍選取前100 s時(shí)間段的數(shù)據(jù)[16]。但該取值范圍并沒(méi)有較好的理論支持，且不同測(cè)溫點(diǎn)其溫升變化范圍不同，一致選取前100 s時(shí)間段的數(shù)據(jù)顯然不夠科學(xué)。所以需要定量分析不同擬合時(shí)間段的選取對(duì)SAR分布的影響。

具體方法為，以10 s為步長(zhǎng)，分別選取原始數(shù)據(jù)前40 s、50 s…120 s時(shí)間段求取斜率，從而分析擬合結(jié)果。在求取斜率過(guò)程中發(fā)現(xiàn)，當(dāng)時(shí)間段取到前120 s及以上時(shí)，很多測(cè)溫點(diǎn)的斜率出現(xiàn)負(fù)值，不滿足實(shí)驗(yàn)處理要求。因此只需對(duì)前40 s、50 s…110 s時(shí)間段的溫度數(shù)據(jù)進(jìn)行分析即可。

（3）擬合預(yù)處理數(shù)據(jù)。由于各測(cè)溫點(diǎn)的SAR值即為其初始階段的溫升斜率，所以根據(jù)SAR在3個(gè)坐標(biāo)軸的分布函數(shù)，即可求得SAR的三維空間分布[16]。其中r方向?yàn)橹笖?shù)衰減，z方向?yàn)?次多項(xiàng)式衰減。對(duì)SAR分布進(jìn)行擬合的傳統(tǒng)方法為最小二乘法[4]，但由于不同測(cè)溫點(diǎn)斜率值數(shù)量級(jí)相差大，擬合過(guò)程中較小數(shù)量級(jí)的SAR值常被忽略，導(dǎo)致取值動(dòng)態(tài)范圍大，因此得到的擬合結(jié)果非常不穩(wěn)定。為了更好地解決該問(wèn)題，我們采用加權(quán)最小二乘法。

在該方法中，通過(guò)歸一法為所有SAR值引入權(quán)值，使得他們?cè)跀M合中的作用相同；并調(diào)用Matlab R2012a中的Curve Fitting Tool工具箱，通過(guò)設(shè)置最大迭代次數(shù)、誤差閾值和有限差分梯度等擬合參數(shù)實(shí)現(xiàn)擬合結(jié)果的最優(yōu)化。

2 實(shí)驗(yàn)結(jié)果

2.1 實(shí)驗(yàn)中溫度誤差分析

擬合結(jié)果的偏差隨所加溫度誤差的增大而增大，且當(dāng)誤差的方差在3.50 ℃以內(nèi)時(shí)，對(duì)實(shí)驗(yàn)結(jié)果的影響不大。以54 ℃為邊界的等SAR值包絡(luò)線，見(jiàn)圖3。等SAR包絡(luò)線是指對(duì)擬合后的SAR分布進(jìn)行計(jì)算，提取以54 ℃為求解條件的SAR邊界。其中0點(diǎn)即為微波發(fā)射點(diǎn)，x軸為實(shí)驗(yàn)中z方向長(zhǎng)度， y軸為實(shí)驗(yàn)中r方向長(zhǎng)度， 單位為 cm。可以看出，當(dāng)加入溫度誤差的方差為0.01～3.50 ℃時(shí)，與實(shí)驗(yàn)平均值結(jié)果相差不大；當(dāng)加入溫度誤差的方差為4.00～5.00℃時(shí)，較實(shí)驗(yàn)平均值結(jié)果相比出現(xiàn)了一定的偏差，且影響非常明顯。

把實(shí)驗(yàn)結(jié)果與加入不同溫度誤差的SAR分布代入到Ansys微波熱療溫度場(chǎng)仿真模型中，模擬結(jié)果見(jiàn)圖4。圖4中從左到右依次為正常、方差為1 ℃和5 ℃的溫度誤差的消融區(qū)范圍。對(duì)消融區(qū)進(jìn)行測(cè)量，得到不同溫度誤差作用下各消融區(qū)長(zhǎng)軸、短軸的長(zhǎng)度匯總見(jiàn)表1?？梢钥闯?，在溫度場(chǎng)仿真過(guò)程中SAR分布對(duì)消融結(jié)果有很大的影響。當(dāng)所加溫度誤差的方差在3.5 ℃以內(nèi)時(shí)，消融區(qū)各長(zhǎng)軸、短軸長(zhǎng)度較為統(tǒng)一，因此該范圍誤差對(duì)仿真結(jié)果影響不大；而當(dāng)所加溫度誤差的方差為5 ℃時(shí)，對(duì)仿真結(jié)果有較大影響。實(shí)驗(yàn)中最大的溫度誤差為2.00 ℃，所以該溫度誤差對(duì)擬合結(jié)果不會(huì)有太大的影響。

2.2 最佳擬合時(shí)間段的選取

實(shí)驗(yàn)數(shù)據(jù)處理中，擬合時(shí)間的選取對(duì)仿真結(jié)果有影響，且當(dāng)選取前60 s、70 s…110 s時(shí)間段數(shù)據(jù)進(jìn)行處理時(shí)，溫度場(chǎng)仿真結(jié)果較理想。以54 ℃為邊界的等SAR值包絡(luò)線，見(jiàn)圖5?？梢钥闯?，當(dāng)時(shí)間選取前60 s、70 s…110 s時(shí)間段時(shí)，包絡(luò)線內(nèi)面積相差不大，而當(dāng)選取前40 s或50 s時(shí)間段時(shí)，其結(jié)果較正常結(jié)果相比明顯變小。

把選取不同時(shí)間所得SAR分布，利用Ansys進(jìn)行微波熱療溫度場(chǎng)仿真，得到消融區(qū)范圍，見(jiàn)圖6。圖中左右兩個(gè)消融區(qū)分別為取前40 s和100 s時(shí)間段的仿真結(jié)果。對(duì)消融區(qū)進(jìn)行測(cè)量，長(zhǎng)軸、短軸的長(zhǎng)度匯總見(jiàn)表2?？梢钥闯觯?dāng)選取前40 s時(shí)間段進(jìn)行溫度場(chǎng)仿真時(shí)，其結(jié)果較前100 s時(shí)間段的擬合結(jié)果有明顯不同，消融區(qū)的長(zhǎng)軸和短軸明顯變小。因此，當(dāng)微波熱療加熱時(shí)間為10 min時(shí)，取前100 s時(shí)間段進(jìn)行分析是可行的，不會(huì)對(duì)溫度場(chǎng)仿真造成誤差影響。

2.3 實(shí)驗(yàn)數(shù)據(jù)擬合

用加權(quán)最小二乘法對(duì)SAR分布進(jìn)行擬合時(shí)，結(jié)果穩(wěn)定性有明顯提高。用加權(quán)最小二乘法做5次SAR分布擬合，SAR分布中6個(gè)參數(shù)擬合結(jié)果的方差，見(jiàn)表3。可以看到，擬合結(jié)果的波動(dòng)性普遍＜0.01，且最大波動(dòng)＜0.40。說(shuō)明用加權(quán)最小二乘法對(duì)SAR分布進(jìn)行擬合時(shí)結(jié)果波動(dòng)性小，穩(wěn)定性高。

對(duì)同一組實(shí)驗(yàn)測(cè)量數(shù)據(jù)分別使用最小二乘法和加權(quán)最小二乘法進(jìn)行4次擬合，并且求得擬合結(jié)果的Ansys溫度場(chǎng)仿真，見(jiàn)圖7。上面4個(gè)消融區(qū)所使用的擬合方法為最小二乘法，消融區(qū)范圍的變化非常明顯；下面4個(gè)消融區(qū)所使用的擬合方法為加權(quán)最小二乘法，消融區(qū)面積非常穩(wěn)定，說(shuō)明加權(quán)最小二乘法結(jié)果可靠有效，適應(yīng)于臨床研究對(duì)仿真結(jié)果穩(wěn)定性的要求。

3 結(jié)論

通過(guò)本研究得出如下3個(gè)結(jié)論：

（1）實(shí)驗(yàn)中測(cè)量溫度誤差的控制決定了微波熱療溫度場(chǎng)模擬的精度,在SAR的實(shí)驗(yàn)測(cè)量中應(yīng)盡可能精確定位測(cè)溫針的位置，減少溫度測(cè)量誤差。

（2）最佳擬合時(shí)間段的選取對(duì)溫度場(chǎng)仿真中SAR分布的求解有著極其重要的意義，實(shí)驗(yàn)結(jié)果表明最佳SAR擬合時(shí)間段為實(shí)驗(yàn)前60 s、70 s…110 s。

（3）SAR實(shí)驗(yàn)中，用加權(quán)最小二乘法擬合SAR分布較傳統(tǒng)方法更加穩(wěn)定可靠。

[1] Alexakis N,Halloran C,Raraty M,et al.Current standards of surgery for pancreatic cancer[J].Br J Surg,2004,91(11):1410-1427.

[2] Jason Chiang,Kieran Hynes,Christopher L Brace,et al.Flowdependent vascular heat transfer during microwave thermal ablation[J].Annual international conference of the IEEE EMBS,2012,34:5582-5585.

[3] LIANG Ping,DONG Bao-wei,YU Xiao-ling,et al.Computeraided dynamic simulation of microwave-induced thermal distribution in coagulation of liver cancer [J].IEEE transaction on biomedical engineering,2001,48(7):821-829.

[4] ZHU Liang,XU Lisa X,Norbert Chencinski,et al.Quantification of the 3-D electromagnetic power absorption rate in tissue during transurethral prostatic microwave thermotherapy using heat transfer model[J]. IEEE transactions on biomedical engineering,1998,45(9):1163-1172.

[5] Wessapan T,Srisawatdhisukul S,Rattanadecho P,et al.The effect of dielectric shield on specific absorption rate and heat transfer in the human body exposed on leakage microwave energy[J].Int Commum Heat Mass Transfer,2011,38(2):255-262.

[6] Grant P Steven,LI Qing,ZHANG Zhong-pu.The constract of two kinds of microwave antenna SAR simulation[J].Mechanics and Materials,2014,553:379-383.

[7] Salloum M,Ma R,Zhu L.An in-vivo experimental study of temperature elevations in animal tissue during magnetic nanoparticle hyperthermia[J].Int J Hyperthermia,2008,24(7):589-601.

[8] Tunc M,Parmaksizoglu C,Cikrikci S.The bio-heat transfer equation and its applications in hyperthermia treatments[J]. Engineering Computations,2006,23(4):451-463.

[9] YANG Lei,HAO Dong-mei,WU Shui-cai,et al.SAR and temperature distribution in the rat head model exposed to electromagnetic field radiation by 900 MHz dipole antenna[J]. Australas Phy Eng Sci Med,2013,36(2):251-257.

[10] Gentili GB,Gori F,Leoncini M.Electromagnetic and thermal models of a water-cooled dipole radiating in a biological tissue[J]. IEEE Trans Biomed.Eng,1991,38(1):98-103.

[11] 翟偉明,盛林,宋亦旭,等.基于影像引導(dǎo)的計(jì)算機(jī)輔助肝癌微波消融[J].計(jì)算機(jī)研究與發(fā)展,2011,48(2):281-288.

[12] 趙磊,繩秀君,吳水才,等.腫瘤熱療中組織熱物性參數(shù)對(duì)熱場(chǎng)分布的影響[J].中國(guó)醫(yī)療設(shè)備,2009,24(4):15-22.

[13] 江漢保,郝晉,李愛(ài)華,等.微波體模[J].中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào),1992,11:199-203.

[14 ]Pennes HH.Analysis of tissue and arterial blood temperatures in the resting human forearm[J].Journal of Applied Physiology,1998,85(1):5-34.

[15] 翟偉明.影像引導(dǎo)下計(jì)算機(jī)輔助介入手術(shù)導(dǎo)航關(guān)鍵技術(shù)的研究[D].北京:清華大學(xué),2010.

[16] 梁萍,董寶瑋,于曉玲,等.超聲引導(dǎo)下植入式微波熱場(chǎng)在體模中的研究[J].中華超聲影像學(xué)雜志,2002,11(4):236-239.

Research on Inf l uence Factors of SAR Distribution in Temperature Field Simulation of Microwave Hyperthermia

ZHANG Man, WU Shui-cai, HU Rui, LIN Lan, GAO Hong-jian
College of Life Science and Biological Engineering, Beijing University of Technology, Beijing 100124, China

In order to improve the accuracy of temperature fi eld simulation of microwave hyperthermia, this study analyzed the influence factors of SAR distribution in terms of quantitative calculation of temperature error, determination of effective time points of SAR fi tting and introduction of optimum fi tting method of SAR distribution. The results of massive calculation and simulation comparison showed that the weighted least squares fi tting method would be more stable and effective than the traditional method when the temperature error was controlled within 3.5 ℃ and the time points of SAR fi tting were choosed as 60 s、70 s…110 s.

microwave hyperthermia; temperature fi eld simulation; SAR fi tting; weighted least square method

R318.04；R318.6

A

10.3969/j.issn.1674-1633.2014.11.007

1674-1633(2014)11-0025-04

2014-09-23

國(guó)家自然科學(xué)基金專項(xiàng)基金項(xiàng)目資助（81127006）。

作者郵箱：284134211@qq.com