Fast and Accurate Prediction of Antenna Reflection Coefficients in Planar Layered Media Environment via Generalized Scattering Matrix

Parameter measurement of planar layered media is a topic of considerable interest, with wide-ranging applications in pavement thickness assessment [1, 2, 3], moisture detection [4, 5, 6], density evaluation [7, 8], underground coal mining [9], and beyond. These applications often involve solving inverse problems—inferring medium properties from the reflection coefficients of antennas. Modern approaches to inverse problems typically rely on optimization algorithms, which iteratively solve the forward problem to approximate the inverse solution. This makes the computational efficiency of the forward solver a critical factor in the overall performance.

The forward modeling of antenna behavior near planar layered media is commonly based on the spatial-domain dyadic Green’s function for layered media (DGLM), expressed via Sommerfeld integrals [10, 11, 12, 13, 14]. This formulation has been adopted in commercial electromagnetic simulation tools, such as Altair FEKO [15], which employ surface integral equations based on discretized antenna geometries to accurately compute reflection coefficients. However, the high computational cost of a single full-wave simulation significantly limits its usability in optimization-based scenarios that require numerous forward evaluations.

To address this limitation, Lambot et al. pioneered a method that models the antenna as an equivalent point source [16], using its far-field radiation to approximate the reflection response from the layered medium—still based on spatial-domain DGLM. Although evaluating spatial-domain DGLM remains challenging due to the presence of infinite-limit Sommerfeld integrals, the point-source approximation significantly reduces the computational burden by requiring evaluation at only a single location. A key limitation, however, is the assumption that the medium is placed in the far field of the antenna. This was later extended in [17], where a multi-point source model was introduced to enable accurate reflection prediction in near-field scenarios. While this approach enhances accuracy, it also increases the number of DGLM evaluations and requires an additional optimization procedure to determine source positions, introducing nontrivial computational overhead.

More recently, Lambot and colleagues proposed a plane-wave spectrum based technique for estimating antenna reflection coefficients [18]. However, their method neglects evanescent wave contributions and multiple reflections between the antenna and the medium, restricting its applicability to scenarios involving weakly scattering antennas.

Motivated by these limitations, this paper presents a novel method for evaluating the electromagnetic behavior of antennas in the presence of planar layered media, with an emphasis on supporting frequency-domain applications that require repeated reflection coefficient evaluations. These might cover most ground-penetrating radar (GPR) operations and electromagnetic material characterization techniques based solely on reflection data [19, 20]. Our approach recharacterizes antenna transmission, reception, and scattering using the generalized scattering matrix (GSM) framework, expressed in terms of discrete spherical vector wave functions (SVWFs), as opposed to continuous plane-wave spectrum expansions. The interaction between the antenna and the layered medium is modeled algebraically by transforming between SVWFs and planar vector wave functions (PVWFs) [21, 22, 23].

Although our method involves a seemingly complex spherical-planar-spherical waveform transformation, the actual computational effort is modest. While the transformation matrix still requires Sommerfeld integrals, the integration limits are very finite—effectively reducing the complexity to that of standard definite integrals. Compared to the spectral-domain method in [18], which involves double integrals, our approach only requires single integrals. Furthermore, the matrix’s inherent symmetry and sparsity significantly improve numerical efficiency.

As this paper further demonstrates, the SVWF-based framework allows for a compact problem size and makes the inclusion of multiple reflections computationally inexpensive. This enables high-fidelity modeling without compromising accuracy for the sake of speed. Ignoring the one-time preprocessing cost of the free-space antenna characterization, our method achieves several orders of magnitude speedup over full-wave simulations in FEKO, while maintaining near-identical accuracy. This substantial advancement in forward problem efficiency could greatly facilitate the adoption of modern computational techniques in inverse problems involving planar layered media.

A core foundation of this paper is the use of PVWFs and SVWFs, which are extensively utilized in [21, 22, 23, 24]. We adopt the time convention $e^{\mathrm{j}\omega t}$ and define the PVWF $\bm{\phi}_{i}\left(\alpha,\beta,\bm{r}\right)=\bm{\phi}_{i}\left(\hat{\gamma},\bm{r}\right)$ as follows:

Here, $i=1$ represents the TE wave, and $i=2$ represents the TM wave. The unit vectors $\hat{\alpha},\hat{\beta},\hat{\gamma}$ are related to the elevation angle $\alpha$ and azimuth angle $\beta$ as follows:

The azimuth angle $\beta$ ranges from $[0,2\pi]$, while the elevation angle $\alpha$ spans $C_{+}\cup C_{-}$. The contours $C_{+}$ and $C_{-}$ are illustrated in Fig. 1a. For propagating waves, $\alpha$ lies along the real axis between 0 and $\pi$.

