A Hybrid Method to Incrementally Extract Road Networks Using Spatio-Temporal Trajectory Data

Abstract

With the rapid development of urban traffic, accurate and up-to-date road maps are in crucial demand for daily human life and urban traffic control. Recently, with the emergence of crowdsourced mapping, a surge in academic attention has been paid to generating road networks from spatio-temporal trajectory data. However, most existing methods do not explore changing road patterns contained in multi-temporal trajectory data and it is still difficult to satisfy the precision and efficiency demands of road information extraction. Hence, in this paper, we propose a hybrid method to incrementally extract urban road networks from spatio-temporal trajectory data. First, raw trajectory data were partitioned into K time slices and were used to initialize K-temporal road networks by a mathematical morphology method. Then, the K-temporal road networks were adjusted according to a gravitation force model so as to amend their geometric inconsistencies. Finally, road networks were geometrically delineated using the k-segment fitting algorithm, and the associated road attributes (e.g., road width and driving rule) were inferred. Several case studies were examined to demonstrate that our method can effectively improve the efficiency and precision of road extraction and can make a significant attempt to mine the incremental change patterns in road networks from spatio-temporal trajectory data to help with road map renewal.

1. Introduction

With the ongoing rapid urbanization, road networks have become increasingly complicated and heavily dense. Generating timely, detailed, and accurate road maps is the foundational demand of intelligent transportation systems and smart city management. Traditional road data acquisition methods, such as map digitization and field surveys, require expensive equipment and labor-intensive indoor processing. Presently, some collaborative mapping programs (e.g., OpenStreetMap, Wikimapia) have made a significant contribution to commercial road maps, but they also depend on the manual processing work of volunteers and are confronted with submissions of varying quality [1]. Widely used positioning devices (e.g., mobile phones) make it possible to create massive trajectory data for vehicles and people, which contain an abundance of information on road structures, traffic conditions, and driving behaviours. Hence, a surge in academic attention has, in recent years, been paid to automatically extracting road-related information from spatio-temporal trajectory data [2,3,4].

Due to the multipath effect of GNSS positioning and the inequality of human travelling activities, spatial–temporal heterogeneities probably occur in vehicle trajectory data (especially low-frequency data) [5]. On the one hand, the successive track points may pass through several street blocks due to a low sample rate or even satellite signal loss. On the other hand, the track points are intensive in some traffic hubs and high-level roads, but are sparse in other low-level roads. This brings about large challenges for accurately and efficiently extracting highly detailed road networks from heterogeneous trajectory data. Additionally, to the best of our knowledge, most existing studies are devoted to constructing a road map version using all the gathered trajectory data [6,7,8,9]. The version-based methods may ignore the road’s changing patterns, implied in multi-temporal trajectory data, and scarcely provide a time-sustainable data update mechanism. Hence, in this paper we explored the crucial cue of changes in road patterns in multi-temporal trajectory data and propose a hybrid method to incrementally extract road networks.

The remainder of this paper is organized as follows. The related works about road extraction from trajectory data are reviewed in Section 2. Section 3 presents the hybrid method for incrementally extracting road networks from spatio-temporal trajectory data. The experimental results and evaluation analysis are presented in Section 4 to demonstrate the efficiency and robustness of the proposed method. Finally, the conclusions and future works are discussed at the end of the paper.

2. Literature Review

In general, the existing road extraction methods comprise three categories: image-based [10], LiDAR-based [11], and trajectory-based methods [12]. Compared to trajectory-based methods, LiDAR-based and image-based methods are probably more precise, but are costly [13]. Recently, ubiquitous track data have become a dominant and area-wide data resource to monitor urban road environments [14]. On the one hand, the increasingly collected trajectory data can provide a good geometric delineation of the road network structure [15,16]. On the other hand, the recorded geographic information (e.g., time, position, direction, and speed) provides a valuable cue to infer semantic rules associated with road networks (e.g., traffic mode, turning restrictions, and speed limit) [17,18]. Moreover, continuously observed trajectory data can also be used for timely monitoring of dynamic traffic flow and changing road patterns [19,20]. This paper mainly concentrates on extracting structural road networks from trajectory data, and thus we undertake a brief review of trajectory-based road extraction models.

