Dispersive FDTD Modeling of Human Body

Ha, Sang-Gyu; Jung, Kyung-Young

doi:10.5515/KJKIEES.2020.31.3.205

J. Korean Inst. Electromagn. Eng. Sci. 2020; 31(3):205-215

pISSN: 1226-3133, eISSN: 2288-226X

DOI: https://doi.org/10.5515/KJKIEES.2020.31.3.205

튜토리얼 및 리뷰/TUTORIAL & REVIEW

인체 분산 FDTD 모델링

하상규¹

, 정경영¹^,^*

Dispersive FDTD Modeling of Human Body

Sang-Gyu Ha¹

, Kyung-Young Jung¹^,^*

Author Information & Copyright ▼

¹한양대학교 전자컴퓨터통신공학과

¹Department of Electronics Computer Engineering, Hanyang University

^*Corresponding Author: Kyung-Young Jung (e-mail: kyjung3@hanyang.ac.kr)

© Copyright 2019 The Korean Institute of Electromagnetic Engineering and Science. This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Mar 07, 2020; Revised: Mar 25, 2020; Accepted: Mar 30, 2020

Published Online: Mar 31, 2020

요약

본 논문에서는 인체와 전자기파의 상호작용을 해석하기 위한 인체 분산 FDTD(Finite-Difference Time-Domain) 모델링에 대해 검토한다. 주파수에 따라 전기적 특성이 변하는 78개의 인체조직 분산 특성을 표현하는 Debye 모델, Lorentz 모델, Cole-Cole 모델, Davidson-Cole 모델, Havriliak-Negami 모델, 2차 복소분수함수 모델, 4차 복소분수함수 모델에 대해 논의하고, 분산 모델을 FDTD 알고리즘에 적용하는 방법을 검토한다. 또한, 인체의 기하학적 구조를 모델링하기 위한 방법으로 Virtual Family Tool의 인체팬텀과 미디어 기법에 대해 논의한다.

Abstract

Herein, we review the dispersive finite-difference time-domain(FDTD) modeling of the human body for analyzing its interaction with electromagnetic waves. We discuss the dispersive modeling of 78 human tissues using the Debye model, Lorentz model, Cole-Cole model, Davidson-Cole model, Havriliak-Negami model, two-pole complex rational function(CRF) model, and four-pole CRF model. Subsequently, we review a method that applies the dispersion model to the FDTD algorithm. Moreover, we discuss the computational human phantom of the Virtual Family tool and the media method for the geometric structure of the human body.

Keywords: Finite-Difference Time-Domain(FDTD); Dispersion Model of Human Body

Ⅰ. 서 론

전자기파 수치해석 기술은 안테나, 레이다 및 마이크로파 시스템을 포함한 다양한 전자파 문제의 해를 매우 정확히 예측할 수 있게 빠르게 발전하였다. 다양한 수치해석법 중 FDTD(Finite-Difference Time-Domain)기법은 전자기파 산란 문제를 해결하기 위해 K. S. Yee에 의해 1966년 처음으로 도입되었고, 정확성과 견고성으로 인해 RF에서부터 광학 주파수에 이르는 전자기 현상을 해석하는데 널리 사용되고 있다^[1]～^[3]. 주파수 영역 해석기법인 모멘트법(MOM: moment of method)^[4]과 유한 요소법(FEM: finite element method)^[5]과 달리 시간 영역 해석기법인 FDTD는 한 번의 해석으로 광대역 주파수 응답을 얻을 수 있으며, 행렬계산을 요하지 않기 때문에 병렬처리에 용이하다. 또한, 맥스웰 회전 방정식(Maxwell Curl Equation)을 직접 이산화하기 때문에 복잡하고 다양한 경계조건을 해석하기에 매우 용이하다. 현재까지 FDTD 기법은 지속적으로 개발되고 있으며, 정보통신기술, 생명공학기술, 나노기술^[6], 그리고 에너지 기술^[7]에 이르기까지 다양한 응용분야에서 활용도가 매우 높다. 한편, 정보통신 기술과 생명공학 기술을 융합으로 전자기파를 이용한 유방암 진단, 전자파 온열치료, 인체 삽입 안테나 및 무선 캡슐 내시경 등과 같이 검진과 치료를 위해 인체와 전자기파의 상호작용에 대한 활발한 연구가 진행되고 있다.

한편, 심장, 위, 근육, 지방, 피부 등과 같이 다양한 인체조직은 주파수에 따라 복소 유전율이 변하는 분산 특성을 갖는다. 따라서 인체와 전자기파의 상호작용을 정확히 해석하기 위해 분산 특성을 고려한 인체 분산 FDTD 모델링과 인체조직의 기하학적 구조 모델링이 요구된다. 최근 인체조직의 분산 특성을 표현하는 모델들에는 Debye 모델^[8], Lorentz 모델^[9], Cole-Cole 모델^[10], Davidson-Cole 모델^[11], Havriliak-Negami 모델^[12], 2차 복소분수함수 모델^[13], 4차 복소분수함수 모델^[14]이 있다. 따라서, 본 논문에서는 인체조직의 분산 특성을 표현하는 분산 모델에 대해 논의하고, 분산 FDTD 알고리즘에 적용하는 방법에 대해 검토한다. 또한, 인체팬텀과 미디어 기법을 이용한 인체의 기하학적 구조 모델링 방법에 대해 논의한다.

Ⅱ. 수치해석 기법

2-1 Yee 알고리즘

FDTD 수치해석 기법은 유한 체적에 대한 맥스웰 회전 방정식을 시간과 공간 영역으로 이산화하고, 중심차분 근사법을 이용하여 차분방정식을 도출한다. 이러한 FDTD 기법은 그림 1과 같이 서로 다른 위치에 존재하는 이산화된 전계 및 자계의 공간격자를 사용하며, 양함수(explicit) 계산 형태로 전계는 이전 시간의 전계 및 자계에 의해 순차적인 형태로 계산된다^[1],^[2].

그림 1. | Fig. 1. Yee 공간격자에서 전계 및 자계 성분 | Electric and magnetic field of Yee grid.