In this paper, we also employ SVWFs $\bm{u}_{n}^{(p)}(\bm{r})$, derived from real-valued spherical harmonics [25, Appendix A]. Here, $p=1,4$ denote regular and outgoing waves, respectively. The subscript $n\rightarrow\tau\sigma ml$ indicates polarization type $\tau={1,2}$ (TE, TM), parity $\sigma={e,o}$ (even, odd), angular degree $l={1,2,\dots,L_{\max}}$, and order $m={0,1,\dots,l}$.

Both SVWFs and PVWFs satisfy the homogeneous vector wave equation, and therefore, they can be converted into one another [24]:

Here, $z>0$ corresponds to $C_{+}$, and $z<0$ corresponds to $C_{-}$. The integral notation is defined as $\int_{C}{\mathrm{d}\hat{\gamma}}=\int_{C}{\sin\alpha\mathrm{d}\alpha\int_{0}^{2\pi}{\mathrm{d}\beta}}$. The function $B_{ni}(\hat{\gamma})$ is given by:

where

Here, the upper (lower) form of $A_{ni}(\beta)$ applies for even (odd) parity $\sigma=e(o)$. The $\delta$ is the Kronecker delta. A bar over $i$ swaps the TE/TM types of PVWFs (i.e., $\bar{i}=3-i$). $B_{ni}^{\dagger}(\hat{\gamma})$ represents the same expression as $B_{ni}(\hat{\gamma})$ but with all explicit $\mathrm{j}$ replaced by $-\mathrm{j}$. The auxiliary functions

where $\tilde{P}_{l}^{m}(u)$ represents the normalized associated Legendre functions [26]. An efficient algorithm for evaluating these functions for all $l$ less than $L_{\max}$ is provided in [26]. This can also be extended to evaluate $\tilde{P}_{l}^{m}(u)$ for imaginary $u$.

A key integral identity used later this work is:

where $I_{\tau\sigma,\tau^{\prime}\sigma^{\prime}}^{i}$ is nonzero only in the following cases:

The GSM for antennas in free space (the source-scattering formulation) is defined as [27, 28]:

where $\mathbf{\Gamma}$ represents the port S-parameters of the antenna (for a single-mode and single-port antenna, this corresponds to the reflection coefficient), $\mathbf{R}$ is the receiving matrix, $\mathbf{T}$ is the transmitting matrix, and $\mathbf{S}$ is the SVWFs scattering matrix. This paper, building upon previous work [29], uses the method of moments (MoM) to determine the GSM of the antenna.

The GSM relates the incident and response fields of the antenna as:

where $\mathbf{v}$ and $\mathbf{w}$ are the expansion coefficients for the incoming and outgoing feeding modes at the antenna ports, respectively. $\mathbf{a}$ represents the expansion coefficients for the regular spherical wave expansion of the external incident wave $\bm{E}^{\mathrm{inc}}$, and $\mathbf{f}$ are the expansion coefficients for the outgoing spherical waves of the antenna’s radiated (scattered) field $\bm{E}^{s}$, specifically:

Consider a planar interface located at $z=z_{\mathrm{I}}$ (with $z_{\mathrm{I}}<0$), as shown in Fig. 2. Further assume that the antenna is actively fed only by wave ports, such that the incident wave entering the antenna from the external space is exclusively the reflected field of $\bm{E}^{s}$, i.e., $\bm{E}^{\mathrm{inc}}=\bm{E}^{s\mathrm{R}}$. To determine $\bm{E}^{s\mathrm{R}}$, we must calculate how each SVWF component $\bm{u}_{n}^{(4)}(\bm{r})$ in (9b) is reflected by the interface. To do so, we transform $\bm{u}_{n}^{(4)}(\bm{r})$ into downward PVWFs using the relation in (3):

The reflected wave of the PVWF $\bm{\phi}_{i}(\alpha,\beta,\bm{r})$, generated by the interface, is given by

where $\rho_{i}^{\text{Ref}}\left(\alpha\right)=\rho_{i}^{1}(\alpha)e^{-2\mathrm{j}k\hat{\gamma}\cdot\bm{d}}$, and $\rho_{i}^{1}(\alpha)$ is the well-known Fresnel reflection coefficient [see Appendix A]. The phase shift, referenced from the origin, is given by $e^{-2\mathrm{j}k\hat{\gamma}\cdot\bm{d}}$, where $\hat{\gamma}\cdot\bm{d}=z_{\mathrm{I}}\cos\alpha$. Thus, the reflected wave of $\bm{u}_{n}^{(4)}(\bm{r})$ can be written as

Equation (3) also provides the SVWF representations for $\bm{\phi}_{i}\left(\pi-\alpha,\beta,\bm{r}\right)$:

where $a^{\prime}=l^{\prime}+m^{\prime}+\tau^{\prime}+i+1$. Consequently, the reflected field $\bm{u}_{n}^{\mathrm{R}}\left(\bm{r}\right)$ can be expressed as

