孫劼,張璞,李成偉,劉文麗
中國計量科學研究院 醫學與生物計量研究所,北京 100013
上世紀70年代,Lauterbur提出磁共振成像(Magnetic Resonance Imaging,MRI)的方法和實驗裝置[1],經過近50年的發展,醫用MRI設備(以下簡稱“MRI設備”)已比較普及。為提高圖像信噪比,MRI設備正從低場(≤1.5 T)向高場、超高場(≥7 T)方向發展[2-3]。隨著MRI設備磁場強度的增強,其射頻場電磁波與人體組織相互作用也在加強,使受檢者受到熱損傷的風險提高,由此引發的安全問題不容忽視。
射頻(Radio Frequency,RF)是高頻電流通過導體產生的電磁波,頻率范圍在300 kHz~30 GHz之間[4]。RF與人體組織相互作用,攜帶的能量被組織吸收,導致局部溫度升高的現象稱為生物熱效應。特定能量吸收率(Speci fi c Absorption Rate,SAR)是定量研究射頻生物熱效應的指標,其定義為生物組織單位時間、單位質量吸收射頻的能量,單位為W/kg[5]。SAR可分為全身SAR和局部SAR,全身SAR指人體全身組織SAR的平均值,局部SAR指人體任意10 g組織SAR的平均值[5]。
為保證受檢者的安全,相關國際組織[6-7]起草了相應的標準、國際建議及指南,對MRI設備中SAR的安全限值做出了規定。其中,國際電工委員會(International Electrotechnical Commission,IEC)在IEC 60601-2-33中規定MRI設備射頻信號在受檢者體內產生SAR的安全限值(表1)[6]。美國食品藥品監督管理局在其指南中規定“使用MRI設備掃描人體時,受檢者頭部和軀干任意1 g組織SAR的安全限值為8 W/kg”[7]。
使用MRI設備掃描人體過程中,若要評估受檢者體內不同組織SAR的分布情況,需先計算出設備附載人體后內部電磁場的分布狀態。這需要用到電磁場數值計算方法并配合Mimics等專用軟件,在給定初值與邊界條件下求出Maxwell方程組的數值解,即得到人體組織內電磁場分布情況,從而得到相應區域內任一一點的SAR值。

表1 IEC 60601-2-33中規定醫用磁共振成像設備SAR限值(W/kg)
自20世紀50年代有學者以Maxwell方程為研究對象,求出近似解析解開始,隨著數學理論及計算機硬件不斷升級,電磁場數值計算方法也在不斷地提高和發展[8]。目前,在基本算法基礎上,發展出了無單元法、小波分析法、混合算法等一系列數值計算方法[9]。伴隨高性能GPU、云計算的出現,基于各種硬件技術的加速算法也成為電磁場數值計算方法的發展方向。
2.1.1時域有限差分法
時域有限差分法(Finite Difference Method,FDTD)由Yee提出[10]。計算時將Maxwell微分方程按照三維Yee元胞結構(圖1)抽樣離散,并在時域中對磁場和電場交替抽樣。通過以上方法將Maxwell方程離散后轉化為顯式差分方程,并在時域內進行迭代求解[11]。FDTD是目前廣泛應用的一種數值計算方法。

