## 1. Introduction

Engineers using microwave radio links have to address the effects of multipath propagation that arise when several rays arrive at the receiver after travelling along different paths from the transmitter. Rays, along their propagation, undergo reflections at the earth’s surface or at variations in the refractive index or its gradient.

Since the ‘60s, standard ray theories for radio wave applications use effective earth radius concepts, ray bending based on Snell’s law in a layered spherical atmosphere, or analytical method, limited to simple refractive index profiles (Du Castel, 1966; Livingston, 1970). At that time, scientists were interested in statistics of fading, distinguishing between *fast* or *slow fading* and the typical results were either nomograms of attenuation as a function of distance for a set of common frequency and height values or cumulative distribution functions representing the fraction of time that the signal was expected to be received at or below a given level (CCIR, 1978). Several techniques were employed to calculate the field strength, such as Cornu's spiral derived by the Huygens' principle together with the introduction of the specular and the diffuse reflection contribution (Hall, 1979). In that approach, the field amplitude was determined assuming the principle of conservation of energy, i.e. the flux of the energy along a ray is the same at every cross section of a narrow tube of rays (Keller, 1957).

Nowadays, radio communications live a renaissance due to the development of mobile and wireless communications and the increasing demand of wideband services. The need of greater rate of transmission is accomplished using higher frequencies for which effects like the tropospheric scatter can affect dramatically the quality of the received signal. Thus, aspects such as analysis of delay power spectra or of potential intersystem interferences require a more realistic modelling. As a matter of fact, local fluctuations in the refractive index *n* can cause scatter, while its abrupt changes with the height can cause reflection and originate ducting layers, in which the rays are reflected and refracted back in such a way that the field is trapped inside these layers (Bean & Dutton, 1968).

Although the goal of radio engineering is to transmit reliably from the source to the receiver, this is of no use if the received signal is unintelligible due to interference or if the transmission causes interference with other systems. So an analysis of the channel behaviour is mandatory to meet both link design and system requirements. Besides, the availability of geographic databases, digital elevation models and the increasing computing power let us face more realistic description to model the interaction of the propagating wave with the atmosphere and the surrounding scenario (Driessen, 2000; Kurner & Cichon, 1993; Lebherz et al., 1992).

In more recent developments, ray modelling of wave propagation addresses the dispersive effects of perturbed atmosphere on the performance of high-capacity digital radio channel (Akbarpour & Webster, 2005; Sevgi & Felsen, 1998).

Besides telecommunications applications, a proper modelling of atmospheric propagation is of concern in atmospheric sciences with particular focus on radio occultation data analysis for assimilation with numerical weather prediction models (Pany et al., 2001). Also, the decrease of propagation speed related to the atmospheric refractivity causes tropospheric delay, which influences applications of the Global Navigation Satellite System (Eresmaa et al., 2008).

The author presents a widely diffuse technique in the domain of seismic studies that accommodates lateral variations in the medium properties (Farra, 1993). As an asymptotic technique, based on high frequency approximation, it permits fast computation but provides a local solution of propagation problem (Červenŷ et al., 1977). The technique was adapted for application with electromagnetic waves and specifically tailored for signals travelling in the atmosphere. The hypothesis of horizontal uniformity can be removed and no stratification is needed to calculate the ray trajectory. The terrain profile coordinates are mapped over a Cartesian reference through analytical transformation. Refractive index values are given in the same Cartesian reference. Hence, propagation is modelled in a two dimensional range-height scenario over irregular terrain through non-homogeneous atmosphere.

The field amplitude is evaluated by means of a perturbation technique using paraxial rays to observe the wave front structure along the path from the transmitter towards the receiver.

This approach allows to analyze system performance and the channel impulse response in presence of any kind of atmosphere, characterized by local values of refractive index *n* in the bi-dimensional panel including the antennae and the terrain profile. The variations of the refractive index *n* along the third dimension are taken equal to zero; nonetheless, the field amplitude is calculated correctly because of the paraxial approximation which takes place in a 3D domain. Median power strength of the numerical results was compared with predictions given by Friis' Formula (Balanis, 1996).

## 2. Modelling