Download Original Figure

선형, 등방성, 비분산 매질에 대한 맥스웰 회전 방정식은 다음과 같다.

∇ × H → = ϵ ∂ H → ∂ t

(1)

∇ × E → = − μ ∂ H → ∂ t

(2)

여기서 ϵ은 유전율, μ 투자율을 나타낸다. 식 (1)과 식 (2)의 벡터 방정식은 직각 좌표계에서 6개의 스칼라 방정식으로 나타낼 수 있다.

∂ E x ∂ t = 1 ϵ (∂ H z ∂ y − ∂ H y ∂ z)

(3)

∂ E y ∂ t = 1 ϵ (∂ H x ∂ z − ∂ H z ∂ x)

(4)

∂ E z ∂ t = 1 ϵ (∂ H y ∂ x − ∂ H x ∂ y)

(5)

∂ H x ∂ t = 1 μ (∂ E y ∂ z − ∂ E z ∂ y)

(6)

∂ H y ∂ t = 1 μ (∂ E z ∂ x − ∂ E x ∂ z)

(7)

∂ H z ∂ t = 1 μ (∂ E x ∂ y − ∂ E y ∂ x)

(8)

Yee 공간격자에 나타나 있는 전계 및 자계의 각 성분을 시간과 공간에 대해 중심차분 근사하면 암페어 법칙(Ampere’s law)과 패러데이 법칙(Faraday’s law)을 이산 수식으로 표현할 수 있다.

E x | i + 1 / 2, j, k n + 1 = E x | i + 1 / 2, j, k n + Δ t ϵ Δ y (H z | i + 1 / 2, j + 1 / 2, k n + 1 / 2 − H z | i + 1 / 2, j − 1 / 2, k n + 1 / 2) − Δ t ϵ Δ z (H y | i + 1 / 2, j, k + 1 / 2 n + 1 / 2 − H y | i + 1 / 2, j, k − 1 / 2 n + 1 / 2)

(9)

E y | i, j + 1 / 2, k n + 1 = E y | i, j + 1 / 2, k n + Δ t ϵ Δ z (H x | i, j + 1 / 2, k + 1 / 2 n + 1 / 2 − H x | i, j + 1 / 2, k − 1 / 2 n + 1 / 2) − Δ t ϵ Δ x (H z | i + 1 / 2, j + 1 / 2, k n + 1 / 2 − H z | i − 1 / 2, j + 1 / 2, k n + 1 / 2)

(10)

E z | i, j, k + 1 / 2 n + 1 = E z | i, j, k + 1 / 2 n + Δ t ϵ Δ x (H y | i + 1 / 2, j, k + 1 / 2 n + 1 / 2 − H y | i − 1 / 2, j, k + 1 / 2 n + 1 / 2) − Δ t ϵ Δ y (H x | i, j + 1, k + 1 / 2 n + 1 / 2 − H x | i, j − 1 / 2, k + 1 / 2 n + 1 / 2)

(11)

H x | i + 1 / 2, j, k + 1 / 2 n + 1 / 2 = H x | i, j + 1 / 2, k + 1 / 2 n − 1 / 2 + Δ t μ Δ z (E y | i, j + 1 / 2, j, k + 1 n − E y | i, j + 1 / 2, k n) − Δ t μ Δ y (E x | i, j + 1, k + 1 / 2 n − E z | i, j, k + 1 / 2 n)

(12)

H y | i + 1 / 2, j, k + 1 / 2 n + 1 / 2 = H y | i + 1 / 2, j, k + 1 / 2 n − 1 / 2 + Δ t μ Δ x (E z | i + 1, j, j, k + 1 / 2 n − E z | i, j, k + 1 / 2 n) − Δ t μ Δ z (E x | i + 1 / 2, j, k + 1 n − E x | i + 1 / 2, j, k n)

(13)

H z | i + 1 / 2, j + 1 / 2, k n + 1 / 2 = H z | i + 1 / 2, j + 1 / 2, k n − 1 / 2 + Δ t μ Δ y (E x | i + 1 / 2, j + 1, k n − E x | i + 1 / 2, j, k n) − Δ t μ Δ x (E y | i + 1, j + 1 / 2, k n − E y | i, j + 1 / 2, k n)

(14)

여기서 i, j, k 는 해석 영역에 대한 직각 좌표계의 x, y, z 공간 지표, Δx, Δy, Δz는 셀 간격, Δt는 이산 시간간격, n은 nΔt의 시간을 나타낸다. FDTD 알고리즘에서 전계와 자계는 식 (9)～식 (14)를 이용해 매 시간 갱신된다. 한편, FDTD 계산은 중심차분 근사로 이산화 했기 때문에 수치적인 분산 특성을 갖는다. 중심차분 근사로 발생하는 시간과 공간에서 오차는 Δx와 Δt가 작아짐에 따라 (Δx)²와 (Δt)² 만큼 감소하는 2차의 정확도(O[(Δx)²], O[(Δt)²])를 갖는다. 따라서 수치적인 분산 특성을 최소화하기 위해 Δx, Δy, Δz는 펄스의 최대 주파수에서 파장의 최소 λ/10보다 작아야 한다. 한편, FDTD 알고리즘은 양함수 형태를 가지고 있으므로, Δt는 다음의 안정조건을 만족해야 한다^[2].

Δ t ≤ 1 c 1 (Δ x) 2 + 1 (Δ y) 2 + 1 (Δ z) 2

(15)

여기서 c는 해석 공간에서 최대 전파 속도를 나타낸다. CFL(Courant-Friedrichs-Lewy) 상수는 다음과 같이 표현된다.

C F L ≤ c Δ t 1 (Δ x) 2 + 1 (Δ y) 2 + 1 (Δ z) 2 ≤ 1

(16)

CFL 상수는 FDTD의 수치적인 안정성을 위해 1보다 작거나 같아야 한다.

2-2 분산 FDTD 알고리즘

선형, 등방성, 주파수 의존 매질에 대한 맥스웰 구성 관계(constitutive relation) 방정식은 다음과 같다.

