Electromagnetic Wave Propagation Modeling for Finding Antenna Specifications and Positions in Tunnels of Arbitrary Cross-Section

Electromagnetic wave propagation in confined environments, such as mine, road or railway tunnels, has inadmissible quality due to multi-path and high scattering effects of receivable signals. The ever-increasing work difficulties and the large demand for systems operating in these scenarios, requires more and more innovative wireless systems to guarantee a reliable communication. From the radiating element point of view, some investigations have been given over to the development of novel technologies to provide radio communication in these environments. In 1956, a paper [46] disclosed what is now generally known as the leaky-feeder principle for propagation of VHF or UHF signals through a tunnel. Leaky coaxial cables (LCXs) have been employed since then [18, 38]. Nevertheless, high cost, maintenance and installation difficulties are the main disadvantages of leaky feeders. These problems are easily solved in short tunnels, but in the case of longer tunnels the leaky feeder solution becomes too expensive. For these reasons, solutions that are based on the use of distributed antennas are becoming more interesting [6].


Introduction
Electromagnetic wave propagation in confined environments, such as mine, road or railway tunnels, has inadmissible quality due to multi-path and high scattering effects of receivable signals.The ever-increasing work difficulties and the large demand for systems operating in these scenarios, requires more and more innovative wireless systems to guarantee a reliable communication.From the radiating element point of view, some investigations have been given over to the development of novel technologies to provide radio communication in these environments.In 1956, a paper [46] disclosed what is now generally known as the leaky-feeder principle for propagation of VHF or UHF signals through a tunnel.Leaky coaxial cables (LCXs) have been employed since then [18,38].Nevertheless, high cost, maintenance and installation difficulties are the main disadvantages of leaky feeders.These problems are easily solved in short tunnels, but in the case of longer tunnels the leaky feeder solution becomes too expensive.For these reasons, solutions that are based on the use of distributed antennas are becoming more interesting [6].
Antenna design for these scenarios is a challenge : Current tunnel cross-sections can have arbitrary shapes and antennas have to operate close to tunnel walls.These antennas have to be discreet and small and propose higher operating frequencies.Moreover, possible interference among systems may be suffered due to the growing number of wireless systems working simultaneously.Contrary to free-space applications, the specifications and parameters involved for antennas working in these environments have not yet been thoroughly established.The understanding and modeling of radio-wave propagation are essential steps toward deriving optimum radio communication systems and deployment in these scenarios.Various studies have been done considering several models to find or analyze radiation characteristics and positioning of antennas in tunnels [21,27,28] or simply to model the wave-propagation in these environments [13,15,19,34,41,57].The models used to study the wave-propagation in these environments are limited due to the intrinsic principles, applicability and assumptions on which they are based.Asymptotic methods are mostly used [55].However, asymptotic methods cannot predict well the fields for antennas associated with systems which are strongly affected by the surrounding environment.Near-field considerations have to be accounted for and, thus full-wave methods should be used.
The problem of antenna operation in tunnel environments has been addressed very briefly.From the positioning point of view, the U.K. guidance (GK/GN06602) on train rooftop antenna positioning as well as some manufacturers provide some obstruction, separation, installation and maintenance specifications for train antennas.Alternative solutions offer improvements in the channel capacity by means of the transmit, receive or spatial diversity, i.e., by using multiple antennas in the transmission/reception located at different positions [44,45,56].A modal interpretation of the principles of these solutions are studied in [37,43].The determination of field specifications to excite correct modes to avoid causing interference between systems, can be also analyzed by using the modal interpretation.Due to multi-modal propagation in these environments, different positions and types of excitation are involved.A modal approach is considered to provide a straightforward description of the wave propagation and the degrees of freedom in these environments [1,54].These modes may be combined in some way to produce a desired effect, e.g., the maximization of the transfered power by the modes.Thus, modal-based analysis and optimization techniques for mode-weight adjustment are used to determine correct co-excitation of these modes, and fulfill the desired specifications.This approach is often used in adaptive array theory [26,47].However, to the best of the authors' knowledge, no such treatment has been reported to deal with multi-modal propagation in tunnel environments.This chapter is organized as follows : Section II introduces the modal approach for guiding structures.It is based on a full-wave method, namely the Transmission Line Matrix (TLM) method.These methods has been hampered by their large computational time when compared to asymptotic methods when large size environments are considered.Thus, a suitable 2.5 D TLM implementation to reduce the computational time and to include lossy dielectric walls of tunnels is briefly presented [2].The computation cost is reduced compared to typical solutions by using the concept of Surface Impedance Boundary Condition (SIBC).Section III is devoted to the description of a methodology for the determination of antenna field specifications and positioning in operational scenarios at high frequencies.Section IV presents the validation of this methodology for a simple canonical case.Lastly, section V describes the analysis and results for a real scenario representative of tunnel environments.Finally, discussions and conclusions are developed.

Modeling approach for guiding structures
The high complexity of today's relevant EM problems has popularized the use of numerical approaches to approximate exact solutions [30].The main limitation of these methods is the high computation time when electrical-large structures are considered.With the increasing performance of computers, research in electromagnetic modeling through numerical methods has been continuously advancing.New techniques have merged and the capability of the existing methods has been improved.A modal approach based in the Transmission Line Matrix (TLM) method, is considered to analyze the wave propagation in tunnels.A reduced TLM formulation, is employed to find the mode field distribution of the uniform guiding structure.Adequate boundary conditions are additionally used to to limit the computation volume.The modeling approach can be applied for a tunnel with arbitrary cross section, including lossy walls.

