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

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

Info

Publication number
JPWO2006126303A1
JPWO2006126303A1 JP2007517722A JP2007517722A JPWO2006126303A1 JP WO2006126303 A1 JPWO2006126303 A1 JP WO2006126303A1 JP 2007517722 A JP2007517722 A JP 2007517722A JP 2007517722 A JP2007517722 A JP 2007517722A JP WO2006126303 A1 JPWO2006126303 A1 JP WO2006126303A1
Authority
JP
Japan
Prior art keywords
order
speed
expansion coefficient
complex
legendre
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.)
Granted
Application number
JP2007517722A
Other languages
English (en)
Other versions
JP5108512B2 (ja
Inventor
梅村 晋一郎
晋一郎 梅村
東 隆
隆 東
哲矢 林
哲矢 林
村山 直之
直之 村山
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hitachi Healthcare Manufacturing Ltd
Original Assignee
Hitachi Medical Corp
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hitachi Medical Corp filed Critical Hitachi Medical Corp
Priority to JP2007517722A priority Critical patent/JP5108512B2/ja
Publication of JPWO2006126303A1 publication Critical patent/JPWO2006126303A1/ja
Application granted granted Critical
Publication of JP5108512B2 publication Critical patent/JP5108512B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88Sonar systems specially adapted for specific applications
    • G01S15/89Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8979Combined Doppler and pulse-echo imaging systems
    • G01S15/8981Discriminating between fixed and moving objects or between objects moving at different speeds, e.g. wall clutter filter
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/06Measuring blood flow
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/48Diagnostic techniques
    • A61B8/488Diagnostic techniques involving Doppler signals
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88Sonar systems specially adapted for specific applications
    • G01S15/89Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8977Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Molecular Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Acoustics & Sound (AREA)
  • Biophysics (AREA)
  • Hematology (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

N個の時系列信号を0次から(N−1)次までの離散的ルジャンドル関数を基底として展開するステップ(S4)と、(2n−1)次の展開係数と(2n+1)次の展開係数との線形結合に虚数単位を乗じて、2n次の展開係数と線形結合した2n次の複素展開係数と、(2n+1)次の展開係数に虚数単位を乗じた後、2n次の展開係数と(2n+2)次の展開係数と線形結合した(2n+1)次の複素展開係数とを算出するステップ(S5)と、各複素展開係数のうち最大の絶対値を持つ係数の次数mを決定する次数決定ステップ(S4)と、次数mに対応した各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出するステップ(S8)とからなる速度計測方法及びこの方法を実施する速度計測装置。

Description

本発明は、パルス状の反射波を用いた速度計測方法及びこれを用いた速度計測装置に関する。
パルス・ドップラ速度計では、複数のパルス波を運動反射物に送信し、図1のような受信された複数のエコー信号が解析される。すなわち、各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べた時系列信号が解析され、運動反射物に関する速度などの情報が得られる。最も一般的な信号処理方法としては、これらの受信信号を直交検波し、複素時系列信号として解析する。このようなパルス・ドップラ速度計は、生体内部の血流などを検出・描画する超音波診断装置、雨雲などを検出・描画する気象レーダー、飛行物体を検出する航空レーダーなどとして広く用いられている。
もし、反射物が1つしか存在しなければ、前記の時系列信号の位相回転速度 ΔΦ/Δtすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度vを、次式のように容易に求めることができる。
v=λΔΦ/Δt/2 (1)
ここで、λは波長、Δtはパルス送信の時間間隔、ΔΦは位相回転角である。運動速度vの符号は、反射物が近づく場合には、運動速度vの符号が正となり、遠ざかる場合には負となる。
ところで、前記のような実際の受信エコー中では、静止反射物によるエコー信号、いわゆるクラッタ信号の強度が、着目する運動反射物によるエコー信号の強度よりも、一般に、数桁以上大きい。また、静止反射物によるエコー信号は、実際には、時間軸上で完全に静止しているわけではなく、途中の媒質のゆらぎや、医療用超音波診断装置の場合には、血流以外の静止臓器の遅い動きのために、ドリフトしている。このため、静止反射物によるエコー信号の周波数成分は、位相回転速度ゼロの直流成分だけでなく、位相回転速度がゼロでない低周波成分を有している。
このため、前記時系列信号に単純に (1)式の処理を適用しただけでは、運動反射物の速度を検出することはできない。そこで、実際のパルス・ドップラ速度計は、静止反射物によるエコー信号を抑圧し、運動反射物によるエコー信号を相対的に強調するMTI(Moving Target Indicator) 処理とよばれる信号処理を行った後に、速度の検出あるいは分析を行うように構成されている。
MTI処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点がある。
1) N1個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出・分析処理部に入力されるデータ点が (N1−1) 個目減りする。
2) 急峻な遮断特性のフィルタを得にくい。
パルス・ドップラ速度計の場合は、運動反射物によるエコー信号を保ちながら、静止反射物のドリフト成分を除去する急峻な遮断特性が必要であり、特に、2)が問題となる。また、N回の送受信により得られるN点の時系列信号から速度検出・分析処理をするとき、データ点が目減りする1)の問題は、実際に速度検出・分析演算に用い得るデータ点が(N−N1+1) 点に減ってしまうことを意味する。これは、医療診断用超音波血流描画装置のように実時間性が重要な応用では、望ましいことではない。
この問題を解決する処理としては、非特許文献1に記載されているようにPolynomial Regression Filter (多項式フィット・フィルタ)が提案されている。これは、時系列信号に、0次式(定数)、1次式、2次式、・・・・、M次式を順次最小2乗フィットし、フィットした成分を引き去ることにより、元の時系列信号のもつドリフト成分を除去する処理である。この処理は、N点の入力時系列信号にN×Nの正方行列を掛けることで表現される。このため、出力信号として同じN点の時系列信号が得られ、
1) データ点の目減りがない。
また、その低域遮断特性は、フィットする最大次数Mの大きさによるが、
2) 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 42., 927-937頁 (1995年)
ところが、この多項式フィット・フィルタによる処理の後段において、v=λΔΦ/Δt/2式を用いて運動物体の速度の検出・分析を行うと、遮断周波数の近くで速度演算誤差が、低域遮断フィルタと比較しても大きくなってしまうという問題点がある。このため、多項式フィット・フィルタが有する急峻な遮断特性を十分に使うことができない。
そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。
本発明は、速度計測対象物に対してN個のパルス波を所定間隔で送信し、受信エコー信号を送信時刻の順番に並べたN個の時系列信号を0次から(N−1)次までの離散的ルジャンドル関数を基底として展開し、この展開係数を用いて速度を算出する。特に、1次以上(N−1)次以下の(2n−1)次の展開係数と(2n+1)次の展開係数との線形結合に虚数単位を乗じて、2n次の展開係数と線形結合した2n次の複素展開係数と、(2n+1)次の展開係数に虚数単位を乗じた後、2n次の展開係数と(2n+2)次の展開係数と線形結合した(2n+1)次の複素展開係数とを算出し、各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数mを決定し、次数mに対応した前記各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出する。
パルス波が照射された運動反射物からは位相回転速度(角周波数)ΔΦ/Δtを有したエコー信号が得られ、ドリフト振動を伴った静止反射物からは低周波成分を有したクラッタ信号が得られる。このクラッタ信号の強度は、エコー信号の強度よりも数桁大きい。エコー信号を次数が大きいほど周波数原点付近で急速にゼロに収束する特性を有する離散的ルジャンドル関数を基底として展開し、各展開係数の横軸を位相回転速度(角周波数)に対応させることによって、運動反射物のエコー信号と静止対象物の低周波成分とが区別される。つまり、静止反射物と運動反射物との反射信号を区別することができる。
本発明によれば、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減できる。
N回のパルス送信によって反射した受信信号を並べた図である。 離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。 離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。 離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。 離散的フーリエ展開を用いた速度分析による速度スペクトル係数の例を示した図である。 従来の低域遮断型MTIフィルタと多項式フィット・フィルタのもつ通過振幅特性の典型例を示す図である。 従来の低域遮断型MTIフィルタと多項式フィット・フィルタをそれぞれ通した後、離散的フーリエ展開を用いた速度分析を行った場合に得られる速度スペクトル係数の例を示す図である。 一連の8点からなる信号に関する離散的ルジャンドル関数の偶数次を示す図である。 一連の8点からなる信号に関する離散的ルジャンドル関数の奇数次を示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の自乗和比である。 ルジャンドル係数の展開係数の自乗和比より求めた速度校正曲線である。 本実施形態の速度計測装置を用いた超音波診断装置の構成を示すブロック図である。 本実施形態のドップラ速度検出アルゴリズムの例を示す図である。 本実施形態による速度検出結果の例を示す図である。 従来方式による速度検出結果の例を示す図である。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。 本実施形態による速度検出結果の例を示す図である。 従来方式による速度検出結果の例を示す図。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。 本実施形態による速度検出結果の例を示す図である。 従来方式による速度検出結果の例を示す図である。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。 本実施形態のドップラ速度検出方式による検出結果の一の例を示す図である。 本実施形態のドップラ速度検出方式による検出結果の他の例を示す図である。 従来方式による速度検出結果の例を示す図である。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。
(フィルタの原理)
各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べたN点の離散的時系列信号から速度分析を行う場合を考える。このような、パルス・ドップラ速度分析用として典型的な時系列複素信号は、受信されたエコー信号に90°位相の異なる2つの搬送周波数信号を乗じて得られる1組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる(図1参照)。
このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、n= 1,2,・・・, N とならんだ時系列信号を、周波数を表す指数k=−N/2,・・,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個となり、これらは、直交関数系を形成している。
反射物が1つしか存在せず、それがナイキスト限界範囲内の一定速度で運動している場合を想定し、運動速度の関数として、前記展開により得られる離散的フーリエ展開係数を計算し、N=8の場合について、その絶対値を図2にプロットした。横軸は、運動速度に対応する位相回転速度を、パルス送信の時間間隔の逆数 Pulse Repetition Frequency (PRF) の2π倍を単位として表示した。図2では展開係数を表す8個の関数が重なって見づらいので、そのうち4個ずつを見分けやすい線により図3及び図4に示した。
送信時間間隔のN倍の時間で、位相が整数回回転するような運動速度では、ゼロでない展開係数はただ1つであり、他の展開係数はすべてゼロとなり、きれいに速度分析されている。ところが、それ以外の位相回転が整数回とならない一般の運動速度では、それぞれの展開係数が、横軸上ピークからはるかに離れた運動速度においても、ピークの1/10程度の大きさの絶対値をもっている。これは、速度分析計として−20dB程度の遠隔的クロストークをもつことを意味する。検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号がドリフトするとき、これはクラッタ信号と運動反射物によるエコー信号とを区別できないという重大な問題を生じさせる。
この問題は、N個の時系列信号に、ハニング関数のような緩やかに立ち上がって、下がる窓関数による重みをかけることにより、かなり抑圧できるが、本質的には解決しない。また、この方法は、時系列信号の数が大幅に目減りしたのと同様の、望ましくない効果を生じ、特に、時系列信号の数が元々少ない場合には適さない方法である。
具体例として、ナイキスト限界の2/3の一定速度で運動する検出・分析対象である運動反射物があり、その1000倍のエコー信号振幅をもつクラッタ信号が、ナイキスト限界速度の1/100でドリフトしている場合を考える。図5には、このときの前記速度分析計の出力である速度スペクトルを示す。検出・分析対象である運動反射物のスペクトルは、速度ゼロ付近にピークをもつクラッタ信号からのクロストーク成分に覆いかくされてしまい、この出力スペクトルの中では検出不可能である。このことが、まさに、クラッタ信号を遮断するMTIフィルタが必要な最も本質的な理由である。
次に、MTI処理として最もよく知られ、時間領域において畳み込み積和で表現される通常の低域遮断フィルタについて、その通過振幅特性の典型例を、図6中に実線で示す。このようなフィルタにn1 = 1,・・・, (N+M−1) とならんだ時系列信号S1(n1) を、入力すると、出力としてn0=1,・・・, N とならんだ時系列信号S0(n0)が得られる。
S0(n0)=ΣF(m)S1(n0+m−1) (3)
ここで、Σは m=1,・・・,Mの和を表す。図6では、低域遮断フィルタとして最も簡単なM=3,F(1)=−1, F(2)=2, F(3)=−1の場合について実線で示した。図5と同様の入力信号を、この低域遮断フィルタを通した後、図5と同様に速度分析すると、図7中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が1/2000倍程度に抑圧された結果、ナイキスト限界の0.3倍から0.4倍の速度にピークをもつ運動反射物のスペクトルが見えるようになっている。
同じ目的のMTIフィルタとして、前記の多項式フィット・フィルタを設計したときの通過振幅特性の典型例を、図6中に破線で示す。この例では、0次から3次式までをフィットして引き去った。前記の畳み込み型フィルタよりも急峻な低域遮断特性が見てとれる。また、図5と同様の入力信号を、この低域遮断フィルタを通した後、図5と同様に速度分析した結果を、図7中に破線で示した。低域遮断特性の急峻さのために、クラッタ信号が完全に抑圧されている。多項式フィット・フィルタ特有の急峻な低域遮断特性は以下のように理解することができる。
連続座標系における狭義ルジャンドル多項式については、岩波数学公式III特殊関数82−85頁 に記述されており、例えば、8次まで記すと
0(x)=1
1(x)=x=cosθ
2(x)=(1/2)(3x2−1)=(1/4)(3cos2θ+1)
3(x)=(1/2)(5x3−3x)=(1/8)(5cos3θ+3cosθ)
4(x)=(1/8)(35x4−30x2+3)=(1/64)(35cos4θ+20cos2θ+9)
5(x)=(1/8)(63x5−70x3+15x)=(1/128)(63cos5θ+35cos3θ+30cosθ)
6(x)=(1/16)(231x6−315x4+105x2−5)=(1/512)(231cos6θ+126cos4θ+105cos2θ+50)
7(x)=(1/16)(429x7−639x5+315x3−35x)=(1/210)(429cos7θ+231cos5θ+189cos3θ+175cosθ)
8(x)=(1/128)(6435x8−12012x6+6930x4−1260x2+35)=(1/214)(6435cos8θ+3432cos6θ+2772cos4θ+2520cos2θ+1225)
である。離散的座標系についても使えるように、一般化すると、偶数(2n)次ルジャンドル多項式は、2nより低い次数のルジャンドル多項式すべてと直交する最高次数2nの偶関数と定義することができ、奇数(2n+1)次ルジャンドル多項式は、(2n+1)より低い次数のルジャンドル多項式すべてと直交する最高次数(2n+1)の奇関数と定義することができる。
例として、一連の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(1),A(6)を破線、A(2),A(5)を1点鎖線、A(3),A(4)を点線にて示した。
図10から図12を、図2から図4と同様、入力複素正弦波に対する周波数応答として眺めてみる。ピークを主応答、それ以外を副応答とみる。ルジャンドル関数は、図2から図4のフーリエ展開系に比べ、主応答より高周波側の副応答は大きいが、低周波側には副応答をもたず、次数が大きいほど周波数原点付近で急速にゼロに収束する。ここでは解析的証明は省略するが、周波数の関数としてのm次ルジャンドル関数は、周波数原点において0次から(m−1)次導関数まですべてが0となる。
前記の多項式フィット・フィルタは、入力信号を、この離散的ルジャンドル関数を基底として展開し、0からM次までを除き、(M+1)次以上の成分の和を出力信号とするフィルタであるということができる。したがって、出力信号は、周波数原点付近で周波数の(M+1)乗に比例する挙動でゼロに収束していく。図6及び図7中に破線で示された例のように、多項式フィット・フィルタの低域遮断特性が優れているのは、このようなルジャンドル関数の本質的な特性に由来するものである。このように考えると、クラッタという周波数原点付近に桁違いに大きな振幅をもつ妨害成分の影響を抑えて、着目する運動物体を検出したり、その速度を分析したりしなければならないドップラ速度検出には、多項式フィット・フィルタやルジャンドル関数展開を利用する方式は、本質的に適しているといえる。
フーリエ展開系では、複素正弦波を基底とすることにより、図2から図5や図7に示したように符号つきの速度分析が可能であることは広く知られており、そのような分析方法が広く用いられている。ところが、ルジャンドル展開系は、前記のように、ドップラ速度検出に極めて適した性質をもつものの、このまま用いたのでは、図10から図20より見てとれるように、絶対値が等しいが符号の異なる速度に対して等しい応答を示すので、符号つきの速度分析は不可能である。
そこで、次数が1つだけ異なる偶数次と奇数次のルジャンドル関数を虚数単位 j を係数として線型結合して、複素ルジャンドル関数を形成し、これにより符号つきの速度分析を可能とする。これは、偶数次のルジャンドル関数が余弦波に、奇数次のルジャンドル関数が正弦波に似た変化を示すことに着目し、余弦関数と正弦関数を虚数単位 j を係数として線型結合することにより、位相角をj倍した指数関数、すなわち複素正弦関数が得られ、この位相角の増減から符号つき速度が得られるのと同様に、複素ルジャンドル関数から符号つき速度を得ようとするものである。
具体的には、一連のN点からなる時系列信号を、まず、0次から (N−1)次の離散的ルジャンドル関数を基底として展開して、A(0)からA(N−1)までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数
C(±(2n−3/2))=A(2n−2)±jA(2n−1) (1≦2n−1≦N-1)
C(±(2n−1/2))=A(2n)±j A(2n−1) (2≦2n≦N−1;n:自然数, 複号同順) (4)
を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。
図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の範囲で構築できることが確認されている。
図13から図15に示した複素ルジャンドル係数の応答は、ルジャンドル展開係数の特徴である位相回転速度原点すなわちドップラ周波数原点のまわりの極めて優れたカットオフ特性を示している反面、図2から図4に示したフーリエ展開係数の応答と比較して、周波数サイドローブレベルが高めである。すなわち、周波数サイドローブレベルは、フーリエ展開係数の図2から図4の例では、0.2程度であるのに対し、複素ルジャンドル係数の図14の例では、0.6近くに達している。周波数サイドローブレベルが高いことは、S/N比すなわちノイズレベルに対する信号レベルが低い場合に、速度検出誤差を生む可能性が高いことを意味するので、周波数サイドローブレベルは抑圧することが望ましい。
そこで、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−1))=α(2n−1)×A(2n−2) +β(2n−1)×A(2n)±j A(2n−1) (1≦2n−1≦N−1)
C(±2n)=A(2n)±j [α(2n)×A(2n−1)+ β(2n)×A(2n+1) ] (2≦2n≦N−1;n:自然数,複号同順;α及びβ:正の実数 ) (5)
となる。
図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程度に抑えられる。
図21には、mを自然数とするとき、(m−1)次から(N−1)次までのルジャンドル係数の自乗和とm次から(N−1)次までのルジャンドル係数の自乗和との比
R(m)=[A(m)2+・・・+A(N−1)2] / [A(m−1)2+A(m)2+・・・+A(N−1)2] (1≦m≦N−1) (6)
を、前記の例と同様N=8の場合について示した。これらの曲線は、ドップラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップルする。その単調増加部分を用いることにより、図22に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドップラ速度を算出すればよい。
検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号が存在し、それがドリフトするときには、そのドリフトの程度が大きくなるに応じて、C(±1)、C(±2)、……と、次数の小さな係数を順次除いた複素ルジャンドル係数系を用いて、速度検出・分析を行えばよい。このようにルジャンドル係数の遮断次数を制御することにより、クラッタ信号の影響を抑圧しながら、運動反射物の速度検出・分析を行うことができる。
(第1実施形態)
以下、本発明の実施形態について説明する。
図23は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子1を構成する各素子は、切り替えスイッチ群2を介して、送波ビームフォーマ3と受波ビームフォーマ10に接続されている。送波ビームフォーマ3では、送受信シークエンス制御部6による制御に従って、送信波形メモリ5から送信波形選択部4により選択されて読み出された波形を用い、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子1の各素子により超音波パルスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子1に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ10では、送受信シークエンス制御部6による制御に従って、指向性を持つ受信感度を生成すべく、各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークエンス制御部6による制御に従って、受波メモリ選択部11により選択された受波信号メモリ12中の1つのバンクへ一旦書き込まれ、ドップラ信号分析すべきN個の時系列信号がそろったのちに読み出されて、位相回転検出器13,15、ダウン・ミキシング器14、血流信号検出・分析部16において速度を検出・分析すべく信号処理される。
読み出されたN個の時系列受信信号は、各々を得るために行われたパルス送信の時刻を基準とする時相の等しい部分を、送信時刻の順番に並べたのちに処理される。この速度検出アルゴリズムについて、図24のフローチャート及び図23を参照して説明する。まず、位相回転検出器13では、クラッタ信号を含んだままの時系列信号から位相回転を求める。位相回転を検出する最も典型的な方法はn=1,・・・, Nとならんだ隣接する時系列複素信号 S(n),S(n+1)をもとに、
P(n)=S(n+1)S(n)* /||S(n+1)||/||S(n)|| (n=1,・・・,N−1) (7)
を計算し、P(n)の平均値Paの位相から平均位相回転速度を求める。
ここで、S(n)*はS(n)の複素共役、||S(n)|| はS(n)の絶対値である。一般に、クラッタ・エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、Paは、概ねクラッタ・エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される(図24のS1)。また、時系列複素信号S(n)は、受信されたエコー信号に90°位相の異なる2つの搬送周波数信号を乗じて得られる。
ダウン・ミキシング器14では、ここで求めた平均値Pa又はPaの周辺領域における空間平均値を用い、エコー信号の位相回転平均値がゼロとなるようにミキシング処理が行われる。いいかえれば、ダウン・ミキシングによりクラッタ平均速度による位相回転が打ち消される(図24のS2参照)。すなわち、S(n)とPaより
Sd(n)=S(n)Pa* n (n=1,・・・, N) (8)
を得る。このダウン・ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。
得られた隣接サンプルである時系列複素信号Sd(n),Sd(n+1)をもとに、位相回転検出器15では、位相回転検出器13と同様に、
Pd(n)=Sd(n+1)Sd(n)* / ||Sd(n+1)|| /||Sd(n)|| (n=1,・・・, N−1) (9)
を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される(図24のS3)。この位相回転最大値の大きさに応じて、血流信号検出・分析部16の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される(図24のS4)。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数Mがより高く設定され、クラッタ成分が効果的に抑圧される。
ここで、本実施形態の特徴部である血流信号検出・分析部16の動作について、もう少し詳しく述べる。入力された時系列信号 Sd(1) ,・・・, Sd(N) は、まず、図8及び図9のような0から(N−1)次の離散的ルジャンドル関数に展開され、その展開係数A(0) ,・・・, A(N−1) が求められる。この計算は、次のような行列演算により容易に行うことができる。n次ルジャンドル関数を行ベクトル L(n) で表し、これをn=0から(N−1)までならべてN×N行列LLを生成すると、Sd(1) ,・・・, Sd(N)を要素とする列ベクトルSdからA(0) ,・・・,A(N−1)を要素とする列ベクトルAを得るための行列は、
F=(L・tLL)-1LL (10)
と求めることができる。ここで、tLLはLLの転置行列、LL-1はLLの逆行列を表す。この行列を予め準備しておけば、行列演算
A=F・Sd (11)
によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数A(0),・・・, A(N−1)から、式(5)に従って、C(±1) ,・・・, C(±(N−1))の複素ルジャンドル係数が生成される(図24のS5)。いいかえれば、n次ルジャンドル関数を行ベクトルL(n)で表し、これをn=0から(N−1)まで列方向に並べたN×N行列LLを転置した転置行列tLLを生成する過程と、転置行列tLLにN×N行列LLを前から乗算したものの逆行列(LL・tLL)-1を算出する過程と、N×N行列LLに逆行列(LL・tLL)-1を前から乗算する過程とによって算出された行列FFを要素が時系列信号である列ベクトルSd(n)に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数(遮断次数M)が決定され(図24のS4)、C(±1) から C(±(M+1)) までの低次数複素ルジャンドル係数が棄却される。
次に、棄却されないで残された正次数の複素ルジャンドル係数の自乗和と負次数の複素ルジャンドル係数の自乗和とを比較して、速度符号が決定される(図24のS6参照)。具体的に、大きな方の次数の符号を速度の符号とする。さらに、速度符号次数の複素ルジャンドル係数のうち、絶対値最大の係数が検出される(図24のS7参照)。具体的に、正なら正の次数、負なら負の次数の複素ルジャンドル係数と、最大次数(N−1)のルジャンドル係数のうち、絶対値最大の係数の次数mが検出される。
そして、検出したルジャンドル係数に対応する校正曲線を用いて血流速度を算出する(図24のS8参照)。すなわち、次数mに対応する図22に示したような速度校正曲線を用いて、式(6)のように(m−1)次から(N−1)次までのルジャンドル係数の自乗和とm次から(N−1)次までのルジャンドル係数の自乗和との比からドップラ速度を求める。
得られた血流速度信号は、エコー振幅検出器17及びエコー振幅圧縮器18により得られる静止臓器からのエコー信号と共に、スキャンコンバータ19に入力される。スキャンコンバータでは、入力された複数の信号を適宜重畳して表示器20にて2次元ないし3次元表示すべく、信号の生成・制御を行う。
次に、本実施形態のドップラ速度計測装置の動作例を従来方式と比較して示す。クラッタ・エコー信号の振幅は、血流エコーの300倍とし、クラッタ速度については、初期値がゼロで一定の加速度で立ち上がる場合を想定した。時系列信号数N=8の場合について、数値計算シミュレーションを行って比較した。従来方式としては、MTIフィルタとして図6の場合の低域遮断フィルタを用い、その後段において、式(5)を用いて位相回転量平均値Paを得たのと同様の方法で、速度を算出する処理を血流信号検出・分析部16において行う例を選んだ。さらに、MTIフィルタを、やはり図5に例を示したような多項式フィット・フィルタに置き換えた場合とも、比較した。この場合、遮断次数Mは、本発明の方法と同様、位相回転検出器15の出力信号に応じて制御した。
図25から図27には、クラッタの最終到達速度がナイキスト限界速度の0.8%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、及び多項式フィット・フィルタ方式による速度検出結果を示した。横軸に入力としての血流速度、縦軸に出力として検出された速度をとり、有効速度検出範囲を実線、範囲外を点線によりプロットした。また、理想的場合を一点鎖線により示した。このように、クラッタ速度が低い場合には、いずれの方式も、ほぼ正確な速度分析結果を与えているが、本実施形態の方式による誤差が最も小さい。
図28から図30には、同様に、クラッタの最終到達速度がナイキスト限界速度の3%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。このクラッタ速度域では、低域遮断フィルタ方式による速度検出は、ほぼ破綻に至っている。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差は大きくないものの、血流速度がそれより小さい場合の誤差が大きく、何らかの対策を加えない限り、使うことができない。これに対し、本実施形態の方式は、有効速度検出範囲において正確な速度検出結果を与えるだけでなく、血流速度がそれより小さい場合、素直に検出速度ゼロを出力している。このように、有効速度検出範囲において正確な速度検出結果を出力し、かつ、血流速度がそれより小さい場合、素直に、検出速度ゼロを出力するか、あるいは速度検出不能であることを表示することは、本実施形態の速度計測装置の特徴であるということができる。この特徴を用いれば、着目した速度計測装置について、装置内部の構成を仔細に検査するに及ぶことなく、本実施形態の実施の有無を確認することができる。
図31から図33には、さらに、クラッタの最終到達速度がナイキスト限界速度の20%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。ここまでクラッタ速度が大きくなると、従来方式による速度検出は、完全に破綻している。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差が大きくなると同時に、血流速度がそれより小さい場合の誤差も著しく大きくなってしまっている。これに対し、本実施形態の方式は、有効速度検出範囲がやや狭くなるものの、血流速度がナイキスト限界速度の60%以上の範囲で、ほぼ正確な速度検出結果を与えるだけでなく、血流速度がそれより未満の場合の動作が、やはり素直である。
図34から図37には、血流速度自身がナイキスト限界速度の8%に対応するランダムな位相誤差をもつ場合について、それぞれ、式(5)を用いる本実施形態の方式、式(4)を用いる方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。クラッタの最終到達速度については、ナイキスト限界速度の10%とした。血流速度自身が位相誤差をもつために、負側の低中速領域において、本実施形態の方式も破綻しかかっているが、式(4)を用いる方式では、速度の符号が誤って検出されているのに対し、式(5)を用いる方式では、速度符号については、誤ることなく検出されている。これは、周波数サイドローブレベルを抑圧するよう構築された式(5)を用いる方式の効果である。多項式フィット・フィルタ方式も、同じ負側の低中速領域において、速度検出が破綻しており、その他の速度域においても、かろうじて速度符号が検出できたにとどまっている。また、クラッタ速度が大きいため、従来方式による速度検出は完全に破綻している。
以上述べたように本実施形態によれば、血流を、臓器の運動と峻別して描画することができる。この特徴をさらに生かして、図23の超音波診断装置では、位相回転検出器13により検出されたクラッタ信号すなわち臓器エコーの位相回転速度の信号をスキャンコンバータ19に入力し、臓器の運動速度あるいはその空間微分の分布を示す像を血流像と重畳あるいは並べて表示可能としている。この構成の有用性を、肝腫瘍の場合を例として述べる。肝腫瘍の辺縁部は新生血管が発達し、血流動態が周辺の正常肝とは異なる。また、周辺の正常肝とは固さが異なるため組織の運動も周辺とは異なる。したがって、血流像に加えて臓器の運動速度の場所による変化を表示することは、肝腫瘍の診断にきわめて有用な画像を提供することになる。
また、生体内部の血流など運動する反射物の速度を、そのエコー信号に比べ桁違いに振幅が大きいクラッタ・エコー信号の影響を除いて、正確に検出することができる。具体的には、臓器の超音波探触子への方向の運動速度が、1mm/s変化するなか、超音波探触子への方向の速度成分 3mm/s以上の血流をリアルタイムに描画できる。このように、本実施形態によれば、的確な医療診断を可能とする血流検出・描画機能つき超音波診断装置を提供することが可能となる。すなわち、本実施形態を実施した装置の医用診断上の有用性はきわめて大きく、従って、医用診断を支える工業における本実施形態の意義も、また、大きい。さらに、本実施形態の方法は、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置など、他のパルス・ドップラ装置の運動反射体検出能力をも飛躍的に向上させることができ、本実施形態の工業的ならびに社会的意義は、その意味でもきわめて大である。
(変形例)
本実施形態は前記した実施形態に限定されるものではなく、例えば以下のような種々の変形が可能である。
(1)前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数mを求め、次数mを境界とした展開係数間の自乗和の比を用いてドップラ速度を算出したが、複素展開係数間の自乗和の比を用いることもできる。
(2)前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置に用いることができる。
本発明は、パルス状の反射波を用いた速度計測方法及びこれを用いた速度計測装置に関する。
パルス・ドップラ速度計では、複数のパルス波を運動反射物に送信し、図1のような受信された複数のエコー信号が解析される。すなわち、各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べた時系列信号が解析され、運動反射物に関する速度などの情報が得られる。最も一般的な信号処理方法としては、これらの受信信号を直交検波し、複素時系列信号として解析する。このようなパルス・ドップラ速度計は、生体内部の血流などを検出・描画する超音波診断装置、雨雲などを検出・描画する気象レーダー、飛行物体を検出する航空レーダーなどとして広く用いられている。
もし、反射物が1つしか存在しなければ、前記の時系列信号の位相回転速度 ΔΦ/Δtすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度vを、次式のように容易に求めることができる。
v=λΔΦ/Δt/2 (1)
ここで、λは波長、Δtはパルス送信の時間間隔、ΔΦは位相回転角である。運動速度vの符号は、反射物が近づく場合には、運動速度vの符号が正となり、遠ざかる場合には負となる。
ところで、前記のような実際の受信エコー中では、静止反射物によるエコー信号、いわゆるクラッタ信号の強度が、着目する運動反射物によるエコー信号の強度よりも、一般に、数桁以上大きい。また、静止反射物によるエコー信号は、実際には、時間軸上で完全に静止しているわけではなく、途中の媒質のゆらぎや、医療用超音波診断装置の場合には、血流以外の静止臓器の遅い動きのために、ドリフトしている。このため、静止反射物によるエコー信号の周波数成分は、位相回転速度ゼロの直流成分だけでなく、位相回転速度がゼロでない低周波成分を有している。
このため、前記時系列信号に単純に (1)式の処理を適用しただけでは、運動反射物の速度を検出することはできない。そこで、実際のパルス・ドップラ速度計は、静止反射物によるエコー信号を抑圧し、運動反射物によるエコー信号を相対的に強調するMTI(Moving Target Indicator) 処理とよばれる信号処理を行った後に、速度の検出あるいは分析を行うように構成されている。
MTI処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点がある。
1) N1個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出・分析処理部に入力されるデータ点が (N1−1) 個目減りする。
2) 急峻な遮断特性のフィルタを得にくい。
パルス・ドップラ速度計の場合は、運動反射物によるエコー信号を保ちながら、静止反射物のドリフト成分を除去する急峻な遮断特性が必要であり、特に、2)が問題となる。また、N回の送受信により得られるN点の時系列信号から速度検出・分析処理をするとき、データ点が目減りする1)の問題は、実際に速度検出・分析演算に用い得るデータ点が(N−N1+1) 点に減ってしまうことを意味する。これは、医療診断用超音波血流描画装置のように実時間性が重要な応用では、望ましいことではない。
この問題を解決する処理としては、非特許文献1に記載されているようにPolynomial Regression Filter (多項式フィット・フィルタ)が提案されている。これは、時系列信号に、0次式(定数)、1次式、2次式、・・・・、M次式を順次最小2乗フィットし、フィットした成分を引き去ることにより、元の時系列信号のもつドリフト成分を除去する処理である。この処理は、N点の入力時系列信号にN×Nの正方行列を掛けることで表現される。このため、出力信号として同じN点の時系列信号が得られ、
1) データ点の目減りがない。
また、その低域遮断特性は、フィットする最大次数Mの大きさによるが、
2) 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 42., 927-937頁 (1995年)
ところが、この多項式フィット・フィルタによる処理の後段において、v=λΔΦ/Δt/2式を用いて運動物体の速度の検出・分析を行うと、遮断周波数の近くで速度演算誤差が、低域遮断フィルタと比較しても大きくなってしまうという問題点がある。このため、多項式フィット・フィルタが有する急峻な遮断特性を十分に使うことができない。
そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。
本発明は、速度計測対象物に対してN個のパルス波を所定間隔で送信し、受信エコー信号を送信時刻の順番に並べたN個の時系列信号を0次から(N−1)次までの離散的ルジャンドル関数を基底として展開し、この展開係数を用いて速度を算出する。特に、1次以上(N−1)次以下の(2n−1)次の展開係数と(2n+1)次の展開係数との線形結合に虚数単位を乗じて、2n次の展開係数と線形結合した2n次の複素展開係数と、(2n+1)次の展開係数に虚数単位を乗じた後、2n次の展開係数と(2n+2)次の展開係数と線形結合した(2n+1)次の複素展開係数とを算出し、各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数mを決定し、次数mに対応した前記各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出する。
パルス波が照射された運動反射物からは位相回転速度(角周波数)ΔΦ/Δtを有したエコー信号が得られ、ドリフト振動を伴った静止反射物からは低周波成分を有したクラッタ信号が得られる。このクラッタ信号の強度は、エコー信号の強度よりも数桁大きい。エコー信号を次数が大きいほど周波数原点付近で急速にゼロに収束する特性を有する離散的ルジャンドル関数を基底として展開し、各展開係数の横軸を位相回転速度(角周波数)に対応させることによって、運動反射物のエコー信号と静止対象物の低周波成分とが区別される。つまり、静止反射物と運動反射物との反射信号を区別することができる。
本発明によれば、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減できる。
(フィルタの原理)
各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べたN点の離散的時系列信号から速度分析を行う場合を考える。このような、パルス・ドップラ速度分析用として典型的な時系列複素信号は、受信されたエコー信号に90°位相の異なる2つの搬送周波数信号を乗じて得られる1組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる(図1参照)。
このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、n= 1,2,・・・, N とならんだ時系列信号を、周波数を表す指数k=−N/2,・・,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個となり、これらは、直交関数系を形成している。
反射物が1つしか存在せず、それがナイキスト限界範囲内の一定速度で運動している場合を想定し、運動速度の関数として、前記展開により得られる離散的フーリエ展開係数を計算し、N=8の場合について、その絶対値を図2にプロットした。横軸は、運動速度に対応する位相回転速度を、パルス送信の時間間隔の逆数 Pulse Repetition Frequency (PRF) の2π倍を単位として表示した。図2では展開係数を表す8個の関数が重なって見づらいので、そのうち4個ずつを見分けやすい線により図3及び図4に示した。
送信時間間隔のN倍の時間で、位相が整数回回転するような運動速度では、ゼロでない展開係数はただ1つであり、他の展開係数はすべてゼロとなり、きれいに速度分析されている。ところが、それ以外の位相回転が整数回とならない一般の運動速度では、それぞれの展開係数が、横軸上ピークからはるかに離れた運動速度においても、ピークの1/10程度の大きさの絶対値をもっている。これは、速度分析計として−20dB程度の遠隔的クロストークをもつことを意味する。検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号がドリフトするとき、これはクラッタ信号と運動反射物によるエコー信号とを区別できないという重大な問題を生じさせる。
この問題は、N個の時系列信号に、ハニング関数のような緩やかに立ち上がって、下がる窓関数による重みをかけることにより、かなり抑圧できるが、本質的には解決しない。また、この方法は、時系列信号の数が大幅に目減りしたのと同様の、望ましくない効果を生じ、特に、時系列信号の数が元々少ない場合には適さない方法である。
具体例として、ナイキスト限界の2/3の一定速度で運動する検出・分析対象である運動反射物があり、その1000倍のエコー信号振幅をもつクラッタ信号が、ナイキスト限界速度の1/100でドリフトしている場合を考える。図5には、このときの前記速度分析計の出力である速度スペクトルを示す。検出・分析対象である運動反射物のスペクトルは、速度ゼロ付近にピークをもつクラッタ信号からのクロストーク成分に覆いかくされてしまい、この出力スペクトルの中では検出不可能である。このことが、まさに、クラッタ信号を遮断するMTIフィルタが必要な最も本質的な理由である。
次に、MTI処理として最もよく知られ、時間領域において畳み込み積和で表現される通常の低域遮断フィルタについて、その通過振幅特性の典型例を、図6中に実線で示す。このようなフィルタにn1 = 1,・・・, (N+M−1) とならんだ時系列信号S1(n1) を、入力すると、出力としてn0=1,・・・, N とならんだ時系列信号S0(n0)が得られる。
S0(n0)=ΣF(m)S1(n0+m−1) (3)
ここで、Σは m=1,・・・,Mの和を表す。図6では、低域遮断フィルタとして最も簡単なM=3,F(1)=−1, F(2)=2, F(3)=−1の場合について実線で示した。図5と同様の入力信号を、この低域遮断フィルタを通した後、図5と同様に速度分析すると、図7中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が1/2000倍程度に抑圧された結果、ナイキスト限界の0.3倍から0.4倍の速度にピークをもつ運動反射物のスペクトルが見えるようになっている。
じ目的のMTIフィルタとして、前記の多項式フィット・フィルタを設計したときの通過振幅特性の典型例を、図6中に破線で示す。この例では、0次から3次式までをフィットして引き去った。前記の畳み込み型フィルタよりも急峻な低域遮断特性が見てとれる。また、図5と同様の入力信号を、この多項式フィット・フィルタを通した後、図5と同様に速度分析した結果を、図7中に破線で示した。低域遮断特性の急峻さのために、クラッタ信号が完全に抑圧されている。多項式フィット・フィルタ特有の急峻な低域遮断特性は以下のように理解することができる。
連続座標系における狭義ルジャンドル多項式については、岩波数学公式III特殊関数82−85頁 に記述されており、例えば、8次まで記すと
0(x)=1
1(x)=x=cosθ
2(x)=(1/2)(3x2−1)=(1/4)(3cos2θ+1)
3(x)=(1/2)(5x3−3x)=(1/8)(5cos3θ+3cosθ)
4(x)=(1/8)(35x4−30x2+3)=(1/64)(35cos4θ+20cos2θ+9)
5(x)=(1/8)(63x5−70x3+15x)=(1/128)(63cos5θ+35cos3θ+30cosθ)
6(x)=(1/16)(231x6−315x4+105x2−5)=(1/512)(231cos6θ+126cos4θ+105cos2θ+50)
7(x)=(1/16)(429x7−639x5+315x3−35x)=(1/210)(429cos7θ+231cos5θ+189cos3θ+175cosθ)
8(x)=(1/128)(6435x8−12012x6+6930x4−1260x2+35)=(1/214)(6435cos8θ+3432cos6θ+2772cos4θ+2520cos2θ+1225)
である。離散的座標系についても使えるように、一般化すると、偶数(2n)次ルジャンドル多項式は、2nより低い次数のルジャンドル多項式すべてと直交する最高次数2nの偶関数と定義することができ、奇数(2n+1)次ルジャンドル多項式は、(2n+1)より低い次数のルジャンドル多項式すべてと直交する最高次数(2n+1)の奇関数と定義することができる。
例として、一連の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(1),A(6)を破線、A(2),A(5)を1点鎖線、A(3),A(4)を点線にて示した。
図10から図12を、図2から図4と同様、入力複素正弦波に対する周波数応答として眺めてみる。ピークを主応答、それ以外を副応答とみる。ルジャンドル関数は、図2から図4のフーリエ展開系に比べ、主応答より高周波側の副応答は大きいが、低周波側には副応答をもたず、次数が大きいほど周波数原点付近で急速にゼロに収束する。ここでは解析的証明は省略するが、周波数の関数としてのm次ルジャンドル関数は、周波数原点において0次から(m−1)次導関数まですべてが0となる。
前記の多項式フィット・フィルタは、入力信号を、この離散的ルジャンドル関数を基底として展開し、0からM次までを除き、(M+1)次以上の成分の和を出力信号とするフィルタであるということができる。したがって、出力信号は、周波数原点付近で周波数の(M+1)乗に比例する挙動でゼロに収束していく。図6及び図7中に破線で示された例のように、多項式フィット・フィルタの低域遮断特性が優れているのは、このようなルジャンドル関数の本質的な特性に由来するものである。このように考えると、クラッタという周波数原点付近に桁違いに大きな振幅をもつ妨害成分の影響を抑えて、着目する運動物体を検出したり、その速度を分析したりしなければならないドップラ速度検出には、多項式フィット・フィルタやルジャンドル関数展開を利用する方式は、本質的に適しているといえる。
フーリエ展開系では、複素正弦波を基底とすることにより、図2から図5や図7に示したように符号つきの速度分析が可能であることは広く知られており、そのような分析方法が広く用いられている。ところが、ルジャンドル展開系は、前記のように、ドップラ速度検出に極めて適した性質をもつものの、このまま用いたのでは、図10から図20より見てとれるように、絶対値が等しいが符号の異なる速度に対して等しい応答を示すので、符号つきの速度分析は不可能である。
そこで、次数が1つだけ異なる偶数次と奇数次のルジャンドル関数を虚数単位 j を係数として線型結合して、複素ルジャンドル関数を形成し、これにより符号つきの速度分析を可能とする。これは、偶数次のルジャンドル関数が余弦波に、奇数次のルジャンドル関数が正弦波に似た変化を示すことに着目し、余弦関数と正弦関数を虚数単位 j を係数として線型結合することにより、位相角をj倍した指数関数、すなわち複素正弦関数が得られ、この位相角の増減から符号つき速度が得られるのと同様に、複素ルジャンドル関数から符号つき速度を得ようとするものである。
具体的には、一連のN点からなる時系列信号を、まず、0次から (N−1)次の離散的ルジャンドル関数を基底として展開して、A(0)からA(N−1)までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数
C(±(2n−3/2))=A(2n−2)±jA(2n−1) (1≦2n−1≦N-1)
C(±(2n−1/2))=A(2n)±j A(2n−1) (2≦2n≦N−1;n:自然数, 複号同順) (4)
を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。
図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の範囲で構築できることが確認されている。
図13から図15に示した複素ルジャンドル係数の応答は、ルジャンドル展開係数の特徴である位相回転速度原点すなわちドップラ周波数原点のまわりの極めて優れたカットオフ特性を示している反面、図2から図4に示したフーリエ展開係数の応答と比較して、周波数サイドローブレベルが高めである。すなわち、周波数サイドローブレベルは、フーリエ展開係数の図2から図4の例では、0.2程度であるのに対し、複素ルジャンドル係数の図14の例では、0.6近くに達している。周波数サイドローブレベルが高いことは、/N比すなわちノイズレベルに対する信号レベルが低い場合に、速度検出誤差を生む可能性が高いことを意味するので、周波数サイドローブレベルは抑圧することが望ましい。
そこで、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−1))=α(2n−1)×A(2n−2) +β(2n−1)×A(2n)±j A(2n−1) (1≦2n−1≦N−1)
C(±2n)=A(2n)±j [α(2n)×A(2n−1)+ β(2n)×A(2n+1) ] (2≦2n≦N−1;n:自然数,複号同順;α及びβ:正の実数 ) (5)
となる。
図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程度に抑えられる。
図21には、mを自然数とするとき、(m−1)次から(N−1)次までのルジャンドル係数の自乗和とm次から(N−1)次までのルジャンドル係数の自乗和との比
R(m)=[A(m)2+・・・+A(N−1)2] / [A(m−1)2+A(m)2+・・・+A(N−1)2] (1≦m≦N−1) (6)
を、前記の例と同様N=8の場合について示した。これらの曲線は、ドップラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップルする。その単調増加部分を用いることにより、図22に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドップラ速度を算出すればよい。
検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号が存在し、それがドリフトするときには、そのドリフトの程度が大きくなるに応じて、C(±1)、C(±2)、……と、次数の小さな係数を順次除いた複素ルジャンドル係数系を用いて、速度検出・分析を行えばよい。このようにルジャンドル係数の遮断次数を制御することにより、クラッタ信号の影響を抑圧しながら、運動反射物の速度検出・分析を行うことができる。
(第1実施形態)
以下、本発明の実施形態について説明する。
図23は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子1を構成する各素子は、切り替えスイッチ群2を介して、送波ビームフォーマ3と受波ビームフォーマ10に接続されている。送波ビームフォーマ3では、送受信シークエンス制御部6による制御に従って、送信波形メモリ5から送信波形選択部4により選択されて読み出された波形を用い、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子1の各素子により超音波パルスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子1に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ10では、送受信シークエンス制御部6による制御に従って、指向性を持つ受信感度を生成すべく、各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークエンス制御部6による制御に従って、受波メモリ選択部11により選択された受波信号メモリ12中の1つのバンクへ一旦書き込まれ、ドップラ信号分析すべきN個の時系列信号がそろったのちに読み出されて、位相回転検出器13,15、ダウン・ミキシング器14、血流信号検出・分析部16において速度を検出・分析すべく信号処理される。
読み出されたN個の時系列受信信号は、各々を得るために行われたパルス送信の時刻を基準とする時相の等しい部分を、送信時刻の順番に並べたのちに処理される。この速度検出アルゴリズムについて、図24のフローチャート及び図23を参照して説明する。まず、位相回転検出器13では、クラッタ信号を含んだままの時系列信号から位相回転を求める。位相回転を検出する最も典型的な方法はn=1,・・・, Nとならんだ隣接する時系列複素信号 S(n),S(n+1)をもとに、
P(n)=S(n+1)S(n)* /||S(n+1)||/||S(n)|| (n=1,・・・,N−1) (7)
を計算し、P(n)の平均値Paの位相から平均位相回転速度を求める。
ここで、S(n)*はS(n)の複素共役、||S(n)|| はS(n)の絶対値である。一般に、クラッタ・エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、Paは、概ねクラッタ・エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される(図24のS1)。また、時系列複素信号S(n)は、受信されたエコー信号に90°位相の異なる2つの搬送周波数信号を乗じて得られる。
ダウン・ミキシング器14では、ここで求めた平均値Pa又はPaの周辺領域における空間平均値を用い、エコー信号の位相回転平均値がゼロとなるようにミキシング処理が行われる。いいかえれば、ダウン・ミキシングによりクラッタ平均速度による位相回転が打ち消される(図24のS2参照)。すなわち、S(n)とPaより
Sd(n)=S(n)Pa* n (n=1,・・・, N) (8)
を得る。このダウン・ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。
得られた隣接サンプルである時系列複素信号Sd(n),Sd(n+1)をもとに、位相回転検出器15では、位相回転検出器13と同様に、
Pd(n)=Sd(n+1)Sd(n)* / ||Sd(n+1)|| /||Sd(n)|| (n=1,・・・, N−1) (9)
を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される(図24のS3)。この位相回転最大値の大きさに応じて、血流信号検出・分析部16の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される(図24のS4)。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数Mがより高く設定され、クラッタ成分が効果的に抑圧される。
ここで、本実施形態の特徴部である血流信号検出・分析部16の動作について、もう少し詳しく述べる。入力された時系列信号 Sd(1) ,・・・, Sd(N) は、まず、図8及び図9のような0から(N−1)次の離散的ルジャンドル関数に展開され、その展開係数A(0) ,・・・, A(N−1) が求められる。この計算は、次のような行列演算により容易に行うことができる。n次ルジャンドル関数を行ベクトル L(n) で表し、これをn=0から(N−1)までならべてN×N行列LLを生成すると、Sd(1) ,・・・, Sd(N)を要素とする列ベクトルSdからA(0) ,・・・,A(N−1)を要素とする列ベクトルAを得るための行列は、
F=(L・tLL)-1LL (10)
と求めることができる。ここで、tLLはLLの転置行列、LL-1はLLの逆行列を表す。この行列を予め準備しておけば、行列演算
A=F・Sd (11)
によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数A(0),・・・, A(N−1)から、式(5)に従って、C(±1) ,・・・, C(±(N−1))の複素ルジャンドル係数が生成される(図24のS5)。いいかえれば、n次ルジャンドル関数を行ベクトルL(n)で表し、これをn=0から(N−1)まで列方向に並べたN×N行列LLを転置した転置行列 LLを生成する過程と、転置行列 LLにN×N行列LLを前から乗算したものの逆行列(LL・ LL)-1を算出する過程と、N×N行列LLに逆行列(LL・ LL)-1を前から乗算する過程とによって算出された行列FFを要素が時系列信号である列ベクトルSd(n)に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数(遮断次数M)が決定され(図24のS4)、C(±1) から C(±(M+1)) までの低次数複素ルジャンドル係数が棄却される。
次に、棄却されないで残された正次数の複素ルジャンドル係数の自乗和と負次数の複素ルジャンドル係数の自乗和とを比較して、速度符号が決定される(図24のS6参照)。具体的に、大きな方の次数の符号を速度の符号とする。さらに、速度符号次数の複素ルジャンドル係数のうち、絶対値最大の係数が検出される(図24のS7参照)。具体的に、正なら正の次数、負なら負の次数の複素ルジャンドル係数と、最大次数(N−1)のルジャンドル係数のうち、絶対値最大の係数の次数mが検出される。
そして、検出したルジャンドル係数に対応する校正曲線を用いて血流速度を算出する(図24のS8参照)。すなわち、次数mに対応する図22に示したような速度校正曲線を用いて、式(6)のように(m−1)次から(N−1)次までのルジャンドル係数の自乗和とm次から(N−1)次までのルジャンドル係数の自乗和との比からドップラ速度を求める。
得られた血流速度信号は、エコー振幅検出器17及びエコー振幅圧縮器18により得られる静止臓器からのエコー信号と共に、スキャンコンバータ19に入力される。スキャンコンバータでは、入力された複数の信号を適宜重畳して表示器20にて2次元ないし3次元表示すべく、信号の生成・制御を行う。
次に、本実施形態のドップラ速度計測装置の動作例を従来方式と比較して示す。クラッタ・エコー信号の振幅は、血流エコーの300倍とし、クラッタ速度については、初期値がゼロで一定の加速度で立ち上がる場合を想定した。時系列信号数N=8の場合について、数値計算シミュレーションを行って比較した。従来方式としては、MTIフィルタとして図6の場合の低域遮断フィルタを用い、その後段において、式(5)を用いて位相回転量平均値Paを得たのと同様の方法で、速度を算出する処理を血流信号検出・分析部16において行う例を選んだ。さらに、MTIフィルタを、やはり図5に例を示したような多項式フィット・フィルタに置き換えた場合とも、比較した。この場合、遮断次数Mは、本発明の方法と同様、位相回転検出器15の出力信号に応じて制御した。
図25から図27には、クラッタの最終到達速度がナイキスト限界速度の0.8%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、及び多項式フィット・フィルタ方式による速度検出結果を示した。横軸に入力としての血流速度、縦軸に出力として検出された速度をとり、有効速度検出範囲を実線、範囲外を点線によりプロットした。また、理想的場合を一点鎖線により示した。このように、クラッタ速度が低い場合には、いずれの方式も、ほぼ正確な速度分析結果を与えているが、本実施形態の方式による誤差が最も小さい。
図28から図30には、同様に、クラッタの最終到達速度がナイキスト限界速度の3%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。このクラッタ速度域では、低域遮断フィルタ方式による速度検出は、ほぼ破綻に至っている。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差は大きくないものの、血流速度がそれより小さい場合の誤差が大きく、何らかの対策を加えない限り、使うことができない。これに対し、本実施形態の方式は、有効速度検出範囲において正確な速度検出結果を与えるだけでなく、血流速度がそれより小さい場合、素直に検出速度ゼロを出力している。このように、有効速度検出範囲において正確な速度検出結果を出力し、かつ、血流速度がそれより小さい場合、素直に、検出速度ゼロを出力するか、あるいは速度検出不能であることを表示することは、本実施形態の速度計測装置の特徴であるということができる。この特徴を用いれば、着目した速度計測装置について、装置内部の構成を仔細に検査するに及ぶことなく、本実施形態の実施の有無を確認することができる。
図31から図33には、さらに、クラッタの最終到達速度がナイキスト限界速度の20%である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。ここまでクラッタ速度が大きくなると、従来方式による速度検出は、完全に破綻している。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差が大きくなると同時に、血流速度がそれより小さい場合の誤差も著しく大きくなってしまっている。これに対し、本実施形態の方式は、有効速度検出範囲がやや狭くなるものの、血流速度がナイキスト限界速度の60%以上の範囲で、ほぼ正確な速度検出結果を与えるだけでなく、血流速度がそれより未満の場合の動作が、やはり素直である。
図34から図37には、血流速度自身がナイキスト限界速度の8%に対応するランダムな位相誤差をもつ場合について、それぞれ、式(5)を用いる本実施形態の方式、式(4)を用いる方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。クラッタの最終到達速度については、ナイキスト限界速度の10%とした。血流速度自身が位相誤差をもつために、負側の低中速領域において、本実施形態の方式も破綻しかかっているが、式(4)を用いる方式では、速度の符号が誤って検出されているのに対し、式(5)を用いる方式では、速度符号については、誤ることなく検出されている。これは、周波数サイドローブレベルを抑圧するよう構築された式(5)を用いる方式の効果である。多項式フィット・フィルタ方式も、同じ負側の低中速領域において、速度検出が破綻しており、その他の速度域においても、かろうじて速度符号が検出できたにとどまっている。また、クラッタ速度が大きいため、従来方式による速度検出は完全に破綻している。
以上述べたように本実施形態によれば、血流を、臓器の運動と峻別して描画することができる。この特徴をさらに生かして、図23の超音波診断装置では、位相回転検出器13により検出されたクラッタ信号すなわち臓器エコーの位相回転速度の信号をスキャンコンバータ19に入力し、臓器の運動速度あるいはその空間微分の分布を示す像を血流像と重畳あるいは並べて表示可能としている。この構成の有用性を、肝腫瘍の場合を例として述べる。肝腫瘍の辺縁部は新生血管が発達し、血流動態が周辺の正常肝とは異なる。また、周辺の正常肝とは固さが異なるため組織の運動も周辺とは異なる。したがって、血流像に加えて臓器の運動速度の場所による変化を表示することは、肝腫瘍の診断にきわめて有用な画像を提供することになる。
また、生体内部の血流など運動する反射物の速度を、そのエコー信号に比べ桁違いに振幅が大きいクラッタ・エコー信号の影響を除いて、正確に検出することができる。具体的には、臓器の超音波探触子への方向の運動速度が、1mm/s変化するなか、超音波探触子への方向の速度成分 3mm/s以上の血流をリアルタイムに描画できる。このように、本実施形態によれば、的確な医療診断を可能とする血流検出・描画機能つき超音波診断装置を提供することが可能となる。すなわち、本実施形態を実施した装置の医用診断上の有用性はきわめて大きく、従って、医用診断を支える工業における本実施形態の意義も、また、大きい。さらに、本実施形態の方法は、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置など、他のパルス・ドップラ装置の運動反射体検出能力をも飛躍的に向上させることができ、本実施形態の工業的ならびに社会的意義は、その意味でもきわめて大である。
(変形例)
本実施形態は前記した実施形態に限定されるものではなく、例えば以下のような種々の変形が可能である。
(1)前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数mを求め、次数mを境界とした展開係数間の自乗和の比を用いてドップラ速度を算出したが、複素展開係数間の自乗和の比を用いることもできる。
(2)前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置に用いることができる。
N回のパルス送信によって反射した受信信号を並べた図である。 離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。 離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。 離散的フーリエ展開係数の絶対値を、速度に対応する符号つき周波数の関数として示した図である。 離散的フーリエ展開を用いた速度分析による速度スペクトル係数の例を示した図である。 従来の低域遮断型MTIフィルタと多項式フィット・フィルタのもつ通過振幅特性の典型例を示す図である。 従来の低域遮断型MTIフィルタと多項式フィット・フィルタをそれぞれ通した後、離散的フーリエ展開を用いた速度分析を行った場合に得られる速度スペクトル係数の例を示す図である。 一連の8点からなる信号に関する離散的ルジャンドル関数の偶数次を示す図である。 一連の8点からなる信号に関する離散的ルジャンドル関数の奇数次を示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 複素ルジャンドル係数の規格化絶対値を、速度に対応する符号つき周波数の関数として示す図である。 ルジャンドル展開係数の自乗和比である。 ルジャンドル係数の展開係数の自乗和比より求めた速度校正曲線である。 本実施形態の速度計測装置を用いた超音波診断装置の構成を示すブロック図である。 本実施形態のドップラ速度検出アルゴリズムの例を示す図である。 本実施形態による速度検出結果の例を示す図である。 従来方式による速度検出結果の例を示す図である。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。 本実施形態による速度検出結果の例を示す図である。 従来方式による速度検出結果の例を示す図。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。 本実施形態による速度検出結果の例を示す図である。 従来方式による速度検出結果の例を示す図である。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。 本実施形態のドップラ速度検出方式による検出結果の一の例を示す図である。 本実施形態のドップラ速度検出方式による検出結果の他の例を示す図である。 従来方式による速度検出結果の例を示す図である。 多項式フィット・フィルタ方式による速度検出結果の例を示す図である。

