1. Introduction
In recent years, sustainability measures for aging infrastructure facilities have become important for realizing sustainable development goals [1]. One example is the occurrence of road subsidence accidents because of the leakage of pollutants attributed to the cracks and deterioration of buried pipes in underground structures in urban areas. Early detection through periodic inspection and maintenance is essential for preventing such accidents; however, the former is time consuming and labor intensive, and therefore, nondestructive investigations [2], [3] are required for maintenance and exploration.
Ground-penetrating radar (GPR) [4], [5], which uses an electromagnetic wave, has recently attracted considerable research attention as a method for geological, land mine, and archaeological site surveys, and it can be applied to overcome the environmental problems. In the previous studies, the pulse response of structures was analyzed using square or slanted cavities in a dispersive medium by combining the fast inverse Laplace transform (FILT) [6]-[8], Fourier series expansion method (FSEM) [9], and multilayer division method (MDM) [9]. The effects of these cavities on the reflection response waveforms were investigated [9].
In numerical technique for previous paper [9], it has been analyzed the multiple matrix product for multilayer media by using MDM. However, there were unstable problem when the number of layers increased. The continued fraction expanded method (CFEM) [10]-[13] was used to analyze multilayered structures for resolving this problem. CFEM was also applied for determining the electromagnetic wave reflection properties of layered media, and the obtained properties were analyzed and compared with those of the matrix method (MM) in homogeneous media.
In this study, we propose a transient analysis technique that combines FILT and CFEM for analyzing multilayered dispersive media [14]. The numerical results were obtained using the time response, internal time response waveforms, and electric field distribution [15] of the reflection component. Moreover, we verify the calculation accuracy of the two types of FILT methods, namely, (PFILT) and Cesaro FILT (the so-called CFILT method) [7], [16], using a convergence test.
2. Method of Analysis
Figure 1 shows the multilayered dispersive media considered in this study. The structure shown in Fig. 1 comprises \(L\) layers of dispersive media (complex permittivity \(\varepsilon_{l}\)) with uniform thickness \(d_{\alpha}\) (\(\alpha = 1 \sim L\)) in the \(x\)-direction for each region. The permeability was assumed to be \(\mu_{0}\) in all regions, and the waveform of the incident pulse at \(x=0\) was assumed to be a sine pulse in the time domain. It can be expressed as [12], [13]
| \[\begin{align} E_{0}^{(i)}(t) = [u(t) - u(t - t_{w})] \sin \{2\pi t/t_{w}\}, \tag{1} \end{align}\] |
where \(t_{w}\) (\(:=1/f_{0}\), \(f_{0}\): center frequency) denotes the pulse width and \(u(t)\) represents a unit step function. The image function \(E_{0}^{(i)}(s)\) for the complex frequency domain is expressed as
| \[\begin{align} E_{0}^{(i)}(s) = \frac{(2\pi/t_{w})} {s^{2} + (2\pi/t_{w})^{2}} (1-e^{-st_{w}}). \tag{2} \end{align}\] |
The electromagnetic field in each region can be expressed by
| \[\begin{align} & E_{z}^{(l)}(s,x) = A_{l} e^{-k_{l} x} + B_{l} e^{k_{l} x}, (l = 0 \sim L), \tag{3} \\ & H_{y}^{(l)}(s,x) = \frac{1}{s\mu_{0}} \frac{\partial E_{z}(s,x)}{\partial x}, \nonumber\\ &\hphantom{H_{y}^{(l)}(s,x)}=\frac{-k_{l}}{s\mu_{0}}(A_{l} e^{-k_{l}x} - B_{l} e^{k_{l}x}), \tag{4} \\ & k_{0} \triangleq s\sqrt{\varepsilon_{0}\mu_{0}} = s/c_{0},\ k_{l}(s)(:=s\sqrt{\varepsilon_{l}\mu_{0}}), \nonumber \end{align}\] |
where \(k_{0}\) and \(c_{0}\) represent the wave number in free space and the speed of light, respectively; further, \(A_{0} = E_{0}^{(i)}(s)\), and \(B_{0}=R(s)\).
 |
