相関と波数ベクトルフィルタリングを用いた

平成２３年度 修 士 論 文

RF 相関と波数ベクトルフィルタリングを用いた

高精度ずり弾性波映像化

指導教員 山越 芳樹 教授

群馬大学大学院工学研究科

電気電子工学専攻

吉原 由貴

第 第 第

第 １ １ １ １章 章 章 章 序論 序論 序論 序論

近年、日本のがん死亡率は増加し、1981年には脳卒中と入れかわって死亡原因の第1位 となっている。日本におけるがん死亡数の増加の主な原因は、人口構成の高齢化、高齢者 人口の増加であると言われており、今後もますます増加すると考えられる。Fig.1-1に部位 別がんの死亡率を示す。がんの中でも特に女性では乳がん、男性では前立腺がんによる死 亡率が増加しており、生体表面のがんの定量的な診断が求められている。しかし、現在、

乳がん検査として用いられるマンモグラフィは、X線を用いるため、X線による被ばくなど 安全面で問題がある。また、得られる画像の読影が難しく正確な読影は医師や検査技師の 経験に頼る部分が大きいため、誤って診断されるケースも多い。そのため、安全かつ定量 的ながんの診断法を確立するために数々の研究がなされており、なかでも、正常な組織に 比べてがん組織が硬いという特徴を利用した組織弾性計測が近年注目を浴びている。

生体軟組織などの比較的柔らかい物体の表面から周波数1kHz程度までの低周波振動を 加えると、その放射エネルギーの大部分は生体中を横波として伝搬し、その伝搬速度や減 衰係数は、ずり粘弾性係数などのずれ粘弾性パラメータと関連があることが知られている。

また、生体組織のずり粘弾性特性は、組織を触ったときの硬さや触感と密接に関係してい る。そのため、生体組織について低周波振動の伝搬速度や減衰などが測定できれば、画像 などの視覚的な判断のみに頼ることなく、疾病の進行度の定量的な評価や早期発見が期待 でき、これらは組織の特性化のために有用である。しかし、生体組織の機械的構造は非常 に複雑であり、組織境界等で反射や屈折が生じるため、これが時として測定精度に影響を 与えてしまい、肝臓などの比較的一様な組織でしか伝搬速度を精度良く測定できないのが 現状であるといえる。さらにサンプルボリューム内の超音波散乱体構造の変化により推定 誤差が生じ、スパイク状のノイズが発生し画像が劣化する問題がある。そのため、非一様 かつ複雑な境界面をもつ組織においても、精度良く組織内部の粘弾性を測定できるシステ ムが求められている。さらに臨床においては測定結果の明快かつ信頼たる画像化が求めら れている。

そこで本研究では測定対象を体表近くの組織、特に乳腺および前立腺とし、本論文では 生体内ずり弾性波の映像化としてＲＦ相関によるずり弾性波変位計測、波数ベクトルフィ ルタリングを用いたずり弾性波映像系を提案した。さらにシミュレーションおよび寒天フ ァントムを用いた実験により提案手法の評価を行った。

本論文の構成は以下のようになっている。

第２章では、生体組織内部における低周波振動の伝搬について示す。

第３章では、ずり弾性波映像化の基本原理について示す。

第４章では、シミュレーションによる提案手法の評価について示す。

第５章では、寒天ファントムを用いた実験による提案手法の評価について示す。

第６章では、結論と今後の課題について示す。

Fig.1-1 部位別がん死亡率

前立腺 前立腺 前立腺 前立腺

乳房 乳房

乳房 乳房

第 第 第

第 ２ ２ ２章 ２ 章 章 章 生体軟組織内部 生体軟組織内部 生体軟組織内部 生体軟組織内部における における 低周波振動 における における 低周波振動 低周波振動 低周波振動 の の の伝搬 の 伝搬 伝搬 伝搬

本章では、生体軟組織内部における低周波振動の伝搬について示す。

生体組織の粘弾性パラメータと低周波振動の伝搬速度および減衰の関係について以下に 示す。

外部から媒質に振動を与えると、その振動は一般に縦波・横波として伝搬するが、生体 のような粘弾性媒質中では、Hookeの法則が成り立つVoigtモデルと仮定することによりこ の縦波・横波の伝播速度および減衰係数は次式で与えられる。[Ref.1]

① 縦波

伝搬速度 ：

[ ]

Re

v

v

l

g

= ω

^（2-1-1）

減衰係数 ：

α

_l

= − Im[ ] g

^（2-1-2）

ただし、

( )

1

2 2

2

g ρω

v

µ λ

 

 

=  

+

 

 

（2-1-3）

② 横波

伝搬速度 ：

[ ]

Re

v

v

t

h

= ω

^（2-1-4）

減衰係数 ：

α

_t

= − Im[ ] h

^（2-1-5）

ただし、

1

2 2

h ρω

v

µ

 

