number=N, year=YYYY, firstpage=xx, lastpage=yy, doi=10.5311/JOSIS.YYYY.II.NNN,
received=, returned=, revised=, accepted=,

# A Bayesian approach to location estimation of mobile devices from mobile
network operator data

Martijn Tennekes  Department of Methodology, Statistics Netherlands, The
Netherlands  Yvonne A.P.M. Gootzen  Department of Methodology, Statistics
Netherlands, The Netherlands

###### Abstract

Mobile network operator (MNO) data are a rich data source for official
statistics, such as present population, mobility, migration, and tourism.
Estimating the geographic location of mobile devices is an essential step for
statistical inference. Most studies use the Voronoi tessellation for this,
which is based on the assumption that mobile devices are always connected to
the nearest radio cell.

This paper uses a modular Bayesian approach, allowing for different modules of
prior knowledge about where devices are expected to be, and different modules
for the likelihood of connection given a geographic location. We discuss and
compare the use of several prior modules, including one that is based on land
use. We show that the Voronoi tessellation can be used as a likelihood module.
Alternatively, we propose a signal strength model using radio cell properties
such as antenna height, propagation direction, and power. Using Bayes’ rule,
we derive a posterior probability distribution that is an estimate for the
geographic location, which can be used for further statistical inference.

We describe the method and provide illustrations of a fictional example that
resembles a real-world situation. The method has been implemented in the R
packages mobloc and mobvis, which are briefly described.

###### keywords:

mobile network operator data, geographic location, present population,
Bayesian statistics

##  1 Introduction

Mobile network operator (MNO) data have shown to be a rich potential source
for official statistics, in particular on present population [7, 2, 6, 18, 43,
35], mobility [21, 13, 3, 8, 30, 41, 15, 44, 9], migration [24, 42, 23], and
tourism [7]. Such statistics can be used to support a wide range of policy
issues, for instance regarding mobility restrictions in order to delay and
contain the COVID-19 pandemic [12, 28, 4].

Two types of MNO data are needed for statistical inference: cell plan data,
which contains information about the mobile communication network
infrastructure, and event data, which are data from transaction events between
the network and mobile devices.

A mobile communication network is also called a cellular network, where each
cell enables mobile communication for a specific land area. The cell plan
specifies the geographical location of the cells and their physical
properties. Two types of cells can be distinguished [29]: a macro cell that is
placed in a cell tower or on top of a roof, which has a theoretical range of
about 30 kilometers, and a small cell that is used to enable mobile
communication inside buildings and in dense urban areas, and has a theoretical
range of two kilometers.111A small cell is a general name for micro, pico, and
femto cells, which have theoretical ranges of respectively 2000, 200, and 20
meters. An important cell property is the propagation direction. A cell can be
directional or omnidirectional. In practice, most real-world deployed macro
cells are directional, often covering an angle of about 120 degrees whereas
small cells are usually omnidirectional.

There are generally two types of event data, depending on which transaction
events are contained. Event data that are used to calculate the costs in order
to bill customers are called Call Detail Records (CDR), which contain events
related to active mobile phone use, such as initiating a call 222CDR data
contain records about calls (initiating and receiving), SMSes (sending and
receiving), and mobile data usage. Note that in several studies the term CDR
is used for data that only contains call and SMS events, and alternative terms
Data Detail Records (DDR) or Event Data Records (EDR) are used for data that
also include mobile data usage events.. Event data that also contains passive
events such as location updates, are called signalling data. These data are
used by the MNO for network analysis and optimization. Signalling data are
usually much richer than CDR data and are therefore recommended for
statistical inference.

The geographic location of devices is important for statistical inference on
MNO data. Advanced location estimation techniques are used to facilitate
mobile communication, but these require radio measurement data, such as the
Received Signal Strength (RSS), the Signal to Interference plus Noise Ratio
(SINR), and Timing Advance [16, 17, 45, 5]. However, these data are usually
not stored in event data. The only location related variable that is always
available in event data is the identification number of the serving cell.

The majority of studies on mobile network data use Voronoi tessellation [7] to
distribute the geographic location of logged events. The geographic area is
divided into Voronoi regions such that each Voronoi region corresponds to the
geographic location of a cell tower and each point in that region is closer to
that cell tower than to any other cell tower.

There are a few downsides to using Voronoi tessellation to estimate the
geographic location of devices. First of all, it assumes that all cells are
omnidirectional. As described above, large range cells are often directional.
The second downside of Voronoi tessellation is that it does not take other
cell properties into account, such as power, height, and tilt. Third, the
coverage areas of cells overlap in reality, especially in urban areas. This is
because of load balancing; if a cell has reached full capacity, neighbouring
cells that also have coverage are able to take over communication with mobile
devices. This means that a mobile phone is not always connected to the nearest
cell nor to the cell with best signal. In urban areas, a mobile phone switches
frequently between cells333There are several smart phone apps that show where
the connected cell is located, e.g. Network Cell Info [25].. The fourth and
last downside of using Voronoi tessellation we would like to mention is that
it does not take into account where people are expected to be.

A couple of variations of the Voronoi algorithm have been proposed to overcome
some of these limitations [11, 6, 34]. One improvement is to shift the
locations of the Voronoi points from the cell tower locations towards the
direction of propagation. Alternatively, when the _coverage area_ is known for
each cell, i.e. the area which is served by the cell, the location of the
centroids of these coverage areas can be used as Voronoi points. Another
improvement is to create a Voronoi tessellation for cells with a large range,
and subsequently assign each small cell to the Voronoi region they are located
in. The Voronoi method can be extended with auxiliary data sources, such as
land use, to improve the geographic location of devices [14].

An alternative to Voronoi tessellation is a Bayesian approach. Previously
proposed Bayesian approaches assume availability of SINR data or Timing
Advance data [45, 27]. We propose a modular Bayesian framework to estimate the
geographic location of devices, which consists of two main modules, namely the
_location prior module_ and the _connection module_. The former employs a
priori information about where devices are expected to be to produce a
_location prior_ , and the latter uses network information to estimate the
probability of connection to cells to produce a _connection likelihood_. The
model is generic in the sense that various data sources and methods can be
used for both of these modules. Timing Advance data can taken into account in
our modular approach, but it is not required.

Other Bayesian approaches often assume a uniform distribution for the location
prior [17, 45]. While a uniform prior generally speaking represents an
"objective" prior, we argue this is not true for location estimation of mobile
phones. A main contribution of this paper is the exploration of the effects of
various prior distributions. The most important one is the usage of land use
data. This is arguably the most straightforward option, since more devices are
expected to be in certain land use categories, such as urban areas, than in
other land use categories, for instance grasslands. However, any data source
that contains information about where devices are expected to be can be used.

The connection likelihood describes the estimated probability that a device is
connected to a certain cell given its actual location, taking potential
overlap between cell coverage areas into account. For this component, we
propose a signal strength model which models the propagation per cell using
physical properties of the cells, such as height, direction, and tilt.
However, the Voronoi method and each of the aforementioned variations can also
be used here.

The Bayesian framework can be extended iteratively with modules based on
available auxiliary data. One such module is the Timing Advance module, which
incorporates a variable from which the distance between a mobile device and
its serving cell can be estimated.

The outcome of the Bayesian framework is the _location posterior_ , which is
an estimate of the location of a device given that it is connected to a
certain cell. This location is not a single point, but a geospatial
distribution. For implementation, we recommend to overlay the area of interest
with a grid of square tiles, typically 100 by 100 meters. The location
posterior then defines the probability that a device is located in a grid tile
given that it is connected to a certain cell. This location posterior can be
used for further statistical inference.

A similar Bayesian location estimation approach has been developed
simultaneously with our work [27]. Though their overall approach seems similar
at first glance, there are a two significant differences. First, the proposed
connection likelihood relies on aggregated Timing Advance data, which are
often unavailable. The connection likelihood method that we propose uses cell
plan data which are always available because these data are needed to build
and maintain the mobile communication network infrastructure. When
(aggregated) Timing Advance data are available, they can be used in our Timing
Advance module. Second, the authors of [27] advise against using a priori
information such as land use, raising the argument that it is not objective.
We argue that there is no such thing as an objective prior and will show that
a non-uniform prior has added value for estimating the geographic location of
devices.

The outline of the paper is as follows. In Section 2 we describe our modular
Bayesian framework. The location prior, connection likelihood, and timing
advance module will be discussed in Sections 2.1, 2.2, and 2.3 respectively. A
small example of the framework will be provided in Section 2.4. How this
framework can be used in statistical inference will be discussed in Section
2.5. We introduce the signal strength model that can be used for the
connection likelihood in Section 3. In Section 4 we illustrate our approach
with an fictional application. The implementation is described in Section 5.
We conclude with a couple of remarks in Section 6.

##  2 A modular Bayesian framework

In the proposed method, we will overlay the geographic area of interest with a
grid. The main advantage of using grid tiles is that different geospatial
datasets can be combined without the need to calculate spatial intersections,
which is a time consuming operation. Moreover, the mathematics described below
is easier since all grid tiles have the same area.

The key to the proposed localisation method is Bayes’ formula, which is used
in the following way:

$\centering\mathbb{P}(g\mid a)\propto\mathbb{P}(g)\mathbb{P}(a\mid
g),\@add@centering$ (1)

where $g$ represents a grid tile and $a$ a cell. The probability
$\mathbb{P}(g)$ that a device is located in grid tile $g$ without any
connection knowledge represents the location prior about the relative
frequency of events at grid tile $g$. The connection likelihood
$\mathbb{P}(a\mid g)$ is the probability that a device is connected to cell
$a$ given that the device is located in grid tile $g$. The location posterior
$\mathbb{P}(g\mid a)$ represents the probability that a device is located in
grid tile $g$ given that the device is connected to cell $a$.

As a technical aside, we note that the expression $\mathbb{P}(a\mid g)$ plays
two distinct roles in our discussion. When we model $\mathbb{P}(g\mid a)$ via
Eq. (1), the cell $a$ is considered fixed. Hence $\mathbb{P}(a\mid g)$ is not
a probability distribution, but a likelihood function with parameter $g$. On
the other hand, when we will propose a model for $\mathbb{P}(a\mid g)$ in
Section 2.2, it will be considered as a probability distribution over $a$ for
fixed $g$. Although this distinction is not critical for the reader’s
understanding of our model, the intended interpretation should be clear from
context. For simplicity, we will refer to $\mathbb{P}(a\mid g)$ as the
connection likelihood henceforth.

