Cellular Simulation for Distributed Sensing over Complex Terrains

Abstract

Long-range radio transmissions open new sensor application fields, in particular for environment monitoring. For example, the LoRa radio protocol enables connecting remote sensors at a distance as long as ten kilometers in a line-of-sight. However, the large area covered also brings several difficulties, such as the placement of sensing devices in regards to topology in geography, or the variability of communication latency. Sensing the environment also carries constraints related to the interest of sensing points in relation to a physical phenomenon. Thus, criteria for designs are evolving a lot from the existing methods, especially in complex terrains. This article describes simulation techniques based on geography analysis to compute long-range radio coverages and radio characteristics in these situations. As radio propagation is just a particular case of physical phenomena, it is shown how a unified approach also allows for characterizing the behavior of potential physical risks. The case of heavy rainfall and flooding is investigated. Geography analysis is achieved using segmentation tools to produce cellular systems which are in turn translated into code for high-performance computations. The paper provides results from practical complex terrain experiments using LoRa, which confirm the accuracy of the simulation, and scheduling characteristics for sample networks. Performance tables are produced for these simulations on current Graphics Processing Units (GPUs).

1. Introduction

Climate change and natural evolution seriously impact several social and economic aspects, including living conditions, human health, and development. Autonomous observation is a key point for understanding evolution and the management of territories. This is achieved using small stations equipped with processors, sensors, transceivers, and power supplies. Distributed sensors sample physical behavior. They enable analyzing environmental changes and to predict short-term and medium-term transformations.

Wireless Sensor Networks (WSNs) offered several frameworks to connect sensors. Radio standards such as 802.15.4 [1] propose short range mesh connected topologies suitable for indoor or local deployments. Recent innovations include low power long-range radio systems, such as LoRa [2] or Sigfox [3]. Long range means covering a large surface, more risk during transmissions, long signal latencies, and more difficulties coming from the ground topology. Star organization, where a central node addresses remote sensors is the more natural technique to collect data on surfaces within a range of ten thousand square meters.

Complex geographic areas such as shores or islands, hills, valleys oppose physical difficulties to the radio propagation. At the sensor level, the coverage also defines zones where information is accessible, whatever it may be, wildlife including insects, physical elements such as water or gas. Therefore, a precise coverage computation is critical for monitoring efficiency which matches radio connectivity objectives, and sensing accuracy in regards to physical phenomena. This article describes principles of a general method based on geo-localized cellular systems. It explains how a radio coverage or a physical coverage can be computed from the same framework, allowing for obtaining a better match between observation engines and observed phenomena.

A cellular approach was chosen to address the variety of entities appearing in physical situations. Geography variation is one of them, with elevation accidents, rivers, lakes, forest, and shores. By fragmenting geography into cells, it becomes possible to isolate different behaviors, to simulate, then to reassemble these independent behaviors using higher level simulation paradigms, as shown in Figure 1 and [4].

The approach of fine grain cellular structure and coarse grain partitioning appears to be flexible, allowing different components to coordinate. It is also compatible with geographic information since each cell is localized.

The introduction firstly gives an overview on principles behind three components in the methodology: (i) the terrain analysis (Section 1.2), (ii) the signal propagation and coverage (Section 1.3), and (iii) the representation of physical phenomena behaviors (Section 1.4). Generic algorithms for sensing deployments can be built in relation to metrics for these components. This section will refer to a grid of points having geometric and geographic coordinates. Preliminary examples present a large region with mountains and creeks along the seashore, likely to be equipped with sensors. This level of explanation matches application needs at an engineering level.

Further sections will detail tool principles, with the presentation of the cellular methodology (Section 2), then algorithms which simulate the physical case of flooding (Section 3), and long-range radio propagation (Section 4). Appendix A will present validation results obtained by these algorithms.

1.1. Related Work

Use of a cellular automata approach for modeling of natural phenomena in terms of terrain complexity has attracted the interest of several scientists. Climate change impacts and sea level rise, leading to severe weather, salinity intrusion and natural disasters are described in [6,7,8,9,10]. From these studies, proper plans and forethought can support adaptation and resilience. In addition, other publications have investigated hydrological modeling with respect to homogeneous of landform types, for example, surface water flow simulation, complex river system modeling, and flash flood forecasting [11,12,13]. Solutions for problems related to land use, land erosion, and landslide have been proposed by [10,14,15]. Another related research topic is environmental pollution such as oil spill, air pollution, and so on [16,17]. Many natural disasters (e.g., volcano, earthquake, tsunami, snow avalanches, etc.) could be forecasted through computer programs based on cellular technology [18,19,20,21]. Modeling wildlife propagation and urban in relation to food, environmental variables are also research fields for CA [22,23,24].

1.2. A Characterization of Terrain Complexity

