7649263 磁気共鳴イメージング装置、画像解析装置及び流体の解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2025-03-11

(45)【発行日】2025-03-19

(54)【発明の名称】磁気共鳴イメージング装置、画像解析装置及び流体の解析方法

(51)【国際特許分類】

   A61B 5/055 20060101AFI20250312BHJP

【ＦＩ】

A61B5/055 382

A61B5/055 380

【請求項の数】 12

(21)【出願番号】P 2022003914

(22)【出願日】2022-01-13

(65)【公開番号】P2023103067

(43)【公開日】2023-07-26

【審査請求日】2024-04-01

(73)【特許権者】

【識別番号】306037311

【氏名又は名称】富士フイルム株式会社

(74)【代理人】

【識別番号】110000888

【氏名又は名称】弁理士法人山王坂特許事務所

(72)【発明者】

【氏名】尾藤 良孝

(72)【発明者】

【氏名】越智 久晃

【審査官】佐野 浩樹

(56)【参考文献】

【文献】特表２０１９－５３４０８６（ＪＰ，Ａ）

【文献】特開２００５－０１１０４３（ＪＰ，Ａ）

【文献】特許第６４６７３４１（ＪＰ，Ｂ２）

【文献】特開２０１６－１３１８３６（ＪＰ，Ａ）

【文献】特開２００８－０５４７２１（ＪＰ，Ａ）

【文献】特開平０５－０３１０９９（ＪＰ，Ａ）

【文献】特表２００５－５２５２０７（ＪＰ，Ａ）

【文献】Yoshitaka Bito et al.，Low b-value diffusion tensor imaging for measuring pseudorandom flow of cerebrospinal fluid，Magnetic Resonance in Medicine，2021年，Vol. 86，p. 1369-1382，<検索日：2024.10.31>, < DOI: 10.1002/mrm.28806>, <URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.28806>

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ａ６１Ｂ ５／０５５

(57)【特許請求の範囲】

【請求項1】

ＭＰＧパルスを含む拡散強調イメージングのパルスシーケンスを実行し計測データを取得する計測部と、
前記計測部が取得した計測データを用いて、検査対象内の流体の疑似拡散テンソルを算出する演算部と、を備え、
前記演算部は、前記疑似拡散テンソルと、流速のボクセル内分布の推定モデルとを用いて、前記疑似拡散テンソル以外の流体パラメータとして、壁せん断応力（ＷＳＳ）、運動エネルギー（ＫＥ）或いはエネルギーロス（ＥＬ）、振動せん断インデックス（ＯＳＩ）、圧力勾配（ＰＧ）のいずれか１以上を含む流体パラメータを算出する流体パラメータ算出部を含み、
前記流体パラメータ算出部は、前記疑似拡散テンソルを用いて流速分布の分散を算出するとともに、前記推定モデルとして、ボクセル内の流速分布を矩形分布又はガウス分布と仮定した推定モデルを用い、算出した流速分布の分散と前記推定モデルとを用いて、流速微分ベクトルを算出し、当該流速微分ベクトルに基づき前記流体パラメータを算出することを特徴とする磁気共鳴イメージング装置。

【請求項2】

請求項１に記載の磁気共鳴イメージング装置であって、
前記演算部は、前記流体パラメータの値を画素値とする画像を生成する画像生成部をさらに備え、前記画像を表示装置に表示させることを特徴とする磁気共鳴イメージング装置。

【請求項3】

請求項１に記載の磁気共鳴イメージング装置であって、
前記推定モデルは流体における線形流からの乖離を示す指標を含み、
前記演算部は、前記疑似拡散テンソルから当該指標を算出することを特徴とする磁気共鳴イメージング装置。

【請求項4】

請求項３に記載の磁気共鳴イメージング装置であって、
前記演算部は、複数の拡散時間で算出された前記疑似拡散テンソルから前記指標を算出することを特徴とする磁気共鳴イメージング装置。

【請求項5】

請求項１に記載の磁気共鳴イメージング装置であって、
前記計測部は、前記拡散強調イメージングを検査対象の周期的体動と同期して計測データを収集し、
前記流体パラメータ算出部が、前記周期的体動の情報を加えて前記流体パラメータを算出することを特徴とする磁気共鳴イメージング装置。

【請求項6】

請求項１に記載の磁気共鳴イメージング装置であって、
前記計測部は、前記拡散強調イメージングのパルスシーケンスを実行して複素数値の計測データを取得し、
前記演算部は、前記複素数値の計測データの位相情報を用いて、検査対象内の流体の流速を算出し、
前記流体パラメータ算出部は、前記疑似拡散テンソルと、前記流速と、前記推定モデルとを用いて、前記流体パラメータを算出することを特徴とする磁気共鳴イメージング装置。

【請求項7】

請求項１に記載の磁気共鳴イメージング装置であって、
前記計測部は、前記拡散強調イメージングのパルスシーケンスを用いた第１の計測と、フェイズコントラスト法のパルスシーケンスを用いた第２の計測とを行い、
前記演算部は、前記第２の計測で得た計測データを用いて前記流体パラメータを算出する第２流体パラメータ算出部をさらに備えることを特徴とする磁気共鳴イメージング装置。

【請求項8】

請求項７に記載の磁気共鳴イメージング装置であって、
前記演算部は、前記流体パラメータ算出部が算出した流体パラメータと前記第２流体パラメータ算出部が算出した流体パラメータとを統合するパラメータ統合部をさらに備えることを特徴とする磁気共鳴イメージング装置。

【請求項9】

請求項７に記載の磁気共鳴イメージング装置であって、
計測対象である流体に関する条件を受け付ける入力部をさらに備え、
前記計測部は、前記入力部が受け付けた条件に応じて、前記拡散強調イメージングのパルスシーケンス及び前記フェイズコントラスト法のパルスシーケンスを択一的に行うことを特徴とする磁気共鳴イメージング装置。

【請求項10】

