1992年[ 技術開発研究助成 ] 成果報告 ： 年報第06号
心臓疾患の音響的精密診断のための心音の計測技術・時系列分析法の開発に関する研究
概要
Ⅰ.まえがき
我々はすでに,短い信号間のスペクトルの遷移パターンを,低いSN比の下でも高精度に推定できる新しいスペクトル推定法を提案した[1,2]。心臓の動きに伴って出力される心音の一拍分の中で,心房から心室へ血液が流れ込む直後に出力されるIV音は,心室壁の振動が含まれるため,壁心筋の弾性的特性に関する情報があると考えられるが,振幅が小さくSN比が悪い上に継続時間も短いため,通常のスペクトル解析を適用できない。しかし,上記の手法を,STRESSTESTで得られた複数のIV音に適用した結果,心筋梗塞の患者と正常者十数名において,得られたスペクトル遷移パターンに差異が認められた。こうした結果は,弾性的・力学的特徴に基づく心臓疾患の音響的診断の可能性を示すものである。しかし,これらの実験においては,心音は食道内の加速度センサで検出しており,患者に対する苦痛を無視できない。また,心臓各部分の振動の積分の結果として得られる心音に対して,心臓の局所的な振動の検出と詳細な解析によって,新たな知見が得られる可能性もある。
以上のように,超音波ドプラによって,胸壁上から心臓の各部位の振動を計測するための新しい診断システムの開発が望まれる。しかし,心臓壁の振動は,数Hzで振幅±8～15mmの拍動の上に重畳しているため,反射波の戻って来るまでの時間は,数μs～20μsの範囲で刻々と変化する。従って,従来の診断装置で用いられているようなアナログ回路によるゲートやサンプルホールドでは,タイミングがずれてしまい,大振幅の拍動にのっているある1点の微小振動を,一拍全体にわたって計測することができないという問題がある。
そこで,我々は,直交検波回路の出力を,標本化周波数1MHz程度で高速にA/D変換し,得られたディジタル信号に対して,反射波の戻って来るまでの時間を各パルスごとに求めるなどの処理を詳細に行って,大振幅で動いている対象の微小振動の計測を行う手法を提案し[3],計算機シミュレーションと水槽を用いた評価実験,大動脈弁近傍の大動脈壁の振動の計測を行った[4,5]。
II.内容
1,振動速度計測の原理
図1のように,超音波トランスジューサからの送信出力である,角周波数嫡,パルス幅ToのRFパルス波形Vin(t)を次式で表わす。
u(t>は単位ステップ関数,ん。は振幅を表わす。時刻tにおいて,対象物体がトランスジューサとの距離x(t)の位置をトランスジューサから遠ざかる方向を正方向として速度v(t)で振動しているとする。このとき,時刻tにおける往路/復路の各々の片道の伝搬時間ztは,音速をCとすると,Zt=x(t)/oで表される。往路と復路の伝播による位相遅れは,ドプラ角周波数偏移を△ω、とおくと,各々e-/moz,,e一ノ(ω・+△ω・)τであるから,再びトランスジューサへ戻ってきたときの波形Vout(t)は次式で表される。
ここで,A。utは振幅を表わす。参照波Vref(t)=e一酬をこのVout(t>に掛け合わせることによって得られる波形の低周波成分(直交検波出力)Vm(t)は,次式で表される。
もし,このVm(t)が超音波周波数f0=ω0/2πに近い周波数でA/D変換される場合には,その結果得られるディジタル信号から,パルスの戻ってきた時刻の同定が可能となる。微小振動が大振幅の拍動上にのっている場合にも,各々のパルス位置の検出を行うことができる。その反射パルスの検出された時点でのVm(t)の位相をBtと表わす。
さらに瀾係式△ωt≒-2ω0(t)/cとx(t)-/v(t)∫v(t)dtを用いることによって,位相Btは,次のように整理できる。
通常の超音波ドプラでは,トランスジューサと対象問の距離x(t)が,一定である場合のみを扱っている。そのため,△T後に送信される次のパルスの反射波に関して位相e.△Tを求め,θ之の差をとれば,(5)式の第1項に相当する部分は相殺される。第2項のτ濁一定であり,従って,△ωこのみを考えればよいことになる。
しかし,心臓における大振幅の拍動上の微小振動検出に関しては,第1項の変化を無視できない。しかも,通常△ωfは,ω。に対して十分小さいため,第2項が無視できる。従って,位相Btは,次のように表される。
この位相Btと,次の反射波の出力の位相Br.△rとの差は,次式で表される。
ここで,右辺の伝搬時間の差(Zt+△,一τ∂は,対象までの距離の差x(t+△t)-x(t>と音速Cで割った量である。さらに,時刻t+△tとtの中間の時刻における,対象の速度v(t+△t/2)を,{x(t+△t)-x(t)}/△tで定義する。これらから,対象の速度v(t)は,次のように推定できる。
すなわち,パルスの受信間隔△死その間での直交検波出力の位相の差(Br.△,一θ∂によって,対象の振動速度を推定することができる。
変位関数をx(t)=aX。sin(2nfX。t)とおくと,速度の最大値vmaxは2παx。fX。で,△T間での位置の変化の最大値△xmax=vm。X△T,△T間での伝搬時間の差△τ=Zt+△T一銑の最大値△Zm。Xは2vmax△T/Cとなる。従って(5)式右辺の中で支配的な項であるω。による位相変化の最大値1△Blm。xは2ωo△Tmaxで表される。
例えば,fX。=1.5Hz,aX。=±3.5mmのとき,「△BlmaX=2.76rad=158degと大きな値となる。これは検波出力には,対象位置の変位による伝搬時間の変化にf(=3.5MHz)が掛け合わさって,検波出力の位相変化がおこることを意味しており,対象の位置の変化の取扱いが非常に重要であることがわかる。
また,本手法による推定限界は,パルス送信繰り返し周期△Tの中での検波出力の位相差1θ,.△T-Btuがπ未満であることが最低必要となる。従って(8)式より,推定できる瞬時速度の最大値は,次の不等式を満足する必要がある。
右辺の値は,△T=500μs,f=3.5MHz,o=1500m/sのときの値である。
2.シミュレーション実験による精度評価
以上の処理によって推定される速度の精度の評価を行っために,シミュレーション実験を行った。推定された系列v(t>と,真の速度波形v(t)との差の二乗和をv(t)のパワーで正規化して得た値(正規化誤差パワー)ρ(ax。,fX。)を次式で計算する。
その際の変位関数には,x(t)=aX。sin(2砿。t)を用い,その振幅aX。と周波fX。の2つのパラメータを変化させ,各々の場合に関して,正規化誤差パワーによって評価した結果を図2に示す。ここでは,(3)式のVmfit)を計算機内部で12bit整数により表現した場合に,推定された振動速度に関する正規化誤差パワーを示す。図2(a)(b)は各々,用いるA/D変換器が標本化周波数100kHz,1MHzであるとして,Vm(t)を10μs,1μ∫間隔で生成した場合に対応する。
これらいずれの場合にも,(9)式のvm。x=0.107[m/s]を,推定の上限としている。また,図2(b)に対して,図2(a)のVm(t)が10μs間隔でA/D変換された場合には,ごく狭い条件範囲でしか,推定できないことがわかる。空間分解能を上げる点以外に,A/D変換の精度が有限であることによって,位相θご一上に発生する測定誤差を図2(b)のように改善できる点からも,高い周波数でA/D変換することの有用性がわかる。
3.水槽実験
実際に低周波大振幅の運動に,100Hz帯域での微小振動を重畳させられる水槽を作製した。図3はその概観を示す。モータと偏心カムによって,低周波大振幅の変位(振幅±3.5mm,1.5Hzのほぼ正弦的変位)を発生させる。その偏心カム上の加振器によって,微小振動を発生させ,
水槽内のトランスジューサを用いて,それら2成分の和として運動するゴム板の振動速度を推定した。実験の条件としては,RF周波数f=3.5[MHz],パルス間隔△T=500[μs]である。
図4は,振幅±80μm,30～100Hz帯域の白色雑音で,加振器を微小振動させた場合の推定結果を示している。図4(c)は,反射パルスの戻って来るまでの時間2×rt,図4(b)は,本手法による速度の推定値v(t)とそのスペクトルである。本手法で得られた波形とそのパワースペクトルには,偏心カムによる低周波成分も検出されている。
図(b-1)の波形の後半部分にカットオフ周波数30HzのHIGHPASSFILTERをかけて得られた波形を図(c)に示す。また,図(d)と(e)には各々,レーザ変位計による変位測定値を差分して得た速度波形vd(t),加振器の入力信号の差分波形び,n(t)を示す。図(d)のレーザ変位計による変位測定の際には,偏心カムによる低周波変動成分をHIGHPASSFILTERで除いている。
また,図4(c)(d)の2つの信号間のコヒーレンス関数γ2(f)を計算し,その結果を図4(f)に示す。図4(c)(d)(e)の波形が互いにほぼ対応がとれ,図4(f)のコヒーレンス関数も30～100Hzの帯域でほぼ1であることから,大変位に重畳する振動の計測に,本方法が有効であることが確認された。
4.大動脈弁付近の大動脈壁の振動の検出例
市販の超音波診断装置によって,心臓の断面映像を調べ,トランスジューサの方向を1方向に固定し,その方向に関して,RFパルス送信・受信して得られた信号に関して,上述の本手法を適用した結果を図5に示す。20代前半の健康な男性の大動脈弁付近の大動脈壁からの反射パルスに関して処理した例である。図のように,弁の開閉に伴って発生する振動が計測されている。なお,超音波ビームは,この中隔面にほぼ垂直になるように実験を行っている。また,パルス送信間隔△Tは,250[μs]に設定している。
III.成果
心臓壁の振動計測を想定し,低周波大振幅変位上の微小振動の計測のための新しい手法を提案し,計算機シミュレーションと実際の水槽実験から手法の評価を行って,心臓の振動を検出する際に必要な精度が得られることを確認した。
IV.まとめ
心音の分析を通じて出てきた問題を解決するために,心筋の局所的振動計測法を開発し,精度評価と人体への適応実験を行った。今後,計測システムとして整備するとともに,心臓の各部位の振動の計測とその結果得られた振動波形に論文[1,2]の手法を適用して心疾患の診断法の開発を行う予定である。
我々はすでに,短い信号間のスペクトルの遷移パターンを,低いSN比の下でも高精度に推定できる新しいスペクトル推定法を提案した[1,2]。心臓の動きに伴って出力される心音の一拍分の中で,心房から心室へ血液が流れ込む直後に出力されるIV音は,心室壁の振動が含まれるため,壁心筋の弾性的特性に関する情報があると考えられるが,振幅が小さくSN比が悪い上に継続時間も短いため,通常のスペクトル解析を適用できない。しかし,上記の手法を,STRESSTESTで得られた複数のIV音に適用した結果,心筋梗塞の患者と正常者十数名において,得られたスペクトル遷移パターンに差異が認められた。こうした結果は,弾性的・力学的特徴に基づく心臓疾患の音響的診断の可能性を示すものである。しかし,これらの実験においては,心音は食道内の加速度センサで検出しており,患者に対する苦痛を無視できない。また,心臓各部分の振動の積分の結果として得られる心音に対して,心臓の局所的な振動の検出と詳細な解析によって,新たな知見が得られる可能性もある。
以上のように,超音波ドプラによって,胸壁上から心臓の各部位の振動を計測するための新しい診断システムの開発が望まれる。しかし,心臓壁の振動は,数Hzで振幅±8～15mmの拍動の上に重畳しているため,反射波の戻って来るまでの時間は,数μs～20μsの範囲で刻々と変化する。従って,従来の診断装置で用いられているようなアナログ回路によるゲートやサンプルホールドでは,タイミングがずれてしまい,大振幅の拍動にのっているある1点の微小振動を,一拍全体にわたって計測することができないという問題がある。
そこで,我々は,直交検波回路の出力を,標本化周波数1MHz程度で高速にA/D変換し,得られたディジタル信号に対して,反射波の戻って来るまでの時間を各パルスごとに求めるなどの処理を詳細に行って,大振幅で動いている対象の微小振動の計測を行う手法を提案し[3],計算機シミュレーションと水槽を用いた評価実験,大動脈弁近傍の大動脈壁の振動の計測を行った[4,5]。
II.内容
1,振動速度計測の原理
図1のように,超音波トランスジューサからの送信出力である,角周波数嫡,パルス幅ToのRFパルス波形Vin(t)を次式で表わす。
u(t>は単位ステップ関数,ん。は振幅を表わす。時刻tにおいて,対象物体がトランスジューサとの距離x(t)の位置をトランスジューサから遠ざかる方向を正方向として速度v(t)で振動しているとする。このとき,時刻tにおける往路/復路の各々の片道の伝搬時間ztは,音速をCとすると,Zt=x(t)/oで表される。往路と復路の伝播による位相遅れは,ドプラ角周波数偏移を△ω、とおくと,各々e-/moz,,e一ノ(ω・+△ω・)τであるから,再びトランスジューサへ戻ってきたときの波形Vout(t)は次式で表される。
ここで,A。utは振幅を表わす。参照波Vref(t)=e一酬をこのVout(t>に掛け合わせることによって得られる波形の低周波成分(直交検波出力)Vm(t)は,次式で表される。
もし,このVm(t)が超音波周波数f0=ω0/2πに近い周波数でA/D変換される場合には,その結果得られるディジタル信号から,パルスの戻ってきた時刻の同定が可能となる。微小振動が大振幅の拍動上にのっている場合にも,各々のパルス位置の検出を行うことができる。その反射パルスの検出された時点でのVm(t)の位相をBtと表わす。
さらに瀾係式△ωt≒-2ω0(t)/cとx(t)-/v(t)∫v(t)dtを用いることによって,位相Btは,次のように整理できる。
通常の超音波ドプラでは,トランスジューサと対象問の距離x(t)が,一定である場合のみを扱っている。そのため,△T後に送信される次のパルスの反射波に関して位相e.△Tを求め,θ之の差をとれば,(5)式の第1項に相当する部分は相殺される。第2項のτ濁一定であり,従って,△ωこのみを考えればよいことになる。
しかし,心臓における大振幅の拍動上の微小振動検出に関しては,第1項の変化を無視できない。しかも,通常△ωfは,ω。に対して十分小さいため,第2項が無視できる。従って,位相Btは,次のように表される。
この位相Btと,次の反射波の出力の位相Br.△rとの差は,次式で表される。
ここで,右辺の伝搬時間の差(Zt+△,一τ∂は,対象までの距離の差x(t+△t)-x(t>と音速Cで割った量である。さらに,時刻t+△tとtの中間の時刻における,対象の速度v(t+△t/2)を,{x(t+△t)-x(t)}/△tで定義する。これらから,対象の速度v(t)は,次のように推定できる。
すなわち,パルスの受信間隔△死その間での直交検波出力の位相の差(Br.△,一θ∂によって,対象の振動速度を推定することができる。
変位関数をx(t)=aX。sin(2nfX。t)とおくと,速度の最大値vmaxは2παx。fX。で,△T間での位置の変化の最大値△xmax=vm。X△T,△T間での伝搬時間の差△τ=Zt+△T一銑の最大値△Zm。Xは2vmax△T/Cとなる。従って(5)式右辺の中で支配的な項であるω。による位相変化の最大値1△Blm。xは2ωo△Tmaxで表される。
例えば,fX。=1.5Hz,aX。=±3.5mmのとき,「△BlmaX=2.76rad=158degと大きな値となる。これは検波出力には,対象位置の変位による伝搬時間の変化にf(=3.5MHz)が掛け合わさって,検波出力の位相変化がおこることを意味しており,対象の位置の変化の取扱いが非常に重要であることがわかる。
また,本手法による推定限界は,パルス送信繰り返し周期△Tの中での検波出力の位相差1θ,.△T-Btuがπ未満であることが最低必要となる。従って(8)式より,推定できる瞬時速度の最大値は,次の不等式を満足する必要がある。
右辺の値は,△T=500μs,f=3.5MHz,o=1500m/sのときの値である。
2.シミュレーション実験による精度評価
以上の処理によって推定される速度の精度の評価を行っために,シミュレーション実験を行った。推定された系列v(t>と,真の速度波形v(t)との差の二乗和をv(t)のパワーで正規化して得た値(正規化誤差パワー)ρ(ax。,fX。)を次式で計算する。
その際の変位関数には,x(t)=aX。sin(2砿。t)を用い,その振幅aX。と周波fX。の2つのパラメータを変化させ,各々の場合に関して,正規化誤差パワーによって評価した結果を図2に示す。ここでは,(3)式のVmfit)を計算機内部で12bit整数により表現した場合に,推定された振動速度に関する正規化誤差パワーを示す。図2(a)(b)は各々,用いるA/D変換器が標本化周波数100kHz,1MHzであるとして,Vm(t)を10μs,1μ∫間隔で生成した場合に対応する。
これらいずれの場合にも,(9)式のvm。x=0.107[m/s]を,推定の上限としている。また,図2(b)に対して,図2(a)のVm(t)が10μs間隔でA/D変換された場合には,ごく狭い条件範囲でしか,推定できないことがわかる。空間分解能を上げる点以外に,A/D変換の精度が有限であることによって,位相θご一上に発生する測定誤差を図2(b)のように改善できる点からも,高い周波数でA/D変換することの有用性がわかる。
3.水槽実験
実際に低周波大振幅の運動に,100Hz帯域での微小振動を重畳させられる水槽を作製した。図3はその概観を示す。モータと偏心カムによって,低周波大振幅の変位(振幅±3.5mm,1.5Hzのほぼ正弦的変位)を発生させる。その偏心カム上の加振器によって,微小振動を発生させ,
水槽内のトランスジューサを用いて,それら2成分の和として運動するゴム板の振動速度を推定した。実験の条件としては,RF周波数f=3.5[MHz],パルス間隔△T=500[μs]である。
図4は,振幅±80μm,30～100Hz帯域の白色雑音で,加振器を微小振動させた場合の推定結果を示している。図4(c)は,反射パルスの戻って来るまでの時間2×rt,図4(b)は,本手法による速度の推定値v(t)とそのスペクトルである。本手法で得られた波形とそのパワースペクトルには,偏心カムによる低周波成分も検出されている。
図(b-1)の波形の後半部分にカットオフ周波数30HzのHIGHPASSFILTERをかけて得られた波形を図(c)に示す。また,図(d)と(e)には各々,レーザ変位計による変位測定値を差分して得た速度波形vd(t),加振器の入力信号の差分波形び,n(t)を示す。図(d)のレーザ変位計による変位測定の際には,偏心カムによる低周波変動成分をHIGHPASSFILTERで除いている。
また,図4(c)(d)の2つの信号間のコヒーレンス関数γ2(f)を計算し,その結果を図4(f)に示す。図4(c)(d)(e)の波形が互いにほぼ対応がとれ,図4(f)のコヒーレンス関数も30～100Hzの帯域でほぼ1であることから,大変位に重畳する振動の計測に,本方法が有効であることが確認された。
4.大動脈弁付近の大動脈壁の振動の検出例
市販の超音波診断装置によって,心臓の断面映像を調べ,トランスジューサの方向を1方向に固定し,その方向に関して,RFパルス送信・受信して得られた信号に関して,上述の本手法を適用した結果を図5に示す。20代前半の健康な男性の大動脈弁付近の大動脈壁からの反射パルスに関して処理した例である。図のように,弁の開閉に伴って発生する振動が計測されている。なお,超音波ビームは,この中隔面にほぼ垂直になるように実験を行っている。また,パルス送信間隔△Tは,250[μs]に設定している。
III.成果
心臓壁の振動計測を想定し,低周波大振幅変位上の微小振動の計測のための新しい手法を提案し,計算機シミュレーションと実際の水槽実験から手法の評価を行って,心臓の振動を検出する際に必要な精度が得られることを確認した。
IV.まとめ
心音の分析を通じて出てきた問題を解決するために,心筋の局所的振動計測法を開発し,精度評価と人体への適応実験を行った。今後,計測システムとして整備するとともに,心臓の各部位の振動の計測とその結果得られた振動波形に論文[1,2]の手法を適用して心疾患の診断法の開発を行う予定である。