The position of Eq. (1) in our modular Bayesian framework is illustrated in
Fig. 1. The location prior and connection likelihood are estimates produced
from models we call the location prior module and connection module,
respectively. The former module may use land use data as input, while we
propose to use cell plan data. The location prior is then updated by the
connection likelihood to form the location posterior. The localisation method
can at this point be considered as complete and this posterior may be used for
further statistical inference. Alternatively, the posterior can be seen as a
new prior and updated again with other likelihoods if the necessary
information is available. Fig. 1 shows how a Timing Advance module might be
used to construct an updated location posterior.

The rest of this section is devoted to proposals for, and elaborations of
these modules.

Figure 1: A modular framework for modelling location posteriors.

###  2.1 Location prior

We define the _location prior_ $\mathbb{P}(g)$ as the probability that a
device is present in grid tile $g$, such that

$\sum_{g\in\mathcal{G}}\mathbb{P}(g)=1,$ (2)

where $\mathcal{G}$ is the set of all possible location estimates in our
model, in other words, the whole grid.

The definition of the location prior function will be based on assumptions
about where a device is expected to be. In this section, we propose four
options: the uniform prior, the land use prior, the network prior, and the
composite priors.

####  2.1.1 Uniform prior

When we use the _uniform prior_ , we assume the probability of a device being
in any grid tile is the same value for every grid tile:

$\mathbb{P}_{\text{uniform}}(g):=\frac{1}{|\mathcal{G}|}.$ (3)

A uniform prior is sometimes viewed as uninformative. In the case of mobile
phone data, however, the implicit assumption that any grid tile is as likely
as the next can lead to an underestimation of devices in urban areas and an
overestimation of devices in rural areas. We therefore advise against using
the uniform prior as a default prior without consciously assessing the
plausibility of the underlying assumption.

####  2.1.2 Land use prior

An alternative to the assumption of uniformity is to use administrative
sources on land use for the location prior. One would for example expect more
devices in an urban area than in a meadow. Hence, our second proposed prior is
called the _land use prior_ and it is based on a proportional expectation of
the number of devices $n(g)$ in grid tile $g$ such that:

$n(g)\propto\mathbb{E}[\text{number of devices in $g$}]\quad\text{for all
}g\in\mathcal{G}.$ (4)

The land use prior is then defined as:

$\mathbb{P}_{\text{land
use}}(g):=\frac{n(g)}{\sum_{g^{\prime}\in\mathcal{G}}n(g^{\prime})}.$ (5)

Due to the normalisation in its definition, the land use prior does require an
explicit estimate of the number of devices per grid tile. Any proportional
measure has the same effect. One way to utilise this is when information is
available on land use classes of the grid, such as levels of urbanisation. Let
there be $K$ land use classes, each with their own relative expected number of
devices: $u_{1},u_{2},\ldots,u_{K}$. Let
$w_{1}(g),w_{2}(g),\ldots,w_{K}(g)\in[0,1]$ be the proportion of the grid tile
that is covered by each class respectively, such that

$\sum_{\mathclap{k=1,\ldots,K}}w_{k}(g)=1\quad\text{for all }g\in\mathcal{G}.$
(6)

Then $n(g)$ can be modelled as

$n(g):=\sum_{\mathclap{k=1,\ldots,K}}u_{k}\cdot w_{k}(g).$ (7)

An example of a simple land use classification and the relative expected
numbers is shown in Table 1.

Table 1: An example of land use classes and their relative expected number of
devices.

Land use class | $u_{k}$
---|---
Urban | 1.0
Main roads | 0.5
Other land | 0.1
Water | 0.0

One of the downsides of the land use prior is that the assumptions based on
administrative sources are less flexible in the case of major events. A
festival in a location which would ordinarily be a quiet meadow, but suddenly
contains many devices, is not accounted for in the land use prior. Such events
can be recognized by the positioning and setup of the cells. For instance,
extra small cells are often used to compensate for large numbers of devices
[40, 39].

It can also be worthwhile to let the land use prior depend on the time and
day. For instance, the expected number of devices in industrial areas might be
smaller during nighttime and weekends compared to daytime and working days.

####  2.1.3 Network prior

The following prior, which we call the _network prior_ , is defined as

$\mathbb{P}(g):=\frac{\sum_{a\in\mathcal{A}}s(g,a)}{\sum_{a\in\mathcal{A}}\sum_{g^{\prime}\in\mathcal{G}}s(g^{\prime},a)},$
(8)

where $s(g,a)$ for a grid tile $g$ and cell $a$ is called _signal dominance_ ,
and $\mathcal{A}$ denotes the set of all cells in the network. In Section 2.2,
a general notion of _signal dominance_ will be introduced to serve as input
for our model for the connection likelihood. A specific instance of it, the
_Voronoi signal dominance_ $s_{\text{Vor}}(g,a)$, will be defined in that
section as well, and another instance $s_{\text{strength}}(g,a)$ will be
elaborated on in Section 3. Any of these instances can be used for Eq. (8).

This prior, together with our model for the connection likelihood
$\mathbb{P}(a\mid g)$, which is discussed in Section 2.2, simplifies Eq. (1)
to

$\mathbb{P}(g\mid a)\propto s(g,a).$ (9)

The interpretation of the prior (8) depends on the instance of the signal
dominance $s(g,a)$ one has chosen. Let $F\colon\mathcal{G}\to\mathcal{A}$ be
an arbitrary function that maps grid tiles to cells. Then, by setting

$s_{F}(g,a):=\begin{cases}1&\text{if $a=F(g)$,}\\\
0&\text{otherwise},\end{cases}$ (10)

the prior in Eq. (8) becomes uniform. Since $s_{\text{Vor}}$ is a specific
instance of the function $s_{F}$, the network prior with $s=s_{\text{Vor}}$
also simplifies to the uniform prior.

When using the term _network prior_ from here on, we refer to Eq. (8) with
signal dominance based on the instance $s_{\text{strength}}(g,a)$.

Basically, the network prior reflects the distribution of the total signal
over all the grid tiles. This prior contains implicit knowledge about where an
MNO is expecting people. The placement of cells is not without reason;
generally, more cells are placed in crowded areas, such as city centers, than
in quiet rural areas. Note that we could have defined the network prior using
the cell density. However, since the network capacity also depends on the type
and configuration of the cells and on the environment (buildings and trees
will generally have a negative effect on the propagation of the signal) we use
the signal dominance, through which these aspects are taken into account.

There are two aspects to be aware of when using the network prior. First, the
placement of cells is based on estimated peak traffic rather than the average
expected number of devices. MNOs normally provide better network coverage in
railway stations than in residential areas, since the estimated peak traffic
is higher; people typically use their phone more actively in railway stations
and moreover, the expected number of devices fluctuates more over time. The
second aspect to be aware of is that MNOs might place extra (partially
overlapping) cells in an area for reasons other than an expected increase of
the number of devices. They might do this, for example, when they detect that
the quality of the network connection is insufficient on certain sections of
land. In summary, the total signal strength of the network does not always
reflect the estimated number of devices.

From a Bayesian perspective, it may seem odd to use the same input, i.e. the
signal strength, in both the connection likelihood distribution and the
location prior. However, we use it in two different, complementary, ways. The
example in Section 2.4 includes calculations in which the signal strength
model is used for both the connection likelihood and location prior.

####  2.1.4 Composite priors

Our fourth and final proposed prior is less theoretically substantiated and
more driven by practical considerations. One can combine all three priors
described earlier as follows:

$\displaystyle\mathbb{P}_{\text{composite}}(g)$
$\displaystyle:=\pi_{\text{uniform}}\cdot\mathbb{P}_{\text{uniform}}(g)+{}$
(11) $\displaystyle\phantom{{}:={}}\pi_{\text{land
use}}\cdot\mathbb{P}_{\text{land use}}(g)+{}$
$\displaystyle\phantom{{}:={}}\pi_{\text{network}}\cdot\mathbb{P}_{\text{network}}(g),$

where $\pi:=(\pi_{\text{uniform}},\pi_{\text{land use}},\pi_{\text{network}})$
is any vector in $\mathbb{R}^{3}$ such that $0\leq\pi_{\text{uniform}}\leq 1$,
$0\leq\pi_{\text{land use}}\leq 1$, $0\leq\pi_{\text{network}}\leq 1$, and
$\pi_{\text{uniform}}+\pi_{\text{land use}}+\pi_{\text{network}}=1$.

The components of $\pi$ represent the contributions of the three previously
defined priors to the final _composite prior_
$\mathbb{P}_{\text{composite}}(g)$. Hence both the advantages and
disadvantages of all three priors are mixed. The optimal choice of $\pi$
depends on the situation at hand and must be derived experimentally. A
composite prior is harder to interpret from a theoretical point of view
compared to the three priors discussed so far, which can be seen as an
additional disadvantage.

###  2.2 Connection likelihood

We define the _connection likelihood_ $\mathbb{P}(a\mid g)$ for a cell $a$ and
a grid tile $g$ to be the probability that when a device located in grid tile
$g$ generates an event at some cell, it does so at $a$. We model this
probability as

$\mathbb{P}(a\mid
g):=\frac{s(g,a)}{\sum_{a^{\prime}\in\mathcal{A}}s(g,a^{\prime})},$ (12)

where $\mathcal{A}$ is the set of all cells in the MNOs network and
$s(g,a)\in[0,\infty)$ stands for the _signal dominance_ (an umbrella term
introduced by ourselves) received in grid tile $g$ from cell $a$. That is, the
connection likelihood is the ratio of the signal dominance received from cell
$a$ to the total value of signal dominance received from all cells. Different
choices for modelling the signal dominance are possible, and any choice
defines the _connection module_ in Fig. 1. Note that $\mathbb{P}(a\mid g)$ is
independent of rescaling the function $s(g,a)$ by a constant, and our
convention is that $s(g,a)$ should be defined so as to take on values in the
interval $[0,1]$.

One simple method to define the connection module is via Voronoi tessellation.
In this case, $s(g,a)$ is set to be

$s_{\text{Vor}}(g,a):=\begin{cases}1&\text{if }g\in\text{Vor}(a),\\\
0&\text{otherwise,}\end{cases}$ (13)

where $\text{Vor}(a)$ is the set of grid tiles of which the centroids lie in
the Voronoi region surrounding cell $a$. Denoting (12) in this case by
$\mathbb{P}_{\text{Vor}}(a\mid g)$ then gives

$\mathbb{P}_{\text{Vor}}(a\mid g)=s_{\text{Vor}}(g,a),$ (14)

and combining this connection likelihood with the uniform prior
$\mathbb{P}_{\text{uniform}}(g)$ results in the location posterior

