C-SWE SWE法による高精度針位置可視化システム

令和元年度 修 士 論 文

C-SWE 法による高精度

針位置可視化システム

指導教員 山越 芳樹 教授

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

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

太田 聖人

1 C-SWE 法による高精度針位置可視化システム －目次－ ページ 第 1 章 序章 3 第 2 章 せん断波計測について 5 2-1 せん断波とは 5 2-2 生体内組織における低周波振動の伝播 6 2-3 せん断波計測で期待されるパラメータと臨床的有用性 8 第 3 章 実時間せん断波映像法の原理 9 3-1 超音波パルスドプラ法による組織内振動伝播計測 9 3-2 カラーフロー映像系（CFI）の流速推定アルゴリズム 12 3-3 CFI の流速推定アルゴリズムによるせん断波波面検出 14 3-4 定量的なせん断波画像の構成法 22 3-5 波数ベクトルフィルタリング 23 第４章 せん断波映像による高精度穿刺針・鍼の位置可視化の基本原理 24 4-1 長軸方向針位置推定法 24 4-2 短軸方向針位置推定法 28 4-3 Volterra Filter を用いた画像補正 32 4-4 CFI を用いた実時間映像法 34 第５章 シミュレーションによる針位置可視化システムの検証 36 5-1 速度・方向図のシミュレーション 36 5-2 d 推定シミュレーション 39 5-3 CFI を用いた d 推定シミュレーション 41 第 6 章 ファントムを用いた高精度針位置推定の評価 43 6-1 実験系の構成と実験方法 43 6-2 ファントムを用いた実験結果 46 6-3 実験結果への 3 次元位置推定の定量性評価 50 6-4 血管模擬したファントムの実験結果 54

2 第７章 CFI を用いた実時間映像法の評価 56 7-1 実験系の構成と実験方法 56 7-2 CFI を用いたリアルタイム実験結果 59 7-3 CFI を用いたリアルタイム 3 次元位置推定の定量性評価 60 7-4 鶏胸肉を用いた実験結果 63 第８章 結論 66 8-1 結論 66 8-2 今後の課題 67 謝辞・参考文献 68

3 第 1 章 序章 抗がん剤など薬剤の体内への連続投与, 経口摂取ができない患者への栄養補給など, 長い 期間の間, 穿刺針を中心静脈（主に上大静脈と下大静脈）に留置する中心静脈留置（中心 静脈カテーテル）が行われることが多いが, この手技に関して医療事故が少なからず発生 している. このため超音波ガイド下での穿刺（超音波ガイド法）が推奨されているが, 超 音波ガイド下であっても依然として医療事故が発生している. 超音波ガイド法は従来の解 剖学的な構造を推察して穿刺を行う方法（ランドマーク法）に比べ穿刺成功率が高く合併 症の発生が低いとされているが, 超音波ガイド下であっても多くの医療事故が発生してい る. この大きな理由として, 留置針を静脈内に正しく留置することが難しいことが挙げら れる. これは, ①留置針自体が非常に細く, 超音波画像の輝度が十分に得られない場合があ る, ②超音波画像は 2 次元画像であるので, 留置針の一部が超音波画像面に達するまで留置 針の位置がわからず盲目状態で穿刺を行わなければならない, という問題があるためであ る. また, 痛み治療のための鍼治療では, ファシアなどの対象部位に正確に鍼を刺すことが 求められる. そのため, 穿刺針と同様にエコーガイド下刺鍼が行われている. しかし, 鍼の 直径は 0.18[mm]程度と非常に細く, B モード像から鍼位置を見つけ出すことは困難であり, かつ鍼が超音波映像面に達しないと鍼の可視化ができないという問題点が挙げられる. こ こで, Fig.1-1-1 に実際の穿刺針や鍼を挿入したときの B モード像の例を示す. Fig.1-1-1(b) の鍼より, 鍼が超音波映像面であっても B モード像から鍼を観測することができないこと が確認できる. これらの問題を解決し, 安全, 正確に穿刺・刺鍼を行うには, 3 次元的に針の位置を把握 する必要がある. このようなナビゲーションシステムにおいても, 使える技術は 2 次元の 超音波映像を応用したものになるが, 2 次元の超音波映像においていかに 3 次元的な針の位 置を計測し可視化するか, 通常の超音波映像では超音波反射強度が低いために針が見えづ らいという課題に対して, いかに針をコントラスト高く明瞭に可視化するかが必要になる. また, 実際の臨床利用を踏まえると, 付加的装置なしで簡便に測定が行えるシステムにす ることが求められる. そこで本研究では針自体を加振すると針自体から機械的振動波であるせん断波が放射さ れ, せん断波は生体組織内を伝播する. せん断波の伝播速度の映像を超音波プローブから 送信される超音波で測定し, せん断波の伝播図（せん断波の波面の伝播を表す画像）やせ ん断波の速度図（せん断波の伝播速度の画像）, せん断波の方向図（せん断波の伝播方向 の画像）を再生することで, 針の位置の 3 次元的位置の可視化を行う. このとき, C-SWE 法を通じて, せん断波の伝播図や速度図, 方向図を得ることで, 3 次元に針の位置の可視化 を行い, 針の 3 次元位置可視化を行うシステムを実現する. さらに本研究では, リアルタイム性良く簡便な測定システムを構築する. 針自体を加振

4

