TOP > 国内特許検索 > パルスドプラ計測装置、その方法及びそのプログラム > 明細書

明細書 :パルスドプラ計測装置、その方法及びそのプログラム

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4729765号 (P4729765)
公開番号 特開2007-215816 (P2007-215816A)
登録日 平成23年4月28日(2011.4.28)
発行日 平成23年7月20日(2011.7.20)
公開日 平成19年8月30日(2007.8.30)
発明の名称または考案の名称 パルスドプラ計測装置、その方法及びそのプログラム
国際特許分類 A61B   8/06        (2006.01)
G01S  15/50        (2006.01)
FI A61B 8/06
G01S 15/50
請求項の数または発明の数 5
全頁数 20
出願番号 特願2006-040587 (P2006-040587)
出願日 平成18年2月17日(2006.2.17)
審査請求日 平成20年10月24日(2008.10.24)
特許権者または実用新案権者 【識別番号】599016431
【氏名又は名称】学校法人 芝浦工業大学
発明者または考案者 【氏名】田中 直彦
個別代理人の代理人 【識別番号】100064414、【弁理士】、【氏名又は名称】磯野 道造
【識別番号】100111545、【弁理士】、【氏名又は名称】多田 悦夫
審査官 【審査官】宮川 哲伸
参考文献・文献 特開平2-203849(JP,A)
特許第2953083(JP,B2)
特開平3-228751(JP,A)
特開平2-215445(JP,A)
特開2005-46527(JP,A)
特開平4-197249(JP,A)
調査した分野 A61B 8/06
G01S 15/50
特許請求の範囲 【請求項1】
対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測装置であって、
前記超音波を前記対象物に送波する送波手段と、
前記対象物からの前記超音波に対する反射波を受波する受波手段と、
この受波手段で受波した反射波から、周波数毎の位相を検出するスペクトル解析手段と、
このスペクトル解析手段で検出された周波数毎の位相について、前記周期毎に位相差を演算する位相差演算手段と、
前記周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで2πの整数倍の位相を加算又は減算することで前記位相差を修正する位相差修正手段と、
前記周波数と当該周波数における位相差とを座標点とする座標系において、前記複数の位相修正パターン毎に、前記周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する直線近似手段と、
この直線近似手段で演算された前記誤差が最小となる直線の傾きを、前記移動速度に対応する位相差の傾きとして選択する速度特定手段と、
を備えていることを特徴とするパルスドプラ計測装置。
【請求項2】
前記位相差演算手段は、前記周波数毎の位相差を、予め定めた送受波回数毎に平均化することを特徴とする請求項1に記載のパルスドプラ計測装置。
【請求項3】
前記スペクトル解析手段で検出された位相に基づいて、クラッタ成分を除去するフィルタリング手段を備えていることを特徴とする請求項1又は請求項2に記載のパルスドプラ計測装置。
【請求項4】
対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測方法であって、
前記超音波を前記対象物に送波し、前記対象物からの反射波を受波する送受波ステップと、
前記反射波から、周波数毎の位相を検出するスペクトル解析ステップと、
このスペクトル解析手段で検出された周波数毎の位相について、前記周期毎に位相差を演算する位相差演算ステップと、
前記周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで2πの整数倍の位相を加算又は減算することで前記位相差を修正する位相差修正ステップと、
前記周波数と当該周波数における位相差とを座標点とする座標系において、前記複数の位相修正パターン毎に、前記周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する直線近似ステップと、
この直線近似ステップで演算された前記誤差が最小となる直線の傾きを、前記移動速度に対応する位相差の傾きとして選択する速度特定ステップと、
を含んでいることを特徴とするパルスドプラ計測方法。
【請求項5】
対象物に対して一定の周期で送波された超音波に対する反射波に基づいて、前記対象物の移動速度を計測するために、コンピュータを、
前記反射波から、周波数毎の位相を検出するスペクトル解析手段、
このスペクトル解析手段で検出された周波数毎の位相について、前記周期毎に位相差を演算する位相差演算手段、
前記周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで2πの整数倍の位相を加算又は減算することで前記位相差を修正する位相差修正手段、
前記周波数と当該周波数における位相差とを座標点とする座標系において、前記複数の位相修正パターン毎に、前記周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する直線近似手段、
この直線近似手段で演算された前記誤差が最小となる直線の傾きを、前記移動速度に対応する位相差の傾きとして選択する速度特定手段、
として機能させることを特徴とするパルスドプラ計測プログラム。
発明の詳細な説明 【技術分野】
【0001】
本発明は、対象物に送波した超音波に対する反射波に基づいて、移動する対象物の速度を計測するパルスドプラ計測装置、その方法及びそのプログラムに関する。
【背景技術】
【0002】
従来、人体の診断を行う装置として、生体内の血流の速度を計測するパルスドプラ法を用いた診断装置が種々開発され、臨床の場で広く利用されている。これらの診断装置は、生体内の血流に対して超音波のパルスを一定周期で繰り返し送波し、計測対象である赤血球が動くことにより生じる反射波の位相変化を検出することで、血流の速度を計測している。
ここで、図9を参照して、従来のパルスドプラ法について、その概念を説明する。
パルスドプラ法は、図9(a)に示すように、送波器Sから、一定周期で超音波(送波パルス)を、動きを伴う対象物(例えば、血流)Bに送波し、受波器Rで、対象物Bからの反射波を受波する。このとき、1番目の送波と2番目の送波との間に対象物BがΔxだけ移動したとすると、図9(b)に示すように、伝播時間差Δtだけ受波波形がシフトする。
この移動量Δxと伝播時間差Δtとの関係は、伝播媒質中の音速をc〔m/s〕とすると、以下の(1)式で表される。
【0003】
【数1】
JP0004729765B2_000002t.gif

