Trajectory data road network construction method based on multistage grid characteristics

Abstract

The invention discloses a trajectory data road network construction method based on multilevel grid characteristics, which comprises the following steps: firstly, preprocessing original trajectory data to remove noise and abnormal values, then constructing a multilevel grid based on cleaned data, and calculating internal and neighborhood characteristics of the multilevel grid to mine semantic information of the trajectory data; a key grid is automatically identified by using supervised learning methods such as a random forest; and finally, realizing accurate extraction of the road network by combining an improved four-neighborhood refinement algorithm with abnormal structure processing. According to the method, high-precision key grid classification and recognition are realized through a multi-level semantic feature mining technology and an improved four-neighborhood refining algorithm in combination with a random forest model, so that the integrity and precision of road network extraction are effectively improved, and a road network abnormal structure is successfully eliminated; and a digital road network construction solution with strong adaptability, high accuracy and good robustness is provided for walking and vehicle trajectory data.

Description

Technical Field

[0001] The present invention belongs to the technical field of intelligent transportation, and particularly relates to a method for constructing a road network of trajectory data based on multi-level grid features. Background Art

[0002] With the acceleration of the urbanization process and the continuous expansion of population growth, urban road data is of great significance for promoting the construction of an urban comprehensive road traffic system and improving the connectivity and penetration of the infrastructure in the metropolitan area. The road network is a basic element of geographic information systems and intelligent transportation, and has a positive significance for improving the quality of community life and promoting urban and rural development and construction. The road network can reflect the development trend of a city. An accurate road network is the basis for wide applications and an important basis for urban structure and road traffic system planning. Therefore, constructing a high-precision and high-quality road network is a key research content in the field of traffic geographic information.

[0003] Crowdsourced trajectory data has richer social attributes, semantic information, and temporal information, and at the same time has a larger data volume. YANG X, Zhang Yunfei, etc. extracted the road network by analyzing the semantic information contained in GNSS trajectory data. Using GNSS trajectory data, based on the Delaunay triangulation theory, Yang Wei and Ai Tinghua respectively extracted the road centerline and road boundary, which can adapt to trajectory data with large density differences. Zheng Tianjing used pedestrian trajectory data, constructed a trajectory data density map, combined with Morse theory, and constructed a pedestrian road network by extracting the "ridge line" in the density map. Guo YJ adjusted the trajectory to construct a more compact density to extract the pedestrian road network. On this basis, Zhou Baoding fully explored the physiological characteristics of pedestrian gait and proposed pedestrian dead reckoning (PDR) to construct an indoor pedestrian road network. YANG L, AI M L constructed a human-flow probability field (HFPF) and could effectively extract the main pedestrian paths and smaller branches by combining with the method of hydrological analysis.

[0004] However, the above methods for constructing a road network by mining trajectory data mainly construct the road network based on the method of trajectory point density, and rarely consider the multi-level semantic features of trajectory data, such as the direction of trajectory points and the semantic features between trajectory segments. The extraction effect of the road network in areas with sparse trajectory point distribution is greatly affected, and the characteristics of crowdsourced trajectory data, such as different qualities, uneven coverage, and data noise, all affect the final extraction result of the road network. Summary of the Invention

[0005] The purpose of the embodiment of the present invention is to provide a method for constructing a road network of trajectory data based on multi-level grid features. By matching trajectory points into the constructed grid, semantic information including features such as density and direction in the multi-level grid is mined. At the same time, in combination with the context relationship, by identifying key grids, the road network is extracted using the key grids.

[0006] To solve the above technical problems, the technical solution adopted by the present invention is a method for constructing a road network of trajectory data based on multi-level grid features, which is specifically carried out according to the following steps:

[0007] S1. Preprocess the original trajectory;

[0008] S2. Construct a multi-level grid based on the preprocessed original trajectory, determine the internal and neighborhood features of the grid to mine the semantic information of the trajectory data;

[0009] S3. Detect the key grids based on supervised learning to determine the actual road network image;

[0010] S4. Accurately extract the road network based on the improved four-neighborhood refinement algorithm and abnormal structure processing.

[0011] Further, the specific steps of S1 are as follows:

[0012] S101. Define each GPS trajectory segment as T j ={P1,P2,…,P i}, where j ∈ (1,…,J), J is the number of trajectory segments, i ∈ (1,…,I), I is the number of trajectory points included in the trajectory segment, and P i is the i-th point of the j-th trajectory segment, where:

[0013] P i =(t i ,x i ,y i ) (1)

[0014] Among them, t i is the timestamp, and x i , y i are longitude and latitude respectively;

[0015] S102. Delete the duplicate and redundant points of the trajectory segment, and at the same time split the trajectory segments with abnormal distances and abnormal speeds;

[0016] S103. Determine the direction vector of the trajectory point, calculate according to the temporal characteristics of the pedestrian trajectory points in the trajectory segment, and the direction feature vector d i of P i is:

[0017] d i =(xi+1 -x i ,y i+1 -y i ) (2)

[0018] Integrate the trajectory points, construct and compile grid index information, and define the set of all grids as G = {g1, g2, …, g n}, where n is the total number of grids, g n is the divided grid, and define the trajectory points in g n as where n ∈ (1, …, N), N is the total number of grids constructed, and a, b, j ∈ (1, …, J), where a, b, j are trajectory segment indices.

[0019] Furthermore, the specific steps of S2 are as follows:

[0020] S201. Compile grid indices for the trajectory points and determine the internal characteristic indices of multi-level grids based on these indices;

[0021] S202. Determine the grid neighborhood characteristic indices.

