Discrepancies of Measured SAR between Traditional and Fast Measuring Systems

Abstract

Human exposure to mobile devices is traditionally measured by a system in which the human body (or head) is modelled by a phantom and the energy absorbed from the device is estimated based on the electric fields measured with a single probe. Such a system suffers from low efficiency due to repeated volumetric scanning within the phantom needed to capture the absorbed energy throughout the volume. To speed up the measurement, fast SAR (specific absorption rate) measuring systems have been developed. However, discrepancies of measured results are observed between traditional and fast measuring systems. In this paper, the discrepancies in terms of post-processing procedures after the measurement of electric field (or its amplitude) are investigated. Here, the concerned fast measuring system estimates SAR based on the reconstructed field of the region of interest while the amplitude and phase of the electric field are measured on a single plane with a probe array. The numerical results presented indicate that the fast SAR measuring system has the potential to yield more accurate estimations than the traditional system, but no conclusion can be made on which kind of system is superior without knowledge of the field-reconstruction algorithms and the emitting source.

1. Introduction

Human exposure [1] has drawn much public attention recently due to the wide usage of wireless communication equipment. Ensuring the absorbed energy lies within the safe range [2] requires the accurate quantification of the so-called specific absorption rate (SAR). Efforts towards such quantification have been made by researchers both in simulation tools [3,4,5] and measuring systems [6,7,8,9]. For simulation techniques, when the distance is large (relative to the wavelength of the electric field), the front of incident wave upon the body can be treated as planar [10,11], while near fields need to be analyzed with a high-resolution human module with a numerical code (FDTD [12] or FEM [13]) when the human body is close to sources [14,15,16]. In measurements, the exposure to base stations is usually estimated by measuring the intensity of the electric field at various locations and through the construction of a path loss model it is possible to have a statistical view [17]. The exposure to user equipment is traditionally measured by moving a probe inside a liquid-filled phantom, which simulates the composition of human body or head, to estimate the whole-body and/or local SAR [9].

The traditional measuring system [9] when applied to user equipment is the concern of this paper. The peak spatial-average ( 1 $\mathrm{g}$ or 10 $\mathrm{g}$) SAR is often of interest in practice. The probe carries out the so-called area scan (two-dimensional scanning over a coarse grid to find the location of maximum SAR) and zoom scan (three-dimensional scanning with a finer grid to determine the peak mass-averaged SAR using interpolation and extrapolation techniques). A complete measurement consumes tens of minutes while more than 100 measurements have to be performed to check full compliance (with different frequency bands, working modes and device positions) of a product. Considering that with today’s massive production the compliance of millions (even billions) of products needs to be assessed, the efficiency of such a measuring system is intolerable and in great demand to be improved.

Research has been carried out to speed up the process by reducing the number of measurement points and deducing the peak spatial-average SAR based on parametric models [18] or empirical observations [19]. However, since the method is model-dependent, the estimation accuracy for newly emerged devices (e.g., equipped with MIMO terminals) is not guaranteed. The estimation approach based on the technique of plane-wave expansion [20] has been proposed in which only the electromagnetic properties of the medium are required. Such an approach has recently been applied to the computation of power density of millimetre waves [8].

As described above, various approaches exist to obtain improvements in SAR measurement. The problem arises when discrepant results are observed from fast measuring systems that are developed by different manufacturers [21]. Due to reasons of commercial security, it is hard to justify which system is more accurate than others and also not easy to conclude that the fast measuring system is biased when the generated results challenge the traditional system, which is usually considered as the reference.

Efforts are made in this paper to reveal the reason for estimation discrepancies and try to answer the following questions:

• Why do discrepancies appear for the estimation of SAR by different fast measuring systems?
• Can we say fast measuring systems generate biased estimations if they differ appreciably from the traditional SAR measuring system?
• Which of the traditional measuring system and the fast measuring system is the more accurate?

Despite various methodologies for the fast SAR measuring system, the concerned system measures the electric field inside the phantom by a vector-probe on a plane. Then, the field in other positions of interest is reconstructed by algorithms (e.g., plane-wave expansion), which can be carried out with high computational efficiency. The study is based on analytical functions, which simulate the wave propagation of electric field inside the phantom. To avoid too complex an analysis, the estimation discrepancies due to the post-processing procedures are mainly investigated and factors like scattering from the phantom shell will not be considered.

