携帯電話に対する人体ドシメトリの計算機シミュレーション
王建青藤原修名古屋工業大学 電気情報工学科

１．まえがき携帯電話はその利便性のゆえに爆発的に普及しており，その勢いは留まることを知さららない．一方，携帯電話使用時には枢要器官の頭部は局所的に強い電磁界に曝されることになり，その安全利用を目的とする局所吸収指針が世界各国において相次いで制定されている．一般に電波と呼ばれる高周波電磁波の人体影響は，内部組織の温度上昇に起因するものとされ，携帯電話を対象とした上述の安全指針もこれを共通のベースとしている．そのとき，評価尺度としては単位体重当たりの吸収電力，すなわち SAR（Specific Absorption Rate，単位 W/kg）が用いられ，SAR を定量することを「ドシメトリ（Dosimetry）」という．しかしながら，人体内部の SAR は直接測定が不可能なため，ドシメトリ評価には数値モデルによる計算機シミュレーションかファントムと呼ばれる擬似人体モデルでの測定に頼らざるを得ない．携帯電話のドシメトリの計算評価には人体頭部を解剖学的に詳細に模擬した数値モデルが用いられ，そのドシメトリのシミュレーションも数値モデルの分解能の向上とともに大規模となり，大型計算機の利用が欠かせないものとなっている．
本文では，携帯電話に対する人体頭部ドシメトリの計算手法の概要を述べ，東北大学情報シナジーセンターのスーパーコンピュータを利用した筆者らのシミュレーション例[1]を紹介する．また，プログラムのベクトル化と並列化の向上手法についても述
べる．

２．計算法
ドシメトリの指標となる SAR は，生体組織内の電界を E，導電率をσ，密度をρとす
れば，SAR= σ E / 2 ρ で計算される．従って，携帯電話による頭部ドシメトリの計算
2

評 価 は ， 結 局 生 体 組 織 内 の 電 界 E を 求 め る こ と に 帰 着 し ， 一 般 に FDTD
（Finite-Difference Time-Domain）法[2]が用いられる．FDTD 法とは，電界 E，磁界H に関するマクスウェルの方程式を時間領域と空間領域とで差分化し，その差分式を時間領域で逐次計算することで計算領域内の電磁界を数値的に求める手法をいう．電界 E，磁界 H に関するマクスウェルの方程式について FDTD 定式化を行うと，例えば
電界 E と磁界 H のｚ成分はそれぞれ H yn +1/ 2 (i, j , k ) − H yn +1/ 2 (i − 1, j , k ) ∆t / ε 
1 − σ∆t / 2ε n∆xEzn +1 (i, j , k ) =Ez (i, j , k ) +1 + σ∆t / 2ε
1 + σ∆t / 2ε  H xn +1/ 2 (i, j , k ) − H xn +1/ 2 (i, j − 1, k ) −∆y E yn (i + 1, j , k ) − E yn (i, j , k ) ∆t ∆xn +1/ 2n −1/ 2H z (i, j , k ) = H z (i, j , k ) +µ  Exn (i, j + 1, k ) − Exn (i, j , k ) −∆y
(2)で表せる．ここで，式(1)，式(2)のµ，ε，σはそれぞれ媒質の透磁率，誘電率，導電率である．空間領域における差分化として，計算対象を三次元の微小立方体か直方体（セル）に分割し，各セルに生体組織と対応する電気定数（誘電率εと導電率σ）を割り付ける．
∆X∆Z
HxHyEzEyHz∆YCALL SETUP
T=0
DO N=1, N_ TI MESTEPCALL E_ FI ELD_ CALCUL ATI ONCALL E_ FI ELD_ PML
T=T+DT/ 2
CALL H_ FI ELD_ CALCUL ATI ON
CALL H_ FI ELD_ PML
T=T+DT/ 2
END DOEx

図１

FDTD 単位セ ル内の 電磁界 配置

図２

FDTDプ ログラ ム例

