領域別フィルタリングによるC-SWE法の高分解能化

令和元年度

修 士 論 文

領域別フィルタリングによる

C-SWE 法の高分解能化

指導教員 山越 芳樹 教授

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

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

T181D009 伊藤 拓海

1

領域別フィルタリングによる C-SWE 法の高分解能化

目次 第１章 序論 第２章 せん断波による硬さ計測 2-1 せん断波とは 2-2 生体内組織におけるせん断波の伝播 2-3 せん断波計測で期待されるパラメータと有用性 2-4 高分解能化への要求 第３章 連続せん断波映像法 [Continuous SWE] 3-1 超音波パルスドプラ法による組織内振動伝搬計測 3-2 カラーフロー映像系（CFI）の流速推定アルゴリズム 3-3 CFI の流速推定アルゴリズムによるせん断波の波面検出 3-4 定量的なせん断波画像の構成法 3-5 波数ベクトルフィルタリング 第４章 連続せん断波映像法の実験系と特徴 4-1 提案手法の実験系と特徴 4-2 測定の課題と対策 4-2-1 貼り付け型加振器の導入

4-2-2 せん断波伝播方向指標 Shear Wave Propagation Direction Index の導入 4-3 提案手法の信頼性 4-3-1 ARFI 法との比較 4-3-2 測定における信頼性 第５章 領域別フィルタの導入 5-1 要求分解能と高分解能化の手法 5-2 領域別フィルタのアルゴリズム 5-2-1 全体フロー 5-2-2 各フィルタの概要 5-3 先験的情報別の分割結果 5-4 分割結果まとめ

2 第６章 領域別フィルタによる分解能評価 6-1 シミュレーションによる分解能評価 6-2 二層ファントム測定による分解能評価 6-3 分解能評価まとめ 第７章 運動負荷を加えた際の僧帽筋硬度の経時的な変化計測 7-1 実験系・実験プロトコル 7-2 高分解能化によるせん断波伝播 7-3 実験結果 7-3-1 運動負荷による伝播速度変化量 7-3-2 Control 伝播速度への依存性 第8 章 結論 謝辞・参考文献

3

第

1 章 序論

スポーツ庁の「スポーツ実地状況等に関する世論調査」によると，日本のスポーツ人口は 年々増えてきており，スポーツ医学への関心も高まってきている．近年，整形外科等のスポ ーツ医学の分野では，筋肉の硬さ「筋硬度」を数値で評価する手法が求められているが，筋 肉の硬さや損傷の評価基準は医師や患者の主観によるものが多く，定量的に筋肉の緊張や 損傷の度合いを表すパラメータが必要である．筋硬度を定量的に測定することで，リハビリ テーションや筋肉トレーニング，高齢者医療を客観的な指標で行うことが可能となる．近年， 筋硬度を定量的に測定する様々な方法が研究されているが，新たな方法として超音波エラ ストグラフィによる組織弾性計測がある． 周波数1 [kHz]程度までの振動を，生体表演から生体組織などの比較的柔らかい物体に加 えると，その放射エネルギーの大部分は組織中を横波として伝播する．そして，その伝播速 度や減衰係数は，せん断粘弾性パラメータと関連していることが知られている．また，生体 組織のせん断粘弾性特性は，生体組織を触った時の硬さや感触と密接に関係しているため， 生体組織に低周波振動を加えた時の伝搬速度や減衰などを測定することができれば，疾病 の進行度の定量的な評価や癌等の病気の早期発見が期待できる．しかし，生体組織内ではせ ん断波が組織境界で屈折，反射するため伝播が非常に複雑になるため，測定精度に影響を及 ぼす．そのため，伝搬速度を精度良く測定できるのは肝臓などの比較的一様な組織に限定さ れてしまう．しかし，複雑な境界面を持つ非一様な組織においても，精度よく組織内部の粘 弾性を測定できるシステムが必要である．さらに，臨床においては信頼でき，分かりやすい 測定結果の画像化等も求められている． 現在，超音波エラストグラフィとして広く用いられているものは，ARFI 法，SWE であ る．それらの方法は，機械的振動波であるパルスせん断波を組織に伝播させ，その伝播速度 を超音波映像装置で測定し組織の弾性を評価している．しかし，この手法が搭載された超音 波映像装置は非常に高価であり使用できる医療機関が限られてしまうこと，また，骨付近で 測定を行う際に生体組織の温度上昇が懸念されること，伝播速度にバラつきがあるという 問題点がある．そのため，より安価で安全，定量的な測定方法を開発する必要がある． そこで本研究では，測定対象を体表近くの組織をターゲットとした新規のせん断波エラ ストグラフィとして，連続せん断波映像法（Continuous Shear Wave Elastography :C-SWE 法）を提案する．本稿では，連続せん断波を組織に伝播させ，その伝播速度を計測する連続 せん断波エラストグラフィを用いて，生体模擬ファントムや骨格筋を測定し評価を行った． 特に僧帽筋をターゲットとして，運動負荷を加えた時の筋弾性の経時的な変化を測定した． また，骨格筋を測定していく中で，腱や筋膜等の微細な組織の測定では要求分解能が高く， 分解能の向上が必要であることがわかった．そこで、先見的情報から弾性が異なると思われ る複数領域に分割し，分割された領域に対して個別にフィルタをかけることで，高分解能化 を図ったことについて報告する．

4

第

2 章 せん断波による硬さ計測

本章ではせん断波の特徴と工学的な研究課題，生体軟組織内部における低周波振動の伝播 について示す．さらに，せん断波計測により期待される臨床意義や目的について示す．

2-1 せん断波とは

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

5

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-3) 𝜇 = 𝜇1+ jω𝑣𝜇2𝜆 = 𝜆1+ 𝑗𝜔𝑣𝜆2 𝜇1 : せん断弾性係数𝜆1∶ 体積弾性係数 𝜇2 : せん断弾性係数𝜆2∶ 体積弾性係数 ρ ∶ 密度 ω𝑣: 振動周波数 𝑅𝑒[ ]，𝐼𝑚[ ]：[]内の複素数の実数部，虚数部 また，これら縦波や横波の他に生体の表面付近では表面波が存在するが，この伝播速度は ほぼ横波の伝播速度に等しいことが知られている．上記の波動の中で，縦波は圧縮性の波 であり，媒質を圧縮することにより伝播する．一方，横波は非圧縮性の波であり，媒質を 等体積のまま，横方向に挟み切るように変形させながら伝播していくため，せん断波とも 呼ばれている．ここで，周波数が1[kHz]程度以下の低周波振動であると，外部から与えら れた振動のエネルギーはそのほとんどが横波に変換されると考えられている[2]．

