%~Mouliné par MaN_auto v.0.27.3 2024-01-08 13:13:07
\documentclass[CRPHYS, Unicode, XML, screen, Thematic, published]{cedram}




\usepackage{rotating}



%\newcommand{\bbR}{\mathbb{R}}
%\newcommand{\bbC}{\mathbb{C}}
%\newcommand{\bbZ}{\mathbb{Z}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\RubricEF{URSI-France 2023 Workshop}{Journées URSI-France 2023}
%\TopicEF{Section information missing}{Intervention en colloque}


\title{A Simulation Method Suited for the Whole French Territory Electromagnetic Waves Exposure}

\alttitle{Une méthode de simulation adaptée à l'exposition aux ondes électromagnétiques de l'ensemble du territoire français}

\author[N. Noé]{\firstname{Nicolas} \lastname{Noé}\IsCorresp}
\address{Centre Scientifique et Technique du Bâtiment, Nantes}
\email[N. Noé]{nicolas.noe@cstb.fr}

\author[L. Sefsouf]{\firstname{Lydia} \lastname{Sefsouf}}
\address{Agence nationale des fréquences, Maisons-Alfort}
\email[L. Sefsouf]{lydia.sefsouf@anfr.fr}

\author[J.-B. Dufour]{\firstname{Jean-Benoit} \lastname{Dufour}}
\address{Geomod, Saint-Didier-au-Mont-d'Or}
\email[J.-B. Dufour]{jean-benoit.dufour@geomod.fr}

\author[S. Carré]{\firstname{Samuel} \lastname{Carré}}
\addressSameAs{1}{Centre Scientifique et Technique du Bâtiment, Nantes}

\author[E. Conil]{\firstname{Emmanuelle} \lastname{Conil}}
\addressSameAs{2}{Agence nationale des fréquences, Maisons-Alfort}

\author[N. Bounoua]{\firstname{Nabila} \lastname{Bounoua}}
\addressSameAs{3}{Geomod, Saint-Didier-au-Mont-d'Or}

\author[J.-B. Agnani]{\firstname{Jean-Benoit} \lastname{Agnani}}
\addressSameAs{2}{Agence nationale des fréquences, Maisons-Alfort}

%\shortrun

\keywords{\kwd{EMF exposure} \kwd{mobile telephony} \kwd{numerical simulation}}

\altkeywords{\kwd{Exposition aux champs électromagnétiques}
\kwd{Téléphonie mobile}
\kwd{Simulation numérique}}

\begin{abstract}
As part of the process for monitoring public exposure to electromagnetic waves, the national frequency Agency is carrying out, in collaboration with the Ministry of Ecological Transition and Territorial Cohesion, the Scientific and Technical Center for Building (CSTB) and Geomod, a numerical modeling of the electromagnetic wave exposure levels emitted by mobile telephony base stations on a national scale. This paper presents a dedicated simulation method for the numerical modeling of the whole French territory's exposure to EMF. This method accounts for EMF exposure everywhere (outdoors and inside buildings), while performing fast enough to fulfill operational constraints. The simulation method relies on a ray-based 2.5D and 3D mixed approach that takes into account computation areas depending on the radiated power pattern of the antennas. The method is yet deployed on pilot areas before a full deployment.
\end{abstract}
\editornote{This article follows the URSI-France workshop held on 21 and 22 March 2023 at Paris-Saclay.}

\begin{altabstract}
Dans le cadre de son plan de surveillance de l'exposition du public aux ondes, l'Agence nationale des fréquences réalise, en collaboration avec le ministère de la Transition écologique et de la Cohésion des territoires, le Centre scientifique et technique du bâtiment et Geomod, une modélisation numérique des niveaux d'exposition aux ondes électromagnétiques émises par les antennes relais de téléphonie mobile à l'échelle du territoire national. Cette modélisation nécessite la mise en \oe{}uvre d'une méthode de calcul spécifique. La méthode présentée ici permet de restituer un niveau d'exposition en tout lieu (à l'extérieur et à l'intérieur des bâtiments) tout en conservant des temps de calcul compatibles avec les contraintes opérationnelles. La méthode est basée sur une approche à rayons combinée 2,5D et 3D et sur la prise en compte de zones de calcul dépendant des caractéristiques d'émission des antennes. La méthode est aujourd'hui mise en \oe{}uvre sur des zones d'expérimentation avant son déploiement complet.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
\dateposted{2024-03-05}
\begin{document}
\maketitle

\section{Context and main goal}

Engagement 8 of the fourth health and environment national plan~\cite[\emph{plan national santé environnement, PNSE4}]{pnse4} aims at containing exposure to electromagnetic waves. It focuses mainly on the knowledge of the electric field generated by mobile telephony base station antennas. Within this plan, the French Frequency Agency (\emph{Agence nationale des fréquences, ANFR}) has been tasked with the Ministry of Ecological Transition (\emph{Ministère de la Transition \'Ecologique et de la Cohésion des territoires}), the Scientific and Technical Center for Building (\emph{Centre Scientifique et Technique du Bâtiment, CSTB}) and the company Geomod, with creating a numerical modelling of electromagnetic wave exposure levels induced by mobile telephony base station antennas at a national scale.

Figure~\ref{Fig_carte} was generated using a beta release of the web service that exposes those simulation results. The left part of the figure is a color map of the exposure level~1.5m above the ground. On the right part of the figure, each building outline (in the horizontal plane) is colored using the maximum indoor exposure level on the frontage, along a vertical line.

One of the specificities of this numerical modelling is the extensive computation area to consider. As a matter of fact, up to this day, most exposure simulations were only performed in the vicinity of the antennas. These simulations are usually dedicated to assessing compliance with exposure standards or to enforcing local regulations. For instance in France, mandatory records based on simulation have to be submitted when a new antenna is deployed~\cite[\emph{dossiers d'information mairie}]{dim}. Similarly, environmental permits in the Bruxelles Capitale region in Belgium~\cite{brux} also require simulation results. These regulations either (France) conform with the\linebreak European Commission recommendation on the limitation of exposure of the general public to electromagnetic fields (0~Hz to 300~GHz), based on the guidelines published by the International Commission on Non Ionising Radiation Protection~\cite[ICNIRP]{icnirp} or can be more restrictive (Bruxelles Capitale region). Nevertheless they are always limited to the vicinity of the antennas. As far as coverage maps are concerned, they span very large areas, but with a low spatial resolution. The aim pursued here is to estimate exposure levels anywhere in space, even far from antennas, along with a precise knowledge of the exposure levels close to the antennas.

Such exposure maps could then be used for statistical analysis (to compare exposure levels between different areas for instance) and to monitor temporal evolution of exposure level (due to the ever changing radioelectric sources) or even to test future scenarios before they are implemented.



\section{Implementation}

Achieving the aforementioned goal implies being able to produce this whole France map in a reasonable time and to be able to update it on a regular basis. Figure~\ref{Fig_schema1a} gives an overview of the full simulation process, starting from databases (stage~0) to the final exposure maps (stages~6 and~7).

An adapted computation method is needed, that handles each base station (a batch of antennas) independently, as illustrated on Figure~\ref{Fig_schema1b}, with stages~1 to~5. Results on batches of antennas are merged at local scale, using regular tiles (stage~6) to deliver a cumulated exposure level everywhere in France (stage~7). This cumulated level might be identical to zero, meaning it should be below a measurement threshold (typically 0.05~V/m).

\begin{figure}[!hp]
\includegraphics[width=0.48\textwidth]{carte.png}
\includegraphics[width=0.45\textwidth]{cartev.png}
\caption{Example of exposure maps : electric field 1.5m above the ground (left) and on the contours of buildings (right).}\label{Fig_carte}
\end{figure}

\begin{figure}[!hp]
\includegraphics[width=0.9\textwidth,page=1,clip,trim=0cm 0.6cm 0cm 0.3cm]{schema1}
\vspace*{-5pt}
\caption{Overview of the full process : from databases to web services for end users.}\label{Fig_schema1a}
\end{figure}

\begin{figure}[h]
\includegraphics[width=0.9\textwidth,page=2,clip,trim=0cm 0.6cm 0cm 0.3cm]{schema1}
\caption{Overview of the simulation process for a single antenna batch.}\label{Fig_schema1b}
\end{figure}

\subsection{Input data}

\subsubsection{Antennas}

Antennas descriptions come from ANFR database STATIONS (an internal database filled with information by telephone operators). From the computation side, an antenna is associated to a 3D position and an orientation (aiming given as two angles : azimuth in the horizontal plane and mechanical tilt in the vertical plane). An antenna corresponds to a physical ``box'', and each antenna is made of one or more transmitters (one transmitter per technology - 2G, 3G, 4G, 5G and frequency band). Each transmitter is associated with an EIRP (Effective Isotropic Radiated Power) and a far field radiation pattern (including electrical tilt). This radiation pattern (from manufacturer databases) is usually described by two cut-planes of the gain (one in the horizontal plane and one in the vertical plane). A full 3D radiation pattern is recomposed from these two 2D patterns (see Figure~\ref{Fig_MSI}) in order to be able to estimate gain in any 3D direction. \cite{msi3d} details a recomposition method and references others but for base station antennas, it appears that the method used by radio planning software~\cite{atoll} is the most suited.

\begin{figure}[h]
\includegraphics[width=0.35\textwidth,clip,trim=0cm 0.55cm 0cm 0.4cm]{msilin2.png}
\includegraphics[width=0.35\textwidth,clip,trim=0cm 0.55cm 0cm 0.4cm]{msilin3.png}
\caption{Radiation pattern of an antenna (linear) : on the left, the horizontal and the vertical cut-planes (input data), on the right the recomposed 3D pattern used for simulation.}\label{Fig_MSI}
\end{figure}

The STATIONS database does not give a precise location of antennas. There is only a single approximate location of each base station (even if the base station is made of antennas with different azimuths at different corners of a building). This single location (that could be interpreted as the barycenter of all antennas) has also uncertainty due to numerical rounding of latitude and longitude in degree / minute / second natural integers), leading to a potential 20~meters error. Hence, the antennas have to be relocated using a dedicated algorithm. Nevertheless, the type of antenna support (mast or building) is known. For the latter support, the algorithm relocates antennas to the center of the building, while not moving them more than 30~meters from their original position.

\begin{figure}[h]
\includegraphics[width=0.9\textwidth,page=3,clip,trim=0cm 5.4cm 0cm 2.05cm]{schema1}
\caption{Relocation algorithm.}\label{Fig_recalage}
\end{figure}

The relocation algorithm (stage~3 on Figure~\ref{Fig_schema1b}) is illustrated on Figure~\ref{Fig_recalage}. While more realistic algorithms (dispatching different azimuths to different corners of buildings~\cite{amsterdam}) could be used, this choice is a mutual agreement of operators and regulator. The choice of the algorithm is unrelated to the simulation method and precise positions of antennas could be available in a near future or in local databases. Albeit, it must be noted that this introduces extra uncertainty in the simulation results to be discussed later.

\subsubsection{Geometrical model}

The geometrical model (see Figure~\ref{Fig_strasbourg}) mostly comes from the BDTOPO geometrical database from the French mapmaking agency ({\em Institut Géographique National}, IGN). In this case, it is made of 2.5D buildings, i.e. 3D polygons representing the outline of the roof of the building, and the building height. As a consequence, no precise shape of the roof is known. Nevertheless, models with more precise modelling (of the roof) can also be used as no hypothesis is made on the shape of the buildings (that can be true 3D buildings). The digital terrain model also comes from IGN databases (BDALTI).

\begin{figure}[h]
\includegraphics[width=0.6\textwidth]{Fig6.png}
\caption{Sample of the geometrical model used for simulation : BDTOPO IGN (on the left), 3D model provided by {\em Eurométropole de Strasbourg} with more detailed roofs (on the right).}\label{Fig_strasbourg}
\end{figure}

The shape of the roofs and the precise height of buildings is an important parameter in the forecoming simulations, as propagation paths will be computed in a 3D environment, and changing from LOS to NLOS (or reciprocally) due to some different masking effect can lead to important changes in exposure levels.

\subsection{Preprocessing}

\subsubsection{Computation areas}

All horizontally colocated antennas (within a five meter circle) are gathered in so-called batches (stage~0 on Figure~\ref{Fig_schema1a}). For each batch, a computation area surrounding the batch center in the horizontal plane is then precomputed (stage~1 on Figure~\ref{Fig_schema1b}). A minimal electric field threshold is chosen (the minimal level of exposure to be simulated later), and for each azimuth the maximum distance above which the electric field falls below the chosen threshold is computed using a heuristic method (extended Hata model~\cite{hatax}). The heuristic method is instantaneous and only uses the geographical area type (urban, rural, \dots), the power and the horizontal pattern of each transmitter in the batch. This approach avoids relying on a general maximum propagation distance and ensures an homogeneous behaviour of the simulation for antennas with heterogeneous powers and is illustrated in Figure~\ref{Fig_zoneH}.

\begin{figure}[h]
\includegraphics[width=0.3\textwidth]{zone.png}
\hspace{-1ex}
\fbox{\includegraphics[width=0.6\textwidth]{maillage.png}}
\hspace{-30ex}
\begin{tabular}{p{4cm}}\scriptsize computation points on the ground (magenta) and on buildings (red)\\\hphantom{x}
\end{tabular}
\caption{Precomputation of the computation area around an antenna base station made of three antennas with three azimuths (on the left) and computation points (on the right).}\label{Fig_zoneH}
\end{figure}

The overall process for a full French department (Bas-Rhin) is illustrated on Figure~\ref{Fig_zones}, for a threshold of -40~dBV/m (i.e 0.01~V/m). The geographical area type is translated from the population density $P$ (as inhabitants per square meter) in the city or city district at the antennas location using the formula in Table~\ref{eqdensite}.

\begin{table}[h]
\caption{From population density to environment type for extended Hata precomputation}\label{eqdensite}
\centering
\resizebox{1\textwidth}{!}{
\begin{tabular}{|l||l|l|l|l|l|} \hline
population density $P$& $P \le 25$& $25 < P \le 300$& $300 < P \le 1500$& $1500 < P \le 5000$& $P > 5000$\\ \hline environment type& rural& suburban& urban (small)& urban (median)& urban (large)\\ \hline
\end{tabular}
}
\end{table}

There are 5705~antennas with 21681~transmitters, that are divided into 334~computation batches. On one hand, it can be seen that the areas are larger in rural areas and that they are even parts of the territory where nothing will be computed (no exposure). On the other hand, in dense urban areas, most areas are overlapping, while they are smaller. This precomputation is mandatory for fast, exposure-level oriented simulations.

\begin{figure}[h]
\includegraphics[width=0.333\textwidth]{zones67a.png}
\includegraphics[width=0.41\textwidth]{zones67b.png}
\caption{Computation areas for the French Bas-Rhin department. On the left, the full department (within the magenta outline) and on the right Strasbourg city area only.}\label{Fig_zones}
\end{figure}

\subsubsection{Computation model}

After sizing the computation area, the ground and buildings are extracted from geographical databases. The ground is transformed into a TIN (triangular irregular network), with an interpolated normal to ensure continuous reflection. The buildings are made of 3D polygons. The computation area is then populated with computation points (receivers), on the ground and on the building frontages. Calculation point density varies with distance to the center of the area (see Figure~\ref{Fig_zoneH}). The horizontal step ranges from 2~meters to 50 meters at 1~km. Receivers on the front of buildings are created as vertical columns of multiple points, at different heights (with a vertical step three times smaller than the horizontal one). Far from the antennas, only one receiver per building frontage is created and for far and low individual houses, a single receiver is created atop the roof.

This varying receiver density aims at speeding up computations, but the method used is generic and does not rely on this adaptive receiver mesh.

\subsection{Geometrical computations}

A ray-tracing model is used for simulation, in order to find propagation paths connecting antennas and computation points. These propagation paths take into account reflection and diffraction by ground and buildings. Asymptotic methods (geometrical optics and uniform theory of diffraction) are well suited to the mobile telephony frequency range (700~MHz - 3.5~GHz, i.e. wavelengths ranging from 43~cm to 8~cm) and to the level of detail and precision of the geometrical models (approximately 1~m). Each computation is a two stages process, with a full 3D computation and a 2.5D computation.

\subsubsection{Geometrical computation : 3D stage}

A full 3D computation in done within a horizontal circular area centered on the antennas, with a typical 200~meter radius (see Figure~\ref{Fig_23D}, a subset of the model of Figure~\ref{Fig_zoneH}). This 3D computation takes into account reflection and diffraction by buildings (frontages and roofs) and ground, and is based on an adaptive 3D beam tracing algorithm~\cite{icare}. No assumption is made on the orientation of the reflecting surfaces and diffracting edges, that can be slanted or not. Computations are limited to two reflections (typically one on the buildings, one on the ground) and one diffraction (by building edges).

\begin{figure}[h]
\includegraphics[scale=0.5]{trajets2d.png}
\includegraphics[scale=0.5]{trajets3d.png}
\caption{On the left: examples of 2.5D computed paths (ground profile as a navy blue line, and 3D paths including ground reflection as light blue lines). On the right: examples of 3D computed paths.}\label{Fig_23D}
\end{figure}

The purpose of this 3D computation stage is to get more accurate results in the vicinity of the antennas without having to simplify the buildings' shape. Masking effect by different roof shapes (two-pitches, flat, \dots) and reflection not only on vertical planes are then enabled.

\subsubsection{Geometrical computation: 2.5D stage}

In addition to the 3D computation, a 2.5D computation is done on the whole area. A double loop spans all pairs of horizontal points (antennas and receivers at the same horizontal location). For each of these pairs, an algorithm (based on a 2D ray-tracer with all non-vertical building and ground edges) extracts the ground profile as seen from atop. Then, using antenna and receiver heights above the ground, an always existing 3D path is computed, using the convex envelop method~\cite{combeau2d}. This 3D path can either be direct or made of isolated and / or successive diffractions on the buildings and on the ground. Finally, ground reflections~(in a plane perpendicular to the profile plane) are added, with a single reflection between each pair of transmitter / isolated diffraction point / receiver as illustrated in Figure~\ref{Fig_23D}.

The purpose of this 2.5D computation stage is to have a result at every computation point as such a path always exists. This stage is also very fast as it's a 2D + 1D algorithm and thus can be done on larger areas. Furthermore, no building reflection is computed and ground reflections and diffractions are approximated compared to a dedicated 3D computation.

\subsubsection{Aggregation of 3D and 2.5D geometrical results}

A final geometrical stage aggregates 3D and 2.5D geometrical results. This is mandatory in order not to compute a similar 3D path in both stages, as this could happen (for direct paths, ground reflected path and others). This aggregation stage uses the path ``history'' (number and type of interaction - reflection, diffraction \dots, and on which objects). In the case of two overlapping paths, the priority is given to the one computed with the 3D algorithm, and the one computed with the 2.5D algorithm is discarded.

Finally, for each antenna - receiver pair, the path difference between the two main paths (with shortest lengths) is also stored. This information will be used later to interpolate the electric field at any location.

\subsection{Physical computation stage}

\subsubsection{Electric field computation}

The physical stage computes for each 3D path (whether it originates from the 2.5D or the 3D computation stage), the associated electric field, generated at the receiver location, for each transmitter associated with the path. The electric field is a complex 3D vector and is computed using the radiation pattern of the transmitter, its polarization, and the reflection and diffraction coefficients along the path. The underlying physical model can handle successive diffractions, whether they are isolated diffractions on edges or creeping waves (by crawling over the top of a building or over a hill), without limitation of the number of diffractions.

While from a theoretical point of view the electromagnetic properties could be different for each building, no such information exists. As a consequence, a global approach is used with unique building and ground materials. The materials are defined by relative permittivity $\epsilon_r$ (-) and conductivity $\sigma$ (S/m). The buildings' material is a concrete like one ($\epsilon_r = 6$ and $\sigma = 0.03$, leading to an average reflection loss of $\approx$ 6.7~dB in the considered frequency range). The ground material is a very dry ground, taken from an ITU recommendation~\cite{itumat}, leading to an average reflection loss of $\approx$ 9.2~dB.

As far as reflection coefficients are concerned, most construction materials (concrete, brick, wood, \dots) have not-so-different reflection loss ($\approx$ 6~dB), so the influence of the materials (aside from metallic ones) is limited.

On the contrary there is as large diversity of the soil type in different areas or even in the same area. Nevertheless it is very likely that no ground reflected path exists without an equivalent path minus this ground reflection, with a higher level, and the diffraction model used does not take into account ground material. Computations without any ground reflections also exhibited that having no ground material downgrades results (compared to measurements).

The electric field is computed 1.5~meters above the ground but also inside buildings, just after the first wall. This is done by using a heuristic transmission model (from~\cite[COST 231]{COST}) that handles incident angle of incoming waves on the outer wall. Once again, a unique transmission loss factor is used, a 3,8~dB normal (typical of a double glazing).

The band associated with each transmitter is sampled in 16~frequencies, each of them carrying a fraction of the power. The electric field is computed at each of these frequencies before a final quadratic summation to get a scalar electric field generated by the transmitter at the receiver location. This scalar value, along with the path difference, is stored in the output files.

\subsubsection{Assessment of first Fresnel ellipsoid clearance}

In the case of a direct path, the first Fresnel ellipsoid is built around it and its clearance $p \in [0;1]$ is computed using ray-tracing, with discretized curved rays as illustrated in Figure~\ref{Fig_ell}. For a knife-edge like obstacle~\cite{itudeygout} this clearance can be related to the $v$ parameter with $v = 2 \sqrt{|\delta| / \lambda} = (2 p - 1) / \sqrt{2}$, where $\delta$ is the marching difference. An extra path loss $L = 12 p$~dB ($p = 0$ for a total clearance, $p \rightarrow 1$ for a full obstruction) is then applied to the electric field of this direct path as an approximation of the extra path loss.

\begin{figure}[h]
\includegraphics[width=0.31\textwidth]{ell1.png}
\includegraphics[width=0.21\textwidth]{ell2.png}
\includegraphics[width=0.44\textwidth,clip,trim=0.2cm 0.12cm 0.2cm 0.25cm]{deygout.png}
\caption{First Fresnel ellipsoid for a direct path, with partial obstruction by two buildings (on the left), sampled curved rays used for the clearance computation (on the middle) and extra path loss approximation (on the right).}\label{Fig_ell}
\end{figure}

One can also notice the successive improvements of the simulation method in Figure~\ref{figpdf75x}, from the pure 2.5D method, the 2.5D + 3D method (with 3D reflections on buildings) and the 2.5D + 3D method with first Fresnel ellipsoid clearance, gives the better results up to day.

\begin{figure}[h]
\begin{center}
%\resizebox{0.8\width}{!}{
%\begin{turn}{90}\scriptsize \hspace{2ex}cumulated density function (\%)
%\end{turn}
\includegraphics[scale=1]{Fig11.png}
\caption{Cumulated density function of total electric field in Paris: evolution of the simulation method.}\label{figpdf75x}
\end{center}
\end{figure}


\subsection{Post-processing}

Some post-processing is done to reduce results variability and to generate global exposure maps.

\subsubsection{Hotspots refinement and vertical smoothing}

The main goal of the simulation is to evaluate the electric field everywhere (exposure maps). But it should also be possible to estimate the higher levels (``hotspots''), close to the antennas. One must notice that there is a strong vertical spatial variability of the electric field in front of the antennas. This is due to the narrow vertical beamwidth of the radiation pattern (usually 6 to 9~degrees, while the horizontal beamwidth is $\approx 60^\circ$). As a consequence, the choice of the (discrete) computation points can lead to high variation of the maximum electric field on a vertical line, as illustrated in Figure~\ref{Fig_chaud}, and also miss the higher exposure point. To fix this problem, in the case of a line-of-sight from the antenna to the receiver, a new computation point is added.

\begin{figure}[h]
\includegraphics[trim=0.25cm 0 0 0,clip,width=7cm]{lissage.png}
\resizebox{0.62\width}{!}{\includegraphics[trim=0 0.1cm 0 0.9cm,clip,width=10cm,page=4]{max+lisse.pdf}}
\caption{Effect of the vertical sampling on the maximum exposure point computation in front of antennas and simulated electric field on a column of receivers with regular vertical sampling (without and with smoothing), with added hotspot.}\label{Fig_chaud}
\end{figure}

These extra computation points are dynamically added during computation by the software. For each column of receivers (receivers at the same horizontal location, with different heights), receivers in LOS of an antenna are tagged. The precise height of the maximum exposure point is determined using the radiation pattern and the possible height ranges (between two or more adjacent LOS points) with a free-field formulation.


The thin violet crosses in Figure~\ref{Fig_chaud} are the input computation points on a vertical in front of an antenna. An extra hotspot point (fat red cross) is dynamically added during computation at roughly 15.5~m high. This hotspot captures the maximum exposure level that would otherwise have been underestimated. The interpolation between the points is done using a cubic positive interpolation that allows to estimate the electric field at any height (violet curve).

While this hotspot extension ensures a very stable maximum electric field, it is also quite unrealistic because it is a purely theoretical value located at a very precise point. Hence, a vertical smoothing is applied afterward. Using the previous interpolation, a sliding average with a 1.5 meter window is applied. This smoothing mimics the French measurement protocol where measurement points are averaged over three different heights.
Simulation on Paris showed that the average electric field does not change that much (2\% to 3\% lower with smoothing), while the peak values have more sensitive changes (up to 15\% lower).

\subsubsection{Interpolation of results}

As the electric field is computed on a per-base station basis, the computation points are different from one computation batch to another, and there is some overlapping between the computation areas too. Consequently, the electric field from a computation batch has to be interpolated everywhere in space in order to be able to cumulate electric field from different batches at the same location. This interpolation is done using both the scalar electric field and the path difference to reproduce fading and interferences on exposure maps.

The final rendering is done on a 2 m $\times$ 2 m rasterized image and with square tiles covering the whole French territory.

\section{Results}

\subsection{Computation time}

The validation tests are done with a 14 cores Intel Xeon Gold 6132 CPU machine. The computation on a whole French department (Bas-Rhin) is equivalent to 76 hours on a monothreaded machine but there is a strong dispersion between sites (see Table~\ref{tabtemps}).

\begin{table}[h]
\caption{Computation time per site for full Bas-Rhin department}\label{tabtemps}
\centering
%\resizebox{0.7\textwidth}{!}{
\begin{tabular}{|r|r|r|r|r|r|r|r|} \hline
average& rms& min& max& P50& P90& P95& P99 \\ \hline 03mn 04s& 03mn 17s& 00mn 04s& 45mn 44s& 02mn 24s& 05mn 20s& 06mn 58s& 13mn 02s \\ \hline
\end{tabular}
%}
\end{table}

\section{Comparison with measurements}

\subsection{Methodology}

The electric field has been measured in different (outdoor) locations following ANFR protocol~\cite{protocol}. The measurements evaluate the electric field per band and operator averaged over 6 minutes. The simulations are done with EIRP of transmitters (declared authorized power) and reduction factors are applied to the results to account for the real load. This reduction factor is 4 dB for each technology apart from beamsteering antennas.

The measurement points cover a lot of different situations (LOS, NLOS, rural, urban, \dots). There is still uncertainty in the location of the points and the comparisons are carried ``as is''. Three indicators are used to compare measurements : Pearson correlation $p_{\textrm{\scriptsize lin}}$ (on measured and simulated electric field in V/m), signed mean error $\epsilon_{1,\textrm{\scriptsize lin}}$ (of difference between measured and simulated electric field in V/m) and quadratic mean error $\epsilon_{2,\textrm{\scriptsize log}}$ (of difference between measured and simulated electric field in dB).

If $E_{\textrm{\scriptsize mes},i}$ (respectively $E_{\textrm{\scriptsize sim},i}$) is the measured (resp. simulated) full band electric field for point $i, i \in \{1, 2, \dots n\}$, we have:
\begin{align}
p_{\mathrm{lin}} &= \frac {\sum_{i=1}^n (E_{\mathrm{sim},i} - \bar{E}_{\mathrm{sim}}) (E_{\mathrm{mes},i} - \bar{E}_{\mathrm{mes}})} {\sqrt{\sum_{i=1}^n (E_{\mathrm{sim},i} - \bar{E}_{\mathrm{sim}})^2 \sum_{i=1}^n (E_{\mathrm{mes},i} - \bar{E}_{\mathrm{mes}})^2}}
\left\{
\begin{array}{ll} \bar{E}_{\mathrm{sim}} = \frac{1}{n} \sum_{i=1}^n E_{\mathrm{sim},i} \\ \bar{E}_{\mathrm{mes}} = \frac{1}{n} \sum_{i=1}^n E_{\mathrm{mes},i}
\end{array} \right.\\
\epsilon_{1,\mathrm{lin}} &= \frac{1}{n} \sum_{i=1}^n \left(E_{\mathrm{sim},i} - E_{\mathrm{mes},i} \right) \mbox{(V/m)}
\hspace{2ex}\epsilon_{2,\mathrm{log}} = \sqrt{\frac{1}{n} \sum_{i=1}^n \left(L_{\mathrm{sim},i} - L_{\mathrm{mes},i} \right)^2} \mbox{(dB)}
\end{align}

\subsection{Comparison on a city}

The first comparisons were made in the city of Strasbourg (Bas-Rhin, France), located within a 4~km $\times$ 4~km square. There are 575~antennas and 2910~transmitters. 48~measurement points are scattered within the city. Two datasets were considered: one based on generic data (called B, representative of future simulations) and one based on more precise data (called A):
\begin{itemize}
\item The B dataset uses IGN BDTOPO buildings and ANFR STATIONS database with automatic relocation.
\item The A dataset uses 3D buildings with roofs and improved antennas locations (1/3 of manually positioned antennas using Google Street View, 1/3 of antennas using positions from the city database and the last 1/3 with automatic relocation, see Figure~\ref{Fig_strasbourg2}).
\item The BA dataset is a mixed dataset, using A antennas locations and B dataset buildings.
\end{itemize}

It was not possible to achieve a real ``reference'' dataset as far as antennas positions are concerned, mainly due to different time stamps for measurements and existing antenna databases.

\begin{figure}[h]
\includegraphics[scale=1]{strasbourg_zone.png}
\caption{City of Strasbourg antennas (green circles for precise checked locations and red circles for automatically relocated ones) and measurement points (yellow pins).}\label{Fig_strasbourg2}
\end{figure}

The simulated indicators are shown in Table~\ref{Tab_strasbourg} and the corresponding scatter plot is in Figure~\ref{Fig_strasbourg_nuages}. Because of the uncertainties in the input data (antenna locations and radiation patterns, constant reduction factor for load, receivers locations and geometry~\cite{amsterdam}), the results are on par with previous work~\cite{mithra} and also literature~\cite{munich2d}, acknowledging that those works relied on a controlled emitter (constant power and well-known radiation pattern). It must also be noted that since the comparisons were carried on isolated points, no spatial averaging could be done on measurements, whereas final exposure maps will have a spatial averaging. The scatter plots show a better agreement with the A dataset than with the B one. Scatter plots for free-field computation (overestimating) and the Hata model (underestimating) prove that simulation can handle heterogeneous exposure cases (LOS, NLOS, both occurring from different antennas in a lot of exposure cases).

\begin{table}[h]
\caption{Comparison between measurements and simulation in Strasbourg city (48 points) for A, B and BA datasets.}\label{Tab_strasbourg}
\centering
%\resizebox{0.45\textwidth}{!}{
\begin{tabular}{|l||r|r|r|} \hline
dataset / indicator&$p_{\textrm{lin}}$ (-)&$\epsilon_{1,\textrm{lin}}$ (V/m)&$\epsilon_{2,\textrm{log}}$ (dB)\\ \hline\hline A&0.800&0.094&7.6\\ \hline B&0.537&-0.046&6.3\\ \hline BA&0.776&0.081&7.2\\ \hline
\end{tabular}
%}
\end{table}

\begin{figure}[h]
%\resizebox{0.8\textwidth}{!}{
\includegraphics[scale=0.6]{Campagne PNSE4 Strasbourg AAy.png}
\includegraphics[scale=0.6]{Campagne PNSE4 Strasbourg BBy.png}
\includegraphics[scale=0.6]{Campagne PNSE4 Strasbourg DA+HA.png}
%}
\caption{Scatter plot of measurements (horizontal axis, log scale in V/m) and simulations (vertical axis, log scale in V/m) : A dataset on the left, B dataset at the center, free-field (blue) and Hata (orange) on the right.}\label{Fig_strasbourg_nuages}
\end{figure}


\subsection{Comparison on a full French department}

\subsubsection{Overview}

There are 200~measurement points in the French Bas-Rhin department (67) including 48~measurement points in Strasbourg, and 300~measurement points in the Paris department (75). See Figure~\ref{Fig_67_75}. The first result is that for Paris, all measured points are in computation areas of several antennas (because of the dense urban area), i.e. no measured point has a null simulated electric field. To the contrary, in Bas-Rhin there are 17 out of computation areas. For 15 of them, they were labeled as below the sensitivity of measurement (0.05~V/m). The other points are close to the French - German border, and foreign transmitters are not taken into account. These results validate the precomputation stage.

\begin{figure}[h]
\includegraphics[width=0.290\textwidth]{strasbourg_tout.png}
\includegraphics[width=0.600\textwidth]{paris_tout.png}
\caption{Measured points in Bas-Rhin (left, 200 points) and Paris (right, 300 points).}\label{Fig_67_75}
\end{figure}

\begin{table}[!h]
\caption{Palette index of electric field value for colormap rendering.}\label{tabpalette}
\vspace*{-15pt}
\begin{center}
\resizebox{1\textwidth}{!}{
\begin{tabular}{|c|p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}p{0.9cm}|} \hline
E (V/)&$[0;0.4[$&$[0.4;1[$&$[1;2[$&$[2;3[$ &$[3;4[$ &$[4;5[$ &$[5;6[$ &$[6;9[$ & $[9;12[$& $[12;\infty[$ \\ \hline index&0&1&2&3&4&5&6&7& 8& 9 \\ \hline color& \multicolumn{10}{|l|}{
\includegraphics[width=12.4cm,height=0.25cm]{palette.png}} \\ \hline
\end{tabular}
}
\end{center}
\end{table}


There are 4959~antennas, with 21681~transmitters dispatched in 334~computation batches for Bas-Rhin in a $\approx 4800~\text{km}^2$ territory. There are 8547~antennas, with 41940~transmitters dispatched in 1458~computation batches for Paris in a $\approx 105~\text{km}^2$ territory.

The final exposure maps (see Figure~\ref{Fig_carte} for instance) are rendered using a \emph{viridis} palette, with ten indexes given in Table~\ref{tabpalette}. These indexes are a good trade-off between a linear and a logarithmic scale for comparing results, as it does not emphasize errors on high (respectively low) electric field values like a linear (resp. logarithmic or dB) scale does.



\subsubsection{Distribution of errors}

The distribution of errors between simulation and measurements is illustrated in Figure~\ref{figerr67}. The average error is close to being null, while the distribution shape is close to a normal one. There are a few points where error is quite high. It has been investigated and this could be explained by missing geometry in the case of overestimation (missing large and high cemetery walls, with measurements behind the walls) and by the use of a unique load factor for all antennas. The remaining errors might also be due to incorrect input data (only the overall tilt is known, for instance, and the choice of electrical or mechanical tilt can influence the shape of sidelobes).

In addition to the distribution of relative error in V/m, there is a palette index error. A palette error of zero means that the measured and simulated values are equal. A palette error of one means, for instance, that one of the values is at the lower bound of a class and that the other is at the upper bound, but it could also cover two classes (if both values are in the middle of the class for instance). For Bas-Rhin,
more than 80\% of points have an index error between 0 and 1 (75\% for Paris).

\begin{figure}[h]
\begin{center}
\includegraphics[scale=0.53]{67_6_doit2.pdf}
\includegraphics[scale=0.53]{67_6_doit2.pdf}
\vspace*{-8pt}
\caption{Error distribution in Bas-Rhin (electric field on the left, palette index on the right).}\label{figerr67}
\end{center}
\end{figure}

\subsubsection{Cumulated density functions of electric field}\label{secres}

It is also interesting to observe the cumulated density functions (cdf) of the electric field. The cdf for measurements and simulation in Paris is illustrated in Figure~\ref{figpdf75_1} for the global electric field (including all bands) and on Figure~\ref{figpdf75_2} for specific bands. The 3.5 GHz band (5G only) is not available yet as measurements are too sparse and the simulation method for beamforming antennas is not finalized yet.

The shapes of the distribution are very similar, aside from the maximum value itself.

\begin{figure}[h]
\centering
%\resizebox{0.8\width}{!}{
%\begin{turn}{90}\scriptsize \hspace{2ex}cumulated density function (\%)
%\end{turn}
%\includegraphics[page=1,width=0.48\textwidth,clip,trim=0.7cm 0.5cm 0.2cm 0.9cm]{75_6_doit1.pdf}} \\
%\vspace*{-0.2cm}\resizebox{0.8\width}{!}{{\scriptsize electric field (V/m)}}
%\vspace*{-12pt}
\includegraphics[scale=1]{figures/fig17.png} 
\caption{Cumulated density function of total electric field in Paris (including 5G).}\label{figpdf75_1}
\end{figure}

\begin{figure}[h]
\begin{center}

\includegraphics[scale=1.14]{figures/fig18-a.png} 
\includegraphics[scale=1.14]{figures/fig18-b.png} 
\includegraphics[scale=1.14]{figures/fig18-c.png} 

\vspace*{-8pt}
\caption{Cumulated density function for the 700 MHz / 800 MHz / 900 MHz / 1800 MHz / 2100 MHz / 2600 MHz bands.}\label{figpdf75_2}
\end{center}
\end{figure}


\section{Conclusion and future works}

The method presented in this paper allows to tackle the computation of mobile telephony exposure maps at a full country scale. The method is fast and partial recomputations can be done on a regular basis. The first comparisons with measurements give acceptable accuracy given the uncertainties of input data, the most prominent being the precise antenna locations (as shown by the very similar results between A and B datasets). The results are expected to improve in the future as the quality of input data will increase (precise antenna location from operators database, more detailed buildings database - including roofs, and case dependant load factors).

A large part of the work needed to provide a full France exposure map service has not been addressed in this paper as we only focused on the simulation method. Indeed, the full computation process is automatized and it is an application platform that handles databases, computations and web services for result queries (whether they are electric field at a given location, exposure maps or statistics in cities or countries).

Additional work could be done to improve simulation results. For instance, the Fresnel ellipsoid approach could be extended to one-reflection paths (on the ground or on buildings) and the relocation algorithm will be improved with feedback from initial computations. New validations are carried out on foot tests and drive tests, allowing spatial / time averaging on measurements that could be replicated in simulations.

Finally, for 5G beamforming transmitters, the current simulations are not yet conclusive. The approach used (an envelop of all beams with a load factor, as recommended in~\cite{lignesdirectrices}) seems not to reflect the current behaviour of antennas. One solution to explore could be taking into account local visibility of 5G antennas to determine which beams are likely to be enabled because users (outside on the street or inside buildings) could stand in these beams, and which beams are likely to be disabled, in a simplified approach like the one in~\cite{cstb5g}.

\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

%\nocite{*}


\bibliographystyle{crunsrt}
\bibliography{crphys20230644}
\end{document}
