バイオイメージングレポート課題（超音波イメージング、東担当分） 2016夏 日本語版 バイオイメージングレポート課題（超音波イメージング、東担当分）　2016夏 日本語版 　超音波撮像の原理を理解するには、手を動かして計算する必要があるが、講義では演習時間を設けられなかったので、以下の２題をレポート課題とする。数値計算の経験が無い人にも解けるように、appendixにエクセルでの計算例を説明するので、参考にしてください。勿論、matlab, mathematica, Cなど親しんでいるプログラムがある人はそれを使っても良い。（難しい場合は、機械、電気、精密系など数値計算になじみのある学生もしくは、私宛（azuma@fel.t.u-tokyo.ac.jp）に質問して下さい。） ・〆切：8月12日　提出方法：上記メールアドレスまでメール提出（提出後2日以内には受け取ったことを返信メールするので、それが無い場合は御手数ですが再度確認のメールをください） 課題１：ビーム計算 　二つの音源（A,B）から放射された超音波が形成する波面の様子を図１に示す。図中で赤と青の実線はそれぞれ音源A及びBから波長の整数倍の位置を示し、赤と青の破線はそれぞれ音源A及びBから半波長の奇数倍の位置を示す。これらは、音源AとBから同一周波数f[Hz]、同振幅、同位相でcos ωt の波が放射されている場合の、放射開始から周期Tの整数倍の時間が経過したときの波面（実線がa圧力が1となる位置、破線が圧力が-1となる位置）に等しい。ここで周期T=1/f、角周波数ω=2πfである。 　次に、二つの音源からの合成波の振幅分布を考えるために、代表的な点として、二つの音源からの距離が共に波長の整数倍の位置（黒丸）、共に半波長の奇数倍（緑丸）、片方からの距離が波長の整数倍、もう片方からの距離が半波長の奇数倍（白丸）の位置を示す。前述のように実線と破線でそれぞれ振幅が1と-1であるから、黒丸、緑丸、白丸の位置それぞれで、合成波の振幅は2,-2,0となる。この結果を、観測面上での位置と音圧のプロットを行うと図２のようになる。 音源B 音源A nλ (n+1)λ (n+0.5)λ 観測点　1, 2, 3, 4, 5 (n+1.5)λ 観測面 観測点　1 2 3 4 5 図２：合成波の圧力分布（観測面上） 図１：二つの音源から放射されている超音波の波面

A（x1） x1 A（x2） x2 波を複素数で考える手法の導入： 図２では、ある時刻における圧力分布を説明した。実際には下記に示すように、波は伝搬と共に、圧力は振動しており、ある瞬時の圧力を観測したのみでは、その点における真の超音波のエネルギーを調べたことにはならない。（波の圧力P（x,t）＝A（x）cos（wt＋θ）としたときに、場所毎の波の強度分布はA（x）であるが、例えばwt＋θ＝π/2となる時刻tにPを観測したら、A（x）に関わらずP＝０となってしまう。） 時刻ｔ 観測点１での圧力 の時間変動 観測点２での 圧力の時間変動 A（x1） A（x2） x1 x2 ここで、複素数を導入する。 複素平面 実部 虚部 時刻ｔ A（x）cos（wt＋θ） 実部のみ考えると 波の圧力は時間変動しているが 実部と虚部、共に考えると 複素平面 実部 虚部 時刻ｔ A（x）cos（wt＋θ） A（x）sin（wt＋θ） {cos（wt＋θ）}2＋{sin（wt＋θ）}2は常に１になるので、実部と虚部の自乗和は{A（ｘ）}2となる。 このように、波を実部に加えて、位相が90度ずれた虚部も加えて考えると、時間で振動する 効果を考えずに、圧力の空間分布を調べることが可能となる。