The remaining part of this paper is organized as follows. Section 2 and Section 3 introduce the methodology of SAR estimation in the traditional system and the concerned fast measuring system, respectively. Numerical results with the flat phantom and 11 emitting sources in Section 4 indicate the accuracy and stability of the fast estimation and present the estimation discrepancies between two kinds of systems. Conclusions are made in Section 5.

2. Traditional SAR Measuring System

The definition of the spatial-average specific absorption rate (sSAR) is

$$\mathrm{sSAR}=\frac{1}{V}{\int }_{\mathbb{V}}\frac{{\left|\mathbf{E}\left(\mathbf{r}\right)\right|}^2\sigma \left(\mathbf{r}\right)}{\rho \left(\mathbf{r}\right)}\mathrm{d}\mathbf{r},$$

(1)

where V is the volume of the region of interest $\mathbb{V}$, $\left|\mathbf{E}\right(\mathbf{r}\left)\right|$ is the root-mean-square of the electric field at the location $\mathbf{r}$, and $\sigma$ and $\rho$ denote conductivity and density, respectively. $\mathbb{V}$ is a cube with side length 1 mm for 1 $\mathrm{g}$ SAR and $21.5$ mm for 10 $\mathrm{g}$ SAR when $\rho =1\mathrm{g}{/}^3$. The value of $\sigma$ depends on the wave frequency and the reference value can be found in the literature [22,23]. As seen from (1), the phase of the electric field is not required and thus only the amplitude is measured by the probe in the traditional measuring system, the composition of which is sketched in Figure 1, where a probe is moved by a robot arm to measure the amplitude of electric fields inside a phantom due to the emitting device under test. The integration in (1) is carried out numerically (e.g., by the trapezoidal rule) and the amplitude of the electric field is required at dense sampling points inside $\mathbb{V}$.

Rather than making intensive measurements inside the whole phantom, area scan and zoom scan are performed sequentially and followed by interpolations and extrapolations to provide the amplitude of the electric fields at the desired locations. The procedure is summarized below.

• Area scan: measure fields according to a two-dimensional coarse grid, the distance of which to the phantom surface is fixed, to locate the local maxima of the amplitude of electric fields.
• Zoom scan: a three-dimensional scanning within cubes centered at the location of local maxima, the grid step being smaller than that in the area scan.
• Interpolation and extrapolation: linear interpolation and cubic spline interpolation (and extrapolation) are used as necessary to deduce the amplitude at the points in a finer grid.
• Peak spatial-average SAR: obtained by performing numerically the integration in (1) based on the interpolated and extrapolated amplitude.

Documentary standards are available for the specific requirements on the above measurements and post-processing. Here, the standard [9] is followed, which includes the requirements in

Table 1. Configurations of traditional sSAR measuring system.

Area scan	maximum grid spacing	20 mm if $f<$ 3 GHz and $60/f$ mm otherwise
	maximum distance between probe and surface of phantom	5 mm if $f<$ 3 GHz and $\delta ln2/2$ mm otherwise
Zoom scan	horizontal grid spacing	$\le min\{24/f,8\}$ mm
	minimum scan size	$30\mathrm{mm}\times 30\mathrm{mm}\times 30\mathrm{mm}$ if $f<$ 3 GHz and $22\mathrm{mm}\times 22\mathrm{mm}\times 22\mathrm{mm}$ otherwise
	maximum distance between probe and surface of phantom	5 mm if $f<$ 3 Ghz and $\delta ln2/2$ mm otherwise

“ln” denotes natural logarithm, f wave frequency in GHz, $\delta$ plane-wave skin depth.

. For the interpolation and extrapolation, no specific algorithms are required or recommended. Here, the method of linear and cubic spline interpolation are applied. Note that the sSAR drift is also tested for practical instruments, but not considered during the simulation.

3. Fast SAR Measuring System Based on Field Reconstruction

