Mesh-Wise Prediction of Demographic Composition from Satellite Images
  Using Multi-Head Convolutional Neural Network
Population aging is one of the most serious problems in certain countries. In
order to implement its countermeasures, understanding its rapid progress is of
urgency with a granular resolution. However, a detailed and rigorous survey
with high frequency is not feasible due to the constraints of financial and
human resources. Nowadays, Deep Learning is prevalent for pattern recognition
with significant accuracy, with its application to remote sensing. This paper
proposes a multi-head Convolutional Neural Network model with transfer learning
from pre-trained ResNet50 for estimating mesh-wise demographics of Japan as one
of the most aged countries in the world, with satellite images from
Landsat-8/OLI and Suomi NPP/VIIRS-DNS as inputs and census demographics as
labels. The trained model was performed on a testing dataset with a test score
of at least 0.8914 in $\text{R}^2$ for all the demographic composition groups,
and the estimated demographic composition was generated and visualised for 2022
as a non-census year.
Yuta Sato
Mesh-Wise Prediction of Demographic Composition from Satellite Images Using Multi-Head Convolutional Neural Network

###### Abstract

Population aging is one of the most serious problems in certain countries. In order to implement its countermeasures, understanding its rapid progress is of urgency with a granular resolution. However, a detailed and rigorous survey with high frequency is not feasible due to the constraints of financial and human resources. Nowadays, Deep Learning is prevalent for pattern recognition with significant accuracy, with its application to remote sensing. This paper proposes a multi-head Convolutional Neural Network model with transfer learning from pre-trained ResNet50 for estimating mesh-wise demographics of Japan as one of the most aged countries in the world, with satellite images from Landsat-8/OLI and Suomi NPP/VIIRS-DNS as inputs and census demographics as labels. The trained model was performed on a testing dataset with a test score of at least 0.8914 in R\({}^{2}\) for all the demographic composition groups, and the estimated demographic composition was generated and visualised for 2022 as a non-census year.

 demographic composition estimation, population aging, remote sensing, deep learning, convolutional neural network, transfer learning

## I Introduction

While the total population of human beings has increased radically, some countries and regions are being faced with unprecedented aging of society. As a forecast, 1 in 6 people in the world will be more than 60 years old, whose number will be doubled from 1 billion in 2020 to 2.1 billion in 2050 [1]. The radical change of demographic composition in countries and regions would affect the feasibility of social welfare and urban planning policies of governments such as pensions, medical care, transportation, and other physical and social infrastructure [2, 3]. In order to adjust to the transition of societies, obtaining up-to-date demographic composition is of urgent importance. As seen in the case of Italy [4], a low fertility rate stems from multiple factors including identity towards family, employability, and composition of family, but the importance of each factor highly varies depending on sub-regions. In terms of making data-driven decisions on the population aging with high geographic heterogeneity, it is not sufficient to only obtain the total population, but instead crucial to gain metrics of the demographic composition. One of the most accurate and detailed data sources for demographic composition is a census of each country. Although the granularity varies depending on countries and regions, a census provides precise attributes of residents in each tract or mesh at some timing of "census year", such as population by age, number of households, and average income in the area. However, due to its thorough methodology and financial or administrative constraints, a census is usually conducted every several years, such as ten years in the U.S. and the U.K. [5, 6]. Thus, it is not feasible to take place similar national surveys more frequently, such as on an annual basis, and approximation is required by estimation from other data sources.

Among others, Japan is now being suffered from the unprecedented population aging in the history of human beings. In estimation [7], the ratio of people with ages over 64 years old occupied 30.0 percent of the whole population in Japan, which is the second highest among all the other countries next to Monaco. In parallel, the polarisation of the population is gaining momentum with Tokyo Metropolitan Area enjoying the population influx from the other regions since the 1970s [8]. Under the circumstances, there are an increasing number of "marginal settlements" (_genkai sharaku_) in the rural areas, where lives became infeasible due to the lack of young labours and/or infrastructure and the residents would have to abandon their housings to migrate to the urban area [9]. Considering the fast pace of population aging with an estimation of the population ratio over 64 years old as 38.8% by 2050 [10], it is crucial to understand the up-to-date demographic composition in a granular manner, from the perspective of urban and regional planning. Although the Statistics Bureau of Japan is providing the mesh-wise demographic composition as the result of the census, the interval of its provision is every five years [11]. To fulfill the gap between the census years, precise estimation is of importance.

