WO2006126303A1

WO2006126303A1 - 速度計測方法及びこれを用いた速度計測装置

Info

Publication number: WO2006126303A1
Application number: PCT/JP2006/300081
Authority: WO
Inventors: Shinichiro Umemura; Takashi Azuma; Tetsuya Hayashi; Naoyuki Murayama
Original assignee: Hitachi Medical Corporation
Priority date: 2005-05-27
Filing date: 2006-01-06
Publication date: 2006-11-30
Also published as: EP1884196A4; US20090177091A1; CN101184443A; JP5108512B2; EP1884196A1; CN101184443B; US7946992B2; JPWO2006126303A1

Abstract

　Ｎ個の時系列信号を０次から（Ｎ－１）次までの離散的ルジャンドル関数を基底として展開するステップ（Ｓ４）と、（２ｎ－１）次の展開係数と（２ｎ＋１）次の展開係数との線形結合に虚数単位を乗じて、２ｎ次の展開係数と線形結合した２ｎ次の複素展開係数と、（２ｎ＋１）次の展開係数に虚数単位を乗じた後、２ｎ次の展開係数と（２ｎ＋２）次の展開係数と線形結合した（２ｎ＋１）次の複素展開係数とを算出するステップ（Ｓ５）と、各複素展開係数のうち最大の絶対値を持つ係数の次数ｍを決定する次数決定ステップ（Ｓ４）と、次数ｍに対応した各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出するステップ（Ｓ８）とからなる速度計測方法及びこの方法を実施する速度計測装置。

Description

明細書

速度計測方法及びこれを用いた速度計測装置

技術分野

[0001] 本発明は、パルス状の反射波を用いた速度計測方法及びこれを用いた速度計測装置に関する。

背景技術

[0002] パルス 'ドッブラ速度計では、複数のパルス波を運動反射物に送信し、図 1のような受信された複数のエコー信号が解析される。すなわち、各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べた時系列信号が解祈され、運動反射物に関する速度などの情報が得られる。最も一般的な信号処理方法としては、これらの受信信号を直交検波し、複素時系列信号として解析する。このようなパルス 'ドップラ速度計は、生体内部の血流などを検出 ·描画する超音波診断装置、雨雲などを検出'描画する気象レーダー、飛行物体を検出する航空レーダーなどとして広く用いられている。

もし、反射物が 1つしか存在しなければ、前記の時系列信号の位相回転速度 Δ Φ Z A tすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度 Vを、次式のように容易に求めることができる。

ν= λ Δ Φ / Μ/2 (1)

ここで、 λは波長、 A tはパルス送信の時間間隔、 Δ Φは位相回転角である。運動速度 Vの符号は、反射物が近づく場合には、運動速度 Vの符号が正となり、遠ざかる場合には負となる。

[0003] ところで、前記のような実際の受信エコー中では、静止反射物によるエコー信号、いわゆるクラッタ信号の強度が、着目する運動反射物によるエコー信号の強度よりも、一般に、数桁以上大きい。また、静止反射物によるエコー信号は、実際には、時間軸上で完全に静止しているわけではなぐ途中の媒質のゆらぎや、医療用超音波診断装置の場合には、血流以外の静止臓器の遅い動きのために、ドリフトしている。このため、静止反射物によるエコー信号の周波数成分は、位相回転速度ゼロの直流成分だけでな位相回転速度がゼロでない低周波成分を有している。

[0004] このため、前記時系列信号に単純に（1)式の処理を適用しただけでは、運動反射物の速度を検出することはできない。そこで、実際のノ^レス'ドッブラ速度計は、静止反射物によるエコー信号を抑圧し、運動反射物によるエコー信号を相対的に強調する MTI(Moving Target Indicator)処理とよばれる信号処理を行った後に、速度の検出あるいは分析を行うように構成されてレ、る。

MTI処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点、がある。

1) N1個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出'分析処理部に入力されるデータ点が（N1— 1)個目減りする。

2) 急峻な遮断特性のフィルタを得にくい。

[0005] パルス 'ドッブラ速度計の場合は、運動反射物によるエコー信号を保ちながら、静止反射物のドリフト成分を除去する急峻な遮断特性が必要であり、特に、 2)が問題となる。また、 N回の送受信により得られる N点の時系列信号から速度検出 ·分析処理をするとき、データ点が目減りする 1)の問題は、実際に速度検出 ·分析演算に用い得るデータ点が (N— N1 + 1)点に減ってしまうことを意味する。これは、医療診断用超音波血流描画装置のように実時間性が重要な応用では、望ましレ、ことではなレ、。

[0006] この問題を解決する処理としては、非特許文献 1に記載されているように Polynomial

Regression Filter (多項式フィット 'フィルタ）が提案されている。これは、時系列信号に、 0次式（定数）、 1次式、 2次式、 · · · ·、 M次式を順次最小 2乗フィットし、フィットした成分を引き去ることにより、元の時系列信号のもつドリフト成分を除去する処理である。この処理は、 N点の入力時系列信号に N X Nの正方行列を掛けることで表現される。このため、出力信号として同じ N点の時系列信号が得られ、

1) データ点の目減りがない。

また、その低域遮断特性は、フィットする最大次数 Mの大きさによるが、

2) 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。非特許文献 1： IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Co ntrol, Vol. 42., 927- 937頁 (1995年)

発明の開示