Benefiting from ubiquitous track data, many targeted models have been proposed to extract road networks from different levels, covering from low-detailed road information (central line, incline value, and surface roughness) to high-detailed road structures (e.g., intersection models, lane-level structure) [4,16,21,22]. In brief, the existing methods of trajectory-based road map reconstruction can be divided into four categories: rasterization methods, clustering methods, incremental methods, and other combined methods. The first rasterization methods convert vector trajectories to raster images and then extract the road’s central lines based on morphological thinning or skeletonization algorithms [23,24,25]. Rasterization methods can quickly obtain a road skeleton map, but are confronted with two critical problems: one is the probability of missing the topological relations between road segments, and the other is the difficulty in setting the optimal grid sizes in vector-to-grid conversation. The second clustering methods aggregate similar trajectory points or lines into different clusters associated with road lanes, segments, or intersections, and then delineate their geometric structures and turning rules by shape fitting algorithms and statistical analysis [26,27,28,29,30]. The clustering methods can efficiently recognize the topological relations and turning rules but urgently need to improve the insensitivity of handling abnormal trajectories and density disparities. The third incremental methods mainly adopt an iterative strategy, which match all the obtained trajectories in some time period with a reference road network (or an empty map) and then iteratively merge the target trajectories to renew the reference road networks [31,32,33,34]. The incremental methods make a good attempt at continuously updating road networks, but it is fairly time consuming to process the whole long-term observed trajectories [35]. Hence, some researchers proposed combined methods to generate structural graph models of road networks, such as the divide-and-conquer method [36,37], the structure learning method [38], the topological Morse theory [39], and the structure-aware filtering method [40].

Generally, previous studies carried out important exploratory research but also engendered some drawbacks to be solved, including vulnerability to spatio-temporal heterogeneities of crowdsourcing trajectory data as well as not fully exploring the potentially changing patterns implied in different temporal trajectories. Hence, the aim of this paper was to propose a time-sustainable model to incrementally extract urban road networks from continuously observed trajectory data. The contributions of this paper are mainly as follows: (1) the obtained trajectory data are separately partitioned into K time series for initializing K-temporal road networks, decreasing the time consumption of processing from all gathered trajectories; (2) the gravitational model was applied to geometrically adjust the K-temporal road networks, providing a changing road pattern discovery approach for road data renewal; and (3) curve fitting and rule mining models were developed to generate a structural road graph, improving the robustness of our method to accurately delineate various road shapes and topologies.

3. Method

According to the literature review in [12], traditional road inference methods directly deal with all gathered trajectory data to achieve a complete road map, which may ignore the road changing patterns implied in the trajectory data and may simultaneously require huge geometric computation to process all trajectory data. This paper thus proposes a hybrid method to incrementally extract urban road networks, which attempts to partition the original trajectory data into K temporal trajectory datasets and analyses the spatial–semantic inconsistencies between different time series. Generally, the proposed approach consists of three main steps, outlined as follows:

• Partitioning original trajectory data into K time series and initializing a road network version for each time series trajectory data based on morphological rasterization operations;
• Adjusting the K-temporal road networks, initialized above, according to a gravitation force model to amend the geometric deviations between different temporal road networks;
• Constructing urban road networks by a k-segment fitting algorithm and then, inferring the associated road parameters with a statistical trend analysis.

In general, a trajectory point records a variety of information, such as longitude, latitude, time, speed, heading, and license plate, which can be formally expressed as P_i = (lon_i, lat_i, time_i, speed_i, head_i, car_i,…). Compared to ordinary GPS data, DGPS data have the advantage of high-precision and dense sampling, but they also have a high cost and a long cycle of data collection. Therefore, ordinary GPS data are more conducive to rapid updates of urban road networks [32]. In view of the spatio-temporal heterogeneities of raw trajectory data (Figure 1a), trajectory filtering pre-processing is implemented as follows: (1) removing the track points with an abnormally large heading direction (e.g., 360º heading) and breaking the track lines with a few nearly stationary points; (2) removing the track lines with three or more large bending points (e.g., 50º turning); and (3) deleting the track lines with two or more abnormal distance points (e.g., more than 1.3km). The GPS points after pro-processing are shown in Figure 1b.

