TR2RM: an urban road network generation model based on multisource big data

ABSTRACT

Road networks are an important part of transportation infrastructure through which people experience a city. The existing methods of vector map data generation mainly depend on a single data source, e.g. images, trajectories, or existing raster maps, which are limited by information fragmentation due to incomplete data. This study proposes an urban road network extraction framework named trajectory and remote-sensing image to RoadMap (TR2RM) based on deep learning technology by combining high-resolution remote sensing images with big trajectory data; this framework is composed of three components. The first component focuses on feature map generation by fusing remote sensing images with trajectories. The second component is composed of a novel neural network architecture denoted as AD-LinkNet, which is used to identify roads from the fused dataset of the first component. The last component is a postprocessing step that aims to generate the vector map accurately. Taking Rome, Beijing, and Wuhan as examples, we conduct extensive experiments to verify the effectiveness of the TR2RM. The results showed that the correctness of both the topology and geometry of the generated road network based on the TR2RM in Rome, Beijing, and Wuhan was 83.86% and 88.27%, 74.72% and 80.36%, and 73.83% and 77.7%, respectively.

KEYWORDS:

• Road network generation
• semantic segmentation
• remote sensing
• deep learning
• big trajectories

1. Introduction

Road network data play a crucial role in many fields, such as urban planning, intelligent transportation, and navigation (Li et al. 2021a). Approaches to generating road network data from a single data source have been widely discussed in recent years. One of the most common approaches for obtaining large-scale road network data is to use remote sensing images as the data source. The technology for identifying roads from remote sensing images involves two stages: ‘obtaining road information using classic machine learning methods’ and ‘identifying roads based on artificial intelligence (AI) technology’. In the first stage, road elements at the pixel level were recognized from images based on classic extraction algorithms, such as digital morphology (Lin et al. 2008), tensor voting (Medioni et al. 2008), conditional random fields (CRFs) (Wegner, Montoya-Zegarra, and Schindler 2013), or templates (Miao et al. 2014). These methods face the issues of a long processing cycle and poor generalization capability, which promote the development of the second stage. In the second stage, AI techniques are used to improve the road identification performance. Approaches in this stage can be classified into two types: constructing road network data based on generative artificial intelligence (AI) techniques (Juhász et al. 2023) and extracting road elements using deep neural networks (Sun et al. 2019a). The former pays much attention to generating nonexistent but realistic images based on models such as CycleGAN (Zhu et al. 2017), StyleGAN (Karras, Laine, and Aila 2019), or ChatGPT (Zhang et al. 2023). In comparison, the second type of AI model focuses on improving the road information extraction performance from the original images; examples include the deep convolutional encoder-decoder network (DCED) (Panboonyuen et al. 2017), D-LinkNet (Zhou, Zhang, and Wu 2018), HsgNet (Xie et al. 2019), and D-CrossLinkNet (Huang et al. 2020). Although these efforts have indeed made considerable progress in recent years, the quality of the image itself, caused by overexposure or sheltering by trees or buildings, reduces the accuracy of information extraction and leads to information fragmentation (Liu et al. 2018). As shown in Figure 1a, roads highlighted by yellow lines could be easily missed during identification because of sheltering.

Figure 1. A comparison between remote sensing images and trajectories: (a) road information captured by images and (b) road information recorded by trajectories.

Trajectories collected by thousands of volunteers worldwide (e.g. tracking data shared in the OpenStreetMap (OSM) platform) or taxis and car-hailing companies (e.g. DiDi) are another option for constructing road network data (Karimi and Kasemsuppakorn, 2013; Ruan et al. 2018; Yang et al. 2022a; Zhang et al. 2020). Compared with remote sensing images, trajectories have better topological connectivity, which is beneficial for generating vector road network data at different levels, such as the road centerline level (Deng et al. 2019; Gao et al. 2021) and the traffic lane level (Mahmoud et al. 2021). Three kinds of methods have emerged for generating road network data based on trajectories: incremental insertion, clustering, and machine vision-based methods (Ahmed et al. 2015; Yang et al. 2020). The first obtains road networks by incrementally inserting tracking points into an initially empty map based on map-matching ideas (Gao et al. 2021). The second type used a set of tracking points as input and clustered them in different ways to acquire road segments or intersections and construct road network data, e.g. Mahmoud et al. (2021), Deng et al. (2019), and Yang et al. (2018a). The last category extracts road networks by transforming trajectories into discretized images and recognizing road segments via computer vision algorithms, such as road information generation based on the DeepMG (Ruan et al. 2020) and road intersection identification via the Mask-RCNN model (Yang et al. 2020). Regardless of the specific technique used for road network generation, trajectories, whether collected by taxis, online car-hailing, or other people, still face the issues of unbalanced distributions and noise (Abdollahi et al. 2020; Sun et al. 2019a), which can lead to missed or misidentified roads.

By analyzing studies of road network data generation using remote sensing images or Global Navigation Satellite System (GNSS) tracking data, issues caused by the data itself always exist. Relying only on extraction techniques is not enough to eliminate their shortcomings entirely. Therefore, how to obtain road information through fostering the strength of multisource data and circumventing their weaknesses has been the focus of research. For instance, Sun et al. (2019b) proposed extracting road information from aerial images first and then applying crowdsourced trajectories to repaint the identified road segments. Zhang et al. (2020) rasterized the GPS trajectory of a floating car into a raster map and used the processed raster map to label the satellite image to obtain the road extraction sample set to reduce the manual labeling cost. Li et al. (2021b) demonstrated the use of a multisource fusion network (MTMSAF) to simultaneously extract roads and detect intersections by combining remote sensing images with GPS tracking data. Qin et al. (2022) proposed an incremental road network updating method that combines GPS trajectory data and UAV remote sensing images. This method uses the hidden Markov model (HMM) to identify problematic road sections and uses deep learning technology to automatically extract road sections from remote sensing images to update problematic road sections. However, previous works failed to address the features that represent road information from trajectories; instead, they simply superimposed them on remote sensing images for road recognition. This may increase the probability of misidentification due to data inconsistencies and the problems of road object fractures and poor connectivity not being well improved.