Brief overview
The modal approach is a rigorous technique, allowing us to have an appropriate physical insight of the propagation in tunnel environments.It allows us to express the fields in terms of modes and generalize the applicability of our algorithm for tunnels with arbitrary cross section.This theory is based on functional analysis.The total fields in a guiding structure, traveling in the ±z direction, can be expanded by the sum of a set of linearly-independent functions Ψ.Their mutual orthogonality is not necessarily required.The space spanned by these functions is given by: where E guided and H guided are the total guided-electric and magnetic fields, respectively, N is the total number of modes.Ên ( r) and Ĥn (x, y) are the normalized mode profiles, w n ( r) are the mode weight coefficients.Finally, α and β are the attenuation and phase constants, respectively.Considering a particular guiding structure, these modes are solutions for the fields which satisfy both Maxwell's equations and the boundary conditions of the waveguide.
To fully characterize the modes, field profiles and the attenuation and phase constants have to be determined.A modal approach based on the TLM method is employed to this end.
The TLM is a time-domain explicit algorithm based on local wave decomposition.The computational domain is simulated by a network of interconnected transmission-lines.All fields components are computed at the center of the nodes (connexions) by linear combination of impulse voltages traveling in the transmission-lines at discrete time Δt.The explicit scheme means that voltages at a later time from the system are calculated at the current time.The method is carried out basically through two processes: scattering and connection.In the scattering process, voltage pulses k V i are incident upon the node from each of the link-lines (halfway between two nodes) at each time step k.These pulses are then scattered to produce a set of scattered voltages, k V r , which become incident on adjacent nodes at the next time step k + 1.In the connection process, pulses are simply exchanged among immediate neighbors.
A simplified approach based on this method and by using the procedure described in [39,58], is briefly introduced in the next subsection.

EM Modeling of arbitrary cross-section over-sized waveguides: The 2.5-D TLM approach
The simplification of guiding problems in TLM was first introduced in [29] to obtain the dispersion properties of guiding structures.They proposed the reduction of the calculation region by introducing the field dependence along the propagation direction z, which can be expressed by the terms e −jβz , where β is the phase constant.Two points along the longitudinal 263 Electromagnetic Wave Propagation Modeling for Finding Antenna Specifi cations and Positions in Tunnels of Arbitrary Cross-Section distance of the guide can be characterized by the phase difference β (z 2 − z 1 ).Then, it is possible to find a relationship between the reflected voltages of the node at the time kΔt and the incident ones at the time (k + 1)Δt.The Symmetrical Condensed Node (SCN), shown in Fig. 1, constitutes the most widely used TLM node for 3D structures.The expressions relating to these voltages for this node are given by equations (2,5), where k V r 4 , k V r 2 , k V r 9 , k V r 8 are the reflected voltages at the time kΔt and k+1 V i 4 , k+1 V i 2 , k+1 V i 9 , k+1 V i 8 are the incident voltages at the time (k + 1)Δt respectively, on lines 4, 2 9 and 8.
Note that the characteristic parameter β must be imposed at the beginning of each calculation.Then, the resulting fields correspond to the solution of the fields for this value.Two traveling waves are injected in opposite directions in some of the nodes forming a standing wave.A time-domain signal, such as a Dirac or a Gaussian pulse is usually employed to analyze the structure in a wide-frequency range.After a considerable number of time step iterations, the response is taken at some nodes and correspond to the superposition of the modal fields.By performing a Fast Fourier Transform (FFT), the response in frequency domain can be obtained.The peaks in the spectrum correspond to the modes.Several β values have to be imposed if the response of the modes for different phase constants is desired.Based on this idea, a 2.5 Dimensional SCN for guiding structures formulation was employed to determine the fields [39,58].The term 2.5D is used as the 3D computational domain is reduced to a 2D mesh in the guide cross-section.However, cells are 3D ones and account for all 6 electromagnetic field components.
In addition to the reduced node for guiding structures, the concept of Surface Impedance Boundary Condition (SIBC) is also considered to efficiently model the tunnel walls, while reducing the computational domain, as seen in Fig. 2 and explained below.

SIBC boundary conditions to model the tunnel walls
Tunnels are made with dielectric walls and may contain other materials (trains, cars, objects etc.).For volumic methods, such as TLM, Finite Difference Time Domain (FDTD) or Finite Element Method (FEM), the limitation of the computational domain within dielectric materials is usually achieved by using Perfect Matched Layer (PML) as boundary conditions.However, this technique is too expensive in terms of complexity and consequently computation time, when dealing with very large electrical objects.Since we are interested only in the EM fields within the hollow region of the tunnel, an appropriate boundary conditions, namely, the Surface Impedance Boundary Conditions (SIBC), constitutes the mathematical artifice to limit the computation volume.For detailed description of its formulation, validation and results obtained with the TLM method, the reader is referred to [2,3,58].Here, an introduction of this concept is briefly presented.
Boundary conditions in TLM are simulated by introducing a relationship between incident and reflected voltages at the boundary.Its expression depends on the nature of the materials on both sides.An interface of air and a lossy dielectric, as in the case of tunnel walls, can be simulated by introducing a reflection coefficient Γ, given by [10]: Γ is in general complex and would alter the shape of the excitation pulses, which normally cannot be accounted for in the TLM method [10].V i armj and V r armj represent the incident and reflected voltages of the j-arm of the node, respectively, and Zs correspond to the Surface Impedance of a lossy dielectric [50].Thus, the SIBC concept allows the efficient modeling of lossy dielectrics, replacing one medium by a local reflection coefficient at the interface.This avoids meshing and calculating fields beyond the interface.For tunnel wall modeling applications, these field values are not required for mode parameter computation.However, the presence of the the dielectric beyond the interface is accounted for by the SIBC.This approach is applied in the next sections to study the radio wave propagation for waveguide with arbitrary cross-section.However, it is limited to walls of permittivity e r >> 1 but without any restrictions on the conductivity σ.However, the problem of finding antenna specifications is examined and thoroughly discussed first.

Justification of the choice
One of the most important properties of an antenna is that it radiates or receives radio waves in a way related to field distribution.In the simplest scenario of a wireless propagation channel, radio waves propagate in free space, expand spherically from the TX to the RX antenna.Antennas are then designed according to a desired radiation pattern and considering far fields.A more complicated situation occurs if the radio waves propagate through more complex environments, such as tunnels.Waves find several paths through a complex environment with a variety of several scattering obstacles.Antenna specifications for these environments are needed.
Modal theory is usually employed to study the correct excitation in guiding structures [11,40,51].In the particular case of radio-wave propagation in tunnel environments, the modal theory has been employed in the past.Some of these approaches have considered canonical geometries and excitation sources [15][16][17]43].Thus, in the next subsections, a methodology to find the field specifications and positioning in these scenarios, based on some recent research, is discussed.

Basic principles
In this subsection, essential principles used with this methodology to obtain the field specifications and antenna positioning in confined environments, are described.These principles make use of field modal expansion, the mode orthogonality.They will allow the optimum antenna positioning for mode excitation.

