変位推定による鍼置のその場可視化システム

令 和 二 年 度 修 士 論 文

変位推定による鍼位置のその場可視化システム

指導教員 山越 芳樹 教授

群馬大学大学院理工学府 理工学専攻

電子情報・数理教育プログラム

小川 智也

1 －目次－ 第 1 章 序章 3 第 2 章 せん断波計測について 5 2-1 せん断波とは 2-2 生体内組織における低周波振動の伝播 2-3 せん断波計測で期待されるパラメータと臨床的有用性 第 3 章 実時間せん断波映像法の原理 9 3-1 超音波パルスドプラ法による組織内振動伝播計測 3-2 カラーフロー映像系（CFI）の流速推定アルゴリズム 3-3 CFI の流速推定アルゴリズムによるせん断波波面検出 3-4 定量的なせん断波画像の構成法 3-5 波数ベクトルフィルタリング 第 4 章 測定機構について 24 4-1 タブレットエコー装置について 4-2 鍼と加振器の一体化仕様について 第 5 章 変位推定法の原理 28 5-1 鍼位置推定の基本的なアイディア 5-2 微分変位法について 5-3 Loupas 法について 5-4 Loupas 法での測定誤差改善 5-5 振幅図の s/n 値向上のためのフィルタ設計 第 6 章 変位振幅による鍼の位置可視化の基本原理 39 6-1 理論構成のための基本的なアイディア 6-2 長軸方向鍼位置推定法 6-3 短軸方向鍼位置推定法 第 7 章 シミュレーションによる鍼位置可視化システムの検証 43 7-1 長軸 d 推定シミュレーション 7-2 短軸 d 推定シミュレーション

2 第 8 章 ファントムを用いた鍼位置推定の評価 49 8-1 実験系の構成と実験方法 8-2 ファントムを用いた実験結果 8-3 実験結果への 3 次元位置推定の定量性評価 8-4 鶏むね肉を用いた実験結果 8-5 当研究室で提案した C-SWE 法との比較 第 9 章 結論 67 9-1 結論 9-2 今後の課題 謝辞・参考文献 68

3 第 1 章 序論 現在、痛み治療のためにファシア・ハイドロリリース治療のような鍼治療が注目されてい る。鍼治療では, ファシアなどの対象部位に正確に刺鍼することが求められる。そのた め、穿刺針と同様にエコーガイド下刺鍼が行われている。しかし、Fig.1-1 に示すように鍼 の直径は穿刺針と比較すると、0.18[mm]と非常に細く、B モード画像から鍼位置を見つけ 出すことは困難であり、 かつ鍼が超音波映像面に達したとしても鍼の可視化ができない という問題点が挙げられる。これらの問題を解決し, 安全, 正確に穿刺・刺鍼を行うには, 3 次元的に鍼の位置を把握し、変化する鍼位置をリアルタイムで追従する必要がある。 これまで、本研究室では C-SWE 法を用いて 3 次元に穿刺針の位置の可視化を行い, 穿 刺針の 3 次元位置可視化を行うシステムを実現してきた.C-SWE 法とは.連続せん断波映像 法とも呼び Fig.1-1 に示すように CFI（＝カラーフロー画像）や速度図,方向図を得る手法 である.しかし,この手法で作成される速度図と方向図は CFI を複数フレーム用いて作成さ れるため、作成に時間を要するため実時間映像化が困難である.そのため,昨年度には鍼自 体を加振した時に得られる CFI の１フレームの情報から,せん断波の波面や振動振幅に着目 して,3 次元的位置の可視化を実現させリアルタイム性良く簡便な測定システムを構築した. しかし,CFI は刺鍼中途の状態を観察すると鍼位置の視認性が悪い.つまり時間変化への応答 性が悪く,鍼位置を常に把握できない. そこで,本研究ではせん断波の振幅を１フレームごとに算出する変位振幅推定を用いるこ とで,実時間映像で鍼位置を常に把握できる新しい鍼位置可視化システムを構築する. また,タブレットエコー装置を導入し臨床利用を考慮した簡便な測定システムを新たに構築 し,シミュレーションおよび基礎実験でその有用性を確認する. Fig.1-1 穿刺針と鍼の比較・鍼挿入時の B モード画像例

4

5 第 2 章 せん断波計測について 本章ではせん断波の特徴と工学的な研究課題, 生体軟組織内部における低周波振動の伝 播について示す. さらに, せん断波計測により期待される臨床意義や目的について示す. 2-1 せん断波とは ここでは, せん断波の特徴とそこから考えられる工学的な課題について示す. せん断波の特徴 1. 波長 波長は数ミリメートルであるため, 高分解能測定が求められている. 2. 振幅 振幅は数十ミクロン以下であり, 高精度超音波計測技術が求められている. 3. 周波数 主にせん断波の減衰により制限され, 現在利用できるのは 100[Hz]～数[kHz]である. 工学的な研究課題 1. せん断波の波動としての性質 伝播方向が一様ではなく多重反射や回折, 減衰の問題がある. 2. せん断波の励起方法 振幅を得ようとすると加振器のサイズが大きく重くなる. また, 効率の問題もあり加振器 の発熱の問題がある. 3. パラメータ推定法, その物理, 臨床的意味づけ

6 2-2 生体内部の組織における低周波振動の伝播 生体組織の粘弾性パラメータと低周波振動の伝播速度および減衰の関係について以下に 示す. 外部から媒質に振動を伝えると, その振動は一般的に縦波・横波として伝播する. 生体の ような粘弾性媒質中では, Hooke の法則が成り立つ Voigt モデルと仮定することにより, 組 織の粘弾性を推定する手法が提案[1]されている. この縦波・横波の伝播速度および減衰係数 は次式で与えられる. ① 縦波 伝播速度 ： 𝑣_𝑙 = 𝜔𝑣 𝑅𝑒[𝑔] (2-2-1) 減衰係数 ： 𝛼𝑙= −𝐼𝑚[𝑔] (2-2-2) ただし，g = { 𝜌𝜔𝑣2 (2𝜇+λ)} 1 2 (2-2-3) ② 横波 伝播速度 ： 𝑣_𝑡 = 𝜔𝑣 𝑅𝑒[ℎ] (2-2-4) 減衰係数 ： 𝛼𝑡= −𝐼𝑚[ℎ] (2-2-5) ただし，ℎ = {𝜌𝜔𝑣2 𝜇 } 1 2 (2-2-6) 𝜇 = 𝜇1+ jω𝑣𝜇2 𝜆 = 𝜆1+ 𝑗𝜔𝑣𝜆2 𝜇1 : せん断弾性係数 𝜆1∶ 体積弾性係数 𝜇2 : せん断弾性係数 𝜆2∶ 体積弾性係数 ρ ∶ 密度 ω𝑣: 振動周波数 𝑅𝑒[ ]，𝐼𝑚[ ]：[]内の複素数の実数部，虚数部 また, これら縦波や横波の他に生体の表面付近では表面波が存在するが, この伝播速度は ほぼ横波の伝播速度に等しいことが知られている. 上記の波動の中で, 縦波は圧縮性の波 であり, 媒質を圧縮することにより伝播する. 一方, 横波は非圧縮性の波であり, 媒質を等 体積のまま, 横方向に挟み切るように変形させながら伝播していくため, せん断波とも呼 ばれている. ここで, 周波数が 1[kHz]程度以下の低周波振動であると, 外部から与えられ た振動のエネルギーはそのほとんどが横波に変換されると考えられている[2]. ここで, (2-2-4)式, (2-2-5)式で与えられる横波の伝播速度と減衰係数を, 粘弾性パラメータ

7 を用いて書くと, 𝑣𝑡= √ 2(𝜇12+𝜔𝑣2𝜇22) ρ(𝜇1+√𝜇12+𝜔𝑣2𝜇22) (2-2-7) 𝛼𝑡= √ 𝜌𝜔𝑣2(𝜇1+√𝜇12+𝜔𝑣2𝜇22) 2(𝜇₁2_+𝜔 𝑣2𝜇22) (2-2-8) となる. したがって, もし, 媒質の弾性が粘性にまさり, 𝜇1≫ ω𝑣𝜇1の関係が成り立つときには, 𝑣𝑡1≅ √ 𝜇1 𝜌 (2-2-9) 𝛼𝑡1≅ 0 (2-2-10) となり, 伝播速度は, 単にせん断弾性係数と媒質の密度のみの関数となる. このとき, 𝜇1が 大きいということは, 媒質が硬いということであり, 硬い媒質ほど伝播速度は速くなる. 一方, 媒質の粘性が弾性にまさり𝜇1≪ ω𝑣𝜇1の関係が成り立つときには, 𝑣𝑡2≅ √ 2𝜔𝑣𝜇2 𝜌 (2-2-11) 𝛼𝑡2≅ √ 𝜌𝜔𝑣 2𝜇2 (2-2-12) となり, 𝑣𝑡2・𝛼𝑡2とも粘性係数と密度の関数になり, この場合𝑣𝑡2・𝛼𝑡2の周波数依存性（分 散性）が現れてくる. Fig.2-2-1 に弾性体と粘弾性体の周波数別伝播速度を示す. Fig.2-2-1 弾性体と粘弾性体の周波数別伝播速度 0 1 2 3 4 5 6 100 200 300 400 500 600 700 800 900 1000 伝播速度 [m /s ec ] 加振周波数[Hz] 弾性率2.26kPa, 粘性率2.38Pa・s 弾性のみの場 合（2.26kPa）