Fig. 1 Structure and coordinate system for multilayered dispersive media |
Moreover, \(A_{l}\), \(B_{l}\) represent unknown coefficients determined from the boundary conditions and symbol \(l\) corresponds to each region.
The dielectric constant of the dispersion medium \(\varepsilon_{l}(s)\) is expressed as a Sellmeier trinomial and orientational polarization considering the water ratio, which includes the soil.
| \[\begin{align} \frac{\varepsilon_{l}(s)} {\varepsilon_{0}} \triangleq 1 + \sum_{m=1}^3 \frac{\Omega_{m,l}^{2}} {s^{2} + g_{m,l}s + \omega_{m,l}^{2}} + \frac{\beta_{2,l}} {1+s\beta_{1,l}}. \tag{5} \end{align}\] |
Substituting Eqs. (3) and (4) into the boundary conditions (\(x = d_{t}^{(l)}\)) yields the relationship between the unknown coefficients \(A_{l}\), \(B_{l}\) for each layer as
| \[\begin{align} x = d_{t}^{(l)} : \left. \begin{array}{l} E_{z}^{(l)}(s,x) = E_{z}^{(l+1)}(s,x)\\ H_{y}^{(l)}(s,x) = H_{y}^{(l+1)}(s,x) \end{array} \right\}. \tag{6} \end{align}\] |
Further, Eq. (6) can be expanded as
| \[\begin{align} & A_{l} e^{-k_{l} d_{t}^{(l)}} + B_{l} e^{k_{l} d_{t}^{(l)}} = A_{l+1} e^{-k_{l+1} d_{t}^{(l)}} + B_{l+1} e^{k_{l+1} d_{t}^{(l)}}, \tag{7} \\ & A_{l} e^{-k_{l} d_{t}^{(l)}} - B_{l} e^{k_{l} d_{t}^{(l)}} = q_{l(l+1)} \{A_{l+1} e^{-k_{l+1} d_{t}^{(l)}} - B_{l+1} e^{k_{l+1} d_{t}^{(l)}}\}. \tag{8} \end{align}\] |
The reflection coefficients are derived from Eqs. (7) and (8) using the CFEM. We add Eqs. (7) and (8) to obtain
| \[\begin{align} A_{l} e^{-jk_{l} d_{t}^{(l)}} = {\rm P}_{l(l+1)} \{A_{l+1} e^{-k_{l+1} d_{t}^{(l)}} + \eta_{l(l+1)} B_{l+1} e^{k_{l+1} d_{t}^{(l)}}\}. \tag{9} \end{align}\] |
Considering the differences between Eqs. (7) and (8), we obtain
| \[\begin{align} B_{l} e^{k_{l} d_{t}^{(l)}} = {\rm P}_{l(l+1)} \{\eta_{l(l+1)} A_{l+1} e^{-k_{l+1} d_{t}^{(l)}} + B_{l+1} e^{k_{l+1} d_{t}^{(l)}}\}, \tag{10} \end{align}\] |
where,
| \[\begin{aligned} & {\rm P}_{l(l+1)} \triangleq \frac{1+q_{l(l+1)}}{2},\quad q_{l(l+1)} \triangleq \frac{\mu_{l} k_{x}^{(l+1)}}{\mu_{l+1} k_{x}^{(l)}} = \frac{1}{q_{(l+1)l}}, \\ & \eta_{l(l+1)} \triangleq \frac{1-q_{l(l+1)}}{1+q_{l(l+1)}} = -\eta_{(l+1)l}. \end{aligned}\] |
Thus, the relational equation for this ratio can be obtained as
| \[\begin{align} \frac{B_{l}}{A_{l}} = \frac{{\rm P}_{l(l+1)} \{\eta_{l(l+1)} A_{l+1} e^{-k_{l} d_{t}^{(l)}} + B_{l+1} e^{k_{l} d_{t}^{(l)}}\}} {{\rm P}_{l(l+1)} \{A_{l+1} e^{-k_{l+1} d_{t}^{(l)}} + \eta_{l(l+1)} B_{l+1} e^{k_{l} d_{t}^{(l)}}\}} e^{-2k_{l} d_{t}^{(l)}}. \tag{11} \end{align}\] |
Equation (11) can be summarized to obtain the final equation for the CFEM, which is
| \[\begin{align} \frac{B_{l}}{A_{l}} = \frac{e^{-2k_{l} d_{t}^{(l)}}}{\eta_{l(l+1)}} + \frac{ \{1 - (1/\eta_{l(l+1)})^{2}\}^{-2(k_{l} + k_{l+1})d_{t}^{(l)}} }{ (1/{\eta_{l(l+1)}}) e^{-2k_{l+1} d_{t}^{(l)}} + (B_{l+1}/A_{l+1}) }. \tag{12} \end{align}\] |
The amplitude ratio (reflection coefficient) in the complex frequency domain can be obtained numerically by substituting the unknown coefficient ratio \(B_{l+1}/A_{l+1}\) of denominator into the second term in Eq. (12).
Therefore, the reflection coefficient \(R\) obtained using CFEM was multiplied by Eq. (2) for numerically obtaining the reflected electric field in the normalized time domain using the FILT method as
| \[\begin{align} E_{z}^{(r)}(T) &:= \frac{1}{2\pi j} \int^{\gamma + j\infty}_{\gamma - j\infty} E_{z}^{(r)}(S,X)e^{ST} dS, \nonumber\\ &\hphantom{:} = \frac{e^{a}}{T}\sum_{m=0}^n (-1)^{m} \mathop{\rm Im} [E_{z}^{(r)}(S,X)], \tag{13} \\ &\hphantom{:} = \frac{e^{a}}{T}\sum_{m=0}^n C_{\kappa} (N,m)(-1)^{m} \mathop{\rm Im}[E_{z}^{(r)}(S,X)], \tag{14} \end{align}\] |
where,
| \[\begin{aligned} & X \triangleq x/c_{0} t_{w},\quad T = t/t_{w},\quad S \triangleq \frac{a+j(m+0.5)\pi}{T},\\ & C_{\kappa}(N,m) \triangleq \frac{N!(\kappa + N - m + 1)!} {(\kappa + N)!(N - m + 1)!}. \end{aligned}\] |
Equations (13) and (14) represent the PFILT and CFILT methods, respectively. In addition, \(N\), \(\kappa\), \(S(:=st_{w})\), \(T(:=t/{t_{w}})\), and \(X(:=x/(c_{0} t_{w}))\) represent the truncation mode number of the FILT method, number of terms in the Cesaro summation [16], normalized complex frequency, normalized time, and normalized coordinate, respectively.
3. Numerical Results
In the following analysis, the numerical calculations were set as \(a=4\), \(N=100\), \(f_{0}=1\) GHz, \(L=4\), and normalized thickness \(D_{\alpha}(:=d_{\alpha}/(c_{0} t_{w}))=0.25\). Further, we employed complex dielectric constants \(\varepsilon_{1}(s)\) and \(\varepsilon_{2}(s)\) with a water ratio of 5% and 10% for the dispersion medium in regions, respectively.
Figure 2 shows the comparison of exact solution [14] and CFEM by combined with PFILT method under the above conditions. The results are in good agreement, thereby confirming the validity of the proposed method.
 |