D → (ω) = ϵ (ω) E → (ω)

(17)

B → (ω) = μ (ω) H → (ω)

(18)

여기서 D는 전속 밀도, B는 자속밀도를 나타내며, 인체는 비자성 매질이기 때문에 μ(ω)는 μ₀로 정의된다. 식 (17)의 처리 방법에 따라 매질의 분산 특성 알고리즘은 PLRC(Piecewise Linear Recursive Convolution)^[15],^[16], RC(Recursive Convolution)^[17]～^[19]방법과 ADE(Auxiliary Differential Equation)^[20],^[21]방법으로 나뉜다. 먼저, 전계성분을 근사하기 위해 PLRC 방법은 컨볼루션 적분식을 하나의 구간 [mΔt, (m+1)Δt]에 대해 선형적으로 근사하는 반면, RC 방법은 계단형태로 근사한다. 따라서 근사의 정확도는 PLRC 방법이 RC 방법보다 우수하다. PLRC와 RC 방법은 다른 방법에 비해 적은 메모리를 요구하며, 다양한 분산 매질들을 하나의 알고리즘으로 표현할 수 있으나, 유도과정이 복잡하며 컨볼루션 적분식의 고유한 선형성 때문에, 비선형 분산 매질의 모델링이 불가능하다. 한편, ADE 방법은 적분식 대신 보조 미분 방정식을 사용하며, 푸리에 역변환을 이용해 주파수 영역으로 표현된 분산 모델을 간단하게 시간 영역으로 이산화 할 수 있다. ADE의 근사 정확도는 RC 방법보다는 우수하고, PLRC 방법과는 동일한 2차 정확도를 가진다. ADE 방법의 가장 큰 장점은 비선형 분산 매질을 쉽게 모델링할 수 있다는 것이다.

Ⅲ. 인체 분산 FDTD

3-1 인체조직 분산 데이터

인체조직의 유전율 특성에 대한 광범위한 데이터는 전자기 분야와 생물학적 시스템의 상호작용을 연구하는 응용분야에서 필요성이 크다. 1996년 Gabriel은 10 Hz～100 GHz의 주파수 범위에서 인체조직의 분산 특성을 4차 Cole-Cole 모델^[22]을 이용해 표현하였다. 주파수에 따른 도전율, 상대 유전율, 손실 탄젠트, 파장, 그리고 침투깊이가 제공되며, Aorta(대동맥)부터 Vitreous Humor(안구 유리체)까지 총 78개의 인체조직(표 1)에 대하여 데이터화 하였다.

표 1. | Table 1. Gabriel의 인체조직 데이터 | Gabriel’s human tissue.

No.	Human tissue
1	Air
2	Aorta
3	Bladder
4	Blood
5	Blood vessel
.	.
77	Uterus
78	Vitreous humor

Download Excel Table

3-2 Debye 분산 모델

도체 손실을 포함하는 Debye 분산 모델^[8]은 식 (19)로 표현된다.

ϵ r (ω) = ϵ ∞ + ∑ p = 1 P ϵ s, p − ϵ ∞ 1 + j ω τ p + σ 1 + j ω ϵ 0

(19)

여기서 ω는 각주파수, ε_s,p는 DC에서 유전상수, ε_∞는 무한 주파수에서 유전상수, ϵ₀는 진공에서 유전상수, τ_p는 시상수, σ는 도전율을 나타낸다. UWB(3.1~10.6 GHz)에서 인체조직의 분산 특성을 잘 표현하는^[23] 1차 Debye 모델은 변화율이 매우 큰 0.4～3 GHz 대역에서 정확성이 현저히 떨어진다. DC～1.5 GHz 대역까지 모델링 가능한 2차 Debye 모델^[24]은 도체 손실을 고려하지 않아 정확성이 떨어진다. 한편, 뉴턴기법과 최소 제급 기법을 이용한 비선형 최적화 문제를 풀어야 하는 도체 손실을 포함한 2차 Debye 모델은 부적절한 초기값을 선택하면 정확한 모델링을 할 수 없다^[25].

3-3 Lorentz 분산 모델링

Lorentz 분산 모델^[9],^[26]은 식 (20)으로 표현된다.

ϵ r (ω) = ϵ ∞ + ∑ p = 1 P (ϵ s, p − ϵ ∞) ω p 2 ω p 2 + 2 j ω δ p − ω 2

(20)

여기서 ω_p는 매질의 공진 주파수, δ_p는 감쇄계수를 나타낸다. Debye 모델과 마찬가지로, Lorentz 모델계수를 찾기 위해 복잡한 비선형 최적화 기법을 사용해야 한다. 각 계수의 크기는 보통 10²⁰배 이상 차이가나기 때문에 초기값 선택이 적절하지 않으면 정확한 인체 분산 모델링이 불가능하다.

3-4 Cole-Cole 분산 모델링

Gabriel의 연구로 데이터화된 인체조직의 분산 특성은 4차 Cole-Cole 분산 모델^[10]로 표현된다.

ϵ r (ω) = ϵ ∞ + ∑ p = 1 4 ϵ s, p − ϵ ∞ 1 + (j ω τ p) 1 − α + σ 1 + j ω ϵ 0

(21)

여기서 α는 0과 1사이의 계수이며, 이완현상에서 감쇄 요소로 유전율 곡선의 모양을 결정한다. 일반적으로 α는 정수가 아니므로 식 (21)은 감마함수(Γ)를 포함한 복잡한 적분으로 표현된다. 따라서 시간영역에서 인체 분산 특성을 표현하기 매우 어렵다. 또한, 인체 분산 FDTD 알고리즘 구현 시 전계성분의 업데이트 수식이 증가하여 해석 효율이 떨어지는 단점이 있다. 한편, Cole-Cole 모델 관계식에서 α가 0인 특수한 경우 식 (21)은 Debye 모델과 같다.

3-5 Davidson-Cole 분산 모델링

Davidson-Cole 분산 모델^[11]은 식 (22)로 표현된다.