したときに得られる CFI の 1 フレームの情報のみから, せん断波の波面や振動振幅に着目 して, 3 次元的位置の可視化を実現させる. また, 針を加振させる加振器を超小型のものを 導入し, 臨床利用を考慮した簡便な測定システムを新たに構築させる.

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 第４章 せん断波映像による穿刺針・鍼の映像化法の基本原理 本研究では, 穿刺針や鍼の映像法として, ・手法Ⅰ：せん断波の速度図を用いた針位置可視化及び 3 次元の高精度な位置推定を実現 した「せん断波速度・方向図を用いた高精度映像法」 ・手法Ⅱ：せん断波の波動伝播や振動振幅に着目しリアルタイム性の良い測定系を実現した 「Color Flow 画像を用いた実時間映像法」 の 2 つの手法を提案する. 以降, 手法Ⅰについて 4-1～4-3 節, 手法Ⅱについて 4-4 節にて基本原理を述べる. 4-1 長軸方向針位置推定法 穿刺針は超音波映像面に対する長軸方向に穿刺されるとする. 穿刺針自体を加振し, 穿 刺針から放射されるせん断波の伝播速度像には, 穿刺針が速度の高い部分から超音波映像 面に投影した穿刺針の 3 次元での位置が, 伝播速度の高い部分の幅から映像面と穿刺針と の距離が測定でき, 両者を併せることで穿刺針自体が超音波映像面内になくても, 超音波 映像面に対する穿刺針の 3 次元位置が測定できる. (1)モデル(長軸での穿刺針の挿入) 穿刺針と超音波映像面との位置関係を Fig. 4-2-1 に示す. 2 次元(x,y)平面を考え, 穿刺針 は図に対して原点に垂直方向に穿刺され, それから X 軸方向に距離 d だけ離れた位置に超 音波映像面があるとする. Fig. 4-2-1 長軸で挿入する穿刺針のモデル

25 (2)穿刺針の 3 次元位置推定の原理 Fig. 4-2-2 に定式化のための変数の定義を示す. Fig. 4-2-2 穿刺針と超音波映像面 Fig. 4-2-2 において, 超音波映像面上Δｙだけ離れた点 Q における距離差Δd は, ∆d = √𝑑2_{+ ∆𝑦}2_{− d ≅ d (1 +}1 2( ∆𝑦 𝑑) 2 ) − 𝑑 =1₂∆𝑦_𝑑2 (4-2-1) ここで d は穿刺針と超音波映像面までの距離である. 今, 半径𝑟0の穿刺針の中心から, せん断波が球面状に伝播しているとすると, 点 P に対す る点 Q のせん断波の位相差Φは, 𝛷 = 𝑘𝑇∆𝑑 = 𝑘𝑇 2 ∆𝑦2 𝑑 (4-2-2) ここで, 𝑘𝑇はせん断波の波数であり, 𝑘𝑇 = 2𝜋 𝜆𝑇= 2𝜋𝑓𝑏 𝑣𝑇 (4-2-3) また, 𝜆𝑇：せん断波の波長 𝑣𝑇：せん断波の伝播速度 𝑓𝑏 ：せん断波の周波数 である. 点 Q における位相差Φのために, 超音波の映像面上 y 軸方向には, 位相の変化による偽 りのせん断波の伝播が生じているように見える, この偽りのせん断波の伝播速度𝑣𝑃は, 𝑣𝑃= 2𝜋𝑓𝑏 𝜕𝛷 𝜕𝑦 = 2𝜋𝑓𝑏( 𝜕𝛷 𝜕𝑦) −1 _(4-2-4)

26 （4-2-2）式より, 𝜕𝛷 𝜕𝑦

=

𝑘𝑇∆𝑦 𝑑 であることを考慮すると, 𝑣𝑃= 2𝜋𝑓𝑏( 𝑑 𝑘𝑇∆𝑦) (4-2-5) （4-2-3）式を代入すると, 𝑣𝑃 = 2𝜋𝑓𝑏( 𝑑 ∆𝑦) 𝑣_𝑇 2𝜋𝑓𝑏= 𝑑 ∆𝑦𝑣𝑇 (4-2-6) （4-2-6）式より, 穿刺針と映像面との距離 d が大きいほど超音波の映像面の y 軸上に現 れる偽りの伝播速度𝑣𝑃は, 媒質の伝播速度𝑣𝑇に対して高くなる. 今, 超音波映像面上で, せん断波の伝播速度が𝑣𝑃 > 𝑣𝑇𝐻となる y 軸方向の幅 W を考え る. 𝑊 =２∆𝑦，𝑣𝑝= 𝑣𝑇𝐻とおくと, （4-2-6）式より, W = 2 𝑣𝑇 𝑣_𝑇𝐻𝑑 (4-2-7) 今, 穿刺針と超音波映像面までの距離 d は, 𝑑 =_2𝑣1 𝑇𝑣𝑇𝐻𝑊 (4-2-8) 穿刺針の半径は𝑟0であるので, 映像面から穿刺針の端面までの距離𝑑𝑒は, 𝑑𝑒 = 1 2𝑣_𝑇𝑣𝑇𝐻𝑊 − 𝑟0 (4-2-9) (3)精度を高める手法 Fig. 4-2-3 において, 複数の閾値𝑣𝑇𝐻を定め, 各々の𝑊𝑣_𝑇𝐻 を求めることで, 複数の𝑣𝑇𝐻∙ 𝑊か ら平均値を求め, （4-2-9）式より, 𝑑𝑐𝑎𝑙𝑐を求める. Fig. 4-2-3 複数の閾値による評価 𝑑𝑐𝑎𝑙𝑐 = ( 1 2𝑣𝑇) [𝑣̅̅̅̅̅̅̅̅̅̅̅̅̅] − 𝑟𝑇𝐻∙ 𝑊𝑣𝑇𝐻 0 (4-2-9)

27

Fig. 4-2-4 において, 複数点のデータを用い, 平均値を求め, 信頼性が高まる評価をするこ とが可能となる.

