JP3842607B2 - Spectral estimation method and magnetic resonance imaging apparatus - Google Patents
Spectral estimation method and magnetic resonance imaging apparatus Download PDFInfo
- Publication number
- JP3842607B2 JP3842607B2 JP2001311917A JP2001311917A JP3842607B2 JP 3842607 B2 JP3842607 B2 JP 3842607B2 JP 2001311917 A JP2001311917 A JP 2001311917A JP 2001311917 A JP2001311917 A JP 2001311917A JP 3842607 B2 JP3842607 B2 JP 3842607B2
- Authority
- JP
- Japan
- Prior art keywords
- value
- waveform
- spectrum
- data
- magnetic resonance
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、いわゆる自己回帰モデルを用いたスペクトル推定方法に関し、特に、Ｓ／Ｎ比の低い観測データに関して精度よくスペクトル推定をおこなうスペクトル推定方法、スペクトル推定方法を用いた磁気共鳴映像撮像方法、スペクトル推定方法を用いた磁気共鳴映像撮像装置、およびこれらの方法をコンピュータ上で実行するプログラムに関する。
【０００２】
【従来の技術】
従来、撮像対象物を磁場雰囲気中に配置し、核磁気共鳴現象に起因して撮像対象物から放出されるＮＭＲ信号を観測することにより波形データとして取り込み、その波形データを、周波数を変数とした信号強度分布であるスペクトルに関するデータに変換し、撮像対象物の内部構造を撮像する磁気共鳴映像撮像（以下、「ＭＲＩ」と言う）装置が知られている。ＭＲＩ装置によって撮像された画像の解像度を向上させるためには、ＮＭＲ信号からなる波形データをスペクトルに変換する工程において、ピーク周波数の特定を精度よくおこなうことが必要である。
【０００３】
観測によって得られた波形データについて、波形データがどのようなスペクトルを有するかを知るためにおこなわれる手法として、高速フーリエ変換が一般に用いられている。高速フーリエ変換とは、観測された波形データのうち、間隔Δｔで離散化したＮ個の時刻ｔ1、ｔ2、ｔ3、・・・、ｔi、・・・、ｔNにおける強度ｘ（ｔi）をピックアップした上で、離散フーリエ変換をおこなうことを基礎とするアルゴリズムである。高速フーリエ変換を用いてスペクトルの推定をおこなった場合、観測された波形データから観測波の周波数分布を容易に求めることができ、波形データに周期性があればスペクトルにおいて、ある一定の周波数についてピークが現れ、波形データを構成するそれぞれの波の周波数を求めることができる。
【０００４】
しかし、高速フーリエ変換によってスペクトル推定をおこなった場合にはピーク周波数を推定する精度に限界があるため、一定の範囲Δωにわたって誤差が生ずるという問題を有する。具体的には、波形データからピックアップするデータ相互の間の時間間隔Δｔおよび観測点数Ｎを用いて、
Δω＝１／（ＮΔｔ）・・・・・・・（１）
と表される。式（１）は、高速フーリエ変換によってスペクトル推定をおこなう限り、どのようなアルゴリズムに基づいても成立する。したがって、高速フーリエ変換によるスペクトル推定は、Δωだけ周波数に対して不確定性を有するという問題点を有する。
【０００５】
このような周波数誤差を生じない方法として、自己回帰モデルを用いたスペクトル推定の方法が知られている（S. Haykin. Ed, Nonlinear Methods of Spectral Analysis, Springer-Verlag, 1979 など）。自己回帰モデルとは、時刻ｔnにおけるデータｘnをｔnよりも前の時刻ｔn-2、ｔn-3、ｔn-4、・・・・におけるデータの線形結合で表現することをいう。すなわち、自己回帰モデルを用いたスペクトル推定とは、時刻ｔnにおけるデータｘnが、
ｘn＝ａ1ｘn-1＋ａ2ｘn-2＋ａ3ｘn-3＋・・・＋ａMｘn-M ・・・（２）
の形で表されることを基礎として、（２）式で表されたｘnについて、任意の数ｚを用いてｚ変換をおこない、周波数を変数とした信号強度の分布に関するデータを得ることによりスペクトル推定をおこなう方法である。
【０００６】
図８（ａ）は、ある一定の現象における１秒間の間の信号強度を示す時系列情報である。そして、図８（ｂ）は、（ａ）の時系列情報を高速フーリエ変換することにより得られたスペクトルである。これに対し、図８（ｃ）は、（ａ）の時系列情報を自己回帰モデルを用いたスペクトル推定により得られたスペクトルである。図８（ｂ）と（ｃ）を比較すると、高速フーリエ変換では全体的にピークの半値幅が大きく、３１．５Ｈｚのピークと、３２．５Ｈｚのピークが分離できていない。本来２つのピークからなる部分が１つのピークで表されているため、図８（ｂ）から得られる周波数が、実際のものと異なったものとなる。これに対して、図８（ｃ）は、ピークの半値幅も狭く、各ピークの分離も良好になされているため、周波数も正確な値を得ることができる。
【０００７】
【発明が解決しようとする課題】
一方で、自己回帰モデルを用いたスペクトル推定についても問題点は存在する。自己回帰モデルにおいては、（２）式の右辺の項数（以下「モデル次数」と言う）Ｍをどのように決定するかがスペクトル推定において重要となる。一般には自己回帰モデルで過去の値ｘn-iから現在の値ｘnを予測した時の予測誤差の２乗の期待値であるＦＰＥが最小になるＭを自己回帰モデルの次数とする方法が知られている。具体的には、ＦＰＥは以下の式で表される。
ＦＰＥ＝（１＋（Ｍ＋１）／Ｎ）（１−（Ｍ＋１）／Ｎ）-1（Σ（ｘn−Σａiｘn-i））・・・・（３）
ここで、Σａiｘn-iは、１≦ｉ≦Ｍについて和をとり、Σ（ｘn−Σａiｘn-i）は、０≦ｎ≦Ｎについて和をとるものとする。
【０００８】
しかし、この方法も十分ではない。なぜならば、ＦＰＥはあくまで統計的に最適なモデル次数Ｍを与える規範に過ぎず、これによって実際の観測データに対して必ず良好なスペクトルを導き出すものではないためである。したがって、経験上は、有用な周波数データを得るためのモデル次数Ｍは、データ点数Ｎの１／３〜１／２の値とする手法が知られており、データのＳ／Ｎ比が十分大きい場合にはこのような方法で決定した次数Ｍを用いることで図８（ｃ）に示す通りのスペクトルに関するグラフを得ることができる。
【０００９】
しかし、この経験則についても問題点が存在する。すなわち、十分なＳ／Ｎ比を得ることのできない波形データについては、上記の方法を適用しても正確なスペクトルを得ることはできない。特に、磁気共鳴映像の撮像においては、人体から生ずるノイズ等を無視することはできず、また、患者の負担を考慮して測定時間を長くとることができない。したがって、高いＳ／Ｎ比を有する波形データを得ることは困難であり、上記のモデル次数Ｍの決定方法をＭＲＩに用いることは適切ではない。
【００１０】
また、波形データのＳ／Ｎ比が小さい場合には、上記の問題点とは別の問題点も存在する。すなわち、自己回帰モデルを用いてスペクトル推定をおこなった場合には、有意な波形データだけでなく、ノイズを形成する波についてもスペクトル推定をおこなってしまう。そのため、得られたスペクトルは、波形データに起因する周波数のピークがノイズ波に起因するピーク中に埋もれてしまい、有用なデータとしての活用が難しくなる。
【００１１】
さらに、Ｓ／Ｎ比の小さい波形データからのスペクトル推定は、ＭＲＩ装置以外にも、電子回路や機械システムの特性解析など、さまざまなケースにおいて必要となる。これらの場合も、従来の手法によって自己回帰モデルを用いたスペクトル推定をおこなっても、有効なデータを得ることができない。
【００１２】
本発明は、上記問題点を解決するためになされたものであり、自己回帰モデルを用いたスペクトル推定において、Ｓ／Ｎ比が小さい波形データでも、分解能のスペクトル推定がおこなえる方法、この方法を用いた磁気共鳴映像撮像方法、磁気共鳴映像撮像装置、およびこの方法をコンピュータ上で実行するプログラムを提供することを目的とする。
【００１３】
【課題を解決するための手段】
上述した課題を解決し、目的を達成するため、第１の観点にかかる発明は、時間を変数とする波形データの所定の時刻における信号強度と、該所定の時刻までの離散時刻における所定数の信号強度とに基づいて周波数を変数とする信号強度波形を推定するスペクトル推定方法であって、観測データにノイズ除去フィルタを適用して前記波形データを導出する波形データ導出工程と、前記波形データから周波数を変数とする第１の信号強度波形を取得し、該取得した第１の信号強度波形の信号強度と雑音強度から前記所定数を決定する決定工程とを有することを特徴とする。
【００１４】
この第１の観点にかかる発明によれば、波形データの導出に関しノイズ除去フィルタを適用し、所定数の決定を信号強度と雑音強度から求めることとしているため、Ｓ／Ｎ比が低い観測データについても、分解能の高いスペクトル推定をおこなうことができる。
【００１５】
また、第２の観点にかかる発明は、第１の観点にかかる発明において、前記時間を変数とする波形データの所定の時刻における信号強度を前記所定の時刻までの離散時刻における所定数の信号強度の線形結合により表現した自己回帰モデルにより前記信号強度波形を推定する際に、前記線形結合における各項の係数を最小２乗法により決定する係数決定工程をさらに含んだことを特徴とする。
【００１６】
この第２の観点にかかる発明によれば、線形結合における各項の係数を最小２乗法により決定することとしたため、雑音を有する観測データに関しても、正確に各項の係数を決定することができる。
【００１７】
また、第３の観点にかかる発明は、第１または第２の観点にかかる発明において、前記波形データ導出工程は、−１と、前記観測データをフーリエ変換して得られた周波数を変数とする第２の信号強度波形における最大ピークの半値幅から導出される係数と、時間に関する変数との積を指数成分として持つ指数関数を含むノイズ除去フィルタ関数を前記観測データに乗算して前記波形データを導出することを特徴とする。
【００１８】
この第３の観点にかかる発明によれば、ノイズ除去フィルタ関数が指数関数を含み、その指数成分に、ピークの半値幅から導出される係数を有し、さらにノイズ除去フィルタ関数を観測データに乗算して波形データを導出することとしたため、観測データが有するノイズ成分を的確に除去し、導出される波形データを減衰振動関数とすることができる。
【００１９】
また、第４の観点にかかる発明は、第３の観点にかかる発明において、前記波形データ導出工程は、前記最大ピークの半値幅が０．５〜３．０Ｈｚとなる場合には該最大ピークの半値幅の値を前記係数の値とし、前記最大ピークの半値幅が０．５Ｈｚ以下の値となる場合には０．５を前記係数の値とし、前記最大ピークの半値幅が３．０Ｈｚ以上の値となる場合には３．０を前記係数の値としたことを特徴とする。
【００２０】
この第４の観点にかかる発明によれば、具体的にノイズ除去フィルタ関数の指数成分を決定したため、より的確にノイズ成分を除去することができる。
【００２１】
また、第５の観点にかかる発明は、第１〜第４のいずれかの観点にかかる発明において、前記決定工程は、前記第１の信号強度波形における信号強度の最大値と雑音強度値との比から前記所定数を決定することを特徴とする。
【００２２】
この第５の観点にかかる発明によれば、所定数を、前記スペクトルにおける信号強度の最大値と雑音強度値との比から求めることとしたため、ＦＰＥを用いた場合よりも高分解能のスペクトル推定をおこなうことができる。
【００２３】
また、第６の観点にかかる発明は、第１〜第５のいずれかの観点にかかる発明において、前記決定工程は、前記第１の信号強度波形における信号強度の最大値と雑音強度値との比の（−１／２）乗の値に比例定数を掛けた値の整数部分を前記所定数としたことを特徴とする。
【００２４】
この第６の観点にかかる発明によれば、所定数の決定をより具体的におこなうことで、さらに高分解能のスペクトル推定をおこなうことができる。
【００２５】
また、第７の観点にかかる発明は、撮像対象物を磁場雰囲気中に配置し、該配置した撮像対象物に対して前記磁場雰囲気の磁場強度に比例した中心周波数を有する電磁波を送信し、核磁気共鳴現象に基づいて前記撮像対象物から得られるＮＭＲ信号を受信し、受信した前記ＮＭＲ信号の波形データに基づいて自己回帰モデルによるスペクトルの推定をおこない、推定された前記スペクトルを位置に対する輝度情報に置き換えることにより撮像対象物の内部構造の映像を出力する磁気共鳴映像撮像方法において、前記スペクトルの推定は、観測データにノイズ除去フィルタを適用して前記波形データを導出する波形データ導出工程と、前記波形データから周波数を変数とする第１の信号強度波形を取得し、該取得した第１の信号強度波形の信号強度と雑音強度から前記所定数を決定する決定工程とを有することを特徴とする。
【００２６】
この第７の観点にかかる発明によれば、スペクトルの推定において、ノイズ除去フィルタの適用をおこない、信号強度と雑音強度から所定数を決定するため、Ｓ／Ｎ比の低い観測データから高分解能のスペクトルが推定され、そのスペクトルに基づいて高解像度の磁気共鳴映像を形成することができる。
【００２７】
また、第８の観点にかかる発明は、第７の観点にかかる発明において、前記波形データ導出工程は、前記観測データを高速フーリエ変換して得られた第２の信号強度波形において前記撮像対象物内部に存在する水分子を構成する水素原子に起因したピークの半値幅を指数成分の一部として持つ指数関数を含むノイズ除去フィルタ関数を前記観測データに乗算して前記波形データを導出することを特徴とする。
【００２８】
この第８の観点にかかる発明によれば、ノイズ除去フィルタ関数に含まれる指数関数における指数成分について、撮像対象物内部に存在する水分子を構成する水素原子に起因したピークの半値幅を含むこととしたため、有用なノイズ除去フィルタ関数を構成することができ、観測データにおける雑音成分を効果的に除去することができる。したがって、高い解像度を有する磁気共鳴映像を撮像することができる。
【００２９】
また、第９の観点にかかる発明は、第７または第８の観点にかかる発明において、前記決定工程は、前記第１の信号強度波形における信号強度の最大値と雑音強度値との比の（−１／２）乗の値に比例定数を掛けた値の整数部分を前記所定数としたことを特徴とする。
【００３０】
この第９の観点にかかる発明によれば、所定数の決定を具体的に規定することでさらにＦＰＥを用いた場合に比較してＳ／Ｎ比の低い観測データから高分解能のスペクトル推定をおこない、高解像度の磁気共鳴映像を撮像することができる。
【００３１】
また、第１０の観点にかかる発明は、第７〜第９のいずれかの観点にかかる発明において、前記スペクトル推定工程により得られた信号強度の分布から、複数のピーク間の周波数差の変動に基づいて前記撮像対象物の内部の温度変動を検知する検知工程をさらに含んだことを特徴とする。
【００３２】
この第１０の観点にかかる発明によれば、代謝物を構成する原子に起因するピークは、温度変化によって共鳴周波数も変動するため、その変動幅を測定することで撮像対象物内部の温度変動を検出することができる。
【００３３】
また、第１１の観点にかかる発明は、第１〜第１０の観点にかかる発明のいずれか一つの方法をコンピュータ上で実行するプログラムとしたことで、第１〜第１０のいずれかの観点における発明にかかる方法をコンピュータによって実行することができる。
【００３４】
また、第１２の観点にかかる発明は、均一な磁場を発生する静磁場発生磁石と、位置および時間に応じて変動する磁場を発生する勾配磁場発生部と、撮像対象物に電磁波を照射する送信部と、前記撮像対象物からのＮＭＲ信号を含む観測データを受信する受信部と、磁気共鳴画像を形成するための画像形成処理部とを有し、前記均一な磁場および前記変動する磁場の強度に比例した周波数の電磁波の送信およびＮＭＲ信号の受信により磁気共鳴映像を撮像する磁気共鳴映像撮像装置であって、前記画像形成処理部は、前記観測データとノイズ除去フィルタ関数とに基づいて時系列情報からなる波形データを導出し、導出された波形データを高速フーリエ変換して得られた第１の信号強度波形における信号強度と雑音強度から多項式の項数を決定しつつ自己回帰モデルによるスペクトル推定をおこない、周波数を変数とするスペクトルを推定することにより磁気共鳴映像を形成する磁気共鳴映像撮像装置であることを特徴とする。
【００３５】
この第１２の観点にかかる発明によれば、画像形成処理部が自己回帰モデルを用いたスペクトルの推定をおこない、スペクトルの推定において、ノイズ除去フィルタ関数を用いて波形データを導出し、信号強度と雑音強度によって多項式の項数を決定する磁気共鳴映像撮像装置としたため、Ｓ／Ｎ比の低い観測データであっても自己回帰モデルを用いて高分解能のスペクトル推定をおこない、高解像度の磁気共鳴映像を撮像可能な装置を提供することができる。
【００３６】
また、第１３の観点にかかる発明は、第１２の観点にかかる発明において、前記ノイズ除去フィルタ関数は、前記観測データを高速フーリエ変換して得られた周波数を変数とする第２の信号強度波形において前記撮像対象物内部に存在する水分子を構成する水素原子に起因したピークの半値幅を指数成分の一部として持つ指数関数を含んだことを特徴とする。
【００３７】
この第１３の観点にかかる発明によれば、ノイズ除去フィルタ関数を水素原子に起因したピークの半値幅を指数成分に含む指数関数を含むこととしたため、画像形成処理部は、観測データから効果的にノイズ成分を除去することができ、高解像度の磁気共鳴映像を撮影できる磁気共鳴映像撮像装置を提供できる。
【００３８】
また、第１４の観点にかかる発明は、第１２または第１３の観点にかかる発明において、前記多項式の項数は、前記第１の信号強度波形における前記信号強度と前記雑音強度との比を（−１／２）乗した値に比例定数を乗算した値の整数部分であることを特徴とする。
【００３９】
この第１４の観点にかかる発明によれば、画像形成処理部が、ＦＰＥを用いてモデル次数を決定した場合と比較して、Ｓ／Ｎ比の低い観測データに対してより適切にモデル次数を提供することができる磁気共鳴映像撮像装置を提供することができる。
【００４０】
また、第１５の観点にかかる発明は、第１２〜第１４のいずれかの観点にかかる発明において、前記画像形成処理部は、前記スペクトル推定により得られたスペクトルにおいて、複数のピーク間の周波数差の変動に基づいて撮像対象物の温度変動を検知することを特徴とする。
【００４１】
この第１５の観点にかかる発明によれば、処理手段が、複数のピーク間の周波数の差を検知することで撮像対象物の内部における温度変動を検知する磁気共鳴映像撮像装置としたため、直接測定することのできない物体内部の温度について、正確に検知することができる。
【００４２】
【発明の実施の形態】
実施の形態１．
以下、図面を参照して実施の形態１にかかるＭＲＩ装置について説明する。図１は、実施の形態１にかかるＭＲＩ装置の概要を示す模式図である。実施の形態１にかかるＭＲＩ装置は、静磁場発生部１と、勾配磁場を印加する勾配コイル部２と、電磁波の照射およびＮＭＲ信号を受信するＲＦコイル部３とを有する。また、勾配コイル部２は、勾配駆動部４によって制御され、ＲＦコイル部３は、送信部５および受信部６の制御を受ける。また、ＲＦコイル部３は、アナログ・ディジタル変換部７を介してコンピュータ８に接続されている。さらに、勾配駆動部４、送信部５、受信部６、アナログ・ディジタル変換部７はそれぞれ制御部９の制御を受け、また、コンピュータ８は操作部１０および表示部１１を有し、制御部９を制御する。
【００４３】
静磁場発生部１は、撮像対象物１２が配置された空間に時間的および空間的に均一な磁場を印加するためのものである。静磁場発生部１は、永久磁石または超電導磁石などからなり、１．０Ｔ程度の強磁場を発生する機能を有する。
【００４４】
勾配コイル部２は、勾配磁場を発生するためのものである。上述のように、ＭＲＩ装置は、撮像対象物１２を構成する特定の核種（主に水素原子）に対して核磁気共鳴現象を利用するが、位置情報を調べるために静磁場とは別に空間的および時間的に変動する勾配磁場を撮像対象物１２に対して印加する必要があるためである。
【００４５】
ＲＦコイル部３は、撮像対象物１２に対して電磁波を照射する電磁波送信手段としての機能と、核磁気共鳴現象によって撮像対象物１２を構成する原子から放出されるＮＭＲ信号を受信するＮＭＲ信号受信手段としての機能を有する。なお、ＲＦコイル部３によって照射される電磁波の周波数は、静磁場発生部１により印加される磁場の強度が１．０Ｔの場合には、数ＭＨｚ〜数十Ｍｚであることが望ましい。
【００４６】
勾配駆動部４は、勾配コイル部２を制御するためのものである。勾配コイル部２によって印加される勾配磁場は、その大きさが空間的および時間的に変動させる必要があるため、勾配駆動部４によって勾配コイル部２を制御して勾配磁場を変動させている。
【００４７】
送信部５、受信部６は、それぞれＲＦコイル部３を制御するためのものである。上述の通り、ＲＦコイル部３は、電磁波送信手段としての機能と、ＮＭＲ信号受信手段としての機能を有する。具体的には、送信部５から送られた電流に基づいてＲＦパルスを発生して電磁波を送信し、受信したＮＭＲ信号を受信部６に出力する。受信部６はＲＦコイル部３から受け取ったＮＭＲ信号をアナログ・ディジタル変換部７に伝える。
【００４８】
アナログ・ディジタル変換部７は、ＲＦコイル部３が受信したＮＭＲ信号に関するアナログデータをディジタルデータに変換するためのものである。撮像対象物１２から放出され、コンピュータ８上でスペクトル解析をするためには観測データをディジタルデータの形式に変換しておく必要があるためである。
【００４９】
コンピュータ８は、ＲＦコイル部３で受信され、アナログ・ディジタル変換部７でディジタルデータに変換された観測データから、スペクトル推定をおこなうためのものである。コンピュータ８は、スペクトル推定をおこない、スペクトル推定の結果より撮像対象物１２の内部構造を分析し、磁気共鳴映像を形成する機能を有する。
【００５０】
また、コンピュータ８は、制御部９を制御するためのものでもある。コンピュータ８に接続された操作部１０に入力された情報に基づいて、コンピュータ８は、磁気共鳴映像の撮像の開始や、スライス幅の決定等をおこない、制御部９に指令を送る機能を有する。
【００５１】
表示部１１は、コンピュータ８で形成された磁気共鳴映像を表示するためのものである。表示部１１の構造としては、ＣＲＴディスプレイや、ＴＦＴ液晶画面のように直接映像を表示するものでもよいし、プリンタのように映像を他の媒体に出力するものでもよい。
【００５２】
次に、実施の形態１にかかるＭＲＩ装置を用いて、実際に磁気共鳴映像を撮像する方法について、説明する。実施の形態１における磁気共鳴映像の撮像は、大きく分けて磁場の印加および電磁波の照射、受信による観測データの取得工程と、取得された観測データをスペクトルに変換し、磁気共鳴映像を形成する工程とからなる。
【００５３】
まず、観測データの取得工程について、図２を参照して、説明する。最初に、撮像対象物１２をＭＲＩ装置の内部に搬入し、静磁場発生部１による発生した静磁場雰囲気中に撮像対象物１２を配置する（ステップＳ１０１）。なお、静磁場発生部１が超電導磁石からなる場合は、あらかじめ超電導磁石を液体ヘリウムなどで冷却し、超伝導状態を作り出しておく必要がある。
【００５４】
そして、勾配コイル部２がスライス方向に勾配磁場を印加するとともに、ＲＦコイル部３は撮像対象物１２に対して電磁波を照射する（ステップＳ１０２）。スライス方向に勾配磁場を印加することにより、撮像対象物１２は、共鳴周波数が異なり、法線がスライス方向と一致する多数のスライスに分割される。そして、任意のスライスについての共鳴周波数と同一周波数の高周波パルスからなる電磁波を撮像対象物１２に照射することで、所望のスライスに含まれる原子のスピンのみを励起することができ、原子の位置に関して、スライス方向の座標について特定することができる。
【００５５】
そして、勾配コイル部２は、撮像対象物１２に対して位相勾配磁場を印加する（ステップＳ１０３）。位相勾配磁場を印加することにより、ステップＳ１０２において選択されたスライス中に存在する原子から放出される電磁波の位相が、位相勾配磁場の分布にしたがって異なるものとなる。具体的には、スライスは、位相方向を短手方向とした複数の短冊に分割され、異なる短冊に属する原子からは、異なる位相の電磁波が放出される。したがって、本ステップ終了時点において電磁波を放出する原子の位置は、スライス方向の座標および位相方向の座標について特定できることとなる。
【００５６】
そして、勾配コイル部２によって周波数勾配磁場を印加しながら、撮像対象物１２を構成する原子から放出されるＮＭＲ信号を、ＲＦコイル部３によって受信する（ステップＳ１０４）。周波数方向に勾配磁場を印加することにより、ＲＦコイル部３は、周波数方向に沿って異なる周波数のＮＭＲ信号を受信することができる。したがって、ＮＭＲ信号を放出する原子に関して、周波数方向の座標を特定することができ、ステップＳ１０２〜Ｓ１０３とあわせて、ＮＭＲ信号を放出する原子の位置が完全に特定される。そして、ステップＳ１０２〜Ｓ１０４の走査を位相方向のマトリックス数だけ繰り返し実施し、２次元のデータセットを得る。
【００５７】
以上ステップＳ１０１からステップＳ１０４の工程を経ることによって、撮像対象物１２を構成する原子がどの位置に、どれだけ存在するかを示す時系列情報である観測データを取得することができる。なお、これらの工程における勾配磁場の印加および電磁波の照射、ＮＭＲ信号の受信はそれぞれ直接には勾配駆動部４、送信部５、受信部６によっておこなわれるが、これらの動作は相互に連携させておこなう必要があるため、制御部９は勾配駆動部４、送信部５、受信部６にそれぞれ指令を出し、適切な同期をとっている。
【００５８】
なお、上記観測データの取得工程では、撮像対象物１２の内部器官を構成する代謝物と、撮像対象物１２の内部に含まれる水についておこなわれる。元来、撮像対象物１２の内部構造を知るためには内部器官についての情報のみがあれば足りるが、本実施の形態１においては、後のスペクトル推定において水に関する情報についても必要となるためである。具体的には、次のようにして代謝物および水についての観測データを取得する。すなわち、内部器官を構成する代謝物の分子に含まれる水素原子と、水分子に含まれる水素原子とでは、分子構造の相違から共鳴周波数が異なる。したがって、ＲＦコイル部３から照射される電磁波の周波数を変化させることで、内部器官に関する観測データと、水に関する観測データを取得することができる。
【００５９】
次に、上記の工程で取得された観測データをスペクトルに変換し、磁気共鳴映像を形成する工程について、図３〜図６を参照して説明する。なお、図４〜図６に示した図は、本工程の内容理解の補助を目的とした模式的な図であり、実際の観測データなどとは必ずしも一致しない。
【００６０】
まず、水に関する観測データによって、代謝物の観測データの位相を補正する（ステップＳ２０１）。観測データの取得においては、時間に応じて変化する勾配磁場を印加することによりおこなうため撮像対象物１２には誘導起電力が生じ、渦電流が発生する。この渦電流によって観測データの値に乱れが生じ、この乱れを補正するために、本ステップはおこなわれる。
【００６１】
そして、水に関する観測データを高速フーリエ変換する（ステップＳ２０２）。この際の高速フーリエ変換は従来の手法によっておこなうことが可能で、本ステップの結果、水に関するスペクトルを得ることができる。その後、ステップＳ２０２で得られたスペクトルについて、水に関するピークの半値幅を測定する（ステップＳ２０３）。
【００６２】
そして、ステップＳ２０３で測定した水に関するピークの半値幅からノイズ除去フィルタ関数を決定する（ステップＳ２０４）。具体的には、ノイズ除去フィルタ関数を指数関数とし、指数成分を−１と、ピークの半値幅と、時間変数ｔとの積からなるとする。ここで、半値幅の値が０．５Ｈｚよりも小さい場合には、ピークの半値幅にかえて０．５を代入し、半値幅の値が３．０Ｈｚよりも大きい場合には、ピークの半値幅にかえて３．０を代入する。
【００６３】
その後、代謝物に関する観測データに対して、ステップＳ２０４で決定したノイズ除去フィルタ関数を適応する（ステップＳ２０５）。具体的には、図４（ａ）に示すように、時間を変数とする観測データｌ1と、同じく時間を変数とするノイズ除去フィルタ関数ｌ2とを掛け合わせ、新たに図４（ｂ）に示す波形データｌ3を得る。この波形データｌ3は、ノイズ除去フィルタ関数を掛け合わせることにより、減衰振動波形となっていることが、図４（ｂ）より明らかである。本ステップにより、ノイズ成分を除去した波形データが得られる。
【００６４】
そして、ステップＳ２０５で得られた波形データに対して、高速フーリエ変換を行う（ステップＳ２０６）。高速フーリエ変換をおこなうことで、図５に示すような波形データのスペクトルを得ることができる。
【００６５】
その後、波形データのスペクトルにおいて、図５に示すピーク強度の最大値ＰSと、ノイズ成分の大きさＰNを測定し、強度比であるＳ／Ｎ比を導出する（ステップＳ２０７）。そして、Ｓ／Ｎ比に基づいてモデル次数Ｍを決定する（ステップＳ２０８）。ここで、Ｓ／Ｎ比が２５０より小さい場合は、以後Ｓ／Ｎ比を２５０として計算し、Ｓ／Ｎ比が４００よりも大きい場合は、以後Ｓ／Ｎ比を４００として計算をおこなう。
【００６６】
具体的には、Ｓ／Ｎ比の−１／２乗の値に対して比例定数αをかけた値をモデル次数Ｍとする。なお、比例定数αは、測定データのデータポイント数と、データ取得の時間間隔によって決定される値であり、たとえば、データポイント数が２０４８の場合には、比例定数αの値は５０となる。
【００６７】
そして、ステップＳ２０８で決定したモデル次数Ｍに基づいて、Ｍ個の項数を有する多項式を作成する（ステップＳ２０９）。具体的には、図６で示すように、相互の間隔がΔｔの離散時刻ｔ0、ｔ1、ｔ2、ｔ3、・・・・・、ｔM、・・・・・、ｔNについて、強度ｘ0、ｘ1、ｘ2、ｘ3、・・・・・、ｘM、・・・・、ｘNを波形データから読みとり、これらについて以下の多項式を作成する。
ｘM+i＝ａ1ｘM+i-1＋ａ2ｘM+i-2＋ａ3ｘM+i-3＋・・・・・＋ａM-1ｘM+i-(M-1)＋ａMｘM+i-M・・・・・（４）
ここで、ｉ＝０、１、２、・・・・・、Ｎであるとし、Ｎ＝２Ｍとする。したがって（４）式はＭ個の多項式からなる。また、ｘjは波形データから読みとられた具体的な数値であることから、（４）式における未知数はａ1、ａ2、・・・・・、ａMである。したがって（４）式はＭ個の未知数に対するＭ個の方程式であり、これらの方程式を解くことで自己回帰モデルにおける係数を導出する。
【００６８】
その後、（４）式のうちｉ＝０の式について、ｚ変換を行う（ステップＳ２１０）。（４）式のうちｉ＝０の式をｚ変換したＸ（ｚ）は、次の式であらわされる。
Ｘ（ｚ）＝（１−Σａiｚ-1）-1 ・・・・・（５）
ここで、Σａiｚ-1は、ｉについて１からＭまでの和をとるものとする。
【００６９】
そして、（５）式に対して変数変換をおこなう（ステップＳ２１１）。具体的には、変数ｚを、（ｉｆΔｔ／２π）に置き換える。ここで、ｉは虚数単位であり、ｆは周波数であり、Δｔは波形データを読みとった時刻の間隔である。具体的には、変数変換の結果、（６）式は、次の式に変形される。
Ｘ（ｆ）＝（１−Σａi（ｉｆΔｔ／２π）-1）-1 ・・・（６）
ここで、ａiは（５）式から具体的に求められており、Δｔ、２π、ｉは定数であることから、（７）式は、周波数ｆ（Ｈｚ）を変数とする関数であることが分かる。
【００７０】
最後に、（６）式の絶対値を求める（ステップＳ２１２）。（６）式は虚数ｉを含む式であるため、スペクトルを知るためには（７）式を自然数からなる関数に変換する必要があるためである。具体的には、（６）式の右辺に共役複素数を掛け合わせ、その積の平方根を求める。以上で、取得された観測データをスペクトルに変換する工程が終了する。
【００７１】
本工程によって得られたスペクトルについて、図７（ｂ）にグラフを示す。なお、図７（ｂ）においては、ノイズ除去フィルタ関数の指数成分であるピーク半値幅を１．８Ｈｚ、モデル次数Ｍを３００としてスペクトル推定をおこなっている。同一の観測データをもとに高速フーリエ変換によって得られたスペクトルである図７（ａ）と比較して、スペクトルのピークの半値幅が明らかに狭い。また、副次的な効果としてノイズ除去フィルタ関数を観測データに掛け合わせたことにより、結果のスペクトルにおけるノイズ成分の大きさも図７（ｂ）のグラフのほうが明らかに小さくなっていることが分かる。
【００７２】
その後、上記の工程で得られたスペクトルに基づいて各スライス面に対する２次元の磁気共鳴画像がコンピュータ上で作成され、表示部１１に出力されることにより、本実施の形態１にかかるＭＲＩ装置を用いた磁気共鳴映像の撮像が終了する。
【００７３】
ステップＳ２０４において、決定されるノイズ除去フィルタ関数を指数関数としたのは、次の理由である。画像データ取得工程におけるＲＦコイル部３による電磁波の照射は短時間におこなわれるため、これに対応して撮像対象物１２を構成する原子から放出されるＮＭＲ信号の強度は、ＲＦコイル部３による電磁波照射終了のあとには指数関数的に減少し、減衰振動波形を形成するはずである。しかし、観測データでは、図４の直線ｌ1に示すように、必ずしも指数関数的に減少しない。これは、ノイズ成分の強度が大きいために、核磁気共鳴現象に基づいて放出されたＮＭＲ信号の減少分がキャンセルされてしまうためである。そのため、指数減少関数を観測データに掛け合わせることにより、観測データからノイズ成分が除去され、核磁気共鳴現象に基づくＮＭＲ信号成分のみを波形データとして抽出することができる。
【００７４】
また、ステップＳ２０４におけるノイズ除去フィルタ関数の指数成分の決定について、水に関するピークの半値幅に比例するものとしたのは、次の理由による。すなわち、各スペクトルの半値幅は、時系列データである観測データの減衰振動成分に対応しており、スペクトルにおけるピークの半値幅を用いることで効果的にノイズ成分を除去できる。また、半値幅を測定するピークは、ピーク強度が大きいものを用いることが望ましく、水分子を構成する水素原子に起因するピークを選択している。
【００７５】
このように、本実施の形態１にかかるＭＲＩ装置は、従来の高速フーリエ変換を用いた磁気共鳴映像の撮像と比較して、スペクトル推定を正確におこなうことが可能である。その結果、高い解像度を有する磁気共鳴映像を撮像することが可能である。
【００７６】
また、本実施の形態１におけるスペクトル推定により、Ｓ／Ｎ比の小さい観測データに関しても自己回帰モデルを用いて分解能の高いスペクトル推定をおこなうことができる。したがって、ＭＲＩ装置による磁気共鳴映像の撮像以外にも、たとえば、気象現象の解析、電子回路や機械システムの特性解析などにおいて、Ｓ／Ｎ比の小さい時系列情報から、分解能の高いスペクトルを得ることができる。
【００７７】
なお、観測データの取得工程において、ステップＳ１０２、Ｓ１０３、Ｓ１０４は、上述したように、撮像対象物１２中の原子の位置を特定するためのステップである。これに対して、水に関する観測データは、あくまで代謝物に関する観測データをスペクトル変換する際に用いられるためのものにすぎないことから、水に関する観測データを取得する工程においては、勾配磁場の印加を省略してもよい。同様の理由から、水に関する観測データの取得に際しては、２次元のデータセットを得る必要がなく、ステップＳ１０２〜Ｓ１０４を繰り返しおこなうことを省略してもよい。
【００７８】
また、波形データは、ノイズ除去フィルタ関数によってノイズを低減しているものの、ノイズを完全に除去したものではない。したがって、ステップＳ２０９において、厳密には（４）式について左辺と右辺との間の等号は成立しない場合がある。したがって、より正確には最小２乗法によってａ1、ａ2、・・・・・、ａMを導出することが望ましい。具体的には、まず、ｉ＝０、１、２、・・・・、Ｎ−Ｍに関する（４）式を行列式で表現する。ここで、（ｘM、ｘM+1、・・・・、ｘN）からなる列ベクトルをｄと定義し、（ａ1、ａ2、ａ3、・・・・、ａM）からなる列ベクトルをａとする。また、Ｍkl＝ｘM-1+k-lを要素とするＭ×Ｎ行列をＨとする。これらｄ、ａ、Ｈにより（４）式は、
ｄ＝Ｈ・ａ ・・・・・・（７）
の形で表現することができる。したがって、最小２乗法によりａは、
ａ＝ＨtＨ-1Ｈtｄ ・・・・・（８）
と表現することができ、この（８）式をとくことでａが具体的に求められる。最小２乗法を用いることにより、より正確にａ1、ａ2、・・・・、ａMを求めることができる。また、これ以外にも、Ｙｕｌｅ−Ｗａｌｋｅｒ法や、Ｂｕｒｇ法によってａ1、ａ2、ａ3、・・・・、ａMの決定をおこなうことも有効である。
【００７９】
また、図１においてＲＦコイル部３は、撮像対象物１２である患者の頭部のみを覆う構造となっているが、本実施の形態１にかかるＭＲＩ装置はこの構造に限定されない。たとえば、患者の腹部を撮像するためにＲＦコイル部３の形状を変更してもよいし、患者の全身を撮像可能とするような形状としてもよい。さらに、図１に示すＭＲＩ装置は、水平方向に整磁場を印加する、いわゆる水平磁場型のＭＲＩ装置であるが、これ以外にも、たとえば、垂直方向に静磁場を印加する、いわゆる垂直磁場型の構造からなるＭＲＩ装置であってもよい。
【００８０】
さらに、本実施の形態１においては、図１に示したように、ＭＲＩの構造に関して、勾配駆動部４、送信部５など、装置に対して制御をおこなう部分について、それぞれ別個独立のものとしている。これは、本実施の形態１にかかるＭＲＩ装置の構造の説明を容易にする目的によるものであって、必ず各部分が分離した構造となることを示すものではない。したがって、たとえば、制御部９をコンピュータ８に内蔵する構造としてもよいし、アナログ・ディジタル変換部７をコンピュータ８に内蔵して、アナログデータを直接コンピュータ８に入力する構造としてもよい。
【００８１】
また、本実施の形態１における観測データの取得工程に関してスライス方向は、通常は撮像対象物１２である患者の頭部からつま先に向かう方向とするが、これに限定されるものではない。撮像する磁気共鳴映像の領域にあわせてスライス方向は変化させることが可能である。また、周波数方向、位相方向、スライス方向は直交座標系を形成するのが通常であるが、本実施の形態１における各方向は必ずしも直交座標系に限定されるものではない。
【００８２】
また、本実施の形態１におけるスペクトル推定について、操作部１０を通じてオペレータが手動でスペクトル推定をおこなってもよいが、より望ましい形態としては上記のスペクトル推定をおこなうプログラムをコンピュータ８上で実行することが望ましい。プログラムをコンピュータ８上で実行する形態とすることで、自動的にスペクトル推定をおこなうことができる。なお、スペクトル推定のみをプログラムで実行してもよいが、ＭＲＩ装置の制御、並びに推定されたスペクトルから磁気共鳴映像を形成する工程についてもプログラムで実行させる形態としてもよい。
【００８３】
実施の形態２．
次に、実施の形態２にかかるＭＲＩ装置について、説明する。実施の形態２にかかるＭＲＩ装置は、磁気共鳴映像を撮像する機能と共に、撮像対象物の内部の温度変化を検出する機能を有する。実施の形態２にかかるＭＲＩ装置の構造は、図１に示す実施の形態１にかかるＭＲＩ装置の構造と同様である。
【００８４】
また、本実施の形態２における観測データの取得工程およびスペクトルの推定は実施の形態１における工程とほぼ同様におこなわれるものとする。そして、これらの工程によって得られたスペクトルから、撮像対象物１２の内部における温度変動を観測する。
【００８５】
具体的には、実施の形態１と同様の方法でスペクトル推定をおこなった後、コンピュータ８上において、水を構成する水素原子に起因するピークの周波数と、代謝物を構成する水素原子に起因するピークの周波数との差について、推定されたスペクトルからの読みとりをおこなう。
【００８６】
そして、この読みとった周波数差を基準となるデータと比較することで、撮像対象物１２の内部温度がどのように変化したのか検知する。たとえば、あらかじめ３６度における周波数差に関するデータをあらかじめ取得しておき、本実施の形態２において読みとった周波数差を既知のデータと対比することにより、内部温度が３６度からどれだけ変動したのかを調べる。
【００８７】
一般に、代謝物の共鳴周波数は温度に依存して変動することが知られている。一方、水分子を構成する水素原子の共鳴周波数は温度に依存しないため、代謝物を構成する水素原子の共鳴周波数と、水分子を構成する水素原子の共鳴周波数との差を求めることにより、温度変動を検知することが可能である。代謝物を構成する水素原子の共鳴周波数の変動はわずかな値であったため、従来の高速フーリエ変換によるスペクトル推定では検出が困難であった。しかし、本実施の形態１におけるスペクトル推定では、精度よくスペクトル推定をおこなうことが可能となることから、温度変化に対する共鳴周波数の変動を十分検出することができる。すなわち、本実施の形態２にかかるＭＲＩ装置は、直接に温度測定をおこなうことのできない物体の内部について、温度変動を検知することができる。
【００８８】
このように、撮像対象物１２の内部の温度変動を測定することで、次の利点が生ずる。たとえば、本実施の形態２にかかるＭＲＩ装置を人間の脳に関する診断装置として利用することができる。
【００８９】
脳挫傷等の脳疾患において、脳の温度を知ることはきわめて重要である。すなわち、脳の温度が上昇した場合には、脳浮腫を併発して脳の体積が増大し、患者が死亡する場合がある。このような症例において、本実施の形態２にかかるＭＲＩ装置で脳の温度の上昇を検知できれば、適切な治療、たとえば低体温療法を施すことによって患者を治療することができる。
【００９０】
なお、本実施の形態２においては、磁気共鳴画像の撮像そのものは目的としていないため、観測データの取得工程において、勾配磁場を印加して撮像する工程を省略し、静磁場のみを印加して、ＲＦコイル部３による電磁波の送信をおこなってもよい。撮像対象物１２の内部の温度変動を検知するためには、代謝物の共鳴周波数と、水の共鳴周波数が分かれば十分なためである。逆に、観測データの取得工程を実施の形態１の場合と同様におこない、温度変動の測定と共に撮像対象物１２の磁気共鳴映像を撮像する構成としてもよい。
【００９１】
なお、本実施の形態１および２においては、自己回帰モデルを用いたスペクトル推定方法について説明したが、自己回帰モデルを用いたスペクトル推定方法は、いわゆるＭＥＭ（Maximum Entropy Method）と同一のスペクトル推定方法であることが証明されている(S. L. Marple, Digital Spectral Analysis, Prentice-Hall, 1987)。したがって、本発明は、このＭＥＭについても適用できることは当然である。
【００９２】
【発明の効果】
上述してきたように、第１の観点にかかる発明によれば、波形データの導出に関しノイズ除去フィルタを適用し、所定数の決定を信号強度と雑音強度から求めることとしているため、Ｓ／Ｎ比が低い観測データについても、分解能の高いスペクトル推定をおこなうことができるという効果を奏する。
【００９３】
また、第２の観点にかかる発明によれば、線形結合における各項の係数を最小２乗法により決定することとしたため、雑音を有する観測データに関しても、正確に各項の係数を決定することができるという効果を奏する。
【００９４】
また、第３の観点にかかる発明によれば、ノイズ除去フィルタ関数が指数関数を含み、その指数成分に、ピークの半値幅から導出される係数を有し、さらにノイズ除去フィルタ関数を観測データに乗算して波形データを導出することとしたため、観測データが有するノイズ成分を的確に除去し、導出される波形データを減衰振動関数とすることができるという効果を奏する。
【００９５】
また、第４の観点にかかる発明によれば、具体的にノイズ除去フィルタ関数の指数成分を決定したため、より的確にノイズ成分を除去することができるという効果を奏する。
【００９６】
また、第５の観点にかかる発明によれば、所定数を、前記スペクトルにおける信号強度の最大値と雑音強度値との比から求めることとしたため、ＦＰＥを用いた場合よりも高分解能のスペクトル推定をおこなうことができるという効果を奏する。
【００９７】
また、第６の観点にかかる発明によれば、所定の項数の決定をより具体的におこなう構成としたことで、さらに高分解能のスペクトル推定をおこなうことができるという効果を奏する。
【００９８】
また、第７の観点にかかる発明によれば、スペクトルの推定において、ノイズ除去フィルタの適用をおこない、信号強度と雑音強度から所定数を決定するため、Ｓ／Ｎ比の低い観測データから高分解能のスペクトルが推定され、そのスペクトルに基づいて高解像度の磁気共鳴映像を形成することができるという効果を奏する。
【００９９】
また、第８の観点にかかる発明によれば、ノイズ除去フィルタ関数に含まれる指数関数における指数成分について、撮像対象物内部に存在する水分子を構成する水素原子に起因したピークの半値幅を含むこととしたため、有用なノイズ除去フィルタ関数を構成することができ、観測データにおける雑音成分を効果的に除去することができる。したがって、高い解像度を有する磁気共鳴映像を撮像することができるという効果を奏する。
【０１００】
また、第９の観点にかかる発明によれば、所定数の決定を具体的に規定することでさらにＦＰＥを用いた場合に比較してＳ／Ｎ比の低い観測データから高分解能のスペクトル推定をおこない、高解像度の磁気共鳴映像を撮像することができるという効果を奏する。
【０１０１】
また、第１０の観点にかかる発明によれば、代謝物を構成する原子に起因するピークは、温度変化によって共鳴周波数も変動するため、その変動幅を測定することで撮像対象物内部の温度変動を検出することができるという効果を奏する。
【０１０２】
また、第１１の観点にかかる発明によれば、第１〜第１０の観点にかかる発明のいずれか一つの方法をコンピュータ上で実行するプログラムとしたことで、第１〜第１０のいずれかの観点における発明にかかる方法をコンピュータによって実行することができるという効果を奏する。
【０１０３】
また、第１２の観点にかかる発明によれば、画像形成処理部が自己回帰モデルを用いたスペクトルの推定をおこない、スペクトルの推定において、ノイズ除去フィルタ関数を用いて波形データを導出し、信号強度と雑音強度によって多項式の項数を決定する磁気共鳴映像撮像装置としたため、Ｓ／Ｎ比の低い観測データであっても自己回帰モデルを用いて高分解能のスペクトル推定をおこない、高解像度の磁気共鳴映像を撮像可能な装置を提供することができるという効果を奏する。
【０１０４】
また、第１３の観点にかかる発明によれば、ノイズ除去フィルタ関数を水素原子に起因したピークの半値幅を指数成分に含む指数関数を含むこととしたため、画像形成処理部は、観測データから効果的にノイズ成分を除去することができ、高解像度の磁気共鳴映像を撮影できる磁気共鳴映像撮像装置を提供できるという効果を奏する。
【０１０５】
また、第１４の観点にかかる発明によれば、画像形成処理部が、ＦＰＥを用いてモデル次数を決定した場合と比較して、Ｓ／Ｎ比の低い観測データに対してより適切にモデル次数を提供することができる磁気共鳴映像撮像装置を提供することができるという効果を奏する。
【０１０６】
また、第１５の観点にかかる発明によれば、処理手段が、複数のピーク間の周波数の差を検知することで撮像対象物の内部における温度変動を検知する磁気共鳴映像撮像装置としたため、直接測定することのできない物体内部の温度について、正確に検知することができるという効果を奏する。
【図面の簡単な説明】
【図１】実施の形態１にかかるＭＲＩ装置の構造を示す模式図である。
【図２】実施の形態１における観測データ取得の工程を示すフローチャートである。
【図３】実施の形態１におけるスペクトル推定の工程を示すフローチャートである。
【図４】（ａ）は、実施の形態１におけるスペクトル推定において、ノイズ除去フィルタ関数を観測データに掛け合わせる模式図である。（ｂ）は、掛け合わせた結果として導出された波形データを示す模式図である。
【図５】実施の形態１におけるスペクトル推定において、モデル次数Ｍを決定する際のピーク強度およびノイズ強度を示す模式図である。
【図６】実施の形態１におけるスペクトル推定において、観測データから自己回帰モデルによる多項式を導出する際の模式図である。
【図７】（ａ）は、実施の形態１にかかるＭＲＩ装置で得られた観測データを高速フーリエ変換することにより、得られたスペクトル図である。（ｂ）は、同一の観測データを本実施の形態１におけるスペクトル推定方法をもちいて導出したスペクトル図である。
【図８】（ａ）は、時系列情報である観測データを示すグラフ図であり、（ｂ）は（ａ）のデータを高速フーリエ変換することによって得られたスペクトル図である。また、（ｃ）は、（ａ）のデータから自己回帰モデルを用いてスペクトル推定をおこなったことによるスペクトル図である。
【符号の説明】
１ 静磁場発生部
２ 勾配コイル部
３ ＲＦコイル部
４ 勾配駆動部
５ 送信部
６ 受信部
７ アナログ・ディジタル変換部
８ コンピュータ
９ 制御部
１０ 操作部
１１ 表示部[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a spectrum estimation method using a so-called autoregressive model, and in particular, a spectrum estimation method for accurately estimating spectrum with respect to observation data having a low S / N ratio, a magnetic resonance imaging method using the spectrum estimation method, and a spectrum The present invention relates to a magnetic resonance imaging apparatus using an estimation method, and a program for executing these methods on a computer.
[0002]
[Prior art]
Conventionally, an imaging object is placed in a magnetic field atmosphere and captured as waveform data by observing an NMR signal emitted from the imaging object due to a nuclear magnetic resonance phenomenon, and the waveform data is used as a variable. 2. Description of the Related Art A magnetic resonance imaging (hereinafter referred to as “MRI”) apparatus that converts data related to a spectrum that is a signal intensity distribution and images an internal structure of an imaging target is known. In order to improve the resolution of the image captured by the MRI apparatus, it is necessary to specify the peak frequency with high accuracy in the process of converting the waveform data composed of NMR signals into a spectrum.
[0003]
For the waveform data obtained by observation, a fast Fourier transform is generally used as a technique performed to know what spectrum the waveform data has. Fast Fourier transform is N times t discretized at an interval Δt in the observed waveform data. 1 , T 2 , T Three , ..., t i , ..., t N Intensity x (t i ), And an algorithm based on performing a discrete Fourier transform. When the spectrum is estimated using the fast Fourier transform, the frequency distribution of the observed wave can be easily obtained from the observed waveform data. If the waveform data has periodicity, the spectrum has a peak at a certain frequency. Appears, and the frequency of each wave constituting the waveform data can be obtained.
[0004]
However, when spectrum estimation is performed by Fast Fourier Transform, there is a problem that an error occurs over a certain range Δω because there is a limit to the accuracy of estimating the peak frequency. Specifically, using the time interval Δt between the data picked up from the waveform data and the number N of observation points,
Δω = 1 / (NΔt) (1)
It is expressed. Expression (1) is established based on any algorithm as long as spectrum estimation is performed by fast Fourier transform. Therefore, spectrum estimation by fast Fourier transform has a problem that it has uncertainty with respect to frequency by Δω.
[0005]
As a method that does not cause such a frequency error, a spectrum estimation method using an autoregressive model is known (S. Haykin. Ed, Nonlinear Methods of Spectral Analysis, Springer-Verlag, 1979, etc.). Autoregressive model is time t n Data x in n T n T before n-2 , T n-3 , T n-4 It is expressed by linear combination of data in. That is, spectrum estimation using the autoregressive model is the time t n Data x in n But,
x n = A 1 x n-1 + A 2 x n-2 + A Three x n-3 + ... + a M x nM ... (2)
X expressed by the formula (2) on the basis of being expressed in the form of n The spectrum estimation is performed by performing z-transform using an arbitrary number z and obtaining data on the signal intensity distribution with the frequency as a variable.
[0006]
FIG. 8A shows time-series information indicating the signal intensity for one second in a certain phenomenon. FIG. 8B shows a spectrum obtained by performing a fast Fourier transform on the time series information of FIG. On the other hand, FIG. 8C shows a spectrum obtained by estimating the time series information of FIG. 8A using the autoregressive model. Comparing FIGS. 8B and 8C, the full width at half maximum of the peak is large as a whole in the fast Fourier transform, and the peak at 31.5 Hz and the peak at 32.5 Hz cannot be separated. Since the part consisting of two peaks is represented by one peak, the frequency obtained from FIG. 8B is different from the actual one. On the other hand, in FIG. 8C, the half-value width of the peak is narrow and each peak is well separated, so that an accurate value can be obtained for the frequency.
[0007]
[Problems to be solved by the invention]
On the other hand, there is a problem with spectrum estimation using an autoregressive model. In the autoregressive model, how to determine the number of terms (hereinafter referred to as “model order”) M on the right side of the equation (2) is important in spectrum estimation. In general, past value x in autoregressive model ni To current value x n There is known a method in which M that minimizes FPE, which is the expected value of the square of the prediction error when predicting, is the order of the autoregressive model. Specifically, FPE is represented by the following formula.
FPE = (1+ (M + 1) / N) (1- (M + 1) / N) -1 (Σ (x n -Σa i x ni )) ... (3)
Where Σa i x ni Is the sum of 1 ≦ i ≦ M and Σ (x n -Σa i x ni ) Shall take the sum for 0 ≦ n ≦ N.
[0008]
However, this method is not sufficient. This is because FPE is merely a norm that gives a statistically optimal model order M, and does not necessarily lead to a good spectrum for actual observation data. Therefore, from experience, it is known that the model order M for obtaining useful frequency data is a value that is 1/3 to 1/2 of the number N of data points, and the S / N ratio of the data is sufficiently large. In some cases, the order relating to the spectrum as shown in FIG. 8C can be obtained by using the order M determined by such a method.
[0009]
However, there are problems with this rule of thumb. That is, for waveform data for which a sufficient S / N ratio cannot be obtained, an accurate spectrum cannot be obtained even if the above method is applied. In particular, in magnetic resonance imaging, noise generated from the human body cannot be ignored, and the measurement time cannot be increased in consideration of the burden on the patient. Therefore, it is difficult to obtain waveform data having a high S / N ratio, and it is not appropriate to use the above-described method for determining the model order M for MRI.
[0010]
Further, when the S / N ratio of the waveform data is small, there is another problem different from the above problems. That is, when spectrum estimation is performed using an autoregressive model, spectrum estimation is performed not only on significant waveform data but also on waves forming noise. For this reason, in the obtained spectrum, the peak of the frequency due to the waveform data is buried in the peak due to the noise wave, making it difficult to utilize it as useful data.
[0011]
Furthermore, spectrum estimation from waveform data with a small S / N ratio is required in various cases such as characteristic analysis of electronic circuits and mechanical systems in addition to MRI apparatuses. Even in these cases, effective data cannot be obtained even if spectrum estimation using an autoregressive model is performed by a conventional method.
[0012]
The present invention has been made in order to solve the above-described problems. In spectrum estimation using an autoregressive model, a method capable of estimating resolution spectrum even with waveform data having a small S / N ratio, and this method are used. It is an object of the present invention to provide a magnetic resonance imaging method, a magnetic resonance imaging apparatus, and a program for executing the method on a computer.
[0013]
[Means for Solving the Problems]
In order to solve the above-described problems and achieve the object, the invention according to the first aspect is directed to a signal intensity at a predetermined time of waveform data having time as a variable, and a predetermined number of discrete data up to the predetermined time. A spectrum estimation method for estimating a signal intensity waveform having a frequency as a variable based on a signal intensity, wherein a waveform data deriving step of deriving the waveform data by applying a noise removal filter to observation data, and the waveform data And a determination step of acquiring a first signal intensity waveform having a frequency as a variable and determining the predetermined number from the signal intensity and noise intensity of the acquired first signal intensity waveform.
[0014]
According to the first aspect of the invention, a noise removal filter is applied to derive waveform data, and a predetermined number of decisions are obtained from signal intensity and noise intensity. Also, it is possible to perform spectrum estimation with high resolution.
[0015]
The invention according to the second aspect is the invention according to the first aspect, wherein the signal intensity at a predetermined time of the waveform data having the time as a variable is a predetermined number of signal intensities at discrete times up to the predetermined time. When the signal intensity waveform is estimated by the autoregressive model expressed by the linear combination, a coefficient determination step of determining the coefficient of each term in the linear combination by the least square method is further included.
[0016]
According to the invention according to the second aspect, since the coefficient of each term in the linear combination is determined by the least square method, the coefficient of each term can be accurately determined even for observation data having noise. .
[0017]
The invention according to a third aspect is the invention according to the first or second aspect, wherein the waveform data deriving step uses −1 and a frequency obtained by Fourier transforming the observation data as a variable. The waveform data is obtained by multiplying the observed data by a noise removal filter function including an exponential function having an exponential component that is a product of a coefficient derived from the half-width of the maximum peak in the second signal intensity waveform and a variable related to time. It is derived.
[0018]
According to the third aspect of the invention, the noise removal filter function includes an exponential function, the exponential component has a coefficient derived from the half width of the peak, and the noise removal filter function is multiplied by the observation data. Thus, since the waveform data is derived, the noise component of the observation data can be accurately removed, and the derived waveform data can be used as a damped oscillation function.
[0019]
The invention according to a fourth aspect is the invention according to the third aspect, wherein the waveform data deriving step is configured such that the half-value width of the maximum peak is 0.5 to 3.0 Hz. When the half-value width is the value of the coefficient, and the half-value width of the maximum peak is 0.5 Hz or less, 0.5 is the value of the coefficient, and the half-value width of the maximum peak is 3.0 Hz or more. In this case, the coefficient value is set to 3.0.
[0020]
According to the fourth aspect of the invention, since the exponent component of the noise removal filter function is specifically determined, the noise component can be removed more accurately.
[0021]
Further, the invention according to a fifth aspect is the invention according to any one of the first to fourth aspects, wherein the determining step is a step of calculating a signal strength maximum value and a noise strength value in the first signal strength waveform. The predetermined number is determined from the ratio.
[0022]
According to the fifth aspect of the invention, since the predetermined number is obtained from the ratio of the maximum value of the signal intensity and the noise intensity value in the spectrum, spectrum estimation with higher resolution than that using FPE is performed. Can be done.
[0023]
The invention according to a sixth aspect is the invention according to any one of the first to fifth aspects, wherein the determining step is a step of calculating a signal intensity maximum value and a noise intensity value in the first signal intensity waveform. An integer part of a value obtained by multiplying the value of the ratio to the (-1/2) power by a proportionality constant is defined as the predetermined number.
[0024]
According to the sixth aspect of the invention, it is possible to perform spectrum estimation with higher resolution by more specifically determining the predetermined number.
[0025]
According to a seventh aspect of the present invention, an imaging object is arranged in a magnetic field atmosphere, an electromagnetic wave having a center frequency proportional to the magnetic field intensity of the magnetic field atmosphere is transmitted to the arranged imaging object, Receiving an NMR signal obtained from the object to be imaged based on a magnetic resonance phenomenon, estimating a spectrum using an autoregressive model based on the waveform data of the received NMR signal, and using the estimated spectrum to obtain luminance information for a position In the magnetic resonance imaging method for outputting an image of the internal structure of the imaging object by replacing with, the estimation of the spectrum is a waveform data derivation step of deriving the waveform data by applying a noise removal filter to the observation data; A first signal intensity waveform having a frequency as a variable is acquired from the waveform data, and the acquired signal of the first signal intensity waveform Characterized in that the degree and the noise intensity and a determination step of determining the predetermined number.
[0026]
According to the seventh aspect of the invention, in the spectrum estimation, a noise removal filter is applied, and the predetermined number is determined from the signal intensity and the noise intensity. A spectrum is estimated, and a high-resolution magnetic resonance image can be formed based on the spectrum.
[0027]
The invention according to an eighth aspect is the invention according to the seventh aspect, wherein the waveform data deriving step uses the imaging object in a second signal intensity waveform obtained by fast Fourier transform of the observation data. Deriving the waveform data by multiplying the observed data by a denoising filter function including an exponential function having a half-width of a peak caused by hydrogen atoms constituting water molecules existing therein as a part of an exponential component. Features.
[0028]
According to the eighth aspect of the invention, the exponential component in the exponential function included in the noise removal filter function includes the half width of the peak due to the hydrogen atoms constituting the water molecules present in the imaging target. Therefore, a useful noise removal filter function can be configured, and noise components in the observation data can be effectively removed. Therefore, a magnetic resonance image having high resolution can be taken.
[0029]
Further, the invention according to a ninth aspect is the invention according to the seventh or eighth aspect, wherein the determining step is a ratio of a maximum value of the signal strength to a noise strength value in the first signal strength waveform ( The integer part of the value obtained by multiplying the value of -1/2) to the proportional constant is defined as the predetermined number.
[0030]
According to the ninth aspect of the invention, by speci? Cally defining a predetermined number of determinations, high-resolution spectrum estimation is performed from observation data having a loWer S / N ratio than when using FPE. High resolution magnetic resonance images can be taken.
[0031]
Further, the invention according to the tenth aspect is the invention according to any one of the seventh to ninth aspects, wherein the frequency difference variation between a plurality of peaks is changed from the signal intensity distribution obtained by the spectrum estimation step. The method further includes a detection step of detecting a temperature fluctuation inside the imaging object based on the detection target.
[0032]
According to the invention according to the tenth aspect, the peak due to the atoms constituting the metabolite also changes the resonance frequency due to the temperature change. Therefore, by measuring the fluctuation range, the temperature fluctuation inside the imaging target can be reduced. Can be detected.
[0033]
In addition, the invention according to the eleventh aspect is the program according to any one of the first to tenth aspects, in which any one of the methods according to the first to tenth aspects is executed on a computer. The method according to the invention can be executed by a computer.
[0034]
The invention according to the twelfth aspect includes a static magnetic field generating magnet that generates a uniform magnetic field, a gradient magnetic field generating unit that generates a magnetic field that varies according to position and time, and transmission that irradiates an imaging object with electromagnetic waves. Unit, a receiving unit that receives observation data including NMR signals from the imaging object, and an image formation processing unit for forming a magnetic resonance image, and the intensity of the uniform magnetic field and the varying magnetic field A magnetic resonance imaging apparatus for imaging a magnetic resonance image by transmitting an electromagnetic wave having a frequency proportional to the frequency and receiving an NMR signal, wherein the image formation processing unit is time-series based on the observation data and a noise removal filter function Determining the number of polynomial terms from the signal strength and noise strength in the first signal strength waveform obtained by deriving waveform data consisting of information and fast Fourier transforming the derived waveform data Characterized in that performs a spectrum estimation by autoregressive model, a magnetic resonance image pickup device for forming a magnetic resonance image by estimating the spectrum of the frequency as a variable with.
[0035]
According to the twelfth aspect of the invention, the image formation processing unit estimates a spectrum using an autoregressive model, and in the spectrum estimation, derives waveform data using a noise removal filter function, Since the magnetic resonance imaging apparatus determines the number of polynomial terms based on noise intensity, high-resolution magnetic resonance imaging is performed by performing high-resolution spectral estimation using an autoregressive model even for observation data with a low S / N ratio. Can be provided.
[0036]
The invention according to a thirteenth aspect is the invention according to the twelfth aspect, wherein the noise removal filter function is a second signal intensity waveform whose variable is a frequency obtained by fast Fourier transform of the observation data. And an exponential function having a half width of a peak caused by hydrogen atoms constituting water molecules existing inside the imaging object as part of an exponential component.
[0037]
According to the thirteenth aspect of the invention, since the noise removal filter function includes the exponential function including the half width of the peak due to the hydrogen atom in the exponential component, the image formation processing unit is effective from the observation data. In addition, it is possible to provide a magnetic resonance imaging apparatus capable of removing a noise component and capturing a high-resolution magnetic resonance image.
[0038]
The invention according to a fourteenth aspect is the invention according to the twelfth or thirteenth aspect, wherein the number of terms of the polynomial is a ratio of the signal strength and the noise strength in the first signal strength waveform ( -1/2) It is an integer part of a value obtained by multiplying a value multiplied by a proportionality constant.
[0039]
According to the fourteenth aspect of the present invention, the image formation processing unit more appropriately sets the model order for the observation data having a low S / N ratio compared to the case where the model order is determined using the FPE. A magnetic resonance imaging apparatus that can be provided can be provided.
[0040]
The invention according to a fifteenth aspect is the invention according to any one of the twelfth to fourteenth aspects, wherein the image formation processing unit has a frequency difference between a plurality of peaks in the spectrum obtained by the spectrum estimation. It is characterized in that the temperature fluctuation of the object to be imaged is detected on the basis of the fluctuation.
[0041]
According to the fifteenth aspect of the invention, since the processing means is a magnetic resonance imaging apparatus that detects temperature fluctuations inside the imaging object by detecting a difference in frequency between a plurality of peaks, direct measurement is performed. The temperature inside the object that cannot be detected can be accurately detected.
[0042]
DETAILED DESCRIPTION OF THE INVENTION
The MRI apparatus according to the first embodiment will be described below with reference to the drawings. FIG. 1 is a schematic diagram illustrating an outline of the MRI apparatus according to the first embodiment. The MRI apparatus according to the first embodiment includes a static magnetic
[0043]
The static
[0044]
The
[0045]
The RF coil unit 3 functions as an electromagnetic wave transmitting means for irradiating the imaging object 12 with an electromagnetic wave, and receives an NMR signal that receives an NMR signal emitted from atoms constituting the imaging object 12 by a nuclear magnetic resonance phenomenon. It has a function as a means. The frequency of the electromagnetic wave irradiated by the RF coil unit 3 is preferably several MHz to several tens of Mz when the strength of the magnetic field applied by the static magnetic
[0046]
The
[0047]
The
[0048]
The analog /
[0049]
The
[0050]
The
[0051]
The display unit 11 is for displaying a magnetic resonance image formed by the
[0052]
Next, a method for actually imaging a magnetic resonance image using the MRI apparatus according to the first embodiment will be described. The imaging of magnetic resonance images in the first embodiment is roughly divided into a step of acquiring observation data by applying a magnetic field, irradiating and receiving electromagnetic waves, and a step of converting the acquired observation data into a spectrum to form a magnetic resonance image It consists of.
[0053]
First, the observation data acquisition process will be described with reference to FIG. First, the imaging object 12 is carried into the MRI apparatus, and the imaging object 12 is placed in the static magnetic field atmosphere generated by the static magnetic field generation unit 1 (step S101). In addition, when the static magnetic
[0054]
The
[0055]
And the
[0056]
Then, while applying the frequency gradient magnetic field by the
[0057]
As described above, through the steps S101 to S104, it is possible to obtain observation data that is time-series information indicating at what position and how many atoms make up the imaging target 12 exist. The application of the gradient magnetic field, the irradiation of the electromagnetic wave, and the reception of the NMR signal in these steps are performed directly by the
[0058]
Note that the observation data acquisition step is performed for metabolites constituting the internal organs of the imaging target 12 and water contained in the imaging target 12. Originally, in order to know the internal structure of the object 12 to be imaged, it is only necessary to have information on the internal organs. However, in the first embodiment, information on water is also required in the subsequent spectrum estimation. is there. Specifically, observation data about metabolites and water is acquired as follows. That is, the resonance frequency differs between the hydrogen atom contained in the metabolite molecule constituting the internal organ and the hydrogen atom contained in the water molecule due to the difference in molecular structure. Therefore, the observation data regarding the internal organs and the observation data regarding water can be acquired by changing the frequency of the electromagnetic wave irradiated from the RF coil unit 3.
[0059]
Next, the process of converting the observation data acquired in the above process into a spectrum and forming a magnetic resonance image will be described with reference to FIGS. Note that the diagrams shown in FIGS. 4 to 6 are schematic diagrams for the purpose of assisting in understanding the contents of this process, and do not necessarily match actual observation data.
[0060]
First, the phase of the observation data of the metabolite is corrected with the observation data regarding water (step S201). Acquisition of observation data is performed by applying a gradient magnetic field that changes with time, so that an induced electromotive force is generated in the imaging target 12 and an eddy current is generated. This eddy current causes disturbance in the value of the observation data, and this step is performed to correct this disturbance.
[0061]
And the observation data regarding water are fast Fourier-transformed (step S202). The fast Fourier transform at this time can be performed by a conventional method, and as a result of this step, a spectrum relating to water can be obtained. Thereafter, the half width of the peak related to water is measured for the spectrum obtained in step S202 (step S203).
[0062]
And a noise removal filter function is determined from the half value width of the peak about water measured at Step S203 (Step S204). Specifically, it is assumed that the noise removal filter function is an exponential function, the exponent component is −1, the half width of the peak, and the time variable t. Here, when the half width is less than 0.5 Hz, 0.5 is substituted for the half width of the peak, and when the half width is greater than 3.0 Hz, the half of the peak Substitute 3.0 for the price range.
[0063]
Thereafter, the noise removal filter function determined in step S204 is applied to the observation data regarding the metabolite (step S205). Specifically, as shown in FIG. 4 (a), observation data l with time as a variable 1 And a noise removal filter function l with time as a variable 2 And the waveform data l shown in FIG. Three Get. This waveform data l Three FIG. 4B clearly shows that a damped oscillation waveform is obtained by multiplying the noise removal filter function. By this step, waveform data from which noise components have been removed is obtained.
[0064]
Then, fast Fourier transform is performed on the waveform data obtained in step S205 (step S206). By performing the fast Fourier transform, a spectrum of waveform data as shown in FIG. 5 can be obtained.
[0065]
Thereafter, in the spectrum of the waveform data, the maximum value P of the peak intensity shown in FIG. S And the magnitude P of the noise component N And an S / N ratio that is an intensity ratio is derived (step S207). Then, the model order M is determined based on the S / N ratio (step S208). Here, when the S / N ratio is smaller than 250, the S / N ratio is calculated as 250, and when the S / N ratio is larger than 400, the S / N ratio is calculated as 400.
[0066]
Specifically, the model order M is a value obtained by multiplying the value of the S / N ratio to the -1/2 power by a proportionality constant α. Note that the proportionality constant α is a value determined by the number of data points of measurement data and the data acquisition time interval. For example, when the number of data points is 2048, the value of the proportionality constant α is 50.
[0067]
Then, a polynomial having M number of terms is created based on the model order M determined in step S208 (step S209). Specifically, as shown in FIG. 6, a discrete time t with a mutual interval Δt. 0 , T 1 , T 2 , T Three , ..., t M , ..., t N About strength x 0 , X 1 , X 2 , X Three , ..., x M , ..., x N Is read from the waveform data, and the following polynomial is created for these.
x M + i = A 1 x M + i-1 + A 2 x M + i-2 + A Three x M + i-3 + ... + a M-1 x M + i- (M-1) + A M x M + iM (4)
Here, i = 0, 1, 2,..., N, and N = 2M. Therefore, equation (4) consists of M polynomials. X j Is a specific numerical value read from the waveform data, the unknown in equation (4) is a 1 , A 2 , ..., a M It is. Therefore, Equation (4) is M equations for M unknowns, and the coefficients in the autoregressive model are derived by solving these equations.
[0068]
Thereafter, z conversion is performed for the equation of i = 0 in the equation (4) (step S210). (4) X (z) obtained by z-converting the equation of i = 0 is expressed by the following equation.
X (z) = (1-Σa i z -1 ) -1 (5)
Where Σa i z -1 Is the sum of 1 to M for i.
[0069]
Then, variable conversion is performed on the equation (5) (step S211). Specifically, the variable z is replaced with (ifΔt / 2π). Here, i is an imaginary unit, f is a frequency, and Δt is an interval of time when waveform data is read. Specifically, as a result of the variable conversion, equation (6) is transformed into the following equation.
X (f) = (1-Σa i (IfΔt / 2π) -1 ) -1 ... (6)
Where a i Is specifically obtained from the equation (5), and Δt, 2π, and i are constants. Therefore, it can be understood that the equation (7) is a function having the frequency f (Hz) as a variable.
[0070]
Finally, the absolute value of equation (6) is obtained (step S212). This is because the equation (6) is an equation including an imaginary number i, so that it is necessary to convert the equation (7) into a function composed of natural numbers in order to know the spectrum. Specifically, the right side of equation (6) is multiplied by a conjugate complex number to obtain the square root of the product. This completes the process of converting the acquired observation data into a spectrum.
[0071]
About the spectrum obtained by this process, a graph is shown in FIG.7 (b). In FIG. 7B, spectrum estimation is performed by setting the peak half-value width, which is an exponential component of the noise removal filter function, to 1.8 Hz and the model order M to 300. Compared to FIG. 7A, which is a spectrum obtained by fast Fourier transform based on the same observation data, the half width of the spectrum peak is clearly narrow. Further, as a secondary effect, it can be seen that the magnitude of the noise component in the resulting spectrum is clearly smaller in the graph of FIG. 7B by multiplying the noise removal filter function by the observation data.
[0072]
Thereafter, a two-dimensional magnetic resonance image for each slice plane is created on a computer based on the spectrum obtained in the above process and output to the display unit 11, whereby the MRI apparatus according to the first embodiment is changed. Imaging of the used magnetic resonance image is completed.
[0073]
The reason why the noise removal filter function determined in step S204 is an exponential function is as follows. Since the irradiation of the electromagnetic wave by the RF coil unit 3 in the image data acquisition process is performed in a short time, the intensity of the NMR signal emitted from the atoms constituting the imaging target 12 corresponding to this is the electromagnetic wave by the RF coil unit 3 After the end of irradiation, it should decrease exponentially and form a damped oscillation waveform. However, in the observation data, the straight line l in FIG. 1 As shown, it does not necessarily decrease exponentially. This is because the decrease in the NMR signal emitted based on the nuclear magnetic resonance phenomenon is canceled because the intensity of the noise component is large. Therefore, by multiplying the observation data by the exponential decrease function, the noise component is removed from the observation data, and only the NMR signal component based on the nuclear magnetic resonance phenomenon can be extracted as the waveform data.
[0074]
The reason for determining the exponential component of the noise removal filter function in step S204 is proportional to the half-value width of the peak related to water for the following reason. That is, the half-value width of each spectrum corresponds to the damped oscillation component of the observation data that is time-series data, and the noise component can be effectively removed by using the half-value width of the peak in the spectrum. Moreover, it is desirable to use a peak with a high peak intensity as a peak for measuring the half width, and a peak due to a hydrogen atom constituting a water molecule is selected.
[0075]
As described above, the MRI apparatus according to the first embodiment can accurately perform spectrum estimation as compared with imaging of a magnetic resonance image using conventional fast Fourier transform. As a result, it is possible to capture a magnetic resonance image having a high resolution.
[0076]
In addition, spectrum estimation in the first embodiment can perform spectrum estimation with high resolution using the autoregressive model even for observation data with a small S / N ratio. Therefore, in addition to magnetic resonance imaging using an MRI apparatus, for example, obtaining a spectrum with a high resolution from time-series information with a low S / N ratio in the analysis of weather phenomena and the characteristics of electronic circuits and mechanical systems. Can do.
[0077]
In the observation data acquisition process, steps S102, S103, and S104 are steps for specifying the position of the atom in the imaging target 12, as described above. On the other hand, the observation data related to water is only for use in spectral conversion of the observation data related to metabolites. It may be omitted. For the same reason, it is not necessary to obtain a two-dimensional data set when acquiring observation data relating to water, and repeating steps S102 to S104 may be omitted.
[0078]
In addition, although the waveform data is reduced in noise by the noise removal filter function, it is not completely removed. Therefore, strictly speaking, in step S209, the equal sign between the left side and the right side may not be established in the expression (4). Therefore, more precisely, a 1 , A 2 , ..., a M It is desirable to derive Specifically, first, the equation (4) regarding i = 0, 1, 2,..., NM is expressed by a determinant. Where (x M , X M + 1 , ..., x N ) Is defined as d, and (a 1 , A 2 , A Three ... a M ) Is a column vector. M kl = X M-1 + kl Let H be an M × N matrix whose elements are. With these d, a, and H, the equation (4) is
d = H · a (7)
It can be expressed in the form of Therefore, by the least square method, a is
a = H t H -1 H t d (8)
A can be specifically obtained by taking the equation (8). By using the least squares method, a 1 , A 2 ... a M Can be requested. Other than this, the Yule-Walker method and Burg method can be used to 1 , A 2 , A Three ... a M It is also effective to make a decision.
[0079]
In FIG. 1, the RF coil unit 3 is configured to cover only the patient's head, which is the imaging object 12, but the MRI apparatus according to the first embodiment is not limited to this structure. For example, the shape of the RF coil unit 3 may be changed in order to capture an image of the patient's abdomen, or the shape may be such that the entire body of the patient can be imaged. Further, the MRI apparatus shown in FIG. 1 is a so-called horizontal magnetic field type MRI apparatus that applies a rectifying magnetic field in the horizontal direction, but other than this, for example, a so-called vertical magnetic field type that applies a static magnetic field in the vertical direction. The MRI apparatus which consists of these structures may be sufficient.
[0080]
Furthermore, in the first embodiment, as shown in FIG. 1, regarding the MRI structure, the parts that control the apparatus, such as the
[0081]
In the observation data acquisition process according to the first embodiment, the slice direction is normally the direction from the patient's head, which is the imaging object 12, toward the toes, but is not limited thereto. The slice direction can be changed according to the region of the magnetic resonance image to be imaged. In addition, the frequency direction, the phase direction, and the slice direction usually form an orthogonal coordinate system, but each direction in the first embodiment is not necessarily limited to the orthogonal coordinate system.
[0082]
In addition, regarding the spectrum estimation in the first embodiment, the operator may manually perform spectrum estimation through the
[0083]
Next, an MRI apparatus according to the second embodiment will be described. The MRI apparatus according to the second embodiment has a function of detecting a temperature change inside the imaging target as well as a function of imaging a magnetic resonance image. The structure of the MRI apparatus according to the second embodiment is the same as the structure of the MRI apparatus according to the first embodiment shown in FIG.
[0084]
In addition, the observation data acquisition process and spectrum estimation in the second embodiment are performed in substantially the same manner as in the first embodiment. And the temperature fluctuation in the inside of the imaging target 12 is observed from the spectrum obtained by these processes.
[0085]
Specifically, after spectrum estimation is performed in the same manner as in the first embodiment, the peak frequency attributed to hydrogen atoms constituting water and the hydrogen atoms constituting metabolites are generated on the
[0086]
Then, by comparing the read frequency difference with reference data, it is detected how the internal temperature of the imaging object 12 has changed. For example, data relating to a frequency difference at 36 degrees is acquired in advance, and the frequency difference read in the second embodiment is compared with known data to examine how much the internal temperature has changed from 36 degrees. .
[0087]
In general, it is known that the resonance frequency of a metabolite varies depending on temperature. On the other hand, since the resonance frequency of the hydrogen atom constituting the water molecule does not depend on the temperature, the difference between the resonance frequency of the hydrogen atom constituting the metabolite and the resonance frequency of the hydrogen atom constituting the water molecule is obtained by calculating the temperature. It is possible to detect fluctuations. Since the fluctuation of the resonance frequency of the hydrogen atom constituting the metabolite was a slight value, it was difficult to detect by spectrum estimation by conventional fast Fourier transform. However, in the spectrum estimation according to the first embodiment, the spectrum estimation can be performed with high accuracy, so that the fluctuation of the resonance frequency with respect to the temperature change can be sufficiently detected. That is, the MRI apparatus according to the second embodiment can detect temperature fluctuations in the inside of an object that cannot be directly measured for temperature.
[0088]
Thus, the following advantage arises by measuring the temperature fluctuation inside the imaging object 12. For example, the MRI apparatus according to the second embodiment can be used as a diagnostic apparatus related to the human brain.
[0089]
In brain diseases such as brain contusions, it is extremely important to know the temperature of the brain. That is, when the temperature of the brain rises, the brain volume increases due to brain edema, and the patient may die. In such a case, if an increase in brain temperature can be detected by the MRI apparatus according to the second embodiment, the patient can be treated by performing appropriate treatment, for example, hypothermia therapy.
[0090]
In the second embodiment, since the imaging of the magnetic resonance image itself is not intended, in the observation data acquisition step, the step of applying the gradient magnetic field and imaging is omitted, and only the static magnetic field is applied, The RF coil unit 3 may transmit electromagnetic waves. This is because it is sufficient to know the resonance frequency of the metabolite and the resonance frequency of water in order to detect the temperature fluctuation inside the imaging object 12. On the contrary, the observation data acquisition process may be performed in the same manner as in the first embodiment, and the magnetic resonance image of the imaging object 12 may be captured together with the measurement of the temperature fluctuation.
[0091]
Although the spectrum estimation method using the autoregressive model has been described in the first and second embodiments, the spectrum estimation method using the autoregressive model is the same as the so-called MEM (Maximum Entropy Method) spectrum estimation method. (SL Marple, Digital Spectral Analysis, Prentice-Hall, 1987). Therefore, the present invention is naturally applicable to this MEM.
[0092]
【The invention's effect】
As described above, according to the invention according to the first aspect, the noise reduction filter is applied for the derivation of the waveform data, and the predetermined number is determined from the signal intensity and the noise intensity. Even for observation data with low resolution, it is possible to perform spectrum estimation with high resolution.
[0093]
Further, according to the invention according to the second aspect, since the coefficient of each term in the linear combination is determined by the least square method, the coefficient of each term can be accurately determined even for observation data having noise. There is an effect that can be done.
[0094]
According to the third aspect of the invention, the noise removal filter function includes an exponential function, the exponential component has a coefficient derived from the half width of the peak, and the noise removal filter function is used as observation data. Since the waveform data is derived by multiplication, the noise component of the observation data can be accurately removed, and the derived waveform data can be used as a damped oscillation function.
[0095]
In addition, according to the fourth aspect of the invention, since the exponent component of the noise removal filter function is specifically determined, there is an effect that the noise component can be removed more accurately.
[0096]
Further, according to the fifth aspect of the invention, since the predetermined number is obtained from the ratio of the maximum value of the signal intensity and the noise intensity value in the spectrum, spectrum estimation with higher resolution than when FPE is used. There is an effect that can be performed.
[0097]
In addition, according to the sixth aspect of the invention, the configuration in which the predetermined number of terms is determined more specifically provides an effect that spectrum estimation with higher resolution can be performed.
[0098]
Further, according to the seventh aspect of the invention, in the spectrum estimation, a noise removal filter is applied and a predetermined number is determined from the signal intensity and the noise intensity, so that the high resolution can be obtained from the observation data having a low S / N ratio. The spectrum is estimated and a high-resolution magnetic resonance image can be formed based on the spectrum.
[0099]
Further, according to the eighth aspect of the invention, the exponential component in the exponential function included in the noise removal filter function includes the half width of the peak due to the hydrogen atoms constituting the water molecules present inside the imaging object. Therefore, a useful noise removal filter function can be configured, and noise components in the observation data can be effectively removed. Therefore, it is possible to capture a magnetic resonance image having a high resolution.
[0100]
In addition, according to the ninth aspect of the invention, by specifying the predetermined number of decisions specifically, high-resolution spectrum estimation can be performed from observation data having a lower S / N ratio than when FPE is used. This produces an effect that a high-resolution magnetic resonance image can be taken.
[0101]
Also, according to the tenth aspect of the invention, since the peak due to the atoms constituting the metabolite also changes the resonance frequency due to temperature change, the temperature fluctuation inside the imaging object is measured by measuring the fluctuation range. There is an effect that can be detected.
[0102]
According to the eleventh aspect of the invention, any one of the first to tenth aspects of the invention is a program for executing the method on a computer. The method according to the aspect of the invention can be executed by a computer.
[0103]
According to the twelfth aspect of the invention, the image formation processing unit estimates a spectrum using an autoregressive model, and in the spectrum estimation, waveform data is derived using a noise removal filter function, and the signal intensity Since the magnetic resonance imaging device determines the number of polynomial terms based on the noise intensity, high-resolution magnetic resonance is estimated using an autoregressive model even for observation data with a low S / N ratio. There is an effect that a device capable of capturing an image can be provided.
[0104]
According to the thirteenth aspect of the invention, since the noise removal filter function includes the exponential function including the half width of the peak due to the hydrogen atom in the exponential component, the image formation processing unit can obtain the effect from the observation data. In particular, the noise component can be removed and a magnetic resonance imaging apparatus capable of capturing a high-resolution magnetic resonance image can be provided.
[0105]
According to the fourteenth aspect of the present invention, the model order is more appropriately applied to observation data having a low S / N ratio compared to the case where the image forming processing unit determines the model order using FPE. It is possible to provide a magnetic resonance imaging apparatus capable of providing
[0106]
According to the fifteenth aspect of the invention, since the processing means is a magnetic resonance imaging apparatus that detects temperature fluctuations inside the imaging object by detecting a frequency difference between a plurality of peaks, There is an effect that the temperature inside the object that cannot be measured can be accurately detected.
[Brief description of the drawings]
FIG. 1 is a schematic diagram showing a structure of an MRI apparatus according to a first embodiment;
FIG. 2 is a flowchart showing a process of obtaining observation data in the first embodiment.
FIG. 3 is a flowchart showing a spectrum estimation process in the first embodiment.
4A is a schematic diagram of multiplying observation data with a noise removal filter function in spectrum estimation in
5 is a schematic diagram showing a peak intensity and a noise intensity when determining a model order M in spectrum estimation according to
6 is a schematic diagram when deriving a polynomial by an autoregressive model from observed data in spectrum estimation in
FIG. 7A is a spectrum diagram obtained by performing fast Fourier transform on observation data obtained by the MRI apparatus according to the first embodiment; (B) is the spectrum figure which derived | led-out the same observation data using the spectrum estimation method in this
8A is a graph showing observation data that is time-series information, and FIG. 8B is a spectrum diagram obtained by performing a fast Fourier transform on the data in FIG. 8A. Further, (c) is a spectrum diagram obtained by performing spectrum estimation from the data of (a) using an autoregressive model.
[Explanation of symbols]
1 Static magnetic field generator
2 Gradient coil
3 RF coil section
4 Gradient drive
5 Transmitter
6 Receiver
7 Analog / digital converter
8 Computer
9 Control unit
10 Operation part
11 Display
Claims (9)
観測データにノイズ除去フィルタを適用して前記波形データを導出する波形データ導出工程と、
前記波形データから周波数を変数とする第１の信号強度波形を取得し、該取得した第１の信号強度波形における信号強度の最大値と雑音強度値との比から前記所定数を決定する決定工程と、を有することを特徴とするスペクトル推定方法。A spectrum estimation method for estimating a signal strength waveform having a frequency as a variable based on signal strength at a predetermined time of waveform data having time as a variable and a predetermined number of signal strengths at discrete times up to the predetermined time. And
A waveform data deriving step of deriving the waveform data by applying a noise removal filter to the observation data;
A determination step of acquiring a first signal intensity waveform having a frequency as a variable from the waveform data and determining the predetermined number from a ratio between a maximum value of signal intensity and a noise intensity value in the acquired first signal intensity waveform. And a spectrum estimation method characterized by comprising:
前記画像形成処理部は、前記観測データとノイズ除去フィルタ関数とに基づいて時系列情報からなる波形データを導出し、導出された波形データを高速フーリエ変換して得られた第１の信号強度波形における信号強度と雑音強度から多項式の項数を決定しつつ自己回帰モデルによるスペクトル推定を行い、周波数を変数とするスペクトルを推定することにより磁気共鳴映像を形成し、
前記多項式の項数は、前記第１の信号強度波形における前記信号強度の最大値と雑音強度値との比から決定することを特徴とする磁気共鳴映像撮像装置。A static magnetic field generating magnet that generates a uniform magnetic field, a gradient magnetic field generating unit that generates a magnetic field that varies according to position and time, a transmission unit that irradiates an electromagnetic wave to the imaging target, and an NMR signal from the imaging target A reception unit for receiving observation data including the image data, and an image formation processing unit for forming a magnetic resonance image, and transmitting an electromagnetic wave having a frequency proportional to the intensity of the uniform magnetic field and the varying magnetic field and an NMR signal A magnetic resonance imaging apparatus for imaging magnetic resonance imaging by receiving
The image formation processing unit derives waveform data composed of time-series information based on the observation data and a noise removal filter function, and a first signal intensity waveform obtained by fast Fourier transform of the derived waveform data The spectrum is estimated by the autoregressive model while determining the number of terms of the polynomial from the signal intensity and noise intensity at, and the magnetic resonance image is formed by estimating the spectrum with the frequency as a variable.
The number of terms of the polynomial is determined from a ratio between a maximum value of the signal intensity and a noise intensity value in the first signal intensity waveform.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2001311917A JP3842607B2 (en) | 2001-10-09 | 2001-10-09 | Spectral estimation method and magnetic resonance imaging apparatus |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2001311917A JP3842607B2 (en) | 2001-10-09 | 2001-10-09 | Spectral estimation method and magnetic resonance imaging apparatus |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2003116811A JP2003116811A (en) | 2003-04-22 |
| JP3842607B2 true JP3842607B2 (en) | 2006-11-08 |
Family
ID=19130658
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2001311917A Expired - Fee Related JP3842607B2 (en) | 2001-10-09 | 2001-10-09 | Spectral estimation method and magnetic resonance imaging apparatus |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3842607B2 (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2012135466A (en) * | 2010-12-27 | 2012-07-19 | Ge Medical Systems Global Technology Co Llc | Magnetic resonance imaging apparatus |
| JP6160366B2 (en) * | 2013-08-26 | 2017-07-12 | 日本電子株式会社 | NMR signal processing system |
| CN111398755B (en) * | 2020-04-21 | 2022-07-12 | 武汉朕泰智能科技有限公司 | Cable partial discharge waveform extraction method based on short-time FFT (fast Fourier transform) segmentation technology |
-
2001
- 2001-10-09 JP JP2001311917A patent/JP3842607B2/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JP2003116811A (en) | 2003-04-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP5198859B2 (en) | MRI temperature measurement involving phase mapping and reference medium used as phase reference | |
| US9720066B2 (en) | Magnetic resonance imaging apparatus and control method thereof | |
| US20140125333A1 (en) | Magnetic resonance imaging apparatus and gradient magnetic field waverform estimation method | |
| CN102859387B (en) | Adopt the magnetic resonance imaging of the chemical species of spectrum model | |
| WO2001017428A1 (en) | Magnetic resonance imaging device and method therefor | |
| KR20070042875A (en) | Magnetic resonance imaging apparatus | |
| US8049497B2 (en) | MRI RF encoding using multiple transmit coils | |
| JP2005270285A (en) | Magnetic resonance imaging apparatus and magnetic resonance imaging data processing method | |
| JP5259715B2 (en) | Magnetic resonance imaging apparatus and execution order determination method | |
| JP6038654B2 (en) | Magnetic resonance imaging apparatus and vibration error magnetic field reduction method | |
| US10698060B2 (en) | Magnetic resonance imaging apparatus and medical image processing apparatus | |
| JP2019195588A (en) | Magnetic resonance imaging apparatus and multi-slice imaging method | |
| JP4781120B2 (en) | Magnetic resonance imaging apparatus and magnetic resonance spectrum measuring method | |
| JPH064066B2 (en) | Magnetic resonance imaging equipment | |
| JP3842607B2 (en) | Spectral estimation method and magnetic resonance imaging apparatus | |
| JP5465565B2 (en) | Magnetic resonance imaging system | |
| CN110074786B (en) | Nuclear magnetic resonance shimming method and device, computing equipment and nuclear magnetic resonance imaging system | |
| JP4912802B2 (en) | Magnetic resonance imaging apparatus, transmission sensitivity distribution measuring apparatus, and transmission sensitivity distribution measuring method | |
| JP2014045775A (en) | Magnetic resonance imaging apparatus | |
| JP5167556B2 (en) | Magnetic resonance imaging system | |
| JPH07265277A (en) | Mr bloodstream information collecting method and mri device | |
| JP6157976B2 (en) | Magnetic resonance imaging apparatus and method | |
| JP3884282B2 (en) | MRI equipment | |
| JP6727154B2 (en) | Magnetic resonance imaging equipment | |
| JP2015019813A (en) | Magnetic resonance imaging apparatus and method for calculating correcting b1 map |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A625 | Written request for application examination (by other person) |
Free format text: JAPANESE INTERMEDIATE CODE: A625 Effective date: 20040406 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20060124 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060131 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060407 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20060711 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20060810 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100818 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100818 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100818 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110818 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110818 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120818 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120818 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120818 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120818 Year of fee payment: 6 |
|
| LAPS | Cancellation because of no payment of annual fees |