JP2007215816A - パルスドプラ計測装置、その方法及びそのプログラム - Google Patents

Abstract

【課題】対象物の速度を超音波パルスの反射波の位相差により計測する際に、位相差における位相の折り返しによる誤計測を防止することを可能にしたパルスドプラ計測装置を提供する。
【解決手段】パルスドプラ計測装置１は、反射波から、周波数毎の位相を検出するフーリエ変換手段３１と、周波数毎の位相について位相差を演算する位相差演算手段３３と、その位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算する位相差修正手段３４１と、周波数と当該周波数における位相差とを座標点とする座標系において、位相修正パターン毎に周波数毎の座標点を直線に近似することで、その直線の傾きとその誤差とを演算する直線近似手段３４２と、その誤差が最小となる直線の傾きを、移動速度に対応する位相差の傾きとして選択する速度特定手段３６と、を備えることを特徴とする。
【選択図】図１

Description

本発明は、対象物に送波した超音波に対する反射波に基づいて、移動する対象物の速度を計測するパルスドプラ計測装置、その方法及びそのプログラムに関する。

従来、人体の診断を行う装置として、生体内の血流の速度を計測するパルスドプラ法を用いた診断装置が種々開発され、臨床の場で広く利用されている。これらの診断装置は、生体内の血流に対して超音波のパルスを一定周期で繰り返し送波し、計測対象である赤血球が動くことにより生じる反射波の位相変化を検出することで、血流の速度を計測している。
ここで、図９を参照して、従来のパルスドプラ法について、その概念を説明する。
パルスドプラ法は、図９（ａ）に示すように、送波器Ｓから、一定周期で超音波（送波パルス）を、動きを伴う対象物（例えば、血流）Ｂに送波し、受波器Ｒで、対象物Ｂからの反射波を受波する。このとき、１番目の送波と２番目の送波との間に対象物ＢがΔｘだけ移動したとすると、図９（ｂ）に示すように、伝播時間差Δｔだけ受波波形がシフトする。
この移動量Δｘと伝播時間差Δｔとの関係は、伝播媒質中の音速をｃ〔ｍ／ｓ〕とすると、以下の（１）式で表される。

ここで、送波パルスの繰り返し周期をＴ〔ｓ〕とし、対象物Ｂの速度をｖ〔ｍ／ｓ〕とすると、前記（１）式は、以下の（２）式となる。

この伝播時間差Δｔは、図９（ｃ）に示すように、シフトした受波のスペクトルにおいて位相のずれとなって現れる。すなわち、伝播時間差Δｔは、図９（ｄ）に示すように、シフトした受波のスペクトルの位相差（位相回転）に比例する。
ここで、送波パルスの中心周波数をω_ｃ、ω_ｃにおける位相差をΔθ_ｃとすると、図９（ｄ）の直線の傾きはΔθ_ｃ／ω_ｃとなる。これが伝播時間差Δｔに相当するため、前記（２）式より、以下の（３）式を得ることができる。

このように、対象物Ｂの速度ｖは、送波パルスの中心周波数ω_ｃと位相差Δθ_ｃとにより求めることができる。
そして、このパルスドプラ法を用いた診断装置は、血流の速度を表示装置上にカラー表示させることにより、心臓、肝臓等の診断、血行動態のリアルタイムの診断等、多くの診断に利用されている。
しかし、パルスドプラ法を用いた診断装置では、前記（３）式に示すように、血流計測を反射波の位相変化により行っているため、図９（ｅ）に示すように、血流の速さによっては、位相の折り返し（いわゆるエイリアシング）が発生し、正しい計測結果が得られないという問題がある。

この問題を解決する手法として、種々の提案がなされている。例えば、２つの異なる周波数をもつパルスを予め定めた周期で対象物に送波し、それぞれ送波したパルスと反射波との位相差の差をとることで、位相差が±πを超えエイリアシングが発生する場合でも、２つの位相差の差が±π以内でありエイリアシングが発生しないことを利用して、対象物の速度を計測する手法（以下、２周波法という）が提案されている（特許文献１、特許文献２参照）。
また、別の手法として、パルスを異なる周期で対象物に送波し、それぞれ送波したパルスと反射波との位相差をとり、周期の差により発生する位相差の差を利用することで、対象物の速度を計測する手法（以下、２周期法という）が提案されている（特許文献３参照）。

また、これらの従来の手法では、送受波を複数回行い位相の平均化を行うことで、計測速度の精度を高めている。さらに、従来の手法では、送受波するパルスの中心周波数により位相の検出を行っているが、パルスの帯域を狭くする方が、Ｓ／Ｎの点で有利であると考え、比較的パルス長の長いバースト波を用いている。
また、従来の手法では、実際の血流計測において、血管壁、骨等の動きのない反射体からの反射波を拾ってしまうクラッタ成分を除去するため、反射波の位相の差に基づいて、動きのない反射体からの反射波の信号を除去するＭＴＩ（Moving Target Indication）フィルタを使用している。
特許第２９５３０８３号公報（段落０００９） 米国特許第４５３４３５７号明細書 特開平４−１９７２４９号公報（第３−５頁、第３−５図）

前記した従来の手法（２周波法、２周期法）では、Ｓ／Ｎ改善の観点から、受波した反射波の位相の平均化を行ったり、比較的パルス長の長いバースト波を使用したりすることで、計測の精度を高めようとしていた。
しかし、実際の生体内の血流計測において、精度を高めることができない原因は、各々の赤血球からの反射波の干渉による位相の変動であることを、本発明者は解明した。

