# An interactive city choice model and its application for measuring the
intercity interaction

Xiang-Yu Jia  Key Laboratory of Transport Industry of Big Data Application
Technologies for Comprehensive Transport, Ministry of Transport, Beijing
Jiaotong University, Beijing 100044, China  Er-Jian Liu<EMAIL_ADDRESS>Institute of Transportation System Science and Engineering, Beijing Jiaotong
University, Beijing 100044, China  Chun-Yan Chen  Traffic Control Technology
Co.,Ltd, Beijing 100875, China  Xiao-Yong Yan<EMAIL_ADDRESS>Institute of
Transportation System Science and Engineering, Beijing Jiaotong University,
Beijing 100044, China

###### Abstract

Measuring the interaction between cities is an important research topic in
many disciplines, such as sociology, geography, economics and transportation
science. The traditional and most widely used spatial interaction model is the
gravity model, but it requires the parameters to be artificially set. In this
paper, we propose a parameter-free interactive city choice (ICC) model that
measures intercity interaction from the perspective of individual choice
behavior. The ICC model assumes that the probability of an individual choosing
to interact with a city is proportional to the number of opportunities in the
destination city and inversely proportional to the number of intervening
opportunities between the origin city and the destination city, calculated
using the travel time in the transportation network. The intercity interaction
intensity can be obtained by calculating the product of this probability and
the origin city’s population. We apply the ICC model to measure the
interaction intensity among 339 cities in China and analyze the impact of
changes in the Chinese land transportation network from 2005 to 2018 on the
intercity and city interaction intensity. Compared with the previously widely
used spatial interaction models, the measurement results of the ICC model are
more consistent with the actual situation.

## Introduction

The rapid development of cities worldwide and the acceleration of urbanization
have led to more than half of the world’s population living in
citiesRitchieRoser-482 ,and thus cities have become the main location for
human activities in today’s societyBarthelemy-550 . The connections between
cities through the transportation network promotes the flow of people, goods,
information, money and skills among cities; such flow between cities is called
intercity interactionOrd-552 ; Besag-549 . Understanding and predicting
intercity interaction patterns has long been an important research topic in
sociology, geography, economics, transportation science and many other
disciplinesWilson-560 ; najem2018debye . It also has great significance in the
rational formulation of urban development strategiesBatty-486 ;
paldino2015urban , the promotion of regional sustainable
developmentGatelyHutyra-556 , communicable disease
controlHufnagelBrockmann-490 ; EubankGuclu-491 ; BalcanColizza-492  and other
fields. As the intercity interaction intensity increases, cities are no longer
regarded as isolated individuals but as interdependent urban systemsBatten-553
. Therefore, understanding the intercity interaction and establishing a model
that can accurately measure the interaction between cities are of great value
for optimizing the spatial structure of urban agglomerationsBarthelemy-554 ;
LeNEChet-555 .

The gravity model was the first model proposed to measure intercity
interactionStewart-557 . The model assumes that the intensity of interaction
between two cities is proportional to the product of their sizes (e.g.,
population, GDP) and inversely proportional to a power law function of their
distance. The gravity model is simple in form and is widely used to predict
intercity interactions, such as intercity travelJungWang-499 , commuting
tripsViboudBjORnstad-501 ; lenormand2016systematic , population
migrationTobler-487  and international tradeFagiolo-567 . However, this model
is based on analogy with Newton’s law of universal gravitation and does not
involve individual spatial choice behaviorWilson-560 ; YanZhou-508 .
Furthermore, the parameter of the gravity model’s power law distance function
is artificially defined. For example, some researchers set the parameter to
1GoldenbergLevy-564 ; backstrom2010find ; rahman2009australia , while others
set it to 2KringsCalabrese-534 ; Keum-577 ; AlamUddin-578 . Therefore, it
would be a valuable contribution to establish a parameter-free model to
measure intercity interactions from the perspective of individual spatial
choice behavior.

Simini et al. took an important step forward in spatial interaction modeling
by establishing a parameter-free model named the radiation
modelSiminiGonzALez-509  to predict commuting trips between counties in the
U.S. This model assumes that the individual will consider the employment
opportunities provided by the work location and the benefits that the
opportunities may bring to him/her when choosing a work location. He/she will
choose the work location nearest to his/her home that offers a benefit greater
than the best offer available in his/her home county. Some researchers improve
the radiation model and propose various commuting prediction models, such as
the radiation model with selectionsimini2013human  and the flow and jump
modelvarga2018commuting . Recently, many researchers have applied the
radiation model or improved radiation models to measure intercity interaction
intensityLiFeng-561 ; Tian-562 ; ZhengKuang-563 . However, the radiation model
assumes that the individual will only choose the nearest location with a
higher benefit than his/her home, which reflects a cautious tendency of
individual choice behavior. It can predict commuting trips but is not suitable
for predicting general travelLiuYan-566  because travelers may choose not only
the closest location with a higher benefit than the origin but also other
locations with higher benefits than the origin and intervening destinations .
To solve this problem, Liu and Yan proposed another parameter-free model,
named the opportunity priority selection (OPS) modelLiuYan-519 , that adopts
the perspective of individual destination choice behavior. The OPS model
assumes that when the individual chooses a destination, he/she will choose a
location with a higher benefit than the benefit of the origin and the benefits
of the intervening opportunitiesStouffer-498 . This reflects an exploratory
tendency in individual choice behavior and can accurately predict human
mobility within and between cities. Compared with the radiation model, the OPS
model can better describe individual destination choice behavior between
cities, which implies that the OPS model is more suitable for measuring the
intercity interaction intensity. However, applications of the OPS model to
measure intercity interactions is still lacking.

In this paper, we establish an intercity interaction measurement model named
the interactive city choice (ICC) model by improving the OPS model. We apply
this model to measure the intercity interaction intensity in China and analyze
the impact of the change in China’s land transportation network from 2005 to
2018 on the city interaction intensity. Compared with the earlier gravity
model and radiation model, the results obtained by our model are more
consistent with reality, implying that it captures the underlying mechanism of
individual behavior in choosing to interact with cities.

## Interactive city choice (ICC) model

The OPS modelLiuYan-519  assumes that when an individual chooses a
destination, similar to the classic radiation modelSiminiGonzALez-509  and the
population-weighted opportunities modelYanZhao-512 ; YanWang-515 , he/she
first evaluates the benefit of the opportunities in each location, in which
the number of opportunities in a location is proportional to the location’s
population, and the benefit of opportunities is a random variable with a
distribution of $p(z)$ ($p(z)$ can be any continuous distribution). After
evaluating the benefit, the individual will select a location that presents
higher benefits than the origin and any intervening opportunities. According
to the above assumption, when an individual at location $i$ makes a choice for
location $j$, the probability of location $j$ being selected is

$Q_{ij}=\int_{0}^{\infty}{\mathrm{Pr}_{m_{i}+s_{ij}}(z)}\mathrm{Pr}_{m_{j}}(>z)\mathrm{d}z,$
(1)

where $m_{i}$ is the number of opportunities at location $i$, and $s_{ij}$ is
the number of intervening opportunities (i.e., the sum of the number of
opportunities at all locations at a shorter distance to $i$ than
$j$Stouffer-498 ; see Fig. 1(a)). ${\mathrm{Pr}_{m_{i}+s_{ij}}(z)}$ is the
probability that the maximum benefit obtained after $m_{i}+s_{ij}$ sampling is
exactly $z$, and ${\mathrm{Pr}_{m_{j}}(>z)}$ is the probability that the
maximum benefit obtained after $m_{j}$ samplings is greater than z.

Since $\mathrm{Pr}_{x}(<z)=p(<z)^{x}$, we can obtain

$\mathrm{Pr}_{x}(z)=\frac{\mathrm{d}\mathrm{Pr}_{x}(<z)}{\mathrm{d}z}=xp(<z)^{x-1}\frac{\mathrm{d}p(<z)}{\mathrm{d}z}.$
(2)

Substituting Eq. (2) in Eq.(1), we obtain

$\displaystyle Q_{ij}$
$\displaystyle=\int_{0}^{\infty}(m_{i}+s_{ij})p(<z)^{m_{i}+s_{ij}-1}\frac{\mathrm{d}p(<z)}{\mathrm{d}z}[1-p(<z)^{m_{j}}]\mathrm{d}z$
(3)
$\displaystyle=(m_{i}+s_{ij})\int_{0}^{1}{[p(<z)^{s_{ij}+m_{i}-1}-p(<z)^{m_{j}+s_{ij}+m_{i}-1}]}\mathrm{d}p(<z)$
$\displaystyle=(m_{i}+s_{ij})(\frac{p(<z)^{s_{ij}+m_{i}}}{s_{ij}+m_{i}}\Big{|}_{0}^{1}-\frac{p(<z)^{m_{j}+s_{ij}+m_{i}}}{m_{j}+s_{ij}+m_{i}}\Big{|}_{0}^{1})$
$\displaystyle=(m_{i}+s_{ij})(\frac{1}{s_{ij}+m_{i}}-\frac{1}{m_{j}+s_{ij}+m_{i}})$
$\displaystyle=\frac{m_{j}}{m_{i}+s_{ij}+m_{j}}.$

Figure 1: Schematic diagram of the calculation of intervening opportunities.
Each dot represents a location. Location $i$ (blue dot) is the origin, and
location $j$ (green dot) is the destination. (a) In the OPS model, the
intervening opportunity $s_{ij}$ refers to the sum of the number of
opportunities at all locations (i.e., locations in the orange circle, except
location $i$ and location $j$) whose distance $d$ from location $i$ is less
than the distance $d_{ij}$ from location $i$ to location $j$. (b) In the ICC
model, intervening opportunity$s_{ij}$ refers to the sum of the number of
opportunities at all locations (i.e., locations in the orange area, except
location $i$ and location $j$) whose travel time $t$ from location $i$ is less
than the travel time $t_{ij}$ from location $i$ to location $j$, where the
travel time is obtained by calculating the shortest travel time path through
the transportation network.

From Eq. (3), we can see that the OPS model can calculate the probability of
an individual choosing a destination without any adjustable parameters.
However, there are two problems in applying the OPS model to measure intercity
interactions. One is that the intervening opportunities $s_{ij}$ in the OPS
model are calculated by the geographic distance between two locations, as
shown in Fig. 1(a). However, in reality, locations are connected by a
transportation network. Most individuals compare which locations are easier to
reach by taking travel time, instead of geographic distance, as the main
factor. Therefore, the intervening opportunities $s_{ij}$ should be calculated
by the travel time between two locationsRenErcsey-Ravasz-516 , as shown in
Fig. 1 (b). The other is that the OPS model assumes that the number of
location opportunities is proportional to its population. However, the number
of opportunities provided by each city is not directly proportional to the
population but is more related to the city’s industrial scale, GDP or other
economic indicatorsReillyOthers-568 , among which the most commonly used
indicator is GDPFagiolo-567 . Therefore, it is more reasonable to use GDP to
reflect the number of opportunities. To solve these two problems, we establish
a probability model for individuals to choose to interact with a city. We
assume that the intervening opportunities $s_{ij}$ are calculated by the
travel time and that the number of location opportunities is proportional to
the city’s GDP. Given the GDP of each city, we can calculate the probability
that an individual located in city $i$ selects city $j$ as

$Q_{ij}=\frac{m_{j}}{m_{i}+s_{ij}+m_{j}},$ (4)