磁気共鳴イメージング装置において計測された計測データまたは当該計測データから計算された疑似拡散テンソルの値を受け付ける受付部と、
前記受付部が受け付けた計測データまたは疑似拡散テンソルの値を用いて、流体パラメータを算出する演算部と、を備え
前記演算部は、前記疑似拡散テンソルと、流速のボクセル内分布の推定モデルとを用いて、前記疑似拡散テンソル以外の流体パラメータとして、壁せん断応力（ＷＳＳ）、運動エネルギー（ＫＥ）或いはエネルギーロス（ＥＬ）、振動せん断インデックス（ＯＳＩ）、圧力勾配（ＰＧ）のいずれか１以上を含む流体パラメータを算出する流体パラメータ算出部を含み、
前記流体パラメータ算出部は、前記疑似拡散テンソルを用いて流速分布の分散を算出するとともに、前記推定モデルとして、ボクセル内の流速分布を矩形分布又はガウス分布と仮定した推定モデルを用い、算出した流速分布の分散と前記推定モデルとを用いて、流速微分ベクトルを算出し、当該流速微分ベクトルに基づき前記流体パラメータを算出することを特徴とする画像解析装置。

【請求項11】

請求項１０に記載の画像解析装置であって、
前記流体パラメータ算出部は、前記疑似拡散テンソルから、流速分布の分散を算出する分散算出部と、算出した流速分布の分散と前記推定モデルとを用いて、流速微分ベクトルを算出する微分算出部と、をさらに備えることを特徴とする画像解析装置。

【請求項12】

流体を含む検査対象について磁気共鳴イメージング装置が取得した拡散テンソル画像を用いた流体の解析方法であって、
前記拡散テンソル画像を用いて流体の流速分布の分散を算出するステップと、
前記分散とボクセル内の流速の分布形状とを用いて流速の微分値を算出するステップと、
前記微分値を用いて、流体の流れの特徴を示す流体パラメータを算出するステップと、
を含み、
前記流速の微分値を算出するステップは、前記ボクセル内の流速の分布形状として矩形分布又はガウス分布を用い、
前記流体パラメータを算出するステップは、壁せん断応力（ＷＳＳ）、運動エネルギー（ＫＥ）或いはエネルギーロス（ＥＬ）、振動せん断インデックス（ＯＳＩ）、圧力勾配（ＰＧ）のいずれか１以上を含む流体パラメータを算出することを特徴とする流体の解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、磁気共鳴イメージング（以下、ＭＲＩと略す）装置及びＭＲＩ装置で取得した計測データの解析技術に関し、特に拡散テンソル画像の解析技術に関する。

【背景技術】

【0002】

人の脳や体内を流れる流体の流れを解析することにより、種々の病変の診断に有効な指標が得られる。例えば、壁せん断応力（ＷＳＳ：wall shear stress）や流体の運動エネルギー（ＫＥ: kinetic energy）或いは運動エネルギーロス（ＥＬ）など流体パラメータは、血液や脳脊髄液の流れの異常を検出するために有効な指標であり、それを計測・解析する技術が種々提案されている。

【0003】

その一つとして、特許文献１では、ＭＲＩのフェイズコントラスト（ＰＣ）法を用いて、得られた計測データから運動エネルギーロス等を解析する手法を開示している。ＰＣ法では、ＮＭＲ信号の位相をディフェイズする傾斜磁場パルスとリフェイズする傾斜磁場パルスとを組み合わせることで、静止部分からの信号と流れのある部分からの信号との位相を異ならせて、所望の領域内を流れる血液などの流体の速度を算出する。特許文献１に記載された技術では、ＰＣ法で得られた流体の速度ベクトル値のマップ（速度ベクトル値を画素値とする血流ベクトル画像）から、隣接ベクトルの差分情報を算出し、この差分を用いてＥＬ等を算出している。

【0004】

また、流体を観察できるＭＲＩの撮像方法として、拡散強調イメージング（ＤＷＩ）や拡散テンソルイメージング（ＤＴＩ）がある。ＤＴＩは、少なくとも６方向のＭＰＧパルスを用いたＤＷＩを行って得た計測データから、ボクセル毎に拡散の異方性を表す情報を拡散テンソルで示す手法であり、脳脊髄液の疑似ランダム流などの観察に有効であることが、非特許文献１に報告されている。

【先行技術文献】

【特許文献】

【0005】

【文献】特許第６４６７３４１号明細書

【非特許文献】

【0006】

【文献】Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ ｉｎ Ｍｅｄｉｃｉｎｅ ２０２１；８６、１３６９―１３８２

【発明の開示】

【発明が解決しようとする課題】

【0007】

ＰＣ法では、位相を拡散するディフェイズ傾斜磁場パルスに対し、所定の時間間隔を置いてリフェイズ傾斜磁場パルスを印加し、信号を取得するパルスシーケンスを繰り返し実施する必要があり、画像の空間分解能が比較的低い。空間分解能が低い血流ベクトル画像において隣接ベクトルの差分を取った場合、流速の微分値は、実際の微分値よりも低く計算されることになり、それをもとに算出されるＥＬ等の血流パラメータ値の誤差も増大するおそれがある。特に流体の速度が遅い場合には、誤差は拡大しやすい。

【0008】

本発明は、ＭＲＩを用いて流体パラメータを精度よく算出する技術を提供することを課題とする。特に速度が遅い流体について、流体パラメータ算出精度を向上することを課題とする。

【課題を解決するための手段】

【0009】

上記課題を解決するため、本発明は拡散テンソルイメージングで得られる流速情報と、流速分布のモデルとを用いて、流速の微分ベクトルを算出し、流体パラメータを算出する。

【0010】

即ち本発明のＭＲＩ装置は、ＭＰＧパルスを含む拡散強調イメージングのパルスシーケンスを実行し計測データを取得する計測部と、計測部が取得した計測データを用いて、検査対象内の流体の流速と拡散の情報を有する疑似拡散テンソルを算出する演算部と、を備え、演算部は、疑似拡散テンソルと、流速のボクセル内分布の推定モデルとを用いて、疑似拡散テンソル以外の流体パラメータを算出する流体パラメータ算出部を含む。