This study proposes extracting road networks by fusing remote sensing images with big trajectories based on the novel framework TR2RM, which consists of three components. The first component focuses on merging remote sensing images and trajectories into feature maps. In the second component, we design a novel model, AD-LinkNet, to identify road segments using feature maps as the input. The last component is a postprocessing algorithm used to correct the errors of the generated road networks of the second component. We verify the proposed TR2RM by using three kinds of datasets collected in Rome, Beijing, and Wuhan in various scenarios. The main contributions of this study are as follows.

1. We construct a new framework, TR2RM, to automatically extract road networks by combining remote sensing images with trajectories.

2. A novel model, AD-LinkNet, is designed based on the original version of LinkNet to increase the suitability of road network generation from the fused dataset.

3. We present a postprocessing method to solve the problems of small, isolated areas and incorrect topological relations in the generated results.

2. Methodology

The overall structure of TR2RM is shown in Figure 2. The first component of TR2RM is used to fuse tracking data with the corresponding remote sensing images. Then, we used the fused data to train the newly designed AD-LinkNet model to further identify roads. The third component is applied to recognize isolated small areas or road breakages and correct them to obtain a complete vector road network.

Figure 2. The structure of TR2RM.

2.1. Feature extraction and data fusion

In this study, we converted the coordinate systems of remote sensing images, trajectories, and ground truths into WGS84 data to avoid position deviation. Here, ground truths were obtained from the OSM platform based on areas covered by remote sensing images and modified manually against high-resolution remote sensing images from Google Earth to maintain reliability. The crowdsourced big trajectory data are composed of a set of tracking points, and each of them records spatiotemporal information about the moving object, including the longitude, latitude, elevation, time stamps, and speed. To facilitate trajectory data computation by using a deep learning model and eliminate errors between remote sensing images and trajectory data, we rasterized the data based on the pixel size of the corresponding remote sensing image. Specifically, we obtained corresponding remote sensing images from Google Earth based on the area covered by trajectories and the collection time. In Figure 3, four vertices record the coverage scope of trajectories in the experimental area. The parameters Lat_max and Lat_min and Lon_max and Lon_min indicate the maximum and minimum values of latitude and longitude, respectively. The total number of pixels in both the horizontal direction and vertical direction in the corresponding remote sensing image is denoted as I and J, respectively. Each pixel can be identified based on its column number in the grid image. Based on the pixel size of the corresponding remote sensing image, we can generate a new blank image. The intervals for both the longitude (denoted as Interval_lon) and latitude (denoted as Interval_lat) of each pixel in this blank image can be computed based on Equation (1). For a tracking point p_i(lat_i, lon_i), we can calculate which pixel grid in the blank image it is in based on Equation (2). In Equation (2), parameters x and y represent the column number of the pixel in the grid image. The function ‘Math. round’ is used to round a numerical value.

Figure 3. The conversion of trajectories to raster based on the pixel size of the remote sensing images.

(1) $\mathrm{Interva}{\mathrm{l}}_{\mathrm{lon}}=\frac{(\mathrm{Lo}{n}_{max}-\mathrm{Lo}{n}_{min})}{I},\mathrm{Interva}{\mathrm{l}}_{\mathrm{lat}}=\frac{(\mathrm{La}{t}_{max}-\mathrm{La}{t}_{min})}{J}$(1) (2) $x=\frac{\mathrm{Math}.\mathrm{round}(p_{\mathrm{lon}}-\mathrm{Lo}{n}_{min})}{\mathrm{Interva}{l}_{\mathrm{lon}}},y=\frac{\mathrm{Math}.\mathrm{round}(p_{\mathrm{lat}}-\mathrm{La}{t}_{min})}{\mathrm{Interva}{l}_{\mathrm{lat}}}$(2)

The gray value of each pixel in the blank image is determined by the features of the tracking points falling into the pixel. According to the existing work, the higher the tracking point density is, the more likely it is to be a road (Yang et al. 2020). Moreover, Ruan et al. (2020) indicated that the speeds of moving objects on roads or nonroads (e.g. parking lots) were different. Thus, we obtain the gray value of each pixel based on the point density and the speed at which the tracking point falls into it. Moreover, since each tracking point represents the geographic location of the moving object at time stamp t_i, the trajectory composed of these tracking points also represents the approximate geometric structure of the road. Therefore, ‘Distance’ is the last feature to participate in gray value computation. The specific steps for obtaining the gray value of each converted pixel of trajectories are illustrated as follows.

1. Point density. Let the number of trajectory points falling in a grid unit (also called a grid cell or pixel cell) T be n, and we use the threshold value N to confirm the gray value. That is, for N = 1, the gray value of a pixel cell is set to 0 when n is smaller than N. The gray value of a pixel cell is set to 128 if n is equal to N; otherwise, the gray value is set to 255.

2. Speed. To reduce the influence of outliers, we calculate the standard deviation of the speed of trajectories in a pixel to confirm its gray value. The velocity of each tracking point t is assumed to be denoted as v_t. Equation (3) shows the computation of the average speed v_T of the tracking points. The parameter n represents the total number of tracking points in T. The standard deviation σ_T of the current grid cell can be computed based on Equation (4). (3) $v_T={\displaystyle \frac{\sum _{i=1}^n{v}_i}{n}}$(3) (4) ${\sigma }_T=\sqrt{{\sigma }^2}=\sqrt{\frac{\sum _{i=1}^N{(v_i-v_T)}^2}{N}}$(4) where N is the number of grid cells surrounding the current grid cell T whose speed is not equal to 0 and $v_i$ is the velocity of grid cell $i$. The current grid cell is considered to be isolated if the speeds of its eight neighbors are all equal to 0. The gray value of the isolated grid cell is set to 0. The gray values of the other grid cells are confirmed based on a trisection strategy. That is, the gray value of a grid cell is set to 0, 128, or 255 when σ_T> 0.8, 0.5<σ_T< 0.8, or 0<σ_T< 0.5, respectively.

3. Distance. In this study, we use the distance between two adjacent tracking points to estimate the gray value of the current grid cell regardless of whether the points are in the same grid cell to reduce the influence of sparse sampling. The gray value of a grid cell is set to 255 if both adjacent tracking points fall in it. One tracking point is set to 128 if it falls inside the grid cell and its adjacent point falls within the neighbors of the grid cell. In other cases, it is set to 0.