【0004】
ここで、送波パルスの繰り返し周期をT〔s〕とし、対象物Bの速度をv〔m/s〕とすると、前記(1)式は、以下の(2)式となる。
【0005】
【数2】
JP0004729765B2_000003t.gif

【0006】
この伝播時間差Δtは、図9(c)に示すように、シフトした受波のスペクトルにおいて位相のずれとなって現れる。すなわち、伝播時間差Δtは、図9(d)に示すように、シフトした受波のスペクトルの位相差(位相回転)に比例する。
ここで、送波パルスの中心周波数をω、ωにおける位相差をΔθとすると、図9(d)の直線の傾きはΔθ/ωとなる。これが伝播時間差Δtに相当するため、前記(2)式より、以下の(3)式を得ることができる。
【0007】
【数3】
JP0004729765B2_000004t.gif

【0008】
このように、対象物Bの速度vは、送波パルスの中心周波数ωと位相差Δθとにより求めることができる。
そして、このパルスドプラ法を用いた診断装置は、血流の速度を表示装置上にカラー表示させることにより、心臓、肝臓等の診断、血行動態のリアルタイムの診断等、多くの診断に利用されている。
しかし、パルスドプラ法を用いた診断装置では、前記(3)式に示すように、血流計測を反射波の位相変化により行っているため、図9(e)に示すように、血流の速さによっては、位相の折り返し(いわゆるエイリアシング)が発生し、正しい計測結果が得られないという問題がある。
【0009】
この問題を解決する手法として、種々の提案がなされている。例えば、2つの異なる周波数をもつパルスを予め定めた周期で対象物に送波し、それぞれ送波したパルスと反射波との位相差の差をとることで、位相差が±πを超えエイリアシングが発生する場合でも、2つの位相差の差が±π以内でありエイリアシングが発生しないことを利用して、対象物の速度を計測する手法(以下、2周波法という)が提案されている(特許文献1、特許文献2参照)。
また、別の手法として、パルスを異なる周期で対象物に送波し、それぞれ送波したパルスと反射波との位相差をとり、周期の差により発生する位相差の差を利用することで、対象物の速度を計測する手法(以下、2周期法という)が提案されている(特許文献3参照)。
【0010】
また、これらの従来の手法では、送受波を複数回行い位相の平均化を行うことで、計測速度の精度を高めている。さらに、従来の手法では、送受波するパルスの中心周波数により位相の検出を行っているが、パルスの帯域を狭くする方が、S/Nの点で有利であると考え、比較的パルス長の長いバースト波を用いている。
また、従来の手法では、実際の血流計測において、血管壁、骨等の動きのない反射体からの反射波を拾ってしまうクラッタ成分を除去するため、反射波の位相の差に基づいて、動きのない反射体からの反射波の信号を除去するMTI(Moving Target Indication)フィルタを使用している。