The proposed ray tracing technique is widely used in seismic for subsoil investigation and it follows an approach based upon the integration of the Eikonal equation with Hamiltonian-Jacobi technique (Kružkov, 1975). According to this formulation, rays are defined by the vector *y(τ)=(x(τ), p(τ))*, where *x(τ)* and *p(τ)* are the position vector and the slowness (inverse of phase velocity) vector along the ray, both function of the integration variable *τ* and of the initial conditions (launching point and direction). The slowness vector *p* is defined as *y*(*τ*) satisfies the Hamilton differential equations, whose solutions describe the wave propagation in the medium, under the asymptotic assumption (Abramovitz & Stegun, 1970):

where _{P} and x represent the gradient computed versus the slowness vector and the position vector, respectively. H is the Hamiltonian function describing the wave propagation in the considered medium, i.e. the atmosphere, chosen as follows:

being v (x) the propagation speed of the medium at the position vector x.

The advantage of such an approach is twofold: it takes into account the medium inhomogeneities that originate multipath propagation and wavefronts folding (caustics), and it permits to consider both vertical and long range variations of the atmospheric model, without any approximation like flat-earth model (Hall, 1979). Also, it permits to keep separated the different field contributions due to the different interactions between the rays and the surrounding scenario, which impacts on the phase calculation of the overall received field.

### 2.1. Ray tracing

A suitable discussion for ray tracing can be found in (Farra, 1993); nevertheless, in this section some basic concepts and analytical details are reported.

Substituting (2) in (1), one obtains the equations for the vector y(τ):

x and p define the 6D phase-space, where the solution of (3), i.e. the rays, are the characteristic lines of the Eikonal equation, thus enabling the interpolation without ambiguity on a single fold (Vinje et al., 1993). For this reason rays are also called bi-characteristics lines of the wave equation. By choosing as integration variable the ray propagation delay or travel time T, system (3) becomes:

The solution of system (4) describes, under the asymptotic assumption, the propagation in the atmosphere in terms of ray trajectories and time delay.

### 2.2. Amplitude calculation

Anomalous conditions such as medium inhomogeneity or velocity model variations can affect dramatically the wavefront propagation. Under these circumstances the computation based upon the spherical divergence assumption results in erroneous evaluations.

One possible choice is to compute the amplitude through paraxial rays (first order perturbation theory), used to determine the flow tube, whose deformations, together with the theorem of the conservation of the flux along the ray path, allow the calculation of the amplitude A at the time τ.

Let us consider a reference ray with characteristic vector y_{0}(τ)=(x_{0}(τ), p_{0}(τ)). A paraxial ray is obtained from the reference one by applying the first order perturbation theory, so paraxial rays coordinates are defined by:

Tracing paraxial rays consists in finding in the phase-space the canonical perturbation vector δy(τ)=(δx(τ),δp(τ)). These perturbations in the trajectory are due to small changes in the initial conditions of the ray.

The linear system for the calculation of paraxial rays is obtained by inserting the perturbation vector given by (5) in system (1) and developing to the first order:

where

is a 6 6 matrix whose elements are the derivatives of the Hamiltonian function calculated on the reference ray.

Paraxial rays allow to define the flow tube, schematically represented in figure 1.

The deformation of the flow tube along the ray path

where τ is the sampling parameter, θ the elevation angle and ϕ the azimuth angle. The amplitude of the plane wave associated to the ray is computed considering both the deformation of the flow tube and the conservation of power density flow law along the ray path. The volume element of the flow tube can be expressed as

From the conservation of the power density flow along the ray (Keller, 1957) and known the amplitude A_{i} at τ_{i}, it follows that the amplitude A_{i+1} at τ_{i+1} is:

Equation (10) shows that the ray amplitude depends on the velocity model. Hence, paraxial rays take into account anomalous propagation conditions.

The singularity of the Jacobian J is a singularity, known as caustic, of the plane wave associated to the ray (Kravtsov & Orlov, 1990). Caustics arise when the ray field folds. Caustics are of first order when the flow tube looses one dimension and of the second order when two are the dimensions lost. The wave front folding affects the phase of the field carried by the wave with a phase shift of

### 2.3. Parameterization

The modelling is performed in a 2D panel, as far as the terrain profile and the atmosphere characteristics are concerned, while the ray trajectories are computed in a 3D space.

This assumption is based on the hypothesis that lateral variations in the propagation medium and in the terrain profile are negligible in the third dimension, which does not strictly holds.

