超音波による動脈壁の弾性特性と組織性状の非侵襲イメージングに関する研究
長谷川英之

金井

浩

医工学研究科 医工学専攻
計測・診断医工学講座 准教授
工学研究科 電子工学専攻
電子制御工学講座 准教授
E-mail: hasegawa@us.ecei.tohoku.ac.jp

工学研究科 電子工学専攻
電子制御工学講座 教授
医工学研究科 医工学専攻
計測・診断医工学講座 教授
E-mail: hkanai@ecei.tohoku.ac.jp

１．はじめに近年，心筋梗塞や脳梗塞などの循環器系疾患が我が国でも大きな問題となっている．これら循
環器疾患の大きな原因の 1 つが動脈硬化症である
ことは広く知られており，重篤な循環器系疾患の発症を予防するためには，動脈硬化症の早期診断が重要である．動脈硬化症の進展により動脈壁の弾性特性は大きく変化するため [1,2]，脈波速度法やスティフネスパラメータ法など動脈壁の弾性特性の計測法について古くから研究開発が行われてきた [3-6]．これらの手法では，脈波伝搬速度を測定した 2 点間，もしくは血圧変化による動脈直径の変化を計測した部位の円周全体，それぞれの平均的な弾性特性を評価している．したがって，動脈壁の局所弾性特性，特に局所的に発生する動脈硬化病変の弾性特性の計測は困難である．動脈壁
局所の弾性特性を評価するために本研究グループ
では，心拍による動脈壁の径方向変位の径方向空間分布を計測することにより壁の径方向ひずみを推定し，血圧との関係から弾性特性を評価する手法を開発した [7-11]．この手法では，内圧変化による円筒管の変位は径方向のみであると仮定できることから，動脈壁の径方向のみ変位を，超音波
ビーム方向の変位として推定しているが，近年，動脈が，心臓の収縮や圧力波の伝搬などの影響によりその長軸方向 (超音波ビームと直交する方向:
ラテラル方向) にも変位していることが確認された[12]．したがって，動脈壁の動態をより高精度に計測するためには，高い時間分解能での計測を実現するともに，動脈壁をその長軸方向にもトラッキングする必要がある．本研究では，動脈壁の高速超音波イメージング方法について検討を行うとともに，長軸方向変位も含め，動脈壁の変位計測法に関する検討を行った [13,14]．

２．原理
２.１ 動脈の高速超音波イメージング法
図1 に示されるように，従来のリニア走査では送信・受信ともにビームを形成するため，超音
波送信繰り返し周波数を fPRF，ビーム本数をNbmとするとフレームレートは f PRF / N bm となる．

図1. Illustration of beam forming. Left: Conventional
linear scan. Right: Parallel beam forming.
一方，parallel beam forming (PBF) [15] では，広いビームを送信 (本報告では平面波) し，そのビーム内に複数の受信ビームを形成するため，リニ
ア走査と同じビーム本数を得るために必要な送信回数を減少させることができる．ただし，最終的なビームの指向性は，送信ビームと受信ビームの指向性の積となるため [16]，従来のリニア走査に比べ方位分解能は劣化する．
本報告では，N t = 96 チャネル同時送受信を行うことにより，parallel beam forming を試みた．1
つの受信ビームを形成するために使用する素子の数をN e とし，第 m ( m = 1, 2, ..., M ) 回目の送受信に使用する素子の番号を (N t − N e ) (m − 1) + i ( i= 1, 2, ..., N t) とすると，一回の送受信で(N t − N e)
本のビームを形成することができる．本報告では，
送信回数 M = 3，受信ビーム形成に使用する素子
数 N e = 72 とすることにより，M (N t − N e) = 72 本のビームを形成した．f PRF = 10416 Hz (観察可能深度: 7.4 cm) であるから，以上の処理によりフレー
ムレート f PRF / M = 約3500 Hz を実現できる．
２.２ 超音波ビーム方向の変位計測法
本研究では，動脈壁の径方向ひずみ ( = 壁厚の変化) の高精度計測を目的としている．頸動脈壁の厚みは約1 mm，径方向ひずみは数%であるため，壁厚の変化は数十ミクロンである．本研究において，動脈に対して送受信された超音波RF信号の標
本化周波数は40 MHzであり，この場合超音波ビー
ム方向の標本点間隔は約20 μmとなる．したがって，数十ミクロンの厚み変化を精度良く計測するためには，標本点間隔以下の変位を推定する必要があ
る．本研究では，このように粗く標本化されたデータを用いて，標本点の間隔以下のビーム方向変位を推定する手法について検討を行った．
２.２.１ 受信超音波の位相を用いた変位計測法
動脈壁に対して送受信され標本化された超音
波RF信号を用いて，標本化間隔よりも小さなビー
ム方向変位を高精度に推定するために，本研究グ
ループでは，受信超音波RF信号の位相を用いた変位推定法“位相差トラッキング法”を開発した [7]．
超音波診断装置の超音波プローブから送信された
超音波パルスは，微小振動している動脈壁で反射されて超音波プローブに戻るまでに，伝搬距離による位相遅れが生じる．著者らが開発した位相差
トラッキング法では，深さ (ビーム方向位置) z からの反射波の直交検波信号 g(z; t) = gI (z; t) + jgQ (z; t) から，連続する2つのフレームにおいて送受信された