概略すると、生体内の血流を計測する場合、超音波の反射体となるのは赤血球である。この赤血球は、超音波の一波長内に極めて多数かつランダムに存在し、多数散乱体として振舞う。このような多数散乱体に対して、超音波のパルスを送波すると、各々の赤血球からの反射波は互いに重なり合い干渉することになる。その結果、血流が一定の速度で流れていたとしても、反射波の位相は不規則に変動することになる。
これに対し、従来の手法（２周波法、２周期法）は、干渉による位相変動の影響を考慮していないため、計測の精度を高めることができなかった。すなわち、従来の手法においては、位相差の差を求めるため、位相の変動による影響を強く受けることになり、エイリアシングの問題は解決できても、計測の精度を逆に劣化させることとなっていた。

さらに、実際の血流計測において、従来の手法は、ＭＴＩフィルタによりクラッタ成分を除去している。しかし、従来の手法は、送受波するパルスの中心周波数により位相の検出を行っているため、送波周期の整数倍に相当する距離だけ移動した反射体については、位相差がゼロになってしまい、実際には対象物が動いているにも拘らず、クラッタ成分として除去されてしまう。このように、従来の手法では、求められない速度（以下、ブラインド速度という）が存在するという問題があった。

本発明は、以上のような課題を解決するためになされたものであり、対象物の速度を超音波パルスの反射波の位相差により計測する際に、位相差における位相の折り返し（エイリアシング）による誤計測を防止するとともに、対象物の速度の計測精度を高め、さらに、ブラインド速度をなくすことを可能にしたパルスドプラ計測装置、その方法及びそのプログラムを提供することを目的とする。

本発明は、前記目的を達成するために創案されたものであり、まず、請求項１に記載のパルスドプラ計測装置は、対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測装置であって、送波手段と、受波手段と、スペクトル解析手段と、位相差演算手段と、位相差修正手段と、直線近似手段と、速度特定手段と、を備える構成とした。

かかる構成において、パルスドプラ計測装置は、送波手段によって、超音波のパルスを、速度を計測する対象となる対象物に送波し、受波手段によって、対象物から反射された反射波を受波する。このとき、対象物が移動していると反射波の位相がずれることになる。
また、パルスドプラ計測装置は、スペクトル解析手段によって、反射波から、周波数毎の位相を検出する。すなわち、スペクトル解析手段は、送波したパルスにおける１つの周波数だけでなく、予め定めた複数の周波数における位相を検出する。

さらに、パルスドプラ計測装置は、位相差演算手段によって、スペクトル解析手段で検出された周波数毎の位相について、周期毎に位相差を演算する。この位相差は、対象物が移動することによる移動量に相当するため、対象物の速度に対応することとなる。
そして、パルスドプラ計測装置は、位相差修正手段によって、周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算することで位相差を修正する。この位相修正パターンは、どの周波数の位相差に対して２πの整数倍の位相を加算又は減算するかを示すものである。なお、この位相修正パターンは、最大周波数から順に２πの整数倍を加算又は減算するパターンが望ましい。これによって、位相の折り返し（エイリアシング）が発生する場合であっても、その折り返しを補正する方向に位相差が修正されることになる。

また、パルスドプラ計測装置は、直線近似手段によって、周波数と当該周波数における位相差とを座標点とする座標系において、予め定めた複数の位相修正パターン毎に、周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する。なお、この直線近似は、最小二乗法を用いて演算することができる。また、この誤差が最も小さい場合、周波数と当該周波数における位相差との座標点が、ほぼ直線上にプロットされることになり、エイリアシングにより発生した位相差の位相の折り返しが修正されたことになる。
そして、パルスドプラ計測装置は、速度特定手段によって、直線近似手段で演算された誤差が最小となる直線の傾きを、対象物の移動速度に対応する位相差の傾きとして選択する。

また、請求項２に記載のパルスドプラ計測装置は、請求項１に記載のパルスドプラ計測装置であって、前記位相差演算手段が、周波数毎の位相差を、予め定めた送受波回数毎に平均化することを特徴とする。

かかる構成において、パルスドプラ計測装置は、位相差演算手段によって、反射波の変動を平均化することで、ノイズや、反射波の干渉による位相の変動を抑圧する。

さらに、請求項３に記載のパルスドプラ計測装置は、請求項１又は請求項２に記載のパルスドプラ計測装置において、スペクトル解析手段で検出された位相に基づいて、クラッタ成分を除去するフィルタリング手段を備える構成とした。

かかる構成において、パルスドプラ計測装置は、フィルタリング手段によって、スペクトル解析手段で検出された位相に基づいて、クラッタ成分を除去する。すなわち、フィルタリング手段によって、予め定めた基準よりも遅い位相回転のスペクトル成分（低周波スペクトル変動成分）を除去し、速い位相回転のスペクトル成分を抽出する。これによって、遅い動きによる信号成分（クラッタ成分）が除去される。

また、請求項４に記載のパルスドプラ計測方法は、対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測方法であって、送受波ステップと、スペクトル解析ステップと、位相差演算ステップと、位相差修正ステップと、直線近似ステップと、速度特定ステップと、を含んでいることを特徴とする。

