溫控射頻消融溫度場仿真模型參數(shù)敏感性分析

王笑茹，高宏建，吳水才，白燕萍

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

溫控射頻消融溫度場仿真模型參數(shù)敏感性分析

王笑茹，高宏建，吳水才，白燕萍

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

目的研究溫控射頻消融溫度場仿真模型中特性參數(shù)的敏感性，分析特性參數(shù)對于消融區(qū)溫度分布的影響。方法利用有限元技術(shù)建立溫度場仿真模型，基于方差分析理論研究基本特性參數(shù)變化對于消融區(qū)尺寸和不同點(diǎn)溫度的方差貢獻(xiàn)率（SS%）和主效應(yīng)。結(jié)果研究表明，在近場區(qū)域中，消融前期電導(dǎo)率（Sigma）和組織等效電阻（R）的SS%約為22%，消融后期SS%約為47.5%；在遠(yuǎn)場區(qū)域中，消融前期Sigma和R的SS%約為26%，消融后期SS%約為43%；對于消融區(qū)尺寸而言，消融前期Sigma和R的SS%約為36%，消融后期SS%約為42%；根據(jù)主效應(yīng)分析，參數(shù)Sigma和R與溫度成正相關(guān)，比熱容和密度與溫度成負(fù)相關(guān)，導(dǎo)熱率與溫度的相關(guān)性隨時(shí)間和位置而有所變化。結(jié)論在射頻熱消融過程中，相比于導(dǎo)熱率、密度、比熱容和介電常數(shù)，電導(dǎo)率和組織等效電阻具有顯著敏感性。

射頻消融；溫度場仿真；模型參數(shù)；敏感性分析

引言

射頻熱消融（Radiofrequency Ablation，RFA）技術(shù)是將交流電場施加到腫瘤組織，通過電阻加熱提升靶組織的溫度，使之產(chǎn)生凝固性壞死，從而實(shí)現(xiàn)對腫瘤的原位滅活[1]。該技術(shù)因具有微創(chuàng)、治療效果顯著等優(yōu)點(diǎn)而在肝腫瘤治療中得到了廣泛應(yīng)用，但熱消融的質(zhì)量仍取決于臨床醫(yī)生的經(jīng)驗(yàn)。研究表明[2-4]，熱消融效果主要取決于消融區(qū)的范圍、形狀及其內(nèi)部的溫度場分布。因此，熱消融溫度場的精確表征顯得尤為重要。

計(jì)算機(jī)建模仿真技術(shù)不僅能夠觀測消融區(qū)溫度變化，還有助于醫(yī)生術(shù)前制定合理手術(shù)計(jì)劃。張曼[5]基于Pennes和Hyperbolic傳熱方程，討論了不同傳熱模型對于溫度場的影響。聶曉慧[6]研究了脊柱腫瘤射頻消融適形治療的溫度場。Jin等[7]建立了基于MRI的甲狀腺射頻消融有限元模型。但這些文獻(xiàn)中將溫度場模型中的熱參數(shù)和電參數(shù)取為固定值，未考慮生物組織熱物性參數(shù)的特異性變化，從而使得消融結(jié)果存在較大的誤差。黃維等[8]基于反饋溫度的比例積分控制方法調(diào)節(jié)電極針施加電壓，構(gòu)建溫控射頻消融仿真模型，但未涉及特性參數(shù)對于中心電壓的影響。針對上述問題，可通過溫度反饋來調(diào)節(jié)組織參數(shù)，減少仿真誤差。由于溫度場仿真模型中具有多個(gè)特性參數(shù)，因此反饋調(diào)節(jié)需要重點(diǎn)關(guān)注對誤差影響顯著的因素，即敏感性參數(shù)。為了解決這一問題，本文通過方差分析對溫度場模型中的特性參數(shù)進(jìn)行敏感性分析，獲得對模型具有顯著性影響的參數(shù)，由此為特性參數(shù)的精確表征和溫度場的反饋調(diào)節(jié)提供科學(xué)依據(jù)。

