Spectral Properties of Bistatic Radar Signals using the Ray Tracing Technique and a Facet Approach

1. Introduction

Remote sensing using bistatic radar configurations, where the transmitter and receiver are spatially separated, uses scattering off a surface to measure surface characteristics such as roughness, reflectivity, and dielectric constants. Bistatic radar for the study of statistical surface characteristics was first applied to the lunar surface [1] in the 1960s where the spacecraft (S/C) and a terrestrial ground station served as the radar transmitter and receiver, respectively. In addition to the moon, bistatic radar on space missions has been used to study the surfaces of extra-terrestrial bodies in the Solar System [2,3,4]. Bistatic radar used to study planetary surfaces includes not only orbital-Earth configurations but also the use of multiple terrestrial ground stations. Terrestrial bistatic radar observations have provided crucial information for space mission planning [5]. Additionally, a bistatic configuration between a lander and a spacecraft has been used to provide surface roughness data by the Huygens lander on Titan for the NASA/ASI/ESA Cassini-Huygens mission [6] and by the Philea lander on comet 67P/Churyumov-Gerasimenko (67P/C-G) for the ESA Rosetta mission [7].

The relationship between the characteristics of the received radar signal and the statistics of a rough surface, such as the root-mean-square (rms) slope, is used to constrain the roughness of the surface at scales of the incident wavelength (cm). Polarization of the bistatic forward scattering can be used to determine the dielectric constant of the surfaces as well detect the presence of water-ice on lunar and planetary surfaces using same sense circular polarization of the signal. The roughness of geological surfaces at microwave wavelengths is important to planetary geology and enables identification of geological unit types, geomorphological changes over time, geophysical processes, and the presence of volatiles on the Moon [8,9], Venus [10,11], Mars [12,13,14], Titan [15,16], Pluto [17], asteroid Vesta [18] and comet 67P/C-G [19,20] among others. Characterizing surface roughness also provides crucial data for planning space missions that include landers and rovers as well as surface sampling. Future missions such as ESA’s Jupiter Icy Moons Explorer (JUICE) will utilize bistatic radar to study the surface properties on Ganymede [21]. In addition, future bistatic experiments could include micro satellite networks instead of a S/C in an orbital bistatic configuration as well as a possible unmanned aerial vehicle (UAV) to rover configuration like the Mars Helicopter demonstration (Ingenuity) currently being tested with the Perseverance rover on Mars [22,23,24,25].

Typically, bistatic radar data uses a downlink configuration between the S/C and Earth, where the transmitter is on-board the S/C as part of the telecommunications system and the terrestrial ground station is the receiver. However, for the New Horizon mission to Pluto, an uplink configuration, where the receiver is on-board the S/C and the terrestrial ground station transmits commands after the scattered signal is received, provided the necessary higher signal-to-noise ratio given the extreme distance involved as well as the high speed of the S/C during its flyby at Pluto [26].

2. Background and Goals

Orbital bistatic radar experiments transmit continuous sinusoidal signals from the S/C moving with a certain velocity with respect to the surface at a particular incident angle. As a first approximation, the surface of a planet, which at microwave wavelengths can be assumed to behave as a gently undulating surface, is considered to scatter the transmitted radar signal primarily in a specular manner [2]. The angle of incidence equals the angle of reflection, as shown in Figure 1 assuming the azimuthal angle between the plane of incidence and observation is zero. The surface is broken down into various planar scattering elements, each tangential to the existing surface, and the received signal is considered to be a summation of specular reflections from individual scattering elements or facets that are oriented to receive the direct energy of the transmitted signal. This model produces quasi-specular scattering due to randomly oriented perfectly flat facets. Diffuse scattering due to wavelength-scale surface features and sub-surface scattering can be part of the total scattering from a planetary surface, but are neglected in this analysis. The quasi-specular model is based on the Kirchhoff or tangent-plane approximation and the stationary phase or high-frequency approximation to describe the surface (e.g., [27,28]). Kirchhoff’s approximation assumes that the spatial scale of the surface structure is larger than the incident wavelength and the high-frequency approximation assumes that only those surface facets whose slope enables a specular reflection of the transmitted ray contribute to the radar cross-section. The size of the surface facet should approximate at least one wavelength, thus reducing diffraction effects, while the dimension of the effective Fresnel zone is the upper limit [29,30,31]. In this case, where the electromagnetic wavelength is much smaller than the surface radii of curvature, geometrical optics approximations are valid.

The typical incidence angles for planetary bistatic orbital configurations have ranged between 40 and 70 degrees. While analytical models relying on quasi-specular scattering are not considered reliable for high and low (grazing) incidence angles [32,33] due to shadowing, diffraction, and multiple scattering effects, in Ref [34] Gunnarsdottir found data produced in multiple high incidence angle bistatic radar experiments for Mars in good agreement with MOLA data from the Mars Odyssey mission.

The power of the received signal emitted from the surface facet reflects the characteristics of the transmitter and receiver antennas and of the surface facets as shown in Equation (1)