Fig. 2 Comparison of continued fraction expanded method and exact solution |
Figures 3 (a) and (b) show the reflection response waveforms for a comparison of PFILT and CFILT methods, respectively. From in Fig. 3 (a), we can see clearly that the amplitude of the time response becomes smaller in the range of \(0.5 \leq T \leq 1.0\) and \(4.0 \leq T \leq 6.0\) when compared with that for the PFILT and CFILT methods. The peak value tends to become smaller with an increase in the value of parameter \(\kappa\). Figure 3 (b) show reflection response for magnified view of \(10 \leq T \leq 20\) when we computed \(N=50\).
 |
Fig. 3 Reflection response of PFILT and CFILT methods (\(L=4\)) |
Figure 4 shows the convergence of the truncation mode number with respect to the reflected electric field \(E_{z}^{(r)}(T)\), with a fixed \(T=2.5\), under the conditions of Fig. 2. From Fig. 4, we can see the following feature: It can be obtained that the relative error to the extrapolated true value in \(|E_{z}^{(r)}(T)|\) is about 1.6% and 1% for PFILT and CFILT methods, respectively, when \(N \geq 100\). In addition, the convergence speed increases with an increase in parameter \(\kappa\).
 |
Fig. 4 Convergence of the truncation mode number \(1/N\) |
The internal reflected response waveform is analyzed for investigating the physical phenomena of the reflection response waveform obtained in Fig. 3. The reflected wave component in Eq. (3) can be expressed as follows to efficiently utilize the CFEM amplitude ratio for calculating the internally reflected electric field.
| \[\begin{align} \frac{E_{z}^{(\alpha)}(s,x)} {A_{\alpha}(s)} = \frac{B_{\alpha}(s)} {A_{\alpha}(s)} \exp(k_{\alpha}x). \tag{15} \end{align}\] |
Figure 5 shows the reflection response waveforms at the midpoint (\(X_{\alpha} = (2\alpha - 1)d_{\alpha}/2\), \(\alpha = 1 \sim 4\)) of the thickness of each layer under the conditions shown in Fig. 3. In Fig. 5, we can see that the internal response waveforms increase in sequence with a change in observation points. This is considered to be the reflected wave from each layer until the rise time at the fourth layer, after which the amplitude of the reflected wave from the conductor plate is considered to be amplified by the reflected wave at each layer. And also, the peak values of PFILT and CFILT methods are different in the internal response at observation point of the fourth layer. As this reason, we can understand that the number of terms in the denominator of \(C_{\kappa}(N,m)\) in Eq. (14) is larger than that of the numerator.
 |