かかる手順において、パルスドプラ計測方法は、送受波ステップで、超音波のパルスを、速度を計測する対象となる対象物に送波し、対象物からの反射された反射波を受波する。また、パルスドプラ計測方法は、スペクトル解析ステップで、反射波から、周波数毎の位相を検出する。
さらに、パルスドプラ計測方法は、位相差演算ステップで、スペクトル解析ステップで検出された周波数毎の位相について、周期毎に位相差を演算する。この位相差は、対象物が移動することによる移動量に相当するため、対象物の速度に対応することとなる。

そして、パルスドプラ計測方法は、位相差修正ステップで、周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算することで位相差を修正する。
また、パルスドプラ計測方法は、直線近似ステップで、周波数と当該周波数における位相差とを座標点とする座標系において、予め定めた複数の位相修正パターン毎に、周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する。
そして、パルスドプラ計測方法は、速度特定ステップで、直線近似ステップで演算された誤差が最小となる直線の傾きを、対象物の移動速度に対応する位相差の傾きとして選択する。

さらに、請求項５に記載のパルスドプラ計測プログラムは、対象物に対して一定の周期で送波された超音波に対する反射波に基づいて、前記対象物の移動速度を計測するために、コンピュータを、スペクトル解析手段、位相差演算手段、位相差修正手段、直線近似手段、速度特定手段、として機能させることを特徴とする。

かかる構成において、パルスドプラ計測プログラムは、スペクトル解析手段によって、反射波から、周波数毎の位相を検出する。
さらに、パルスドプラ計測プログラムは、位相差演算手段によって、スペクトル解析手段で検出された周波数毎の位相について、周期毎に位相差を演算する。
そして、パルスドプラ計測プログラムは、位相差修正手段によって、周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算することで位相差を修正する。

また、パルスドプラ計測プログラムは、直線近似手段によって、周波数と当該周波数における位相差とを座標点とする座標系において、予め定めた複数の位相修正パターン毎に、周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する。
そして、パルスドプラ計測プログラムは、速度特定手段によって、直線近似手段で演算された誤差が最小となる直線の傾きを、対象物の移動速度に対応する位相差の傾きとして選択する。

本発明は、以下に示す優れた効果を奏するものである。
請求項１、請求項４又は請求項５に記載の発明によれば、反射波の位相差にエイリアシングが発生する場合であっても、その発生したエイリアシングを複数の周波数における位相差によって、修正することができる。これによって、本発明は、対象物の速度を計測する際の誤計測を防止することができる。

請求項２に記載の発明によれば、反射波の変動を平均化することで、ノイズや、反射波の干渉による位相の変動を抑圧するため、例えば、血流のような多数散乱体を計測の対象物とする場合、反射波が重なり合い干渉する場合であっても、位相変動の影響を小さくすることができ、精度よく対象物の計測を行うことができる。

請求項３に記載の発明によれば、フィルタリング手段が、複数の周波数においてクラッタ成分を除去するため、ある周波数において位相差が所定値未満となり速度が検出できない場合であっても、他の周波数において位相差が発生しているため、動いている対象物を動いていない物体であると認識することがない。このため、本発明は、速度検出ができない、いわゆるブラインド速度をなくすことができる。

以下、本発明の実施の形態について図面を参照して説明する。
［パルスドプラ計測方法（広帯域ドプラ法）の基本原理］
最初に、図３及び図４を参照して、本発明に係るパルスドプラ計測方法の基本原理について説明する。なお、ここでは、本発明に係るパルスドプラ計測方法を「広帯域ドプラ法」と呼び、一般的なパルスドプラ計測方法を「従来のパルスドプラ法」と呼ぶこととする。

図３は、従来のパルスドプラ法と本発明に係る広帯域ドプラ法との違いを説明するためのグラフ図であって、（ａ）は従来のパルスドプラ法で送受波するパルスの振幅スペクトルと、伝播時間差における振幅スペクトルの位相差を示し、（ｂ）は広帯域ドプラ法で送受波するパルスの振幅スペクトルと、伝播時間差における振幅スペクトルの位相差を示している。なお、図３の（ａ−１）及び（ｂ−１）では、横軸を周波数、縦軸を振幅としている。また、図３の（ａ−２）及び（ｂ−２）では、横軸を周波数、縦軸を位相差としている。
図４は、広帯域ドプラ法におけるエイリアシングが発生した場合の位相差の修正を説明するためのグラフ図である。なお、図４では、横軸を周波数、縦軸を位相差としている。

通常のパルスドプラ法は、図３（ａ）の（ａ−１）に示すパルスＰ_Ｎの中心周波数ω_ｃにおいて、（ａ−２）に示すように振幅スペクトルの位相差の傾きを求めることで、速度を得ている。この場合、位相差が±π以下となる低速で移動する対象物に対しては、正常に速度を得ることができるが、位相差が±πを超えるような高速で移動する対象物に対しては、位相差の位相の折り返し（エイリアシング）が発生し、位相差の傾きによる速度を正しく得ることができない。なお、従来のパルスドプラ法において、このエイリアシングの問題を解決する手法が、背景技術で説明した２周波法や２周期法である。