$$\mathit{d}{\mathit{P}}_{\mathit{R}}=\left(\frac{{\mathit{P}}_{\mathit{T}}{\mathit{G}}_{\mathit{T}}}{4\mathit{\pi }{\mathit{R}}_{\mathit{T}}^2}\right)\mathit{\sigma }\left(\frac{{\mathit{A}}_{\mathit{R}}}{4\mathit{\pi }{\mathit{R}}_{\mathit{R}}^2}\right)$$

(1)

where dP_R is the incremental received signal power, P_T is the transmitted power, G_T is the gain of the transmitter antenna, A_R is the effective area of the receiver antenna, while R_T and R_R are the distances to and from the surface facet from the transmitter and the receiver, respectively. Finally, σ is the bistatic radar cross section (RCS) of the individual surface facet and is written as:

$$\sigma ={\sigma }_{0}dS$$

(2)

and is a function of the slope or tilt angle of each element. The slope indicates how closely the facet aligns with the case of specular reflection. The smaller the slope, the more closely the facet aligns to specular reflection and the greater its effect on the received signal.

Scattering laws based on Gaussian and exponential forms have been used to investigate scattering off a quasi-specular surface. However, the most commonly accepted scattering law for the bistatic RCS is Hagfor’s law which is based on a Gaussian height distribution and an exponential autocorrelation function for surface length and frequency components. Hagfor’s law is also used with other surface models, such as fractal surfaces, where it is considered a special case, along with Gaussian and exponential scattering laws, in fractal-based scattering laws [35].

The S/C produces a radar swath of data as it moves across the surface and causes the received signal to be spectrally broadened due to variations in the Doppler shift at different points within the footprint and along the surface track. The surface facets, which are also moving with a planetary velocity, produce a bistatic Doppler shift relative to the transmitter and receiver. The degree to which the received signal is spectrally broadened is proportional to the surface roughness:

$$B=4\left(ln2{)}^{\frac{1}{2}}\right(V\mathsf{\zeta }/)\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\phi }$$

(3)

where B is the half-power bandwidth of the received signal, V is the velocity of the specular point across the surface, ζ is the rms surface roughness and φ is the incidence angle [1]. Thus, a wider received signal bandwidth is a result of a rougher surface and/or a higher velocity S/C. However, if the surface becomes too rough, some of the facet scattering will be outside of the illumination pattern of the transmitting antenna thus reducing the power of the received signal (beam-limitation).

In this paper, we investigate how the variability of surface roughness affects the spectral broadening of the received radar signal using a quasi-specular model with both a physical optics approach and a geometrical optics approach by means of ray-tracing. The surface is modeled as a collection of local facets with Gaussian surface statistics and includes both isotropic and anisotropic surfaces. While numerical methods such as the Integral Equations Model (IEM) and the Methods of Moments (MoM) can solve Maxwell’s equations numerically and provide an exact solution, they are computationally expensive. Ray-tracing provides an approximate solution to Maxwell’s equations in the high-frequency regime. Physical optics is used to generate the surface model and to obtain Doppler shifts and subsequent changes to the received signal bandwidth. We will also investigate the effect of minimum ray density on the results of this model.

This numerical model has been developed using the MATLAB interface for SPICE (Spacecraft, Planet, Instrument, “C-matrix”, Events) open access code developed by NASA’s NAIF (Navigation and Ancillary Information Facility). The numerical model can be seen as consisting of three basic steps; model parameter generation, ray propagation and signal analysis. Model generation includes the generation of basic parameters necessary to determine the received signal characteristics such as positions of the transmitting and receiving antenna (satellite and terrestrial ground station) and planetary target, antenna footprint coverage, beam angle as well as generation of the rough faceted Gaussian surface. This data can be self-generated or else pulled from SPICE mission databases. Electromagnetic wave propagation in this model between the satellite, planetary surface and receiving antenna will be approximated by ray-tracing technique. After this, the Doppler shift of the signal and the signal bandwidth can be calculated so that the relationship between the surface roughness and the spectral broadening can be ascertained.

3. Surface Properties from Quasi-Specular Radar Backscatter

We present the model used to generate the scattered field from a rough surface illuminated by a transmitting antenna simulating wave propagation using ray tracing algorithms. The discussion will first focus on the ray tracing algorithm used to determine the position of the reflected ray on a smooth surface and then a rough Gaussian surface. To determine the direction of the reflected ray, the rough surface will be decomposed into triangular micro facets. The receipt of the signals will be determined using both circular region-based, and angle-based receiving conditions. Finally, the characteristics of the signal will be examined. The parameters used in this model are presented in Table 1. The positions of the components of this model (satellite transmitting antenna, planetary surface, and receiving antenna) can be self-determined or else downloaded from the SPICE libraries for a specific space mission. The planetary surface considered here is a perfectly conducting spherical one with a Gaussian height distribution.