敏感性分析是一種定量描述輸入?yún)?shù)對輸出響應(yīng)程度的方法[9]，其量化指標(biāo)為誤差的方差貢獻(xiàn)率。方差貢獻(xiàn)率越大，則參數(shù)對模型響應(yīng)的影響也就越大[10]。本研究的目的在于通過敏感性分析，獲得各特性參數(shù)的影響顯著性，在反饋應(yīng)用中去除敏感性較低的參數(shù)，重點(diǎn)考慮敏感性較高的參數(shù)，由此有效地降低模型的復(fù)雜度，減少數(shù)據(jù)運(yùn)算，進(jìn)一步提高仿真模型的精度。

在臨床射頻熱消融手術(shù)中，常用的是溫控射頻消融儀，其根據(jù)所設(shè)置的參數(shù)（中心溫度、溫升速率等）實(shí)時(shí)地調(diào)節(jié)輸出功率，通過功率補(bǔ)償獲得恒定的中心溫度，最終獲得所需的熱凝固效果。本研究基于恒溫?zé)嵯趦xRFA-I（Blade Co.，Ltd.，Beijing，China）和仿肝組織體模建立溫控射頻熱消融溫度場模型，建模軟件為COMSOL Multiphysics軟件（COMSOL Inc.，Palo Alto，CA，USA）；然后利用Minitab（State College，PA，USA）中的因子分析來研究各特性參數(shù)的敏感性，為后續(xù)反饋工作提供可靠依據(jù)。

1 熱消融溫度場的建模

熱消融仿真模型的構(gòu)建主要基于電磁技術(shù)和生物傳熱技術(shù)，通過電磁技術(shù)獲得各個(gè)點(diǎn)的熱源，并且通過生物傳熱技術(shù)獲得各個(gè)點(diǎn)的溫度分布。

1.1 電磁源和生物傳熱建模

在射頻消融仿真過程中，電流場作為熱源，是焦耳熱產(chǎn)生的主要原因，其控制方程為拉普拉斯方程(1)～(3)：

其中J表示電功，σ表示電導(dǎo)率，E表示電場強(qiáng)度，V表示電壓源，由功率公式可知：

其中R（單位：Ω）為組織的等效電阻值，P（單位：W）根據(jù)實(shí)驗(yàn)中的溫控射頻消融儀獲得，在本研究中為時(shí)間的函數(shù)：

本研究中采用的生物傳熱方程為經(jīng)典的Pennes生物傳熱方程（PBE），具體描述如下：

其中T為溫度（℃），t為時(shí)間（s），ρ為體積質(zhì)量密度（kg/m3），c為比熱容（J/kg·℃），k為導(dǎo)熱率[J/(s·m·℃)]，w為體積血液灌注率[kg/(m3·s)]，Q為代謝熱生成速率[J/(m3·s)]，S為比吸收率SAR（W/kg），下標(biāo)b表示組織的血液特性。

1.2 基于有限元技術(shù)的溫度場模型求解

模型的建立和求解在COMSOL Multiphysics軟件中完成，其中采用的技術(shù)為有限元方法（Finite Element Methods，F(xiàn)EM），其通過計(jì)算機(jī)采用分片近似，逼近整體的研究思想來求解物理問題[11]。本文采用基于多極射頻電極的溫控射頻消融溫度場模型，該模型的尺寸與實(shí)驗(yàn)體模（50 mm×50 mm×70 mm）的尺寸相同?？紤]到射頻電極分布的對稱性，并且為了減少仿真計(jì)算時(shí)間，僅建立1/2幾何模型，大小為50 mm×25 mm×70 mm，見圖1。多極射頻電極形狀，見圖2。電極周圍的長方體為仿肝組織，6個(gè)子電極均勻分布在幾何模型的一側(cè)，兩個(gè)子電極間的夾角為55°，每個(gè)子電極的撓曲度不同，相鄰兩個(gè)子電極的撓曲度為104°和57°。

圖1 多極射頻消融溫度場仿真的幾何模型

圖2 多極射頻消融電極的實(shí)際模型