28 4-2 短軸方向針位置推定法 穿刺針は超音波映像面に対する短軸方向に穿刺されるとする. 穿刺針自身を加振し, 穿 刺針から放射されるせん断波を C-SWE 法を用いて可視化すると, C-SWE 法で得られるせ ん断波の伝播速度図には, 穿刺針が速度の高い部分として映像化されるだけでなく, 速度 の高い部分の幅から映像面と穿刺針との距離が測定できる. (1)モデル(短軸での穿刺針の挿入) 穿刺針と超音波映像面との位置関係を Fig. 4-3-1 に示す. 2 次元(x,y)平面を考え, 穿刺針 は図に対して水平に置かれ, 穿刺針の先端（原点）から X 軸方向に距離𝑙だけ離れた位置に 超音波映像面が Y 軸に対する角度𝜃で斜めにあるとする. Fig. 4-3-1 短軸で挿入する穿刺針のモデル (2)穿刺針の 3 次元位置推定の原理 Fig. 4-3-2 に定式化のための変数の定義を示す.

29

Fig. 4-3-2 穿刺針と超音波映像面

Fig. 4-3-2 において, 超音波映像面状Δｙだけ離れた点 Q における距離差Δ𝑙は, ∆𝑙 = √∆y2_{+ (𝑑 cos 𝜃)}2_{− 𝑑 cos 𝜃}

= 𝑑 cos 𝜃 √1 + (_{𝑑 cos 𝜃}∆y )2− 𝑑 cos 𝜃 (4-3-1) 𝑑 cos 𝜃 ≫ ∆yであるから, テイラー展開すると,

≅ 𝑑 cos 𝜃 (1 +1₂(_{𝑑 cos 𝜃}∆y )2) − 𝑑 cos 𝜃 =1_{2𝑑 cos 𝜃}∆y2 (4-3-2)

ここで𝑑は穿刺針と超音波映像面までの距離である. 今, 穿刺針の先端から, せん断波が球面状に伝播しているとすると, 点 P に対する点 Q の せん断波の位相差Φは, 𝛷 = 𝑘𝑇∆𝑥 = 𝑘𝑇 2 ∆𝑦2 𝑑 cos 𝜃 (4-3-3) ここで, 𝑘𝑇はせん断波の波数であり, 𝑘𝑇= 2𝜋 𝜆𝑇= 2𝜋𝑓𝑏 𝑣𝑇 (4-3-4) また, 𝜆𝑇：せん断波の波長 𝑣𝑇：せん断波の伝播速度 𝑓𝑏 ：せん断波の周波数

30 である. 点 Q における位相差Φのために, 超音波の映像面上には, 位相の変化による偽りのせん 断波の伝播が生じているように見える. この偽りのせん断波の伝播速度𝑣𝑃は, 𝑣𝑃 = 2𝜋𝑓𝑏 𝜕𝛷 𝜕𝑦 = 2𝜋𝑓𝑏( 𝜕𝛷 𝜕𝑦) −1 _(4-3-5) （4-3-2）式より, 𝜕𝛷 𝜕𝑦

=

𝑘𝑇∆𝑦 𝒅 𝐜𝐨𝐬 𝜃であることを考慮すると, 𝑣𝑃 = 2𝜋𝑓𝑏( 𝑑 cos 𝜃 𝑘𝑇∆𝑦) (4-3-6) （4-3-3）式を代入すると, 𝑣𝑃= 𝑣_𝑇𝑑 cos 𝜃 ∆𝑦 (4-3-7) （4-3-7）式より, 穿刺針先端と映像面との距離𝑙が大きいほど超音波の映像面上に現れる 偽りの伝播速度𝑣𝑃は, 媒質の伝播速度𝑣𝑇に対して高くなる. 今, 超音波映像面上で, せん断波の伝播速度が𝑣𝑃 > 𝑣𝑇𝐻となる y 軸方向の幅 W を考え る. 𝑊 =２∆𝑦，𝑣𝑝= 𝑣𝑇𝐻とおくと, （4-3-7）式より, 𝑊 =2𝑣𝑇𝑑 cos 𝜃 𝑣𝑇𝐻 (4-3-8) 今, 穿刺針先端から超音波映像面までの距離 d は, 𝑑 =_{2 𝑐𝑜𝑠 𝜃𝑣}1 𝑇𝑣𝑇𝐻𝑊 (4-3-9) (3)精度を高める手法 Fig. 4-3-3 において, 複数の閾値を定め, 各々の𝑊𝑣_𝑇𝐻 を求めることで, 複数の𝑣𝑇𝐻∙ 𝑊から 平均値を求め, （4-3-9）式より, 𝑑𝑐𝑎𝑙𝑐を求める. Fig. 4-3-3 複数の閾値による評価

31 𝑑𝑐𝑎𝑙𝑐= ( 𝟏 𝟐 𝒄𝒐𝒔 𝜽𝒗𝑻) [𝑣̅̅̅̅̅̅̅̅̅̅̅̅̅] (4-3-9) 𝑇𝐻∙ 𝑊𝑣𝑇𝐻 Fig. 4-3-4 において, 4 方向（①X 軸方向, ②左斜方向, ③Z 軸方向, ④右斜方向）のデータ を用い, 平均値を求め, 信頼性が高まる評価をすることが可能となる. Fig. 4-2-4 多方向による評価

