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

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

北京工業大學 生命科學與生物工程學院，北京 100124

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

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

北京工業大學 生命科學與生物工程學院，北京 100124

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

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

引言

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

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

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

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

1 熱消融溫度場的建模

熱消融仿真模型的構建主要基于電磁技術和生物傳熱技術，通過電磁技術獲得各個點的熱源，并且通過生物傳熱技術獲得各個點的溫度分布。

1.1 電磁源和生物傳熱建模

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

其中J表示電功，σ表示電導率，E表示電場強度，V表示電壓源，由功率公式可知：

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

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

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

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

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

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

圖2 多極射頻消融電極的實際模型

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

2 特性參數的敏感性分析

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

表1 特性參數和取值范圍

表2 方差分析的32實驗安排

2.1 遠場點參數敏感性分析

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

圖3 各參數在t=50 s（a）和t=360 s（b）時對遠場點溫度的主效應圖

圖4 仿真過程中各參數對遠場點方差貢獻率趨勢圖

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

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

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

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

2.2 近場點參數敏感性分析

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

表3 t=50 s以及t=360 s 時的遠場點方差分析

圖6 各參數在t=50 s（a）和t=360 s（b）時對近場點溫度的主效應圖

圖7 仿真過程中各參數對近場點方差貢獻率趨勢圖

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

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

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

表4 t=50 s以及t=360 s時的近場點方差分析

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

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

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

3 結論

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

圖9 仿真過程中各參數對40℃等溫面包容體積方差貢獻率趨勢圖

圖10 仿真過程中各參數對80℃等溫面包容體積方差貢獻率趨勢圖

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

[1] 周茂勝,常欣.射頻消融原理及應用[J].當代醫學,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] 張曼.腫瘤熱消融治療手術規劃系統及關鍵技術研究[D].北京:北京工業大學,2016.

[6] 聶曉慧.脊柱腫瘤射頻消融適形治療的溫度場研究[D].北京:北京工業大學,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] 黃維,羅洪艷,潘進洪,等.構建基于肝臟CT圖像的溫控射頻消融仿真模型[J].中國介入影像與治療學,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].北京師范大學學報(自然科學版),2008,44(1):1.

[11] 李進霞,張偉杰,張素香.有限元法及應用狀況[J].科技創新導報,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] 包家立.電磁醫療設備的生物物理基礎與應用[J].中國醫療設備,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] 孫登華,錢鋒,孫光,等.射頻消融技術的臨床應用進展[J].吉林醫學,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

國家自然科學基金資助項目（71661167001）；北京市自然科學基金資助項目（7174279）。

吳水才，教授，主要研究方向為生物醫學信息檢測與處理、生物醫學電子與醫療儀器，生物醫學圖像處理。

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