為了得到有限元模型的定解，設(shè)定仿真模型的邊界條件和初始條件。對于電流場，多極射頻電極的邊界設(shè)定為熱源，肝組織邊界設(shè)定為地電極，電勢為0 V；對于生物傳熱物理場，模型初始溫度為28℃，模型外部熱邊界條件為28℃；模型加熱時(shí)間為360 s；忽略組織血液灌注率w對仿真溫度的影響。

2 特性參數(shù)的敏感性分析

為研究熱物性參數(shù)和電特性參數(shù)對熱消融溫度場的影響，利用Minitab軟件設(shè)計(jì)32組仿真實(shí)驗(yàn)，具體設(shè)置見表1。其中各參數(shù)的取值范圍源自已有報(bào)道[12-13]，見表2。為了研究特性參數(shù)對溫度分布的影響，選取包括遠(yuǎn)場點(diǎn)（1.5 cm，0 cm）和近場點(diǎn)（0.5 cm，0 cm）的多個(gè)點(diǎn)以及消融區(qū)等溫面進(jìn)行敏感性分析。

表1 特性參數(shù)和取值范圍

表2 方差分析的32實(shí)驗(yàn)安排

2.1 遠(yuǎn)場點(diǎn)參數(shù)敏感性分析

根據(jù)表2中的32組實(shí)驗(yàn)安排，分別在Comsol軟件中進(jìn)行溫度場模型仿真，每組實(shí)驗(yàn)消融時(shí)間為360 s。參數(shù)R、Cp、k、Sigma、rho和Epsilon在t=50 s和t=360 s時(shí)對遠(yuǎn)場點(diǎn)的主效應(yīng)圖，見圖3。該主效應(yīng)圖反映了各參數(shù)對消融溫度的影響趨勢。各參數(shù)在整個(gè)仿真過程中遠(yuǎn)場點(diǎn)方差貢獻(xiàn)率的趨勢圖，見圖4，其中Cp和rho變化趨勢相同，在圖中重合；R和Sigma的變化趨勢十分相似。由圖3和圖4可知，隨著仿真時(shí)間的增加，Cp和rho對遠(yuǎn)場點(diǎn)消融溫度的影響越來越小，且呈負(fù)相關(guān)；k對遠(yuǎn)場點(diǎn)消融溫度的影響一直很弱；R和Sigma對遠(yuǎn)場點(diǎn)消融溫度的影響一直很大，為顯著性影響因素；Epsilon對消融溫度結(jié)果無影響。

圖3 各參數(shù)在t=50 s（a）和t=360 s（b）時(shí)對遠(yuǎn)場點(diǎn)溫度的主效應(yīng)圖

圖4 仿真過程中各參數(shù)對遠(yuǎn)場點(diǎn)方差貢獻(xiàn)率趨勢圖

為了得到各參數(shù)對消融溫度場影響的定量分析，進(jìn)一步對各組實(shí)驗(yàn)結(jié)果進(jìn)行方差分析，結(jié)果見表3。F值是衡量數(shù)據(jù)均值的方差，該值越大則參數(shù)影響越顯著；P值為用于判定假設(shè)檢驗(yàn)結(jié)果的參數(shù)[9]，根據(jù)方差分析，P值＜0.05說明該參數(shù)具有顯著性；Adj_SS(SS)表示參數(shù)的方差和；Adj_MS(MS)表示參數(shù)方差的平均值；貢獻(xiàn)率SS%表示方差偏離均值的平方和，其大小可定量地反應(yīng)各參數(shù)對消融溫度場的影響顯著性，如公式(7)所示：

式中，i取A、B、C、D、E和F；Adj_SSi表示參數(shù)i的方差和。

表3中的SS%結(jié)果和圖4的趨勢圖表明，消融前期各參數(shù)敏感性為：Cp（rho）＞Sigma＞R＞k＞兩因子交互作用＞Epsilon（Epsilon=0，故表3中未列出）；在消融后期各參數(shù)的敏感性為：Sigma＞R＞Cp（rho）＞k＞兩因子交互作用＞Epsilon；其中Epsilon對消融溫度始終沒有影響，因此在后續(xù)分析中可忽略參數(shù)Epsilon。