[0022] Furthermore, the internal characteristic indices of the multi-level grids in S201 include the convex hull area and perimeter, the point density within the convex hull, the number of grid trajectory points, the area of the minimum circumscribed circle of the grid trajectory points, the perimeter of the minimum circumscribed circle of the grid trajectory points, the area of the minimum circumscribed rectangle of the grid trajectory points, the HMR index, the HMC index, the number of grid trajectory point direction clusters, and the similarity of the trajectory segments.

[0023] Furthermore, the specific method for determining the internal characteristic indices of the multi-level grids is as follows:

[0024] S2011. Construct a convex hull based on the trajectory points within the grid. Traverse the sorted points from left to right, connect the points that satisfy the cross product ≤ 0 in sequence, remove the intermediate points that do not meet the conditions to construct the lower convex hull. Subsequently, traverse the sorted points from right to left, connect the points that satisfy the cross product ≤ 0 in sequence, and also remove the intermediate points to construct the upper convex hull. Merge the lower convex hull and the upper convex hull, and remove the duplicate points at the beginning and end to obtain the set of trajectory points that finally form the convex hull. Then determine the convex hull area S hull (g n ) and the convex hull perimeter L hull (g n ); and determine the point density Den hull (g n ) within the convex hull based on the convex hull area S hull (g n ):

[0025]

[0026] Among them, count(P gn ) is the number of trajectory points in grid g n ;

[0027] S2012. Construct the minimum circumscribed circle of the grid trajectory points, and determine the area S mincir (g n ) and the perimeter L mincir (g n ) of the minimum circumscribed circle;

[0028] S2013. Based on the convex hull constructed in S2011, use the rotating calipers algorithm to traverse all possible rotation angles, determine the rectangle with the smallest area, and thus obtain the minimum circumscribed rectangle of the grid trajectory points, and determine the area S minretc (g n ) of the minimum circumscribed rectangle of the grid trajectory points;

[0029] S2014. Determine the HMR index and HMC index based on the convex hull area S hull (g n ), the area S mincir (g n ) of the minimum circumscribed circle of the grid trajectory points, and the area S minretc (g n ) of the minimum circumscribed rectangle of the grid trajectory points:

[0030]

[0031] S2015. Cluster the direction vectors of the grid trajectory points: For the data set P gn , if there are at least MinPts samples in the ε-neighborhood of the sample P j , then the sample P j is a core point:

[0032]

[0033] Among them, P i ∈P gn , P gn is the set of trajectory points located in the grid with index n, N ε (P j ) is the number of sub-sample sets in the ε-neighborhood that contain samples in the sample set P gn and have an angle difference not greater than ε with P j , ε is the angle difference threshold. If N ε (P j )≥MinPts, then P j is a core object and is thus clustered into a cluster; if N ε (P j) < MinPts, then P j is a non-core object. Proceed to the next step of judgment. If P j is within the ε-neighborhood of a certain core point, then assign P j to the cluster where the core point is located. If P j is not within the ε-neighborhood of any arbitrary core point, then mark P j as a noise point. angle(P i , P j ) is the function for calculating the angle difference of trajectory points:

[0034]

[0035] are respectively the polar angles after converting the sample trajectory point and other trajectory points from the rectangular coordinate system to the polar coordinate system:

[0036]

[0037] where y Pi , x Pi , y Pj , x Pj are respectively the longitude and latitude coordinates of point P i , P j ; Determine the number of grid trajectory point direction clusters by controlling the angle threshold and the number of samples;

[0038] S2016. Determine the similarity of trajectory segments in the grid based on the dynamic time warping algorithm:

[0039] First, determine two ordered trajectory segments T1 = {p1, p2,..., p m}, T2 = {q1, q2,..., q n}, where m and n are the number of trajectory points in the two trajectory segments, and m, n > 1. The similarity metric formula for the two trajectory segments is:

[0040] DTW(T1, T2) = f(m, n) (10)

[0041] where f(m, n) is the distance accumulation formula:

[0042]

[0043] Starting from the end of the trajectory segment, use the Euclidean distance between two points as the accumulation of the warping path distance. ||·|| is the Euclidean distance between two points. ||p m - q n || is the Euclidean distance between the m-th point p m of T1 and the n-th point p nThe Euclidean distance between them is then used, and the minimum cumulative distance among the three possible paths that reach the current distance is selected through the minimum function min{}.

[0044] Furthermore, the specific steps of S202 are as follows:

[0045] S2021: Based on the density of the grid's own trajectory points, calculate the density difference between the grid and the grids in the eight adjacent directions. If there are no adjacent grids, assign a density difference of 0 to that direction.

[0046] S2022: Taking the grid as the center, calculate the absolute value of the density difference between two adjacent grids that are on the same straight line as the central grid, then add them up and take the average value as the direction density difference.

[0047] Furthermore, in the automatic detection of key grids based on supervised learning in S3, a random forest is selected as the training model; when constructing the feature engineering, each grid is represented by the grid center, and the internal and neighborhood features of the grids determined in S2 are input into the random forest model; in the classification, after calculating the multi-level grid features, the grids are defined as key grids and non-key grids, and are assigned 1 and 0 respectively. The key grids and non-key grids are defined by judging whether the grid center is on the road buffer zone, and the actual pedestrian road network image is determined by using the OSM road network and the method of visual manual interpretation of high-precision remote sensing images.

[0048] Furthermore, the specific steps of S4 are as follows:

[0049] S401: Preprocess the road network image obtained in S3: By calculating the sum of the connected component values and removing the pixel points whose connected component values are less than a specific threshold.