Modal expansion of the fields
The modal expansion (1) is only valid for closed source-free regions.Inside the source region, besides the discrete mode series expansion, an additional term is needed to account for the fields of the sources [4].Additionally, for open waveguides, there is a continuous part of the spectrum related to radiation modes [4,59].Thus, the total fields for wave-guiding structures can be generally expressed by (7): where E sources (z) and E radiated ( r) are orthogonal complements of the functions Ψ E and Ψ H .An analogous expression can be written for the magnetic fields.Inside the source region, the weight coefficients w n ( r) in (1) represent the mode amplitude and take into account the z-dependence of the modes and the exciting sources, if any.They depend on the coordinate z as a result of the source actions, modeling the transverse field components of the sources.However, for source-free regions and assuming that modes are orthogonal, they are constants w n ( r) = w n .The longitudinal components of the source are modeled by the term E sources ( r) [4].
From the point of view of the fields, modes can be excited and, hence, treated separately.However, from the point of view of the total power carried by the modes, they may be coupled.This point will be discussed in the next subsection by means of the quasi-orthogonality relationship in hollow waveguides.

Quasi-orthogonality relation for hollow-lossy waveguides
Mode orthogonality is directly related to the properties of functions defined on Hilbert spaces H.The elements of these spaces are defined with an inner product.Consider a set of functions Ψ 1 , Ψ 2 , Ψ 3 ...Ψ N in a Hilbert space H.The inner product of Ψ k (r) and Ψ l (r) for k not necessarily equal to l is defined by: where * indicates the conjugate.The functions Ψ k (r) and Ψ l (r) are said to be orthogonal if their inner product for k = l is zero.Otherwise they are said to be non-orthogonal.
The existence of the orthogonality property allows us to simplify the analysis of a given waveguide.This property plays a key role in solving problems of excitation and scattering in waveguides [40].The orthogonality between modes means that each mode carries its own power independently of other modes.In lossless waveguides with perfectly electric or magnetic walls, the orthogonality is an inherent property.In waveguides with lossy dielectric walls, such as tunnels, orthogonality no longer holds and so-called quasi-orthogonality can prevail.A complete treatment of this problem can be found in [4].To better understand this consequence, consider the k-th and l-th modes propagating in a source-free region.The curl Maxwell equations applied to these modes are given by: where σ m and σ e are the electric and magnetic conductivities, respectively.Using Maxwell's equations, one finds: Application of the two-dimensional divergence theorem: into (10) yields: where S is a closed surface with interface contour L, Z s denotes the surface impedance, H denotes the tangential component of H at the boundaries and: where γ is the propagation constant and ξ is equal to either k or l.Substitution of (13) into equation ( 12) yields to: 267 Electromagnetic Wave Propagation Modeling for Finding Antenna Specifi cations and Positions in Tunnels of Arbitrary Cross-Section The real (time-averaged) total power transmitted P (z) and dissipated Q (z) can be calculated by equations (15).
Every single mode carries a real self-power flow P k = P kk and self-power loss If self powers are coupled to each other, i.e.Q kl = 0, it leads to the non-orthogonality between modes.Then, equation ( 18) allows us to analyze the cases where the non-orthogonality between modes occurs.It is worth noting that, for a hollow-dielectric waveguide, such a tunnel, the surface integral in Q (z) vanishes, and the expression reduces to: where L represent a closed contour around the surface of the hollow-waveguide.In this case, modes are non-orthogonal due to the fact that fields do not vanish at the boundaries leading to Q kl = 0. Another interesting case occurs when a waveguide with perfectly electric walls and loaded with a lossy dielectric, is considered.In this case, the expression (18) reduces to: The field components are obtained by multiplying the expressions for the fields for an unloaded waveguide with e αz , where α is the attenuation constant [7].Evaluating expression (20) for this case leads to zero.Thus, the modes remains orthogonal even when losses are present in the waveguide.However, the integration surface S can be chosen such that the integral in (20) does not vanish, obtaining the non-orthogonality between modes.This simple case will be considered later to evaluate our methodology for waveguides with non-orthogonal modes.
The direct consequence of the non-orthogonality of modes is that the total power is not additive, i.e., it is incorrect to state that the total power loss is the sum of the power loss of the modes propagating independently [4,59].Thus, the powers of the modes in lossy-hollow-dielectric waveguides, such as tunnels, are in general, mutually coupled and this fact has to be considered.

Optimum positioning for mode excitation
In the most general case, a wave-guiding structure can be excited by a current source or by an external field.Several modes can be excited and the main objective consists in deriving the weight coefficients due to a given source.For lossless waveguides, the orthogonality property implies that there is no power exchange between modes.The formalism for determining these coefficients is well-known and is usually done by considering current sheets [40,51].In practice, such currents are approximated by probes or loops.The orthogonality in this case simplifies the analysis and modes can be selectively excited or excluded.
Nevertheless, in lossy-dielectric waveguides, energy exchange between modes can take place.
The theory that describes the interaction between them is known as the Coupled Mode Theory.This theory allows us to determine the mode amplitudes.The problem of the weight calculation is rather complicated and lies outside the scope of this chapter.The reader is referred on this point to the works described in [4,35,59].The mutual coupling between modes causes the mode amplitudes to be coupled and the resulting coefficients are no longer constant but become dependent on the propagation direction.The desired set of the equations of mode excitation for k = 1, 2, ..., can be written in the form (21). where: S s and C s correspond to the surface and contour for the source, respectively.J e a , J e b , J m a , J m b are electric and magnetic currents to model the source.The subscripts a and b denotes the external bulk and surface sources.The ideal excitation position can be studied by considering two modes, the resulting expression for the mode amplitude of one of the modes is given by: where N 12 and N 21 are the coupling factors between the modes, γ 1 is the propagation constant of the first mode, R 1 and R 2 correspond to the excitation terms (23) for the first two modes.
In [43], it was demonstrated that these values can be relatively small.By comparing ( 24) and ( 23), it can be noted that the ideal excitation for a given mode is predominantly dominated by R 1 .In fact, the dot product in (23) indicates that the excitation should coincide with the mode profile.As a result, to maximize the excitation of a given eigenfunction, the source should be located in regions where it is predominant and avoid those where it vanishes.In a more general case, the eigenfunction of the electric and magnetic field should be maximized.Thus, a proper excitation should be placed in regions where the power is maximum.This analysis may be extended for the expression of the weight coefficient for all modes, which is given by: where w(z) is a vector containing the mode amplitudes, e −γz is a vector containing the exponential terms of the mode propagation constant, N is a matrix of the coupling factors 269 Electromagnetic Wave Propagation Modeling for Finding Antenna Specifi cations and Positions in Tunnels of Arbitrary Cross-Section given by (22) and R is a matrix with the source terms given by (23).Appropriate mode excitation should be located in regions where the sum of the mode field profiles is maximum and avoid those where it vanishes.More generally speaking, locations where the total power carried by the modes maximizes, correspond to those where the electric and magnetic fields are simultaneously better excited.