Using Kirchhoff’s and the high-frequency assumptions, we assume the local points on the planetary surface to be a tangent plane approximation to the actual surface and we assume the roughness is smaller than the wavelength thus the scattering of the local surface is equivalent to the specular reflection of a smooth surface. In this case, the direction of the reflected wave will obey Huygens’ principle such that the incident and reflected (in this case scattered) angle and points are equal and the incident ray, reflected ray, and surface normal all lie in the same plane as seen in Figure 1. The magnitude of the reflected ray can be determined using Fresnel’s reflection coefficients. For this model, we will assume all the Fresnel reflection coefficients are unity. Thus, to determine the reflected ray, the incident direction, the incident surface location, and the normal to the surface at the interception point are necessary. As shown in Figure 2, the position of the satellite at an epoch time $t_i$, is designated as $P_{t_i}^S$, the position of the receiver antenna at epoch time is $P_{t_i}^R$, and $P_{t_i}^P$ is the centre position of the planetary surface. The position of the central ray intercepting the surface at epoch time is $P_{t_i}^O$, which can be obtained by referring to $P_{t_i}^S$, $P_{t_i}^R$ and $P_{t_i}^P$ or else by using the SPICE software functions.

While illuminating the surface with a large number of radio rays would lead to a higher resolution it would also lead to slower ray-tracing computation due to the increased need for ray interception checking. To limit the ray tracing computational time, the model initially sets the number of rays equal to the number of surface elements, ${(N+1)}^2$, with N being equal to the number of divisions along each side of the half-power beamwidth, which generates a square grid of $(N+1)\times (N+1)$ elements.

3.1. The Simulation of Reflected Signals

To trace the ray through the model, we assume that the satellite transmits a bundle of rays at a certain half-beam angle, α. For the sake of the algorithm, the beam footprint can be assumed to consist of a central curve along with a left and right curve. This algorithm processes the ray interception points along the central curve first and then the lateral points going left to right along the footprint. We use the Law of Reflection which states that the incident ray, reflected ray, and surface normal all exist in the same plane and that the incident angle equals the reflected angle. In short, a point on the central curve exists within a plane defined by that point and the positions of the satellite and planet and this plane has a normal which goes through the point and the center planet position as seen in Figure 2. Thus, the reflected ray position and angle can be found rotating the point around the normal. While in reality the reflected ray position equals the surface interception ray position, for the sake of this algorithm it is defined as another point. This new point, similar to the initial point, is defined as existing within a plane defined by itself and the positions of the satellite and planet with a normal going through itself and the planet’s position. At the local ray interception position, on the central curve, $P_{t_i}^O$, the reflected ray position and direction can be obtained by the rotation of $P_{t_i}^O$ around the normal to the plane on the central curve determined by the surface interception point, $P_{t_i}^O$ along with $P_{t_i}^S$ and $P_{t_i}^P$, denoted as $\mathsf{\Gamma }(P_{t_i}^S,P_{t_i}^P,P_{t_i}^O)$. The normalized vector, denoted as ${\rm Z}(P_{t_i}^S,P_{t_i}^P,P_{t_i}^O)$, that $P_{t_i,O}$ rotates around, is the normal to this plane, $\mathsf{\Gamma }(P_{t_i}^S,P_{t_i}^P,P_{t_i}^O)$, and goes through point $P_{t_i}^S$. The range of the rotation angle, $\beta$, is between $[-\alpha ,\alpha ]$, and the $P_{t_i}^O$ rotated by $\beta$ is denoted as $P_{t_i,\beta }^O$, shown in Figure 2; this point represents the reflected/scattered point. This point, similar to the initial ray interception point, $P_{t_i}^O$, is also defined as existing within a plane defined by the reflected point, $P_{t_i,\beta }^O$ along with $P_{t_i}^S$ and $P_{t_i}^P$. The normal to $\mathsf{\Gamma }(P_{t_i}^S,P_{t_i}^P,P_{t_i}^O)$ at $P_{t_i,\beta }^O$ is denoted as ${\rm Z}(P_{t_i}^S,P_{t_i}^P,P_{t_i,\beta }^O)$. The vector from $P_{t_i}^O$ to $P_{t_i}^S$ , denoted as $V\left({P_{t_i}^O,P}_{t_i}^S\right)$, is perpendicular to the normal at $P_{t_i,\beta }^O$, ${\rm Z}(P_{t_i}^S,P_{t_i}^P,P_{t_i,\beta }^O)$.