[0050] S402: Perform morphological closing operations on the preprocessed road network image: Use the structural element to perform dilation and erosion operations on the target image in turn.

[0051] S403: Refine the road network image after the morphological closing operation using an improved four-neighborhood binary image thinning algorithm to obtain a single-pixel road skeleton line.

[0052] S404: Use Kalman filtering to perform smooth fitting processing on the pixel center points of the single-pixel road skeleton line, and then connect the smooth fitting points of adjacent pixels to generate the final road network diagram.

[0053] Furthermore, the improved four-neighborhood binary image thinning algorithm in S403 is specifically as follows:

[0054] (1) For square regions and continuous square regions, the square regions are reduced in dimension to connected triangles, and then the connected triangles are further processed. First, find the regions that are squares globally. Using the method of neighborhood perception, by judging whether the sum of the pixel neighborhood values is less than 1, the pixels are then removed, reducing the dimension to connected triangles, and then the connected triangles are processed;

[0055] (2) For connected triangles of different structures, first judge whether they are four different structures of connected triangles that meet the conditions. If they meet the conditions, check the outer neighborhood pixels of the pixels that make up the connected triangle. If the sum of the outer neighborhood pixel values is less than 1, then assign the pixels that make up this connected triangle to 0 and remove them from the final road skeleton line;

[0056] Finally, extract the single-pixel road network from the remaining connected triangles through the four-neighborhood binary image thinning algorithm.

[0057] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0058] The present invention realizes the effective analysis of the deep information of trajectory data through multi-level semantic feature mining technologies such as trajectory similarity analysis and grid trajectory point density calculation, and combines the random forest model to perform high-precision classification and recognition of key grids (the accuracy rate of the road area under the 200×200 grid reaches 80% and the recall rate reaches 79.5%). The method of the present invention has strong migration ability and can adapt to pedestrian and vehicle trajectory data; the method of the present invention replaces the traditional buffer method through the analysis of grid context relationship, significantly improving the integrity and accuracy of road network extraction, and the road length and precision rate indicators are better than the traditional methods; the present invention improves the four-neighborhood thinning algorithm, successfully eliminating abnormal road network structures such as connected triangles and square pixels, ensuring the topological quality. The present invention forms a technical closed-loop in multiple links such as multi-source trajectory semantic mining, key area recognition, and road network optimization, providing a solution that takes into account adaptability, accuracy, and robustness for the construction of high-precision digital road networks. Brief Description of the Drawings

[0059] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0060] Figure 1 It is a schematic diagram of the direction clustering of trajectory points within the grid in this embodiment;

[0061] Figure 2 It is a visualization diagram of the DTW principle in this embodiment;

[0062] Figure 3 is the calculation diagram of the neighborhood feature index of this embodiment;

[0063] Figure 4 is the schematic diagram of the random forest principle of this embodiment;

[0064] Figure 5 is the schematic diagram of the abnormal structure of this embodiment;

[0065] Figure 6 is the determination diagram of four different structure connected triangles and neighborhood pixels of this embodiment;

[0066] Figure 7 is the feature visualization diagram of this embodiment; among them, (a) is the number diagram of trajectory points, (b) is the HMC index diagram, (c) is the convex hull area diagram, (d) is the number diagram of trajectory point clusters, among which, (e) is the trajectory point density diagram, (f) is the DTW index diagram;

[0067] Figure 8 is the diagram of the original trajectory points and the prediction results of this embodiment; among them, (a) is the diagram of the original trajectory points, (b) is the prediction result diagram;

[0068] Figure 9 is the comparison diagram before and after the processing result of the isolated pixel points of this embodiment; among them, (a) is before processing, (b) is after processing;

[0069] Figure 10 is the comparison diagram before and after the closing operation processing of this embodiment; among them, (a) is before processing, (b) is after processing;

[0070] Figure 11 is the partial enlarged view of the extraction result of this embodiment. Detailed implementation manners

[0071] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0072] This embodiment provides a method for constructing a road network of trajectory data based on multi-level grid features. By constructing a multi-level grid feature index, deeply mining the deep information contained therein, and using a random forest model, inputting the relevant feature information into the random forest model for model training to identify key grids, then using the trajectory data of other regions for model verification, and finally based on the training results, using morphological methods to construct the road network.

[0073] In some specific embodiments, the specific process of the road network construction method for trajectory data based on multi-level grid features is as follows:

[0074] S1. Preprocessing of original trajectories

[0075] Pedestrian trajectory data is the trajectory of a pedestrian moving along a pedestrian path recorded by a GPS receiver device, recording the position information of the pedestrian's walking state, including longitude, latitude, and timestamp. Due to the delay in starting the GPS device and weak signals in complex situations, etc., in this embodiment, the original trajectory data is preprocessed through the following process:

[0076] Define each GPS trajectory segment as T j ={P1, P2, …, P i}), where j ∈ (1, …, J), J is the number of trajectory segments, i ∈ (1, …, I), I is the number of trajectory points included in the trajectory segment, and P i is the i-th point of the j-th trajectory segment. Among them, P i is defined as:

[0077] P i =(t i , x i , y i ) (1)

[0078] Among them, t i is the timestamp, and x i , y i are longitude and latitude respectively.

[0079] (1) Delete the duplicate and redundant points of the trajectory segment, and at the same time split the trajectory segments with abnormal distances and abnormal speeds. In some specific embodiments, the duplicate and redundant points and the trajectory segments with abnormal distances and abnormal speeds are screened by constructing an Euclidean distance matrix and a speed matrix.