Sensor investigation in a zone almost certainly starts with a glance at a geographic map or an aerial photography. If the zone includes mountains, hills, rivers, or shores, two questions will arise, about the sensing point locations, the reachability of an infrastructure network, and the adequacy with the sensing objectives. Design tools will help to take into account the geographic characteristics and possible obstacles to signal propagation, while keeping track of these objectives. Sensing systems are application specific, point to point, and differ a lot from mobile radio communication based on regular coverages.

Terrain complexity metrics allow for measuring ground irregularity, and offer the possibility to split and isolate zones of similar characteristics. Figure 2 displays complexity for a complex zone (Soummam Valley, leading to Bejaia city, the north of Algeria) based on Digital Elevation Models (DEMs) analysis.

Several metrics are known for this. Terrain Ruggedness Index (TRI) and Topographic Position Index (TPI) for DEM provide a quantitative measure of surface turbulence [25]. TRI is the square root of the summed squared deviation in elevation between a cell and its eight neighbor grid cells [26]. TPI is the difference in elevation value of a center cell and the mean of eight adjacent neighbor cells [27]. Both TRIs and TPIs grids illustrate the distribution of terrain heterogeneity and may be displayed in the form of a map as shown in Figure 2. Figure 3 display analysis results for four cases, also showing how terrain characteristics allow separating concerns about sensor layouts. Partitioning the ground according to terrain complexity (as shown in Figure 4) is a preliminary operation for network layout and selection of the pair to pair radio transmissions, to be obtained through computer programs.

1.3. Designing for Long-Distance Radio Coverage

Setup of long-distance radio communication links is difficult if there is an obstacle in the geographical topology.

In smooth terrain areas, radio propagation can be considered as near line-of-sight (NLoS). In this case, the communication distance is expressed in relation to power reduction along the path. Free space path loss (FSPL) equation was proposed for this aim (see Section 5 and Equation (9) in Section 5).

Concerning complex terrains, the choice of a radio position can be application context dependent, topology-dependent, or it can come from a layout algorithm. Thus, it is critical to obtain a description of the zones covered by a signal, plus other characteristics such as propagation delay from the source, and quality of the signal. Section 5 will explain how this can be done using cellular algorithms which mimic the spatial physical propagation over the terrain.

Figure 5 illustrates the results of a coverage computation, from an expected source in the middle of the figure. The irregular yellow shape reflects the difficulty to choose source and destination given that any position between them can interrupt the signal propagation. Deploying a network of several nodes on long distances is not tractable by hand.

1.4. Observing Physical Phenomena

Another component for observation methodology is the necessary focus on the phenomenon to monitor. The nature of sensing varies a lot, depending on application fields. Sensing can be related to biology, climate, geology, human activities, or composition of parameters. Complex terrains will bring a complex behavior, as it is the case for radio signal propagation.

To illustrate the observation problem, let us consider heavy rain in the same region as in Section 1.3. The physical fact to be monitored is water streaming from the mountain to the river, then to the sea. Rain can be predicted from the meteorology services simulations, or just observed using radars, satellites, or ground devices. Causality for water flooding is the nature and shape of the terrain, sharing characteristics with terrain complexity and signal propagation.

Cellular simulation presented in Section 2 will compute propagation of water given an expected, or existing rain history taken from a public service (this case), or extracted from radar networks [28] or satellite observations [29]. Figure 5 displays the water streaming from the mountain down to the river in this particular case of a heavy storm whose profile is shown in Figure 6.

In the resulting maps, high values represent flooding and possible risks. This analysis demonstrates that simulation allows for establishing dependencies between risks and positions. Choice of real sensing positions can be produced by computer programs. For example, the control that appeared in Figure 7 can connect on radio coverage zones and obtain measures from critical positions in the rivers. The difference in communication latencies is another problem in the long range. This is also a place where computer tools can help to produce a schedule of sensor communication automatically.

The flexibility of cells embedding local parameters makes it possible to represent many physical phenomena such as ocean or river wave effects, rain effects and flooding, air or water pollution, species behavior, sound, and alerts reachable zones, etc. Parallel algorithms were designed to simulate rain flooding (Section 3) and propagation of the LoRa radio protocol (Section 5.1), the objective being to obtain high performances on these kinds of problems.

Now, the article explains the geographic space structuration based on cellular systems.

2. Mapping Geographic Space into Cell Systems

2.1. Geographic Map and User Interfaces

Several domains such as physics, medicine, or biology need sensing systems and geometric references. The case of geography includes necessary projections of a spherical surface into flat maps, with many options depending on applications, precision and considered location. Software tools exist for conversions [31]. Designing observation systems implies considering distances and geometry relating to applications, and therefore use of one or several reference systems to manage physical phenomena and their perceptions coherently. This section proposes an explanation of map building and displays supporting the meaning of map contents, and data abstraction for geo-localized information additionally fetched from other sources, including simulations.