=  

 

（2-1-6）

1

j

v 2

µ µ = + ω µ λ λ = +

1

j ω λ

v 2

µ

1^{：ずれ弾性係数}

λ

1^{：体積弾性係数}

µ

2^{：ずれ粘性係数}

λ

2^{：体積粘性係数}

ρ

^：密度

ω

v^{：振動角周波数}

Re[]

^、

Im[]

^：

[]

内の複素数の実数部、虚数部

また、生体表面近くにはこれら縦波や横波の他に表面波が存在するが、この伝搬速度はほ ぼ横波の伝搬速度に等しいことが知られている。上記の波動の中で、縦波は圧縮性の波で あり、媒質を圧縮することにより伝搬する。一方、横波は非圧縮性の波であり、媒質を等 体積のまま、横方向に変形させながら伝搬していくため、ずり波とも呼ばれている。ここ で、周波数が1 kHz程度以下の低周波振動であると、外部から与えた振動のエネルギーは、

そのほとんどが横波に変換されると考えられている。[Ref.2]

ここで、（2-1-4）式、（2-1-5）式で与えられる横波の伝搬速度と減衰係数をずり粘弾性パラ

メータを用いて書くと、

( )

( )

2 2 2

1 2

2 2 2

1 1 2

2

_v

t

v

v µ ω µ

ρ µ µ ω µ

= +

+ +

（2-1-7）

( )

( )

2 2 2 2

1 2 1

2 2 2

1 2

2

v v

t

v

ρω µ ω µ µ

α µ ω µ

+ −

= +

^（2-1-8）

となる。

したがって、もし、媒質の弾性が粘性にまさり、

1 v 2

µ >> ω µ

の関係が成り立つときには、

1 1

v

t

µ

≅ ρ

^（2-1-9）

1

0

α

t

≅

^（2-1-10）

となり、伝播速度は、単にずれ弾性係数と媒質の密度のみの関数となる。このとき、

µ

1^が

大きいということは、媒質が硬いということであり、硬い媒質ほど伝搬速度は速くなる。

一方、媒質の粘性が弾性にまさり、

1 v 2

µ << ω µ

の関係が成り立つときには、

2 2

2

_v

v

t

ω µ

≅ ρ

^（2-1-11）

2

2

2 v t

α ρω

≅ µ

^（2-1-12）

となり、

v

t2^、

α

_t2ともずれ粘性係数と密度の関数になり、この場合

2

v

t ^、

α

_t2^{の周波数依存}

性（分散性）が現れてくる。

以下に弾性体と粘弾性体の周波数別伝搬速度を示す。

Fig.2-1 弾性体と粘弾性体の周波数別伝搬速度

0 1 2 3 4 5 6

100 200 300 400 500 600 700 800 900 1000

弾性率2.26kPa,粘 性率2.38Pa・s 弾性のみの場合

（2.26kPa）

伝搬速度[m/sec]

加振周波数[Hz]

第 第 第

第 ３ ３ ３章 ３ 章 章 章 ずり ずり弾性波映像化 ずり ずり 弾性波映像化 弾性波映像化 の 弾性波映像化 の の の基本原理 基本原理 基本原理 基本原理

本章では、３－１で超音波パルスドプラ法による組織内振動伝搬計測の基本原理につい て示す。さらに

３－２ RF相関による微小変位計測の基本原理

３－３ 波数ベクトルによる伝搬速度推定法の基本原理

３－４ 波数ベクトルフィルタリングによるずり弾性波映像化の基本原理 ３－５ 信号処理フロー

を示す。

３－１ 超音波パルスドプラ法による組織内振動伝搬計測の基本原理

組織内振動伝搬計測は、組織表面から振動を印加することで組織内に振動を励起させ、

内部を伝搬する振動を超音波で計測するものである。これは、組織内部を多数の超音波散 乱体と考えると、組織内部に超音波を送波し、超音波散乱体から反射してくる超音波がド ップラー効果によって周波数変調を受けていることに着目したものである。したがって、

超音波散乱体から反射した超音波を直交検波することで得られるドップラー信号から組織 内部を伝搬する振動を推定することができる。

今、Fig.3-1-1に示すような超音波トランスデューサに遠ざかる方向に、周波数

f

v^、速度

v t ( )

で振動する超音波散乱体に対して超音波パルスを送波する場合を考える。

Fig.3-1-1 計測モデル

散乱体の運動ξ(t)は次式で表すことができる。

( ) t

0

cos(2 f t

_v

)

ξ = ξ π − ϕ

^{（3-１-1）}

ただし

ξ

0^{：振動振幅}

ϕ

^{：初期位相}

f

0

( )

v t f

v

超音波 トランスデューサ

超音波散乱体

この時、超音波散乱体に反射した超音波の周波数

f

^は

( )

0

c v t

f f

c

= −

^(3-1-2)

f

