A Hierarchical Matching Method for Vectorial Road Networks Using Delaunay Triangulation

: The primary objective of vectorial road network matching is to identify homonymous roads from two di ﬀ erent data sources. Previous methods usually focus on matching road networks with the same coordinate system but rarely with di ﬀ erent or unknown coordinate systems, which may lead to nontrivial and nonsystematic deviations (e.g., rotation angle) between homonymous objects. To ﬁll this gap, this study proposes a novel hierarchical road network matching method based on Delaunay triangulation (DTRM). First, the entire urban road network is divided into three levels (L1, L2, L3) by using the principle of stroke. Then, the triangular meshes are constructed from L2, and the minimum matching unit (MMU) in the triangular mesh is used instead of the traditional “node-arc” unit to measure the similarity for the matching of L2. Lastly, a hierarchical matching solution integrating the probabilistic relaxation method and MMU similarity is yielded to identify the matching relationships of the three-level road network. Experiments conducted in Wuhan, China, and Auckland, New Zealand, show that the MMU similarity metrics can e ﬀ ectively calculate the similarity value with di ﬀ erent rotation angles, and DTRM has higher precision than the benchmark probability-relaxation-matching method (PRM) and can correctly identify the most matching-relationships with an average accuracy of 89.63%. This study provides a matching framework for road networks with di ﬀ erent or even unknown coordinate systems and contributes to the integration and updating of urban road networks. represented by the sum of the similarity values divided by the number of triangles. The correspondence between the nodes of the two MMUs is also obtained. comparison curve of DR between DRTM and PRM by using the 10 pairs of nodes. The average level of DTRM is signiﬁcantly higher than that of PRM. The average DR value of DTRM is 0.64, whereas that of PRM is 0.48. This result shows that the similarity values of the MMU are higher than those of PRM. This condition is conducive to the identiﬁcation of correct pairs of elements with the same name. above 90% of the candidate strokes have obvious geometrical di ﬀ erences with the target matched stroke. The HitRate value approaches 0 in the Dr range of 0–0.5, indicating there are very few candidate strokes that are very similar to the target matched stroke. In other words, the di ﬀ erence rate between the target matched strokes and the candidate strokes is signiﬁcant when the proportions are 5% and 10%, which is beneﬁcial to obtain correct matching results. When the values of Con p are larger than 10%, the HitRate value decreases in the Dr range of 0.5–1 and increases in the Dr range of 0–0.5. When Con b is set to 300 m, the distribution trend is similar to that when Con b is set to 50 m, indicating that the bu ﬀ er size has little inﬂuence on HitRate . Therefore, 10% is suggested as the selection proportion of the ﬁrst level L1 to ensure the robustness and accuracy of the matching.


Introduction
Geospatial data matching that connects different data sets and makes them interoperable has a wide range of applications in the spatial analysis due to the need to utilize information from different data sources [1]. On the one hand, large amounts of various geospatial data can be easily obtained from federal, state, and local governments or private companies. On the other hand, identifying corresponding objects is an essential step in change detection and continuous incremental updating. Thus, new features can easily be detected and heterogeneous data sets can be integrated into an enriched product by improving positional or semantic accuracies [2].
Vectorial road networks are digital representations of road maps. The main task of matching road networks is building the corresponding relationships of road-object pairs that represent the same segment of a real-world road in heterogeneous road maps [3,4]. It is an important prerequisite for road network integration, change detection, and data updating. The demand for matching road networks The remainder of this paper is organized as follows. Section 2 elaborates the proposed method for synthesizing the hierarchical matching strategy and similarity matching based on the MMU. Section 3 describes the experiments conducted using three different real datasets, analyzes the results, and discusses factors that may affect the performance of the proposed method. Section 4 concludes this study.

Methodological Framework
As shown in Figure 1, the proposed framework synthesizes the hierarchical matching strategy and MMU similarity matching; it is a four-step process. First, the framework divides the entire urban road network (Ω) into three levels: the skeleton roads (L1: Υ), the frame roads of subregions (L2: Φ), and the remaining part of the subregion roads excluding the frame roads (L3: χ = Ω − Υ − Φ). Then, L2 roads are converted into triangular meshes constrained by road nodes and segments by using constrained Delaunay triangulation (CDT). Subsequently, the similarity of MMU in the meshes is calculated to identify the matching relationship of L2. Lastly, a hierarchical strategy is applied to identify the matching relationship of the three levels of the road network, especially by using the matched relationship in L2 to derive the relationship of the element pairs in L3. Characterizing the geometric and topological properties of roads is crucial for matching processes [21]. In this study, the MMU in the triangular meshes is introduced and interpreted as the node-area structure to tackle the problem of nonsystematic bias generated by different coordinate systems. Meanwhile, to consider different scales, the MMU similarity is embedded into the hierarchical matching framework, which then drives the matching applied at the nonskeletal level.
The following subsections provide detailed descriptions in accordance with the four steps listed in this framework. The following subsections provide detailed descriptions in accordance with the four steps listed in this framework.

Hierarchical Generation of Road Network
The hierarchical generation process for road network consists of three steps. First, by simulating the human visual cognitive mechanism, the skeleton of the urban road network is extracted using the principle of stroke [22]. Here, stroke is defined as the natural functional units of a network [23]. It is constructed by aggregating road segments in accordance with different properties, such as street names or the angle among neighboring road segments, which is known as the principle of stroke. The longest k strokes are selected as L1. The value of k is determined by a proportion of 10% and how the proportion is determined is discussed in Section 3.7.
The roads with the node connectivity of 1 are selected as nonskeletal roads (L3). Let NS be the number of strokes, LS be the length of a stroke, and T be the threshold of stroke length. ∆NS, the difference between the number of strokes in the source and target datasets whose length is greater than T, is defined as where NS OT denotes the number of strokes with a length greater than T in the source dataset, NS DT denotes the number of strokes with a length greater than T in the target dataset, and S T is a collection of all road lengths in the source and target datasets. The minimum threshold of the stroke length is represented as C and formulated as The strokes whose lengths are less than C are also classified as L3, whereas the remaining roads are classified as L2. The difference between the number of strokes with a length greater than C in the source and target datasets should be 0, or the derivative of the difference should be 0. Therefore, the numbers of strokes obtained from L2 in the source and target datasets are similar. This condition is a prerequisite for matching.
As shown in Figure 2, according to the principle of hierarchical generation, an example of a road network is first divided into subregions which form L1. L1 represents the skeleton of the road network. Then, for the upper-right subregion in L1, the road network of L2 and L3 generated from the source road network and the target road network are illustrated, respectively. The difference between the source road network and the target road network is subtle for L1, moderate for L2, and obvious for L3.

Delaunay Triangulation Construction Constrained by Road Nodes and Segments
In this subsection, L2 that is composed of nodes and road segments is converted into a triangular mesh. As shown in Figure 3a, the upper-left corner shows the node-arc data to be triangulated, including the six road nodes of A, B, C, D, E, and F, and a road segment AB. In the absence of segment AB, the triangulation is constructed as shown in Figure 3b. In this study, the

Delaunay Triangulation Construction Constrained by Road Nodes and Segments
In this subsection, L2 that is composed of nodes and road segments is converted into a triangular mesh. As shown in Figure 3a, the upper-left corner shows the node-arc data to be triangulated, including the six road nodes of A, B, C, D, E, and F, and a road segment AB. In the absence of segment AB, the triangulation is constructed as shown in Figure 3b. In this study, the triangulation needs to be performed under the constraints of nodes and segments, as shown in Figure 3c.

