Theory of Microwave Remote Sensing of Vegetation Effects, SoOp and Rough Soil Surface Backscattering

Abstract

In this paper, we provide updates on our recent work on the theory of microwave remote sensing for applications in remote sensing of soil moisture and snow water equivalent (SWE). The three topics are the following. (i) For the effects of forests and vegetation, we developed the hybrid method of NMM3D full-wave simulations over the vegetation field and forest canopies. In the hybrid method, we combined the use of commercial off-the-shelf software and wave multiple scattering theory (W-MST). The results showed much larger transmission than classical radiative transfer theory. (ii) In signals of opportunity at L-band and P-band, which are radar bistatic scattering in the vicinity of the specular direction, we developed the Analytical Kirchhoff solution (AKS) and Numerical Kirchhoff approach (NKA) in the calculations of coherent waves and incoherent waves. We also took into account of the effects of topographical elevations and slopes which have strong influences. (iii) In rough surface radar backscattering, we used the volume integral equation approach for NMM3D full-wave simulations for soil surfaces with kh up to 15. The simulations were calculated for the X-band and Ku-band and the results showed saturation effects. The simulation results can be applied to microwave remote sensing of SWE at these two frequencies.

1. Introduction

Soil moisture and snow water equivalent are two geophysical variables that are monitored globally and continually. Information of soil moisture has applications in weather forecasting, modeling of climate change, agricultural productivity, water resources management, drought prediction, flood area mapping, and the ecosystem. Seasonal snow on land is responsible for processes and feedbacks that affect the global climate system, freshwater availability to billions of people, biogeochemical activity including exchanges of carbon dioxide and trace gases, and ecosystem services.

Presently, the ESA SMOS satellite and NASA SMAP satellite have been using L-band radiometry in measuring soil moisture. The NISAR satellite that will be launched in 2023 will use both L-band and S-band radar backscattering. The ESA’s Sentinel 1 radar at C-band has also been used to map soil moisture. The Copernicus Imaging Microwave Radiometer (CIMR) will be launched by ESA. There are five frequencies in CIMR of which the L-band, C-band, and X-band will be useful for the monitoring of soil moisture. Recently, P-band reflectometry is also proposed for the remote sensing of soil moisture. In microwave remote sensing of soil moisture, a major challenge is the ability of microwaves to penetrate through the vegetation field and forest canopies above the soil. The original proposal of SMAP at L-band was that the sensitivity to soil moisture will be limited by the vegetation field and forest canopies with the upper limit of VWC (Vegetation Water Content) at 5 kg/m². This upper limit is based on the models of radiative transfer equation (RTE) and the distorted Born approximation (DBA) using the discrete scatterer model. The RTE and DBA approaches have been used since the 1980s and are extensions of the water cloud model [1]. Based on the two models, attenuation through vegetation will increase with frequency, making it even more difficult for remote sensing of soil moisture at the higher frequencies of S-, C-, and X-bands. However, in recent years, we have been successful in performing full-wave simulations of Numerical solution of Maxwell equations in 3D vegetation based on the hybrid method [2,3,4,5]. The results of full-wave simulations showed much larger transmission than that of the results of RTE/DBA.

In P-band and L-band Signals of Opportunity (SoOp), the transmitters on existing satellites are utilized. Satellites with receivers are launched to measure the reflected signals. Global Navigation Satellite System Reflectometry (GNSS-R) is an application of SoOp at the L-band. The operating GNSS-R missions include the Techdemosat-1 (TDS-1) [6] launched by UK in 2014, the Cyclone Global Navigation Satellite System (CYGNSS) [7] launched by NASA in 2016, and Bufeng launched by China in 2019. The GNSS-R data are collected in the form of Delay Doppler Maps (DDMs). Researchers have shown the potential of soil moisture retrieval by the GNSS-R data [8,9,10]. Recently there are also interests in using the phase of the P-band signals of opportunity at 370 MHz. In Reference [11], it was proposed to infer snow water equivalent using the phase of the reflected waves. With new SAR technology [12] implemented for Signals of Opportunity with much higher spatial resolutions, the use of SoOp will gain importance for remote sensing in land applications. The distinction between SoOp and radar backscattering is that in SoOp, the bistatic direction is in the vicinity of the specular direction within a few degrees. The signals are influenced strongly by topography. In theoretical modeling for signals of opportunity, we developed the Analytical Kirchhoff solution (AKS) and Numerical Kirchhoff approach (NKA) in first-principles calculations of coherent waves and incoherent waves [13,14,15,16].

In microwave remote sensing of snow water equivalent (SWE), the X- and Ku-band radar backscattering measurements provide the means to produce SWE information at the temporal and spatial scales necessary to advance water resource management across the northern hemisphere [17]. The ESA Cold Regions Hydrology High-Resolution Observatory (CoReH2O) mission (dual-frequency X- and Ku-band; completed Phase A at ESA in 2013) was a major impetus [18]. The physical basis for estimating SWE from X- and Ku-band radar measurements is volume scattering by millimeter-scale snow grains. Significant progress was made over the past decade in understanding the X-band and Ku-band radar response to variations in SWE, snow microstructure, and snow wet/dry state. Driven by the success of recent air-borne measurements, ground-based measurements, [19,20,21] and improved physical models and validations of retrieval algorithms [22,23,24], a Terrestrial Snow Mass Mission (TSMM) in phase 0 has been supported by the Canadian Space Agency. In addition to volume scattering, another contribution to the X-band and Ku-band radar backscattering signal is the rough surface backscattering from the soil surface below the snow. The backscattering is the sum of the volume scattering contribution and the rough surface scattering contribution. The rough-surface scattering is a “nuisance effect” as it is not correlated with SWE and needs to be subtracted to obtain the volume scattering contribution. With the advent of computers and computation methods, full-wave simulations of numerical solutions of Maxwell equations for rough soil-surface scattering have been carried out for L-bands up to kh = 3 where k is the wavenumber and h is the rms height [25,26,27]. In mountainous regions, the rms height of rough surfaces can be 6 cm so that kh = 25 at the Ku band frequency of 17.2 GHz

Recently, we performed full-wave simulations to apply to rough-surface scattering at X-band and Ku-band, up to kh = 15 using the volume integral equation (VIE) approach. The results showed rough-surface scattering saturation with increase of rms height and frequency.

The focus of this paper was to present the exciting recent results from our group on the theory of microwave remote sensing. In the three topics, (i) the full-wave simulations of NMM3D are recent works that can provide different results from RTE. It is significant because the models of RTE and a similar model of DBA have been used for several decades. (ii) For the signals of opportunity, the SoOp topic is new as SoOp missions were only recently launched. The new microwave theory is that the scattered direction is in the vicinity of the specular direction. The theory shows that, unlike radar backscattering, specular reflection is affected by topography strongly and there are both coherent waves and incoherent waves. (iii) For rough-surface scattering, the new theory is for kh up to 15 while old models are limited to kh = 3. The new theory means that the results can be applied up to the Ku-band.

In Section 2, we study the effects of forests and vegetation in microwave remote sensing of soil moisture. We describe the hybrid method of NMM3D full-wave simulations over the vegetation scene and forests scene. In Section 3, we describe the model of signals of opportunity. We describe the Analytical Kirchhoff solution (AKS) and Numerical Kirchhoff approach (NKA) in first-principles calculations of coherent waves and incoherent waves in the vicinity of the specular scattered direction. In addition to microwave centimeter roughness, we also took into accounts of the effects of topographical elevations and slopes. The results of NKA and AKS were indistinguishable. Comparisons were also made with the results of two geometric optics methods. In Section 4, we describe NMM3D full-wave simulations of rough soil-surface scattering for soil surfaces with kh up to 15. The simulations were applied for the X-band and Ku-band to calculate rough-surface scattering of the soil surface with and without a layer of snow above the soil. We also illustrate the retrieval of rough surface rms height and soil moisture using UAVSAR L-band data. Section 5 is the Conclusion.

2. Vegetation and Forest Effects in Microwave Remote Sensing of Soil Moisture

In microwave remote sensing of soil moisture, a major challenge is the capability of microwaves to penetrate through the vegetation and forest canopies. The accurate modeling of wave propagation at multiple frequencies in the vegetation layer above the soil is also crucial for evaluating the feasibility of the multifrequency approach of remote sensing of soil moisture.

Since the 1980s, the discrete scatterer approach has been implemented using Radiative Transfer Equation (RTE) and the Distorted Born Approximation (DBA). In these two approaches, the vegetation canopy is considered as a layer of randomly distributed scatterers. Specifically, the branches and leaves are treated as single scatterers and are modeled, respectively, by cylinders and disks. The positions of branches and leaves are assumed to be random and are assumed to be statistically homogeneous in space. The scattering and absorption cross sections are calculated by adding the scattering cross sections and absorption cross sections of each branch and each leaf. The orientations of disks and cylinders are characterized by probability distribution over the Euler angles of rotations [28].

In RTE and DBA, the vegetation or forest canopy is approximated as a layered medium with an effective attenuation rate calculated based on the independent scattering approximation. The effective attenuation rate is basically an assumption of homogenization treating the vegetation canopy as effectively homogeneous. In this model of backscattering, there are three contributions which are rough-surface scattering, volume scattering, and double bounce [29,30]. The models of first-order RTE and DBA are essentially the same as they have the same formulas of attenuation, optical thickness tau, and bistatic volume scattering. The only difference between RTE and DBA is a factor of 2 in the double bounce term, which is due to backscattering enhancement [22,31,32,33]. The RTE and DBA have several inherent approximations when applied to vegetation and forests.

(1) The positions of the scatterers are assumed to be uniformly random. However, such a description is not consistent with the clustering of scatterers in a tree or plant. For example, in pine trees, there are hundreds of needle leaves attached to a branch. In addition, the assumption of uniformly random positions does not account for the gaps among plants and trees. Unlike ray optics, the abilities of microwave to pass through the gaps depend on the sizes of the gaps relative to the microwave wavelength.

(2) The scatterers such as leaves and branches are assumed to be independent so that the scattered intensities are added incoherently. The frequency dependence of scattering is that of the individual cylinder and/or individual disk, which leads to a strong increase of scattering and absorption by vegetation. The collective scattering and absorption effects due to clustering are ignored.

(3) Far-field assumptions are made throughout the RTE and DBA models. This means far field between the scatterers and between the scatterers and the underlying soils. Since the far field distance is size squared divided by wavelength, the far field assumptions are not obeyed in vegetation/forests.

Recently, we were successful in performing full-wave solutions of Numerical Maxwell Model of 3D Simulations. The method is a hybrid method combining Computational Electromagnetics (CEM) of a single object with Wave Multiple Scattering Theory (W-MST). We added a word “Wave” in front of MST to distinguish from the multiple scattering theory of the Radiative Transfer Equations. The MST in RTE is intensity addition rather than wave interactions with amplitudes and phases. In the following subsections, we first describe the T matrix of a single scatterer that can be calculated from CEM. This is followed by the W-MST based on Foldy-Lax equations.

In addition to remote sensing, the MST model has been extensively used in solid state physics, photonics, and acoustics.

2.1. T Matrix of a Plant or a Tree

The W-MST formulation is a Foldy-Lax Multiple Scattering equation using the T matrices of single scatterers. In the hybrid method, we treat a single plant such as a corn plant, a wheat plant, or a soya bean plant, or a single tree as a single scatterer. This is different from the RTE and DBA in that a single branch or a single leaf is treated as a single scatterer. The single scatterer T matrix is the scattering T matrix of an isolated scatterer in the absence of all other scatterers.

Consider cylindrical coordinates $\left(\rho ,\varphi ,z\right)$. The Vector Cylindrical Waves (VCW) that are outgoing are

$${\overline{M}}_n\left(k_z,\overline{r}\right)=\left[\widehat{\rho }\frac{in}{\rho }{H}_n^{\left(1\right)}\left(k_{\rho }\rho \right)-\widehat{\varphi }{k}_{\rho }{H}_n^{{\left(1\right)}^{\prime }}\left(k_{\rho }\rho \right)\right]exp\left(i{k}_{z}z+in\varphi \right)$$

(1)

$${\overline{N}}_n\left(k_z,\overline{r}\right)=\left[\widehat{\rho }\frac{i{k}_{\rho }{k}_z}{k}{H}_n^{{\left(1\right)}^{\prime }}\left(k_{\rho }\rho \right)-\widehat{\varphi }\frac{n{k}_z}{k\rho }{H}_n^{\left(1\right)}\left(k_{\rho }\rho \right)+\widehat{z}\frac{k_{\rho }^2}{k}{H}_n^{\left(1\right)}\left(k_{\rho }\rho \right)\right]exp\left(i{k}_{z}z+in\varphi \right)$$

(2)

where $k$ is the wavenumber, $k_{\rho }=\sqrt{k^2-k_z^2},-\infty <k_z<\infty ,n=0,\pm 1,\pm 2,\cdots ,$ and $H_n^{\left(1\right)}\left(k_{\rho }\rho \right)$ are Hankel functions of the first kind and of order $n$. There are two indices, $n$ and $k_z$, for the wave functions. The regular wave functions are the expressions above which Hankel functions are replaced by Bessel functions $J_n\left(k_{\rho }\rho \right)$.

$$Rg{\overline{M}}_n\left(k_z,\overline{r}\right)=\left[\widehat{\rho }\frac{in}{\rho }{J}_n\left(k_{\rho }\rho \right)-\widehat{\varphi }{k}_{\rho }{J}_n^{\prime }\left(k_{\rho }\rho \right)\right]exp\left(i{k}_{z}z+in\varphi \right)$$

(3)

$$Rg{\overline{N}}_n\left(k_z,\overline{r}\right)=\left[\widehat{\rho }\frac{i{k}_{\rho }{k}_z}{k}{J}_n^{\prime }\left(k_{\rho }\rho \right)-\widehat{\varphi }\frac{n{k}_z}{k\rho }{J}_n\left(k_{\rho }\rho \right)+\widehat{z}\frac{k_{\rho }^2}{k}{J}_n\left(k_{\rho }\rho \right)\right]exp\left(i{k}_{z}z+in\varphi \right)$$

(4)

In Figure 1, we show a corn plant with height 1.25 m, and the extent of the leaves is to 0.4 m. For the calculation of T matrices from commercial software (such as FEKO or HFSS), we enclose the plant in a cylindrical volume of radius $R_C=0.45$ m and height 1.3 m as shown in Figure 1. The T matrices results are for the use in Foldy–Lax multiple scattering equations with other corn plants. Other corn plants are enclosed by similar cylinders, and they lie outside the enclosing cylinders of each other.

The T matrices in vector cylindrical waves are for the region $\rho >R_C$ to treat wave interactions among plants. The fields above the plant and below the plant for $\rho \le R_C$ do not enter Foldy–Lax multiple scattering equations. After Foldy–Lax equations are solved, the fields in these “inner” regions, $\rho \le R_C$, can be calculated by Huygen’s principles. In Figure 1, the $\beta$ is the orientation angle and $\alpha$ is the azimuthal angle for the leaf. They are labeled to characterize the orientations of the stalks and the leaves. These are used traditionally in RTE as each stalk is treated as a single scatterer and each leaf is treated as a single scatterer. Both $\beta$ and $\alpha$ of leaves and stalks obey orientation distributions. Generally, $\alpha$ is assumed to have uniform distribution between 0 and $2\pi$. In the hybrid method, we perform orientations differently as we treat a plant as a single scatterer. Each plant is vertical. Thus, the plant itself does not have any $\beta$. The plant does have the azimuthal angle which is also assumed to be uniform between 0 and $2\pi$. We construct a corn field of many plants such that the orientations of leaves and stalks are described by the same distributions as done traditionally. We have several plants such that the distributions of $\beta$ angles of the stalks and leaves agree with traditional characterizations. The T matrices of these several plants are calculated with commercial software. The $\alpha$ angles are next introduced. The T matrices with $\alpha$’s are just rotations of the T matrices without the need of calculations of new T matrices.

The plant or tree is within a “virtual” cylinder of radius $R_C$ (Figure 1) based on Wave Multiple Scattering Theory (W-MST), the electric field outside the enclosing cylinder, that is, $\rho \ge R_C$ is the sum of ${\overline{E}}^{ex}$, the “final” exciting field, and ${\overline{E}}^s$, the “final” scattered field. Both ${\overline{E}}^{ex}$ and ${\overline{E}}^s$ refer to that of the “single scatterer“, which is inside the enclosing cylinder.

The final exciting field and final scattered field are respectively given by

$${\overline{E}}^{ex}\left(\overline{r}\right)={\displaystyle \sum }_n{{\displaystyle \int }}_{-\infty }^{\infty }d{k}_z\left[a_n^{EM}\left(k_z\right)Rg{\overline{M}}_n\left(k_z,\overline{r}\right)+a_n^{EN}\left(k_z\right)Rg{\overline{N}}_n\left(k_z,\overline{r}\right)\right]$$

(5)

$${\overline{E}}^s\left(\overline{r}\right)={\displaystyle \sum }_n{{\displaystyle \int }}_{-\infty }^{\infty }d{k}_z\left[a_n^{SM}\left(k_z\right){\overline{M}}_n\left(k_z,\overline{r}\right)+a_n^{SN}\left(k_z\right){\overline{N}}_n\left(k_z,\overline{r}\right)\right]$$

(6)

In W-MST, the T matrix is that of a single scatterer, which is the scattering of the scatterer in the absence of all other scatterers. The scattered field coefficients are expressed in terms of the T matrices.

$$a_n^{SM}\left(k_z\right)={\displaystyle \sum }_{n^{\prime }}{{\displaystyle \int }}_{-\infty }^{\infty }d{k}_z^{\prime }\left[T_{n{n}^{\prime }}^{\left(M,M\right)}\left(k_z,k_z^{\prime }\right)a_{n^{\prime }}^{EM}\left(k_z^{\prime }\right)+T_{n{n}^{\prime }}^{\left(M,N\right)}\left(k_z,k_z^{\prime }\right)a_{n^{\prime }}^{EN}\left(k_z^{\prime }\right)\right]$$

(7)

$$a_n^{SN}\left(k_z\right)={\displaystyle \sum }_{n^{\prime }}{{\displaystyle \int }}_{-\infty }^{\infty }d{k}_z^{\prime }\left[T_{n{n}^{\prime }}^{\left(N,M\right)}\left(k_z,k_z^{\prime }\right)a_{n^{\prime }}^{EM}\left(k_z^{\prime }\right)+T_{n{n}^{\prime }}^{\left(N,N\right)}\left(k_z,k_z^{\prime }\right)a_{n^{\prime }}^{EN}\left(k_z^{\prime }\right)\right]$$