【特許文献1】特許第2953083号公報(段落0009)
【特許文献2】米国特許第4534357号明細書
【特許文献3】特開平4-197249号公報(第3-5頁、第3-5図)
【発明の開示】
【発明が解決しようとする課題】
【0011】
前記した従来の手法(2周波法、2周期法)では、S/N改善の観点から、受波した反射波の位相の平均化を行ったり、比較的パルス長の長いバースト波を使用したりすることで、計測の精度を高めようとしていた。
しかし、実際の生体内の血流計測において、精度を高めることができない原因は、各々の赤血球からの反射波の干渉による位相の変動であることを、本発明者は解明した。
【0012】
概略すると、生体内の血流を計測する場合、超音波の反射体となるのは赤血球である。この赤血球は、超音波の一波長内に極めて多数かつランダムに存在し、多数散乱体として振舞う。このような多数散乱体に対して、超音波のパルスを送波すると、各々の赤血球からの反射波は互いに重なり合い干渉することになる。その結果、血流が一定の速度で流れていたとしても、反射波の位相は不規則に変動することになる。
これに対し、従来の手法(2周波法、2周期法)は、干渉による位相変動の影響を考慮していないため、計測の精度を高めることができなかった。すなわち、従来の手法においては、位相差の差を求めるため、位相の変動による影響を強く受けることになり、エイリアシングの問題は解決できても、計測の精度を逆に劣化させることとなっていた。
【0013】
さらに、実際の血流計測において、従来の手法は、MTIフィルタによりクラッタ成分を除去している。しかし、従来の手法は、送受波するパルスの中心周波数により位相の検出を行っているため、送波周期の整数倍に相当する距離だけ移動した反射体については、位相差がゼロになってしまい、実際には対象物が動いているにも拘らず、クラッタ成分として除去されてしまう。このように、従来の手法では、求められない速度(以下、ブラインド速度という)が存在するという問題があった。
【0014】
本発明は、以上のような課題を解決するためになされたものであり、対象物の速度を超音波パルスの反射波の位相差により計測する際に、位相差における位相の折り返し(エイリアシング)による誤計測を防止するとともに、対象物の速度の計測精度を高め、さらに、ブラインド速度をなくすことを可能にしたパルスドプラ計測装置、その方法及びそのプログラムを提供することを目的とする。
【課題を解決するための手段】
【0015】
本発明は、前記目的を達成するために創案されたものであり、まず、請求項1に記載のパルスドプラ計測装置は、対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測装置であって、送波手段と、受波手段と、スペクトル解析手段と、位相差演算手段と、位相差修正手段と、直線近似手段と、速度特定手段と、を備える構成とした。
【0016】
かかる構成において、パルスドプラ計測装置は、送波手段によって、超音波のパルスを、速度を計測する対象となる対象物に送波し、受波手段によって、対象物から反射された反射波を受波する。このとき、対象物が移動していると反射波の位相がずれることになる。
また、パルスドプラ計測装置は、スペクトル解析手段によって、反射波から、周波数毎の位相を検出する。すなわち、スペクトル解析手段は、送波したパルスにおける1つの周波数だけでなく、予め定めた複数の周波数における位相を検出する。
【0017】
さらに、パルスドプラ計測装置は、位相差演算手段によって、スペクトル解析手段で検出された周波数毎の位相について、周期毎に位相差を演算する。この位相差は、対象物が移動することによる移動量に相当するため、対象物の速度に対応することとなる。
そして、パルスドプラ計測装置は、位相差修正手段によって、周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで2πの整数倍の位相を加算又は減算することで位相差を修正する。この位相修正パターンは、どの周波数の位相差に対して2πの整数倍の位相を加算又は減算するかを示すものである。なお、この位相修正パターンは、最大周波数から順に2πの整数倍を加算又は減算するパターンが望ましい。これによって、位相の折り返し(エイリアシング)が発生する場合であっても、その折り返しを補正する方向に位相差が修正されることになる。
【0018】
また、パルスドプラ計測装置は、直線近似手段によって、周波数と当該周波数における位相差とを座標点とする座標系において、予め定めた複数の位相修正パターン毎に、周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する。なお、この直線近似は、最小二乗法を用いて演算することができる。また、この誤差が最も小さい場合、周波数と当該周波数における位相差との座標点が、ほぼ直線上にプロットされることになり、エイリアシングにより発生した位相差の位相の折り返しが修正されたことになる。
そして、パルスドプラ計測装置は、速度特定手段によって、直線近似手段で演算された誤差が最小となる直線の傾きを、対象物の移動速度に対応する位相差の傾きとして選択する。
【0019】
また、請求項2に記載のパルスドプラ計測装置は、請求項1に記載のパルスドプラ計測装置であって、前記位相差演算手段が、周波数毎の位相差を、予め定めた送受波回数毎に平均化することを特徴とする。
【0020】
かかる構成において、パルスドプラ計測装置は、位相差演算手段によって、反射波の変動を平均化することで、ノイズや、反射波の干渉による位相の変動を抑圧する。
【0021】
さらに、請求項3に記載のパルスドプラ計測装置は、請求項1又は請求項2に記載のパルスドプラ計測装置において、スペクトル解析手段で検出された位相に基づいて、クラッタ成分を除去するフィルタリング手段を備える構成とした。
【0022】
かかる構成において、パルスドプラ計測装置は、フィルタリング手段によって、スペクトル解析手段で検出された位相に基づいて、クラッタ成分を除去する。すなわち、フィルタリング手段によって、予め定めた基準よりも遅い位相回転のスペクトル成分(低周波スペクトル変動成分)を除去し、速い位相回転のスペクトル成分を抽出する。これによって、遅い動きによる信号成分(クラッタ成分)が除去される。
【0023】
また、請求項4に記載のパルスドプラ計測方法は、対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測方法であって、送受波ステップと、スペクトル解析ステップと、位相差演算ステップと、位相差修正ステップと、直線近似ステップと、速度特定ステップと、を含んでいることを特徴とする。
【0024】
かかる手順において、パルスドプラ計測方法は、送受波ステップで、超音波のパルスを、速度を計測する対象となる対象物に送波し、対象物からの反射された反射波を受波する。また、パルスドプラ計測方法は、スペクトル解析ステップで、反射波から、周波数毎の位相を検出する。
さらに、パルスドプラ計測方法は、位相差演算ステップで、スペクトル解析ステップで検出された周波数毎の位相について、周期毎に位相差を演算する。この位相差は、対象物が移動することによる移動量に相当するため、対象物の速度に対応することとなる。
【0025】
そして、パルスドプラ計測方法は、位相差修正ステップで、周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで2πの整数倍の位相を加算又は減算することで位相差を修正する。
また、パルスドプラ計測方法は、直線近似ステップで、周波数と当該周波数における位相差とを座標点とする座標系において、予め定めた複数の位相修正パターン毎に、周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する。
そして、パルスドプラ計測方法は、速度特定ステップで、直線近似ステップで演算された誤差が最小となる直線の傾きを、対象物の移動速度に対応する位相差の傾きとして選択する。
【0026】
さらに、請求項5に記載のパルスドプラ計測プログラムは、対象物に対して一定の周期で送波された超音波に対する反射波に基づいて、前記対象物の移動速度を計測するために、コンピュータを、スペクトル解析手段、位相差演算手段、位相差修正手段、直線近似手段、速度特定手段、として機能させることを特徴とする。
【0027】
かかる構成において、パルスドプラ計測プログラムは、スペクトル解析手段によって、反射波から、周波数毎の位相を検出する。
さらに、パルスドプラ計測プログラムは、位相差演算手段によって、スペクトル解析手段で検出された周波数毎の位相について、周期毎に位相差を演算する。
そして、パルスドプラ計測プログラムは、位相差修正手段によって、周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで2πの整数倍の位相を加算又は減算することで位相差を修正する。
【0028】
また、パルスドプラ計測プログラムは、直線近似手段によって、周波数と当該周波数における位相差とを座標点とする座標系において、予め定めた複数の位相修正パターン毎に、周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する。
そして、パルスドプラ計測プログラムは、速度特定手段によって、直線近似手段で演算された誤差が最小となる直線の傾きを、対象物の移動速度に対応する位相差の傾きとして選択する。
【発明の効果】
【0029】
本発明は、以下に示す優れた効果を奏するものである。
請求項1、請求項4又は請求項5に記載の発明によれば、反射波の位相差にエイリアシングが発生する場合であっても、その発生したエイリアシングを複数の周波数における位相差によって、修正することができる。これによって、本発明は、対象物の速度を計測する際の誤計測を防止することができる。
【0030】
請求項2に記載の発明によれば、反射波の変動を平均化することで、ノイズや、反射波の干渉による位相の変動を抑圧するため、例えば、血流のような多数散乱体を計測の対象物とする場合、反射波が重なり合い干渉する場合であっても、位相変動の影響を小さくすることができ、精度よく対象物の計測を行うことができる。
【0031】
請求項3に記載の発明によれば、フィルタリング手段が、複数の周波数においてクラッタ成分を除去するため、ある周波数において位相差が所定値未満となり速度が検出できない場合であっても、他の周波数において位相差が発生しているため、動いている対象物を動いていない物体であると認識することがない。このため、本発明は、速度検出ができない、いわゆるブラインド速度をなくすことができる。
【発明を実施するための最良の形態】
【0032】
以下、本発明の実施の形態について図面を参照して説明する。
[パルスドプラ計測方法(広帯域ドプラ法)の基本原理]
最初に、図3及び図4を参照して、本発明に係るパルスドプラ計測方法の基本原理について説明する。なお、ここでは、本発明に係るパルスドプラ計測方法を「広帯域ドプラ法」と呼び、一般的なパルスドプラ計測方法を「従来のパルスドプラ法」と呼ぶこととする。
【0033】
図3は、従来のパルスドプラ法と本発明に係る広帯域ドプラ法との違いを説明するためのグラフ図であって、(a)は従来のパルスドプラ法で送受波するパルスの振幅スペクトルと、伝播時間差における振幅スペクトルの位相差を示し、(b)は広帯域ドプラ法で送受波するパルスの振幅スペクトルと、伝播時間差における振幅スペクトルの位相差を示している。なお、図3の(a-1)及び(b-1)では、横軸を周波数、縦軸を振幅としている。また、図3の(a-2)及び(b-2)では、横軸を周波数、縦軸を位相差としている。
図4は、広帯域ドプラ法におけるエイリアシングが発生した場合の位相差の修正を説明するためのグラフ図である。なお、図4では、横軸を周波数、縦軸を位相差としている。
【0034】
通常のパルスドプラ法は、図3(a)の(a-1)に示すパルスPの中心周波数ωにおいて、(a-2)に示すように振幅スペクトルの位相差の傾きを求めることで、速度を得ている。この場合、位相差が±π以下となる低速で移動する対象物に対しては、正常に速度を得ることができるが、位相差が±πを超えるような高速で移動する対象物に対しては、位相差の位相の折り返し(エイリアシング)が発生し、位相差の傾きによる速度を正しく得ることができない。なお、従来のパルスドプラ法において、このエイリアシングの問題を解決する手法が、背景技術で説明した2周波法や2周期法である。
【0035】
一方、広帯域ドプラ法は、図3(b)の(b-1)に示すパルスPのスペクトル全域において、(b-2)に示すように複数の周波数(ω~ω)において位相差を求め、その複数の位相差により傾きを求めることで、速度を得ている。この場合であっても、位相差が±πを超えるような高速で移動する反射体に対しては、位相差の位相の折り返し(エイリアシング)が発生する。
【0036】
そこで、広帯域ドプラ法は、図4に示すように、位相の折り返しがあった周波数以上(ここではω以上)の位相差について、2πの位相を加えることで、折り返しの生じた位相差を修正する。なお、位相の折り返しがパルスPのスペクトル内で複数回発生する場合は、2πの整数倍の位相を加える。そして、広帯域ドプラ法は、位相差が修正された直線Lの傾きから速度を求める。なお、この直線の求め方については、後で詳細に説明を行うこととする。
すなわち、従来のパルスドプラ法は、中心周波数ωのみを使用するため、比較的パルス長の長い狭帯域なバースト波を用いるのに対し、広帯域ドプラ法は、逆に、複数の周波数を使用するため、広帯域のパルス波を用いることとしている。
なお、ここでは、位相差に2πを加算する例を示しているが、対象物B(図1参照)の動きが逆の場合は、位相差の傾きが負となるため、位相差から2πを減算ことにより位相差の修正を行うこととする。
【0037】
このように、広帯域ドプラ法は、複数の周波数に基づいて、エイリアシングが発生する位相差の傾きを修正するため、エイリアシングが発生するような速度で移動する移動体についても、その速度を計測することができる。
以下、この広帯域ドプラ法を実現するためのパルスドプラ計測装置の構成とその動作について詳細に説明を行う。
【0038】
[パルスドプラ計測装置の構成]
まず、図1を参照して、本発明に係るパルスドプラ計測装置の構成について説明する。図1は、パルスドプラ計測装置の構成を示すブロック図である。
図1に示したパルスドプラ計測装置は、対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、対象物の移動速度を計測するものである。なお、ここでは、パルスドプラ計測装置1は、生体内の血流の速度を計測するものとして機能させることとする。
図1に示すように、パルスドプラ計測装置1は、送受波装置2と、制御装置3とを備えている。
【0039】
送受波装置2は、対象物に対して超音波を送波し、対象物からの反射波を受波するものである。この送受波装置2は、送波手段21と受波手段22とを備えている。
【0040】
送波手段21は、超音波を一定の周期で対象物(ここでは、血流〔具体的には、赤血球〕)Bに送波するものである。この送波手段21は、従来のパルスドプラ計測装置で使用される送波パルスのパルス長よりも短い広帯域な超音波を送波する。
例えば、従来使用されていたパルス長が2.4μs(マイクロ秒)に対し、1.0μsのパルス長の送波パルスを使用する。このパルス長は、その帯域内に複数の周波数(図3(b)参照)が設定可能な長さであればよい。
【0041】
受波手段22は、送波手段21から送波された超音波に対して、対象物(血流)Bから反射される反射波を受波するものである。この受波手段22で受波された反射波は、デジタル信号に変換されて制御装置3に出力される。
【0042】
制御装置3は、送受波装置2で受波した反射波をデジタル信号として入力し、反射波の位相の変化に基づいて、対象物Bの移動速度を計測するものである。ここでは、制御装置3は、フーリエ変換手段31と、MTIフィルタ32と、位相差演算手段33と、位相差傾き演算手段34と、修正パターン記憶手段35と、速度特定手段36とを備えている。
【0043】
フーリエ変換手段(スペクトル解析手段)31は、送受波装置2の受波手段22で受波した反射波から、周波数毎の位相(スペクトル)を検出するものである。例えば、フーリエ変換手段31には、一般的なFFT(Fast Fourier Transform:高速フーリエ変換)を用いることができる。このときの周波数の標本化間隔は、検出すべき速度の上限によって決める。すなわち、速度の計測範囲において、周波数の標本化間隔内の位相回転が±π以内となるように予め標本化間隔を定めておく。
ここで、時刻tにおける(k-1)番目の反射波をfk-1(t)、フーリエ変換後の周波数ωにおけるスペクトル値をFk-1(ω)とする。また、時刻(t+T)におけるk番目の反射波をf(t)=fk-1(t+Δt)とする。ここでTは送波パルスの時間間隔で、Δtは対象物Bの移動により生じた伝搬時間の変化量である。すると、k番目の反射波f(t)をフーリエ変換した後の周波数ωにおけるスペクトル値F(ω)は、以下の(4)式で表される。
【0044】
【数4】
JP0004729765B2_000005t.gif