Delaunay Triangulation Construction Constrained by Road Nodes and Segments
In this subsection, L2 that is composed of nodes and road segments is converted into a triangular mesh. As shown in Figure 3a, the upper-left corner shows the node-arc data to be triangulated, including the six road nodes of A, B, C, D, E, and F, and a road segment AB. In the absence of segment AB, the triangulation is constructed as shown in Figure 3b. In this study, the triangulation needs to be performed under the constraints of nodes and segments, as shown in Figure 3c. When many concave roads exist in the boundary of the street block, as depicted in Figure 4a, the generated triangular mesh may have some narrow and elongated triangles, as shown in Figure 4b. In Figure 5a, some nodes inside the boundary have a triangular mesh on only one side because the boundary serves as an external constraint. This biased mesh exerts a great influence on the subsequent matching process. To mitigate the boundary effect, the triangular meshes are optimized as follows: (1) The convex hull (M) of all nodes in the frame roads of the source and target datasets is calculated.
(2) The convex hull M is extended using a given distance threshold to obtain a new convex hull N that is parallel to M. (3) The nodes of all line segments on the convex hull N and the nodes and segments of the frame roads in the road network dataset are triangulated together.
The optimized triangulations by applying the above steps are illustrated in Figure 5b. Hereafter, the mesh structure for road network matching is constructed. When many concave roads exist in the boundary of the street block, as depicted in Figure 4a, the generated triangular mesh may have some narrow and elongated triangles, as shown in Figure  4b. In Figure 5a, some nodes inside the boundary have a triangular mesh on only one side because the boundary serves as an external constraint. This biased mesh exerts a great influence on the subsequent matching process. To mitigate the boundary effect, the triangular meshes are optimized as follows: (1) The convex hull (M) of all nodes in the frame roads of the source and target datasets is calculated. (2) The convex hull M is extended using a given distance threshold to obtain a new convex hull N that is parallel to M. (3) The nodes of all line segments on the convex hull N and the nodes and segments of the frame roads in the road network dataset are triangulated together.
The optimized triangulations by applying the above steps are illustrated in Figure 5b. Hereafter, the mesh structure for road network matching is constructed.

Similarity Metrics based on MMU
In this subsection, an MMU is constructed as the node-area structure. The similarity among MMUs depends on the similarity between basic triangle pairs. The triangular structure is continuous in space, and its similarity is unaffected by the rotation angle. It can efficiently characterize the spatial structure of the road network.
The MMU is defined as a set of triangles with the same vertices and is formulated as

Similarity Metrics Based on MMU
In this subsection, an MMU is constructed as the node-area structure. The similarity among MMUs depends on the similarity between basic triangle pairs. The triangular structure is continuous in space, and its similarity is unaffected by the rotation angle. It can efficiently characterize the spatial structure of the road network.
The MMU is defined as a set of triangles with the same vertices and is formulated as where CP represents the center node of the MMU. Let M be a mesh model, M = {V, E, T}, where V is a set of vertices, E is a set of edges, and T is a set of triangles. A triangle t, t = {v i , e i ; i = 0, 1, 2}, comprises three vertices and three edges. An edge e i , e i = {v p , v q , T i }, consists of endpoints v p and v q and the set of triangles T i that is bounded by it. A vertex v i , v i = {x, y, E i , T i }, is denoted with its coordinates and the edges E i and triangles T i that contain it. TS represents the collection of triangles T i originated from CP, which is formulated as TS = {T I , CP ∈ T I and i = 1, 2, . . . , n} where each triangle T i has its central node CP, and these n triangles are stored in a clockwise arrangement. For each triangle T, we first store the central node CP and then its two vertices in a clockwise order; it is represented as Figure 6 shows an MMU, in which the set of triangles {a, b, c, d, e} centered at point P are surrounded by points A, B, C, D, and E.
TSM P represents the set of MMUs generated using all nodes in the MMU with P as the center node and is defined as Formula (6). S P represents the set of all points of MMU P .
The similarity between MMU P and MMU O is calculated as ISPRS Int. J. Geo-Inf. 2020, 9, 509 8 of 31 For triangle T i (∆ABC) in TS O and triangle T j (∆DEF) in TS P , vertices A and D, vertices B and E, and vertices C and F correspond to each other. The similarity between the two angles ∠A (set to a) and ∠D (set to x) is , b = aP/3, and P is generally set to 50%. For triangles T i and T j , their similarity Sim ij is calculated as follows: Let P be the node in the source dataset and O be the candidate that matches P in the target dataset, then the similarity between MMU P and MMU O is calculated via the following five steps: (1) MMU P with P as the center node and TSM O with O as the center node are prepared.
(2) MMU P centering on P is moved to O to obtain a new MMU, namely, MMU P .
(3) The distance between all nodes in MMU P and TSM O is calculated to form a distance matrix.
We use the distance matrix to select the nodes in TSM O , in which the sum of the distance between all the nodes in MMU P and them is the smallest, which is called CPS O . The selected node cannot be repeated.
where each triangle Ti has its central node CP, and these n triangles are stored in a clockwise arrangement. For each triangle T, we first store the central node CP and then its two vertices in a clockwise order; it is represented as  TSMP represents the set of MMUs generated using all nodes in the MMU with P as the center node and is defined as Formula (6). SP represents the set of all points of MMUP.
The similarity between MMUP and MMUO is calculated as For triangle Ti (Δ ) in and triangle Tj (Δ ) in , vertices A and D, vertices B and E,   The correspondence between the nodes of the two MMUs is also obtained. Figure 7 shows an example of TSMP and TSMO, in which the red lines represent MMUP and MMUO, respectively. After moving point P to point O, MMUP also moves with P. The distance matrix of all nodes in MMUP and all nodes in TSMO is calculated, and a set of nodes with the smallest sum of distances from all nodes in MMUP is obtained.  Figure 8. Then, the similarity between all triangles in MMUP and MMUO' is calculated to acquire a triangle similarity matrix, as shown in Table 1. The sum of similarity values on all diagonals in the similarity matrix is calculated to generate the largest one. We characterize the average similarity of the MMU by using the sum of diagonal similarities of the matrix. For example, the sum of similarity of the diagonal line (a-v, b-z, c-y, d-x, e-w) is calculated as 4.12, then the similarity between the two MMUs of O and P is 4.12/5 = 0.824.   Figure 8. Then, the similarity between all triangles in MMU P and MMU O is calculated to acquire a triangle similarity matrix, as shown in Table 1. The sum of similarity values on all diagonals in the similarity matrix is calculated to generate the largest one. We characterize the average similarity of the MMU by using the sum of diagonal similarities of the matrix. For example, the sum of similarity of the diagonal line (a-v, b-z, c-y, d-x, e-w) is calculated as 4.12, then the similarity between the two MMUs of O and P is 4.12/5 = 0.824. MMUO, respectively. After moving point P to point O, MMUP also moves with P. The distance matrix of all nodes in MMUP and all nodes in TSMO is calculated, and a set of nodes with the smallest sum of distances from all nodes in MMUP is obtained.  Figure 8. Then, the similarity between all triangles in MMUP and MMUO' is calculated to acquire a triangle similarity matrix, as shown in Table 1. The sum of similarity values on all diagonals in the similarity matrix is calculated to generate the largest one. We characterize the average similarity of the MMU by using the sum of diagonal similarities of the matrix. For example, the sum of similarity of the diagonal line (a-v, b-z, c-y, d-x, e-w) is calculated as 4.12, then the similarity between the two MMUs of O and P is 4.12/5 = 0.824.   Based on the similarity among MMUs, the candidates matched to each MMU can be determined. Figure 9 shows that after the previous calculations, the correspondence between each node, that is, the correspondence between (A, S), (B, T), (C, M), (D, N), and (E, R), and the similarity between their MMUs can be obtained. Therefore, the similarity between the neighboring MMUs of O and P is utilized to optimize the similarity of O and P and obtain a similarity calculation result that is globally optimal.
Probability relaxation iteration is performed on the MMU similarity. In each iteration, the matching probability of MMU O and MMU P needs to be adjusted using the matching probability (i.e., similarity) of the triangle pairs of its subordinates. The similarity between the pairs of triangles <i,j> in the MMUs in the t-th iteration Sim ij can be updated using the weighted sum of the probability of this iteration and the average matching probability of the triangle pairs under it. Thereby, the similarity between O and P is updated as where n is the number of triangle pairs in the two MMUs, in which the value is an integer; the larger the value is, the faster the convergence will be. Variables i and j correspond to the sequence numbers in a clockwise order of the triangles in MMU O and MMU P , respectively. The above probability matrix iterates until convergence. The pair with the highest matching probability is the final matching pair.
between their MMUs can be obtained. Therefore, the similarity between the neighboring MMUs of O and P is utilized to optimize the similarity of O and P and obtain a similarity calculation result that is globally optimal. Probability relaxation iteration is performed on the MMU similarity. In each iteration, the matching probability of MMUO and MMUP needs to be adjusted using the matching probability (i.e., similarity) of the triangle pairs of its subordinates. The similarity between the pairs of triangles <i,j> in the MMUs in the t-th iteration Simij can be updated using the weighted sum of the probability of this iteration and the average matching probability of the triangle pairs under it. Thereby, the similarity between O and P is updated as where n is the number of triangle pairs in the two MMUs, in which the value is an integer; the larger the value is, the faster the convergence will be. Variables i and j correspond to the sequence numbers in a clockwise order of the triangles in MMUO and MMUP, respectively. The above probability matrix iterates until convergence. The pair with the highest matching probability is the final matching pair.
(a) (b) The stop conditions of the iteration are twofold: (1) The number of iterations reaches the specified threshold (usually 20); (2) The difference of the current similarity value and the last iteration similarity value is less than the specified threshold (usually 0.05). During the iterative process, the similarity of correct matches will decrease slowly, whereas the similarity of error matches will decrease rapidly. A clear gap exists in the similarity between correctly and wrongly matched ones, then the final matching relationship among MMUs is determined using the optimized similarity. Figure 10 shows a certain rotation offset between two MMUs with O and P as the central nodes. To handle the matching tasks with the rotation angle, the following processing should be performed:  2) The difference of the current similarity value and the last iteration similarity value is less than the specified threshold (usually 0.05). During the iterative process, the similarity of correct matches will decrease slowly, whereas the similarity of error matches will decrease rapidly. A clear gap exists in the similarity between correctly and wrongly matched ones, then the final matching relationship among MMUs is determined using the optimized similarity. Figure 10 shows a certain rotation offset between two MMUs with O and P as the central nodes. To handle the matching tasks with the rotation angle, the following processing should be performed: (1) The similarity of all triangle pairs in the two MMUs is calculated to obtain a similarity matrix M.
(2) The two largest similarity values on each diagonal of the similarity matrix, M ij , (3) M k ij , are selected, and the sum of the two values, which is denoted as Sim ij , is calculated. (4) The largest similarity values Sim ij , which indicate the similarity values on the i-th diagonal, are selected as the best indicator for characterizing the matching relationships between the triangles in the two MMUs. (5) The maximum similarity value M ij on the i-th diagonal is chosen, and the deflection angle δ of the two triangles is calculated in accordance with the correspondence between the triangle pair <i,j>. Then, the MMU centered at P is rotated using the angle δ to generate a new MMU, as shown in Figure 11. Lastly, the similarities between the two MMUs are calculated.
(5) The maximum similarity value on the i-th diagonal is chosen, and the deflection angle δ of the two triangles is calculated in accordance with the correspondence between the triangle pair <i,j>. Then, the MMU centered at P is rotated using the angle δ to generate a new MMU, as shown in Figure 11. Lastly, the similarities between the two MMUs are calculated.