Nowadays, pattern recognition of non-linear relationships between different types of datasets is prevalent with potential application to problems in geography, due to the increase of data volume and new machine learning algorithms proposed [12]. In particular, Deep Learning models as extensions of ANN have shown remarkable performances in the area of remote sensing, utilising satellite images (i.e. land-use classification [13]). Since some well-known satellites such as Landsat-8 [14], Sentinel-2 [15], and Suomi NPP [16] periodically observe the surface of the Earth globally, they constantly provide a rich amount of data source possibly enough to approximate the other socio-economic variables. When it comes to the estimation of population, previous studies have already tried to estimate the population density from satellite images using Convolutional Neural Network (CNN) such as
in India [17] and China [18, 19]. However, the estimation of demographic composition rather than the total population is understudied, especially in the context of population aging. This is mainly because the estimation of the population tended to focus on capturing the fast pace of urbanisation around cities in 'developing' countries, where rigid data source like census is radically not available. Moreover, some models from the previous studies relied on a complex mixture of multiple data sources, including POI data (i.e. [18]). Although better prediction scores could be expected as the variety and amount of input data are increased, this approach implicitly lowers the generalizability of models because of data availability depending on the area of interest. Thus, it is important to develop prediction models by intentionally limiting the input data source only to those available commonly among countries and regions which are (and will be) suffering from population aging.

Being aware of the issues mentioned above, this paper proposes a CNN model of predicting the demographic composition of Japan on an annual basis, from the satellite images of Landsat-8/OLI and Suomi NPP/VIIRS-DNS, which are publicly available at any location on the Earth. Since the mesh-wise census demographic composition is available in Japan every five years, this paper employs demographic composition from the census conducted in 2015 and 2020 as labels for training, validation, and testing. At the end of the paper, a new dataset is synthesised from the satellite images in 2022, which is not a census year. Due to the public data availability in any country and region, this paper shall be significant in the sense that it has generalisability of granular estimation of demographic composition not only in Japan but also the other countries and regions.

## II Related Works

### _Remote sensing for demographics estimation before Deep Learning_

In the field of remote sensing, satellite images have been utilised for a variety of purposes such as tracking descertification, aerosol pollution, and urbanisation [20]. Since its first launch in 1972, Landsat Missions operated by the U.S. Geological Survey (USGS) and National Aeronautics and Space Administration (NASA) have been employed as a major source of surface reflectance of the whole globe [14]. Although the history is relatively shallower, Sentinel Mission from European Space Agency (ESA) followed Landsat as another source of satellite images since the first launch in 2014 with several applications including cropland monitoring, forest ecology, and urbanisation [15, 21]. When it comes to demographics estimation, the urban population of Ohio in the U.S. was estimated from the impervious surface pixel value of Landsat-7 Enhanced Thematic Mapper (ETM+) images, by Ordinary Least Square (OLS) regression on a pixel basis and Spatial Autoregressive model (SAR) on a zonal basis [22]. However, as the authors stated, the accuracy of the zonal estimation was "unacceptable" due to the lack of capturing non-linearity in the pattern of pixels with a linear regression model.

While Landsat and Sentinel are observing the earth's surface reflectance in the daytime, other satellites focus on night-time light (NTL). With the global coverage, the Defence Meteorological Satellite Program's Operational Linescan System (DMSP-OLS) by the U.S. Department of Defence and Suomi National Polar-Orbiting Partnership's Visible Infrared Imaging Radiometer Suite Day-Night Band (Suomi-NPP/VIIRS-DNB) are ones of the major sources for NTL [23]. Previous studies pointed out that NTL could approximate socio-economic indicators such as land-use boundary [24] and GDP [25, 26]. For the estimation of population, the correlations between NTL from DMSP-OLS and 1km mesh-wise population density were clarified in Hokkaido Prefecture, Japan by geographically weighted regression (GWR), with a prediction score of 0.8833 in \(R^{2}\)[27]. In another study [28], 100m grid-wise population density was estimated in the extent of China using multiple sources including NTL from Suomi-NPP/VIIRS-DNB, Normalized Difference Vegetation Index (NDVI) from Moderate Resolution Imaging Spectroradiometer (MODIS), Digital Elevation Model (DEM), and the other POI data by Random Forest regression, whose accuracy score varied \(R^{2}\) between 0.3 and 0.8 depending on provinces.

### _Deep Learning and transfer learning with application to remote sensing_