[0007] ところが、この多項式フィット 'フィルタによる処理の後段において、 ν= λ Δ Φ/ A t /2式を用いて運動物体の速度の検出'分析を行うと、遮断周波数の近くで速度演算誤差が、低域遮断フィルタと比較しても大きくなつてしまうという問題点がある。このため、多項式フィット ·フィルタが有する急峻な遮断特性を十分に使うことができなレ、。そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。

[0008] 本発明は、速度計測対象物に対して Ν個のパルス波を所定間隔で送信し、受信ェコー信号を送信時刻の順番に並べた Ν個の時系列信号を 0次から（Ν— 1)次までの離散的ルジャンドル関数を基底として展開し、この展開係数を用いて速度を算出する。特に、 1次以上 (N—1)次以下の（2η_ 1)次の展開係数と（2η+ 1)次の展開係数との線形結合に虚数単位を乗じて、 2η次の展開係数と線形結合した 2η次の複素展開係数と、（2η+ 1)次の展開係数に虚数単位を乗じた後、 2η次の展開係数と（2η+ 2)次の展開係数と線形結合した（2η+ 1)次の複素展開係数とを算出し、各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数 mを決定し、次数 mに対応した前記各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出する。

[0009] パルス波が照射された運動反射物からは位相回転速度（角周波数） Δ Φ/ Δ ΐを有したエコー信号が得られ、ドリフト振動を伴った静止反射物からは低周波成分を有したクラッタ信号が得られる。このクラッタ信号の強度は、エコー信号の強度よりも数桁大きい。エコー信号を次数が大きいほど周波数原点付近で急速にゼロに収束する特性を有する離散的ルジャンドル関数を基底として展開し、各展開係数の横軸を位相回転速度（角周波数）に対応させることによって、運動反射物のエコー信号と静止対象物の低周波成分とが区別される。つまり、静止反射物と運動反射物との反射信号を区別することができる。 [0010] 本発明によれば、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減できる。

図面の簡単な説明

[0011] [図 1]N回のパルス送信によって反射した受信信号を並べた図である。

[図 2]離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。

[図 3]離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。

[図 4]離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。

[図 5]離散的フーリエ展開を用いた速度分析による速度スペクトル係数の例を示した図である。

[図 6]従来の低域遮断型 MTIフィルタと多項式フィット 'フィルタのもつ通過振幅特性の典型例を示す図である。

[図 7]従来の低域遮断型 MTIフィルタと多項式フィット 'フィルタをそれぞれ通した後、離散的フーリエ展開を用いた速度分析を行った場合に得られる速度スペクトル係数の例を示す図である。

[図 8]—連の 8点からなる信号に関する離散的ルジャンドル関数の偶数次を示す図である。

[図 9]一連の 8点からなる信号に関する離散的ルジャンドル関数の奇数次を示す図である。

[図 10]ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 11]ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 12]ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 13]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 14]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 15]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 16]ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 17]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 18]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 19]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 20]複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。

[図 21]ルジャンドル展開係数の自乗和比である。

[図 22]ルジャンドル係数の展開係数の自乗和比より求めた速度校正曲線である。

[図 23]本実施形態の速度計測装置を用いた超音波診断装置の構成を示すブロック図である。

[図 24]本実施形態のドッブラ速度検出アルゴリズムの例を示す図である。

[図 25]本実施形態による速度検出結果の例を示す図である。

[図 26]従来方式による速度検出結果の例を示す図である。

園 27]多項式フィット 'フィルタ方式による速度検出結果の例を示す図である。

[図 28]本実施形態による速度検出結果の例を示す図である。

[図 29]従来方式による速度検出結果の例を示す図。

園 30]多項式フィット 'フィルタ方式による速度検出結果の例を示す図である。

[図 31]本実施形態による速度検出結果の例を示す図である。

[図 32]従来方式による速度検出結果の例を示す図である。 [図 33]多項式フィット 'フィルタ方式による速度検出結果の例を示す図である。

[図 34]本実施形態のドッブラ速度検出方式による検出結果の一の例を示す図である

[図 35]本実施形態のドッブラ速度検出方式による検出結果の他の例を示す図である

[図 36]従来方式による速度検出結果の例を示す図である。

[図 37]多項式フィット 'フィルタ方式による速度検出結果の例を示す図である。

発明を実施するための最良の形態

[0012] (フィルタの原理）

各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べた N点の離散的時系列信号力速度分析を行う場合を考える。このような、パルス 'ドッブラ速度分析用として典型的な時系列複素信号は、受信されたェコ一信号に 90° 位相の異なる 2つの搬送周波数信号を乗じて得られる 1組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる（図 1参照）。

このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、 n= 1， 2, · · · , Nとならんだ時系列信号を、周波数を表す指数 k= _NZ2， · ·， 0, · · , N/2の複素正弦波系

Cs(n， k) = cos [ k π (2n_N_ 1)/N ] + j sin [ k π (2n_N_ 1)/N ] (2) で展開する。 jは虚数単位である。

負の kは反射物が遠ざかる向きの速度に対応し、正の kは近づく向きの速度に対応する。 k= ±N/2はナイキスト限界であり、正負が区別されないので、（2)式により表される（N + 1)個の複素正弦波関数のうち独立なものは N個となり、これらは、直交関数系を形成している。