x L 以下、複数の音源からの波面の合成により、超音波ビームを計算する手法を説明する。 原点 　以下、複数の音源からの波面の合成により、超音波ビームを計算する手法を説明する。 ①音源を半径Lの円弧上に配列し、その円の中心を原点とする。音源から円の中心までは等距離になるので、円の中心が幾何学的な焦点位置となる。（以下、ある一つの音源の位置を(Xs,Ys）とする。 ②N個の観測点に関して観測点番号に対応するx座標の 位置を計算。 ③一つの音源から各観測点への距離を計算。 ④上記の求めた距離を、波長（=音速/周波数）で割り、位相を計算。波の実部（cosθ）、虚部(sinθ)を計算。 ⑤異なる音源に関しても、上記の③と④のステップを計算。 ⑥上記で計算した各音源からの波の実部と虚部それぞれの合成を行い、合成波の実部と虚部を計算。 ⑦上記の合成波の絶対値（実部と虚部の二乗和の平方根）を計算し、ビームパターンをプロットする。 原点 観測点：N点 dx dθ 音源配列 （半径Lの円弧上） 口径 図３：音源と観測点 ヒント：上記のXs,Ysは、音源の角度間隔をdθとすると、中心からN番目の音源ではXs = L sin((N-1)dθ), Ys = L sin((N-1)dθ) 以下、具体的な計算課題を示す。 A:音源数2と5のそれぞれ場合に関して、（周波数f:3MHz：、音速v：1540m/s、幾何学焦点の距離L：0.1m、観測点数N:32、観測点間距離dx：0.1mm（=1e-4m）、音源間隔は8°）の条件にて計算を行い、ビームがどのように変化するかを、計算結果を用いながら説明せよ。 B:素子数が一定の場合、①間隔を広げて全体の口径（図参照）を広げた場合と、②全体の口径は狭いが間隔が密に配列した場合に、それぞれ観測面上での超音波ビームがどのようになるか議論せよ。 C:今回の検討では、素子幅は十分に狭いものとして、また素子からの伝搬距離に応じた減衰や回折の考慮もしていない。これらの項を考慮すると、ビームはどのように変化すると考えられるか議論せよ。（グラフをプロットする必要は無い） D:もし伝搬経路において、音速が部分的に不均一であった場合にどのような影響が生じるか。考えられることを議論せよ。（グラフをプロットする必要は無い）

課題２：生体内での多重反射によるアーチファクトと、深部の反射体からの信号の識別に関する問題 図１に示すようなイメージングデバイスで、撮像対象内のある断面での撮像を行う場合、もし図２に示すように、撮像面のみを切り出して体外で撮像が可能であれば、反射体の反射強度が大きければ大きい程、イメージングが容易である。 　超音波エコー撮像の場合、音速と密度の積である、音響インピーダンスZが不連続となる界面で反射が生じる。では、生体内で音響インピーダンスの違いが大きいほど、撮像が容易なのか、以下の問題を解きながら考察を行おう。 実際の撮像では図２のように切り出すことが出来ないので、体外から超音波を照射し、撮像面まで音を伝搬させる必要がある。 A 撮像対象 反射体 イメージングデバイス 図1 B A B 撮像面 図2 時間 エコー強度 図３ 反射源１ 反射源２ 送信/受信器 エコー2 エコー1 エコー3 エコー4 そこで、図３に示すように音の伝搬経路に沿った、一次元の問題を考える。撮像面内に音響インピーダンスに不均一がある（＝反射体がある）と考える場合、当然、伝搬経路上にも音響インピーダンスに不均一があると考える方が妥当である。 　そこで図の左から右側に超音波を送信し、対象物からのエコー信号を、送信と同じ左側で受信する場合に、二つの音響インピーダンスが異なる領域（図中の反射源１と反射源２）からのエコーに関して検討する。反射源１に関しては、左側の境界で反射するエコー１、左側の境界を透過、右側で反射、再び左側を透過するエコー２、反射源内を2往復してから出てくるエコー３などが存在する。エコー１と２はイメージング上有用であるが、エコー３はアーチファクトとなる。エコー３とエコー４の強度比が、実際のイメージングにおいて重要となるので、以下このエコー３とエコー４の強度比を検討する。 反射率Rに対して、透過率T＝1-Rとして、伝搬中の減衰や散乱を考慮しないとエコー３の強度ははT2R3に比例し、エコー４の強度はT4Rに比例する。 課題：①反射率を0~1の範囲で変えて、エコー３と４の強度をそれぞれプロットせよ。この結果からエコーでのイメージングに適した反射率の範囲を議論しなさい。 ②音響インピーダンスz1とz2の界面での反射率はR＝|z1-z2|/(z1+z2)である。脂肪、水などの生体を構成する物質の密度と音速を調べ、その積から音響インピーダンス、水が多い領域と脂肪が多い領域の界面での反射率を計算し、上記①の結果と比べて議論しなさい。 ③実際のイメージングと上記の検討の違いに関して検討し、議論を行いなさい。

Appendix： Excelでの計算ステップの例： １．観測点番号を配列として定義 ②左記の①で作成したセルの一つ下のセルを 選び（例ではC列3行）、ここに「=C2+1」と入力する。 もしくは図中の赤丸の部分に式を入力しても良い。 注：上記C2は①で作成したセルに対応させる。 （選択したいセルをクリックしても良い。）この状態でエンターキーを入力すると、下図のようになる ③前の②のステップで入力したセル の右下の角をマウスでつかみ、 角を掴んだまま 下におろすと、 ④観測点の番号が上から順に入れながら、 必要な観測点の数だけセルを作る。

対応するB列の説明に記述しているだけで、ここではA列とB列の間に対応関係はリンクさせていない） ２．観測点のｘ座標の計算例 ①定数の設定：　B列2（N=32）や、B列3(dx=0.2mm=2e-4m)に示すように計算で繰り返し使う定数をセルにいれておく。（A列2や、A列3は、 対応するB列の説明に記述しているだけで、ここではA列とB列の間に対応関係はリンクさせていない） ②観測点のｘ座標の計算：　観測点番号Nに対して X = (N-32/2-0.5) × dx　を図中のセルD列2行目のように =(C2-B2/2-0.5)*B3　と入力する。 このままセルの右下角をマウスでつかんで下に引っ張ると 観測点の番号（C列）に応じて、x座標（D列）の値が連動して 計算される。 ③しかし、同時に観測点数（N=32）と、観測点間の 距離（dx=0.2mm）が格納されたセル番号も変化してしまう。そこで連動させたくない変数はB2を$B$2のように変更する。（同様にB3も$B$3）その上でセルD2の右下角を掴んで、下に引っ張れば、観測点の番号が入ったC列の値 だけが連動して変化して、D列の計算が出来る。 ２．観測点と、音源間の距離、位相、実部、虚部の計算例 ①他の定数、周波数、音速、波長、幾何学焦点の距離、音源間隔などを 設定し、音源の位置座標を計算する。 （例えば、このように音源のｘ座標=B8*SIN(B9/180*3.1415) 、ｙ座標 =B8*COS(B9/180*3.1415) ）のように計算。円周率πはPI()としても良い。 ②音源と観測点の座標が計算されたら、二点(x1,y1)と(x2,y2)間の距離 を計算（E列）、波長で割って位相θを求め （F3=E3/$B$7*2*3.1415 ）、実部（G3=cos(F3)）、虚部（H2=sin(F3)）を計算して、 これまでと同様にそれぞれのセルの右下を掴んで下に引っ張れば、 観測点ごとの、音源１からの波の実部、虚部が計算できる。 (𝑥1−𝑥2) 2 + (𝑦1−𝑦2) 2 ③音源一つ毎に、E,F,G,H列の 計算をしたのと、同様の操作を 行う。

３．波の合成の計算例 ①各音源からの波の合成は、実部、虚部それぞれに関して、図のO（O3=G3+L3）、P列（P3=H3+M3） のように計算し、絶対値（実部と虚部の二乗の和の平方根）を計算する。 最後にQ列、観測点番号1~32までをプロットすると、ビーム計算が完成する。