The foundational concept of Artificial Neural Network (ANN) was initially proposed as a model of mimicking the human brain's biological neurons as multi-layer perceptron (MLP) in the sense that input information propagates in a non-linear manner for pattern recognition of objects such as images, sounds, and sentences [29]. ANN introduced an activation function for each neuron, which enables the model to capture the non-linearity of input data sources. ANN is composed of multiple layers of neurons with weights, and each of the weights is trained through a chain of partial derivatives on a loss function called backpropagation. Although the concept was established earlier, it was computationally expensive at the moment when Graphics Processing Unit (GPU) was not adjusted for the computation of backpropagation [30]. With the advance of GPU, an ANN model called AlexNet was proposed for image recognition tasks with the introduction of a convolution layer, which captures patterns of neighbouring pixels efficiently by sharing parameters as filters [31]. This 8-layer CNN model scored a remarkable top-5 classification error of 15.3% on the ImageNet dataset consisting of 1.2 million 224 \(\times\) 224 pixels with red-green-blue (RGB) colours with 1,000 image categories. AlexNet also introduced a new activation function called Rectified Linear Unit (ReLU) as \(\sigma(x)=\max(0,x)\) for inputs \(x\). Compared with previous activations functions such as sigmoid and tanh, ReLU performed better for capturing non-linearity in the CNN model. Since the success of AlexNet, a variety of CNN models have been proposed with increased accuracy and the number of convolution layers such as GoogleNet with 22 layers [32] and VGGNet with up to 19 layers [33]. However, the increase of layer numbers did not simply reflect on better prediction, as the gradient in backpropagation would be unstable to be explosive or vanished. To solve this problem and make the layers "deeper", ResNet was proposed by adopting residual layers
which introduced skip connections from previous convolution layers [34]. With this improvement, ResNet could augment layers up to 152 layers, with 3.57% top-5 classification error on the ImageNet dataset, which outperformed the error of human beings with 5.1% [35].

As the accuracy of CNN models was improved, well-known pre-trained CNN models were fine-tuned and re-utilised for multiple domains other than the image classification of the ImageNet dataset, namely transfer learning [36]. While default CNN models require a huge amount of datasets (i.e. 1,000 categories classification of 1.2 million images from ImageNet), transferred ones can achieve a considerable level of prediction with fewer samples if the task domain is similar. One of the simplest ways of transfer learning is to append or modify the last layer of models with holding the others the same at the beginning of training. For instance, ImageNet-based CNN models including AlexNet, GoogLeNet, and VGGNet could be transferred via fine-tuning into two medical image detection tasks from CT-scanned slices [37].

In the field of remote sensing, CNN models have been mainly used for image segmentation tasks, where the pixel-wise classification output size is the same as the one from input images, such as cloud detection [38], deforestation monitoring [39], and land-use classification [40, 41]. In the case of demographics estimation, previous studies proposed CNN models for regression tasks of population density. As an example, ImageNet-pre-trained CNN models including VGGNet and ResNet were employed for estimating the total population in 30 arc-second grids (approximately 1 km on the equator) in the U.S., using Sentinel-2 images in RGB and near-infrared (NIR) [19]. The test error of logged population in the Metro Dallas area scored \(R^{2}\) of 0.875 and 0.906 for VGGNet and ResNet, respectively. Others developed an original 6-layer CNN model with feature extraction of MLP estimated the 100m mesh-wise population density in Shenzhen, China, using multiple data sources including NTL from VIIRS-DNS and NDVI from the Chinese Academy of Sciences, and social sensing POI data [18]. Although daytime multispectral images like Landsat and Sentinel were not used, the prediction on the testing dataset counted 0.77 in \(R^{2}\). While these previous studies showed the success of CNN models in the task of population estimation, its more detailed breakdown into the demographic composition is understudied. Also, some of the data sources employed included local dataset which was only available in specific countries or regions, causing the limitation in generalisability of trained models into other geographic locations.

## III Datasets and study areas

