基于快速體積分方程和磁共振的人體電磁特性反演方法
【專利摘要】基于快速體積分方程和磁共振的人體電磁特性反演方法，涉及核磁共振成像。結(jié)合玻恩迭代算法或變分玻恩迭代算法或變形玻恩迭代算法、快速傅里葉變換、穩(wěn)定雙共軛梯度算法以及共軛梯度算法，進(jìn)行反復(fù)的正演和反演迭代至結(jié)果收斂，求解人體電磁特性參數(shù)，實(shí)現(xiàn)人體電磁特性參數(shù)磁共振成像的方法，得到人體組織電磁特性參數(shù)分布的二維或三維圖像，呈現(xiàn)組織內(nèi)各區(qū)域電磁特性參數(shù)差異的同時(shí)，也可對其進(jìn)行定量研究，可用于醫(yī)學(xué)上的疾病研究和指導(dǎo)臨床診斷、治療。
【專利說明】
基于快速體積分方程和磁共振的人體電磁特性反演方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明設(shè)及核磁共振成像，尤其是設(shè)及一種基于快速體積分方程和磁共振的人體 電磁特性反演方法。
【背景技術(shù)】
[0002] 1991年化acke第一次提出了基于核磁共振的電特性成像化PT)。電特性成像與核 磁共振一樣都是非侵入式的成像方法，它利用核磁共振的數(shù)據(jù)反演出組織的電特性參數(shù)分 布。磁共振檢測伴隨著熱效應(yīng)，運(yùn)可能給人體帶來傷害，尤其是未來要發(fā)展的高場磁共振檢 巧。，熱效應(yīng)與人體的電特性參數(shù)分布有關(guān)，了解電特性參數(shù)的分布可W提前評估出磁共振 檢測給人體帶來的熱效應(yīng)及其影響，運(yùn)對于核磁共振成像的發(fā)展有很大的意義。癌癥由于 其高死亡率和低治愈率而受到人們的廣泛關(guān)注，有研究表明，癌癥細(xì)胞比正常細(xì)胞的電特 性參數(shù)高很多，比如乳癌細(xì)胞比正常細(xì)胞大200%，而膀脫癌細(xì)胞比正常細(xì)胞大100%，利用 電特性成像就可W更清晰明了地分辨出腫瘤和正常組織，運(yùn)有利于癌癥的早期發(fā)現(xiàn)和及時(shí) 理療。目前，基于核磁共振的電特性成像方法已經(jīng)有多種，比如正交鳥籠線圈法、多通道傳 輸接收法、局部麥克斯韋成像法等等。但是運(yùn)些方法都有一定的局限性，它們都基于麥克斯 韋等式，大多都假設(shè)電特性參數(shù)的分布是局部均勻的，沒有將電磁場邊界條件考慮在內(nèi)，雖 然簡化了計(jì)算，但使得電特性參數(shù)在組織的邊界處會變得不可靠，且對于噪聲也更加敏感。 Balidemaj等人雖提出了利用核磁共振系統(tǒng)測量的數(shù)據(jù)進(jìn)行迭代反演的對比源反演方法， 有效地抑制邊界出現(xiàn)的錯(cuò)誤，但作者只重建了二維模型，且沒有將射頻屏蔽考慮在內(nèi)。
[0003] 21世紀(jì)初定量磁化率成像(QSM)被提出，利用MRI方法獲取磁化率定量圖像，對磁 化率做定量分析，可W為生物醫(yī)學(xué)研究提供一種非創(chuàng)傷性手段，而且對組織內(nèi)部順磁性鐵 含量的測量，有利于對腦部血管疾病和神經(jīng)病變疾病的診斷和治療。QSM需要對相位信息進(jìn) 行解纏繞及去除背景場的預(yù)處理來獲得反映局部磁場變化的場圖，再結(jié)合重建算法重建出 磁化率圖像。對于QSM方法，在病態(tài)逆問題的求解情況下，背景場去除效果不理想的話會對 磁化率的求解產(chǎn)生較大影響，甚至?xí)?dǎo)致最終的磁化率成像結(jié)果不可用，目前兩種去除背 景場較好的方法是復(fù)雜諧波偽影去除法(SHARP)和偶極場投影法(PDF)。由場圖信息重建磁 化率圖像是一個(gè)不適定逆問題，穩(wěn)定求解和精確定量是關(guān)鍵問題。常用的QSM重建方法有多 方向采樣磁化率計(jì)算方法(COSMOS)、貝葉斯正則化方法、k空間加權(quán)微分法(W邸)等。由于成 像儀器腔體限制和病人舒適角度問題，很難獲得多個(gè)不同擺放方向的腦部成像，運(yùn)極大限 審IJ 了 COSMOS方法的臨床應(yīng)用。貝葉斯正則化方法，隨著匹配噪聲參數(shù)增大，磁化率重建圖像 清晰度變好，顯示出很好的去噪能力。WH)方法成像能顯示出腦組織局部細(xì)節(jié)，且沒有出現(xiàn) 明顯的條狀偽影，但成像參數(shù)(如場強(qiáng)、回波時(shí)間、翻轉(zhuǎn)角等)會對相位處理和磁化率分布圖 的精度產(chǎn)生影響。