32 4-3 Volterra Filter を用いた画像補正 (1) 画像補正の基本原理 提案手法により得られたせん断波速度図において, 雑音と信号が混在することにより, 速度図上で穿刺針以外の領域においてアーティファクトが生じ, 穿刺針の正確な可視化が 困難となる. そこで Volterra Filter を適用させ, 速度図の速度値の大きい領域を強調しながら, 速度値 の小さい領域における強調処理を平滑化する画像補正方法を設計した[3]. 本フィルタでは 速度図の各ピクセルにおいて計算を行うが, Fig.4-4-1 に示すように, 計算するピクセルとそ の周囲を定義する. Fig. 4-4-1 画像の一部 本研究では使用している 2 次の Volterra Filter は以下の通りである. 𝑌𝑖= 𝐶1𝑋𝑖2− 𝐶2𝑋8∙ 𝑋1− 𝐶3𝑋6∙ 𝑋3− 𝐶4𝑋7∙ 𝑋2− 𝐶5𝑋5∙ 𝑋4 (4-3-1) 𝐶1, 𝐶2, 𝐶3, 𝐶4, 𝐶5は Volterra Filter のフィルタ係数である. 本研究では𝐶1= 3, 𝐶2= −1, 𝐶3= −1, 𝐶4= −0.5, 𝐶5= −0.5に定めることにする. 各ピクセルに出力する信号𝑌𝑖0は, 𝑌𝑖0= 𝑋𝑖+ 𝑌𝑖∙ ρ (4-3-2) である. (2) 画像補正の手順 Fig. 4-4-2 の流れのように, 画像補正を行う. ① 入力する速度図の信号は Volterra Filter を通じて, 新たな信号𝑌𝑖が得られる. ② 新たな信号𝑌𝑖は定数ρで乗算する. ③ ②の結果を元の入力信号に加え, 補正した出力信号𝑌𝑖0が得られる.

33

Fig. 4-4-2 Volterra Filter のブロック図

Fig.4-4-2 のような画像補正を施すことによって, 穿刺針以外で生じるアーティファクトを 抑えながら, 針位置を強調することができ, より針位置を正確に可視化することを可能に する.

34 4-4 CFI を用いた実時間映像法 前節までで述べてきたせん断波速度・方向図を用いた高精度映像法は, 速度図や方向図を 基にしているため, 高精度な 3 次元位置推定が可能であるが, リアルタイム性に乏しいとい う問題点がある. 実際に, CFI は 1[s]間に 10 フレーム取得することができるが, 速度図, 方 向図の再構成には, 最小でも 16 フレーム程度は必要であり, 2[s]程度の時間を要する. すな わち, 速度図や方向図のもととなる CFI を用いて針位置推定や d 推定, 先端推定を行う手 法を提案することで, リアルタイム性の良い測定システムとなり, 実際の臨床利用を考え たときに好ましい. そこで, 本研究ではせん断波波面の特徴に着目し, CFI を用いた針の実時間映像化を実現 した. Fig.4-6-1 に測定時のイメージを示す. Fig.4-6-1 測定イメージ ここで, 超音波映像面と針との距離 d が大きくなるとせん断波の特徴として, ・せん断波振幅が小さくなる ・せん断波位相の空間的な勾配が小さくなる 従って, それぞれの特徴から CFI では, ・CFI 上での測定推定領域が狭まる ・CFI 上に現れるせん断波波面が変化する が観測される. この CFI 上に現れるせん断波波面の変化を穿刺針や鍼の実時間映像化に使 用する. (1) 長軸挿入時 長軸挿入時には, せん断波波面の伝播に着目して d 推定を行う. Fig.4-6-2 にイメージを示 す.

35 Fig.4-6-2 CFI を用いた実時間映像法（長軸）イメージ 伝播図 針を微小振動させた場合, 媒質内部では針上の各々の点を波源とみなしことができ, 全体 として線波源として考えることができる. この時, d が大きくなった時を考えると, Fig.4-6-2 よりせん断波波面の鍼面に対する角度が大きくなっていることがわかる. このように, 長 軸挿入時には, せん断の伝播に着目して d 推定を行う. (2) 短軸挿入時 短軸挿入時には, せん断波の振動振幅に着目して d 推定を行う. Fig.4-6-3 にイメージを示 す. Fig.4-6-3 CFI を用いた実時間映像法（短軸）イメージ 振幅図 d が大きくなるに伴い, せん断波は 1/d で振動振幅が減衰する. すなわち, d が大きくなる とドプラ信号振幅が減少し, CFI 画像が現れにくくなる. このように, 短軸挿入時には, せ ん断波の振幅に着目して d 推定を行う.

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

37 Fig.5-1-1 速度図シミュレーション結果(長軸) Fig5-1-2 方向図シミュレーション結果(長軸) 以上の結果より, Fig.5-1-1 の速度図の結果より, どの d であっても針の位置を可視化できて おり, また, d が大きくなるほど速度の高い部分の幅が大きくなっていることがわかる. こ れを数値化することで詳細な d 推定が可能であることが示唆される. これについては, 次節 でシミュレーションする. 次に, 短軸方向に針を挿入した時の速度図と方向図のシミュレーション結果をそれぞれ Fig.5-1-3, Fig.5-1-4 に示す.

38 Fig.5-1-3 速度図シミュレーション結果(短軸) Fig5-1-4 方向図シミュレーション結果(短軸) 以上の結果より, 長軸と同様に Fig5-1-3 の速度図の結果より, どの d であっても針の位置 を可視化できている. また, d が大きくなるほど速度の高い円状の領域が大きくなっている ことも確認できる.

39 5-2 d 推定シミュレーション (1) シミュレーション概要 本シミュレーションでは, 取得した速度図や方向図より d 推定や先端推定を行う手法に ついて第 4 章でその理論について述べたが, それについて数値シミュレーションを行い, 正 当性を確認する. (2) シミュレーション方法 本シミュレーションのモデル条件については 5-1 節のものと同等のものを設定する. 本節 で新たに設定する条件を以下に述べる. ・伝播速度の閾値𝑣𝑇𝐻= 5.0[𝑚 𝑠⁄ ] (3)シミュレーション結果 まず, 長軸において d 推定を行った結果を Fig.5-2-1 に示す. なお, 取得した速度図は, Fig.5-1-1 であり, その速度図をもとに d 推定を行っている. ノイズ耐性や信頼性を高める ため, 複数点で d 推定(5 点)を行った平均を示す. Fig.5-2-1 d 推定シミュレーション結果(長軸) Fig.5-2-1 の 結果 よ り , 設 定 し た 距離 d と 推 定値 と の 最 大誤 差は d=2.5[mm] の時 の 0.114[mm]であり, 平均して 0.075[mm]の誤差であった. 次に, 短軸において d 推定を行った結果を Fig.5-2-2 に示す. なお, 取得した速度図は, Fig.5-1-3 であり, その速度図をもとに d 推定を行っている. ノイズ耐性や信頼性を高める ため, 複数方向で d 推定(4 方向)を行った平均を示す.