The gray value of each pixel in the conversion of trajectories to raster data includes 0, 128, and 255, representing nonroad pixels, potential road pixels, and road pixels, respectively. The reason why we converted the point density, speed, and distance of trajectories into three levels is that two levels increase the probability of false recognition due to the uncertainty of trajectories. Moreover, four levels increase the computational complexity of the deep learning model from an exponential multiplication of 3–4. Additionally, to reduce the uncertainty caused by noise and outliers, we smoothed the raster map converted by trajectories by counting pixel types in an 8-neighborhood window. Specifically, for a pixel, if its original type is a ‘road pixel’ or ‘potential road pixel’, then it will be referred to as a ‘nonroad pixel’ when the total number of pixels without road information in its 8-neighborhood field is greater than H. In contrast, if the original type of a pixel is a ‘nonroad pixel’, then it will be referred to as a ‘potential road pixel’ when the total number of pixels with road information in its 8-neighborhood field is greater than H. Based on this smoothing strategy, we can obtain a three-dimensional vector in each grid cell and express it as $X_{\mathrm{traj}}\in {\mathbb{R}}^{I\times J\times 3}$. Moreover, the remote sensing images are expressed as $X_{\mathrm{rs}}\in {\mathbb{R}}^{I\times J\times 3}$. As shown in Figure 1 in the Appendix, we connect three types of feature maps generated from trajectories with the corresponding remote sensing images based on the channel latitude to ensure that the size of the input dataset of the subsequent deep learning model is not changed (see Eq. 5). (5) $X_{\mathrm{merge}}\in {\mathbb{R}}^{I\times J\times 6}=\mathrm{concatenate}(X_{\mathrm{traj}}\in {\mathbb{R}}^{I\times J\times 3},X_{\mathrm{rs}}\in {\mathbb{R}}^{I\times J\times 3})$(5) where$X_{\mathrm{merge}}\in {\mathbb{R}}^{I\times J\times 6}$ is a new feature map obtained by integrating the raster grids converted by trajectories with remote sensing images.

2.2. Road detection based on a newly designed AD-LinkNet model

Figure 4 shows the newly designed AD-LinkNet model obtained by introducing a new AD-block to replace the middle layer cascaded by the multilayer dilated convolution of D-LinkNet. We use a 1 × 1 convolution to reduce the dimension of feature maps from a 6-dimensional channel to a 3-dimensional channel in the beginning part of AD-LinkNet, highlighted by the red dotted box in Figure 4, to ensure that the features of two different datasets can be fully integrated.

Figure 4. The structure of the AD-LinkNet.

The newly designed AD-Block in the AD-LinkNet model is used to overcome two limitations of D-LinkNet. First, the dilated convolution structure of the D-LinkNet model causes a serious gridding phenomenon during training and prevents many pixels from participating in the convolution operation (Wang et al. 2018). Second, D-LinkNet obtains only local features and ignores the rich context information of images, even though it increases the receptive field of the convolution kernel. We embed the newly designed AD-Block in the four Res-blocks to expand the diversity of the receptive fields in each stage and further increase the semantic accuracy at the pixel level. Moreover, the AD-Block can fuse the features of the decoder component with the features of the following component to avoid bottom feature and high-level feature loss and improve the quality of stored images of detailed road features by preserving global information. Figure 5 shows that the newly designed AD-Block is composed of three criss-cross attention (CCA) modules and a hybrid dilated convolution (HDC) module with dilation rates of 1, 2, and 5.

Figure 5. The structure of the AD-Block.

The CCA module (Huang et al. 2019) embedded in the AD-Block is used to enhance the ability of global information learning by AD-LinkNet without increasing the amount of calculation and memory usage by considering its advantages in context information collection in both the horizontal and vertical directions. The steps of the CCA module for collecting the pixel-level context information from the horizontal and vertical directions of the feature map are as follows. Taking a downsampled feature map $X\in {\mathbb{R}}^{C\times H\times W}$ as the input data, two 1 × 1 convolution kernels are used to process the data. After the convolution calculation, we obtain two new feature maps $\{Q,K\}\in {\mathbb{R}}^{C^{\mathrm{\prime }}\times H\times W}$where $C^{\mathrm{\prime }}$ indicates the number of channels. Since Q and K are essentially the same feature map, we can obtain a set ${\mathrm{\Omega }}_u\in {\mathbb{R}}^{C^{\mathrm{\prime }}\times (H+W-1)}$ by gathering a certain location u in K and its feature points in the same row and column. The parameter ${\mathrm{\Omega }}_{i,u}\in {\mathbb{R}}^{C^{\mathrm{\prime }}}$ represents the i-th element in ${\mathrm{\Omega }}_u$. To obtain the correlation between points in the set ${\mathrm{\Omega }}_u\in {\mathbb{R}}^{C^{\mathrm{\prime }}\times (H+W-1)}$, the similarity computation shown in Equation (6) is used. (6) $\mathrm{Scor}{\mathrm{e}}_{i,u}=Q_u{\mathrm{\Omega }}_{i,u}^T$(6) where $\mathrm{Scor}{\mathrm{e}}_{i,u}$ is the correlation score of point u and its feature points in the same row and column. Then, the channel compression method is applied to obtain the required attention map $A\in {\mathbb{R}}^{C\times H\times W}$. That is, we convolve $X\in {\mathbb{R}}^{C\times H\times W}$with a 1 × 1 convolution kernel to obtain a new feature map $V\in {\mathbb{R}}^{C\times H\times W}$. Based on the location u in the set of V and its feature points in the same column, we can obtain a set ${\mathrm{\Phi }}_u\in {\mathbb{R}}^{C\times (H+W-1)}$. We use an aggregation operation to further obtain the feature map weighted according to positional relationships, as shown in Equation (7). (7) $X_u^{\mathrm{\prime }}=\sum A_{i,u}{\mathrm{\Phi }}_{i,u}+X_u$(7) where $X_u^{\mathrm{\prime }}$ is the feature vector of the output feature map $X^{\mathrm{\prime }}\in {\mathbb{R}}^{C\times H\times W}$ at position $u$ and $A_{i,u}$ is the scalar value of the attention map $A$ at channel $i$ and position $u$. To further obtain the contextual information from all the pixels, we connected several CCA modules in series and embedded them in the AD-Block. The number of CCA modules embedded in the AD-Block is denoted as n.