$\mathbb{P}_{\text{Vor}}(g\mid
a)=\begin{cases}\bigl{\lvert}\\{g^{\prime}\in\mathcal{G}\mid
g^{\prime}\in\text{Vor}(a)\\}\bigr{\rvert}^{-1}&\text{if
$g\in\text{Vor}(a)$,}\\\ 0&\text{otherwise.}\end{cases}$ (15)

In Section 3 we propose a different, more advanced definition of the
connection module by first approximating the signal strength $S(g,a)$ measured
in $\mathrm{dBm}$, and then applying a transformation to it to obtain a signal
dominance $s_{\text{strength}}(g,a)$.

###  2.3 Incorporating Timing Advance data

Some MNOs include in their signalling data a so called _Timing Advance_
variable [1]. Timing advance is also known as round-trip time [45], or trip-
time bins [27]. The value of such a variable represents a time duration, and
it is used to estimate and adjust for communication delays. However, in
combination with the speed of radio waves it can alternatively be used to
estimate the distance between device and cell [20]. In the case of 4G
signalling data the Timing Advance variable takes on values in the discrete
set $\\{0,1,\ldots,1282\\}$. If an event then contains such a value $\tau$ for
this variable, the associated device is located approximately in an annulus
centered around the antenna of width $78\text{\,}\mathrm{m}$ and an inner
circle of radius $\tau\cdot$78\text{\,}\mathrm{m}$$.

One should be aware, though, of errors present in this approximation. As the
authors of [31] observed, the Timing Advance variable may take on different
values while the device has the same distance to the cell, and even when the
device is at the same location at different times.

Knowledge of $\tau$ may be used to improve the location posterior
$\mathbb{P}(g\mid a)$ through Bayesian updating:

$\mathbb{P}(g\mid a,\tau)\propto\mathbb{P}(g\mid a)\mathbb{P}(\tau\mid a,g),$
(16)

where the _Timing Advance likelihood_ $\mathbb{P}(\tau\mid a,g)$ can be
modelled as the fraction of the grid tile $g$ which lies in the annulus around
the cell $a$ specified by $\tau$. Since the maximum value of $\tau$ is $1282$,
this new location posterior $\mathbb{P}(g\mid a,\tau)$ equals $0$ for grid
tiles further than approximately $100\text{\,}\mathrm{km}$ from the cell.

Computing the likelihoods $\mathbb{P}(\tau\mid a,g)$ for all Timing Advance
annuli around all cells $a$ in the network and all tiles $g$ in the grid that
is used might prove to be too expensive computationally if calculated in the
way we suggest above. One could therefore model $\mathbb{P}(\tau\mid a,g)$
more coarsely as being $1$ if the centroid of $g$ lies in the annulus
specified by $\tau$, and $0$ otherwise. If the grid tiles used are
substantially larger than $78\text{\,}\mathrm{m}$, though, such as the $100$
by $100$ meter tiles we propose, this coarser Timing Advance likelihood model
could increase the location estimation error for devices located in the
smaller annuli. These errors can be mitigated by merging adjacent annuli. For
example, one could model

$\mathbb{P}(\tau\mid a,g):=\begin{cases}1&\text{if }g\in\text{Ann}(\tau,b),\\\
0&\text{otherwise,}\end{cases}$ (17)

where $\text{Ann}(\tau,b)$ specifies the grid tiles whose centroids lie in the
annuli specified by $\\{\tau-b,\ldots,\tau+b\\}$, where $b$ is a globally
defined integer, independent of $\tau$, $a$ and $g$, that determines how many
annuli are merged on both sides to the annulus corresponding to $\tau$.

Non-binary distributions can be constructed from Timing Advance data [27]. The
proposed connection likelihood relies on location of antenna, propagation
direction and a PRACH curve, which is the total number of connections per
month, per cell, and per Timing Advance distance band. PRACH curves can be
incorporated in our Timing Advance module if $\mathbb{P}(\tau\mid a,g)$ is
modelled as the observed ratio between devices that connect to Timing Advance
annulus $\tau$ and all devices connected to antenna $a$.

###  2.4 Example

To illustrate the computations involved in the model, we consider a small
fictional island of 1 by 3 kilometers. All numerical values mentioned in this
example can be found in Table 2. The island can be divided into three grid
tiles of equal size, $g_{1}$, $g_{2}$ and $g_{3}$. Note that for a more
realistic example we would use much smaller grid tiles, but for simplicity
each tile in this example measures 1 by 1 kilometer.

Figure 2: A schematic top view of an island of 1 by 3 kilometers.

There are two towns A and B on the island, of which the latter is about three
times as large as the former. Two cells, $a_{1}$ and $a_{2}$, have been
installed. These are illustrated in Fig. 2. Cell $a_{1}$ has perfect signal
dominance in $g_{1}$ and $g_{2}$, but no signal in $g_{3}$, while cell $a_{2}$
has perfect signal dominance in $g_{3}$ and $g_{2}$, but no signal in $g_{1}$.
Perfect signal dominance and no signal are expressed by values $1$ and $0$,
respectively. The signal dominance can be interpreted as the outcome of the
signal strength model from Section 3.

We calculate the four priors as defined in Section 2.1 based on the above
information. The connection likelihood for each cell is calculated according
to Eq. 12 based on the signal dominance values. Finally, the connection
likelihood and priors are combined to location posteriors for all combinations
of cells and priors.

Table 2: The corresponding numbers of the fictional example where the composite prior is based on $\pi:=(\pi_{\text{uniform}},\pi_{\text{land use}},\pi_{\text{network}})=(0,\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}})$. |  |  Grid tile $g$
---|---|---
|  | $g_{1}$ | $g_{2}$ | $g_{3}$
Signal dominance | $s(g,a_{1})$ | 1 | 1 | 0
| $s(g,a_{2})$ | 0 | 1 | 1
Location priors | $\mathbb{P}_{\text{uniform}}(g)$ | $\nicefrac{{1}}{{3}}$ | $\nicefrac{{1}}{{3}}$ | $\nicefrac{{1}}{{3}}$
| $\mathbb{P}_{\text{land use}}(g)$ | $\nicefrac{{1}}{{4}}$ | $\nicefrac{{3}}{{4}}$ | 0
| $\mathbb{P}_{\text{network}}(g)$ | $\nicefrac{{1}}{{4}}$ | $\nicefrac{{2}}{{4}}$ | $\nicefrac{{1}}{{4}}$
| $\mathbb{P}_{\text{composite}}(g)$ | $\nicefrac{{1}}{{4}}$ | $\nicefrac{{5}}{{8}}$ | $\nicefrac{{1}}{{8}}$
Connection likelihood | $\mathbb{P}(a_{1}\mid g)$ | 1 | $\nicefrac{{1}}{{2}}$ | 0
| $\mathbb{P}(a_{2}\mid g)$ | 0 | $\nicefrac{{1}}{{2}}$ | 1
Location posterior | $\mathbb{P}_{\text{uniform}}(g\mid a_{1})$ | $\nicefrac{{2}}{{3}}$ | $\nicefrac{{1}}{{3}}$ | 0
| $\mathbb{P}_{\text{uniform}}(g\mid a_{2})$ | 0 | $\nicefrac{{1}}{{3}}$ | $\nicefrac{{2}}{{3}}$
| $\mathbb{P}_{\text{land use}}(g\mid a_{1})$ | $\nicefrac{{2}}{{5}}$ | $\nicefrac{{3}}{{5}}$ | 0
| $\mathbb{P}_{\text{land use}}(g\mid a_{2})$ | 0 | 1 | 0
| $\mathbb{P}_{\text{network}}(g\mid a_{1})$ | $\nicefrac{{1}}{{2}}$ | $\nicefrac{{1}}{{2}}$ | 0
| $\mathbb{P}_{\text{network}}(g\mid a_{2})$ | 0 | $\nicefrac{{1}}{{2}}$ | $\nicefrac{{1}}{{2}}$
| $\mathbb{P}_{\text{composite}}(g\mid a_{1})$ | $\nicefrac{{4}}{{9}}$ | $\nicefrac{{5}}{{9}}$ | 0
| $\mathbb{P}_{\text{composite}}(g\mid a_{2})$ | 0 | $\nicefrac{{5}}{{7}}$ | $\nicefrac{{2}}{{7}}$

The location posterior allows for all kinds of further calculations such as
modelling the distribution of devices or even persons over grid tiles, when
connection events are counted for each cell during a time interval. This
process is complicated enough to be viewed as a separate research topic,
especially when one is interested in probability distributions rather than
point estimates. We consider it part of the _further statistical inference_ in
the framework from Fig. 1.

###  2.5 Statistical inference

The outcome of the modular system described in this paper is the location
posterior $\mathbb{P}(g\mid a)$, which specifies the probability that a device
is located in grid tile $g$, given that it is connected to cell $a$. This can
be used to calculate the total number of devices that are present at a
specific location during a specific time interval, or the number of devices
that move from one city to another. However, many applications in official
statistics are about numbers of people, for instance the number of visitors of
a tourist destination during holidays, or the number of people who commute
between two cities. Additional methods and auxiliary data are needed to
translate estimates of devices to estimates of people.

A generic framework has been proposed to organize the production process
needed for statistical inference on MNO data [33, 35]. According to this
framework, the production process runs through three distinct layers. The
bottom layer is called the data- or D-layer and consists of the processing of
raw mobile network data, which takes place at the MNO. The processing methods
that take place in this layer are dependent on the mobile network technology
used. The statistics- or S-layer is the top layer in which the processed
mobile phone data is used for statistical purposes. The convergence- or
C-layer connects these two layers with processing mobile network data sources
into data that can be used for statistical purposes. This intermediate layer
is needed since mobile network technology is complex and constantly changing.
The output of the C-layer should be a stable source for the S-layer, in which
this is used in combination with other data sources to produce statistics.

Our framework takes place in the D-layer, since mobile network data is
processed for constructing the connection likelihood. Note that it does not
matter which method is used for this process, since all described methods use
mobile network data, e.g. the Voronoi method uses cell tower locations. The
output of our framework, i.e. the location posterior, belongs in the C-layer,
since this does not depend on technology, and hence can be used directly for
statistical purposes. Using prior information could be theoretically be placed
in the S-layer. Note that the process should ideally be run at the MNO due to
potential privacy issues.

##  3 Signal strength model