【0045】
この(4)式で演算されたスペクトル値F(ω)は、順次、MTIフィルタ32に出力される。
【0046】
MTIフィルタ(フィルタリング手段)32は、フーリエ変換手段31によって求められたスペクトルから、クラッタ成分を除去するものである。なお、ここでフィルタリングを行う基準となる閾値は、計測対象となる対象物の速度範囲により予め定めておく。
このMTIフィルタ32でフィルタリングされたフィルタリングされた信号は、位相差演算手段33に出力される。
【0047】
ここで、図5を参照して、MTIフィルタの原理について説明する。図5は、MTIフィルタの原理を説明するための説明図であって、(a)は従来のパルスドプラ法におけるMTIフィルタの原理、(b)は広帯域ドプラ法におけるMTIフィルタの原理を示している。
【0048】
通常、MTIフィルタは、1回目の受波と2回目の受波との差に基づいて、その差が所定値よりも小さいものを、血管壁、骨等の動きのない対象物からの反射波(クラッタ成分)として除去している。このように、MTIフィルタは、受波の差をとるため、対象物の動きが送波周期の整数倍に相当する距離だけ移動した場合、動きのない対象物とみなし、クラッタ成分として除去してしまう。すなわち、MTIフィルタにおける伝達特性は、図5(a)に示すような櫛形となり、計測することができない速度領域(ブラインド速度)が発生する。
そこで、従来のパルスドプラ法では、MTIフィルタにおいて、図5(a)に示すように、ブラインド速度が発生しない速度範囲(例えば、図中、速度範囲W)でフィルタリングを行い、速度の計測を行っている。
【0049】
一方、広帯域ドプラ法では、基本的なMTI処理については、従来と同様であるが、図5(b)に示すように、処理する周波数帯域が広く、複数の周波数(ω,ω,ω,…)において、フィルタリングを行うため、個々の周波数においてブラインド速度が発生しても、計測する周波数領域において、ブラインド速度が発生しなくなり、推定可能な速度領域を拡大させることができる。
図1に戻って、パルスドプラ計測装置の構成について、説明を続ける。
【0050】
位相差演算手段33は、MTIフィルタ32によってフィルタリングされた信号において、隣接する周期間における位相差を周波数毎に演算するものである。ここでは、位相差演算手段33は、位相差として、平均クロススペクトルの位相を用いることとする。平均クロススペクトルとは、2つのスペクトルの、ある周波数成分どうしを掛け合わせて平均化したものである。
すなわち、位相差演算手段33は、送受波をa回行う場合、以下の(5)式により平均クロススペクトルC(ω)を演算する。ただし、*は複素共役を示している。
【0051】
【数5】
JP0004729765B2_000006t.gif