(8)

where the T matrix coefficients are $T_{n{n}^{\prime }}^{\left(M,M\right)}\left(k_z,k_z^{\prime }\right)$, $T_{n{n}^{\prime }}^{\left(M,N\right)}\left(k_z,k_z^{\prime }\right)$, $T_{n{n}^{\prime }}^{\left(N,M\right)}\left(k_z,k_z^{\prime }\right)$, and $T_{n{n}^{\prime }}^{\left(N,N\right)}\left(k_z,k_z^{\prime }\right)$. The T matrix coefficients have dimensions of length for this case of cylindrical waves.

The self-consistent equations of W-MST allow ${\overline{E}}^{ex}$ and ${\overline{E}}^s$ to be calculated in a self-consistent manner using the single scatterer T matrix. Since $-\infty <k_z<\infty$, and $k_{\rho }=\sqrt{k^2-k_z^2}$, the waves are propagating for $\left|k_z\right|\le k$ and evanescent in the $\widehat{\rho }$ direction for $\left|k_z\right|>k$. Evanescent waves decay exponentially. At $1/4$ wavelengths, the evanescent waves become small. At the SMAP frequency of $1.41 \mathrm{GHz}$, this $1/4$ wavelength corresponds to $5.3 \mathrm{cm}$. This means if the neighboring circumscribing cylinder is separated by more than 5 cm, the evanescent waves are negligible. We shall truncate with $\left|k_z\right|\le k$. The angular direction $\theta$ can be used with

$$k_z=k\cos \theta ;{{\displaystyle \int }}_{-k}^{k}d{k}_z=k{{\displaystyle \int }}_0^{\pi }d\theta \sin \theta$$

(9)

In (9) there is a factor, $\sin \theta$, that corresponds to the decrease of the solid angle when $\theta \to 0$, decreasing the contributions of $\theta \to 0$.

The cylindrical wave expansions are analytic wave functions. The use of Bessel functions and Hankel functions in the wave functions ${\overline{M}}_n$ and ${\overline{N}}_n$ mean that they are applicable to all distance ranges of $\rho$ from near field to intermediate distance and to far field. On the other hand, the coefficients of field expansions, which include exciting field coefficients, the scattered field coefficients, and the T matrix coefficients are independent of distance. This means we can obtain these coefficients at any distance including the far field, the intermediate field, and the near field. Once these coefficients are obtained, they can be substituted in the analytic field expansions and the field expansions are valid for all distance ranges. We use far field to extract the T matrix coefficients because commercial software provides far-field scattering amplitudes.

In the far field, we related the T matrix of vector spherical waves to that of far field scattering amplitudes [28]. Below, we relate the T matrix of vector cylindrical waves (VCW) to the far field scattering amplitudes. In the far field of $\rho \to \infty$

$${\overline{M}}_n\left(k_z,\overline{r}\right)\to -\widehat{\varphi }{k}_{\rho }\left(i\sqrt{\frac{2}{\pi k_{\rho }\rho }}exp\left(i\left(k_{\rho }\rho -\frac{n\pi }{2}-\frac{\pi }{4}\right)\right)\right)exp\left(i{k}_{z}z\right)exp\left(in\varphi \right)$$

(10)

$${\overline{N}}_n\left(k_z,\overline{r}\right)\to \left[-\widehat{\rho }{k}_z+\widehat{z}{k}_{\rho }\right]\frac{k_{\rho }}{k}\sqrt{\frac{2}{\pi k_{\rho }\rho }}exp\left(i\left(k_{\rho }\rho -\frac{n\pi }{2}-\frac{\pi }{4}\right)\right)exp\left(i{k}_{z}z\right)exp\left(in\varphi \right)$$

(11)

Substitute into Equation (6) and apply the method of stationary phase. The spherical coordinates are $\left(r,\theta ,\varphi \right).$ The stationary phase points are at $k_z=k\cos \theta$, $k_{\rho }=k\sin \theta$. Since $\widehat{\theta }=\cos \theta \widehat{\rho }-\sin \theta \widehat{z}$, then

$$\begin{gathered} {\overline{E}}^s\left(\overline{r}\right)\to -\widehat{\varphi }k\sin \theta \frac{2exp\left(ikr\right)}{r}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)a_n^{SM}\left(k\cos \theta \right) \\ -\widehat{\theta }k\sin \theta \frac{2exp\left(ikr\right)}{ir}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)a_n^{SN}\left(k\cos \theta \right) \end{gathered}$$

(12)

Let the exciting field be a plane wave in the direction ${\widehat{k}}_i=\sin {\theta }_i\cos {\varphi }_i\widehat{x}+\sin {\theta }_i\sin {\varphi }_i\widehat{y}+\cos {\theta }_i\widehat{z}$ with polarizations vectors ${\widehat{v}}_i={\widehat{\theta }}_i$ of vertical polarization and ${\widehat{h}}_i={\widehat{\varphi }}_i$ of horizontal polarization. Let the amplitudes be $E_{vi}$ and $E_{hi}$. Then, from Tsang et al., volume 1 [28]

$$\begin{gathered} {\overline{E}}^{ex}\left(\overline{r}\right)=\left({\widehat{v}}_i{E}_{vi}+{\widehat{h}}_i{E}_{hi}\right)exp\left(i{\overline{k}}_i\cdot \overline{r}\right) \\ ={\displaystyle \sum }_n\frac{i^{n}exp\left(-in{\varphi }_i\right)}{k_{i\rho }}\left[i{E}_{hi}Rg{\overline{M}}_n\left(k_{i\rho },k_{iz},\overline{r}\right)-E_{vi}Rg{\overline{N}}_n\left(k_{i\rho },k_{iz},\overline{r}\right)\right] \end{gathered}$$

(13)

The exciting field coefficients are

$$a_{n^{\prime }}^{EM}\left(k_z^{\prime }\right)=\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{hi}\delta \left(k_z^{\prime }-k_{iz}\right)$$

(14)

$$a_{n^{\prime }}^{EN}\left(k_z^{\prime }\right)=\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{vi}\delta \left(k_z^{\prime }-k_{iz}\right)$$

(15)

where $\delta \left(k_z^{\prime }-k_{iz}\right)$ is the delta function indicating that $k_z^{\prime }=k_{iz}.$ Substitute (14) and (15) into Equations (7) and (8), and we obtain the scattered field coefficients as

$$\begin{gathered} a_n^{SM}\left(k\cos \theta \right)={\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(M,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{hi}\right) \\ +{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(M,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(-\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{vi}\right) \end{gathered}$$

(16)

$$\begin{gathered} a_n^{SN}\left(k\cos \theta \right)={\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(N,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{hi}\right) \\ +{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(N,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(-\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{vi}\right) \end{gathered}$$

(17)

The far field scattered field is

$$\begin{gathered} \begin{array}{c}{\overline{E}}^s\left(\overline{r}\right)\hfill \\ \hfill =-i\widehat{\varphi }k\sin \theta \frac{2}{i}\frac{exp\left(ikr\right)}{r}{\displaystyle \sum }_{n}exp(-i\frac{n\pi }{2})exp\left(in\varphi \right)[{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(M,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{hi}\right)\hfill \\ \hfill \\ +{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(M,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)(-\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{vi})]\hfill \\ \hfill \\ -\widehat{\theta }k\sin \theta \frac{2}{i}\frac{exp\left(ikr\right)}{r}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)[{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(N,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{hi}\right)\hfill \\ \hfill \\ +{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(N,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)(-\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}{E}_{vi})]\hfill \end{array} \end{gathered}$$

(18)

The far field scattered field is written as

$${\overline{E}}^s\left(\overline{r}\right)=\frac{exp\left(ikr\right)}{r}\left({\widehat{v}}_s{E}_{vs}+{\widehat{h}}_s{E}_{hs}\right)$$

(19)

where ${\widehat{v}}_s={\widehat{\theta }}_s$, ${\widehat{h}}_s={\widehat{\varphi }}_s$. The scattered field components $E_{vs}$ and $E_{hs}$ have dimensions of length and are in terms of scattering amplitudes, $f_{vv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$, $f_{vh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$, $f_{hv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$, and $f_{hh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right).$

$$E_{vs}=f_{vv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)E_{vi}+f_{vh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)E_{hi}$$

(20)

$$E_{hs}=f_{hv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)E_{vi}+f_{hh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)E_{hi}$$

(21)

Equations (18)–(21), we obtain

$$f_{vv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)=-k\sin \theta \frac{2}{i}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)\left[{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(N,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)(-\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }})\right]$$

(22)

$$f_{vh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)=-k\sin \theta \frac{2}{i}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)\left[{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(N,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}\right)\right]$$

(23)

$$f_{hv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)=-ik\sin \theta \frac{2}{i}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)\left[{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(M,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)(-\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }})\right]$$

(24)

$$f_{hh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)=-ik\sin \theta \frac{2}{i}{\displaystyle \sum }_{n}exp\left(-i\frac{n\pi }{2}\right)exp\left(in\varphi \right)\left[{\displaystyle \sum }_{n^{\prime }}{T}_{n{n}^{\prime }}^{\left(M,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)\left(\frac{i^{n^{\prime }+1}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}\right)\right]$$

(25)

The above expressions indicate that the T matrix coefficients, $T_{nn\prime }$ are Fourier coefficients of the double Fourier series of the far-field scattering amplitudes. Using Fourier coefficients integration, we have

$$T_{n{n}^{\prime }}^{\left(N,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)=\frac{i\sin {\theta }_i}{8{\pi }^2\sin \theta }{i}^{n-n^{\prime }}{{\displaystyle \int }}_0^{2\pi }d\varphi {{\displaystyle \int }}_0^{2\pi }d{\varphi }_{i}exp\left(-in\varphi +i{n}^{\prime }{\varphi }_i\right)f_{vv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$$

(26)

$$T_{n{n}^{\prime }}^{\left(N,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)=\frac{-\sin {\theta }_i}{8{\pi }^2\sin \theta }{i}^{n-n^{\prime }}{{\displaystyle \int }}_0^{2\pi }d\varphi {{\displaystyle \int }}_0^{2\pi }d{\varphi }_{i}exp\left(-in\varphi +i{n}^{\prime }{\varphi }_i\right)f_{vh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$$

(27)

$$T_{n{n}^{\prime }}^{\left(M,N\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)=\frac{\sin {\theta }_i}{8{\pi }^2\sin \theta }{i}^{n-n^{\prime }}{{\displaystyle \int }}_0^{2\pi }d\varphi {{\displaystyle \int }}_0^{2\pi }d{\varphi }_{i}exp\left(-in\varphi +i{n}^{\prime }{\varphi }_i\right)f_{hv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$$

(28)

$$T_{n{n}^{\prime }}^{\left(M,M\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)=\frac{i\sin {\theta }_i}{8{\pi }^2\sin \theta }{i}^{n-n^{\prime }}{{\displaystyle \int }}_0^{2\pi }d\varphi {{\displaystyle \int }}_0^{2\pi }d{\varphi }_{i}exp\left(-in\varphi +i{n}^{\prime }{\varphi }_i\right)f_{hh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$$

(29)

The far-field scattering amplitudes for incident plane waves are computed in commercial software such as FEKO and HFSS.

2.2. Wave Multiple Scattering Theory (W-MST)

Let there be N scatterers. All the scatterers are within enclosing cylinders. The T matrices are used to describe the scattered fields outside the circumscribing cylinders at all distance ranges. Consider an incident wave with incident field ${\overline{E}}^{inc}\left(\overline{r}\right)$. For W-MST, the final exciting field on one scatterer, ′$q$′, is the sum of the incident field and the scattered fields from all other scatterers except itself. The Foldy–Lax wave multiple scattering equations are

$${\overline{E}}_q^{ex}\left(r\right)={\overline{E}}^{inc}\left(\overline{r}\right)+{\displaystyle \sum }_{p=1,p\ne q}^N{\overline{E}}_{qp}^s\left(\overline{r}\right)$$

(30)

where $q=1,2,\dots ,N$. In (30) ${\overline{E}}_q^{ex}\left(r\right)$ is the “final” exciting field of scatterer $q$, and ${\overline{E}}_{qp}^s\left(\overline{r}\right)$ is the “final” scattered field from scatterer $p$ to scatterer $q$. In (30), the summation of $p$ is over all scatterers, $p=1,2,\dots ,N$, with $p\ne q$. Using cylindrical coordinate system, $\overline{r}=\overline{\rho }+z\widehat{z}$, the center location of the $qth$ scatterer is ${\overline{r}}_q={\overline{\rho }}_q+z_q\widehat{z}=x_q\widehat{x}+y_q\widehat{y}+z_q\widehat{z}$. Since the scatterer is inside a circumscribing cylinder, the ${\overline{\rho }}_q$ coordinate coincides with the horizontal center of the cylinder. Since each scatterer is enclosed by a cylindrical surface, the final excitation fields can be expanded using the incoming VCW. For scatterer $q$ centered at ${\overline{r}}_q$, ${\overline{E}}_q^{ex}\left(r\right)$ is expanded as

$${\overline{E}}_q^{ex}\left(\overline{r}\right)=\sum _m\int d{k}_z\left[Rg{\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)a_{m,q}^{EM}\left(k_z\right)+Rg\overline{N}\left(k_z,{\overline{rr}}_q\right)a_{m,q}^{EN}\left(k_z\right)\right]$$

(31)

In (31), ${\overline{rr}}_q$ =$\overline{r}-{\overline{r}}_q$, $Rg{\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)$ and $Rg\overline{N}\left(k_z,{\overline{rr}}_q\right)$ are VCW-centered at ${\overline{r}}_q$. The incident plane wave is in direction ${\overline{k}}_i=k_{i\rho }\cos {\varphi }_i\widehat{x}+k_{i\rho }\sin {\varphi }_i\widehat{y}+k_{iz}\widehat{z}$ with $k_{iz}=k\cos {\theta }_i$ and $k_{i\rho }=k\sin {\theta }_i$. The electric field of the incident plane wave is

$${\overline{E}}^{inc}(\overline{r})=(\hat{v}{E}_{vi}+\hat{h}{E}_{hi})exp(i{\overline{k}}_i·\overline{r})$$

(32)

To balance the Foldy–Lax Equation (30), all terms need to have VCW with the center at the $q$th scatterer, ${\overline{r}}_q$. To achieve that, we make use of wave transformation. For the incident wave, we set $exp\left(i{\overline{k}}_i\cdot \overline{r}\right)=$ $exp\left(i{\overline{k}}_i\cdot {\overline{r}}_q\right)exp\left(i{\overline{k}}_i\cdot \left(\overline{r}-{\overline{r}}_q\right)\right).$ A phase shift results when we use the VCW expansion in Equation (13) with center at ${\overline{r}}_q,$

$$\begin{array}{cc}\hfill {\overline{E}}_q^{inc}\left(\overline{r}\right)=& \int d{k}_{z}exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\sum _m\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}\hfill \\ \hfill & \times \left[i{E}_{hi}Rg{\overline{M}}_m\left(k_{\rho },k_z,{\overline{rr}}_q\right)-E_{vi}Rg{\overline{N}}_m\left(k_{\rho },k_z,{\overline{rr}}_q\right)\right]\delta \left(k_z-k_{iz}\right)\hfill \end{array}$$

(33)

The Dirac delta function $\delta \left(k_z-k_{iz}\right)$ indicates that the incident wave is with $k_z=k_{iz}$ in the integration $\int d{k}_z$. Since the scattered field from $p$ to $q$ originates from scatterer $p$, the outgoing VCW from scatterer $p$ are initially expressed with origin at ${\overline{r}}_p$. The expansion is

$${\overline{E}}_{qp}^s\left(\overline{r}\right)={\displaystyle \sum }_n{\int }_{-\infty }^{\infty }d{k}_z\left[{\overline{M}}_n\left(k_z,{\overline{rr}}_p\right)a_{n,p}^{SM}\left(k_z\right)+{\overline{N}}_n\left(k_z,{\overline{rr}}_p\right)a_{n,p}^{SN}\left(k_z\right)\right]$$

(34)

To translate ${\overline{E}}_{qp}^s\left(\overline{r}\right)$ to origin, ${\overline{r}}_q$ of $q$ scatterer, we use the translational addition theorem. Let ${\overline{r}}_p={\overline{\rho }}_p+z_p\widehat{z}=x_p\widehat{x}+y_p\widehat{y}+z_p\widehat{z}$, $\overline{{\rho }_p{\rho }_q}={\overline{\rho }}_p-{\overline{\rho }}_q$. The quantity ${\varphi }_{\overline{{\rho }_p{\rho }_q}}=\mathrm{t}\mathrm{a}{\mathrm{n}}^{-1}\frac{y_p-y_q}{x_p-x_q}$ is the angle that ${\overline{\rho }}_p-{\overline{\rho }}_q$ makes with the $x$ axis. The translational addition theorems are (Tsang et al. volume 2)

$${\overline{M}}_n\left(k_z,{\overline{rr}}_p\right)={\displaystyle \sum }_{m}Rg{\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)H_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i\left(m-n\right){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)$$

(35)

$${\overline{N}}_n\left(k_z,{\overline{rr}}_p\right)={\displaystyle \sum }_{m}Rg{N}_m\left(k_z,{\overline{rr}}_q\right)H_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i\left(m-n\right){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)$$

(36)

Then we have

$$\begin{gathered} {\overline{E}}_{qp}^s\left(\overline{r}\right)=\sum _n{\int }_{-\infty }^{\infty }d{k}_z\sum _{m}Rg{\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)H_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i\left(m-n\right){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)a_{n,p}^{SM}\left(k_z\right) \\ \begin{array}{c}+{\displaystyle \sum _n}{\int }_{-\infty }^{\infty }d{k}_z{\displaystyle \sum _m}Rg{\overline{N}}_m\left(k_z,{\overline{rr}}_q\right)H_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i\left(m-n\right){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)a_{n,p}^{SN}\left(k_z\right)\end{array} \end{gathered}$$

(37)

Since all terms now have the dependence Rg${\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)$ or Rg${\overline{N}}_m\left(k_z,{\overline{rr}}_q\right)$, we balance the coefficients of each term. Balancing the terms with Rg${\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)$, we have

$$\begin{array}{c}{a}_{m,q}^{EM}\left(k_z\right)=exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}i{E}_{hi}\delta \left(k_z-k_{iz}\right)+{\displaystyle \sum _{p\ne q}}{\displaystyle \sum _n}{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)\\ \times exp\left(-i(m-n){\varphi }_{\overline{{\rho }_{ppq}}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)a_{n,p}^{SM}\left(k_z\right)\end{array}$$

(38)

Balancing Rg${\overline{N}}_m\left(k_z,{\overline{rr}}_q\right)$, we have

$$\begin{array}{c}{a}_{m,q}^{EN}\left(k_z\right)=exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}\left[-E_{vi}\right]\delta \left(k_z-k_{iz}\right)\\ +{\displaystyle \sum _{p\ne q}}{\displaystyle \sum _n}{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)\\ \left.\times exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)a_{n,p}^{SN}\left(k_z\right)\right)\end{array}$$

(39)

The multiple scattering equations in (38) and (39) relate exciting field coefficients $a_{m,q}^{EM}\left(k_z\right),a_{m,q}^{EN}\left(k_z\right)$ and scattering field coefficients $a_{n,p}^{SM}\left(k_z\right)$, $a_{n,p}^{SN}\left(k_z\right)$ among the scatterers. To obtain Foldy–Lax equations with the exciting field coefficients $a_{m,q}^{EM}\left(k_z\right)$, $a_{m,q}^{EN}\left(k_z\right)$ as unknowns, we use the T matrix relations

$$a_{n,p}^{SM}\left(k_z\right)=\sum _{n^{\prime }}\int d{k}_z{}^{\prime }\left[T_{n{n}^{\prime },p}^{(M,M)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EM}\left(k_z{}^{\prime }\right)+T_{n{n}^{\prime },p}^{(M,N)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EN}\left(k_z{}^{\prime }\right)\right]$$

(40)

$$a_{n,p}^{SN}\left(k_z\right)=\sum _{n^{\prime }}\int d{k}_z{}^{\prime }\left[T_{n{n}^{\prime },p}^{(N,M)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EM}\left(k_z{}^{\prime }\right)+T_{n{n}^{\prime },p}^{(N,N)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EN}\left(k_z{}^{\prime }\right)\right]$$

(41)

Then, the Foldy–Lax MST equations with exciting field coefficients as field unknowns are

$$\begin{array}{cc}\hfill & a_{m,q}^{EM}\left(k_z\right)=exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}i{E}_{hi}\delta \left(k_z-k_{iz}\right)\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \times \sum _{n^{\prime }}\int d{k}_z{}^{\prime }\left[T_{n{n}^{\prime },p}^{(M,M)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EM}\left(k_z{}^{\prime }\right)+T_{n{n}^{\prime },p}^{(M,N)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EN}\left(k_z{}^{\prime }\right)\right]\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill a_{m,q}^{EN}\left(k_z\right)=& exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}\left[-E_{vi}\right]\delta \left(k_z-k_{iz}\right)\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \times \sum _{n^{\prime }}\int d{k}_z{}^{\prime }\left[T_{n{n}^{\prime },p}^{(N,M)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EM}\left(k_z{}^{\prime }\right)+T_{n{n}^{\prime },p}^{(N,N)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EN}\left(k_z{}^{\prime }\right)\right]\hfill \end{array}$$

(43)

In solving the Foldy–Lax equations numerically, it is sometimes convenient to use scattered field coefficients as field unknowns. To obtain scattered field coefficients, MST, multiply (38) by $\sum _m\int d{k}_z{T}_{m^{\prime }m,q}^{\left(M,M\right)}\left({k_z}^{″},k_z\right)$ and multiply (39) by $\sum _m\int d{k}_z{T}_{m^{\prime }m,q}^{\left(M,N\right)}\left({k_z}^{″},k_z\right)$. Adding the two equations gives

$$\begin{array}{c}\hfill a_{m^{\prime },q}^{SM}\left(k_z^{″}\right)=\sum _{m}exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}\left[T_{m^{\prime }m,q}^{(M,M)}\left(k_z^{″},k_{iz}\right)i{E}_{hi}-T_{m^{\prime }m,q}^{(M,N)}\left(k_z^{″},k_{iz}\right)\left(E_{vi}\right)\right]\\ \hfill +\sum _m\int d{k}_z\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\\ \hfill \left[T_{m^{\prime }m,q}^{(M,M)}\left(k_z^{″},k_z\right)a_{n,p}^{SM}\left(k_z\right)+T_{m^{\prime }m,q}^{(M,N)}\left(k_z^{″},k_z\right)a_{n,p}^{SM}\left(k_z\right)\right]\end{array}$$

(44)

Next, multiply (38) by $\sum _m\int d{k}_z{T}_{m^{\prime }m,q}^{\left(N,M\right)}\left({k_z}^{″},k_z\right)$ and multiply (39) by $\sum _m\int d{k}_z{T}_{m^{\prime }m,q}^{\left(N,N\right)}\left({k_z}^{″},k_z\right)$. Add the two equations gives

$$\begin{array}{c}\hfill a_{m^{\prime },q}^{SN}\left(k_z^{″}\right)=\sum _{m}exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}\left[T_{m^{\prime }m,q}^{(N,M)}\left(k_z^{\prime \prime },k_{iz}\right)i{E}_{hi}-T_{m^{\prime }{m}_{,q}}^{(N,N)}\left(k_z^{″},k_{iz}\right)\left(E_{vi}\right)\right]\\ \hfill +\sum _m\int d{k}_z\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\\ \hfill \left[T_{m^{\prime }m,q}^{(N,M)}\left(k_z^{″},k_z\right)a_{n,p}^{SM}\left(k_z\right)+T_{m^{\prime }m,q}^{(N,N)}\left(k_z^{″},k_z\right)a_{n,p}^{SM}\left(k_z\right)\right]\end{array}$$

(45)

Equations (44) and (45) are the Foldy–Lax W-MST equations with scattered field coefficients as field unknowns. They are also called the self-consistent MST equations. We can also have exciting field Foldy–Lax MST equations by separating out the incident field Dirac delta functions in (42) and (43). Let

$$a_{m,q}^{EM}\left(k_z\right)=exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}i{E}_{hi}\delta \left(k_z-k_{iz}\right)+a_{m,q}^{EOM}\left(k_z\right)$$

(46)

$$a_{m,q}^{EN}\left(k_z\right)=exp\left(i{\overline{k}}_i·{\overline{r}}_q\right)\frac{i^{m}exp\left(-im{\varphi }_i\right)}{k_{i\rho }}\left[-E_{vi}\right]\delta \left(k_z-k_{iz}\right)+a_{m,q}^{EON}\left(k_z\right)$$

(47)

where “O” in the superscript $a_{m,q}^{EOM}\left(k_z\right)$ and $a_{m,q}^{EON}\left(k_z\right)$ stands for from “other” scatterers. Then the integral equations for $a_{m,q}^{EOM}\left(k_z\right)$ and $a_{m,q}^{EON}\left(k_z\right)$ are

$$\begin{array}{cc}\hfill & a_{m,q}^{EOM}\left(k_z\right)=\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}{T}_{n{n}^{\prime },p}^{(M,M)}\left(k_z,k_{iz}\right)exp\left(i{\overline{k}}_i·{\overline{r}}_p\right)\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}i{E}_{hi}\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}\int d{k}_z{}^{\prime }{T}_{n{n}^{\prime },p}^{(M,M)}\left(k_z,k_z{}^{\prime }\right)a_{n^{\prime },p}^{EOM}\left({k_z}^{\prime }\right)\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}{T}_{n{n}^{\prime },p}^{(M,N)}\left(k_z,k_{iz}\right)exp\left(i{\overline{k}}_i·{\overline{r}}_p\right)\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}\left[-E_{vi}\right]\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}\int d{k_z}^{\prime }{T}_{n{n}^{\prime },p}^{(M,N)}\left(k_z,{k_z}^{\prime }\right)a_{n^{\prime },p}^{EON}\left({k_z}^{\prime }\right)\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill a_{m,q}^{EON}\left(k_z\right)=& \sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}{T}_{n{n}^{\prime },p}^{(N,M)}\left(k_z,k_{iz}\right)exp\left(i{\overline{k}}_i·{\overline{r}}_p\right)\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}i{E}_{hi}\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}\int d{k}_z{}^{\prime }{T}_{n{n}^{\prime },p}^{(N,M)}\left(k_z,{k_z}^{\prime }\right)a_{n^{\prime },p}^{EOM}\left({k_z}^{\prime }\right)\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|\overline{{\rho }_p{\rho }_q}\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}{T}_{n{n}^{\prime },p}^{(N,N)}\left(k_z,k_{iz}\right)exp\left(i{\overline{k}}_i·{\overline{p}}_p\right)\frac{i^{n^{\prime }}exp\left(-i{n}^{\prime }{\varphi }_i\right)}{k_{i\rho }}\left[-E_{vi}\right]\hfill \\ \hfill & +\sum _{p\ne q}\sum _n{H}_{m-n}^{\left(1\right)}\left(k_{\rho }\left|{\overline{\rho }}_p{\rho }_q\right|\right)exp\left(-i(m-n){\varphi }_{\overline{{\rho }_p{\rho }_q}}\right)exp\left(-i{k}_z\left(z_p-z_q\right)\right)\hfill \\ \hfill & \sum _{n^{\prime }}\int d{k_z}^{\prime }{T}_{n{n}^{\prime },p}^{(N,N)}\left(k_z,{k_z}^{\prime }\right)a_{n^{\prime },p}^{EON}\left({k_z}^{\prime }\right)\hfill \end{array}$$