ϵ r (ω) = ϵ ∞ + ∑ p = 1 p ϵ s, p − ϵ ∞ (1 + j ω τ p) β

(22)

여기서 β는 0과 1사이의 계수이며, 이완현상에서 비대칭 요소로 유전율 곡선의 너비를 결정한다. 일반적으로 β는 정수가 아니다. 따라서 1+jω는 분수 차수를 포함하므로 감마함수를 포함한 시간 컨볼루션 적분 때문에 인체 분산 FDTD 알고리즘 구현이 매우 어렵다^[27]. 한편, β가 1인 경우 도전율이 없는 Debye 모델과 같다.

3-6 Havriliak-Negami 분산 모델링

Havriliak-Negami 분산 모델^[12]은 아래와 같이 표현된다.

ϵ r (ω) = ϵ ∞ + ∑ p = 1 p ϵ s, p − ϵ ∞ [(1 + j ω τ p) 1 − α] β

(23)

여기서 α와 β는 0과 1 사이의 계수이다. α와 β가 1인 경우 Debye 분산 모델로 표현된다. β가 1인 경우 Cole-Cole 분산 모델로, α가 1인 경우, Davidson-Cole 분산 모델로 표현된다. 한편, Cole-Cole 분산 모델, Davidson-Cole 분산 모델 그리고 Harvriliak-Negami 분산 모델들은 더욱 복잡한 분산 매질 모델링에 용이하지만, 분수 차수 미분으로 인체 분산 FDTD 알고리즘 구현이 어려우며, 구현이 되더라도 계산효율이 매우 떨어지는 단점이 있다.

3-7 2차 복소분수함수 분산 모델링

2차 복소분수함수 분산 모델은 식 (24)로 표현된다^[13].

ϵ r (ω) = A 0 + A 1 (j ω) + A 2 (j ω) 2 1 + B 1 (j ω) + B 2 (j ω) 2

(24)

여기서 A₀, A₁, A₂, B₁, B₂는 계수를 나타낸다. 계수는 복소수 커브피팅법^[28]을 적용하여 행렬 계산으로 추출할 수 있으며, 초기값을 요하지 않는다.

[A 0 A 1 A 2 B 1 B 2] = [λ 0 0 − λ 2 T 1 S 2 0 λ 2 0 − S 2 T 3 λ 2 0 − λ 4 T 3 S 4 T 1 − S 2 − T 3 U 2 0 S 2 T 3 − S 4 0 U 4] = [S 0 T 1 S 2 0 U 2]

(25)

λ h = ∑ k = 0 m (ω k) h

(26)

S h = ∑ k = 0 m (ω k) h R k

(27)

T h = ∑ k = 0 m (− ω k) h I k

(28)

U h = ∑ k = 0 m (ω k) h [(R k) 2 + (I k) 2]

(29)

여기서 ω_k는 샘플링 주파수, k는 인덱스, m은 인덱스의 최대값, R_k는 샘플링 주파수에서 $ε r ′$ 의 측정값, I_k는 샘플링 주파수에서의 $ε r ″$ 측정값을 나타낸다. 한편, 특정 주파수 대역에서 변화율이 큰 복소 유전율은 정확한 분산 특성을 표현하기 어렵다. 따라서 정확한 분산 특성을 표현하기 위해 식 (26)～식 (29)에 주파수에 따른 가중치 함수를 적용할 수 있으며, 최적 가중치(WF: weighting factor)는 PSO(particle swarm optimization) 기법^[29]을 적용하여 도출하고, 제곱평균 오차를 비교하여 정확성을 확인한다.

E r m s = ∑ i N | ε r (ω) − ε r, d a t a (ω i) | 2 ∑ i N | ε r, d a t a (ω i) | 2

(30)

여기서 ε_r,data(ω_i)는 Gabriel 인체조직의 상대 유전율을 나타낸다. 주파수에 따른 인체조직의 상대 유전율은 그림 2와 그림 3에 도시하였으며, 실선은 분산 모델을, 심볼은 Gabriel 측정 데이터를 나타낸다. 그림에서와 같이 측정 데이터와 2차 복소분수함수 모델링이 매우 일치함을 확인할 수 있다.

그림 2. | Fig. 2. 2차 복소분수함수 상대 유전율의 실수부^[13] | Real part of the relative permittivity for 2-pole CRF.

Download Original Figure

그림 3. | Fig. 3. 2차 복소분수함수 상대 유전율의 허수부^[13] | Imaginary part of the relative permittivity for 2-pole CRF.

Download Original Figure

2차 복소분수함수 분산 모델에 대한 업데이트 수식을 도출하기 위해 맥스웰 구성 관계 방정식을 이용한다.

D → (ω) = ϵ 0 [A 0 + A 1 (j ω) + A 2 (j ω) 2 1 + B 1 (j ω) + B 2 (j ω) 2] E → (ω)

(31)

식 (31)에 푸리에 역변환과 중심차분 근사를 적용해 정리하면 전계에 대한 업데이트 수식을 아래와 같이 얻을 수 있다.

E → n + 1 = C a E → n + C b E → n − 1 + C c D → n + 1 + C d D → n + C e D → n − 1

(32)

여기서 α₀=A₀Δt², α₁=2A₁Δt, α₂=4A₂, β₀=Δt²/ϵ₀, β₁=2B₁Δt/ϵ₀, β₂=4B₂/ϵ₀로 정의하면 계수 C_a, C_b, C_c, C_d, C_e는 아래와 같다.

C a = − 2 (α 0 − α 2) / (α 0 + α 1 + α 2) C b = − (α 0 − α 1 + α 2) / (α 0 + α 1 + α 2) C c = (β 0 + β 1 + β 2) / (α 0 + α 1 + α 2) C d = 2 (β 0 − β 2) / (α 0 + α 1 + α 2) C e = (β 0 − β 1 + β 2) / (α 0 + α 1 + α 2)

한편, 식 (33)과 같이 상태 공간(state-space) 신호처리^[30] 기법을 이용한 인체 분산 FDTD 알고리즘은 식 (32)에 비해 20% 개선된 메모리 효율을 가진다.