一方、広帯域ドプラ法は、図３（ｂ）の（ｂ−１）に示すパルスＰ_Ｗのスペクトル全域において、（ｂ−２）に示すように複数の周波数（ω_１〜ω_ｎ）において位相差を求め、その複数の位相差により傾きを求めることで、速度を得ている。この場合であっても、位相差が±πを超えるような高速で移動する反射体に対しては、位相差の位相の折り返し（エイリアシング）が発生する。

そこで、広帯域ドプラ法は、図４に示すように、位相の折り返しがあった周波数以上（ここではω_７以上）の位相差について、２πの位相を加えることで、折り返しの生じた位相差を修正する。なお、位相の折り返しがパルスＰ_Ｗのスペクトル内で複数回発生する場合は、２πの整数倍の位相を加える。そして、広帯域ドプラ法は、位相差が修正された直線Ｌの傾きから速度を求める。なお、この直線の求め方については、後で詳細に説明を行うこととする。
すなわち、従来のパルスドプラ法は、中心周波数ω_ｃのみを使用するため、比較的パルス長の長い狭帯域なバースト波を用いるのに対し、広帯域ドプラ法は、逆に、複数の周波数を使用するため、広帯域のパルス波を用いることとしている。
なお、ここでは、位相差に２πを加算する例を示しているが、対象物Ｂ（図１参照）の動きが逆の場合は、位相差の傾きが負となるため、位相差から２πを減算ことにより位相差の修正を行うこととする。

このように、広帯域ドプラ法は、複数の周波数に基づいて、エイリアシングが発生する位相差の傾きを修正するため、エイリアシングが発生するような速度で移動する移動体についても、その速度を計測することができる。
以下、この広帯域ドプラ法を実現するためのパルスドプラ計測装置の構成とその動作について詳細に説明を行う。

［パルスドプラ計測装置の構成］
まず、図１を参照して、本発明に係るパルスドプラ計測装置の構成について説明する。図１は、パルスドプラ計測装置の構成を示すブロック図である。
図１に示したパルスドプラ計測装置は、対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、対象物の移動速度を計測するものである。なお、ここでは、パルスドプラ計測装置１は、生体内の血流の速度を計測するものとして機能させることとする。
図１に示すように、パルスドプラ計測装置１は、送受波装置２と、制御装置３とを備えている。

送受波装置２は、対象物に対して超音波を送波し、対象物からの反射波を受波するものである。この送受波装置２は、送波手段２１と受波手段２２とを備えている。

送波手段２１は、超音波を一定の周期で対象物（ここでは、血流〔具体的には、赤血球〕）Ｂに送波するものである。この送波手段２１は、従来のパルスドプラ計測装置で使用される送波パルスのパルス長よりも短い広帯域な超音波を送波する。
例えば、従来使用されていたパルス長が２．４μｓ（マイクロ秒）に対し、１．０μｓのパルス長の送波パルスを使用する。このパルス長は、その帯域内に複数の周波数（図３（ｂ）参照）が設定可能な長さであればよい。

受波手段２２は、送波手段２１から送波された超音波に対して、対象物（血流）Ｂから反射される反射波を受波するものである。この受波手段２２で受波された反射波は、デジタル信号に変換されて制御装置３に出力される。

制御装置３は、送受波装置２で受波した反射波をデジタル信号として入力し、反射波の位相の変化に基づいて、対象物Ｂの移動速度を計測するものである。ここでは、制御装置３は、フーリエ変換手段３１と、ＭＴＩフィルタ３２と、位相差演算手段３３と、位相差傾き演算手段３４と、修正パターン記憶手段３５と、速度特定手段３６とを備えている。

フーリエ変換手段（スペクトル解析手段）３１は、送受波装置２の受波手段２２で受波した反射波から、周波数毎の位相（スペクトル）を検出するものである。例えば、フーリエ変換手段３１には、一般的なＦＦＴ（Fast Fourier Transform：高速フーリエ変換）を用いることができる。このときの周波数の標本化間隔は、検出すべき速度の上限によって決める。すなわち、速度の計測範囲において、周波数の標本化間隔内の位相回転が±π以内となるように予め標本化間隔を定めておく。
ここで、時刻ｔにおける（ｋ−１）番目の反射波をｆ_ｋ−１（ｔ）、フーリエ変換後の周波数ωにおけるスペクトル値をＦ_ｋ−１（ω）とする。また、時刻（ｔ＋Ｔ）におけるｋ番目の反射波をｆ_ｋ（ｔ）＝ｆ_ｋ−１（ｔ＋Δｔ）とする。ここでＴは送波パルスの時間間隔で、Δｔは対象物Ｂの移動により生じた伝搬時間の変化量である。すると、ｋ番目の反射波ｆ_ｋ（ｔ）をフーリエ変換した後の周波数ωにおけるスペクトル値Ｆ_ｋ（ω）は、以下の（４）式で表される。

この（４）式で演算されたスペクトル値Ｆ_ｋ（ω）は、順次、ＭＴＩフィルタ３２に出力される。

ＭＴＩフィルタ（フィルタリング手段）３２は、フーリエ変換手段３１によって求められたスペクトルから、クラッタ成分を除去するものである。なお、ここでフィルタリングを行う基準となる閾値は、計測対象となる対象物の速度範囲により予め定めておく。
このＭＴＩフィルタ３２でフィルタリングされたフィルタリングされた信号は、位相差演算手段３３に出力される。