(49)

It should be noted that the Foldy–Lax equations are entirely in terms of the wave interactions outside the circumscribing cylinders. Thus, the wave interactions are completely described by the T matrices of the scatterers.

2.3. Final Fields

After solving the Foldy–Lax MST equations, the final fields both outside and inside the circumscribing cylinders can be calculated by the following procedures.

2.3.1. Outside the Enclosing Cylinders

We obtained the scattered field coefficients from solutions of the Foldy–Lax W-MST equations. Outside the enclosing cylinder, the electric field is equal to the sum of the incident fields and the scattered fields from all the scatterers

$$\overline{E}\left(\overline{r}\right)={\overline{E}}^{inc}\left(\overline{r}\right)+{\displaystyle \sum }_{q=1}^N{\overline{E}}_q^s\left(\overline{r}\right)$$

(50)

The final scattered fields are calculated by

$${\overline{E}}_q^s\left(\overline{r}\right)={\displaystyle \sum }_n\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{k}_z\left[{\overline{M}}_n\left(k_z,{\overline{rr}}_q\right)a_{n,q}^{SM}\left(k_z\right)+{\overline{N}}_n\left(k_z,{\overline{rr}}_q\right)a_{n,q}^{SN}\left(k_z\right)\right]$$

(51)

2.3.2. Inside the Enclosing Cylinders

By solving the Foldy–Lax MST equations, we obtained the final exciting field coefficients $a_{n,q}^{EM}\left(k_z\right)$, $a_{n,q}^{EN}\left(k_z\right)$, $q=1,2,\cdots ,N$.

Since commercial software were used to solve the scattering problem of a single object for plane wave excitations, we express the exciting field as plane waves by using transformations of cylindrical waves in terms of plane waves. The final exciting field of the $q$-th scatterer is

$${\overline{E}}_q^{ex}\left(\overline{r}\right)={\displaystyle \sum }_m\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{k}_z\left[Rg{\overline{M}}_m\left(k_z,{\overline{rr}}_q\right)a_{m,q}^{EM}\left(k_z\right)+Rg{\overline{N}}_m\left(k_z,{\overline{rr}}_q\right)a_{m,q}^{EN}\left(k_z\right)\right]$$

(52)

From Tsang et al. volume 1, let $k_x=k_{\rho }\cos {\varphi }_k$, $k_y=k_{\rho }\sin {\varphi }_k$,

$$exp\left(i{k}_{x}x+i{k}_{y}y\right)=exp\left(i{k}_{\rho }\rho \cos \left(\varphi -{\varphi }_k\right)\right)={\displaystyle \sum }_n{J}_n\left(k_{\rho }\rho \right)exp\left(in\left(\varphi -{\varphi }_k\right)+in\frac{\pi }{2}\right)$$

(53)

Multiply both sides by ${\int }_0^{2\pi }d{\varphi }_{k}exp\left(i{n}^{\prime }{\varphi }_k\right)$ and integrate

$$J_n\left(k_{\rho }\rho \right)exp\left(in\varphi \right)=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }_{k}exp\left(in\left({\varphi }_k-\frac{\pi }{2}\right)\right)exp\left(i{k}_{x}x+i{k}_{y}y\right)$$

(54)

The scalar cylindrical waves

$$Rg{\psi }_n\left(k_z,\overline{r}\right)=J_n\left(k_{\rho }\rho \right)exp\left(in\varphi \right)exp\left(i{k}_{z}z\right)$$

(55)

Then

$$Rg{\psi }_n\left(k_z,\overline{r}\right)=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }_{k}exp\left(in\left({\varphi }_k-\frac{\pi }{2}\right)\right)exp\left(i{k}_{x}x+i{k}_{y}y+i{k}_{z}z\right)$$

(56)

Equation (56) expresses scalar cylindrical waves in terms of scalar plane waves. For VCW, we then have the transformation to vector plane waves

$$Rg{\overline{M}}_n\left(k_z,\overline{r}\right)=\frac{\left(-ik\sin {\theta }_k\right)}{2\pi }{\int }_0^{2\pi }d{\varphi }_{k}exp\left(in\left({\varphi }_k-\frac{\pi }{2}\right)\right)\left(\widehat{h}\left({\theta }_k,{\varphi }_k\right)exp\left(i{k}_{x}x+i{k}_{y}y+i{k}_{z}z\right)\right)$$

(57)

$$Rg{\overline{N}}_n\left(k_z,\overline{r}\right)=\frac{\left(-k\sin {\theta }_k\right)}{2\pi }{\int }_0^{2\pi }d{\varphi }_{k}exp\left(in\left({\varphi }_k-\frac{\pi }{2}\right)\right)\left(\widehat{v}\left({\theta }_k,{\varphi }_k\right)exp\left(i{k}_{x}x+i{k}_{y}y+i{k}_{z}z\right)\right)$$

(58)

where $k_z=k\cos {\theta }_k$, $k_{\rho }=k\sin {\theta }_k$ and $\widehat{v}\left({\theta }_k,{\varphi }_k\right)$ and $\widehat{h}\left({\theta }_k,{\varphi }_k\right)$ are vertical and horizontal polarizations

$$\widehat{v}\left({\theta }_k,{\varphi }_k\right)=\cos {\theta }_k\cos {\varphi }_k\widehat{x}+\cos {\theta }_k\sin {\varphi }_k\widehat{y}-\sin {\theta }_k\widehat{z}$$

(59)

$$\widehat{h}\left({\theta }_k,{\varphi }_k\right)=-\sin {\varphi }_k\widehat{x}+\cos {\varphi }_k\widehat{y}$$

(60)

Then

$$\begin{gathered} {\overline{E}}_q^{ex}\left(\overline{r}\right)=\frac{-ik}{2\pi }\sum _n\int d{k}_z{a}_{n,q}^{EM}\left(k_z\right)\sin {\theta }_k{\int }_0^{2\pi }d{\varphi }_{k}exp\left(in\left({\varphi }_k-\frac{\pi }{2}\right)\right)\widehat{h}\left({\theta }_k,{\varphi }_k\right)exp\left(i\overline{k}\cdot \left(\overline{r}-{\overline{r}}_q\right)\right) \\ \begin{array}{c}-\frac{k}{2\pi }{\displaystyle \sum _n}\int d{k}_z{a}_{n,q}^{EN}\left(k_z\right)\sin {\theta }_k{\int }_0^{2\pi }d{\varphi }_{k}exp\left(in\left({\varphi }_k-\frac{\pi }{2}\right)\right)\widehat{v}\left({\theta }_k,{\varphi }_k\right)exp\left(i\overline{k}\cdot \left(\overline{r}-{\overline{r}}_q\right)\right)\end{array} \end{gathered}$$

(61)

The above provides a spectrum of plane waves exciting (incident on) scatterer $q$. Using commercial software, the fields inside the enclosing cylinder $q$ can be calculated.

2.4. Rotation and Efficient Use of Re-usable T Matrices