This section describes the propagation of signal strength originating from a
single cell. We distinguish two types of cells: omnidirectional and
directional, resulting in two different propagation models. Omnidirectional
cells have no aimed beam and their coverage area can be thought of as a
circular disk. Directional cells point in a certain direction and their
coverage area can be thought of as an oval with one axis of symmetry. In
practice, small cells are usually omnidirectional and normal cells (i.e.
attached to cell towers or placed on rooftops) are often directional [19].

###  3.1 Omnidirectional cells

For omnidirectional cells, propagation of the signal strength $S(g,a)$ is
modelled as

$S(g,a):=S_{0}-S_{\text{dist}}(r_{g,a}),$ (18)

where $S_{0}$ is the signal strength at $r_{0}=1$ meter distance from the cell
in $\mathrm{dBm}$ and $r_{g,a}$ is the distance between the middle point of
grid tile $g$ and cell $a$ in meters. The value of $S_{0}$ can be different
for every cell and is assumed to be a known property. In cell plan data, it is
common to list the power $P$ of a cell in Watt, rather than the signal
strength in $\mathrm{dBm}$. The value of $S_{0}$ can be calculated from $P$
using the conversion between Watt and $\mathrm{dBm}$ [10]:

$S_{0}=30+10\log_{10}(P).$ (19)

The function $S_{\text{dist}}(r)$ returns the loss of signal strength as a
function of distance $r$:

$S_{\text{dist}}(r):=10\log_{10}(r^{\gamma})=10\gamma\log_{10}(r),$ (20)

where $\gamma$ is the path loss exponent, which resembles the reduction of
propagation due to reflection, diffraction and scattering caused by objects
such as buildings and trees [36]. In free space, $\gamma$ equals 2, but in
practice higher values should be used. As a rule of thumb, 4 can be used for
urban areas and 6 for indoor environments. Special situations, such as
tunnels, could improve the propagation such that a value of less than 2 is
applicable. The path loss exponent can be approximated by using the land use
register.

In Fig. 3, the signal loss as a function of the distance is shown for a cell
with $10\text{\,}\mathrm{W}$ power that is standing in an urban environment
($\gamma=4)$.

Figure 3: Signal loss as a function of the distance for a specific cell.

###  3.2 Directional cells

A directional cell is a cell that is aimed at a specific angle. Along this
angle, the signal strength is received at its best. However, the signal can
also be strong in other directions. It is comparable to a speaker producing
sound in a specific direction. The sound is audible in many directions, but is
much weaker at the sides and the back of the speaker. We specify the beam of a
directional cell $a$ by four parameters:

  * •

The azimuth angle $\varphi_{a}$ is the angle from the top view between the
north and the direction in which the cell is pointed, such that
$\varphi_{a}\in[0,360)$ degrees. Note that cell towers and rooftop cells often
contain three cells with 120 degrees in between.

  * •

The elevation angle $\theta_{a}$ is the angle between the horizon plane and
the tilt of the cell. Note that this angle is often very small, typically only
four degrees. The plane that is tilt along this angle is called the elevation
plane.

  * •

The horizontal beam width $\alpha_{a}$ specifies in which angular difference
from the azimuth angle in the elevation plane the signal loss is
$3\text{\,}\mathrm{dB}$ or less. At $3\text{\,}\mathrm{dB}$, the power of the
signal is halved. The angles in the elevation plane for which the signal loss
is $3\text{\,}\mathrm{dB}$ correspond to $\varphi_{a}\pm\alpha_{a}/2$. In
practice, these angles are around 65 degrees.

  * •

The vertical beam width $\beta_{a}$ specifies the angular difference from
$\theta_{a}$ in the vertical plane orthogonal to $\varphi_{a}$ in which the
signal loss is $3\text{\,}\mathrm{dB}$. The angles in which the signal loss is
$3\text{\,}\mathrm{dB}$ correspond to $\theta_{a}\pm\beta_{a}/2$. In practice,
these angles are around 9 degrees.

Let $\delta_{g,a}$ be the angle in the elevation plane between the azimuth
angle $\varphi_{a}$ and the orthogonal projection on the elevation plane of
the line between the center of cell $a$ and the center of grid tile $g$.
Similarly, let $\varepsilon_{g,a}$ be the angle from the side view between the
line along the elevation angle $\theta_{a}$ and the line between the center of
cell $a$ and the center of grid tile $g$. Note that $\varepsilon_{g,a}$
depends on the cell property of the installation height above ground level. We
model the signal strength for directional cells as

$S(g,a):=S_{0}-S_{\text{dist}}(r_{g,a})-S_{\text{azi}}(\delta_{g,a},\alpha_{a})-S_{\text{elev}}(\varepsilon_{g,a},\beta_{a}),$
(21)

where $S_{0}$ is the signal strength at $r_{0}=1$ meter distance from the
cell, in the direction of the beam so that $\delta=0$ and $\varepsilon=0$. The
signal loss due to distance to the cell, azimuth angle difference and
elevation angle difference is specified by $S_{\text{dist}}$, $S_{\text{azi}}$
and $S_{\text{elev}}$, respectively. The definition of $S_{\text{dist}}$ is
similar to the omnidirectional cell and can be found in Eq. (20).

Each cell type has its own signal strength pattern for both the azimuth and
elevation angles. These patterns define the relation between signal loss and
the offset angles, i.e., $\delta_{g,a}$ for the azimuth and
$\varepsilon_{g,a}$ for the elevation angles. We model the radiation pattern
for both $S_{\text{azi}}$ and $S_{elev}$ by a linear transformation of the
Gaussian formula, each with different values for parameters $c$ and $\sigma$.
Let

$f(\varphi)=c-ce^{-\frac{\varphi^{2}}{2\sigma^{2}}},$ (22)

where $c$ and $\sigma^{2}$ are constants, whose value is determined by
numerically solving equations for a set of constraints. These constraints are
different for $S_{\text{azi}}$ and $S_{\text{elev}}$ and depend on cell
properties.

The resulting patterns are shown in Fig. 4. The black line shows the relation
between signal loss and angle in the azimuth plane (left) and elevation plane
(right). The grey circles correspond to the signal loss; the outer circle
means $0\text{\,}\mathrm{dB}$ loss (which is only achieved in the main
direction), the next circle corresponds to $5\text{\,}\mathrm{dB}$ loss, an so
forth. The red lines denote the angles corresponding to
$3\text{\,}\mathrm{dB}$ loss. The angle between the red lines is $2\alpha_{a}$
in the azimuth plane and $2\beta_{a}$ in the elevation plane.

Figure 4: Radiation patterns for the azimuth and elevation planes.

Although these models approximate the general curve of real radiation
patterns, the radiation patterns are more complex in reality, e.g. they often
contain local spikes caused by so-called side and back lobes.

Fig. 5 (top row) illustrates the signal strength at the ground level from
above for a specific cell. In this case, the cell is placed at $x=0$, $y=0$ at
55 meters above ground level in an urban environment ($\gamma=4$), has a power
of $10\text{\,}\mathrm{W}$, and is directed eastwards with an elevation angle
(tilt) of 5 degrees, a horizontal beam width of 65 degrees and a vertical beam
width of 9 degrees. Notice that the signal strength close to the cell, which
on ground level translates to almost under the cell, is lower than at a couple
of hundred meters distance. This is caused by relatively large $\varepsilon$
angles at grid tiles nearby the cell.

Figure 5: Signal strength (top row) and signal dominance (bottom row) at ground level. Table 3: Indication of quality for signal strength in 4G networks. Signal strength ($\mathrm{dBm}$) | Quality
---|---
$-70$ or higher | Excellent
$-90$ to $-70$ | Good
$-100$ to $-90$ | Fair
$-110$ to $-100$ | Poor
$-110$ or less | Bad or no signal

###  3.3 Signal dominance

The assignment of a cell to a mobile device does not only depend on received
signal strength, but also on the capacity of the cells. The process of
assigning devices to cells while taking into account the capacity of the cells
is also called _load balancing_.

Our model allows for two phenomena that we feel should not be overlooked. The
first is the switching of a device when it is receiving a bad signal to a cell
with a better signal. Table 3 describes how the signal strength can be
interpreted in terms of quality for 4G networks [19]. The second phenomenon is
the switching between cells that is influenced by some decision making system
in the network that tries to optimize the load balancing within the network.
The specifics of this system are considered unknown.

We assume that a better signal leads to a higher chance of connection. When a
device has multiple cells available with a signal strength above a certain
threshold, say $-90\text{\,}\mathrm{dBm}$, the signal strengths are both more
than good enough and the cell with the highest capacity is selected rather
than the cell with the best signal strength. When the choice is between cells
with a lower signal strength, one can imagine that their relative differences
play a more important role in the connection process. However, when there are
multiple cells available with a poor signal strength, it can be assumed that
the signal strength value is less important than having capacity. In short, we
assume that signal strength plays a more important role in load balancing when
it is in the middle range instead of in the high quality or low quality
ranges.

To model this take on the load balancing mechanism, we use a logistic function
to translate the signal strength $S(g,a)$ to the more interpretable signal
dominance measure $s_{\text{strength}}(g,a)$, which is then used to define the
connection likelihood (12). Let us define

$s_{\text{strength}}(g,a):=\frac{1}{1+\exp\bigl{(}-S_{\text{steep}}\left(S(g,a)-S_{\text{mid}}\right)\bigr{)}},$
(23)

where $S_{\text{mid}}$ and $S_{\text{steep}}$ are parameters that define the
midpoint and the steepness of the curve respectively. Fig. 6 shows an example
of Eq. (23).

Figure 6: Logistic relation between signal strength ($\mathrm{dBm}$) and
signal dominance, where $S_{\text{mid}}$ and $S_{\text{steep}}$ are set to
$-92.5\text{\,}\mathrm{dBm}$ and $0.2\text{\,}\mathrm{dBm}$ respectively to
resemble Table 3.

The signal dominance at ground level is shown in Fig. 5 (bottom row). The
values that are shown are normalized by the sum of all values over all grid
tiles, such that the normalized values form a probability distribution.
Compared to the signal strength shown in Fig. 5 (top row), the signal
dominance puts more emphasis on the geographic area that is in the range of
the cell. Whether these signal dominance values resemble reality, should be
validated by field tests.

##  4 Application

A fictional example that resembles a real-world situation is illustrated in
Fig. 7(a), where a base map is shown with artificially placed cells. The
northwestern part of this area is urbanised, while the rest is mostly rural
with some small villages. The black dots represent directional macro cells,
which are deployed in triplets, with propagation angles about 120 degrees
apart from each other. The white dots represent omnidirectional small cells.

Figure 7: The application of our signal strength model.