【0052】
このように、予め定めた送受波回数毎に位相差を平均化することで、多数散乱体に対して、超音波のパルスを送波し、反射波が互いに重なり合い干渉する場合であっても、位相変動の影響を抑圧することができる。
この位相差演算手段33で演算された位相差(平均クロススペクトルの位相)は、位相差傾き演算手段34に出力される。
なお、位相差は、通常、以下の(6)式で表されるため、位相差演算手段33は、送受波をa回行う場合、以下の(7)式により平均の位相差Δθを求めることとしてもよい。
【0053】
【数6】
JP0004729765B2_000007t.gif

【0054】
【数7】
JP0004729765B2_000008t.gif

【0055】
位相差傾き演算手段34は、周波数毎の位相差に基づいて、周波数と当該周波数における位相差とを座標点とする座標系において、複数の位相修正パターン毎に、座標点を最小二乗法により直線に近似することで、その直線の傾き(位相差の傾き)と、その誤差とを演算するものである。
【0056】
ここで、図6を参照して、位相差傾き演算手段34において行う処理についてその概要を説明する。図6は、位相差における位相の折り返しが発生した場合における位相の修正の概念を示す概念図である。なお、図6では、横軸に周波数ω、縦軸に位相差(平均クロススペクトルの位相)を示し、使用する周波数をω~ω10としている。
【0057】
図6(a)は、周波数ω以上において、位相差における位相の折り返しが発生している状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似する。この場合、求めた直線Lの傾きは、位相の折り返し後(周波数ω以上)における位相の影響により、誤差が大きくなってしまう。
【0058】
図6(b)は、図6(a)の状態から、周波数ω及びω10の位相差に2πを加えた状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似すると、求めた直線Lの傾きの誤差は、図6(a)と比べ小さくなっているが、まだ、誤差を含んでいる。
【0059】
図6(c)は、図6(b)の状態から、周波数ω及びωの位相差に2πを加えた状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似すると、求めた直線Lの傾きの誤差が最も小さくなる。
【0060】
図6(d)は、図6(c)の状態から、周波数ω及びωの位相差に2πを加えた状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似すると、求めた直線Lの傾きの誤差は、逆に増加してしまう。
【0061】
このように、周波数ω~ω10の最大周波数から順に位相差に2πを加算することで、最小二乗法における誤差が最小となる直線を求めることができる。なお、ここでは、図の理解を容易にするため、2つの周波数毎に2πの位相を加算しているが、実際は、最大周波数から順に1つの周波数毎に2πを加算する。
なお、ここでは、最大周波数から順に位相差に2πを加算する例を示しているが、対象物B(図1参照)の動きの逆の場合は、2πを減算する。
【0062】
図1に戻って、パルスドプラ計測装置の構成について、説明を続ける。
図6で説明した位相の傾きを修正するため、ここでは、位相差傾き演算手段34は、位相差修正手段341と、直線近似手段342とを備えている。
【0063】
位相差修正手段341は、所定数の周波数の位相差(平均クロススペクトルの位相)に対して、当該周波数のうち最大周波数から順に2πの整数倍の位相を加算又は減算することで位相差を修正するものである。ここでは、位相差修正手段341は、後記する修正パターン記憶手段35により予め定めた複数の位相修正パターンで示される位相を加算又は減算するため、その位相修正パターンの数だけ備えられている(位相差修正手段341~341)。
ここでは、位相差修正手段341~341が、周波数毎に位相差に2πを加算する処理を並列処理する。これによって、例えば、図6(b)~(d)で示した位相差に2πを加算する処理が、個々の位相差修正手段341~341によって、並列に処理される。
この位相差修正手段341で修正された位相差(平均クロススペクトルの位相)の値は、当該位相差修正手段341に対応する直線近似手段342に出力される。
【0064】
直線近似手段342は、周波数と当該周波数における位相差とを座標点とする座標系において、複数の位相修正パターン毎に、座標点を最小二乗法により原点を通る直線に近似するものである。ここでは、位相差修正手段341は、位相差修正手段341に対応して、位相修正パターンの数だけ備えられている(直線近似手段342~342)。
そして、直線近似手段342は、近似した直線の傾き(位相差の傾き)と、最小二乗法における誤差とを速度特定手段36に出力する。
【0065】
ここで、数式により、直線近似手段342が、直線の傾きとその誤差とを求める手法について説明する。
まず、平均クロススペクトルの位相θを以下の(8)式、振幅qを以下の(9)式とする。
【0066】
【数8】
JP0004729765B2_000009t.gif