超音波パルス間の位相差 Δθz(t) を複素相互相関関数
により検出する．
KejΔθˆ z ( t )=) g ( z + kΔZ ; t )∑ g ( z + kΔ Z ; t + Δ T ・
*

k =−K
K

(1)

) g ( z + kΔZ ; t )∑ g ( z + kΔ Z ; t + Δ T ・
*

k =−K

ここで，ΔT とΔZ は超音波パルスの送信繰り返し周期と標本点の深さ方向の間隔，K は複素相関関数を算出するために使用するビーム方向の標本点の数 (片側) である．また，*は複素共役を示す． (1)式により検出された位相差 をもとに，超音波ビー
ム方向の動脈壁の微小振動速度 vz(t) は，連続する2つのフレームの中間時点での値として，次のように求められる.
c Δθˆz (t )(2)vˆz (t ) = − 0
4π f 0 ΔT
ここで， f0は送信超音波の周波数，c0は音速である．
(2)式により推定された速度を時間積分することにより変位が得られる．以上のような処理を，各深さに適用することにより，超音波ビーム方向 (動脈径方向) の変位分布 が得られる．
２.２.２ Gradient-Based オプティカルフロー法
超音波RF信号の補間等を必要としない高速な標本化間隔以下の変位の推定法の1つとしてgradient-basedオプティカルフロー法がある [17]．この手法は，移動前後で超音波画像の輝度値・パターンが変化しないものと仮定し，画像中の任意の点の変位を画像の輝度勾配から推定するものである．
時刻tにおける画像上の関心点P(x, z)の輝度値をI(x, z; t)とし，微小時間 ΔT 後に点P が P’(x + Δx, z +Δz) に移動したとする．移動前後で輝度値が変わらないと仮定すると，(3)式の関係が成り立つ．
I ( x , z ; t ) = I ( x + Δx , z + Δz ; t + ΔT )(3)

右辺をTaylor展開すると(4)式が得られる．
I ( x, z; t ) = I ( x, z; t ) + Δx

∂I
∂I
∂I
+ Δz + ΔT
+L
∂z
∂t
∂x

(4)
2次以上の高次の項を無視して両辺をΔT で割り，
ΔT → 0とすると(5)式の関係が得られる．
vx

∂I
∂I ∂I
+ vz
+
=0
∂x
∂z ∂t

(5)
ここで，vx ( = vx(t))とvz ( = vz(t))はぞれぞれ，ラテラ
ル方向とビーム方向の速度である．

関心点Pの周りに関心領域RPを設定してRP内の点が全て同じ速度で移動すると仮定し，(6)式の二乗誤差 α を評価する．
∂I ∂I ⎞
⎛ ∂I
+ vz
+ ⎟
wx , z ⎜ vx∑∂
∂z ∂t ⎠x⎝x , z∈R P

α=

∂I ∂I

x , z ∈R P

(6)
ここで，wx,z は誤差を評価する際の重みであり，関心点Pの輝度値を用いた．αを最小化するために，αをvxとvzで偏微分して零とおいて連立方程式を解く
ことにより，ビーム方向の速度vz が(7)式により得
られる．

x, z

x, z

x , z∈ R P⎛ ∂I ⎞
wx , z ⎜ ⎟∑⎝ ∂z ⎠
x , z∈ R P

vˆx =

2

∂I ∂I

∑ w ∂z ∂t ∑ w ∂x ∂z
2− ∑ wx , z
x , z∈ R P⎛∂I ∂I ⎞
⎜ ∑ wx , z⎟
2⎜∂x ∂z ⎟⎠⎛ ∂I ⎞ ⎝ x , z∈R P∑ wx, z ⎜⎝ ∂x ⎟⎠ −2
x , z∈ R P⎛ ∂I ⎞w⎜⎟∑ x, z ⎝ ∂z ⎠x , z∈ R P

∂I ∂I
∂x ∂t

2

(8)
vˆ x (t ) を時間積分することによりラテラル方向の変位 Δxˆ (t ) が得られる．
２.３.２ 2次元相関法

∂I ∂I

∂I ∂I

∑ w ∂x ∂t ∑ w ∂x ∂z

x , z∈ R P

x, z

x , z∈R P⎛ ∂I ⎞vˆz =∑ w ⎜⎝ ∂x ⎟⎠x , z∈R P

x, z

2− ∑ wx , z
x , z∈ R P

∂I ∂I
∂z ∂t

x, z

⎛
∂I ∂I ⎞
⎜ ∑ wx , z⎟
2⎜∂x ∂z ⎟⎠⎛ ∂I ⎞ ⎝ x , z∈R P∑ wx, z ⎜⎝ ∂z ⎟⎠ −2
x z∈R P⎛ ∂I ⎞w⎜⎟∑ x, z ⎝ ∂x ⎠x , z∈R P

2(7)
この場合も，推定された速度 vˆ z (t ) を時間積分する
ことによりビーム方向の変位 Δzˆ (t ) が得られる．

超音波RF信号の2次元相互相関関数を用いた変位推定法についても検討を行った．時刻 t における関心点P(x, z)の超音波RF信号をs(x, z; t)と定義すると，時刻 t と t + ΔT における超音波RF信号間の相互相関関数 r(mx, mz; t) は(9)式で表される．
r ( mx , m z ; t ) =σ (t ) =2∑ s ( x , z ; t ) ⋅ s ( x + m ΔX , z + m ΔZ ; t + ΔT )x

x , z∈ R Pz

σ (t )σ (t + ΔT )∑ s ( x, z ; t )(9)2

x , z∈ R P

２.３ 超音波ビームと直交する方向 (ラテラル方向) の変位計測法
本研究におけるラテラル方向変位を計測する
目的は，動脈壁がその長軸方向に変位する場合に，その長軸方向変位に対応して関心領域を追随させ，常に同一領域の径方向ひずみを推定することである．したがって，変位計測の分解能がラテラル方向のデータ間隔 ( = 超音波ビームの走査間隔) に依存しても良い．
２.３.１ Gradient-Based オプティカルフロー法
2次元画像に対するgradient-basedオプティカル
フロー法では，輝度の時間方向，空間方向の勾配(空間微分) をもとに2次元の速度 vx，vz を同時に推定することができる．超音波ビーム方向の速度・変位推定の場合と同様に，(6)式をvx とvz について偏微分して零と置くことにより得られる連立方程式をラテラル方向速度vx について解くと(8)式のようにvx が推定できる．

σ 2 (t + ΔT ) = ∑ s ( x + m x ΔX , z + m z ΔZ ; t + ΔT )

2

(10)
(11)

x , z∈R p

ここで，mx とmz はラテラル方向とビーム方向のラグ，ΔX は超音波ビームの間隔である．相互相関関数 r(mx, mz; t) を最大とするラグ m̂ x ， m̂z を探索する
ことにより，ラテラル方向の速度vx(t)が(12)式により得られる．
vˆx (t ) =mˆ x ΔX
ΔT(12)(12)式から分かるように，速度推定の分解能は，超音波ビームの走査間隔に依存する．(12)式の速度
vx(t)を時間積分することによりラテラル方向の変位
Δx(t)が得られる．
２.３.３ Lateral Modulation 法
Jensenらは，受信ビーム形成時にラテラル方向に音場を変調した場合，その変調成分の位相変化が対象物のラテラル方向の変位に依存することを利用したラテラル方向変位推定方法を開発した[18]．この手法では，ラテラル方向の変調成分の位相を用いているため，2次元相関法のように超音波
ビーム間隔に依存することなくラテラル方向の変位を推定することができる．

本研究では，送信時には平面波 (ステア角φ = 0度)，受信時には2つの異なるステア角度φ1，φ2の平面波が関心深さz = 20 mmにおいて重なるよう，それぞれビームフォーミングを行った．受信時のビームフォーミングにおいて，2つの平面波のステア
角度をφ1 = -18度とφ2 = 18度に設定した場合と，ステア角度をφ1 = 18度とφ2 = -18度に設定した場合のラテラル方向の変調成分の位相が90度異なるため，それぞれのステア角度の組み合わせで得られた送受信信号 fI (x, z; t)，fQ (x, z; t)を元に，時刻 t とt + ΔT間の対象物のラテラル方向変位による関心点P(x, z)の送受信信号の位相変化Δθx(t)は(13)式で表される．
LejΔθˆ ( t )x=が見られるものの，この程度の分解能の犠牲により，フレームレートは約 60 Hz から約 3500 Hz にまで飛躍的に向上している．

) f ( x + lΔX , z ; t )∑ f ( x + lΔX , z ; t + ΔT ・
*

l =−L
L

) f ( x + lΔX , z ; t )∑ f ( x + lΔX , z ; t + ΔT ・*

(13)
ここで，f (x, z; t) = fI (x, z; t) + j fQ (x, z; t)であり，L は(13)式の複素相互相関関数を算出する際に使用する
ラテラル方向の標本点の数 (片側) である．
(13)式により推定された位相変化Δθx(t)を元に，
ラテラル方向の速度vx(t)が(14)式により算出される．
l =−Lvˆx (t ) =λx Δθˆx (t )
2πΔT

(14)
ここで，λx は，ラテラル方向に変調した音場の1周期の距離である．(14)式により得られた速度 vˆ x (t ) を
時間積分することにより，ラテラル方向の変位
Δxˆ (t ) が得られる．
３．動脈の高速超音波イメージングに関する基礎検討
３.１ ファントムに対する送受信結果
本報告では，生体組織を模擬したゲル体中に張ったナイロンワイヤ (直径 0.1 mm) に対して送受信を行うことにより parallel beam forming の方位分解能を検討した．
図 2(a) と 2(b) はそれぞれ，通常のリニア走査および parallel beam forming により得られたファン
トムの B モード像である．4 つの異なる深さにワイヤが張ってあるが，いずれの深さにおいても通常のリニア走査の方が parallel beam forming に比べ方位分解能が若干優れていることが分かる．また，
parallel beam forming では，ワイヤとバックグラウ
ンドの輝度の差が減少しており，コントラストの劣化も見られる．図 2(c) は，図 2(a) と 2(b) の水色の線上の振幅分布を示す．parallel beam forming の方が通常のリニア走査に比べ若干の分解能の低下

図 2. Ultrasonic B-mode images obtained by (a)conventional linear scanning and (b) parallel beam
forming. (c) Amplitude profiles along the cyan lines in(a) and (b).
３.２ 理論的考察
超音波プローブ表面の開口を，点音源の集合と仮定し，図 3 に模式的に示した開口から超音波を放射した際の，距離 z 上の点 P’(x, z) における音
圧 R(x) を考える．以下では z を固定し，x 方向の音圧分布について検討するため，z は定数として扱う．まず，図 3 の点 P(ξ, 0) にある点音源の音圧
r(ξ) の，P’ 点における音圧への寄与を考える．PP’
間の距離 d(ξ) = ((x − ξ)2 + z2)1/2 であり，点 P の音源
は点音源であるから，点 P の音源による点 P’ の音圧への寄与 a(ξ) は次式で表される．
⎧ 2π⎫exp⎨ jd (ξ )⎬
λ⎩⎭a(ξ ) = r (ξ ) ⋅d (ξ )(15)

ここで，λ は超音波の波長である．
(15) 式で示される r(ξ) は伝搬距離 d(ξ) により減衰するが，距離 d(ξ) の最大値と最小値の差は，

伝搬距離による減衰 (1/d(ξ)) に対しては十分に小さいと仮定すると，(15) 式は(16) 式のように減衰(1/d(ξ)) を無視することができる．
⎧ 2π
⎫
(16)
a (ξ ) = r (ξ ) exp⎨ jd (ξ )⎬⎩ λ⎭
点 P’ における音圧 R(x) は，(16) 式を開口幅 Lで積分したものであり，(17) 式で表される．
R( x) = ∫L / 2− x− L / 2− x

⎧ 2πr (ξ ) exp⎨ j⎩ λ⎫( x − ξ )2 + z 2 ⎬ d ξ⎭(17)

実際の超音波プローブは配列型探触子であるため，
音源 r(ξ) は離散的に存在する．(17) 式を離散系で表すと(18) 式となる．
N e / 2−n

R ( n) =∑i =− Ne / 2−n⎧ 2πr (i ) exp ⎨ j⎩ λ⎫(i − n ) 2 Δ x 2 + z 2 ⎬⎭(18)

ここで，R(n) ≡ R(nΔx)，r(i) ≡ r(iΔx) であり，Δx と
Ne はそれぞれ，素子間隔と開口の素子数である．

図 3. Assigned coordinates.

(21) 式を離散系で表すと(22) 式となる．
⎧ π

Ne −n

⎫

(22)

∑ r (i) exp⎨⎩ j λz (i − 2n)iΔx ⎬⎭R ( n) =2i=− Ne −nさらに，πξ2/(zλ) << 1 であり，πξ2/(zλ) による位相の変化は無視できると仮定すると，(21) 式は(23)式のように近似できる(Fraunhofer 近似) [18],[19]．
R( x) = ∫L / 2− x− L / 2− x

⎧ 2π ⎫r (ξ ) exp ⎨− jxξ ⎬ d ξ
λz ⎭
⎩

(23)

(23) 式は，1/(λz) でスケール化されているが r(ξ) のフーリエ変換に対応しており，音源の音圧分布と音場分布の対応付けが容易になる．例えば，r(ξ) =cos(πξ /L) とすることは (アポダイゼーション)，フーリエ変換において矩形窓ではなくハニング窓を用いることに対応し，メインローブの幅は広がるがサイドローブを抑える効果があることなどが容易に推察できる．ただし，ξ2/λ ≈ 103 m の場合において z >> 1600 m という条件が必要である[19]．本
報告における条件 ξ2/λ ≈ 60 m について考えると，z>> 96 m となるため，本報告の条件下においてFraunhofer 近似の精度は良くない．
本報告における parallel beam forming では，送信時には焦点を設けないため，送信音場は(17) 式
により表現できる．しかし，受信時に点(x0, z) に関して焦点を形成した場合には，受信音場 R’(x) は(24) 式で与えられる．
R' ( x) = ∫L / 2− x− L / 2− xr ' (ξ )

(

(24)

)

(17) 式により R(x) が求められるが，良く使用される近似について述べておく．伝搬距離 d(ξ) は⎧ 2π⎫d (ξ ) − (ξ − x0 ) 2 + z 2 + z ⎬ d ξ× exp⎨ j
λ⎩⎭(x − ξ )z2

ここで，r’(ξ) はプローブ表面での受信感度分布である．離散系における(24) 式の表現は，(24) 式のexp の乗数第 3 項は一定であるから無視すると(25)式で示される．

d (ξ ) = z 1 +

2

(19)

であり，(x− ξ)2/z2 が 1 よりも十分小さい場合には，(19) 式は(20) 式のように近似できる(para-axial 近似) [19]．
d (ξ ) = z +(x − ξ )2
2z

(20)

2

(20) 式における z + x /(2z) は，(17) 式の積分においては定数であるから無視すると，(17) 式は(21) 式のように変形できる．
R( x) = ∫L / 2− x− L / 2− x

⎧ π⎫r (ξ ) exp ⎨ j (ξ − 2 x)ξ ⎬ d ξ
⎩ λz
⎭

(21)

R ' ( n) =⎧ 2π
r ' (i ) exp⎨ j⎩ λi = − N e / 2− n
Ne / 2−n∑− (iΔx − x0 ) 2 + z 2( (i − n) Δx + z
2

2

2

(25)

)}

