脈沖梯度磁場由梯度線圈產生，廣泛用于磁共振成像設備的信號定位系統中。但該脈沖也會同時在線圈周圍導體結構中感應出渦流場，從而導致總磁場改變，使圖像出現變形或偽影，影響臨床診斷。本文通過實驗測量了 1.5 T 超導磁共振成像設備中脈沖梯度導致的渦流場；根據電感-電阻串聯模型，分析了其數學模型，該模型表明渦流磁場可以由獨立項和交叉項組成；然后本文定量測定了每項指數函數的振幅和時間常數等關鍵參數，為設計數字濾波器以進行預加重和波形整形提供了依據。本文涉及的測量是在設計的專用工裝上進行，它由基于中心對稱的六個方向上的水模和各自的射頻接收線圈構成，通過加載測試序列可以得到有無渦流的原始一維數據，經計算可將其轉化為有意義的采樣數據。根據本文研究提及的數學模型，將采樣數據用最小二乘法擬合后表明，數據符合電感-電阻串聯模型，渦流效應確實存在。通過對比實驗同時也進一步證實，用數字濾波器對脈沖進行預加重后，渦流場被有效減小。本文建立了系統函數的數學模型，驗證了預加重在減小渦流上的實用價值，提出了改進辦法，為磁共振設備減輕渦流效應提供了參考。
引用本文： 何汶靜, 祝元仲, 王文周, 鄒凱, 張凱, 賀超. 磁共振梯度渦流場的定量實驗與分析. 生物醫學工程學雜志, 2017, 34(2): 220-226. doi: 10.7507/1001-5515.201506059 復制
引言
隨著磁共振成像（magnetic resonance imaging，MRI）設備在國內醫院的逐漸普及，通過 MRI 在盡量短的時間內獲得高分辨率、高信噪比、低偽影的高質量圖像成為未來醫學診斷領域的發展方向，其中脈沖梯度磁場技術是空間編碼、彌散加權成像的關鍵技術，同時也是制約縮短成像時間的主要問題之一。脈沖梯度磁場由置于 MRI 設備中的梯度線圈產生，應用在幾乎所有的成像序列中[1]。由于在產生脈沖梯度的同時也會感應出渦流（eddy cur-rents，EC），從而會導致回波信號的缺失和變形，引起偽影。具體來說，當梯度線圈由一個梯形波驅動時，會在其反方向感應出 EC，也就是在梯形的上升沿或下降沿產生相反的效應[2]，使上升或下降時間持續幾百微秒。雖然理論上可采用延長時間并等待 EC 自行消失的辦法以減輕偽影，但由于成像過程中需要幾百次的梯度線圈開關，從而使得實際上等待 EC 自行消失的辦法操作性較低，而某些需要短脈沖梯度磁場的序列更是無法實現。
由于 EC 是由 MRI 設備中多個零部件共同引起的，所以找到 EC 產生的具體部分并不具有可操作性，即從理論上找出源頭、建模、解決問題的方式不合適。因此采用直接定量測定 EC 的大小，并找到具有操作性的辦法從工程上來說比較實際。為此，國內外提出了各種測量和降低 EC 的辦法，例如改變梯度序列[3]，改變射頻脈沖傳輸方式[4]，或者預加重波形補償等辦法[5-7]。目前尚未見關于提取實驗數據和關于如何實際應用這些數據獲得波形的報道。基于已有研究成果，本文應用先驅的實驗方法在超導 MRI 設備的梯度線圈上進行實驗，獲得實驗數據；進一步推導了 EC 產生的系統函數，計算出預加重波形的形狀；并提出利用迭代方法，可進一步減輕偽影。在本文的具體測試過程中，首先根據數學模型編寫測量序列，然后根據在 MRI 設備上的測試實驗完成數據采樣，再將采樣數據在 Matlab 平臺中處理，計算出關鍵參數和校正波形，最后將校正波形在譜儀中實現，并驗證實驗結果。
本文運用預加重波形補償的辦法，詳細描述了實驗過程中如何在設備上提取實驗數據和如何應用這些數據獲得預加重波形以有效減小渦流，并在原有基礎上提出預加重波形的推導過程和改進辦法，為國內磁共振設備研究企業和機構提供一點參考。
1 理論基礎
1.1 不均勻磁場的求解
測試從推導求解不均勻磁場理論分布的表達式開始，尋找緩解 EC 造成的偽影的方法。考慮一個由梯度線圈產生的梯度磁場，EC 將在其上升或下降沿產生，因此采用定量 EC 大小的方法來處理。我們通過在梯度線圈上施加一個較長持續時間的階躍波形，當梯度關閉的時候，再施加一個 90° 射頻脈沖用以在水模中產生自由衰減信號（free induction decay signal，FID）[6]。然后立即用接收線圈接收，經模擬/數字（analogue to digital，A/D）轉換獲得一維采樣數據，該 FID 只在 EC 背景下采得，因此可以反映出 EC 的大小。若假設上升沿和下降沿的 EC 相等，則上升沿情況可同理推出。
在下降沿，Bloch 方程描述了時變函數磁矢量 M 和加載的磁場強度 B 的關系[8]：
| $\frac{{{\rm{d}}{{M}}}}{{{\rm{d}}t}} \!=\! {\rm{\gamma }}{{M}} \! \times \! {{B}} \!-\! \frac{{{M_x}{{X}} \! + \! {M_y}{{Y}}}}{{{{\rm{T}}_2}}} \!-\! \frac{{({M_z} \!-\! \left| {{{{M}}_0}} \right|){{Z}}}}{{{{\rm{T}}_1}}}$ |
其中 M 0 是僅有主磁場 B 0 時沿 Z 軸的縱向磁矢量，理論上只有 Z 方向的分量，T1 和 T2 分別表示橫向和縱向弛豫時間。
假設旋轉坐標頻率等于拉莫爾進動頻率，則式（1）右邊第一項變為
| $\left\{ \begin{array}{l}\frac{{{\rm{d}}{M_z}}}{{{\rm{d}}t}} = - \frac{{({M_z} - \left| {{{{M}}_0}} \right|)}}{{{{\rm{T}}_1}}}\\\frac{{{\rm{d}}{{{M}}_{xy}}}}{{{\rm{d}}t}} = - \frac{{{{{M}}_{xy}}}}{{{{\rm{T}}_2}}}\end{array} \right.$ |
由 90° 射頻脈沖產生的橫向磁化矢量的 FID 信號可由式（2）解出
| $\left\{ \begin{array}{l}{{{M}}_{xy}} = {M_x} + {{i}}{M_y} = \left| {{{{M}}_0}} \right|{{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_{\rm{2}}}}}}}{{\rm{e}}^{ - {{i}} { \text{ω}} {}_0t}}\\{M_z} = \left| {{{{M}}_0}} \right|(1 - {{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_{\rm{1}}}}}}})\end{array} \right.$ |
假設 B （ r ）是實驗室坐標下 r 處 EC 的感應磁場，它感應出的接收線圈中的電壓V（t）為
| $\begin{aligned}& V(t) = \! - \! \frac{\partial }{{\partial t}}\int_{{\rm{object}}}^{} \! {{{B}}({{r}}) \! \cdot \! {{M}}({{r}},t){\rm{d}}r} \\ &= - \! \int_{{\rm{object}}} \!\!\! {[{B_{r,x}}({{r}}) \!\! \frac{{\partial {M_x}({{r}},t)}}{{\partial t}} \!\! + \! {B_{r,y}}({{r}}) \! \frac{{\partial {M_y}({{r}},t)}}{{\partial t}}]} {\rm{d}}r\end{aligned}$ |
其中，Br，x（ r ）、Br，y（ r ）是 B （ r ）在 r 處的 X 、 Y 方向的分量，Mx （ r ，t）、My （ r ，t）是 M （ r ，t）在 r 處的 X 、 Y 方向的分量。由于Mz （ r ，t）相對于 |Mxy （ r ，t）| 變化緩慢，所以Br，z（ r ）Mz （ r ，t）被忽略。將式（3）代入式（4）中，簡化后得
| $\begin{aligned}{V_{\sin }}(t) = & - \int_{{\rm{object}}} { {{ω}} ({{r}}){B_{r,xy}}({{r}})\left| {{{{M}}_0}} \right|} \cdot \\ &{{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_2}}}}}\sin [ {{ω}} ({{r}}) + {{ \text{φ}}_0}]{\rm{d}}{{r}}\end{aligned}$ |
或
| $\begin{aligned}{V_{\cos }}(t) = & \int_{{\rm{object}}} { {{ω}} ({{r}}){B_{r,xy}}({{r}})\left| {{{M}}_0} \right|} \cdot \\ &{{\rm e}^{ - \frac{t}{{{T_2}}}}}\cos [ {{ω}} ({{r}}) + {{ \text{φ}}_0}]{\rm{d}}{{r}}\end{aligned}$ |
其中用ω（ r ）代替了式（3）中的ω0，表示此處ω 是位置的函數。最后用相位檢測器（phase-sensitive detection，PSD）獲得信號 S （t），其中實部SR （t） 等于Vsin（t），虛部SI （t） 等于Vcos（t）：
| ${{{S}}(t) = {S_R}(t) + i{S_I}(t)$ |
如果不均勻磁場 ΔB（τ）是和時間相關的函數，則式（5）、式（6）中的相位可以寫為
| $\phi = {\text{ω}} ({{{r}}}) + {{φ}_0} = {\rm{\gamma }}\int_0^t {\Delta B({ \text{τ}} ){\rm{d}}{\text{τ}} } + \Delta {\text{ω}} t$ |
將式（5）和式（6）代入式（7），可以得到
| ${{S}}(\!t\!) \!=\!\!\! \int_{{\rm{object}}}^{} \!\!\!\!{ { \text{ω}} (\!{{r}}\!){{{B}}_{r,xy}}(\!{{r}}\!) \! \left| {{{{M}}_0}} \right| \! {{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_{\rm{2}}}}}}}{{\rm{e}}^{ - {{i}}[{\rm{\gamma }}\int_0^t {\Delta B(\! { \text{τ}} \!){\rm d}{ \text{τ}} + \Delta { \text{ω}} t} ]}}{\rm{d}}{{r}}} $ |
式中不均勻磁場 ΔB（τ）是需要求解的量，包含在信號 S （t）的相位信息中，所以在后續的計算中僅需要關注相位信息，即式（5）、式（6），又式（5）、式（6）中 φ0 是初相位，為常數，所以相位 的變化量僅由不均勻磁場 ΔB（τ）積分造成，考察 S （t）的相位變化量即可求得 ΔB（τ）， S （t）可根據式（7）用 PSD 測量獲得。出于習慣，將變量τ 用t 代換，假設 ΔB（t）由以下兩項構成
| $\Delta B(t) = {B_0}(t) + G(t)x$ |
其中B0（t）是偽空間場，而G（t）是由 EC 引起的偽梯度場，x 是在 X 軸上的位置。
由于同一磁體系統的B0（t）、G（t）僅為時間函數，在 X 軸不同位置x1、x2 處上放置水模，即可獲得不同的 ΔB1（t）和 ΔB2（t），其表達式為：
| $\left\{ \begin{array}{l}\Delta {B_1}(t) = {B_0}(t) + G(t){x_1}\\\Delta {B_2}(t) = {B_0}(t) + G(t){x_2}\end{array} \right.$ |
式（11）兩式相減，即可計算G（t）的表達式。兩式分別乘以x2、x1，再相減，可計算B0（t）的表達式：
| $\left\{ \begin{array}{l}G(t) = \frac{{\Delta {B_1}(t) - \Delta {B_2}(t)}}{{{x_1} - {x_2}}}\\{B_0}(t) = \frac{{{x_1}\Delta {B_2}(t) - {x_2}\Delta {B_1}(t)}}{{{x_1} - {x_2}}}\end{array} \right.$ |
為了求解B0（t）、G（t），x1、x2 可直接測量水模位置得出；ΔB1（t）和 ΔB2（t）可在得到采樣信號 S （t）后，根據式（7），信號相位 對時間的微分得出。
1.2 測試和擬合
根據以上理論推導，可得為解出G（t）和B0（t）的關鍵是要得到信號相位 ，而 是通過式（7）反正切計算得出： ，為此測試序列至少應包含：① X 方向梯度線圈加載階躍電流信號i（t）=u（t）；② 通過體線圈輸出 90° 射頻脈沖；③ 在回波時間（echo time，TE）后，從接收線圈接收信號，獲得式（7）中的一組一維數據 S （t）；但是，由于 EC 的持續時間為 105 s 量級，即需要考察的 S （t）應持續至少 1 s，而單次 FID 由于衰減較快，采樣能獲得的有效信號時間長度僅為 10–2 s 量級，所以必須采用多次采集數據，最后將各段數據拼接起來的方式獲得整個 S （t）[6]，所以測量的第 ③ 步應該改為：多次改變 TE 的值，獲得多組數據 Si （t），拼接獲得 S （t）。
獲得 S （t）后，即可根據式（12）計算出B0（t）、G（t），至此，B0（t）、G（t）的原始數據已全部獲得。但由于其中包含大量隨機噪聲，還需用最小二乘法進行曲線擬合。根據電感—電阻串聯（induc-tance-resistance，L-R）模型擬合G（t）-t 曲線。由 L-R 模型，G（t）-t 曲線可以描述為多指數項的和[8]，如式（13）所示：
| $G(t) = \sum\limits_{k = 1}^n {{{\rm{a}}_k}{{\rm{e}}^{ - \frac{t}{{{{\rm{{ \text{τ}} }}_k}}}}}} $ |
其中 ak 是常系數，τk 是時間常數，n 是指數項的個數。最小二乘法可用作擬合算法[9]。
同理可推出B0（t）-t 曲線。
應當指出的是，此處得到的G（t）和B0（t）都是 X 方向的梯度Gx（t）用 X 方向放置的水模測出的 X 方向 EC 的感應磁場，因此應將其寫為Gx2x（t）和B0x2x（t）。實際上為了更準確的測量，還應測量交叉項，即Gx2y（t）和Bx2y（t）、Gx2z（t）和Bx2z（t）、Gy2z（t）和By2z（t）等。
1.3 波形校正
將形成 EC 效應的整體考慮為一個系統H，如圖 1 所示，輸入為x（t），輸出為y（t）。
圖1
EC 效應的系統
Figure1.
System of EC
x（t）經傅里葉變換為X（jω），y（t）經傅里葉變換為Y（jω），則如圖 1 所示系統函數可表示為：
| $Y({\rm{j}} { \text{ω}} ) = X({\rm{j}} { \text{ω}} ) \cdot H({\rm{j}} { \text{ω}} )$ |
根據要求，波形校正后輸出信號y（t）是理想的方波信號 r（t），傅里葉變換為 R（jω），則如果已知系統函數H（jω），便可根據式（14）得到輸入函數x（t）的傅里葉變換X（jω），再進行反變換，即可得到輸入函數x（t），實現波形校正。
| $X({\rm{j}} { \text{ω}} ) = \frac{1}{{H({\rm{j}} { \text{ω}} )}} \cdot {\rm{R}}({\rm{j}} { \text{ω}} )$ |
問題變為求解系統函數H（jω）。為此，根據小節 1.2 中介紹的測試和擬合方法，可得一對已知條件：輸入x（t）為階躍信號 u（t），輸出y（t）為式（13）計算后得到的擬合曲線G（t），由式（14）可知，將x（t）、y（t）分別進行傅里葉變換，可最后求得系統函數。
同理，根據階躍函數和矩形脈沖函數的關系，也可以推導另一對已知條件，同樣可求出系統函數：如果輸入x（t）為矩形波 r2T（t），r2T是高度為 1 的矩形函數，寬度為（–T，T），則在假設矩形波上升沿和下降沿產生的 EC 大小相等，方向相反的情況下，由 L-R 模型，矩形波 r2T（t）引起的系統響應G（t）可以表示為
| $\begin{aligned}G(t) = & (E - \sum\limits_{k = 1}^n {{{\rm{a}}_k}{{\rm{e}}^{ - \frac{{t + T}}{{{{ \text{τ}}_k}}}}}} ){{\rm{r}}_{2T}} + \\ &E \sum\limits_{k = 1}^n {{{\rm{a}}_k}{{\rm{e}}^{ - \frac{{t - T}}{{{{ \text{τ}}_k}}}}}} {\rm{u}}(t - T)\end{aligned}$ |
其中E 表示輸入方波x（t）的最大值；u（t）為高度為 1 的階躍函數。
2 實驗
本文實驗在一臺醫用 1.5 T、710 mm 孔徑的超導 MRI（四川成都，奧泰醫療）上進行。施加的梯度為 12.5 mT/m。測試工裝在 X 、 Y 、 Z 三個軸的正負對稱六個位置分別放置一個直徑為 2 mm 的球形水模，并將其置于磁體孔徑的中心位置。每個水模距離中心位置的距離為 10 cm，且包裹接收線圈。
脈沖序列設計為梯度線圈加載的階躍電流，信號持續時間為 1 s，每一個重復周期（repeat time，TR）顛倒梯度方向，即改變階躍電流方向；階躍信號停止后經等待時間（time delay，TD）施加 90° 射頻脈沖，TD 的變化范圍為 1 ms 到 2 s，其中第一個周期等待 1 ms，此后每個 TR 增加 8 ms；90° 射頻脈沖后立即采樣獲得 Si （t），采樣過程持續 10.24 ms。
3 實驗結果
實驗首先測量了 X 方向的梯度線圈在 X 、 Y 、 Z 各方向產生的 EC 感應磁場G（t），每次測量需做兩次，以完成梯度波形有預加重、無預加重的對比。其中時間軸在階躍電流信號關閉并開始接收信號時重置為 0。每個 TR 獲得一段數據，每段數據包含 128 個數據（如前所述采樣時間為 10.24 ms，每 0.08 ms 采一個數據，共采得 128 個數據），經平臺 Matlab 后處理時將所有段數據（共 256 段）拼接起來，如圖 2 所示，即灰色部分，獲得約 2 s 的 EC 磁場衰減情況。由于其中包含大量背景噪聲，還需用最小二乘法進行曲線擬合，如圖中黑色曲線所示。由于篇幅限制，圖 2 僅展示了 X 方向的梯度線圈在 X 、 Z 方向感應出的 EC 感應磁場。
在進行渦流補償（eddy current compensation，ECC）以前， X 方向梯度線圈感應出 X 方向的 EC 的采樣數據可根據式（12）算出，并拼接得到。采樣數據的曲線擬合采用 e 指數模型，若擬合為三項指數模型，則可以寫為
| ${G_{x{\rm{2}}x}}(t) = {{\rm{a}}_1}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_1}}}}} + {{\rm{a}}_2}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_2}}}}} + {{\rm{a}}_3}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_3}}}}}$ |
其中 a1=—0.023 675，τ1=1.643 1×102；a2=—1.249 9×10–4，τ2=2.649 4×102；a3=—3.906 2×10–6，τ3=2.690 3×102。采樣數據和擬合曲線如圖 2 左上圖所示。
圖2
ECC 前后感應磁場的比較
Figure2.
Induced magnetic field compensation before and after ECC
而圖中展示的 X 方向梯度線圈在 Z 方向的 EC 感應磁場，可以寫為
| ${G_{x{\rm{2}}z}}(t) = {{\rm{a}}_1}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_1}}}}} + {{\rm{a}}_2}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_2}}}}} + {{\rm{a}}_3}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_3}}}}}$ |
其中 a1=6.210 9×10–4，τ1=2.632 7×102；a2=—3.906 2×10–6，τ2=2.690 3×102；a3=—8.007 8×10–5，τ3=2.740 9×102。由于此處a3>a2，所以此處的二指數項已經足以描述當前精度。
同理， X 方向梯度線圈在 Y 方向的 EC 感應磁場也可以測出：
| ${G_{x{\rm{2}}y}}(t) = {{\rm{a}}_1}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_1}}}}} + {{\rm{a}}_2}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_2}}}}}$ |
其中 a1=1.909 2×10–4，τ1=2.707 2×102; a2=3.906 2×10–6，τ2=2.715 2×102.
圖中還展示了使用預加重脈沖波形的結果，即進行 ECC 后的 EC 效應，在相同條件下，用專用工裝上的水模測出 EC 效應明顯下降，在t=0 時，EC 的感應磁場下降至 10–5 Gs/cm 數量級。
以 X 方向校正波形為例，施加梯度 12.5 mT/m 換算后得 1.25 Gs/cm，從圖 2 可看出 EC 的感應磁場在 10–3 Gs/cm 數量級，將E=1.25 Gs/cm，及 a1～a3，τ1～τ3 帶入式（16）中，可得出當輸入為方波 r2T 時，輸出響應為G（t）。如果要求輸出為方波，則根據式（15），可求出校正后的輸入脈沖的傅里葉變換X（jω），再進行反變換，即可得到輸入脈沖x（t）。如圖 3 所示。其余校正波形可同理得到。
圖3
感應脈沖和校正脈沖
Figure3.
Induced magnetic field and correction pulse
如圖 4 所示，比較了在 ECC 前后用 19 cm 直徑球形水模成像的模型。補償前的成像圖在邊緣處表現出變形和偽影，在圖像中測出形變最大處直徑僅為 15.83 cm，而補償后成像直徑在 18.7～19.1 cm 范圍內，各方向直徑基本相等。除此之外，還可以采用迭代進一步校正輸入脈沖，即使用校正脈沖后再次測試序列，測試校正后 EC 的感應磁場，以達到進一步降低 EC 的作用。肉眼觀察迭代后的水模圖像和補償前區別不大。人體試驗如圖 5 所示，在醫用 1.5 T、710 mm 孔徑的超導 MRI 上進行的試驗：受試者為女性志愿者，32 歲，試驗采用頭線圈采樣數據，彌散加權成像（diffusion weighted ima-ging，DWI）掃描序列。掃描參數為 TR：5 500 ms；回波時間：89.1 ms；翻轉角 1：90°；翻轉角 2：160°；相位視野（field of view，FOV）：240 mm；讀出 FOV：240 mm；層厚 6 mm；總掃描層數：18 層；DWI 方向：3；最大擴散敏感系數 b：1 000 s/mm2。如圖 5 所示，比較了 1 次 ECC 和 3 次 ECC 后 DWI 成像的區別，圖中節選兩個不同的斷層圖片，其中左圖中箭頭標注的是經 1 次 ECC 補償圖像，邊緣相對模糊，而右圖經 3 次 ECC 補償后圖像邊緣更清晰。這說明測試結果可以準確地指導數字濾波器實現脈沖波形的預加重，而成像結果也驗證了 EC 測試過程和測試結果的正確性。
圖4
ECC 前后球形水模（
=19 cm）成像
Figure4.
Spherical phantom（
=19 cm）images before and after ECC
圖5
1 次 ECC 和 3 次 ECC 后 DWI 成像 箭頭標記出模糊的邊緣
Figure5.
Comparison of DWI images ECC between one ECC and three ECCs The differences are marked by arrows.
4 結論
EC 可能引起 MRI 圖像的幾何畸變，也可能產生偽影，降低成像速度。為了降低 EC 引起的不均勻場，理論上可采取預加重梯度線圈脈沖波形的方式，將脈沖波形由矩形脈沖更改為適合當前磁共振設備 EC 的脈沖波形，為此，必須首先測試設備產生的不均勻磁場強度，再根據測試結果，推導預加重梯度線圈脈沖波形，設計數字濾波器實現補償。由于脈沖波形是根據測試結果設計的，這說明對每一臺磁共振設備，即使型號相同，由于所用零件批次可能不同，裝機環境不同，人為因素等原因，也應在裝機后專門進行 EC 校正。
測試中，如果假設設備產生的 EC 在矩形脈沖的上升階段和下降階段互為逆過程，并假設脈沖持續時間足夠長，以至于上升階段和下降階段不會互相影響，則僅需測試上升和下降其中之一即可。所以在實驗過程中應盡量使上面兩個假設得到滿足。本文的測試中，測試了梯度線圈在加載足夠持續時間的階躍脈沖時感應產生的 EC，即僅測試矩形脈沖的上升階段的系統響應，下降過程可相應導出，從而得到整個矩形脈沖的系統響應。
測試需用專有工裝，共有 9 個參數需要測試。即當 X 、 Y 、 Z 方向梯度線圈分別加載階躍脈沖時，對 X 、 Y 、 Z 方向 EC 的感應磁場，用這 9 個參數指導設計校正脈沖。
另外，多次迭代可進一步降低 EC，因此本文通過對比實驗數據（如圖 2 所示）發現，1 次 ECC 后的 EC 效應即可得到明顯壓制，EC 的感應磁場下降了 2 個數量級，至 10–5 Gs/cm，而多次重復可使 EC 效應進一步降低，成像質量得到改善。MRI 的測試序列繁多，而本文主要用 DWI 序列完成人體試驗，其它需用到梯度短脈沖的序列還需要進一步實驗驗證。
引言
隨著磁共振成像（magnetic resonance imaging，MRI）設備在國內醫院的逐漸普及，通過 MRI 在盡量短的時間內獲得高分辨率、高信噪比、低偽影的高質量圖像成為未來醫學診斷領域的發展方向，其中脈沖梯度磁場技術是空間編碼、彌散加權成像的關鍵技術，同時也是制約縮短成像時間的主要問題之一。脈沖梯度磁場由置于 MRI 設備中的梯度線圈產生，應用在幾乎所有的成像序列中[1]。由于在產生脈沖梯度的同時也會感應出渦流（eddy cur-rents，EC），從而會導致回波信號的缺失和變形，引起偽影。具體來說，當梯度線圈由一個梯形波驅動時，會在其反方向感應出 EC，也就是在梯形的上升沿或下降沿產生相反的效應[2]，使上升或下降時間持續幾百微秒。雖然理論上可采用延長時間并等待 EC 自行消失的辦法以減輕偽影，但由于成像過程中需要幾百次的梯度線圈開關，從而使得實際上等待 EC 自行消失的辦法操作性較低，而某些需要短脈沖梯度磁場的序列更是無法實現。
由于 EC 是由 MRI 設備中多個零部件共同引起的，所以找到 EC 產生的具體部分并不具有可操作性，即從理論上找出源頭、建模、解決問題的方式不合適。因此采用直接定量測定 EC 的大小，并找到具有操作性的辦法從工程上來說比較實際。為此，國內外提出了各種測量和降低 EC 的辦法，例如改變梯度序列[3]，改變射頻脈沖傳輸方式[4]，或者預加重波形補償等辦法[5-7]。目前尚未見關于提取實驗數據和關于如何實際應用這些數據獲得波形的報道。基于已有研究成果，本文應用先驅的實驗方法在超導 MRI 設備的梯度線圈上進行實驗，獲得實驗數據；進一步推導了 EC 產生的系統函數，計算出預加重波形的形狀；并提出利用迭代方法，可進一步減輕偽影。在本文的具體測試過程中，首先根據數學模型編寫測量序列，然后根據在 MRI 設備上的測試實驗完成數據采樣，再將采樣數據在 Matlab 平臺中處理，計算出關鍵參數和校正波形，最后將校正波形在譜儀中實現，并驗證實驗結果。
本文運用預加重波形補償的辦法，詳細描述了實驗過程中如何在設備上提取實驗數據和如何應用這些數據獲得預加重波形以有效減小渦流，并在原有基礎上提出預加重波形的推導過程和改進辦法，為國內磁共振設備研究企業和機構提供一點參考。
1 理論基礎
1.1 不均勻磁場的求解
測試從推導求解不均勻磁場理論分布的表達式開始，尋找緩解 EC 造成的偽影的方法。考慮一個由梯度線圈產生的梯度磁場，EC 將在其上升或下降沿產生，因此采用定量 EC 大小的方法來處理。我們通過在梯度線圈上施加一個較長持續時間的階躍波形，當梯度關閉的時候，再施加一個 90° 射頻脈沖用以在水模中產生自由衰減信號（free induction decay signal，FID）[6]。然后立即用接收線圈接收，經模擬/數字（analogue to digital，A/D）轉換獲得一維采樣數據，該 FID 只在 EC 背景下采得，因此可以反映出 EC 的大小。若假設上升沿和下降沿的 EC 相等，則上升沿情況可同理推出。
在下降沿，Bloch 方程描述了時變函數磁矢量 M 和加載的磁場強度 B 的關系[8]：
| $\frac{{{\rm{d}}{{M}}}}{{{\rm{d}}t}} \!=\! {\rm{\gamma }}{{M}} \! \times \! {{B}} \!-\! \frac{{{M_x}{{X}} \! + \! {M_y}{{Y}}}}{{{{\rm{T}}_2}}} \!-\! \frac{{({M_z} \!-\! \left| {{{{M}}_0}} \right|){{Z}}}}{{{{\rm{T}}_1}}}$ |
其中 M 0 是僅有主磁場 B 0 時沿 Z 軸的縱向磁矢量，理論上只有 Z 方向的分量，T1 和 T2 分別表示橫向和縱向弛豫時間。
假設旋轉坐標頻率等于拉莫爾進動頻率，則式（1）右邊第一項變為
| $\left\{ \begin{array}{l}\frac{{{\rm{d}}{M_z}}}{{{\rm{d}}t}} = - \frac{{({M_z} - \left| {{{{M}}_0}} \right|)}}{{{{\rm{T}}_1}}}\\\frac{{{\rm{d}}{{{M}}_{xy}}}}{{{\rm{d}}t}} = - \frac{{{{{M}}_{xy}}}}{{{{\rm{T}}_2}}}\end{array} \right.$ |
由 90° 射頻脈沖產生的橫向磁化矢量的 FID 信號可由式（2）解出
| $\left\{ \begin{array}{l}{{{M}}_{xy}} = {M_x} + {{i}}{M_y} = \left| {{{{M}}_0}} \right|{{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_{\rm{2}}}}}}}{{\rm{e}}^{ - {{i}} { \text{ω}} {}_0t}}\\{M_z} = \left| {{{{M}}_0}} \right|(1 - {{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_{\rm{1}}}}}}})\end{array} \right.$ |
假設 B （ r ）是實驗室坐標下 r 處 EC 的感應磁場，它感應出的接收線圈中的電壓V（t）為
| $\begin{aligned}& V(t) = \! - \! \frac{\partial }{{\partial t}}\int_{{\rm{object}}}^{} \! {{{B}}({{r}}) \! \cdot \! {{M}}({{r}},t){\rm{d}}r} \\ &= - \! \int_{{\rm{object}}} \!\!\! {[{B_{r,x}}({{r}}) \!\! \frac{{\partial {M_x}({{r}},t)}}{{\partial t}} \!\! + \! {B_{r,y}}({{r}}) \! \frac{{\partial {M_y}({{r}},t)}}{{\partial t}}]} {\rm{d}}r\end{aligned}$ |
其中，Br，x（ r ）、Br，y（ r ）是 B （ r ）在 r 處的 X 、 Y 方向的分量，Mx （ r ，t）、My （ r ，t）是 M （ r ，t）在 r 處的 X 、 Y 方向的分量。由于Mz （ r ，t）相對于 |Mxy （ r ，t）| 變化緩慢，所以Br，z（ r ）Mz （ r ，t）被忽略。將式（3）代入式（4）中，簡化后得
| $\begin{aligned}{V_{\sin }}(t) = & - \int_{{\rm{object}}} { {{ω}} ({{r}}){B_{r,xy}}({{r}})\left| {{{{M}}_0}} \right|} \cdot \\ &{{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_2}}}}}\sin [ {{ω}} ({{r}}) + {{ \text{φ}}_0}]{\rm{d}}{{r}}\end{aligned}$ |
或
| $\begin{aligned}{V_{\cos }}(t) = & \int_{{\rm{object}}} { {{ω}} ({{r}}){B_{r,xy}}({{r}})\left| {{{M}}_0} \right|} \cdot \\ &{{\rm e}^{ - \frac{t}{{{T_2}}}}}\cos [ {{ω}} ({{r}}) + {{ \text{φ}}_0}]{\rm{d}}{{r}}\end{aligned}$ |
其中用ω（ r ）代替了式（3）中的ω0，表示此處ω 是位置的函數。最后用相位檢測器（phase-sensitive detection，PSD）獲得信號 S （t），其中實部SR （t） 等于Vsin（t），虛部SI （t） 等于Vcos（t）：
| ${{{S}}(t) = {S_R}(t) + i{S_I}(t)$ |
如果不均勻磁場 ΔB（τ）是和時間相關的函數，則式（5）、式（6）中的相位可以寫為
| $\phi = {\text{ω}} ({{{r}}}) + {{φ}_0} = {\rm{\gamma }}\int_0^t {\Delta B({ \text{τ}} ){\rm{d}}{\text{τ}} } + \Delta {\text{ω}} t$ |
將式（5）和式（6）代入式（7），可以得到
| ${{S}}(\!t\!) \!=\!\!\! \int_{{\rm{object}}}^{} \!\!\!\!{ { \text{ω}} (\!{{r}}\!){{{B}}_{r,xy}}(\!{{r}}\!) \! \left| {{{{M}}_0}} \right| \! {{\rm{e}}^{ - \frac{t}{{{{\rm{T}}_{\rm{2}}}}}}}{{\rm{e}}^{ - {{i}}[{\rm{\gamma }}\int_0^t {\Delta B(\! { \text{τ}} \!){\rm d}{ \text{τ}} + \Delta { \text{ω}} t} ]}}{\rm{d}}{{r}}} $ |
式中不均勻磁場 ΔB（τ）是需要求解的量，包含在信號 S （t）的相位信息中，所以在后續的計算中僅需要關注相位信息，即式（5）、式（6），又式（5）、式（6）中 φ0 是初相位，為常數，所以相位 的變化量僅由不均勻磁場 ΔB（τ）積分造成，考察 S （t）的相位變化量即可求得 ΔB（τ）， S （t）可根據式（7）用 PSD 測量獲得。出于習慣，將變量τ 用t 代換，假設 ΔB（t）由以下兩項構成
| $\Delta B(t) = {B_0}(t) + G(t)x$ |
其中B0（t）是偽空間場，而G（t）是由 EC 引起的偽梯度場，x 是在 X 軸上的位置。
由于同一磁體系統的B0（t）、G（t）僅為時間函數，在 X 軸不同位置x1、x2 處上放置水模，即可獲得不同的 ΔB1（t）和 ΔB2（t），其表達式為：
| $\left\{ \begin{array}{l}\Delta {B_1}(t) = {B_0}(t) + G(t){x_1}\\\Delta {B_2}(t) = {B_0}(t) + G(t){x_2}\end{array} \right.$ |
式（11）兩式相減，即可計算G（t）的表達式。兩式分別乘以x2、x1，再相減，可計算B0（t）的表達式：
| $\left\{ \begin{array}{l}G(t) = \frac{{\Delta {B_1}(t) - \Delta {B_2}(t)}}{{{x_1} - {x_2}}}\\{B_0}(t) = \frac{{{x_1}\Delta {B_2}(t) - {x_2}\Delta {B_1}(t)}}{{{x_1} - {x_2}}}\end{array} \right.$ |
為了求解B0（t）、G（t），x1、x2 可直接測量水模位置得出；ΔB1（t）和 ΔB2（t）可在得到采樣信號 S （t）后，根據式（7），信號相位 對時間的微分得出。
1.2 測試和擬合
根據以上理論推導，可得為解出G（t）和B0（t）的關鍵是要得到信號相位 ，而 是通過式（7）反正切計算得出： ，為此測試序列至少應包含：① X 方向梯度線圈加載階躍電流信號i（t）=u（t）；② 通過體線圈輸出 90° 射頻脈沖；③ 在回波時間（echo time，TE）后，從接收線圈接收信號，獲得式（7）中的一組一維數據 S （t）；但是，由于 EC 的持續時間為 105 s 量級，即需要考察的 S （t）應持續至少 1 s，而單次 FID 由于衰減較快，采樣能獲得的有效信號時間長度僅為 10–2 s 量級，所以必須采用多次采集數據，最后將各段數據拼接起來的方式獲得整個 S （t）[6]，所以測量的第 ③ 步應該改為：多次改變 TE 的值，獲得多組數據 Si （t），拼接獲得 S （t）。
獲得 S （t）后，即可根據式（12）計算出B0（t）、G（t），至此，B0（t）、G（t）的原始數據已全部獲得。但由于其中包含大量隨機噪聲，還需用最小二乘法進行曲線擬合。根據電感—電阻串聯（induc-tance-resistance，L-R）模型擬合G（t）-t 曲線。由 L-R 模型，G（t）-t 曲線可以描述為多指數項的和[8]，如式（13）所示：
| $G(t) = \sum\limits_{k = 1}^n {{{\rm{a}}_k}{{\rm{e}}^{ - \frac{t}{{{{\rm{{ \text{τ}} }}_k}}}}}} $ |
其中 ak 是常系數，τk 是時間常數，n 是指數項的個數。最小二乘法可用作擬合算法[9]。
同理可推出B0（t）-t 曲線。
應當指出的是，此處得到的G（t）和B0（t）都是 X 方向的梯度Gx（t）用 X 方向放置的水模測出的 X 方向 EC 的感應磁場，因此應將其寫為Gx2x（t）和B0x2x（t）。實際上為了更準確的測量，還應測量交叉項，即Gx2y（t）和Bx2y（t）、Gx2z（t）和Bx2z（t）、Gy2z（t）和By2z（t）等。
1.3 波形校正
將形成 EC 效應的整體考慮為一個系統H，如圖 1 所示，輸入為x（t），輸出為y（t）。
圖1
EC 效應的系統
Figure1.
System of EC
x（t）經傅里葉變換為X（jω），y（t）經傅里葉變換為Y（jω），則如圖 1 所示系統函數可表示為：
| $Y({\rm{j}} { \text{ω}} ) = X({\rm{j}} { \text{ω}} ) \cdot H({\rm{j}} { \text{ω}} )$ |
根據要求，波形校正后輸出信號y（t）是理想的方波信號 r（t），傅里葉變換為 R（jω），則如果已知系統函數H（jω），便可根據式（14）得到輸入函數x（t）的傅里葉變換X（jω），再進行反變換，即可得到輸入函數x（t），實現波形校正。
| $X({\rm{j}} { \text{ω}} ) = \frac{1}{{H({\rm{j}} { \text{ω}} )}} \cdot {\rm{R}}({\rm{j}} { \text{ω}} )$ |
問題變為求解系統函數H（jω）。為此，根據小節 1.2 中介紹的測試和擬合方法，可得一對已知條件：輸入x（t）為階躍信號 u（t），輸出y（t）為式（13）計算后得到的擬合曲線G（t），由式（14）可知，將x（t）、y（t）分別進行傅里葉變換，可最后求得系統函數。
同理，根據階躍函數和矩形脈沖函數的關系，也可以推導另一對已知條件，同樣可求出系統函數：如果輸入x（t）為矩形波 r2T（t），r2T是高度為 1 的矩形函數，寬度為（–T，T），則在假設矩形波上升沿和下降沿產生的 EC 大小相等，方向相反的情況下，由 L-R 模型，矩形波 r2T（t）引起的系統響應G（t）可以表示為
| $\begin{aligned}G(t) = & (E - \sum\limits_{k = 1}^n {{{\rm{a}}_k}{{\rm{e}}^{ - \frac{{t + T}}{{{{ \text{τ}}_k}}}}}} ){{\rm{r}}_{2T}} + \\ &E \sum\limits_{k = 1}^n {{{\rm{a}}_k}{{\rm{e}}^{ - \frac{{t - T}}{{{{ \text{τ}}_k}}}}}} {\rm{u}}(t - T)\end{aligned}$ |
其中E 表示輸入方波x（t）的最大值；u（t）為高度為 1 的階躍函數。
2 實驗
本文實驗在一臺醫用 1.5 T、710 mm 孔徑的超導 MRI（四川成都，奧泰醫療）上進行。施加的梯度為 12.5 mT/m。測試工裝在 X 、 Y 、 Z 三個軸的正負對稱六個位置分別放置一個直徑為 2 mm 的球形水模，并將其置于磁體孔徑的中心位置。每個水模距離中心位置的距離為 10 cm，且包裹接收線圈。
脈沖序列設計為梯度線圈加載的階躍電流，信號持續時間為 1 s，每一個重復周期（repeat time，TR）顛倒梯度方向，即改變階躍電流方向；階躍信號停止后經等待時間（time delay，TD）施加 90° 射頻脈沖，TD 的變化范圍為 1 ms 到 2 s，其中第一個周期等待 1 ms，此后每個 TR 增加 8 ms；90° 射頻脈沖后立即采樣獲得 Si （t），采樣過程持續 10.24 ms。
3 實驗結果
實驗首先測量了 X 方向的梯度線圈在 X 、 Y 、 Z 各方向產生的 EC 感應磁場G（t），每次測量需做兩次，以完成梯度波形有預加重、無預加重的對比。其中時間軸在階躍電流信號關閉并開始接收信號時重置為 0。每個 TR 獲得一段數據，每段數據包含 128 個數據（如前所述采樣時間為 10.24 ms，每 0.08 ms 采一個數據，共采得 128 個數據），經平臺 Matlab 后處理時將所有段數據（共 256 段）拼接起來，如圖 2 所示，即灰色部分，獲得約 2 s 的 EC 磁場衰減情況。由于其中包含大量背景噪聲，還需用最小二乘法進行曲線擬合，如圖中黑色曲線所示。由于篇幅限制，圖 2 僅展示了 X 方向的梯度線圈在 X 、 Z 方向感應出的 EC 感應磁場。
在進行渦流補償（eddy current compensation，ECC）以前， X 方向梯度線圈感應出 X 方向的 EC 的采樣數據可根據式（12）算出，并拼接得到。采樣數據的曲線擬合采用 e 指數模型，若擬合為三項指數模型，則可以寫為
| ${G_{x{\rm{2}}x}}(t) = {{\rm{a}}_1}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_1}}}}} + {{\rm{a}}_2}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_2}}}}} + {{\rm{a}}_3}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_3}}}}}$ |
其中 a1=—0.023 675，τ1=1.643 1×102；a2=—1.249 9×10–4，τ2=2.649 4×102；a3=—3.906 2×10–6，τ3=2.690 3×102。采樣數據和擬合曲線如圖 2 左上圖所示。
圖2
ECC 前后感應磁場的比較
Figure2.
Induced magnetic field compensation before and after ECC
而圖中展示的 X 方向梯度線圈在 Z 方向的 EC 感應磁場，可以寫為
| ${G_{x{\rm{2}}z}}(t) = {{\rm{a}}_1}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_1}}}}} + {{\rm{a}}_2}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_2}}}}} + {{\rm{a}}_3}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_3}}}}}$ |
其中 a1=6.210 9×10–4，τ1=2.632 7×102；a2=—3.906 2×10–6，τ2=2.690 3×102；a3=—8.007 8×10–5，τ3=2.740 9×102。由于此處a3>a2，所以此處的二指數項已經足以描述當前精度。
同理， X 方向梯度線圈在 Y 方向的 EC 感應磁場也可以測出：
| ${G_{x{\rm{2}}y}}(t) = {{\rm{a}}_1}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_1}}}}} + {{\rm{a}}_2}{{\rm{e}}^{ - \frac{t}{{{{ \text{τ}}_2}}}}}$ |
其中 a1=1.909 2×10–4，τ1=2.707 2×102; a2=3.906 2×10–6，τ2=2.715 2×102.
圖中還展示了使用預加重脈沖波形的結果，即進行 ECC 后的 EC 效應，在相同條件下，用專用工裝上的水模測出 EC 效應明顯下降，在t=0 時，EC 的感應磁場下降至 10–5 Gs/cm 數量級。
以 X 方向校正波形為例，施加梯度 12.5 mT/m 換算后得 1.25 Gs/cm，從圖 2 可看出 EC 的感應磁場在 10–3 Gs/cm 數量級，將E=1.25 Gs/cm，及 a1～a3，τ1～τ3 帶入式（16）中，可得出當輸入為方波 r2T 時，輸出響應為G（t）。如果要求輸出為方波，則根據式（15），可求出校正后的輸入脈沖的傅里葉變換X（jω），再進行反變換，即可得到輸入脈沖x（t）。如圖 3 所示。其余校正波形可同理得到。
圖3
感應脈沖和校正脈沖
Figure3.
Induced magnetic field and correction pulse
如圖 4 所示，比較了在 ECC 前后用 19 cm 直徑球形水模成像的模型。補償前的成像圖在邊緣處表現出變形和偽影，在圖像中測出形變最大處直徑僅為 15.83 cm，而補償后成像直徑在 18.7～19.1 cm 范圍內，各方向直徑基本相等。除此之外，還可以采用迭代進一步校正輸入脈沖，即使用校正脈沖后再次測試序列，測試校正后 EC 的感應磁場，以達到進一步降低 EC 的作用。肉眼觀察迭代后的水模圖像和補償前區別不大。人體試驗如圖 5 所示，在醫用 1.5 T、710 mm 孔徑的超導 MRI 上進行的試驗：受試者為女性志愿者，32 歲，試驗采用頭線圈采樣數據，彌散加權成像（diffusion weighted ima-ging，DWI）掃描序列。掃描參數為 TR：5 500 ms；回波時間：89.1 ms；翻轉角 1：90°；翻轉角 2：160°；相位視野（field of view，FOV）：240 mm；讀出 FOV：240 mm；層厚 6 mm；總掃描層數：18 層；DWI 方向：3；最大擴散敏感系數 b：1 000 s/mm2。如圖 5 所示，比較了 1 次 ECC 和 3 次 ECC 后 DWI 成像的區別，圖中節選兩個不同的斷層圖片，其中左圖中箭頭標注的是經 1 次 ECC 補償圖像，邊緣相對模糊，而右圖經 3 次 ECC 補償后圖像邊緣更清晰。這說明測試結果可以準確地指導數字濾波器實現脈沖波形的預加重，而成像結果也驗證了 EC 測試過程和測試結果的正確性。
圖4
ECC 前后球形水模（
=19 cm）成像
Figure4.
Spherical phantom（
=19 cm）images before and after ECC
圖5
1 次 ECC 和 3 次 ECC 后 DWI 成像 箭頭標記出模糊的邊緣
Figure5.
Comparison of DWI images ECC between one ECC and three ECCs The differences are marked by arrows.
4 結論
EC 可能引起 MRI 圖像的幾何畸變，也可能產生偽影，降低成像速度。為了降低 EC 引起的不均勻場，理論上可采取預加重梯度線圈脈沖波形的方式，將脈沖波形由矩形脈沖更改為適合當前磁共振設備 EC 的脈沖波形，為此，必須首先測試設備產生的不均勻磁場強度，再根據測試結果，推導預加重梯度線圈脈沖波形，設計數字濾波器實現補償。由于脈沖波形是根據測試結果設計的，這說明對每一臺磁共振設備，即使型號相同，由于所用零件批次可能不同，裝機環境不同，人為因素等原因，也應在裝機后專門進行 EC 校正。
測試中，如果假設設備產生的 EC 在矩形脈沖的上升階段和下降階段互為逆過程，并假設脈沖持續時間足夠長，以至于上升階段和下降階段不會互相影響，則僅需測試上升和下降其中之一即可。所以在實驗過程中應盡量使上面兩個假設得到滿足。本文的測試中，測試了梯度線圈在加載足夠持續時間的階躍脈沖時感應產生的 EC，即僅測試矩形脈沖的上升階段的系統響應，下降過程可相應導出，從而得到整個矩形脈沖的系統響應。
測試需用專有工裝，共有 9 個參數需要測試。即當 X 、 Y 、 Z 方向梯度線圈分別加載階躍脈沖時，對 X 、 Y 、 Z 方向 EC 的感應磁場，用這 9 個參數指導設計校正脈沖。
另外，多次迭代可進一步降低 EC，因此本文通過對比實驗數據（如圖 2 所示）發現，1 次 ECC 后的 EC 效應即可得到明顯壓制，EC 的感應磁場下降了 2 個數量級，至 10–5 Gs/cm，而多次重復可使 EC 效應進一步降低，成像質量得到改善。MRI 的測試序列繁多，而本文主要用 DWI 序列完成人體試驗，其它需用到梯度短脈沖的序列還需要進一步實驗驗證。