図１は単位セル内の電磁界各成分の空間配置を示す．基本的には電界はセルの各辺に沿って，磁界は面の中心に垂直に割り当てられ，電界と磁界は空間的に相互に配置される．これは，電界の回転が磁界を，磁界の回転が電界を作るというマクスウェルの方程式を満たすような配置となっている．一方，時間軸においては電界と磁界は時(1)間的に相互に配置されることになる．例えば，電界を t=n∆t (∆t: 時間ステップ)の整数次の時刻に，磁界を t=(n+1/2) ∆t の半整数次の時刻にそれぞれ割当て，t=(n-1) ∆t の電界 En −1と t=(n-1/2) ∆t の磁界 H

n −1/ 2nとから E を， Hn−1/ 2nと E とから H

n +1/ 2をそれぞれ計算する，というように電界，磁界が順次計算される．特に，式(1)，(2)からもわかるように，(i, j , k ) での電界を計算するためにはその前の∆t/2 時刻の隣接磁界，同様に (i, j , k ) での磁界を計算するためにはその前の∆t/2 時刻の隣接電界が用いられるた
め，FDTD 法は大型計算機の得意とする大規模並列計算に特に適しているといえる．
なお，計算機で取り扱える解析領域は有限であるため，解析領域を仮想な境界で閉じておく必要がある．この仮想的な境界を吸収境界といい，その条件を吸収境界条件という．代表的な吸収境界条件は，境界に仮想的な吸収媒質を置いて入射波を減衰させる PML（Perfectly Matched Layer）である．
FDTD 計算のプログラム例は図 2 に示す．

３．頭部内 SAR のシミュレーション例

図３

頭 部数値 モデル

図３は日本人成人男性の頭部 MRI(Magnetic Resonance Imaging)データから筆者らが製作した頭部数値モデルを示す．MRI データは，各水平断面において 256ｘ256ピクセル（約 1 mm 四方）の空間分解能，9 ビットグレースケールの濃淡分解能を有する．この MRI データを基に，各ピクセルを放射線医師の指導の下で，17 種類のRGB(Red-Green-Blue)コードのいずれかに指定することで皮膚，脂肪，筋肉，骨，脳など 17 種類の組織を同定した．こうして得られた頭部数値モデルは，一辺 2mm の立方体セルを約 53 万個集積して構成されている．
図 4 は携帯電話機による頭部垂直断面内の SAR 空間分布を示す．なお，図(a)は携帯電話機のアンテナが伸張時，図(b)は収納時のものである．図から，SAR の最大値は携帯電話機側の頭部表面耳付近で生じ，そこから遠ざかるに従って SAR は減衰し，頭部内部にはホットスポットは形成されないことがわかる．また，アンテナ収納時には短いアンテナ部及び筐体上方に電流が集中した結果，高 SAR 領域は耳付近に集中する
ば く ろが，アンテナ引き出し時にはそれが頭部上方に分散され，曝露領域が広くなるもののピーク SAR 値は低くなる．
このように，頭部内ＳＡＲ分布をシミュレーションすることで，人体に対する安全性評価を行うことが可能となる．

2 n [W/k g]n=5n=−7(a)(b)

図４ 携 帯電話 による 頭部垂 直断面 内のS AR空間 分布(a)アンテ ナ伸張 時；(b)アンテ ナ収納 時

４．FDTD コードの性能向上
前述ドシメトリで用いた筆者らの FDTD コードの性能向上にあたっては，東北大学情報シナジーセンターから有益な助言と指導を頂いた．この FDTD コードは当初
Fortran 77 で作成されたもので，PC Unix 上で使用することを想定し，演算は全て単精度で行われている．しかし，スーパーコンピュータ(NEC SX-4)は演算を倍精度で行うために，型宣言が単精度の場合には，単精度⇔倍精度間の変換が必要になり，単精度版のままでは演算が遅くなる．これに対処するために，プログラムの先頭に全ての変数を倍精度に指定するように
IMPLICIT REAL*8 (A-H,O-Z)を追加することにした．
FDTD 計算における電界と磁界の計算部は，式(1)と(2)からもわかるように高ベクトル化率と高並列化率が容易であるが，PML 吸収境界条件の部分は，一般に PML の層
K(1,NY,K)