ここで，x0 は焦点の x 軸方向の位置である．
(24) 式に近似を適用すると，(23) 式に対応して，
R' ( x) = ∫L / 2− x− L / 2− x

⎧ 2π⎫r ' (ξ ) exp ⎨− j
( x − x0 )ξ ⎬ d ξ
λz
⎩
⎭

(26)

が得られる．(26) 式は para-axial 近似のみから導かれる．(22)式などと同様に，(26) 式の離散形は(27)式で示される．
R ' ( n) =2π⎧⎫
r ' (i ) exp ⎨− j
(i Δ x − x 0 )i Δ x ⎬λz⎩⎭i =− N e / 2−nN e / 2− n

∑

(27)

図 4 は，使用した装置と同じパラメータ λ =0.15 mm，Δx = 0.2 mm，Ne = 72，r(ξ) = 1，r’(ξ) =cos(πξ/L) として，z = 15 mm と z = 50 mm における音場を計算したものである(図 2(c) は z = 13 mm)．
図 4(a) には焦点を形成しない場合の音場 |R(x)| を，
近似なしおよび para-axial 近似あり，の場合についてそれぞれ算出したものを示してある．

れるように，parallel beam forming によるビーム径の拡大は主にピークの裾付近に見られ，図 2 の実
験結果と同様の傾向であった．また，図 4 では，(1) と(2) に z = 15 mm と z = 50 mm における分布をそれぞれ示してあるが，アポダイゼーションなしの場合 (|R(x)|) には z = 15 mm の近い距離での paraaxial 近似の精度があまり良くないことが分かる．
アポダイゼーション有の場合 (|R’(x)|) に近似精度が向上するのは，(19) 式において(x − ξ)2/z2 が大きい領域にある音源の寄与が小さくなるからである．
４．動脈壁変位計測法に関する基礎実験結果
４.１ 基礎実験システム
散乱体 (グラファイト粉末) を重量比 5%で混入させたシリコーンゴムで作製された模擬血管を用いて評価実験を行った．ファントムの内・外直径はそれぞれ 8 mm，10 mm であり，壁の弾性率は750 kPa である．循環系を模擬した水槽実験系において，拍動流ポンプを用いて模擬血管内圧の変化を発生させた．
また，測定は超音波診断装置の 10 MHz リニア型プローブを用いて行い，受信超音波 RF 信号を
40 MHz，16 bit で A/D 変換した．フレームレート
は 286 Hz である．内圧は図 5 に示すように圧力センサで計測した．