圖1 三維Yee元胞結構圖
2.1.2矩量法
矩量法(Method Of Moments,MOM)由Harrington在專著中提出,其原理是將頻域上Maxwell積分方程組轉化為近似有限求和的形式,并以此為基礎建立代數方程組,通過求出代數方程組的數值進而得到電磁場中每一點磁場、電場分布的數值計算方法[12]。矩量法具有求解過程簡單、步驟統一、應用便捷等優點,同時也存在著計算量大、占用大量計算機資源的缺點。目前矩量法主要應用于天線、電磁波散射等方面問題的計算。
2.1.3有限元法
有限元法(Finite Element Method,FEM)于20世紀60年代提出,Silvester最早將該方法應用于電磁場問題計算。有限元法是依據變分原理將復雜問題離散化、分解為有限多個簡單問題后再分別進行求解的一種近似數值計算方法[13]。有限元法具有求解容易、收斂性好、計算程序通用性強等優點,但對無邊界區域卻難以求解。
2.1.4無單元法
無單元法(Element-Free Method,EFM)可突破單元域瓶頸求解復雜邊界條件下的邊值問題,是一種解決工程電磁場計算問題的有效新方法。無單元法由Lancaster提出,是基于滑動最小二乘法將計算場域離散成若干獨立的點構造場函數,并利用初始邊界條件求出數值解的方法。國外有學者使用該方法解決二維力學場邊值問題[14-15]。劉素貞等[16]又將無單元法應用于電磁場數值計算研究,王立鵬等[17]嘗試了基于伽遼金型無單元多重網格法,進一步解決了無單元法計算電磁場時效率不高的問題,并取得一定成效。
2.1.5小波分析法
法國工程師Morlet率先提出了小波分析的概念[18]。與傳統的分析方法相比,小波分析采用的局部變換方法可對函數或信號進行多尺度細化分析,從而處理Fourier變換難以解決的問題[19]。小波分析法在電磁場數值計算中的應用主要有:小波變換和小波展開。Baharav等[20]通過對經小波變換后的阻抗矩陣方程進行壓縮提高方程的稀疏性,降低運算量并提高電磁場數值計算效率。Wang[21]將小波分析中的小波展開與邊界元算法相結合,提出基于多尺度小波分解的電磁場微分方程解法,用于分析多導體傳輸中的電磁場問題。楊仕友等[22]將一般緊支集正交小波應用于電磁場數值計算,提出一種適合處理通用邊界條件的新方法,并將其成功應用于實際工程電磁場的數值計算。目前,小波基函數求解偏微分方程時邊界處理問題、小波變換電磁場快速算法問題等是小波分析應用于電磁場工程領域的研究熱點。
2.1.6混合算法
利用不同算法各自的優勢,將不同算法應用于求解同一電磁場問題的混合算法也是電磁場數值計算的發展方向。對此,國內外學者進行了一系列研究,Braess等[23]嘗試在FEM法的框架中引入多重網格法以提高問題求解效率。王海彬等[24]嘗試將MOM與FDTD法混合應用于計算脈沖照射電磁場中人體頭部不同部位SAR值。李莉等[25]將FEM與邊界元法混合應用于求解包含運動導體的二維電磁場計算。
2.1.7基于硬件技術的加速算法
基于GPU(Graphic Processing Unit,GPU)的加速算法圖形處理器是用于圖像信息處理與生成的專用處理器。與CPU相比,GPU更適用于流處理計算。GPU性能不斷增強,其內部的可編程浮點單元運算能力提高,使其適于處理具有“計算密集型”特點的問題[26]。在電磁學領域,為提高電磁場計算效率,開發各種基于GPU的加速算法已成為研究熱點。Krakiwsky等[27]嘗試利用GPU對FDTD算法進行加速。我國學者楊正龍等[28]在物理光學和物理繞射理論基礎上,開發基于GPU的射線追蹤算法,該方法可用于快速計算具有凹腔結構目標的電磁散射問題。
基于云計算的加速算法云計算是隨著互聯網發展產生的新型商業化計算服務模式。云計算具有使用維護成本低、計算資源利用率高等優點。基于云計算的加速算法已應用于電磁場計算中,楊亮[29]將云計算與并行算法相結合,將電磁場空間劃分為多個子空間,使用多個CPU并行計算各個子空間內電磁場,將全部計算結果合并得到整個空間的電磁場分布。金亮[30-31]將云計算技術引入電磁場耦合問題計算中,搭建了基于OpenMPI方法的云計算平臺,并在該平臺上實現電磁場高效并行計算。
2.2.1人體模型研究進展
上世紀90年代起,已有Slavin等[32]開始研制人體解剖學數字化模型并將其應用于醫學仿真。早期的人體數字化模型大都基于志愿者和病人的CT掃描圖像建立,主要用于醫用輻射物理學的相關研究。由于這類模型對軟組織分辨率不高,因此,不適于對人體不同介電性質的軟組織進行仿真。為此,Nagaoka等[33]通過磁共振成像建立了日本成年男性和成年女性的高分辨率全身模型,模型包括51個不同組織,可用于射頻電磁場下的人體仿真。Lee等[34]統計了年齡在18~24歲之間韓國成年男性的平均體格數據,并選擇一名男性志愿者通過其磁共振圖像建立包含30種不同組織結構的人體模型。Christ等[35]基于高分辨率磁共振圖像建立了一組稱為“Virtual Family”的解剖學人體模型(其中包括一名成年男性、一名成年女性以及兩名兒童),這組模型具有80多種不同組織結構,可用于各類電磁環境暴露仿真評估。Szczerba等[36]基于幾何圖像拓撲改進“Virtual Family”人體模型解決其圖形質量和一致性問題。Wu等[37]研究建立了可用于電磁環境仿真的基于中國人體質的解剖學人體模型,包括一名成年男性(含87種組織結構)和一名成年女性(含90種組織結構)。隨后,Li等[38]利用磁共振成像建立了嬰兒(12個月)模型(含32種組織結構)。
2.2.2仿真實驗研究進展
為研究不同磁場強度、不同類型線圈、不同掃描序列條件下,人體組織中SAR分布情況,研究人員采用前文所介紹的電磁場數值計算方法,結合計算機技術利用Mimics等專用軟件,對線圈和人體組織進行計算機建模仿真研究。Wolf等[39]通過仿真實驗,給出可用于高場強MRI設備SAR仿真的人體模型簡化方法。Cornelis等[40]用FDTD算法仿真計算出3 T MRI設備上使用體線圈掃描人體時B1場分布,并將仿真結果與實際測量值比較發現兩者之間存在相關性,驗證該仿真方法可行。Ibrahim等[41]采用FDTD法對鳥籠線圈下人體頭部的SAR值進行仿真計算,結果表明SAR值分布與人體組織特性有關。Neufeld等[42]通過仿真“Virtual Family”人體模型在3 T場強MRI設備上,分別使用多源發射體線圈和傳統體線圈掃描時,模型內部各組織SAR分布情況,研究多源發射體線圈的安全性。仿真結果表明目前應用于強場MRI設備的多源發射體線圈可導致更多熱量沉積,存在更大的安全風險。
曾雁冰等[43]基于仿真方法建立鳥籠線圈與女性盆腔模型,利用FDTD算法計算該模型在64 Hz射頻條件下的SAR值分布,找出局部熱損傷風險較高的組織。黃綺華等[44]利用FDTD算法對磁場強度分別為1.5 T、3 T和7 T的MRI設備掃描人體時,B1場均勻性和SAR分布進行仿真計算,研究人體內部B1場均勻性、SAR分布隨MRI設備磁場強度的變化規律。仿真結果表明:當MRI設備場強從1.5 T分別增加到3 T和7 T時,B1場的均勻性逐漸降低,人體各組織SAR的平均值逐漸增高。
隨著我國經濟發展和醫療服務需求的持續增長,各類醫療機構中MRI設備裝機量不斷增加。MRI設備數量的增多、新設備的投入以及大量使用多年的在用設備,這些都對MRI設備使用中安全風險防護提出更高的要求。如何準確計量SAR并對MRI設備的安全性進行評估,降低受檢者的熱損傷幾率,是安全風險防護中的重點,也是醫生和受檢者普遍關注、需要解決的問題。對MRI設備進行科學有效的計量溯源是解決上述問題的有效方法。
國務院2014年6月實施的《醫療器械監督管理條例》及藥監總局2016年2月實施的《醫療器械使用質量監督管理辦法》中均規定對MRI設備這類需要定期檢查、檢驗、校準、保養、維護的醫療器械,“應當按照產品說明書的要求進行檢查、檢驗、校準、保養、維護”,這是對醫療機構中在用MRI設備開展周期性計量校準的政策依據。
醫療設備質控以規避醫療設備風險為出發點,以醫療設備全生命周期的質量保障為目標,以技術性檢測為基礎手段,以完整的管理流程為執行標準,以數據收集和分析為持續改進方向和管理工作,是醫院質量管理的重要組成部分。MRI設備屬于醫院醫療設備質控重點關注的大型設備,對其進行周期性計量校準是日常質控的重要技術手段。
(1)磁體系統性能計量需求。磁體系統由主磁體、梯度線圈和射頻線圈組成,作為MRI設備中主要部件其磁場強度數值直接影響設備掃描圖像的真實性和有效性,需準確計量。
(2)成像質量計量需求。MRI設備作為影像設備其信噪比、均勻性、空間線性、空間分辨力、低對比度分辨力等成像質量指標均需計量校準。
(1)技術規范需求。國內部分經濟發達省份針對醫用MRI設備起草了地方計量檢定規程或校準規范,規定了MRI設備磁場強度、成像質量的計量要求及各參數的計量方法。但上述規范中對MRI設備重要安全指標SAR的計量性能未作規定,也未給出科學有效的SAR計量方法,這影響到SAR計量校準工作的開展。
(2) 計量設備需求。SAR計量設備包括通用模型 (受體)及精確測量SAR值的專用測量裝置。通用模型由模擬人體不同部位(頭部、軀干)輪廓的殼體及內部仿真溶液構成,用于模擬實際掃描時受檢者人體組織。SAR值專用測量裝置是測量MRI設備空間中每一點上SAR值的專用測量設備。目前國內缺乏符合中國人生理特征的SAR計量通用模型及SAR測量裝置,無法對受檢者受檢時體內能量沉積和升溫的情況進行有效評估,這影響到SAR計量校準工作的開展。研制MRI設備SAR計量通用模型,搭建SAR測量裝置,實現對MRI設備SAR值的精確計量是目前亟待解決的問題,需要進一步深入研究。
SAR作為評價MRI設備射頻生物熱效應的量化指標成為研究熱點。基于電磁場數值計算方法的計算機仿真實驗是研究SAR的有效方法。對MRI設備SAR的計算仿真、并以有效計量溯源保證仿真數據準確,有助于進一步量化研究相關設備射頻在人體組織內生物熱效應,有效降低受檢者的安全風險。
[參考文獻]
[1] Lauterbur PC.Image formation by induced local interactions:Examples of employing nuclear magnetic resonance[J].Nature,1973:242(5394):190-191.
[2] Edelstein WA,Glover GH,Hardy CJ,et al.The intrinsic signal to noise ratio in NMR imaging[J].Magn Reson Med,1986,3(4):604.
[3] 楊保聯.超高場磁共振人體成像應用研究和醫學前景[J].波譜學雜志,2015,32(4):707-714.
[4] 姜曉強,張鵬高,王成發.射頻技術發展與應用現狀[J].信息通信,2012,(3):206.
[5] 畢帆,王龍辰,李斌.MR射頻能量特定吸收率的計算及其臨床應用[J].中國醫療器械雜志,2014,38(6):423-426.
[6] IEC 60601-2-33,Particular Requirements for the Basic Safety and Essential Performance of Magnetic Resonance Equipment for Medical Diagnosis[S].Intemational Electrotechnical Commission,2013.
[7] FDA CDRH Guidance 793,Guidance for Industry and FDA Staff: Criteria for Signi fi cant Risk Investigations of Magnetic Resonance Diagnostic Devices[S].FDA,2003.
[8] 關海爽,馬文閣.電磁場數值計算新方法的研究[J].遼寧工業大學學報(自然科學版),2006,26(4):230-233.
[9] 趙云芳,王紅霞,竹有章,等.電磁場數值計算方法[A].陜西省物理學會2008年學術年會[C].2008.
[10] Yee KS.Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media[J].IEEE T Antenn Propag,1966,14(3):302-307.
[11] 葛德彪,閏玉波.電磁波時域有限差分方法[M].西安:西安電子科技大學出版社,2005.
[12] Gibson WC.The Method of Moments in Electromagnetics[M].Pmladelphia:Taylor &Francis,2007.
[13] Jin J.The Finite Element Method in Electromagnetics[M].New York:Wiley,1993.
[14] Nayroles B,Touzot G,Villon P.Generalizing the fi nite element method: diffuse approximation and diffuse elements[J].Comput Mech,1992,10(5):307-318.
[15] Lei Y,Friswell MI,Adhikari S.A Galerkin method for distributed systems with non-local damping[J].Int J Solids Struct,2006,43(11-12):3381-3400.
[16] 劉素貞,楊慶新,陳海燕,等.無單元法在電磁場數值計算中的應用研究[J].電工技術學報,2001,16(2):30-33.
[17] 王立鵬,王欣彥,唐任遠.EFG-MG法及其在電磁場數值計算中的應用[J].電工技術學報,2010,25(1):1-5.
[18] 鄧東皋,彭立中.小波分析[J].數學進展,1991,20(3):294-310.
[19] 彭玉華,董曉龍.小波變換在電磁場數值計算中的應用[J].電子學報,1996,24(12):46-52.
[20] Baharav Z,Leviatan Y.Impedance matrix compression using adaptively constructed basis functions[J].IEEE T Antenn Propag,1996,44(9):1231-1238.
[21] Wang G.On the utilization of periodic wavelet expansions in the moment methods[J].IEEE T Antenn Propag,1995,43(10):2495-2498.
[22] 楊仕友,倪米正.小波—伽遼金有限元法及其在電磁場數值計算中的應用[J].中國電機工程學報,1999,(1):56-61.
[23] Braess D,Verfurth R.Multigrid methods for nonconforming fi nite element methods[J].Siam J Numer Anal,1990,27(4):979-986.
[24] 王海彬,牛中奇.MOMTD/FDTD混合法分析線天線與近場媒質的相互作用[J].微波學報,2007,23(3):5-9.
[25] 李莉,潘杰,朱翠霞,等.有限元—邊界元耦合法在二維電磁場數值計算中的應用及優化措施[J].山東師范大學學報(自然科學版),2013,28(3):59-65.
[26] 吳恩華,柳有權.基于圖形處理器(GPU)的通用計算[J].計算機輔助設計與圖形學學報,2004,16(5):601-612.
[27] Krakiwsky SE,Turner LE,Okoniewski MM.Graphics processor unit (GPU) acceleration of finite-difference time-domain(FDTD) algorithm[A].International Symposium on Circuits and Systems[C].New York:IEEE,2004:265-268.
[28] 楊正龍,金林,李蔚清.基于GPU的圖形電磁計算加速算法[J].電子學報,2007,35(6):1056-1060.
[29] 楊亮.云計算與電磁仿真結合的可行性分析[J].電子技術與軟件工程,2014,(2):29-30.
[30] 金亮,邱運濤,楊慶新,等.基于云計算的電磁場有限元計算方法研究[J].中國科技論文,2016,11(11):1310-1314.
[31] 金亮,邱運濤,楊慶新,等.基于云計算的電磁問題并行計算方法[J].電工技術學報,2016,31(22):5-11.
[32] Slavin KV.The visible human project[J].Sur Neur,1997,48(6):638.
[33] Nagaoka T,Watanabe S,Sakurai K,et al.Development of realistic high-resolution whole-body voxel models of Japanese adult males and females of average height and weight, and application of models to radio-frequency electromagnetic- fi eld dosimetry[J].Phys Med Biol,2004,49(1):1-15.
[34] Lee,Ae-kyoung,Choi WY,et al. Development of korean male body model for computational dosimetry[J].ETRI J,2006,28(1):107-110.
[35] Christ A,Kainz W,Hahn EG,et al.The Virtual familydevelopment of surface-based anatomical models of two adults and two children for dosimetric simulations[J].Phys Med Biol,2010,55(2):N23.
[36] Szczerba D,Neufeld E,Zefferer M,et al.Unstructured mesh generation from the Virtual Family, models for whole body biomedical simulations[J].Pro Comput Sci,2010,1(1):837-844.
[37] Wu T,Tan L,Shao Q,et al.Chinese adult anatomical models and the application in evaluation of RF exposures[J].Phys Med Biol,2011,56(7):2075-2089.
[38] Li C,Chen Z,Yang L,et al.Generation of infant anatomical models for evaluating electromagnetic field exposures[J].Bioelectromagnetics,2015,36(1):10-26.
[39] Wolf S,Diehl D,Gebhardt M,et al.SAR simulations for highfi eld MRI: how much detail, effort, and accuracy is needed[J].Magnet Reson Med,2013,69(4):1157-1168.
[40] Cornelis VDB,Bartels LW,Van BB,et al.The use of MR B+1 imaging for validation of FDTD electromagnetic simulations of human anatomies[J].Phys Med Biol,2006,51(19):4735-4746.
[41] Ibrahim TS,Lee R,Baertlein BA,et al.B1 field homogeneity and SAR calculations for the birdcage coil[J].Phys Med Biol,2001,46(2):609.
[42] Neufeld E,Gosselin MC,Murbach M,et al.Analysis of the local worst-case SAR exposure caused by an MRI multi-transmit body coil in anatomical models of the human body[J].Phys Med Biol,2011,56(15):4649.
[43] 曾雁冰,周宇佳,黃綺華,等.MR射頻能量特定吸收率在女性盆腔不同組織中分布規律的研究[J].中國醫學物理學雜志,2012,29(5):3632-3637.
[44] 黃綺華,高勇,辛學剛.高場和超高場MR下人體內B_1場均勻性及SAR隨場強變化規律的研究[J].中國生物醫學工程學報,2013,32(1):21-27.
本文編輯 王婷