Optimum field excitation for mode excitation
Field specifications for correct mode excitation has mainly been studied in the literature for applications in optical communications.The selective mode excitation in a hollow waveguide is studied in [42,48,49].A linearly polarized input light is converted to a TE 01 -like profile to excite this mode in a circular waveguide.Several papers on the design of mode converters can be consulted by the interested reader in [24,36].A Gaussian beam launched axially into a circular fiber is used in [20], for exciting low-order modes.Finally, only canonical cases have been reported for the integration of antennas in guiding structures [9,12,14,22,23,51].
An optimum source excitation in dielectric waveguides, such as tunnels, is still to be determined and its characteristics are determined by the nature of the weight coefficients given by (1).As mentioned before, multi-modal propagation in these environments can cause detrimental effects to the received power.A similar problem arises in adaptive array theory where the array output power has to be minimized by choosing adequate weight coefficients of an N-element array.Optimization techniques are employed to find a proper set of coefficients.A complete treatment of this problem can be found in [26,47].Optimum mode coefficients for correct mode excitation in tunnels can be determined by following a similar optimization procedure as in adaptive array theory.The total power carried out by the modes, given by (16), must be optimized.The power must be restricted to a given range, and only finite values of the weight coefficients should be considered.The problem can be stated into different ways.First, by maximizing the transmitted power through M receivable modes in Rx, such that its value is bounded.Its formulation can be expressed mathematically as the maximization of: P (x, y, z) = w H P ΔTx w (26) subject to the constrains: The second choice consist of maximizing the received power through M transmittable modes in Tx, such that its value is higher than a given threshold.The problem can be reduced to the maximization of: subject to the constrains: Equation ( 26) and ( 28) are the matrix form of ( 16) in Tx and Rx.The vector w contains the weight coefficients and w H denotes the conjugate transpose of w.P is a hermitian N × N matrix (i.e.P = P H ) formed from the cross-Poyinting powers of the N modes in the structure, P is a non-diagonal matrix, the (l, k) element of this matrix is given by: Hence, P = P ΔTx or P = P ΔRx correspond to the power densities in Wm −2 for an element of area ΔΩ and represent the evaluation of (30) in Tx or Rx.The terms δ max Tx = P max Tx /Ae Tx and δ min Rx = P min Rx /Ae Rx are the maximum and minimum powers densities in Tx and Rx, respectively.P max Tx and P min Rx correspond to the maximum and minimum powers that can be transmitted and received in an effective area Ae Tx in Tx and Ae Rx in Rx.These power densities will be clear later on.The first constraint in (27) ensures that the power density (or simply the power) is bounded, whereas in (29), it ensures that the power is higher than a given threshold.
In general, only one of these problem has to be solved, however the treatment of both will be considered here for better understanding.The third constraint avoids the solution w = 0 insuring that the weight coefficient vector has a nonzero magnitude and that c is necessarily a real constant.Finding w optimum to satisfy ( 27) can be accomplished by the method of Lagrange Multipliers.The partial derivative of the unconstrained problem (31) with respect to each variable yields a system of equations to find the values of the weight coefficients.However, w is a complex vector and it is not clear how the derivative operates on the real and imaginary parts.A complex gradient operator defined in [5] is rather employed to this end.This operator or its complex conjugate are suitable for determining stationary points of a real function, such as P. The necessary Kuhn-Tucker conditions [53], to maximize the following cost function: are given by: where ∇ ν denotes the complex gradient operator with respect to ν, Ω = Tx or Rx and δ = δ max Tx or δ = −δ min Rx .
The first, third and fifth equation in (32) are satisfied for P ΔΩ w = λw, ρ = 0. Thus, the optimum vector is an eigenvector of the Poyinting matrix P ΔΩ .The eigenvectors are usually scaled so that the norm of each is 1, satisfying the second equation in (32), i.e. w H w = c = 1.Lastly, the fourth equation in (32) are satisfied for δ max Tx − w (x, y, z Tx ) H P ΔTx w (x, y, z Tx ) = δ max Tx − λ ≥ 0 or w (x, y, z Rx ) H P ΔRx w (x, y, z Rx ) − δ min Rx = λ − δ min Rx ≥ 0. Thus, at each position on the propagation direction z, the optimum solution corresponds to the eigenvector(s) associated to δ min Rx ≤ λ or λ ≤ δ max Tx : The sufficient condition for (26) to have a constrained relative maximum, is that the determinant of the Hessian matrix formed by the second partial derivatives of (31), is positive definite.The application of this condition in this case leads to zero.Thus, the critical point is called a saddle point.The characteristic of a saddle point is that it corresponds to a relative minimum or maximum.By considering the solution of the problem ( 26) or ( 28) subject only to the second constraint in (27) or (29), it is possible to discriminate further.In this case, the optimum solution vector can be found using (33) and the determinant of the Hessian matrix is given by (34).Thus, the eigenvectors of P ΔΩ correspond to the maxima of (26).
The additional constraint can be included in the solution by examining the norm of the eigenvectors of P ΔΩ .It is noteworthy that the equality in the first constraint in ( 27) or ( 29) can be reached by evaluating the weight vectors (35) or (36) at z Tx or z Rx , respectively.
where λ max (z) and w max (z) correspond to the z-dependent maximum eigenvalues and eigenvectors of (33).The first term on the right-hand side of ( 35) and ( 36) normalizes the eigenvalues to the maximum or minimum desired powers.These magnitudes are modified so that w H Tx w Tx = δ max Tx /λ max Tx and w H Rx w Rx = δ min Rx /λ max Rx , which still satisfies the second equation in (32).In other words, optimal weight vectors in (35) and (36) ensure that the maximum and minimum powers in (33) coincide with the desired output power densities δ max Tx or δ min Rx at Tx and Rx, by modifying their norm.Thus, these vectors correspond to the maxima of (26).It should be noted that all the eigenvectors of (33) are solutions and here, only the maximum one was considered.Finally, the optimum fields at z Tx and z Rx can be found by using (1).They can be written in the matrix form:

Assumptions
In this subsection, the assumptions considered in this methodology are examined.First, to apply the modal theory, fields in these structures are assumed to vary as e −(jωt±jβ)z .Electromagnetic waves propagate along a uniform infinitely long structure and field solutions must satisfy the boundary conditions.Tunnel walls are supposed infinitely thick and hence fields are considered to be at negligible levels when they reach the end of the wall.As a result, the term in (7) corresponding to the radiated fields, is assumed to be zero.The geometry of a hollow dielectric waveguide with arbitrary cross-section is shown in Fig. 3. Tunnels are assumed to be straight.Corners and objects can be simulated as long as the tunnel cross section remains constant.
To calculate the mode parameters, a few points over the cross-section should be excited by considering the 2.5-Dimensional TLM approach for guiding structures.Then, the source term in ( 7) should be considered.These sources are employed to obtain all the possible modes that can be excited, but their description is not needed.Since this term depends on the z-field component of the source, it can be neglected by exciting the structure only by the transversal components.This means that the guide is excited on a slice of the tunnel and remains applicable for Transversal-Electric (TE), Transversal-Magnetic (TM) and hybrid mode excitation.It is due to the fact that at least the value of one transversal component is non-zero.
The modal expansion (1), allow us to obtain all the modes that propagate in the structure.However, assuming that the objective is to design the transmitter antenna, in a more realistic scenario, only the modes excited by this antenna (source) and captured by the receiver should be considered.This can be accomplished by exciting the structure in the transmitter region, as shown in Fig. 3. Ae Tx is the area where the power would be radiated if the source is placed somewhere in Tx.For an unloaded tunnel it corresponds to the whole tunnel cross-section.The excitation is simulated by an electric field probe concentrated in a differential of area ΔS in Tx.The power density of the i-th element has to be calculated over ΔS to apply the restrictions given by (27).
Finally, the orthogonality property for the modal functions in (1) is not mandatory.Its existence guarantees that the total power is the addition of individual modal power.However, as explained before, in general for lossy hollow-dielectric waveguides, such as tunnels, modes are non-orthogonal and the weight coefficients depend on z.Thus, these cross-power interactions have to be considered.

Algorithm flow chart
The antenna synthesis problem can be stated as the inverse of the analysis problem, i.e., given a set of design specifications, such as required fields, type of excitation and positioning, determine an optimum antenna.Optimization techniques can be employed to solve a constrained problem and fulfill the design requirements.The flow chart in Fig. 4 illustrates the process for finding the optimum field specifications and positioning in tunnel environments.
It is divided into five steps.The first one, concerns the definition of the inputs to simulate the tunnel.The characteristic parameters of the modes α, β, E (x, y) and H (x, y), are calculated at this point.In the second step, the modes are discriminated by their power carried through the tunnel, so that only modes among those with the highest power, are considered.Then, in the third step, an optimization procedure on the weight coefficients is carried out to obtain an optimal vector that maximizes the transmitted power by the modes.In the fourth step, the optimum weight coefficients to maximize the mean power along the propagation direction are obtained.Finally, the optimum fields specifications, type of excitation and source positioning are determined by using the optimum weights.

Description
Finding antenna design specifications is a synthesis problem.It is of the utmost importance for correct antenna integration in over-sized guiding structures considered here.The fundamental problem facing antenna engineers is to improve the power reception in a certain desired region by using all possible parameters, so that the communication system meets the required specifications.An optimization algorithm is proposed and described below.It is divided into five steps.
The starting point is to solve this electromagnetic problem with the 2.5-Dimensional TLM for guiding structures.It consists in defining the geometry of the structure and constitutive parameters of the involved materials.The solution region is divided into a number of elementary nodes and their maximum size is defined by the minimum signal wavelength (or maximum frequency).Its value must be much smaller than the minimum wavelength to avoid the so-called numerical dispersion errors; a rule of thumb is to consider this value as one-tenth of the wavelength at the maximum frequency.The analysis is restricted to regions where the transmitter (Tx) and receiver (Rx) are usually located, and only the fields in these regions should be sampled.The sampling step is defined by the operation frequency and it should be chosen small enough to guarantee that fields are almost constant, for instance by considering a value of one-tenth of the wavelength at the operation frequency.Then, by using the previous entries in the 2.5-D TLM for guiding structures, the characteristic parameters for the modes can be calculated.
It is worth to remember that this study is meant to define the main factors affecting radio communication in tunnel environments, and to establish some criteria to diminish the undesired effects.Multi-modal propagation is one of the most detrimental factors degrading communication in tunnel environments.The problem of mode excitation has already been studied by means of modal theory [15].High-order modes are responsible for rapid fluctuations in power, and their effects should be mitigated as much as possible.Moreover, a high number of modes leads to a considerable calculation time.Thus, in the second step of the proposed methodology, the number of considered modes is reduced.The strategy is straightforward: Assuming that the transmitter Tx is placed at z = 0, modes are classified by their power contribution in Rx located at z along the longitudinal distance.Since the calculation domain has been discretized, the (l, k) element of the matrix of crossed powers with dimensions N × N and given by (30), is approximated by using: where ΔΩ is the approximation of the differential element of area in Ω and (x, y) correspond to the coordinates in the region for the receiver where the power has to be optimized.A singular-value decomposition (SVD) of (38), P Ω = UΣV, is done at each point along the z-direction.U and V are unitary matrices composed of eigenvectors of P H Ω P Ω , and Σ is a diagonal matrix containing the singular values of P Ω .The choice of the number of selected modes M is made at this stage.We look at the ratio of the various singular values with respect to the largest one.Consider the ratio for the i-th singular value ς i , given by: A threshold value for the parameter u, indicating the number of significant digits that pertains to the i-th mode, is established.The first M modes with singular values along z above this threshold, where M ≤ N, are considered.Hence, those have some significant contribution to the response at Ω, while others should be neglected.This technique is usually employed to best estimate the order of significant poles in a function that can be expressed as a sum of complex exponentials, as in (1).Note that the procedure accounts for the interaction among all modes.A similar parameter is used in MIMO systems for confined environments as a figure of merit to evaluate its performances [56].
In the third step, the set of optimum weight coefficients to excite the modes is defined by using the optimization procedure presented in subsection 3.2.4.In practice, the maximum power Electromagnetic Wave Propagation Modeling for Finding Antenna Specifi cations and Positions in Tunnels of Arbitrary Cross-Section density that can be transmitted by Tx and/or the minimum power density that can be detected at Rx, denoted as δ max Tx and δ min Rx , respectively, are defined by the industrial specifications.The optimum weight coefficients for all the distances 0 ≤ z ≤ Z max are calculated by using (33).The approximation is done in a differential element of area in Tx or Rx where the excitation is concentrated and (x, y) correspond to the coordinates where it is located.Equation ( 33) is evaluated for all the points in Ω.A set of the optimum weight coefficients for each point in the (x, y)-plane, are obtained.The upper limit for the weight amplitudes can be found by scaling the power density using (35) or (36), so that the extreme values in Tx or Rx are satisfied.For Ω = Tx, the weight amplitudes in (35) have to be scaled with respect to the maximum eigenvalue in Tx, and for Ω = Rx, the weights in (36) have to be scaled with respect to the maximum eigenvalue in Rx.It should be noticed that, according to the previous step, the number of contributing modes is M and only the most contributing eigenvalues and eigenvectors are considered, so that the amplitudes of the remaining eigenvectors can be adjusted such that Pmin Rx ≤ w H Rx P ΔRx w Rx in (27) or w H Tx P ΔTx w Tx ≤ Pmax Tx for (29).Hitherto, the optimum weight coefficients were calculated at each point in Ω, obtaining a set of values that depends on the (x, y, z) position.The z-dependence of these coefficients is due to the non-orthogonality between modes.The (x, y)-dependence indicates which coefficients should be used if the structure is excited in the (x, y)-plane.Thus, the fourth step consists in finding the optimum position in the (x, y)-plane to calculate these coefficients.In doing so, the mean of the total power is calculated at each point (x, y) in Ω for the distances z along the axis of the tunnel.This procedure is repeated for all the points in Ω, obtaining a matrix containing the means for all the points (40).Next, the criterion is similar to that already explained.It consists in finding the positions where the mean of the power in Ω versus distance is maximum.P(x, y) = Lastly, in the fifth step, the field specifications to excite the tunnel and the best transmitter location are determined.The total fields are calculated by using the optimum weight coefficients obtained in the third step at the position where the mean ( 40) is maximum and by using the mode characteristic parameters obtained in the first step.Expression (37) gives the optimum fields at any cross section z ≥ 0. This expression can also be employed to define the privileged polarization by evaluating the total field components separately, and observing the predominant one.Finally, points where the matrix ( 40) is maximized belong to the best excitation points for the M considered modes according to subsection 3.2.3 and ,thus, the best source location.