To determine a position laterally offset from the central curve, denoted as $P_{t_i,\beta ,\phi }^O$, the offset length between the corresponding point on the central curve, ${\mathrm{P}}_{{\mathrm{t}}_{\mathrm{i}},\mathsf{\beta }}^{\mathrm{O}}$, and ${\mathrm{P}}_{{\mathrm{t}}_{\mathrm{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathrm{O}}$ can be determined using the geometrical relationships between the length of the vector from $P_{t_i}^O$ to $P_{t_i}^S$, denoted as $V\left({P_{t_i}^O,P}_{t_i}^S\right)$, which is perpendicular to the normal at ${\mathrm{P}}_{{\mathrm{t}}_{\mathrm{i}},\mathsf{\beta }}^{\mathrm{O}}$, $(P_{t_i}^S,P_{t_i}^P,P_{t_i,\beta }^O)$ and the angle, $\phi$, formed between the vectors $V(P_{t_i}^S,{\mathrm{P}}_{{\mathrm{t}}_{\mathrm{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathrm{O}})$ and $V(P_{t_i}^S,{\mathrm{P}}_{{\mathrm{t}}_{\mathrm{i}},\mathsf{\beta }}^{\mathrm{O}})$ as shown in Figure 2. The distance between $P_{t_i,\beta }^O$ and $P_{t_i,\beta ,\phi }^O$, denoted as $\lambda (P_{t_i,\beta }^O,P_{t_i,\beta ,\phi }^O)$, can be obtained by

$$\lambda \left(P_{t_i,\beta }^O,P_{t_i,\beta ,\phi }^O\right)=\mathrm{t}\mathrm{a}\mathrm{n}\left(\phi \right)\times \lambda (P_{t_i}^S,P_{t_i,\beta }^O)$$

(4)

The position of the offset point is thus:

$$P_{t_i,\beta ,\phi }^O=P_{t_i,\beta }^O+\lambda \left(P_{t_i,\beta }^O,P_{t_i,\beta ,\phi }^O\right)\times NV\left(P_{t_i,\beta }^O,P_{t_i,\beta ,\phi }^O\right)$$

(5)

where $\phi \in [-\alpha ,\alpha ]$, $NV\left(P_{t_i,\beta }^O,P_{t_i,\beta ,\phi }^O\right)$ is normalized $V\left(P_{t_i,\beta }^O,P_{t_i,\beta ,\phi }^O\right)$.

Using a ray-tracing algorithm we can determine the position and direction of scattered/reflected rays. However, these reflection points are generated by assuming a smooth surface. In practice, the surface is rough, and each illuminated surface element can be described as having a height profile with a particular root mean square height distribution, denoted here as $\mathrm{r}\mathrm{m}\mathrm{s}\mathrm{H}$, and a particular correlation length which determines the spatial frequencies over the surface in the x and y directions, denoted herein as $\mathrm{C}\mathrm{o}\mathrm{X}$ and $\mathrm{C}\mathrm{o}\mathrm{Y}$. In this model, we consider the local illuminated surface to be a square area with a local surface coordinate system $\mathrm{L}\mathrm{X}-\mathrm{L}\mathrm{Y}-\mathrm{L}\mathrm{Z}$. The origin of coordinate system is denoted as $\mathrm{L}\mathrm{O}$. We use the spacecraft as a reference point and maintain that the directions from north to south and from left to right are positive. The characteristics of the $\mathrm{L}\mathrm{X}-\mathrm{L}\mathrm{Y}-\mathrm{L}\mathrm{Z}$ system are shown in Figure 3.

Overlaid upon this coordinate system is an auxiliary square grid, with Gaussian elevation associated to the nodes of the grid, used to simulate the rough surface; thus, the number of auxiliary elements in the $\mathrm{L}\mathrm{X}$ or $\mathrm{L}\mathrm{Y}$ direction is $N+1$, resulting in ${(N+1)}^2$ total elements. The number $N+1$ must conform to the Nyquist theorem in the spatial domain such that the number of auxiliary elements per unit length is inversely related to the correlation length for its lower limit. This relationship also provides constraints for the side length $SL$ of the square facet (N/SL > 2/CoX).

The element side length visible to the transmitter at a certain time epoch is approximately:

$$SL=2\times \mathrm{t}\mathrm{a}\mathrm{n}\left(\alpha \right)\times \lambda (P_{t_i}^S,P_{t_i}^O)$$

(6)

To simplify the description, the generated auxiliary elements overlaid on the $\mathrm{L}\mathrm{X}-\mathrm{L}\mathrm{Y}-\mathrm{L}\mathrm{Z}$ coordinate system are numbered as ${\kappa }_{a,b}$ where $aϵ[1,N+1]$ and $bϵ[1,N+1]$ are row number and column number respectively. The elevation data of ${\kappa }_{a,b}$ is initialized by $H_{a,b}$.

Figure 4a shows an example of an $N\times N$ grid of auxiliary elements. To generate a random matrix with a Gaussian distribution, as seen in Figure 4b, we use MATLAB’s random number generator function $H=rmsH\times randn(N,N)$.

To produce a spatially correlated surface, this surface is convolved with a Gaussian filter using MATLAB’s discrete Fast Fourier Transform (FFT) function. When $\mathit{C}\mathit{o}\mathit{X}=\mathit{C}\mathit{o}\mathit{Y}$, the surface is considered isotropic, meaning the roughness is the same in the $\mathit{L}\mathit{X}$ and $\mathit{L}\mathit{Y}$ directions and thus scatters light evenly in all directions. To achieve a spatially correlated isotropic surface, the Gaussian filter with a low pass frequency response is

$$\mathit{G}\mathit{f}={\mathit{e}\mathit{x}\mathit{p}}^{-\left(\frac{\left|{\mathit{X}(\mathit{\kappa }}_{\mathit{a},\mathit{b}})\right|+\left|{\mathit{Y}(\mathit{\kappa }}_{\mathit{a},\mathit{b}})\right|}{\mathit{C}\mathit{o}\mathit{X}/2}\right)}$$

(7)

Where ${\mathit{X}(\mathit{\kappa }}_{\mathit{a},\mathit{b}})$ and ${\mathit{Y}(\mathit{\kappa }}_{\mathit{a},\mathit{b}})$ are coordinates of ${\mathit{\kappa }}_{\mathit{a},\mathit{b}}$ on $\mathit{L}\mathit{X}$ and $\mathit{L}\mathit{Y}$ respectively. Using this Gaussian filter in addition to the inverse $\mathit{F}\mathit{F}\mathit{T}\left(\mathit{I}\mathit{F}\mathit{F}\mathit{T}\right)$, to change from the spectral domain to the spatial one, we can obtain the rmsH for the facet,

$${\mathit{H}}_{\mathit{a},\mathit{b}}=2\times \mathit{X}\mathit{L}\times \mathit{I}\mathit{F}\mathit{F}\mathit{T}\left(\mathit{F}\mathit{F}\mathit{T}\left({\mathit{H}}_{\mathit{a},\mathit{b}}\right)\mathbf{*}\mathit{F}\mathit{F}\mathit{T}\right(\mathit{G}\mathit{f}\left)\right)/(\mathit{N}\times \mathit{C}\mathit{o}\mathit{X})$$

(8)

where * is the dot product between two vectors. The spatially correlated surface seen in Figure 4c, is the result of the convolution of H and the Gaussian filter for an isotropic surface. On the other hand, When $\mathit{C}\mathit{o}\mathit{X}$ and $\mathit{C}\mathit{o}\mathit{Y}$ are different, the surface is considered anisotropic, meaning the roughness varies in the $\mathit{L}\mathit{X}$ and $\mathit{L}\mathit{Y}$ directions thus scattering light in different directions with different intensities. To achieve a spatially correlated surface for an anisotropic, the Gaussian filter with a low pass frequency response is

$$\mathit{G}\mathit{f}={\mathit{e}\mathit{x}\mathit{p}}^{-\left(\frac{\left|{\mathit{X}(\mathit{\kappa }}_{\mathit{a},\mathit{b}})\right|}{\raisebox{1ex}{$\mathit{C}\mathit{o}\mathit{X}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+\frac{\left|{\mathit{Y}(\mathit{\kappa }}_{\mathit{a},\mathit{b}})\right|}{\raisebox{1ex}{$\mathit{C}\mathit{o}\mathit{Y}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)}$$

(9)

With this auxiliary grid we can correct the initial reflection point by adding the rms height ($\mathit{r}\mathit{m}\mathit{s}\mathit{H}$)

$${\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}={\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}\times \frac{\mathit{\lambda }({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{P}},{\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}})+{\mathit{H}}_{\mathit{a},\mathit{b}}}{\mathit{\lambda }({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{P}},{\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}})}$$

(10)

Now that the reflection points on the rough surface elements can be calculated there is the possibility that an overlapping or shadowed area possibly exists between two adjacent illuminated surface elements when the moving step of the spacecraft during the simulation is short. Considering that the elevation data of each illuminated surface is independently generated, it is necessary to correct the elevation data of the latter illuminated surface located in the overlapping part. Supposing the moving step is denoted as $\mathit{\varrho }$, the point belonging to the latter illuminated surface and overlapping part can be denoted as ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho },\mathit{\beta },\mathit{\phi }}^{\mathit{O}}$. After projecting the overlapping part on the $\mathit{L}\mathit{X}-\mathit{L}\mathit{Y}-\mathit{L}\mathit{Z}$ system, it is easy to find which element belonging to the former illuminated surface that ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho },\mathit{\beta },\mathit{\phi }}^{\mathit{O}}$ will be located in. Supposing this element is ${\mathit{\xi }}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho },\mathit{\beta },\mathit{\phi }}$, ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho },\mathit{\beta },\mathit{\phi }}^{\mathit{O}}$ can be corrected and replaced by the intersection of vector $\mathit{V}({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho }}^{\mathit{P}},{\mathbf{P}}_{{\mathbf{t}}_{\mathbf{i}}+\mathit{\varrho },\mathsf{\beta },\mathsf{\phi }}^{\mathbf{O}})$ and element ${\mathit{\xi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}$.

Once the reflection positions are determined, the direction vector of the reflected ray must be determined. The rough surface is meshed using a collection of non-overlapping triangular facets approximating the surface shape. In this way, the intersection of a ray with each triangular facet is calculated. The points selected to create the surface element are ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}$, ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }+\frac{2\times \mathit{\alpha }}{\mathit{N}}}^{\mathit{O}}$ and ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta }+\frac{2\times \mathit{\alpha }}{\mathit{N}},\mathit{\phi }+\frac{2\times \mathit{\alpha }}{\mathit{N}}}^{\mathit{O}}$, where ${\mathbf{P}}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathbf{O}}$ is selected as the reflection point of this surface element. An arbitrary facet denoted as ${\mathit{\xi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}$ is seen in Figure 5 projected above the local surface coordinate system.

To determine the normal to the triangular facet plane, ${\mathsf{\chi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}=\mathit{{\rm Z}}({\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta }+\frac{2\times \mathit{\alpha }}{\mathit{N}},\mathit{\phi }+\frac{2\times \mathit{\alpha }}{\mathit{N}}}^{\mathit{O}},{\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }+\frac{2\times \mathit{\alpha }}{\mathit{N}}}^{\mathit{O}},{\mathit{P}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}})$, the cross product of two edges is calculated. To find the vector of the reflected ray ${\mathit{R}}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathit{R}}$, the component of the incident ray is reversed parallel to the triangle facet normal such that