In modeling vegetation and forests, we model a plant or a tree as a single object. We use commercial software to calculate the scattering of a plant or a tree. The three common methods of commercial software are the Finite Element Method (FEM), Finite Difference Time Domain method (FDTD), and the Method of Moment (MoM). Initially, we used the commercial software HFSS (3D Hugh Frequency Simulation Software; version R15; Ansys, Canonsburg, PA, USA), which is based on FEM. Recently, we also used FEKO (Feldberechnung für körper mit beliebiger oberfläche), which can be translated as “field calculations involving bodies of arbitrary shape”. FEKO is MoM based and accelerated with Multilevel Fast Multipole (MLFMA). Using either HFSS or FEKO, we used incident plane waves and calculate the far field scattering amplitudes $f_{vv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$, $f_{vh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$, $f_{hv}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$ and $f_{hh}\left(\theta ,\varphi ;{\theta }_i,{\varphi }_i\right)$. The options of far-field scattering amplitudes are available in commercial software. Using the above transformation formulates, we obtain the T matrix coefficients of vector cylindrical waves. In the above, we neglected evanescent waves. To include evanescent waves, field solutions can be obtained in the surface of the enclosing cylinder. Then, Huygens’s principles can be used to accurately extract the T matrix coefficients. We have adopted the Huygens principles previously [3,4,5].

The T matrices are reusable and can be used to generate the T matrices of sech as 100 plants in a vegetation field. In traditional modeling of vegetations, the branches and leaves obey orientation distribution as described by the Euler angles. Let the axis of symmetry be

$${\widehat{z}}_b=\sin \beta \cos \alpha \widehat{x}+\sin \beta \sin \alpha \widehat{y}+\cos \beta \widehat{z}$$

(62)

where $\beta$ and $\alpha$ are the orientation angles of rotation. The probability distributions are $p_{\beta }\left(\beta \right)$ and $p_{\alpha }\left(\alpha \right)$ for $\beta$ and $\alpha$ respectively such that

$${\int }_0^{\pi }d\beta p_{\beta }\left(\beta \right){\int }_0^{2\pi }d\alpha p_{\alpha }\left(\alpha \right)=1$$

(63)

In the past, $p_{\beta }\left(\beta \right)$ and $p_{\alpha }\left(\alpha \right)$ are prescribed for branches and leaves. The probability density function is usually assumed to be uniform between 0 and $2\pi$.

We can generate a large number such as 100 different plants in the following manner. We first generate 5 plants that obey $p_{\beta }\left(\beta \right)$ that are prescribed for branches and leaves. The T matrices of these 5 plants are calculated by the procedure earlier. In the hybrid method, the plant is a single scatterer. Thus, the plant has a prescribed $p_{\alpha }\left(\alpha \right)$. Then the 5 plants are rotated by $\alpha$ such that $p_{\alpha }\left(\alpha \right)$ is obeyed. For the entire plant as a single scatterer ${\widehat{z}}_b=\widehat{z}$. For a rotation by $\alpha$, the T matrix follows a simple rotation relation as described below. If we have 20 rotations of $\alpha$, we then have 100 different plants in the field. We can also have 100 rotations of $\alpha$ to have 100 plants. Thus, in this manner, we only need to calculate the 5 T matrices of plants for the 100 plants. The 100 plants can be placed randomly. They can also be placed in a quasi-periodic manner with random displacement from the periodic positions. The T matrix elements are stored and are of multiple uses for Monte Carlo simulation. The computations of T matrix are a set-up step. The 5 T matrices can also be stored for future use for various densities of plants and for many realizations of random displacements of plants.

Below, we derive the relations of a rotation of T matrix for vector cylindrical waves. The body axes of the plant by arbitrary choice are ${\widehat{x}}_b$ and ${\widehat{y}}_b$. So that the relation between body axes and principal axes are as follows

$${\widehat{x}}_b=\cos \alpha \widehat{x}+\sin \alpha \widehat{y}$$

(64)

$${\widehat{y}}_b=-\sin \alpha \widehat{x}+\cos \alpha \widehat{y}$$

(65)

We calculate the T matrix coefficients in the body frame by using the far field scattering amplitude in the body frame. For example,

$$T_{n{n}^{\prime }}^{\left(N,N\right)\left(b\right)}\left(k\cos \theta ,k\cos {\theta }_i\right)=\frac{i\sin {\theta }_i}{8{\pi }^2\sin \theta }{i}^{n-n^{\prime }}{\int }_0^{2\pi }d{\varphi }_i^{\left(b\right)}{\int }_0^{2\pi }d{\varphi }_{i}exp\left(-in{\varphi }^{\left(b\right)}+i{n}^{\prime }{\varphi }_i^{\left(b\right)}\right)f_{vv}^{\left(b\right)}\left({\theta }_s,{\varphi }_s^{\left(b\right)};{\theta }_i,{\varphi }_i^{\left(b\right)}\right)$$

(66)

where superscript ($b$) stands for body frame. Next, suppose the body frame is rotated by an angle $\alpha$ about the vertical $z$ axis.

Then the scattering amplitude in the principal frame f_vv and the body frame $f_{vv}^{\left(b\right)}$ are related by the relationship

$$f_{vv}\left({\theta }_s,{\varphi }_s;{\theta }_i,{\varphi }_i\right)=f_{vv}^{\left(b\right)}\left({\theta }_s,{\varphi }_s^{\left(b\right)};{\theta }_i,{\varphi }_i^{\left(b\right)}\right)=f_{vv}^{\left(b\right)}\left({\theta }_s,{\varphi }_s-\alpha ;{\theta }_i,{\varphi }_i-\alpha \right)$$

(67)

Substituting in, we obtain the following relations for the T matrices

$$T_{n{n}^{\prime }}^{\left(j,l\right)}\left(k\cos {\theta }_s,k\cos {\theta }_i\right)=exp\left(i\left(n^{\prime }-n\right)\alpha \right)T_{n{n}^{\prime }}^{\left(j,l\right)\left(b\right)}\left(k\cos {\theta }_s,k\cos {\theta }_i\right)$$

(68)

where $j,l=M,N$.

2.5. Computational Efficiency of the Hybrid Method for Statistical Moments of Fields

The hybrid method is adapted to solving wave MST for a large number of scatterers of vegetation and forests.

In Figure 2, we have 2 plants with different $\beta$. With the 2 plants we rotate 9 different $\alpha$ to have 9 corn plants. Each of the 9 plants has an $\alpha$. We first place them periodically. Next, we displace the corn plants from their periodic positions by a defined amount. That generates a corn field with quasi-periodic arrangement of corn.

We make use of several key elements of such Monte Carlo simulations.

(i) Each plant or each tree is a single scatterer, the leaves are aggregated on a branch, and branches are aggregated on a plant or a tree. Neither branches nor leaves are single scatterers.

(ii) Each plant or each tree is a distinct random scatterer with a pre-computed T matrix. The plant is rotated by an angle $\alpha$ and the corresponding T matrix is obtained. This provides a different plant because the plant does not have azimuthal symmetry.

(iii) The positions of plants or trees are random. They can be placed in a quasi-periodic manner such as in a row structure with random displacements from the periodic positions.

(iv) A field of many distinct plants such as 100 plants can be generated with a few, such as 5, precomputed and stored T matrices.

In random scattering problems, the scattered field solutions have fluctuations of amplitudes and phases. Thus, averaging results over realizations are needed to obtain useful solutions for remote sensing applications.

The computational objective is to calculate the Statistical Moments of Fields in Vegetation Fields and Forests. At present, the primary interests are the first and statistical moments. This objective is distinctly different from CEM, the objective of which is to calculate the “exact” field solution of each realization.

For each realization, we iterate Foldy–Lax multiple scattering equations with the T matrices of the scatterers to obtain the multiple orders of scattering by scatterers. The iterations provide the first-order solution, the second-order solution, the third-order solution, etc. The iteration, unlike traditional methods such as conjugate gradient or biconjugate gradient in which iterations are carried out until the convergence of the field solution, is achieved for each realization. We iterate and then take realization averaging, for example, after the solution is up to second order, we take realization averaging. Then, after the solution is up to the 4th order, we take the realization averaging. This iteration method in Monte Carlo simulations is analogous to the analytical perturbation method in random scattering. For example, in the second-order small perturbation method for rough surfaces, an averaging is taken after the second-order solution. In the fourth-order perturbation method of rough-surface scattering, averaging of realizations is taken over the solution is obtained up to the fourth order. The reasons for taking averages after even order are to ensure energy conservation up to that order.

For the memory requirement, such as for 100 plants, it is only required to store the 5 T matrices for the different ${\beta }^{\prime }s$. The rotation of different ${\alpha }^{\prime }s$ is merely rotating the 5 T matrices. For the translational addition theorem, the order dependence of the Hankel function is on $m-n$, the difference between the mth order and the nth order. Let the order be truncated at $N_{max}$ with $n,m=0,\pm 1,\pm 2,\dots ,N_{max}$. The memory requirement for the translation theorem is $4N_{max}$+1 and not ${\left(2N_{max}+1\right)}^2$. This reduction of memory is significant because the enclosing cylinder is of moderate to large radius, requiring moderate to large values of $N_{max}$. In the matrix times column vector, the dependence on $m-n$ means that FFT can be applied. In

Table 1. Comparison between CEM method and hybrid method.

	CEM Method	Hybrid Method	Comments
Full Wave	Entire problem of Np number of plants	Single plant is an object. T matrix is obtained for the plant	Each plant is an object In RT: each leaf is an object. Each branch is an object
Field Solutions	N = number of unknowns in field solutions, in millions	Multiple scattering of Np plants Np moderate number, e.g., 100	NA
Reusable	Not reusable, solve N field unknowns for each realization	T matrices of a few plants plus azimuthal α rotations, reused for (i) configurations, random, quasi-periodic (ii) realizations (iii) future use	T matrices put on shelves for future use T matrices Portable
Iteration Solution of Each Realization	large number of iterations (e.g., conjugate gradient) for large number of N field unknowns to reach convergence of “exact” field solution	Iterate Foldy–Lax to obtain multiple scattering order solutions, second order, fourth order, tenth order, even orders of solutions	Significant wave iterations within a plant which are included in T matrix of a plant, less wave interactions between plants For vegetation fields and forests, no more than 10 multiple orders of solutions
Averaging over Realizations to Calculate Statistical Moments	Averaging “exact” solutions over Nr realizations	Averaging Nr realizations after second order, fourth order, sixth order, until statistical moments converge	Averaged second order solutions, fourth-order solutions, are analogous to analytical random media theory SPM in rough surfaces and iterative solutions of radiative transfer equation

, we compare the methodology of computational electromagnetic method (CEM) method of full-wave simulations and the hybrid method. In the table, $N_P=$ number of plants = 100, $N=$ number of field unknowns in full-wave simulations, and $N_r=$ number of realizations = 30 to 100.

2.6. Calculations and Validation of T Matrices of A Single Corn Plant Using Commercial Software

We next illustrate the calculations of T matrix of a single corn plant and the validation of the calculations. We first used FEKO to calculate the far-field scattering amplitudes. Figure 3 shows the real and imaginary parts of $f_{vv}$ as a function of ${\varphi }_s$ for ${\theta }_i={140}^o$, ${\varphi }_i=0^o$; ${\theta }_s={140}^o$.

We used the Fourier coefficients equations in (26)–(29) to calculate the T matrix coefficients from the far field scattering amplitudes. We used $N_{max}=6$. In Figure 4, we plot the T matrix coefficients $T_{m{n}^{\prime }}^{NN}\left({\theta }_s,{\theta }_i\right)$ as a function of $n$ for $n^{\prime }=0$.

There are 2 incident polarizations, $v$ and $h$. For each pair of $\left({\theta }_i,{\varphi }_i\right)$, we can have as many values of $\left({\theta }_s,{\varphi }_s\right)$ for $f_{vv}\left({\theta }_s,{\varphi }_s,{\theta }_i,{\varphi }_i\right)$, $f_{hv}\left({\theta }_s,{\varphi }_s,{\theta }_i,{\varphi }_i\right)$, $f_{vh}\left({\theta }_s,{\varphi }_s,{\theta }_i,{\varphi }_i\right)$, $f_{hh}\left({\theta }_s,{\varphi }_s,{\theta }_i,{\varphi }_i\right)$. Using the far-field scattering amplitudes, we calculated the $T_{n{n}^{\prime }}^{MM}\left({\theta }_s,{\theta }_i\right)$, $T_{n{n}^{\prime }}^{MN}\left({\theta }_s,{\theta }_i\right)$, $T_{n{n}^{\prime }}^{NM}\left({\theta }_s,{\theta }_i\right)$, and $T_{n{n}^{\prime }}^{NN}\left({\theta }_s,{\theta }_i\right)$. The choice of $N_{max}$ depends on ${\theta }_s$ and ${\theta }_i$. The dimensions of T matrix are in

Table 2. The choice of $N_{max}$ for different $\theta$ or ${\theta }_i$.

$\mathit{\theta }$$\mathbf{or} {\mathit{\theta }}_{\mathit{i}}$	${\mathit{N}}_{\mathit{m}\mathit{a}\mathit{x}}$	$2{\mathit{N}}_{\mathit{m}\mathit{a}\mathit{x}}+1$	$\mathit{\theta }$$\mathbf{or} {\mathit{\theta }}_{\mathit{i}}$	${\mathit{N}}_{\mathit{m}\mathit{a}\mathit{x}}$	$2{\mathit{N}}_{\mathit{m}\mathit{a}\mathit{x}}+1$
2°	2	5	95.18°	6	13
7.18°	2	5	100.35°	5	11
12.35°	2	5	105.53°	5	11
17.53°	2	5	110.71°	5	11
22.71°	3	7	115.88°	5	11
27.88°	4	9	121.06°	4	9
33.06°	4	9	126.24°	4	9
38.24°	4	9	131.41°	4	9
43.41°	4	9	136.59°	4	9
48.59°	4	9	141.76°	4	9
53.76°	4	9	146.94°	4	9
58.94°	4	9	152.12°	4	9
64.12°	5	11	157.29°	3	7
69.29°	5	11	162.47°	2	5
74.47°	5	11	167.65°	2	5
79.65°	5	11	172.82°	2	5
84.82°	6	13	178°	2	5
90°	6	13

. The total dimension is $614\times 614$.

Next, we validated the accuracy and the use of T matrix coefficients of calculating the scattered field. We use the T matrix of dimension $614\times 614$ and the scattered field Equations of (7) and (8) with the incident field coefficients of Equations (14) and (15). We calculate the scattered field as a function of a position for the positions as shown in Figure 5a. We calculated the scattered field in 2 ways: using FEKO directly and using the expression of scattered field in terms of T matrices. In Figure 5b, we show the comparisons. The results were in very good agreement, validating the accuracies and the use of T matrix $T_{n{n}^{\prime }}^{NN}$ from the far-field calculated from Equation (29).

2.7. Numerical Results of Hybrid Method of Vegetation Field and Forests

Results of the Hybrid Method can be found in found in refs. [2,4,5,34]. Below we illustrate two examples.

Simulations were performed for the transmission through a simulated forest [2] consisting of 196 cylinders of 20 m height and 12 cm diameter arranged as shown in Figure 6. The results are tabulated in

Table 3. Transmission coefficient from RTE based on distorted Born approximation (RTE/DBA) and the hybrid method from Figure 6. The table is adapted with permission from Ref. [2]. © 2019 IEEE.

	RTE/DBA	Hybrid Method
Transmission	0.35	0.66

. The results showed that the transmission is 1.89 times that of the results of RTE.

For the transmission through vegetation and forests, there have been papers using C-band Sentinel 1 to study soil moisture. Such studies coincide with the results in this paper, i.e., the C-band penetration through vegetation is much larger than predicted by RTE. For forests, the SMOS has the retrieval algorithm applied to forests meaning that the L-band can penetrate through forests with VWC much larger than 5 kg per square meter. SMAP has been conducting a ground and airborne campaign to verify the penetration through forests at L-Band. Thus, although experimental comparison is lacking, there is support from the community that there can be much larger penetration through vegetation fields and forests from the L-band to the C-band.

To consider frequency dependence, we next show an example of the transmission through a field consisting of 196 wheat plants (Figure 7a) at L-, S-, and C-bands [4,5] as a function of VWC. Results of Figure 7b are compared with RTE. Firstly, the results showed the transmission of full-wave simulations are much larger than RTE. Secondly, the transmission at C-band was only slightly less than that at S-band, showing the frequency dependence is weak between S-band and C-band. On the other hand, RTE showed a big drop in transmission from the S- to the C-band, indicating a strong increase of attenuation with frequency from the S-band to the C-band predicted by RTE. The full-wave simulation results using the hybrid method had little difference between the S-band and C-band, showing saturation with frequency between the S-band and C-band.

Our current work constitutes (i) developing fast computation techniques for calculating the solutions of Foldy–Lax equations in the VCW representation and (ii) calculation of solutions for the cases of forests with trees up to 20 m tall. For optical scattering work, readers should refer to papers by Ping Yang from Texas A&M University and Karri Muinonen from the University of Helsinki.

3. Signals of Opportunity

3.1. GNSS-R and SoOp Introduction

The operating GNSS-R missions include the Techdemosat-1 (TDS-1) [6], launched by UK in 2014, the Cyclone Global Navigation Satellite System (CYGNSS) [7], launched by NASA in 2016, and Bufeng-1, launched by China in 2019 [35]. The GNSS-R data are collected in the form of Delay Doppler Maps (DDMs), which have been applied to the retrieval of ocean wind speed [36], sea ice thickness [37], and monitoring the wetland changes [38].

For the SoOp CYGNSS mission at the L-band, experimental results and comparisons were presented in the IGARSS 2021 paper by Campbell et al. There has been a ground campaign conducted for CYGNSS since 2018.

There are major differences between the traditional rough-surface bistatic scattering formula and the reflected signals in SoOp. The comparison between GNSS-R and radar back scattering is shown in

Table 4. Comparison between radar backscattering and GNSS-R.

	Radar Backscattering with 40 Degrees Incident Angle	GNSS-R Observation Close to Specular Direction
Scattering	Large angle from specular	Small angle from specular
Kirchhoff integral	Not Valid as Kirchhoff predicts VV is comparable to HH	Accurate near specular direction
Roughness	Microwave roughness Topography have small effects	Topography strong influence +microwave roughness
Mean field intensity/Covariance of field	Covariance of fields only	Mean field intensity and Covariance of fields
Gamma/sigma0	−25 dB to 0 dB	Much Larger values 10 dB to 30 dB

. In the usual rough-surface bistatic scattering, the formula is an extension of radar backscattering [39]. In radar rough-surface backscattering from soil surfaces, contributions to ${\sigma }_0$ come from the microwave roughness which have the ratios of correlation length to rms height ratio around 10. In extension to bistatic scattering, the formula is merely the radar backscattering with the scattered direction changed to that of the bistatic direction and thus only the contributions from microwave roughness are included. However, in reflection from signals of opportunity including GNSS-R, the scattered direction is in the vicinity of the specular direction within a few degrees. The first difference of this close-to-specular scattering is that the bistatic scattering has both the mean field intensity and the covariance of the field. Thus, both the mean field intensity and covariance of field need to be included in the ${\sigma }_0$ [39]. The second difference is that even if the mean fields are ignored, the remaining bistatic covariance of field are strongly dependent on topographical elevations and slopes which can affect the ${\sigma }_0$ by many decibels. The topographical slopes are much less than that of microwave roughness. On the other hand, for L-band microwave remote sensing of soil moisture, the models of SMAP (soil moisture active and passive) and NISAR (NASA-ISRO SAR) primarily include the effects of microwave roughness with topography playing a lesser role.

In the past, there have been three popular models. One is the extension of the radar backscattering of microwave roughness to the bistatic direction as described earlier. The second model is the coherent model which uses the image theorem of a point source to obtain the reflected signal that captures the Fresnel zone effects. The results are multiplied by an attenuation factor $exp\left(-4k^2{h}_1^2{\cos }^2{\theta }_i\right)$, where $h_1$ is the root mean square (rms) height of microwave roughness, $k$ is the wave number, and ${\theta }_i$ is the incident angle. The attenuation factor is quite significant. For microwave roughness of $h_1=3$ cm, ${\theta }_i={40}^{\circ }$, the attenuation factor is −10 dB. The coherent model is calculating the mean field from small $h_1$ surfaces. The third model is the incoherent model based on the assumption that because of the topographical large elevation changes in land surfaces, the received signal is incoherent. The incoherent models are based on Geometric Optics (GO) model. The GO model is also an approximation of Kirchhoff integral using the method of stationary phase. In the GO model, the scattered intensity is proportional to the probability density function $p_{\mathrm{pdf}}\left(p,q\right)$ where $p$ and $q$ are the slopes in the horizontal $x$ and $y$ directions, respectively. The choice of pdf is Gaussian, so that $p_{\mathrm{pdf}}\left(p,q\right)=exp\left[-\left(p^2+q^2\right)/\left(2s^2\right)\right]/\left(2\pi s^2\right)$, where $s$ is the rms slope and a small number is usually used. It was first used in GNSS-R for ocean reflection [40] by truncating the ocean spectrum to eliminate the small roughness so as to derive the pdf of slopes. For GNSS-R land applications, it was also used by refs. [41,42]. However, a second version of geometric optics in ocean problem was proposed in ref. [43] in which the attenuation factor $exp\left(-4k^2{h}_1^2{\cos }^2{\theta }_i\right)$ is attached to account for the effects of microwave small roughness. It is labeled as the “Improved Geometric Optics Model (IGOM)”, which is GO-Att (Geometric Optics with Attenuation). The differences between the two GO models are significant even for $h_1=3$ cm. The two GO models have been applied to land surfaces [39,42].

We applied Kirchhoff models to signals of opportunity for both L-band and P-band. In our models, the surface profiles are composed of a summation of three kinds of roughness/topography (Figure 8):

$$z=f_1\left(x,y\right)+f_2\left(x,y\right)+f_3\left(x,y\right)$$

(69)

where $f_1\left(x,y\right)$ is the microwave roughness with rms height of 6 cm or less, and $f_2$ and $f_3$ are for the topography. In the CYGNSS investigations, extensive measurements have been taken to measure the rms heights and correlation lengths of the microwave roughness at San Luis Valley. The topography $f_3\left(x,y\right)$ is the coarse scale topography as given by DEM. It is labeled as “coarse” because the DEM is of horizontal resolution of 30 m. A linear interpolation is used to obtain $f_3\left(x,y\right)$, so that $f_3\left(x,y\right)$ corresponds to tilted planar patches with 30 m scale. The $f_2\left(x,y\right)$ is labeled as “fine-scale topography” that is in-between the coarse topography and the microwave roughness. The fine scale topography $f_2\left(x,y\right)$ will have rms height of 5 cm and above and horizontal correlations of 5 m to 10 m. The ratio of correlation length to rms height of fine-scale topography $f_2\left(x,y\right)$ is of the order of 100, which is 10 times that of the microwave roughness of $f_1\left(x,y\right)$. Recently, lidar measurements have been taken [44].

In our approach, we have used the Kirchhoff theory. Although the vector Kirchhoff theory is not valid for radar backscattering particularly for the VV polarization, it is quite accurate in the vicinity of the specular direction as shown by numerical solutions of Maxwell equations. In

Table 5. Summary of the two Kirchhoff approaches.

Models	Numerical Kirchhoff Approach (NKA) [13]	Analytical Kirchhoff Solution (AKS) [16]
Discretization	2 cm	30-m DEM patch
Monte Carlo simulations	Monte Carlo Speckle fluctuations	Analytical No Monte Carlo No fluctuations
CPU time for one DDM pixel of 15 km	Intensive 1 week for one DMM	Fast 1 h for one DDM f1 and f2 constant
Validation	Accurate benchmark based on brute force calculations	Validated by NKA
DEM Coarse f3	Planar with slope, deterministic	Planar with slope, deterministic
Fine scale f2: random	Monte Carlo average	Analytical average
Microwave f1: random	Monte Carlo average	Analytical average
Combining roughness	$f=f_1+f_2+f_3$	$f_{12}=f_1+f_2$ combined dividing line not needed
Spectrum W(k)	Can directly use W(k)	Can directly use W(k)
Histogram statistics of amplitude and phase	Yes	No

, we summarize our two approaches when applied to L-band signals of opportunity. The first approach, which we initially used, is a numerical Kirchhoff approach [13]. We used a patch size of 2 cm by 2 cm in the numerical discretization of calculating the Kirchhoff integral. The surface is characterized by the height function $f=f_1+f_2+f_3$. Because of fluctuations due to $f_1$ and $f_2$, Monte Carlo simulations were performed and averages were taken over realizations over $f_1$ and $f_2$. The results are treated as having benchmark accuracy as no approximations are used aside from the Kirchhoff integral. An intensive CPU with high-performance computation is required.

In the second approach, the Analytical Kirchhoff Solution (AKS) [16], we performed analytical averaging for the random surface. In this approach, we used $f_{12}=f_1+f_2$ combined. The roughness $f_{12}$ was super imposed on the $f_3$ planar patches of DEM. Then rough-surface scattering theory was applied to $f_{12}\left(x,y\right)$ with analytical solutions derived for the coherent waves and incoherent waves. The salient features of the AKS model are as follows. (i) Analytical expressions are obtained for both coherent and incoherent waves by taking analytical averaging over the random $f_{12}$; (ii) Monte Carlo simulations are not required, making the AKS computationally efficient; and (iii) the analytical solutions are expressed in terms of the spectrum so that the dividing line between microwave roughness and the fine scale topography is not required, and the rough surface spectrum derived from lidar elevation measurements can be incorporated directly. The results of AKS and Numerical Kirchhoff Approach (NKA) are indistinguishable for both the coherent waves and the incoherent waves. The agreements validate the AKS as NKA is a brute force accurate method using 2 cm discretization and high-performance computers. In Table 5, we compare the two approaches. Since the two approaches provide indistinguishable results, we shall label the results as the Kirchhoff approach. The numerical results of the Kirchhoff approach have been compared with that of GO [40] and GO-Att (IGOM) [42,43].

This section is organized as follows. In Section 3.2, the geometry of the SoOp problem is discussed. In Section 3.3, NKA is discussed. In Section 3.4, we discuss AKS. Numerical results generated with the Lidar data surface profile are shown in Section 3.5 and the CPU time comparison between AKS and NKA is presented in Section 3.6.

3.2. Coherent and Incoherent Models

The first two models initially proposed to the GNSS-R or SoOp application are the coherent model and the incoherent model, which only account for the coherent wave or the incoherent wave only, but not both. The geometry considered by the two models is shown in Figure 9.

The coherent model is:

$$\frac{P_r}{P_t}=\frac{G_t{G}_r{\lambda }^2}{{\left(4\pi \right)}^2{\left(R_{pt}+R_{pr}\right)}^2}{\left|R_{lr}\right|}^2{e}^{-4k^2{h}^2{\cos }^2\theta }$$

(70)

where $G_t$ and $G_r$ are the gain of transmitter and receiver, and $R_{pt}$ and $R_{pr}$ are the distance from transmitter to specular point and specular point to receiver. $\lambda$ is the wavelength, and $R_{lr}$ is the reflection coefficient of circular polarized wave given as:

$$R_{lr}=R_v+i{R}_h$$

(71)

k is the wave number, h is the rms height of microwave roughness, and $\theta$ is the incident angle at the specular point.

The coherent model is from a past model of rough-surface scattering and assumes a single elevation and the contributions arise from Fresnel zones. The attenuation caused by the surface roughness in the exponential term. Rigorous derivation for the coherent model can be found in the appendix of ref. [13].

The incoherent model is also based on a past scattering model and treats the GNSS-R as a special case of bistatic radar scattering with the scattering angle in the specular direction. In the bistatic radar equation, the incoherent model is:

$$\frac{P_r}{P_t}=\frac{G_t}{4\pi R_t{}^2}\frac{1}{4\pi R_r{}^2}\frac{G_r{\lambda }^2}{4\pi }\int \gamma dA$$

(72)

where $\gamma$ is the bistatic scattering coefficient of incoherent wave. Equation (72) is proportional to surface area. Based on geometric optics approximation, the bistatic scattering coefficient is given as:

$$\gamma =\frac{{\left|R_{lr}\right|}^2}{2s^2}$$

(73)

where s is the rms slope of the surface, and $R_{lr}$ is the reflection coefficient of circular polarized wave given in Equation (71). In the case of Gaussian correlation function, $s=\sqrt{2}h/l$, where $l$ is the correlation length of the rough surface.

Both the coherent and incoherent model consider only a single-scale roughness. The effects of topography elevations and slopes are not included in the coherent model or the incoherent model.

3.3. Geometric Descriptions of SoOp: Topography and Rough Surface

We consider the surface height $f\left(x,y\right)$ as composed of a summation of three kinds of roughness/topography (Figure 8).

$$f\left(x,y\right)=f_1\left(x,y\right)+f_2\left(x,y\right)+f_3\left(x,y\right)$$

(74)

We define

$$f_{12}\left(x,y\right)=f_1\left(x,y\right)+f_2\left(x,y\right)$$

(75)

where $f_1\left(x,y\right)$ is the microwave roughness with rms height of 0.25 wavelength or less. The $f_3\left(x,y\right)$ is the DEM topography at 30 m scale. The intermediate scale $f_2\left(x,y\right)$ is labeled as fine-scale topography between microwave roughness and DEM topography. Previously, profiles of $f_2\left(x,y\right)$ were largely unknown. However, recent lidar measurements confirm the existence of non-zero $f_2\left(x,y\right)$. A major difference between $f_1\left(x,y\right)$ and $f_2\left(x,y\right)$ is that the ratio of correlation length to rms height of $f_2\left(x,y\right)$ is 10 times larger than that of $f_1\left(x,y\right)$. We shall assume that $f_1\left(x,y\right)$ and $f_2\left(x,y\right)$ are stochastic Gaussian processes that are statistical homogeneous. The $f_3\left(x,y\right)$, derived from DEM, is deterministic and will be assumed to be consisting of planar patches of 30 m with slopes $p_3$ and $q_3$ in the $x$ and $y$ directions respectively. This means that the second-order derivatives of $f_3\left(x,y\right)$ are equal to zero. Thus, the geometry consists of DEM 30 m planar patches with slope and stochastic roughness $f_1\left(x,y\right)$ and $f_2\left(x,y\right)$ superimposed on the planar patch.

The descriptions of $f_1\left(x,y\right)$, $f_2\left(x,y\right)$, and $f_{12}\left(x,y\right)$ are shown in

Table 6. Descriptions of $f_1\left(x,y\right)$, $f_2\left(x,y\right)$, and $f_{12}\left(x,y\right)$.

Scales	$\mathbf{Microwave} \mathbf{Roughness} {\mathit{f}}_1\left(\mathit{x},\mathit{y}\right)$	$\mathbf{Fine}-\mathbf{Scale} \mathbf{Topography} {\mathit{f}}_2\left(\mathit{x},\mathit{y}\right)$	$\mathbf{Combined} \mathbf{Profile} {\mathit{f}}_{12}\left(\mathit{x},\mathit{y}\right)$
Correlation Function	$h_1^2{C}_1\left(x,y\right)$	$h_2^2{C}_2\left(x,y\right)$	$h^{2}C\left(x,y\right)$
Spectrum	$W_1\left(k_x,k_y\right)$	$W_2\left(k_x,k_y\right)$	$W\left(k_x,k_y\right)$

. For real-life profiles, it can be difficult to form a dividing line between $f_1$ and $f_2$. Thus, the last column $f_{12}\left(x,y\right)$ is a combination of microwave roughness and fine-scale topography without a dividing line. In Table 6, we also list the correlation functions and the spectral densities. The combined case of $f_{12}\left(x,y\right)$ is the general case with correlation function $h^{2}C\left(x,y\right)$ and spectral density of $W\left(k_x,k_y\right)$. The decompositions of $f_{12}\left(x,y\right)$ separately into $f_1\left(x,y\right)$ and $f_2\left(x,y\right)$ are special cases of $f_{12}\left(x,y\right)$.

Since $f_{12}\left(x,y\right)$ is statistically homogeneous, the correlation function is:

$$h^{2}C\left(x,y\right)=\langle f_{12}\left(x^{\prime },y^{\prime }\right)f_{12}\left(x^{\prime }+x,y^{\prime }+y\right)\rangle ={\iint }_{-\infty }^{+\infty }d{k}_{x}d{k}_{y}W\left(k_x,k_y\right)e^{i{k}_{x}x+i{k}_{y}y}$$

(76)

with $C\left(0,0\right)=1$. Then we have

$$h^2={\iint }_{-\infty }^{+\infty }d{k}_{x}d{k}_{y}W\left(k_x,k_y\right)$$

(77)

For the case of isotropic rough surface:

$$h^{2}C\left(x,y\right)=h^{2}C\left(\rho \right)=2\pi {\int }_{-\infty }^{+\infty }d{k}_{\rho }{k}_{\rho }W\left(k_{\rho }\right)J_0\left(k_{\rho }\rho \right)$$

(78)

where $J_0$ is the Bessel function of zero-th order. Let the GNSS-R transmitter and receiver be located in the $xz$ plane, as shown in Figure 10. Note that the transmitter is located at:

$$T_x=\left(x_t,0,h_t\right)$$

(79)

Let the receiver be located at:

$$R_x=\left(x_r,0,h_r\right)$$

(80)

Since $x_t$ is negative, let

$$x_s=-x_t$$

(81)

Then

$$x_r+x_s=d; x_s=-x_t=\frac{d{h}_t}{h_r+h_t}; x_r=\frac{d{h}_r}{h_r+h_t}$$

(82)

In GNSS-R applications, the receivers are in the specular directions, which means d can be expressed in terms of the incident angle ${\theta }_i$:

$$\frac{d}{h_r+h_t}=\tan {\theta }_i=\tan {\theta }_s$$

(83)

Thus, we have

$$x_s=-x_t=h_t\tan {\theta }_i;x_r=h_r\tan {\theta }_s$$

(84)

3.4. Numerical Kirchhoff Approach (NKA)

Let a point on the observed surface have position $\overline{r}=\left(x,y,z\right)$ where $z=f_1\left(x,y\right)+f_2\left(x,y\right)+f_3\left(x,y\right)$. The stochastic roughness and fine-scale topography are superimposed on the planar patch of $f_3\left(x,y\right)$:

$$f_3\left(x,y\right)=z_n+p_{3n}\left(x-x_n\right)+q_{3n}\left(y-y_n\right)$$

(85)

where $p_{3n}$ and $q_{3n}$ are the coarse topography slopes in x and y directions. The scattered field from the land surface is given by the following equation under Kirchhoff approximation [39]:

$${\overline{E}}_s=\frac{ik}{4\pi }\sqrt{\frac{P_t{G}_t{\eta }_0}{2\pi }}{\iint }_{S}d\overline{r}\frac{e^{R_t+R_r}}{R_t{R}_r}\left(\stackrel{=}{I}-{\widehat{k}}_s{\widehat{k}}_s\right)·\overline{F}\left(\alpha ,\beta \right)$$

(86)

where the $\overline{F}\left(\alpha ,\beta \right)$ accounts for the polarization of the scattered wave, which is given as:

$$\begin{gathered} \begin{array}{c}\overline{F}\left(\alpha ,\beta \right)=\sqrt{1+{\alpha }^2+{\beta }^2}[\left(-1+R_h\right)({\widehat{e}}_i·{\widehat{k}}_i){\widehat{q}}_i+\left(1+R_v\right)\left({\widehat{e}}_i·{\widehat{p}}_i\right)\widehat{n}\times {\widehat{q}}_i+{\widehat{k}}_s\times \left[\left(1+R_h\right)\left({\widehat{e}}_i·{\widehat{q}}_i\right)\right]\widehat{n}\times {\widehat{q}}_i\\ \\ +\left(1-R_v\right)\left({\widehat{e}}_i·{\widehat{p}}_i\right)(\widehat{n}·{\widehat{k}}_i){\widehat{q}}_i]\end{array} \end{gathered}$$

(87)

In the expression above, $\alpha$ and $\beta$ are the surface slopes in x and y directions for $f_1\left(x,y\right)+f_2\left(x,y\right)+f_3\left(x,y\right)$, which is given as:

$$\alpha =\frac{d}{dx}f\left(x,y\right)=\frac{d}{dx}{f}_1\left(x,y\right)+\frac{d}{dx}{f}_2\left(x,y\right)+p_{n_3}$$

(88)

$$\beta =\frac{d}{dy}f\left(x,y\right)=\frac{d}{dy}{f}_1\left(x,y\right)+\frac{d}{dy}{f}_2\left(x,y\right)+q_{n_3}$$

(89)

${\widehat{e}}_i$ is the incident wave polarization vector; in the GNSS-R application, the incident wave is circular polarized with

$${\widehat{e}}_i={\widehat{v}}_i+i{\widehat{h}}_i$$

(90)

The incident wave direction is given by ${\widehat{k}}_i$ and the scattering direction is given by ${\widehat{k}}_s$, which are given by the following equations:

$$\begin{array}{c}{\hat{k}}_i=\frac{{\overline{r}}^{\prime }-{\overline{R}}_T}{R_t}=\frac{(x^{\prime }-x_t)\hat{x}+(y^{\prime }-y_t)\hat{y}+(z^{\prime }-z_t)\hat{z}}{\sqrt{{(x^{\prime }-x_t)}^2+{(y^{\prime }-y_t)}^2+{(z^{\prime }-z_t)}^2}}\\ {\hat{k}}_s=\frac{{\overline{R}}_r-{\overline{r}}^{\prime }}{R_r}=\frac{(x_r-x^{\prime })\hat{x}+(y_r-y^{\prime })\hat{y}+(z_r-z^{\prime })\hat{z}}{\sqrt{{(x_r-x^{\prime })}^2+{(y_r-y^{\prime })}^2+{(z_r-z^{\prime })}^2}}\end{array}$$

(91)

The distance between $\overline{r}$ and receiver is given by $R_r$. Due to the high elevation in the orbit of the transmitter satellite, the incident direction is the same for the observation area. The receiver orbit is much lower and the scattering direction changes according to the position of the location point (x, y, z) on the height profile. Let $\widehat{n}$ be the unit normal at the point (x, y, z):

$$\widehat{n}=\frac{-\frac{df\left(x,y\right)}{dx}\widehat{x}-\frac{df\left(x,y\right)}{dy}\widehat{y}+\widehat{z}}{\sqrt{1+{\left(\frac{df\left(x,y\right)}{dx}\right)}^2+{\left(\frac{df\left(x,y\right)}{dy}\right)}^2}}$$

(92)

${\widehat{p}}_i$ and ${\widehat{q}}_i$ are local polarization vectors defined at the point:

$$\begin{array}{c}{\widehat{q}}_i=\frac{{\widehat{k}}_i\times \widehat{n}}{\left|{\widehat{k}}_i\times \widehat{n}\right|}\\ {\widehat{p}}_i={\widehat{q}}_i\times {\widehat{k}}_i\end{array}$$

(93)

$R_v$ and $R_h$ are the Fresnel reflection coefficients with local incidence angle:

$${\theta }_{il}=-a\cos (\widehat{n}·{\widehat{k}}_i)$$

(94)

$$\begin{array}{c}{R}_h=\frac{k\cos {\theta }_{il}-\sqrt{k_1^2-k^2{\sin }^2{\theta }_{il}}}{k\cos {\theta }_{il}+\sqrt{k_1^2-k^2{\sin }^2{\theta }_{il}}}\\ R_v=\frac{{\epsilon }_{1}k\cos {\theta }_{il}-{\epsilon }_0\sqrt{k_1^2-k^2{\sin }^2{\theta }_{il}}}{{\epsilon }_{1}k\cos {\theta }_{il}+\nu \sqrt{k_1^2-k^2{\sin }^2{\theta }_{il}}}\end{array}$$

(95)

The integral in (86) is performed with brute force using a high-performance computer with 2 cm discretization at the L-band. To obtain the coherent wave and incoherent intensity, we use Monte Carlo averaging, and thus the mean field intensity and the intensity of covariance of fields are given as:

$${|\langle {\overline{E}}_s^{mean field}\rangle |}^2=|\frac{1}{N}\sum _1^N{{\overline{E}}_{sn}|}^2$$

(96)

$${|\langle {\overline{E}}_s^{Cov}\rangle |}^2=\frac{1}{N}\sum _1^N{|{\overline{E}}_{sn}-\langle {\overline{E}}_s\rangle |}^2$$

(97)

After obtaining the intensity for mean field and covariance of field, we next obtain the bistatic scattering coefficients by:

$${\gamma }_{mean field/Cov}=\frac{{\left(4\pi \right)}^2{R}_t^2{R}_r^2}{2\eta P_t{G}_{t}A\cos {\theta }_i}{I}_s^{mean field/Cov}$$

(98)

3.5. Analytical Kirchhoff Solutions (AKS)

In the analytical Kirchhoff solutions, we divided the land surface into patches and apply far-field approximations to calculate the scattering from each patch. We first considered the mean field and Covariance of fields from a single patch of size L. Let the center of the patch be given by

$${\overline{r}}_n=x_n\widehat{x}+y_n\widehat{y}+z_n\widehat{z}$$

(99)

On this patch, $z_n=f_3\left(x_n,y_n\right)$, and the slopes of this patch are

$$\begin{array}{c}\frac{d{f}_3\left(x_n,y_n\right)}{dx}=p_{3n}\\ \frac{d{f}_3\left(x_n,y_n\right)}{dy}=q_{3n}\end{array}$$

(100)

The distance between the patch center and the transmitter and the receiver are respectively $R_{nt}$ and $R_{nr}$:

$$\begin{array}{c}{R}_{nt}=\sqrt{{\left(x_n-x_t\right)}^2+y_n^2+{\left(z_n-h_t\right)}^2}\\ R_{nr}=\sqrt{{\left(x_n-x_r\right)}^2+y_n^2+{\left(z_n-h_r\right)}^2}\end{array}$$

(101)

In order to apply the far field approximation to the patch, we need to have:

$$\begin{array}{c}{R}_{nt},R_{nr}\gg \frac{L^2}{\lambda }\\ L\ll \sqrt{\lambda R_{nt}},\sqrt{\lambda R_{nr}}\end{array}$$

(102)

In GNSS-R, the elevation of the transmitter is much larger than that of the receiver, which is $R_r=500 \mathrm{km}$. For the P band, we have $freq=370$ MHz, $\lambda =0.81 \mathrm{m}$. For the L band, $freq=1.575$ GHz, $\lambda =0.19 \mathrm{m}$. Thus,

$$\begin{array}{c}L\ll 308 \mathrm{m}, L-band\\ L\ll 636 \mathrm{m}, P-band\end{array}$$

(103)

is required over the domain of integration. If the discretization is less than that requirement, then the variations of the scattered field direction and the variations of the phase factor need to be variables over the domain the integration. With respect to the center of the patch, the wave vector for the incident wave is given by:

$$\begin{array}{c}\begin{array}{c}{\overline{k}}_{in}=k_{inx}\widehat{x}+k_{iny}\widehat{y}-k_{inz}\widehat{z}\\ k_{inx}=\frac{k\left(x_n-x_t\right)}{R_{nt}}\end{array}\\ \begin{array}{c}{k}_{iny}=\frac{k{y}_n}{R_{nt}} \\ k_{inz}=\frac{k\left(h_t-z_n\right)}{R_{nt}}\end{array}\end{array}$$

(104)

The incident polarization vectors on the patch are:

$$\begin{array}{c}{\widehat{h}}_{in}=\frac{\widehat{z}\times {\overline{k}}_{in}}{\left|\widehat{z}\times {\overline{k}}_{in}\right|}\\ {\widehat{v}}_{in}={\widehat{h}}_{in}\times {\overline{k}}_{in}\end{array}$$

(105)

For the scattering wave, we have the wave vector defined as follows:

$$\begin{array}{c}\begin{array}{c}{\overline{k}}_{sn}=k_{snx}\widehat{x}+k_{sny}\widehat{y}+k_{snz}\widehat{z}\\ k_{snx}=\frac{k\left(x_r-x_n\right)}{R_{nr}}\end{array}\\ \begin{array}{c}{k}_{sny}=-\frac{k{y}_n}{R_{nr}}\\ k_{snz}=\frac{k\left(h_r-z_n\right)}{R_{nr}}\end{array}\end{array}$$

(106)

The corresponding polarization vectors are:

$$\begin{array}{c}{\widehat{h}}_{sn}=\frac{\widehat{z}\times {\overline{k}}_{sn}}{\left|\widehat{z}\times {\overline{k}}_{sn}\right|}\\ {\widehat{v}}_{sn}={\widehat{h}}_{sn}\times {\overline{k}}_{sn}\end{array}$$

(107)

The local incident angle on the patch is given by:

$$\cos {\theta }_{in}=-\widehat{n}·{\overline{k}}_{in}$$

(108)

The mean field intensity and covariance of fields correspond to the average and the variance of the integral $I_n^{\left(N\right)}$. They are evaluated using methods described in page 80 of Tsang and Kong volume 3 [39]. The average of $I_n^{\left(N\right)}$, denoted by $I_n^{\left(N\right)}$, is expressed as:

$$\langle I_n^{\left(N\right)}\rangle =kL{e}^{-\frac{k_{dnz}^2{h}^2}{2}}sinc\left[\left(\frac{k_{dnx}}{k_{dnz}}+p_{3n}\right)\frac{k_{dnz}L}{2}\right]sinc\left[\left(\frac{k_{dny}}{k_{dnz}}+q_{3n}\right)\frac{k_{dnz}L}{2}\right]$$

(109)

where $k$ is the wavenumber, ${\overline{k}}_{dn}={\overline{k}}_{in}-{\overline{k}}_{sn}=k_{dnx}\widehat{x}+k_{dny}\widehat{y}+k_{dnz}\widehat{z}$, and $L=30$ m is the planar patch size of DEM $f_3\left(x,y\right)$. The sinc functions represent the scattering pattern of coherent of the patch. $k_{dnx}$ and $k_{dny}$ represent the separation of the patch from the specular point, and $p_{3n}$ and $q_{3n}$ are the slopes of the patch for x and y directions.

The variance of $I_n^{\left(N\right)}$ is, assuming isotropic roughness and isotropic fine scale topography of $f_{12}\left(x,y\right)$:

$$D_{I_n}^{\left(N\right)}=\langle {|I_n^{\left(N\right)}|}^2\rangle -|\langle {I_n^{\left(N\right)}}^2\rangle |=2\pi k^2{\int }_0^{\infty }d\rho \rho J_0\left(\rho \sqrt{{\left(k_{dnx}+p_{3n}{k}_{dnz}\right)}^2+{\left(k_{dny}+q_{3n}{k}_{dnz}\right)}^2}\right)\left[e^{-k_{dnz}^2{h}^2\left(1-C\left(\rho \right)\right)}-e^{-k_{dnz}^2{h}^2}\right]$$

(110)

3.5.1. Multiple DEM Patches

For land surfaces which are divided into multiple DEM patches with size of $L$ (e.g., $L=30$ m), the mean fields were obtained by field summation from each patch. The covariance of fields was obtained by incoherent summation from each patch.

3.5.2. Mean Field

The bistatic scattering coefficient for mean field intensity of an area with N patches is given as the following equation.

$${\gamma }_{mean field}=\frac{R_t^2{R}_r^2\cos {\theta }_i}{N\pi }{\left|{\displaystyle \sum }_{n=1}^N\frac{e^{ik\left(R_{nt}+R_{nr}\right)}}{R_{nt}{R}_{nr}}\frac{{\widehat{v}}_{sn}{R}_v\left({\theta }_{in}\right)+i{\widehat{h}}_{sn}{R}_h\left({\theta }_{in}\right)}{\sqrt{2}}\langle I_n^{\left(N\right)}\rangle \right|}^2$$

(111)

Small slope approximation is applied to $f_3$ to simply the expression. Reference [16] provides more details. To calculate the ${\gamma }_{mean field}$ within a certain area, we need first to obtain the $\langle I_n^{\left(N\right)}\rangle$ from each patch using Equation (109). The exponential term $\frac{e^{ik\left(R_{nt}+R_{nr}\right)}}{R_{nt}{R}_{nr}}$ keeps track of the phase of mean field for each patch. $R_{nt}$ and $R_{nr}$ are the distance between the specular point and transmitter and receiver respectively. The term $\frac{{\widehat{v}}_{sn}{R}_v\left({\theta }_{in}\right)+i{\widehat{h}}_{sn}{R}_h\left({\theta }_{in}\right)}{\sqrt{2}}$ represents the reflection of circular polarized wave on the patch, and the Fresnel reflection coefficients with local incident angle ${\theta }_{in}$ are given by:

$$\begin{array}{c}{R}_{hn}=\frac{k\cos {\theta }_{il}-\sqrt{k_1^2-k^2{\sin }^2{\theta }_{in}}}{k\cos {\theta }_{il}+\sqrt{k_1^2-k^2{\sin }^2{\theta }_{in}}}\\ R_{vn}=\frac{{\epsilon }_{1}k\cos {\theta }_{in}-{\epsilon }_0\sqrt{k_1^2-k^2{\sin }^2{\theta }_{in}}}{{\epsilon }_{1}k\cos {\theta }_{in}+\nu \sqrt{k_1^2-k^2{\sin }^2{\theta }_{in}}}\end{array}$$

(112)

3.5.3. Covariance of Fields

$${\gamma }_{Co var}=\frac{R_t^2{R}_r^2\cos {\theta }_i}{N\pi }{\displaystyle \sum }_{n=1}^N\frac{1}{R_{nt}^2{R}_{nr}^2}\frac{{\left|R_v\left({\theta }_{in}\right)-R_h\left({\theta }_{in}\right)\right|}^2}{4}{D}_{I_n}^{\left(N\right)}$$

(113)

To calculate the ${\gamma }_{Co var}$ within a certain area, the $D_I^{\left(N\right)}$ for each patch is calculated by Equation (113). The reflection coefficient for the circular polarized wave is calculated by $\frac{{\left|R_v\left({\theta }_{in}\right)-R_h\left({\theta }_{in}\right)\right|}^2}{4}$. Equation (45) is obtained by small slope approximation to $f_3$. The detailed derivation can be found in ref. [16].

3.6. Two Geometric Optics Approaches

The geometric optics consists of making a high-frequency approximation of $k\to \infty$. Since microwave roughness does not fall in that category, $f_1\left(x,y\right)$ is ignored in the traditional geometric optics approach. The stationary phase approximation is then applied to the Kirchhoff integral to the phase term with $f_{23}\left(x,y\right)=f_2\left(x,y\right)+f_3\left(x,y\right)$. The results are independent of frequency. The GO was first applied to GNSS-R for ocean [40]. However, to account for the effects of small roughness of $f_1\left(x,y\right)$, an improved geometric optics model (IGOM) was proposed in ref. [43], which was recently adopted in CYGNSS [42].

In IGOM, an attenuation term $e^{-4k^2{h}_1^2{\cos }^2{\theta }_i}$ is added. The approach requires the decomposition of $f_{12}$ separately into $f_1+f_2$.

The bistatic scattering coefficient for GO is:

$${\gamma }^{GO}=\frac{R_t^2{R}_r^2}{N\cos {\theta }_i}{\displaystyle \sum }_{n=1}^N\frac{1}{R_{nt}^2{R}_{nr}^2}\frac{{\left|R_v\left({\theta }_{in}\right)-R_h\left({\theta }_{in}\right)\right|}^2}{4}\frac{1}{2s_{2n}^2}\frac{{\left|k_{dn}\right|}^4}{{\left|k_{dnz}\right|}^4}{e}^{-\frac{{\left(\frac{k_{dnx}}{k_{dnz}}+p_{3n}\right)}^2+{\left(\frac{k_{dny}}{k_{dnz}}+q_{3n}\right)}^2}{2s_{2n}^2}}$$

(114)

where N is the number of patches in the observation area, $\frac{1}{2}\left(R_v+R_h\right)$ is the reflection coefficient of circular polarized wave. $s_2=\sqrt{2}\frac{h_2}{l_2}$ is the rms slope of $f_2$ assuming a Gaussian correlated surface.

The IGOM is given by

$${\gamma }^{IGOM}={\gamma }^{GO}{e}^{-4k^2{h}_{1n}^2{\cos }^2{\theta }_{in}}$$

(115)

Notice that both GO and IGOM are derived from the Kirchhoff integral based on stationary phase approximation. The differences between GO and IGOM are small for small $h_1$. In

Table 7. Typical values of the attenuation factor in GO-Att. $k=33 m^{-1},{\theta }_i={40}^{°}$.

${\mathit{h}}_1\left(\mathbf{c}\mathbf{m}\right)$	$-4{\mathit{k}}^2{\mathit{h}}_1^2{\mathbf{cos}}^2{\mathit{\theta }}_{\mathit{i}}$	$\mathit{e}\mathit{x}\mathit{p}\left(-4{\mathit{k}}^2{\mathit{h}}_1^2{\mathbf{cos}}^2{\mathit{\theta }}_{\mathit{i}}\right)\left(\mathbf{d}\mathbf{B}\right)$
$1$	$0.26$	$-1.11$
$2$	$1.02$	$-4.44$
$3$	$2.30$	$-10.0$
$6$	$6.39$	$-40.0$

, we tabulate $e^{-4k^2{h}_{1n}^2{\cos }^2{\theta }_i}$ at the L-band. As shown in Table 7, with $h_1=3 \mathrm{cm}$, the difference between GO and IGOM is 10 dB. In mountainous regions, when $h_1=6 \mathrm{cm}$, the difference between GO and IGOM is 40 dB.

3.7. Numerical Results for L-Band and P-Band

In this section, we present the results for L- and P-band simulations. In Figure 11, we show the bistatic scattering coefficient near the specular direction of AKS, GO, and IGOM in both L- and P-bands. It is shown that the GO stayed the same for L- and P-band and IGOM had a 2 dB difference while AKS showed a difference about 8 dB in the specular direction. Additionally, notice that at the L-band, the AKS solution overlaps with the GO. This is due to the large $k{h}_2$ value in the L-band since the rms height is already 1 quarter of the wave length in the L-band, in which case the Geometric optic limit is satisfied. However, in the P-band, the rms height is only 1/16 of the wavelength.

The mean field bistatic coefficient and covariance are shown in Figure 12 for the P-band and L-band. The observation area was centered at $\left(0,0,0\right)$, which is the specular point. The incident angle and scattering angle with respect to the specular point was fixed at ${\theta }_i=40°$. The transmitter and receiver positions are given by Equations (79) and (80). Equations (111) and (113) were used to calculate the bistatic scattering coefficients for mean field intensity and Covariance of fields. The rms height and correlation length for $f_2\left(x,y\right)$ were 5.78 cm and 3.75 m respectively for the entire 5 km by 5 km area. The land topography data were obtained from SRTM [45] centered at at 37.2116° N, 105.9756° W. The land surface was interpolated into 60 m × 60 m grids, and the slopes for coarse topography was obtained by taking the different between adjacent data points.

The bistatic scattering coefficient for the mean field at P-band was much larger than the covariance of fields for an area as large as 5 km. For the L-band, the mean field was much smaller than the covariance of fields. This is due to the large rms height of 5.78 cm for $f_2\left(x,y\right)$. Monte Carlo was performed to obtain the mean and covariance of fields for NKA. For this case, the number of realizations $N=300$ was used. Notice that for the L band, the NKA does not agree with AKS for mean field in larger areas. This is due to the fact that when the area increases, there will be more additions of complex numbers. To obtain convergent results, more realizations are required for NKA.

3.8. L-Band: Track-Wise Comparison of DDM with CYGNSS Data

We compared simulated Delay Doppler Maps (DDM) with CYGNSS data. The track-wise comparison was performed with CYGNSS v3.0 data from Physical Oceanography Distributed Active Archive Center (PODAAC). We used the track data collected near the Z1 cal/val site (${37}^{°}{11}^{\prime }{26.26}^{″} \mathrm{N}$,${105}^{°}{59}^{\prime }{31.64}^{″} \mathrm{W}$) in San Luis Valley, CO, USA. The data were collected on day 301 of 2019 by craft 02 in channel 3. The time index is from 101,225 to 101,246. The specular points of the track of data we picked start from (${37.1817}^{°} \mathrm{N}$,${106.1732}^{°} \mathrm{W}$), and end at (${37.2451}^{°} \mathrm{N}$,${105.4590}^{°} \mathrm{W}$).

In Figure 13, a track-wise comparison for the peak value bistatic radar cross section (BRCS) between CYGNSS data and AKS is shown. The soil permittivity was selected as $\epsilon =3.293+0.198i$ based on the ground soil moisture measurement. We kept $h_1=1 \mathrm{cm}$, $h_2=5 \mathrm{cm}$, $l_1=10 \mathrm{cm}$, and changed $l_2$ from $50h_1$ to $130h_1$ for mountainous and planar areas respectively. The table is shown in

Table 8. Correlation length $l_2$ and coarse topography slope.

Longitude	${\mathit{l}}_2$$(\mathbf{Times} \mathbf{of} {\mathit{h}}_2)$	$\langle \sqrt{{\mathit{p}}_3^2+{\mathit{q}}_3^2}\rangle$
<−105.90	125	/
−105.86	69	0.0599
−105.83	53	0.0785
−105.79	52	0.0797
−105.76	50	0.0837
−105.73	61	0.0679
−105.69	83	0.0503
−105.66	74	0.0569
−105.62	93	0.0448
−105.59	126	0.0331
>−105.56	125	/

.

The size of the area that contributes to the Delay–Doppler along the track is 45,000 km². The $p_{3n}$ and $q_{3n}$ of the 30 m patches were derived from DEM over this area. As is shown in Figure 13, the simulation results of AKS are in good agreement with the data. We also show the mean slope using the scale on the right vertical axis. The result of Bistatic coefficient is opposite to that of the topographical slope. When the topographical slope is small, the bistatic coefficient is large. When the slope is large, the bistatic coefficient is small. The result clearly shows the importance of topography in reflectometry. One limitation of this comparison is that $h_1=1 \mathrm{cm}$ is assumed to be small. The small $h_1$ means that the difference between GO and IGOM is small.

With the incorporation of SAR technology in SoOP, we are studying the specular reflections with much finer spatial resolutions. We are also studying the vegetation and forest effects using the full-wave simulations methods.

4. Rough Surface

In this section, we review our recent method and results of calculating the scattering of the snow/soil interface at L-, C-, X-, and Ku-bands [46,47]. Two distinct assessments of rough-surface scattering simulations are (i) what is the largest $kh$ simulated? (ii) Are the simulations using exponential correlation functions or Gaussian correlation functions? Physical models of rough-surface scattering have been studied with the two classical methods of the small perturbation method (SPM) and the Kirchhoff approach [39,48]. Advanced analytical methods include the advanced integral equation model [49] and small slope approximation and its extensions [50,51]. A description of roughness is through $kh$, which is the product of the EM wavenumber of the medium above the rough surface and the surface rms height $h$. For soil-surface scattering, the previous results of analytical models and numerical simulations have been limited to $kh<3$ because of the past focus on L-band. Using a times series of backscatter measurements of SMAP from HH and VV polarizations, both soil moisture and surface rms height were retrieved at 3 km resolution for the 13 April–7 July 2015 period of SMAP radar operations [52]. The product shows that the global median rough surface heights are 2 cm and rms heights are 5 cm in the mountain regions. At Ku band, for the case of rms height of 5 cm, $kh$ is 18 for air/soil interface and 21.6 for snow/soil interface. The larger $kh$ for the snow/soil interface is because of the larger wavenumber in snow than in air. Thus, the past limit of $kh<3$ has limited applications to the L-band and C-band. Recently, we performed full-wave simulations with $kh$ up to 15, thus widening the applicability of full-wave simulations up to the Ku-band.

The soil surface can be described by a stationary Gaussian random process, so that knowledge of its covariance function is sufficient to describe its properties. The covariance function is further parametrized in terms of its rms height and correlation length. Ground measurements have been made of these properties [53,54,55]. The ground measurements of correlation lengths have a maximum of 10 cm. We label these roughness measurements as “limited correlation length up to 10 cm”. However, in global retrieval of soil moisture using six months of the NASA Soil Moisture Active Passive (SMAP) radar data at the L-band, the roughness was modeled as having a constant ratio of correlation length to rms height [52,56]. The ratios range from 5 to 20. This means that if the rms height is 5 cm, then the correlation length 100 cm if the ratio is 20. The two assumptions of “limited correlation length” and “constant ratio” are widely different for rms heights beyond 2 cm. The rms heights in mountainous regions are large and can be up to 5 cm. Based on our recent simulations, we found that the constant ratio approach provides more acceptable results. Surface scattering also depends on the soil permittivity, which in turn depends on soil moisture (which describes the volume of water present per unit volume of soil) and texture (which describes soil composition). Given these parameters, empirical models are available to calculate the soil permittivity [57,58]. In addition to soil properties, land cover including litter and vegetation as well as rock outcrops impact the spatial variability of surface permittivity.

We applied the Method of Moment (MoM) in full-wave simulations of numerical solutions of Maxwell’s equations (NMM3D) to L-band radar backscatter for the SMAP mission [25,26]. The full-wave simulations were used to generate look-up table (LUT) [27]. The LUT were initially used for air/soil interfaces. However, the LUT are based on incident angles of the upper medium and the relative dielectric constant between the two media on the two sides of the rough surfaces. By adjusting the relative dielectric constants and the incidence angle using Snell’s law, the NMM3D LUT are also applicable for all combinations of relative dielectric constants such as rough interfaces of snow/soil, air/snow, and snow/permafrost, etc. Both the AIEM model and the NMM3D have been previously limited to $kh<3$.

Recently, NMM3D calculations were performed with $kh$ up to 15 ($h=4.16 \mathrm{cm}$ for 17.2 GHz) [46,47].

Table 9. kh of 5-cm rms height for air/soil and snow/soil interface from L- to Ku-band.

$\mathbf{rms} \mathbf{Height}, \mathit{h} = 5 \mathbf{c}\mathbf{m}$	L Band (1.26 GHz)	S Band (2.5 GHz)	C Band (5.4 GHz)	X Band (9.6 GHz)	Low Ku Band (13.6 GHz)	High Ku Band (17.2 GHz)
kh, air/soil interface	1.32	2.62	5.66	10.06	14.25	18.02
kh, snow/soil interface	1.58	3.14	6.79	12.07	17.10	21.63

lists how $kh$ increases for $h$ = 5 cm from the L-band to the Ku-band for the air/soil interface and snow/soil interface. With $kh=15$, we are able to reach the Ku-band for the calculation for rms height up to 5 cm for most soil surfaces in the world. In this section, we describe how we calculate the rough soil-surface scattering using the volume integral equation (VIE) method with periodic boundary conditions.

Three-dimensional full-wave simulations of rough-surface scattering have drawn significant interest in recent years because of the advance of Computational electromagnetics (CEM), High-Performance Computing (HPC), and commercial software. We have applied MoM accelerated by SMCG Other methods and commercial software have also been used [26,27]. Among user-friendly commercial software are HFSS, which are based on the Finite element method (FEM), FEKO, which is based on the Method of Moments, accelerated by Multilevel Fast Multipole (MLFMA), and CST, which is based on the Finite Difference Time Domain method (FDTD). The advancement of these methods in rough-surface scattering can be assessed by the largest kh value calculated and the ability to deal with rough surfaces with exponential correlation functions. Soil surfaces are well described by an exponential correlation function as concluded in the paper by Oh et al. 1992 [59]. Surfaces with exponential correlation functions have the fine-scale features that make them more challenging than surfaces with Gaussian correlation functions. For HFSS, it was reported in paper by H. Lawrence et al. [60]. The paper showed results of exponential correlation functions with kh = 1. For CST, it was reported that rough surface simulations were performed for Gaussian correlation functions with h = 0.1 wavelength, which means kh = 0.63 [61]. For FEKO, it was reported in the paper by Y. W. Wei et al. that rough-surface scattering was calculated with Gaussian correlation functions for h = 0.06 wavelength, which means kh = 0.418 [62]. The extended boundary condition technique has been used for exponential correlation function with h at 0.044 wavelength corresponding to kh = 0.27 [63]. In our previous work, we reported the SMCG method in MoM with RWG basis functions with kh up to 4 and exponential correlation functions [27,64]. At this moment, the kh values of the reported full-wave simulations are much smaller than the kh = 15 reported in this section. To summarize, we previously performed NMM3D simulations using MoM–SMCG for rough soil surfaces with exponential correlation functions up to kh = 4. Other methods including FEM, FDTD and Extended Boundary condition (EBC), HFSS, FEKO, and CST, etc. have been limited to even smaller kh: less than 1.

4.1. Formulation of VIE Using Periodic Boundary Conditions and Periodic Half Space Dyadic Green’s Function

We applied the VIE approach with periodic boundary conditions to calculate scattering by a rough soil surface. Since we only consider rough surfaces, the formulation of combined volume and surface scattering in Tan 2016 is simplified [65].

Consider an incident plane wave incident onto a rough surface. The media above and below the rough surface are labeled medium 0 and medium 1 respectively. The permittivities are $\epsilon$ and ${\epsilon }_1$ respectively. The wavenumbers are $k$ and $k_1$ respectively. The medium, 0, can be air or snow with $\epsilon$ =${\epsilon }_0$ for air and $\epsilon$ =${\epsilon }_{snow}$ for snow. The ${\epsilon }_1$ in this paper is ${\epsilon }_1$ =${\epsilon }_{soil}$ for soil. For a single realization of a rough surface (Figure 14) with maximum height $f_{max}$ and minimum height $f_{min}$, we set $\mathrm{z}=0$ at $f_{\min }$ and $z$ =$d$ at $f_{\max }$ so that

$$d=f_{max}-f_{min}$$

(116)

Then, the region $0\le z\le d$ is an inhomogeneous region. The inhomogeneous permittivity in that region is ${\epsilon }_m\left(\overline{r}\right)$ is

$${\epsilon }_m\left(\overline{r}\right)=\{\begin{array}{c}\epsilon \mathrm{for} \overline{r} \mathrm{in} \mathrm{medium} 0\\ {\epsilon }_1 \mathrm{for} \overline{r} \mathrm{in} \mathrm{medium} 1\end{array}$$

(117)

The region is inhomogeneous, and the problem can be cast into VIE with the surface scattering problem converted into a volume scattering problem in the inhomogeneous region. We discretize the inhomogeneous regions into small cubes that are much smaller than a wavelength. Each small cube has permittivity either $\epsilon$ or ${\epsilon }_1$. Consider an incident plane wave launched from medium 0.

$${\overline{E}}_{inc}\left(\overline{r}\right)={\overline{E}}_0{e}^{i{k}_{ix}x+i{k}_{iy}y+i{k}_{iz}z}$$

(118)

where ${\overline{E}}_0$ is vertical polarized (TM) or horizontal polarized (TE) and

$$\begin{array}{c}{k}_{ix}=k\sin {\theta }_i\cos {\varphi }_i\\ k_{iy}=k\sin {\theta }_i\sin {\varphi }_i\\ k_{iz}=k\cos {\theta }_i\end{array}$$

(119)

Let ${\overline{k}}_{i\rho }=k_{ix}\widehat{x}+k_{iy}\widehat{y}$ be the horizontal component of the incident wave vector. Because we use the half space dyadic Green’s function as the background with the rough surface creating the inhomogeneous medium, there is a specular reflected wave when the rough surface is planar, which comes from reflection of a planar surface between medium 0 and medium 1 due to the planar surface at $z=0$.

$${\overline{E}}_{ref}\left(\overline{r}\right)={\overline{E}}_{0}R\left(\theta \right)e^{i{k}_{ix}x+i{k}_{iy}y+i{k}_{iz}z}$$

(120)

where $R\left(\theta \right)$ is the Fresnel reflection coefficient at angle ${\theta }_i$. The reflection coefficient is $R\left({\theta }_i\right)=R^{TM}\left({\theta }_i\right)$ for TM and $R\left({\theta }_i\right)=R^{TE}\left({\theta }_i\right)$ for TE. By periodic boundary conditions, we mean that the horizontal domain is divided into periodic regions and the permittivities are repeated periodically. The (0,0) period is

$$\begin{array}{c}-\frac{L_x}{2}\le x\le \frac{L_x}{2}\\ -\frac{L_y}{2}\le y\le \frac{L_y}{2}\end{array}$$

(121)

Then, the cells are labeled by integers $p,q$ with p, q = 0, ±1, ±2, … Let

$${\overline{\rho }}_{pq}=p{L}_x\widehat{x}+q{L}_y\widehat{y}$$

(122)

The periodic permittivities are

$${\epsilon }_m\left(\overline{r}+{\overline{\rho }}_{pq}\right)={\epsilon }_m\left(\overline{r}\right)$$

(123)

The boundary value problem is 3D electromagnetic wave scattering, and the periodicity is 2D (labeled as a 3D2D problem). The electric field obeys the Bloch condition

$$\overline{E}\left(\overline{r}+{\overline{\rho }}_{pq}\right)=\overline{E}\left(\overline{r}\right)e^{i{\overline{k}}_{i\rho }·{\overline{\rho }}_{pq}}$$

(124)

By making use of periodic Green’s functions, The VIE is converted into one period (0,0) for the inhomogeneous region

$$\begin{gathered} \overline{E}\left(\overline{r}\right)={\overline{E}}_{inc}\left(\overline{r}\right)+k^2{\int }_0^{d}d{z}^{\prime }{\int }_{-\frac{L_x}{2}}^{\frac{L_x}{2}}d{x}^{\prime }{\int }_{-\frac{L_y}{2}}^{\frac{L_y}{2}}d{y}^{\prime }{\stackrel{=}{G}}_P\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right) \left(\frac{{\epsilon }_m\left(\overline{r}\right)}{\epsilon }-1\right)\overline{E}\left(\overline{r}\right) \\ ;-\frac{L_x}{2}\le x\le \frac{L_x}{2},-\frac{L_y}{2}\le y\le \frac{L_y}{2},0\le z\le d \end{gathered}$$

(125)

where $\overline{r}$ and ${\overline{r}}^{\prime }$ are in the inhomogeneous region of period (0,0). The ${\stackrel{=}{G}}_P\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)$ is the half space periodic Dyadic Green’s function. It is decomposed into the free space part ${\stackrel{=}{G}}_P^{\left(0\right)}$ and the reflection part ${\stackrel{=}{G}}_P^{\left(R\right)}$.

$${\stackrel{=}{G}}_P\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)={\stackrel{=}{G}}_P^{\left(0\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)+{\stackrel{=}{G}}_P^{\left(R\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)$$

(126)

By applying the discrete dipole approximation (DDA) with small cube size $\Delta V$ and taking into account the singularity of the dyadic Green’s function of the free space response ${\stackrel{=}{G}}_P^{\left(0\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right),$ we obtain the matrix equation.

$${\overline{p}}_i={\stackrel{=}{\alpha }}_i\left[{\overline{E}}_{inc}\left({\overline{r}}_i\right)+{\overline{E}}_{ref}\left({\overline{r}}_i\right)+\frac{k^2}{\epsilon }{\displaystyle \sum }_{j\ne i} {\stackrel{=}{G}}_P^{\left(0\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)·{\overline{p}}_j+\frac{k^2}{\epsilon }{\displaystyle \sum }_j {\stackrel{=}{G}}_P^{\left(R\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)·{\overline{p}}_j\right]$$

(127)

where

$${\overline{p}}_i=\Delta V\epsilon \left(\frac{{\epsilon }_m\left({\overline{r}}_i\right)}{\epsilon }-1\right)\overline{E}\left({\overline{r}}_i\right)$$

(128)

is the polarization of the cube.

$${\stackrel{=}{\alpha }}_i=\Delta V\epsilon \left(\frac{{\epsilon }_m\left({\overline{r}}_i\right)}{\epsilon }-1\right){\left[\stackrel{=}{I}-k^2\left(\frac{{\epsilon }_m\left({\overline{r}}_i\right)}{\epsilon }-1\right)\stackrel{=}{S}\right]}^{-1}$$

(129)

is the polarizability of the i-th cube and $\stackrel{=}{S}$ is the singular integral over the self cube.

$$\stackrel{=}{S}={\int }_{V_C^{\left(0\right)}}^{}{\stackrel{=}{G}}_P^{\left(0\right)}\left(\overline{0},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)d{\overline{r}}^{\prime }$$

(130)

The expression for half space dyadic Green’s functions are the non-periodic ones in Tsang et al. volume 1 with the wave vectors replaced by Bloch vectors [28]. The free space part is

$${\stackrel{=}{G}}_P^{\left(0\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)=\frac{i}{2\mathsf{\Omega }}{\displaystyle \sum }_{\alpha }\frac{1}{k_{z\alpha }}\left[{\widehat{e}}_{\alpha }\left(k_{z\alpha }\right){\widehat{e}}_{\alpha }\left(k_{z\alpha }\right)+{\widehat{h}}_{\alpha }\left(k_{z\alpha }\right){\widehat{h}}_{\alpha }\left(k_{z\alpha }\right)\right]e^{i{\overline{k}}_{\rho }·\left(\overline{r}-{\overline{r}}^{\prime }\right)} for z>z^{\prime }$$

(131)

$${\stackrel{=}{G}}_P^{\left(0\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)=\frac{i}{2\mathsf{\Omega }}{\displaystyle \sum }_{\alpha }\frac{1}{k_{z\alpha }}\left[{\widehat{e}}_{\alpha }\left(-k_{z\alpha }\right){\widehat{e}}_{\alpha }\left(-k_{z\alpha }\right)+{\widehat{h}}_{\alpha }\left(-k_{z\alpha }\right){\widehat{h}}_{\alpha }\left(-k_{z\alpha }\right)\right]e^{i{\overline{K}}_{\rho }·\left(\overline{r}-{\overline{r}}^{\prime }\right)} for z<z^{\prime }$$

(132)

and the reflected part is

$${\stackrel{=}{G}}_P^{\left(R\right)}\left(\overline{r},{\overline{r}}^{\prime },{\overline{k}}_{i\rho }\right)=\frac{i}{2\mathsf{\Omega }}{\displaystyle \sum }_{\alpha }\frac{1}{k_{z\alpha }}\left[R^{TE}{\widehat{e}}_{\alpha }\left(k_{z\alpha }\right){\widehat{e}}_{\alpha }\left(-k_{z\alpha }\right)+R^{TM}{\widehat{h}}_{\alpha }\left(k_{z\alpha }\right){\widehat{h}}_{\alpha }\left(-k_{z\alpha }\right)\right]e^{i{\overline{k}}_{\rho }·\left(\overline{r}\right)}{e}^{-i{\overline{K}}_{\rho }·{\overline{r}}^{\prime }} for z<z^{\prime }$$

(133)

where ${\overline{k}}_{\rho }$ and ${\overline{K}}_{\rho }$ are respectively the upward Bloch vector and downward Bloch vector.

For 2D periodicity,$\mathsf{\Omega }=L_x{L}_y$, the reciprocal lattice vectors are

$${\overline{K}}_{mn}=m\frac{2\pi }{L_x}\widehat{x}+n\frac{2\pi }{L_y}\widehat{y}$$

(134)

where $m,n=0,\pm 1,\pm 2,\dots$. The Bloch vectors are

$${\overline{k}}_{imn}={\overline{k}}_{i\rho }+{\overline{K}}_{mn}$$

(135)

so that

$$k_{zmn}=\sqrt{k^2-{\left|{\overline{k}}_{imn}\right|}^2}$$

(136)

We also use $\alpha$ to represent the double index $\left(m,n\right)$ so that ${\sum }_{\alpha }={\sum }_{m,n}$, $k_{z\alpha }=k_{zmn}$.

The upward propagation Bloch vector is

$${\overline{k}}_{\alpha }={\overline{k}}_{imn}+{\overline{k}}_{zmn}\widehat{z}$$

(137)

The downward propagation Bloch vector is

$${\overline{K}}_{\alpha }={\overline{k}}_{imn}-{\overline{k}}_{zmn}\widehat{z}$$

(138)

The definitions of TE and TM polarization vector $\widehat{e}$ and $\widehat{h}$ follow that of Tsang et al. volume 1 and are

$${\widehat{e}}_{\alpha }\left(k_{z\alpha }\right)={\widehat{e}}_{\alpha }\left(-k_{z\alpha }\right)=\frac{1}{\left|{\overline{k}}_{imn}\right|}\left(k_{y\alpha }\widehat{x}-k_{x\alpha }\widehat{y}\right)$$

(139)

$${\widehat{h}}_{\alpha }\left(k_{z\alpha }\right)=-\frac{k_{z\alpha }}{k\left|{\overline{k}}_{imn}\right|}\left(k_{x\alpha }\widehat{x}+k_{y\alpha }\widehat{y}\right)+\frac{\left|{\overline{k}}_{imn}\right|}{k}\widehat{z}$$

(140)

$${\widehat{h}}_{\alpha }\left(-k_{z\alpha }\right)=\frac{k_{z\alpha }}{k\left|{\overline{k}}_{imn}\right|}\left(k_{x\alpha }\widehat{x}+k_{y\alpha }\widehat{y}\right)+\frac{\left|{\overline{k}}_{imn}\right|}{k}\widehat{z}$$

(141)

4.2. Modeling and Estimation of the Roughness

In NASA’s Soil Moisture Active Passive (SMAP) mission, the objective is to retrieve soil moisture globally using radar and radiometer observations. For the radar retrieval part, the lookup table (LUT) is applied according to landcover and vegetation type [52]. Kim et al. applied a time-series approach to retrieve soil moisture [29]. For a given pixel, the rms height is assumed constant, so it is not only retrieving soil moisture but also rms height [29]. The retrieved soil moisture and rms height are provided in SMAP L3 Radar Global Daily 3 km EASE-Grid Soil Moisture products through the National Snow & Ice Data Center (NSIDC) [66]. Although SMAP’s radar was only active for 3 months, the amount of retrieved global rms height map provides the basis for our analysis here. Figure 15a plots the rms heights obtained over North America based on the average of rms height maps from the 3 months products. The Figure shows most locations to have rms heights of 2 cm or less while mountainous regions show rms heights from 4 to 6 cm. A histogram of these rms heights for North America is further provided in Figure 15b. In rough-surface scattering at the Ku-band, a upper limit of $kh$ = 2 or 3 was used, which missed the fact that the $kh$ upper limit is actually 20, which is ten times larger than $kh$ = 2. Next, we consider the scattering problem taking into account those fine features at higher frequency.

Consider an incident plane wave on a rough surface as shown in Figure 16. The direction of the incident wave is ${\overline{k}}_i=k{\widehat{k}}_i=k_{ix}\widehat{x}+k_{iy}\widehat{y}-k_{iz}\widehat{z}$, where $k_{ix}=k \sin {\theta }_i\cos {\varphi }_i$, $k_{iy}=k\sin {\theta }_i\sin {\varphi }_i$, and $k_{iz}=k\cos {\theta }_i$. ${\theta }_i$ and ${\varphi }_i$ are the incident elevation and azimuth angle, respectively. The scattered wave ${\overline{E}}_s$ is with scattered direction of ${\overline{k}}_s=k{\widehat{k}}_s=k_{sx}\widehat{x}+k_{sy}\widehat{y}+k_{sz}\widehat{z}$, where $k_{sx}=k\sin {\theta }_s\cos {\varphi }_s$, $k_{sy}=k\sin {\theta }_s\sin {\varphi }_s$, and $k_{sz}=k\cos {\theta }_s$. ${\theta }_s$ and ${\varphi }_s$ are the scattering elevation and azimuth angle, respectively.

The roughness profile $f\left(x,y\right)$ is given by

$$\begin{array}{c}f\left(x,y\right)=f_1\left(x,y\right)+f_2\left(x,y\right)\\ h=\sqrt{h_1^2+h_2^2}\\ h^{2}C\left(\rho \right)=h_1^2{C}_1\left(\rho \right)+h_2^2{C}_2\left(\rho \right)\end{array}$$

(142)

where $f_1\left(x,y\right)$ is the finer roughness in millimeter scale with rms height $h_1$, correlation length $l_1$, and correlation function $C_1\left(\rho \right)$. The $f_2\left(x,y\right)$ is roughness in centimeters with rms height $h_2$, correlation length $l_2$, and correlation function $C_2\left(\rho \right)$. For the combined $f_1\left(x,y\right)$ and $f_2\left(x,y\right)$ in $f\left(x,y\right)$, the rms height $h$ and correlation function $C(\rho$) are given by Equation (141). $l$ is the correlation length of $f\left(x,y\right)$ that $C\left(l\right)=exp\left(-1\right)$.

Figure 17 plots the exponential term of the Kirchoff integral in the backscattering direction as a function of $\rho$ from the L- to the Ku-band. Results indicate that the exponential term has effects on backscatter for $\rho$ more than 10 cm for the L-band. For frequencies above the X-band, the exponential term has effects for $\rho$ less than 3 cm, which is less than the correlation length of $f_2$ but comparable with the correlation length of $f_1$. The millimeter-scale roughness $f_1$ has effects for high frequency. To investigate the scattering of rough surfaces at the X- and Ku-bands, it is important to include fine-scale roughness $f_1$ superimposed on $f_2$, providing a two-scale roughness for the present simulations.

In the VIE-NMM3D approach, we modeled the rough-surface scattering using a Monte Carlo simulation in which the scattering from each surface realization was computed using a discrete dipole approximation (DDA) to the volume integral [67,68,69].

The rough surface realizations generated are stationary Gaussian random processes having horizontal dimensions $L\times L$. For each realization, the minimum of the rough surface profile, $f_{min}=MIN\left\{f\left(x,y\right)\right\}$, is defined as $z=0$, and the “depth” of the simulation domain must satisfy $d>f_{max}-f_{min}$, where $f_{max}=MAX\left\{f\left(x,y\right)\right\}$. Note this definition of the coordinate system results in a mean surface height > 0 and results in phase offsets for each realization. This is immaterial for the purposes of this work because only the incoherent backscattering coefficient $\gamma ({\theta }_s$, ${\varphi }_s)$ at oblique incidence is examined. The DDA simulation domain of dimension $L\times L\times d$ is then discretized into small cubes whose dimensions are small compared with the wavelengths of interest and to the roughness scales of interest. Cubes within the soil region are modeled as having the permittivity of soil while those in free space are modeled as having the permittivity of free space. Because these cubes exist in the presence of the soil half space z < 0, a half space Green’s function was used to compute the radiation from each cube with $z=0$ representing the half space interface. Periodic boundary conditions were also applied in both horizontal directions to reduce the effect of the horizontal boundaries. This has the effect of discretizing the scattered field into a set of Floquet mode plane waves. However, it should be expected that surface periods $L$ much greater than the wavelengths of interest will produce Floquet mode amplitudes that approximate the scattering in those directions of an infinite continuous rough surface when observed in the same scattering direction.

4.3. The Frequency and Roughness Responses of the Backscattering from the Rough Surface

In this section, we provide a few simulation results [46]. More simulation results can be found in the PhD thesis of Jiyue Zhu [47]. In Figure 18, the backscattering is plotted as a function of kh up to kh = 6. The ratio of correlation length to rms height was 7. The simulated frequency was 9.6 GHz. Figure 18 shows that the DDA simulations are in good agreement with previous MoM full-wave simulations LUT up to kh = 1.5 [70]. The backscattering of DDA simulations first increased with kh and then saturated around kh = 2.5. The physical basis behind such phenomena is that the magnitude of waves reflected back to the backscattering direction increase with the roughness and then finally reaches a maximum value (which is corresponding to the geometric optics (GO) realm [39]). Saturation was also postulated by Oh’s empirical model without experimental data in the large kh regime [59].

In Figure 19, the backscattering is plotted as a function of $kh$ for rough surfaces with “limited correlation lengths”. The rms heights (0.3 cm, 0.6 cm, 0.9 cm, 1.2 cm, 1.5 cm, and 1.8 cm) and correlation lengths (2.64 cm, 5.04 cm, 5.85 cm, 6.3 cm, 6.75 cm, and 7.2 cm) are from field measurements in San Luis Valley (SLV), CO, USA. A laser range finder mounted on a horizontal bubble level supported by a tripod at each end was used to measure small-scale surface roughness along the 1 m baseline of the level at multiple locations [42]. The simulated frequency was 13.6 GHz. We considered a snow-covered, rough soil surface. Backscattering first increased with $kh$ to reach the peak around $kh=2.1$ and then decreased. Such a decrease is caused by the decrease of the ratio of correlation length and rms height. In the GO limit, smaller ratio yields a weaker backscatter.

In Figure 20, we show the VV backscattering as a function of rms heights at C-, X-, and Ku-bands. The figures show saturation effects, meaning that rough-surface scattering saturates at large rms heights (~3–6 cm) and at high frequencies. The new results of $kh$ up to 15 are useful for studying rough surface radar backscattering at the X- and Ku- bands for snow/soil interfaces.

4.4. Using the Retrieved rms Height from UAVSAR L-Band Data to Simulate Backscattering at X- and Ku-Bands

Interaction of Radar Waves with the Ground Surface Beneath the Snowpack

Figure 21 shows radar backscattering from a snow layer. The interaction of radar with snow overlying a rough soil surface involves refraction and scattering from the air/snow interface. In Figure 21, we show the incident plane wave at incident angle ${\theta }_i$. Based on Snell’s law, the transmitted angle is ${\theta }_t$ so that the incident angle on the rough soil surface is ${\theta }_t$. With snow relative permittivity of approximately 1.44, the refracted angle ${\theta }_t$ was 32.4 degrees, which corresponds to the 40 degrees incidence angle of ${\theta }_i$ in air. Because of the magnitude of dielectric contrasts, the contributions of rough-surface scattering are from the snow–soil rough interface and not from the air/snow interface. The air/snow interface has stronger scattering contribution for wet snow which is outside the domain of the X-band and Ku-band volume scattering approach for SWE retrieval.

The volume scattering method at the X-band the and Ku-band has been proposed to retrieve SWE for dry snow. For such a case, the rough surface contribution arises from the rough snow–soil surface below the snow layer. It is proposed in this paper that L-band and C-band radar signals can be used to retrieve both the rms height and the moisture of the soil. Since microwave scattering depends on permittivity. The retrieval of soil moisture covers the various soil conditions of soil texture, freeze–thaw, etc. For wet snow, a rough-surface contribution arises from the air/snow surface. However, the volume scattering method is not proposed for wet snow.

Rough-surface scattering from the snow/soil interface contributes to radar observations as represented by ${\sigma }_{rs}exp\left(-2\tau sec{\theta }_t\right)$. The term is affected by the rough soil-surface scattering ${\sigma }_{rs}$ and the attenuation through snow factor of $exp\left(-2\tau sec{\theta }_t\right)$. $\sigma$ refers to either VV or HH. $rs$ stands for rough surface and $\tau$ is optical thickness. The rough soil-surface scattering contribution is not related to SWE and it should be removed when retrieving SWE. The subtraction of surface scattering has been used to improve the accuracy of SWE retrieval [17,24]. The approach of removing ${\sigma }_{rs}exp\left(-2\tau sec{\theta }_t\right)$ consists of using a combination of data and electromagnetic models. The removal of surface scattering is a significant part of the retrieval algorithm that is discussed later in this section.

To estimate rough-surface scattering at the X-band and Ku-band, there are two approaches, (a) and (b). Approach (a) is as described in ref. [18]. In approach (a), radar observations at the X- and Ku-bands at a specific location before the snowfall are used to estimate the surface backscattering [18]. Such an approach neglects changes in soil properties and background land-cover during the snow-on season. In approach (b), the estimation of rough-surface scattering, ${\sigma }_{rs}$, at the X-band and Ku-band consists of using a combination of measurement data and electromagnetic models. The measurement data include backscattering data at the L-band, C-band, X-band and Ku-band before and after snow. The electromagnetic model of rough-surface scattering results is based on using look-up table (LUT) results of Numerical Maxwell 3D model of full-wave simulations of surface scattering from the L-band to the Ku-band. In this section, we illustrate approach (b) using L band UAVSAR data. We used a two-step procedure. In the first step of approach (b), we used co-polarized radar time series observations at the L- and C-bands, which have much larger surface scattering than volume scattering, to estimate soil permittivity and surface roughness. The use of the C-band time series together with the L-band time series enhances the existing L-band algorithm in retrieving rms height and soil moisture. Using L-band and C-band data before and after snow fall together with LUT of Maxwell equations for the L-band and C-band, we retrieve the rough surface rms heights and soil moistures. The retrieval is carried out for both before snow fall and after snow fall. After snow fall, Snell’s law is used to adjust the incident angle at the snow–soil interface to account for the refraction angle at the air/snow interface. In step 2 of approach (b), the rms heights and soil permittivities retrieved in step 1 are used in the NMM3D LUT to obtain the model backscattering results of ${\sigma }_{rs}$ at the X- and Ku-bands for the same roughness and soil moisture parameters. Approach (b) is applicable even when there are changes of rms heights and soil moistures during the snow season. Approach (b) assumes the availability of matchup spaceborne L- and/or C-band SAR observations with revisits of 10 days. The revisits provided by the Sentinel-1 and NISAR systems, as well as future proposed continuation missions, suggest that such datasets are likely to be available during the time frame of a future snow observing mission.

We next give an example of approach (b) using UAVSAR L band data to demonstrate the retrieval of rms heights and soil permittivity. We consider, as in step 1, using L-band UAVSAR radar full polarization observations under snow-on conditions using the NMM3D LUT [27]. The UAVSAR dataset was collected from February to March in 2017, in the SnowEx 2017 campaign using five flights over the Grand Mesa region in Colorado, United States. In situ soil moisture measurements were also collected throughout 2017 from an installed meteorological observation station. We apply the time series retrieval algorithm developed for the SMAP mission at the station location, that is, at the point-scale [52,56]. From the retrieved soil permittivity, the soil moisture was derived using Mironov’s empirical model [57]. The comparison of retrieved and measured soil moisture is shown in Figure 22a. In addition to retrieving soil moisture, the rms height was also retrieved from the time series and was at a single value of 1.9 cm. The cost function for retrieving soil moisture and rms height is as follows. By using the L-band rough surface lookup table, we used VV and HH to find the minimum cost function against the observed backscatters. The rms height and time series soil moisture were retrieved at the same time. The retrieved rms height from L-band is assumed to be invariant with time and subsequently can be applied for rough-surface scattering at higher frequencies.

$$C_{min}=min\left\{{\displaystyle \sum }_n\left[{\left({\sigma }_{VV}^{LUT}\left(h,{ϵ}_n\right)-{\sigma }_{VV}^{obs}\right)}^2+{\left({\sigma }_{HH}^{LUT}\left(h,{ϵ}_n\right)-{\sigma }_{HH}^{obs}\right)}^2\right]\right\}$$

(143)

The above description is the rationale of the Kim et al., 2017, algorithm showing that soil moisture changes but not rms height. Wet and frozen soils are assumed to have the same rms height [52].

The retrieval soil moisture is in good agreement, as shown in Figure 22a, with the measured in situ soil moisture for a period of 8 weeks from 6 February 2017 to 31 March 2017. The agreement achieved has a root mean square error (RMSE) of 0.047 ${\mathrm{m}}^3/{\mathrm{m}}^3$, a correlation of 0.95, and a bias of 0.039 ${\mathrm{m}}^3/{\mathrm{m}}^3$. Based on the measurements, the soil moistures at this site were relatively constant during the dry snow season (6 February to 8 March 2017).

To complete approach (b), we carried out step 2 and applied the retrieved rms height of 1.9 cm, and soil properties from Figure 22a to calculate the surface scattering contributions with snow attenuation, ${\sigma }_{rs}exp\left(-2\tau sec{\theta }_t\right)$, at the X band of 9.6 GHz, the low Ku band of 13.4 Ghz, and the high Ku-band of 17.2 Ghz. These three frequencies were considered for the proposed satellite snow water equivalent (SWE) missions. Figure 22b shows surface scattering including snow attenuation. The results in Figure 22b show that the rough soil-surface scattering contribution, including snow attenuation, was around −12 dB at the X band and decreased to −14 dB at the high Ku band of 17.2 GHz. Higher frequencies such as at the Ku band typically experience higher volume scattering and greater attenuation from the snow layer of the surface scattering contributions. Continued studies are required to improve and validate the approach (b). These studies include extending NMM3D surface modeling studies and LUT into cases with rms heights of four wavelengths or more so that the LUTs can be applied to the Ku-band of 17.2 GHz up to 7 cm of rms heights. Unlike volume scattering of snow, rough-surface scattering has a stronger dependence on incidence angle. Thus, the effects of topographical slopes that cause changes in incidence angles, particularly in mountainous regions, should be included in the retrieval. The sample size in Figure 22a,b is small. However, with the availability of L band and C band data, full-wave simulations of rough soil-surface scattering from the L band to the Ku band up to $kh=15$ allow retrievals of rms heights of the soil surfaces and the soil-surface scattering at the X band and Ku band can be determined.

Our current work is computing look-up tables for rough soil-surface scattering for kh up to 15. We are also continuing investigating various CEM techniques of rough-surface scattering in the large kh domain.

5. Conclusions

In this paper, we reviewed our recent work in three topics of the theory of microwave remote sensing. The effects of vegetation and forests are important topics of microwave remote sensing of soil moisture and the snow water equivalent. The hybrid method of full-wave simulations showed results that are significantly different from that of RTE and DBA. As remote sensing of soil moisture extends from the P-band to the X-band while remote sensing of snow water equivalent covers the P-band to the Ku-band, such studies using full-wave simulations will provide new solutions different from that of RTE/DBA. Signals of opportunity offer new efficient ways of remote sensing at P-band and L-band. The theory of microwave remote sensing at bistatic scattering close to the specular direction requires modeling techniques different from radar backscattering and radiometry. For remote sensing of land surfaces at high frequencies, such as the X-band Ku-band, the theory of microwave remote sensing full-wave simulations needs to be carried out to much larger kh values than before. The study of multiscale surfaces and their effects at large kh continue to be a challenge.