為了直觀地觀測對溫度場影響顯著的參數(shù)Sigma和R之間的交互作用，進(jìn)一步繪制出不同參數(shù)和消融溫度結(jié)果的等值線圖，該等值線圖反映了兩個(gè)參數(shù)變量和響應(yīng)之間的關(guān)系，見圖5。該圖表示Sigma和R同時(shí)增大時(shí)，消融溫度會略有增加，如圖5中右上角區(qū)域所示。另外結(jié)合因子交互作用的SS%值（＜0.2%），可忽略不計(jì)參數(shù)間交互作用對消融區(qū)的影響。

2.2 近場點(diǎn)參數(shù)敏感性分析

圖5 Sigma和R在t=50 s（a）和t=360 s（b）的遠(yuǎn)場等溫線圖

表3 t=50 s以及t=360 s 時(shí)的遠(yuǎn)場點(diǎn)方差分析

圖6 各參數(shù)在t=50 s（a）和t=360 s（b）時(shí)對近場點(diǎn)溫度的主效應(yīng)圖

圖7 仿真過程中各參數(shù)對近場點(diǎn)方差貢獻(xiàn)率趨勢圖

該分析采用與遠(yuǎn)場點(diǎn)參數(shù)敏感性分析相同的模型仿真條件。參數(shù)R、Cp、k、Sigma、rho和Epsilon在t=50 s、和t=360 s時(shí)對近場點(diǎn)的主效應(yīng)圖，見圖6。各參數(shù)在整個(gè)仿真過程中近場點(diǎn)方差貢獻(xiàn)率的趨勢圖，見圖7，其中Cp和rho變化趨勢相同，在圖中重合；R和Sigma的變化趨勢十分相似。由圖6和圖7可知，隨著仿真時(shí)間的增加，Cp和rho對近場點(diǎn)消融溫度的影響越來越小，呈負(fù)相關(guān)；k對近場點(diǎn)消融溫度的影響一直很弱；R和Sigma對遠(yuǎn)場點(diǎn)消融溫度的影響一直很大，為顯著性影響因素；Epsilon對近場點(diǎn)溫度始終無影響。

t=50 s以及t=360 s時(shí)的近場點(diǎn)方差分析結(jié)果，見表4。表4中的SS%結(jié)果和圖7的趨勢圖表明，消融前期各參數(shù)的敏感性大小為：Cp（rho）＞Sigma＞R＞k＞交互作用＞Epsilon；消融后期各參數(shù)的敏感性為：Sigma＞R＞Cp（rho）＞k＞交互作用＞Epsilon。

參數(shù)Sigma和R的交互作用對消融溫度的影響，見圖8。隨著Sigma和R的增大近場點(diǎn)消融溫度升高，如圖中右上角區(qū)域所示。結(jié)合因子交互作用的SS%值（＜0.2%），故可忽略不計(jì)參數(shù)間交互作用對消融區(qū)的影響。

表4 t=50 s以及t=360 s時(shí)的近場點(diǎn)方差分析

圖8 Sigma和R在t=50 s（a）和t=360 s（b）的近場等溫線

2.3 各參數(shù)對40℃和80℃等溫面包容體積的敏感性分析

選取不同位置的點(diǎn)研究特性參數(shù)的敏感性，具有特例性，因此本研究進(jìn)一步討論了各參數(shù)對不同等溫面包容體積的敏感性。在熱消融手術(shù)中，通常選用50℃等溫面[14-15]和60℃等溫面[16-17]作為消融區(qū)邊界，即可認(rèn)為50℃等溫面包容體積或60℃等溫面包容體積為熱損傷組織體積。因此本研究中分別選取40℃等溫面和80℃等溫面作為近場點(diǎn)和遠(yuǎn)場點(diǎn)的綜合影響效果。各參數(shù)在不同時(shí)刻對40℃和80℃等溫面包容體積的方差貢獻(xiàn)率趨勢圖，見圖9～10。如圖所示，其中Cp和rho變化趨勢相同，在圖中重合；R和Sigma的變化趨勢十分相近；消融期間參數(shù)R和Sigma對等溫面包容體積影響一直很顯著；Epsilon對等溫面包容體積無影響；參數(shù)Cp和rho對等溫面包容體積影響敏感性隨著仿真時(shí)間的增加先略有所增加然后變小，整體影響趨勢是減小的；這些結(jié)果與單點(diǎn)敏感性分析結(jié)果具有較好的一致性。