The Voronoi map for this example, which we will use as a reference for our
signal strength model, is shown in Fig. 7(b). We apply the Voronoi algorithm
only to macro cells, where we shift the locations of the Voronoi points
towards the propagation directions. Subsequently, the grid tile areas in which
the small cells are located are subtracted from the Voronoi tessellation, and
assigned to the small cells.

Fig. 7(c) shows the _best server map_ based on our signal strength model. The
best server map is a tessellation where each geographic location is assigned
to the cell for which the signal strength is the strongest. This map provides
insights in the results of our signal strength model. Observe for instance
that some areas are relatively large and a few areas even consist of multiple
disjoint patches. This is caused by a combination of factors, for instance the
power, height, and elevation angle (tilt) of the cells, and the land
elevation.

Figs. 7(d) and (e) show the signal strength (in $\mathrm{dBm}$) and signal
dominance ($s$) of one specific cell respectively, which is colored red. The
propagation direction of this is south-west. Signal strength measurement tools
can be used to validate the modelled signal strength, and if needed, to adjust
the model parameters. The signal dominance is harder to validate, since it is
a model of the deployed load balancing mechanism. A deep analysis of MNO data
using ground truth location data, for instance GPS data, can be used for this
validation task.

Fig. 8 shows the results of our modular framework for the same region, shown
in (a). The rest of this figure is arranged as a cross table where the
location priors are shown as rows in (c), the connection likelihoods as
columns in (b), and the location posteriors as table cells in (d). For
instance, the bottom right map of this figure is the result of using the
composite prior (which consists of half the land use prior and half the
network prior) and the signal strength likelihood.

In the maps shown in Fig. 8(b), (c), and (d), sequential color palettes have
been applied where yellow corresponds to relatively low values and
respectively dark green, blue, and brown correspond to relatively high values.
For instance, the dark blue areas in the land use prior map correspond to the
areas with buildings and roads, and therefore, where devices are expected to
be. The blue areas in the network prior map highlight the areas that have a
high network coverage.

Two different connection likelihoods for the selected cell are shown in Fig.
8(b). Observe that the Voronoi likelihood has the same shape as the Voronoi
area of this cell shown in Fig. 7(b) in orange, which follows directly from
Eq. (14) which states that a likelihood value equals 1 if and only if the
corresponding grid tile is located in the Voronoi area. The high values for
the signal strength likelihood are located in the south/south-west direction
of the cell, whereas the propagation direction is south-west. This can be
explained by the fact that there is less overlap with other cells in the south
direction in comparison to the west direction.

The combination of the location priors and the connection likelihoods result
in the location posteriors shown in Fig. 7(d). The land use prior is seen to
have a strong effect on the posterior distributions. The network prior places
more weight on the areas with better network coverage. Note that, by applying
Eq. (9), the combination of network prior and signal strength likelihood
results in a posterior distribution that is a rescaled version of the signal
dominance, shown in Fig. 7(e). Finally, the combination of the composite prior
and signal strength likelihood takes into account both land use and network
coverage, and would therefore be a good candidate to be used for further
statistical inference.

Figure 8: Example showing how three location priors and two connection
likelihoods are combined into six location posteriors.

The final step in our modular framework is to update the location posterior by
Timing Advance. Recall from Section 2.3 that Timing Advance is a variable in
MNO data from which the distance between the device and the cell can be
estimated. Fig. 9 illustrates the location posterior shown in the bottom right
map of Fig. 8 (composite prior, signal strength likelihood) updated with
$\tau$ values of 15, 25, and 35. The Timing Advance likelihood defined by Eq.
(17) has been used with parameter $b$ set to 1, which means three annuli have
been merged, in Fig. 9(a) the annuli that correspond to $\tau$ values of 14,
15, and 16. In this application, the width of each annulus is 78.12 meter, so
this means that only grid tiles whose centroids are between 1094 to 1328 away
from the selected cell have a positive probability in the updated posterior.
For Fig. 9(b) and (c) those distance intervals are [1878, 2109] and [2656,
2890] respectively.

Figure 9: Updated location posterior with Timing Advance.

##  5 Implementation

The methods described in this paper have been implemented in mobloc [37] and
mobvis [38], packages for the programming language R. The former package is
used for the calculations, the latter for the graphical user interfaces and
visualizations. In this section we present a general work flow for mobloc and
illustrate how mobvis is used in this work flow. Details and reproducible
examples can be found in the documentation of the packages.

The general work flow of mobloc is depicted in Fig. 10. The first step in the
work flow is to collect all relevant datasets. The most important dataset is
is the cell plan. The geographical locations of the cells is required, whereas
other physical properties, such as height, propagation direction, and power,
are recommended in order to run the signal strength model. Other datasets that
can be used as input are elevation data, land use data, and administrative
region boundaries.

The collected data need to be in the correct format in order to be processed.
Information about which object classes are supported are described in the
package documentation [37]. The preprocessing stage is also used to make sure
that the spatial objects have one common Coordinate Reference System (CRS). It
is strongly recommended to use a CRS that preserves areas and distances for
the region of interest [26]. The function validate_cellplan is used to
validate the cell plan and to check whether it is consistent with the other
input data objects.

Figure 10: Work flow of mobloc.

The third step in the work flow is to run the signal strength model. This step
is only required if the network prior is used or if the signal strength model
is used in the connection likelihood module. The parameters for the signal
strength model can be configured with a graphical user interface started with
the function setup_sig_strength_model from the mobvis package. A screenshot of
this tool is shown in Fig. 11. The function compute_sig_strength from the
mobloc package is used to compute the signal strength and the corresponding
signal dominance.

Figure 11: Screenshot of the signal strength model setup tool. Table 4: An
overview of mobloc functions.

Module | Option | Function
---|---|---
Location prior | Uniform prior | create_uniform_prior
| Network prior | create_network_prior
| Auxiliary data (e.g. land use) prior | create_prior
Connection likelihood | Signal strength | create_strength_llh
| Voronoi | create_voronoi_llh
Bayesian update |  Posterior $\mathbb{P}(g\mid a)$, see Eq. (1) | calculate_posterior
|  Posterior $\mathbb{P}(g\mid a,\tau)$, updated | update_posterior_TA
|  with Timing Advance, see Eq. (16) |

The fourth and final step is to run the whole modular framework that is
depicted in Fig. 1. The functions that correspond to the modules are listed in
Table 4.

The mobvis package contain functions to visualize the results. All maps from
Figs. 7, 8, and 9 have been created with mobvis. This package also contains an
interactive tool, started with explore_mobloc, in which the results are
visualized. Via this tool, it is possible to select different cells, and
compare different module options. A screenshot is shown in Fig. 12, in which
the same fictional cell plan data is used as described in Section 4.

Figure 12: Screenshot of the mobloc exploration tool.

The output of mobloc is a data table with the following variables: cell
identification number, tile identification number, optionally a Timing Advance
value, and the corresponding posterior probability. This table can be used for
further statistical processing and analysis with MNO data [33, 35].

##  6 Concluding remarks

We have proposed a Bayesian approach to estimate the location of mobile
devices using mobile operator network data. The methods described in this
paper are modular, in the sense that one method can easily be replaced by
another. If, for instance, a better propagation model exists, which for
example uses a 3d model of the environment, this can be used together with the
other methods. The same applies for the signal dominance function, connection
likelihood, and location prior.

An MNO often facilitates mobile communication via multiple generations of
networks (e.g. 3G and 4G). For each generation, an MNO maintains a cellular
network. We currently assume that a mobile device only connects to cells from
one network generation, that is, the latest generation which the device
supports. When this assumption holds, the methods can be independently applied
for each generation. The networks can be viewed as independent because a
device will only connect to cells from one generation. In reality however,
this assumption might not hold for reasons such as coverage gaps, capacity and
network-specific optimal communication mode such as text, voice and internet.
More research on switching between cells from different generations is needed.

For each generation, cells will serve at different frequencies. For instance
an MNO may have a network of 4G cells that serve at 900, 1800 and
$2100\text{\,}\mathrm{Hz}$. As with handling of cells from multiple
generations, more research is needed on the process of switching between
frequencies.

Our propagation model can be calibrated using other data sources. Field
measurements of received signal strength could provide insights in several
parameters, for instance the power of the cells and the path loss exponent.
The calibration process should ideally be executed for each mobile phone
network, since they may be configured in different ways.

Alternative estimation methods that build on our methods have been proposed
recently. One method is the maximum likelihood estimator, which uses the
connection likelihood [34, 22]. Our signal strength model can be used here for
the connection likelihood, but this is not necessary. More recently, another
method has been proposed that uses the connection likelihood to select a range
of candidate estimates, and subsequently selects the estimate that best
confirms prior information [32]. Although the results of both methods are
promising when applied to fictitious examples, they rely more on the
connection likelihood than our method does. In scenarios where the connection
likelihood is of lesser quality due to incomplete cell plan data or where
prior data are available that have good predictive value, our location
posterior estimation method may be preferable. Further research is needed to
compare these methods by applying them with real data and validate them with
external sources.

MNOs often have additional data sources that can be used for validation, such
as coverage maps and best server maps, which are usually created with
propagation models and measurements. One could compare these maps to those
resulting from location estimation methods. However, be aware that these maps
are often created for other purposes, such as network analysis and
optimization. More research is needed on how to interpret the
(dis)similarities between such maps.

A more thorough validation approach is to compare logged location data from a
sample of devices, e.g. with GPS, and compare this to the location estimates
that are based on MNO data from the same sample of devices. This would require
consent both from the MNO as well as from the device owners.

The work of this paper has been focused solely on static situations; we
estimate the location posterior probabilities of a device given the cell at
which a certain event has been logged. However, when we consider the dynamic
behaviour of a device and take into account the sequence of logged events,
much more information can be extracted by interrelating the corresponding
location posterior probabilities, for instance by using a Hidden Markov Model
[35].