To proceed with the development of a CNN model for estimating demographic composition in Japan, this paper employed Landsat-8/OLI Collection 2 Level-2 [42] and Suomi NPP/VIIRS-DNS [16] as inputs, and mesh-wise population census data provided by Statistics Bureau of Japan via e-Stat database [43] as ground truth labels. Samples were obtained from 2015 and 2020 with labels of demographic composition and 2022 without labels as a non-census year for generating estimated labels. The geographic extent of the samples was defined as the whole meshes available in Japan either from the census of 2015 or 2020. Each of the mesh cells was geographically defined as "Basic Grid Square" with a geodesic system of WGS 84, which spans 30 arc-seconds and 45 arc-seconds (approximately 1km square in the territory of Japan) for latitude and longitude, respectively [44]. While the population census provides the number of residents in the mesh cell by age as of October 1st in the census year, aggregated number of three demographic composition groups of residents below 14 years old, between 15 and 64 years old, and over 64 years old were used as labels for a CNN model. This paper utilised satellite images both from daytime (Landsat-8/OLI) and night-time (Suomi NPP/VIIRS-DNS) so that the input could include more detailed information about land-use and human activities. Assumedly, the existence of only either building structures or NTL may not necessarily indicate the settlement of people (i.e. abandoned housings in rural areas; 24-hour factories in industrial areas). Thus, stacking multiple images as a single tensor could support a CNN model to recognise the demographic composition in the area from multiple aspects in the inputs. Landsat-8/OLI consists of nine multispectral images depending on the wavelength of light. Out of the nine bands, band 4, 3, and 2 were used as RGB images of natural colour. In addition to it, band 7, 6, and 5 were also used for another RGB image as false colour, so that wavelengths out of human eyes could also be recognised by a CNN model. Moreover, three synthesised indices were appended as the third RGB image from the original Landsat/OLI multispectral images as follows; NDVI, Normalized Difference Built-up Index (NDBI), and Normalized Difference Water Index (NDWI). While NDVI as an index for vegetation was already used for population estimation (i.e. [28]), this paper additionally employed the indices for buildings and water to explicitly delineate the different characteristics of lands. By referring to [45, 46], these three indices were calculated as follows:

\[\text{NDBI}=\frac{\text{MIR}-\text{NIR}}{\text{MIR}+\text{NIR}}=\frac{\text {Band 6}-\text{Band 5}}{\text{Band 6}+\text{Band 5}} \tag{1}\]

\[\text{NDVI}=\frac{\text{NIR}-\text{Red}}{\text{NIR}+\text{Red}}=\frac{\text {Band 5}-\text{Band 4}}{\text{Band 5}+\text{Band 4}} \tag{2}\]

\[\text{NDWI}=\frac{\text{Green}-\text{NIR}}{\text{Green}+\text{NIR}}=\frac{ \text{Band 3}-\text{Band 5}}{\text{Band 3}+\text{Band 5}} \tag{3}\]

where NIR and MIR stand for near-infra-red and middle-infra-red, respectively.

Figure 1 illustrates the overview of processing images from the three data sources for a CNN model. Since satellite images are periodically collected from the space as snapshots of the Earth, it is inevitable to encounter cloud coverage in the images. To avoid biases from cloud coverage, this paper employed adjusted images where each pixel is an annual average of those without clouds. For Landsat-8/OLI, Google Earth Engine API was used for obtaining the adjusted images of "USGS Landsat 8 Level 2, Collection 2, Tier 1" [47]. In the case of Suomi NPP/VIIRS-DNS, adjusted raster
files were acquired from Annual VNL V2.1 maintained by the Earth Observation Group [48, 49, 50]. The spatial resolutions of Landsat-8/OLI and Suomi NPP/VIIRS-DNS are approximately 30m and 750m, respectively. To fulfill the gap, the pixels of Suomi NPP/VIIRS-DNS were adjusted to those from Landsat-8/OLI. Due to the coarse spatial resolution of satellite images for the size of original meshes from the census, mesh cells were firstly aggregated from 30 arc-seconds \(\times\) 45 arc-seconds (1km square) to 5 arc-minutes \(\times\) 7 arc-minutes and 30 arc-seconds (10km square). After the cloud coverage removal by annual averaging, two types of satellite images were cropped by the aggregated mesh cells. The cropped images were pre-processed for resizing the shape into 334 \(\times\) 334 pixels square so that the following CNN model can deal with the images as inputs with 30m spatial resolution for the 10km square extent. In the pre-processing, the values of Suomi NPP/VIIRS-DNS were re-scaled with a natural logarithm. The output of this flow was the stacked images of Band 2 to 7 from Landsat-8/OLI, their synthesised three indices (NDBI, NDVI, NDWI), and the NTL from Suomi NPP/VIIRS-DNS.

Based on the processing per sample explained in Figure 1, Figure 2 illustrates the flow of sample management as a framework of a supervised learning algorithm. In order to train a CNN model for the task of demographic composition

Fig. 1: Flow chart of data processing per one single sample.