JI(1,NY,NZ)(1,1,K)

L=3(b)(NX,1,NZ)(NX,NY,K)
L=2

L=1(NX,NY,NZ)(1,1,NZ)

L=4(NX,1,K)

PM L層(NX,NY,1)(1,J,NZ)(1,1,1)

L=6(NX,J,NZ)(NX,1,1)
L=2L=1(a)(1,J,1)

L=5(NX,J,1)(c)
図５ FDTD 解析空 間とPML 各層の 配置(a) 全FDTD 解析空 間；(b) I J平面断 面図； (c) I K平面 断面図

数 M が 10 前後（筆者らのプログラムでは M=12）しかなく，ベクトル化と並列化においては工夫を要する．図５に FDTD 解析領域と PML 媒質の配置を示す．PML 媒質層 は 解 析 領 域 を 囲 ん で L=1 か ら ６ ま で 計 ６ 面 が あ る ． 図 ６ に 電 界 に 対 す る Ｅ_FIELD_PML の例を示す．このプログラムは文献２で紹介されたもので，広く用いられているようである．PML 媒質層（図中グレーの領域）が解析領域を囲んで６面があるので，各 PML 媒質の始点の座標，終点の座標を LPMLII(L,1)，LPMLJJ(L,1)，．．．
，LPMLKK(L,2)によって定めている．ただし，L=1 は，I=1 に接する PML 層，L=2 は
I=NX に接する層，同様に，例えば，L=6 は K=NZ に接する層としている．このプログラムは見かけ上非常にコンパクトに書かれているが，次のような問題点がある．
（１） ベクトル化率を向上させるためには，最内側のループの長さが最大になるようなループ構成にする必要があるが，図６のプログラムでは L=5 と 6 のときに最内側の K に関するループは 11 か 12 の長さしかなく，ベクトル長は極めて短い．例えば，E_FIELD_PML の平均ベクトル長は理想値 256 に対して僅か 31であった．
（２） 並列化の効率を向上させるためには，できるだけ外側のループで並列化する必要があるが，図６のプログラムでは，最も外側の L のループは長さが６と短い．
これでは，例え 16 や 32 個の CPU で並列演算を行わせても並列化の効果は実質的に発揮されない（後述の［備考］参照）．

DO L=1, 6
I 0=LPML I I ( L, 1 )I 1=LPML I I ( L, 2 )J 0=LPML J J ( L, 1 )J 1=LPML J J ( L, 2 )K0=LPML KK( L, 1 )K1=LPML KK( L, 2 )L1=LPML ST( L)
DO I =I 0 , I 1- 1
DO J = J 0+1, J 1- 1
DO K=K0+1 , K1- 1
EX_FI EL D_CALCULATI ON
L 1=L1+1
END DOEND DOEND DO
L2=LPML ST( L)
DO I =I 0 +1, I 1- 1
DO J = J 0, J 1- 1
DO K=K0+1 , K1- 1
EY_FI EL D_CALCULATI ON
L 2=L2+1
END DOEND DOEND DO
L3=LPML ST( L)
DO I =I 0 +1, I 1- 1
DO J = J 0+1, J 1- 1
DO K=K0, K1- 1
EZ_FI EL D_CALCULATI ON
L 3=L3+1
END DOEND DOEND DOEND DO

図６

E_FIELD_PMLの プログ ラム例(文献２より)

DO K=2, M
DO J =2, NY- 1
DO I = 1, NX- 1
EX_ FI ELD_ CALCUL ATI ON
END DOEND DOEND DO
DO K=M+1, NZ- M
DO J =2, M
DO I = 1, NX- 1
EX_ FI ELD_ CALCUL ATI ON
END DOEND DODO I =1, M
DO J = M+1, NY- M
EX_ FI ELD_ CALCUL ATI ON
END DOEND DODO I =NX- M, NX- 1
DO J = M+1, NY- M
EX_ FI ELD_ CALCUL ATI ON
END DOEND DODO J =NY- M+1, NY- 1
DO I = 1, NX- 1
EX_ FI ELD_ CALCUL ATI ON
END DOEND DOEND DO
DO K=NZ- M+1, NZ- 1
DO J =2, NY- 1
DO I = 1, NX- 1
EX_ FI ELD_ CALCUL ATI ON
END DOEND DOEND DO
以下EY，EZ同様