【0067】
【数9】
JP0004729765B2_000010t.gif

【0068】
ここで、Δωはクロススペクトルの標本化間隔、nは整数である。ここで、レンジゲート長をTとすると、Δω=2π/Tと表される。
また、直線の傾き(位相差の傾き)kは、以下の(10)式で表される。
【0069】
【数10】
JP0004729765B2_000011t.gif

【0070】
なお、この(10)式において、eは誤差項である。
このとき、平均クロススペクトルの振幅(振幅スペクトル)qで重み付けられた誤差の二乗和、すなわち、eQeを最小にする傾きkの推定量kは、以下の(11)式で求められる。
【0071】
【数11】
JP0004729765B2_000012t.gif

【0072】
また、推定誤差eQeは、以下の(12)式で求められる。
【0073】
【数12】
JP0004729765B2_000013t.gif

【0074】
このように求められた直線の傾き(推定量k)とその誤差(推定誤差eQe)は、速度特定手段36に出力される。
【0075】
修正パターン記憶手段35は、位相差修正手段341で加算する位相を位相修正パターンとして複数記憶するものであって、ハードディスク等の一般的な記憶装置である。
ここでは、修正パターン記憶手段35には、位相修正の対象となる周波数分の位相修正パターンが記憶されている。なお、この位相修正パターンには、位相を加算する対象となる周波数を示す位置と、その加算量を記述しておく。
【0076】
ここで、図7を参照(適宜図1参照)して、位相修正パターンの内容について具体的に説明する。図7は、位相修正パターンの一例を示すパターン図である。ここでは、位相を加算する対象となる周波数をω~ω15の15個とし、位相修正パターンの数をID~ID24の25個とした例を示している。ここで、ω15が最大周波数を示している。
この位相修正パターンにおいて、「0」は位相の加算を行わないことを示す。また、「1」、「2」、「3」、「4」はそれぞれ「1×2π」、「2×2π」、「3×2π」、「4×2π」の位相の加算を行うことを示している。
【0077】
また、ID~ID24は、位相差傾き演算手段34において、複数の位相差修正手段341のそれぞれが、固有のIDにより位相修正パターンを読み出し、位相差の修正を行う。これによって、位相差傾き演算手段34が並列処理により位相差の修正を行うことが可能になる。
【0078】
なお、ここでは、2πの整数倍の位相を加算する周波数を15個、位相修正パターンの数を25種類としているが、この数に限定されるものではない。この位相修正パターンを増やすことで、広範囲な速度を推定することが可能になる。
また、ここでは、2πの整数倍の位相を加算する位相修正パターンのみを示しているが、この位相修正パターンは、対象物の動きが逆の場合は、逆に2πを減算するパターンを示すこととする。また、対象物の動きの方向が予め定まっていない場合は、2πを加算するパターンと、減算するパターンとを合わせて設定しておくことが望ましい。
図1に戻って、パルスドプラ計測装置の構成について、説明を続ける。
【0079】
速度特定手段36は、直線近似手段342で演算された誤差が最小となる直線の傾きを、対象物Bの移動速度に対応する位相差の傾きとして選択するものである。
なお、この誤差が最小となる直線は、位相差の位相の折り返し(エイリアシング)を適正に修正したものであるため、対象物Bの動きにより位相の折り返しが発生する場合であっても、正しく対象物Bの速度が特定されることになる。
【0080】
以上説明した制御装置3は、CPU、RAM、ROM等を備えた一般的なコンピュータにより実現することができる。すなわち、制御装置3は、コンピュータを前記した各手段として機能させるパルスドプラ計測プログラムにより動作させることができる。
【0081】
以上、パルスドプラ計測装置1の構成について説明したが、本発明はこの構成に限定されるものではない。
ここでは、位相差修正手段341と直線近似手段342とを複数備え、並列に位相差の傾き(直線の傾き)を求めることとしたが、位相差修正手段341と直線近似手段342とをそれぞれ1つの構成とし、位相差の傾き(直線の傾き)をシリアルに順次演算することとしてもよい。
【0082】
また、ここでは、パルスドプラ計測装置1は、MTIフィルタ32を備えているが、測定地点から超音波を送波する方向に対象物以外の物体が存在しない環境において測定を行う場合、MTIフィルタ32を省略することとしてもよい。例えば、測定対象が測定地点に対して上空を飛行する飛行機の場合は、必ずしもMTIフィルタ32を備える必要はない。
なお、パルスドプラ計測装置1は、生体内の血流を計測する以外にも、前記した飛行機の速度測定や、雲の動きを測定する気象レーダ等、動きのある物体の速度計測に利用することができる。
【0083】
[パルスドプラ計測装置の動作]
次に、図2を参照(適宜図1参照)して、本発明に係るパルスドプラ計測装置の動作について説明する。図2は、パルスドプラ計測装置の動作を示すフローチャートである。
【0084】
(送受波ステップ)
まず、パルスドプラ計測装置1は、送受波装置2の送波手段21によって、超音波を一定の周期で対象物Bに送波し(ステップS1)、受波手段22によって、対象物Bから反射される反射波を受波する(ステップS2)。この受波した反射波は、デジタル信号に変換されて制御装置3に出力される。
【0085】
(スペクトル解析ステップ)
そして、パルスドプラ計測装置1は、制御装置3のフーリエ変換手段31によって、反射波に対して高速フーリエ変換(FFT)を行うことで、スペクトルを検出する(ステップS3)。
(フィルタリングステップ)
そして、パルスドプラ計測装置1は、MTIフィルタ32によって、ステップS3で検出されたスペクトルの位相回転の大きさによって、クラッタ成分を除去する(ステップS4)。
【0086】
(位相差演算ステップ)
そして、パルスドプラ計測装置1は、位相差演算手段33によって、隣接する周期間における位相差を周波数毎に演算する(ステップS5)。なお、ここでは、位相差として、前記(5)式に示すように、平均クロススペクトルの位相を演算する。
【0087】
その後、パルスドプラ計測装置1は、位相差傾き演算手段34によって、周波数と当該周波数における位相差とを座標点とする座標系において、複数の位相修正パターン毎に、座標点を最小二乗法により直線に近似することで、その直線の傾き(位相差の傾き)と、その誤差とを演算する。なお、本動作は、ここでは、位相修正パターン毎に、並列に動作する。
すなわち、パルスドプラ計測装置1は、位相差修正手段341によって、当該位相差修正手段341に対応する位相修正パターンを修正パターン記憶手段35から読み出す(ステップS6)。
【0088】
(位相差修正ステップ)
そして、位相差修正手段341は、周波数毎の位相差(平均クロススペクトルの位相)に、位相修正パターンに対応して、2πの整数倍の位相を加算する(ステップS7)。
【0089】
(直線近似ステップ)
その後、パルスドプラ計測装置1は、直線近似手段342によって、周波数と当該周波数における位相差とを座標点とする座標系において、各座標点を、最小二乗法により原点を通る直線に近似することで、その直線の傾き(前記(11)式参照)と、その誤差(前記(12)式参照)とを演算する(ステップS8)。
ステップS6~ステップS8、…、ステップS6~ステップS8の動作については、ステップS6~ステップS8の動作において、位相修正パターンが異なるだけであるため説明を省略する。
なお、ここでは、位相差傾き演算手段34の動作を、位相修正パターン毎に並列に動作させることとしたが、単一の位相差修正手段341と直線近似手段342とにより、ステップS6~ステップS8、ステップS6~ステップS8、…、ステップS6~ステップS8の動作を順次行うこととしてもよい。
【0090】
(速度特定ステップ)
そして、パルスドプラ計測装置1は、速度特定手段36によって、ステップS8、S8、…、S8で演算された誤差が最小となる直線の傾きを、対象物Bの移動速度に対応する位相差の傾きとして選択する(ステップS9)。
以上の動作によって、パルスドプラ計測装置1は、広帯域の周波数を用いて計測を行うため、対象物Bの動きにより位相差において位相の折り返しが発生する場合であっても、正しく対象物Bの速度を特定することができる。また、パルスドプラ計測装置1は、個々の周波数においてブラインド速度が発生しても、計測する周波数領域において、ブラインド速度が発生しないため、推定可能な速度領域を従来よりも拡大させることができる。
【0091】
[シミュレーション結果]
最後に、従来の手法(従来のパルスドプラ法、2周波法)と、本発明に係るパルスドプラ計測装置1による広帯域ドプラ法との速度推定の評価結果について説明する。
ここでは、一次元多数散乱体モデルを用いて、速度推定のシミュレーションを行った。シミュレーションの条件は、以下の表に示す通りである。
【0092】
【表1】
JP0004729765B2_000014t.gif