[0013] 反射物が 1つしか存在せず、それがナイキスト限界範囲内の一定速度で運動している場合を想定し、運動速度の関数として、前記展開により得られる離散的フーリエ展開係数を計算し、 N = 8の場合について、その絶対値を図 2にプロットした。横軸は、運動速度に対応する位相回転速度を、パルス送信の時間間隔の逆数 Pulse Repet ition Frequency (PRF)の 2 π倍を単位として表示した。図 2では展開係数を表す 8個の関数が重なって見づらいので、そのうち 4個ずつを見分けやすい線により図 3及び図 4に示した。

[0014] 送信時間間隔の Ν倍の時間で、位相が整数回回転するような運動速度では、ゼロでない展開係数はただ 1つであり、他の展開係数はすべてゼロとなり、きれいに速度分析されている。ところが、それ以外の位相回転が整数回とならない一般の運動速度では、それぞれの展開係数が、横軸上ピークからはるかに離れた運動速度においても、ピークの 1/10程度の大きさの絶対値をもっている。これは、速度分析計として— 20dB程度の遠隔的クロストークをもつことを意味する。検出'分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号がドリフトするとき、これはクラッタ信号と運動反射物によるエコー信号とを区別できないという重大な問題を生じさせる。

[0015] この問題は、 N個の時系列信号に、ハニング関数のような緩やかに立ち上がって、下がる窓関数による重みをかけることにより、力なり抑圧できる力本質的には解決しない。また、この方法は、時系列信号の数が大幅に目減りしたのと同様の、望ましくない効果を生じ、特に、時系列信号の数が元々少ない場合には適さない方法である。

[0016] 具体例として、ナイキスト限界の 2/3の一定速度で運動する検出 ·分析対象である運動反射物があり、その 1000倍のエコー信号振幅をもつクラッタ信号が、ナイキスト限界速度の 1/100でドリフトしている場合を考える。図 5には、このときの前記速度分析計の出力である速度スペクトルを示す。検出'分析対象である運動反射物のスベクトルは、速度ゼロ付近にピークをもつクラッタ信号からのクロストーク成分に覆い力べされてしまレ、、この出力スペクトルの中では検出不可能である。このこと力まさに、クラッタ信号を遮断する MTIフィルタが必要な最も本質的な理由である。

[0017] 次に、 MTI処理として最もよく知られ、時間領域において畳み込み積和で表現される通常の低域遮断フィルタについて、その通過振幅特性の典型例を、図 6中に実線で示す。このようなフィルタに nl = 1 , · · · , (N + M- 1) とならんだ時系列信号 Sl(n 1)を、入力すると、出カとして₁10 = 1, ' ' '， Nとならんだ時系列信号 SO(nO)が得られる。

S0(n0)=∑ F (m) Sl(nO+m- 1) (3)

ここで、∑は m=l， ···， Mの和を表す。図 6では、低域遮断フィルタとして最も簡単な M = 3， F(l)=-1, F(2)=2, F (3) = _ 1の場合について実線で示した。図 5 と同様の入力信号を、この低域遮断フィルタを通した後、図 5と同様に速度分析すると、図 7中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が 1 Z2000倍程度に抑圧された結果、ナイキスト限界の 0.3倍から 0.4倍の速度にピークをもつ運動反射物のスペクトルが見えるようになってレ、る。

[0018] 同じ目的の MTIフィルタとして、前記の多項式フィット 'フィルタを設計したときの通過振幅特性の典型例を、図 6中に破線で示す。この例では、 0次から 3次式までをフィットして引き去った。前記の畳み込み型フィルタよりも急峻な低域遮断特性が見てとれる。また、図 5と同様の入力信号を、この低域遮断フィルタを通した後、図 5と同様に速度分析した結果を、図 7中に破線で示した。低域遮断特性の急峻さのために、クラッタ信号が完全に抑圧されている。多項式フィット 'フィルタ特有の急峻な低域遮断特性は以下のように理解することができる。

[0019] 連続座標系における狭義ルジャンドル多項式については、岩波数学公式 III特殊関数 82— 85頁に記述されており、例えば、 8次まで記すと

Ρ j =x = cos θ

1

P (x) = (l/2) (3x²-l) = (l/4) (3cos2 θ +1)

2

Ρ (χ) = (1/2) (5χ³-3χ) = (1/8) (5cos3 θ +3cos θ )

3

Ρ (χ) = (1/8) (35χ⁴-30χ + 3) = (1/64) (35cos4 θ +20cos2 θ +9)

4

Ρ (χ) = (ΐ/8) (63χ⁵-70χ +15χ) = (1/128) (63cos5 θ +35cos3 θ +30c

5

OS θ )

Ρ (Χ) = (1/16) (231χ⁶-315χ⁴+105χ²-5) = (1/512) (231cos6 θ +126

6

cos4 θ +105cos2 θ +50)

Ρ (χ) = (1/16) (429χ^?-639χ⁵ + 315χ³-35χ) = (1/2¹⁰) (429cos7 θ +23 lcos5 θ +189cos3 θ +175cos θ ) P (x) = (1/128) (6435x⁸- 12012x⁶ + 6930x⁴- 1260x² + 35) = (l/2¹⁴) (6

8

435cos8 Θ + 3432cos6 Θ + 2772cos4 Θ + 2520cos2 Θ + 1225)