図 4. Theoretical sound fields at (1) z = 15 mm and (2) z= 50 mm. (a) Unfocused fields |R(x)| with and withoutpara-axial approximation. (b) Focused fields |R’(x)|(focus: x = 0 mm and z = 15 or 50 mm) with and withoutpara-axial approximation. Sound fields |R’(x)|2 and|R(x)||R’(x)| generated by conventional linear scan andparallel beam forming without approximation (c) andwith para-axial approximation (d).
図 4(b-1) と 4(b-2) はそれぞれ，x0 = 0，z = 15
mm および x0 = 0，z = 50 mm に焦点を形成した場
合の音場分布|R’(x)|である．図 4(a) と 4(b) の音場
分布から，通常のリニア走査の送受信音場|R’(x)|2および本報告における parallel beam forming の送受
信音場|R(x)||R’(x)| を算出した結果を図 4(c)に示す(図 4(d) は para-axial 近似による)．図 4(c-1) に示さ

図 5. Experimental system for basic experiments.
４.２ 超音波ビーム方向の変位計測結果
４.２.１ 位相を用いた手法による計測結果
図 6 は，模擬血管の長軸 B モード断層像である．各超音波ビーム位置で得られた超音波 RF 信号に 2.2.1 節の手法を適用して後壁内の変位分布を推定した．(1)式の K で決定される複素相関関数を算出する窓の長さは，使用した超音波パルスの-20 dB
幅を参照し，0.25 mm に設定した．また，超音波ビーム 3 本 (±1 本) の複素相関関数を加算平均している [11]．