Claims (9)

  1. 速度計測対象物に対してN個のパルス波を所定間隔で送信し、前記パルス波の受信エコー信号を送信時刻の順番に並べた時系列信号を用いて前記速度計測対象物を構成する運動反射物の速度を計測する速度計測方法であって、
    前記N個の時系列信号を0次から(N−1)次までの離散的ルジャンドル関数を基底として展開する展開ステップと、
    自然数nとするとき、前記展開ステップで展開した1次以上(N−1)次以下の(2n−1)次の展開係数と(2n+1)次の展開係数との線形結合に虚数単位を乗じて、2n次の展開係数と線形結合した2n次の複素展開係数と、(2n+1)次の展開係数に虚数単位を乗じた後、2n次の展開係数と(2n+2)次の展開係数と線形結合した(2n+1)次の複素展開係数とを算出する複素展開係数算出ステップと、
    前記各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数mを決定する次数決定ステップと、
    前記次数mに対応した前記各展開係数間あるいは前記各複素展開係数間の自乗和の比から前記運動反射物に関する符号付き速度信号を算出する速度信号算出ステップと
    からなる速度計測方法。
  2. 前記展開ステップは、n次ルジャンドル関数を行ベクトルL(n)で表し、これをn=0から(N−1)まで列方向に並べたN×N行列LLを転置した転置行列tLLを生成する過程と、前記転置行列tLLに前記N×N行列LLを前から乗算したものの逆行列(LL・tLL)-1を算出する過程と、前記N×N行列LLに前記逆行列(LL・tLL)-1を前から乗算する過程とによって算出された行列FFを要素が前記時系列信号である列ベクトルSd(n)に前から乗算するステップであることを特徴とする請求の範囲第1項に記載の速度計測方法。
  3. 静止反射物に対する前記受信エコーの大きさの変動あるいはドリフトの程度によって、前記複素展開係数を小さな次数から徐々に無視することを特徴とする請求の範囲第2項に記載の速度計測方法。
  4. 前記運動反射物の速度が有効速度検出範囲よりも小さい場合に検出速度ゼロを出力するか否かによって請求の範囲第3項に記載の速度計測方法の実行を判定するステップを備えたことを特徴とする速度計測方法。
  5. 速度計測対象物に対してN個のパルス波を所定間隔で送信し、前記パルス波の受信エコー信号を送信時刻の順番に並べた時系列信号を用いて前記速度計測対象物を構成する運動反射物の速度を計測する速度計測装置であって、
    前記N個の時系列信号を0次から(N−1)次までの離散的ルジャンドル関数を基底として展開する展開手段と、
    自然数nとするとき、前記展開ステップで展開した1次以上(N−1)次以下の(2n−1)次の展開係数と(2n+1)次の展開係数との線形結合に虚数単位を乗じて、2n次の展開係数と線形結合した2n次の複素展開係数と、(2n+1)次の展開係数に虚数単位を乗じた後、2n次の展開係数と(2n+2)次の展開係数と線形結合した(2n+1)次の複素展開係数とを算出する複素展開係数算出手段と、
    前記各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数mを決定する次数決定手段と、
    前記次数mに対応した前記各展開係数間あるいは前記各複素展開係数間の自乗和の比から前記運動反射物に関する符号付き速度信号を算出する速度信号算出手段と
    からなる速度計測装置。
  6. 前記展開手段は、n次ルジャンドル関数を行ベクトルL(n)で表し、これをn=0から(N−1)まで列方向に並べたN×N行列LLを転置した転置行列tLLを生成する過程と、前記転置行列tLLに前記N×N行列LLを前から乗算したものの逆行列(LL・tLL)-1を算出する過程と、前記N×N行列LLに前記逆行列(LL・tLL)-1を前から乗算する過程とによって算出された行列FFを要素が前記時系列信号である列ベクトルSd(n)に前から乗算する手段であることを特徴とする請求の範囲第5項に記載の速度計測装置。
  7. 静止反射物に対する前記受信エコーの大きさの変動あるいはドリフトの程度によって、前記複素展開係数を小さな次数から徐々に無視することを特徴とする請求の範囲第6項に記載の速度計測装置。
  8. 前記パルス波は、超音波であり、
    前記運動反射物は、血流であることを特徴とする請求の範囲第5項ないし第7項に記載の速度計測装置。
  9. 前記パルス波は、超音波であり、
    前記速度計測対象物は、肝臓であり、
    前記肝臓の速度の空間微分の分布及び血流像を表示する表示手段を備えたことを特徴とする請求の範囲第5項ないし7項に記載の速度計測装置。