0^{：超音波の中心周波数}

c

^：音速

この反射波が超音波トランスデューサで受信される時の周波数

f '

^は

' ( )

f c f

c v t

= +

^(3-1-3)

(3-1-2),(3-1-3)式より

( ) ( ) ( )

( )

0 0

c v t c v t

f c f f

c v t c c v t

− −

′ = × =

+ +

^(3-1-4)

したがって、反射超音波のドプラ周波数シフト

∆ f

^は

( ) ( ) ( )

( )

0 0 0 0

2

c v t v t

f f f f f f

c v t c v t

− −

∆ = ′ − = − =

+ +

^(3-1-5)

となる。

超音波ドプラ法で組織内の速度を観測する場合、組織内での音速は約1500m/secであり、

そ れ と 比 較 し て 観 測 し よ う と す る 組 織 内 の 速 度 は 1~十 数 m/sec と 微 小 で あ る の で 、

( )

c v t

^となり、(3-1-5)式は次式のように近似することができる。

( )

0

f 2v t f

∆ − c

^(3-1-6)

この時、超音波の位相変化

∆ φ

^は

( ) ( )

0

0

2 4

4 ( )

f dt f v t dt c

f t

c φ π

π π ξ

∆ = ∆

= −

= −

∫

∫

^(3-1-7)

となるので、この散乱体からの受信信号（RF信号）

y t ( )

^は

0

0 0

0

( ) ( ) sin(2 2 )

( ) sin(2 4 ( ) 2 )

( ) sin{2 ( 2 ( ) ) 2 }

u

u

u

y t A t f t k z

A t f t f t k z

c

A t f t t k z

c

π φ

π π ξ

π ξ

= + ∆ −

= − −

= − −

(3-1-8)

ただし

( ) A t

^：振幅

k

u ^{：超音波パルスの波数}

z

：トランスデューサ、散乱体間の距離

となる。

次にRF信号に、位相が互いに90度異なる超音波周波数成分を畳み込み積分し低域通過フ ィルターをかけ、QI信号を得る。

(ⅰ)I信号

RF信号にキャリア信号を乗算すると

0 0

0 0

0

2 ( )

'( ) ( ) sin{2 ( ) 2 }sin(2 )

4 ( ) 4 ( )

( ) {cos(4 2 ) cos( 2 )}

2

u

u u

I t A t f t t k z f t

c

f t f t

A t f t k z k z

c c

π ξ π

π ξ π ξ

π

= − −

= − − − − − −

(3-1-9)

となる。ここで

2 ω

0付近の信号を低域通過フィルターで除くと、

4

0

( )

( ) ( ) cos( 2 )}

2

^u

f t

I t A t k z

c

= − π ξ −

(3-1-10)

となりI信号を得る。

（ⅱ）Q信号

（ⅰ）と90度異なるキャリア信号を乗算すると

0 0

0 0

0

2 ( )

'( ) ( ) sin{2 ( ) 2 }cos(2 )

4 ( ) 4 ( )

( ) {sin(4 2 ) sin( 2 )}

2

u

u u

Q t A t f t t k z f t

c

f t f t

A t f t k z k z

c c

π ξ π

π ξ π ξ

π

= − −

= − − + − −

(3-1-11)

となる。（ⅰ）と同様に低域通過フィルターを用いると

4

0

( )

( ) ( ) sin{ 2 }

2

^u

f t

Q t A t k z

c

= − π ξ −

(3-1-12)

となりQ信号を得る。

従来法であるArc-tan法を用いた変位推定はQI信号より

1 0

( ) tan ( ) 2

4 ( )

c Q t

t n

f I t

ξ π

π

⁻

   

 

= −

  

±



   

  ^(3-1-13)

となる。

ただし

tan

⁻1関数の主値の範囲を考慮し、

I

^{の符号の変化時に}

± 2n π

のオフセットを加える。

３－２ RF相関による微小変位計測の基本原理

ずり弾性波により散乱体を振動させると、(3-1-8)式より、送波パルスごとのRF信号にわ ずかな位相差が生じる。この位相差は散乱体の移動量を表したものである。よってパルス 間のRF信号で相関をとり、位相差を計測することで散乱体の移動量を推定する。

（１） RF相関による微小変位計測

今、ある散乱体からのRF信号模式図をFig.3-2-1に示す。

Fig.3-2-1 RF信号模式図

このRF信号をパルスごとに区切り二次元的に表すと、(3-1-8)式より

0