E → n + 1 = C c D → n + 1 + W → 1 n W → 1 n = C d D → n + C a E → n + W → 2 n W → 2 n = C e D → n − 1 + C b E → n

(33)

그림 4는 주파수에 따른 인체조직의 효과를 나타낸 것이다. 3차원 해석공간은 100×100×100 셀, 공간간격 Δx=Δy=Δz=0.2 mm, 시간간격 Δt=0.3675 ps로 설정하였으며, 인체조직은 20×20×20 셀로 해석영역 중앙에 위치한다. 입력 여기파는 변조된 가우시안 펄스로 S(20, 50, 50)의 위치에서 인가되고, 관찰점은 인체조직의 중심 O(50, 50, 50)로 설정하였다. 그림에서와 같이 2차 복소분수함수 분산 모델 기반 FDTD 알고리즘의 해석 결과가 Yee FDTD 알고리즘의 해석 결과와 매우 일치함을 확인할 수 있다. 한편, 실선은 2차 복소분수함수 분산 FDTD 해석 결과를, 심볼은 Yee FDTD 해석 결과를 나타낸다.

그림 4. | Fig. 4. 인체조직의 효과^[13] | Effects of human tissues.

Download Original Figure

3-8 4차 복소분수함수 분산 모델링

2차 복소분수함수 분산 모델은 D와 E의 구성 관계를 멀티폴(multipole)로 표현할 수 없으므로, 광대역 인체 분산 모델링에는 적합하지 않다^[14],^[28]. 따라서 인체조직의 광대역 분산 특성은 4차 복소분수함수 분산 모델을 이용해 표현할 수 있다. 주파수에 따른 인체조직의 상대 유전율은 그림 5와 그림 6에 도시하였다. 그림에서와 같이 4차 복소분수함수 분산 모델링이 측정 데이터와 매우 일치함을 확인할 수 있다. 한편, 실선은 분산 모델을, 심볼은 Gabriel 측정 데이터를 나타낸다.

ϵ r (ω) = A 0 + A 1 (j ω) + A 2 (j ω) 2 + A 3 (j ω) 3 + A 4 (j ω) 4 1 + B 1 (j ω) + B 2 (j ω) 2 + B 3 (j ω) 3 + B 4 (j ω) 4

(34)

그림 5. | Fig. 5. 4차 복소분수함수 상대유전율의 실수부^[14] | Real part of the relative permittivity for 4-pole CRF.

Download Original Figure

그림 6. | Fig. 6. 4차 복소분수함수 상대유전율의 허수부^[14] | Imaginary part of the relative permittivity for 4-pole CRF^[14].

Download Original Figure

인체 분산 FDTD 알고리즘은 2차 복소분수함수 분산 모델과 유사한 절차로 구현된다. 상태 공간 신호처리 기법을 이용한 전계의 업데이트 식은 식 (35)와 같으며, 메모리 효율은 기존 대비 33% 개선된다.

E → n + 1 = C e D → n + 1 + W → 1 n W → 1 n = C f D → n + C a E → n + W → 2 n W → 2 n = C g D → n − 1 + C b E → n − 1 + W → 3 n W → 3 n = C h D → n − 2 + C b E → n − 2 + W → 4 n W → 4 n = C i D → n − 3 + C d E → n − 3

(35)

여기서 α₀=A₀Δt⁴, α₁=A₁Δt³, α₂=A₂Δt², α₃=A₃Δt, α₄=A₄, β₀=B₀Δt⁴, β₁=B₁Δt³, β₂=B₂Δt², β₃=B₃Δt, β₄=B₄로 정의하면 계수 C_a, C_b, C_c, C_d, C_e, C_f, C_g, C_h, C_i는 아래와 같다.

C a = (α 1 + 2 α 2 − 2 α 3 − 8 α 4) / (α 3 + 2 α 4) C b = − 2 (α 0 − 2 α 2 + 6 α 4) / (α 3 + 2 α 4) C c = − (− α 1 + 2 α 2 + 2 α 3 − 8 α 4) / (α 3 + 2 α 4) C d = − (− α 3 + 2 α 4) / (α 3 + 2 α 4) C e = (β 3 + 2 β 4) / ϵ 0 (α 3 + 2 α 4) C f = (β 1 + 2 β 2 − 2 β 3 − 8 β 4) / ϵ 0 (α 3 + 2 α 4) C g = 2 (β 0 − 2 β 2 + 6 β 4) / ϵ 0 (α 3 + 2 α 4) C h = (− β 1 + 2 β 2 + 2 β 3 − 8 β 4) / ϵ 0 (α 3 + 2 α 4) C i = (− β 3 + 2 β 4) / ϵ 0 (α 3 + 2 α 4) .

그림 7과 그림 8은 인체 분산 FDTD 알고리즘의 정확성을 확인하기 위해 1차원 불균질(inhomogeneous) 구조에서 반사계수와 투과계수를 이론값과 비교한 것이다. 입력 여기파는 양의 방향으로 자유공간과 인체조직의 경계(z=0)에 입사된다. 반사계수는 z=−0.2 mm에서 계산되고, 투과계수는 z=0.2 mm에서 계산된다. 그림에 서와 같이 4차 복소분수함수 분산 모델 기반 FDTD 알고리즘의 해석 결과가 이론값과 매우 일치함을 확인할 수 있다. 여기서 실선은 4차 복소분수함수 분산 FDTD 해석을, 심볼은 이론값을 나타낸다. 지금까지 인체 분산 FDTD 모델들에 대해 검토하였다. 표 2와 표 3은 인체 분산 특성을 표현하는 7개의 모델들을 정리한 것이다.

그림 7. | Fig. 7. 반사계수와 투과계수의 크기^[14] | Amplitudes of reflection (Γ) and transmission (τ) coefficients^[14].

Download Original Figure

그림 8. | Fig. 8. 반사계수와 투과계수의 각도^[14] | Angles of reflection (Γ) and transmission (τ) coefficients^[14].

Download Original Figure