JP2007517722A 2005-05-27 2006-01-06 速度計測方法及びこれを用いた速度計測装置 Expired - Fee Related JP5108512B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2007517722A JP5108512B2 (ja) 2005-05-27 2006-01-06 速度計測方法及びこれを用いた速度計測装置

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2005154844 2005-05-27
JP2005154844 2005-05-27
JP2007517722A JP5108512B2 (ja) 2005-05-27 2006-01-06 速度計測方法及びこれを用いた速度計測装置
PCT/JP2006/300081 WO2006126303A1 (ja) 2005-05-27 2006-01-06 速度計測方法及びこれを用いた速度計測装置

Publications (2)

Publication Number Publication Date
JPWO2006126303A1 true JPWO2006126303A1 (ja) 2008-12-25
JP5108512B2 JP5108512B2 (ja) 2012-12-26

Family

ID=37451736

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2007517722A Expired - Fee Related JP5108512B2 (ja) 2005-05-27 2006-01-06 速度計測方法及びこれを用いた速度計測装置

Country Status (5)

Country Link
US (1) US7946992B2 (ja)
EP (1) EP1884196A4 (ja)
JP (1) JP5108512B2 (ja)
CN (1) CN101184443B (ja)
WO (1) WO2006126303A1 (ja)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8585733B2 (en) 2006-04-19 2013-11-19 Vibrynt, Inc Devices, tools and methods for performing minimally invasive abdominal surgical procedures
JP5872323B2 (ja) * 2011-03-29 2016-03-01 株式会社東芝 X線ct装置及び画像処理方法
US9636076B2 (en) 2011-03-29 2017-05-02 Toshiba Medical Systems Corporation X-ray CT apparatus and image processing method
WO2013056160A2 (en) 2011-10-13 2013-04-18 Masimo Corporation Medical monitoring hub
CN102503138A (zh) * 2011-10-18 2012-06-20 武汉理工大学 一种烧结法建筑装饰微晶玻璃的着色方法
US20130172734A1 (en) * 2011-12-30 2013-07-04 General Electric Company Flow measurement with time-resolved data
JP6253360B2 (ja) * 2013-02-13 2017-12-27 キヤノン株式会社 被検体情報取得装置、被検体情報取得方法、及びプログラム
TWI507709B (zh) * 2013-03-20 2015-11-11 Univ Nat Taiwan 藉由格雷編碼激發之超音波都卜勒偵測方法
KR20160046669A (ko) * 2014-10-21 2016-04-29 알피니언메디칼시스템 주식회사 빔포밍 장치, 초음파 이미징 장치 및 빔포밍 방법
US10342509B2 (en) * 2015-03-27 2019-07-09 Alpinion Medical Systems Co., Ltd. Beamforming device, ultrasonic imaging device, and beamforming method allowing simple spatial smoothing operation
DE102015207894A1 (de) * 2015-04-29 2016-11-03 Bayer Pharma Aktiengesellschaft Ermitteln der Geschwindigkeit eines Fluids mit Hilfe eines bildgebenden Verfahrens
WO2019216086A1 (ja) * 2018-05-09 2019-11-14 古野電気株式会社 気象レーダ装置、気象観測方法、および、気象観測プログラム

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS6480351A (en) * 1987-09-24 1989-03-27 Hitachi Medical Corp Ultrasonic doppler meter
JPH0484953A (ja) * 1990-07-30 1992-03-18 Matsushita Electric Ind Co Ltd 超音波ドプラ映像装置
JP2003521341A (ja) * 2000-01-31 2003-07-15 アー.ヤー. アンゲルセン、ビョルン 医療用超音波イメージングにおける位相面収差およびパルス残響の補正

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1003768B (zh) * 1985-04-01 1989-04-05 复旦大学 一种双向定量测量血流绝对速度的超声多普勒方法与仪器
US5228009A (en) * 1992-04-10 1993-07-13 Diasonics, Inc. Parametric clutter elimination

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS6480351A (en) * 1987-09-24 1989-03-27 Hitachi Medical Corp Ultrasonic doppler meter
JPH0484953A (ja) * 1990-07-30 1992-03-18 Matsushita Electric Ind Co Ltd 超音波ドプラ映像装置
JP2003521341A (ja) * 2000-01-31 2003-07-15 アー.ヤー. アンゲルセン、ビョルン 医療用超音波イメージングにおける位相面収差およびパルス残響の補正

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
CSNC200901611046, 複素ルジャンドル法によるクラッタに強いドプラ速度検出, "一般演題抄録(口演)その1", 超音波医学 Vol.32, Supplement JAPANESE JOURNAL OF MEDICAL ULTRASONICS, 2005, 第32巻, p.281, JP, (社)日本超音波医学会 *
JPN6009038695, Hans Torp, "Clutter Rejection Filters in Color Flow Imaging: A Theoretical Approach", IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, 199703, vol.44, no.2, pp.417−424 *