【0011】

また本発明の画像解析装置は、上記ＭＲＩ装置の演算部の機能を持つ独立した解析装置であり、ＭＲＩ装置において計測された計測データまたは当該計測データから計算された疑似拡散テンソルの値を受け付ける受付部と、受付部が受け付けた計測データまたは疑似拡散テンソルの値を用いて、流体パラメータを算出する演算部と、を備え、演算部は、疑似拡散テンソルと、流速の分布の推定モデルとを用いて、流体の流れの特徴を示す指標となる流体パラメータを算出する流体パラメータ算出部を含む。

【0012】

さらに本発明は、流体を含む検査対象について磁気共鳴イメージング装置が取得した拡散テンソル画像を用いた流体の解析方法であって、拡散テンソル画像を用いて流体の流速の分散を算出するステップと、分散とボクセル内の流速の分布形状とを用いて流速の微分ベクトルを算出するステップと、微分ベクトルを用いて、流体の流れの特徴を示す流体パラメータを算出するステップと、を含む流体の解析方法を提供する。

【発明の効果】

【0013】

本発明は、ボクセル内の流体について推定分布を用いることで、ＤＴＩから流体パラメータを算出するという新規な手法を提供する。本発明は、ボクセル内の平均化した流速とそのボクセル間での差分を用いるＰＣ法の手法に比べ、ボクセル内の流速分布を利用しており、且つ低分解能に伴う算出誤差を回避できるので、精度良く流体パラメータを算出することができる。特に速度の遅い流体について、空間分解の低下に伴う算出精度の低下を抑制することができる。

【図面の簡単な説明】

【0014】

【図1】ＭＲＩ装置の全体概要を示す図。

【図2】演算部の機能ブロック図。

【図3】流体パラメータ算出部の処理の流れを示すフロー。

【図4】拡散テンソルイメージングに用いるパルスシーケンスの一例を示す図。

【図5】各種流れと、それに対応する拡散テンソル（ＤＴ）との関係を示す図。

【図6】脳画像及びその拡散テンソル画像の一例を示す図。

【図7】線形流と拡散傾斜磁場との関係を模式的に示す図。

【図8】流速分布の推定モデルの例を示す図。

【図9】線形流と拡散とが混じった流れを模式的に示す図。

【図10】実施形態２の演算部の一例を示す機能ブロック図。

【図11】実施形態３の演算部の機能ブロック図。

【図12】（A）、（B）はそれぞれ実施形態３のＧＵＩの例を示す図。

【発明を実施するための形態】

【0015】

以下、本発明のＭＲＩ装置の実施形態を、図面を参照して説明する。

【0016】

まず本発明が適用されるＭＲＩ装置の構成を説明する。ＭＲＩ装置１の構成は、図１に示すように、一般的なＭＲＩ装置の構成と同様であり、被検体が置かれる検査空間に均一な静磁場を発生するマグネット１１と、マグネット１１が発生する静磁場に対し磁場勾配を与える傾斜磁場コイル１２と、被検体にパルス状の高周波磁場を印加し被検体の組織を構成する原子の原子核に核磁気共鳴を起こさせる送信コイル及び被検体から発生する核磁気共鳴信号を受信する受信コイルを備えたプローブ１３と、受信コイルに接続された受信器１４と、送信コイルが接続された高周波磁場発生器１５と、傾斜磁場コイル１２が接続された傾斜磁場電源１６と、受信器１４、高周波磁場発生器１５及び傾斜磁場電源１６を所定のパルスシーケンスに従って制御するシーケンサ１７と、計算機２００とを備えている。上述した要素のうち計算機２００を除く要素を包括して、ここでは計測部１０という。

【0017】

計測部１０の受信器１４が受信した核磁気共鳴信号は、デジタル化され計測データとして計算機２００に渡される。

【0018】

計測部１０を構成する各部の構造、機能などは、公知のＭＲＩ装置と同様であり、また本発明は公知の種々のタイプのＭＲＩ装置や要素に適用することができるので、ここでは計測部１０の詳細な説明は省略する。

【0019】

計算機２００は、ＭＲＩ装置全体を制御する制御系及び計測部１０の計測データに対し画像生成を含む種々の演算を行う演算部２０として機能するもので、ＣＰＵやメモリを備えた汎用の計算機やワークステーションで構成することができる。また計算機２００は、ＭＲＩ装置とデータの授受が可能な、ＭＲＩ装置とは独立した計算機やワークステーションでもよく、その場合、演算部２０の機能の全部または一部を独立した装置（画像解析装置）で行うことができる。このような画像解析装置は、ＭＲＩ装置からのデータを受け付ける受付部を有し、以下説明する計算機２００と同様の画像解析機能を実現する。

【0020】

また計測データの処理について、ＰＬＣ（Programable logic controller）やＡＳＩＣ、ＦＰＧＡ等のハードウェアを用いることもでき、それらも広い意味で演算部２０に包含される。

【0021】

本実施形態のＭＲＩ装置は、計測部１０が、拡散強調イメージング（ＤＷＩ）のパルスシーケンスを用いた計測を行い、演算部２０がＤＷＩで得られた計測データを用いて、被検体の体内を流れる血液や脳脊髄液などの流体の流れの異常を検出するための流体パラメータを算出することが特徴であり、そのため、計算機２００は、シーケンサ１７（計測制御部）に指令を送り、計測部１０が入力装置（入力部）３０等を介してユーザが指定した撮像条件（撮像パラメータ）でＤＷＩシーケンスを実行するように制御するとともに、演算部２０は計測データを用いて上述した流体パラメータを算出するための演算を行う。また必要に応じて、算出したパラメータの値を画素値とするパラメータ画像を生成し、表示装置４０に表示させたり、記憶媒体５０に格納或いは医療情報データベース等の別の記憶装置等に転送する。