## References

  * [1] 3GPP.  TS Evolved Universal Terrestrial Radio Access (E-UTRA); Physical layer procedures, 3 2019.
  * [2] Ahas, R., Aasa, A., Yuan, Y., Raubal, M., Smoreda, Z., Liu, Y., Ziemlicki, C., Tiru, M., and Zook, M. Everyday space–time geographies: using mobile phone-based sensor data to monitor urban activity in Harbin, Paris, and Tallinn.  International Journal of Geographical Information Science 29, 11 (2015), 2017–2039.
  * [3] Alexander, L., Jiang, S., Murga, M., and González, M. C. Origin–destination trips by purpose and time of day inferred from mobile phone data.  Transportation Research Part C: Emerging Technologies 58 (2015), 240 – 250.  https://doi.org/10.1016/j.trc.2015.02.018.
  * [4] Badr, H., Du, H., Marshall, M., Dong, E., Squire, M., and Gardner, L. Association between mobility patterns and covid-19 transmission in the usa: a mathematical modelling study.  The Lancet Infectious Diseases 20, 11 (Nov. 2020).
  * [5] Calabrese, F., Ferrari, L., and Blondel, V. Urban sensing using mobile phone network data: A survey of research.  ACM Computing Surveys 47, 2, 1–20.
  * [6] De Meersman, F., Seynaeve, G., Debusschere, M., Lusyne, P., Dewitte, P., Baeyens, Y., Wirthmann, A., Demunter, C., Reis, F., and Reuter, H. Assessing the quality of mobile phone data as a source of statistics.  In European Conference on Quality in Official Statistics (2016), Eurostat.
  * [7] Deville, P., Linard, C., Martin, S., Gilbert, M., Stevens, F. R., Gaughan, A. E., Blondel, V. D., and Tatem, A. J. Dynamic population mapping using mobile phone data.  Proceedings of the National Academy of Sciences 111, 45 (2014), 15888–15893.  10.1073/pnas.1408439111.
  * [8] Diao, M., Zhu, Y., Joseph Ferreira, J., and Ratti, C. Inferring individual daily activities from mobile phone traces: A boston example.  Environment and Planning B: Planning and Design 43, 5 (2016), 920–940.  10.1177/0265813515600896.
  * [9] Dujardin, S., Jacques, D., Steele, J., and Linard, C. Mobile phone data for urban climate change adaptation: Reviewing applications, opportunities and key challenges.  Sustainability 12, 4 (02 2020), 1501.
  * [10] Figueiras, J. a., and Frattasi, S. Mobile Positioning and Tracking: From Conventional to Cooperative Techniques.  John Wiley and Sons, Ltd., 2010.
  * [11] Graells-Garrido, E., Peredo, O. F., and García, J. Sensing urban patterns with antenna mappings: The case of Santiago, Chile.  Sensors 16, 7 (2016), 1098.
  * [12] Grantz, K., Meredith, H., Cummings, D., Metcalf, C., Grenfell, B., Giles, J., Mehta, S., Solomon, S., Labrique, A., Kishore, N., Buckee, C., and Wesolowski, A. The use of mobile phone data to inform analysis of covid-19 pandemic epidemiology.  Nature Communications 11, 1 (Dec. 2020).
  * [13] Iqbal, M. S., Choudhury, C., Wang, P., and Gonzalez, M. C. Development of origin–destination matrices using mobile phone call data.  Transportation Research Part C: Emerging Technologies 40 (03 2014), 63–74.  10.1016/j.trc.2014.01.002.
  * [14] Järv, O., Tenkanen, H., and Toivonen, T. Enhancing spatial accuracy of mobile phone data using multi-temporal dasymetric interpolation.  International Journal of Geographical Information Science 31, 8 (2017), 1630–1651.
  * [15] Jiang, S., Yang, Y., Gupta, S., Veneziano, D., Athavale, S., and González, M. C. The timegeo modeling framework for urban motility without travel surveys.  Proceedings of the National Academy of Sciences 113, 37 (2016).
  * [16] Khalaf-Allah, M., and Kyamakya, K. Bayesian mobile location in cellular networks.  In 2006 14th European Signal Processing Conference (2006), pp. 1–5.
  * [17] Khalaf-Allah, M., and Kyamakya, K. Mobile location in gsm networks using database correlation with bayesian estimation.  In 11th IEEE Symposium on Computers and Communications (ISCC’06) (2006), pp. 289–293.  10.1109/ISCC.2006.103.
  * [18] Kondor, D., Grauwin, S., Kallus, Z., Gódor, I., Sobolevsky, S., and Ratti, C. Prediction limits of mobile phone activity modelling.  Royal Society open science 4, 2 (2017).
  * [19] Kora, A. D., Elono Ongbwa, B. A., Cances, J.-P., and Meghdadi, V. Accurate radio coverage assessment methods investigation for 3G/4G networks.  Computer Networks 107, P2 (Oct. 2016), 246–257.
  * [20] Kreher, R., and Gaenger, K. LTE Signaling, Troubleshooting and Optimization, 1 ed.  John Wiley and Sons, Ltd., 2011.
  * [21] Kung, K., Sobolevsky, S., and Ratti, C. Exploring universal patterns in human home-work commuting from mobile phone data.  PloS one 9 (11 2013).
  * [22] Laan, D. v. d., and Jonge, E. d. Maximum likelihood reconstruction of population densities from mobile signalling data.  In Proceedings of the NetMob 2019 Conference (2019).
  * [23] Lai, S., Erbach-Schoenberg, E., Pezzulo, C., Ruktanonchai, N., Sorichetta, A., Steele, J., Li, T., Dooley, C., and Tatem, A. Exploring the use of mobile phone data for national migration statistics.  Palgrave Communications 5, 34 (2019).
  * [24] Lu, X., Wrathall, D. J., Sundsøy, P. R., Nadiruzzaman, M., Wetter, E., Iqbal, A., Qureshi, T., Tatem, A., Canright, G., Engø-Monsen, K., and Bengtsson, L. Unveiling hidden migration and mobility patterns in climate stressed regions: A longitudinal study of six million anonymous mobile phone users in bangladesh.  Global Environmental Change 38 (2016), 1 – 7.  10.1016/j.gloenvcha.2016.02.002.
  * [25] M2Catalyst, LLC. Network Cell Info app, 2020.
  * [26] Maling, D. H. Coordinate Systems and Map Projections. Oxford: Pergamon Press., 1992.
  * [27] Ogulenko, A., Benenson, I., Omer, I., and Alon, B. Probabilistic positioning in mobile phone network and its consequences for the privacy of mobility data.  Computers Environment and Urban Systems 85 (10 2020), 101550.  10.1016/j.compenvurbsys.2020.101550.
  * [28] Oliver, N., Lepri, B., Sterly, H., Lambiotte, R., Deletaille, S., De Nadai, M., Letouzé, E., Salah, A. A., Benjamins, R., Cattuto, C., Colizza, V., de Cordes, N., Fraiberger, S. P., Koebe, T., Lehmann, S., Murillo, J., Pentland, A., Pham, P. N., Pivetta, F., Saramäki, J., Scarpino, S. V., Tizzoni, M., Verhulst, S., and Vinck, P. Mobile phone data for informing public health actions across the covid-19 pandemic life cycle.  Science Advances 6, 23 (2020).
  * [29] Panwar, N., Sharma, S., and Singh, A. K. A survey on 5g: The next generation of mobile communication.  Physical Communication 18 (2016), 64 – 84.  Special Issue on Radio Access Network Architectures and Resource Management for 5G.
  * [30] Pucci, P., Manfredini, F., and Tagliolato, P. Mapping urban practices through mobile phone data.  PoliMI SpringerBriefs Series, 02 2015.
  * [31] Raitoharju, M., Ali-Löytty, S., and Wirola, L. Estimation of base station position using timing advance measurements.  In Proceedings of SPIE – The International Society for Optical Engineering (12 2010), vol. 8285.
  * [32] Ricciato, F., and Coluccia, A. On the estimation of spatial density from mobile network operator data, 2020.
  * [33] Ricciato, F., Lanzieri, G., Wirthmann, A., and Seynaeve, G. Towards a methodological framework for estimating present population density from mobile network operator data.  Pervasive and Mobile Computing 6 (2020), 101263.
  * [34] Ricciato, F., Widhalm, P., Craglia, M., and Pantisano, F. Beyond the “single-operator, CDR-only" paradigm: An interoperable framework for mobile phone network data analyses and population density estimation.  Pervasive and Mobile Computing 35, February 2017 (2017), 65–82.
  * [35] Salgado, D., Sanguiao, L., Oancea, B., Barragán, S., and Necula, M. An end-to-end statistical process with mobile network data for official statistics, 2020.  Submitted to EPJ Data Science.
  * [36] Srinivasa, S., and Haenggi, M. Path loss exponent estimation in large wireless networks.  In 2009 Information Theory and Applications Workshop (Feb 2009), pp. 124–129.
  * [37] Tennekes, M. mobloc: Mobile phone location algorithms and tools, 2020.
  * [38] Tennekes, M. mobvis: Visualization of mobile phone location algorithm results, 2020.
  * [39] Tolstrup, M. Indoor Radio Planning: A Practical Guide for 2G, 3G and 4G. Wiley, 2015.
  * [40] Wang, S., Zhao, W., and Wang, C. Budgeted cell planning for cellular networks with small cells.  Vehicular Technology, IEEE Transactions on 64 (10 2015), 4797–4806.
  * [41] Widhalm, P., Yang, Y., Ulm, M., Athavale, S., and González, M. C. Discovering urban activity patterns in cell phone data.  Transportation 42, 4 (2015), 597–623.
  * [42] Wilson, R., zu Erbach-Schoenberg, E., Albert, M., Power, D., Tudge, S., Gonzalez, M., Guthrie, S., Chamberlain, H., Brooks, C., Hughes, C., Pitonakova, L., Buckee, C., Lu, X., Wetter, E., Tatem, A., and Bengtsson, L. Rapid and near real-time assessments of population displacement using mobile phone data following disasters: The 2015 Nepal earthquake.  PLOS Currents Disasters (Feb 2016).
  * [43] Xu, Y., Belyi, A., Bojic, I., and Ratti, C. Human mobility and socioeconomic status: Analysis of Singapore and Boston.  Computers, Environment and Urban Systems 72 (2018), 51–67.
  * [44] Zagatti, G. A., Gonzalez, M., Avner, P., Lozano-Gracia, N., Brooks, C. J., Albert, M., Gray, J., Antos, S. E., Burci, P., zu Erbach-Schoenberg, E., Tatem, A. J., Wetter, E., and Bengtsson, L. A trip to work: estimation of origin and destination of commuting patterns in the main metropolitan regions of Haiti using CDR.  Development Engineering 3 (2018), 133–165.
  * [45] Zang, H., Baccelli, F., and Bolot, J. Bayesian inference for localization in cellular networks.  In Proceedings of the 29th Conference on Information Communications (2010), INFOCOM’10, IEEE Press, p. 1963–1971.


# End-to-End Language Identification using Multi-Head Self-Attention and 1D
Convolutional Neural Networks

###### Abstract