where the $nn^{\prime}$-th element of the matrix $\mathcal{W}$ is given by:

The integral over $\beta$ in (15) can be computed analytically, resulting in the identity in (6). By substituting $u=\cos\alpha$, we arrive at:

The Sommerfeld integral in (16) traditionally starts from $\mathrm{j}\infty$ on the contour $C_{-}^{\prime}$ in Fig. 1b. However, recent research suggests that the starting point of this integral should be truncated by a very finite value $\mathrm{j}\tilde{\kappa}_{m}(\tilde{\kappa}_{m}>1)$, rather than a large value approaching infinite. This modification significantly reduces the computational complexity of (16) while ensuring convergence everywhere above the planar interface. An empirical formula for determining $\tilde{\kappa}_{m}$ is from [30] as:

which holds for $0.5\leq kR_{\min}\leq 10,L_{\max}\leq 20$. The maximum degree number $L_{\max}$ is usually chosen as [31]:

where $R_{\min}$ is the minimal radius of the circumscribing sphere around the antenna, and $\iota\in[2,7]$. One can prove the symmetry property $\mathcal{W}_{nn^{\prime}}=\mathcal{W}_{n^{\prime}n}$, which further reduces the computational burden for evaluating $\mathcal{W}$.

With the above results, we can express $\bm{E}^{s\mathrm{R}}(\bm{r})$ as

or more compactly:

where

Because GSM connects $\mathbf{v},\mathbf{w}$ and $\mathbf{a},\mathbf{f}$ as

we substitute $\mathbf{a}=\mathcal{W}\mathbf{f}$ into the second equation of (22) to obtain

Solving for $\mathbf{f}$ gives

where $\mathbf{G}\equiv\frac{1}{2}\mathcal{W}$. Substituting this expression into the first equation of (22) leads to the following relation:

This implies that the S-parameters matrix of the antenna near the layered planar medium is

The Neumann series expansion of $\left[\mathbf{1}-\left(\mathbf{S}-\mathbf{1}\right)\mathbf{G}\right]^{-1}$ reads

Each term in this expansion represents interactions of different orders with the layered planar medium. This equation allows for an explicit separation of the contributions from various orders of ground echoes. If multiple reflections between the antenna and the interface are neglected, equation (27) can be approximated by taking only the first term (i.e., the identity matrix $\mathbf{1}$). However, this approximation holds only when the antenna’s scattering is small ($\mathbf{S}\to\mathbf{1}$). For greater precision, it is usually necessary to consider higher-order terms in the multiple reflections or to directly calculate $\left[\mathbf{1}-\left(\mathbf{S}-\mathbf{1}\right)\mathbf{G}\right]^{-1}$.

An interesting observation about equation (27) is that $\mathbf{\Gamma},\mathbf{R},\mathbf{T},\mathbf{S}$ are matrices associated with the antenna’s behavior in free space. In contrast, $\mathbf{G}$ is a matrix that depends on the layered planar medium and the position of the interface $z_{\mathrm{I}}$, but is independent of the antenna’s characteristics. Therefore, when either $z_{\mathrm{I}}$ or the medium parameters change, only $\mathbf{G}$ needs to be reevaluated. This process is computationally efficient and can typically be completed within milliseconds, as further detailed in Section IV-D.

This section demonstrates the numerical performance of the proposed formulas in various planar layered environments. A multi-mode horn antenna, as described in [29], is used as the transmitter. The coordinate origin (geometric center of the antenna) is positioned 80 mm from the radiation aperture. Over the frequency range from 3.2 GHz to 3.8 GHz, the antenna operates in five modes: $\text{TE}_{10}$, $\text{TE}_{20}$, $\text{TE}_{01}$, $\text{TE}_{11}$, and $\text{TM}_{11}$. In the following examples, unless otherwise stated, the maximum degree of SVWFs used is $L_{\max}=12$, and the integration truncation number is $\tilde{\kappa}_{m}=1.5$.

This paper introduces a novel numerical method for computing the reflection coefficient of antennas near planar layered media. The method models the transmitting, receiving, and scattering behavior of antennas through a GSM defined by spherical vector wave functions. By expanding SVWFs into PVWFs and applying Fresnel reflection coefficients for the layered medium, the interaction between the antenna and the layered structure is evaluated. Importantly, no approximations are made in this process, ensuring that the computational accuracy matches that of FEKO.

The proposed method benefits from low matrix dimensionality and reduced algebraic complexity. When excluding the one-time preprocessing step to compute the antenna’s GSM in free space, the computation time for each layered media configuration is on the order of milliseconds—vastly outpacing full-wave simulations in FEKO. This makes the method particularly suitable for real-time applications and presents a promising new approach for solving inverse problems involving layered media.