where $m_{i}$ is the GDP of city $i$, and $s_{ij}$ is the sum of the GDP of
all cities whose travel time from city $i$ is shorter than that of city $j$
(see the orange area in Fig. 1 (b)). Furthermore, if we know the total
population $n_{i}$ of city $i$, we can calculate the interaction intensity
from city $i$ to city $j$ as

$T_{ij}={n_{i}}{Q_{ij}}=\frac{{n_{i}}{m_{j}}}{m_{i}+s_{ij}+m_{j}}.$ (5)

We name Eq. (5) the interactive city choice (ICC) model. It should be noted
that the spatial interaction intensity $T_{ij}$ is not an actual flow volume
but a dimensionless value. Furthermore, according to the spatial interaction
intensity $T_{ij}$, we can calculate the interaction intensity of the city as

$\displaystyle
T_{i}=\frac{\sum\limits_{{j}\neq{i}}T_{ij}+\sum\limits_{{k}\neq{i}}T_{ki}}{2},$
(6)

where $\sum\limits_{{j}\neq{i}}T_{ij}$ is the sum of the outgoing interaction
intensity and $\sum\limits_{{k}\neq{i}}T_{ki}$ is that of the incoming
interaction intensitySimYaliraki-513 .

## Application of the ICC model to measuring intercity interaction intensity

In this section, we apply the ICC model to measure the interaction intensity
between cities in China and analyze the impact of the change in the Chinese
land transportation network from 2005 to 2018 on city interaction intensity.
It should be noted that China started to build high-speed railways in 2005;
thus, we select 2005 as the starting year. Because we can only download
Chinese economic and demographic data up to 2018, when we started this work,
we select 2018 as the ending year.

### Data and processing methods

We select 339 Chinese cities, including 333 prefecture-level cities, four
municipalities (Beijing, Tianjin, Shanghai, Chongqing) and two special
administrative regions (Hong Kong and Macao), as the research objects. We
download the population and GDP data of the 339 Chinese cities in 2005 and
2018 from the official website of the National Bureau of Statistics of China
and the data of the cities’ central points, Chinese road networks and railway
networks in 2018 from the OpenStreetMap website. The reason for selecting the
road and railway network data is that the total annual transportation volume
of these two land transportation modes accounts for more than 84% of the total
annual transportation volume of all intercity transportation modes (including
railway, road, waterway and airway) in both 2005 and 2018, as shown in Table
1. In these two types of data, roads include highways, national roads,
provincial roads, county roads and township roads; railways include high-speed
railways and normal railways. Since travel between any two cities can be
realized through national roads, we select the national road data as the basic
land transportation network data. Furthermore, we add three other types of
data (i.e., highways, normal railways and high-speed railways) that are
designed to provide faster travel than national roads in the land
transportation network. We establish the 2018 land transportation network, in
which the edges represent highways, national roads, normal railways or high-
speed railways. We add the city central point to the land transportation
network by connecting it to the nearest road within the urban area. We also
connect it to the nearest railway if there is a railway station within the
urban area. We also need to assign the travel time value to each edge to
calculate the intercity travel time in the land transportation network. We
know the length of each edge in the land transportation network, so we only
need to set the speeds of these four transportation modes (i.e., national
road, highway, normal railway and high-speed railway) to calculate the travel
time. According to the standards, including Code for Design of Railway Line
(TB 10098-2017) and Design Specification for Highway Alignment (JTG D20-2017),
the design speed range of high-speed railway is from 250 km/h to 350 km/h, of
the normal railway is from 80 km/h to 200 km/h, of the highway is from 80 km/h
to 120 km/h and of the national road is from 60 km/h to 100 km/h. For
simplicity, we use the median value of the speed range, i.e., 300 km/h as the
assumed speed for the high-speed railway, 140 km/h for the normal railway, 100
km/h for the highway and 80 km/h for the national road. We then calculate the
travel time of each edge by dividing its length by its assumed speed. We
obtain the 2005 land transportation network by deleting roads and railways
built after 2005 in the 2018 land transportation network according to the 2005
Chinese road map and 2005 Chinese railway map, as shown in Fig. 2.

Table 1: The proportion of annual passengers and freight transportation volume for various transportation modes in China Transportation mode |  2005 Passenger transportation volume(%) |  2005 freight transportation volume( %) |  2018 Passenger transportation volume(%) |  2018 freight transportation volume(%)
---|---|---|---|---
railway |  6.26 |  14.7 |  18.81 |  7.83
road |  91.9 |  72.35 |  76.17 |  76.73
waterway |  1.1 |  11.49 |  1.56 |  13.58
airway |  0.75 |  0.02 |  3.4 |  0.01
Figure 2: Chinese land transportation network and selected important cities.
(a) 2005. (b) 2018. Each circle represents a provincial capital city or
municipality. The size and color of the circle represent the population and
GDP of the city, respectively. Solid lines represent high-speed railways,
dashed-dotted lines represent normal railways, dashed lines represent highways
and dotted lines represent national roads.

### Calculation of the intercity interaction intensity

Figure 3: Distribution of interaction intensity among Chinese cities. (a)
2005. (b) 2018. The black dashed line is the Heihe-Tengchong line. The
thickness and color of the other lines indicate the intercity interaction
intensity.

We apply the ICC model to calculate the intercity interaction intensity in
2005 and 2018. We first calculate the travel time between cities by finding
the shortest intercity travel time path in both the 2005 and 2018 land
transportation networks. According to the intercity travel time, we can obtain
$s_{ij}$ by summing the GDP of all cities whose travel time from city $i$ is
less than the travel time from city $i$ to city $j$. We then calculate the
intercity interaction intensity in 2005 and 2018 according to Eq. (5). The
results are shown in Fig. 3, from which we can see that the intercity
interaction intensity in the east of the Heihe-Tengchong lineHuanyong-586  is
higher than that in the west in both 2005 and 2018. Considering Fig. 2 and
Fig. 3 comprehensively, we can see that the interaction intensity between
large cities in 2018 is significantly higher than that in 2005. This increase
is due to the more intensive construction of high-speed railways and highways
between large cities during these 13 years. Additionally, the intercity travel
time has been greatly shortened with the development of the land
transportation network. We can see from Fig. 4 that the proportion of short-
time interactions increases and the proportion of long-time interactions
decreases from 2005 to 2018. The longest travel time was shortened from 54 h
in 2005 to 43 h in 2018. Furthermore, as shown in the subgraph of Fig. 4, it
only takes 7.5 h to fall within 95% of the total national interaction
intensity in the land transportation network in 2018, while it took 13 h to
reach that in 2005.

Figure 4: Travel time distribution. The square dot and circular dot in the
figure represent the proportion of the intercity interaction intensity with a
travel time of $t$ hours in 2005 and 2018, respectively. The solid line and
dotted lines in the subgraph represent the cumulative probability distribution
of travel time in 2005 and 2018, respectively.

### Calculation of the city interaction intensity

Figure 5: City interaction intensity and the difference between city
interaction intensity in 2005 and that in 2018. (a) 2005 and (b) 2018 Chinese
city interaction intensity. The dotted line is the Heihe-Tengchong line and
each circle represents a city. The color and size of the circle indicate the
city interaction intensity. (c) The difference between city interaction
intensity in 2005 and that in 2018. Each circle represents a city and the size
of the circle represents the absolute value of the difference value.

We calculate the interaction intensity of Chinese cities in 2005 and 2018
according to Eq. (6). The results are shown in Fig. 5(a-b), from which we can
see that the two years show a similar distribution of city interaction
intensity, i.e., the cities with high interaction intensity are mainly
concentrated east of the Chinese Heihe-Tengchong Line. This is mainly because
these cities have a more developed economy, more opportunities and a more
intensive surrounding land transportation network, as shown in Fig. 2.
Furthermore, we calculate the difference in city interaction intensity between
2005 and 2018, as shown in Fig. 5 (c). We can see that the interaction
intensity of many cities, e.g., Wuhan, significantly improves, while the
interaction intensity of some cities, e.g., Hong Kong, decreases. The reasons
for this phenomenon are the changes in these cities’ GDP and the development
of land transportation networks (especially high-speed railway networks) from
2005 to 2018. For example, Wuhan’s GDP increased from 223.823 billion yuan in
2005 to 1484.729 billion yuan in 2018, with a growth rate of 663.35%. In
addition, the two national high-speed railway arteries (i.e., Beijing-
Guangzhou and Shanghai-Hanrong high-speed railways) established after 2005
both pass through Wuhan. The development of Wuhan’s GDP and transportation
infrastructure have led to a rapid increase in its attractiveness for
interaction, so the interaction intensity of Wuhan significantly improves. In
contrast, Hong Kong’s GDP increased from 1384.5 billion yuan in 2005 to
2400.098 billion yuan in 2018, with a growth rate of 73.35%. It has the lowest
GDP growth rate compared with other cities, which reduces its attraction for
interaction. In addition, we can see from Fig. 5(c) that the interaction
intensity of most cities along the high-speed railway increases, while that of
cities far from the high-speed railway generally decreases. This is mainly
because the travel time between cities along the high-speed railway and other
cities has been significantly reduced with the construction and rapid
development of the high-speed railway. Cities along the high-speed railway
will be chosen with higher probability by those choosing to interact with a
city.

Figure 6: Ranking of city interaction intensity. (a) 2005. (b) 2018. The
horizontal axis represents the ranking of cities, and the vertical axis
represents the city interaction intensity. Figure 7: The distribution of the
intercity interaction intensity calculated by different models. (a) Radiation
model. (b) Gravity model. The thickness and color of the lines indicate the
intercity interaction intensity.

We further rank the 339 cities according to their interaction intensity in
2005 and 2018. The results are shown in Fig. 6, from which we can see that the
ranking of city interaction intensity changes greatly from 2005 to 2018. For
example, Guangzhou’s ranking rises from fourth place to third place,
Shenzhen’s ranking rises from ninth place to seventh place, and Hong Kong’s
ranking drops from third place to fifteenth place. Although these three cities
are geographically close and their 2018 GDPs are similar (as shown in Table
2), they differ greatly in their city interaction intensity due to their
different positions in the land transportation network. Guangzhou is one of
three national comprehensive transportation hubs with three national road
arteries (i.e., Beijing-Guangzhou Line, Guiyang-Guangzhou Line and Nan-
Guangzhou Line). Shenzhen is also a comprehensive transportation hub that
connects Hong Kong, Macao and mainland China.