In this work, we propose a new approach for language identification using
multi-head self-attention combined with raw waveform based 1D convolutional
neural networks for Indian languages. Our approach uses an encoder, multi-head
self-attention, and a statistics pooling layer. The encoder learns features
directly from raw waveforms using 1D convolution kernels and an LSTM layer.
The LSTM layer captures temporal information between the features extracted by
the 1D convolutional layer. The multi-head self-attention layer takes outputs
of the LSTM layer and applies self-attention mechanisms on these features with
M different heads. This process helps the model give more weightage to the
more useful features and less weightage to the less relevant features.
Finally, the frame-level features are combined using a statistics pooling
layer to extract the utterance-level feature vector label prediction. We
conduct all our experiments on the 373 hrs of audio data for eight different
Indian languages. Our experiments show that our approach outperforms the
baseline model by an absolute 3.69% improvement in F1-score and achieves the
best F1-score of 95.90%. Our approach also shows that using raw waveform
models gets a 1.7% improvement in performance compared to the models built
using handcrafted features.

Index Terms—  multi-head self-attention,language identification, 1D-CNNs

##  1 INTRODUCTION

A recent development in the area of deep neural networks has shown tremendous
improvements in speech systems, including speech recognition[1,2,3], emotion
recognition[20], speaker identification [13]. Previously, the language
identification field was dominated by i-vectors [5], which uses Gaussian
mixture models traditionally. Even today, i-vectors are considered to be the
best model in the case of low data regime. However, recent developments in the
field of deep learning show that deep neural networks are one of the
dominating approaches in language identification. Previously, W. Geng et al.
[4] proposed to use deep features extracted from a neural network trained for
speech recognition and showed that deep neural models are capable of obtaining
excellent performance over the classical systems[7,9,10,15,19]. Recently,
Time-delay neural networks(TDNN) have shown excellent performance for speech
recognition tasks [21]. X-vector [13] built for speaker identification was
used for language identification [14], and its shown to be one of the best
methods for language identification. The recent trend in sequence to sequence
mapping problems [18] involves the attention mechanism [16]. The attention
mechanism is one of the very well known techniques being used in sequence
mapping problems [17], and today’s state of the art speech recognition models
are built using the attention-based models. These models process sequential
inputs by iteratively selecting relevant features using attention technique.
Due to the efficiency of the attention technique for sequence mapping
problems, A. Vaswani et al. [22] propose to use stacks of attention layers
alone. They have shown remarkable results in machine translation. Motivated by
the work [22], this paper proposes to use multi-head self-attention in
combination with 1D convolutional neural network front-end processing for
language identification. Ours proposed model takes raw waveform as input
directly and extracts features that are useful for the LID task. The model
consists of a sequence of Residual blocks [23] of 1D convolutional layers to
extract features from the raw audio. Since the convolutional layer does not
have the capability to capture temporal information, we use an LSTM layer on
top of the convolutional network to capture time-varying information from the
input. Finally, the LSTM output feature sequence is fed into a multi-head
self-attention block consisting of multiple attention heads to select
important features from different parts of the input feature sequence using
attention weighting. Finally, an utterance level feature vector is generating
using a statistics pooling layer, and we classify the utterance level feature
to predict the class label.The organization of the paper is as follows. In
section 2, we explain our proposed approach in detail. In section 3, we give a
detailed analysis of the dataset collection and curation process, and in
section 4, we explain our experimental setup in detail. Finally, in section 5,
we describe our results.

##  2 Proposed method

In this section, we explain our proposed approach in detail. The detailed
model architecture is shown in Figure 1. Our model consists of 3 main stages,
1) An Encoder layer, which includes multiple 1D convolutional layers with
residual connections, an LSTM layer, 2) Multi-head self-attention layer to
select important features for language identification using attention
weighting, and 3)statistic pooling layer to obtain utterance level feature
vector for classification. The model takes raw audio waveform as input and
applies initial 1D convolution operations along with 1D max-pooling, as shown
in Figure 1. The initial convolutional layer features will go through a series
of there 1D Residual blocks followed by an LSTM layer. We then use multi-head
self-attention to extract relevant features from different parts of the input.
The statistics pooling layer generates an utterance level feature vector
containing language discriminative properties. The output of the statistics
pooling layer gives us a single feature vector called the utterance level
feature vector. This utterance level feature vector is fed into a projection
layer followed by a softmax layer to predict the class label. We explain the
details of each of these blocks in the following section.

Fig. 1: Proposed model architecture

###  2.1 Encoder

The encoder of our model consists of a series of three residual blocks
combined with a single LSTM layer, as shown in Figure 1. The encoder takes a
raw waveform signal and applies an initial 1D convolution operation consisting
of 64 filters of 1x7 filter size followed by a max-pooling operation. The max-
pooling is applied with a kernel size of 1x3 with stride 2. After the initial
convolution and max-pooling operation, we send the output a sequence of
residual blocks. The details of a single residual block is shown in Figure 1.
Each residual blocks operate with 1D convolution kernels of size 1x3. Residual
Block -1 consist of 2*64 convolution kernels of size 1x3, Residual Block -2
consist of 2*128 convolution kernels of size 1x3 and Residual Block -3
contains 2*256 convolution kernels of size 1x3. The output of Residual Block
-2 and Residual Block -3 will go through a 1D max-pooling operation. Since
LSTMs are known to capture long-range dependencies between frames in speech,
we use LSTM at the end of the final residual block to capture temporal
information. The output of the final residual block Residual Block -3 after
the max-pooling operation is sent to a single unidirectional LSTM layer with a
hidden size of 256. Let $\boldsymbol{X_{n}=[x_{1},x_{2}..x_{n},...x_{N}]}$ be
raw audio sequence with N samples.

$\boldsymbol{H^{A}}=\textsf{Encoder}(\boldsymbol{X_{n}})$ (1)

Where, Encoder is a mapping function which consists of initial 1D
convolutional layer and max-pooling operation, sequence of there residual
blocks Residual Block-1,Residual Block-2 and Residual Block-3 along with LSTM
layer . After this operation, we obtain a feature sequence
$\boldsymbol{H^{A}=[h_{1},h_{2}.....h_{T}]}$ of length T (T$<<$N). Typically
after the convolution operation, $\boldsymbol{H^{A}}$ can be looked at as a
feature matrix whose x-axis is a time dimension, and the y-axis is a feature
dimension. The feature dimension in our case is 256, as the hidden layer size
of the LSTM is the same.

###  2.2 Multi-Head Self-Attention

In this section, we describe multi-head self-attention in detail. It consists
of 3 different linear blocks, one for query, one for key, and another for
value. Each linear block consists of M independent linear layers. Where M is
the number of heads. The multi-head attention block takes features
$\boldsymbol{H^{A}=[h_{1},h_{2}.....h_{T}]}$ from LSTM and applies linear
transformation to create $\boldsymbol{Q_{i}}$, $\boldsymbol{K_{i}}$ and
$\boldsymbol{V_{i}}$ using $\boldsymbol{i^{\text{th}}}$ linear layers where,
$\boldsymbol{i=[1,2.....M]}$ and M is the total number of attention heads. The
$\boldsymbol{Q_{i}}$, $\boldsymbol{K_{i}}$ and $\boldsymbol{V_{i}}$ are fed
into scaled dot product attention layer. The scaled dot product attention
$\boldsymbol{A_{i}}$ for $\boldsymbol{i^{\text{th}}}$ head is defined as
follows.

$\boldsymbol{A_{i}}=\textsf{Softmax}(\boldsymbol{\frac{Q_{i}K_{i}}{d_{q}}})\boldsymbol{V_{i}}$
(2)

Where $\boldsymbol{d_{q}}$ is the dimension of the query vector. We combine
the attention output from all the heads using simple concatenation and feed
into the feed-forward layer.

$\boldsymbol{A}=\textsf{Concat}(\boldsymbol{A_{1},A_{2},A_{3}...A_{i}.....A_{M}}){W_{0}}$
(3)

Where, $\boldsymbol{A_{i}}$ is a $\boldsymbol{d_{q}}$x$\boldsymbol{T}$
dimensional matrix. Since the Concat operation is applied to the feature
dimension of all the matrices, the final output attention matrix
$\boldsymbol{A}$ from the multi-head attention block will have
$\boldsymbol{Md_{q}}$x$\boldsymbol{T}$ matrix dimensions.

The multi-head attention layer helps in finding features that are more
relevant for language identification. The scaled dot product attention
achieves this by giving more weighting to the more relevant features and less
weighting to less relevant features. Due to the presence of multiple heads in
the attention layer, this process selects features from different parts of the
input and helps in obtaining better language classification performance.

###  2.3 Statistics pooling

The idea of the statistics pooling layer is similar to max pooling. In the
case of statistics pooling, we compute the mean and standard deviation from
frame-level features. The mean and standard deviation features are
concatenated to create the utterance level feature vector, as described in the
equation below. Let $\boldsymbol{A=[a_{1},a_{2}.....a_{T}]}$ is the output
from multi-head attention block.

$\boldsymbol{P}=\textsf{Concat}({\textsf{mean}(\boldsymbol{A})},{\textsf{std}(\boldsymbol{A})})$
(4)

Where, $\boldsymbol{a_{i}}$ is a feature vector of dimension
$\boldsymbol{M*d_{q}}$ and $\boldsymbol{P}$ is final pooled feature vector
using statistics pooling layer. Since the dimension of the utterance level
feature vector is $\boldsymbol{P}$ become bigger when M is large, we add a
projection layer on top to the statistics pooling layer (Figure 1) in order to
reduce the dimension of $\boldsymbol{P}$. We take the output from this
projection layer to visualize the utterance level embeddings for different
languages.

##  3 Dataset

In this section, we describe our data collection process. We collect and
curate videos from Youtube111www.youtube.com using manual labeling. We ask
annotators to look for videos for eight languages in Youtube and manually
verify it make sure the video does not contain multiple Indian languages. Most
of these videos contain background noise or music signals. Sometime the video
may contain a mix of English and other Indian languages due to code-mixing. We
use an in-house speech v/s non-speech detection model to detect only the
speech segments. We clip only the speech segments from every video and discard
the non-speech part of the video. After preprocessing, our total dataset
contains 373.27 hrs of audio data for 8 Indian languages. Our dataset includes
Hindi, English, Kannada,Tamil, Telugu, Malayalam, Gujarati and Marathi. These
languages are officially spoken in the North and South regions of India. We
split the dataset into training and evaluation part, and the statistics of the
training and evaluation parts are shown in Table 1.