$${\mathit{R}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{R}}={\mathit{R}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{I}}-2\mathbf{*}\left({\mathit{R}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{I}}\mathbf{*}{\mathit{\chi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}\right)\mathbf{*}{\mathit{\chi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}$$

(11)

where * is the dot product, ${\mathit{R}}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathit{I}}$ is the incident ray.

3.2. Modeling the Reception of Reflected Signals

Compared with the reality, the simulation of radio waves in terms of quantity is usually under-estimated. Hence if the receiving condition for reflection signals is strictly set, just several reflected signals can be captured by the receiver and the analysis result will also be underestimated. As shown in Figure 6, a circular receiving range $\mathit{\Psi }$ is utilized to define the receiving ability of receiver. This circle is perpendicular to the vector from ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{O}}$ to ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{R}}$, denoted as $\mathit{V}({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{O}},{\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{R}})$, and the center of this circle is ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{R}}$. However, as has mentioned before, the underestimated quantity of simulated signals will cause negative influence on the analysis result. So, this receiving plane should be put at a nearer position to the explored planet, which can be denoted as ${\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{R}\mathit{\text{'}}}$. The intersection point of ${\mathit{R}}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathit{R}}$ and receiving plane is denoted as ${\mathit{\delta }}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}$. The reflected signal is regarded to be received when meeting that

$$\mathsf{\lambda }({\mathit{\delta }}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }},{\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{R}})\le \mathit{\Psi }$$

(12)

3.3. The Doppler Shift of the Received Signals

Figure 7 shows the relative spatial positions of transmitter, surface and receiver, as well as several incident and reflected signals, where the sum of ${\mathit{L}}_{1\mathit{x}}$ and ${\mathit{L}}_{2\mathit{x}}$ are the total path length of the ray. Assume the transmitting signal is denoted as

$${\mathit{R}}^{\mathit{I}}=\mathit{A}\times \mathit{c}\mathit{o}\mathit{s}(2\times \mathit{\pi }\times {\mathit{f}}_{\mathit{T}}\times \mathit{t}+{\mathit{\varphi }}_0)$$

(13)

Assuming that each reflected signal reaching the receiver have the same field intensity. Hence, the reflected signal reaching receiver can be denoted as

$${\mathit{R}}^{\mathit{R}}={\mathit{A}}^{\prime }\times \mathit{c}\mathit{o}\mathit{s}(2\times \mathit{\pi }\times {\mathit{f}}_{\mathit{T}}\times (\mathit{t}+\mathit{\Delta }\mathit{t})+{\mathit{\varphi }}_0)$$

(14)

where ${\mathit{f}}_{\mathit{T}}$ is the frequency of transmitting signal and $\mathit{\Delta }\mathit{t}=\frac{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}{\mathit{v}}$. $\mathit{v}$ is the velocity of radio wave, with the value of 3 ∗ ${10}^8$m/s. Equation 19 can be further expressed by

$${\mathit{R}}^{\mathit{R}}={\mathit{A}}^{\prime }\times \mathit{c}\mathit{o}\mathit{s}(2\times \mathit{\pi }\times {\mathit{f}}_{\mathit{T}}\times (\mathit{t}+\frac{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}{\mathit{v}})+{\mathit{\varphi }}_0)$$

(15)

The equation of signal via element ${\mathit{\xi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}$ is denoted as

$${\mathit{R}}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{R}}={\mathit{A}}^{\prime }\times \mathbf{cos}[2\times \mathit{\pi }\times {\mathit{f}}_{\mathit{T}}\times (\mathit{t}+\frac{\mathit{\lambda }({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{S}},{\mathbf{P}}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathbf{O}})+\mathit{\lambda }({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}}^{\mathit{R}},{\mathbf{P}}_{{\mathbf{t}}_{\mathbf{i}},\mathsf{\beta },\mathsf{\phi }}^{\mathbf{O}})}{\mathit{v}})+{\mathit{\varphi }}_0]$$