Fig. 2: Flow chart of data splitting for the supervised learning framework.
estimation, all the processed images and labels from the census year of 2015 and 2020 were split into those for training, validation, and testing. In order to ensure even distribution among the dataset groups, stratified sampling was employed by mesh cell ID and total population. In other words, unique values of strings for mesh cell IDs were extracted from all the samples, and averaged values of the total population between 2015 and 2020 were obtained. It is noteworthy to mention that there was no geographical duplication between the training, validation, and testing dataset since it is already rare to observe a radical urbanisation process in Japan and geographical duplication might cause data leakage problems [51] in the supervised learning framework. The mesh cell IDs were assigned to 10 partitions by ten deciles of the averaged total population in the histogram. Out of the 10 partitions, samples were evenly split into training, validation, and testing with a ratio of 8:1:1. The split samples were imported to PyTorch Dataloader objects [52] for mini-batch learning of a CNN model. Images from 2022 were directly imported to PyTorch Dataloader objects without labels. Figure 3 shows the spatial distribution of samples from the training, validation, and testing dataset from the years 2015 and 2020. For the mesh cells where there was no data was omitted from the output since it does not necessarily mean the non-existence of residents, but also implies the lack of conducting a population census (i.e. Japanese government is claiming the sovereignty of the Kuril Islands without conducting a population census since they are under the effective control of Russian government). Table I describes the summary statistics of samples by the type of dataset and sample year. As seen on the left in Figure 4, the distribution of the population is heavily skewed. On the other hand, the logarithm with base 10 was applied for the original number of the population, which is shown on the right in Figure 4. Thus, for the labels of the CNN model, this paper employs the logarithm of the original population with base 10, However, some samples contain 0 for certain demographic composition groups, where the logarithm is not defined. To avoid errors in the following computation, 0 was replaced with 1 (\(10^{0}\)) in advance. Figure 5 shows three randomly sampled images and labels. As seen in the figure, the images with RGB of Band 4, 3, and 2 were close to natural colour, while those with Band 7, 6, 5 and NDBI, NDVI, and NDWI were expressed in false colour. The combination of NDBI, NDVI, and NDWI could explicitly delineate water, forest, and others.

Fig. 3: Spatial distribution of training, validation, and test dataset by two census years. JGD2011 was applied for the coordinate reference system. The map is based on OpenStreetMap.
NTL from Suomi NPP/VIIRS-DNS highlighted the location of human activities at night.

## IV Methodology

Let the CNN model \(f:=\mathcal{X}\rightarrow\mathcal{Y}\), where \(\mathcal{X}:=\{\mathbf{X}_{i}|\mathbf{X}_{i}\in\mathbb{R}^{334\times 334\times 1 2}\}\) a feature tensor space consisting of four types of RGB-based images and \(\mathcal{Y}:=\{\mathbf{y}_{i}|\mathbf{y}_{i}\in\mathbb{R}^{3}\}\) a label vector space consisting of three number of demographic composition groups. Then, the predicted demographic composition vector \(\hat{\mathbf{y}}_{it}\in\mathcal{Y}\) is expressed as following:

\[\hat{\mathbf{y}}_{it}=f(\mathbf{X}_{i};\mathbf{\Theta}_{t}) \tag{4}\]

where \(\mathbf{\Theta}_{t}\) is an element of \(\mathbf{\Theta}\) as a set of layer weights in \(f\) such that \(\mathbf{\Theta}=\{\mathbf{\Theta}_{t}\}\) and \(t\in\mathbb{N}\) denotes the order of training. The CNN model \(f\) was trained through the gradient descent method as backpropagation such that

\[\boldsymbol{\theta}_{t}^{l}=\boldsymbol{\theta}_{t-1}^{l}-\text{ Optim}(\alpha,\nabla\boldsymbol{\theta}_{t-1}^{l}) \tag{5}\]

\[\nabla\boldsymbol{\theta}_{t-1}^{l}=\frac{\partial\mathcal{L}(\hat{\mathbf{y} }_{it-1},\mathbf{y}_{i})}{\partial\boldsymbol{\theta}_{t-1}^{l}} \tag{6}\]

where \(\boldsymbol{\theta}_{t}^{l}\in\mathbf{\Theta}_{t}\) and \(\text{Optim}(\alpha,\nabla\theta_{t}^{l})\) denote the \(l\)-th layer of \(\mathbf{\Theta}_{t}\) and the optimiser with the learning rate \(\alpha\in\mathbb{R}^{+}\) and gradient \(\nabla\theta_{t}^{l}\), respectively. For the loss function \(\mathcal{L}(\hat{\mathbf{y}}_{it},\mathbf{y}_{i})\), this paper adopts the formula of following:

\[\mathcal{L}(\hat{\mathbf{y}}_{it},\mathbf{y}_{i})=\frac{1}{N}\sum_{i}^{N}|| \hat{\mathbf{y}}_{it}-\mathbf{y}_{i}|| \tag{7}\]

where \(N\in\mathbb{N}\) denotes the number of samples in a dataset. For the optimiser, Adam [53] was applied with following:

\[\text{Optim}(\alpha,\nabla\boldsymbol{\theta}_{t}^{l};\beta_{1},\beta_{2})= \alpha\frac{\hat{m}_{t}}{\sqrt{v_{t}+\epsilon}}\nabla\boldsymbol{\theta}_{t-1} ^{l} \tag{8}\]

\[m_{t}=\beta_{1}m_{t-1}-(1-\beta_{1})\nabla\boldsymbol{\theta}_{t-1}^{l} \tag{9}\]

\[v_{t}=\beta_{2}v_{t-1}-(1-\beta_{2})(\nabla\boldsymbol{\theta}_{t-1}^{l})^{2} \tag{10}\]

\[\hat{m}_{t}=\frac{m_{t}}{1-\beta_{1}^{l}} \tag{11}\]

\[\hat{v}_{t}=\frac{v_{t}}{1-\beta_{2}^{l}} \tag{12}\]

where \(\beta_{1}\) and \(\beta_{2}\) denote the hyperparameters of the first and second moment vectors. After every training of all samples in the training dataset as one epoch, the trained CNN model \(f\) was validated with samples in the validation dataset by calculating \(\mathcal{L}\). This one cycle of training and validation was iterated over 250 epochs, and a model with the least validation error was selected as the best model in terms of generalisability. The chosen best model was tested by samples in the testing dataset

Fig. 4: Distribution of population in all the samples from the training, validation, and testing dataset by demographic composition groups. Scales on the y-axis are original and logarithm with base 10 for the left and right figure, respectively. The blue lines across ages indicate that the values were from the same sample in a mesh cell either in 2015 or 2020.
and \(R^{2}_{(g)}\) was calculated for all the demographic composition groups \(g\) as follows:

\[R^{2}_{(g)}=\frac{\sum(\hat{y}^{(g)}_{i}-\bar{y}^{(g)}_{i})^{2}}{\sum(y^{(g)}_{i} -\bar{y}^{(g)}_{i})^{2}} \tag{13}\]

Figure 6 displays the architecture of the CNN model \(f\). The model employed multi-head CNN based on ResNet50 [34], so that multi-layered images could be processed as multiple RGB-coloured tensors. Four images were set as inputs of RGB-coloured tensors, where NTL was duplicately stacked for each colour. These inputs were passed through the first convolution layer and concatenated along the dimension of channels. The concatenated tensor was convoluted with 1 \(\times\) 1 filter and aggregated by max pooling layer for down-sizing, which only keeps a maximum value in the window. The pooled tensor was further convoluted 48 times with skip connections. The output from all the convolutions was pooled by 1 \(\times\) 1 average pooling layer as an operation of taking an average in the window, leading to 1 \(\times\) 1 \(\times\) 2048 tensor. The tensor was flattened to a 2048-dimension vector and passed through two fully-connected layers with the final output of a 3-dimension vector \(\hat{\mathbf{y}}_{it}\), estimating the demographic composition in the sample. The output size of each convolution and pooling layer was described in Figure 6 by the formula of following:

\[n_{\text{out}}=\lfloor\frac{n_{\text{in}}+2p-k}{s}\rfloor+1 \tag{14}\]