ここで、図５を参照して、ＭＴＩフィルタの原理について説明する。図５は、ＭＴＩフィルタの原理を説明するための説明図であって、（ａ）は従来のパルスドプラ法におけるＭＴＩフィルタの原理、（ｂ）は広帯域ドプラ法におけるＭＴＩフィルタの原理を示している。

通常、ＭＴＩフィルタは、１回目の受波と２回目の受波との差に基づいて、その差が所定値よりも小さいものを、血管壁、骨等の動きのない対象物からの反射波（クラッタ成分）として除去している。このように、ＭＴＩフィルタは、受波の差をとるため、対象物の動きが送波周期の整数倍に相当する距離だけ移動した場合、動きのない対象物とみなし、クラッタ成分として除去してしまう。すなわち、ＭＴＩフィルタにおける伝達特性は、図５（ａ）に示すような櫛形となり、計測することができない速度領域（ブラインド速度）が発生する。
そこで、従来のパルスドプラ法では、ＭＴＩフィルタにおいて、図５（ａ）に示すように、ブラインド速度が発生しない速度範囲（例えば、図中、速度範囲Ｗ）でフィルタリングを行い、速度の計測を行っている。

一方、広帯域ドプラ法では、基本的なＭＴＩ処理については、従来と同様であるが、図５（ｂ）に示すように、処理する周波数帯域が広く、複数の周波数（ω_１，ω_２，ω_３，…）において、フィルタリングを行うため、個々の周波数においてブラインド速度が発生しても、計測する周波数領域において、ブラインド速度が発生しなくなり、推定可能な速度領域を拡大させることができる。
図１に戻って、パルスドプラ計測装置の構成について、説明を続ける。

位相差演算手段３３は、ＭＴＩフィルタ３２によってフィルタリングされた信号において、隣接する周期間における位相差を周波数毎に演算するものである。ここでは、位相差演算手段３３は、位相差として、平均クロススペクトルの位相を用いることとする。平均クロススペクトルとは、２つのスペクトルの、ある周波数成分どうしを掛け合わせて平均化したものである。
すなわち、位相差演算手段３３は、送受波をａ回行う場合、以下の（５）式により平均クロススペクトルＣ（ω）を演算する。ただし、＊は複素共役を示している。

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

位相差傾き演算手段３４は、周波数毎の位相差に基づいて、周波数と当該周波数における位相差とを座標点とする座標系において、複数の位相修正パターン毎に、座標点を最小二乗法により直線に近似することで、その直線の傾き（位相差の傾き）と、その誤差とを演算するものである。

ここで、図６を参照して、位相差傾き演算手段３４において行う処理についてその概要を説明する。図６は、位相差における位相の折り返しが発生した場合における位相の修正の概念を示す概念図である。なお、図６では、横軸に周波数ω、縦軸に位相差（平均クロススペクトルの位相）を示し、使用する周波数をω_１〜ω_１０としている。

図６（ａ）は、周波数ω_７以上において、位相差における位相の折り返しが発生している状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似する。この場合、求めた直線Ｌ_ａの傾きは、位相の折り返し後（周波数ω_７以上）における位相の影響により、誤差が大きくなってしまう。

図６（ｂ）は、図６（ａ）の状態から、周波数ω_９及びω_１０の位相差に２πを加えた状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似すると、求めた直線Ｌ_ｂの傾きの誤差は、図６（ａ）と比べ小さくなっているが、まだ、誤差を含んでいる。

図６（ｃ）は、図６（ｂ）の状態から、周波数ω_７及びω_８の位相差に２πを加えた状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似すると、求めた直線Ｌ_ｃの傾きの誤差が最も小さくなる。

図６（ｄ）は、図６（ｃ）の状態から、周波数ω_５及びω_６の位相差に２πを加えた状態を示している。ここで、周波数毎の位相差の座標を、最小二乗法により原点を通る直線に近似すると、求めた直線Ｌ_ｄの傾きの誤差は、逆に増加してしまう。

このように、周波数ω_１〜ω_１０の最大周波数から順に位相差に２πを加算することで、最小二乗法における誤差が最小となる直線を求めることができる。なお、ここでは、図の理解を容易にするため、２つの周波数毎に２πの位相を加算しているが、実際は、最大周波数から順に１つの周波数毎に２πを加算する。
なお、ここでは、最大周波数から順に位相差に２πを加算する例を示しているが、対象物Ｂ（図１参照）の動きの逆の場合は、２πを減算する。

図１に戻って、パルスドプラ計測装置の構成について、説明を続ける。
図６で説明した位相の傾きを修正するため、ここでは、位相差傾き演算手段３４は、位相差修正手段３４１と、直線近似手段３４２とを備えている。

位相差修正手段３４１は、所定数の周波数の位相差（平均クロススペクトルの位相）に対して、当該周波数のうち最大周波数から順に２πの整数倍の位相を加算又は減算することで位相差を修正するものである。ここでは、位相差修正手段３４１は、後記する修正パターン記憶手段３５により予め定めた複数の位相修正パターンで示される位相を加算又は減算するため、その位相修正パターンの数だけ備えられている（位相差修正手段３４１_１〜３４１_ｎ）。
ここでは、位相差修正手段３４１_１〜３４１_ｎが、周波数毎に位相差に２πを加算する処理を並列処理する。これによって、例えば、図６（ｂ）〜（ｄ）で示した位相差に２πを加算する処理が、個々の位相差修正手段３４１_１〜３４１_ｎによって、並列に処理される。
この位相差修正手段３４１で修正された位相差（平均クロススペクトルの位相）の値は、当該位相差修正手段３４１に対応する直線近似手段３４２に出力される。