6 ここで，(2-2-4)式，(2-2-5)式で与えられる横波の伝播速度と減衰係数を，粘弾性パラメー タを用いて書くと， 𝑣𝑡= √ 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）

7

2-3 せん断波計測で期待されるパラメータと有用性

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

8

2-4 高分解能化への要求

SWE において，適用範囲における要求分解能を Fig.2-4-1 に示す． Fig.2-4-1 高分解能化への要求 現在のC-SWE 法の分解能は 5[mm]程度である． Fig.2-4-1 より，骨格筋の硬さ測定においては現在の C-SWE 法で十分に測定可能である といえる．しかし，乳がんの診断や腱の測定，そしてファッシアや穿刺針の測定においては 現在のC-SWE では測定が困難であるといえる． それらを測定するためにC-SWE 法の高分解能化が必要である．

9

第

3 章 連続せん断波映像法[Continuous SWE]

3-1 超音波パルスドプラ法による組織内振動伝播計測

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

ξ(𝑡) = 𝜉

₀

𝑠𝑖𝑛(2π𝑓

_𝑏

𝑡 + 𝜙

_𝑏

)

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

𝜙

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

𝑓 =

𝑐+𝑣(𝑡) 𝑐

𝑓

0 (3-1-2) 𝑐：音速

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-11) となり，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𝜋𝑓₀ 𝑐 2𝑣 𝑖 Δ𝑡) (3-2-3) 𝑄𝑖= 𝑎 sin (𝜙0+ 2𝜋𝑓0 𝑐 2𝑣 𝑖 Δ𝑡) (3-2-3)式は，(3-2-4)式のように書くこともできる． Q⃗⃗ 𝑖= 𝑎 𝑒𝑥𝑝( 𝑗(𝜙0+ 2𝜋𝑓₀ 𝑐 2𝑣 𝑖 Δ𝑡)) (3-2-4) ここで，第i 番目の超音波パルスの位相と，第 i+1 番目の超音波パルスの位相の差Δ𝜙𝑖を考 える．これは， Δ𝜙𝑖= 𝑎𝑟𝑔 (Q⃗⃗ 𝑖+1Q⃗⃗ 𝑖 ∗ ) (3-2-5)

13 と推定できるので，(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 番目の受信超音波パルスの位相𝜙𝑖は，

𝜙

_𝑖

= 𝜙

₀

+

2𝜋𝑓0 𝑐

2𝜉

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

𝜔

_𝑏

=

2𝜋 4Δ𝑡 (3-3-4) つまり，せん断波の周波数であらわすと， 𝑓𝑏= 1 4Δ𝑡 (3-3-5) さらに，振動の初期位相として 𝜙0= 0 (3-3-6) が満たされるとする．

15 上記条件（(3-3-5)式および(3-3-6)式）は，せん断波の伝播による組織の変位振動の周期が 超音波の４パルスに等しく，かつ初期位相が 0 の条件であり，これを変位振幅として図に 表すと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) となる．

𝛥𝑡

2𝛥𝑡

𝜉

3𝛥𝑡

4𝛥𝑡

0

𝑡

16 ここで，i=0,1,2,3 について，直交検波器の出力を求めてみると， 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 となる．

17 (3-3-11)-(3-3-15)式の関係をベクトル図であらわすと ① 0 ≤ 𝜉0≤ 𝜆 8 の場合 Fig.3-3-2 に示すように，すべてのベクトルは第一象限と第四象限にある． 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 のようになる． Fig.3-3-7 CFI における流速導出アルゴリズムを使ったせん断波の波面再生