図 7(a)は，46 本の超音波ビームそれぞれに沿った最大変位量の分布を算出し，それらの平均値と標準偏差を深さごとに点と縦棒で示したものである．図 7(a)のように得られた壁内変位分布を，超音波ビーム方向 (ファントム径方向) に空間微分することにより，図 7(b)の壁内径方向ひずみ分布を得た．実線は理論値であり，別途計測したファ
ントム壁の弾性率 E = 750 kPa と圧力センサで計測された内圧変化 Δp から(28)式により算出した [19]．
ε r (r ) = −ri 2 ro2
3Δp2
2 (ro − ri 2 )r 2 E(28)
計測した径方向のひずみの平均値の，理論値からの誤差は 12%，標準偏差は 14.1%であった．

(7)式に基づき超音波ビーム方向の速度 vz(t)を推定する際の関心領域 RP のビーム方向の大きさは，(1)式の位相を用いた場合と同様，使用した超音波
パルスの-20 dB 幅を参照し，0.25 mm に設定した．
図 8(1)は，複数の異なる関心領域 RP のラテラル方
向の幅を用いて推定したビーム方向変位Δx(t)である．図 8(1-a)における関心領域 RP のラテラル方向の幅 (±0.2 mm)は，図 7(a)の変位分布を推定した場合と同じであるが，図 7(a)に比べ推定した変位分
布の標準偏差が大きいことが分かる．変位分布の標準偏差を十分に減少させるためには，超音波ビームの-20 dB 幅のほぼ 2 倍に対応する±0.8 mm 程
度以上の関心領域 RP のラテラル方向の幅が必要であった．また，図 8(2)に示す径方向ひずみ分布の標準偏差を，位相を用いた推定結果と同程度に抑えるためには，±1.6 mm の関心領域 RP のラテラル方向の幅が必要であった．その場合でも，推定されたひずみの大きさは全体的には理論値と同程度であったものの，管外側の方が推定されたひずみ量が大きく，プロファイル的には理論値と異なる結果となった．