Tiled web maps are used to represent geographic data and to display information as flat graphics. This appears in proprietary software, or free access map systems such as OpenStreetMap. In Tiled web maps, the Earth is accessed by a projection rule called Web Mercator [32], and most often by a 3D designation mechanism XYZ used by web browsers and tile servers. The Z parameter provides a zoom factor from the whole Earth to more and more detailed view of the Earth fragments. At a given zoom factor Z, XY gives the 2D index of a tile inside the rectangle of tiles of this level Z. Based on this, tiles of graphical information are composed and delivered by servers to clients, usually in the form of $256\times 256$ pixel images (Figure 8 presents a map composed from such tiles). These images are composed by a server rendering engine that extracts objects from a database and draws them according to some style (see Figure 9).

The browser Quickmap [33] supports standard map tiles, including OpenStreetMap, and a variety of other items, either for maps or aerial images. This tool also allows for describing sensor systems, communication links, and mobile trajectories, for example, Low Earth Orbit (LEO) satellites [34].

At a map browser level, users specify or observe geometric points relative to a zoom level. Real points are geographical positions, with coordinates specified according to the current projection system. In the case of Quickmap, they are double precision floating point numbers for latitude and longitude. Another concern is the necessity to use real distances in meters, over the Earth surface, for a task such as coverage computation. In the case of Web Mercator transformations between addressing mechanisms for (lat, long) absolute geographical specifications (x, y) distance in meters, (x, y) pixels, and $XYZ$ tile access operations can be managed with formula from [35].

2.2. Definition of Cells

Geographical cells can now be defined as objects grouping local data, local system behaviors, and graphic representation. This is described as an object-oriented class, as variables, and methods that are strictly local to the cells (Figure 10):

• Cell locations are bounded to geometric locations on maps, in relation to tile containers, and when possible to geographical locations.
• Basic cell contents are extracted from the map or image fragment as supported by the browser tool.
• Cell content extensions are obtained from external databases. Most of the cases are digital elevations, and also climate, or weather characteristics and historical data.
• Cell size is chosen to match a particular physical phenomenon or sensing requirements.
• Cell behaviors are local procedures operating on a cell state. They need to be programmed to produce simulation data that in turn can be sent back to databases, or displayed on tiles.

Examples in Section 1 illustrate how geographical abstraction, representation, and simulation can interact (Figure 2, Figure 5 and Figure 7).

Given a partitioning of a geographical zone into a cell system, we need to relate the observation system and physical behavior. The possibility to simulate physical processes will be of great help to correctly characterize a layout of sensors in connection with physical evolution.

Cellular automata (CA) allow for modeling physical activities as processes that exchange information and evolve according to transition rules. The cell concept binds geographical fragments to such processes whose assembly is generated automatically.

2.3. Cellular Automata Principles

2.3.1. Synchronous Systems

Cellular automata were invented by John Von Neumann and colleagues with the aim to build a self-reproducing machine abstraction [36]. To support this goal, a two-dimensional space was configured with automata governed by a small set of states. This representation of space was used in several scientific domains and bound to physical behavior having similar properties [37]. CA can be described, and specified as a discrete space which associates cells.

In synchronous cellular automata, cells evolve, step by step, following a discrete time. Such CA have been described in a variety of languages and executed on specific machines [38] following four simple patterns:

• The cellular space is represented by an assembly of similar cells. A common notion of neighborhoods defines local communications following observed physical dependencies. The spatial organization can be either regular or irregular, possibly with disconnected subsystems as shown in Figure 11, item 3.
• The evolution of each cell is defined in a set of states as observed in a real system (quantities, colors, boolean). Change of states are operated by procedures associated with representing transition rules from step to step: $Stat{e}_t\to Stat{e}_{t+1}$.
• The neighborhood represents physical dependencies, for example, signal propagation, or downward flooding. These dependencies are connectivities from cell to neighbor cells. Common neighborhoods are Von Neumann and Moore with four and eight cardinal directions, respectively. Item “Process architecture” in Figure 11 illustrates a Moore neighborhood.
• The transition rule defines the behavior of each cell evolution under influence of its neighborhood and local sensed influences. The state of the whole cell system synchronously changes, time step by time step.

Non deterministic behaviors include random variation at the physical level. This was used in lattice gas simulation [39] and to represent life cycle and species evolutions [40].

2.3.2. Variability in Large Systems and Asynchronism

Considering large systems, synchronism and massive data parallel execution can be an obstacle, both in behavior modeling and performances because of sparse data spaces.