Table 2: Ranking of population, GDP and city interaction intensity calculated by the ICC model, radiation model and gravity model Rank |  Population |  GDP |  ICC model |  Radiation model |  Gravity model($\alpha$=1) |  Gravity model($\alpha$=2)
---|---|---|---|---|---|---
1 |  Shanghai |  Shanghai |  Beijing |  Beijing |  Shanghai |  Shenzhen
2 |  Beijing |  Beijing |  Shanghai |  Guangzhou |  Beijing |  Hong Kong
3 |  Chengdu |  Shenzhen |  Guangzhou |  Zhengzhou |  Shenzhen |  Guangzhou
4 |  Tianjin |  Hong Kong |  Wuhan |  Tianjin |  Guangzhou |  Foshan
5 |  Guangzhou |  Guangzhou |  Tianjin |  Chengdu |  Tianjin |  Shanghai
6 |  Shenzhen |  Tianjin |  Chengdu |  Shanghai |  Hong Kong |  Beijing
7 |  Baoding |  Suzhou |  Zhengzhou |  Wuhan |  Suzhou |  Hangzhou
8 |  Wuhan |  Chengdu |  Shenzhen |  Hangzhou |  Wuhan |  Tianjin
9 |  Shijiazhuang |  Wuhan |  Nanjing |  Xi’an |  Hangzhou |  Suzhou
10 |  Suzhou |  Hangzhou |  Changsha |  Shenzhen |  Nanjing |  Nanjing
11 |  Linyi |  Nanjing |  Xi’an |  Changsha |  Zhengzhou |  Xi’an
12 |  Zhengzhou |  Qingdao |  Suzhou |  Shijiazhuang |  Chengdu |  Zhengzhou
13 |  Nanyang |  Wuxi |  Hangzhou |  Nanjing |  Foshan |  Dongguan
14 |  Xi’an |  Changsha |  Xuzhou |  Xuzhou |  Changsha |  Wuhan
15 |  Hangzhou |  Ningbo |  Hong Kong |  Shenyang |  Wuxi |  Shaoxing
16 |  Handan |  Zhengzhou |  Shijiazhuang |  Suzhou |  Hefei |  Wuxi
17 |  Harbin |  Foshan |  Hefei |  Harbin |  Xuzhou |  Xianyang
18 |  Qingdao |  Quanzhou |  Chungking |  Jinan |  Ningbo |  Chengdu
19 |  Weifang |  Nantong |  Qingdao |  Changchun |  Shijiazhuang |  Jiaxing
20 |  Wenzhou |  Chungking |  Foshan |  Handan |  Dongguan |  Changsha

We next apply the radiation modelSiminiGonzALez-509  and gravity
modelJiaoTang-580  to calculate the intercity interaction intensity (see Fig.
7) and the interaction intensity of each city in 2018. We further list the
ranking of the interaction intensity of the top 20 Chinese cities calculated
by the ICC model, radiation model and gravity model (see Table 2). Compared
with the radiation model and gravity model, the result of the ICC model is
more consistent with the actual situation. This is because the radiation model
assumes that when the individual chooses an interactive city, he/she only
chooses the nearest city with a higher GDP than his/her homeSiminiGonzALez-509
; the gravity model assumes that the intercity interaction intensity is
proportional to the product of the two cities’ size (population or
GDP)JiaoTang-580  and inversely proportional to a power function of the travel
time with a power exponent of 1 or 2Odlyzko-581 . Neither the radiation model
nor the gravity model can accurately capture the mechanism of individuals
choosing to interact with cities, so they cannot reasonably reflect the city
interaction intensity. Unlike the radiation model and gravity model, the ICC
model’s ranking of the cities’ interactive intensity suggests that our model
can better describe the essential mechanism behind the individual choice to
interact with a city.

## Conclusion and discussion

The measurement of intercity interaction has important research significance.
In this paper, we develop the ICC model from the perspective of individual
choice behavior, which assumes that the probability of an individual choosing
to interact with a city is proportional to the number of opportunities, as
expressed by the GDP of the destination city, and inversely proportional to
the number of intervening opportunities, calculated by the shortest travel
time in the land transportation network. Multiplying this probability by the
origin city’s population, one can obtain the intercity interaction intensity.
To demonstrate the advantage of the ICC model, we apply the ICC model to
measure the interaction intensity among 339 cities in China. After collecting
and processing the big data related to intercity interaction, we analyze the
impact of the change in the land transportation network from 2005 to 2018 on
the intercity and city interaction intensity. We find that the travel time
between cities has decreased and the interaction intensity between large
cities has increased due to the development of land transportation. In
particular, the interaction intensity of cities along high-speed railways has
greatly increased. Compared with the previously widely used spatial
interaction models (i.e., the gravity model and radiation model), the results
obtained by the ICC model are more consistent with the actual situation.

The proposed ICC model not only helps us better measure the intercity
interaction intensity but also offers potential additional applications. For
example, the ICC model provides a new perspective for identifying suburbs,
which is a hot topic in geographical research. The traditional suburban
identification method usually refers to population density and the nature of
residential landTian-562 . The ICC model introduces the spatial interaction
intensity between city districts, which can improve the method of suburban
identification. In addition, the ICC model can calculate the interaction
intensity within and between urban agglomerations, providing valuable
indicators for a comprehensive evaluation of the degree of urban
agglomerationLiFeng-561 ; ShiYan-572 , which is of great significance for
urban agglomeration sustainable development.

Although the ICC model can obtain more realistic results when measuring
intercity interaction intensity, it still has room for expansion in practical
applications. In this paper, we use GDP, which is a key factor affecting the
number of opportunities, to reflect the number of opportunities. In reality,
there are many other factors, e.g., urban population, industrial size and
industrial structure, that also affect a city’s opportunities. Therefore, we
can use multiple factors to calculate the number of opportunities in future
applications. In addition, we only use the travel time calculated by the
shortest time path algorithm in the land transportation networks, including
roads and railways, to measure the interaction intensity among 339 Chinese
cities. However, the importance of various transportation modes is different
in different countries or regions. For example, airways are an important mode
of passenger transportation between the U.S. citieszhang2017energy , and
waterways are the main mode of freight transport between European
citiesbenga2019assesment . Therefore, future research can consider extending
the land transportation network to a more comprehensive three-dimensional
transportation network including roads, railways, airways and waterways to
make more reasonable measurements of intercity interactions.

## Acknowledgements

We thank Professor Kai Liu and her postgraduate student Wei-Hua Zhu at Beijing
Normal University for their help with our early research.

## Funding

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 71822102, 71621001, 71871010).

## References

  * (1) Ritchie, H., Roser, M.: Urbanization. Our world in data (2019). https://ourworldindata.org/urbanization
  * (2) Barthelemy, M.: The statistical physics of cities. Nat. Rev. Phys. 1(6), 406–415 (2019)
  * (3) Ord, K.: Estimation methods for models of spatial interaction. J. Am. Stat. Assoc. 70(349), 120–126 (1975)
  * (4) Besag, J.: Spatial interaction and the statistical analysis of lattice systems. J. R. Stat. Soc. Ser. B-Stat. Methodol. 36(2), 192–225 (1974)
  * (5) Wilson, A.G.: A family of spatial interaction models, and associated developments. Environ. Plan. A 3(1), 1–32 (1971)
  * (6) Najem, S., Faour, G.: Debye–hückel theory for refugees’ migration. EPJ Data Sci. 7(1), 1–9 (2018)
  * (7) Batty, M.: The size, scale, and shape of cities. Science 319(5864), 769–771 (2008)
  * (8) Paldino, S., Bojic, I., Sobolevsky, S., Ratti, C., González, M.C.: Urban magnetism through the lens of geo-tagged photography. EPJ Data Sci. 4(1), 5 (2015)
  * (9) Gately, C.K., Hutyra, L.R., Wing, I.S.: Cities, traffic, and co2: A multidecadal assessment of trends, drivers, and scaling relationships. Proc. Natl. Acad. Sci. U.S.A. 112(16), 4999–5004 (2015)
  * (10) Hufnagel, L., Brockmann, D., Geisel, T.: Forecast and control of epidemics in a globalized world. Proc. Natl. Acad. Sci. U.S.A. 101(42), 15124–15129 (2004)
  * (11) Eubank, S., Guclu, H., Kumar, V.A., Marathe, M.V., Srinivasan, A., Toroczkai, Z., Wang, N.: Modelling disease outbreaks in realistic urban social networks. Nature 429(6988), 180 (2004)
  * (12) Balcan, D., Colizza, V., Gon C C Alves, B., Hu, H., Ramasco, J.E.J., Vespignani, A.: Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. U.S.A. 106(51), 21484–21489 (2009)
  * (13) Batten, D.F.: Network cities: creative urban agglomerations for the 21st century. Urban Stud. 32(2), 313–327 (1995)
  * (14) Barthelemy, M.: The Structure and Dynamics of Cities. Cambridge University Press, Cambridge (2016)
  * (15) Le N E Chet, F.: Urban spatial structure, daily mobility and energy consumption: a study of 34 european cities. Cybergeo (2012)
  * (16) Stewart, J.Q.: An inverse distance variation for certain social influences. Science 93(2404), 89–90 (1941)
  * (17) Jung, W.-S., Wang, F., Stanley, H.E.: Gravity model in the korean highway. Europhys. Lett. 81(4), 48005 (2008)
  * (18) Viboud, C.E.C., Bj O Rnstad, O.N., Smith, D.L., Simonsen, L., Miller, M.A., Grenfell, B.T.: Synchrony, waves, and spatial hierarchies in the spread of influenza. Science 312(5772), 447–451 (2006)
  * (19) Lenormand, M., Bassolas, A., Ramasco, J.J.: Systematic comparison of trip distribution laws and models. J. Transp. Geogr. 51, 158–169 (2016)
  * (20) Tobler, W.: Migration: Ravenstein, thornthwaite, and beyond. Urban Geogr. 16(4), 327–343 (1995)
  * (21) Fagiolo, G.: The international-trade network: gravity equations and topological properties. J. Econ. Interact. Coord. 5(1), 1–25 (2010)
  * (22) Yan, X.-Y., Zhou, T.: Destination choice game: A spatial interaction theory on human mobility. Sci. Rep. 9(1), 9466 (2019)
  * (23) Goldenberg, J., Levy, M.: Distance is not dead: Social interaction and geographical distance in the internet era. arXiv:0906.3202 (2009)
  * (24) Backstrom, L., Sun, E., Marlow, C.: Find me if you can: improving geographical prediction with social and spatial proximity. In: Proceedings of the 19th International Conference on World Wide Web, pp. 61–70 (2010)
  * (25) Rahman, M.M.: Australia’s global trade potential: evidence from the gravity model analysis. In: Proceedings of the 2009 Oxford Business and Economics Conference (OBEC 2009), pp. 1–41 (2009)
  * (26) Krings, G., Calabrese, F., Ratti, C., Blondel, V.D.: Urban gravity: a model for inter-city telecommunication flows. J. Stat. Mech. Theory Exp. 2009(07), 07003 (2009)
  * (27) Keum, K.: Tourism flows and trade theory: a panel data analysis with the gravity model. Ann. Reg. Sci. 44(3), 541–557 (2010)
  * (28) Alam, M., Uddin, G., Taufique, K., Others: Import inflows of bangladesh: the gravity model approach. Int. J. Econ. Finance 1(1), 131–139 (2009)
  * (29) Simini, F., Gonz A Lez, M.C., Maritan, A., Barab A Si, A.-L.A.S.O.: A universal model for mobility and migration patterns. Nature 484(7392), 96 (2012)
  * (30) Simini, F., Maritan, A., Néda, Z.: Human mobility in a continuum approach. PLoS One 8(3), 60069 (2013)
  * (31) Varga, L., Tóth, G., Néda, Z.: Commuting patterns: the flow and jump model and supporting data. EPJ Data Sci. 7(1), 37 (2018)
  * (32) Li, F., Feng, Z., Li, P., You, Z.: Measuring directional urban spatial interaction in china: A migration perspective. PLoS One 12(1), 0171107 (2017)
  * (33) Tian, Y.: Mapping suburbs based on spatial interactions and effect analysis on ecological landscape change: A case study of jiangsu province from 1998 to 2018, eastern china. Land 9(5), 159 (2020)
  * (34) Zheng, W., Kuang, A., Wang, X., Chen, J.: Measuring network configuration of the yangtze river middle reaches urban agglomeration: Based on modified radiation model. Chin. Geogr. Sci. 30(4), 677–694 (2020)
  * (35) Liu, E.-J., Yan, X.-Y.: A universal opportunity model for human mobility. Sci. Rep. 10(1), 1–10 (2020)
  * (36) Liu, E., Yan, X.: New parameter-free mobility model: Opportunity priority selection model. Physica A 526, 121023 (2019)
  * (37) Stouffer, S.A.: Intervening opportunities: a theory relating mobility and distance. Am. Sociol. Rev. 5(6), 845–867 (1940)
  * (38) Yan, X.-Y., Zhao, C., Fan, Y., Di, Z., Wang, W.-X.: Universal predictability of mobility patterns in cities. J. R. Soc. Interface 11(100), 20140834 (2014)
  * (39) Yan, X.-Y., Wang, W.-X., Gao, Z.-Y., Lai, Y.-C.: Universal model of individual and population mobility on diverse spatial scales. Nat. Commun. 8(1), 1639 (2017)
  * (40) Ren, Y., Ercsey-Ravasz, M.A.R., Wang, P., Gonz A Lez, M.C., Toroczkai, Z.A.N.: Predicting commuter flows in spatial networks using a radiation model based on temporal ranges. Nat. Commun. 5, 5347 (2014)
  * (41) Reilly, W.J., Others: Methods for the study of retail relationships (1929)
  * (42) Sim, A., Yaliraki, S.N., Barahona, M., Stumpf, M.P.: Great cities look small. J. R. Soc. Interface 12(109), 20150315 (2015)
  * (43) Huanyong, H.: Distribution of china’s population: Accompanying charts and density map(in chinese). Acta Geogr. Sin. 2(2), 33–74 (1935)
  * (44) Jiao, L.M., Tang, X., Liu, X.P.: Spatial linkage and urban expansion: An urban agglomeration perspective. Prog. Geogr. 35(10), 1177–1185 (2016)
  * (45) Odlyzko, A.: The forgotten discovery of gravity models and the inefficiency of early railway networks. Œconomia (5-2) (2015)
  * (46) Shi, Q., Yan, X., Jia, B., Gao, Z.: Freight data-driven research on evaluation indexes for urban agglomeration development degree. Sustainability 12(11), 4589 (2020)
  * (47) Zhang, J., Fang, H., Wang, H., Jia, M., Wu, J., Fang, S.: Energy efficiency of airlines and its influencing factors: A comparison between china and the united states. Resour. Conserv. Recycl. 125, 1–8 (2017)
  * (48) Benga, G.C., Savu, D., Savu, S.V., Olei, A., Iacobici, R.I.: Assesment of trends in inland waterway transport within european union. In: Advanced Engineering Forum, pp. 247–254 (2019)