【發(fā)明內(nèi)容】

[0004] 本發(fā)明的目的在于針對現(xiàn)有技術(shù)中使用的電特性成像化PT)算法在組織邊界處值 不可靠，對于噪聲更加敏感，導(dǎo)致成像分辨率較差，W及磁共振(MR)成像中幅值信息和相位 信息分開成像，且利用相位信息的磁化率成像需要解纏繞和去除背景場處理等問題，提供 一種基于快速體積分方程和磁共振的人體電磁特性反演方法。
[0005] 本發(fā)明包括W下步驟：
[0006] 1)測量實(shí)驗(yàn)數(shù)據(jù)，利用磁共振的影像技術(shù)測量人體，得到人體的Bi+場，所述Bi+場 是磁共振射頻場的正旋場；
[0007] 2)進(jìn)行迭代計(jì)算，具體步驟為：
[0008] ①假定初始的電特性對比度xE;W(r)和磁化率xH;W(r)，利用方程（1)和(2)，進(jìn)行 正演計(jì)算，解出E^(r)和HW(r);
[0009] 所述方程(1)和(2)如下：
[0012]其中，Dinv代表反演區(qū)域，r和r'為反演區(qū)域中的位置，kb為背景介質(zhì)的波數(shù)，ω為 角頻率，j為虛數(shù)符號；分別是背景介質(zhì)的復(fù)介電常數(shù)和磁導(dǎo)率;EinK(r)和E(r)分別為 入射電場和總電場，HinK(r)和H(r)分別為入射磁場和總磁場;g(r，r'）是標(biāo)量格林函，xE(r) 是電特性對比度，xH(r)是磁化率，定義為W下方程：
[0016] 其中，方程(4)和(5)中的和μ0-)分別是人體復(fù)介電常數(shù)和磁導(dǎo)率，而復(fù)介電 常數(shù)又可W用方程(6)表示：
[0017]
[0018] 其中，e(r)，〇(r)分別是人體的介電常數(shù)，電導(dǎo)率；
[0019] ②將第n-1次迭代的總電場E(r)和總磁場H(r)，代入方程(7)和(8)，進(jìn)行反演計(jì) 算，解出第η次迭代的xEiW(r)和xHsWw;
[0020] 所述方程(7)和(8)如下：
[0021]
[0022]
(8)
[0023] 其中，Hsea(r)散射磁場;Hl + isea(r)為正旋散射磁場，定義為化+ isea(r)二化sea(r) + 冊sea(r)]/2;HxSea(r)和HySea(r)分別為散射磁場Hsea(r)的X分量與y分量;Bl+(r)為磁共振正 旋場的磁通密度，BiWne(r)為背景介質(zhì)下磁共振正旋場的磁通密度；
[0024] ③將步驟②解出的χΕ=ω0-)和χΗ=ω0-)，代入方程（1)和（2)，進(jìn)行正演計(jì)算，解出 第η次迭代的Eb>(r)和Η^>0-);
[0025] ④重復(fù)步驟②和③，直到利用乂6=^的，乂11=(。)的代入方程(7)和(8)算出的81+的 與測量值的相對殘差在10% W內(nèi)，則算法收斂，結(jié)束迭代。
[0026] 所述迭代的算法可采用玻恩迭代算法(ΒΙΜ)、變分玻恩迭代算法(VBIM)或變形玻 恩迭代算法(DBIM)等。
[0027] 在步驟2)第①和③步驟中，所述正演計(jì)算為己知電特性對比度xE(r)和磁化率χΗ (r)，利用方程（1)和(2)，結(jié)合穩(wěn)定雙共輛梯度算法(BCGS)和快速傅里葉變換(FFT)，計(jì)算出 總電場E(;r)和總磁場H(;r)。
[0028] 在步驟2)第②步驟中，所述反演計(jì)算為己知總電場E(r)和總磁場H(r)，利用方程 (7)和(8)，結(jié)合共輛梯度算法(CG)和快速傅里葉變換(FFT)，計(jì)算出電特性對比度χΕ(r)和 磁化率xH(r)。
[0029] 3)利用方程(4)，巧），（6)結(jié)合迭代結(jié)果中的電特性對比度xE(r)和磁化率xH(r)，計(jì) 算出e(r)，〇(r)，y(r)(或磁化率xH(r))的分布，輸出結(jié)果，計(jì)算完畢。
[0030] 本發(fā)明具體設(shè)及電和磁體積分方程，結(jié)合玻恩迭代算法(BIM)或變分玻恩迭代算 法(VBIM)或變形玻恩迭代算法(DBIM)、快速傅里葉變換(FFT)、穩(wěn)定雙共輛梯度算法(BCGS) W及共輛梯度算法(CG)，進(jìn)行反復(fù)的正演和反演迭代至結(jié)果收斂，求解人體電磁特性參數(shù) (包括電導(dǎo)率0,介電常數(shù)ε和磁導(dǎo)率μ)，實(shí)現(xiàn)人體電磁特性參數(shù)磁共振成像的方法，得到人 體組織電磁特性參數(shù)分布的二維或Ξ維圖像，呈現(xiàn)組織內(nèi)各區(qū)域電磁特性參數(shù)差異的同 時(shí)，也可對其進(jìn)行定量研究，可用于醫(yī)學(xué)上的疾病研究和指導(dǎo)臨床診斷、治療。
[0031] 本發(fā)明能有效抑制邊界出現(xiàn)的錯(cuò)誤、減小噪聲影響、提高成像分辨率，結(jié)合MR信號 的幅值信息和相位信息，利用快速傅里葉變換(FFT)加快重建速度、減少重建時(shí)間，結(jié)合玻 恩迭代算法(BIM)或變分玻恩迭代算法(VBIM)或變形玻恩迭代算法(DBIM)，穩(wěn)定雙共輛梯 度算法(BCGS)，共輛梯度算法(CG)，同時(shí)反演電特性參數(shù)(電導(dǎo)率0,介電常數(shù)ε)和磁特性參 數(shù)(磁導(dǎo)率μ)。
[0032] 本發(fā)明的有益技術(shù)效果如下：
[0033] 提出一種基于快速體積分方程算法和磁共振的人體電磁特性反演方法，對于電特 性成像，能夠有效抑制邊界出現(xiàn)的錯(cuò)誤、減小噪聲影響、提高成像分辨率，利用快速傅里葉 變換(FFT)加快重建速度、減少重建時(shí)間，結(jié)合變分玻恩迭代算法(ΒΙΜ)或變分玻恩迭代算 法(VBIM)或變形玻恩迭代算法(DBIM)，穩(wěn)定雙共輛梯度算法(BCGS)，共輛梯度算法(CG)，同 時(shí)反演電特性參數(shù)(電導(dǎo)率0、介電常數(shù)0和磁特性參數(shù)(磁導(dǎo)率μ)。
【附圖說明】
[0034] 圖1是反演出的人腦相對介電常數(shù)Ξ維分布。
[0035] 圖2是反演出的人腦電導(dǎo)率Ξ維分布。
[0036] 圖3是反演出的人腦磁化率Ξ維分布。
[0037] 圖4是反演出的人腦相對介電常數(shù)分布Χ-Υ截面。
[0038] 圖5是反演出的人腦電導(dǎo)率分布Χ-Υ截面。
[0039] 圖6是反演出的人腦磁化率分布Χ-Υ截面。
【具體實(shí)施方式】
[0040] 本發(fā)明利用快速體積分方程算法和磁共振數(shù)據(jù)進(jìn)行人體的二維或Ξ維電磁特性 反演，其中快速體積分方程算法中結(jié)合了玻恩迭代算法(ΒΙΜ)或變分玻恩迭代算法(VBIM) 或變形玻恩迭代算法(DBIM)，快速傅里葉變換(FFT)，穩(wěn)定雙共輛梯度算法(BCGS)，共輛梯 度算法(CG)。
[0041 ]本實(shí)施例W變分玻恩迭代算法(VBIM)為例子進(jìn)行解釋?！揪唧w實(shí)施方式】如下：
[0042] 1)測量實(shí)驗(yàn)數(shù)據(jù)。利用磁共振的影像技術(shù)測量人體，得到Βι+場。Βι+場是磁共振射 頻場的正旋場。
[0043] 2)進(jìn)行迭代計(jì)算，包括W下步驟：
[0044] ①假定初始的電特性對比度xEiW(r)和磁化率xHi(〇)(r)，進(jìn)行正演計(jì)算，解出eW (r)和 H(〇)(r)。
[0045] ②利用第n-1次迭代的總電場E(r)和總磁場H(r)，進(jìn)行反演計(jì)算，解出第η次迭代 的 xE;(n)(r)和 xH;(n)(r)。
[0046] ③再由步驟②解出的χΕ=ω0-)和χΗ=ω(')，進(jìn)行正演計(jì)算，解出第η次迭代的E(n> (r)和 H(n)(r)。
[0047] ④重復(fù)步驟②和③，直到xEiW(r)，xHiW(r)代入方程（14)和（15)算出的Bi+(r)與 測量值的相對殘差在10% W內(nèi)，則算法收斂，結(jié)束迭代。
[0048] 步驟①和③所述正演計(jì)算為己知電特性對比度xE(r)和磁化率xH(r)，利用W下方 程計(jì)算總電場E(r)和總磁場H(r):
[0053]其中，Dinv代表反演區(qū)域，r和r'為反演區(qū)域中的位置，kb為背景介質(zhì)的波數(shù)，ω為 角頻率，j為虛數(shù)符號。分別是背景介質(zhì)的復(fù)介電常數(shù)和磁導(dǎo)率。Eine(r)和E(r)分別為 入射電場和總電場，HinK(r)和H(r)分別為入射磁場和總磁場。A(r)和F(r)分別為磁勢矢量 和電勢矢量。g(r，r'）是標(biāo)量格林函數(shù)，xE(r)是電特性對比度，xH(r)是磁化率，定義為W下 方程：
[0057]其中，方程(4)和(5)中的e^r)和μ0-)分別是人體復(fù)介電常數(shù)和磁導(dǎo)率，而復(fù)介電 常數(shù)又可W用方程(6)來表示：
[0化引
[0059] 其中，e(r)，〇(r)分別是人體的介電常數(shù)，電導(dǎo)率。
[0060] 便于計(jì)算，利用沖激函數(shù)作為基函數(shù)和測試函數(shù)，則方程（1)，（2)，（3)和(4)被離 散化為下列方程：
[0067] je[i，j]，ke[i，K]，le[i，L]
[006引 A X、Δ y和Δ Z分別為單元格在x、y和Z方向的大小。j，k，l分別為離散后X方向，y方 向和Z方向的位置。J，K，L為X方向，y方向和Z方向的離散點(diǎn)數(shù)。
[0069 ] 而方程(11)和(12)中的Δ V為單元格體積，定義為Αν=ΔχΑγΔζ。
[0070] 由于，標(biāo)量格林函數(shù)g(r，r'）具有平移不變性，所W方程（11)和（12)可W利用快速 傅里葉變換(FFT)進(jìn)行快速計(jì)算。此方法可W節(jié)省大量時(shí)間W及內(nèi)存，使計(jì)算更加快速。
[0071] 離散化后，可W將方程寫成W下形式：
[0075] L代表的是線性操作。
[0076] 直接解方程（13)會花大量時(shí)間和內(nèi)存，所W在本發(fā)明當(dāng)中利用穩(wěn)定雙共輛梯度算 法(BCGS)和快速傅里葉變換(FFT)計(jì)算方程(13)得出總電場E(r)和總磁場H(r)。
[0077] 步驟②所述反演計(jì)算為己知總電場E(r)和總磁場H(r)，利用W下方程計(jì)算電特性 對比度xE(r)和磁化率xH(r):
[0080] 其中，Hsea(r)散射磁場。Hl + ;sea(r)為正旋散射磁場，定義為Hl + ;sea(r) =化sea(r) + 冊sKa(r)]/2,其中，Η嚴(yán)(r)和HySKa(r)分別為散射磁場Hsc^a(r)的X分量與y分量。Bi+(r)為磁 共振正旋場的磁通密度，BiWne(r)為背景介質(zhì)下磁共振正旋場的磁通密度。GbE(r，r'）和GbH (r，r'）分別為電并矢格林函數(shù)和磁并格林函數(shù)。
[0081] 本發(fā)明中利用變分玻恩迭代來計(jì)算物體的電磁特性分布。所W方程（14)寫成W下 形式：
[0082]
[0083] 其中，
[0084] 5Hsca(r)=Hsca(r)-Hsca;(n-i)(r) (17)
[0085] 代表的是測量值與n-1次迭代的散射磁場的差值。
[0088] 代表的是第η次迭代與第n-1次迭代計(jì)算的電特性對比度xE(r)和磁化率xH(r)的差 值。
[0089] 為了便于計(jì)算，利用沖激函數(shù)作為基函數(shù)和測試函數(shù)，則方程(16)寫為W下形式：
[0090]
[0091] 其中，電并矢格林函數(shù)GbE(r，r'）和磁并格林函數(shù)GbH(r，r'）都具有平移不變性，所 W也可W利用快速傅里葉變換(FFT)來進(jìn)行快速計(jì)算。
[0092] 當(dāng)利用方程(20)計(jì)算δχΕ0-)和δχΗ0-)時(shí)，如果直接用矩陣方式進(jìn)行計(jì)算的話，花費(fèi) 的時(shí)間與內(nèi)存過于龐大，所W本發(fā)明利用共輛梯度算法(CG)，結(jié)合快速傅里葉變換來計(jì)算δ xE(r)和 δχΗ0-)，然后利用方程（18)和（19)算出 xEiW(r)和 χΗ=ω^)。
[0093] 3)利用方程(6)，（7)，（8)結(jié)合迭代結(jié)果中的電特性對比度xE(r)和磁化率xH(r)，計(jì) 算出e(r)，〇(r)，y(r)(或磁化率xH(r))的分布，輸出結(jié)果，計(jì)算完畢。
[0094] 表 1
[0095]