40 Fig.5-2-2 d 推定シミュレーション結果(短軸) Fig.5-2-2 の 結果 よ り , 設 定 し た 距離 d と 推 定値 と の 最 大誤 差は d=1.0[mm] の時 の 0.146[mm]であり, 平均して 0.081[mm]の誤差であった. よって, 本理論を適用させること によって, 長軸, 短軸ともに高精度に d 推定が可能であることが, シミュレーションにより 明らかとなった.

41 5-3 CFI を用いた d 推定シミュレーション (1) シミュレーション概要 本節では, 4-4 節において CFI の単一フレームから d 推定を行うための理論について述べ たがここでは, シミュレーションを通して検証を行う. (2) シミュレーション方法 本シミュレーションのモデル条件については 5-1 節のものと同等のものを設定する. な お, 本節では, CFI の波面の様子を観測するために CFI から再構成した後の伝播図をシミ ュレーションし, 伝播する波面の様子を観測する. 長軸は, 取得した伝播図を基に針から放射されたせん断波の波面の角度を計測する. 短軸は, 各 d における振動振幅をシミュレーションし, CFI 上で振動振幅が低い領域は波 面が現れない性質を再現するため, 一定の閾値以下の振幅を表示しない伝播図を出力し, 比較を行う. (3) シミュレーション結果 まず, 長軸時の伝播図のシミュレーション結果を Fig.5-3-1 に示す. d の範囲は 0.0～ 3.0[mm]の範囲で 1.0[mm]間隔で行う. Fig.5-3-1 伝播図シミュレーション結果(長軸) ここで, Fig.5-3-1 中に示す角度θをそれぞれ求めた結果を Fig.5-3-2 に示す. Fig.5-3-2 長軸角度θ測定結果（シミュレーション）

42 Fig.5-3-2 の結果より, 角度θは d が大きくなるほど直線的に大きくなっていることがわか る. よって, 本手法を d 推定に用いることが可能であることが明らかとなった. 次に, 短軸時の振動振幅のシミュレーション結果を Fig.5-3-3 に示す. d の範囲は 0.0～ 6.0[mm]の範囲で 2.0[mm]間隔で行う. Fig.5-3-3 振動振幅図シミュレーション結果（短軸） Fig.5-3-3 の結果より, d が大きくなるほど針から放射されるせん断波の振動振幅は小さくな っていることが確認できる. この時, 一定の閾値を設定し, それ以下の振動振幅の時には出 力しない伝播図をシミュレーションした. この結果を Fig.5-3-4 に示す. Fig.5-3-4 伝播図シミュレーション結果（短軸）振動振幅考慮 Fig.5-3-4 の結果より, 振動振幅を考慮した伝播図は d が大きくなるほど表示される領域が 小さくなっていることがわかる. この傾向より, 本手法を d 推定に用いることが可能である と明らかとなった.

43 第 6 章 ファントムを用いた高精度針位置推定の評価 6-1 実験系の構成と実験方法 (1)実験概要 前章までにおいて, 理論及びシミュレーションによって, 取得した CFI よりせん断波の 速度図や方向図などの各種マップによって高精度に針位置推定ができることが明らかとな った. これが実際に成り立っているのかを実験を行い検証する. 本章での提案手法として は, 手法Ⅰの「せん断波速度・方向図を用いた高精度映像法」にあたる. 6-2 節では CFI よ り再構成した各種マップの出力結果と速度図を基にした d 推定の結果および評価を行う. また, 6-3 節では実際に血管模擬したファントムによって実際に穿刺することを再現した実 験結果例を示す. (2)実験系の構成 Fig.6-1-1 に実験系の構成図を示す. Fig.6-1-1 実験系の構成 また, 実験では加振器は Fig.6-1-2 に示すリニア駆動アクチュエータを用いた. 各アクチュ エータの特性は Tab.6-1-1 のようになる. リニア駆動アクチュエータは共振系であるため, 振動周波数は共振周波数（約 300[Hz]）付近に限定される. Tab.6-1-1 アクチュエータの特性 アクチュエータ 動作電圧 [V] 共振周波数 [Hz] 最大振幅 [μm] リニア振動アクチュエータ 5 314 1000

44 Fig. 6-1-2 リニア振動アクチュエータ Panasonic の電動歯ブラシより摘出 加振器の外形を Fig.6-1-3 に示す. Fig. 6-1-3 使用加振器の外形 以下, 実験条件を述べる. ・超音波映像装置 EUB8500(Hitachi 製) ・加振周波数 276.5[Hz] ・穿刺針の振動振幅 11[μm] ・穿刺針の挿入角度 35[deg] ・せん断波の伝播速度 5.1[m/s] ・加振対象 アルギン酸ファントム (3)実験方法 1.加振器先端に取り付けた 22[G]の穿刺針を穿刺し, 穿刺針を振動させ, ファントム内部に せん断波を励起させる. 2.長軸実験時には, 超音波映像面につながれた超音波プローブをせん断波伝播方向と平行 にあてる. 短軸実験時には, 超音波映像面につながれた超音波プローブをせん断波伝播方向と垂直 にあてる. 3.Fig.6-1-4 に示すように, 図中矢印の方向に 0.5[mm]ずつ動かしていく. 4.各 d における CFI を取得し, 画像の再構成を行い, 伝播図, 速度図, 方向図など各種せん 断波伝播に関する図を取得する.