# LiDAR guided Small obstacle Segmentation

Aasheesh Singh∗1, Aditya Kamireddypalli∗1, Vineet Gandhi1 and K Madhava
Krishna1 * Equal Contribution1 KCIS, International Institute of Information
Technology Hyderabad, India<EMAIL_ADDRESS>
###### Abstract

Detecting small obstacles on the road is critical for autonomous driving. In
this paper, we present a method to reliably detect such obstacles through a
multi-modal framework of sparse LiDAR(VLP-16) and Monocular vision. LiDAR is
employed to provide additional context in the form of confidence maps to
monocular segmentation networks. We show significant performance gains when
the context is fed as an additional input to monocular semantic segmentation
frameworks. We further present a new semantic segmentation dataset to the
community, comprising of over 3000 image frames with corresponding LiDAR
observations. The images come with pixel-wise annotations of three classes
off-road, road, and small obstacle. We stress that precise calibration between
LiDAR and camera is crucial for this task and thus propose a novel Hausdorff
distance based calibration refinement method over extrinsic parameters. As a
first benchmark over this dataset, we report our results with 73 % instance
detection up to a distance of 50 meters on challenging scenarios.
Qualitatively by showcasing accurate segmentation of obstacles less than 15
cms at 50m depth and quantitatively through favourable comparisons vis a vis
prior art, we vindicate the method’s efficacy. Our project and dataset is
hosted at _https://small-obstacle-dataset.github.io/_

##  I INTRODUCTION

Within the realm of perception for autonomous cars, small obstacle detection
is a critical problem. Small obstacles fall precariously on the border of
being classified as drivable space or obstacles. It is prudent for the
planning module of an autonomous car to be informed of the small obstacles in
its environment.

A straightforward extension of semantic segmentation methods that label
semantically contiguous regions accurately, while occupying a tiny part of the
overall image space has proven difficult [1]. Since such small area classes
occupy an insignificant portion in the predicted label space, they accrue
meager error costs that result in diminished gradient magnitudes, providing
negligible updates to the model parameters. To overcome these, methods such as
[2, 3] propose cost terms sensitive to small area class labels (median
frequency balancing), which enjoy limited success. However, such methods have
not showcased their results on obstacles lying on the road but small obstacle
classes like poles and traffic signs seen from afar. Obstacles lying on the
road pose a greater challenge because such instances are rare in carefully
curated autonomous driving datasets. Additionally, they may not be easy to
spot if the appearance is similar to the road (an example is shown in Figure
1, where large stones lie on a cemented road). In this paper, we develop an
alternate paradigm - one that fuses LiDAR and Image inputs effectively for on-
road small obstacle detection.

In our method, LiDAR scans, available from the sparse Velodyne Puck 16, are
used to detect possible obstacle regions in a scan through a discontinuity
detector. A contiguous set of 3D points between two discontinuities form a
candidate region where a small obstacle may reside. These candidate regions
are projected onto the image plane through LiDAR-Camera extrinsics and are
further blurred through the Gaussian blur operator. Such blurred regions in
the image are assigned probability scores that represent the confidence of a
small obstacle and henceforth called as confidence maps. The confidence map
along with the input RGB image forms the four-channel input to a deep
convolutional network architecture. The network - supervised on three classes,
namely the road class, off-road class, and small obstacle class, - can
accurately detect small obstacle regions hitherto difficult for a purely
image-based classification framework (a motivating example is shown in Figure
1). We show significant improvement in detection and classification accuracy
of small obstacle regions by such fusion of LiDAR and image-based sensing
modalities vis a vis image only classifier architectures. Formally our paper
makes the following contributions:

|
---|---
(a) | (b)
|
(c) | (d)
Figure 1: Example detections: (a) Original image with two small obstacles
(green rectangles). (b) Monocular RGB baseline fails to detect both the
obstacles. (c) The context (confidence map) generated using LiDAR point cloud.
(d) The proposed method (combining monocular image with LiDAR context)
successfully detects both the obstacles.  Figure 2: The proposed small
obstacle detection pipeline. The confidence module generates a confidence map
based on the point cloud. The result is concatenated with the current image
and used to predict a segmentation mask.

  1. 1.

A novel LiDAR-Camera small obstacle segmentation dataset that comprises over
3000 frames with per pixel annotation of 3 classes, on-road, off-road, and
small obstacle along with LiDAR-Camera calibration extrinsics. This dataset
can also be viewed as an image only or LiDAR only dataset for small obstacle
segmentation.

  2. 2.

A novel pipeline that combines confidence map information obtained from LiDAR
point clouds along with RGB channel inputs to obtain significant boost over
purely image-based architectures such as [4, 5]. It is important to note that
sparse LiDAR scans can miss small obstacles entirely, and hence it is not
prudent to rely on the current LiDAR input alone. This is overcome by
propagating prior LiDAR confidence maps to the current frame, such that
continuous per frame small obstacle detection is deemed possible even when
LiDAR scans miss the small obstacle completely. We portray these results in
section V-B. The use of the sparse Puck 16 LiDAR over its denser 64 scan-line
counterpart, which is more amenable to point cloud segmentation, is yet
another cornerstone of this work.

  3. 3.

Detailed ablation studies across alternative methods for various modules that
constitute the overall pipeline, such as LiDAR-Camera calibration, confidence
map generation, and temporal propagation of confidence maps.

  4. 4.

A novel Haussdorf distance-based markerless LiDAR-Camera calibration
refinement that significantly improves training time and accuracy as a result
of improved extrinsics.

  5. 5.

And finally, to the best of our knowledge, this is the first such method that
combines sparse LiDAR data with a monocular image to segment small obstacles,
unlike most segmentation approaches that operate exclusively on only image or
dense point cloud data.

##  II RELATED WORK

The field of semantic segmentation has seen tremendous success over the past
few years. The success can be attributed to the applicability of CNN based
encoder-decoder architectures [3, 4, 5] and the availability of large datasets
which allow efficient training [6, 7]. However, handling class imbalance has
always remained a challenge as small objects usually contribute less to the
segmentation loss. Many semantic segmentation works follow a relatively simple
cost-sensitive approach via an inverse frequency rebalancing scheme, e.g., [8]
or median frequency re-weighting [2]. The small classes, such as pole and
signboards, in autonomous driving datasets like Cityscapes are unable to
supervise a network to detect small-obstacles lyingon the road. Apart from
class imbalance (and class mismatch), monocular small obstacle detection also
faces the challenge of weaker appearance cues (Figure 1). To highlight these
concerns, we use two state of the art image segmentation networks DeepLab
[4]and HRNet [5] as our baselines. Both networks were pretrained using the
Cityscape dataset and were fine-tuned on the proposed IIIT Small obstacle
detection dataset.

Specific efforts have been made in the area of on-road small obstacle
detection. The carefully curated Lost and Found dataset [9] has been pivotal
for the progress in this direction. Most of the approaches augment the
appearance cues with a depth map obtained using a stereo camera. The work by
Pinggera _et al._ [9] performs statistical hypothesis tests in disparity space
directly on stereo image data, assessing free space and obstacle hypotheses on
independent local patches. Ramos _et al._ [10] combines deep learning
approaches with geometric cues. MergeNet [11] proposes a novel architecture
combining local structural cues obtained from a Stripe based network with a
full image output. MergeNet takes inspiration from Stixel based processing
[12] and demonstrates the ability to be trained from only a few images (as low
as 135). The major limitation of these approaches is the reliance on
inaccurate depth maps, especially in low texture scenarios or in cases where
the obstacle appearance is similar to that of the road.

In this paper, we employ a LiDAR + RGB pipeline for on-road small obstacle
detection. LiDAR data has been successfully used for problems like road object
segmentation [13], curb detection [14] etc. Combining RGB data has shown to
bring improvements in both the tasks of semantic segmentation [15] and object
detection [16] (in contrast to pure LiDAR methods). The fusion of the two
modalities particularly shows improvements for detecting objects at a long-
range (>40m-70m) over their RGB and LiDAR only counterparts [15]. Our work is
in a similar spirit. Furthermore, the small obstacle detection problem is
difficult to be solved using LiDAR alone because (a) the small obstacle may
not intersect with LiDAR rings, especially when using sparse LiDAR (as used in
our work) and (b) even if the discontinuities are detected the overall extent
of the object will be unknown in the image space. We propose a method fusing
sparse 16 channel LiDAR with monocular data and provide pixel-level
segmentation in the image space.

Our work is also related to the literature in LiDAR and RGB calibration [17,
18]. We show that minor calibration errors can be extremely problematic for
the task of small obstacle detection. Consequently, we propose a method for
fine refinement of extrinsic parameters for LiDAR-RGB calibration. The
intuition of the approach is to align the detected LiDAR discontinuities on
the road with nearby ground truth small obstacle segmentation labels in image
space. The refinement algorithm leverages Hausdorff distance for the task and
is agnostic to the original calibration algorithm.