Rather than a single probe, a probe array is used in the fast SAR measuring system, which is sketched in Figure 2. The probes measure the amplitude and phase of electric fields on a single plane. The solution to fields at other positions of interest is obtained by a field-reconstruction algorithm. Note that the algorithms used in commercial products are usually inaccessible for reasons of security protection. Here, the technique of plane-wave expansion (PWE), which is popularly utilized when measurements are taken on plane(s), is applied. To study the effects of different field-reconstruction algorithms, the PWE method, which is found inaccurate for some cases, is performed with different settings.

3.1. Plane-Wave Expansion (PWE)

Based on the PWE theory, the electric field can be represented as an integral of planar waves,

$$E(x,y,z)=\frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\tilde{E}(k_x,k_y)e^{i(k_{x}x+k_{y}y-k_{z}z)}\mathrm{d}{k}_x\mathrm{d}{k}_y,$$

(2)

where $(x,y,z)$ is the Cartesian coordinates of the observation point, $k_x$, $k_y$ and $k_z$ are respectively the x, y, z component of the wavenumber of the liquid k with the identity $k_x^2+k_y^2+k_z^2=k^2$. The imaginary part of $k_z$ is non-negative so that the Sommerfeld radiation condition [24] is satisfied. Note that (2) is valid for the x, y and z component of the electric fields and the three components can be independently treated. Thus, the following derivations apply to each component.

From (2), we see that $E(x,y,z)$ is a two-dimensional inverse Fourier transform of $\tilde{E}(k_x,k_y)e^{-i{k}_{z}z}$, i.e., $E(x,y,z)$ can be reconstructed if the corresponding spectrum $\tilde{E}(k_x,k_y)e^{-i{k}_{z}z}$ is known. Set $\tilde{E}(k_x,k_y,z)$ $=\tilde{E}(k_x,k_y)e^{-i{k}_{z}z}$. It follows that

$$E(x,y,z_{\mathrm{rec}})=\frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\tilde{E}(k_x,k_y,z_{\mathrm{mea}})e^{-i{k}_z(z_{\mathrm{rec}}-z_{\mathrm{mea}})}{e}^{i(k_{x}x+k_{y}y)}\mathrm{d}{k}_x\mathrm{d}{k}_y,$$

(3)

where $z_{\mathrm{mea}}$ and $z_{\mathrm{rec}}$ denote the z-coordinate of the measurement and reconstruction plane, respectively, and the multiplicand $e^{-i{k}_z(z_{\mathrm{rec}}-z_{\mathrm{mea}})}$ is the so called spectral propagator. From (2), we see $\tilde{E}(k_x,k_y,z_{\mathrm{mea}})$ is the two-dimensional Fourier transform of $E(x,y,z_{\mathrm{mea}})$, i.e.,

$$\tilde{E}(k_x,k_y,z_{\mathrm{mea}})={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }E(x,y,z_{\mathrm{mea}})e^{-i(k_{x}x+k_{y}y)}\mathrm{d}x\mathrm{d}y,$$

(4)

which is to be estimated based on the measured fields.

Electric fields are measured on a plane with probes uniformly distributed in a specific domain, whose centre is usually within the main lobe of the radiated field. Denote by $∆x$ and $∆y$ the minimum interval between probes along the x and y axis, respectively. The coordinates of the probes are $(x_p^{n_x},y_p^{n_y})=(x_p^{\mathrm{c}},y_p^{\mathrm{c}})+n_x\widehat{x}+n_y\widehat{y}$, where $(x_p^{\mathrm{c}},y_p^{\mathrm{c}})$ is the centre of the measurement domain, $n_x$, $n_y$ are integers in the intervals $[-(N_x-1)/2,(N_x-1)/2]$ and $[-(N_y-1)/2,(N_y-1)/2]$, respectively. $N_x{N}_y$ is the number of probes and assumed odd. $\widehat{x}$ and $\widehat{y}$ are unit vectors along the x and y axis, respectively.

