Propagation of Radio Waves in a Realistic Environment using a Parabolic Equation Approach

(1)

Linköping University | Department of Physics, Chemistry and Biology Master’s Thesis, 30 hp | Applied Physics and Electrical Engineering Spring term 2019 | LITH-IFM-A-EX--19/3636--SE

Propaga on of Radio Waves in a

Realis c Environment using a

Parabolic Equa on Approach

Jonas Ehn

Examiner: Peter Münger Supervisor: Tobias Hansson

(2)

(3)

Avdelning, Institution Division, Department

Division of Theoretical Physics

Department of Physics, Chemistry and Biology Linköping University

SE-581 83 Linköping, Sweden

Datum Date 2019-06-10 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport

URL för elektronisk version —

ISBN — ISRN

LITH-IFM-A-EX--19/3636--SE Serietitel och serienummer Title of series, numbering

ISSN —

Titel Title

Utbredning av radiovågor i en realistisk miljö genom paraboliska ekvationer

Propagation of Radio Waves in a Realistic Environment using a Parabolic Equation Approach

Författare Author

Jonas Ehn

Sammanfattning Abstract

Radars are used for range estimation of distant objects. They operate on the principle of sending electromagnetic pulses that are reflected off a target. This leads to the propagation of electromagnetic waves over large distances. As the waves propagate, they are affected by several aspects that decrease the performance of the radar system. In this master thesis, we derive a mathematical model that describes electromagnetic propagation in the tropo-sphere. The model developed is based on a parabolic equation and uses the split-step Fourier method for its numerical solution. Using the model, we estimate the influence of a varying, complex, refractive index of the atmosphere, different lossy materials in the ground, terrain, and oceans. The terrain is described using a piecewise linear shift map method. The mod-elling of the ocean is done using a novel model which is a combination of terrain for large swells and Miller surface roughness for smaller waves, both based on a Pierson-Moskowitz sea spectrum. The model is validated and found to agree very well, with results found in the literature.

Nyckelord

Keywords Electromagnetic fields, wave propagation, radar, parabolic equation, split-step Fourier method

(4)

(5)

(6)

(7)

Abstract

Radars are used for range estimation of distant objects. They operate on the prin-ciple of sending electromagnetic pulses that are reflected off a target. This leads to the propagation of electromagnetic waves over large distances. As the waves propagate, they are affected by several aspects that decrease the performance of the radar system. In this master thesis, we derive a mathematical model that de-scribes electromagnetic propagation in the troposphere. The model developed is based on a parabolic equation and uses the split-step Fourier method for its nu-merical solution. Using the model, we estimate the influence of a varying, com-plex, refractive index of the atmosphere, different lossy materials in the ground, terrain, and oceans. The terrain is described using a piecewise linear shift map method. The modelling of the ocean is done using a novel model which is a combination of terrain for large swells and Miller surface roughness for smaller waves, both based on a Pierson-Moskowitz sea spectrum. The model is validated and found to agree very well, with results found in the literature.

(8)

(9)

Acknowledgements

First, I would like to thank my supervisor Robert Jonsson and Saab Dynamics for the possibility to do this master’s thesis. Thank you, Robert, for all discussions, great advice and support over the course of these weeks. Your help has been in-dispensable. I would also like to thank the rest of the staff at Saab Dynamics that I have come into contact with during this thesis. Thank you for being so patient when answering all my questions. I am also very grateful for all the thoughtful feedback that I have received from Tobias Hansson, my supervisor at the Depart-ment of Physics, Chemistry and Biology, both on the report and on numerical simulation in general.

I would also like to thank Peter Holm at FOI for very inspiring discussions about pe-modelling. My understanding would have been greatly reduced without your help.

This thesis marks the end of six years of studies in Linköping. I want to take this opportunity to thank all my friends who have made these years an unforgettable time. It is thanks to all of you that these years have been truly awesome.

Finally, I would like to thank my family for always supporting and believing in me. You have played a larger part in the making of this thesis than you know.

Linköping, June 2019 Jonas Ehn

(10)

(11)

Notation

Physical Quantities

Symbol Quantity Unit

ω Angular frequency [rad s−₁

]

ϵr Complex relative permittivity [1] x Distance along the surface of the Earth [m] J Electric current density [A m−2] D Electric displacement field [C m−2] E Electric field [V m−1_]

pz Fourier space variable [1]

f Frequency [Hz]

γ Grazing angle of the em-waves [rad]

z Height above the surface of the Earth [m]

ha Height of the antenna [m]

H Magnetic field [A m−1_]

B Magnetic flux density [T]

M Modified refractivity [1]

µ Permeability [H m−₁

]

ϵ Permittivity [F m−1]

F Propagation factor [1]

σc Radar cross section [m2] ae Radius of the Earth [m]

u Reduced wave function [A m−1_{] or [V m}−1_]

n Refractive index [1]

N Refractivity [1]

ρ Roughness reduction factor [1]

ρc Charge density [C m−3]

Te Temperature [K]

ψ Wave function [A m−1] or [V m−1]

k Wave number of the em-waves [m−1]

κ Wave number of the sea waves [m−₁

]

(14)

Mathematical Notation Symbol Meaning

[A, B] Commutator between A and B, [A, B] = AB − BA

∂ Differential operator

C _{Fourier cosine transform} F _{Fourier transform} S _{Fourier sine transform}

j Imaginary unit

I0 Modified Bessel function of the first kind of order 0

∇ _{Nabla operator}

A Phasor, A = ℜ[Ae−jωt]

u˜(pz) Transform of u(z)

tri(x) Triangle function, tri(x) = max (1 − |x|, 0) ea Unit vector in direction of a

A Vector

∇2 _{Vector Laplace operator}

Abbreviations

Abbreviation Meaning

dmft _{Discrete Mixed Fourier Transform} em _{Electromagnetic}

fft _{Fast Fourier Transform}

itu International Telecommunication Union mft Mixed Fourier Transform

pe Parabolic Equation pec Perfect Electric Conductor ssfm Split-Step Fourier Method

(15)

1

Introduction

This is a master thesis that deals with electromagnetic (em) fields that are prop-agating in the atmosphere. The context is of that of radar but other applications, such as communications, are also applicable. The thesis work has been performed at Saab Dynamics in Linköping in cooperation with Linköping University. This chapter includes a brief background, the aim of the study and some limitations. Last in the chapter is an outline of the thesis.

1.1 Background

Radar is a tool that is used to detect and determine the distance to distant objects. A radar builds upon the principle that an antenna emits an em-pulse, and then listens for the echo that is returned when the pulse is reflected off a target. In the case of a ground-based radar, the em-waves propagate through the atmosphere in close proximity to the Earth, which means that they are heavily influenced by both the atmosphere and reflections off the surface of the Earth.

How em-waves propagate has been known for a long time [1]. What makes this a complicated problem is that the refractive index of the atmosphere depends on the conditions of the atmosphere, such as altitude and current weather [2]. The electromagnetic properties of the Earth are highly dependent on the type of surface, such as dry soil, snow, water, etc. [3, 4]. These dependencies mean that the performance of a radar system varies depending on where, and under which conditions, it is deployed. A radar system on a ship might behave very differently close to shore than compared to over an open ocean. A commander of a ship might therefore be interested in finding out how the radar system of the ship behaves under the current conditions. Another example can be where to best place a radar or communications installation due to the varying geographical