【0022】

上記機能を実現する演算部２０の機能ブロック図の一例を図２に示す。図示するように、演算部２０は、ＤＷＩの計測データを処理し、疑似拡散テンソル（流速と拡散の混合した情報を表すテンソル）等を算出するＤＴ算出部２１、流体ベクトルを用いて流体パラメータを算出する流体パラメータ算出部２３、及び画像生成部２５を備える。流体パラメータ算出部２３は、画素（ＤＴ画像のボクセル）毎に、流速の分散を算出する分散算出部２３１と、分散算出部２３１が算出した分散と、算出した分散と流速の分布（ボクセル内分布の形状）の推定モデルとを用いて流速の微分ベクトルを算出する微分算出部２３３と、微分ベクトルを用いて個々の流体パラメータを算出する各種パラメータ算出部とを有する。流体パラメータには、例えば、ＷＳＳ、ＥＬのほか、ＰＧ（圧勾配：pressure gradient）、ＯＳＩ（振動せん断インデックス：Oscillatory shear index）などがあるが、図２では代表として、ＷＳＳ及びＥＬを算出するＷＳＳ算出部２３５、ＥＬ算出部２３７のみを示している。

【0023】

ＤＷＩの計測データから疑似拡散テンソルを算出するＤＴＩ算出部２１の機能は、従来のＤＴＩと同様であり、よく知られているが、本実施形態のＭＲＩ装置は、流体パラメータ算出部２３が、ＤＴＩで得られる疑似拡散テンソルを利用して、ＰＣ法では得られない精度の高い流体パラメータを算出することが特徴である。具体的な流体パラメータ算出については後述の実施形態で述べることとし、以下では、図３を参照して各実施形態に共通する処理の概略を説明する。

【0024】

まずシーケンサ１７の制御のもとで、計測部１０がＤＷＩのパルスシーケンスを実行し、計測データを収集する（Ｓ１）。ＤＷＩのパルスシーケンスとしては、一般的にＥＰＩシーケンスが用いられるが、マルチショットＥＰＩやＳＳ（シングルショット）ＦＳＥなどのパルスシーケンスを用いてもよい。図４に典型的なＤＷＩのパルスシーケンスを示す。図中、401、403はＲＦパルス、402、404はスライス傾斜磁場（Gs）、405は位相エンコード傾斜磁場（Gp）、406はリードアウト傾斜磁場（Gr）、407はNMR信号である。

【0025】

図示するように、ＤＷＩのパルスシーケンスは、励起用ＲＦパルス401で励起した磁化を反転させる180度ＲＦパルス403の印加前後に、流体の位相拡散を促す傾斜磁場パルス（ＭＰＧパルス）408、409を印加することが特徴であり、このＭＰＧパルス408、409の強度を異ならせることで、ボクセル内での流体の速度の違い（速度分散）に応じた位相分散を生じさせ、結果としてボクセルの信号強度を低下させる。言い換えれば、この信号強度の低下によって、ボクセル内での流体の速度分散（疑似拡散）の情報を得ることができる。すなわち、流体を対象にした場合、ＤＷＩにより流速分散の情報を得ることができる。ＭＰＧパルスの強度の指標はｂ値と呼ばれ、速度の早い流体について流体パラメータを求める場合には、比較的低いｂ値とすることが好ましい。ｂ値は、予め所定の範囲内で複数のｂ値を設定しておいてもよいし、ユーザが特定のｂ値を設定することも可能である。

【0026】

またＭＰＧパルスを印加する軸を異ならせることで拡散および疑似拡散の異方性に関する情報が得られる。図４は、リードアウト方向（Ｇｒ軸）にＭＰＧパルスを印加する例を示しているが、拡散テンソルイメージング（ＤＴＩ）では、ＤＴを算出するために６方向以上のＭＰＧの印加が必要であり、予め設定された複数の軸について図４と同様に計測を行う。

【0027】

ＥＰＩでは、1回の励起用ＲＦパルス後に、位相エンコード傾斜磁場パルス405及び反転するリードアウト傾斜磁場パルス406を印加することで、1枚の画像を構成する複数のＮＭＲ信号407を得ることができる。ＮＭＲ信号407は、印加された位相エンコード量とリードアウト傾斜磁場で決まるｋ空間の位置に配置され、デジタルデータ（計測データ）として、計算機２００に渡され、ここで演算部２０が種々の演算を行う。

【0028】

演算部２０は、まずＤＴ算出部２１において、ボクセル毎（画素毎）の拡散テンソルを算出する（Ｓ２）。

【0029】

拡散テンソルＤＴの算出方法は周知であるが、簡単に説明すると拡散テンソルＤＴは、ｂベクトルを変化させて撮影した複数のＤＷＩの信号強度から、次式（１）で表される対応する直線の勾配を求めることで決定することができ、３×３の行列で表される。

【0030】

【数1】

式（１）中、Ｓｉは所定ｂベクトルの際に得られた信号強度、ｂはｂベクトルの絶対値、Ｇｉはｂベクトルの傾斜磁場印加方向を表す単位ベクトルである。なお、ｂ値を多次元の数値、すなわちｂベクトルと同義とする記述方法もある。

【0031】

この拡散テンソルＤＴの対角化を行うことで固有値と固有ベクトルが得られる。固有値と固有ベクトルに変換することで、ボクセルにおける拡散テンソルの主軸とそれに直交する２軸が求められることになる。また固有ベクトル対角要素を用いてみかけの拡散係数ＡＤＣ等を算出することができる。

【0032】