Assume that the amplitude of electric fields outside the measurement domain is quite small. The sampling density satisfies the Nyquist sampling theorem, i.e., $2\pi /∆x\ge k_x^{\max }$, $2\pi /∆y\ge k_y^{\max }$. $k_x^{\max }$ and $k_y^{\max }$ denote the maximum value of $k_x$ and $k_y$ respectively with the indication that the spectrum $\tilde{E}(k_x,k_y)$ is approximately zero outside the intervals $-k_x^{\max }/2\le k_x\le k_x^{\max }/2$, $-k_y^{\max }/2\le k_y\le k_y^{\max }/2$. Then (4) can be well approximated by the discrete form

$$\tilde{E}(k_x^{m_x},k_y^{m_y},z_{\mathrm{mea}})=\sum _{n_x=-(N_x-1)/2}^{(N_x-1)/2}\sum _{n_y=-(N_y-1)/2}^{(N_y-1)/2}E(x^{n_x},y^{n_y},z_{\mathrm{mea}})e^{-i(k_x^{m_x}{x}^{n_x}+k_y^{m_y}{y}^{n_y})},$$

(5)

where $k_x^{m_x}=m_x{k}_x^{\max }/M_x$, $m_x=-(M_x-1)/2,\dots ,(M_x-1)/2$ and $k_y^{m_y}=m_y{k}_x^{\max }/M_y$, $m_y=-(M_y-1)/2,\dots ,(M_y-1)/2$, $M_x$, $M_y$ being the number of sampled spatial-frequency points, which are set as odd numbers. Setting the centre $(x_p^c,y_p^c)$ of the measurement domain as the origin of the x and y axis, i.e., $x_p^c=0$, $y_p^c=0$, (5) is rewritten as the standard form of the discrete Fourier transform (DFT),

$$E(m_x,m_y,z_{\mathrm{mea}})=\sum _{n_x^{\prime }=0}^{N_x-1}\sum _{n_y^{\prime }=0}^{N_y-1}E(n_x,n_y,z_{\mathrm{mea}})e^{-i2\pi (m_x{n}_x^{\prime }/M_x+m_y{n}_y^{\prime }/M_y)},$$

(6)

with the notation $n_x^{\prime }=n_x+(N_x-1)/2$, $n_y^{\prime }=n_y+(N_y-1)/2$.

Then $\tilde{E}(k_x,k_y,z_{\mathrm{rec}})$ is computed as $\tilde{E}(k_x,k_y,z_{\mathrm{mea}})e^{-i{k}_z(z_{\mathrm{rec}}-z_{\mathrm{mea}})}$ and the field solution in the spatial domain is obtained after the inverse DFT, i.e.,

$$E(n_x,n_y,z_{\mathrm{rec}})=\frac{1}{M_x{M}_y}\sum _{m_x=0}^{M_x-1}\sum _{m_y=0}^{M_y-1}E(m_x,m_y,z_{\mathrm{rec}})e^{i2\pi (m_x{n}_x^{\prime }/M_x+m_y{n}_y^{\prime }/M_y)}.$$

(7)

Set $\mathbf{F}=\left[e^{-i2\pi (m_x{n}_x^{\prime }/M_x+m_y{n}_y^{\prime }/M_y)}\right]$, $\mathbf{invF}=[e^{i2\pi (m_x{n}_x^{\prime }/M_x+m_y{n}_y^{\prime }/M_y)}/\left(M_x{M}_y\right)]$, and $\mathbf{P}(z_{\mathrm{rec}}-z_{\mathrm{mea}})=\mathrm{diag}\left\{e^{-i{k}_z(z_{\mathrm{rec}}-z_{\mathrm{mea}})}\right\}$. The above reconstruction procedures are expressed in matrix form as

$${\mathbf{E}}_{\mathrm{rec}}=\mathbf{invF}·\mathbf{P}·\mathbf{F}·{\mathbf{E}}_{\mathrm{mea}},$$

(8)

where ${\mathbf{E}}_{\mathrm{rec}}$, ${\mathbf{E}}_{\mathrm{mea}}$ are column vectors composed of $N_x{N}_y$ sampled electric fields. In summary, the solution to the field on the reconstruction plane is an inverse DFT (denoted by $\mathbf{invF}$) of the spectrum, which is obtained by multiplying the spectrum of measured fields ${\mathbf{E}}_{\mathrm{mea}}$ by the propagator $\mathbf{F}$.