3.1. Initialization of K-temporal Road Networks Based on Morphological Rasterization

Spatio-temporal trajectory data provide a tireless and detailed perception of urban road networks. To mine the road change patterns implied in spatio-temporal trajectory data, we aggregated the entire trajectories into K temporal cliques and then extracted K-temporal road features from these trajectory cliques based on the morphological rasterization method.

In theory, morphological operations contain two basic operators (dilation and erosion), as well as some other derived operators (e.g., opening, closing, and thinning) [41]. As illustrated in Figure 2a, the dilation operator of binary image A by structuring element B₁ is defined as A⊕B₁ = {x, y|(B₁)_{x, y}∩A ≠ ∅}. Like that shown in Figure 2b, the erosion operator of binary image A by structuring element B₂ is defined as AΘB₂ = {x, y|(B₂)_{x, y}⊆A}. It can be seen that the dilation operator aims to expand the shapes of binary image A, and, conversely, the erosion operator aims to thin the shapes of binary image A and remove the crack edges.

The vector trajectory data in each temporal clique were first converted into binary raster images with two specified grid and density thresholds. After that, road features were then extracted by a dilation–opening–thinning operator strategy. First, to decrease the impact of coarse sampling, a dilation operation of structuring element B₁ was necessarily performed to expand the shape of the road regions and to connect the fractured road segments. Then, we executed an opening operator on the binary trajectory image after dilation. The opening operator is a derived operator which executes the basic erosion and dilation operators in order with the same structuring element B₂. As illustrated in Figure 3, Figure 3b is the result of directly performing an opening operator on the original trajectory image in Figure 3a, and Figure 3c is the result of performing the opening operator after a dilation operator on the original image. It can be seen that compared to directly performing the opening operator (Figure 3b), the geometric characteristics of the road segments are better preserved by adopting the dilation–opening operators (Figure 3c).

After that, a sophisticated thinning algorithm was utilized to generate road central lines from the trajectory images processed by the above dilation and opening operations. Image thinning, or skeletonization, is the process of obtaining the skeleton graph of objects from binary raster images. There are two representative thinning methods, i.e., iterative methods and non-iterative methods. In view of image noise immunity and preserving road completeness, we adopted the sophisticated iterative method, based on mathematical morphology in this paper. In detail, the target image A was iteratively thinned by a set of structuring elements {B₁,…,B_N} according to Equation (1). The morphological thinning operator was to set the pixel values as 0 when they are hit by the structuring elements B₁,…,B_N. The single-temporal road network can be obtained by inversely converting the thinned images to the vector skeleton map. Other temporal road networks are similarly extracted by the above mathematical morphology method.

$$A\otimes \left\{B_{1\dots N}\right\}=\left(\dots \left(\left(A\otimes B_1\right)\otimes B_2\right)\dots \right)\otimes B_N$$

(1)

3.2. Adjustment of K-temporal Road Networks Based on a Gravitation Force Model

Due to unstable fluctuations in positional accuracy and sampling rate, the road networks obtained from different temporal trajectory data show uncertain geometric deviations and potential road change information. This paper applied the gravitation force model to adjust the geometric deviation between the K-temporal road networks to achieve a complete and up-to-date road structure. The key idea of the gravitation force model is to iteratively compute the positional adjustment values by balancing the gravitation forces between K-temporal road central lines and the spring forces of moving road central lines.

We first utilized a simple buffering analysis to recognize similar road segments that are matched to an identical road entity, which are called K-temporal corresponding road segments. For example, in Figure 4.

R₁…R_i…R_j…R_K is one representative set of K-temporal corresponding road segments. According to the physical gravitation force model, the gravitation force of other road segments {R₁…R_j…R_K} (j ≠ i) to one road segment R_i at its vertex t, i.e., f₁(t, R_j), is defined by the following formula:

$$F_1\left(t\right)={\displaystyle {\displaystyle \sum }_{j=1,j\ne i}^K}-\frac{M\ast \left(K-1\right)}{{\delta }^3\sqrt{2\pi }}\cdot exp\left(-\frac{{\left|t-R_j\right|}^2}{2{\delta }^2}\right)\left(t-R_j\right),t\in R_i$$

(2)

where (t−R_j) indicates the projecting vector consisting of the vertex t and the closest point of R_j to t, and |t−R_j| measures the norm of the projecting vector (t−R_j). K is the total number of K-temporal corresponding road segments, and M and σ are the predefined parameters, set as M = 1 and σ = 10, respectively, according to [36].

In addition, the spring force of the moving the vertex t to a new position x is defined as Equation (3), where μ is the predefined spring force coefficient, set to 0.005 by default.

$$F_2\left(t\right)=\mu \cdot \left(x-t\right)$$

(3)

According to the physical gravitation force model, the optimal new position of each vertex t in one road segment R_i is calculated and updated iteratively by making the sum of the gravitation forces of other road segments {R₁…R_j…R_N} (j ≠ i) equal to the spring force of moving the vertex t to a new position x, i.e., F₁(t) = F₂(t). The positional adjustment process of that road segment R_i will terminate when the positional differences of all vertices between two iterations are smaller than a specified value. All K-temporal corresponding road segments are spatially adjusted by the same gravitation-force-based fusion process. Figure 5a shows some examples of K-temporal corresponding road segments extracted from K partitioned trajectory datasets, and Figure 5b shows those road segments after positional adjustment, which indicates that the geometric deviation between different temporal road segments is efficiently improved to generate a geometry-consistent road network.

3.3. Construction of Urban Road Networks Based on k-segment Fitting and Statistical Analysis

3.3.1. Geometry Construction of Urban Road Networks Based on the K-segment Fitting Algorithm

After the positional adjustment of K-temporal corresponding road networks, the same road entity still has multiple geometric representations. These K-temporal road geometries should be fitted as a principal curve to construct a final geometry representation for each road entity. The principal curve is defined as a smooth and self-consistent curve that passes through the distribution centre of all the observed points [42]. We developed a k-segment principal curve fitting algorithm to generate the final road networks in order for it to have a good performance in delineating various point distributions [4]. The k-segment fitting algorithm was designed to iteratively partition all observed points into different Voronoi sets and to find a smooth curve of k segments that satisfy the objective function in Equation (4):

$$\begin{array}{c}{\displaystyle \underset{f=\left\{s_1\dots s_k\right\}}{\mathrm{argmin}}}\left\{G\left(f\right)={{\displaystyle \sum }}_{i=1\dots k}{\displaystyle {\displaystyle \sum }_{p\in V_i}}d\left(p,s_i\right)\right\}\\ \mathbf{s}.\mathbf{t}.V_i=\left\{p\in P_n|i={\displaystyle \underset{j=1\dots k}{\mathrm{argmin}}}d\left(p,s_j\right)\right\}\end{array}$$

(4)

where P_n is a set of n observed points, p $\in$ P_n. In this paper, P_n indicates all the vertex points of one K-temporal corresponding road segment set. The curve f, consisting of k segments s₁…s_k, is the principal curve fitted by P_n. Additionally, d(p, s_j) indicates the minimum distance from point p to any points in the segment s_j. V₁…V_i…V_k denote the Voronoi partitions of the observed points, which means that the points in V_i must be closer to s_i than they are to any other segment.

It can be seen from Equation (4) that a minimum sum of distances is achieved between the observed points and the fitted curve f. As described in detail in the literature [4], a first principal line is initialized as f = {s₁}, k = 1; then, a new segment is iteratively inserted into f with k = k + 1, based on Equation (4), until the sum of the distances between observed point set P_n and the k segments reaches a steady value and finally, the k segments are linked as a smooth curve by the greedy Hamiltonian path search strategy. In this paper, the K-temporal road segments after positional adjustment were reclassified to different road entities and the vertex points in each K-temporal corresponding road segment set were regarded as an independent observed point set to fit a principal curve for geometrically delineating the corresponding road entity. It should be noted that the k in the principal curve fitting (Equation (4)) is different from the former K. K is the total number of trajectory cliques and k is the total segment number of the principal curve for delineating one road geometry.