Validation
The first three stages are essentially the core of the methodology presented in the previous subsection: The modal approach, the non-orthogonality between modes and an optimization technique were employed to determine the optimum weight coefficients.To validate these steps, the study of a simple theoretical "reference solution" for a canonical geometry was considered.A metallic-rectangular waveguide was studied in this subsection to help understand the treatment of a realistic scenario.
A dielectric-loaded WR-90 rectangular waveguide was considered as an exact field solution exists.The loading lossy material has a conductivity of σ = 0.01Sm −1 and a relative dielectric constant equals to the unity.The waveguide has a width w = 2.286 cm and height h = 4 cm and the operation frequency was chosen at f = 14 GHz.The objective consists in maximizing the power conducted by the first two modes in a region on the left side over the cross-section of the guide.Figure 5 illustrates the geometry of the waveguide and the cutoff frequencies of the first three mode.For the sake of simplicity, the gray region represent the possible locations for the transmitter and receiver, and consequently, where the power has to be maximized.As it was explained at the end of section 3.2.2, it should be pointed out that because only a part of the cross-section of the waveguide is considered, the orthogonality property is destroyed, so that this example serves as validation case of a real tunnel environment.Regarding the first step of the proposed methodology, the exact expression for the fields can be found in [51] for this case.So, the 2.5-D TLM is not required to determine the mode characteristic parameters.For the second step, it can be demonstrated that the matrix of crossed powers is given by: : ω is the angular frequency, μ = μ 0 is the permeability inside the waveguide, w and h are the waveguide dimensions and γ 10 = α 10 + jβ 10 and γ 20 = α 20 + jβ 20 are the propagation constant of the first and the second modes, respectively; α 10 , α 20 and β 10 , β 20 are their corresponding attenuation and phase constants.Lastly, the P l,k -element in (41) can be calculated by using (30).For reasons that will become clear later, one can consider the first two singular values of P. Additionally, it is worth noting that the term P 3,3 is zero for frequencies below 14.8 GHz.In that case, the matrix (41) has only two non-zero singular values.Thus, the number of significant modes is set to M = 2 and only the terms P 1,1 , P 1,2 , P 1,1 and P 2,2 in (41), are considered.
The third step is to obtain the solution for the constrained problem (26).The maximum power density that can be transmitted by Tx or the minimum detectable power density in Rx location, may be specified.Suppose that the maximum and minimum desired power densities are δ max Tx = 1 Wm −2 and δ min Rx = 10 μWm −2 , respectively.Equations ( 35) and ( 36) are employed to calculate the optimum weight coefficients and the results are shown in Fig. 6.These coefficients fulfill the power constraints at the maximum and minimum distances.Any choice for the weight amplitudes between the solutions for the maximum and minimum powers could be considered.By observing these figures, one can intuitively suppose that they can be approximated by a sum of complex exponentials due to their dependence with the modes.The Matrix Pencil method is a very popular technique to estimate the parameters of a sum of complex exponentials [25].By using this method, the weight coefficients can be determined and are listed in (43) for the case of δ max Tx .This result demonstrates that some coupling between non-orthogonal modes exists.Similar results can be obtained for δ min Rx .
Figure 7 illustrates the amplitudes of the inner terms of the matrix P as a function of the longitudinal distance.The fluctuations observed in the crossed terms coincide with those of the weight amplitudes, confirming the mutual dependence between modes, as shown in (43).Finally, the black curves in Fig. 7 illustrate the constrained solutions for the maximum power at z = 0 and minimum one at z = z zmax = 50 cm.It is interesting to observe that the influence of the crossed terms P 12 and P 12 is higher at closer distances to the source, i.e. at z = 0.This is due to the strong coupling close to Tx, as can be seen in Fig. 6.
To verify this result, the intersection of the curves for the total power (26) and the constraints (27) has to be found.Different weight amplitudes w were considered in the power function at the maximum and minimum distances and the points for which the constraints were achieved were plotted.It is important to note that, thanks to the fact that only two modes were considered, a 3D plot of the power versus the amplitudes of both modes can be obtained, as shown in Fig. 8.This result confirms that the solution given by equations ( 35) and ( 36) maximize ( 27) the power at Tx and Rx and their values correspond to Pmax Tx and Pmin Rx .