(16)

Therefore, the frequency of received signal can be expressed by

$$\begin{gathered} {\mathit{f}}_{\mathit{R}}=\frac{1}{2\times \mathit{\pi }}\times \frac{\mathit{d}\left[2\times \mathit{\pi }\times \mathit{f}\times \left(\mathit{t}+\frac{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}{\mathit{v}}\right)+{\mathit{\varphi }}_0\right]}{\mathit{d}\mathit{t}} \\ =\frac{1}{2\times \mathit{\pi }}\times \left[\frac{\mathit{d}\left(2\times \mathit{\pi }\times \mathit{f}\times \mathit{t}\right)}{\mathit{d}\mathit{t}}+\frac{\left(2\times \mathit{\pi }\times \mathit{f}\times \frac{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}{\mathit{v}}\right)}{\mathit{d}\mathit{t}}+\frac{{\mathit{\varphi }}_0}{\mathit{d}\mathit{t}}\right] \\ =\frac{1}{2\times \mathit{\pi }}\times \left[2\times \mathit{\pi }\times \mathit{f}+\frac{2\times \mathit{\pi }\times \mathit{f}}{\mathit{v}}\times \left(\dot{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}\right)\right] \\ =\mathit{f}+\mathit{f}\times \frac{\dot{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}}{\mathit{v}} \\ =\mathit{f}\times (1+\frac{\dot{{\mathit{L}}_{1\mathit{x}}+{\mathit{L}}_{2\mathit{x}}}}{\mathit{v}}) \end{gathered}$$