45 (a) 長軸穿刺時 (b) 短軸穿刺時 Fig. 6-1-4 実験方法 Fig.6-1-5 に実際に実験を行う際の写真（長軸実験時）を示す. なお, 設定距離 d を正確に設 定するため, 線形ステージを使用した. Fig.6-1-5 実験写真（穿刺針・長軸実験時）

46 6-2 ファントムを用いた実験結果 (1)長軸実験結果 ファントムを用いて長軸穿刺にて測定を行った後, 得られた測定結果を示す. Fig.6-2-1 に穿 刺針に振動を励起させていないときの B モード画像, 励起したときの CFI, CFI をもとに再 構成した伝播図, 速度図, 方向図の結果を d=0.0 と 2.0[mm]の時のそれぞれのものを示す. (a) d=0.0[mm] (b)d=2.0[mm] Fig.6-2-1 測定結果（長軸実験時） Fig.6-2-1 より, B モード画像では針が超音波理像面上に水平に存在しない場合 (d=2.0[mm])では, 針位置は観測することは困難であるが, 速度図では観測することが可 能である. また d=0.0[mm]の時の B モード像と速度図を比較しても明らかなように速度図 の方が針位置と他領域のコントラストの差が大きい画像であり, より明確に針位置の推定 が可能であることがわかる. 次に各 d における(d=0.0,1.0,2.0,3.0[mm])速度図及び方向図の結果を Fig.6-2-2 に示す.

47 (a)速度図 (b)方向図 Fig.6-2-2 各 d における速度図・方向図の測定結果（長軸実験時） まず, Fig.6-2-2(a)の各 d における速度図より d が大きくなるほど速度値の高い領域の幅 が大きくなっていることがわかる. さらに, (b)の方向図よりいずれの d においても穿刺針 位置の上下では, 針上が加振点のため, せん断波の伝播方向が大きく変化しており, 針先端 位置になると前方に放射していることがわかる. これらは, 理論およびシミュレーション で明らかになっており, 実際の実験でも同様のことがいえることがわかる. (2)短軸実験結果 長軸実験と同様にして, ファントムを用いて得られた測定結果を示す. 2-3 に Fig.6-2-1 と同様に各種図を示す.

48 (a) d=0.0[mm] (b)d=2.0[mm] Fig.6-2-1 測定結果（短軸実験時） 長軸と同様に, Fig.6-2-3 より針先端が超音波映像面上になくても(d=2.0[mm])速度図で先 端位置を示せており, コントラスト高い画像で針位置を示せている. また, (a3)や(b3)のせ ん断波伝播図をみると, 針の先端位置から同心円状にせん断波が伝播していることが確認 できる. 次に各 d における(d=0.0,1.0,2.0,3.0[mm])速度図及び方向図の結果を Fig.6-2-4 に示す.

49 (a)速度図 (b)方向図 Fig.6-2-4 各 d における速度図・方向図の測定結果（短軸実験時） まず, Fig.6-2-4(a)の各 d における速度図より d が大きくなるほど値の高い円状の領域の直 径が大きくなっていることがわかる. さらに, (b)の方向図よりいずれの d においても針先 端位置では 360 度全方向に放射状に伝播していることがわかる. また, せん断位置について は, d=0.0[mm]以外の d=1.0～3.0[mm]においては, 最終的に先端が超音波映像面に達する d=0.0[mm]となる位置を示せていることが確認できた.

50 6-３ 実験結果への 3 次元位置推定の定量性評価 本節では, 超音波映像面と針との距離 d を算出する理論を基に, 長軸, 短軸ともに, d 推定 の精度の評価を行う. また, 長軸実験では, 方向図を用いた先端推定の精度評価を行う. (1) 長軸実験 d 推定結果 前節の測定結果 Fig.6-2-2(a)の d=0.0～3.0[mm]の各 d の速度図の測定結果を基に d 推定 を行った結果を Fig.6-3-1 に示す. Fig.6-3-1 d 推定結果（長軸） Fig.6-3-1 の推定結果の誤差の評価を行うと, 設定距離との最大の誤差は d=1.0[mm]で 0.18[mm], d=1.0～3.0[mm]での平均誤差は 0.08[mm]であった. この結果より高精度に d の推定が行えており, これにより 2 次元画像から 3 次元的な位置の推定が高精度に行えて いることが実際の実験より明らかとなった. なお, d=0.0～1.0[mm]においては穿刺針自体 の太さの影響を受けて誤差が生じたことが考えられる. (2) 短軸実験 d 推定結果 次に長軸実験と同様にして, 前節の測定結果 Fig.6-2-4(a)の d=0.0～3.0[mm]の各 d の速 度図の測定結果をもとに d 推定を行った結果を Fig.6-3-2 に示す.

51 Fig.6-3-2 d 推定結果（短軸） Fig.6-3-2 の推定結果の誤差の評価を行うと, d=1.0～3.0[mm]での設定距離との最大誤差は d=2.0[mm]で 0.38[mm], 平均誤差は 0.18[mm]であった.この結果より, 長軸実験と同様に 高精度に d 推定が行えていることがわかる. (3) 長軸実験での先端推定精度評価 ここでは, d=0.0[mm]での B モード画像と方向図を用いた先端推定の精度について評価 する. まず, 先端推定を行う際に, 速度図ではなく方向図を用いる理由を述べる. 針の先端推定 を行う際, 速度図を用いると針位置以外のアーティファクトの影響を強く受け, 正確さに 欠けるという問題点がある. Volterra Filter を通した場合でも完全に取り除くことができな いケースもある. そこで, 方向図を基にした新たな手法を提案する. では, 実際に推定する方法について述べる. 長軸実験時, 針位置の上下では, 針上が加振 点となっているため, せん断波の伝播方向が大きく変化している. また, 針の先端位置では, せん断波が 360 度全方向に伝播する特異性を持っている. この変化点を先端推定に用いる. Fig.6-3-3 にその流れを示す.