である。離散的座標系についても使えるように、一般化すると、偶数（2n)次ルジヤンドル多項式は、 2nより低レ、次数のルジャンドル多項式すべてと直交する最高次数 2n の偶関数と定義することができ、奇数（2n+ l)次ルジャンドル多項式は、（2n+ l)より低い次数のルジャンドル多項式すべてと直交する最高次数 (2n+ l)の奇関数と定義すること力 Sできる。

[0020] 例として、一連の 8点からなる時系列信号について、 0次から 7次までの離散的ルジヤンドル関数を計算し、偶数次を図 8に、奇数次を図 9に示した。 0次 · 1次を実線、 2 次 · 3次を破線、 4次 · 5次を 1点鎖線、 6次 · 7次を点線にて示した。振幅は、 RMS (2 乗平均の平方根)値が 1となるよう規格化した。符号は、最低次項の符号が正となるよう、すなわち、偶数次多項式については定数項、奇数次多項式については 1次項の符号が正となるよう定めた。図 10には、その周波数スペクトルをプロットした。周波数成分の振幅の絶対値を、個々のルジャンドル展開係数の最大値で規格化して示した。図 10では 8個の関数が重なってしまって見づらいので、そのうち 4個ずつを見分けやすい線により図 11及び図 12に示した。 A (0) , A (7)を実線、 A(l) , A (6)を破線、 A (2) , A (5)を 1点鎖線、 A (3)， A (4)を点線にて示した。

[0021] 図 10から図 12を、図 2から図 4と同様、入力複素正弦波に対する周波数応答として眺めてみる。ピークを主応答、それ以外を副応答とみる。ルジャンドル関数は、図 2から図 4のフーリエ展開系に比べ、主応答より高周波側の副応答は大きいが、低周波側には副応答をもたず、次数が大きいほど周波数原点付近で急速にゼロに収束する。ここでは解析的証明は省略する力周波数の関数としての m次ルジャンドル関数は、周波数原点において 0次から (m_ l)次導関数まですべてが 0となる。

[0022] 前記の多項式フィット 'フィルタは、入力信号を、この離散的ルジャンドル関数を基底として展開し、 0から M次までを除き、（M + 1)次以上の成分の和を出力信号とするフィルタであるということができる。したがって、出力信号は、周波数原点付近で周波数の (M + 1)乗に比例する挙動でゼロに収束していく。図 6及び図 7中に破線で示された例のように、多項式フィット 'フィルタの低域遮断特性が優れているのは、このようなルジャンドル関数の本質的な特性に由来するものである。このように考えると、クラッタという周波数原点付近に桁違いに大きな振幅をもつ妨害成分の影響を抑えて、着目する運動物体を検出したり、その速度を分析したりしなければならないドッブラ速度検出には、多項式フィット 'フィルタやルジャンドル関数展開を利用する方式は、本質的に適しているといえる。

[0023] フーリエ展開系では、複素正弦波を基底とすることにより、図 2から図 5や図 7に示したように符号つきの速度分析が可能であることは広く知られており、そのような分析方法が広く用いられている。ところ力ルジャンドル展開系は、前記のように、ドッブラ速度検出に極めて適した性質をもつものの、このまま用いたのでは、図 10力ら図 20より見てとれるように、絶対値が等しいが符号の異なる速度に対して等しい応答を示すので、符号つきの速度分析は不可能である。

[0024] そこで、次数が 1つだけ異なる偶数次と奇数次のルジャンドル関数を虚数単位 jを係数として線型結合して、複素ルジャンドル関数を形成し、これにより符号つきの速度分析を可能とする。これは、偶数次のルジャンドル関数が余弦波に、奇数次のルジヤンドル関数が正弦波に似た変化を示すことに着目し、余弦関数と正弦関数を虚数単位 jを係数として線型結合することにより、位相角を j倍した指数関数、すなわち複素正弦関数が得られ、この位相角の増減から符号つき速度が得られるのと同様に、複素ルジャンドル関数から符号つき速度を得ようとするものである。

具体的には、一連の N点からなる時系列信号を、まず、 0次から（N—1)次の離散的ルジャンドル関数を基底として展開して、 A(0)から A (N_ 1)までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数

C (±(2n_ 3Z2)) =A (2n- 2) ±jA (2n_ l) (l≤2n_ l≤N_l)

C (±(2n_ lZ2)) =A (2n) ±j A (2n_ l) (2≤2n≤N_ l;n:自然数,複号同順） (4)

を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。

[0025] 図 13力図 15には、図 10から図 12と同様 N = 8の場合について、複素正弦波入力に対する出力としての複素ルジャンドル係数の絶対値を、符号つき速度に対応する位相回転速度の関数として示した。複素ルジャンドル係数の絶対値は、それぞれの最大値で規格化してプロットした。図 13中の実線は C(6.5)、破線は C(_6.5)、図 14中の実線は C (5.5)、破線は C(_5.5)、図 15中の実線は C (4.5)、破線は C(-4.5)である。これらの図から、複素ルジャンドル係数が、速度に対し正負を区別して応答していることがわかる。なお、このような速度に対する応答をもつ複素ルジヤンドル系が N≤ 35の範囲で構築できることが確認されてレ、る。

[0026] 図 13から図 15に示した複素ルジャンドル係数の応答は、ルジャンドル展開係数の特徴である位相回転速度原点すなわちドッブラ周波数原点のまわりの極めて優れたカットオフ特性を示している反面、図 2から図 4に示したフーリエ展開係数の応答と比較して、周波数サイドローブレベルが高めである。すなわち、周波数サイドローブレべルは、フーリエ展開係数の図 2から図 4の例では、 0.2程度であるのに対し、複素ルジヤンドル係数の図 14の例では、 0.6近くに達してレヽる。周波数サイドローブレベルが高いことは、 S/N比すなわちノイズレベルに対する信号レベルが低い場合に、速度検出誤差を生む可能性が高いことを意味するので、周波数サイドローブレベルは抑圧することが望ましい。

[0027] そこで、 A(2n)をもとに C(±2n)を生成する際、 A(2n)と同じ周波数において最大値をとるよう A(2n—1)と A(2n—1)とを線型結合して用レ、、また、 A(2n— 1)をもとに C(±(2n— 1))を生成する際、 A(2n— 1)と同じ周波数において最大値をとるよう A( 2n-2)と A(2n)とを線型結合して用いる。

すなわち、

C(±(2n_l)) =ひ（2n_l)XA(2n_2) + j3 (2n- 1) X A(2n) ±j A(2n_l) (l≤2n-l≤N-l)