図 6. Longitudinal ultrasonic B-mode image of aphantom.
４.２.２ Gradient-Based オプティカルフローによる計測結果
2 次元画像に対するオプティカルフロー法は，位相を用いた方法と異なり 2 次元変位Δx(t)，Δz(t)が推定可能であるため，超音波ビーム方向の変位計測精度が十分であれば，オプティカルフロー法
のみで径方向ひずみ分布と血管長軸方向のトラッキングを同時に実現できる．本節では，標本点間隔以下の変位を推定できる gradient-based オプティ
カルフロー法の超音波ビーム方向の変位計測精度について検討を行った．

図 7. (a) Displacements relative to that of the luminalinterface and (b) radial strains of the posterior wall thephantom estimated using phases of ultrasonic RF echoes.
Plots and vertical bars show means and standard
deviations of the estimates obtained at all ultrasonicbeams.
４.３ ラテラル方向の変位計測結果

本実験では，ファントムには内圧変化を与えるのみで長軸方向には固定されている．非常にプリミティブな検討ではあるが本節では，ファント
ムに径方向変位・ひずみのみを与えた場合に，ラテラル方向変位の推定結果が零となるかどうかを評価することにより，各推定法の安定性について検討を行った．

出されており，この手法によるラテラル方向変位の推定が難しいことが分かる．
４.３.２ 2 次元相関法による推定結果
超音波 RF 信号間の 2 次元相互相関を用いた変位の推定結果を図 10 に示す．変位を計測した点の数は，(全ビーム本数) / (関心領域 RP に含まれるビーム本数) = 2 である．関心領域 RP のラテラル方向の幅は，gradient-based オプティカルフロー法を用いた場合に変位推定結果の標準偏差が十分に低下した±1.6 mm に設定した．RP のビーム方向の幅は，
図 10(1)と 10(2)においてそれぞれ 0.25 mm と 0.75
mm に設定した．図 10(a)と 10(b)はそれぞれ， 内圧上昇前後の B モード断層像と関心点を示している．
図 10(1)に示されるように，関心領域 RP のビー
ム方向の幅が 0.25 mm (使用した超音波パルスの-20dB 幅程度)の場合には，ファントムはその長軸方向には静止しているにも関わらず，ラテラル方向 ( =長軸方向) の変位が検出されていることが分かる(赤線が計測点)．図 10(2)に示されるように，関心
領域 RP のビーム方向の幅を 0.75 mm に設定することにより，関心領域内のパターンがより unique となり，安定した結果が得られた．