( ', )

( ', ) ( ', ) sin{2 ( 2 t ) 2

_u

( ')}

y t A t f k z t

c

τ = τ π τ − ξ τ − + ϕ

(3-2-1)

ただし

'

t

^{：パルストリガ時間}

τ

^{：超音波遅延時間}

( ') t

ϕ

^{：初期位相}

1

^st

pulse 2

^nd

pulse 3

^rd

pulse 4

^th

pulse ….

t

パルストリガ時間[msec]

µ

いま、基本となる超音波パルス

i

番目のパルストリガ時間を

t '

_i、次のパルストリガ時間を

'

_i 1

t

₊ ^{とし、深さ}

z

_l を中心に振動する散乱体からのRF信号を次式のように表す。

0

0

( ' , )

( ' , ) ( ' , ) sin{2 ( 2 ) 2 ( ')}

( ' , ) sin{2 ( ' , )}

i l

i l i l l u l

i l l i l

y t A t f t k z t

c

A t f C t

ξ τ

τ τ π τ ϕ

τ π τ τ

= − − +

= −

(3-2-2)

1

1 1 0

1 0 1

( ' , )

( ' , ) ( ' , ) sin{2 ( 2 ) 2 ( ')}

( ' , ) sin{2 ( ' , )}

i l

i l i l l u l

i l l i l

y t A t f t k z t

c

A t f C t z

ξ τ

τ τ π τ ϕ

τ π τ

+ + +

+ +

= − − +

= −

(3-2-3)

ただし

2

_l

l

z

τ = c

(3-2-4)

0

( ', )

( ', ) 4 t 2

_u

( ')

C t f k z t

c

τ = π ξ τ + + ϕ

(3-2-5)

とおく。

( ' , )

_i 1

y t

₊

τ

^{の位相をある遅延時間}

τ

shift^{分ずらした}

y t ( ' ,

_i+1

τ τ +

_shift

)

^と

y t ( ' , )

_i

τ

^{の二乗誤差を、}

τ

lを中心に求めると、二乗誤差

2

( ' ,

_{i shift}

, )

_l

e t τ τ

^は

/ 2

2 2

1 / 2

/ 2

0 / 2

2

1 0 1

( ' , , ) { ( ' , ) ( ' , )}

[ ( ' , ) sin{2 ( ' , )}

( ' , ) sin{2 ( ) ( ' , )}]

l

l l

l

nT

i shift l i i shift

nT nT

i i l

nT

i shift shift i l

e t y t y t d

A t f C t z

A t f C t z d

τ

τ τ

τ

τ τ τ τ τ τ

τ π τ

τ τ π τ τ τ

+

+

− +

−

+ +

= − +

= −

− + + −

∫

∫

(3-2-6)

となる。

ただし

2 2

n n

shift

τ τ τ

− ≤ ≤

τ

n^{：RF信号のシフト幅}

T

^{：超音波周期}

n

^{：任意の相関長さ}

ここで

τ

n^は1パルス間における最大移動量以上とする。

－

－

－

－ :::: y t ( ' , )

_i

τ

－ － － －

_:

_{y t ( ' , )}

_i+1

τ

--- --- --- ---: y t ( ' ,

_i+1

τ τ +

_shift

)

Fig.3-2-3

i

^番目のRF信号とシフトさせた

i + 1

^番目のRF信号

さらに

2

( ' ,

_{i shift}

, )

_l

e t τ τ

^{を正規化すると}

2 2

/ 2 2 / 2

( ' , , ) ( ' , , )

( ' , )

l

l

i shift l

normal i shift l nT

i nT

e t e t

y t d

τ τ

τ τ τ τ

+

τ τ

−

=

∫

(3-2-7)

となる。

ま た 、 あ る 時 間 差

τ '

^{に 対 し 、}

y t ( ' , )

_i

τ = y t ( ' ,

_i+1

τ τ + ')

が 成 り 立 つ 場 合 、 二 乗 誤 差

2

( ' , , )

normal i shift l

e t τ τ

^は

τ

shift ^{の 二 次 式 で 表 さ れ る 。 パ ル ス 間 の RF 信 号 は}

( ' , )

_i

( ' ,

_i 1

')

y t τ ≈ y t

₊

τ τ +

^{となるため、}

e

_normal²

( ' , t

_i

τ

_shift

, ) τ

_l ^{は次式で表される。}

2 2

, , ,

( ' , , )

normal i shift l i l shift i l shift i l

e t τ τ ≈ a τ + b τ + c

^(3-2-8)

ただし

,

,

,

,

, i l i l i l

a b c

^：定数

τ

Fig.3-2-4 RF信号の二乗誤差

(3-2-8)式の解

'

' ,

i l

t τ

τ

^は

ⁱ

^番目と

ⁱ + ¹

^番目のRF信号の位相差に相当する。よって(3-2-8)式 を解くことで

i

^番目と

i + 1

番目の位相差が得られる。

(3-2-8)式を二階微分すると

2

, ,

2

( ' , , )

i shift l

2

i l shift i l shift

e t τ τ a b

τ τ

∂ = +

∂

^(3-2-9)

となり、解

'

' ,

i l

t τ

τ

^は

, ' , ,

, ' ,

,

2 ' 0

' 2

i l

i l

i l t i l

i l t

i l

a b

b a

τ τ

τ τ

+ =

= −

(3-2-10)

と表される。

2 τn

−

shift

τ '

' ,

i l

t τ

τ

2

( ' , , )

normal i shift l

e t τ τ

2 τn

ここで

y t ( ' , )

_i

τ

_l ^と

y t ( ' , )

_i+1

τ

_l ^{の位相差は(3-2-2)、(3-2-3)式より}

1 1

1 0 1 0

sin { ( ' , )} sin { ( ' , )} {2

⁻

y t

_i+

τ

_l

−

⁻

y t

_i

τ

_l

= π τ f

_l

− C t ( ' , )} {2

_i+

z

_l

− π τ f

_l

− C t ( ' , )}

_i

z

_l

( ' , )

_i 1 _l

( ' , )

_{i l}

C t

₊

z C t z

= − +

(3-2-11)

(ただし

( ' , )

_{i l}

( ' , )

_i 1 _l

A t τ A t

₊

τ

⁾

(3-2-5)式より、(3-2-11)式は散乱体の振動変位量を表す。

よって深さ

z

l^{にある散乱体の}

t '

_i^、

t '

_i+1^間の変位

∆ ξ ( ' , ) t

_i

z

_l ^は

( ' , ) '

' ,

i l

i l t

t z

_τ

ξ τ

∆ =

^(3-2-12)

となる。

これより散乱体振動は

' '

' 0

' '

' ,

' 0

( ', ) ( ' , ) '

' '

i

i i

i l i

t t

l i l i

t

t t

t i

t

t z t z dt

τ

dt

ξ ξ

τ

=

=

=

=

= ∆

=

∫

∫

(3-2-13)

となる。

（２）適応的変位推定

（１）ではRF信号の相関による微小変位計測の原理を述べた。計測対象の散乱体が超音 波パルス間で

∆ ξ ( ' , ) t

_i

z

_l ^{振動した場合、(3-2-1)式、(3-2-12)式より、RF信号は}

1 ,

( ' , ) ( ' , ' )

i l

i i t

y t τ = y t

₊

τ τ +

_τ (3-2-14)

と表せる。

しかし実際にはサンプルボリューム内の散乱体が３次元的に変化することなどから初期位 相

ϕ ( ') t

が変化し、パルス間でのRF信号はわずかに変化する。この時の

' , '

_{i i} 1

t t

₊ ^{間の二乗誤} 差を

2

( ' , , )

normal i shift l

e t τ τ

^{とする。また、}

e

normal²

^{( ' ,} t

i

τ

shift

^{, )} τ

l ^{のモデル関数を}

e

true²

^{( ' ,} t

i

τ

shift

^{, )} τ

l

とし、次式で表す。

2 2

' , ' , ' ,

( ' , , )

i l i l i l

true i shift l t shift t shift t

e t τ τ = a

_τ

τ + b

_τ

τ + c

_τ (3-2-15)

ここで計測した RF 信号の二乗誤差

2

( ' , , )

normal i shift l

e t τ τ

^{とモデル関数}

e

true²

^{( ' ,} t

i

τ

shift

^{, )} τ

l ^の二

乗誤差をとることで、RF信号が変化し推定変位の誤差が大きい個所を特定することができ る。

2

( ' , , )

normal i shift l

e t τ τ

^と

e

true²

^{( ' ,} t

i

τ

shift

^{, )} τ

l の二乗誤差より推定誤差

2

( ' ,

_{i shift}

, )

_l

E t τ τ

^は

/ 2

2 2 2

2 / 2

/ 2

2 / 2

{ ( ' , , ) ( ' , , )}

( ' , , )

( ' , , )

n

n

n

n

normao i shift l true i shift l shift

i shift l

normao i shift l shift

e t e t d

E t

e t d

τ

τ

τ τ

τ τ τ τ τ

τ τ

τ τ τ

−

−

−

=

∫

∫

(3-2-16)

と表される。

推定誤差により推定個所を適応的に選択することで高精度な変位推定が可能となる。以 後、推定誤差が大きい部分を抜き出し、深さ方向に線形補間を行う変位推定法を適応的RF 相関法と呼ぶ。

３－３ 波数ベクトルによる伝播速度推定法の基本原理

今、加振器を用いて生体表面から生体内に振動を励起することを考える。計測モデルの

断面図をFig 3-3-1に、上面図をFig3-3-2に示す。加振源を原点とし、加振源からN種類の

周波数で加振した場合を考える。測定領域の微小区間において波は平面波で入射し、反射、

屈折はないものと仮定する。また、波はある周波数

f

v^に対して

L

^{波伝播するとし、}z軸に対 する各伝播方向の天位角を

θ

_l

( l = 1, 2,..., ) L

^{、方位角を}

φ

_l

^{( l} = ^{1, 2... ) L}

^{とする推定に用い} る領域の大きさを(開口長)をd、開口の中心を(X,Z,Y)とする。

Fig.3-3-1 計測モデル断面図

スキャン

θ

l

d

d X

…

加振器 トランスデューサ

Z

Fig.3-3-2 計測モデル上面図

また、加振周波数

f

v^の第l波の振動によるx方向の波数

k

lx

′

^、z方向の波数

k′

lz^、y方向の波 数

k

ly

′

^{はそれぞれ}

sin sin cos sin cos

lx l l l

lz l l

ly l l l

k k

k k

k k

θ φ

θ

θ φ

′ =

′ =

′ =

(3-3-1)

ただし

l

:

v

k 加振周波数 f ,第 l 波の波数

で表されるので式(3-3-1)より、ある測定点

( , , ) x z y

^{での振動振幅}

ξ ( , , , ) x z y t

^{は次式で表さ} れる。

1

1

( , ) cos{2 ( ' ' ' )}

cos{2 ( ' )}

L

l v lx lz ly l

l L

l v l l

l

t f t k x k z k y

f t

ξ δ π ϕ

δ π ϕ

=

=

= − + + +

= − ⋅ +

∑

∑

p

k p

(3-3-2)

ただし

: :

' ( , , ) :

( , , ) :

:

v

l v

l lx lz ky v

l v

f

f l

k k k f l

x z y

f l

δ

ϕ

′ ′ ′

=

= k p

　　　　　　　　　 加振周波数

　　　　　　　　 加振周波数 ,第 波の振動振幅 　　 加振周波数 ,第 波の波数ベクトル 　　　　　 測定位置の位置ベクトル

　　　　　　　 　加振周波数 ,第 波の初期位相

x

y 加振器

トランスジューサー

スキャン

φ

l

…

である。また、計測時間

T

を加振周波数の整数倍とし、

ξ ( , , , ) x z y t

^{を加振周波数帯域のみ} のフーリエ変換を行うことによって得られる複素振幅

^

( , f

_v

)

ξ

^{p は次式で表される。}

^

0

0 1

0 1

( , ) 2 ( , ) exp( 2 )

2 cos{2 ( ' )}exp( 2 )

1 [exp{ {2 ( ' )}} exp{ {2 ( ' )}]exp( 2 )

T

v v

T L

v l l v

l T L

v l l v l l v

l

f t j f t dt

T

f t j f t dt

T

j f t j f t j f t dt

T

ξ ξ π

π ϕ π

π ϕ π ϕ π

=

=

= −

= − ⋅ + −

− ⋅ + + − − ⋅ + −

∫

∫ ∑

∫ ∑

p p

k p

k p k p

　 =

1

exp{ ( ' )}

L

l l

l

j ϕ

=

= ∑ −

^k

⋅ +

^p (3-3-3)

また開口の中心位置の位置ベクトルをP

= ( , , ) X Z Y

^{とし、複素変位}

^

( , f

_v

)

ξ

^{p をx方向、z方}

向、ｙ方向の3次元フーリエ変換を行うことにより、次式を得る。

{ }

/ 2 / 2 / 2^

3 / 2 / 2 / 2

( , ,

_v

) 1

^{Y d Z d X d}

( ,

_v

) exp ( )

Y d Z d X d

P f f j dxdzdy

d

₋⁻ ₋⁻ ₋⁺

ξ

= ∫ ∫ ∫ − ⋅

k P p k p

' ' '

1

exp{ ( ' ) }sin sin sin

2 2 2

L

y ly

x lx z lz

l l

l

k k

k k k k

j c d c d c d

δ

=



−





−

 

−

  

= −

     

 

     

∑

^{k k P (3-3-4)}

ただし

:

d

開口長

( , ,

_v

)

P k P f

^は

k

_x

= k

_lx

'

^、

k

_z

= k

_lz

'

^、

k

_y

= k

_ly

'

においてピークを持ち、その点から離れるに従 ってsinc 関数的に振動しながら減衰していく。よって

P ( , , k P f

_v

)

のピーク位置よりずり弾 性波の波数ベクトル

k

_l

'

が推定される。さらに開口長の中心座標

( , , ) X Y Z

^{を移動させるこ} とにより、ずり弾性波の伝搬速度

v ( , P f l

_v

, )

^{は次式で表される。}

( ,

_v

, ) 2

^v

l

v f l ₌ π f P ′

k

^(3-3-5)

また方位角

φ ( , ^P f l

_v

, )

^、天位角

θ ^{( , P} f l

v

^{, )}

^は

1

'

( , , ) tan ( ) '

lx v

ly

f l k

φ P =

⁻

k

(3-3-6)

2 2

1

' '

( , , ) tan ( )

'

ly lx

l

lz

k k

f l k

θ P =

⁻

⁺

(3-3-7)

となる。

３－４ 波数ベクトルフィルタリングによる弾性波映像化の基本原理

ここまで超音波パルスドプラ法を用いてずり弾性波の瞬時変位を推定し、そこからずり 弾性波の伝播速度を推定するまでを説明してきたが、ここで空間領域と波数領域でずり弾 性波の伝播の様子を観察することはそれぞれ利点があるので下にその利点を示す。

空間領域

・ずり弾性波の伝播の様子が振幅、位相として直接観察できる。

・反射、屈折、回折の様子が可視化できる。

波数領域

・ずり弾性波の伝播方向、伝播速度の違いで分離できる。

・振幅の違いでマッピングできる。

という利点がある。そこでずり弾性波の複素変位を一度波数領域へ変換し、複数ずり弾性 波を分離した後、ある伝播方向、伝播速度のずり弾性波のみ取り出すようなフィルターを かけ複素変位空間へ戻すことで、ある伝播方向、伝播速度の波のみ可視化することができ る。反射波が存在するような場所やその強さなど直接動画像として可視化することができ る。

（１） 波数ベクトルフィルタリングの基本原理

今ある位置

P = ( , , ) X Z Y

^{での開口幅}dにおける波数スペクトラムは式(3-3-4)より

1

( , , )

exp{ ( ' ) }sin sin sin

2 2 2

v L

y ly

x lx z lz

l l

l

P f

k k

k k k k

j c d c d c d

δ

=

− ′

′ ′  

− −

   

= −      

     

∑

k P

k k P

で表される。ここである伝播速度を取り出すようなフィルター

F k P

_v

( , , ) f

^は

2 2

( , ) 1

^{v v}

v v low high

f f

F f ₌ π _≥ v v _≥ π

k k k

　　 　　　 　かつ

(3-4-1)

0 その他

ただし

: :

low high

v v

低域カット速度 高域カット速度

またある方向のみ取り出すようなフィルター

,

( , , ) F

_{θ φ} k P

f

1 1

,

2 2 2 2

1 1

( ) 1 tan ( ) tan ( )

tan ( ) tan ( )

x x

low high

y y

y x y x

low high

z z

k k

F k k

k k k k

k k

θ φ

φ φ

θ θ

− −

− −

= ≥ ≥

+ +

≥ ≥

k 　 　　　　　　 　かつ　　 　　

かつ 　かつ　　 　

0 その他 (3-4-2) ただし、

: : :

:

low high low high

φ φ θ θ

低域カット方位角 高域カット方位角 低域カット天位角 高域カット天位角

である。ここでフィルター後の波数スペクトラム

P ′ ( , , )

k P

f

^は

( , ,

_v

) ( , ,

_v

)

_v

( ,

_v

)

,

( )

P ′

k P

f = P

k P

f F

k

f F

_{θ φ} k (3-4-3)

となり、

P ′ ( , , )

k P

f

を逆フーリエ変換で空間領域へ戻すと、

/ / /

^

/ / /

( , ) ( , , ) exp( )

d d d

n v v x z y

d d d

f P f j dk dk dk

π π π

π π π

ξ

− − −

′ ^p = ∫ ∫ ∫ ′ ^{k P k p} ⋅

^(3-4-4)

となり、ある伝播方向、伝播速度を持つずり弾性波を取り出すことができる。

また開口がオーバーラップするようにの中心座標P

= ( , , ) x z y

^{を移動させていき、ある位}

置p

= ( , , ) x y z

^{での推定結果が}N個あるとすると、それらを平均処理することにより。

1

ˆ ( , ) 1 ˆ ( , )

N

v n v

n

f f

ξ N ξ

=

′

^p

= ∑ ′

^{p (3-4-5)}

となり、空間的にずり弾性波が伝播する様子が確認できる。

（２）波数ベクトルフィルタリングによるずり弾性波映像化

（１）で波数ベクトルフィルタリングにより特定の伝搬速度のずり弾性波を取り出すこ とが可能と述べた。しかし、速度の異なる２つの物体が存在し、一方の伝搬速度のみを取 り出すフィルターをかけた場合でも物体の境界を画像で判別することは難しい。これはず り弾性波が組織中で減衰してしまうためである。加振源に近い領域では波数スペクトラム の振幅が大きいため、フィルター範囲外の波動のサイドローブが大きく影響してしまう。

そこで新たな波数ベクトルフィルタリングを用いて特定のずり弾性波映像化を行う。

いま計測領域内に伝播速度が

1

,

2

v v

^{の異なる２つの組織}

O O

₁

,

₂が存在するとする。この時 速度は

v

1

> v

2^{とし、周波数}

f

_v^{で解析する。}

Fig.3-4-1 計測領域内の非一様物体モデル

今ある位置P

= ( , , ) X Z Y

^{での開口幅}dにおける波数スペクトラムは式(3-3-4)より

1

( , )

sin sin sin

2 2 2

L

y ly

x lx z lz

l l

P

k k

k k k k

c d c d c d

δ

=

− ′

′ ′  

− −

   

=      

     

∑

k P

で表される。

2

,

2

O v = v

1

,

1

O v = v

x

z

トランスデューサ

計測領域

d

d

組織

O

1中を伝播する波動のみを取り出すフィルター

( ,

_thr1

)

F

k

v

^は

1 1

( ,

_thr

) 1 2 f

^v _thr

F v ₌ π _≥ v

k k

　　 　　　 　

(3-4-6)

0 その他

となる。ただし

2 thr1 1

v < v < v

^(3-4-7)

ここでフィルター後の波数スペクトラム

'( , ,

_thr1

)

P

k P

v

^{は次式で表される。}

1 1

( , ,

_thr

) ( , ) ( ,

_thr

)

P ′

k P

v = P

k P

F

k

v

(3-4-8)

( , ,

_thr1

)

P ′

k P

v

を逆フーリエ変換し空間領域へ戻すと

/ / /

^

1 1

/ / /

( , ) ( , , ) exp( )

d d d

thr thr x z y

d d d

v P v j dk dk dk

π π π

π π π

ξ

− − −

′ ^p = ∫ ∫ ∫ ′ ^{k P k p} ⋅

^(3-4-9)

となる。

次に組織

O O

1

,

2中を伝搬する波動を取り出すフィルター

( ,

_thr2

)

F

k

v

^は

2 2

( ,

_thr

) 1 2 f

^v _thr

F v ₌ π _≥ v

k k

　　 　　　 　

(3-4-10)

0 その他

となる。ただし

2 2

v

thr

< v

^(3-4-11)

フィルター後の波数スペクトラム

'( , ,

_thr1

)

P

k P

v

^は

2 2

( , ,

_thr

) ( , ) ( ,

_thr

)

P ′

k P

v = P

k P

F

k

v

(3-4-12)

となり、これを空間領域へ戻すと次式で表される。

/ / /

^

2 2

/ / /

( , ) ( , , ) exp( )

d d d

thr thr x z y

d d d

v P v j dk dk dk

π π π

π π π

ξ

− − −

′ ^p = ∫ ∫ ∫ ′ ^{k P k p} ⋅

^(3-4-13)

組織

O

1^{中を伝搬する波動}

^

( , v

_thr1

)

ξ ^′

^{p を組織}

O O

1

,

2^{中を伝搬する波動}

^

( , v

_thr2

)

ξ ^′

^{p で規格化す}

ると、

^

^

1

^

2

( , ) ( )

( , )

thr s

thr

v v ξ ξ

ξ

′ = ′

′

p p

p

(3-4-14)

となる。

^ s

( )

ξ ^′

^{p は伝搬速度}

v

thr2^から

v

_thr1の成分を抑える事ができ、

1

v

thr ^{以上で伝搬する組織}

を強調するような画像を得る事ができる。

^

( , v

_thr2

)

ξ ^′

^p で規格化することにより減衰や複雑な 伝搬の波動がある場合でも測定が可能となる。

また、式(3-4-5)と同様に開口がオーバーラップするよう中心座標を移動させ、

^ s

( )

ξ ^′

^{p の移動}

平均を求めることで、空間的に

1

v

thr 以上で伝搬する組織を強調する画像を得る事ができる。

３－５ 信号処理フロー

（１） 信号処理のフロー

以下に提案手法の処理手順を示す。

（Ⅰ）ハードウェア

① 低周波のずり弾性波を組織内部に励起する

② 超音波パルスを送波し、RF信号を得る

③ 直交検波を行い、QI信号を計算する

④ A/D変換したQI信号をPCへ転送する

（Ⅱ）ソフトウェア

① QI信号をサンプリング周波数1GHzにアップサンプリングする

② アップサンプリングされた QI 信号にそれぞれ位相が90 度ずれたキャリア信号を乗 算

③ 上記の信号を和算し、RF信号を復元

④ RF相関により組織内の振動変位を計算

⑤ 振動変位を加振周波数帯のフーリエ変換を行い、複素変位を得る

⑥ 微小区間での２次元フーリエ変換を行う

⑦ 波数ベクトルフィルタリングによる映像化

ピークスペクトルからずり弾性波の伝搬方向、伝搬速度を推定

（２） RF信号の復元

変位測定を行うためQI信号からRF信号を復元する必要がある。

数MHzのQI信号を1GHzにアップサンプリングし、低域通過フィルターを用いて雑音を 低減する。アップサンプリングした後のQI信号は(3-1-10)式、(3-1-12)式より

4

0

( ) ( , )

( , ) cos( ')}

2

^u

f t

I t A t k z

c τ π ξ

τ = − −

4

0

( ) ( , )

( , ) sin{ '}

2

^u

f t

Q t A t k z

c τ π ξ

τ = − −

と表す。