[0080] (2) Calculate the direction vector of the trajectory points, which is calculated according to the temporal sequence characteristics of the pedestrian trajectory points in the trajectory segment. The calculation formula for the direction feature vector of P i is as follows:

[0081] d i =(Δx i , Δy i )=(x i+1 -x i , y i+1 -y i ) (2)

[0082] Integrate the trajectory points, construct and compile grid index information, and define the set of all grids as G = {g1, g2, …, g n}, where n is the total number of grids, and gn For the divided grid, define the trajectory points in g n as where n ∈ (1, …, N), N is the total number of grids constructed, a, b, j ∈ (1, …, J), J is the number of trajectory segments, and a, b, j are the indices of the trajectory ends.

[0083] S2. Calculate based on multi - level grid features

[0084] Trajectory data contains rich semantic information, and the semantic information contained in trajectory data is multi - level. There are correlations between trajectory points, and there are also interaction correlations between trajectory segments. Therefore, this embodiment adopts a calculation method based on multi - level grid features. Through multi - level grid design, on the basis of ensuring the global nature of trajectory data, internal grid features and grid neighborhood features are proposed to fully explore the semantic information contained in trajectory data.

[0085] S201. Compile a grid index for the trajectory points to match the trajectory points P n in g gn and calculate the multi - level internal grid feature index based on this.

[0086] To match the trajectory points P n in g gn , first compile a grid index for the trajectory points. Then, based on these indices, calculate the multi - level grid feature index.

[0087] (1) Construct a convex hull according to the trajectory points in the grid. Traverse the sorted points from left to right, and successively connect the points that meet the condition of "not turning right" (i.e., cross - product ≤ 0), and remove the intermediate points that do not meet the conditions to construct the lower convex hull. Subsequently, traverse the sorted points from right to left, successively connect the points that meet the condition of "not turning right", and also remove the intermediate points to construct the upper convex hull. Merge the lower convex hull and the upper convex hull, and remove the repeated points at the beginning and end to obtain the set of trajectory points that finally form the convex hull. The area S hull (g n ) and the perimeter L hull (g n ) of the convex hull are used as feature indices. At the same time, calculate the point density Den hull (g n ) of the trajectory points within the convex hull as a feature index. The point density calculation formula based on the convex hull is as follows:

[0088]

[0089] where count(P gn ) is the number of trajectory points in the grid g n , and at the same time count(P gn) It is also one of the characteristic indices of the grid.

[0090] (2) Construct the minimum circumscribed circle of the grid trajectory points, and calculate the area S of the minimum circumscribed circle of the grid trajectory points mincir (g n ), in this embodiment, the Welzl algorithm is used to construct the minimum circumscribed circle. Its core idea is to gradually construct the minimum circle through recursion, combine geometric properties and randomization strategies, and efficiently solve the problem of the minimum circumscribed circle. Construct the minimum circumscribed circle of the trajectory points in the grid, where the area S of the minimum circumscribed circle mincir (g n ) and the perimeter L of the minimum circumscribed circle mincir (g n ) are used as one of the characteristic indices.

[0091] (3) Through the constructed convex hull, use the rotating calipers algorithm to traverse all possible rotation angles to find the rectangle with the smallest area, so as to obtain the minimum circumscribed rectangle of the grid trajectory points, that is, traverse each side of the convex hull as the direction of one side of the rectangle. There is only one minimum rectangle that can contain all the scattered points through each side. By comparing the areas of the rectangles constructed by each side of the convex hull, the rectangle that meets the minimum area is the minimum circumscribed rectangle of the point set. In this embodiment, the area of the convex hull formed by the trajectory points is used as the actual area of the target object. After constructing the minimum circumscribed rectangle of the trajectory points in the grid, the area S of the minimum circumscribed rectangle of the grid trajectory points minretc (g n ) is used as one of the characteristic indices.

[0092] (4) Based on the convex hull area S hull (g n ), the area S of the minimum circumscribed circle of the grid trajectory points mincir (g n ) and the area S of the minimum circumscribed rectangle of the grid trajectory points minretc (g n ), construct two characteristic indices, named HMR index and HMC index respectively. Their calculation formulas are as follows:

[0093]

[0094] (5) Through the improved DBSCAN clustering algorithm, cluster the direction vectors of the trajectory points. The general DBSCAN clustering algorithm directly uses the Euclidean distance between two points for neighborhood comparison. In this embodiment, since it is clustering the direction of the trajectory points, the angle difference (in degrees) between two direction vectors is used for clustering. For a certain data set P gn , if there are at least MinPts samples in the ε-neighborhood of the sample P j , then the sample P j is called a core point, that is:

[0095]

[0096] where P i ∈P gn , P gn is the set of trajectory points located in the grid index n, N ε (P j ) is the number of subsets of samples in the ε-neighborhood containing the sample set P gn whose angular difference from P j is no greater than ε, denoted as |N ε (P j )|, and ε is the angular difference threshold. If N ε (P j ) ≥ MinPts, then P j is a core object and is thus clustered into a cluster. If N ε (P j ) < MinPts, then P j is a non-core object and needs to be further judged. If P j is located within the ε-neighborhood of a certain core point, then P j is assigned to the cluster where the core point is located. If P j is not located within the ε-neighborhood of any arbitrary core point, then P j is marked as a noise point. angle(P i , P j ) is a function for calculating the angular difference of trajectory points, and the calculation principle is as follows:

[0097]

[0098] and are the polar angles after converting the sample trajectory point and other trajectory points from the rectangular coordinate system to the polar coordinate system, respectively. The polar angle calculation formula is as follows, and the unit is converted from radians to degrees:

[0099]