図 8. Axial (radial) displacements (1) and strains (2) of
the phantom estimated by gradient-based optical flow.
Axial size of the kernel was set at 0.25 mm. Lateral sizesof kernels were set at (a) ± 0.2 mm, (b) ± 0.4 mm, (c) ±0.8 mm, and (d) ± 1.6 mm.
４.３.１ Gradient-Based オプティカルフロー法による推定結果図 9 は，gradient-based オプティカルフロー法
によるラテラル方向変位 Δx(t)の推定結果である．
プロットと縦棒はそれぞれ，46 本の超音波ビームにおいて推定されたラテラル方向最大変位の平均値と標準偏差を示している．関心領域 RP のラテラ
ル方向の幅を大きくすることにより，標準偏差を減少させることはできるが，ラテラル方向幅±1.6mm の場合でも標準偏差は大きく，また，ファント
ムはラテラル方向には静止 (ラテラル方向の変位零) しているにも関わらず 0.10 mm 程度の変位が検図 9. Lateral displacements (1) and strains (radial) (2) of
the phantom estimated by gradient-based optical flow.
Axial size of the kernel was set at 0.25 mm. Lateral sizesof kernels were set at (a) ± 0.2 mm, (b) ± 0.4 mm, (c) ±0.8 mm, and (d) ± 1.6 mm.
４.３.３ Lateral Modulation 法による推定結果
まず，点散乱体 (0.1 mm のワイヤ 4 本) に対して送受信を行い，送受信音場の評価を行った．図11 はその結果である．図 11(a)と 11(b)はそれぞれ，
ラテラル方向変調の in-phase 成分と quadrature 成分

を示しており，散乱体の位置に見られるパターンのラテラル方向の変調成分の位相が図 11(a)と 11(b)で 90 度異なっていることが分かる．しかし，点拡がり関数としては通常のフォーカシングに比べ大きくなってしまっており，空間分解能が劣化してしまっていることが分かる．
図 12 は，Lateral Modulation 法におけるビーム
フォーミングにより得られたファントムの B モード像である．図 12 から分かるように，イメージン
グのための十分な空間分解能が得られたとは言い難く，この手法については今後の検討が必要である．

図 11. Ultrasonic echoes from 4 nylon wires (diameter:
0.1 mm) embedded in agar obtained by the lateralmodulation method. Echo signals from (a) in-phase and(b) quadrature beamformers.

図 12. B-mode image of the phantom obtained by thelateral modulation method.
５．in vivo 計測結果

図 10. Axial and lateral displacements of the phantom
estimated by 2-D correlation of RF echoes. Axial sizes ofkernels were set at (1) 0.25 mm and (2) 0.75 mm. Lateralsize of the kernel was set at ± 1.6 mm. (a) At minimum
internal pressure. (b) At maximum internal pressure.

基礎実験による検討をもとに，ラテラル方向のトラッキングには 2 次元相関法を，超音波ビー
ム方向の変位計測には位相を用いた手法を使用し，ヒト頸動脈における in vivo 計測を行った．図 13(a)と 13(b)は，心電図 R 波のタイミングおよび心電図
R 波から 0.3 秒後の B モード断層像上に，関心点を赤線で示したものである．図 13(a)と 13(b)において，動脈壁は収縮期に 0.5 mm 程度頭側に移動していることが B モード断層像から見て取れるが，設定した関心点はその動きを良く追従していることが分かる．図 13(c)は，受信超音波 RF 信号の位相を用いて動脈壁のひずみ分布を推定したものである．

図 13(c)から，心拍による血圧の上昇により，-10%程度の径方向ひずみ (壁厚が減少) が発生していることが分かる．また，動脈壁が長軸方向に 0.6 mm
程度伸張していることが確認された．

