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

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

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

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

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

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

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

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

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

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

1 原理與方法

1.1 實驗原理

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

體模中SAR根據(jù)Pennes生物傳熱方程計算[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均為零，且在微波開啟瞬間導(dǎo)熱項可忽略不計，則此方程可簡化為位組織吸收的微波輻射SAR等于局部穩(wěn)定變化函數(shù)。因此，只需得到實驗溫升初始階段斜率即可得到該點的SAR值。

1.2 研究方法

用實驗測量方法計算SAR分布，具體求解流程見圖2。

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

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

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

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

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

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

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

2 實驗結(jié)果

2.1 實驗中溫度誤差分析

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

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

2.2 最佳擬合時間段的選取

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

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

2.3 實驗數(shù)據(jù)擬合

用加權(quán)最小二乘法對SAR分布進(jìn)行擬合時，結(jié)果穩(wěn)定性有明顯提高。用加權(quán)最小二乘法做5次SAR分布擬合，SAR分布中6個參數(shù)擬合結(jié)果的方差，見表3?？梢钥吹?，擬合結(jié)果的波動性普遍＜0.01，且最大波動＜0.40。說明用加權(quán)最小二乘法對SAR分布進(jìn)行擬合時結(jié)果波動性小，穩(wěn)定性高。

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

3 結(jié)論

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

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

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

（3）SAR實驗中，用加權(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].計算機研究與發(fā)展,2011,48(2):281-288.

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

[13] 江漢保,郝晉,李愛華,等.微波體模[J].中國生物醫(yī)學(xué)工程學(xué)報,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)下計算機輔助介入手術(shù)導(dǎo)航關(guān)鍵技術(shù)的研究[D].北京:清華大學(xué),2010.

[16] 梁萍,董寶瑋,于曉玲,等.超聲引導(dǎo)下植入式微波熱場在體模中的研究[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

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

作者郵箱：284134211@qq.com