표 2. | Table 2. 인체 분산 모델 수식 | Dispersion model equation for human tissue.

Dispersion model	Equation
Debye^[8]	$ϵ r (ω) = ϵ ∞ + ∑ p = 1 P ϵ s, p − ϵ ∞ 1 + j ω τ p + σ 1 + j ω ϵ 0$
Lorentz^[9]	$ϵ r (ω) = ϵ ∞ + ∑ p = 1 P (ϵ s, p − ϵ ∞) ω p 2 ω p 2 + 2 j ω δ p − ω 2$
Cole-Cole^[10]	$ϵ r (ω) = ϵ ∞ + ∑ p = 1 4 ϵ s, p − ϵ ∞ 1 + (j ω τ p) 1 − α + σ 1 + j ω ϵ 0$
Davidson-Cole^[11]	$ϵ r (ω) = ϵ ∞ + ∑ p = 1 p ϵ s, p − ϵ ∞ (1 + j ω τ p) β$
Havriliak-Negami^[12]	$ϵ r (ω) = ϵ ∞ + ∑ p = 1 p ϵ s, p − ϵ ∞ [(1 + j ω τ p) 1 − α] β$
Two-pole CRF^[13]	$ϵ r (ω) = A 0 + A 1 (j ω) + A 2 (j ω) 2 1 + B 1 (j ω) + B 2 (j ω) 2$
Four-pole CRF^[14]	$ϵ r (ω) = A 0 + A 1 (j ω) … + A 4 (j ω) 4 1 + B 1 (j ω) + … + B 4 (j ω) 4$

Download Excel Table

표 3. | Table 3. 인체 분산 모델 특징 | Dispersion model characteristics for human tissue.

Dispersion model	Characteristics
Debye^[8]	They should use a nonlinear optimization. The accuracy is sensitive to initial values.
Lorentz^[9]
Cole-Cole^[10]	They are suitable for complicated dispersive media. FDTD implementation is difficult due to the fractional derivatives.
Davidson-Cole^[11]
Havriliak-Negami^[12]
Two-pole CRF^[13]	No complicated optimization is needed. No initial values are needed.
Four-pole CRF^[14]

Download Excel Table

한편, 4차 복소분수함수 분산 모델 기반 FDTD 알고 리즘의 수치적 안정도는 맥스웰 방정식과 구성 관계식에 중심차분 근사를 적용하여 안정도 다항식을 유도하였다^[14].

3-9 인체 팬텀 모델링

인체의 경우, 부신(adrenal gland), 동맥(artery), 방광(bla-dder) 등과 같이 78개의 다양한 조직으로 구성되어 있어 각각의 조직에 대한 복합 유전율을 식 (36)과 같이 고려해야 한다.

E x | i, j + 1 / 2, k + 1 / 2 n + 1 / 2 = C a | i, j + 1 / 2, k + 1 / 2 E x | i, j + 1 / 2, k + 1 / 2 n − 1 / 2 + C b y | i, j + 1 / 2, k + 1 / 2 (H z | i, j + 1, k + 1 / 2 n − H z | i, j, k + 1 / 2 n) − C b z | i, j + 1 / 2, k + 1 / 2 (H y | i, j + 1 / 2, k + 1 n − H y (| i, j + 1 / 2, k n))

(36)

여기서 계수 C_a, C_bs는 아래와 같이 정의된다.

C a | i, j + 1 / 2, k + 1 / 2 = (1 − σ i, j + 1 / 2, k + 1 / 2 Δ t 2 ϵ i, j + 1 / 2, k + 1 / 2) / (1 + σ i, j + 1 / 2, k + 1 / 2 Δ t 2 ϵ i, j + 1 / 2, k + 1 / 2) C b s | i, j + 1 / 2, k + 1 / 2 = (σ i, j + 1 / 2, k + 1 / 2 Δ t ϵ i, j + 1 / 2, k + 1 / 2 Δ s) / (1 + σ i, j + 1 / 2, k + 1 / 2 Δ t 2 ϵ i, j + 1 / 2, k + 1 / 2)

공간상의 위치 (i, j, k)에 따라 변경되는 복소 유전율을 표현하기 위해 FDTD 알고리즘 계수 C_a, C_bs, ε, σ는 3차원 실수 배열로 선언되기 때문에 해석 시 많은 메모리 공간을 요한다. 따라서 인체 분산 FDTD 알고리즘 구현 시 다양한 인체조직의 복소 유전율을 효율적으로 표현하기 위해 미디어 기법(Media Method)^[2]를 이용한다. 미디어 기법은 조직에 따라 공간을 나누고, 계수를 미디어 번호로 설정하는 기법이다. 미디어 기법을 적용한 전계에 대한 업데이트 식은 식 (37)과 같다.

m = M E D I A i, j + 1 / 2, k + 1 / 2 E x | i, j + 1 / 2, k + 1 / 2 n + 1 / 2 = C a | m E x | i, j + 1 / 2, k + 1 / 2 n − 1 / 2 + C b y | m (H z | i, j + 1 / 2, k + 1 / 2 n − H z | i, j, k + 1 / 2 n) − C b z | m (H y | i, j + 1 / 2, k + 1 n − H y | i, j + 1 / 2, k n)

(37)

인체조직의 복소 유전율에 대한 정보는 Gabriel의 측정 데이터^[22]를 이용하고, 인체구조에 대한 정보는 Virtual Family Tool^[31]에서 복셀(voxel) 데이터를 추출해 이용한다.

Ⅳ. 흡수 경계 조건