(16)

and meteorological conditions at different prospective sites. It is therefore a field that has been thoroughly studied from the beginning of radio communications and the invention of the radar [2, 3].

Today, there exist simulation software dedicated to doing these kinds of simula-tions [5, 6]. However, there are some flaws with this software. The Advanced Propagation Model [5] is a very advanced simulation software that has many interesting features. But the fact that the source code is not publicly available means that it is impossible to actually determine how it works. The software tool petool _{by Özgün et al. [6] has a source code that is publicly available, but its} description of a sea surface is rather rudimentary.

1.2 Aim

The aim of this study is to develop a model that can be used to calculate the propagation of radio waves in a complex environment. This model shall include that. . .

• . . . the em-waves are diffracted after being reflected off the Earth

• . . . the refractive index of the atmosphere is not constant, but a function of the altitude and current weather conditions

• . . . the em-waves are refracted in the atmosphere, and can therefore not be expected to travel in a straight path

• . . . the Earth is curved

• . . . it is possible that the waves propagate over a sea surface that is rough with water-waves. The model shall be able to handle up to sea state 6, Code table 3700 in [7].

The domain of computation shall be up to some hundreds of kilometres in dis-tance and up to a kilometre in altitude. The model shall be able to operate for a frequency range of 1–20 GHz. The model shall be implemented as a matlab program that can be run on a personal computer.

If possible, the model should also be able to model the propagation in a domain that contains some terrain.

The aim is to construct a model that includes support for all these features. A literature study shall be conducted to find a suitable approach to these features and their integration. If there are several methods present in the literature, the most suitable shall be chosen.

1.3 Delimitations

The properties of the radar system are not evaluated. This is not in the scope of this thesis. The results from the model should therefore be independent of the radar system used. The radar system is only specified as working at a certain frequency and with the beam pattern of the antenna in the far field.

(17)

1.4 Layout 3 The fact that the radar system is not evaluated means that no exotic signals are considered. This would lead to signal processing considerations that are outside of the scope of this thesis. The em-field emitted by the antenna is supposed to be sufficiently well behaved. That is, there are no extreme em-fields that require consideration of the near-field, non-linearities etc.

The simulations take place in a 2D-environment. This is an approximation that introduces some error in the model since scattering around objects in the hori-zontal plane is neglected. But this error is supposed to be of little importance in most real radar applications. The model can be applied in each direction in the horizontal plane from a point of interest to get the correct behaviour in all directions.

1.4 Layout

Chapter 2gives the theory that is used to build up the model. It describes each feature in section 1.2, above and how they can be approached in a working model. It contains details on the em-propagation and refractive conditions. It also gives a thorough derivation of the governing equations.

Chapter 3contains a brief overview of the work that has been done in the field of em-propagation at large ranges. It contains a review of some of the parts that have to be investigated in greater detail in order to get a good model of long-range propagation. It is concluded with a discussion of some other software comparable to the one developed in this study.

Chapter 4gives an outline of the current model and how it relates to the theory in the theory chapter. It also gives an overview of which simulations that are performed to give the results in the results and discussion chapter.

Chapter 5presents an overview of each of the features of the model, such as varying refractive index, terrain etc., that are derived in this study. The features are presented one by one with an analysis. Each feature is compared to the results presented by other authors to validate the model developed here. It also contains some novel results.

Chapter 6concludes the report by commenting on how the model answers to the goals set up in section 1.2. The chapter also contains some suggestions for future research.

(18)

(19)

2

Theory

This chapter introduces some of the key elements of the theoretical foundations of the present work to the reader. The chapter starts with relating the radar sys-tem to the field of electromagnetics and how the parabolic equation is introduced. Much of the chapter is then dedicated to the description of the split-step Fourier method, ssfm. Afterwards follow descriptions of the refractive conditions in the atmosphere and boundary conditions.

2.1 Radar

A radar system is used for the detection, and the determination of the distance to, distant objects. The name radar is an acronym of the words radio detection and ranging. The working principle of a radar system is that it sends out an em-pulse towards a target and then listens for the echo that is reflected off the target. This can be seen in the schematic drawing in Figure 2.1. A radar system consists of a transmitting antenna that sends out a signal, a signal generator, a receiving antenna and equipment for signal processing. It is common that the transmitting and receiving antenna are the same, such a radar is called a monostatic radar [3].

As the electromagnetic pulse travels to the target, it interacts with what is in its way, usually the atmosphere and the underlying terrain. The em-waves are subject to reflection, refraction and absorption. This interaction can be described by electromagnetic theory and the electrical properties of the terrain and the atmosphere.

(20)

Transmitted wave Re ected wave Radar antenna Radar target

Figure 2.1: A schematic drawing of a radar system, with the radar antenna to the left and a target, in this case an aircraft, up right.

The performance of a monostatic radar system is most commonly described by the radar equation [3]

Pr =

PtGAeσc

(4π)2_R4F

4_, _(2.1)

where Pris the power received at the antenna in Watts, Ptis the power transmitted

by the antenna in Watts, G is the antenna gain in decibel which describes how much of the signal that is radiated in a given direction, Aeis the effective antenna

aperture in m2which is a measure of how large area the antenna uses to capture the signal, σcis the radar cross section of the target in m2 which describes how

much of the incident radiation that is reflected off the target, R is the distance between the antenna and the target, and F is the propagation factor. The equation describes two-way propagation.

The propagation factor, F, can be said to describe everything in the environment that affects the propagation of the waves, so it includes atmospheric effects, shad-owing and reflections off objects etc. It is defined as [2, 3]

F := _||E|

E₀|, (2.2)

where E is the actual field, in the presence of an atmosphere, terrain etc. and E0

is the corresponding field in free space, i.e. in the total absence of atmospheric effects and influence of any objects. It is therefore common to use the propagation factor to describe how the environment affects the system. If the propagation factor is known, it is possible to obtain the actual field up to a phase factor by multiplying the propagation factor with the corresponding free space field. It is also possible to define the propagation factor in terms of path loss, which

(21)

2.2 Electromagnetic waves 7 is the loss along the propagation path, as is done in [6, 8]. Keeping the notation consistent, this gives

F2= Lp

Lf p

, (2.3)

where Lp is the actual path loss and Lf p is the free space path loss. The square

comes from the fact that the path loss is related to the power rather than the field intensity. The path loss is a ratio of how much of the power radiated in one point that arrives at another point.

One common approximation regarding the propagation factor is to consider the propagation to be reciprocal, that is, the propagation effects are assumed to be the same for both the transmitted and the reflected pulse. This means that it is pos-sible to calculate F for only the transmitted pulse, to get a description of one-way propagation, and then just square it to get a description of two-way propagation. Reciprocity is not true in general, but it is a convenient approximation consider-ing how much it simplifies the calculations.

2.2 Electromagnetic waves

The electromagnetic field is a vector field that consists of the two fields, E and H, which are the electric and magnetic fields, respectively. These fields vary over both space and time, so at a certain point r at time t, the fields are E(r, t) and H(r, t). The interaction of these fields in space and time are described by Maxwell’s equations [9] ∇ ·_{D(r, t) = ρ}_c_{(r, t)} _(2.4a) ∇ ·_{B(r, t) = 0} _(2.4b) ∇ ×_{E(r, t) = −}∂B(r, t) ∂t (2.4c) ∇ ×_{H(r, t) = J(r, t) +}∂D(r, t) ∂t , (2.4d)