52 Fig.4-5-2 方向情報を用いた先端推定法(長軸) 針位置以外の箇所では Fig.4-5-2 に示すようにピクセル間の差が小さく, 針位置ではその差 が大きいという特徴がある. この差を用いて方向差図を求める. その後, 閾値𝐷𝑖𝑟𝑑𝑖𝑓_𝑡ℎを設 定し, 図中 x 軸の右から比較を行い, 最初に閾値以上となる点を先端と推定する. 次に, 実際に穿刺針をファントムに挿入した際の結果(d=0.0[mm])を基にして, 幾何学的 に求めた先端位置との誤差を測定し, 本手法の評価を行う. まず, B モード画像と先端位置 の拡大図を Fig.4-5-3 に示す. Fig.4-5-3 先端推定（B モード画像） Fig.4-5-3 より, 先端位置と他の位置とでのコントラストの差が小さく, また先端位置で白 く表示されている領域があるが正確に推定することはできない. よって, B モード画像での 高精度な先端位置推定は困難であることがわかる. 次に, 速度図と先端位置の拡大図を Fig.4-5-4 に示す.

53 Fig.4-5-4 先端推定（速度図） Fig.4-5-4 より, 速度図での先端推定は, 針位置以外でも速度の高い領域があるため, 先端 位置を 1 点に決めることができない. よって, 速度図より高精度な先端位置推定は困難であ ることがわかる. さらに, 方向図を用いた先端位置の結果を Fig.4-5-4 に示す. Fig.4-5-5 先端推定（方向図） 今回, 方向図より推定した先端位置は, 方向図上でせん断はが 360 度全方向に伝播する点 を先端と推定しており, Fig.4-5-4 では黄色十字がぶつかる点(A 点)である. 幾何学的に求め た先端位置と A 点の誤差はおよそ 0.49[mm]に抑えることができている. したがって, せん 断波の方向図を用いることで先端位置を 1 点に定めることができ, また, 高精度に推定す ることが可能である.

54 6-4 血管模擬したファントムの実験結果 穿刺針で穿刺する際には, 血管などといった対象物に向けて針を挿入していく. 本節で は, その状況を模擬し, 正確に位置ナビゲーションを行うことができるか評価を行う. (1)実験方法 今回実験を行う際に作成したファントムを Fig.6-4-1 に示す.血管を模擬するために直径 8[mm]の穴をあけたファントムを使用する. Fig.6-4-1 血管模擬ファントム 実際に穿刺している様子を Fig.6-4-2 に示す. 今回実験を行う際の実験系は 6-1 節に従う. Fig.6-4-2 血管模擬実験穿刺の様子 (2)実験結果 実験を行って得られた結果を Fig.6-4-3 に示す. なお,(a)が穿刺針の先端が血管に到達する 前, (b)が到達した後を模擬している. それぞれ B モード像, CFI, 伝播図, 速度図, 方向図を 示す.

55 (a) 穿刺針先端位置：血管外 (b) 穿刺針血管位置：血管内 Fig.6-4-3 血管模擬実験測定結果(d=0.0[mm]) Fig.6-4-3 より B モード画像に比べ速度図や方向図を用いることでコントラスト高い画 像で針位置を示すことができている. また, 針先端位置も正確に追従可能であり, 針が 血管到着前後での違いは明確であり, 対象物に向けて穿刺する際の正確な位置ナビゲー ションとして, 適していることがわかる.

56 第７章 CFI を用いた実時間映像法の評価 7-1 実験系の構成と実験方法 (1)実験概要 第 6 章では、手法Ⅰのせん断波の速度・方向図を用いた高精度映像法についての実験結 果及び定量評価について述べたが, 本章では, 手法Ⅱの CFI を用いた実時間映像法につい ての実験結果について述べる. (2)実験系の構成 Fig.7-1-1 に実験系の構成図を示す. Fig.7-1-1 実験系の構成 また, 実験では Fig.7-1-2 に示す超小型リニア駆動アクチュエータを用いた. Fig.7-1-2 使用加振器（リニア振動アクチュエータ）

57 本加振器の性能について, 以下に示す. ・共振周波数：170[Hz] ・駆動電圧：0.1～2.5[V] ・駆動電流：170[mA](定格) ・外形：円型（直径 10[cm]） 今回の実験は鍼を挿入することとし, 使用する鍼を Fig.7-1-3 に示す. なお, 鍼の直径は 0.08[mm]で, 鍼長は 40[mm]であるものを用いた. Fig.7-1-3 使用鍼 以下, 実験条件を述べる. ・超音波映像装置 EUB8500(Hitachi 製) ・加振周波数 276.5[Hz] ・加振器の振動振幅 60[μm] ・せん断波の伝播速度 2.6[m/s] ・加振対象 アルギン酸ファントム (3)実験方法 1.超小型加振器をゴム手袋内に装着させる.(Fig.7-1-3(a)) 2.超小型加振器が装着されているゴム手袋と鍼の鍼柄を Fig.7-1-3(b)のように接触させる. 3.長軸実験時には, 超音波映像面につながれた超音波プローブをせん断波伝播方向と平行 にあてる. 短軸実験時には, 超音波映像面につながれた超音波プローブをせん断波伝播方向と垂直 にあてる. 4.超音波映像面と鍼との距離 d を 1.0[mm]ずつ変化させながら, CFI を取得する.

58 (a) 加振器装着の様子 (b) 鍼所持の様子 Fig.7-1-3 超小型加振器装着方法 Fig.7-1-4 に実際に実験を行う際の写真（長軸実験時）を示す. Fig.7-1-4 実験写真（鍼・長軸実験時）