3.3.2. Attribute Inferencing of Urban Road Networks Based on a Statistical Analysis

In addition to road geometries, road attributes, such as road width and turning directions, are also indispensable information for navigation and route planning applications. We reconstructed the linkage relations between target trajectory points and the fitted road segments by an adaptive buffer analysis, and further determined the road width and one-/two-way driving rules, according to the point distributions and heading direction statistics of the associated trajectory points. The semantic inference process consists of the following three steps.

Road Width Computation: For each road segment, we started to build its buffer region with a radius of ε₀ and compute the ratio r₀ of the trajectory points that fall in the buffer region buffer (ε₀); then, we rebuild the buffer region with an increasing radius of ε_i = ε₀ + i*s (where i is the iterative number and s = 0.1 m is the iterative step) to compute a new ratio r_i of trajectory points falling in the buffer (ε_i); the above buffer analysis continues until the trajectory point ratio r_i reaches 96%, and the buffer radius after iteratively buffering is inferred as the road width W (shown in Figure 6).

Driving rule identification: For each road segment, we computed the angles between the road segment and the associated trajectory lines in the buffer (W), which were used to classify all trajectory lines into three types according to Equation (5); then, the ratios of Type = 1 or 2 were calculated to judge the driving rules of that road segment. In this paper, when more than 90%, or less than 10%, of associated trajectory lines belong to Type = 1 or Type = 2, the road segment was regarded as a one-way road; otherwise, it was regarded as a two-way road.

$$Type=\{\begin{array}{c}\hfill 1 (same direction),\theta \in \left[0, \mathsf{\Phi }\right]\hfill \\ \hfill 2 (opposite direction), \theta \in \left[180-\mathsf{\Phi }, 180+\mathsf{\Phi }\right]\hfill \\ \hfill 0 (abnormal direction), otherwise\hfill \end{array}$$

(5)

Two-way road reconstruction: For the identified two-way roads, the geometry representations of dual carriageways are computed according to Equation (6):

$$P_i{}^{\prime }=P_i\pm \frac{\overrightarrow{W}}{6}$$

(6)

where P_i = (x_i, y_i) is the vertex coordinates of road central lines; $\overrightarrow{W}$ is a vector, of which the length equals the road width W and the direction is perpendicular to the road segment; and P_i′ = (x_i′, y_i′) is the corresponding vertex coordinates of the dual-carriageway representation.

Due to coarse sampling or two-way-road reconstruction, some generated roads may be disconnected or misconnected with other roads. In view of preserving road connectivity, a post-processing of topology reconstruction was necessarily implemented. First, to distinguish the dead-end roads that are generally far from other roads, we used a very small search radius (e.g., 1 m) to find the connected-end roads and to repair their topological relations with their neighbouring roads. As illustrated in Figure 7a, if only one neighbouring road at one node is found, the road node of L₁ will be extended straight to the unique neighbouring road (Figure 7b); otherwise, that road node of L₂ will be extended by inserting a new segment consisting of the nearest node and the road node itself (Figure 7b). Certainly, a fixed search radius may not be enough to judge dead-end roads and connected-end roads. It could be better to combine the spatial proximity of road segments and the topological connectivity of the associated trajectory lines to deal with these singleton segments.

4. Results

In the experimental analysis, one taxi trajectory dataset in Shenzhen, China, from the 1st to 21st January 2012, was used to verify our proposed method. The trajectory dataset has a sample rate of about 60–100 s, has a positional precision of about 10–15 m, and contains about 1.2 million trajectory points and 16,000 trajectories per day. The parameter settings for incrementally extracting road networks from the trajectory data are listed in

Table 1. Parameter settings for incrementally extracting urban road networks from trajectory data.

Parameter	B1	B2	M	δ	μ	kmax	ε0	s	Φ
Value	$\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 1& 1\end{array}\right]$	$\left[\begin{array}{ccc}0& 1& 0\\ 1& 1& 1\\ 0& 1& 0\end{array}\right]$	1.0	10.0	0.005	24	6 m	0.1 m	30°