Identification of the Matching Relationship by Using the Hierarchical Matching Strategy
The spatial relationship of neighboring roads is relatively stable during stratification. The road structure that constitutes the urban skeleton can be used as a reliable hierarchical matching context for the local matching units. The matching relationships are identified from high level to low level.
First, the matching relationship of L1 is identified. The common part of the skeletons is constructed to determine certain matches of L1.
Second, the matching relationship of L2 is identified. Let Va {ai, i∈{1,2,3,...,m}} and Vb {bi, i∈{1,2,3,...,n}} be the node sets of L2 in the source and target datasets. Va has m nodes. Vb has n nodes. ai represents the i-th node in Va, and bi represents the i-th node in Vb (bi, i∈{1,2,3,...,n}). On the basis of the MMU similarity metrics in Subsections 2.2-2.4, the steps for matching L2 are as follows: (1) The convex hull M of all nodes in Va and Vb is calculated and extended using a certain distance (here, the distance value is set to 90 m) to obtain a new convex hull N parallel to M.  (5) The maximum similarity value on the i-th diagonal is chosen, and the deflection angle δ of the two triangles is calculated in accordance with the correspondence between the triangle pair <i,j>. Then, the MMU centered at P is rotated using the angle δ to generate a new MMU, as shown in Figure 11. Lastly, the similarities between the two MMUs are calculated.

Identification of the Matching Relationship by Using the Hierarchical Matching Strategy
The spatial relationship of neighboring roads is relatively stable during stratification. The road structure that constitutes the urban skeleton can be used as a reliable hierarchical matching context for the local matching units. The matching relationships are identified from high level to low level.
First, the matching relationship of L1 is identified. The common part of the skeletons is constructed to determine certain matches of L1.
Second, the matching relationship of L2 is identified. Let Va {ai, i∈{1,2,3,...,m}} and Vb {bi, i∈{1,2,3,...,n}} be the node sets of L2 in the source and target datasets. Va has m nodes. Vb has n nodes. ai represents the i-th node in Va, and bi represents the i-th node in Vb (bi, i∈{1,2,3,...,n}). On the basis of the MMU similarity metrics in Subsections 2.2-2.4, the steps for matching L2 are as follows: (1) The convex hull M of all nodes in Va and Vb is calculated and extended using a certain distance (here, the distance value is set to 90 m) to obtain a new convex hull N parallel to M.

Identification of the Matching Relationship by Using the Hierarchical Matching Strategy
The spatial relationship of neighboring roads is relatively stable during stratification. The road structure that constitutes the urban skeleton can be used as a reliable hierarchical matching context for the local matching units. The matching relationships are identified from high level to low level.
First, the matching relationship of L1 is identified. The common part of the skeletons is constructed to determine certain matches of L1.
Second, the matching relationship of L2 is identified. Let V a {a i , i∈{1,2,3, . . . ,m}} and V b {b i , i∈{1,2,3, . . . ,n}} be the node sets of L2 in the source and target datasets. V a has m nodes. V b has n nodes. a i represents the i-th node in V a , and b i represents the i-th node in V b (b i , i∈{1,2,3, . . . ,n}). On the basis of the MMU similarity metrics in Sections 2.2-2.4, the steps for matching L2 are as follows: (1) The convex hull M of all nodes in V a and V b is calculated and extended using a certain distance (here, the distance value is set to 90 m) to obtain a new convex hull N parallel to M. Third, the matching relationship of L3 is identified. L3 is usually an extension of L2; hence, the connections between L2 and L3 are used as the context to derive the matching of L3 on the basis of the matching of L2. Let L i be the i-th road of the source dataset in L2 and L j be the j-th road of the target dataset in L2. The steps for matching L3 are as follows: (1) Each L i has two endpoints, Vi 1 and Vi 2 . The connected road in L3 is denoted as CL i with two endpoints CVi 1 , CVi 2 , and CVi 1 is on L i . Similarly, each L j has two endpoints, Vj 1 and Vj 2 . (2) In accordance with the matching results of L2, if the matching pairs of (Vi 1 , Vj 1 ) and (Vi 2 , Vj 2 ) are obtained, then the matching relationship between (L i , L j ) can be determined. The road connected to L j is denoted as CLj (j = 1,2, 3,...,k). CVj 1 , CVj 2 are the two endpoints of L j, and CVj 1 is on L j . (3) If CL i and CL j meet the following two conditions: 1) CL i and CL j are on the same side of L i and L j , respectively; 2) The relative position of CVi 1 in L i is closest to the relative position of CVj 1 in L j , and the difference should not be greater than 20%, then CVi 1 matches CVj 1 . (4) If CVi 2 and CVj 2 are not connected to other roads, then CVi 2 matches CVj 2 , and the matched nodes are placed into the same road layer. (5) Steps 1-4 are repeated until no new matches exist.
As illustrated in Figure 12, road AB is from L2 of the source dataset, and road QW is from L2 of the target dataset. GF, CD, and EH are from L3 of the source dataset, and TR is from L3 of the target dataset. Given that (A, Q) and (B, W) are known matching nodes, AB matches QW. Three roads GF, EH, and CD are connected to road AB, where nodes F, H, and D are on road AB. Road TR is connected to road QW, including R. The relative positions of points F, H, and D on road AB are 25%, 30%, and 76%, respectively; and R is located at 28% of road QW. Roads GF and CD are located on the left side of road AB, road EH is located on the right side of road AB, and road TR is located on the left side of road QW. GF and TR satisfy the two conditions in step 3. Therefore, F matches R. In accordance with step 4, G matches T for G and T are not connected to other roads.
(1) Each Li has two endpoints, Vi1 and Vi2. The connected road in L3 is denoted as CLi with two endpoints CVi1, CVi2, and CVi1 is on Li. Similarly, each Lj has two endpoints, Vj1 and Vj2. (2) In accordance with the matching results of L2, if the matching pairs of (Vi1, Vj1) and (Vi2, Vj2) are obtained, then the matching relationship between (Li, Lj) can be determined. The road connected to Lj is denoted as CLj (j = 1,2, 3,...,k). CVj1, CVj2 are the two endpoints of Lj, and CVj1 is on Lj. As illustrated in Figure 12, road AB is from L2 of the source dataset, and road QW is from L2 of the target dataset. GF, CD, and EH are from L3 of the source dataset, and TR is from L3 of the target dataset. Given that (A, Q) and (B, W) are known matching nodes, AB matches QW. Three roads GF, EH, and CD are connected to road AB, where nodes F, H, and D are on road AB. Road TR is connected to road QW, including R. The relative positions of points F, H, and D on road AB are 25%, 30%, and 76%, respectively; and R is located at 28% of road QW. Roads GF and CD are located on the left side of road AB, road EH is located on the right side of road AB, and road TR is located on the left side of road QW. GF and TR satisfy the two conditions in step 3. Therefore, F matches R. In accordance with step 4, G matches T for G and T are not connected to other roads.