Also Published As

Publication number Publication date
EP1884196A4 (en) 2009-12-23
US7946992B2 (en) 2011-05-24
US20090177091A1 (en) 2009-07-09
CN101184443A (zh) 2008-05-21
JP5108512B2 (ja) 2012-12-26
CN101184443B (zh) 2011-12-28
WO2006126303A1 (ja) 2006-11-30
EP1884196A1 (en) 2008-02-06

Similar Documents

Publication Publication Date Title
JP5108512B2 (ja) 速度計測方法及びこれを用いた速度計測装置
Ebrahimkhanlou et al. Damage localization in metallic plate structures using edge-reflected lamb waves
Giusti et al. ISAR image resolution enhancement: Compressive sensing versus state-of-the-art super-resolution techniques
Aderson Multi-dimensional velocity estimation with ultrasound using spatial quadrature
Stankovic et al. Inverse radon transform–based micro-Doppler analysis from a reduced set of observations
Zhao et al. Mode identification and extraction of broadband ultrasonic guided waves
JP4369427B2 (ja) ドプラ速度検出装置及びそれを用いた超音波診断装置
Jia et al. Rigid and elastic acoustic scattering signal separation for underwater target
US11719613B2 (en) Method and device for quantifying viscoelasticity of a medium
Noumura et al. Feasibility of low-frequency directive sound source with high range resolution using pulse compression technique
Mozumi et al. Anti-aliasing method for ultrasonic 2D phase-sensitive motion estimator
JP2007215816A (ja) パルスドプラ計測装置、その方法及びそのプログラム
Mozumi et al. Impact of spacing of ultrasound receiving beams on estimation of 2D motion velocity
CN107907591B (zh) 多组分固液两相混合物组分浓度的超声检测系统和方法
Azouz et al. Design and implementation of an enhanced matched filter for sidelobe reduction of pulsed linear frequency modulation radar
Sando et al. High-frequency ultrasonic airborne Doppler method for noncontact elasticity measurements of living tissues
Tian et al. Parametric modeling in food package defect imaging
Madhavanunni et al. Triangulation based vector flow imaging with non-steered plane waves for transverse flows
JPH05146438A (ja) 超音波装置およびそれを用いたベクトルドプラ法
Xiaoming et al. Application of Hilbert–Huang transform to laser Doppler velocimeter
Liu et al. Order analysis for measurement of Lamb wave propagation distance
Latif et al. Dispersion analysis of acoustic circumferential waves using time-frequency representations
Nam et al. Polar coordinate-based guided wave beamforming imaging using scanning LDV
Kuga et al. Alias-free estimation of blood velocity using phase difference and amplitude correlation of chirp echo spectrum
Li Enhanced Ultrasonic Imaging and Non-Destructive Evaluation Using Angle-of-Arrival Estimation in Dimension-Reduced Beam Space

Legal Events

Date Code Title Description
A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20101109

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20110111

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20110308

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20110509

A02 Decision of refusal

Free format text: JAPANESE INTERMEDIATE CODE: A02

Effective date: 20120313

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20120611

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20120615

A911 Transfer to examiner for re-examination before appeal (zenchi)

Free format text: JAPANESE INTERMEDIATE CODE: A911

Effective date: 20120705

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: 20121002

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20121005

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: 20151012

Year of fee payment: 3

LAPS Cancellation because of no payment of annual fees