We add the HDC module to the AD-Block to improve the receptive field and avoid the occurrence of the grid phenomenon (Wang et al. 2018). The dilation rate of the HDC module in the AD-Block is set to 1, 2, or 5 based on the standardized extended convolution block principle of HDC proposed in Wang et al. (2018) (see Figure 6).

Figure 6. The effect diagram of HDC.

In the HDC module, the equivalent size of the ordinary convolution kernel after the expansion operation can be computed based on Equation (8). (8) $k^{\mathrm{\prime }}=d\times (k-1)+1$(8) where k is the size of the input convolution kernel, d represents the expansion coefficient and can be artificially set as a hyperparameter, and $k^{\mathrm{\prime }}$is the equivalent size of the ordinary convolution kernel after the expansion operation. According to Equation (8), the receptive field of the HDC module is equivalent to the receptive field of standard convolution kernel sizes of 3, 5, and 11. In Figure 6, the receptive field of the HDC module covers the whole area, eliminating the gridding effect and obtaining more complete contextual information. After the expansion operation, the size of the feature map outputted from the HDC module can be calculated according to Equation (9). (9) $o=\lfloor {\displaystyle \frac{i+2p-k^{\mathrm{\prime }}}{s}}\rfloor +1$(9) where $i$ is the size of the input feature map, $o$ is the size of the output feature map, $p$ is the size of the pooling operation, $s$ represents the step size, and $\lfloor \dots \rfloor$ indicates a rounding-down operation. It should be noted that the convolution kernel size of the expansion convolution in this study is set to 3 × 3 to ensure that the sizes of the input feature map and output feature map are equal when the values of s and p are equal to 1 and d, respectively. In addition, by analyzing the similarities between ResNet with skip connections and the dilated convolution block, we add the skip connection to the HDC module to enhance and supplement the original features. That is, we connect the output feature map of the HDC module with the original input feature map to strengthen it. Then, we apply a 1 × 1 convolution to fuse their features and further reduce the computational overhead in the subsequent operations.

2.3. Postprocessing

At present, CNN-related methods mainly focus on object identification and cannot solve the issues of road breakage and small isolated areas caused by road feature loss or poor-quality images during road information extraction (Li et al. 2021c). In this study, we propose a postprocessing method to address the problems of fracture of detected road objects, poor connectivity, and small isolated areas generated after semantic segmentation based on AD-LinkNet. The proposed postprocessing procedure includes four steps: (1) isolated small area removal, (2) breakpoint detection, (3) breakpoint reconnection, and (4) road network refinement. The execution details for each step are described below.

1. Isolated area removal. We removed isolated areas from the outputted images of AD-LinkNet based on the principle of road connectivity. That is, a road segment should connect with one or more road segments in the real world. The isolated small areas are usually incorrectly detected and need to be removed. Therefore, we extract the contour information of each independent road segment in the image and set it as a vector set $V$, where each element retains a set vector v of points constructed by continuous edge points; see Equation (10). Each point set is an independent connected domain outline. We calculate the area of each connected domain based on the obtained contour information (see Eq. 11). Then, the connected domains are regarded as the isolated area and removed if the area is smaller than the threshold $\delta$. (10) $V=\{v_1,v_2,\dots ,v_n\},v_i=\{p_1,p_2,\dots ,p_m\}$(10) (11) $\mathrm{are}{\mathrm{a}}_i=\mathrm{contourArea}\left(v_i\right)$(11)

2. Breakpoint detection. The output dataset of AD-LinkNet is composed of binary images that contain only two kinds of objects, namely, road points and nonroad points. Based on the observations that the gray value of a point will extend around itself if it is not a breakpoint; otherwise, its gray value will change substantially compared with its surrounding points. According to this observation, we detect the breakpoint from the output result of AD-LinkNet by the following steps. First, we use the Zhang – Suen skeleton extraction algorithm (Zhang and Suen 1984) to extract skeleton lines from road areas for subsequent topology reconstruction and breakpoint connection operations. As shown in Figure 7, each road area is regarded as an independent domain and denoted A, B, C, D, and E (see Figure 7a). By using the Zhang – Suen skeleton algorithm, we obtain the refined road skeleton line, as shown in Figure 7b. Then, the Shi-Tomasi Corner detector (Shi 1994) is applied to identify potential breakpoints, which can also be called corner points.

Figure 7. Breakpoint detection (a) Output results of AD-LinkNet (b) road skeleton line (c) corner point and (d) breakpoint.

Given a point $(x,y)$, its pixel value is denoted as $I(x,y)$. $w(x,y)$ represents the corresponding window function, and its value is set to 1. Let this point move along the X- and Y-axes with tiny increments of u and v, respectively. The change in pixel gray value obtained by moving the window in all directions can be calculated based on Equation (12). (12) $E(u,v)=\sum _{(x,y)}w(x,y)\times {\left[I\right(x+u,y+v)-I(x,y\left)\right]}^2$(12) Based on the computation shown in Equation (12), the point (x, y) is regarded as a corner point if the value of E(u,v) is very large. In addition, we use Taylor expansion to simplify the computation of Equation (12); see Equations (13) and (14). (13) $E(u,v)\approx [u,v]M\left[\begin{array}{c}u\\ v\end{array}\right]$(13) (14) $M=\sum _{(x,y)}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\to R^{-1}\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]R$(14) where $M$ is a real symmetric matrix, and its eigenvalues ${\lambda }_1$ and ${\lambda }_2$ can be extracted after diagonalization. Based on these computations, a corner detection function RF is obtained to find the corner points; see Equation (15). (15) $\mathrm{RF}=min({\lambda }_1,{\lambda }_2)$(15) As shown in Figure 7c, the corner points of each skeleton line can be detected based on the above computation. Since not all corner points belong to breakpoints, we use an eight-neighborhood method to further confirm breakpoints. A corner point p is considered an isolated point and is abandoned when there is no road point in its neighborhood $\delta \left(p\right)$. A road point is confirmed if there are only two road points in $\delta \left(p\right)$. For other situations, road point p is considered a bifurcation point.