##  III METHODOLOGY

###  III-A Overview

Figure 3: Break point identification (for point $R$). The points $P$, $Q$ and
$R$ are consecutive points that lie on the same LiDAR ring. $O$ represents the
centre of the LiDAR sensor. $OP$, $OQ$, $OR$ and $OR^{\prime}$ corresponds to
$d_{i}$, $d_{i+1}$, $d_{i+2}$ and $d_{p}$ respectively.

The proposed pipeline is illustrated in Figure 2. At a given time, the input
is an RGB image and the corresponding LiDAR point cloud. The output is a
three-class (road, small obstacle, and off-road) semantic segmentation map. At
first, the LiDAR point cloud is used to detect discontinuities (breakpoints)
on the road using geometric reasoning (the detected discontinuities are
illustrated with green rectangles in Figure 2). The detected discontinuities
are then projected onto the image space and are used to generate a pixel level
confidence map. Since the LiDAR used in our experiments is sparse, it may miss
out on some of the small obstacles (in cases where no LiDAR point falls on
them). The confidence map is further augmented using temporal propagation of
discontinuities detected in previous frames. The confidence map with the RGB
image (4D input) and the corresponding ground truth maps are then used for
training a deep convolutional network for semantic segmentation. We now
describe each of these steps in detail.

###  III-B Obstacle Confidence Maps from LiDAR

Breakpoint Detection: We identify geometric break points within each ring of
the Point Cloud, where a ring is defined as one complete 360∘ scan for a given
channel of the LiDAR. Our sensor VLP-16 gives 16 such rings within a vertical
resolution of 30∘ (-15 to 15). We use the approach discussed in [19] to
isolate points of depth discontinuity or break points within each such ring.
As illustrated in Figure 3, for a given triplet of consecutive points,
$P=\\{p_{i},p_{i+1},p_{i+2}\\}$ having distances
$D=\\{d_{i},d_{i+1},d_{i+2}\\}$ from the LiDAR origin, we utilise the LiDAR’s
horizontal angular resolution angle $\theta$, and the measured distances
$d_{i},d_{i+1}$ at $p_{i}$ and $p_{i+1}$ respectively to predict the distance
$\mathit{d_{p}}$ at $p_{i+2}$.

$\mathit{d_{p}}=\frac{d_{i}d_{i+1}}{2d_{i}\cos{\theta}-d_{i+1}}$ (1)

The point $p_{i+2}$ is categorised as a break point if the difference between
the predicted distance $\mathit{d_{p}}$ and measured distance $d_{i+2}$ is
beyond a certain threshold $\mathbf{d_{th}}$.

$|d_{i+2}-\mathit{d_{p}}|\geq d_{th}$ (2)

We then filter out the breakpoints detected outside the road. For sequences
where road curbs and pedestrian pathways are elevated from the road plane, we
use our detected breakpoints to isolate the road segments, ring wise within
the Point Cloud, in a similar fashion to what has been discussed in [19]. We
use a pre-trained SqueezeSeg [13] for road segmentation when the curbs are not
available. Our method is agnostic to the choice of road segmentation
algorithm, and any pre-trained algorithm either in image or LiDAR space can be
used for the task.

Figure 4: LiDAR points projected onto the image space. Points within each pair
of breakpoints (denoted in Red) represent a small obstacle segment (AB).
Figure 5: Temporal Aggregation. (A) Current Confidence Map generated by break
point detection and segment isolation. Note that the two discontinuities in
the Point Cloud (depicted within green rectangles) correspond to current
Confidence Map. A third obstacle is missed by the current LiDAR scan. (B) A
template of the third small-obstacle is retrieved from contexts given in the
past frames. (C) Using the template extracted and the current image, temporal
aggregation is performed to generate an aggregated temporal confidence mask

Small Obstacle Segment Isolation: We define gradient direction at a breakpoint
A as $G(A)=sign(d_{i+2}-\mathit{d_{p}})$. A consecutive pair of breakpoints A
and B are said to belong to an obstacle segment if $G(A)$ is negative and
$G(B)$ is positive. Small obstacle segments are then filtered out based on the
horizontal/azimuthal angular resolution between the breakpoints A and B
(example detected small obstacle segments as a pair of breakpoints are
illustrated in Figure 4). This spread/resolution differentiates big obstacles
such as cars from their smaller counterparts. For this dataset, we define this
threshold as $\mathbf{2^{\circ}}$ and set the value of $\mathbf{d_{th}}$ as
0.4. This allows us to detect an obstacle as small as 15 cms in height from
about 50 meters away.

Confidence Maps: To generate a pixel level confidence map, we project the
detected small obstacle segments onto the image plane. The projection is
performed using the mapping $\phi:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}$:

$\mathbf{x}=\phi(\mathit{X},\mathit{K},\mathit{T}).$ (3)

where $\mathit{X}$ is the set of 3-D LIDAR points belonging to a small
obstacle segment and $\phi$ is the projection operator. $\phi$ is parametrized
by camera intrinsic matrix ($\mathit{K}$) and extrinsics matrix
($\mathit{T}$). The set of projected points, $\mathit{x}$, then serve as
anchors of our confidence regions, where each point in this set is assigned a
confidence value of 1 and the confidence values in the neighbourhood of this
point follow a _Gaussian_ distribution (with mean centered on the point and
variance parameterized by $\sigma$).

###  III-C Temporal Propagation of Confidences

We further augment the confidence maps using the temporal propagation of the
previously detected small obstacles. This step helps to gather context for
obstacles which are missed by the sparse LiDAR in the current frame. The
updated confidence map is then used as the input to the segmentation network.
The propagation is performed in a two step process: (a) we select a Region of
Interest (ROI) in image space around the projections of the previously LiDAR
detected small obstacle segments and (b) we find the corresponding ROI in the
current frame using template matching in RGB space. We reduce the search space
during template matching using the coarse odometry estimates between the past
and the current frame. An example is illustrated in Figure 5.

##  IV SMALL OBSTACLE DATASET

We introduce the Small Obstacle Dataset to the community considering the
importance of this problem for on-road navigation. Occurrences of small
obstacles like bricks, granite rocks, stones may not be too uncommon in some
parts of the world. Such objects are often on roads adjacent to or bordering
large scale construction sites. Sometimes it is possible that sentient beings
such as cats and dogs may be asleep on the roads. To this end, we curate a
novel IIIT Small Obstacle Dataset encompassing significant variations in terms
of obstacles, road types and lighting.

###  IV-A Description

The dataset includes monocular RGB images and synchronized LiDAR scans with
Odometry information. It comprises of 15 sequences in total and utilises a
varied set of small obstacle instances. The details in terms of total images
and train/val/test splits are given in Table I We used different roads and
different set of obstacles while recording the train, val and test sequences.
Test split is kept to be most challenging in terms of turns, occlusions and
shadows to better evaluate the generalizability. We kept the definition of a
small obstacle as an object whose longest dimension is less than the diameter
of a standard car wheel - which is around 21 inches.

####  IV-A1 Sensor and Software Setup

The Sensor and Software setup for recording the data is as follows

  * •

ZED Stereo camera (only left feed).

  * •

Velodyne Puck (VLP-16).

  * •

Vehicle - Mahindra E2O electric car.

  * •

ROS Kinetic

Split | # seq. | # Images | with curbs | without curbs
---|---|---|---|---
Train | 9 | 1937 | 6 | 3
Validation | 4 | 530 | 2 | 2
Test | 2 | 460 | 1 | 1
Total | $\mathbf{15}$ | $\mathbf{2927}$ | 9 | 6
TABLE I: Description of the Small Obstacle Dataset

####  IV-A2 Annotations

We provide accurate pixel-level semantic annotations specific to the task of
small obstacle segmentation for a selected set of images for all the sequences
(sampling more densely around small obstacles). The total number of annotated
images is 2927(Table I). We also release all the raw sequences. The dataset is
segregated into three classes, namely, Road, Off-Road and Small-Obstacle.
Everything except road and small obstacles is labelled as off-road class
(buildings, cars, pedestrians etc.). The Small-Obstacle class is further
annotated with nine different types of sub-classes (stones, bricks, plastic,
dogs etc.).

###  IV-B Improving LiDAR-Camera Calibration

LiDAR and Camera provide complementary streams of information in depth and
colour and have become an indispensable part of the perception module for
autonomous driving. To make vital environmental inferences, it is essential
that they must be well calibrated so that their mutual information can be
fused together. The problem hence involves finding the 6-DoF rigid body
transformation matrix which allows for flexible transformation between their
respective co-ordinate frames. This is done by solving for point-
correspondences either between 2D-3D points i.e pixel co-ordinates in Camera
frame and 3-D points in Point-Cloud or through 3D-3D associations. In the
latter, specific setups like AR based ArUCo markers are used, which facilitate
retrieving the 3-D coordinates in the camera frame through encoded patterns.
We use this method for our setup utilising ArUco markers as detailed in [18].

While the method is specifically adapted to work with sparse 16 channel LiDAR
over other methods - which use a denser 64 channel LIDAR(for eg. used in KITTI
dataset [20]), it still has its limitations. It relies on recognition of the
patterns in ArUco marker to calculate the 3-D correspondences for the Camera
frame, therefore the markers are kept within a few metres of the setup during
the optimisation process. We observed that the calibration matrix obtained
using point correspondences, estimated using this method(i.e corners of the
marker), were not accurate at greater distances( 50m), posing a significant
problem for our task (an example is illustrated in Figure 6). We therefore
propose a method to do fine refinement over the extrinsics parameters
specifically for tasks where precise calibration is essential for sensor
fusion. Below we explain our methodology and show Quantitative results in V of
improvement over the original task with this refined calibration method.

Figure 6: (a) Point Cloud projection using calibration matrix obtained through
[18]. (b) Improved calibration using our method. We can observe how projection
error becomes more significant for farther obstacles.

2D-3D Point Correspondence: The main intuition of our approach is that if the
calibration is accurate, the small obstacle segments detected in 3D LiDAR
point cloud when projected onto 2D image space should intersect with ground
truth annotations of the small obstacle class. Formally, let $\mathit{S_{1}}$
denote the set of 3D points falling on multiple small obstacles and
$\mathit{P_{2}}$ denote the set of 2D points belonging to the corresponding
annotated pixel label in Image. We optimise over the extrinsic calibration
parameters such that the projection set $\mathit{P_{1}\in\mathbb{R}^{2}}$ of
$\mathit{S_{1}}$ obtained through (4), has the following property:
$\mathbf{P_{1}}\subseteq\mathbf{P_{2}}$, i.e. all detected small obstacle
segments in the point cloud lie within their respective pixel label bounds.
The transformation operation $\phi:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}$ is
defined as:

$\begin{pmatrix}x\\\ y\end{pmatrix}=\pi\begin{pmatrix}R\begin{pmatrix}X\\\
Y\\\ Z\end{pmatrix}+t\end{pmatrix}$ (4)

where $\pi$ subsumes the camera intrinsics ($f_{x}$,$f_{y}$,$c_{x}$,$c_{y}$)
and R, t are the rotation and translation matrix in SO(3) respectively.