The PWE approach suffers from being an ill-conditioned problem [25] when reconstructing high spatial-frequency (mentioned as frequency in the remaining part) components. When $k_z$ has a large imaginary part, i.e., $k_x^2+k_y^2\gg k^2$, the value of the propagator $e^{-i{k}_z(z_{\mathrm{rec}}-z_{\mathrm{mea}})}$ would be very large when $z_{\mathrm{rec}}>z_{\mathrm{mea}}$. As a result, the effects of approximation errors (e.g., due to invalid assumptions for (5)) or measurement noise are amplified by the propagator and the estimation would be highly biased and with a large variance.

Considering that the energy of high-frequency components is usually small, a stable reconstruction is often reached, without losing much accuracy, by only considering the spectrum at low frequencies. With lossless cases, only propagating plane-wave components, i.e., $k_x^2+k_y^2<k^2$, are considered. In the concerned cases, since the equivalent liquid is lossy (wavenumber k is complex), the integrand stands for evanescent waves. However, a similar constraint $k_x^2+k_y^2\le {\left|k\right|}^2$ is applied.

Note that in practice, only the x and y component of the electric field, denoted by $E_x$ and $E_y$ respectively, are required to be measured or reconstructed, since the z component can be obtained based on the identity

$${\tilde{E}}_x(k_x,k_y)k_x+{\tilde{E}}_y(k_x,k_y)k_y+{\tilde{E}}_z(k_x,k_y)k_z=0.$$

(9)

3.2. Field Reconstruction Making Use of More High-Frequency Components

When the energy of high-frequency components is small and can be neglected, reconstructing only the low-frequency spectrum is reasonable. Otherwise, since a part of the energy is not taken into account in the algorithm, the SAR value tends to be underestimated, as shown by the numerical results in Section 4. Thus, high-frequency components need to be reconstructed. Replacing the constraint of the spectrum by $\sqrt{k_x^2+k_y^2}<\epsilon$, a higher-frequency spectrum would be considered as $\epsilon$ increases. For the commonly utilized approach in Section 3.1, the threshold $\epsilon$ is set as $\left|k\right|$. To quantify the additional spectrum considered, $\epsilon$ is expressed by $\left|k\right|+(\sqrt{{\left(k_x^{max}\right)}^2+{\left(k_y^{max}\right)}^2-\left|k\right|)}\delta$. The parameter $\delta \in [0,1]$ quantifies how much high-frequency spectrum is reconstructed. For the traditional PWE method, $\delta =0$. The reconstruction of the complete spectrum is achieved by setting $\delta =1$.

4. Numerical Results

The hyperparameters in the described PWE method include the number of frequency points $M_x$, $M_y$ and $\delta$. $M_x$ and $M_y$ are set as $2N_x+1$ and $2N_y+1$, respectively, while no significant effect is observed with larger values. Different values are given to $\delta$ to study the effects of reconstructing the high-frequency spectrum. Uncertainty due to, e.g., measurement errors and coupling effects is not considered in the following analysis except in Section 4.4.

4.1. Configurations

The studied fast measuring system is sketched in Figure 2, where the device under test is placed 5 mm away from the flat phantom and the vector-probe array at the plane with $z=19.25$. An array of $29\times 29$ probes is uniformly distributed in the surface defined by $-10 \mathrm{cm} \le x,y \le 10 \mathrm{cm}$ with interval 7 mm. For analysis purposes, probes are represented as points and measurements of electric fields are assumed to be accurate. The measured electric field inside the phantom is obtained based on the analytical function

$$E\left(\mathbf{r}\right)=\sum _{i=1}^{80}\sum _{j=1}^{80}\left({\nabla }^2{G}_k(\mathbf{r}-{\mathbf{d}}_{i,j}){\mathbf{p}}_{i,j}+k^2{G}_k(\mathbf{r}-{\mathbf{d}}_{i,j}){\mathbf{p}}_{i,j}\right).$$

(10)

The scatterings by different emitting sources are simulated by varying the setting of dipole position ${\mathbf{d}}_{i,j}\in {\mathbb{R}}^3$ and moment ${\mathbf{p}}_{i,j}\in {\mathbb{R}}^3$. $\mathbf{r}$ denotes the position of the observation point, $G_k\left(\mathbf{r}\right)=e^{ik\left|\mathbf{r}\right|}/4\pi \left|\mathbf{r}\right|$ and k is the wavenumber of the phantom. The 11 cases tested are with the flat phantom and configured according to