Results for a rectangular tunnel
The determination of optimum-field specifications and antenna positioning in tunnel environments following the basic principles and procedure outlined in the previous subsections is now straightforward.Multi-modal propagation is experimented in these scenarios and the correct excitation of these modes determines the efficiency of a given source.Tunnels are made with dielectric walls and contain other materials (trains, cars, objects etc.).The use of a full-wave approach for mode parameters calculation now becomes evident due to the necessity of modeling the tunnels walls, near-fields and arbitrary cross sections.The calculation time is the main disadvantage of full-wave methods compared with commonly used ray tracing techniques.The reduced 2.5-D TLM node and the SIBC concept avoid large region meshing and hence, reduce the required computer cost.In this section, the source optimization procedure is applied in the case of a real tunnel environment.
First of all, the wave propagation in a rectangular tunnel with transverse dimension w = 7m × h = 5m was considered.The tunnel walls were modeled by a lossy dielectric with r = 12 and σ = 0.02Sm −1 .The central operating frequency and the maximum frequency to simulate for the 2.5 D TLM for guiding structures were established at 2.4 GHz and 3 GHz, respectively.Thus, for 2.5-D TLM simulation the mesh size was set to Δl = 0.009 m.The regions for the Tx and Rx were discretized by taking samples at each λ/10 at the operation frequency, i.e., ΔS = c 0 / f op .The tunnel transverse section and the zones in which the transmitter (Tx) and receiver (Rx) can be located, are shown in (Fig. 9).This tunnel was simulated using the 2.5-D TLM for guiding structures and the time-domain SIBC implementation.The transversal field components of the electric field were excited in Tx in order to calculate the field specifications for the transmitter.The attenuation and phase constants, and the field profiles of the modes were determined at the operation frequency.
Figure 10 shows the attenuation and phase constant up to 3 GHz and Fig. 11 shows the field profiles of the first two modes.Slight discrepancies are observed for the mode profiles and the phase constant β.The discrepancies in the attenuation constant α for low frequencies are explained by two facts: First, the comparison was made with a commonly used theory, namely Marcatilli's theory, which neglects the fields on the corners.Secondly, because the calculation of the attenuation constant α is derived from β.One can show that errors on β can produce large calculation error on α.However at the operation frequency, this error can be considered to be negligible.
In the next step, the matrix of crossed powers P Ω of dimensions N × N was calculated for the first N =32 modes and a singular value decomposition of P Ω was carried out.So, the region of the receiver was sampled at one-tenth of the wavelength of the operation frequency and the (l, k)-element of (30) was approximated by the Riemmann sum (38).The choice of the M best modes to carry the power was made at this stage.The normalized amplitudes of the singular values were calculated and only those accurate up to 4 significant digits were considered.The amplitudes were obtained for the first 1,000 m in the tunnel, which constitute the maximum distance between Tx and Rx.The singular values for P ΔRx are shown in Fig. 12, the value M was set to 11.The remaining singular values for which the ratio in equation ( 39) is below 10 −u = 10 −4 are neglected and are not be used in the calculation of the reduced matrix P of dimension M × M. One can note that the propagation in the tunnel is mainly dominated by the first two modes.[31,33] and the procedure in this paper [3].Next, the matrix P ΔΩ of dimensions M × M was calculated at each point on the cross section of the tunnel and all distances along the z-direction.The optimum weight coefficients that 281 Electromagnetic Wave Propagation Modeling for Finding Antenna Specifi cations and Positions in Tunnels of Arbitrary Cross-Section maximize the cost functions ( 26) and ( 28) subject to the constrains ( 27) and ( 29), respectively, were calculated for all points.The power density specifications were defined to be 0 dB Wm −2 at z Tx = 0 and -40 dBWm −2 at z Rx = 1, 000 m.This can be accomplished by scaling the weight coefficients, as explained in subsection 3.2.4.Then, for each pair of coordinates (x, y), the field evolution was calculated along z and the mean power was obtained by using (40).Figure 13 shows the mean power for different positions in Tx and Rx.It is interesting to note that the mean power is maximum for positions close to the center for both Tx and Rx.The evolution of the normalized weight coefficients as a function of the longitudinal distance for points where the mean power is maximized in Tx and Rx, were computed (see Fig. 14).To satisfy the condition at z Tx = 0 and z Rx = 1, 000 m, the levels should be increased or decreased by 26.78 dB or -35.53 dB, respectively.Using these coefficients, the optimum power in Tx and Rx were calculated.The best and worst case for the mean power were also computed and compared.As expected, the required specifications were satisfied, as illustrated in Fig. 15.Moreover, from this figure, one can observe that, by correctly locating the source, improvements in the received power of the order of the order of 10 dB, can be obtained.
Finally, the optimum electric fields were obtained by using the optimum weight coefficients in (37).By way of illustration, the three field components of the electric field that satisfy the power density specifications in Tx were calculated.They were computed for the overall tunnel cross section at z = 0, as shown in Fig. 16.Intuitively, one could expect that optimum fields are most influenced by the lowest attenuated modes.This result can be confirmed by observing that the points where the x-component of the electric field in Fig. 16 is minimum at corners, and the z-component is minimum at half of the width as it occurs for the E x 11 mode.This observation can also be confirmed as the x-component is maximum at one-third, as well as one-quarter of the width as for both E x 21 and E x 31 modes.It is important to clarify the fact that the weight coefficients vary along the longitudinal distance is due to the non-orthogonality between modes.Note that this variation does not imply that the tunnel must be excited with different weights at each point.Therefore, coefficients should be given only at z = 0. Finally, the optimum positions to locate the transmitter are those where the mean power is maximum, i.e., for the points close to the white regions in Fig. 13.