Asynchronous cellular automata were proposed to represent reactive systems where events are propagated in a way similar to Communicating Sequential Processes [41]. However, in the case of physical simulation, sampling and time references are often mandatory. Variability of computation frequency can be obtained by isolation and simulation of sub-systems. This was implemented on GPUs as related in [42].

The High Level Architecture (HLA) also allows for sequencing cellular sub-systems at a different speed and to provide data exchanges between them. For example, Ref. [4] shows a countryside simulation with forest, forest fire, sensors and river pollution composed together.

Thus, we can consider the grouping of cells into synchronous subsystems as an efficient way to manage variability. Composing simulations at a high level is shown in Figure 1, and was technically discussed in [4].

2.4. Cellular Automata Parallel Execution Models

A sequential execution model will read states in an array $A_t$ for time t, and loop over, writing a similar array $A_{t+1}$. Once the step is completed, the two arrays are exchanged, and a new turn begins. CA is easy to parallelize preserving this behavior, grouping operations together and observing turns completion. Parallelism is mandatory because problems are large, if not huge, a lot of small cells can be critical for simulation precision, and some applications require examining very large geographical regions.

• The synchronous distributed messaging model [43] can support parallel computation by associating cells with communicating processes. In this case, process progress by locked steps, based on messages being sent and received to or from neighbor nodes. The steps are split into two phases, one for communications with the neighborhood, the other to execute the transition rule. Figure 12 shows an internal node representation and the outside connection with three input and output links. This model does not need to specify the relative speed of processes, and can therefore be used for multi-cores or supercomputers.
• Another way to take advantage of parallelism is to use data parallel Single Instruction Multiple Data (SIMD) processors to execute a group of processes simultaneously. Current graphics accelerators propose solutions up to two thousand processors working concurrently, and exchanging data synchronously in shared memory. State spaces must be copied to the accelerator memory. Then, a loop of steps can be run completely on the acceleration, which is very efficient.

Production of code from cell systems to this schema was done for two targets:

• Asynchronized Occam communicating processes [44,45] targeting KRoC compiler and multi-cores [46],
• CUDA code production [47] targeting NVIDIA tools and accelerators.

2.5. Cellular Systems Workflow Organization

Most of the work described in this article is supported by a set of dedicated tools from the University of Bretagne Occidentale (UBO). The design flow can be summarized as follows:

• Zone selection is done by moving a graphical window anywhere, with any level of zoom. The tool extracts graphic tiles, and displays the contents.
• Cell segmentation is obtained by splitting the view into rectangles of a given size specified as a width × height value. Cells will carry an image and a geographical location from the underlying image.
• Binding cell together and producing a cell system implies the choice of a connectivity (Moore, Von Neumann), and possibly filtering cells by colors or elevation. This step also injects external values from a variety of sources, including elevation.
• Adding behavior programs the cellular system at the local level, given a cell system architecture that step 3 (Figure 11) automatically produces.

Steps 1 to 3 are interactive and can be achieved in minutes. Step 4 is a concurrent programming activity specific to a concurrent platform which requires elaboration of CA transition rule and coding for the target platform.

Physical phenomena simulation can reveal places of interest for sensing. In the case of flooding, simulation track accumulation and circulation of water.

3. Heavy Rain Simulation on a Complex Terrain

To illustrate cell system interactions and behavior, we use the example of a zone receiving and flooding rain, as discussed in Section 1.4.

3.1. Managing Space

Space is defined by geographical coordinate bounds, and a possible selection operated on cells. Selected cells have common properties, such as an elevation above or under a threshold level (Figure 11, step C), some color characteristics, or some signature produced from observed parameters. Thus, the cell system can be rectangular or have arbitrary shapes and isolated subsystems.

Practically, cell systems first appear as arrays of objects that carry geometric coordinates, geographic coordinates, and a pixel array extracted from the original image.

Besides the zone structure, designers need to specify neighborhood, representing local physical influences and producing the necessary connectivity between processes. We can just admit that the space animation will be obtained by messages exchanged between neighbors. Execution software and hardware will allow adaptation to the intermediate communicating process model.

The case of rain reveals several interactions:

• Millimeters of water falling on the ground. For this simulation, we admit that the quantity can vary over time, but will remain uniform.
• Water disappearing locally for reasons such as absorption or evaporation.
• Water passed locally from cell to cell according to elevation differences.

Other cases include different specific problems for signal propagation (Section 5), sound propagation [48], or insect swarms behavior [40]. For flood modeling, Von Neumann neighborhood is convenient. Figure 13 shows such a neighborhood, with the North cell lost, and some elevation differences appearing above and below a center cell.

3.2. Transition Rule