Table 2. Physical parameters and reference values of peak 10 $\mathrm{g}$ SAR.

Index	1	2	3	4	5	6	7	8	9	10	11
f (MHz)	850	1800	1900	2450	5500	5800	750	1950	750	835	1750
${ϵ}_r$	42.23	40.45	40.28	39.37	33.30	32.64	42.47	40.20	42.47	42.26	40.53
$\sigma$ (S/m)	0.89	1.39	1.45	1.87	5.18	5.55	0.85	1.49	0.85	0.88	1.35
10g SAR	0.58	0.48	0.48	0.43	0.29	0.28	0.28	0.41	0.66	0.65	0.52

. Note that the reference value of peak 1 $\mathrm{g}$ spatial-average SAR (mentioned as 1 $\mathrm{g}$ SAR in the following part) equals 1 for all cases.

4.2. Verification of Post-Processing Procedures

The post-processing includes the field reconstruction of regions of interest and the followed computation of peak spatial-average SAR. The former is verified by observing the reconstructed field on the plane of $z=0$ mm and comparing with the reference as shown in Figure 3 for the fourth case, where the maps of field intensity in the first row are the reference, the second row gives the reconstructed results, and the last row shows the absolute error of reconstruction. The resolution of the map is 2 mm. As seen, the amplitude of the reconstructed fields is a little different from the reference, the absolute errors being small. Setting the resolution as 1 mm, with the reconstructed field of the cube bounded by $-10 \mathrm{cm} \le x,y \le 10 \mathrm{cm}$, $0 \mathrm{cm} \le z \le 3 \mathrm{cm}$, the computed peak spatial-average SAR equals $0.993$ for 1 $\mathrm{g}$ and $0.427$ for 10 $\mathrm{g}$, which are quite close to the reference values 1 and $0.431$.

The computed peak 1 $\mathrm{g}$ and 10 $\mathrm{g}$ SAR of the 11 cases based on the reconstructed field are shown in Figure 4. For most cases, with the commonly used PWE approach, i.e., setting $\delta =0$, the peak sSAR value can be accurately estimated. However, e.g., with respect to the first, seventh, and eighth case, increasing the value of $\delta$ significantly improves the estimation accuracy since a larger part of high spatial-frequency spectrum is considered. Note that the value of $\delta$ should not be too large due to the involved ill-conditioned problem. For instance, when $\delta$ is increased to 0.03, the estimated peak sSAR value changes from closely approximated solutions to a highly biased one. An insight into the ill-conditioned problem is given through the investigation of the field reconstruction problem with respect to the seventh case and the uncertainty analysis.

4.3. Problem in Field Reconstructions

While the development of advanced field-reconstruction algorithms is not the concern of this paper, the challenges encountered in the field reconstruction can be presented. With respect to the seventh case, where the value of peak sSAR is underestimated, the amplitude of the reconstructed $x,y$ component of electric fields is shown in Figure 5a. The colour limits are set to be the same for each component. Note that since the z component is determined by the $x,y$ counterparts, only the x and y components of the electric field are shown. As observed, the electric fields with a large amplitude are distributed intensively in the spatial domain. Consequently, as shown by Figure 5b, the energy of the spectrum spreads widely and the high-frequency spectrum cannot be neglected. If following the commonly utilized approach and only reconstructing the spectrum constrained by $k_x^2+k_y^2\le {\left|k\right|}^2$, the smooth region of the electric field is well reconstructed but the small regions, which are centered at the peak values and influential to the estimation of the peak spatial-average SAR, are poorly reconstructed. Increasing $\delta$ by 0.03 to make more use of the spectrum, the reconstruction accuracy is improved. However, setting $\delta$ to 1 to consider the complete spectrum, very poor reconstructions are observed. From the corresponding spectrum, it is seen that the low-frequency spectrum is well retrieved but the high-frequency part suffers from large deviations. That is due to the ill-conditioning mentioned in Section 3.1. The effects of approximation errors in the PWE approach (e.g., due to insufficiently small sampling density) are amplified after the back-propagation stage of the reconstruction algorithm.