3 結(jié)論

本研究基于溫控射頻熱消融仿真模型，分析了包括熱參數(shù)和電參數(shù)在內(nèi)的特性參數(shù)對近場點(diǎn)、遠(yuǎn)場點(diǎn)和不同等溫面包容體積的影響。結(jié)果表明，參數(shù)Sigma和R在消融期內(nèi)對各點(diǎn)消融溫度的敏感性始終很顯著，且成正相關(guān)；參數(shù)Cp和rho對各點(diǎn)消融溫度的敏感性隨著仿真時(shí)間的增加而減小，且成負(fù)相關(guān)；參數(shù)Epsilon在消融期間內(nèi)對各點(diǎn)消融溫度沒有影響。k的變化沒有較好的一致性：對于近場點(diǎn)而言，開始階段通過熱傳導(dǎo)吸收熱量，所以k增加時(shí)該點(diǎn)消融溫度會升高，成正相關(guān)；消融后期通過傳熱傳導(dǎo)釋放熱量，因此與k成負(fù)相關(guān)。對于遠(yuǎn)場點(diǎn)而言，主要通過熱傳導(dǎo)吸收熱量，因此與k成正相關(guān)。k對不同等溫面包容體積的敏感性在消融前期略有不同，在消融中后期變化趨勢基本一致。k的SS%始終很小，可忽略不計(jì)，在今后的研究中可將k取為最優(yōu)固定值。由于Sigma和R在消融期內(nèi)對消融溫度的敏感性很顯著，因此在反饋調(diào)節(jié)中可重點(diǎn)考慮反饋參數(shù)Sigma和R。對敏感性不顯著的參數(shù)如Cp、rho、k和Epsilon可取為最優(yōu)固定值。

圖9 仿真過程中各參數(shù)對40℃等溫面包容體積方差貢獻(xiàn)率趨勢圖

圖10 仿真過程中各參數(shù)對80℃等溫面包容體積方差貢獻(xiàn)率趨勢圖

綜上所述，本研究通過敏感性分析獲得了溫控射頻熱消融仿真模型中的顯著性影響因素Sigma和R，為敏感性參數(shù)的反饋調(diào)節(jié)中提供了科學(xué)依據(jù)。

[1] 周茂勝,常欣.射頻消融原理及應(yīng)用[J].當(dāng)代醫(yī)學(xué),2009,15(21): 18-19.

[2] Clasen S,Rempp H,Schmidt D,et al.Multipolar radiofrequency ablation using internally cooled electrodes in ex vivo bovine liver: correlation between volume of coagulation and amount ofapplied energy[J].Eur J Radiol,2012,81(1):111-113.

[3] Berjano EJ,Hornero F.Thermal-electrical modeling for epicardial atrial radiofrequency ablation[J].IEEE Trans Biomed Eng,2004,51(8):1348-1357.

[4] Matschek J,Bullinger E,Haeseler FV,et al.Mathematical 3D modelling and sensitivity analysis of multipolar radiofrequency ablation in the spine[J].Math Biosci,2017,284:51-60.

[5] 張曼.腫瘤熱消融治療手術(shù)規(guī)劃系統(tǒng)及關(guān)鍵技術(shù)研究[D].北京:北京工業(yè)大學(xué),2016.

[6] 聶曉慧.脊柱腫瘤射頻消融適形治療的溫度場研究[D].北京:北京工業(yè)大學(xué),2016.

[7] Jin C,He Z,Liu J.MRI-based finite element simulationon radiofrequency ablation of thyroid cancer[J].Comput Methods Programs Biomed,2014,113(2):529-538.

[8] 黃維,羅洪艷,潘進(jìn)洪,等.構(gòu)建基于肝臟CT圖像的溫控射頻消融仿真模型[J].中國介入影像與治療學(xué),2014,11(8):532-536.