59 7-2 CFI を用いたリアルタイム実験結果 (1)長軸実験結果 ファントムを用いて,長軸刺鍼にて測定を行った後, 得られた測定結果を示す. なお, d の値 を 0.0～3.0[mm]の範囲で 1.0[mm]ずつ変化させた際の結果を Fig.7-2-1 に示す. Fig.7-2-1 鍼測定結果（長軸実験時） Fig.7-2-1 中に示す先端位置付近に着目すると, せん断波の波面の伝播方向が鍼上と比べて 広がっていることが確認できる. すなわち, せん断波の伝播の様子を観測することで先端 位置を推定することができる. 本実験で用いた鍼(直径 0.08[mm])は B モードでは観測する ことが困難であるが, 鍼を微小振動させ, CFI を取得することで鍼位置を可視化することを 可能にした. また, CFI を用いることで CFI のフレームレートは 10[fps]であり, リアルタイ ム性良く測定することができる. (2)短軸実験結果 ファントムを用いて, 短軸刺鍼にて測定を行った後, 得られた測定結果を示す. なお, d の値 を 0.0～3.0[mm]の範囲で 1.0[mm]ずつ変化させた際の結果を Fig.7-2-2 に示す. Fig.7-2-2 鍼測定結果（短軸実験時） Fig.7-2-2 中に示す先端位置付近に着目すると, その点からせん断波が放射的に伝播してい ることがわかる. すなわち, その点を先端と推定することができる. (1)(2)の結果より, CFI を用いることでリアルタイム性良く鍼位置を可視化することができ ることが, 実際の実験においても明らかとなった.

60 7-3 CFI を用いたリアルタイム 3 次元位置推定の定量性評価 本節では, 7-2 節で示した CFI での測定結果より, CFI から d 推定を行うための理論を基 に実際の実験結果で評価を行ったので, その結果を示す. (1)長軸実験時 長軸実験時には, せん断波の伝播に着目して, CFI での針から放射される伝播する波面の 角度を計測することで, d 推定が可能である. よって, 実際に波面の角度の計測を行う例を Fig.7-1-3 に示す. Fig7-3-1. 波面の角度計測例 この時, CFI はせん断波の波面は位相が変化するのに伴い, 常に Fig.7-3-1 に示す波面とな らない. そのため, 各 d での 10 フレーム分の角度θを計測し, その平均をグラフ化したも のを Fig.7-3-2 に示す. Fig.7-3-2 波面角度計測結果（長軸実験時） Fig.7-3-2 より, d の値が大きくなるほど波面の角度𝜃𝑎𝑣𝑒が直線的に大きくなっていること が確認できる. Fig.7-3-2 中に示す点線は近似直線であり, その直線との d 推定の最大誤差 は d=1.0[mm]の時 0.08[mm]で, 0.0～3.0[mm]の範囲の平均誤差は 0.05[mm]であった. よ って, せん断波波面の角度を計測することで, CFI の単一フレームでリアルタイムに d 推定 を行えることが明らかとなった.

61 (2)短軸実験時 短軸実験時には, せん断波の振動振幅に着目して, CFI の ROI 内に現れるせん断波波面 の面積を求めることで, d 推定が可能である. まず, 今回測定した CFI を Fig7-3-3 に示す. Fig.7-3-3 CFI 結果例(短軸実験時) この得られた CFI に対し, せん断波の波面が現れる領域に色付けを行う二値化処理を行っ た結果を Fig.7-3-4 に示す. Fig.7-3-4 波面表示領域例(Quality 図) ここで, 出力された図の中で色付けされたピクセルの数を𝑄𝑢𝑎𝑙𝑖𝑡𝑦と定義する. 各 d での 10 フレームの𝑄𝑢𝑎𝑙𝑖𝑡𝑦の平均𝑄𝑢𝑎𝑙𝑖𝑡𝑦𝑎𝑣𝑒をグラフ化したものを Fig.7-3-5 に示す. なお, CFI 内 の最大のピクセル数は 31,450 である.

62 Fig.7-3-5 Quality 結果（短軸実験時） Fig.7-3-5 より, d が大きくなると𝑄𝑢𝑎𝑙𝑖𝑡𝑦 値が直線的に小さくなっており, ROI 内に現れる せん断波の領域は小さくなっていることが確認できる. これは, Fig.7-3-4 の図を見ても ROI の下部領域が顕著にその傾向が表れていることがわかる. また, Fig.7-3-5 中の点線は近似直 線であり, その直線との d 推定の最大誤差は, d=4.0[mm]の時 0.29[mm]であった. さらに, d=0.0～6.0[mm]の範囲の平均誤差は 0.19[mm]であった. よって, せん断波の振幅をせん 断波の波面の領域を求めることで, CFI の単一フレームでリアルタイムに d 推定を行えるこ とが明らかとなった.

63 7-4 鶏胸肉を用いた実験結果 生体内部への刺鍼は, 内部の繊維などの影響を受け, CFI にノイズを含み明瞭な二値化画 像とならないこともある. 本節では, 鶏胸肉を対象に刺鍼し, 測定を行い, CFI の単一フレ ームから画像処理を行い, ノイズの除去された二値化画像を生成した結果を示す. (1)実験方法 実験条件や実験系は 7-1 節に従うものとする. 実際に刺鍼している様子を Fig.7-4-1 に示す. Fig.7-4-1 鶏胸肉実験刺鍼の様子 なお, 今回二値化画像を生成する際の画像処理の流れを Fig.7-4-2 に示す. 本画像処理によ って, CFI の流速推定値が最大となる領域に色付けがなされる二値化画像が生成される. Fig.7-4-2 二値化画像処理の流れ