4.4. Uncertainty of Factors

Assume the concerned factors follow the distributions described in

Table 3. Description and distribution of input variables. $\mathcal{U}(a,b)$ denotes the uniform distribution with limits a and b, and $\mathcal{N}(\mu ,\tau )$ denotes the normal distribution with mean $\mu$ and standard deviation $\tau$.

Variable	Description	Distribution
$x_p,y_p,z_p$ (mm)	Cartesian coordinates of the probe position	$a_p^{\mathrm{Ref}}+\mathcal{U}(-0.1,0.1)$, a being x, y, or z
${ϵ}_r$	relative permittivity	${ϵ}_r^{\mathrm{Ref}}+{ϵ}_r^{\mathrm{Ref}}\mathcal{U}(-0.1,0.1)$
$\sigma$ (S/m)	conductivity	${\sigma }^{\mathrm{Ref}}+{\sigma }^{\mathrm{Ref}}\mathcal{U}(-0.1,0.1)$
c (dB)	coupling coefficient	$c^{\mathrm{Ref}}+\mathcal{U}(-2,2)$
$\left|E\right|$	amplitude of electric field	${\left|E\right|}^{\mathrm{Ref}}+{\left|E\right|}^{\mathrm{Ref}}\mathcal{N}(0,0.025)$
$\angle E$ (radian)	phase angle of electric field	$\angle E^{\mathrm{Ref}}+\angle E^{\mathrm{Ref}}\mathcal{N}(0,0.025)$

, where the superscript “Ref” denotes the reference or recorded value. The mutual inductance among probes yields the so-called coupling effects. In practice, since the composition of probes may not be exactly the same, the coupling effects can differ for different pairs of probes and are hard to be analyzed due to many involved coupling coefficients. Here, the coupling effects of each probe are analyzed by only taking into account the mutual inductance from its eight neighbour probes. As a result, the measured noisy field is obtained by

$$\mathbf{E}\left(\mathbf{r}\right)={\mathbf{E}}^{\mathrm{Ref}}\left(\mathbf{r}\right)+\sum _{p_x=-1}^1\sum _{p_y=-1}^1\left[\begin{array}{cc}{c}_{p_x,p_y}^{x,x}& c_{p_x,p_y}^{x,y}\\ c_{p_x,p_y}^{y,x}& c_{p_x,p_y}^{y,y}\end{array}\right]{\mathbf{E}}_{p_x,p_y}^{\mathrm{Ref}}\left(\mathbf{r}\right),$$

(11)

where the subscripts $p_x$ and $p_y$ indicate the position of the neighboured probes while ${\mathbf{E}}_{0,0}^{\mathrm{Ref}}$ denote the electric field measured by the concerned probe. Note that the reference of the measured field is computed by (10) when the uncertainty of the probe position and the phantom electromagnetic properties has been taken into account. Since the z-component is usually computed (rather than measured) according to (9), the electric field in (11) is a column vector composed of the x and y component.

To simplify the analysis further, the matrix of coupling coefficients with respect to different probes is assumed to be the same. The reference value of coupling coefficients applied is given by

$${\mathbf{c}}^{x,x}=\left[\begin{array}{ccc}-0.29-0.27\mathrm{i}& -0.58-0.73\mathrm{i}& -0.46-0.33\mathrm{i}\\ 0.46-1.48\mathrm{i}& 100& 0.37-1.18\mathrm{i}\\ -0.33-0.24\mathrm{i}& -0.65-0.83\mathrm{i}& -0.24-0.21\mathrm{i}\end{array}\right]\times {10}^{-2},$$

(12a)

$${\mathbf{c}}^{x,y}=\left[\begin{array}{ccc}-0.3-0.3\mathrm{i}& -0.5-0.8\mathrm{i}& -0.2+0.2\mathrm{i}\\ 1.5-7.0\mathrm{i}& 2.5-3.8\mathrm{i}& 1.6-0.6\mathrm{i}\\ -0.5-1.0\mathrm{i}& -0.5-0.7\mathrm{i}& -0.1+0.3\mathrm{i}\end{array}\right]\times {10}^{-3},$$