図７

E_FIELD_PMLの プログ ラム例(改良版)実際に，SX-4 でのトレース解析結果によれば，FDTD 全計算時間の 6 割が PML の部分で費やされていることがわかった．この問題を対処する手段として，まず，変数 Lをなくし，従来の L，I，J，K に関する４重ループを I，J，K に関する３重ループに変更した．その結果，プログラム上では I，J，K に関する３重ループが従来の３つから 18 になる．つぎに，ベクトル化率を向上させるために，最内側のループ変数をループ長の長いものにした．即ち，L=1 と 2 に対応する PML 層については J を最内側のループ変数，L=3〜6 に対応する PML 層については I を最内側のループ変数にし，ループ長が PML の層数 M しかないような変数を最内側のループにしない．さらに，並列化の効率を上げるために，最外側のループをできるだけループ長の長いものに変更した．このような考え方に基づいた電界に関する PML のプログラム例を図７に示す．
なお，磁界に関する PML も全く同様である．
表１にベクトル化と並列化の向上効果を示す．このとき NX=NY=NZ=251，時間ステップ数 N_TIMESTEP=100 とした．表から，プログラム全体の平均ベクトル長は 75.8から 231.2 に３倍以上（PML だけの平均ベクトル長が 242 まで）向上したこと，またこのベクトル化の向上で演算時間が 51.1%に低減されたことがわかる．さらに，従来のプログラムでは L が最大６なので，１6 並列しても８並列時と同様な演算時間を要したが，プログラムの改良で並列化の効果が向上し，改良前に比べて８並列では 74%，16 並列では 82％の演算時間の短縮ができた．
表１

実 行時間 と平均 ベクト ル長

並列せず

従来
8並列

16並列並列せず

改良後
8並列

16並列

実行時間（相対値）

100.0

27.6

27.7

51.1

7.2

4.9

平均ベクトル長

75.8

75.8

75.8

231.2

231.2

231.2

［備考］
図６のプログラムは，自動並列化のオプション-Pauto だけではＩのループでの並列化となり，粒度（並列実行される部分の時間）が小さくて並列化の効率は悪い．そこで，つぎの並列化指示行を「DO L=1,6」の直前に指定して，Ｌのループで並列化し，粒度を大きくすることもできる．ただし，この場合でも，Ｌの長さは６なので，６個の CPU しか使われない．
!cdir paralleldo private(I0,I1,J0,J1,K0,K1,L1,L2,L3)

ここで，「private( )」は並列化に際して CPU ごとに確保すべき変数を指定したものである．

５．むすび
スーパーコンピュータを利用した携帯電話に対する頭部ドシメトリに関する筆者らの FDTD 計算例を紹介し，PML 吸収境界条件のベクトル化と並列化の向上手法を述べた．FDTD 法は極めて並列化に向くアルゴリズムであり，スーパーコンピュータ上でその威力が一層発揮される．人体を解析対象とする研究分野では，人体数値モデルの分解能や組織数は年を追うごと高精度かつ大規模になっており，スーパーコンピュータの利用価値は一層高まるものと予想される．
謝辞
本研究での計算の一部は，東北大学情報シナジーセンターのスーパーコンピュータを使用して行われたものである．計算コードのベクトル化・並列化にあたっては，情報シナジーセンターから有益な指導と助言をいただいた．
文献[1] J. Wang and O. Fujiwara, “Dosimetry in the human head for portable
telephones,” in The Review of Radio Science 1999-2002, Edited by W.R. Stone,
Wiley-Interscience, 2002, pp.51-63.
[2] 宇野亨，FDTD 法による電磁界およびアンテナ解析，コロナ社，1998.