拡散テンソルは、図５（下側）に示すような楕円体モデルで表され、ボクセル内の流れが、例えば、同図上側の左端に示すような同じ速度で流れる（１）プラグ流の場合には分子拡散のみを反映した小さな真球状となり、中央の３つのように（２）、（３）層流の場合には、異方性を持つ楕円体状となり、等方的にランダムな流れがある場合（４）には大きな真球状となる。ＤＴＩは、図５の上側に示すような流体の信号値から、拡散テンソルを求めて下側のような画像を得るものであるが、ＤＴＩから流れの状態を示す指標である流体パラメータを直接算出することはできない。

【0033】

本実施形態の流体パラメータ算出部２３は、ボクセル内の流速の分布に仮定を置くことで、その分散と流速微分ベクトルを算出し、流速微分ベクトルを用いた流体パラメータの算出を可能にするものである。

【0034】

このため流体パラメータ算出部２３は、まず、分散算出部２３１が、ＤＴ算出部２１が算出した各ボクセルの拡散テンソルＤＴを用いて流速分布ρの分散を算出する（Ｓ３）。次に仮定した流速分布ρの形状と、算出した分散とを用いて流速の微分値を算出する（Ｓ４）。次に個別のパラメータ算出部（２３５、２３７）が、流速の微分値を用いて、ＷＳＳやＥＬを算出する（Ｓ５）。

【0035】

画像生成部２５は、ＤＷＩで得られる形態画像やＤＴＩ（拡散テンソル画像）とともに、算出したパラメータの値を示す表示画像を生成する（Ｓ６）。表示されるＤＴＩの一例を図６に示す。図６の上側は脳画像のサジタル断面を示す図で、下側は、脳画像のうち四角で囲った部分のＤＴＩを示す図である。図示するように、ボクセル単位でＤＴＩが表示されることで、脳脊髄液が流れる壁面が描出され、また壁せん断応力を算出すべき法線方向を提示或いは把握することができる。

【0036】

以上概略を説明したように、本実施形態のＭＲＩ装置は、ボクセル内の流体ベクトルの分布の形状に推定した形状を用いることで、ＤＷＩ計測により得た拡散ベクトルからベクトルの微分値を算出することができ、この微分値を用いて精度よく流体パラメータを算出することができる。

【0037】

以下、演算部２０による処理の具体的な実施形態を説明する。

【0038】

＜実施形態１＞
本実施形態では、流体の流れが線形流であることを前提とし、また流速の分布のモデルとして矩形分布を用いる。

【0039】

線形流とは、流体の輸送距離が時間経過に比例する流れである。この線形流に、流体の位相を拡散させる一対の傾斜磁場（図４のＭＰＧパルス408,409：これらは反転180度403を挟んでいるので正負反転している）を印加すると、図７に示すように、拡散時間τｄは、パルスの印加時間をδ、1つ目のパルスから２つ目のパルスまでの時間をΔとすると、
τｄ＝Δ-δ/３ （２）
で表される。

【0040】

分散算出部２３１は、この拡散時間τｄを用いて、各ボクセルの流速分布ρ（ｖ、ｘi）の分散Cov(ρ)を算出する。式（１）に示したＤＴ算出部２１が算出したＤＴ（３×３の行列）について、ｂを０に近づけたときの極限をとったときのテンソルをＰとすると、Ｐは、流速分布ρの分散Cov(ρ)と次式の関係があることがわかっている（非特許文献１：式（３））。

【0041】

【数3】

【0042】

さらに条件によっては次のように単純化できる。
（１）Cov(ρ)がＤよりも十分大きい場合

【数31】

【0043】

（２）さらにΔ>>δのときにはτ_ｄ＝Δの近似が可能なため

【数32】

【0044】

分散算出部２３１は、この式（３）或いは（３１）又は（３２）を用いて、各ボクセルの流速分布ρ（ｖ、ｘi）の分散Cov(ρ)を算出する。なお、式（３）を用いる場合には疑似拡散テンソルＰから拡散係数（テンソル）Ｄを減算する処理を行うが、先見情報を用いる方法と追加計測データを用いる方法の二通りがある。先見情報を用いる方法では、予め知られている情報である自由水の拡散係数（例えば生体温度３７℃における自由水の拡散係数3×10＾(-9) m^2/s）をＰから減算する。追加計測データを用いる方法には、さらに二通りがある。一つは、τｄを変更しながら複数のＰを計測し、それをフィッティングすることでCov(ρ)とＤを同時に求める方法である。もう一つは、速度成分を打ち消すようなバイポーラーパルスの組で構成されるＭＰＧを用いることでＤのみを計測しＰから減算する。分散Cov(ρ)算出の具体的な流れについては以下のようになる。ここでは簡単のため主に２次元で説明する。

【0045】

図５上側（２）に示す層流（平行）の場合、一般性を与えるために、適当な回転座標変換によりｘ方向の層流と考えると、
∇Vxx = ∇Vyx = ∇Vyy = 0，∇Vxy = a （aは所定の値、∇Vijは流速Vのi方向の流速成分をj方向で微分した値を表す）
及び
Cov(ρ_xρ_x) = a^2, Cov(ρ_xρ_y) = Cov(ρ_yρ_x) = Cov(ρ_yρ_y) = 0
が成立する。

【0046】

なお３次元の場合には、y方向とz方向が不定となるため、その割合を決定する必要がある。例えば、近傍のCovの変化率を算出し、その比率に応じてｙ方向及びｚ方向に値を割り振る。またＤＴＩで位相情報が得られている場合には、位相から各ボクセルの平均速度を計算し、近傍との速度差の二乗で割り振る、などの手法が取りえる。なお、位相から各ボクセルの平均速度を計算する方法は、ＰＣ法と同じである。また、後述のＰＣ法を併用している場合にも、各ボクセルの平均速度を計算し、近傍との速度差の二乗で割り振るなどの手法が取りえる。また通常のＤＴＩから計算される流速微分ベクトルでは正負の符号が判らないが、ＤＴＩもしくは併用したＰＣの符号を適用することで、正負の符号を付与することが可能である。

【0047】