Let the 1 x 6 vector $\xi$ = $(\nu,\omega)\in$ se(3) denote our initial coarse
calibration extrinsics where [$\nu_{x},\nu_{y},\nu_{z}$] is the translation
vector and [$\omega_{x},\omega_{y},\omega_{z}$] is the rotational vector along
their respective axis. The rotational vector $\omega$ $\in$ so(3) can be
converted to R $\in$ SO(3) using the exponential map exp: so(3) $\to$ SO(3);
$\omega$ $\to$ $e^{\omega}$.

$R=I+\frac{w}{||w||}\sin{||w||}+\frac{w^{2}}{||w||^{2}}(1-\cos{||w||^{2}})$
(5)

According to our chosen point correspondences stated above, we now define the
projection loss through Hausdorff distance.

Hausdorff Projection Loss: The directed Hausdorff distance between two sets
$\mathbf{P_{1}}$ and $\mathbf{P_{2}}$ is defined as:

$d(P_{1},P_{2})=\underset{x\in P_{1}}{\sup}\underset{y\in P_{2}}{\inf}||x-y||$
(6)

The loss function thus gives a measure of the distance the set $P_{1}$ has to
travel in Image space so as to be completely contained inside of $P_{2}$.
Gradients with respect to each of the 6 elements of the se(3) transformation
vector $\xi$ can then be back-propagated while minimising the loss with any
gradient based optimizer. For the experiment, we sample a few frames from each
sequence of the Training split of the Dataset and use Adam [21] optimizer with
a learning rate of $1\text{e-}5$.

##  V EXPERIMENTS AND DISCUSSION

###  V-A Evaluation metrics

We evaluate our model’s performance on both instance-level and pixel-wise
metrics:

Instance Detection Rate (IDR): Instance-level detection rate (IDR) is defined
as the fraction of obstacle instances, taken across the dataset, which are
detected by the network. For this metric, an instance is marked correctly
detected if more than 20% of the pixels of the predicted obstacle overlaps
with the ground truth of that instance. For extracting instances from pixel-
level predictions we use a similar approach to [11]. The IDR metric is
formally calculated as:

$\mathit{IDR}=\frac{TP_{obstacle}}{TI_{obstacle}}$ (7)

Where a true positive, $TP_{obstacle}$, represents a set of predicted
instances that have an overlap greater than a threshold of 0.2 with the
ground-truth. $TI_{obstacle}$ is the total instances of the obstacle class,
taken across the entire dataset.

instance False Detection Rate (iFDR) This metric is defined as the fraction of
instances that have been incorrectly detected as small obstacles:

$\mathit{iFDR}=\frac{FP_{obstacle}}{PRED_{obstacle}}$ (8)

Where a $FP_{obstacle}$ is defined as any predicted small obstacle prediction
instance that has no overlap with ground truth and $PRED_{obstacle}$ are the
total number of instances predicted as small obstacles.

Pixel Detection Rate (PDR) $\mathit{PDR}$ is the fraction of small obstacle
pixels that have been correctly predicted.

$\mathit{PDR_{obstacle}}=\frac{TPX_{obstacle}}{GTX_{obstacle}}$ (9)

where $TPX_{obstacle}$ denotes the total number of correctly predicted small
obstacle pixels and $GTX_{obstacle}$ is the total number of small obstacle
pixels in the image.

mean Intersection over Union (mIoU) The mean Intersection over Union (mIoU) is
a commonly used metric for semantic segmentation. It measures the average
Intersection over Union across all classes.

$\mathit{IoU=\frac{TPX}{TPX+FPX+FNX}}$ (10)

Where, $TPX$, $FPX$ and $FNX$ are pixel level True Positives, False Positives
and False Negatives for the small obstacle class.

Network | Method | Instance-level | Pixel-level
---|---|---|---
IDR | iFDR | PDR | mIoU
DeepLab-V3+[4] | Image | 0.39 | 0.28 | 0.37 | 0.73
Image + CM | 0.50 | 0.26 | 0.45 | 0.74
Image + CM + TP | 0.60 | 0.18 | 0.60 | 0.76
HRNet[5] | Image | 0.44 | 0.25 | 0.27 | 0.70
Image + CM | 0.47 | 0.25 | 0.32 | 0.70
Image + CM + TP | 0.63 | 0.21 | 0.51 | 0.73
TABLE II: Performance comparison between various inputs to the semantic
segmentation architectures.

###  V-B Experiments

We compare our method with two monocular semantic segmentation baselines. We
present thorough ablations to motivate each component of the proposed
pipeline. The evaluations are performed on the test set of the proposed IIIT
Small Obstacle Dataset. Each experiment is trained on 2 Nvidia GTX-1080ti GPUs
for 15 epochs. A batch size of 6 is kept during training.

As to the choice of the segmentation network, we picked two state of the art
networks DeepLabV3+ [4] and HRNet [5] on CityScapes benchmark. They also made
a suitable choice due to availability of pre-trained model/weights and well
documented reproducible code base. We perform three set of experiments
considering pre-trained DeepLabV3+ and HRNet as base networks. The pre-
training is performed on the CityScapes dataset.

Image baseline:  Our baseline for comparison is the performance of base
network on the Small Obstacle Dataset, when trained on images only. We fine-
tune a pre-trained base network on the proposed dataset to give a 3 class per
pixel prediction. We supervise the network using the Cross-Entropy loss with
inverse frequency rebalancing.

Image + LiDAR confidence masks: Next, we modify the base-network by
initialising a new channel for the input layer so as to accept a 4D input
tensor. We then train the network on images concatenated with confidence masks
(generated using LiDAR projections) along the channel dimension. This
experiment only considers the current LiDAR for generating the confidence
maps.

Image + LiDAR confidence masks + TP: This experiment is similar to (Image +
LiDAR) in network configuration. However, the temporal context (4 previous
frames) is considered to generate the confidence maps (to compensate of
obstacles missed in the current frame).

###  V-C Results

In Table II we evaluate the performance of our framework. All metrics are
reported for the test set. The DeepLab-V3+ image-only baseline achieves an
$IDR$ of 0.39. Our proposed method of temporally aggregated confidence maps
shows a 53% improvement to an $IDR$ of 0.60, while reducing the $iFDR$ by 35%
to a score of 0.18. Furthermore, this improvement comes with a 4% increase to
0.76 in the mIoU metric.

A similar trend is seen when our framework is trained with the HRNet
segmentation network. The HRNet image-only baseline achieves a higher $IDR$
than DeepLab-V3+ at 0.44. When trained with our proposed temporally aggregated
confidence masks along with images, performance improves by 43% to a score of
0.63, while the $iFDR$ drops by 16% to a score of 0.21.

We further note the increase in the $IDR$ metric between utilising obstacle
confidences only from the current frame (Image + CM), and temporally
propagating such confidences from previous frames (Image + CM + TP). As can be
seen in Table II, there is a 20% corresponding increase from 0.50 to 0.60 by
augmenting the confidence map using temporal observations. This reinforces the
importance of temporal confidence propagation.

We present a few qualitative results across our test and validation splits in
Fig. 7. It can be observed that even when there is a high confidence value on
legitimate regions within the image (in Fig. 7, legs of a man, wheel of a
car), the joint representation has learnt to classify these regions as non-
obstacle(road/off-road).

Method | IDR | PDR
---|---|---
RANSAC | 0.424 | 0.43
Conv-1D | 0.421 | 0.373
Geometric (ours) | 0.50 | 0.45
TABLE III: A comparison of various small obstacle detection techniques on
Point Cloud.

##  VI ABLATION STUDY

We show the following ablation experiments across our three main
contributions: (a) Detection of small obstacles in the LiDAR Point Cloud, (b)
Temporal propagation of confidences and (c) improved extrinsics calibration.
All of the results presented below use DeepLab-V3 as the base network
architecture.

####  VI-1 Detection of Obstacle Segments

We experiment with multiple approaches for detecting a small obstacle within a
Point Cloud ring. In Table-III we show our results with three such methods for
the Image + Confidence Map Baseline. We first experiment with a plane-fitting
method where points falling on small obstacles are classified as outliers
through Random sample consensus [22]. Through this we can successfully detect
bigger small obstacle instances however it fails on small instances such as
rocks,bricks etc. We then study a learning based method where a small 1-D
convolutional network is trained to classify each point within the curb
boundaries as a road or an obstacle point. The best results were obtained
using geometric approach (Sec III-B) and was thus used in our final pipeline.

####  VI-2 Temporal Propagation of Confidences

We tried two different methods for temporal propagation and observed their
impact on the final output (Table IV). In Forward Projection, we directly
project obstacle detections in previous 4 frames onto the current Image frame
using Odometry information calculated through LiDAR Odometry and Mapping[23].
However due to less than desired accuracy in the Odometry estimate this
doesn’t localize the confidence map to the true small obstacle instance on the
Image and shows minor improvement. The proposed temporal aggregation (Sec
III-C) gives much improved results.

####  VI-3 Improved Extrinsics Calibration

Here we study the effect of better calibration on our method comparing the
results on Image + Confidence Map baseline. In Table-V it can be seen how
better calibration significantly affects the detection of small obstacle
instances. It is to be noted that increasing the spatial spread of the
Confidence region/map with coarse extrinsics (Pre-Hausdorff) does not
translate to a comparable performance with a better set of extrinsics obtained
using our method. This highlights the need for calibration refinement in our
pipeline.

Method | IDR | iFDR | PDR
---|---|---|---
Forward Projection | 0.52 | 0.23 | 0.43
Temporal Aggregation | 0.60 | 0.18 | 0.60
TABLE IV: A comparison of methods for temporal propagation of detected obstacles. Extrinsics Parameters | $\sigma$ | IDR | iFDR | PDR
---|---|---|---|---
Pre-Hausdorff | 5 | 0.412 | 0.263 | 0.363
7 | 0.426 | 0.217 | 0.36
Post-Hausdorff | 5 | 0.50 | 0.26 | 0.45
7 | 0.49 | 0.22 | 0.464
TABLE V: Comparison on the detection performance between Coarse and Fine
calibration obtained using our method. $\sigma$ denotes the variance of the
Gaussian blur.

##  VII CONCLUSIONS

Our work focuses on the problem of detecting small obstacles lying on road by
fusing monocular RGB and LiDAR data. Our experiments show that IDR increases
by 40-50$\%$ when a given image only segmentation network is augmented with
Confidence Maps provided by LiDAR. Thorough experiments using two state of the
art base networks DeepLABV3+ and HRNet demonstrate the efficacy of our
approach. We also show that small obstacle ground truth annotations can be
exploited to improve the extrinsic calibration, which in turn improves the IDR
(Table- V). We further present thorough ablation studies to justify the design
choices. Overall, our method is able to detect 73% of the small obstacles (as
small as 15cm high) within the range of 50 meters.

|  |  |
---|---|---|---
|  |  |
|  |  |
|  |  |
(a) | (b) | (c) | (d)
Figure 7: (a) Image. (b) Confidence Map (c) Prediction. (d) Ground-truth