(17)

In the frequency domain, the received overall signal will show a Gaussian shape because of the Doppler shift. To get the frequency of overall received signals, a simple method is fast Fourier transform (FFT). However, because usually the signal frequency ${\mathit{f}}_{\mathit{T}}$ configured for spacecraft is high, and the sampling frequency conducting FFT should be higher than ${\mathit{f}}_{\mathit{T}}$, the performance of FFT will cause unacceptable computation consumption. Considering this drawback, the method employed here is to record the frequency of all the received signals and compute the largest frequency difference. This frequency difference can be approximated to be $3\mathit{\sigma }$ confidence interval of the Guassian shape and is crucial to obtain the half-power width of the frequency-domain signal.

3.4. The Computation of the Rms Slope

For radar signature analysis, the rms slope ${\mathit{\tau }}_0$ rather than the rms height, $\mathit{r}\mathit{m}\mathit{s}\mathit{H}$, is typically used. The rms slope is the derivative of the local surface tilt in the z direction as a function of the length of the local facet. The computation of ${\mathit{\tau }}_0$ is conducted on the statistics of all surface elements. As shown in Figure 8, the slope for element ${\mathit{\xi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}$, denoted as ${\mathit{\tau }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}$, can be obtained by

$${\mathit{\tau }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}=\mathbf{a}\mathbf{c}\mathbf{o}\mathbf{s}(\frac{{\mathsf{\chi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}\mathbf{*}\mathit{V}({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho }}^{\mathit{P}},{\mathbf{P}}_{{\mathbf{t}}_{\mathbf{i}}+\mathit{\varrho },\mathsf{\beta },\mathsf{\phi }}^{\mathbf{O}})}{{|\mathsf{\chi }}_{{\mathit{t}}_{\mathit{i}},\mathit{\beta },\mathit{\phi }}^{\mathit{O}}\left|\mathbf{*}\right|\mathit{V}({\mathit{P}}_{{\mathit{t}}_{\mathit{i}}+\mathit{\varrho }}^{\mathit{P}},{\mathbf{P}}_{{\mathbf{t}}_{\mathbf{i}}+\mathit{\varrho },\mathsf{\beta },\mathsf{\phi }}^{\mathbf{O}})|})$$

(18)

The rms slope ${\mathit{\tau }}_0$ of the illuminated surface is then obtained by

$${\mathit{\tau }}_0=\mathbf{r}\mathbf{m}\mathbf{s}\left(\left[\begin{array}{ccc}{\mathit{\tau }}_{{\mathit{t}}_{\mathit{i}},-\mathit{\alpha },-\mathit{\alpha }}& \cdots & {\mathit{\tau }}_{{\mathit{t}}_{\mathit{i}},-\mathit{\alpha },\mathit{\alpha }}\\ ⋮& \cdots & ⋮\\ {\mathit{\tau }}_{{\mathit{t}}_{\mathit{i}},\mathit{\alpha },-\mathit{\alpha }}& \cdots & {\mathit{\tau }}_{{\mathit{t}}_{\mathit{i}},\mathit{\alpha },\mathit{\alpha }}\end{array}\right]\right)$$

(19)

4. Numerical Simulation Experiments

To investigate the reliability and numerical stability of BIRDS, in this section we study the influence of the number of incidence waves on the backscattering surface of the planetary body. This is useful to quantify the computational burden of each simulation as a function of the achievable accuracy in the estimation of the surface roughness.

4.1. Case 1: Mars Bistatic Radar Experiments

We start by simulating a BSR experiment scenario at Mars, remembering that a usual problem encountered in the real Mars environment is that, with a typical S/C and orbital design, the beam is limited, and cannot work efficiently on Mars’ high-rough surface. To prevent this problem, in this simulated Mars BSR experiment, the half-beam width of the onboard HGA is fixed to ${15}^{\mathbf{o}}$. Moreover, to maximize the number of amount of received signals reflected from the surface, we fixed a receiving mask cone of 10° aperture, which is an unrealistically large value, but this avoids simulating a huge number of RF beams, so to speed up the numerical simulation. We decided to consider, as the only independent parameter, the angular separation between two adjacent incidence waves, which can be denoted as $∆\mathbf{A}$. In our simulation, we varied $∆\mathbf{A}$ between 0.155° to 0.075° (with steps of 0.02°), which means the total number of incident waves, simulated in the experiment, varied between 29387 and 160000. For each numerical experiment carried out using a different value of $∆\mathbf{A}$ we computed the half-power bandwidth of the total simulated received RF signal (which is the sum of all the individual simulated waves) to compare it to the analytical prediction given by Eq. (3). The left panels of Figure 8 depict this comparison for a variety of different BIRDS-simulated RMS slopes of the illuminated area on the surface of Mars. From the top row to the bottom row the $∆\mathbf{A}$ angle is gradually decreased and, while the prediction of Eq. (3) is shown as a solid line, the results obtained from 5 independent runs are given by box plot. The right panels of Figure 8 also show the number of the received waves in the five different values of the $∆\mathbf{A}$ angle. For each RMS slope, we show our results as follows: the red horizontal line represents the average value of the received signals (for the same five runs of BIRDS described above), the blue box represents the standard deviation of the same fine runs, and the black dashed lines represent the minimum and maximum limits of the five runs. These plots show that while increasing the number of the received signals (via decreasing the $∆\mathbf{A}$ angle) the half-power bandwidth of the received echo from BIRDS gradually approaches the result of Eq. (3) and the vertical bars become smaller, see left panels of Figure 8. In the current range of RMS slope for the illuminated surface, the number of received waves increases as a quadratic function, see right panels of Figure 8, and the panel in the last row (plot (j)) offers the best match to a quadratic function, with the lowest variabilities.

The left panel of Figure 9 depicts the average value of the error between the numerical simulations and the analytical prediction given by Eq. (3) by generating different number of incident signals. It shows that when the number of received waves is larger than 40000, the error can be limited to about 50 Hz. When the number of incident waves is 160000, the error can be limited to about 25 Hz. It indicates that increasing the quantity of incident waves can effectively improve the system’s stability. The standard deviation of the error between the numerical simulations and the analytical prediction given by Eq. (3), for 5 runs of the BIRDS simulator, is also shown in the right panel of Figure 9. This shows that the standard deviation of the error in the half-power bandwidth for a small quantity of incident signals (29387, in particular, in Figure 9) is totally unreliable, as it is remarkably small (even smaller than the average error), an artifact of the few received signals. On the other hand, when the quantity of the received waves reaches 90000, the standard deviation of the error can be limited to about 35 Hz.

In Table 1, the rms slope of the generated surface and the one computed from Eq. (3) are compared. The rms slope computed using Eq. (3) relies on the rms height and irregularity horizontal length. The rms heights (in km) are given in the first column, and the rms slopes computed using the numerical simulations carried out under the different hypothesis of grid points are shown in the second to sixth columns, while the rms slopes computed via Eq. (3) are given in the last column. An interesting phenomenon can be observed is that by referring to the results shown in Figure 8 and in Table 1, the computed bandwidth is relevant to the precision of the reconstructed rms slope. In order to lower the error in the rms slope estimation, the computed bandwidth must closely reach the result of Eq (3).

4.2. Case 2: JUICE

The performance of BIRDS was also assessed using as a reference scenario a trajectory designed for ESA’s JUICE mission [36] during the so-called Ganymede orbital phase (Jupiter Icy moons Explorer Consolidated Report on Mission Analysis v.5, ESA Document [37]). Here, we explore the capability of BIRDS to cope with a dynamic environment where the S/C is moving and flying over a terrain with varying values of surface roughness.

We set the scene in the time range from 2032-09-01-00:00:00 to 2032-09-01-00:04:00 with a period of four minutes. The simulated S/C radio frequency is at X-band (around 8.4 GHz), the beamwidth is 2 degrees, the irregularity factor for generating the simulated planetary surface is 0.25 km, the angular separation between two adjacent incidence waves is 0.1 degrees, and the same reception mask cone of 10° described in Section 4.1 was used. We selected eight epochs, from 2032-09-01-00:00:03 to 2032-09-01-00:00:17, with an increasing step of 2 seconds. For each illumination site, the rms heights are randomly generated as 0.003 km, 0.002 km, 0.002 km, 0.003 km, 0.0025 km, 0.0025 km, 0.003 km and 0.004 km, respectively. The rms slopes of these eight illumination surface are ${1.1679}^{\mathrm{o}}$, ${0.7786}^{\mathrm{o}}$, ${0.7786}^{\mathrm{o}}$, ${1.1679}^{\mathrm{o}}$, ${0.97325}^{\mathrm{o}}$, ${0.97325}^{\mathrm{o}}$, ${1.1679}^{\mathrm{o}}$ and ${1.5572}^{\mathrm{o}}$, respectively. In each run of BIRDS, the surface is randomly generated, so the obtained bandwidth is slightly different. In order to avoid numerical instability problems, each case is run independently 5 times, and the bandwidths of the received waves are presented in box plot in Figure 10. These results also show good fitting to the results obtained from equation (3).

5. Conclusions

We presented a 3D model to characterize the scattered field of a rough surface on an extra-terrestrial body, aimed at simulating the expected performance of future orbital bistatic radar experiments. In particular, we have investigated how the variability of surface roughness impacts the spectral broadening of the received signal using physical optics approximations and ray tracing on a surface model using a facet approach with Gaussian properties. We have assessed the limits and the effectiveness of numerical simulations, focusing on the minimum number of rays impacting the planetary surface and then reaching the receiver, in order to obtain statistically meaningful results. The predicted analytical linear relationship between the spectral broadening of the signal and the surface roughness was confirmed in our numerical simulation tool and was applied to two scientifically compelling scenarios of practical interest at Mars and at Jupiter’s moon Ganymede. For the latter case, the simulator could be used to carry out a preliminary assessment of JUICE BSR in the future; using realistic values of half-power beamwidth, more facets, and more computational burden, we can better understand the occurrence of beam-limitation and explain the shape of real echoes collected using JUICE 3GM X or Ka band signals.