Therefore, the propagation occurs in the so called Earth system (s,h), where s is the range measured along the idealized spherical earth and h is the altitude taken along radial direction passing through the Earth centre, respectively. Instead, for sake of simplicity, all calculation are developed in a Cartesian reference system (x,z) related to the Earth by the following transformation relationships:

where R_{0} is the Earth radius. The transformation from cartographic coordinates (s,h) to absolute coordinates (x,z) is given by the inverse transformation:

being

This way to proceed allows to take into account the actual geometry of the problem without introducing approximations like equivalent Earth radius. Also, it keeps separated the domain in which the radio link characteristics are defined and the computational one, where rays are traced.

The characteristics of the atmosphere are taken into account in terms of its propagation velocity using several quantities that are generally function of the position vector x:

where c is the speed of the light and n(x) is the refractive index of the atmosphere, which usually is expressed in terms of the refractivity N:

N is a dimensionless quantity but in literature it is often measured in N-units or in parts per million (ppm). Its value depends upon the altitude and the range locally according to the atmospheric pressure P, the vapour pressure P_{V} and the temperature T as in the following (Bean and Dutton, 1968) :

Besides the variations in the refractive index, or its gradient, rays along their propagation undergo reflections and diffraction at the earth’s surface. At the present stage of the modelling, neither the diffraction mechanisms nor any scattering features are included except reflection.

The trajectory of the reflected ray is given by the Snell’s law under the hypothesis that locally the ground acts as a perfect smooth surface. The initial conditions with which the reflected ray is traced depend both on the geometric characteristics of the terrain profile and on the ray direction of incidence. The amplitude and the phase of the rays reflected from the ground are computed according to the chosen ground parameterization, i.e. the electric soil properties. These determine the value of the ground impedance η, which is a function of both soil permittivity ε_{r} and conductivity σ (Ulaby, 1999):

being ε_{0} and ω the free space permittivity and the angular frequency respectively.

Once that the rays trajectories are computed, together with their complex amplitudes - weighted by the antenna pattern - and the travel times, they are classified into different wave fronts and then interpolated on the locations where the field is desired. The total field is obtained by addition of the single contribution of each wave front.

## 3. Received field estimation

The field reconstruction is described by means of a pilot example of a 4.5 GHz point-to-point radio link 80 km long affected by multipath. The aim of the technique is to evaluate the vertical electric field (or the received power) intensity at the receiver range. According to the considered geometry, this leads to the prediction of the vertical profile of the electric field E(h) (or P(h)), where h is the receiver height ranging between

The modelled atmospheric condition is indicated in the panel of figure 4, while figure 5 shows the vertical profile of N at few chosen distances. The N value is comprised between 280 and 360 N-units and it shows variations in both dimensions. This atmospheric model -taken from (Bean and Dutton, 1968), p.325 - let us focus on several kind of interactions and their effects on the propagating signal that could possible occur in actual conditions.

To appreciate the role of the variations in the refractivity, and without loss of generality, the ground profile was chosen as a plane surface characterized by the values of permittivity ε_{r} = 3 and conductivity σ = 0.001 [S/m] in the frequency band of interest.

The top image in figure 6 shows the trajectories of the rays belonging to a chosen elevation angular sector (-0.5 ÷ 0.5 with respect to the local horizon) traced in the cartographic system. The transmitter is in the origin of abscissas axis at 100 m above the ground. Besides the wavefront that carries the direct arrival of the signal, one may observe that several wavefronts are involved, resulting from the inhomogeneities of the atmosphere and ground reflection.

The bottom image of figure 6 shows the rays pertaining to the 30 wavefronts that results from the ray tracing process. The first step consists in grouping rays in different wave fronts, each characterized by the same history such as same number and location of caustics or reflections. Wavefronts are numbered and rays belonging to the same wavefront can be used to interpolate the field. As a matter of fact, rays arrive at the receiver range at discrete values so that the vertical profiling of the field intensity involves three steps: wavefronts separation; interpolation of the ray comb belonging to the same wavefront; coherent addition of the wavefronts themselves.

### 3.1. Wave front separation

The first criterion of wavefront separation that could come in mind is that of the travel time arrivals. Figure 7a shows the multi-valued behaviour of the wavefront delay of arrivals at different heights, while figure 7b shows how the parameterization of the wavefronts in the ray launching angle resolves the ambiguity: working in the angle domain seems the natural way for unfolding multi-valued ray fields (Operto et al., 2000).