[0100] where y Pi , x Pi , y Pj , x Pj are the longitude and latitude coordinates of points P i and P j , respectively.

[0101] Take the number of clusters of the grid trajectory point directions as one of the feature indices. The clustering of the feature point direction vectors is as Figure 1As shown, the number of grid trajectory point direction clusters is determined by controlling the angle threshold and the number of samples. In this embodiment, the angle threshold is set to 30°, and the sample number threshold is 5.

[0102] (6) Trajectory segment similarity, which is used to describe the degree of chaos of the trajectory segments within the grid. Since each trajectory segment in the grid is divided into different lengths and the number of trajectory points that make up the trajectory segment also varies, therefore, in this embodiment, based on the Dynamic Time Warping (DTW) algorithm, to adapt to the problem of different lengths and numbers of trajectory points of the trajectory segments, the sum of the distances between similar points in two trajectory segments is calculated through DTW, which is called the warping path distance, to measure the similarity between the trajectory segments. The smaller the warping path distance, the more similar the two trajectory segments are. Conversely, the two trajectories are not similar, as Figure 2 shown. In the figure, the difference in values between two points is used to replace the distance.

[0103] The Dynamic Time Warping algorithm is based on two ordered trajectory segments T1 = {p1, p2,..., p m} and T2 = {q1, q2,..., q n}, where m and n are the number of trajectory points in the two trajectory segments, and m, n > 1. The similarity metric formula for the two trajectory segments is:

[0104] DTW(T1, T2) = f(m, n) (10)

[0105] In the formula, f(m, n) is the distance accumulation formula, and the calculation principle is as follows:

[0106]

[0107] Starting from the end of the trajectory segment, in this embodiment, the Euclidean distance between two points is used as the accumulation of the warping path distance. ||·|| is the two-norm of the two-point coordinates, that is, the Euclidean distance. In the formula, ||p m -q n || is the Euclidean distance between the m-th point p m of T1 and the n-th point p n of T2. Subsequently, the minimum cumulative distance among the three possible paths to reach the current distance is selected through the minimum function min{}. After calculating the similarity between all pairs of trajectory segments in the grid, the degree of chaos of the trajectory segments within the grid is measured by calculating the percentage of the warping path distances between all pairs of trajectory segments in the grid that are greater than a certain threshold.

[0108] S202. Determine the grid neighborhood feature index

[0109] Since the road network has global connectivity as a whole, as Figure 3, after constructing the internal feature index of the grid, this embodiment proposes a density-based grid neighborhood feature index in combination with the grid context relationship.

[0110] (1) Neighborhood density difference: Based on the density of the grid's own trajectory points, calculate the density difference between the grid and the grids in the eight adjacent directions. If there is no adjacent grid, assign the density difference in this direction as 0, which serves as the eight feature indices of the grid.

[0111] (2) Direction density difference: Taking the grid as the center, calculate the absolute value of the density difference between two adjacent grids that are on the same straight line as the central grid, and then add them up and take the average value as the direction density difference.

[0112] This embodiment mines the spatial semantic information of trajectory data through grid features: Inside the grid, the complexity of the activity range is evaluated through the convex hull area and perimeter (a large convex hull area reflects dispersed activities, and a small area corresponds to dense aggregation), the trajectory point density characterizes the regional activity level (high density indicates busy roads or public places), the area of the minimum circumscribed circle / rectangle reveals the distribution pattern characteristics (such as the combination of a large circumscribed circle and a small convex hull implies multiple hot spot distributions), the HMR / HMC composite index depicts the spatial distribution pattern, the direction vector clustering identifies the dominant movement trend (such as pedestrian / vehicle direction preferences), and the DTW trajectory similarity quantifies the movement consistency (a high score indicates a complex path, and a low score reflects an orderly flow); at the neighborhood level, the eight-direction neighborhood density difference detects the functional boundary and the hot spot transition zone, and the direction density difference analyzes the spatial trend characteristics of the main roads. By integrating the above-mentioned internal geometric features, motion pattern parameters, and neighborhood correlation relationships of the grid, the system constructs an interpretation framework from micro motion behaviors to macro path structures, accurately identifies semantic entities such as high-frequency activity areas and core traffic corridors, and provides quantitative support for the refined road network extraction and urban spatial function analysis based on trajectory data.

[0113] S3. Key grid detection based on supervised learning to establish the actual road network

[0114] In the selection of the training model for supervised learning, use Figure 4 Random Forests as the training model. As an ensemble learning method, Random Forests has good adaptability to complex data and noise and has good generalization ability.

[0115] When constructing the feature engineering, each grid is represented by the grid center, and a total of 18 feature indices such as the corresponding grid direction class number, grid point density, trajectory similarity, 8 neighborhood density differences, and 4 direction density differences are input into the Random Forest model.

[0116] In classification, after calculating the multi-level grid features, the grids are defined as key grids and non-key grids, and are assigned values of 1 and 0 respectively. The key grids and non-key grids are defined by judging whether the grid center is on the road buffer. The OSM road network and the method of visually interpreting high-precision remote sensing images are used to determine the actual pedestrian road network image.

[0117] S4. Road network connection based on morphology

[0118] The prior art usually refines based on the four-neighborhood of pixels, which has certain limitations. Therefore, this embodiment proposes a binary image refinement algorithm based on the improved four-neighborhood for extracting the single-pixel road network.

[0119] This embodiment first preprocesses the original image, screens and removes the isolated pixel points in the original image. Utilizing the characteristics of the continuity and integrity of the road network, this embodiment uses the method of the sum of connected component values for screening, and eliminates the pixel points whose sum of connected component values is less than a certain threshold. The sum of connected component values is the sum of all binary image values of this part.