直線近似手段３４２は、周波数と当該周波数における位相差とを座標点とする座標系において、複数の位相修正パターン毎に、座標点を最小二乗法により原点を通る直線に近似するものである。ここでは、位相差修正手段３４１は、位相差修正手段３４１に対応して、位相修正パターンの数だけ備えられている（直線近似手段３４２_１〜３４２_ｎ）。
そして、直線近似手段３４２は、近似した直線の傾き（位相差の傾き）と、最小二乗法における誤差とを速度特定手段３６に出力する。

ここで、数式により、直線近似手段３４２が、直線の傾きとその誤差とを求める手法について説明する。
まず、平均クロススペクトルの位相θ_ｎを以下の（８）式、振幅ｑ_ｎを以下の（９）式とする。

ここで、Δωはクロススペクトルの標本化間隔、ｎは整数である。ここで、レンジゲート長をＴ_Ｒとすると、Δω＝２π／Ｔ_Ｒと表される。
また、直線の傾き（位相差の傾き）ｋは、以下の（１０）式で表される。

なお、この（１０）式において、ｅは誤差項である。
このとき、平均クロススペクトルの振幅（振幅スペクトル）ｑ_ｎで重み付けられた誤差の二乗和、すなわち、ｅ^ＴＱｅを最小にする傾きｋの推定量ｋ_ｅは、以下の（１１）式で求められる。

また、推定誤差ｅ^ＴＱｅは、以下の（１２）式で求められる。

このように求められた直線の傾き（推定量ｋ_ｅ）とその誤差（推定誤差ｅ^ＴＱｅ）は、速度特定手段３６に出力される。

修正パターン記憶手段３５は、位相差修正手段３４１で加算する位相を位相修正パターンとして複数記憶するものであって、ハードディスク等の一般的な記憶装置である。
ここでは、修正パターン記憶手段３５には、位相修正の対象となる周波数分の位相修正パターンが記憶されている。なお、この位相修正パターンには、位相を加算する対象となる周波数を示す位置と、その加算量を記述しておく。

ここで、図７を参照（適宜図１参照）して、位相修正パターンの内容について具体的に説明する。図７は、位相修正パターンの一例を示すパターン図である。ここでは、位相を加算する対象となる周波数をω_１〜ω_１５の１５個とし、位相修正パターンの数をＩＤ_０〜ＩＤ_２４の２５個とした例を示している。ここで、ω_１５が最大周波数を示している。
この位相修正パターンにおいて、「０」は位相の加算を行わないことを示す。また、「１」、「２」、「３」、「４」はそれぞれ「１×２π」、「２×２π」、「３×２π」、「４×２π」の位相の加算を行うことを示している。

また、ＩＤ_０〜ＩＤ_２４は、位相差傾き演算手段３４において、複数の位相差修正手段３４１のそれぞれが、固有のＩＤにより位相修正パターンを読み出し、位相差の修正を行う。これによって、位相差傾き演算手段３４が並列処理により位相差の修正を行うことが可能になる。

なお、ここでは、２πの整数倍の位相を加算する周波数を１５個、位相修正パターンの数を２５種類としているが、この数に限定されるものではない。この位相修正パターンを増やすことで、広範囲な速度を推定することが可能になる。
また、ここでは、２πの整数倍の位相を加算する位相修正パターンのみを示しているが、この位相修正パターンは、対象物の動きが逆の場合は、逆に２πを減算するパターンを示すこととする。また、対象物の動きの方向が予め定まっていない場合は、２πを加算するパターンと、減算するパターンとを合わせて設定しておくことが望ましい。
図１に戻って、パルスドプラ計測装置の構成について、説明を続ける。

速度特定手段３６は、直線近似手段３４２で演算された誤差が最小となる直線の傾きを、対象物Ｂの移動速度に対応する位相差の傾きとして選択するものである。
なお、この誤差が最小となる直線は、位相差の位相の折り返し（エイリアシング）を適正に修正したものであるため、対象物Ｂの動きにより位相の折り返しが発生する場合であっても、正しく対象物Ｂの速度が特定されることになる。

以上説明した制御装置３は、ＣＰＵ、ＲＡＭ、ＲＯＭ等を備えた一般的なコンピュータにより実現することができる。すなわち、制御装置３は、コンピュータを前記した各手段として機能させるパルスドプラ計測プログラムにより動作させることができる。

以上、パルスドプラ計測装置１の構成について説明したが、本発明はこの構成に限定されるものではない。
ここでは、位相差修正手段３４１と直線近似手段３４２とを複数備え、並列に位相差の傾き（直線の傾き）を求めることとしたが、位相差修正手段３４１と直線近似手段３４２とをそれぞれ１つの構成とし、位相差の傾き（直線の傾き）をシリアルに順次演算することとしてもよい。

また、ここでは、パルスドプラ計測装置１は、ＭＴＩフィルタ３２を備えているが、測定地点から超音波を送波する方向に対象物以外の物体が存在しない環境において測定を行う場合、ＭＴＩフィルタ３２を省略することとしてもよい。例えば、測定対象が測定地点に対して上空を飛行する飛行機の場合は、必ずしもＭＴＩフィルタ３２を備える必要はない。
なお、パルスドプラ計測装置１は、生体内の血流を計測する以外にも、前記した飛行機の速度測定や、雲の動きを測定する気象レーダ等、動きのある物体の速度計測に利用することができる。