C(±2n) =A(2n) ±j [ひ（2n)XA(2n— 1)+ j3 (2n) X A(2n+ 1) ] (2≤2n≤ N— l;n:自然数，複号同順；ひ及び /3：正の実数）（5)

となる。

[0028] 図 16から図 17には、図 13から図 15と同様 N = 8の場合について、ルジャンドル係数 A (4)、 A (5)、 A (6)から複素ルジャンドル係数 C(± 5)を生成する例を示した。図 16中の実線は A(5)の絶対値のドッブラ周波数応答である。一方、点線は、その A( 5)と同じ周波数において最大値をとるように A (4)と A(6)とを線型結合して生成した係数の絶対値のドッブラ周波数応答である。図 17は、その A (4)と A(6)との線型結合を実部、 A (5)を虚部として得られる複素ルジャンドル係数 C (± 5)の絶対値のドッブラ周波数応答である。図中実線は C (5)であり、破線は C (_ 5)である。周波数サイドローブレベルは、 0. 4程度に抑えられている。これは、図 14に示された C (± 5. 5) や図 15に示された C ( ±4. 5)の周波数サイドローブレベル 0. 6程度に比べ、大きく改善されている。同様に生成した複素ルジャンドル係数 C (± 6)、 C ( ±4)、 C ( ± 3) の絶対値のドッブラ周波数応答を図 18から図 20に示した。これらは図 17と同様に、実線が C (6)、 C (4)、 C (3)であり、破線が C (一 6)、 C (一 4)、 C (一 3)である。周波数サイドローブレベルは、やはり 0· 4程度に抑えられている。すなわち、 A (2n— 1)と同じ角周波数において最大値をとるよう A (2n— 2)と A (2n)とを線型結合することによって、サイドローブレベルが 0. 6力ら 0. 4程度に抑えられる。

[0029] 図 21には、 mを自然数とするとき、（m— 1)次から (N— 1)次までのルジャンドル係数の自乗和と m次から (N— 1)次までのルジャンドル係数の自乗和との比

R (m) = [A(m) ²+ · · · +A (N - 1)²] / [A(m— 1) ²+A (m) ²+ · · · +A(N - 1)² ] (l≤m≤N- l) (6)

を、前記の例と同様 N = 8の場合について示した。これらの曲線は、ドッブラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップノレする。その単調増加部分を用いることにより、図 22に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドッブラ速度を算出すればよい。

[0030] 検出 ·分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号が存在し、それがドリフトするときには、そのドリフトの程度が大きくなるに応じて、 C (± l)、 C ( ± 2)、……と、次数の小さな係数を順次除いた複素ルジャンドル係数系を用いて、速度検出'分析を行えばよい。このようにルジャンドル係数の遮断次数を制御することにより、クラッタ信号の影響を抑圧しながら、運動反射物の速度検出- 分析を行うことができる。

[0031] (第 1実施形態）以下、本発明の実施形態について説明する。

図 23は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子 1を構成する各素子は、切り替えスィッチ群 2を介して、送波ビームフォーマ 3と受波ビームフォーマ 10に接続されている。送波ビームフォーマ 3では、送受信シークェンス制御部 6による制御に従って、送信波形メモリ 5から送信波形選択部 4により選択されて読み出された波形を用レ、、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子 1の各素子により超音波ノ^レスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子 1に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ 10では、送受信シークェンス制御部 6による制御に従って、指向性を持つ受信感度を生成すベぐ各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークェンス制御部 6による制御に従つて、受波メモリ選択部 11により選択された受波信号メモリ 12中の 1つのバンクへ一且書き込まれ、ドッブラ信号分析すべき N個の時系列信号がそろったのちに読み出されて、位相回転検出器 13， 15、ダウン'ミキシング器 14、血流信号検出'分析部 1 6において速度を検出'分析すべく信号処理される。

読み出された N個の時系列受信信号は、各々を得るために行われたパルス送信の時刻を基準とする時相の等しい部分を、送信時刻の順番に並べたのちに処理される。この速度検出アルゴリズムについて、図 24のフローチャート及び図 23を参照して説明する。まず、位相回転検出器 13では、クラッタ信号を含んだままの時系列信号から位相回転を求める。位相回転を検出する最も典型的な方法は η= 1，· · ·， Νとならんだ隣接する時系列複素信号 S (η)， S (η+ 1)をもとに、

P (n) = S (n+ l) S (n) * /||S (n+ 1) ||/||S (n) || (n= 1， · . · ,N_ 1)

(7)

を計算し、 P (n)の平均値 Paの位相から平均位相回転速度を求める。

ここで、 S (n) *は S (n)の複素共役、 ||S (n) IIは S (n)の絶対値である。一般に、クラッタ 'エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、 p_aは、概ねクラッタ 'エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される（図 24の Sl)。また、時系列複素信号 S (n)は、受信されたエコー信号に 90° 位相の異なる 2つの搬送周波数信号を乗じて得られる。