where \(n_{\text{in}}\) and \(n_{\text{out}}\) denote the size of input and output, respectively. \(p\), \(k\), and \(s\) stand for the size of zero-padding, filter, and strides, respectively. \(\lfloor.\rfloor\) operation indicates obtaining an integer number from the float number. For instance, the output size of the first convolution for Landsat-8/OLI Band 4,3, and 2 was \(\lfloor\frac{334+2*3-7}{2}\rfloor+1=167\). As shown in Figure 6, all the parameters of this CNN model were transferred from ResNet50 being pre-trained by the ImageNet dataset, except the second convolution layer and the last fully-connected layer. Moreover, batch normalisation was implemented after each of the convolution layers, where each batch in the forward propagation was normalised with the average and variance, so that the training process could be stabilised [54]. ReLU was used as an activation function for all the convolution layers after batch normalisation. In addition to batch normalisation, dropout [31] was implemented for avoiding overfitting. Dropout is an operation for removing a certain amount of neurons within a layer during training so that the CNN model would be robust to the heterogeneity in the dataset as if it behaves as ensemble learning. However, as a previous study reported [55], the co-existence of dropout and batch normalisation in the convolution layer would worsen the performance of CNN models due to variance shifts. Thus, dropout was applied only to the last two fully-connected layers, where batch normalisation was not implemented in the original ResNet50. Since all the labels were biased to \(10^{0}\) and their values were always non-negative, the last layer was also passed through ReLU. For the hyperparameters of this machine learning framework, the learning rate \(\alpha\) and coefficients of the first and second moment \(\beta_{1}\) and \(\beta_{2}\) for Adam were set as 0.001, 0.9, and 0.999, respectively. The ratio of dropout was set as 0.25 for both of the two fully-connected layers. As an environment of model training, AWS Sagemaker Studio was used with an instance type of g4dn.xlarge consisting of 4 vCPUs, 16 GB memory, and 1 GPU of NVIDIA Tesla T4 [56].

## V Result

As a result of the training and validation process, Figure 7 illustrates the training and validation loss function over 250

Fig. 5: Images and labels from three random samples from the training dataset.
epochs. In the early epochs, the validation loss function had a volatility of up to 1.2220. While the training loss steadily decreased throughout the epochs, validation loss reached its lower bound of 0.35 at around epoch 50 and little improvement was observed in the following epochs. The training loss recorded 0.3186 as the minima at epoch 238, and the model at the epoch was employed as the one with the most generalisability for the testing and dataset generation.

With the trained model chosen above, the testing dataset was utilised for calculating the testing score in \(R^{2}_{(g)}\), as well as the loss function. To do so, all the images in the testing dataset were passed through the CNN model as forward propagation, and the outputs were used for calculating the loss function as a test error. As seen in Table I, the loss function in the testing dataset was 0.3215. This value was close to 0.3186 in the validation dataset, while 0.1032 in the training dataset

Fig. 6: Architecture of ResNet50-based multi-head CNN model.

Fig. 7: Training and validation progress of the CNN model. The grey horizontal and vertical dotted lines indicate the lowest value of the loss function in the validation dataset throughout the epochs and its epoch, respectively. The intersection of the dotted lines indicates the loss function of the validation dataset at epoch 238.
was relatively smaller than the other two. The outputs obtained were matched with the labels in the testing dataset and \(R^{2}_{(g)}\) was calculated for all the demographic composition groups \(g\). As Table II and Figure 8 illustrate, the \(R^{2}_{(g)}\) reached at least 0.8914 for all the demographic composition groups. Although mesh cells with a high population had slight overestimation, the regression lines were aligned along the diagonal of references. Based on the trained CNN model, Figure 9 shows the estimated ratio of the population over 64 years old in 2022, which is not a census year. The mesh cells were employed from those where the population was counted either in 2015 or 2020.

## VI Discussion

As population aging is in progress among an increasing number of countries and regions, prompt understanding of demographic composition is of urgent importance. This paper demonstrated the mesh-wise estimation of demographic composition in Japan as one of the countries with the most aged society in the world, from the dataset of population census in 2015 and 2020 as labels and satellite images of Landsat-8/OLI and Suomi NPP/VIIRS-DNS. The chosen ResNet50-based CNN model after the training and validation process was tested by the testing dataset with metrics of \(R^{2}_{(g)}\) between original and estimated values for each of the demographic composition groups. In the end, the trained model was used for estimating the demographic composition of Japan in 2022 as a non-census year, and the aged population rate was visualised. The result showed a remarkable test score of \(R^{2}_{(g)}\) with at least 0.8914 for all the demographic composition groups, indicating three contributions of following; Firstly, the result displays the effective usage of pre-trained ResNet50 [34] from the ImageNet dataset consisting of more than 1.2 million images as a 1,000-objects classification task for the demographic composition estimation task with only 6,724 training samples. This result corresponds to previous studies (i.e. [19]) in the sense that the remote sensing tasks share a similar domain with the ImageNet object classification task as image recognition. This performance was seen even with the introduction of multi-head inputs in the first convolution layer, where multiple RGB-colour images were accommodated. Secondly, this paper proposes the feasibility of estimating demographic composition rather than total population from satellite images, so that the progress of population aging could be observed without conducting formal surveys such as census. As the values of the loss function both for the validation and testing dataset were not deviated from around 0.32, it is reasonable to confirm that the model held generalisability enough to estimate other geographical locations which were not included in this dataset. Lastly, the proposed CNN model only employed publicly available satellite images as inputs for enduring generalisability in any other geographic location. Since Landsat-8/OLI and Suomi NPP/VIIRS-DNS cover the surface of the earth holistically and periodically, the trained CNN model could be further utilised in other countries in any year. Considering the performance and data availability, the proposed model could help fulfilling the gap between census years not only in Japan but also the other countries and regions for the purposes of urban planning policymaking, especially where rigorous population surveys are limited due to a lack of financial and human resources.