［パルスドプラ計測装置の動作］
次に、図２を参照（適宜図１参照）して、本発明に係るパルスドプラ計測装置の動作について説明する。図２は、パルスドプラ計測装置の動作を示すフローチャートである。

（送受波ステップ）
まず、パルスドプラ計測装置１は、送受波装置２の送波手段２１によって、超音波を一定の周期で対象物Ｂに送波し（ステップＳ１）、受波手段２２によって、対象物Ｂから反射される反射波を受波する（ステップＳ２）。この受波した反射波は、デジタル信号に変換されて制御装置３に出力される。

（スペクトル解析ステップ）
そして、パルスドプラ計測装置１は、制御装置３のフーリエ変換手段３１によって、反射波に対して高速フーリエ変換（ＦＦＴ）を行うことで、スペクトルを検出する（ステップＳ３）。
（フィルタリングステップ）
そして、パルスドプラ計測装置１は、ＭＴＩフィルタ３２によって、ステップＳ３で検出されたスペクトルの位相回転の大きさによって、クラッタ成分を除去する（ステップＳ４）。

（位相差演算ステップ）
そして、パルスドプラ計測装置１は、位相差演算手段３３によって、隣接する周期間における位相差を周波数毎に演算する（ステップＳ５）。なお、ここでは、位相差として、前記（５）式に示すように、平均クロススペクトルの位相を演算する。

その後、パルスドプラ計測装置１は、位相差傾き演算手段３４によって、周波数と当該周波数における位相差とを座標点とする座標系において、複数の位相修正パターン毎に、座標点を最小二乗法により直線に近似することで、その直線の傾き（位相差の傾き）と、その誤差とを演算する。なお、本動作は、ここでは、位相修正パターン毎に、並列に動作する。
すなわち、パルスドプラ計測装置１は、位相差修正手段３４１_１によって、当該位相差修正手段３４１_１に対応する位相修正パターンを修正パターン記憶手段３５から読み出す（ステップＳ６_１）。

（位相差修正ステップ）
そして、位相差修正手段３４１_１は、周波数毎の位相差（平均クロススペクトルの位相）に、位相修正パターンに対応して、２πの整数倍の位相を加算する（ステップＳ７_１）。

（直線近似ステップ）
その後、パルスドプラ計測装置１は、直線近似手段３４２_１によって、周波数と当該周波数における位相差とを座標点とする座標系において、各座標点を、最小二乗法により原点を通る直線に近似することで、その直線の傾き（前記（１１）式参照）と、その誤差（前記（１２）式参照）とを演算する（ステップＳ８_１）。
ステップＳ６_２〜ステップＳ８_２、…、ステップＳ６_ｎ〜ステップＳ８_ｎの動作については、ステップＳ６_１〜ステップＳ８_１の動作において、位相修正パターンが異なるだけであるため説明を省略する。
なお、ここでは、位相差傾き演算手段３４の動作を、位相修正パターン毎に並列に動作させることとしたが、単一の位相差修正手段３４１と直線近似手段３４２とにより、ステップＳ６_１〜ステップＳ８_１、ステップＳ６_２〜ステップＳ８_２、…、ステップＳ６_ｎ〜ステップＳ８_ｎの動作を順次行うこととしてもよい。

（速度特定ステップ）
そして、パルスドプラ計測装置１は、速度特定手段３６によって、ステップＳ８_１、Ｓ８_２、…、Ｓ８_ｎで演算された誤差が最小となる直線の傾きを、対象物Ｂの移動速度に対応する位相差の傾きとして選択する（ステップＳ９）。
以上の動作によって、パルスドプラ計測装置１は、広帯域の周波数を用いて計測を行うため、対象物Ｂの動きにより位相差において位相の折り返しが発生する場合であっても、正しく対象物Ｂの速度を特定することができる。また、パルスドプラ計測装置１は、個々の周波数においてブラインド速度が発生しても、計測する周波数領域において、ブラインド速度が発生しないため、推定可能な速度領域を従来よりも拡大させることができる。

［シミュレーション結果］
最後に、従来の手法（従来のパルスドプラ法、２周波法）と、本発明に係るパルスドプラ計測装置１による広帯域ドプラ法との速度推定の評価結果について説明する。
ここでは、一次元多数散乱体モデルを用いて、速度推定のシミュレーションを行った。シミュレーションの条件は、以下の表に示す通りである。

なお、本シミュレーションにおいては、従来の手法と広帯域ドプラ法とでは、使用する送波パルスのパルス長のみを変え、従来の手法においては２．４μｓ、広帯域ドプラ法では、１．０μｓとした。また、広帯域ドプラ法において、８１種類の位相修正パターンを使用した。

図８は、各手法におけるシミュレーション結果をグラフ化したものであり、（ａ）は従来のパルスドプラ法、（ｂ）は従来の２周波法、（ｃ）は本発明に係る広帯域ドプラ法による結果を示している。なお、ここでは２周波法は、特許文献２に記載されている手法を用いている。また、図８の各グラフでは、横軸に速度測定対象である散乱体の予め設定した速度、縦軸にシミュレーションにより推定した速度を示している。理想的には、予め設定した速度と推定した速度とが一致し、プロット位置（図中「＋」）は、傾きが「１」の直線上に乗ることが望ましい。