where J is the electric current density and ρc is charge density. The fields D

and B are related to the electric and magnetic fields, respectively, through the dependence of the material in which the em-field is propagating and are called electric displacement and magnetic flux density, respectively. These equations describe all classical electromagnetic phenomena. They can be solved together with suitable boundary conditions to get a unique solution [10].

This thesis has the delimitation that it only considers em-fields that are suffi-ciently well-behaved. This means that it is, in general, possible to make the as-sumption that the em-field is time harmonic, e.g. that it has a time dependence of cos (ωt). Then it is possible to define the phasors corresponding to the fields as E = ℜ[Ee−jωt], where E is known as the phasor. Using the phasor notation for all field variables, Maxwell’s equations for a time-harmonic field can be written

(22)

as [9] ∇ ·_D_{= ρ}_c _(2.5a) ∇ ·_B_{= 0} _(2.5b) ∇ ×_E_{= jωB} _(2.5c) ∇ ×_H_{= J − jωD.} _(2.5d)

2.2.1 Em-propagation in vacuum

The simplest case of propagation is the case of propagation in vacuum, since there are no currents nor charges to take into account, so that J = 0 and ρc = 0.

This means that the constitutive relations become D = ϵ0Eand B = µ0H. This

makes it possible to write the Maxwell’s equations in vacuum as [10]

∇ ·_E_{= 0} _(2.6a)

∇ ·_H_{= 0} _(2.6b)

∇ ×_E_{= jωµ}₀_H _(2.6c) ∇ ×_H_{= −jωϵ}₀_E_. _(2.6d)

It is usually quite difficult to solve Maxwell’s equations in general, but in many cases, it is possible to simplify the problem. Since em-waves are waves, they satisfy the Helmholtz equation. The Helmholtz equation can be obtained from Maxwell’s equations via a few simple steps [10]. Taking the rotation of Equa-tion 2.6d and using EquaEqua-tion 2.6c on the right hand side gives

∇ ×_{(∇ × H) = ∇ × (−jωϵ}₀_{E) = −jωϵ}₀_{(∇ × E)} = −jωϵ0(jωµ0H) = ω2ϵ0µ0H.

(2.7) For the left hand side, it is possible to use the vector identity ∇ × ∇ × V = ∇(∇ · V) − ∇2V, for any vector V. This gives

∇ ×_{(∇ × H) = ∇ (∇ · H) − ∇}2_H_{= −∇}2_H _(2.8) due to Equation 2.6b. If we define k = ω√ϵ0µ0, we get the standard form of the

Helmholtz equation. The same calculations can also be performed to obtain the Helmholtz equation for the E field. Then the two Helmholtz equations are [9]

(︂

∇2_{+ k}2)︂_H_{= 0} _(2.9a) (︂

∇2_{+ k}2)︂_E_{= 0.} _(2.9b) These are known as the homogeneous vector Helmholtz equations [9]. The inter-pretation of k is as the magnitude of the wave vector k. These equations have the transverse plane wave solution [9]

E(r) = E0ejk·res (2.10a)

H(r) = H0ejk·rek×s (2.10b)

(23)

2.2 Electromagnetic waves 9 the E-field. The polarization vector es is always perpendicular to the direction

of propagation k. This means that E always is perpendicular to H and both the fields to the direction of propagation k. So, the three vectors, E, H and k are all perpendicular to each other [9].

The magnitude of the E and H-fields are related through the intrinsic impedance of the material, η. The impedance is defined as the ratio of the magnitudes of the electric field and its corresponding magnetic field [4]

η := E0 H0

=√︃ µ0

ϵ0

≈_120π _[Ω]. _(2.11) Here it can be seen that the fields only differ by a real factor. The fields are therefore oscillating in phase with each other, but with different amplitudes.

2.2.2 _{Em-propagation in lossy media}

All the actual materials in this thesis are modelled to be simple media, that is, they are supposed to be linear, isotropic and homogeneous. The air is actually not homogeneous from an electromagnetic perspective, but the deviations are so small that they are of no concern to the general discussion here. These as-sumptions lead to the constitutive relations being D = ϵE and B = µH, for scalar

ϵ = ϵrϵ0 and µ = µrµ0. For non-magnetic materials µ = µ0, which applies for

most materials in this thesis. We also make the assumption that currents, J = σ E, may be present in these materials. This gives Maxwell’s equations as [9]

∇ ·_E_{= 0} _(2.12a)

∇ ·_H_{= 0} _(2.12b)

∇ ×_E_{= jωµH} _(2.12c)

∇ ×_H_{= J − jωϵE = (σ − jωϵ) E.} _(2.12d) The last equation looks quite different from the vacuum case, the rest is basically exchanging ϵ0and µ0with ϵ and µ. The last equation can be brought back to the

form we are used to via the introduction of a complex permittivity ϵc[9] ϵc:= ϵ + j

σ

ω. (2.13)

Using the complex permittivity, it is possible to write the last of the Maxwell’s equations as

∇ ×_H_{= −jωϵ}_c_E_, _(2.14) which is the same form as we had in the case of vacuum, [9].

The principle of duality says that Maxwell’s equations are invariant to the kind of linear shifts done in Equation 2.14 in a simple medium [9]. This means that it is possible to obtain the Helmholtz equations in the same manner as above, with the only difference being that the permittivity is now complex. To not change our previous definition of k, we define the refractive index n as n :=√ϵrµr which

(24)

the Helmholtz equations as (︂

∇2_{+ k}2_n2)︂_H_{= 0} _(2.15a) (︂

∇2_{+ k}2_n2)︂_E_{= 0} _(2.15b) and their solutions as

E(r) = E0ejnk·res (2.16a)

H(r) = H0ejnk·rek×s. (2.16b)

The difference here is the fact that the quantity nk is complex. It is common to divide this number into its real and imaginary part to help the discussion, as

kn = ar+ jbi. (2.17)

Then the exponent in the solution of the Helmholtz equation can be split up using

ar and bi. So, the solution can be written using this split as

E(r) = E0ejarek ·_r

e−biek·r_e

s (2.18)

where the term e−biek·r_{is an attenuation term since it decreases with increasing}

r. It is common to define some metric of how fast the waves attenuate in the material. This is usually done with the help of the skin depth which is defined as [10] δs:= 1 b = √︄ 2 ω2_ϵµ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ √︃ 1 + (︃ _σ ϵω )︃2 −₁ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ −1 2 . (2.19)

The skin depth is the distance in meters the field has to propagate in a material before its amplitude has decreased with a factor 1/e.

The intrinsic impedance of a lossy medium is defined as above, but now the ratio becomes a complex number

η =√︃ µ ϵc = √︃ µ ϵ + jσ_ω, (2.20)

due to ϵc [9]. This means that the ratio between the magnitudes of the electric

and magnetic field can be seen as a real part that is a scaling between the two fields, but also an imaginary part that can be seen as a phase difference. This means that the electric and magnetic fields will no longer oscillate in phase [9].

2.2.3 Boundary conditions

At the interface between two different materials, Maxwell’s equations give some boundary conditions that have to be fulfilled by the field on both sides of the boundary. The situation at the interface can be seen in Figure 2.2. The boundary

(25)