[0033] ダウン'ミキシング器 14では、ここで求めた平均値 Pa又は Paの周辺領域における空間平均値を用い、エコー信号の位相回転平均値がゼロとなるようにミキシング処理が行われる。いいかえれば、ダウン'ミキシングによりクラッタ平均速度による位相回転が打ち消される（図 24の S2参照）。すなわち、 S(n)と Paより

Sd(n)=S(n)Pa*ⁿ(n=l,---, N) (8)

を得る。このダウン'ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。

[0034] 得られた隣接サンプルである時系列複素信号 Sd(n) , Sd(n+ 1)をもとに、位相回転検出器 15では、位相回転検出器 13と同様に、

Pd(n)=Sd(n+l)Sd(n)V ||Sd(n+l)|| /||Sd(n)|| (η=1,···, Ν— 1)

(9)

を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される（図 24の S3)。この位相回転最大値の大きさに応じて、血流信号検出'分析部 16の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される（図 24の S4)。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数 Mがより高く設定され、クラッタ成分が効果的に抑圧される。

[0035] ここで、本実施形態の特徴部である血流信号検出'分析部 16の動作について、もう少し詳しく述べる。入力された時系列信号 Sd(l)，···， Sd(N) は、まず、図 8及び図 9のような 0から (N_ 1)次の離散的ルジャンドル関数に展開され、その展開係数 A(O) ,···, A(N-l)が求められる。この計算は、次のような行列演算により容易に行うことができる。 n次ルジヤンドノレ関数を行べクトノレ L(n)で表し、これを n=0から (N—1) までならベて NX N行歹 1JLLを生成すると、 Sd(l) ,···, Sd(N)を要素とする列ベクトル Sd力 Α(0)，· · ·，Α(Ν— 1)を要素とする列ベクトル Aを得るための行列は、

と求めることができる。ここで、は LLの転置行歹 IJ、 LL— ¹は LLの逆行歹を表す。この行列を予め準備しておけば、行列演算

A=F-Sd (11)

によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数 Α(0),·· ·, Α(Ν_1)から、式 (5)に従って、 C(±l)，···， C (土 (N—1))の複素ルジヤンドノレ係数が生成される（図 24の S5)。いいかえれば、 n次ルジャンドル関数を行ベクトル L (n)で表し、これを n = 0から（N— 1)まで列方向に並べた N X N行列 LLを転置した転置行列 tLLを生成する過程と、転置行列 tLLに N X N行列 LLを前から乗算したものの逆行列（LL · tLU— ¹を算出する過程と、 N X N行列 LLに逆行列（LL · tLL)— ¹を前力乗算する過程とによって算出された行歹 IJFFを要素が時系列信号である列べタトル Sd(n)に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数 (遮断次数 M )が決定され（図 24の S4)、 C(±l)から C (土 (M+1))までの低次数複素ルジヤンドル係数が棄却される。

[0036] 次に、棄却されなレ、で残された正次数の複素ルジャンドル係数の自乗和と負次数の複素ルジャンドル係数の自乗和とを比較して、速度符号が決定される（図 24の S6 参照）。具体的に、大きな方の次数の符号を速度の符号とする。さらに、速度符号次数の複素ルジャンドル係数のうち、絶対値最大の係数が検出される（図 24の S7参照 )。具体的に、正なら正の次数、負なら負の次数の複素ルジャンドル係数と、最大次数 (N— 1)のルジャンドル係数のうち、絶対値最大の係数の次数 mが検出される。

[0037] そして、検出したルジャンドル係数に対応する校正曲線を用いて血流速度を算出する（図 24の S8参照）。すなわち、次数 mに対応する図 22に示したような速度校正曲線を用いて、式 (6)のように (m_l)次から (N_l)次までのルジャンドル係数の自乗和と m次から (N_ 1)次までのルジャンドル係数の自乗和との比力ドッブラ速度を求める。

[0038] 得られた血流速度信号は、エコー振幅検出器 17及びエコー振幅圧縮器 18により得られる静止臓器からのエコー信号と共に、スキャンコンバータ 19に入力される。スキャンコンバータでは、入力された複数の信号を適宜重畳して表示器 20にて 2次元ないし 3次元表示すベぐ信号の生成 ·制御を行う。

[0039] 次に、本実施形態のドッブラ速度計測装置の動作例を従来方式と比較して示す。

クラッタ 'エコー信号の振幅は、血流エコーの 300倍とし、クラッタ速度については、初期値がゼロで一定の加速度で立ち上がる場合を想定した。時系列信号数 N = 8の場合について、数値計算シミュレーションを行って比較した。従来方式としては、 MT Iフィルタとして図 6の場合の低域遮断フィルタを用レ、、その後段において、式 (5)を用レ、て位相回転量平均値 Paを得たのと同様の方法で、速度を算出する処理を血流信号検出 ·分析部 16において行う例を選んだ。さらに、 MTIフィルタを、やはり図 5に例を示したような多項式フィット 'フィルタに置き換えた場合とも、比較した。この場合、遮断次数 Mは、本発明の方法と同様、位相回転検出器 15の出力信号に応じて制御した。