[0096] 人腦設(shè)置參數(shù)如表1，頻率為63.87MHz。
【主權(quán)項(xiàng)】
1.基于快速體積分方程和磁共振的人體電磁特性反演方法，其特征在于包括以下步 驟： 1) 測量實(shí)驗(yàn)數(shù)據(jù)，利用磁共振的影像技術(shù)測量人體，得到人體的B1+場，所述B1+場是磁共 振射頻場的正旋場； 2) 進(jìn)行迭代計(jì)算，具體步驟為： ① 假定初始的電特性對比度xM())(r)和磁化率xH;(())(r)，利用方程（1)和(2)，進(jìn)行正演 計(jì)算，解出 E(())(rWPHw(r); 所述方程(1)和(2)如下：其中，Dinv代表反演區(qū)域，r和r'為反演區(qū)域中的位置，kb為背景介質(zhì)的波數(shù)，ω為角頻 率，j為虛數(shù)符號；分別是背景介質(zhì)的復(fù)介電常數(shù)和磁導(dǎo)率;Eine(r)和E(r)分別為入射 電場和總電場，H inIr)和H(r)分別為入射磁場和總磁場;g(r，r'）是標(biāo)量格林函，xE(r)是電 特性對比度，x H(r)是磁化率，定義為以下方程：其中，方程(4)和(5)中的^(r)和μ(Γ)分別是人體復(fù)介電常數(shù)和磁導(dǎo)率，而復(fù)介電常數(shù) 又可用方程(6)表示：其中，Kr)，〇(r)分別是人體的介電常數(shù)，電導(dǎo)率； ② 將第n-1次迭代的總電場E(r)和總磁場H(r)，代入方程(7)和(8)，進(jìn)行反演計(jì)算，解 出第η次迭代的xE;(n)(r)和x H;(n)(r); 所述方程(7)和(8)如下：(?其中，Hsea(r)散射磁場;^+^(r)為正旋散射磁場，定義為出+;^(〇 =[舊。^)+犯/。3 (r) ]/2;Hxs°a(r)和Hys°a(r)分別為散射磁場H s°a(r)的X分量與y分量;Bi+(r)為磁共振正旋場 的磁通密度，B廣ine(r)為背景介質(zhì)下磁共振正旋場的磁通密度； ③ 將步驟②解出的xE;(n)(r)和xH;(n)(r)，代入方程（1)和(2)，進(jìn)行正演計(jì)算，解出第η次 迭代的 E(n)(r)和 H(n)(r); ④ 重復(fù)步驟②和③，直到利用乂￡;(")&)^11;(11)(〇代入方程(7)和(8)算出的8 1+(〇與測 量值的相對殘差在1 〇 %以內(nèi)，貝1J算法收斂，結(jié)束迭代； 所述迭代的算法采用玻恩迭代算法、變分玻恩迭代算法或變形玻恩迭代算法； 3)利用方程(4)，（5)，（6)結(jié)合迭代結(jié)果中的電特性對比度xE(r)和磁化率xH(r)，計(jì)算出 ε(Γ)，σ(Γ)，μ0-)(或磁化率x H(r))的分布，輸出結(jié)果，計(jì)算完畢。2. 如權(quán)利要求1所述基于快速體積分方程和磁共振的人體電磁特性反演方法，其特征 在于在步驟2)第①和③步驟中，所述正演計(jì)算為己知電特性對比度x E(r)和磁化率xH(r)，利 用方程（1)和(2)，結(jié)合穩(wěn)定雙共輒梯度算法和快速傅里葉變換，計(jì)算出總電場E(r)和總磁 場H(r)。3. 如權(quán)利要求1所述基于快速體積分方程和磁共振的人體電磁特性反演方法，其特征 在于在步驟2)第②步驟中，所述反演計(jì)算為己知總電場E(r)和總磁場H(r)，利用方程(7)和 (8)，結(jié)合共輒梯度算法和快速傅里葉變換，計(jì)算出電特性對比度x E(r)和磁化率xH(r)。
【文檔編號】A61B5/055GK105877747SQ201610190456
【公開日】2016年8月24日
【申請日】2016年3月30日
【發(fā)明人】柳清伙, 洪榮漢, 李勝男, 劉娜, 張建華
【申請人】廈門大學(xué)