2.2 Electromagnetic waves 11 conditions are [10] ϵ1E ⊥ 1 = ϵ2E ⊥ 2 (2.21a) E∥₁ = E∥₂ (2.21b) µ1H ⊥ 1 = µ2H ⊥ 2 (2.21c) H∥₁ = H∥₂ (2.21d) where E⊥is the electric field component perpendicular to the interface, E∥is the electric field parallel to the interface, H⊥ is the magnetic field component per-pendicular to the interface and H∥is the magnetic field parallel to the interface. Several things follow from this. One of the most important is Snell’s law [9, 10]

sin θt sin θi = n1 n2 = η2 η1 =√︃ ϵr1 ϵr2 . (2.22)

which says that the angle of refraction, θt, is related to the ratio between the

refractive indices of the two materials. It is assumed that both materials are non-magnetic so that µ1 = µ2 = µ0. It is Snell’s law that causes refraction in the

atmosphere when the refractive index of the air varies, see section 2.5. Another thing that follows are the Fresnel equations of reflection [9]

Γ⊥=

E0r

E0i =

η2cos θi−η1cos θt

η2cos θi+ η1cos θt (2.23a)

Γ∥= E0r E0i = η2cos θt−η1cos θi η2cos θt+ η1cos θi . (2.23b)

where Ei and Erare the magnitudes of the incident field and the reflective field,

medium 1,η1 medium 2,η2 e_y ex ez E_i θi γ E_t θt E_r

Figure 2.2: Wave reflection and transmission at an interface between two media.

(26)

respectively as indicated in Figure 2.2. The quantity Γ relates the reflected field to the incident field. The difference between Γ⊥ and Γ∥ is that they correspond

to the cases where the incident field is polarized perpendicular or parallel to the plane of the interface, respectively. So, in the first case the polarization vector es⊥is perpendicular to the normal of the interface en = ezin Figure 2.2, so that

es⊥·en= 0. In the case of parallel polarization, es∥·en≠ 0. The Fresnel equations of reflection are not of much use to the current study in their present form. It is preferably to be able to express them in terms of the grazing angle, γ, which is the complementary angle to the incident angle θi (see Figure 2.2) and the medium

impedance η. One then finds that [2, 11]

R+0 = Γ∥=

Y sin γ −√︁Y − cos2_γ

Y sin γ +√︁Y − cos2_γ (2.24a)

R−₀ = Γ⊥=

sin γ −√︁Y − cos2_γ

sin γ +√︁Y − cos2_γ (2.24b)

where Y = η₁2/η₂2 = n₂2/n2₁ = ϵr1/ϵr2, and where for the second equality it is

as-sumed that µ1= µ2= µ0for a non-magnetic media [9]. The reflection coefficients

can be used to calculate the field that is reflected off the interface and will be used later in the boundary conditions and the implementation of ground in the model. The explicit notation for phasors, E, is dropped after this section for notational simplicity. All quantities related to the fields E and B are still considered to be phasors unless explicitly noted.

2.3 The parabolic equation

Solving the Helmholtz equation is, in general, an easier problem than solving Maxwell’s equations, but a quite difficult problem nevertheless. To make this easier, it is common to do some additional approximations and solve a simplified problem instead. In the field of propagation of radio waves, the most common is to use a parabolic equation to describe the propagation.

The parabolic approximation was first described by Leontovich and Fock [12]. It did not become popular until an algorithm, the ssfm, to efficiently solve it was presented by Hardin and Tappert [13] in 1973 and the computational power in computers had advanced sufficiently.

This section presents the derivation of the parabolic equation for both a flat and a spherical Earth. The situation that we want to calculate can be seen in Figure 2.3.

2.3.1 Flat Earth

In the case of a flat Earth, the coordinate system seen in Figure 2.3 is an ordinary Cartesian one, since the coordinate exalong the surface of the Earth is a straight

line. The derivation of the parabolic equation starts with the Helmholtz equation obtained above