Experimental Area and Data
The proposed method is implemented using the C++ programming language. Experiments are conducted on a computer equipped with an NVIDIA GeForce GTX 960 GPU and an Intel(R) Pentium(R) D 3.00GHz with 3.5 GB RAM. The software-development environment is composed of Visual Studio 2010 and MapGIS 10 (Zondy, Wuhan, Hubei, China).

Experimental Area and Data
The proposed method is implemented using the C++ programming language. Experiments are conducted on a computer equipped with an NVIDIA GeForce GTX 960 GPU and an Intel(R) Pentium(R) D 3.00GHz with 3.5 GB RAM. The software-development environment is composed of Visual Studio 2010 and MapGIS 10 (Zondy, Wuhan, Hubei, China).
To evaluate the effectiveness and performance of our method, three pairs of road networks are utilized. Areas 1 and 2 are from Wuhan, China and Area 3 belongs to Auckland, New Zealand. Each road network pair comes from different data producers. The source map represents the road network to be matched, and the target map represents the road network that is used to match the source map.
First, the selection of experimental data takes into account the influence of scale. The source and target maps of Area 1 have the same scales, whereas the maps for Area 2 and Area 3 have different scales. The three areas are illustrated in Figures 13-15 in a way of spatial overlay. Figure 13 shows the selected partial road network of Wuhan with the same scale, in which the difference between the source and target data is small. Figures 14 and 15 show the urban road data of Wuhan or Auckland with different scales. The source and target datasets are consistent in the global structure; in local areas, the target dataset has smaller roads. That is, a large portion of 1:0 relationships exist, and the differences in these data may cause a remarkable difference in the basic matching unit-MMU and certain challenges to the matching algorithm. The information on road nodes, road segments, and strokes in the three experimental areas are shown in Table 2. utilized. Areas 1 and 2 are from Wuhan, China and Area 3 belongs to Auckland, New Zealand. Each road network pair comes from different data producers. The source map represents the road network to be matched, and the target map represents the road network that is used to match the source map.
First, the selection of experimental data takes into account the influence of scale. The source and target maps of Area 1 have the same scales, whereas the maps for Area 2 and Area 3 have different scales. The three areas are illustrated in Figures 13-15 in a way of spatial overlay. Figure 13 shows the selected partial road network of Wuhan with the same scale, in which the difference between the source and target data is small. Figures 14 and 15 show the urban road data of Wuhan or Auckland with different scales. The source and target datasets are consistent in the global structure; in local areas, the target dataset has smaller roads. That is, a large portion of 1:0 relationships exist, and the differences in these data may cause a remarkable difference in the basic matching unit-MMU and certain challenges to the matching algorithm. The information on road nodes, road segments, and strokes in the three experimental areas are shown in Table 2.   road network pair comes from different data producers. The source map represents the road network to be matched, and the target map represents the road network that is used to match the source map.
First, the selection of experimental data takes into account the influence of scale. The source and target maps of Area 1 have the same scales, whereas the maps for Area 2 and Area 3 have different scales. The three areas are illustrated in Figures 13-15 in a way of spatial overlay. Figure 13 shows the selected partial road network of Wuhan with the same scale, in which the difference between the source and target data is small. Figures 14 and 15 show the urban road data of Wuhan or Auckland with different scales. The source and target datasets are consistent in the global structure; in local areas, the target dataset has smaller roads. That is, a large portion of 1:0 relationships exist, and the differences in these data may cause a remarkable difference in the basic matching unit-MMU and certain challenges to the matching algorithm. The information on road nodes, road segments, and strokes in the three experimental areas are shown in Table 2.       Area1_Src  107  68  31  Area1_Des  109  69  32  Area2_Src  715  431  148  Area2_Des  2083  1486  463  Area3_Src  608  386  134  Area3_Des  2805  1574  764 Second, street patterns are also considered when choosing experimental data. The patterns of urban streets are the result of many interacting phenomena over time, such as geography, natural setting, and socio-economic transformation. In this process, several typical street patterns are gradually formed: grid-like (planned) street networks, irregular (self-evolved) street networks, and hybrid street networks (planned and self-evolved structures coexist) [24,25]. Hierarchical partitioning should be conducted in the three experimental areas by using our proposed method. For Area 1, hierarchical partitioning is skipped and the matching calculation is directly performed because of its small size. Area 1 is numbered as R1. Area 2 is divided into four subregions, numbered R2, R3, R4, and R5. Area 3 is divided into three subregions, numbered R6, R7, and R8. The eight subregions are shown in Figure 16. R1, R2, R3, and R4 belong to irregular (self-evolved) street networks; R5 and R7 are mainly grid-like (planned) street networks; R6 and R8 are mostly hybrid street networks. The eight subregions cover three typical patterns of urban streets; therefore, the effectiveness of our proposed model is examined with these eight subregions.   Area1_Src  107  68  31  Area1_Des  109  69  32  Area2_Src  715  431  148  Area2_Des  2,083  1,486  463  Area3_Src  608  386  134  Area3_Des 2,805 1,574 764 Second, street patterns are also considered when choosing experimental data. The patterns of urban streets are the result of many interacting phenomena over time, such as geography, natural setting, and socio-economic transformation. In this process, several typical street patterns are gradually formed: grid-like (planned) street networks, irregular (self-evolved) street networks, and hybrid street networks (planned and self-evolved structures coexist) [24,25]. Hierarchical partitioning should be conducted in the three experimental areas by using our proposed method. For Area 1, hierarchical partitioning is skipped and the matching calculation is directly performed because of its small size. Area 1 is numbered as R1. Area 2 is divided into four subregions, numbered R2, R3, R4, and R5. Area 3 is divided into three subregions, numbered R6, R7, and R8. The eight subregions are shown in Figure 16. R1, R2, R3, and R4 belong to irregular (self-evolved) street networks; R5 and R7 are mainly grid-like (planned) street networks; R6 and R8 are mostly hybrid street networks. The eight subregions cover three typical patterns of urban streets; therefore, the effectiveness of our proposed model is examined with these eight subregions.