また図５中央（３）に示す層流（線形拡大）の場合は、流れをある点から湧き出す流体の一部（半径r，角度θ）と捉えて、極座標系に変換する。非圧縮性の液体では、流量と断面積の関係から、次の式が成立する。
∇Vrr = (1/(2πr))’ = -1/(2πr²) （2次元の場合）
∇Vrr = (1/(4πr²)’ = - 1/(2πr³) （3次元の場合）
∇Vrθ = ∇Vθr = ∇Vθθ = 0
Cov(ρ_rρ_r) = 1/(4π²r⁴)（2次元の場合）
Cov(ρ_rρ_r) = 1/(4π²r⁶)（3次元の場合）
Cov(ρ_θρ_θ) = r²θ², Cov(ρ_rρ_θ) = Cov(ρ_θρ_r) = 0

【0048】

3次元の場合には、層流（平行）の場合と同様に、y方向とz方向が不定となるので、その割合を決定する。例えば、近傍のCovの変化率を算出し、その比率に応じて割り振る、ＤＴＩで位相情報が得られている場合には、位相から各ボクセルの平均速度を計算し、近傍との速度差の二乗で割り振るなどである。また、後述のＰＣ法を併用している場合にも、各ボクセルの平均速度を計算し、近傍との速度差の二乗で割り振るなどの手法が取りえる。

【0049】

次に微分算出部２３３が、流速分布の形状について推定を置き、算出した分散Cov(ρ)を用いて流速微分ベクトルを算出する。ここで、流速分布ρの形状が矩形分布である、即ち図８に示すように、ボクセル内で空間的に線形な分布であると仮定すると、矩形分布の分散と流速微分ベクトル∇Ｖは、次式（４）の関係となり、流速微分ベクトル∇Ｖは式（５）で表される。

【0050】

Cov(ρ) = (vmax-vmin)²/12 = (∇V・Δx)²/12 （４）
∇V ＝ √(12 Cov(ρ))/Δx （５）
式（４）、（５）において、Δxはボクセルサイズである。

【0051】

微分算出部２３３は式（５）により流速微分ベクトルを算出する。ここで、bが十分小さい場合に計測したＤＴは、Ｐとみなすことができるので、低ｂ値で得たＤＴを式（３）のＰとして、式（５）に代入すると、式（５）は
∇V = √(24 (Ｐ-Ｄ)/τd)/Δx （５１）
となる。式（３１）で示したように、ＰがＤよりも十分大きいときにはＤを無視ことができる。算出される∇Vの符号は判らないが、∇Vを用いることで符号に依存しない流体パラメータを算出することができる。

【0052】

流体パラメータ算出部２３は、こうして算出した流速微分ベクトルを用いて各種の流体パラメータを算出する。

【0053】

例えば、ＷＳＳ算出部２３５は、まず、図６に示したように、流体が流れる壁面を抽出する。壁面の抽出には、ＤＴＩのＭＤ（平均拡散係数：mean diffusivity）に閾値を設定する方法、ＭＲＡ（ＭＲアンジオグラフィ）や強度のＴ２Ｗ撮像など液体を強調する画像を別途計測しておき輝度に閾値を設定する方法が考えられる。次いで、壁面に対する法線ベクトルｎを計算し、法線方向における流速微分値を計算する。ＷＳＳは、この法線方向の流速微分値∇V・nを用いて、次式により算出することができる。

【0054】

【数6】

式（６）において、ＷはＷＳＳ、ｗは流速微分値∇V・nを表す。μは流体の粘度（固定値）（例えば、血液であれば、μ＝0.004Pa*s）である(以下、同じ)。

【0055】

またＥＬ算出部２３７は、対象とする領域Ｖを指定し、次式（７）により、エネルギーロス（ＥＬ）を算出する。

【数7】

【0056】

またＥＬ算出部２３７は、上記式（７）と数学的に等価な別手法を用いることも可能である。すなわち、ＤＴＩの各要素、もしくは固有値を直接使用し、流速微分ベクトルを暗に使用することが可能である。例えば、上記式（７）を対角化を行う適切な座標変換を施したのち展開して得られる一要素∇Vi,j^2はＤＴＩの一要素Pi,jと定数倍を除いて等価となる。ここで定数倍は、速度分布として矩形分布やガウス分布などを仮定すれば一意に定まる。このように、ＤＴＩで得られた情報から、流速微分ベクトルを明に求めずに、ＥＬを計算してもよい。また図２では、ＷＳＳ算出部とＥＬ算出部のみを例示したが、流体パラメータとして、ユーザが指定した方向nに関するＰＧ（圧勾配）を算出してもよい。ＰＧの算出方法は∇V・nで表される。

【0057】

こうして算出した流体パラメータは、図６に示したＤＴＩとともに表示装置４０に表示するなどユーザに提示する。これによりユーザは、ＤＴＩとともに流体の異常に関する指標を得ることができる。

【0058】

本実施形態によれば、ＤＴＩのデータ（疑似拡散テンソル）と、流体の流速分布についての推定モデルを用いることで、流速の微分ベクトルを求めることができ、それを用いて流体パラメータを算出することができる。ここで、微分ベクトルを算出するための流速分布の分散は、ＤＴＩにより直接計測可能なため、ＰＣ法のような空間分解能の影響を受けにくい。従って、他方（例えばＰＣ法）と同じ空間分解能であれば精度を向上することができ、空間分解能を下げられるため計測時間を短縮することが可能となる。

【0059】

＜実施形態１の変形例１＞
上述した実施形態１では、流体パラメータとして、ＷＳＳ及びＥＬを用いる場合を説明したが、ＤＷＩ計測時に体動同期撮影を行い、体動に関する情報を得ている場合には、その情報を加えて、さらに別の流体パラメータを算出することも可能である。

【0060】

体動同期撮影には、よく知られているように、心臓の拍動もしくは心電図に同期して計測を行う心拍同期撮像もしくは心電同期撮像や、呼吸動に同期した呼吸同期撮像、さらには両方の体動のタイミングを用いる同期撮像があり、そのいずれにも対応できる。また同期撮像の手法として、同期信号をトリガーとして計測を行う方法と、連続して計測したデータに対し、同期信号を用いて事後的に解析するレトロスペクティブ解析方法とがあるが、同期信号の情報が得られるものであれば、特に限定されない。

【0061】

同期信号を用いた流体パラメータとしては、例えば、ＴＡＷＳＳ（時間平均ＷＳＳ）やＯＳＩ（振動せん断インデックス）があり、ＴＡＷＳＳは、各時相で計算されたＷＳＳをW(t)とおくと、次式（８）により算出することができる。

【数8】

【0062】

またＯＳＩは次式（９）より算出することができる。

【数9】

【0063】

＜実施形態１の変形例２＞
実施形態１では、流速の分布の推定モデルを矩形分布としてが、推定モデルは矩形分布に限らず、流体が流れる臓器の形状やサイズを考慮して変形が可能であり、流速の分布を示す関係式を変更することで、実施形態１と同様に流体パラメータを算出することができる。

【0064】

一例として、図８（ｂ）に示すような正規分布を用いる場合を説明する。流速確率密度関数ρ（ｖ）を正規分布±3σで切り捨て)で近似し、ボクセル内で空間的に線形に分布と仮定すると同分布の分散と流速微分は次の関係をもつ。
Cov(ρ) = σ² = (∇V・Δx)²/6 （１０）
ここでΔxはボクセルサイズである。

