C-SWE法による痙縮腕の定量評価

令和元年度 修 士 論 文

C-SWE 法による痙縮腕の定量評価

指導教員 山越 芳樹 教授

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

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

堀口 悠希

1 C-SWE 法による痙縮腕の定量評価 目次 第１章 序論 第２章 せん断波計測について 2-1 せん断波とは 2-2 生体内組織における低周波振動の伝搬 2-3 せん断波計測で期待されるパラメータと臨床的有用性 第３章 実時間せん断波映像法の原理 3-1 超音波パルスドプラ法による組織内振動伝搬計測 3-2 カラーフロー映像系（CFI）の流速推定アルゴリズム 3-3 CFI の流速推定アルゴリズムによるせん断波の波面検出 3-4 定量的なせん断波画像の構成法 3-5 波数ベクトルフィルタリング 第４章 痙縮筋の弾性計測の現状と課題 4-1 測定装置の現状 4-2 せん断波の異方性による測定の影響 4-3 せん断波の加振方向と実際に見えるせん断波について 第５章 2 次元 TE_FDTD 法を波面伝播に適用 第６章 FDTD 法による波面伝播評価と基礎実験 6-1 せん断波伝播方向に影響する波動現象 6-2 波面伝播の評価方法 6-3 第１層の有無による伝播への影響 6-4 加振源―測定 ROI 間の距離の変化による伝播への影響 6-5 ファントムでの検証実験 第７章 実時間せん断波映像法を適用した痙縮腕の評価 7-1 測定概要 7-2 痙縮・非痙縮の上腕二頭筋測定例 7-3 外れ値除去を行った ROI 内の最大値で評価 7-4 せん断波の伝播方向のバラつきの評価 第８章 結論 8-1 結論

2 8-2 今後の課題

3 第1章

序論

今日の日本は社会の高齢化が進んでいる．平均寿命が90 歳まで延びることが想定されて いる中，自立した生活を送れる期間「健康寿命」が寿命にどれだけ近づくことができるかで 人生の豊かさが大きく左右される．健康寿命を延ばす要素は様々であるが，身体的自由を保 障するために筋疾患が大きく関わる．手足が動かしにくくなったり，勝手に動いてしまう痙 縮は脳梗塞，脳出血，くも膜下出血の後遺症で発症することがある．これらの疾患は合わせ て年間約29 万人が発症する[1]．関節拘縮による運動機能障害を引き起こす筋ジストロフィ ーで最も発症率が高いデュシンヌ型筋ジストロフィーは10 万人に 6 人の出生男児が発症す る[2]．基礎疾患である糖尿病から神経障害を併発し筋委縮が起こることがある．「糖尿病が 強く疑われる者」は日本に1000 万人いると推計されており[3]，筋肉が動かなくなることに よる運動機能障害は一般的になっていると言える．筋委縮は治療薬が開発されたり，その他 の運動機能障害においてもリハビリテーションによる症状緩和を実施しているが，症状の 度合いの定量的な測定は臨床で普及されていない．筋弾性の定量評価として，パルス状のせ ん断波の組織内を伝播する伝播速度を測定する ARFI 法が普及されているが，装置が大型 で携帯不可，高価であることや測定できる伝播速度に上限があることが筋疾患を測定する 際に利便性が低いと考えられる．山越芳樹教授考案の Continuous-SWE 法（C-SWE 法） は組織を伝播する連続せん断波を可視化して組織の弾性を評価することができる．周波数 76 [Hz]程度の低周波振動を，生体組織などの比較的柔らかい物体の表面に加えると，その 放射エネルギーの大部分は物体中を横波として伝播する. そして，その伝播速度や減衰係数 は，せん断粘弾性パラメータと関連していることが知られている．また，生体組織のせん断 粘弾性特性は，生体組織を触った時の硬さや感触と密接に関係しているため，生体組織につ いて低周波振動での伝播速度や減衰などが測定できれば，疾病の進行度の定量的な評価や 早期発見が期待でき，これらは組織の特性化のために有用である． 本研究では，痙縮を罹患している患者を定量的に測定できるシステムを構築することを 目標としている．痙縮は，脳梗塞や交通事故などの後遺症で体を動かす神経が麻痺して体を 動かすことができなくなり，筋肉が萎縮して硬化する症状である．現在，痙縮を回復する方 法は確立されておらず，薬剤投与やリハビリテーションでの経過観察が一般的であるが，筋 硬度を測定する手法を導入していないことが多く，効果を定量的に評価できていない．痙縮 筋の硬度測定では，症状の度合いが様々なことから患者による測定値のばらつきが大きく， 現在普及されている測定技術では信頼性を確保できないことが理由として挙げられる．そ こで C-SWE 法のリアルタイム性，せん断波の可視化，装置が廉価で携帯可である特徴か ら，症状の度合いを定量的に測定する方法が確立することで，薬剤投与やリハビリテーショ ンの効果を確認することができ，今後の治療方法の発展が期待できる． 本研究では，痙縮筋のターゲットを上腕二頭筋とした．患者ごとに測定値がばらつく原因 を，FDTD(時間領域差分法)を用いた波面伝播シミュレーションから考察をしたうえで痙縮 患者を測定し，痙縮腕と非痙縮腕で有意差が認められるか検討した．