1. Breakpoint reconnection. To ensure the completeness of the road network in terms of topology, we need to connect the detected breakpoints. The reconnection rules of breakpoints in this study include two points. First, the closer two nodes are, the more similar the two nodes are, and the more likely they are to belong to the same road. Second, the change in the direction trend between the current road point and its adjacent road points is not very large if they are on the same road. Based on these two rules, we design Algorithm 1 for breakpoint reconnections.

Algorithm 1: ConnectBreakPoints

Table

Input: Road networks detected from RS + traj $D_i\in R^{I\times J\times 3}$

Output: Enhanced road network $D_i^{\mathrm{\prime }}\in R^{I\times J\times 3}$

1: data=[length,width] // Store the road type represented by the pixel

2: flag=[length,width] // Flagbit that determines whether the current pixel has been accessed

3: radis // Radius of buffer

4: length // Maximum length

5: TypeImg($D_i$,data)

6: for int x=0 to length do // i rows

7: for int y=0 to width do // j columns

8: if data[x][y] is breakpoint then

9: neighbor(data[x,y])

10: for int i=0 to 8 do

11: if data[i]! =0 then

12: circle(data[x][y], data[i],radis) // Extend the radius in the direction of point i

13: if another breakpoint in circle, then

14: t_x=breakpoint_lat

15: t_y=breakpoint_lng

16: connect(data[x][y],t_x,t_y)

17: end if

18: if no breakpoint in circle, then

19: Traverse each pixel to obtain the nearest road point t(x,y) from data[x][y]

20: if distance(x,y,t_x,t_y)<length then

21: connect(data[x][y], t_x,t_y)

22: end for

23: end for

24: return $D^{\mathrm{\prime }}$// result

In Algorithm 1, the direction D of the road from breakpoint point p to road point $p^{\mathrm{\prime }}$ is calculated. By moving the grid unit r forward in the direction D to search for breakpoints within the buffer of r, the breakpoint q and the breakpoint p will be connected if their Euclidean distance is the lowest among the breakpoints in the buffer. The breakpoint p will be connected to a road point with the lowest Euclidean distance if there is no other breakpoint in the buffer from the start to the end of the search. In addition, breakpoint p is considered a termination point with no connection to other road points if there are only nonroad points in the buffer from the start to the end of the search. The parameter r is a hyperparameter and is discussed in the experiment.

1. Road segment restoration. We use the region growing algorithm (Adams and Bischof 1994) to restore the width of the detected road. The size of the growth area is calculated based on the average width of two connected road sections, as shown in Equation (14). (14) $m_{i,j}={\displaystyle \frac{1}{2}\left({\displaystyle \frac{V_i}{T_i}+{\displaystyle \frac{V_j}{T_j}}}\right)}$(14) where m_i,j indicates the size of the growth area, i and j represent the nodes of two road sections that are going to be connected, V_i and V_j are the numbers of pixels in the connected domain where i and j are located before the thinning operation, and T_i and T_j are the numbers of pixels in the connected domain where i and j are located before after the thinning operation. A road network map was constructed by superimposing the image generated with the pixel growth algorithm, and the results were isolated area removal.

3. Experiment

3.1. Datasets and preprocessing

We used three datasets collected in Rome, Beijing, and Wuhan to verify the proposed method. In Rome, trajectories were collected by more than 300 taxis in one month in 2014 and were provided by Ahmed et al. (2015). In Wuhan, trajectories were gathered by 83 pedestrians within a month in 2016 and given by Yang et al. (2020). Based on the collection times and coverage areas of trajectories in these two areas, we obtained corresponding remote-sensing images from Google Earth. Trajectories and images in Beijing were obtained from Sun et al. (2019b); these data were collected from 28,000 taxis within one week and obtained from the Gaode Map in 2018. Table 1 shows the details about the datasets used in this study.

Table 1. Dataset details.

Items	Rome	Beijing	Wuhan
Area of the experimental region (km^2)	399.96	44.17	6.48
Data Source (RS GPS)	(Goole Earth)	Sun et al. (2019b)	(Goole Earth)
Data Collection Time	2014	2018	2016
Total number of tracking points	2482340	1630796	696318
Resolution of remote sensing images (cm/pixel)	44.41	50	51.43
Image Resolution	40960 × 30464	1024 × 1024	4864 × 4352
Data Form (RS GPS)	(TIF CSV)	(TIF PKL)	(TIF CSV)
Training Set	1728	278	106
Validation Set	432	70	26
Testing Set	432	70	26
Total number of vertices of the ground truth	3051	3166	513
Total number of edges of the ground truth	3050	3165	512