8 2-3 せん断波計測で期待されるパラメータと臨床的有用性 せん断波の伝播速度は, 臨床的な有用性が明らかにされているが, せん断波計測によっ て得られる情報としては, この他にも Fig.2-3-1 に示すような情報も得られると考えられ る. 測定量 物理パラメータ 臨床意義 計測時の問題点 伝播速度 せん断弾性係数 組織の硬さ 多重反射，減衰 減衰係数 粘性係数 粘性評価 多重反射，屈折， 反射 伝播速度の 周波数依存性 粘性評価， 測定の定量性向上 多重反射，減衰， 空間分解能 共振現象 せん断弾性係数 組織のボリュームの 大きさ 減衰，空間分解能 非線形性 初期応力， 媒質の非線形性 組織非線形性評価 振動振幅の減衰 異方性 伝播速度の方向性 繊維方向，繊維化 三次元伝播方向 Fig.2-3-1 せん断弾性波によって得られる情報

9 第 3 章 実時間せん断波映像法について(C-SWE 法) 3-1 超音波パルスドプラ法による組織内振動伝播計測 組織内振動伝播計測は, 組織表面から振動を印加することで組織内に振動を励起させ, 内部を伝播する振動を超音波で計測するものである. これは, 組織内部を多数の超音波散 乱体と考えると, 組織内部に超音波を送波し, 超音波散乱体から反射してくる超音波がド ップラー効果によって周波数変調を受けていることに着目したものである. したがって, 超音波散乱体から反射した超音波を直交検波することで得られるドップラー信号から組織 内部を伝播する振動を推定することができる. 今, Fig.3-1-1 に示すような超音波トランスデューサに近づく方向に, 周波数𝑓𝑏, 速度v(𝑡) で振動する超音波散乱に対して超音波パルスを送波する場合を考える. Fig.3-1-1 計測モデル 散乱体の運動ξ(𝑡)は次式で表すことができる.

ξ(𝑡) = 𝜉

₀

𝑠𝑖𝑛(2π𝑓

_𝑏

𝑡 + 𝜙

_𝑏

)

(3-1-1) ただし, 𝜉0：振動振幅

𝜙

_𝑏：初期位相 この時, 超音波散乱体に反射した超音波の周波数 𝑓 は,

𝑓 =

𝑐+𝑣(𝑡) 𝑐

𝑓

0 (3-1-2) 𝑓0：超音波の中心周波数 𝑐 ：音速

10 この反射波が超音波トランスデューサで受信されるときの周波数𝑓′_は,

𝑓

′

=

𝑐 𝑐−𝑣(𝑡)

𝑓

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

𝑓

′

=

𝑐 𝑐−𝑣(𝑡)

×

𝑐+𝑣(𝑡) 𝑐

𝑓

0

=

𝑐+𝑣(𝑡) 𝑐−𝑣(𝑡)

𝑓

0

(3-1-4) したがって, 超音波のドプラ周波数シフト∆𝑓は,

∆𝑓 = 𝑓

′

− 𝑓

0

=

𝑐+𝑣(𝑡) 𝑐−𝑣(𝑡)

𝑓

0

− 𝑓

0

=

2𝑣(𝑡) 𝑐−𝑣(𝑡)

𝑓

0

(3-1-5) となる. 超音波ドプラ法で組織内の速度を観測する場合, 組織内での音速は約 1500[m/sec]であ り, それと比較して観測しようとする組織内の速度は 1~10 数[m/sec]と微小であるので, c ≫ v(𝑡)となり, (3-1-5)式は次式のように近似することができる.

∆𝑓 ≅

2𝑣(𝑡) 𝑐

𝑓

0 (3-1-6) この時, 超音波の位相変化∆𝜙は,

∆𝜙 = 2π ∫(∆𝑓)𝑑𝑡

=

4𝜋𝑓0 𝑐

∫ 𝑣(𝑡)𝑑𝑡

=

4𝜋𝑓0 𝑐

𝜉(𝑡)

(3-1-7) となるので, この散乱体からの受信信号𝑟(𝑡)は,

𝑟(𝑡) = 𝐴(𝑡)𝑠𝑖𝑛(2π𝑓

₀

𝑡 + ∆𝜙 − 2𝑘

_𝑢

𝑍)

= 𝐴(𝑡)𝑠𝑖𝑛 (2π𝑓

₀

𝑡 +

4𝜋𝑓0 𝑐

𝜉(𝑡) − 2𝑘

𝑢

𝑍)

= 𝐴(𝑡)𝑠𝑖𝑛 {2π𝑓

0

(𝑡 + 2

𝜉(𝑡) 𝑐

) − 2𝑘

𝑢

𝑍}

(3-1-8) ただし, 𝐴(𝑡)：振幅 𝑘𝑢 ：超音波パルスの波数 𝑍 ：トランスデューサ，散乱体間の距離

11 となる. よって超音波パルス間で微小変位𝜉(∆𝑡)による位相ずれが生じる. 次に RF 信号に, 位相が互いに 90 度異なる超音波周波数成分を畳み込み積分し低域通 過フィルターをかけ, IQ 信号を得る. (ⅰ) I 信号 RF 信号にキャリア信号を乗算すると,

𝐼

′

(𝑡) = 𝐴(𝑡)𝑠𝑖𝑛 {2π𝑓

0

(𝑡 + 2

𝜉(𝑡) 𝑐

) − 2𝑘

𝑢

𝑍} 𝑠𝑖𝑛(2π𝑓

0

)

= 𝐴(𝑡) 2 {𝑐𝑜𝑠 (4π𝑓0𝑡 + 4π𝑓0𝜉(𝑡) 𝑐 − 2𝑘𝑢𝑍) − 𝑐𝑜𝑠 ( 4π𝑓0𝜉(𝑡) 𝑐 − 2𝑘𝑢𝑍)}

(3-1-9) となる. ここで 2ω0付近の信号を低域通過フィルターで除くと,

𝐼(𝑡) =

𝐴(𝑡) 2

𝑐𝑜𝑠 (

4π𝑓0𝜉(𝑡) 𝑐

− 2𝑘

𝑢

𝑍)

(3-1-10) となり I 信号を得る. (ⅱ) Q信号 (ⅰ)と 90 度異なるキャリア信号を乗算すると,

𝑄

′

(𝑡) = 𝐴(𝑡)𝑠𝑖𝑛 {2π𝑓

0

(𝑡 + 2

𝜉(𝑡) 𝑐

) − 2𝑘

𝑢

𝑍} 𝑐𝑜𝑠(2π𝑓

0

)

= 𝐴(𝑡) 2 {𝑠𝑖𝑛 (4π𝑓0𝑡 + 4π𝑓0𝜉(𝑡) 𝑐 − 2𝑘𝑢𝑍) − 𝑠𝑖𝑛 ( 4π𝑓0𝜉(𝑡) 𝑐 − 2𝑘𝑢𝑍)} (3-1-11) となる. (ⅰ)と同様に低域通過フィルターを用いると,

𝑄(𝑡) =

𝐴(𝑡) 2

𝑠𝑖𝑛 (

4π𝑓0𝜉(𝑡) 𝑐

− 2𝑘

𝑢

𝑍)

(3-1-12) となり,Q 信号を得る.