일반적인 문제에서 전자기파는 무한하게 퍼져 나가기 때문에 무한 공간에 대한 계산이 요구된다. 그러나 수치적으로 무한 공간을 계산하는 것은 불가능하므로 해석 공간을 제한하는 흡수 경계 조건(ABCs: absorbing boundary conditions)이 필요하다^[2]. 흡수 경계 조건은 1980년대 이후에 많은 종류가 개발되었으며, 크게 해석적인 방법과 가상의 흡수 매질을 이용하는 방법이 있다. Mur^[32], Liao^[33], Higdon^[34],^[35]에 의해 제안된 해석적인 방법은 파동 방정식에서 최외각 경계면의 전자계의 성분을 이용해 최외각 전계와 자계의 값을 구하고, 이를 통해 전파하는 것처럼 구현한다. 한편, 가상의 흡수 매질을 이용하는 방법은 해석 공간의 외각 영역에 전자기파 흡수 매질을 둠으로써 해석 공간으로 반사파가 발생하지 않게 하는 경계 조건이다. 1994년 Berenger에 의해 제안된 PML(Perfectly Matched Layer)^[36]은 하나의 필드 성분을 두 개의 수직 성분으로 나누어 표현하는 Split-Field로 구현해 주파수와 입사각에 무관하게 이론적으로 완전한 흡수 성능을 가진다. 따라서 PML은 분산 특성이 강한 복합 매질 해석에 매우 유용하게 사용된다^[37]. 한편, Kuzuoglu와 Mittra는 기존 PML의 단점을 극복하는 CFS(Complex Frequency Shifted) PML을 제안하였다^[38]. CFS-PML은 기존 PML과 달리, 저주파수 영역 및 감쇄파에 대해 좋은 흡수성능을 나타내며, 매우 긴 해석이 요구되는 문제에 대해 late-time 수치발산이 일어나지 않는 장점이 있다.

Ⅴ. 결 론

본 논문에서는 인체와 전자기파의 상호작용을 해석하기 위한 인체 분산 FDTD 모델링에 대해 검토하였다. 이를 위해, 먼저 인체조직의 분산 특성을 표현하는 Debye 모델, Lorentz 모델, Cole-Cole 모델, Davidson-Cole 모델, Havriliak-Negami 모델, 2차 복소분수함수 모델, 4차 복소분수함수 모델들에 대해 논의하고, 분산 모델들을 FDTD 알고리즘에 적용하는 방법에 대해 검토하였다. 또한, 분산 매질을 FDTD 알고리즘으로 해석하는 PLRC, RC, ADE 방법에 대해 검토하였다. 마지막으로, 인체의 기하학적 구조를 정확하고, 효율적으로 모델링하기 위한 방법으로 인체팬텀과 미디어 기법에 대해 논의하고, 분산 특성이 강한 복합 매질 해석에 용이한 흡수 경계 조건에 대해 검토하였다.

References

[1].

K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302-307, May 1966.

[2].

A. Taflove, S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. Boston, MA, Artech House, 2005.

[3].

S. D. Gedney, Introduction to the Finite-Difference Time-Domain(FDTD) Method for Electromagnetics, London, Morgan & Claypool, 2011.

[4].

R. F. Harrington, Field Computation by Moment Methods, New York, NY, Maclillan, 1968.

[5].

J. Jin, The Finite Element Method in Electromagnetics, New York, NY, John Wiley & Sons, 1993.

[6].

K. Y. Jung, F. L. Teixeira, and R. M. Rean, “Au/SiO₂ nanoring plasmon waveguides at optical communication band,” Journal of Lightwave Technology, vol. 25, no. 9, pp. 2757-2765, Sep. 2007.

[7].

M. A. Sefunc, A. K. Okyay, and H. V. Demir, “Plasmonic backcontact grating for P3HT:PCBM organic solar sells enabling strong optical absorption increased in all polarizations,” Optics Express, vol. 19, no. 15, pp. 14200-14209, Jul. 2011.

[8].

H. K. Rouf, F. Costen, and M. Fujii, “Modelling EM wave interaction with human body in frequency dependent Crank Nicolson method,” Journal of Electromagnetic Waves and Applications, vol. 25, no. 17-18, pp. 2429-2441, 2011.

[9].

M. A. Alsunaidi, A. A. Al-Jabr, “A general ADE-FDTD algorithm for the simulation of dispersive structures,” IEEE Photonics Technology Letters, vol. 21, no. 12, pp. 817-819, Jun. 2009.

[10].

I. T. Rekanos, T. G. Papadopoulos, “An auxiliary differential equation method for FDTD modeling of wave propagation in Cole-Cole dispersive media,” IEEE Transactions on Antennas and Propagation, vol. 58, no. 11, pp. 3666-3674, Nov. 2010.

[11].

I. T. Rekanos, “FDTD schemes for wave propagation in Davidson-Cole dispersive media using auxiliary differential equations,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 3, pp. 1467-1478, Mar. 2012.

[12].

I. T. Rekanos, “FDTD modeling of Havriliak-Negami media,” IEEE Microwave and Wireless Components Letters, vol. 22, no. 2, pp. 49-51, Feb. 2012.

[13].

S. G. Ha, J. Cho, J. Choi, H. Kim, and K. Y. Jung, “FDTD dispersive modeling of human tissues based on quadratic rational function,” IEEE Transactions on Antennas and Propagation, vol. 61, no. 2, pp. 996-999, Feb. 2013.

[14].

S. G. Ha, J. Cho, E. K. Kim, Y. B. Park, and K. Y. Jung, “FDTD dispersive modeling with high-order rational constitutive parameters,” IEEE Transactions on Antennas and Propagation, vol. 63, no. 9, pp. 4233-4238, Sep. 2015.

[15].

D. F. Kelley, R. J. Luebbers, “Piecewise linear recursive convolution for dispersive media using FDTD,” IEEE Transactions on Antennas and Propagation, vol. 44, no. 6, pp. 792-797, Jun. 1996.

[16].

G. X. Fan, Q. H. Liu, “An FDTD algorithm with perfectly matched layers for general dispersive media,” IEEE Transactions on Antennas and Propagation, vol. 48, no. 5, pp. 637-646, May 2000.

[17].

R. Luebbers, F. R. Hunsberger, K. S. Kunz, R. B. Standler, and M. Schneider, “A frequency-dependent finite difference time-domain formulation for dispersive materials,” IEEE Transactions on Electromagnetic Compatibility, vol. 32, no. 3, pp. 222-227, Aug. 1990.

[18].