To increase the certainty of the trajectories, we first removed outliers from the original datasets using the method proposed by Yang et al. (2018b). Then, we filtered the tracking points again based on their speed. That is, the tracking point p_i is removed if its velocity from p_i−1 to p_i is greater than the maximum velocity threshold v_max. Here, the maximum speeds v_max of driving and walking are set to 20 and 2.83`m/s, respectively, according to previous work (Yang et al. 2020). The trajectories covered in Beijing, Rome, and Wuhan were converted into grids based on the resolution of the corresponding images and used to construct a feature map based on the method illustrated in Section 2.1. According to existing work (Sun et al. 2019b), images of Beijing were clipped into a series of sub-images based on a size of 1024 × 1024. To maintain consistency, datasets, including the corresponding feature maps in Beijing and images and feature maps in Rome, were also clipped to the same size. In Wuhan, since the experimental area was far smaller than that in Rome and Beijing, we set the input size of the images and feature maps in TR2RM to 512 × 512 according to the conventions of image clipping (Sun et al. 2019a). The ground truth was used to quantitatively evaluate the generated road network (Table 1).

3.2. Experimental setting and evaluation indicators

Codes of the proposed TR2RM in this study were written using Python 3.6, and PyCharm2022.3.7 was used as an integrated development environment. The newly designed AD-LinkNet in TR2RM was developed based on PyTorch 1.4.0 and trained on a server equipped with two NVIDIA GeForce RTX 2080Ti graphics cards. The Room dataset contains 2160 images, the Beijing dataset contains 348 images, and the Wuhan dataset contains 132 images. The training, validation, and test sets were divided at a 6:2:2 ratio for these datasets in this study (Huang et al. 2020). The specific samples in each category during model practice were automatically assigned from the dataset based on their percentage. The parameters, including the initial learning rate, momentum, and weight attenuation coefficient, were set to 0.001, 0.95, and 0.0004, respectively. Moreover, the loss function of the proposed AD-LinkNet applies a weighted combination of Binary Cross Entropy (BCE) and Dice loss (see Equation 15). (15) $\mathrm{L}=\lambda L_{\mathrm{BCE}}+(1-\lambda )L_{\mathrm{Dice}}$(15) where $\lambda$ is a hyperparameter used to control the proportion of two different loss functions. In the experiment, this value was set to 0.5. The batch size for model training and testing was 4. The number of training rounds on each dataset is 200, thus achieving the best performance of the model on two GPUs and SyncBN. Additionally, we applied random scaling and flipping during model training and testing. The specific steps include horizontal flip, vertical flip, rotation, random operation selection, and random noise operation. Specifically, images in the original dataset were horizontally flipped with a probability of 50% first. Then, images obtained in the previous step and from the original dataset were vertically flipped with a probability of 50% and processed with a clockwise rotation of 180 degrees with a probability of 50%. The next step involved translating and scaling the data to a probability of 20%. The last step involved adding a certain amount of noise to all the images with a probability of 20%. By using these steps, we enhanced the learning ability of the neural network model, which avoids the risk of overfitting and makes the model more robust and generalizable.

The value of parameter r during postprocessing was computed based on the width of the road, as shown in Equation (16). (16) $r={\displaystyle \frac{R_{\mathrm{width}}}{I}}$(16) where R_width indicates the width of the road and parameter I denotes the area represented by a single grid cell. According to road construction standards in China, the widths of urban trunk roads and roads or streets in residential areas range from 30 to 40 m and from 7 to 9 m, respectively. In this study, most of the roads in the experimental area in Rome and Beijing were urban trunk roads. In Wuhan, the type of roads included roads or streets in a residential area. Therefore, the values of r in Rome, Beijing, and Wuhan were set to 80, 80, and 20, respectively.

The performance of the proposed AD-LinkNet in TR2RM was estimated based on the following indicators: accuracy, precision, recall, overall and individual mIoU, and F1-score. The methods used to compute these indicators were the same as those used in existing studies (Yang et al. 2022a). Moreover, the correctness of the geometry and topology is used to estimate the performance of the TR2RM. The calculation steps for geometric correctness and topological correctness are illustrated in the Appendix based on existing studies (Kasemsuppakorn and Karimi 2013). According to the urban road construction standard of China, the buffer distance used to match the road with the benchmark is set to 10 m.

3.3. Result evaluation and comparison

3.3.1. Estimation of road generation using a single dataset and fused datasets

To improve this limitation, by using a single dataset to extract road networks, we propose generating road networks by fusing images with big trajectory data based on the designed AD-LinkNet model. To verify the performance of road extraction by using the fused datasets, we used three kinds of datasets to estimate the generation results of the road network. These three datasets included a trajectory dataset, remote sensing images, and a combined dataset of trajectories and images from the cities of Rome, Beijing, and Wuhan. These three different datasets are denoted as Traj, RS, and RS + Traj, respectively. Table 2 shows the estimation results of road generation by using these three datasets. Here, the best result of each evaluating indicator is displayed in bold.

Table 2. Road network generation results from different data sources.

Area	Dataset	Accuracy	Precision	Recall	mIoU	F1
Rome	Traj	87.90	50.86	76.72	44.06	61.17
	RS	93.77	69.15	90.06	64.24	78.23
	RS + Traj	94.39	71.42	91.33	66.89	80.16
	Traj	91.49	61.89	77.67	52.54	68.89
Beijing	RS	93.22	67.43	85.21	60.37	75.29
	RS + Traj	93.91	70.34	86.07	63.15	77.41
	Traj	91.42	59.76	82.35	52.98	69.26
Wuhan	RS	79.20	22.83	48.55	18.38	31.06
	RS + Traj	92.83	66.56	82.83	58.49	73.81

Based on the statistics shown in Table 2, we can find that the performance of road identification in the cities of Rome, Beijing, and Wuhan using RS + Traj is better than that of using a single dataset such as Traj or RS. Specifically, the average values of accuracy, precision, recall, mIoU, and F1 score for road network generation by using the RS + Traj dataset in Rome, Beijing, and Wuhan were approximately 93.71%, 69.44%, 86.74%, 62.84%, and 77.13%, respectively. In the city of Rome, the mIoU and F1 score of the generated road networks based on the fused dataset increased by 2.65% and 1.93%, respectively, compared with those of the single-dataset RS. Similarly, the mIoU and F1 score of road identification also increase by nearly 3% and 2%, respectively, when using the fused dataset in Beijing compared with when using only the RS dataset. Moreover, the estimation results shown in Table 2 also indicate that the road extraction performance based on the fused dataset is obviously better than that based only on the Traj dataset in the cities of Rome and Beijing. In Wuhan, the mIoU and F1 score of road identification based on the RS + Traj dataset increase by 5.51% and 4.55%, respectively, compared with those of the test based on the Traj dataset. Interestingly, we also found that the performance of the RS model was superior in all the cases except for Wuhan, where the accuracy, precision, recall, mIoU, and F1 score were significantly lower. Moreover, the statistics illustrated that the road extraction performance in Wuhan using the Traj dataset was better than that using the RS dataset in both Beijing and Rome. This is a result of the differences in data quality and completeness in these three regions. For the RS dataset, manual inspection revealed that the issue of road sheltering by trees in images in Wuhan was more serious than that in images from the other two regions (see Figure 8a). This leads to lower values of all elevation indicators of road identification using RS in Wuhan, e.g. accuracy and precision, than those of Beijing and Rome. In addition, the imbalance and incomplete coverage of trajectories in Beijing and Rome are more serious than those in Wuhan, which results in roads with fewer tracking points or even no misclassified or missed tracking points.

Figure 8. Visualization of road network generation results by using different datasets from the cities of Rome Beijing and Wuhan; (a) RS dataset (b) Traj dataset (c) ground truth (d) road network generation using the RS dataset (d) road network generation using the Traj dataset and (f) road network generation using the RS + Traj dataset.

Figure 8 shows several typical identification examples from different datasets in Beijing, Rome, and Wuhan. As shown in Figure 8d, some road segments in Rome are misidentified using RS due to sheltering. This situation is also obvious in Beijing and Wuhan. In Wuhan, trees and buildings sheltered most of the roads in the experimental area cannot be identified based only on remote sensing images, as shown in the last two rows of the images in Figure 8. Compared with road generation results using RS, although road extraction using the Traj dataset can avoid losing some road segments, there is still the problem of incomplete information, especially in Beijing and Wuhan. In Figure 8e, we can see that some road segments are missed because there are only a small amount of tracking data or even no tracking data. By using the fused datasets, the errors shown in Figures 8d and e are improved; see Figure 8f.

3.3.2. Performance estimation for AD-LinkNet

The proposed AD-LinkNet is compared with 7 other typical models, namely, UNet (Ronneberger, Fischer, and Brox 2015), SegNet (Badrinarayanan, Kendall, and Cipolla 2017), DenseNet (Huang et al. 2017), Deeplabv3+ (Chen et al. 2017a, 2017b), LinkNet (Chaurasia and Culurciello 2017), D-LinkNet (Zhou, Zhang, and Wu 2018), and NL-LinkNet (Wang, Seo, and Jeon 2021), to evaluate its performance on road segmentation based on the Traj + RS datasets in the cities of Rome, Beijing, and Wuhan. It should be noted that we did not add an ablation test because the only difference between AD-LinkNet and D-LinkNet was that the former had a new module, AD-Block. Therefore, the comparison between AD-LinkNet and D-LinkNet is equivalent to an ablation test. Table 3 lists the quantitative results of each model by using the dataset collected in the city of Rome. In Table 3, the proposed AD-LinkNet achieves the best performance on the Rome dataset when compared with other models, in which the mIoU and F1 scores reach 66.89% and 80.16%, respectively. Moreover, the mIoU and F1 score obtained by using AD-LinkNet increased by nearly 2% and 1%, respectively, compared with those obtained by using NL-LinkNet.

Table 3. Comparison of segmentation performance of commonly used models on the Roman dataset.

Method	Accuracy	Precision	Recall	mIoU	F1
Unet	89.88	56.07	85.59	51.23	67.75
SegNet	81.17	60.18	85.48	54.60	70.63
DenseNet	91.13	59.38	90.43	55.87	71.69
Deeplabv3+	91.85	61.72	90.50	57.96	73.39
LinkNet	93.78	69.26	89.74	64.18	78.18
D-LinkNet	93.83	69.10	91.08	64.72	78.58
NL-LinkNet	93.89	69.32	91.08	64.91	78.72
AD-LinkNet	94.39	71.42	91.33	66.89	80.16

We also conducted the same comparison experiment by using the Beijing dataset. According to the results shown in Table 4, the mIoU and F1 score for road network generation based on AD-LinkNet increased by 1.62% and 1.22%, respectively, compared with those of the best performing baseline.

Table 4. Comparison of segmentation performance of commonly used models on the Beijing dataset.

Method	Accuracy	Precision	Recall	mIoU	F1
Unet	92.30	63.27	86.91	57.77	73.23
SegNet	92.86	65.66	86.30	59.46	74.57
DenseNet	93.18	66.81	86.92	60.71	75.55
Deeplabv3+	93.17	67.25	85.13	60.18	75.14
LinkNet	92.86	65.08	88.77	60.13	75.10
D-LinkNet	93.39	67.70	87.07	61.52	76.18
NL-LinkNet	93.38	67.56	87.33	61.53	76.19
AD-LinkNet	93.91	70.34	86.07	63.15	77.41

In contrast to those in the first two datasets, the trajectories contained in the Wuhan dataset were collected by pedestrians and exhibited random and irregular characteristics (Yang et al. 2022b). The statistics shown in Table 5 indicate that both the mIoU and F1 score of road network generation based on AD-LinkNet increase by nearly 2% compared with those of NL-LinkNet. Compared with those of D-LinkNet, the mIoU and F1 values increased by 5.5% and 4.53%, respectively. As shown in Table 5, there is an approximately 8% difference between our study and the UNet model. By analyzing the identification results, we found that the Unet model recognized more nonroad pixels as road points. This will increase the values of false-positives (FPs) and true positives (TPs) while decreasing the values of true negatives (TNs) and false-negatives (FNs). In comparison, AD-LinkNet focused on improving the precision, mIoU, and F1 score, even if some road pixels were missing. Moreover, the results also indicate that the performance of AD-LinkNet on the Wuhan dataset is not as good as that on the first two regions. This may be because the randomness and irregularity of vehicle trajectories are much lower than those of pedestrian trajectories due to the limitations of traffic rules and relatively straight motor roads. Thus, future work should focus on how to detect multiple types of roads in urban areas using fused datasets.

Table 5. Comparison of segmentation performance of commonly used models on the Wuhan dataset.

Method	Accuracy	Precision	Recall	mIoU	F1
Unet	88.84	49.27	90.99	46.98	63.92
SegNet	89.15	50.70	82.09	45.65	62.69
DenseNet	90.32	54.20	89.73	51.04	67.58
Deeplabv3+	91.45	61.30	76.81	51.72	68.18
LinkNet	91.26	58.41	85.54	53.17	69.42
D-LinkNet	91.10	56.90	88.53	52.99	69.28
NL-LinkNet	92.13	62.03	86.35	56.49	72.19
AD-LinkNet	92.83	66.56	82.83	58.49	73.81

3.3.3. Estimation for the postprocessing of TR2RM

Postprocessing is the second component of TR2RM and is used mainly to modify the errors caused by road segmentation via the AD-LinkNet model, such as broken road connection points, blurred road edges, and misidentified areas. To estimate its effectiveness, comparative performance tests were conducted using the fused datasets from Rome, Beijing, and Wuhan. Table 6 shows the results of road generation with and without the postprocessing method.

Table 6. Performance of postprocessing algorithms on various datasets.

Dataset	Data Source	mIoU	F1	Geometrical correctness (%)	Topological correctness (%)
Roman	Original	66.89	80.16	84.18	81.75
	Postprocessing	67.22	80.40	88.27	83.86
Beijing	Original	63.15	77.41	77.20	72.16
	Postprocessing	63.26	77.51	80.36	74.72
Wuhan	Original	58.49	73.81	73.24	71.03
	Postprocessing	58.77	74.04	77.70	73.83

In Table 6, the mIoU and F1 score of road network generation after postprocessing in Rome, Beijing, and Wuhan increased by (0.33%, 0.24%), (0.11%, 0.1%) and (0.28%, 0.23%), respectively. Moreover, the values of geometrical correctness and topological correctness of road extraction after postprocessing in Rome, Beijing, and Wuhan also increased by (4.09%, 2.11%), (3.16%, 2.56%), and (4.46%, 2.8%), respectively. By comparison, we find that the improvement in the correctness of both the geometry and topology after postprocessing is more significant than that in the mIoU and F1 score. This is because postprocessing improves the result by removing several pixels identified as isolated small areas or adding several pixels to modify the breakpoints and blurred road edges; see Figure 9.

Figure 9. Road network extraction results after postprocessing; (a) (b) (c) and (d) show cases of breakpoint reconnection; (d) isolated small area removal.

In Figure 9, the first four columns and the last column mainly indicate the comparisons of breakpoint reconnection and isolated small area removal using the postprocessing method, respectively. The comparison of the results in Figures 9a, c, and d indicates that the postprocessing technique can connect the breakpoints very well, whether there are breakpoints near two parallel roads (see Figure 9a and d) or on intersection roads (see Figure 9c). In addition, the postprocessing method can also connect the breakpoint on the zigzag road, as shown in Figure 9b. The effectiveness of road segment restoration and isolated area removal using the postprocessing method is also verified in Figure 9c and e, respectively. These modifications are highly important for accurately estimating the geometry and topology of generated road vector maps. However, the number of pixels modified during postprocessing is limited, thereby leading to improvements in the mIoU and F1 score that are not obvious.

3.3.4. Visualization of the whole experimental area

This paper provides a panoramic view of road network generation in Rome, Beijing, and Wuhan. The pink mask in the picture shows the road network generated by the framework proposed in this paper, and the white mask overlayed by the pink mask is the ground truth. From the three panoramic views, the road network generated by the framework proposed in this paper is only slightly different in the nuances of roads (Figure 10).

Figure 10. Panoramic view of road network generation in Rome Beijing and Wuhan: (a) road network in Rome (b) road network in Beijing and (c) road network in Wuhan.

4. Conclusion

We propose a novel framework, TR2RM, to automatically extract road networks by integrating high-resolution remote sensing images with big trajectory data. There are three components in TR2RM, namely, feature extraction and fusion, proposed AD-LinkNet model training, and road network postprocessing. Taking Rome, Beijing, and Wuhan as examples, we carried out experiments to estimate the performance of the TR2RM. The results show that the average values of accuracy, precision, recall, mIoU, and F1 score of road network generation using the RS + Traj dataset for Rome, Beijing, and Wuhan were approximately 93.71%, 69.44%, 86.74%, 62.84%, and 77.13%, respectively, which are generally higher than those obtained using the RS and Traj datasets. We compared the proposed AD-LinkNet with 7 other typical models. According to the statistics, compared with those of the best-performing baselines, the mIoU and F1 score of AD-LinkNet on the basis of the fused datasets from Beijing, Rome, and Wuhan increased by nearly (2% and 1%), (2% and 1%), and (2% and 2%), respectively. Additionally, we evaluated the correctness of the generated road networks in three experimental regions from the aspects of geometry and topology before and after postprocessing. The results showed that the correctness of both the geometry and topology in Rome, Beijing, and Wuhan increased by (4.09%, 2.11%), (3.16%, 2.56%), and (4.46%, 2.8%), respectively.

This study has several limitations. First, the experimental area used for validation was still too small because of the limited access to data from both public trajectories and high-resolution remote sensing images. Then, the proposed method extracts road networks at the centerline level. Therefore, one direction of future work should focus on detailed road information detection from fused images and trajectories, such as traffic lanes, turnings, and road infrastructures, and extending the experimental area. Moreover, the implications of merging two datasets with differing temporal information should be explored in future work, such as road change detection, despite our efforts to maintain temporal consistency during data collection. In addition, the postprocessing of TR2RMs should be improved from the perspectives of threshold confirmation and optimization strategy establishment with additional road features in the future.

Acknowledgments

The authors would like to sincerely thank the anonymous reviewers for their constructive comments and valuable suggestions to improve the quality of this article. The authors sincerely thank all the organizations and scholars who provided the data and technical support for this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data and code that support the findings of this study are openly available at https://figshare.com/s/1779f7508e2c5974f792.

Additional information

Funding

This work was jointly supported by the National Natural Science Foundation of China [grant number: 42271449].

Appendix

The calculation method of the correctness both in geometry and topology: the geometric correctness (denoted as δ₁) is defined as the percentage of the generated road network compared to the ground truth and calculated by Equation (A1). The topological correctness (denoted as δ₂) indicates the percentage of correctly connected intersections in the generated road network compared to the number of ground truth which can be computed based on Equation (A2). (A1) ${\delta }_1=n^m/n^g$(A1) (A2) ${\delta }_2=x_m/x_g$(A2) where n^m is the length of road segments who can be matched with the ground truth and n^g is the length of all generated road segments. The parameters x_m and x_g represent the number of correct intersections and the total number of all detected intersections, respectively.

The figure of trajectories and remote sensing images fusing:

Figure A1. Feature extraction and fusion.