Each cell will receive rainwater, dispatch part of this water to neighbor cells with lower elevation, and absorb another part. A complete study will take into account the ground specific characteristics, and current weather (ECOCLIMAP [50] supporting French meteorology AROME model provides more than 20 parameters for cells of 1 km × 1 km). While real dependencies are water leaking from cell to cell, the abstract behavior is represented by messages sent and received to and from neighbors. Vertical water behavior is directly encoded in the transition rule that can be specified as follows:

• t: time t represented by a step number,
• $Q_t$: water quantity in a center cell at the beginning of time step t,
• ${\alpha }_t$: rainfall at time t,
• $\beta$: percentage of water remaining on each cell after each step,
• $cel{l}_i.elevation$: elevation value of $cel{l}_i$,
• $cel{l}_c.elevation$: elevation value of center cell,
• ${\delta }_i$: the difference of elevation between neighbor $cel{l}_i$ and the center cell,
• $\Delta$: sum of all elevation differences,
• $out{F}_i$: amount of water out coming from center cell to neighbor $cel{l}_i$,
• $in{F}_i$: amount of water in coming to center cell from neighbor $cel{l}_i$.

The water quantity of center cell at time $t+1$ after local absorption is

$$remaining=Q_t\times \beta ,$$

(1)

but the center cell receives rainfall within its area the following:

$$rain={\alpha }_t.$$

(2)

The quantity of water coming from neighbors is taken into account

$$received=\sum _{i=1}^{n}in{F}_i,$$

(3)

where n is a number of neighbors of the center cell.

A center cell also has to distribute water to surrounding neighbors in the proportion of elevation differences. This is calculated as follows:

$${\delta }_i=cel{l}_c.elevation-cel{l}_i.elevation.$$

(4)

It is obvious that the water of a center cell only flows to neighbor cells with lower elevations so that only positive values of $delt{a}_i$ are used to compute the total proportion of all differences

$$elevation Above=\Delta =\sum {\delta }_i,\forall {\delta }_i>0.$$

(5)

Equation (6) figures out the amount of water that a cell distributes to its lower neighbors. For higher elevation neighbors, with ${\delta }_i<0$, there is no water flowing out, $out{F}_i=0$:

$$sen{t}_i=out{F}_i=Q_t\times ({\delta }_i/\Delta ),\forall {\delta }_i>0.$$

(6)

Eventually, the center cell updates its own quantity of water by adding and subtracting. The whole system executes a synchronous transition from t to $t+1$:

$$Q_{t+1}:=remaining+rain+received-\sum _{i=1}^{n}sen{t}_i.$$

(7)

Notice that this transition rule needs to know the neighbor relative elevations $cel{l}_i.elevation$. In reality, physical rules will discover these relations naturally. In the case of simulations, a preliminary procedure called neighborhood discovery is executed to establish initial knowledge such as the number and the elevations of cells around.

4. Parallel Algorithms for Cellular Long-Range Coverage Computations

4.1. General Idea

Once a geographic space has been selected, a concern is to compute reachable cells in the line-of-sight from an emitting position. According to Section 2.4, we admit that each cell is represented by a process, and simulation is achieved by cycling on a synchronous parallel program: communication with neighbors, updating the local state, and preparing next cycle communications. The line-of-sight (LoS) is a ray broadcast in any direction from a root emitter. Propagation will be stopped or modified by ground obstacles such as hills, valley, etc. The simulation parallel algorithm mimics the physical behavior, by propagating the signal inside a spanning tree rooted at the emitter cell, and covering progressively all the space in concentric “circles” as shown in Figure 14.

Each new step in the algorithm adds a new circle, and the computation finishes in $2\times log\left(n\right)$ steps where n is the number of cells. During ray propagation, the ground profile is collected into a route. A route profile is shown in Figure 15. Routes are completed progressively based on positions and elevations. Thus, each cell can decide if the emitter is visible or not by comparing its slope to root to previous ones in the received profile, as shown in Algorithm 1.

Algorithm 1 Setup visible nodes from emitter based on the LoS condition.

4.2. Vertical Model

Signal propagation progresses in steps according to time. It builds progressively a route from the emitter to any other cell, as an array holding traversed locations and signal parameters (Figure 15). As it is the case for most of the cellular simulation algorithms, nodes start with a discovery stage, where geometric coordinates are exchanged to bind communication channels to cardinal directions. This is mandatory to deal with the irregularity of cell system shapes as shown in Figure 11, step (3), and also to maintain horizontal signal direction, as discussed in Section 4.3.

To start the building a propagation, an emitter has to initialize a tracing route with its own identity and geographic location. Other cells start their processes with an empty tracing route. Note that the tracing route, in this system, is the payload of the message for communications within the cell system. The root cell then sends its message to all neighbors via eight directions from the Moore neighborhood. Each cell accepts only one message and then assigns the corresponding owner of the accepted message as its parent. This cell must send an acknowledgment message back to the parent cell to confirm the dependency relationship between them. The next step in the transition rule is to insert its own location into the tracing route, then to pack routes and spread out to neighbors. If side is the maximum between width and height, then, at most, side rounds are necessary to cover an entire zone.