図 13. In vivo experimental results. (a) Initial positions ofpoints of interest shown by the red lines. (b) Positions of
the points of interest at 0.3 s after the R-wave ofelectrocardiogram. (c) Spatial distribution of radial(axial) strains of the arterial wall estimated using phasesof ultrasonic echoes.
６．まとめ
本研究では，動脈壁の変位計測法に関する検討を行った．心臓の収縮や圧力波の伝搬による影響などにより，動脈がその長軸方向に移動する場
合もあることから，長軸方向(ラテラル方向)の動き
をトラッキングしながら径方向変位の分布を高精度に計測し，壁の径方向ひずみ分布を推定する方法を検討した．その結果，ラテラル方向のトラッキングには超音波 RF 信号の 2 次元相関法，ビーム方向の変位計測には位相を用いた変位推定法を用いた場合に最も良い結果が得られた．ただし，ラテラル方向の変位計測に関する検討は，ラテラル
方向には静止しているファントムについて行ったのみであり，今後さらなる検討が必要である．

参考文献[1] Lee RT, Grodzinsky AJ, Frank E, Kamm RD, andSchoen FJ. Structure-dependent dynamic behaviorof fibrous caps from human atherosclerotic plaques.
Cirdulation 83, 1764-1770, 1991.
[2] Loree HM, Grodzinsky AJ, Park SY, Gibson LJ, and
Lee RT. Static circumferential tangential modulusof human atherosclerotic tissue. J Biomech 27, 195204, 1994.
[3] Benthin M, Dahl P, Ruzicka R, and Lindström K.
Calculation of pulse-wave velocity using crosscorrelation -Effects of reflexes in the arterial tree.
Ultrasound Med Biol 17, 461-469, 1991.

[4] Kanai H, Kawabe K, Takano M, Murata R,
Chubachi N, and Koiwa Y. New method forevaluating local pulse wave velocity by measuring
vibrations on aortic wall. Electron Lett 30, 534-536,
1993.
[5] Hoeks APG, Ruissen CJ, Hick P, and Reneman RS.
Transcutaneous detection of relative changes inartery diameter. Ultrasound Med Biol 11, 51-59,
1985.
[6] Länne T, Stale H, Bengtsson H, Gustafsson D,
Bergqvist D, Sonesson B, Lecerof H, and Dahl P.
Noninvasive measurement of diameter changes inthe distal abdominal aorta in man. Ultrasound Med
Biol 18, 451-457, 1992.
[7] Kanai H, Sato M, Koiwa Y, and Chubachi N.
Transcutaneous measurement and spectrum analysisof heart wall vibrations. IEEE Trans Ultrason
Ferroelect Freq Contr 43, 791-810, 1996.
[8] Hasegawa H, Kanai H, Koiwa Y, and Chubachi N.
Noninvasive evaluation of Poisson’s ratio of arterial
wall using ultrasound. Electron Lett 33, 340-342,
1997.
[9] Hasegawa H, Kanai H, Hoshimiya N, and Koiwa Y.
Evaluating the regional elastic modulus of acylindrical shell with nonuniform wall thickness. J
Med Ultrason 31, 81-90, 2004.
[10] Kanai H, Hasegawa H, Ichiki M, Tezuka F, andKoiwa Y. Elasticity imaging of atheroma withtranscutaneous ultrasound -preliminary study-.
Circulation 107, 3018-3021, 2003.
[11] Hasegawa H and Kanai H. Strain imaging ofarterial wall with translational motion compensationand error correction. IEEE 2007 Intern’l Ultrason
Symp Proc 860-863, 2007.
[12] Cinthio M, Ahlgren AR, Jansson T, Eriksson A,
Persson HW, and Lindström K. Evaluation of an
ultrasonic echo-tracking method for measurements ofarterial wall movements in two dimensions. IEEE TransUltrason Ferroelect Freq Contr 52, 1300-1311, 2005.
[13] Hasegawa H and Kanai H. Simultaneous imaging ofartery-wall strain and blood flow by high frame rateacquisition of RF signals. IEEE Trans Ultrason
Ferroelect Freq Contr 55, 2626-2639, 2008.
[14] 長谷川英之, 金井 浩. 超音波による動脈壁の変位計測に関する検討. 電子情報通信学会技術研究
報告 108, 19-24, 2008.
[15] Tanter M, Bercoff J, Sandrin L, and Fink M. Ultrafast
compound imaging for 2-D motion vector estimation:
application to transient elastography. IEEE TransUltrason Ferroelect Freq Contr 49, 1363-1374, 2002.
[16] Mahafza BR. Introduction to radar analysis. CRC
Press, Boca Raton, 1998.
[17] Mikic I, Krucinski S, and Thomas JD.
Segmentation and tracking in echocardiographic
sequences: active contours guided by optical flow
estimates. IEEE Trans Med Imaging 17, 274-284,
1998.
[18] Jensen JA and Munk P. A new method for
estimation of velocity vectors. IEEE Trans Ultrason
Ferroelect Freq Contr 45, 837-851, 1998.
[19] Goodman JW. Introduction to Fourier optics.
McGraw Hill, New York, 1968.
[20] Timoshenko SP and Goodier JN, Theory of
Elasticity. New York, McGraw Hill, 1970.