[0120] Then, the initial road network is extracted by using rasterization methods such as closing operation, erosion, and thinning. First, considering the possible breakpoint phenomenon, the closing operation is directly performed on the original image to obtain the result. The morphological closing operation is to perform dilation and erosion operations on the target image in turn using the structuring element. Dilation is to perform a logical AND operation on each pixel in the original image using the structuring element. If all are 0, otherwise it is 1. The erosion operation is similar to the dilation operation. In the stage of performing the logical AND operation, if all are 1, the image is 1, otherwise it is 0.

[0121] In the refinement stage of the road network, the general parallel iterative refinement algorithm based on the four-neighborhood quickly screens candidate points through the four-neighborhood and the ray length, and combines the 0-1 pattern number to verify the connectivity to prevent breakpoints. However, for road networks with different abnormal structures (such as Figure 5 ), such as connected triangles, square grids, continuous square grids, etc. To solve this phenomenon, a binary image refinement algorithm based on the improved four-neighborhood is proposed. On the original basis, the following two methods are proposed:

[0122] (1) For square regions and continuous square regions, through dimensionality reduction processing, that is, reducing the square region to a connected triangle, and then further processing the connected triangle. First, find the regions that are squares globally, and use the neighborhood perception method. By judging whether the sum of the pixel neighborhood values is less than 1, the pixels are then removed to reduce them to connected triangles, and then the connected triangles are processed. Among them, for the top-left pixel of the square region, its neighborhood pixels are the left, upper pixels, and the top-left pixel; for the top-right pixel of the square region, its neighborhood pixels are the right, upper pixels, and the top-right pixel; for the bottom-left pixel of the square region, its neighborhood pixels are the left, lower pixels, and the bottom-left pixel; for the bottom-right pixel of the square region, its neighborhood pixels are the right, lower pixels, and the bottom-right pixel. (2) For connected triangles with different structures, first judge whether they are four different structures of connected triangles that meet the conditions. If they meet, check the outer neighborhood pixels of the pixels that make up the connected triangle. If the sum of the outer neighborhood pixel values is less than 1, it means that this connected triangle may not be part of the real road network, so it is assigned a value of 0, that is, it is removed from the final road skeleton line. For the outer neighborhood pixels, different constructed connected triangles have different outer neighborhood pixels. Taking the connected triangle with the right-angle pixel at the bottom-left as an example, for the right-angle pixel, its outer neighborhood pixels are the left, lower, and bottom-left pixels; for the pixel above the right-angle pixel, its outer neighborhood pixels are the left, upper pixels, and the top-left and top-right pixels; for the pixel to the right of the right-angle pixel, its outer neighborhood pixels are the right, lower pixels, and the bottom-right and top-right pixels. As Figure 6 shown.

[0123] Finally, a single-pixel road skeleton line is obtained. Since the connection between pixels does not represent the real road network, smoothing fitting processing is required. Taking the center point of the pixel as the operation object, Kalman filtering is used for smoothing fitting operation. Connect the points obtained after smoothing fitting by Kalman filtering in the way of connecting adjacent pixels to obtain the final road network.

[0124] The experimental data used in this embodiment include the walking trajectory data on the campus of the Yuehai Campus of Shenzhen University, including longitude coordinates, latitude coordinates, and timestamps, the high-resolution remote sensing images of the Yuehai Campus of Shenzhen University, and at the same time, the vehicle trajectory data in Wuhan and the OSM urban road network in Wuhan are selected.

[0125] In this embodiment, the walking trajectory data collected from the Yuehai Campus of Shenzhen University is used for method training and testing. By calculating the multi-level grid feature index, the passing characteristics and density aggregation degree of the trajectory points in the grid are measured, and the visualization of its characteristics is as Figure 7 (a) - (f) shown.

[0126] The grid sizes are set to 100×100, 200×200, 300×300, 400×400, and 600×600 respectively for classification tests with different granularity divisions. Through visual interpretation of high-resolution remote sensing and combined with the road network of the Yuehai Campus of Shenzhen University on Tianditu, the training samples are labeled with category tags. Based on the above grid divisions of different sizes, the random forest algorithm is used for model training, and the ratio of the training set to the test set is set to 8:2, obtaining the classification results in Table 1.

[0127] Table 1 Classification accuracy statistics under different grid divisions

[0128]

[0129]

[0130] Experiments show that in this embodiment, when the grid is divided into 200×200, the average accuracy rate of the road area calculated reaches 80%, and the average recall rate reaches 79.5%. At this time, the effect of the model in this embodiment is the best.

[0131] To verify the transferability of the method in this embodiment, vehicle trajectory data in different regions are selected for experiments. As Figure 8 (a) shows, Wuhan is selected as the experimental area, and the vehicle trajectory data of Wuhan are input into the trained model. The key grids are identified by the model, and the identification results are shown in 8(b), where the red ones are the grids identified as key grids by the model, and the green ones indicate the grids identified as non-key grids by the model.

[0132] After completing the road network area identification, morphological processing is carried out, and then the road network is extracted. The correctly identified road grid rasters are converted into binary images.

[0133] First, the original image is preprocessed. For the isolated pixel points in the original image, screening and removal are carried out. In this embodiment, the method of the sum of connected component values is used for screening, and the pixel points with the sum of connected component values less than a certain threshold are removed. The sum of connected component values is the sum of all binary image values in this part, as Figure 9 (a)(b) show.