4.3. Horizontal Model and Directed Breadth-First Search

The original Breadth-First search (BFS) algorithm builds a tree covering a whole network starting from a root position. For radio propagation, this tree would allow reaching any cell in a minimum number of steps, propagating a route to this cell without managing signal horizontal directions. In the case of radio signals, there is the necessity to guide the route horizontally to maintain the best approximation of a straight line. A sequential algorithm that reduces horizontal errors to a minimum was suggested by Bresenham to draw graphic lines on bitmaps [51]. In our case, it was preferred to find a parallel distributed approach based on Directed BFS (DBFS) that could also match particular geographic or electromagnetic dynamic considerations. Figure 16 shows DBFS routes resulting from this algorithm according to an emitting node and some target receivers.

At the difference of a BFS, the cell transition function makes a decision on receipt of route messages based on the source position and its own source position. It adds its data to the route, and it forwards the route to neighbors of interest. The neighborhood was discovered at the initialization stage, associating link indexes and cardinal directions.

5. Radio Signal Propagation on Complex Terrains

Radio propagation models are used to describe the qualifications and reliability of links using radio frequency. These models were based on the physics of the diffusion of electromagnetic waves with respect to both constructive and destructive interference of environment. They can be practically applied to parameters inside routes transmitted from cell to cell, according to Section 4.

Free space path loss (FSPL) is the simplest model based on the Friis transmission equation. By assuming unobstructed along the transmission path, FSPL represents the line-of-sight decay of an electromagnetic wave as a function for distance only. This model exposes significant limits for complex terrain areas. Empirical path loss models were proposed to cope with the effects of turbulent terrains. An advantage is that these models attempt to estimate power loss as a function of distance and radio frequency taking into account the terrain heterogeneity as well as the effects of real transmission environment [52,53].

5.1. Free Space Path Loss Model

In terms of radio communication, a Free Space Condition is a region where there is no obstacle along the propagation path of radio waves. A radio signal is emitted by a transmitter and then it propagates in any direction at the speed of light. In this case, the signal energy can be received by an antenna in inverse proportion of the distance away from the signal source. In addition, the received power depends on the transmitted power, and gains of both transmit and receive antennas. Fading margin and other loss such as system loss, cables, and connectors are also considered. The power of the signal lost on the path is called the Free Space Loss represented by

$$FSPL=20lo{g}_{10}\left(\frac{4\Pi d}{\lambda }\right).$$

(8)

Algorithm 2 is derived from Equation (8) in order to describe how to calculate the attenuation of signal power due to free space path loss. As a consequence, the free space power received is given by Equation (9), the Friis free space equation [54]:

$$P_r\left(d\right)=\frac{P_t{G}_t{G}_r{\left(\lambda \right)}^2}{{\left(4\pi \right)}^2{d}^{2}L},$$

(9)

where

• $P_t$: transmitted power,
• $P_r\left(d\right)$: received power,
• $G_t$: transmitter antenna gain,
• $G_r$: receiver antenna gain,
• d: distance between transmitter and receiver in meters,
• L: system loss factor,
• $\lambda$: wavelength in meters.

Algorithm 2 Received signal power in the Free Space Path Loss Model.

1: procedure SignalPower

2:   Input: txPower, waveLength, distance

3:   Output: receivedPower

4:   for each neighbor of $cel{l}_i$ do

5:      if neighbor.visible then

6:      receivedPower ← txPower × (SQR(waveLen /(4.0 $\times valuePi$ ))× SQR(1/distance))

7:      end if

8:   end for

9: end procedure

5.2. Single Knife-Edge Diffraction Model

Diffraction appears around objects such as buildings, vegetation, etc. The losses caused by these obstacles can be described by a Knife Edge Diffraction model. This is considered as a function of the path difference around the obstacles and could be explained by Fresnel zones. Diffraction loss of signal after a knife edge is calculated as a function of the Fresnel parameter v:

$$v=h\left(\sqrt{\frac{2}{\lambda }(\frac{1}{t_t}+\frac{1}{d_r})}\right),$$

(10)

where h is the height of obstacle. $d_t$, $d_r$ are the distances from the obstacle to the emitter and the receiver, respectively. Algorithm 3 demonstrates the computation of received power within a shadowing area because of a single knife edge obstacle:

$$L\left(v\right)=6.9+20log(\sqrt{{(v-0.1)}^2+1}+v-0.1).$$

(11)

Algorithm 3 Received signal power in the Single Knife-Edge Diffraction Model.