(12b)

and $c_{p_x,p_y}^{y,x}=c_{p_x,-p_y}^{x,y}$, $c_{p_x,p_y}^{y,y}=c_{p_x,-p_y}^{x,x}$ are assumed for the matrices ${\mathbf{c}}^{y,x}$, ${\mathbf{c}}^{y,y}$. This assumption is based on the geometric symmetry of the neighboured probes and is valid when the composition of all probes is exactly the same. Here, the self inductance is not considered and thus the corresponding coefficient is set to 1. Transforming the unit into decibel (dB), the uncertainty of the coupling effects is considered by allowing the truly applied coefficient to have a deviation up to 2 dB.

With 500 Monte Carlo trials, the estimated peak 1 $\mathrm{g}$ SAR of the fourth case is plotted in Figure 6 indicating the uncertainty of all factors described in Table 3. The values of each factor are shown as a boxplot, a standardized way of portraying quartiles of a distribution of values [26]. As shown in Section 4.2, with no uncertainty, the 1 $\mathrm{g}$ and peak 10 $\mathrm{g}$ SAR are estimated with a high accuracy for the 4th case. With uncertainties, unbiased estimations are observed when $\delta =0$, but an overestimation is likely to be obtained when $\delta =0.01$, 0.02 or 0.03. The extent of the bias increases with $\delta$. Besides, due to the ill-conditioning, the variance of the estimation increases with $\delta$ as a higher spatial-frequency spectrum is reconstructed. From the variance, one identifies the most influential factors. When $\delta =0$ or 0.01, the uncertainty of relative permittivity and conductivity contributes most to the variation of the estimation. Despite the large bias and variance of the estimation, the influence factors are changed as the z coordinate of the probe and the measurement (amplitude and phase) accuracy of the electric field when $\delta =0.02$ or 0.03. For validation purposes, a further 500 Monte Carlo trials were independently carried out; the obtained results led to the same conclusion.

4.5. Comparison between the Traditional and Fast Measuring Systems

The estimated peak sSAR value of the concerned 11 cases by the traditional system and the measuring system applying the PWE field reconstructions is presented in Figure 7. Both linear and spline interpolation (and extrapolation) are used in the post-processing of the traditional measuring system. While for most cases all approaches yield well approximated estimations, the technique of spline interpolation seems superior to the linear one, most probably because of its greater approximation power. However, for the fifth, sixth, seventh and eighth cases, the estimation by the traditional system is largely biased, especially for the peak 1 $\mathrm{g}$ SAR. In contrast, although the sampling density is smaller (7 mm for the fast system and $\le 5 \mathrm{mm}$ for the traditional one), the approach based on field reconstructions yields more accurate estimations (except in the seventh case in the estimation of peak 1 $\mathrm{g}$ SAR).

5. Conclusions

Based on the analytical function and 11 body phantom cases, the traditional and the fast SAR measuring system are studied by simulating the measurement process. While both systems generate well approximated estimations for most of the cases, the approach based on the PWE field reconstruction seems to have the potential to achieve more accurate estimations than the approach based on interpolation and extrapolation. However, the reconstruction algorithm suffers from ill-conditioning, which leads to a trade-off between the reconstruction accuracy and reliability. Reconstructing a high spatial-frequency spectrum gives the possibility of yielding more accurate solutions, but simultaneously the estimation always suffers from a higher variance, and vice versa. Moreover, the estimation accuracy may vary with emitting sources. Therefore, without the knowledge of the field-reconstruction algorithm and the source, it is hard to conclude which kind of measuring system is superior to the other.

Another challenge follows the above conclusion. Due to security and commercial interests, the supplier of a measuring system in general will not provide the code or the details of the reconstruction (for the concerned fast system) or the interpolation (for traditional system) algorithm. As shown by the given results, the traditional system is not the gold standard. Therefore, it is necessary to develop a methodology to quantify the measurement accuracy barely based on the product. The involved challenges may include the consideration of various antennas in practice when designing reference antennas, the determination of the reference SAR value, and the uncertainty quantification.