However, there are certain limitations in the proposed CNN model when it comes to enduring generalisability for further usage. First and foremost, biases in the dataset could result in overfitting in other datasets out of this paper. In terms of the spatial biases stemming from the dataset, the proposed CNN model was trained only by the satellite images in the realm of the population census held by the Japanese government. Since this paper did not conduct testing by a dataset from other countries and regions, the generalisability of the trained CNN model is still uncertain if the condition of morphological landscape in satellite images is different from the one in Japan. To consider the temporal biases in the proposed CNN model, the dataset was obtained from the years 2015 and 2020, where each sample was independently treated and did not explicitly consider the temporal relationship in the same geographical location. Also, the satellite images were pre-processed by annually averaging pixel values from those without clouds, which lacked temporal transition within a sample year. From the model architecture point of view, one way to enhance the proposed CNN model for overcoming these problems is

Fig. 8: Prediction performance of the CNN model in the testing dataset by demographic composition groups. The real lines were drawn from the linear regression of the predicted values on the original ones. The dotted lines were drawn as references which indicate the symmetric normal lines which could have been realised if all the samples would have been predicted precisely.
to implement spatio-temporal CNN with the introduction of Long Short-Term Memory (LSTM) blocks [57] or Attention mechanism [58]. Moreover, rather simply and naturally, data argumentation by expanding the geographical locations and sample years could improve the proposed CNN model. Another limitation is the precision of the proposed CNN model. Although the test score reached 0.8914 for the \(R_{(g)}^{2}\) of all the demographic composition groups \(g\) in the testing dataset, Figure 8 indicates "heteroscedasty" in the correlation between the original and estimated values. In particular, samples with small population tended to be overestimated with poor precision and higher variance than the others. This could be due to the spatial resolution of the satellite images in the input, where 30m for Landsat-8/OLI and 750m for Suomi NPP/VIIRS-DNS were not granular enough to detect the shape of individual buildings in the extent. Not surprisingly, replacement of Landsat-8/OLI and Suomi NPP/VIIRS-DNS with finer satellite images and stacking of other data sources like spatially interpolated POI data would solve this problem. However, this approach would limit the data availability for geographical locations where the data source cannot be obtained. One potential method without employing new data sources is to synthesise satellite images with finer spatial resolution using Generative Adversarial Network (GAN) for super-resolution [59]. In the future study, the proposed CNN model could be applied to other countries and regions with super-resolution to enhance the generalisability and precision for further trustworthy estimation.

## VII Conclusion

This paper proposed a multi-head CNN model based on pre-trained ResNet50 to estimate demographic composition in Japan as one of the most aged societies from the dataset of satellite images Landsat-8/OLI as daytime and Suomi NPP/VIIRS-DNS as night-time in the year 2015 and 2020. The trained model showed at least 0.8914 in \(R_{(g)}^{2}\) for all the demographic composition groups \(g\) consisting of people below 14 years old, between 15 and 64 years old, and over 64 years old in the testing dataset. The trained model estimated the demographic composition in 2022 as a fulfillment of a dataset in one of the non-census years of Japan. Since the proposed CNN model requires inputs only from publicly available satellite images at any location on the earth, it holds generalisability to other countries and regions in terms of data availability. Although further studies are expected for evaluating the performance in different geographical locations and enhancing the coarse spatial resolution in the inputs, this paper shall contribute to policymaking for solving problems amongst population aging without conducting a rigorous population survey such as a census.

## Acknowledgment

The author would like to thank Dr. Elisabetta Pietrostefani as his supervisor in his Master's program at the London School of Economics and Political Science for her advice provision

Fig. 9: Estimated rate of aged population in Japan as of 2022. JGD2011 was applied for the coordinate reference system. The map is based on OpenStreetMap.
to this paper as the author's dissertation, including the data collection and visualisation methods.