1: procedure diffractionLoss    2:   Input: fsplPower, h, $d_1$, $d_2$    3:   Output: diffPower    4:   v = h * SQRT((valuePi/2)×((1/$d_1$) + (1/$d_2$)))    5:   if v < 0 then    6:      diffLoss ← 0    7:   else    8:      if v < 2.4 then    9:      diffLoss ← 6 + (9 × v) + (1.27 × SQR(v))   10:      else   11:      diffLoss ← 13 + 20×$lo{g}_{10}\left(v\right)$   12:      end if   13:   end if   14:   diffPower = fsplPower - diffLoss   15: end procedure

In Algorithm 3, receivedPower is defined as the signal strength in dBm obtained at a receiver. In addition, sensibilityPower is the minimum value of signal strength in which condition a receiver can detect and demodulate a message successfully.

5.3. Okumura–Hata Model

This is a dedicated radio propagation model to predict the path loss regarding geographic environments such as urban, suburban, and open areas [53,54]. The formula for the Okumura–Hata model is given by

$$L=69.55+26.16lo{g}_{10}\left(f\right)-13.82lo{g}_{10}\left(h_b\right)-a\left(h_m\right)+(44.9-6.55lo{g}_{10}\left(h_b\right))lo{g}_{10}\left(d\right)-K,$$

(12)

where

• $f_{MHz}$: center carrier frequency of transmission band in MHz.
• $h_b$: antenna height of base station in meter,
• $h_m$: antenna height of mobile node in meter,
• $d_{km}$: distance in kilometers,
• $a\left(h_m\right)$ for each type of area and K, see

  Table 2. List of $a\left(h_m\right)$ and K values for each type of area [54].

  Type of Area	a$\left({\mathit{h}}_{\mathit{m}}\right)$	K
  Open		$4.78{\left[\mathit{l}\mathit{o}{\mathit{g}}_{10}\left(\mathit{f}\right)\right]}^2-18.33\mathit{l}\mathit{o}{\mathit{g}}_{10}\left(\mathit{f}\right)+40.94$
  Suburban	$[1.1\mathit{l}\mathit{o}{\mathit{g}}_{10}\left(\mathit{f}\right)-0.7]{\mathit{h}}_{\mathit{m}}-[1.56\mathit{l}\mathit{o}{\mathit{g}}_{10}\left(\mathit{f}\right)-0.8]$	$2{\left[\mathit{l}\mathit{o}{\mathit{g}}_{10}(\mathit{f}/28)\right]}^2+5.4$
  Small city		0
  Large city	$3.2{\left[\mathit{l}\mathit{o}{\mathit{g}}_{10}(11.75{\mathit{h}}_{\mathit{m}})\right]}^2+4.97$	0

  .

The model is valid for radios using carrier frequencies from 150 to 1500 MHz with an effective height of antennas from 30 m to 1000 m and distances ranging from 1 km to 100 km. Table 2 presents list of $a\left(h_m\right)$ and K values for different terrain types.

5.4. Correctness of the Different Radio Communication Models

The terrain was a mountain area investigated by applying the three models during simulations. The computation results are shown in

Table 3. Comparison in number of points found in coverage from a total of 96,431 points.

Model	Points
Free space path loss	19,537
Single knife-edge diffraction	22,305
Okumura-Hata	17,621

. Figure A5 displays the coverage obtained by the model of single knife-edge diffraction, and effective coverage measured by a mobile.

6. Conclusions

To optimize communication ranges, two major concerns are radio waves blocked by obstacles and propagation loss in proportion with a point-to-point distance. Terrains to be covered can be classified and partitioned as smooth, or rough according to metrics (Section 1.2). In smooth terrain areas, radio transmission can be considered as a near line-of-sight (NLoS). As a result, the issue of topologies can be almost neglected and the communication distance is expressed in relation to power reduction along the path. The free space path loss (FSPL) equation was proposed for this aim (Section 5). For example, using LoRa technology for a point-to-point connection in NLoS conditions, the distance for transceiving data is expected to reach hundreds of kilometers [55].

In complex terrain areas, heterogeneity of topologies seriously impacts the quality of radio links, especially for low power and long-range communication networks. Different propagation models allow for representing the median of the expected path loss such as the Longley–Rice model, the ITU model, and the Okumura–Hata model, as described in Section 5.

Appendices give details on experimental measurements, with the communication distance checked to be around 5 km in complex urban area and up to 20 km in rural and shore ones. Experiment results shown in

Table A1. Three areas with different complexity terrains. Resolution and Actual size columns give the size in pixels associated with actual size in meters for each cell in regular grid data. The maximum communication in each experiment is shown by the column Range. The Match and Mismatch column provides a number of visible points which are right matched and mismatched, respectively, between simulation results and obtained values in real measurements.