Datset Train Eval Duration Files Duration Files English 44.21 15963 11.28 4074
Kannada 35.95 12988 8.90 3216 Gujarati 30.24 10933 7.312 2642 Hindi 38.79
14004 9.72 3510 Malayalam 34.98 12636 8.94 3228 Tamil 63.15 22774 15.95 5753
Telugu 35.07 12666 8.54 3087 Marathi 16.23 5873 4.013 1449

Table 1: Train and evaluation splits for different languages (Duration is in
Hrs)

##  4 Experiments

We conduct all our experiments on in house dataset collected for 8 Indian
languages. Our proposed model consists of an encoder, multi-head self-
attention block, and statistics pooling layer followed by projection and a
softmax layer. We randomly select a 4sec audio signal from each audio file
during training. Since our data has a sampling rate of 16KHz, we get 64000
samples from every file during training. We feed a 1x64000 dimensional signal
into our encoder. We conduct multiple experiments to see the effectiveness of
the multi-head self-attention module for language identification. We first
train a standalone 1D convolutional neural network model as the first baseline
model. We refer this system as ResNet. We also train a 1D convolutional neural
network in combination with a unidirectional LSTM as a second baseline, and we
refer to it as ResNet-LSTM. Finally, our proposed model is built using a 1D
convolutional neural network, LSTM, and multi-head self-attention. We refer to
it as ResNet-LSTM-MHA-Raw.

We conduct multiple experiments to see the effectiveness of our model on the
duration of the audio during training. We train 3 different models ResNet-
LSTM-MHA-2Sec, ResNet-LSTM-MHA-3Sec, ResNet-LSTM-MHA-4Sec which takes
2sec,3sec and 4sec audio data respectively during training. Our final
experiments study the effectiveness of using raw waveform methods instead of
handcrafted features. We set up an experiment to train the model using MFCC
features as inputs instead of the raw waveform to our model. We extract 13
dimensional MFCC (with delta and double-delta) feature for every 25ms using a
10ms frameshift for this experiment. The MFCC based model is referred to as
ResNet-LSTM-MHA-MFCC while the raw waveform based model is referred to as
ResNet-LSTM-MHA-Raw in this paper. We use Adam [24] optimizer to train all our
models with a learning rate of 0.001 for up to 25 epochs. We use a batch size
of 64 during training. We train all our models using
Pytorch222https://pytorch.org/ toolkit.

##  5 RESULTS

In this section, we describe the evaluation of different models and their
performances. We train 2 baseline models ResNet and ResNet-LSTM. The first
baseline model ResNet consists fo a sequence of 3 Residual blocks made up of
1D convolution kernels. We can think of this model as a ResNet [23] with an
average pooling layer replaced by a statistics pooling layer. This model takes
4sec raw audio data and predicts the language label. The ResNet model has an
F1-score of 88.67% on the test dataset. The second baseline model ResNet-LSTM
consists of the same setting as baseline-1, but it has an extra LSTM layer on
top of CNN in order to capture long-range temporal information. The
performance of this model is 92.21%, as shown in Table 2. We compare our
baseline models with our proposed model ResNet-LSTM-MHA-RAW, which contains a
multi-head attention layer and operates with the raw waveform as input. Table
2 shows that our model gets 3.69% absolute improvement in F1-score compared to
the second baseline model. We also create a model that takes MFCC features as
input instead of raw audio refer to as ResNet-LSTM-MHA-MFCC. We show that raw
waveform based models can get 1.7% improvement over handcrafted feature-based
models.

System F1-Score ResNet (baseline-1) 89.67% ResNet-LSTM (baseline-2) 92.21%
ResNet-LSTM-MHA-MFCC (ours) 94.22% ResNet-LSTM-MHA-RAW (ours) 95.90%

Table 2: Comparison of different architectures for language identification.
Bold indicates the best performance

In order to see the effect of the input length during training, we conduct an
experiment to train the model using 2Sec, 3Sec, and 4Sec audio data, and we
refer to these models as ResNet-LSTM-MHA-2Sec, ResNet-LSTM-MHA-3Sec and
ResNet-LSTM-MHA-4Sec respectively. The results of these experiments are shown
in Table 3. It shows that longer audio data tends to improve the F1-score on
the test data due to longer context signals.

Table 3: Comparison of models trained with different segment duration. Bold
indicates the best performance

System F1-Score ResNet-LSTM-MHA-2Sec 92.64% ResNet-LSTM-MHA-3Sec 94.40%
ResNet-LSTM-MHA-4Sec 95.90%

Finally, we visualize the utterance level embeddings extracted from the
projection layer for all the languages. We extract embeddings for 6500
randomly selected test utterances for t-sne visualization. Each embedding has
dimension. We reduce the dimension of the embeddings to 2 using the t-sne
technique. The t-sne plot of the 2-D embeddings is shown in Figure 2. It can
be clearly seen that the proposed model learns very good language
discriminative features at the segment level.

Fig. 2: t-sne plot of utterence level embeddings

##  6 Conclusion

In this work, we propose a new architecture for language identification using
multi-head self-attention and 1D convolutional neural networks. We propose to
use raw waveform directly as input instead of handcrafted features to learn
language discriminative feature using 1D convolution operations. Our model
uses multi-head self-attention to learn and select more important features for
language identification task. We finally use a statistics pooling approach to
extract utterance level language representation from frame-level features. We
collect and curate 373hrs audio data for 8 Indian languages Hindi, English,
Kannada,Tamil, Telugu, Malayalam, Gujarati and Marathi. Our experiments show
that multi-head self-attention in combination with raw waveform based 1D
convolutional neural network model obtains the best performance on our
evaluation dataset. We extract the utterance level embeddings for our
evaluation data and visualize the clustering effect using t-sne. The
visualization clearly shows that the model learns very good language
discriminative features.

## References

  * [1] W. Chan, N. Jaitly, Q. V. Le, and O. Vinyals, “Listen, Attend and Spell: A Neural Network for Large Vocabulary Conversational Speech Recognition” in ICASSP, 2016
  * [2] S. Kim, T. Hori, and S. Watanabe, “Joint CTC-attention based endto-end speech recognition using multi-task learning,” in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2017, pp. 4835–4839.
  * [3] Ronan Collobert, Christian Puhrsch, and Gabriel Synnaeve, “Wav2letter: an end-to-end convnet based speech recognition system,” CoRR, vol.abs/1609.03193, 2016.
  * [4] W. Geng, J. Li, S. Zhang, X. Cai, and B. Xu, “Multilingual tandem bottleneck feature for language identification,” in Sixteenth Annual Conference of the International Speech Communication Association, 2015\.
  * [5] N. Dehak, P. Kenny, R. Dehak, P. Dumouchel, and P. Ouellet, “Front-end factor analysis for speaker verification,” IEEE Transactions on Audio, Speech,and Language Processing vol. 19, no. 4,pp. 788–798, 2011.
  * [6] Y. Song, B. Jiang, Y. Bao, S. Wei, and L.-R. Dai, “I-vector representation based on bottleneck features for language identification,” Electronics Letters, vol. 49, no. 24, pp. 1569–1570, 2013
  * [7] I. Lopez-Moreno, J. Gonzalez-Dominguez, O. Plchot, D. Martinez, J. Gonzalez-Rodriguez, and P. Moreno, “Automatic language identification using deep neural networks”, in Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on. IEEE, 2014, pp. 5337–5341
  * [8] S. Hochreiter and J. Schmidhuber, “Long short term memory”, Neural computation, 1997.
  * [9] J. Gonzalez-Dominguez, I. Lopez-Moreno, H. Sak, J. Gonzalez Rodriguez, and P. J. Moreno, “Automatic language identification using long short-term memory recurrent neural networks.” in INTERSPEECH, 2014, pp. 2155–2159.
  * [10] S. Ganapathy, K. J. Han, S. Thomas, M. K. Omar, M. Van Segbroeck, and S. S. Narayanan, “Robust language identification using convolutional neural network features.” ISCA INTERSPEECH, 2014\.
  * [11] Lozano Diez, A. Zazo Candil, R. Gonzalez Dominguez, J., Toledano, D.T., and Gonzalez Rodrıguez, J. (2015), “An end-to-end approach to language identification in short utterances using convolutional neural networks”, INTERSPEECH 2015
  * [12] Garcia-Romero, D., McCree, A. “Stacked Long-Term TDNN for Spoken Language Recognition,” Proc. Interspeech 2016, 3226- 3230\.
  * [13] D. Snyder, D. Garcia-Romero, G. Sell, D. Povey, and S. Khudanpur,“X-vectors: Robust dnn embeddings for speaker recognition,” in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2018
  * [14] D. Snyder, D. Garcia-Romero, A. McCree, G. Sell, D. Povey, and S. Khudanpur, “Spoken language recognition using x-vectors,” in Odyssey: The Speaker and Language Recognition Workshop, 2018
  * [15] C Bartz, T Herold, H Yang, and C Meinel ,“Language identification using deep convolutional recurrent neural networks”, CoRR,abs/1708.04811, 2017
  * [16] D. Bahdanau, K. Cho, and Y. Bengio. “Neural machine translation by jointly learning to align and translate”, arXiv preprint arXiv:1409.0473, 2014
  * [17] J. Chorowski, D. Bahdanau, K. Cho, and Y. Bengio, “End-to-end continuous speech recognition using attention-based recurrent nn: first results,” arXiv preprint arXiv:1412.1602, 2014
  * [18] I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to sequence learning with neural networks,” in Advances in neural information processing systems, 2014, pp. 3104–3112.
  * [19] M. Ravanelli and Y.Bengio, “Speaker Recognition from raw waveform with SincNet,” in Proc. of SLT, 2018
  * [20] M. Sarma, P. Ghahremani, D. Povey, N. K. Goel, K. K. Sarma, and N. Dehak, “Emotion identication from raw speech signals using dnns”, Proc. Interspeech 2018, pp. 3097–3101, 2018.
  * [21] V. Peddinti, D. Povey, and S. Khudanpur, “A time delay neural network architecture for efficient modeling of long temporal contexts,” in Proceedings of INTERSPEECH, 2015.
  * [22] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in Advances in Neural Information Processing Systems, 2017, pp. 5998–6008
  * [23] K. He, X. Zhang, S. Ren, and J. Sun. “Deep residual learning for image recognition,” Computer Vision and Pattern Recognition (CVPR), 2016
  * [24] Diederik P. Kingma and Jimmy Ba, “Adam: A Method for Stochastic Optimization”, In Proceedings of the International Conference on Learning Representations (ICLR), 2014
  * [25] Okabe, Koji and Koshinaka, Takafumi and Shinoda, Koichi,“Attentive Statistics Pooling for Deep Speaker Embedding”, In Interspeech 2018