Model Evaluation Indices
Three evaluation indices are introduced to evaluate the matching results of the models quantitatively. They are precision (P), recall (R), and F1-score (F1).
Precision represents the matching accuracy in object matching. It is defined as the percentage of correctly matched pairs concerning the total number of matched pairs and represented as where TP (true positive) stands for the number of road pairs that are correctly matched, and FP (false positive) stands for the number of road pairs that are incorrectly matched. The closer the P-value is to 1, the more accurate the matching features identified by the algorithm.
Recall is the percentage of correctly matched pairs concerning the number of real matching pairs and is represented as where FN (false negative) stands for the number of actually corresponding pairs in two data sets that are not detected. F1-score combines precision and recall to provide a single metric to show the real success of our model in comparison with PRM and is expressed as follows: The effectiveness of the proposed algorithm is evaluated by comparing it with a benchmark. In recent years, the probability-relaxation-matching (PRM) algorithm has been employed in road network matching studies [2,7,[12][13][14]. Thus, as one of the representative algorithms, the matching method (henceforth PRM) proposed by [7] is adopted as the benchmark and implemented under the same experimental environment in this study. In the following experiments, for convenience, our model is called DTRM. The ending condition for the iteration of PRM is set on the basis of the experience value (0.0005) [26].

Evaluation of Algorithm Performance
We evaluate the algorithm performance with the three experimental areas (Areas 1-3). Specifically, our proposed algorithm is used to match the eight subregions and compare the result with the manually identified matching relationship to evaluate the algorithm precision. Tables 3-5 present the statistical results.  Figure 17 delineates the F1-score of DTRM and the benchmark PRM with different road networks. The number on the horizontal axis corresponds to the eight subregions (from R1 to R8). The F1-score of DTRM is greater than 85%, and more than half of the F1-score is greater than 93%, peaking at 98.68% for R1.  Figure 17 delineates the F1-score of DTRM and the benchmark PRM with different road networks. The number on the horizontal axis corresponds to the eight subregions (from R1 to R8). The F1-score of DTRM is greater than 85%, and more than half of the F1-score is greater than 93%, peaking at 98.68% for R1. Figure 17. F1-score of DTRM and PRM in the eight areas (areas 1-8). Table 3, the F1-score of DTRM for Area 1 (R1) reaches 98.68%, which shows that DTRM performs well in small-scale road data at the same scale and can accurately obtain the matching results. In Table 4, the F1-score of DTRM for Area 2 (partial roads in Wuhan) is above 93%, (except for R2 which is slightly lower than 88.23%). Compared with PRM, the F1-score of the three subregions of R2, R3, and R4 is higher than that of PRM, and the F1-score of R5 is slightly lower than that of PRM. In Table 5, the F1-score of DTRM is approximately 86% in the three subregions of Area 3 (Auckland). It indicates that the performance of DTRM in the three subregions is almost the same, without much fluctuation. By contrast, the performance of PRM in the three subregions is inconsistent. The F1-score in R8 is only 82.96%, which is much lower than that of DTRM.

As listed in
Overall, the F1-score of DTRM varies from region to region, mainly following the same trend as the benchmark. The majority of subregions (75%) have higher F1-scores when using DTRM compared with the situation when using PRM. Hence, the performance of DTRM has been improved to a certain extent. Figure 18 shows the difference between the F1-scores of DTRM and PRM. In the eight experimental subregions, the average F1-score difference between DTRM and PRM is 1.25 (the maximum value is 3.91% and the minimum value is 1.07%). The F1-score difference of six of these subregions is positive, with an average value of 2.11. The F1-score difference between the two subregions (R5 and R7) is negative. For the subregions with negative values, the absolute values are less than 0.02, which is in an acceptable error range. Figure A2d and Figure A3b show that the error matches (red lines) are mainly distributed in the non-grid area of narrow areas. The reason for the precision decrease of the two cases could be explained as follows: For such narrow areas, the matching relationship established in L1 is limited. In the complex hybrid road matching, there are As listed in Table 3, the F1-score of DTRM for Area 1 (R1) reaches 98.68%, which shows that DTRM performs well in small-scale road data at the same scale and can accurately obtain the matching results. In Table 4, the F1-score of DTRM for Area 2 (partial roads in Wuhan) is above 93%, (except for R2 which is slightly lower than 88.23%). Compared with PRM, the F1-score of the three subregions of R2, R3, and R4 is higher than that of PRM, and the F1-score of R5 is slightly lower than that of PRM. In Table 5, the F1-score of DTRM is approximately 86% in the three subregions of Area 3 (Auckland). It indicates that the performance of DTRM in the three subregions is almost the same, without much fluctuation. By contrast, the performance of PRM in the three subregions is inconsistent. The F1-score in R8 is only 82.96%, which is much lower than that of DTRM.
Overall, the F1-score of DTRM varies from region to region, mainly following the same trend as the benchmark. The majority of subregions (75%) have higher F1-scores when using DTRM compared with the situation when using PRM. Hence, the performance of DTRM has been improved to a certain extent. Figure 18 shows the difference between the F1-scores of DTRM and PRM. In the eight experimental subregions, the average F1-score difference between DTRM and PRM is 1.25 (the maximum value is 3.91% and the minimum value is 1.07%). The F1-score difference of six of these subregions is positive, with an average value of 2.11. The F1-score difference between the two subregions (R5 and R7) is negative. For the subregions with negative values, the absolute values are less than 0.02, which is in an acceptable error range. Figures A2d and A3b show that the error matches (red lines) are mainly distributed in the non-grid area of narrow areas. The reason for the precision decrease of the two cases could be explained as follows: For such narrow areas, the matching relationship established in L1 is limited. In the complex hybrid road matching, there are gaps between the pattern parameters; while the hybrid street pattern needs to be handled at the boundary between different street patterns. Such operation may affect the accuracy of matching and is one of the limitations of the proposed method in this study. During this type of operation, the trade-off between the upper level matched relationship (i.e., L2) and the similarity of this layer (i.e., L3) needs to be further studied. boundary between different street patterns. Such operation may affect the accuracy of matching and is one of the limitations of the proposed method in this study. During this type of operation, the trade-off between the upper level matched relationship (i.e., L2) and the similarity of this layer (i.e., L3) needs to be further studied. On the whole, when faced with different patterns of road networks, DTRM can achieve better performance than PRM and correctly identify the most matching relationships (89.63% on average). Figures A1-A3 show the matching results of the eight subregions for Areas 1-3, respectively, where the orange line indicates the correctly matched nodes, and the red line indicates the incorrectly matched nodes (Appendix).

Rotation Test on MMU
To verify the correctness and effectiveness of our proposed algorithm, the robustness of the MMU represented by a node-area structure in the face of rotation offset should be examined first. The experimental subregion R4 is selected to conduct the rotation test in this subsection. The similarity measures of the MMU (node-area unit) in our proposed method and the node-line unit in PRM are compared under different rotation offsets.
First, DR is defined as the degree of difference between the source element and the target matched element and that between the source element and a nonmatched candidate matching set, as shown in Formula (14). The greater DR is, the better the similarity indicator based on the MMU unit will be, and the easier it is to identify the correct matching relationship.
where represents the similarity value between the source node and its correctly matched node, represents the similarity value between the source node and the i-th candidate node, and n represents the number of candidates with the largest n matching similarity values to the source node (except for the correctly matched node).
Ten pairs of nodes in the region R3 (the location distribution is shown in Figure 19) are selected, and the difference in DR between DRTM and PRM without rotation is calculated.

Rotation Test on MMU
To verify the correctness and effectiveness of our proposed algorithm, the robustness of the MMU represented by a node-area structure in the face of rotation offset should be examined first. The experimental subregion R4 is selected to conduct the rotation test in this subsection. The similarity measures of the MMU (node-area unit) in our proposed method and the node-line unit in PRM are compared under different rotation offsets.
First, DR is defined as the degree of difference between the source element and the target matched element and that between the source element and a nonmatched candidate matching set, as shown in Formula (14). The greater DR is, the better the similarity indicator based on the MMU unit will be, and the easier it is to identify the correct matching relationship.
where sim r represents the similarity value between the source node and its correctly matched node, sim oi represents the similarity value between the source node and the i-th candidate node, and n represents the number of candidates with the largest n matching similarity values to the source node (except for the correctly matched node). Ten pairs of nodes in the region R3 (the location distribution is shown in Figure 19) are selected, and the difference in DR between DRTM and PRM without rotation is calculated. Figure 20 shows the comparison curve of DR between DRTM and PRM by using the 10 pairs of nodes. The average level of DTRM is significantly higher than that of PRM. The average DR value of DTRM is 0.64, whereas that of PRM is 0.48. This result shows that the similarity values of the MMU are higher than those of PRM. This condition is conducive to the identification of correct pairs of elements with the same name.    Figure 21 demonstrates that when the rotation angle is 0°, which means that no rotation occurs, the DR of PRM is approximately 0.5. That is, the similarity value between the node to be matched and the correctly matched one greatly differs from the similarity value between the feature to be matched and the nonmatched candidate, which is beneficial to obtaining correct matches. When the rotation angle is 15°, the DR is approximately 0.1. The similarity between the node to be matched and the correctly matched node is not as obvious as the similarity between the node to be matched and the nonmatched nodes. As the rotation angle increases (45°-180°), the DR value reaches only 0.01, and some negative values exist. This result shows that after a certain rotation, the similarity value between the node to be matched and the correct matching node is considerably low     Figure 21 demonstrates that when the rotation angle is 0°, which means that no rotation occurs, the DR of PRM is approximately 0.5. That is, the similarity value between the node to be matched and the correctly matched one greatly differs from the similarity value between the feature to be matched and the nonmatched candidate, which is beneficial to obtaining correct matches. When the rotation angle is 15°, the DR is approximately 0.1. The similarity between the node to be matched and the correctly matched node is not as obvious as the similarity between the node to be matched and the nonmatched nodes. As the rotation angle increases (45°-180°), the DR value reaches only 0.01, and some negative values exist. This result shows that after a certain rotation, the similarity value between the node to be matched and the correct matching node is considerably low  Figure 21 demonstrates that when the rotation angle is 0 • , which means that no rotation occurs, the DR of PRM is approximately 0.5. That is, the similarity value between the node to be matched and the correctly matched one greatly differs from the similarity value between the feature to be matched and the nonmatched candidate, which is beneficial to obtaining correct matches. When the rotation angle is 15 • , the DR is approximately 0.1. The similarity between the node to be matched and the correctly matched node is not as obvious as the similarity between the node to be matched and the nonmatched nodes. As the rotation angle increases (45 • -180 • ), the DR value reaches only 0.01, and some negative values exist. This result shows that after a certain rotation, the similarity value between the node to be matched and the correct matching node is considerably low compared with the similarity value between the node to be matched and the nonmatched node, and PRM has been unable to obtain a correct match.
As shown in Figure 22, the DR values of the 10 pairs of elements are above 0.5 regardless of the rotation angle, and the DR values of some element pairs are above 0.7 and 0.8. This shows that when using our proposed DTRM algorithm, compared with the similarity value of the node to be matched, the similarity value of the node matched with the correctly matched node always maintains a significant difference under different rotation angles. It shows that DTRM can identify the correct matching pairs efficiently from the perspective of the similarity measurement based on the MMU unit and has strong robustness. Furthermore, it indicates that the MMU and its similarity metrics can effectively calculate the similarity value under different rotation angles and overcome the limitation of the node-arc structure.
ISPRS Int. J. Geo-Inf. 2020, 7, x FOR PEER REVIEW 19 of 31 compared with the similarity value between the node to be matched and the nonmatched node, and PRM has been unable to obtain a correct match. As shown in Figure 22, the DR values of the 10 pairs of elements are above 0.5 regardless of the rotation angle, and the DR values of some element pairs are above 0.7 and 0.8. This shows that when using our proposed DTRM algorithm, compared with the similarity value of the node to be matched, the similarity value of the node matched with the correctly matched node always maintains a significant difference under different rotation angles. It shows that DTRM can identify the correct matching pairs efficiently from the perspective of the similarity measurement based on the MMU unit and has strong robustness. Furthermore, it indicates that the MMU and its similarity metrics can effectively calculate the similarity value under different rotation angles and overcome the limitation of the node-arc structure. compared with the similarity value between the node to be matched and the nonmatched node, and PRM has been unable to obtain a correct match. As shown in Figure 22, the DR values of the 10 pairs of elements are above 0.5 regardless of the rotation angle, and the DR values of some element pairs are above 0.7 and 0.8. This shows that when using our proposed DTRM algorithm, compared with the similarity value of the node to be matched, the similarity value of the node matched with the correctly matched node always maintains a significant difference under different rotation angles. It shows that DTRM can identify the correct matching pairs efficiently from the perspective of the similarity measurement based on the MMU unit and has strong robustness. Furthermore, it indicates that the MMU and its similarity metrics can effectively calculate the similarity value under different rotation angles and overcome the limitation of the node-arc structure.  To further understand the changes of the similarity values in terms of the basic matching unit DTRM and PRM with different rotation angles, the first pair of nodes and six nodes with the highest similarity values in the candidate matching set are selected. The results of PRM and DTRM are shown in Figures 23 and 24. The horizontal axis represents the index of nodes in the candidate set, the vertical axis represents the similarity value of the node to be matched with the node in the candidate set, and Node 3 is the node that correctly matches. To further understand the changes of the similarity values in terms of the basic matching unit DTRM and PRM with different rotation angles, the first pair of nodes and six nodes with the highest similarity values in the candidate matching set are selected. The results of PRM and DTRM are shown in Figures 23 and 24. The horizontal axis represents the index of nodes in the candidate set, the vertical axis represents the similarity value of the node to be matched with the node in the candidate set, and Node 3 is the node that correctly matches. As can be seen from Figure 23, when the road network does not have a rotational offset (the rotation angle is 0°), PRM correctly matches the node with the highest similarity value after the probability relaxation iteration. Nonetheless, the similarity value on Node 3 differs from those on other nodes greatly. When the road network rotates slightly (the rotation angle is 15°), although the similarity value on Node 3 is still the highest, the difference between Node 3 and other nodes is inconsiderable. When the rotation angle continues to increase, PRM can no longer correctly obtain the correct matching nodes, and the difference between Node 3 and other nodes tends to be gentle. This shows that in the PRM method when facing road network data with rotational offset, the similarity measure among basic matching units becomes invalid.  As can be seen from Figure 23, when the road network does not have a rotational offset (the rotation angle is 0°), PRM correctly matches the node with the highest similarity value after the probability relaxation iteration. Nonetheless, the similarity value on Node 3 differs from those on other nodes greatly. When the road network rotates slightly (the rotation angle is 15°), although the similarity value on Node 3 is still the highest, the difference between Node 3 and other nodes is inconsiderable. When the rotation angle continues to increase, PRM can no longer correctly obtain the correct matching nodes, and the difference between Node 3 and other nodes tends to be gentle. This shows that in the PRM method when facing road network data with rotational offset, the similarity measure among basic matching units becomes invalid. As can be seen from Figure 23, when the road network does not have a rotational offset (the rotation angle is 0 • ), PRM correctly matches the node with the highest similarity value after the probability relaxation iteration. Nonetheless, the similarity value on Node 3 differs from those on other nodes greatly. When the road network rotates slightly (the rotation angle is 15 • ), although the similarity value on Node 3 is still the highest, the difference between Node 3 and other nodes is inconsiderable. When the rotation angle continues to increase, PRM can no longer correctly obtain the correct matching nodes, and the difference between Node 3 and other nodes tends to be gentle. This shows that in the PRM method when facing road network data with rotational offset, the similarity measure among basic matching units becomes invalid.
In Figure 24, the six curves are almost the same with different rotation angles. The similarity value calculated on the basis of the MMU on Node 3 is the highest, and a significant difference exists between Node 3 and other candidate nodes. This result demonstrates that the MMU of the node-area structure constructed using our proposed algorithm and the similarity metric of the MMU is effective.

Rotation Test on the Road Network
In this subsection, we verify the validity and correctness of the DTRM proposed in this paper under the rotation condition based on the road network by comparing it with PRM. We perform angle rotations (e.g., 15 • , 30 • , 45 • , 90 • , and 180 • ) by using the experimental subregion of R4. The obtained results are shown in Figure 25. When the rotation angle is inconsiderably large (below 15 • ), PRM can still obtain a good matching result. As the rotation angle increases, the result of PRM drops sharply. At 90 • , the precision and recall are less than 5%. This result is consistent with the result of the previous similarity measure function. DTRM can show good performance with different rotation degrees. The precision and recall decrease only slightly when the rotation angle increases, but they are still approximately 90%. Thus, DTRM performs well when matching data with different rotation angles. In Figure 24, the six curves are almost the same with different rotation angles. The similarity value calculated on the basis of the MMU on Node 3 is the highest, and a significant difference exists between Node 3 and other candidate nodes. This result demonstrates that the MMU of the node-area structure constructed using our proposed algorithm and the similarity metric of the MMU is effective.

Rotation Test on the Road Network
In this subsection, we verify the validity and correctness of the DTRM proposed in this paper under the rotation condition based on the road network by comparing it with PRM. We perform angle rotations (e.g., 15°, 30°, 45°, 90°, and 180°) by using the experimental subregion of R4. The obtained results are shown in Figure 25. When the rotation angle is inconsiderably large (below 15°), PRM can still obtain a good matching result. As the rotation angle increases, the result of PRM drops sharply. At 90°, the precision and recall are less than 5%. This result is consistent with the result of the previous similarity measure function. DTRM can show good performance with different rotation degrees. The precision and recall decrease only slightly when the rotation angle increases, but they are still approximately 90%. Thus, DTRM performs well when matching data with different rotation angles. The matching result of subregion R3 with a rotation angle of 30° is selected to conduct further analysis. Figure 26 shows the mismatches of DTRM. The number of mismatches is increased only by one compared with the number of mismatches with no rotation ( Figure A1 in the Appendix A). This shows that DTRM can robustly handle the road network matching with rotation offset. Figure 27 shows the error matches of PRM for subregion R4 with a rotation angle of 30°. It can be found that the nodes that are incorrectly matched using PRM are concentrated near the boundary of the road network. The similarity measurement function of PRM depends on the spatial distance and topological connections of the element pairs. When the spatial distance of the elements with the same name deviates greatly due to the rotation offset, it is difficult for PRM to obtain the correct result. The matching result of subregion R3 with a rotation angle of 30 • is selected to conduct further analysis. Figure 26 shows the mismatches of DTRM. The number of mismatches is increased only by one compared with the number of mismatches with no rotation ( Figure A1 in the Appendix A). This shows that DTRM can robustly handle the road network matching with rotation offset. Figure 27 shows the error matches of PRM for subregion R4 with a rotation angle of 30 • . It can be found that the nodes that are incorrectly matched using PRM are concentrated near the boundary of the road network. The similarity measurement function of PRM depends on the spatial distance and topological connections of the element pairs. When the spatial distance of the elements with the same name deviates greatly due to the rotation offset, it is difficult for PRM to obtain the correct result.

Sensitivity Analysis of the Buffer Threshold
The Euclidean distance between features with the same name also changes with the rotation. To retrieve potential matching candidates more accurately, the buffer size needs to be increased appropriately. In this section, the sensitivity of the buffer threshold is examined by gradually increasing the buffer size as the rotation angle increases. Then the precision and recall of the two methods are calculated to reveal the impact of the buffer size on the experimental results.
As can be seen from Figure 28, for different rotation angles (e.g., 15°, 30°, 45°, 60°, 75°, and 90°), when the buffer radius is relatively small (i.e., 200 m), the precision and recall of the algorithm are relatively low. This may be when the buffer is excessively small; the elements with the same name in the target dataset are not in the buffer, resulting in a low matching rate. However, when the buffer size reaches a certain value and continues to increase, regardless of using DTRM or PRM, no evident fluctuation in accuracy and recall rates exists. This shows that the two algorithms are insensitive to the buffer size and have good robustness.

Sensitivity Analysis of the Buffer Threshold
The Euclidean distance between features with the same name also changes with the rotation. To retrieve potential matching candidates more accurately, the buffer size needs to be increased appropriately. In this section, the sensitivity of the buffer threshold is examined by gradually increasing the buffer size as the rotation angle increases. Then the precision and recall of the two methods are calculated to reveal the impact of the buffer size on the experimental results.
As can be seen from Figure 28, for different rotation angles (e.g., 15°, 30°, 45°, 60°, 75°, and 90°), when the buffer radius is relatively small (i.e., 200 m), the precision and recall of the algorithm are relatively low. This may be when the buffer is excessively small; the elements with the same name in the target dataset are not in the buffer, resulting in a low matching rate. However, when the buffer size reaches a certain value and continues to increase, regardless of using DTRM or PRM, no evident fluctuation in accuracy and recall rates exists. This shows that the two algorithms are insensitive to the buffer size and have good robustness.

Sensitivity Analysis of the Buffer Threshold
The Euclidean distance between features with the same name also changes with the rotation. To retrieve potential matching candidates more accurately, the buffer size needs to be increased appropriately. In this section, the sensitivity of the buffer threshold is examined by gradually increasing the buffer size as the rotation angle increases. Then the precision and recall of the two methods are calculated to reveal the impact of the buffer size on the experimental results.
As can be seen from Figure 28, for different rotation angles (e.g., 15 • , 30 • , 45 • , 60 • , 75 • , and 90 • ), when the buffer radius is relatively small (i.e., 200 m), the precision and recall of the algorithm are relatively low. This may be when the buffer is excessively small; the elements with the same name in the target dataset are not in the buffer, resulting in a low matching rate. However, when the buffer size reaches a certain value and continues to increase, regardless of using DTRM or PRM, no evident fluctuation in accuracy and recall rates exists. This shows that the two algorithms are insensitive to the buffer size and have good robustness. rotation angle, using a fixed buffer radius may miss some candidate matching pairs and the buffer radius for selecting candidate matching pairs should be increased. A properly increased buffer radius can effectively alleviate this problem. On the contrary, an excessively enlarged buffer radius may introduce incorrect matching pairs that are similar in structure and reduce the matching precision.

Impact of the Hierarchical Classification Threshold
As generating a hierarchical road network is a starting point of the proposed method, it is necessary to confirm the rationality of the stratification. Thus, this subsection explains how the stratification threshold is selected. In previous studies, head/tail breaks were used to stratify different levels of the road network for there are more small things than large things in living structures [22]. The stratification threshold was set to different values in accordance with the downstream application requirements, such as cartographic generalization, traffic flow predicting, and so on. In this study, the proportion of strokes selected in L1 is crucial in the hierarchical generation process. In accordance with an appropriate proportion, k is determined and the longest k strokes are selected as L1. Then the hierarchy of road networks could be stratified according to the steps in Section 2.2. Further, the effect of buffer size on the matching performance of PRM and DTRM under different rotation angles is compared. The precision of PRM decreases significantly as the rotation angle increases. For the same rotation angle, the precision improves slightly as the buffer distance increases, but it still does not reach the precision of PRM with no rotation angle. As the rotation angle increases, the precision of DTRM slowly decreases, while still maintaining a high level (i.e., approximately 80%). As shown in Figure 28c, at a certain angle of 45 • , when the buffer radius increases slightly, the precision of DTRM gradually rises to a relatively stable level. The rotation angle generally exerts different effects on the two matching methods. DTRM shows good resistance to the rotation angle. When the rotation angle increases, the influence of the rotation angle on the precision can be effectively suppressed by appropriately increasing the buffer radius. Due to the rotation angle, using a fixed buffer radius may miss some candidate matching pairs and the buffer radius for selecting candidate matching pairs should be increased. A properly increased buffer radius can effectively alleviate this problem. On the contrary, an excessively enlarged buffer radius may introduce incorrect matching pairs that are similar in structure and reduce the matching precision.

Impact of the Hierarchical Classification Threshold
As generating a hierarchical road network is a starting point of the proposed method, it is necessary to confirm the rationality of the stratification. Thus, this subsection explains how the stratification threshold is selected. In previous studies, head/tail breaks were used to stratify different levels of the road network for there are more small things than large things in living structures [22]. The stratification threshold was set to different values in accordance with the downstream application requirements, such as cartographic generalization, traffic flow predicting, and so on. In this study, the proportion of strokes selected in L1 is crucial in the hierarchical generation process. In accordance with an appropriate proportion, k is determined and the longest k strokes are selected as L1. Then the hierarchy of road networks could be stratified according to the steps in Section 2.2.
To determine a proper proportion for L1, Dr and HitRate are used to evaluate the different distribution of the target matched strokes and candidate strokes under different proportions. Dr is defined as the degree of difference in similarity between a source stroke and its matched stroke and that between the source stroke and a non-matched candidate stroke, which is represented by Formula (15).
where Max is the similarity between the matched stroke and its source stroke, which is denoted as reference value; V is the similarity between the non-matched candidate stroke and its source stroke. The larger the Dr, the more obvious the difference between the matched stroke and the non-matched one is, and the higher the matching accuracy is. In other words, a proper proportion should keep the value of Dr as large as possible.
HitRate is defined as the ratio of the candidate strokes that fall in the range of a specified Dr value to all candidate strokes under the specified conditions, which is represented in Formula (16). There are three conditions in Formula (16). Con b denotes the condition of a buffer size; Con p denotes the condition of a proportion of the longest strokes for L1; Con dr denotes the condition of a specified Dr value. M represents the number of candidate strokes that satisfy the conditions of both Con b and Con p . N represents the number of candidate strokes that satisfy the conditions of Con b , Con p, and Con dr .
Area 2 is taken as the test area. The buffer size in Con b is set to 50 m and 300 m. The proportion in Con p for L1 is set to 5%, 10%, 15%, 20%, 25%, and 30%, respectively. The distribution relationship between HitRate and Dr is illustrated in Figure 29. For example, B50_P5 represents the scenario of two combined conditions (the buffer size of Con b is 50 m and the proportion of Con p is 5%). When Con b is set to 50 m and the values of Con p are 5% or 10%, the mean value of HitRate reaches 0.92 when Dr is in the range of (0.5-1]; indicating above 90% of the candidate strokes have obvious geometrical differences with the target matched stroke. The HitRate value approaches 0 in the Dr range of 0-0.5, indicating there are very few candidate strokes that are very similar to the target matched stroke. In other words, the difference rate between the target matched strokes and the candidate strokes is significant when the proportions are 5% and 10%, which is beneficial to obtain correct matching results. When the values of Con p are larger than 10%, the HitRate value decreases in the Dr range of 0.5-1 and increases in the Dr range of 0-0.5. When Con b is set to 300 m, the distribution trend is similar to that when Con b is set to 50 m, indicating that the buffer size has little influence on HitRate. Therefore, 10% is suggested as the selection proportion of the first level L1 to ensure the robustness and accuracy of the matching. distribution of the target matched strokes and candidate strokes under different proportions. Dr is defined as the degree of difference in similarity between a source stroke and its matched stroke and that between the source stroke and a non-matched candidate stroke, which is represented by Formula (15).
where Max is the similarity between the matched stroke and its source stroke, which is denoted as reference value; V is the similarity between the non-matched candidate stroke and its source stroke. The larger the Dr, the more obvious the difference between the matched stroke and the non-matched one is, and the higher the matching accuracy is. In other words, a proper proportion should keep the value of Dr as large as possible.
HitRate is defined as the ratio of the candidate strokes that fall in the range of a specified Dr value to all candidate strokes under the specified conditions, which is represented in Formula (16). There are three conditions in Formula (16). Conb denotes the condition of a buffer size; Conp denotes the condition of a proportion of the longest strokes for L1; Condr denotes the condition of a specified Dr value. M represents the number of candidate strokes that satisfy the conditions of both Conb and Conp. N represents the number of candidate strokes that satisfy the conditions of Conb, Conp, and Condr.
Area 2 is taken as the test area. The buffer size in Conb is set to 50 m and 300 m. The proportion in Conp for L1 is set to 5%, 10%, 15%, 20%, 25%, and 30%, respectively. The distribution relationship between HitRate and Dr is illustrated in Figure 29. For example, B50_P5 represents the scenario of two combined conditions (the buffer size of Conb is 50 m and the proportion of Conp is 5%). When Conb is set to 50 m and the values of Conp are 5% or 10%, the mean value of HitRate reaches 0.92 when Dr is in the range of (0.5-1]; indicating above 90% of the candidate strokes have obvious geometrical differences with the target matched stroke. The HitRate value approaches 0 in the Dr range of 0-0.5, indicating there are very few candidate strokes that are very similar to the target matched stroke. In other words, the difference rate between the target matched strokes and the candidate strokes is significant when the proportions are 5% and 10%, which is beneficial to obtain correct matching results. When the values of Conp are larger than 10%, the HitRate value decreases in the Dr range of 0.5-1 and increases in the Dr range of 0-0.5. When Conb is set to 300 m, the distribution trend is similar to that when Conb is set to 50 m, indicating that the buffer size has little influence on HitRate. Therefore, 10% is suggested as the selection proportion of the first level L1 to ensure the robustness and accuracy of the matching.

Conclusions
This study has established a novel methodology for hierarchical road network matching based on Delaunay triangulation (DTRM) and a new problem formulation to conquer the problem of road network matching with nonsystematic bias in different coordinate systems. In accordance with the human visual hierarchy and regionalized cognitive mechanism, the entire road network is divided into three levels. The upstream matched relationship provides a stable contextual reference for the matching analysis of downstream road entities, thereby improving the matching accuracy. In particular, in such a road network matching task, a new MMU is designed instead of the traditional "node-arc" unit. On the basis of the rotation-invariant characteristics of triangles, the similarity metrics based on MMU help solve the problem of the rotation angle of nonskeletal road network matching. Compared with the traditional matching methods based on local features, DTRM is not limited to the specific quantization of road geometry and topological properties and always maintains better global hierarchical features throughout the matching process. This study opens up a new avenue for hierarchical road network-matching methods along different tracks.
Experiments are conducted using real datasets of Wuhan, China and Auckland, New Zealand to verify the efficacy of our approach (DTRM). We compare our proposed method with the PRM method. The results are summarized as follows: (1) When faced with different patterns of the road network (e.g., sparse or dense and grid or radial), DTRM can achieve higher precision than the benchmark PRM and correctly identify the most matching relationships (89.63% on average). (2) In terms of DR indicators, the average DR value of DTRM is significantly higher than that of PRM. In the case of rotation, the similarity metrics based on MMU have better discriminating power than the traditional index in PRM and are more conducive to the identification of correct pairs of elements with the same name. The results also show that the structure of MMU and its corresponding similarity metric can effectively overcome the rotation problem and have strong robustness. (3) Comparison of P and R indices of DRTM and PRM shows that PRM can obtain good matching when the rotation angle is inconsiderably large (below 15 • ). As the rotation angle increases, the P and R values of PRM drop sharply. By contrast, the accuracy and recall of DTRM decrease only slightly, although both remain at approximately 90%. DTRM can effectively solve the problem of road matching with different rotation angles, considering that the similarity index of MMU is insensitive to increased rotation angle.
The model developed in this study is suitable for the integration and update of urban road networks and their applications in navigation systems. The methodology also has the potential to be applied to other urban sectors for infrastructure applications with sector-oriented modifications, such as public management [27], urban planning [28,29] and transport modeling [30].
In future work, several issues need to be studied in depth. First, at the junction of irregular and regular areas, the robustness based on MMU structure matching needs further consideration. Second, the time cost is noticed less than the accuracy of the road network matching methods. The performance issue becomes increasingly prominent as the volume of matching data increases. Third, introducing a learning method based on basic matching pairs is an interesting topic to solve the road network matching problem.