12 3-2 カラーフロー映像系(CFI）の流速推定アルゴリズム 今, 超音波パルスを同一方向に N パルス送波すると, i 番目の超音波パルスに対する受信 超音波の位相𝜙𝑖は, 𝜙_𝑖= 𝜙0+ 2𝜋𝑓0 𝑐 2𝑣 𝑖 Δ𝑡 (3-2-1) ここで, 𝜙0: 初期位相 𝑓0: 超音波の中心周波数 𝑐: 音速 𝑣: 流速 Δ𝑡: 超音波パルス間の時間間隔 (3-2-1)式より, i 番目の受信 RF 信号𝑟𝑖は, 𝑟𝑖= 𝑟0 sin (2𝜋 𝑓0 𝑡 + 𝜙𝑖) = 𝑟0 sin (2𝜋 𝑓0 𝑡 + 𝜙0 + 2𝜋𝑓0 𝑐 2𝑣 𝑖 Δ𝑡) (3-2-2) この受信 RF 信号を直交検波器で直交検波すると, その複素直交検波出力𝑄⃗ 𝑖, および𝑄⃗ 𝑖の実 部信号および虚部信号である In phase 信号𝐼𝑖と, Quadrature 信号𝑄𝑖は, Q⃗⃗ _{𝑖= 𝐼𝑖+ 𝑗𝑄𝑖} 𝐼𝑖= 𝑎 cos (𝜙0+ 2𝜋𝑓0 𝑐 2𝑣 𝑖 Δ𝑡) (3-2-3) 𝑄𝑖= 𝑎 sin (𝜙0+ 2𝜋𝑓0 𝑐 2𝑣 𝑖 Δ𝑡) (3-2-3)式は, (3-2-4)式のように書くこともできる. Q⃗⃗ 𝑖 = 𝑎 𝑒𝑥𝑝( 𝑗(𝜙0+ 2𝜋𝑓₀ 𝑐 2𝑣 𝑖 Δ𝑡)) (3-2-4)

13 ここで, 第 i 番目の超音波パルスの位相と, 第 i+1 番目の超音波パルスの位相の差Δ𝜙𝑖を考 える. これは, Δ𝜙𝑖= 𝑎𝑟𝑔 (Q⃗⃗ 𝑖+1Q⃗⃗ 𝑖 ∗ ) (3-2-5) と推定できるので, (3-2-4)式を代入すると, Δ𝜙𝑖= 𝑎𝑟𝑔 ( 𝑎2 𝑒𝑥𝑝( 𝑗 2𝜋𝑓0 𝑐 2𝑣 Δ𝑡)) = 2𝜋𝑓0 𝑐 2𝑣 Δ𝑡 (3-2-6) よって流速𝑣は, 次式で求められる.

𝑣 =

𝑐 2𝜋𝑓0∙2Δ𝑡

Δ𝜙

𝑖

=

𝑐 2𝜋𝑓0∙2Δ𝑡

𝑎𝑟𝑔 (Q

⃗⃗

_𝑖+1

_Q

⃗⃗

_𝑖∗

₎

(3-2-7) (3-2-7)式のカッコ内は, IQ 信号を使うと, Q⃗⃗ _𝑖+1Q⃗⃗ _𝑖∗_{= (𝐼𝑖+1+ 𝑗𝑄𝑖+1) (𝐼𝑖+ 𝑗𝑄𝑖)}∗ = (𝐼_𝑖+1+ 𝑗𝑄𝑖+1) (𝐼𝑖− 𝑗𝑄𝑖) = 𝐼𝑖+1𝐼𝑖+ 𝑄𝑖+1𝑄𝑖+ 𝑗(𝐼𝑖𝑄𝑖+1− 𝐼𝑖+1𝑄𝑖) (3-2-8) と書けることより, 流速の推定式として,

𝑣 =

𝑐 2𝜋𝑓0∙2Δ𝑡

𝑎𝑟𝑐𝑡𝑎𝑛 (

𝐼𝑖𝑄𝑖+1−𝐼𝑖+1𝑄𝑖 𝐼𝑖+1𝐼𝑖+𝑄𝑖+1𝑄𝑖

)

(3-2-9) CFI では, S/N を向上させるために, 連続した超音波 N パルスから得た直交検波出力信号を 用いて以下の式で流速を推定している.

𝑣 =

𝑐 2𝜋𝑓0∙2Δ𝑡

𝑎𝑟𝑐𝑡𝑎𝑛 (

𝐸𝑈 𝐸𝐿

)

(3-2-10) 𝐸𝑈= ∑𝑁𝑖=1𝐼𝑖𝑄𝑖+1− 𝐼𝑖+1𝑄𝑖 𝐸𝐿= ∑ 𝐼𝑖+1𝐼𝑖+ 𝑄𝑖+1𝑄𝑖 𝑁 𝑖=1

14 3-3 CFI の流速推定アルゴリズムによるせん断波の波面検出 今, CFI の流速推定アルゴリズムをせん断波により反射体が正弦的に振動している場合に 適用する. せん断波が伝播して組織が正弦的に変動すると, 組織変位𝜉は次式のように表すことがで きる.

𝜉 = 𝜉

₀

sin (𝜔

_𝑏

𝑡 + 𝜙

₀

)

(3-3-1) 𝜔𝑏 : 振動角周波数 𝜙0 : 初期位相 このとき, i 番目の受信超音波パルスの位相𝜙𝑖は,

𝜙

𝑖

= 𝜙

0

+

2𝜋𝑓0 𝑐

2𝜉

(3-3-2) 直交検波器の出力は, (3-2-3)式と同様に, 𝐼𝑖= 𝑎 cos (𝜙0+ 2𝜋𝑓0 𝑐 2𝜉) (3-3-3) 𝑄𝑖 = 𝑎 sin (𝜙0+ 2𝜋𝑓0 𝑐 2𝜉) となる. ここでせん断波の角周波数に対して, 下記の条件(周波数条件)が成り立つ場合を考える.

𝜔

_𝑏

=

2𝜋 4Δ𝑡

(3-3-4) つまり, せん断波の周波数であらわすと, 𝑓𝑏= 1 4Δ𝑡 (3-3-5) さらに, 振動の初期位相として, 𝜙0= 0 (3-3-6) が満たされるとする. 上記条件（(3-3-5)式および(3-3-6)式）は, せん断波の伝播による組織の変位振動の周期 が超音波の４パルスに等しく, かつ初期位相が 0 の条件であり, これを変位振幅として図に

15 表すと Fig.3-3-1 にようになる. Fig.3-3-1 仮定した変位振幅 Fig.3-3-1 と同じ振動振幅は, せん断波の振動周波数が高く, エイリアジングにより低い周 波数に折り返す場合にも生じるが, この時の振動周波数は, ｍを整数として,

𝑓

_𝑏

=

1 2

(𝑚 +

1 2

)

1 Δ𝑡 (3-3-7) として表される. このため, 以下の議論は, (3-3-7)式が成り立つ場合にも成立するので,せ ん断波の周波数として(3-3-7)式が成り立てばよい（ CFI でせん断波を映像化するときの周波 数条件）. この時, 変位𝜉は, 𝜉 = 𝜉0 sin (2𝜋 𝑓𝑏 𝑖 Δ𝑡) (3-3-8) と表される. この時, 直交検波器の出力信号である I,Q 信号は, 𝐼𝑖= 𝑎 cos ( 4𝜋𝑓₀ 𝑐 𝜉) 𝑄𝑖= 𝑎 sin ( 4𝜋𝑓₀ 𝑐 𝜉) (3-3-9) となる. ここで, i=0,1,2,3 について, 直交検波器の出力を求めてみると,

16 i = 0の場合 {𝐼𝑖 = 𝑎 𝑄𝑖= 0 (3-3-10) i = 1の場合 𝐼𝑖= 𝑎 cos ( 4𝜋𝑓0 𝑐 𝜉0) (3-3-11) ただしλを超音波の波長とすると, ① 0 ≤ 𝜉0≤ 𝜆 8 の場合 { 𝐼𝑖 ≥ 0 𝑄𝑖≥ 0 (3-3-12) ② 𝜆 8 ≤ 𝜉0≤ 3𝜆 8 の場合 { 𝐼𝑖 ≤ 0 𝑄𝑖≥ 0 i = 2の場合 {𝐼𝑖 = 𝑎 𝑄𝑖= 0 (3-3-13) i = 3 の場合 𝐼𝑖= 𝑎 cos ( 4𝜋𝑓₀ 𝑐 𝜉0) (3-3-14) ただし, ① 0 ≤ 𝜉0≤ 𝜆 8 の場合 { 𝐼𝑖 ≥ 0 𝑄𝑖≤ 0 (3-3-15) ② 𝜆 8 ≤ 𝜉0≤ 3𝜆 8 の場合 { 𝐼𝑖 ≤ 0 𝑄𝑖≤ 0 となる. (3-3-11)-(3-3-15)式の関係をベクトル図であらわすと, ① 0 ≤ 𝜉0≤ 𝜆 8 の場合 Fig.3-3-2 に示すように,すべてのベクトルは第一象限と第四象限にある.

17 Fig.3-3-2 𝟎 ≤ 𝝃𝟎≤ 𝝀 𝟖 での直交検波器の出力信号 ③ 𝜆 8 ≤ 𝜉0≤ 3𝜆 8 の場合 ｉ=1 とｉ=3 の時のベクトルは第二象限と第三象限にある. Fig.3-3-3 𝝀 𝟖 ≤ 𝝃𝟎≤ 𝟑𝝀 𝟖 での直交検波器の出力信号

18 これらを Tab.3-3-1 にまとめる. Tab.3-3-1 直交検波器の出力信号 𝑖 𝐼𝑖 𝑄𝑖 0 𝑎 0 1 𝐼𝑎 * 𝑄𝑎 (正) 2 𝑎 0 3 𝐼𝑎 * −𝑄𝑎 (負) * 0 ≤ 𝜉0≤ 𝜆 8 のとき𝐼𝑎≥ 0、 𝜆 8 ≤ 𝜉0≤ 3𝜆 8 のとき𝐼𝑎≤ 0 次に, この IQ 信号のパターンに対して, CFI による速度推定値を求めてみる. まず(3-2-10)式で示される, 流速導出アルゴリズムは次の 2 つの基本演算からなる. Fig.3-3-4 流速導出の基本演算

19 ここで超音波パルスの送受信数 N=11 の場合に, CFI による流速導出アルゴリズムを図式化 すると, Fig.3-3-5 CFI における流速導出アルゴリズム (3-3-4), (3-3-6)式の 2 つの条件がともに満たされているとき, CFI における流速推定は Fig.3-3-6 のようになる.

20 Fig.3-3-7 CFI における流速導出アルゴリズムを使ったせん断波の波面再生 ここで, {0 ≤ 𝜉0≤ 𝜆 8 のとき 𝐸𝐿≥ 0 𝜆 8 ≤ 𝜉0≤ 3𝜆 8 のとき 𝐸𝐿 ≤ 0 となるが, ともに EU=0 であるので, 実軸を EL, 虚軸を EUとするベクトルは, ELが正の場合 は実軸上の正の方向を向くベクトルとなり, 流速推定値は 0 になる. 一方, ELが負の場合は 実軸上の負の方向を向くベクトルとなり, 流速推定値は正の最大値, または負の最大値（ナ イキスト周波数で決まる最大の流速値）になる. つまり, ① ELが正になる条件（せん断波による振動振幅が0 ≤ 𝜉0≤ 𝜆 8 の場合） 流速 0 になる. ② ELが負になる条件（せん断波による振動振幅が𝜆 8 ≤ 𝜉0≤ 3𝜆 8の場合） 振動振幅の位相が 0 度, および 180 度になる位置にＣＦＩ画像には流速最大の部分が 現れる. この条件は, せん断波の振幅により, せん断波による振動位相が０または 180 度の時に, 特 異なパターンが CFI 画像に現れることを示しており, これを振幅条件と呼ぶ.

21 せん断波が組織中を伝播しているとき, CFI 画像の中から上記に示したような特徴ある部 分を抽出することにより, せん断波の位相（0 度または 180 度）が推定できることになる. せん断波が等位相になる部分はせん断波の波面を再現することに相当するので, この方法 により, CFI 画像からせん断波の波面を再現できることになる. この方法は, 周波数条件（3-3-7 式）が成り立つときに, CFI の推定アルゴリズムが, せん 断波の 0 度と 180 度の位相を検出するディジタルフィルターになっていることに着目した, せん断波の映像化法である. 横軸を初期位相𝜙𝑏, 縦軸を振動振幅𝜉0として, 以下の条件で, 流速推定の数値シミュレーションをおこなった結果を Fig.3-3-7 に示す. [シミュレーション条件] 超音波中心周波数 𝑓0 6.5𝑀𝐻𝑧 超音波伝播速度 𝑐 _{1500 𝑚 𝑠}⁄ パルス繰り返し周波数 1 𝑑𝑡⁄ 365𝐻𝑧 パルス本数 𝑁 11 加振周波数 𝑓𝑏 91.25𝐻𝑧 Fig.3-3-7 数値シミュレーション結果 周波数条件は, 理論的には（3-3-7）式であらわされるが, 実際には, せん断波の周波数が この条件に近いときでも, 流速の最大値または流速０の部分が CFI 画像上に現れる. その ため, せん断波の周波数が周波数条件に近いときにも, せん断波の波面が再現できる. 𝜙𝑏

22 3-4 定量的なせん断波画像の構成法 周波数条件が成立する時, せん断波の位相 0 度および 180 度付近の 2 か所で流速推定値 𝐹𝑉𝑀が最大値または 0 の値を示す. そのため, カラードプラ像のフレームレートによって 見えるせん断波の偽りの角周波数を𝜔𝑎𝑙𝑖𝑎𝑠とすると, 𝐹𝑉𝑀は角周波数𝜔𝑝= 2𝜔𝑎𝑙𝑖𝑎𝑠の矩形波 となる. ここで, この矩形波の基本波のスペクトラム成分は, フーリエ変換を用いて 𝐹𝐹𝑉𝑀(𝑥, 𝑧, 𝜔𝑝) = ∫ 𝐹𝑉𝑀(𝑥, 𝑧, 𝑡) 𝑇𝐶𝐹𝐼 0 𝑒𝑥𝑝(𝑗𝜔𝑝𝑡)𝑑𝑡 (3-4-1) と表せる. また, その位相スペクトラム成分𝜃𝐹𝑉𝑀(𝑥, 𝑧)は, 以下のように導出される. 𝜃𝐹𝑉𝑀(𝑥, 𝑧) = 𝑎𝑟𝑔 (𝐹𝐹𝑉𝑀(𝑥, 𝑧, 𝜔𝑝)) (3-4-2) zx 平面を伝播する平面波の波数の x 成分𝑘𝑥と z 成分𝑘𝑧とすると, CFI で観測される波面の位 相は, せん断波の位相𝜙の二倍変化するため, 𝑥方向の単位長さあたりの超音波照射時間遅 れ∆𝑇𝑝を考慮すると, 𝑘𝑥(𝑥, 𝑧) = 1 2 𝜕 𝜃𝐹𝑉𝑀(𝑥,𝑧) 𝜕𝑥 + 𝜔𝑏∆𝑇𝑝 (3-4-3) 𝑘𝑧(𝑥, 𝑧) = 1 2 𝜕 𝜃𝐹𝑉𝑀(𝑥,𝑧) 𝜕𝑧 (3-4-4) また, |𝑘⃗ |は次式で表される. |𝑘⃗ | = √𝑘𝑥2+ 𝑘𝑧2= 2𝜋 𝜆 = 2𝜋 𝑣𝑏 𝑓𝑏 (3-4-5) よって, せん断波の伝播速度𝑣𝑏は 𝑣𝑏(𝑥, 𝑧) = 2𝜋𝑓𝑏 √𝑘_𝑥2_+𝑘 𝑧2 = 2𝜋𝑓𝑏 √(1 2 𝑑 𝜃𝐹𝑉𝑀(𝑥,𝑧) 𝑑𝑥 +𝜔𝑏∆𝑇𝑝)2+( 1 2 𝑑𝜃𝐹𝑉𝑀(𝑥,𝑧) 𝑑𝑧 )2 (3-4-6)

23 3-5 波数ベクトルフィルタリング 反射波により定在波が存在すると, せん断波の位相は空間的に変調されてしまう. せん 断波の複素振幅マップを二次元フーリエ変換して空間周波数上のスペクトルを計算し, 波 数ベクトルフィルタを適用してフーリエ逆変換することで任意の成分を抽出することがで きる. この方法を C-SWE 法で取得したせん断波位相マップに適用することを考える. 前方へ伝播する入射波と, 後方へ伝播する反射波が観測された一次元について考えると, せん断波の複素振幅𝑆(𝑥)は 𝑆(𝑥) = 𝐴𝐹𝑒𝑥𝑝[𝑗(𝑘𝑝𝑥 + 𝜑𝐹)] + 𝐴𝐵 𝑒𝑥𝑝[𝑗(−𝑘𝑝𝑥 + 𝜑𝐵)] (3-5-1) ただし, 𝐴𝐹:入射波の振幅 𝐴𝐵:反射波の振幅 𝜑𝐹:入射波初期位相 𝑘𝑝:せん断波の波数 𝜑𝐵:反射波初期位相である. 𝐴𝐵≪ 𝐴𝐹のとき, 入射波と𝑆(𝑥)との最大の位相差∆𝜃は ∆𝜃 = 𝑎𝑟𝑐𝑡𝑎𝑛 (𝐴𝐵 𝐴𝐹) (3-5-2) したがって, CFI から得られる複素信号𝑆𝐶𝐹𝐼(𝑥)は振幅情報が失われるので, 𝜑𝐹, 𝜑𝐵を無視 すると, 𝑆𝐶𝐹𝐼(𝑥) = 𝑒𝑥𝑝[𝑗 𝑎𝑟𝑔(𝑆(𝑥))] = 𝑒𝑥𝑝[𝑗{𝑘𝑝𝑥 + ∆𝜃 𝑠𝑖𝑛(−2𝑘𝑝𝑥)}] (3-5-3) 上式について, フーリエ級数展開すると 𝑆𝐶𝐹𝐼(𝑥) = ∑∞𝑛=−∞𝐽𝑛(∆𝜃)exp(𝑗𝑘𝑝𝑥) 𝑒𝑥𝑝[−2𝑗𝑛𝑘𝑝𝑥] (3-5-4) ここで, 𝐽𝑛(𝑥):n 次のベッセル関数である. 𝑘𝑝周りのスペクトラム成分のみを抽出するフィルタを適用することで, 入射波のせん断波 位相マップ𝜃𝐹𝑃𝑊は以下のように導出される. 𝑆𝐶𝐹𝐼′(𝑥) = 𝐽0(∆𝜃)𝑒𝑥𝑝(𝑗𝑘𝑝𝑥) (3-5-5) 𝜃𝐹𝑃𝑊(𝑥) = 𝑎𝑟𝑔(𝑆𝐶𝐹𝐼′(𝑥)) (3-5-6)

24 第 4 章 測定機構について 4-1 タブレットエコー装置について 今回使用したタブレットエコー装置を Fig.4-1-1 に示す. Fig.4-1-1 タブレットエコー装置本体 また,以下にタブレットエコー装置の性能について示します. サイズ：幅約 70×長さ 150×厚み 35[mm] 本体のみ 重量 ：300[g]以下(ケーブル込み) ケーブル長：1.5[m]

25 4-2 鍼と加振器の一体化仕様について (1) 使用する鍼と加振器 今回,実験に使用した鍼を Fig.4-2-1, 加振器について Fig.4-2-2 に示す.超小型リニア駆動ア クチュエータを用いた.なお、鍼の直径は 0.18[mm]で,鍼長は 40[mm]であるものを用いた. Fig.4-2-1 使用鍼 Fig.4-2-2 使用する超小型加振器（超小型リニア振動アクチュエータ） 本加振器の性能について, 以下に示す. ・共振周波数：170[Hz] ・駆動電圧：0.1～2.5[V] ・駆動電流：170[mA](定格) ・外形：円型（直径 10[mm]） : 厚さ：4.0[mm] コンパクトで簡便な加振方法で,鍼を加振させるためには,リニア振動アクチュエータを直 接鍼に当てて,加振することが有効です.

26 (2) 実験において再現性の高い機構 鍼位置可視化の実験を行うにあたって,鍼を再現性のある条件で刺鍼できるようにしたい. そのためには,鍼と超小型加振器を一体型となるように固定して鍼に加振を伝えることがで きる固定具が必要です.今回,3D プリンタを使用して PLA 樹脂製の加振器と鍼の固定具を制 作しました.その時の固定具の設計図を Fig.4-2-3,完成品を Fig.4-2-4 に示します。 Fig.4-2-3, 固定器具設計図 Fig.4-2-4, 固定器具完成品 鍼と超小型加振器,固定具を合体させたものを Fig.4-2-5 に示します。

27 Fig.4-2-5 鍼と加振器固定図 Fig.4-2-5 を実験では、基礎実験の際に用いる. (3)実際の治療に近い機構 (2)では,基礎実験のような定量性のある結果を求められるような加振機構であったが,実 際の痛み治療での現場では,医師が自らの手で患者に刺鍼を行っています.そのため,実際の 医療現場と同様に,指で鍼を直接持ち加振させる方法が必要になります. そのためには, Fig, 4-2-5(a)のようにゴム手袋内の親指に装着する.その後、超小型加振器 が装着されているゴム手袋と鍼の鍼柄を Fig.4-2-6 (b)のように接触させる. (a) 加振器装着の様子 (b) 鍼所持の様子 Fig. 4-2-6 超小型加振器装着方法 Fig. 4-2-6 に示すように,超小型加振器を装着することで,実際の治療に近い機構を再現した. この固定機構は,リアルタイム性実験において使用しました.

28 第 5 章 変位振幅推定法の原理 5-1 鍼位置推定の基本的なアイディア 鍼位置推定の基本的なアイディアについて Fig.5-1-1 を用いて説明します. Fig.5-1-1 鍼位置推定の基本的アイディア まず初めに,鍼に加振器を用いて微小振動を加えます.すると,鍼からせん断波が生体内部に 放射され伝播していきます.その後,発生したせん断波がタブレットエコー装置に到達する と,CFI を構成する IQ データを記録します.記録した IQ データから,微分変位法を用いてせ ん断波の振幅を出力します.これによって,せん断波の振幅の大きさから鍼の 3 次元位置を可 視化します.このアイディアのメリットとして,せん断波の振幅は時間変化への応答性が良 く鍼の 3 次元位置の可視化が可能です.また,振幅の大きさの変化で鍼位置を観察できる,つ まり視認性が高いこともメリットとして挙げられる.

29 5-2 微分変位法について 鍼位置推定の基本的アイディアに使用している「微分変位法」について説明します. 微分変位法とは,IQ 信号から生体変位(変位振幅)を推定し,クオリティとして評価する手法 です.この手法は,変位が正弦波である場合にその微分値も変位振幅に応じた大きさとなる ことを利用しています. 今、エコー装置から取得した IQ 信号が複素数である場合を考える. IQ 信号𝑄̃(t)は,次式で表すことができる.

𝑄̃ = Ae

jϕ(t) (5-2-1) ただし, A：振幅 ϕ(t)：位相 今,微小な時間∆ｔだけ離れて採取した 2 つの複素 IQ 信号から次式を考える.

∆𝑄

̃ =𝑄̃ (t + ∆ｔ) ∗ 𝑄̃(𝑡)

∗ _(5-2-2) ただし、

𝑄̃(𝑡)

∗：

𝑄̃(𝑡)

の位相共役 (4-2-2)式に(4-2-1)式を代入して

∆Q = Aexp (jϕ (t + ∆ｔ)) ∗ Aexp(−jϕ(t))

= 𝐴

2

_{exp (jϕ (t + ∆ｔ) − jϕ(t))}

_(5-2-3)

jϕ (t + ∆ｔ) − jϕ(t)

は微小時間の変化を表しているので

∆Q ≅ 𝐴

2

_{exp (𝑗∆ｔ}

dϕ(t)

𝑑𝑡

)

30 (5-2-3)式より

Δ𝑄

̃

の実数部と虚数部を用いて

∆ｔ

での

∆ϕ(t)

(変位の微分)を求める. ∆ϕ(t) = 𝑎𝑡𝑎𝑛 {𝐼𝑚(Δ𝑄̃ ) 𝑅𝑒(Δ𝑄̃ )} = 𝑎𝑡𝑎𝑛 { 𝐴2sin (∆ｔdϕ(t)_𝑑𝑡 ) 𝐴2_{cos (∆ｔ}dϕ(t) 𝑑𝑡 ) } ∆ϕ(t) = 𝑎𝑡𝑎𝑛 {tan (∆ｔdϕ(t) 𝑑𝑡 )} = ∆ｔ dϕ(t) 𝑑𝑡 (5-2-4) 今,IQ 信号の位相が正弦的な振動で引き起こされたとすると

ϕ(t) = 𝐾𝜉

₀

sin(𝜔

_𝑏

𝑡 + 𝜙

_𝑏

)

と表せる. (5-2-5) ただし,

_{K =}

4𝜋𝑓0 𝑐

𝑐

: 光速

𝑓

₀: 超音波の中心周波数

𝜉

₀ : 変位振幅

𝜔

_𝑏 : 変位周波数

𝜙

_𝑏 : 変位周波数 (5-2-4)式と(5-2-5)式より ∆ϕ(t) = ∆ｔdϕ(t) 𝑑𝑡 = ∆ｔ d 𝑑𝑡ϕ(t) = ∆ｔ d 𝑑𝑡{

𝐾𝜉

0

sin

(

𝜔

𝑏

𝑡 + 𝜙

𝑏)} = ∆ｔ

𝐾𝜉

₀ d 𝑑𝑡{

sin

(

𝜔

𝑏

𝑡 + 𝜙

𝑏)} =

𝐾𝜉

₀∆ｔ

𝜔

𝑏

cos

(

𝜔

𝑏

𝑡 + 𝜙

_𝑏) (5-2-6)

31 (5-2-6)式を振動周期 T の有限区間で加振周波数𝜔𝑏でフーリエ解析すると、

F(𝜔

𝑏

) = ∫

∆ϕ

(

t

)

𝑇 2 −𝑇₂

exp (−𝑗𝜔

𝑏

𝑡)𝑑𝑡

(5-2-7)

=

1

𝑇

∫

∆ϕ

(

t

)

𝑇 0

exp (−𝑗𝜔

𝑏

𝑡)𝑑𝑡

=

1

𝑇

∫ 𝐾𝜉

0∆ｔ

𝜔

𝑏

cos(𝜔

𝑏

𝑡 + 𝜙

𝑏

)

𝑇 0

exp (−𝑗𝜔

𝑏

𝑡)𝑑𝑡

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

∫ cos(𝜔

𝑏

𝑡 + 𝜙

𝑏

)

𝑇 0

exp (−𝑗𝜔

𝑏

𝑡)𝑑𝑡

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

∫ cos(𝜔

𝑏

𝑡 + 𝜙

𝑏

)

𝑇 0

exp (−𝑗𝜔

𝑏

𝑡)𝑑𝑡

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

1

2

∫ cos(2𝜔

𝑏

𝑡 + 𝜙

𝑏

) + cos 𝜙

𝑏 𝑇 0

− 𝑗 sin(2𝜔

𝑏

𝑡 + 𝜙

𝑏

) + 𝑗𝑡 sin 𝜙

𝑏

𝑑𝑡

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

1

2

[

1

2𝑤

_𝑏

sin(2𝜔

𝑏

𝑡 + 𝜙

𝑏

) + 𝑡 cos 𝜙

𝑏

+ 𝑗

1

2𝑤

_𝑏

cos(2𝜔

𝑏

𝑡 + 𝜙

𝑏

)

+ 𝑗𝑡 sin 𝜙

𝑏

]

𝑇

0

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

1

2

{(

1

2𝑤

_𝑏

sin(2𝜔

𝑏

𝑇 + 𝜙

𝑏

)

+ 𝑇 cos 𝜙

_𝑏

+ 𝑗

1

2𝑤

_𝑏

cos(2𝜔

𝑏

𝑇 + 𝜙

𝑏

) + 𝑗𝑇 sin 𝜙

𝑏

)

− (

1

2𝑤

_𝑏

sin 𝜙

𝑏

+ 𝑗

1

2𝑤

_𝑏

cos 𝜙

𝑏

)}

ここで,

𝜙

_𝑏

= 0

とすると

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

1

2

{

1

2𝑤

_𝑏

sin 2𝜔

𝑏

𝑇 + 𝑇 + 𝑗

1

2𝑤

_𝑏

cos 𝜙

𝑏

− 𝑗

1

2𝑤

_𝑏

}

32

=

1 𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏 1 2

{

1 2𝑤𝑏

sin 2𝜔

𝑏

𝑇 − 𝑗

1 𝑤𝑏

sin

2

_𝜔

𝑏

𝑇 + 𝑇}

(5-2-8) さらに,ここで

_{T =}

𝑛𝜋 𝜔_𝑏 (n=1,2,3…)とすると

F(𝜔

_𝑏

) =

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

1

2

(0 − 𝑗0 + 𝑇)

=

1

𝑇

𝐾𝜉

0∆ｔ

𝜔

𝑏

1

2

𝑇

F(𝜔

𝑏

)

=

1 2

𝐾𝜉

0

∆ｔ

𝜔

𝑏

(5-2-9) これによって,加振周波数

_𝜔

_𝑏におけるスペクトラムを求めた. (5-2-9)式より

F(𝜔

_𝑏

)

は複素数であるため,絶対値をとり, スペクトラム振幅を求める. この時の加振周波数のスペクトラム振幅は,加振振幅相当値である.

𝜉

₀

=

2 𝐾𝜉0∆ｔ𝜔𝑏

| F(𝜔

𝑏

)|

(5-2-10)

33 5-3 Loupas 法について 5-2 の「微分変位法」において,従来であれば基準信号と微小な時間（

∆ｔ

）だけ離れた シフト信号によって(5-2-4)式つまり“Arctangent 法”から変位の微分値

∆ϕ(t)

を求めてい ます.しかし.今回は”Loupas 法”を用いて変位の微分値

∆ϕ(t)

を求めています. それでは“Loupas 法”について説明していきます. そもそも、“Loupas 法”とは 2 次元自己相関として組織の動き（変位）を追従されるために 使用される計算方法です. ここで,“Loupas 法”の計算式について以下に示します.[2] U : 求める変位

𝑐

: 光速

𝑓

_𝑐: 超音波の中心周波数 m,n : それぞれの変位が計算される慈雨方向の範囲と時間的な範囲 分子(atan 部分)：微小時間方向の平均位相シフト[n 方向] 分母(atan 部分)：深さ方向の平均ドプラシフト[m 方向] “Loupas 法”は,分子部分の「微小時間方向の平均位相シフト」を分母部分の「深さ方向の平 均ドプラシフト」によって補正を行うことで変位を求めている. 変位が位相シフト推定問題 解いて定式化されている場合,自己相関アルゴリズムとして一般的に使用されています. そもそも微分変位法に使用されていた“Arctngent 法”や”Looupas 法”の分子のみの値,1 次元 自己相関で知られる“河西法”はすべて,「中心周波数に対する平均位相シフトの値」を表し ています.この平均位相のシフトの計算は,IQ 信号の基となる RF 信号の中心周波数（超音波 の中心周波数）が常に一定であると仮定して計算が行われています.しかし,実際には対象内 部の散乱体の振動によって中心周波数は変動し、正確な変位を求めることができません.そ のため,分母部分の平均ドプラシフトによる補正をすることで,より正確に変位を求められ ます.

34 メリット (1)超音波の中心周波数の局所的な変化に対して S/N の向上が見込め,変位の精度が高い (2)変位推定値は軸方向の寸法に沿って滑らかに変化する. (3)アルゴリズムの計算速度が速い (4)正規化された相互相関と同様の性能を発揮する デメリット (1)微小時間方向の広がりが小さく,低 S/N であるいくつかの場合においては,アルゴリズム の分母が物理的でないほど小さくなるため,変位を

−

𝜆 2

~

𝜆 2 に制限するフィルタが必要にな ります. (2)中心周波数からの変動が大きすぎる信号がある場合に,補正が上手くいかないため測地 誤差が生じやすい.

35 5-4 Loupas 法での測定誤差改善 前節では“Loupas 法”でのデメリットを挙げたが,この節では特にデメリット(2)の改善案 を発表します.改善案としまして,周波数の変動幅に補正可能範囲を設けて範囲を超えたデ ータを取り除く処理を行います. 短軸結果を例に,説明します. (a)短軸振幅図 (b)短軸ヒストグラム Fig.5-4-1 短軸の振幅図とヒストグラム Fig.5-4-1(a)のように,短軸の振幅図を表示して,黒丸内の鍼位置とみられる振幅の最大値を 中心に走査線 7 本分の値から,Fig.5-4-1(b)のようにヒストグラムを作成しました.このヒス トグラムから,補正可能範囲を-0.1~0.1 までとしました。 補正可能範囲を設定した場合の結果と設定しない場合の結果の比較を Fig.5-4-2 に示します.

36 (a)範囲の制限なし (b)-0.1 から 0.1 までの範囲 Fig.5-4-2 補正可能範囲の適応比較 Fig.5-4-2 に結果から補正可能範囲を設定することで中心周波数の大きすぎる変動によって 発生する振幅図内の特異なデータを除去することができました.これによって,振幅図の S/N が向上していることが分かります.

37 5-5 振幅図の s/n 値向上のためのフィルタ設計 鍼位置の 3 次元可視化を目的とする上で,S/N 値の向上は必須です. そのため,前述の周波数変動幅の補正可能範囲の設定に追加して,“Noise_Reduction Filter” を導入しました.一般的なノイズ除去手段として Median フィルタが存在します.しかし,この フィルタでは隣り合う値の差が大きいほど,高い値を補正してしまうため,ノイズ以外の値 も一緒に補正をしてしまう欠点を持ちます. そのため,“Noise_Reduction Filter”では, 特定の点の振幅とその上下を比較し,差が大きい時 のみ振幅を修正する処理を行う事で,インパルス状のノイズを除去する働きを持たせること で,振幅のノイズの特徴である急峻な変化のみを補正することができます.

ここで, “Noise_Reduction Filter”の働きを Fig.5-5-1 とフィルタの効果を Fig.5-5-2 に示 します.

Fig.5-5-1 “Noise_Reduction Filter”の働き

X2、X3 に値が周りに比べて急峻に高いノイズが存在しています.この高いノイズに“pct”(設 定値)を乗算し,その値が X1、X4 の値と比較し大きかった時のみ比較した値に補正します. “pct”の設定は,シミュレーション結果より算出し pct = 0.6 となっています.

38

Fig.5-5-2 フィルタの効果

フィルタを適応する前と後では,明らかにノイズが提言されていることが分かる.

これによって,この“Noise_Reduction Filter”は振幅図の S/N 値向上に非常に有用であるこ とが分かる.

39 第 6 章 変位振幅による鍼の位置可視化の基本原理 本研究では,超音波映像面に対して並行方向から刺鍼を行うことを「長軸刺鍼」,直交方向か ら刺鍼を行うことを「短軸刺鍼」と呼称します. 6-1 理論構成のための基本的なアイディア 鍼の位置可視化には,「せん断波振幅」に注目します.せん断波振幅は,波の伝播によって減 衰していきます.鍼と超音波映像面との距離ｄが大きくなると,その値に応じてせん断波振 幅は減少していく性質を利用します.これによって,振幅の大きさの変化で鍼位置を観察で きるため視認性が高い.またその変化は時間変化への応答性が良いため,鍼が刺鍼の途中で あっても正確に鍼位置を推定できると考えています. 6-2 長軸方向鍼位置推定法 鍼は超音波映像面に対して並行方向から刺鍼されるとする.鍼自体を加振して,対象に刺 鍼を行うと,せん断波は鍼から円筒波として放射されます,また,せん断波振幅の最大値は鍼 の形に沿って現れます. 鍼と超音波映像面との距離ｄの値が大きくなるに従って振幅の値は小さくなる.そのため, 振幅の値が大きい位置だけを取り出して表示できれば,その図は鍼の位置を示すものになる. 長軸での鍼位置推定の流れについて,Fig.6-2-1 に示す. Fig.6-2-1 長軸における鍼位置推定の流れ

40 長軸での鍼位置推定では,まず振幅図の最大値を見つけます. その後、最大値に対してしきい値を決めてしきい値を以上の振幅値を持つ部分だけを抽出 して二値化図を作成する.またこの時に,しきい値以上の振幅値の平均振幅を算出します. つまり,振幅の二値化図では,その先端部分から鍼先端部を推定できます.また,平均振幅から 鍼と超音波映像面との距離ｄを推定できます. これら二つの推定から鍼の 3 次元位置可視化を行います.

41 6-3 短軸方向鍼位置推定法 鍼は超音波映像面に対して直交方向から刺鍼されるとする.鍼自体を加振して,対象に刺 鍼を行うと,せん断波は鍼先から球面波として放射されます,また,せん断波振幅の最大値は 鍼先端の一点に現れます. 鍼と超音波映像面との距離ｄの値が大きくなるに従って振幅の値は小さくなります.その ため,振幅の最大値は常に鍼先端がその時の最大値であるため,振幅の最大値を探索すれば 鍼位置推定は可能となる. 短軸での鍼位置推定の流れについて,Fig.6-3-1 に示す. Fig.6-3-1 短軸における鍼位置推定の流れ

42 短軸での鍼位置推定も,まず振幅図の最大値を見つけます. 次にその値と位置を記録して,記録した座標を中心に特定の範囲をとる. その範囲内の平均振幅を算出します. つまり,鍼先端が振幅の最大値であることを利用して、範囲内に必ず鍼先端が存在する.これ によって鍼先先端部を推定できます.また,変位内の平均振幅から鍼と超音波映像面との距 離ｄを推定できます. これら二つの推定から鍼の 3 次元位置可視化を行います. 「特定の範囲の設定について」 今回,使用した範囲は縦方向に 0.7[mm],横方向に 0.12[mm]となりました. なぜこのような範囲になったのかというと,タブレットエコー装置の解析におけるデータ間 隔が関係しています. タブレットエコー装置の解析におけるデータ間隔 縦方向は 0.0966[mm],横方向は 0.597[mm] 縦方向のデータ点を 7 点とると,縦横同程度の長さと長さになるが,横方向のデータが 1 点のみとなり,範囲を作ることができないため中心から前後１点ずつデータを取得した.

43 第 7 章 シミュレーションによる鍼位置可視化システムの検証 本章では,第 6 章で述べた鍼位置可視化の基本原理について,シミュレーションを行った結 果を示す. 7-1 長軸 d 推定シミュレーション (1)シミュレーション概要 本シミュレーションは,鍼からせん断波が媒質内部に伝播し,これにより鍼と超音波映像 面との距離ｄを変化させたとき,せん断波の振幅図はどのようになるのか,理論に沿ってい るのかを明らかにすることを目的としている. (2)シミュレーション方法 シミュレーションモデルとして, 針長を 20[mm]とし, 針上の各々の点を波源とする. そ の各点から加振周波数の正弦波が球面的に伝播するとする. 任意の位置で各々の波源から の複素振幅の和を求めて図示します. 以下に本シミュレーションのモデル条件を述べる. ・加振周波数：156.0[Hz] ・鍼の挿入角度：35[deg] ・媒質の伝播速度：2.5[m/s] 以上の条件の時,鍼と超音波映像面までの距離ｄを 0.0～3.0[mm]の範囲で 1.0[mm]の間隔 で変化させたときのせん断波の振幅図,およびしきい値を用いた二値化図をシミュレーショ ンしていく. (3)シミュレーション結果 まず, 鍼と超音波映像面までの距離ｄを 0.0～3.0[mm]の範囲での振幅図のシミュレーシ ョン結果を Fig.7-1-1 に示し,鍼位置を追加した結果を Fig.7-1-2 に示す. Fig.7-1-1 振幅図シミュレーション結果(長軸)

44 Fig.7-1-2 振幅図シミュレーション結果(長軸)鍼位置追加 Fig.7-1-1 や Fig.7-1-2 の結果から d が大きくなるにつれて振幅の大きさが小さくなってい ることが分かります.このことから,振幅の大きさの変化は視認性が良く鍼位置を観察でき ます. 続いて,振幅図にしきい値を加えその値以上の振幅値を持つ部分だけを抽出し二値化したも のを Fig.7-1-3 に示します.また,Fig.7-1-2 と同様に鍼位置を追加した結果を Fig.7-1-4 に示 します.このときのしきい値は振幅最大値の 45%となっています.このしきい値を設定した 理由としては,振幅での鍼位置推定のシミュレーション条件として鍼と超音波映像面までの 距離ｄを 0.0～3.0[mm]の範囲で解析行うため, ｄ=3.0[mm]において最大値から鍼位置の 振幅が消えない最大の値の比率をしきい値として設定しました.

45 Fig.7-1-4 振幅二値化図シミュレーション結果(長軸)鍼位置追加 Fig.7-1-3, Fig.7-1-4 の結果から,振幅図の二値化図からどのｄであっても鍼先の位置を可視 化することができます.また,二値化図においてｄが大きくなると二値化図の範囲が大きく なっていることが分かります.これを数値化することで詳細な d 推定が可能であることが示 唆されます. 最後にしきい値以上の振幅値の平均振幅をシミュレーションより算出します. 各ｄにおける平均振幅を算出して,相対値として図示すると Fig.7-1-5 になります. Fig.7-1-5 各 d における平均振幅(相対値)のシミュレーション結果(長軸) Fig.7-1-5 の結果から、鍼と超音波映像面との距離ｄが大きくなると振幅の値が小さくなる ことが確認できた.このため,平均振幅を求めることでｄ推定が可能であることがシミュレ ーションから明らかになった. 上記のことから,長軸ではしきい値による振幅図の二値化から鍼先端部を推定し,しきい 値以上の振幅値の平均振幅から d 推定を行うことで鍼の 3 次元位置可視化を実現できる.

46 7-2 短軸 d 推定シミュレーション (1)シミュレーション概要 本シミュレーションは,鍼からせん断波が媒質内部に伝播し,これにより鍼と超音波映像 面との距離ｄを変化させたとき,せん断波の振幅図はどのようになるのか,理論に沿ってい るのかを明らかにすることを目的としている. (2)シミュレーション方法 本シミュレーションのモデル条件については 7-1 節のものと同等のものを設定する. 本節 では鍼と超音波映像面までの距離ｄを 0.0～2.0[mm]の範囲で 0.5[mm]の間隔で変化させた (3)シミュレーション結果 まず, 鍼と超音波映像面までの距離ｄを 0.0～3.0[mm]の範囲での振幅図のシミュレーシ ョン結果を Fig.7-2-1 に示し,鍼位置を追加した結果を Fig.7-2-2 に示す. Fig.7-2-1 振幅図シミュレーション結果(短軸) Fig.7-2-2 振幅図シミュレーション結果(短軸)鍼位置追加 Fig.7-2-1 や Fig.7-2-2 の結果から d が大きくなるにつれて振幅の大きさが小さくなってい ることが分かります.このことから,長軸シミュレーション結果同様に振幅の大きさの変化 は視認性が良く鍼位置を観察できます. 続いて,振幅図内の最大振幅を探査してその位置を中心に縦 0.7[mm],横 1.2[mm]の範囲 を作成して,振幅図に追加する. Fig.7-2-3 に設定範囲を追加した結果と Fig.7-2-4 には,設定 範囲だけではなく鍼位置を追加した結果を示します.

47 Fig.7-2-3 振幅図シミュレーション結果(短軸)設定範囲追加 Fig.7-2-4 振幅図シミュレーション結果(短軸)設定範囲と鍼位置追加 Fig.7-2-3 と Fig.7-2-4 の結果からどのｄであっても設定範囲が鍼先の位置とほぼ同一であ ることから鍼位置可視化することができます.また,設定範囲内においてｄが大きくなると 範囲内の振幅が小さくなっていることが分かります.これを数値化することで詳細な d 推定 が可能であることが示唆されます.

48 最後に設定範囲内の平均振幅をシミュレーションより算出します. 各ｄにおける平均振幅を算出して,相対値として図示すると Fig.7-2-5 になります. Fig.7-2-5 各 d における平均振幅(相対値)のシミュレーション結果(短軸) Fig.7-2-5 の結果から,鍼と超音波映像面との距離ｄが大きくなると振幅の値が小さくなる ことが確認できた.このため,平均振幅を求めることでｄ推定が可能であることがシミュレ ーションから明らかになった. 上記のことから,短軸では設定範囲から鍼先端部を推定し,設定範囲内の平均振幅から d 推定 を行うことで鍼の 3 次元位置可視化を実現できる.

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.5

1.0

1.5

2.0

平均振幅

(相対値

)

設定距離d[mm]

49 第 8 章 ファントムを用いた鍼位置推定の評価 8-1 実験系の構成と実験方法 (1)実験概要 前章までにおいて, 理論及びシミュレーションによって, 取得した振幅図を二値化や範 囲設定より高精度に鍼位置推定ができることが明らかとなった. これが実際に成り立って いるのかを実験を行い検証する. (2)実験系の構成 Fig8-1-1 に実験系の構成図を示す. Fig8-1-1 実験系の構成

50 今回の実験では,痛み治療に用いられる鍼を刺鍼することとして,使用する鍼を Fig8-1-2 に 示します.鍼の直径は 0.18[mm]で鍼長は 40[mm]であるものを使用した. Fig8-1-2 使用鍼 また,今回 Fig8-1-3 に示す超小型加振器(超小型リニア駆動アクチュエータ)を使用した. 本加振器の性能について, 以下に示す. ・共振周波数：170[Hz] ・駆動電圧：0.1～2.5[V] ・駆動電流：170[mA](定格) ・外形：円型（直径 10[cm]） Fig8-1-3 使用加振器(超小型リニア振動アクチュエータ) 以下に実験条件を述べる. ・超音波映像装置 タブレットエコー装置 ・加振周波数 156.0[Hz] ・加振器の振動振幅 33[μm] ・せん断波の伝播速度 2.6[m/s] ・加振対象 アルギン酸ファントム[鍼位置推定評価実験] こんにゃく[リアルタイム性実験]

51 (3)実験方法 実験方法について今回は, [リアルタイム性実験]と[鍼位置推定評価実験]を行いました. 長軸刺鍼の場合には,超音波映像面に対して並行方向から刺鍼を行う. Fig.8-1-4(a) 短軸刺鍼の場合には,超音波映像面に対して直交方向から刺鍼を行う Fig.8-1-4(b) [リアルタイム性実験] 1.超小型加振器を(Fig.8-1-6)のようにをゴム手袋内に装着させる. 2.指で鍼をつまみ加振器と鍼を接触させる 3.鍼と超音波映像面との距離ｄを (長軸)1.0[mm]ずつ変化させながら振幅図を取得する (短軸)0.5[mm]ずつ変化させながら振幅図を取得する [鍼位置推定評価実験] 1.鍼と超小型加振器を Fig.8-1-5 に示すように固定具によって一体化させる 2.鍼と超音波映像面との距離ｄを (長軸)1.0[mm]ずつ変化させながら振幅図を取得する. (短軸)0.5[mm]ずつ変化させながら振幅図を取得する. (a) 長軸刺鍼 (b) 短軸刺鍼 Fig.8-1-4 長軸刺鍼と短軸刺鍼の手法

52 Fig.8-1-5 鍼と超小型加振器固定の様子 Fig.8-1-6 超小型加振器を装着した様子 Fig.8-1-7(a)に[鍼位置推定評価実験]の様子を示し Fig.8-1-7(b)に[リアルタイム性実験]の様子を示した. (a) リアルタイム性実験 (b) 鍼位置推定評価実験 Fig.8-1-7 実験の様子

53 8-2 ファントムを用いた実験結果 [リアルタイム性実験] (1)長軸実験結果 こんにゃくを用いて長軸刺鍼を行い,得られた結果を Fig.8-2-1 に示す. なお d = 0.0[mm]の時の結果である. Fig.8-2-1 リアルタイム計測(長軸結果) Fig.8-2-1 の結果から,鍼の先端と推定される部分の振幅は大きく,かつ鍼が刺鍼している最 中にも振幅の大きい部分は鍼先を追従するように変化していることが分かる. これによって, 長軸での振幅の大きさの変化は視認性が良いことが確認できます.

54 (2)短軸実験結果 こんにゃくを用いて短軸刺鍼を行い,得られた結果を Fig.8-2-2 に示す. なお d = 0.0[mm]の時の結果である. Fig.8-2-2 リアルタイム計測(短軸結果) Fig.8-2-2 の結果から,鍼の先端と推定される部分の振幅は鍼と超音波映像面との距離ｄが 小さくなると振幅が大きくなり,かつ振幅が大きくなる過程を視覚的にとらえやすい結果と なった.これによって, 短軸での振幅の大きさの変化は視認性が良いことが確認できます. リアルタイム性実験を通じて,長軸,短軸共に振幅の変化への視認性が高く,振幅の変化は時 間変化への応答性が良いことが明らかとなった.

55 [鍼位置推定評価実験] (1)長軸実験結果 ファントムを用いて,長軸刺鍼にて測定を行った後, 得られた測定結果を示す. なお, d の 値を 0.0～3.0[mm]の範囲で 1.0[mm]ずつ変化させた際の結果を Fig.8-2-3 に示す. Fig.8-2-3 鍼測定結果（長軸） Fig.8-2-3 を見ると,鍼と超音波映像面までの距離ｄが大きくなると鍼位置を思われる位置 の振幅が小さくなっていることが分かります。この振幅図の結果は,シミュレーション結果 と同様であり振幅の大きさの変化は視認性が良いことが確認でいます.ここで,しきい値以 上の振幅値を持つ部分だけを二値化したものを Fig.8-2-4 に示した. Fig.8-2-4 鍼測定結果二値化図(長軸) 二値化図において,ｄが大きくなると二値化された時の領域が広がっていることが分かりま す. ここで,しきい値以上の振幅の平均値を求めます. 各ｄにおける平均振幅を算出して,相対値として図示すると Fig.8-2-5 になります.

56

Fig.8-2-5 各ｄにおける平均振幅（長軸）

Fig.8-2-5 の結果から,シミュレーション結果と同様に,鍼と超音波映像面との距離ｄが大き くなると振幅の値が小さくなることが確認できた.このため,平均振幅を求めることでｄ推 定が可能であることが明らかになった.

57 (2)短軸実験結果 ファントムを用いて,短軸刺鍼にて測定を行った後, 得られた測定結果を示す. なお, d の 値を 0.0～2.0[mm]の範囲で 0.5[mm]ずつ変化させた際の結果を Fig.8-2-6 に示す. Fig.8-2-6 鍼測定結果（短軸） Fig.8-2-6 を見ると,鍼と超音波映像面までの距離ｄが大きくなると鍼位置を思われる位置 の振幅が小さくなっていることが分かります。シミュレーション結果と同様であり振幅の 大きさの変化は視認性が良いことが確認でいます. ここで,最大値を探査してその位置を中心に縦 0.68[mm],横 1.21[mm]の範囲を Fig.8-2-7 に表示します.さらに各ｄにおける範囲内の振幅の平均値を求め相対値として図示すると Fig.8-2-8 となる.

58

Fig.8-2-7 鍼測定結果範囲表示（短軸）

59

Fig.8-2-8 の結果から,シミュレーション結果と同様に,鍼と超音波映像面との距離ｄが大き くなると振幅の値が小さくなることが確認できた.このため,平均振幅を求めることでｄ推 定が可能であることが明らかになった.

60 8-3 実験結果への３次元位置推定の定量性評価 (1) 長軸実験結果 ここで,長軸での d=0.0[mm]の時の二値化図の鍼位置精度を調べます. 今回の結果を Fig.8-3-1 に示し,別のデータでも精度を求めた時の結果を Fig.8-3-2 ます. Fig.8-3-1 今回の結果による二値化図の鍼位置精度(d=0.0[mm])

61 Fig.8-3-2 異なる結果での二値化図の鍼位置精度(d=0.0[mm]) Fig.8-3-1 と Fig.8-3-2 はともに白い破線が,実験結果から幾何学的に求めた鍼先先端位置に なります.また,二値化画像の鍼先端位置は図内の最も右下に存在する赤い部分の右端とし た. Fig.8-3-1 での二値化図と幾何学的に鍼位置を求めた時の誤差は,0.12[mm]となった. また, Fig.8-3-2 での二値化図と幾何学的に鍼位置を求めた時の誤差は,0.27 [mm]となった. よって,最大値のしきい値以上を持つ振幅を抽出して二値化図を作成することで,鍼先の位 置を高精度に推定することができることが明らかになった. また,Fig.8-2-5 の長軸の各ｄにおける平均振幅（相対値）の推移と, Fig.7-1-5 の長軸のシミ ュレーション結果を比較するために同じグラフにまとめたものを Fig.8-3-3 として図示する と以下のようになります.

62 Fig.8-3-3 長軸 d 推定のシミュレーションとの比較 Fig.8-3-3 の結果を見ると,実験結果はシミュレーション結果とほぼ同程度の推移をするこ とが分かった.最大誤差はｄ＝3.0[mm]のときで,誤差が 0.31[mm]となった. このことから,最大値のしきい値でその値を超える振幅値の平均振幅相対値を求めると,ｄ 推定が可能であることが明らかとなった. Fig.8-3-2 の鍼先の誤差精度結果と,リアルタイム性実験の結果とを合わせることで, 長軸においてリアルタイム性の高い高精度な鍼位置可視化が実現できることが分かる.

63 (2) 短軸実験結果 Fig.8-2-8 の短軸の各ｄにおける範囲内の平均振幅(相対値)の推移と,Fig.7-2-5 の短軸シミ ュレーション結果を比較するためにグラフにまとめたものを Fig.8-3-4 として図示すると以 下のようになります. Fig.8-3-4 短軸 d 推定のシミュレーションとの比較 Fig.8-3-4 の結果を見ると, 実験結果はシミュレーション結果とほぼ同程度の推移をするこ とが分かった.最大誤差は,ｄ＝2.0[mm]の時であり,その誤差は 0.29[mm]であった. このことから,振幅図内の最大値を探査し, 縦 0.68[mm],横 1.21[mm]の範囲で平均振幅を 算出して,その相対値を用いてことでｄ推定が可能であることが明らかになった. よって,リアルタイム性実験の結果と合わせることで,短軸においてリアルタイム性の高い 高精度な鍼位置可視化が実現できることが分かる.

64 8-4 鶏むね肉を用いた実験結果 鶏むね肉(疑似生体)内部への刺鍼を行い,その時の振幅の伝播を調べ理論に沿った結果が 導出できるか検証する (1) 実験方法 実験条件や実験系は 8-1 節に従うものとする.実際に刺鍼している様子を Fig.8-4-1 に示す. Fig.8-4-1 鶏むね肉刺鍼の様子

65 (2)実験結果

実験を行って得た結果を Fig.8-4-3 及び Fig.8-4-4 に示す. Fig.8-4-3 は,長軸刺鍼 また Fig.8-4-4 は短軸刺鍼での実験結果です. Fig.8-4-3 鶏むね肉長軸結果 Fig.8-4-4 鶏むね肉短軸結果 長軸短軸共に B モード画像では,鍼の位置を推定することは困難であった. 長軸の(b)CFI や(c)二値化図では,鍼の位置を推定することができます. 短軸の(b)(c)の振幅図において振幅の最大値が一点になっているため,鍼先を推定すること ができます.

66 8-5 当研究室で提案した C-SWE 法との比較 最後に,本研究では変位振幅を用いた実時間映像法を提案しました. ここで,これまで鍼の 3 次元可視化で使用されてきた C-SWE 法と本研究について特徴をま とめたものを Fig.8-5-1 に示します. Fig.8-5-1 C-SWE 法と本研究の提案手法の特徴

67 第 9 章 結論 9-1 結論 せん断波振幅の性質を利用することで,振幅の情報のみで長軸・短軸共に鍼の 3 次元位置と 先端位置の推定が可能になった.長軸では,二値化画像から鍼位置を推定でき,その誤差は 0.12[mm]や 0.27[mm]と非常に制度が良い.短軸では縦 0.68[mm]横 1.21[mm]の範囲内に 必ず鍼先が存在するため長軸同様に位置精度が良い. 振幅の大きさを可視化することで,リアルタイム性の高い鍼の 3 次元位置可視化が可能とな った.長軸では,最大値のしきい値以上の振幅の平均振幅(相対値)を算出することで鍼と超音 波映像面との距離ｄを高精度で推定することができた .また,短軸では縦 0.68[mm]横 1.21[mm]の範囲内の平均振幅（相対値）を算出することで長軸同様に精度よくｄ推定が可 能となりました. これらによって,変位振幅推定法によって実時間映像で鍼位置を常に把握できる新しい鍼の 3 次元位置の可視化を実現した. 9-2 今後の課題 せん断波の振幅によって,鍼の 3 次元位置可視化が可能になったが, 臨床利用の実現を 想定し, 生体内部の複雑な構造であっても, 鍼位置可視化が行えるか検証を行う必要があ る. そのため, 今回鶏胸肉を対象に刺鍼したが, 他の対象にも刺鍼を行い, さらに検証を行 っていく必要がある.

68 謝辞 本研究を進めるにあたり, 終始適切なご指導をいただいた群馬大学大学院理工学府 山越 芳樹教授に深く感謝申し上げます. 日頃の測定においてご支援いただいた荻野毅技官に感 謝申し上げます. 研究を共にし, 日々の実験や解析にご協力いただいた学部 4 年竹澤颯人氏 に心より感謝いたします. 最後に, 研究室での学生生活においてお世話になりました山越 研究室の皆様に感謝の意を表します. 参考文献 [1]砂川和宏， 金井浩. "動脈壁組織性状診断を目的としたずり弾性波伝播の計測とずり粘 弾性推定の検討. " 超音波医学 33. 1 (2006): 65-74.

[2] Xuan Ding, Debaditya Dutta, Ahmed M.Mahmoud, Bryan Tillman, Steven A.Leers, Kang Kim 2015 January