マイクロ波マンモグラフィのための高精度境界抽出法と逆散乱解析法の融合

マイクロ波マンモグラフィのための

高精度境界抽出法と逆散乱解析法の融合

学 籍 番 号

1731122

氏

名

則武 和輝

情報・ネットワーク工学専攻

電子情報システムコース

指 導 教 員

木寺 正平 准教授

副指導教員

野村 英之 准教授

提

出

日

平成

31

年

3

月

12

日

乳癌は，国内の女性の癌罹患率が最も高い癌であり，早期に治療を受けることで 10 年生存率が 90 ％程度となるため，定期的な検診による早期発見が重要である． しかしながら，既存の検診技術である X 線マンモグラフィは被曝や乳房圧迫等の 身体的負担が大きいため，国内の受診率は 40%程度に留まっている．これに対し， 近年マイクロ波マンモグラフィが注目されている．マイクロ波による乳房計測は， 低コストかつ被曝や痛みを伴わないため，検査頻度の高い簡易スクリーニング技術 として有望である．先行研究により，マイクロ波帯の乳癌細胞の複素誘電率は，脂 肪に対し約 10 倍，乳腺組織に対し約 1.2 倍のコントラストを有することが報告さ れており，この差異を識別する高精度な癌識別が期待されている．同コントラスト を利用するマイクロ波画像化法は，主にレーダ方式と逆散乱解析法に分けられる． 何れの手法においても，画像化精度は乳房境界抽出精度に大きく依存することが分 かっている． 本稿ではまず，レーダ方程式の一種である Envelope 法を用いた高精度境界抽出 法を導入する．一般に，マイクロ波マンモグラフィは近傍界計測であり，相互結合 によって乳房からの散乱波は，距離計測で利用される参照波形に対して歪むことが 分かっている．同波形の不整合による距離推定誤差が，乳房境界抽出精度に大きく 影響する．そこで，FDTD(Finite-difference time-domain) 法を用いて近傍界効果を 再現した参照波形を導入することで，高精度境界抽出を実現する． 次に，レーダ方式の一種である DAS(Delay-and-Sum) に着目し，高精度境界抽 出法の有用性を検証する．同方式では，表面反射波の抑圧が重要となっている．一 方、表面反射波を効果的に抑圧するためには，乳房境界を正確に求める必要がある． このため，前述した高精度境界抽出法によって乳房境界を推定し，FDTD 法によ り表面反射波を生成させる．これにより癌細胞の応答を抑圧せずに表面反射波を抑 圧することができるため，DAS 画像化精度を改善させることができる．表面反射 波抑圧性能をさらに改善するために，部分微分波形最適化による抑圧を提案する．

また，逆散乱解析法式の一種である DBIM(Distorted Born Iterative Method) に よる誘電率分布推定において，前述の境界抽出法を導入し，精度改善を図る．逆散 乱解析法は，画像化領域の推定精度に画像化精度が依存するため，高精度境界抽出 による精度改善が期待される．乳房 MRI 画像より抽出される精緻な乳房モデルを 用いた 2 次元及び 3 次元 FDTD 解析により，各手法の特性を評価する．

1 序論 1 1.1 研究背景 . . . 1 1.2 研究目的 . . . 2 2 マイクロ波マンモグラフィ 4 2.1 乳房組織の複素誘電率特性 . . . 6 2.2 レーダ方式 . . . 9 2.2.1 ビームフォーミング法 . . . 9 2.2.2 時間反転法 . . . 10 2.2.3 RPM 法 . . . 12 2.2.4 乳房境界抽出法 . . . 13 2.3 逆散乱解析法 . . . 16 2.3.1 積分方程式による問題の定式化 . . . 17 2.3.2 DBIM . . . 19 2.3.3 CSI 法 . . . 20 2.3.4 FBTS 法 . . . 22 2.3.5 各種逆散乱解析法の比較 . . . 24 2.4 各種表面反射波抑圧手法 . . . 25 3 FDTD 法と Envelope 法による高精度乳房境界推定法 28 3.1 システムモデル . . . 29 3.2 波形歪補正に基づく推定距離更新 . . . 30

4.1 システムモデル . . . 32 4.2 部分微分波形最適化による表面反射波抑圧 . . . 32 4.3 処理手順 . . . 33 4.4 2 次元 FDTD 解析による性能評価 . . . 34 4.4.1 計算モデル . . . 34 4.4.2 表面反射波抑圧結果 . . . 37 4.4.3 DAS による画像化性能の評価 . . . 39 5 高精度乳房境界推定法と逆散乱解析法の統合 41 5.1 システムモデル . . . 41 5.2 高精度境界抽出法と DBIM . . . 41 5.3 処理手順 . . . 42 5.4 2 次元 FDTD 解析による性能評価 . . . 44 5.4.1 計算モデル . . . 44 5.4.2 ROI 推定結果 . . . 46 5.4.3 誘電率再構成結果の検証 . . . 47 5.5 3 次元 FDTD 解析による性能評価 . . . 51 5.5.1 計算モデル . . . 51 5.5.2 ROI 推定誤差評価 . . . 54 5.5.3 誘電率再構成結果の検証 . . . 56 5.6 精度改善のための各種検討 . . . 58 5.6.1 行列方程式における係数調整 . . . 59 5.6.2 L2 ノルム正則化の導入 . . . 60 5.6.3 ROI 真値での DBIM の結果の確認 . . . 61 5.6.4 TV 拘束の導入 . . . 64 5.6.5 高勾配 ROI での DBIM による精度改善の検証 . . . 66 5.7 3 次元 DBIM による画像化の課題 . . . 68 6 結論 73

序論

1.1

研究背景

乳癌は，乳房組織である乳腺にできる悪性腫瘍であり，国内外の女性の癌罹患率 が最も高い癌である [1][2]．乳癌の 5 年生存率は治療開始病期が，I 期の場合 100%， II 期では 96%，III 期で 80.8%，IV 期で 37.1%となっている [3]．よって，早期発見， 特に大きさ 2cm 以下でリンパ節や他の臓器への転移のない I 期での早期発見が非 常に重要となっている．既存の検診技術として，X 線マンモグラフィや超音波によ る画像診断が行われている．X 線マンモグラフィは高い空間分解能を有するが，X 線被曝及び乳房圧迫による被験者への身体的負担が大きい問題を有する．超音波診 断は安全かつ非圧迫による検診が可能であるが，接触計測が必須であり，癌細胞の 識別はプローブを走査する検者の技能に依存する問題を有する．特に，日本人女性 の乳房は欧米女性に比べ比較的に小さく乳腺組織が密である．そのため，X 線マン モグラフィの身体的負担が大きく，日本人の受診率は 41%に留まっている [4]．こ れは，欧米女性の受診率が 70%を超えているのに対して低い水準である [4]．前述 のように，乳癌に関する社会課題として，早期発見が重要であるが受診率が低いと いう問題がある． こうした中，マイクロ波マンモグラフィが既存の検診技術と比較して，低コスト かつ被曝や非接触計測で痛みを伴わないため，より検査頻度の高い簡易スクリー ニング技術として注目されている．マイクロ波マンモグラフィに関する研究は 20 年程前から盛んであり，生体組織の複素誘電率を調査する研究，アンテナや観測シ

ステムなどのハードウェアに関する研究，画像化手法に関する研究に分けることが できる [5]．癌細胞の複素誘電率と，正常細胞である脂肪，乳腺組織の複素誘電率 の間には，それぞれ 10 倍，1.2 倍のコントラスト比を有することが報告されている [6]．このコントラストを基にした画像化手法は主に，レーダ方式と逆散乱解析法 (トモグラフィ方式) に分けられる．レーダ方式 [7] は，コントラストによって生じ る後方散乱波に対して結像処理を行い，癌細胞位置を特定する手法である．逆散乱 解析法 [8] は，コントラストによって生じる散乱波を用いて，散乱を引き起こした 複素誘電率分布を，積分方程式を基に逆問題を解いて直接再構成する手法である． レーダ方式においては，結像処理のための焦点距離算出と表面反射不要波の抑圧の ために乳房境界の情報を事前情報として利用する．また，逆散乱解析法では，基と なる積分方程式の積分領域かつ画像化領域 (ROI:Region of Intrested) を事前情報と して与える必要がある．このように，非接触計測を仮定する場合のマイクロ波マン モグラフィでの画像化精度は，いずれの手法においても乳房境界情報の推定精度に 依存する．

1.2

研究目的

