(57)【要約】      （修正有）【課題】物質内に流入する電界などの物の量を正確に把握することのできる物流入量算出方法を提供する。【解決手段】スカラーポテンシャルを未知数とする連立１次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われることを特徴とする。【選択図】なし【特許請求の範囲】
【請求項１】
スカラーポテンシャルを未知数とする連立1次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、
上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、
上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われる
ことを特徴とする物体内への物流入量算出方法。
【請求項２】
請求項１記載の方法であって、
上記の所定の物が電流である方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物体内への物流入量算出方法に関するものである。
【背景技術】
【０００２】
近年、中間周波数帯 (300 Hz &#8211; 10 MHz)の電磁界を利用した機器の利用拡大に伴い、機器から発生する電磁界による生体への影響に関心が寄せられている。準静的近似が有効な周波数帯における電界ばく露時の人体内誘導電界を計算する手法の一つとして SPFD (Scalar Potential Finite Difference)法が挙げられる（非特許文献１及び２参照）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】T.W. Dawson, J. De Moerloose, and M.A. Stuchly, “Com-parison of magnetically induced ELF fields in humans com-puted by FDTD and scalar potential FD codes,” Appl.Comput. Electromagn. Soc. J., vol.11, no.3, pp.63&#8211;71, 1996.
【非特許文献２】T.W. Dawson, K. Caputa, and M.A. Stuchly, “High-resolution organ dosimetry for human exposure to low-frequency electric fields,” IEEE Trans. Power Deliv., vol.13,no.2, pp.366&#8211;373, 1998.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上述の手法は、スカラーポテンシャルを未知数とする大規模な連立一次方程式を解くことに帰着されるため、連立一次方程式を高速に解ける数値計算ライブラリを用いて解くことで、解析時間の短縮が期待できるものの、SPFD 法の連立一次方程式は、右辺ベクトル b の全要素の和が 0 という条件を満たさない場合、解が存在しないという問題がある。この右辺ベクトルの値は人体表面から流れ込む電流から計算され、物理的に考えるとその総和は 0 になるはずである。しかし、この値を外部電界計算によって求める場合は、数値誤差やモデルの離散化による誤差などが含まれ、電流の総和は 0 にならないことが多い。
そのため、誤差が最も小さくなる妥協解を求めることが必要である。しかし、右辺ベクトルの要素和が 0 でないとき、解法によっては解が収束しないことや大きな誤差が出てしまうことがあり、数値計算ライブラリを適用出来ない場合がある。
従って本発明の目的は、物質内に流入する電界などの物の量を正確に把握することのできる物流入量算出方法を提供することにある。
【課題を解決するための手段】
【０００５】
本発明は、以下の各発明を提供するものである。
１．スカラーポテンシャルを未知数とする連立1次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、
上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、
上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われる
ことを特徴とする物体内への物流入量算出方法。
この方法において上記外部環境データは、実測データでも計算値でもよい。また上記の所定の物の具体例としては、電流、気体又は液体等が挙げられる。
２．１記載の方法であって、上記の所定の物が電流である方法。
【発明の効果】
【０００６】
本発明によれば、物質内に流入する電界などの物の量を正確に把握することができる。
【図面の簡単な説明】
【０００７】
【図１】図１は、ボクセルの節点配置に関する説明図である。
【図２】図２は、回路網でモデル化された人体表面を示す概略図である。
【図３】図３は、回路モデルを示す概要図である。
【図４】図４は、人体モデルと微小ダイ大ポールとを示す模式図である。
【図５】図５は、誘電電界強度分布の結果を示すチャートである。
【図６】図６は、本発明の方法と他の方法との計算時間を示すチャートである。
【発明を実施するための形態】
【０００８】
以下、本発明についてさらに詳細に説明する。
本研究では、右辺ベクトルの総和が 0 でない場合におけるSPFD 法の連立一次方程式の計算方法について検討を行う。まず、様々な数値解法に適用でき、誤差が最小となる解が得られるような右辺ベクトルの補正方法を提案する。そして、補正方法の妥当性を検討するとともに、実際に補正を用いて誘導電界解析を行い、様々な数値解法への適用可能性を検討する。
2. 1 SPFD 法 [5]
SPFD 法を用いて人体内誘導電界を計算する場合、2 段階の解析を行う。はじめに、電界ばく露時に人体表面に誘起される表面電荷を計算するため、人体外部の電界分布を計算する。次に、得られた表面電荷を境界条件として用い、SPFD 法で体内誘導電界を計算する。SPFD 法では、人体をボクセルでモデル化し、生体内のボクセルに含まれる全ての節点について、スカラーポテンシャルを未知数とする以下の方程式を立てる。
式 (1) は図１に示す節点0とその節点に隣接する6つの接点r(r =1&#8764;6)のスカラーポテンシャルからなる方程式となる。
なお、Srは節点0&#8722;r間のアドミタンス、ψr は節点 r におけるスカラーポテンシャルを表している。すなわち、SPFD法は図２のような回路網にモデル化される。また、境界条件は外部電界計算によって求められた表面電荷密度分布ρs を用いて、電流連続の式により式 (2) と表される。
なお、&#710;nはモデル表面に垂直な法線単位ベクトルJ は電流密度を意味する。この境界条件は図２の回路網では電流源として表される。
【０００９】
2. 2 SPFD 法の連立一次方程式
SPFD 法の連立一次方程式を行列形式Ａｘ＝ｂで表わせば、数値計算ライブラリを用いた計算が可能となる。なお、この時の係数行列Ａはアドミタンス、右辺ベクトル b は表面から流れ込む電流から成り立つ。ここで、 SPFD 法における係数行列 A の行ベクトルの和及び列ベクトルの和は 0 であり、rank-deficientであるため、非正則である。係数行列が非正則である連立一次方程式では、解 x が一意に決まらない。また、この時の連立一次方程式は右辺ベクトルが係数行列 A の列空間に存在するときのみ解を持ち、右辺ベクトルがこの条件を満たさない時、解は存在しない。SPFD 法の連立一次方程式では、右辺ベクトルの全要素の総和が0であれば、右辺ベクトルは列空間に存在し、解は存在するようになる。行ベクトルの和が0であることからも、右辺ベクトルの全要素の和が0でなければ、連立一次方程式Ax =bを満たす解が存在しないことがわかる。また、SPFD 法における右辺ベクトルの全要素の和は表面から流れ込む電流の総和であり、理論的には電流連続の式を満たし、この総和は0になるはずである。しかし、右辺ベクトルの値は外部電界計算によって求めるため、数値誤差や離散化による誤差などが含まれ、多くの場合、この総和は 0 にはならない。
このような連立一次方程式に対して、反復法を用いて解こうとすると解法によっては解が発散してしまうことがある。1 行1 列削除し、係数行列を正則にして解けば解は収束するが、電流の総和が 0 からかけ離れている場合は得られる解に大きな誤差が含まれる。
そのため、このような問題に対して、連立一次方程式 Ax = bを満たす解ではなく、残差ノルム ‖r‖L2 = ‖b   Ax‖L2 を最小とするような解を妥協解として得る手法が考えられる。
【００１０】
2. 3 残差ノルム最小の解を得る方法
残差ノルム最小の解を得る手法として、Moore-Penrose の擬似逆行列 [6] や最小化残差法 [7] などがあるが、ここでは、解法にかかわらず右辺ベクトルの補正をすることで残差ノルム最小の解を得る方法を提案する。
連立一次方程式 Ax = b の残差ノルムが最小である時の残差ベクトルを rmin とし、その時の解ベクトルを x&#8727;とすると、これらは式 (3) を満たす。式 (3) を整理すると、式 (4) となる。
式 (4) は係数行列 A、右辺ベクトル b &#8722; rmin の連立一次方程式が残差ノルムを最小とする解 x&#8727;を持つことを表しているので、これを解けば残差ノルム最小の解を得ることができる。元の連立一次方程式 Ax = b から補正後の右辺ベクトル&#732;b を得るのに必要な演算は式 (5) のみなので、様々な解法に適用可能であると予想される。
2. 4 残差ノルムの最小値及び残差ベクトル
SPFD 法における N 元連立一次方程式 Ax = b の残差ノルム ‖r‖L2 の最小値を考える。残差ノルムの定義は式 (6) であり、残差ベクトルの i 番目の要素 ri は右辺ベクトル b、解ベクトル x の i 番目の要素 bi 、 xi と係数行列 A の i 行 j 列の要素aij を用いて、式 (7) のように表される。
式 (6) の最小値を求めるには式 (8) の関数 f (r) の最小値を求めれば良い。
次に右辺ベクトルの全要素の和が式 (9) であるときの SPFD 法の連立一次方程式における残差ベクトル r が受ける制約条件を考える。
【００１１】
【数１】
2019124595.tif 000002
【００１２】
ここで残差ベクトルの全要素の和は式 (10) のように表せる。
第 1 項は右辺ベクトルの要素和であり、今回は式 (9) としたので、Itotal となる。また、SPFD 法の係数行列は任意の行で式(11) を満たす。さらに、対称行列であるので、式 (12) も同時に満たす。
式 (9)、(12) から式 (10) は式 (13) となる。
次に、式 (13) の制約のもとで式 (8) の関数 f (r) の最小値及び、その時の残差ベクトル r をラグランジュの未定乗数法を用いて求める。このときのラグランジュ関数は式 (14) となる。
最小となる r を求めるためには式 (15)、(16) の連立方程式を解けばよい。
これを解くと、ri 、 λ はそれぞれ以下のように求まる。
以上より、残差ノルム最小となる残差ベクトル rmin、残差ノルム ‖rmin‖L2 はそれぞれ式 (19)、(20) となることが示された。
実際に計算を行う際は、右辺ベクトルの総和 Itotal と要素数 Nから残差ノルム最小となる残差ベクトル rmin を求め、右辺ベクトルを式 (5) のように補正すればよい。
【００１３】
2. 5 補正方法の妥当性評価
前節で提案した補正方法を図 3 の回路モデルに用いて実際に計算を行い。その妥当性を評価した。この回路はアドミタンス S1 &#8764; S6、電流源 Ia、 Ib、 Ic からなり、未知数を節点電位 V1 &#8764; V 4 とする。それぞれのアドミタンスをS1 = 1 S/m、 S2 = 2 S/m、 S3 = 3 S/m、 S4 = 4 S/m、 S5 =5 S/m、 S6 = 6 S/m とし、それぞれの電流値は総和が 0 でなくなるように、Ia = 1 A、 Ib = &#8722;2 A、 Ic = &#8722;1 A とした。この時の連立一次方程式は式 (21) となる。
【００１４】
【数２】
2019124595.tif 000003
【００１５】
【数３】
2019124595.tif 000004
【００１６】
このときの電流の総和は &#8722;2.0 であるので、式 (19)、(20) から残差ノルム最小となる残差ベクトルは rmin = &#8722;0.5[1 1 1 1]t、 残差ノルムは ‖rmin‖L2 = 1.0 となる。この連立一次方程式の右辺ベクトルに対して補正を行なった上で共役勾配法 (CG: ConjugateGradient Method) を用いて解く。また、 Moore-Penrose の擬似逆行列を用いて同様の連立一次方程式を解く。Moore-Penroseの擬似逆行列 A+と右辺ベクトルとの積 x = A+b によって得られる解 x は残差ノルム最小の解であるという特徴があるため、その時の残差ベクトル及び残差ノルムはそれぞれ、式 (19)、(20)を満たすことが予想される。この条件で計算を行った時、それぞれの手法によって得られる解 V1 &#8764; V4、電界 E1 &#8764; E6 を比較する。なお、電界値は節点間の電位差から求める。
表 1&#8211;3 に各手法で得られたスカラーポテンシャル、電界、残差ベクトル及び残差ノルムを示す。表 2、 3 から補正を用いて得られた電界、 残差ベクトル、 残差ノルムの値は Moore-Penroseの擬似逆行列を用いて得られた値とそれぞれ一致することが確認できた。したがって、補正により残差ノルム最小の解を求められることが確認された。
【００１７】
【表１】
2019124595.tif 000005
【００１８】
3. 数値人体モデルを用いた解析
3. 1 ばく露条件
図 4 のように数値人体モデル前方に波源として微小ダイポールを配置した場合の内部誘導電界計算を提案の補正方法とともに行った。波源の周波数は 85 kHz とし、数値人体モデルには、情報通信研究機構によって開発された日本人成人男性モデルTARO モデル [8] を用いた。本モデルは空間分解能 2 mm で51 種類の組織で構成されている。外部電界計算時の TARO モデルの材質は完全導体とし、内部誘導電界計算時の人体組織の電気定数は文献 [9] より引用した。外部電界計算には AET 社の有限積分法に基づく電磁解析シミュレーションソフト CSTSTUDIO SUITE Low Frequency Solver を用いた。このときに得られた電流分布の総和及び最小残差ノルムなどの値を表 4に示す。
3. 2 残差ノルム最小の解の妥当性
補正によって求められた残差ノルム最小の解の妥当性を検討した。外部電界計算によって得られた電流値に対して、式 (22)を行い、電流の総和 Itotal = 0 である電流分布を模擬した。
この電流分布を用いて得られる誘導電界を正しい値であるとし、補正によって求められる誘導電界と比較する。
図 5 に SPFD 法によって求められた内部誘導電界強度分布を示す。また、内部誘導電界の統計量は表 5 のようになっており、補正によって得られた値は概ね一致していた。
【００１９】
【表２】
2019124595.tif 000006
【００２０】
3. 3 様々な解法を用いた計算
数値計算ライブラリに含まれる連立一次方程式の様々な数値解法を用いて計算を行い、その時の計算時間を測定した。解法として、
&#8226; 共役勾配（CG: Conjugate Gradient）法
&#8226; 不完全コレスキー分解付き CG（ICCG: Incomplete  Cholesky decomposition CG）法
&#8226; 不完全 LU 分解付き CG （ILUCG: Incomplete LU CG）法
&#8226; 最小残差（MINRES: Minimize Residual）法
&#8226; 対称 LQ（SYMMLQ: Symmetric LQ）法
&#8226; 代数的多重格子法付き CG （AMG-CG: Alegebraic Multi  Grid CG）法
を用いる。CG 法、ICCG 法、ILUCG 法、MINRES 法、SYMMLQ 法は Mathworks 社の MATLAB 2016a の線形代数ライブラリを用いた。AMG-CG 法は代数的多重格子 (AMG)法 [10] を CG 法に適用した解法であり、みずほ情報総研によって開発された数値計算ライブラリを用いた。本ライブラリは、CPU による計算と GPU（Graphics Processor Unit）による並列計算の両方が可能であり、今回は CPU、GPU、両方の計算をそれぞれ行った。なお、CPU の計算では、計算の並列化は行わなかった。これらの計算を表 6 の計算環境で実行した。
収束条件は相対残差ノルム ‖b   Ax‖L2/‖b‖L2 < 1.0 × 10&#8722;5とし、各解法が収束条件を満たすまでにかかった時間を計測する。なお、計算時間の計測は前処理を含む各数値解法の解を求めるルーチンのみとした。
また、これらの解法に加え比較対象として、従来広く用いられている逐次過緩和（SOR: Successive Over Relaxation）法の計算も行った。SOR 法は、異なる収束基準を用いていることや収束性が緩和係数の値に依存することなどから、単純に上記の解法と計算時間の比較をすることは出来ないが今回は計算の一例として、緩和係数を 1.5 とし、50000 反復にかかった計算時間を比較対象とする。なお、50000 反復目の相対残差ノルムは3.75 × 10&#8722;2であったため、‖b   Ax‖L2/‖b‖L2 < 1.0 × 10&#8722;5上記の収束条件を満たすためにはさらに反復計算を行う必要がある。なお、SOR 法の計算は C 言語で記述し、bash 上で実行した。
それぞれの数値解法ごとの計算時間を図 6 に示す。AMG-CG法は一番計算時間が短く GPU による並列計算を行った場合は3.2 秒、並列化を行わなかった場合は 96.7 秒であった。次に計算時間が短かったのは前処理付き共役勾配法の ILUCG 法、ICCG 法の 2 つで両者とも計算時間は 1200 秒程度であった。これらの解法を SPFD 法の連立一次方程式のとしてを用いた場合、従来用いられてきた CG 法や SOR 法を用いた場合に比べて、より高速な解析が可能であることがわかった。
【００２１】
4. 結 論
本研究では、右辺ベクトルの全要素の和が 0 でない矛盾したSPFD 法の連立一次方程式に対する補正方法の提案と有効性の検討を行なった。提案した右辺ベクトルの補正方法を用いれば、残差ノルム最小の解を得ることができ、妥当な結果が得られることを確認した。また、右辺ベクトルの要素和が大きい時、解法によっては大きな誤差が出たり、発散する場合があったが、補正を連立一次方程式を解く時の前処理として行うことで残差ノルム最小の解を数値計算ライブラリに含まれる様々な数値解法で得られることが確認できた。また、AMG-CG 法、ILUCG法、ICCG 法などの解法を用いれば、従来手法より高速な解析が可能であることがわかった。
【図１】
2019124595.tif 000007
【図２】
2019124595.tif 000008
【図３】
2019124595.tif 000009
【図４】
2019124595.tif 000010
【図５】
2019124595.tif 000011
【図６】
2019124595.tif 000012