Area Name	Description	Resolution	Actual Size	Range	Match	Mismatch	% Error
Albert 1er	Urban area	$5\times 5$ pixels	24 m	3 km	64	13	20.21%
Plougastel	River and its banks	$5\times 5$ pixels	96 m	9 km	129	4	3.10%
Trevezel	mountain area	$5\times 5$ pixels	191 m	11 km	45	8	17.78%

confirm the interest of a computer-aided approach. It may be concluded that communication in a LoRa network is strongly dependent on the considered environment. Hence, the actual coverage prediction must take topographic complexity into account as a critical factor. An assessment of physical simulation was also conducted for flash flooding problems (Section 3). This algorithm was applied to the practical case of an intense rain found on 3 June 2018, and compared with flooding observations. The rain occurred in a rough terrain urban area and the simulation was able to retrieve major flooding places and level of water. There is a major interest in these kinds of physical simulations that reveal places of interest for sensing and establish causality between events.

The software developed in this project is easy to use for cellular system generation. Implementing cellular transition rules, tuning and verifying these rules necessitate in-depth investigations. Several domains were investigated, from sound propagation to insect behaviors. Thus, the methodology appears very general and flexible (see

Table 4. A characterization of space and signal relations.

Datagrid	Model	Resolution (m)
SRTM30	2D/3D	30 (https://lta.cr.usgs.gov/SRTM1Arc)
SRTM90	2D/3D	90 (https://blogs.esri.com/esri/arcgis/2015/06/26/terrain-3d-now-with-global-srtm-30-meter-content/)
Weather broadcast	2D/3D	3000 (https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/numerical-weather-prediction)
LiDAR	3D	5 (https://content.meteoblue.com/en/research-development/data-sources/nmm-modelling/model-domain)
Coastal ocean	2D/3D	1200 (https://pdfs.semanticscholar.org/55be/a487827c28aaaf713017c499e4f33aed62fd.pdf)
Earth magnetic	2D/3D	3700 (http://0-www-sciencedirect-com.brum.beds.ac.uk/science/article/pii/S0377042798002465)
Ocean acoustics	2D/3D	2500 (https://www.ngdc.noaa.gov/geomag/emag2.html)
		0.5 (Peter Wille. Chapter 5: The Sea Floor—Natural Formations.
Road traffic noise	2D/3D	In Sound Images of the Ocean: in Research and Monitoring.
		Springer-Verlag Berlin Heidelberg, 2005. ISBN 978-3-540-27910-5)
Pickcell	2D/2.5D	30

). Parallel programs implement described algorithms very efficiently (

Table 5. Execution times for three processes. There are 58,725 cells in a regular squared grid with a cell size 3 × 3 pixels corresponding to actual area 115 × 115 m.

Processes	Number of Rounds	Time (s)
Compute ruggedness index	1	0.026
Flood simulation	32	0.103
Coverage prediction	261	0.211

and

Table 6. Execution times for flood simulations which are performed on Linux Ubuntu PCs with Intel$\circledR$ Core$™$ i3-4005U CPU @ 1.70GHz × 4, 8 GiB DDRAM, card NVIDIA GeForce 820M (96 CUDA Cores), Intel$\circledR$ Core$™$ i7 CPU 920 @ 2.67GHz × 8, 4 GiB DDRAM, card NVIDIA GeForce GTX 680 (1536 CUDA Cores) and Intel$\circledR$ Core$™$ i7-7700K CPU @ 4.20GHz × 8, 16 GiB DDRAM, card NVIDIA GeForce GTX 1070 (1920 CUDA Cores).

			Execution Time
Resolution (Pixels)	Number of Cells	820M	GTX 680	PGTX 1070
3 × 3	58,725	28.168 (ms)	6.5269 (ms)	1.3373 (ms)
5 × 5	21,060	10.411 (ms)	2.1728 (ms)	518.64 ($\mu$s)
10 × 10	5226	2.5234 (ms)	454.59 ($\mu$s)	83.696 ($\mu$s)
15 × 15	2340	1.1155 (ms)	181.48 ($\mu$s)	65.916 ($\mu$s)
20 × 20	1287	602.47 ($\mu$s)	171.19 ($\mu$s)	62.085 ($\mu$s)

). Efficiency is mandatory if space exploration strategies are to be developed.

Several research opportunities are created by this work:

• Synthesizing long-range mesh radio networks would allow for covering very large surfaces and require optimization of relay placement. Powerful parallel algorithms can help to overcome the coverage problem complexity.
• Line-of-sight propagation and water spreading are examples of common spatial behaviors found in nature. Many biological, mechanical or physical processes have similar properties. Sensing of environment could take benefits from generic frameworks helping to develop and compose such simulations.
• Applicability for large problems available for the public can be obtained by publishing web interfaces to services dependent on two separated concerns: zone selection and cell resolution, and cellular system libraries.