[0134] Then, rasterization methods such as closing operation, erosion, and thinning are used to extract the initial road network. First, considering the possible break point phenomenon, the closing operation is directly used on the original image to obtain the result. The morphological closing operation is to perform dilation and erosion operations on the target image in turn using the structural element. Dilation is to perform a logical AND operation on each pixel in the original image using the structural element. In this embodiment, the closing operation uses a cross-shaped structural element with a length of 9.

[0135] As shown in the figure, where Figure 10 (a) is the original image, Figure 10(b) is the image after performing closing operation.

[0136] Finally, perform a thinning operation on the image. Use a binary image thinning algorithm based on the improved four-neighborhood to thin the graph, and at the same time use Kalman filtering to smooth it.

[0137] Figure 11 It is a partial enlarged view of the extraction result in Wuhan, showing the extracted roads and the original trajectory points. For the Wuhan area, the extracted roads have a better fitting effect compared to the original trajectory.

[0138] In order to evaluate the road extraction effect and the accuracy of the extracted roads, in this embodiment, a 7m buffer is established based on the real roads, and the extracted roads falling within this buffer are the correctly extracted roads. The real roads are obtained through OSM, and the extraction quality is evaluated according to the following three indicators: (1) Precision P, that is, the ratio of the length of the correctly extracted roads to the total length of the extracted roads; (2) Recall R, the ratio of the number of road points identified by the random forest within 4 meters of the real road network buffer to the number of road points identified by the random forest; (3) The total length of the correctly extracted roads.

[0139] Taking the OSM road data as a reference, the comparison of the evaluation indicators of the road extraction by the model in this embodiment and the rasterization method of JIANG Y J et al. is shown in Table 2.

[0140] Table 2 Comparison of evaluation indicators of different road extraction methods

[0141]

[0142] From the above road extraction comparison experiment results and the comparison results of evaluation indicators, it can be concluded that when extracting the road network from the same trajectory data, the method adopted in this embodiment has improved the accuracy of road extraction. At the same time, for the total length of the correctly extracted roads, compared with the rasterization method of JIANG Y J et al. which directly performs buffer analysis on the trajectory data and then directly converts it into raster data, the method adopted in this embodiment matches the trajectory points to the grid, analyzes the detailed information contained in the trajectory points within the grid, and analyzes the context relationship of the grid, so as to perform more effective analysis of the road network information and thus extract a more complete road network.

[0143] Each embodiment in this specification is described in a related manner. For the same or similar parts between each embodiment, reference can be made to each other. Each embodiment focuses on the differences from other embodiments. In particular, for the system embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and for the related parts, reference can be made to the partial description of the method embodiment.

[0144] The above are only the preferred embodiments of the present invention and are not intended to limit the protection scope of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention are all included within the protection scope of the present invention.

Claims

1. A method for constructing a road network of trajectory data based on multi-level grid features, characterized in that, The specific steps are as follows: S1. Preprocess the original trajectory; S2. Construct a multi-level grid based on the preprocessed original trajectory, determine the internal and neighborhood features of the grid to mine the semantic information of the trajectory data; S3. Detect the key grids based on supervised learning to determine the actual road network image; S4. Accurately extract the road network based on the improved four-neighborhood refinement algorithm and abnormal structure processing.

2. The road network construction method according to claim 1, wherein The specific steps of S1 are as follows: S101. Define each GPS track segment as T j = {P1, P2, …, P i}, where j ∈ (1, …, J), J is the number of track segments, i ∈ (1, …, I), I is the number of track points included in the track segment, and P i is the i-th point of the j-th track segment, where: P i = (t i , x i , y i ) (1) where t i is the timestamp, and x i , y i are the longitude and latitude respectively; S102. Delete the duplicate redundant points of the trajectory segment, and at the same time split the trajectory segments with abnormal distances and abnormal speeds; S103. Determine the direction vector of the trajectory point, which is calculated based on the temporal characteristics of the pedestrian trajectory points in the trajectory segment. The direction feature vector d of P i is i as follows: d i = (x i+1 - x i , y i+1 - y i ) (2) Integrate the trajectory points, construct and compile grid index information, and define the set of all grids as G = {g1, g2, …, g n}, where n is the total number of grids, and g n is the divided grid. Define the trajectory points in g n as where n ∈ (1, …, N), N is the total number of grids constructed, and a, b, j ∈ (1, …, J), and a, b, j are trajectory segment indices.

3. The road network construction method according to claim 1, characterized in that The specific steps of S2 are as follows: S201. Compile grid indexes for the trajectory points, and determine the internal feature indexes of the multi-level grid based on these indexes; S202. Determine the grid neighborhood feature index.

4. The road network construction method according to claim 3, characterized in that The internal feature indexes of the multi-level grid in S201 include the convex hull area and perimeter, the point density within the convex hull, the number of grid trajectory points, the area of the minimum circumscribed circle of the grid trajectory points, the perimeter of the minimum circumscribed circle of the grid trajectory points, the area of the minimum circumscribed rectangle of the grid trajectory points, the HMR index, the HMC index, the number of direction clusters of the grid trajectory points, and the similarity of the trajectory segments.