【0093】
なお、本シミュレーションにおいては、従来の手法と広帯域ドプラ法とでは、使用する送波パルスのパルス長のみを変え、従来の手法においては2.4μs、広帯域ドプラ法では、1.0μsとした。また、広帯域ドプラ法において、81種類の位相修正パターンを使用した。
【0094】
図8は、各手法におけるシミュレーション結果をグラフ化したものであり、(a)は従来のパルスドプラ法、(b)は従来の2周波法、(c)は本発明に係る広帯域ドプラ法による結果を示している。なお、ここでは2周波法は、特許文献2に記載されている手法を用いている。また、図8の各グラフでは、横軸に速度測定対象である散乱体の予め設定した速度、縦軸にシミュレーションにより推定した速度を示している。理想的には、予め設定した速度と推定した速度とが一致し、プロット位置(図中「+」)は、傾きが「1」の直線上に乗ることが望ましい。
【0095】
図8(a)に示すように、従来のパルスドプラ法では、ナイキスト限界が0.42m/sとなるため、これを超える速度ではエイリアシングが発生している。また、散乱体からの反射波の干渉により推定した速度にばらつきが見られる。
また、図8(b)に示すように、従来の2周波法では、エイリアシングは発生していないものの、干渉による位相変動の影響を受けるため、推定した速度に大きなばらつきが発生している。なお、このシミュレーション結果は、従来の2周期法にも同様のことが言える。
しかし、図8(c)に示すように、本発明に係る広帯域ドプラ法では、エイリアシングが発生せず、さらに図8(a)の従来のパルスドプラ法に比べ、推定速度のばらつきが小さく、±1.4m/sの範囲で精度よく速度推定がなされている。
このように、本発明に係る広帯域ドプラ法は、エイリアシングが発生せず、散乱体からの反射波の干渉の影響が少なく、精度よく速度推定を行うことができる。
【図面の簡単な説明】
【0096】
【図1】本発明に係るパルスドプラ計測装置の構成を示すブロック図である。
【図2】本発明に係るパルスドプラ計測装置の動作を示すフローチャートである。
【図3】従来のパルスドプラ法と本発明に係る広帯域ドプラ法との違いを説明するためのグラフ図である。
【図4】広帯域ドプラ法におけるエイリアシングが発生した場合の位相差の修正を説明するためのグラフ図である。
【図5】MTIフィルタの原理を説明するための説明図である。
【図6】位相の折り返しが発生した場合における位相の修正の概念を示す概念図である。
【図7】位相修正パターンの一例を示すパターン図である。
【図8】シミュレーション結果を示すグラフ図である。
【図9】従来のパルスドプラ法の概念を説明するための説明図である。
【符号の説明】
【0097】
1 パルスドプラ計測装置
2 送受波装置
21 送波手段
22 受波手段
3 制御装置
31 フーリエ変換手段(スペクトル解析手段)
32 MTIフィルタ(フィルタリング手段)
33 位相差演算手段
34 位相差傾き演算手段
341 位相差修正手段
342 直線近似手段
35 修正パターン記憶手段
36 速度特定手段
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8