Analysis and discussion
A procedure for determining the optimum antenna field specification and location in guiding structures was developed and presented in details.The understanding and modeling of radio-wave propagation through the 2.5 D TLM for guiding structures, were essential steps toward deriving this approach.Radio wave propagation inside tunnels was traditionally modeled using ray tracing and ray launching techniques based on geometrical optics (GO) [13,15,19,34,41,57].However, the analysis can be rather complicated at long ranges due to the fast growing number of contributing rays and breaks down completely at and near caustics [32].Moreover, the main disadvantage consist in that near-field effects cannot be modeled, and as a result, antenna performances cannot be guaranteed by using these methods.This explains why full-wave techniques are preferred.In the past, full-wave techniques, such us Finite Element Method (FEM), were employed to analyze the wave propagation in these structures [8].However, these solutions were very limited due to the complexity of the resulting Electromagnetic problem and computation burden.The fundamental strength of the proposed approach is the simplification of the calculation to deal with this electromagnetic problem.Nowadays other recent attempts to deal with these structures are emerging and they are still under development [52,55].
The modal representation of the fields and the non-orthogonality between modes were the main assumptions involved in developing this methodology.A metallic-rectangular waveguide was considered to validate this procedure.The validation was done by comparing the power density obtained with the optimum weight coefficients obtained with our procedure and comparing with those obtained for different combinations of the weight amplitudes.The results confirm that this approach allows us to obtain optimal weight coefficients to excite the modes.Lastly, the study of a rectangular tunnel was considered.The characteristic parameters of the modes were computed and the most contributing modes in Ω were determined.The resulting modes are coupled to each other and an isolated excitation of them is not possible, as observed in Fig. 14.This fact can also be seen in Fig. 6 for the validation case.In Fig. 14, one can observe that the amplitude of the first mode is almost constant.Thus it can be considered quasi-orthogonal due to its independence with the other modes.Large fluctuations are observed for higher order modes, which is in agreement with the theory.Analysis of the elements of the weight vector w with the corresponding mode functions Ψ E ( r) and Ψ H ( r) in (1) reveals that, for a correct excitation of the 23 and E x 13 must be considered.Regarding this figure, it can be noted that the wave-propagation is predominantly dominated by the first three modes.This result is in agreement with the works in [43] for a rectangular tunnel.The weighted sum of these coupled modes constitute the optimum fields to excite the structure.On the one hand, from these results of the optimum fields, one can conclude that, for this case, an excitation by the x-component of the electric field is preferred.In [15], it was also verified that, for rectangular tunnels, the vertically polarized modes are much more affected than the horizontally polarized ones.In practice, these fields are difficult to generate and they are usually approximated with dipoles or loops.Consequently, a deeper study to link these solutions with typical antennas are currently under investigation.On the other hand, from the result for the optimum positions, one can conclude that there is a significant influence of the first mode and the transmitting source should be located as close as possible to the center.
Generalization of these concepts to analyze arbitrary-shaped tunnels may be possible by using the procedure outlined in this chapter.Full characterization of propagation in different railway scenarios (several tunnel cross-section and dimensions, loaded and unloaded tunnels) may be analyzed with the 2.5-Dimensional TLM for guiding structures.Useful information, such as field distribution, mode propagation characteristics (attenuation and phase constants) and depolarization effects for various scenarios, can be obtained.The treatment of a more realistic scenario where different cross sections are present, can be studied by establishing a common region to calculate the source specifications.This region should correspond to the biggest possible common area for all the cases considered.A well-adapted antenna for these environments should match as well as possible all results obtained for each case.Loaded or unloaded tunnels with different cross sections may be then analyzed with the proposed procedure.

Conclusions
A methodology for finding optimum antenna field specifications and positioning in tunnel environments has been presented and discussed.A TLM approach suitable for the analysis of radio-wave propagation in tunnels was employed and briefly presented.By this approach, the phase constant, attenuation constant and field profiles can be determined.The concept of Surface Impedance Boundary Condition (SIBC) was additionally employed to reduce the TLM computational domain.Combined techniques can be applied to electrically large uniform guiding structures such as lossy tunnels with arbitrary cross-section.
The optimality criterion of the proposed methodology is the maximization of the power density in a given region, as required in industrial applications.The optimum field specifications and positioning in a guiding structure are obtained by means of the modal fields and the optimum weight coefficients satisfying the power density requirements.A metallic waveguide was used to validate this approach and results for a rectangular tunnel with lossy dielectric walls were obtained.Results in terms of field distribution and positioning allow us to analyze the relevant modes and the coupling between them, as well as the privileged polarizations and enhancements in the performances by using a proper excitation.
Concerning future work, results of this approach were presented for a Single Input and Single Output (SISO) system.However, the procedure can be extended for Multiple Input Multiple Output (MIMO) systems to further improve communication system in tunnels.

Figure 3 .
Figure 3. Geometry of an uniform section infinitely long hollow waveguide.

Figure 4 .
Figure 4. Flowchart of the process for determining the field specifications and positioning in tunnel environments.

Figure 5 .
Figure 5. a) Geometry of the rectangular waveguide and region where the power has to be maximized and b) Cutoff frequencies of the first three modes.

Figure 6 .w
Figure 6.Optimum weight coefficients satisfying the maximum and the minimum required power densities.

Figure 7 .
Figure 7. Optimum power density conducted by the modes and terms of the matrix P versus distance.

Figure 8 .
Figure 8. Optimum weight coefficients satisfying the maximum and the minimum required power densities.

Figure 9 .
Figure 9. Schema of the rectangular tunnel with w = 7m and b = 5m.Rx denotes the position for the receiver (fixed) and Tx position of the receiver (to be optimized).

Figure 10 .
Figure 10.a) Attenuation constant in dB/km, and b) Dispersion curves for the E y 11 and E x11 modes, calculated in the frequency range from 0.4 GHz to 3 GHz with Marcatilli's theory[31,33] and the procedure in this paper[3].

Figure 11 .
Figure 11.Field configuration in dBV/m of the E y 11 and E x 11 modes: a) Calculated with Marcatilli's theory and b) Calculated with the procedure in this paper [3].

Figure 12 .
Figure 12.Normalized amplitudes of the first eleven singular values versus distance.

Figure 13 .
Figure 13.Mean of the power at z = 0 for Tx and z = z max for Rx.

Figure 14 .
Figure 14.Normalized optimum weight coefficients for the non-zero valued modes in dB versus distance to satisfy the required specifications in Tx and Rx.

Figure 15 .
Figure 15.Power versus distance for the best and worst points in Tx and Rx.

Figure 16 .
Figure 16.Optimum electric fields at z = 0 for optimum excitation of the rectangular tunnel.