図８（ａ）に示すように、従来のパルスドプラ法では、ナイキスト限界が０．４２ｍ／ｓとなるため、これを超える速度ではエイリアシングが発生している。また、散乱体からの反射波の干渉により推定した速度にばらつきが見られる。
また、図８（ｂ）に示すように、従来の２周波法では、エイリアシングは発生していないものの、干渉による位相変動の影響を受けるため、推定した速度に大きなばらつきが発生している。なお、このシミュレーション結果は、従来の２周期法にも同様のことが言える。
しかし、図８（ｃ）に示すように、本発明に係る広帯域ドプラ法では、エイリアシングが発生せず、さらに図８（ａ）の従来のパルスドプラ法に比べ、推定速度のばらつきが小さく、±１．４ｍ／ｓの範囲で精度よく速度推定がなされている。
このように、本発明に係る広帯域ドプラ法は、エイリアシングが発生せず、散乱体からの反射波の干渉の影響が少なく、精度よく速度推定を行うことができる。

本発明に係るパルスドプラ計測装置の構成を示すブロック図である。 本発明に係るパルスドプラ計測装置の動作を示すフローチャートである。 従来のパルスドプラ法と本発明に係る広帯域ドプラ法との違いを説明するためのグラフ図である。 広帯域ドプラ法におけるエイリアシングが発生した場合の位相差の修正を説明するためのグラフ図である。 ＭＴＩフィルタの原理を説明するための説明図である。 位相の折り返しが発生した場合における位相の修正の概念を示す概念図である。 位相修正パターンの一例を示すパターン図である。 シミュレーション結果を示すグラフ図である。 従来のパルスドプラ法の概念を説明するための説明図である。

符号の説明

１ パルスドプラ計測装置
２ 送受波装置
２１ 送波手段
２２ 受波手段
３ 制御装置
３１ フーリエ変換手段（スペクトル解析手段）
３２ ＭＴＩフィルタ（フィルタリング手段）
３３ 位相差演算手段
３４ 位相差傾き演算手段
３４１ 位相差修正手段
３４２ 直線近似手段
３５ 修正パターン記憶手段
３６ 速度特定手段

Claims (5)

1. 対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測装置であって、
   前記超音波を前記対象物に送波する送波手段と、
   前記対象物からの前記超音波に対する反射波を受波する受波手段と、
   この受波手段で受波した反射波から、周波数毎の位相を検出するスペクトル解析手段と、
   このスペクトル解析手段で検出された周波数毎の位相について、前記周期毎に位相差を演算する位相差演算手段と、
   前記周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算することで前記位相差を修正する位相差修正手段と、
   前記周波数と当該周波数における位相差とを座標点とする座標系において、前記複数の位相修正パターン毎に、前記周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する直線近似手段と、
   この直線近似手段で演算された前記誤差が最小となる直線の傾きを、前記移動速度に対応する位相差の傾きとして選択する速度特定手段と、
   を備えていることを特徴とするパルスドプラ計測装置。
2. 前記位相差演算手段は、前記周波数毎の位相差を、予め定めた送受波回数毎に平均化することを特徴とする請求項１に記載のパルスドプラ計測装置。
3. 前記スペクトル解析手段で検出された位相に基づいて、クラッタ成分を除去するフィルタリング手段を備えていることを特徴とする請求項１又は請求項２に記載のパルスドプラ計測装置。
4. 対象物に対して超音波を一定の周期で送波し、その反射波に基づいて、前記対象物の移動速度を計測するパルスドプラ計測方法であって、
   前記超音波を前記対象物に送波し、前記対象物からの反射波を受波する送受波ステップと、
   前記反射波から、周波数毎の位相を検出するスペクトル解析ステップと、
   このスペクトル解析手段で検出された周波数毎の位相について、前記周期毎に位相差を演算する位相差演算ステップと、
   前記周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算することで前記位相差を修正する位相差修正ステップと、
   前記周波数と当該周波数における位相差とを座標点とする座標系において、前記複数の位相修正パターン毎に、前記周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する直線近似ステップと、
   この直線近似ステップで演算された前記誤差が最小となる直線の傾きを、前記移動速度に対応する位相差の傾きとして選択する速度特定ステップと、
   を含んでいることを特徴とするパルスドプラ計測方法。
5. 対象物に対して一定の周期で送波された超音波に対する反射波に基づいて、前記対象物の移動速度を計測するために、コンピュータを、
   前記反射波から、周波数毎の位相を検出するスペクトル解析手段、
   このスペクトル解析手段で検出された周波数毎の位相について、前記周期毎に位相差を演算する位相差演算手段、
   前記周波数毎の位相差に対して、予め定めた異なる複数の位相修正パターンで２πの整数倍の位相を加算又は減算することで前記位相差を修正する位相差修正手段、
   前記周波数と当該周波数における位相差とを座標点とする座標系において、前記複数の位相修正パターン毎に、前記周波数毎の座標点を直線に近似することで、その直線の傾きと、その誤差とを演算する直線近似手段、
   この直線近似手段で演算された前記誤差が最小となる直線の傾きを、前記移動速度に対応する位相差の傾きとして選択する速度特定手段、
   として機能させることを特徴とするパルスドプラ計測プログラム。

Images