【0065】

従って、上述した式（５）は、次式（５２）となる。
∇V = √(6 Cov(ρ))/Δx （５２）
式（５２）は、矩形分布の場合と同様に、ＤＴＩで計測したＰ（＝ｂが十分小さいときのＤＴ）を用いて、次式（５３）となり、
∇V = √(12 (Ｐ-Ｄ)/τd)/Δx （５３）
微分値∇Ｖを算出することができる。本変形例でも、ＰがＤよりも十分大きいときにはＤを無視してよい。

【0066】

その後、算出した微分ベクトルを用いてＷＳＳやＥＬなどの流体パラメータを算出することは実施形態１と同様である。

【0067】

＜実施形態１の変形例３＞
本変形例では、拡散強調イメージングのパルスシーケンス自体を利用して、フェイズコントラスト（速度に依存する位相シフト）を取得し、フェイズコントラストから求めた流速を低ｂ値のＤＴＩと併用もしくは代替として用いる。

【0068】

即ち、計測部１０は、実施形態１と同様にＤＷＩのパルスシーケンスを実行して複素数値の計測データを取得する。演算部２０は、複素数値の計測データの位相情報を用いて、検査対象内の各ボクセルにおける流体の平均流速を算出する。流体パラメータ算出部２３は、ＤＴ算出部２１が算出した疑似拡散テンソルと、フェイズコントラストから算出した流速と、流体の流速分布の推定モデルとを用いて、流速微分ベクトル∇Vを算出する。具体的には二つの方法がある。

【0069】

一つの方法は、前述したように疑似拡散テンソルだけでは不定になる流速微分ベクトルの方向割合を、フェイズコントラストから算出したボクセル毎の流速から計算される周囲との差分で決定する方法である。例えば、流速の差分の二乗で割合を決定する。もう一つの方法は、フェイズコントラストから算出した流速を用いて、周囲ボクセルとの差分を∇ＶＰとし、疑似拡散テンソルから上述の方法で求められた流速微分ベクトルを∇ＶＤとし、両者の平均や加重平均、もしくはフェイズコントラストで求められる流速の閾値でどちらかを選択する。なお、後者の方法については、別途フェイズコントラスト用の計測を併用する実施形態３で詳細に説明する。

【0070】

その後、流速微分ベクトルを用いて、流体パラメータを算出することは実施形態１と同様である。

【0071】

＜実施形態２＞
実施形態１では、線形流を前提として、流速の微分ベクトルを算出したが、本実施形態では、線形流と拡散との中間的な複雑な流れの流体パラメータを算出する。線形流と拡散とが混じった流れは、図９に示すように、例えば、ボクセル内で渦がある場合や、流れが複雑に屈曲している場合などで、輸送能が低下し、輸送距離と時間とが比例しない場合である。

【0072】

本実施形態では、異常拡散に倣って、次式（３３）により、時間変化をモデル化した指標αを用いて、拡散時間（τｄ）のα乗「τｄ^α」（０≦α≦１）と分散Cov(ρ)を算出する。
Ｐ＝τｄ^α/２*（Cov（ρ））+Ｄ （３３）

【0073】

ここでα＝０の場合は拡散、α＝１の場合は線形流を意味する。

【0074】

式（３３）をもとに算出した分散を用いて、微分∇Ｖを算出すること、さらに流体パラメータを算出すること、及び、微分の算出において、流速分布の推定形状として矩形形状或いはガウス分布形状を用いてもよいことは、実施形態１或いはその変形例と同様である。

【0075】

指標αの計測および計算手法を説明する。この場合の流体パラメータ算出部の機能ブロック図を図１０に示す。図示するように、本実施形態では、拡散指標算出部２３９が追加される。その他の要素は、図２と同様であり、同じ符号で示す。

【0076】

拡散要素αを推定する一つの手法として、拡散傾斜磁場（ＭＰＧ）の間隔を異ならせた複数回（例えば３回）の計測を行う。１回目の間隔をΔ１、２回目の間隔をΔ２、３回目の間隔をΔ３・・としたとき、それぞれで求められる流体の流速が間隔に対し線形に変化する場合には、線形流とみなすことができ、線形性からの外れ方によりαを算出する。

【0077】

その後、式（３）を式（３３）に置き換える以外は実施形態１と同じである。なお、本説明では指標αを計測により求める場合を説明したが、生体モデルとして指標αを固定して流体パラメータを求めても良い。典型的には、α＝１と固定した場合が線形流のモデルであり、α＝０と固定した場合が拡散のみのモデルである。