Fig. 5 Reflected electric field distribution for different observation points |
Finally, we analyze the reflected electric field distribution at \(T=2.5\) and 3.022 under the conditions as shown in Fig. 5 for further investigating the above results. We compared the PFILT and CFILT methods, and the spatial distribution was normalized to the maximum value of the analysis domain.
Figures 6 (a) and 6 (c) were observed that there is no difference between the PFILT and CFILT results, which is also observed in Figs. 6 (b) and 6 (d). We can see clearly that reflections from the conductor plate are observed in Figs. 6 (a) and 6 (c), and reflection waves in the fourth layer propagate to the third layer in Figs. 6 (b) and 6 (d).
 |
Fig. 6 Distribution of reflected electric field |
Consequently, Fig. 6 shows that, once normalized by the maximum value, the difference in the peak value is irrelevant, and similar results can be obtained for the PFILT and CFILT methods in the distribution of the reflected electric field.
4. Conclusions
In this paper, we propose the transient analysis technique to analyze the multilayered dispersive media by using the combination of the FILT method and the CFEM. Numerical results are given by the reflection response, inside time response waveforms, and the electric field distribution of reflection component. In addition, we verified the calculation accuracy for FILT method of the two types from convergence test.
As a result, we were able to obtain that it is possible to analyze small truncation mode number of FILT method combined with CFEM. In our proposed method, it is effective and powerful, because of indicating same result those of PFILT in distribution of reflected electric field.
We want to apply the CFILT combined with CFEM method to multilayer structures with inhomogeneous distributions. In the future, we will have investigated time response analysis and electromagnetic field distributions with above mentioned structure.
Acknowledgments
This study was partially supported by JSPS KAKENHI grant number JP21K04239.
References
[1] M. Sato, “Social Contributions by Subsurface Electromagnetic Survey,” IEICE Trans. Electron., vol.J103-C, no.3, pp.186-193, March 2020.
[2] Y. Yoshihiro and S. Kidera, “Inverse Scattering Enhanced Imaging Method of Object Buried in Multi-layered Media for Ground Penetrating Radar,” IEICE Tech. Report, AP2022-47, pp.35-40, July 2022.
[3] N. Nishimura and S. Hirose, “Inverse problems and non-destructive testing,” Journal of structural Eng., A, vol.37A, pp.425-434, 1991.
[4] R. Persico, “Introduction to Ground Penetrating Radar ― Inverse Scattering and Data Processing,” Wiley, 2014.
CrossRef
[5] H.M. Jol, “Ground Penetrating Radar Theory and Application,” Elsevier, 2009.
CrossRef
[6] T. Hosono, “Numerical inversion of Laplace transform and some applications to wave optics,” Radio Science, vol.16, no.6, pp.1015-1019, 1981.
CrossRef
[7] T. Hosono, “International series of monographs on advanced electromagnetic: FILT,” vol.2, Science House Inc., 2013.
[8] T. Hosono, “Error Analysis and Improvement of FILT,” IEICE Trans. Electron., vol.J81-C-I, no.4, pp.215-221, April 1998.
[9] R. Ozaki and T. Yamasaki, “Transient Response Analysis of a Dispersive Periodic Grating After Deformation,” URSI Radio Science Letters, vol.4, pp.1-4, 2022.
CrossRef
[10] J.A. Kong, “Electromagnetic Wave Theory,” Wiley, pp.384-389, 1990.
[11] K. Itaya, R. Ozaki, and T. Yamasaki, “Analysis of Reflection Properties for Electromagnetic Waves by a Layered Homogeneous Media,” IEE Japan, 1-016, p.22, March 2023.
[12] K. Itaya, R. Ozaki, and T. Yamasaki, “Fundamental Study on Time Response Analysis of Multilayered Dispersive Media Structure,” IEICE Tech. Report, EMT2023-46, pp.187-192, July 2023.
[13] K. Itaya, R. Ozaki, and T. Yamasaki, “Analysis of Transient Response for Multilayered Dispersive Media Structure Using Continued Fraction Expanded Method,” IEICE Society Conf. 2023, C-1-3, p.3, Sept. 2023.
[14] R. Ozaki and T. Yamasaki, “Transient Scattering Analysis of Electromagnetic Waves for Dispersion Media,” IEICE Tech. Report, EMT2022-58, pp.79-84, Nov. 2022.
[15] J. Sonoda, T. Kon, M. Sato, and Y. Abe, “Characteristics of Detection for Cavity under Reinforced Concrete Using Ground Penetrating Radar by FDTD Method,” IEICE Trans. Electron., vol.J100-C, no.8, pp.302-309, 2017.
[16] T. Hosono, “Cesaro’s Method of Summation and Its Application to Cherenkov Radiation,” IEICE Trans. Electron., vol.J82-C-I, no.3, pp.101-108, March 1999.