The wavefronts separation is based on a broad and on a narrow selection process applied consecutively. The broad selection process involves the ray history: as for example, two rays belong to two different wavefronts if they are reflected a different number of times. The narrow selection process uses as criterion the angular distance: if two rays are reflected the same number of times they belong to the different wavefronts when their elevation starting angles difference is greater than δθ. This angular distance has the role of guarantee the proper accuracy in the interpolation of the ray field and it depends on the velocity model of the medium (Sun, 1992).

The results of the wavefront separation process are organized in a database whose structure is schematically represented in figure 8. Wavefronts and rays are organized in a record frame that gives for the k^{th} wavefront the propagation delay τ_{k}, the amplitude A_{k} and the phase ϕ_{k} obtained by interpolation of the travel time T_{ik}, a_{ik} and φ_{ik} of the rays at their arrival height, where i pertains to the individual rays in the wavefront k.

Each ray record refers is composed of the following parameters: the height of arrival, the travel time, which depend on the ray history; the amplitude and the phase, which depend on the ray history, on the reflection coefficients and on the system parameters.

The various quantities calculated by the ray tracing are referred to as Green’s function attributes. Figure 9 shows among the attributes those that directly contribute to the field calculation.

### 3.2. Wave front interpolation

The field calculation along the vertical profile at the receiver range is performed through the interpolation of the different wavefronts contribution. Let K be the number of the wavefronts, each including N_{k} arrivals, the electric field associated to the i^{th} arrival of the k^{th} wavefront is:

where f is the frequency of operation of the link. The contribution of the k^{th} wavefront to the total field is obtained through interpolation of the travel times, of the amplitudes and of the phases in the following way:

where g_{k} is a generic interpolating function and h is the height. As result, the interpolation of local values returns the field intensity in the target zone. Finally, the overall field is computed adding coherently all the wavefronts contributions:

The vertical profile of the received power is obtained by applying locally the Poynting's theorem to a spherical surface with radius R:

being η the medium impedance and R the receiver range.

The solid line in the leftmost panel of figure 10 shows the power profile due to the direct arrival, while the dotted line represents the theoretical values predicted by Friis’ formula (Balanis, 1996). Here, the vertical power profile is shown along the entire computational domain, from 0 up to 300 m, in order to show the effects of the tropospheric multipath on the received signal. The panel in the middle and the rightmost one of figure 10 show the power profile computed adding all the wavefront contribution at vertical (V) and horizontal (H) polarization.

## 4. Radio channel analysis

### 4.1. Time domain approach

The aim of time domain approach is to estimate the radio channel response to a transmitted pulse with a given waveform

where H represents the Hilbert transform operator. The corresponding pass band signal, which is transmitted on radio channel, is obtained with a frequency shift:

The channel pulse response can be written as a sequence of ideal shifted pulses of the transmitted waveform with their own amplitude and phase:

where the amplitudes, the phases and the propagation delays have been calculated, respectively, through interpolation of a_{ik}, φ_{ik} and T_{ik} at height h. These procedure has low computational cost as the ray trajectories and propagation delay τ_{k} do not depend on the frequency and they can be computed once for all. In equation (23) some quantities such as amplitudes A_{k} and phases ϕ_{k} depends on the frequency of operation but using a narrow band approximation one can assume that amplitude and phase variations are negligible within the signal bandwidth. Under these hypothesis, the pass band received signal

The received base band signal

whose real part gives the received signal

The panel in figure 11 reports the estimated field and the time delay as a function of height ranging between _{0}= 125 m and Δh = 25 m observed at the receiver range. This represents the impulse response of the system computed at vertical polarization. The input signal is a Nyquist wavelet and it is drawn in the left image of figure 11 as a dotted line. Also, figure 11 shows the received signal (solid line) at four given height of the panel. Two in-phase but slightly delayed arrivals are evident in panel a). In panel b) there are 2 arrivals, with one of them with a π/4 phase shift, while in panel d) the arrivals are of opposite phase.

### 4.2. Frequency domain approach