. The structuring elements B₁ and B₂ in road network initialization (Section 3.1) and the buffering parameters ε₀ and s in road attribute inference (Section 3.3.2) were empirically determined for good performance to approach different temporal trajectory data. The parameters M, δ, and μ in the gravitation force model (Section 3.2) and kmax in k-segment fitting (Section 3.3.1) were uniformly set by referring to the literature in [36] and [4], respectively.

4.1. Experimental Results

Figure 8 illustrates the results of road network initialization, based on the morphological rasterization method. Figure 8a–c shows the binary images converted from three-temporal trajectory data and Figure 8d–f shows the road networks roughly extracted from the aforementioned three-temporal binary images. It can be seen from the red ellipses in Figure 8 that spatial inconsistencies and complementarities coexist in different temporal road networks because of spatio-temporal heterogeneities between trajectory data. The proposed K-temporal extraction strategy can efficiently discover the spatio-temporal road change characteristics that are implied in continuously observed trajectory data, providing significant clues for road change detection and self-updating. Moreover, the K-temporal road network initialization method can cause efficiency promotion, in comparison to directly extracting the road network from the whole trajectory dataset. In practice, it is a feasible solution to utilize parallel computation to further decrease the time consumption of K-temporal road network initialization.

Figure 9a,c depicts the results of overlaying the K-temporal road networks extracted by the rasterization method, and Figure 9b,d elaborates on the spatial adjustment results based on the gravitation force model. As shown in the blue rectangles in Figure 9, the geometric deviations between different temporal road networks are efficiently improved by the force-based fusion model, improving the geometric accuracy of the generated road networks.

After that, all K-temporal road segments were reclassified by a simple buffer search to make sure that each road segment cluster was associated with its identical road entity. As shown in Figure 10a, the road segments dotted in the same colour belong to one road segment cluster and are associated with the same road entity. It can be seen that each road segment cluster is one similarly shaped but data-compressed version of the raw trajectories, and the vertex points in each cluster are also one similar-shaped, but data-compressed, version of the raw trajectory points. In this way, we utilized the vertex points of each road segment cluster to generate the corresponding road central lines, based on the k-segment fitting algorithm. As shown in Figure 10b, the fitted road segments (especially some curve road segments) are structure-complete and geometry-consistent, providing a good delineation of the whole road skeleton. This demonstrates the good performance of our method in generating a detailed and up-to-date road map.

Based on the heading distribution of the associated trajectory points, the driving rule of each individual road segment can be identified as of two types, i.e., one-way and two-way. Figure 11a displays the identified one-way (grey line) and two-way (black line) road segments, indicating that a large proportion of road segments are two-way roads, and the one-way roads are mostly low-level streets or road intersection segments. Figure 11b shows the results of two-way road reconstruction. The experimental results illustrate that the trajectory data provide a rich and timely source to mine both the geometry and semantic information of urban road networks.

4.2. Result Analysis and Evaluation

To qualitatively evaluate the extraction precision of the proposed method, we overlaid the road networks extracted by our method and another three methods from Davies, Ahmed, and Wang together [27,34,39]. As depicted in the enlarged views of rectangles A and B in Figure 12, compared to the blue, green, and yellow road networks extracted by Davies, Ahmed, and Wang, higher spatial consistencies are clearly seen between our extracted road networks (red lines) and the background image. In particular, a geometric deviation of about 2–15 m can be seen between our extracted road segments and the reference remote sensing image (Google Earth).

In addition, an authorized road network was used as benchmark data to quantitatively evaluate the spatial consistency of the extracted road networks. In detail, a set of buffer regions of the benchmark roads was constructed using different buffer distances (e.g., 1–15 m). The respective length percentages of extracted roads that fell in the different buffer regions were calculated. As listed in