(︂

(27)

2.3 The parabolic equation 13

ey ex

ez

paraxial direction

Figure 2.3:The general calculation domain. An antenna is placed to the left in the domain and it emits electromagnetic radiation that travels to the right. The domain is limited by the ground below which is indicated with the curvy line in the figure. Theez-axis is perpendicular to the surface of the Earth, and theex-axis follows the average surface of the Earth. The coordinate system thus corresponds to a Cartesian system in the case of a flat Earth, but not in the case of a spherical Earth.

where k = |k| is the magnitude of the wave vector, n is the refractive index and ψ is a scalar wave function. The scalar wave function ψ corresponds to the compo-nent Eyof the electric field E in the case of horizontal polarization, since it is the

only non-zero component of the electric field in that case. In the case of vertical polarization, ψ corresponds to the component Hy since it is the only non-zero

component of the magnetic field H in that case [8]. The implication of this is that we are always considering a field that is oscillating in the xy-plane.

It is assumed that it is only necessary to describe the propagation in the xz-plane, thus we separate the perpendicular field dependence and drop the y-term to get

(︄ ∂2 ∂x2 + ∂2 ∂z2 + k 2_n2 )︄ ψ = 0. (2.26)

The approximation can be done if we assume that the field ψ is independent, or at most linearly dependent, of y, which is not true but quite close to how a real radar works by looking in one direction at a time. This approximation is quite a large one since it means that we will neglect any scattering around objects and only consider scattering that occurs over objects. It also reduces the domain of calculation from 3d to 2d. The 3d domain can be reconstructed by doing this calculation in all directions and then recombining them to get the volume solution. However, the effects due to scattering from around objects will still be missing.

Equation 2.26 is satisfied by the field component ψ if the refractive index n is constant over the entire domain. Unfortunately, n is often varying with both altitude and range. (See section 2.5 for a complete description of n.) However, the equation is still a good approximation if n varies slowly with respect to x and

(28)

We introduce the reduced wave function

u(x, z) = ψ(x, z)e−jkx. (2.27)

This is done by assuming that the energy is mainly propagating in the x-direction, and, by doing this change of variable, rapid oscillations due to the carrier wave of the field are factored away [14]. This gives a field u that is varying slowly with respect to x [8].

Substituting u, defined in Equation 2.27, into Equation 2.26 and performing the differentiations on the exponential gives

[︄ ∂2 ∂z2 + ∂2 ∂x2 + 2jk ∂ ∂x+ k 2(︂ n2−₁)︂ ]︄ u(x, z) = 0. (2.28)

To further simplify the expression, it is assumed that ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ∂2_u ∂x2 ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ≪ ⃓ ⃓ ⃓ ⃓ ⃓2jk ∂u ∂x ⃓ ⃓ ⃓ ⃓ ⃓ (2.29)

which is known as the paraxial or the parabolic approximation. This approxi-mation is part of the method to simplify the expression. The approxiapproxi-mation is done on the grounds that the reduced function is a slowly varying function in x [15]. The approximation means that the solution will be limited to fields that are propagating close to the so-called paraxial direction [14]. The paraxial direction is the main direction of the propagation of the waves, in this case, the x-direction, see Figure 2.3.

The paraxial approximation lets us remove the small second derivative with re-spect to x, and the standard parabolic equation

[︄ ∂2 ∂z2 + 2jk ∂ ∂x+ k 2(︂ n2−₁)︂ ]︄ u(x, z) = 0. (2.30)

is obtained which is often used as the starting point of many authors, e.g. [16, 17]. At this stage it is noted that this equation is a first order differential equation with respect to x, something that will become crucial when attempting to solve it.

2.3.2 Round Earth

In the case of a spherical Earth, things are slightly more complicated since the e_x-axis in Figure 2.3 no longer corresponds to a straight line in reality. Therefore, it is necessary to do a transformation to obtain a rectangular domain. The deriva-tion follows the same steps, except for a coordinate transform to account for the spherical Earth.

We start by imagining the Earth as a sphere with radius ae, on which we place the

radar antenna at the top. If we introduce an ordinary spherical coordinate system (r, θ, φ) with the origin at the centre of the Earth, the antenna will be placed in the point (ae+ ha, 0, φ). The origin of the domain in Figure 2.3 corresponds to

the point (ae, 0, φ) in this coordinate system. It is possible to set up Maxwell’s

equations for this situation, and from there derive the Helmholtz equation if we assume azimuthal independence. However, we want to represent the Helmholtz

(29)

2.3 The parabolic equation 15 equation in the coordinate system in Figure 2.3 where the coordinates have an easy interpretation. This can be done using a conformal transformation if we write the planes xz and rθ as complex variables using

Ξ= x + jz, ζ = (ae+ h)(sin (θ) + j cos (θ)). (2.31)

The Möbius transformation corresponding to the relation between the two coor-dinate systems is then [18]

Ξ= 2ae ae+ jζ ζ + jae

, (2.32)

where aeis the radius of the Earth.

It is then possible to approximate the Möbius transformation using these rela-tions between the coordinate systems [8]

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ x = aeθ z = aeln (︄ 1 + h ae )︄ ≈_h (2.33)

where ae is the radius of the Earth and h the height above the surface of the

Earth. The coordinate θ comes from the spherical coordinate system. The last approximation is valid for the region where z, x ≪ ae[18].

This gives that the Helmholtz equation, Equation 2.26, can be transformed and written in the new coordinate system, seen in Figure 2.3, as [8, 18]

(︄ ∂2 ∂x2 + ∂2 ∂z2 + ⃓ ⃓ ⃓ ⃓ ⃓ dζ dΞ ⃓ ⃓ ⃓ ⃓ ⃓ 2 k2n(x, z)2 )︄ ψ = 0 (2.34)

where the scalar field ψ relates to the electric and magnetic fields as [8, 18]

ψh= √︄ kaesin (︄ x ae )︄ exp (︄ z 2ae )︄

Eφ, Horizontal polarization (2.35a)

ψv = √︃ ϵ0 ϵ √︄ kaesin (︄ x ae )︄ exp (︄ z 2ae )︄ Hφ, Vertical polarization. (2.35b)

We introduce a reduced wave function u here as well to obtain [︄ ∂2 ∂z2 + ∂2 ∂x2 + 2jk ∂ ∂x+ k 2 (︄⃓ ⃓ ⃓ ⃓ ⃓ dζ dΞ ⃓ ⃓ ⃓ ⃓ ⃓ 2 n2−₁ )︄]︄ u(x, z) = 0, (2.36)

which is very similar to Equation 2.28 in the previous subsection. As before, the next step is to do the parabolic approximation, but here we also do the approxi-mation that [18], k2 (︄⃓ ⃓ ⃓ ⃓ ⃓ dζ dΞ ⃓ ⃓ ⃓ ⃓ ⃓ 2 n2−₁ )︄ ≈_k2 (︄ n2−_{1 +}2z ae )︄ =: k2(m2−_1), _(2.37) where we have introduced a modified index of refraction m2 = n2+ (2z)/aethat

(30)

includes a term that accounts for the curvature of the Earth. The exact require-ments for this approximation can be found in [18]. Using the notation m gives the parabolic equation for the case of the spherical Earth on the same form as for the flat Earth

[︄ ∂2 ∂z2 + 2jk ∂ ∂x + k 2(︂ m2−₁)︂ ]︄ u(x, z) = 0, (2.38)

but with a different interpretation of u due to the new interpretation of the scalar wave function ψ in Equation 2.35. This is the version of the parabolic equation that will be used for the rest of the thesis.

2.4 The split-step Fourier method

The split-step Fourier method was presented as an algorithm intended to solve the parabolic wave equation by Hardin and Tappert [13] in 1973. The method is a marching method which means that it successively calculates the next step from the previous one. I.e., given the field at an initial point, u(x = x0, z), the ssfm

enables the calculation of the entire domain by calculating first u(x = x1, z) and

then u(x = x2, z) etc., up to some maximum x-value, given that the steps xi+1−xi

are small enough.

What makes this method special is that half of the marching takes place in the spatial domain, while the other part takes place in the Fourier domain. So the solution is advanced in half steps. This will become clearer at the end of this section.

The following section is dedicated to the derivation of a propagator, or an evolu-tion operator, that is used to march the soluevolu-tion forward. This secevolu-tion contains two different methods to get to such a propagator under different approxima-tions. The different methods used to obtain them causes the propagators to be slightly different and to have different precision.

2.4.1 Method 1

The method to obtain a propagator follows the outline presented by Agrawal [19]. We start with the parabolic equation Equation 2.38. The parabolic equation can be rewritten in a more convenient form

∂u ∂x = [︄ j 2k ∂2 ∂z2 + jk 2 (︂ m(x, z)2−₁)︂ ]︄ u(x, z), (2.39)

which highlights the fact that it is a partial differential equation of the first order with respect to x.

This version of the parabolic equation can be restated in a more compact form as

∂u

∂x =

[︂

Lˆ + Nˆ]︂u(x, z), (2.40)

with the help of the two operators Lˆ and Nˆ , defined as

Lˆ := j 2k ∂2 ∂z2 Nˆ := jk 2 (︂ m(x, z)2−₁)︂_. _(2.41)

(31)

2.4 The split-step Fourier method 17 Equation 2.40 has a solution that can be expressed using the operators Nˆ and Lˆ. The solution is formally

u(x + ∆x, z) = e

(︂

Lˆ+Nˆ)︂∆x

u(x, z), (2.42)

for a ∆x that is small enough and under the assumption that the operators are both independent of x. This assumption is not entirely true since m in general is a function of both x and z. But for small variations of m, this is approximately true. We assume that m is almost constant, or at least that it varies slowly with respect to x.

The factorization of the exponential operator can be approximated using a Suzuki-Trotter decomposition [20] which can be written

eδAeδB= eδ(A+B)+12δ2[A,B]+O(δ3)_. _(2.43)

We perform the decomposition with the same designations as in [8]

A = 1 k2 ∂2 ∂z2, B = m(x, z) 2₋_1, _{δ =} jk∆x 2 , (2.44) so A and B are Lˆ and Nˆ with the term jk/2 factored out, jkA/2 = Lˆ, jkB/2 = Lˆ. Doing this gives the approximate solution u(x + ∆x, z) as

u(x + ∆x, z) ≈ eδAeδBu(x, z) = e∆xLˆe∆xNˆu(x, z), (2.45) where we have disregarded the quadratic error term exp (0.5δ2_{[A, B]). The error}

introduced when doing this obviously scales with δ2_{, and thereby with (∆x)}2_{, but}

the term δ is not small in general. It is common to have ∆x on the order of 100 m and k ∝ 10 m−1_{for normal radar frequencies. However, it is assumed that A and}

B almost commute, so that the error term becomes small due to the small

com-mutator. The next error term includes nested commutators, [15], and is assumed to be smaller still. It is also possible to calculate the decomposition using a sym-metric Suzuki-Trotter decomposition eδ(A+B) = e(δA)/2eδBe(δA)/2+ O(δ3) [20]. This is more exact, but is for simplicity not used here.

The next step is to calculate the solutions for the operators Nˆ and Lˆ separately. We assume that we have two different regions. One region is in which the wave is propagating through vacuum so that the operator Nˆ becomes 0. The other region is where the refractive effects are taken into account, and there the diffraction operator Lˆ is taken to be zero. So, each marching step is made up of two steps: one with free space propagation and one with only refraction. This can be seen as a system where the wave propagates through a series of thin lenses placed in vacuum [21].

Inside the lens, Lˆ is taken to be zero, then the parabolic equation, Equation 2.40, can be written as ∂u ∂x = Nˆ u(x, z) = jk 2 (︂ m(x, z)2−₁)︂_{u(x, z).} _(2.46)

(32)

This version of the parabolic equation has the solution u(x + ∆x, z) = exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x+∆x ∫︂ x Nˆ dx′ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ u(x, z) = exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ jk 2 x+∆x ∫︂ x m2(x′, z) − 1 dx′ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ u(x, z) ≈_exp(︄ jk∆x 2 [︃ m2 (︃ x +∆x 2 , z )︃ −₁ ]︃)︄ u(x, z), (2.47)

where the integral is due to the fact that Nˆ is dependent of x. The integral should therefore actually be present in both Equation 2.42 and Equation 2.45 as well to account for this. In the third step, the integral has been approximated with the midpoint rule. This approximation should not introduce any large errors even though it is a crude approximation in most cases. It has already been said that

m is almost independent of x since it is almost constant. If it is constant, or if m2would vary linearly, then the midpoint rule is exact, and it can therefore be expected to be reasonably accurate for something that is almost constant. For the case of the free space propagation, where Nˆ is zero and we only consider the operator Lˆ, the parabolic equation, Equation 2.40, reduces to

∂u ∂x = Lˆu(x, z) = j 2k ∂2 ∂z2u(x, z). (2.48)

This equation is most easily solved in the Fourier domain. Taking the Fourier transform with respect to z yields

∂u˜

∂x =

−_jp2_z

2k u˜(x, pz), (2.49)

where u˜ is the wave function in the Fourier domain corresponding to the wave function u in the spatial domain. The equation has the exact solution

u˜(x + ∆x, pz) = exp (︄ −j∆xp 2 z 2k )︄ u˜(x, pz). (2.50)

Now, both these solutions, Equation 2.47 and Equation 2.50, that are valid in each part of space are combined together by the help of Equation 2.45 to get a propagator that takes the solution from one step x to the next x + ∆x. This combined solution is expressed as

u(x + ∆x, z) = exp(︄ jk∆x 2 [︃ m2 (︃ x +∆x 2 , z )︃ −₁ ]︃)︄ F−1 {︄ exp (︄ −j∆xp 2 z 2k )︄ F {_{u(x, z)}} }︄ . (2.51)

(33)

2.4 The split-step Fourier method 19 [8], even though she uses a different derivation.

All the transformations back and forth to the Fourier domain is due to the fact that we only have one propagator, Equation 2.50, in the Fourier domain. The idea is that it is faster to move back and forth than trying to solve the entire parabolic equation in the spatial domain. The fastest way to move back and forth between the Fourier domain and the spatial domain is by using the Fast Fourier Transform, fft. The fft is an algorithm to reduce the number of operations in the discrete version of the Fourier transform, thereby making it faster to compute for larger vector sizes [22].

2.4.2 Method 2

This section presents another method of obtaining a propagator. This method follows the outline of the analysis by Ryan [15] and involves some other approx-imations than the previous one. This results in a propagator that has a different validity and a slightly different expression than Equation 2.51.

The first difference is that we aim to solve the Helmholtz equation for the reduced wave function, Equation 2.28. The first step is to transform the entire equation to the Fourier domain with respect to the z-coordinate. This gives

[︄ −_p2_z₊ ∂ 2 ∂x2 + 2jk ∂ ∂x+ k 2_(m_˜2₋₁₎ ]︄ u˜(x, pz) = 0. (2.52)

Next, we define an operator Wˆ such that

Wˆ (x, pz) =

√︂

−_p2_z _{+ k}2_m_{˜ (x, p}

z)2. (2.53)

It is then possible to rewrite the Helmholtz equation for the reduced wave func-tion, Equation 2.52, with this operator as

(︄ ∂ ∂x+ jk + jWˆ (x, pz) )︄ (︄ ∂ ∂x+ jk − jWˆ (x, pz) )︄ u˜(x, pz) +j[︄ ∂ ∂x, Wˆ (x, pz) ]︄ u˜(x, pz) = 0. (2.54)

The commutator comes from the fact that ∂/∂x and Wˆ does not commute since Wˆ is a function of x through m(x, pz). Assuming an atmosphere that is homogeneous

in range, the commutator becomes zero since m = m(pz) in that case. In reality, m

does vary with range but the variation is considered to be small enough for this to be a reasonable approximation, just as in the previous method, Method 1. It is also clear to see that the left-hand side of this equation consists of one wave propagating in the positive direction and one propagating in the negative x-direction. Considering only the forward propagating wave, and disregarding the commutator, gives

(︄ ∂

∂x+ jk − jWˆ (x, pz)

)︄

(34)

We do the same rearrangement of the equation to highlight it being a first order partial differential equation with respect to x as we did in the previous section.

∂u˜(x, pz)

∂x =

(︂

−_{jk + jW}ˆ (x, p_z₎)︂_{u˜(x, p}_z_{) = jQ}ˆ (x, p_z_{)u˜(x, p}_z_), _(2.56) where we have defined an operator Qˆ

Qˆ (x, pz) := Wˆ (x, pz) − k =

√︂

−_p2_z_{+ k}2_m_{˜ (x, p}

z)2−k. (2.57)

To get a solution that can be marched forward, we want to define an evolution operator Uˆ (x, x0, pz) such that it is possible to calculate the next step given the

previous step as

uˆ(x, pz) = Uˆ (x, x0, pz)u˜(x0, pz), (2.58)

with x0being some value of x smaller than x, x0 ≤x. Substitution of this

expres-sion into Equation 2.56 gives that the evolution operator Uˆ (x, x0, pz) must fulfil

the relation

∂Uˆ (x, x0, pz)

∂x = jQˆ (x, pz)Uˆ (x, x0, pz). (2.59)

This partial differential equation for the evolution operator has the general solu-tion Uˆ (x, x0, pz) = eΩ(x,x0,pz), Ω(x, x0, pz) = ∞ ∑︂ k=1 Ω_k(x, x0, pz), (2.60)

where Ω(x, x0, pz) is a Magnus expansion. The Magnus expansion is an infinite

series, and the first three terms in the sum are [23] Ω1 = j x ∫︂ x0 Qˆ (x1, pz) dx1 (2.61a) Ω2 = −1 2 x ∫︂ x0 dx1 x1 ∫︂ x0 dx2 [︂ Qˆ (x1, pz), Qˆ (x2, pz) ]︂ (2.61b) Ω₃ = −j 6 x ∫︂ x0 dx1 x₁ ∫︂ x0 dx2 x₂ ∫︂ x0 dx3 (︂[︂ Qˆ1, [︂ Qˆ2, Qˆ3 ]︂]︂ +[︂Qˆ3, [︂ Qˆ2, Qˆ1 ]︂]︂)︂ , (2.61c)

where in the last term, the notation Qˆ1 = Qˆ (x1, pz), Qˆ2 = Qˆ (x2, pz) and Qˆ3 =

Qˆ (x3, pz) has been used for brevity. The infinite series is too cumbersome to

han-dle, so it is approximated by its first term. This is a good approximation in most cases since it can be noted that this is no approximation at all if the operator Qˆ commutes with itself for different values of x. This should hold true for Qˆ since Qˆ can be evaluated to be a scalar at each point x in the Fourier domain, and scalars

(35)

2.4 The split-step Fourier method 21 always commute. Applying this assumption gives a short formulation of Ω as

Ω(x, x0, pz) = Ω1(x, x0, pz) = j x

∫︂

x0

Qˆ (x1, pz) dx1. (2.62)

Having expressed Ω in a closed form enables the writing of the solution of the evolution operator Uˆ . The necessary form Uˆ must have to fulfil the requirement in Equation 2.59 is Uˆ (x, x0, pz) = exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ j x ∫︂ x0 Qˆ (x1, pz) dx1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (2.63)

To simplify the evaluation of the integral, we split Qˆ in two parts: Qˆ(x, pz) ≈ Aˆ (pz) + Bˆ(x, pz). This means that the radical in Qˆ has to be approximated

some-how. We choose another method than the one presented by Ryan [15] to do this. The operator Qˆ can be written as

Qˆ = k ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ √︃ −_p2_z k2 + m˜ 2_{(x, p} z) − 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = k(︂√︁1 + Z + ξ − 1)︂ (2.64) with Z = −p 2 z k2, ξ = m˜ 2_{(x, p} z) − 1. (2.65)

Doing a first order Taylor series around ξ = 0 gives √︁ 1 + Z + ξ = √ 1 + Z + ξ 2 √ 1 + Z + O(ξ 2_), _(2.66)

and then doing another Taylor series expansion around Z = 0 in the second term gives √︁ 1 + Z + ξ = √ 1 + Z +ξ 2+ O(Zξ) (2.67) This is the approximation used in [24]. Then we can put Aˆ (pz) =

√

1 + Z = √︂

k2−_p2_z −_{k, and Bˆ(x, p}_z_{) = (k/2)(m}_˜2_{(x, p}_z_{) − 1). Compare these two operators}

Aˆ and Bˆ with the operators Lˆ and Nˆ that were defined in the previous section.

Using these operators, the integral in the exponent in Equation 2.63 can be sepa-rated as x ∫︂ x0 Qˆ (x1, pz) dx1≈(x − x0)Aˆ (pz) + x ∫︂ x0 Bˆ(x1, pz) dx1. (2.68)

(36)

Thus, the evolution operator becomes (with ∆x := x − x0) Uˆ (x, x0, pz) = exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ j∆xAˆ (pz) + j x ∫︂ x0 Bˆ(x1, pz) dx1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (2.69)

We want to split the sum in the exponent to get the right form of the propaga-tor to fit with the split-step algorithm. This is done by the same Suzuki-Trotter decomposition [20] as was done in the first method

exp(︂j∆xAˆ (pz) )︂ exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ j x ∫︂ x0 Bˆ(x1, pz) dx1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≈ ≈_exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ j∆xAˆ (pz) + j x ∫︂ x0 Bˆ(x1, pz) dx1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (2.70)

The propagator Bˆ is easier to calculate in the spatial domain, since n has a mean-ing there. So, we return Bˆ to the spatial domain via an inverse Fourier trans-form. It is assumed that such a transform exists since n represents a physical property and can therefore be considered to be sufficiently well-behaved. There,

Bˆ(x1, z) = (k/2)(m(x1, z)2−1). This causes half of the propagator to be calculated

in the spatial domain, and half in the Fourier domain.

So, inserting the expressions of Aˆ and Bˆ into Equation 2.70, and switching the order of the exponents (that now only are scalars), the full propagator becomes

U (x, x0, z) = exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ jk 2 x ∫︂ x0 (︂ m2(x1, z) − 1 )︂ dx1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ F−1 {︄ exp (︄ j∆x (︄√︂ k2−_p2_z−_k )︄)︄}︄ . (2.71)

It is possible to use the midpoint rule here as well to approximate the integral. This puts the two propagators on the same form, making them easier to compare. Doing this yields a propagator that gives the next step as:

u(x + ∆z, z) = exp(︄ jk∆x 2 (︃ m2(x0+ ∆x 2 , z) − 1 )︃)︄ F−1 {︄ exp (︄ j∆x (︄√︂ k2−_p2_z−_k )︄)︄ F {_{u(z, x)}} }︄ . (2.72)

This propagator is very similar to what is referred to as the wide-angle propagator by Levy [8]. The difference is in the environment term where she obtains m − 1 instead of the current (m2−1)/2. This comes from different approximations of the radical in Qˆ , Equation 2.64.

(37)

2.4 The split-step Fourier method 23 Both propagators that are obtained by the two methods are presented in Table 2.1 for easier comparison. The propagator derived in this section, Equation 2.72, has the same environmental term as the one defined in the previous method, Equation 2.51, what differs between them is the free space term. It can also be seen that the free space term in Equation 2.51 is a Taylor expansion of the rational to the first order of the free space term in Equation 2.72. This is reassuring, since it means that both propagators should describe the same phenomena.

2.4.3 Sampling distance

One of the main steps in these methods is the Fourier transform. This means that the signal will be moved to the Fourier domain and then back. The signal will thus have to be sampled in accordance with the sampling theorem to avoid aliasing and distortion of the signal [25]. This is very important since the recon-struction of the signal occurs several times in each step, meaning that a small reconstruction error each time will accumulate to become a very large recon-struction error in the end.

A function x(t) that has a Fourier transform that exists, and equals zero for |ω| > 2πB is said to be band-limited, where B is its bandwidth. Such a function can be uniquely reconstructed by sampling the signal with a sampling distance of 1/2B, the sampling theorem [26]. The reduced wave function, u(z), is assumed to be band-limited with the bandwidth pmax, since the Fourier transform of u is

assumed to be zero above this value, u˜(pz) = 0 for |pz|> pmax[18]. The sampling

theorem then says that this band-limited function can be uniquely reconstructed from its transform for all points z if the sample points are spaced π/pmaxz-units

apart, or closer. The spacing in z is equal to ∆z = zmax/Nz, where zmax is the

largest z of interest, the top of the domain, and Nz is the number of samples

in the z-direction. But the sample spacing is also, from the sampling theorem,

∆z ≤ π/pmax. So, we get that the minimum number of sample points in the

z-direction is Nz≥(zmaxpmax)/π.

Table 2.1:A comparison between the narrow angle propagator obtained us-ing method 1 and the wide angle propagator usus-ing method 2.

Environmental term:

Narrow angle: exp(︄ jk∆x 2 [︃ m2 (︃ x +∆x 2 , z )︃ −₁ ]︃)︄

Wide angle: exp(︄ jk∆x 2 [︃ m2 (︃ x +∆x 2 , z )︃ −₁ ]︃)︄ Free space term:

Narrow angle: exp (︄ −j∆xp 2 z 2k )︄

Wide angle: exp ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ j∆xk ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ √︃ 1 −p 2 z k2 −1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(38)

The variable pz is interesting. Within the interval −pmax < pz < pmax it can be

interpreted as pz = k sin α, where k is the magnitude of the wave vector and α

is the angle from the x-axis. Then the maximum value pmax can be calculated

as pmax = k sin αmax where αmax is the maximum angle of interest [18]. In this

context, it is also possible to see αmaxas the angle corresponding to the wave with

the highest elevation that can be correctly reconstructed. Using this in the requirement on Nzgives that

Nz ≥ zmaxpmax

π =

zmaxk sin αmax

π =

2zmaxsin αmax

λ , (2.73)

if we use that k = (2π)/λ, with λ being the wavelength of the signal. This is a condition that is easier to fulfil than the earlier expression since | sin α| always is smaller than one regardless of the angle α. This limitation can also be expressed as a requirement on ∆z

∆z ≤ λ

2 sin αmax

(2.74) which is almost always larger than λ in a realistic case. In the parabolic approx-imation, it was assumed that the waves were propagating forward with a small angle relative to the surface. This means that the angle αmaxdoes not have to be

very large to include all the waves, and thus the spacing ∆z can be quite large. An example of 20 GHz and an angle of 5° gives that ∆z should be at least equal to 8.6 cm, which is larger than the wavelength, 1.5 cm.

2.4.4 Initial field

The initial field is of great importance since the split-step Fourier method marches a solution forward. So, the initial field will be the basis for the entire solution. The initial field used in this study is the one described by Levy [8]. It represents a Gaussian beam pattern propagated to the far-field, so that the antenna is outside of the domain, and is written as

uf s(0, z) = kβ 2 √ 2π ln 2exp (−jkθ0z) exp (︄ − β 2 8 ln 2k 2_{(z − h} a)2 )︄ (2.75) where θ0is the elevation angle of the antenna compared with the surface at the

site of the antenna, β is the half-power beam width of the antenna, ha is the

altitude of the antenna over the ground and k is the wave vector of the signal. We want the initial field to satisfy the boundary conditions. This is done by using image theory and assuming a perfectly conducting ground. Then the initial field at x = 0 can be written as u(x = 0, z) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ uf s(0, z) − uf s(0, −z) Horizontal polarization uf s(0, z) + uf s(0, −z) Vertical polarization (2.76) This field is then what is marched out over the entire domain.

(39)

2.5 The atmosphere 25

2.5 The atmosphere

Since the propagation of the em-waves takes place in the atmosphere, the electro-magnetic properties of the atmosphere become important. The most important electromagnetic properties of the atmosphere can be summarized in the index of refraction, n. The index of refraction of a material is a ratio that describes how fast light propagates in that material compared to vacuum. This ratio is given as [27]

n = c0

v , (2.77)

where v is the speed of light in the current medium and c0 is the speed of light

in vacuum. For air, the index of refraction is almost 1, it seldom exceeds 1 with more than a fraction on the order of 1 × 10−4_{. It is usually around 1.0003 for air}

[28]. This makes n quite a cumbersome unit, so the refractivity is introduced as [29]

N := (n − 1) × 106, (2.78)

for a real index of refraction. The refractivity is measured in N-units, [30]. The refractivity of the atmosphere is a function of the temperature, pressure and humidity as described by [28, 29] N = 77.6P Te + 3.73 × 105ew Te2 , (2.79)

where P is the atmospheric pressure in millibars, Teis the temperature in Kelvin

and ew is the partial pressure of water vapour in millibars that is calculated as

[29] ew= esHR= 6.1 exp (︄ 25.22Te−273 Te −_{5.31 ln} (︃ _T e 273 )︃)︄ , (2.80) where es is the saturated water pressure in millibars and HRis the relative

hu-midity in percent.

Most of the time, it is possible to use the approximation of a standard atmosphere rather than having to know all the parameters in the equation above. For a stan-dard atmosphere, the refractivity is a function of the altitude. It decreases with altitude according to [30]

N (z) = 315e−0.136z, (2.81)

where z is the altitude in km above the sea surface. This can be seen in the left part of Figure 2.4. This causes the em-waves to bend downwards in a standard at-mosphere, in accordance with Snell’s law, which is explained in subsection 2.2.3. So, a ray that is launched parallel to the surface of the Earth will bend down slightly. It will bend down with a curvature corresponding to a circle with a ra-dius approximately4_/₃_{of that of the Earth. This is why}4_/₃_{times the radius of}

the Earth, together with straight rays can be used for calculations, since that in-corporates the effects of the atmosphere [3]. Using an Earth radius4_/₃_{times the}

(40)

very common approximation, but it is seldom explained why it is done.

Figure 2.4:VerticalM-profiles, with the standard atmosphere to the left, an evaporation duct in the middle and an elevated duct to the right.

The Earth is curved and it is usually most interesting to see how the em-waves propagate in relation to the curvature of the Earth. It is common to include the curvature of the Earth in the refractivity. The refractivity is modified in such a way that the em-waves are bent upwards with as much as the Earth should have bent away downward. The modified refractivity, M, is defined as [30]

M := (n + z

ae

−_{1) × 10}6_{= N +} z

ae

×₁₀6_, _(2.82) Where z is the height above the local surface of the Earth and ae is the radius

of the Earth (6360 km [27]). This has the same effect as performing a coordinate transformation from a global coordinate system to one that follows the curvature of the Earth. It is more intuitive to have distance as length along the surface of the Earth and altitude as height above the surface of the Earth. This is introduced in the model by adding a factor2z_/_a_e_{to the refractive index term in the}

environ-mental part of the propagator above so that instead of having n2−_{1, we have}

n2−_{1 + 2z/a}_e_{= m}2−_{1 [18].}

The most interesting aspect of M is how it changes with altitude, M′ = dM/dz, since that is what defines the behaviour of the em-waves [3]. Four different cat-egories can be defined. These are standard, sub-refractive, super-refractive and ducting. The corresponding values of M′ can be seen in Table 2.2 and exam-ples of the different beam propagation paths in these conditions can be seen in Figure 2.5.

A sub-refractive condition means that the rays will either bend down less than they do in a standard atmosphere, go straight or even bend upwards. This is said to be a rare condition [3] but can occur when cooler air resides over warmer water

Propagation of Radio Waves in a Realistic Environment using a Parabolic Equation Approach

Propaga on of Radio Waves in a

Realis c Environment using a

Parabolic Equa on Approach

Jonas Ehn

Abstract

Acknowledgements

Contents

Notation

1

Introduction

1.1

Background

1.2

Aim

1.3

Delimitations

1.4

Layout

2

Theory

2.1

Radar

2.2

Electromagnetic waves

2.2.1

Em-propagation in vacuum

2.2.2

Em-propagation in lossy media

2.2.3

Boundary conditions

2.3

The parabolic equation

2.3.1

Flat Earth

2.3.2

Round Earth

2.4

The split-step Fourier method

2.4.1

Method 1

2.4.2

Method 2

2.4.3

Sampling distance

2.4.4

Initial field

2.5

The atmosphere

_{Em-propagation in lossy media}