本稿ではレーダ方式と逆散乱解析法の画像化精度の向上を目的とし，各々の画像 化手法に高精度境界抽出法を統合する手法を提案する．マイクロ波を用いた境界抽 出法は，何れの手法においてもアンテナと皮膚間の距離を推定し，同距離を基に乳 房境界を推定する．一方，マイクロ波マンモグラフィにおける計測は，一般的に， アンテナと乳房の距離が波長以下であるため，近傍界計測となる．近傍界計測では 波形が歪み，アンテナ皮膚間の距離の推定に誤差が生じるため，乳房境界推定精度 が劣化する問題があった．先行研究では，アンテナと皮膚間の相互結合による精度 劣化を FDTD 法で補正する手法を提案した． 本稿ではまず始めに，レーダ方式と高精度境界抽出を統合した手法の提案を行 う．レーダ方式では，表面反射波抑圧と焦点距離抽出の精度が，事前情報として与 える境界情報に依存する．従来の Envelope 法を用いた境界抽出では，近傍界計測 より ROI 推定誤差が大きい．そこで，高精度境界抽出法によって ROI を推定し， 事前情報としてレーダ方式に与える．また上記の手法に部分微分による表面反射波

抑圧を導入することで更なる画像化精度改善を図る． 次に，逆散乱解析法と高精度境界抽出法を統合した手法の提案を行う．逆散乱解 析法は，画像化領域を事前情報として与える必要があり，同領域の推定精度に画 像化精度が依存する．従来の Envelope 法を用いた境界抽出では，近傍界計測より ROI 推定誤差が大きい．そこで，高精度境界抽出法によって ROI を推定し，事前 情報として逆散乱解析法に与え，画像化精度の向上を図る． それぞれの画像化手法で，高精度境界抽出法を適用したことによって最終的な画 質が向上したことを確認するため，統計的な調査と MRI 画像より得られた精緻な 数値乳房モデルを用いた 2 次元及び 3 次元 FDTD 解析により性能評価を行う．

マイクロ波マンモグラフィ

本章では，始めに既存の乳癌検診技術について述べる．次に既存の検診技術に対 するマイクロ波マンモグラフィの利点について述べる．その後でマイクロ波帯での 乳癌組織の複素誘電率特性について述べ，同特性を用いた各種画像化について述 べる． X 線マンモグラフィは，最も一般的な検診技術である．X 線マンモグラフィは直 進する X 線を乳房に照射し，吸収率を画像化する．検出器は µm 間隔で数十万個並 べられ [9]，同検出器毎に得られた吸収率を処理することで高精細な空間分解能を 持つ画像が得られる．また，近年，照射点を走査し断層的な画像を得るトモシンセ シスと呼ばれる手法が製品化されている．実際の X 線マンモグラフィとトモシン セシスでの画像化結果を図 2.1 に示す．最新のトモシンセシス技術によってより高 精細かつ組織の重なりを確認することが可能となった．一方，X 線マンモグラフィ やトモシンセシスでは 1 回の照射で 1 枚の 2 次元画像を取得するため，低被曝化， 正常細胞と癌細胞の重なりの減少，体動抑圧，画質向上の観点から乳房を圧迫し計 測を行う．同圧迫の身体的負担が大きいことが X 線を利用する手法の問題点であ る．また，石灰化前の初期乳癌では画像化の基になる線吸収係数のコントラスト が非常に低くなる．線吸収係数は，乳腺組織で 0.80cm−1_{であるのに対し，乳癌は} 0.85cm−1で石灰化後は 12.5cm−1である [10]．石灰化後のコントラストは 15.6 倍程 度であるのに対し 1.06 倍程度のコントラストであるため，初期乳癌の検出は困難 であり誤警報率が高い問題を有する．