Table 2. Length percentage of extracted roads falling in different buffer regions of the reference road network.

Buffer Distance (Meters)	1	2	4	6	8	10	12	15
The proposed method	15.11%	30.02%	54.92%	70.61%	80.66%	87.32%	91.71%	96.32%
Davies’s Method [27]	23.26%	30.66%	43.93%	51.78%	56.08%	58.09%	59.08%	60.18%
Ahmed’s method [34]	05.89%	11.81%	22.26%	30.70%	37.46%	43.11%	47.94%	54.03%
Wang’s Method [39]	14.46%	16.26%	19.78%	23.16%	25.87%	28.42%	30.99%	35.71%

, more than 50% of the roads extracted by our method fell in the 4 m buffer regions of the reference roads, whereas only 43.93%, 22%, and 19.78% of the roads extracted by the methods in [27,34,39] fell in the same buffer regions. It can be seen that almost all roads extracted (96%) by our method are very close to the benchmark roads, within 15 m, showing a high spatial consistency with the ground truth.

Furthermore, to quantitively evaluate the precision of driving rules inferences, we cross-checked between the driving rules inferred by our method and the driving rules provided by official or commercial map providers (e.g., Google Maps and Baidu Maps). For the experimental case, we extracted 450 roads in total, of which 361 roads were inferred as two-way roads and only 89 roads were inferred as one-way roads. By comparing with the professional map services, 327 of the 361 extracted two-way roads were correctly defined as two-way roads and 65 of the 89 extracted one-way roads were correctly defined as one-way roads. The results indicate that 87.11% of the extracted roads are consistent with professional map services, demonstrating high robustness and accuracy of road detection and semantic extraction achieved by the proposed method.

In order to verify the efficiency increase of the K-temporal road extraction strategy, we calculated the elapsed time of both the force-based adjustment and k segment curve fitting based on raw trajectory data and K-temporal road networks initialized by the rasterization method. As listed in

Table 3. Elapsed time statistics of force-based adjustment and k segment curve fitting based on raw trajectory data and multi-temporal road networks.

	Force-Based Adjustment	k-Segment Curve Fitting
Raw trajectory data	6.02 s	1.48 s
K-temporal road networks	>5 h	7.68 s

, it took more than 5 h to accomplish the spatial adjustment of raw trajectory data, but 6.02 s was enough for the K-temporal road network adjustment. An obvious efficiency increase was also found in the k-segment curve fitting of different temporal road networks. This demonstrates that the extraction method based on K-temporal trajectory partitioning achieves a higher time efficiency than directly dealing with all trajectory data.

5. Conclusions

Ubiquitous trajectory data are a dominant and area-wide resource for road network extraction and updating. Presently, spatio-temporal heterogeneities in trajectory data are a significant obstacle to generating highly precise road maps. Since existing research has often ignored the potential road change patterns implied in different time-series trajectories, in this paper, we proposed a hybrid method to incrementally extract urban road networks from spatio-temporal trajectory data. The proposed method partitioned raw trajectory data into K time series and initially extracted K-temporal road networks from different trajectory cliques based on the rasterization method. The geometrical inconsistencies between the K-temporal road networks were then adjusted by a gravitation force model. After that, the final geometry delineation and associated road attributes were constructed by k-segment fitting and statistical analysis, respectively. The experimental results show that the proposed method can effectively process spatio-temporally heterogeneous trajectory data in order to generate a geometry-accurate and semantic-rich road network. Moreover, all trajectory data were partitioned into multiple time slices, which is a useful attempt at exploring road change patterns from area-wide and time-continuous trajectory data.

However, due to the large sampling heterogeneity, the proposed method still has some limitations in discovering low-level streets and complicated intersections. Secondly, the proposed method possibly detects some fake change information because of traffic flow fluctuations. Hence, future works will focus on improving the reliabilities of parameter settings to adaptively extract various road structures (e.g., low-level streets and road junctions). More significantly, a road fusion model incorporating traffic flow analysis and multi-temporal change validation needs to be developed to automatically filter fake change information and incrementally update road networks.