20 ここで， { 0 ≤ 𝜉0≤ 𝜆 8 のとき 𝐸𝐿≥ 0 𝜆 8 ≤ 𝜉0≤ 3𝜆 8 のとき 𝐸𝐿≤ 0 となるが，ともにEU=0 であるので，実軸を EL，虚軸をEUとするベクトルは，ELが正の 場合は実軸上の正の方向を向くベクトルとなり，流速推定値は0 になる．一方，ELが負の 場合は実軸上の負の方向を向くベクトルとなり，流速推定値は正の最大値，または負の最大 値（ナイキスト周波数で決まる最大の流速値）になる． つまり， ① ELが正になる条件（せん断波による振動振幅が0 ≤ 𝜉₀≤𝜆 8 の場合） 流速0 になる． ② ELが負になる条件（せん断波による振動振幅が𝜆 8 ≤ 𝜉0≤ 3𝜆 8の場合） 振動振幅の位相が0 度，および 180 度になる位置で CFI 画像には流速最大の部分が 現れる． この条件は，せん断波の振幅により，せん断波による振動位相が0 または 180 度の時に， 特異なパターンがCFI 画像に現れることを示しており，これを振幅条件と呼ぶ． せん断波が組織中を伝播しているとき，CFI 画像の中から上記に示したような特徴ある 部分を抽出することにより，せん断波の位相（0 度または 180 度）が推定できることにな る．せん断波が等位相になる部分はせん断波の波面を再現することに相当するので，この方 法により，CFI 画像からせん断波の波面を再現できることになる． この方法は，周波数条件（3-3-7 式）が成り立つときに，CFI の推定アルゴリズムが，せ ん断波の0 度と 180 度の位相を検出するディジタルフィルターになっていることに着目し た，せん断波の映像化法である．横軸を初期位相𝜙𝑏，縦軸を振動振幅𝜉0として，以下の条件 で，流速推定の数値シミュレーションをおこなった結果をFig.3-3-7 に示す．

21 [シミュレーション条件] 超音波中心周波数 𝑓0 6.5𝑀𝐻𝑧 超音波伝播速度 𝑐 _{1500 𝑚 𝑠}⁄ パルス繰り返し周波数 1 𝑑𝑡⁄ 365𝐻𝑧 パルス本数 𝑁 11 加振周波数 𝑓𝑏 91.25𝐻𝑧 Fig.3-3-7 数値シミュレーション結果 周波数条件は，理論的には（3-3-7）式で表されるが，実際には，せん断波の周波数がこ の条件に近いときでも，流速の最大値または流速0 の部分が 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 波数ベクトルフィルタリング

反射波により定在波が存在すると，せん断波の位相は空間的に変調されてしまう．せん断 波の複素振幅マップを二次元フーリエ変換して空間周波数上のスペクトルを計算し，波数 ベクトルフィルタを適用してフーリエ逆変換することで任意の成分を抽出することができ る．この方法をCD SWI 法で取得したせん断波位相マップに適用することを考える． 前方へ伝播する入射波と，後方へ伝播する反射波が観測された一次元について考えると， せん断波の複素振幅𝑆(𝑥)は 𝑆(𝑥) = 𝐴𝐹𝑒𝑥𝑝[𝑗(𝑘𝑝𝑥 + 𝜑𝐹)] + 𝐴𝐵 𝑒𝑥𝑝[𝑗(−𝑘𝑝𝑥 + 𝜑𝐵)] (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 提案手法の実験系と特徴

せん断波励起部では，発振器からの低周波信号を加振器から印加し，第三章で示した周波 数条件を満たした１ｋHz 程度以下の連続的な振動を骨格筋端部に加え，せん断波を筋繊維 方向に伝搬させる．超音波プローブは骨格筋に並行になるように当てて筋組織を描画する． このとき血流の映像化に用いられる超音波カラーフロー画像を採取するが，この画像上に は第三章で示した特定条件下で，せん断波伝搬に起因した波状パターンが現れるので，これ をPC 内に画像インターフェースを介して実時間で取り込み，せん断波波面を再現する． 本研究では，超音波映像装置にEUB-8500(日立メディコ)，EUB7500(日立メディコ)を用 いた．提案手法の測定系をFig.4-1-1 に示す． Fig.4-1-1 測定系 超音波装置の画像を画像処理用PC へ取り込むために使用した，コンバーター・キャプ

25 Fig.4-1-2 超音波装置画像取得装置 提案手法の特徴として，以下のものが挙げられる。 ① 汎用の超音波映像装置のカラードプラ機能である流速検出アルゴリズムをせ ん断波の波面検出に使うため，超音波装置本体の改造を一切必要としない．ま た，超音波映像装置のビデオ出力を画像処理することで，連続的なせん断波を 使う組織弾性の映像系が構成できる． 汎用の超音波装置に，加振源と画像処理用のPC，専用の画像処理ソフトを付 けることで，簡単に定量性の高い組織弾性映像系が構成できる．（組織弾性の映 像系が簡便な方法で得られる．） ② 超音波映像装置の流速検出アルゴリズムをせん断波の波面検出に活用してい るので，せん断波映像を得るために必要な信号処理能力は一般のPC で十分で あり，このため実時間でせん断波の波面が組織中を伝播していく様子が動画像 で確認できる．（実時間で波面の動きが再生されるため，波面の動きや伝播方向 の乱れから組織の機械的なマクロ構造が観察できる．） ③ 連続的なせん断波（周波数 1[kHz]程度以下の生体表面からの振動の印加で生 体組織中に励起される）を使っているので，超音波の放射圧を使う方法（シーメ ンス社Virtual Touch 等）に比べて生体への高い安全性を有する．また静圧を生 体表面から印加しそのときの組織ひずみを映像化する方法（日立，エラストグラ ィ等）に比べて，せん断波を用いた計測により定量性が高い． ④ 連続的なせん断波を用いるため，反射波等により定在波が発生してしまうと速 度推定に誤差を生じてしまう．しかし，波数ベクトルフィルタによる反射波の除 去や，実時間での波面観測による加振点の選定により，改善することが可能．

26

4-2 測定の課題と対策

骨格筋に対して SWE を適用した例は多いが，その中で代表的な文献からせん断波伝播速 度の計測値をまとめ，せん断波伝播速度計測における必要精度を検討する． 疾患の有無や，ある動作をした際の骨格筋のせん断波伝播速度をまとめたものを Fig.4-2-1 に示す． Fig.4-2-1 必要測定精度

Fig.4-2-1 より Neck Pain の有無による僧帽筋の伝播速度の差は 12％であり，変化が小

さいことがわかる．これより，骨格筋の評価の際には誤差5％程度の測定精度が必要である と考えられる． また，測定における課題とその対策をFig.4-2-2 に示す． Fig.4-2-2 本手法の課題と対策

4-2-1 貼り付け型加振器の導入

従来，測定者が手に持ち加振する加振器を使用していたが，この方法では測定ごとに加 振位置が異なってしまう可能性から位置の再現性が低い．また，長時間の測定では測定者 の持ち手が疲れてしまう．などといった課題があった．そこで本手法では生体に貼り付け ることが可能である小型加振器を導入した．

27 今回使用する加振器はオーディオエキサイタ(Tectonic Elements 製)を使用した．オーデ ィオエキサイタの詳細を以下に示す． 定格出力 0.8W 直流抵抗 7.8Ω インダクタンス 0.13mH 重量 18.6g 以上より、定格出力で2.5V の出力電圧を得ることができる． 従来の加振器，新加振器をFig.4-2-1-1 に示す． Fig.4-2-1-1 加振器 また，導入した加振器を生体に貼り付け使用している例をFig.4-2-1-2 に示す． 生体に貼り付ける際は，テーピングテープを使用して張り付けた． Fig.4-2-1-2 生体貼り付け例

28 この新加振器は重さ 18.6g と計量であるため，生体に貼り付け使用しても被検者の負担 にはならない．この貼り付け型加振器により，測定ごとの位置ずれを抑え位置の再現性を高 めることが可能である． また，加振器用の制御回路(Fig.4-2-1-3)には発振器が一体となっているものを用いた． Fig.4-2-1-3 ドライバ加振器用の制御回路 これらを使用した際の振動振幅(Fig.4-2-1-4)，温度特性(Fig.4-2-1-5)を示す． Fig.4-2-1-4 加振器の出力電圧と振幅の関係 Fig.4-2-1-4 より定格出力で 210um の振幅を得られるため，振幅条件を満たすことを確 認できる．

29

Fig.4-2-1-5 定格駆動における加振器の温度変化

Fig.4-2-1-5 より 2.5V 駆動で 30 分使用時の最高温度は 35.9°C であるため，生体で使用 しても問題ないことを確認できた．

4-2-2 せん断波伝播方向指標 Shear Wave Propagation Direction Index の導

入

せん断波は生体内で 3 次元的な複雑な伝播をしているため，速度推定において推定誤差

が生じてしまう．そこで，せん断波の伝播の一様性を評価する指標として SWDI(Shear

Wave Propagation Direction Index)を導入した．SWDI は式(4-2-2-1)で表させる．

(4-2-2-1)

ただし，

M : x 方向のピクセル数 N : z 方向のピクセル数

30

Fig4-2-2-1 θ の基準 SWDI とせん断波の伝播角度との関係を Fig.4-2-2-2 に示す．

Fig.4-2-2-2 SWDI とせん断波伝播角度の関係

31 Fig.4-2-2-3 SWDI と伝播図 この時のSWDI とせん断波伝播速度の関係を Fig.4-2-2-4 に示す． Fig.4-2-2-4 SWDI と伝播速度 Fig.4-2-2-3，Fig.4-2-2-4 より，SWDI が高い値を示すほどせん断波の伝播は一様であり， せん断波伝播速度の変動が小さいことがわかる．SWDI が 95 点以上のデータを選択した場 合，変動係数を5％以内に抑えることが出来ると推察される． Fig.4-2-2-2 より，SWDI を 95 点以上にするためには，せん断波の伝播角度を約 15 度以内 に収める必要があるとわかる． これらのことより， SWDI を利用することにより，精度高く測定ができるようになると考 えられる．

32

4-3 提案手法の信頼性

本節では，提案手法で得られた伝播速度の信頼性について記述する．

4-3-1 ARFI 法との比較

まず初めに，提案手法とARFI 法の比較を行う． ARFI 法は，超音波による生体組織の硬さ評価に使われる別の手法である． ARFI 法（音響放射圧法） 機械的振動波であるせん断波を，比較的強力な超音波の収束により発生させる．組織中に 伝播するせん断波の伝播速度を超音波で測定して組織の弾性を評価する方法．高フレーム レートの高価なエコー装置が必要で，測定ごとにクーリング時間が必要である． Fig.4-3-1-1 に提案手法と ARFI 法の比較を示す． Fig.4-3-1-1 提案手法と ARFI 法の比較 ARFI 法との比較として，提案手法は連続波を使用しているため，測定可能なせん断波の 伝播速度に上限がないこと，また，伝播方向等も解析的に導出することが可能であるとこが 挙げられる．

33

次に，生体模擬ファントムを使用して提案手法とARFI 法の速度測定値の比較を行った．

[実験条件]

超音波映像装置：ARFI 法 ACUSON S3000 (Siemens) C-SWE 法 EUB7500(HITACHI) 加振周波数：76.0Hz

使用ファントム：OST 製ファントム（No.4，No.8）Fig.4-3-1-2 に示す．

Fig.4-3-1-2 OST 製ファントム（No.4，No.8） [実験方法]

1. 加振器でファントム表面を振動させ，ファントム内部にせん断波を励起させる． 2. 超音波映像装置につながれた超音波プローブをせん断波伝搬方向と平行に当て

る．

3. ファントムを変えて１，２を繰り返す．

Fig.4-3-1-3 に ARFI 法で得た速度図，C-SWE 法で得た伝播図と速度図を示す．

34

Fig.4-3-1-3 二つの方法でのファントム測定結果

Fig.4-3-1-4 得られた速度値

Fig.4-3-1-3 の速度図より，どちらの手法でも一様な伝播得られたことがわかる．また， C-SWE 法では伝播図より，せん断波の伝播の様子も確認できる．

Fig.4-3-1-4 では，C-SWE 法と ARFI 法を PhantomNo.8 でキャリブレーションした後， No.4 のファントムで速度比較をした結果を示し，二つの方法での伝播速度はほぼ同じであ ることを確認できる．

35

4-3-2 測定における信頼性

次に，提案手法での測定における信頼性について記述する． 提案手法の測定の信頼性を評価するため実験を行った． 再現性は，検者間信頼性と検者内信頼性の二つで評価を行う． まずは，検者内信頼性について記述する． 1 人の検者が 1 人の被験者に対して測定をすることで，測定の信頼性を評価した． [実験条件] 超音波映像装置：EUB7500(HITACHI) 加振周波数：76.0Hz 測定対象：僧帽筋 [実験内容] １人の検者が１人の被験者に対して二回測定を行った． 加振器は張り付けたまま固定し，プローブは生体から離す．しかし，ペンで生体にマーク を付けることで，同一部位を測定できるようにした．測定では，SWDI 値が９５点以上のも のを使用し，３データを１測定分とて2 測定間で有意差検定を行った． 測定写真をFig.4-3-2-1，結果を Fig.4-3-2-2 に示す． Fig.4-3-2-1 僧帽筋の測定

36 Fig.4-3-2-2 検者内再現性結果 Fig.4-3-2-2 より，二回の測定ともに，同じような伝播図が得られたことが確認できる． せん断波平均伝播速度値は，1 測定目：3.20[m/s]，2 測定目：3.11[m/s]となり，二測定の間 の有意確率は0.213 と有意差は認められなかった． 以上のことから，検者内信頼性において，本測定方法で信頼性ある測定ができているとい える．

37 次に，検者間信頼性について記述する． ２人の検者が1 人の被験者に対して測定を行った際の信頼性を評価した． [実験条件] 超音波映像装置：EUB7500(HITACHI) 加振周波数：76.0Hz 測定対象：僧帽筋 [実験内容] 2 人の検者が１人の被験者に対して 1 回ずつ測定を行った． 加振器は張り付けたまま固定し，プローブは生体から離す．しかし，ペンで生体にマーク を付けることで，同一部位を測定できるようにした．測定では，SWDI 値が９５点以上のも のを使用し，３データを１測定分とて2 測定間で有意差検定を行った． 結果をFig.4-3-2-3 に示す．

38 Fig.4-3-2-3 検者間再現性結果 Fig.4-3-2-3 より，2 人の検者の測定で，同じような伝播図が得られたことが確認できる． せん断波平均伝播速度値は，検者1：3.20[m/s]，検者 2：3.13[m/s]となり，二測定の間の有 意確率は0.399 と有意差は認められなかった． 以上のことから，検者間信頼性において，本測定方法で信頼性ある測定ができているとい える． このように，提案手法では，検者内信頼性，検者内信頼性共に信頼性のある測定が行えてい ると考えられる．

39

第

5 章 領域別フィルタの導入

第2 章第 4 節で高分解能化への要求について記述した． 本章では，領域別フィルタを導入し高分解能化するまでのアルゴリズムについて記載する．

5-1 要求分解能と高分解能化の手法

第2 章 4 節でも記載したが，筋の計測，そして腱，癌などの計測を行っていくために必要 であろう分解能をFig.5-1-1 に示す． Fig.5-1-1 要求分解能 現在のC-SWE の分解能は約 5[mm]程度であり，現在本研究室で開発中のリハビリテーシ ョンなどで使用できる小型のタブレットエコーの分解能は約8mm から 10ｍｍ程度である． Fig.5-1-1 より， アキレス腱などの腱の硬さ測定の要求分解能は約3-5mm， 痛み治療に使われるファシアハイドロリリースなどにおける要求分解能は1-2mm 針灸などによるエコーガイドでは，0.2-0.4mm の分解能が求められている． 現在の測定法での分解能では測定が難しいため，C-SWE の適用範囲の拡大するには高分解 能化が必要である．

40 高分解能にするためには，せん断波の周波数（加振周波数）を今の周波数よりも高い周波数 にすることが最も簡単な手法である．しかし，せん断波の周波数を高くすると，せん断波の 減衰が大きくなってしまい侵入深さが制限されてしまう． 高周波伝播の例として，半腱様筋で加振周波数を245.8[Hz]と 73.6[Hz]に設定し計測した例 をFig.-5-1-2 に示す． Fig.5-1-2 加振周波数でのせん断波伝播の違い Fig.5-1-2 より，加振周波数 245.8[Hz]よりも 73.6[Hz]の低い周波数のせん断波を伝播させ るほうがせん断波の伝播の減衰を抑えることができ，侵入深さを確保することができるこ とが確認できる． このことから，高分解能化にするためにせん断波の周波数を上げることは良い方法ではな いことがわかる． そこで，私は領域別フィルタを使った高分解能化への手法を提案する． 領域別フィルタを用いた提案手法は，「B モード画像」「せん断波伝播位相図」「せん断波伝 播方向図」等の先験的に得られる情報を使用して，臓器境界を算出．そのあと領域毎に個別 のフィルタを適用し伝播速度推定を行うことで，高分解能を実現する． 先験的情報としてなぜそれらが使えるのかを以下に記述する． ・B モード 筋膜で区切られた別々の領域に分割が可能である． ・せん断波伝播位相図，せん断波伝播方向図 せん断波の伝播特性の変化を利用し別々の領域に分割が可能である．

41

5-2 領域別フィルタのアルゴリズム

前節で，領域別フィルタを使用した高分解能化への手法を提案した． 本節では，領域別フィルタのアルゴリズムについて記述する． まず，領域別フィルタを用いた高分解能化への手法をFig.5-2-1 に示す． Fig.5-2-1 領域分割フィルタを用いた高分解能化への手法 Fig.5-2-1 に示したように提案手法は，２つのステップに分かれている． まず，ステップ１として，解析したいデータに対して，領域別フィルタにより境界を算出し 複数領域に分割を行う．ステップ2 として，領域分割を行った後，分割領域毎に波数ベクト ル上で適応的ウィナーフィルタを演算し,全領域にごく弱い低域通過フィルタをかけること により領域間の不連続性を取り除く処理を行う． 次に，領域別フィルタのアルゴリズムを記述する．

42

5-2-1 全体フロー

領域別フィルタのフローチャートをFig.5-2-1-1 に示す．

43 ここからはフローチャートに沿って説明していく． まず，先験的情報の選択から連結領域へのラベリングの部分について記述する． 解析したいデータに対して，先験的情報別に処理が異なる． B モード使用時にはｸﾞﾚｰｽｹｰﾙ化，背景減算を行い，位相図・方向図の場合は上下ピクセルの 差分を取ることで変化量を算出する． 次に，傾き抽出フィルタをかけることで，ノイズの低減と筋膜の強調を行う．傾き抽出フィ ルタの細かい概要については後述する． 次に適応的二値化を行う．適応的二値化は画像ごとに適応的に閾値を算出し二値化を行う． これによって，画像ごとに適した閾値で筋膜の描出を行うことができる．適応敵二値化につ いても細かい概要については後述する． 次に，連結領域へのラベリングを行い，抽出した筋膜に対してラベリングすることによって 区別する．連結領域へのラベリングについても後述する． フローに応じた画像をFig.5-2-1-2 に示す．

44 Fig.5-2-1-2 フローに応じた画像 次に，連結領域の検出の後から固定閾値二値化までの部分について記述する． 連結領域へのラベリングの後に，ラベルの長さの計算を行い，算出したラベルの長さをもと に，細かい筋繊維の削除や描出の補足を行う． Fig.5-2-1-2 の連結領域へのラベリング画像を見てわかるように，二値化して描出した筋膜 には，細かい筋繊維などの雑音や描出できなかった筋膜が混在している． これにたいして，ラベル毎の長さを算出するために，筋膜がどうかの判断を行う．

45 閾値としてROI の横幅の長さを設定し， ラベルの長さが閾値と同じ、つまり筋膜である場合はなにも行わない． 閾値の半分以上，つまりROI の半分の長さ以上のものは，筋膜とみなし直線で補足． 閾値の半分未満は，細かい筋繊維とみなして削除．の処理を行う． Fig.5-2-1-3 に計算処理を行った後の画像を記載する． Fig.5-2-1-3 ラベルの長さによる筋膜の判別結果 続いて，固定閾値による二値化から，連結領域へのラベリングについて記述する． 固定閾値による二値化を行い，筋膜で分割された筋を描出する．そして，その描出した筋に 対して，連結領域へのラベリングを行い，筋をラベルによって区別する． 固定閾値による二値化，連結領域へのラベリングを行ったあとの画像をFig.5-2-1-4 に示す． Fig.5-2-1-4 固定閾値による二値化，連結領域へのラベリング結果

46 Fig.5-2-1-4 はラベル毎に色分けして表示した図で，筋を筋膜で分割された領域毎にラベル 付けできていることがわかる．これによって，ラベル付けを行った領域別にフィルタをかけ ることが可能となる． ここまでが前節で示したstep1 の処理となる． つぎに，step2 の処理である，分割された領域に対してのフィルタの部分について記述する． Fig.5-2-1-1 のフローの「ラベル毎にフィルタ」「全領域にフィルタ」の部分に該当する． ラベル毎のフィルタでは，連結領域へのラベリングによりラベル付けを行った筋に対して， 個別に方向性フィルタやウィナーフィルタをかることで，クラッタノイズの除去，反射波成 分の除去，進行波成分の抽出を行う． ラベル毎のフィルタ（分割領域への個別のフィルタ）をかけるとこで，伝播図において領域 間に不連続性が生じてしまうため，これを取り除くために全領域に対してごく弱い低域通 過フィルタをかける． ラベル毎のフィルタ処理，全領域へのフィルタをかけた際の伝播図を Fig.5-2-1-5 に示す． Fig.5-2-1-5 ラベル毎のフィルタ，全領域へのフィルタでの伝播図 以上のフローチャートに沿って，領域分割を行い，筋膜を介した複数領域に分割を行うこと で，領域別フィルタによる高分解能化を図る．

47

5-2-2 各フィルタの概要

ここでは，前節で紹介した「傾き抽出フィルタ」，「二値化」，「連結領域へのラベリング」の 3 つについて記述していく． まず，傾き抽出フィルタについて記述する． 傾き抽出フィルタは，二次元波数ベクトル上で特定の範囲の傾きを持つ線分を抽出する フィルタである．今回は，-5°から 15°の線分を強調するように設定しているため，Fig.5-2-2-1 のようなフィルタとなる． Fig.5-2-2-1 -5°から 15°を通過域とするフィルタ これにより，細かい筋繊維の細分化や筋膜を残しながらノイズの低減を行うことが出来る． 以下に例を示す． 傾き抽出フィルタなしでB モードを二値化した結果について Fig.5-2-2-2 に示す． Fig.5-2-2-2 傾き抽出フィルタなしで二値化した場合

48 Fig.5-2-2-2 のように，B モード上で筋膜の輝度値が低く，明瞭に確認できない筋膜がある 場合，筋膜の一部が描出できない場合がある． 次に，傾き抽出フィルタを使用した場合の二値化結果をFig.5-2-2-3 に示す． Fig.5-2-2-3 傾き抽出フィルタを使用し二値化した場合 Fig.5-2-2-3 より，傾き抽出フィルタを使用することによって，筋膜部分が強調されて，輝 度値が低くて描出できなかった筋膜を描出できていることがわかる． 次に，二値化の手法について記述する． 二値化には，閾値を画像ごとに適応的に決定する「適応的二値化」を使用している． 例として，閾値を固定（閾値𝑡 = 75.0）で二値化したものを Fig.5-2-2-4 に示す。 Fig.5-2-2-4 閾値固定（𝑡 = 75.0）で二値化

49 画像３は閾値が低いため，筋膜下部の細かい筋繊維を検出してしまう．一方で、画像４は閾 値が高いため，検出したい筋膜が途切れて検出してしまう． このように，すべての画像に対して同じ閾値で二値化してしまうと，検出できる筋膜にば らつきが生じてしまい，検出したいものが検出できなかったり，検出したくないものが検出 できなかったりしてしまう． よって，Fig.5-2-2-5 に示すようなアルゴリズムで画像ごとに閾値を決定する適応的二値 化を行っている． Fig.5-2-2-5 二値化決定フローチャート Fig.5-2-2-5 のアルゴリズムにより決めた閾値を使って二値化（適応的二値化）を行った 画像をFig.5-2-2-6 に示す． Fig.5-2-2-6 適応的二値化 Fig.5-2-3-6 のように，検出したくないノイズ成分は検出せず，検出したい筋膜を検出で きていることがわかる．これにより，幅広い B モード画像で有効的に二値化を行うことが

50 出来る． 次に，連結領域の抽出について記述する． 連結領域の抽出は，ある閾値以上の値を持っていてつながっている領域に対しラベル付け を行うことで，複数領域に分離し領域分けを行う手法である． これにより，ラベルを使った分類やラベル毎の処理が可能となる． 例として，Fig.5-2-2-7 の二値化画像に対して，連結領域へのラベリングを行ったものを Fig.5-2-2-8 に示す． Fig.5-2-2-7 連結領域抽出前画像 Fig.5-2-2-8 連結領域の抽出 Fig.5-2-2-8 のように，筋膜で分割された領域にラベルを振り分けることで，筋を分割する ことが出来る． このラベリングにより，分割された領域に個別にフィルタをかけることが可能となる．

51

5-3 先験的情報別の分割結果

本章で，領域分割をするにあたり使用する先験的情報として， 「B モード画像」，「せん断波位相図」，「せん断波方向図」の三つを提案した． 本節では，三つの情報で領域分割を行った結果を示す． 領域分割を行うにあたり，三つのデータを例として使用した． 分割結果をFig.5-3-1 に示す． Fig.5-3-1 先験的情報別の分割結果 Fig.5-3-1 より， B モード画像を利用した場合の領域分割は，筋膜が明瞭に確認できる場合において最も有 効であることを確認した． 位相図，方向図を利用した場合の領域分割は，筋膜による位相ズレがある場合に有効であ り，一様に伝播している組織に対しては適用が困難であることがわかる． 僧帽筋などの骨格筋を測定する場合は，B モードに筋膜が明瞭に確認できるため，B モー ド画像を利用した領域分割が有効であると考えられる

52

5-4 分割結果まとめ

今回提案した先験的情報を用いた領域分割を， 群馬大学医学部付属病院整形外科で2018 年 9 月から 2019 年 1 月に測定したデータに対 して適用した．結果をFig.5-3-1，Fig.5-3-2，Fig.5-3-3 に示す． Fig.5-3-1 領域分割フィルタ適用結果 適用するにあたり， 最初に，骨格筋に対して最も有効であると考えられるB モードを利用した領域分割を適用 し，分割できなかった例に対して位相・方向図を利用した領域分割を適用した． 結果を三つのパターンとして示す． パターンA：B モードを利用した領域分割が有効であった．(Fig.5-3-2) パターンB：B モードを利用した領域分割が有効でなく，位相・方向図を利用した領域分 割が有効であった．(Fig.5-3-3) パターンC：領域分割が有効でなかった．

53 Fig.5-3-2 パターン A 結果 Fig.5-3-3 パターン B 結果 全117 データの分割結果として， パターンA が 117 例中 113 例，パターン B が 117 例中 4 例，パターン C は 117 例中 0 例 となった． このことから，全117 データについて領域分割を用いた高分解能化は適用可能であること がわかる．

54

第６章 領域別フィルタによる分解能評価

前章で領域分割フィルタを用いた分割の結果について記述した． 本章では，領域別フィルタを用いた高分解能C-SWE 法の分解能の評価を行ったので報告す る．

6-1 シミュレーションによる分解能評価

領域別フィルタを用いた高分解能C-SWE 法での分解能を評価するためにシミュレーシ ョンを行った． 汎用の超音波映像装置をモデルとして，1 ピクセルの大きさを 0.13[mm]，ROI を縦 185 ピクセル，横170 ピクセルとした．ROI 内には速度の異なる上下 2 層の物体を仮定し，2 層 間の速度遷移領域の幅から分解能を求める．分解能は，伝播速度の異なる二層の物体を測定 したと仮定して，z 方向の速度変化をとったときに，推定した速度の誤差が 10％に収まる 範囲を分解能として定義した． 推定した速度の誤差が 10%以内に収まる幅を分解能 R として定義した．分解能の定義を Fig.6-1-1 に示す． Fig.6-1-1 分解能の定義 仮定した物体(Fig.6-1-2)は，上の層の伝播速度を 2.0[m/s]，下の層の伝播速度を 3.0[m/s] とした．

55 Fig.6-1-2 仮定した物体 この物体に対して，従来法シミュレーションにて出力した速度図と伝播図をFig.6-1-3 に 示す． Fig.6-1-3 従来法での速度図と伝播図 また，仮定した物体の中心のラインで取得した速度変化を表すグラフをFig.6-1-4 に示 す． Fig.6-1-4 従来法での速度変化

56 Fig.6-1-4 より，従来法での分解能 R は 3.68[mm]と算出された． 次に，提案手法である領域別フィルタを用いた高分解能C-SWE によるシミュレーション 結果を示す．速度図と伝播図をFig.6-1-5 に，同じラインでの速度変化を Fig.6-1-6 に示 す． Fig.6-1-5 提案手法での速度図と伝播図 Fig.6-1-6 提案手法での速度変化 Fig.6-1-6 より，提案手法での分解能 R は 2.11[mm]と算出された． このことから，シミュレーションにおいて 従来法の分解能：3.68[mm]，提案手法の分解能：2.11[mm]となり，提案手法では分解能 を従来法の約1.7 倍程度に向上できることが確認できた．

57

6-2 二層ファントム測定による分解能評価

次に，弾性の異なる二層ファントムを作成し，それを用いて実際に分解能を評価した． ファントムは，上の層の伝播速度が約3.0[m/s]，下の層の伝播速度が約 3.4[m/s]のファント ムを作成した． この二層ファントムについて，従来法で解析を行った時の分解能と提案手法で解析を行っ た時の分解能を算出し比較を行った．解析には 5 測定分のデータを使用し，その平均値と 標準偏差を算出した．例として1 例紹介する． 測定した，二層ファントムのB モードを Fig.6-2-1 に示す． Fig.6-2-1 二層ファントム B モード 従来法で二層ファントムを測定したときのせん断波伝播図，速度図をFig.6-2-2 に示す． Fig.6-2-2 二層ファントム従来法結果 また，提案手法で解析したときのせん断波伝播図，速度図をFig.6-2-3 に示す．

58 二層ファントムでは，B モード上で二層の境界が筋膜のように明瞭に現れていないため，領 域分割には位相図を用いた． Fig.6-2-3 二層ファントム提案手法結果 次に，Fig.6-2-2，Fig.6-2-3 で示した速度図の黄点線ライン上（中央のライン）での速度変 化をFig.6-2-4，Fig.6-2-5 に示す． Fig.6-2-4 従来法での速度変化

59 Fig.6-2-5 提案手法での速度変化 Fig.6-2-4，Fig.6-2-5 より分解能をそれぞれ算出すると， 従来法3.81[mm] ，高分解能 C-SWE2.50[mm]となった． 5 測定分の平均値と標準編差を算出すると， 従来法：4.21±0.47[mm] 高分解能C-SWE：2.74±0.19[mm] となり，高分解能C-SWE では分解能を従来法の約 1.5 倍程度に向上できることが確認でき た．

60

6-3 分解能評価まとめ

本章では，シミュレーションによる分解能評価と，二層ファントムを用いた分解能評価を 行った． 二つの方法で，算出した分解能をFig.6-3-1 に示す．二層ファントムを測定して算出した 実測値は5 測定の平均値を表示している． Fig.6-3-1 シミュレーションと実測値による分解能 Fig.6-3-1 より，シミュレーション，実測値それぞれにおいて提案手法では分解能が向上し ていることが確認できる． 今回のシミュレーションではノイズを加えていないため，分解能は実測値の方が低く現れ ていると考えられる．

61

第

7 章 運動負荷を加えた際の僧帽筋硬度の経時的な変化計測

本章では，提案手法で僧帽筋に運動負荷を加えた際の経時的な変化を計測したのでその結 果について記述する．

7-1 実験系・実験プロトコル

[実験条件] 超音波映像装置 EUB7500(HITACHI) 加振周波数 76.0Hz 測定部位 僧帽筋 被検者 13 人 [実験方法] プローブの位置は肩甲骨の内上角と肩峰後角を結ぶ線分の垂直二等分を中心とする位置 とし，運動負荷前後で同じ位置にプローブを置くためにプローブをあてる位置を皮膚表面 に予めマーキングした.加振源の位置を運動負荷前後で一定とするために小型加振器を生体 表面にテーピングテープで貼り付けた.加振源の振動振幅はおよそ 170μm で加振した． 運動負荷は重さ 5kg のダンベルを利き腕で上げ下げするダンベルシュラッグ(Fig.7-1-1) とし，ダンベルシュラッグ10 回を 1 セットとして，5 セット（計 50 回）行った. 計測は，Fig.7-1-2 に示す測定プロトコルに従い，負荷前（control）と負荷中 1 セット毎, さらに負荷終了後10 分，30 分，60 分の 3 計測を加えた計 9 回行った. Fig.7-1-3 が実験の様子であり，計測は座位にて行い，利き腕を計測した. Fig.7-1-1 ダンベルシュラッグの概要

62

Fig.7-1-2 測定プロトコル

63

7-2 高分解能化によるせん断波伝播

次に，測定で得られたせん断波伝播図について考察した． 従来法で得られた伝播図と提案手法で得られた伝播図で比較を行った． まず，負荷によるせん断波伝播の変化について，負荷前と負荷10 回ごとに測定した伝播 図の内，従来法で解析したものをFig.7-2-1，提案手法で解析したものを Fig.7-2-2 に示 す． Fig.7-2-1 負荷中被験者 A 従来法でのせん断波伝播図 Fig.7-2-2 負荷中被験者 A 提案手法でのせん断波伝播図 Fig.7-2-1，Fig.7-2-2 より 従来法では確認できなかった筋膜による位相ズレが，提案手法では確認できることがわか る． 負荷前では筋膜部分での位相ズレが起きてないことが確認でき，10 回，20 回負荷を行っ た後の伝播図では，従来法では位相ズレが確認できないが，提案手法では筋膜の部分で位 相ズレが起き始めていることが確認できる．このことより，運動により筋膜に変化が起き ていることが高分解能化することにより確認できる． 続いて，30 回，40 回，50 回負荷を行ったときは，従来法では，若干の位相ズレが確認で き，提案手法では，位相ズレがより大きくなっていることが確認できる． また，提案手法では，50 回負荷を行ったときに筋膜を介した深浅の筋肉での速度の差を明

64 瞭に確認することができる． このような負荷による筋繊維での位相ズレは，被験者13 名中 7 名で確認できた．次に，被 験者別の変化の違いについて確認した． 被検者B での負荷前と負荷 10 回ごとに測定した伝播図を Fig.7-2-3 に示す． Fig.7-2-3 負荷中被験者 B 提案手法でのせん断波伝播図 Fig.7-2-2，Fig.7-2-3 より，どちらの被験者も，筋膜による位相ズレが確認できていること がわかる．被検者A と比べ被検者 B は，運動初期の段階から，筋膜による位相ズレが生じ ており，筋膜を介した深浅の筋肉の伝播速度の差が発生していることがわかる．また，被験 者B は被験者 A よりも筋膜での位相ズレが大きいことが確認できる． このように，負荷による変化は被験者毎にことなり，せん断波伝播図を使用することで筋膜 の癒着等を評価することが出来る可能性が示唆された． 次に，休憩によるせん断波伝播の変化について考察した． まず，被験者A における従来法での伝播図と提案手法での伝播図を Fig.7-2-1-4，Fig.7-2-5 に示す． Fig.7-2-1-4 休憩中被験者 A 従来法でのせん断波伝播図

65 Fig.7-2-5 休憩中被験者 A 提案手法でのせん断波伝播図 Fig.7-2-1-4，Fig.7-2-5 より，提案手法において，休憩後も筋膜による位相ズレが起きてい ることが確認できる．また，提案手法において負荷50 回の時の筋膜を介した深浅の筋肉で の伝播速度の違いが，休憩しているときも引き続き現れていることがわかる． このことから，1 時間の休憩では負荷による筋肉への変化は回復できていないことが推察で きる． 同様に，被験者B の休憩中の伝播図を Fig.7-2-6 に示す． Fig.7-2-6 休憩中被験者 B 提案手法でのせん断波伝播図 Fig.7-2-6 より， 被検者A と同様に，被験者 B も筋膜休憩中において筋膜での位相ズレが確認できる． 被検者A，B 共通の変化として，休憩中には筋膜を介した筋肉での伝播速度の違いが明瞭 に確認できる． このように領域別フィルタリングによる高分解能化によって，負荷 休憩による筋繊維の 微細な変化が確認できるようになった．

66

7-3 実験結果

本実験で得られた僧帽筋の測定データの例を Fig.7-3-1 に示す． Fig.7-3-1 僧帽筋での測定例 また，7-1 節で示した実験プロトコルに従って，僧帽筋運動負荷時の経時的変化をまと めたものをFig.7-3-2 に示す． Fig.7-3-2 僧帽筋運動負荷時の経時的変化 Fig.7-3-2 から，運動負荷を加えることによって，せん断波の伝播速度が変化しているこ とを確認した． 次に，得られた結果について考察を行う．

67

7-3-1 運動負荷による伝播速度変化量

ここでは，トレーニングによるせん断波の伝播速度変化量，すなわち筋弾性変化のしやす さについて考察を行う． トレーニングする前の伝播速度𝑉𝑐と、トレーニング中の伝播速度変化の最大値𝑉𝑚𝑎𝑥，伝播 速度変化の最小値𝑉𝑚𝑖𝑛、そしてその差で表される∆Vに着目した。定義を Fig.7-3-1-1 に示す． Fig.7-3-1-1 ∆Vの定義 すべての被験者について，𝑉𝑐と∆Vの関係を図示したものを Fig.7-3-1-2 に示す． 図示するにあたり，横軸に運動前の伝播速度𝑉𝑐をとり，縦軸に変化量∆Vをとった． つまり，Fig7-3-1-2 の右側は運動前の僧帽筋の硬さが硬く，左側は柔らかい． 上側は負荷によって筋弾性変化しやすく，下側は変化しにくい． と考えられる．