[0040] 図 25から図 27には、クラッタの最終到達速度がナイキスト限界速度の 0. 8%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、及び多項式フィット 'フィルタ方式による速度検出結果を示した。横軸に入力としての血流速度、縦軸に出力として検出された速度をとり、有効速度検出範囲を実線、範囲外を点線によりプロットした。また、理想的場合を一点鎖線により示した。このように、クラッタ速度が低い場合には、いずれの方式も、ほぼ正確な速度分析結果を与えているが、本実施形態の方式による誤差が最も小さレ、。

[0041] 図 28から図 30には、同様に、クラッタの最終到達速度がナイキスト限界速度の 3% である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット 'フィルタ方式による速度検出結果を示した。このクラッタ速度域では、低域遮断フィルタ方式による速度検出は、ほぼ破綻に至っている。また、多項式フィット' フィルタ方式は、有効速度検出範囲における誤差は大きくないものの、血流速度がそれより小さい場合の誤差が大き何らかの対策を加えない限り、使うことができなレ、。これに対し、本実施形態の方式は、有効速度検出範囲において正確な速度検出結果を与えるだけでなぐ血流速度がそれより小さい場合、素直に検出速度ゼロを出力している。このように、有効速度検出範囲において正確な速度検出結果を出力し、かつ、血流速度がそれより小さい場合、素直に、検出速度ゼロを出力する力、あるいは速度検出不能であることを表示することは、本実施形態の速度計測装置の特徴であるということができる。この特徴を用いれば、着目した速度計測装置について、装置内部の構成を仔細に検査するに及ぶことなぐ本実施形態の実施の有無を確認すること力 Sできる。

[0042] 図 31から図 33には、さらに、クラッタの最終到達速度がナイキスト限界速度の 20% である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット 'フィルタ方式による速度検出結果を示した。ここまでクラッタ速度が大きくなると、従来方式による速度検出は、完全に破綻している。また、多項式フィット'フィノレタ方式は、有効速度検出範囲における誤差が大きくなると同時に、血流速度がそれより小さい場合の誤差も著しく大きくなつてしまっている。これに対し、本実施形態の方式は、有効速度検出範囲がやや狭くなるものの、血流速度がナイキスト限界速度の 60%以上の範囲で、ほぼ正確な速度検出結果を与えるだけでなぐ血流速度がそれより未満の場合の動作が、やはり素直である。

[0043] 図 34から図 37には、血流速度自身がナイキスト限界速度の 8%に対応するランダムな位相誤差をもつ場合について、それぞれ、式 (5)を用いる本実施形態の方式、式 (4)を用いる方式、低域遮断フィルタ方式、多項式フィット 'フィルタ方式による速度検出結果を示した。クラッタの最終到達速度については、ナイキスト限界速度の 10%とした。血流速度自身が位相誤差をもっために、負側の低中速領域において、本実施形態の方式も破綻し力かっている力式 (4)を用いる方式では、速度の符号が誤って検出されているのに対し、式 (5)を用いる方式では、速度符号については、誤ることなく検出されている。これは、周波数サイドローブレベルを抑圧するよう構築された式 (5 )を用いる方式の効果である。多項式フィット 'フィルタ方式も、同じ負側の低中速領域において、速度検出が破綻しており、その他の速度域においても、かろうじて速度符号が検出できたにとどまつている。また、クラッタ速度が大きいため、従来方式による速度検出は完全に破綻している。

[0044] 以上述べたように本実施形態によれば、血流を、臓器の運動と峻別して描画することができる。この特徴をさらに生力て、図 23の超音波診断装置では、位相回転検出器 13により検出されたクラッタ信号すなわち臓器エコーの位相回転速度の信号をスキャンコンバータ 19に入力し、臓器の運動速度あるいはその空間微分の分布を示す像を血流像と重畳あるいは並べて表示可能としている。この構成の有用性を、肝腫瘍の場合を例として述べる。肝腫瘍の辺縁部は新生血管が発達し、血流動態が周辺の正常肝とは異なる。また、周辺の正常肝とは固さが異なるため組織の運動も周辺とは異なる。したがって、血流像に加えて臓器の運動速度の場所による変化を表示することは、肝腫瘍の診断にきわめて有用な画像を提供することになる。

[0045] また、生体内部の血流など運動する反射物の速度を、そのエコー信号に比べ桁違いに振幅が大きレ、クラッタ'エコー信号の影響を除レ、て、正確に検出することができる。具体的には、臓器の超音波探触子への方向の運動速度が、 lmm/s変化するなか、超音波探触子への方向の速度成分 3mm/s以上の血流をリアルタイムに描画できる。このように、本実施形態によれば、的確な医療診断を可能とする血流検出 · 描画機能つき超音波診断装置を提供することが可能となる。すなわち、本実施形態を実施した装置の医用診断上の有用性はきわめて大きぐ従って、医用診断を支える工業における本実施形態の意義も、また、大きい。さらに、本実施形態の方法は、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス ·ドップラ ·レーダー装置など、他のパルス ·ドッブラ装置の運動反射体検出能力をも飛躍的に向上させることができ、本実施形態の工業的ならびに社会的意義は、その意味でもきわめて大である。