## References

  * [1] D. Guo, L. Zhu, Y. Lu, H. Yu, and S. Wang, “Small object sensitive segmentation of urban street scene with spatial adjacency between object classes,” _IEEE Transactions on Image Processing_ , vol. 28, no. 6, pp. 2643–2653, 2018.
  * [2] D. Eigen and R. Fergus, “Predicting depth, surface normals and semantic labels with a common multi-scale convolutional architecture,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 2650–2658.
  * [3] V. Badrinarayanan, A. Kendall, and R. Cipolla, “Segnet: A deep convolutional encoder-decoder architecture for image segmentation,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 39, no. 12, pp. 2481–2495, 2017.
  * [4] L. Chen, G. Papandreou, F. Schroff, and H. Adam, “Rethinking atrous convolution for semantic image segmentation,” _CoRR_ , vol. abs/1706.05587, 2017. [Online]. Available: http://arxiv.org/abs/1706.05587
  * [5] J. Wang, K. Sun, T. Cheng, B. Jiang, C. Deng, Y. Zhao, D. Liu, Y. Mu, M. Tan, X. Wang, W. Liu, and B. Xiao, “Deep high-resolution representation learning for visual recognition,” _CoRR_ , vol. abs/1908.07919, 2019.
  * [6] M. Cordts, M. Omran, S. Ramos, T. Rehfeld, M. Enzweiler, R. Benenson, U. Franke, S. Roth, and B. Schiele, “The cityscapes dataset for semantic urban scene understanding,” in _Proc. of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016.
  * [7] M. Everingham, S. A. Eslami, L. Van Gool, C. K. Williams, J. Winn, and A. Zisserman, “The pascal visual object classes challenge: A retrospective,” _International journal of computer vision_ , vol. 111, no. 1, pp. 98–136, 2015.
  * [8] J. Xu, A. G. Schwing, and R. Urtasun, “Tell me what you see and i will show you where it is,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2014, pp. 3190–3197.
  * [9] P. Pinggera, S. Ramos, S. Gehrig, U. Franke, C. Rother, and R. Mester, “Lost and found: detecting small road hazards for self-driving vehicles,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 1099–1106.
  * [10] S. Ramos, S. Gehrig, P. Pinggera, U. Franke, and C. Rother, “Detecting unexpected obstacles for self-driving cars: Fusing deep learning and geometric modeling,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2017, pp. 1025–1032.
  * [11] K. Gupta, S. A. Javed, V. Gandhi, and K. M. Krishna, “Mergenet: A deep net architecture for small obstacle discovery,” _CoRR_ , vol. abs/1803.06508, 2018. [Online]. Available: http://arxiv.org/abs/1803.06508
  * [12] D. Pfeiffer and U. Franke, “Towards a global optimal multi-layer stixel representation of dense 3d data.” in _BMVC_ , vol. 11, 2011, pp. 51–1.
  * [13] B. Wu, A. Wan, X. Yue, and K. Keutzer, “Squeezeseg: Convolutional neural nets with recurrent crf for real-time road-object segmentation from 3d lidar point cloud,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 1887–1893.
  * [14] G. Zhao and J. Yuan, “Curb detection and tracking using 3d-lidar scanner,” in _2012 19th IEEE International Conference on Image Processing_. IEEE, 2012, pp. 437–440.
  * [15] G. P. Meyer, J. Charland, D. Hegde, A. Laddha, and C. Vallespi-Gonzalez, “Sensor fusion for joint 3d object detection and semantic segmentation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2019, pp. 0–0.
  * [16] C. R. Qi, W. Liu, C. Wu, H. Su, and L. J. Guibas, “Frustum pointnets for 3d object detection from rgb-d data,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 918–927.
  * [17] F. M. Mirzaei, D. G. Kottas, and S. I. Roumeliotis, “3d lidar–camera intrinsic and extrinsic calibration: Identifiability and analytical least-squares-based initialization,” _The International Journal of Robotics Research_ , vol. 31, no. 4, pp. 452–467, 2012.
  * [18] A. Dhall, K. Chelani, V. Radhakrishnan, and K. M. Krishna, “Lidar-camera calibration using 3d-3d point correspondences,” _arXiv preprint arXiv:1705.09785_ , 2017.
  * [19] W. S. Wijesoma, K. S. Kodagoda, and A. P. Balasuriya, “Road-boundary detection and tracking using ladar sensing,” _IEEE Transactions on robotics and automation_ , vol. 20, no. 3, pp. 456–464, 2004.
  * [20] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _International Journal of Robotics Research (IJRR)_ , 2013\.
  * [21] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [22] M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,” _Communications of the ACM_ , vol. 24, no. 6, pp. 381–395, 1981.
  * [23] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time.” in _Robotics: Science and Systems_ , vol. 2, no. 9, 2014.


# FOUR-DIMENSIONAL POLARISATION-RING-SWITCHING FOR DISPERSION-MANAGED OPTICAL
FIBRE SYSTEMS

Bin Chen1,2 Chigo Okonkwo1 Hartmut Hafermann3 Alex Alvarado1 1Department of
Electrical Engineering, Eindhoven University of Technology, The Netherlands.
<EMAIL_ADDRESS>
2School of Computer and Information, Hefei University of Technology, China.
3Mathematical and Algorithmic Sciences Lab, Paris Research Center, Huawei
Technologies France SASU, France.

###### Abstract

The recently introduced 4D 64-ary polarisation-ring-switching format is
investigated in dispersion-managed systems. Numerical simulations show a reach
increase of $25\%$ with respect to PM-8QAM. This gain is achieved from the
nonlinear tolerance of the format and a 4D demapper using correlated noise
assumptions.

###### keywords:

CODED MODULATION, GENERALISED MUTUAL INFORMATION, DISPERSION MANAGED SYSTEM.

##  1 Introduction

Fibre nonlinearities are considered to be one of the limiting factors for
achieving higher information rates in coherent optical fibre transmission
systems EssiambreJLT2010 ; SecondiniJLT2013 . Multiple digital signal
processing (DSP) algorithms have been investigated to reduce the impact of
intra-channel and inter-channel nonlinearities Cartledge:17 . One appealing
alternative is to properly design multi-dimensional (MD) modulation formats,
which can minimize the signal power variations, and thus, reduce the nonlinear
interference (NLIN) power. This property is especially suitable for
dispersion-managed (DM) links with a strong nonlinear cross-phase modulation
(XPM) Chagnon:13 ; ReimerOFC2016 ; Kojima2017JLT .

We recently proposed in BinChenJLT2019 the 4D 64-ary polarisation-ring-
switching (4D-64PRS) format, which was obtained by jointly optimising the
coordinates and labeling of the constellation. 4D-64PRS is a nonlinearity-
tolerant format which was shown to outperform other formats at the same
spectral efficiency (6 bit/4D-sym) in both the linear and nonlinear channels.
All the reported results in BinChenJLT2019 were obtained for dispersion
unmanaged systems only. The first contribution of this paper is to extend
those results to DM links.

In this paper, we use generalised mutual information (GMI) as a figure of
merit. GMI is very popular because it directly shows the net data rates
obtained by state-of-the-art soft-decision forward error correction (SD-FEC)
AlvaradoJLT2018 . GMI calculations can include the demapper structure at the
receiver. In this case, the GMI is calculated using log-likelihood ratios
(LLRs, or soft bits), which are in turn calculated using a particular
assumption on the channel.

Typically, circularly symmetric Gaussian statistics are assumed in the LLR
calculation. This assumption is accurate for uncompensated long-haul fibre
systems ErikssonJLT2016 ; KeykhosraviARXIV2018 . However, in systems with
strong phase noise (e.g., due to residual phase noise) or large nonlinear
fibre effects (like in a DM link), this assumption is no longer valid.
Transmission experiments have shown a decrease in achievable information rates
and coding gain if there is a strong mismatch between the true channel and the
channel assumption made by the demapper ZhaoOFC2013 ; ErikssonJLT2016 . The
second contribution of this paper is to show that relatively large GMI gains
are obtained when 4D-64PRS is combined with a demapper that lifts the
circularly symmetric Gaussian assumption.

In this paper, we show the nonlinearity tolerance of 4D-64PRS over PM-8QAM in
DM links and a demapper that treats the four dimensions jointly and allow
correlations between them. It is shown that 4D-64PRS outperforms PM-8QAM by up
to $0.55$ bit/4D-sym in terms of GMI ($12\%$ rate increase) for a 2400 km DM
transmission link. 4D-64PRS also yields 600 km ($25\%$) reach increase.
Finally, we also show the gain increases as the number of WDM channels
increase. This leads us to believe that 4D-64PRS is well-suited for ultra
wideband fibre systems.

##  2 The 4D-64PRS Modulation Format

Figure 1: 2D-projections of the designed 4D-64PRS modulation and associated
bits/polarisation bits mapping. The rings are given by
$R_{1}^{2}={\nu_{1}^{2}+\nu_{3}^{2}}$ and $R_{2}^{2}=2\nu_{2}^{2}$.

Fig. 1 shows the coordinates and the corresponding binary labeling of the
4D-64PRS with $M=2^{m}=64$ points (6 bit/4D-sym) for two 2D-projections. The
coordinates of 4D-64PRS are such

$\displaystyle\vspace{-1.5em}\bm{s}_{i}\in$
$\displaystyle\\{[\pm\nu_{1},\pm\nu_{3},\pm\nu_{2},\pm\nu_{2}],[\pm\nu_{3},\pm\nu_{1},\pm\nu_{2},\pm\nu_{2}],$
$\displaystyle[\pm\nu_{2},\pm\nu_{2},\pm\nu_{1},\pm\nu_{3}],[\pm\nu_{2},\pm\nu_{2},\pm\nu_{3},\pm\nu_{1}]\\},\vspace{-1em}$

where $i\in\\{1,2,\cdots,64\\}$,
$\nu_{2}=\operatorname{Re}\left(R_{1}e^{j(\phi_{1})}\right)=\sqrt{1/2}R_{2}$,
$\nu_{1}=\operatorname{Re}\left(R_{1}e^{j(\phi_{2})}\right)$ and
$\nu_{3}=\operatorname{Im}\left(R_{1}e^{j(\phi_{2})}\right)$, $R_{j}$ and
$\phi_{j}$ indicate the radius and phase of inner ($j=1$) and outer ring
($j=2$). $\theta=\phi_{1}-\phi_{2}>0$ is the relative phase between inner and
outer ring points. We can observe that the 4D-64PRS is highly symmetric. The
quadrants in each polarisation (or orthants in 4D space) are defined by 4 out
of 6 bits, namely, $[b_{1},b_{2}]$ and $[b_{4},b_{5}]$ for the first and
second polarisation, respectively. These $4$ bits can be seen as the ones
determining which of the 16 orthants is to be chosen. The remaining 2 bits
[$b_{3},b_{6}$] determine the symbol from the 4 possible constellation points
in the same orthant. In other words, [$b_{3},b_{6}$] determine the color of
the transmitted points while [$b_{1},b_{2},b_{4},b_{5}$] determine the
coordinate of the point in the same color through
$[(-1)^{b_{2}}s_{i,1},(-1)^{b_{1}}s_{i,2},(-1)^{b_{4}}s_{i,3},(-1)^{b_{3}}s_{i,4}]$
with $i\in\\{1,2,\ldots,64\\}$.

##  3 Nonlinear Channel Performance

The simulation setup is as follows: a dual-polarisation multi-span WDM system
with 11 co-propagating channels are generated at the transmitter at a symbol
rate of 45 GBaud with a root-raised cosine (RRC) pulse shape with a roll-off
factor of 0.1. Each WDM channel carries $2^{16}$ symbols in two polarisations
at the same power. The WDM optical signals are launched into a multi-span
optical link. Each span consists of a 80 km large effective area fibre (LEAF)
with 0.219 dB/km attenuation, 4.255 ps/nm/km dispersion and 1.464
(w$\cdot$km)-1 nonlinear coefficient. Each span is followed by an ideal
$100\%$ in-line chromatic dispersion compensation fibre and an erbium-doped
fibre amplifier (EDFA) with a noise figure of 5 dB. The transmission link is
simulated employing a split-step Fourier solution of the Manakov equation with
a step size of 0.1 km to guarantee the accuracy. polarisation mode dispersion
(PMD) is not considered.