【0078】

＜実施形態３＞
一般に拡散テンソルイメージングは、流速が遅い流体については精度高くＤＴＩ値（疑似拡散テンソル）を求められることが知られているが、流速が早い流体では拡散により信号が減衰しきってしまうため、精度が低下する。一方、ＰＣ法は流速が遅い場合には、計測で求められる位相差が小さいため、精度が低下する。

【0079】

本実施形態では、両者の特性を生かし、２つの計測を択一的或いは合わせて行うことで、広い速度範囲で流体パラメータ算出の精度を向上する。

【0080】

本実施形態の演算部の機能ブロックを図１１に示す。図１１において、図２に示す要素と同じ機能を持つ要素については同じ符号で示し、重複する説明は省略する。

【0081】

図１１に示すように、本実施形態の演算部２０は、ＤＴＩを用いた流体パラメータ算出部２３とは別に、ＰＣ法で得られた計測データを用いて流体パラメータを算出する第２の流体パラメータ算出部２３Ｂを備えている。また２つの流体パラメータ算出部の結果を統合するパラメータ統合部２９を備えた構成としてもよい。

【0082】

第２の流体パラメータ算出部２３Ｂは、図示しないが、ＰＣ法の計測データを用いて流速を算出する流速算出部２３１Ｂ、流速の微分ベクトルを算出する微分算出部２３３Ｂなどを備え、微分算出部２３３Ｂが算出した流速の微分ベクトルを用いて、上述した式（６）、（７）等を用いてＷＳＳやＥＬ等の流体パラメータを算出する。

【0083】

流体パラメータ算出部２３の処理とするか、第２の流体パラメータ算出部２３Ｂの処理とするかは、流体の速度について所定の閾値を予め設定しておく或いは対象とする流体が血液か脳せき髄液か、動脈流か静脈流か、などに応じて予め設定しておいてもよいし、ユーザが表示装置４０に表示されるＧＵＩなどを介して選択し、演算部２０に設定してもよい。

【0084】

このようなＧＵＩの例を図１２（A）、（B）に示す。図１２に示す例では、撮像条件を設定する画面1200において、流体パラメータを算出することがデフォルトで設定されているか、流体パラメータの算出の「要」とのユーザ指示を受け付けると、ＰＣ法による流体パラメータ算出とする流体の部位、種類、流速などの情報を択一的に受け付けるＧＵＩが表示される。

【0085】

このＧＵＩにより、例えば所定の流速が設定されると、計測部１０は、流速が予め定めた流速より早い場合には、ＤＷＩ撮像ではなくＰＣ法の撮像を行い、設定された流速が予め定めた流速より遅い場合には、ＤＷＩ撮像を行う。

【0086】

計測部１０がＰＣ法の撮像を行った場合には、第２の流体パラメータ算出部２３Ｂは、従来法（例えば特許文献１などに記載される技術）により流体パラメータを算出する。また計測部１０がＤＷＩ撮像を行った場合には、実施形態１や実施形態２と同様に、ＤＴ算出部２１がＤＴＩを作成し、流体パラメータ算出部２３がＤＴＩの情報を用いて流体パラメータを算出する。

【0087】

以上は、ＤＷＩ撮像かＰＣ撮像かを択一的に行う例であるが、両者を併用し、その結果を重み付け加算してもよい。その場合には、計測部１０はＤＷＩ撮像とＰＣ撮像とをそれぞれ実行する。撮像の順序は限定されない。演算部２０は、両撮像で得たデータを、それぞれ、ＤＴ算出部２１と第２の流体パラメータ算出部２３Ｂに渡し、流体パラメータ算出部２３による処理と、第２の流体パラメータ算出部２３Ｂによる処理とを平行して行う。

【0088】

流体パラメータの算出において、通常のＤＴＩから計算される流速微分ベクトルでは正負の符号が判らないが、ＰＣの符号を適用する。これにより算出のために正負の符号が必要な流体パラメータをも算出することができる。

【0089】

その後、パラメータ統合部２９は、両方の流体パラメータ算出部２３、２３Ｂの処理結果を次式のように重み付け加算する。
ＷＳＳ＝ω1・ＷＳＳ_ＤＷ＋ω2・ＷＳＳ_ＰＣ
式中、ＷＳＳ_ＤＷは流体パラメータ算出部２３の算出結果、ＷＳＳ_ＰＣは流体パラメータ算出部２３Ｂの算出結果である。ここではＷＳＳの例を示したが、他の流体パラメータについても同様である。

【0090】

重みω1、ω2は、流体の種類や速度に関わらず、０．５、０．５等の定数にしておいてもよいが、流体の種類や速度に応じて予め所定の重みを決めておき、例えば、ＧＵＩを介してユーザが入力した流体の種類や速度、或いは流速算出部２３１Ｂが算出した流体の速度に対応する重みを採用してもよい。

【0091】

本実施形態によれば、流体パラメータを算出するための基礎データの計測方法として、特徴の異なる二つの撮像を採用し、それぞれの算出結果を組み合わせることで、幅広い速度や部位に対応して精度のよい流体パラメータを算出することができる。

【符号の説明】

【0092】

１：ＭＲＩ装置、２：画像処理装置、１０：計測部、１１：静磁場発生部（マグネット）、１２：プローブ（送受信ＲＦコイル）、１３：傾斜磁場コイル、１４：受信器、１５：高周波磁場発生器、1６：傾斜磁場電源、１７：シーケンサ（計測制御部）、２０：演算部、２１：ＤＴ算出部、２３：流体パラメータ算出部、２３Ｂ：第２流体パラメータ算出部、２５：画像生成部、３０：入力装置（入力部）、４０：表示装置、５０：記憶媒体、２００：計算機、２３１：分散算出部、２３３：微分算出部、２３５：ＷＳＳ算出部、２３７：ＥＬ算出部、２３９：拡散指標算出部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】