[0046] (変形例）

本実施形態は前記した実施形態に限定されるものではなぐ例えば以下のような種々の変形が可能である。

(1)前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数 mを求め、次数 mを境界とした展開係数間の自乗和の比を用いてドッブラ速度を算出した力複素展開係数間の自乗和の比を用いることもできる。

(2)前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レー一などのパルス ·ドッブラ ·レーダー装置に用いることができる。

Claims

請求の範囲

[1] 速度計測対象物に対して N個のパルス波を所定間隔で送信し、前記パルス波の受信エコー信号を送信時刻の順番に並べた時系列信号を用いて前記速度計測対象物を構成する運動反射物の速度を計測する速度計測方法であって、

前記 N個の時系列信号を 0次から (N— 1)次までの離散的ルジャンドル関数を基底として展開する展開ステップと、

自然数 nとするとき、前記展開ステップで展開した 1次以上 (N— 1)次以下の（2n— 1)次の展開係数と（2n+ 1)次の展開係数との線形結合に虚数単位を乗じて、 2n次の展開係数と線形結合した 2n次の複素展開係数と、 (2n+ l )次の展開係数に虚数単位を乗じた後、 2n次の展開係数と（2n+ 2)次の展開係数と線形結合した（2n+ 1 )次の複素展開係数とを算出する複素展開係数算出ステップと、

前記各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数 mを決定する次数決定ステップと、

前記次数 mに対応した前記各展開係数間あるいは前記各複素展開係数間の自乗和の比から前記運動反射物に関する符号付き速度信号を算出する速度信号算出ステツプと

からなる速度計測方法。

[2] 前記展開ステップは、 n次ルジャンドル関数を行ベクトル L (n)で表し、これを n=0 から (N_ 1)まで列方向に並べた N X N行列 LLを転置した転置行列を生成する過程と、前記転置行列に前記 N X N行列 LLを前から乗算したものの逆行列（LL • D— ¹を算出する過程と、前記 N X N行歹 IJLLに前記逆行列（LI^LL)— ¹を前から乗算する過程とによって算出された行歹 IJFFを要素が前記時系列信号である列べクトノレ Sd (n)に前から乗算するステップであることを特徴とする請求の範囲第 1項に記載の速度計測方法。

[3] 静止反射物に対する前記受信エコーの大きさの変動あるいはドリフトの程度によつて、前記複素展開係数を小さな次数から徐々に無視することを特徴とする請求の範囲第 2項に記載の速度計測方法。

[4] 前記運動反射物の速度が有効速度検出範囲よりも小さい場合に検出速度ゼロを出力するか否かによって請求の範囲第 3項に記載の速度計測方法の実行を判定するステップを備えたことを特徴とする速度計測方法。

[5] 速度計測対象物に対して N個のパルス波を所定間隔で送信し、前記パルス波の受信エコー信号を送信時刻の順番に並べた時系列信号を用いて前記速度計測対象物を構成する運動反射物の速度を計測する速度計測装置であって、

前記 N個の時系列信号を 0次から (N— 1)次までの離散的ルジャンドル関数を基底として展開する展開手段と、

自然数 nとするとき、前記展開ステップで展開した 1次以上 (N— 1)次以下の（2n— 1)次の展開係数と（2n+ 1)次の展開係数との線形結合に虚数単位を乗じて、 2n次の展開係数と線形結合した 2n次の複素展開係数と、 (2n+ l)次の展開係数に虚数単位を乗じた後、 2n次の展開係数と（2n+ 2)次の展開係数と線形結合した（2n+ 1 )次の複素展開係数とを算出する複素展開係数算出手段と、

前記各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数 mを決定する次数決定手段と、

前記次数 mに対応した前記各展開係数間あるいは前記各複素展開係数間の自乗和の比から前記運動反射物に関する符号付き速度信号を算出する速度信号算出手段と

からなる速度計測装置。

[6] 前記展開手段は、 n次ルジャンドル関数を行ベクトル L (n)で表し、これを n= 0から

(N— 1)まで列方向に並べた N X N行列 LLを転置した転置行列を生成する過程と、前記転置行列 ⁴LLに前記 N X N行列 LLを前から乗算したものの逆行列（LL ^LL を算出する過程と、前記 N X N行歹 IJLLに前記逆行列 (LL 'tLL)— ¹を前力乗算する過程とによって算出された行歹 IJFFを要素が前記時系列信号である列ベクトル Sd ( n)に前から乗算する手段であることを特徴とする請求の範囲第 5項に記載の速度計測装置。

[7] 静止反射物に対する前記受信エコーの大きさの変動あるいはドリフトの程度によつて、前記複素展開係数を小さな次数から徐々に無視することを特徴とする請求の範囲第 6項に記載の速度計測装置。

[8] 前記パルス波は、超音波であり、

前記運動反射物は、血流であることを特徴とする請求の範囲第 5項ないし第 7項に記載の速度計測装置。

[9] 前記パルス波は、超音波であり、

前記速度計測対象物は、肝臓であり、

前記肝臓の速度の空間微分の分布及び血流像を表示する表示手段を備えたことを特徴とする請求の範囲第 5項ないし 7項に記載の速度計測装置。