4

第2章 せん断波計測について

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

5 2-2 生体内部の組織における低周波振動の伝播 生体組織の粘弾性パラメータと低周波振動の伝播速度および減衰の関係について以下に 示す． 外部から媒質に振動を伝えると，その振動は一般的に縦波・横波として伝播する．生体の ような粘弾性媒質中では，Hooke の法則が成り立つ Voigt モデルと仮定することにより， 組織の粘弾性を推定する手法が提案[4]されている．この縦波・横波の伝播速度および減衰 係数は次式で与えられる． ① 縦波 伝播速度 ： 𝑣𝑙 = 𝜔𝑣 𝑅𝑒[𝑔] (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]程度以下の低周波振動であると，外部から与えら れた振動のエネルギーはそのほとんどが横波に変換されると考えられている[5]．

6 ここで，(2-2-4)式，(2-2-5)式で与えられる横波の伝播速度と減衰係数を，粘弾性パラメー タを用いて書くと， 𝑣𝑡= √ 2(𝜇12+𝜔𝑣2𝜇22) ρ(𝜇₁+√𝜇₁2_+𝜔 𝑣2𝜇22) (2-2-7) 𝛼𝑡= √ 𝜌𝜔_𝑣2_(𝜇 1+√𝜇12+𝜔𝑣2𝜇22) 2(𝜇12+𝜔𝑣2𝜇22) (2-2-8) となる． したがって，もし，媒質の弾性が粘性にまさり，𝜇1≫ ω𝑣𝜇1の関係が成り立つときには， 𝑣𝑡1≅ √ 𝜇1 𝜌 (2-2-9) 𝛼𝑡1≅ 0 (2-2-10) となり，伝播速度は，単にせん断弾性係数と媒質の密度のみの関数となる．このとき，𝜇1 が大きいということは，媒質が硬いということであり，硬い媒質ほど伝播速度は速くな る． 一方，媒質の粘性が弾性にまさり𝜇1≪ ω𝑣𝜇1の関係が成り立つときには， 𝑣𝑡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/ se c] 加振周波数[Hz] 弾性率2.26kPa, 粘性率2.38Pa・s 弾性のみの場 合（2.26kPa）

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

8

1. 実時間せん断波映像法の原理

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

ξ(𝑡) = 𝜉

₀

𝑠𝑖𝑛(2π𝑓

_𝑏

𝑡 + 𝜙

_𝑏

)

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

𝜙

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

𝑓 =

𝑐+𝑣(𝑡) 𝑐

𝑓

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

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

𝑓

′

₌

𝑐 𝑐−𝑣(𝑡)

𝑓

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

𝑓

′

=

𝑐 𝑐−𝑣(𝑡)

×

𝑐+𝑣(𝑡) 𝑐

𝑓

0

=

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

𝑓

0

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

∆𝑓 = 𝑓

′

− 𝑓

₀

=

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

𝑓

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π𝑓

0

𝑡 + ∆𝜙 − 2𝑘

𝑢

𝑍)

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

₀

𝑡 +

4𝜋𝑓

0

𝑐

𝜉(𝑡) − 2𝑘

𝑢

𝑍)

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

0

(𝑡 + 2

𝜉(𝑡) 𝑐

) − 2𝑘

𝑢

𝑍}

(3-1-8) ただし， 𝐴(𝑡)：振幅 𝑘𝑢 ：超音波パルスの波数 𝑍 ：トランスデューサ，散乱体間の距離 となる．よって超音波パルス間で微小変位𝜉(∆𝑡)による位相ずれが生じる．

10 次にRF 信号に，位相が互いに 90 度異なる超音波周波数成分を畳み込み積分し低域通 過フィルターをかけ，QI 信号を得る． (ⅰ) I 信号 RF 信号にキャリア信号を乗算すると

𝐼

′

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

₀

(𝑡 + 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 信号を得る．

11 3-2 カラーフロー映像系(CFI)の流速推定アルゴリズム いま，超音波パルスを同一方向にN パルス送波すると，i 番目の超音波パルスに対する受 信超音波の位相𝜙𝑖は， 𝜙𝑖= 𝜙0+ 2𝜋𝑓₀ 𝑐 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𝜋𝑓0 𝑐 2𝑣 𝑖 Δ𝑡)) (3-2-4) ここで，第i 番目の超音波パルスの位相と，第 i+1 番目の超音波パルスの位相の差Δ𝜙𝑖を考 える．これは， Δ𝜙𝑖= 𝑎𝑟𝑔 (Q⃗⃗ 𝑖+1Q⃗⃗ 𝑖 ∗ ) (3-2-5) と推定できるので，(3-2-4)式を代入すると，

12 Δ𝜙𝑖= 𝑎𝑟𝑔 ( 𝑎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

13 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) が満たされるとする．

14 上記条件（(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𝜋𝑓0 𝑐 𝜉) 𝑄𝑖= 𝑎 sin ( 4𝜋𝑓0 𝑐 𝜉) (3-3-9) となる．

𝛥𝑡

2𝛥𝑡

𝜉

3𝛥𝑡

4𝛥𝑡

0

𝑡