R. Luebbers, F. Hunsberger, and K. Kunz, “A frequency-dependent finite difference time-domain formulation for transient propagation in plasma,” IEEE Transactions on Antennas and Propagation, vol. 39, no. 1, pp. 29-34, Jan. 1991.

[19].

R. Luebbers, “Lossy dielectric in FDTD,” IEEE Transactions on Antennas and Propagation, vol. 41, no. 11, pp. 1586-1588, Nov. 1993.

[20].

T. Kashiwa, I. Fukai, “A treatment by FD-TD method of dispersive characteristics associated with electronic polarization,” Microwave and Optics Technology Letters, vol. 3, no. 6, pp. 203-205, Jun. 1990.

[21].

R. M. Joseph, S. C. Hagness, and A. Taflove, “Direct time integration of Maxwell's equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses,” Optics Letters, vol. 16, no. 18, pp. 1412-1414, Sep. 1991.

[22].

C. Gabriel, “Compilation of the dielectric properties of body tissue at RF and microwave frequency,” Brooks Air Force, San Antonio, TX, AL/OE-TR-1996-0037, 1996.

[23].

J. Wang, D. Su, “Design of an ultra wideband system for in-body wireless communications,” in the 2006 4th Asia-Pacific Conference on Environmental Electromagnetics, Dalian, 2006, pp. 565-568.

[24].

O. P. Gandhi, C. M. Furse, “Current induced in the human body for exposure to ultrawideband electromagnetic pulses,” IEEE Transactions on Electromagnetic Compatibility, vol. 39, no. 2, pp. 174-180, May 1997.

[25].

T. Wuren, T. Takai, M. Fujii, and I. Sakagami, “Effective 2-Debye-pole FDTD model of electromagnetic interaction between whole human body and UWB radiation,” IEEE Microwave and Wireless Components Letters, vol. 17, no. 7, pp. 483-485, Jul. 2007.

[26].

M. A. Alsunaidi, A. A. Al-Jabr, “General ADE-FDTD algorithm for the simulation of dispersive structures,” IEEE Photonics Technology Letters, vol. 21, no. 12, pp. 817-819, Jun. 2009.

[27].

K. B. Oldham, J. Spanier, The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order, Doven, NY, Mineola, 2006.

[28].

E. C. Levi, “Complex-curve fitting,” IRE Transactions on Automatic Control, vol. 4, no. 1, pp. 37-44, May 1959.

[29].

H. Chung, S. G. Ha, J. Choi, and K. Y. Jung, “Accurate FDTD modeling for dispersive media using rational function and particle swarm optimization,” International Journal of Electronics, vol. 102, no. 7, pp. 1218-1228, Oct. 2014.

[30].

B. C. Kuo, Automatic Control Systems, Englewood Cliffs, NJ, Prentice-Hall, 1975.

[31].

A. Christ, W. Kainz, E. G. Hahn, K. Honegger, M. Zefferer, and E. Neufeld, et al., “The virtual family-development of surface-based anatomical CAD models of two adults and two children for dosimetric simulations,” Physics in Medicine and Biology, vol. 55, no. 2, pp. 23-38, Dec. 2010.

[32].

G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Transactions on Electromagnetic Compatibility, vol. EMC-23, pp. 377-382, Nov. 1981.

[33].

Z. F. Liao, K. L. Huang, B. P. Yang, and Y. F. Yuan, “A transmitting boundary for transient wave analysis,” Science in China Series A-Mathematics, Physics, Astronomy & Technological Science, vol. 27, no. 10, pp. 1063-1076, Oct. 1984.

[34].

R. L. Higdon, “Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation,” Mathematics of Computation, vol. 47, no. 176, pp. 437-459, Oct. 1986.

[35].

R. L. Higdon, “Numerical absorbing boundary conditions for the wave equation,” Mathematics of Computation, vol. 49, pp. 65-90, Jul. 1987.

[36].

J. P. Berenger, “A perfectly matched layer for the absorption of electromagntic waves,” Journal of Computational Physics, vol 144, no. 2, pp. 185-200, Oct. 1994.

[37].

K. Y. Jung, S. Ju, and F. L. Teixeira, “Two-stage perfectly matched layer for the analysis of plasmonic structures,” IEICE Transactions on Electronics, vol. E93-C, no. 8, pp. 1371-1374, Aug. 2010.

[38].

M. Kuzuoglu, R. Mittra, “Frequency dependence of the constitutive parameters of causal perfectly matched absorbers,” IEEE Microwave and Guided wave letters, vol. 6, no. 12, pp. 447-449, Dec. 1996.

Author Information

하 상 규[다쏘시스템/솔루션컨설턴트]

jkiees-31-3-205-i1

https://orcid.org/0000-0002-3826-1025
2011 2월: 아주대학교 전자공학과 (공학사)
2017년 8월 한양대학교 전자컴퓨터통신공학과 (공학박사)
2017년 7월～2019년 8월: LG전자 선임연구원
2019년 9월～현재: 다쏘시스템코리아 솔루션 컨설턴트
[주 관심분야] 전자파 수치해석, 인체 전자파 모델링, 비포스터 안테나, 5G 위상 배열 안테나

정 경 영[한양대학교/부교수]

jkiees-31-3-205-i2

https://orcid.org/0000-0002-7960-3650
1996년 2월: 한양대학교 전파공학과 (공학사)
1998년 2월: 한양대학교 전자통신공학과 (공학석사)
2008년 8월: 미국 Ohio State University 전기컴퓨터공학 (공학박사)
1998년 1월~2001년 4월: 현대전자 전임연구원
2001년 5월～2004년 5월: 팬택앤큐리텔 선임연구원
2008년 8월～2009년 2월: 미국 Ohio State University 전기컴퓨터공학 Post-Doctoral Researcher
2009년 3월～2011년 2월: 아주대학교 전자공학부 전임강사
2011년 3월～2016년 2월: 한양대학교 융합전자공학부 조교수
2016년 3월～현재: 한양대학교 융합전자공학부 부교수
[주 관심분야] 전자파 수치해석, 나노 전자파, 인체 전자파