図 2.1: X 線 マ ン モ グ ラ フィ ，ト モ シ ン セ シ ス で の 画 像 化 例 (出 典:国 立 が ん 研 究 セ ン タ ー, 放射線技術部, 乳房 X 線検査 (マンモグラフィ)，URL： https://www.ncc.go.jp/jp/ncch/division/radiological technology /radiological diagnosis/xsenkensa/020.html) X 線マンモグラフィ以外の既存の検診技術として，超音波診断がある．超音波診 断では，超音波をプローブから乳房に放射し音響インピーダンスの違いを画像化す る．超音波診断は，被曝がなく安全かつ圧迫が不要な利点を有する．そのため，乳 腺が密な乳房や低年齢層へ適用されることが多い．一方で，同診断の有効性は確立 されておらず，厚生労働省指定の研究が行われている．また，超音波診断の識別は 検査者の技能に依存する問題があり，これに伴って乳癌診断のための技能を有する 超音波検査士不足が問題となっている [11]．

その他の検診技術としては，PET(Positron Emission Tomography) による画像 診断がある．PET-CT は，癌細胞がブドウ糖を多く消費する特異性を利用し，造影 剤を投入して CT によって画像化する手法である．造影剤が癌細胞に集まるため， 初期癌や癌の悪性度まで測ることができる．一方，検査コストが高いことや造影剤 を注射するため高血圧，糖尿病患者で識別率が下がる問題がある． こうした中，マイクロ波マンモグラフィが注目されている．X 線マンモグラフィ に対して安全かつ非圧迫である利点を持ち，超音波に対しては簡便という利点を持

つ．また，PET-CT に対しては安全かつ低コストである利点を持つ．マイクロ波マ ンモグラフィは 2000 年頃から盛んに研究されてきた分野であり次の 3 つに大別さ れる． • 乳房組織の複素誘電率の周波数特性の調査を行う研究 • 計測システムやアンテナなどのハードウェアに関する研究 • 観測された電気信号から画像化を行う研究 乳房組織の複素誘電率の周波数特性の調査によって，乳癌細胞と正常組織である脂 肪と乳腺との間には，それぞれ 10 倍，1.2 倍程度のコントラストを有することが報 告されている [6]．このコントラストによって電磁散乱界が生じ，同散乱界を計測 して画像化する研究が行われている．特に，画像化アルゴリズムに関する研究は 多岐にわたり，大別してレーダ方式と逆散乱解析法がある．レーダ方式は癌細胞境 界での後方散乱波を用いて結像処理を行い，同癌細胞位置を特定する手法である． 同方式の処理は比較的単純で計算コストは低いという利点を有する．一方で，画像 が反射波強度分布となることや，癌細胞が乳腺の中に埋まっている低コントラスト な状況で画像化精度が著しく劣化する欠点を有する．逆散乱解析法は散乱波を用い て，ヘルムホルツ型積分方程式を基にした逆問題を解くことによって複素誘電率分 布を直接再構成する手法である．逆散乱解析法は，乳房の複素誘電率分布を直接再 構成する利点を有する．一方で，逆問題を解く際に計算コストが高く悪条件下での 解析となることや，事前情報に画像化精度が依存する問題を有する．

2.1

乳房組織の複素誘電率特性

本節では [6] で調査されたマイクロ波帯での乳房組織の複素誘電率特性について 述べる．提案手法の性能評価で用いた数値乳房モデルは [6] を基に作成されたもので ある．University of Wisconsin(UW) と University of Calgary(UC) 大学病院で摘出 された乳房に対し，誘電分光法によって特性が調査された．図 2.2 に摘出された乳 房サンプルの一例を示す．同調査のサンプル数や患者データの概要を表 2.1 に示す．

図 2.2: 乳房組織サンプル，白組織：癌細胞，橙組織：脂肪 (出典:[6]“A large-scale study of the ultrawideband microwave dielectric properties of normal, benign and malignant breast tissues obtained from cancer surgeries”)

表 2.1: 乳房組織誘電特性調査 [6] の概要 総患者数 196 総サンプル数 319 患者年齢 35-87 歳 (UW) 19-90 歳 (UC) 摘出後から観測までの時間 20-221 min (UW) 11-240 min (UC) 観測中の組織温度 18.0 度-25.7 度 (UW) 19.9 度-27.2 度 (UC)

表 2.1 のサンプルに対し，VNA(vector network analyzer) を用いて 0.5GH ｚ-20GHz の特性を計測した．計測の際，サンプルを台に置きプローブを押し付けて S パラメータ (S11) を観測し特性を得る．図 2.3 に測定結果の 1 部を示す．ここで，正 常細胞サンプルの脂肪の割合が 0-30%であるグループが group1，31-84%が group2， 85-100%が group3 である．同図より，脂肪 (group1) と癌細胞では 10 倍程度のコン トラストがあり，乳房組織 (group3) と癌細胞では 1.2 倍程度のコントラストがある ことが確認できる．但し，同図の結果は 319 サンプルの各 group，癌細胞で得られ

図 2.3: 測定結果の Cole-Cole モデル (実線) によるフィッチング，(a)(b):group1-3(正 常細胞)，(c)(d)：癌細胞，(a)(c)：比誘電率，(b)(d)：導電率 (出典:[6]“A large-scale study of the ultrawideband microwave dielectric properties of normal, benign and malignant breast tissues obtained from cancer surgeries”)

た特性の中央値であり，個人差を有することに注意すべきである．また，in-vivo(生 体上) ではなく切除されたサンプルであり，温度が低く血が流れていないことにも 注意すべきである．[6] での調査結果をまとめると以下となる． •  乳房生体組織のマイクロ波帯での複素誘電率周波数特性：分散性モデルの 一種である Cole-Cole モデルによって再現される • 脂肪組織-癌細胞のコントラスト：10 倍程度 • 乳腺組織-癌細胞のコントラスト：1.2 倍程度 同調査結果を基に次節から述べる各種画像化が行われている．

2.2

レーダ方式

レーダ方式は，癌細胞と正常細胞のコントラストによって生じる後方散乱波を用 いて結像処理を行い，癌細胞の位置を検出する手法である．癌細胞と正常細胞のコ ントラストで生じる後方散乱波以外は不要波となり，同不要波の除去が大きな課題 となっている． 2.2.1 ビームフォーミング法 ビームフォーミング法の一種である MIST-BF 法 [7] は，生体組織で生じる周波 数分散を考慮して結像処理を行う手法である．始めに，周波数分散性を考慮しない DAS(Delay-and-Sum) 法について述べる．アンテナ素子 i での不要波を完全に除去 した癌細胞応答を si(t) とする．該当画像化ピクセル位置 r と各素子位置から算出 される伝搬経路から，遅延時間 τ (r)iを算出し結像処理を行う DAS 法は次の式で 表される． G(r) = I ∑ i=1 si(τ (r)i) (2.1) ここで，G(r) は結像結果，I は素子数であり，P (r) = G(r)2_{として電力分布 [dB]} を得ることができる．一方，MIST-BF 法では周波数空間で結像処理を行う．この 時に，周波数分散性を考慮した関数 Wi(r, ωj) を重みづけし，次式のように総和を とる． G(r) = I ∑ i=1 J ∑ j=1 Wi(r, ωj)Si(ωj)e−jωjτ (r)i (2.2) ここで，Si(ω) = F T [si(t)]，J は使用周波数領域である．また重みづけ関数 Wi(r, ωj) は各該当画像化ピクセル位置 r，アンテナ素子，周波数 ω から最適化により算出す る ([7]：3 章参照)．MIST-BF 法による画像化結果を図 2.4 に示す．癌細胞と正常細 胞のコントラストが大きく，アンテナ素子 i での不要波に埋もれていない癌細胞応 答 si(t) と遅延時間 τ (r)iが正確に得らた場合は図 2.4 のように癌細胞の検知が可能 である．しかし，表面反射波を完全に抑圧し正確な遅延時間 τ (r)iを求めることは， 特に，非接触計測の場合その両方が困難である．また，遅延時間 τ (r)i算出性能は 乳房境界抽出性能に依存していることが重要な点である．

図 2.4: MIST ビームフォーミングでの画像化例 (出典:Microwave imaging via space-time beamforming for early   detection of breast cancer)

2.2.2 時間反転法 時間反転法 (TR 法：Time Reversal)[12] は，癌細胞からの後方散乱波の応答時間 を基に画像化を行うレーダ方式の一種である．TR 法では，癌細胞からの応答時間， 式 2.12.2 での τ (r)iが得られていること仮定している点に注意する必要がある (前 述のように，本来は皮膚，乳腺からの不要波が大きく癌細胞応答を得ることが困難 である)．レーダ方式によって t = 0s から t = T s > 0 まで観測を行い，アンテナ素 子 i ごとに癌細胞応答 si(t) と，同応答時間 τ (r)iを得たとする．次に，時間逆伝搬 FDTD 法等で t = T s から t = 0s まで時間逆伝搬する．図 2.5 に [12] で用いられる 時間逆伝搬 FDTD を示す．時間逆伝搬中に，癌細胞応答時間 τ (r)iのタイミングで 得られた応答 si(t) をアンテナ素子位置に与える．これによって，各アンテナから 逆伝搬された応答が t = 0s で癌細胞位置に集中し，同位置の特定を行う．また背景 媒質は未知であるため，乳房の場合代表的な値で均質媒質を仮定し逆伝搬を行う． [12] では，均質媒質による逆伝搬誤差を抑圧するために，t = 0s で逆伝搬を止める のでなく，エントロピーが最小となる時間で止めることで癌細胞位置に鋭い応答を 得ることができる．[12] の TR 法による乳癌の画像化結果を図 2.6 に示す．ここで 同図の結果では，癌細胞と正常細胞のコントラストは [6] で報告されたものよりも 非常に大きいことを仮定している．また，乳房境界情報と皮膚の厚さ，誘電特性は 既知であることを仮定している．同図より，TR 法によって癌細胞位置が特定でき ていることが確認できる．一方，高コントラストである仮定が [6] で報告された乳 房の特性と異なっているため現状での画像化は困難であると考えられる．

図 2.5: 時間逆伝搬 FDTD の例，t = T 0 に時間反転 (左上→右上→左下→右下)(出 典 [12]:Time Reversal With the FDTD Method for Microwave Breast Cancer De-tection )

図 2.6: TR 法による画像化の例，上段：観測モデルの誘電率分布，下段：TR 法に よる画像化結果 (出典 [12]:Time Reversal With the FDTD Method for Microwave Breast Cancer Detection )

2.2.3 RPM 法

(a)波面 (b)赤点：候補点群，赤線：癌細胞境界 図 2.7: 背景媒質上の FDTD 法によって生成した波面と候補点群

マイクロ波マンモグラフィのための拡張 RPM(Range Point Migration) 法 [13] は，

背景媒質が既知であることを仮定し得られた癌細胞応答時間 τ (r)iを用いて FDTD 法によって波面を生成し，RPM 法 [14] の処理を適用する手法である．まず始めに， 得られた癌細胞応答時間 τ (r)i,jを用いて，既知であることを仮定した背景媒質上 で FDTD を実行し波面を生成する．ここで，i は素子を表し，j は癌細胞応答の第 j 到来波を表す．次に，同波面上に一定間隔で候補点群 qi,j(θe) を生成する．ここで， e = 1, .., E は各素子と到来波ごとに生成する候補点群数である．図 2.7 に生成した 波面と候補点群の例を示す．各素子ごとに図 2.7(b) の候補点群を生成し，RPM 法 による評価を行い推定点群を得る．但し，生成した波面の振幅は乳房の場合の不均 質が影響しているため，RPM 法 [14] での評価の内，振幅値による重みづけは行わ ない．図 2.8 にマイクロ波マンモグラフィのための拡張 RPM 法による画像化結果 を示す．同図より，癌細胞境界の推定が行えていることが確認できる．但し，背景 媒質が既知でという厳しい仮定である点に注意すべきである．

20 40 60 80 100120 x [mm] 20 40 60 80 100 120 140 y [mm] 5 10 15 20 ǫ 図 2.8: 赤点：マイクロ波マンモグラフィのための拡張 RPM 法による画像化結果 例，カラーバー：誘電率，中央癌細胞 (a) (b) 図 2.9: (a): 観測モデル，(b):レーザーによる推定結果，赤線：真値，黒線：推定 Voxel 位置 (出典：[15]，“Laser Surface Estimation for Microwave Breast Imaging Systems”)

2.2.4 乳房境界抽出法

本章では，マイクロ波レーダ方式による乳房境界抽出の意義と先行研究について 述べる．

࢘௖ ࢘_௧ ࢘_௥ ࢖௞ǡ௟ ࢖௞ǡ௟ାଵ ࢘_௞ ܴ෨ሺ࢘௧ǡ ࢘௥ሻ ࡿ Dielectric material 図 2.10: Envelope 法の原理 計測が可能であることがある．非接触計測において乳房境界情報は，事前情報とし て何れの画像化手法においても非常に重要である．レーダ方式においては，結像 処理の際の焦点距離 (遅延時間) 算出や表面反射波抑圧時の窓関数の精度が境界情 報に依存する．また，逆散乱解析法においては基になる積分方程式の積分領域 (乳 房存在位置) を事前情報として与える必要があり，同事前情報に逆問題の精度が依 存する．このような中，様々な境界抽出法が提案されている．レーザーセンサーを 用いた手法 [15] での乳房境界情報抽出結果を図 2.9 に示す．直進性のあるレーザー による推定によって，高精度乳房境界抽出が成されていることが確認される．しか し，マイクロ波マンモグラフィ装置にレーザーセンサー装置の機能を持たせるには ハードウェア的な複雑さが発生する．そのため，マイクロ波マンモグラフィにおい てはマイクロ波によって乳房境界抽出を行うことが望ましい．そのため，レーザー 等に対して直進性が低く距離分解能が低いマイクロ波での乳房境界抽出法が各種 提案されている [16][17]．マイクロ波による境界抽出では，何れの手法においても， 始めにアンテナ素子と乳房表面との距離を算出し，同距離を基に推定を行う． 2.2.4.1 Envelope 法 Envelope 法 [16] は，焦点を送受素子位置 ((rt, rr))，長軸半径を R(rt, rr) として 得られる楕円の包絡線を用いて目標境界を推定する手法である．ここで，推定した 距離を ˜R(rt, rr) とし，同距離は受信波形の整合フィルタ出力のピーク応答時間か

図 2.11: Envelope 法による乳房境界抽出, 破線：真値, 実線：推定値 (出典：[18],”Breast surface reconstruction algorithm for a multi-static radar-based breast imaging sys-tem”) ら算出される．図 2.10 に Envelope 法の概略を示す．まず，ある面領域 S(2 次元の 場合線領域) を K サンプル点に分け，同点を rk(k = 1, ..., K) とする．次に，任意 の点 rc_{から r} kまでの線分と，送受アンテナ (rt, rr) と距離 ˜R(rt, rr) から求まる楕 円との交点 pk,lを求める．複数ある交点 pk,lの中から，次式で推定点を決定する． pk = arg min pk,l √ |rc− p k,l| 2 (2.3) 領域 S 上のサンプル点 rkごとに上式で推定点を求めることで，包絡線上の推定点 群を得ることができる．ここで，rt, rrを，rt= rrとすると Mono-Sratic 型，隣り 合う素子で組合せると Bi-Static 型，複数組合せると Multi-Static 型となり，処理 の複雑さや状況に合わせて選択が可能である．Envelope 法は，表面形状が滑らか な場合，単純かつ低計算コストで高精度な境界情報を抽出することが可能である． 一方，同手法の精度は推定距離 ˜R(rt, rr) の精度に依存する．Envelope 法による乳 房形状推定 [18] の結果を図 2.11 に示す．同図から，距離推定が正しく行えている 場合，高精度に乳房境界情報を抽出できていることが分かる． 2.2.4.2 BSID 法

Envelope 法と似た手法である BSID(breast surface identification) 法 [17] は，素子

図 2.12: BSID による乳房境界抽出結果，グレースケールバー：誤差 (出典：[17], ” Estimating the breast surface using UWB microwave monostatic backscatter mea-surements”)

ライン (thin plate spline) によって曲面近似を行って境界情報を抽出する手法であ

る．まず始めに，素子位置 (rt, rr) と推定距離 ˜R(rt, rr) から IBST 法 [19] によって散 乱点 pm(r) を算出する．ここで，Envelope 法が推定距離群の積分処理であるのに対 し，IBST 法は微分処理となっている．次に，M 点ある散乱点 pm(r) を，N >> M 点に補間する．補間後のサンプル点 pn(r) に対し，乳房の領域を上下に分けて薄板 スプラインによって曲面近似を行い境界情報を得る．図 2.12 に [17] での乳房境界 抽出結果を示す．同図より高精度に乳房境界情報を抽出できていることが分かる． 但し，誤差が大きい部分では 2mm 程度の誤差を有することが確認できる．

2.3

逆散乱解析法

逆散乱解析法は，散乱波を用いて散乱を引き起こした複素誘電率分布を直接再構 成する手法である．電磁散乱界現象を表すヘルムホルツ型積分方程式を基に逆問題 を解くことで対象の電気的特性の直接再構成が可能である．同方式は散乱波を用い るため，後方散乱波のみを用いるレーダ方式に比べて活用できるデータが多い利点 を有する．一方，逆散乱解析法は複素誘電率分布を直接再構成するため未知数が多

く，画像化精度が初期設定の精度に依存することや計算コストが高いといった問題 を有する．また未知数が多いことは，悪条件下で逆問題になるため局所解に陥る問 題を有する．本章ではまず始めに，基となる積分方程式について述べ，その後で各 種逆散乱解析法について述べる． 2.3.1 積分方程式による問題の定式化 電磁界を表すマクスウェルの方程式から得られる，自由空間かつ電荷や電流が存 在しない場合の同次ヘルムホルツ方程式は次式となる． ( ∇2_{+ k}2 0 ) ET(r) = 0 (2.4) k0 = ω√ϵ0µ0 (2.5)

ここで，ET _{は全電界 (Total Field)，r は位置ベクトル，ϵ}

0, µ0は真空の誘電率，透 磁率である．式 2.5 において，電荷や電流は存在せず，比誘電率 ϵr(r) の散乱体が ある場合の非同次ヘルムホルツ方程式は次式となる． ( ∇2_{+ k}2 0 ) ET(r) =−k2₀ET(r)(ϵr(r)− 1) (2.6) 次に，式 2.6 の微分方程式の解をグリーン関数を用いた積分方程式にによって表 す．まず，一般的な微分方程式を微分演算子_{L を用いて表すと，} Lf(x) = −g(x) (2.7) となる．この時，同式の解がグリーン関数 G(x, x′) を用いて f (x) = f0(x) + ∫ ν G(x, x′)g(x′)dx′ (2.8) で求められる．ここでグリーン関数は， L[G(x, x′_{)] = −δ(x − x}′_{) (2.9)} を満たすものである．このようなグリーン関数を用いた式 2.8 を微分方程式_{Lf(x) =} −g(x) に代入すると，Lf0(x) = 0 より， ∫ ν L[G(x, x′_)]g(x′_)dx′ _{=−g(x) (2.10)} ∫ ν {L[G(x, x′_{)] + δ(x− x}′_{)} g(x}′_)dx′ _{= 0 (2.11)}

となり，同式が任意の x について成り立つため，式 2.8 によって解が得られる．こ のように，微分方程式_{Lf(x) = −g(x) を積分方程式 f(x) = −L}−1_{g(x) で解く利点} は近似解の構成に有利である点や，非同次項，非同次境界条件に依らない点である． 微分方程式の解をグリーン関数用いた積分方程式で解く式 2.8 を，微分方程式で ある非同次ヘルムホルツ方程式 2.6 に適用すると次の積分方程式が得られる． ES(r) = ET_(r)− EI_{(r) (2.12)} = ω2_µ∫ ΩG0(r, r′)E T_(r′_)o(r′_)dr′ _(2.13)

ここで，EI_{(r) は入射電界と呼ばれる自由空間での電場，E}S_{(r) は散乱電界，o(r) =}

ϵ(r)−ϵ0(r) は複素誘電率のコントラスト，G0(r, r′) は自由空間での r, r′間のグリー ン関数 (伝達関数) である．また，Ω は積分領域であり，同領域内に散乱体 ϵ(r) が 存在する．ヘルムホルツ型積分方程式の自由空間でのグリーン関数は式 2.9 の条件 や境界条件から解析的に求められる以下が一般的に用いられる． 3D : G0(r, r′) = 1 4π|r − r′|e ±ik|r−r′_| (2.14) 2D : G0(r, r′) =± i 4H (2) 0 (k|r − r′|) (2.15) ここで H(2) 0 は第 2 種ハンケル関数である． 式 2.13 での散乱電界とそれを表す積分方程式は，自由空間 ϵ0(r) と目標散乱体 ϵ(r) で得られる観測位置 r での電界の差分である．ここで，ある任意の背景媒質 ϵb(r) で得られた電界と目標散乱体 ϵ(r) との間の散乱電界は次式の積分方程式で表 すことができる． ∆ET(r) = ET_(r)− ET b (r) (2.16) = ω2µ∫_ΩGb(r, r′)ET(r′)∆o(r′)dr′ (2.17) ここで，Gb(r, r′) は背景媒質でのグリーン関数であり，∆o(r) = o(r)− ob(r) は背 景媒質とのコントラスト関数である．逆散乱解析法は，積分方程式である式 2.13， もしくは式 2.17 を基に逆問題を解く手法である．何れの手法においても次節では， 種々ある逆散乱解析法の中で代表的な 3 つの手法について，簡単な原理と利点，欠 点について述べる．また，何れの手法においても基となる積分方程式内の積分領域 Ω の決定が重要となる．

2.3.2 DBIM

DBIM(Distorted Born Iterative Method) は，順問題と Born 近似を行った線形 逆問題を交互に繰り返し解き推定解を求める手法である [8]．背景媒質の更新と共 に Born 近似とグリーン関数を更新することで，高コントラストな媒質であっても 妥当な解を得ることができる手法である．同手法は，数学的にはガウスニュートン 法と同じである．同手法の簡単な原理を以下で述べる．積分方程式 2.17 は，未知 数が ET_(r′_{) と ∆o(r}′_{) であり，非線形問題となっている．そこで，線形化のために} ET(r′) を背景媒質より解析的に求められる Eb(r′) にする Born 近似を行う．Born 近似を施した積分方程式は次式となる． ∆ET(r) = ET_{(r)− E}T b (r) (2.18) ≃ ω2_µ∫ ΩGb(r, r′)E b_(r′_)∆o(r′_)dr′ _(2.19)

上式を用いて ∆o(r′) について線形逆問題を解き，on+1(r′) = on(r′) + ∆o(r′) とす

ることで背景媒質を更新する．次に，更新された背景媒質上で順問題を解き，背景 媒質全電界 Eb_(r′_{) とグリーン関数 G} b(r, r′) を更新する．その後，再度 Born 近似 を施した式 2.19 を構成し線形逆問題を解いて更新量 ∆o(r′) を求める．上記の処理 をコスト関数 ∆ET(r) ₂が小さくなるように ∆o(r′) を繰り返し更新し，終了条件 を満たした最終的な複素誘電率分布を得る．DBIM では，背景媒質を更新によって Born 近似 Et_(r′_{)≃ E}b_(r′_{) の誤差が低減するため，最終的な推定精度が改善するた} め，高いコントラストを有する不均質な散乱体であっても高精度な推定が可能であ る．一方，同手法の欠点は，順問題と逆問題を交互に解く必要があり，計算コスト が高い点である．また，順問題と逆問題を解く手法は種々あるため，状況にあった 解法を選択する難しさもある．図 2.13 に，[8] での DBIM による画像化結果を示す． 乳房のような高コントラストかつ分散媒質でも，初期設定や各種パラメータを適切 に設定することである程度の画像が得られている．但し，[8] では乳房境界と皮膚 の厚さは真値であることを仮定している．

図 2.13: DBIM 画像化例,(a)-(c):Class3 乳房モデル，(d)-(f):DBIM による推定結果， (a)(d):x-y 平面，(b)(e):y-z 平面，(c)(f):x-z 平面 (出典:[8]“Microwave imaging via space-time beamforming for early detection of breast cancer”)

2.3.3 CSI 法

CSI(Contrast Source Inversion) 法 [20] は，式 2.13 において未知数となっている

全電界 (Source)ET_(r′_{), o(r}′_{) とコントラストを一つの未知数，コントラストソース} として求める手法である．CSI 法では，物理方程式とデータ方程式と呼ばれる 2 つ の行列形式積分方程式を用いる．同 2 式は式 2.13 から導かれ，以下のように表さ れる． 物体方程式 (領域 D 全体の電界のふるまいに関する式) uk = uInck + G D kχkuk = uInck + G D kwk (2.20) データ方程式 (取得データに関する式) fk= uk− uInck = G S kχkuk = GSkwk (2.21)

ここで，uk, uInck は領域 D 全体の全電界，入射電界，S はアンテナが位置する面領 域 (2 次元：線分領域)，fkは S 上で得られる観測 (散乱) 電界，χkは複素誘電率分 布と自由空間のコントラスト，k = 1, ..., K は未知数であり，コントラストソース は次式で表される． wk = χkuk (2.22) 積分方程式では未知数が χi, uiの 2 つであったのに対し，上式では未知数が wkの 1 つになり線形化されている．具体的な解き方は，まず始めに適当な初期設定とし て χk,n=0, uk,n=0を与える．ここで，n は更新回数である．次に，以下のコスト関数 を最小化にするようなコントラストソース wkを逆問題を解き求める． F = ∑ k fk− GSkwk 2 S ∑ k∥fk∥ 2 S + ∑ k χkuInck − wk+ χkGDkwk 2 D ∑ k∥χku Inc k ∥ 2 D (2.23) wk = arg min wk F (2.24) 次に，上記で求まったコントラストソース wkを用いて式 2.20 で全電界 uk,nをアッ プデートし，さらに式 2.22 でコントラスト χk,nを更新する．更新されたコントラ ストを用いて再度式 2.24 でコントラストソースを求める．上記の繰り返しを行う ことで最終的なコントラスト χk,nを得る． CSI 法の利点は上記の処理手順からわかるように順問題を解く必要がないこと である．一方，繰り返し回数は一般に多く計算コストは高い問題点を有する．ま た，コントラストソース wkは数学的，物理的に自明解を持たないことが分かって いる．しかし，多くの問題において有効性が確認されており，様々な研究がなされ ている．その中に CSI 法に全変量によるエッジ保存のための拘束を加えた MR-CSI 法 (Multiplicative Regularized Contrast Source Inversion)[21] があり，CSI 法より も空間分解能が高い結果を得ている手法もある．MR-CSI 法と DBIM による脚の 画像化結果比較を図 2.14 に示す．同図より，正則化を導入した MR-CSI 法によっ て高精度に脚部の画像化が行えていることが確認できる．一方で，DBIM では複数 の周波数を組合せて画像再構成を行うことができるが，CSI 法では各周波数毎にコ ントラストを求める必要があり，分散性媒質に対する画像化では比較検討が必要で ある．

図 2.14: 画像化例,(a)(b):脚部複素誘電率，(c)(d):MR-CSI 法による画像化結果， (e)(f):DBIM による画像化結果，(a)(c)(e):複素誘電率実部，(b)(d)(f):複素誘電率虚 部 (出典:[22]”Comparison of an Enhanced Distorted Born Iterative Method and the Multiplicative-Regularized Contrast Source Inversion method”)

2.3.4 FBTS 法 FBTS(Forward-Backward Time-Stepping) 法 [23] は時間空間において，観測電界 と，推定複素誘電率分布で得られる推定電界との差分を最小にするように，複素誘 電率分布を最適化する手法である．推定に用いる評価汎関数は次式となる． Q(p) = ∫ T 0 M ∑ m=1 N ∑ n=1 ∥vm(p; rn, t)− ˜vm(ptrue, rn, t)∥ 2 (2.25) ここで，p は複素誘電率のパラメータ (e.g. 単極デバイモデル)，˜vm(ptrue, rn, t) は， 送信 m の時の rnの位置での全電界の真値，vm(p; rn, t) は p を仮定した時の全電界， M, N, T はそれぞれ送信点数，未知数，計測時間である．上記の評価汎関数のフレ シェ微分による勾配法をより解を得る．簡単な処理手順を以下に述べる．まず始め に，適当な分布 p を仮定し t = 0 から Forward time-stepping によって vm(p; rn, t)

を計算する．計算された全電界より，残差等価波源は次式となる． um(r, t) = N ∑ n=1 [vm(p; rn, t)− ˜vm(ptrue, rn, t)] δ(r− rn) (2.26) 同残差等価波源を用いて t = T から Backward time-stepping を行い，次式の随伴 方程式より随伴解 wmを得る． Lwm= um (2.27) ここで，L は随伴微分演算子である．同随伴解を用いて汎関数の勾配を次式で得る． g(r) = 2 M ∑ m=1 ∫ T 0 wm(p; rn, t) ∂ ∂(ct)vm(p; rn, t) (2.28) 同勾配から勾配法によってパラメータ p を評価汎関数 2.25 が小さくなるように更 新する．更新された p を用いて再度 Forward time-stepping によって vm(p; rn, t) を 計算し，上記の処理を収束条件を満たすまで繰り返す． FBTS 法は，時空間での評価関数を最小化するため，実験装置とシミュレーショ ンとの間の周波数特性誤差に対して堅牢である利点を有する．また，時間領域デー タを基にするためデータが多く，非線形性が減少している利点を有する．一方，1 回の更新ごとに time-stepping，つまり FDTD による解析が必須であるため，高計 算コストである欠点を有する．図 2.15(e) に，FBTS 法 [23] による画像化結果を示 す．癌細胞と正常細胞のコントラストが図 2.15(a) のように高い場合，癌細胞の画 像化が行えている． 図 2.15: FBTS 画像化例,(a):Class3 乳房モデル，(e):FBTS による推定結果 (出 典:[23]“Advances in the 3-D Forward-Backward Time-Stepping (FBTS) Inverse Scattering Technique for Breast Cancer Detection”)

2.3.5 各種逆散乱解析法の比較 前述の逆散乱解析法のまとめと比較を表 2.2 に示す．なお，FBTS 法の利点につ いては総合大会 2018 で長崎大学の田中先生に伺った内容も含まれている点に注意 すべきである． 表 2.2: ROI 推定誤差 Err_Ωˆ 比較 DBIM CSI 法 FBTS 法 原理 順問題解析と Born 近 似による線形逆問題を 交互に繰り返し解く 物体方程式とデータ方 程式の逆問題を交互に 解く 順問題解析 (FT) と随 伴界解析 (BT) によっ て勾配を計算し推定値 を更新 利点 高コントラストかつ不 均質目標であっても妥 当な解を得る 順問題解析が不要であ る 時間領域のコスト関数 のため，データが多く 非線形性が低い 順問題と逆問題の解析 手法が数多あり，状況 に合わせ解析手法を選 択できる 計算コストが他の 2 手 法よりも少ない 時間領域のコスト関数 のため，実験データに 対する堅牢性が高い 欠点 高計算コスト 単一周波数ごとの推定 であるため分散媒質に 対し精度劣化が考えら れる 高計算コスト (DBIM より高い) 各逆問題で求める解が 更新量であり正則化な どを行いにくい 評価関数の自由度が低 い アルゴリズムが複雑

2.4

各種表面反射波抑圧手法

一般に，空気と皮膚，皮膚と脂肪で生じる表面反射波は非常に大きく，同不要波 に癌細胞応答が埋もれてしまう．前述の通り，レーダ方式において必要なのは癌細 胞応答のみであり，表面反射波は抑圧の対象となる．そのため，レーダ方式での画 像化は前処理として表面反射波抑圧を行う．レーダ方式による画像化精度は，同表 面反射波抑圧性能に依存している．こうした中，マイクロ波マンモグラフィにおけ る表面反射波抑圧手法が各種提案されている．代表的な平均波抑圧，Wiener フィ ルタ，エントロピー法についてレビュー論文 [24] を参考に以下にまとめる．また以 下では，bi(t) を観測信号，si(t) を抑圧後の信号，i を素子番号，I を素子数とする． 2.4.0.1 平均波形抑圧 平均波形抑圧は，表面で発生する波形の形状が各素子で相似していること仮定し 以下の式で抑圧を行う． si(t) = bi(t)− 1 I− 1 I ∑ j=1,i̸=j bj(t) (2.29) 平均波形抑圧では，該当信号素子以外の癌細胞応答が含まれる信号を利用するた め，癌細胞応答まで抑圧してしまうリスクを有する． 2.4.0.2 Wiener フィルタ Wiener フィルタは，信号処理で広く適用されるフィルタ処理であり，最適化に よりフィルタの重みづけ関数を決定する．マイクロ波マンモグラフィでは，該当素 子の受信信号を，該当素子以外の受信信号を用いて最小にするような重みづけ関 数 qiを，表面反射波があるであろう応答時間 ts ≤ t ≤ teの間で最小化するように 求める．この最適化問題を解く手法を種々あるが，代表的なものは最初二乗法であ る．抑圧信号は次式で表される． si(t) = bi(t)− qiTDi(t) (2.30)

ここで，Di(t) は Di(t) = [ dT₁, ..., di− 1T, di+ 1T, ..., dTI ] (2.31) dk = [bk(t′ − M), ..., bk(t′), ...bk(t′+ M ), ] (2.32) であり，M は窓関数の幅である．該当素子以外での受信信号を，t′_{を中心とした幅} を持った波形情報で構成された行列 Di(t) を，巧く組み合わせて波形の最適化を行 うための重みづけ関数 qiを次式で求める． qi = arg min qi te ∑ t=ts ∥bi(t)− qiTDi(t)∥2 (2.33) 同式で得られた重みづけ関数 qi を用いて式 2.30 によって表面反射波抑圧を行う． Wiener フィルタによる表面反射波抑圧は性能が非常に良いが，パラメータが複数 あることや，平均波形抑圧と同様に癌細胞応答を含む波形を用いるため，癌細胞応 答まで抑圧してしまうリスクを有する． 2.4.0.3 エントロピー法 エントロピー法は，平均波形抑圧と同様に各素子での波形形状が相似しているこ とを仮定し，エントロピーが高い，つまり波形形状が一致している時間を窓関数に よって完全に 0 にする手法である．反対に，癌細胞応答は各素子ごとに異なるため エントロピーが低く，窓関数によって抑圧されずそのまま通過する．窓関数の具体 的な算出方法を以下に示す．まず，各受信信号を確率変数 (p_{≥ 0,}∑p = 1) として 扱うために以下のように正規化する． pi(t) = ∥b i(t)∥2 ∑I j=1∥bj(t)∥2 (2.34) 次に，各時間 t での情報量としてのエントロピーを次式で求める． Hα(t) = 1 1− αlog { _I ∑ i=1 [pi(t)]α } (2.35) 同エントロピーは，ノイズ等の影響を受けやすいため，次の平滑化エントロピーを 求める． H_αs(t) = 1 M k=t+M∑ k=t Hα(t) (2.36)

(a)観測モデル(非接触) (b)平均波形抑圧 (c) Wienerフィルタ 図 2.16: [25] での表面反射波抑圧手法の比較 ここで，M は平滑化幅である．平滑化エントロピーをを用いて以下のように窓関 数を構成する． W (t) =      0, eHαs(t) _{> N} 0 1, otherwise (2.37) ここで，1 < N0 < I は経験的に決定するパラメータである．この窓関数を用いて 以下の式で抑圧信号を得る． si(t) = W (t)bi(t) (2.38) エントロピー法では，窓関数によって完全に表面反射波を除去できる一方，癌細胞 位置が皮膚に近い場合は癌細胞応答も除去してしまうリスクを有する． 図 2.16 に [25] での表面反射波抑圧と Wiener フィルタによる抑圧の例を示す．同 図 (a) に示すように非接触計測であることを想定している．また，同図 (b)(c) で示 す結果は [25] の Scenario III の皮膚の厚さが未知であるを想定した上での結果であ る．同図より平均波形抑圧と Wiener フィルタによって表面反射波抑圧が行えてい ることが確認できる．但し，平均波形抑圧では癌，皮膚に近い領域で抑圧性能が劣 化している．Wiener フィルタでは，乳房モデル中央付近に虚像が生まれている． 先行研究では，上記の何れの手法において問題となっている癌細胞応答抑圧のリ スクを避けた再現波形による表面反射波抑圧法を提案した．本稿 4 章では，先行研 究での提案手法同様に癌細胞応答抑圧を避けたまま性能をより向上させるために， 部分微分による表面反射波抑圧を提案する．

FDTD

法と

Envelope

法による高精度乳房境

界推定法

2 章で述べたように，何れの画像化手法においても乳房境界情報は非常に重要であ る．同時に，マイクロ波マンモグラフィでは装置等の観点からマイクロ波による境 界抽出が望まれる．こうした中，マイクロ波による乳房境界抽出法が各種提案さ れており，何れの手法もアンテナ素子と乳房表面の距離を基にしている．一方で， マイクロ波マンモグラフィの観測は，一般にアンテナ素子を乳房境界から近接した 場所に配置するため波長以下の近傍界計測となる．そのため，近傍界と相互結合の 影響によって観測波形が歪み，放射電界のみの参照波形との間に波形の不整合が生 じる．図 3.1 に，参照波形と近傍界計測の影響が大きく歪んだ波形と，同影響が小 さく歪の少ない波形を示す．同図からわかるように，観測距離が半波長以下の観測 波形では参照波形と異なる波形となることが確認できる．同波形不一致によって， アンテナ素子と乳房間の距離推定に誤差が生じ，これに伴い乳房境界抽出精度が劣 化する．[26] では，既知の形状の水のカップの反射波を参照波形としている．しか し，アンテナ素子と乳房間の距離は素子によって異なり，近傍界計測の影響も素子 によって異なる．本章では，近傍界計測の影響を FDTD 法によって再現し，再現 波形との時間シフトを抽出することによって推定距離を更新する高精度境界抽法に ついて述べる．

3.1

システムモデル

図 3.2 にシステムモデルを示す．乳房のモデルは，損失，等方，分散性媒質を仮 定する．アンテナアレーを乳房近傍に配置し，送信点を切り替えと全アンテナでの 受信を繰り返し，全ての送受アンテナの組合せでデータを取得する．送信位置 rt， 受信位置 rt，時間 t での電界変化を Eobs(rt, rr; t) とする．ある送受信素子の組み合 わせ (rt, rr) で得られる受信波形と，参照波形 Eref(t) との相互相関関数のピーク値 から τ (rt, rr) を抽出し，距離 R(rt, rr) = cτ (rt, rr)/2 を得る．但し c は真空中の光 速である． 0 0.5 1 1.5 2 2.5 3 Time [nsec.] -2 -1.5 -1 -0.5 0 0.5 1 Norma Amplitude 図 3.1: マイクロ波マンモグラフィにおける受信波形と参照波形の形状不一致の例， λ:波長 7[ 5[ ƌĞĂƐƚŵĞĚŝĂ ;ĂĚŝƉŽƐĞ͕ĨŝďƌŽŐůĂŶĚƵůĂƌͿ 6NLQ 5[ 5[ 5[ 5[ 5[ 5[ 5[ 図 3.2: システムモデル

ཧ↷Ἴᙧ
_{௧ǡ ௥} ൌ ሺ௧ǡ ௥ሻȀʹ dĂƌŐĞƚ ほ  ほ Ἴᙧ Ἴᙧ䛾୙ᩚྜ ㊥㞳᥎ᐃ ՜̃_{ሺ௧ǡ ௥ሻ}
 _{௧ǡ ௥} ൌ ̃_{ሺ௧ǡ ௥ሻȀʹ} ƐƚŝŵͲ ĂƚĞĚ ȟ &d ;^ŝŵƵůĂƚŝŽŶͿ ෌⌧Ἴᙧ Ἴᙧ⿵ṇ ほ Ἴᙧ ෌⌧Ἴᙧ ο࣎෤ ᫬㛫䝅䝣䝖 㔞οᢳฟ
_{௧ǡ ௥} hƉĚĂƚĞĚ ൌ
 ௧ǡ ௥ ൅ οͬϮ ƚƌƵĞ ㊥㞳᭦᪂ 図 3.3: 波形歪補正に基づく距離点更新の概略図

3.2

波形歪補正に基づく推定距離更新

図 3.3 に波形歪補正に基づく距離点更新の概略図を示す．素子-乳房間距離 R(rt, rr) = cτ (rt, rr)/2 の状況で観測を行い，受信信号 Eobs(rt, rr; t) を得る．受信信号 Eobs(rt, rr; t) と参照波形 Eref_{(t) を用いて ˜τ (r} t, rr) を算出し距離 ˜R(rt, rr) を推定する．この時， 近傍界計測での相互結合作用より波形が歪むため，抽出した時間シフト量は次式の ように誤差を含む． ˜ τ (rt, rr) = τ (rt, rr) + ∆τ (rt, rr) (3.1) ここで，∆τ (rt, rr) は時間シフト以外の波形不一致によって生じる推定誤差である． 次に，推定距離群 ˜R(rt, rr)|t = 1, .., M, r = 1, .., N を用いて Envelope 法を適用し 誤差 ∆τ を含んだ境界情報 ˜Ω を推定する． ˜Ω と代表的な乳房組織のパラメータを 用いて再現乳房数値モデルを生成し，同モデルを用いて FDTD 法により再現波形 Eest_(r t, rr; t) を得る．再現波形 Eest(rt, rr; t) は，相互結合作用を含んでおり，放射 電界のみの参照波形よりも，観測波形に形状が近い．同再現波形を用いて次式で時

ཧ↷Ἴᙧܧ୰ୣ୤ 䛸ほ Ἴᙧܧ୭ୠୱ 䜢⏝䛔䛯㊥㞳ܴ෨䛾᥎ᐃ ㊥㞳ܴ෨䜢⏝䛔䛯ŶǀĞůŽƉĞἲ䛻䜘䜛ቃ⏺ȳ෩᥎ᐃ ቃ⏺ȳ෩䜢⏝䛔䛯෌⌧ஙᡣᩘ್䝰䝕䝹䛾⏕ᡂ ෌⌧䝰䝕䝹䜢⏝䛔䛯&dἲ䛻䜘䜛෌⌧Ἴᙧܧୣୱ୲_䛾⏕ᡂ ᭱⤊ⓗ䛺ܴ෠䜢⏝䛔䛯ŶǀĞůŽƉĞἲ䛻䜘䜛ቃ⏺ȳ෡⟬ฟ ほ Ἴᙧܧ୭ୠୱ 䛸෌⌧Ἴᙧܧୣୱ୲ 䜢⏝䛔䛯㊥㞳䛾᭦᪂ܴ෠ zĞƐ EŽ ཰᮰ุᐃ 図 3.4: FDTD による波形補正に基づく高精度境界抽出法のフローチャート 間シフト量 ∆˜τ (rt, rr) を求める． ∆˜τ (rt, rr) = arg min τ [ Eobs(rt, rr; t) ⋆ Eest(rt, rr; t) ] (τ ) (3.2) ここで，⋆ は相互相関演算子である．同式で算出された時間シフト量 ∆˜τ (rt, rr) は， 観測波形と再現波形の形状が近いため，乳房までの距離推定誤差によって生じた値 に近い (∆˜τ (rt, rr) ≃ ∆τ(rt, rr))．次に，算出された時間シフト量 ∆˜τ (rt, rr) を用 いて推定距離を更新する． ˆ R(rt, rr) = ˜R(rt, rr) + c∆˜τ (rt, rr)/2 (3.3) 上記の処理を ˆR が収束するまで行い，最終的に得られた推定距離に Envelope 法を 適用し高精度な境界 ˆΩ を得る． 上記の高精度境界抽出法のフローチャートを図 3.4 に示す．

FDTD

法による表面波抑圧法とレーダ画像化

法の統合

先行研究において，2.4 で紹介した何れの手法において問題となっている癌細胞応 答抑圧のリスクを避けた再現波形による表面反射波抑圧法を提案した．本稿では， 先行研究の性能を向上させるために，部分微分波形最適化による表面反射波抑圧を 提案する．

4.1

システムモデル

システムモデルは 3 章と同様の図 3.2 である．但し，モノスタティック (rt = rr) で観測された電界 Eobs_(r i; t) を用いて画像化を行う．ここで，i = 1, ..., M は送受素 子数とする．これに伴い，Envelope 法も Mono-Sratic 型 (rt= rr) を仮定する．

4.2

部分微分波形最適化による表面反射波抑圧

先行研究で提案した手法では，癌細胞応答の抑圧を避けるため高精度に推定され た境界抽出情報を用いて，表面反射波再現波 Eest_(r i; t) を FDTD で計算し，次式で 抑圧波形を得た． ˜ sconv(ri; t) = Eobs(ri; t)− maxtEobs(ri; t) maxtEest(ri; t) Eest(ri; t) (4.1)

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Time [nsec] -1 -0.5 0 0.5 1 Normalized Amplitude α_=-0.4 α=-0.2 α=0.0 α=0.2 α=0.4 図 4.1: 部分微分係数による波形変化の一例 この時，Eest_(r i; t) を算出する際の乳房再現モデルと，実際の乳房との複素誘電率 分布の差異より，抑圧後の波形 ˜sconv_(r i; t) には表面反射波残差が存在する．そこで， 本稿では次式の部分微分波形最適化による表面反射波抑圧を行う． ˜

sprop(ri; t) = Eobs(ri; t)− ˜AEest(ri; t− ˜τ; ˜α) (4.2) ここで，最適化パラメータは次式より求める． ( ˜ A, ˜τ , ˜α ) = arg min (A,τ,α) ∫ Tr+Tw Tr ∥Eobs_(r i; t)− AEest(ri; t− τ; α)∥2dt (4.3) ここで，Tr ≤ t ≤ Tr+ Twは表面反射波応答時間，−1 ≤ α ≤ 1 は部分微分係数で 部分微分は次式で定義される． E(ri; t; α) =F−1[(jω) α_{F [E(r} i; t)]] (4.4) 図 4.1 に部分微分係数 α を変更した際の波形の変化を示す．同図より，部分微分に よって波形形状が変化していることが確認できる．

4.3

処理手順

本稿での提案手法による画像化の処理手順を以下に示す．ここで，本稿で新たに 導入したのは手順 3 である． 手順

1. 高精度境界抽出法 (3 章) を適用し境界情報 ˆΩ を抽出 2. ˆΩ と代表的な乳房のパラメータを用いて FDTD 法により再現波形 Eest_(r i; t) を算出 3. 部分微分波形最適化による表面反射波抑圧→ ˜sprop(ri; t) 4. ˜sprop_(r i; t) を用いた DAS 法による画像化→ P (r) 4.3.0.1 比較手法 提案手法による表面反射波抑圧の性能評価比較を行うため，以下の 2 つの従来手 法と提案手法の比較を行った． 従来手法 1:境界情報：波形補正なし，式 4.1 による抑圧 従来手法 2:境界情報：波形補正なし，MS-MIST[27] での調整平均波形抑圧   非接触計測時に対応するため以下の式による抑圧を行う． ˜ sconv(ri; t) = Eobs(ri; t)− 1 M − 1 M ∑ j=1,j̸=i Eobs(rj; t− ˜τj) (4.5) ここで，˜τjは各ピーク応答時を揃える遅延時間であり，観測波形同士の相互 相関によって求める． 提案手法:境界情報：波形補正あり，式 4.2 による抑圧

4.4

2

次元

FDTD

解析による性能評価

4.4.1 計算モデル 中心周波数 2.45GHz，帯域幅 2.7GHz レイズコサイン変調パルスを送信波形とす る．2 次元問題，TE 波 (Transverse Electric) を仮定する．Mono-Sratic 型送受信素

子数は 31 とする．単極デバイモデル (ϵ(ω) = ϵ∞+ ∆ϵ 1 + jωτP + σs jωϵ0 , τP : 緩和時間) を仮定した分散性 FDTD 法を用いて散乱データを取得する．FDTD でのグリッド サイズは 0.5mm，時間ステップは 1.18ps である．

20 40 60 80 100 120 x[mm] 20 40 60 80 100 120 y[mm] 5 10 15 20 Relative permittivity( ǫ ( ∞ )) (a) 20 40 60 80 100 120 x[mm] 20 40 60 80 100 120 y[mm] 0 5 10 15 20 25 30 35 Relative permittivity( ∆ ǫ = ǫ0 -ǫ∞ ) (b) 20 40 60 80 100 120 x[mm] 20 40 60 80 100 120 y[mm] 0 0.2 0.4 0.6 0.8 Conductivity[S/m]( σ ) (c) 図 4.2: CASE1 乳房均質モデル (a):ϵ_∞，(b):∆ϵ，σs，白点：素子位置 20 40 60 80 100 120 x[mm] 20 40 60 80 100 120 y[mm] 5 10 15 20 Relative permittivity( ǫ ( ∞ )) (a) 20 40 60 80 100 120 x[mm] 20 40 60 80 100 120 y[mm] 0 5 10 15 20 25 30 35 Relative permittivity( ∆ ǫ = ǫ0 -ǫ∞ ) (b) 20 40 60 80 100 120 x[mm] 20 40 60 80 100 120 y[mm] 0 0.2 0.4 0.6 0.8 1 Conductivity[S/m]( σ ) (c) 図 4.3: CASE2 乳房不均質モデル (a):ϵ_∞，(b):∆ϵ，σs，白点：素子位置 図 4.2 と図 4.3 に FDTD で数値計算を行う CASE1 と CASE2 の乳房モデルと白点 で素子位置を示す．図 4.3 の CASE2 モデルは MRI 画像から統計的に複素誘電率分 布を決定したものであり，脂肪が支配的なモデル (Class2) である [28]．また，デバ イモデルにおける緩和時間 τPは 1.5×10−11[s] とする．CASE1 のモデルは，CASE2 のモデルの皮膚を除く乳房で平均をとり均質媒質にしたものである．合わせて両 CASE の (x,y)=(41mm,48mm) の位置に 4mm の癌細胞を付加する．同癌細胞のデ バイパラメータは，(ϵ_∞, ∆ϵ, σs) = (20.0, 38.0, 0.8) とする．皮膚近くに癌細胞を配 置し，従来手法 2 における癌細胞応答抑圧のリスクに対する提案手法の有効性を確 認する． 観測散乱データにホワイトガウシアンノイズを付加する．同ノイズの SNR は所 望波である癌細胞応答をシグナルとした場合，16dB である．この時，波源での直 接波をシグナルとした場合の SNR は 70dB であり，マイクロ波マンモグラフィに

0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -5 0 5 Amplitude ×10-5 (a)観測波形 0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -5 0 5 Amplitude ×10-5 (b)抑圧後の波形：従来手法1 0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -5 0 5 Amplitude ×10-5 (c) 抑圧後の波形：従来手法2 0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -5 0 5 Amplitude ×10-5 (d)抑圧後の波形：提案手法 図 4.4: CASE1：波形データ比較 おいて妥当なノイズである [29]． 次に，提案手法における各種パラメータの設定，処理の内容を以下に述べる． 提案手法式 4.3 における Trの値は各素子ごとに Tr,i= 2 ˆR(ri)/c とし，Twは送信波 形パルス幅の 1.5 倍である Tw = 0.75ns とする．式 4.3 の最適化では非線形最適化 法の Nelder-Mead 法を用いる．Nelder-Mead 法では，初期値 (A0_{, τ}0_{, α}0_{) への依存} があるため 2 度最適化を行う．1 度目の最適化では，(A0_{, τ}0_{, α}0_{) = (1, 0, 0) とし最} 適化を行い各素子 i でパラメータ(_A˜_{i, ˜τi, ˜αi})_{を得る．同パラメータの中には初期} 値依存から，局所解に陥った解が存在するため，全素子での平均値(_{A, ¯}¯˜ _{τ, ¯˜ α˜})を新 たな初期値に設定し再度 NelderMead 法を適用する．

0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -6 -4 -2 0 2 4 6 Amplitude ×10-5 (a)観測波形 0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -6 -4 -2 0 2 4 6 Amplitude ×10-5 (b)抑圧後の波形：従来手法1 0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -6 -4 -2 0 2 4 6 Amplitude ×10-5 (c) 抑圧後の波形：従来手法2 0 0.5 1 1.5 2 Time [nsec] 0 1 2 3 4 5 6 θ [rad] -6 -4 -2 0 2 4 6 Amplitude ×10-5 (d)抑圧後の波形：提案手法 図 4.5: CASE2：波形データ比較 4.4.2 表面反射波抑圧結果 図 4.4 に CASE1 での観測波形と各手法での抑圧後の波形を示す．同図の横軸は 時間，縦軸は乳房中心からの角度 θ にあるアンテナで得れた波形を示す．同図 (a) の 0.5ns 付近に見える応答が表面反射波である．従来手法 1，2 に対して，提案手 法による表面反射波抑圧性能が改善していることが確認される．また，従来手法 2 では癌細胞応答が若干抑圧されているのに対し，提案手法では癌細胞応答が残って いることが確認できる．一方，提案手法の θ = 1, 3 付近の表面反射波が，他の素子 に比べて残差が存在することが確認できる．これは，非線形問題である式 4.3 での 解が同素子位置において局所解に陥った可能性が考えられる．同様に，CASE2 で の波形データの比較を図 4.5 に示す．均質モデル同様，提案手法による癌細胞応答 の抑圧を避けた表面反射波抑圧性能の向上が確認できる．また，CASE2 の提案手

表 4.1: 各 CASE，手法での定量評価 従来手法 1 従来手法 2 提案手法 RMSE[mm] CASE 1 2.1 0.6 CASE 2 2.6 0.5 PPRR[dB] CASE 1 -28.5 -38.1 -42.8 CASE 2 -21.8 -32.6 -37.8 法の部分微分波形最適化では局所解に陥らず，何れの素子でも精度よく抑圧を行え ていることが確認できる．次に，乳房境界推定精度と表面反射波抑圧性能の定量 評価として算出した RMSE(root mean square errors) と PPRR(The Peak-to-peak Response Ratio) を表 4.1 に示す．RMSE は，境界抽出精度を評価するもので，次 式で定義する． RMSE = v u u t 1 M M ∑ m=1 min pest ||p est,m− ptrue,m||2 (4.6)

ここで，pest,m, ptrue,mはそれぞれ推定点と真値，M は Envelope 法による推定点数

である．PPRR は，表面反射波抑圧性能を評価するもので次式で定義する． PPRR = 1 N N ∑ i=1 log₁₀ ( max ˜s(ri; t) max Eobs_(r i; t) |Tr≤t≤Tw ) (4.7) 提案手法によって高精度な乳房境界情報が推定でき，部分微分波形最適化にって抑 圧性能が 10dB 程度改善していることが確認できる．