(57)【特許請求の範囲】
【請求項１】
静磁場発生手段、傾斜磁場発生手段、高周波磁場発生手段および高周波磁場検出手段を備えた撮像手段と、撮像シーケンスに基き前記撮像手段の動作を制御する制御手段と、前記高周波磁場検出手段が検出した核磁気共鳴信号を用いて、画像再構成を含む演算を行う演算手段と、を備えた磁気共鳴イメージング装置であって、
前記演算手段は、前記傾斜磁場発生手段に入力される入力波形と当該入力波形によって前記傾斜磁場発生手段が発生した出力波形とを用いて、前記出力波形に影響を与える複数の要素の応答関数の合算である応答関数を算出する応答関数算出手段と、前記応答関数計算手段が計算した応答関数を用いて、前記撮像シーケンスに設定された傾斜磁場パルスの入力波形から出力波形を算出する傾斜磁場算出手段と、前記傾斜磁場算出手段が算出した出力波形の計算値を用いて、前記撮像シーケンスに設定された高周波磁場パルスの修正および/または前記撮像シーケンスで得られた核磁気共鳴信号のk空間座標の修正を行う補償手段とを備えたことを特徴とする磁気共鳴イメージング装置。
【請求項２】
請求項1に記載の磁気共鳴イメージング装置であって、
前記応答関数算出手段は、前記出力波形を変化させる複数の要素毎に応答関数を算出し、前記傾斜磁場算出手段は、前記複数の要素の応答関数の合算である応答関数を用いて出力波形の算出を行うことを特徴とする磁気共鳴イメージング装置。
【請求項３】
請求項2に記載の磁気共鳴イメージング装置であって、
前記応答関数算出手段は、前記要素毎に応答関数モデルを設定するモデル設定手段と、前記モデル設定手段が設定した応答関数モデルによって算出した出力波形の計算値と、出力波形の実測値との差が最小となるように応答関数モデルを適正化するモデル適正化手段とを備え、前記モデル適正化手段が適正化した応答関数モデルを各要素の応答関数とすることを特徴とする磁気共鳴イメージング装置。
【請求項４】
請求項1に記載の磁気共鳴イメージング装置であって、
前記応答関数算出手段は、前記入力波形として複数の入力値を用い、複数の入力値とそれに対応する実測値との組合せを用いて、入力値と出力波形の実測値との関係から、直接、応答関数を算出することを特徴とする磁気共鳴イメージング装置。
【請求項５】
請求項4に記載の磁気共鳴イメージング装置であって、
前記応答関数算出手段は、算出後の応答関数にフィルタをかけるフィルタ手段を備えることを特徴とする磁気共鳴イメージング装置。
【請求項６】
請求項5に記載の磁気共鳴イメージング装置であって、
前記フィルタ手段は、ウィナーフィルタ、メディアンフィルタを含むことを特徴とする磁気共鳴イメージング装置。
【請求項７】
請求項1に記載の磁気共鳴イメージング装置であって、
前記傾斜磁場発生手段は、発生した傾斜磁場の出力波形に基き入力波形を制御する制御回路を備え、
前記出力波形を変化させる複数の要素は、前記制御回路による制御に伴う前記出力波形の歪みを含むことを特徴とする磁気共鳴イメージング装置。
【請求項８】
請求項1に記載の磁気共鳴イメージング装置であって、
前記応答関数算出手段は、前記静磁場発生手段が発生する静磁場の中心に対する位置が異なる複数の傾斜磁場について、それぞれ、応答関数を算出することを特徴とする磁気共鳴イメージング装置。
【請求項９】
請求項1に記載の磁気共鳴イメージング装置であって、
前記傾斜磁場算出手段が算出する傾斜磁場パルスは、前記撮像シーケンスにおけるスライス選択傾斜磁場であり、前記補償手段は、前記スライス選択傾斜磁場の計算値に基き、前記撮像シーケンスに設定された高周波パルスを補正する高周波パルス補正手段を備えたことを特徴とする磁気共鳴イメージング装置。
【請求項１０】
請求項9に記載の磁気共鳴イメージング装置であって、
前記高周波パルス補正手段は、補正前後の高周波パルスの照射時間の変化を算出する照射時間計算手段を備え、前記照射時間計算手段が算出した照射時間の変化が所定の増加率を超えるときに、前記傾斜磁場算出手段の入力値を変更し、前記傾斜磁場算出手段が前記変更された入力値から推定した出力波形の推定値を用いて、前記高周波パルスの補正を行うことを特徴とする磁気共鳴イメージング装置。
【請求項１１】
請求項1に記載の磁気共鳴イメージング装置であって、
前記傾斜磁場算出手段が算出する傾斜磁場パルスは、前記撮像シーケンスにおける読出し傾斜磁場であり、前記補償手段は、前記読出し傾斜磁場の計算値に基き、前記撮像シーケンスで得た核磁気共鳴信号が配置されるk空間座標を補正するk空間座標補正手段を備えたことを特徴とする磁気共鳴イメージング装置。
【請求項１２】
請求項11に記載の磁気共鳴イメージング装置であって、
前記k空間座標補正手段は、補正前後のk空間座標の大きさの変化を算出するk空間サイズ計算手段を備え、前記k空間サイズ計算手段が計算したk空間の大きさの変化が所定の縮小率を超えるときに、前記傾斜磁場算出手段の入力値を変更し、前記傾斜磁場算出手段が変更後の入力値から算出した出力波形の再計算値を用いて、前記k空間座標の補正を行うことを特徴とする磁気共鳴イメージング装置。
【請求項１３】
磁気共鳴イメージング装置の傾斜磁場発生手段が発生する出力傾斜磁場波形を推定する方法であって、
前記傾斜磁場発生手段の入力値およびそれによって前記傾斜磁場発生手段が発生した出力傾斜磁場波形の実測値を入力するステップと、
前記出力傾斜磁場波形に影響を与える複数の因子について、それぞれ、応答関数モデルを設定するステップと、
前記入力した傾斜磁場波形の入力値および実測値を用いて、各応答関数モデルを適正化するステップと、
所望の傾斜磁場の入力値に対し、前記適正化した応答関数モデルを用いて、出力傾斜磁場波形の推定値を算出するステップと、
を含む傾斜磁場波形推定方法。
【請求項１４】
請求項13に記載の傾斜磁場波形推定方法であって、
前記複数の因子は、傾斜磁場により発生する渦電流と、前記傾斜磁場発生手段に備えられた制御回路の出力とを含むことを特徴とする傾斜磁場波形推定方法。
【請求項１５】
請求項14に記載の傾斜磁場波形推定方法であって、
前記実測値は、前記傾斜磁場発生手段の制御回路を無効にした時の第1の実測値と、前記制御回路を有効にした第2の実測値を含み、
前記応答関数モデルを適正化するステップは、前記第１の実測値および第2の実測値を用いて、制御回路による出力傾斜磁場波形の応答関数を求めることを特徴とする傾斜磁場波形推定方法。
【請求項１６】
磁気共鳴イメージング装置における出力傾斜磁場波形を推定する方法であって、
傾斜磁場発生手段の少なくとも2種の入力値およびそれらによって前記傾斜磁場発生手段が発生した出力傾斜磁場波形の実測値を入力するステップと、
前記少なくとも2種の入力値およびそれに対応する実測値を用いて、前記出力傾斜磁場波形に影響を与える複数の因子の影響を含む前記傾斜磁場発生手段の応答関数を算出するステップと、
所望の傾斜磁場の入力値に対し、前記応答関数を用いて、出力傾斜磁場波形の推定値を算出するステップと、を含む傾斜磁場波形推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、被検体中の水素や燐等からの核磁気共鳴信号(以下、NMR信号と呼ぶ)を測定し、核の密度分布や緩和時間分布等を画像化する磁気共鳴イメージング(以下、MRIと呼ぶ)装置に関し、特に、傾斜磁場の歪みに起因する画質劣化を補償する技術に関する。
【背景技術】
【０００２】
MRI装置は、被検体、特に人体の組織を構成する原子核スピンが発生するNMR信号を計測し、その頭部、腹部、四肢等の形態や機能を2次元的に或いは3次元的に画像化する装置である。撮影においては、被検体を静磁場(分極磁場Ｂ0)内に配置した上で、特定の領域を選択励起するためにスライス選択傾斜磁場パルスと共に高周波磁場パルスを印加し、その後、位相エンコード傾斜磁場パルスや読出し傾斜磁場パルスを印加することにより、励起範囲内をエンコードし、位置情報を付与する。【０００３】
傾斜磁場パルスは、直交する三軸方向に線形の傾斜磁場を発生する複数のコイルに、傾斜磁場電源からパルス状に変化する電流(入力傾斜磁場波形という)を供給することにより、所望の傾斜磁場パルス(出力傾斜磁場波形という)を出力するようにしている。理想的には、入力傾斜磁場波形とそれによってMRI装置に発生する傾斜磁場波形は一致しなければならないが、種々の要因により誤差が生じる。【０００４】
誤差を生じる原因の一つは、傾斜磁場パルスによる急激な磁界の変化により、MRI装置を構成する磁性体部品あるいは電気回路に発生する渦電流による影響である。渦電流は、傾斜磁場パルスによる磁界の変化とは逆方向の磁界の変化を時間及び空間に依存して与え、出力傾斜磁場波形を歪ませる。【０００５】
また一般に、傾斜磁場コイルの出力傾斜磁場波形は、傾斜磁場コイルの特性(Ｑ値)に依存する応答遅れ等の歪みを含むが、これを補償するためにフィードバック制御などを行う回路が、傾斜磁場電源を含む傾斜磁場発生系には含まれている場合が多い。このような制御回路は、渦電流の補償の役割も果たすが、入力傾斜磁場波形と異なる出力傾斜磁場波形を生成する一因にもなり得る。【０００６】
これら渦電流及び傾斜磁場電源の制御回路の影響により生じる出力傾斜磁場波形の歪みは、入力傾斜磁場波形の形状に依存して変化し、その傾斜磁場の印加軸に応じて種々の問題を生じる。【０００７】
例えば、スライス選択傾斜磁場パルスの歪みは、励起プロファイル及び励起位置に誤差を生じさせる。特に傾斜磁場パルスの印加強度を変化させながら高周波磁場パルスを印加するVERSE(variable rate selective excitation)法などの場合には影響が大きい。【０００８】
読み出し傾斜磁場パルスの歪みは、画像に歪みやゴースト等のアーチファクトを生じさせる。特に、読み出し傾斜磁場パルスの印加極性を反転させながら計測を行うEPI(Echo Planar Imaging)計測や、k空間をらせん状に走査するスパイラル計測、k空間中心から走査を開始する超短TE計測などの場合に歪みの影響が顕著である。【０００９】
上記問題を解決するためには、入力傾斜磁場波形に依存して変化する出力傾斜磁場波形の歪みを補償する必要がある。既存の補償技術として、以下に示すような手法が提案されている。【００１０】
一つは、シムコイルや傾斜磁場コイルなど用いて渦電流を補償する磁場を印加するというものであり、傾斜磁場パルスを印加した後に発生する渦電流により生じる磁場を時間的・空間的に計測し、その計測結果から得られる渦電流の出力特性を打ち消す傾斜磁場パルスの波形を出力して出力傾斜磁場波形の歪みを補償する(特許文献1、特許文献2)。【００１１】
またスライス選択傾斜磁場パルスの歪みの補償に特化した手法として、出力傾斜磁場波形に応じて高周波磁場パルスの照射タイミングを変更する方法が提案されている(特許文献3)。この方法では、傾斜磁場パルスの歪みを測定し、パルス面積の重心から遅延時間を算出し、遅延時間をもとに高周波磁場パルスと傾斜磁場パルスのタイミングを変更することで、励起プロファイルの劣化を抑制する。また、固定の入力傾斜磁場波形を用いることで、入力傾斜磁場波形に応じた出力傾斜磁場波形の変化を排除する方法も提案されている(非特許文献1)。この方法では、高周波磁場パルスの励起幅を変更する場合は、入力傾斜磁場波形を変化させるのではなく、高周波磁場パルスの振幅を調整している。【００１２】
読出し傾斜磁場パルスの歪みの補償に特化した手法としては、システム応答をRLC回路でモデル化し、入力傾斜磁場波形に対する出力傾斜磁場波形を推定し、その結果を用いてNMR信号のk空間における座標位置を算出し、画像の歪みを補償する方法が提案されている(非特許文献2)。この方法では、モデル化した数式の係数は、画像の見た目を基準にして決定している。また上記手法のモデル化の代わりに、撮像終了後に傾斜磁場波形を実測した結果からNMR信号のk空間における座標位置を算出し、画像の歪みを補償する方法も提案されている(非特許文献3)。
【先行技術文献】
【特許文献】
【００１３】
【特許文献１】特開平10-272120号公報
【特許文献２】特開2010-42275号公報
【特許文献３】特開2006-149564号公報
【非特許文献】
【００１４】
【非特許文献１】A.M.Takahashi,Reduction of Slice Select Artifacts in Half Pulse Excitations used in Ultrashort TE(UTE)Imaging,ISMRM2010-4955
【非特許文献２】S.H.Cho et al,Compensation of eddy current by an R-L-C circuit model of the gradient system,Proc.Intl.Soc.Mag.Reson.Med.16:1156(2008)
【非特許文献３】Ethan K.Brodsky et al,Characterizing and Correcting Gradient Errors in Non-Cartesian Imaging:Are Gradient Errors Linear Time-Invariant(LTI)?,Magnetic Resonance in Medicine 62:1466-1476(2009)
【発明の概要】
【発明が解決しようとする課題】
【００１５】
しかしながら、上述した出力傾斜磁場波形の補償技術にはそれぞれ以下のような問題点がある。
まず、渦電流を補償する傾斜磁場パルスを出力する方法(特許文献1、特許文献2)では、時間的・空間的に高次の複雑な変化を示す渦電流について、渦電流により生じる磁場を完全に補正する渦電流補正コイルを実現することは困難かつ高コストであるという問題がある。また、傾斜磁場電源の制御回路の振る舞いについては考慮されておらず、出力傾斜磁場波形の歪みを完全に除去できない。【００１６】
次に、高周波磁場パルスの照射タイミングを出力傾斜磁場波形に応じて変更する方法(特許文献3)では、高周波磁場パルスの照射タイミングのみを変更しているため出力傾斜磁場波形の非線形な歪みについては補償できていない。入力傾斜磁場波形を固定する方法(非特許文献1)では、特定の入力傾斜磁場波形のみを使用して高周波磁場パルスの振幅を調整することは、SAR(Specific Absorption Rate：被検体内へのエネルギ蓄積量を示す基準値)の増加を招くおそれがある。【００１７】
さらにシステム応答をRLC回路でモデル化する方法(非特許文献2)では、システムの応答がRLC回路によるモデル式に当てはまらない場合、出力傾斜磁場波形の歪みを表現しきれず、画像の歪みやゴースト等のアーチファクトを残すことになる。これは、傾斜磁場電源によりフィードバック制御などの複雑な命令がなされている場合などに起こりえる。また、モデル式の各係数が画像の見た目を基準にして決定されるため、係数を決定する作業者に応じてモデル式の近似精度がばらつくという問題がある。【００１８】
撮像後に傾斜磁場波形を実測して再構成に用いる方法(非特許文献3)では、実測の傾斜磁場波形をk空間座標の算出に用いているため、ノイズによる揺らぎが重畳することになる。ノイズを低減するためにフィルタなどの処理を施した場合には、フィルタ処理による傾斜磁場波形の歪みが問題となる。【００１９】
そこで本発明は、入力傾斜磁場波形に応じて変化する出力傾斜磁場波形の歪みを精度よく算出し、上記算出した出力傾斜磁場波形を用いて各種補正を行うことを目的とする。
【課題を解決するための手段】
【００２０】
上記課題を解決する本発明は、入力傾斜磁場波形と当該入力傾斜磁場波形に対応する出力傾斜磁場波形を算出し、入力傾斜磁場波形と出力傾斜磁場波形とを用いて、該出力傾斜磁場波形に影響を与える複数の要素の応答関数の合算である応答関数を算出し、応答関数を用いて、撮像シーケンスに設定された傾斜磁場パルスの入力傾斜磁場波形から出力傾斜磁場波形を算出し、算出した出力傾斜磁場波形の計算値を用いて、各種補正を行う。補正の例としては、撮像シーケンスに設定された高周波磁場パルスの修正、及び/又は、撮像シーケンスで得られた核磁気共鳴信号のk空間座標の修正を行う。
【発明の効果】
【００２１】
本発明によれば、傾斜磁場の出力波形に影響を与える複数の要素、例えば、渦電流や傾斜磁場電源の制御回路のそれぞれの応答関数を加味した応答関数を算出し、それによって傾斜磁場の出力波形を算出するので、渦電流のみならず制御回路に起因する傾斜磁場の歪みも補償することができ、高コストな補償用傾斜磁場コイル等を用いることなく、従来技術では取り除くことができなかった画像の歪みやゴースト等のアーチファクトを抑制することができる。
【図面の簡単な説明】
【００２２】
【図１】MRI装置の全体構成を示す図。
【図２】ディジタル信号処理装置の機能ブロック図。
【図３】第一実施形態による全体的な処理手順を示すフローチャート。
【図４】出力傾斜磁場波形の空間別測定を説明する図。
【図５】渦電流により発生する磁場による出力傾斜磁場波形の歪みを説明する図。
【図６】第一実施形態における渦電流の応答関数を求める処理手順を示すフローチャート。
【図７】傾斜磁場電源の制御回路により出力傾斜磁場波形に加えられるオーバシュートもしくはアンダーシュートを説明する図。
【図８】第一実施形態における傾斜磁場電源の制御回路の応答関数を求める処理手順を示すフローチャート。
【図９】第一実施形態の応答関数を用いて推定した出力傾斜磁場波形と、入力磁場波形及び実際の出力傾斜磁場波形とを示すグラフ。
【図１０】渦電流のみを考慮した応答関数を用いて推定した出力傾斜磁場波形と、入力磁場波形及び実際の出力傾斜磁場波形とを示すグラフ。
【図１１】第一実施形態における出力傾斜磁場波形の推定処理のフローチャート。
【図１２】第一実施形態における入力傾斜磁場波形及び高周波磁場パルスの再計算処理手順を示すフローチャート。
【図１３】再計算前後の高周波磁場パルスをそれぞれ用いて得られた画像を示す図。
【図１４】第一実施形態における入力傾斜磁場波形及びk空間座標の再計算処理手順を示すフローチャート。
【図１５】k空間座標の再計算処理前後の再構成画像を示す図。
【図１６】第二実施形態における全体的な処理手順を示すフローチャート。
【図１７】第二実施形態における応答関数算出の処理手順を示すフローチャート。
【図１８】第三実施形態における全体的な処理手順を示すフローチャート。
【図１９】第三実施形態の応答関数の算出処理に用いる入力傾斜磁場波形の例を示す図。
【図２０】図19に示す台形波のラプラス変換波形を示す図。
【発明を実施するための形態】
【００２３】
以下、添付図面に従って本発明のMRI装置の好ましい実施形態について詳説する。なお、発明の実施形態を説明するための全図において、同一機能を有するものは同一符号を付け、その繰り返しの説明は省略する。
【００２４】
最初に、本発明に係るMRI装置の一例の全体概要を図1に基づいて説明する。図1は、本発明に係るMRI装置の一実施形態の全体構成を示すブロック図である。このMRI装置は、NMR現象を利用して被検体の断層画像を得るものであり、図1に示すように、静磁場発生系2と、傾斜磁場発生系3と、送信系5と、受信系6と、信号処理系7と、シーケンサ4とを備えて構成される。
【００２５】
静磁場発生系2は、図示していないが、永久磁石、常電導磁石あるいは超電導磁石の静磁場発生源を含み、被検体1が置かれる空間に均一な静磁場を発生させるものである。静磁場の方向によって、垂直磁場方式、水平磁場方式などがあり、垂直磁場方式であれば、被検体1の周りの空間にその体軸と直交する方向に、水平磁場方式であれば、体軸方向に均一な静磁場を発生させる。
【００２６】
傾斜磁場発生系3は、MRI装置の座標系(静止座標系)であるX、Y、Zの3軸方向に巻かれた傾斜磁場コイル9と、それぞれの傾斜磁場コイルを駆動する傾斜磁場電源10とから成り、後述のシ-ケンサ4からの命令に従ってそれぞれのコイルの傾斜磁場電源10を駆動することにより、X、Y、Zの3軸方向に傾斜磁場Gx、Gy、Gzを印加する。撮影時には、スライス面(撮影断面)に直交する方向にスライス方向傾斜磁場パルスを印加して被検体1に対するスライス面を設定し、そのスライス面に直交して且つ互いに直交する残りの2つの方向に位相エンコード方向傾斜磁場パルスと周波数エンコード方向傾斜磁場パルスを印加して、エコー信号にそれぞれの方向の位置情報をエンコードする。
【００２７】
シーケンサ4は、高周波磁場パルスと傾斜磁場パルスをある所定のパルスシーケンスで繰り返し印加する制御手段で、信号処理系7(ディジタル信号処理装置8)の制御で動作し、被検体1の断層画像のデータ収集に必要な種々の命令を送信系5、傾斜磁場発生系3、および受信系6に送る。パルスシーケンスは、撮像方法によって異なる種々のものがあり、それらは、予めプログラムとして磁気ディスク18に格納されている。
【００２８】
送信系5は、被検体1の生体組織を構成する原子の原子核スピンに核磁気共鳴を起こさせるために、被検体1に高周波磁場パルスを照射するもので、高周波発振器11と変調器12と高周波増幅器13と送信側の高周波コイル(送信コイル)14aとから成る。高周波発振器11から出力された高周波磁場パルスをシーケンサ4からの指令によるタイミングで変調器12により振幅変調し、この振幅変調された高周波磁場パルスを高周波増幅器13で増幅した後に被検体1に近接して配置された高周波コイル14aに供給することにより、高周波磁場パルスが被検体1に照射される。
【００２９】
受信系6は、被検体1の生体組織を構成する原子核スピンの核磁気共鳴により放出されるエコー信号(NMR信号)を検出するもので、受信側の高周波コイル(受信コイル)14bと信号増幅器15と直交位相検波器16と、A/D変換器17とから成る。送信側の高周波コイル14aから照射された電磁波によって誘起された被検体1の応答のNMR信号が被検体1に近接して配置された高周波コイル14bで検出され、信号増幅器15で増幅された後、シーケンサ4からの指令によるタイミングで直交位相検波器16により直交する二系統の信号に分割され、それぞれがA/D変換器17でディジタル量に変換されて、信号処理系7に送られる。
【００３０】
信号処理系7は、各種データ処理と処理結果の表示及び保存等を行うもので、ディジタル信号処理装置8と、磁気ディスク18、光ディスク19等の外部記憶装置と、ディスプレイ20と、ディジタル信号処理装置8による処理結果やディジタル信号処理装置8が使用するデータ等を保存するROM21、RAM22などの内部メモリを有し、受信系6からのデータがディジタル信号処理装置8に入力されると、ディジタル信号処理装置8が信号処理、画像再構成等の処理を実行し、その結果である被検体1の断層画像をディスプレイ20に表示すると共に、外部記憶装置の磁気ディスク18等に記録する。ディジタル信号処理装置8は、上述した画像再構成等の処理のほかに、事前計測データを用いた装置特性に関する計算や、磁気ディスク18に格納されたパルスシーケンスの実行に必要なパラメータの計算を行い、シーケンサ4の動作を制御する。
【００３１】
操作部25は、MRI装置の各種制御情報や上記信号処理系7で行う処理の制御情報を入力するもので、トラックボール又はマウス23、及び、キーボード24から成る。この操作部25はディスプレイ20に近接して配置され、操作者がディスプレイ20を見ながら操作部25を通してインタラクティブにMRI装置の各種処理を制御する。
【００３２】
なお、図1において、送信側の高周波コイル14aと傾斜磁場コイル9は、被検体1が挿入される静磁場発生系2の静磁場空間内に、垂直磁場方式であれば被検体1に対向して、水平磁場方式であれば被検体1を取り囲むようにして設置されている。また、受信側の高周波コイル14bは、被検体1に対向して、或いは取り囲むように設置されている。
【００３３】
本実施形態のMRI装置は、さらに信号処理系7の機能として、傾斜磁場コイル9への入力波形(入力傾斜磁場波形)と、それによって静磁場空間に発生する傾斜磁場波形(出力傾斜磁場波形)との関係を表す応答関数を算出する手段と、この応答関数を用いて出力傾斜磁場波形を算出する傾斜磁場波形算出手段とを備えている。これら手段を備えたディジタル信号処理装置8の詳細を図2に示す。
【００３４】
図示するように、ディジタル信号処理装置8は、傾斜磁場の応答関数を計算する応答関数算出部81、算出された応答関数を用いて、実際の撮像シーケンスで使用される傾斜磁場の出力傾斜磁場波形を算出する出力傾斜磁場波形算出部82、出力傾斜磁場波形算出部82が算出した出力傾斜磁場波形を用いて撮像シーケンスで使用する高周波パルス等を補正するシーケンス計算部83、出力傾斜磁場波形算出部82が算出した出力傾斜磁場波形を用いて、NMR信号が配置されるk空間の座標位置を補正するk空間座標算出部84、NMR信号を用いて画像再構成演算を行って画像データを作成する画像再構成部85を備えている。
【００３５】
以下、上記構成における本実施形態のMRI装置の動作、特にディジタル信号処理装置8で行われる処理を中心に説明する。
<第一の実施形態>
本実施形態のMRI装置のディジタル信号処理装置8で行われる処理は、大きく分けて、4つの処理からなる。4つの処理の関係を図2のブロック図を用いて説明する。1つめの処理101は、事前測定として応答関数算出部81で行われる処理であり、渦電流及び傾斜磁場電源の制御回路の応答関数1112を算出する処理である。処理102〜104は、実際の撮像において行われる処理である。処理102は、処理101で求めた応答関数1112と、撮像シーケンス毎に設定される入力傾斜磁場波形1113、1114を用いて、それぞれの出力傾斜磁場波形を推定する処理であり、出力傾斜磁場波形算出部82が行う。処理103は、処理102の結果を用いて、撮像シーケンス高周波磁場パルスを再計算する処理であり、シーケンス計算部83の高周波磁場パルス算出部831が行う。処理104は、処理102の結果を用いて、計測されるNMR信号1115の各サンプリング点に対応するk空間座標を再計算する処理であり、k空間座標算出部84が行う。以下、それぞれの処理について詳説する。
【００３６】
<<応答関数算出処理101>>
処理101の渦電流及び傾斜磁場電源の制御回路の応答関数算出処理について図3に示すフローチャートを用いて全体の流れを説明する。
【００３７】
ステップ301において、傾斜磁場電源の制御パラメータを無効化し、フィードバック制御やフィードフォワード制御などの動作を傾斜磁場電源が行わないように設定する。傾斜磁場電源の制御を無効にした状態で、ステップ302において、出力傾斜磁場波形を空間別に測定する。出力傾斜磁場の測定は、公知技術(例えば、特許文献3に記載される技術)を用いて行うことができる。空間別の測定は、例えば図4に示すように、X、Y、Zの各軸方向に50mm間隔で測定する。この測定結果である出力傾斜磁場波形1111及び測定時に用いた入力傾斜磁場波形1110を、信号処理系のメモリRAM22(図1)に格納する。
【００３８】
ステップ303において、RAM22に格納された入力傾斜磁場波形及び計測された出力傾斜磁場波形を読み出し、渦電流の応答関数を算出する。渦電流の応答関数を算出方法については後述する。次にステップ304において、傾斜磁場電源の制御パラメータを有効化し、フィードバック制御やフィードフォワード制御などが動作する状態にする。その状態において、出力傾斜磁場波形を測定する(ステップ305)。この測定結果もまた、RAM22に格納する。ステップ306において、傾斜磁場電源の制御回路の応答関数を算出する。傾斜磁場電源の制御回路の応答関数を算出方法についても後述する。以上の処理は、X、Y、Zの3軸方向に配された傾斜磁場コイル毎に行う。
【００３９】
まずステップ303の渦電流の応答関数の算出について説明する。渦電流により生じる磁場は、図5に例として示すように、入力傾斜磁場波形501に対して出力傾斜磁場波形502を歪ませる。入力傾斜磁場波形と出力傾斜磁場波形の関係は例えば、次のような式で表現できる。
0005916240.tif 000002
【００４０】
ここで、sはラプラス変数、Gout(s)は出力傾斜磁場波形のラプラス変換形、Gin(s)は入力傾斜磁場波形のラプラス変換形、H(s)は渦電流の伝達特性である。なお、出力傾斜磁場波形は、図4に示したように空間別に測定された個々のデータである。式(1)は一例であり、例えば下記式のような別の表現に置き換えることもできる。
0005916240.tif 000003
【００４１】
以下の説明では、式(1)の表現を用いることとする。入力傾斜磁場波形Gin(s)と出力傾斜磁場波形Gout(s)が既知であることから、理論的には渦電流の伝達特性H(s)は次式から算出できる。
0005916240.tif 000004
渦電流の応答関数は、式(3.1)を逆ラプラス変換することで得られる。
0005916240.tif 000005
【００４２】
ここで、ILT[]は逆ラプラス変換を表す。ただし、式(3.1)の計算では、Gin(s)の値が零となる場合にはH(s)の値が発散し、正しい値が求められない。このことを防ぐため、本実施形態では、渦電流の応答関数をモデル式に当てはめ、モデル式中の係数を検索して求める。渦電流の応答関数のモデル式の一例を下記式(4)に示す。この式(4)の応答関数は渦電流により発生する磁場の過渡応答を表現している。
0005916240.tif 000006
【００４３】
ここで、tは時間、α1〜αnは振幅ゲイン、τ1〜τnは時定数である。式(4)中の指数関数の数nは、MRI装置を構成する部品が持つ時定数の種類で決まる。発明者らの経験では、nは10あれば十分に渦電流により発生する磁場の過渡応答を表現することができる。
【００４４】
渦電流の応答関数を求めるために、入力傾斜磁場波形に式(4)のh(t)を畳みこんだ結果が、最も出力傾斜磁場波形に近い傾斜磁場波形を与えるα1〜αn及びτ1〜τnを検索する。検索手順を図6のフローチャートを用いて説明する。
【００４５】
まずステップ601で、振幅ゲインおよび時定数の初期値を設定する。初期値として、ここでは、振幅ゲインαの各係数を全てゼロに設定し、時定数τにはゼロ以外の正の数値を設定する。ループステップ602(繰り返し処理)では処理のループカウンタとしてiを設定する。このiは、αとτの添字の番号に対応する。すなわち、i=1の場合、α1とτ1の最適値を検索することになる。ステップ603では、αiのみを変化させ、最適なαiの値を検索する。このとき、最適値は、下記式の評価量e1が最小になる値である。
0005916240.tif 000007
【００４６】
ここで、gout(t)は出力傾斜磁場波形、gin(t)は入力傾斜磁場波形、h(t)は式(4)に示した応答関数である。式(5)は、出力傾斜磁場波形と入力傾斜磁場波形に応答関数を畳み込んだ結果がどの程度類似しているかについて数値化しているものである。従って、式(5)以外の評価式を用いても本質的に変わりはない。-1<αi<1の範囲内を0.001刻みでαiを変更してe1を計算し、e1が最小になる値になるαiを最適値とする。ここで、黄金分割法などに代表される検索アルゴリズムを用いてもよい。
【００４７】
検索した最適値αiをh(t)に設定した上で、ステップ604において、今度はτiのみを変化させ、評価量e1が最小になるτiの値を検索する。検索方法についてはαiと同様である。αiとτiを交互に最適値を検索し、その値を更新し続けた結果、αiとτiの値が収束しているかについて分岐ステップ605で判断する。αiとτiの収束条件は例えば以下のように定義できる。
0005916240.tif 000008
ここで、αi_prevとτi_prevは前回検索時に検出されたαiとτiの最適値である。すなわち、最低2回の繰り替えし検索が必要となる。
【００４８】
上記収束条件を満たさなければ、再度αiとτiの検索を行う。収束条件を満たせば、iがn未満であればiをインクリメントし、次のαiとτiの組み合わせの検索に移る。iがn以上であれば終了する。なお、αは正だけではなく負の値も取り得る。αは正の値のみになることが渦電流の発生原理からは正しいが、α1〜αnを順番に求め、後で求めるα及びτの成分を無視して検索しているが故に負の値が生じ得る。以上の手順で渦電流の応答関数を求める。
【００４９】
次に図3のステップ306の傾斜磁場電源の制御回路の応答関数の算出について説明する。傾斜磁場電源の制御回路は、図7に例として示すように入力傾斜磁場波形701に対して出力傾斜磁場波形702にオーバシュートもしくはアンダーシュートなどの挙動を含ませることになる。傾斜磁場電源が出力する傾斜磁場パルスは、入力波形にこの制御回路からの成分を加えたものであり、下記式のように表現できる。
0005916240.tif 000009
ここで、gamp(t)は傾斜磁場電源が出力する傾斜磁場パルス、gin(t)は入力傾斜磁場波形、u(t)は傾斜磁場電源の制御成分である。
【００５０】
この制御成分についても、渦電流の場合と同様にモデル式で定義する。制御成分のモデル式は、制御回路が採用している制御手法に応じて適宜定義することができる。制御手法としては、PID制御、それより簡素なP制御やPD制御、PI制御、もしくはより複雑なニューラルネットワークなどの知的制御を想定することができ、本発明を適用するMRI装置に搭載される傾斜磁場電源の仕様に合わせて制御理論を選択すればよい。本実施形態では傾斜磁場電源の制御回路がPID制御を行っている場合を想定し、u(t)を下記のモデル式で定義する。
0005916240.tif 000010
【００５１】
ここで、KPは比例ゲイン、KIは積分ゲイン、KDは微分ゲインである。また、式(9)のDiff(t)は下記式で定義される。
0005916240.tif 000011
【００５２】
ここで、delayはPID制御の入力を算出する際の時間遅れである。gamp(t)と出力傾斜磁場波形の関係は以下である。
0005916240.tif 000012
【００５３】
ここで、gout(t)は図3のステップ305でRAM22に格納された出力傾斜磁場波形である。h(t)は式(4)に示した渦電流の応答関数であり、上述した処理により算出され、既知である。そこで、入力傾斜磁場波形と傾斜磁場電源の制御成分u(t)との合計であるgamp(t)に、ステップ303で算出したh(t)を畳みこんだときに、最も出力傾斜磁場波形gout(t)に近い傾斜磁場波形を与えるu(t)を求める。u(t)に含まれるdelay、KP、KI、KDを検索する手順を図8のフローチャートを用いて説明する。ステップ801では、delay、KP、KI、KDを全てゼロに初期設定する。ステップ802では、delayのみを変化させ最適なdelayを検索する。このとき、最適値は、下記式の評価量e2が最小になる値である。
0005916240.tif 000013
【００５４】
評価量e2が最小になる値になるdelayの検索は、上述した応答関数の係数α1〜αn及びτ1〜τnの検索処理と同様に行う。すなわち、特定の範囲内で全点検索を行ってもよいし、黄金分割法などに代表される検索アルゴリズムを用いてもよい。delayを求めることで、Diff(t)が求まることから、次にステップ803において最適なKPを検索する。ここでもやはりKPが最適となる条件は、式(12)の評価量e2が最小になる値である。同様の処理で、ステップ804においてKIを検索し、ステップ805においてKDを検索する。最終的に求められる傾斜磁場電源の制御回路の応答関数は下記式で表わされる。
0005916240.tif 000014
【００５５】
以上の処理により、渦電流の応答関数h(t)と制御回路の応答関数(制御成分u(t))が算出されたので、入力傾斜磁場波形が与えられれば、式(11)より出力傾斜磁場波形を計算できることになる。式(11)を 解くための上記応答関数h(t)(式(4))およびu(t) のパラメータは、図4に示したような空間別に測定された個々のデータについてそれぞれ求められ、信号処理系7のRAM22に保存され、実際の撮像においてパルスシーケンスの計算や画像再構成に使用される。
【００５６】
式(11)を用いて出力傾斜磁場波形を推測した結果例を図9に示す。図中、波形901は入力傾斜磁場波形、波形902は実測して得られた出力傾斜磁場波形、そして波形903は式(11)を用いて推定した出力傾斜磁場波形である。波形902と波形903はよく一致していることから、式(11)により出力傾斜磁場波形を精度良く推定できることが分かる。一方、図10は、傾斜磁場電源の制御成分u(t)を加えない場合の出力傾斜磁場波形の推定結果例である。
実測して得られた出力傾斜磁場波形902と、u(t)をゼロとした式(11)から推定した出力傾斜磁場波形904との間に矢印で示した部分に偏差が生じている。このことから、傾斜磁場電源の制御成分を加味しなければ、高精度に出力傾斜磁場波形を推定できないことが分かる。
【００５７】
<<出力傾斜磁場波形算出処理102>>
次に、図2のブロック図に戻り、出力傾斜磁場波形算出部82で行われる処理102について説明する。実際の撮像にあたっては、撮像目的に応じてシーケンス計算部83により撮像パルスシーケンスが計算され、計算された撮像パルスシーケンスに含まれる高周波パルス(振幅、エンベロープ)、スライス選択傾斜磁場パルス(入力傾斜磁場波形1113)、読出し傾斜磁場パルス(入力傾斜磁場波形1114)などがRAM22に保存される。
【００５８】
処理102では、応答関数算出手段81が算出した応答関数1112と、設定された入力傾斜磁場波形1113、1114を用いて、撮像パルスシーケンスの出力傾斜磁場波形を算出する。図11に、出力傾斜磁場波形の算出処理のフローチャートを示す。
【００５９】
まずステップ1101において、撮像視野の座標位置からX、Y、Zの3軸方向における視野中心位置を取得する。視野中心位置は、下式から算出される。
0005916240.tif 000015
ここで、AxisはX、Y、Zの各軸を表し(以下、同じ)、FOV(Axis,min)とFOV(Axis,max)は撮像条件に応じて指定される各軸の撮像視野の最小座標位置と最大座標位置である。
【００６０】
ステップ1102において、RAM22に保存された空間別の応答関数h(t)及び傾斜磁場電源の制御成分u(t)から、ステップ1101で取得した各軸の視野中心位置に対応する応答関数h(t)および制御成分u(t)を選択し、取得する。ここで、応答関数h(t)及び傾斜磁場電源の制御成分u(t)の取得は、視野中心位置に最も近い位置で算出された結果を取得してもよいし、空間別に算出された結果を補間して応答関数h(t)及び傾斜磁場電源の制御成分u(t)を求めてもよい。補間する場合は、下記式に従う。
0005916240.tif 000016
ここで、h(Axis,t)とu(Axis,t)は各軸の視野中心位置における応答関数と傾斜磁場電源回路の制御成分である。また、P1およびP2は、視野中心位置の前後方向でそれぞれ最も近い応答関数および傾斜磁場電源回路の制御成分の算出位置を表す。
【００６１】
次に、ステップ1103において、X、Y、Zの3軸方向に対する入力傾斜磁場波形を取得する。入力傾斜磁場波形は、撮像パルスシーケンスに応じて設定され、RAM22に保存されているものであり、図2に示される、スライス傾斜磁場パルス1113もしくは読み出し傾斜磁場パルス1114である。どちらの場合であっても、処理の内容は同じである。
【００６２】
ステップ1104において、次式で表される3軸の入力傾斜磁場波形を合成した波形を算出する。
0005916240.tif 000017
ここで、gin(Axis,t)は取得した各軸の入力傾斜磁場波形である。
【００６３】
ステップ1105において、軸毎の入力傾斜磁場波形に対して応答関数h(t)及び傾斜磁場電源の制御成分u(t)を用いて出力傾斜磁場波形を推定する。出力傾斜磁場波形の推定は、式(11)を3軸に拡張した形として下記式で表わされる。
0005916240.tif 000018
【００６４】
ここで、gest(Axis,t)は計算により求めた各軸の出力傾斜磁場波形の推定結果である。
各軸の出力傾斜磁場波形を算出した後に、ステップ1106において、下記式に従い各軸の波形を合成する。
0005916240.tif 000019
以上のように推定されたスライス傾斜磁場及び読出し傾斜磁場の各出力傾斜磁場波形は、それぞれ、高周波磁場パルスの再計算(図2：103)およびNMR信号のk空間座標位置の補正(図2：処理104)に用いられる。
【００６５】
<<高周波磁場パルス算出処理103>>
処理103では、出力傾斜磁場波形算出部82で求めたスライス傾斜磁場の出力傾斜磁場波形gest(t)と入力傾斜磁場波形gin(t)を用いて、高周波磁場パルスの再計算を行う。高周波磁場パルスの再計算は下記式に従い、高周波パルスのエンベロープおよび周波数を計算する。
0005916240.tif 000020
ここで、RFenv'(t)は再計算後の高周波磁場パルスのエンベロープ、RFenv(t)は再計算前の高周波磁場パルスのエンベロープ、RFfreq'(t)は再計算後の高周波磁場パルスの周波数、RFfreq(t)は再計算前の高周波磁場パルスの周波数である。
【００６６】
これらの式(19)〜(21)からわかるように、ここでは、スライス傾斜磁場の出力傾斜磁場波形の強度に比例して高周波磁場パルスのエンベロープの振幅、周波数、照射時間を変更している。これら式を基本として、装置の制約や撮像条件を考慮した変更を行うことができる。
【００６７】
変更を要する一例として、高周波磁場パルスの周波数を時間に応じて変化させることが装置の構成上困難な場合がある。その場合には、高周波磁場パルスの位相を時間と共に変化させ、結果として周波数の変調を実現してもよい。周波数変調を実現するための位相の算出は次式に従う。
0005916240.tif 000021
ここで、g_standardは高周波磁場パルスの周波数計算の基準となる傾斜磁場強度[T/m]、distanceは磁場中心から励起位置までの距離[m]、γは磁気回転比[Hz/T]、Δtは高周波磁場パルスの位相変更間隔時間[s]である。
【００６８】
ただし、位相を変更できる範囲は0〜360[deg]であるため、distanceの範囲もまた限られる。従って、distanceの設定可能な範囲を広くするためには、位相の変化は少ない方が好ましい。そこで、位相の変化を最小とするために、g_standard及び高周波磁場パルスの周波数を下記式に従って計算する。
0005916240.tif 000022
ここで、Average[  ]は高周波磁場パルス印加区間内での平均値計算処理、RFfreq'は時間的変化のない再計算後の高周波磁場パルスの周波数、RFfreqは時間的変化のない再計算前の高周波磁場パルスの周波数である。式(23)及び式(24)によりg_standardとRFfreq'を算出することで、位相の変化が最小に抑えられ、distanceの上限値を大きく取ることが可能となる。
【００６９】
変更を要する他の例として、高周波磁場パルスの照射時間の抑制がある。上記式(19)〜式(21)を用いた高周波磁場パルスの再計算の結果として、高周波磁場パルスの照射時間は変化するが、パルスシーケンスによっては、照射時間の変化、特に照射時間の延長が好ましくない場合がある。照射時間を変化させない高周波磁場パルスの再計算手順を図12に示すフローチャートを用いて説明する。
【００７０】
ステップ1021と1022において、軸毎の出力傾斜磁場波形の算出、各軸の出力傾斜磁場波形の合成を行う。これらの処理は、出力傾斜磁場波形算出部82で行われるものであり、図11のステップ1103とステップ1104で行う処理と同じである。次に、ステップ1031において、式(19)〜式(21)に従い高周波磁場パルスを計算する。その結果を受け、ステップ1032において、照射時間の増加率βRFを算出する。増加率βRFの算出は次式に従う。
0005916240.tif 000023
ここで、gin(t)は入力傾斜磁場波形、gest(t)は式(18)より求めた出力傾斜磁場波形の推定結果、StartTimeは高周波磁場パルス照射開始時刻、EndTimeは高周波磁場パルス照射終了時刻を表す。
【００７１】
分岐ステップ1033において、増加率βRFが1.0を越えるか否かを判断する。増加率βRFが1.0を越える場合、計算後の高周波磁場パルスの照射時間が延長していることになる。
その場合にはステップ1034に移り、スライス傾斜磁場の入力傾斜磁場波形の再作成を行う。入力傾斜磁場波形の再作成は次式に従う。
0005916240.tif 000024
照射時間の増加率βRFを入力傾斜磁場波形に掛けることで、入力傾斜磁場波形の強度を増加させ、それにより高周波磁場パルスの印加時間を短縮させる。この時、入力傾斜磁場波形の全体を一様に増加させるだけであるため、特別なシムコイルや傾斜磁場コイルなどは必要としない。
【００７２】
入力傾斜磁場波形の再作成後、その結果を入力として処理102のステップ1021に戻し、出力傾斜磁場波形の算出(ステップ1021、1022)、高周波磁場パルスの再計算(ステップ1031)、照射時間の増加率βRFの算出(ステップ1032)を行う。再び算出した増加率βRFが1.0以下であれば処理を終了し、1.0以下になっていなければ再度、入力傾斜磁場波形を作成する。このとき、式(26)のgin(t)は直前に再作成された入力傾斜磁場波形である。
【００７３】
このようにスライス傾斜磁場の入力傾斜磁場波形の再作成を繰り返すことで、最終的に高周波磁場パルスの照射時間は、再計算前の高周波磁場パルスの照射時間に比べて延長しないようにすることができる。これにより、再計算後の高周波磁場パルスを使用する際に、スライス傾斜磁場パルス以外のパルスシーケンスを変更する必要がなくなる。
【００７４】
図13は、再計算前の高周波磁場パルスを用いて得られた画像1301と、再計算後の高周波磁場パルスを用いて得られた画像1302である。どちらの画像も、VERSE法を用いた超短TE計測での撮像結果である。再計算前の高周波磁場パルスを用いて得られた画像1301に比べて、再計算後の高周波磁場パルスを用いて得られた画像1302は、矢印で示す背景信号(撮影断面以外からの信号)の混入が少ないことから、良好な励起プロファイルが得られていることが分かる。
【００７５】
<<k空間座標再計算処理104>>
続いて、図2のk空間座標算出部84が行うk空間座標再計算処理104について説明する。k空間座標再計算処理104では、出力傾斜磁場波形算出部82で求めた読出し傾斜磁場の出力傾斜磁場波形gest(t)と入力傾斜磁場波形gin(t)を用いてNMR信号が配置されたk空間座標の再計算を行う。図14に、k空間座標の再計算手順を示すフローチャートの一例を示す。この処理手順も、図12に示す高周波磁場パルスの再計算のフローチャートと同様の処理の流れで、k空間座標の再計算結果によって読出し傾斜磁場の入力傾斜磁場波形を再作成する。
【００７６】
まずステップ1041において、処理102において求めた、読出し傾斜磁場の入力傾斜磁場波形gin(t)と出力傾斜磁場波形gest(t)を用いて、k空間座標の算出を行う。出力傾斜磁場波形gest(t)を用いたk空間座標の再計算は下記式に従う。
0005916240.tif 000025
ここで、γは磁気回転比[Hz/T]である。この式は、一般的な入力傾斜磁場波形gin(t)を用いたk空間座標の算出の式(次式(28))におけるgin(j)をgest(j)に置き換えたものである。
0005916240.tif 000026
【００７７】
次に、再計算によって算出した座標の範囲が縮小しているか否かを調べる。再計算したk空間座標k’(t)は、もともとのk空間座標k(t)に対して、その範囲が縮小する可能性がある。その場合、k空間内の想定された範囲内に信号が充填されないことになり、画像の分解能を低下させる。そこで、ステップ1042において、k空間の縮小率βkを算出する。
縮小率βkの算出は次式に従う。
0005916240.tif 000027
ここで、EndTimeは読み出し傾斜磁場パルスの印加終了時刻(もしくはNMR信号の計測終了時刻)を表す。
【００７８】
分岐処理ステップ1043において、縮小率βkが1.0未満であるか否かを判断する。縮小率βkが1.0未満である場合、再計算後のk空間座標k’(t)は再計算前のk空間座標k(t)に対して縮小していることになる。その場合にはステップ1044に移り、読み出し傾斜磁場の入力傾斜磁場波形の再作成を行う。入力傾斜磁場波形の再作成は次式に従う。
0005916240.tif 000028
【００７９】
照射時間の縮小率βkの逆数を入力傾斜磁場波形に掛けることで、入力傾斜磁場波形の強度を増加させ、それによりk空間座標の範囲を広げる。この時、式(30)からわかるように、入力傾斜磁場波形は全体を一様に増加させられるだけであるため、特別なシムコイルや傾斜磁場コイルなどは必要としない。入力傾斜磁場波形の再作成後、その結果を入力として処理102に戻し、出力傾斜磁場波形の算出(ステップ1021、1022)、k空間座標の再計算(ステップ1041)、縮小率βkの算出(ステップ1042)を行う。
【００８０】
分岐処理ステップ1043において、再び算出した縮小率βkが1.0以上であれば処理を終了し、1.0以上になっていなければ再度、入力傾斜磁場波形を作成する。このとき、式(30)のgin(t)は直前に再作成された入力傾斜磁場波形である。
【００８１】
このように入力傾斜磁場波形の再作成を繰り返すことで、最終的にk空間座標の範囲は縮小せず、画像の分解能を保つことができる。図15は、再計算前のk空間座標を用いて得られた画像1501と、再計算後のk空間座標を用いて得られた画像1502である。どちらの画像も、EPI計測での撮像結果であり、読み出し傾斜磁場は画像における横方向に印加されている。再計算前のk空間座標を用いて得られた画像1501に比べて、再計算後のk空間座標を用いて得られた画像1502は、矢印で示す画像の歪みが低減されていることから、読み出し傾斜磁場の歪みの補償が適切に行われていることが分かる。
【００８２】
撮像は、上記処理103で修正された高周波磁場パルスを用いて行われ、また上記処理104で修正されたk空間座標を用いて、取得したNMR信号の画像再構成を行う。これにより選択したスライスへの外部信号の混入がなく、且つ画像の歪みやアーチファクトの抑制された画像を得ることができる。
【００８３】
以上、説明した本実施形態の、従来技術に対する利点は、次のとおりである。
まず、渦電流を補償する傾斜磁場パルスを出力する方法(特許文献1、特許文献2)では、時間的・空間的に高次の複雑な変化示す渦電流による磁場を完全に補正することは困難かつ高コストであるのに対し、本実施形態では特別な磁場補正用コイルを必要とせず、渦電流による磁場を含めた高周波磁場パルスとk空間座標の再計算を行うことで、励起プロファイルの劣化と画像の歪みやゴースト等のアーチファクトを抑制できる。
【００８４】
また、高周波磁場パルスの照射タイミングを出力傾斜磁場波形に応じて変更する方法(特許文献3)では、高周波磁場パルスの照射タイミングのみを変更しているのに対して、本実施形態では渦電流及び傾斜磁場電源の制御回路の応答関数を用いて出力傾斜磁場波形を推定するため、出力傾斜磁場波形の非線形な歪みについても補償する事が可能である。
【００８５】
固定の入力傾斜磁場波形を用いる方法(非特許文献1)では、特定のパルスシーケンスのみにしか適用できないのに対して、本実施形態ではあらゆる入力傾斜磁場波形に対する出力傾斜磁場波形を推定できるため、固定の入力傾斜磁場波形を使用しなければならないといった制限がない。
【００８６】
さらに、システム応答をRLC回路でモデル化する方法(非特許文献2)では、傾斜磁場電源の制御回路の振る舞いについては考慮されていないが、本実施例では傾斜磁場電源の制御回路の振る舞いも考慮するため、より高精度に高周波磁場パルス及びk空間座標を再計算できる。加えて、システム応答をRCL回路でモデル化する方法(非特許文献2)では、画像の見た目を判断基準としてモデル式の各係数を求めているが、本実施例では出力傾斜磁場波形に基づいてモデル式の係数を算出するため、人間の判断に依る必要がなく、精度の高いモデル式を算出することができる。
【００８７】
また、撮像後に傾斜磁場波形を実測して再構成に用いる方法(非特許文献3)では、ノイズによる揺らぎが重畳することになるが、本実施例ではモデル式を用意して渦電流と傾斜磁場電源回路の制御成分を表現しているため、推定する出力傾斜磁場波形にノイズが重畳しない。
【００８８】
<第二の実施形態>
本実施形態においても、図2に示すディジタル信号処理部8の構成と、ディジタル信号処理部8で4つの特徴的な処理101〜104が行われることは、第一実施形態と同様である。
本実施形態は、渦電流及び傾斜磁場電源の制御回路の応答関数を算出する処理101を、第一実施形態より簡易な手続きで求めることを特徴としている。すなわち、第一実施形態では、傾斜磁場電源回路の制御を無効にした状態と有効にした状態でそれぞれ2回の出力傾斜磁場波形の測定を行ったが、本実施形態では、1回の出力傾斜磁場波形の測定のみ行う。本実施形態の処理の流れを図16のフローチャートを用いて説明する。
【００８９】
まず、ステップ1601において、出力傾斜磁場波形を空間別に測定する。この時、傾斜磁場電源の制御パラメータには変更を加えず、装置の通常状態としておく。すなわち、傾斜磁場電源の制御パラメータが有効になっている状態とする。出力傾斜磁場波形の測定結果は、RAM22に格納する。この処理は、第一実施形態のフローを示す図3のステップ302、305で行われる処理と同様である。
【００９０】
次に、ステップ1602において、RAM22に格納された入力傾斜磁場波形及び計測された出力傾斜磁場波形を読み出し、渦電流と傾斜磁場電源の制御回路の応答関数を算出する。
【００９１】
渦電流と傾斜磁場電源の制御回路の応答関数の算出について、図17のフローチャートを用いて説明する。ステップ1701において渦電流の応答関数を算出する。この処理は上述した第一実施形態でのステップ303の処理と同様である。すなわち、次式(31)の評価量e1が最小になる応答関数h(t)を検索する。
0005916240.tif 000029
ここで、gout(t)はステップ1601でRAM22に格納された出力傾斜磁場波形であり、gin(t)は入力傾斜磁場波形である。h(t)は式(4)で定義される渦電流の応答関数であり、第一実施形態でのステップ303と同様に式(4)の係数α1〜αn、τ1〜τnを検索する。
【００９２】
式(31)は、第一実施形態の式(5)と同形状であるが、式(5)の出力傾斜磁場波形gout(t)は傾斜磁場電源の制御回路を無効にして測定したものであるのに対し、式(31)の出力傾斜磁場波形gout(t)には、傾斜磁場電源の制御回路による歪みが含まれている。そのため、算出される渦電流の応答関数h(t)には、渦電流だけではなく傾斜磁場電源の制御回路の特性も含まれることになる。
【００９３】
次に、ステップ1702において傾斜磁場電源の制御回路の応答関数を算出する。この処理は第一実施形態のステップ306の処理と同じである。すなわち式(32)の評価量e2が最小になるgamp(t)を検索する。式(32)のgamp(t)は、式(33)で定義される。
0005916240.tif 000030
ここでu(t)は、式(9)で定義される傾斜磁場電源の制御成分であるが、ステップ1701で求めた渦電流の応答関数h(t)の中には、傾斜磁場電源の制御回路の応答成分も含まれているため、渦電流の応答関数のモデル(式(4))では表現しきれなかった成分を式(9)に当てはめて算出することになる。
【００９４】
ステップ1702で傾斜磁場電源の制御回路の応答関数u(t)を求めた後に、ステップ1703において渦電流の応答関数h(t)を再計算する。これは、傾斜磁場電源の制御回路の応答関数u(t)を定義することで、電流の応答関数h(t)(式(4))の各係数(α1〜αn及びτ1〜τn)の最適値が変化するためである。応答関数h(t)を再計算は、式(31)のgin(t)を、ステップ1702で求めたgamp(t)に置き換えて、その評価量e1を最小とするh(t) の各係数を検索することにより行う。以上の処理により、渦電流の応答関数と傾斜磁場電源の制御成分を加味した応答関数が求められる。応答関数は、空間別に測定された個々のデータについてそれぞれ求められ、信号処理系7のRAM22に保存される。
【００９５】
撮像に際し、RAM22に保存した応答関数を用いて、出力傾斜磁場波形の推定する処理102、及び高周波磁場パルスやk空間座標を再計算する処理103、104の手順は第一実施形態と同様であり、説明を省略する。
【００９６】
本実施形態によれば、応答関数の算出回数が増えるものの、出力傾斜磁場波形の測定回数を減らすことができ、傾斜磁場電源の制御パラメータを変更する作業も省くことができ、より簡素な手順で励起プロファイルを再計算するために必要な応答関数を算出できる。
その他の効果は、第一実施形態と同様である。
【００９７】
<第三実施形態>
本実施形態においても、図2に示すディジタル信号処理部8の構成と、ディジタル信号処理部8で4つの処理101〜104が行われることは、第一実施形態と同様である。本実施形態は、渦電流及び傾斜磁場電源の制御回路の応答関数を特定のモデル式に当てはめずに求めることを特徴としている。本実施形態の処理の流れを図18に示す。図18のステップ1801、1802は、第二実施形態で説明した図16のステップ1601、1602と同様であるが、本実施形態の特有の処理は、ステップ1801において、複数種類の入力傾斜磁場波形を用いて出力傾斜磁場波形を測定する処理、ステップ1802において、モデル式を用いずに渦電流と傾斜磁場電源の制御回路の応答関数を直接導出する処理、ステップ1803において、導出した応答関数を補正する処理である。
【００９８】
既に説明したように、入力傾斜磁場波形Gin(s)(ラプラス変換形)と出力傾斜磁場波形Gout(s) (ラプラス変換形)が既知であれば、理論的には渦電流の伝達特性H(s)は式(3.1)から算出できる。式(3.1)を再掲する。
0005916240.tif 000031
しかし、Gin(s)の値が零となる場合にはH(s)の値が発散し、正しい値が求められない。
これを防ぐため、第一および第二実施形態では、渦電流の応答関数をモデル式に当てはめて、最適な評価値を与える応答関数を導出した。
【００９９】
本実施形態では、モデル式を用いることなく、式(3.1)或いはそれを逆ラプラス変換した式(3.2)から算出される応答関数h(t)を用いる。この際、零点の影響を排除するために、渦電流と傾斜磁場電源制御回路の応答関数の算出に用いる入力傾斜磁場波形を複数種類用意する。
【０１００】
具体的には、ステップ1801において、複数種類の入力傾斜磁場波形を用いて、それぞれ出力傾斜磁場波形を測定し、RAM22に保存する。入力傾斜磁場波形なるべく急峻な形状を持ったものがよい。なだらかな変化の波形では、高周波域が含まれず、高域の応答を算出できないためである。理想的には、インパルス波を用いることが望ましいが、傾斜磁場電源の出力性能などの制約によりインパルス波を出力することが難しいため、本実施形態では、図19に示すような複数の台形波を用いる。
【０１０１】
ステップ1802では、式(3.1)により伝達特性H(s)を算出する。その際、一つの入力傾斜磁場波形(台形波)のスペクトル(ラプラス変換形)に現れる零点を、他の入力傾斜磁場波形(台形波)のスペクトルデータを用いて置き換えることにより、零点を回避する。具体的には、次式に従う。
0005916240.tif 000032
ここで、Gin1(s)は、台形波1901のラプラス変換形、Gout1(s)は、入力傾斜磁場波形として台形波1901を用いた場合の出力傾斜磁場波形のラプラス変換形、Gin2(s)は、台形波1902のラプラス変換形、Gout2(s)は、入力傾斜磁場波形として台形波1902を用いた場合の出力傾斜磁場波形のラプラス変換形である。
【０１０２】
図19に示す台形波群をラプラス変換して得られる振幅スペクトルの概略を図20に図示する。図20に示すように、台形波1901、1902のスペクトル2001、2002に零点が現れる位置は、スペクトル2001とスペクトル2002とで異なるので、式(34)に示すように、値が大きい方のスペクトルを用いて計算することにより、零点を回避することができる。なお、図19では台形波を2つ示しているが、3つ以上の台形波を用意してもよい。その場合には、式(34)において、複数の入力傾斜磁場波形のうち、スペクトルの値が最も大きい(零点ではない)値となる入力傾斜磁場波形とそれによる出力傾斜磁場波形の組合せを選んで伝達特性の算出に用いればよい。
【０１０３】
次に、ステップ1803において、ステップ1802で算出した応答関数h(t)からノイズの影響を排除する処理を行う。これは、Gout1(s)及びGout2(s)は実測データであるため、ノイズが重畳しており、正しい伝達特性が求められないケースが存在するためである。本実施形態では、このようなノイズを、ラプラス空間で2段階のフィルタ処理を行うことにより除去する。
【０１０４】
まず式(34)で求めた伝達特性に、次式の最適フィルタを掛け算し、ノイズの影響を抑制する。
0005916240.tif 000033
ここで、H'(s)は最適フィルタ適用後の伝達特性、Φ(s)は最適フィルタである。最適フィルタΦ(s)には、公知の技術であるウィナーフィルタを用いる。ウィナーフィルタの定義式は下記である。
0005916240.tif 000034
ここで、N(s)はGout1(s)のノイズ成分、S(s)はGout1(s)のノイズを除いた信号成分である。N(s)及びS(s)の求め方は、既知の技術に従い、Gout1(s)のパワースペクトルの最も高周波の信号値をN(s)とし、Gout1(s)のパワースペクトルからN(s)を減算したものをS(s)とする。
【０１０５】
このような最適フィルタによりある程度ノイズの影響は低減できるが、実用的には最適フィルタ適用後にもノイズの影響が問題になることがある。そこで、最適フィルタ適用後の伝達特性に対してノイズ除去フィルタを実行する。本実施例ではノイズ除去フィルタとしてメディアンフィルタを用いる。
0005916240.tif 000035
ここで、H”(s)はメディアンフィルタ適用後の伝達特性である。
【０１０６】
最終的に、渦電流と傾斜磁場電源制御回路を合わせた応答関数は、式(37)を逆ラプラス変換することで得られる。
0005916240.tif 000036
ここで、ILT[]は逆ラプラス変換を表す。
【０１０７】
以上のようにして求めた応答関数を用いて、出力傾斜磁場波形を推定し、高周波磁場パルス・k空間座標を再計算する手順は第一実施形態と同様である。
【０１０８】
このように本実施形態では、実測データに含まれるノイズをラプラス空間における2つのフィルタ処理により低減することで、応答関数を特定のモデル式に当てはめることなく求めることができ、様々なハードウェア構成(傾斜磁場コイル等)やソフトウェア構成(傾斜磁場電源の制御方式)に対して柔軟に適用することができる。
【０１０９】
以上、本発明のMRI装置の動作、特に傾斜磁場の応答関数算出部の動作の各実施形態を説明したが、本発明は、上述した実施形態に限定されず、種々の変更を行うことが可能である。例えば、渦電流による出力傾斜磁場波形の歪みや傾斜磁場電源の制御回路による出力傾斜磁場波形の歪みを表現するモデル式は、実施形態で例示したものに限らず、他のモデル式を用いることも可能である。また上述した実施形態では、出力傾斜磁場波形を歪ませる原因として、渦電流と傾斜磁場電源の制御回路を取り上げ、それらによる応答関数を出力傾斜磁場波形の算出に用いる場合を説明したが、その他の要素、例えば、永久磁石を用いたMRI装置で見られるヒステリシスを伴う残留磁場といったような他の成分を出力傾斜磁場波形の算出に加えることもできる。また、各実施形態の説明のために示したフローチャートの処理手順は一例であり、省略できる処理もあるし、また必要に応じて別の処理を追加することも可能である。
【０１１０】
本発明は、スライス傾斜磁場について算出した出力傾斜磁場波形により高周波磁場パルスを計算する場合や、読出し傾斜磁場について算出した出力傾斜磁場波形により、NMR信号が配置されるk空間座標を補正する場合に好適に適用することができるが、位相エンコード傾斜磁場についても本発明の適用を妨げるものではない。また読出し傾斜磁場には、読出し傾斜磁場と位相エンコード傾斜磁場の区別のないラジアルスキャンやスパイラルスキャンのようなパルスシーケンスにおける2方向の傾斜磁場も含まれる。
【０１１１】
さらに上記実施形態では、スライス傾斜磁場及び読出し傾斜磁場の両方について出力傾斜磁場波形を推定し、高周波磁場パルスの再計算およびk空間座標の再計算をともに行う場合を説明したが、いずれか一方のみを行う場合も本発明に含まれる。
【０１１２】
以上の各実施形態の説明で明らかになる本願発明の概要は以下の通りである。即ち、
本発明のMRI装置は、
静磁場発生手段、傾斜磁場発生手段、高周波磁場発生手段および高周波磁場検出手段を備えた撮像手段と、撮像シーケンスに基き前記撮像手段の動作を制御する制御手段と、前記高周波磁場検出手段が検出した核磁気共鳴信号を用いて、画像再構成を含む演算を行う演算手段と、を備え、
前記演算手段は、前記傾斜磁場発生手段に入力される入力波形と当該入力波形によって前記傾斜磁場発生手段が発生した出力波形とを用いて、前記出力波形に影響を与える複数の要素の応答関数の合算である応答関数を算出する応答関数算出手段と、前記応答関数計算手段が計算した応答関数を用いて、前記撮像シーケンスに設定された傾斜磁場パルスの入力波形から出力波形を算出する傾斜磁場算出手段と、前記傾斜磁場算出手段が算出した出力波形の計算値を用いて、前記撮像シーケンスに設定された高周波磁場パルスの修正および／または前記撮像シーケンスで得られた核磁気共鳴信号のk空間座標の修正を行う補償手段とを備えたことを特徴とする。
【０１１３】
好ましくは、
前記応答関数算出手段は、前記出力波形を変化させる複数の要素毎に応答関数を算出し、前記傾斜磁場算出手段は、前記複数の要素の応答関数の合算である応答関数を用いて出力波形の算出を行う。
【０１１４】
また、好ましくは、
前記応答関数算出手段は、前記要素毎に応答関数モデルを設定するモデル設定手段と、前記モデル設定手段が設定した応答関数モデルによって算出した出力波形の計算値と、出力波形の実測値との差が最小となるように応答関数モデルを適正化するモデル適正化手段とを備え、前記モデル適正化手段が適正化した応答関数モデルを各要素の応答関数とする。
【０１１５】
また、好ましくは、
前記応答関数算出手段は、前記入力値として複数の入力値を用い、複数の入力値とそれに対応する実測値との組合せを用いて、入力値と出力波形の実測値との関係から、直接、応答関数を算出する。
【０１１６】
また、好ましくは、
前記応答関数算出手段は、算出後の応答関数にフィルタをかけるフィルタ手段を備える。
【０１１７】
また、好ましくは、
前記フィルタ手段は、ウィナーフィルタ、メディアンフィルタを含む。
【０１１８】
また、好ましくは、
前記傾斜磁場発生手段は、発生した傾斜磁場の出力波形に基き入力波形を制御する制御回路を備え、
前記出力波形を変化させる複数の要素は、前記制御回路による制御に伴う前記出力波形の歪みを含む。
【０１１９】
また、好ましくは、
前記応答関数算出手段は、前記静磁場発生手段が発生する静磁場の中心に対する位置が異なる複数の傾斜磁場について、それぞれ、応答関数を算出する。
【０１２０】
また、好ましくは、
前記傾斜磁場算出手段が算出する傾斜磁場パルスは、前記撮像シーケンスにおけるスライス選択傾斜磁場であり、前記補償手段は、前記スライス選択傾斜磁場の計算値に基き、前記撮像シーケンスに設定された高周波パルスを補正する高周波パルス補正手段を備える。
【０１２１】
また、好ましくは、
前記高周波パルス補正手段は、補正前後の高周波パルスの照射時間の変化を算出する照射時間計算手段を備え、前記照射時間計算手段が算出した照射時間の変化が所定の増加率を超えるときに、前記傾斜磁場算出手段の入力値を変更し、前記傾斜磁場算出手段が前記変更された入力値から推定した出力波形の推定値を用いて、前記高周波パルスの補正を行う。
【０１２２】
また、好ましくは、
前記傾斜磁場算出手段が算出する傾斜磁場パルスは、前記撮像シーケンスにおける読出し傾斜磁場であり、前記補償手段は、前記読出し選択傾斜磁場の計算値に基き、前記撮像シーケンスで得た核磁気共鳴信号が配置されるk空間座標を補正するk空間座標補正手段を備る。
【０１２３】
また、好ましくは、
前記補正手段は、補正前後のk空間座標の大きさの変化を算出するk空間サイズ計算手段を備え、前記k空間サイズ計算手段が計算したk空間の大きさの変化が所定の縮小率を超えるときに、前記傾斜磁場算出手段の入力値を変更し、前記傾斜磁場算出手段が変更後の入力値から算出した出力波形の再計算値を用いて、前記k空間座標の補正を行う。
【０１２４】
本発明の磁気共鳴イメージング装置の傾斜磁場発生手段が発生する出力傾斜磁場波形を推定する方法は、
前記傾斜磁場発生手段の入力値およびそれによって前記傾斜磁場発生手段が発生した出力傾斜磁場波形の実測値を入力するステップと、
前記出力傾斜磁場波形に影響を与える複数の因子について、それぞれ、応答関数モデルを設定するステップと、
前記入力した傾斜磁場波形の入力値および実測値を用いて、各応答関数モデルを適正化するステップと、
所望の傾斜磁場の入力値に対し、前記適正化した応答関数モデルを用いて、出力傾斜磁場波形の推定値を算出するステップと、
を含むことを特徴とする。
【０１２５】
好ましくは、
前記複数の因子は、傾斜磁場により発生する渦電流と、前記傾斜磁場発生手段に備えられた制御回路の出力とを含む。
【０１２６】
また、好ましくは、
前記実測値は、前記傾斜磁場発生手段の制御回路を無効にした時の第1の実測値と、前記制御回路を有効にした第2の実測値を含み、
前記応答関数モデルを適正化するステップは、前記第1の実測値および第2の実測値を用いて、制御回路による出力傾斜磁場波形の応答関数を求める。
【０１２７】
本発明の磁気共鳴イメージング装置の傾斜磁場発生手段が発生する出力傾斜磁場波形を推定する別の方法は、
前記傾斜磁場発生手段の少なくとも2種の入力値およびそれらによって前記傾斜磁場発生手段が発生した出力傾斜磁場波形の実測値を入力するステップと、
前記少なくとも2種の入力値およびそれに対応する実測値を用いて、前記出力傾斜磁場波形に影響を与える複数の因子の影響を含む前記傾斜磁場発生手段の応答関数を算出するステップと、
所望の傾斜磁場の入力値に対し、前記応答関数を用いて、出力傾斜磁場波形の推定値を算出するステップと、を含む。
【産業上の利用可能性】
【０１２８】
本発明によれば、MRI装置において、多大な負担を伴うことなく、種々の要因によって生じる出力傾斜磁場波形の歪みを高精度に補償することができ、これにより、画像の歪みやゴースト等のアーチファクトのない信頼性の高い画像を提供することができる。
【符号の説明】
【０１２９】
1  被検体、2  静磁場発生系、3  傾斜磁場発生系、4  シーケンサ、5  送信系、6  受信系、7  信号処理系、8  ディジタル信号処理部(CPU)、9  傾斜磁場コイル、10  傾斜磁場電源、11  高周波発信器、12  変調器、13  高周波増幅器、14a  高周波コイル(送信コイル)、14b  高周波コイル(受信コイル)、15  信号増幅器、16  直交位相検波器、17  A/D変換器、18  磁気ディスク、19  光ディスク、20  ディスプレイ、21  ROM、22  RAM、23  トラックボール又はマウス、24  キーボード、25  操作部、81  応答関数算出部、82  傾斜磁場波形算出部、83  シーケンス計算部、84  k空間座標算出部、85  画像再構成部。
【図１】
0005916240.tif 000037
【図２】
0005916240.tif 000038
【図３】
0005916240.tif 000039
【図４】
0005916240.tif 000040
【図５】
0005916240.tif 000041
【図６】
0005916240.tif 000042
【図７】
0005916240.tif 000043
【図８】
0005916240.tif 000044
【図９】
0005916240.tif 000045
【図１０】
0005916240.tif 000046
【図１１】
0005916240.tif 000047
【図１２】
0005916240.tif 000048
【図１３】
0005916240.tif 000049
【図１４】
0005916240.tif 000050
【図１５】
0005916240.tif 000051
【図１６】
0005916240.tif 000052
【図１７】
0005916240.tif 000053
【図１８】
0005916240.tif 000054
【図１９】
0005916240.tif 000055
【図２０】
0005916240.tif 000056