5. The road network construction method according to claim 4, characterized in that The specific method for determining the internal feature indexes of the multi-level grid is as follows: S2011. Construct a convex hull based on the trajectory points within the grid. Traverse the sorted points from left to right, and successively connect the points that satisfy the cross product ≤ 0, removing the intermediate points that do not meet the conditions to construct the lower convex hull. Subsequently, traverse the sorted points from right to left, successively connect the points that satisfy the cross product ≤ 0, and also remove the intermediate points to construct the upper convex hull. Merge the lower convex hull and the upper convex hull, removing the repeated points at the beginning and end to obtain the set of trajectory points that finally form the convex hull. Subsequently, determine the convex hull area S hull (g n ) and the convex hull perimeter L hull (g n );And based on the convex hull area S hull (g n ) determine the point density Den hull (g n ) within the convex hull range: where count(P gn ) is the number of trajectory points in grid g n ; S2012. Construct the minimum circumscribed circle of the grid trajectory points and determine the area S of the minimum circumscribed circle of the grid trajectory points mincir (g n ) and the perimeter L of the minimum circumscribed circle mincir (g n ); S2013. Based on the convex hull constructed in S2011, use the rotating calipers algorithm to traverse all possible rotation angles, determine the rectangle with the smallest area, thereby obtaining the minimum bounding rectangle of the grid trajectory points, and determine the area S of the minimum bounding rectangle of the grid trajectory points minretc (g n ); S2014, based on the convex hull area S hull (g n ) the area S of the minimum circumscribed circle of the grid trajectory points mincir (g n ) and the area S of the minimum circumscribed rectangle of the grid trajectory points minretc (g n ) to determine the HMR index and the HMC index: S2015. Cluster the grid trajectory point direction vectors: For the dataset P gn , if there are at least MinPts samples in the ε-neighborhood of the sample P j , then the sample P j is a core point: Among them, P i ∈P gn , P gn is the set of trajectory points located in the grid index n, N ε (P j ) is the number of subsets of samples in the ε-neighborhood containing the sample set P gn whose angular difference from P j is no greater than ε. ε is the angular difference threshold. If N ε (P j ) ≥ MinPts, then P j is a core object and is thus clustered into a cluster; if N ε (P j ) < MinPts, then P j is a non-core object and the next judgment is made. If P j is located within the ε-neighborhood of a certain core point, then P j is assigned to the cluster where the core point is located. If P j is not located within the ε-neighborhood of any arbitrary core point, then P j is marked as a noise point. angle(P i , P j ) is the function for calculating the angular difference of trajectory points: They are respectively the polar angles after the sample trajectory points and other trajectory points are converted from the rectangular coordinate system to the polar coordinate system: Among them, are the longitude and latitude coordinates of points P i and P j respectively; the number of grid trajectory point direction clusters is determined by controlling the angle threshold and the number of samples. S2016. Determine the similarity of the trajectory segments in the grid through the dynamic time warping algorithm: First, determine two ordered trajectory segments T1 = {p1, p2, …, p m}, T2 = {q1, q2, …, q n}, where m and n are the number of trajectory points in the two trajectory segments, and m, n > 1. The similarity metric formula for the two trajectory segments is as follows: DTW(T1,T2) = f(m,n) (10) where f(m,n) is the distance accumulation formula: Starting from the end of the trajectory segment, the Euclidean distance between two points is used as the accumulation of the distance of the rectified path. ||·|| is the Euclidean distance between two points, and ||p m -q n || is the Euclidean distance between the m-th point p m of T1 and the n-th point p n of T2. Subsequently, the minimum cumulative distance among the three possible paths to reach the current distance is selected through the minimum value function min{}.

6. The road network construction method according to claim 4, wherein The specific steps of S202 are as follows: S2021. Calculate the density difference between the grid and the eight adjacent grids in different directions based on the density of the grid's own trajectory points. If there is no adjacent grid, assign the density difference in this direction as 0; S2022. Take the grid as the center, calculate the absolute value of the density difference between the two adjacent grids on the same straight line as the center grid, then add them up and take the average value as the direction density difference.

7. The road network construction method according to claim 1, characterized in that In the automatic detection of key grids based on supervised learning in S3, the random forest is selected as the training model; when constructing the feature engineering, each grid is represented by the grid center, and the internal and neighborhood features of the grid determined in S2 are input into the random forest model; in the classification, after calculating the features of the multi-level grid, the grid is defined as a key grid and a non-key grid, and are assigned 1 and 0 respectively. Determine the actual pedestrian road network image by judging whether the grid center is on the road buffer zone and using the OSM road network and visually interpreted high-precision remote sensing images.

8. The road network construction method according to claim 1, characterized in that The specific steps of S4 are as follows: S401. Preprocess the road network image obtained in S3: Calculate the sum of the connected component values and remove the pixel points with the sum of the connected component values less than a specific threshold; S402. Perform morphological closing operation on the preprocessed road network image: Use the structural element to perform dilation and erosion operations on the target image in turn; S403. Refine the road network image after the morphological closing operation using the binary image refinement algorithm based on the improved four-neighborhood to obtain the single-pixel road skeleton line; S404. Use Kalman filtering to perform smoothing fitting on the pixel center points of the single-pixel road skeleton line, and then connect the smoothed fitting points of adjacent pixels to generate the final road network diagram.

9. The road network construction method according to claim 8, wherein The four-neighborhood binary image thinning algorithm described in S403 is specifically as follows: (1) For the square region and the continuous square region, reduce the square region to a connected triangle, and then further process the connected triangle. First, find the region that is a square globally, and use the neighborhood-aware method. By judging whether the sum of the pixel neighborhood values is less than 1, the pixel is then removed, reducing it to a connected triangle, and then the connected triangle is processed; (2) For connected triangles with different structures, first judge whether they are four different types of connected triangles that meet the conditions. If they do, check the outer neighborhood pixels of the pixels that make up the connected triangle. If the sum of the outer neighborhood pixel values is less than 1, then assign the pixels that make up this connected triangle to 0 and remove them from the final road skeleton line; Finally, extract the single-pixel road network from the remaining connected triangles through the four-neighborhood binary image thinning algorithm.

Images