After match filter and ideal phase compensation, performance of the
constellation with cardinality $M=2^{m}$ is evaluated via the GMI as
AlvaradoJLT2018

$\displaystyle\vspace{-1em}\text{GMI}\approx
m-\frac{1}{N_{s}}\sum_{k=1}^{m}\sum_{j=1}^{N_{s}}\log_{2}\left(1+e^{(-1)^{b_{k,j}L_{k,j}}}\right),\vspace{-1em}$
(1)

where $L_{k,j}$ is the log-likelihood ratios (LLRs) with $k$ denoting the bit
position and $j$ denoting the $j$th received symbol with $j=1,2,\cdots,N_{s}$.
For the general case of an $N$-dimensional demapper, the LLRs are calculated
by the demapper using the channel law $f_{\bm{Y}|\bm{X}}$.

$L_{k,j}=\log\frac{\sum_{\bm{s}_{i}\in\mathcal{X}^{0}_{k}}f_{\bm{Y}|\bm{X}}(\bm{y}|\bm{s}_{i})}{\sum_{\bm{s}_{i}\in\mathcal{X}^{1}_{k}}f_{\bm{Y}|\bm{X}}(\bm{y}|\bm{s}_{i})},$
(2)

where $y$ is received symbols, $\bm{s}_{i}$ denote the $i$th constellation
point, $\mathcal{X}^{b}_{k}$ is the set of constellation symbols labeled by a
bit $b\in\\{0,1\\}$ at bit position $k$.

In this paper, we assume the channel low
$f_{\bm{Y}|\bm{X}}(\bm{y}|\bm{s}_{i})$ to be an $N$-dimensional Gaussian
distribution in ErikssonJLT2016 ; SillekensJLT2016 ,

$\vspace{-0.5em}f_{\bm{Y}|\bm{X}}(\bm{y}|\bm{s}_{i})=\frac{1}{\sqrt{(2\pi)^{N}|\sum_{i}|}}e^{-\frac{1}{2}(\bm{y}_{j}-\bm{s}_{i})^{T}\bm{\sum}_{i}^{-1}(\bm{y}_{j}-\bm{s}_{i})},$
(3)

where $|\bm{\sum}_{i}|$ is the determinant of the covariance matrix
$\bm{\sum}_{i}$ for $i$th constellation point. In this paper, we consider two
different assumptions of Gaussian distribution. In the first case, we assume
independent and identically distributed Gaussian (4D-iidG model) random
variables, and thus, $\bm{\sum}_{i}$ is a $4\times 4$ identity matrix
multiplied with the average noise variance $\sigma^{2}$, which is identical
for all $M$ constellation points. In the second case, we do not make an
assumption on the covariance matrix, and thus, we assume correlated
$N$-dimensional Gaussian noise (4D-CG model). Note that LLRs mismatched to the
channel or computed with some approximation result in a rate loss
AlvaradoJLT2018 . Different LLR correction strategies can be used to further
improve the rate as well as the performance of the FEC decoder, i.e., via LLR
post-processing AlvaradoECOC2016 .

The transmission performance is shown in Fig. 2 for three modulation formats
(PM-8QAM, 6b4D-2A8PSK Kojima2017JLT and 4D-64PRS) and two distribution models
over 2400 km. In the linear regime, 4D-64PRS offers a GMI gain of 0.3
bit/4D-sym over PM-8QAM. This gain increases to 0.4 bit/4D-sym at the optimum
launch power in the nonlinear regime. By using the 4D-CG model, the GMI
improves by additional 0.15 bit/4D-sym compared to the 4D-iidG model. The
largest obtained gain is therefore 0.55 bit/4D-sym gain ($12\%$ net data rate
increase). Fig. 2 also shows that 4D-64PRS is more tolerant to nonlinearities,
which is shown by the higher optimum launch power, and by the larger gain in
the highly nonlinear regime.

Figure 2: Measured GMI versus launch power after 2400 km.

(a) PM-8QAM

(b) 4D-64PRS

Figure 3: Recovered constellations: (a) PM-8QAM and (b) 4D-64PRS for two
polarisations after 10 spans. The x-polarisation is displayed as a 2D
projection, whilst the y-polarisation is shown as a 3D projection including
the probabilities.

Fig. 3 shows the recovered constellations for PM-8QAM and 4D-64PRS at 800 km
and optimum launch power (2D projection). The phase noise nature of NLIN is
visible in the case of PM-8QAM, but its relative significance reduces for
4D-64PRS. In addition, the noise distribution of 4D-64PRS is closer to a
circularly-symmetric distributions. From Fig. 3, it is clear that the NLIN
partly depends on the modulation format and can be mitigated by the constant-
modulus property of 4D-64PRS. Compared to uniform PM-8QAM, Fig. 3 also shows
that 4D-64PRS induces a nonuniform distribution when projected onto 2D.

The GMI of PM-8QAM, 6b4D-2A8PSK and 4D-64PRS are compared in Fig. 4. GMI of
targeted three data rates shown in the dotted lines, which correspond to 250
Gb/s, 225 Gb/s and 200 Gb/s. For the 4D-iidG model under consideration, the
4D-64PRS is capable of outperforming PM-8QAM and 6b4D-2A8PSK at all three
rates. After applying 4D-CG model on the 4D-64PRS to achieve the same GMI of
4.55 bit/4D-sym, 4D-64PRS with 4D-CG model yields a reach increase of 600 km
($25\%$) and 500 km ($20\%$) with respect to PM-8QAM with 4D-iidG model and
4D-CG model, respectively. The transmission distances achieved by 4D-64PRS at
the three rates are 1900 km, 2500 km and 3000 km, respectively.

Figure 4: The measured GMI with respect to the transmission distances for 11
WDM channels at 250 Gbps, 225 Gbps and 200 Gbps achievable data rates.

Figure 5: Dependence of the information rate at optimum launch power on the
number of channels by using 4D-CG model. The solid/dashed lines are a linear
fit of the simulation results.

To conclude, Fig. 5 shows the dependence of the information rates (in Gb/s) on
the number of channels at optimum launch power. For 2400 km, 4D-64PRS yields
an increase in information rate of 10 Gb/s ($4\%$) compared to PM-8QAM for the
5 channels case. This advantage increases with the number of channels, which
provides a gain of 20 Gb/s ($12\%$) for 41 channels. Similar gains are also
observed at a distance of 2000 km. This clearly shows the effectiveness of
4D-64PRS in suppressing nonlinearity.

##  4 Conclusions

We investigated the performance of the recently introduced 4D-64PRS modulation
format in dispersion managed links. The nonlinearity tolerance of 4D-64PRS was
shown to yield large reach and data rate increases with respect to regular
PM-8QAM. These gains were shown to becomes even higher when the number of WDM
channels increases. In this paper we also highlighted the importance of using
the correct channel model assumption in the demapper, which resulted in
additional gains. We believe 4D-64PRS is a promising candidate for upgrading
installed legacy low-dispersion fibres with high nonlinearity. The design of
higher dimensional modulation formats and low-complexity demappers are left
for future research.

Acknowledgements: This work was supported by Huawei France through the NLCAP
project. The work of B. Chen is partially supported by the National Natural
Science Foundation of China (NSFC) under Grant 61701155. The work of A.
Alvarado is supported by the Netherlands Organisation for Scientific Research
(NWO) via the VIDI Grant ICONIC (project number 15685). The authors would like
to thank Tobias Fehenberger and Eric Sillekens for useful discussion.

## References

## References

  * [1] Essiambre, R., Kramer, G., Winzer, P. J., et al.: ‘Capacity limits of optical fiber networks’, Journal of Lightwave Technology, 2010, 28, (4), pp. 662–701, 10.1109/JLT.2009.2039464
  * [2] Secondini, M., Forestieri, E., Prati, G.: ‘Achievable information rate in nonlinear WDM fiber-optic systems With arbitrary modulation formats and dispersion maps’, Journal of Lightwave Technology, 2013, 31, (23), pp. 3839–3852, 10.1109/JLT.2013.2288677
  * [3] Cartledge, J. C., Guiomar, F. P., Kschischang, F. R., et al.: ‘Digital signal processing for fiber nonlinearities’, Opt. Express, 2017, 25, (3), pp. 1916–1936, 10.1364/OE.25.001916
  * [4] Chagnon, M., Osman, M., Zhuge, Q., et al.: ‘Analysis and experimental demonstration of novel 8PolSK-QPSK modulation at 5 bits/symbol for passive mitigation of nonlinear impairments’, Opt. Express, 2013, 21, (25), pp. 30204–30220, 10.1364/OE.21.030204
  * [5] Reimer, M., Gharan, S. O., Shiner, A. D., et al.: ‘Optimized 4 and 8 dimensional modulation formats for variable capacity in optical networks’, Proc. 2016 Optical Fiber Communications Conference and Exhibition (OFC), Mar. 2016
  * [6] Kojima, K., Yoshida, T., Koike-Akino, T., et al.: ‘Nonlinearity-tolerant four-dimensional 2A8PSK family for 5-7 bits/symbol spectral efficiency’, Journal Lightwave Technology, 2017, 35, (8), pp. 1383–1391, 10.1109/JLT.2017.2662942
  * [7] Chen, B., Chigo, O., Hafermann, H., et al.: ‘Polarization-ring-switching for nonlinearity-tolerant geometrically-shaped four-dimensional formats maximizing generalized mutual information’, 2019, 37, (14), pp. 3579–3591, 10.1109/JLT.2019.2918072
  * [8] Alvarado, A., Fehenberger, T., Chen, B., et al.: ‘Achievable information rates for fiber optics: applications and computations’, Journal of Lightwave Technology, 2018, 36, (2), pp. 424–439
  * [9] Eriksson, T. A., Fehenberger, T., Andrekson, P. A., et al.: ‘Impact of 4D channel distribution on the achievable rates in coherent optical Communication experiments’, Journal of Lightwave Technology, 2016, 34, (9), pp. 2256–2266, 10.1109/JLT.2016.2528550
  * [10] Keykhosravi, K., Secondini, M., Durisi, G., et al.: ‘How to increase the cchievable information rate by per-channel dispersion compensation’, arXiv e-prints arXiv: 1812.03556, 2018
  * [11] Zhao, Y., Qi, J., Hauske, F. N., et al.: ‘Beyond 100G optical channel noise modeling for optimized soft-decision FEC performance’, Proc. 2013 Optical Fiber Communications Conference and Exhibition (OFC), Mar. 2012
  * [12] Sillekens, E., Alvarado, A., Okonkwo, C. M., et al.: ‘An experimental comparison of coded modulation strategies for 100 Gb/s Transceivers’, Journal of Lightwave Technology, 2016, 34, (24), pp. 5689–5697, 10.1109/JLT.2016.2628101
  * [13] Alvarado, A., Szczecinski, L., Fehenberger, T., et al.: ‘Improved soft-decision forward error correction via post-processing of mismatched log-likelihood ratios’, Proc. 2016 European Conference on Optical Communication (ECOC), Sep. 2016