Detecting global urban expansion over the last three decades using a fully convolutional network

Our study area is located from 65°S–75°N latitude worldwide, identical to the areal coverage of the NTL data (figure 1). This zone contains 0.13 billion km² of land surface area and accounts for 89.26% of the global land surface area. It includes all of Asia, Africa, Latin America and the Caribbean, and Oceania, and most of Europe and Northern America (figure 1).

Zoom In Zoom Out Reset image size

This study makes use of several data sources with global coverage, specifically the NTL, NDVI, and LST data mentioned above, to detect urban expansion. The 500 m Moderate Resolution Imaging Spectroradiometer (MODIS) Urban Land Cover (MOD500) data are used to select the training samples. Urban population statistics, Landsat images, and existing GHS data and the ESACCI global land cover data were also obtained.

The NTL data were obtained from the NOAA/NGDC (http://ngdc.noaa.gov/eog, accessed 20 September 2017). These data were collected by the Defense Meteorological Satellite Program's Operational Linescan System (DMSP-OLS) and include uncalibrated NTL data acquired in 1992 and radiance-calibrated NTL data from 1996, 2000, 2006, and 2010 at a resolution of 1 km. The monthly NTL data from 2016 at a resolution of 742 m acquired by the Suomi National Polar-orbiting Partnership's Visible Infrared Imaging Radiometer Suite (NPP-VIIRS) were also obtained. The approach introduced by Elvidge et al (2009) is applied to intercalibrate the uncalibrated NTL data collected in 1992 to improve their continuity and comparability, and the approach developed by Elvidge et al (2013) is performed to produce the annual mean value composite of the 2016 VIIRS NTL data. To reduce the influence of gas flares, we follow Zhou et al (2014, 2015) to exclude gas flare pixels from the NTL data using the gas flare mask obtained from NOAA/NGDC.

The NDVI data include the MODIS NDVI data for 2000, 2006, 2010, and 2016 and the NOAA-AVHRR NDVI data for 1992 and 1996 at a resolution of 1 km. The MODIS NDVI data were obtained from the MODIS 16 d composite products (http://ladsweb.nascom.nasa.gov, accessed 20 September 2017). The NOAA-AVHRR NDVI data were obtained from the USGS website (https://lta.cr.usgs.gov/NDVI, accessed 20 September 2017). Using the maximum value composite (MVC) approach introduced by Holben (1986), an annual MVC of the NDVI data for each year is produced to reduce cloud contamination.

The LST data in 2000, 2006, 2010, and 2016 at a resolution of 1 km are derived from the MODIS eight-day composite product (MOD11A2) (http://ladsweb.nascom.nasa.gov, accessed 20 September 2017). The nighttime LST data in this product are used because the LST in nighttime generally has better performance for distinguishing between urban and non-urban land (Buyantuyev and Wu 2010, Zakšek and Oštir 2012). Consistent with Mildrexler et al (2009), an annual MVC of the LST data for each year is also produced.

The MOD500 data used to select training samples were obtained from the Center for Sustainability and the Global Environment at the University of Wisconsin-Madison (http://www.sage.wisc.edu/, accessed 20 September 2017). The MOD500 data are produced through supervised classification of MODIS multi-spectral data circa 2001 (Schneider et al 2009, 2010). The data have an OA of 93%, which is a relatively accurate representation of the global urban land area (Schneider et al 2009, 2010). All of the NTL, NDVI, LST, and MOD500 data are resampled to a spatial resolution of 1 km. In addition, all of the data layers are geometrically corrected with reference to the MOD500 data.

In addition, the statistical data on the global urban population used for the accuracy assessment were obtained from the World Urbanization Prospects published by the United Nations Department of Economic and Social Affairs (2015) (https://esa.un.org/unpd/wup/, accessed 20 September 2017). The Landsat images were obtained from the USGS (http://glovis.usgs.gov.com, accessed 20 September 2017). The built-up grid and settlement grid of the GHS Layer were obtained from the European Commission Joint Research Center (http://ghslsys.jrc.ec.europa.eu/, accessed 20 September 2017). Finally, the ESACCI data were obtained from the European Space Agency (http://maps.elie.ucl.ac.be/CCI/viewer/index.php, accessed 20 September 2017).

Following the basic FCN framework presented by Long et al (2015), we develop a new FCN for detecting global urban expansion (figure 2). The FCN mainly includes one input layer, one pooling layer, three convolutional layers, one concatenation layer, and one output layer, and each convolutional layer follows an activation layer (figure 2(b)). The input layer is used to integrate the NTL, NDVI, and LST data to provide both socioeconomic and physical information. Longitudes and latitudes are also entered to provide information on the location of urban land. The convolutional and pooling layers are used to obtain feature information from the data at the pixel, neighborhood and urban region scales. The concatenation layer is used to integrate the feature information at the three scales to identify urban land (figure 2(b)).

Zoom In Zoom Out Reset image size

When extracting global urban land in 2000 and after, the input layers include five bands, specifically NTL, NDVI, LST, longitude, and latitude. However, when extracting urban land for 1992 and 1996, the input layer includes only four bands, specifically NTL, NDVI, longitude, and latitude. That is, the MODIS LST data are unavailable.

The three convolutional layers and the pooling layer employ different numbers and sizes of convolutional kernels to obtain information on urban land at three scales of pixel, neighborhood and urban region. Both the numbers and sizes of the convolutional kernels are determined following previous studies and trial-and-error procedures (LeCun et al 1990, Szegedy et al 2014, Long et al 2015).

The first convolutional layer contains 80 convolutional kernels with a size of 1 × 1 pixel and is used to obtain information at the pixel scale. Its relationship with the input layer can be represented by the equation (1):

$\begin{eqnarray}&&Con{v}_{i,j}^{1}=\displaystyle \sum _{I}^{{N}_{0}}({W}_{I,i}^{1}\cdot Inpu{t}_{I,j})+Bia{s}_{i}^{1},\end{eqnarray} \tag{ 1 }$

where $Con{v}_{i,j}^{1}$ and $Inpu{t}_{I,j}$ denote the DN values in the first convolution layer and the input layer respectively; i and I indicate the band; j indicates the pixel; ${W}_{I,i}^{1}$ is the weight of this convolution layer; $Bia{s}_{i}^{1}\,$ is the bias; and ${N}_{0}$ is the number of bands.

The second convolutional layer contains 40 convolutional kernels with a size of 3 × 3 pixels and is used to obtain information at the neighborhood scale. Its relationship with the input layer can be represented by the equation (2):

$\begin{eqnarray}&&Con{v}_{i,j}^{2}=\displaystyle \sum _{I}^{{N}_{0}}\displaystyle \sum _{d}^{{K}_{1}\cdot {K}_{1}}({W}_{I,i,d}^{2}\cdot Inpu{t}_{I,j,d})+Bia{s}_{i}^{2},\end{eqnarray} \tag{ 2 }$

where $Con{v}_{i,j}^{2}$ and $Inpu{t}_{I,j,d}$ denote the DN values in the second convolution layer and the input layer respectively; d indicates the pixel in the neighborhood of the jth pixel; ${W}_{I,i,d}^{2}$ is the weight of this convolution layer; $Bia{s}_{i}^{2}$ is the bias; and ${K}_{1}$ is the size of the neighborhood.

The third convolutional layer (which contains 20 convolutional kernels with a size of 1 × 1 pixel) and a pooling layer with a kernel size of 51 × 51 pixels are combined to obtain information at the scale of urban regions. The kernel size of 51 × 51 pixels used in the pooling layer is determined based on the largest urban region in the world (Angel et al 2016). The relationship between the third convolutional layer and the pooling layer can be represented by the equations (3), (4):

$\begin{eqnarray}&&Con{v}_{i,j}^{3}=\displaystyle \sum _{I}^{{N}_{0}}({W}_{I,i}^{3}\cdot Poo{l}_{I,j})+Bia{s}_{i}^{3},\end{eqnarray} \tag{ 3 }$

where $Con{v}_{i,j}^{3}$ and $Poo{l}_{I,j}$ denote the DN values in the third convolution layer and the pooling layer respectively; ${W}_{I,i}^{3}$ is the weight of this convolution layer; $Bia{s}_{i}^{3}$ is the bias

$\begin{eqnarray}&&Poo{l}_{I,j}=\displaystyle \sum _{g}^{{K}_{2}\cdot {K}_{2}}Inpu{t}_{I,j,g}/({K}_{2}\cdot {K}_{2}),\end{eqnarray} \tag{ 4 }$

where $Inpu{t}_{I,j,g}$ denotes the DN value of the gth pixel within the urban region of the jth pixel in the Ith band in the input layer; and ${K}_{2}$ denotes the size of the urban region.

The concatenation layer is used to integrate the information obtained by the three convolutional layers, and their relationships can be represented by the equation (5):

$\begin{eqnarray}&&Con{c}_{i,j}=\left\{\begin{array}{c}f(Con{v}_{i,j}^{1}),\,i\leqslant 80\\ f\left(Con{v}_{i-80,j}^{2}\right),80\lt i\leqslant 120\,\\ f(Con{v}_{i-120,j}^{3}),120\lt i\leqslant 140\end{array}\right.,\end{eqnarray} \tag{ 5 }$

where $Con{c}_{i,j}$ denotes the DN value in the concatenation layer, and $f(x)\,$ denotes the activation function, which can be represented by the equation (6):

$\begin{eqnarray}&&f\left(x\right)=\,\max (0,x).\end{eqnarray} \tag{ 6 }$

The relationship between the concatenation layer and the output layer can be represented by the equation (7):

$\begin{eqnarray}&&Urba{n}_{j}=\left\{\begin{array}{c}1,\,Outpu{t}_{j}^{urban}\gt Outpu{t}_{j}^{non \mbox{-} urban}\\ 0,\,Otherwise\,\end{array}\right.,\end{eqnarray} \tag{ 7 }$

where $Urba{n}_{j}$ denotes the class value in the output layer, in which 1 represents urban land whereas 0 represents non-urban land. $Outpu{t}_{j}^{urban}$ and $Outpu{t}_{j}^{non \mbox{-} urban}$ denote the probability that the jth pixel represents urban land or non-urban land respectively. These probabilities can be calculated using the following equations:

$\begin{eqnarray}\begin{array}{lll}Outpu{t}_{j}^{urban} & = & f\left(\displaystyle \sum _{i}^{140}\left({W}_{i}^{urban}\cdot Con{c}_{i,j}\right)\right.\\ & & +\,\left.Bia{s}^{urban}\Space{0ex}{3.0ex}{0ex}\right),\end{array}\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray}\begin{array}{lll}Outpu{t}_{j}^{non \mbox{-} urban} & = & f\left(\displaystyle \sum _{i}^{140}\left({W}_{i}^{non \mbox{-} urban}\cdot Con{c}_{i,j}\right)\right.\\ & & +\,\left.Bia{s}^{non \mbox{-} urban}\Space{0ex}{3.0ex}{0ex}\right),\end{array}\end{eqnarray} \tag{ 9 }$

where ${W}_{i}^{urban}$ and ${W}_{i}^{non \mbox{-} urban}$ are weights, and $Bia{s}^{urban}$ and $Bia{s}^{non \mbox{-} urban}$are biases. The FCN is developed using the Convolutional Architecture for Fast Feature Embedding (Caffe) software platform (Jia et al 2014).

The objective of calibration is to obtain the weights used in the FCN. Essentially, calibration involves iteratively estimating these weights using optimization algorithms with the goal of minimizing the loss function. This process includes three main steps: selecting training samples, determining the loss function, and setting the optimization algorithm and iteration parameters. First, we select training samples that represent urban and non-urban land from the MOD500 data and NTL data that cover the same period. The selection of training samples can be calculated using the equation (10):

$\begin{eqnarray}&&Trai{n}_{i}=\left\{\begin{array}{c}\,1,\,Mo{d}_{i}=1\,\\ 0,\,Mo{d}_{i}=0\,and\,PU{L}_{i}=1\\ -1,\,Otherwise\,\end{array}\right.,\end{eqnarray} \tag{ 10 }$

where $Trai{n}_{i}$ and $Mo{d}_{i}$ denote the class of the ith pixel in the training sample data and the MOD500 data respectively. A value of 1 denotes urban land, a value of 0 denotes non-urban land, and a value of −1 denotes the background, which is not included in the calibration to improve efficiency. $PU{L}_{i}$ denotes the class value in a potential urban land data, which takes values of 1 (potential urban land) or 0 (not potential urban land). The pixels with DN values greater than 1 nW cm⁻² sr⁻¹ in the NTL data were recognized as potential urban land (Dou et al 2017).

Further, we calculate the loss function after the output layer based on the equation (11) (Jia et al 2014):

$\begin{eqnarray}&&l\left(y,z\right)=\,\mathrm{log}\,\displaystyle {\sum }_{j}{e}^{{z}_{j}}-{z}_{y},\end{eqnarray} \tag{ 11 }$

where $z$ denotes the inputs to the loss function; $l(y,z)$ denotes the loss function of $z;$ ${z}_{j}$ denotes the result of estimating the jth class; and $y$ denotes the class of $z$ in the training samples.

Finally, following Long et al (2015), we utilize the stochastic gradient descent algorithm to perform the optimization. The iteration step size is 1 × 10^–7, and the maximum number of iterations is 1 million. The parameters of gamma and momentum both are 0.99. The weight decay is 0.0005.

First, we calibrate four FCNs with different remote sensing data using the Caffe deep learning platform so that the urban land at different times can be detected effectively (Jia et al 2014). Specifically, the first FCN applied for detecting urban land in 1992, is calibrated using the intercalibrated DMSP-OLS NTL data. The second FCN applied for detecting urban land in 1996, is calibrated on the basis of the radiance-calibrated DMSP-OLS NTL data. The third FCN applied for detecting urban land in 2000, 2006, and 2010, is calibrated on the basis of the radiance-calibrated DMSP-OLS NTL and LST data. The fourth FCN applied for detecting urban land in 2016, is calibrated on the basis of the NPP-VIIRS NTL and LST data. All of the FCNs involve NDVI, longitude, and latitude. We then use the calibrated FCNs to extract the global urban land for 1992, 1996, 2000, 2006, 2010, and 2016, respectively. We follow this step with the post-classification processing described below to ensure that the extracted global urban lands in different years are comparable. Because urban land has relatively high population densities, stronger NTL, and lower NDVI values, the areas of urban expansion generally show increased population densities and NTL and decreased NDVI values (Elvidge et al 2007, He et al 2014). In contrast, areas that have been transformed from urban land to non-urban land usually show decreased population densities and NTL and increased NDVI values. Therefore, we correct the extracted global urban land using the NTL, NDVI, and population densities from the HYDE dataset. The procedure can be represented by equation (12):

$\begin{eqnarray}&&Urba{n}_{i,j}^{Cor}=\left\{\begin{array}{c}1,\,\left(U{D}_{i,j}^{Pre}\ne -1{\rm{}}and{\rm{}}Urba{n}_{i,j}^{Pre}=1\right)\,or\,\left(U{D}_{i,j}^{Post}\ne 1\,and\,Urba{n}_{i,j}^{Post}=1\right)\\ Urba{n}_{i,j},\,Otherwise\,\end{array}\right.,\end{eqnarray} \tag{ 12 }$

where $Urba{n}_{i,j}^{Cor}$ denotes the value in the corrected urban land data; i indicates the year; j indicates the pixel; $Urba{n}_{i,j}^{Pre}$ and $Urba{n}_{i,j}^{Post}$ denote the values in the urban land data in the years before and after the ith year, respectively. The value of 1 denotes urban, while 0 denotes non-urban. $U{D}_{i,j}^{Pre}$ and $U{D}_{i,j}^{Post}$ denote the possible conversion to urban land of the jth pixel in the previous period and the latter period of the ith year. For these quantities, the value of 1 denotes the conversion from non-urban land to urban land; the value of −1 denotes the transformation of urban land into non-urban land; and the value of 0 denotes an uncertain conversion. $U{D}_{i,j}^{Pre}$ can be calculated using the following equation:

$\begin{eqnarray}&&U{D}_{i,j}^{Pre}=\left\{\begin{array}{c}\,1,\,Po{p}_{i,j}\gt Po{p}_{i,j}^{Pre}\,and\,{\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}\gt {\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}^{Pre}\,and\,{\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}\lt {\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}^{Pre}\\ -1,\,Po{p}_{i,j}\lt Po{p}_{i,j}^{Pre}\,and\,{\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}\lt {\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}^{Pre}\,and\,{\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}\gt {\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}^{Pre}\\ 0,\,Otherwise\,\end{array}\right.,\end{eqnarray} \tag{ 13 }$

where $Po{p}_{i,j},$ ${\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j},$ and ${\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}$ denote the population density, NTL, and NDVI corresponding to the ith year, respectively. $Po{p}_{i,j}^{Pre},$ ${\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}^{Pre},$ and ${\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}^{Pre}$ denote the population density, NTL, and NDVI in the year before the ith year, respectively. $U{D}_{i,j}^{Post}$ can be calculated using the following equation:

$\begin{eqnarray}&&U{D}_{i,j}^{Post}=\left\{\begin{array}{c}\,1,Po{p}_{i,j}\lt Po{p}_{i,j}^{Post}\,and\,{\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}\lt {\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}^{Post}\,and\,{\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}\gt {\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}^{Post}\\ -1,\,Po{p}_{i,j}\gt Po{p}_{i,j}^{Post}\,and\,{\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}\gt {\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}^{Post}\,and\,{\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}\lt {\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}^{Post}\\ 0,\,Otherwise\,\end{array}\right.,\end{eqnarray} \tag{ 14 }$

where $Po{p}_{i,j}^{Post},$ ${\rm{N}}{\rm{T}}{{\rm{L}}}_{i,j}^{Post},$ and ${\rm{N}}{\rm{D}}{\rm{V}}{{\rm{I}}}_{i,j}^{Post}$ denote the population density, NTL, and NDVI in the year after the ith year, respectively. Finally, we obtain information on the global urban expansion from 1992–2016.

We first evaluate the estimated urban expansion area resulting from this study at the national level using population census data (Sutton 2003, Zhang and Seto 2011). Using census data provided by the United Nations, we analyze the relationships between the estimated growth of urban land and the growth of urban population at the national level for 1992–2000, 2000–2010, 2010–2016, and 1992–2016. The estimated urban land expansion is strongly correlated (R² ≥ 0.5) with the census data at a statistically significant level of 0.01 (figure 3).

Zoom In Zoom Out Reset image size

Following the procedure presented by Potere et al (2009), we then use finer-resolution remote sensing data (Landsat imagery with a 30 m resolution) to evaluate the estimated urban expansion (figure 4). During this process, we select a total of 34 cities as samples that represent various levels of socioeconomic conditions in different regions using the stratified sampling approach introduced by Angel et al (2016) based on urban population and location (figure 1). The selected cities include two cities from group 1 (with populations greater than 10 million), two cities from group 2 (with populations between 1 million–10 million), and two cities from group 3 (with populations less than 1 million) located in each of the five following regions: Asia, Africa, Europe, Northern America, and Latin America and the Caribbean. We also selected two cities from group 2 and two cities from group 3 in Oceania, where no city has a population greater than 10 million (figure 1). We then extract the areas of urban expansion between 1992–2016 for the selected cities through visual interpretation of the Landsat images as reference data. In comparing the estimated urban expansion data obtained using the FCN and the reference data derived from the Landsat imagery, we calculate the OA, the kappa coefficient, the commission error (CE), and the omission error (OE) according to Olofsson et al (2014) (table 2, figure 4). The assessment results in the average values of 90.9%, 0.47%, 43.0%, and 49.9% for OA, kappa, CE, and OE, respectively (table 2, figure 4). Among the selected cities, Kigali in Rwanda displays the highest kappa value of 0.78 and is followed by Valledupar in Columbia and Malatya in Turkey, which display kappa values between 0.70 and 0.75 (table 2). Six cities, including Ho Chi Minh City in Vietnam, Medan in Indonesia, and Guadalajara in Mexico, display kappa values between 0.6–0.7 (table 2). Seven cities, including Mumbai in India; Portland, Oregon in the United States; and Cairo in Egypt, display kappa values between 0.5–0.6 (table 2).

Zoom In Zoom Out Reset image size

In addition, we compare the accuracy of our results with that of existing global urban expansion data with reference to Schneider et al (2010) and Zhou et al (2015). Among the six existing datasets that contain global urban expansion data (table 1), we first select the GHS built-up data, the GHS SMOD data, and the ESACCI data for comparison; these datasets have spatial resolutions and temporal coverage that are similar to those used in this study (tables 1–2, figure 4). We then assess the accuracy of these datasets using the reference data extracted from Landsat imagery for the selected cities and compare them with our results. This comparison shows that the accuracy of the results produced by this study is much higher than those of the existing datasets (table 2, figure 4). Specifically, the GHS built-up data display the average values of 90.1%, 0.20%, 49.7%, and 82.5% for OA, kappa, CE, and OE, respectively (table 2). The GHS SMOD data display the average values of 87.4%, 0.15%, 67.0%, and 83.3% for OA, kappa, CE, and OE, respectively (table 2). The ESACCI data display the average values of 90.5%, 0.19%, 49.7%, and 84.5% for OA, kappa, CE, and OE, respectively (table 2). The average OA and kappa of the results from this study are 0.4%–3.5% and 0.27–0.32 higher than those of these three global urban expansion datasets, respectively, whereas the average CE and OE obtained in this study are 6.7%–24.0% and 32.7%–34.6% lower (table 2).

The world has experienced rapid urban expansion. The global urban land increased from 274.7 thousand km² to 621.1 thousand km²; these values correspond to an urban expansion area of 346.4 thousand km² and an annual growth rate of 3.5% (figure 5, table 3).

Zoom In Zoom Out Reset image size

Asia displays the largest urban expansion area. Its urban land area has increased from 65.6 thousand km²–209.8 thousand km². These values correspond to an increase in urban area of 144.2 thousand km², which accounts for 41.6% of the global urban expansion area (figure 5(b), table 3). Among the countries in Asia, China displays the largest urban expansion area of 72.2 thousand km²; this value is 20.8% of the global urban expansion area. The urban areas in India and Turkey increased by 14.4 thousand km² and 7.4 thousand km², and these values account for 4.2% and 2.2% of the global urban expansion area, respectively. The urban expansion areas of five countries, Iran, Saudi Arabia, Vietnam, Indonesia, and Malaysia, contribute between 1%–2% of the global urban expansion area.

In Northern America, the urban land area increased by 69.7 thousand km² (from 78.9 thousand km²–148.5 thousand km²), and this increase accounts for 20.1% of the global urban expansion area (figure 5(b), table 3). In particular, the United States increased its urban land area by 66.4 thousand km², which accounts for 19.2% of the global urban expansion area. In Latin America and the Caribbean, the urban land area increased from 36.2 thousand km²–99.4 thousand km². This change corresponds to a net gain of 63.2 thousand km², which accounts for 18.2% of the global urban expansion area (figure 5(b), table 3). Within this region, Brazil, Argentina, and Mexico increased their urban land areas by 28.8 thousand km², 11.4 thousand km², and 7.9 thousand km² and account for 8.3%, 3.3%, and 2.3% of the global urban expansion area, respectively.

In Europe, the urban land area increased from 81.5 thousand km²–121.3 thousand km². This increase corresponds to an urban expansion area of 39.9 thousand km², which accounts for 11.5% of the worldwide growth (figure 5(b), table 3). Three countries in Europe, Italy, France, and Spain, increased their urban areas by between 5.0 thousand km²–5.5 thousand km², taking up approximately 1.5% of global urban expansion area. Russia expanded its urban areas by 3.6 thousand km², taking up 1.1% of global urban expansion area.

The areas of urban expansion in Africa and Oceania account for less than 10% of the worldwide growth (figure 5(b), table 3). In Africa, the urban land area grew from 7.7 thousand km² to 27.6 thousand km², an additional 19.9 thousand km², which accounts for 5.7% of the global growth. Within this region, South Africa increased its urban land area by 5.2 thousand km², which accounts for 1.5% of the global expansion. In Oceania, the urban land area increased from 4.8 thousand km²–14.5 thousand km². This change corresponds to an urban expansion area of 9.7 thousand km², which accounts for 2.8% of the worldwide growth. Within this region, Australia increased its urban land area by 8.4 thousand km², which corresponds to 2.4% of the global expansion.

Presently, many classification approaches are utilized to facilitate the detection of global urban expansion (Chen et al 2015, Pesaresi et al 2016, UCL-Geomatics 2017). For example, Chen et al (2015) detected global urban expansion between 2000–2010 using pixel- and object-based methods with knowledge (POK). Pesaresi et al (2013, 2016) used a symbolic machine learning method to detect global urban expansion from 1975–2015. UCL-Geomatics (2017) applied an unsupervised classification method to study global urban expansion from 1992–2015. However, effective methods to support the study of global urban expansion are still lacking. For example, the POK method requires large amounts of manual work and time. Unsupervised classification methods and symbolic machine learning methods ignore the features of urban land over multiple scales, such as shape, texture, and background information, resulting in the confusion of urban land and barren land. In addition, the lack of training samples of historical urban land areas around the world has led to relatively low accuracy and inconsistencies among the global urban land areas extracted using classification methods.

The proposed FCN has advantages in integrating remote sensing data from multiple sources, combining features over multiple scales, and making up for the lack of historical training samples. First, FCNs inherit the abilities of artificial neural networks and permits the adaptive processing of massive data from multiple sources (LeCun et al 2015). The FCN used in this study employs large numbers of neurons in the convolutional layers to obtain information from different data sources, including NTL, NDVI, LST, longitude, and latitude. The fusion of various types of data is made possible through the use of a nonlinear transformation function in the activation layers. The method then selects appropriate information to effectively characterize the socioeconomic and natural attributes of urban land and provides abundant and reliable information for urban land extraction (figures 2(b) and 6(a)). Second, the FCN uses three convolutional layers and one pooling layer to obtain urban land information at the pixel, neighborhood, and urban region scales. It effectively distinguishes urban built-up areas from barren land (figures 2(b), 6(b)). Third, when calibrated using training samples at one time, the FCN can be applied repeatedly due to its repetitive utilization of weights, which it inherits from deep learning models (Long et al 2015). This advantage solves the problems associated with a lack of historical training samples of urban land around the world but also improves the continuity and comparability of the extracted urban land at different times (figure 6(c)). Consequently, the proposed FCN can accurately extract urban land areas at one time and has the ability to integrate data from multiple sources and features over multiple scales (figures 6(a), (b)). For instance, the extraction of urban land in Beijing in 2016 performed using the FCN is associated with a kappa index of 0.79 and an OA of 95.1% (figures 6(a), (b)). Furthermore, we can effectively identify urban land at different times by repeatedly using the FCN. The detected urban expansion in Beijing from 1992 to 2016 is associated with a kappa index of 0.62 and an OA of 93.8% (figure 6(c)). In addition, the proposed FCN can automatically extract urban land areas from remote sensing data, which save large amounts of manual work and time. Thus, it provides an effective approach for detecting global urban expansion.

Zoom In Zoom Out Reset image size

Several long-term time series of surface reflectance data are available for the study of land use dynamics. These datasets include those collected by MERIS, AVHRR, PROBA-V, and SPOT-VGT and the high-resolution radar data obtained by TerraSAR-X, TanDEM-X, and Sentinel 1, as well as crowd-sourced data, such as OpenStreetMap (Esch et al 2017, UCL-Geomatics 2017). These data are obtained from multiple sources and have different spatial and temporal resolutions. The proposed FCN method can easily adapt to process these datasets. Thus, it has great potential to be widely used to detect global urban land expansion using time series data.

In this study, we produced a new global urban expansion dataset using the proposed FCN. This dataset has several merits over existing datasets (table 1). First, it provides adequate spatial information with 1 km resolution, which is much more precise than the HYDE dataset. Second, our dataset includes continuously comparable urban land data from 1992–2016, filling gaps in existing datasets of MGIS, GLC30 and IMPSA. Third, our dataset reveals higher accuracy than the GHS dataset, the FROM-GLC dataset and the ESACCI dataset.

In addition, the inconsistent definitions of urban land have been recognized as a primary reason for the discrepancies of existing global urban land datasets (Liu et al 2014). In this study, we followed Schneider et al (2009, 2010) and defined urban land as built-up area, i.e. the area dominated (more than 50% in cover) by non-vegetated, human-constructed elements, such as roads, buildings, runways, and industrial facilities (Liu et al 2014). Such definition makes our dataset comparable with most existing datasets with similar resolution (table 1).

However, our study also has some limitations. For example, the blooming effects of NTL data and the discontinuity between the DMSP-OLS data and the NPP-VIIRS data still result in relatively high OE and CE in some cities (e.g. Los Angeles, New York, and Victoria) (table 2).

In the future, we will attempt to adjust the number of convolutional layers and the size of the convolutional kernels of our FCN using high-resolution data (Long et al 2015, Fu et al 2017) to enable the detection of global urban expansion at resolutions finer than 1 km. We will also quantify the spatiotemporal patterns of global urbanization based on the results (Seto et al 2011, Liu et al 2016a). In addition, we will use the ecological models of Carnegie-Ames-Stanford approach, Integrated Valuation of Ecosystem Services and Tradeoffs, and Service Path Attribution Networks to evaluate the impacts of global urban expansion on ecosystems (McDonald et al 2008, 2013, Bren d'Amour et al 2017).