The computation of the electric field in the target zone for a given frequency f_{0} as reported in §3.2 strictly holds under the hypothesis of the transmission of a sinusoidal signal, which has an infinite duration and an infinitesimal band. As radio systems transmit signals with finite duration and band, there is the need to evaluate the response of the radio channel to the effects of frequency variation inside the system bandwidth of operation.

Let us assume that the radio link is transmitting a signal with bandwidth B around the carrier frequency f_{0}. Therefore, given the frequency f and using the equation (17) principle, the received electric field expression is:

Applying equation (27) for each frequency in the interval spanning from _{k} and the phase terms ϕ_{k}, which both depends (directly or indirectly) on the frequency of operation f.

Frequency domain approach can lead also to a characterization of radio channel in the time domain applying Fourier transform to the radio channel transfer function provided by equation (27). However, when the bandwidth of the transmitted signal is narrow, that is when f_{0} >> B, this procedure leads to a waste of computational resources. A more efficient way to manage narrow band channels is to face the problem in the time domain.

The transfer function of the radio channel is shown in figure 12 together with the frequency behaviour computed at the same four height levels highlighted in the time domain analysis. Again the target zone is comprised between 100 and 150 m and the signal has vertical polarization. One can notices how the transfer function may change from a smooth behaviour, panel c) to a sharper one with a notch deepness of about 30 dB, panel a).

## 5. Numerical validation

The ray tracing technique was developed to set up a prediction tool that is operative both as a MATLAB-based package and a C++ executable code. Unfortunately, beside an industrial application whose results are not available for publication at this stage, the computed field and power vertical profiles were not validated against measurements. Nevertheless the quality of the prediction performance can be observed when compared against canonical scenarios or other authors reference solutions.

Usually, when developing solver codes, one starts with the comparison of the predicted results with analytical or assessed solutions to check the correctness of the implementation. If this is the proper choice at that stage, it could lead to a flat and formal exposition of the methodology under study when describing the single steps towards the overall result. That is why the particular example chosen to illustrate the technique throughout §3 and §4 is far complicated, as it appears from the number and variety of wavefronts and from the aspects of the radio link characteristics both in time and in frequency domain. As the complexity could conceal errors or misevaluations, the author chose as reference solution an example in which a standard atmosphere model describes the velocity of the propagation medium over a flat perfect electric conductor ground. Even in such a simple example, the importance of the wave front discrimination could induce in errors when reconstructing the field.

### 5.1. Two ray example

The transmitter, in the origin of the abscissa axis, is at 150 m above the ground and the receiver range is 50 km apart. The atmospheric model is a standard exponential one with a gradient of 40 ppm in the first 1000m. Figure 13 shows the trajectories together with the two wavefront, one associated with the direct arrival and another with the reflected one. The ground is a perfect conductor with reflection coefficient equal to -1.

The two arrivals, direct and reflected, are shown in figure 14: on the leftmost panel the wavefront time delay dispersion results in 2 distinct contribution at the same instant. As in §3.1, the parameterization in the angular domain nicely unfolds the integration path (rightmost panel).

The predicted vertical received power profile is reported in figure 15, together with the direct arrival contributions given by Friis’ formula. One can notice the typical pattern of zero and maxima of reception due to the rotating in-phase and counter-phase summation of the two arrivals.

Finally, figure 16 shows the impulse response panel calculated for the two ray model and vertical polarization. Again, the received signal shows different behaviour due to the atmospheric multipath as a function of the height.

## 6. Conclusion

The ray tracing technique based on the resolution of the Eikonal equation with Hamilton-Jacobi technique is well suited to model the dispersive effects of perturbed atmosphere on the performance of high-capacity digital radio channel.

The advantage of such an approach is the possibility of describing the medium with its variation in the lateral domain in presence of atmospheric multipath and of non uniform terrain profile.

The field profiling at the receiving site is obtained through 2,5D propagation as, notwithstanding all variations in the third dimension are taken equal to zero, the field amplitude is calculated correctly because of the paraxial approximation which takes place in a 3D domain.

At the time of writing the author is refining the code towards a full 3D implementation, which requires a reformulation in spherical coordinates, and the modelling of the mechanism of ground scattering and diffraction.

## Acknowledgments

The author would like to thank Professor Giuseppe Drufuca (Politecnico di Milano) for his useful comments and suggestions and Simone Re (former research associate, Politecnico di Milano) for his contribution to the code development and calculations.