15 ここで，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) ② 𝜆₈ ≤ 𝜉0≤ 3𝜆 8 の場合 { 𝐼𝑖≤ 0 𝑄𝑖≥ 0 i = 2の場合 {𝐼𝑖 = 𝑎 𝑄𝑖= 0 (3-3-13) i = 3 の場合 𝐼𝑖= 𝑎 cos ( 4𝜋𝑓0 𝑐 𝜉0) (3-3-14) ただし， ① 0 ≤ 𝜉0≤ 𝜆 8 の場合 { 𝐼𝑖 ≥ 0 𝑄𝑖≤ 0 (3-3-15) ② 𝜆₈ ≤ 𝜉0≤ 3𝜆 8 の場合 { 𝐼𝑖≤ 0 𝑄𝑖≤ 0 となる．

16 (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 _𝟖𝝀 ≤ 𝝃𝟎≤ 𝟑𝝀 𝟖 での直交検波器の出力信号

17 これらを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 流速導出の基本演算 ここで超音波パルスの送受信数 N=11 の場合に，CFI による流速導出アルゴリズムを図式 化すると

18

Fig.3-3-5 CFI における流速導出アルゴリズム

(3-3-4),(3-3-6)式の 2 つの条件がともに満たされているとき，CFI における流速推定は Fig.3-3-6 のようになる．

19 ここで， { 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 に示す． [シミュレーション条件] 超音波中心周波数 𝑓0 6.5𝑀𝐻𝑧 超音波伝播速度 𝑐 _{1500 𝑚 𝑠}⁄ パルス繰り返し周波数 1 𝑑𝑡⁄ 365𝐻𝑧 パルス本数 𝑁 11 加振周波数 𝑓𝑏 91.25𝐻𝑧

20 Fig.3-3-7 数値シミュレーション結果 周波数条件は，理論的には（3-3-7）式で表されるが，実際には，せん断波の周波数がこ の条件に近いときでも，流速の最大値または流速0 の部分が CFI 画像上に現れる．そのた め，せん断波の周波数が周波数条件に近いときにも，せん断波の波面が再現できる． 𝜙𝑏

21 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) 𝜃𝑒𝑠𝑡= tan−1( 𝑘𝑧(𝑥,𝑧) 𝑘𝑥(𝑥,𝑧)) (3-4-6) ここで，x 軸方向の波数𝑘𝑥を実際の伝播図の見ばえに合った値に補正し，𝑣𝑏を計算する． 始めに𝑘𝑥から補正前の伝播速度𝑣𝑜𝑙𝑑を計算する． 𝑣𝑜𝑙𝑑(𝑥, 𝑧) = 2𝜋𝑓𝑏 √𝑘𝑥2+𝑘𝑧2 (3-4-7) 𝑣𝑛𝑒𝑤(𝑥, 𝑧) = 𝛼𝑣𝑜𝑙𝑑(𝑥, 𝑧) 𝛼：補正値 (3-4-8) 𝑘𝑥′(𝑥, 𝑧) = 2𝜋𝑓𝑏 𝑣𝑛𝑒𝑤(𝑥,𝑧) (3-4-9) よって，せん断波の伝播速度𝑣𝑏は 𝑣𝑏(𝑥, 𝑧) = 𝜃𝑒𝑠𝑡(𝑥, 𝑧) × 2𝜋𝑓𝑏 𝑘𝑥′(𝑥,𝑧) (3-4-10)

22 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)

23

第

4 章 痙縮筋の弾性計測の現状と課題

4-1 測定装置の現状

Fig4-1-1 に ARFI 法（Acuson S3000 HELX(Siemens)）で健常者と痙縮患者の上腕二 頭筋を測定したときの平均伝播速度を示す[6]．この実験では健常者７人，痙縮患者８人で 測定を行っている．痙縮患者の非痙縮側の伝播速度は健常者のもの差が小さい．しかし，痙 縮患者の痙縮側の伝播速度はばらつきが大きく，非痙縮側や健常者のものと分布領域が被 っている．このことから，現状の痙縮筋の測定では被検者による測定値のばらつきが大きく， 測定値の信頼性を確保できないと言える．また，従来のSWE 装置では Time of Flight 法で 計測しているため，測定できるせん断波の伝播速度の最大値に上限がある．健常者の筋肉よ りも硬いことが想定される筋疾患の筋肉を測定する場合，せん断波の伝播速度を上限以上 に表示することができず，標準偏差が小さくなり測定データの精度が低下することが考え られる．フレームレートを上げて測定できるせん断波の伝播速度の上限を上げた場合は，B モード画像の画質が低下してしまう．また，いずれの装置も大型で非常に高価であることか ら，筋疾患のリハビリテーションの現場に実装する際に不便である． 本研究で使用するC-SWE は測定できるせん断波の伝播速度に上限がない．また，超音波 装置を改良する必要がなく，測定環境を整えるための費用が廉価である．また，タブレット エコー装置を使用すれば携帯も容易であり，臨床の現場で実装が現実的であると言える． Fig4-1-1 健常者と痙縮患者の上腕二頭筋を測定したときの平均伝播速度

24 4-2 せん断波の異方性による測定の影響 Fig4-2-1 に大腿直筋の測定位置と B モード画像,Fig4-2-2 に大腿直筋の測定位置とプローブ の向きによるせん断波伝播速度の変化について示す[7]．Fig4-2-2 からわかるように，プロ ーブを配置する向きによってせん断波の伝播速度は50％程度異なる場合がある．構造が均 一な骨格筋で，せん断波が伝播する方向によって伝播速度が変化する異方性があるため，疾 患がある骨格筋を測定するときも異方性によって複雑に伝播することが考えられる．その ため，疾患を持つ筋肉に対して高精度な弾性計測を実現するためには，せん断波が筋繊維方 向にまっすぐ伝播する条件を探すことが重要である． Fig4-2-1 大腿直筋の測定位置と B モード画像 Fig4-2-2 大腿直筋の測定位置とプローブの向きによるせん断波伝播速度の変化

25 4-3 せん断波の加振方向と実際に見えるせん断波について Fig4-3-1 にせん断波の加振方向と実際に見えるせん断波の関係について示す． 加振器は接している組織に対して浅部から深部の方向に球面波を励起する．それにも関わ らず，実際に観測される生体を伝播するせん断波は，筋繊維に沿った（プローブ方向と一 致した）平面波に近い波面である．この現象について検証を行ったことが今までなく，原 因が明かされていない状態であった． 4-2 の「せん断波が筋繊維方向にまっすぐ伝播する条件」，4-3 の「加振器は接している組 織に対して浅部から深部の方向に球面波を励起するにも関わらず，実際に観測される生体 を伝播するせん断波は，筋繊維に沿った平面波に近い波面である原因」を考察するために FD-TD 法を用いたせん断波伝播評価を行った． Fig4-3-1 せん断波の加振方向と実際に見えるせん断波の関係

26

第

5 章

2 次元 TE_FDTD 法を波面伝播に適用

波源も物体の構造もｙ軸方向で変化のない状態であると仮定した2 次元 TE_FDTD 法では， 電磁界は以下の様になる． 𝑬𝒙𝑛(𝑖 + 1 2, 𝑘) = 𝐶𝑥(𝑖 + 1 2, 𝑘) 𝑬𝒙 𝑛−1 (𝑖 +1₂, 𝑘) +𝐶𝑥𝑦(𝑖 + 1 2, 𝑘) {𝑯𝒚 𝑛−1₂ (𝑖 +1₂, 𝑘 +1 2) − 𝑯𝒚 𝑛−1₂ (𝑖 +1₂, 𝑘 −1 2)} (5-1a) 𝑬𝒛𝑛(𝑖, 𝑘 + 1 2) = 𝐶𝑧(𝑖, 𝑘 + 1 2) 𝑬𝒛 𝑛−1 (𝑖, 𝑘 +1₂) −𝐶𝑧𝑦(𝑖, 𝑘 + 1 2) {𝑯𝒚 𝑛−1₂ (𝑖 +1₂, 𝑘 +1₂) − 𝑯𝒚𝑛− 1 2(𝑖 −1 2, 𝑘 + 1 2)} (5-1b) ただし 𝐶𝑥(𝑖 + 1 2, 𝑘) = 1 −𝜎 (𝑖 + 1 2 , 𝑘) ∆𝑡 2𝜀 (𝑖 +1_{2 , 𝑘)} 1 +𝜎 (𝑖 + 1 2 , 𝑘) ∆𝑡 2𝜀 (𝑖 +1_{2 , 𝑘)} 𝐶𝑥𝑦(𝑖 + 1 2, 𝑘) = ∆𝑡 𝜀 (𝑖 +1_{2 , 𝑘)} 1 +𝜎 (𝑖 + 1 2 , 𝑘) ∆𝑡 2𝜀 (𝑖 +1_{2 , 𝑘)} 1 ∆𝑧 𝐶𝑧(𝑖, 𝑘 + 1 2) = 1 −𝜎 (𝑖, 𝑘 + 1 2) ∆𝑡 2𝜀 (𝑖, 𝑘 +1₂₎ 1 +𝜎 (𝑖, 𝑘 + 1 2) ∆𝑡 2𝜀 (𝑖, 𝑘 +1₂₎ 𝐶𝑧𝑦(𝑖, 𝑘 + 1 2) = ∆𝑡 𝜀 (𝑖, 𝑘 +1₂₎ 1 +𝜎 (𝑖, 𝑘 + 1 2) ∆𝑡 2𝜀 (𝑖, 𝑘 +1₂₎ 1 ∆𝑥

27 である． 𝑯𝒚𝑛+ 1 2_{(𝑖 +}1 2, 𝑘 + 1 2) = 𝑯𝒚 𝑛−1₂ (𝑖 +1 2, 𝑘 + 1 2) −𝐶𝑦𝑥(𝑖 + 1 2, 𝑘 + 1 2) {𝑬𝒛 𝑛_{(𝑖 + 1, 𝑘 +}1 2) − 𝑬𝒛 𝑛_{(𝑖, 𝑘 −}1 2)} +𝐶𝑦𝑧(𝑖 + 1 2, 𝑘 + 1 2) {𝑬𝒙 𝑛 (𝑖 +1₂, 𝑘 + 1) − 𝑬𝒙𝑛(𝑖 + 1 2, 𝑘)} (5-2) ただし 𝐶𝑦𝑥(𝑖 + 1 2, 𝑘 + 1 2) = ∆𝑡 𝜇 (𝑖 +1_{2 , 𝑘 +}1₂₎ 1 ∆𝑥 𝐶𝑦𝑧(𝑖 + 1 2, 𝑘 + 1 2) = ∆𝑡 𝜇 (𝑖 +1_{2 , 𝑘 +}1₂₎ 1 ∆𝑧 である． せん断波波面伝播シミュレーションでは，磁界に該当するｙ軸方向に伝播するせん断波を 𝑉𝑦，電界に該当する成分を𝑃𝑥, 𝑃𝑧とする．加振源から発振されるせん断波を単一パルス波と して，以下の式を用いる． 𝑉𝑦= 30−(𝑓𝑡−1) 2 ∙ sin(2𝜋𝑓𝑡) (5-3) ここで，𝑓は周波数で，𝑓, 𝑡は時刻[𝑠𝑒𝑐] である． 微小ステップ時間∆t, ステップ数を nとするとt = n∆tとなる．また，単一セルを正方形とす ることで∆𝑑 = ∆𝑥 = ∆𝑧とする． 同一媒質中は弾性が均一でせん断波の伝播速度も均一であるとする．媒質1,2 の伝播速度 が𝑣1, 𝑣2(𝑣1> 𝑣2)のとき， 𝑣𝑝𝑙= 𝑣1∙ ∆𝑡 ∆𝑑 (5-4) cs = 1 + (𝑣1 𝑣2− 1) ∙ 𝑎 (5-5) 𝑄𝑠𝑝= 𝑉𝑝𝑙∙ (𝑐𝑠 ∙ 𝑣2 𝑣1) 2 (5-6) を定義する．式(5-25)の a は変数で(0, 1 11, 2 11, 3 11⋯ 11 11)という 12 個の数値のどれかが入る． これは，媒質1,2 の境界の伝播速度を𝑣1から𝑣2に緩やかに変化させるために用いる

28 （例 媒質1 の伝播速度を適用するときは a=0,媒質 2 の伝播速度を適用するときは a=11 11）． とする．式(5-19)(5-22)に(5-24)(5-25)(5-26)を代入すると 𝑷𝒙𝑛(𝑖 + 1 2, 𝑘) = 𝑷𝒙 𝑛−1 (𝑖 +1₂, 𝑘) +𝑄𝑠𝑝{𝑽𝒚𝑛− 1 2(𝑖 +1 2, 𝑘 + 1 2) − 𝑽𝒚 𝑛−1 2(𝑖 +1 2, 𝑘 − 1 2)} (5-7a) 𝑷𝒛𝑛(𝑖, 𝑘 + 1 2) = 𝑷𝒛 𝑛−1 (𝑖, 𝑘 +1₂) −𝑄𝑠𝑝{𝑽𝒚𝑛− 1 2(𝑖 +1 2, 𝑘 + 1 2) − 𝑽𝒚 𝑛−1 2(𝑖 −1 2, 𝑘 + 1 2)} (5-7b) 𝑽𝒚𝑛+ 1 2(𝑖 +1 2, 𝑘 + 1 2) = 𝑽𝒚 𝑛−1₂ (𝑖 +1₂, 𝑘 +1₂) −𝑣𝑝𝑙{𝑷𝒛𝑛(𝑖 + 1, 𝑘 + 1 2) − 𝑷𝒛 𝑛 (𝑖, 𝑘 −1₂)} +𝑣𝑝𝑙{𝑷𝒙𝑛(𝑖 + 1 2, 𝑘 + 1) − 𝑷𝒙 𝑛_{(𝑖 +}1 2, 𝑘)} (5-8)

29

第

6 章 FD-TD 法による波面伝播評価と基礎実験

6-1 せん断波伝播方向に影響する波動現象 せん断波の伝播方向に影響する波動現象を列挙し，この現象によってせん断波の伝播方向 が変化する条件を探索する． 1. 波の伝播距離（円の半径）による曲率の変化 Fig6-1-1 に球面波の半径による波の伝播方向の変化のイメージを示す．測定 ROI を模した 赤い四角の中に球面波があるとき，ROI の辺に接している 2 点を通る円の半径が曲率半径 R となり，これは球面波の半径と同じである．ROI 内の曲線の曲がり具合を示す曲率 k は k =1 𝑅で表すことができる．Fig6-1-1 で，𝑅1> 𝑅2のとき，それぞれの曲率𝑘1, 𝑘2は 𝑘1> 𝑘2が成立し，球面波の伝播距離が長くなるほどROI 内で観測できる波の曲率は小さく なると言える．これを生体のせん断波伝播の状況に用いると「せん断波の伝播距離（加振源 ―ROI 間距離）が長くなるほど観測できるせん断波の曲率が小さくなる」と言える．これを 検証する伝播評価を6-3 に示す． Fig6-1-1 球面波の半径による波の伝播方向の変化のイメージ 2. 波の屈折 Fig6-1-2 に 2 層構造のアルギン酸カルシウムファントムを測定したときのせん断波の屈折 の様子，Fig6-1-3 に波の屈折のイメージ，入射媒質での波の伝播速度を𝑣1，入射角を𝜃1，屈 折媒質での波の伝播速度を𝑣2，屈折角を𝜃2とすると，スネルの法則は， sin 𝜃1 sin 𝜃2= 𝑣1 𝑣2 （6-1） が成立する．これをsin 𝜃2について整理すると， sin 𝜃2= 𝑣2 𝑣1sin 𝜃1 （6-2）

30 となる．Fig6-1-4 に式（6-2）を𝑣2 𝑣1ごとに計算した結果を示す．Fig6-1-4 から， 𝑣2 𝑣1> 1.0を満 たしていれば波の屈折によって𝜃2> 𝜃1が成立し，屈折後の波は x 軸方向に対して水平にな ると言える．また，𝜃1が一定の場合， 𝑣2 𝑣1が大きくなるほど𝜃2 が大きくなる．これを生体のせ ん断波伝播の状況に用いると「脂肪層の下の筋肉を測定するとき，せん断波の屈折によって， 筋肉層を伝播するせん断波は筋繊維に沿って水平に伝播しやすくなる」と言える．波の屈折 現象を利用した波面伝播評価を6-4 に示す． Fig6-1-2 2 層構造のアルギン酸カルシウムファントムで観測したせん断波の屈折 Fig6-1-3 に波の屈折のイメージ

31

32 6-2 波面伝播の評価方法 本研究では群馬大学三輪空司准教授が作成した、2 次元 FD-TD 法を用いた波面伝播シミュ レーションプログラムを使用した． せん断波が筋組織をモデリングしたシミュレーションを行うことで，せん断波が筋組織に 水平に伝播する条件を考察する． ある時刻 t での𝑉𝑦を log スケールの平面に表示する．せん断波伝播の変化を評価するため に，波面の角度を計算する．Fig6-1 に波面伝播シミュレーションの一例と波面の先端の角 度の計算の概要を示す．ある閾値の振幅の領域を「波面」として赤線で描画する．波面の赤 線の隣り合うプロット点を(𝑥1, 𝑧1), (𝑥1, 𝑧1)とすると，波面の先端の角度はθで表すことがで きるため，θは θ = tan−1₍𝑥2−𝑥1 𝑧2−𝑧1) (6-1) となる．超音波装置を想定した測定 ROI を配置して，ｚ軸方向 20mm 分のθの平均値を 「ROI 内角度𝜃𝑎𝑣」として，評価の指標とした． シミュレーションでは加振源からROI 間の距離 L，測定対象の層構造および第１層の伝播 速度を比較対象とした． Fig6-1 波面伝播シミュレーションの一例と波面先端の角度の計算

33 6-3 第１層の有無による伝播への影響 人体の骨格筋では，筋肉の上に皮膚や脂肪層があるため，筋肉層の上に別の層があること によって筋肉層を伝播するせん断波の波面に影響があるか検証した． 第１層が無い均一媒質を伝播するせん断波の伝播速度を3.9[m/s]としてシミュレーション した結果をFig6-2-2-1,に示す．第２層を伝播するせん断波の伝播速度を 3.9[m/s]，厚さ 8mm の第１層を伝播するせん断波の伝播速度を 1.0[m/s]，2.0[m/s]，3.0[m/s]としてシミ ュレーションを行った結果をFig6-2-2-2，Fig6-2-2-3 Fig6-2-2-4 に示す． Fig6-3-1 １層構造のシミュレーション結果 Fig6-3-2 2 層構造で第 1 層の伝播速度が 1.0[m/s]のシミュレーション結果

34

Fig6-3-3 2 層構造で第 1 層の伝播速度が 2.0[m/s]のシミュレーション結果

35

6-4 加振源―測定 ROI 間の距離の変化による伝播への影響

加振源からROI までの距離 L を変化させた場合のせん断波の角度への影響を検証した．

一様媒質を伝播するせん断波の伝播速度を 3.9[m/s]としてシミュレーションした結果を， Fig6-2-1, Fig 6-2-2, Fig 6-2-3, Fig 6-2-4 に示す．

Fig6-2-1 L=14.9mm の場合のシミュレーション結果

36 Fig6-2-3 L=36.3mm の場合のシミュレーション結果 Fig6-2-4 L=45.2mm の場合のシミュレーション結果 L が大きくなるにつれて𝜃𝑎𝑣が小さくなった． 波の屈折と伝播距離による曲率の変化の相乗効果で波の伝播方向が水平に近づくと考えら れる．6-3，6-4 のシミュレーション結果をまとめたグラフを Fig6-4-5 に示す．

37 Fig6-4-5 シミュレーション結果 以上より， 1. 曲率の変化により，L が約 10mm から 40mm まで変化すると𝜃𝑎𝑣は約50%になった 2. せん断波の屈折により，1 層構造に対して 2 層構造の𝜃𝑎𝑣が約50%になった 文献値によると，せん断波の伝播方向が90[deg]変化するとせん断波の伝播速度が最大で 1[m/s]程度変化すると言われている．骨格筋において，脂肪速度は約 2[m/s]と言われてい るため，第1 層が 2[m/s]のファントムでの𝜃𝑎𝑣の最小値，約5[deg]では伝播速度が 0.06[m/s]変化すると推測できる．波の曲率の変化と屈折の効果で，異方性が測定誤差とし て支障が無いレベルまで小さくなったと言える． 実際の測定でも同様の結果が得られるか，ファントムを用いて実験した．

38 6-5 ファントムでの検証実験 6-2,6-3 のシミュレーションから推察された，せん断波の波面伝播の傾向が実測定でも見ら れるか検証した． 1. 実験のために以下のアルギン酸カルシウムファントムを作成した． ⅰ濃度が7.56w%均一なファントム ⅱ上部約5mm が濃度 4.35w%,下部約 20mm が濃度 7.56w%のファントム 2. 各ファントムを C-SWE で測定したところ，ⅰの伝播速度が 3.3[m/s],ⅱの上部の伝播 速度が2.4[m/s]，下部の伝播速度が 2.9[m/s]であった 3. ⅰのファントムで，加振源-プローブ間距離（L）が 20mm，50mm の状態でそれぞれ 測定した． 4. ⅰ，ⅱのファントムで L を 20mm とし，同じ深さの位置に ROI を配置して測定を行 った． Fig6-4-1 に加振源の図を示す．Fig6-4-2 の a に加振源プローブ間距離 L が 20mm，b に L が50mm のときの実験風景を示す．加振器の加振源は加振器の端から 20mm の位置にあ るため，プローブから30mm 離れた位置に加振器を置くことで，加振源とプローブは 50mm 離れている状態である． Fig6-4-1 加振器内の加振源の配置

39

Fig6-4-2 実験風景

Fig6-4-3 に 1 層構造，L=20mm のときの測定結果，Fig6-4-4 に 1 層構造，L=50mm のと きの測定結果を示す．

40

Fig6-4-4 1 層構造，L=50mm のときの測定結果

L が大きくなることによってせん断波の角度が小さくなり，シミュレーションと同様の結 果が得られた．

41 Fig6-4-5 2 層構造，L=20mm のとき測定した結果 Fig6-4-3 と比較して，2 層構造，L=20mm のファントムの方がせん断波の角度が小さくな り，シミュレーションと同様の傾向が得られた． Fig6-4-7 に同一ファントムで 5 回測定した結果を示す． 左のグラフが加振源―プローブ間距離で比較したもの．右が第1 層の有無を比較したグラ フである．

42 Fig6-4-7 ファントム実験結果の 5 データ統計 Fig6-4-7 から， 1. L が 20mm から 50mm まで変化すると，せん断波の角度が約 50%になった． 2. L が 20mm から 50mm まで変化すると，せん断波の角度が約 50%になった． と言える．せん断波の曲率の変化と屈折の相乗効果でせん断波の伝播方向が水平に近づく と考えられる．

43

第

7 章 実時間せん断波映像法を適用した痙縮腕の評価

本実験は自治医科大学の倫理委員会の承認のもと行われた． 痙縮患者の上腕二頭筋の測定を行い，痙縮側と非痙縮側で有意差が得られるか実験した． 7-1 測定概要 [実験条件] 超音波映像装置 GE LOGIQ7 加振周波数 73.8Hz 測定部位 上腕二頭筋 検者 3 人 被検者 １人 [実験方法] 1. 患者はベッドに仰臥位になり，腕を伸ばした．(Fig7-1) 2.左腕の測定箇所を決定後，加振器をテーピングテープで固定 3.20 データ測定する． 4.3 人の検者で同様に測定を行う． 5.右腕でも 1~3 の手順で測定を行う． 患者左腕で痙縮の症状が見られた．以下では，痙縮腕を「患側」，非痙縮腕を「健側」と表 記する． Fig7-1 ダンベルシュラッグの概要

44 7-2 測定例 本実験で得られた上腕二頭筋の測定データの例を Fig7-2 に示す。 Fig7-2 健側と患側の伝播図での比較 健側はせん断波が水平に伝播していることに対して，患側ではせん断波が歪み，波面が複雑 である．このことから，患側は，筋の萎縮を起こしているが，筋肉が硬くなった部分と柔ら かいままの部分が複雑に分布していると考えられる． 従来は，ROI 内の平均伝播速度を弾性の指標としていたが，一様に伝播していない痙縮腕 を評価することに不適であると考える．次項では，複雑に伝播したせん断波を評価する方法 を提案する． 7-3 新評価方法の提案 患側は健側と比較して伝播図の波面傾きから定義される伝播方向のバラつきが大きかった 一方で，波面の間隔から求められる伝播速度の最大値は患側の方が大きかった．このことか ら，観測された画像の特徴の違いを明確に数値で表すために次の2 点の手法を導入した． 1. 各ピクセルの伝播速度の最大値𝑉𝑚𝑎𝑥で評価(7-3-1) 2. せん断波の伝播方向のバラつき𝜃𝑠𝑑で評価(7-3-2)

45 7-3-1 外れ値除去を行った ROI 内の最大値𝑉𝑚𝑎𝑥で評価 せん断波の伝播速度の最大値で評価をした場合，せん断波の波面に歪みがあるデータでは 最大値が特異的に高くなることがあった．そこで，四分位数を用いた外れ値除去を行ったあ と伝播速度の最大値を𝑉𝑚𝑎𝑥として評価を行った．四分位数を用いた外れ値除去の方法を Fig7-3-1-1 に示す。 1. 全ピクセルの伝播速度を数値の昇順に並べ替える． 2. 全データのうち 25%の位置にある数値を第一四分位数（𝑉𝑄1），75%の位置にある数値を 第三四分位数（𝑉𝑄3）からIQR を計算する． 3. 最大値の外れ値𝑉𝑜を求める．𝑉𝑜より大きい値は異常値とみなして除外し，除外後の最大 値を𝑉𝑚𝑎𝑥とする． Fig7-3-1-1 四分位数を用いた外れ値除去の方法 検者3 人で患側と健側それぞれにおいて𝑉𝑚𝑎𝑥を5 回の計測データで評価した結果を Fig7-3-1-2 に示す．

46

Fig7-3-1-2 𝑉𝑚𝑎𝑥を用いた健側・患側の有意差評価

２検者で患側・健側間の有意差が認められた．検者 C については，測定データのバラつ

47 7-3-2 せん断波の伝播方向のバラつき𝜃𝑠𝑑の評価 痙縮患者の患側は，せん断波が一様に伝播しなかったため，ROI 内のせん断波の伝播方向 の指標である，せん断波の角度の標準偏差𝜃𝑠𝑑を用いて評価を行った．𝜃𝑠𝑑の一例を Fig7-3-2-1 に示す． せん断波がx 軸方向に水平に伝播する状態を 0[deg]とする． Fig7-3-2-1 𝜃𝑠𝑑の例 検者3 人で患側と健側それぞれにおいて ROI 内の角度の標準偏差𝜃𝑠𝑑を5 回の計測デー タで評価した結果をFig7-3-2-2 に示す．

48

Fig7-3-2-2 𝜃𝑠𝑑を用いた健側・患側の有意差評価

全ての検者で健側・患側間で有意差が認められた．

健側の𝜃𝑠𝑑が検者間での差が大きいことについて考察をする．

49

Fig.7-3-2-3 に検者ごとの CFI と伝播図

検者A ではプローブと体表の間のゲル層に気泡があり，B モード画像が不明瞭である．ま た，検者 A,B,C を比較すると，筋肉の厚さが異なることから，検者間でプローブ配置，角 度の再現性に問題があったと考えられる．

50

第

8 章 結論

8-1 結論 本研究では，FDTD シミュレーションを用いてせん断波が筋繊維方向に伝播する条件を 評価し，被検者の状態が測定に与える影響について考察した．また，これを踏まえて痙縮患 者を測定し，痙縮側・非痙縮側で有意差が認められるか実験をした． 1. FDTD シミュレーションによるせん断波が筋繊維方向に伝播する条件評価  加振源―ROI 間距離が 30mm~50mm 程度離れていること  測定する組織は 1 層構造より 2 層構造の方が望ましい． 2. 痙縮腕の評価実験 痙縮腕は部分的に波面が歪み，筋組織内で硬いところと柔らかいところが複雑に分布 していることが推察された． 外れ値除去を行ったROI 内の最大値𝑉𝑚𝑎𝑥と，せん断波伝播方向のROI 内バラつき𝜃𝑠𝑑 をともに評価することで，健側・患側間の有意差が認められた． 8-2 今後の課題 1. このほかのせん断波が筋繊維方向に伝播条件を考察し，シミュレーション及び実機での 評価（3 次元の FDTD シミュレーションを用いることが望ましい）． 2. 痙縮腕をモデルとしたシミュレーション 3. より多くの痙縮患者で実験を行い，痙縮腕の評価方法，測定プロセスを明確化 4. がんの可視化への応用

51

謝辞

本研究を進めるにあたり，終始適切なご指導を頂いた群馬大学大学院理工学府 山越芳樹 教授に深く感謝申し上げます．また，本研究の臨床的有用性の評価は共同研究に基づいたも のであり，共同研究者である自治医科大学金谷裕司先生、倉品渉先生、紺野啓先生に深く感 謝いたします．日ごろの測定においてご支援いただいた遠坂俊明客員教授，荻野毅技官に感 謝申し上げます．研究を共にし，日々の実験や解析にご協力いただいた修士1 年 小久保大 輔氏に心より感謝いたします．最後に，研究室での学生生活においてお世話になりました山 越研究室の皆様に感謝の意を表します．

52

参考文献

[1] Naoyuki Takashima et al. Circulation Journal Vol. 81, no. 11, 1636-1646 (2017)

[2] Murry Aitken et al. Understanding Neuromuscular Disease Care. IQVIA Institute. Parsippany, NJ. (2018).

[3] 厚生労働省 平成 28 年「国民健康・栄養調査」

[4] 砂川和宏， 金井浩. "動脈壁組織性状診断を目的としたずり弾性波伝搬の計測とずり 粘弾性推定の検討. " 超音波医学 33. 1 (2006): 65-74.

[5] Oestreicher HL. "Field and impedance of an oscillating sphere in a viscoelastic medium with an application to biophysics. " The Journal of the Acoustical Society of America， 23. 6 (1951): 707-714.

[6] JING GAO et al.” QUANTITATIVE ULTRASOUND IMAGING TO ASSESS THE BICEPS BRACHII MUSCLE IN CHRONIC POST-STROKE SPASTICITY” Ultrasound in Med. & Biol., Vol. 44, No. 9, pp. 19311940, 2018