[9] Bueno C,Hauschild MZ,Rossignolo JA,et al.Sensitivity analysis of the use of life cycle impact assessment methods:a case study on building materials[J].J Clean Prod,2016,112(1):2208-2220.

[10] 蔡毅,邢巖,胡丹.敏感性分析綜述[J].北京師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2008,44(1):1.

[11] 李進(jìn)霞,張偉杰,張素香.有限元法及應(yīng)用狀況[J].科技創(chuàng)新導(dǎo)報(bào),2012,(31):120-121.

[12] Rossmanna C,Haemmerich D.Review of temperature dependence of thermal properties, dielectric properties,and perfusion of biological tissues at hyperthermic and ablation temperatures[J].Crit Rev Biomed Eng,2014,42(6):467-492.

[13] Cavagnaro M,Pinto R,Lopresto V.Numerical models of microwave thermal ablation procedures[A].Microwave Conference (EuMC),2014 44thEuropean[C].New York:IEEE,2014:480-483.

[14] 包家立.電磁醫(yī)療設(shè)備的生物物理基礎(chǔ)與應(yīng)用[J].中國醫(yī)療設(shè)備,2016,31(4):6-13.

[15] Peng T,O’Neill D,Payne S.Mathematical study of the effects of different intrahepatic cooling on thermal ablation zones[A]. International Conference of the IEEE Engineering in Medicine and Biology Society[C].New York:IEEE,2011:6866-6869.

[16] 孫登華,錢鋒,孫光,等.射頻消融技術(shù)的臨床應(yīng)用進(jìn)展[J].吉林醫(yī)學(xué),2012,33(13):2823-2825.

[17] Haase S,Süss P,Schwientek J,et al.Radiofrequency ablation planning: An application of semi-infinite modelling techniques[J].Eur J Oper Res,2012,218(3):856-864.

本文編輯 袁雋玲

Parametric Sensitivity Analysis of Simulation Model for Temperature-Controlled Radiofrequency Ablation Temperature Field

WANG Xiaoru, GAO Hongjian, WU Shuicai, BAI Yanping
College of Life Science & Bio-Engineering, Beijing University of Technology, Beijing 100124, China

ObjectiveTo study the sensitivity of the characteristic parameters in the temperature-controlled radiofrequency ablation (RFA) temperature field model and to analyze the signi ficance of these parameters for the temperature distribution.MethodsThe finite element method was used to establish the temperature field simulation model. The variance contribution rate (SS%) and the main effect of the characteristic parameters on the size and the temperature of the ablation area were then discussed based on the variance analysis theory.ResultsThe study showed that in the near field region, the SS% of the electrical conductivity (Sigma) and the tissue equivalent resistance (R) were about 22% in the early ablation session and about 47.5% in the posterior ablation session. In the far field region, the SS% of Sigma and R were about 26.9% and 43% respectively in the early and posterior ablation stage; with regard to ablation area size, the SS% of Sigma and R were about 36% and 46.5% in the early and posterior ablation stage, respectively. According to the main effect analysis, Sigma and R were directly correlated with temperature, speci fic heat capacity and density were inversely-related to temperature, and the correlation between thermal conductivity and temperature varied with time and position.ConclusionDuring in the RFA thermal ablation, compared with the thermal conductivity, density, speci fic heat capacity and dielectric constant, electrical conductivity and tissue equivalent resistance have signi ficant sensitivity.

radiofrequency ablation; temperature field simulation; model parameters; sensitivity analysis

R318

A

10.3969/j.issn.1674-1633.2017.09.006

1674-1633(2017)09-0023-06

2017-05-16

2017-05-23

國家自然科學(xué)基金資助項(xiàng)目（71661167001）；北京市自然科學(xué)基金資助項(xiàng)目（7174279）。

吳水才，教授，主要研究方向?yàn)樯镝t(yī)學(xué)信息檢測與處理、生物醫(yī)學(xué)電子與醫(yī)療儀器，生物醫(yī)學(xué)圖像處理。

通訊作者郵箱：wushuicai@bjut.edu.cn